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

(1.浙江水利水電專科學(xué)校水利工程系,浙江 杭州 310018;2.信息產(chǎn)業(yè)電子第十一設(shè)計(jì)研究院科技工程股份有限公司杭州分院,浙江杭州 310020;3.浙江華東工程安全技術(shù)有限公司,浙江杭州 310014;4.浙江省錢塘江管理局勘測(cè)設(shè)計(jì)院,浙江杭州 310016)

土坡穩(wěn)定分析是估計(jì)天然及人工邊坡安全系數(shù)的一個(gè)常用工具,一般進(jìn)行邊坡穩(wěn)定分析時(shí),由于測(cè)量負(fù)孔隙水壓力和孔隙氣壓力并將其納入穩(wěn)定分析中比較困難,對(duì)于地下水位以上非飽和區(qū)中由負(fù)孔隙水壓力和孔隙氣壓力所提供的部分抗剪強(qiáng)度通常予以忽略不計(jì)。近年來,隨著人們對(duì)多相流的研究越來越多,在土坡穩(wěn)定分析中全面考慮孔隙水壓力和孔隙氣壓力成為可能。例如,Dijke等[1]給出了一個(gè)較完整的多相流模型,Matthew等[2]提出二相流模型的一種空間和時(shí)間離散方法,White等[3]分析了多孔介質(zhì)飽和度和相對(duì)滲透系數(shù)關(guān)系,Laroche等[4]對(duì)多相流模型中毛細(xì)壓力和飽和度的關(guān)系做了研究。

本文建立水-氣二相流模型來模擬土體邊坡在穩(wěn)定滲流情況和降雨情況下的水相和氣相滲流場(chǎng),全面考慮土體中的水相和氣相流動(dòng)以及水相和氣相成分的相互轉(zhuǎn)換,利用二相流計(jì)算結(jié)果求得邊坡安全系數(shù),分析孔隙氣壓力和負(fù)孔隙水壓力在土坡穩(wěn)定中的作用。

1 水-氣二相流模型

1.1 定義以及基本假設(shè)

非飽和土是四相混合體,包括固相、氣相、液相以及稱為收縮膜的水氣分界面[5]。在土體中可能有多種液相,本文只考慮水相。非飽和土中的滲流為包括氣相和液相流動(dòng)的二相流,氣相和液相都包括水和空氣兩種組分。液相和氣相通過蒸發(fā)與溶解相互轉(zhuǎn)換[5],相態(tài)及其相互轉(zhuǎn)換關(guān)系見圖1。

圖1 相態(tài)及相態(tài)轉(zhuǎn)換

模型假設(shè)土體骨架不變形,即孔隙率為常數(shù),孔隙流體不可壓縮,忽略化學(xué)生物學(xué)反應(yīng),恒溫,不考慮本構(gòu)關(guān)系中滲透系數(shù)與飽和度之間的滯后作用,氣相中所有組分都遵守理想氣體定律,氣相中水組分分壓力等于此溫度下的飽和蒸汽壓力。

1.2 模型控制方程及其求解

模型考慮氣相和液相的流動(dòng)以及各組分之間的質(zhì)量傳遞,模型的控制方程即組分κ的質(zhì)量守恒方程[6]:

土體固有滲透系數(shù):

式中:下標(biāo) α和κ分別代表相態(tài)和組分,取g時(shí)為氣相或空氣組分,取w時(shí)為水相或水組分;φ為土體孔隙率;Sα為α相的飽和度;ρm,α=nα/Vα,nα為 α相的摩爾質(zhì)量,Vα為nα對(duì)應(yīng)的體積,m3;ρα為α相的密度;xα,κ為α相中組分κ的摩爾分?jǐn)?shù);t為時(shí)間;ksw為土的飽和滲透系數(shù);krα為α相的相對(duì)滲透系數(shù);pα為α相的壓力;μα為α相的絕對(duì)黏滯性,溫度一定的情況下為常數(shù);g為重力加速度;Qκ為κ組分的源匯項(xiàng),kg/s。

