CN107657075B

CN107657075B - 模拟地下水介质交界面处达西速度的区域分解有限元法

Info

Publication number: CN107657075B
Application number: CN201710732804.4A
Authority: CN
Inventors: 谢一凡; 吴吉春; 薛禹群; 谢春红
Original assignee: Nanjing University
Current assignee: Nanjing University
Priority date: 2017-08-24
Filing date: 2017-08-24
Publication date: 2019-09-03
Anticipated expiration: 2037-08-24
Also published as: CN107657075A

Abstract

本发明公开了一种模拟地下水介质交界面处达西速度的区域分解有限元法，包括：运用伽辽金法对地下水流问题变分，剖分研究区域，应用有限元法获得水头；根据研究区域的介质组成，应用不同介质的交界面将研究区域分解为若干单一介质子区域，应用区域分解法将研究区域上的地下水达西速度求解问题分解为子区域上的子问题；选取一个子问题，运用伽辽金法变分，应用Yeh的伽辽金模型获得达西速度，结合折射定律获得该子区域和其他子区域交界面上另一侧的达西速度，作为相邻子问题的第一类边界条件；选取下一个子问题求解，重复这一过程直到所有子问题求解完毕。该方法通过区域分解法降低了达西速度的计算消耗，应用折射定律保证交界面处达西速度的精度。

Description

模拟地下水介质交界面处达西速度的区域分解有限元法

技术领域

本发明属于水力学技术领域，具体涉及一种模拟地下水介质交界面处达西速度的区域分解有限元法。

背景技术

地下水是水资源的重要组成部分，是人类赖以生存的重要资源。地下水的流速和流量能够精确刻画地下水的运动状态，对于地下水问题特别是溶质运移问题具有重要意义。因此，开发精确、高效的地下水达西速度算法，对于考察地下水数值模拟具有非常重要的意义。

自然界中大多数地下水介质具有非均质性，并且地下水介质交界面、裂隙、透镜体等因素均会影响地下水介质的非均质性。在模拟许多自然界中的地下水问题时，可以根据所研究区域的介质组成，通过不同介质的交界面将研究区分为若干个具有单一介质的子区域。然而，在模拟子区域交界面处的达西速度时，包括Yeh的伽辽金有限元模型在内的许多经典算法难以保证达西速度满足折射定律这一基本法则，即无法保证在不同介质交界面处的法向达西速度连续，切向达西速度不连续。为了解决这一问题，Zhou等在2001年应用折射定律改进了Yeh的伽辽金有限元模型，通过迭代保证了达西定律在不同介质交界面处满足折射定律，能够获得比Yeh的伽辽金有限元模型更加精确的达西速度分布。

随着经济、科技的快速发展，大尺度地下水问题备受关注。地下水数值模拟时涉及的时间、空间尺度也越来越大，导致计算消耗的急剧上升。同时，在模拟此类问题时，Zhou的模型由于在处理交界面时需要迭代，降低了整体的计算效率。在计算大尺度地下水问题时存在困难。

针对这一问题，本发明将采用区域分解法提高达西速度的计算效率，降低计算消耗，并保证交界面处达西速度的精度。本发明借鉴了Yeh在1981年于“On the computationof Darcian velocity and mass balance in the finite element modeling ofgroundwater flow”一文中提出有限元模型，和Zhou等在2001年于“Accurate calculationof specific discharge in heterogeneous porous media”一文中提出的改进模型。由于区域分解法、Yeh的伽辽金模型、Zhou的改进模型在现有研究中展现了很好的有效性和稳定性。因此，本发明具有很强的可行性。

发明内容

有鉴于此，本发明的目的在于提供一种模拟地下水介质交界面处达西速度的区域分解有限元法，以解决现有技术中计算达西速度时效率较低和不同介质交界面处达西速度不满足折射定律的问题，本发明的方法通过区域分解法降低了达西速度的计算消耗，并应用折射定律保证交界面处达西速度的精度。

