微极Navier-Stokes方程的一种二阶时间步算法

A Second-Order Time-Step Algorithm for the Micropolar Navier-Stokes Equations

摘要: 针对微极Navier-Stokes方程(MNSE)，我们提出了一种新的二阶时间步算法。对MNSE中的非线性项进行了线性化处理，并在线速度、压力和角速度的离散解中加入了“曲率稳定”项，旨在改进常用的“速率稳定”。该方法不仅克服了阻力产生的数值不稳定，并且在不增加计算复杂性的情况下将解的精度从一阶提高到二阶。然后我们给出了算法的无条件稳定性。最后，通过数值实验验证了预测的收敛速度。

Abstract: We propose a new second-order time-step algorithm for the micropolar Navier-Stokes equations (MNSE). The nonlinear term in MNSE is linearized, and the term “curvature stabilization” is added to the discrete solutions of online velocity, pressure and angular velocity, which aims to improve the commonly used “rate stabilization”. This method not only overcomes the numerical instability caused by resistance, but also improves the accuracy of the solution from the first order to the second order without increasing the computational complexity. And then we give the unconditional stability of the algorithm. Finally, the predicted convergence rate is verified by numerical experiments.

1. 引言

微极Navier-Stokes方程(MNSE)是一种重要的物理模型，该方程一般由Navier-Stokes方程耦合角动量方程而成，可以描述具有微观结构的流体运动现象，如材料质点具有平移和旋转自由度时的不可压缩流体的运动等。针对具有微观结构的流体模型的新型高性能的数值算法研究具有重要的理论和现实价值。1964年，Eringen提出了一种通用的微极流体理论来描述这类流体的运动，同时考虑了内部子结构的特性。文 [1] 证明了磁微极流体运动强解的存在性和唯一性。事实上，大量的研究人员使用数值方法和技术来讨论这类问题。例如：Nochetto等人 [2] 提出了线速度和角速度解耦的MNSE一阶半隐式全离散有限元方法；Salgado研究了微极Navier-Stokes方程(MNSE)的全离散分数阶时间步长投影有限元法 [3] ；文 [4] 针对Navier-Stokes方程提出了一种新的、二阶精度的数值正则化方法；文 [5] 提出并分析了微极Navier-Stokes方程(MNSE)的两种全离散投影数值格式；文 [6] 将SAV方法与压力投影方法相结合，提出了求解微极Navier-Stokes方程的无条件能量稳定数值格式；文 [7] 提出了不可压缩微极性流体流动控制方程的二阶Gauge-Uzawa格式；文 [8] 针对二维/三维非定常不可压缩热–微极流体方程，提出了三种分裂有限元格式，分别为标准分裂格式、一致分裂格式。文 [9] 将对流项的标量辅助变量(SAV)方法和耦合项的隐–显(IMEX)处理方法相结合，提出了求解微极Navier-Stokes方程的解耦、线性和无条件能量稳定数值格式。主要研究内容是针对微极Navier-Stokes方程提出一种新的、二阶无条件稳定的时间步方法，给出该格式的无条件稳定性。最后用数值实验验证该格式的收敛速度。

2. 符号和预备知识

微极Navier-Stokes方程可以描述如下：对于有界正则域 $\Omega \subset {ℝ}^d\left(d=2or3\right)$ ，给定最终时间 $T$ ，找到 $u:\Omega \times \left(0,d\right]\to R^d，p:\Omega \times \left(0,d\right]\to R$ 和 $w:\Omega \times \left(0,d\right]\to R^d$ 满足：

$\begin{array}{l}{u}_t-\left(s+s_r\right)\Delta u+\left(u\cdot \nabla \right)u+\nabla p=2s_r\nabla \times w+f\\ \nabla \cdot u=0\\ j{w}_t-\left(c_a+c_d\right)\Delta w+j\left(u\cdot \nabla \right)w-\left(c_0+c_d-c_a\right)\nabla \nabla \cdot w+4s_{r}w=2s_r\nabla \times u+g\end{array}$ (1)

边界条件为： $u\left(x,t\right)=0,w\left(x,t\right)=0,\forall \left(x,t\right)\in \partial \Omega \times \left(0,T\right]$ ，

初始条件为： $u\left(x,0\right)=u_0,w\left(x,0\right)=w_0,\forall x\in \Omega .$ 。

