近年来, 不可压缩流与多孔介质流之间的耦合问题受到人们的广泛关注. 该耦合问题能够有效地模拟地表水与地下水之间的运移问题, 石油开采问题, 工业中与流体过滤有关的问题以及生物流体力学中的血液循环问题等. Navier-Stokes/Darcy 模型和 Stokes/Darcy 模型是研究该耦合问题的重要模型, 研究者已经针对这两个模型进行了深入的研究. 本文重点讨论 Stokes/Darcy 模型。

目前, 针对定常 Stokes/Darcy 模型的研究已经相对成熟, 其数值方法包括有限元法^{[1⇓⇓⇓-5]}, 间断伽辽金法^[6⇓-8], 界面松弛法^[9-10], 拉格朗日乘子法^[11-12], 区域分解法^[13⇓-15], 两重网格或多重网格法^{[16⇓⇓-19]}等. 现在, 学者们越来越重视非定常 Stokes/Darcy 模型的研究, 并通过不同的方法获得其时间半离散格式的相关结论. 如整体相同时间步长方法^[20⇓-22], 不同区域不同时间步长方法^{[23⇓⇓-26]}. 然而, 现存的方法仍然存在一些不足. 例如, 某类数值方法需要消耗大量的计算时间才能得到相对较高的数值精度; 某类算法需要对时间步长给予特定的限制条件其稳定性和误差估计才能得到证明; 某类算法的稳定性需要增加一个稳定化项才能得到证明. 因此, 研究更高效, 适应性更强的算法是十分必要的.

多步法是求解常微分方程初值问题的基本数值解法. 多步法的主要思想是保留和使用前面多步数值来预测当前时刻数值使算法尽可能获得较高精度. 本文的研究基于一阶线性多步法的 $\theta$ -格式, 其中 $\theta$ 是满足 $0<\theta<1/2$ 的参数. 特别地, 当 $\theta=0$ 时, 该算法为向后欧拉方法. 与向后欧拉方法相比较, 线性多步法更具一般性, 其适应性更强. Time filter 算法^[27⇓-29]在程序中易实现, 并且能够有效地提高算法的收敛阶和计算效率. 该算法是模块化的, 只需增加额外的两行代码就能够将线性多步法的精确度从一阶提高到二阶. 结合线性多步法和 time filter 算法的优势, 本文提出线性多步法加 time filter 算法, 并理论分析了耦合和解耦的线性多步法加 time filter 算法的稳定性和误差估计. 最后, 本文进行了三个数值实验. 第一个实验模拟了耦合和解耦算法的有效性. 第二个实验体现出耦合和解耦算法的收敛阶从一阶提高到二阶, 并且通过比较可以得出结论: 解耦算法的计算效率相对较高. 第三个实验说明相比于向后欧拉方法, 线性多步法几乎在每个时刻的误差相对较小, 并且随着$\theta$的不断增大, 线性多步法的误差会减小. 类似地, 与向后欧拉方法加 time filter 算法相比, 线性多步法加 time filter 算法也具有同样的特点.

本文的其余部分组织如下: 第 2 节介绍具有 Beavers-Joseph Saffman (BJS) 界面条件的 Stokes/Darcy 模型, 并给出其弱形式; 第 3 节提出耦合和解耦的线性多步法加 time filter 算法并分析其稳定性; 第 4 节分析耦合和解耦算法的误差估计; 第 5 节利用数值试验进一步佐证理论分析的结果.

在有界光滑区域 $\Omega\subset{\bf R}^d$, $d=2,3$ 内, 我们讨论自由流体流动和多孔介质流动的耦合问题. 其中 $\Omega=\Omega_f\cup\Gamma\cup\Omega_p$, $\Omega_f$ 和 $\Omega_p$ 分别表示两个不相交, 连续且有界的自由流体区域和多孔介质区域. $\Gamma=\overline{\Omega}_f\cap\overline{\Omega}_p$ 是两区域之间的交界面, $\Gamma_f=\partial\Omega_f\cap\partial\Omega$, $\Gamma_p=\partial\Omega_p\cap\partial\Omega$. 为了方便起见, 假设 $\partial\Omega_f$, $\partial\Omega_p$ 都足够光滑, ${\bf n}_f$ 和 ${\bf n}_p$ 分别表示 $\partial\Omega_f$ 和 $\partial\Omega_p$ 上的单位外法向向量. 如图 1 所示.

