含裂隙土质隧道降雨入渗双通道渗流模型

邱军领; 贾玎; 赖金星; 唐琨杰; 强磊

doi:10.11779/CJGE20231215

含裂隙土质隧道降雨入渗双通道渗流模型 English Version

长安大学公路学院, 陕西西安 710064

基金项目:

国家自然科学基金青年基金项目 52208386

国家自然科学基金面上项目 51978066

国家自然科学基金面上项目 52278393

长安大学中央高校基本科研业务费专项资金项目 300102213202

详细信息

作者简介:
邱军领(1989—)，男，博士，副教授，主要从事隧道及地下工程方面的研究工作。E-mail：junlingqiu@chd.edu.cn

通讯作者:
赖金星, E-mail：laijinxing@chd.edu.cn

中图分类号: U451
计量
- 文章访问数: 314
- HTML全文浏览量: 27
- PDF下载量: 79
出版历程
- 收稿日期: 2023-12-10
- 网络出版日期: 2024-05-29
- 刊出日期: 2025-02-28

Dual-channel seepage model for tunnels with fissured soil under rainfall infiltration

School of Highway, Chang'an University, Xi'an 710064, China

摘要

摘要: 通过含单一裂隙土层局部渗流性态分析中得到的现象和结果，提出双通道渗流模型理论。基于阶跃函数推导了降雨诱发渗流的边界条件及其转换公式，以既有文献中的降雨入渗模型为基础，讨论了阶跃函数过渡区长度的设定，结果发现过渡区长度为4或8个单位时，计算结果较为准确。应用上述结果，通过COMSOL Multiphysics多物理场耦合方法模拟了双通道渗流效应，并探究降雨入渗情况下含裂隙土体的隧道围岩变形特征。同时根据可能影响隧道降雨入渗变形特征的因素包括裂隙位置、降雨强度、裂隙宽度进行多因素影响性分析。结果表明：随着降雨时间的增长，近裂隙一侧隧道拱顶沉降和上部水平位移值明显高于远裂隙侧；裂隙宽度越大，降雨入渗扩展越快，隧道产生的位移越大。对比2 mm裂隙宽度的拱顶沉降和上部水平位移，隧道上方存在8 mm宽度裂隙时，拱顶沉降最大值达328 mm，增幅达6.5%，上部水平位移最大值为26 mm，增幅达44%；降雨强度则对隧道围岩的影响较小。
- 含裂隙土质隧道 /
- 降雨入渗 /
- 双通道渗流效应 /
- 阶跃函数 /
- 变形特征
Abstract: The theory of a dual-channel seepage model is proposed based on the phenomena and results obtained from the analysis of local seepage behavior in soil layers with a single fissure. By using the step function, the boundary conditions and transformation formula for rainfall-induced seepage are derived. Based on the existing rainfall infiltration models in literatures, the setting of the transition zone length for the step function is discussed. It is found that the calculated results are more accurate when the transition zone length is 4 or 8 units. Applying the above results, the dual-channel seepage effects are simulated using the COMSOL Multiphysics, and the deformation characteristics of the surrounding rock of tunnel with fractured soil under rainfall infiltration are explored. At the same time, a multi-factor impact analysis is conducted based on the factors that may affect the deformation characteristics of rainfall infiltration of the tunnel, including fissure location, rainfall intensity and fissure width. The results show that with the increase of rainfall time, the settlements of the arch and upper horizontal displacement of the tunnel at the side near the fissure are significantly higher than those at side the far from the fissure. The larger the fissure width, the faster the rainfall infiltration and expansion, and the greater the displacement generated by the tunnel. Comparing the settlement of the arch crown with a fissure width of 2 mm and the upper horizontal displacement, when there is a fissure width of 8 mm above the tunnel, the maximum settlement of the arch crown reaches 328 mm, an increase of 6.5%, and the maximum horizontal displacement of the upper part is 26 mm, an increase of 44%. The intensity of rainfall has a relatively small impact on the surrounding rock of the tunnel.
- tunnel with fissured soil /
- rainfall infiltration /
- dual-channel seepage effect /
- step function /
- deformation characteristic

HTML全文

0. 引言

边坡是人类活动中极为重要的自然地质环境之一。如何保证边坡稳定成为了岩土工程一个重要的研究领域^[1]。在进行边坡稳定分析时，最主要的判断依据就是稳定安全系数。根据安全系数的大小来评判边坡是否稳定，这就涉及安全系数的定义问题。由于各种分析方法所依据的安全系数定义的不同^[2]，对于同一个计算对象，其计算出来的安全系数值一般也有所不同。但现有的一些安全系数的定义式中，很少考虑到滑动的矢量特征，且对边坡真实状态引入了过多的人为假定。为了克服这方面的局限性，葛修润院士从滑动的基本概念出发，指出滑动本身就是一个矢量变化的过程，并从矢量的角度入手，可求出滑裂面的下滑力与抗滑力，从而进一步求出边坡的安全系数。这一求解方法，称之为矢量和分析方法（vector sum method, VSM）^[3-4]。矢量和分析方法是基于力是矢量这一基本概念，安全系数定义为坡体在受各种荷载（包括自重）作用下，潜在滑面上提供的抵抗坡体滑动的合力矢与作用在潜在滑动面上贡献滑动力的合力矢在滑体的整体滑动趋势方向上投影的比值。该方法具有明确的物理意义，考虑了滑动方向，能较好地反映边坡的整体稳定性。

矢量和分析方法是根据力的矢量特性提出的一种边坡稳定性分析方法，该方法在求到了滑体的合力后，通过假定主滑方向，认为滑体沿主滑动方向平面滑动，目前常用如下4种方法来确定主滑方向：①基于潜在滑面各点真实剪应力合力矢单位方向；②基于极限下滑力的合力矢方向；③基于边坡潜在滑面剪出点指向剪入点的位置矢量方向；④基于最小势能原理。然而，在滑坡发生时，滑动面通常呈现弧形、平面、折线或曲面，这意味着坡体滑动不仅涉及平动问题，还需考虑坡体的转动，因此本文针对坡体转动提出了一种新的边坡稳定性分析方法：力矩和法（torque sum method, TSM）。此方法定义边坡稳定安全系数F为：坡体在受各种荷载（包括自重）作用条件下，潜在滑动面所能提供的抵抗坡体绕某点转动的极限强度所提供的力矩之和与滑动体所受各种荷载（包括自重）所引发坡体绕这点转动的力矩之和的比值。力矩和法可以直接采用显式计算，可以对任意滑动面上的安全系数进行求解计算，具有清晰的物理含义，概念清楚，可以直接根据滑动面的实际应力分布以及滑动面所能提供的抗滑能力进行计算，弥补了矢量和方法不能考虑坡体可能发生转动失稳的不足。本文给出了力矩和法的定义表达式，适用于一般的边坡岩土体力学分析。运用本方法计算了两个典型的边坡算例，并与传统分析方法的计算结果进行了对比分析。

