0 引言

我国滑坡地质灾害频发，工程中针对不同成因的滑坡常采取不同的治理措施。已有研究^[1-2]表明，绝大多数滑坡都是直接或间接地由坡内地下水渗流场变化而引发的，因此截排水措施几乎在所有的滑坡治理中都会用到。现有常用的截排水措施包括：地表排水沟、地下排水洞、排水盲沟、水平排水洞等，这些排水措施已经在滑坡防治中得到了较为成熟的应用，但是这些措施都是利用重力势能进行排水，而新近提出的充气截排水技术^[3-5]通过在坡体后缘水位线以下充入有压气体，排开土体孔隙水形成非饱和土截水帷幕来达到截渗减排的目的。这一过程涉及到复杂的固、液、气三相耦合效应，而已有的关于充气截排水的研究只考虑了水-气两相流相互作用的的过程。如：康剑伟等^[5]在不考虑充气过程中土体变形的条件下，基于水-气二相流理论，模拟了充气阻渗过程；杜丽丽等^[6]通过数值模拟耦合气-水两相，探究了渗透系数、充气气压、孔隙率对截排水效果的影响；余文飞等^[7]建立了二维长边坡数值模型，探究了充气时坡体内不同阶段水-气二相流具体的运动特征。另外已有的物理模型试验^{[3, 8]}也证明了充气截排水的有效性，并研究了相关的影响因素，但以上研究均未考虑土体骨架可能发生的变形。本文首先通过数值模拟求得充气过程中瞬时非饱和渗流场的解答，再将得到的结果应用到土体应力应变计算中，通过非完全耦合计算分析最终得到充气引起土体变形的规律。

1 水-气运移与土体变形机理分析

充气截排水技术的关键在于往坡体后缘冲入有压气体，形成非饱和截水帷幕，从而减少地下水入渗到下游潜在滑坡区。对于非饱和截水帷幕形成过程的研究是十分重要的，这个过程涉及到土体中水-气相互作用并不断发生运移，同时由于孔隙气压力与水压力的作用，在宏观上土体骨架在不同阶段也表现出相应的变形特征。余文飞等^[7]结合数值分析结果得到充气截排水过程大致可分为3个阶段，即充气非饱和区形成与扩展阶段、不稳定两相流阶段和充气截排水工作阶段。起初，充气形成的非饱和区在坡体内部地下水位线以下持续扩大，但气体未逸出坡体；紧接着，非饱和区突破坡体地下水位线，使得原本被封存在非饱和区内的有压气体在短时间内快速释放到坡体外，导致非饱和区内气压迅速下降，水-气二相流的平衡被打破，水-气界面的形状和运移处于非稳定状态，但非饱和区范围仍在不断扩大，下游坡体内地下水位由于非饱和区拦截的作用而逐渐降低，截排水效果随时间具有波动性；最后，随着充气的进行，从充气点到坡面形成相对固定的气体逸出通道，坡体内水的渗流路径也重新趋于稳定，此时非饱和区形态处于相对固定状态。以上过程很好地说明了充气截排水的工作机理，但是并没有考虑充气过程中土体的变形。由于向土体中充气，刚开始气压力会对土体有一个托顶的作用使土体发生隆起变形，而当气体逸出后土体又会发生沉降变形。此过程类似于浅层气的形成和释放^[9-11]。浅层气作为封闭的压力气体形成后可使地层微量膨胀，气体与土颗粒一起承担上覆的荷重，当气体在一定条件下逸出时就会发生以垂直沉降为主的变形。但当气体在无控制条件下释放时，含气层气压力急剧下降，快速的气流对土体产生冲刷，带走大量的土、砂颗粒，并且会严重扰动含气土层，降低土体强度。因此上述相关研究对于同样要在土层中形成含气非饱和区的充气截排水而言具有很强的借鉴意义。

2 控制方程与模型建立

