Multiscale finite element method–triple grid model for simulation of groundwater flows
-
摘要: 传统有限元法在模拟地下水流问题时常需要精细剖分描述含水介质的非均质性以保证精度，导致计算消耗过高。传统多尺度有限元法（MSFEM）能缓解这一问题，但在处理高计算量问题时仍需较高消耗来构造基函数。提出了一种用于模拟地下水流的三重尺度有限元模型（MSFEM-T），该方法在MSFEM的粗、细两种尺度网格之间引入中网格，从而可以在粗网格单元内基于中、细两种尺度网格应用MSFEM本身替代有限元法构造基函数，能够显著降低基函数的构造消耗以提高整体计算效率。此外，MSFEM-T还提出了一种基于粗、中、细三重网格的超样本技术，可以进一步提升其计算精度。数值算例结果显示MSFEM-T的精度与MSFEM和精细剖分有限元法（LFEM-F）的精度相近，且计算效率更高。Abstract: The traditional finite element method often requires fine element grids to describle the heterogeneity of medium to ensure the accuracy for numerical modeling of groundwater, which leads to a large amount of calculation consumption. The multiscale finite element method can alleviate this problem, but it still needs a high cost to formulate the basis function when dealing with high computational complexity. A multiscale finite element method-triple grid model (MSFEM-T) is proposed for the simulation of groundwater flows. The MSFEM-T introduces an intermediate grid between the coarse grid and the fine grid, so that the basis function in the coarse grid can be established using the MSFEM instead of the FEM based on the intermediate and fine grids, therefore reducing the construction consumption of the basis function and improving the overall calculation efficiency. Moreover, the MSFEM-T uses an over-sampling method based on the coarse, intermediate and fine grids, which can further improve its calculation accuracy. The results show that the accuracy of the MSFEM-T is similar to that of the MSFEM and the finite element method of fine elements (LFEM-F), but the computational efficiency is much higher.
-
0. 引言
地下水作为水资源的重要组成部分，不合理的利用地下水资源会引起地面沉降、海（咸）水入侵及地下水污染等一系列问题[1-2]。有限元法作为求解地下水问题有效的数值方法[3]，但是应用有限元法求解大尺度非均质地下水问题时，由于网格不能骑跨岩层分界线，常需要精细剖分以保证精度，导致未知项过多而产生大量的计算消耗[4]。在模拟大尺度、长时间和非均质等复杂特性地下水问题时，有限元法需要巨大的计算成本，对硬件也有较高的需求[5]。因此，开发高效、精确的数值模拟方法对于中国岩土、水利工程中土壤及地下水问题的分析具有十分重要的意义[6]。
多尺度有限元法（Multiscale Finite Element Method，缩写MSFEM）是一种新型地下水数值模拟方法，已被广泛应用于断裂力学、结构健康监测、地下水流模拟、材料科学等领域[7-9]。MSFEM通过在粗网格上满足局部椭圆型微分算子来构造粗尺度基函数，能够有效获取水头的宏观特性，从而提升尺度以降低模拟所需的计算消耗。近年来许多科学家们致力于MSFEM的研究与改进工作，以提高该方法的效果与应用性，发展了广义多尺度有限元法[10]，混合多尺度有限元法[11]，多尺度有限体积法[12-14]，多尺度有限元降基法[15]，多尺度概率法[16]，多尺度有限元–有限元模型[17]等多种算法。同时，国内学者也应用MSFEM成功进行了上海和苏锡常地区的地面沉降的模拟[18]。然而，MSFEM在模拟高计算量问题时仍需大量计算消耗用于构造基函数[19-20]。随着经济发展，我国对地下水问题的时、空尺度和模拟精度要求与日俱增，MSFEM基函数较低的构造效率将限制其在水文地质领域的发展与应用。
2004年，Ye等[8]的研究证明MSFEM在模拟地下水流问题时，解的精度更依赖基函数边界条件而非基函数本身。本文的前期工作中[20]，谢一凡等基于三重网格应用区域分解技术替换有限元法进行了MSFEM的基函数构造，提高了构造效率。然而，该工作[20]需要先将基函数的构造问题分解为若干子问题，并需对粗网格单元进行相应的区域分解，而且要预先求解各子问题的定解条件，步骤较繁琐；该网格剖分的灵活度也因此受限。此外，该工作[20]也没有设计相应超样本技术来处理物理尺度和网格尺度之间因共振效应引起的误差，限制了解的精度。
本文提出了一种新型的三重尺度有限元模型(Multiscale Finite Element Method-Triple Grid Model，缩写MSFEM-T)来模拟地下水流运动问题，具有极高的精度与效率。与前期工作[20]不同，MSFEM-T在粗网格单元上使用MSFEM本身构造基函数，可直接求解退化的椭圆问题，无需为基函数预先规划子问题的区域和定解条件，十分便捷简单。同时，MSFEM-T在粗网格单元上应用了非常灵活的多尺度网格，中、细尺度网格的数目和位置均不受限制。此外，本文提出了基于三重网格的超样本技术，使得MSFEM-T能够满足高精度模拟的需求。数值模拟结果显示MSFEM-T的精度与MSFEM、LFEM-F相近，且具有更高的计算效率；在使用超样本技术后，MSFEM-T能够获得比LFEM-F更高的精度，但计算成本极低。
1. 原理及算法
1.1 MSFEM-T网格剖分方式
MSFEM-T的研究区网格剖分与MSFEM相同（图 1），但具有不同的粗网格单元的剖分方式。MSFEM在剖分三角形粗网格单元时直接将其剖分为细网格单元（图 2），而MSFEM-T粗网格的细分过程分为两个步骤：①将每个三角形粗网格单元划分为中网格单元；②再将每个中网格单元划分为细网格单元（图 3）。设研究区Ω被剖分为三角形粗网格单元Δijk。MSFEM-T的中网格单元为三角形单元Δabc，细网格单元为三角形单元Δaabbcc。
![]()
图 2 MSFEM粗网格剖分方式Figure 2. Sketch map of standard subdivision of coarse element of MSFEM![]()
图 3 MSFEM-T粗网格剖分方式Figure 3. Sketch map of standard subdivision of coarse element of MSFEM-T（1）步骤1：将粗网格单元Δijk剖分为中网格单元Δabc
如图 3（a）所示，MSFEM-T将三角形粗网格单元Δijk以规则的剖分方式划分为中网格单元，即等分Δijk各边界并连接形成中网格单元。
因此，在需要将粗网格单元划分为相同数目的细网格时，MSFEM-T和MSFEM将粗网格剖分的份数不同。设三角形粗网格Δijk需要被划分为144个三角形细网格单元。在图 2中，MSFEM在三角形粗网格单元Δijk上直接使用12×12的标准剖分方式来获得这144个三角形细网格单元。对于MSFEM-T，则先将Δijk划分为中网格，再将中网格划分为细网格，以应用MSFEM算法构造基函数。在图 3（a）中，MSFEM-T通过4×4的标准三角形单元剖分方式将粗网格分为16个三角形中网格单元。
（2）步骤2：将中网格单元Δabc划分为细网格单元Δaabbcc
在第一步划分Δijk为中网格后，MSFEM-T需要再将各中网格以规则剖分方式形成细网格。在图 3（a）中，Δijk已被划分为16个中网格。为了获得与MSFEM相同数目的144个细网格，MSFEM-T再将每个中网格划分为9个大小相同的细网格（图 3（b））。MSFEM-T的网格包含了多个尺度且具有很高的灵活性，可在应用MSFEM算法构造基函数时有效降低构造矩阵的阶数，提高基函数构造效率。
由于MSFEM-T与MSFEM均为规则剖分，可编写自动剖分程序，二者网格剖分的前处理的工作量和便捷性十分相近。
1.2 MSFEM-T基函数
与MSFEM不同的是，MSFEM-T的基函数不是直接使用有限元法求解，而是应用MSFEM算法进行求解。MSFEM-T的粗网格基函数能够高效提取各尺度网格的信息，并反映到粗尺度上，能够显著提高计算效率。MSFEM-T粗网格基函数的构建可以分为两步。
（1）步骤1：构造中网格单元基函数
首先在中网格Δabc上求解退化的椭圆型方程来构造中网格单元基函数φ。以Δabc上的基函数φa为例，退化的椭圆型方程为
\nabla\cdot \boldsymbol{K}(x,y)\nabla {\varphi _a} = 0 {\text{((}}x,y{\text{)}} \in {\Delta _{abc}})\text{，} (1) 式中， \boldsymbol{K} 为渗透系数张量， {\varphi _a} 为中网格顶点a的基函数。 {\varphi _a} 的边界条件可由基函数的线性、振荡型边界条件确定[5]。
应用伽辽金变换于式（1），可得
{J_{M{M_\tau }}} = \iint_{{\Delta _{abc}}} {\boldsymbol{K}(x,y)}\nabla {\varphi _a} \cdot \nabla {N_{M{M_\tau }}}{\text{d}}x{\text{d}}y = 0\\ \;\;\;\;\;\;\left(\tau=1,2, \cdots, p_{\mathrm{p}}\right)\text{，} (2) 式中，{N_{M{M_\tau }}}为未知节点M{M_\tau }处的线性基函数值，pp为未知节点总数。
根据MSFEM理论[5]，{\varphi _a}在Δabc内的每个细网格Δaabbcc上可被表示为
\begin{aligned} \varphi_a(x, y)= & \varphi_a(a a) N_{a a}+\varphi_a(b b) N_{b b}+\varphi_a(c c) N_{c c} \\ & \left((x, y) \in \Delta_{{aabbcc }}\right), \end{aligned} (3) 式中，{N_{aa}}，{N_{bb}}，{N_{cc}}分别为细网格Δaabbcc3个顶点aa，bb，cc处的有限元线性基函数。
将式（3）代入式（2），并离散到细网格单元上，经过化简、移项可得方程组{A_1}{\varphi _a} = {f_1}。式中，{A_1}为中网格上对称正定的刚度矩阵，{f_1}为右端自由项。求解可得到中网格单元基函数{\varphi _a}在中网格Δabc内所有结点的值。同理，{\varphi _b}, {\varphi _c}求解过程类似。
（2）步骤2：应用MSFEM算法构造粗网格基函数
在粗网格Δijk上求解退化的椭圆型方程：
\nabla \cdot \boldsymbol{K}(x, y) \nabla \mathit{\Psi}_i=0,(x, y) \in \triangle_{i j k}, (4) 式中， {\mathit{\Psi} _i} 为i点的粗网格基函数， {\mathit{\Psi} _i} 的边界条件可由基函数的线性、振荡型边界条件确定。
应用伽辽金变换于式（4），并离散到中网格单元上，可得
\begin{gathered} J_{M_\tau}=\sum \iint_{\Delta_{a b c}} K(x, y) \nabla \mathit{\Psi}_i \cdot \nabla \varphi_{M_\tau} \mathrm{d} x \mathrm{d} y=0 \\ (\tau=1,2, \cdots, p), \end{gathered} (5) 式中，{\varphi _{{M_\tau }}}为未知结点{M_\tau }处的中网格基函数，p为未知节点总数。
{\mathit{\Psi} _i}在每个中网格Δabc内可以被中网格基函数表示，{\mathit{\Psi} _i}可被表示为
\begin{gathered} \mathit{\Psi}_i(x, y)=\mathit{\Psi}_i(a) \varphi_a+\mathit{\Psi}_i(b) \varphi_b+\mathit{\Psi}_i(c) \varphi_c \\ \left((x, y) \in \triangle_{a b c}\right) \end{gathered}, (6) 式中，{\varphi _a}，{\varphi _b}，{\varphi _c}分别为中网格Δabc内3个顶点a，b，c处的中网格单元基函数。
将式（6）代入式（5），再离散到细网格单元上，并结合式（3）。经过化简、移项可得方程组{A_2}{\mathit{\Psi} _i} = {f_2}。式中{A_2}为粗网格上对称正定的刚度矩阵；{f_2}为右端自由项。求解可得到{\mathit{\Psi} _i}在粗网格Δijk内的所有结点的值。同理，{\mathit{\Psi} _j}, {\mathit{\Psi} _k}求解过程类似。
值得说明的是，MSFEM-T的剖分方法会减少边界节点的数目，会轻微弱化振荡边界条件的效果。本文的算例显示：在一般非均质条件下，MSFEM-T的剖分对于精度降低的程度很小，所得结果依然非常精确，已经可以适用于一般的地下水问题。此外，当K的高度振荡时，MSFEM-T可以调整中、细尺度单元的比例，来增加中尺度边界节点数量，从而解决这一问题。
1.3 MSFEM和MSFEM-T基函数的构造消耗
设粗网格Δijk剖分为{n_0}^2个三角形细网格单元。MSFEM需要在每个粗网格有\frac{{{n_0}^2 - 3{n_0} + 2}}{2}个内部节点，因此需要求解一个\frac{{{n_0}^2 - 3{n_0} + 2}}{2} \times \frac{{{n_0}^2 - 3{n_0} + 2}}{2}的方程组来构建一个粗网格上的基函数。
在MSFEM-T中，设需要将粗网格Δijk剖分为{N^2}个中网格单元，再将每个中网格单元剖分为N{N^2}个细网格单元（其中{n_0} = N \times NN），MSFEM-T需要在每个中网格中求解{N^2}个\frac{{N{N^2} - 3NN + 2}}{2} \times \frac{{N{N^2} - 3NN + 2}}{2}阶的方程组来构造中尺度基函数，再求解一个\frac{{{N^2} - 3N + 2}}{2} \times \frac{{{N^2} - 3N + 2}}{2}的方程组来构造粗尺度基函数。以1.1节的剖分为例，如图 2，3所示，当MSFEM和MSFEM-T均将粗网格剖分144个细网格时（ {n}_{0}=12，N=4，{N}_{\text{N}}=3 ），构造粗尺度基函数计算复杂度方面，MSFEM需求解1个55×55阶的方程组，而MSFEM-T仅需求解16个1×1的方程和1个3×3的方程组。
由于解方程组的运算数一般为方程阶数n的指数级，故MSFEM-T可以显著降低计算运算数和工作量。
1.4 超样本技术
模拟地下水流运动时，多孔介质的物理尺度和剖分的网格大小相近会导致共振效应，引发谐振误差[5]。本文改进了Hou等基于MSFEM的工作[5]，提出了适用于MSFEM-T三重网格的超样本技术。MSFEM-T在应用超样本技术时，需要先将原粗网格单元Δijk放大，再在放大的单元（超样本单元）上应用1.2节的过程构造临时的粗网格基函数{\mathit{\Phi} _i},{\mathit{\Phi} _j}和{\mathit{\Phi} _k}，最后用临时基函数表示原粗网格基函数，即
\left.\begin{array}{l}{\mathit{\Psi} }_{i}={C}_{11}{\mathit{\Phi} }_{i}+{C}_{12}{\mathit{\Phi} }_{j}+{C}_{13}{\mathit{\Phi} }_{k}\text{，}\\ {\mathit{\Psi} }_{j}={C}_{21}{\mathit{\Phi} }_{i}+{C}_{22}{\mathit{\Phi} }_{j}+{C}_{23}{\mathit{\Phi} }_{k}\text{，}\\ {\mathit{\Psi} }_{k}={C}_{31}{\mathit{\Phi} }_{i}+{C}_{32}{\mathit{\Phi} }_{j}+{C}_{33}{\mathit{\Phi} }_{k}\text{，}\end{array}\right\} (7) 式中，{C_{\alpha \beta }}为超样本系数(\alpha ,\beta= 1,2,3)。以{\mathit{\Psi} _i}为例，{C_{\alpha \beta }}满足
\left.\begin{array}{l}{\mathit{\Psi}}_{i}({x}_{i},{y}_{i})={C}_{11}{\mathit{\Phi} }_{i}({x}_{i},{y}_{i})+{C}_{12}{\mathit{\Phi} }_{j}({x}_{i},{y}_{i})+{C}_{13}{\mathit{\Phi} }_{k}({x}_{i},{y}_{i})=1\text{，}\\ {\mathit{\Psi} }_{i}({x}_{j},{y}_{j})={C}_{11}{\mathit{\Phi} }_{i}({x}_{j},{y}_{j})+{C}_{12}{\mathit{\Phi} }_{j}({x}_{j},{y}_{j})+{C}_{13}{\mathit{\Phi} }_{k}({x}_{k},{y}_{k})=0\text{，}\\ {\mathit{\Psi} }_{i}({x}_{k},{y}_{k})={C}_{11}{\mathit{\Phi} }_{i}({x}_{k},{y}_{k})+{C}_{12}{\mathit{\Phi} }_{j}({x}_{k},{y}_{k})+{C}_{13}{\mathit{\Phi} }_{k}({x}_{k},{y}_{k})=0。\end{array}\right\} (8) 求解式（8）并结合式（7）得到原粗网格基函数{\mathit{\Psi} _i}的值。{\mathit{\Psi} _j}，{\mathit{\Psi}_k}求解过程类似，这里不再赘述。
1.5 运用MSFEM-T模拟二维地下水流问题
以地下水二维非稳定流为例来说明，控制方程为
\begin{gathered} \frac{\partial}{\partial x}\left(K_x \frac{\partial H}{\partial x}\right)+\frac{\partial}{\partial y}\left(K_y \frac{\partial H}{\partial y}\right)+W=S_{\mathrm{s}} \frac{\partial H}{\partial t} \\ ((x, y) \in \mathit{\Omega}), \end{gathered} (9) 式中，\mathit{\Omega} 为研究区域，W为源汇项，{K_x}，{K_y}为坐标轴x，y方向的渗透系数分量，{S_{\text{s}}}为贮水率。
对式（9）进行伽辽金变换，并离散到每个粗网格单元Δijk上：
{\displaystyle {\iint }_{\mathit{\Omega} }\left[\frac{\partial }{\partial x}\left({K}_{x}\frac{\partial H}{\partial x}\right)+\frac{\partial }{\partial y}\left({K}_{y}\frac{\partial H}{\partial y}\right)+W-{S}_{\text{s}}\frac{\partial H}{\partial t}\right]}\cdot {\mathit{\Psi}}_{i}\text{d}x\text{d}y=0。 (10) 根据多尺度有限元法的原理，在粗网格Δijk上，{\mathit{\Psi} _i}可被线性表示为
H(x,y) = {H_i}{\mathit{\Psi} _i}(x,y) + {H_j}{\mathit{\Psi} _j}(x,y) + {H_k}{\mathit{\Psi} _k}(x,y)\text{，} (11) 式中，{H_i}，{H_j}，{H_k}分别是Δijk顶点i，j，k处的水头值。
结合前面得到的基函数{\mathit{\Psi} _i}，{\mathit{\Psi} _j}和{\mathit{\Psi} _k}，分别结合式（3）、（6）、（11），将式（10）逐层离散到粗、中、细网格上，可以得到式（10）的具体表达式。结合问题的定解条件，可获得Δijk上的单元刚度矩阵和右端自由项，联立所有粗网格可以得到研究区上的总方程组，使用Cholesky分解法进行求解，即可获得研究区的水头值。
2. 算例分析
采用如下缩写形式：LFEM表示采用线性边界条件的有限单元法；LFEM-F表示精细剖分的LFEM；MSFEM-L、MSFEM-O分别表示采用线性边界条件、振荡边界条件的MSFEM；MSFEM-T-L、MSFEM-T-O分别表示采用线性边界条件、振荡边界条件的MSFEM-T；MSFEM-T-OS-n表示采用超样本技术且放大倍数为n的MSFEM-T。
本节应用MSFEM-T模拟了渗透系数连续变化的二维地下水稳定流问题、渗透系数渐变的二维地下水非稳定流问题和渗透系数随机变化的二维地下水稳定流问题，并将模拟的结果与LFEM-F、MSFEM、LFEM进行比较。MSFEM-T的细网格总数与LFEM-F、MSFEM相同。同时，所有数值方法均使用C++编写，方程组解法相同，均未采用并行计算技术，并在相同条件下运行。
2.1 参数连续变化的二维地下水稳定流问题
二维地下水稳定流方程为
\frac{\partial }{\partial x}\left(K\frac{\partial H}{\partial x}\right)+\frac{\partial }{\partial y}\left(K\frac{\partial H}{\partial y}\right)=0\text{ ((}x,y\text{)}\in \mathit{\Omega} \text{) }\text{，} (12) 式中，研究区\mathit{\Omega} =[50 m，150 m]×[50 m，150 m]，渗透系数为连续变化的K(x,y) = {x^2} m/d，此方程具有解析解H = 3{y^2} - {x^2} + 20000 m，第一类边界条件由解析解给出。本案例来源于Ye等[8]的工作，主要考察MSFEM-T处理连续变化的渗透系数的能力。
采用LFEM、LFEM-F、MSFEM-L、MSFEM-O、MSFEM-T-L和MSFEM-T-O。在第一套网格中，LFEM、MSFEM、MSFEM-T将研究区剖分为1800个粗网格单元。MSFEM将每个粗网格直接剖分为36个细网格；MSFEM-T先将每个粗网格剖分为9个中网格，再将每个中网格剖分为4个细网格。LFEM-F直接将研究区剖分为与MSFEM、MSFEM-T相同的64800细网格。
图 4给出了在y=100 m处各数值方法求解的水头绝对误差图。由图可知，LFEM误差最大；MSFEM-T-L的水头精度与MSFEM-L十分相近，且优于LFEM；MSFEM-O、MSFEM-T-O和LFEM-F的误差最小。其中，MSFEM-T-O的水头精度略低于MSFEM-O但相差不大，这是由于MSFEM-T的网格剖分会略微弱化边界条件效果。由图 4各方法的精度关系可以看出，MSFEM-T能够获得与LFEM-F、MSFEM相近的效果，使用基函数的振荡边界条件能够显著提高MSFEM-T和MSFEM的计算精度，结果显示MSFEM-T-O能够较好的处理非均质地下水问题。
![]()
图 4 y=100 m处各数值方法水头绝对误差图Figure 4. Absolute errors of hydraulic heads of different methods, with y=100 m此外，MSFEM-T还可以使用超样本技术降低谐振误差，进一步提高精度，甚至可以超越LFEM-F，这是由于放大的超样本单元能够减少MSFEM-T误差中参数物理尺度和网格尺度的比例项部分（谐振误差）[7]，且原网格边界（现落在超样本单元内部）上的基函数从一维解变成更加精确的二维解，从而提高了计算精度。图 5展示了应用不同倍数超样本技术时的MSFEM-T-O-OS的水头绝对误差，并与LFEM-F结果进行对比，随着放大倍数的增大，MSFEM-T-O-OS的结果更加精确。该结果表明，使用超样本技术可以显著提高MSFEM-T的计算精度，且放大倍数越大，精度越高。当MSFEM-T-O-OS放大倍数为1.14倍时，其结果甚至可以优于LFEM-F，显示MSFEM-T在模拟地下水流问题时可以一定程度上替代LFEM-F。此外，MSFEM-T和MSFEM-T-O-OS的计算本例用时分别为5，6 s，相比于LFEM-F的14301 s，能够大幅降低计算消耗。
![]()
图 5 不同放大倍数时水头绝对误差图Figure 5. Absolute errors of hydraulic heads with different magnifications为了进一步比较MSFEM和MSFEM-T的计算效率，MSFEM-T和MSFEM使用第二套网格将研究区剖分为3200个粗网格单元，每个粗网格单元剖分为81个细网格单元，即共259200个细网格单元。图 6展示了MSFEM和MSFEM-T剖分为不同单元数时的计算用时和计算绝对误差，其中阴影和白色柱状图分别表示不同网格数目时的计算用时，实线和虚线则分别对应不同网格数目时的计算绝对误差。
![]()
图 6 不同单元数时的计算用时和计算绝对误差图Figure 6. CPU time and absolute errors of methods with different element numbers由阴影柱状图和实线可知，当研究区剖分为64800个单元时，MSFEM-T的CPU计算用时约为MSFEM的20%。由白色柱状图和虚线可知，当研究区剖分为更多的259200个单元时，MSFEM-T的计算时间仅为MSFEM的16%。对比两种数目网格下的MSFEM-T和MSFEM的结果，网格的加密会令基函数的构造消耗占总计算消耗的比率上升，可令MSFEM-T与MSFEM之间的CPU用时比率进一步下降。结果表明，相比于MSFEM，MSFEM-T在计算精度相近的情况下，具有更高的计算效率。在计算大尺度问题时，MSFEM-T效率优势将更加明显。
2.2 参数渐变的二维地下水非稳定流问题
该案例是基于实际情况山前冲积平原水流问题的模拟[6]，控制方程为式（9）。其中研究区\mathit{\Omega} =[0 m, 10000 m]×[0 m，10000 m]；渗透系数自左向右从1 m/d缓慢升至251 m/d，即K(x,y) = (40 + x)/40m/d；研究区左、右边界为定水头边界，左右边界水头分别10，0 m，上下为隔水边界；在（5000，5000 m）处有一口流量为1000 m3/d的抽水井；含水层厚度为10 m；贮水率{S_{\text{S}}} = 0.00001 - 0.000009x/1000；初始时刻水头{H_0}(x,y) = 10 - x/1000；时间步长1 d，总时间为5 d；本算例无解析解，故本例将精细剖分为64800个单元的LFEM-F数值解作为“标准解”。
本例中各方法的网格剖分数目均和2.1节相同。图 7展示了各数值方法在y=5000 m时的水头值，各方法均使用第一套网格。由图可知，LFEM计算误差最大，MSFEM-L和MSFEM-T-L其次，采用振荡边界条件的MSFEM-O和MSFEM-T-O曲线最接近LFEM-F的结果，误差最小。使用超样本技术后，当MSFEM-T-O的放大倍数为1.1时，获得了最高的计算精度。值得注意的是，在（5000 m，5000 m）处因为有抽水井的影响，各方法都会出现一定程度的误差，此处精度从高到低依次为MSFEM-T-OS-1.1、MSFEM-O、MSFEM-T-O、MSFEM-L、MSFEM-T-L、LFEM。该结果也表明：MSFEM-T采用振荡边界条件具有比线性边界条件更精确的计算结果；和2.1相同，MSFEM-T的网格剖分会略微弱化边界条件效果（MSFEM-T精度略低于MSFEM但基本相同）；使用超样本技术可以提高MSFEM-T处理抽水井处的计算精度；MSFEM-T和MSFEM二者都具有较好的处理地下水非稳定流的能力。
![]()
图 7 各数值方法在y=5000 m时的水头值Figure 7. Calculated heads of different methods, with y=5000 m图 8比较了各数值方法剖分为不同单元数（64800，259200个网格，剖分同例2.1节）时的计算时间和计算绝对误差图，其中柱状图和折线图代表信息和上例相同。在非稳定流问题中，MSFEM-T和MSFEM仅需在首个时间步构造基函数，而在后面时间步可直接使用此基函数。因此，除基函数计算消耗外，同样粗网格上两种方法的计算时间几乎相同。同时，由于时间步的增加，MSFEM-T和MSFEM处理非稳定流时的总计算用时会略高于稳定流问题。图 9中还对比了MSFEM-L和MSFEM-T-L构造基函数计算用时：MSFEM-T基函数的CPU计算用时明显低于MSFEM，具有更高的基函数构造效率。
![]()
图 8 不同单元数时的计算用时和计算绝对误差图Figure 8. CPU time and absolute errors of methods with different element numbers![]()
图 9 不同单元数时的基函数计算用时图Figure 9. CPU time of basics function of MSFEM-L and MSFEM-T-L with different element numbers2.3 参数随机变化的二维地下水稳定流问题
稳定流控制方程为式（12）。其中研究区\Omega =[0 m，10000 m]×[0 m，10000 m]，原点位于(0 m，0 m)；左右边界为定水头边界，左边界水头为16 m，右边界水头为11 m，上下为隔水边界；无源汇项。本案例使用GSlib序贯高斯模拟器在规则的具有400×400个结点的网格上随机生成渗透系数场[21]，InK相关长度均为100 m，均值为0，变异函数为指数型，方差{\sigma ^2}分别取为2，4。图 10显示了两种不同方差时的InK分布图。
使用LFEM，LFEM-F，MSFEM和MSFEM-T求解本案例，其中MSFEM和MSFEM-T均采用了振荡边界条件。LFEM-F将研究区剖分为45000个单元，LFEM，MSFEM和MSFEM-T均将研究区剖分为1250个粗网格单元，其中MSFEM将每个粗网格直接剖分为36个细网格；MSFEM-T先将每个粗网格剖分为9个中网格，再将每个中网格剖分为4个细网格。因此MSFEM和MSFEM-T共计将研究区剖分为45000个单元，与LFEM-F一致。
图 11（a），（b）分别显示了方差{\sigma ^2}取2，4时各数值方法在y=5000 m处的水头相对误差。从图 11（a）可看出，LFEM误差最大，MSFEM和MSFEM-T其次且非常接近LFEM-F，LFEM-F曲线波动幅度最小，精度最高。图 11（b）显示，当方差增大时，各数值方法的计算误差都相应增大，但依旧保持之前的趋势。该案例进一步说明，MSFEM和MSFEM-T求解随机非均质多孔介质的能力十分相近，均能有效的处理随机非均质的地下水问题。与例2.1节、2.2节类似，在方差{\sigma ^2}为2，4时，MSFEM-T用时仍低于MSFEM，并远低于LFEM-F的计算时间，展示了MSFEM-T具有极高的效率。
![]()
图 11 各数值方法在y=5000 m处的水头相对误差图Figure 11. Relative errors of hydraulic heads of different methods, with y=5000 m3. 结论
本文提出了一种模拟地下水流运动的三重尺度有限元模型（MSFEM-T），该方法通过改进粗网格的剖分方式提高了基函数的计算效率，能够高效、精确求解渗透系数连续变化、渐变和随机变化的二维稳定、非稳定流问题。
（1）在使用相同数量网格情况下，相比于MSFEM和LFEM-F等传统方法，MSFEM-T具有更高的计算效率和相近的计算精度，特别是在大尺度问题中MSFEM-T具有更明显的优势。
（2）与MSFEM结果类似，MSFEM-T对基函数的边界条件敏感，使用振荡边界条件的计算结果比线性边界条件更准确；MSFEM-T的剖分方式会微量降低边界条件的效果，在一般非均质情况下对精度影响较小。
（3）在使用超样本技术后，MSFEM-T的计算精度能够进一步提高，放大倍数增大到一定程度时，精度甚至可以超过精细剖分的有限元法。
-
![]()
图 2 MSFEM粗网格剖分方式
Figure 2. Sketch map of standard subdivision of coarse element of MSFEM
![]()
图 3 MSFEM-T粗网格剖分方式
Figure 3. Sketch map of standard subdivision of coarse element of MSFEM-T
![]()
图 4 y=100 m处各数值方法水头绝对误差图
Figure 4. Absolute errors of hydraulic heads of different methods, with y=100 m
![]()
图 5 不同放大倍数时水头绝对误差图
Figure 5. Absolute errors of hydraulic heads with different magnifications
![]()
图 6 不同单元数时的计算用时和计算绝对误差图
Figure 6. CPU time and absolute errors of methods with different element numbers
![]()
图 7 各数值方法在y=5000 m时的水头值
Figure 7. Calculated heads of different methods, with y=5000 m
![]()
图 8 不同单元数时的计算用时和计算绝对误差图
Figure 8. CPU time and absolute errors of methods with different element numbers
![]()
图 9 不同单元数时的基函数计算用时图
Figure 9. CPU time of basics function of MSFEM-L and MSFEM-T-L with different element numbers
-
[1] 陈雄, 张岩, 王艺伟, 等. 苏北沿海三市三维地下水流数值模拟[J]. 吉林大学学报(地球科学版), 2018, 48(5): 1434–1450. https://www.cnki.com.cn/Article/CJFDTOTAL-CCDZ201805036.htm CHEN Xiong, ZHANG Yan, WANG Yi-wei, et al. Numerical simulation of three dimensional groundwater flow in three coastal cities of north Jiangsu[J]. Journal of Jilin University (Earth Science Edition), 2018, 48(5): 1434–1450. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-CCDZ201805036.htm
[2] YE S J, LUO Y, WU J C, et al. Three-dimensional numerical modeling of land subsidence in Shanghai, China[J]. Hydrogeology Journal, 2016, 24(3): 695–709. doi: 10.1007/s10040-016-1382-2
[3] 李宁, 杨敏, 李国锋. 再论岩土工程有限元方法的应用问题[J]. 岩土力学, 2019, 40(3): 1140–1148, 1157. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201903035.htm LI Ning, YANG Min, LI Guo-feng. Revisiting the application of finite element method in geotechnical engineering[J]. Rock and Soil Mechanics, 2019, 40(3): 1140–1148, 1157. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201903035.htm
[4] 薛禹群, 谢春红. 地下水数值模拟[M]. 北京: 科学出版社, 2007. XUE Yu-qun, XIE Chun-hong. Numerical Simulation for Groundwater[M]. Beijing: Science Press, 2007. (in Chinese)
[5] 周杰, 汪德爟. 有限元水流计算中内存和运行效率初探[J]. 水科学进展, 2004, 15(5): 593–597. doi: 10.3321/j.issn:1001-6791.2004.05.010 ZHOU Jie, WANG De-guan. Exploration on memory requirement and operation efficiency of finite element method in flow calculation[J]. Advances in Water Science, 2004, 15(5): 593–597. (in Chinese) doi: 10.3321/j.issn:1001-6791.2004.05.010
[6] 张丙印, 朱京义, 王昆泰. 非饱和土水气两相渗流有限元数值模型[J]. 岩土工程学报, 2002, 24(6): 701–705. doi: 10.3321/j.issn:1000-4548.2002.06.006 ZHANG Bing-yin, ZHU Jing-yi, WANG Kun-tai. Finite element modeling of two-phase seepage in unsaturated soil[J]. Chinese Journal of Geotechnical Engineering, 2002, 24(6): 701–705. (in Chinese) doi: 10.3321/j.issn:1000-4548.2002.06.006
[7] HOU T Y, WU X H. A multiscale finite element method for elliptic problems in composite materials and porous media[J]. Journal of Computational Physics, 1997, 134(1): 169–189. doi: 10.1006/jcph.1997.5682
[8] YE S J, XUE Y Q, XIE C H. Application of the multiscale finite element method to flow in heterogeneous porous media[J]. Water Resources Research, 2004, 40(9): 337–348.
[9] 范颖, 王磊, 章青. 多尺度有限元法及其应用研究进展[J]. 水利水电科技进展, 2012, 32(3): 90–94. https://www.cnki.com.cn/Article/CJFDTOTAL-SLSD201203025.htm FAN Ying, WANG Lei, ZHANG Qing. Research progress and application of multiscale finite element method[J]. Advances in Science and Technology of Water Resources, 2012, 32(3): 90–94. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLSD201203025.htm
[10] CHEN J, CHUNG E T, HE Z K, et al. Generalized multiscale approximation of mixed finite elements with velocity elimination for subsurface flow[J]. Journal of Computational Physics, 2020, 404: 109133. doi: 10.1016/j.jcp.2019.109133
[11] CHEN Z M, HOU T. A mixed multiscale finite element method for elliptic problems with oscillating coefficients[J]. Mathematics of Computation, 2003, 72(242): 541–576.
[12] JENNY P, LEE S H, TCHELEPI H A. Multi-scale finite-volume method for elliptic problems in subsurface flow simulation[J]. Journal of Computational Physics, 2003, 187(1): 47–67. doi: 10.1016/S0021-9991(03)00075-5
[13] 贺新光, 任理. 求解非均质多孔介质中非饱和水流问题的一种自适应多尺度有限元方法: Ⅰ. 数值格式[J]. 水利学报, 2009, 40(1): 38–45, 51. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200901007.htm HE Xin-guang, REN Li. Adaptive multi-scale finite element method for unsaturated flow in heterogeneous porous media I. Numerical scheme[J]. Journal of Hydraulic Engineering, 2009, 40(1): 38–45, 51. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200901007.htm
[14] 贺新光, 任理. 求解非均质多孔介质中非饱和水流问题的一种自适应多尺度有限元方法: Ⅱ. 数值结果[J]. 水利学报, 2009, 40(2): 138–144. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200902003.htm HE Xin-guang, REN Li. Adaptive multi-scale finite element method for unsaturated flow in heterogeneous porous media Ⅱ. Numerical results[J]. Journal of Hydraulic Engineering, 2009, 40(2): 138–144. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200902003.htm
[15] 黄梦杰, 贺新光. 求解非均质多孔介质中随机水流问题的多尺度有限元降基方法[J]. 水资源与水工程学报, 2019, 30(6): 86–95. https://www.cnki.com.cn/Article/CJFDTOTAL-XBSZ201906014.htm HUANG Meng-jie, HE Xin-guang. Reduced multiscale finite element basis method for solving the stochastic water flow problems in heterogeneous porous media[J]. Journal of Water Resources and Water Engineering, 2019, 30(6): 86–95. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-XBSZ201906014.htm
[16] SHI L S, ZHANG D X, LIN L, et al. A multiscale probabilistic collocation method for subsurface flow in heterogeneous media[J]. Water Resources Research, 2010, 46(11): 386.
[17] 谢一凡, 吴吉春, 王益, 等. 一种模拟节点达西流速的多尺度有限元-有限元模型[J]. 岩土工程学报, 2022, 44(1): 107–114, 202. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202201010.htm XIE Yi-fan, WU Ji-chun, WANG Yi, et al. Multiscale finite element-finite element model for simulating nodal Darcy velocity[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(1): 107–114, 202. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202201010.htm
[18] SHI X Q, WU J C, YE S J, et al. Regional land subsidence simulation in Su-Xi-Chang area and Shanghai City, China[J]. Engineering Geology, 2008, 100(1/2): 27–42.
[19] XIE Y F, WU J C, XUE Y Q, et al. Modified multiscale finite-element method for solving groundwater flow problem in heterogeneous porous media[J]. Journal of Hydrologic Engineering, 2014, 19(8): 04014004.
[20] XIE Y F, WU J C, XUE Y Q, et al. Efficient triple-grid multiscale finite element method for solving groundwater flow problems in heterogeneous porous media[J]. Transport in Porous Media, 2016, 112(2): 361–380.
[21] DEUTSCH C V. GSLIB: Geostatistical Software Library and User's Guide[M]. 2nd ed. Oxford: Oxford University Press, 1998.
-
期刊类型引用(2)
1. 江博. 长源水电站大坝坝肩破坏成因分析研究. 陕西水利. 2025(03): 11-14 . 
百度学术
2. 王顿，裴丽欣，张礼中，樊连杰，卢丽，李习文，刘潇桐，李俊楠，梁林德，白雪冬. 基于文献计量学近20年国内外地下水数值模拟研究进展及展望. 环境工程. 2023(S1): 240-247 . 百度学术
其他类型引用(6)
下载:
计量
- 文章访问数: 156
- HTML全文浏览量: 42
- PDF下载量: 19
- 被引次数: 8
