基于混合格式Richards方程的非均質(zhì)土體滲流模擬

何玲麗 ，田東方

(1.水利與環(huán)境國家級實(shí)驗(yàn)教學(xué)示范中心，湖北 宜昌 443002； 2.三峽大學(xué) 水利與環(huán)境學(xué)院，湖北 宜昌 443002)

許多巖土工程問題涉及到非飽和滲流過程，例如邊坡工程中的降雨入滲[1]，環(huán)境工程中地下水污染物運(yùn)移等。非飽和滲流過程可由Richards方程描述[2]，該方程有3種基本格式：體積含水率格式（θ-form）、壓力水頭格式（h-form）和混合格式（mixed form）[3-4]。在數(shù)值模擬求解Richards方程時，體積含水率格式的優(yōu)點(diǎn)是允許使用較大時間步長且具有質(zhì)量守恒特性。但由于體積含水率在材料界面不連續(xù)，一般來說該格式適用于均質(zhì)土體；經(jīng)過Hills等[5-7]的發(fā)展，也可用于模擬非均質(zhì)土體，但只適用于非飽和滲流模擬。壓力水頭格式適用于飽和-非飽和滲流模擬，但在模擬干土入滲時，該格式需要采用非常小的時間步長和細(xì)密網(wǎng)格，才能保證質(zhì)量守恒。為此，Phoon等[8-9]提出使用混合格式，同時在迭代過程中將本次迭代的含水率在上次迭代的含水率處按泰勒級數(shù)展開為壓力水頭形式。

由于混合格式具有質(zhì)量守恒特性，因而也被用于主變量轉(zhuǎn)換法（primary switching technique），目的在于結(jié)合體積含水率格式和壓力水頭格式的優(yōu)點(diǎn)，即模擬時對非飽和區(qū)域采用體積含水率格式，飽和區(qū)域采用壓力水頭格式[10]。隨后，Hao等[11-12]提出對非飽和區(qū)域采用混合格式，而飽和區(qū)域采用壓力水頭格式的方法。He等[13]進(jìn)一步指出，混合格式同樣適用于飽和-非飽和滲流模擬。

綜上所述，基于混合格式的數(shù)值模擬具有質(zhì)量守恒特性且能模擬飽和-非飽和滲流的優(yōu)點(diǎn)。但是采用有限元模擬時，該格式需知節(jié)點(diǎn)處的體積含水率，由于體積含水率在材料分界面不連續(xù)，使得模擬非均質(zhì)土體滲流時存在一定困難。目前，針對混合格式模擬時體積含水率不連續(xù)問題的詳細(xì)探討鮮有發(fā)現(xiàn)。因此，本研究目的在于建立一種基于混合格式Richards方程的非均質(zhì)土體飽和-非飽和滲流有限元模擬方法。

1 基本理論

1.1 非飽和滲流控制方程

以忽略源匯項(xiàng)的二維各向同性問題為例，混合格式Richards方程[1]為：

式中：θ為體積含水率[L-1]；t為時間[T]；H為總水頭[L]，H=y+h，y為位置水頭[L]，h為壓力水頭[L]；K=KrKs為滲透系數(shù)[L/T]，Kr為相對滲透系數(shù)，Ks為飽和滲透系數(shù)[L/T]；x、y為坐標(biāo)[L]，y軸豎直向上為正。L和T分別表示長度、時間量綱，下同。

在求解域 Ω 內(nèi)初始條件為：

邊界條件為，在 Γh上為本質(zhì)邊界條件：

在 Γq上為自然邊界條件：

式中：H0、、ˉ為 已知函數(shù)；nx、ny分別為 Γq外法線單位向量在x、y軸的分量。

1.2 Gardner模型

Richards方程的求解還需補(bǔ)充 θ與h的關(guān)系，即土水特征曲線（SWCC），以及滲透性函數(shù)Kr的表達(dá)式，本文采用Gardner模型[14]：

式中：Se為有效飽和度；θr為殘余體積含水率； θs為飽和體積含水率；α為擬合參數(shù)[L-1]。

2 混合格式Richards方程有限元模型

2.1 滲流控制方程的變分提法

