模拟二维水流运动的三次样条多尺度有限元方法
【技术领域】
[0001] 本发明属于水力学技术领域，具体涉及一种模拟多孔介质中二维水流运动的三次 样条多尺度有限元方法。
【背景技术】
[0002] 地下水是十分重要的水体，是水资源的重要组成部分。世界各国有很多城市的大 部分用水都取自地下水。同时，在地质工程活动中，地下水的分布和运移是一项重要因素。 因此，研究地下水位和渗流速度的计算方法和数值模拟，对于考察地下水分布规律及预测 具有十分重要的意义。
[0003] 多尺度有限单元法（Hou and Wu 1997)是由科学工作者提出的一种求解非均质介 质中地下水流问题的方法。多尺度有限单元法的核心思想是在网格单元上令基函数满足局 部微分算子，从而通过基函数抓住介质的非均质性质。和常规的有限单元方法或者有限差 分方法不同，多尺度有限单元法不要求单元内部的渗透系数必须近似为常数。在处理非均 质问题时，多尺度有限单元法可以使用较少的网格单元数，无需对研究区精细剖分。因此， 多尺度有限单元法比有限单元方法或有限差分方法的计算效率更高，在进行大尺度研究区 域中的水流模拟时能够节约大量计算成本。同时，许多科学工作者已经从数学推导和数值 模拟上证明了多尺度有限单元法能够很好的求解椭圆和抛物型问题，并且高效、收敛、精确 (Hou and Wu 1997, Hou et al. 1999, Deng et al. 2008, Ye et al. 2004)。然而，和有限单 元法一样，多尺度有限单元法无法获得连续的节点水头导数，导致其无法通过达西定律获 得连续、精确的达西渗流速度场，不能对地下水运动状态进行精确描述。
【发明内容】
[0004] 针对于上述问题，本发明的目的在于提供一种模拟二维水流运动的三次样条多尺 度有限元方法，该方法将三次样条法（Zhang et al. 1994)与多尺度有限单元法有机结合， 以解决现有技术中无法通过达西定律获得连续、精确的达西渗流速度场，不能对地下水运 动状态进行精确描述的问题。
[0005] 为达到上述目的，本发明的模拟二维水流运动的三次样条多尺度有限元方法，包 括步骤如下：
[0006] (1)根据所要模拟的研究区域确定边界条件，设定粗网格单元尺度h，剖分该研究 区域，得到粗网格单元，对每一个粗网格单元进行细剖分得到细网格单元；
[0007] (2)根据渗透系数K以及基函数的边界条件，求解退化的椭圆型问题，确定基函 数，形成有限元空间；
[0008] (3)计算各粗网格单元的刚度矩阵，相加得总刚度矩阵；根据研究区域的边界条 件、源汇项，计算右端项，形成有限元方程；
[0009] (4)提供有限元方程的有效解法，求得研究区域上每个节点的水头；
[0010] (5)在研究区每条网格单元线上运用三次样条技术求解水头的一阶导数值；
[0011] (6)通过达西定律求解研究区节点上的达西渗流速度Γ;
[0012] (7)通过研究区域每个节点的水头及基函数得到每一个粗网格单元上每一个节点 水头；
[0013] (8)在粗网格每条单元线上运用三次样条技术求解水头的一阶导数值；
[0014] (9)通过达西定律求解每一粗网格单元上的达西渗流速度Vf。
[0015] 进一步地，所述的步骤（1)中，采用矩形单元剖分研究区域，以形成粗网格单元。
[0016] 进一步地，所述的步骤（1)中，采用三角形单元剖分粗网格单元，以形成细网格单 J L· 〇
[0017] 进一步地，所述的步骤（2)、（3)中，细网格单元上的渗透系数K、源汇项值近似取 这个单元的所有顶点上的渗透系数、源汇项的平均值。
[0018] 进一步地，所述的步骤（5)、（6)与步骤（7)、（8)、（9)相互独立。
[0019] 本发明的有益效果：
[0020] 1.本发明令多尺度有限单元法可以获得其研究区及粗网格单元节点上的达西渗 流速度，并保证它们的连续性；
[0021] 2.与现有计算达西渗流速度的常用方法相比，在研究区网格节点数目相同时，本 发明获得的精度更高；
[0022] 3.本发明可以将计算整个研究区达西渗流速度的问题转化为求解各个粗网格单 元节点的达西渗流速度的问题，从而节约大量计算成本；
[0023] 4.本发明能有效求解非均质多孔介质地下水流场和达西渗流速度场，且原理简 单，尚效精确；
[0024] 5.在求解大范围，长时间等高计算量问题时，本发明的效率更高。
【附图说明】
[0025] 图1绘示本发明中三次样条法示意图。
[0026] 图2绘示二维连续介质模型，在研究区的每边被剖分为10、20、30等份时各数值法 获得的G的平均相对误差。
[0027] 图3绘示二维高振荡水头模型，各数值法获得的水头在截面y = 0. 7处的相对误 差。
[0028] 图4绘示二维高振荡水头模型，各数值法获得的「在截面y = 0. 7处的相对误差。
【具体实施方式】
[0029] 为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说 明，实施方式提及的内容并非对本发明的限定。
[0030] 本发明的模拟二维水流运动的三次样条多尺度有限元方法求解粗网格单元上的 达西渗流速度Vf的过程如下：
[0031] 先通过多尺度有限单元法得到粗网格单元顶点水头，再通过基函数的插值公式获 得所有粗网格单元上的节点（即细尺度节点）上的水头。假设[a，b]是粗网格单元口 1]H上的一条网格线（参照图1所示），有n+1个节点，使用样条函数S近似估计水头H，样条函 数S为三次多项式，S在[a，b]上任一节点的值等于该点的水头值；S的二阶导数在[a，b] 上连续，S在[a, b]上某区间[χτ χτ]的表达式如下：
[0033] 由于S的二阶导数在[a，b]上连续，可以得到一个关于S-阶导数的方程组，求解 可以得到S的一阶导数，该导数在节点上是连续的。利用该导数替代水头的一阶导数，再通 过达西定律，可以得到粗网格单元节点上的连续的达西渗流速度vf。若需要求解研究区域 网格节点（即粗尺度节点）上的达西渗流速度V%只需利用多尺度有限单元法求得的所有 研究区网格节点上的水头值，再通过三次样条函数得到连续的水头导数，通过达西定律求 解即可。具体步骤和上述求解，的步骤类似。通过对多孔介质下的二维连续介质模型，二 维高振荡水头模型，二维振荡介质模型，二维多尺度介质模型（两种不同尺度），二维渐变 介质模型，二维稳定流潜水介质模型的数值模拟，发现本发明与解析解吻合的很好；和常用 求解达西渗流速度的方法相比，在研究区网格节点数目相同时，本发明获得的精度更高；在 总节点数目（包括细尺度节点）相同，本发明精度和现有常用方法相近，并可以将计算整个 研究区达西渗流速度的问题转化为求解各个粗网格单元节点的达西渗流速度的问题，节约 大量计算成本。
[0034] 下面结合具体实施例对本发明做进一步的解释，其中涉及一些简写符号，以下为 注解：
[0035] M:研究区边界被等分份数
[0036] Vx:X方向上的达西渗流速度
[0037] 巧:研究区网格节点上的达西渗流速度Vx，即粗尺度达西渗流速度
[0038] F/ :粗单元网格节点上的达西渗流速度Vx，即细尺度达西渗流速度
[0039] LFEM:传统有限单元法，采用三角形网格单元
[0040] LFEM-F:精细剖分的传统有限单元法，采用三角形网格单元
[0041] CMSFEM :三次样条多尺度有限单元法，采用矩形网格单元
[0042] CST-AS :CMSFEM应用解析的水头求解细尺度达西渗流速度
[0043] CMSFEM-tr :三次样条多尺度有限单元法，采用三角形网格单元
[0044] Method_Yeh:Yeh的线性伽辽金有限元模型，采用LFEM求解水头
[0045] Method-Yeh-F: Yeh的线性伽辽金有限元模型（精细剖分），采用LFEM-F求解水头
[0046] Method-Zhang:三次样条法，采用LFEM求解水头
[0047] CMSFEM和CMSFEM-tr的基函数均使用振荡边界条件，由于三次样条技术要求应用 区域为矩形，CMSFEM-tr仅能够计算粗尺度上的达西渗流速度。在求解时，CMSFEM是在 研究区网格上求解的，而计算F/时，CMSFEM是在每个粗网格单元上求解的。在计算人的相 对误差时，若在某节点Vx，且数值解Vx< 10 则将该节点的相对误差记为0%，否则记为 100%〇
[0048] 实施例1 :二维连续介质模型
[0049] 研究区为一正方形单元：Ω = [0, lm] X [0, lm]，渗透系数 K(x，y) = (1+x) (1+y) m/d，水流方程为：
[0051] 边界条件为定水头边界条件
^此模型有解析解：H = Xy(l-X) (l_y)〇
[0052] 求解粗尺度达西渗流速度：
[0053] 采用CMSFEM、Method-Yeh和Method-Zhang求解此模型。为了保证研究区上网 格节点数目相同，CMSFEM、Method-Yeh和Method-Zhang均将研究区的每边剖分为相同份 数。三种方法均将研究区每边剖分为1〇、20、30等份（M= 10、20、30)。Method-Yeh和 Method-Zhang采用三角形剖分，则将研究区剖分为200、800、1800个三角形单元；CMSFEM采 用正方形剖分，将研究区剖分为1〇〇、400、900个矩形粗网格单元，并将每个粗网格单元剖 分为50个三角形单元。
[0054] Method-Yeh和Method-Zhang通过有限单元法（LFEM)求解水头，当研究区的每边 被分为10、20、30等份时，水头的平均相对误差分别为2. 29 %、0. 58 %和0. 26%，CPU时间 分别为1秒、1秒、2秒。CMSFEM采用多尺度有限单元法（MSFEM)求解水头，当研究区的每 边被分为10、20、30等份时，水头的