模拟多孔介质中二维水流运动的高效多尺度有限元方法
【技术领域】
[0001] 本发明属于水力学技术领域，具体涉及一种模拟多孔介质中二维水流运动的高效 多尺度有限元方法。
【背景技术】
[0002] 地下水资源是水资源的重要组成部分，是工业、农业和城市用水的重要来源之一。 在水文地质学中地下水位能够反映地下水所具有的机械能的大小。地下水的分布情况与工 程的实施方案、施工方式、施工时间、工程资金等因素关系密切;所以，研究地下水位相关的 数值计算方法，对于分析地下水分布、运动情况是非常必要，具有重要的研究价值。
[0003] 传统有限单元法是目前常用的地下水数值计算方法之一，在水文地质领域应用十 分广泛。然而，该方法受到单元内部岩性必须相同的限制，在处理非均质多孔介质中的地下 水问题时，必须通过精细剖分的方式保证精度，在研究区域较大时，需要大量的计算时间和 空间，效率较低。因此，国外数学家提出了多尺度有限单元法(Hou and Wu 1997)用于解决 这一问题，是一种既可以减少剖分单元数又能保证计算结果精度的高效方法。该方法突破 了有限单元法的单元内部岩性必须相同的限制，在求解非均质多孔介质中的地下水问题 时，该方法可以使用较大的单元。同时，该方法还能够通过基函数满足粗网格单元的局部椭 圆问题来抓住细尺度的岩性性质，从而节约计算消耗并保证精度(Hou and Wu 1997，Ye et al.2004，Xie et al2014)。鉴于上述传统有限单元法的局限性和多尺度有限单元法的优越 性，开展多尺度有限单元法的改进算法具有重要的理论和实际意义。近年来，经济、科技发 展迅速，人们对于地下水问题越来越关注，也谋求求解更加复杂的地下水问题，相应研究面 积、水文周期也越来越大，如地面沉降问题、海水入侵问题等。在求解此类问题时，多尺度有 限单元法的基函数构造消耗过高，需要进一步提高。
【发明内容】
[0004] 针对于上述现有技术的不足，本发明的目的在于提供一种模拟多孔介质中二维水 流运动的高效多尺度有限元方法，该方法运用区域分解技术改进了基函数的构造算法和剖 分方式，将基函数的构造问题分解为若干个子问题，分批求解未知项，能够大幅降低构造基 函数所需的计算消耗，以解决现有技术中求解复杂的地下水问题时基函数构造消耗过高等 问题。
[0005] 为达到上述目的，本发明的一种模拟多孔介质中二维水流运动的高效多尺度有限 元方法，包括步骤如下：
[0006] (1)根据所要模拟的研究区域确定边界条件，设定粗网格单元尺度，剖分该研究区 域，得到粗网格单元；
[0007] (2)设定中网格单元尺度，剖分上述粗网格单元，得到中网格单元；
[0008] (3)设定细网格单元尺度，剖分上述中网格单元，得到细网格单元；
[0009] (4)根据渗透系数K以及基函数的边界条件，以中网格单元为最小子单元，在粗网 格单元上求解退化的椭圆型问题，确定所有中网格单元顶点处的基函数值；
[0010] (5)运用区域分解技术将上述粗网格上的局部椭圆问题分解为每个中网格单元上 的子问题；
[0011] (6)根据渗透系数K、中网格单元顶点处的基函数值以及改进的基函数边界条件得 到所有子问题的边界条件，以细网格单元为最小子单元，在每个中网格单元上求解子问题 得到基函数在每个中网格单元中所有节点上的值；
[0012] (7)计算各粗网格单元的刚度矩阵，相加得总刚度矩阵；根据研究区域的边界条 件、源汇项，计算右端项，形成有限元方程；
[0013 ] (8)采用cho 1 e shy分解法，求得研究区域上每个节点的水头。
[0014] 优选地，上述的步骤(1)中，采用三角形单元剖分研究区域，以形成粗网格单元。
[0015] 优选地，上述的步骤(2)中，采用三角形单元剖分粗网格单元，以形成中网格单元。
[0016] 优选地，上述的步骤(3)中，采用放射状的三角形单元剖分方式剖分中网格单元， 以形成细网格单元。
[0017]优选地，上述的步骤(4)中，所述的研究区域最小子单元上的渗透系数K取这个单 元的所有顶点上的渗透系数平均值。
[0018] 优选地，上述的步骤(6)中，所述的研究区域最小子单元上的渗透系数K取这个单 元的所有顶点上的渗透系数平均值。
[0019] 优选地，上述的步骤(7)中，细网格单元上的源汇项值取这个单元的所有顶点上的 源汇项的平均值。
[0020] 优选地，上述的步骤(7)中，粗网格单元上的源汇项值取这个单元中所有细网格单 元的源汇项的平均值。
[0021] 本发明的有益效果：
[0022] 1.令多尺度有限单元法可以将其每个粗网格单元上的基函数构造问题转换为若 干个子问题，从而分批求解未知项，大幅降低基函数的构造消耗；
[0023] 2.可以将整个研究区剖分为粗、中、细三重网格单元，单元灵活度更高，单元的抗 畸形能力较强；
[0024] 3.改进了基函数的边界条件，不要求边界端点必须分别为0和1;
[0025] 4.在研究区粗网格数目相同且粗网格的细网格数目相同时，本发明获得的精度和 传统多尺度有限单元法相近，计算时间大幅降低；
[0026] 5.能够有效求解非均质多孔介质地下水稳定流和非稳定流问题，如二维稳定流的 连续介质问题、二维非稳定流的渐变介质问题，二维稳定流的突变介质问题，二维稳定流的 多尺度介质模型(两种不同尺度），二维潜水流模型的数值模拟，且原理简单，高效精确；
[0027] 6.在求解大尺度，长时间等高计算量问题时，本发明的效率更高。
【附图说明】
[0028] 图1为高效多尺度有限单元法的粗网格单元剖分方式示意图。
[0029] 图2为高效多尺度有限单元法的中网格单元剖分方式示意图。
[0030]图3为二维稳定流的连续介质模型，子例1.1中各数值方法在y=100m截面处的水 头绝对误差值的示意图。
[0031] 图4为二维非稳定流的渐变介质模型，子例2.1中各数值方法在y = 5200m截面处的 水头曲线的不意图。
[0032] 图5为二维非稳定流的渐变介质模型，子例2.2中各数值方法所计算的在y = 5200m 截面处的水头曲线的示意图。
[0033] 图6为二维潜水流模型，各数值方法在y = 0.5截面处水头的相对误差值的示意图。
【具体实施方式】
[0034] 为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说 明，实施方式提及的内容并非对本发明的限定。
[0035] 参照图1所示，本发明的一种模拟多孔介质中二维水流运动的高效多尺度有限元 方法，其基函数构造过程如下：
[0036]在粗网格单元Δ,诉上考虑构造方程，即退化的椭圆方程：
[0037] V · Α'νψ, = 0,
[0038] 其中Κ为渗透系数，Ψi为粗网格单元Δ ijk(图1)在i点的基函数。
[0039] 步骤一:在粗网格单元△#(图1)上，以中网格单元为最小单元，求解构造方程，可 以得到基函数在0点（即Mo点）的值;通过基函数的边界条件可以得到基函数在其余中网格 顶点（即i，j，k，e，f，p，q)上的值。
[0040]步骤二:运用区域分解技术将Δ uk上的构造问题转换为中网格单元Δ Qgl、Δ _、 Δ。#、AQfg、Δ<^、AqkgApfe、Apej上的子问题；中网格单元的某一边界ξτι上的改进的基函数 线性、振荡边界条件如下：
[0043] 步骤三:在每个中网格单元上（图2)，以细网格单元为最小单元，求解子问题即可 得到基函数在内点Μι_Μ8的值。
[0044] 其中：
[0045] 在步骤一中，在通过基函数的边界条件获得基函数在中网格顶点上的值时，可以 对边界加密提升所获精度。
[0046] 在步骤二中，中网格单元数目并不仅限于图1所示的8个，可以根据需要增多。
[0047] 完成上述基函数的构造后，通过变分形式即可求得单元刚度矩阵，相加得总刚度 矩阵，运用cholesky分解法求解总刚度矩阵和右端项的方程组即可获得水头值。
[0048] 下面结合具体实施例对本发明做进一步的解释，其中涉及一些简写符号，以下为 注解:LFEM:有限单元法；
[0049] LFEM-F:有限单元法(精细剖分）；
[0050] MSFEM-L:传统多尺度有限单元法，使用线性边界条件；
[0051 ] MSFEM-0:传统多尺度有限单元法，使用振荡边界条件；
[0052] EMSFEM-L:高效多尺度有限单元法，使用线性边界条件；
[0053] EMSFEM-0:高效多尺度有限单元法，使用振荡边界条件。
[0054]本发明所有实施例在基函数构造的步骤一中，利用粗网格的基函数边界条件得到 中网格顶点基函数值时，将粗网格每边分为8份，并使用振荡边界条件。
[0055] 实施例1:二维稳定流的连续介质模型
[0056] 研究区为正方形区域：Ω =[50m，150m] X [50m，150m]，渗透系数K(x，y)=x2m/d， 研究方程为稳定流方程：
[0058] 边界条件为定水头边界条件//|iu = V -3/ <源汇项为〇，此模型有解析解:!1 = ￥- 3y2〇
[0059] 子例 1 · 1:采用 LFEM，LFEM-F，MSFEM-L，MSFEM-0，EMSFEM-L 和 EMSFEM-0 求解。LFEM-F 将研究区剖分为88200个单元，其他