令${\Omega}$是${{\Bbb R}} ^2$上的一个有界区域,它是由两个子区域${\Omega}_{1}$和${\Omega}_{2}$通过一个共享界面$I$耦合而成,并在每个界面的边界$\partial {\Omega}_{i}$上有其外法线向量${\bf n}_{i} (i = 1, 2)$.接着,考虑该区域的流体-流体相互作用模型,该模型已进行了初步研究^[13],如下所示:给定粘性系数$\nu_{i} > 0,$摩擦系数$\kappa\in {{\Bbb R}},$外力项${\bf f}_{i}:[0, T]\rightarrow H^{1}(\Omega_{i})^2$和初始值${\bf u}_{i}(0)\in H^{1}(\Omega_{i})^2$,找到${\bf u}_{i}:[0, T]\times \Omega_{i}\rightarrow{{\Bbb R}} ^{2}$和$p_{i}:[0, T]\times \Omega_{i}\rightarrow{{\Bbb R}}$满足

其中$|\cdot|$代表欧几里得范数, ${\bf u}_{i}$和$p_{i}$分别代表了速度和压力.除此之外, $\tau$是任意切向量且满足$\tau \cdot {\bf n}_{i} = 0.$值得注意的是该流体-流体相互作用模型应用了界面$I$上的非线性界面条件,该界面$I$用于对线段进行建模.

众所周知,在许多重要的科学,工程和工业应用中需要一种流体与另一种流体的多领域和多物理耦合问题的近似解(参见文献[3-4, 10, 32]).事实上,流体-流体相互作用模型就是其中之一,它用来模拟两种不同粘性流体流动的动态特征,如不均匀的血液流动^[14]、海洋-大气界面^{[3-4, 23-25]}等等.虽然该模型可以将海洋-大气问题的动态核心简化为最简单的形式,但其仍然保留了模拟海洋-大气问题中的一些难点.此外,在数值方法方面,它仍然是一类重要且具有挑战性的问题.主要难点在于它是压力项、不可压缩条件、非线性项和具有某些非线性界面项相互作用条件构成的复杂耦合系统.因此,人们在研究该模型有效的数值方法上投入了大量的精力.

Bresch和Koko^[7]采用算子分裂法和基于优化的非重叠区域分解方法,分别求解了一个耦合退化的Stokes问题和一个耦合线性的Advection-diffusion问题.此外, Connors等人^[13]提出了两种解耦的时间步径方法,即隐显式法(implicit-explicit method)和几何平均法(geometric-averaging method),并证明了几何平均法是无条件稳定的.相比,隐显式方法虽然是最简单,最自然的解耦方法^[12],但在一定条件下是稳定的,这已被Zhang等人^[31]证明.由Connors和Howell^[11]提出了另一种无条件稳定方法,该方法采用了解耦子问题和不同时间步长策略求解流体-流体相互作用模型.最近, Qian等人^[26]针对流体-流体相互作用问题提出了一种局部投影稳定和特征解耦的时间步径方法.其采用几何平均法处理非线性界面条件,得到了一种无条件稳定的分块方法.这种处理非线性界面的方法也被用于其他的一些算法^[21-22].除此之外, Aggul等人^[1]提出了一种预测-校正类型方法,该方法也是无条件稳定的,且具有二阶时间精度. Conors^[9]提出了一个统计湍流模型,用于在平面界面上耦合两种流体的集成计算.

众所周知,关于两步有限元法的文献非常丰富. Chorin^[8]和Temam^[30]设计了一种投影算法,该算法经历了一些演变并得到了进一步的发展,如压力校正方法^[28]和矩阵分解方法^[27]. Chorin和Temam的算法是基于一个中间速度场在空间上的投影向量场,将求解鞍点系统线性解的困难转化为两个较简单的线性解,并可以通过在粘性步中显化压力和在投影步中修正压力来改进.在投影步骤中加入一个压力修正项,根据精确解的一些正则性假设,改进了Chorin/Temam投影法对中间速度和压力的估计.然而,投影法的步末速度不满足精确的边界条件.因此,根据Chorin/Temam投影算法的思想,提出了分步法^{[16, 29]},其中在每一步中只需求解一个泊松方程得到压力解,大大降低了计算成本.