为达到上述目的，本发明的一种模拟地下水介质交界面处达西速度的区域分解有限元法，包括步骤如下：

(1)确定所要模拟的地下水达西速度问题所在研究区域上的地下水流问题的定解条件，设定网格尺度，剖分研究区域，得到研究区域的网格单元；

(2)应用有限元法获得研究区域上的地下水流问题的方程组，即运用伽辽金法对地下水流问题变分，在步骤(1)中的网格剖分下，构造有限元基函数，结合渗透系数K获得研究区域上的每一单元上的关于水头的单元刚度矩阵，相加得总刚度矩阵，根据研究区域的边界条件和源汇项计算右端项，形成方程组；

(3)采用cholesky分解法计算步骤(2)中获得的方程组，求得研究区域上每个网格节点的水头；

(4)根据研究区域的介质组成，运用不同介质的交界面将研究区域分解为若干个互不重叠的具有单一介质的子区域，各子区域直接应用步骤(1)中的剖分获得网格单元；

(5)确定研究区域上所要模拟的地下水达西速度问题的定解条件，即确定研究区域上的达西定律方程的定解条件；应用区域分解法将研究区域上的地下水达西速度问题分解为子区域上子问题；

(6)选取第一个需要求解的子问题，应用上述步骤(3)中的水头获得应用Yeh的伽辽金模型解该子问题时所需的水头条件；

(7)应用Yeh的伽辽金模型获得步骤(6)中选取的子问题的方程组，即在该子问题所在子区域上运用伽辽金法对该子问题变分，构造有限元基函数，结合渗透系数K获得该子问题所在子区域上的每一单元上的关于达西速度的单元刚度矩阵，相加得总刚度矩阵，根据该子问题所在子区域的边界条件和源汇项计算右端项，形成方程组；

(8)采用cholesky分解法计算步骤(7)中获得的方程组，求得步骤(6)中选取的子问题所在子区域的每个节点上的达西速度；

(9)结合折射定律获得步骤(6)中选取的子问题所在子区域和相邻子区域交界面上另一侧的达西速度值，该达西速度值将作为相邻子区域上的子问题的第一类边界条件；

(10)选取下一个需要求解的子问题，应用上述步骤(3)中的水头获得应用Yeh的伽辽金模型解该子问题时所需的水头条件；若该子问题所在子区域某边界上的达西速度值已通过相邻子问题结合折射定律得到，则作为第一类边界条件；

(11)应用步骤(7)-(8)获得步骤(10)中选取的子问题所在子区域上所有的达西速度值，并结合折射定律计算该子区域和相邻子区域交界面上另一侧的达西速度值，该达西速度值将作为相邻子区域上的子问题的第一类边界条件；

(12)重复步骤(10)-(11)，直到步骤(5)中的所有子问题求解完毕。

优选地，上述的步骤(1)中，采用三角形单元剖分研究区域。

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

优选地，上述的步骤(7)中，所述的渗透系数K取这个单元的所有顶点上的渗透系数平均值。

优选地，上述的步骤(2)中，所述的源汇项取这个单元的所有顶点上的源汇项平均值。

优选地，上述的步骤(7)中，所述的源汇项取这个单元的所有顶点上的源汇项平均值。

优选地，上述的步骤(11)中，所计算的交界面上另一侧的达西速度值应为之前未进行求解的。

本发明的有益效果：

1、应用区域分解法将地下水达西速度的求解问题转化为若干子问题，将整个研究区域上的高阶速度方程组转化为若干低阶方程组，大幅提高了计算效率；

2、能够应用折射定律保证不同介质交界面处的达西速度的精度；

3、在研究区域剖分相同时，本发明模拟地下水达西速度问题的精度和Zhou的模型相近，且本发明将整个研究区域上的高阶速度方程组转化为若干低阶方程组并无需迭代，计算时间更少；

