2018, Vol. 39
滑坡是灾难性大、突发性强的地质灾害, 而降雨是诱发滑坡地质灾害最主要的原因, 且降雨积累量越大发生边坡失稳的概率越高.习念念等[1]从流固耦合模型、水气耦合模型对非饱和土边坡的降雨入渗过程进行分析, 但由于非饱和土复杂的水利特性, 导致其固-液-气三相耦合密不可分.Lam等[2]将非饱和土的固结理论和水体运移理论相结合, 广泛应用于岩土工程实践的饱和-非饱和渗流控制方程中.熊勇林等[3]采用Fortran语言自主编制开发了一套水-土-气三相渗流-变形的耦合程序来进行有限元分析.本文基于非饱和土力学理论, 结合降雨入渗过程, 建立了考虑孔隙介质可压缩性的二维非饱和土坡固-液-气三相耦合控制方程, 利用COMSOL multiphysics软件模拟对比分析, 验证了本文所建立的三相耦合控制方程的正确性, 并得到不同因素对渗流场和位移场的影响规律, 因此，在岩土力学多耦合领域，研究所有对固-液-气动态耦合体系的影响具有重要的理论意义和工程实践价值.
1 非饱和土渗流基本理论 1.1 土-水特征曲线Van Genuchten[4]提出的VG土/水特征曲线模型是目前比较经典的模型之一, 体积含水率θ的表达式为
|
(1) |
式中:θ为土体中水的体积分数; θr为残余土体中水的体积分数; θs为土体中水的饱和体积分数; α1, m, n均为经验拟合参数; h为基质压力水头.
1.2 渗透性函数模型Van Genuchten土-水特征曲线公式:
|
(2) |
式中, Se为有效饱和度.采用Mualem[5]模型, 令m=1-1/n, 得
|
(3) |
其中, kr为相对渗透系数.
2 固-液-气三相耦合模型 2.1 两相流耦合模型流体非饱和带多孔系统的渗流满足质量守恒通式[6]:
|
(4) |
式中:α是液相(w)或气相(a); ρα是流体密度(kg/m3); Se，α是两种流体的饱和度; qα是流体的单位流量(m3/s); t是时间变量(s); Qα是汇源项(kg/m3·s).
由于存在压力梯度, 多孔介质均符合达西定律的线性渗流运动, 其表达式为
|
(5) |
其中:vα是流体达西流速(m/s); kint是绝对渗透率; krα是相对渗透系数; pα是流体的压强(Pa); g表示重力加速度(m/s2); ηα是流体的动力黏度(Pa·s).
定义多孔介质中单位面积的基质吸力pc为
|
(6) |
本文假定二维土坡的孔隙完全被两种流体填充:
|
(7) |
将饱和度随基质吸力变化用比水容量C表示:
|
(8) |
当忽略孔隙介质的可压缩特性时, 将式(4)~式(8)五个方程耦合, 即可得将流体压缩系数引入两相流耦合的控制方程:
|
(9) |
其中，βα为介质的压缩系数.
式(9)即是水气两相流耦合模型渗流场的控制方程.
2.2 力学平衡方程在二维平面体系中的非饱和单元土体力学平衡方程表达式为
|
(10) |
其中:σ为土体总应力; 
因此得到了非饱和土固-液-气三相的渗流-变形耦合控制方程式(9)、式(10).
3 固-液-气三相耦合模型验证本文采用Liakopoulos砂柱排水试验来验证上述推导模型的正确性[7-9].
3.1 计算模型此试验采用Del Monte砂柱, 高1 m, 直径0.1 m, 底部设有透水石, 侧部采用圆筒隔水.假定此问题是轴对称问题, 给予一定的边界及初始条件, 建立二维的砂柱排水试验计算模型.计算参数如表 1所示.
|
表 1 砂柱模型的物理参数 Table 1 Physical parameters of sand column model |
图 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 Change of pore water pressure of Ⅰ-Ⅰ cross section with height |
2) 气体迁移过程.图 3为砂柱中心截面Ⅰ-Ⅰ不同时刻孔隙气压随高度变化的分布规律, 由图可知:相对的低气压区首先位于砂柱的顶部位置, 随着时间的推移气体由砂柱顶部向内部逐渐入渗, 低压区逐渐向中部移动, 使上部气压值逐渐增大, 最后整体砂柱内的孔隙气压渐趋于标准大气压.
|
图 3 Ⅰ-Ⅰ截面孔隙气压随高度的变化 Fig.3 Change of pore air pressure of Ⅰ-Ⅰ cross section with height |
3) 变形过程.图 4为砂柱中心截面Ⅰ-Ⅰ不同时刻竖向位移随高度变化的分布规律.由图可知:砂柱随孔隙水的排出发生沉降变形, 竖向位移的峰值在砂柱顶部出现, 并且逐渐增大, 前20 min砂柱竖向位移变化的速率较快, 随后速度逐渐减慢, 120 min时其沉降值渐趋稳定.
|
图 4 Ⅰ-Ⅰ截面竖向位移随高度的变化 Fig.4 Change of vertical displacement of Ⅰ-Ⅰ cross section with height |
以上三相渗流-变形耦合模型数值分析结果与砂柱排水试验结果基本吻合, 证明本文所建立的三相渗流-变形耦合控制方程的正确性及其应用价值.
4 降雨下非饱和土坡因素分析本文利用COMSOL multiphysics软件分析渗透率ks及进气值系数α对三相耦合状态下二维边坡计算结果的影响.
4.1 数值模型及参数条件图 5是二维边坡几何模型, 尺寸如图所示.上部是降雨边界, 底边是地下水位边界, 并在x=50 m处设置计算截面Ⅰ-Ⅰ, 表 2是计算物理参数.图 6是二维网格划分图, 共划分为705个域单元, 77个边和6个顶点单元.
|
图 5 二维边坡模型图 Fig.5 2D slope model |
|
|
表 2 二维土坡模型的计算物理参数 Table 2 Parameters of two-dimensional soil slope model |
|
图 6 二维边坡网格划分 Fig.6 Two-dimensional slope grids |
初始条件:pw=-ρwgy, pa=105 Pa.
边界条件:边坡底部ux=uy=0, pw=0, qa=0;侧边ux=0, qw=0, qa=0;上部qw=q, 整个入渗边界为透气边界, 设置pa=105 Pa.
4.2 影响因素分析 4.2.1 渗透率ks的影响渗透率与孔隙率及体积应变的关系式为
|
(11) |
式中:nt为孔隙率; εv为体积应变; ks0为初始饱和土体渗透率.
图 7为在Ⅰ-Ⅰ截面处考虑ks耦合变化及ks不变的条件下孔隙水压力随高度变化的情况.t=24 h时, 这两种条件下坡顶处的孔隙水压力值较接近; t=48 h时, 坡顶处的孔隙水压力差值增大.表明随着降雨入渗的进行, 当两场耦合时, 土体变形且土粒孔隙逐渐被压缩, ks也随之逐渐减小, 从而对降雨渗流过程的推进起到阻滞作用.因此同一时刻ks耦合变化的条件下孔隙水压力的增长速率比ks不变的条件下低.
|
图 7 不同时刻Ⅰ-Ⅰ截面的孔隙水压力分布图 Fig.7 Pore water pressure of Ⅰ-Ⅰ cross section at different 10:09:21 |
图 8考虑ks耦合变化及ks不变条件下的坡面顶点竖向位移随时间变化的分布情况.两种条件下坡面顶点竖向位移变化的趋势基本相同, 降雨初期两种条件下的竖向位移曲线基本重合, 随着降雨入渗的进行竖向位移变化的差值随之增大, 最后最大差值趋于稳定.因此考虑ks耦合变化对降雨入渗条件下的非饱和土坡变形存在一定影响, 但总体看来此影响较弱.
|
图 8 坡面顶点的竖向位移分布图 Fig.8 Vertical displacement of peak point of slope |
非饱和土体的进气值(AEV)定义为空气刚进入土体时其对应基质吸力的值, λ≈1/AEV, 称为进气值系数.
图 9表示在Ⅰ-Ⅰ截面处不同λ的条件下孔隙水压力随高度变化的分布情况.λ为0.1时孔隙水压力增长速率最慢, 而λ为0.02时其速率最快.表明λ对非饱和土坡的降雨入渗有一定的影响, 同一时刻进气值系数λ越小, 基质吸力消散速度越快.
|
图 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 Vertical displacement of the peak point of the slope |
1) 基于非饱和土力学的基本理论建立了考虑孔隙介质可压缩性的二维非饱和土坡三相耦合控制方程, 经典Liakopoulos砂柱排水试验结果证明该控制方程的正确性.
2) 随着砂柱底部不断排水, 其上部先形成非饱和区且孔隙水压力随之降低, 最后达到稳定状态.同时气体从砂柱顶部不断向内部入渗, 孔隙气压力也逐步增大, 最终逐渐趋向于标准大气压.砂柱的竖向沉降也随着降雨时长的增加逐渐稳定.
3) 渗透率ks对降雨渗流过程的推进起阻滞作用, 同时也延缓了非饱和土坡的变形, 但影响较弱.进气值系数λ对非饱和土坡变形沉降的过程有一定影响, 但总体影响较弱.
| [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 |