基于连续-非连续耦合方法的降雨滑坡数值模拟研究

陈闻潇， 石 崇， 李汪洋， 张一平

（1.河海大学岩石力学与堤坝工程教育部重点实验室，南京 210098； 2.河海大学岩土工程科学研究所，南京 210098）

滑坡过程分析是岩土工程评价与灾害评估的重要内容［1］，工程师通过分析滑坡过程中出现的裂隙位置、分布、孕育发展过程判断滑坡破坏模式，进一步评估潜在滑坡灾害的规模和范围，探讨合理的防灾减灾措施［2-5］.

滑坡数值仿真是进行滑坡过程分析的重要手段. 对于降雨滑坡的模拟分析，李宁等［6］提出了改进的Mein-Larson 降雨入渗模型；蒋中明等［7］建立了基于FLAC有限差分的非饱和渗流计算方法；谢强等［8］分析了降雨影响下边坡的安全系数变化. 但是这些学者提出的降雨计算模型都是基于连续介质力学，不能分析降雨滑坡造成的大变形破坏、滑坡运动过程、堆积形态与灾害影响范围等.

离散元方法被广泛运用到计算边坡的滑坡过程研究中［9-10］，与连续力学方法相比，其具有能模拟大变形的优势. 在离散元方法中，考虑降雨影响的主要方法是参数折减或计算等效入渗深度［11］，或采用流固耦合［12］的方法，但是这些方法中假定土体饱和，不满足土的饱和-非饱和特性. 离散元中采用的简化方法势必会造成结果的偏差. 另外，滑坡往往仅限于一定范围，如果整个模型都采用离散单元方法，很容易造成颗粒或块体数目过多，导致计算工作量超出电脑承受能力.

因此部分学者建议采用连续-非连续耦合计算方法分析滑坡等大变形问题［13-14］，该方法将可能发生破坏的土体用离散元求解，而变形量小的周边岩土体用连续介质方法计算，能极大地提高计算速度与精度［15］. 由于离散元进行饱和-非饱和降雨入渗模拟相对困难，然而在连续数值模拟方法中很容易做到，如果将两者结合，进行连续-非连续耦合计算，那么不仅能模拟降雨入渗，同时可以模拟降雨诱导的滑坡大变形破坏.

本文基于连续-非连续耦合分析方法，利用边界墙耦合传递方法，提出了降雨条件下滑坡的模拟分析方法，借助案例分析了该方法的可行性，并探讨了耦合滑坡数值仿真过程中的控制因素及影响规律，研究结果可为连续-非连续数值仿真在岩土工程中的应用提供依据.

1 连续-非连续耦合分析实现方法

目前连续-非连续数值模拟实现方法主要有两种：基于边界控制颗粒［16］和基于边界控制墙体［17］的耦合方法. 本文采用的连续-非连续耦合分析算法以有限差分算法（FLAC）和颗粒离散元方法（PFC）原理为基础，通过基于边界控制墙体的方法来进行信息交互与耦合计算，属于边界控制墙体方法. 边界墙由三角形面组成，其顶点速度和位置为时间的函数. 耦合逻辑的工作原理是：在每一个耦合分析步中，获取非连续模型与边界墙接触的力和力矩，通过等效力系统，将其施加于连续模型的节点上，在连续模型中计算后，将节点的速度和位置传递给非连续颗粒，通过非连续模型中力-位移法则，计算得出新时间步的力和力矩，如此就完成了一个耦合循环.

在耦合模型中当离散颗粒与三角形的边界墙相接触时，在墙面上的接触点为CP. 如图1所示，三角形边界墙顶点位置为xi（i=1，2，3），将三角形的顶点与点CP连接，得到三个区域的面积为Ai（i=1，2，3），三角形的总面积为A. 定义每个顶点的权重因子wi（i=1，2，3）为一个顶点对边的三角形面积除以三角形的总面积，即wi=Ai/A. 定义ri（i=1，2，3）为从点CP指向各自三角形顶点的向量，有ri=xi-CP.

非连续模型每一时步将力和力矩传递至连续模型，其实现需要通过等效力方法. 作用于接触点CP处的接触力为F，黏接产生的接触力矩为M. 通过等效力方法计算后的作用力位于边界墙的每个顶点xi（i=1，2，3），其力的大小为Fi（i=1，2，3），通过下式计算：

图1 连续-非连续耦合边界墙Fig.1 Continuous-discontinuous coupling boundary wall

