边坡的渗流–应力耦合数学模型与稳定性分析

Mathematical Model of Seepage-Stress Coupling and Stability Analysis of Slope

摘要: 针对松散堆积体边坡在降雨与开挖耦合作用下的稳定性问题，本文建立了饱和–非饱和渗流与应力耦合的偏微分方程模型。基于虚功原理与有效应力原理，推导了渗流–应力耦合控制方程，采用FLAC3D平台进行三维数值模拟，并结合强度折减法构建稳定性判据。通过对比单一开挖、单一降雨及耦合工况下的数值试验，揭示了耦合效应对塑性区扩展与安全系数的协同放大效应。结果表明，耦合作用下边坡安全系数显著降低至1.18，揭示了渗流–应力反馈机制在边坡灾变过程中的数学本质。

Abstract: To investigate the stability of loose accumulation slopes under the coupled effects of rainfall and excavation, this study establishes a set of partial differential equations describing the coupled saturated-unsaturated seepage and stress fields. Based on the principles of virtual work and effective stress, the governing equations for seepage-stress coupling are derived. Three-dimensional numerical simulations are performed using the FLAC3D platform, and a stability criterion is developed through the strength-reduction method. By comparing numerical experiments under excavation only, rainfall only, and coupled conditions, the study reveals a synergistic amplification effect of the coupling process on the expansion of the plastic zone and the reduction of the safety factor. The results indicate that the slope safety factor decreases significantly to 1.18 under coupled conditions, highlighting the mathematical essence of the seepage-stress feedback mechanism in slope failure evolution.

1. 引言

工程勘察表明，松散堆积体边坡是山区线路建设中的关键挑战。该类边坡由古滑坡体、崩坡积体等结构性差的物质构成，其松散多孔的特征使其易受降雨入渗软化，同时自稳能力差，对开挖卸荷等扰动极为敏感。因此，描述渗流–应力耦合作用的偏微分方程组成为表征并预测此类边坡失稳的关键数学模型。典型案例如箭丰尾公路，在持续强降雨条件下进行切坡开挖，最终触发大规模滑坡，造成重大社会经济损失[1]。由此可见，建立精确的数学模型以探究降雨条件下开挖边坡的稳定性演化规律，对灾害防控具有重要的科学与实践价值。

隧道开挖引发的卸荷效应直接导致岩体变形与力学性质劣化，显著影响边坡工程的稳定性与安全[2]。为了探究隧道开挖作用下影响边坡稳定性的因素，魏燃等[3]利用FLAC3D软件对软弱围岩浅埋段隧道施工进行数值模拟，揭示了边坡渐进式失稳的演化机制。邓东平等[4]将极限平衡推力法与迭代计算策略相结合，实现了安全系数与滑面强度参数的动态交互分析。穆成林等[5]则依据岩体结构类型，筛选定性与定量指标，引入未知测度理论构建了新的边坡稳定性预测体系。葛虎胜等[6]结合典型工程案例，系统阐明了高陡边坡在开挖过程中的位移演化规律及稳定状态变化。在降雨入渗方面，水分迁移会显著削弱土体的黏聚力与内摩擦角，降低抗剪强度，并引起孔隙水压力上升，从而加剧边坡失稳风险[7]。姚茂宏等[8]在经典Green-Ampt模型的基础上，引入非均匀初始含水率条件，建立了考虑饱和区渗流与气压作用的入渗模型，并推导出对应的边坡稳定性系数解析表达式。顾玉明等[9]基于降雨入渗理论与Geo-Studio平台，通过设定不同渗透系数比值，探讨了渗透各向异性对边坡渗流场及稳定状态的响应规律。任东兴等[10]借助ABAQUS模拟了降雨入渗区内土体强度参数的弱化过程，重点分析了暴雨条件下边坡的应力–变形特征及其稳定性演变。陈景松等[11]基于饱和–非饱和渗流理论与强度折减有限元方法，研究了不同降雨类型(均匀型、间歇型)及岩土渗透性对混合花岗岩路堑边坡渗流特性与稳定性的影响机制。

开挖卸荷与降雨入渗的多场耦合作用是诱发边坡非线性失稳演化的重要机制，其过程机理复杂，是当前岩土工程领域的难点。张治国等[12]基于分层Green-Ampt模型模拟降雨入渗过程，结合复变函数理论对雨浸地层进行分区映射，提出了降雨影响下隧道开挖引起的地层位移场与应力场理论解析方法。熊齐欢等[13]通过室内直剪试验，系统分析了不同饱和度条件下坡积粉质黏土与全风化凝灰熔岩的抗剪强度变化，建立了黏聚力与内摩擦角随饱和度变化的数学关系。陈亚烽等[14]基于瑞典条分法与一阶线性应变软化模型，提出了一种边坡渐进破坏分析方法，推导了应变软化型边坡的极限平衡表达式，并获得了不同破坏阶段的安全系数。

