地震作用下边坡及滑坡稳定性及破坏过程分析

李兆宇, 马宗源, 魏睿真, 焦 凯

(1.中国水利水电第三工程局有限公司, 陕西 西安 710024; 2.西安理工大学 岩土工程研究所, 陕西 西安 710048)

地震是影响边坡稳定性的主要外部因素之一,而强震引发的山体滑坡是地震诱发的主要震害类型之一,对地震后期的灾后重建工作造成巨大影响。通过总结分析2008年汶川地震中形成的滑坡特征及分布情况,发现地震诱发滑坡的破坏力及危害性极大,有些甚至超过了地震本身所造成的损失[1]。

滑坡稳定性分析方法主要由土力学中的边坡稳定分析方法改进发展而来。上世纪初提出的圆弧滑动面的极限平衡方法,用来分析静力情况下边坡稳定问题,后来基于土条间力平衡假定提出了条分方法[2-3]。我国学者基于二维情况下Spencer条分方法提出了一种新的三维极限平衡边坡稳定分析方法且具有一定实用性[4]。上世纪70年代研究人员提出利用强度折减的方式结合有限元方法分析边坡稳定,后来称为强度折减法[5-6]。与传统条分方法相比,有限元强度折减法不需要假定滑动面及土条间的作用力,可自动搜索边坡的滑动面并计算安全系数[7]。

目前拟静力方法是边坡地震稳定性问题的主要分析方法。该方法主要基于静力条件下边坡稳定性分析法基本思路,即通过极限平衡方法或有限元强度折减法确定边坡的地震安全系数[8-9]。该方法虽然简化了计算,但是却忽略了地震产生的惯性力影响以及边坡的动力反应过程,所以无法反映地震过程中边坡的动力响应及边坡失稳后的破坏过程。此外,拟静力方法使用恒定的体力荷载考虑地震影响,使得边坡安全系数的计算结果过小,往往偏于保守。Newmark基于刚塑性滑块假定提出了地震条件下边坡永久位移的计算公式并广泛应用于边坡地震稳定性问题的计算分析[10-11]。该方法能够合理考虑地震效应的影响,但仍然基于小变形及极限平衡假定进行计算,无法考虑土体的动力特性及分析地震过程中边坡的失稳破坏过程。此外研究人员采用分型方法考虑岩土体力学特性、边坡形态特征及地震烈度等因素来综合评价边坡的地震稳定性[12]。

本文基于显式有限元及动力学大变形方法提出一种适合边坡地震稳定性及其破坏过程的非线性分析方法[13-14],分析了地震情况下滑坡的稳定性及滑坡失稳后的破坏过程。

1 理论及方法

1.1 土动力学基本理论

采用土动力学模型分析地震情况下边坡稳定性问题。选取Hardin-Drnevich骨干加载曲线模型预测土的动模量和阻尼比随动应变的衰减关系。动三轴试验数据说明Hardin-Drnevich模型更适于预测饱和黄土的动应力应变关系[10]。根据Hardin-Drnevich模型,土的动剪切模量和阻尼比可写为动剪应变的函数[15-17]:

G/Gmax=1/(1+γ/γref)

(1)

D/Dmax=(γ/γref)/(1+γ/γref)

(2)

式中:G为土体的动剪切模量;D为土体的阻尼比;Gmax为土体的初始剪切模量;Dmax为土体的最大阻尼比;γ为土体的动剪应变幅值;γref为参考剪应变(γref=τmax/Gmax,即G/Gmax=0.5时动剪切模量对应的动剪应变幅值γ)。γref越大,土的动剪切模量随动应变衰减越快,土的阻尼越大。

三维情况下一般将γ取为广义剪应变形式,即动力过程中土体的剪切模量和阻尼比与广义剪应变相关。但是广义剪应变是恒为正的,不能反映循环加载时剪应变的正反变化,即不能准确反映反向加载时的应力路径。本文采用将循环加载过程中卸载到反向加载阶段的广义剪应变增量△γ乘以一个折减系数a,以抵消反向加载过程中广义剪应变不能为负的影响。根据黏弹性理论构建阻尼和弹簧并联型的黏弹性本构模型:

