0.   引言

含裂隙土体在自然界及实际工程中普遍存在。土体中发育的裂隙是许多工程灾害的诱因，例如隧道失稳、边坡滑塌、水土流失等^[1]。而在降雨条件下，裂隙的存在为雨水入渗提供了优势通道^[2]，对于隧道工程来说，降雨对于含裂隙土层隧道的不利影响主要体现在两方面：一是降雨导致渗流场变化，使得更多的土体转向饱和状态，从而引起作用在土体上的水荷载增大和土体抗剪强度降低；二是在持续降雨作用下，裂隙中的渗流速率明显高于基质中，一旦入渗到隧道影响范围内，便会对隧道稳定性产生影响^[3-6]。

目前有许多学者针对降雨入渗条件下的隧道或边坡稳定性进行了研究。Frenelus等^[7]根据非线性弹塑性理论，提出了渗流条件下软岩隧道周围应力、位移和塑性半径的解析解，同时基于应变软化模型，研究了隧道围岩的流固耦合特性。Chen等^[8]通过实地调查、理论分析和现场监测等手段，分析了降雨条件下导致云南省某隧道塌方的潜在因素，结果表明降雨会降低土体剪切强度，导致隧道坍塌。Si等^[9]通过建立地表径流、地下水流和土力学的耦合模型，并应用到隧道工程中，提出了一种评估降雨条件下隧道洞口稳定性的可靠方法；韩同春等^[10]基于伽辽金法得到了理查兹方程离散化的解，并应用到边坡工程中分析了边坡稳定性；开迪尔丁·吾拉木等^[11]基于非饱和渗流理论建立了滑坡稳定性计算模型，提供了一种计算降雨型黄土滑坡降雨阈值方法；曾铃等^[12]基于有限元模拟开展了一维渗流及二维渗流计算，得到了边坡内部的渗流场特征，并研究了多因素影响下的渗流场变化特征；邱军领等^[13]通过模型试验和数值模拟相对比，提出并分析了突发高压渗流作用下的黄土地铁隧道水压阶跃效应；张万志等^[14]通过FLAC的二次开发，分析了降雨条件下膨胀性黄土隧道围岩的破坏演化；高阳等^[15]通过离散元软件，分析了强降雨条件下，岩质隧道洞口段围岩裂隙及渗流的变化特征；谢强等^[16]建立了DDA瞬态整体渗流矩阵，并通过相关工程案例进行了数值模拟验证了该模型的正确性，为处理非饱和瞬态流固耦合问题提供了新思路；黄阜等^[17]基于三维离散型破坏机制和极限分析上限定理，提出了一种考虑土体非均质性和孔隙水压力耦合作用下的算法，并通过数值模拟进行了验证。上述研究对流固耦合的模拟可分为3类：①含裂隙土层的边坡降雨入渗；②降雨入渗条件下的隧道稳定性分析；③流固耦合的数值模型研究。而目前对降雨入渗的模拟大多只设定固定流量或固定压力的边界，没有考虑边界转换的过程，同时考虑裂隙土层中不同裂隙形式对隧道稳定性的影响研究也亟需补充。

鉴于此，本文以含单一裂隙土层局部渗流试验中得到的现象为引，提出双通道渗流理论，并基于阶跃函数法，提出降雨诱发渗流的边界条件，同时讨论了阶跃函数过渡区长度的合理设置。应用有限元软件COMSOL Multiphysics对含裂隙土体的非饱和渗流进行数值模拟，并以裂隙为切入点，探究降雨入渗作用下不同裂隙形式对隧道变形特性的影响。

1.   含单一裂隙土层局部渗流性态分析

1.1   试验方案

模型试验装置由模型箱、进水系统、排水系统、传感器和孔隙水压力计组成。模型箱尺寸为1.5 m×0.9 m×1.5 m，外部为钢架结构，内部为透明钢化玻璃板制成的箱型结构（见图 1）。模型箱上部设置进水管，底部设置排水系统。模型试验中地层为黄土地层，取自西安地铁8号线长鸣路车辆段。假定土层中某一点为出水点，该出水点下方有纵向裂隙。以模型箱侧面为主视角，沿纵向在出水点前后等距离选取两个断面，每个断面距出水点50 cm，每一横断面上布设6条传感器，选择裂隙正下方同一断面上的3点作为土层位移测点，通过传感器监测进水后的土层位移。在出水点同一横断面内，于出水点下方10 cm裂隙处和裂隙左右两侧等距各布设1个孔隙水压力测点，共3个测点，测点与测点之间距离15 cm。

