高等土力学高等土力学 (8).pdf

《高等土力学高等土力学 (8).pdf》由会员分享，可在线阅读，更多相关《高等土力学高等土力学 (8).pdf（7页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第 33 卷增刊 1 岩 土 力 学 Vol.33 Supp.1 2012 年 9 月 Rock and Soil Mechanics Sep.2012 收稿日期：2002-05-28 基金项目：国家自然科学基金资助项目（No.50879039）；973 计划课题资助（No.2010CB732103）；中央高校基本科研业务费专项资金资助（No.2011JBM070）。第一作者简介：陈曦，男，1977 年生，博士，副教授，主要从事岩土工程数值计算方法等方面的研究。E-mail: 文章编号：文章编号：10007598(2012)增 1023707 非饱和渗流非饱和渗流 Richards 方程数值求

2、解的欠松弛方法方程数值求解的欠松弛方法 摘摘 要：要：非饱和土渗流理论是岩土工程问题的基础理论，在土石坝渗流、污染物传输、冻土渗流相变和边坡稳定分析等领域有着广泛的应用。非饱和土渗流 Richards 方程的数值求解过程中，某些参数如水力传导系数计算不当可能引起非线性方法，如 Picard 方法或 Newton 方法的迭代收敛震荡，从而导致非线性迭代方法收敛缓慢和精度降低。为了消除或降低迭代收敛震荡对求解精度和计算性能的影响，目前主要采用欠松弛方法。通过一维入渗算例和二维非均质土坝渗流算例演示已有欠松弛方法的局限性，进而提出新的短项混合欠松弛法，并对其实用性和可靠性进行验证。关关 键键 词：词

3、：非饱和渗流；Richards 方程；有限元；迭代收敛震荡；欠松弛法 中图分类号：中图分类号：TU 443 文献标识码：文献标识码：A Under-relaxation methods for numerical solution of Richards equation of variably saturated flow Abstract:Variably saturated flow theory is the fundamental theory of geotechnical engineering;and it has wide applications to seepage of

4、 earth rockfill dam,contaminant transport,seepage and phase transformation of frozen soil,and slope stability analysis.In numerical solution of Richards equation for variably saturated flow,inappropriate computation of some parameters,particularly for the hydraulic conductivity,may lead to convergen

5、ce oscillation of nonlinear iterative methods such as Picard method or Newton method,and consequently result in slow convergence and inaccurate solution.To eliminate or reduce the influence of iterative convergence oscillation,the under-relaxation method is often employed.By using 1D problem and 2D

6、heterogeneous earth dam,the limitations of available under-relaxation methods are demonstrated.Furthermore,a novel short-term mixed under-relaxation method is proposed with verified practicability and reliability.Key words:variably saturated flow;Richards equation;finite elements;iterative convergen

7、ce oscillation;under-relaxation 1 引 言 非饱和土渗流理论是非饱和土力学理论体系的一个重要组成部分，Richards 方程是非饱和土渗流理论的基本方程，在土石坝渗流、污染物传输、冻土渗流相变和边坡稳定性分析等领域有着广泛的应用，例如，污染物在饱和与非饱和土中迁移的基本方程可以通过 Richards 方程和对流扩散方程导 出1。对于渗流引起的边坡失稳问题，目前主要采用渗流与变形非耦合的分析方法2，研究结果表 明，非饱和渗流场和侵润面的计算对边坡的稳定性和安全系数等结果具有显著的影响34，开展Richards 方程数值求解方法的研究具有十分重要的意义。针对 Ric

8、hards 方程数值求解过程中的各种问题，研究人员已经开展了相关的研究，如吴梦喜等5研究了非饱和渗流求解过程中的数值弥散现象（即土体饱和度比较低处孔隙负压计算的结果分布规律 238 岩 土 力 学 2012 年 较乱，出现与实际物理表现不符的现象），为消除这一现象提出了变坐标的特征有限元法。Richards 方程空间离散和时间差分后，一般可采用非线性迭代方法来求解每一时间步所对应的方程，常用的方法有 Picard 迭代法和 Newton 迭代法，Mehl6通过比较 Picard 迭代法和 Newton 迭代法，得出与 Newton迭代法相比，Picard 的结论通常更加简单有效。非饱和土的水力

