非饱和渗流问题的自适应欠松弛变量变换方法

程勇刚，常晓林，李典庆，陈 曦

（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.