王述红, 杨天娇, 刘欢, 佟可蕙
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2017-09-11
基金项目：国家自然科学基金资助项目(U1602232，51474050)；地质灾害防治与地质环境保护国家重点实验室资助项目(SKLGP2014K011)；辽宁省高等学校优秀人才支持计划项目(LN2014006)。
作者简介：王述红(1969-)，男，江苏泰州人，东北大学教授，博士生导师。

摘要：基于非饱和土力学的基本理论, 建立了考虑孔隙介质可压缩性的二维非饱和土坡固-液-气三相耦合控制方程, 利用COMSOL multiphysics软件的PDE平台, 将所建立的耦合控制方程的数值模拟结果与经典Liakopoulos砂柱排水的试验结果进行比较, 验证了所建立的三相耦合控制方程的正确性, 并讨论不同因素对渗流场和位移场的影响.模拟结果表明:渗透率k_s对降雨渗流过程的推进起阻滞作用, 同时延缓了非饱和土坡的变形, 但影响较弱.进气值系数对非饱和土坡变形沉降的过程有一定影响, 但总体影响也较弱.
关键词：非饱和土固-液-气三相耦合渗透率进气值系数降雨
Verification and Analysis of Three-Phase Coupling Model for Soil Slope Under Rainfall Conditions
WANG Shu-hong, YANG Tian-jiao, LIU Huan, TONG Ke-hui
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: WANG Shu-hong, E-mail: shwangneu@126.com
Abstract: Based on the basic theory of unsaturated soil mechanics, a solid-liquid-gas three-phase coupling control equation for two-dimensional unsaturated soil slopes was established by considering the compressibility of porous media. The numerical simulation results of the established coupling control equations based on the PDE platform of COMSOL multiphysics software was compared with the experimental results of the classical Liakopoulos sand column. The established three-phase coupling control equation was verified, and the influences of different factors on the seepage field and displacement field were discussed. The simulation results showed that the permeability can slightly retard the rainfall seepage process and delay the deformation of unsaturated soil slope. The air intake coefficient has a slight impact on the settling process of unsaturated soil slope.
Key words: unsaturated soilsolid-liquid-gas three-phase couplingpermeabilityair intake coefficientrainfall
滑坡是灾难性大、突发性强的地质灾害, 而降雨是诱发滑坡地质灾害最主要的原因, 且降雨积累量越大发生边坡失稳的概率越高.习念念等^[1]从流固耦合模型、水气耦合模型对非饱和土边坡的降雨入渗过程进行分析, 但由于非饱和土复杂的水利特性, 导致其固-液-气三相耦合密不可分.Lam等^[2]将非饱和土的固结理论和水体运移理论相结合, 广泛应用于岩土工程实践的饱和-非饱和渗流控制方程中.熊勇林等^[3]采用Fortran语言自主编制开发了一套水-土-气三相渗流-变形的耦合程序来进行有限元分析.本文基于非饱和土力学理论, 结合降雨入渗过程, 建立了考虑孔隙介质可压缩性的二维非饱和土坡固-液-气三相耦合控制方程, 利用COMSOL multiphysics软件模拟对比分析, 验证了本文所建立的三相耦合控制方程的正确性, 并得到不同因素对渗流场和位移场的影响规律, 因此，在岩土力学多耦合领域，研究所有对固-液-气动态耦合体系的影响具有重要的理论意义和工程实践价值.
1 非饱和土渗流基本理论1.1 土-水特征曲线Van Genuchten^[4]提出的VG土/水特征曲线模型是目前比较经典的模型之一, 体积含水率θ的表达式为

(1)

式中:θ为土体中水的体积分数; θ_r为残余土体中水的体积分数; θ_s为土体中水的饱和体积分数; α₁, m, n均为经验拟合参数; h为基质压力水头.
1.2 渗透性函数模型Van Genuchten土-水特征曲线公式:

(2)

式中, S_e为有效饱和度.采用Mualem^[5]模型, 令m=1-1/n, 得

(3)

其中, k_r为相对渗透系数.
2 固-液-气三相耦合模型2.1 两相流耦合模型流体非饱和带多孔系统的渗流满足质量守恒通式^[6]:

(4)

