基于多场耦合的膨胀土边坡非饱和降雨入渗分析

韦秉旭，陈亮胜，肖罗明，郑 威，白玉祥

(1.长沙理工大学 交通运输工程学院，长沙 410114； 2.湖南路桥建设集团有限责任公司，长沙 410004)

1 研究背景

降雨诱发膨胀土边坡产生滑坡、崩塌等地质灾害，严重威胁公路工程建设和运营安全，造成了大量的人员伤亡、财产损失和环境破坏[1]。降雨入渗引起土坡湿度场的变化是一个典型的非饱和渗流过程，随着降雨的持续，大气影响深度内坡体的孔隙水压力、含水率和重度显著提高，导致基质吸力消散和抗剪强度降低，逐渐形成暂态饱和区[2]。此外，降雨入渗会引起膨胀土膨胀变形和强度软化等不良性质的出现，加剧边坡地质灾害的形成[3]。因此，掌握降雨入渗对膨胀土边坡非饱和渗流影响的时间和空间分布规律，能为膨胀土边坡防治提供理论基础与技术支持。

数值模拟因其可复现性、求解速度快、数据易采集、节约成本等优点，成为广大学者研究膨胀土边坡非饱和渗流的重要手段。文献[4-6]采用渗流与应力的非耦合数值方法，直接从孔隙水压力、含水率、渗透力和暂态饱和区的时空分布规律来阐述膨胀土边坡内部非饱和渗流的过程，并说明了暂态饱和区与渗流发展速度、面积、时间等诸多因素有关，但目前尚未建立一个描述不同降雨历时期间暂态饱和区变化规律的定量指标，且此方法忽略了土体变形对土骨架和孔隙水的相互作用，导致其渗透特性、水分运移通道和渗流过程与实际存在差异。有限差分软件FLAC3D具有强大的非线性应力应变求解功能，能有效处理非饱和土渗流场与应力场相互耦合的难题。蒋中明等[7]基于内嵌FISH语言的二次开发平台，修正并完善了非饱和渗流的降雨入渗计算功能。汪华斌等[8]、谢强等[9]着重处理入渗边界计算方法与非饱和区流量的两大关键问题，充分验证了有限差分软件FLAC3D实现非饱和土流-固耦合数值模拟的可靠性。实际上，膨胀土的非饱和降雨入渗过程伴随着膨胀应变和应力应变的产生，是典型的膨胀-渗流-应力耦合问题，如何将膨胀土的膨胀应变引入非饱和土的流-固耦合模型成为膨胀土渗流耦合分析的另一难点。

鉴于此，基于有限差分渗流理论、降雨入渗模块的二次开发和膨胀应变的施加，采用C++和FISH语言编制程序，提出一种综合考虑应力应变和膨胀应变的非饱和膨胀土多场耦合分析法。以某膨胀土路堑边坡为研究对象，运用FLAC3D实现了降雨入渗条件下膨胀土边坡非饱和渗流的数值模拟，并探讨了不同降雨时间下膨胀应变对饱和度和暂态面积比变化规律的影响。

2 非饱和膨胀土多场耦合分析法

2.1 有限差分渗流理论

准静态比奥理论是有限差分FLAC3D计算非饱和渗流的方法，它将多孔介质结构的饱和度、孔隙压力和体积应变建立联系[10]，假设非饱和渗流中不考虑气相的影响，其流体连续性方程的表达式为

(1)

式中：s为饱和度；uw为孔隙水压力；εv为体积应变；M为比奥模量；ζ为孔隙水运移过程中发生的流体体积改变量；n为孔隙率；t为时间；α为比奥系数。

当膨胀土边坡经降雨入渗处于饱和-非饱和渗流阶段时，因含水率的变化胀缩现象明显，总应变更新为应力应变和膨胀应变之和[11]，即

(2)

将式(2)代入连续性方程式(1)中，得到考虑膨胀应变的膨胀土连续性方程为

在降雨入渗条件下，非饱和土多孔介质水分运移的过程遵循达西定律和质量守恒定律，其达西流动方程的表达式为

qi=-Kilk(s)(uw-ρfxjgj)il。

(4)

式中：qi为单位流量向量;ρf为流体密度;Kil为渗透系数张量；Kil=ks/γw，ks为达西定律中单元的饱和渗透系数，γw为流体重度;k(s)为相对渗透系数，满足k(s)=s2(3-2s)∈[0，1]；xj为笛卡尔坐标分量;gj为重力加速度分量。

流体质量平衡方程的表达式为

∂ζ/∂t=-qi,i+qv。

(5)

式中qv为流体源强度。

