0.   引言

边坡失稳在大多数情况下不是瞬间发生整体破坏，而是由局部滑移破坏逐步发展至整体贯通破坏的过程^[1]。由于成因不同导致滑坡的表现形式各异，人工开挖或库水位骤降一般会引起由前缘启动向后缘发展的滑坡，称为牵引式滑坡^[2-4]；坡顶堆载或降雨可能引发由后缘启动向前缘发展的滑坡，称为推移式滑坡^[5-6]。相比于牵引式滑坡，推移式滑坡在破坏前坡顶会有明显的预警变形，“塑性”破坏的特点更加明显^[7]。

传统的极限平衡法（LEM）只能分析边坡的安全系数而无法计算边坡整体的应力变形情况，且需要假设滑裂面的形状和位置；弹塑性有限元方法（FEM）在变形过大的情况下引起网格畸变，也难以准确刻画边坡的滑裂面发展规律。

近年来，很多学者提出了多种更贴合实际的分析方法。王庚荪等^[8]提出了一种含剪切带的单元模型并将其应用于边坡破坏的过程分析，是国内较早对边坡渐近破坏过程的模拟研究；唐芬等^[9]针对强度折减法采用单一安全系数对强度参数进行折减的不足，提出了对c，φ采用不同的折减系数的双安全系数方法，但是仍未能在行业内达到普遍共识；Conte等^[10]在有限元法中采用非局部弹黏塑性本构模型用以克服不收敛和网格依赖性等问题，进而模拟应变软化边坡剪切带的发展和边坡破坏过程；Mohammadi等^[11]采用大变形有限元法模拟了边坡破坏过程中的变形情况；卢应发等^[12-13]建立了一种全新的条块力-位移分析方法，并提出几种有针对性的边坡稳定性系数计算方法；沈华章等^[14]、薛海斌等^[15]提出了一种基于矢量和法并考虑应变软化的边坡破坏分析方法。然而上述方法受到传统有限元框架所限，难免引入各种假定，而且模拟效果仍不够直观。

扩展有限元法（XFEM）是近些年发展出的不需要网格重构而可以模拟裂缝、剪切带等不连续场的有效方法，正得到广泛的应用^[16-19]。滑坡问题中包含不连续位移场，用XFEM进行研究分析与传统方法相比不需要假定滑裂面的形状和位置，程序可自行搜索和模拟滑裂面的发生发展过程，而且克服了网格畸变问题。

本文基于XFEM构造模拟土质边坡滑裂面发生发展过程的方法，并讨论材料参数对均质土坡推移式滑动的滑裂面深度和滑裂面形态的影响，在此基础上研究成层土坡推移式滑动的滑裂面扩展规律。

1.   XFEM方法介绍

1.1   不连续位移场的构造

XFEM通过对有限元计算域内部分结点富集额外自由度来构造非连续场。对于非连续位移场，利用常规自由度插值构造的连续位移场${\boldsymbol{u}^{{\text{cont}}}}$叠合富集自由度构造的不连续位移场${\boldsymbol{u}^\text{disc}}$来描述：

式中ψ为域内所有结点的集合，ξ为域内富集结点的集合，$\xi \subset \psi$；$\boldsymbol{N}(\boldsymbol{x})$为形函数；u为常规结点位移向量；q为富集结点位移向量；H为Heaviside阶跃函数，且

其中，${\varOmega _{\text{A}}}$和${\varOmega _{\text{B}}}$分别为单元内不连续面切割而成的两个子域，如图 1所示。虽然在整个单元上由于存在不连续面，其位移场和其他场函数也都是不连续的，但是在两个子域上这些场函数是连续的，因而可以将单元积分在两个子域内分别进行^[20]。

图  1  被不连续界面贯穿单元及其积分示意图

Figure  1.  Schematic diagram of an element cut by discontinuous interface and its integration

下载: 全尺寸图片 幻灯片

1.2   控制方程

