模拟地下水溶质运移的新型有限体积多尺度有限元方法与流程

本发明属于水力学技术领域，具体涉及一种模拟地下水溶质运移的新型有限体积多尺度有限元方法(nfvmsfem)。

背景技术：

地下水是地球上最主要、分布最广泛的天然水资源之一，是人类主要引用水来源。地下水资源的滥用和过度利用，地下水污染日益严峻。对流-弥散方程被广泛的应用于描述地下水溶质运移模拟问题。然而，在模拟中，对流项和弥散项的结合导致了数值频散和振荡，经典的有限元方法和有限差分法都会产生较大的误差。同时，在对流占优的情况，有限元法需要精细的网格来减少数值离散和振荡，这导致了较高的计算成本。此外，潜在的大尺度非均质地下水问题可能使这种情况更差，并需要更多的成本。最后，在经典的有限元方法中，简单的差分格式不足以处理瞬态问题中的时间项，精细的时间步长需要更大的计算消耗。另一方面，弥散速度在评价弥散项对溶质运移过程的贡献方面起着至关重要的作用，能够同时计算浓缩速度和弥散速度的模型比较少见，而分别计算这两个未知量需要较高的计算成本。为了解决上述问题，本发明提出了对流-弥散方程的新型有限体积多尺度有限元方法。新型有限体积多尺度有限元方法是基于多尺度有限单元法(houandwu1997)、有限体积法和yeh有限元方法，应用了有限体积多尺度有限元法的框架。

多尺度有限单元法(houandwu1997)是一种基于有限单元法的高效方法，其通过在粗网格单元上构造多尺度基函数抓住介质的细尺度信息，并直接在粗尺度上求解水头，能够大幅降低计算消耗。随后，科学工作者在“finitevolumemultiscalefiniteelementmethodforsolvingthegroundwaterflowproblemsinheterogeneousporousmedia”一文中提出了有限体积多尺度有限元法(heandren2005)，将有限体积法和多尺度有限单元法有机结合，可以通过有限体积理论保证水流的局部质量守恒，具有比多尺度有限单元法更高的计算精度。然而，科学工作者只给出了有限体积多尺度有限元法的水流方程的解法，如何应用有限体积多尺度有限元法计算溶质运移还有待研究，特别是上述解溶质运移问题的相关问题无法解决。

另一方面，大部分算法无法保证非均质地下水问题中达西速度的连续性，精度受到了限制。yeh有限元方法(“onthecomputationofdarcianvelocityandmassbalanceinthefiniteelementmodelingofgroundwaterflow”)能够解决达西连续性问题，但没有溶质运移问题中弥散速度相应算法。同时，由于其有限元框架的限制，该方法效率较低。

技术实现要素：

针对于上述现有技术的不足，本发明的目的在于提供一种模拟地下水溶质运移的新型有限体积多尺度有限元方法，解决现有技术中求解溶质运移问题的计算效率过低、求解对流占优情况下溶质运移问题的精度较低及无法估计弥散速度并保证其连续性的问题。本发明方法针对溶质运移问题进行了改进，将能够高效精确的解对流扩散方程，同时获得浓度和弥散速度两种未知量，并可以对流占优情况下获得精确的解。

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

本发明的一种模拟地下水溶质运移的新型有限体积多尺度有限元方法，步骤如下：

(1)根据所要模拟的地下水溶质运移问题确定研究区域的边界条件，设定粗、细网格单元的尺度，将研究区域剖分为粗网格单元，将每一粗网格单元剖分为细网格单元，形成多尺度网格；

(2)以步骤(1)中多尺度网格的粗网格单元上的每一未知节点为有限体积网格中的矩形控制体积的基点，连接各有限体积网格中的矩形控制体积的基点周围粗网格单元的中心点，形成以各基点为中心的矩形控制体积，由各个矩形控制体积构成有限体积网格；

(3)在步骤(1)中多尺度网格的每一粗网格单元上，根据弥散系数d、多尺度基函数的边界条件公式，求解基于弥散系数的退化椭圆方程，获得多尺度基函数；