与标准的投影方法不同,粘性分离法是一种分步法,它对粘分项进行了分割,允许对所有子步骤的速度执行完整的狄利克雷边界条件.此外,该方法是对对流效应与不可压缩性的解耦.因此,近年来粘性分离法得到了广泛的关注. Blasco等^[6]提出了求解三维的非定常不可压缩Navier-Stokes方程的一阶时间离散粘性分离方法.在连续解的正则性假设下,得到了时间离散格式的一些误差估计^[5].此外, Guillén-González等^[17]得到了时间离散方案的最优误差估计,其中必须包含初始时间步长的权重,以推导出压力的最优误差估计.再者,他们还利用时间离散格式作为辅助问题来研究完全离散的有限元格式,获得速度和压力的最佳一阶近似(参见文献[18]).此外,粘性分离算法还被应用于海洋的原始方程^[19], Oldroyd模型^[35]和Boussinesq问题^[34]. Zhang等^[33]针对粘弹性流问题,提出了一种大时间步径的粘性分离方法,并在稳定条件下得到了时间半离散格式的误差估计.上述基于粘性分离技巧的算法均采用空间连续有限元法.对于不连续有限元法,在Blasco等^[6]介绍的粘性分离算法的基础上, Girault等人^[15]采用对称内罚Galerkin法和非对称内罚Galerkin法求解了非定常的不可压缩Navier-Stokes方程.

本文建立并研究了全离散粘性分离有限元法求解流体-流体相互作用模型.与许多模拟流体-流体相互作用模型方法一样,该方法也是无条件稳定的.此外,该方法还采用了粘性在时间上的分解技术和空间上的混合有限元法,其中时间项包括两个步骤.该方法允许对所有子步骤的速度施加完整的狄利克雷边界条件,并将对流影响与不可压缩性解耦.论文的其余部分安排如下:在第2章,引入了一些数学基础,并给出模型(1.1)的变分形式.在第3章中,列出流体-流体相互作用模型的全离散粘性分离有限元方法.第4章列出了该方法的稳定性和误差估计.最后通过数值实验对全离散化方案的理论结果进行了验证.

在本章,我们回顾一些必要的符号和引理,稍后将使用它们.分别用$\|\cdot\|_{0}$和$(\cdot, \cdot)$来表示常见的$L^{2}(\Omega)$范数及其内积,用$L^{p}(\Omega)$范数和$W^{m, p}(\Omega)$范数代表$\|\cdot\|_{L^{p}(\Omega)}$和$\|\cdot\|_{W^{m, p}(\Omega)}$对于$m\in\mathbb{N^{+}}$, $1\leq p\leq \infty$.特别的, $H^{m}(\Omega)$被用来表示Sobolev空间$W^{m, 2}(\Omega)$. $\|\cdot\|_{m}$用于表示在$H^{m}(\Omega)$的范数. $X$是$\Omega$中的赋范函数空间, $L^{p}(0, T; X)$表示定义在$[0, T]\times \Omega$上所有函数的空间,在此函数空间的范数表示为

且其范数是有界的.当$p = \infty$,表示在时间上取最大的范数.

对于流体-流体相互作用模型(1.1),我们引入以下Hilbert空间:对于$i = 1, 2$,有

和

${\bf f}_{i}$是定义在$X_{i}$对偶空间$X_{i}'$的函数,其范数表示为

基于上述函数空间的定义,问题(1.1)的变分形式如下:找到${\bf u}_{i}:[0, T]\rightarrow X_{i}$和$p_{i}:[0, T]\rightarrow M_{i}$ ($i, j = 1, 2, i\neq j$),使得对于任意的$({\bf v}_{i}, q_{i})\in (X_{i}, M_{i})$,满足

