引言

近年来我国的地震频发,在地震力的作用下极易触发边坡滑坡。在地震中由于滑坡造成的次生灾害十分常见,尤其是大型和巨型滑坡造成的灾害非常严重。在发生过地震的灾区存在着数量众多的碎石土堆积体边坡,其坡面形态多种多样,多数处于亚稳定状态,当余震发生时极可能再次滑坡(黄润秋,2011)。因此,对碎石土边坡稳定性的研究具有重要意义。

地震波输入的入射角对边坡的动力响应影响较大(邓龙胜等,2012)。一般在计算边坡的地震反应时,地震波输入通常采用垂直向上入射的剪切波方式,这对于震源距边坡较远时是合理的;然而当震源距边坡较近时,地震波并非垂直入射,而是以一定的倾斜角度向上传播,这时地震动呈现出空间变化特性(苑举卫等,2010),其中一个有效的证据是利用强震记录可统计出低频段S波的入射角均值为56º左右(尤红兵等,2009)。

目前,很多学者对地震作用下边坡的稳定性进行了相关研究。如:黄博等(2013)基于半无限弹性空间地震波的传播理论,对斜入射地震波在土体中的动应力路径进行了分析;李山有等(2005)应用显示有限元法,对倾斜断层台阶在地震体波斜入射下的地震反应进行了数值模拟;何刘等(2014)采用自制的单向振动试验台,对不同坡面形态的土质边坡失稳破坏过程进行了模型试验研究;卢坤林等(2014)从失稳边坡资料统计分析、模型试验和理论计算三个方面,研究了坡面形态对边坡稳定的影响。但是,对斜入射地震波作用下的碎石土边坡以及坡面形态对边坡稳定性影响的研究还较少。

本文在前人的大量研究基础之上,基于有限元差分软件FLAC3D展开了多种坡面形态的对比分析,其中包括不同坡度坡面的对比,不同凹凸度坡面的对比和不同台阶数坡面的对比。本文得出的坡面形态对边坡稳定性的影响规律及关系,可为碎石土边坡在斜入射地震波作用下的稳定性研究提供基础数据,也为类似边坡的防治提供了相关的技术依据。

1 方案分析
1.1 数值模拟理论基础

FLAC3D采用完全非线性分析方法并基于显式差分法,使用由周围区域真实密度得出的网格节点集中质量,求解全部运动方程。采用摩尔-库伦模型,程序可以自动计算永久变形,并同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响(陈育民等,2009)。因此,FLAC3D可以模拟边坡在地震动作用下的完全非线性响应。

1.2 斜入射地震作用的实现

理论上,任何一种地震波的入射方式都可以分解为垂直和水平两个方向上的震动(陈晓利等,2011)。首先根据波动理论,将静态边界上两个方向的波动场分别分解为自由波场(已分解得到的两个方向的震动波)和散射波场(由局部场地效应引起)两部分(廖振鹏,2002)。然后,将地震波动的问题转化为波源的问题,即将地震波转化为直接作用在静态边界上法向和切向的等效荷载。

基于人工边界等效荷载法(张如林等,2010),总波场ut可分解为自由波场uf和散射波场us

ut=uf+us
(1)

由位移连续条件和力学平衡条件,静态人工边界上任意一点B的运动方程可写为:

(2)
ì
í
î
MBTu••tB(t)+CBTutB(t)+KBTutB(t)=FfBT(t)+FsBT(t)
MBNv••tB(t)+CBNvtB(t)+KBNvtB(t)=FfBN(t)+FsBN(t)

式中,MBTMBN为质量矩阵;CBTCBN分别为静态人工边界节点上平面内切向和法向的阻尼器系数;KBTMBN为刚度矩阵;F fBT(t)F fBN(t)分别为模拟自由波场需施加在静态人工边界节点上的切向和法向荷载;F sBT(t)F sBN(t)分别为模拟散射波场需施加在静态人工边界点上的切向和法向荷载,可由静态人工边界提供,故无需在边界节点上额外施加。

自由波场的模拟,需要首先设定已知的原自由波场为u0(xyt)和ν0(xyt)。其中,在静态人工边界上B点产生的切向位移和法向位移分别为u0(xByBt)和ν0(xByBt);切向应力和法向应力分别为ν0(xByBt)和σ0(xByBt),xByB为静态人工边界上节点B的坐标。

