模拟多孔介质中水流达西速度的Yeh‑多尺度有限元方法

摘要

本发明公开了一种模拟多孔介质中水流达西速度的Yeh‑多尺度有限元方法，运用伽辽金法将求解问题变分；将研究区剖分为粗网格单元，对所有粗网格单元进行剖分，得到细网格单元；在每一粗网格单元上求解退化椭圆方程，得到基函数；运用基函数求解变分形式，得到总刚度矩阵；根据研究区的源汇项与边界条件得到右端项；联立得到水头方程组；运用有效的数值方法求解该方程组得到研究区的节点水头；结合Yeh的伽辽金有限元模型，运用所构造的基函数和研究区的水头值，在研究区域直接求解达西方程，得到粗尺度节点上连续的达西渗透流速，再利用基函数线性表示细尺度达西渗透流速。与现有技术相比，该方法具有相近的精度和更高的效率。

说明书

技术领域

本发明属于水力学技术领域，具体涉及一种模拟多孔介质中二维水流达西速度的Yeh-多尺度有限元方法。

背景技术

地下水主要分布在多孔介质中，是水资源的重要组成部分，在模拟地下水的流动和溶质运移时，精确模拟地下水的流速和流量十分重要。因此，对地下水的达西渗透流速的计算方法的研究和数值模拟，对于考察地下水分布及溶质运移状态具有十分重要的意义。

Yeh的伽辽金有限单元模型[Yeh 1981]是一种经典求解地下水达西渗透流速的方法，该方法直接在研究区上运用有限单元法求解达西方程，具有较高的计算精度[Zhang 1994]，应用十分广泛[Boufadel et al.1999；Jang et al 2007；Park et al.2007；Park et al.2008]。许多研究证明了该方法能够保证达西渗透流速在节点上的连续性，从而保证通过截面的流入流出量相等。然而，在求解大尺度或非均质介质中的地下水问题时，该方法的有限元性质要求单元内部渗透系数是常数，必须对研究区域采用精细剖分，求解达西渗透流速时需要大量的计算时间和空间消耗，效率较低。

上世纪末，为了提升有限单元法求解非均质多孔介质中水流问题的计算效率，数学和水文研究工作者提出了多尺度有限单元法[Hou,T.Y.,and X.H.Wu(1997)]。该方法通过在粗网格单元上令基函数满足局部微分算子，抓住细尺度信息。因此，该方法单元内部渗透系数可变，无需精细剖分即可求解大尺度或非均质介质中的地下水问题，具有很高的计算效率。许多数学、水文工作者已经证明了多尺度有限单元法比有限单元法的计算效率更高，并且收敛、精确[Hou,T.Y.,and X.H.Wu(1997),X.H.Wu,and Z.Cai(1999),W.Deng et al.(2008),Ye,S.,Y,Xue,and C.Xie(2004)]。然而，多尺度有限单元法缺乏求解达西渗透流速的手段，无法保证速度和流量在节点上的连续性。

为了解决上述问题，本发明将Yeh的伽辽金模型与多尺度有限单元法有机结合，综合了两种方法的优点，提出了一种全新的计算达西渗透流速的方式。

发明内容

针对于上述现有技术的不足，本发明的目的在于提供一种模拟多孔介质中水流达西速度的Yeh-多尺度有限元方法，通过数值模拟得到的结果与解析解吻合，在相同单元数目的情况，计算达西渗透流速效率高于现有技术的方法。

为达到上述目的，本发明的一种模拟多孔介质中水流达西速度的Yeh-多尺度有限元方法，包括步骤如下：

(1)根据所要模拟的研究区域确定边界条件，设定网格单元尺度，剖分研究区域，得到粗网格单元，对所有粗网格单元剖分得到细网格单元；

(2)根据渗透系数以及基函数的边界条件，在每一粗网格单元上求解退化的椭圆型问题，确定基函数，形成有限元空间；

(3)计算粗网格单元的单元刚度矩阵，累加得总刚度矩阵；根据研究区域的边界、源汇项，计算右端项，形成多尺度有限元的水头方程组；

(4)使用cholesky分解法，求得研究区域上粗尺度节点的水头，结合基函数线性表示细尺度水头；

(5)结合Yeh的伽辽金模型，运用上述步骤(2)中所构造的基函数和步骤(4)中所获研究区的水头值；在所研究区域直接运用多尺度有限单元法求解达西方程；得到粗尺度节点上的达西渗透流速的方程组；

(6)使用cholesky分解法，求得研究区粗尺度节点上的达西渗透流速；