式中:α是液相(w)或气相(a); ρ_α是流体密度(kg/m³); S_e，α是两种流体的饱和度; q_α是流体的单位流量(m³/s); t是时间变量(s); Q^α是汇源项(kg/m³·s).
由于存在压力梯度, 多孔介质均符合达西定律的线性渗流运动, 其表达式为

(5)

其中:v_α是流体达西流速(m/s); k_int是绝对渗透率; k_rα是相对渗透系数; p_α是流体的压强(Pa); g表示重力加速度(m/s²); η_α是流体的动力黏度(Pa·s).
定义多孔介质中单位面积的基质吸力p_c为

(6)

本文假定二维土坡的孔隙完全被两种流体填充:

(7)

将饱和度随基质吸力变化用比水容量C表示:

(8)

当忽略孔隙介质的可压缩特性时, 将式(4)~式(8)五个方程耦合, 即可得将流体压缩系数引入两相流耦合的控制方程:

(9)

其中，β_α为介质的压缩系数.
式(9)即是水气两相流耦合模型渗流场的控制方程.
2.2 力学平衡方程在二维平面体系中的非饱和单元土体力学平衡方程表达式为

(10)

其中:σ为土体总应力; 为多孔介质三相的总密度, =nS_{e, w}ρ_w+nS_{e, a}ρ_a+(1-n)ρ_s，n为土体孔隙率.
因此得到了非饱和土固-液-气三相的渗流-变形耦合控制方程式(9)、式(10).
3 固-液-气三相耦合模型验证本文采用Liakopoulos砂柱排水试验来验证上述推导模型的正确性^[7-9].
3.1 计算模型此试验采用Del Monte砂柱, 高1 m, 直径0.1 m, 底部设有透水石, 侧部采用圆筒隔水.假定此问题是轴对称问题, 给予一定的边界及初始条件, 建立二维的砂柱排水试验计算模型.计算参数如表 1所示.
表 1(Table 1)

表 1 砂柱模型的物理参数Table 1 Physical parameters of sand column model ks/kPa ρs/(kg·m-3) n0 k0/m2 1E-9 2 850 0.297 5 4.5E-13 kw/kPa ρw0/(kg·m-3) ηw/(Pa·s) ρa0/(kg·m-3) 2E+6 1 000 1.1E-2 1.25 ηa/(Pa·s) p/kPa 1.8E-5 101

表 1 砂柱模型的物理参数 Table 1 Physical parameters of sand column model

3.2 初始条件和边界条件图 1是砂柱初始及边界条件示意图.图中Ⅰ-Ⅰ是对称轴中心截面.
图 1(Fig. 1)

图 1 砂柱初始及边界条件Fig.1 Initial and boundary conditions of sand column

用COMSOL multiphysics建立网格划分图, 此模型包括1个域, 4个边界及4个端点单元.共划分为608个域单元及110个边单元.
3.3 结果分析1) 土体非饱和渗流过程.图 2为不同时刻的孔隙水压力随高度变化的分布规律, 由图可知:在实验初期耦合模型所得的孔隙水压力比砂柱试验数据稍小, 孔隙水压力的变化速率比试验数据大.这是由于在试验初期砂柱的渗流过程及由此渗流过程所引起的变形还未完全稳定, 砂柱中水的排出受到重力作用的同时, 还受到由初始沉降而引起的砂土“挤压”变形的影响.20 min后沉降逐渐完成, 试验数据与模拟结果越来越接近.
图 2(Fig. 2)

图 2 Ⅰ-Ⅰ截面孔隙水压力随高度的变化Fig.2 Change of pore water pressure of Ⅰ-Ⅰ cross section with height

2) 气体迁移过程.图 3为砂柱中心截面Ⅰ-Ⅰ不同时刻孔隙气压随高度变化的分布规律, 由图可知:相对的低气压区首先位于砂柱的顶部位置, 随着时间的推移气体由砂柱顶部向内部逐渐入渗, 低压区逐渐向中部移动, 使上部气压值逐渐增大, 最后整体砂柱内的孔隙气压渐趋于标准大气压.
图 3(Fig. 3)

图 3 Ⅰ-Ⅰ截面孔隙气压随高度的变化Fig.3 Change of pore air pressure of Ⅰ-Ⅰ cross section with height