图  1  模型箱示意图

Figure  1.  Schematic diagram of model box

下载: 全尺寸图片 幻灯片

1.2   试验结果分析与讨论

由图 2可知，测点1的孔隙水压力值要明显高于测点2和测点3。在进水后的4 h内，3个测点的孔隙水压力值都迅速增加，进入快速增长阶段。测点1的增速要明显高于测点2和测点3，测点1的孔隙水压力最大值达到140 kPa，远超测点2和测点3的最大值70 kPa。在进水后4~20 h内，3个测点的孔隙水压力值缓慢增长，进入缓慢增长阶段。进水后21 h裂隙被水不断侵蚀扩张，孔隙水压力值进入水压突增阶段，产生阶段式的突变。土层位移图与孔隙水压力图的规律大致相同，但只有裂隙下方土体的变形有明显的阶段性变化（见图 3）。随着进水时间的增加，处在裂隙下方土层的位移量明显高于裂隙两侧的土层，经历变形快速增长阶段后进入缓慢增长阶段。在时间达到21 h后，裂隙下方土层位移值也产生突变现象。

图  2  孔隙水压力变化规律

Figure  2.  Variation laws of pore water pressure

下载: 全尺寸图片 幻灯片

图  3  土层位移变化规律

Figure  3.  Variation laws of soil displacement

下载: 全尺寸图片 幻灯片

随着进水时间的增加，水分不仅会沿着裂隙产生垂直渗流，还会沿着裂隙的法向产生水平方向的渗流，导致裂隙水平截面上的孔隙水压力增大、土层的位移增大。虽然试验假定了土体中只含有单一裂隙，但实际渗流过程中，水分不仅仅只沿裂隙进行单一通道入渗，还会沿着裂隙的法向水平扩展，因此将该现象称作双通道渗流现象，即裂隙-基质渗流现象。

2.   含裂隙土层双通道渗流模型构建

土体裂隙中的渗流速度要明显高于基质中，这就是所谓的优先流现象。土体中的优先流通常使用双通道渗流模型来描述^[18]，如图 4所示。双通道渗流理论认为土体由基质域和大孔隙域组成，发生渗流时还会在两域之间产生水分交换。

图  4  双通道渗流示意图

Figure  4.  Schematic diagram of dual-channel seepage

下载: 全尺寸图片 幻灯片

双通道渗流模型的控制方程如下：

其中，

式中：${C_{\text{f}}}$，${C_{\text{m}}}$为大孔隙域、基质域对应的比容水度；${S_{\text{e}}}$为有效饱和度；${S_{\text{s}}}$为储水率；${h_{\text{f}}}$，${h_{\text{m}}}$为大孔隙域、基质域对应的压力水头；$\nabla$为微分算子；${K_{\text{f}}}$，${K_{\text{m}}}$为大孔隙域、基质域对应的不饱和导水率；${\varGamma _{\rm{\mathsf{ω}}}}$为水交换项；${\alpha _{{\rm{\mathsf{ω}}}}}$为一阶传质系数；$\beta$为土壤团粒几何形状因子，球形一般取15，矩形一般取3；$d$为团粒块中心到孔隙边缘的距离，通常取0.01 m；${\varUpsilon _{{\rm{\mathsf{ω}}} }}$为经验参数，一般取0.4；${K_{\text{a}}}$为大孔隙域与基质域界面处的不饱和导水率；${\omega _{\text{f}}}$，${\omega _{\text{m}}}$分别为大孔隙域与基质域占土体的体积分数。

Richards方程可用于描述双通道渗流模型中的大孔隙域与基质域：

将$\frac{{\partial {h_{\text{p}}}}}{{\partial t}}$展开，得到

代入式（6）中，得到关于孔隙水压力$p$的Richards方程：

“裂隙流”理论定义了大孔隙域的边界^[19]，使流体能够沿着大孔隙域边界的法向流动，可近似模拟大孔隙域与基质域的水分交换（如图 5所示）。裂隙流的控制方程为