4、在研究区域剖分相同时，本发明模拟地下水达西速度问题时能够保证速度满足折射定律，而Yeh的模型不能，且本发明将整个研究区域上的高阶速度方程组转化为若干低阶方程组，计算时间更少；

5、本发明能够有效处理平行于坐标轴、不平行于坐标轴、交叉的介质交界面。

附图说明

图1为区域分解有限元法的流程图；

图2为含有平行于坐标轴的交界面的研究区域的示意图；

图3为实施例1各方法在垂直交界面的截面x＝0.5处的达西速度V_x值的示意图；

图4为实施例1各方法在交界面AB处的达西速度V_x值的示意图；

图5为含有不平行于坐标轴的交界面的研究区域的示意图；

图6为实施例2各方法在交界面AB处的达西速度V_x值的示意图；

图7为实施例2各方法在交界面CD处的达西速度V_x值的示意图；

图8为当zone 1和zone 3的渗透系数为100m/d时，各方法在x＝80m处的达西速度V_x值的示意图；

图9为当zone 1和zone 3的渗透系数为10m/d时，各方法在x＝80m处的达西速度V_x值的示意图；

图10为当zone 1和zone 3的渗透系数为1000m/d时，各方法在x＝80m处的达西速度V_x值的示意图；

图11为含有交叉交界面的研究区域的示意图。

具体实施方式

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

本发明的一种模拟地下水介质交界面处达西速度的区域分解有限元法，包括步骤如下：

(12)重复步骤(10)-(11)，直到步骤(5)中的所有子问题求解完毕。

交界面折射定律：

在两种不同介质交界面处，达西速度应满足折射定律：

V⁺·n＝V^-·n

H⁺＝H^- (1)

其中V是达西速度，H为水头，n为交界面法向单位向量，“+”、“-”号分别代表交界面两侧的渗透系数。

由于水头在所有的交界面节点是连续的，其梯度J满足:

J⁺·s＝J^-·s (2)

其中s为交界面切向单位向量。

根据达西定律，可以得到下式：

根据公式(1)-(3)，折射定律可以被表示为：

V⁺·n＝V^-·n

[K^-]^-1·V^-·s＝[K⁺]^-1·V⁺·s (4)

区域分解有限元法(DDFEM)：

以二维地下水稳定流问题的求解过程为例(图1)，介绍区域分解有限元法的实施过程。设法向单位向量n＝[n_x,n_y]^T，切向单位向量s＝[n_y,-n_x]^T，由折射定律，我们可以得到交界面两侧达西速度的差异，即DDFEM的Jump函数：

DDFEM需要先获得研究区域所有节点上的水头值，然后再通过水头求解达西速度，需在研究区域Ω上求解如下方程：

应用有限元法求解上述公式(6)获得水头，即运用伽辽金法对地下水流问题变分，在步骤(1)中的网格单元剖分下，构造有限元基函数，根据渗透系数K获得每一单元上的关于水头的单元刚度矩阵，相加得总刚度矩阵；根据研究区域的边界条件和源汇项计算右端项，形成方程组，应用cholesky分解法，求得研究区域上每个节点的水头。

在获得水头分布后DDFEM可以开始达西速度的求解，以x方向的达西速度求解过程为例，在研究区域Ω上考虑如下达西定律方程：

假设研究区域由两种不同介质构成，参照图2所示，可由介质交界面AB分成zone1,zone 2两个区域。应用区域分解法，可以将(7)式的求解问题分解为zone 1和zone 2两个子区域上的(7)式的求解问题，DDFEM将应用Yeh的伽辽金模型分别进行zone 1和zone 2上的(7)式的求解，所有子区域的剖分应用的是在解水头过程中使用的网格。

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

其中n_n为zone 1上的总节点数目。

根据Yeh的伽辽金有限元模型，在任意网格单元Δ_ijk上，水头和达西速度均可由有限元线性基函数表示：

