轴对称承压井瞬态渗流问题

关键词： 轴对称 瞬态 线性 渗流 向后差分

本节以轴对称的承压井瞬态渗流问题为例，讲解如何使用 FEtch 系统对柱坐标系下的线性抛物型方程进行求解。

控制方程

微分形式

取柱坐标系 $Or\theta z$，对称轴为 $z$ 轴，径向为 $r$ 轴，环向为 $\theta$ 轴。对于轴对称问题，未知量与 $\theta$ 无关，只与 $r$ 和 $z$ 有关，因此可只研究 $rz$ 平面上的截面部分。

对于承压井瞬态渗流问题，满足如下偏微分方程 \(s \frac{\partial u}{\partial t}-\frac{1}{r}\frac{\partial}{\partial r}\left(k r\frac{\partial u}{\partial r}\right)-\frac{\partial}{\partial z}\left(k \frac{\partial u}{\partial z}\right)=0 \ \left(\text{in } \Omega\right) \tag{1}\label{eq1}\) 其中，$s$ 为材料的储水系数，$u$ 表示压力水头，$k$ 是渗透系数。根据需要，结合初始条件和边界条件，$u$ 同时可以作为增量形式的降深（$du=u_0-u$）来使用。

这里考虑两类边界条件:

第一类边界条件 \(u=\bar{u} \ \left(\text{on } \Gamma_1\right)\)

第二类边界条件 \(-k\frac{\partial u}{\partial \mathbf{n}} \cdot \mathbf{n}=-\bar{q}\ \left(\text{on } \Gamma_2\right)\)

其中，$\bar{u}$ 和 $\bar{q}$ 为边界上的水头和流量值。$\mathbf{n}$ 为边界上的单位外法向量。水流 $\bar{q}$ 以进入区域为正。

特别地，对于大范围求解域问题，井径可以视为无限小，因此经常被假设为线状源汇，此时第二类边界条件退化为 \(-k\frac{\partial u}{\partial r}r\bigg|_{r \rightarrow 0}= -\bar{Q}\ \left(\text{on } \Gamma_3\right)\) 其中， $\bar{Q}=\frac{Q}{2\pi H}$ ，$Q$ 为井的流量，$H$ 为井深。

弱形式

运用虚位移原理求水头，由微分方程 $\eqref{eq1}$ 两端乘以压力水头的变分，并在全域上积分，得 \(2\pi\int_{\Omega}\left[s \frac{\partial u}{\partial t}-\frac{1}{r}\frac{\partial}{\partial r}\left(k r\frac{\partial u}{\partial r}\right)-\frac{\partial}{\partial z}\left(k \frac{\partial u}{\partial z}\right)\right] \delta u\ rdrdz=0\) 对于二阶导数的部分，利用分部积分将其中的一阶导数转化到压力水头的变分上，得 \(\int_{\Omega} s \frac{\partial u}{\partial t} \delta u\ rdrdz+\int_{\Omega}\left(k \frac{\partial u}{\partial r} \frac{\partial \delta u}{\partial r}+k \frac{\partial u}{\partial z} \frac{\partial \delta u}{\partial z}\right) rdrdz=\int_{\Gamma} k \frac{\partial u}{\partial \mathbf{n}} \cdot \mathbf{n} \delta u \ rd \Gamma\) 对于分部积分得到的边界积分项，代入边界条件，得 \(\int_{\Omega} s \frac{\partial u}{\partial t} \delta u\ rdrdz+\int_{\Omega}\left(k \frac{\partial u}{\partial r} \frac{\partial \delta u}{\partial r}+k \frac{\partial u}{\partial z} \frac{\partial \delta u}{\partial z}\right) rdrdz\\=\int_{\Gamma_2} \bar{q} \delta u \ rd \Gamma+\int_{\Gamma_3} \bar{Q} \delta u d \Gamma \tag{2}\label{eq2}\) 这就是压力水头场的最终的弱解形式。

如果与直角坐标系下的弱解形式进行对比，可以发现唯一的不同点在于柱坐标系下的积分项与径向坐标 $r$ 相关，因此，程序编制的关键是如何将径向坐标作为参数引入到单元矩阵的计算中。

