算例-非饱和渗流

非饱和渗流问题

关键词：非饱和 渗流 瞬态 非线性 向后差分 迭代 变时间步长

许多岩土工程问题会涉及到非饱和渗流过程，如降雨入渗或地下水变化时土质边坡与堤坝的稳定性评价；垃圾填埋场内部污染物质的运移模拟；冻土中相变发生时的渗流过程分析；高放核废料的地质深埋处理等。因此，有效地模拟和分析非饱和渗流过程有着重要的实际意义，一直是岩土工程、水利工程、环境工程中的一项热门课题。

控制方程

微分形式

Richards 方程是非饱和渗流理论的基本方程。以孔隙水压力水头 $h$ 作为基本变量，$h =p/γ_w$，$p$ 为孔隙水压力，$γ_w$ 为水的容重。以 $z$ 表示竖直方向的坐标。如果考虑重力作用，则 $z$ 还代表了位置水头。此时，总水头为 $h+z$。二维坐标系下的 Richards 方程可表示为 \(n\frac{\partial S_{r}}{\partial h} \frac{\partial h}{\partial t}-\frac{\partial}{\partial x}\left(K \frac{\partial h}{\partial x}\right)-\frac{\partial}{\partial z}\left(K \frac{\partial h}{\partial z}\right)-\frac{\partial K}{\partial z}=0 \ \left(\text{in } \Omega\right) \tag{1}\label{eq1}\) 式中，$n$ 为多孔介质的孔隙率，$S_r$ 为饱和度，$K$ 为渗透系数，$x$ 为水平方向的坐标，$t$ 为时间。如果不考虑重力作用，则上式左侧最后一项 ${\partial K}/{\partial z}$ 为 $0$。

在初始时刻，压力水头为 $h_0$，即初始条件为 \(h=h_0 \ \left(\text{in } \Omega\right)\)

在边界 $\Gamma_h$ 上压力水头为定值 $\bar{h}$，在边界 $\Gamma_q$ 上流量为定值 $\bar{q}$ (流入为正)，则边界条件可表示为：

1. 强制边界条件 \(h=\bar{h}\ \left(\text{on } \Gamma_{h}\right)\)
2. 自然边界条件 \(-K\frac{\partial(h+z)}{\partial \mathbf{n}} \cdot \mathbf{n}=-\bar{q} \ \left(\text {on } \Gamma_q\right)\)

其中，$\mathbf{n}$ 为多孔介质表面的单位法向量。

弱形式

运用变分原理，由方程 $\eqref{eq1}$ 得 \(\int_{\Omega}\left[n\frac{\partial S_{r}}{\partial h} \frac{\partial h}{\partial t}-\frac{\partial}{\partial x}\left(K \frac{\partial h}{\partial x}\right)-\frac{\partial}{\partial z}\left(K \frac{\partial h}{\partial z}\right)-\frac{\partial K}{\partial z}\right] \delta u d \Omega=0\) 将其转化为弱形式，可得 \(\int_{\Omega} n\frac{\partial S_{r}}{\partial h} \frac{\partial h}{\partial t} \delta h d \Omega+\int_{\Omega}\left(K \frac{\partial h}{\partial x} \frac{\partial \delta h}{\partial x}+K \frac{\partial h}{\partial z} \frac{\partial \delta h}{\partial z}\right) d \Omega \\=-\int_{\Omega} K \frac{\partial \delta h}{\partial z} d \Omega+\int_{\Gamma} K \frac{\partial (h+z)}{\partial \mathbf{n}}\cdot \mathbf{n} \delta h d \Gamma\) 代入边界条件，可得 \(\int_{\Omega} n\frac{\partial S_{r}}{\partial h} \frac{\partial h}{\partial t} \delta h d \Omega+\int_{\Omega}\left(K \frac{\partial h}{\partial x} \frac{\partial \delta h}{\partial x}+K \frac{\partial h}{\partial z} \frac{\partial \delta h}{\partial z}\right) d \Omega \\=-\int_{\Omega}K \frac{\partial \delta h}{\partial z} d \Omega+\int_{\Gamma_q} \bar{q} \delta h d \Gamma \tag{2}\label{eq2}\)