(3)

使用Visual Fortran语言编制黄土动本构关系计算程序,由VUMAT子程序接口导入ABAQUS有限元计算软件。采用一个四节点减缩积分单元计算单元的剪应力及剪应变关系,并与黄土动三轴试验测试数据进行对比。对该单元的顶部约束水平自由度,单元底部水平方向上加载正弦加速度时程曲线,正弦时程函数见式(4),峰值加速度为0.04 g。

(4)

式中:A为加速度;t为振动时间。

图1为不同类型土的动剪切模量、阻尼比与动应变的关系,其中陕北黄土数据取自文献[16],砂土和黏土数据取自文献[18]和[19]。黄土动力学参数为: 最大剪切刚度Gmax=29.2MPa, 最大阻尼比Dmax=0.159,参考剪应变γref=0.03。可以看出根据Hardin-Drnevich骨干加载曲线建立黄土动力本构模型能够很好模拟黄土加卸载过程中的滞回曲线形状的应力应变关系。

图1 土的动剪切模量、阻尼比与动剪应变的关系Fig.1 Relationship between dynamic shear modulus, damping ratio and dynamic shear strain of soil

图2显示了输入正弦幅值不同系数a取值情况下广义剪应变的时程变化曲线。如a=1.0情况下动剪切模量和阻尼比完全由广义剪应变确定,a=0情况下动剪切模量和阻尼比只与广义剪应变正向加载幅值有关,与卸载及反向加载过程无关。图3黄土的动应力与动应变关系可以看出a=1.0情况下在反向加载过程中动应力应变曲线切线斜率明显偏大,而a=0情况下的动应力应变曲线呈椭圆形状,a=0.5时动应力应变曲线与试验结果相近。

图2 广义剪应变时程曲线Fig.2 Generalized shear strain time history curve

图3 黄土动应力应变关系理论预测与试验数据的对比Fig.3 Comparison of prediction and experimental data on dynamic stress-strain relationship of loess

1.2 计算方法

使用强度折减法计算边坡的安全系数(FOS),边坡抗剪强度参数折减方式可写为[7,20]:

(5)

式中:c和φ分别为边坡土体原始的黏聚力及内摩擦角;cf和φf分别为经过折减之后的边坡土体的黏聚力及内摩擦角。边坡位移突然增大(折减系数与边坡最大位移关系曲线拐点)时刻的折减系数确定为边坡安全系数。

显式有限元方法和隐式有限元方法最大的区别在于求解非线性问题是否迭代,是否所有待求物理量在同一时刻获得解答。隐式方法处理非线性问题时一般无法直接求解,需要采用迭代方法求解雅可比矩阵,对于高度非线性问题迭代求解可能会不收敛。显式方法依靠时间积分求解控制方程,无需迭代直接进行求解,求解高度非线性问题具有一定优势。显式有限元计算过程存在波动引起求解误差,需要控制时间步长保证求解稳定性。根据节点力的平衡方程,加速度可写为:

(6)

式中:M为节点质量矩阵;P为外力;I为单元内力。位移表示的平衡方程的显式积分形式可写为:

(7)

式中:t为时间;Δt为时间增量。保证求解稳定性的时间步长由系统最高频率及系统阻尼确定,稳定步长计算如下式:

(8)

式中:ωmax为时间;ξ为系统阻尼。

选用ABAQUS有限元软件显式分析模块及大变形模式,采用二维平面应变问题假定进行边坡稳定性问题的计算分析。采用C3D8R减缩积分单元对边坡模型进行网格划分。以刚性基底的边坡算例计算静力情况下边坡的安全系数对本文方法进行对比验证。边坡坡高30 m,坡比1∶2。分别使用显式及隐式有限元方法计算边坡的安全系数,其中显式有限元方法采用大变形,隐式有限元方法采用小变形模式进行计算。边坡土体计算参数为:弹性模量E=100.0MPa,Poisson比υ=0.25,重度γ=20 kN/m3,黏聚力为c=30 kPa,内摩擦角φ=20°,剪胀角ψ=0(不考虑剪胀性的影响),重力加速度g=9.81 m/s2。图4(a)为边坡计算网格划分图,图4(b)～(c)为不同方法确定的边坡极限状态时网格变形计算结果对比。图5为不同方法确定的边坡最大位移与强度折减系数SRF的关系,给出了本文使用显式有限元方法(EFEM)及使用隐式有限元(IFEM)和显式有限差分方法(EFDM)的计算结果对比,其中隐式有限元方法采用SLOPE64边坡稳定分析程序进行计算[21],显式有限差分方法使用FLAC软件进行计算[22]。