为准确实现波动输入,还需要在静态人工边界上施加连续分布的等效荷载,以使静态人工边界上产生的位移和应力与原自由波场中的位移与应力相同,即满足以下条件:

(3)
ì
í
î
u(xByBt)=u0(xByBt)
v(xByBt)=v0(xByBt)
(4)
ì
í
î
τ(xByBt)=τ0(xByBt)
σ(xByBt)=σ0(xByBt)

式中,u(xByBt)和ν(xByBt)、τ(xByBt)和σ(xByBt)分别为施加的等效荷载在静态人工边界上B点产生的位移和应力。
图 1 静态人工边界及其脱离体示意图 Fig. 1Static artificial boundary and its free body

如图1所示,根据一般力学分析中的脱离体概念,可将静态人工边界与其附加的阻尼元件脱离,则静态人工边界上B点的应力为:

τ(xByBt) A=F fBT(t)-fBT(t)
(5)
σ(xByBt) A=F fBN(t)-fBN(t)
(6)

式中,A为边界点B所代表的边界面积;F fBT(t)F fBN(t)分别为阻尼元件和静态人工边界联接处的切向荷载和法向荷载。

将式(3)代入到式(6)可得:

fBT(t)=CBTu(xByBt)=CBTu0(xByBt)
(7)
fBN(t)=CBNv(xByBt)=CBNv0(xByBt)
(8)

联立式(4)、式(7)和式(8),根据已知的自由波场求出人工边界上的速度和应力,进而分别计算出静态人工边界上B点需要施加的切向和法向等效荷载:

F fBT(t)=CBTu0(xByBt)+τ0(xByBt) A
(9)
F fBN(t)=CBNv0(xByBt)+σ0(xBBt) A
(10)

这样便实现了地震波转化为直接作用在静态边界上法向和切向的等效荷载,就可以在FLAC3D里实现地震波斜入射条件下碎石土边坡动力响应的数值模拟。

1.3 坡面形态的选择

碎石土边坡的坡面形态复杂多样,自然风化形成的边坡和地震波作用后形成的边坡均可分为四种基本的类型:凸面、凹面、顺直面和台阶面,以及它们的相互组合形成的组合形式。为了便于研究不同坡面形态对斜入射地震波作用下碎石土边坡稳定性的影响,笔者选择了以下简化对比分析:

(1)顺直型坡面的坡度对碎石土边坡稳定性的影响;

(2)顺直型(即不凹不凸)、凸面型和凹面型坡面的凹凸度,对碎石土边坡稳定性的影响;

(3)台阶型坡面的台阶数对碎石土边坡稳定性的影响。

2 算例验证及分析
2.1 算例模型

为验证前面介绍的在FLAC3D中实现斜入射地震波作用的数值模拟方法,本文的算例模型采用二维计算模型,如图2所示。算例中模型的各种介质参数为:剪切模量G为16.7MPa;弹性模量E为40MPa;泊松比V为0.2;容重D为2000kg/m3。

图 2 计算模型及监测点图 Fig. 2Computational model and monitoring site
图 3 输入位移图 Fig. 3Displacement time history of input wave

在地震反应分析时,输入的地震波通常采用的是SV波,这里仅进行SV波输入下的验证。入射波采用正弦波位移,位移波时程如图3所示,幅值为1m,持续时间为2s,波速Vs为1400m/s,分别进行0°垂直向上入射和30°斜入射。在计算模型中选取两个监测点A和B,其中A为自由表面节点,B为模型底部节点,如图2所示。

2.2 算例结果及分析

在入射角分别为0°和30°的SV波作用下,通过FLAC3D计算可得到位移数值解,根据地震波反射原理推导出位移理论解,如图4—图7所示。

图 4 SV波0°垂直入射下A点位移时程 Fig. 4Displacement time history of point A under 0° vertically incident SV wave
图 5 SV波0°垂直入射下B点位移时程 Fig. 5Displacement time history of point B under 0° vertically incident SV wave
图 6 SV波30°垂直入射下A点位移时程 Fig. 6Displacement time history of point A under 30° vertically incident SV wave
图 7 SV波30°垂直入射下B点位移时程 Fig. 7Displacement time history of point B under 30° vertically incident SV wave

