非飽和滲流Richards方程數(shù)值求解的欠松弛方法

陳 曦，于玉貞，程勇剛

（1. 北京交通大學 土木建筑工程學院，北京 100044；2. 清華大學 水沙科學與水利水電工程國家重點實驗室，北京 100084;3. 武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072）

1 引 言

非飽和土滲流理論是非飽和土力學理論體系的一個重要組成部分，Richards方程是非飽和土滲流理論的基本方程，在土石壩滲流、污染物傳輸、凍土滲流相變和邊坡穩(wěn)定性分析等領域有著廣泛的應用，例如，污染物在飽和與非飽和土中遷移的基本方程可以通過 Richards方程和對流擴散方程導出[1]。對于滲流引起的邊坡失穩(wěn)問題，目前主要采用滲流與變形非耦合的分析方法[2]，研究結(jié)果表明，非飽和滲流場和侵潤面的計算對邊坡的穩(wěn)定性和安全系數(shù)等結(jié)果具有顯著的影響[3－4]，開展Richards方程數(shù)值求解方法的研究具有十分重要的意義。

針對 Richards方程數(shù)值求解過程中的各種問題，研究人員已經(jīng)開展了相關的研究，如吳夢喜等[5]研究了非飽和滲流求解過程中的數(shù)值彌散現(xiàn)象（即土體飽和度比較低處孔隙負壓計算的結(jié)果分布規(guī)律較亂，出現(xiàn)與實際物理表現(xiàn)不符的現(xiàn)象），為消除這一現(xiàn)象提出了變坐標的特征有限元法。Richards方程空間離散和時間差分后，一般可采用非線性迭代方法來求解每一時間步所對應的方程，常用的方法有 Picard迭代法和 Newton迭代法，Mehl[6]通過比較Picard迭代法和Newton迭代法，得出與Newton迭代法相比，Picard的結(jié)論通常更加簡單有效。非飽和土的水力傳導系數(shù)（即滲透系數(shù)）具有較強的非線性特征，是導致常規(guī)Richards有限元方程的收斂性較差的主要原因之一，一些研究人員嘗試采用變量變換方法來改善Richards有限元方程的求解性能。Williams等[7]認為，非飽和滲流變量（如壓力水頭）在較小空間和較短時間范圍內(nèi)的快速變化，是引起Richards有限元方程非線性較強和求解困難的主要原因，基于常規(guī)的方程格式和求解方法未必有效，建議針對變量進行變換來改善Richards有限元方程的非線性和收斂性，并研究不同的變換方法。Miller等[8]提出基于誤差控制的自適應策略，將其同時應用于Richards方程的空間離散和時間差分，提出空間和時間自適應求解方案。

在Richards方程的數(shù)值求解過程中，某些參數(shù)尤其是水力傳導系數(shù)的計算，需要采用欠松弛（under-relaxation）法,不同的欠松弛法對非飽和滲流的數(shù)值求解精度和計算效率具有顯著的影響[9－10]。變量變換技術的性能也是依賴于欠松弛方法[11]，從基本變量進行欠松弛方法的研究尤為重要。本文通過數(shù)值實例指出了現(xiàn)有欠松弛法的局限性，提出一種新的混合欠松弛方法，并對其實用性和可靠性進行了驗證。

2 非飽和滲流Richards有限元方程

經(jīng)過多年的發(fā)展，非飽和土滲流Richards方程已衍生出3種基本格式，即壓力水頭格式（h-form）、混合格式（mixed form）和體積含水率格式（θ-form）。盡管混合格式的Richards方程被證實具有嚴格質(zhì)量守恒的特性，壓力水頭格式的Richards方程在實際工程中的應用仍十分廣泛。

混合格式的Richards方程表達式為