(7)通过上述步骤(2)中所构造的基函数和上述步骤(6)中所获粗尺度达西渗透流速线性表示细尺度节点上的达西渗透流速。

优选地，上述步骤(1)中形成粗网格单元的剖分采用的是三角形单元剖分。

优选地，上述步骤(1)中形成细网格单元的剖分采用的是三角形单元剖分。

优选地，上述步骤(2)、(5)中细网格单元上的渗透系数值取这个单元的所有顶点上的渗透系数的平均值。

优选地，上述步骤(3)中细网格单元上的渗透系数源汇项值取这个单元的所有顶点上的源汇项的平均值。

优选地，所述的Yeh-多尺度有限元方法具体还包括：先通过多尺度有限单元法求解研究区域上的水头数值，再运用多尺度有限单元法所构造的基函数和所获水头值在研究区上直接求解达西方程，假设K_xy＝K_yx＝0，则x方向的达西方程为：

$V_x=-K_x\frac{\partial H}{\partial x}---\left(1\right)$

假设Ψ_I为对应于I点的多尺度有限单元法的基函数，乘以(1)式两端可得：

$\int {\int }_{\Omega }{V}_x{\Psi }_{I}dxdy=-\int {\int }_{\Omega }{K}_x\frac{\partial H}{\partial x}{\Psi }_{I}dxdy,I=1,2,\mathrm{...}{N}_n---\left(2\right)$

其中N_n是研究区上的总节点数目；

根据Yeh的伽辽金模型和多尺度有限单元法的基本理论，在任意粗网格单元Δ_ijk上，达西渗透流速可以被基函数线性表示：

V_x＝V_x(i)Ψ_i+V_x(j)Ψ_j+V_x(k)Ψ_k>

将(3)式代入(2)式，再将(2)式离散到粗网格单元上，得到每个粗网格单元上的表达式，由于多尺度有限单元法的基函数可以在每个细网格单元上的线性基函数表示，再将每个粗网格单元上的表达式离散到细网格上，可以得到粗网格上关于V_x的单元刚度矩阵和右端项，相加得达西渗透流速的总方程组，求解可以得到研究区粗尺度节点上的达西渗透流速值，再通过(3)式即可得到细尺度达西渗透流速值。

本发明的有益效果：

本发明提出了一种计算达西渗透流速的方式，运用多尺度有限单元法提高Yeh的伽辽金有限元模型计算达西渗透流速的效率，并继承Yeh的模型能够保证达西渗透流速和流量的连续性的特点。与现有经典算法相比，在研究区网格剖分相同时，本发明具有更高的计算效率、更高精度；与精细剖分的Yeh的伽辽金有限元相比，在细网格单元数目相同时，两方法精度相近，但计算消耗不到Yeh方法的1％。本发明能高效、精确求解多种情况下的非均质多孔介质地下水达西渗透流速场，在处理非稳定流和非线性等高计算问题时，该方法的效率优势更加明显。

通过对多孔介质下的二维振荡介质稳定流模型、二维渐变介质非稳定流模型、二维振荡水头稳定流模型、二维冲积平原稳定流模型、二维稳定流潜水介质模型的数值模拟，本发明所获结果与解析解吻合的很好。与几种经典求解达西渗透流速的方法相比，在研究区网格节点数目相同时，本发明获得的粗尺度达西速度的精度更高；在总节点数目(包括粗网格上的节点)相等时，本发明计算细尺度达西速度的效率更高。

附图说明

图1为Yeh-多尺度有限元方法的研究区剖分示意图。

图2为Yeh-多尺度有限元方法的粗网格单元剖分示意图。

图3为二维振荡介质稳定流模型，各数值法在截面y＝0.6米处的水头值。

图4为二维振荡介质稳定流模型，各数值法获得的在截面y＝0.6米处的相对误差示意图。

图5为二维渐变介质非稳定流模型，各数值法获得的在截面y＝0.6米处的绝对误差示意图。

图6为二维渐变介质非稳定流模型，各数值法获得的在截面y＝0.6米处的绝对误差示意图。

图7为本发明方法流程图。

具体实施方式

为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。

参照图7所示，本发明的一种模拟多孔介质中水流达西速度的Yeh-多尺度有限元方法，包括步骤如下：

(1)根据所要模拟的研究区域确定边界条件，设定网格单元尺度，剖分研究区域，得到粗网格单元，对所有粗网格单元剖分得到细网格单元；

(2)根据渗透系数以及基函数的边界条件，在每一粗网格单元上求解退化的椭圆型问题，确定基函数，形成有限元空间；

