A Subgrid Stabilizing Postprocessed Mixed Finite Element Method for the Navier-Stokes Equations

定义1  设Ω是${{\mathbb{R}}^{d}}$(d = 2，3)中具有齐次Dirichlet边界条件的有界区域，我们考虑下面的Navier-Stokes方程：

其中：u：Ω→${{\mathbb{R}}^{2}}$表示速度矢量，p：Ω→$\mathbb{R}$表示压力，f：Ω→${{\mathbb{R}}^{2}}$表示流体驱动的体积力，ν＞0为流体粘性系数，u₀是使得▽·u₀ = 0的初始速度，并且u_t = ∂u/∂t，定常数T表示最后时刻.

下面针对定义1提出的数学问题，我们引进两个Hilbert空间：

对H，V分别赋予L²(Ω)^d空间和H₀¹(Ω)^d空间的内积.用‖·‖_l表示Sobolev空间W^l，q(Ω)^d的范数，‖·‖_-l表示其对偶空间的范数.后面还会涉及到商空间H^l(Ω)/$\mathbb{R}$，其范数定义为：

设$\mathscr{T}$_h = (τ_i^h，ϕ_i^h)_i∈Ih(h＞0)是区域Ω_h的准均匀网格剖分，其中h是单元τ_i^h∈${{\mathscr{T}}_{\mathit{h}}}$的最大直径，ϕ_i^h：τ₀→τ_i^h为单射.当r≥2时，基于网格剖分${{\mathscr{T}}_{\mathit{h}}}$，我们定义下列有限元空间：

其中P_r-1(τ₀)表示定义在τ₀上的至多r-1次多项式空间.

接下来，我们考虑Hood-Taylor元(X_h，r，Q_h，r-1)，其中

并且这个有限元对(X_h，r，Q_h，r-1)满足inf-sup条件^[6].

根据非线性项b(·，·，·)的定义^[9]，我们具体给出离散的非线性项b_h(·，·，·)为：

亚格子模型是基于一个椭圆型投影算子${{\mathit{\Pi }}_{h}}:{{H}_{0}}^{1}{{\left(\mathit{\Omega } \right)}^{d}}\to {{R}_{1}} = \left\{ v\in {{H}_{0}}^{1}{{\left(\mathit{\Omega } \right)}^{d}}:v\left| _{\pi _{i}^{h}} \right.\in {{({{p}_{1}})}^{d}}, \forall {{\tau }_{i}}{{~}^{h}}\in {{T}_{h}} \right\}$所定义，其中算子Π_h的定义为：

下面根据投影算子Π_h，定义亚格子稳定项为：

其中α＞0是稳定化参数，通常α = o(h^s)(s是与所使用有限元空间的阶数有关的参数).

下面在T时刻逼近问题(1)-(4).基于3种后处理的有限元方法^[6-8]和亚格子稳定化方法^[9]，我们的亚格子稳定化的后处理有限元方法具体如下：

算法1  基于Stokes问题的稳定化后处理有限元方法.

1) 给定u₀^h逼近u(0)，寻找u_h：[0，T]→X_h，r和Q_h，r-1使得：

其中${{{\mathit{\boldsymbol{\dot{u}}}}}_{h}} = \partial {{\mathit{\boldsymbol{u}}}_{h}}/\partial t$，并且通常取u₀^h = S_h(u₀)，S_h是离散Stokes投影^[6].

2) 对离散的速度和压力(u_h(T)，p_h(T))后处理，寻找$({{{\mathit{\boldsymbol{\tilde{u}}}}}_{{\tilde{h}}}}, \text{ }{{{\mathit{\boldsymbol{\tilde{p}}}}}_{{\tilde{h}}}})\in \left(\tilde{X}, \tilde{Q} \right)$使得：

其中后处理的有限元$\left(\tilde{X}, \tilde{Q} \right)$有以下两种取法：

方法1  更细网格上的同阶Hood-Taylor元$\left(\tilde{X}, \tilde{Q} \right) = ({{X}_{\tilde{h}, r}}, {{Q}_{\tilde{h}, \text{ }r-1}}), \ r\ge 3, \ \tilde{h} < h$.

方法2 相同网格上的高阶Hood-Taylor元$\left(\tilde{X}, \tilde{Q} \right) = ({{X}_{\tilde{h}, r+1}}, {{Q}_{\tilde{h}, \text{ }r}}), \ r\ge 3, \ \tilde{h} = h$.

算法2 基于Newton型问题的稳定化后处理有限元方法.

1) 给定u^h(0)逼近u(0)，寻找u_h：[0，T]→X_h，r和p_h：[0，T]→Q_h，r-1，使得(u_h，p_h)满足(5)-(6)式；

2) 对离散的速度和压力(u_h(T)，p_h(T))后处理，寻找$({{{\mathit{\boldsymbol{\tilde{u}}}}}_{{\tilde{h}}}}, \text{ }{{{\mathit{\boldsymbol{\tilde{p}}}}}_{{\tilde{h}}}})\in \left(\tilde{X}, \tilde{Q} \right)$使得：