边坡失稳产生滑裂面的过程可视为一种具有内部演化界面的初边值问题，不失一般性，如图 2所示。研究域为$\varOmega$，b为域内体力。${\varGamma _t}$为力边界，其上作用$\boldsymbol{\bar t}$，${\varGamma _u}$为位移边界，满足$\boldsymbol {u = \bar u}$，且${\varGamma _t} \cap {\varGamma _u} = \emptyset$。${\varGamma _{\rm{d}} }$表示裂缝，为内部边界，分为$\varGamma _{\rm{d}} ^ +$和$\varGamma _{\rm{d}} ^ -$两个界面，即${\varGamma _{\rm{d}} } = \varGamma _{\rm{d}} ^ + \cup \varGamma _{\rm{d}} ^ -$，${\boldsymbol{n}^ + }$，${\boldsymbol{w}^ + }$，${\boldsymbol{t}^ + }$和${\boldsymbol{n}^ - }$，${\boldsymbol{w}^ - }$，${\boldsymbol{t}^ - }$分别对应界面上的法向量、位移向量和应力向量。

图  2  边界非线性初边值问题示意图

Figure  2.  Schematic diagram of nonlinear initial boundary value problem

下载: 全尺寸图片 幻灯片

在保守力系作用下，由虚功原理可得到准静态动量方程的等效积分弱形式：

式中，σ和ε分别为应力张量和应变张量。

裂缝面上的本构方程为

式中，${\boldsymbol{D}_\varGamma }$为局部坐标下的裂缝本构矩阵，Q为局部坐标与全局坐标的转换矩阵。

根据虚位移原理，分别对式（3）中u和q取变分，可整理得到控制方程的有限元格式：

式中，${\boldsymbol{K}}_{uu}={\displaystyle {\int }_{\varOmega }{\boldsymbol{B}}^{\rm{T}}\boldsymbol{DB}\text{d}\varOmega }，{\boldsymbol{K}}_{qu}\boldsymbol{=}{\boldsymbol{K}}_{uq}{}^{\rm{T}}={\displaystyle {\int }_{\varOmega }{\boldsymbol{B}}^{\rm{T}}\boldsymbol{DB}(\boldsymbol{H}-}$${\boldsymbol{H}}_{I})\text{d}\varOmega ，{\boldsymbol{K}}_{qq}={\displaystyle {\int }_{\varOmega }{(\boldsymbol{H}-{\boldsymbol{H}}_{I})}^{\rm{T}}{\boldsymbol{B}}^{\rm{T}}\boldsymbol{DB}(\boldsymbol{H}-{\boldsymbol{H}}_{I})\text{d}\varOmega }，{\boldsymbol{K}}_{qq}^{{\varGamma }_{d}}=$$\int_{\varGamma_{\mathrm{d}}} \boldsymbol{N}^{\mathrm{T}} \boldsymbol{Q}^{\mathrm{T}} \boldsymbol{D}_{\Gamma} \boldsymbol{Q} \boldsymbol{N} \mathrm{d} \varGamma_{\mathrm{d}}, \quad \boldsymbol{F}_u=\int_{\varOmega} \boldsymbol{N}^{\mathrm{T}} \boldsymbol{b} \mathrm{d} \varOmega+\int_{\varGamma_t} \boldsymbol{N}^{\mathrm{T}} \overline{\boldsymbol{t}} \mathrm{d} \varGamma, \quad \boldsymbol{F}_q=$${\displaystyle {\int }_{\varOmega }{(\boldsymbol{H}-{\boldsymbol{H}}_{I})}^{\rm{T}}{\boldsymbol{N}}^{\rm{T}}\boldsymbol{b}}\mathrm{d}\varOmega +{\displaystyle {\int }_{{\varGamma }_{t}}{\left(\boldsymbol{H}-{\boldsymbol{H}}_{I}\right)}^{\rm{T}}{\boldsymbol{N}}^{\rm{T}}\boldsymbol{\bar{t} }}\mathrm{d}\varGamma 。$

