基于PFC-Geostudio耦合的边坡降雨变形破坏研究

2022-12-30朱志男李慧亭刘太辉成昊凌

河北工程大学学报(自然科学版) 2022年4期

朱志男，李慧亭，陈羽，刘太辉，成昊凌

(1.中国建筑第八工程局有限公司，甘肃兰州 730030；2.西安交通大学人居环境与建筑工程学院，陕西西安 710049)

滑坡灾害对人类的生命财产安全造成了极大的威胁，根据国内外的滑坡记录，许多滑坡都是由降雨导致[1]。2011年“9.16”特大暴雨诱发了四川省南江县数以千计的群发性缓倾角浅层土质滑坡[2]，2013年7月，在持续二十多天的降雨作用下，陕西省延安市全市有7 594处发生滑塌，造成了42人遇难，经济损失达66.15亿元[3-4]。可见，研究降雨诱发滑坡的机理、分析边坡滑坡过程对提出边坡防灾减灾的对策具有重要意义。

目前国内外学者对降雨诱发滑坡做了大量的研究，研究方法主要为理论分析、室内试验与数值模拟[5]，由于数值模拟成本低、效果显著的特点被广泛应用[6]。许方领等通过有限元软件ABAQUS软件分析了降雨入渗后的渗流场及位移和安全系数的变化规律[7]，王磊等利用有限元软件Geostudio软件分析了降雨条件下坡体含水率和孔压响应情况，并探讨了裂缝对陡坡稳定性的影响[8]。通过有限元方法模拟渗流比较简单高效，但有限元方法难以考虑大变形，无法反映滑坡过程中裂隙的延伸过程，可以通过颗粒流方法来分析边坡的滑坡过程，近年来众多学者开始利用颗粒离散元的方法来研究边坡的稳定性，汪华安通过强度折减的方法借助颗粒流数值模拟研究了烟墩岭边坡在降雨工况下的变形破坏模式[9]，张家勇以贵阳省开阳县鱼鳅坡滑坡为研究对象，通过PFC3D分析了滑坡的破坏过程[10]；吴谦通过PFC3D对黄土边坡降雨冲刷过程进行了流固耦合模拟，研究了降雨过程中边坡遭受侵蚀程度及水流侵蚀能力的分布规律[11]；但颗粒流方法缺乏有效的本构模型来反应饱和-非饱和应力应变关系，现有的PFC-CFD耦合方法计算效率较低，且流体网格形状较为规则单调，难以应用于工程实际。

为解决有限元方法无法计算大变形以及颗粒流方法中流体-颗粒耦合计算较为复杂的问题，本文将有限元与离散元相结合，提出了一种PFC-Geostudio耦合的方法，将有限元软件Geostudio的渗流计算结果导入到颗粒离散元方法中，通过等效方法来研究边坡在降雨条件下的变形情况以及滑坡过程。

1 PFC-Geostudio边坡降雨模拟方法

降雨对边坡土体有三个方面的影响：(1)土体由天然状态转化为饱和/非饱和状态，重度增加，增大了下滑力；(2)土体强度参数减小，如摩擦角、内聚力等，增大滑坡风险；(3)土体内发生渗流，土颗粒受到渗流力的作用，当降雨强度较大时，坡表还会发生径流，对土颗粒产生拖曳力作用[12]。因此，可通过Geostudio来计算渗流场，然后将渗流计算的结果导入至PFC方法中，根据颗粒所在位置的含水率，改变其重度与接触参数，实现降雨导致土体重度和强度参数的变化模拟；对渗流区域施加渗流力，当降雨强度较大时，额外考虑拖曳力，以实现降雨对土体的力学作用的模拟；通过以上方法即可在颗粒流中等效实现降雨效果，随后可进行边坡变形与滑坡分析，具体流程可见图1。

图1 Geostudio渗流计算模型图Fig.1 Geostudio-PFC calculation method

1.1 Geostudio渗流计算

有限元软件Geostudio内置的SEEP/W模块，通过非饱和土力学计算方法，可以快速模拟边坡在降雨条件下渗流场变化情况，并得到边坡内浸润线分布以及各点的含水量。为研究在降雨条件下边坡内的渗流场，建立如图 2所示的有限元模型，模型长度为50 m，高度为20 m，单元网格长度设置为0.3 m，整个边坡模型划分为8 676个单元，8 879个节点，将模型的上边界设置为降雨边界，两侧及底部设置为零流量边界，渗流计算中所使用的参数如表1所示。

图2 Geostudio渗流计算模型图Fig.2 Geostudio seepage calculation model diagram

