基于亏量校正求解自然对流问题的后验误差估计子

Recovery-Based Error Estimator for the Natural-Convection Problem Based on Defect-Correction Method

摘要:

Abstract: 本文给出NC方程基于亏量校正方法的恢复型误差估计子。我们在亏量步求解人工黏性系数的稳定非线性问题；在校正步通过线性化方程来校正这个残量，在亏量和校正步我们都采用Oceen迭代。考虑应力张量 ，若我们采用最低阶协调速度和压力有限元对 求解离散NC问题，应力张量的有限元逼近 是不连续的分片常数。我们利用超收敛片恢复技巧构造连续空间的量 ，得到了局部恢复型误差估计子 。最后，通过数值实验验证了基于亏量校正方法的恢复型误差估计子是有效的。 This article solves an artificial viscosity stabilized nonlinear problem in the defect step, and corrects the residual by linearized equations in the correction step for a few steps. In both the defect and correction steps, we use the Oseen iterative scheme to solve the discrete nonlinear equations, considering the stress tensor . For the speed and pressure of NC problem, we use the lowest finite element pair . The discrete finite element stress tensor approximation is discontinuous piecewise constant. We obtain the local recovery type error estimator by the convergence recovery technique to construct in continuous space. Finally, the stability, accuracy and efficiency of the proposed method are confirmed by several numerical investigations.

1. 引言

Natural convection (NC)自然对流问题，即浮力驱动流体运动，在实际问题中经常出现。定常NC问题是粘性流体不可压缩流动和传热过程的耦合方程，其中不可压缩流体是Boussinesq’s近似，在大气动力学中，它是一个重要的强迫耗散非线性系统，它不仅包含速度矢量场和压力场，而且包含温度场，从热力学角度讲，我们知道，流体的运动粘度会产生一定的热量，因此，流体的运动必须伴随着温度的相互转换，速度和压力。因此，对这一非线性系统的研究具有重要意义。对于自然对流方程的数值分析和数值结果，Christie和Mitchell [1] ，Boland和Layton [2] 等人进行了很多研究。此外，Ҫɪbɪk和Kaya [3] 提出了基于投影的稳定化有限元方法求解稳态自然对流问题。

我们提出一种基于亏量校正方法后验恢复的误差估计方法来解决NC问题。亏量校正方法的主要思想是：我们在亏量步求解人工黏性系数的稳定非线性问题；在校正步通过线性化方程来校正这个残量，在亏量和校正步我们都采用Oceen迭代。众所周知，缺陷校正方法是求解非线性稳态问题的一种迭代改进技术。该方法在不需要网格的情况下，提高了解的精度，对对流占优问题有效。它被用来解决Layton等人的Navier-Stokes方程 [4] 。作为一个非线性问题，NC方程的近似误差由两部分组成：一个是由有限元离散化产生的，另一个是由非线性方程迭代产生的。由Oceen迭代方案缺陷校正技术相结合，提出了一种新的算法。在亏量及校正的步骤中，我们利用Oceen方案求解离散方程，这种新方法继承了每个算法的最佳特征，它是无条件稳定的。

本文中，我们采用最低阶协调速度和压力有限元对 $P_1-P_0-P_1$ 离散NC问题。应力张量 $\sigma =\nabla u-pI+\nabla T$ 的有限元逼近 ${\sigma }_h=\nabla u_h-p_{h}I+\nabla T_h$ 是不连续的分片常数。我们利用超收敛片恢复技巧 [5] 构造连续空间的量 $G\left({\sigma }_h\right)$ ，得到了局部恢复型误差估计子 ${\Vert {\sigma }_h-G\left({\sigma }_h\right)\Vert }_K$。这样计算出的近似解是比一般有限元方法求解的近似解更加逼近真解。所以，我们结合亏量校正方法和恢复型误差估计子求解稳定的NC方程，不仅能有效地解决对流占有问题，还能得到更加精确的近似解。

本文主要分为4个部分，首先介绍NC方程的主要性质和经典结果。第二部分提出亏量校正方法，并对其稳定性和误差估计进行简单分析。第三部分，结合恢复型误差估计子求解NC方程。最后数值实验验证了我们方法的有效性和可靠性。

2. 预备知识

2.1. 定常NC方程

假设Ω是一个有界、凸开区域，其边界 $\partial \Omega$ 是Lipschitz连续的，令 ${\Gamma }_T=\partial \Omega /{\Gamma }_B$ ，其中 ${\Gamma }_B$ 是 $\partial \Omega$ 的正则开子集。我们考虑下面的NC方程：