控制方程的變量分為主要變量和次要變量?jī)深?次要變量為主要變量的函數(shù)。采用有限差分法對(duì)控制方程進(jìn)行空間和時(shí)間的離散[7],用Newton-Raphson法[6]求解模型控制方程。將式(1)簡(jiǎn)化為函數(shù)形式:

式中 x為主要變量。泰勒級(jí)數(shù)展開式(3),忽略高階項(xiàng),在第m+1次迭代步、k+1次時(shí)間步:

令F(xk+1,m+1)=0,將式(4)改寫為

式中:B為系數(shù)矩陣,B=?F/?x;u為主要變量增量,uj=xk+1,m+1-xk+1,m;F(xk+1,m)為第 k+1時(shí)間步、m迭代步的誤差值。系數(shù)矩陣B用數(shù)值差分法來求,系數(shù) bij值為式中:Δxj為一微小增量;n為主要變量數(shù)乘以單元數(shù)的值。

Newton-Raphson法具體計(jì)算步驟如下:第 k+1時(shí)間步,初始值 xk+1,0為上一時(shí)間步的解,當(dāng)‖F(xiàn)(xk+1,m)‖2＞ε時(shí),求 解 B(xk+1,m)u=-F(xk+1,m)[8],令 xk+1,m+1=xk+1,m+η uj(η＜1),為迭代的步長(zhǎng)系數(shù)(減幅系數(shù)),進(jìn)行下一迭代步計(jì)算,直到‖F(xiàn)(xk+1,m)‖2＜ε,進(jìn)行下一時(shí)間步計(jì)算。

1.3 主要變量及相態(tài)轉(zhuǎn)換

在給定的控制體中,相態(tài)并非固定不變,對(duì)應(yīng)不同相態(tài)需采用不同的主要變量,在計(jì)算中單元的相態(tài)和主要變量需在每一迭代步后重新調(diào)整。根據(jù)Gibbsian相態(tài)定律[6],一個(gè)多組分多相流系統(tǒng)中自由度的數(shù)量為F=C-P+2+(P-1)=C+1,其中C為組分?jǐn)?shù),P為相數(shù),本模型中,由于溫度為常量,只需另外選擇兩個(gè)獨(dú)立的主要變量。只有氣相時(shí),xg,w和pg為主要變量;只有液相時(shí),xw,g和 pg為主要變量;水相氣相同時(shí)存在時(shí),Sw和pg為主要變量。相態(tài)轉(zhuǎn)換的執(zhí)行程序見圖2。

圖2 相態(tài)轉(zhuǎn)換執(zhí)行程序

1.4 次要變量

1.4.1摩爾分?jǐn)?shù)和密度

α相組分κ的摩爾分?jǐn)?shù)xα,κ滿足:

其中氣相中水組分的摩爾分?jǐn)?shù):

式中pw,sat為此溫度下的飽和蒸汽壓力。

根據(jù)Henry定律和Dalton分壓定律,得到水相中空氣組分的摩爾分?jǐn)?shù)[9]:

式中:KH為Henry系數(shù),KH=(0.894 2+1.47exp(-0.04394T))×10-10Pa-1;T為絕對(duì)溫度。

忽略壓力影響,在常溫下,水相密度ρw取常數(shù)1000kg/m3。由理想氣體定律可得

式中:R為氣體常數(shù);Mκ為組分κ的分子質(zhì)量。

1.4.2毛細(xì)壓力和相對(duì)滲透系數(shù)

收縮膜兩側(cè)承受的水壓力和空氣壓力之差稱為毛細(xì)壓力,即pc=pw-pg≤0Pa,毛細(xì)壓力和相對(duì)滲透系數(shù)都為飽和度的函數(shù),采用Genuchten提出的近似公式[10]:

