改进的Picard法在非饱和土渗流中的应用研究

朱帅润; 李绍红; 何博; 吴礼舟

doi:10.11779/CJGE202204014

改进的Picard法在非饱和土渗流中的应用研究 English Version

朱帅润^{1, 2,},
李绍红³,
何博^{1, 2},
吴礼舟^1, ,

1.
重庆交通大学山区桥梁及隧道工程国家重点实验室, 重庆 400074
2.
成都理工大学环境与土木工程学院, 四川成都 610059
3.
西南交通大学地质工程系, 四川成都 610031

基金项目:

国家重点研发计划项目 2018YFC1504702

国家自然科学基金项目 41790432

地质灾害防治与地质环境保护国家重点实验室自主课题 SKLGP2019Z010

详细信息

作者简介:
朱帅润(1992—)，男，博士研究生，主要从事岩土工程数值计算方法等方面的研究工作。E-mail：zhushuairun@163.com

通讯作者:
吴礼舟, E-mail: lzwu@cqjtu.edu.cn

中图分类号: TU411.4
计量
- 文章访问数: 299
- HTML全文浏览量: 90
- PDF下载量: 187
出版历程
- 收稿日期: 2021-06-23
- 网络出版日期: 2022-09-22
- 刊出日期: 2022-03-31

Application of improved Picard method in unsaturated seepage

ZHU Shuai-run^{1, 2,},
LI Shao-hong³,
HE Bo^{1, 2},
WU Li-zhou^1, ,

1.
State Key Laboratory of Mountain Bridge and Tunnel Engineering, Chongqing Jiaotong University, Chongqing 400074, China
2.
College of Environment and Civil Engineering, Chengdu University of Technology, Chengdu 610059, China
3.
Department of Geological Engineering, Southwest Jiaotong University, Chengdu 610031, China

摘要

摘要: Richards方程广泛应用于非饱和渗流的数值模拟以及其他相关领域。在数值求解中，可以采用有限体积法进行数值离散，进而采用Picard方法进行迭代求解。然而，为了获得可靠准确的数值解，通常均匀网格的空间步长是很小的，特别是在一些不利的数值条件下，比如降雨入渗于干燥土壤中，这往往使得迭代过程耗时，甚至不收敛。因此，结合Chebyshev形式的非均匀网格，提出了一种基于非均匀网格的两层网格校正法的改进Picard迭代方法（NTG-PI）。通过3个非饱和渗流算例，并与传统方法和解析解对比，对改进方法的数值精度、收敛率和加速效果进行了验证。结果表明，相对于传统的Picard方法和自适应松弛Picard方法，提出的方法NTG-PI可以在较少的离散节点数下获得较高的数值精度，以及较高的计算效率。该方法可以对非饱和渗流的数值模拟提供一定参考。
- 非饱和渗流 /
- Picard迭代 /
- 数值精度 /
- 非均匀网格 /
- 两层网格校正方法
Abstract: Richards' equation is widely used in the simulation of unsaturated seepage and related fields. In the numerical solution, the finite volume method can be used for numerical discretization, and then the Picard method can be used for iterative solution. However, in order to obtain reliable and accurate numerical solutions, the spatial step size of a uniform grid is usually very small, especially under some unfavorable numerical conditions, such as rainfall infiltration into dry soil, which often makes the iterative process time-consuming or even unable to converge. Therefore, combined with the non-uniform grid in the form of Chebyshev, an improved Picard iteration method (NTG-PI) is proposed based on the non-uniform two-grid correction scheme. Through three examples of unsaturated seepage, the numerical accuracy, convergence rate and acceleration effect of the improved algorithm are validated by comparing the traditional methods and analytical solutions. The results show that compared with the traditional Picard and the adaptive relaxation Picard methods, the proposed method NTG-PI can obtain higher numerical accuracy with a smaller number of discrete nodes, and also has higher computational efficiency. The proposed method can provide a reference for the numerical simulation of unsaturated seepage.
- unsaturated seepage /
- Picard iteration /
- numerical accuracy /
- non-uniform grid /
- two-grid correction scheme

HTML全文

0. 引言

在不同应力、排水条件和边界条件等因素的影响下，土体颗粒、孔隙或组构的排列在不同方向上有明显差异，天然土和重塑土表现出一定程度的各向异性（初始各向异性）。在加载或变形过程中，土体的组构随着塑性变形而发生改变，从而产生新的各向异性（诱导各向异性）。土体的各向异性微观表现为颗粒排列和孔隙结构的不同^[1-2]。宏观表现为土体的力学性质取决于应力张量和应变张量的方向；屈服面的主方向偏离应力主轴方向；屈服面的形状随各向异性的变化而变化^[3]。倾斜屈服面和旋转硬化法则可以统一描述重塑土的初始各向异性和诱导各向异性。

在临界状态理论框架下，采用类似于剑桥模型的推导方式，通过假定塑性剪切应变和塑性体积应变耦合的功耗散方程和相关联流动法则，得到各向异性模型屈服面和塑性势面表达式^[4]，能够较好地描述屈服面主方向对应力主轴的偏离。许多各向异性模型^[5-7]采用该表达式为各向异性模型的屈服面。但该式对屈服面的描述不准确，运用该式计算的不排水三轴试验的应力路径和预测的水平位移也不准确^[3]。这可能是由于各向异性土屈服面表现出不同的形状。在各向同性模型中，许多学者采用了不同形态的屈服面^[8-10]。例如，剑桥模型的子弹型、修正剑桥模型的椭圆形以及Chen等^[10]提出的其他形状。Lagioia等^[11]基于剪胀和应力比之间的数学关系，提出了一个适用于各向同性土通用的屈服和塑性势面表达式，该表达式可以描述多种形状的屈服面。这些模型，未考虑到组构对屈服面形状的影响。因此可以尝试发展一个形式相对灵活、屈服轨迹形状可以变换的各向异性土屈服面表达式，以实现模型预测功能的改进。

土体是由不规则的土颗粒组成的散粒体结构，在应力作用下表现出各向异性的演化特性。在持续应力作用下，土颗粒通过滑动或滚动调节其排列形式，形成新的各向异性结构以适应应力的变化。这也是各向异性模型建模的核心问题。在目前提出的各种旋转硬化法则中，主要区别在于对不同加载路径土体组构达到的平衡状态的定义。Newson等^[12]和Pestana等^[6]假设在平衡状态时，各向异性土的屈服面倾角与施加的恒定应力比加载路径一致。Dafalias等^[13]认为屈服面倾斜程度与恒定应力比加载路径偏离静水压力轴的程度成正比，比例系数为小于1的常数。Wheeler等^[3]、Dafalias等^[14]则发现，在应力作用下土体各向异性达到平衡状态时，采用非线性描述屈服面倾角与施加的恒定应力比关系时，本构模型的预测效果更好。

以上模型在土体达到临界状态失稳破坏时对土体组构的描述也存在很大的争议。理论上，对临界状态的试样进行不同应力路径的卸载再加载试验，可以测得临界状态土体的屈服面和流动法则。但由于对试样临界状态的预估和保持非常困难，难以通过试验手段实现试样临界状态特性的测量。Wheeler等^[3]、Yang等^[15]和Li等^[16]假设在临界状态下土体达到恒定的各向异性组构；Dafalias等^[14]推测临界状态下各向异性或各向同性组构存在的可能性；Zhao等^[17]和Fu等^[18]采用离散元数值方法对临界状态条件下的土体组构进行研究。然而，这些研究并未得出一致的结论。

