堰塞体形成全过程的连续离散耦合数值模拟

王 叶，周 伟，马 刚，陈 远，常晓林

(1. 武汉大学水资源与水电工程科学国家重点实验室，武汉 430072；2. 武汉大学水工岩石力学教育部重点实验室，武汉 430072)

自然界中堰塞体的形成因素很多，根据统计，有90%以上是地震诱发滑坡形成的[1]，滑坡体堆积河谷形成堰塞体，一旦失稳将引发下游洪水灾害，因此深入研究地震作用下的堰塞体形成机制，对于评估堰塞体的稳定性、堰塞湖潜在次生灾害可能性及其损失具有重要意义。

目前对于堰塞体形成机制的研究方法主要有有限元法(FEM)[2,3]、离散元法(DEM)[4-7]、非连续变形分析方法(DDA)[8,9]，连续-离散耦合分析方法(FDEM)[11-17]等。徐文杰等[2]在对肖家桥滑坡区详细的地质调查基础上，结合动力有限元分析进行，再现滑坡三维空间失稳过程；刘传正等[3]通过有限差分方法研究了红石岩山体在地震的动力响应特征和崩塌的形成机理。在强震作用下，岩土体发生局部大变形甚至形成强间断面(岩石拉断、沿结构面滑移)，此时基于连续介质力学的传统有限元方法和有限差分法受小变形假设的限制而很难适用。曹琰波等[4]采用离散单元法，对唐家山滑坡由变形累积到破坏滑动的全过程进行了模拟，研究了地震作用下顺层岩质滑坡全过程；Wu等[5]采用离散元方法分析一个地震诱发的大规模滑坡，研究滑坡体滑动过程和滑坡后的形态等；焦玉勇等[6]采用离散单元法模拟了边坡失稳过程；邬爱清等[8]采用非连续变形分析方法，以唐家山滑坡形成的堰塞坝形态和位置作为目标函数，对唐家山滑坡过程进行复演，探讨地震荷载下高速滑坡的形成机制；Zhang等[9]提出附加接触模型的改进的DDA方法，研究有黏聚力材料的边坡模型，评估边坡节理面处的抗剪强度；Huang等[10,11]采用考虑节理动摩擦退化的改善非连续变形分析方法(IDDA)研究地震诱发的滑坡体的速度、位移等动力学行为。

DEM和DDA能模拟边坡失稳及滑动全过程，但是它们对裂纹、应力、应变的表征不够直观，因此有学者开始采用连续离散耦合分析方法来研究边坡失稳及滑动全过程。Barla等[12]采用FDEM对位于意大利Alpetto矿山的2处边坡进行数值模拟，与早前的研究进行比较，验证了FDEM模拟边坡失稳的可行性；Antolini等[13]采用FDEM方法对Torgiovannetto di Assisi滑坡体进行数值模拟，研究滑坡体的启动机制和演化过程，为滑坡体的风险评估提供参考。Piovano等[14]研究了意大利Aosta山谷的Beauregard滑坡，表明运用FDEM方法能够预测深层滑坡的典型不稳定面和破坏机理；常晓林等[15]采用FDEM模拟了在暴雨等情况下边坡失稳滑移的过程，较好地模拟了抗剪断参数降低时边坡从小变形-大变形-局部滑动-整体滑移的全过程失效；Zhou等[16]基于FDEM方法研究了节理倾角对边坡破坏模式、运动过程、破坏后的堆积形态等的影响。O. K. Mahabadi等[17]应用基于FDEM的Y-Geo代码较好地模拟了山崖侵蚀的过程。但是FDEM进行滑坡全过程模拟大多限于二维，对三维滑坡体研究甚少。

在连续-离散耦合分析方法中引入界面单元，采用内聚力模型能有效地捕捉加载过程中裂纹的扩展路径，模拟岩体渐进破坏过程，充分发挥连续介质力学方法和非连续介质力学方法各自的优势[15-18]，实现岩体从连续状态到非连续状态的转化过程。本文采用FDEM，结合红石岩堰塞体，对三维滑坡由变形破坏到整体滑动的全过程进行模拟，以研究地震作用下岩质滑坡的失稳破坏过程，直观地捕捉滑坡体从开裂破坏到最终的堆积状态、滑坡体失稳过程中的速度分布情况等。

1 基于内聚力模型的岩石开裂扩展模拟