式中：θ 為體積含水率，θ = nSf，n為多孔介質(zhì)的孔隙率，Sf為流體的飽和度；s為源匯項；K =KsKr(Θ)，Kr(Θ)為相對滲透張量，為有效飽和度Θ的函數(shù)；Ks為飽和滲透系數(shù)張量；(x,t)∈Ω ×(0,T ]。非飽和滲流數(shù)值分析通常采用 Mualem[12]定義的水力傳導系數(shù)，即

式中：ks為飽和水力傳導系數(shù)；有效飽和度或標準化含水率Θ可由van Genuchten[13]定義的土-水特征曲線計算，即

式中：h為壓力水頭；θ、θs、θr分別為體積含水率、飽和含水率和殘余含水率；av、nv、mv分別為經(jīng)驗系數(shù)，mv通常簡化為mv= 1-1/nv。

根據(jù)比存儲量概念（即土-水特征曲線的導數(shù)）：

式中：mw為土-水特征曲線的斜率；γw為水的單位重度。

利用式（4），將式（1）表示為h-form的Richards方程，即

對式（5）應用Galerkin加權(quán)殘量法，經(jīng)過空間離散和時間差分后可得h-form的Richards有限元方程：

式中：下標n表示時間步指標；M、C分別為等價質(zhì)量矩陣和水力傳導特征矩陣，由相應單元矩陣Me、Ce組裝而成：

式中： N、B分別為形函數(shù)矩陣及其導數(shù)矩陣。右端矢量為b=s+g+f，其中：

式中：z為節(jié)點高程矢量；n為邊界單位法線矢量。

式（6）的增量迭代形式又稱為修正Picard（這里記為mPicard）迭代，表示如下：

當式（9）中第一式中 JP替換為準確雅克比矩陣（JN= J時），mPicard迭代法成為Newton迭代法。

3 非線性迭代收斂震蕩與欠松弛方法

在迭代過程中，式（2）中第一式水力傳導系數(shù)的計算需要輸入適當?shù)膲毫λ^：

根據(jù)Tan等[9]的研究，當=hn+1,m時，水力傳導系數(shù)采用當前時間步、當前迭代步的壓力水頭進行計算，這種未使用欠松弛的方法記為UR0。然而，直接使用當前迭代步的壓力水頭來計算水力傳導系數(shù)可能會引起如圖 1所示迭代收斂震蕩。Phoon等[10]對這種迭代收斂震蕩現(xiàn)象進行了分析，認為上述震蕩現(xiàn)象是由水力傳導系數(shù)相對滲流量計算不協(xié)調(diào)造成。為了降低迭代收斂震蕩對求解精度和計算效率的影響，可采用適當?shù)那匪沙诜椒?。當采用前一時間步結(jié)束時的壓力水頭與當前時間步當前迭代步的壓力水頭的均值來計算水力傳導系數(shù)時，即

這種方法稱為UR1方法。由于UR1方案簡單有效，商業(yè)軟件GeoStudio中的SEEP/W模塊采用了這一方法。Tan等[9]和 Phoon等[10]分別研究了 3種不同的欠松弛方法，即上述的UR0、UR1和下面的UR2方法：

這里，UR2方法采用當前時間步最近兩次迭代步的壓力水頭的均值來計算水力傳導系數(shù)。他們通過一維算例得出“UR2方案在網(wǎng)格較粗的情況下能夠給出更加精確的壓力水頭場，而UR1盡管收斂較快，結(jié)果卻明顯偏離正確解”的結(jié)論[9-10]。采用同一算例對1 m厚、初始條件為干燥的土層進行有表面雨水入滲的模擬，具體參數(shù)見章節(jié) 5.1算例。圖 1給出 10單元網(wǎng)格一維非飽和滲流的計算結(jié)果，顯示了 UR1只需少量迭代步數(shù)便可達到收斂，但其收斂值（約為-8.0 m）明顯偏離了解析解-0.0216 m。UR2和UR0都能收斂到較為精確的解，但UR2相對 UR0需要較少的迭代步數(shù)。圖 1的局部放大圖顯示UR2和UR0欠松弛方法作用下迭代法仍存在圍繞正確解的迭代收斂震蕩現(xiàn)象，正是由于這種迭代收斂震蕩導致了Picard或Newton等非線性迭代方法收斂緩慢和數(shù)值求解精度降低等問題。