定义n为指向三角形法线方向的单位向量，则沿三角形边界墙表面的切向力矢量为：

切向单位矢量为：

力的求解通过建立局部坐标系，使得局部坐标系中x分量与n方向一致，y分量与s方向一致，此时每个ri的x分量为0，即ri,x=0 . 此外，将重心加权wi（i=1，2，3）施加在三角形平面上最大接触力的方向. 这种简化可以直接求解局部坐标系统中顶点的力和力矩.

连续模型每一时步将节点的速度传递至非连续模型. 节点传递至边界墙顶点xi（i=1，2，3）的速度分别为Vi（i=1，2，3），假定速度场在三角形中线性变化，则可以通过重心插值方式得到接触点CP处的速度值，并将该点的速度值作为非连续颗粒的速度值，其计算公式如下：

2 降雨荷载条件考虑与计算

2.1 降雨入渗的饱和-非饱和计算

本文提出的降雨滑坡的连续-非连续耦合分析流程如图2所示，模型的模拟分为降雨过程的模拟和滑坡过程的模拟. 降雨模拟通过连续力学方法计算，由于降雨是长时间尺度的过程，其持续时间通常以小时和天为单位，在分析时通常不考虑变形［7，18］，故固定变形以节省效率，采用基于饱和-非饱和渗流理论对降雨进行计算，并对孔压场和渗流场进行记录. 滑坡过程的模拟将允许模型自由变形，将可能发生破坏的区域进行离散化，基于耦合边界墙的方法生成连续-非连续模型，并将孔压场和渗流场通过等效力和土体弱化的方法叠加进模型，进而对降雨滑坡进行连续-非连续耦合计算分析.

图2 耦合方法流程图Fig.2 Flowchart of coupling method

降雨过程模拟采用饱和-非饱和渗流理论. 非饱和土渗流分析中考虑土-水特征线和水力传导方程，其中土-水特征线用来描述含水率与吸力的关系，水力传导方程用来描述渗透系数与吸力的关系. 本文采用常用的Van Genuchten（VG）模型［19］来模拟非饱和土的水力特性. VG模型中土-水特征曲线的表达式为：

式中：θ为土的体积含水率；θr为残余体积含水率；θs为饱和体积含水率；P 为土体的孔隙水压力；ɑ、n′、m′为拟合参数. 由于体积含水率与饱和度s 关系为θ=n·s，通过上式可得到土体饱和度与负孔隙水压力的关系式为：

式中：Sr为残余饱和度，非饱和土的渗透系数k( s) 与饱和渗透系数K的关系为：

为了实现降雨边坡的模拟，需要对边坡表面的降雨边界条件进行实时调整修正，即动态边界条件. 当降雨强度小于土体的入渗率时，边坡表面不会产生积水，入渗量的值取降雨强度值. 当降雨强度大于土体的入渗率时，边坡表面会产生积水，入渗量的值取土体的最大入渗率. 在计算的每一个时间步都会通过自定义FISH函数，将降雨强度与土体的入渗率的大小进行对比，从而实现动态边界条件的模拟.

2.2 降雨影响的连续-非连续等效力方法

在通过饱和-非饱和渗流计算，得到边坡的渗流场与孔压场之后，需要通过等效简化的方式，将降雨影响施加到非连续力学的颗粒中. 首先必须考虑由于雨水入渗而引起的土体重度的增加，受降雨影响的土体重度ρ 可按下式计算：

式中：ρd为土的干重度；θ为体积含水率；ρw为水的重度. 在饱和土体中由于水头压力差的存在会对土体产生渗流力，渗流力的方向与渗流方向一致［20］，且忽略非饱和土体内渗流力的影响. 在饱和区域内施加的等效渗透作用力Ff为：

式中：J为渗透力；V为单元网格或颗粒单元的体积；γW为水的重度（N/m3），i为水力梯度，其大小值为土体中两点水势之差与其渗透距离的比值. 颗粒在饱和作用下受到的浮力的大小可以由以下公式计算得出：

此外，土体在饱和状态下，由于软化作用使其物理力学参数存在一定的折减，本文根据工程经验将饱和区内的土体参数折减15%. 在非连续模型中，降雨荷载的影响考虑由公式（8）～（10）计算的土重增加、渗流力和浮力及参数折减. 在连续模型中，由于不是滑坡发生的主要区域，仅考虑公式（8）计算的土重增加和参数折减，以保证合理的计算效率.