(4)定义任一矩形控制体积的基点处的浓度粗尺度解为：浓度在该矩形控制体积上的积分除以该矩形控制体积的面积；

(5)在步骤(1)中多尺度网格的每一粗网格单元上，根据弥散系数d、多尺度基函数，及x方向上的fick定律方程，应用yeh有限元方法，将弥散速度项放置在根据fick定律方程形成的方程组左端，将浓度项放置在该方程组右端，得到关于弥散速度的方程组，将弥散速度的系数矩阵求逆，并与方程组右端的浓度项的系数向量相乘，得到x方向弥散速度矩阵；

(6)在步骤(1)中多尺度网格的每一粗网格单元上，根据弥散系数d、多尺度基函数，及y方向上的fick定律方程，应用yeh有限元方法，将弥散速度项放置在根据fick定律方程形成的方程组左端，将浓度项放置在该方程组右端，得到关于弥散速度的方程组，将弥散速度的系数矩阵求逆，并与方程组右端的浓度项的系数向量相乘，从而得到y方向弥散速度矩阵；

(7)在步骤(2)中有限体积网格的每一个矩形控制体积上将描述溶质运移问题的对流扩散方程积分，将浓度粗尺度解对时间的偏微分代入对流扩散方程，应用散度定理变换对流扩散方程；

(8)将步骤(7)中经过散度定理变换的对流扩散方程离散到与构成该矩形控制体积的相关粗网格单元上，应用相关粗网格单元的多尺度基函数、弥散速度矩阵将弥散项、对流项离散，应用crank-nicolson格式处理时间项，获得该矩形控制体积关于浓度粗尺度解的方程的具体细尺度形式；

(9)联立研究区域内所有矩形控制体积的关于浓度粗尺度解的方程，应用qr分解法进行求解，获得研究区域内各个节点的浓度粗尺度解；

(10)基于每个粗网格单元上的x，y方向的弥散速度矩阵分别获得该粗网格单元上的x，y方向的细尺度弥散速度；

(11)在步骤(1)中多尺度网格的节点上，分别平均该节点的x，y方向的细尺度弥散速度以分别获得该节点的x，y方向的粗尺度弥散速度。

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

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

优选的，所述步骤(3)中构造多尺度基函数具体包括：以与i点相关的多尺度基函数ψi为例，在矩形粗网格单元□ijkl上应用有限单元法求解如下基于弥散系数的退化椭圆方程：

其中，d为渗透系数，ψi的边界条件可以选择振荡边界条件或线性边界条件。

优选的，所述步骤(4)中定义浓度的粗尺度解具体为：以i点为例，在i的浓度粗尺度解为浓度c在矩形控制体积ii积分并除以ii的面积sii：

优选的，所述步骤(8)中矩形控制体积的相关粗网格单元为包含该矩形控制体积的最小粗网格单元子集，即上述步骤(2)中有限体积网格中的矩形控制体积的基点周围的粗网格单元。

本发明的有益效果：

1、本发明的方法，能够有效降低对流占优引发的数值弥散和振荡。

2、本发明在解决非稳定溶质运移问题时能够在对流占优条件下使用较大的时间步长，从而显著降低计算消耗。

3、本发明在模拟溶质运移问题时能够在一次计算过程中获得浓度和x，y方向的细尺度弥散速度，无需额外计算成本获得弥散速度。

4、本发明可以保证弥散速度的连续性和守恒性，从而获得较高的计算精度。

5、本发明能够获得和精细剖分的有限元法相近的浓度精度，但计算成本远小于该方法。

6、本发明在对流占优的条件下的浓度精度、计算效率高于同条件下的多尺度有限单元法，但计算成本几乎一致。

7、本发明能够获得和精细剖分的yeh的有限元模型相近的弥散速度的精度，但仅使用远低于其的计算成本。

8、本发明能够高效求解稳定和非稳定溶质运移问题，可以有效处理多种变化状态下的弥散系数。

附图说明