圖1 3種欠弛松方法作用下Picard方法的迭代收斂行為Fig.1 Convergence behaviors of Picard iteration

4 一種混合欠松弛方法

理論上，上述的欠松弛方法可以統(tǒng)一在下面的廣義欠松弛方法的框架內(nèi)，其表達形式為

式（14）表示水力傳導系數(shù)計算所采用的壓力水頭為前一時間步結(jié)束時的壓力水頭和當前時間步所有迭代步已知壓力水頭的線性組合。很明顯，采用廣義欠松弛方法來構(gòu)造近似壓力水頭主要存在以下兩個問題：①隨著迭代步數(shù)的增加，廣義欠松弛方法需要存儲前面所有迭代步的壓力水頭；②公式右端各項壓力水頭的系數(shù)很難確定。針對上述問題，可只保留式（14）右端的2～3項。根據(jù)對圖1中欠松弛方法的觀察，做出如下推測：①前一時間步結(jié)束時的壓力水頭hn在欠松弛迭代法中具有阻尼收斂震蕩的作用，如UR1；②為了獲得較為精確的收斂結(jié)果，應盡可能保留最近的迭代結(jié)果，即hn+1,m，如UR0和UR2；③增加前幾步迭代的結(jié)果也有阻尼收斂震蕩的作用，如UR2中的hn+1,m-1，但阻尼效果較弱。

基于上面的理論分析、數(shù)值觀察和推測，提出一種混合欠松弛方法，即

式中：α∈[0,1]和β∈[0,1]為欠松弛經(jīng)驗參數(shù)，由于hn的阻尼作用，α也稱為阻尼系數(shù)。

此短項欠松弛法可以看作 UR0、UR1和 UR2的混合方法，記為 mUR(α,β)，而 UR0、UR1和UR2法也可以看作是這種混合欠松弛方法的特例，即當α= 0，β= 0，mUR退化為UR0；當α= 0.5，β= 0，mUR退化為UR1；當α= 0，β= 0.5，mUR成為UR2。根據(jù)前面的Richards有限元方程和欠松弛方法編制了非飽和滲流計算程序和可視化界面，通過具體實例指出現(xiàn)有欠松弛方法的局限性，并對所提出的混合欠松弛方法的優(yōu)越性進行驗證。

5 數(shù)值評價

5.1 一維均質(zhì)土的非飽和滲流

算例1. 對一維均質(zhì)土的非飽和非穩(wěn)態(tài)滲流問題進行模擬[10]。模型飽和體積含水率θs和殘余體積含水率θr分別為0.363 和0.186，van Genuchten模型中的參數(shù)av、nv分別為1 m-1和1.53，飽和滲透系數(shù)為ks= 1×106m/s。模型中，1 m厚土層劃分為10單元有限元網(wǎng)格（即Δz = 0.1 m）。在均勻干燥的初始條件下（設為h(z,0) = -8.0 m），上表面有水逐漸滲入時的邊界條件設為h(0,t) = -8.0 m，h(1.0,t)=0.0。計算t = 46800 s的壓力水頭分布時，隨網(wǎng)格和時間步長逐漸加密，壓力水頭逐漸收斂，如圖2所示。采用UR2欠松弛方法，當網(wǎng)格密度達到100單元時，壓力水頭分布已經(jīng)足夠精確，與1000個單元網(wǎng)格對應的曲線非常接近。采用1000個單元時，在0.8 m的位置處，壓力水頭的計算結(jié)果為-0.0208 m，接近解析解為-0.0216 m[10]。

