CN111507026A

CN111507026A - 一种模拟节点达西渗透流速的双重网格多尺度有限单元法

Info

Publication number: CN111507026A
Application number: CN201910828194.7A
Authority: CN
Inventors: 谢一凡; 赵文凤; 吴吉春; 鲁春辉; 叶逾; 谢春红; 谢镇泽
Original assignee: Nanjing University; Hohai University HHU
Current assignee: Nanjing University; Hohai University HHU
Priority date: 2019-09-03
Filing date: 2019-09-03
Publication date: 2020-08-07

Abstract

本发明公开了一种模拟节点达西渗透流速的双重网格多尺度有限单元法，包括：将研究区剖分成粗网格单元，再将粗网格单元剖分成细网格单元；运用伽辽金法得到所求问题的变分形式；在每一粗网格单元上构造基函数；运用基函数求解变分形式，得到总刚度矩阵；根据研究区的源汇项与边界条件，得到水头的方程组；运用有效的数值方法求解，得到研究区每一节点水头值；将研究区网格沿所求达西渗透流速方向或反方向平移一段极小距离，得到平移后的网格；在平移后的网格上再次进行上述水头求解步骤，获得水头；根据每一点平移前后的水头差和位移差，得到连续的水头一阶导数，获得连续、精确的达西渗透流速；应用插值可直接获得细尺度节点的达西渗透流速。

Description

一种模拟节点达西渗透流速的双重网格多尺度有限单元法

技术领域

本发明涉及水力学领域，具体指代一种模拟节点达西渗透流速的双重网格多尺度有限单元法。

背景技术

在地下水含水层中的溶质运移的数值计算中，精确、连续的达西渗透流速能够显著提高溶质运移方程的计算精度。传统达西渗透流速有限元算法因受有限元框架限制，必须精细剖分以描述非均质介质，导致计算成本很高。因此，科学工作者提出多尺度有限单元法(MSFEM)[Hou,T.Y.,and X.H.Wu(1997)]弥补这一缺陷。该方法能够通过基函数抓住细尺度信息，无需精细剖分即可获得精确的解，具有极高的计算效率。近年来，我们证明了应用MSFEM能够改进传统有限元类达西渗透流速算法，显著提高计算效率并保持精度。

Vedat Batu的双重网格技术(D-FEM)是一种有效的有限元类达西渗透流速算法，具有原理简单，技术便于操作的特点。该方法仅仅需要进行水流方程的求解，无需专门进行达西方程的计算，即可获得连续的达西渗透流速。其主要步骤是：将原有剖分网格上的所有内部节点沿着所需求解的达西渗透流速的正向或者反向平移一段极小距离，通过平移前后的水头差和位移差，再结合达西定律即可求解出该点达西渗透流速。由于平移距离极小，平移前后的点可以视为同一点，故获得的达西渗透流速具有连续性。然而，该方法依然被有限元框架所限制，在模拟大尺度、非均质地下水问题时需要较高的计算成本。

发明内容

针对于上述现有技术的不足，本发明的目的在于提供一种模拟节点达西渗透流速的双重网格多尺度有限单元法(D-MSFEM)，以解决现有技术中传统有限元类方法为保证解的精确度需要对研究区进行精细剖分，从而产生大量计算消耗，导致计算效率较低的问题。

为达到上述目的，本发明采用的技术方案如下：

本发明的一种模拟节点达西渗透流速的双重网格多尺度有限单元法，包括以下步骤：

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

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

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

(4)提供多尺度有限元的水头方程组的有效解法，求得研究区域上各节点的水头；

(5)将研究区网格沿所需求解的达西渗透流速方向或反方向平移一段距离，得到平移后的网格，在平移后的网格上再次应用步骤(1)-(4)，求解得到平移之后各节点上的水头；

(6)利用节点平移前后的水头差、位移差，根据达西定律，结合双重网格技术得到粗网格节点的达西渗透流速；

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

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

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

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

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

优选地，所述步骤(5)中的一段距离需小于粗网格尺度的百分之一。