图  5  基本工况模型示意图

Figure  5.  Schematic diagram of model for basic working conditions

下载: 全尺寸图片 幻灯片

式中：${q_{\text{f}}}$为裂隙中每单位长度的体积流量；${K_{\text{f}}}$为渗透率；$\mu$为流体的动力黏度；${d_{\text{f}}}$为孔径；${\nabla _{\text{T}}}$为限制在裂隙切向平面上的梯度算子；$p$为压力；$\rho$为流体密度；$g$为重力加速度；$D$为纵坐标。其物理意义是使用切向导数来定义沿边界的流动。

将裂隙流方程与土体的材料特性联系起来，会产生一个新的压力方程：

式中：${\varepsilon _{\text{f}}}$为孔隙率；${Q_{\text{m}}}$为质量源项。

3.   含裂隙土质隧道的降雨入渗数值模拟

3.1   数值模型

本文使用COMSOL Multiphysics多场耦合软件来模拟降雨入渗条件下含裂隙土层中的优先流现象，探究隧道变形特征。模型的长度为50 m，高度为40 m。模型共设置两个土层，从地表向下分别为砂性土层（7 m）和黏性土层（33 m）。在模型中部（x=20 m处）存在一条垂直裂隙，宽度为0.002 m，贯穿整个模型。地下水位设置为25 m处，降雨强度取0.05 m/d。土层的本构模型采用Drucker-Prager（DP）准则，基质的非饱和水力参数由Van-Genuchten（VG）方程确定，裂隙的非饱和水力参数由Brooks-Corey（BC）方程确定。各土层物理力学参数及非饱和水力参数见表 1，2。

表  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.2   基于阶跃函数的降雨入渗边界模拟

含裂隙土的降雨入渗模拟有两个关键点：一是流量边界和压力边界的模拟（见图 6）。降雨初期，土体入渗能力强，此时为流量边界；当降雨强度超过土体的饱和导水率时，土体入渗能力不断降低，土体表面会覆盖一层厚度较小的水膜，边界条件转换为压力边界，两种边界方程如下所示：

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

Figure  6.  Schematic diagram of rainfall infiltration mode conversion

下载: 全尺寸图片 幻灯片

式中：$n$为边界的外法线向量；${v_0}$为降雨强度；${R_{\text{b}}}$为流体从空气单元到基体的传导率；${D_{\text{b}}}$为外部高程；$D$为纵坐标。

二是混合入渗的模拟（见图 6）。土体表面孔隙水压力值由正转负的过程即为混合入渗。引入${c_1}$，${c_2}$作为降雨入渗边界条件的控制系数：

将控制系数与两种边界方程结合起来，可以得到降雨入渗边界条件的转换方程：

当两种边界进行转换时，需要对控制系数设置阶跃函数，如图 7所示。在阶跃函数设置过渡区不仅可以近似模拟混合入渗模式，还能保证计算的收敛性。既有文献中尚还缺乏对于阶跃函数过渡区的讨论，本文通过改变过渡区长度讨论其合理的取值。

图  7  阶跃函数曲线

Figure  7.  Curve of step function

下载: 全尺寸图片 幻灯片

阶跃函数的x轴为孔隙水压力，y轴为条件控制系数，将孔隙水压力零点附近范围内设为过渡区，为使过渡区函数形式平滑，采用二阶连续可导函数，改进后的控制系数如下：

式中：$p = \left( {\frac{{ - 2}}{{{L^2}}}} \right){p^2} - \left( {\frac{2}{L}} \right)p + \frac{1}{2}$，$L$为过渡区长度。

将过渡区长度$L$分别设置为2，4，8，16，32个孔隙水压力单位长度，部分模型参数参考年庚乾等^[20]建立的降雨入渗模型，各参数取值见表 3。

表  3  模型参数

Table  3.  Model parameters

Ks/(m·d-1)	v0/(m·d-1)	θs/%	θr/%	α/m-1	n	H0/m
1.39×10-5	4Ks	40	4	2.5	2.1	-0.4

下载: 导出CSV 

| 显示表格