1. 基于力矩和法的安全系数分析方法

1.1 安全系数的定义

对于可能发生坡体转动的情况，参照矢量和方法，边坡稳定安全系数F可以通过计算坡体应力状态获得。力矩和方法的安全系数F定义为潜在滑动面所能提供的绕某转动中心的极限抗转动力矩的总和∑M_r与作用在滑动面上计算应力绕转动中心的滑动力矩的总和∑M_s之比，即

$F=\frac{\sum \boldsymbol{M}_{\text{r}}}{\sum \boldsymbol{M}_{\text{s}}} 。$

(1)

而极限抗转动力矩是由滑动面上的极限应力强度提供。由于力是一个矢量，不但有作用的大小，而且有作用的方向，所以式(1)中求和是考虑了力的矢量性质，遵循矢量求和运算法则合成的一个力矩和。

1.2 基本假定

对于一种稳定性分析方法，若要将其所涉及到的所有因素都考虑进去，问题将会是复杂和繁琐的^[5]。如果通过适当的假设，可以让分析方法更加的有针对性，并把握住问题的主要矛盾，进而降低问题的复杂性。对力矩和法，将建立在如下假定的基础之上：

（1）对某一特定的边坡进行稳定性分析的先决条件，是通过地质勘察，确定边坡的基本物理力学参数和边坡的载荷以及边坡的力学边界条件，而坡体内应力分布状态可由其它应力分析方法（如有限元法、离散元法、条分法等）计算获得。

（2）在二维边坡问题中，确定边坡所构成的计算范围为 $\mathit{Ω}$ ，假设在该计算范围内已知边坡的潜在滑动面L，构成滑动体的区域为S。

（3）在边坡潜在滑动面上，假设岩土体的材料强度特性服从Mohr-Coulomb准则，抗剪强度为

${\tau _{\text{f}}} = c + \sigma \tan \varphi 。$

(2)

（4）对绕某点转动的力矩和法安全系数F(x_p, y_p)定义为：转动点p坐标为(x_p, y_p)，滑动面上所提供抵抗坡体转动的各力沿转动点的力矩代数和∑M_r(x_p, y_p)与滑动面上计算应力沿转动点的力矩代数和∑M_s(x_p, y_p)的比值：

$F({x_{\text{p}}},{y_{\text{p}}}) = \frac{{\sum {{M_{\text{r}}}({x_{\text{p}}},{y_{\text{p}}})} }}{{\sum {{M_{\text{s}}}({x_{\text{p}}},{y_{\text{p}}})} }} 。$

(3)

（5）滑动面是滑动体与母体脱离发生绕某点转动的面，可以是任意形状。对于滑动（转动）的滑体，滑动面是约束边界。做力矩和计算时，滑动体假定为刚体。

1.3 滑动体转动中心的确定

当分析某一具体的潜在滑动体时，也即确定了其对应的潜在滑动面。作为平面刚体，刚体最基本的运动形式为沿某直线的平动和绕某点的转动，如果不考虑其后的运动轨迹，从静止时的运动趋势考虑，可以认为刚体只有平动和转动，而且其启动加速度在理论上是可以计算的。

对于边坡滑体，以往的知识认为滑体是沿滑动面滑动（或转动），笔者认为滑动面是边坡滑体失稳后，产生运动形成的轨迹，与边坡初始失稳时形成的失稳面（或破坏面）是有区别的，如果将滑动面的形式与位移或运动的形式相关联，则无法解释现有的分析滑动面为非圆弧或非直线的方法中，滑动体作为刚体时如何发生。当边坡初始失稳时，无论滑动面形式如何，滑动面只是滑动体的一个边界，其相互作用力可以认为是边界上的力边界条件。

因为刚体的基本位移是平动和转动，分析滑动体时，应该按照其基本的位移特征进行分析。矢量和方法考虑了滑动体的平动位移，认为滑动体将沿某一直线运动，并通过力矢量求和等相应的假设，获得其安全系数。作为对矢量和方法的一个补充，本文将只考虑滑动体的转动。对于同时考虑平动和转动的复杂情况，将在往后的工作中开展。

值得强调的是，这里所指的滑动体的平动和转动，均与滑动面无关，仅与滑动体本身和其受力有关，平动和转动也只是指滑动体在初始失稳时的位移趋势，与滑动体失稳后的位移和滑动面的形成无直接关系。实质上，边坡的失稳是指的滑动体在初始位移趋势上的失稳，并在考虑初始位移趋势上计算其失稳安全系数。

如何获得滑动体的初始位移趋势或初始转动中心，可以按照边坡稳定性分析的普遍方法，也即求解边坡最小安全系数的方法来获得。定义其初始位移趋势（或转动中心）是使滑动体具有最小安全系数时对应的转动中心。

仅考虑滑动体转动时，如图 1为边坡稳定性分析时的安全系数求解示意图。中可以通过有限元法或其它方法，计算得到潜在滑动面上任何一点i的正应力 ${\sigma _i}$ 和切应力 ${\tau _i}$ 。设滑动面上点i处岩土体内聚力为 ${c_i}$ ，内摩擦角为 ${\varphi _i}$ ，中 ${\alpha _i}$ 为滑动面在i点上的切线相对于坐标系x轴的正向夹角。

图 1 转动点搜索示意图

Figure 1. Schematic diagram of rotating point searching

下载: 全尺寸图片幻灯片

设点p(x_p, y_p)为转动点，为了寻找所计算滑动体具有最小安全系数的转动点，本文采用了较直观的网格（Grid）搜索法，转动点假设位于设定的一矩形搜索区域，矩形区域均布网格节点，作为试算的转动点。每一节点只需计算一次安全系数，就可通过比较获得区域内具有安全系数最小的节点。如果获得的节点位于矩形区域边界上，则需调整区域。该方法与圆弧滑动面用Grid方法搜索滑动圆心类似（图 2）。

图 2 安全系数求解示意图

Figure 2. Solution of factor of safety

下载: 全尺寸图片幻灯片

1.4 力矩和法的安全系数计算

