基于水氣二相流的邊坡降雨入滲有限元分析

習(xí)念念，童富果，劉 剛，程慶超

(1.三峽大學(xué) 水利與環(huán)境學(xué)院，湖北 宜昌 443002； 2.湖北省水利水電規(guī)劃勘測(cè)設(shè)計(jì)院，武漢 430064)

?

基于水氣二相流的邊坡降雨入滲有限元分析

習(xí)念念1，2，童富果1，劉 剛1，程慶超1

(1.三峽大學(xué) 水利與環(huán)境學(xué)院，湖北 宜昌 443002； 2.湖北省水利水電規(guī)劃勘測(cè)設(shè)計(jì)院，武漢 430064)

為揭示邊坡降雨入滲基本規(guī)律，基于水氣二相流理論，以三峽庫(kù)區(qū)大坪滑坡為例，采用有限單元法分析邊坡降雨入滲水氣運(yùn)動(dòng)規(guī)律。結(jié)果表明：根據(jù)水氣流入、流出的方向不同，坡表可分為吸入?yún)^(qū)和溢出區(qū)，通常坡體上部為水、氣吸入?yún)^(qū)，下部為水、氣溢出區(qū)；水、氣進(jìn)出坡表的總速率受降雨特性、坡面產(chǎn)流情況、氣溫變化等外在因素的影響較小，具有相對(duì)穩(wěn)定性；在坡表吸入?yún)^(qū)，當(dāng)降雨強(qiáng)度大于水氣入滲總速率時(shí)，入滲強(qiáng)度等于水氣入滲總速率，吸入?yún)^(qū)產(chǎn)流；反之，入滲強(qiáng)度等于降雨強(qiáng)度，吸入?yún)^(qū)不產(chǎn)流。此外，坡表飽和區(qū)水壓力可通過(guò)孔隙氣傳遞到地下水面，增大地下水的壓力水頭，不利于坡體穩(wěn)定。

水氣二相流；降雨入滲；有限單元法；非飽和滲流；三峽庫(kù)區(qū)邊坡

1 研究背景

邊坡降雨入滲是一個(gè)復(fù)雜的非飽和滲流過(guò)程，其入滲強(qiáng)度主要取決于水、氣兩相流體在土體中的耦合驅(qū)動(dòng)效應(yīng)[1]，并與土體滲透性、水氣黏滯性以及坡面邊界條件等因素有關(guān)。以往研究[2-4]大多基于Richard方程來(lái)描述降雨入滲，因其假定孔隙氣壓恒定且處處相等(通常設(shè)為大氣壓)，實(shí)質(zhì)上忽略了水、氣耦合作用，僅考慮了水的運(yùn)動(dòng)，在特定情況下可能會(huì)與工程實(shí)際情況出入較大。大量研究成果表明，在邊界封閉條件下，孔隙氣對(duì)液相(水)流動(dòng)的影響不可忽略[1,5-6]。雨水入滲過(guò)程中，土體中空氣受液相(水)的擠壓而產(chǎn)生與水體運(yùn)動(dòng)方向相反的頂托力，進(jìn)而雨水入滲速率會(huì)因水壓力梯度的減小而降低。唐海行等[7]，彭勝等[8]，孫冬梅等[9]學(xué)者考慮氣相影響，實(shí)現(xiàn)了邊坡降雨的水、氣入滲過(guò)程，認(rèn)為采用水氣二相流理論探究邊坡降雨入滲問(wèn)題更符合實(shí)際。目前國(guó)內(nèi)關(guān)于二相流的研究主要探究入滲過(guò)程中壓力等參數(shù)變化，驗(yàn)證了氣相的影響，較少涉及降雨入滲規(guī)律的一般歸納。本文基于水氣二相流的基本理論及方法，對(duì)降雨入滲條件下的三峽庫(kù)區(qū)大坪滑坡進(jìn)行模擬分析，以研究邊坡降雨入滲的一般規(guī)律，是對(duì)現(xiàn)有降雨入滲規(guī)律的補(bǔ)充，為工程實(shí)例提供更多依據(jù)。

2 水氣二相流控制微分方程

水氣二相流控制微分方程包含液相(水)流動(dòng)方程和氣相流動(dòng)方程，兩者通過(guò)飽和度、孔壓、基質(zhì)吸力、孔隙率等參變量的關(guān)聯(lián)性來(lái)實(shí)現(xiàn)水氣二相流耦合描述。水、氣在土體孔隙中的流動(dòng)需各自滿足壓力驅(qū)動(dòng)下的質(zhì)量守恒方程，同時(shí)考慮土體變形、水氣的可壓縮性、黏滯性以及飽和度變化等因素，即可得水氣二相流控制微分方程[10]