9、传导系数（即渗透系数）具有较强的非线性特征，是导致常规 Richards 有限元方程的收敛性较差的主要原因之一，一些研究人员尝试采用变量变换方法来改善 Richards 有限元方程的求解性能。Williams 等7认为，非饱和渗流变量（如压力水头）在较小空间和较短时间范围内的快速变化，是引起 Richards 有限元方程非线性较强和求解困难的主要原因，基于常规的方程格式和求解方法未必有效，建议针对变量进行变换来改善 Richards 有限元方程的非线性和收敛性，并研究不同的变换方 法。Miller 等8提出基于误差控制的自适应策略，将其同时应用于 Richards 方程的空间离散和时间差分，提

10、出空间和时间自适应求解方案。在 Richards 方程的数值求解过程中，某些参数尤其是水力传导系数的计算，需要采用欠松弛（under-relaxation）法,不同的欠松弛法对非饱和 渗流的数值求解精度和计算效率具有显著的影 响910。变量变换技术的性能也是依赖于欠松弛方法11，从基本变量进行欠松弛方法的研究尤为重要。本文通过数值实例指出了现有欠松弛法的局限 性，提出一种新的混合欠松弛方法，并对其实用性和可靠性进行了验证。2 非饱和渗流 Richards 有限元方程 经过多年的发展，非饱和土渗流 Richards 方程已衍生出 3 种基本格式，即压力水头格式（h-form）、混合格式（mixe