图4 显式有限元、隐式有限元及显式有限差分法计算出的边坡安全系数及网格变形结果对比Fig.4 Comparison of safety factor and mesh deformation calculated by different methods

图5 边坡最大位移与强度折减系数SRF的关系Fig.5 Relationship between maximum displacement of slope and strength reduction factor SRF

从图5可以看出边坡最大位移随强度折减系数SRF的增加逐渐增大,当位移突然增大, 曲线出现明显拐点时确定为边坡的临界破坏状态。隐式有限元方法确定的边坡安全系数FOS=1.36,显式有限差分方法确定的边坡安全系数FOS=1.43。对比两种方法可看出,使用显式有限元及强度折减方法同样可以确定边坡安全系数,但由于显式有限元方法采用大变形模式,因此边坡位移比隐式方法要大。相比隐式方法,显式有限元方法采用显式时间积分方案求解动力学方程,材料及几何非线性问题不需要进行迭代,因此更适用于动力学及大变形问题的计算分析[23]。

2 边坡地震稳定性及破坏过程分析

2.1 计算条件

计算设置同上节,土的计算参数取值详见表1。在折减强度参数的同时将边坡土体的动力学参数也进行折减,如下式所示[17]:

(9) 表1 计算参数取值表Tab.1 Computation list of the parameter value

地震动输入选取天然地震加速度时程,分别来自日本Kobe、美国Imperial Valley、Loma Prieta、Northridge以及土耳其Kocaeli地震[24]。将边坡模型底面设置为加速度边界,为了比较地震动特性对边坡稳定性的影响,将地震输入的水平及竖向分量选取同一种加速度时程,其中竖向分量为水平向分量的2/3。边坡两侧设置零加速度边界(在该处加速度强制为零),显式计算中约束加速度同样也会约束位移自由度,因此边坡两侧边界尽量设置远离边坡坡体区域。地震动具体信息详见表2,地震加速度时程见图6。图7为5%阻尼情况下五种地震加速度时程的反应谱曲线,可以看出五个地震波中美国的Northridge地震波峰值反应加速度最大,土耳其Kocaeli地震波特征周期Tg最大。

表2 地震动时程特征信息Tab. 2 Ground motion time history characteristic information

图6 五个地震加速度时程Fig.6 Five seismic acceleration time histories

图7 五个地震时程的反应谱曲线Fig.7 Response spectrum curves of five seismic time histories

2.2 计算结果

分别采用五种地震加速度时程进行计算,得到地震之后边坡最大位移与折减系数的关系以及最终确定的安全系数量值(见图8)。从图8结果可以看出加载Kocaeli地震加速度时程确定的边坡的安全系数的数值最小(FOS=1.20),其他四种地震情况确定的边坡安全系数结果相近。从五个地震加速度时程对应的反应谱曲线可以看出,Kocaeli地震的峰值反应加速度最小,但其特征周期Tg最大,即地震峰值能量持时较其他四个地震时间长,因此在其他条件不变的情况下Kocaeli地震确定的边坡安全系数是最小的。

图8 五个地震得到的边坡最大位移与折减系数的关系Fig.8 Relationship between maximum displacement and reduction coefficient obtained from five earthquakes

图9为不同地震过程中边坡坡顶水平向位移随时间的变化,当折减系数SRF=1.21,地震时间到20 s时边坡坡体出现持续滑动破坏。

图9 地震过程中边坡坡顶水平向位移时程变化曲线Fig.9 The time-history curve of horizontal displacement of slope top during earthquake

