Page 102 - 2021年第52卷第8期

-∇∙ ( K∇Ψ i ) = 0 (x，y )∈ Δ ijk （4）

在设定多尺度基函数的边界条件后 ［14］ ，式（4）适定。
在粗网格单元 Δ 中的任意细网格单元 Δ 内，Ψ 可被线性表示为 ［14］ ：
ijk abc i
)
Ψ (x，y = Ψ ( )N + Ψ ( ) b N + Ψ ( ) c N c (x，y )∈ Δ abc （5）
a
i
i
i
b
i
a
a
式中： N 、 N 、 N 为有限元线性基函数； Ψ ( )、Ψ ( ) b 、Ψ ( ) c 为 Δ abc 顶点处的 Ψ 的值。对式
i
b
a
i
i
c
i
（4）进行伽辽金变换，再将式（5）代入，可以得到关于 Ψ 的方程组，求解可得 Ψ 在 Δ 中所有节点的
i
i
ijk
值。其他多尺度基函数的构造过程与 Ψ 类似，这里不再赘述。
i
2.5 运用 MSFEM-D 模拟二维地下水问题
2.5.1 运用 MSFEM-D 模拟地下水水头 以地下水二维非稳定流问题为例，考虑如下方程：
S ∂H - ∂ æ ç K ∂H ö ÷ - ∂ æ ç K ∂H ö ÷ = W (x，y )∈ Ω （6）
∂t ∂x è x ∂x ø ∂y è y ∂y ø
式中：W 为源汇项；Ω为研究区；S 为贮水系数。在地下水模拟问题的定解条件确定后，式（6）适定。
在水头模拟时，无需运用交界面进行分区，可以直接剖分研究区。在式（6）两边乘以多尺度基函
数Ψ ，进行伽辽金变换，再离散到每一粗网格单元 Δ 上，可得：
i
ijk
é æ ∂H ö ∂Ψ i æ ∂H ö ∂Ψ i ∂H ù
 ê ê ç K x ∂x ÷ + ç K y ∂y ÷ ∂y + S ∂t Ψ dxdy =  WΨ dxdy （7）
ú ú i
i
Δ ë è ø ∂x è ø û Δ
ijk ijk
根据文献［14］，在粗网格单元 Δ 上，Ψ 可被线性表示为：
ijk
i
)
)
)
H (x，y = H Ψ (x，y + H Ψ (x，y + H Ψ (x，y ) （8）
i
j
j
i
k
k
式中 H 、H 、H 分别为水头 H 在粗网格单元 Δ 顶点 i、j、k 处的值。
i
j
k
ijk
将式（7）离散到细网格上，并代入式（5）、式（8）可得式（7）的具体表达式。再运用 Ψ 、Ψ 乘以式
j k
（6）两端，并重复上述过程，可得粗网格 Δ 上关于水头 H 的单元刚度矩阵。联立所有粗网格上的单
ijk
元刚度矩阵，结合源汇项 W、初始条件、边界条件即可得到右端项，从而得到研究区上的总方程
组，结合 Cholesky 分解法等矩阵求解方法即可获得研究区上的水头在研究区粗尺度节点上的值。
2.5.2 运用 MSFEM-D 模拟达西流速 在获得研究区的水头分布后，MSFEM-D 可以运用这些已知的
［14］ ［8］ ［3］
水头值，结合区域分解技术、MSFEM 、Yeh 的伽辽金有限元法 以及 Zhou 的 Jump 函数 进行达西
流速的求解。其中，水头值可以提供研究区的全局信息，区域分解技术可以分离介质避免互相影
响、MSFEM 基函数可以提高模拟效率、Yeh 的伽辽金法可以保证达西流速的连续性。以 V 的求解过
x
程为例，在图 2 所示的研究区Ω上考虑达西方程：
V = -K x ∂H （9）
∂x
x
如 2.2.2 节所述，MSFEM-D 已经将研究区Ω划分为 Zone 1、Zone 2 两个子区域。MSFEM-D 运用区
域分解技术将研究区上的达西流速求解问题划分为 Zone 1、Zone 2 这两个子区域的子问题。设先求解
Zone 1 上的子问题，则需在 Zone 1 上考虑式（9）。设子区域 Zone 1 上的粗尺度网格数目为γ ，粗尺度
1
节点数目为 Nu ，运用Ψ ，I = 1，2，…，Nu 乘以达西方程的两端，并在 Zone 1 上积分，并将 Zone 1
1 I 1
离散为γ 个粗网格 Δ ，可得：
1 ijk
γ γ
1 1 ∂H
V Ψ dxdy = - å K Ψ dxdy (x，y （10）
å x I x ∂x I )∈ Zone 1
1 1
在任意粗网格 Δ 上，根据现有文献中的理论 ［8，14］ ，V 可以被多尺度基函数线性表示为：
ijk
x
)
)
)
)
V (x，y = V ( x ，y Ψ (x，y + V (x ，y )Ψ (x，y + V (x ，y )Ψ (x，y ) （11）
x x i i i x j j j x k k k
)
)
式中 V ( x ，y 、V ( x ，y 、V ( x ，y k ) 为达西流速 V 在粗网格 Δ 顶点 i、j、k 上的值。
x
i
i
j
x
j
x
ijk
x
k
— 982 —