Page 31 - 2024年第55卷第12期

了比较和分析。其中，ＭＳＦＥＭ－ Ｊ的粗网格单元数目远小于其他两种方法的网格单元数，而其细网格
与 Ｚｈｏｕ － Ｆ、Ｙｅｈ － Ｆ的网格形状、大小、总数均相同，且三种方法均使用了线性基函数边界条件。此
外，所有程序均使用 Ｃ ＋＋ 编写，并在同一软、硬件条件下运行。三种方法的程序均已进行了全面的测
试，能够与解析解精确吻合，具有较高的可靠性。
３．１ 含水平介质交界面的非均质各向异性二维地下水稳定流
问题 本算例是由 Ｃａｉｎｅｌｌｉ等 ［９］ 工作中的无量纲算例修改而
成的，控制方程为式（ ３）。研究区为［０，１］ × ［０，１］的正方
形（图 ３），含有一水平交界面 ＡＢ（ｙ ＝ ０．５ ），交界面下侧区域
为 Ｚｏｎｅ１，上侧区域为 Ｚｏｎｅ２。Ｚｏｎｅ１和 Ｚｏｎｅ２的渗透系数
具有各向异性：Ｚｏｎｅ１渗透系数为 Ｋ ＝ １０（１ ＋ ｘ）（１ ＋ ｙ），Ｋ ＝
ｙ
ｘ
５０ （１ ＋ ｘ）（１ ＋ ｙ），Ｚｏｎｅ２的渗透系数 Ｋ ＝ １００ （１ ＋ ｘ）（１ ＋ ｙ），Ｋ ＝
ｘ
ｙ
５００ （１ ＋ ｘ）（１ ＋ ｙ）。本算例具有解析解，如式（１６）。研究区四
周的第一类边界条件以及算例的源汇项均基于水头解析解设
定 ［９］ 。达西流速解析解则可以根据达西定律，由水头解析解
图 ３ 研究区粗网格剖分示意图
和渗透系数得到。

]( )
[ （２ｘ ＋ ９ ）ｙ － ２ ｘ ＋ ９ ３ ４
ｓｉｎ ２ ２０ ｙ＋ １ Ｚｏｎｅ１
Ｈ（ｘ，ｙ） ＝ ３ ４ （１６）
２ｙ ＋ ９
２ｘ ＋ ９
[ （２ｘ ＋ ９ ）（２ｙ ＋ ９ ） ]( )( )
ｓｉｎ ４０ － ２０ ２０ ＋ １ Ｚｏｎｅ２
本算例中，ＭＳＦＥＭ － Ｊ将研究区划分为 ３２００（４０ × ４０ × ２）个三角形粗网格单元，每个粗网格单元划分为
４（２ × ２）个三角形细网格单元；Ｚｈｏｕ － Ｆ和 Ｙｅｈ － Ｆ将研究区划分为 １２８００（８０ × ８０ × ２）个三角形粗网格单
元。ＭＳＦＥＭ－ Ｊ采用 ＭＳＦＥＭ计算水头，Ｚｈｏｕ － Ｆ和 Ｙｅｈ － Ｆ均是通过精细剖分的有限元方法求解的水头。
虽然 ＭＳＦＥＭ－ Ｊ的粗网格数远小于其他两者的单元数，但它们所模拟的水头精度十分接近，ＭＳＦＥＭ－ Ｊ、
Ｚｈｏｕ － Ｆ、Ｙｅｈ － Ｆ三种方法与解析解的相对误差分别为 ０．０２％、０．０１％、０．０１％，说明 ＭＳＦＥＭ－ Ｊ模拟的
水头与精细剖分的有限元法基本相同，具有很高的精度。同时，虽然 ＭＳＦＥＭ－ Ｊ仅计算了粗尺度节点
上的水头值，但细尺度水头值可由式（ ５）直接获得，故其所获的水头解数量与另两种方法相同。
图 ４展示了 ＭＳＦＥＭ－ Ｊ、Ｚｈｏｕ － Ｆ和 Ｙｅｈ － Ｆ所模拟的达西流速 Ｖ，图 ４（ａ）—（ｃ）为各方法整个研究
ｘ
区的全局达西流速；图 ４（ｄ）（ｅ）分别为交界面 ＡＢ和截面 ｘ ＝ ０．５ 上的 Ｖ，图 ４（ｆ）则为 Ｖ在交界面 ＡＢ
ｘ
ｘ
上的相对误差。其中，ＭＳＦＥＭ－ Ｊ结合了多尺度基函数和 Ｙｅｈ的伽辽金有限元模型模拟达西流速 ［１０，１２］ ，
并通过 ＪＵＭＰ向量修正 ［５］ ；Ｚｈｏｕ － Ｆ和 Ｙｅｈ － Ｆ均使用了 Ｙｅｈ的伽辽金有限元模型模拟达西流速 ［１０］ ，但
Ｚｈｏｕ － Ｆ采用了 ＪＵＭＰ向量修正。从图 ４（ａ）—（ｃ）可以看出，ＭＳＦＥＭ－ Ｊ可通过多尺度基函数捕获粗网
格的细尺度特征，在非交界面区域能够得到与后两者精度相当的达西流速。ＭＳＦＥＭ－ Ｊ和 Ｚｈｏｕ － Ｆ通
过 ＪＵＭＰ向量修正了介质交 界面 处的达西 流速，没有 像 Ｙｅｈ － Ｆ一 样 在介 质 交 界面 附近 出现 较大的
误差。
图 ４（ｄ）—（ｆ）将三种方法所模拟的达西流速 Ｖ和解析解进行了进一步比较。根据折射定律，介质
ｘ
交界面的切向达西流速需按渗透系数成比例。因此，图 ４（ｄ）中的 ＭＳＦＥＭ－ Ｊ和 Ｚｈｏｕ － Ｆ以及解析解均
在 ＡＢ处具有两个不同数值，且 ＭＳＦＥＭ－ Ｊ和 Ｚｈｏｕ － Ｆ的曲线均与解析解几乎重合。由于 Ｙｅｈ － Ｆ未考虑
折射定律，其在 ｘ ＝ ０．５处仅有一个数值且与解析解相差巨大。同时，如图 ４（ｅ）所示，ＭＳＦＥＭ－ Ｊ和
Ｚｈｏｕ － Ｆ在非交界面区域依然获得了精确的数值；图 ４（ｆ）展示了三种方法在交界面处的相对误差，显
示 Ｙｅｈ － Ｆ的相对误差和 ＭＳＦＥＭ－ Ｊ和 Ｚｈｏｕ － Ｆ具有数量级的差异。由于 Ｙｅｈ － Ｆ是采用平均的方法计算
达西流速的 ［１０］ ，其交界面附近的节点的解对交界面处的解的误差进行了 “补偿”，因而也产生了较大
误差。ＭＳＦＥＭ － Ｊ和 Ｚｈｏｕ － Ｆ在交界面 ＡＢ上下侧的 Ｖ总体平均相对误差分别为 ０．０８％和 ０．０７％，而 Ｙｅｈ － Ｆ
ｘ
在 ＡＢ界面上侧的 Ｖ平均相对误差为 １００％，下侧也在 ４５％左右。
ｘ
图 ５统计了三种方法在本算例模拟水头和两个方向的达西流速所使用的总方程阶数和所消耗的 ＣＰＵ
４
— １ ４ ３ —