在岩体连续-离散耦合分析方法中[12-18]，岩石被离散为由实体单元和无厚度界面单元组成的系统，每个实体单元均为单独的离散块体。界面单元失效之前，离散块体之间通过界面单元“黏结”成一个整体，以模拟岩石材料的连续变形。实体单元只发生弹性变形，损伤和断裂仅发生在界面单元上。界面单元的破坏准则为带拉断的Mohr-Coulomb准则，界面单元的应力状态满足破坏准则后，采用基于断裂能的线性损伤演化模型，损伤因子达到1.0后完全失效。界面单元失效后从模型中删除，两侧实体单元发生接触关系，采用罚函数法进行接触分析。

采用内聚力模型模拟岩石材料的界面开裂行为，通过界面单元的失效模拟裂纹的萌生、扩展、交汇和贯通。在张拉应力状态下，简化的内聚力区应力分布见图1。真实裂缝的尖端应力为零，内聚力区尖端应力为材料的抗拉强度，从真实裂缝尖端到内聚力区尖端应力逐渐提高。二维界面单元具有2个积分点，均位于单元的中平面上。采用共节点的方式实现界面单元与相邻实体单元之间力、位移的传递。图2给出局部模型的有限元网格拓扑结构图，裂缝可从布置界面单位的位置扩展。

图1 准脆性材料的内聚力区应力分布Fig.1 Illustration of cohesive zone for quasi brittle material

图2 网格拓扑图(实体单元面积收缩便于显示界面单元，实际界面单元厚度为0)Fig.2 Mesh topology (the elastic elements are shrunk for illustration purpose, the actually element thickness is zero)

在界面单元出现损伤之前，假定其本构关系是线弹性：

(1)

式中：tn、ts分别为界面单元的法向应力、切向应力；Kij为刚度矩阵，本文计算时未考虑应力、应变的法向分量与切向分量之间的耦合关系，因此Kns=0 ；εn、εs分别为界面单元的法向应变、切向应变。

法向应变、切向应变定义为：

(2)

式中：δn、δt分别为界面单元沿着法向和切向产生的位移；T0为界面单元的本构厚度。

本文采用无厚度的界面单元模拟材料开裂，倘若取T0为几何厚度，将会导致应变奇异。为了消除应变奇异，在进行界面单元的应变和应力计算时采用的是单元的本构厚度。通常定义T0=1，则界面单元的应变与相应方向上的位移大小相等。界面单元的线弹性本构方程可写成：

(3)

式中：kn、ks分别为界面单元的法向刚度、切向刚度。

为了模拟裂纹的萌生、扩展过程，需要定义合理的裂缝起裂准则，考虑岩石等准脆性材料的破坏往往是由拉应力和剪应力共同作用导致的，本文采用二次应力准则判断界面单元是否开裂，采用带拉断的Mohr-Coulomb准则计算岩石材料的抗剪强度：

(4)

式中：c为岩石材料的黏聚力；φ为岩石材料的内摩擦角；t0n为抗拉强度；t0s为抗剪强度。

二次应力准则：当各个方向的应力与其相应的临界应力的比值的平方和达到1时，界面单元即开始出现损伤：

(5)

在损伤演化阶段，界面单元的本构方程为：

(6)

ts=(1-D)ksδs

式中：D为无量纲的损伤因子，当D=0时，表明界面单元未出现损伤；当D=1时，表明界面单元完全失效，失去承载能力。

损伤因子的表达式如下：

(7)

式中：t0eff为达到破坏准则时的等效应力；δ0m为达到破坏准则时的等效位移；δmaxm为加载历史中的最大等效位移。

在大多数情况下，界面单元处于在混合加载状态下，同时发生法向和切向的变形，因此需要定义界面单元的等效应力teff和等效位移δm：

(8)

通过定义界面单元在损伤过程中耗散的能量，即可控制界面单元的损伤演化过程。在本文中，耗散能量也就是通常所说的断裂能，它等于应力-分离曲线下面所包含的面积大小。本文模拟时采用基于线性软化的Benzeggagh-Kenane准则，界面单元的应力-分离曲线见图3。

图3 界面单元的应力-分离曲线(Ⅰ、Ⅱ)Fig.3 Constitutive relations of the cohesive elements

Benzeggagh-Kenane准则的表达式如下：

GT=Gs+GⅠ

(9)

式中：Gshear为界面单元在剪切荷载作用下的断裂能；GT为界面单元在混合荷载作用下的断裂能；η为材料常数，通常由弯曲试验测得，本文取2。

