基于FLAC3D的倾倒变形岩质滑坡地震反应分析

郭 翔，孙 旭

山东电力工程咨询院有限公司，山东 济南 250013

0 引言

地震荷载作用下岩质滑坡的稳定性分析是岩土工程界的重要研究方向之一，涉及工程地质、岩石力学、地震力学等多个学科，目前主要的研究方法有拟静力法、动力分析法、试验法等[1-2]。

动力分析法于20世纪60年代开始研究并发展起来，目前通用的软件是FLAC3D，它基于拉格朗日有限差分法，由于其强大的计算功能和模拟功能，在岩土工程界得到了广泛的应用[3-4]。

我国西南高山峡谷地区地质作用强烈，岩石变形破坏严重，发育了大量的倾倒变形岩体，岩体在自重及构造应力等多重因素的影响下，发生弯曲或弯折[5-7]。同时，西南地区是我国地震活动频发的地区，近年来随着西南地区大型水电站的建设，地震作用下边(滑)坡的稳定性问题日益突出[8]，倾倒变形岩体在地震动荷载作用下的稳定性变化及二次变形的发展趋势目前鲜有研究。基于此，论文以澜沧江上游某在建水电站一处倾倒变形岩质滑坡为例，用FLAC3D对其在地震荷载下的反应进行了动力模拟和分析，该成果对完善地震引发滑坡灾害的研究具有一定的理论和实际意义[9]。

1 滑坡工程地质概况

澜沧江上游某在建水电站位于青藏高原东部边缘地带，此区域为典型的高山峡谷地貌，区域地质构造活动明显，倾倒变形岩体广泛发育[10]。其中澜沧江右岸发育有H6滑坡。

滑坡形态明显，坡面平缓，坡度多在30°～40°之间，滑坡堆积体周界以外大片基岩裸露，岩性为侏罗系板岩，倾倒变形发育强烈(见图1)。滑坡厚度约10～60 m，初估方量约720万m3。受公路开挖的影响，公路岸边小规模的次级滑坡发育，长约100 m。坡上植被较稀疏，岸坡结构为反向坡。

图1 滑坡上游边界侧冲沟Fig.1 Gully at upstream boundary of landslide

滑坡的前缘高程为1 408 m，后缘高程为1 650 m。滑坡前缘临界正常蓄水1 408 m，滑坡后缘高于正常蓄水1 408水位线242 m。在堆积体的前缘，县级公路穿过，路面高程约为1 407～1 409 m，公路后坡基岩大范围出露，倾倒变形强烈。

2 力学参数及地震加速度

H6滑坡可分为4个部分。其中表层部分为倾倒溃屈岩体，岩体强烈倾倒折断、坠覆，无法识别岩层层状结构，岩体风化非常严重，分布范围有限；其下为倾倒松动岩体，岩体强烈倾倒折断，整体张裂松弛，局部架空，张裂隙张开度几厘米至十厘米左右，张裂隙有碎石或黏土填充，岩体风化较为严重，厚度约5～10 m；再下部为倾倒松弛岩体，岩体倾倒较为强烈，岩体不仅发生层内张裂变形，局部也发生切层剪张破裂，张开宽度几毫米左右，岩体风化较为严重，厚度约8～15 m；底部为倾倒蠕变岩体，岩体倾倒变形较弱，层内错动带剪切错位，层内岩体微量张裂变形，风化程度比较弱，此类型倾倒体发育深度大(4种类型的倾倒岩体见图2)。

图2 倾倒变形岩体的类型Fig.2 Types of toppling deformation rock mass

其相应物理力学参数见表1：

表1 岩体物理力学参数Table.1 Physical and mechanical parameters of rock mass

区域范围内地震活动强烈，其中，大于7级的地震发生过6次，6～6.9级地震30余次。工程区50年超越概率10%的地震动峰值加速度为131.9 gal，相应地震基本烈度为Ⅶ度。50年超越概率5%、100年超越概率2%的地震动峰值加速度分别为166.3 gal、257.8 gal。

3 滑坡模型的建立

3.1 模型及参数设定

模型以边坡走向NE方向作为X轴，NW方向作为Y轴，竖直方向为Z轴。以实际高程坐标为竖直方向坐标，建立直角坐标系。根据现场的地质资料，H6模型采用了13个剖面。模型地面高程1 320 m，顶面高程最高1 720 m。边坡地面高程1 385 m，坡顶高程最高1 720 m，边坡沿走向方向长度为577.1 m，沿倾向坡长551.4 m。模型自下而上分为4层，共有7 082个单元，1 856个节点。

3.2 边界条件设置