离散化

对式 $\eqref{eq2}$ 进行空间离散化后，可写成矩阵形式 \(M(U) \dot{U}+S(U)U=F\) 其中 $U$ 为压力水头向量，$\dot{U}$ 为压力水头对时间的导数向量，$M(U)$ 为质量矩阵， $S(U)$ 为刚度矩阵，$F$ 为载荷向量。采用向后差分法进行时域离散，可将上式改写为 \(M(U) \frac{U^{t+\Delta t}-U^{t}}{\Delta t}+S(U) U^{t+\Delta t}=F^{t+\Delta t}\) 整理得 \(\left(M(U)+S(U) \Delta t\right) U^{t+\Delta t}=F^{t+\Delta t} \Delta t+M(U) U^{t}\) 其中，上标 $t+\Delta t$ 代表当前时间步，$t$ 代表上一时间步。$M(U)$ 和 $S(U)$ 均为 $U$ 的函数，因此，上式是一个非线性方程组，无法直接求解，需要引入迭代过程。

最常用的处理方法是采用直接迭代法（也称 Picard 法），进行线性化，得到 \(\left(M(U^i)+S(U^i) \Delta t\right) U^{i+1,t+\Delta t}=F^{t+\Delta t} \Delta t+M(U^i) U^{t}\)

其中，上标 $i$ 表示 $t+\Delta t$ 时间步下的迭代步数。上式在每个时间步下都需要反复地迭代计算，直至误差 $\Vert U^{i+1}-U^i\Vert$ 小于某一设定的阈值。直接迭代法可以通过使用松弛因子，大幅提高收敛速度。松弛迭代法已被写入 NFE 算法模板库， 是 FEtch 系统默认的非线性迭代算法。

材料本构模型

非线性问题远比线性问题要复杂。由于非线性参数的函数形式千变万化，所以通常需要针对具体的问题专门地设计程序。对于 FEtch 系统而言，由于脚本文件是完全自主可控的，只要基本框架搭好，后面需要改动的地方极少。

孔隙水压力和相对渗透系数是饱和度的函数，它们之间的关系通常需要根据试验数据进行拟合，目前已有多种通用模型可以描述，如广泛使用的 van Genuchten 模型、Brooks-Corey 模型等。为了后续算例的计算需要，这里考虑以下两种模型。

1.自然指数模型

这里将持水特征曲线和渗透系数函数假设为压力水头 $h$ 的自然指数形式，主要应用于求解析解。式中，$S_{r}^{irr}$ 为残余饱和度，$\beta$ 为与土性相关的模型参数，$K_s$ 为饱和渗透系数。

2.van Genuchten 模型

式中，有效饱和度 $S_{e}=\left(S_{r}-S_{r}^{irr}\right) /\left(1-S_{r}^{irr}\right)$，$\alpha_v$ 、$n_v$ 和 $m_v$是与土性相关的模型参数，通常取 $m_v=1-1/n_v$。

后续程序开发的重点是定义和使用函数 $S_r(h)$ 和 $K(h)$ 。用户需认真体会有限元语言是如何实现的。

脚本编写

一个具体问题的有限元语言表述，通常需要包含 5 种类型的文件：PDE、FBC、MDI、GCN 和 NFE。下面依次给出详细介绍。

PDE 和 FBC 文件

在这个例子中，只有一个物理场，需要给出孔压场对应的 PDE 文件和 FBC 文件。根据有限元语言的语法规则，将方程 $\eqref{eq2}$ 的体积分项写入 PDE 文件如下。

这里将 shap 行和 gaus 行的参数设为待定，参数值 %1、%2 和 %3 将从 MDI 文件中取得。

flg 参数用来标记是否考虑重力场。考虑时 flg 取值为1.0， 否则为 0.0。

