基于水-气二相流模型的土坡稳定性分析

张晓悦,王 栋,张晓乐,沈跃军

(1.浙江水利水电专科学校水利工程系,浙江 杭州 310018;2.信息产业电子第十一设计研究院科技工程股份有限公司杭州分院,浙江杭州 310020;3.浙江华东工程安全技术有限公司,浙江杭州 310014;4.浙江省钱塘江管理局勘测设计院,浙江杭州 310016)

土坡稳定分析是估计天然及人工边坡安全系数的一个常用工具,一般进行边坡稳定分析时,由于测量负孔隙水压力和孔隙气压力并将其纳入稳定分析中比较困难,对于地下水位以上非饱和区中由负孔隙水压力和孔隙气压力所提供的部分抗剪强度通常予以忽略不计。近年来,随着人们对多相流的研究越来越多,在土坡稳定分析中全面考虑孔隙水压力和孔隙气压力成为可能。例如,Dijke等[1]给出了一个较完整的多相流模型,Matthew等[2]提出二相流模型的一种空间和时间离散方法,White等[3]分析了多孔介质饱和度和相对渗透系数关系,Laroche等[4]对多相流模型中毛细压力和饱和度的关系做了研究。

本文建立水-气二相流模型来模拟土体边坡在稳定渗流情况和降雨情况下的水相和气相渗流场,全面考虑土体中的水相和气相流动以及水相和气相成分的相互转换,利用二相流计算结果求得边坡安全系数,分析孔隙气压力和负孔隙水压力在土坡稳定中的作用。

1 水-气二相流模型

1.1 定义以及基本假设

非饱和土是四相混合体,包括固相、气相、液相以及称为收缩膜的水气分界面[5]。在土体中可能有多种液相,本文只考虑水相。非饱和土中的渗流为包括气相和液相流动的二相流,气相和液相都包括水和空气两种组分。液相和气相通过蒸发与溶解相互转换[5],相态及其相互转换关系见图1。

图1 相态及相态转换

模型假设土体骨架不变形,即孔隙率为常数,孔隙流体不可压缩,忽略化学生物学反应,恒温,不考虑本构关系中渗透系数与饱和度之间的滞后作用,气相中所有组分都遵守理想气体定律,气相中水组分分压力等于此温度下的饱和蒸汽压力。

1.2 模型控制方程及其求解

模型考虑气相和液相的流动以及各组分之间的质量传递,模型的控制方程即组分κ的质量守恒方程[6]:

土体固有渗透系数:

式中:下标 α和κ分别代表相态和组分,取g时为气相或空气组分,取w时为水相或水组分;φ为土体孔隙率;Sα为α相的饱和度;ρm,α=nα/Vα,nα为 α相的摩尔质量,Vα为nα对应的体积,m3;ρα为α相的密度;xα,κ为α相中组分κ的摩尔分数;t为时间;ksw为土的饱和渗透系数;krα为α相的相对渗透系数;pα为α相的压力;μα为α相的绝对黏滞性,温度一定的情况下为常数;g为重力加速度;Qκ为κ组分的源汇项,kg/s。

控制方程的变量分为主要变量和次要变量两类,次要变量为主要变量的函数。采用有限差分法对控制方程进行空间和时间的离散[7],用Newton-Raphson法[6]求解模型控制方程。将式(1)简化为函数形式:

式中 x为主要变量。泰勒级数展开式(3),忽略高阶项,在第m+1次迭代步、k+1次时间步:

令F(xk+1,m+1)=0,将式(4)改写为

式中:B为系数矩阵,B=∂F/∂x;u为主要变量增量,uj=xk+1,m+1-xk+1,m;F(xk+1,m)为第 k+1时间步、m迭代步的误差值。系数矩阵B用数值差分法来求,系数 bij值为式中:Δxj为一微小增量;n为主要变量数乘以单元数的值。

