CN112347678B - 一种同时模拟地下水流和达西速度的新型多尺度有限元法 - Google Patents

Abstract

本发明公开了一种同时模拟地下水流和达西速度的新型多尺度有限元法，包括步骤：设定粗、细网格单元的尺度，将研究区剖分为粗网格单元，将粗单元剖分为细网格单元，获得多尺度网格；在粗网格单元上求解退化的椭圆方程，构造基函数；在粗网格单元上求解达西方程，速度矩阵；运用伽辽金法和格林公式得到问题的变分形式，并离散到粗网格上；将达西定律代入各粗网格上变分形式的分量，将水头偏导项转换为速度项，应用速度矩阵将达西速度项应用水头的粗尺度解线性表示，获得该粗网格的单元刚度矩阵，相加得水头的总方程；用有效的矩阵解法求解总方程，同时通过速度矩阵获得达西速度值。与多种经典方法相比，本发明的新型多尺度有限元法具有更高的效率。

Description

一种同时模拟地下水流和达西速度的新型多尺度有限元法

技术领域

本发明属于水力学技术领域，具体涉及一种同时模拟地下水流和达西速度的新型多尺度有限元法。

背景技术

地下水是水资源必不可少的组成部分，具有水质良好、水量稳定等优点，是人类生活与经济发展的十分重要的来源。由于我国生产力布局和地下水资源分布的不匹配，不仅会引起地下水资源的储量减少、水质污染等现象，还会出现海水入侵、地面沉降、土壤恶化、工程地质灾害等地质、环境问题。数值模拟是地下水研究的重要技术手段，它能够预测地下水资源和环境状态及其变化趋势。因此，建立精确高效的水流模型综合分析水头、达西速度等参数具有重要意义。

地下水系统具有非均质性，并常跨越很多尺度，应用传统有限元法需要保证网格内部渗透系数可为常数，直接进行非均质地下水问题水头求解时需要非常大的计算消耗。因此，科研工作者提出了多尺度有限单元法(Hou and Wu 1997)改进了传统有限元法，通过在粗单元上构造多尺度基函数抓住介质的细尺度信息，并直接在粗尺度上求解水头，能够大幅降低计算消耗。

然而，科学工作者主要研究多尺度有限元法的水流方程的解法，而多尺度有限元法达西速度计算方法研究很少。同时，多尺度有限元算法无法保证非均质地下水问题中达西速度的连续性，精度受到了限制。Yeh在1981年于“On the computation of Darcianvelocity and mass balance in the finite element modeling of groundwater flow”一文中提出的伽辽金有限元模型可以解决达西速度连续性的问题，具有比同类方法更高的精度和应用范围。然而，该方法的有限元框架依然会在非均质地下水的模拟时产生大量的计算消耗。另一方面，在模拟地下水综合问题时，水头和达西速度都需要计算，大多数算法需要求解不同的方程进行计算，步骤繁琐，计算消耗高。

发明内容

发明目的：为了克服背景技术的不足，本发明公开了一种同时模拟地下水流和达西速度的新型多尺度有限元法，本方法组合了多尺度有限元法和Yeh的伽辽金有限元方法，通过多尺度基函数提高水头和达西速度的计算效率，通过Yeh的伽辽金有限元方法构造速度矩阵保证达西速度的连续性，通过将该速度矩阵嵌入多尺度有限元法的算法框架实现达西速度和水头的同时计算，能够保证达西速度和水头的质量守恒性，解决现有技术中求解非均质地下水头和连续达西速度的计算效率过低的问题，并实现了求解一个方程即可获得水头和达西速度两项参数。

技术方案：本发明的同时模拟地下水流和达西速度的新型多尺度有限元法，包括以下步骤：

S1、根据研究区域确定所要模拟的地下水问题的定解条件，设定粗网格尺度，对该研究区域进行网格剖分，得到粗单元，粗网格单元的顶点为粗尺度节点；