图1为本发明的新型有限体积多尺度有限元法的研究区剖分示意图；其中，实线为多尺度有限元网格，虚线为有限体积法网格。

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

图3为本发明以i为基点生成的矩形控制体积ii(□abcd)及其相关的粗网格单元eⅰeⅱeⅲeⅳ示意图。

图4为二维稳定流溶质运移模型，各方法在情形一不同d时的粗尺度浓度平均相对误差值示意图。

图5a为二维稳定流溶质运移模型，method-yeh-f方法在情形二计算的粗尺度弥散速度绝对误差的示意图。

图5b为二维稳定流溶质运移模型，nfvmsfem方法在情形二计算的粗尺度弥散速度绝对误差的示意图。

图5c为二维稳定流溶质运移模型，method-yeh方法在情形二计算的粗尺度弥散速度绝对误差的示意图。

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

具体实施方式

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

本发明的一种模拟地下水溶质运移的新型有限体积多尺度有限元方法(nfvmsfem)，步骤如下：

步骤(1)：根据所要模拟的地下水溶质运移问题确定研究区域的边界条件，设定粗、细网格单元的尺度，将研究区域剖分为粗网格单元，将每一粗网格单元剖分为细网格单元，形成多尺度网格；

步骤(2)：以步骤(1)中多尺度网格的粗网格单元上的每一未知节点为有限体积网格中的矩形控制体积的基点，连接各有限体积网格中的矩形控制体积的基点周围粗网格单元的中心点，形成以各基点为中心的矩形控制体积，由各个矩形控制体积构成有限体积网格；

设研究区ω为矩形区域(如图1所示)，如步骤(1)所述，按照图1中的实线将研究区剖分为粗网格单元，按照图2将每一个粗网格单元剖分为细网格单元(图2)，从而形成多尺度网格；如步骤(2)所述，以多尺度网格的粗网格单元上的某一未知节点i为例，以i点为矩形控制体积的基点，如图3所示，连接其周围粗网格单元eⅰeⅱeⅲeⅳ的中心a、b、c、d，获得矩形控制体积ii(□abcd)，图2虚线为研究区ω的有限体积网格。

步骤(3)：构造多尺度基函数；以与i点相关的多尺度基函数ψi为例，在矩形粗网格单元□ijkl(图2)上应用有限单元法求解如下基于弥散系数的退化椭圆方程：

其中，d为渗透系数，ψi的边界条件可以选择振荡边界条件或线性边界条件，本示例中选择振荡边界条件。

求解过程为：设m＝1，2，...，p为□ijkl的内点，在式(1)左右乘以与m＝1，2，...，p相关的线性基函数，再进行伽辽金有限单元变分，可以得到式(1)的变分形式；在此基础上，将该变分形式离散到各个细网格单元(图2中的三角形单元，如△abc)上，并通过线性基函数将ψi展开，可以获得关于ψi的关于未知节点m＝1，2，...，p的方程组，应用cholesky分解法求解该方程组，即可获得ψi在矩形粗网格单元□ijkl上的值。

对于研究区每个粗网格单元的每个顶点的相关多尺度基函数进行上述求解过程，即可获得所有多尺度基函数的值。

步骤(4)：定义浓度的粗尺度解；以i点为例，在i的浓度粗尺度解为浓度c在矩形控制体积ii(□abcd)积分并除以ii的面积sii：

步骤(5)：构造x、y方向的弥散速度矩阵；设h＝x,y，在矩形粗网格单元□ijkl上考虑fick定律方程：

式中，vdh为在方向h的弥散速度；dh为在方向h的弥散系数。

应用有限单元法求解式(3)，应用矩形粗网格单元□ijkl所有节点相关的线性基函数：nτ,τ＝1,2,..，nf分别乘以达西定律方程的两端，积分得：

式中，nf为□ijkl的总结点数目；

在矩形粗网格单元□ijkl上浓度被多尺度基函数表示为：

c＝φiψi+φjψj+φkψk+φlψl(5)

基于yeh的有限元模型(yeh1981)，在矩形粗网格单元□ijkl的任意子单元δabc上，有：

