满足速度压力 $\inf$-$\sup$ 条件的多种混合有限元算法^{[1, 2]}已经被深入研究.$\inf$-$\sup$ 条件使得更理想的有限元空间不能使用,比如低阶 $P_1/P_1$ 有限元对.实际科学计算中,低次等阶 $P_1/P_1$ 有限元对由于其数据结构简单、计算量小且精度高等优点更受到欢迎而被广泛使用.但核心问题是 $P_1/P_1$ 有限元对不满足离散的 $\inf$-$\sup$ 条件会导致压力震荡.

为了解决这些困难,科学家们致力于研究稳定化技巧. 因此近些年里,很多稳定化技巧(参见文献[3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27])被用于补偿相容性条件.近期研究工作主要集中在几种稳定化技巧上,包括罚方法^{[24, 25]}、正则化方法^{[27, 28, 29]}、局部 Gauss 积分法^{[10, 11, 12]} 和多尺度增量方法^{[19, 20, 21]}.%(参见文献[19, 20, 21]).

虽然求解 Navier-Stokes 方程的稳定化有限元算法已经有很多理论结果,但关于这些算法之间的比较及数值解行为的研究工作有限.这些算法在求解 Stokes 方程^{[30, 31]}和二阶偏微分方程^{[32, 33]} 时的数值行为研究工作已经进行.本文的目的在于研究 Navier-Stokes 方程四种关于等阶 $P_1/P_1$ 有限元对稳定化算法的数值行为,这四种算法分别包括罚方法、正则化方法、局部 Gauss 积分法和多尺度增量方法. 对这四种算法的比较工作从四个方面入手,分别是稳定化参数、收敛速度、耗时和相容性. 数值试验表明在计算大雷诺数问题时多尺度稳定化算法具有好的稳定性、精度及计算效率,因此对求解大雷诺数 Navier-Stokes 方程时多尺度稳定化有限元算法是一种不错的选择.

本文的结构安排如下: 第二部分介绍 Sobolev 函数空间和 Navier-Stokes 方程.接下来在第三部分列出求解 Navier-Stokes 方程的几种常见稳定化有限元算法,详细介绍了各个算法的稳定技巧. 在第四部分给出了数值算例,通过试验客观比较各个算法的数值行为. 最后一部分对文章的主要研究结果进行总结.

假设 $\Omega$ 是 ${\Bbb R}^{2}$ 中一个具有 Lipschitz -连续边界 $\partial\Omega$ 的有界凸区域.本文研究如下具有 Dirichlet 边界条件的定常 Navier-Stokes 方程