2 红石岩滑坡体FDEM三维模拟

2.1 工程概况

红石岩堰塞坝形成于2014年发生的云南鲁甸“8·03”地震。“8·03”地震发生后，左岸滑坡堆积物表层松动并向河床滑动，右岸山体产生大规模滑坡，在河床形成堰塞湖。根据现场调查，堰塞体顶部呈马鞍形，顶部左岸高，右岸低，右岸边缘为滑坡岩石堆积体。堰塞体组成松散，最大厚度约103 m。顶部横河向最低高程点1 222 m，堰塞体左岸最高点为1 270 m，上游迎水面平均坡比约1∶2.5，下游面平均坡比约1∶5.5。顺河向底宽约910 m，沿1 222 m高程坝轴线长度约307 m。右岸滑坡后，地形发生了较大变化，上部滑床后缘为陡崖，高度估计约150～200 m，中部形成一个横河向近水平且倾向下游的斜面地形，在此斜面以下为陡崖。左岸坡虽有表面松散并向下滚落，但体量不大，总的地形未有大的改变。堰塞体物质组成较为复杂，主要为右岸边坡崩滑后形成的崩塌堆积物。

2.2 FDEM模拟所需参数率定

采用连续-离散耦合分析方法进行数值模拟时，需要对所用参数进行率定。采用二维平面应变模型进行双轴数值压缩试验来率定模拟参数，通过不同围压下的偏应力应变曲线得到应力莫尔圆从而计算出岩土体的宏观抗剪强度参数。用于参数率定的数值试样尺寸为80 mm×160 mm(见图4)。刚性加载板与试样之间定义接触关系，采用位移控制方式加载施加于上、下刚性加载板，围压施加在左右加载板上。刚性加载板与数值试样之间的摩擦系数取为0.1。

图4 双轴压缩试验试样尺寸及加载示意图Fig.4 Size of numerical experiments of biaxial compression and illustration of loading

红石岩滑坡体各土层的物理力学参数见表1。在FDEM数值模拟中，界面单元的参数包括刚度、抗拉强度、黏聚力c、内摩擦角φ、断裂能。通过反复试算，使得数值试验得到的力学参数与表1所示的室内试验值接近。表2对比了数值试验得到的各土层宏观黏聚力、内摩擦角与室内试验各土层抗剪强度，相对误差都在10%以内。最终确定的FDEM模拟参数见表3。

表1 红石岩岩体室内试验物理力学参数Tab.1 The physical and mechanical parameters of laboratory tests for Hongshiyan rocks

表2 红石岩岩体力学参数与数值试样计算力学参数结果对比Tab.2 The comparison of mechanical parameters ofrocks for Hongshiyan rocks between laboratorytests and numerical experiments

表3 红石岩岩体数值模拟界面单元细观参数Tab.3 The mesoscopic parameters of cohesive interface elements for numerical simulation of Hongshiyan rocks

2.3 计算模型

根据震前边坡原地形图进行滑坡前的还原。由于地震影响范围较广，为减小边界反射效应，将模型边界范围取得足够大。三维计算模型中，x轴方向为横河向，由左岸指向右岸为正；y轴方向为顺河向，向下游为正；z轴方向为竖直向，向上为正。横河向长1 427.1 m，纵河向1 620.5 m，竖直向970 m。

地震荷载作用下，只考虑边坡底部的弹性作用，消除边坡底部对地震的放大作用，采用有刚度无质量方案进行分析，故底部约束为固定约束，侧向边界为法向约束。在尽量考虑实际地形地貌和地质结构的前提下，同时考虑计算效率，三维整体模型共剖分为341 651 个单元和376 002 个节点，其中实体单元为215 271个，界面单元为126 380 个。本文重点关注滑坡形成堰塞体的整个过程，所以只在滑坡体部分插入界面单元。整个FDEM模拟分2步进行，首先进行地应力平衡，得到初始应力场，再施加地震加速度荷载，进行地震分析。

在计算模型底部输入3个方向的地震加速度，模拟滑坡体在地震作用下的启动、滑动到最终堆积直至形成堰塞体的全过程。加速度时程曲线见图5(a)、5(b)、5(c)，地震作用过程历时20 s(10～30 s)，总的作用时间为100 s。

图5 地震加速度动力时程曲线Fig.5 Time-history curves of earthquake acceleration