1.3   界面接触算法

在形成滑裂面后，滑裂面两侧位移不连续，可以相互错动或者分离，即压剪模式或张剪模式；滑裂面两侧的应力联系通过接触算法实现，本研究在裂缝界面处于张剪状态时采用黏聚裂纹模型^[21]，在处于压剪状态时采用Willner的摩擦接触理论^[22]。相对位移可以分为沿着滑裂面法向的正相对位移${w_{\text{n}}}$和沿着滑裂面切向的切相对位移${w_{\text{s}}}$。

图 3给出了接触算法流程。在每一迭代步中，首先判断相对位移模式：若${w_{\text{n}}}$＞0，表示接触点上的裂缝面相互嵌入，采用压剪模式对其进行计算；若${w_{\text{n}}}$＜0，表示接触点上的裂缝面相互远离，采用张剪模式对其进行计算。接触算法在集成总刚时，会给出${\boldsymbol{D}_\varGamma }$矩阵，在应力积分时，会给出接触应力${t_\text{n}}$和${t_\text{s}}$，进而参与残差平衡计算。

图  3  接触算法流程图

Figure  3.  Flow chart of contact algorithm

下载: 全尺寸图片 幻灯片

1.4   土体的滑裂判别

土体滑裂面尖端和裂缝尖端类似，会产生应力集中和应力主轴偏转。然而土体材料由于其碎散性和弹塑性，一般应力集中程度不高，因此可以利用尖端的一个局部区域中的特征应力状态来判断滑裂面的扩展条件，称这个局部区域为扩展控制域。基于此设置扇形扩展控制域如图 4所示，C扇面位于缝尖的正前方以判定张拉裂缝为主，圆心角为90°，A，B扇面对称于裂缝两侧以判断剪切裂缝为主，圆心角均为135°。每个扇形控制域的特征应力状态由域内所有积分点的应力加权平均得到，其中权函数为

图  4  扇形扩展控制域示意图

Figure  4.  Schematic diagram of fan-shaped extended control domain

下载: 全尺寸图片 幻灯片

式中，r为控制区域内某积分点到缝尖O点的距离，R为扩展控制域的半径，建议取为1.5倍~3倍的网格特征长度。

对于张拉破坏采用最大拉应力准则；对于剪切破坏采用Mohr-Coulomb准则。分别计算尖端A、B、C扩展控制域的拉应力水平T_L和剪应力水平S_L：

式中，σ₁为大主应力，σ₃为小主应力，f_t为土体的抗拉强度，φ为土体的内摩擦角，c为土体的黏聚力。当任意扩展控制域内的任意一种应力水平达到1，则当前控制域内发生相应类型的破坏，这个控制域称为“优势控制域”。当${\sigma _3}$＜0且T_L=1时，出现张拉破坏面，扩展方向为垂直于小主应力的方向；当S_L=1时，出现与小主应力方向成φ/2+45°角的剪切破坏面，具体扩展方向需要根据已有滑裂面进一步判断^[16]。此处引入两个设定条件：①惯性条件，即宏观裂缝向前扩展时总是趋向于沿着与已有裂缝最接近的方向进行；②区域占优假定，即裂缝扩展方向优先落在“优势控制域”内。当两个假定不能同时满足时，优先满足区域占优假定，但不应出现大于90°的转角。

2.   推移式滑坡的滑裂面扩展规律研究

土质边坡的滑裂破坏除了受到结构形式和荷载条件等外在因素影响外，土体自身的变形特性（模量）和强度特性（内摩擦角和黏聚力）也至关重要，可以归结为内在因素。本文基于人工设计理想边坡，对土体材料采用Mohr-Coulomb理想弹塑性本构模型，分析弹性模量和内摩擦角对土质边坡滑裂破坏过程的影响。