$\begin{equation}\left\{\begin{array}{lll}-\nu\Delta +(u\cdot\nabla)u+\nabla p=f,&\ \ \ \mbox{在$ \Omega$中,}\\\nabla \cdot u=0,&\ \ \ \mbox{在$ \Omega$中,}\\u=0,&\ \ \ \mbox{在$\partial\Omega$上,}\end{array}\right.\end{equation}$

(2.1)

其中 $\nu>0$ 是粘性系数,$u=(u_{1}(x_{1},x_{2}),u_{2}(x_{1},x_{2}))$ 表示速度,$p=p(x_{1},x_{2})$ 为压力,$f=(f_{1}(x_{1},x_{2}),f_{2}(x_{1},x_{2}))$ 则是右端体积力.

令 \begin{eqnarray*}X\triangleq(H^{1}_{0}(\Omega))^{2} ,\ \ \ Y\triangleq(L^{2}(\Omega))^{2} ,\ \ \ M\triangleq L^{2}_{0}(\Omega)=\{q\in L^{2}(\Omega): \int_{\Omega}q {\rm d}x=0\}.\end{eqnarray*} 它们皆是 Sobolev 函数空间. 进一步,$(L^{2}(\Omega))^{m}$ ($m=1,2,4$ ) 具有 $L^{2}$ -内积 $(.,.)$ 和 $L^{2}$ -范数$\|\cdot\|_{0}$;$X$ 赋予如下通常的内积和范数

$((u,v))\triangleq {{((\nabla u,\nabla \text{v}))}_{\Omega }},\ \ \ |u{{|}_{1}}\triangleq {{((u,u))}^{\frac{1}{2}}}=(\nabla u,\nabla \text{v})_{\Omega }^{\frac{1}{2}}.$

同时,我们给出一般的 Sobolev 函数空间 $W^{m,p}(\Omega)$,其具有全范数 $\|\cdot\|_{m,p}$ 和半范 $|\cdot|_{m,p}$ ($m,p\geq0$).一般约定,用 $H^{m}(\Omega)$ 表示 $W^{m,2}(\Omega)$,$\|\cdot\|_{m}$ 简记 $\|\cdot\|_{m,2}$,$|\cdot|_{m}$ 简记 $|\cdot|_{m,2}$.

另外,我们引入如下三线性项

$b(u;v,w)=((u\cdot\nabla)v,w)_{\Omega},\ \ \ \forall u,v,w\in X,$

(2.2)

该三线性形式 $b(\cdot;\cdot,\cdot)$ 具有如下性质(参见文献 ^{[1, 34]}):

$|b(u;v,w)|\leq C|u|_{1}|v|_{1}|w|_{1},\ \ \ \forall u,v,w\in X,$

(2.3)

$b(u;v,w)=-{{((u\cdot \nabla )w,v)}_{\Omega }}-{{(\nabla \cdot u,v\cdot w)}_{\Omega }},\ \ \ \forall u,v,w\in X,$

(2.4)

$b(u;v,v)=-\frac{1}{2}{{(\nabla \cdot u,v\cdot v)}_{\Omega }},\ \ \ \forall u,v,w\in X,$

(2.5)

其中 $C$ 为仅依赖于区域 $\Omega$ 的正常数.

此时方程(2.1) 的弱形式表示为: 找 $(u,p)\in X\times M$ 使得

$a(u,v)-d(p,v)+d(q,u)+b(u;u,v)=(f,v)_{\Omega},\ \ \ \ \forall (v,q)\in X\times M,$

(2.6)

其中 $a(u,v)\triangleq\nu (\nabla u,\nabla v)_{\Omega}$,$\ d(p,v)\triangleq(p,\nabla\cdot v)_{\Omega}.$且双线性形式 $d(\cdot,\cdot)$ 满足如下的 $\inf-\sup$ 条件(参见文献 ^{[1, 34]})

$\sup_{v\in X}\frac{d(q,v)}{\|\nabla \text{v}\|_0}\geq\beta\|q\|_0,\ \ \ \ \forall q\in M.$

(2.7)

这一小节里我们给出方程 (2.1) 的有限元逼近形式. 对区域 $\Omega$ 进行三角剖分 (记为 $\tau_{h}$),每一个剖分单元由具有边界 $\partial K$ 的闭三角单元 $K$ 构成.$\{p_{j},j=1,2,\cdots ,N\}$ 表示剖分 $\tau_{h}$ 的所有节点的集合.$\varepsilon_{h}$ 标记该剖分 $\tau_{h}$ 所有边的集合,而 $\varepsilon^{int}_{h}$ 标记与 $\partial\Omega$ 没有交集的所有内边的集合.进一步,$h_{K}$ 为单元 $K$ 的直径,$h=\max\limits_{K\in\tau_{h}}{h_{K}}$ 是网格参数.假设剖分 $\tau_{h}$ 是正则且拟一致的剖分,即存在 $C>0$ 使得

$\frac{{{h}_{K}}}{{{\rho }_{K}}}\le C,\ \ \forall \ K\in {{\tau }_{h}},$

(3.1)

其中 $\varrho_{K}$ 表示含于单元 $K$ 内最大球的直径.$X_{h}$ 和 $M_{h}$ 分别表示基于剖分 $\tau_{h}$ 的速度和压力有限元逼近空间. 其定义如下 $$X_{h}\triangleq\{v\in X\cap (C^{0}(\overline{\Omega}))^{2}:v|_{K}\in P_{1}(K)^{2},~~ \forall\ K \in \tau_{h}\}, $$ $$M_{h}\triangleq\{q \in M: q|_{K}\in P_{1}(K),~~ \forall\ K \in \tau_{h}\}, $$ 其中 $P_{1}(K)$ 表示单元 $K$ 上线性或者双线性函数空间.

此时,可以给出 Navier-Stokes 方程经典有限元 Galerkin 变分形式如下:寻找 $(u_{h},\ p_{h})\in X_{h}\times M_{h}$ 使得

$a(u_{h},v_{h})-d(p_{h},v_{h})+d(q_{h},u_{h})+b(u_{h};u_{h},v_{h})=(f,v_{h})_{\Omega},~~\forall \ (v_{h},q_{h})\in X_{h}\times M_{h}.$

(3.2)

接下来,用 $U$ 表示速度在节点处值形成的数组,$P$ 表示压力在节点处值形成的数组且 $F$ 表示右端体积力在节点处取值形成的数组.这时 (3.2) 式可以写为如下的矩阵方程形式

$\left(\begin{array}{cc}A+N(x)~~ & B\\B^T~~ & 0\end{array}\right)\left(\begin{array}{c}UP \end{array}\right)=\left(\begin{array}{c}F0 \end{array}\right),$

其中矩阵 $A$ 和 $B$ 是运用 $X_h$ 和 $M_h$ 两个空间的基,分别从双线性形式 $a(\cdot,\cdot)$ 及 $d(\cdot,\cdot)$ 推导出来.$B^T$ 表示矩阵 $B$ 的转置. 矩阵 $N(x)$ 是由非线性项 $b(\cdot;\cdot,\cdot)$ 推导所得,且为关于变量 $x=(x_1,x_2)$ 的函数. 再注意到等阶 $P_1/P_1$ 有限元对不满足如下离散的 $\inf$-$\sup$ 条件^[35]

$\sup_{v_{h}\in X_h}\frac{d(v_{h},q_h)}{\|\nabla \text{v}_{h}\|_0}\geq \beta_1\|q_h\|_0,%\mbox{or} \sup_{v_{h}\in NC_h}\frac{d_h(v_{h},q_h)}{\|\nabla \text{v}_{h}\|_{0,h}}\geq \beta_1\|q_h\|_0 \forall q_h\in M_h,$

其中 $\beta_1$ 依赖于 $h$. 为了克服此困难,好多种稳定化技巧被用于稳定等阶 $P_1/P_1$ 有限元配对.

方法 1 (加罚方法^{[24, 25]}) 寻找 $(u_{h},\ p_{h})\in X_{h}\times M_{h}$ 使得

$a(u_{h},v_{h})-d(p_{h},v_{h})+d(q_{h},u_{h})+\frac{\varepsilon}{\nu}(p_h,q_h)_{\Omega}+b(u_{h};u_{h},v_{h})=(f,v_{h})_{\Omega},$

(3.3)

其中对任意的 $(v_{h},q_{h})\in X_{h}\times M_{h}$,$\varepsilon>0$ 是罚参数.加罚方法相容性弱使得其计算精度不高,不能很快的收敛,达不到理论分析结果^{[1, 36, 37]} 中所期望达到的最优收敛阶.此外,该方法数值解行为明显依赖于罚参数 $\varepsilon$ 的选取. 其变分形式 (3.3) 可以表示为如下的矩阵形式

$\left( \begin{matrix} A+N(x)~~ & B \\ {{B}^{T}}~~ & \frac{\varepsilon }{\nu }D \\ \end{matrix} \right)\left( \begin{matrix} UP \\ \end{matrix} \right)=\left( \begin{matrix} F0 \\ \end{matrix} \right),$

其中矩阵 $D$ 是由运用 $M_h$ 的基,从稳定项 $(p_h,q_h)$ 推导得来.

方法 2 (正则性方法^{[28, 29]}) 寻找 $(u_{h},\ p_{h})\in X_{h}\times M_{h}$ 使得

$\begin{align} & a({{u}_{h}},{{v}_{h}})-d({{p}_{h}},{{v}_{h}})+d({{q}_{h}},{{u}_{h}})+{{\tau }_{1}} \\ & \sum\limits_{K\in {{\tau }_{h}}}{\frac{h_{k}^{2}}{\nu }}{{(\nabla {{p}_{h}},\nabla {{q}_{h}})}_{K}}+b({{u}_{h}};{{u}_{h}},{{v}_{h}})={{(f,{{v}_{h}})}_{\Omega }}+ \\ & {{\tau }_{1}}\sum\limits_{K\in {{\tau }_{h}}}{\frac{h_{k}^{2}}{\nu }}{{(f,\nabla {{q}_{h}})}_{K}},\forall \ ({{v}_{h}},{{q}_{h}})\in {{X}_{h}}\times {{M}_{h}},\\ \end{align}$

(3.4)

其中 $\tau_1$ 是依赖于网格尺寸 $h$ 的稳定化参数. (3.4) 式为不相容的,其矩阵形式

$\left(\begin{array}{cc}A+N(x)~~ & B\\B^T ~~& \tau_1\frac{h^2_k}{\nu}E\end{array}\right)\left(\begin{array}{c}UP \end{array}\right)=\left(\begin{array}{c}F \tau_1\frac{h^2_k}{\nu}F_1\end{array}\right),$

其中矩阵 $E$ 为利用 $M_h$ 的基,从稳定项 $\sum\limits_{K\in \tau_h}(\nabla p_h,\nabla q_h)_K$ 推导所得.矩阵 $F_1$ 是由稳定化项 $\sum\limits_{K\in \tau_h}(f,\nabla q_h)_K$ 推导得出.

方法 3 (局部 Gauss 积分法^{[10, 11, 12]}) 寻找 $(u_{h},\ p_{h})\in X_{h}\times M_{h}$ 使得

$a(u_{h},v_{h})-d(p_{h},v_{h})+d(q_{h},u_{h})+G(p_h,q_h)+b(u_{h};u_{h},v_{h})=(f,v_{h})_{\Omega},$

(3.5)

其中 $\forall (v_{h},q_{h})\in X_{h}\times M_{h}$,$G(p_h,q_h)$ 定义为

$\begin{eqnarray*}G(p_h,q_h)=\sum\limits_{K\in \tau_h}\bigg\{\int_{K,2}p_hq_h {\rm d}x-\int_{K,1}p_hq_h {\rm d}x\bigg\},\forall p_h,q_h\in M_h,\end{eqnarray*}$

这里 $\int_{K,i}q(x){\rm d}x$ 表示单元 $K$ 上的一个局部 Gauss 积分,当多项式的阶数 $j$ 取值 $j=1,2$ 时该积分可以取到精确解.其中 $q(x)=p_hq_h$ 为阶数不高于2的多项式. (3.5)式矩阵形式为

$\left(\begin{array}{cc}A+N(x) & B\\B^T & G\end{array}\right)\left(\begin{array}{c}UP \end{array}\right)=\left(\begin{array}{c}F0 \end{array}\right),$

其中矩阵 $G$ 为利用 $M_h$ 基,由稳定化项 $G(p_h,q_h)$ 推导得到.

方法 4 (多尺度增量方法^[19]) 寻找 $(u_{h},\ p_{h})\in X_{h}\times M_{h}$ 使得

$\begin{align} & a({{u}_{h}},{{v}_{h}})-d({{p}_{h}},{{v}_{h}})+d({{q}_{h}},{{u}_{h}})+b({{u}_{h}};{{u}_{h}},{{v}_{h}}) \\ & +\sum\limits_{K\in {{\tau }_{h}}}{{{\tau }^{RFB}}}{{(-\nu \nabla {{\text{u}}_{h}}+{{u}_{h}}\nabla {{u}_{h}}+\nabla {{p}_{h}},\nu \nabla {{\text{v}}_{h}}+\nabla {{\text{v}}_{h}}{{u}_{h}}+\nabla {{q}_{h}})}_{K}} \\ & +{{\tau }_{2}}\sum\limits_{E\in \varepsilon _{h}^{int}}{\frac{{{h}_{E}}}{\nu }}{{({{[[\nu {{\partial }_{n}}{{u}_{h}}]]}_{E}},{{[[\nu {{\partial }_{n}}{{v}_{h}}]]}_{E}})}_{E}}= \\ & {{(f,{{v}_{h}})}_{\Omega }}+\sum\limits_{K\in {{\tau }_{h}}}{{{\tau }^{RFB}}}{{(f,\nu \nabla {{\text{v}}_{h}}+\nabla {{\text{v}}_{h}}{{u}_{h}}+\nabla {{q}_{h}})}_{K}},\\ \end{align}$

(3.6)

这里$\forall {{\text{v}}_{h}}\in {{X}_{h}}$,$\ q_{h}\in M_{h}$. 稳定化参数由如下式子给出

$\tau^{RFB}=\left\{ \begin{array}{ll} \frac{h^{2}_{K}}{12\nu},&\mbox{若 }~ 0\leq P_{e_{K}}<1,\\[3mm] \frac{h_{K}}{2|u_{h}|_{K}},~~ &\mbox{若 }~ P_{e_{K}}\geq1,\end{array} \right.$

(3.7)

其中 $P_{e_{K}}=\frac{|u_{h}|_{K}h_{K}}{6\nu}$ 为 Peclet 数,$|u_{h}|_{K}:=\frac{\|u_{h}\|_{0,K}}{|K|^{\frac{1}{2}}}$.$\tau_2$ 为不依赖于网格尺寸 $h$ 的稳定化参数. (3.6) 式保持相容性,其矩阵形式如下

$\left(\begin{array}{cc}A+N(x) & B\\B^T+\tau_2B_2+\tau^{RFB}N_1(x)+\tau^{RFB}N_2(x) & \tau^{RFB}H+\tau^{RFB}H_1\end{array}\right)\left(\begin{array}{c}UP \end{array}\right)=\left(\begin{array}{c}F+\tau^{RFB}I(x)\\\tau^{RFB}I\end{array}\right),$

其中矩阵 $B_2,N_1(x)$ 和 $N_2(x)$ 分别从项 \begin{eqnarray*}\sum\limits_{E\in \varepsilon^{int}_{h}}\frac{h_{E}}{\nu}([[\nu\partial_{n}u_{h}]]_{E},[[\nu \partial_{n}v_{h}]]_{E})_{E},\sum\limits_{K\in \tau_{h}}(u_{h}\nabla u_{h},\nabla \text{v}_{h}u_{h})_K,\sum\limits_{K\in \tau_{h}}(u_{h}\nabla u_{h},\nabla q_h)_K,\end{eqnarray*}推导所得; $H$ 和 $H_1$ 分别由项 $\sum\limits_{K\in \tau_{h}}(\nabla p_h,\nabla q_h)_K$ 和 $\sum\limits_{K\in \tau_{h}}(\nabla p_h,\nabla \text{v}_{h}u_{h})_K$ 得到;而矩阵 $I(x)$ 则由项 $\sum\limits_{K\in \tau_{h}}(f,\nabla \text{v}_{h}u_{h})_K$ 推导得到;最后 $I$ 由从项 $\sum\limits_{K\in \tau_{h}}(f,\nabla q_h)_K$ 所得.

这一小节里,我们对上一节所给的四种稳定化算法给出数值算例进行比较.所有的非线性项采用 Oseen 迭代,各个算法皆在同一平台下,由有限元软件 Freefem++^[38]实现.

在区域 $\Omega=[0,1]\times[0,1],\nu=0.1$ 条件下,考虑 Navier-Stokes 方程 (2.1),且右端体积力 $f$ 通过原方程代入如下真解求得 $$u_{1}(x_{1},x_{2})=10x_{1}^{2}(x_{1}-1)^2x_{2}(x_{2}-1)(2x_{2}-1),p(x_{1},x_{2})= 10(2x_{1}-1)(2x_{2}-1), $$ $$ u_{2}(x_{1},x_{2})=-10x_{1}(x_{1}-1)(2x_{1}-1)x_{2}^{2}(x_{2}-1)^2. $$加罚方法的稳定化项依赖于罚参数 $\epsilon$.为使加罚法达到期望的精度,由文献 ^[39] 知罚参数需满足 $\epsilon\leq1.0e-4$,因此这里我们选取 $\epsilon=1.0e-5$. 进一步,正则性方法和多尺度增量法的稳定化项由稳定化参数 $\tau_1,\tau_2$,$\tau^{RFB}$,$\tau^{RFB}$,其中 $\tau^{RFB}$ 由式 (3.7) 决定. $\tau_1,\tau_2$ 的取值未知.局部 Gauss 积分方法其稳定化项不依赖任何与网格尺寸有关的参数.

这里我们采用一致的三角网格剖分,选用等阶 $P_1/P_1$ 有限元对,分别用加罚法、正则性法、局部 Gauss 积分法和多尺度增量方法求解 Navier-Stokes 方程.估计关于速度的 $L^2$ -范数和 $H^1$ -范数,压力的 $L^2$ -范数,相应的结果由表 1-4 展示. 理论上这些稳定化算法(参见文献 ^{[11, 19, 20, 21, 25, 27, 28, 29]}),其速度的 $L^2$ -范数和 $H^1$ -范数收敛阶分别应达到 $O(h^2)$ 和 $O(h)$,压力的 $L^2$ -范数收连阶应达到 $O(h)$. 这里我们选取 $\tau_1=1/8$ 和 $\tau_2=1/12$,因为经验得知这样选取可以得到比较好的精度. 从表 1-4 可知,除了加罚法外其余方法的速度收敛阶都达到理论分析的结果.加罚法的速度 $L^2$ -范数收敛阶随着网格的加密没有达到 $O(h^2)$.正则性方法和多尺度增量法其压力的 $L^2$ -范数渐近收敛速度近乎达到了$O(h^{2})$,局部 Gauss 积分方法其压力的 $L^2$ -范数渐近收敛速度近乎达到了 $O(h^{1.5})$,都比理论分析的预期结果好些.

另外注意到表 1-4中列出的 CPU 时间. 我们的所有计算都在同一计算平台下进行.由表 1-4 可知随着网格的加密,多尺度增量方法耗时最多,同时局部 Gauss 积分方法耗时最少.然而由表 3和表 4知,多尺度增量算法关于压力的收敛速度快些.关于此我们将在后面给出方腔驱动流算例进行更强有力的说明.

这里给出正方形区域上的方腔驱动流,该区域上具有无滑移的边界条件,仅在上边界 $\{(x,1):0<x<1\}$ 处速度满足 $u=(1,0)$,具体可见图 1.我们的目的是用这几种稳定化有限元算法 (加罚方法、正则性方法、局部 Gauss 积分法和多尺度增量法) 求解大雷诺数下的方腔驱动流,客观比较各个算法的稳定性和效率.具体分析由下面几段文字予以说明.

首先,四种稳定化有限元算法在不同雷诺数 $Re=10,6000,8000$ 下的压力等高线图和速度流线图由图 2-9 分别给出.

其次,我们对四种算法分别进行详细分析,具体如下:

由图 2,我们可以看到在雷诺数 $Re$ 取值10时,加罚法关于压力的结果失去了稳定性.同时,由于图 3中 $Re=10$ 时方腔左边底部拐角处出现了小的震荡,因此该算法对速度的计算也不好,也失去了稳定性.更糟糕的是,随着雷诺数的增大,该加罚方法完全失去了稳定性和有效性.

关于局部 Gauss 积分法,图 4和图 5分别给出了其压力等高线图和速度流线图. 从图 4可以看到,在雷诺数 $Re$ 取值10时,压力出现了震荡. 随着雷诺数的增大,该方法同样失去了稳定性和有效性.

在图 6-7中分别给出了正则性方法的压力等高线图和速度流线图.在 $Re=6000$ 时,当迭代次数 $m=439$ 时我们求得了其逼近解.进而当 $Re=8000$ 时,正则性算法失去了其有效性,结果发散.

对于多尺度增量有限元算法,图 8-9分别给出了当雷诺数 $Re$ 等于10,6000,8000时,其压力等高线图和速度流线图. 另外图 10中,我们特别给出了当 $Re=10000$ 时的压力等高线图和速度流线图.随着雷诺数的增大,第一个主漩涡向方腔的中心移动; 而第二个漩涡可能会出现在方腔的右下角底部拐角处,第三个漩涡出现在左下角底部拐角处. 当雷诺数 $Re=10000$ 时会有五个漩涡被扑捉到,这一结果与文献 ^[40] 中的结果一致. 另外,速度流线图和压力等高线图的结果也与文献 ^[41] 中的结果契合.

最后我们需要指出的是,在雷诺数 $Re=6000$ 时,采用多尺度增量有限元算法求解方腔流迭代次数 $m=89$ 时便得到了好的逼近结果.进一步,当雷诺数 $Re$ 等于8000时,由上述分析可知正则性算法已经失去了稳定性.因此,在雷诺数取值很大 $(\geq5000)$ 时方法4(多尺度增量算法)比方法2(正则性算法)相对有效些.总之,求解大雷诺数问题时方法4(多尺度增量算法)相对于其余三种算法具有好的稳定性和有效性,不失为求解 Navier-Stokes 方程的一种不错的方法.

上面的两个算例说明了多尺度增量有限算法在求解大雷诺数问题时好的稳定性和有效性.为了进一步加强我们这一结论可靠性,我们给出如下的两个经典算例再次进行说明.

这一部分我们给出经典的后台阶流问题,该问题有一个正则性的拐角.该问题的边界条件和几何特性参见图 11^[42]. 这里上边界、下边界以及入流部分具有一致的自由边界条件. 出流边界具有无外力条件.选取雷诺数 $Re=150$,最大入流速度 $u_{x\max}=1$ 及台阶高度 $l=0.5$.三角剖分选取三种不同网格尺寸,每个所含的三角单元个数分别为 234,932 和 2989.图 12-13分别给出了逼近解的速度流线图和压力等高线图.在粗网格上,多尺度增量有限元算法有一些小震荡.而当网格加密,可以看出压力仅有很轻微的震荡. 这些结论与文献 ^{[42, 43]}中结论一致.由此可知多尺度增量有限元算法对于求解后台阶流问题是一种不错的方法.

这一部分给出第四个算例,即研究圆柱绕流问题^[44].该问题的边界条件和几何特性见图 14. 上下计算边界及入流区域具有一致的自由边界条件.在出流边界施加无外力边界条件. 该圆柱绕流问题受雷诺数的影响,雷诺数这里由自由流速度和圆柱直径 $d$ 决定,即 $Re=\frac{u_{\infty}d}{\nu}$.

当选取雷诺数 $Re=0.01,10,26$ 时该问题的速度流线图和压力等高线图分别由图 15-17给出.从这三个图,可知圆柱后面出现了两个附加的涡流,且它们随着雷诺数的增大而逐渐变大. 此现象与文献 ^[42] 中的结果一致.因此可知该多尺度增量有限元算法可以较好的计算经典的圆柱绕流问题.

在这篇文章中,我们首先总结了采用 $P_1/P_1$ 有限元对求解定常 Navier-Stokes 方程的几个稳定化有限元算法.随后我们给出数值试验将几种稳定化算法进行比较. 通过对试验数据的客观分析我们得到了如下的结论:

(1) 加罚方法的稳定性和有效性依赖于罚参数的选取.当罚参数取值很小时,加罚方法的稳定性很好.但是这并不意味着罚参数可以无限取小,当罚参数非常小时系统矩阵的条件数变大以至于发散.在求解大雷诺数问题时该加罚方法失去了稳定性.

(2) 局部 Gauss积分方法不依赖于稳定性参数,且耗时最短.但是该方法在求解大雷诺数问题时不稳定,方法失效.

(3) 对正则性方法和多尺度增量法而言,两种方法皆可以求解大雷诺数问题.然后由比较可知,多尺度增量方法相比正则性方法在计算大雷诺数问题 $(\geq5000)$ 时更有效些.