3 堰塞体形成全过程分析

图6为堰塞体最大剖面对比图(对比范围相同)，实际堰塞体最大堆积厚度约为103 m，模拟的对应纵剖面最大堆积厚度为101.97 m，由于模型的简化和模拟范围限制，2者的堆积形态稍有不同，但大致相似。

图6 堰塞体最大纵剖面对比Fig.6 Comparison for longitudinal section of landside dam

定义堰塞体堆积最大宽度为Lmax，堆积最大高度为Hmax，宽高比Lmax/Hmax，用这3个参数表征堆积形态。表4统计了4个典型剖面堆积形态实测值和FDEM计算值及2者的相对误差(百分制)。由于模型简化和网格划分关系，表4中个别计算值与实测值有较大差别，故其相对误差较大。图7对比了4个典型横剖面实测与FDEM模拟的堆积形态。从表4和图7可以看出模拟得到的横剖面堆积形态与实测堆积形态基本一致，进一步验证了本文FDEM模拟边坡失稳、滑动形成堰塞体的合理性。

表4 典型横剖面堆积形态实测值与计算值对比Tab.4 The comparison of accumulation state of typical sections between measured value and calculated value

滑坡体不同时刻的滑动状态见图8。选取4个时刻滑坡典型剖面g-g′的形态[见图9(a)～图9(d)] 研究地震诱发滑坡失稳过程运动特征。从图8和图9可以看出，在地震初期由于受地震作用，表层岩体首先出现裂缝然后贯通形成整体滑裂面；随着地震作用，滑坡体内部损伤不断累积，t=20 s左右时，整个滑坡体与下覆基岩出现明显的运动分离[见图9(a)]；然后在地震和重力共同作用下滑坡体高速下滑，且在下滑过程中滑坡体不断破碎、解体，受对岸山体的阻挡，与之发生高速撞击，进一步破碎、解体、堆积、堵塞河道，在60 s后滑坡体逐渐稳定，逐渐堆积密实，100 s时基本静止，形成堰塞体。

图7 典型横剖面实测与计算堆积形态对比Fig.7 The comparison of accumulation state of typical sections between actual measurement and calculation

图8 滑坡体不同时刻滑动形态Fig.8 Sliding patterns of landslide at different time

图9 滑坡体典型剖面g-g′不同时刻滑坡形态Fig.9 Sliding patterns for a typical section named “g-g′” at different time

图10为滑坡体典型剖面g-g′不同时刻速度云图。t=20 s左右滑坡体速度开始增大，滑块底部首先达到8 m/s左右，在地震和重力作用下，滑坡体的速度持续增大，t=30 s左右高速下滑，达到20 m/s左右，速度的分布受边坡坡度的影响，由于上部坡度较陡而中部坡度缓，下部与中部坡度相近，但下部是临空面，而中部有下部岩体阻挡，故岩体失稳时滑动速度中部较小，上部和下部速度较大；t=40 s滑坡体已经解体，滑坡体在下滑过程中积聚很大的动能，故此时分布有最大滑坡速度，撞击对岸岩体后速度迅速改变方向，不断堆积河谷，直至逐渐稳定，60 s左右速度逐渐减小为零。

图10 滑坡体典型剖面g-g′不同时刻速度云图Fig.10 Cloud pictures of velocity of a typical section named “g-g′” at different time

图11显示了系统的能量演化曲线(E为能量值，Emax为相应能量的最大值)，包括系统总动能历时曲线、系统摩擦耗散能历时曲线和破碎耗散能历时曲线。在历时20 s左右，系统的总动能和摩擦耗散能出现突变，边坡开始失稳，滑坡体高速下滑，在40 s左右由于撞击对岸岩体，动能开始迅速减小。在t=60 s，各能量基本保持不变，滑坡堆积体达到稳定状态。从摩擦耗散能和破碎耗散能历时曲线可知，随着滑坡体的失稳和滑动，失效界面单元不断增加，导致滑坡体由连续体向非连续体转换，并最终形成松散的堆积体。

图11 系统能量演化曲线Fig.11 Evolution curves of system energy

4 结 论

(1)将FDEM应用于红石岩的三维边坡失稳模拟，重现了堰塞体形成的全过程。在地震作用10 s左右，系统动能出现突变，边坡开始失稳，地震总的作用时间20 s后，滑坡体高速下滑，进而受到对岸撞击，不断解体并堆积在河谷。通过FDEM分析能够直观地捕捉滑坡体各个时刻的滑动状态、速度分布和滑坡体的最终堆积形态。