这里的 $u,p$ 和 $w$ 分别是线速度、压力和角速度。 $f$ 和 $g$ 的函数分别表示平滑的外部施加的力和力矩。我们假设所有材料常数 $j,s,s_r,c_0,c_a$ 和 $c_d$ 的运动粘度均为常数且为正数。我们设 $s_0=s+s_r,c_1=c_a+c_d,c_2=c_0+c_a-c_d$ ，并且 $c_1,c_2>0$ 。特别的，当 $d=2$ 时，我们令 $u=\left(u_1,u_2,0\right),w=\left(0,0,w\right)$ ，那么有

$\nabla \times u=\left(0，0，\frac{\partial u_2}{\partial x}-\frac{\partial u_1}{\partial y}\right),\nabla \times w=\left(\frac{\partial w}{\partial y},-\frac{\partial w}{\partial x},0\right)$

我们令 $L^p\left(\Omega \right)$ 空间中的范数表示为 ${\Vert \cdot \Vert }_{L^p}$ ， $W_k^p\left(\Omega \right)$ 空间中的范数表示为 ${\Vert \cdot \Vert }_{W_k^p}$ 。特别的， $L^2{\left(\Omega \right)}^d$ 空间中的范数表示为 $\Vert \cdot \Vert$ ， $L^2{\left(\Omega \right)}^d$ 空间中的范数表示为 ${\Vert \cdot \Vert }_{\infty }$ 。接下来，我们定义一些Sobolev空间：

$\begin{array}{l}X=H_0^1{\left(\Omega \right)}^d=\left\{v\in H^1{\left(\Omega \right)}^d:v|{}_{\partial \Omega }=0\right\}\\ Q=L_0^2\left(\Omega \right)=\left\{q\in L^2\left(\Omega \right):{\displaystyle {\int }_{\Omega }q}dx=0\right\}\\ Y=\left\{z\in H^1{\left(\Omega \right)}^d:z\cdot n|{}_{\partial \Omega }=0\right\}\end{array}$

此外，无散函数的空间由下式给出：

$V=\left\{v\in X:\left(\nabla \cdot v,q\right)=0,\forall q\in Q\right\}$

对于在整个时间间隔 $\left(\text{0,T}\right)$ 上定义的函数 $u\left(x,t\right)$ ，我们定义以下范数：

${\Vert u\Vert }_{L^{\infty }\left(0,T;X\right)}=\text{ess}\underset{0<t<T}{\text{sup}}{\Vert u\left(t\right)\Vert }_X,{\Vert u\Vert }_{L^p\left(0,T;X\right)}={\left({\displaystyle {\int }_0^T{\Vert u\left(t\right)\Vert }_X^p}dt\right)}^{\frac{1}{p}},1\le p<\infty$

令 $V^{-1}$ 表示 $V$ 的对偶空间， $V^{-1}$ 的范数由下式给出

${\Vert f\Vert }_{-1}=\underset{0\ne v\in V}{\sup }\frac{\left(f,v\right)}{\Vert \nabla v\Vert }$

我们定义离散无散空间：

$V_h=\left\{v_h\in X_h:\left(\nabla \cdot v_h,q_h\right)=0,\forall q_h\in Q_h\right\}$

通常的显式偏斜对称三线性形式定义为：

$b^{*}\left(u,v,w\right)=\frac{1}{2}\left(u\cdot \nabla v,w\right)-\frac{1}{2}\left(u\cdot \nabla w,v\right)$

我们可以将(1)的弱形式表示为：找到 $u:\left(0,T\right]\to X,p:\left(0,T\right]\to Q$ 和 $w:\left(0,T\right]\to X$ ，在时间 $t\in \left(0,T\right]$ 上满足：

$\begin{array}{l}\left(u_t,v\right)+s_0\left(\nabla u,\nabla v\right)+b^{*}\left(u,u,v\right)-\left(p,\nabla \cdot v\right)=2s_r\left(\nabla \times w,v\right)+\left(f,v\right)\\ \left(\nabla \cdot u,q\right)=0\\ j\left(w_t,z\right)+c_1\left(\nabla w,\nabla z\right)+j{b}^{*}\left(u,w,z\right)+c_2\left(\nabla \cdot w,\nabla \cdot z\right)+4s_r\left(w,z\right)\\ =2s_r\left(\nabla \times u,z\right)+\left(g,z\right)\end{array}$