本文研究由 Stokes 方程和 Darcy 方程组成的 Stokes/Darcy 模型. 设 $T>0$ 是一个有限的时间. 自由流体区域 $\Omega_f$ 流动由 Stokes 方程控制: 求流体速度 ${\bf u}_f:\Omega_f\times[T]\rightarrow{\bf R}^d$ 和压力 $p_f:\Omega_f\times[T]\rightarrow{\bf R}$ 满足

其中

分别是应力张量和应变张量, $\nu>0$ 是粘性系数, ${\bf g}_f({\bf x},t)$ 是外力项.

多孔介质区域 $\Omega_p$ 流动由以下方程控制

其中, 第一个方程是饱和流动模型, $S$ 是比质量储能系数, $\phi_p=z+\frac{p_p}{\rho g}$ 是测压水头, $p_p$ 表示流体在 $\Omega_p$ 中的压力, $z$ 是从参考平面算起的高度, $\rho$ 和 $g$ 分别表示密度和重力加速度. 第二个方程是 Darcy 定律, ${\bf u}_p$ 是多孔介质区域 $\Omega_p$ 的流体速度, 与 $\phi_p$ 的梯度成正比关系. ${\bf K}$ 表示水利传导率张量, 它是一个最小特征值大于零的对称正定的矩阵, 且允许在空间中变化. ${\bf g}_p({\bf x},t)$ 是源项.

将方程(2.4)代入方程(2.3)可得 Darcy 方程: 求测压水头 $\phi_p:\Omega_p\times[T]\rightarrow{\bf R}^d$ 满足

下面考虑 Stokes/Darcy 模型满足以下初值条件

以及齐次 Dirichlet 边界条件

交界面条件是 Stokes/Darcy 模型的重要组成部分, 包括质量守恒条件, 法向力的平衡条件和 Beavers-Joseph-Saffman(BJS) 条件. 具体表示如下: 在 $\Gamma$ 上满足

其中 ${\boldsymbol{\tau}}_i$, $i=1,2,\cdots,d-1$, 是沿着 $\Gamma$ 的正交切向单位向量, $\alpha$ 是一个确定的实验参数, ${\bf \Pi}$ 表示渗透率, 与水利传导率张量存在关系: ${\bf K}=\frac{{\bf \Pi}g}{\nu}$.

首先定义 Hilbert 空间

区域 $D$ 中, $(\cdot,\cdot)_D$ 表示 $D=\Omega_f$ 或 $\Omega_p$ 中的内积. 并对空间 ${\bf U}$ 配备以下范数: $\forall$$\underline{\bf u}=({\bf u}_f,\phi_p)\in {\bf U}$,有

用 $\Vert\cdot\Vert_{{\bf H}_f/{\bf H}_p}$, $\Vert\cdot\Vert_{f/p}$ 和 $\Vert\cdot\Vert_{\Gamma}$ 分别表示范数 ${\bf H}^1(\Omega_{f/p})$, $L^2(\Omega_{f/p})$, $L^2(\Gamma)$. 同时定义下列范数

对函数 $v(x,t)$ 定义范数

其次, 非定常Stokes/Darcy 耦合模型(2.1)-(2.10)式的弱形式描述如下: 对于 ${\bf g}_f\in L^2(0,T; {\bf L}^2(\Omega_f))$, ${\bf g}_p\in L^2(0,T; {\bf L}^2(\Omega_p))$, 求 $\underline{\bf u}=(u_f, \phi_p) \in (L^2(0,T; {\bf H}_f)\cap L^\infty(0,T; L^2(\Omega_f)) \times L^2(0,T; {\bf H}_p)\cap L^\infty(0,T; L^2(\Omega_p)) )$ 和 $p_f\in L^2(0,T; Q_f)$ 使得 $\forall$$(\underline{\bf v},q_f) \in{\bf U}\times Q_f$ 满足

其中

${\bf U}'$ 是 ${\bf U}$ 的对偶空间, $P_{\tau}(\cdot)$ 是局部切平面上的投影, 表示为 $P_{\tau}({\bf v}_f)={\bf v}_f-({\bf v}_f\cdot{\bf n}_f){\bf n}_f$.

双线性形式 $a(\cdot,\cdot)$ 满足连续性和强制性: $\forall$$\underline{\bf u}, \underline{\bf v}\in {\bf U},$${\bf u}_f\in {\bf H}_f,$$\phi_p \in {\bf H}_p,$有