土坡模型和计算网格如图 5所示，地基厚2 m，长24 m；土坡底部长20 m，顶部长10 m，高10 m；在土坡顶部的边缘有一长2 m、高1 m的刚性块，用于加载制造土坡滑裂破坏，计算采用四结点四边形平面应变单元。对地基的底部施加全位移约束，对地基的左右两侧和土坡的右侧施加水平位移约束。

图  5  土坡模型和网格划分

Figure  5.  Soil slope model and meshing

下载: 全尺寸图片 幻灯片

土体基本材料参数为：弹性模量E为10 MPa，泊松比$\nu$为0.4，重度$\gamma$为20.0 kN/m³，内摩擦角φ为10°，黏聚力c为30 kPa，抗拉强度f_t为20 kPa。计算过程分为两个步骤：步骤1对整个模型施加重力荷载进行自重平衡计算；步骤2首先进行位移清零，然后在刚性块的中心施加竖直向下的位移荷载u= 0.5 m，分成100个增量步进行。

基本材料参数的模拟结果如图 6所示。滑裂面首先从顶部加载后缘产生，然后逐步向下发展，直至贯通坡面，形成整体滑动。在第22个加载步完成后滑裂面已经贯通滑体，此时滑坡体中上部表现为以沉降变形为主，而滑坡体中下部表现为以横向挤出变形为主。之后滑体将发生大滑移变形，边坡完全失稳破坏。

图  6  滑裂面形态和土坡变形

Figure  6.  Morphology of sliding surface and deformation of soil slope

下载: 全尺寸图片 幻灯片

2.1   均质土坡推移式滑动的滑裂面扩展规律

（1）弹性模量对均质土坡滑裂面的影响

分析不同弹性模量下均质土坡推移式滑裂面的滑裂深度和贯通所需增量步数的变化规律。选取弹性模量分别为5，10，15，20，30，50 MPa，其他参数不变，结果汇总如图 7所示。纵坐标为滑裂深度h（滑裂面最低点到土坡顶部的竖直距离）和滑裂面贯穿时的加载增量步数n，横坐标为弹性模量E。结果显示，弹性模量对均质土坡推滑破坏的滑裂面深度影响很小，不同的弹性模量，其他条件相同的均质土坡得到几乎相同的滑裂深度。但弹性模量越大，滑裂面贯穿时的增量步数越少，即总的加荷量越小。这一规律说明均质土坡的推移破坏是由应力和强度控制，而非变形控制，也说明了本文的强度控制区域占优判别的扩展分析方法符合问题的实际物理机制。

图  7  不同弹性模量下滑裂面深度与位移荷载的变化

Figure  7.  Variation of depth of sliding surface and displacement load with elastic modulus

下载: 全尺寸图片 幻灯片

（2）内摩擦角对均质土坡滑裂面的影响

分析不同内摩擦角下均质土坡推移式滑裂面的滑裂深度和贯通所需增量步数的变化规律，选取内摩擦角的变化区间为5°~26°，其他参数不变。

计算结果汇总如图 8所示，其中纵坐标分别为滑裂面深度h和滑裂面贯穿时的增量步数n，横坐标为内摩擦角φ。由图 8可见，均质土坡推移式滑动的滑裂面深度随着土体内摩擦角的增大而近似线性增加，且滑裂面贯穿需要的位移荷载也近似线性增加，两者的斜率分别为0.28和1.30。

图  8  滑裂面深度与位移荷载随内摩擦角的变化

Figure  8.  Variation of depth of sliding surface and displacement load with internal friction angle

下载: 全尺寸图片 幻灯片

部分算例的具体滑裂面形态如图 9所示。由图 9可以看出，较大的内摩擦角会产生更深的滑裂路径，提高边坡抗滑力和整体稳定性，结果符合一般认知规律。

图  9  不同内摩擦角下的滑裂面形态

Figure  9.  Variation of morphology of sliding surface with internal friction angle

下载: 全尺寸图片 幻灯片

（3）滑裂面扩展过程中反力的变化