$\{\begin{array}{l}-Pr\Delta u+\left(u\cdot \nabla \right)u+\nabla p=PrRajT,\text{in}\Omega ,\\ \nabla \cdot u=0,\text{in}\Omega ,\\ -\kappa \Delta T+u\cdot \nabla T=\gamma ,\text{in}\Omega ,\\ u=0,\text{in}\partial \Omega T=0,\text{in}{\Gamma }_T,\partial T/\partial n=0,\text{in}{\Gamma }_B,\end{array}$ (2.1)

这里 $u=\left(u_1\left(x\right),u_2\left(x\right)\right)$ 是流体速度， $p=p\left(x\right)$ 是流体压力， $T\left(x\right)$ 是温度， $\gamma$ 是外力方程， $Pr,Ra>0$ 是Prandtl数和Rayleigh数， $\kappa >0$ 是热传导系数， $j=\left(0,1\right)$ 是一个二维向量，n是 ${\Gamma }_B$ 的外法线。

接下来，介绍一些Hilbert空间：

$X=H_0^1{\left(\Omega \right)}^2,W=\left\{s\in H^1\left(\Omega \right):s=0\text{on}\Gamma \right\}$ ，

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

$L^2\left(\Omega \right)$ 中定义L²-标量内积 $\left(\cdot ,\cdot \right)$ 和L²-范数 ${\Vert \cdot \Vert }_0$。我们在空间X和W中同样定义内积 $\left(\nabla u,\nabla v\right)$ 和范数 ${\Vert \nabla u\Vert }_0$。标准定义是在Sobolev空间 $H^{m,p}\left(\Omega \right)$ 中，其范数为 ${\Vert \cdot \Vert }_{m,p}$ ， $m,p\ge 0$。我们将 $H^{m,2}\left(\Omega \right)$ 和其范数 ${\Vert \cdot \Vert }_m$ 写作和 ${\Vert \cdot \Vert }_{m,2}$。此外，定义 $H_0^1\left(\Omega \right)$ 的对偶空间为 $H^{-1}\left(\Omega \right)$ ，对偶空间中的范数定义为
${\Vert f\Vert }_{-1}=\underset{v\in H_0^1}{\sup }\frac{\left|\left(f,v\right)\right|}{{\Vert \nabla v\Vert }_0}$ ，且 $X^{-1}$ 是X的对偶空间。

方便起见，我们假设u，v属于同一个有限元空间X。下面分别定义两个连续双线性形式 $a\left(\cdot ,\cdot \right)$ 和 $d\left(\cdot ,\cdot \right)$ 和一个三线性形式 $b\left(u;v,w\right)$ ：

$a\left(u,v\right)=\left(\nabla u,\nabla v\right),\forall u,v\in X,$

$d\left(v,q\right)=\left(q,\nabla \cdot v\right),\forall v\in X,\forall q\in M,$ (2.2)

$b\left(u;v,w\right)=\left(\left(u\cdot \nabla \right)v,w\right)+\frac{1}{2}\left(\left(\nabla \cdot u\right)v,w\right)=\frac{1}{2}\left(\left(u\cdot \nabla \right)w,v\right)-\frac{1}{2}\left(\left(u\cdot \nabla \right)v,w\right),\forall u,v,w\in X.$ (2.3)

此外，我们需要在定义两个双线性形式 $\bar{a}\left(\cdot ,\cdot \right)$ 和一个三线性形式 $\bar{b}\left(\cdot ;\cdot ,\cdot \right)$ 在 $W\times W$ 和 $X\times W\times W$ 中，有

$\bar{a}\left(T,s\right)=\left(\nabla T,\nabla s\right),\forall T,s\in W,$

$\bar{b}\left(u;T,s\right)=\left(\left(u\cdot \nabla \right)T,s\right)+\frac{1}{2}\left(\left(\nabla \cdot u\right)T,s\right)$ $=\frac{1}{2}\left(\left(u\cdot \nabla \right)T,s\right)-\frac{1}{2}\left(\left(u\cdot \nabla \right)s,T\right),\forall u\in X,\forall T,s\in W.$

根据上面定义的双线性形式和三线性形式，有以下两个已知结论( [6] [7] [8] )：

$b\left(u;v,w\right)=-b\left(u;w,v\right),\left|b\left(u;v,w\right)\right|\le N{\Vert \nabla u\Vert }_0{\Vert \nabla v\Vert }_0{\Vert \nabla w\Vert }_0,\forall u,v,w\in X,$ (2.4)

$\bar{b}\left(u;T,s\right)=-\bar{b}\left(u;s,T\right),\left|\bar{b}\left(u;T,s\right)\right|\le \bar{N}{\Vert \nabla u\Vert }_0{\Vert \nabla T\Vert }_0{\Vert \nabla s\Vert }_0,\forall u,T,s\in X,$ (2.5)

其中

$N=\underset{u,v,w\in X}{\sup }\frac{\left|b\left(u;v,w\right)\right|}{{\Vert \nabla u\Vert }_0{\Vert \nabla v\Vert }_0{\Vert \nabla w\Vert }_0},\bar{N}=\underset{u\in X,T,s\in W}{\sup }\frac{\left|\bar{b}\left(u;T,s\right)\right|}{{\Vert \nabla u\Vert }_0{\Vert \nabla T\Vert }_0{\Vert \nabla s\Vert }_0},$

是两个只依赖于 $\Omega$ 的定常数。

现在给出一个定义在 $\left(X,M\right)\times \left(X,M\right)$ 上的广义双线性形式：

$B\left(\left(u,p\right);\left(v,q\right)\right)=Pra\left(u,v\right)-d\left(v,p\right)+d\left(u,q\right),\forall \left(u,p\right),\left(v,q\right)\in \left(X,M\right).$

它满足连续性和 $inf-sup$ 条件：

$\left|B\left(\left(u,p\right);\left(v,q\right)\right)\right|\le c\left({\Vert \nabla u\Vert }_0+{\Vert p\Vert }_0\right)\left({\Vert \nabla v\Vert }_0+{\Vert q\Vert }_0\right),\forall \left(u,p\right),\left(v,q\right)\in \left(X,M\right),$

$\underset{v,q\in \left(X,M\right)}{\sup }\frac{\left|B\left(\left(u,p\right);\left(v,q\right)\right)\right|}{{\Vert \nabla v\Vert }_0+{\Vert q\Vert }_0}\ge {\beta }_1\left({\Vert \nabla u\Vert }_0+{\Vert p\Vert }_0\right),\forall \left(u,p\right),\left(v,q\right)\in \left(X,M\right),$

其中， $c>0$ 是常数。它是与网格尺寸无关的，但是依赖于X和本文中的其他参数， ${\beta }_1>0$ 依赖于 $\Omega$ 和v。

NC方程(2.1)的变分形式如下：求解 $\left(u,p,T\right)\in X\times M\times W$ 对于任意的 $\left(v,q,s\right)\in X\times M\times W$ ，

$\{\begin{array}{l}B\left(\left(u,p\right);\left(v,q\right)\right)+b\left(u;u,v\right)=PrRa\left(jT,v\right),\\ \kappa \bar{a}\left(T,s\right)+\bar{b}\left(u;T,s\right)=\left(\gamma ,s\right).\end{array}$ (2.6)

下面给出原方程的解的存在唯一性的经典结论：

存在至少一个解对 $\left(u,p,T\right)\in X\times M\times W$ 满足(2.6)且

${\Vert \nabla T\Vert }_0\le {\kappa }^{-1}{\Vert \gamma \Vert }_{-1},$

${\Vert \nabla u\Vert }_0\le Ra{\kappa }^{-1}{\Vert \gamma \Vert }_{-1}.$ (2.7)

此外，若 $Pr,Ra,\kappa ,\gamma$ 满足下面的唯一性条件：

$0<\delta =\left(Pr-Ra{\kappa }^{-1}{\Vert \gamma \Vert }_{-1}\left(N+Pr\bar{N}{\kappa }^{-1}\right)\right)<Pr,$ (2.8)

则 $\left(u,p,T\right)$ 是方程(2.6)的唯一解。

2.2. 亏量校正方法

首先，我们给出有限元空间 $\left(X_h,M_h,W_h\right)$ ，

$X_h=\left\{v_h\in C^0\left(\Omega \right)\cap X;{v_h|}_T\in P_1{\left(T\right)}^2,\forall T\in {\tau }_h\right\},$

$M_h=\left\{q_h\in L^2{\left(\Omega \right)}^2\cap M;{q_h|}_T\in P_0{\left(T\right)}^2,\forall T\in {\tau }_h\right\},$

$W_h=\left\{v_h\in C^0\left(\Omega \right)\cap W;{v_h|}_T\in P_1\left(T\right),\forall T\in {\tau }_h\right\}.$

其中 $P_1\left(T\right)$ 表示T上的分片线性多项式， $P_0\left(T\right)$ 表示T的分片常数多项式。

我们给出 $L^2$ 正交投影的定义：

$\left({\rho }_{h}q,q_h\right)=\left(q,q_h\right),\forall q\in M,q_h\in M_h.$

显然有限元空间对 $\left(X_h,M_h,W_h\right)$ 是不满足离散的 $inf-sup$ 条件的，但是满足下面的性质A。

存在映射 $r_h:H^2{\left(\Omega \right)}^2\cap X\to X_h$ ， ${\bar{r}}_h:H^2\left(\Omega \right)\cap W\to W_h$ 和 ${\rho }_h:M\to M_h$ 满足

$\begin{array}{l}{\Vert \nabla \left(u-r_{h}u\right)\Vert }_0\le Ch{\Vert \Delta u\Vert }_0,\forall u\in H^2{\left(\Omega \right)}^2\cap X,\\ {\Vert p-{\rho }_{h}p\Vert }_0\le Ch{\Vert p\Vert }_1,\forall p\in H^1\left(\Omega \right)\cap M,\\ {\Vert \nabla \left(T-{\bar{r}}_{h}T\right)\Vert }_0\le Ch{\Vert \Delta T\Vert }_0,\forall T\in H^2\left(\Omega \right)\cap W,\end{array}$ (2.9)

根据上述定义给出原方程的离散形式：求解 $\left(u_h,p_h,T_h\right)\in X_h\times M_h\times W_h$ 满足

$\{\begin{array}{l}B\left(\left(u_h,p_h\right);\left(v_h,q_h\right)\right)+b\left(u_h;u_h,v_h\right)=PrRa\left(j{T}_h,v_h\right),\\ \kappa \bar{a}\left(T_h,s_h\right)+\bar{b}\left(u_h;T_h,s_h\right)=\left(\gamma ,s_h\right).\end{array}$ (2.10)

由于我们选取的是 $P_1-P_0-P_1$ 为不稳定元，所以考虑下面的罚跳跃稳定化方法：

$S_h\left(p_h,q_h\right)={\beta }_0{\displaystyle \underset{e\in {\Gamma }_T}{\sum }{h}_e{\displaystyle {\int }_e{\left[p_h\right]}_e{\left[q_h\right]}_e}},\forall p_h,q_h\in M_h.$

其中 ${\Gamma }_T$ 为 ${\tau }_h$ 的单元内边的集合， ${\left[q\right]}_e$ 为q在边e上的跳跃， ${\beta }_0$ 为给定的稳定参数。

$\{\begin{array}{l}{B}_h\left(\left(u,p\right);\left(v,q\right)\right)+b\left(u;u,v\right)=PrRa\left(jT,v\right),\\ \kappa \bar{a}\left(T,s\right)+\bar{b}\left(u;T,s\right)=\left(\gamma ,s\right).\end{array}$ (2.11)

这里的 $B_h\left(\left(u,p\right);\left(v,q\right)\right)=B\left(\left(u,p\right);\left(v,q\right)\right)+S_h\left(p_h,q_h\right)$。

接下来，我们将介绍亏量校正方法 [9] ，它包括：在亏量步求解人工粘性稳定非线性问题，然后用线性化方程校正几个步骤的亏量解。在亏量及校正的步骤，我们使用最稳定的Oseen迭代求解离散方程。

新方法描述如下：

考虑下面稳定问题：求解 $\left(u_h,p_h,T_h\right)\in \left(X_h,M_h,W_h\right)$ ，对任意的 $\left(v_h,q_h,s_h\right)\in \left(X_h,M_h,W_h\right)$ 有

$\{\begin{array}{l}{B}_h\left(\left(u_h,p_h\right);\left(v_h,q_h\right)\right)+b\left(u_h,u_h,v_h\right)+\alpha a\left(u_h,v_h\right)=PrRa\left(j{T}_h,v_h\right),\\ \left(\kappa +\alpha \right)\bar{a}\left(T_h,s_h\right)+\bar{b}\left(u_h;T_h,s_h\right)=\left(\gamma ,s_h\right),\end{array}$ (2.12)

其中 $\gamma >0$ 是一个稳定的人工粘性系数。我们定义它的迭代解 $\left(u_h^m,p_h^m,T_h^m\right)\in \left(X_h,M_h,W_h\right)$ ，对任意的 $\left(v_h,q_h,s_h\right)\in \left(X_h,M_h,W_h\right)$ 满足

$\{\begin{array}{l}{B}_h\left(\left(u_h^m,p_h^m\right);\left(v_h,q_h\right)\right)+b\left(u_h^{m-1},u_h^m,v_h\right)+\alpha a\left(u_h^m,v_h\right)=PrRa\left(j{T}_h^m,v_h\right),\\ \left(\kappa +\alpha \right)\bar{a}\left(T_h^m,s_h\right)+\bar{b}\left(u_h^{m-1};T_h^m,s_h\right)=\left(\gamma ,s_h\right),\end{array}$ (2.13)

这里的 $m=1,2,3,4,\cdots$。

离散方程的亏量 $D\left(u_h^m,p_h^m,T_h^m\right)$ 定义为

$\begin{array}{c}\left(D\left(u_h^m,p_h^m,T_h^m\right),v_h\right)=PrRa\left(j{T}_h,v_h\right)+\left(\gamma ,s_h\right)-B_h\left(\left(u_h,p_h\right);\left(v_h,q_h\right)\right)\\ -b\left(u_h,u_h,v_h\right)-\kappa \bar{a}\left(T_h^m,s_h\right)-\bar{b}\left(u_h^{m-1};T_h^m,s_h\right).\end{array}$ (2.14)

根据下列方程求解校正量 $\left({\epsilon }_h^0,{\theta }_h^0,{\xi }_h^0\right)$

$\{\begin{array}{l}{B}_h\left(\left(u_h^m,p_h^m\right);\left(v_h,q_h\right)\right)+b\left(u_h^{m-1},u_h^m,v_h\right)+\alpha a\left(u_h^m,v_h\right)=PrRa\left(j{T}_h^m,v_h\right),\\ \left(\kappa +\alpha \right)\bar{a}\left(T_h^m,s_h\right)+\bar{b}\left(u_h^{m-1};T_h^m,s_h\right)=\left(\gamma ,s_h\right)+\alpha \bar{a}\left(T_h^{m-1},s_h\right).\end{array}$ (2.15)

这样，离散方程的亏量校正方法可以描述如下：

步骤1：求解稳定有限元迭代解 $\left(u_h^m,p_h^m,T_h^m\right)\in \left(X_h,M_h,W_h\right)$ ，对 $\forall \left(v_h,q_h,s_h\right)\in \left(X_h,M_h,W_h\right)$ 有

$\{\begin{array}{l}{B}_h\left(\left(u_h^m,p_h^m\right);\left(v_h,q_h\right)\right)+b\left(u_h^{m-1},u_h^m,v_h\right)+\alpha a\left(u_h^m,v_h\right)=PrRa\left(j{T}_h^m,v_h\right),\\ \left(\kappa +\alpha \right)\bar{a}\left(T_h^m,s_h\right)+\bar{b}\left(u_h^{m-1};T_h^m,s_h\right)=\left(\gamma ,s_h\right),\end{array}$ (2.16)

其中 $m=1,2,3,\cdots$。

步骤2：求解校正解 $\left(u_{n+1}^h,p_{n+1}^h,T_{n+1}^h\right)\in \left(X_h,M_h,W_h\right)$ 满足

$\{\begin{array}{l}{B}_h\left(\left(u_{n+1}^h,p_{n+1}^h\right);\left(v_h,q_h\right)\right)+b\left(u_n^h,u_{n+1}^h,v_h\right)+\alpha a\left(u_{n+1}^h,v_h\right)=PrRa\left(j{T}_{n+1}^h,v_h\right)+\alpha a\left(u_n^h,v_h\right),\\ \left(\kappa +\alpha \right)\bar{a}\left(T_{n+1}^h,s_h\right)+\bar{b}\left(u_n^h;T_{n+1}^h,s_h\right)=\left(\gamma ,s_h\right)+\alpha \bar{a}\left(T_n^h,s_h\right),\end{array}$ (2.17)

其中 $\left(u_0^h,p_0^h,T_0^h\right)=\left(u_h^m,p_h^m,T_h^m\right),n=0,1,2,3,\cdots$。

流程图如下：

2.3. 稳定性和误差估计

我们给出上述方法的主要稳定性和误差估计结果。

定理2.1：假设 $\left(u_n^h,p_n^h,T_n^h\right)\in \left(X_h,M_h,W_h\right)$ 是方程(2.17)的解，满足：

${\Vert \nabla u_n^h\Vert }_0\le \frac{PrRa\left(1+\alpha {\kappa }^{-1}\right)}{\left(\kappa +\alpha \right)\left(Pr+\alpha \right)}{\Vert \gamma \Vert }_{-1},$

${\Vert \nabla T_n^h\Vert }_0\le \frac{1+\alpha {\kappa }^{-1}}{\kappa +\alpha }{\Vert \gamma \Vert }_{-1}.$

根据(2.4)式和(2.5)式，(2.8)式和(2.13)式可以很容易得到定理2.1的证明。

定理2.2：假设 $\left(u,p,T\right),\left(u_n^h,p_n^h,T_n^h\right)\in \left(X_h,M_h,W_h\right)$ 分别是方程(2.6)与(2.17)的解，C为正常数，则有下面估计：

${\Vert \nabla \left(u-u_n^h\right)\Vert }_0+{\Vert p-p_n^h\Vert }_0+{\Vert \nabla \left(T-T_n^h\right)\Vert }_0\le Ch+CPrRa\frac{1+\alpha {\kappa }^{-1}}{\kappa +\alpha }{\Vert \gamma \Vert }_{-1}+C\alpha \frac{PrRa\left(1+\alpha {\kappa }^{-1}\right)}{\left(\kappa +\alpha \right)\left(Pr+\alpha \right)}{\Vert \gamma \Vert }_{-1}.$

证明：令 $\left(e,\eta ,\xi \right)=\left(r_{h}u-u_n^h,{\rho }_{h}p-p_n^h,{\bar{r}}_{h}T-T_n^h\right)$。用(2.6)减去(2.17)且 $\left(v,q,s\right)=\left(v_h,q_h,\xi \right)$ 得

$\begin{array}{c}\kappa {\Vert \nabla \xi \Vert }_0\le \left(\kappa +\bar{N}{\Vert \nabla u_n^h\Vert }_0\right){\Vert \nabla \left({\bar{r}}_{h}T-T\right)\Vert }_0+\bar{N}{\Vert \nabla \left(r_{h}u-u\right)\Vert }_0{\Vert \nabla T\Vert }_0+\bar{N}{\Vert \nabla e\Vert }_0{\Vert \nabla T\Vert }_0+{\Vert \nabla T_n^h\Vert }_0+{\Vert \nabla T_{n+1}^h\Vert }_0\\ \le \left(\kappa +\bar{N}\frac{PrRa\left(1+\alpha {\kappa }^{-1}\right)}{\left(\kappa +\alpha \right)\left(Pr+\alpha \right)}{\Vert \gamma \Vert }_{-1}\right){\Vert \nabla \left({\bar{r}}_{h}T-T\right)\Vert }_0+\bar{N}\frac{1+\alpha {\kappa }^{-1}}{\kappa +\alpha }{\Vert \gamma \Vert }_{-1}{\Vert \nabla e\Vert }_0+C_0\frac{1+\alpha {\kappa }^{-1}}{\kappa +\alpha }{\Vert \gamma \Vert }_{-1}.\end{array}$

和

$\begin{array}{l}Pra\left(e,e\right)-d\left(e,p-{\rho }_h{p}_n^h\right)+d\left(u-r_h{u}_n^h,\eta \right)+b\left(e,u,e\right)+b\left(u_n^h,e,e\right)\\ =Pra\left(r_{h}u-u,e\right)+b\left(r_{h}u-u,u,e\right)+b\left(u_n^h,r_{h}u-u,e\right)\\ +PrRa\left(\left(j\xi ,e\right)+\left(j\left(T-{\bar{r}}_{h}T\right),e\right)\right)+\alpha a\left(u_{n+1}^h,e\right).\end{array}$

化简得

$\begin{array}{l}\left(Pr-P{r}^{-1}NRa{\kappa }^{-1}{\Vert \gamma \Vert }_{-1}-N\frac{PrRa\left(1+\alpha {\kappa }^{-1}\right)}{\left(\kappa +\alpha \right)\left(Pr+\alpha \right)}{\Vert \gamma \Vert }_{-1}\right){\Vert \nabla e\Vert }_0^2\\ \le \left(Pr-C_{0}N\frac{PrRa\left(1+\alpha {\kappa }^{-1}\right)}{\left(\kappa +\alpha \right)\left(Pr+\alpha \right)}{\Vert \gamma \Vert }_{-1}\right){\Vert \nabla \left(r_{h}u-u\right)\Vert }_0{\Vert \nabla e\Vert }_0\\ +PrRa{\Vert \nabla \xi \Vert }_0{\Vert \nabla e\Vert }_0+{\Vert \nabla \left(T-{\bar{r}}_{h}T\right)\Vert }_0{\Vert \nabla e\Vert }_0\\ +\alpha \frac{PrRa\left(1+\alpha {\kappa }^{-1}\right)}{\left(\kappa +\alpha \right)\left(Pr+\alpha \right)}{\Vert \gamma \Vert }_{-1}{\Vert \nabla e\Vert }_0.\end{array}$

根据Young不等式可以得到

${\Vert \nabla e\Vert }_0^2+{\Vert \nabla \eta \Vert }_0^2+{\Vert \nabla \xi \Vert }_0^2\le Ch+C\alpha \frac{PrRa\left(1+\alpha {\kappa }^{-1}\right)}{\left(\kappa +\alpha \right)\left(Pr+\alpha \right)}{\Vert \gamma \Vert }_{-1}+CPrRa\frac{1+\alpha {\kappa }^{-1}}{\kappa +\alpha }{\Vert \gamma \Vert }_{-1}.$

3. 基于亏量校正方法的恢复型误差估计子

3.1. 恢复型误差估计子

本小节中，我们将构造出一个恢复型误差估计子并分析它的性质。假设 $\left(u_h,p_h,T_h\right)$ 数值结果，我们考虑应力张量 $\sigma :=\nabla u-pI+\nabla T$ ，和它的有限元近似 ${\sigma }_n^h:=\nabla u_n^h-p_n^{h}I+\nabla T_n^h$。

其中I是 $2\times 2$ 的单位矩阵。Zienkienwicz-Zhu估计子的主要思想是将间断有限元梯度进行恢复得到一个连续的恢复项G (参考 [10] [11] )。

在 ${\tau }_h$ 中分别定义 $N,N_h,{\varphi }_v$ 和 ${\omega }_v$ 为：在 ${\tau }_h$ 中所有顶点的集合，在 ${\tau }_h$ 中的在 $\Omega$ 内部所有顶点的集合，v ( $\forall v\in N$ )的基函数共享一个顶点v的所有单元的集合 $\left({\omega }_v:=\text{supp}{\varphi }_z\right)$。

结合文献 [11] 中的超收敛片恢复技巧将分片常数张量 ${\sigma }_n^h$ 恢复到连续的从而获得一个恢复应力张量 $G\left({\sigma }_n^h\right)$。

对任意顶点 $v\in N$ 和它的片 ${\omega }_v$ ，定义

$G\left({\sigma }_n^h\right)\left(v\right)={\displaystyle {\sum }_{T\in {\omega }_v}\left(\left|T\right|{\sigma }_n^h\right)/\left|{\omega }_v\right|},$

其中 $\left|T\right|$ 是三角剖分T的面积。细节可以参考文献 [12] 。

构造一个恢复型误差估计子。定义局部误差和全局误差：

${\eta }_T:={\Vert {\sigma }_n^h-G\left({\sigma }_n^h\right)\Vert }_T,\eta :={\left({\displaystyle \underset{T\in {\tau }_h}{\sum }{\eta }_T^2}\right)}^{\frac{1}{2}}.$

恢复型误差估计子 ${\eta }_T$ 有以下几个重要性质。证明方式与文献( [12] [13] )相似，不再赘述。

引理3.1：根据上述 ${\sigma }_n^h$ 和 $G\left({\sigma }_n^h\right)$ 的定义，存在两个与网格无关的常数 $C_1$ 和 $C_2$ 有

$C_1\left({\Vert \left[{\sigma }_n^h\cdot n_e\right]\Vert }_{S_h}+{\Vert {\left[p_n^h\right]}_e\Vert }_{S_h}+{\Vert {\left[\nabla T_n^h\right]}_e\Vert }_{S_h}\right)n_e$ $\le \Vert {\sigma }_n^h-G\left({\sigma }_n^h\right)\Vert \le C_2\left({\Vert \left[{\sigma }_n^h\cdot n_e\right]\Vert }_{S_h}+{\Vert {\left[p_n^h\right]}_e\Vert }_{S_h}+{\Vert {\left[\nabla T_n^h\right]}_e\Vert }_{S_h}\right).$

引理3.2：根据上述 ${\sigma }_n^h$ 和 $G\left({\sigma }_n^h\right)$ 的定义，存在一个与网格无关的常数C，可以得到以下 ${\Vert \left[\nabla u_h\cdot n_e\right]\Vert }_{S_h}$ 和 ${\Vert {\left[\nabla T_h\right]}_e\Vert }_{S_h}$ 误差估计

$C\left({\Vert \left[\nabla u_n^h\cdot n_e\right]\Vert }_{S_h}+{\Vert {\left[p_n^h\right]}_e\Vert }_{S_h}+{\Vert {\left[\nabla T_n^h\right]}_e\Vert }_{S_h}\right)\le {\Vert \left[\nabla {\sigma }_n^h\cdot n_e\right]\Vert }_{S_h}+{\Vert {\left[p_n^h\right]}_e\Vert }_{S_h}$ .

3.2. 可靠性和有效性分析

定理3.1：假设 $\left(u,p,T\right)$ 和 $\left(u_n^h,p_n^h,T_n^h\right)$ 分别是(2.6)和(2.17)的解，存在一个与网格无关的正常数C有

${\Vert \nabla \left(u-u_n^h\right)\Vert }_0+{\Vert p-p_n^h\Vert }_0+{\Vert \nabla \left(T-T_n^h\right)\Vert }_0\le C\left(\eta +{\left|\log h_{T,\min }\right|}^{\frac{1}{2}}{\eta }^2\right).$

证明：存在一个强制算子 $Q_h:X\times M\times W\to X_h\times M_h\times W_h$ 有

$Q_h\left(v,q,s\right)=\left(I_{h}v,0,I_{h}s\right).$

因此，根据引理3.1，引理3.2，用Cauchy-Schwarz不等式，有

$\begin{array}{l}{\langle B_h\left(\left(u_n^h,p_n^h,T_n^h\right)\right),\left(I_{h}d-Q_h\right)\left(v,q,s\right)\rangle }_{X_h\times M_h\times W_h}\\ =Pra\left(u_n^h,\left(I_{h}d-Q_h\right)v\right)-d\left(\left(I_{h}d-Q_h\right)v,p_n^h\right)\\ +d\left(u_n^h,\left(I_{h}d-Q_h\right)q\right)+b\left(u_n^h;u_n^h,\left(I_{h}d-Q_h\right)v\right)\\ +\kappa \bar{a}\left(T_n^h,\left(I_{h}d-Q_h\right)s\right)+\bar{b}\left(u_n^h;T_n^h,\left(I_{h}d-Q_h\right)s\right)\\ -PrRa\left(j{T}_n^h,\left(I_{h}d-Q_h\right)v\right)-\left(\gamma ,\left(I_{h}d-Q_h\right)s\right)\end{array}$

$\begin{array}{l}\le {\displaystyle \underset{e\in S_h}{\sum }{\left(\left[{\sigma }_h\cdot n_e\right],v-I_{h}v\right)}_e}+C\eta {\Vert q\Vert }_0+C\eta {\Vert v\Vert }_1+b\left(u_n^h;u_n^h,\left(I_{h}d-Q_h\right)v\right)\\ +\bar{b}\left(u_n^h;T_n^h,\left(I_{h}d-Q_h\right)s\right)+\left(\gamma ,\left(Id-Q_h\right)s\right)\\ \le C\left(\eta +{\left|{\log }_{h_T,\min }\right|}^{\frac{1}{2}}{\eta }^2\right){\Vert \nabla v\Vert }_0+{\Vert q\Vert }_0+{\Vert \nabla s\Vert }_0.\end{array}$

根据 $B_h$ 的定义，近似解 $\left(u_n^h,p_n^h,T_n^h\right)$ 满足 ${\langle B_h\left(\left(u_n^h,p_n^h,T_n^h\right)\right),Q_h\left(v,q,s\right)\rangle }_{X_h\times M_h\times W_h}=0$。

同时满足 ${\langle B\left(\left(u_n^h,p_n^h,T_n^h\right)\right),Q_h\left(v,q,s\right)\rangle }_{X\times M\times W}=0$。所以，结合引理，可得一下估计

${\Vert \nabla \left(u-u_n^h\right)\Vert }_0+{\Vert p-p_n^h\Vert }_0+{\Vert \nabla \left(T-T_n^h\right)\Vert }_0\le C\left(\eta +{\left|\log h_{T,\min }\right|}^{\frac{1}{2}}{\eta }^2\right).$

定理3.1得证。

从引理3.1，可得 $\eta \le C\left({\Vert \left[{\sigma }_n^h\cdot n_e\right]\Vert }_{S_h}+{\Vert {\left[p_n^h\right]}_e\Vert }_{S_h}+{\Vert {\left[\nabla T_n^h\right]}_e\Vert }_{S_h}\right)$。

为了证明恢复型误差估计子 $\eta$ 的有效性，首先需要估计：

${\Vert \left[{\sigma }_n^h\cdot n_e\right]\Vert }_{S_h},{\Vert {\left[p_n^h\right]}_e\Vert }_{S_h}$ 和 ${\Vert {\left[\nabla T_n^h\right]}_e\Vert }_{S_h}.$

引理3.3：假设 $\left(u,p,T\right)$ 和 $\left(u_n^h,p_n^h,T_n^h\right)$ 分别是(2.6)和(2.17)的解。存在一个与网格无关的正常数C有

${\Vert \left[{\sigma }_n^h\cdot n_e\right]\Vert }_{S_h}\le C{\Vert \nabla \left(u-u_n^h\right)\Vert }_0+{\Vert p-p_n^h\Vert }_0+{\Vert \nabla \left(T-T_n^h\right)\Vert }_0.$

证明：从 ${\sigma }_n^h$ 的定义出发，可得

${\Vert \left[{\sigma }_n^h\cdot n_e\right]\Vert }_{S_h}\le C{\Vert \left[\left(\nabla u_n^h-p_n^{h}I\right)\cdot n_e\right]\Vert }_{S_h}+C{\Vert \left[\nabla T_n^h\cdot n_e\right]\Vert }_{S_h}\le C{\Vert \nabla \left(u-u_n^h\right)\Vert }_0$ $+{\Vert p-p_n^h\Vert }_0+{\Vert \nabla \left(T-T_n^h\right)\Vert }_0.$

下面定义插分算子 $\Pi :L^2\left(\Omega \right)\to R_{1,h}\left(\Omega \right)$。用迹不等式有

${\Vert {\left[\nabla T_n^h\right]}_e\Vert }_{S_h}\le {\Vert {\left[\nabla T_n^h-\Pi \nabla T\right]}_e\Vert }_{S_h}\le C(\Vert \nabla T-\nabla T_n^h\Vert$ $+\Vert \nabla T-\Pi \nabla T\Vert )\le C{\Vert \nabla \left(u-u_n^h\right)\Vert }_0+{\Vert p-p_n^h\Vert }_0+{\Vert \nabla \left(T-T_n^h\right)\Vert }_0.$

${\Vert {\left[p_n^h\right]}_e\Vert }_{S_h}\le {\Vert {\left[p_n^h-\Pi p\right]}_e\Vert }_{S_h}$ $\le C\left(\Vert p-p_n^h\Vert +\Vert p-\Pi p\Vert \right)\le C{\Vert \nabla \left(u-u_n^h\right)\Vert }_0+{\Vert p-p_n^h\Vert }_0+{\Vert \nabla \left(T-T_n^h\right)\Vert }_0.$

定理3.2：假设 $\left(u,p,T\right)$ 和 $\left(u_n^h,p_n^h,T_n^h\right)$ 分别是(2.6)和(2.17)的解，存在一个与网格无关的常数C有

$\eta \le C{\Vert \nabla \left(u-u_n^h\right)\Vert }_0+{\Vert p-p_n^h\Vert }_0+{\Vert \nabla \left(T-T_n^h\right)\Vert }_0.$

证明：结合上面的估计，用引理3.1和引理3.2，可得

$\eta \le C\left({\Vert \left[{\sigma }_n^h\cdot n_e\right]\Vert }_{S_h}+{\Vert {\left[p_n^h\right]}_e\Vert }_{S_h}+{\Vert {\left[\nabla T_n^h\right]}_e\Vert }_{S_h}\right)$ $\le C{\Vert \nabla \left(u-u_n^h\right)\Vert }_0+{\Vert p-p_n^h\Vert }_0+{\Vert \nabla \left(T-T_n^h\right)\Vert }_0.$

完成了有效性的证明。

4. 数值实验

在本小节中，结合自适应方法给出两个个数值算例来验证我们所提方法的有效性。对于处理非线性部分，我们均采用Ocien迭代。

为了方便起见，下面给出几个定义：

(i) $DOF:=$ 表是三角剖 ${\tau }_h^j$ 的自由度个数，

(ii) $e_r:={\left({\Vert u-u_h\Vert }_1^2+{\Vert p-p_h\Vert }_0^2+{\Vert T-T_h\Vert }_1^2\right)}^{\frac{1}{2}}/{\left({\Vert u\Vert }_1^2+{\Vert p\Vert }_0^2+{\Vert T\Vert }_1^2\right)}^{\frac{1}{2}},$

(iii) ${\eta }_r:=\eta /{\left({\Vert u\Vert }_1^2+{\Vert p\Vert }_0^2+{\Vert T\Vert }_1^2\right)}^{\frac{1}{2}},$

(iv) #Math_205#

4.1. 光滑解问题

第一个算例我们求解 $\Omega ={\left[0,1\right]}^2$ 区域上的一个光滑真解问题，从而验证我们的方法，对于光滑真解是有效的。我们给出真解 $u=\left(u_1,u_2\right)$ ，压力p和温度T如下：

$u_2\left(x_1,x_2\right)=-10x_1\left(x_1-1\right)\left(2x_1-1\right)x_2^2{\left(x_2-1\right)}^2,$

$p\left(x_1,x_2\right)=10\left(2x_1-1\right)\left(2x_2-1\right),$

$T\left(x_1,x_2\right)=u_1\left(x_1,x_2\right)+u_2\left(x_1,x_2\right).$

从表1中可以看出，基于亏量校正方法的后验误差估计是适用的。 $I_{eff}$ 趋于1表示我们的估计子与真实误差是有效渐进的，说明我们的方法是有效的。图1中(a)是初始网格，(b)是自适应加密后的网格，从图中可以看出，加密是有效的。

4.2. 方腔流

第二个数值实验我们考虑没有真解的方腔流问题。定义计算区域 $\Omega ={\left[0,1\right]}^2$。取 $Pr=0.71$ ， $\kappa =1$ ， $\alpha =h^2/m$ ，边界定义如下：左边界和底部取 $T=0$ ，底边取 $\partial T/\partial n=0$ ，其余取 $T=4y\left(1-y\right)$ ，且速度上

Table 1. The result of smooth solution with Pr = 0.71, κ = 1, Ra = 1000

表1. 光滑解结果，其中Pr = 0.71，κ = 1，Ra = 1000

(a) (b)

Figure 1. From left to right: original meth, refined meth

图1. 从左到右分别是初始网格和网格加密结果

(a) (b)

Figure 2. From left to right: original meth, refined meth

图2. 从左到右分别是初始网格和网格加密结果

(a) (b) (c)

Figure 3. The numerical streamlines with Ra = 10, Ra = 1000, Ra = 1e5

图3. 数值结果的速度流线图：Ra = 10，Ra = 1000，Ra = 1e5

(a) (b) (c)

Figure 4. The numerical isobars with Ra = 10, Ra = 1000, Ra = 1e5

图4. 数值结果的等压图：Ra = 10，Ra = 1000，Ra = 1e5

为0的Dirichlet条件。图1表示了自适应方法的初始网格和自适应后的网格。

图2表示初始网格和自适应加密后的网格，表示自适应方法在捕捉大误差处是有效的，加密结果是我们想要的。图3，图4分别表示不同Rayleigh数10，10³，10⁵的速度流线图和等温图。容易看出，我们的方法在求解没有真解的方腔流问题同样有效，其次在求解大Rayleigh方面比一般有限元更加有效。

5. 结论与展望

本文提出了基于亏量校正方法的恢复型误差估计子，通过数值实验，验证了我们方法的有效性和可靠性。同时，在求解方腔流大Rayleigh数问题时，由于在亏量步添加了人工黏性系数，使我们避免了对流占优情况的发生，而后的校正步与后验误差估计子使近似解更加地逼近真解，这样得出的结果是优于一般有限元方法的。此方法可以推广到其他流体的问题。

基金项目

不可压缩磁流体动力学方程的非协调自适应有限元高效算法研究(11701493)。

参考文献