$a_\Gamma(\cdot,\cdot)$ 满足反对称性

下面给出Poincar$\acute{e}$, Trace, Sobolev 不等式: 存在常数 $C_p$, $C_t$ 和 $C_s$ 仅仅依赖于区域 $\Omega_f$, $\tilde{C}_p$, $\tilde{C}_t$ 和 $\tilde{C}_s$ 仅仅依赖于区域 $\Omega_p$, 使得 $\forall$${\bf v}_f\in{\bf H}_f$ 和 $\psi_p\in{\bf H}_p$,有

定常的耦合Stokes/Darcy 模型的适定性已被证明. 这里, 假定在非稳态的情况下也同样成立. 本文主要研究其数值解.

对任意给定的小参数 $h>0$, 构造 $\Omega$, $\Omega_f$ 和 $\Omega_p$ 的正则三角剖分 ${\cal T}_h$, ${\cal T}_{fh}$ 和 ${\cal T}_{ph}$. 为了简便, 设 $\Omega_f$ 和 $\Omega_p$ 是足够光滑的区域. 通过 MINI 元 (P1b-P1) 离散的 ${\bf H}_{fh}\subset{\bf H}_f$, $Q_{fh}\subset Q_f$ 和通过线性拉格朗日元 (P1) 离散的 ${\bf H}_{ph}\subset {\bf H}_p$ 是有限元空间, 使得空间对 $({\bf H}_{fh},Q_{fh})$ 满足离散的 LBB 条件

并定义 ${\bf u}_h={\bf H}_{fh}\times {\bf H}_{ph}$.

定义线性投影算子(见文献[31]): $\forall$$\underline{\bf v}_h\in{\bf u}_h$, $q_{fh}\in Q_{fh}$ 和 $t\in(0,T]$, $P_h:(\underline{\bf u}(t),p_f(t)) \in {\bf U}\times Q_f\to(P_h^{\underline{\bf u}}\underline{\bf u}(t), P_h^{p}p_f(t)) \in {\bf u}_h\times Q_{fh}$ 满足

假设 $(\underline{\bf u}(t),p_f(t))$ 足够光滑, $(\underline{\bf u}(t),p_f(t))$ 的投影 $(P_h^{\underline{\bf u}}\underline{\bf u}(t)$, $P_h^{p}p_f(t))$ 满足性质

Inverse 不等式: 存在常数 $C_I$, $\tilde{C}_I$ 依赖于有限元网格中的角度, 使得 $\forall$${\bf v}_{fh}\in{\bf H}_{fh}$, $\psi_{ph}\in {\bf H}_{ph}$,

理论分析时将用到下面两个引理.

${\bf引理3.1}$ (文献[30,第4章引理 3.3]) 设 $\delta =\beta_2-\frac{\alpha_2}{2}>0$. 系数 $\alpha_i$ 和 $\beta_i$ 满足下面关系

${\bf引理3.2}$ (文献[引理 1]) $\forall$${\bf v}_f\in {\bf H}_f$, $\phi_p\in{\bf H}_p$, 存在 $C_k>0$, 使得 $\forall$$\varepsilon>0$,有

另外, $\forall$${\bf v}_{fh}\in {\bf H}_{fh}$, $\phi_{ph}\in{\bf H}_{ph}$, 存在 $\tilde{C}_k>0$, 使得 $\forall$$\tilde{\varepsilon}>0$, 有

这里引入相关的记号

本节其余部分将给出非定常 Stokes/Darcy 模型耦合和解耦的线性多步法加 time filter 算法以及其稳定性. 将时间段 $[T]$ 进行均匀剖分, 并令 $t_m=m\Delta t$, $m=0,1,\cdots,N$, 其中 $\Delta t=\frac{T}{N}$ 是时间步长. $(\underline{\bf u}_h^{m+1},p_{fh}^{m+1})=({\bf u}_{fh}^{m+1},p_{fh}^{m+1},\phi_{ph}^{m+1})$ 表示 $(\underline{\bf u}_h(t_{m+1}), p_{fh}(t_{m+1}))=({\bf u}_{fh}(t_{m+1}),p_{fh}(t_{m+1}),\phi_{ph}(t_{m+1}))$ 的全离散格式的数值解.

耦合的线性多步法加 time filter 算法描述如下.

$\blacktriangledown$ 线性多步法 (一阶收敛)