(1)

(2)

水氣二相流控制微分方程是基于時(shí)間和空間的非線性偏微分方程組，求解時(shí)需先進(jìn)行時(shí)間和空間離散，然后通過(guò)迭代求解。本計(jì)算空間離散化采用Garlerkin有限單元法，時(shí)間離散采用差分法。計(jì)算時(shí)孔隙氣壓pg和飽和度S為未知量，通過(guò)對(duì)方程(1)、方程(2)循環(huán)計(jì)算以實(shí)現(xiàn)迭代求解。

3 數(shù)值計(jì)算模型

3.1 工程概況及幾何模型

本研究以三峽庫(kù)區(qū)大坪滑坡為案例進(jìn)行計(jì)算，大坪滑坡位于巴東縣境內(nèi)長(zhǎng)江北岸，上距巴東縣城3km，下距三峽工程壩址61km。該滑坡體物質(zhì)主要為松散巖石碎屑層和黃土層，透水性好，有利于地下水的下滲，滑坡體底部為千枚巖，透水性相對(duì)較弱，容易形成相對(duì)隔水層。同時(shí)滑體物質(zhì)厚度較大，為滑坡提供了物質(zhì)基礎(chǔ)。大坪滑坡體前緣寬1 100m，后緣寬800m，縱向(順坡向)長(zhǎng)550m，后緣分布在高程280～338m一線，總體呈弧形。取其主滑動(dòng)剖面進(jìn)行數(shù)值建模，剖面水平長(zhǎng)度為600m，豎直向高度為350m。為較好地模擬降雨入滲，網(wǎng)格剖分時(shí)滑坡體表層采用厚度較小單元，其厚度沿坡體深度方向按0.25，0.5，0.75m逐漸增加。剖分所得有限元計(jì)算網(wǎng)格如圖1所示，共計(jì)3 499個(gè)單元，3 427個(gè)節(jié)點(diǎn)。

圖1 大坪滑坡有限元分析網(wǎng)格

3.2 初始及邊界條件

為盡可能減小初始狀態(tài)不確定性對(duì)研究帶來(lái)的影響，本分析結(jié)合三峽庫(kù)區(qū)多年降雨實(shí)測(cè)資料(由國(guó)土資源部三峽庫(kù)區(qū)地質(zhì)災(zāi)害防治工作指揮部提供)，以坡體全飽和態(tài)起算，計(jì)算模擬至坡體滲流場(chǎng)相對(duì)穩(wěn)定，取其穩(wěn)定后的滲流場(chǎng)為計(jì)算的初始條件。

模型后緣側(cè)面及底部為不透水邊界，該邊界水位及其變化情況可由計(jì)算得到。庫(kù)水位以下部分坡表面及前緣側(cè)面為已知水壓邊界，其數(shù)值大小取決于庫(kù)水位高程(庫(kù)水位高程為145m)。庫(kù)水位以上坡表為已知?dú)鈮毫吔?，鑒于坡面徑流水頭相較于大氣壓力相對(duì)很小，可忽略，本計(jì)算中氣壓力邊界上的氣壓力均設(shè)為大氣壓。與傳統(tǒng)單相流計(jì)算相比，坡表水、氣入滲率可由計(jì)算結(jié)果得出，不需要作為流量邊界而事先給定。

3.3 主要本構(gòu)關(guān)系和參數(shù)

本文土水特征曲線采用常用的VanGenuchten(1980)模型[11]，基質(zhì)吸力(pc≡pg-pl)與飽和度的關(guān)系可表示為

pc=-p0(Se-1/λ-1)1-λ。

(3)

(4)

(5)

式(1)、式(2)中涉及的其他主要參數(shù)及其數(shù)值分別為n=0.15，k=4.0×10-12m2，g=9.8N/kg，ρg=1.29kg/m3，ρl=1×103kg/m3，μl=1.0×10-3N·s/m2，μg=1.0×10-5N·s/m2。

4 邊坡降雨入滲分析