图 8为不同过渡区长度情况下，模型表面入渗率曲线图。本文数值模拟中应用了“空气单元法^[21]”，当入渗边界为流量边界时，入渗率峰值变化体现在了过渡区长度的变化上；当入渗边界为压力边界时，积水层囊括了空气单元，导致计算结果与文献中略有不同，但规律大致相同。在降雨初期，即流量边界情况下，当过渡区长度L取2，16，32个单位时，入渗率峰值偏大，均超过7.5×10^-5 m/s；过渡区长度取4，8个单位时，入渗率峰值与文献中接近，约为6.9×10^-5 m/s。在降雨中期和后期，即压力边界情况下，入渗率结果无明显变化，均趋于稳定值。为了保证计算结果的准确性，后续将过渡区长度统一设为4个单位。

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

Figure  8.  Variation laws of surface infiltration rate of model

下载: 全尺寸图片 幻灯片

3.3   多因素隧道变形特征影响分析

将模型的裂隙形式和降雨强度作变化，进行多因素影响性分析。按裂隙与隧道之间的相对位置划分为3种工况，分为裂隙在隧道左侧（X=20 m）、隧道右侧（X=30 m）和隧道上方（X=25 m），然后改变降雨强度和裂隙宽度。为统一变量，多工况分析中降雨强度统一取基本工况中的2倍，裂隙宽度取基本工况中的4倍，共设5种工况（见表 4）。

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

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

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

下载: 导出CSV 

| 显示表格

4.   隧道变形特征分析与讨论

4.1   基本工况

图 9为基本工况下的土体饱和度随时间变化云图。可以看出，随着降雨时间增长，土体表层饱和度不断增大，降雨持续12 h后土体中出现优先流现象，裂隙附近饱和度明显高于基质，且产生裂隙与基质的水分交换，符合双通道渗流模型的规律。

图  9  土层饱和度分布

Figure  9.  Distribution of soil saturation

下载: 全尺寸图片 幻灯片

基本工况中隧道拱顶沉降和上部水平位移规律如图 10~12所示。其中隧道水平位移负值表示沿水平向左的位移，正值代表沿水平向右的位移。由图 12（a）可以看出，降雨持续过程中，隧道拱顶沉降逐渐增大，且拱顶沉降曲线逐渐向近裂隙一侧偏移，沉降最大值达到306 mm。由图 12（b）可知，隧道上部水平位移先达到稳定值，后由于裂隙优先流现象，隧道上部产生的最大水平位移超过了20 mm。随着降雨持续时间增长，位移曲线不断向裂隙处扩展，由大致对称的斜“S”型曲线逐渐变为斜“C”型曲线。

图  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

下载: 全尺寸图片 幻灯片

4.2   降雨强度影响性分析

图 13，14为工况1和2的隧道拱顶沉降和水平位移，工况1和2在增加了一倍降雨强度数值的情况下，拱顶沉降最大值与基本工况相比变化不大，上部水平位移的变化则与裂隙与隧道间相对位置有关。当裂隙在隧道两侧时，隧道近裂隙一侧的水平位移变化较大，最大超过20 mm；当裂隙在隧道上方时，隧道水平位移变化较小，最大不超过4 mm。

图  13  工况1隧道变形规律

Figure  13.  Deformation laws of tunnel under condition 1

下载: 全尺寸图片 幻灯片

图  14  工况2隧道变形规律

Figure  14.  Deformation laws of tunnel under condition 2

下载: 全尺寸图片 幻灯片

针对该现象，对该模型入渗边界的入渗率和隧道上方的裂隙底部的边界通量进行分析，图 15为模型降雨强度分别为0.05，0.1，0.2 m/d的入渗率结果图，图 16为裂隙底部边界通量结果图。可以看出，随着降雨强度的增加，模型表面的入渗率在降雨初期不断增大，最大值分别达到了2.5×10^-5，3.2×10^-5，3.3×10^-5 m/s，之后随着降雨时间的增加，3种工况表面入渗率逐渐降低，每种工况入渗率均达到稳定值0.7×10^-5 m/s。而降雨强度为0.05，0.1 m/d时，裂隙底部的边界通量变化不大，均小于0.0003 kg/(m²·s)，而当降雨强度达到0.2 m/d时，边界通量在降雨48h后逐渐增大，最大值达到0.009 kg/(m²·s)，后在降雨63 h后逐渐减小。

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