11、d form）和体积含水率格式（-form）。尽管混合格式的 Richards 方程被证实具有严格质量守恒的特性，压力水头格式的 Richards 方程在实际工程中的应用仍十分广泛。混合格式的 Richards 方程表达式为()(),hzst=+&Kx （1）式中：为体积含水率，=nSf，n为多孔介质的孔隙率，Sf为流体的饱和度；s为源汇项；K=KsKr()，Kr()为相对渗透张量，为有效饱和度的函数；Ks为饱和渗透系数张量；(x,t)(0,T。非饱和渗流数值分析通常采用Mualem12定义的水力传导系数，即()()()()vvsr21r11 mmlkk kk=（2）式中：ks为饱和水力传导系

12、数；有效饱和度或标准化含水率可由van Genuchten13定义的土-水特征曲线计算，即()()()()vvrvsr1,0 1,0mnhha hhhh=+=（）（）（3）式中：h为压力水头；、s、r分别为体积含水率、饱和含水率和残余含水率；av、nv、mv分别为经验系数，mv通常简化为mv=11/nv。根据比存储量概念（即土-水特征曲线的导数）：()wwddhhm=（4）式中：mw为土-水特征曲线的斜率；w为水的单位重度。利用式（4），将式（1）表示为h-form的Richards方程，即()()(),h hhzst=+&Kx （5）对式（5）应用Galerkin加权残量法，经过空间离散和时

13、间差分后可得h-form的Richards有限元方程：()1111111nnnnnnnntt+=+MChMhb（6）式中：下标n表示时间步指标；M、C 分别为等价质量矩阵和水力传导特征矩阵，由相应单元矩阵Me、Ce组装而成：()TTddeeeeeeeeh=MMN NCCB KB （7）式中：N、B分别为形函数矩阵及其导数矩阵。右端矢量为b=s+g+f，其中：()()TTTTdddeeeheeeeeeeeeszhz=+ssNggB KC zffNK n（8）式中：z为节点高程矢量；n为边界单位法线矢量。增刊 1 陈 曦等：非饱和渗流 Richards 方程数值求解的欠松弛方法 239 式（6）的

14、增量迭代形式又称为修正Picard（这里记为mPicard）迭代，表示如下：()1,1,11,1,11,1,1Pnmnmnmnmnmnm+=+Jhhrhhh （9）式中：()1,1,11,1,1,1,11,PnmnmnnmPnmnmnmnnnmtt+=+=JhMCrJ hMhb（10）当式（9）中第一式中PJ替换为准确雅克比矩阵（JN=J时），mPicard迭代法成为Newton迭代 法。3 非线性迭代收敛震荡与欠松弛方法 在迭代过程中，式（2）中第一式水力传导系数的计算需要输入适当的压力水头：()1,s1,nmrnmkk kh+=（11）根据Tan等9的研究，当1,1,nmnmhh+=时，水

15、力传导系数采用当前时间步、当前迭代步的压力水头进行计算，这种未使用欠松弛的方法记为UR0。然而，直接使用当前迭代步的压力水头来计算水力传导系数可能会引起如图1所示迭代收敛震荡。Phoon等10对这种迭代收敛震荡现象进行了分析，认为上述震荡现象是由水力传导系数相对渗流量计算不协调造成。为了降低迭代收敛震荡对求解精度和计算效率的影响，可采用适当的欠松弛方法。当采用前一时间步结束时的压力水头inh与当前时间步当前迭代步的压力水头的均值1,inmh+来计算水力传导系数时，即()1,1,2iiinmnnmhhh+=+（12）这种方法称为UR1方法。由于UR1方案简单有效，商业软件GeoStudio中的S

16、EEP/W模块采用了这一方法。Tan等9和Phoon等10分别研究了3种不同的欠松弛方法，即上述的UR0、UR1和下面的UR2方法：()1,1,11,2iiinmnmnmhhh+=+（13）这里，UR2方法采用当前时间步最近两次迭代步的压力水头的均值来计算水力传导系数。他们通过一维算例得出“UR2方案在网格较粗的情况下能够给出更加精确的压力水头场，而UR1尽管收敛较快，结果却明显偏离正确解”的结论9-10。采用同一算例对1 m厚、初始条件为干燥的土层进行有表面雨水入渗的模拟，具体参数见章节5.1算例。图1给出10单元网格一维非饱和渗流的计算结果，显 示了UR1只需少量迭代步数便可达到收敛，但其

17、 收敛值（约为-8.0 m）明显偏离了解析解-0.021 6 m。UR2和UR0都能收敛到较为精确的解，但UR2相对UR0需要较少的迭代步数。图1的局部放大图 显示UR2和UR0欠松弛方法作用下迭代法仍存在围绕正确解的迭代收敛震荡现象，正是由于这种迭代收敛震荡导致了Picard或Newton等非线性迭代方法收敛缓慢和数值求解精度降低等问题。图图 1 3 种欠弛松方法作用下种欠弛松方法作用下 Picard 方法的迭代收敛行为方法的迭代收敛行为 Fig.1 Convergence behaviors of Picard iteration 4 一种混合欠松弛方法 理论上，上述的欠松弛方法可以统一在

18、下面的广义欠松弛方法的框架内，其表达形式为 1,011,11,1,iiiiinmnnrnrmnmhhhhh+=+LL （14）式中：r为非线性迭代步数。公式右端系数应满足条件11.0mrr=。式（14）表示水力传导系数计算所采用的压力水头为前一时间步结束时的压力水头和当前时间步所有迭代步已知压力水头的线性组合。很明显，采用广义欠松弛方法来构造近似压力水头主要存在以下两个问题：随着迭代步数的增加，广义欠松弛方法需要存储前面所有迭代步的压力水头；公式右端各项压力水头的系数很难确定。针对上述问 题，可只保留式（14）右端的23项。根据对图1中欠松弛方法的观察，做出如下推测：前一时间050100150

19、-9.0-8.0-7.0-6.0-5.0-4.0-3.0-2.0-1.00.01.00.0z=0.8 m 处的压力水头值/m 迭代步数/步 UR0 UR1 UR2 240 岩 土 力 学 2012 年 步结束时的压力水头hn在欠松弛迭代法中具有阻尼收敛震荡的作用，如UR1；为了获得较为精确的收敛结果，应尽可能保留最近的迭代结果，即 hn+1,m，如UR0和UR2；增加前几步迭代的结果也有阻尼收敛震荡的作用，如UR2中的hn+1,m-1，但阻尼效果较弱。基于上面的理论分析、数值观察和推测，提出一种混合欠松弛方法，即()()1,1,11,11iiiinmnnmnmhhhh+=+（15）式中：0,1

20、和0,1为欠松弛经验参数，由于hn的阻尼作用，也称为阻尼系数。此短项欠松弛法可以看作UR0、UR1和UR2的混合方法，记为mUR(,)，而UR0、UR1和UR2法也可以看作是这种混合欠松弛方法的特例，即当=0，=0，mUR退化为UR0；当=0.5，=0，mUR退化为UR1；当=0，=0.5，mUR成为UR2。根据前面的Richards有限元方程和欠松弛方法编制了非饱和渗流计算程序和可视化界面，通过具体实例指出现有欠松弛方法的局限性，并对所提出的混合欠松弛方法的优越性进行验证。5 数值评价 5.1 一维均质土的非饱和渗流一维均质土的非饱和渗流 算例1.对一维均质土的非饱和非稳态渗流问题进行模拟1

21、0。模型饱和体积含水率s和残余体积含水率r分别为0.363 和0.186，van Genuchten模型中的参数av、nv分别为1 m-1和1.53，饱和渗透系数为ks=1106 m/s。模型中，1 m厚土层划分为10单元有限元网格（即z=0.1 m）。在均匀干燥的初始条件下（设为h(z,0)=-8.0 m），上表面有水逐渐渗入时的边界条件设为h(0,t)=-8.0 m，h(1.0,t)=0.0。计算t=46 800 s的压力水头分布时，随网格和时间步长逐渐加密，压力水头逐渐收敛，如图2所示。采用UR2欠松弛方法，当网格密度达到100单元时，压力水头分布已经足够精确，与1 000个单元网格对应

22、的曲线非常接近。采用1 000个单元时，在0.8 m的位置处，压力水头的计算结果为0.020 8 m，接近解析解为0.021 6 m10。仍采用10单元有限元网格，令混合欠松弛方法中参数。取固定值0.5，观察变化对求解精度和效率的影响。图3为mUR(,0.5)作用下mPicard 迭代法在t=46 800 s时的压力水头分布。这里，mUR(0,0.5)变为UR2欠松弛方法，其结果用十字标记表示。由图可见，随着的增大，压力水头分布逐渐偏离红色虚线表示的密网格解。当=0.5时，压力水头分布与三角形标记表示的UR1方法给出的结果几乎一致，意味着参数影响甚微。图图 2 压力水头分布随网格和时间步长加密

23、的收敛情况压力水头分布随网格和时间步长加密的收敛情况 Fig.2 Convergence of pressure head distribution with the refinement of finite element mesh and time step 图图 3 mUR 中参数中参数 对压力水头分布计算精度的影响对压力水头分布计算精度的影响 Fig.3 Effects of in mUR on solution accuracy 通过上述分析可知，当网格划分较粗时，应足够小（例如(0,0.01）才可确保收敛到较为精确的解。相应地，将mUR(,0.5)作用下mPicard迭代法在z=0

24、.8 m高程处解的迭代收敛过程绘制在图4中，这里选取部分具有代表性的值。有趣的是，随着值的增加，迭代求解的精度逐渐下降，mPicard迭代求解的速度却在增加，且当值足够小时，计算精度下降不是很明显。但=0.01与=0（即UR2方法）相比，mPicard所需的迭代数分别为68步和112步，迭代步减少约为40%；当较-10.0-8.0-6.0-4.0-2.00.00.00.10.20.30.40.50.60.70.80.91.0z=0.1 m,t=46 800 s z=0.05 m,t=11 700 s z=0.02 m,t=1 872 s z=0.01 m,t=468 s z=0.001 m,t

25、=4.68 s 压力水头/m 高程 z/m-10-8-6-4-2020.00.10.20.30.40.50.60.70.80.91.0=0.9 =0.5 =0.3 =0.1 =0.01 =0.005 =0.001 =0.0 压力水头/m 高程 z/m =0.7 增刊 1 陈 曦等：非饱和渗流 Richards 方程数值求解的欠松弛方法 241 大时，迭代收敛尽管很快，结果却偏离正确解较大，这一点与UR1方法相似，印证了前面的推测的合理性，证实了混合欠松弛方法mUR(,)中hn对迭代震荡具有较好的阻尼作用，且此算例中参数益取一个较小的值。图图 4 mUR 中参数中参数 对对 mPicard 迭代

26、性能的影响迭代性能的影响 Fig.4 Effects of in mUR on the performance of mPicard 确定mUR(,)中的参数后（设定为=0.01），还需要进一步研究和评价变化对求解精度和计算效率的影响。根据数值结果，mUR(,)中变化时所得到的压力水头分布与图3中=0.01对应的曲线基本上一致（因此这里没有绘出），表明此算例中求解精度主要由控制。图5显示参数对mPicard的计算性能的影响。由图可见，与UR2的112步（图中红色虚线标定）相比，取0.5以下一定范围的值都可以导致约50%的迭代步数节省，表明混合欠松弛方法在保证精度的同时，可有效提高迭代法的计算效

27、率。图图 5 mUR 中参数中参数 对对 mPicard 迭代性能的影响迭代性能的影响 Fig.5 Effect of in mUR on the performance of mPicard 5.2 二维非均质坝体的非饱和渗流二维非均质坝体的非饱和渗流 算例2.参照于玉贞等的文献4，坝体和2层堤基的非饱和土参数分别为：堤身土，(s,r,av,nv,ks)=(0.441,0.106,1.9 m1，1.31,6.710-7 m/s)；堤基黏土层：(s,r,av,nv,ks)=(0.490,0.123,0.8 m1,1.09,1.810-7 m/s)；堤基砂土层：(s,r,av,nv,ks)=(0

28、.411,0.049,7.5 m1,1.89,210-5 m/s)。采用3组 有限元网格（记为A，B和C），分别含有160、320和2 008个4节点四边形流体单元。对于最密的网格划分C，时间步长设为t=900 s（288个时间步），为了提高计算效率，这里采用稀疏格式的对称连续超松弛（SSOR）预处理的共轭梯度（CG）迭代法并结合Eisenstat技巧14来求解每一非线性迭代步 对称正定线性方程组。坝体一侧初始水位为17 m 标高，另一侧保持恒定水位为10 m标高，在此水位差条件下达稳态渗流。以此稳态渗流状态为初始条件，坝体一侧水位从17 m标高处以1 m d1的速度开始下降，经过3 d，水位

29、下降3 m后坝体内非稳态渗流压力水头分布如图6 所示。图图 6 水位下降水位下降 3 m 后的压力水头分布（单位：后的压力水头分布（单位：m）Fig.6 Distribution of pressure head after 3 m water level drawdown(unit:m)图6根据网格C和UR1方法求出的非稳态渗流压力水头场，在网格C条件下，采用UR1、UR2和mUR几种欠松弛方法所得到的压力水头分布基本一致，因此，可将其用于评价粗网格的计算精度。图7为根据网格A和网格B所获得的0压力水头线，时间步长分别设为t=10 800 s（24个时间步）和t=5 400 s（48个时间步

30、）。由图可见，对于2个粗网格，蓝色实线对应UR1方法所获得的计算精 度（其累计迭代数分别为3 612和8 207）要明显好于红色实线表示的UR2方法（其累计迭代数分别为8 579和18 179），而且随着网格加密，UR1的计算结果较快收敛于正确解（这里指虚线对应的加密网格获得的0压力水头线）。需要强调的是，比较上述两个算例可知，欠松弛方法的性能会因问题不同而呈现显著差别，Tan等9和Phoon等10的数值研究并不全面。参照UR1方法（即=0），固定=0，从0.6变化到1.0，计算得到的0压力水头线与图7中UR1给出的0压力水头线基本一致，表明上述参数y/m x/m 02040 6080 100

31、120-10.0-8.0-6.0-4.0-2.00.02.0=0.1 =0.3 =0.5 迭代步数/步 =0.0 =0.001 =0.005 =0.01 z=0.8m 处的压力水头/m 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.05060708090100110120总迭代步数/步 242 岩 土 力 学 2012 年 情况下计算精度变化甚微。图8为相应的累计迭代数变化曲线，括号中的数字为mPicard未收敛的累计时间步数，这里mPicard采用的收敛允许值为110-5，最大允许迭代数为500。如图所示，保持=0，从0.5变化到1.0时，累计迭代数基本呈下降

32、趋势，表明在这种情况下，采用较大的可以显著减少计算量。对于网格A和网格B，mUR(1.0,0)作用下的累计迭代数分别是2 772和5 770，与UR1作用下的累计迭代数（分别为3 612和8 207）相比，计算量分别节省了23%和30%。(a)网格 A (b)网格 B 图图 7 水位下降水位下降 3 m 后的后的 0 压力水头线压力水头线 Fig.7 The 0 pressure head curves after 3 m water level 图图 8 mUR 中参数中参数 对计算量的影响对计算量的影响 Fig.8 Effects of in mUR on computing cost 6

33、 结 论（1）非饱和渗流数值求解时，水力传导系数相对渗流量计算不协调会造成迭代法收敛震荡，欠松弛方法的主要作用是对这种迭代震荡施加一种阻尼。通过数值试验可知，对于迭代法的收敛震荡，前一时间步结束时的压力水头hn的阻尼抑制效果明显。（2）算例演示了现有两种常用欠松弛方法计算性能的显著差异，表明单一欠松弛方法所存在的局限性。在一维粗网格入渗算例中，UR1方法收敛较快，但收敛结果却偏离正确解较大，UR2方法作用下非线性迭代法可收敛到较为正确的解，但收敛速度相对较慢。二维堤坝渗流算例与一维入渗算例结论不同，UR1方法收敛较快且计算精度也明显好于UR2。（3）基于广义欠松弛方法，发展出的混合欠松弛方法m

34、UR可以看作是UR1和UR2的泛化方法。通过mUR(,)中的系数可调整阻尼程度和欠松弛程度。数值试验表明文中所提出的混合欠松弛方法克服了单一欠松弛方法的局限，具有更好的适用性和计算性能。参参 考考 文文 献献 1 李锡夔.饱和非饱和土壤中污染物运移过程的数值模拟J.力学学报,1998,30(3):321331.LI Xi-kui.Numerical modelling of pollutant transport in unsaturated/saturated soilsJ.Acta Mechanica Sinica,1998,30(3):321332.2 陈曦,刘春杰.逐步完善的非饱和土边

35、坡稳定性有限元分析方法J.岩土工程学报,2011,33(增刊 1):380384.CHEN Xi,LIU Chun-jie.Staged development of finite element methods for stability of unsaturated soil slopesJ.Chinese Journal of Geotechnical Engineering,2011,33(Supp.1):380384.3 黄润秋,戚国庆.非饱和渗流基质吸力对边坡稳定性的影响J.工程地质学报,2002,10(4):343348.HUANG Run-qiu,QI Guo-qing.The

36、 effect of unsaturated soil suction on slope stabilityJ.Journal of Engineering Geology,2002,10(4):343348.4 于玉贞,林鸿州，李荣建,等.非稳定渗流条件下非饱和土边坡稳定分析J.岩土力学,2008,29(11):28922898.UR2 UR1 y/m x/m UR2 UR1 y/m x/m 0.5 0.6 0.7 0.8 0.9 1.02 0003 0004 0005 0006 0007 0008 0009 000网络 A 网络 B(16)(12)(10)(10)(11)(10)(7)(6

37、)(5)(5)(5)(5)累计迭代数/步 增刊 1 陈 曦等：非饱和渗流 Richards 方程数值求解的欠松弛方法 243 YU Yu-zhen,LIN Hung-chou,LI Rong-jian,et al.Stability analysis of unsaturated soil slope under transient seepage flow stateJ.Rock and Soil Mechanics,2008,29(11):28932897.5 吴梦喜,高蓮士.饱和-非饱和土体非稳定渗流数值分析J.水利学报,1999,30(12):3842.WU Meng-xi,GAO L

38、ian-shi.Saturated unsaturated unsteady seepage numerical analysisJ.Journal of Hydraulic Engineering,1999,30(12):3842.6 MEHL S.Use of Picard and Newton iteration for solving nonlinear ground water flow equationsJ.Ground Water,2006,44(4):583594.7 WILLIAMS G A,MILLER C T,KELLEY C T.Transformation appro

39、aches for simulating flow in variably saturated porous mediaJ.Water Resources Research,2000,36(4):923934.8 MILLER C T,ABHISHEK C,FARTHING M W.A spatially and temporally adaptive solution of Richards equationJ.Advances in Water Resources,2006,29(4):525545.9 TAN T S,PHOON K K,CHONG P C.Numerical study

40、 of finite element method based solutions for propagation of wetting frontsJ.Journal of Geotechnical and Geoenvironmental Engineering,ASCE,2004,130(3):254263.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

41、J.Geotechnical and Geological Engineering,2007,25(5):525541.11 CHENG Y G,PHOON K K,TAN T S.Unsaturated soil seepage analysis using a rational transformation method with under-relaxationJ.International Journal of Geomechanics,ASCE,2008,8(3):207212.12 MUALEM Y.A new model for predicting the hydraulic

42、conductivity of unsaturated porous mediaJ.Water Resources Research,1976,12(3):513522.13 VAN GENUCHTEN M T.A closed-form equation for predicting the hydraulic conductivity of unsaturated soilsJ.Soil Science Society of America,1980,44(5):892898.14 CHEN X,TOH K C,PHOON K K.A modified SSOR preconditioner for sparse symmetric indefinite linear systems of equationsJ.International Journal for Numerical Methods in Engineering,2006,65(6):785807.