vdh＝vdh(a)na+vdh(b)nb+vdh(c)nc(6)

根据上式(4)、(5)、(6)得到关于弥散速度vdh的方程组：

[α^h][vdh]＝[β^h][φ^c](7)

将弥散速度的系数矩阵求逆，并与右端浓度项的系数向量相乘，得到弥散矩阵，即：

[γ^h]＝[α^h][β^h](8)

式中，[γ^h]为h方向的弥散系数矩阵。

步骤(6)：以描述非稳定溶质运移的对流扩散方程为例：

式中，ux、uy分别为x,y方向上的水流速度，nw为源汇项。

将上式在每个矩形控制体积上积分，并应用散度定理；以矩形控制体积ii上的求解过程为例：

结合fick定律，可以得到：

步骤(7)：将式(5)、(6)、(8)代入式(11)，并离散到与构成矩形控制体积ii的相关粗网格单元上，得到式(11)右端的详细形式，记为mi；

设时间步长为dt，记tk，tk+1时刻的mi为mi(tk+1)，应用crank-nicolson格式，可以得到ii上的关于浓度粗尺度解的方程：

所述步骤(7)中矩形控制体积的相关粗网格单元为包含该矩形控制体积的最小粗网格单元子集，即上述步骤(2)中有限体积网格中的矩形控制体积的基点周围的粗网格单元。如图2所示eⅰ、eⅱ、eⅲ、eⅳ为矩形控制体积□abcd相关的粗网格单元。

步骤(8)：联立研究区域上所有矩形控制体积上的方程(12)，可以得到关于浓度粗尺度解的总方程，应用qr分解法求解，得到浓度粗尺度解在研究区任意节点上的值；

步骤(9)：根据式(8)可以得到细尺度弥散速度vdh。

步骤(10)：在步骤(1)中多尺度网格的节点上，分别平均该节点的x,y方向的细尺度弥散速度以分别获得该节点的x,y方向的粗尺度弥散速度，以i点为例：

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

lfem：有限单元法；

lfem-f：精细剖分的有限单元法；

method-yeh：yeh的有限单元模型；

method-yeh-f：精细剖分的yeh的有限单元模型；

nfvmsfem：新型有限体积多尺度有限元法；

peclet数：

courant数：

式中，u为水流速度，l为单元尺度，d为弥散系数，dt为时间步长。

实施例1：二维稳定流溶质运移模型：

研究区域ω＝[0,1]]x[0,1]，ux是x方向的水流速度，dx＝dy＝d；式(13)的第一类边界条件和源汇项由解析解c＝xy-xy²-x²y+x²y²给出。

采用nfvmsfem与msfem、lfem-f和lfem进行浓度求解；nfvmsfem、msfem和lfem将ω的每个边划分为相同的数字(nc)等份，nfvmsfem和msfem再将每个粗网格单元划分为32个三角形细网格单元；lfem将ω划分为2nc²三角形单元。lfem-f将ω划分为32nc²个三角形单元。nfvmsfem和msfem有两层网格剖分，因此有两种peclet数。基于粗网格单元尺度和细网格单元尺度的peclet数分别定义为pec和pef。根据上面的网格划分，可以得到一些关系：pec＝pelfem，pef＝pelfem-f和pef＝0.25pec。

nfvmsfem能够同时获得弥散速度，本例采nfvmsfem、method-yeh-f和method-yeh进行弥散速度求解。method-yeh-f、method-yeh需要lfem-f、lfem来获得浓度分布，因而他们的剖分分别与lfem-f、lfem相同。

情形一：

令d＝0.001，0.0001，0.00001，0.000001，u_x＝1，nc＝30，来测试对流占优时nfvmfem处理高peclet数的能力，相应的pec分别为33.3，333.3，3333.3，33333.3。