其中$({\bf u}_{i, t}, {\bf v}_{i}) = \int_{\Omega_{i}}\frac{\partial {\bf u}_{i}}{\partial t}{\bf v}_{i}{\rm d}{\Omega}_i, i = 1, 2,$连续的双线性项$a(\cdot, \cdot)$和$d(\cdot, \cdot)$分别被定义在$X_{i}\times X_{i}$和$X_{i} \times M_{i}$,表示为

对于函数${\bf u}_{i},$${\bf v}_{i},$${\bf w}_{i}\in X_i$,非线性项$b(\cdot, \cdot, \cdot)$定义为

非线性项的上界和一些性质将在下一个数值分析中用到,在下面的引理中给出.

引理2.1^[13] 给定${\bf u}_{i},$${\bf v}_{i},$${\bf w}_{i}\in X_{i},$对于$i = 1, 2$,有

此外,如果${\bf v}_{i}\in H^{2}(\Omega_{i})^2,$则

其中$C > 0$是一个依赖于$\Omega_{i}$的正常数.

此后的$C$将被记为一个正常数,它最多取决于$\kappa$, $\Omega_{i}$, $T$, $\nu_{i}$, ${\bf u}_{i}$和$p_{i}$,其在不同的情况下代表不同的值.在下面的数值分析中,将使用以下界面上的有界性.

引理2.2^[13] 假设$\alpha_{i},$$\alpha_{j},$$\beta_{i},$$\beta_{j}$是任意正常数,则对于任意的$({\bf u}_{i}, {\bf v}_{j})\in (X_{i}, X_{j}),$和$({\bf w}_{i}, {\bf w}_{j})\in (X_{i}, X_{j})$,满足

其中$\|\cdot\|_{I}$表示$L^{3}(I)$,且其可以被$H^1$ -范数约束.

众所周知,离散Gronwall引理在收敛性分析中将扮演重要的角色,因此我们将其列在下面的引理中.

引理2.3^[20] 假设$k$和$a_{n},$$b_{n},$$d_{n},$及整数$n_{1}\leq n$都是非负数且满足

当$s = m-1$时,有

首先,对于每个正整数$N$,令$\{t_n = n\Delta{t}\}_{n = 0}^{N}$是对$[0, T]$均匀划分得到的时间步长,其中时间步长$\Delta{t} = \frac TN$.其次,我们定义$X_{i}^{h}\subset X_{i}$和$M_{i}^{h}\subset M_{i} (i = 1, 2)$是速度和压力的有限元空间,子区域$\Omega_{i}$的三角形剖分定义为$\pi^{h}_{i} = \left\{K_{i}\right\}$,其中$K_{i}$是剖分后的三角形,且其直径定义为$h_{i, K_{i}}$,记${ } h_{i} = \max_{K_{i}\in\pi_{i}^{h}}h_{i, K_{i}}$是$\pi^{h}_{i}$中最大的直径.此外,有限元空间$X_{i}^{h}\times M_{i}^{h}$假设满足inf-sup条件或者为了离散压强的稳定性满足${\rm LBB}^{h}$条件

其中$\beta_{i}^{h}$是一个不依赖于$h_{i}$的正常数.更进一步,定义$h = \max\{h_1, h_2\}$是$\Omega$三角形剖分中最大的直径.此外,我们介绍$X_i^{h}$的子空间$V_h$,定义为

本文在离散空间中采用${\rm P_{1}+bubble/P_{1}}$ (MINI-element^[2])有限元对.最后,记$({\bf u}_{i}^{n}, p_{i}^{n})$是模型(1.1)在时间$t_{n}$的有限元近似解.

全离散的粘性分离有限元方法包括两个步骤.第一步:给定${\bf u}_{i}^{n-1},$${\bf u}_{i}^{n}\in X_{i}^{h}$和${\bf u}_{j}^{n-\frac12}\in X_{j}^{h}$,对$n = 0, 1, \cdots, N-1$,找到${\bf u}_{i}^{n+\frac{1}{2}}\in X_{i}^{h}$使得

