非饱和渗流问题的有限元求解

学术科技 2024-10-11 15:54 山东

点击上方「蓝字」关注我们

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

控制方程

微分形式

Richards 方程是非饱和渗流理论的基本方程。我们以孔隙水压力水头作为基本变量，，为孔隙水压力，为水的容重；以表示竖直方向的坐标。如果考虑重力作用，则还代表了位置水头。此时，总水头为。在二维直角坐标系下，Richards 方程可以表示为

\begin{aligned} 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) \end{aligned}

式中，为多孔介质的孔隙率，为饱和度，为渗透系数，为水平方向的坐标，为时间。如果不考虑重力作用，则上式左侧最后一项为。

在初始时刻，压力水头为，即初始条件为

h=h_0 \ \left(\text{in } \Omega\right)

假设在边界上压力水头为定值，在边界上流量为定值 (流入为正)，则边界条件可表示为：

强制边界条件

h=\bar{h}\ \left(\text{on } \Gamma_{h}\right)

自然边界条件

-K\frac{\partial(h+z)}{\partial \mathbf{n}} \cdot \mathbf{n}=-\bar{q} \ \left(\text {on } \Gamma_q\right)

其中，为多孔介质表面的单位法向量。

弱形式

运用变分原理，由以上方程得

\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

将其转化为弱形式，可得

\begin{aligned} &\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}} \delta h d \Gamma \end{aligned}

代入边界条件，可得

\begin{aligned} &\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 \end{aligned}

离散化

对上式进行空间离散化后，可写成矩阵形式

M(U) \dot{U}+S(U)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}

其中，上标代表当前时间步，代表上一时间步。和均为的函数，因此，上式是一个非线性方程组，无法直接求解，需要引入迭代过程。

最常用的处理方法是采用直接迭代法（也称 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}

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

材料本构模型

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

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

自然指数模型

S_{r}(h)=S_{r}^{irr}+\left(1-S_{r}^{irr}\right) e^{\beta h}

K(h) = K_s e^{\beta h}

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

van Genuchten 模型

S_{r}(h)=S_{r}^{irr}+\left(1-S_{r}^{i r r}\right)\left[1+\left(-\alpha_{v} h\right)^{n_{v}}\right]^{-m_{v}}

K(h)=K_s \sqrt{S_{e}}\left[1-\left(1-S_{e}^{1 / m_{v}}\right)^{m_{v}}\right]^{2}

式中，有效饱和度，、和是与土性相关的模型参数，通常取。

算例1

一维均质土的非饱和渗流

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

前处理

进入 GiD 前处理，建立一个的长方形区域，它的两个对角的角点坐标分别为（0, 0）和（0.1, 1）。网格剖分采用 4 节点等参单元。为了在保证计算精度的同时尽量减少网格数量，从下到上按 1.1:1 的比例进行网格局部加密，共使用了 40 个单元，82 个节点。

计算结果

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

压力水头的演化过程

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

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

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

算例2

二维砂槽的非饱和渗流

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

前处理

进入 GiD 前处理，建立一个的长方形区域。采用 4 节点等参单元，将该区域剖分为 750 个双线性四边形等参单元，共 806 个结点。

计算结果

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

压力水头的演化过程

部分等值线的演化过程

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

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

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

选取入渗前沿的点 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) 作为观察点，压力水头随时间变化的情况如下图所示。

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

往期推荐

一图读懂 FEtch 的工作流程

谈谈有限元程序开发，兼谈 FEtch 入门

通用前后处理软件 GiD 简介

非饱和渗流问题的有限元求解

控制方程

微分形式

弱形式

离散化

材料本构模型

算例1

前处理

计算结果

算例2

前处理

计算结果

推荐阅读