用H+和H-表示在材料分界面 Γm兩側(cè)不同材料區(qū)域內(nèi)的總水頭函數(shù)（圖1），則當(dāng) (x,y)∈Γm時，保證總水頭函數(shù)在材料分界面上連續(xù)的條件可寫為：

圖1 離散后的單元及節(jié)點(diǎn)示意Fig.1 Sketch of elements and nodes

設(shè)非飽和滲流問題的泛函記為 π，并假定其滿足本質(zhì)邊界條件（式（3）），則：

式中：c= ?θ/?h，為容水度[L-1]；λ為拉格朗日乘子，是坐標(biāo)的函數(shù)。

可證明泛函的變分 δπ =0與式（1）等價，且包含自然邊界條件（式（4））和連續(xù)性條件（式（7））。證明如下：假設(shè)求解域 Ω 包含若干種不同材料區(qū)域，其中一個區(qū)域用 Ωk表示，的本質(zhì)邊界；表示 Ωk的自然邊界；表示 Ωk的材料分界面；則的變分為：

式（10）中，q按式（4）計(jì)算，q+和q-按下式計(jì)算：

式中：上標(biāo)“+”和“-”表示相應(yīng)物理量屬于兩種不同材料時在邊界的取值。

若 δπ =0，則根據(jù)式（12）等號右邊的第1項(xiàng)可得Richards控制方程（式（1））；第2項(xiàng)可得控制方程需要的自然邊界條件（式（4））；第3項(xiàng)可得總水頭函數(shù)連續(xù)性條件（式（7））；第4項(xiàng)可得 λ的物理意義為邊界流量。

2.2 求解域離散及近似函數(shù)

本文選用四邊形四節(jié)點(diǎn)單元，滲流和拉格朗日乘子兩套網(wǎng)格的拓?fù)湫畔⒋_定如下：

滲流網(wǎng)格為適應(yīng)體積含水率在材料分界面的不連續(xù)，賦予材料分界面上的滲流節(jié)點(diǎn)兩個編號（指整體編號），如圖1中的點(diǎn)p1，在單元中編號為，而在單元中編號為。用數(shù)組Lφ記錄單元的節(jié)點(diǎn)編號，例如

拉格朗日乘子網(wǎng)格中的節(jié)點(diǎn)坐標(biāo)同滲流節(jié)點(diǎn)（位于材料分界面的節(jié)點(diǎn)），但從1開始獨(dú)立編號（指整體編號），如圖1中點(diǎn)p1和p2，編號分別為、。用Lλ記錄單元節(jié)點(diǎn)的整體編號，如Lλ(kλ,1)=。用數(shù)組L記錄單元節(jié)點(diǎn)對應(yīng)的滲流節(jié)點(diǎn)整體編號，例如

式中：Nr為四邊形四節(jié)點(diǎn)單元形函數(shù)；Hr為單元節(jié)點(diǎn)處總水頭值；r=1,2,3,4。

總水頭函數(shù)在求解域 Ω 的近似為：

式中：wr為一維兩節(jié)點(diǎn)單元形函數(shù)； λr為單元節(jié)點(diǎn)處拉格朗日乘子值；r=1,2。

拉格朗日乘子函數(shù)在求解域 Γm的近似為：

2.3 有限元離散格式

將式（14）和（16）代入式（8），并令 δπ =0，可得式（1）的有限元離散格式，寫成矩陣形式：

式（17） 對時間采用向后差分離散后可得：

當(dāng)前，自然資源確權(quán)登記就是將相對完整的生態(tài)功能區(qū)域作為一個自然資源登記單元，自然資源統(tǒng)一確權(quán)登記將各類自然資源的質(zhì)量、數(shù)量和保護(hù)要求全面摸清，并通過登記的法律手段予以公示明確，落實(shí)到每一個產(chǎn)權(quán)人或者使用權(quán)人，有助于充分掌握自然資源家底，并根據(jù)自然資源容量和承載力進(jìn)行分類開發(fā)和保護(hù)，做到自然資源分類施策。

式中：m為 時間步數(shù)；k為第m+1時步內(nèi)的迭代次數(shù)； Δt為時間步長。

為保證求解過程中質(zhì)量守恒，按Celia等[2]的建議，將節(jié)點(diǎn)處體積含水率 θ按下式展開：

式中：C為對角矩陣，c(i,i)= ?θi/?hi，i為節(jié)點(diǎn)號。