力矩和法安全系数的计算包括转动力矩和计算与抗转动力矩和计算两方面：

（1）转动力矩和计算

滑动体上任一微段 ${{\Delta }}{l_i}$ 的正应力 ${\sigma _i}$ 和切应力 ${\tau _i}$ 是由滑体自重和作用在滑体上的外载荷所产生的，理论上滑动面（也是边界）上的所有正应力 ${\sigma _i}$ 和切应力 ${\tau _i}$ 与滑动体的自重和外荷载是一个静力平衡系。正应力 ${\sigma _i}$ 和切应力 ${\tau _i}$ 对转动中心求力矩并沿滑动面求和，就是滑体的转动力矩和。同时，岩床反作用在滑体上的力为 ${\sigma '_i}$ 和 ${\tau '_i}$ ，它们数值相等，方向相反。如果能够求到滑动体的重心，转动力矩和的求解可以直接求重心对转动中心的力矩与滑体上的外荷载对转动中心的力矩的总和。这与沿滑动面由 ${\sigma _i}$ 和 ${\tau _i}$ 求力矩和是一致的。

将滑动面上的应力对转动点力矩沿滑动面积分（或求和）即为式（3）中的分母项。

滑动面上任何一点i的切应力 ${\tau _i}$ 的斜率为 ${k_{\tau i}}$ ，所在直线方程为

$y={k}_{\tau i}x+{b}_{i} 。$

(4)

式中： ${k_{\tau i}} = \tan {\alpha _i}$ ， ${b_i} = {y_i} - {k_{\tau i}}{x_i}$ ，(x_i, y_i)为i点坐标。

转动点(x_p, y_p)到 ${\tau _i}$ 所在直线的距离为

${d}_{\tau i}=\frac{\mid {k}_{\tau i}{x}_{\text{p}}+\text{(}-1\text{)}{y}_{\text{p}}+{b}_{i}\mid }{\sqrt{{k}_{\tau i}{}^{2}+{\text{(}-1\text{)}}^{2}}} 。$

(5)

切应力 ${\tau _i}$ 提供转动力矩，力矩和为 $\sum _{i=1}^{n}{\tau }_{i}\text{Δ}{l}_{i}{d}_{\tau i}$ 。

滑动面上任何一点i的正应力 ${\sigma _i}$ 的斜率为 ${k_{\sigma i}}$ ，所在直线方程为

$y={k}_{\sigma i}x+{B}_{i} 。$

(6)

其中， ${k_{\sigma i}} = - \frac{1}{{{k_{\tau i}}}} = - \frac{1}{{\tan {\alpha _i}}}$ ， ${B_i} = {y_i} - {k_{\sigma i}}{x_i}$ ，(x_i, y_i)为i点坐标。

转动点(x_p, y_p)到 ${\sigma _i}$ 所在直线的距离为

${d}_{\sigma i}=\frac{\mid {k}_{\sigma i}{x}_{\text{p}}+(-1){y}_{\text{p}}+{B}_{i}\mid }{\sqrt{{k}_{\sigma i}{}^{2}+{(-1)}^{2}}} 。$

(7)

正应力 ${\sigma _i}$ 提供转动力矩，力矩和为 $\sum\limits_{i=1}^{n}{\sigma }_{i}\text{Δ}{l}_{i}{d}_{\sigma i}$ 。

滑动体的总的转动力矩和为

$\sum\limits_{i=1}^{n}{M}_{\text{s}}({x}_{\text{p}},{y}_{\text{p}}) ={\displaystyle \sum _{i=1}^{n}{\tau }_{i}\text{Δ}{l}_{i}{d}_{\tau i}} +\sum\limits_{i=1}^{n}{\sigma }_{i}\text{Δ}{l}_{i}{d}_{\sigma i} 。$

(8)

由于转动点的位置不同，滑动面的形状可能不规则，某一弧段的 ${\sigma _i}$ 可能提供下滑力矩也可能提供抗滑力矩，体现其计算力矩时可正可负（如果定义沿坡面倾向方向转动为正），但作为力矩和方法，不需要区分此时的应力是否是产生下滑力矩或抗滑力矩，只需将抗滑力矩设为负值，进行力矩和计算即可。在计算切应力 ${\tau _i}$ 时也如此，不需要区分，只需根据实际发生的计算应力进行计算即可。不需要区分任意弧段的 ${\sigma _i}$ 和 ${\tau _i}$ 产生的力矩的正负，更符合“和”方法的本意，而且在力学上，滑动面上的 ${\sigma _i}$ 和 ${\tau _i}$ ，等效于滑动体的重力和作用于滑动体上的全部外力，求和时结果是一致的。

（2）抗转动力矩和计算