(3)计算粗网格单元的单元刚度矩阵，累加得总刚度矩阵；根据研究区域的边界、源汇项，计算右端项，形成多尺度有限元的水头方程组；

(4)使用cholesky分解法，求得研究区域上粗尺度节点的水头，结合基函数即可线性表示细尺度水头；

(5)结合Yeh的伽辽金模型，运用上述步骤(2)中所构造的基函数和步骤(4)中所获研究区的水头值；在所研究区域直接运用多尺度有限单元法求解达西方程；得到粗尺度节点上的达西渗透流速的方程组；

(6)使用cholesky分解法，求得研究区粗尺度节点上的达西渗透流速；

(7)通过上述步骤(2)中所构造的基函数和上述步骤(6)中所获粗尺度达西渗透流速线性表示细尺度节点上的达西渗透流速。

其中，上述步骤(1)中形成粗网格单元的剖分采用的是三角形单元剖分。

其中，上述步骤(1)中形成细网格单元的剖分采用的是三角形单元剖分。

其中，上述步骤(2)、(5)中细网格单元上的渗透系数值取这个单元的所有顶点上的渗透系数的平均值。

其中，上述步骤(3)中细网格单元上的渗透系数源汇项值取这个单元的所有顶点上的源汇项的平均值。

本发明先通过多尺度有限单元法求解研究区域上的水头数值，再运用多尺度有限单元法所构造的基函数和所获水头值在研究区上直接求解达西方程，假设K_xy＝K_yx＝0，则x方向的达西方程为：

$V_x=-K_x\frac{\partial H}{\partial x}---\left(1\right)$

假设Ψ_I为对应于I点的多尺度有限单元法的基函数，乘以(1)式两端可得：

$\int {\int }_{\Omega }{V}_x{\Psi }_{I}dxdy=-\int {\int }_{\Omega }{K}_x\frac{\partial H}{\partial x}{\Psi }_{I}dxdy,I=1,2,\mathrm{...}{N}_n---\left(2\right)$

其中N_n是研究区上的总节点数目；

根据Yeh的伽辽金模型和多尺度有限单元法的基本理论，在任意粗网格单元Δ_ijk上，达西渗透流速可以被基函数线性表示：

V_x＝V_x(i)Ψ_i+V_x(j)Ψ_j+V_x(k)Ψ_k>

将(3)式代入(2)式，再将(2)式离散到粗网格单元上，得到每个粗网格单元上的表达式，由于多尺度有限单元法的基函数可以在每个细网格单元上的线性基函数表示，再将每个粗网格单元上的表达式离散到细网格上，可以得到粗网格上关于V_x的单元刚度矩阵和右端项，相加得达西渗透流速的总方程组，求解可以得到研究区粗尺度节点上的达西渗透流速值，再通过(3)式即可得到细尺度达西渗透流速值。

下面结合具体实施例对本发明做进一步的解释，其中涉及一些简写符号，以下为注解：

N：研究区边界被等分份数；

V_x：x方向上的达西渗透流速；

研究区网格节点上的达西渗透流速V_x，即粗尺度达西渗透流速；

粗单元网格节点上的达西渗透流速V_x，即细尺度达西渗透流速；

FEM：传统有限单元法；

FEM-F：精细剖分的传统有限单元法；

MSFEM-Y：Yeh-多尺度有限元方法；

Method-Yeh：Yeh的线性伽辽金模型，采用FEM求解水头；

Method-Yeh-F：Yeh的线性伽辽金模型(精细剖分)，采用FEM-F求解水头；

Method-Zhang：三次样条法，采用LFEM求解水头；

MSFEM-Y的基函数使用振荡边界条件。此外，在计算V_x的相对误差时，若在某节点V_x，且数值解V_x＜10^-10，则将该节点的相对误差记为0％；否则记为100％。

实施例1:二维振荡介质模型

研究区为一正方形单元：Ω＝[0,1m]×[0,1m]，水流方程为：

$-\frac{\partial }{\partial x}\left(K\frac{\partial H}{\partial x}\right)-\frac{\partial }{\partial y}\left(K\frac{\partial H}{\partial y}\right)=W,---\left(4\right)$

渗透系数为：

$K=\frac{1}{2+P_1\sin [\pi (x+y)·P_2]}---\left(5\right)$

其中P₁＝1.965,P₂＝1.99，源汇项W及定水头边界由解析解H＝xy(1-x)(1-y)确定。