式中:ξ,λ為經(jīng)驗(yàn)參數(shù),取 ξ=1-1/λ;P0為土的進(jìn)氣值;Se為有效飽和度,Se=(Sw-Srw)/(Ssw-Srw),其中Srw為剩余水飽和度,Ssw為最大水飽和度。

1.5 定解條件

a.邊界條件的實(shí)現(xiàn):在滲流場(chǎng)邊界外圍加一圈薄層單元,稱之為邊界條件單元,通過賦初值給邊界條件單元來設(shè)定滲流場(chǎng)的邊界條件,假設(shè)其體積為一個(gè)極大值,如1050m3,以保證在計(jì)算過程中邊界條件不發(fā)生改變。

b.穩(wěn)定滲流計(jì)算的初始條件:假設(shè)邊坡內(nèi)部土體單元水相飽和,主要變量為 xw,g和pg,令 xw,g=10-10,pg等于標(biāo)準(zhǔn)大氣壓;上下游水位以上的邊界條件單元?dú)庀囡柡?主要變量為 xg,g和pg,令xg,g=0.9999,pg等于標(biāo)準(zhǔn)大氣壓;上下游水位以下的邊界條件單元水相飽和,主要變量為 xw,g和pg,令 xw,g=10-10,pg等于該邊界單元形心點(diǎn)處的水壓力值。降雨情況的初始條件為穩(wěn)定滲流計(jì)算的結(jié)果。

c.兩種邊界條件:Dirichlet邊界條件,即已知?dú)鈮毫吔鐥l件和已知水壓力邊界條件,通過給邊界條件單元賦初值實(shí)現(xiàn);Neumann邊界條件,即流量邊界條件,包括水相流量和氣相流量邊界,給邊界條件單元施加一個(gè)對(duì)應(yīng)于組分質(zhì)量守恒方程(式(1))的源匯項(xiàng)。在穩(wěn)定滲流情況下源匯項(xiàng)為0,在降雨情況下,水相源匯項(xiàng)計(jì)算公式為[5]

式中:A為與降雨方向垂直的土體表層有效面積;qt為t時(shí)刻的降雨強(qiáng)度。

2 土坡穩(wěn)定分析公式

采用圓弧滑動(dòng)面,根據(jù)簡(jiǎn)化Bishop法,假設(shè)XRXL=0,其中XR與XL分別為土條右側(cè)和左側(cè)的條間豎向剪力[11],已知水-氣二相流模型計(jì)算出的單元孔隙壓力值,要將該計(jì)算結(jié)果應(yīng)用于邊坡穩(wěn)定分析中,就要根據(jù)水-氣二相流模型網(wǎng)格節(jié)點(diǎn)的孔隙壓力插值得到滑動(dòng)面上各點(diǎn)的孔隙壓力值,可采用等參單元的有限元逆變換來實(shí)現(xiàn)[12],本文采用Taylor級(jí)數(shù)展開的線性項(xiàng)進(jìn)行逆變換,插值求得各土條底邊中點(diǎn)的孔隙壓力值,然后根據(jù)力矩平衡來計(jì)算邊坡安全系數(shù)。

根據(jù)豎向力平衡,作用于第 i土條底面上的總法向力Ni計(jì)算公式為

其中

式中:βi為i土條底面斜向長(zhǎng)度;c′為有效黏聚力;φ為有效內(nèi)摩擦角,而參數(shù) φb代表抗剪強(qiáng)度隨毛細(xì)壓力的增加情況,當(dāng)土條底面位于飽和區(qū)內(nèi)時(shí),φb采用 φ值;αi為i土條底面中點(diǎn)的切線與水平面的夾角;F為安全系數(shù);Wi為i土條重,浸潤(rùn)線以下按飽和密度計(jì)算;pw,i,pg,i分別為i土條底面中點(diǎn)的孔隙水壓力和孔隙氣壓力,為相對(duì)壓力。