岩土体的材料强度服从Mohr-Coulomb准则，滑动体沿滑动面发生剪切，其抗转动力矩则由滑动面上的切向抗剪强度和法向反力提供。对于滑动面上的任一微段 ${\Delta }{l_i}$ ，岩土体的抗剪强度 ${\tau _{{\text{f}}i}}$ 由莫尔库仑强度理论进行计算，方向与位移趋势方向相反，法向反力 ${\sigma '_i}$ = ${\sigma _i}$ ，方向相反。

${\tau _{{\text{f}}i}} = {c_i} + {\sigma _i}\tan {\varphi _i}。$

(9)

滑动面上切向抗剪强度和法向反力对转动点的力矩积分（或求和）构成了抗转动力矩和，即为式（3）中的分子项。

${\tau _{{\text{f}}i}}$ 所在直线与 ${\tau _i}$ 一致，故转动点(x_p, y_p)到 ${\tau _{{\text{f}}i}}$ 所在直线的距离也为 ${d_{\tau i}}$ ，提供抗转动力矩之和为 $\sum _{i=1}^{n}{\tau }_{\text{f}i}\text{Δ}{l}_{i}{d}_{\tau i}$ 。

转动点(x_p, y_p)到 ${\sigma '_i}$ 所在直线的距离也为 ${d_{\sigma i}}$ ，提供抗转动力矩之和为 $\sum\limits_{i=1}^{n}{{\sigma }^{\prime }}_{i}\text{Δ}{l}_{i}{d}_{\sigma i}$ 。

抗转动力矩和即为

$\sum\limits_{i=1}^n M_{\mathrm{r}}\left(x_{\mathrm{p}}, y_{\mathrm{p}}\right)=\sum\limits_{i=1}^n \tau_{\mathrm{f}i} \Delta l_i d_{\tau i}+\sum\limits_{i=1}^n \sigma_i^{\prime} \Delta l_i d_{\sigma i}。$

(10)

（3）安全系数的计算

由式（8），（10），并代入式（3），得到力矩和法安全系数的表达式（1）或（3）的计算式。

$F({x_{\text{p}}},{y_{\text{p}}}) = \frac{{\sum ({c_i} + {\sigma _i}\tan {\varphi _i}){\Delta }{l_i}{d_{\tau i}} + \sum {{\sigma '}_i}{\Delta }{l_i}{d_{\sigma i}}}}{{\sum {\tau _i}{\Delta }{l_i}{d_{\tau i}} + \sum {\sigma _i}{\Delta }{l_i}{d_{\sigma i}}}}\text{，}$

(11)

写成积分形式为

$F({x_{\text{p}}},{y_{\text{p}}}) = \frac{{\int {\left[ {({c_i} + {\sigma _i}\tan {\varphi _i}){d_{\tau i}} + {{\sigma '}_i}{d_{\sigma i}}} \right]{\text{d}}l} }}{{\int {({\tau _i}{d_{\tau i}} + {\sigma _i}{d_{\sigma i}}){\text{d}}l} }}。$

(12)

1.5 力矩和方法求解步骤

力矩和方法的基本计算步骤如下：

（1）根据边坡地质勘查资料，确定边坡的潜在滑动面和对应的滑动体。

（2）建立考虑所分析边坡和滑动面的力学（应力）分析模型（如有限元、离散元、条分法等可以获得滑动面应力状态的方法均可），根据边坡的载荷、边界条件和岩土体物理力学参数，得到潜在滑动面的计算应力分布状态。

（3）将潜在滑动面分为若干段，从分析模型中获得每段编号、位置坐标和几何参数，从计算所得应力分布状态获得每滑面段的应力。

（4）根据滑动面位置，拟定试算的转动点位置区域，依次以区域内节点为转动点，根据式（8），（10）计算滑动面上的转动力矩和与抗转动力矩和，再按式（11）或（12）计算安全系数。

（5）对所试算区域的所有转动节点，通过计算找到最小安全系数，如果最小安全系数对应节点位于区域内部，可以认为此时的最小安全系数为所分析滑动面的基于力矩和法分析的安全系数，及其对应的转动节点。若最小安全系数对应节点位于区域边界，则回到第（4）步骤，调整试算的转动点位置区域。

2. 基于条分法模型的力矩和方法求解

力矩和法具有物理含义明确、概念清晰的特点，它可以直接采用显式计算，在已经确定了滑动面位置前提下，只要得到边坡在荷载和自重下的计算应力场，就可以通过力矩和法计算边坡的稳定性安全系数。为了获得滑动面上的应力，力矩和法是需要结合其它分析方法的。在边坡稳定性分析方法中，条分法是目前应用最广泛的分析方法之一。为了获得滑动面的应力，本文将采用条分法计算模型，通过条分法的计算获得滑动面上的应力，从而进行力矩和方法的计算。基于条分法的滑动面应力求解，其一便于方法之间的对比特别是与广泛应用的条分法的结果之间的对比分析；其二是更易于实现，避免了有限元计算时对模型和程序的较多依赖。

条分法中，Morgenstern-Price方法（M-P法）是一种同时满足力和力矩平衡条件的通用条分法^[6]，它是边坡极限平衡条分法^[7-9]中最具普遍性和假设条件的一种分析方法。下面将介绍M-P法与力矩和法相结合，求解安全系数。

如图 3所示，将滑体划分成n个竖直条块，取任意条块i，其中条块i的底面长度是 ${l_i}$ ，底边与X轴的夹角是 ${\alpha _i}$ ，转动点坐标是(x_p, y_p)，条块i底面左下脚坐标是(x_i1, y_i1)，条块i底面右下脚坐标是(x_i, y_i)。W_i是条块i的重力，E_i1，E_i分别是条块i左右两侧的水平条间力，X_i1，X_i分别是条块i左右两侧的条间剪力，h_i-1为作用在条块i左侧水平条间力E_i-1离基底的距离，h_i为作用在条块i右侧水平条间力E_i离基底的距离，H是坡高，β是坡角，P_i是底部的法向力, T_i是底部的切向剪力, u_i是底部的孔隙水压， ${\varphi '_i}$ 是滑动表面的有效内摩擦角，有效黏聚力为 ${c'_i}$ ，绕点(x_p, y_p)转动的安全系数F(x_p, y_p)。

图 3 条分法安全系数求解示意图

Figure 3. Schematic diagram of solving factor of safety by slices method

下载: 全尺寸图片幻灯片

M-P法假定条间力函数关系为 ${X_i} = \lambda f(x){E_i}$ 。式中 $\lambda$ 为任意常数， $f(x)$ 可取半正弦函数、常数等。 $f(x) =$ ${\text{sin}}\left[ {(x - a)/(b - a)} \right]{\text{π }}$ （正弦曲线），其中a，b分别为滑动面两端的x坐标。如果 $f(x)$ =1，则条间力( ${X_i}$ 和 ${E_i}$ 的合力)方向平行，即为斯宾塞法对条间力的假设。当力与力矩的平衡被严格满足时，在一个合理的范围内，对条间力的不同假设对其安全系数的影响不大^[10]。

根据莫尔-库仑准则及安全系数F的定义，建立条块的水平和垂直方向力平衡、滑体整体水平向力平衡以及将作用在条块上的力对土条条块底中点取矩建立力矩平衡。为了便于计算力矩和安全系数，得到 ${P_i}$ 和 ${T_i}$ 的表达式，有如下公式：

${E}_{i}={E}_{i-1}+ {P}_{i}\mathrm{sin}{\alpha }_{i}-{T}_{i}\mathrm{cos}{\alpha }_{i} \text{，}$

(13)

${P}_{i}=\frac{{W}_{i}-({X}_{i}-{X}_{i-1})-\frac{1}{F}({{c}^{\prime }}_{i}{l}_{i}-{u}_{i}{l}_{i}\mathrm{tan}{{\varphi }^{\prime }}_{i})}{\mathrm{cos}{\alpha }_{i}+\frac{1}{F}(\mathrm{sin}{\alpha }_{i}\mathrm{tan}{{\varphi }^{\prime }}_{i})} \text{，}$

(14)

${T}_{i}=\frac{{{c}^{\prime }}_{i}{l}_{i}+({P}_{i}-{u}_{i}{l}_{i})\mathrm{tan}{{\varphi}^{\prime }}_{i}}{F} \text{，}$

(15)

$\sum {({P_i}\sin {\alpha _i} - {T_i}\cos {\alpha _i})} = 0 \text{，}$

(16)

${E_{i - 1}}\left( {{h_{i - 1}} + {y_{i - 1}} - \frac{{{y_{i - 1}} + {y_i}}}{2}} \right) - {E_i}({h_i} - \frac{{{y_{i - 1}} + {y_i}}}{2} + {y_i}) -\\ ~~~~~~~~~~ {X}_{i-1}\left(\frac{{x}_{i-1}+{x}_{i}}{2}-{x}_{i-1}\right)-{X}_{i}({x}_{i}-\frac{{x}_{i-1}+{x}_{i}}{2})\text{=}0\text{ }。$

(17)

岩土体的抗剪强度 ${\tau _{{\text{f}}i}}$ 提供抗滑力，其值由莫尔库仑强度准则进行计算，滑动面上条块i的 ${\tau _{{\text{f}}i}}$ 的斜率为 ${k_{\tau {\text{f}}i}}$ ，所在直线方程为

$y={k}_{\tau \text{f}i}x+{b}_{i} \text{，}$

(18)

式中： ${k_{\tau {\text{f}}i}} = \frac{{{y_i} - {y_{i - 1}}}}{{{x_i} - {x_{i - 1}}}}$ ； ${b_i} = {y_i} - \frac{{{y_i} - {y_{i - 1}}}}{{{x_i} - {x_{i - 1}}}}{x_i}$ 。

转动点(x_p, y_p)到 ${\tau _{{\text{f}}i}}$ 所在直线的距离为

${d}_{\tau \text{f}i}=\frac{\mid {k}_{\tau \text{f}i}{x}_{\text{p}}+(-1){y}_{\text{p}}+{b}_{i}\mid }{\sqrt{{k}_{\tau \text{f}i}{}^{2}+{(-1)}^{2}}} 。$

(19)

抗剪强度 ${\tau _{{\text{f}}i}}$ 提供抗转动力矩，之和为

$\sum {\tau _{{\text{f}}i}}{l_i}{d_{\tau {\text{f}}i}} = \sum [{c'_i}{l_i} + ({P_i} - {u_i}{l_i})\tan {\varphi '_i}]{d_{\tau {\text{f}}i}} 。$

(20)

式中： ${l_i}$ 为条块i的基底长度。

滑动面上条块i的基底反力 ${P}_{i}$ 的斜率为 ${k_{Pi}}$ ，所在直线方程为

$y={k}_{Pi}x+{B}_{i} 。$

(21)

式中： ${k_{Pi}} = - \frac{1}{{{k_{\tau {\text{f}}i}}}} = - \frac{{{x_i} - {x_{i - 1}}}}{{{y_i} - {y_{i - 1}}}}$ ； ${B_i} = \frac{{{y_i} + {y_{{\text{i}} - 1}}}}{2} - {k_{Pi}} \cdot$ $\frac{{{x_i} + {x_{i - 1}}}}{2}$ 。

转动点(x_p, y_p)到 ${P}_{i}$ 所在直线的距离为

${d}_{Pi}=\frac{\mid {k}_{Pi}{x}_{\text{p}}+(-1){y}_{\text{p}}+{B}_{i}\mid }{\sqrt{{k}_{Pi}{}^{2}+{(-1)}^{2}}} \text{，}$

(22)

${P}_{i}$ 提供的抗转动力矩和为 $\sum {{P_i}{d_{Pi}}}$ 。

抗转动力矩和即为

$\sum {M}_{\text{r}}\text{(}{x}_{\text{p}},{y}_{\text{p}}\text{)}=\sum {\tau }_{\text{f}i}{l}_{i}{d}_{\tau \text{f}i}+\sum {P}_{i}{d}_{Pi}。$

(23)

${W_i}$ 提供转动力矩，之和为 $\sum {{W_i}{d_{Wi}}}$ 。其中 ${d_{Wi}} = {x_{\text{p}}} -$ ${x_i}$ 。

力矩和安全系数表达式如下：

$F{\text{(}}{x_{\text{p}}},{y_{\text{p}}}{\text{)}} = \frac{{\sum {M_{\text{r}}}{\text{(}}{x_{\text{p}}},{y_{\text{p}}}{\text{)}}}}{{\sum {M_{\text{s}}}{\text{(}}{x_{\text{p}}},{y_{\text{p}}}{\text{)}}}} = \frac{{\sum {\tau _{{\text{f}}i}}{l_i}{d_{\tau fi}} + \sum {P_i}{d_{Pi}}}}{{\sum {W_i}{d_{Wi}}}}。$

(24)

3. 算例

3.1 圆弧滑动面边坡算例

选用ACADS考题1中EX1(c)，此边坡为非均质边坡，由3层土组成，其计算模型见图 4，材料参数见表 1。

图 4 ACADS考题EX1(c)计算模型

Figure 4. Model of problem of EX1(c) of ACADS

下载: 全尺寸图片幻灯片

表 1 ACADS考题EX1(c)材料参数

Table 1. Material parameters of EX1(c) of ACADS

土层	c/kPa	φ/(°)	γ/(kN·m^-3)
1	0	38	19.5
2	5.3	23	19.5
3	7.2	20	19.5

下载: 导出CSV

| 显示表格

本文暂时不对边坡稳定力矩和法的临界滑动面的搜索问题进行讨论，图 4中的滑动面为边坡稳定性分析软件slide求得的临界滑动面，本算例用此滑动面计算边坡稳定的力矩和法安全系数，其中考题EX1(c)的滑动面同时穿过3种土层。

表 2为条分法计算的条块数据。表 3为多种方法计算结果。计算均采用同一个滑动面，从表中可以看出，本文方法计算得到的安全系数最小，其值为1.043，且对于圆弧滑动面，当转动点为圆心时，本文方法计算结果为1.378，与裁判答案基本一致，误差小于2%。

表 2 条方法计算过程

Table 2. Calculation process by slices method

条块i	基底长l_i/ m	条重W_i/ (kN·m^-3)	抗滑力 ${\tau _{{\text{f}}i}}{l_i}$ /(kN·m^-1)	基底反力P_i/ (kN·m^-1)
1	5.23	110.11	75.25	96.31
2	3.97	232.52	143.43	183.58
3	3.30	261.32	110.98	220.29
4	3.35	284.60	119.98	263.44
5	3.19	243.01	117.80	260.58
6	3.13	167.48	95.96	201.65
7	3.17	59.87	50.78	76.88

下载: 导出CSV

| 显示表格

表 3 安全系数计算结果

Table 3. Calculation results of safety factor

计算方法	圆心坐标/m	转动点坐标/m	安全系数
费伦纽斯法	(35.77, 23.54)		1.263
毕肖普法	(35.77, 23.54)		1.372
简布法	(35.77, 23.54)		1.296
M-P法	(35.77, 23.54)		1.379
矢量和法			1.375
力矩和法		(33.5, 4.6)	1.043

下载: 导出CSV

| 显示表格

3.2 非圆弧滑面边坡算例

非圆弧滑面边坡算例选用ACADS提供的考题3中EX3(b)，EX3(b)是一个含软弱夹层的非均质边坡考核题目，边坡计算模型见图 5，岩土参数见表 4。

图 5 ACADS考题EX3(b)计算模型

Figure 5. Model of problem of EX3(b) of ACADS

下载: 全尺寸图片幻灯片

表 4 ACADS考题EX3(b)材料参数

Table 4. Material parameters of EX3(b) of ACADS

土层	c/kPa	φ/(°)	γ/(kN·m^-3)
1	28.5	20	18.84
2	0	10	18.84

下载: 导出CSV

| 显示表格

对于具有软弱夹层的边坡，其滑面必定有一段会沿软弱夹层滑动，此时若依然采用圆弧滑动面，得到的结果往往是不正确的，故这里采用指定的折线滑动面进行计算边坡稳定的力矩和法安全系数。图 5给出剖面几何特征和指定滑动面位置，将滑动体划分为7个条块，滑动面同时穿过2种土层，底部一段沿软弱夹层滑动，滑动面控制点坐标见表 5。

表 5 滑动面控制点坐标

Table 5. Coordinates of control points for sliding surface

点号	X/m	Y/m
1	10.69	20.00
2	20.50	7.00
3	40.00	6.50
4	42.15	7.75

下载: 导出CSV

| 显示表格

表 6为条分法计算的条块数据。表 7为多种方法计算结果。

表 6 条分法计算过程

Table 6. Calculation process by slices method

条块i	基底长l_i/m	条重W_i/ (kN·m^-3)	抗滑力 ${\tau _{{\text{f}}i}}{l_i}$ / (kN·m^-1)	基底反力P_i/ (kN·m^-1)
1	8.14	300.33	299.76	185.96
2	8.14	825.64	491.09	711.65
3	6.50	1158.27	207.78	1178.36
4	6.50	780.68	144.07	817.07
5	6.50	403.10	75.24	426.69
6	0.99	20.82	5.24	29.71
7	1.49	9.21	56.75	39.08

下载: 导出CSV

| 显示表格

表 7 安全系数计算结果

Table 7. Calculated results of factor of safety

计算方法		转动点坐标	安全系数
简布法			1.195
陆军工程师团法			1.364
M-P法			1.286
矢量和法			1.221
力矩和法	条间力函数为常数	(40, 6.5)	1.051
力矩和法	条间力函数为半正弦函数	(40, 6.5)	1.049

下载: 导出CSV

| 显示表格

从表 7可以看出，在非圆弧滑面的情况下，本文方法得到的安全系数较小，而不同的条间力假设对力矩和安全系数的影响并不明显，而对于含软弱夹层的非均质边坡，滑面偏离圆弧较大，本文方法的优越性是能够处理各种非圆弧滑面。

4. 结论

（1）提出了分析边坡稳定性的新方法即力矩和方法。该方法对一确定的滑坡体和滑动面，假定滑坡体绕某点转动，求出滑坡体在自重和外荷载作用下的转动力矩和，同时求出滑动面上基岩法向反力和剪切强度对滑动体的抗转动力矩和，定义抗转动力矩和与转动力矩和的比值为对应的安全系数。本文中，转动点由最小安全系数决定，并通过试算搜索得到。

（2）力矩和方法是对矢量和方法的一个补充。滑动体视为刚体时，它的位移应该包括刚体的平动和转动两部分（只有外力合力通过重心时没有外力引起的转动），矢量和方法认为滑动体是沿某一直线方向平动，不考虑其转动，而力矩和方法考虑了转动，两个方法互为补充。当然，若能同时考虑滑体的平动和转动，将会更合理，也是本文后续将延续的工作。

（3）本文在求解滑动面应力时，采用了更简便的条分法计算（M-P法）。并结合M-P法模型，推导了力矩和方法求解公式。

（4）算例1结果显示，对于所分析的滑动面，由于所选择的转动点不同，力矩和方法所得到的安全系数较其它几种方法更小。这主要是因为最终所得最小安全系数对应的转动点不是所分析滑动面圆弧的圆心点。滑动面圆弧的圆心点只是力矩和方法所分析转动点的一个，其对应的安全系数不一定是滑动体绕某点发生转动时的最小安全系数。

这一结果也提示一点，具有圆弧滑动面的滑坡体，其初始失稳转动时所对应的转动点（或圆心）不一定是圆弧滑动面对应的圆心。进一步可认为，任意形状的滑动面，对于滑动体，滑动面只是其约束边界条件，与滑动体初始失稳或破坏的位移趋势无直接关系。这与人们习惯的认知是不同的，值得后续进一步的研究。

算例1中，当转动点取为圆心时，本文方法求得的安全系数与裁判答案基本一致，说明力矩和方法计算安全系数是可行的。

（5）力矩和方法可以分析任意滑动面。滑动体所对应的滑动面，只是滑动体的边界条件，其形状不会直接影响滑动体作为刚体自身的位移趋势，而是通过边界提供反力影响滑动体的后续位移轨迹。算例2给出了本文方法分析非圆弧滑动面的情况，并可以搜索到比原结果较小的安全系数，说明假定滑动体绕某点转动的位移趋势是值得对边坡进行稳定性分析时考虑。

图 1 模型箱示意图

Figure 1. Schematic diagram of model box

下载: 全尺寸图片幻灯片

图 2 孔隙水压力变化规律

Figure 2. Variation laws of pore water pressure

下载: 全尺寸图片幻灯片

图 3 土层位移变化规律

Figure 3. Variation laws of soil displacement

下载: 全尺寸图片幻灯片

图 4 双通道渗流示意图

Figure 4. Schematic diagram of dual-channel seepage

下载: 全尺寸图片幻灯片

图 5 基本工况模型示意图

Figure 5. Schematic diagram of model for basic working conditions

下载: 全尺寸图片幻灯片

图 6 降雨入渗模式转换示意图

Figure 6. Schematic diagram of rainfall infiltration mode conversion

下载: 全尺寸图片幻灯片

图 7 阶跃函数曲线

Figure 7. Curve of step function

下载: 全尺寸图片幻灯片

图 8 模型表面入渗率变化规律

Figure 8. Variation laws of surface infiltration rate of model

下载: 全尺寸图片幻灯片

图 9 土层饱和度分布

Figure 9. Distribution of soil saturation

下载: 全尺寸图片幻灯片

图 10 基本工况隧道拱顶沉降云图

Figure 10. Nephogram of vault settlement of tunnel under basic conditions

下载: 全尺寸图片幻灯片

图 11 基本工况隧道水平位移云图

Figure 11. Nephogram of horizontal displacement of tunnel under basic conditions

下载: 全尺寸图片幻灯片

图 12 基本工况隧道变形规律

Figure 12. Deformation laws of tunnel under basic condition

下载: 全尺寸图片幻灯片

图 13 工况1隧道变形规律

Figure 13. Deformation laws of tunnel under condition 1

下载: 全尺寸图片幻灯片

图 14 工况2隧道变形规律

Figure 14. Deformation laws of tunnel under condition 2

下载: 全尺寸图片幻灯片

图 15 模型表面入渗率变化规律

Figure 15. Variation laws of surface infiltration rate of model

下载: 全尺寸图片幻灯片

图 16 裂隙底部边界通量变化规律

Figure 16. Variation laws of flux at bottom boundary of fissures

下载: 全尺寸图片幻灯片

图 17 工况3隧道变形规律

Figure 17. Deformation laws of tunnel under working condition 3

下载: 全尺寸图片幻灯片

图 18 工况4隧道变形规律

Figure 18. Deformation laws of tunnel under working condition 4

下载: 全尺寸图片幻灯片

图 19 工况5隧道变形规律

Figure 19. Deformation laws of tunnel under working condition 5

下载: 全尺寸图片幻灯片

图 20 不同工况下的隧道拱顶沉降

Figure 20. Settlement of tunnel vault under different operating conditions

下载: 全尺寸图片幻灯片

图 21 不同工况下的隧道水平位移

Figure 21. Horizontal displacements of tunnel under different operating conditions

下载: 全尺寸图片幻灯片

表 1 各土层物理力学参数

Table 1 Physical and mechanical parameters of soil layers

土层	密度/(kg∙m^-3)	含水率/%	孔隙率	弹性模量/MPa	泊松比	黏聚力/kPa	摩擦角/(°)	饱和渗透系数/(m∙s^-1)
砂质土层	1380	17.9	0.412	18	0.3	21	27	2.1e^-5
黏性土层	1500	16.7	0.385	25	0.3	25	30	1.3e^-5

下载: 导出CSV

表 2 各区域非饱和水力参数

Table 2 Unsaturated hydraulic parameters of various regions

区域	模型	α/m^-1	n	l	θ_r
基质	VG	2.0	1.5	0.5	0.04
裂隙	BC	2.68	0.131	0.5	0.04

下载: 导出CSV

表 3 模型参数

Table 3 Model parameters

K_s/(m·d^-1)	v₀/(m·d^-1)	θ_s/%	θ_r/%	α/m^-1	n	H₀/m
1.39×10^-5	4K_s	40	4	2.5	2.1	-0.4

下载: 导出CSV

表 4 多工况分析中各参数变化值

Table 4 Change of various parameters in multi-operating condition analysis

工况	参数变化
1	裂隙位置X=30 m，裂隙宽度u=0.002 m，降雨强度V₀=0.1 m·d^-1
2	裂隙位置X=25 m，裂隙宽度u=0.002 m，降雨强度V₀=0.1 m·d^-1
3	裂隙位置X=20 m，裂隙宽度u=0.008 m，降雨强度V₀=0.1 m·d^-1
4	裂隙位置X=30 m，裂隙宽度u=0.008 m，降雨强度V₀=0.1 m·d^-1
5	裂隙位置X=25 m，裂隙宽度u=0.008 m，降雨强度V₀=0.1 m·d^-1

下载: 导出CSV

参考文献(21)

[1]	蔡国庆, 韩博文, 韦靖威, 等. 复杂水-力路径下非饱和砂质黄土增湿变形特性[J]. 岩石力学与工程学报, 2022, 41(增刊1): 3073-3080. CAI Guoqing, HAN Bowen, WEI Jingwei, et al. Wetting deformation characteristics of unsaturated sandy loess under complex hydro-mechanical paths[J]. Chinese Journal of Rock Mechanics and Engineering, 2022, 41(S1): 3073-3080. (in Chinese)
[2]	李锦辉, 郭凌波, 张利民. 考虑裂隙动态变化时裂隙土土水特征曲线的预测方法研究[J]. 岩石力学与工程学报, 2013, 32(6): 1288-1296. LI Jinhui, GUO Lingbo, ZHANG Limin. Prediction of swcc for unsaturated cracked soil considering development process of cracks[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(6): 1288-1296. (in Chinese)
[3]	VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(5): 892-898. doi: 10.2136/sssaj1980.03615995004400050002x
[4]	殷宗泽, 徐彬. 反映裂隙影响的膨胀土边坡稳定性分析[J]. 岩土工程学报, 2011, 33(3): 454-459. http://cge.nhri.cn/article/id/13962 YIN Zongze, XU Bin. Slope stability of expansive soil under fissure influence[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(3): 454-459. (in Chinese) http://cge.nhri.cn/article/id/13962
[5]	邵生俊, 杨春鸣, 焦阳阳, 等. 湿陷性黄土隧道的工程性质分析[J]. 岩土工程学报, 2013, 35(9): 1580-1590. http://cge.nhri.cn/article/id/15269 SHAO Shengjun, YANG Chunming, JIAO Yangyang, et al. Engineering properties of collapsible loess tunnel[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(9): 1580-1590. (in Chinese) http://cge.nhri.cn/article/id/15269
[6]	赖金星, 樊浩博, 来弘鹏, 等. 软弱黄土隧道变形规律现场测试与分析[J]. 岩土力学, 2015, 36(7): 2003-2012, 2020. LAI Jinxing, FAN Haobo, LAI Hongpeng, et al. In-situ monitoring and analysis of tunnel deformation law in weak loess[J]. Rock and Soil Mechanics, 2015, 36(7): 2003-2012, 2020. (in Chinese)
[7]	FRENELUS W, PENG H, ZHANG J Y. Seepage actions and their consequences on the support scheme of deep-buried tunnels constructed in soft rock strata[J]. Infrastructures, 2024, 9(1): 13. doi: 10.3390/infrastructures9010013
[8]	CHEN L L, WANG Z F, WANG Y Q. Failure analysis and treatments of tunnel entrance collapse due to sustained rainfall: a case study[J]. Water, 2022, 14(16): 2486. doi: 10.3390/w14162486
[9]	SI J L, LIU S P, ZHANG H J, et al. Failure investigation and treatments of tunnel entrance collapse in weak diatomaceous soil induced by heavy rainfall through coupling surface and groundwater flows[J]. Engineering Failure Analysis, 2023, 150: 107337. doi: 10.1016/j.engfailanal.2023.107337
[10]	韩同春, 何露, 林博文, 等. 一种裂隙土的双重入渗模型及对边坡稳定的影响[J]. 华南理工大学学报(自然科学版), 2019, 47(5): 123-129, 138. HAN Tongchun, HE Lu, LIN Bowen, et al. Double infiltration model of fractured soil and its influence on slope stability[J]. Journal of South China University of Technology (Natural Science Edition), 2019, 47(5): 123-129, 138. (in Chinese)
[11]	开迪尔丁·吾拉木, 张紫昭, 张艳阳, 等. 基于COMSOL Multiphysics的降雨型滑坡临界降雨阈值计算模型研究: 以新疆新源县喀拉海依苏滑坡隐患体为例[J]. 工程地质学报, 2023, 31(4): 1364-1374. KAIDIERDING Wuramu, ZHANG Zizhao, ZHANG Yanyang, et al. Comsol multiphysic based calculation model of critical rainfall threshold for rainfall-induced landslide: a case study of karahayisu landslide in Xinyuan County, Xinjiang[J]. Journal of Engineering Geology, 2023, 31(4): 1364-1374. (in Chinese)
[12]	曾铃, 史振宁, 付宏渊, 等. 降雨入渗对边坡暂态饱和区分布特征的影响[J]. 中国公路学报, 2017, 30(1): 25-34. ZENG Ling, SHI Zhenning, FU Hongyuan, et al. Influence of rainfall infiltration on distribution characteristics of slope transient saturated zone[J]. China Journal of Highway and Transport, 2017, 30(1): 25-34. (in Chinese)
[13]	邱军领, 秦祎文, 赖金星, 等. 突发高压渗流作用下黄土地铁隧道水压阶跃效应分析[J]. 岩土工程学报, 2023, 45(4): 758-767. doi: 10.11779/CJGE20220062 QIU Junling, QIN Yiwen, LAI Jinxing, et al. Step effects of hydraulic pressure of metro tunnels in loess under sudden high-pressure seepage[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(4): 758-767. (in Chinese) doi: 10.11779/CJGE20220062
[14]	张万志, 徐帮树, 曾仲毅, 等. 降雨入渗下膨胀性黄土隧道围岩破坏演化[J]. 东南大学学报(自然科学版), 2018, 48(4): 736-744. ZHANG Wanzhi, XU Bangshu, ZENG Zhongyi, et al. Research on failure evolution process of surrounding rock of swelling loess tunnel under rainfall infiltration[J]. Journal of Southeast University (Natural Science Edition), 2018, 48(4): 736-744. (in Chinese)
[15]	高阳, 孙浩凯, 刘德军, 等. 强降雨影响下破碎复理岩地层隧道洞口段失稳机理[J]. 中南大学学报(自然科学版), 2019, 50(9): 2295-2303. GAO Yang, SUN Haokai, LIU Dejun, et al. Collapse mechanism of tunnel portal sectionin broken flysch under influence of heavy rainfall[J]. Journal of Central South University (Science and Technology), 2019, 50(9): 2295-2303. (in Chinese)
[16]	谢强, 陈昱成, 傅翔, 等. 非饱和瞬态渗流的DDA流固耦合模型研究[J]. 岩土工程学报, 2024, 46(2): 299-306. doi: 10.11779/CJGE20221026 XIE Qiang, CHEN Yucheng, FU Xiang, et al. Fluid-solid coupling model for discontinuous deformation analysis of unsaturated transient seepage[J]. Chinese Journal of Geotechnical Engineering, 2024, 46(2): 299-306. (in Chinese) doi: 10.11779/CJGE20221026
[17]	黄阜, 纪恒博, 王子钦, 等. 考虑土体非均质性和孔隙水压力耦合作用的地下连续墙槽壁稳定性研究[J]. 岩土工程学报, 2024, 46(3): 539-548. doi: 10.11779/CJGE20221401 HUANG Fu, JI Hengbo, WANG Ziqin, et al. Stability of slurry trench walls considering coupling effects of soil heterogeneity and pore water pressure[J]. Chinese Journal of Geotechnical Engineering, 2024, 46(3): 539-548. (in Chinese) doi: 10.11779/CJGE20221401
[18]	SHAO W, BOGAARD T, BAKKER M. How to use COMSOL multiphysics for coupled dual-permeability hydrological and slope stability modeling[J]. Procedia Earth and Planetary Science, 2014, 9: 83-90. doi: 10.1016/j.proeps.2014.06.018
[19]	马慧. COMSOL Multiphysics基本操作指南和常见问题解答[M]. 人民交通出版社, 2009. MA Hui. COMSOL Multiphysics Basic Operation Guide and Frequently Asked Questions[M]. Beijing: China Communications Press, 2009. (in Chinese)
[20]	年庚乾, 陈忠辉, 张凌凡, 等. 边坡降雨入渗问题中两种边界条件的处理及应用[J]. 岩土力学, 2020, 41(12): 4105-4115. NIAN Gengqian, CHEN Zhonghui, ZHANG Lingfan, et al. Treatment of two boundary conditions for rainfall infiltration in slope and its application[J]. Rock and Soil Mechanics, 2020, 41(12): 4105-4115. (in Chinese)
[21]	侯晓萍, 樊恒辉. 基于COMSOL Multiphysics的非饱和裂隙土降雨入渗特性研究[J]. 岩土力学, 2022, 43(2): 563-572. HOU Xiaoping, FAN Henghui. Study on rainfall infiltration characteristics of unsaturated fractured soil based on COMSOL Multiphysics[J]. Rock and Soil Mechanics, 2022, 43(2): 563-572. (in Chinese)