S2、设定细网格尺度，对上述每一粗网格单元进行网格剖分，得到细单元，细网格单元的顶点为细尺度节点；

S3、在S2中的粗网格的剖分下，根据渗透系数K、粗单元顶点上的多尺度基函数值以及多尺度基函数边界条件公式，在每一粗网格单元上求解退化的椭圆型水流方程以构造多尺度基函数，获得多尺度基函数在每一粗单元上所有细尺度节点上的值；

S4、在S2中的粗网格的剖分下，根据渗透系数K、粗单元上的基函数，结合Yeh的伽辽金有限元理论，在每一粗网格单元上分别求解x及y方向上的达西定律方程，分别获得每一粗网格单元上x及y方向上的速度矩阵；

S5、在研究区域上，运用伽辽金法和格林公式得到所求问题的变分形式，并离散到S1所获的粗网格单元上，获得每个粗网格单元上该变分形式的分量；

S6、在每一粗网格单元上，将x、y方向的达西定律方程代入每个粗网格上变分形式的分量，将水头对x、y的一阶偏导项分别转换为x、y方向的速度项；应用S4中所获的相应粗网格单元上x和y方向速度矩阵，将该达西速度项应用水头的粗尺度解线性表示，获得该粗网格上的单元刚度矩阵；

S7、将S6中所获的每一粗网格单元上的单元刚矩阵相加，得到研究区总刚度矩阵；根据研究区域的边界条件、源汇项，计算右端项，形成水头的粗尺度解的方程，采用cholesky分解法，求得研究区域上每个节点的水头；

S8、在每个粗单元上，应用S4中构造的每一粗网格单元的x和y方向速度矩阵，获得相应粗网格单元上的细尺度达西速度；

S9、在S1中的研究区剖分的网格线上，在该网格线上的S1中研究区剖分所获的每一粗尺度节点上，分别平均从该节点的周围的粗单元中按照S8获得的该节点的x、y方向细尺度达西速度值，来获得该节点的粗尺度x、y方向的达西速度值；

S10、在S1中的研究区剖分的网格线上，在该网格线上的S2中的粗单元剖分所获的细尺度节点上，分别平均从该节点的相关粗网格中按照S8获得的该节点的x、y方向细尺度达西速度值，来获得该节点的x、y方向粗尺度达西速度值。

其中，S1中采用矩形单元网格剖分研究区域，以形成矩形粗网格单元。

进一步的，S2中采用直角三角形单元网格剖分粗网格单元，以形成直角三角形细网格单元。

进一步的，S10中粗网格线上的每一细尺度节点不包括S9中已经计算过的节点。

有益效果：与现有技术相比，本发明的优点为：首先，本发明求解一次矩阵即可获得地下水水头和x、y方向上的粗、细尺度达西速度，不需要额外的计算消耗计算达西速度；其次，可以保证x、y方向上达西速度的连续性和守恒性；同时，本发明能够获得和精细剖分的有限元法、多尺度有限单元法同阶的地下水水头精度，计算时间和多尺度有限元法一致，远小于精细剖分的有限元法；再而，本发明能够获得和精细剖分的Yeh的伽辽金有限元方法相近的地下水达西速度的精度，但计算消耗远小于该方法；最后，本发明能够精确高效求解地下水稳定流和非稳定流问题，可以有效处理多种非均质介质。

附图说明

图1为本发明新型多尺度有限元法的研究区剖分示意图；

图2为本发明新型多尺度有限元法的粗单元剖分示意图；

图3为二维连续介质稳定流模型，各方法在y＝0.7截面处的水头绝对误差示意图；

图4为二维连续介质稳定流模型，各方法在y＝0.7截面处的粗尺度达西速度绝对误差示意图；

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

图6为二维渐变介质稳定流模型，各方法在y＝6000m截面处的粗尺度x方向达西速度的绝对误差示意图；

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

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

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

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

一种同时模拟地下水流和达西速度的新型多尺度有限元法，包括以下步骤：