给定 $({\underline{\bf u}}_h^0,{p}_{fh}^0)$ 和 $({\underline{\bf u}}_h^1,{p}_{fh}^1)$, 求 $\hat{\underline{\bf u}}_h^{m+1}=(\hat{\bf u}_{fh}^{m+1}$, $\hat{\phi}_{ph}^{m+1})\in{\bf u}_h$, $\hat{p}_{fh}^{m+1}\in Q_{fh}$, $m=1,2,$$\cdots,N-1$, 使得 $\forall$$\underline{\bf v}_h\in{\bf u}_h$ 和 $q_{fh}\in Q_{fh}$,

$\blacktriangledown$ Time Filter 算法 (二阶收敛)

通过 time filter 算法更新线性多步法的解 $(\hat{\underline{\bf u}}_h^{m+1},\hat{p}_{fh}^{m+1})$,

${\bf注3.1}$ 无论压力 $\hat{p}_{fh}^{m+1}$ 是否经 time filter 更新, 最终的结论几乎一致.

下面给出耦合的线性多步法加 time filter 算法的稳定性分析.

${\bf定理3.1}$ (耦合算法的稳定性) 设 $\underline{\bf u}_h^{m+1}$ 是线性多步法加 time filter 算法的解. 对于 $N\geq 2$, 有

其中 $C$ 是一个常数, 与 $h$, $\Delta t$ 以及其他参数无关, 并且在不同的情况下可能取不同的值.

${\bf证}$ 由 (3.8) 式可得

将上述两式代入 (3.7) 式, $\forall\ \underline{\bf v}_h\in {\bf u}_h$, $q_{fh}\in Q_{fh}$,有

假设 $\underline{\bf v}_h=2\Delta t B(\underline{\bf u}_h^{m+1},\theta)$, $q_{fh}=2\Delta t B(p_{fh}^{m+1},\theta)$, 下面对 (3.8) 式的每一项进行分析. 首先, 由于$\frac{(1-\theta)(3-2\theta)}{2}-\frac{1}{2}\frac{3-2\theta}{2}>0$, 因此借助引理3.1, 有

其次, 利用双线性形式 $a(\cdot,\cdot)$ 的强制性可得

由无散度特性可知, 左边第三项等于 0. 最后, 根据 Young 不等式和 Cauchy Schwarz 不等式, 有

结合 (3.10)-(3.12) 式, 并将 (3.9) 式从 $m=1$ 直到 $m=N-1$ 求和, 整理得

由于

因此

从而得到耦合算法的稳定性结论.证毕.

解耦的线性多步法加 time filter 算法表述如下.

$\blacktriangledown$ 线性多步法 (一阶收敛)

给定 $({\bf u}_{fh}^0,p_{fh}^0)$ 和 $({\bf u}_{fh}^1,p_{fh}^1)$, 求 $ (\hat{\bf u}_{fh}^{m+1},\hat{p}_{fh}^{m+1})\in({\bf H}_{fh},Q_{fh})$, $m=1,2,\cdots,N-1$, 使得 $\forall\ {\bf v}_{fh} \in{\bf H} _{fh}$, $q_{fh}\in Q_{fh}$,有

给定 $\phi_{ph}^0$ 和 $\phi_{ph}^1$, 求 $\hat{\phi}_{ph}^{m+1}\in{\bf H}_{ph}$, $m=1,2,\cdots,N-1$, 使得 $\forall$$\psi_{ph} \in{\bf H}_{ph}$,有

$\blacktriangledown$ Time Filter 算法 (二阶收敛)

通过 time filter 算法更新线性多步法的解 $(\hat{\bf u}_{fh}^{m+1},\hat{p}_{fh}^{m+1},\hat{\phi}_{ph}^{m+1})$,

${\bf注3.2}$ 交界面项利用二阶外推方法近似 $\hat{\bf u}_{fh}^{m+1}$ 和 $\hat{\phi}_{ph}^{m+1}$, 即 $\hat{\bf u}_{fh}^{m+1}\approx 2{\bf u}_{fh}^{m}-{\bf u}_{fh}^{m-1}$, $\hat{\phi}_{ph}^{m+1}\approx 2\phi_{ph}^{m}-\phi_{ph}^{m-1}$. 类似地, 无论压力 $\hat{p}_{fh}^{m+1}$ 是否经 time filter 更新, 最终的结果几乎一致.

下面给出解耦的线性多步法加 time filter 算法的稳定性分析.

