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

張曉悅,王 棟,張曉樂,沈躍軍

(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.