本文先利用岩土软件GEO-STUDIO中的SEEP/W和AIR/W模块耦合，模拟得到坡体充气情况下的瞬时非饱和渗流场；然后将得到的孔压改变导入SIGMA/W模块中的应力-应变计算的每个荷载步，得到相应孔压改变的体积位移变化。上述过程就是非完全耦合分析。相比于完全耦合分析^[12]，非完全耦合分析可以先获取对非饱和渗流域的完全了解，有助于得到正确的水利边界条件、恰当的材料性质和合适的时间序列，更易于控制且数值上优于同时求解渗流-变形的耦合方程。

2.1 充气时水-气两相运移控制方程

根据质量守恒定律，在等温条件下，忽略非饱和土体内的化学、相变等因素的影响，建立同时包含孔隙气压力和孔隙水压力的质量守恒方程。假设水为不可压缩流体，则密度不变，所以其质量守恒方程可以用体积守恒方程(式(1))来表示；而考虑气体可压缩性的质量守恒方程可用式(2)^[13]表示。

式中：m_w为土水特征曲线在某一特定孔隙水压力处的斜率；γ_w为水的重度；t为时间；K_w为随含水率变化而变化的渗透系数；H_w为渗流过程中的总水头；p_a为孔隙气压力；Q_w为水的单宽流量；θ_a为空气的单位体积分数；T为温度常数；ρ_a为气体密度；K_a为随含水率变化而变化的透气系数；γ_oa为标准状态下气体重度；ρ_oa为标准状态下气体密度；R为气体常数，对于干燥空气，R=287 J/(kg·K)。

在非饱和土中，基质吸力是由存在于水气交界面的收缩膜而产生的，其值等于孔隙气压力和孔隙水压力之差，且与水的体积分数有着密切的关系，它们之间的关系通过土水特征曲线来反映。本文通过Van Genuchten模型^[14]来定义土水特征曲线，见式(3)；并在此基础上确定水、气渗透率曲线，见式(4)、(5)。

式中：θ_w为水体积分数；θ_s为饱和水体积分数；θ_r为残余水体积分数；Ψ为基质吸力；a，n，m均为拟合经验参数，n=1/(1－m)；K_s为饱和土的渗透系数；K_g为透气系数；K_dg为干土透气系数；s_w为有效饱和度；q为系数，q=2.9。

式(3)—(5)联立质量守恒方程共同构成了描述坡体中水-气二相运移数学模型的控制方程。

2.2 渗流-变形耦合方程

本文假设所用土体材料的本构模型为线弹性模型，对于二维问题，非饱和土介质的增量应变-应力关系^[15]可以写成

式中：σ_x，σ_y，σ_z分别为x，y，z方向正应力；τ_xy为垂直于x方向沿y方向切应力；ε_x, ε_y, ε_z分别为x，y，z方向正应变；γ_xy为垂直于x方向沿y方向切应变；u_a为孔隙气压力；u_w为孔隙水压力; u_a-u_w为基质吸力；μ为泊松比；E为土的刚度模量，可以是常数或上覆土层有效应力的函数。

2.3 数值计算模型及参数选取

采用二维边坡进行模拟分析，具体模型及网格划分如图 1所示。滑坡体位于基岩上方，滑面位置固定，坡体厚12 m、长100 m、坡度21°。充气点位置坐标为(25，31)。

首先进行水-气二相流的模拟。自然条件下未充气时的边界条件：坡体后缘AB水头为49 m，坡体表面AD为透水透气边界，坡体右侧CD为潜在渗流面边界但不透气，滑面以下基岩面BC设为不透水边界。充气时的边界条件：水相边界条件同上；气相边界条件为，充气位置选取合适的充气气压值，坡体表面AD为透气边界，气压大小等于大气压的值，坡体两侧AB、CD和底部BC为不透气边界。采用粉砂作为模型材料，Van Genuchten模型的拟合参数取值为a=7，n=1.86，m=0.462，非饱和土渗透系数函数和透气系数函数只需给定饱和土的渗透系数K_s和干土的透气系数K_dg即可确定，结合工程地质手册^[16]及工程经验，取K_dg=10K_s=3.0×10^-3 cm/s，θ_s=0.352 m³/m³，θ_r=0.040 m³/m³。之后将水-气二相流模拟结果导入SIGMA/W模块进行计算，假定坡体AB、CD边界只能沿竖向移动，BC边界水平和竖向都固定，选取材料本构模型为线弹性模型，各参数取值如表 2所示。