3. 算法及其稳定性

提出了一种二阶全离散时间步进算法，该算法将非线性项的线性处理与离散弱解在时间上的稳定项相结合，并给出了该方法的无条件稳定性。

3.1. 二阶时间步算法

我们定义最终时间 $T$ 、时间步数 $N$ 、时间步 $\Delta t=\frac{T}{N}$ 和 $t_n=n\Delta t,n=1,\mathrm{2...}N$ 。 $t=t_{n+1}$ 处的插值由下式给出：

${\mathcal{J}}_{n+1}^{ϵ}\left(a\right):=\frac{\nu +ϵ}{\nu }{a}^{n+1}+\left(1-\frac{\nu +2ϵ}{\nu }\right)a^n+\frac{ϵ}{\nu }{a}^{n-1}$

$t=t_{n+1}$ 处的外推由下式给出：

${\mathcal{H}}_{n+1}\left(a\right):=2a^n-a^{n-1}$

接下来我们给出二阶全离散时间步算法。

算法3.1 给定 $u_h^{n-1},u_h^n,w_h^{n-1},w_h^n$ 和 $p_h^{n-1},p_h^n$ ，找到 $\left(u_h^{n+1},p_h^{n+1},w_h^{n+1}\right)\in \left(X_h,Q_h,X_h\right),n=0,1,\dots ,N-1$ 满足：

$\begin{array}{l}\left(\frac{3u_h^{n+1}-4u_h^n+u_h^{n-1}}{2\Delta t},v_h\right)+s_0\left({\mathcal{J}}_{n+1}^{ϵ}\left(\nabla u_h\right),\nabla v_h\right)\\ +b^{\ast }\left({\mathcal{H}}_{n+1}\left(u_h\right),{\mathcal{J}}_{n+1}^{ϵ}\left(u_h\right),v_h\right)-\left({\mathcal{J}}_{n+1}^{ϵ}\left(p_h\right),\nabla \cdot v_h\right)\\ =2s_r\left({\mathcal{J}}_{n+1}^{ϵ}\left(\nabla \times w_h\right),v_h\right)+\left(f^{n+1},v_h\right)\\ \left(\nabla \cdot u^{n+1},q_h\right)=0\\ j\left(\frac{3w_h^{n+1}-4w_h^n+w_h^{n-1}}{2\Delta t},z_h\right)+j{b}^{\ast }\left({\mathcal{H}}_{n+1}\left(u_h\right),{\mathcal{J}}_{n+1}^{ϵ}\left(w_h\right),z_h\right)\\ +c_1\left({\mathcal{J}}_{n+1}^{ϵ}\left(\nabla w_h\right),\nabla z_h\right)+c_2\left({\mathcal{J}}_{n+1}^{ϵ}\left(\nabla \cdot w_h\right),\nabla \cdot z_h\right)+4s_r\left({\mathcal{J}}_{n+1}^{ϵ}\left(w_h\right),z_h\right)\\ =2s_r\left({\mathcal{J}}_{n+1}^{ϵ}\left(\nabla \times u_h\right),z_h\right)+\left(g^{n+1},z_h\right)\end{array}$ (2)

3.2. 稳定性分析

为了分析稳定性，我们通过 $F\in {ℝ}^{n\times n}$ 引入对称矩阵 $F=I+\frac{4ϵ}{\nu }I$ ，以及对称矩阵 $G\in {ℝ}^{2n\times 2n}$ 如下：

$G=\left(\begin{array}{cc}\frac{5\nu +2ϵ}{4\nu }I& -\frac{\nu +ϵ}{2\nu }I\\ -\frac{\nu +ϵ}{2\nu }I& \frac{\nu +2ϵ}{4\nu }I\end{array}\right)$

我们定义 $F$ 范数和 $G$ 分别为：

${\Vert u\Vert }_F=\left(u,Fu\right),{\Vert \left[\begin{array}{c}u\\ w\end{array}\right]\Vert }_G^2=\left(\left[\begin{array}{c}u\\ w\end{array}\right],G\left[\begin{array}{c}u\\ w\end{array}\right]\right),\forall u,w\in {ℝ}^{n\times n}$

定理3.2算法3.1无条件稳定并且满足以下界限：