由力矩平衡可得安全系數(shù)為

聯(lián)合式(16)和(17),迭代求解安全系數(shù),若不考慮土體的負(fù)孔隙水壓力pw,i和負(fù)孔隙氣壓力pg,i,則式(16)和式(17)變?yōu)檫吰路€(wěn)定性分析的飽和土計(jì)算公式。

3 土坡穩(wěn)定分析算例

取一均質(zhì)土坡進(jìn)行計(jì)算,左側(cè)水位11m,右側(cè)水位16m,坡度27°,底邊不透水,邊坡截面見圖3。取兩個(gè)滑動(dòng)面進(jìn)行研究,滑動(dòng)面1為深層滑動(dòng)面,大多數(shù)土條底面位于飽和區(qū),滑動(dòng)圓心為(16.00m,25.00m),滑動(dòng)半徑為16.16m;滑動(dòng)面2為淺層滑動(dòng)面,大多數(shù)土條底面位于非飽和區(qū),滑動(dòng)圓心為(20.00m,27.00m),滑動(dòng)半徑為13.74m。

圖3 邊坡剖面(單位:m)

邊坡土體材料參數(shù):飽和滲透系數(shù)ksw=1.0×10-6m/s,經(jīng)驗(yàn)參數(shù) ξ=0.457,進(jìn)氣值 P0=20kPa,孔隙率φ=0.421;剩余水飽和度Srw=0.15,最大水飽和度Ssw=1.0,有效黏聚力 c′=10.1kPa,有效內(nèi)摩擦角 φ=42.6°,參數(shù) φb=35.0°,溫度 T=15℃,土體密度為1930kg/m3,土體飽和密度為2000kg/m3。

利用水-氣二相流模型進(jìn)行模擬,得到土坡在穩(wěn)定滲流情況及降雨強(qiáng)度為7.2mm/h、降雨時(shí)間為5h的降雨結(jié)束瞬時(shí)的孔壓及飽和度如圖4和圖5所示(圖中孔隙氣壓力和孔隙水壓力分別用相應(yīng)的壓力水頭(pg-p0)/ρwg 和(pw-p0)/ρwg表示,其中 p0為大氣壓力。圖6和圖7同)。

圖4 穩(wěn)定滲流情況水-氣二相流模型計(jì)算結(jié)果

圖5 降雨情況水-氣二相流模型計(jì)算結(jié)果

由圖4和圖5可知,邊坡浸潤(rùn)線以下的飽和區(qū),孔隙氣壓力等于孔隙水壓力。穩(wěn)定滲流情況下,浸潤(rùn)線以上的非飽和區(qū),孔隙氣壓力基本為0Pa,存在負(fù)孔隙水壓力;降雨后,土體表層形成暫態(tài)飽和區(qū),非飽和區(qū)的水相飽和度增加,使得負(fù)孔隙水壓力減小,由于雨水入滲,孔隙中空氣尚未全部消散,空氣所產(chǎn)生的阻力阻礙水相滲入,造成孔隙氣壓力明顯增大,在浸潤(rùn)線附近,孔隙氣壓力形成一條基本為0Pa的過渡帶。

僅考慮水相流動(dòng)用有限元法計(jì)算[13],得到土坡在穩(wěn)定滲流和降雨后的孔隙水壓力等值線如圖6和圖7所示,所以通過對(duì)比水-氣二相流和單相流計(jì)算的結(jié)果,可驗(yàn)證水-氣二相流模型的正確性。

圖6 單相穩(wěn)定滲流情況孔隙水壓力(單位:m)

圖7 單相降雨情況孔隙水壓力(單位:m)