本发明的有益效果：

本发明运用多尺度有限单元法提高了Batu的D-FEM计算达西渗透流速的效率，并继承了D-FEM能够保证达西渗透流速和流量的连续性特点。与现有传统有限元类方法相比，在研究区剖分网格相同时，本发明具有更高的计算效率、计算精度；与精细剖分的D-FEM、Yeh伽辽金模型相比，在细网格单元数目相同时，三种方法精度相近，但计算消耗远小于Batu方法和Yeh伽辽金模型。

本发明能高效、精确地模拟多种情况下的地下水含水层中的达西渗透流速场。通过对地下水二维均质介质稳定流模型、地下水二维高振荡水头稳定流模型以及地下水二维渐变介质非稳定流模型的数值模拟，本发明所获结果能很好契合解析解。与多种达西渗透流速的传统方法相比，在研究区网格节点(粗网格上的节点)数目相等时，本发明获得的粗尺度节点达西渗透速度的精度更高；在总节点数目(包括粗网格和细网格上的节点)相同时，本发明计算细尺度达西渗透速度的效率更高。

附图说明

图1为本发明方法的原理图。

图2为双重网格多尺度有限元法的研究区第一重(原)网格示意图。

图3为双重网格多尺度有限元法的研究区第二重(平移后)网格示意图。

图4为二维均质介质模型中各数值方法所计算的水头平均相对误差与粗网格单元尺度关系示意图。

图5为二维均质介质模型中各数值方法所需计算时间与粗网格单元尺度关系示意图。

图6为二维均质介质模型中各数值方法所计算的Vx平均相对误差与粗网格单元尺度关系示意图。

图7为二维高振荡水头模型中渗透系数K分布图。

图8为二维高振荡水头模型中在截面y＝0.6处各数值方法所计算达西渗透流速的绝对误差示意图。

图9为二维渐变介质非稳定流模型中各数值方法结果对比示意图。

具体实施方式

参照图1所示，本发明的一种模拟节点达西渗透流速的双重网格多尺度有限单元法，是将双重网格技术与多尺度有限单元法有机结合提出的高效达西渗透流速算法。D-MSFEM的主要思想是：先应用多尺度网格剖分研究区域，运用MSFEM高效地求解各节点的水头值；再沿所求达西渗透流速方向(或反方向)的坐标轴将研究区内部网格平移一段极小距离，并运用MSFEM在平移后的网格上再次求解水头，获得平移前后的每个网格节点的位移差和水头差；最后应用双重网格技术获得连续的水头一阶导数，从而通过达西定律获得连续、精确的达西渗透流速。在实施例中，本发明将D-MSFEM与D-FEM以及Yeh的伽辽金有限元模型等传统方法进行对比，得到结论为：D-MSFEM能够仅应用和传统方法相近的时间即可获得精细剖分网格上的所有水头值和达西渗透流速值，具有极高的计算效率和精度。

以求解粗网格单元□_ijkl(□_ijkl表示一个四顶点分别为i,j,k,l的四边形单元)上k点在x方向上的达西渗透流速为例，图2展示了D-MSFEM的第一重网格构造，将研究区剖分成N×N个粗网格单元，设研究区域原点为(0，0)，用粗实线表示。每个粗网格单元再剖分成8个细网格单元，用细虚线表示。图2只给出了边界上的单元剖分，中间单元省略，用点虚线表示。其中□_ijkl为一示例粗网格单元。在D-MSFEM的第一重网格基础上，运用多尺度有限单元法求解出点k的水头值为H₁。