仍采用10單元有限元網(wǎng)格，令混合欠松弛方法中參數(shù)。β取固定值 0.5，觀察α變化對求解精度和效率的影響。圖3為mUR(α,0.5)作用下mPicard迭代法在t = 46800 s時的壓力水頭分布。這里，mUR(0,0.5)變?yōu)閁R2欠松弛方法，其結(jié)果用十字標記表示。由圖可見，隨著α的增大，壓力水頭分布逐漸偏離紅色虛線表示的密網(wǎng)格解。當α= 0.5時，壓力水頭分布與三角形標記表示的 UR1方法給出的結(jié)果幾乎一致，意味著β參數(shù)影響甚微。

圖2 壓力水頭分布隨網(wǎng)格和時間步長加密的收斂情況Fig.2 Convergence of pressure head distribution with the refinement of finite element mesh and time step

圖3 mUR中參數(shù)α 對壓力水頭分布計算精度的影響Fig.3 Effects of α in mUR on solution accuracy

通過上述分析可知，當網(wǎng)格劃分較粗時，α應足夠?。ɡ绂痢?0,0.01]）才可確保收斂到較為精確的解。相應地，將 mUR(α,0.5)作用下 mPicard迭代法在z = 0.8 m高程處解的迭代收斂過程繪制在圖4中，這里選取部分具有代表性的α值。有趣的是，隨著α值的增加，迭代求解的精度逐漸下降，mPicard迭代求解的速度卻在增加，且當α值足夠小時，計算精度下降不是很明顯。但α= 0.01與α=0（即UR2方法）相比，mPicard所需的迭代數(shù)分別為68步和112步，迭代步減少約為40%；當α較大時，迭代收斂盡管很快，結(jié)果卻偏離正確解較大，這一點與UR1方法相似，印證了前面的推測①的合理性，證實了混合欠松弛方法 mUR(α,β)中 hn對迭代震蕩具有較好的阻尼作用，且此算例中參數(shù)α益取一個較小的值。

圖4 mUR中參數(shù)α 對mPicard迭代性能的影響Fig.4 Effects of α in mUR on the performance of mPicard

確定 mUR(α,β)中的參數(shù)α后（設定為α=0.01），還需要進一步研究和評價β變化對求解精度和計算效率的影響。根據(jù)數(shù)值結(jié)果，mUR(α,β)中β變化時所得到的壓力水頭分布與圖3中α= 0.01對應的曲線基本上一致（因此這里沒有繪出），表明此算例中求解精度主要由α控制。圖5顯示β參數(shù)對mPicard的計算性能的影響。由圖可見，與UR2的112步（圖中紅色虛線標定）相比，β取0.5以下一定范圍的值都可以導致約 50%的迭代步數(shù)節(jié)省，表明混合欠松弛方法在保證精度的同時，可有效提高迭代法的計算效率。

圖5 mUR中參數(shù) β 對mPicard迭代性能的影響Fig.5 Effect of β in mUR on the performance of mPicard

5.2 二維非均質(zhì)壩體的非飽和滲流

算例2. 參照于玉貞等的文獻[4]，壩體和2層堤基的非飽和土參數(shù)分別為：堤身土，(θs,θr,av,nv,ks) = (0.441,0.106,1.9 m-1，1.31,6.7×10-7m/s)；堤基黏土層：(θs,θr,av,nv,ks) = (0.490,0.123,0.8 m-1,1.09,1.8×10-7m/s)；堤基砂土層：(θs,θr,av,nv,ks) =(0.411,0.049,7.5 m-1,1.89,2×10-5m/s)。采用 3 組有限元網(wǎng)格（記為A，B和C），分別含有160、320和2008個4節(jié)點四邊形流體單元。對于最密的網(wǎng)格劃分C，時間步長設為Δt = 900 s（288個時間步），為了提高計算效率，這里采用稀疏格式的對稱連續(xù)超松弛（SSOR）預處理的共軛梯度（CG）迭代法并結(jié)合 Eisenstat技巧[14]來求解每一非線性迭代步對稱正定線性方程組。壩體一側(cè)初始水位為 17 m標高，另一側(cè)保持恒定水位為10 m標高，在此水位差條件下達穩(wěn)態(tài)滲流。以此穩(wěn)態(tài)滲流狀態(tài)為初始條件，壩體一側(cè)水位從17 m標高處以1 m d-1的速度開始下降，經(jīng)過3 d，水位下降3 m后壩體內(nèi)非穩(wěn)態(tài)滲流壓力水頭分布如圖6所示。

