非飽和滲流問題的自適應欠松弛變量變換方法

程勇剛，常曉林，李典慶，陳 曦

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

1 引 言

非飽和土滲流問題一直是巖土工程、水利工程中的一項熱門研究課題。例如，在邊坡穩定問題的分析中，工程界已普遍認識到非飽和土中的負孔隙水壓力，即基質吸力，對于邊坡的穩定有著正面的貢獻。而隨著降雨的入滲，土體中的基質吸力將逐漸消散，導致其抗剪強度的降低，邊坡的穩定性也將降低，甚至引起失穩。因此，準確地模擬和分析非飽和土滲流問題具有重要的實際意義。

然而，在使用數值模擬方法，如有限元法等求解非飽和土滲流問題時，由于土-水特征曲線和滲透率函數的強烈非線性，有限元計算中經常出現迭代不收斂、計算誤差大等問題。為了提高數值計算的效率和精度，國內外的研究者已經進行了許多有意義的工作。例如，Celia等[1]提出采用混合型Richards方程，依據孔壓變化直接計算含水率變化的改進Picard迭代方法來確保質量平衡，從而提高迭代效率。彭華[2]、周桂云[3]等提出容水度的修正公式，提高了收斂速度。吳夢喜[4]采用積分法處理孔隙水壓力對時間的導數項以提高收斂效率。Pan等[5]提出采用變量變換算法，將滲流控制方程中的未知數（通常為壓力水頭）變換為其他變量，降低方程的非線性程度，從而提高非線性迭代的收斂速度和計算精度。利用這些方法，非飽和土滲流問題數值計算的效率和精度都得到了某種程度的提高。然而，對于一些更為復雜、非線性程度更高的工程實際問題，計算效率和精度的問題仍未得到很好的解決。

本文基于變量變換的思想，結合時間步長自適應技術提出了一種求解非飽和滲流問題的新方法——欠松弛RFT變換方法（ATUR1），并通過兩個二維入滲問題的數值算例驗證了 ATUR1方法在提高計算效率和精度上的效果。

2 欠松弛RFT變換方法

描述非飽和滲流問題的Richards方程可用下式表示：

式中：H為總水頭，等于壓力水頭h與位置水頭z的和；K為孔隙水滲透系數，是壓力水頭h的函數；Q為邊界流量；θ為土壤含水率，也是壓力水頭 h的函數，該函數用土-水特征曲線描述；λ為容水度；t為時間。

由于K、θ函數的強烈非線性，非飽和滲流方程中壓力水頭h的解在空間和時間上的分布也具有強烈的非線性。圖1給出了一個一維非飽和入滲問題中壓力水頭h在空間和時間上的分布的例子。可以看出，在滲透前沿附近，隨著水的入滲，土體含水率提高，壓力水頭h在空間和時間上都有一個劇烈的變化。這種強烈的非線性經常造成數值計算中出現迭代不收斂、計算誤差大等問題。針對這個問題，采用RFT變量變化算法[5]，可表示為

式中：p為變換水頭；h為原壓力水頭；β為變換參數，通常為負數。

通過這種合理變換，將壓力水頭h變換為另一個變量p，可以大大降低Richards方程中未知數在空間和時間上的非線性程度，從而改善這種非線性所帶來的計算收斂困難和精度差等問題。圖1給出了RFT變換的效果。可以看出，與原壓力水頭h的分布相比，變換水頭p在空間和時間上的變化都要平緩得多。

圖1 變量變換算法的效果Fig.1 Effects of variable transformation method

通過RFT變量變換，Richards方程變換為

其中：

相應的有限元格式可表示為

式中：｛p｝為節點變換水頭向量；｛p｝,t為節點變換水頭向量時間偏導；[k ]*為變換后的滲透系數矩陣。d為單元厚度。

采用隱式Euler方法處理時間偏導，得到

式中：｛pn｝為時間步 n時的變換水頭；Δt為時間步長。