在描述各向异性土体变形时，运用单一形式的屈服面和塑性势面预测的三轴不排水试验的应力路径不准确^[19]，许多模型通过采用非关联流法则^{[7, 13, 15]}实现了对模型预测结果的改进。但是由于增加了需要标定的塑性势面参数，使模型比较复杂。可以选择合适屈服面形态来避免使用非相关联流动法则，发展简单的各向异性模型。

基于以上论述，在假定相关联流动法则的前提下，本文通过引入参数n，建立了一个可以描述屈服面形态变化的率无关各向异性模型。在该模型中，假定当剪切变形达到临界状态时，土的各向异性变量为0，即临界状态时，屈服面会旋转至与静水压力轴同轴，并保持不变。该模型在n=1且不考虑屈服面旋转时，可以退化为经典的修正剑桥模型。运用显示积分算法实现模型的数值化并通过模型预测结果与恒定应力比加载试验、变换应力路径试验、排水三轴压缩试验、不排水三轴压缩试验和不排水三轴伸长试验的试验结果对比，验证模型的准确性。

1. 各向异性模型

1.1 通用的各向异性土屈服面表达式

（1）恒定应力比固结阶段，各向异性土的压缩特性

恒定应力比固结即在固结过程中等比例增加平均有效应力p′和偏应力q到目标围压，并保持η=q/p′（K= ${\sigma '_3}$ / ${\sigma '_1}$ =(3-η)/(2η+3)）不变，其中， ${\sigma '_1}$ ， ${\sigma '_2}$ 和 ${\sigma '_3}$ 为有效应力分量。在恒定应力比固结路径下，可以产生稳定的初始各向异性组构。试验结果表明，如图 3（c），在比体积-平均有效应力半对数ν–lnp′空间，恒定应力比固结的土体变形行为可以通过一组与正常固结线斜率相同，近似平行的压缩曲线簇表示^[14]，其斜率与应力比大小无关，仅截距受应力比影响。因此可假设各向异性土在恒定应力比固结作用下的比体积（ν=1+e，e为孔隙比）变化量为

${\text{d}}e = {\text{d}}\nu = - \lambda {\text{d}}(\ln p') {, }$

(1)

式中，λ为土的压缩指数，通过压缩曲线的斜率获得。

（2）恒定应力比固结后，保持 $p'$ 不变的排水剪切阶段，各向异性土的体积变形特性

对于正常固结土或弱超固结土，土体只发生剪切体缩。在平均净应力p′相同时，恒定应力比η=α固结结束时，土体的孔隙比e_α最大；各向异性土剪切到临界状态时，保持体应变不变，设临界状态时土体的孔隙比为e_M，其值最小。因此等平均有效应力剪切过程中土体的孔隙比可采用插值形式表示为

$e = {e_\alpha } - ({e_\alpha } - {e_M})\psi {, }$

(2)

式中，ψ为关于η和α的函数，且满足0≤ψ≤1，即当η=α时，ψ=0；当η=M时，ψ=1。在修正剑桥模型中，ψ=ln(1+η²/M²)/ln2；在Ohno等^[20]给出的各向同性土本构模型中，ψ=ln(1+ηⁿ/Mⁿ)/ln2，实质上是采用了对数型插值函数描述剪胀方程，使各向同性模型可以适应多种形态的屈服面，且当n=2时可退化为修正剑桥模型；在Dafalias^[4]给出的各向异性土的屈服面中，ψ=ln[1+(η-α)²/(M²-α²)]/ln[1+(M-α)/(M+α)]。对于各向异性土，可假设

$\psi = \ln \left( {1{\text{ + }}\frac{{{{({{(\eta - \alpha )}^2})}^n}}}{{{{(M{\text{ + }}\alpha )}^n}{{(M - \alpha )}^n}}}} \right)/\ln \left( {1{\text{ + }}\frac{{{{(M - \alpha )}^n}}}{{{{(M{\text{ + }}\alpha )}^n}}}} \right) {, }$

(3)

式中，M为临界状态应力比，n为控制屈服面形状的参数。当n=1时，ψ=ln[1+(η-α)²/( ${M^2}$ -α²)]/ln[1+(M-α)/(M+α)]，即退化为Dafalias^[4]的模型形式；当α=0时，ψ=ln(1+η²ⁿ/M²ⁿ)/ln2，即为Ohno等^[20]的模型形式；当n=1且α=0时，ψ=ln(1+η²/M²)/ln2，即为修正剑桥模型（MCC）的形式。

对方程（2）求导可得

${\text{d}}e = - ({e_\alpha } - {e_M}){\text{d}}\psi {\rm{。}}$

(4)

（3）任意应力作用下，各向异性土的变形特性

各向异性土在任意应力状态下的体积变形可分解为恒定应力比固结和等p′剪切作用下的体积变形，其固结作用使土体密实并保持稳定的各向异性组构；而剪切作用可以改变土体内部的孔隙、颗粒或团粒排列方向，使颗粒之间发生摩擦滑动，最终导致土体失稳破坏。因此任意应力作用下各向异性土的孔隙比变化量为

${\text{d}}e = - \lambda {\text{d(}}\ln p'{\text{)}} - {\text{(}}{e_\alpha } - {e_M}{\text{)d}}\psi {, }$

(5)

积分得

$e - {e_{\alpha 0}} = - \lambda \ln \frac{{p'}}{{{{p'}_{\alpha 0}}}} - ({e_\alpha } - {e_M})\psi {, }$

(6)

式中，e_α0为形成稳定各向异性组构时的初始孔隙比，相应的平均有效应力为 ${p'_{\alpha 0}}$ 。

因此各向异性土的总体积应变为 ${\varepsilon _{\text{v}}} = \frac{{{e_{\alpha 0}} - e}}{{1 + {e_{\alpha 0}}}}$ ，即

${\varepsilon _{\text{v}}} = \frac{\lambda }{{1 + {e_{\alpha 0}}}}\ln \frac{{p'}}{{{{p'}_{\alpha 0}}}} + \frac{{{e_\alpha } - {e_M}}}{{1 + {e_{\alpha 0}}}}\psi {\rm{。}}$

(7)

（4）弹性变形特性

本文假设各向异性组构对土体弹性变形的影响可以忽略不计，土体的弹性行为是各向同性的，弹性体积应变和弹性剪切应变增量与修正剑桥模型（MCC）相同：