S1、根据研究区域确定所要模拟的地下水问题的定解条件，设定粗网格尺度，如图1所示，对该研究区域进行网格剖分，得到粗单元，粗网格单元的顶点为粗尺度节点；

S2、设定细网格尺度，对上述每一粗网格单元进行网格剖分，得到细单元，细网格单元的顶点为细尺度节点；

S3、在S2中的粗网格的剖分下，根据渗透系数K、粗单元顶点上的多尺度基函数值以及多尺度基函数边界条件公式，在每一粗网格单元上求解退化的椭圆型水流方程以构造多尺度基函数，获得多尺度基函数在每一粗单元上所有细尺度节点上的值；

S4、在S2中的粗网格的剖分下，根据渗透系数K、粗单元上的基函数，结合Yeh的伽辽金有限元理论，在每一粗网格单元上求解x及y方向上的达西定律方程，分别获得每一粗网格单元上x及y方向上的速度矩阵；

S5、在研究区域上，运用伽辽金法和格林公式得到所求问题的变分形式，并离散到S1所获的粗网格单元上，获得每个粗网格单元上该变分形式的分量；

S6、在每一粗网格单元上，将x、y方向的达西定律方程代入每个粗网格上变分形式的分量，将水头对x、y的一阶偏导项分别转换为x、y方向的达西速度项；应用S4中所获的相应粗网格单元上x和y方向速度矩阵，将该达西速度项应用水头的粗尺度解线性表示，获得该粗网格上的单元刚度矩阵；

S7、将S6中所获的每一粗网格单元上的单元刚矩阵相加，得到研究区总刚度矩阵；根据研究区域的边界条件、源汇项，计算右端项，形成水头的粗尺度解的方程，采用cholesky分解法，求得研究区域上每个节点的水头；

S8、在每个粗单元上，应用S4中构造的每一粗网格单元的x和y方向速度矩阵，获得相应粗网格单元上的细尺度达西速度；

S9、在S1中的研究区剖分的网格线上，在该网格线上的S1中研究区剖分所获的每一粗尺度节点上，分别平均从该节点的周围的粗单元(图1中e₁-e₄)中按照S8获得的该节点的x、y方向细尺度达西速度值，来获得该节点的粗尺度x、y达西速度值；

S10、在S1中的研究区剖分的网格线上，在该网格线上的S2中的粗单元剖分所获的细尺度节点上，分别平均从该节点的相关粗网格中按照S8获得的该节点的x、y方向细尺度达西速度值，来获得该节点的x、y方向粗尺度达西速度值。

其中，S1中采用矩形单元网格剖分研究区域，以形成矩形粗网格单元。

S2中采用直角三角形单元网格剖分粗网格单元，以形成直角三角形细网格单元。

S10中粗网格线上的每一细尺度节点不包括S9中已经计算过的节点。

以研究区域Ω为例，S1中的网格剖分为按照图1的细实线对研究区进行的剖分，以获得粗单元；S2中的网格剖分为按照图2的对每一粗单元细分，以获得细单元。

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

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

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

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

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

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

据多尺度有限元法(Hou and Wu 1997)的工作，在粗单元□_ijkl上

H(x,y)＝H_iΨ_i(x,y)+H_jΨ_j(x,y)+H_kΨ_k(x,y)+H_lΨ_l(x,y) (4)

其中Ψ_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 (5)

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

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

[A^x][V_x]＝[B^x][H] (6)

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

[V_x]＝[A^x]^-1[B^x][H] (7)

y方向的速度矩阵的构造过程与x方向的速度矩阵构造过程类似。

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

经过S5，应用伽辽金变换和格林公式，可得：

S6，将达西定律方程代入(9)，可得

将速度矩阵带入式(10)，可以得到□_ijkl上关于水头粗尺度解的单元刚度矩阵。S7，联立研究区上所有粗网格的单元刚度矩阵，即可得到研究区上关于水头粗尺度解的方程组，求解后可以得到研究区上粗尺度网格上每个节点的水头值。在每一粗单元上，通过基函数插值可以得到细尺度水头。S8，在每一粗单元上，根据粗尺度水头和速度矩阵可以获得细尺度的达西速度。再根据S9、S10即可以得到粗尺度达西速度。

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