為了進一步提高非線性迭代計算的效率，采用欠松弛技術來減少迭代過程中的振蕩現象。Tan等[6]指出，目前常用的欠松弛技術有兩種：一種采用每個時間步中點的平均水頭來定義下一迭代步的土體參數，他們定義為UR1。這種欠松弛技術也被GEOSLOPE公司的SEEP/W軟件所采用[7]；另一種欠松弛方法則采用每個迭代步中點的平均水頭來定義下一迭代步的土體參數，他們定義為 UR2。Tan等[6]的研究指出，UR1欠松弛技術可以大大改善每一個時間步中迭代收斂到穩定解的速度，但在時間步長較大的情況下精度比較差；而UR2欠松弛技術可以在較大的時間步長和較為稀疏的網格時獲得更為精確的解，但在每一個時間步的計算中，迭代收斂的速度相對較差。本文提出的自適應欠松弛變換方法基于UR1欠松弛技術，可定義為ATUR1。而作者前期的研究表明，基于UR2欠松弛技術的變換方法效果較差，因此，不再贅述。在本文所提出的欠松弛 RFT變換方法中，欠松弛技術被應用于變換水頭，即

3 時間步長自適應方法

在常見的非飽和滲流數值計算方法中，一般都采用固定空間網格和固定的時間步長。然而，對于復雜的滲流問題，如同時涉及性質迥異的不同土體的情況，這種固定時間步長的方法可能效率不高。

本文采用一種基于后驗誤差估計的時間步長自適應方法[8]。該方法基于隱式Euler方法，計算每個時間步的時間誤差，并根據此誤差值調整下一個時間步的步長。研究表明，該方法可以有效地控制整個計算過程的誤差。

3.1 誤差估算

在一階精度的隱式Euler方法中，有

而根據二階精度的Thomas-Gladwell算法[9]，有

因此，隱式 Euler方法的誤差可由式（16）和式（18）的差來進行估算，表示為

值得注意的是，當應用于本文提出的欠松弛RFT變換方法時，誤差估算仍舊基于壓力水頭值。

3.2 時間步長自適應

在估算完每個時間步的誤差之后，可根據此誤差值評估當前時間步長是否合適，并調整下一個時間步的步長，即如果誤差滿足如下所示的條件，則接受當前時間步長為

式中：τR和τA分別為相對誤差和絕對誤差允許值；i為節點號。具有最大誤差的節點編號記為iCrit。

如果當前時間步長滿足要求，則繼續進行下一個時間步的計算，新的時間步長調整為

如果當前時間步長不滿足誤差要求，則重新計算當前步，時間步長調整為

rmax、rmin、s、EPS均為控制參數，目的是控制時間步長的變化范圍，從而提高算法的穩定性。根據研究，本文采用如下取值：rmax=4.0，rmin=0.1，s=0.9， EPS = 10-10。前期研究表明，控制參數在上述推薦數值附近的一般變動對自適應方法的效率和穩定性并無本質的影響。

4 數值算例

4.1 Kirkland二維入滲問題

首先研究Kirkland二維入滲問題[10]。由于該問題中土壤的初始負孔隙水壓力極高，而且涉及性質迥異的兩種土，使得數值計算難度較大。圖2給出了該二維入滲問題的示意圖。計算區域由砂土和黏土間隔填充。土體的土-水特征曲線采用 van Genuchten模型[11]，非飽和滲透系數函數采用Mualem模型[12]，具體參數見表1。除表面入滲處入滲量q = 5 cm/d 外其他邊界均為不透水邊界。土體初始負孔隙水壓力水頭為-500 m。有限元網格采用4節點四邊形單元，單元大小取為0.05 m×0.05 m。

圖2 Kirkland入滲 (單位: m)Fig.2 Kirkland's infiltration (unit: m)

表1 Kirkland入滲土體參數Table 1 Soil properties for Kirkland's infiltration

根據本文所提出的時間步長自適應欠松弛RFT變換方法（ATUR1），編制了計算程序 THFELA，對該算例進行計算。計算采用兩組不同的誤差允許值，分別為：①τR=0.5，τA=5.0 m；②τR=0.1，τA=1.0 m。

