CN110083853B

CN110083853B - 模拟地下水流运动的有限体积Yeh多尺度有限元法

Info

Publication number: CN110083853B
Application number: CN201811146786.2A
Authority: CN
Inventors: 谢一凡; 吴吉春; 鲁春辉; 叶逾; 薛禹群; 谢春红
Original assignee: Nanjing University; Hohai University HHU
Current assignee: Nanjing University; Hohai University HHU
Priority date: 2018-09-29
Filing date: 2018-09-29
Publication date: 2022-09-20
Anticipated expiration: 2038-09-29
Also published as: CN110083853A

Abstract

本发明公开了一种模拟地下水流运动的有限体积Yeh多尺度有限元法，包括：确定粗、细单元的尺度；将研究区剖分为粗单元；将每一粗单元剖分为细单元；以粗网格上的每一未知节点为中心，连接与该节点相关粗单元的中心，将研究区剖分为互不重叠的体积元；在每一粗单元上构造基函数，构造x和y方向速度矩阵；在每一体积元上将水流方程积分，应用散度定理进行变换，通过速度矩阵将达西速度项应用水头的粗尺度解线性表示，获得该体积元粗尺度水头解的方程，结合差分格式，相加水头的总方程；采用有效的计算方法求解总方程，同时通过速度矩阵获得达西速度值。与多种经典方法相比，本发明的有限体积Yeh多尺度有限元法具有更高的计算效率。

Description

模拟地下水流运动的有限体积Yeh多尺度有限元法

技术领域

本发明属于水力学技术领域，具体涉及一种模拟地下水流运动的有限体积Yeh多尺度有限元法。

背景技术

21世纪以来，地下水污染、海水入侵等问题的重要性与日俱增，建立精确高效的水流模型综合分析水头、达西速度等参数具有重要意义。精确高效的水头、达西速度解法能够精确描述地下水位及运移状态，能够令地下水工作者能更好的了解地下水系统情况，具有重要的研究意义。

地下水系统的非均质性常跨越很多尺度，应用传统有限元法直接进行水头求解非常困难，需要非常大的计算消耗。因此，科研工作者提出了多尺度有限单元法(Hou and Wu1997)改进了传统有限元法，通过在粗单元上构造多尺度基函数抓住介质的细尺度信息，并直接在粗尺度上求解水头，能够大幅降低计算消耗。随后，科学工作者在“Finite volumemultiscale finite element method for solving the groundwater flow problems inheterogeneous porous media”一文中提出了有限体积多尺度有限元法(He and Ren)，将有限体积法和多尺度有限单元法有机结合，可以通过有限体积理论保证水流的局部质量守恒，具有比多尺度有限单元法更高的计算精度。

然而，科学工作者只给出了有限体积多尺度有限元法的水流方程的解法，如何应用有限体积多尺度有限元法计算达西速度还有待研究。另一方面，有限元等算法无法保证非均质地下水问题中达西速度的连续性，精度受到了限制。Yeh在1981年于“On thecomputation of Darcian velocity and mass balance in the finite elementmodeling of groundwater flow”一文中提出的伽辽金有限元模型可以解决达西速度连续性的问题，具有比同类方法更高的精度和应用范围。然而，这种方法是基于有限元法的，在模拟非均质介质中的地下水问题的达西速度时也会出现计算消耗过大的问题。此外，在水头和达西速度都需要计算的情况时，大多数算法需要求解不同的方程进行计算，步骤繁琐。

发明内容

本发明的目的在于提供一种模拟地下水流运动的有限体积Yeh多尺度有限元法，以解决现有技术中求解非均质地下水头和连续达西速度的计算效率过低的问题；本发明方法组合了有限体积多尺度有限元法和Yeh的伽辽金有限元模型，通过多尺度基函数提高水头和达西速度的计算效率，通过Yeh的伽辽金有限元模型构造速度矩阵保证达西速度的连续性，通过将该速度矩阵嵌入有限体积多尺度有限元法的算法框架实现达西速度和水头的同时计算，能够保证达西速度和水头的质量守恒性。

为达到上述目的，本发明的一种模拟地下水流运动的有限体积Yeh多尺度有限元法，包括步骤如下：

(1)根据研究区域确定所要模拟的地下水问题的边界条件，设定粗尺度，对该研究区域进行网格剖分，得到粗单元；

(2)设定细尺度，对上述每一粗单元进行网格剖分，得到细单元；

(3)根据步骤(1)中的粗尺度，以步骤(1)中获得的研究区剖分网格上的每一未知节点为中心，连接与该节点相关粗单元的中心，获得该节点的体积元，将研究区剖分为互不重叠的矩形体积元；

(4)在步骤(2)中的网格剖分下，根据渗透系数K、粗单元顶点上的多尺度基函数值以及多尺度基函数边界条件公式，在粗单元上求解退化的椭圆型问题，获得多尺度基函数在粗单元所有节点上的值；