从图4可以看出,SV波在0°垂直入射时,由于斜坡自由表面的反射作用,入射波和反射波在自由表面处幅值叠加,自由表面节点A的位移幅值为入射位移幅值的1.5倍,数值解和理论解符合很好。从图5可以看出,在底部节点B的位移时程中可看到两个波,其中第一个周期的正弦波为入射波,第二个周期的正弦波为从自由表面反射回来的波,它们在到达底部之后的时间内,其位移值几乎维持在0附近,即在底面静态人工边界上能够全部透射出去。这表明本文采用的静态人工边界和等效荷载输入方法具有较高的计算精度。从图6和图7可以看出,SV波在30°斜入射下恰好垂直于斜坡,入射波和反射波在自由表面处幅值叠加,自由表面节点A的位移幅值为入射位移幅值的1.73倍,底部节点B的位移幅值为入射位移幅值的0.87倍。无论是0°垂直入射,还是30°斜入射,表面A点和底部B点的数值计算结果与理论解都非常符合,这表明采用人工边界和等效荷载波动输入的方式,可以进行碎石土边坡在斜入射地震波作用下的分析。

3 数值模型的建立
3.1 模型尺寸及监测点布置

模型基本尺寸为:长250m;宽60m;高130m。为便于对比分析,可将各种坡面形态的边坡模型进行以下分组。

(1)坡度影响组

建立3个不同坡面坡度的模型,其坡度比分别为:1:1.73、1:1.3和1:1。在3个不同坡度模型的坡面上分别布置坡肩A点、坡腰B点和坡脚C点,如图8所示。

(2)凹凸型影响组

建立5个不同坡面凹凸度的模型,其分别为:凹型2个、不凹不凸型(即顺直型)1个和凸型2个。在5个不同凹凸度模型的坡面上分别布置坡肩A点、坡腰B点和坡脚C点,如图9所示。

图 8 坡度影响组模型尺寸及检测点 Fig. 8The effect of declivity on model size and monitoring site
图 9 凹凸度影响组模型尺寸及检测点 Fig. 9The effect of declivity on model size and monitoring site

(3)台阶级数影响组

建立3个不同坡面台阶数的模型,其分别为:无台阶(即顺直型)、一级台阶和两级台阶。在无台阶模型的坡面上分别布置坡肩A1点、坡腰B1点和坡脚C1点;在一级台阶模型的坡面上分别布置坡肩A2点、坡腰B2点和坡脚C2点;在两级台阶模型的坡面上分别布置坡肩A3点、坡腰B3点和坡脚C3点,如图10所示。

图 10 台阶数影响组模型尺寸及检测点 Fig. 10The effect of slope steps on model size and monitoring site
3.2 模型材料

本文的重点是探究斜入射地震波作用下坡面形态对碎石土边坡变形的影响,因此,假定边坡为均质碎石土边坡,采用吴周明等(2012)给出的物理力学参数取值,如表1所示。

表1 碎石土的物理力学参数表 Table 1 Physical and mechanical parameters of gravel-soil
3.3 模型边界条件与阻尼设置

在模型四周设置自由场边界,对生成的自由场网格进行group赋值,保证在人工边界处的传播特性与在原连续介质中保持一致。这样可较好地模拟无限域的效果,使向上的面波在边界上不会发生扭曲。在模型底部设置静态边界,但输入地震波作用时需将地震波转换为人工边界上的等效荷载输入。

在FLAC3D中有瑞利阻尼、局部阻尼和滞后阻尼。考虑到瑞利阻尼的理论与常规动力分析方法类似,而且实践证明瑞利阻尼计算得到的加速度响应规律比较符合实际。因此,本文采用瑞利阻尼,最小临界阻尼比ξmin取5%,最小中心频率ωmin取16.8。

3.4 地震波的选择与转换

本文选取汶川地震波的SV波,图11为汶川地震波加速度时程,其加速度峰值为0.3g,持时为20s,入射角为56º。按照上述方案分析中的斜入射实现方法,利用数学软件MATLAB计算出汶川地震波的等效荷载,然后编成table文件,再输入到FLAC3D中,从而实现地震波的斜入射。

图 11 汶川地震波加速度时程 Fig. 11Acceleration time history of Wenchuan earthquake
4 各组模型的计算及结果分析
4.1 坡度影响组