Newton-Raphson法具体计算步骤如下:第 k+1时间步,初始值 xk+1,0为上一时间步的解,当‖F(xk+1,m)‖2＞ε时,求 解 B(xk+1,m)u=-F(xk+1,m)[8],令 xk+1,m+1=xk+1,m+η uj(η＜1),为迭代的步长系数(减幅系数),进行下一迭代步计算,直到‖F(xk+1,m)‖2＜ε,进行下一时间步计算。

1.3 主要变量及相态转换

在给定的控制体中,相态并非固定不变,对应不同相态需采用不同的主要变量,在计算中单元的相态和主要变量需在每一迭代步后重新调整。根据Gibbsian相态定律[6],一个多组分多相流系统中自由度的数量为F=C-P+2+(P-1)=C+1,其中C为组分数,P为相数,本模型中,由于温度为常量,只需另外选择两个独立的主要变量。只有气相时,xg,w和pg为主要变量;只有液相时,xw,g和 pg为主要变量;水相气相同时存在时,Sw和pg为主要变量。相态转换的执行程序见图2。

图2 相态转换执行程序

1.4 次要变量

1.4.1摩尔分数和密度

α相组分κ的摩尔分数xα,κ满足:

其中气相中水组分的摩尔分数:

式中pw,sat为此温度下的饱和蒸汽压力。

根据Henry定律和Dalton分压定律,得到水相中空气组分的摩尔分数[9]:

式中:KH为Henry系数,KH=(0.894 2+1.47exp(-0.04394T))×10-10Pa-1;T为绝对温度。

忽略压力影响,在常温下,水相密度ρw取常数1000kg/m3。由理想气体定律可得

式中:R为气体常数;Mκ为组分κ的分子质量。

1.4.2毛细压力和相对渗透系数

收缩膜两侧承受的水压力和空气压力之差称为毛细压力,即pc=pw-pg≤0Pa,毛细压力和相对渗透系数都为饱和度的函数,采用Genuchten提出的近似公式[10]:

式中:ξ,λ为经验参数,取 ξ=1-1/λ;P0为土的进气值;Se为有效饱和度,Se=(Sw-Srw)/(Ssw-Srw),其中Srw为剩余水饱和度,Ssw为最大水饱和度。

1.5 定解条件

a.边界条件的实现:在渗流场边界外围加一圈薄层单元,称之为边界条件单元,通过赋初值给边界条件单元来设定渗流场的边界条件,假设其体积为一个极大值,如1050m3,以保证在计算过程中边界条件不发生改变。

b.稳定渗流计算的初始条件:假设边坡内部土体单元水相饱和,主要变量为 xw,g和pg,令 xw,g=10-10,pg等于标准大气压;上下游水位以上的边界条件单元气相饱和,主要变量为 xg,g和pg,令xg,g=0.9999,pg等于标准大气压;上下游水位以下的边界条件单元水相饱和,主要变量为 xw,g和pg,令 xw,g=10-10,pg等于该边界单元形心点处的水压力值。降雨情况的初始条件为稳定渗流计算的结果。

c.两种边界条件:Dirichlet边界条件,即已知气压力边界条件和已知水压力边界条件,通过给边界条件单元赋初值实现;Neumann边界条件,即流量边界条件,包括水相流量和气相流量边界,给边界条件单元施加一个对应于组分质量守恒方程(式(1))的源汇项。在稳定渗流情况下源汇项为0,在降雨情况下,水相源汇项计算公式为[5]

式中:A为与降雨方向垂直的土体表层有效面积;qt为t时刻的降雨强度。

2 土坡稳定分析公式

采用圆弧滑动面,根据简化Bishop法,假设XRXL=0,其中XR与XL分别为土条右侧和左侧的条间竖向剪力[11],已知水-气二相流模型计算出的单元孔隙压力值,要将该计算结果应用于边坡稳定分析中,就要根据水-气二相流模型网格节点的孔隙压力插值得到滑动面上各点的孔隙压力值,可采用等参单元的有限元逆变换来实现[12],本文采用Taylor级数展开的线性项进行逆变换,插值求得各土条底边中点的孔隙压力值,然后根据力矩平衡来计算边坡安全系数。