综上所述，尽管现有研究在降雨入渗与施工扰动对岩土体力学性能的影响方面已形成较系统的理论框架，但多数成果仍局限于单一因素或特定工况的数值模拟与经验分析，缺乏对开挖与降雨动态耦合作用下边坡失稳演化过程的系统性数学建模与理论分析。为此，本文针对松散堆积体边坡，旨在建立一个渗流–应力全耦合的数学模型，开展开挖与降雨耦合作用下的失稳演化规律研究。采用FLAC3D三维显式有限差分法求解该耦合偏微分方程组，结合强度折减法作为稳定性判据，系统分析开挖扰动与降雨入渗共同作用下的位移发展规律、塑性区与剪应变增量分布特征。从数学上揭示耦合工况下边坡的破坏模式与失稳演化规律，以期为地质灾害的防控提供理论参考。

2. 数学模型与数值方法

2.1. 渗流–应力耦合控制方程

边坡系统的稳定性由渗流场与应力场的耦合作用所控制，该过程可由一组偏微分方程描述[15] [16]。

2.1.1. 应力场平衡方程

假设土体为连续介质，根据虚功原理，应力平衡方程的等效积分形式为：

${\displaystyle {\int }_V\left(\sigma :\delta \epsilon \right)\text{d}V=}{\displaystyle {\int }_V\left(f\cdot \delta u\right)\text{d}V+{\displaystyle {\int }_S\left(p\cdot \delta u\right)\text{d}S}}$ (1)

式中，σ为总应力张量；δε为虚应变张量；δu为虚位移向量；f为单位体积力向量，对于非饱和土为土骨架重与水重之和，即f = (1 − n)ρ_sg + nsρ_wg，其中n为孔隙率，g为重力加速度，ρ_s与ρ_w分别为土颗粒和水的密度，s为饱和度；V为体积积分域；S为边界面积分域。

2.1.2. 有效应力原理与本构关系

采用Bishop有效应力原理来描述非饱和土：

$\sigma ={\sigma }^{\prime }-\chi s{p}_{w}m$ (2)

式中，${\sigma }^{\prime }$ 为有效应力张量；p_w为孔隙水压力；s为饱和度；m为单位张量，对于三维问题，m = [1, 1, 1, 0, 0, 0]^T；χ为有效应力参数。

假设土骨架为线弹性材料，有效应力增量与应变增量满足线性本构关系：

$d{\sigma }^{\prime }=D:d\epsilon$ (3)

式中，D为弹性刚度张量；dε为应变增量张量。

2.1.3. 渗流场方程与耦合项

土体变形时，土颗粒体积不变，孔隙率n随体积应变ε_v演化：

$n=\frac{{\varphi }_0+{\epsilon }_v}{1+{\epsilon }_v}$ (4)

式中，${\varphi }_0$ 为初始孔隙率；ε_v为体积应变。

考虑降雨入渗的质量源项，边坡内水的质量守恒方程为：

$\frac{\partial }{\partial t}\left(ns{\rho }_w\right)+\nabla \cdot \left({\rho }_{w}v\right)=Q$ (5)

式中，n为孔隙率；s为饱和度；ρ_w为水的密度；v为渗流速度向量；t为时间；▽为散度算子；Q为降雨入渗质量源项。

渗流速度与总水头梯度成正比，考虑重力的达西定律为：

$v=-\frac{1}{{\gamma }_w}k\cdot \nabla \left(p_w+{\rho }_{w}gz\right)$ (6)

式中，v为渗流速度向量；k为渗透系数张量；γ_w为水的容重；z为高程水头；g为重力加速度；$\nabla$ 为梯度算子。

关键耦合项出现在孔隙率n上，孔隙率随土体变形而变化：

$n={\varphi }_0+{\epsilon }_v$ (7)

其中：

${\epsilon }_v=\nabla \cdot u$ (8)

将有效应力原理和本构关系代入虚功原理的增量形式，得到含孔隙水压力增量的耦合方程：

${\displaystyle \int {}_V\left(\delta \epsilon :D:d\epsilon \right)\text{d}V}-{\displaystyle \int {}_V\left(\delta \epsilon :\left(sd{p}_{w}m\right)\right)\text{d}V}=d$ (9)