在變量變換算法中，需要選擇合適的變換參數β。Pan等[5]在提出RFT變量變換算法時，也對此進行了初步的研究。他們指出，可以根據式（4）計算得到的 K*曲線的形狀來確定變換參數β的大小。圖 3給出了不同β值對應的砂土 K*曲線。可以看出，隨著β的減小， K*曲線將由單調遞減變為先增大后減小的非單調曲線。Pan等[5]建議，計算時應選擇使 K*曲線仍維持單調遞減時β的最小允許值。從圖3中可以看出，此時 β≈-2 m-1。然而，在本文所提出的自適應欠松弛 RFT變換方法ATUR1中，由于欠松弛技術的影響，Pan等的建議值不再適用。作者對大量算例的數值研究表明，ATUR1中，變換參數可取為Pan等建議值的1/2，即 K*曲線仍維持單調遞減時β最小允許值的1/2。此時計算精度和效率最高。而當計算中需要處理兩種以上的土體時，首先應按照上述的單調要求，根據各個土體材料的參數計算β的最小允許值，然后選擇其中的最大值。本算例中，砂土材料取β=- 1.0 m-1，黏土材料取 β=- 1.6 m-1，因此，最終取β= -1 .0 m-1。

圖3 不同變換參數時砂土材料的K＊曲線Fig.3 K＊ curves for sand with different β values

圖4 Kirkland入滲水頭等值線圖(ATUR1方法)Fig.4 Pressure head contours of Kirkland's infiltration (ATUR1 method)

圖4給出了入滲12.5 d之后的模擬結果。與采用極小的時間步長和極密的有限元網格的計算結果（可視為準確解）相比，ATUR1的結果完全一致。與文獻[10]中的結果（圖5）相比，ATUR1的結果也基本一致，說明了 ATUR1方法的可靠性和準確性。

圖5 Kirkland入滲的水頭等值線圖(文獻[10])Fig.5 Pressure head contours of Kirkland's Infiltration (from Ref. [10])

4.2 Forsyth二維入滲問題

Forsyth二維入滲問題也是文獻中常用的例子之一[13]。由于該問題涉及性質完全不同的4種土體介質，且空間分布較為復雜，使其數值模擬難度更大。圖6給出了Forsyth二維入滲問題的示意圖。計算區域可分為4個部分，分別由4種不同的土填充。與上例一樣，土體的土-水特征曲線采用 van Genuchten模型[11]，非飽和滲透系數函數采用Mualem模型[12]，具體參數見表2。土體初始負孔隙水壓力水頭為-7.34 m。除表面入滲處入滲量 q =2 cm/d外，其他邊界均為不透水邊界。有限元網格采用4節點四邊形單元，共分為1 780個單元和1 890個節點。

圖6 Forsyth入滲 (單位: m)Fig.6 Forsyth's infiltration (unit: m)

表2 Forsyth入滲土壤參數Table 2 Soil properties for Forsyth's infiltration

采用本文所提出的時間步長自適應欠松弛RFT變換方法（ATUR1）進行模擬。根據4.1節中所述的原則，變換參數可取為 β=-2 .0 m-1。計算采用3組不同的誤差允許值，分別為：①τR=1.0，τA=1.0 m；②τR=0.5，τA=0.5 m；③τR=0.1，τA=0.1 m。

圖7給出了ATUR1方法計算得到的入滲30 d之后的飽和度等值線圖。可以看出，不同的誤差允許值給出的計算結果基本類似。而隨著誤差允許值的減小，計算結果逐漸收斂。與采用極小的時間步長和極密的有限元網格的計算結果（可視為準確解）相比，在采用第3組誤差允許值時，ATUR1的結果與準確解完全一致，說明了 ATUR1方法的可靠性和準確性。

圖7 Forsyth入滲飽和度等值線圖（ATUR1方法）Fig.7 Saturation contours of Forsyth's infiltration (ATUR1 method)

圖8給出了文獻[13]中對此問題的解。可以看出，與ATUR1的結果相比，文獻[13]的結果在左下角有著很大的區別。ATUR1給出的滲透前沿更為陡峭，而文獻[13]的結果彌散度更大。仔細比較發現，文獻[13]的計算中采用的總的時間步數極少（整個計算只用了29次迭代），顯然，時間步長會很大，而這樣大的時間步長明顯不足以得到一個準確的解。文獻[14]也發現了文獻[13]計算中的不足。圖9給出了文獻[14]中對此問題的計算結果。可以看出，如果采用文中所述的TBFN方法，計算所需的時間步數也較少（15個時間步），其計算結果與文獻[13]中的結果類似。而如果采用更為準確的PCOSN方法，計算所需的時間步數將大幅增加（199個時間步），但計算所得到的滲透前沿將更為陡峭，其結果與本文提出的 ATUR1方法的結果基本一致，進一步說明了本文提出的 ATUR1方法的可靠性和準確性。同時也說明，在非飽和滲流問題的數值求解中，必須選擇適當的時間步長或采用時間步長自適應技術，采用不同的方法或不同的時間步長參數進行對比計算，以保證計算結果的準確性。