H(x,y)＝H_iN_i+H_jN_j+H_kN_k

V_x(x,y)＝V_x(i)N_i+V_x(j)N_j+V_x(k)N_k (9)

其中V_x(i),V_x(j),V_x(k)分别为i,j,k点的达西速度值，H_i,H_j,H_k分别为i,j,k点的水头。

将(9)式代入(8)式，(8)式的左右两端可以分别被表示为：

和

由(10)、(11)式，可以得到zone 1上的达西速度方程组，应用cholesky分解法，可以得到zone 1所有节点的达西速度值。

在zone 1和zone 2的交界面AB上应用jump函数公式(5)，即可获得靠zone 2一侧的达西速度值，并将作为求解zone 2上达西定律方程时的第一类边界条件。在zone 2上应用Yeh的伽辽金模型解(7)式，进行类似zone 1上的求解过程，即可得到zone 2所有节点的达西速度值。

需要说明的是，DDFEM并不限定有限元法作为水头的求解器，DDFEM也可以应用其他方法进行水头的求解。下面结合具体实施例对本发明做进一步的解释，其中涉及一些简写符号，以下为注解：

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

FEM：传统有限单元法；

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

DDFEM:区域分解有限元法，采用FEM求解水头；

DDFEM-AS:区域分解有限元法，采用解析解水头；

Method-Zhou：Zhou的改进模型(即Zhou2001年工作中的Scheme 1)，采用FEM求解水头；

Method-Zhou-F：Zhou的改进模型(精细剖分)，采用FEM-F求解水头；

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

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

实施例1：二维稳定流模型，含有平行于坐标轴的交界面

研究区如图2所示，地下水流和达西速度的方程分别为(6)、(7)式。研究区域Ω＝[0,1]×[0,1]。研究区域含有两种不同的介质，可以被交界面AB(y＝0.5)分为两个子区域，两个子区域的渗透系数分别为：zone 1:K_1x＝(1+x)(1+y),K_1y＝10(1+x)(1+y)，zone 2:K_2x＝10(1+x)(1+y),K_2y＝100(1+x)(1+y)。本例具有解析解：

水头：

V_x:

本例的源汇项W以及第一类边界条件根据解析解水头给出。

研究区被三角形网格剖分为7200个单元，水头由FEM求解，速度由DDFEM、Method-Zhou、Method-Yeh求解。

各方法在垂直交界面的截面x＝0.5处的达西速度V_x值如图3所示。由折射定律可知，交界面两侧的切向达西速度不连续，点(0.5,0.5)位于交界面AB上，应有两个不同的达西速度值。DDFEM、Method-Zhou和解析解均在点(0.5,0.5)处有2个不同的达西速度值，且两点之间的差距满足折射定律，显示了DDFEM，Method-Zhou能够有效处理交界面。另一方面，Method-Yeh在此截面上的点处处连续，在点(0.5,0.5)处存在误差。同时，Method-Yeh在点(0.5,0.5)附近的点的达西速度值会对点(0.5,0.5)处的达西速度值进行补偿，也出现了一定的误差，显示Method-Yeh在交界面附近的达西速度精度较差。

各方法在交界面AB处的达西速度V_x值如图4所示。图中DDFEM和Method-Zhou的结果十分接近，均和解析解重合，显示了很高的计算精度。和解析解一样，DDFEM和Method-Zhou均在AB的所有节点上具有两个达西速度值，且它们之间的差距满足jump函数公式(5)，证明了两种方法的确能够保证交界面的切向达西速度满足折射定律公式(4)。而Method-Yeh在AB的所有节点上均仅有一个值，其曲线大约处于解析解的两条曲线的中间位置，显示了该方法获得的达西速度值在AB上处处连续，但不满足折射定律。因此，DDFEM和Method-Zhou在处理交界面问题时，具有比经典方法Method-Yeh更高的计算精度。