采用MSFEM-Y、Method-Yeh、Method-Zhang和Method-Yeh-F求解此模型。为了保证粗网格节点数目相同，MSFEM-Y、Method-Yeh和Method-Zhang均将研究区的每边剖分为相同的30份(N＝30)，共1800(2NN)个三角形单元(图1)。MSFEM-Y再将每个单元剖分为25个细网格单元(图2)。为了保证单元数和MSFEM-Y的细网格单元数目相同，Method-Yeh-F将研究区每边剖分为150份(N＝150)，共45000个细网格单元。

粗尺度达西渗透流速：

MFEM-Y、Method-Yeh(Method-Zhang)和Method-Yeh-F的水头值分别由MSFEM,FEM,FEM-F求得，其中各方法在截面y＝0.6米的水头值如图3所示。从图中可以看出，MSFEM和FEM-F的水头均与解析解十分接近，精度高于FEM所获水头。

在求得水头之后，我们将使用上述四种方法求解达西速度，其中MSFEM-Y、Method-Yeh和Method-Zhang计算的在截面y＝0.6米的平均相对误差如图4所示。从图中可以看出，MSFEM-Y的曲线形状与Method-Yeh一致，但精度高于Method-Yeh和Method-Zhang。当x＝0.5米时，解析解的速度值为0，会产生较大相对误差。此时，MSFEM-Y仅在x＝0.5米处有一较高误差，而Method-Yeh和Method-Zhang均具有两处较大误差。本实施例中，Method-Yeh和Method-Zhang计算H和的CPU时间分别为3秒和2秒。

细尺度达西渗透流速

同时，MSFEM-Y不仅可以计算还可以计算如图2所示，每个粗网格单元有21个细尺度节点，则MSFEM-Y需要计算37800个细尺度节点。由于粗网格单元的部分边界互相重合。研究区内仅有22801个不重复的细尺度节点，因此我们在表1中仅仅比较MSFEM-Y和Method-Yeh-F在这22801个细尺度节点的值。

表1

如表1所示，MSFEM-Y的精度略低于Method-Yeh-F。因为MSFEM-Y的是通过(在30×30的网格上求解)和基函数(在5×5网格上求解)获得的，而Method-Yeh-F则是在150×150的网格上求解同时，MSFEM-Y的所需的CPU时间远小于Method-Yeh-F，其计算水头的时间仅为Method-Yeh-F的1％，而计算达西渗透流速的时间仅为Method-Yeh-F的0.04％，显示了MSFEM-Y具有更高的计算效率。

实施例2：二维渐变介质非稳定流模型

研究区为一正方形单元：Ω＝[0,1m]×[0,1m]，水流方程为：

$S\frac{\partial H}{\partial t}-\frac{\partial }{\partial x}\left(K\frac{\partial H}{\partial x}\right)-\frac{\partial }{\partial y}\left(K\frac{\partial H}{\partial }\right)=W,---\left(6\right)$

渗透系数K(x,y)＝(1+x)(1+y)m/d。贮水系数为0.1/m，含水层厚度为1米，时间步长为1天，总时间为5天。源汇项W及定水头边界由解析解：H＝xy(1-x)(1-y)e^-t确定。

采用MSFEM-Y、Method-Yeh、Method-Zhang和Method-Yeh-F求解此模型。剖分与上述实施例一相同，即MSFEM-Y、Method-Yeh和Method-Zhang将研究区剖分为1800个单元，MSFEM-Y将每个粗网格剖分为25个细单元，Method-Yeh-F将研究区剖分45000个单元。

粗尺度达西渗透流速

和上述实施例一相同，MSFEM-Y和Method-Yeh-F获得了比其他方法更精确的水头值，MSFEM-Y、Method-Yeh和Method-Zhang在y＝0.6米处的的绝对误差如图5所示。可以看出，MSFEM-Y的曲线位于Method-Yeh和Method-Zhang的曲线下方，绝对误差要小很多。显示MSFEM-Y比Method-Yeh和Method-Zhang在研究区剖分相同时具有更高的计算精度。

细尺度达西渗透流速

由MSFEM-Y和Method-Yeh-F计算的在截面y＝0.6米处绝对误差如图6所示。可以看出MSFEM-Y和Method-Yeh-F两方法的精度非常接近，且误差值非常低。由于MSFEM的是由得到的，对比图5、图6，可以发现MSFEM-Y在两图的曲线形状相似。本实施例中MSFEM-Y使用246秒计算水头，3秒计算和而Method-Yeh-F则需24020秒计算水头，441秒计算和上述实施例一相同，MSFEM-Y所需的时间不到Method-Yeh-F的1％。这一结果证明了MSFEM-Y具有较高的计算效率和精度。结合上述实施例一，可以发现MSFEM-Y在求解非稳定流问题时优势更加明显。

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。