降雨入滲是一個(gè)非飽和滲流過(guò)程，基于水氣二相流有限元計(jì)算可獲得坡體孔隙氣壓、水壓、水飽和度等隨時(shí)間空間的變化情況，基于分析需要，可繪制任意時(shí)刻等值線分布圖(如圖2)?？紤]氣相影響后，降雨入滲過(guò)程并非傳統(tǒng)中的均勻下滲，降雨初期，孔隙氣壓較小，降雨入滲強(qiáng)度相對(duì)較大，隨著降雨歷時(shí)增加，淺層坡體含水量增加，坡表逐漸飽和，降雨入滲強(qiáng)度會(huì)隨孔隙氣壓增大而逐漸減小。地下水位會(huì)受降雨強(qiáng)度、降雨歷時(shí)、降雨間隔等因素影響，并在時(shí)間上具有一定的滯后性。坡度變幅不大的情況下，飽和區(qū)的孔隙水壓力與距地表深度正相關(guān)，非飽和區(qū)的孔隙水壓力受初始條件、邊界條件等因素影響與距地表深度之間關(guān)系不顯著。

圖2 大坪滑坡某時(shí)刻等值線分布

圖3 坡表(水、氣)吸入、溢出區(qū)

結(jié)合大坪滑坡水氣兩相流有限元計(jì)算結(jié)果，本文先從坡表降雨入滲特性、降雨入滲強(qiáng)度和產(chǎn)流條件等方面闡述邊坡降雨入滲一般規(guī)律，然后就邊坡內(nèi)水氣耦合傳力機(jī)制進(jìn)行說(shuō)明。

4.1 坡表降雨入滲特性

根據(jù)大坪滑坡計(jì)算結(jié)果，物質(zhì)(水、氣)進(jìn)出坡表的方向和總速率呈現(xiàn)一定規(guī)律性。按照物質(zhì)(水、氣)流入、流出方向不同，坡體表面可分為吸入?yún)^(qū)和溢出區(qū)，水氣流速方向大體朝坡內(nèi)則為吸入?yún)^(qū)，反之則為溢出區(qū)(見圖3)。分析表明，坡表吸入?yún)^(qū)和溢出區(qū)分布具有較強(qiáng)的規(guī)律性，水、氣吸入?yún)^(qū)通常位于滑坡體中上部位，而溢出區(qū)一般位于其下部。吸入?yún)^(qū)和溢出區(qū)的分界線通常不是固定的，除與坡體幾何特性、滲透特性有關(guān)外，還會(huì)受降雨歷時(shí)、降雨間隔、降雨強(qiáng)度等降雨特性的影響。隨降雨持續(xù)時(shí)間的延長(zhǎng)，吸入?yún)^(qū)和溢出區(qū)的分界線會(huì)上移，坡表吸入?yún)^(qū)面積會(huì)縮小，溢出區(qū)面積會(huì)增加；隨降雨間隔的增加，分界線會(huì)逐漸下移，坡表吸入?yún)^(qū)面積會(huì)增大，而溢出區(qū)面積會(huì)逐漸減小。

從數(shù)值上分析，水、氣進(jìn)出坡表的總速率具有相對(duì)穩(wěn)定性，其大小主要取決于坡體幾何形態(tài)尺寸和巖土材料傳輸特性，受降雨特性、坡面產(chǎn)流情況、氣溫變化等外在因素的影響相對(duì)較小。水氣進(jìn)出坡表的速率可描述降雨吸入、溢出強(qiáng)度，圖4為大坪滑坡坡表水、氣吸入(溢出)強(qiáng)度與時(shí)間關(guān)系。曲線顯示，隨降雨歷時(shí)增加，大坪滑坡坡面最大吸入、溢出強(qiáng)度均具有一定的穩(wěn)定性，最大吸入強(qiáng)度約為0.4mm/min，最大溢出強(qiáng)度約為0.5mm/min。坡面水、氣平均吸入和溢出強(qiáng)度值總體均呈平緩下降趨勢(shì)，兩者在數(shù)值大小上非常接近，均約為0.1mm/min。

圖4 坡面水氣吸入、溢出強(qiáng)度時(shí)間關(guān)系曲線

4.2 坡表降雨入滲強(qiáng)度和產(chǎn)流條件

坡表降雨入滲強(qiáng)度和產(chǎn)流條件的判斷與降雨強(qiáng)度和坡面(水、氣)總吸入速率的相對(duì)大小關(guān)系密切。在坡表吸入?yún)^(qū)，水氣入滲量主要取決于坡面邊界水、氣的含量。在地表無(wú)降雨的情況下，吸入?yún)^(qū)邊界吸入的是空氣，地表有降雨情況下，降雨強(qiáng)度小于坡面(水、氣)總吸入速率時(shí)吸入?yún)^(qū)坡面不產(chǎn)流，雨水入滲強(qiáng)度等于降雨強(qiáng)度，吸入?yún)^(qū)邊界將吸入不同含量的水、氣混合物。降雨強(qiáng)度大于坡面(水、氣)總吸入速率時(shí)，坡面會(huì)有徑流產(chǎn)生，入滲強(qiáng)度等于水、氣總吸入強(qiáng)度，坡表氣封，吸入?yún)^(qū)邊界吸入的是水。大坪滑坡案例中，降雨強(qiáng)度大于坡面總吸入強(qiáng)度(約為0.1mm/min)將是其坡面產(chǎn)流的必要條件，由于最大吸入?yún)^(qū)位于坡頂，故當(dāng)降雨強(qiáng)度大于坡表最大吸入強(qiáng)度(0.4mm/min)時(shí)，坡面才會(huì)有全面產(chǎn)流的可能。