图10为Kocaeli地震过程中边坡网格变形图,可看出当地震时间为26 s时刻边坡坡体出现滑裂面,坡脚出现破坏之后出现整体滑动破坏。

图10 Kocaeli地震过程中边坡的破坏过程Fig.10 Progressive failure of slope during Kocaeli earthquake

3 滑坡地震稳定性及破坏过程分析

滑坡地处甘肃省甘南州舟曲县,滑坡地层结构自上而下主要为第四系上更新统马兰黄土,中、上更新统滑坡堆积碎石土及残坡积碎石土,中上志留统千枚岩、板岩及中泥盆统灰岩为主。

图11所示为北山上部滑坡体地层剖面及计算模型。将该滑坡视为二维平面应变问题,分析滑坡在地震条件下的稳定性及其破坏过程。

图11 滑坡计算材料分区及单元网格划分图Fig.11 Slide calculation material partition and unit grid division diagram

边坡一般为均质土层构成且高度较低,而滑坡的地层结构更为复杂且高度较大,因此土体的刚度和强度需考虑初始应力的影响,土体最大剪切刚度由下式确定[25]:

Gmax=κpa(σm/pa)n

(10)

式中:κ为初始刚度系数;σm=(σ1+σ2+σ3)/3;pa为标准大气压;n为幂次系数。对于强风化的残坡积层土体强度参数c采用下式计算:

c=c0(σm/pa)n

(11)

式中c0为地表土体的黏聚力,随着深度增加土体强度逐渐增大直至达到基岩的强度。地震荷载分别选取土耳其Kocaeli和美国Northridge地震加速度时程作为为横向及竖向输入分量,模型两侧为零加速度边界条件。假定基岩在地震过程中不会出现塑性变形为完全弹性体。基岩上部风化残积层强度参数为初始黏聚力c0=40kPa,内摩擦角φ=30°。土动力学参数为初始刚度系数κ=800,n=0.3,最大阻尼比Dmax=0.2,参考剪应变γref=0.167。此外假定滑动面的摩擦系数μ=0.3。

图12和图13分别为折减系数SRF达到1.05后,地震过程中滑坡的广义剪应变及位移云图。从计算结果可以看出,滑坡从形成滑动面到滑动是一个完整的破坏过程。该滑坡属于牵引式滑动破坏,即前缘首先滑动破坏卸荷之后带动后缘产生滑动。滑坡体并未完全滑动到坡脚沟谷底部,而是一半停留在山坡中部,只有一部分滑入谷底(最大滑动距离为1.46 km)。通过分析可知,显式有限元及动力学大变形方法的结果可分析地震情况下诱发滑坡的破坏过程。

图13 地震过程中滑坡体的位移云图Fig.13 Contourof landslide displacement during earthquake

4 结 论

本文提出了一种地震条件下边坡及滑坡稳定性分析方法,并对边坡及滑坡的地震稳定性及破坏过程进行了计算分析。

1) 地震情况下边坡的失稳破坏是一个从坡体变形到出现局部破坏直至整体失稳滑动的发展过程,并非是达到极限状态后突然发生的,须基于土动力学理论使用动力学大变形计算方法分析该问题。

2) 通过本文研究证明,显式有限元及大变形计算方法可直接确定边坡的地震安全系数,同时还可分析地震过程中边坡的动力响应及破坏过程。相对于传统隐式分析方法,本文所提出的方法在静力及小变形条件下分析边坡稳定性具有较大的误差,但在研究地震作用下边坡的稳定性及地震诱发滑坡灾害问题具有显著优势。

3) 通过一个滑坡工程实例的分析,说明本文提出的方法同样适用于地震条件下滑坡稳定性问题的计算分析,并可进一步确定滑坡的滑动距离及破坏范围。本文未考虑降雨、地下水及三维地形等因素对边坡及滑坡稳定性的影响,下一步将建立三维地形模型考虑降雨及地下水位的影响对地震诱发滑坡灾害问题进行计算分析工作,以期为防灾减灾工作提供参考依据。