坡度影响组分别计算和记录了3种不同坡度的碎石土边坡沿坡面上各监测点的加速度、速度和位移。为了描述地震动力作用下碎石土边坡坡面加速度的反映规律,笔者设定坡面各监测点的加速度峰值PGA与模型底部输入地震波的加速度峰值的比值为PGA放大系数(徐光兴等,2008)。根据各监测点观测到的数据,处理后可得到不同坡度的碎石土边坡各监测点水平方向PGA放大系数关系图(见图12)和各监测点水平位移图(见图13)。

图 12 各监测点水平方向PGA放大系数 Fig. 12Horizontal PGA amplification coefficient of each monitoring site
图 13 各监测点水平位移 Fig. 13Horizontal displacement of each monitoring site

由图12可知,通过对比坡面各高度监测点水平方向的PGA放大系数后发现,水平方向PGA放大系数呈现出沿坡面随高度增加而增大的规律;通过对比不同坡度的碎石土边坡坡面的各监测点水平方向PGA放大系数后发现,PGA放大系数呈现出随坡度增加而增大的规律。

由图13可知,3个碎石土边坡坡面各监测点均有不同大小的水平位移,坡度比为1:1.73 的碎石土边坡总体水平位移在0.25—0.78m范围内;坡度比为1:1.3的碎石土边坡总体水平位移在0.5—1.75m范围内;坡度比为1:1的碎石土边坡总体水平位移在0.82—2.89m范围内。这表明在斜入射地震波的作用下碎石土边坡均不稳定,但不同坡度的碎石土边坡的不稳定程度是有差异的,它们分别表现为:坡度比为1:1.73的碎石土边坡相对较为稳定;坡度比为1:1.3的碎石土边坡局部滑动,但整体还未破坏;坡度比为1:1的碎石土边坡大范围滑动,边坡整体破坏。从碎石土边坡各监测点水平位移的大小可以得出,坡面的水平位移随高度增加而增大,随坡度增加而增大,这与坡面的PGA放大系数是相吻合的。

4.2 凹凸度影响组

凹凸度影响组分别计算和记录了5种不同凹凸度的碎石土边坡沿坡面上各监测点的加速度、速度和位移,并绘制出了各边坡的剪应变增量云图,如图14和图15所示,其中,凹型2>凹型1;凸型2>凸型1。

图 14 不同凹度边坡的剪应变增量云图 Fig. 14Shear strain increment cloud of slope with different concavity
图 15 不同凸度边坡的剪应变增量云图 Fig. 15Shear strain increment cloud of slope with different convexity

由图14可知,剪应变增量随边坡凹度增加,从顺直型的近似均布逐步沿坡面向上转移到坡肩集中分布,可以推断凹型坡的竖向位移主要集中在坡肩处。由图15可知,剪应变增量随边坡凸度增加,从顺直型的近似均布逐步沿坡面向中下转移到坡腰集中分布,可以推断凸型坡的水平位移主要集中在坡腰处。因此,凹型边坡和凸型边坡的破坏形式是不同的。

单独对比凹型1和凹型2边坡坡肩的竖向位移后发现,凹型2的终时竖向位移为2.65m,大于凹型1的终时竖向位移1.20m,如图16所示。单独对比凸型1和凸型2边坡坡腰的水平位移后发现,凸型2的终时水平位移为2.35m,大于凸型1的终时水平位移1.20m,如图17所示。由此可以得出,凹型边坡的坡肩是潜在滑坡面,凸型边坡的坡腰下方是潜在滑坡面,并且凹度越大,坡肩的竖向位移越大,边坡越不稳定;凸度越大,坡腰的水平位移越大,边坡越不稳定。

图 16 凹型坡面坡肩A点位移 Fig. 16The shift of point A on concavity slope shoulder
图 17 凸型坡面坡腰B点位移 Fig. 17The shift of point B on convexity slope waist
4.3 台阶级数影响组

台阶级数影响组分别计算和记录了3种不同台阶级数的边坡沿坡面上各监测点的加速度、速度和位移,并绘制出了各监测点位移-时间关系图,如图18所示。

图 18 台阶型各坡面监测点的位移 Fig. 18The shift of each monitoring site on each slope surface of step type slop