${\bf定理3.2}$ (解耦算法的稳定性) 设 ${\bf u}_{fh}^{m+1}$, $\phi_{ph}^{m+1}$ 是线性多步法加 time filter 算法的解. 对于 $N\geq 2$, 有

其中 $C(T)=\exp\Big(\sum\limits_{m=1}^{N-1}\max\left\{\frac{8(1-\theta)\Delta t}{(3-2\theta)\hat{C}_{coe}},\frac{8(1-\theta)\Delta t}{(3-2\theta)\check{C}_{coe}} \right\}\Big)$.

${\bf证}$ 由 (3.15) 式可得

将上述三式分别代入 (3.13) 和 (3.14) 式后相加, $\forall\ {\bf v}_{fh}\in{\bf H}_{fh}, \psi_{ph}\in{\bf H}_{ph}, q_{fh}\in Q_{fh}$,有

假设 ${\bf v}_{fh}=2\Delta t B({\bf u}_{fh}^{m+1},\theta)$, $\psi_{ph}=2\Delta t B(\phi_{ph}^{m+1},\theta)$, $q_{fh}=2\Delta t B(p_{fh}^{m+1},\theta)$. 类似定理 3.1 的证明, 有以下不等式成立

这里主要分析交界面项. 借助引理 3.2, 并代入 $\varepsilon_1=\frac{1}{\hat{C}_{coe}}$, $\varepsilon_2=\frac{g}{\check{C}_{coe}}$, 有

其中, 最后一个不等式是在选择合适的常量 $C_1$ 和 $C_2$ 下得到的.

结合 (3.17)-(3.23) 式, 并将 (3.16) 式从 $m=1$ 直到 $m=N-1$ 求和, 整理可得

由于

以及

因此

结合离散的 Gronwall 引理可得解耦算法的稳定性结论. 证毕.

本节将证明耦合和解耦算法的误差估计. 记算法的精确解为 $(\underline{\bf u}^{m+1},p_f^{m+1})=$$(\underline{\bf u}(t_{m+1}),$$p_f(t_{m+1}))$, 其中 $\underline{\bf u}^{m+1}=({\bf u}_f^{m+1},\phi_f^{m+1})$, $\underline{\bf u}(t_{m+1})=({\bf u}_f(t_{m+1}),\phi_f(t_{m+1}))$, 同时定义误差函数

显然地, 下述关系式成立

需注意 $\eta_{\underline{\bf u}}^0=0$, $\eta_{p}^0=0$, $\eta_{{\bf u}_f}^0=0$, $\eta_{\phi_p}^0=0$.

假设所求解满足下述正则性条件

外力项 ${\bf g}_f$ 和 ${\bf g}_p$ 满足

${\bf定理4.1}$ (耦合算法的二阶收敛性) 在 (4.2)和(4.3) 式的假设下, 对于 $N\geq 2$ 有以下估计

${\bf证}$ 在 $t_{m+1}$, $t_{m}$ 和 $t_{m-1}$ 时刻分别对 (2.11) 式乘$\frac{(1-\theta)(3-2\theta)}{2}$, $(4\theta-1-2\theta^2)$, $\frac{(1-\theta)(1-2\theta)}{2}$ 得到三个对应时刻的等式. 从 (3.9) 式中减去三者之和, 结合误差函数的定义和投影算子的性质 (3.2) 式可得, $\forall\ \underline{\bf v}_h\in{\bf u}_h$, $\ q_{fh}\in Q_{fh}$,有

假设 $\underline{\bf v}_h=2\Delta t B(\eta_{\underline{\bf u}}^{m+1},\theta)$, $q_{fh}=2\Delta t B(\eta_{p}^{m+1},\theta)$. 类似定理 3.1 的证明, 对 (4.4) 式的每一项进行分析.

首先, 借助引理 3.1, 有

利用双线性形式 $a(\cdot,\cdot)$ 的强制性可得

由无散度特性可知, 左边第三项的结果为 0.

其次, 利用积分型余项的泰勒展开式,有

代入并整理得

再结合 Cauchy-Schwarz 不等式,有

可以得到

类似地, 借助积分型余项的泰勒展开式,有

代入并整理得

从而

类似地, 根据积分型余项的泰勒展开式,有

代入并整理得

从而

最后, 结合 (4.5)-(4.8) 式, 将 (4.4) 式从 $m=1$ 直到 $m=N-1$ 求和, 类似定理3.1的证明有