(2)将最大纵剖面和4个典型横剖面的数值模拟的堆积形态与现场实际堆积形态进行对比，2者吻合程度较高。对比了表征堆积形态的参数，基本一致，进一步验证了FDEM模拟边坡失稳和滑动堆积形成堰塞体全过程的合理性。

(3)FDEM三维模拟结果表明，红石岩堰塞体形成机制为地震诱发岩土体内部损伤和变形累积、结构面贯通、滑坡体高速下滑并伴随破坏、解体、撞击对面山体，堆积河道、不断密实，形成堰塞体。

□

[1] Costa J E, Schuster R L. The formation and failure of natural dams [J]. Geological Society of America Bulletin, 1988,100(7):1 054-1 068.

[2] 徐文杰，陈祖煜，何秉顺，等. 肖家桥滑坡堵江机制及灾害链效应研究[J]. 岩石力学与工程学报，2010，29 (5)：933-942.

[3] 刘传正，葛永刚，江兴元，等. 鲁甸地震红石岩崩塌触发机理分析[J]. 防灾减灾工程学报，2016，36(4)：601-607.

[4] 曹琰波，戴福初，许 冲，等. 唐家山滑坡变形运动机制的离散元模拟[J]. 岩石力学与工程学报，2011，30(增1)：2 879-2 887.

[5] Wu Jianhong, Lin Jeenshang, Chen Chaoshi. Dynamic discrete analysis of an earthquake-induced large-scale landslide[J]. International Journal of Rock Mechanics & Mining Sciences, 2009,(46):397-407.

[6] 焦玉勇，葛修润，刘泉声,等. 三维离散单元法及其在滑坡分析中的应用[J]. 岩土工程学报，2000，20(1)：101-104.

[7] Zhou Wei, Lai Zhiqiang, Ma Gang. Effect of base roughness on size segregation in dry granular flows[J].Granular Matter, 2016,(18):83.

[8] 邬爱清，林绍忠，马贵生，等. 唐家山堰塞坝形成机制DDA模拟研究[J]. 人民长江，2008，39(22)：91-95.

[9] Zhang Yingbin, Xu Qiang, Chen Guangqi, et al. Extension of discontinuous deformation analysis and application in cohesive-frictional slope analysis[J]. International Journal of Rock Mechanics & Mining Sciences, 2014,(70):533-545.

[10] Huang Da, Song Yixiang, Cen Duofeng, et al. Numerical modeling of earthquake-induced landslide using an improved discontinuous deformation analysis considering dynamic friction degradation of joints[J]. Rock Mech Rock Eng, 2016,49:4 767-4 786.

[11] Barrett P J. The shape of rock particles, a critical review [J]. Sedimentology, 1980,27(3):291-303.

[12] Marco Barla, Giovanna Piovano, Giovanni Grasselli. Rock slide simulation with the combined finite-discrete element method [J]. Geomechanics, 2012,12(6):711-721.

[13] Francesco Antolini, Marco Barla, Andrea Giorgetti, et al. Combined finite-discrete numerical modeling of runout of the torgio-vannetto di assisi rockslide in central italy [J]. International Journal of Geo-mechanics, 2016,16(6):1-16.

[14] Giovanna Piovano, Francesco Antolini, Marco Barla, et al. Continuum-discontinuum modelling of failure and evolution mechanisms of deep seated land-slides[C]∥6th Int. Conf. on Discrete Element Method, Colorado School of Mines, Golden, CO, USA, 2013:295-300.

[15] 常晓林，胡 超，马 刚，等. 模拟岩体失效全过程的连续-非连续变形体离散元方法及应用[J]. 岩石力学与工程学报，2011，30(10)：2 004-2 011.

[16] Zhou Wei, Yuan Wei, Ma Gang, et al. Combined finite-discrete element method modeling of rockslides [J]. Engineering Computations, 2016,33(5):1 530-1 559.

[17] O K Mahabadi, A Lisjak, A Munjiza, et al. Y-Geo: new combined finite-discrete element numerical code for geome-chanical applications[J]. Geomechanics, 2012,12(6):676-688.

[18] 马 刚，周创兵，常晓林，等. 岩石破坏全过程的连续-离散耦合分析方法[J]. 岩石力学与工程学报，2011，30 (12)：2 444-2 454.