DOI: 10.12677/aam.2025.142079, PDF, HTML, XML,  被引量   
作者: 谢富伟, 马晓华^*：长沙理工大学数学与统计学院，湖南 长沙
关键词: Stokes-Darcy耦合问题；混合有限元；统一近似；二阶收敛性；Coupled Stokes-Darcy Problem； Mixed Finite Element； Unified Approximation； Second-Order Convergence

文章引用：谢富伟, 马晓华. 一种具有二阶收敛性的Stokes-Darcy耦合问题的统一混合有限元近似[J]. 应用数学进展, 2025, 14(2): 388-398. https://doi.org/10.12677/aam.2025.142079

1. 引言

由于在环境科学、生理学、生物工程、水文学和一些工业应用方面的重要性，近年来，许多用于解决Stokes-Darcy耦合问题的有效数值方法得到了快速发展。这两个方程在耦合界面处的取值往往受到质量守恒、法向力平衡和BJS (Beavers-Joseph-Saffman)条件的约束[1] [2]。

经过经典变分后，Stokes方程的混合变分格式在任意协调有限元空间中都是正定的。而Darcy方程混合变分格式的正定性与压力离散空间密切相关。因此，适用于解决Stokes问题的有限元空间一般不适用于Darcy问题。有许多研究在两个区域采用不同的有限元空间(例如：[3]-[8])。其中一些文章[7] [8]采用了压力稳健格式，即速度的收敛性不再依赖于压力，使得压力本身的光滑性和低阶压力有限元空间不会影响速度的收敛。

为了在两个区域上使用相同的有限元空间，一些文章引入了惩罚项[9] [10]。另一方面，在[11]中，作者修改了耦合问题的公式，以便他们可以在两个区域上均匀地使用经典的Mini元素或Taylor-Hood元素。为了进一步发展统一的离散化，在[12]中，根据[13]的思想，作者修改了混合公式，使新问题与原问题有相同的解。在此基础上，作者将经典的Mini元素应用于改进的耦合二维Stokes-Darcy问题，使问题变得简单直接地实现。

在[12]中，作者使用的Mini元中压力和速度均基于 $P_1$ 元，但由于在混合变分格式中速度要求梯度以及散度，所以均基于 $P_1$ 元的速度空间会降低压力空间的收敛性。另外，由于 $P_1$ 元在耦合界面上的自由度不足，作者需要在耦合界面上添加额外的稳定函数来保证近似解的稳定性。

本文的其余部分结构如下。在第2节中，我们陈述了经典的Stokes-Darcy耦合问题并对其进行了修改。在第3节中，我们提出了有限元离散化方法。在第4节中，我们构造了一个Fortin算子，在Γ具有特殊形状的条件下，证明了修改的Stokes-Darcy耦合问题的离散稳定条件，并在第5节中将其扩展到一般情况。最后，在第6节中，我们提供了两个数值例子来证实我们的结论。

2. Stokes-Darcy耦合问题的陈述与修改

2.1. Stokes-Darcy耦合问题的陈述

我们考虑一个多边形的有界开区域 $\Omega \subset R^2$ ，其被分为了两个多边形的子区域 ${\Omega }_S$ 和 ${\Omega }_D$ 。我们设 $\bar{\Omega }={\bar{\Omega }}_S\cup {\bar{\Omega }}_D$ ， ${\Omega }_S\cap {\Omega }_D=\varnothing$ ， ${\bar{\Omega }}_S\cap {\bar{\Omega }}_D=\Gamma$ ，其中 $\Gamma$ 表示耦合界面。最后，我们定义 ${\Gamma }_S=\partial {\Omega }_S\backslash \Gamma$ ， ${\Gamma }_D=\partial {\Omega }_D\backslash \Gamma$ (如图1)。

Figure 1. Two-dimensional domain diagram of coupling problem

图1. 耦合问题的二维区域示意图

我们使用 $n_S={\left(n_1^S,n_2^S\right)}^{\text{T}}$ 和 $n_D={\left(n_1^D,n_2^D\right)}^{\text{T}}$ 来代表 $\partial {\Omega }_S$ 和 $\partial {\Omega }_D$ 上的单位外法向量，考虑到两个区域上的函数性质不同，对于定义在 $\Omega$ 上的任意函数，我们定义 $v_S={v|}_{{\Omega }_S}$ 和 $v_D={v|}_{{\Omega }_D}$ 。

在 ${\Omega }_S$ 中，自由流体运动中速度 $u_S$ 和压力 $p_S$ 由Stokes方程控制

