郝 霜，童富果，劉 剛，薛 松

(三峽大學 水利與環(huán)境學院， 湖北 宜昌 443002)

邊界透氣性對降雨入滲的影響

郝 霜，童富果，劉 剛，薛 松

(三峽大學 水利與環(huán)境學院， 湖北 宜昌 443002)

基于水氣二相流理論，采用有限單元法計算了不同透氣邊界條件下的降雨入滲過程，探討了邊界透氣性對穩(wěn)定入滲強度、孔隙氣壓分布和飽和度分布等的影響。邊界透氣條件是滲流計算中的一個重要控制條件，傳統(tǒng)計算多為定性研究，通過量化透氣條件，分析邊界透氣性對降雨入滲的影響。結(jié)果表明：穩(wěn)定入滲強度與邊界透氣性密切相關(guān)，邊界透氣性越小，穩(wěn)定入滲強度受透氣性影響越敏感，其敏感性隨邊界透氣性增大逐漸降低，邊界透氣性對穩(wěn)定入滲強度的影響呈對數(shù)關(guān)系。此外，邊界透氣性越小，降雨入滲過程中孔隙氣壓力越大，濕潤鋒推移速度越慢。

降雨入滲；透氣條件；水氣兩相流；入滲強度

降雨入滲是一個復雜的非飽和水氣二相流問題，相關(guān)的研究在滑坡治理，城市防洪，保護水土流失等領(lǐng)域應(yīng)用廣泛[1-3]。降雨入滲的本質(zhì)是水氣兩種流體的物質(zhì)傳輸，滲流過程是水流下滲驅(qū)替土體孔隙中空氣，孔隙氣受壓又阻礙水流下滲的過程[4-5]。在物質(zhì)傳輸?shù)臄?shù)值模擬中，物質(zhì)的進出條件是重要的邊界控制條件，控制著物質(zhì)進出的強度、總量和難易程度[6]，因此，在入滲過程的模擬中水、氣進出邊界條件也對滲流過程起著關(guān)鍵的作用。

非飽和土體滲流的邊界問題已有學者做了相關(guān)研究，在國內(nèi)吳良驥等[7]對非飽和滲流問題的研究開展較早，應(yīng)用有限差分法考慮勢邊界、Neuman流量邊界和自由面的混合邊界條件對滲流問題進行了模擬。朱岳明等[8]用有限單元法考慮水相滲流邊界條件對三維非飽和滲流場進行了求解。在國外Gabriela Marinoschi[9]基于Richards方程推導了非飽和土體滲流的非線性模型，應(yīng)用了半滲透的自由邊界條件。大多研究都考慮了滲流邊界透水性對降雨入滲的影響，但降雨入滲是水、氣兩種流體在物質(zhì)傳輸過程中相互耦合作用的過程，因此研究滲流邊界透氣性的影響也同樣重要。

在滲流邊界透氣性的研究中大多為定性分析，李援農(nóng)等[10]利用均質(zhì)土壤積水入滲的試驗，研究了入滲過程中考慮氣阻與孔隙氣自由溢出的邊界透氣條件對降雨入滲的影響；張華等[11]采用室內(nèi)物理試驗的方式研究了積水入滲過程中孔隙中封閉氣泡阻止?jié)B透的機制和規(guī)律；Boris I等[12]基于試驗研究的方法探討了滲流過程封閉氣泡含量對穩(wěn)定入滲率的影響規(guī)律。由于試驗研究無法定量給出邊界透氣程度，也無法定量分析邊界氣相流動對降雨入滲的影響[13]，因此本文基于水、氣二相流的理論，采用有限單元法建立降雨入滲的數(shù)值計算模型，量化邊界透氣性來分析流量邊界上不同程度的透氣條件對降雨入滲的影響。

1 水氣二相流數(shù)值模型

水氣二相流控制微分方程是為描述非飽和多孔介質(zhì)中水、氣運動特征而建立的數(shù)學模型[14-16]。其基本假設(shè)認為土體孔隙是均勻且各向同性的，土、水、氣混合物為連續(xù)介質(zhì)，是滿足連續(xù)介質(zhì)守恒定律，流速緩慢認為其流動狀態(tài)為層流且符合達西定律，水氣間質(zhì)量交換忽略不計。根據(jù)質(zhì)量守恒方程與達西定律可得到孔隙中水相非飽和流動方程：

(1)

對于式(1)中第一項，考慮φw=φSr，可展開為

?t(ρwφw)=φSr?tρw+φρw?tSr+ρwSr?tφ

(2)

考慮水不可壓縮，?tρw=0；不考慮孔隙變形，?tφ=0。則式(2)可變?yōu)?/p>

?t(ρwφw)=φρw?tSr

(3)

對式(1)中第二項，可展開為

(4)

(5)

不考慮內(nèi)源，可得到水流動方程微分形式