LFEM:有限单元法；

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

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

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

NMSFEM-O:新型多尺度有限元法,振荡边界条件；

NMSFEM-L:新型多尺度有限元法，线性边界条件；

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

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

实施例1

二维连续介质稳定流模型

稳定流方程为：

研究区域Ω＝[0,1]×[0,1]，含水层厚度为1，考虑两个不同的渗透系数场，K＝K₁＝1和

狄利克雷边界条件和源汇项W由解析水头H＝xy(1-x)(1-y)给出。

本例采用FEM-F，NMSFEM，MSFEM和FEM求解水头场。NMSFEM，MSFEM和FEM将研究区Ω的每个边界划分为30等份，NMSFEM和MSFEM将研究区Ω分为900(30×30)个矩形粗单元，每个粗单元又分为128(8×8×2)个三角形细单元；FEM将研究区Ω分为1800(30×30×2)个三角形细单元；FEM-F将研究区Ω分为240×240×2个单元，以保证其单元总数与NMSFEM和MSFEM的单元总数相同。

如图3所示，在两个不同渗透系数场中沿y＝0.7截面由FEM-F，NMSFEM-L，MSFEM-L和FEM求解的水头绝对误差。可以看到，两个场中的FEM-F，NMSFEM-L，MSFEM-L曲线几乎重合，说明它们的误差很接近，而FEM曲线明显高于其它曲线，其误差最大。子图(a)显示，NMSFEM-L相较于MSFEM-L和FEM-F更精确；子图(b)表明，由于非均匀介质的影响，这几种方法的误差相较于子图(a)更大。综上，FEM-F，NMSFEM和MSFEM的误差变化范围比FEM小，更能反映它们在解决渗透系数为非均质情况下的高效性。在CPU耗时方面，以K₁为例，FEM-F求解57121个网格需57779s，MSFEM-L求解841个粗尺度网格和57121个细尺度网格需6s，而NMSFEM-L同样花6s可以求解841个粗尺度网格，57121个细尺度网格，以及961×2个粗尺度速度和58081×2个细尺度速度，说明了在获得相似精度的情况下，NMSFEM和MSFEM的计算效率比FEM-F更高，并且能在一个计算过程中对粗/细尺度水头、x,y粗/细尺度速度进行计算。

如图4所示，在两个不同渗透系数场中沿y＝0.7截面由Method-Yeh-F,NMSFEM-L和Method-Yeh求解的粗尺度速度绝对误差。由图可知，NMSFEM-L和Method-Yeh-F的误差接近，且明显比Method-Yeh小；同时，渗透系数对速度的影响远大于对水头的影响。当渗透系数场从K₁变为K₂时，速度的绝对误差相较于水头增加得更多；当K＝K₂时，NMSFEM-L曲线和Method-Yeh-F曲线的振幅变化比Method-Yeh小，说明这两种方法具有较强的抗非均质影响的能力。在CPU耗时方面，NMSFEM-L不花额外的计算消耗求解两个方向上的达西速度，求解相同数量未知量的时间仅为Method-Yeh-F的0.00331816％。

实施例2

二维渐变介质稳定流模型

研究区为一正方形区域：Ω＝[0,10km]×[0,10km]，研究方程为方程(11)，研究区内的渗透系数从左到右满足关系：

含水层厚度为1m。

本例中研究区的左右边界为定水头边界，分别为10m和0m，上下边界为隔水边界。