算法3  基于Ossen型问题的稳定化后处理有限元方法.

1) 给定u^h(0)逼近u(0)，寻找u_h：[0，T]→X_h，r和p_h：[0，T]→Q_h，r-1，使得(u_h，p_h)满足(5)-(6)式；

2) 对离散的速度和压力(u_h(T)，p_h(T))后处理，寻找$({{{\mathit{\boldsymbol{\tilde{u}}}}}_{{\tilde{h}}}}, \text{ }{{{\mathit{\boldsymbol{\tilde{p}}}}}_{{\tilde{h}}}})\in \left(\tilde{X}, \tilde{Q} \right)$使得：

本节将利用FreeFem++软件^[10]进行数值实验来证明算法的有效性和比较3种亚格子稳定化后处理混合有限元方法的不同.在实验中，Hood-Taylor元用于空间离散，算法中第一步的非线性系统由Crank-Nicolson Admas-Bashforth方法求解，第二步在相同的网格上用高阶Hood-Taylor元P₃-P₂元后处理.

特别地，对于算法1，我们选择Navier-Stokes方程的准确解为：

其中，求解区域Ω = [0, 1]×[0, 1] ⊂${{\mathbb{R}}^{2}}$，并且取ν = 1，T = 0.1，α = 0.1 h²|log(h)|. 表 1给出了T时刻稳定化的混合有限元逼近解和后处理所得解的第一个分量的L²-范误差，其中‖u(T)-u_N^h‖₀表示稳定化的混合有限元速度解的误差，${{\left\| \mathit{\boldsymbol{u}}\left(T \right)-{{{\mathit{\boldsymbol{\tilde{u}}}}}_{N}}^{h} \right\|}_{0}}$表示稳定化后处理速度解的误差，第二个分量类似可得；表 2为H¹-范误差；表 3为压力解的误差比较，其中‖P(T)-P_h‖₀为稳定化的混合有限元解的压力误差，${{\left\| \mathit{\boldsymbol{P}}\left(T \right)-{{{\mathit{\boldsymbol{\tilde{P}}}}}_{h}} \right\|}_{0}}$表示稳定化后处理压力误差.通过表 1-3可知，稳定化的后处理混合有限元方法确实提高了稳定化的混合有限元解的精确度.

为了对比3种算法，取网格尺寸h = 1/160，时间步长Δt = 1/8 000，分别取ν = 1，0.1，0.01和0.001，进行数值求解.数值结果在表 4-6中给出，通过比较发现：从花费时间的长短来看，除ν = 1以外，在其它情况下稳定化的Newton型后处理花费的时间相对较多，而稳定化的Ossen型后处理花费的时间相对较少.从精确度来看，用3种算法所得速度的L²-范误差基本一致，但由算法2和算法3所得的速度的H¹-范误差和压力的L²-范误差比算法1更准确.

在这个实验中，我们考虑一个定义在Ω = [0, 6]×[0, 1]上的圆柱绕流问题，其中该圆柱的圆心(x，y) = (1，0.5)，半径为0.15.取入流速度为

并且在其它边界上满足无滑动边界条件^[11].在这个问题中，Hood-Taylor元用于空间离散，粘性系数ν = 0.001.设置网格尺寸为h = 1/12，稳定化参数α = 0.1h²|log(h)|，时间步长Δt = 0.001，分别取最后时刻T = 1，2，4，6，8，10. 图 1描述了3种算法关于动能和时间T的关系.从图中可以看出在T = 1，2，4，6时3种算法所得解的动能差距不大，但在T = 8，10时算法2和算法3所得解的动能比算法1的大.

图 2描绘了在T = 1，2，4，6，8，10的情况下，算法2在圆柱周围的速度分量${{{\mathit{\boldsymbol{\tilde{u}}}}}_{\tilde{h}1}}$的等值线；图 3是速度分量${{{\mathit{\boldsymbol{\tilde{u}}}}}_{\tilde{h}2}}$的等值线；而图 4是压力等值线.对于算法1和算法3，可得到类似的结果.

在这篇文章中，我们主要通过数值实验研究了3种基于亚格子模型的后处理混合有限元方法，实验结果表明：在选取适当的稳定化参数和网格尺寸的条件下，3种稳定化的后处理有限元方法提高了稳定化的混合有限元解的精确度，并且收敛阶也明显提高了一阶.但是从计算时间看，除ν = 1以外，在其它情况下稳定化的Newton型后处理花费的时间相对较多，而稳定化的Ossen型后处理花费的时间相对较少.从精确度来看，用3种算法所得的速度的L²-范误差基本一致，但Newton型后处理和Ossen型后处理方法所得的速度的H¹-范误差和压力的L²-范误差比Stokes型后处理方法更有效.