(6)

根據(jù)以上方法，不考慮孔隙變形、氣體壓縮性，可以得到氣流動方程

(7)

式(6)和式(7)為非飽和土水氣二相流控制微分方程組。

在求解水氣二相流方程時，選擇孔隙氣壓力和飽和度作為基本變量。孔隙水壓力和孔隙氣壓力關(guān)系可以表示為pw≡pg-pc，pc為基質(zhì)吸力。方程組的空間離散化采用Garlerkin有限單元法，方便引入邊界條件。時間域離散采用一維差分方法，將時間域分為一系列的時段，在每個時段Δt=tk+1-tk內(nèi)假設(shè)變量隨時間線性變化。

2 數(shù)值模型

2.1 幾何模型

本模型為一維均質(zhì)壤土土柱模型，高3m，計算單元20個，計算網(wǎng)格節(jié)點數(shù)42個，網(wǎng)格上密下疏，最小單元厚度為0.1m，計算網(wǎng)格模型及尺寸見圖1。為控制計算土體邊界的滲透性，模型以6號單元底部邊界為底部邊界，1～5號單元為固定飽和度的土體，用以控制底部邊界的透氣程度，因此將模型分為上下兩部分。假設(shè)計算土體為均質(zhì)各向同性土體，不考慮孔隙變形，假定降雨強度遠大于入滲強度，地表層始終處于飽和狀態(tài)，忽略其地表徑流對土體內(nèi)部入滲過程的影響。

2.2 輔助方程及參數(shù)

圖1 計算區(qū)域網(wǎng)格模型及尺寸

土水特征曲線采用Van Genuchten模型(VG模型)，為公式(8)，模型曲線見圖2。

圖2 土水特性曲線

(8)

(9)

式中：Se為有效飽和度；Srw為殘余水飽和度;n、m為與粒徑均勻性有關(guān)的參數(shù)，m=1-1/n；p0為與材料特性相關(guān)的參數(shù)。根據(jù)相關(guān)資料，確定壤土計算參數(shù)p0=0.278,m=0.359,Srw=0.1。

液相非飽和滲透系數(shù)與飽和度的關(guān)系采用Mualem模型，如式(10)

(10)

氣相的滲透率函數(shù)曲線采用Brooks-Corey模型，為公式(11)

(11)

式中：Srg為殘余氣飽和度，Srw為殘余水飽和度。

土樣物理模型參數(shù)及VG模型參數(shù)見表1，液相與氣相滲透特性曲線見圖3。

表1 土樣物理參數(shù)及VG模型參數(shù)

圖3 土體非飽和滲透特性曲線

2.3 邊界條件

土柱整體以同種初始飽和度起算，地表為積水入滲已知流量邊界，底部自由流出，地表與底部均為大氣壓邊界，側(cè)邊界為不透水不透氣邊界。因為在自然界中，所有邊界的封氣條件都不是絕對的，其與邊界下土體厚度，飽和度、孔隙率等有關(guān)，所以為研究邊界不同透氣程度對降雨入滲過程的影響，采用邊界下單位作用力的氣流速作為統(tǒng)一衡量邊界透氣性的標準，無量綱化為模型下部單位氣相流速與上部飽和氣相流速的比值，上部飽和氣相流速取0.000168m/s，殘余飽和度為0.1，氣相黏滯系數(shù)為1.84×10-5。該透氣性標準與計算土柱底邊界以下的土柱飽和度正相關(guān)，對應(yīng)關(guān)系見表2。

表2 邊界透氣性對應(yīng)表

3 滲流邊界透氣性對降雨入滲的影響

降雨入滲過程是一個涉及水、氣兩相流動的物理力學過程，影響降雨入滲的主要因素有土體性質(zhì)及顆粒級配，土體初始飽和度，水力梯度，滲流邊界條件和溫度等，而其中滲流邊界上氣相的流動一定程度上影響了水相入滲。本節(jié)首先分析了透氣性與穩(wěn)定水相入滲率的相關(guān)關(guān)系，而后研究了不同邊界透氣性對入滲過程中土體的孔隙氣壓力、相對入滲率和不同降雨歷時下飽和度分布的影響。

3.1 邊界透氣性對相對穩(wěn)定入滲強度的影響

穩(wěn)定入滲強度是計算土體滲流量，衡量土體入滲能力的重要指標，穩(wěn)定入滲強度與土體顆粒結(jié)構(gòu)、初始飽和度和透氣性等因素有關(guān)，其中邊界透氣條件是自然界中普遍影響入滲強度的因素。本節(jié)模擬了土體初始飽和度為0.3時，不同邊界透氣程度下的相對穩(wěn)定入滲強度，研究了邊界透氣性對穩(wěn)定降雨入滲強度的影響規(guī)律(見圖4)。

