CN103778298A - 改进的模拟多孔介质中二维水流运动的多尺度有限元方法 - Google Patents

Abstract

本发明公开了一种改进的模拟多孔介质中二维水流运动的多尺度有限元方法。该方法首先将需要求解的问题转换成变分形式；确定边界条件，设定网格单元尺度h，剖分研究区域，得到粗网格单元；对每一粗网格单元进行细剖分；根据渗透系数K以及基函数的边界条件，求解退化的椭圆型问题，确定基函数；根据基函数得到单元刚度矩阵，相加得到总刚度矩阵；根据研究区域边界条件和源汇项得到右端项；采用有效的计算方法求解总刚度矩阵和右端项联立方程组；求得研究区域上每个节点的水头。通过多种模拟试验，得到的结果与解析解相吻合。与现有技术相比，该方法的精度与其相近，但计算时间不到其的10%；在求解大范围，长时间或者复杂问题时，效率大幅度提高。

Description

改进的模拟多孔介质中二维水流运动的多尺度有限元方法

技术领域

本发明涉及水力学领域，具体涉及一种模拟多孔介质中二维水流运动的改进的多尺度有限元方法。

背景技术

水资源问题是当前和人类生存关系至为密切的一个重要问题。世界各国有很多城市的大部分用水都取自地下水。此外，在地质工程活动中，地下水的分布也是必须考虑的因素。因此，研究地下水位的计算方法和模拟，对于测量地下水的分布情况与预报具有非常重要的意义。

地下水流的一般方程由椭圆型方程描述稳定流的分布,其二维形式为：

$-\frac{\∂}{\∂x}(K_{\mathrm{xx}}\frac{\∂H}{\∂x}+K_{\mathrm{xy}}\frac{\∂H}{\∂y})-\frac{\∂}{\∂y}(K_{\mathrm{yx}}\frac{\∂H}{\∂x}+K_{\mathrm{yy}}\frac{\∂H}{\∂y})=W,---\left(1\right),$

地下水流的一般方程由抛物型方程描述稳定流的分布,其二维形式为：

$\frac{\∂H}{\∂t}-\frac{\∂}{\∂x}(K_{\mathrm{xx}}\frac{\∂H}{\∂x}+K_{\mathrm{xy}}\frac{\∂H}{\∂y})-\frac{\∂}{\∂y}(K_{\mathrm{yx}}\frac{\∂H}{\∂x}+K_{\mathrm{yy}}\frac{\∂H}{\∂y})=W,---\left(2\right),$

这里H为水头，K_xx、K_xy、K_yx、K_yy分别为xx方向、xy方向、yx方向、yy方向的的渗透系数，W为源汇项。

地下水流方程可以用常规的有限单元方法或者有限差分方法求解。但是这些方法要求单元网格内的介质是均质的，在求解非均质问题时，则必须精细剖分以保证单元内部的渗透系数变化率不大，可以近似为常数。在进行大面积的研究区域中的水流模拟时，精细剖分会产生非常多的节点，对计算机的存储量要求很大并且需要大量的计算时间。因此，水文研究工作者提出了多尺度有限单元法[Hou,T.Y.,and X.H.Wu(1997)]来解决这一问题。