本例采用FEM-F，NMSFEM-O，NMSFEM-L，MSFEM-O，MSFEM-L和FEM求解水头场。NMSFEM，MSFEM和FEM将研究区Ω的每一边界分为20等份，NMSFEM和MSFEM将研究区Ω分为400(20×20)个矩形粗单元，每个粗单元分为72(6×6×2)个三角形细单元；FEM将研究区Ω分为800(20×20×2)个三角形细单元；FEM-F将研究区Ω分为120×120×2个网格单元，即和NMSFEM，MSFEM同样数目的细单元。由于本实例没有解析解，因此将精细部分的FEM-F的解近似作为解析解。

图5是FEM-F，NMSFEM-O，NMSFEM-L，MSFEM-O，MSFEM-L，FEM在y＝6000m截面求解的水头值。由图可知，NMSFEM-O结果比FEM-F更精确，表明了NMSFEM方法的有效性，FEM-F，MSFEM-O，NMSFEM-L，MSFEM-L的精确性依次递减。NMSFEM-L解优于MSFEM-L解，FEM解最差，这与实施例1的结果一致。在计算时间方面，NMSFEM和MSFEM求解399个粗尺度网格和14399个细尺度网格需2s，而FEM-F求解与NMSFEM相同的14399个细尺度网格需要928s，和上例一样，NMSFEM能在更短时间求解与其它方法一样数量的粗/细尺度网格，计算效率更高。

图6是NMSFEM-O，NMSFEM-L，Method-Yeh-F和Method-Yeh在y＝6000m截面求解的速度的绝对误差。与水头结果相似，NMSFEM-O速度比Method-Yeh-F速度更精确，NMSFEM-O结果明显好于NMSFEM-L，表明振荡边界条件也可以帮助NMSFEM提高速度的精度，Method-Yeh表现最差。与图5比较可以看出，速度的精度与水头的精度有关，水头越精确，速度也越精确。在计算时间方面，NMSFEM-O可以用更少的计算量获得与Method-Yeh-F相近的结果，计算效率明显更高。

实施例3

二维渐变介质非稳定流模型

研究方程为：

研究区为：Ω＝[0,10km]×[0,10km]，研究区内的渗透系数从左到右满足关系：

左右边界为定水头边界，分别为10m和0m，上下边界为隔水边界。单位储水量为10^-6/m，含水层厚度为1m。在点(6000m,6000m)处有一抽水井，恒定流速为100m/day，总泵送时间为3天，时间步长为1天，本例的初始水头为H₀＝0。

本例采用FEM-F，NMSFEM-O，MSFEM-O和FEM求解水头场。NMSFEM，MSFEM和FEM将研究区Ω分为20等份，NMSFEM和MSFEM将研究区Ω分为400(20×20)个矩形粗单元，每个粗单元分为128(8×8×2)个三角形细单元；FEM将研究区Ω分为800(20×20×2)个三角形细单元；FEM-F将研究区Ω分为160×160×2个单元，以保证其单元总数与NMSFEM和MSFEM的单元总数相同。由于本实例没有解析解，因此将精细部分的FEM-F的解近似作为解析解。

图7是“True solution”、FEM-F，NMSFEM-O，MSFEM-O和FEM在y＝6000m截面处的水头值。FEM-F，NMSFEM-O和MSFEM-O的精确度相近，NMSFEM-O的精确度优于MSFEM-O，FEM表现最差。由于抽水井的影响，在点(6000m，6000m)附近的NMSFEM-O和MSFEM-O结果比在其它点处的结果有更大的误差。抽水井附近的水头变化很快，以致于NMSFEM-O和MSFEM-O很难对其精确求解。然而，更精细的网格可以改善抽水井附近的结果，所以FEM-F在抽水井处的精确度最好。另一方面，在远离井的区域，NMSFEM-O的解优于FEM-F和MSFEM-O的解，表明了NMSFEM方法的有效性。在CPU耗时方面，NMSFEM和MSFEM用7s求解了399个粗尺度水头和25599个细尺度水头，而FEM-F求解25599个水头用了15693s，表明在非稳定流动问题中，NMSFEM和MSFEM比在稳定流动问题中节省更多的计算量。