将达西定律代入渗流连续性方程，并结合孔隙率演化方程，得到耦合渗流方程：

$\frac{\partial }{\partial t}\left(ns{\rho }_w\right)=\nabla \cdot \left(\frac{{\rho }_w}{{\gamma }_w}k\cdot \nabla \left(p_w+{\rho }_{w}gz\right)\right)+Q$ (10)

变形控制方程与渗流控制方程共同构成了渗流–应力全耦合的封闭数学模型，清晰刻画双向耦合机制。这一理论框架是定量分析隧道开挖扰动与降雨入渗耦合下，边坡位移场、剪应变演化及塑性区扩展的核心基础，为后续开展边坡稳定性的数值模拟及防控措施制定提供了严谨的理论支撑。

2.2. 强度折减法原理

在边坡稳定性分析中，我们采用强度折减法作为系统失稳的数学判据。该方法由Zienkiewicz等[17]提出并发展。边坡滑动破坏的本质是岩土体或结构面抗剪强度不足引发的强度失效[18]。计算过程中，在恒定重力场条件下，通过连续折减岩土体的粘聚力c和内摩擦角φ参数，实现边坡渐进破坏过程的模拟。具体数学表达式如下：

$c^{\prime }=\frac{c}{F_c}$ (11)

${\phi }^{\prime }=\arctan \left(\frac{\tan \phi }{F_{\phi }}\right)$ (12)

式中，c为岩土体的实际粘聚力；$c^{\prime }$ 为折减后岩土体的粘聚力；φ为岩土体的实际摩擦角；${\phi }^{\prime }$ 为折减后岩土体的摩擦角；F_c和F_φ为强度折减系数。

该方法的数学原理在于，通过逐步增大强度折减系数值，求解耦合系统在折减后参数下的响应。当数值计算不再收敛，或系统位移发生突变时，即表明该系统解趋于发散，对应于边坡的临界破坏状态。此时对应的折减系数被定义为该边坡系统的安全系数，其数学意义是系统保持稳定状态的强度储备下限。

3. 数值模型与实验设计

为求解上述数学模型并进行数值实验，需建立相应的几何模型并设定初边值条件。数值模型设置参考了张鲁渝等[19]提出的建模准则。根据该准则，当坡脚至左侧边界的水平距离取坡高的1.5倍、坡顶至右侧边界的水平距离取坡高的2.5倍，且模型上下边界的总高不小于坡高的2倍时，可获得较为理想的计算精度。针对净长45 m、净宽5 m、净高4 m的隧道模型，建立了如图1(a)所示的理想化边坡数值模型。模型具体流体边界设置如下：CD边界为透水边界，允许地下水自由通过；DE边界为不透水边界，模拟隔水层以阻止水流迁移；AB边界兼具入渗与溢流功能，既可作为降雨入渗通道，也可在孔隙水压力升高时作为水体排出界面；BC边界则设为潜在溢出边界，当地下水位高于该边界高程时，地下水将由此自然排泄。为构建稳定渗流条件，模型左右两侧分别设置固定水头边界，其中左侧水头为13  m，右侧水头为25  m，从而形成自右向左的水力梯度，为系统提供稳定的渗流驱动力。

模型坐标系的定义为：X轴指向水平，Z轴指向竖直，Y轴沿隧道开挖方向。三维模型的网格划分如图1(b)所示，总共包含79,370个单元和14,604个节点。在边界条件的设置上，坡顶采用自由边界以模拟自然地表的受力状态；前、后、左、右四个侧面施加法向位移约束，限制其侧向变形；边坡底部则施加X、Y、Z三个方向的法向位移约束，以确保底部稳定，防止模型发生整体刚体运动。模型中各岩土体的物理力学参数主要依据现场地质勘察成果，并参考相关文献资料综合确定，具体数值列于表1。

Figure 1. Schematic diagram of an ideal slope numerical model

图1. 理想边坡数值模型示意图

Table 1. Geotechnical material parameters

表1. 岩土材料参数

岩土名称	密度/kg∙m^−^3	杨氏模量/Mpa	泊松比	黏聚力/kpa	内摩 擦角/˚	抗拉强度/Mpa	孔隙率/%	渗透系数/cm∙s^−1
黏土	1900	80	7	0.32	19	0.01	0.2	5 × 10^−^5
基岩	2770	50,000	5	0.37	47	8	0.5	3 × 10^−^11

4. 数值模拟结果与分析