3 数值模拟及结果分析 3.1 非饱和区水-气运移过程模拟

首先进行稳态分析，得到未充气时自然状态下的边坡初始渗流场，结果如图 2所示。

然后计算充气状态下的非饱和瞬态渗流场的变化。图 3所示是当充气压力p=100 kPa时，充气引起的非饱和区形态以及渗流场的变化过程，从阶段1到阶段4，图中颜色越深的区域代表饱和度越大，箭头越大代表流速越大。充气刚开始时，气体排开土体孔隙中的水，非饱和区域范围较小，且离充气点越近的地方饱和度越小，同时水流流动方向发生改变，但坡内地下水位线无明显变化(图 3a)；随着充气的进行，非饱和区逐渐扩大，且坡体内的地下水位线有明显的上升，这是由于非饱和区大小达到一定程度后严重阻碍了后缘来水的渗流路径，水流只能向坡面方向流去导致地下水位局部抬升(图 3b)；紧接着非饱和区突破地下水位线，与下游水位线以上的非饱和区域连通，气体开始逸出坡体，此时非饱和区气压力迅速降低，水气界面向气相方向运动挤压气体，当气压力增大到一定程度时水气界面又向液相方向运动，非饱和区以此规律不断发展直到达到水气平衡状态，而下游水位由于非饱和区对后缘来水的拦截逐渐下降(图 3c)；最后到达稳定状态(图 3d)。

另外，当充气压力较小时，虽然非饱和区不断扩大，但是始终无法突破地下水位线，图 4所示是充气压力p=80 kPa时充气稳定后的渗流场。可以推测气体驱替该气压下所能驱动的部分土体孔隙通道中的水体，使得该区域饱和度下降，且由于土体的封气性，形成一个有明显水气分界面的稳定非饱和区，从而在一定程度上拦截后缘来水。

当充气压力过大时，虽然非饱和区从阶段1到阶段4的时间大大缩短，但是由于气压力超过此坡体的进气气压上限值^[4]，水气运移的平衡已被打破，数值上不再收敛，物理上表现为渗流场紊乱。图 5为充气压力p=120 kPa时某一时刻水气相互作用形成的不稳定紊乱渗流场。

3.2 坡体非饱和区变形模拟结果

将上述计算得到非饱和渗流场作为已知条件导入到SIGMA/W模块进行变形模拟分析。选取图 1中1(20，44)、2(25，42)、3(30，40)、4(35，38)4个测点，模拟在不同充气压力下各点竖向位移随时间的变化，探究充气压力大小及渗流场变化对坡体变形的影响。在充气压力p=100 kPa的条件下，4个测点的竖向位移随时间的变化如图 6所示。

从图 6可以看出，所有点的竖向位移都随时间先增大后减小，最后趋于稳定。这是因为有压气体排开土体孔隙水形成一个相对封闭的有压非饱和区后会对周围土体有托顶挤压的作用，从而使地表土体发生隆起变形；而当封闭非饱和区与外界联通气体逸出时，地表土体又会发生沉降变形，所以各测点的竖向位移形成了如图 6所示的变化趋势。而充气稳定后其中测点1和2的竖向位移为正，测点3的竖向位移约为0，测点4的竖向位移为负。这是因为测点1和2的位置相对于充气非饱和区偏上，且位于水气交界面的水相一侧，气体压力对上游土体有顶托作用且阻碍其向下游的变形，故其竖向位移为正；测点3和4的位置相对于充气非饱和区偏下，且位于水气交界面的气相一侧，其中测点3受向上托顶的力产生的正向竖向位移与受沿坡向推挤的力产生的负向竖向位移相当，故竖向位移约为0，而测点4在充气稳定后受沿坡向推挤的力产生的负向竖向位移占优势，所以最后竖向位移为负。