根据气象部门对降雨强度的划分，降雨通常可分为为小雨、中雨、大雨、暴雨、大暴雨和特大暴雨六个等级，为研究在降雨条件下边坡的渗流场情况，设定降雨强度为30 mm/24 h(大雨)，原始水位为地下30 m，连续48 h降雨后，边坡渗流场如图3(a)所示，可见沿雨水入渗方向，形成了一定厚度的饱和带，继续向下为非饱和带，最下侧土体为天然状态，未受雨水入渗影响，提取边坡前缘、边坡底部、边坡中部、边坡顶部、边坡后缘五个位置沿入渗方向的各点的含水率，形成图3(b)所示曲线图，可见对于图2所示的均质土体边坡，降雨入渗规律与土体的所在位置无关，仅与其所在深度相关，土体含水率最大的区域出现在饱和带，最大含水率为40%，饱和带的宽度约为1.6 m，在非饱和带，土体含水率沿入渗方向持续减小，非饱和带的厚度约为1.8 m，含水率最小为15.5%。

表1 Geostudio中的模型参数Tab.1 Model parameters in Geostudio

图3 降雨影响下边坡渗流场情况Fig.3 Slope seepage field under the influence of rainfall

1.2 PFC-Geostudio耦合方法

为实现PFC-Geostudio的耦合，可在PFC中建立相同尺度的模型，将Geostudio中的计算结果，即每个位置的含水率情况导入；在颗粒流模型中，根据各个颗粒的位置，修改其重度，若颗粒处于饱和区域，则颗粒的重度应取饱和重度，颗粒的接触参数也需按照饱和参数进行折减，根据文献[13]的统计结果，饱和参数可按照天然参数的15%折减，若颗粒处于非饱和区域，颗粒的重度应根据其所在位置的含水率进行折减，其接触参数也应进行相应的折减，具体可见下式：

(1)

(2)

式中：γ1—饱和区内颗粒重度，N/m3；γ2—非饱和区颗粒重度，N/m3；γ3—天然区域颗粒重度，N/m3；ω—颗粒所在位置的含水率，可从Geostudio中的计算结果获取；ω1—饱和含水率；ω3—天然含水率；f1—饱和区颗粒的接触参数(此处以摩擦系数为例)；f3—天然区域颗粒的接触参数。

除改变颗粒重度与接触参数外，在饱和区域由于雨水入渗，颗粒还会受到渗流力的作用，颗粒所受到的渗流力可简化等效为[14-16]:

J=γfiV

(3)

式中：J—渗流力，N；方向为沿雨水入渗方向；γf—流体重度，N/m3；i—流体的水力梯度，V—流体体积，m3。

当降雨强度较大，超过土体的入渗能力，则在边坡表面会发生径流。在水的冲刷作用下，坡表土体会受到冲刷作用，则土颗粒会受到拖曳力的作用，其方向为沿雨水冲刷方向，其大小为：

fdrag=f0ε-χ

(4)

式中：f0—单个颗粒所受的拖曳力，N；ε—颗粒所在流体单元的孔隙率，ε-χ为考虑局部孔隙度的经验系数。这个修正项可以使拖曳力同时适用于孔隙度较高和孔隙度较低的系统，而且流体的雷诺数也可以大幅度取值。单个圆形颗粒所受流体拖曳力为[17]：

(5)

式中：Cd—拖曳力系数；ρf—流体密度，kg/ m3；r为颗粒半径，m；u—颗粒速度矢量，m/s；v—流体速度矢量，m/s；不同位置的坡面径流速度各不相同，当水沿坡面流下时，由于重力的作用，水流的速度从坡顶到坡底持续增加，在边坡底部水流速度达到最大，当水流到达坡底时，由于土颗粒的阻碍，水在坡底流动时，其速度又会持续减小。

如图4所示，位于坡面的颗粒A处的径流速度大小与方向为:

图4 不同位置坡面径流速度示意图Fig.4 Schematic diagram of the velocity of slope runoff at different positions

(6)

(7)

(8)

式中：vA—颗粒A处的径流速度，m/s；vmax—坡表径流的最大冲刷速度，m/s；如图 4所示，h1—颗粒A距坡底的高度，m；h2—坡顶距坡底的高度，m；vAx—颗粒A处径流速度沿x方向的分量，m/s；vAy—颗粒A处径流速度沿y方向的分量，m/s；Δx—水平方向的单位增量;Δy—水平方向增量为Δx时坡面线的高程增量。

位于坡底的颗粒B处的径流速度大小为:

(9)

式中：vB—颗粒B处的径流速度，方向为沿径流的冲刷方向，m/s；x1—颗粒B距坡底的水平距离，m；x2—坡底的长度，m。

式(5)中，拖曳力系数Cd为：

(10)

式(5)中，经验系数χ采用下式进行计算：

(11)

式(10)与式(11)中，Rep表示颗粒雷诺数，并采用下式进行计算：

(12)

式中，μ—流体动力粘度，N·s/m2。