3) 变形过程.图 4为砂柱中心截面Ⅰ-Ⅰ不同时刻竖向位移随高度变化的分布规律.由图可知:砂柱随孔隙水的排出发生沉降变形, 竖向位移的峰值在砂柱顶部出现, 并且逐渐增大, 前20 min砂柱竖向位移变化的速率较快, 随后速度逐渐减慢, 120 min时其沉降值渐趋稳定.
图 4(Fig. 4)

图 4 Ⅰ-Ⅰ截面竖向位移随高度的变化Fig.4 Change of vertical displacement of Ⅰ-Ⅰ cross section with height

以上三相渗流-变形耦合模型数值分析结果与砂柱排水试验结果基本吻合, 证明本文所建立的三相渗流-变形耦合控制方程的正确性及其应用价值.
4 降雨下非饱和土坡因素分析本文利用COMSOL multiphysics软件分析渗透率k_s及进气值系数α对三相耦合状态下二维边坡计算结果的影响.
4.1 数值模型及参数条件图 5是二维边坡几何模型, 尺寸如图所示.上部是降雨边界, 底边是地下水位边界, 并在x=50 m处设置计算截面Ⅰ-Ⅰ, 表 2是计算物理参数.图 6是二维网格划分图, 共划分为705个域单元, 77个边和6个顶点单元.
图 5(Fig. 5)

图 5 二维边坡模型图Fig.5 2D slope model

表 2(Table 2)

表 2 二维土坡模型的计算物理参数Table 2 Parameters of two-dimensional soil slope model E/MPa μ n ρs/(kg·m-3) 7.5 0.3 0.25 2 850 ρw/(kg·m-3) ρa/(kg·m-3) ηw/(Pa·s) ηa/(Pa·s) 1 000 1.25 1.005E-3 1.79E-5 βw/Pa-1 βa/Pa-1 k0/m2 p/kPa 5E-10 5E-4 4.5E-13 1E+5 θs θr α1 m n 0.185 0.346 0.38 0.275 1.379

表 2 二维土坡模型的计算物理参数 Table 2 Parameters of two-dimensional soil slope model

图 6(Fig. 6)

图 6 二维边坡网格划分Fig.6 Two-dimensional slope grids

初始条件:p_w=-ρ_wgy, p_a=10⁵ Pa.
边界条件:边坡底部u_x=u_y=0, p_w=0, q_a=0;侧边u_x=0, q_w=0, q_a=0;上部q_w=q, 整个入渗边界为透气边界, 设置p_a=10⁵ Pa.
4.2 影响因素分析4.2.1 渗透率k_s的影响渗透率与孔隙率及体积应变的关系式为

(11)

式中:n_t为孔隙率; ε_v为体积应变; k_s0为初始饱和土体渗透率.
图 7为在Ⅰ-Ⅰ截面处考虑k_s耦合变化及k_s不变的条件下孔隙水压力随高度变化的情况.t=24 h时, 这两种条件下坡顶处的孔隙水压力值较接近; t=48 h时, 坡顶处的孔隙水压力差值增大.表明随着降雨入渗的进行, 当两场耦合时, 土体变形且土粒孔隙逐渐被压缩, k_s也随之逐渐减小, 从而对降雨渗流过程的推进起到阻滞作用.因此同一时刻k_s耦合变化的条件下孔隙水压力的增长速率比k_s不变的条件下低.
图 7(Fig. 7)

图 7 不同时刻Ⅰ-Ⅰ截面的孔隙水压力分布图Fig.7 Pore water pressure of Ⅰ-Ⅰ cross section at different 10:09:21

图 8考虑k_s耦合变化及k_s不变条件下的坡面顶点竖向位移随时间变化的分布情况.两种条件下坡面顶点竖向位移变化的趋势基本相同, 降雨初期两种条件下的竖向位移曲线基本重合, 随着降雨入渗的进行竖向位移变化的差值随之增大, 最后最大差值趋于稳定.因此考虑k_s耦合变化对降雨入渗条件下的非饱和土坡变形存在一定影响, 但总体看来此影响较弱.
图 8(Fig. 8)

图 8 坡面顶点的竖向位移分布图Fig.8 Vertical displacement of peak point of slope