为系统揭示降雨条件下隧道开挖对边坡稳定性的影响规律，本节设置了三种典型工况进行对比分析，分别为单一隧道开挖工况、单一降雨条件以及开挖与降雨耦合工况。在隧道开挖工况中，采用3米作为每一施工步的开挖步距，该尺寸的设定综合考虑了常见隧道施工工艺的实际操作步距，同时也兼顾数值模拟的合理性与计算效率，能够较好地反映隧道分步开挖过程中围岩应力分阶段释放的动态响应特征。降雨条件设定为0.5 mm/h的强度持续6小时，该降雨强度属于中雨范畴，能够代表一种较为常见且具有潜在诱发滑坡风险的降雨事件。通过对比不同条件下边坡的力学响应与变形特征，重点研究在开挖与降雨耦合作用下边坡的位移场分布规律、塑性区扩展特征、最大剪应变增量场的空间演化。在此基础上，结合定性与定量分析结果，综合分析边坡在耦合作用下的稳定性演化规律，为类似地质的隧道工程设计、稳定性评价与灾害防控提供参考。

4.1. 开挖工况下边坡稳定性分析

开挖步距3 m工况下，边坡的数值模拟结果如图2所示。

Figure 2. Numerical simulation results for the excavation condition

图2. 开挖工况数值模拟结果

根据数值模拟结果(图2(a)~(d))可知，隧道开挖引起的扰动主要集中在拱顶上方的软弱黏土层。施工完成后，边坡顶端与隧道结构均出现明显变形，其中拱顶及上覆黏土层位移最显著，表明该区岩体强度低、约束性差；而下伏坚硬基岩整体保持稳定。在3 m开挖步距下，塑性区主要分布于黏土层与板岩层交界处，并延伸至拱顶周围，坡体顶部出现局部拉张单元。剪应力虽在较大范围内集中，但未形成贯通滑移面，说明由于软弱层与坚硬层刚度差异导致的应力重分布已达到新的平衡，边坡整体仍处于稳定状态。同时，最大剪应变增量集中于拱顶区域，表明该处为应力集中和变形发展的关键部位。该工况下的安全系数为1.77，表明在纯力学扰动下，系统具有较高的稳定性。

从力学机制角度分析，隧道开挖导致的卸荷效应引起总应力下降，依据有效应力原理，土体有效应力相应降低，从而削弱其抗剪强度。在该工况中，应力路径主要表现为卸荷回弹特征，塑性区的形成是应力重分布超过土体峰值强度的直接结果，体现出开挖过程中能量释放与岩土体强度弱化的响应过程。

4.2. 降雨工况下边坡稳定性分析

降雨强度为0.5 mm/h、持续时间6小时的工况下，边坡的数值模拟结果如图3所示。

Figure 3. Numerical simulation results for the rainfall condition

图3. 降雨工况数值模拟结果

根据数值模拟结果(图3(a)~(d))可知，边坡在降雨工况下的位移响应表现出明显的岩性控制特征，整体位移极值(36.051 mm)始终分布于力学性质较差的黏土层区域，揭示该土层对降雨入渗具有更高的敏感性。受低降雨强度条件制约，边坡整体位移场未发生显著变化，表明降雨入渗虽在一定程度上激活了黏土层的变形，但尚未对边坡整体稳定性构成主导性影响。在持续6小时降雨条件下，塑性区逐步向深部扩展，但其发育范围仍相对有限，主要集中于黏土层内部，坡脚区域未形成贯通的塑性滑移带。与此同时，最大剪应变增量分布进一步表明潜在滑移面尚未贯通。该工况下边坡稳定性安全系数为1.47，综合表明边坡整体仍处于基本稳定状态。

降雨入渗侵蚀了边坡土体的基质吸力，依据扩展的莫尔–库仑准则，这直接导致土体有效黏聚力与抗剪强度显著劣化。同时，水分迁移在黏土层中诱发应变局部化，促使塑性区向深部扩展。上述机制共同证明，由降雨触发的强度劣化是驱动边坡变形加剧与稳定性下降的内在力学根源。

4.3. 耦合工况下边坡稳定性分析

在降雨作用下，裂隙周边渗流场发生显著改变，水分入渗导致土体基质吸力下降，进而引起强度与刚度劣化，促使破坏区域沿裂隙扩展；随着裂隙持续发育，坡体内部约束逐渐减弱，此时叠加隧道开挖扰动，将进一步加剧变形响应。当破坏范围扩大且变形累积至临界状态时，潜在滑动面将逐步贯通，诱发整体失稳。本节设置降雨强度为0.5 mm/h、持时6 h，并在降雨过程中以3 m步距进行隧道开挖，系统分析边坡在耦合工况下的变形规律与稳定状态演化特征，数值模拟结果如图4所示。