第二步:基于从(3.1)得到的${\bf u}^{n+\frac{1}{2}}_{i}$,找到${\bf u}_{i}^{n+1}\in X_{i}^{h}$和$p^{n+1}_i\in M_{i}^{h}$使得

在上面全离散的有限元方法(3.1)和(3.2)中,采用半隐格式处理非线性项$b(\cdot, \cdot, \cdot)$会得到一个具有变系数矩阵的线性系统,这样可以节省大量的计算量.而且对于非线性流体界面项,我们采用几何平均的方法(geometric averaging method)^[13].与耦合方法相比,使用粘性分离方法(3.1)和(3.2)的优点是对流效应与不可压缩性的解耦.同时,我们在$\Gamma_{i}$上执行相同的齐次狄利克雷边界条件.

在本章中,我们将研究上一章中关于全离散粘性分离有限元方法(3.1)和(3.2)的稳定性和误差估计.

定理4.1 假设${\bf u}_{i}^{m+1}$$(i = 1, 2)$是算法(3.1)和(3.2)的解,则对于$m = 1, 2, \cdots, N-1,$满足

证 一方面,令(3.1)式中${\bf v}_{i} = 2\Delta t{\bf u}_{i}^{n+\frac{1}{2}}$得到

对于上面等式中的界面项,可得

其次,对于(4.1)式的右端项,有

将(4.1)式, (4.2)式和(4.3)式结合得

另一方面,令(3.2)式中${\bf v}_{i} = 2\Delta t{\bf u}_{i}^{n+1}$可得

根据(4.4)式和(4.5)式得

最后,对(4.6)式从$n = 1, 2, \cdots, m$和$i = 1, 2$求和,即完成证明.

在本章接下来的部分中,我们将主要讨论流体-流体相互作用模型(1.1)的全离散粘性分离方法的误差估计.为了证明粘性分离方法的收敛性,我们回忆Stokes投影^[31].

找到$(R_{i}{\bf u}_{i}, T_{i}p_{i})\in(X_{i}^{h}, M_{i}^{h}) (i = 1, 2)$满足

并且满足下列最优逼近性质

其次,将误差分解成如下格式

其中

现在,根据上述定义,我们得到算法(3.1)和(3.2)的误差估计.

定理4.2 若$({\bf u}_{i}(t_{n+1}),$$p_{i}(t_{n+1}))$和$({\bf u}_{i}^{n+1}, p_{i}^{n+1})$分别是(1.1)式和全离散粘性分离算法(3.1)和(3.2)的解,连续解满足正则性假设: ${\bf u}_{i}\in L^{2}(0, T; H^{2}(\Omega_{i})^2),$${\bf u}_{i, t}\in L^2(0, T; X_i),$$p_{i}\in L^{2}(0, T; H^1({\Omega}_i))$, ${\bf u}_{i, tt}\in L^{2}(0, T; L^{2}(\Omega_{i})^2)$,初始速度值${\bf u}_{i}^{0},$${\bf u}_{i}^{1}$和${\bf u}_{i}^{\frac12}$满足以下估计

并且时间步长$\Delta t$满足条件$\Delta t\leq \Delta t_0$,其中$\Delta t_0$是一个正常数,则有以下误差估计成立

证 一方面,在$t_{n+1}$时, (2.1)式减去(3.1)式得

(4.10)式中的第一项可以重新写作

将(4.9)式中的误差形式分开,并将(4.11)式带入到(4.10)式得

另一方面, (3.2)式加减$R_{i}{\bf u}_{i}(t_{n+1})$得

此外,令(4.12)式中${\bf v}_{i} = 2\Delta t{\bf \tilde{e}}_{i}^{n+\frac{1}{2}}\in V_{i}^{h}$,可得