保持研究区边界条件不变的情况下，内部垂向的粗网格沿着x坐标轴方向平移一段极小距离Δx，构成双重网格多尺度有限元法的第二重网格构造(图3)，粗实线代表平移之前的粗网格边界，粗虚线代表平移之后的粗网格边界，细虚线代表对平移后的粗网格精细剖分。图3中，粗网格单元□_ijkl经过平移之后变成□_ij′k′l，点k平移后的变成点k′。主要思想是，保持研究区边界条件不变，内部垂向的粗网格线沿着坐标轴x方向向右平移一段极小距离Δx。左右边界条件不变，平移之后靠左边界的第一列单元的水平尺度扩大Δx，靠右边界的最后一列单元的水平尺度缩小Δx，中间单元大小不变，位置向右平移了Δx。平移之后的研究区内部粗网格(粗虚线)和研究区外边界重新构成了一组N×N的粗网格单元，将每个粗网格剖分成8个细网格单元，再次运用多尺度有限单元法求解出平移之后点k′的水头值为H₂，得到点k平移前后的水头差Δh。根据达西定律可求得点k在x方向上的达西渗透流速为：

通过求解渗透系数均质的二维稳定流问题、渗透系数非均质振荡的二维稳定流问题，以及渗透系数渐变的二维非稳定流问题中的水头和达西渗透流速，发现本发明能够仅应用和传统方法相近的时间即可获得精细剖分网格上的所有水头值和达西渗透流速值，并具有很高的计算精度。

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

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

AS：解析解；

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

D-FEM：Vedat Batu双重网格有限元方法；

D-FEM-F：精细剖分的Vedat Batu双重网格有限元方法；

Yeh：Gour-Tsyh Yeh的伽辽金有限元方法；

Yeh-F：精细剖分的Gour-Tsyh Yeh的伽辽金有限元方法；

D-MSFEM：双重网格多尺度有限单元法；

D-MSFEM(a,b)：用D-MSFEM方法将研究区剖分为a个粗网格单元，每个粗网格单元被剖分为b个细网格单元；

D-MSFEM-L：D-MSFEM使用基函数线性边界条件；

D-MSFEM-O：D-MSFEM使用基函数振荡边界条件。

实施例1:二维均质介质稳定流模型

研究区域Ω为[0，1]×[0，1]，含水层厚度为1，渗透系数K＝1，解析解为H＝xy(1-x)(1-y)，源汇项W及狄利克雷边界条件由解析解给出，地下水的二维稳定流方程可表示为：

粗尺度达西渗透流速

在本例中，D-MSFEM、D-FEM和Yeh方法均将研究区剖分为N×N份，当N＝10，20，30，40时，粗网格单元尺度分别为0.1，0.05，0.03，0.025。用D-FEM和Yeh方法计算时，将研究区域分别剖分成200，800，1800，3200个三角形粗网格单元(N×N×2)与之对应的，用D-MSFEM方法计算时，将研究区域分别剖分成100，400，900，1600个正方形粗网格单元(N×N)，再将每一个粗网格单元剖分成8个三角形单元(2×2×2)。同样，用Yeh-F方法计算时，将研究区域分别剖分成800，3200，7200，12800个三角形粗网格单元，以获得与D-MSFEM相同的细网格单元数目。D-FEM和D-MSFEM方法中4种尺度对应的横向平移距离Δx分别取0.001，0.0002，0.000053，0.000016。

图4展示了水头的平均相对误差和粗网格单元尺度的关系，结果显示：当粗网格单元尺度相同时，各数值方法的水头精度从高到低依次为Yeh-F、D-FEM-F、D-MSFEM、Yeh、D-FEM，其中Yeh-F、D-FEM-F、D-MSFEM三种方法精度相近，且明显比Yeh、D-FEM精度更高。由于Yeh和D-FEM(Yeh-F和D-FEM-F)都是在有限元的基础上计算，所以当剖分份数相同时，水头误差也相同。从这一结果可以看出，D-MSFEM能够获得和精细剖分的传统方法相近的水头结果。

图5展示了不同粗网格尺度下各数值方法计算达西渗透流速所需的时间(包含了水头计算时间)。当粗网格单元尺度相同时，各数值方法的所花时间从高到低为Yeh-F、D-FEM-F、D-MSFEM、Yeh、D-FEM，其中D-MSFEM、Yeh、D-FEM明显比Yeh-F、D-FEM-F所花时间更少。当粗网格尺度为0.025时，Yeh-F和D-FEM-F两种精细剖分方法所用时间最长，分别为287s、246s；D-MSFEM、Yeh、D-FEM用时较短依次为10s、10s、8s。