由三角不等式可得耦合算法误差估计的结论. 证毕.

${\bf定理4.2}$ (解耦算法的二阶收敛性) 在 (4.2) 和 (4.3) 式的假设下, 对于 $N\geq 2$ 有以下估计

${\bf证}$ 在 $t_{m+1}$, $t_{m}$ 和 $t_{m-1}$ 时刻分别对 (2.11) 式乘 $\frac{(1-\theta)(3-2\theta)}{2}$, $(4\theta-1-2\theta^2)$, $\frac{(1-\theta)(1-2\theta)}{2}$ 得到三个对应时刻的等式. 从 (3.16) 式中减去三者之和, 再结合误差函数的定义和投影算子的性质 (3.2) 式可得 $\forall$${\bf v}_{fh}\in {\bf H}_{fh}$, $\psi_{ph}\in {\bf H}_{ph}$, $q_{fh}\in Q_{fh}$,有

假设 ${\bf v}_{fh}=2\Delta tB(\eta_{{\bf u}_f}^{m+1},\theta)$, $\psi_{ph}=2\Delta tB(\eta_{{\phi}_p}^{m+1},\theta)$, $q_{fh}=2\Delta tB(\eta_{{p}_f}^{m+1},\theta)$. 回顾定理 4.1 的证明, 有以下不等式成立

下面主要分析交界面项. 交界面项可以被分解为

由线性性质可知, 上式中前四项之和等于 $0$. 下面四项借助 (2.12) 式及引理 3.2, 并代入 $\varepsilon_3=\varepsilon_5=\frac{3}{\hat{C}_{coe}}$, $\varepsilon_4=\varepsilon_6=\frac{3g}{\check{C}_{coe}}$, 有

其中, 最后一个不等式是在选择合适的常量 $C_3,C_4,C_5$ 和 $C_6$ 下得到的.

利用带有积分型余项的泰勒展开式,有

代入并整理得

类似地, 根据带有积分型余项的泰勒展开式,有

可以得到

利用引理 3.2 并代入 $\varepsilon_7=\frac{3}{\hat{C}_{coe}},\varepsilon_8=\frac{3g}{\check{C}_{coe}}$ 可得

其中, 第二个不等式是在选择合适的常量 $C_7$ 和 $C_8$ 下得到的.

综上, 交界面项可化简为

结合 (4.10)-(4.23) 式, 将 (4.9) 式从 $m=1$ 直到 $m=N-1$ 求和, 类似定理 3.2 的证明, 有

由三角不等式可得解耦算法误差估计的结论.证毕.

这一节将进行3个数值实验. 第1个实验对耦合和解耦算法的有效性进行模拟. 第2个实验说明算法的收敛阶由一阶提高到二阶, 并通过对比体现出解耦算法的高效性. 第3个实验将耦合和解耦的向后欧拉方法和取不同 $\theta$ 值的线性多步法, 以及在此基础上加 time filter 算法进行比较, 说明了线性多步法和线性多步法加 time filter 算法的误差相对较小. 本文的数值模拟是利用 FreeFEM++ 软件实现的. 并不失一般性, 以下实验均假设相应的物理参数 $\nu=1$, $\rho=1$, $g=1$, $S=1$$\alpha=1$ 以及${\bf K}=k{\bf I}$, 其中 $k=1$. 并且所有的初始条件, 边界条件和源项遵循精确解. 为了简便, 在图例, 图注及表头中, CBEM 和 DBEM 分别表示耦合和解耦的向后欧拉方法, CLMM 和 DLMM 分别表示耦合和解耦的线性多步法, CBEMTF 和 DBEMTF 分别表示耦合和解耦的向后欧拉方法加 time filter 算法, CLMMTF 和 DLMMTF 分别表示耦合和解耦的线性多步法加 time filter 算法.

${\bf 实验 1}$ 这里利用文献[37] 中的算例，设全局区域 $\Omega$ 由不可压缩流区域 $\Omega_f=(0,\pi)\times(0,1)$, 多孔介质区域 $\Omega_p=(0,\pi)\times(-1,0)$, 以及交界面 $\Gamma=(0,\pi)\times{0}$ 组成. Taylor-Hood 元 (P2-P1) 和分段二次多项式 (P2) 被用于自由流体流动方程和多孔介质流动方程. 精确解由下式给出