由图18可知,被台阶分割形成的各小坡段内,坡面位移和整坡的坡面位移变化规律大致相同,沿坡面随高度增加而增大;坡度和总高度相同但台阶数不同的3个边坡,它们各自的总体竖向位移有所差异,其中无台阶边坡的总体竖向位移最大,在1.25—0.12m范围内;一级台阶边坡的总体竖向位移次之,在.96—0.07m范围内;两级台阶边坡的总体竖向位移最小,在0.65—0.08m范围内。这表明坡度和高度相同的碎石土边坡,在斜入射地震波的作用下坡面的台阶级数越多边坡的稳定性越好。

综上笔者认为,加速度在沿坡面传递时,台阶对边坡坡面的加速度峰值放大效应有减弱的作用,从而减小了位移量。同时,台阶将一个长又高的坡面分割成了几个小坡面,降低了坡面的高度,从而减少了坡面加速度峰值的放大效应,继而减小了位移量。

5 结论

本文利用FLAC3D软件进行了各种坡面形态的碎石土边坡在斜入射地震波作用下的动力响应模拟实验,通过大量的数值计算和不同模型的对比,对斜入射地震波作用下的碎石土边坡得到了如下认识:

(1)各种坡面形态的碎石土边坡均有较大的位移,表明碎石土边坡在斜入射地震波作用下大多是不稳定的,存在滑坡的安全隐患,相关部门需重视类似碎石土边坡的防治工作。

(2)在斜入射的情况下,坡面凹凸度对碎石土边坡潜在滑坡面的位置有明显影响,且与一般垂直入射的影响不同:随着坡面凹度的增大,剪应力增量向坡肩处集中,坡肩愈加不稳定;随着坡面凸度的增大,剪应力增量向坡腰处集中,坡腰愈加不稳定。台阶对坡面的加速度放大效应有削弱作用,台阶数越多越稳定。

(3)边坡防治的几点建议:减小坡面的坡度;对凹型坡加强坡肩位置处的防治,对凸型坡加强坡腰位置处的防治;对高边坡可设置多级台阶增加稳定性。

参考文献
[1]陈晓利,李杨,洪启宇等,2011.地震作用下边坡动力响应的数值模拟研究.http://www.ysxb.ac.cn/ysxb/ch/index.aspx 岩石学报 27(6):1899—1908.[本文引用:1次]
[2]陈育民,徐鼎平,2009.FLAC/FLAC3D基础与工程实例. 北京:中国水利水电出版社 181—184.[本文引用:1次]
[3]邓龙胜,范文,2012.黄土边坡动力响应的影响效应研究. 工程地质学报 20(4):483—490.[本文引用:1次]
[4]何刘,吴光,赵志明,2014.坡面形态对边坡动力变形破坏影响的模型试验研究. 岩土力学 33(1):111—117.[本文引用:1次]
[5]黄博,凌道盛,丁浩,2013.斜入射地震波在土体中产生的动应力路径及试验模拟. 岩土工程学报 35(2):276—283.[本文引用:1次]
[6]黄润秋,2011.汶川地震地质灾害后效应分析. 工程地质学报 19(2):145—151.[本文引用:1次]
[7]李山有,马强,韦庆海,2005.地震体波斜入射下的断层台阶地震反应分析. 地震研究 28(3):277—281.[本文引用:1次]
[8]廖振鹏,2002.工程波动理论导引. 北京:科学出版社 81—92.[本文引用:1次]
[9]卢坤林,朱大勇,2014.坡面形态对边坡稳定性影响的理论与试验研究. 岩石力学与工程学报 33(1):35—42.[本文引用:1次]
[10]吴周明,傅军键,2012.碎石土边坡稳定性及影响因素的数值模拟研究. 江西理工大学学报 33(1):30—34.[本文引用:1次]
[11]徐光兴,姚令侃,李朝红等,2008.边坡地震动力响应规律及地震动参数影响研究. 岩土工程学报 30(6):918—923.[本文引用:1次]
[12]尤红兵,赵凤新,荣棉水,2009.地震波斜入射时水平层状场地的非线性地震反应. 岩土工程学报 31(2):234—240.[本文引用:1次]
[13]苑举卫,杜成斌,刘志明,2010.基于设计地震动的地震波斜入射波动输入研究. 四川大学学报(工程科学版) 42(5):250—255.[本文引用:1次]
[14]张如林,楼梦麟,2010.基于FLAC3D的斜入射地震波作用的数值模拟方法研究.http://www.tmgcxb.com/ 土木工程学报 43(S1): 22—27.[本文引用:1次]