Figure  15.  Variation laws of surface infiltration rate of model

下载: 全尺寸图片 幻灯片

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

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

下载: 全尺寸图片 幻灯片

结果表明，随着降雨时间增长，模型表面入渗率变化规律相近，只是最大入渗率结果略有差别，裂隙底部边界通量在极端降雨条件下（V₀=0.2 m/d）时会产生突变。

4.3   裂隙宽度影响性分析

图 17~19分别为工况3~5的隧道拱顶沉降和水平位移结果图。裂隙在隧道两侧时，隧道拱顶沉降均会向裂隙处扩展，最大值产生在降雨72 h时，分别为319，320 mm；上部水平位移曲线则会向近裂隙一侧发展，最大水平位移均产生在降雨72 h时，分别为25，28 mm。当裂隙位置在隧道上方时，裂隙宽度对隧道拱顶沉降的影响更为明显，拱顶沉降最大值达到了328 mm。上部水平位移变化规律与基本工况大致相同，最大值为26 mm。

图  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

下载: 全尺寸图片 幻灯片

由上述分析可知，最大拱顶沉降和水平位移均产生在降雨持续72 h时，因此根据图 20，21纵向比较降雨持续72 h时各工况下隧道变形规律。由图 20可知，裂隙宽度的增大会导致拱顶沉降显著增大。裂隙在隧道两侧时，8 mm裂隙宽度下隧道拱顶沉降比2 mm裂隙宽度下增大约3.9%，裂隙在隧道上方时，该数值增大到6.5%。由图 21可知，不同裂隙位置和裂隙宽度都会显著影响隧道的水平位移。裂隙在隧道左右两侧时，8 mm裂隙宽度下隧道最大水平位移较2 mm裂隙宽度下分别增大8.7%和12%，裂隙在隧道上方时，该数值达到了44.4%。

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

Figure  20.  Settlement of tunnel vault under different operating conditions

下载: 全尺寸图片 幻灯片

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

Figure  21.  Horizontal displacements of tunnel under different operating conditions

下载: 全尺寸图片 幻灯片

结果表明，随着降雨时间的增长，裂隙位置的变化和裂隙宽度的增大会显著加剧隧道沉降值和上部水平位移值的增速。裂隙在隧道上方时，比基本工况大4倍宽度的裂隙会使隧道上部水平位移值增大超40%，严重影响隧道的稳定性。

5.   结论

根据含单一裂隙土层的局部渗流性态表征，提出含裂隙土层的双通道渗流理论，并基于阶跃函数提出改进的降雨入渗边界条件方程，讨论过渡区长度并与既有文献中的结果进行对比验证。以隧道工程为例，采用COMSOL Multiphysics有限元软件对降雨条件下裂隙土层隧道的渗流场和围岩变形特性进行数值模拟研究，相关结果可为复杂地质情况下的隧道变形预测作为参考。

（1）以双渗透模型理论和裂隙流理论作为入渗模型的基础，提出基于阶跃函数的降雨入渗边界条件方程，以既有文献中的降雨入渗模型为基础，讨论了阶跃函数过渡区的合理设置，发现过渡区长度为4或8个单位时，计算结果较准确。

（2）将降雨入渗模型应用到隧道工程中，发现随着降雨时间的增长，隧道拱顶沉降和上部水平位移会不断增大，沉降最大值达306 mm，最大水平位移超过20 mm，且随着降雨时间增长，隧道近裂隙侧的拱顶沉降值和水平位移值变化明显高于远裂隙侧。

（3）分析了影响含裂隙土层隧道降雨入渗变形特性的因素，包括降雨强度、裂隙位置和裂隙宽度。结果表明，降雨持续72 h情况下，靠近裂隙处隧道拱顶沉降和上部水平位移较大；裂隙宽度越大，降雨入渗的渗流扩展越快，隧道产生的位移越大，其中裂隙处于隧道两侧时，拱顶沉降最大值较基本工况增大了3.9%，上部水平位移最大值增大了12%。裂隙处于隧道上方时，拱顶沉降最大值较基本工况增大了6.5%，上部水平位移最大值增大了44%；降雨强度则对隧道围岩的影响较小，3种降雨强度下模型表面入渗率在降雨60 h后均趋于稳定。