图6展示了达西渗透流速V_x的相对误差和粗网格单元尺度的关系。从该图可以看出，各数值方法的V_x精度从高到低依次为Yeh-F、D-FEM-F、D-MSFEM、Yeh、D-FEM。当粗网格尺度为0.025时，Yeh-F、D-MSFEM和D-FEM方法的V_x误差分别为0.054％、0.086％、0.356％。根据图4，粗网格尺度为0.025时，Yeh-F、D-MSFEM和D-FEM方法的水头误差分别为0.037％、0.042％、0.149％。精细剖分的Yeh-F结果比其余方法精确，显示了速度误差主要是由水头误差产生的。D-MSFEM的水头精度和Yeh-F十分接近，所以D-MSFEM相比于D-FEM和Yeh，能够得到更为精确的水头值，从而提高速度精度。

综上所述，可见D-MSFEM能够仅应用和传统方法相近的时间即可获得精细剖分网格上的所有水头值和达西渗透流速值，并具有很高的计算精度。

细尺度达西渗透流速

表1展示了粗网格单元尺度为0.025时，用D-MSFEM(1600,50)计算的某一粗网格单元□_ijkl中细尺度节点速度V_x及其解析解，并比较相对误差。所取粗网格单元□_ijkl的顶点坐标为i(0.8,0.6)，j(0.825,0.6)，k(0.825,0.625)，l(0.8,0.625)。

从表1可看出，D-MSFEM的细尺度达西渗透流速值具有较高精度。同时，这些细尺度达西渗透流速是D-MSFEM直接应用粗尺度节点的达西渗透流速值和多尺度基函数直接插值得到的，无需额外的计算消耗。表1如下：

表1

节点坐标	解析解	D-MSFEM	相对误差/％
				(0.805,0.605)	0.149241	0.145775	2.38
(0.81,0.605)	0.1493	0.148165	0.77
				(0.815,0.605)	0.149277	0.150554	0.85
(0.82,0.605)	0.149343	0.152944	2.35
				(0.805,0.61)	0.14853	0.145119	2.35
(0.81,0.61)	0.148604	0.147498	0.75
				(0.815,0.61)	0.148565	0.149877	0.88
(0.82,0.61)	0.148643	0.152256	2.37
				(0.805,0.615)	0.147825	0.144433	2.35
(0.81,0.615)	0.147902	0.146801	0.75
				(0.815,0.615)	0.147863	0.149168	0.87
(0.82,0.615)	0.147937	0.151536	2.37
				(0.805,0.62)	0.147125	0.143716	2.37
(0.81,0.62)	0.147191	0.146072	0.77
				(0.815,0.62)	0.147167	0.148428	0.85
(0.82,0.62)	0.147226	0.150784	2.36

实施例2：二维高振荡水头稳定流模型

渗透系数振荡的二维稳定流方程为(2)，渗透系数

研究区域Ω为[0，1]×[0，1]，含水层厚度为1，解析解为H＝xy(1-x)(1-y)，源汇项W及狄利克雷边界条件由解析解给出。令振幅A＝1.99，相位Q＝1.5，则渗透系数最大值为最小值的400倍，图7为此非均质介质中渗透系数的三维分布示意图。

本例用Yeh-F、D-MSFEM-O、D-MSFEM-L、Yeh以及D-FEM法求解。粗网格尺度取0.025，D-FEM、Yeh将研究区域分别剖分成3200个三角形粗网格单元(N×N×2)。D-MSFEM将研究区域剖分成1600个正方形粗网格单元(N×N)，再将每一个粗网格单元剖分成8个三角形单元。同样，为获得与D-MSFEM相同的细网格单元数目，用Yeh-F将研究区域剖分成12800个三角形粗网格单元。其中D-FEM和D-MSFEM的横向平移距离Δx均取0.0003。