在本例中，Method-Zhou需要CPU时间353s进行水头和达西速度的计算(其中达西速度需要进行15次迭代)，而Method-Yeh的时间为37s。DDFEM获得了与Method-Zhou几乎一致的结果，但仅需要17s计算水头，5s计算达西速度，即22s CPU时间。因此，DDFEM比其他两种方法具有更高的计算效率，同时能够保证交界面处的达西速度具有很高的计算精度。

实施例2：二维稳定流模型，含有不平行于坐标轴的交界面

研究区域如图5所示，地下水流和达西速度的方程分别为(6)、(7)式。研究区域Ω＝[0,120m]×[0,120m]，被介质交界面AB(y＝20m)和CD(x+y＝160m)分为3个子区域，源汇项W为0。研究区域左右边界为定水头边界分别为0m和1m，上下边界为隔水边界。本例的研究区域为各项同性的介质，即x和y方向的渗透系数相等。3个子区域zone 1-zone3的渗透系数分别为：100m/d、10m/d、100m/d。由于本例没有解析解，将使用Method-Zhou-F的解作为标准对照。

Method-Zhou-F将研究区剖分为18432个三角形单元，并应用FEM-F求解水头。DDFEM、Method-Zhou、Method-Yeh将研究区剖分为1152个三角形单元(图5)，应用FEM求解水头。

各方法在交界面AB、CD处的V_x值分别如图6、7所示。图中DDFEM，Method-Zhou在两个交界面上的达西速度不连续，均满足折射定律，且和Method-Zhou-F的值十分接近。Method-Yeh依然只能在交界面处获得一个达西速度，不满足折射定律，精度低于另两种方法。另一方面，本例中DDFEM的CPU时间不到1秒，而Method-Yeh需要1秒，Method-Zhou需要4秒。若DDFEM和Method-Zhou-F一样，使用精细剖分将研究区剖分为18432个三角形单元，则DDFEM需要438s，约为Method-Zhou-F的CPU时间(10119s)的5％。

本例中，zone 1，zone 3的渗透系数等于100m/d，各方法在x＝80m处的V_x值如图8所示；令zone 1，zone 3的渗透系数等于10m/d，V_x值如图9所示；令zone 1，zone 3的渗透系数等于1000m/d，V_x值如图10所示。我们发现在交界面两侧的渗透系数比值上升时，两侧的达西速度的差异会变大。并且Method-Zhou所需对达西速度的迭代次数以及CPU时间也会上升，而DDFEM的CPU时间则保持不变。这一结果显示，DDFEM在交界面两侧的渗透系数差距上升时更加高效。

实施例3：二维稳定流模型，含有交叉的交界面

研究区如图11所示，地下水流和达西速度的方程分别为(6)、(7)式。研究区Ω＝[0,1]×[0,1]。交界面AB(y＝0.5)和CD(x＝0.5)交叉于点o，将研究区分为4个子区域，zone1-zone 4的渗透系数分别为1、10、10、100。本例具有解析解，其中zone 1和zone 3的解析解分别为实施例1的中zone 1和zone 2的解析解。而zone 2和zone 4的解析解水头为：

本例的源汇项、达西速度、第一类边界条件的表达式均可根据渗透系数和解析解水头给出。

本例采用DDFEM、DDFEM-AS、Method-Yeh、Method-Yeh-AS四种方法求解，均把研究区剖分为28800个三角形单元(图11)。DDFEM、Method-Yeh采用FEM计算水头，DDFEM-AS、Method-Yeh-AS采用解析解水头。

和前面两实施例相同，本例中DDFEM、DDFEM-AS在交界面AB、CD上十分精确并满足折射定律。Method-Yeh、Method-Yeh-AS则依然不满足折射定律，在交界面AB上以及AB、CD交叉点o(0.5,0.5)处具有较大误差。