${\text{d}}\varepsilon _{\text{v}}^{\text{e}} = \frac{\kappa }{{1 + {e_{\alpha 0}}}}\frac{{{\text{d}}p'}}{{p'}} = \frac{{{\text{d}}p'}}{K} {, }$

(8)

${\text{d}}\varepsilon _{\text{s}}^{\text{e}} = \frac{{{\text{d}}q}}{{3G}} {, }$

(9)

式中，κ为土体在e–lnp′空间回弹线斜率–回弹指数。 $K = \frac{{{\text{1 + }}{e_{\alpha 0}}}}{\kappa }p'$ 为体积应变模量。G=3K(1-2μ)/(2+2μ)为剪切应变模量，μ为泊松比。

（5）各向异性土屈服面

由式（7），（8）得各向异性土的弹塑性体应变为

$\varepsilon _{\text{v}}^{\text{p}} = \frac{{\lambda - \kappa }}{{1 + {e_{\alpha 0}}}}\ln \frac{{p'}}{{{{p'}_{\alpha 0}}}} + \frac{{{e_\alpha } - {e_M}}}{{1 + {e_{\alpha 0}}}}\psi {\rm{。}}$

(10)

因此各向异性土屈服面为

$f(p', q, \varepsilon _{\text{v}}^{\text{p}}, \alpha ) = \frac{{\lambda - \kappa }}{{1 + {e_{\alpha 0}}}}\ln \frac{{p'}}{{{{p'}_{\alpha 0}}}} + \frac{{{e_\alpha } - {e_M}}}{{1 + {e_{\alpha 0}}}}\psi - \varepsilon _{\text{v}}^{\text{p}} {, }$

(11)

当土体到达临界状态 $\left| \eta \right| = M$ 时， $\dot \varepsilon _{\text{v}}^{\text{p}} = 0$ ，此时，

${\left. {\frac{{\partial f}}{{\partial p'}}} \right|_{\left| \eta \right|{\text{ = }}M}} = \frac{{\lambda - \kappa }}{{1 + {e_{\alpha 0}}}}\frac{1}{{p'}} + \frac{{{e_\alpha } - {e_M}}}{{1 + {e_{\alpha 0}}}}{\left. {\frac{{\partial \psi }}{{\partial p'}}} \right|_{\left| \eta \right|{\text{ = }}M}} = 0 {, }$

(12)

因此，

$\frac{{{e_\alpha } - {e_M}}}{{1 + {e_{\alpha 0}}}} = - \left( {\frac{{\lambda - \kappa }}{{1 + {e_{\alpha 0}}}}} \right)/p'{\left. {\frac{{\partial \psi }}{{\partial p'}}} \right|_{\left| \eta \right|{\text{ = }}M}} {\rm{。}}$

(13)

对式（3）求导并带入临界状态条件 $\left| \eta \right| = M$ ：

$\frac{{\partial \psi }}{{\partial p'}} = \frac{{ - M}}{{p'\ln \left( {1{\text{ + }}\frac{{{{(M - \alpha )}^n}}}{{{{(M{\text{ + }}\alpha )}^n}}}} \right)}}\frac{{2n{{(M - \alpha )}^{2n{\text{ - }}1}}}}{{{{({M^2} - {\alpha ^2})}^n} + {{(M - \alpha )}^{2n}}}} {, }$

(14)

令

$\zeta = \frac{{{{({M^2} - {\alpha ^2})}^n} + {{(M - \alpha )}^{2n}}}}{{2nM{{(M - \alpha )}^{2n{\text{ - }}1}}}} = \frac{{(M - \alpha )}}{{2nM}}\left( {1{\text{ + }}\frac{{{{(M + \alpha )}^n}}}{{{{(M - \alpha )}^n}}}} \right) {, }$

(15)

则

$\frac{{{e_\alpha } - {e_M}}}{{1 + {e_{\alpha 0}}}} = \left( {\frac{{\lambda - \kappa }}{{1 + {e_{\alpha 0}}}}} \right)\zeta \ln \left( {1{\text{ + }}\frac{{{{(M - \alpha )}^n}}}{{{{(M{\text{ + }}\alpha )}^n}}}} \right) {, }$

(16)

因此屈服面为

$\begin{array}{l}f(p', q, \varepsilon _{\text{v}}^{\text{p}}, \alpha )\\ =\frac{\lambda -\kappa }{1+{e}_{\alpha 0}}\mathrm{ln}\frac{{p}^{\prime }}{{{p}^{\prime }}_{\alpha 0}}+\left(\frac{\lambda -\kappa }{1+{e}_{\alpha 0}}\right)\zeta \mathrm{ln}\left(1\text{+}\frac{{({(\eta -\alpha )}^{2})}^{n}}{{({M}^{2}-{\alpha }^{2})}^{n}}\right)-{\epsilon }_{\text{v}}^{\text{p}}\text{ }{, } \end{array}$

(17)

${\rm{令}}\ \ \ \ \ \ {p'_{\alpha {\text{m}}}} = {p'_{\alpha 0}}\exp \left( {\frac{{1 + {e_{\alpha 0}}}}{{\lambda - \kappa }}\varepsilon _{\text{v}}^{\text{p}}} \right) ,$

(18)

则屈服函数可变为

$f(p', q, \alpha , {p'_{\alpha {\text{m}}}}) = {\left( {1{\text{ + }}\frac{{{{({{(\eta - \alpha )}^2})}^n}}}{{{{({M^2} - {\alpha ^2})}^n}}}} \right)^\zeta } - \frac{{{{p'}_{\alpha {\text{m}}}}}}{{p'}} = 0 {, }$

(19)

式中，α为具有应力比类型的变量，满足α < M，控制着屈服面倾角的变化；p′_αm为η=α时的p′值，控制着屈服面的大小。临界状态线通过屈服函数的最高点，即当η=M时， $\partial f/\partial p' = 0$ 。

当n=1时，ζ=1，屈服函数还可表示为

$f(p', q, \alpha , {p'_{\alpha {\text{m}}}}) = \left( {1 + \frac{{{{(\eta - \alpha )}^2}}}{{{M^2} - {\alpha ^2}}}} \right) - \frac{{{{p'}_{a{\text{m}}}}}}{{p'}} {, }$

(20)

式（20）与Dafalias^[4]和Wheeler等^[3]提出的各向异性模型屈服面一致。

当α=0时，ζ=1/n， ${p'_{\alpha {\text{m}}}}$ 即为屈服面与p′轴的交点，表示为 ${p'_0}$ ，各向同性土屈服函数为

$f(p', q, \alpha , {p'_{\alpha {\text{m}}}}) = {\left( {1{\text{ + }}\frac{{{\eta ^{2n}}}}{{{M^{2n}}}}} \right)^{\frac{1}{n}}} - \frac{{{{p'}_0}}}{{p'}} = 0 {\rm{。}}$

(21)

因此，新的屈服面的最大特点就是通过参数n的不同值，改变屈服面的形状。该方法的实质为土体在各向同性体积应变和临界状态体积应变之间选择不同的插值方式，屈服面的形状则不同，如图 1所示，当n < 1时，屈服面为子弹型，当n=1时，屈服面为椭圆形，当n > 1时，屈服面呈现泪滴型，而且屈服面的形状也会因各向异性参数α的变化而变化。

图 1 M=1.1，屈服面的形状

Figure 1. Shapes of yield surface varying with parameter n

下载: 全尺寸图片幻灯片

1.2 各向异性土演化法则

（1）各向同性硬化

各向异性模型的内变量由 ${p'_{\alpha {\text{m}}}}$ 和 $\alpha$ 组成，分别包含了模型的各向同性硬化和旋转硬化。各向同性硬化描述了屈服面形状随着土体密度的变化产生的伸缩变化，可通过e–lnp′空间的正常固结线求得，现有的各向异性模型^{[3, 15, 21]}大多采用与MCC类似的各向同性硬化法则：

${\dot p'_0} = \left\langle L \right\rangle \frac{{1 + e}}{{\lambda - \kappa }}\frac{{\partial g}}{{\partial p}}{p'_0} {, }$

(22)

式中，L为塑性乘子， $\left\langle {} \right\rangle$ 为Macauley括号， $\left\langle L \right\rangle$ 可以通过一致性 $\dot f = 0$ 条件获得， $\left\langle L \right\rangle \partial g{\text{/}}\partial p$ 是塑性体积应变率。

但由于旋转硬化法则的引入，各向异性土的先期固结应力 ${p'_{\alpha {\text{m}}}}$ 随着 $\alpha$ 的变化并不总是沿着p′轴的，仅仅采用式（22）作为各向异性土的各向同性硬化准则而不做变化，计算的误差较大。Leroueil等^[22]利用当前孔隙比和当前应力状态定义了土的无结构的状态，修正的剑桥模型各向同性椭圆表示了这种材料状态，Burland^[23]又将这种材料状态描述为材料的本征状态。本文引入与各向异性方程（19）相对应的各向同性椭圆方程（21）作为各向异性土的本征屈服面，则无需引进新的参数即可建立如下关系：

${p'_0} = {p'_{\alpha {\text{m}}}}{{\text{(1 + }}{\alpha ^{2n}}{\text{/}}{M^{2n}}{\text{)}}^{1/n}} {\rm{。}}$

(23)

${p'_0}$ 的硬化直接采用式（22）。通过式（22），（23）把对 ${p'_{\alpha {\text{m}}}}$ 的硬化转化为对各向同性屈服面的硬化，确保屈服面在旋转时没有过度伸长积累，尤其可以避免即将达到临界状态时屈服面回转到与静水压力轴p′一致时的不断膨胀。

（2）旋转硬化

土体的当前组构依赖于应力历史，控制着土体的力学行为。土体所受应力比的变化可引起土体组构的变化，并产生不可逆的土体变形。塑性各向异性变形引起的组构的变化可用α的变化表示。

α增量可表示为^[5]

$\dot \alpha = \left\langle L \right\rangle c\frac{1}{{{p_{{\text{at}}}}}}\frac{{p'}}{{{{p'}_{\alpha {\text{m}}}}}}\left[ {{\alpha _{\text{e}}}(\eta ) - \alpha } \right] {, }$

(24)

式中，p_at为标准大气压。 ${\alpha _{\text{e}}}$ (η)为恒定应力比加载时α的平衡值，当 $\alpha = {\alpha _{\text{e}}}{\text{(}}\eta {\text{)}}$ ， $\dot \alpha = 0$ 。c为常数，控制α的旋转速率。在当前应力比下，旋转速率依赖于当前组构与平衡组构之间的差值。这个差值越大，旋转速率就越快。

α_e(η)的形式可通过恒定应力比加载试验数据和临界状态破坏时的土体组构获得。尽管临界状态时，土体是否有组构存在争议。微观试验^[24]和数值模拟^[16]表明，在接近临界状态时，土颗粒的片状结构平行于剪切面，这是高度各向异性的；而在一维加载阶段，黏土颗粒的方向逐渐垂直于加载方向。因此，在土体从一维固结到剪切至临界状态破坏时，土体组构张量的主方向并不一致，组构张量发生了主轴旋转。在不考虑土体组构张量主轴旋转的模型中，假设在临界状态时，土体各向异性组构的边界值为0似乎更简单。即假定当应力水平等于或大于临界状态应力水平时，土体结构最终不再是各向异性的。这也被Chen等^[25]证明是热力学允许的。

这意味着临界状态时旋转硬化法则不再发生变化，即当η=M时， $\dot \alpha = 0$ 且α_e(M)=0；当η=0时， ${\alpha _{\text{e}}}(0)=0$ ，即各向同性应力作用下，最终土体组构也是各向同性的。因此可假设：

${\alpha _{\text{e}}}(\eta ) = \eta {m_{\text{L}}}{\left( {{\text{exp}}\left\langle {1 - \frac{{\left| \eta \right|}}{M}} \right\rangle - 1} \right)^{{n_{\text{L}}}}} {, }$

(25)

即当0 < η < M时， ${\alpha _{\text{e}}}(\eta ) = \eta {m_{\text{L}}}{\left[ {{\text{exp}}\left( {1 - \frac{\eta }{M}} \right) - 1} \right]^{{n_{\text{L}}}}}$ ，当-M < η < 0时， ${\alpha _{\text{e}}}(\eta ) = \eta {m_{\text{L}}}{\left[ {{\text{exp}}\left( {1{\text{ + }}\frac{\eta }{M}} \right) - 1} \right]^{{n_{\text{L}}}}}$ ，当|η| > M时， ${\alpha _{\text{e}}}$ (η)=0。式中，m_L > 0和n_L > 0为模型常数，且 ${\alpha _{\text{e}}}$ (0)=0。为Otabiemi clay的 ${\alpha _{\text{e}}}$ 值随η/M的变化图，其中倾斜三角形点为Wheeler等^[3]计算所得，因此所给出的曲线形式能较好地拟合 ${\alpha _{\text{e}}}$ 在压缩侧的变化规律。

图 2 恒定应力比加载时α_e的值

Figure 2. Values of α_e at constant stress ratio

下载: 全尺寸图片幻灯片

1.3 加载指数和塑性模量

采用相关联流动法则，通过一致性条件 $\dot f = 0$ 和式（23）~（25）可得到加载指数 $\left\langle L \right\rangle$ 和塑性模量K_p：

$\left\langle L \right\rangle = \frac{{\frac{{\partial f}}{{\partial p'}}{\text{d}}p' + \frac{{\partial f}}{{\partial q}}{\text{d}}q}}{{{K_{\text{p}}}}} {, }$

(26)

${K}_{\text{p}}=-\left(\frac{\partial f}{\partial {{p}^{\prime }}_{\alpha \text{m}}}\frac{1+e}{\lambda -\kappa }\frac{\partial g}{\partial p}{{p}^{\prime }}_{0}/A+\frac{\partial f}{\partial \alpha }\frac{c}{{p}_{\text{at}}}\frac{{p}^{\prime }}{{{p}^{\prime }}_{\alpha \text{m}}}\left[{\alpha }_{\text{e}}(\eta )-\alpha \right]\right){\rm{。}}$

(27)

1.4 模型参数确定

提出的各向异性本构模型需要8个参数：4个修正剑桥模型参数（λ，κ，G或μ，M），1个用于表示屈服面形状变化的增加参数n，另外3个参数（c，m_L，n_L）描述土体各向异性硬化。这些参数都可以通过常规室内试验得到。例如固结试验可以确定参数（λ，κ，G或μ），三轴剪切试验可以确定土样的抗剪强度参数和M（在压缩侧和伸长侧的M分别表示为M_c，M_e）；完全各向同性固结试样屈服面的最佳拟合或者通过各向同性固结试样的不排水抗剪强度可求得参数n；变换应力比的加载–卸载–再加载试验可以确定参数c，m_L，n_L^[13-14]。模型还包括3个状态参数，初始比体积ν₀、初始固结压力p′_αm和初始倾角α₀，是描述试样初始状态的参数。土体屈服面的大小和初始倾角由先期应力和应变历史决定。

2. 模型验证

为了验证各向异性模型的预测能力，运用matlab软件和Sloan等^[26]提出的带有自动误差控制的显式修正欧拉法(modified Euler method with automatic error control)实现所建立的各向异性模型的数值化。对不同应力路径、排水条件的Lower Cromer Till（简化为LCT）^[27]试验和Kaolin黏土^[28]不排水试验的力学响应进行预测。试验模型参数的取值如表 1。首先，通过模拟结果与原始恒定应力比固结试验对比，验证该模型对固结特性的预测能力。其次，通过改变应力比的固结试验，验证所提出的旋转硬化法则。然后，对比不同固结路径下LCT试样的排水和不排水三轴压缩或伸长试验结果和预测结果。最后，为了说明模型对于黏性土的适用性，对Kaolin黏土不排水试验进行补充验证。图中线代表预测结果，点为试验结果。

表 1 Lower Cromer Till（LCT）和Kaolin土的力学参数

Table 1. Mechanical parameters of LCT and Kaolin clay

试样	M_c	M_e	λ	κ	μ	c	m_L	n_L	n
LCT	1.18	0.95	0.066	0.0077	0.258	80	0.5	0.02	1.8
Kaolin	1.05	0.85	0.140	0.0500	0.300	100	1.0	0.15	1.5

下载: 导出CSV

| 显示表格

2.1 恒定应力比K固结试验

为不同应力比固结试样的 ${\alpha _{\text{e}}}$ 值的拟合结果，图 3（b）为不同应力比固结过程中α值的变化，从图中可以看出固结结束，α值均趋向 ${\alpha _{\text{e}}}$ 值。图 3（c），（e）分别为e–lnp′空间和轴向应变-体积应变空间（ ${\varepsilon _{\text{1}}}$ - ${\varepsilon _{\text{v}}}$ ），恒定应力比取不同值（K为0.4，0.5，0.667，0.8，1）固结到应力p′=233.33 kPa时，固结试验结果和所建模型预测结果对比图。在参数c=0，n=1时，所建立的模型退化为经典的修正剑桥模型MCC，图 3（d），（f）同时给出了修正剑桥模型的预测结果。Gens^[27]所给的固结试验结果表示在w–lnp′空间，其中w为含水率，通过e=wG_s/S_r（G_s=2.65为土颗粒比重，S_r=1为饱和度）将试验结果变换为e–lnp′空间。由图可知，MCC模型明显高估了土体的体积应变和竖向应变，且未能准确预测土体的临界状态；所建模型的预测结果与Gens^[27]的试验数据之间一致性较好，尤其是e–lnp′空间，较好地再现了恒定应力比固结时，土体的体积变形特性，压缩空间沿着临界状态线（η=M，K=0.34）固结的预测也达到了理想的效果；在 ${\varepsilon _{\text{1}}}$ - ${\varepsilon _{\text{v}}}$ 空间，由于采用了相关联流动法则，当K值较小时，所提出的模型的预测结果与试验结果存在较小偏差，但对应变变化趋势的预测效果较好。只有在各向同性固结K=1时，由于土体结构不发生旋转硬化， ${\varepsilon _{\text{1}}}$ - ${\varepsilon _{\text{v}}}$ 是线性的；在各向异性固结K为0.4，0.5，0.667，0.8时， ${\varepsilon _{\text{1}}}$ - ${\varepsilon _{\text{v}}}$ 是向上弯曲的曲线，且K值越小，各向异性程度越大，曲线向上弯曲越多，这与试验结果的变化趋势一致。

图 3 K固结试验结果^[27]和预测结果对比图

Figure 3. Comparison between model predictions and test data for K consolidation tests

下载: 全尺寸图片幻灯片

2.2 应力比变化的再固结试验

变化应力比的再固结试验，是验证旋转硬化法则正确性较直观的方法。Gens^[27]对LCT土进行了大量的变化应力比再固结试验，可分为各向同性固结到各向异性固结、各向异性固结到各向同性固结、各向异性固结到各向异性固结三类变化应力路径的加载方式。在本节所有试验路径中，AB段表示第一阶段加载（初始固结），BC段表示卸载，CD段表示第二阶段再加载（再固结，K₂加载），P点表示屈服点。

（1）各向同性固结到各向异性固结试验

试样IN-4和IN-5的应力路径如图 4（a）所示。图 4（b）及其后面的图中相对体积变化是以卸载结束时（C点）试样的状态为参考计算得到的。在加载的第二阶段由于假定弹性卸载不改变土体组构，因此试验能直观地揭示恒定应力比加载路径下土体的各向异性的演化过程。图 4（b）中可以观察到，预测结果和试验结果的一致性较好，在K₂加载屈服后（P′D′）的初始阶段，曲线发生了一定程度的弯曲，说明所提出的旋转硬化法则能够较好地反映各向同性固结到各向异性固结土体组构的变化，且加载到一定程度时，土体的组构稳定，不再发生演化，此时不同K₂值加载的△V-lnp′近似平行，与所示的恒定应力比固结试验的结果一致。描述了大主应变增量和体积应变增量比 ${\text{d}}{\varepsilon _{\text{1}}}/{\text{d}}{\varepsilon _{\text{v}}}$ 的演化过程。虽然试验测得的应变增量比 ${\text{d}}{\varepsilon _{\text{1}}}/{\text{d}}{\varepsilon _{\text{v}}}$ 具有一点的离散性，模型的预测结果较好地捕捉了应变增量比 ${\text{d}}{\varepsilon _{\text{1}}}/{\text{d}}{\varepsilon _{\text{v}}}$ 的变化趋势。与初始各向同性状态下的 ${\text{d}}{\varepsilon _{\text{1}}}/{\text{d}}{\varepsilon _{\text{v}}} = 1/3$ 相比， ${\text{d}}{\varepsilon _{\text{1}}}/{\text{d}}{\varepsilon _{\text{v}}}$ 的值随着再固结应力比的增大而急剧增大，并逐渐减小趋于稳定。应力比变化越大，应变比增量 ${\text{d}}{\varepsilon _{\text{1}}}/{\text{d}}{\varepsilon _{\text{v}}}$ 的变化就越明显，这是由于施加的应力各向异性使土体的组构向相应的组构平衡状态发生较大的调整。所提出的模型可以较好地再现IN-4和IN-5试样的力学性能。

图 4 IN-4和IN-5预测结果和试验结果^[27]对比

Figure 4. Comparison between model predictions and test data for reconsolidation tests IN-4 and IN-5

下载: 全尺寸图片幻灯片

（2）各向异性固结到各向同性固结试验

Gens^[27]进行的试验K-4、AM-3和K-2的应力路径如图 5（a）所示。上文已分析过K固结试验，这里仅对各向同性再固结阶段进行预测。如图 5（b）所示的体积变形为再固结阶段的初始弹性变形（CP）和不可逆塑性（PD）变形。从弹性状态过渡到弹塑性状态对应着图 5（c）中应变增量比的突然下降，即塑性应变增量方向的迅速变化。在再固结路径结束时，压缩线（5（b））、应变增量比（5（d））均接近渐近状态，基本达到稳定的各向同性状态。在体积变化中观察到的一些离散点（5（b））实际上是由于试样K-4，AM-3和K-2在各向同性再固结过程中弹性区域不同所引起的。在各向同性固结应力作用下，各向异性固结试样颗粒重新排列以适应应力变化，形成新的各向同性组构。因此，所提出的模型能较好地再现K-4、AM-3和K-2试验的再固结过程。

图 5 K-4、AM-3和K-2预测结果和试验结果^[27]对比

Figure 5. Comparison between model predictions and test data for reconsolidation tests K-4, AM-3 and K-2

下载: 全尺寸图片幻灯片

（3）各向异性固结到各向异性固结试验

试验AX-3和AN-2的应力路径如图 6（a）所示。图 6（b）为AX-3和AN-2再固结阶段的体积变化图。图 6（c）为AX-3和AN-2在再固结阶段的应变增量比随平均应力p′的变化图。图 6（d）为AX-3和AN-2在再固结阶段的屈服面倾角随体积应变的变化图。从图中可以看出，与前文分析一致，试样一旦开始屈服，应变速率增量就会急剧变化。再固结比K₂值越小，各向异性程度越强，屈服面和塑性势面向目标倾角连续旋转，土体组构逐渐发生变化直至稳定。再固结至较大应力时，试样在第一阶段固结时形成的初始组构完全发生变化，形成新的稳定组构。对比试验结果和预测结果，所提出的模型能够准确地预测LCT在恒定应力比加载应力路径变化时的力学响应。

图 6 AX-3和AN-2预测结果和试验结果^[27]对比

Figure 6. Comparison between model predictions and test data for reconsolidation tests AX-3 and AN-2

下载: 全尺寸图片幻灯片

2.3 排水试验

图 7对比了恒定应力比K固结下形成的各向异性试样的排水三轴压缩试验的试验结果和模型预测结果。图 7（a）为试验的剪切应力路径，在剪切开始时这些试样均处于正常固结状态，平均有效固结应力为p₀=233.3 kPa，K分别为0.4，0.5，0.667，0.8，1。图 7（b），（d），（f）分别从剪切应力应变、体积应变剪切应变和应力比偏应变的角度对比了所建模型的预测结果。图 7（c），（e），（g）为剑桥模型的预测结果。

图 7 各向异性土排水试验结果^[27]与预测结果对比图

Figure 7. Comparison between model predictions and test data for drained triaxial compression tests with different consolidation stress ratios K for Lower Cromer Till

下载: 全尺寸图片幻灯片

通过对比可知，所建模型较剑桥模型能够更好地预测正常固结的各向异性土的排水强度、体积应变和剪切应变。

2.4 不排水试验

图 8为初始状态与排水试验一致的各向异性LCT试样的不排水三轴压缩和伸长试验的结果图。所建模型较MCC能够更合理地对各向异性不排水应力路径（图 8（a））和剪应力发展（图 8（c））进行预测。所建模型预测的伸长试验应力路径，特别是对于K值较低的固结试样，不排水应力路径初始保持恒定的p′值，表现为剪切启动时的瞬时弹性变形，直到应力点到达屈服面，这与试验结果的差距较大，但较MCC误差较小。这种不一致的结果源于旋转硬化弹塑性理论框架，在该理论框架中，旋转角起“双重”作用，并存在过度的弹性变形^[23]。而且在压缩侧，模型未能预测出低K固结试样的应变软化特性（图 8（c））。这些可能与试验测试技术或模型框架有关，由于在不排水试验中，土体颗粒的平移、滚动和翻转等运动所受的影响因素较多，弹性变形较小且屈服面和塑性势面的差异较大，可以尝试采用非关联流动法则和边界面模型改进模型的预测结果。在伸长试验接近临界状态时，所提出的模型能够表现出试验测得的“钩”型应力路径。

图 8 各向异性Lower Cromer Till不排水试验结果^[27]与预测结果对比图

Figure 8. Comparison between model predictions and test data for undrained triaxial compression/extension tests with different consolidation stress ratios K for Lower Cromer Till

下载: 全尺寸图片幻灯片

Stipho^[28]对Kaolin黏土的各向同性和各向异性正常固结和超固结试样进行了一系列不排水三轴试验。由于超固结土具有应变软化特性，需要在此模型的基础上建立新的边界面模型。为了简便起见，本文只考虑正常固结试样。选择K为0.57，0.667，0.8，1等4种固结试样进行不排水三轴压缩试验，应力路径如图 9（a）所示。对比了Kaolin黏土试样的应力路径、偏应力–应变(q– ${\varepsilon _{\text{1}}}$ )响应和孔隙水压力–偏应变（u_a– ${\varepsilon _{\text{d}}}$ ）关系的试验结果与所建模型和MCC的预测结果。与LCT相比，所提出的模型能更好的模拟Kaolin黏土在不排水三轴应力路径下的力学响应。模型较好地预测了Kaolin黏土在低K固结试验中表现出应变硬化特性。该模型在Kaolin黏土模拟中仍存在伸长侧的不足。

图 9 各向异性Kaolin不排水试验结果^[28]与预测结果对比图

Figure 9. Comparison between model predictions and test data for undrained triaxial compression/extension tests with different consolidation stress ratios K for Kaolin

下载: 全尺寸图片幻灯片

3. 结语

通过合理的假设各向异性土的各向同性体积应变和临界状态体积应变之间的插值函数和恒定应力比加载土体形成的稳定各向异性结构表达式，特别是假定临界状态时土体不再具有各向异性结构，使用相关联流动法则，建立了黏性土的各向异性本构模型。该模型的主要特点是，通过适当调整参数n，可以获得广泛的屈服面形状，屈服面可以呈现泪滴形、子弹头形和椭圆形，且该参数可以通过拟合三轴试验的试验结果获得。在各向同性和临界状态下，旋转硬化法则边界值表达式为0。所提出的模型在n=1且不考虑屈服面旋转时，可以退化为经典的修正剑桥模型。

采用显示积分算法实现了所建立的黏性土各向异性模型的预测功能。验证了LCT在恒定应力比加载试验、变换应力路径试验、排水三轴压缩试验、不排水三轴压缩试验及不排水三轴伸长试验和Kaolin黏土在不排水三轴伸长、压缩试验中的响应特性，结果表明所建立的模型对于黏性土具有较强的适用性。

图 1 基于非均匀网格节点均质土坡的有限体积法和一维入渗模型的示意图

Figure 1. Schematic drawing of control volume method based on non-uniform grid nodes and 1D infiltration model for homogeneous soil slope

下载: 全尺寸图片幻灯片

图 2 非均匀两层网格校正方法的示意图

Figure 2. Schematic drawing of non-uniform two-grid correction scheme

下载: 全尺寸图片幻灯片

图 3 不同迭代方法下获得的数值解和精确解的比较

Figure 3. Comparison between numerical solutions obtained using different iteration methods and exact solutions

下载: 全尺寸图片幻灯片

图 4 比较不同迭代方法的收敛率

Figure 4. Comparison of convergence rates for different iteration methods

下载: 全尺寸图片幻灯片

图 5 SWCC的拟合曲线

Figure 5. Fitting curves of SWCC

下载: 全尺寸图片幻灯片

图 6 不同SWCCs下的数值结果

Figure 6. Numerical results for different SWCCs

下载: 全尺寸图片幻灯片

图 7 两层土质边坡压力水头的计算剖面

Figure 7. Computed profiles of pressure head for two-layer soil slopes

下载: 全尺寸图片幻灯片

图 8 两层土质边坡F_S的计算剖面

Figure 8. Computed profiles of F_S for two-layer soil slopes

下载: 全尺寸图片幻灯片

表 1 t = 6 h时不同方法的数值精度对比

Table 1 Comparison of numerical accuracy of different methods at t = 6 h

数值离散条件		RSE		RE/%		MRE/%
N	$\Delta t$	PI	NTG-PI	PI	NTG-PI	PI	NTG-PI
100	0.10	0.135	1.6×10^-2	5.87	1.330	55.03	2.98
	0.05	0.138	9.1×10^-3	6.70	0.570	55.44	3.25
	0.01	0.141	5.9×10^-3	7.35	0.032	55.76	3.46
200	0.10	0.121	1.5×10^-2	4.67	1.460	56.32	1.87
	0.05	0.123	7.7×10^-3	5.52	0.710	56.75	0.94
	0.01	0.126	2.0×10^-3	6.19	0.100	57.09	0.83
400	0.10	0.090	1.5×10^-2	2.57	1.490	53.15	1.86
	0.05	0.093	7.7×10^-3	3.44	0.740	53.67	0.93
	0.01	0.095	1.6×10^-3	4.14	0.130	54.09	0.18

下载: 导出CSV

表 2 在N= 100时比较改进方法相对于N-PI的加速比

Table 2 Comparison of speed-up ratios of improved schemes relative to N-PI when N= 100

条件	加速比
$\Delta t$	${S_{{\text{N - PI/NAR - PI}}}}$	${S_{{\text{N - PI/NTG - PI}}}}$
0.1 h	3.33	46.69
0.05 h	1.72	25.99
0.01 h	1.07	14.12

下载: 导出CSV

表 3 参数拟合值

Table 3 Fitting values of parameters

模型	${\theta _{\text{r}}}$	${\theta _{\text{s}}}$	$\left\| \alpha \right\|$	$\beta$	n	R²
Gardner	0.02	0.43	0.011	—	—	0.945
van Genuchten	0.02	0.43	—	0.014	1.74	0.953

下载: 导出CSV

表 4 比较不同SWCCs下改进方法相对于N-PI的加速比

Table 4 Comparison of speed-up ratios of improved schemes relative to N-PI under different SWCCs

条件	加速比 ${S_{{\text{N - PI/NTG - PI}}}}$
N	Gardner	van Genuchten
120	1.66	1.38
240	2.39	2.30
480	3.35	4.23

下载: 导出CSV

表 5 两层土质边坡的模型参数

Table 5 Model parameters for two-layer soil slopes

土层参数	${K_{\text{s}}}$ /(m·h^-1)	${\theta _{\text{s}}}$	${\theta _{\text{r}}}$	$\beta$	n
土层1	5.46×10^-3	0.4686	0.106	1.04	1.3954
土层2	2.25×10^-1	0.3658	0.0286	2.8	2.239

下载: 导出CSV

参考文献(31)

[1]	李梦姿, 蔡国庆, 李昊, 等. 考虑抗拉强度剪断的非饱和土无限边坡稳定性分析[J]. 岩土工程学报, 2020, 42(4): 705–713. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202004017.htm LI Meng-zi, CAI Guo-qing, LI Hao, et al. Stability of infinite unsaturated soil slopes with tensile strength cut-off[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(4): 705–713. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202004017.htm
[2]	覃小华, 刘东升, 宋强辉, 等. 强降雨条件下考虑饱和渗透系数变异性的基岩型层状边坡可靠度分析[J]. 岩土工程学报, 2017, 39(6): 1065–1073. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201706014.htm QIN Xiao-hua, LIU Dong-sheng, SONG Qiang-hui, et al. Reliability analysis of bedrock laminar slope stability considering variability of saturated hydraulic conductivity of soil under heavy rainfall[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(6): 1065–1073. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201706014.htm
[3]	WU L Z, SELVADURAI A P S, ZHANG L M, et al. Poro-mechanical coupling influences on potential for rainfall-induced shallow landslides in unsaturated soils[J]. Advances in Water Resources, 2016, 98: 114–121. doi: 10.1016/j.advwatres.2016.10.020
[4]	陈永贵, 蔡叶青, 叶为民, 等. 处置库膨润土胶体吸附迁移性及核素共同迁移特性研究进展[J]. 岩土工程学报, 2021, 43(12): 2149–2158. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202112001.htm CHEN Yong-gui, CAI Ye-qing, YE Wei-min, et al. Progresses in researches on adsorption and migration properties of bentonite colloids and their co-migration with nuclide in repository[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(12): 2149–2158. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202112001.htm
[5]	RICHARDS L A. Capillary conduction of liquids through porous mediums[J]. Physics, 1931, 1(5): 318–333. doi: 10.1063/1.1745010
[6]	WU L Z, HUANG J S, FAN W, et al. Hydro-mechanical coupling in unsaturated soils covering a non-deformable structure[J]. Computers and Geotechnics, 2020, 117: 103287. doi: 10.1016/j.compgeo.2019.103287
[7]	SRIVASTAVA R, YEH T C J. Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils[J]. Water Resources Research, 1991, 27(5): 753–762. doi: 10.1029/90WR02772
[8]	TRACY F T. Clean two- and three-dimensional analytical solutions of Richards' equation for testing numerical solvers[J]. Water Resources Research, 2006, 42(8): 1–11.
[9]	PANICONI C, PUTTI M. Physically based modeling in catchment hydrology at 50: survey and outlook[J]. Water Resources Research, 2015, 51(9): 7090–7129. doi: 10.1002/2015WR017780
[10]	ZHA Y Y, YANG J Z, SHI L S, et al. Simulating one-dimensional unsaturated flow in heterogeneous soils with water content-based richards equation[J]. Vadose Zone Journal, 2013, 12(2): 1–13. doi: 10.2136/vzj2012.0142
[11]	AN H, ICHIKAWA Y, TACHIKAWA Y, et al. Three-dimensional finite difference saturated-unsaturated flow modeling with nonorthogonal grids using a coordinate transformation method[J]. Water Resources Research, 2010, 46(11): 1–18.
[12]	ZHANG Z Y, WANG W K, YEH T C J, et al. Finite analytic method based on mixed-form Richards' equation for simulating water flow in vadose zone[J]. Journal of Hydrology, 2016, 537: 146–156. doi: 10.1016/j.jhydrol.2016.03.035
[13]	HERRERA P A, MASSABÓ M, BECKIE R D. A meshless method to simulate solute transport in heterogeneous porous media[J]. Advances in Water Resources, 2009, 32(3): 413–429. doi: 10.1016/j.advwatres.2008.12.005
[14]	WU L Z, ZHU S R, PENG J B. Application of the Chebyshev spectral method to the simulation of groundwater flow and rainfall-induced landslides[J]. Applied Mathematical Modelling, 2020, 80: 408–425. doi: 10.1016/j.apm.2019.11.043
[15]	LOTT P A, WALKER H F, WOODWARD C S, et al. An accelerated Picard method for nonlinear systems related to variably saturated flow[J]. Advances in Water Resources, 2012, 38: 92–101. doi: 10.1016/j.advwatres.2011.12.013
[16]	BRENNER K, CANCÈS C. Improving Newton's method performance by parametrization: the case of the richards equation[J]. SIAM Journal on Numerical Analysis, 2017, 55(4): 1760–1785. doi: 10.1137/16M1083414
[17]	ZENG J C, ZHA Y Y, YANG J Z. Switching the Richards' equation for modeling soil water movement under unfavorable conditions[J]. Journal of Hydrology, 2018, 563: 942–949. doi: 10.1016/j.jhydrol.2018.06.069
[18]	ZHA Y Y, YANG J Z, YIN L H, et al. A modified Picard iteration scheme for overcoming numerical difficulties of simulating infiltration into dry soil[J]. Journal of Hydrology, 2017, 551: 56–69. doi: 10.1016/j.jhydrol.2017.05.053
[19]	FARTHING M W, OGDEN F L. Numerical solution of richards' equation: a review of advances and challenges[J]. Soil Science Society of America Journal, 2017, 81(6): 1257–1269. doi: 10.2136/sssaj2017.02.0058
[20]	CHÁVEZ-NEGRETE C, DOMÍNGUEZ-MOTA F J, SANTANA-QUINTEROS D. Numerical solution of Richards' equation of water flow by generalized finite differences[J]. Computers and Geotechnics, 2018, 101: 168–175. doi: 10.1016/j.compgeo.2018.05.003
[21]	DOLEJŠÍ V, KURAZ M, SOLIN P. Adaptive higher-order space-time discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[J]. Applied Mathematical Modelling, 2019, 72: 276–305. doi: 10.1016/j.apm.2019.02.037
[22]	CELIA M A, BOULOUTAS E T, ZARBA R L. A general mass-conservative numerical solution for the unsaturated flow equation[J]. Water Resources Research, 1990, 26(7): 1483–1496. doi: 10.1029/WR026i007p01483
[23]	LIST F, RADU F A. A study on iterative methods for solving Richards' equation[J]. Computational Geosciences, 2016, 20(2): 341–353. doi: 10.1007/s10596-016-9566-3
[24]	陈曦, 于玉贞, 程勇刚. 非饱和渗流Richards方程数值求解的欠松弛方法[J]. 岩土力学, 2012, 33(增刊1): 237–243. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2012S1039.htm CHEN Xi, YU Yu-zhen, CHENG Yong-gang. Under-relaxation methods for numerical solution of Richards equation of variably saturated flow[J]. Rock and Soil Mechanics, 2012, 33(S1): 237–243. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2012S1039.htm
[25]	李文涛, 马田田, 韦昌富. 基于自适应松弛Picard法的高效非饱和渗流有限元分析[J]. 岩土力学, 2016, 37(1): 256–262 https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201601031.htm LI Wen-tao, MA Tian-tian, WEI Chang-fu. An efficient finite element procedure for unsaturated flow based on adaptive relaxed Picard method[J]. Rock and Soil Mechanics, 2016, 37(1): 256–262. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201601031.htm
[26]	BENZI M. Preconditioning techniques for large linear systems: a survey[J]. Journal of Computational Physics, 2002, 182(2): 418–477. doi: 10.1006/jcph.2002.7176
[27]	WANG K, ZHANG J. MSP: a class of parallel multistep successive sparse approximate inverse preconditioning strategies[J]. SIAM Journal on Scientific Computing, 2003, 24(4): 1141–1156. doi: 10.1137/S1064827502400832
[28]	LIU C Y, KU C Y, HUANG C C, et al. Numerical solutions for groundwater flow in unsaturated layered soil with extreme physical property contrasts[J]. International Journal of Nonlinear Sciences and Numerical Simulation, 2015, 16(7/8): 325–335.
[29]	BRIGGS W L, HENSON V E, MCCORMICK S F. A Multigrid Tutorial, Second Edition[M]. Philadelphia: Society for Industrial and Applied Mathematics, 2000.
[30]	IVERSON R M. Landslide triggering by rain infiltration[J]. Water Resources Research, 2000, 36(7): 1897–1910. doi: 10.1029/2000WR900090
[31]	GARDNER W R. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[J]. Soil Science, 1958, 85(4): 228–232. doi: 10.1097/00010694-195804000-00006

施引文献(15)

期刊类型引用(13)

1.	严辉，林沛元. 深圳市岩溶地层标准贯入击数神经网络模型. 地质科技通报. 2025(02): 305-321 . 百度学术
2.	李书兆，孙国栋，申辰，李文逵，罗进华，王教龙. 基于深度学习和地震数据的海上风电场CPT预测研究. 工程地球物理学报. 2025(02): 216-226 . 百度学术
3.	崔纪飞，柏林，饶平平，康陈俊杰，张锟. 基于人工智能算法的氯盐侵蚀混凝土预测模型. 硅酸盐通报. 2024(02): 439-447 . 百度学术
4.	段文魁，王来发，晁华俊，明锋. 冻结过程中土体导热系数预测模型. 中国农村水利水电. 2024(05): 47-52 . 百度学术
5.	唐少容，殷磊，杨强，柯德秀. 微胶囊相变材料改良粉砂土的导热系数及预测模型. 中国粉体技术. 2024(03): 112-123 . 百度学术
6.	姚兆明，王洵，齐健. 土体导热系数智能方法预测及影响因素敏感性分析. 工程热物理学报. 2024(05): 1440-1449 . 百度学术
7.	邓志兴，谢康，李泰灃，王武斌，郝哲睿，李佳珅. 基于粗颗粒嵌锁点高铁级配碎石振动压实质量控制新方法. 岩土力学. 2024(06): 1835-1849 . 百度学术
8.	李林，左林龙，胡涛涛，宋博恺. 基于孔压静力触探试验测试数据的原位固结系数物理信息神经网络反演方法. 岩土力学. 2024(10): 2889-2899 . 百度学术
9.	王红旗，李栋伟，钟石明，贾志文，王泽成，陈鑫，秦子鹏. 石灰改良红黏土导热系数影响因素及模型预测. 科学技术与工程. 2023(05): 2084-2092 . 百度学术
10.	王才进，武猛，蔡国军，赵泽宁，刘松玉. 基于多元分布模型预测土体热阻系数. 岩石力学与工程学报. 2023(S1): 3674-3686 . 百度学术
11.	王健翔，任瑞琪. 电学等效的稳态平板导热系数测试实验装置. 电子制作. 2023(11): 105-109 . 百度学术
12.	王才进，武猛，杨洋，蔡国军，刘松玉，何欢，常建新. 基于生物地理优化的人工神经网络模型预测软土的固结系数. 岩土力学. 2023(10): 3022-3030 . 百度学术
13.	徐明，康雅晶，马斯斯，张鹤. 基于贝叶斯优化的XGBoost模型预测路基回弹模量. 公路交通科技. 2023(11): 51-60 . 百度学术