非饱和区的饱和度和渗透系数是随基质吸力动态变化的参数，其值分别通过土-水特征曲线和渗透系数曲线与孔隙水压力建立联系，从而赋值给每一节点或单元。可通过Fredlund-Xing模型[12]确定基质吸力与饱和度的曲线公式，即

(6)

式中：ψ为基质吸力；ψr为残余基质吸力；a为进气值的倒数；m为控制残余含水率的参数，m=3.67ln(1/si)，si为进气值所对应饱和度；n为曲线斜率。

非饱和区单元的渗透系数k为单元饱和渗透系数与相对渗透系数的乘积，其表达式为

k=kss2(3-2s) 。

(7)

2.2 降雨入渗模块二次开发

基于FLAC3D平台的降雨入渗模块二次开发主要涉及初始孔压场的生成、非饱和区饱和度和渗透系数的设置、边界条件的开发三大技术问题。

(1)在FLAC3D中设置流体抗拉强度，能允许负孔隙水压力的存在，这为实现非饱和渗流数值模拟提供了可能。通过内置变量zone.pp(z)可指定节点孔隙水压力，软件内部运用自动插值的方法获取土体的孔隙水压力，即可生成与实际工程适应的初始孔压场。

(2)通过式(6)利用内置变量gp.sat(z)获取非饱和区的节点饱和度，并采用距离加权的方法转化为单元饱和度。土体的渗透系数为单元属性参数，根据相对渗透系数公式，采用内置变量z.fluid.prop (z)计算土体单元的相对渗透系数，并通过式(7)获取非饱和区单元修正后的渗透系数。

(3)FLAC3D计算渗流边界时，通常按降雨强度或饱和渗透系数的大小作为流量边界施加于坡面，事实上，降雨入渗及出渗边界是一个动态变化的边界[13]。采用接触面单元提供的intface功能获取坡面不同位置的渗透系数，调用降雨强度与坡面渗透系数两者之间的小值作为坡面入渗率，并实时监控坡面节点孔隙水压力的大小，如当节点孔隙水压力>0时，节点附近出现积水，将该节点的流量边界调整为压力边界。

根据上述三大技术问题的解决措施编写FISH程序，采用FISHCALL命令调用程序，实现降雨入渗对土质边坡非饱和渗流的模拟。

2.3 膨胀应变的施加

通过开发的降雨入渗模块与修正有效应力后的应力场模块相互耦合，可构成普通土的非饱和流-固耦合模型，但膨胀土非饱和渗流的变化与普通土的区别主要在于膨胀应变的产生[14]。从式(3)膨胀土的流体连续性方程可知，膨胀土是典型的膨胀变形场-非饱和渗流场-应力场多场耦合问题。目前，在FLAC3D内嵌本构模型中，暂无本构模型适用于膨胀土胀缩性这一特殊性质。因此，需修正内嵌本构模型来体现膨胀应变对膨胀土边坡非饱和渗流的影响。

膨胀应变在各向同性弹性体中沿土体的各个方向相同，因此，膨胀应变的主应变方向相同，切应变为0，其膨胀应变分量为

(8)

根据式(2)可知：应力应变和膨胀应变遵循Hooke定律，其应变分量为

(9)

同时将式(9)变换为总应力分量形式，即

(10)

式中：σij和εij分别为总应力分量和总应变分量；σkk和εkk分别为三个方向的主应力之和和主应变之和；E为弹性模量；v为泊松比；δij为Kornecker符号，当i=j时，δij=1，否则取0。

将式(10)改写成有限差分的总应变增量格式，即

(12)

膨胀应变的施加主要在Visual Studio 2010平台中操作，其开发工作基于内嵌本构模型的C++代码，利用上述理论修正头文件(后缀为.h)和源文件(后缀为.cpp)，并根据文献[15]修正屈服准则和流动法则，由此生成一个动态库文件(后缀为.dll)。其中，头文件的修正包括私有变量和成员函数的重新定义，源文件的修正通过实际函数代替基类的成员函数来实现，在成员函数中，主要对Initialize()和Run()函数进行修正，其作用分别为运行初始化模型的中间变量与根据应变增量计算应力增量。

在FLAC3D计算非饱和渗流的流-固耦合分析中，调用修正后的膨胀土本构模型，即可实现非饱和膨胀土的膨胀-渗流-应力多场耦合分析。其原理为：基于有效应力原理建立应力场与非饱和渗流场的联系，可以通过降雨入渗引起非饱和渗流的时空分布控制膨胀应变来影响应力应变关系的变化过程，也可以通过体现应力应变关系的膨胀应变来影响非饱和降雨入渗的变化过程，三者实现相互耦合，简化计算流程如图1所示。

图1 计算流程Fig.1 Computation process

3 工程算例