3.3 不同充气压力下各测点位移变化分析

为进一步探究充气压力大小对坡体渗流场及土体变形的影响，以及渗流场变化与土体变形之间的联系，选取不同的充气压力进行数值模拟，得到各测点在不同充气压力下竖向位移随时间的变化规律。由图 7a、7b可知，测点1和测点2的在不同充气压力下竖向位移的变化规律类似。当充气压力较小(p≤90 kPa)时，在模拟的时间范围内竖向位移随着充气进行缓慢增大；当充气压力较大(90 kPa < p≤110 kPa)时，竖向位移先增大后减小最后趋于一个相对稳定的值；当充气压力继续增大(p>110 kPa)时，竖向位移开始出现波动的不稳定现象。产生上述不同现象的原因是：充气压力较小时，非饱和区在坡体内不断扩大，气压力缓慢向上托顶上覆土层，但是非饱和区始终不能突破坡体表面；充气压力较大时，非饱和区突破坡体表面前气压力逐渐向上托顶上覆土层，非饱和区突破坡体表面后，气体逸出坡体气压力减小，上覆土层出现沉降；充气压力继续增大超过进气气压上限值时，过大的充气压力会造成渗流场紊乱，同时破坏土体结构，此时气体通过裂隙逸出坡体，非饱和区气压力下降，坡体表面土体发生下陷，而当周围的水重新挤压非饱和区气体，孔隙气压力又开始上升，坡体表面土体发生隆起变形；随后气压增大到一定程度又会逸出导致坡体表面土体发生下陷，如此坡体表面土体竖向位移做类似的上下波动变化，直到最后达到一个相对稳定的状态。

由图 7c可知，在充气压力p=120 kPa时，测点3的竖向位移呈现无规律的上下波动变化。该测点在不同充气压力下初期竖向位移所能达到的峰值随充气压力的增大而增大，而最后稳定的竖向位移值均较小(y < 0.003 m)。这是因为测点3受向上托顶的力产生的正向竖向位移与受沿坡向推挤的力产生的负向竖向位移相当，故位移较小。由图 7d可知，测点4在充气后期竖向位移都为负值，当充气压力超过100 kPa后，测点最后的竖向位移大小无规律，即位于非饱和区域下游的测点的位移量不一定随充气压力的增大而相应增大，而是具有一定的不确定性。

4 结论

1) 在一定且合适的充气压力下，充气非饱和区的形成可分为4个阶段，即充气初始阶段、非饱和区扩大水位抬升阶段、不稳定两相流阶段和非饱和区稳定阶段。而当充气压力较小时，充气形成的非饱和区无法突破坡体表面，仅停留在水位抬升阶段；充气压力过大时，无法形成稳定非饱和区达到截排水工作状态。

2) 在某一特定且合适的充气压力下，坡体表面土体的竖向位移均随时间呈现先增大后减少的规律，与土体中非饱和区的发展具有一定协同性，最后稳定时越靠近坡体上游的土体整体竖向位移越大。

3) 不同充气压力下，偏向充气点上游的坡体表面土体竖向位移随充气压力的增大而增大，偏向充气点下游的坡体表面土体竖向位移的变化规律则不明显。当充气压力过大时，无法形成稳定非饱和区，在土体变形上表现为竖向位移变化无规律。

4) 采用弹性应力-应变关系对充气截排水过程中土体变形的预测基本反映了其变形趋势，具有一定的参考价值，为定量描述非饱和土渗流以及变形的耦合分析问题提供了有效途径。