多尺度有限单元法在剖分研究区域时，不要求单元内部的渗透系数必须近似为常数。该方法通过对单元的细剖分，在单元上通过求解简化的椭圆方程来构造基函数，可以很好的通过基函数抓住介质的非均质性质。特别是对于多孔介质，该介质的非均质性一般包含了很多尺度，常常反映在介质的渗透系数的多尺度波动上。如果利用有限单元法求解，为了保证精确度，需要在所有小尺度上求解，这一过程需要非常大的计算量与计算时间。由于多尺度有限单元法的基函数可以抓住介质的尺度的宏观信息，多尺度方法不需要在小尺度上求解，非常适用于工程计算。此外，多尺度有限单元法已经从数学推导和数值模拟上证明了它能够很好的求解椭圆和抛物型问题，并且收敛、精确、高效[Hou,T.Y.,and X.H.Wu(1997),X.H.Wu,and Z.Cai(1999),W.Deng et al.(2008),Ye,S.,Y,Xue,and C.Xie(2004)]。此外，多尺度有限单元法中的超样本技术能够避免网格尺度与物理介质的尺度产生的共振效应，提高收敛速度，减少误差[Hou,T.Y.,and X.H.Wu(1997)]。与有限单元法和有限差分方法相比，多尺度有限单元法对计算机的存储量和计算时间的要求较低，并且能够保证一定的精度。然而，由于传统多尺度有限单元法的三角细剖分法会产生较多的单元内点，导致在求解基函数时需要的计算量较大，如果水流问题的研究区域过于庞大，周期过长，采用传统的多尺度方法需要大量时间与计算量去求解基函数，效率需要提高。

发明内容

本发明针对现有技术中存在的不足，提供了一种模拟多孔介质中二维水流运动的改进的多尺度有限元方法，其可以应用于解决稳定流以及非稳定流的水流问题。该方法模拟得到的结果与解析解非常吻合。在给定细剖分份数、基函数边界条件的情况下，它与传统多尺度有限单元法的精度接近，但需要的计算时间仅有传统方法的10%。

本发明所述改进的模拟多孔介质中二维水流运动的多尺度有限元方法，包括以下步骤：

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

（2）在每一粗网格单元中，以一个或多个内点为中心，采用放射状的三角形单元进行细剖分，得到该粗网格单元的细网格单元；

（3）根据渗透系数K以及基函数的边界条件，求解退化的椭圆型问题，确定基函数，形成有限元空间；

（4）计算各粗网格单元的刚度矩阵，相加得总刚度矩阵；根据研究区域的边界条件、源汇项，计算右端项，形成有限元方程；

（5）提供有限元方程的有效解法。

上述步骤（1）中，所述形成粗网格单元的剖分采用的是三角形单元剖分。

上述步骤（3）中，细网格单元上的渗透系数K、源汇项值近似取这个单元的所有内点上的渗透系数、源汇项的平均值

本发明提出了一种全新的多尺度有限单元法的细剖分方法，通过采用以内点为中心的放射状的细剖分，与传统方法相比，在将单元剖分为同样份数时，本发明需要的内部节点（即未知数）个数更少，从而减少了计算基函数所需的计算量和计算时间。该方法能有效处理非均质多孔介质地下水流问题，且简单易行，高效精确。在求解大范围，长时间或者复杂问题时，该方法的效率要高很多。

以同一三角形Δ_ijk区域为例，图1（a)为用一个内点将其剖分为27份，图1（b)为用三个内点将其剖分为25份，图1(c)为传统多尺度有限单元法的细剖分方法，需要6个内点，将Δ_ijk剖分为25份；从剖分效率上本发明要远大于传统方法。在计算三角单元Δ_ijk的多重尺度基函数时，如采用图1（a)，由于只有一个内点，可以直接以获得一个基函数的表达式；采用图1（b)，由于有3个内点，需要求解一个3×3的矩阵以获得一个基函数的表达式；采用图1（c),由于有6个内点，需要求解一个6×6的矩阵以获得一个基函数的表达式；从计算量上本发明要选小于传统方法。通过对多孔介质下的二维稳定流以及非稳定流的连续介质模型、二维稳定流以及非稳定流的渐变介质模型，二维稳定流以及非稳定流的突变介质模型，二维稳定流潜水介质模型，二维高度非均质介质模型（六种不同尺度）的数值模拟，发现本发明与解析解吻合的很好，精度与传统方法接近，计算时间只有传统方法的10%。结果显示：针对同一个水流问题，采用图1（b)细剖分获得的结果精度高于采用图1（a)获得的结果的精度，但计算时间也略大；采用图1（c)细剖分获得的结果精度略高于与采用图1（b)获得的精度，但所需计算时间远大于采用（b)所需的时间。