在坡表溢出區(qū)，水氣溢出量除與坡面邊界水、氣的含量有關(guān)，還與坡體飽和度緊密相關(guān)。在飽和狀態(tài)下，降雨入滲強(qiáng)度為0，坡表產(chǎn)流，溢出區(qū)邊界流出的是水；非飽和狀態(tài)下，降雨入滲強(qiáng)度等于降雨強(qiáng)度，坡表不產(chǎn)流，溢出區(qū)邊界流出的是水、氣混合物，若在完全干燥的情況下，溢出區(qū)邊界溢出的是孔隙氣。

4.3 水氣耦合傳力機(jī)制

圖5 水氣耦合傳力機(jī)制示意圖

降雨入滲是水、氣兩相流體在土體孔隙中相互驅(qū)動(dòng)耦合作用的過(guò)程，受重力的作用水體向下運(yùn)動(dòng)，進(jìn)而驅(qū)動(dòng)孔隙氣也向下運(yùn)動(dòng)，空氣因?yàn)閿D壓產(chǎn)生與水體運(yùn)動(dòng)方向相反的頂托力。由于空氣的黏滯性非常小(約為水的1%)，坡表降雨飽和區(qū)水壓可通過(guò)孔隙氣快速傳遞到地下水面，使地下水面承壓(如圖5)，表現(xiàn)了氣體作用對(duì)邊坡穩(wěn)定的不利影響。圖6是大坪滑坡某時(shí)刻的飽和度分布云圖，由圖可見，坡體頂部的飽和水壓力會(huì)經(jīng)非飽和區(qū)內(nèi)封閉的孔隙氣傳到滑坡體前端飽和區(qū)，局部抬高坡體前沿阻滑區(qū)的地下水位，形成不利坡面穩(wěn)定的外法向推力。

圖6 大坪滑坡飽和度分布云圖(傳力示意)

5 結(jié) 論

本文基于水氣二相流理論，采用有限元法計(jì)算分析了大坪滑坡體降雨入滲過(guò)程。結(jié)果表明：①根據(jù)水、氣進(jìn)出坡體的方向不同，坡體表面可分為吸入?yún)^(qū)和溢出區(qū)，通常吸入?yún)^(qū)位于坡體表面中上部，溢出區(qū)位于其中下部；②水、氣進(jìn)出坡表的總速率主要取決于坡體材料傳輸特性和幾何特性，受坡面產(chǎn)流情況和徑流深度影響較小，具有一定的穩(wěn)定性；③降雨強(qiáng)度和水、氣總吸入速率的相對(duì)大小關(guān)系可作為坡表降雨入滲強(qiáng)度和產(chǎn)流條件的判斷條件。在坡表吸入?yún)^(qū)，當(dāng)降雨強(qiáng)度小于坡面水、氣總吸入速率時(shí)不產(chǎn)流，入滲強(qiáng)度等于降雨強(qiáng)度；反之，坡面產(chǎn)流，入滲強(qiáng)度等于水、氣總吸入強(qiáng)度；④坡體表層飽和區(qū)水壓力可通過(guò)非飽和區(qū)的孔隙氣傳遞到地下飽和區(qū)，使地下水面承壓，增大了坡腳外法向推力，不利邊坡穩(wěn)定。

本計(jì)算通過(guò)模擬分析降雨入滲一般規(guī)律，為工程實(shí)例提供了更多理論依據(jù)。因影響降雨入滲過(guò)程因素眾多，故對(duì)入滲規(guī)律的定量分析較少，擬在后續(xù)研究中，分析控制各影響因素，深入量化入滲過(guò)程的一般規(guī)律。

[1] 徐 晗，朱以文，蔡元奇，等.降雨入滲條件下非飽和土邊坡穩(wěn)定分析[J].巖土力學(xué)， 2005,(12):1957-1962.

[2] 周家文，徐衛(wèi)亞，鄧俊曄，等.降雨入滲條件下邊坡的穩(wěn)定性分析[J].水利學(xué)報(bào)，2008，39(9):1066-1073.