2 模型建立与实现方法

2.1 颗粒流模型建立及参数标定

在颗粒流方法中，首先导入边坡的边界生成墙体，然后在边界内生成完整的细观颗粒体系。当按照真实土颗粒的粒径生成模型时，会导致颗粒数量过多，严重影响计算效率，但如果颗粒粒径取值过大，生成颗粒太少，模型的精度也会受到影响，因此建立颗粒流模型时，选择一个合适的粒径范围至关重要，根据孙其诚[18]的研究，当模型的特征长度比值L/R≥50时，颗粒粒径对材料的力学性质影响较小，因此选取颗粒粒径为0.04～0.08 m，共生成颗粒83 771个，由于通过初始命令生成的模型比很不均匀，会出现部分区域颗粒之间空隙较大、部分区域颗粒之间应力集中的现象[19]，因此需要运用伺服机制让模型变的更加均匀，伺服后的模型及其力链分布情况如图5所示。

图5 颗粒流初始模型Fig.5 Initial model of particle flow

在颗粒流方法中，通过颗粒(ball)来模拟材料颗粒，通过接触(contact)来模拟颗粒之间的相互作用，通过赋予颗粒不同的细观参数来模拟颗粒之间不同的力学性质，因此选择合适的接触模型以及对接触赋予合适的参数至关重要[20]。

当采用线性接触模型，颗粒间的接触处于粘结状态时，接触模型是线弹性的，直到受力状态超出强度极限以及粘结被破坏才会使界面处于非粘结状态，在非粘结状态下，线性接触模型处于线弹性和摩擦状态，当剪切力超出库伦极限时就会发生滑动，该粘结模型能很好的模拟土体应力应变特性及其破坏特性，因此颗粒间的接触采用线性接触模型。

由于赋予颗粒的细观参数与材料的宏观参数之间的关系并不是一一对应的，因此为确保颗粒流模型细观参数的合理性，需要进行参数标定，即不断调试模型的细观参数，直到模型的宏观响应与我们所需的模型宏观参数相接近时，参数标定方才结束。构建如图 6(a)所示的岩土试样，试样的孔隙率和粒径取值与边坡模型取值相同，粒径取0.04～0.08 m，孔隙率取0.17。对建立的岩土试样模型进行伺服，使试样模型应力分布均匀，随后赋予接触参数，对岩土试样进行不同围压下的双轴压缩试验，提取模型的应力应变曲线(图 6(a))，绘制莫尔应力圆，如图 6(b)所示，即可得试样的宏观参数，不断调整接触的细观参数，使得试样的宏观参数为c=23 kPa，φ=14 °，边坡模型的细观参数最终取值如表2所示。

图6 参数标定图示Fig.6 Parameter calibration diagram

表2 颗粒流模型细观力学参数Tab.2 The meso-mechanical parameters of the calculation model

2.2 颗粒流模拟中降雨实现方法

经参数标定后，对建立的初始边坡模型赋予参数，施加重力，运行至平衡，随后对所有颗粒的速度和位移清零，此时的模型即为天然状态下的边坡模型。设定降雨强度为30 mm/24 h，经Geostudio渗流计算，将各点的含水率导入颗粒流模型中，遍历所有颗粒，根据颗粒所在位置的含水率按照式(1)和式(2)修改颗粒的重度与接触参数，实现效果如图7所示。

当对饱和区域的颗粒施加渗流力时，首先应筛选出受渗流力影响的颗粒。遍历所有颗粒，由于饱和土的含水率为40%，当颗粒所在位置的含水率大于39%时，即可判定其属于饱和区，会受到渗流力的作用，即可按照式(3)对颗粒施加等效渗流力。当降雨强度大于土体的入渗能力时，坡表由于雨水冲刷，坡表颗粒会受到拖曳力的作用，因此在具体施加拖曳力时，首先需要判断降雨强度是否大于土体入渗能力，仅当降雨强度大于土体入渗能力时才调用拖曳力施加函数，在进行拖曳力的具体施加时，遍历所有颗粒，判定出坡表受径流影响的颗粒，如图8(a)所示，由于拖曳力的方向为雨水径流方向，因此当颗粒位于坡底时，径流的方向为沿水平方向，颗粒位于坡面时，径流的方向为沿坡面方向，不同位置的径流速度的大小与方向可按照图4以及式(6)—(9)考虑，再通过Fish函数获取颗粒的粒径、速度、位置等信息，随后即可根据式(4)—(12)来施加拖曳力。当降雨强度为100 mm/24 h时，颗粒流方法中渗流力与拖曳力的施加效果如图8(b)所示。

3 结果分析

3.1 大雨工况下边坡变形分析