3 降雨滑坡的耦合计算分析

本文通过一个案例来分析该连续-非连续耦合降雨滑坡模拟方法的可行性，为避免关注点在于诱发滑坡的降雨阈值问题，本文选取的边坡安全系数为1.1，在给定降雨影响后条分法安全系数小于1.0，即边坡将在降雨影响下产生滑动. 如图3所示建立连续边坡模型，在此基础上开展降雨条件下连续-非连续耦合滑坡分析.为了兼顾问题的三维性，使用假三维模型来进行平面问题的模拟. 模型共有14 636个网格单元和29 796个网格节点，假三维厚度为3 m，边坡土体饱和渗透系数为1.0×10-4cm/s，公式（5）～（7）所需的VG 模型参数采用文献［7］中给出值，以反映渗透系数、饱和度与基质吸力之间的关系.土体内非饱和区的渗透系数由式（5）～（7）计算.如图3所示，初始孔压场通过在模型两侧边界设定孔压边界条件，然后采用饱和-非饱和计算稳定生成，初始地下水位线（黑线）为孔隙水压力为0的等值线，孔隙水压力为负数的区域为非饱和区.

图3 初始计算模型及孔压分布（单位：kPa）Fig.3 Initial calculation model and pore pressure distribution

给定降雨强度为1.0×10-6m/s，连续降雨10 d，通过上文描述的动态边界条件方法，将降雨荷载施加于模型的上部边界. 图4 为边坡在降雨10 d 后的孔隙水压力分布及饱和度等值图，可见此时边坡的坡脚，坡面以及坡顶已经由负孔隙水压力升至正孔压，这表明此区域内的土体已经由非饱和转化为暂态饱和区，且在坡脚处与坡体下部的正孔压区连通，暂态饱和区内孔隙水压力最大不超过50 kPa. 由于边坡坡面及坡顶的水分在重力作用下向坡脚处下渗，所以坡脚处的饱和区范围更大，地下水位上升幅度也较为明显.

图4 降雨计算结果Fig.4 Rainfall calculation results

在降雨影响结果计算完成后，为了继续对降雨导致滑坡过程进行模拟，需通过边界墙耦合方法，建立连续-非连续耦合模型. 将降雨计算的结果保存至单元内另开辟的存储节点内，目的是在耦合模型生成后，通过上文的等效力方法计算降雨荷载并分别施加于连续模型和非连续模型中. 然后通过连续-非连续耦合方法，将滑坡区域的潜在危险部分用离散介质替换. 本案例中认为整个边坡均属于危险区域，在工程应用中，可结合实际滑坡体及其邻近区域设为危险区域［2］. 为了保证潜在滑坡体与连续介质的连续性，防止平滑的薄弱交界面出现，连续-非连续的接触面采用锯齿状分布.

耦合模型如图5 所示，共生成18 153 个半径处于0.15～0.30 m 之间的颗粒. 通过宏观-细观力学参数标定，得到模型宏-细观力学参数如表1、表2所示，局部阻尼值为0.3. 由图5（c）、（d）、（e）可见离散颗粒和连续单元之间通过边界墙来形成连接，两者紧密接触，这是典型的基于边界控制墙的离散-连续耦合方法.图5（a）为连续非连续初始模型的总应力图，可见在将连续介质替换成离散介质后，传递至下方的力满足自然应力分布，模型的总应力符合条件，但由于离散颗粒具有一定的随机性，所以在容许情况下应力存在细微的波动. 图5（b）为模型初始平衡时产生的位移云图，可见连续-非连续模型之间的位移是连续的，这也证明了耦合的可行性. 但是与单一介质模型平衡时产生的位移不同，耦合时因为初始状态生成时，离散颗粒与单元网格之间的接触并不紧密，接触力也并不均匀，所以最大位移产生在竖直耦合交界面的附近，直到模型逐渐平衡，传递均匀的接触力，位移才不再发展，形成初始模型，因此在之后的模型计算中，需先对此位移清零.

表1 计算模型细观力学参数Tab.1 The mesomechanical parameters of calculation model

表2 计算模型宏观力学参数Tab.2 The macromechanical parameters of calculation model

图5 连续-非连续耦合模型Fig.5 Continuous-discontinuous coupling model

图6 滑坡过程Fig.6 Landslide process