选取位于刚性块中心正下方的土体单元作为对象，研究滑裂面扩展过程中反力的变化。设计了3个计算方案：方案1，E = 10 MPa，φ = 10°；方案2，E = 10 MPa，φ = 24°；方案3，E = 50 MPa，φ = 10°。计算结果见图 10，其中n为计算的增量步数（3种方案的总增量步数均为120步，总位移荷载会根据弹性模量确定，以保证每一增量步内的外力增量保持不变），p为当前增量步下的滑裂面深度与总裂缝深度的比值，$\sigma$为对象单元的竖向应力。

图  10  滑裂面扩展过程中反力的变化图

Figure  10.  Variation of counter force during surface expansion of sliding

下载: 全尺寸图片 幻灯片

由图 10可知，在滑裂面开裂前和刚开裂的一段时间内，对象单元的竖向正应力随着增量步数的增加而增加；当滑裂面扩展到一定深度后，滑裂面向前的扩展总是伴随着反力的减小。这说明推移式滑坡演化过程至少分为2个阶段：后缘压缩阶段和加速变形阶段。对于土体模量较小的情况，在加速变形阶段可能会产生多个变形台阶，滑裂面贯通的突发性较弱。对于土体模量较大的情况，加速变形阶段会使滑裂面直接贯穿坡体，边坡失稳的突发性较强。内摩擦角越大，在加速变形阶段越可能产生的变形台阶数量越多，有利于减弱边坡失稳的突发性。

2.2   成层土坡推移式滑动的滑裂面扩展规律

设计材料分界面分别距离土坡顶部2.5，7.5 m的两个成层土坡。计算过程分为两个步骤：①对整个模型施加重力荷载；②分成120个增量步，在刚性块的中心施加竖直向下的位移荷载u = 0.3 m。

（1）弹性模量对成层土坡滑裂面的影响

地基中土体内摩擦角φ为16°。将位于材料分界面上部的材料的弹性模量记为E₁，位于材料分界面下部的材料的弹性模量记为E₂。对上下分界面分别计算“E₁ = 10 MPa，E₂ = 50 MPa”和“E₁ = 50 MPa，E₂ = 10 MPa”两种工况，其他参数不变，得到滑裂面形态如图 11所示。

图  11  两种成层土坡不同弹性模量的计算结果

Figure  11.  Calculated results of different elastic moduli of two stratified soil slopes

下载: 全尺寸图片 幻灯片

比较图 11（a），（b）与图 11（c），（d）可知：如果上部土体的弹性模量比下部土体的小，则滑裂面深度较小；反之，则滑裂面深度较大。比较图 11（a），11（c）及图 11（b），（d）可知：土体为上软下硬的情况时，上部土层越厚，则滑裂面的深度越大；土体为上硬下软的情况时，上部土层越厚，则滑裂面的深度越小。这一规律与刚性基础基底应力分布规律类似，即上软下硬的边坡发生应力集中，破坏集中在较小范围内；而上硬下软的边坡发生应力扩散，破坏范围也相应扩大。

图 11（b）中，在材料分界面处滑裂面路径发生了较大的变化，是由于这一区域发生了局部张拉破坏，且滑裂面并非是由上缘向下单向发展而成的，这将在第3节中详细分析。

（2）内摩擦角对成层土坡滑裂面的影响

将成层土坡上部土体的内摩擦角记为φ₁，下部土体的内摩擦角记为φ₂。对两个成层土坡分别计算“φ₁ = 10°，φ₂ = 26°”和“φ₁ = 26°，φ₂ = 10°”时的推滑破坏过程，共4种工况，其他参数和方案设置不变，计算得到最终贯通后的滑裂面形态如图 12所示。由图 12可见，与均质土坡类似，整体较大的内摩擦角会产生较深的滑裂面，从而提高边坡的整体安全系数。对于成层土坡，内摩擦角大的土层占总体比例越大（相当于整体的加权平均内摩擦角越大），则最终整体的滑裂深度越大，边坡越趋于安全。