离散化

对式 $\eqref{eq2}$ 进行空间离散化后，可进一步写成矩阵形式 \(M \dot{U}+SU=F\) 其中， $U$ 为水头向量，$\dot{U}$ 为水头对时间的导数向量，$M$ 为质量矩阵，$S$ 为刚度矩阵，$F$ 为载荷向量。采用向后差分法进行时域离散，可将上式改写为 \(M \frac{U^{t+\Delta t}-U^{t}}{\Delta t}+S U^{t+\Delta t}=F^{t+\Delta t}\) 整理得 \((M+S \Delta t) U^{t+\Delta t}=F^{t+\Delta t} \Delta t+M U^{t} \tag{3}\label{eq3}\) 上标 $t+\Delta t$ 代表当前时间步，$t$ 代表上一时间步。该式即为最终需要求解的线性方程组。

脚本编写

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

PDE 和 FBC 文件

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

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

同样地，将方程 $\eqref{eq2}$ 的边界积分项写入 FBC 文件如下。

需要注意 FBC 文件中全局坐标变量 gx 的定义和使用。这里通过 coef 段来引入全局坐标。对于柱坐标系（1dr、2drz、3droz）和球坐标系（3dros）的 FBC 文件，系统默认 coef 段的最后 n 项为全局坐标，其中 n 为空间维度。进一步地，对于非线性和耦合问题， coef 段还会用来引入其他参数，这些参数需要全部置于全局坐标变量之前。

MDI 文件

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

注意，同一个场的 PDE 和 FBC 的单元类型和积分方法需要保证匹配，否则会导致计算精度降低，甚至计算失败的情况发生。

MDI 文件声明的单元类型和积分方式将以参数的形式传递给 PDE 和 FBC 文件中的的 %1、%2 和 %3。如果后续想要使用其他的单元类型或者积分方法，只需修改 MDI 文件即可。

GCN 文件

GCN 文件给出计算流程信息，系统规定其文件名字要与 MDI 文件一致。线性抛物方程计算流程比较简单，包含一个时间步的循环过程。只需填写少量代码，对应的 GCN 文件如下。

NFE 文件

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

算例

Theis 非稳定完整井流

Theis 问题是地下水动力学的经典算例，描述了无限含水层中单个定流量井的渗流过程。它的假定条件为：
①含水层是均质、各向同性、等厚且水平分布，含水层假定为弹性体；
②无垂向补给、排泄；
③渗流满足 Darcy 定律；
④完整井，假定流量沿井壁均匀进水；
⑤水头下降引起地下水从储存量中的释放是瞬时完成的；
⑥抽水前水头面是水平的；
⑦井径无限小且定流量抽水；
⑧含水层侧向无限延伸。

对于这一典型的轴对称问题，这里考虑半径 $R=1000\ \mathrm{m}$，深 $H=10\ \mathrm{m}$ 的含水层，研究抽水井引起的地下水位的瞬变下降情况。模型参数为 $s=10^{-3}$，$k=8\ \mathrm{m/d}$，$\bar{Q}=100\ \mathrm{m^2/d}$。水头初始值为 $0$。试求开采 $3\ \mathrm{h}$ 后水头场的分布。

前处理

填写材料参数文件

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

填写时间步长文件

由 seep.gcn 文件可知，用户自定义了时间步长文件 time.dat，给出时间步长和终止时间。文件内容如下：

0.001 0.125

可知，时间步长为 0.001 d，计算终止时间为 0.125 d （=3 h）。

赋单元材料

赋初始条件

施加第二类边界条件

网格剖分

根据所开发程序的计算要求，采用四边形线性单元进行网格剖分。为保证计算精度，对井附近进行了局部加密。

注：由于模型的长宽比例过大，以上所有截图仅着重显示了模型的左侧部分。

计算结果

水头场的演化过程

上图给出了三个不同观察点处水头随时间的演化曲线。其中，实线为解析解，散点为计算结果。可以发现数值解与解析解符合得很好，充分证明了算法和程序的有效性。