圖6 水位下降3 m后的壓力水頭分布（單位：m）Fig.6 Distribution of pressure head after 3 m water level drawdown (unit: m)

圖6根據(jù)網(wǎng)格C和UR1方法求出的非穩(wěn)態(tài)滲流壓力水頭場，在網(wǎng)格C條件下，采用UR1、UR2和 mUR幾種欠松弛方法所得到的壓力水頭分布基本一致，因此，可將其用于評價粗網(wǎng)格的計算精度。圖7為根據(jù)網(wǎng)格A和網(wǎng)格B所獲得的0壓力水頭線，時間步長分別設為Δt = 10800 s（24個時間步）和Δt = 5400 s（48個時間步）。 由圖可見，對于2個粗網(wǎng)格，藍色實線對應 UR1方法所獲得的計算精度（其累計迭代數(shù)分別為3612和8207）要明顯好于紅色實線表示的UR2方法（其累計迭代數(shù)分別為8579和18179），而且隨著網(wǎng)格加密，UR1的計算結(jié)果較快收斂于正確解（這里指虛線對應的加密網(wǎng)格獲得的0壓力水頭線）。需要強調(diào)的是，比較上述兩個算例可知，欠松弛方法的性能會因問題不同而呈現(xiàn)顯著差別，Tan等[9]和 Phoon等[10]的數(shù)值研究并不全面。

參照UR1方法（即β= 0），固定β= 0，α從0.6變化到1.0，計算得到的0壓力水頭線與圖7中UR1給出的0壓力水頭線基本一致，表明上述參數(shù)情況下計算精度變化甚微。圖8為相應的累計迭代數(shù)變化曲線，括號中的數(shù)字為mPicard未收斂的累計時間步數(shù)，這里 mPicard采用的收斂允許值為1×10-5，最大允許迭代數(shù)為500。如圖所示，保持β=0，α從0.5變化到1.0時，累計迭代數(shù)基本呈下降趨勢，表明在這種情況下，采用較大的α可以顯著減少計算量。對于網(wǎng)格 A和網(wǎng)格 B，mUR(1.0,0)作用下的累計迭代數(shù)分別是2772和5770，與UR1作用下的累計迭代數(shù)（分別為3612和8207）相比，計算量分別節(jié)省了23%和30%。

圖7 水位下降3 m后的0壓力水頭線Fig.7 The 0 pressure head curves after 3 m water level

圖8 mUR中參數(shù) α 對計算量的影響Fig.8 Effects of α in mUR on computing cost

6 結(jié) 論

（1）非飽和滲流數(shù)值求解時，水力傳導系數(shù)相對滲流量計算不協(xié)調(diào)會造成迭代法收斂震蕩，欠松弛方法的主要作用是對這種迭代震蕩施加一種阻尼。通過數(shù)值試驗可知，對于迭代法的收斂震蕩，前一時間步結(jié)束時的壓力水頭hn的阻尼抑制效果明顯。

（2）算例演示了現(xiàn)有兩種常用欠松弛方法計算性能的顯著差異，表明單一欠松弛方法所存在的局限性。在一維粗網(wǎng)格入滲算例中，UR1方法收斂較快，但收斂結(jié)果卻偏離正確解較大，UR2方法作用下非線性迭代法可收斂到較為正確的解，但收斂速度相對較慢。二維堤壩滲流算例與一維入滲算例結(jié)論不同，UR1方法收斂較快且計算精度也明顯好于UR2。