根据相关工程资料，考虑南阳某膨胀土路堑边坡为计算模型，经长期气候作用，距坡面深度3 m范围内分布着棕黄色膨胀土(上层土)，以下为褐黄色膨胀土(下层土)。图2所示为路堑边坡的网格模型，一级坡和二级坡的边坡坡率为别为1∶1.5和 1∶1.75，中间设台阶连接，宽度为2 m，设计Ⅰ-Ⅰ、Ⅱ-Ⅱ两个截面作为计算剖面。

图2 网格模型Fig.2 Geometric model

表1 数值模拟参数Table 1 Numerical simulation parameters

计算参数：膨胀土边坡的土-水特征曲线和饱和渗透系数分别采用压力板和双环渗透试验，并分别采用式(6)和式(7)对数据进行拟合，其余计算参数根据室内外试验确定，参数具体取值见表1。实际工程中假设膨胀土的干密度一定，膨胀应变主要与初始含水率和含水率的变化情况有关。因此，膨胀应变随含水率的变化关系采用无荷膨胀试验确定[16]，其曲线关系如图3所示，并通过MatLab确定其多元回归的拟合公式为

式中：εw为膨胀应变；w0和w分别为初始含水率和过程含水率。

图3 膨胀应变随含水率的变化Fig.3 Variation of expansion strain with water content

边界及水文条件：模型底部固定约束，顶部自由，四周水平约束。在降雨入渗数值计算过程中，路面设为不透水边界，坡顶水平面和坡面设置为降雨入渗及出渗边界。值得说明的是，当降雨强度小于坡面入渗能力时，形成无压入渗阶段，此阶段为流量边界；当降雨强度大于坡面入渗能力时，形成有压入渗阶段，此阶段为压力边界。在降雨结束后，根据坡面入渗及出渗的情况重新调整，当水分由于重力作用向坡内入渗时，坡面孔隙水压力转变为负值，此时调整为零流量边界；当水分向坡外溢出时，坡面孔隙水压力转变为正值，此时调整为零压力边界。

根据当地气象资料拟定降雨强度为1.65×10-6m/s，降雨持续3 d，停雨持续5 d，其初始孔压场如图4所示，并采用暂态分析法每隔12 h记录一次计算结果，实时监测边坡模型的变化情况。

图4 初始孔压场Fig.4 Initial pore water pressure

4 计算结果及分析

4.1 孔隙水压力分布

利用本文提出的多场耦合分析法模拟非饱和膨胀土的降雨入渗过程，得到如图5所示不同降雨历时阶段的孔隙水压力变化。从图5可以看出：降雨过程中，孔隙水压力的变化与降雨时间正相关，坡面最先出现暂态水压力，随降雨时间推移，坡内的孔隙水压力也逐渐上升；降雨结束后，水分因重力作用会继续向坡内运移，湿润锋逐渐被均化，其上缘的孔隙水压力减小，下缘的孔隙水压力增加，且随深度增加，土体固结度提高，湿润锋向坡内扩展速率变慢，介于-90～0 kPa之间孔隙水压力的区域增加。还可以看出：由于下层土的渗透能力小于上层土，湿润锋难以透过渗透率低的下层土，导致降雨入渗对非饱和渗流场的影响区域集中在上层土，对下层土孔隙水压力的影响程度低，在土层交界处形成一条显著的自上而下孔隙水压力逐级递减带，具有距坡面位置越远，受降雨入渗影响程度越小的分布特征。其数值模拟结果满足文献[17]和文献[18]的非饱和膨胀土降雨入渗分析，说明本文提出的多场耦合分析法可进一步研究膨胀应变对膨胀土非饱和降雨入渗的变化规律。

图5 孔隙水压力分布Fig.5 Distribution of pore water pressure

4.2 膨胀应变对饱和度的影响规律

根据边坡初始饱和度的状况，其区域可划分为暂态饱和区、非饱和区、初始饱和区。实际降雨入渗过程中，土体孔隙难以全部被水分填充，故将饱和度>0.9的区域划分为饱和区[19]。

图6为边坡截面Ⅰ-Ⅰ和截面Ⅱ-Ⅱ处的饱和度随边坡深度变化的曲线关系。由图6可知:①降雨3 d时，暂态饱和区的范围成悬挂型向坡内扩展，降雨入渗深度和饱和区深度分别处于2.5～3.0 m、1.8～2.3 m之间，且不考虑膨胀变形的降雨入渗深度和暂态饱和区的扩展深度要大于考虑膨胀应变的情况；②停雨5 d时，湿润锋上缘变为非饱和区，降雨入渗深度介于3.5～4.0 m之间，土层交界处存在较大范围的饱和区，暂态饱和区的面积具有显著的时间滞后性；③截面Ⅰ-Ⅰ和Ⅱ-Ⅱ饱和度变化的差异性主要表现在土层交界处，且考虑膨胀应变的饱和度变化幅度小于不考虑膨胀应变的情况。