(5)在步骤(2)中的网格剖分下，根据渗透系数K、粗单元上的基函数，结合Yeh的伽辽金有限元模型框架，在粗单元上求解x方向上的达西定律方程，获得粗单元上关于V_x的速度矩阵；

(6)在步骤(2)中的网格剖分下，根据渗透系数K、粗单元上的基函数，结合Yeh的伽辽金有限元模型框架，在粗单元上求解y方向上的达西定律方程，获得粗单元上关于V_y的速度矩阵；

(7)在步骤(3)中每一个体积元上将水流方程积分，根据有限体积多尺度有限元法中水头的粗尺度解的定义，将水头的粗尺度解对时间偏微分导数代入水流方程，应用散度定理对水流方程进行变换；

(8)将步骤(7)中的经过散度定理变换得到的水流方程离散到与该体积元相关的粗单元上，在每个粗单元上通过速度矩阵将水流方程中的达西速度项应用水头的粗尺度解线性表示，获得该体积元上关于水头的粗尺度解的方程；

(9)结合Crank-Nicolson格式，联立所有体积元上关于水头的粗尺度解的方程，获得关于水头的粗尺度解的总方程，采用cholesky分解法，求得研究区域上每个节点的水头；

(10)在每个粗单元上，应用步骤(5)、(6)中构造的速度矩阵获得细尺度达西速度；

(11)在步骤(1)中的研究区剖分网格的每一节点上，平均从该节点的相关粗单元中按照步骤(10)获得的该节点的细尺度达西速度值，来获得该节点的粗尺度达西速度值；

(12)在步骤(1)中的研究区剖分的网格线上的每一节点上，根据步骤(2)中的粗单元剖分，平均从该节点的相关粗网格中按照步骤(10)获得的该节点的细尺度达西速度值，来获得该节点的粗尺度达西速度值。

优选地，上述的步骤(1)中，采用矩形单元网格剖分研究区域，以形成粗单元。

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

优选地，上述的步骤(12)中，所述的网格线上的每一节点不包括步骤(11)中已经计算过的节点。

优选地，上述的步骤(12)中的网格线上的节点为该网格线两侧的粗单元按步骤(2)的网格剖分而产生在该网格线上的共同节点。

本发明的有益效果：

1.本发明可以同时获得水头和两个方向上的粗、细尺度达西速度，不需要额外的计算消耗计算达西速度。

2.本发明可以保证达西速度的连续性和守恒性。

3.本发明能够获得和精细剖分的有限元法、有限体积多尺度有限元法几乎一致的水头精度，计算时间和有限体积多尺度有限元法一致，远小于精细剖分的有限元法。

4.本发明能够获得和精细剖分的Yeh的伽辽金有限元模型相近的达西速度的精度，但计算消耗远小于该方法。

5.本发明能够高效求解稳定流和非稳定流问题，可以有效处理连续、渐变、随机对数正态分布变化的非均质介质。

附图说明

图1为本发明有限体积Yeh多尺度有限元法的研究区剖分示意图。

图2为本发明有限体积Yeh多尺度有限元法的粗单元剖分示意图。

图3为二维渐变介质非稳定流模型，各方法在y＝6000m截面处的水头值示意图。

图4为二维渐变介质非稳定流模型，各方法在y＝6000m截面处的粗尺度达西速度值示意图。

图5为二维渐变介质非稳定流模型，各方法在y＝6125m截面处的细尺度达西速度值示意图。

图6为二维渐变介质非稳定流模型，各方法计算细尺度达西速度的计算消耗对比示意图。

图7为二维随机对数正态分布介质稳定流模型，各方法在y＝600m截面处的水头值示意图。

图8为二维随机对数正态分布介质稳定流模型，LFEM-F，FVYMSFEM-O(grid 1)，FVMSFEM-O(grid 1)计算水头的计算消耗对比示意图。

具体实施方式

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

本发明的一种模拟地下水流运动的有限体积Yeh多尺度有限元法，包括步骤如下：

以矩形区域Ω为例，步骤(1)中的网格剖分为按照图1的细实线对研究区进行的剖分，以获得粗单元；步骤(2)中的网格剖分为按照图2的对每一粗单元细分，以获得细单元。步骤(3)的体积元剖分以节点i为中心的体积元D_i为例，连接与节点i相关的粗单元中心A、B、C、D，获得体积元D_i。

多尺度基函数的构造过程以Ψ_i为例，在矩形粗单元□_ijkl(图2)上考虑退化的椭圆型方程：