图  12  两种成层土坡不同内摩擦角的计算结果

Figure  12.  Calculated results of different internal friction angles of two stratified soil slopes

下载: 全尺寸图片 幻灯片

3.   成层土坡中张拉裂缝的分析

第2节在不同弹性模量的成层土坡推滑破坏过程中，分界面附近出现了张拉型裂缝，如图 11（b）所示。对此过程中的应力分布变化进行分析，分别绘出此工况下滑裂面扩展过程中5个典型阶段的大主应力和小主应力分布，如图 13，14所示。

图  13  大主应力等值线图和滑裂面发展轨迹（kPa）

Figure  13.  Contours of major principal stress and trajectory of slding surface (kPa)

下载: 全尺寸图片 幻灯片

图  14  小主应力等值线图和滑裂面发展轨迹（kPa）

Figure  14.  Contours of minor principal stress and trajectory of sliding surface (kPa)

下载: 全尺寸图片 幻灯片

在滑裂面出现之前，刚性块施加的位移荷载使土坡不断向下压缩，由于材料分界面上方土体的弹性模量比下方土体的大，所以下方土体的变形会更大，而在分界面处两者要满足变形协调条件，导致上方土体内产生拉应力。当拉应力超过土体的抗拉强度，就会发生张拉裂缝。之后随着进一步加载和滑裂面发展，滑裂面尖端远离分界面，转变为剪切破坏。

综上分析，图 13（a）~（e），14（a）~（e）给出的滑裂面发展过程为：①加载初期，分层界面上下大主应力分布较为均匀，是以自重应力为主，而上层较硬土体的底部因为变形不协调已经开始出现横向拉应力；②上层土体底部附近出现张拉裂缝；③张拉裂缝向两端扩展；④张拉裂缝两端的应力分布发生变化，剪切应力开始起主导作用，裂缝方向也随之发生转折；⑤剪切裂缝继续发展，向上贯通到破顶加载刚性块的角部附近，向下也继续延伸，直至最终形成贯穿整个滑体的滑裂面。

4.   结论

本文基于XFEM开发了能够模拟土体滑裂面扩展过程的程序。通过对人工设计均质土坡和成层土坡的推滑破坏过程进行模拟分析，并比较土体弹性模量和内摩擦角的影响，得到以下6点结论。

（1）坡顶堆载作用下，滑裂面从顶部加载后缘产生，然后逐步向下发展，直至贯通坡面，形成整体滑动。

（2）弹性模量对均质土坡推移式滑动的滑裂面深度的影响很小，但弹性模量越大，滑裂面贯穿的所需的加载量越少，即更容易形成贯通滑裂。

（3）均质土坡推移式滑动的滑裂面深度随着土体内摩擦角的增大而增大，滑裂面贯穿需要的加载量也随之增大。内摩擦角同时影响着滑裂面的扩展角度和开裂时机。

（4）推移式滑坡演化过程可分为后缘压缩和加速变形两个阶段。土体模量较小时，在加速变形阶段可能会产生多个变形台阶，滑裂面贯通的突发性较弱；土体模量较大时，加速变形阶段滑裂面较快贯通，失稳的突发性较强。内摩擦角越大，则在加速变形阶段可能的变形台阶数越多，从而降低失稳的突发性。

（5）对于成层土坡，上软下硬的边坡发生应力集中，破坏集中在较小范围内，滑裂深度较小；而上硬下软的边坡发生应力扩散，破坏范围也相应扩大，滑裂深度较大。成层土坡的整体（加权平均）内摩擦角越大，则滑裂深度越大。

（6）成层土坡在层间模量相差较大的情况产生局部变形不协调而发生局部的张拉型破坏，进而影响整体的边坡滑裂破坏过程。

这些规律均符合一般认知，也验证了本文方法可以准确、直观地模拟滑裂面的发生发展过程，并提供细观机理上的认识。扩展有限元方法在岩土工程问题中值得进一步更为广泛的应用研究。