根据竖向力平衡,作用于第 i土条底面上的总法向力Ni计算公式为

其中

式中:βi为i土条底面斜向长度;c′为有效黏聚力;φ为有效内摩擦角,而参数 φb代表抗剪强度随毛细压力的增加情况,当土条底面位于饱和区内时,φb采用 φ值;αi为i土条底面中点的切线与水平面的夹角;F为安全系数;Wi为i土条重,浸润线以下按饱和密度计算;pw,i,pg,i分别为i土条底面中点的孔隙水压力和孔隙气压力,为相对压力。

由力矩平衡可得安全系数为

联合式(16)和(17),迭代求解安全系数,若不考虑土体的负孔隙水压力pw,i和负孔隙气压力pg,i,则式(16)和式(17)变为边坡稳定性分析的饱和土计算公式。

3 土坡稳定分析算例

取一均质土坡进行计算,左侧水位11m,右侧水位16m,坡度27°,底边不透水,边坡截面见图3。取两个滑动面进行研究,滑动面1为深层滑动面,大多数土条底面位于饱和区,滑动圆心为(16.00m,25.00m),滑动半径为16.16m;滑动面2为浅层滑动面,大多数土条底面位于非饱和区,滑动圆心为(20.00m,27.00m),滑动半径为13.74m。

图3 边坡剖面(单位:m)

边坡土体材料参数:饱和渗透系数ksw=1.0×10-6m/s,经验参数 ξ=0.457,进气值 P0=20kPa,孔隙率φ=0.421;剩余水饱和度Srw=0.15,最大水饱和度Ssw=1.0,有效黏聚力 c′=10.1kPa,有效内摩擦角 φ=42.6°,参数 φb=35.0°,温度 T=15℃,土体密度为1930kg/m3,土体饱和密度为2000kg/m3。

利用水-气二相流模型进行模拟,得到土坡在稳定渗流情况及降雨强度为7.2mm/h、降雨时间为5h的降雨结束瞬时的孔压及饱和度如图4和图5所示(图中孔隙气压力和孔隙水压力分别用相应的压力水头(pg-p0)/ρwg 和(pw-p0)/ρwg表示,其中 p0为大气压力。图6和图7同)。

图4 稳定渗流情况水-气二相流模型计算结果

图5 降雨情况水-气二相流模型计算结果

由图4和图5可知,边坡浸润线以下的饱和区,孔隙气压力等于孔隙水压力。稳定渗流情况下,浸润线以上的非饱和区,孔隙气压力基本为0Pa,存在负孔隙水压力;降雨后,土体表层形成暂态饱和区,非饱和区的水相饱和度增加,使得负孔隙水压力减小,由于雨水入渗,孔隙中空气尚未全部消散,空气所产生的阻力阻碍水相渗入,造成孔隙气压力明显增大,在浸润线附近,孔隙气压力形成一条基本为0Pa的过渡带。

仅考虑水相流动用有限元法计算[13],得到土坡在稳定渗流和降雨后的孔隙水压力等值线如图6和图7所示,所以通过对比水-气二相流和单相流计算的结果,可验证水-气二相流模型的正确性。

图6 单相稳定渗流情况孔隙水压力(单位:m)

图7 单相降雨情况孔隙水压力(单位:m)

根据水-气二相流模型计算出的孔压值,由式(16)和式(17)求出滑动面1和2在以下4种情况下的安全系数:①仅考虑pw＞0Pa部分,即按饱和土公式计算;②同时考虑正孔隙水压力pw＞0Pa和孔隙气压力pg;③同时考虑正负孔隙水压力pw＞0Pa和pw＜0Pa;④同时考虑正负孔隙水压力pw及孔隙气压力 pg。由水-气二相流模型求得的安全系数见表1。情况2和4与情况1相比较,安全系数的增加百分比见表2。