其中，Ψ_i为与i点相关的多尺度基函数，K为渗透系数，上式的边界条件可以通过有限元的线性、振荡边界条件得到。将上式(1)左右乘以每个未知节点上的有限元线性基函数并通过伽辽金变分，再将式(1)离散到细单元上，并应用有限元线性基函数将Ψ_i展开，即可得到关于Ψ_i的方程组，求解可以得到Ψ_i在粗单元□_ijkl上所有节点的值。

根据有限体积多尺度有限元法(He and Ren 2005)的工作，水头的粗尺度解的定义为：

其中，|D_i|为体积元D_i的面积，H为水头，Φ_i＝Φ(x,y,t)为点i的水头的粗尺度解。

速度矩阵的构造过程以V_x为例，在矩形粗单元□_ijkl上考虑达西定律方程：

其中K_x为x方向的渗透系数。

应用有限元线性基函数N_I分别乘以达西定律方程的两端，并积分，可得：

其中，n为粗单元□_ijkl上的总结点数目。

根据有限体积多尺度有限元法(He and Ren 2005)的工作，在粗单元□_ijkl上

H(x,y)＝Φ_iΨ_i+Φ_jΨ_j+Φ_kΨ_k+Φ_lΨ_l (5)

其中，Ψ_i，Ψ_j，Ψ_k，Ψ_l分别为节点i,j,k,l的多尺度基函数。

根据Yeh的伽辽金模型(Yeh 1981)，在□_ijkl的任意子单元Δ_abc上

V_x(x,y)＝V_x(a)N_a+V_x(b)N_b+V_x(c)N_c (6)

其中，N_a，N_b，N_c分别为节点a,b,c的线性基函数。

将式(4)离散到细单元上，并代入式(5)、(6)可以得到：

[A^x][V_x]＝[B^x][Φ] (7)

矩阵[A^x]可逆，因此可以获得速度矩阵[α^x]＝[A^x]^-1[B^x],由此可得：

[V_x]＝[α^x][Φ] (8)

本发明的有限体积Yeh多尺度有限元法的具体应用过程以二维非稳定流为例，考虑如下水流方程：

以未知节点i的体积元D_i上的积分过程为例，将式(9)在体积元D_i上积分，并应用散度定理，结合式(2)可以得到：

其中，

为贮水系数S在D_i上的值，为常数。

将式(10)离散到与i相关的粗单元上，结合达西定律，并将速度矩阵带入式(10)，可以得到D_i上关于水头粗尺度解的方程，根据有限体积法，结合Crank-Nicolson格式，联立研究区上所有体积元的方程，即可得到研究区上关于水头粗尺度解的方程组，求解后可以得到研究区上粗尺度网格上每个节点的水头值。在每一粗单元上，通过基函数插值可以得到细尺度水头。在每一粗单元上，根据粗尺度水头和速度矩阵可以获得细尺度的达西速度。再根据步骤(11)、(12)即可以得到粗尺度达西速度。

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

AS:解析解

LFEM:有限单元法；

LFEM-F:有限单元法(精细剖分)；

Method-Yeh:Yeh的有限单元模型；

Method-Yeh-F:Yeh的有限单元模型(精细剖分)；

FVMSFEM-O:有限体积多尺度有限元法,振荡边界条件；

FVMSFEM-L:有限体积多尺度有限元法，线性边界条件；

FVYMSFEM-O:本发明的有限体积Yeh多尺度有限元法,振荡边界条件；

FVYMSFEM-L:本发明的有限体积Yeh多尺度有限元法，线性边界条件；

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

研究区Ω为10km×10km的区域，渗透系数为

研究地下水非稳定流方程，即：

本例中研究区的左右边界为定水头边界，分别为10m和0m，上下边界为隔水边界，注水系数S＝10^-6/m，含水层厚度为1m。在坐标(6000m,6000m)处有一抽水井，流量为100m³/day，抽水时间为6天，时间步长为1天。

本例采用LFEM-F，FVYMSFEM-O，FVYMSFEM-L，FVMSFEM-O，FVMSFEM-L和LFEM计算。FVYMSFEM，FVMSFEM，LFEM均将研究区每一边界剖分为40等份。FVYMSFEM，FVMSFEM将研究区剖分为1600个矩形粗单元，每个粗单元剖分为32个三角形细单元；LFEM将研究区剖分为3200个三角形粗单元。LFEM-F将研究区剖分为51200个三角形细单元，即和FVYMSFEM、FVMSFEM同样数目的细单元。由于本例没有解析解，LFEM-F将研究区剖分为204800单元来获得“标准解”(Standard Solution)进行参照。