（3）基于廣義欠松弛方法，發(fā)展出的混合欠松弛方法mUR可以看作是UR1和UR2的泛化方法。通過 mUR(α,β)中的系數(shù)可調(diào)整阻尼程度和欠松弛程度。數(shù)值試驗表明文中所提出的混合欠松弛方法克服了單一欠松弛方法的局限，具有更好的適用性和計算性能。

[1]李錫夔. 飽和－非飽和土壤中污染物運移過程的數(shù)值模擬[J]. 力學學報,1998,30(3): 321－331.LI Xi-kui. Numerical modelling of pollutant transport in unsaturated/saturated soils[J]. Acta Mechanica Sinica,1998,30(3): 321－332.

[2]陳曦,劉春杰. 逐步完善的非飽和土邊坡穩(wěn)定性有限元分析方法[J]. 巖土工程學報,2011,33(增刊1): 380－384.CHEN Xi,LIU Chun-jie. Staged development of finite element methods for stability of unsaturated soil slopes[J].Chinese Journal of Geotechnical Engineering,2011,33(Supp. 1): 380－384.

[3]黃潤秋,戚國慶. 非飽和滲流基質(zhì)吸力對邊坡穩(wěn)定性的影響[J]. 工程地質(zhì)學報,2002,10(4): 343－348.HUANG Run-qiu,QI Guo-qing. The effect of unsaturated soil suction on slope stability[J]. Journal of Engineering Geology,2002,10(4): 343－348.

[4]于玉貞,林鴻州，李榮建,等. 非穩(wěn)定滲流條件下非飽和土邊坡穩(wěn)定分析[J]. 巖土力學,2008,29(11): 2892－2898.YU Yu-zhen,LIN Hung-chou,LI Rong-jian,et al.Stability analysis of unsaturated soil slope under transient seepage flow state[J]. Rock and Soil Mechanics,2008,29(11): 2893－2897.

[5]吳夢喜,高蓮士. 飽和-非飽和土體非穩(wěn)定滲流數(shù)值分析[J]. 水利學報,1999,30(12): 38－42.WU Meng-xi,GAO Lian-shi. Saturated unsaturated unsteady seepage numerical analysis[J]. Journal of Hydraulic Engineering,1999,30(12): 38－42.

[6]MEHL S. Use of Picard and Newton iteration for solving nonlinear ground water flow equations[J]. Ground Water,2006,44(4): 583－594.

[7]WILLIAMS G A,MILLER C T,KELLEY C T.Transformation approaches for simulating flow in variably saturated porous media[J]. Water Resources Research,2000,36(4): 923－934.

[8]MILLER C T,ABHISHEK C,FARTHING M W. A spatially and temporally adaptive solution of Richards’equation[J]. Advances in Water Resources,2006,29(4):525－545.

[9]TAN T S,PHOON K K,CHONG P C. Numerical study of finite element method based solutions for propagation of wetting fronts[J]. Journal of Geotechnical andGeoenvironmental Engineering,ASCE,2004,130(3):254－263.

[10]PHOON K K,TAN T S,CHONG P C. Numerical simulation of richards equation in partially saturated porous media: Under-relaxation and mass balance[J].Geotechnical and Geological Engineering,2007,25(5):525－541.

[11]CHENG Y G,PHOON K K,TAN T S. Unsaturated soil seepage analysis using a rational transformation method with under-relaxation[J]. International Journal of Geomechanics,ASCE,2008,8(3): 207－212.

[12]MUALEM Y. A new model for predicting the hydraulic conductivity of unsaturated porous media[J]. Water Resources Research,1976,12(3): 513－522.

[13]VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America,1980,44(5):892－898.

[14]CHEN X,TOH K C,PHOON K K. A modified SSOR preconditioner for sparse symmetric indefinite linear systems of equations[J]. International Journal for Numerical Methods in Engineering,2006,65(6): 785－807.