根據(jù)水-氣二相流模型計(jì)算出的孔壓值,由式(16)和式(17)求出滑動(dòng)面1和2在以下4種情況下的安全系數(shù):①僅考慮pw＞0Pa部分,即按飽和土公式計(jì)算;②同時(shí)考慮正孔隙水壓力pw＞0Pa和孔隙氣壓力pg;③同時(shí)考慮正負(fù)孔隙水壓力pw＞0Pa和pw＜0Pa;④同時(shí)考慮正負(fù)孔隙水壓力pw及孔隙氣壓力 pg。由水-氣二相流模型求得的安全系數(shù)見表1。情況2和4與情況1相比較,安全系數(shù)的增加百分比見表2。

表1 水-氣二相流模型安全系數(shù)計(jì)算結(jié)果

表2 不同計(jì)算情況下安全系數(shù)增加百分比 %

由表1可知,考慮非飽和區(qū)的孔隙氣壓力和負(fù)孔隙水壓力(情況4)之后,邊坡安全系數(shù)比按飽和土公式計(jì)算(情況1)的結(jié)果有了明顯增加;穩(wěn)定滲流情況下,孔隙氣壓力對(duì)安全系數(shù)的影響可以忽略,主要是負(fù)孔隙水壓力對(duì)安全系數(shù)的增加起作用,降雨后,孔隙氣壓力對(duì)安全系數(shù)的影響變大,孔隙氣壓力起減小安全系數(shù)的作用,但與負(fù)孔隙水壓力的作用相比,孔隙氣壓力對(duì)穩(wěn)定的影響依然較小。由表2可知,滑動(dòng)面與地下水位的距離越大,負(fù)孔隙水壓力和孔隙氣壓力對(duì)邊坡穩(wěn)定的影響就越大。

4 結(jié) 論

a.在穩(wěn)定滲流情況下,非飽和區(qū)孔隙氣壓力基本為0Pa,負(fù)孔隙水壓力不為0Pa,安全系數(shù)的計(jì)算公式中孔隙氣壓力的作用可以忽略;降雨情況下,非飽和區(qū)的負(fù)孔隙水壓力減小,由于空氣阻力的作用,孔隙氣壓力增加,此時(shí)負(fù)孔隙水壓力和孔隙氣壓力都對(duì)邊坡的穩(wěn)定有影響。

b.在考慮了非飽和區(qū)孔隙氣壓力與負(fù)孔隙水壓力提供的抗剪強(qiáng)度之后,邊坡安全系數(shù)比按飽和土公式計(jì)算的結(jié)果有了較大提高,其中孔隙氣壓力使安全系數(shù)減小,負(fù)孔隙水壓力使安全系數(shù)增加,滑動(dòng)面與地下水位距離越大,則負(fù)孔隙水壓力與孔隙氣壓力對(duì)邊坡穩(wěn)定的影響越明顯。

c.在邊坡穩(wěn)定分析時(shí),往往忽略孔隙氣壓力和負(fù)孔隙水壓力的作用,由本文分析可知負(fù)孔隙水壓力在提高土坡穩(wěn)定方面的作用不可忽視,而在降雨情況下孔隙氣壓力對(duì)土坡穩(wěn)定的影響也應(yīng)引起注意。

[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]李強(qiáng),封建湖,蔡體敏,等.二相流計(jì)算的一種差分算法[J].應(yīng)用數(shù)學(xué)和力學(xué),2004,25(5):488-496.

[8]韓林,張子明,倪志強(qiáng).應(yīng)用改進(jìn)算法的對(duì)稱逐步超松弛預(yù)處理共軛梯度法進(jìn)行大體積混凝土仿真計(jì)算[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,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]李凱,陳國(guó)榮.基于滑移線場(chǎng)理論的邊坡穩(wěn)定性有限元分析[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2010,38(2):191-195.

[12]徐燕萍,項(xiàng)陽,劉泉聲,等.等參元逆變換插值法的研究及其應(yīng)用[J].巖土力學(xué),2001,22(2):226-228

[13]黃耀英,吳中如,顧璇.基于區(qū)間參數(shù)攝動(dòng)法的黏彈性區(qū)間有限元研究[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2008,36(5):702-706.