为验证算法的有效性, 本文设置网格尺寸为 $h=1/100$, 时间步长为 0.01 进行数值模拟. 图 2-图 4 分别展示了 $\theta=1/6$, $1/4$, $1/3$ 耦合和解耦的线性多步法以及线性多步法加 time filter 算法的速度等值线和速度流线, 通过观察可以发现本文提出的算法能够有效地模拟流体的流动.

${\bf 实验 2}$ 这里借助文献[31]中的算例, 考虑不可压缩流区域 $\Omega_f=(0,1)\times(1,2)$, 多孔介质区域 $\Omega_p=(0,1)\times(0,1)$, 以及交界面 $\Gamma=(0,1)\times{1}$, 并对自由流体流动方程和多孔介质流动方程分别使用 MINI 元 (P1b-P1) 和线性拉格朗日元 (P1). 下面给出精确解

为体现耦合和解耦算法的收敛阶, 本实验考虑 $\theta=1/4$, 并在表1-表8 中分别列出变化时间步长 $\Delta t$ 和网格尺寸 $h$ 对应的误差, 收敛阶及CPU 时间.

首先, 利用与文献[31] 中相同的方法, 通过固定网格尺寸 $h$, 变化时间步长 $\Delta t$ 计算耦合和解耦算法的收敛阶. 这时近似误差主要由时间步长 $\Delta t$ 控制. 定义

其中 $e_{v}(\Delta t)=\|v_h^{\Delta t}-v_h^{\frac{1}{2}\Delta t}\|_{L^2}$, $v={\bf u}_f$, $\phi_p$, $p_f$, $\rho_{v}\approx \frac{4^{\gamma}-2^{\gamma}}{2^{\gamma}-1}$. 固定网格尺寸 $h=\frac{1}{8}$, 时间步长 $\Delta t$ 由 $\frac{1}{20}$ 变化到 $\frac{1}{320}$ 得到表 1 -表 4. 将表1 和表2, 表 3 和表4分别进行比较, 发现耦合和解耦的线性多步法是一阶收敛的, 而本文提出的线性多步法加 time filter 算法在几乎不增加计算量的情况下可以达到二阶收敛. 这体现了 time filter 算法的二阶收敛性和高效性. 同时通过比较表中 CPU 时间可以发现解耦算法的计算效率相对较高.

其次, 在固定时间步长 $\Delta t$, 变化网格尺寸 $h$ 的情况下, 计算耦合和解耦算法的收敛阶. 这时近似误差主要由时间步长 $h$ 控制. 定义

其中 $e_{v}(h)=\|v_{f}-v_{fh}\|_{L^2}$ 是关于网格尺寸 $h$ 的误差. 固定时间步长 $\Delta t=0.01$, 网格尺寸 $h$ 由 $\frac{1}{4}$ 变化到 $\frac{1}{64}$ 得到表 5 -表 8. 将表 5 和 表 7, 表6 和 表 8 中的 CPU 时间分别进行比较, 可以看出解耦算法所用的 CPU 时间较少, 计算效率相较于耦合算法较高.

${\bf 实验 3}$ 这里, 利用文献 [38-39] 中的算例进行数值实验. 同时选择与实验 2 中相同的区域 $\Omega_f$ 和 $\Omega_p$, 以及有限元. 该精确解是

图 5 分别给出了向后欧拉方法和 $\theta$ 取不同值的线性多步法在时间段为 $[T]$ ($T=5$) 且时间步长为$0.05$ 上 ${\bf u}_f$ 和 $\phi_p$ 在每个时刻的绝对误差. 从图 5 中可以看出: 与向后欧拉方法相比较, 线性多步法关于 ${\bf u}_f$ 和 $\phi_p$ 的绝对误差几乎在每个时刻都相对较小. 同时, 当 $\theta$ 取不同值时线性多步法的绝对误差会随着 $\theta$ 的增大而减小. 类似地, 图 6 分别给出了向后欧拉方法加 time filter 算法和 $\theta$ 取不同值的线性多步法加 time filter 算法在时间段为 $[T]$ ($T=5$) 且时间步长为 $0.05$ 上 ${\bf u}_f$ 和 $\phi_p$ 在每个时刻的绝对误差. 由图 6 知: 与向后欧拉方法加 time filter 算法相比, 无论是耦合还是解耦线性多步法加 time filter 算法关于 ${\bf u}_f$ 和 $\phi_p$ 的绝对误差几乎在每个时刻都相对较小, 同时随着 $\theta$ 的增大, 线性多步法加 time filter 的绝对误差减小.