为研究在降雨条件下边坡的变形情况，设定降雨强度为30 mm/24 h，降雨持续时间为48 h，经Geostudio渗流计算，计算结果如图3(a)所示，将各个位置的含水率导入PFC中，通过1.2节所述的耦合方法来等效实现降雨效果。颗粒流边坡在降雨条件下的变形情况如图9(a)所示，可见最大位移区域位于坡顶，最大位移为8.3 cm，变形主要集中在坡顶饱和区域，土体位移方向如图9(b)所示，可见位移方向主要为土体沉降，产生沉降的原因在于降雨导致土体容重增加以及土体强度参数的减弱，此外坡面土体由于左侧无土体支撑，因此会产生一定沿坡面滑动的趋势，若降雨强度持续增加，则降雨入渗深度会继续增大，同时饱和带与非饱和带的厚度也会持续增加，会导致坡面土体稳定性会进一步降低，有潜在滑坡的趋势。

图10 降雨条件下的测点位移Fig.10 Measuring points displacement under rainfall conditions

图11 暴雨工况下边坡滑坡过程及裂隙延伸情况Fig.11 Landslide process and fissure extension of the slope under rainstorm conditions

在坡底、坡面、以及坡顶设置监测点，以进一步了解在计算过程中边坡不同位置土体的变形情况。由于单独监测一个颗粒随机性较强，可能会出现少数颗粒飞出模型边界的情况，因此选择数个颗粒作为一个测点，形成一个颗粒团，取其平均位移以反映测点所在位置土体的位移，可更准确的反映土体的变形情况，且当边坡发生失稳破坏时，测点颗粒可随之运动，依旧能反映测点所在区域边坡的变形情况，如图10(a)所示，沿坡面共设置测点15个，每个测点半径为0.3 m，平均每个测点包含颗粒20个。当降雨强度为30 mm/24 h、降雨持续时间为48 h时，经PFC-Geostudio耦合计算后，提取各测点位移，形成如图10 (b)所示曲线图，可见沿水平方向，坡底和坡顶各测点位移变化不大，但坡顶位移整体大于坡底，坡面各测点的位移变化较大，各测点的位移随高程增加而增加，最大位移出现在坡顶位置。

3.2 暴雨工况下边坡滑坡过程分析

为研究边坡在暴雨工况下的变形破坏情况，设定降雨强度为100 mm/24 h，降雨持续时间为48 h，此时由于降雨强度远大于土体的入渗速度，因此在坡面会产生径流，需额外考虑拖曳力对坡表颗粒的影响，在暴雨工况下，拖曳力的施加效果如图8(b)所示，经PFC-Geostudio耦合计算，边坡的变形破坏过程如图11所示，边坡的滑动过程可划分为3个阶段。

第一阶段如图11(a)—(b)所示，坡面颗粒由于雨水冲刷作用，在坡脚首先发生破坏，底部颗粒在拖曳力的作用下脱离边坡，沿坡面滑动至坡底，同时边坡内饱和区域形成大量裂隙，但并未滑动；

第二阶段如图11(b)—(d)所示，由于坡脚颗粒冲刷滑离边坡，导致上部土体失去支撑，首先在坡脚发生小范围滑坡，随后裂隙拓展延伸，从坡脚已滑动区域的顶部延伸至坡顶，形成大范围滑坡，在滑动过程中，滑体内部出现大量裂隙，滑体逐渐破碎为散体；

第三阶段如图11(d)—(g)所示，边坡后缘土体由于失去前侧土体的支撑，内部出现大量裂隙，在自重作用下失稳，出现了多次小规模的滑动，边坡滑坡的最终状态如图11(h)所示，坡面土体大范围滑动，滑体破碎并堆积至坡底。

为确定在暴雨条件下坡面不同位置土体的滑坡时序，提取在滑坡起始阶段中测点6、测点8、测点10的速度信息，如图12所示，可见由于坡面径流冲刷作用的影响，位于坡脚的测点6速度最先增加，且速度较大，当坡脚颗粒滑动后，上部土体逐步变形，位于坡面中部的测点8的速度开始逐步增加，随后测点10的速度随之增加，但其速度小于测点6，可见测点8附近土体先于测点10附近土体滑动，滑坡过程为明显的渐进式牵引式滑坡。

图12 滑坡过程各测点位移Fig.12 Displacement of each measuring point during the landslide process

文献[21]搭建室内降雨滑坡框架模型，通过在模型顶部设置喷头来模拟降雨，试验中的滑坡过程如图13(a)—(d)所示，在降雨冲刷作用下表面土体流失，从坡脚处开始发生破坏，产生滑坡临空面，逐步诱发了大规模滑坡，并引发了数次小规模滑坡，滑坡过程为典型的牵引式滑坡，与上述采用耦合计算方法得到的滑坡过程(图13(e)—(h))基本一致，可证实本方法的可行性。