略去式（19）中的高階無窮小，將其代入式（18）得：

式（20）即為混合格式Richards方程的有限元求解模型，可用于模擬非均質(zhì)土體飽和-非飽和滲流。該模型將拉格朗日乘子作為未知量與節(jié)點(diǎn)總水頭一起求解，使計(jì)算量有所增加；線性方程組的系數(shù)矩陣對稱但不具備正定性。根據(jù)該式編制Matlab計(jì)算程序，實(shí)現(xiàn)了非飽和滲流數(shù)值模擬。

3 模型驗(yàn)證

第m+1時步的q按下式確定[8]：

3.1 與Srivastava等[15]的解析解對比

計(jì)算模型為一維兩層土土柱，層厚100 cm；計(jì)算域?yàn)閥∈[-100,100]，單位cm。y=-100 cm處為地下水位，即h= 0；頂部y=100 cm處為入滲邊界，以入滲率為q0的穩(wěn)定狀態(tài)為初始條件，當(dāng)t＞ 0時入滲率變?yōu)閝。土水特征曲線和滲透性函數(shù)按式（5）和（6）確定，參數(shù)見圖2。

圖2 入滲與脫濕時壓力水頭模擬與解析解的對比Fig.2 Comparison of analytical and simulated solutions

模擬時計(jì)算網(wǎng)格 Δy=10 cm，允許誤差10-4cm，時間步長 Δt=0.01 h，時長100 h。不同時刻壓力水頭沿y軸分布的模擬結(jié)果與解析解對比見圖2。

從圖2可見本文模型結(jié)果與解析解基本吻合。不同時刻質(zhì)量守恒比率（rmb）見圖3，可見質(zhì)量守恒比率rmb幾乎為1，這表明本文模型所得結(jié)果滿足質(zhì)量守恒。

圖3 質(zhì)量守恒比率隨時間變化Fig.3 Evolution of mass balance ratio

3.2 與Shahrokhabadi等[4]的數(shù)值解對比

計(jì)算模型如圖4所示，坐標(biāo)單位m。底部對應(yīng)地下水位，即h= 0；頂部左側(cè)為入滲邊界，入滲率q=0.7 m/d，其余為不透水邊界。初始條件為H0=0。土水特征曲線和滲透性函數(shù)按式（5）和（6）確定，計(jì)算區(qū)域坐標(biāo)點(diǎn)（坐標(biāo)單位m）和計(jì)算參數(shù)見圖4。模擬時計(jì)算網(wǎng)格 Δx=0.25 m， Δy=0.15 m，允許誤差10-4m，時間步長 Δt=0.001 d，模擬時長 2 d。

圖4 二維層狀土計(jì)算模型（單位：m）Fig.4 Model of layered soil for simulation(unit: m)

選取在0.1 d和2.0 d時刻，對比本文與文獻(xiàn)所得壓力水頭等值線（見圖5），從圖5可知，兩種方法所得結(jié)果差別很小，這表明了本文模型的正確性。

圖5 壓力水頭等值線對比（單位：m）Fig.5 Comparisons of pressure head contours (unit: m)

質(zhì)量守恒比率與時間關(guān)系見圖6，可見質(zhì)量守恒比率rmb幾乎等于1，這表明本文模型所得結(jié)果滿足質(zhì)量守恒。

圖6 質(zhì)量守恒比率隨時間變化Fig.6 Evolution of mass balance ratio

4 結(jié) 語

對材料分界面的節(jié)點(diǎn)賦予兩個編號，將分界面節(jié)點(diǎn)處的體積含水率根據(jù)不同材料按壓力水頭泰勒展開；應(yīng)用拉格朗日乘子法施加總水頭函數(shù)在材料分界面的連續(xù)條件，從而建立一種基于混合格式Richards方程的非均質(zhì)土體飽和-非飽和滲流有限元數(shù)值模型。應(yīng)用所建模型，開展一維層狀土入滲和脫濕過程模擬以及二維層狀土入滲過程模擬，將所得結(jié)果與一維解析解、二維他人文獻(xiàn)結(jié)果進(jìn)行對比，得出本文模型可用于非均質(zhì)土體滲流模擬。質(zhì)量守恒驗(yàn)證表明所建模型滿足質(zhì)量守恒條件。