附图说明

图1：（a)：改进的细剖分方法，用一个内点将一个粗网格单元剖分为27份;(b)改进的细剖分方法，用三个内点将一个粗网格单元剖分为25份;(c)传统剖分方法，将粗网格剖分为25份。

图2：二维稳定流的连续介质模型，各方法在y=100m剖面上水头的绝对误差。

图3：二维非稳定流的渐变介质模型（冲击平原抽水模型），各方法在y=5200m剖面上模拟的水头。

图4：二维非稳定流突变介质模型，各方法在y=5000m剖面上模拟的的水头。

具体实施方式

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

LFEM:古典有限单元法。

LFEM-F:古典有限单元法(精细剖分）。

MSFEM-L:传统多尺度有限单元法,使用线性边界条件。

MSFEM-O:传统多尺度有限单元法,使用振荡边界条件。

MSFEM-os-O:传统多尺度有限单元法,使用振荡边界条件，使用超样本技术。

MMSFEM-p-L:改进的多尺度有限单元法,使用线性边界条件,使用p个内点进行细剖分。

MMSFEM-p-O:改进的多尺度有限单元法,使用振荡边界条件,使用p个内点进行细剖分。

MMSFEM-p-os--O:改进的多尺度有限单元法,使用振荡边界条件,使用p个内点进行细剖分，使用超样本技术。

实施例1:二维稳定流的连续介质模型

研究区为一正方形单元：Ω=[50m,150m]×[50m,150m]，渗透系数K(x,y)＝x²m/d。水流方程为公式（1），边界条件为定水头边界条件

此模型有解析解：H＝3x²+y²。

采用LFEM,LFEM-F,MSFEM-L,MSFEM-O,MMSFEM-1-L,MMSFEM-1-O,MMSFEM-1-os-O求解此模型。其中，LFEM-F将研究区剖分为1800份，其他方法将研究区剖分为200个粗网格单元。MSFEM采用传统三角剖分方法将每一粗网格细剖分为9个单元，MMSFEM采用改进细剖分方法将每一粗网格细剖分为9个单元。超样本技术采用的超样本单元为粗网格单元的1.01倍。

图2为上述方法计算的水头，在y=100m这个剖面的绝对误差。从图2可知，FEM方法的误差是最大的；MSFEM-L与MMSFEM-L的精度非常接近；MMSFEM-1-O,LFEM-F,MSFEM-O都获得非常精确的解，它们的误差非常接近；MMSFEM-1-os-o获得了最好的结果。MMSFEM-1-os-o的结果要好于LFEM-F，这和T.Y.Hou和S.J.Ye在采用MSFEM-os-O模拟水位模型时得到的结果一致[Hou,T.Y.,and X.H.Wu(1997)，Ye,S.,Y,Xue,and C.Xie(2004)]。

实施例2：二维非稳定流的渐变介质模型（冲击平原抽水模型）

研究区为一正方形单元：Ω=[0,10km]×[0m,10km]，渗透系数从研究区的边界左侧到右侧从1m/d增加到250m/d即K(x,y)＝1+x/40m/d。水流方程为公式（2），左右边界定水头边界，左边界为10m,右边界为0m，上下位隔水边界。含水层厚度为10m，贮水系数S＝0.00001-0.000009x/1000/m。在坐标（5200m,5200m）处有一抽水井，流量为1000m³/d，抽水时间为5天，时间步长为1天。初始时刻的水头H₀(x,y)＝10-x/1000m。此模型没有解析解，因此，采用LFEM-F的解作为标准参照。

采用LFEM,LFEM-F,MSFEM-O,MMSFEM-3-O方法求解此模型。其中，LFEM-F将研究区剖分为125000个三角形单元，其他方法将研究区剖分为1250个单元。此模型中，MMSFEM-3-O采用放射状细剖分，MSFEM-O采用传统细剖分均将粗网格单元细剖分为100个三角形单元。图3为各方法y=5200m剖面的水头。MMSFEM与MSFEM的精度非常接近。此模型中，MMSFEM所需的CPU时间为3秒,MSFEM为308秒，MMSFEM具有更高的的效率。

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

研究区为一正方形单元：Ω=[0,10km]×[0m,10km]，渗透系数在x=2480m这个剖面上发生突变，即当x＜2480m,K＝2m/d;x≥2480m,K＝1000m/d。水流方程为公式（2），左右边界定水头边界，左边界为10m,右边界为0m，上下位隔水边界。在坐标（5000m,5000m）处有一抽水井，流量为6000m³/d，抽水时间为3天，时间步长为1天。含水层厚度为10m，贮水系数x＜2480m,S＝0.000002/m;x≥2480m,S＝0.0005/m。初始时刻的水头H₀(x,y)＝10-x/1000m。此模型没有解析解，因此，采用LFEM-F的解作为标准参照。

采用LFEM,LFEM-F,MSFEM-O,MMSFEM-3-O求解此模型。其中，LFEM-F将研究区剖分为80000个三角形单元，其他方法将研究区剖分为400个单元。此模型中，MMSFEM采用放射状细剖分，MSFEM采用传统细剖分均将粗网格单元细剖分为100个三角形单元。图4为各方法y=5000m剖面的水头。MMSFEM与MSFEM的精度非常接近。此模型中，MMSFEM所需的CPU时间为1秒,MSFEM为73秒。

实施例4：二维稳定流潜水流模型（非线性模型）

水流方程为潜水流方程：

-▽·K(x,y,H)▽H＝W,

$K(x,y,H)=\left(\begin{array}{cc}T(H-b)& 0\\ 0& T(H-b)\end{array}\right)$

由于是非线性模型，此方法需要迭代：

-▽·K(x,y,H^(n-1))▽H⁽ⁿ⁾＝W,

并设定的迭代误差为η，即迭代直至|H⁽ⁿ⁾-H^(n-1)|＜η。

此模型研究区为：Ω=[0,1m]×[0,1m]，边界水头为定水头边界且均为0.b=-4m,T＝(1+x)(1+y).η＝0.00001,H₀＝0.此模型有解析解：H＝xy(1+x)(1+y)。

采用MMSFEM-3-L求解此模型。将研究区剖分为800、1800、3200个单元分别进行三次模拟。每次模拟中，均将粗网格单元细剖分为25个单元。三次模拟需要的迭代次数均为4次，所花费的CPU时间均小于1s。

Claims (3)

1.一种改进的模拟多孔介质中二维水流运动的多尺度有限元方法，其特征在于，包括以下步骤：

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

（2）在每一粗网格单元中，以一个或多个内点为中心，采用放射状的三角形单元进行细剖分，得到该粗网格单元的细网格单元；

（3）根据渗透系数K以及基函数的边界条件，求解退化的椭圆型问题，确定基函数，形成有限元空间；

（4）计算各粗网格单元的刚度矩阵，相加得总刚度矩阵；根据研究区域的边界条件、源汇项，计算右端项，形成有限元方程；

（5）提供有限元方程的有效解法，求得研究区域上每个节点的水头。

2.根据权利要求1所述改进的模拟多孔介质中二维水流运动的多尺度有限元方法，其特征在于步骤（1）中，所述形成粗网格单元的剖分采用的是三角形单元剖分。

3.根据权利要求1或2所述改进的模拟多孔介质中二维水流运动的多尺度有限元方法，其特征在于：步骤（3）中，细网格单元上的渗透系数K、源汇项值近似取这个单元的所有内点上的渗透系数、源汇项的平均值。

Images