文件中调用的自定义函数，由另外创建的 Fortran 源代码文件 swcc.f90 给出，具体文件内容如下。

! swcc.f90
! Exp model
real*8 function sat1(h,sir,beta)
implicit real*8 (a-h,o-z)     
    sat1=sir+(1.0-sir)*dexp(beta*h)
end

real*8 function deri1(h,sir,beta)
implicit real*8 (a-h,o-z)
    deri1 = (1.0-sir)*beta*dexp(beta*h)   
end

real*8 function seep1(h,beta)
implicit real*8 (a-h,o-z)
    seep1 = dexp(beta*h)
end

! VG model      
real*8 function sat2(h,sir,ea,en,em)
implicit real*8 (a-h,o-z)
    sat2=sir+(1.0-sir)*(1.0+(-ea*h)**en)**(-em)  
end

real*8 function deri2(h,sir,ea,en,em)
implicit real*8 (a-h,o-z)
    deri = em*en*ea*(1.0-sir)*(-ea*h)**(en-1)
    deri2 = deri*(1.0+(-ea*h)**en)**(-em-1)
end

real*8 function seep2(h,sir,ea,en,em)
implicit real*8 (a-h,o-z)
    sr = sat2(h,sir,ea,en,em)
    se = (sr-sir)/(1.0-sir)
    seep2 = se**0.5*(1.0-(1.0-se**(1.0/em))**em)**2
end

类似 PDE 文件的填写，将方程 $\eqref{eq2}$ 的边界积分项写入 FBC 文件如下。

MDI 文件

MDI 文件负责组织 PDE 和 FBC 文件，声明程序具体使用的坐标系统、单元类型和积分方法。这里文件命名为 unsat.mdi，具体内容如下。

注意，同一个场的 PDE 和 FBC 的单元类型和积分方法需要保证匹配，否则会导致计算精度降低，甚至计算失败的情况发生。MDI 文件声明的单元类型和积分方式将以参数的形式传递给 PDE 和 FBC 文件中的的 %1、%2 和 %3。如果后续想要使用其他的单元类型或者积分方法，只需修改 MDI 文件即可。

GCN 文件

GCN 文件给出计算流程信息，系统规定其文件名字要与 MDI 文件一致。非线性抛物方程的计算流程较为复杂，包含两大循环过程，时间步循环和迭代步循环。时间步长可以根据非线性迭代的收敛速度进行动态调整。对应的 GCN 文件和说明如下。

NFE 文件

GCN 文件使用了 nparb 算法，系统会基于模板文件自动生成以下的 npa.nfe 文件。如果用户有其他需要，可以在此基础上进一步地修改。

算例1

一维均质土的非饱和渗流

考虑无限长的均质土层，设其厚度为 $L = 1.0\ \mathrm{m}$，底部对应地下水位（$ h=0$），顶部为入渗边界，以入渗流量为 0 时的稳定状态为初始状态，时间 $t>0$ 时入渗流量变为 $q_A$。持水特征曲线和渗透系数函数采用自然指数模型 $\eqref{eq3a}$ 和 $\eqref{eq3b}$ 来描述。模拟过程中的具体参数取值如下表所示。

前处理

建模

进入 GiD 前处理，建立一个 $0.1\ \mathrm{m}\times1.0\ \mathrm{m}$ 的长方形，它的两个对角的角点坐标分别为（0, 0）和（0.1, 1）。

填写材料参数文件

系统自动生成了unsat.mat 文件，用于设置材料参数。参数初始值全部为 PDE 和 FBC 文件给出的默认值。可根据具体问题进行修改。对于本算例，修改完成后的文件内容如下。

填写时间步长文件

由 unsat.gcn 文件可知，用户自定义了时间步长文件 time.dat，给出时间步长和终止时间。mat 文件中与时间相关的参数以小时 $\mathrm{h}$ 为单位， time.dat 与其保持一致。文件内容如下：

0.001 32

可知， 初始时间步长为 $0.01\ \mathrm{h}$，结束时间为 $32\ \mathrm{h}$。

赋单元材料