基于不同d的lfem-f、nfvmsfem、msfem和lfem的在ω的平均相对误差如图4所示。由图可知，各个方法的计算精度从高到低依次为lfem-f、nfvmsfem、msfem和lfem。同时，nfvmsfem的曲线几乎是一条水平直线，而其他方法的曲线均受到d变化的影响。当d变小(pec升高)时，除了nfvmsfem的其他方法误差显著变大。总之，算例结果表明，nfvmsfem具有处理高peclet数的优点，比lfem-f、msfem和lfem更适合于对流占优的情况。

情形二：

设d＝0.1，ux＝0.0001，nc＝30。这是一个弥散占优的情况，来测试nfvmsfem弥散速度的模拟情况。本例将nfvmsfem弥散速度通过与method-yeh-f、method-yeh的计算结果进行比较，来验证该方法计算弥散速度时的精度和效率。

如上所述，method-yeh-f、method-yeh采用两个分离的计算过程依次计算浓度和弥散速度，而nfvmsfem可以在单次计算中获得浓度和x、y方向的弥散速度及其精细值，且计算消耗可以忽略不计，十分高效与方便。

method-yeh-f、method-yeh的浓度分别由lfem-f和lfem求解；方法method-yeh-f、nfvmsfem和方法method-yeh的浓度平均相对误差分别为0.0033％、0.0794％和0.1313％。method-yeh-f方法的cpu时间为3208s，nfvmsfem方法的cpu时间为6s。

图5显示了用nfvmsfem、method-yeh-f、method-yeh方法计算的粗尺度弥散速度vdx的绝对误差。从图中可以看出，nfvmsfem的精度与method-yeh-f方法非常接近，且远小于method-yeh方法。同时，nfvmsfem不需要额外的弥散速度成本，其cpu时间仍保持在6s，method-yeh-f方法在弥散速度计算中需要额外的7879s的cpu时间，其总cpu时间为11087s，method-yeh方法的cpu时间为浓度c(3s)和弥散速度(2s)总共5s，略低于nfvmsfem方法(6s)。然而，method-yeh的c和弥散速度解分别比nfvmsfem的精度低得多。总之，nfvmsfem能够比method-yeh-f更高效的获取相近精度弥散速度，且精度远高于method-yeh。

实施例2：二维非稳定流溶质运移模型

研究方程为方程(9)，研究区域ω＝[0,1]]x[0,1]，ux＝1是x方向的水流速度，uy＝0，dx＝dy＝d,时间步长为dt,时间步数为ndt，总模拟时间ts＝dt*ndt。此问题的dirichlet边界条件和源/汇项由解析解c＝(xy-xy²-x²y+x²y²)e^-t和每种情况下的各个参数值指定。初始浓度分布遵循以下函数：c＝(xy-xy²-x²y+x²y²)。

nfvmsfem、msfem和lfem将ω的每个边划分为相同的数字(nc)等分，各方法的详细划分与在例1中使用的划分相同。对于粗网格单元和细网格单元，nfvmsfem和msfem的courant数分别定义为crc和crf(crf＝4crc)。

设d＝0.02，nc＝40，ts＝0.2，ndt＝2、12、16、25、32、40，相应的crc(crlfem)为4,0.666，0.5，0.32，0.25，0.2。lfem的peclet数和pec相同，为1.25和pef＝0.25pec。

基于不同courant数的nfvmsfem、msfem和lfem的ω浓度平均相对误差如图6所示。图中显示的都是在ts＝0.2的结果。从左到右，courant数从0.2增加到4，而ndt从40减少到2。nfvmsfem获得了最佳精度。同时，nfvmsfem的曲线几乎与坐标轴平行，表明courant数的影响不大。msfem的精度次之，但受courant数变化的影响较大。这一结果说明，当courant数较小时，多尺度网格可以帮助其获得较高的精度；当courant数较大时，多尺度细网格难以抵抗courant数的增长。lfem的精度最差，而且受courant数的影响也很大。

虽然高ndt可以降低courant数，但它也意味着更多的时间步数和更多的cpu时间。一般来说，地下水数值模拟只需要在特定时间点的结果。如图6所示，nfvmsfem在较低时间步数，较大时间步长的时候获得的精度远高于其他方法，说明了nfvmsfem在模拟时更具计算效率。

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