Figure 4. Numerical simulation results for the coupled condition

图4. 耦合工况数值模拟结果

隧道开挖引起的应力释放与降雨入渗导致的土体强度弱化共同影响着边坡稳定性。数值模拟结果(图4(a)~(d))显示，耦合工况下边坡位移最大值出现在黏土层、基岩夹层及隧道拱顶区域。塑性区分布表明，滑体上部单元以拉剪复合破坏为主，滑移面附近则以剪切破坏为主导。最大剪应变增量云图显示坡体内未形成贯通剪切带，但结合强度折减法所得安全系数为1.18，表明边坡处于临界稳定状态，存在失稳风险。

耦合效应的力学本质在于应力场与渗流场的交互作用。开挖扰动了原场地应力，在拱顶和坡肩形成应力松弛区，利于雨水汇集和入渗。孔隙水压力的升高(或基质吸力丧失)直接降低了潜在滑动面上的有效法向应力，从而显著削弱其抗剪能力。同时，塑性应变的发展局部改变了土体的结构性和渗透性，加速了水分的迁移与软化的进程。这种“应力重分布–渗流场演进–强度劣化–塑性区扩展”的动态耦合过程，构成了安全系数显著降低至1.18的内在机制。

开挖与降雨的耦合对边坡稳定性的影响远超过单一工况。该过程通常始于隧道开挖引发的卸荷扰动使潜在弱面形成，继而是降雨带来的岩土体饱水软化及有效应力降低。在此过程中，渗流与应力场相互形成正反馈，渗流使应力重分布，应力变化又引导渗流发展，共同加速塑性区与剪应变区的扩展。尽管在耦合作用下塑性区与剪应变区尚未完全贯通，但其范围的显著扩大已导致位移量激增，并使安全系数降至1.18的临界状态，充分凸显开挖–降雨耦合在边坡灾变过程中的关键控制作用，对工程安全预警与风险防控具有重要参考价值。

5. 讨论

本文通过建立渗流–应力耦合的偏微分方程模型，并结合三维数值模型，定量分析了开挖与降雨耦合作用下松散堆积体边坡的稳定性演化规律。然而，本研究仍存在一定的局限性，需要在未来的工作中加以改进：

1) 设土体为连续、均质且各向同性的线弹性材料，并采用了简化的Bishop有效应力公式。实际的松散堆积体通常表现出显著的非均质性、各向异性及非线性弹塑性，未来的模型应考虑更复杂的本构关系以更真实地反映岩土体的力学行为。

2) 管模型边界设置参考了经典准则，但实际工程中的地质与水文地质条件更为复杂。初始地应力场简化为重力场，未考虑构造应力。初始孔隙水压力分布采用稳定渗流场假设，可能与实际情况存在差异。

6. 结论

为深入探讨隧道开挖与降雨作用下松散堆积体边坡稳定性的耦合效应，本文建立了一个渗流–应力全耦合的偏微分方程数学模型，并采用数值方法进行求解与分析，分别对单一降雨入渗、单一隧道开挖及开挖与降雨耦合工况进行了系统的数值模拟。基于本文设定的特定模型参数与工况条件，通过模拟分析，得出了以下主要结论：

1) 建立的渗流–应力耦合数学模型能够有效表征开挖与降雨耦合作用下边坡的稳定性演化。在本研究设定的工况下，耦合作用通过渗流–应力场的相互作用显著削弱边坡稳定性，其安全系数降幅远高于单一工况。因此，在类似地质条件的边坡稳定性评价中，应该综合考量二者的耦合效应，以准确评估风险并制定防控措施。

2) 合效应的数学物理本质在于力学损伤与水力劣化的非线性交互作用。隧道开挖形成的卸荷区与裂隙网络，为雨水入渗提供了优势路径；水分入侵导致基质吸力丧失和孔隙水压力升高，根据有效应力原理，直接降低了有效应力和抗剪强度，加剧了岩土体的软化和塑性变形。这一相互加剧的正反馈循环过程，是驱动边坡失稳性急剧下降的主要原因。

3) 研究通过数学建模与数值实验，定量揭示了开挖–降雨耦合作用的非线性放大机制。强度折减法作为稳定性判据，其数学意义在于识别了耦合系统在参数空间中的临界状态。该框架为边坡稳定性分析提供了基于偏微分方程控制的定量理论工具。