由于倾倒变形岩体存在裂隙，滑坡边界出现的波的反射将会减小，传播到边界上的能量大部分都会被吸收，因此模型底部采用黏性边界，四周采用自由场边界[11]。根据研究，当入射角小于60°时，其吸收效果较为理想。法向黏阻力tn和切向黏阻力ts的计算公式为：

tn=-ρCPvn

(1)

ts=-ρCsvs

(2)

式中：vn、vs为边界法向和切向速度分量；Cn、Cs为P波和S波的传播速度；ρ为介质的密度。

FLAC3D在模型的侧面生成一维和二维的网格以实现自由场条件，通过阻尼器与自由场网格进行耦合，自由场的不平衡力转移到主体网格的边界上，因此向上的面波在边界上不发生扭曲[12]。

3.3 阻尼的选择

对于动力问题的计算，FLAC3D中常用的阻尼形式有3种，分别为瑞利阻尼、滞后阻尼和局部阻尼[13]。大量的模拟结果证明，局部阻尼较为适合简单问题的求解；滞后阻尼使用中限制多；瑞利阻尼的理论研究较多，其计算结果与实际规律较为吻合，因此论文使用瑞利阻尼。对于瑞利阻尼所要使用的最小阻尼比和最小中心频率两个参数，根据经验岩土体的最小阻尼比一般为0.02～0.05，本次计算模型中的倾倒变形体其值取0.04；最小中心频率是根据速度时程曲线，调整fmin使得能量主要集中在fmin～3fmin范围内，此时的fmin即为最小中心频率。

3.4 人工地震波的合成

目前工程计算中常用到的地震波形式主要有直接记录到的地震波、规范要求的地震波和人工合成的地震波等[14]。论文根据工程地震学的相关理论人工合成地震波，根据场地的区域地震活动，100年超越概率2%的地震动峰值加速度为257.8 gal。库区地震持续时间30 s，以0.02 s的时间间隔采样，共计1 500个数据点，如图3所示，人工合成的最大加速度为258 cm/s2，地震峰值约为0.26 g的地震水平加速度。

图3 人工合成100年超越概率2%地震波加速度时程曲线Fig.3 Synthetic acceleration time history curve of seismic wave with 2% exceeding probability in 100 years

图3中的加速度时程需要由数值积分公式转变为速度时程，见图4。

图4 地震波速度时程曲线Fig.4 Seismic wave velocity time history curve

根据速度时程曲线，求得相应的应力时程，在模拟地震波对边坡的作用时输入。应力时程公式可参照文献[15]。

4 模拟结果与分析

用FLAC3D对边坡进行地震动荷载的模拟，需要先进行静力分析，在此基础上再将地震波加入。模拟分为两步，先输入边坡自重应力，再输入地震动[16-18]。图5和图6为边坡在自重应力下的塑形区和水平方向位移云图，由图可知，天然状态下塑形区零星地分布于坡体的浅表层，高程范围为1 410 m～1 480 m与1 600 m～1 650 m之间，浅表层水平位移多在2 cm～5 cm之间，仅极少数地形凸出部分位移较大，这说明边坡发生大规模剪切破坏的可能性很低，即使发生也只在小部分区域，如发生掉块或滚落等，不会危及整个滑坡的稳定性。

图5 天然状态坡体塑性区分布图Fig.5 Distribution of plastic zone of natural slope

图6 天然状态坡体剪切应变分布图Fig.6 Shear strain distribution of natural slope

图7 地震作用下坡体塑性区分布图Fig.7 Distribution of plastic zone of slope under earthquake action

图8 地震作用下坡体剪切应变分布图Fig.8 Shear strain distribution of slope under earthquake action

图7和图8为边坡在地震作用下的塑形区和水平方向位移云图，由图可知，地震作用下塑形区较天然状态下面积增大，深度也增大，高程范围为1 410 m～1 530 m与1 580 m～1 700 m之间，而水平位移显著增大，多为10 cm～15 cm之间，局部可达20 cm以上。这说明在地震作用下，倾倒变形岩体组成的边坡比较容易破坏，且破坏位置多集中于倾倒程度严重的区域，即倾倒溃屈部位，这些部位较为破碎，由于受到地震动荷载的作用，更易产生位移。同时在坡体表面，坡脚和坡肩位置应力最为集中，随着地震作用的持续，极易在坡肩与坡脚之间形成贯通的破裂面，从而导致边坡发生整体失稳，坡脚将被剪切破坏，坡肩则被拉裂破坏。

5 结语

(1)倾倒溃屈较为破碎的岩体结构对边坡的稳定性非常不利。

(2)在地震作用下，倾倒岩体边坡表层位移较天然状态下明显增大，深部岩体无明显的位移，在水电站建设中，需要着重考虑对倾倒变形程度严重的坡体进行防护，防止发生大规模的失稳破坏。