各方法在y＝6000m处的水头如图3所示。从图3中我们可以看到LFEM-F，FVMSFEM-O，FVYMSFEM-O三者的曲线几乎重合，证明了FVYMSFEM-O具有很高的精度。同时，FVMSFEM-O，FVYMSFEM-O的结果要好于FVMSFEM-L，FVYMSFEM-L，说明FVYMSFEM应用振荡的基函数边界条件能够获得更好的结果。在抽水井附近，LFEM-F，FVMSFEM-O，FVYMSFEM-O均取得了很好的结果，其中LFEM-F的精度要略高。在计算时间方面，FVYMSFEM的时间和FVYMSFEM的时间一致，要远小于LFEM-F。同时，FVYMSFEM还能不花任何额外的计算消耗获得两个方向上的达西速度。

本例的达西速度场由Method-Yeh-F，FVYMSFEM-O和Method-Yeh进行求解，其中Method-Yeh应用了LFEM计算水头时的网格，Method-Yeh-F则将研究区剖分为51200个三角形细单元，即和FVYMSFEM、FVMSFEM同样数目的细单元，并将研究区剖分为204800单元来获得“标准解”(Standard Solution)进行参照。图4为各方法在y＝6000m截面处的粗尺度达西速度，可以看到FVYMSFEM-O、Method-Yeh-F的结果和标准解十分接近，精度远高于Method-Yeh。图5为各方法在y＝6125m截面处的细尺度达西速度，虽然Method-Yeh-F的精度略高于FVYMSFEM-O，但两者相差不大。同时，两方法的解均和标准解十分接近，显示了FVYMSFEM-O具有很高的精度。在图6中，我们还将两方法的计算消耗进行了对比。由于达西速度仅在第6天进行计算，因此Method-Yeh-F计算达西速度的时间大约是计算水头的时间的六分之一。FVYMSFEM-O不需要额外的时间进行达西速度的计算，该方法计算水头V_x、V_y的时间总合也仅占了0.27％的比例。这一结果证明了FVYMSFEM具有极高的计算效率。

实施例2：二维随机对数正态分布介质稳定流模型

研究区Ω为1km×1km的区域，原点为(0m,0m)。渗透系数K_x＝K_y＝K，其中K是通过GSLib(Deutsch and Journel,1998)中的序贯高斯模拟法在400×400的网格上生成的随机对数正态分布系数场，lnK的方差为4，相关长度为150m。研究区的左右边界为定水头边界，分别为16m和11m，上下边界隔水，源汇项为0。

本例采用FVYMSFEM-O，VMSFEM-O，LFEM-F进行求解。FVYMSFEM-O，FVMSFEM-O将研究区剖分为三种不同的网格单元：grid 1:625矩形粗单元，每个粗单元剖分为128个细单元；grid 2:100矩形粗单元，每个粗单元剖分为128个细单元；grid 3:100矩形粗单元，每个粗单元剖分为32个细单元。LFEM-F将研究区剖分为80000个三角形细单元，即和grid 1同样数目的细单元。由于本例没有解析解，LFEM-F将研究区剖分为320000单元来获得“标准解”(Standard Solution)进行参照。

图7为各方法在y＝600m截面处的水头值。各方法的精度从高至低为LFEM-F，FVYMSFEM-O(grid 1)，FVMSFEM-O(grid 1)，FVYMSFEM-O(grid 2)，FVMSFEM-O(grid 2)，FVYMSFEM-O(grid 3)，FVMSFEM-O(grid 3)，其中FEM-F，FVYMSFEM-O(grid 1)，FVMSFEM-O(grid 1)三个方法的精度非常接近。这一结果证明了FVYMSFEM具有很高的精度，并且粗、细尺度的网格细分也能提高FVYMSFEM的精度。本例中LFEM-F，FVYMSFEM-O(grid 1)，FVMSFEM-O(grid 1)的计算消耗对比如图8。和上例一样，FVYMSFEM依然取得了和FVMSFEM一样的计算消耗，计算效率远高于LFEM-F。

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

Claims

1.一种模拟地下水流运动的有限体积Yeh多尺度有限元法，其特征在于，包括步骤如下：

(2)设定细尺度，对每一粗单元进行网格剖分，得到细单元；

2.根据权利要求1所述的模拟地下水流运动的有限体积Yeh多尺度有限元法，其特征在于，上述的步骤(1)中，采用矩形单元网格剖分研究区域，以形成粗单元。

3.根据权利要求1所述的模拟地下水流运动的有限体积Yeh多尺度有限元法，其特征在于，上述的步骤(2)中，采用直角三角形单元网格剖分粗单元，以形成细单元。

4.根据权利要求1所述的模拟地下水流运动的有限体积Yeh多尺度有限元法，其特征在于，上述的步骤(12)中，所述的网格线上的每一节点不包括步骤(11)中已经计算过的节点。

5.根据权利要求1所述的模拟地下水流运动的有限体积Yeh多尺度有限元法，其特征在于，上述的步骤(12)中的网格线上的节点为该网格线两侧的粗单元按步骤(2)的网格剖分而产生在该网格线上的共同节点。