图6 饱和度变化Fig.6 Variation of saturation

造成上述现象的原因是：降雨入渗改变了土体含水率的时空分布，引起含水率发生变化的中下部土体膨胀应变变化量显著增加，产生湿润锋的影响范围和水分运移的扩散速度受控于膨胀应变，当考虑膨胀应变时，膨胀土在降雨增湿过程中产生膨胀变形，因吸水膨胀土体被压密，引起土颗粒间的孔隙尺寸减小，水分经渗流通道的连通孔隙变窄，从而造成膨胀土的渗透性能降低，水分扩散速率变缓，抑制了降雨入渗对水分扩散的过程。在降雨结束后，坡面附近的饱和度降低，膨胀土开始出现失水收缩的现象，引起膨胀应变衰减，但其饱和度在降雨入渗过程中始终大于初始饱和度，膨胀应变继续发挥着抑制降雨入渗对水分扩散的作用。

4.3 膨胀应变对暂态面积比的影响规律

为定量分析降雨入渗过程中膨胀土边坡暂态饱和区面积的时程变化，定义暂态面积比η(暂态饱和区面积与初始非饱和区面积之比)来揭示膨胀土边坡非饱和降雨入渗的变化规律，用以弥补描述暂态饱和区时程效应的不足。暂态面积比的表达式为

式中：ST(m2)和S(m2)分别为暂态饱和区面积和初始非饱和区面积；ST(pt)和S(pt)分别为暂态饱和区像素面积和初始非饱和区像素面积。运用图像处理技术，根据两者像素之比可直接获取暂态面积比。

图7所示为考虑膨胀应变与不考虑膨胀应变的暂态面积比时程变化，可以看出，随降雨历时的增加，暂态面积比先上升后下降，呈“凸”型分布，且暂态面积比的“凸”点发生在停雨前期，其极大值为21.5%。这是由于膨胀土边坡在降雨期间和停雨前期，湿润锋呈饱和状态向坡内扩展，悬挂型暂态饱和区面积持续上升；随停雨时间的持续，湿润锋抵达土层交界处，暂态饱和区的面积向下扩展受阻，而坡面水分不断向坡外溢出或向坡内慢渗，湿润锋被均化，暂态面积比呈不断缩小的趋势，且暂态面积比的降低速率明显小于增长速率。在整个降雨历时期间，针对考虑膨胀应变和不考虑膨胀应变的情况，暂态面积比的“凸”点形成时间分别在停雨2 d和1 d；降雨结束5 d后，暂态面积比分别变为13.1%和10.2%，相差2.9%。其原因是膨胀应变对非饱和降雨入渗过程起抑制作用，水分扩散作用变缓，考虑膨胀应变的暂态面积比的变化速率变慢，其暂态面积比的雨后时间滞后效应长于不考虑膨胀应变时的情况。

图7 暂态面积比的时程变化Fig.7 Change in transient area ratio over time

膨胀土边坡非饱和降雨入渗的数值结果验证了本文提出的多场耦合分析法的正确性和合理性。在非饱和降雨入渗过程中，一方面，因降雨增湿产生的膨胀变形降低了土体的渗透率，抑制了降雨入渗对水分的扩散作用；另一方面，膨胀土的膨胀应变受非饱和降雨入渗过程的控制，导致膨胀应变的发展过程不同步，造成膨胀应变变化量的大小差异幅度也不尽相同。因此，采用数值模拟的方法应综合考虑膨胀土的胀缩性对降雨入渗过程的影响，从而进一步分析多场耦合作用下的非饱和膨胀土边坡的破坏机制与变形机理。

5 结 论

(1)根据有限差分渗流理论和增量理论，将膨胀应变引入非饱和土的流-固耦合模型，采用FISH语言和C++语言编制程序，进行了非饱和降雨入渗模块和膨胀土本构模型的二次开发，提出了一种适用于非饱和膨胀土的多场耦合分析法。

(2)湿润锋影响范围集中在边坡浅层，湿润锋随深度增加向坡内扩展速率变慢，雨后介于-90～0 kPa间的孔隙水压力区域增加，且土层交界处形成一条显著的自上而下孔隙水压力逐级递减带。

(3)降雨增湿产生的膨胀变形造成土体渗透性能降低，抑制了降雨入渗对水分的扩散作用，导致膨胀土边坡的降雨入渗深度、湿润锋扩展和均化速度以及饱和区时空分布受膨胀应变的控制。

(4)暂态面积比的时程效应呈“凸”型分布，其“凸”点发生在停雨阶段，且在考虑膨胀应变后，暂态面积比的变化速率会变慢，且雨后时间滞后效应相应延长。