表1 水-气二相流模型安全系数计算结果

表2 不同计算情况下安全系数增加百分比 %

由表1可知,考虑非饱和区的孔隙气压力和负孔隙水压力(情况4)之后,边坡安全系数比按饱和土公式计算(情况1)的结果有了明显增加;稳定渗流情况下,孔隙气压力对安全系数的影响可以忽略,主要是负孔隙水压力对安全系数的增加起作用,降雨后,孔隙气压力对安全系数的影响变大,孔隙气压力起减小安全系数的作用,但与负孔隙水压力的作用相比,孔隙气压力对稳定的影响依然较小。由表2可知,滑动面与地下水位的距离越大,负孔隙水压力和孔隙气压力对边坡稳定的影响就越大。

4 结 论

a.在稳定渗流情况下,非饱和区孔隙气压力基本为0Pa,负孔隙水压力不为0Pa,安全系数的计算公式中孔隙气压力的作用可以忽略;降雨情况下,非饱和区的负孔隙水压力减小,由于空气阻力的作用,孔隙气压力增加,此时负孔隙水压力和孔隙气压力都对边坡的稳定有影响。

b.在考虑了非饱和区孔隙气压力与负孔隙水压力提供的抗剪强度之后,边坡安全系数比按饱和土公式计算的结果有了较大提高,其中孔隙气压力使安全系数减小,负孔隙水压力使安全系数增加,滑动面与地下水位距离越大,则负孔隙水压力与孔隙气压力对边坡稳定的影响越明显。

c.在边坡稳定分析时,往往忽略孔隙气压力和负孔隙水压力的作用,由本文分析可知负孔隙水压力在提高土坡稳定方面的作用不可忽视,而在降雨情况下孔隙气压力对土坡稳定的影响也应引起注意。

[1]van DIJKEM I J,van der ZEE S E A T M,van DUIJN C J.Multi-phase flow modeling of air sparging[J].Advances in Water Resources,1995,18(6):319-333.

[2]MATTHEW W F,CHRISTOPHER E K,CASS T M.Mixed finite element methods and higher order temporal approximations for variably saturated groundwater flow[J].Advances in Water Resources,2003,26:373-394.

[3]WHITE M D,OOSTR OM M,IMHARD R J.Modeling fluid flow and transport in variably saturated porous media with the STOMP simulator[J].Advances in Water Resources,1995,18(6):353-364.

[4]LAR OCHE C,VIZIKA O.Two-phase flow properties prediction from small-scale data using pore-network modeling[J].Transport in PorousMedia,2005,61:77-91.

[5]FREDLUND D G,AHARDIO H.Soilmechanics for unsaturated soils[M].New York:John&Sons,1993.

[6]CLASS H,HELMIG R,BASTIAN P.Numerical simulation of non-isothermal multiphase multicomponent processes in porous media[J].Advances in Water Resources,2002,25(5):533-550.

[7]李强,封建湖,蔡体敏,等.二相流计算的一种差分算法[J].应用数学和力学,2004,25(5):488-496.

[8]韩林,张子明,倪志强.应用改进算法的对称逐步超松弛预处理共轭梯度法进行大体积混凝土仿真计算[J].河海大学学报:自然科学版,2010,38(3):278-283.

[9]KOLDITZ O,de JONGE J.Non-isothermal two-phase flow in low-permeable porous media[J].Computational Mechanics,2004,33:345-364.

[10]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:892-898.

[11]李凯,陈国荣.基于滑移线场理论的边坡稳定性有限元分析[J].河海大学学报:自然科学版,2010,38(2):191-195.

[12]徐燕萍,项阳,刘泉声,等.等参元逆变换插值法的研究及其应用[J].岩土力学,2001,22(2):226-228

[13]黄耀英,吴中如,顾璇.基于区间参数摄动法的黏弹性区间有限元研究[J].河海大学学报:自然科学版,2008,36(5):702-706.