通过加减$b({\bf u}_{i}(t_{n}), {\bf u}_{i}(t_{n+1}), {\bf \tilde{e}}_{i}^{n+\frac{1}{2}}),$三线性项可以重新写作

现在,开始估计上式的每一个右端项.首先,从三线性项的反对称性质,不难看出

然后,使用引理2.1和Young不等式对(4.15)式右端的每一项进行估计,如下

其次,显然有

同样,可以得到

此外,为了简化证明,我们引入了一些符号.令$[{\bf u}^{n+1}] = {\bf u}_{i}^{n+1}-{\bf u}_{j}^{n+1}$, $[{\bf u}(t_{n+1})] = {\bf u}_{i}(t_{n+1})-{\bf u}_{j}(t_{n+1})$, $[R{\bf u}(t_{n+1})] = R_{i}{\bf u}_{i}(t_{n+1})-R_{j}{\bf u}_{j}(t_{n+1}).$最后, (4.14)式中的界面项可以写作

和

记(4.18)式中的右端前四项为$A_1, \cdots, A_4$,记(4.19)式中的右端每一项为$B_1, \cdots, B_7$,则有如下估计成立.为了估计这些项,首先观察到

其中使用了反三角形不等式.根据以上估计,我们有

和

然后,应用(2.2)式, (4.20)式, (4.24)式和Young不等式,得

此外, $A_{3}$和$B_{2}$可以被(2.4)式, (4.21)式和(4.23)式约束, $B_5$可以被(2.3)式约束,如下

同样, $B_{1}$可以被(4.22)式, (2.4)式和(2.2)式约束

对于$B_{4}$,有

这里使用了与$A_1$相同的技巧.

最后,将(4.25)式, (4.26)式, (4.27)式和(4.28)式结合,得

和

这里选取任意数$\varepsilon_{24} = \varepsilon_{9}+\varepsilon_{10}+\varepsilon_{11}+\varepsilon_{15}, \varepsilon_{25} = \varepsilon_{12}+\varepsilon_{13}+\varepsilon_{14}+\varepsilon_{17} +\varepsilon_{19}+\varepsilon_{21}+\varepsilon_{23}$,和$\varepsilon_{26} = \varepsilon_{18}+\varepsilon_{22}$.

令$\varepsilon_{27} = \varepsilon_{24}+\varepsilon_{25},$$\varepsilon_{28} = \max\{\varepsilon_{16}, \varepsilon_{16}+\varepsilon_{26}, \varepsilon_{20}\}$,将(4.29)式与(4.30)式相加,则界面项被约束

其中, $M^{n+1}$, $L^{n+1}$和$P^{n+1}$定义为

现在,对(4.14)式右端的每一项进行估计,应用Cauchy-Schwarz和Young不等式可得

(4.14)式中的压力项可以被约束

其中$d({\bf \tilde{e}}_{i}^{n+\frac{1}{2}}, T _{i}p_{i}(t_{n+1})) = 0$.

将(4.33)式, (4.32)式, (4.31)式, (4.17)式, (4.16)式合并,选取$\varepsilon_{1}+\varepsilon_{2}+\varepsilon_{4}+\varepsilon_{5}+\varepsilon_{7} +\varepsilon_{8}+\varepsilon_{27}+\frac12\varepsilon_{29}+\varepsilon_{30}+\varepsilon_{31} +\varepsilon_{32} = \frac{1}{8},$$\varepsilon_{3}+\varepsilon_{6} = \frac{1}{8},$$\varepsilon_{28} = \frac{1}{64},$ (4.14)式重新写作

现在,我们对$\Delta tL^{n+1}P^{n+1}$进行约束.根据解的正则性假设和定理4.1,我们得到