圖4 相對穩(wěn)定入滲強度與邊界透氣性的關(guān)系曲線

相對穩(wěn)定入滲強度定義為實際水相穩(wěn)定入滲強度與飽和土體水相穩(wěn)定入滲強度的比值，取無量綱化的比值更有相對意義的對比性，水相飽和穩(wěn)定入滲強度取0.01854m/s。由圖4 可知，穩(wěn)定入滲強度隨土體邊界透氣性的增強而增加，穩(wěn)定入滲強度在土體透氣性較小時敏感性較大，隨邊界透氣性增大敏感性逐漸降低，由此邊界透氣性在較小時對穩(wěn)定入滲強度的影響較大。以對數(shù)函數(shù)擬合邊界透氣性與穩(wěn)定入滲強度的關(guān)系，為y=0.0987lnx+1.7329，擬合確定系數(shù)R2為0.9938，由此可以看出對數(shù)函數(shù)的擬合效果良好。降雨入滲的過程就是土體中水、氣兩種流體在孔隙中相互交替的過程，氣相的流動勢必會影響水相的入滲，而在邊界透氣性較差，氣相密封于土孔隙之中時，氣相對入滲強度的阻滯作用明顯。

3.2 邊界透氣性對入滲過程的影響

土體入滲強度是降雨入滲過程中重要的土體滲透參數(shù)，根據(jù)達西滲流定律v=ki，土體入滲強度是由土體飽和滲透系數(shù)k與水力梯度i決定的，而飽和滲透系數(shù)只與土體本身有關(guān)，因此入滲強度基本由水力梯度決定。在降雨入滲初期表層土體飽和度較小，基質(zhì)吸力較大，水力梯度因此較大，使得降雨初期的入滲強度較大，但隨降雨歷時增加，土體飽和度增加，基質(zhì)吸力快速下降，水力梯度也隨之快速下降，因此在降雨入滲初期入滲強度下降梯度較大直到相對穩(wěn)定。在土體表層飽和后，孔隙氣被表層飽和土體密封于孔隙之中，隨水分不斷入滲，孔隙氣受入滲水分擠壓，使得孔隙氣壓力快速上升并對水分入滲產(chǎn)生阻滯作用，見圖5。在不同的邊界透氣程度下，水相擠壓孔隙氣相而產(chǎn)生的孔隙氣壓力不同。邊界透氣性越小孔隙氣壓力增大越快且維持在一個相對較大的水平，氣壓力消散也較緩慢。邊界透氣性越小穩(wěn)定孔隙氣壓力越大，相應(yīng)對水相入滲阻力越大，入滲強度越小，見圖6，這說明邊界透氣性影響了阻礙了水分入滲，降低了水相入滲強度。

圖5 不同透氣程度下0.3 m處氣壓

本節(jié)計算了透氣性為0.85%、19.87%和57.51%的三種情況在入滲511 min和1005 min時的土體垂直深度處飽和度分布，土體初始飽和度都為0.3，結(jié)果見圖7(a)、圖7(b)。土體垂直深度飽和度分布圖顯示，隨著降雨歷時的增加，土體濕潤鋒隨之下移，土體飽和度從表層到深層逐漸增加，入滲面積隨降雨歷時不斷增大。

圖6 不同透氣程度下土體表層相對入滲強度

圖7 飽和度分布圖

邊界透氣性對土體飽和度分布的影響在土體淺層相對明顯，邊界透氣性越小濕潤鋒下降越緩慢，且隨降雨入滲歷時越久邊界透氣性的影響越大，這說明入滲強度的降低使邊界透氣性對累積入滲量的影響隨降雨歷時的增加逐漸加大，相應(yīng)的濕潤鋒也下降緩慢。

4 結(jié) 論

基于水氣二相流理論與有限單元法，本文分析了邊界透氣性在降雨入滲過程中對穩(wěn)定入滲強度、土體孔隙氣壓力、濕潤鋒和飽和度垂直分布等的影響，可以得出以下結(jié)論。穩(wěn)定降雨入滲強度與邊界透氣性的關(guān)系呈指數(shù)分布，函數(shù)擬合關(guān)系為y=0.0987lnx+1.7329，并且擬合度良好，隨邊界透氣性的增大穩(wěn)定入滲強度也隨之增大，邊界透氣性較小時對穩(wěn)定入滲強度的影響較大。邊界透氣性較小時孔隙氣壓力明顯較高，入滲強度較低，土體濕潤鋒下降緩慢。由此，基于水氣二相流理論的降雨入滲計算需考慮滲流邊界的透氣性，以期達到更精準的滲流模擬過程。

[1] 胡 冉,陳益峰,周創(chuàng)兵.降雨入滲過程中土質(zhì)邊坡的固-液-氣三相耦合分析[J].中國科學:技術(shù)科學,2011，41(11):1469-1482.