$\begin{array}{l}{\Vert u_h^N\Vert }^2+j{\Vert w_h^N\Vert }^2+\frac{2\Delta t}{3}{\displaystyle \underset{n=1}{\overset{N-1}{\sum }}\left(s{\Vert {\mathcal{J}}_{n+1}^{ϵ}\left(\nabla u_h\right)\Vert }^2+c_1{\Vert {\mathcal{J}}_{n+1}^{ϵ}\left(\nabla w_h\right)\Vert }^2+2c_2{\Vert {\mathcal{J}}_{n+1}^{ϵ}\left(\nabla \cdot w_h\right)\Vert }^2\right)}\\ +\frac{1}{3}{\displaystyle \underset{n=1}{\overset{N-1}{\sum }}\left({\Vert u_h^{n+1}-2u_h^n+u_h^{n-1}\Vert }_F^2+j{\Vert w_h^{n+1}-2w_h^n+w_h^{n-1}\Vert }_F^2\right)}\\ \le \frac{1}{3^N}\left({\Vert u_h^0\Vert }^2+j{\Vert w_h^0\Vert }^2\right)+\frac{4N}{3}\left({\Vert \left[\begin{array}{c}{u}_h^1\\ u_h^0\end{array}\right]\Vert }_G^2+j{\Vert \left[\begin{array}{c}{w}_h^1\\ w_h^0\end{array}\right]\Vert }_G^2\right)\\ +\frac{2N\Delta t}{3}{\displaystyle \underset{n=1}{\overset{N-1}{\sum }}\left(s^{-1}{\Vert f^{n+\theta }\Vert }_{-1}^2+c_1^{-1}{\Vert g^{n+\theta }\Vert }_{-1}^2\right)}\end{array}$

4. 数值实验

在这一部分中，我们使用了两个数值实验来验证所提出的方法的理论结果。我们用BDF2Reg测试了算法3.1的收敛速度。我们使用FreeFem++完成以下实验。

在这个实验中，我们的目的是验证理论结果的收敛速度。我们考虑域 $\Omega ={\left[0,1\right]}^2$ ，设置其他参数 $j=s=s_r=c_0=c_a=c_d=1$ 。我们假设已知的函数

$\begin{array}{l}u\left(x,y,t\right)=\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}10\left(x^4-2x^3+x^2\right)\left(2y^3-3y^2+y\right)\cos t\\ -10\left(2x^3-3x^2+x\right)\left(y^4-2y^3+y^2\right)\cos t\end{array}\right)\\ p\left(x,y,t\right)=10\left(2x-1\right)\left(2y-1\right)\cos t\\ w\left(x,y,t\right)=u_1+u_2\end{array}$

作为真解。

对于线速度、压力和角速度，我们使用混合有限元 $P_2-P_1-P_2$ 。误差估计定义如下：

${\Vert v-v_h\Vert }_{2,1}={\left\{\Delta t{\displaystyle \underset{n=1}{\overset{N}{\sum }}{\Vert v\left(t^n\right)-v_h^n\Vert }^2}\right\}}^{1/2},\left(v=u,w\right)$

我们通过改变网格宽度 $h$ 和时间步长 $\Delta t$ 来测试二阶全离 $T=0.1,\Delta t=T/16,ϵ=0.1$ 散时间步进算法3.1的精度。对于表1中的空间收敛阶数，我们取最终时间和不同的网格宽度。对于表2中的时间收敛顺序，我们采用最终时间 $T=3,h=T/128,ϵ=0.1$ 和变化的时间步长。我们可以看到，空间收敛和时间收敛的顺序都为2。因此，该实验验证了理论的正确性。

Table 1. Order of spatial convergence

表1. 空间收敛阶数

Table 2. Order of time convergence

表2. 时间收敛阶数

5. 结语

本文提出了微极Navier-Stokes方程的一种新的二阶全离散时间步算法，我们对非线性项进行线性化处理，在时间上增加离散弱解的稳定项，在空间上增加协调有限元。然后给出了该算法的无条件稳定性，并通过数值试验验证了所提出算法的理论结果。本文提出的算法将解的精度从一阶提高到了二阶，并且可以轻松应用于许多其他模型，是研究这些线性化方法和CFD背景下新的稳定性的初步尝试。稳定项如何单独影响稳定性仍需要更多研究。一般来说，需要更广泛的数值实验，特别是对于存在间歇性和大时间变化的问题。