其中,我们使用了Poincaré不等式和$\|{\bf v}_i\|_{I}\leq C\|{\bf v}_i\|_{0}^{\frac{1}{4}}\|\nabla {\bf v}_i\|_{0}^{\frac{3}{4}},$${\bf v}_{i}\in H^1({\Omega}_i)$ (引理2.1^[13]).更进一步,根据投影(4.8)的性质,得$\|\nabla {\bf \hat{e}}_{i}\|_{0}^{2}\leq Ch^{2}.$因此, $P^{n+1}$可以被约束

从(4.36)式和(4.35)式可得

最后,对$i = 1, 2$和$n = 1, 2, \cdots, m$求和,注意到

重写(4.34)式如下

其中, $D^{n+1}$定义为

接下来,令(4.13)式中${\bf v}_{i} = 2\Delta t{\bf \tilde{e}}_{i}^{n+1}$和$q_i = -2\Delta tp_i^{n+1}$,对$n = 1, 2, \cdots, m$和$i = 1, 2$求和,得到

将(4.37)式和(4.38)式合并得到

假设$\Delta t\leq \Delta t_0$,使得$C\Delta t D^{n+1}\leq1$成立,并使用Gronwall引理和(4.8)式,我们完成了证明.

在本章中,通过一些数值试验来验证上一章中粘性分离有限元算法的理论结果.假设$\Omega_{1} = [0, 1]\times[0, 1],$$\Omega_{2} = [0, 1]\times[-1, 0]$和$I = (0, 1)\times \{0\}$.问题(1.1)的解析解^[26]给出如下

右端项${\bf f}_1 = (f_{1, 1}(t, x, y), f_{1, 2}(t, x, y))$和${\bf f}_2 = (f_{2, 1}(t, x, y), f_{2, 2}(t, x, y))$的选取必须满足: $({\bf u}_1, p_1)$和$({\bf u}_2, p_2)$是式(1.1)的解.

在前两个数值测试中,我们将参数值设置为$\nu_1 = 0.005$, $\nu_2 = 0.005$和$\kappa = 100$^[26].误差表示为

为了证明本文中粘性分离算法方法的无条件稳定性,我们计算了不同空间网格、不同时间步长下的$\|\nabla {\bf u}_{1}^{n}\|_{0}$和$\|\nabla {\bf u}_{2}^{n}\|_{0}$的值,最后的迭代时间$T = 1$,结果如表 1和表 2所示.

通过比较表 1和表 2中$\|\nabla {\bf u}_{1}^{n}\| _{0}$和$\|\nabla {\bf u}_ {2}^{n}\|_{0}$在同一空间网格不同时间步长下的值,很容易得到,他们趋向于一些常数,这表明在两个域$\Omega_{i} (i = 1, 2)$不需要时间步长限制.

一方面,我们测试了该方法相对于$h$的收敛速度.为了避免时间步长的影响, $\Delta t$应该选择足够小的值.因此,我们选择$\Delta t = 0.0001$和0.001,其中$T = 0.1$,依次取$h$ = 1/8、1/16、1/32和1/64.

我们在表 3和表 4中分别给出$\Delta t = 0.0001$和0.001的收敛速度.由表 3和表 4可知,收敛阶大于1.

另一方面,我们测试该算法相对于$\Delta t$的收敛速度.在这里,令$t = h$.此外,我们依次取$\Delta t$ = 1/4、1/8、1/16、1/32、1/64.

表 5列出了本文算法求解$T = 1$时的数值结果.根据表 5,得到速度相对于$\Delta t$的收敛阶大于1.

为了验证所提出的方法对非连续运动粘性问题的有效性,我们选择了粘性系数值为$\nu_1 = 0.005$和$\nu_2 = 0.01$.

表 6列出了用粘性分离有限元算法求解的固定时间步长$\Delta t = 0.01$的数值结果.此外, 表 7列举出了$\Delta t = h$时的一些结果.从表 6和7中可以看出,在$\|\cdot\|_{l^2(0, T; H^1(\Omega_i)^2)}$范数意义下,对于不连续运动粘性问题,粘性分离有限元算法的收敛速度均可以达到一阶精度(或更高).