[2] 林興旺,簡文彬,董巖松,等.降雨條件下某排土邊坡的滲流場及穩(wěn)定性時變分析[J].水利與建筑工程學報,2015,13(1):36-41.

[3] 荊周寶,劉保健,解新妍,等.考慮流固耦合的降雨入滲過程對非飽和土邊坡的影響研究[J].水利與建筑工程學報,2015,13(6):165-171.

[4] 孫冬梅,朱岳明,張明進,等.考慮氣相影響的降雨入滲過程分析研究[J].巖土力學,2008,29(9):2307-2313.

[5] 張 杰,韓同春,豆紅強,等.探討變雨強條件下的入滲過程及影響因素[J].巖土力學,2014，35(S1):451-456.

[6] 朱元駿,邵明安.含礫石土壤降雨入滲過程模擬[J].水科學進展,2010,21(6):779-787.

[7] 吳良驥, Bloomaburg G L. 飽和一非飽和區(qū)中滲流問題的數(shù)值模型[J].水利水運科學研究,1985(2):1-12.

[8] 朱岳明,龔道勇.三維飽和非飽和滲流場求解及其逸出面邊界條件處理[J].水科學進展,2003,14(1):67-71.

[9] Marinoschi G. A free boundary problem describing the saturated-unsaturated flow in a porous medium[J]. Abstract and Applied Analysis, 2004(9):729-755.

[10] 李援農(nóng),林性粹.均質(zhì)土壤積水入滲的氣阻變化規(guī)律及其影響[J].水土保持學報,1997，3(3):89-93.

[11] 張 華,吳爭光.封閉氣泡對一維積水入滲影響的試驗研究[J].巖土力學,2009，30(S2):132-137.

[12] Borsi I, Farina A, Primicerio M. A rain water infiltration model with unilateral boundary condition: qualitative analysis and numerical simulations[J]. Mathematical Methods in the Applied Sciences, 2006,29(17):2047-2077.

[13] 劉 剛,童富果,習念念,等.通氣和封氣條件下降雨對粘性土入滲速率、含水率及孔隙壓力的影響試驗[J].水電能源科學,2015，33(12):19-21.

[14] Tong F, Niemi A, Yang Z, et al. A numerical model of tracer transport in a non-isothermal two-phase flow system for CO2geological storage characterization[J]. Transport in Porous Media, 2013,98(1):173-192.

[15] Carrick S, Buchan G, Almond P, et al. Atypical early-time infiltration into a structured soil near field capacity: The dynamic interplay between sorptivity, hydrophobicity, and air encapsulation[J]. Geoderma, 2011,160(3/4):579-589.

[16] 付建新,宋衛(wèi)東,杜建華.考慮二維降雨入滲的非飽和土邊坡瞬態(tài)體積含水率分析[J].工程科學學報,2015，37(4):407-413.

The Impact of Different Gas Boundary Condition on Rainfall Infiltration

HAO Shuang, TONG Fuguo, LIU Gang, XUE Song

(CollegeofHydraulic&EnvironmentalEngineering,ChinaThreeGorgesUniversity,Yichang,Hubei443002,China)

Based on the theory of two-phase (water and gas) flow, this study presents results from the finite-element analysis of soil rainfall infiltration on different gas permeable conditions, and analyzes the influence of gas boundary permeability on the stable infiltration intensity, pore gas pressure and soil saturation distribution. Gas permeable boundary is one of the essential control conditions in the seepage simulation. A quantitative analysis has been proposed to figure out the impacts of gas permeable boundary on the process of rainfall infiltration, while a majority of conventional studies focus on qualitative research. The simulation results indicate that there is a highly relevant relationship between stable infiltration intensity and gas permeability on the boundary. The stable infiltration intensity is sensitive when the gas boundary permeability is at a low stage, and as the gas boundary permeability going up the sensibility of stable infiltration intensity is declined. The equation between stable infiltration intensity and gas boundary permeability is regressed to be logarithmic. Besides, during the process of rainfall infiltration the pressure of pore gas soars up and the movement of wetting front slows down with lower gas boundary permeability.

rainfall infiltration; gas boundary condition; two-phase flow of water and gas; infiltration intensity

10.3969/j.issn.1672-1144.2017.02.010

2016-12-28

2017-01-23

國家自然科學基金項目(51279090， 51679129)

郝 霜(1991—)，女，內(nèi)蒙古呼和浩特人，碩士研究生，研究方向為水工結(jié)構(gòu)。E-mail：shuanghao@ctgu.edu.cn

童富果(1972—)，男，湖北宜昌人，教授，博導，主要從事水工結(jié)構(gòu)和巖土工程領(lǐng)域的多場耦合(THMC)問題及其數(shù)值方法研究工作。E-mail：tfg@ctgu.edu.cn

P641.2

A

1672—1144(2017)02—0055—05