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

（1.武漢大學(xué) 水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢 430072；2.北京交通大學(xué) 土木建筑工程學(xué)院，北京 100044）

1 引 言

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

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

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

2 欠松弛RFT變換方法

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

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

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

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

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

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

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

其中：

相應(yīng)的有限元格式可表示為

式中：｛p｝為節(jié)點(diǎn)變換水頭向量；｛p｝,t為節(jié)點(diǎn)變換水頭向量時(shí)間偏導(dǎo)；[k ]*為變換后的滲透系數(shù)矩陣。d為單元厚度。

采用隱式Euler方法處理時(shí)間偏導(dǎo)，得到

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

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

3 時(shí)間步長自適應(yīng)方法

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

本文采用一種基于后驗(yàn)誤差估計(jì)的時(shí)間步長自適應(yīng)方法[8]。該方法基于隱式Euler方法，計(jì)算每個(gè)時(shí)間步的時(shí)間誤差，并根據(jù)此誤差值調(diào)整下一個(gè)時(shí)間步的步長。研究表明，該方法可以有效地控制整個(gè)計(jì)算過程的誤差。

3.1 誤差估算

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

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

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

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

3.2 時(shí)間步長自適應(yīng)

在估算完每個(gè)時(shí)間步的誤差之后，可根據(jù)此誤差值評(píng)估當(dāng)前時(shí)間步長是否合適，并調(diào)整下一個(gè)時(shí)間步的步長，即如果誤差滿足如下所示的條件，則接受當(dāng)前時(shí)間步長為

式中：τR和τA分別為相對(duì)誤差和絕對(duì)誤差允許值；i為節(jié)點(diǎn)號(hào)。具有最大誤差的節(jié)點(diǎn)編號(hào)記為iCrit。

如果當(dāng)前時(shí)間步長滿足要求，則繼續(xù)進(jìn)行下一個(gè)時(shí)間步的計(jì)算，新的時(shí)間步長調(diào)整為

如果當(dāng)前時(shí)間步長不滿足誤差要求，則重新計(jì)算當(dāng)前步，時(shí)間步長調(diào)整為

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

4 數(shù)值算例

4.1 Kirkland二維入滲問題

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

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

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

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

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

圖3 不同變換參數(shù)時(shí)砂土材料的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之后的模擬結(jié)果。與采用極小的時(shí)間步長和極密的有限元網(wǎng)格的計(jì)算結(jié)果（可視為準(zhǔn)確解）相比，ATUR1的結(jié)果完全一致。與文獻(xiàn)[10]中的結(jié)果（圖5）相比，ATUR1的結(jié)果也基本一致，說明了 ATUR1方法的可靠性和準(zhǔn)確性。

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

4.2 Forsyth二維入滲問題

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

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

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

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

圖7給出了ATUR1方法計(jì)算得到的入滲30 d之后的飽和度等值線圖?？梢钥闯觯煌恼`差允許值給出的計(jì)算結(jié)果基本類似。而隨著誤差允許值的減小，計(jì)算結(jié)果逐漸收斂。與采用極小的時(shí)間步長和極密的有限元網(wǎng)格的計(jì)算結(jié)果（可視為準(zhǔn)確解）相比，在采用第3組誤差允許值時(shí)，ATUR1的結(jié)果與準(zhǔn)確解完全一致，說明了 ATUR1方法的可靠性和準(zhǔn)確性。

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

圖8給出了文獻(xiàn)[13]中對(duì)此問題的解?？梢钥闯?，與ATUR1的結(jié)果相比，文獻(xiàn)[13]的結(jié)果在左下角有著很大的區(qū)別。ATUR1給出的滲透前沿更為陡峭，而文獻(xiàn)[13]的結(jié)果彌散度更大。仔細(xì)比較發(fā)現(xiàn)，文獻(xiàn)[13]的計(jì)算中采用的總的時(shí)間步數(shù)極少（整個(gè)計(jì)算只用了29次迭代），顯然，時(shí)間步長會(huì)很大，而這樣大的時(shí)間步長明顯不足以得到一個(gè)準(zhǔn)確的解。文獻(xiàn)[14]也發(fā)現(xiàn)了文獻(xiàn)[13]計(jì)算中的不足。圖9給出了文獻(xiàn)[14]中對(duì)此問題的計(jì)算結(jié)果?？梢钥闯?，如果采用文中所述的TBFN方法，計(jì)算所需的時(shí)間步數(shù)也較少（15個(gè)時(shí)間步），其計(jì)算結(jié)果與文獻(xiàn)[13]中的結(jié)果類似。而如果采用更為準(zhǔn)確的PCOSN方法，計(jì)算所需的時(shí)間步數(shù)將大幅增加（199個(gè)時(shí)間步），但計(jì)算所得到的滲透前沿將更為陡峭，其結(jié)果與本文提出的 ATUR1方法的結(jié)果基本一致，進(jìn)一步說明了本文提出的 ATUR1方法的可靠性和準(zhǔn)確性。同時(shí)也說明，在非飽和滲流問題的數(shù)值求解中，必須選擇適當(dāng)?shù)臅r(shí)間步長或采用時(shí)間步長自適應(yīng)技術(shù)，采用不同的方法或不同的時(shí)間步長參數(shù)進(jìn)行對(duì)比計(jì)算，以保證計(jì)算結(jié)果的準(zhǔn)確性。

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

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

5 結(jié) 論

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

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

（3）在非飽和滲流問題的數(shù)值求解中，必須選擇適當(dāng)?shù)臅r(shí)間步長，并采用不同的方法進(jìn)行對(duì)比計(jì)算，以保證計(jì)算結(jié)果的準(zhǔn)確性。

[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]彭華，陳勝宏. 飽和-非飽和巖土非穩(wěn)定滲流有限元分析研究[J]. 水動(dòng)力學(xué)研究與進(jìn)展(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]周桂云. 飽和-非飽和非穩(wěn)定滲流有限元分析方法的改進(jìn)[J]. 水利水電科技進(jìn)展, 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]. 水利學(xué)報(bào), 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.