$\{\begin{array}{l}-\mu \Delta u_S+\nabla p_S=f_S\text{, in }{\Omega }_S,\\ \mathrm{div}{u}_S=g_S\text{, in }{\Omega }_S,\\ u_S=0\text{ , in }{\Gamma }_S,\end{array}$ (2.1)

其中 $f_S\in {\left(L^2\left({\Omega }_S\right)\right)}^2$ 代表单位质量的力， $g_S\in L^2\left({\Omega }_S\right)$ ， $\mu >0$ 代表粘度。

在 ${\Omega }_D$ 中，多孔介质流动运动Darcy方程控制，速度为 $u_D$ ，压力为 $p_D$

$\{\begin{array}{l}\frac{\mu }{K}{u}_D+\nabla p_D=f_S\text{, in }{\Omega }_S,\\ div{u}_D=g_D\text{, in }{\Omega }_S,\\ u_D\cdot n_D=0\text{ , }\text{ in }{\Gamma }_S,\end{array}$ (2.2)

其中 $f_D\in {\left(L^2\left({\Omega }_D\right)\right)}^2$ 代表单位质量的力， $g_D\in L^2\left({\Omega }_D\right)$ ， $\mu >0$ 代表粘度，K代表渗透率常数。

在 $\Gamma$ 上，我们考虑以下耦合界面条件[4]

$\{\begin{array}{l}{u}_D\cdot n_D+u_S\cdot n_S=0,\\ p_S{n}_S-\mu \nabla u_S{n}_S-p_D{n}_S-\mu \frac{\alpha }{\sqrt{K}}\left(u_S\cdot t\right)t=0\end{array}$ (2.3)

其中第一个方程表示质量守恒，第二个方程是由于法向力的平衡和BJS条件( $\nabla u=\left(\frac{\partial u_i}{\partial x_j}\right),i,j=1,2$ )构成， $\alpha$ 是通过实验证据确定的参数， $t$ 是 $\Gamma$ 上的切向量。

我们用粗体表示向量值函数及其组成的空间。对于任意的子空间 $E\subset \Omega$ ， $H^m\left(E\right)$ 的范数和半范数记作 ${\Vert \cdot \Vert }_{m,E}$ 和 ${\left|\cdot \right|}_{m,E}$ 并且用 ${\left(\cdot ,\cdot \right)}_E$ 表示 $L^2\left(E\right)$ 或 $L^2\left(E\right)$ 的内积。当 $E=\Omega$ 时，则省去定义域下标。另外，

$H\left(div,\Omega \right)=\left\{v\in L^2\left(\Omega \right):div\begin{array}{c}v\end{array}\in L^2\left(\Omega \right)\right\}$ ， $L_0^2\left(\Omega \right)=\left\{q\in L^2\left(\Omega \right):{\displaystyle {\int }_{\Omega }q\text{d}x}=0\right\}$ 。

结合边界及耦合界面条件，速度空间可以定义为

$V=\left\{v\in H\left(div,\Omega \right):v_S\in H^1\left({\Omega }_S\right),v=0\text{on}{\Gamma }_S,v\cdot n_D=0\text{on}{\Gamma }_D\text{and}{v}_D\cdot n_D+v_S\cdot n_S=0\text{on}\Gamma \right\}$ (2.4)

其范数为 ${\Vert v\Vert }_V={\left({\left|v\right|}_{1,{\Omega }_S}^2+{\Vert v\Vert }_{0,{\Omega }_D}^2+{\Vert \begin{array}{c}div\end{array}v\Vert }_{0,{\Omega }_D}^2\right)}^{\frac{1}{2}}={\left({\left|v\right|}_{1,{\Omega }_S}^2+{\Vert v\Vert }_{H\left(div,{\Omega }_D\right)}^2\right)}^{\frac{1}{2}}$ 。接下来，我们定义压力空间

$Q=L_0^2\left(\Omega \right)=\left\{q\in L^2\left(\Omega \right):{\displaystyle {\int }_{\Omega }q\text{d}x}=0\right\}$ (2.5)

其范数为 ${\Vert q\Vert }_Q\begin{array}{c}=\end{array}{\Vert q\Vert }_0$ 。

然后，耦合问题(2.1)~(2.3)的混合变分公式可以表述为：找到 $\left(u,p\right)\subset V\times Q$ 满足

$\{\begin{array}{l}a\left(u,v\right)+b\left(v,p\right)=F\left(v\right),\forall v\in V\\ b\left(u,q\right)=G\left(q\right)\forall q\in Q\begin{array}{c}.\end{array}\end{array}$ (2.6)

其中

$\begin{array}{l}a\left(u,v\right)=\mu {\left(\nabla u,\nabla v\right)}_{{\Omega }_S}+\mu \frac{\alpha }{\sqrt{K}}{\displaystyle {\int }_{\Gamma }\left(u_S\cdot t\right)}\left(v_S\cdot t\right)\text{d}\sigma +\frac{\mu }{K}{\left(u,v\right)}_{{\Omega }_D}\begin{array}{c},\end{array}\\ b\left(v,q\right)=\left(q,div\begin{array}{c}v\end{array}\right)\begin{array}{c},\end{array}\\ F\left(v\right)={\left(f_S,v\right)}_{{\Omega }_S}+{\left(f_D,v\right)}_{{\Omega }_D}\begin{array}{c},\end{array}G\left(q\right)=-{\left(g_S,q\right)}_{{\Omega }_S}-{\left(g_D,q\right)}_{{\Omega }_D}\begin{array}{c}.\end{array}\end{array}$ (2.7)

2.2. Stokes-Darcy耦合问题的修改

注意到，双线性型 $a\left(\cdot ,\cdot \right)$ 仅在 $V_{0,D}=\left\{v\in V:{\left(q,div\begin{array}{c}v\end{array}\right)}_{{\Omega }_D}=0\begin{array}{c},\begin{array}{c}\forall \end{array}\end{array}q\in Q\right\}$ 上正定，而在 $V$ 上不为正定。这会给我们构建统一的有限元空间带来困难，因此我们将修改耦合问题。

根据散度定理，可得 $\begin{array}{c}div\end{array}V\subset Q$ ，同时对于 $\forall q\subset Q$ ，都存在 $v\in V$ 使得 $q=div\begin{array}{c}v\end{array}$ 。因此，可以看到 $\begin{array}{c}div\end{array}V$ 与Q是等距同构的。因此(2.6)中第二个方程Q中的元素可以用 $\begin{array}{c}div\end{array}V$ 中的元素替换，即

$\left(div\begin{array}{c}u\end{array},div\begin{array}{c}v\end{array}\right)=\left(g_D,div\begin{array}{c}v\end{array}\right)$ (2.8)

我们取其在 ${\Omega }_D$ 上的部分便可得

${\left(div\begin{array}{c}u\end{array},div\begin{array}{c}v\end{array}\right)}_{{\Omega }_D}={\left(g_D,div\begin{array}{c}v\end{array}\right)}_{{\Omega }_D}$ (2.9)

将(2.9)添加到(2.6)中的第一个方程，我们可以得到修改后的耦合问题：找到 $\left(u,p\right)\subset V\times Q$ 满足

$\{\begin{array}{l}\tilde{a}\left(u,v\right)+b\left(v,p\right)=\tilde{F}\left(v\right),\forall v\in V\\ b\left(u,q\right)=G\left(q\right)\forall q\in Q\begin{array}{c}.\end{array}\end{array}$(2.10)

其中

$\begin{array}{l}\tilde{a}\left(u,v\right)=\mu {\left(\nabla u,\nabla v\right)}_{{\Omega }_S}+\mu \frac{\alpha }{\sqrt{K}}{\displaystyle {\int }_{\Gamma }\left(u_S\cdot t\right)}\left(v_S\cdot t\right)\text{d}\sigma +\frac{\mu }{K}{\left(u,v\right)}_{{\Omega }_D}+{\left(div\begin{array}{c}u\end{array},div\begin{array}{c}v\end{array}\right)}_{{\Omega }_D}\begin{array}{c},\end{array}\\ \tilde{F}\left(v\right)={\left(f_S,v\right)}_{{\Omega }_S}+{\left(f_D,v\right)}_{{\Omega }_D}+{\left(g_D,div\begin{array}{c}v\end{array}\right)}_{{\Omega }_D}\begin{array}{c}.\end{array}\end{array}$(2.11)

定理2.1 对称连续双线性型 $\tilde{a}\left(u,v\right)$ 在 $V$ 上是正定的，即存在正常数 $C_1$ 和 $C_2$ 使得

$\begin{array}{l}\left|\tilde{a}\left(u,v\right)\right|\le C_1{\Vert v\Vert }_V\Vert u\Vert {|}_V,\begin{array}{cc}& \end{array}\forall v,u\in V\\ \tilde{a}\left(u,u\right)\ge C_2{\Vert u\Vert }_V^2,\begin{array}{cc}& \end{array}\forall u\in V\end{array}$ (2.12)

证明：根据迹定理和holder不等式，我们可得

$\begin{array}{c}\left|\tilde{a}\left(u,v\right)\right|\le \left|\mu {\left(\nabla u,\nabla v\right)}_{{\Omega }_S}\right|+\left|\frac{\mu }{K}{(}_u\right|+\left|{\left(\begin{array}{c}div\end{array}u,\begin{array}{c}div\end{array}v\right)}_{{\Omega }_D}\right|+\left|\mu \frac{\alpha }{\sqrt{K}}{\displaystyle {\int }_{\Gamma }\left(u_S\cdot t\right)}\left(v_S\cdot t\right)\text{d}\sigma \right|\\ \le \mu {\left|u\right|}_{1,{\Omega }_S}{\left|v\right|}_{1,{\Omega }_S}+\frac{\mu }{K}{\Vert u\Vert }_{0,{\Omega }_D}{\Vert v\Vert }_{0,{\Omega }_D}+{\Vert \begin{array}{c}div\end{array}u\Vert }_{0,{\Omega }_D}{\Vert \begin{array}{c}div\end{array}v\Vert }_{0,{\Omega }_D}+\mu \frac{\alpha }{\sqrt{K}}{\left({\displaystyle {\int }_{\Gamma }{\left(u_S\cdot t\right)}^2}\text{d}\sigma \right)}^{\frac{1}{2}}{\left({\displaystyle {\int }_{\Gamma }{\left(v_S\cdot t\right)}^2}\text{d}\sigma \right)}^{\frac{1}{2}}\\ \le \mu {\left|u\right|}_{1,{\Omega }_S}{\left|v\right|}_{1,{\Omega }_S}+\frac{\mu }{K}{\Vert u\Vert }_{0,{\Omega }_D}{\Vert v\Vert }_{0,{\Omega }_D}+{\Vert \begin{array}{c}div\end{array}u\Vert }_{0,{\Omega }_D}{\Vert \begin{array}{c}div\end{array}v\Vert }_{0,{\Omega }_D}+\mu \frac{\alpha }{\sqrt{K}}{\left|u\right|}_{1,{\Omega }_S}{\left|v\right|}_{1,{\Omega }_S}\\ \le C_1{\Vert v\Vert }_V{\Vert u\Vert }_V.\end{array}$

同理可得

$\tilde{a}\left(u,u\right)=\mu {\left(\nabla u,\nabla u\right)}_{{\Omega }_S}+\frac{\mu }{K}{\left(u,u\right)}_{{\Omega }_D}+{\left(\begin{array}{c}div\end{array}u,\begin{array}{c}div\end{array}u\right)}_{{\Omega }_D}+\mu \frac{\alpha }{\sqrt{K}}{\displaystyle {\int }_{\Gamma }\left(u_S\cdot t\right)}\left(u_S\cdot t\right)\text{d}\sigma \ge C_2{\Vert u\Vert }_V^2.$

3. 修改的Stokes-Darcy耦合问题的有限元近似

令 $\{{S_h\}}_{h\text{>}0}$ 为区域 $\Omega$ 的拟一致三角剖分族，每个单元 $T\in S_h$ 位于 ${\Omega }_S$ 或 ${\Omega }_D$ 中(即 $\Gamma$ 仅穿过单元的边缘，不穿过任何单元内部)。我们将任意 $T\in S_h$ 的直径表示为 $h_T$ ，且令 $h=\max \left(h_T\right)$ 。根据拟一致网格的性质，我们有 $\frac{h_T}{h}\ge C,\forall T\in S_h$ 。从现在开始，我们将用C表示一个通用正常数，它不一定每次出现时都相同，但都与h无关。

对于 $\forall T\in S_h$ ，气泡函数定义为

$b_T\left(x\right)=\{\begin{array}{l}{\lambda }_1\left(x\right){\lambda }_2\left(x\right){\lambda }_3\left(x\right)\begin{array}{cc},& x\in T,\end{array}\\ 0,\begin{array}{ccc}\begin{array}{ccc}\begin{array}{ccc}& & \end{array}& & \end{array}& & x\notin T.\end{array}\end{array}$ (3.1)

其中 ${\lambda }_1\left(x\right)\begin{array}{c},\end{array}{\lambda }_2\left(x\right)\begin{array}{c},\end{array}{\lambda }_3\left(x\right)$ 分别表示 $\forall T\in S_h$ 三个顶点上的一次帽函数。显然， $b_T\left(x\right)\in H_0^1\left(T\right)\subset H_0^1\left(\Omega \right)$ 。

然后我们给出如下的速度有限元空间

$\begin{array}{l}{V}_h=\{v_h^S\in C^0\left({\Omega }_S\right),v_h^D\in C^0\left({\Omega }_D\right):{v_h|}_T\in P_2\left(T\right),\forall T\in S_h,v_h=0\text{on}{\Gamma }_S\begin{array}{c},\end{array}\\ v_h\cdot n_D=0\text{on}{\Gamma }_D\text{and}{v}_h^D\cdot n_D+v_h^S\cdot n_S=0\text{on}\Gamma \}\oplus {\left(\text{span}\left\{b_T\left(x\right)\begin{array}{c};\end{array}\forall T\in S_h\right\}\right)}^2\end{array}$ (3.2)

和压力的有限元空间

$Q_h=\left\{{q_h^S|}_{{\Omega }_S}\in C^0\left({\Omega }_S\right)\begin{array}{c},\end{array}{q_h^D|}_{{\Omega }_D}\in C^0\left({\Omega }_D\right):{q_h|}_T\in P_1\left(T\right)\begin{array}{c},\end{array}\forall T\in S_h\right\}\cap Q$ (3.3)

然后我们从(2.9)得到离散的混合问题：找到 $\left(u_h,p_h\right)\subset V_h\times Q_h$ 满足

$\{\begin{array}{l}\tilde{a}\left(u_h,v_h\right)+b\left(v_h,p_h\right)=\tilde{F}\left(v_h\right),\forall v_h\in V_h\\ b\left(u_h,q_h\right)=G\left(q_h\right)\begin{array}{cc}& \end{array}\begin{array}{cc}& \end{array}\forall q_h\in Q_h\begin{array}{c}.\end{array}\end{array}$(3.4)

4. 具有特殊形状的耦合界面的稳定性分析

在本节中，我们假设耦合界面 $\Gamma$ 满足以下几何性质(我们会在下一章移除该假设)：

假设4.1 $\Gamma$ 是一条穿过 $\Omega$ 的无分支的折线段，其两个端点是与 $\partial \Omega$ 的交点。

4.1. 一些引理

引理4.2 [14] ${\lambda }_i\left(x\right)\begin{array}{c},\end{array}{\lambda }_j\left(x\right)\begin{array}{c},\end{array}{\lambda }_k\left(x\right)$ 为 $\forall T\in S_h$ 的三个顶点上的一次帽函数，则

${\displaystyle {\int }_T{\lambda }_i^{\alpha }}{\lambda }_j^{\beta }{\lambda }_k^{\gamma } \text{d}x=\frac{\alpha !\beta !\gamma !}{\left(\alpha +\beta +\gamma +2\right)!}2\left|T\right|\begin{array}{c},\end{array}\begin{array}{cc}& \end{array}{\Vert {\lambda }_i^{\alpha }{\lambda }_j^{\beta }{\lambda }_k^{\gamma }\Vert }_1\le C$ (4.1)

其中 $\alpha ,\beta ,\gamma$ 是非负整数。另外，在该单元的边 $e_{ij}$ 上有

${\displaystyle {\int }_e{\lambda }_i^{\alpha }}{\lambda }_j^{\beta }\text{d}\sigma =\frac{\alpha !\beta !}{\left(\alpha +\beta +1\right)!}2\left|e\right|$ (4.2)

引理4.3 [14]设 $P_h:H^1\left(E\right)\to P_1\left(E\right)$ 是 $L^2$ 投影算子，则其满足

${\Vert \begin{array}{c}{P}_h\end{array}v\Vert }_{1,E}+h^{-1}{\Vert v-\begin{array}{c}{P}_h\end{array}v\Vert }_{0,E}+{\Vert v-\begin{array}{c}{P}_h\end{array}v\Vert }_{1,E}\le C{\Vert v\Vert }_{1,E}$ (4.3)

引理4.4 [15]存在常数 $C\begin{array}{c}>\end{array}0$ ，使得

${\Vert v\Vert }_{0\begin{array}{c},\end{array}\partial E}\le C{\Vert v\Vert }^{\frac{1}{2}}{}_{0,E}{\Vert v\Vert }^{\frac{1}{2}}{}_{1,E},\forall v\in H^1\left(E\right)$ (4.4)

4.2. 构造Fortin算子

根据假设4.1可得经过三角剖分后， $\Gamma$ 被分成m段，共有 $m+1$ 个节点。其中每段 $e_{i\Gamma }\left(i=1:m\right)$ 为单元T的一条边，每个节点 $i_{\Gamma }\left(i=1:m\right)$ 为单元T的一个顶点。并且在节点 $i_{\Gamma }$ 上的帽函数 ${\lambda }_{i\Gamma }$ 的支集最多含有 $\Gamma$ 上的两条边，对于相邻的两个帽函数 ${\lambda }_{i\Gamma }\begin{array}{c}，\end{array}{\lambda }_{j\Gamma }$ ，他们支集的交集只含有 $\Gamma$ 上的一条边。定义 ${\psi }_0=4{\lambda }_{0\Gamma }{\lambda }_{1\Gamma }{n}_1$ ， ${\psi }_m=4{\lambda }_{m-1\Gamma }{\lambda }_{m\Gamma }{n}_n$ ， ${\psi }_i={\lambda }_{i\Gamma }\left(2{\lambda }_{i\Gamma }-1\right)\left(n_{i-1}+n_i\right)$ ， $\left(i=1:m-1\right)$ 。其中， $n_i\left(i=1:m\right)$ 为边 $e_{i\Gamma }^{}$ 的单位法向量，指向 ${\Omega }_S$ 外侧。根据定义可得， $span\left({\psi }_i,i=0:m\right)\subset V_h$ ， ${Q_h|}_{\Gamma }\subset span\left({\lambda }_{i\Gamma },i=0:m\right)$ 。

定理4.5 存在算子 $r_h:V\to V_h$ 使得

$b\left(v-\begin{array}{c}{r}_h\end{array}v,q_h\right)=0\begin{array}{c},\end{array}\begin{array}{c}\end{array}\forall v\in V\begin{array}{c},\end{array}\begin{array}{c}\end{array}\forall q_h\in Q_h\begin{array}{c}.\end{array}$ (4.5)

证明：定义

$r_{h}v=P_{h}v+{\displaystyle \underset{T\in S_h}{\sum }{c}_T}{b}_T+{\displaystyle \underset{i=0}{\overset{n}{\sum }}{k}_i}{\psi }_i\begin{array}{c},\end{array}\begin{array}{c}\begin{array}{cc}& \end{array}\end{array}\forall v\in V\begin{array}{c}.\end{array}$ (4.6)

利用格林公式可得

$b\left(v-r_{h}v,q_h\right)={\displaystyle \underset{T\in S_h}{\sum }\left({\displaystyle {\int }_T\nabla }{q}_h\cdot \left(v-r_{h}v\right)\text{d}x\right)}-{\displaystyle {\int }_{\Gamma }\left(q_h^S-q_h^D\right)}\left(v_S-r_h{v}_S\right)\cdot n_S\text{d}t\begin{array}{c}.\end{array}$ (4.7)

由于 $\nabla q_h$ 在T上是常向量且 ${Q_h|}_{\Gamma }\subset span\left({\lambda }_{i\Gamma },i=0:m\right)$ ，所以若式(4.5)成立，只需

${\displaystyle {\int }_{\Gamma }{\lambda }_{i\Gamma }\left(v_S-r_h{v}_S\right)\cdot n_S}\text{d}\sigma =0\begin{array}{cc},& i=0:m\end{array},$ (4.8)

${\displaystyle {\int }_{\Omega }v-r_{h}v}=0\begin{array}{cc},& \forall T\in S_h\end{array}\begin{array}{c}.\end{array}$ (4.9)

将(4.6)代入(4.8)可得

${\displaystyle \underset{j=0}{\overset{m}{\sum }}{k}_j{\displaystyle {\int }_{\Gamma }{\lambda }_{i\Gamma }{\psi }_j\cdot n_S}\text{d}\sigma ={\displaystyle {\int }_{\Gamma }{\lambda }_{i\Gamma }\left(v_S-P_h{v}_S\right)\cdot n_S}\text{d}\sigma \begin{array}{cc},& i=0:m\end{array}}.$ (4.10)

利用引理4.2可以将(4.10)转化为下述线性方程组

$A_{\Gamma }K=F_{\Gamma }$ (4.11)

其中

$\begin{array}{l}{A}_{\Gamma }=diag\left(\frac{\left|e_{1\Gamma }\right|}{3},\frac{\left(\left|e_{1\Gamma }\right|+\left|e_{2\Gamma }\right|\right)\left(1+n_1\cdot n_2\right)}{6},\cdots ,\frac{\left(\left|e_{i\Gamma }\right|+\left|e_{i+1\Gamma }\right|\right)\left(1+n_i\cdot n_{i+1}\right)}{6}\cdots \frac{\left(\left|e_{m-1\Gamma }\right|+\left|e_{m\Gamma }\right|\right)\left(1+n_{m-1}\cdot n_m\right)}{6},\frac{\left|e_{m\Gamma }\right|}{3}\right),\\ K={\left(\begin{array}{cc}{k}_0,& \begin{array}{ccc}{k}_1,& \cdots & \begin{array}{cc}{k}_i& \begin{array}{cc}\cdots & \begin{array}{cc}{k}_{m-1},& k_m\end{array}\end{array}\end{array}\end{array}\end{array}\right)}^{\text{T}},\\ F_{\Gamma }=(\begin{array}{cc}{\displaystyle {\int }_{\Gamma }{\lambda }_{0\Gamma }\left(v_S-P_h{v}_S\right)\cdot n_S}\text{d}\sigma ,& \begin{array}{ccc}{\displaystyle {\int }_{\Gamma }\left({\lambda }_{1\Gamma }-{\lambda }_{0\Gamma }\right)\left(v_S-P_h{v}_S\right)\cdot n_S}\text{d}\sigma ,& \cdots & \begin{array}{cc}{\displaystyle {\int }_{\Gamma }{\lambda }_{i\Gamma }\left(v_S-P_h{v}_S\right)\cdot n_S}\text{d}\sigma & \cdots \end{array}\end{array}\end{array}\\ {\begin{array}{ccc}\begin{array}{cc}& \end{array}& \cdots & \begin{array}{cc}{\displaystyle {\int }_{\Gamma }\left({\lambda }_{m-1\Gamma }-{\lambda }_{m\Gamma }\right)\left(v_S-P_h{v}_S\right)\cdot n_S}\text{d}\sigma ,& {\displaystyle {\int }_{\Gamma }{\lambda }_{m\Gamma }\left(v_S-P_h{v}_S\right)\cdot n_S}\text{d}\sigma )\end{array}\end{array}}^{\text{T}}.\end{array}$

根据拟一致剖分的性质可得，存在与h无关的常数C，使得 $1+n_i\cdot n_{i+1}\ge C>0$ 。因此，方程组(4.11)存在唯一解。

接着利用引理4.2并使用同样的方法将(4.6)代入(4.9)可得

$c_T=\frac{60{\displaystyle {\int }_T\left(v-P_{h}v-{\displaystyle \underset{i=0}{\overset{m}{\sum }}{k}_i}{\psi }_i\right)\text{d}x}}{\left|T\right|} ,\forall T\in S_h.$ (4.12)

由此可得 $c_T$ 存在唯一解，定理4.5得证。

定理4.6 对于 $\forall v\in H_0^1\left(\Omega \right)$ ，算子 $r_h$ 满足

${\Vert r_{h}v\Vert }_V\le C{\Vert v\Vert }_1.$ (4.13)

证明：利用holder不等式， $1+n_i\cdot n_{i+1}\ge C>0$ ，引理4.3以及引理4.4可由(4.11)得

$\begin{array}{c}{\displaystyle \underset{i=0}{\overset{m}{\sum }}{\left|k_i\right|}^2}\le {\displaystyle \underset{i=0}{\overset{m}{\sum }}C{h}^{-2}{\left({\displaystyle {\int }_{\Gamma }{\lambda }_{i\Gamma }}\left(v_S-P_h{v}_S\right)\cdot n_S\text{d}\sigma \right)}^2}\\ \le {\displaystyle \underset{i=0}{\overset{m}{\sum }}C{h}^{-2}\left({\displaystyle {\int }_{\Gamma }{\lambda }_{i\Gamma }^2}\text{d}\sigma {\Vert v_S-P_h{v}_S\Vert }_{0,e_{i\Gamma }}^2\right)}\\ \le C{h}^{-1}{\Vert v_S-P_h{v}_S\Vert }_{0,\Gamma }^2\\ \le C{h}^{-1}{\Vert v_S-P_h{v}_S\Vert }_{0,{\Omega }_S}{\Vert v_S-P_h{v}_S\Vert }_{1,{\Omega }_S}\\ \le C{\Vert v\Vert }_1^2.\end{array}$

利用相同的方法可由(4.11)得

$\begin{array}{c}{\displaystyle \underset{T\in S_h}{\sum }{\Vert c_T\Vert }^2}\le C{h}^{-4}{\displaystyle \underset{T\in S_h}{\sum }{\left({\displaystyle {\int }_T\left(v-P_{h}v-{\displaystyle \underset{i=0}{\overset{m}{\sum }}{k}_i}{\psi }_i\right)\text{d}x}\right)}^2}\\ \le C{h}^{-4}{\displaystyle \underset{T\in S_h}{\sum }{\left({\displaystyle {\int }_T\begin{array}{c}1\end{array}}\begin{array}{c}\text{d}x\end{array}\right)}^2}{\Vert v-P_{h}v-{\displaystyle \underset{i=0}{\overset{m}{\sum }}{k}_i}{\psi }_i\Vert }_{0,T}^2\\ \le C{h}^{-2}{\Vert v-P_{h}v-{\displaystyle \underset{i=0}{\overset{m}{\sum }}{k}_i}{\psi }_i\Vert }_0^2\\ \le C{\left({\Vert v\Vert }_1+h^{-1}{\Vert {\displaystyle \underset{i=0}{\overset{m}{\sum }}{k}_i}{\psi }_i\Vert }_0\right)}^2.\end{array}$

结合引理4.2可得

$\begin{array}{c}{\Vert r_{h}v\Vert }_V\le {\Vert r_{h}v\Vert }_1\le {\Vert P_{h}v\Vert }_1+{\Vert {\displaystyle \underset{T\in S_h}{\sum }{c}_T{b}_T}\Vert }_1+{\Vert {\displaystyle \underset{i=0}{\overset{m}{\sum }}{k}_i{\psi }_i}\Vert }_1\\ ={\Vert P_{h}v\Vert }_1+{\left({\displaystyle \underset{T\in S_h}{\sum }{\Vert c_T{b}_T\Vert }_{1,T}^2}\right)}^{\frac{1}{2}}+{\left({\displaystyle \underset{T\in S_h}{\sum }{\Vert {\displaystyle \underset{i\Gamma \in T}{\sum }{k}_i{\psi }_i}\Vert }_{1,T}^2}\right)}^{\frac{1}{2}}\\ \le C{\Vert v\Vert }_1+C{\left({\displaystyle \underset{T\in S_h}{\sum }{\Vert c_T\Vert }^2}\right)}^{\frac{1}{2}}+C{\left({\displaystyle \underset{i=0}{\overset{m}{\sum }}{\Vert k_i{\psi }_i\Vert }_1^2}\right)}^{\frac{1}{2}}\\ \le C{\Vert v\Vert }_1+C{\left({\displaystyle \underset{T\in S_h}{\sum }{\Vert c_T\Vert }^2}\right)}^{\frac{1}{2}}+C{\left({\displaystyle \underset{i=0}{\overset{m}{\sum }}{\left|k_i\right|}^2}\right)}^{\frac{1}{2}}\\ \le C{\Vert v\Vert }_1+C{h}^{-1}{\Vert {\displaystyle \underset{i=0}{\overset{m}{\sum }}{k}_i}{\psi }_i\Vert }_0+C{\left({\displaystyle \underset{i=0}{\overset{m}{\sum }}{\left|k_i\right|}^2}\right)}^{\frac{1}{2}}\\ \le C{\Vert v\Vert }_1.\end{array}$

定理4.7 存在常数 $C>0$ ，使得

$\underset{0\ne v_h\in V_h}{\sup }\frac{b\left(v_h,q_h\right)}{{\Vert v_h\Vert }_V}\ge C{\Vert q_h\Vert }_Q\begin{array}{cc},& \end{array}\forall q_h\in Q_h.$ (4.14)

证明：根据[16]中的经典结论，存在 $v\in H_0^1\left(\Omega \right)$ 使得 $div\begin{array}{c}v\end{array}=-q_h$ 并且 ${\Vert v\Vert }_1\le C{\Vert q_h\Vert }_0$ ， $\forall q_h\in Q_h$ 。结合定理4.5以及定理4.6，可得：对于 $\forall q_h\in Q_h$ ，存在 $v\in H_0^1\left(\Omega \right)$ 使得

$\underset{0\ne v_h\in V_h}{\sup }\frac{b\left(v_h,q_h\right)}{{\Vert v_h\Vert }_V}\ge \frac{b\left(r_{h}v,q_h\right)}{{\Vert r_{h}v\Vert }_V}=\frac{b\left(v,q_h\right)}{{\Vert r_{h}v\Vert }_V}=\frac{{\Vert q_h\Vert }_0^2}{{\Vert r_{h}v\Vert }_V}\ge C\frac{{\Vert q_h\Vert }_0^2}{{\Vert v\Vert }_1}\ge C{\Vert q_h\Vert }_0.$

5. 一般形状的耦合界面的稳定性分析

在本节中，我们考虑一般情况下的稳定性分析，即 $\Gamma$ 是任意折线段，将 $\Omega$ 分成多个不连通区域。

对此，我们可以将 $\Gamma$ 分成M段(与h无关)仅有端点相交的无分支的折线段 ${\Gamma }_j\left(j=1:M\right)$ ，满足 $\Gamma ={\displaystyle \underset{j=1}{\overset{M}{\cup }}{\Gamma }_j}$ 且每段 ${\Gamma }_j$ 至多有两个端点落在 $\partial \Omega$ 上。根据第四节中的方法，我们可以构造如下的Fortin算子

$r_{h}v=P_{h}v+{\displaystyle \underset{T\in S_h}{\sum }{c}_T}{b}_T+{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=0}{\overset{n}{\sum }}{k}_{ij}}{\psi }_{ij}}\begin{array}{c},\end{array}\begin{array}{c}\begin{array}{cc}& \end{array}\end{array}\forall v\in V\begin{array}{c}.\end{array}$ (5.1)

并运用相同的方法，我们依然能够证明定理4.7在一般情况下成立。

定理5.1 若定理2.1和定理4.7成立，且问题(2.9)的解 $\left(u,p\right)$ 满足 $u_S\in H^3\left({\Omega }_S\right)$ ， $u_D\in H^3\left({\Omega }_D\right)$ $p_S\in H^2\left({\Omega }_S\right)$ ， $p_D\in H^2\left({\Omega }_D\right)$ 则离散问题(3.4)的解 $\left(u_h,p_h\right)$ 满足如下误差估计

${\Vert u-u_h\Vert }_V+{\Vert p-p_h\Vert }_Q\le C{h}^2\left({\Vert u\Vert }_{3,{\Omega }_S}+{\Vert u\Vert }_{3,{\Omega }_D}+{\Vert p\Vert }_{2,{\Omega }_S}+{\Vert p\Vert }_{2,{\Omega }_D}\right).$ (5.2)

6. 数值实验

在本节中，我们将使用三个具有不同形状耦合界面的示例来测试理论收敛结果，其中前两个例子中 $\Gamma$ 穿过 $\Omega$ ，第三个例子中 $\Gamma$ 是 $\Omega$ 中的闭环。我们定义以下误差符号

$e_0\left(u_S\right)={\Vert u_S-u_h^S\Vert }_{0,{\Omega }_S}{e}_1\left(u_S\right)={\left|u_S-u_h^S\right|}_{1,{\Omega }_S}{e}_0\left(u_D\right)={\Vert u_D-u_h^D\Vert }_{0,{\Omega }_D}$

$e_0\left(\begin{array}{c}div\end{array}{u}_D\right)={\Vert \begin{array}{c}div\end{array}{u}_D-\begin{array}{c}div\end{array}{u}_h^D\Vert }_{0,{\Omega }_D}{e}_0\left(p\right)={\Vert p-p_h\Vert }_0$

收敛阶符号定义为

$r_i=\frac{\log \left(\frac{e_i}{{e^{\prime }}_i}\right)}{\log \left(\frac{h}{h^{\prime }}\right)},i=0,2$

6.1. 数值实验例1

在第一个数值算例中，我们设 ${\Omega }_S=\left(0,1\right)\times \left(0,1\right)$ ， ${\Omega }_D=\left(1,2\right)\times \left(0,1\right)$ ， $\Gamma =1\times \left(0,1\right)$ 。解析解为：

$\{\begin{array}{l}{u}_S=\left(\begin{array}{c}\sin \left(\pi x\right)y(y-1)\\ x\left(2x-3\right)y\left(y-1\right)\end{array}\right),p_S=-\pi x^{2}y\left(2y-1\right)+\frac{5\pi }{12}\\ u_D=\left(\begin{array}{c}\frac{1}{2}\sin \left(\pi x\right)\cos \left(\pi y\right)\\ \frac{1}{2}\cos \left(\pi x\right)\sin \left(\pi y\right)\end{array}\right),p_D=-\pi x^{2}y+\frac{5\pi }{12}\end{array}$ (6.1)

上述解析解满足Stokes-Darcy耦合问题(2.1)~(2.3)，其误差和相应的收敛阶如表1所示。

Table 1. The errors and the rates of convergence of example 1

表1. 例1的误差和相应的收敛阶

(a)
h	$e_0\left(u_S\right)$	$r_0\left(u_S\right)$	$e_1\left(u_S\right)$	$r_1\left(u_S\right)$	$e_0\left(u_D\right)$
1/4	6.869e−04	-	1.661e−02	-	6.272e−02
1/8	7.394e−05	3.21	4.036e−03	2.04	1.746e−02
1/16	8.176e−06	3.17	1.001e−03	2.01	4.555e−03
1/32	9.496e−07	3.05	2.497e−04	2.00	1.156e−03
(b)
h	$r_0\left(u_D\right)$	$e_0\left(\begin{array}{c}div\end{array}{u}_D\right)$	$r_0\left(\begin{array}{c}div\end{array}{u}_D\right)$	$e_0\left(p\right)$	$r_0\left(p\right)$
1/4	-	6.077e−02	-	6.945e−02	-
1/8	1.84	1.508e−02	2.01	1.557e−02	2.14
1/16	1.94	3.764e−03	2.00	3.831e−03	2.04
1/32	1.98	9.406e−04	2.00	9.501e−04	2.01

6.2. 数值实验例2

Table 2. The errors and the rates of convergence of example 2

表2. 例2的误差和相应的收敛阶

(a)
h	$e_0\left(u_S\right)$	$r_0\left(u_S\right)$	$e_1\left(u_S\right)$	$r_1\left(u_S\right)$	$e_0\left(u_D\right)$
1/4	3.643e−01	-	9.974e−01	-	9.155e−01
1/8	4.219e−02	3.11	2.459e−01	2.02	2.402e−01
1/16	5.025e−03	3.07	6.148e−02	2.00	6.349e−02
1/32	6.151e−03	3.03	1.526e−02	2.01	1.632e−02
(b)
h	$r_0\left(u_D\right)$	$e_0\left(\begin{array}{c}div\end{array}{u}_D\right)$	$r_0\left(\begin{array}{c}div\end{array}{u}_D\right)$	$e_0\left(p\right)$	$r_0\left(p\right)$
1/4	-	5.638e−00	-	3.770e−01	-
1/8	1.93	1.371e−00	2.04	9.041e−02	2.06
1/16	1.92	3.427e−01	2.00	2.229e−02	2.02
1/32	1.96	8.568e−02	2.00	5.534e−03	2.01

在第二个数值算例中，我们设 ${\Omega }_S=\left(0,\pi \right)\times \left(0,\pi \right)$ ， ${\Omega }_D=\left(0,\pi \right)\times \left(-\pi ,0\right)$ ， $\Gamma =\left(0,\pi \right)\times 0$ 。解析解为：

$\{\begin{array}{l}{u}_S=\left(\begin{array}{c}2\sin y\cos y\cos x\\ \left({\sin }^{2}y-2\right)\sin x\end{array}\right),p_S=\sin x\sin y\\ u_D=\left(\begin{array}{c}-\left({\text{e}}^y-{\text{e}}^{-y}\right)\cos x\\ -\left({\text{e}}^y+{\text{e}}^{-y}\right)\sin x\end{array}\right),p_D=\left({\text{e}}^y-{\text{e}}^{-y}\right)\sin x\end{array}$ (6.2)

同样，上述解析解满足Stokes-Darcy耦合问题(2.1)~(2.3)，其误差和相应的收敛阶如表2所示。

6.3. 数值实验例3

Table 3. The errors and the rates of convergence of example 3

表3. 例3的误差和相应的收敛阶

(a)
h	$e_0\left(u_S\right)$	$r_0\left(u_S\right)$	$e_1\left(u_S\right)$	$r_1\left(u_S\right)$	$e_0\left(u_D\right)$
1/4	5.372e−01	-	1.282e−00	-	1.168e−00
1/8	3.838e−02	3.80	2.840e−01	2.17	3.370e−01
1/16	2.965e−03	3.69	6.540e−02	2.11	9.155e−02
1/32	2.791e−04	3.40	1.601e−02	2.03	2.436e−02
(b)
h	$r_0\left(u_D\right)$	$e_0\left(\begin{array}{c}div\end{array}{u}_D\right)$	$r_0\left(\begin{array}{c}div\end{array}{u}_D\right)$	$e_0\left(p\right)$	$r_0\left(p\right)$
1/4	-	1.013e−00	-	1.282e−00	-
1/8	1.80	2.563e−01	1.90	3.889e−01	1.72
1/16	1.88	5.553e−02	2.20	9.672e−02	2.01
1/32	1.91	1.371e−02	2.01	2.393e−02	2.00

在第三个数值算例中，我们设 ${\Omega }_S=\left(-\frac{5\pi }{4},\frac{5\pi }{4}\right)\times \left(-\frac{5\pi }{4},\frac{5\pi }{4}\right)$ ， ${\Omega }_S=\left(-\frac{5\pi }{2},\frac{5\pi }{2}\right)\times \left(-\frac{5\pi }{2},\frac{5\pi }{2}\right)$ 。解析解为：

$\{\begin{array}{l}{u}_S=\left(\begin{array}{c}\sin \left(x\right)\cos \left(y\right)\\ \cos \left(x\right)\sin \left(y\right)\end{array}\right),p_S=\frac{1}{2}\cos \left(x\right)\cos \left(y\right)\\ u_D=\left(\begin{array}{c}\frac{32\sqrt{2}}{375{\pi }^3}\left(x+\frac{5\pi }{2}\right)\left(x-\frac{5\pi }{2}\right)\cos \left(y\right)\\ \frac{32\sqrt{2}}{375{\pi }^3}\left(y+\frac{5\pi }{2}\right)\left(y-\frac{5\pi }{2}\right)\cos \left(x\right)\end{array}\right),p_D=-\frac{1}{2}\cos \left(x\right)\cos \left(y\right)\end{array}$ (6.3)

第三个例子的解析解同样满足Stokes-Darcy耦合问题(2.1)~(2.3)。在该例子中， ${\Omega }_S$ 被包裹在 ${\Omega }_D$ 中，并且 $\Gamma =\partial {\Omega }_S$ 是 $\Omega$ 中的一个闭环。 $u_S$ 在 $\partial \Omega$ 上的值完全由耦合界面条件(2.3)决定，但如同理论分析一样，我们仍能得到良好的收敛性，如表3所示。

7. 结论

本文给出了平面区域中Stokes-Darcy耦合问题的具有二阶收敛精度的统一近似格式。

相比于在两个区域上使用不同有限元空间的数值格式，虽然在两个区域上的网格剖分可以不相同，网格剖分的选择更加灵活。但两个不同有限元空间的函数在耦合界面 $\Gamma$ 上的迹常常无法兼容，需要添加拉格朗日乘子，以弱形式来满足耦合界面条件 $v_h^D\cdot n_D+v_h^S\cdot n_S=0$ 。而本文提出的统一有限元格式，由于在两个区域上的有限元空间是相同的，所以当我们对整个区域 $\Omega$ 进行统一剖分后，只需耦合界面 $\Gamma$ 上的基函数进行适当的合并，可以自然地满足耦合界面条件 $v_h^D\cdot n_D+v_h^S\cdot n_S=0$ 。

NOTES

^*通讯作者。