图8是Method-Yeh-F，NMSFEM-O，Method-Yeh和“True solution”在y＝6000m处求解的粗尺度达西速度。由图可知，NMSFEM-O与Method-Yeh-F有相似的精确度，NMSFEM-O比Method-Yeh的精确度高。在井的附近(x＝5500～6500m)，NMSFEM-O的精确度比Method-Yeh-F低，但仍比Method-Yeh高，在其它点处NMSFEM-O的精确度比Method-Yeh-F高。

图9是Method-Yeh-F，NMSFEM-O和“True Solution”在y＝6125m处求解的细尺度达西速度。NMSFEM-O和Method-Yeh-F都有较高的精确度；在抽水井附近，NMSFEM-O有较大的误差，而Method-Yeh-F仍有较高的精确度。在CPU耗时方面，NMSFEM-O用时7s计算粗、细尺度水头和x,y方向的达西速度，而Method-Yeh-F用时21213s计算水头和x方向的达西速度，NMSFEM-O求解相同数量的水头和x方向的达西速度的时间仅为Method-Yeh-F的0.033％，因此NMSFEM-O有更高的计算效率。

Claims (1)

1.一种同时模拟地下水流和达西速度的新型多尺度有限元法，其特征在于，包括以下步骤：

S1、根据研究区域确定所要模拟的地下水问题的定解条件，设定粗网格尺度，对该研究区域进行网格剖分，得到粗单元，粗网格单元的顶点为粗尺度节点；

S2、设定细网格尺度，对上述每一粗网格单元进行网格剖分，得到细单元，细网格单元的顶点为细尺度节点；

S3、在S2中的粗网格的剖分下，根据渗透系数K、粗单元顶点上的多尺度基函数值以及多尺度基函数边界条件公式，在每一粗网格单元上求解退化的椭圆型水流方程以构造多尺度基函数，获得多尺度基函数在每一粗单元上所有细尺度节点上的值；

S4、在S2中的粗网格的剖分下，根据渗透系数K、粗单元上的基函数，结合Yeh的伽辽金有限元理论，在每一粗网格单元上求解x及y方向上的达西定律方程，分别获得每一粗网格单元上x及y方向上的速度矩阵；

S5、在研究区域上，运用伽辽金法和格林公式得到所求问题的变分形式，并离散到S1所获的粗网格单元上，获得每个粗网格单元上该变分形式的分量；

S6、在每一粗网格单元上，将x、y方向的达西定律方程代入每个粗网格上变分形式的分量，将水头对x、y的一阶偏导项分别转换为x、y方向的达西速度项；应用S4中所获的相应粗网格单元上x和y方向速度矩阵，将该达西速度项应用水头的粗尺度解线性表示，获得该粗网格上的单元刚度矩阵；

S7、将S6中所获的每一粗网格单元上的单元刚矩阵相加，得到研究区总刚度矩阵；根据研究区域的边界条件、源汇项，计算右端项，形成水头的粗尺度解的方程，采用cholesky分解法，求得研究区域上每个节点的水头；

S8、在每个粗单元上，应用S4中构造的每一粗网格单元的x和y方向速度矩阵，获得相应粗网格单元上的细尺度达西速度；

S9、在S1中的研究区剖分的网格线上，在该网格线上的S1中研究区剖分所获的每一粗尺度节点上，分别平均从该节点的周围的粗单元中按照S8获得的该节点的x、y方向细尺度达西速度值，来获得该节点的粗尺度x、y方向的达西速度值；

S10、在S1中的研究区剖分的网格线上，在该网格线上的S2中的粗单元剖分所获的细尺度节点上，分别平均从该节点的相关粗网格中按照S8获得的该节点的x、y方向细尺度达西速度值，来获得该节点的x、y方向粗尺度达西速度值；

S1中采用矩形单元网格剖分研究区域，以形成矩形粗网格单元；S2中采用直角三角形单元网格剖分粗网格单元，以形成直角三角形细网格单元；

S10中粗网格线上的每一细尺度节点不包括S9中已经计算过的节点。

Images