4.2.2 进气值系数λ的影响非饱和土体的进气值(AEV)定义为空气刚进入土体时其对应基质吸力的值, λ≈1/AEV, 称为进气值系数.
图 9表示在Ⅰ-Ⅰ截面处不同λ的条件下孔隙水压力随高度变化的分布情况.λ为0.1时孔隙水压力增长速率最慢, 而λ为0.02时其速率最快.表明λ对非饱和土坡的降雨入渗有一定的影响, 同一时刻进气值系数λ越小, 基质吸力消散速度越快.
图 9(Fig. 9)

图 9 Ⅰ-Ⅰ截面的孔隙水压力分布图Fig.9 Change of pore water pressure of Ⅰ-Ⅰ cross section with height

图 10表示在不同λ的条件下坡面顶点竖向位移随时间变化的分布情况.在10~30 h时, λ为0.02时坡面顶点沉降最大, 速率最快; λ为0.1时坡面顶点沉降最小, 速率也最慢.而在10 h前3种条件下的曲线几乎重合, 在30 h后最终沉降值基本稳定在同一值.因此进气值系数λ的变化对降雨入渗条件下的非饱和土坡变形有一定的影响, 但总体看来此影响较弱.
图 10(Fig. 10)

图 10 坡面顶点的竖向位移分布图Fig.10 Vertical displacement of the peak point of the slope

5 结论1) 基于非饱和土力学的基本理论建立了考虑孔隙介质可压缩性的二维非饱和土坡三相耦合控制方程, 经典Liakopoulos砂柱排水试验结果证明该控制方程的正确性.
2) 随着砂柱底部不断排水, 其上部先形成非饱和区且孔隙水压力随之降低, 最后达到稳定状态.同时气体从砂柱顶部不断向内部入渗, 孔隙气压力也逐步增大, 最终逐渐趋向于标准大气压.砂柱的竖向沉降也随着降雨时长的增加逐渐稳定.
3) 渗透率k_s对降雨渗流过程的推进起阻滞作用, 同时也延缓了非饱和土坡的变形, 但影响较弱.进气值系数λ对非饱和土坡变形沉降的过程有一定影响, 但总体影响较弱.
参考文献

[1]	习念念, 童富果, 刘刚, 等. 基于水气二相流的边坡降雨入渗有限元分析[J].长江科学院院报, 2017, 34(1): 120–123, 141. ( Xi Nian-nian, Tong Fu-guo, Liu Gang, et al. Finite element analysis of slope infiltration based on water-air two-phase flow[J].Journal of Yangtze River Scientific Research Institute, 2017, 34(1): 120–123, 141.)
[2]	Lam L, Fredlund D G. Transient seepage model of saturated unsaturated soil systems[J].Cabadian Geotechnical Journal, 1987, 24(4): 565–580.DOI:10.1139/t87-071
[3]	熊勇林, 朱合华, 叶冠林, 等. 降雨入渗引起非饱和土边坡破坏的土-水-气三相渗合有限元分析[J].岩土力学, 2017, 38(1): 284–290. ( Xiong Yong-lin, Zhu He-hua, Ye Guan-lin, et al. Finite element analysis of soil-water-gas three-phase seepage induced by rainfall infiltration on unsaturated soil slope[J].Rock and Soil Mechanics, 2017, 38(1): 284–290.)
[4]	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
[5]	Mualem Y, Klute A. Hydraulic conductivity of unsaturated soils:prediction and formulas[M]. New York: Springer, 1986.
[6]	Richards L A. Capillary conduction of liquids through porous medium[J].Physics, 1931, 1(5): 318–333.DOI:10.1063/1.1745010
[7]	Schrefler R, Lewis R. The finite element method in the static and dynamic deformation and consolidation of porous media[M]. New York: Wiley, 1998.
[8]	Nagel F, Meschke G. An elastic-plastic three phase model for partially saturated soil for the finite element simulation of compressed air support in tunneling[J].International Journal for Numerical and Analytical Methods in Geomechanics, 2010, 34(6): 605–625.
[9]	Fredlund D, Xing A. Equation for the soil water characteristic curve[J].Canadian Geotechnical Journal, 1994, 31(4): 521–532.DOI:10.1139/t94-061