赋初始条件

注意，初始条件与几何模型的位置是直接相关的。这里的长方形的底边（同时也是水位线位置）在 x 轴上，因此初始压力水头按 -y 来分布。如果模型位置发生变化，初始压力水头要做相应的换算后，再进行施加。

施加第一类边界条件

其中：
u-I 为边界类型，-1 表示约束，1 表示自由
u-D 为相应的约束值或施加节点流量

施加第二类边界条件

第二类边界条件以线单元的形式施加

网格剖分

根据所开发程序的计算要求，网格剖分采用 4 节点等参单元。为了在保证计算精度的同时尽量减少网格数量，从下到上按 1.1:1 的比例进行网格局部加密，共使用了 40 个单元，82 个节点。

计算结果

通过自动改变时间步长，变化范围从 $0.001\ \mathrm{h}$ 到 $0.512\ \mathrm{h}$，总共使用了 183 个时间步。对计算结果稍加整理，效果如下。

压力水头的演化过程

不同时刻压力水头的空间分布

不同高度处压力水头随时间的变化

可以看到与解析解符合得很好。同时无论在时间和空间域内，模拟结果都没有震荡现象发生，充分证明了算法和程序的有效性。

算例2

二维砂槽的非饱和渗流

Vauclin 等进行的二维室内渗流试验也是检验算法的典型算例。流动区域为 $6.00\ \mathrm{m}\times2.00\ \mathrm{m}$ 的厚 $0.05\ \mathrm{m}$ 的砂土，底部为不透水边界，两侧为自由排水边界。土体为颗粒分布很规则的河砂。初始水位为 $0.65\ \mathrm{m}$，试验开始后在土槽顶部中间 $1.00\ \mathrm{m}$ 的范围均匀施加 $0.148\ \mathrm{m/h}$ 的降雨强度，共 $8\ \mathrm{h}$。试验中对流动域内自由水面位置的变化等数据进行了测量。该类砂土的孔隙水压力水头和渗透系数与饱和度之间的关系由 van Genuchten 模型 $\eqref{eq4a}$ 和 $\eqref{eq4b}$ 来描述。模拟过程中的具体参数取值如下表所示。根据对称性，从中线截取区域的右半部分进行计算。

前处理

填写材料参数文件

对于本算例，材料参数文件 unsat.mat 的内容如下。

填写时间步长文件

由 unsat.gcn 文件可知，用户自定义了时间步长文件 time.dat，给出时间步长和终止时间。mat 文件中与时间相关的参数以小时 $\mathrm{h}$ 为单位， time.dat 与其保持一致。文件内容如下：

0.001 8

可知， 初始时间步长为 $0.01\ \mathrm{h}$，结束时间为 $8\ \mathrm{h}$。

赋单元材料

赋初始条件

施加第一类边界条件

其中：
u-I 为边界类型，-1 表示约束，1 表示自由
u-D 为相应的约束值或施加节点流量

施加第二类边界条件

第二类边界条件以线单元的形式施加

网格剖分

剖分为 750 个双线性四边形等参单元，共 806 个结点。

计算结果

通过自动改变时间步长，变化范围从 $0.001\ \mathrm{h}$ 到 $0.256\ \mathrm{h}$，总共使用了 202 个时间步。对计算结果稍加整理，效果如下。

压力水头的演化过程

部分等值线的演化过程

不同时刻地下水位位置模拟结果与试验结果对比

通过不同时刻地下水位位置（$h =0$ 处）模拟结果与试验结果的对比，可以看出两者吻合较好。

不同观察点的压力水头随时间的变化

选取入渗前沿的点 A (0, 2.0)、B (0,1.6)、C (0, 1.2)、D (0.5, 1.6)、E (0.5, 1.2) 和 F (1.0,1.2) 作为观察点，压力水头随时间变化的情况如下图所示。

可以看到随着入渗的进行，孔隙水压力逐渐上升，向正孔压发展，且整个变化过程是平滑和稳定的，进一步证明了算法和程序的有效性。