Multiscale finite element-finite element model for simulating nodal Darcy velocity
-
摘要: 提出了一种能高效模拟节点达西流速并保证其连续性的多尺度有限元-有限元模型(MSFEM-FEM)。该方法先应用多尺度有限元法(MSFEM)框架改进了Yeh的有限元模型的水头模拟部分以提升效率与精度，再将多尺度网格转化为有限元网格，应用Yeh的有限元框架保证流速的连续性。基于多尺度基函数，MSFEM-FEM能够汲取研究区的全局信息并在粗尺度上高效获得精确的水头解。通过将粗尺度网格转换为有限元网格，MSFEM-FEM能够应用Yeh的有限元框架将水头解中的全局信息导入达西流速，提高达西流速的精度并保证其连续性。在获得粗尺度解后，MSFEM-FEM还能应用多尺度基函数对解进行细尺度重构，从而获得研究区内的细尺度水头与流速。数值模拟结果显示MSFEM-FEM能够高效、精确的求解水头，并能够获得连续、精确的达西流速和流量。Abstract: A multi-scale finite element-finite element model (MSFEM-FEM) is proposed, and it can effectively simulate the nodal Darcy velocity and ensure the velocity continuity. The MSFEM-FEM employs the multi-scale finite element method (MSFEM) to replace the head simulation part of the Yeh's finite element model, thus to improve the efficiency and accuracy. Then, the MSFEM-FEM transforms the multi-scale grid to the finite element one, thus, it can directly apply the Yeh's finite element model to obtain continuous Darcy velocity. Based on the multi-scale basis function, the MSFEM-FEM can extract the global information of the study area which allows it to obtain the accurate head solution efficiently on the coarse scale. By transforming the coarse-scale grid into the finite element grid, the MSFEM-FEM can directly employ the Yeh's finite element model to import the global information from the head solution into the Darcy velocity, which can also improve the accuracy of Darcy velocity and ensure the velocity continuity. In addition, the MSFEM-FEM can apply multi-scale basis function to reconstruct the solutions, so as to obtain the fine-scale head and velocity solutions in the study area. The simulated results of two-dimensional groundwater problems show that the MSFEM-FEM can efficiently and accurately solve the head, Darcy velocity and flux, which outperforms the MSFEM and the Yeh's finite element model.
-
0. 引言
在地下水资源评价和地下水污染迁移等研究中，水文工作者常需要将对流-弥散方程与描述地下水运动的水流方程进行耦合，从而对地下水中的溶质运移状态进行准确模拟。地下水的达西流速不仅能够显示地下水流的运动情况，还反映了对流作用对溶质运移过程的影响，是联系两个方程的纽带。因此，如何高效、精确地模拟达西流速与流量对地下水溶质运移的研究具有重要意义。然而，有限元类方法的不能直接获得连续的节点水头导数，这导致其所计算的达西流速不连续且精度较低[1-3]。
Yeh在1981年建立了基于有限元法的达西流速模型[2]，能够保证节点达西流速的连续性，具有较高的计算精度。该模型先使用有限单元法求解水头以获得全局信息和水头分布，再用有限元法直接在研究区上求解达西方程，从而得到连续的达西流速场，并能保证解的质量平衡[2, 4]。和混合有限元类方法[5-6]相比，Yeh的有限元模型更简单直观，适于实际应用[3]。和同类的Batu的双重网格法[7]和Zhang的三次样条法[3]比，Yeh的有限元模型的精度略高，但应用限制更小[3]。然而，在进行大尺度地下水问题的模拟时，有限元法需要精细剖分保证解的精度[1, 8-9]，故Yeh的模型需要大量计算时间获得水头分布。实际上，该模型的水头模拟部分主要是用于获得研究区的介质信息和模拟水头分布，其流速模拟部分才是保证达西流速连续性的关键。
为了解决有限元法计算效率较低的问题，Hou等[9]在1997年提出了多尺度有限单元法(MSFEM)。MSFEM的核心是通过在粗网格单元求解退化椭圆的方程构造多尺度基函数，并通过多尺度基函数获取细尺度信息，从而节约大量的计算成本[1, 9-10]。在计算大尺度地下水问题时，许多工作已经证明了MSFEM比有限元法更加高效[9-14]，并且MSFEM已经被成功应用于中国地面沉降的模拟工作中[15]。然而，和有限元法相同，MSFEM也无法直接获得连续的水头一阶导数，故其也不能获得连续、精确的达西流速。谢一凡等[16]、赵文凤等[17]应用MSFEM改进了多种经典达西流速算法，取得了较好的结果，发现水头精度对于达西流速的精度有重要影响[16-17]。
本文提出了一种模拟连续节点达西流速的多尺度有限元-有限元模型(MSFEM-FEM)。MSFEM-FEM综合了两种方法各自的优势，实现了水头、连续达西流速和流量的高效计算。MSFEM-FEM应用MSFEM改进了Yeh的有限元模型的水头模拟部分，能够更加精确、高效地获取全局信息和水头分布。通过将多尺网格转换为有限元网格，MSFEM-FEM运用Yeh的有限元框架保证了达西流速与流量的连续性。同时，MSFEM-FEM还能运用多尺度基函数在粗网格单元内重构细尺度的水头和达西流速，从而以极低的计算成本获得与精细剖分的Yeh的有限元模型相同数目的解。
在数值试验中，将MSFEM-FEM和Yeh的有限元模型进行了比较，证明了它可以更高效地模拟水头、达西流速和流量。结果显示MSFEM-FEM获得了与精细剖分的Yeh有限元模型同阶的水头精度和达西流速精度。在其中一个算例中，MSFEM-FEM计算时间甚至不到精细剖分的Yeh有限元法的0.07%。本文还将MSFEM-FEM与基于MSFEM的传统流量计算方法进行了比较，结果显示MSFEM-FEM所计算的流量精度更高，并能保证通过含水层截面的流入、流出量相等。
1. 原理和算法
MSFEM-FEM包含地下水水头模拟和达西流速模拟两个部分，分别由MSFEM和Yeh的有限元模型实现，且各部分的剖分方法不同。在获得粗尺度的水头和达西流速后，MSFEM-FEM可以应用多尺度基函数在水头模拟部分的多尺度网格上实现解的细尺度重构。
1.1 构造网格剖分
MSFEM-FEM对研究区有两种不同的剖分：计算水头的多尺度剖分以及计算达西流速的有限元剖分。设研究区为矩形区域，水平、垂直坐标轴分别为x，y轴，具体网格剖分构造方式如下：
在水头模拟部分，MSFEM-FEM先将研究区Ω的边界等分为若干份，以水平、垂直线段连接等分点将研究区剖分为矩形粗网格单元，完成粗尺度剖分，如，图 1中实线将研究区剖分为100个矩形粗网格单元，${\square _{ijkl}}$为一示例粗网格单元。然后，MSFEM-FEM将每个粗网格单元边界等分，以水平、垂直、对角线将粗网格单元剖分为细网格单元，完成细尺度剖分，将示例粗网格${\square _{ijkl}}$剖分为18个细网格单元，如图 2所示。
在达西流速模拟部分，MSFEM-FEM基于水头模拟部分网格的粗尺度节点将粗尺度网格转换为有限元网格，即连接粗网格单元的对角线以获得有限元单元。如，图 1中虚线将每个粗网格单元剖分为2个三角形有限元单元，即研究区总共被剖分为200个有限元单元。
1.2 构造基函数
设示例矩形粗网格单元${\square _{ijkl}}$(图 2)上的顶点按逆时针顺序为i，j，k，l，对应多尺度基函数为${\varPsi _i}$，${\varPsi _j}$，${\varPsi _k}$，${\varPsi _l}$。以构造基函数${\varPsi _i}$的过程为例，考虑如下简化的椭圆方程[9]：
$$ {\text{ - }}\nabla \cdot ({\boldsymbol K}\nabla {\varPsi _i}) = 0{\text{ ((}}x{\text{, }}y{\text{)}} \in {\square _{ijkl}}) \text{，} $$ (1) 式中，${\boldsymbol K}$为渗透系数。在指定多尺度基函数${\varPsi _i}$的边界条件后，方程(1)是适定的。
对式(1)实行伽辽金变换，可得
$$ {J_{{M_\tau }}} = \iint_{{\square _{ijkl}}} {({\boldsymbol K}\nabla {\varPsi _i}}) \cdot \nabla {\varphi _{{M_\tau }}}{\text{d}}x{\text{d}}y = 0{\text{ (}}\tau = 1, 2, \cdot \cdot \cdot , p) \text{，} $$ (2) 式中，${\varphi _{{M_\tau }}}$为在${M_\tau }$点的线性基函数，p为${\square _{ijkl}}$的内点数。
设细网格单元顶点按逆时针顺序为a，b，c，根据多尺度有限元理论[9]，${\varPsi _i}$可被线性表示为
$$ {\varPsi _i}(x, y) = {\varPsi _i}(a){\varphi _a} + {\varPsi _i}(b){\varphi _b} + {\varPsi _i}(c){\varphi _c}\text{，} $$ (3) 式中，${\varphi _a}$，${\varphi _b}$，${\varphi _c}$为与细网格顶点相关的线性基函数。
将式(3)代入式(2)，可以得到一个关于${\varPsi _i}$的内点值的方程组，该方程组的系数是对称正定的。求解这个方程组，可以得到${\varPsi _i}$在粗网格单元内部所有节点上的值。${\varPsi _j}$，${\varPsi _k}$，${\varPsi _l}$的构造过程与${\varPsi _i}$类似。
1.3 运用MSFEM-FEM模拟地下水流和达西流速
(1) 运用MSFEM-FEM模拟地下水头分布
以研究区Ω上的地下水二维稳定流问题为例：
$$ - \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， $$ (4) 式中，${K_x}$，${K_y}$为渗透系数，H为水头，W为源汇项。在确定方程的定解条件后，该问题适定。
设研究区Ω被剖分为$\gamma $个粗网格单元，在式(4)两边乘以多尺度基函数${\varPsi _i}$，由伽辽金变换可得
$$\sum\limits_1^\gamma {\iint {\left[ {\left( {{K_x}\frac{{\partial H}}{{\partial x}}} \right)\frac{{\partial {\varPsi _i}}}{{\partial x}} + \left( {{K_y}\frac{{\partial H}}{{\partial {\text{y}}}}} \right)\frac{{\partial {\varPsi _i}}}{{\partial y}}} \right]}} {\text{d}}x{\text{d}}y \\ ={\displaystyle {\sum }_{1}^{\gamma }{\displaystyle \iint W{\varPsi }_{i}}}\text{d}x\text{d}y\text{ }。 $$ (5) 根据Hou等[9]，在粗网格单元${\square _{ijkl}}$上，水头可以被多尺度基函数线性表示：
$$ H = {H_i}{\varPsi _i} + {H_j}{\varPsi _j} + {H_k}{\varPsi _k} + {H_l}{\varPsi _l}\text{，} $$ (6) 式中，${H_i}$，${H_j}$，${H_k}$，${H_l}$分别为粗网格单元${\square _{ijkl}}$顶点i，j，k，l处的水头值。
将式(6)代入式(5)在${\square _{ijkl}}$上的分量，可以得到一个线性方程：
$$ {B_{ii}}{H_i} + {B_{ij}}{H_j} + {B_{ik}}{H_k} + {B_{il}}{H_l} = {f_i}\text{，} $$ (7) 式中，系数$ {B}_{i\eta }\text{(}\eta =i，j，k，l) $为常数，右端项${f_i} = $$\int_{{\square _{ijkl}}} {W{\varPsi _i}{\text{d}}x{\text{d}}y} $，通过式(3)离散式(7)到细网格上可以得到${B_{i\eta }}$和${f_i}$的具体形式。类似的，可以得到关于${\varPsi _j}$，${\varPsi _k}$，${\varPsi _l}$在${\square _{ijkl}}$上的方程，联立可以获得${\square _{ijkl}}$的单元方程组。联立所有粗网格单元的单元方程组，结合地下水问题的定解条件，可以得到关于粗尺度水头的方程组。
(2) 运用MSFEM-FEM模拟连续达西流速
在研究区Ω上考虑达西方程：
$$ {v_x} = - {K_x}\frac{{\partial H}}{{\partial x}}\text{，} $$ (8) 式中，${v_x}$为坐标轴x方向的达西流速。
设Ω被剖分为${\gamma _2}$(${\gamma _2}$=$2\gamma $)个三角形有限元单元${\Delta _{ijk}}$(如图 1)，共有${N_n}$个节点。在式(8)两端分别乘以基函数${\varphi _I}(I = 1, 2, \cdots , {N_n})$，在整个研究区上积分后，离散到各个有限元单元△ijk上：
$$ \sum\nolimits_1^{{\gamma _2}} {\iint {{V_x}{\varphi _I}}{\text{d}}x{\text{d}}y = - \sum\nolimits_1^{{\gamma _2}} {\iint {{K_x}\frac{{\partial H}}{{\partial x}}{\varphi _{\text{I}}}}} } {\text{d}}x{\text{d}}y。 $$ (9) 基于Yeh的有限元模型[2]，在${\Delta _{ijk}}$上${v_x}$和水头H可以被与其顶点相关的线性基函数分别表示为
$$ {v_x}(x, y) = {v_x}({x_i}, {y_i}){\varphi _i} + {v_x}({x_j}, {y_j}){\varphi _j} + {v_x}({x_k}, {y_k}){\varphi _k}\text{，} $$ (10) $$ H = {H_i}{\varphi _i} + {H_j}{\varphi _j} + {H_k}{\varphi _k}。 $$ (11) 将式(10)，(11)代入式(9)可以得到关于Vx的方程组，求解可以得到${v_x}$在研究区Ω的有限元网格节点(即多尺度网格的粗尺度节点)上的达西流速。
限于篇幅，本节仅给出了MSFEM-FEM变换其粗尺度网格模拟整个研究区达西流速时的方法。若只需获得局部区域的达西流速场，MSFEM-FEM可以使用该局部区域水头网格来变换有限元网格，并在该局部区域基于上述过程进行达西流速模拟。此外，若对精度有更高要求，MSFEM-FEM也可用式(6)获得细尺度水头，再将其细尺度网格变换为精细有限元网格直接模拟细尺度达西流速，但计算效率较低。
(3) 解的细尺度重构
根据多尺度有限元理论[9]，细尺度水头可以在粗网格上被粗网格顶点的水头值和多尺度基函数表示为式(6)。类似的，细尺度达西流速在粗网格内可以被多尺度基函数和粗网格${\square _{ijkl}}$顶点处的达西流速表示为
$$ {v_x}(x, y) = {v_x}({x_i}, {{\text{y}}_i}){\varPsi _i} + {v_x}({x_j}, {y_j}){\varPsi _j} + {v_x}({x_k}, {y_k}){\varPsi _k} + {v_x}({x_l}, {y_l}){\varPsi _l}。 $$ (12) 由于${\square _{ijkl}}$顶点(即粗尺度节点)处的达西流速已由1.3(2)节获得，MSFEM-FEM可以应用式(12)直接获得细尺度达西流速。
(4) 流量模拟过程
首先，简单介绍传统有限元法模拟流量的方式。设$\varGamma $为研究区内部的一垂直截面，可以由若干单元的边ac组成。和传统有限单元法类似，MSFEM需要通过水头计算流量[1]。根据薛禹群等[1]，设通过某一单元边界ac的流量${Q_{ac}}$为
$$ {Q_{ac}} = - K{l_{ac}}M\frac{{\partial H}}{{\partial n}} = - K{l_{ac}}M\left[ {\frac{{\partial H}}{{\partial x}}\cos (n, x) + \frac{{\partial H}}{{\partial y}}\cos (n, y)} \right]\text{，} $$ (13) 式中，M为含水层厚度，n为ac的外法线单位向量，lac为截面ac的长度。
则MSFEM在粗网格${\square _{ijkl}}$的某一细网格单元上(图 2)，结合式(6)，通过截面ac的流量也可以由式(13)得到。然而，有限元法不能获得连续的水头一阶导数，故其通过$\varGamma $左侧单元和右侧单元以式(13)所计算的${Q_{ac}}$不相同，从而造成其所获的$\varGamma $上流入、流出量不等。类似的，MSFEM也会在$\varGamma $上得到两个不同的流量。
如前文所述，MSFEM-FEM能够保证节点达西流速的连续性，即达西流速在每个节点上的值是唯一的，故MSFEM-FEM可以直接用法向节点达西流速乘以截面长度和含水层厚度获得通过该截面的流量，设$\varGamma $上节点为i=1，2，…，n，MSFEM-FEM在$\varGamma $上流量公式为
$$ {Q_\varGamma } = M\sum\nolimits_{i = 1}^{n - 1} {\frac{{{v_x}(i) + {v_x}(i + 1)}}{2}} {l_{i, i + 1}}\text{，} $$ (14) 式中，${l_{i, i + 1}}$为点i，i+1之间的距离。
2. 数值试验
采用如下缩写形式：${v_x}$为x方向上的达西流速。N为研究区每边被等分的份数，AS为解析解，MSFEM-FEM表示本文的多尺度有限元-有限元模型，FEM表示有限元法，FEM-F表示精细剖分的有限元法，Yeh表示Yeh的有限元模型，Yeh-F表示精细剖分的Yeh的有限元模型，AS-Yeh表示Yeh有限元模型利用解析解水头来求解达西流速。
本节应用MSFEM-FEM模拟不同条件下的地下水水头和达西流速，采用多尺度基函数的振荡边界条件[9-10]和矩形粗网格单元。许多工作已经证明Yeh的有限元模型在同类方法中精度最高且质量平衡[2-4]。类似于FEM-F之于水头，具有精细剖分的Yeh-F能够以非常高的精度获得达西流速，可与实际问题较好地吻合。因此，本文将MSFEM和Yeh、Yeh-F进行了比较，并应用具有精细剖分的Yeh-F作为部分算例的参照解。Yeh/Yeh-F的有限元模型采用FEM/FEM-F模拟水头，使用了三角形网格单元。为了比较计算消耗和精度，Yeh对研究区边界的等分数N和MSFEM-FEM粗网格相同，即Yeh的节点数等于MSFEM-FEM的粗尺度节点数，Yeh的单元数是MSFEM-FEM的粗网格数的两倍。Yeh-F对研究区边界的等分数N则和MSFEM-FEM的细网格相同，即Yeh-F的节点数等于MSFEM-FEM的细尺度节点数，Yeh-F的单元数等于MSFEM-FEM的细网格数。此外，在流量的比较过程中，MSFEM直接使用传统流量求解方法式(13)求解流量，故其所获得的流量在含水层截面上不连续。将MSFEM-FEM与MSFEM求解出的流量和解析解进行了比较，证明了MSFEM-FEM所获流量的连续性和精度。
2.1 参数连续的二维地下水稳定流问题
水头和达西流速的控制方程分别为式(4)，(8)，${K_x} = {K_y} = {x^2}$，渗透系数的最大值与最小值的比值为16，研究区Ω为[5 m，20 m]×[5 m，20 m]，边界条件由解析解$H = {x^2} - 3{y^2}$给出，含水层厚度为1 m。
采用AS-Yeh、MSFEM-FEM和Yeh模拟粗尺度达西流速。MSFEM-FEM和Yeh将Ω的边界分为10、20、30等份。故Yeh将Ω分为200，800，1800个三角形单元；MSFEM-FEM将Ω分为100，400，900个矩形粗网格单元，再把每个粗网格单元分为50个细网格单元。
表 1为各方法求解出的水头的平均相对误差，由于AS-Yeh使用的是解析解，其水头误差为0。从表 1可以看出MSFEM-FEM所获的水头始终比Yeh使用FEM获得的水头精度更高。
表 1 例2.1中各方法的水头相对误差Table 1. Relative errors of head calculated by numerical methods in example 2.1(%) N AS-Yeh MSFEM-FEM Yeh 10 0 0.0067 0.161 20 0 0.0032 0.079 30 0 0.0010 0.024 图 3为AS-Yeh、MSFEM-FEM和Yeh求解出的粗尺度达西流速${v_x}$的相对误差，第一到三行分别对应于N=10，20，30。从每行图片的垂直刻度中可以看出，各方法的误差随着N的增加而减少。关注各个子图上部的高度，可以发现MSFEM-FEM和AS-Yeh误差均低于Yeh。同时，根据第一、二列子图的形状，可知MSFEM-FEM和AS-Yeh的精度几乎一致，说明MSFEM-FEM已经基本去除了水头引起的误差。当N=10，20，30时，MSFEM-FEM的CPU时间分别为少于1，1，6 s，Yeh的CPU时间分别为少于1 s，少于1，1 s，AS-Yeh无需计算水头，时间仅为Yeh的一半。然而，Yeh仅仅计算了研究区上(N+1) 2个粗尺度节点上的达西流速，而MSFEM-FEM不仅仅计算了粗尺度节点的达西流速${v_x}$，还能够应用多尺度基函数获得研究区(5N+1)2个细尺度节点上的达西流速${v_x}$。
![]()
图 3 例2.1中各方法模拟的粗尺度达西流速${v_x}$的平均相对误差Figure 3. Relative errors of ${v_x}$ calculated by numerical methods in example 2.1本例将MSFEM-FEM所计算的细尺度达西流速${v_x}$与Yeh-F所计算的达西流速进行了对比。为了和MSFEM-FEM的N=30的结果对比，Yeh-F将研究区$ \varOmega $每边分为150份，即45000个三角形单元，总结点数22801等于MSFEM-FEM所有的细尺度节点数。表 2为MSFEM-FEM和Yeh-F的数值结果对比，可以看出MSFEM-FEM获得了与Yeh-F一样的水头精度，且两种方法的达西流速平均相对误差均低于0.35%。MSFEM-FEM的达西流速的精度略低于Yeh-F，这是由于Yeh-F直接在精细网格上求解，能够比应用多尺度基函数插值重构达西流速的MSFEM-FEM获取更多的细尺度信息。然而，MSFEM-FEM的达西流速平均相对误差仅为0.33%，对于大部分实际模拟问题已经足够。另一方面，MSFEM-FEM的CPU时间不到Yeh-F的0.07%，具有极高的计算效率。
表 2 例2.1中MSFEM-FEM与Yeh-F的数值结果对比Table 2. Numerical results of MSFEM-FEM and Yeh-F in example 2.1方法 水头平均相对误差/% 达西流速平均相对误差/% CPU时间/s MSFEM-FEM 0.001 0.33 6 Yeh-F 0.001 0.02 8530 MSFEM-FEM能够获得连续的达西流速场，故其所计算的通过任意截面的流量均是唯一的。以N=20的结果为例，本例将MSFEM-FEM所计算的流量与MSFEM运用传统算法[1]即式(13)所获得的流量进行了对比。图 4为解析解以及MSFEM-FEM和MSFEM所计算的通过研究区所有垂直截面x=5.5，6.0，…，19.5 m的数值解流量。从图 4中可以看出，MSFEM-FEM和解析解吻合得很好，且通过每个垂直截面的流量值均唯一；而MSFEM通过截面左、右侧单元所计算的流入、流出量不相等，且具有10%左右的误差。
![]()
图 4 例2.1中解析解以及MSFEM-FEM和MSFEM所计算的通过Ω各垂直截面的流量Figure 4. Analytical flux and flux calculated by MSFEM-FEM and MSFEM at vertical sections in Ω in example 2.12.2 基于冲积平原的参数渐变的地下水二维稳定流模型
本例是基于实际山前冲积平原条件的水流模型的模拟[1, 10]。水头和达西流速的控制方程分别为式(4)和式(8)，研究区Ω为[0, 10 km]×[0, 10 km]，${K_x} = $${K_y} = \frac{{40 + x}}{{40}}$ ${\text{m/d}}$，从左至右缓慢增加，是典型的山前冲积平原介质的变化特征，渗透系数的最大值与最小值的比值为251。研究区上下为隔水边界，左右为一类边界，左边界水头为10 m，右边界水头为0 m，源汇项为0。含水层厚度为10 m。由于此例没有解析解，本文将Yeh-F的数值解作为参照解进行对照。
Yeh-F将研究区每边剖分为250份，共125000个三角形单元。Yeh和MSFEM-FEM将研究区每边剖分为25份，则Yeh将研究区剖分为1250个三角形单元，MSFEM-FEM将研究区剖分为625个矩形粗网格单元。MSFEM-FEM将每个粗网格剖分为200个三角形细网格单元以获得与Yeh-F相同的单元数和细尺度节点数。
图 5给出了截面y=5200 m处粗尺度节点上的参照解以及MSFEM-FEM和Yeh所计算的数值解水头。可以看出，MSFEM-FEM的解与参照解几乎重合，精度高于Yeh。这是因为MSFEM-FEM通过多尺度基函数抓住了全局信息，并反映在水头中。
![]()
图 5 例2.2中在截面y=5200 m处的参照解以及MSFEM-FEM和Yeh的数值解水头Figure 5. Reference solution and head calculated by MSFEM-FEM and Yeh at section y=5200 m in example 2.2图 6展示了截面y=5200 m处MSFEM-FEM和Yeh的粗尺度达西流速${v_x}$的相对误差，MSFEM-FEM依然取得了更精确的结果。MSFEM-FEM和Yeh均在第一个节点(x=400 m处)有较大的误差，且MSFEM-FEM的误差更大，这是由参照解和数值解的尺度不同而引起的。研究区左侧水头为10 m，在左侧快速下降到5.7 m左右，然后平稳变化至0 m；MSFEM-FEM和Yeh模拟达西流速的网格尺度为400 m，在第一个节点处([0 m, 400 m]区间)分别将水头从10 m变化到5.71，6.5 m；参照解的网格尺度为40 m，可以用更细的尺度进行水头变化，由于水头变化是先快后慢，其在[360 m, 400 m]区间水头从5.87 m变化到5.57 m。因此，由于网格尺度的原因，参照解在x=400 m处的水力梯度小于其他方法，MSFEM-FEM的水头解更精确，导致其达西渗流误差集中在第一个节点，而Yeh则将该误差平均到了整个水平尺度，故其每个节点均具有较大的误差。此外，相比例2.1，本例中MSFEM- FEM与Yeh之间的达西流速精度的差距更大，表明MSFEM-FEM能更好适应渗透系数变化更大的情况。
![]()
图 6 例2.2中在截面y=5200 m处的MSFEM-FEM和Yeh的数值解达西渗透流速${v_x}$的相对误差Figure 6. ${v_x}$ calculated by MSFEM-FEM and Yeh at section y=5200 m in example 2.2本例中MSFEM-FEM使用了6 s的CPU时间获得了624个粗尺度水头、62499个细尺度水头、676个粗尺度达西流速${v_x}$以及63001个细尺度达西流速，而Yeh虽然只消耗1 s的CPU时间，但其仅计算了624个粗尺度水头和676个粗尺度。Yeh-F获得了与MSFEM-FEM相同数目的解，但其需要253314 s计算水头、250780 s计算达西流速${v_x}$。然而，图 5，6均显示MSFEM-FEM与参照解(Yeh-F)的精度差距非常小，这证明了MSFEM-FEM在模拟大尺度问题时具有极高的计算效率。
2.3 参数高度振荡的二维地下水稳定流问题
本例的水头和达西流速的控制方程分别为式(4)和式(8)，研究区Ω为[0, 1[L]] ×[0, 1[L]]，${K_x} = {K_y} = $$\frac{1}{{2 + 1.995\sin \left( {\frac{{7{\text{π }}(x + y)}}{4}} \right)}}$[LT-1]，渗透系数的最大值与最小值的比值为800，解析解为$H{\text{ = }}xy(1 - x)(1 - y)$[L]，源汇项$W$以及研究区四边的第一类边界条件由解析解给出，含水层厚度为1[L]。
Yeh和MSFEM-FEM将研究区每边剖分为30份(N=30)，则Yeh将研究区剖分为1800个三角形单元；MSFEM-FEM将研究区剖分为900个矩形粗网格单元，每个粗网格剖分为18个三角形细网格单元。Yeh-F将研究区每边剖分为90份(N=90)，即共16200个单元，获得与MSFEM-FEM相同的单元数和细尺度节点数。
图 7为各方法在截面y=0.6 [L]处的水头绝对误差，可以看出MSFEM-FEM和Yeh-F的解十分接近，具有很高的计算精度；而Yeh则受到了高度振荡渗透系数的影响，产生了很大的误差。这一结果证明了MSFEM-FEM和Yeh-F在模拟水头时具有较强的处理非均质系数的能力。图 8为各方法在截面y=0.6 [L]处的达西流速${v_x}$的值，可以看出MSFEM-FEM、Yeh-F、Yeh基本继承了它们在图 7中水头模拟部分的精度关系，这一结果说明了在渗透系数具有较强振荡时水头误差是达西流速误差的主要来源。
![]()
图 7 例2.3中在截面y=0.6 [L]处的Yeh-F、MSFEM-FEM和Yeh的水头绝对误差Figure 7. Absolute errors of head calculated by Yeh-F, MSFEM-FEM and Yeh at section y=0.6 [L] in example 2.3![]()
图 8 例2.3中在截面y=0.6 [L]处的解析解以及Yeh-F、MSFEM-FEM和Yeh的数值解${v_x}$Figure 8. Analytical solution and ${v_x}$calculated by Yeh-F, MSFEM- FEM and Yeh at section y=0.6 [L] in example 2.3同时，MSFEM-FEM获得了与Yeh-F精度相近数量相同的解，说明MSFEM-FEM在模拟水头时也具有较强的处理非均质渗透系数的能力，但其的CPU时间仅为5 s，不到Yeh-F(475 s)的1.1%。Yeh虽然仅需要1 s的CPU时间，但其计算精度和解的数量均与其他两个方法差距巨大。
3. 结论
本文提出了一种模拟节点达西流速的多尺度有限元——有限元模型(MSFEM-FEM)。该方法能够通过MSFEM框架提升计算效率，并通过多尺度基函数抓住研究区的全局信息以反映到水头解中，再通过Yeh的有限元模型将全局信息传递到达西流速中，从而提高水头和达西流速的精度。同时，MSFEM-FEM能够保证达西流速的连续性，故其所计算的流量在研究区内截面上各节点的值是唯一的。应用多种条件下的地下水数值问题验证了MSFEM-FEM的有效性，得到以下4点结论。
(1) MSFEM-FEM能够精确地获得水头和达西流速，解的精度接近于Yeh-F并高于Yeh。
(2) 在本文比较的多种方法中，MSFEM-FEM具有最高的计算效率，能够节约大量的计算消耗。
(3) MSFEM-FEM能够保证达西流速和流量的连续性，精度高于MSFEM应用传统流量求解方式。
(4) MSFEM-FEM适用于多种不同变化形式的渗透系数，且具有很强的处理介质非均质性的能力；在渗透系数高度振荡时，MSFEM-FEM的精度远超Yeh。
总而言之，MSFEM-FEM是一种高效、精确的水头和达西流速的综合模拟方法，能够有效处理地下水问题的介质非均质性，具有广阔的发展前景。
-
![]()
图 3 例2.1中各方法模拟的粗尺度达西流速的平均相对误差
Figure 3. Relative errors of calculated by numerical methods in example 2.1
![]()
图 4 例2.1中解析解以及MSFEM-FEM和MSFEM所计算的通过Ω各垂直截面的流量
Figure 4. Analytical flux and flux calculated by MSFEM-FEM and MSFEM at vertical sections in Ω in example 2.1
![]()
图 5 例2.2中在截面y=5200 m处的参照解以及MSFEM-FEM和Yeh的数值解水头
Figure 5. Reference solution and head calculated by MSFEM-FEM and Yeh at section y=5200 m in example 2.2
![]()
图 6 例2.2中在截面y=5200 m处的MSFEM-FEM和Yeh的数值解达西渗透流速的相对误差
Figure 6. calculated by MSFEM-FEM and Yeh at section y=5200 m in example 2.2
![]()
图 7 例2.3中在截面y=0.6 [L]处的Yeh-F、MSFEM-FEM和Yeh的水头绝对误差
Figure 7. Absolute errors of head calculated by Yeh-F, MSFEM-FEM and Yeh at section y=0.6 [L] in example 2.3
![]()
图 8 例2.3中在截面y=0.6 [L]处的解析解以及Yeh-F、MSFEM-FEM和Yeh的数值解
Figure 8. Analytical solution and calculated by Yeh-F, MSFEM- FEM and Yeh at section y=0.6 [L] in example 2.3
表 1 例2.1中各方法的水头相对误差
Table 1 Relative errors of head calculated by numerical methods in example 2.1
(%) N AS-Yeh MSFEM-FEM Yeh 10 0 0.0067 0.161 20 0 0.0032 0.079 30 0 0.0010 0.024 
下载: 导出CSV
表 2 例2.1中MSFEM-FEM与Yeh-F的数值结果对比
Table 2 Numerical results of MSFEM-FEM and Yeh-F in example 2.1
方法 水头平均相对误差/% 达西流速平均相对误差/% CPU时间/s MSFEM-FEM 0.001 0.33 6 Yeh-F 0.001 0.02 8530
下载: 导出CSV
-
[1] 薛禹群, 谢春红. 地下水数值模拟[M]. 北京: 科学出版社, 2007: 175–178. XUE Yu-qun, XIE Chun-hong. Numerical simulation for groundwater[M]. Beijing: Science Press, 2007: 175–178. (in Chinese)
[2] YEH G T. On the computation of Darcian velocity and mass balance in the finite element modeling of groundwater flow [J]. Water Resources Research, 1981, 17(5): 1529–1534. doi: 10.1029/WR017i005p01529
[3] ZHANG Z H, XUE Y Q, WU J C. A cubic-spline technique to calculate nodal Darcian velocities in aquifers[J]. Water Resources Research, 1994, 30(4): 975–98. doi: 10.1029/93WR03416
[4] PARK C H, ARAL M M. Sensitivity of the solution of the Elder problem to density, velocity and numerical perturbations[J]. Journal of Contaminant Hydrology, 2007, 92(1/2): 33–49.
[5] SRINIVAS C, RAMASWAMY B, WHEELER M F. Mixed finite element methods for flow through unsaturated porous media[M]// RUSSELL T F, EWING R E, BREBBIA C A, et al, ed. Computational Methods in Water Resources: Numerical Methods in Water Resources. Southampton: Elsevier Applied Science, 1992: 239–246.
[6] D'ANGELO C, SCOTTI A. A mixed finite element method for Darcy flow in fractured porous media with non-matching grids[J]. ESAIM: Mathematical Modelling and Numerical Analysis, 2012, 46(2): 465–489. doi: 10.1051/m2an/2011148
[7] BATU V. A finite element dual mesh method to calculate Nodal Darcy velocities in nonhomogeneous and anisotropic aquifers[J]. Water Resources Research, 1984, 20(11): 1705–1717. doi: 10.1029/WR020i011p01705
[8] 周杰, 汪德爟. 有限元水流计算中内存和运行效率初探[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
[9] 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
[10] 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): W09202.
[11] 贺新光, 任理. 求解非均质多孔介质中非饱和水流问题的一种自适应多尺度有限元方法: Ⅰ数值格式[J]. 水利学报, 2009, 40(1): 38–45, 51. doi: 10.3321/j.issn:0559-9350.2009.01.006 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) doi: 10.3321/j.issn:0559-9350.2009.01.006
[13] 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
[14] HE X G, LI Q Q, JIANG L J. A reduced generalized multiscale basis method for parametrized groundwater flow problems in heterogeneous porous media[J]. Water Resources Research, 2019, 55(3): 2390–2406. doi: 10.1029/2018WR023954
[15] 于军, 吴吉春, 叶淑君, 等. 苏锡常地区非线性地面沉降耦合模型研究[J]. 水文地质工程地质, 2007, 34(5): 11–16. doi: 10.3969/j.issn.1000-3665.2007.05.004 YU Jun, WU Ji-chun, YE Shu-jun, et al. Research on nonlinear coupled modeling of land subsidence in Suzhou, Wuxi and Changzhou areas, China[J]. Hydrogeology & Engineering Geology, 2007, 34(5): 11–16. (in Chinese) doi: 10.3969/j.issn.1000-3665.2007.05.004
[16] 谢一凡, 吴吉春, 薛禹群, 等. 一种模拟节点达西渗透流速的三次样条多尺度有限单元法[J]. 岩土工程学报, 2015, 37(9): 1727–1732. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201509030.htm XIE Yi-fan, WU Ji-chun, XUE Yu-qun, et al. Cubic-spline multiscale finite element method for simulation of nodal Darcy velocities in aquifers[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(9): 1727–1732. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201509030.htm
[17] 赵文凤, 谢一凡, 吴吉春. 一种模拟节点达西渗透流速的双重网格多尺度有限单元法[J]. 岩土工程学报, 2020, 42(8): 1474–1481. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202008018.htm ZHAO Wen-feng, XIE Yi-fan, WU Ji-chun. A dual-mesh multiscale finite element method for simulating nodal Darcy velocities in aquifers[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(8): 1474–1481. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202008018.htm
-
期刊类型引用(3)
1. 王汉鹏，张冰，李梦天，赵盛男. 思政教育融入有限元课程的教学创新与实践. 高教学刊. 2024(28): 83-86+90 . 
百度学术
2. 张升，兰鹏，苏晶晶，熊海斌. 基于PINNs算法的地下水渗流模型求解及参数反演. 岩土工程学报. 2023(02): 376-383+443-444 . 本站查看
3. 谢一凡，谢镇泽，吴吉春，张蔚，谢春红，鲁春辉. 模拟地下水流运动的三重尺度有限元模型. 岩土工程学报. 2022(11): 2081-2088 . 本站查看
其他类型引用(6)
计量
- 文章访问数: 293
- HTML全文浏览量: 27
- PDF下载量: 127
- 被引次数: 9