圖8 Forsyth入滲的飽和度等值線圖(文獻[13])Fig.8 Saturation contours of Forsyth's Infiltration (from Ref. [13])

圖9 Forsyth入滲的飽和度等值線圖(文獻[14])Fig.9 Saturation contours of Forsyth's Infiltration (from Ref. [14])

5 結 論

（1）在非飽和滲流問題的計算中，滲透前沿附近壓力水頭 h在空間和時間上都有一個劇烈的變化，這種強烈的非線性造成了數值計算中迭代不收斂、計算誤差大等問題。

（2）本文提出的欠松弛RFT變換方法（ATUR1）通過變量變換，降低了Richards方程中未知數在空間和時間上的非線性程度，從而改善這種非線性所帶來的計算收斂困難和精度差等問題。欠松弛技術的引入減少了迭代過程中的振蕩現象，進一步提高了非線性迭代計算的效率。時間步長自適應技術則有效地控制整個計算過程的誤差。數值算例說明，ATUR1方法可以有效地提高計算效率和精度，是一種準確有效的計算方法。

（3）在非飽和滲流問題的數值求解中，必須選擇適當的時間步長，并采用不同的方法進行對比計算，以保證計算結果的準確性。

[1]CELIA M A, BOULOUTAS E, ZARBA R L. A general mass-conservative numerical solution for the unsaturated flow equation[J]. Water Resources Research, 1990,26(7): 1483－1496.

[2]彭華，陳勝宏. 飽和-非飽和巖土非穩定滲流有限元分析研究[J]. 水動力學研究與進展(A輯), 2002, 17(2):253－259.PENG Hua, CHEN Sheng-hong. Finite element analyses of unstable seepage in saturated-unsaturated rock[J].Journal of Hydrodynamics (Series A), 2002, 17(2): 253－259.

[3]周桂云. 飽和-非飽和非穩定滲流有限元分析方法的改進[J]. 水利水電科技進展, 2009, 29(1): 5－11.ZHOU Gui-yun. Improvement of finite element analysis of unstable saturated-unsaturated seepage[J]. Advances in Science and Technology of Water Resources, 2009,29(1): 5－11.

[4]吳夢喜. 飽和-非飽和土中滲流Richards方程有限元算法[J]. 水利學報, 2009, 40(10): 1274－1279.WU Meng-xi. Finite element algorithm for Richards’equation for saturated-unsaturated seepage flow[J].Journal of Hydraulic Engineering, 2009, 40(10): 1274－1279.

[5]PAN L, WIERENGA P J. A transformed pressure head-based approach to solve Richards’ equation for variably saturated soils[J]. Water Resources Research,1995, 31(4): 925－931.

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

[7]Geo-Slope International Ltd. SEEP/W engineering book:seepage modeling with SEEP/W[M]. Calgary: Geo-Slope International Ltd., 2004.

[8]KAVETSKI D. Temporal integration in numerical models of unsaturated fluid flow: automatic time-stepping for Richards equation[R]. Newcastle: The University of Newcastle, 2002.

[9]THOMAS R M, GLADWELL I. Variable-order variable-step algorithms for second-order systems (Part I):The methods[J]. International Journal for Numerical Methods in Engineering, 1988, 26(1): 39－53.

[10]KIRKLAND M R, HILLS R G, WIERENGA P J.Algorithms for solving Richards’ equation for variably saturated soils[J]. Water Resources Research, 1992,28(8): 2049－2058.

[11]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.

[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]FORSYTH P A, WU Y S, PRUESS K. Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media[J]. Advances in Water Resources, 1995, 18(1): 25－38.

[14]DIERSCH H J G, PERROCHET P. On the primary variable switching technique for simulating unsaturatedsaturated flows[J]. Advances in Water Resources, 1999,23(3): 271－301.