各方法在两界面的交叉点o(0.5,0.5)处的达西速度值V_x如表1所示。我们可以看到DDFEM、DDFEM-AS的解和解析解十分接近，而Method-Yeh、Method-Yeh-AS具有较大误差。我们注意到DDFEM、DDFEM-AS和解析解在zone 1，zone 2的值要小于zone 3，zone 4的。和解析解不同，DDFEM、DDFEM-AS的两个较小的值十分接近，但并不完全一样。这是因为zone 1，zone 4的值是通过解子问题得到的，而zone 2和zone 3的值则是应用折射定律获得的。由于zone 2或zone 3的值既可以通过折射定律和zone 1的值得到，也可以通过折射定律和zone 4的值得到，因此表1中的值为通过这两种方法得到的值的平均值。表1如下：

表1

zone	AS	DDFEM	DDFEM-AS	Method-Yeh	Method-Yeh-AS
						1	4.05259	3.90248	4.08224	21.3105	22.2
2	4.05259	3.89005	4.06594	21.3105	22.2
						3	4.05259	38.9005	40.6594	21.3105	22.2
4	4.05259	38.7762	40.4964	21.3105	22.2

在本例中DDFEM所需的CPU时间仅为1144s，而Method-Yeh则需要2223s，是DDFEM的一倍。由于两种方法均应用FEM计算水头，FEM解水头所需的CPU时间为1070s。因此，DDFEM计算达西速度的时间仅为73s，不到Method-Yeh(1153s)的7％。因此，在水头已知的情况下，DDFEM能够获得比其他方法高得多的计算效率，这一结果证明了，DDFEM不仅能够保证交界面处达西速度的精度，并且具有很高的计算效率。

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

Claims

1.一种模拟地下水介质交界面处达西速度的区域分解有限元法，其特征在于，包括步骤如下：

(9)结合折射定律获得步骤(6)中选取的子问题所在子区域和相邻子区域交界面上另一侧的达西速度值，将该达西速度值作为相邻子区域上的子问题的第一类边界条件；

(10)选取下一个需要求解的子问题，应用上述步骤(3)中的水头获得应用Yeh的伽辽金模型解该子问题时所需的水头条件；若该子问题所在子区域某边界上的达西速度值已通过相邻子问题结合折射定律得到，则将所述达西速度值作为该子问题该边界上的第一类边界条件；

(11)应用步骤(7)-(8)获得步骤(10)中选取的子问题所在子区域上所有的达西速度值，并结合折射定律计算该子区域和相邻子区域交界面上另一侧的达西速度值，将该达西速度值作为相邻子区域上的子问题的第一类边界条件；

(12)重复步骤(10)-(11)，直到步骤(5)中的所有子问题求解完毕。

2.根据权利要求1所述的模拟地下水介质交界面处达西速度的区域分解有限元法，其特征在于，上述的步骤(1)中采用三角形单元剖分研究区域。

3.根据权利要求1所述的模拟地下水介质交界面处达西速度的区域分解有限元法，其特征在于，上述的步骤(2)中，所述的渗透系数K取这个单元的所有顶点上的渗透系数平均值。

4.根据权利要求1所述的模拟地下水介质交界面处达西速度的区域分解有限元法，其特征在于，上述的步骤(7)中，所述的渗透系数K取这个单元的所有顶点上的渗透系数平均值。

5.根据权利要求1所述的模拟地下水介质交界面处达西速度的区域分解有限元法，其特征在于，上述的步骤(2)中，所述的源汇项取这个单元的所有顶点上的源汇项平均值。

6.根据权利要求1所述的模拟地下水介质交界面处达西速度的区域分解有限元法，其特征在于，上述的步骤(7)中，所述的源汇项取这个单元的所有顶点上的源汇项平均值。

7.根据权利要求1所述的模拟地下水介质交界面处达西速度的区域分解有限元法，其特征在于，上述的步骤(11)中，所计算的交界面上另一侧的达西速度值应为之前未进行求解的。