图8展示了各数值方法所计算的截面y＝0.6上各点的达西渗透流速V_x的绝对误差。截面y＝0.6上各点的V_x的平均绝对误差从小到大依次为Yeh-F，0.88；D-MSFEM-O，2.62；D-MSFEM-L，3.04；Yeh，4.51；D-FEM，10.94。此外，D-MSFEM运用基函数的振荡边界条件的D-MSFEM-O达西渗透流速平均绝对误差为2.62，比应用基函数的线性边界条件时的D-MSFEM-L的3.04更精确。在计算时间上，Yeh-F用时265s，D-MSFEM-L和D-MSFEM-O均为18s，Yeh用时7s，D-FEM用时6s。

由此可见D-MSFEM也能高效、精确的求解渗透系数振荡的二维稳定流问题，并且具有极高的效率。同时，D-MSFEM-O能够在参数振荡的条件下获得更加精确的计算结果。

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

非稳定流的一般方程由抛物型方程表示，其二维形式为：

式中，K为渗透系数，H为水头，W为源汇项，S为贮水系数，Ω为研究区。

研究区域为[0，10km]×[0，10km]，边界条件为上下隔水边界，左右为一类边界，左边界水头为10m，右边界水头为0m，源汇项为0。渗透系数从左至右缓慢增加，具体为

m/day，是典型的山前冲积平原介质的变化特征。时间步长为1day，总时间为5day。由于此算例没有解析解，用Yeh-F将研究区精细剖分成200×200份，获得80000个三角形单元，所计算的达西渗透流速值作为解析解，以此求得各数值方法的达西渗透流速误差。

本例采用Yeh、D-FEM、D-MSFEM以及Yeh-F方法计算，将研究区剖分成10×10份。则D-FEM和Yeh将研究区剖分为200个三角形单元；D-MSFEM将研究区剖分为100个粗网格单元，每个粗网格单元再剖分成8个三角形单元，共800个三角形单元；Yeh-F将研究区剖分成800个三角形单元。其中D-FEM、D-MSFEM的横向平移距离Δx取20m。

图9展示了上述数值方法所计算的达西渗透流速的平均绝对误差，精度从高到低依次为Yeh-F，0.87；D-MSFEM，0.92；Yeh，1.61；D-FEM，1.75。在每个时间步长内，D-MSFEM的计算时间依然和传统有限元方法接近，所用时间小于精细剖分的Yeh-F。由此可见，D-MSFEM同样能精确、高效地求解非稳定流问题。

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

Claims

1.一种模拟节点达西渗透流速的双重网格多尺度有限单元法，其特征在于，包括以下步骤：

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

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

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

（4）提供多尺度有限元的水头方程组的有效解法，求得研究区域上各节点的水头；

（5）将研究区网格沿所需求解的达西渗透流速方向或反方向平移一段距离，得到平移后的网格，在平移后的网格上再次应用步骤（1）-（4），求解得到平移之后各节点上的水头；

（6）利用节点平移前后的水头差、位移差，根据达西定律，结合双重网格技术得到粗网格节点的达西渗透流速；

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

2.根据权利要求1所述的模拟节点达西渗透流速的双重网格多尺度有限单元法，其特征在于，所述步骤（1）中形成粗网格单元的剖分采用的是矩形单元剖分。

3.根据权利要求1所述的模拟节点达西渗透流速的双重网格多尺度有限单元法，其特征在于，所述步骤（1）中形成细网格单元的剖分采用的是三角形单元剖分。

4.根据权利要求1所述的模拟节点达西渗透流速的双重网格多尺度有限单元法，其特征在于，所述步骤（2）中，细网格单元上的渗透系数值取其所有顶点上的渗透系数的平均值。

5.根据权利要求1所述的模拟节点达西渗透流速的双重网格多尺度有限单元法，其特征在于，所述步骤（3）中，细网格单元上的渗透系数、源汇项值取这个单元的所有顶点上的渗透系数、源汇项的平均值。

6.根据权利要求1所述的模拟节点达西渗透流速的双重网格多尺度有限单元法，其特征在于，所述步骤（5）中的一段距离需小于粗网格尺度的百分之一。