[3] 榮 冠，張 偉，周創(chuàng)兵，等.降雨入滲條件下邊坡巖體飽和非飽和滲流計(jì)算[J].巖土力學(xué)，2005，26(l0):1544-1550.

[4] 朱 偉，陳學(xué)東，鐘小春.降雨入滲規(guī)律的實(shí)測(cè)與分析[J]. 巖土力學(xué)， 2006， 27(11): 1873-1879.

[5]LALOUIL,KLUBERTANZG,GULLIETL.Solid-liquid-airCouplinginMultiphasePorousMedia[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics, 2003,27:183-206.

[6] 王葉嬌， 曹 玲， 徐永福. 降雨入滲下非飽和土邊坡臨界穩(wěn)定性分析[J].長(zhǎng)江科學(xué)院院報(bào)， 2013， 30(9): 75-79.

[7] 唐海行，蘇逸深.考慮氣壓勢(shì)影響的降雨入滲數(shù)值模擬研究[J].水科學(xué)進(jìn)展， 1996,(7):8-13.

[8] 彭 勝，陳家軍，王金生，等.包氣帶水氣二相流國(guó)外研究綜述[J].水科學(xué)進(jìn)展， 2000,(3):333-338.

[9] 孫冬梅，朱岳明，張明進(jìn).降雨入滲過(guò)程的水-氣二相流模型研究[J].水利學(xué)報(bào)， 2007,(2):150-156.

[10]TONGFG，JINGL，ZIMMERMANRW.AFully-coupledFiniteElementCodeforModelingThermo-hydro-mechanicalProcessesinPorousGeologicalMedia[C]∥Proceedingsofthe43rdUSRockMechanicsSymposium,Asheville,NC,USA,June28-30, 2009: 28-30.

[11]VANGENUCHTENMTH.AClosedFormEquationforPredictingtheHydraulicconductivityofUnsaturatedSoil[J].SoilScienceSocietyofAmericaJournal,1980, 44(5): 892-898.

[12]MUALEMY.ANewModelforPredictingtheHydraulicConductivityofUnsaturatedPorousMedia[J].WaterResourcesResearch, 1976, 12(3): 513-522.

[13]COREYAT.TheInterrelationBetweenGasandOilRelativePermeabilities[J].ProducersMonthly, 1954, 19: 38-41.

(編輯：劉運(yùn)飛)

Finite Element Analysis of Water-air Two-Phase Infiltration in Slope

XI Nian-nian1，2, TONG Fu-guo1, LIU Gang1, CHENG Qing-chao1

(1.College of Hydraulic and Environmental Engineering, China Three Gorges University,Yichang 443002, China； 2.Hubei Provincial Water Resources and Hydropower Planning Survey and Design Institute, Wuhan 430064,China)

On the basis of the theory of two-phase (water and air) flow, the law of water-air movement caused by rainfall infiltration in slope was revealed by means of finite element method. Daping landslide in Three Gorges Reservoir was taken as a case study. The simulation results show that the slope surface can be divided into the inflow area and the outflow area according to the inflow and outflow of water and air. Usually the upper part of slope is inflow area, and the lower part is outflow area. The total infiltration rate of water and air in slope is relatively stable, slightly affected by rainfall characteristics, slope surface runoff and temperature changes. In the inflow area, when rainfall intensity is greater than infiltration rate, water in slope will run off, otherwise water or air will completely infiltrate. The saturation pressure of slope surface could transmit to groundwater surface through pore air, which increases the groundwater pressure head, unfavorable for slope stability.

water-air two-phase flow; rainfall infiltration; finite element method; unsaturated seepage; slope in Three Gorges Reservoir

2015-11-16；

2016-02-02

國(guó)家自然科學(xué)基金項(xiàng)目(51279090)

習(xí)念念(1990-)，女，湖北襄陽(yáng)人，碩士研究生，主要從事水工結(jié)構(gòu)水氣二相流耦合問(wèn)題方面的研究，(電話)15090906272(電子信箱)nn_ctgu@163.com。

童富果(1972-)，男，湖北宜昌人，教授，博士生導(dǎo)師，主要從事水工結(jié)構(gòu)和巖土工程領(lǐng)域多場(chǎng)耦合(THMC)問(wèn)題及其數(shù)值方法方面的研究，(電話)13972560499(電子信箱)tfg@ctgu.edu.cn。

10.11988/ckyyb.20150975

2017,34(1):120-123，141

TV139.1

A

1001-5485(2017)01-0120-04