降雨滑坡耦合模型计算结果如图6所示. 由图6（a）、（b）可见，在降雨入渗形成饱和区的影响下，边坡滑动且形成了明显的圆弧滑动面，滑坡类型主要为浅层滑坡，滑动区域以饱和区为主，圆弧最大深度处超过了暂态饱和区分界线，这说明了非饱和区的土发生了滑坡. 在降雨影响下，滑坡体迅速沿着坡面滑落，并由于局部阻尼的影响，滑坡体最终在坡脚堆积并停止进一步运动. 滑坡的最大位移为29.3 m，滑坡结束后堆积角约为31.6°. 由图6（c）、（d）可见在滑坡启动阶段，各测点的速度迅速增大，而测点3的速度在7 s左右才明显增大，这表明了该滑坡为牵引式滑坡，前缘由于在降雨作用下更容易饱和，饱和区范围较大，土体强度弱化而率先滑动，而处于后缘的土体破坏较迟，其破坏原因受降雨影响和前缘破坏失稳的双重影响. 滑坡体速度在10～25 s左右达到峰值，随后速度逐渐减小，在175 s左右滑坡结束. 图6（e）为滑坡过程中的裂隙数变化曲线，由于降雨作用的影响，在滑坡的初期就产生了一定的裂隙，并在滑坡过程中，出现了几个裂隙陡增的时段，这说明可能发生了小的崩塌. 计算结果表明，采用连续-非连续耦合计算方式，对降雨-滑坡进行耦合分析，能够较好地进行降雨作用下滑坡过程和破坏特征分析.

4 降雨滑坡中的阻尼影响

在滑坡分析中，通常通过设置阻尼来控制能量衰减，防止动能累积过快［21］. 故在降雨滑坡耦合分析过程中阻尼的参数尤为重要，在不同阻尼条件下的滑坡过程分析图7所示. 由图7（a）可见，在滑坡进行的0～30 s内，饱和区内土体的速度和阻尼呈负相关的趋势，阻尼越小，平均速度峰值越大，达到峰值所需的时间也越短. 在滑坡发生30 s之后，平均速度在不同阻尼影响下呈现相同的趋势，即随着滑坡进行不断减小，滑坡体逐渐稳定. 总体上，在不同阻尼下，滑坡趋势均保持一致，即速度先迅速增大至峰值，然后逐渐减小至稳定，位移表现为先迅速增大而后逐步稳定. 由图7（b）可见，暂态饱和区的平均位移受局部阻尼的影响，具体体现在距坡脚的堆积距离上，阻尼越大堆积距离越小，平均位移也就越小.

从图7（c）、（d）可以看出，当阻尼小于0.3时，滑坡体积与堆积距离明显增加，这是因为阻尼较小，能量耗散缓慢，在滑坡中产生的冲击碰撞以及堆积挤压等具有更大的动能，将会使得滑坡造成更进一步的破坏，故滑坡影响区域也越大. 而当阻尼大于0.3时，滑动区域变化不明显，堆积距离也基本一致. 合理的获取阻尼参数，在降雨滑坡耦合分析过程中尤为重要，通过将不同阻尼的计算结果，与真实滑坡中的滑坡征兆、运动过程和堆积状态对比，来确定合理的阻尼参数，进而实现真实情况下的降雨滑坡耦合模拟，并进行更进一步的分析.

图7 局部阻尼影响因素分析Fig.7 Analysis of influencing factors of local damping

5 结论

本文基于连续-非连续数值模拟耦合分析方法，采用边界墙耦合理论和饱和-非饱和降雨入渗分析，通过案例对降雨滑坡过程进行了模拟，分析了该耦合计算方法的可行性，并研究了局部阻尼在降雨滑坡耦合分析中的影响规律，得到主要结论如下：

1）利用连续模型对饱和-非饱和降雨入渗分析，通过边界墙耦合方法生成连续-非连续耦合模型，并将渗流场通过等效力施加于耦合模型中，实现了降雨影响下的滑坡耦合分析，从而建立了基于连续-非连续耦合分析的降雨滑坡分析方法.

2）通过案例分析验证了降雨滑坡耦合分析的可行性. 计算结果表明该方法可以较好地模拟降雨滑坡，并进行滑坡过程和破坏特征的分析，为降雨滑坡的分析提供了新思路.

3）分析了降雨滑坡耦合分析中，局部阻尼影响因素对计算结果的影响. 结果表明局部阻尼对滑坡结果的影响较大，通过实际滑坡的滑坡征兆、运动过程和堆积状态，合理地获取阻尼参数，在降雨滑坡耦合分析过程中尤为重要.