基于Navier-Stokes方程分析飞机抖振问题

A Class of Spectral Method for Solving Unsteady Navier-Stokes Equations

摘要: 针对不定常的Navier-Stokes方程，给出了一种非线性Galerkin-Legendre谱方法(GL方法)，该方法是将勒让德谱方法和非线性Gelerkin方法结合起来，使得速度与压力的混合元空间不需要满足离散的Babuška-Brezzi不等式条件，给出了GL方法的谱格式，并给出该方法的稳定性证明和误差估计。关于飞机抖振问题已经由基于不可压缩N-S方程的方法的研究取代了原来粘、位流分开计算的研究，因此用谱方法研究飞机抖振问题，通过误差分析，减少阻力，达到节约能源的目的。

Abstract: This paper focuses on the issue about using numerical computation to solve unsteady Navier-Stokes equations. We present a scheme called nonlinear Galerkin-Legendre Spectral method, which combines the Legendre spectral method with nonlinear Galerkin-Legendre method. Therefore, it is not necessary to meet Babuska-Brezzi inequality condition in the velocity space and pressure space. In this paper, we give a spectral scheme of nonlinear Galerkin-Legendre methods for solving 2D N-S equations. What’s more, we also provide the error estimation and proof the stability of this scheme. The method based on the separate calculation of viscous flow and potential flow, was used to solve the problem of plane buffeting, and has been replaced by the method based on the incompressible N-S equations. Therefore, it is necessary to use the spectral method to study the chattering problem of aircraft and to reduce the resistance by error analysis and to save energy.

2. 经典Galerkin方法下的结论

2.1. 预备知识

双线性 设 $a\left(u,v\right)$ 是一个双线性泛函 [6] ，则

$a\left(u,v\right)={\displaystyle {\int }_{\Omega }{u}_x{v}_x\text{d}x}(u,v\in H_0^1(\Omega ))$

三线性设 $a\left(u,v,w\right)$ 是一个三线性泛函，则

$a\left(u,v,w\right)=\frac{1}{2}{\displaystyle {\int }_{\Omega }{\displaystyle \underset{i,j=1}{\overset{n}{\sum }}\left[u_i\frac{\partial v_j}{\partial x_i}{w}_j-u_i\frac{\partial w_j}{\partial x_i}{v}_j\right]\text{d}x}}(u,v,w\in H_0^1(\Omega ))$

对于 $a\left(u,v,w\right)$ 有下列性质(参见 [7] )：

$a\left(u,v,w\right)=-a\left(u,v,w\right);a\left(u,v,v\right)=0$

$a\left(u,v,w\right)=C{\Vert u\Vert }_{H_0^1}{\Vert v\Vert }_{H_0^1}{\Vert w\Vert }_{L^2\left(\Omega \right)};\forall u,v\in H_0^1\left(\Omega \right),w\in L^2(\Omega )$

2.2. Galerkin方法的弱形式

设 $\Omega$ 是一个有界区域。考虑下面定常N-S方程 [8] ：

问题I求 $\left(u,p\right)$ 满足

$\{\begin{array}{l}{u}_t-\nu \Delta u+u\cdot \nabla u+\nabla p=f\left(在\Omega 中\right)\\ \nabla \cdot u=0\left(在\Omega 中\right)\\ u=0\left(在\partial \Omega 上\right)\end{array}$ (1)

其中u表示速度，p表示压力，f表示外力， $\nu$ 是雷诺数的倒数，是常数。

如下记号：

$X=H_0^1\left(\Omega \right),M=L_0\left(\Omega \right)=\left\{q\in L^2\left(\Omega \right);{\displaystyle {\int }_{\Omega }q\text{d}x=0}\right\}$

问题I的变分形式为：

问题I^*求 $\left(u,p\right)\in \left(X,M\right)$ 满足

$\left(u_t,v\right)+a\left(u,v\right)+a_1\left(u,u,v\right)-b\left(p,v\right)+b\left(q,u\right)=\left(f,v\right)$

$\forall \left(v,q\right)\in \left(X,M\right)$

其中

$a\left(u,v\right)=\nu {\displaystyle {\int }_{\Omega }\nabla u\cdot \nabla v\text{d}x},b\left(q,v\right)={\displaystyle {\int }_{\Omega }qdivv\text{d}x}$

,

$\left(f,v\right)={\displaystyle {\int }_{\Omega }f\cdot v\text{d}x}$ .

设 $X_N$ 和 $M_N$ 分别是X和M空间中次数不超过N的多项式空间，则(1)的谱格式为：

$\{\begin{array}{l}\left(u_t,v_N\right)+a\left(u_N,v_N\right)+a_1\left(u_N,u_N,v_N\right)-b\left(p_N,v_N\right)=\left(f,v_N\right)\\ \left(\nabla \cdot u_N,q_N\right)=0,\forall v_N\in X_N,q_N\in M_N\end{array}$ (2)

从文献 [7] 可以得到以下结论。

定理1.1 如果Babuška不等式成立，那么

${\Vert u-u_N\Vert }_1+{\beta }_N{\Vert p-p_N\Vert }_{L^2}\le N^{1-m}{\Vert u\Vert }_m+N^{-s}{\Vert u\Vert }_s$

其中N为插值多项式的个数， ${\beta }_N$ 满足Babuška不等式。

3. 非线性Galerkin-Legendre谱方法

3.1. 准备知识

先引进一些记号。设

$S_N=span\left\{L_0\left(x\right),L_1\left(x\right),\cdots ,L_N\left(x\right)\right\}$ ,

$Q_{2N}=span\left\{L_{N-1}\left(x\right),L_N\left(x\right),\cdots ,L_{2N}\left(x\right)\right\}$ ,

$W_{2N}=\left\{w\in Q_{2N}/w\left(\pm 1\right)=0\right\}$ ,

$S_{2N}=Q_{2N}\oplus S_{N-2}=span\left\{L_0\left(x\right),L_1\left(x\right),\cdots ,L_{2N}\left(x\right)\right\}$ ,

$V_{2N}=V_N\oplus W_{2N}=\left\{v\in S_{2N}/v\left(\pm 1\right)=0\right\}$ .

注意到 $V_N$ 和中的元素均在 $x=\pm 1$ 处为零，即在边界上为零。

设 $P_N$ 是 $L^2\left(\Omega \right)\to V_N$ 的正交投影算子，即

$\left(u-P_{N}u,v\right)=0$ , $\forall v\in S_N$ , $u\in L^2(\Omega )$

设 ${\Pi }_N$ 是 $H_0^1\left(\Omega \right)\to V_N$ 的投影算子，即

$a\left(u-{\Pi }_{N}u,v\right)=0$ , $\forall v\in S_N$ , $u\in H_0^1(\Omega )$

3.2. 非线性Galerkin-Legendre方法

问题I的非线性Galerkin-Lengedre格式为 [9]

问题I^**求 $u_N\in V_N$ ，和 $u_{2N}\in W_{2N}$ ，使得

$\{\begin{array}{l}\frac{\text{d}}{\text{d}t}\left(u_N,\varphi \right)+a\left(u_N,\varphi \right)+a_1\left(u_N+u_{2N},u_N+u_{2N},\varphi \right)-b\left(p_N,\varphi \right)\\ =\left(f,\varphi \right)\text{,}\forall \varphi \in V_N\\ a\left(u_{2N},\omega \right)+a_1\left(u_N+u_{2N},u_N,\omega \right)-b\left(p_{2N},\omega \right)=\left(f,\omega \right)\\ \left(q,\nabla \cdot v\right)=0\text{,}\forall \omega \in W_{2N}\end{array}$ (3)

引理2.1根据三线性泛函的性质，有以下结论

${\Vert u_N\Vert }^2\le C\left({\Vert u\Vert }_1^2+{\Vert f\Vert }^2\right)$

以上引理在 [10] 中已经得证。

定理2.1 在引理2.1成立下，对充分大的N，我们有

$\Vert u-u_N-u_{2N}\Vert \le C{N}^{1-s}+K\Delta t^2$ , $\Vert p-p_N-p_{2N}\Vert \le M{N}^{\frac{d-1}{2}-m}$

其中C为与 ${\Vert u\Vert }_s$ 有关的常数，K为与 $\Vert u_t\Vert$ 有关的常数，M为与 ${\Vert p\Vert }_m$ 有关的常数， $d\ge 1$ ，N为插值多项式的个数。

证明：因为 ${\Pi }_{2N}$ 的定义满足

$\{\begin{array}{l}\frac{\text{d}}{\text{d}t}\left(u,\varphi +\omega \right)+\nu \left({\left({\prod }_{2N}u\right)}_x,{\left(\varphi +\omega \right)}_x\right)+a_1\left(u,u,\varphi +\omega \right)\\ =\left(f,\varphi +\omega \right),\forall \varphi \in V_N,\omega \in W_{2N}\\ b\left(q,\nabla \cdot \left(\varphi +\omega \right)\right)=0,\forall \varphi \in V_N,\omega \in W_{2N}\end{array}$

记

$e={\prod }_{2N}u-\left(u_N+u_{2N}\right)$ , $\xi ={\prod }_{N}u-u_N$ , $\eta ={\prod }_{2N}u-u_N-u{}_{2N}$ .

显然 $\xi \in V_N$ ， $\eta \in W_{2N}$ ， $e=\xi +\eta$ ，特别地

${\Vert e\Vert }_1^2={\Vert \xi \Vert }_1^2+{\Vert \eta \Vert }_1^2$

由(1)式与(3)相减

$\begin{array}{l}\frac{\text{d}}{\text{d}t}\left(u,\varphi +\omega \right)-\frac{\text{d}}{\text{d}t}\left(u_N,\varphi \right)+a\left(e,\varphi +\omega \right)+a_1\left(u,u,\varphi \right)\\ -a_1\left(u_N+u_{2N},u_N+u_{2N},\varphi \right)+a_1\left(u,u,\omega \right)-a_1\left(u_N+u_{2N},u_N,\omega \right)=0\end{array}$

整理得到

$\begin{array}{l}a\left(e,\varphi +\omega \right)+\left({\xi }_t,\varphi \right)\\ =\left[a_1\left(u_N+u_{2N},u_N+u_{2N},\varphi \right)-a_1\left(u,u,\varphi \right)\right]\\ +\left[a_1\left(u_N+u_{2N},u_N,\omega \right)-a_1\left(u,u,\omega \right)\right]+\left({\left({\Pi }_{N}u-u\right)}_t,\varphi \right)-\left(u_t,\omega \right)\end{array}$

不妨取 $\varphi =\xi$ ， $\omega =\eta$ ，于是

$\begin{array}{c}{\Vert e\Vert }_1^2+\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert \xi \Vert }^2=\left[a_1\left(u_N+u_{2N},u_N+u_{2N},\xi \right)-a_1\left(u,u,\xi \right)\right]\\ +\left[a_1\left(u_N+u_{2N},u_N+u_{2N},\eta \right)-a_1\left(u,u,\eta \right)\right]\\ +\left[\left({\left({\Pi }_{N}u-u\right)}_t,\xi \right)-\left(u_t,\eta \right)\right]\\ \triangleq I_1+I_2+I_3\end{array}$

其中

$\begin{array}{c}{I}_1=a_1\left(u_N+u_{2N},u_N+u_{2N},\xi \right)-a_1\left(u,u,\xi \right)\\ =-a_1\left(u_N+u_{2N},u-{\prod }_{2N}u,\xi \right)-a_1\left(u_N+u_{2N},e,\xi \right)-a_1\left(u-{\prod }_{2N}u+e,u,\xi \right)\\ \triangleq I_{11}+I_{12}+I_{13}\end{array}$

$\begin{array}{c}{I}_2=a_1\left(u_N+u_{2N},u_N,\eta \right)-a_1\left(u,u,\eta \right)\\ =-a_1\left(u_N+u_{2N},u-{\prod }_{N}u,\eta \right)-a_1\left(u_N+u_{2N},\xi ,\eta \right)-a_1\left(u-{\prod }_{2N}u+e,u,\eta \right)\\ \triangleq I_{21}+I_{22}+I_{23}\end{array}$

注意到 $I_{12}+I_{22}=0$ 。根据Young不等式

$\begin{array}{c}{I}_{13}+I_{23}=-a_1\left(u-{\prod }_{2N}u+e,u,e\right)\\ \le C\left({\Vert u-{\prod }_{2N}u\Vert }_1+{\Vert e\Vert }_1\right){\Vert u\Vert }_1\Vert e\Vert \\ \le \frac{\nu }{8}{\Vert e\Vert }_1^2+M\left({\Vert u-{\prod }_{2N}u\Vert }_1^2+{\Vert e\Vert }^2\right)\end{array}$

因为 ${\Vert {\prod }_N\Vert }_1\le C{\Vert u\Vert }_1$ ，应用Young不等式，Poincare不等式以及 ${\Vert \eta \Vert }_1\le {\Vert e\Vert }_1$ ，对充分大的N，有

$\begin{array}{c}{I}_{11}+I_{21}=a_1\left(e-{\prod }_{2N}u,u-{\prod }_{2N}u,e\right)+a_1\left(e-{\prod }_{N}u,{\prod }_{2N}u-{\prod }_{N}u,\eta \right)\\ \le \frac{\nu }{4}{\Vert e\Vert }_1^2+M\left({\Vert u-{\prod }_{2N}u\Vert }_1^2+{\Vert e\Vert }^2\right)+M{N}^{-2}{\Vert {\prod }_{N}u-{\prod }_{2N}u\Vert }_1^2\end{array}$

因为 $\eta \in W_{2N}$ ，

$\begin{array}{c}{I}_3=\left[\left({\left({\Pi }_{N}u-u\right)}_t,\xi \right)-\left(\left(I-P_{N-2}\right)u_t,\eta \right)\right]\\ \le \frac{\nu }{8}\left({\Vert \xi \Vert }_1^2+{\Vert \eta \Vert }_1^2\right)+C\left({\Vert \left({\Pi }_{N}u-u\right)\Vert }_{-1}^2+{\Vert \left(I-P_{N-2}\right)u_t\Vert }_{-1}^2\right)\\ \le \frac{\nu }{8}{\Vert e\Vert }_1^2+C\left({\Vert \left({\Pi }_{N}u-u\right)\Vert }_{-1}^2+{\Vert \left(I-P_{N-2}\right)u_t\Vert }_{-1}^2\right)\\ \le \frac{\nu }{8}{\Vert e\Vert }_1^2+\Delta t^2{\Vert u_t\Vert }_{-1}^2\end{array}$

综上

${\Vert e\Vert }_1^2+\frac{1}{2}\frac{\text{d}}{\text{d}t}{\Vert \xi \Vert }^2\le M{\Vert e\Vert }^2+M{N}^{-2}{\Vert {\prod }_{N}u-{\prod }_{2N}u\Vert }_1^2+\Delta t^2{\Vert u_t\Vert }_{-1}^2$ .

由Gronwall不等式

$\Vert u-u_N-u_{2N}\Vert \le N^{1-s}{\Vert u\Vert }_s+\Delta t^2{\Vert u_t\Vert }_{-1}^2$

其中 $s\ge 1$ ，N为基函数个数。

现在对压力作误差估计。

记 $V_N=\left\{v_N\in H_0^1\left(\Omega \right)/\left(q_N,\nabla \cdot v_N\right)=0,\forall q_N\in M_N\right\}$ 在文献 [11] 中有以下结果：

引理2.2 如果Babuška条件成立，则

${\beta }_N\Vert p-p_N-p_{2N}\Vert \le \underset{\begin{array}{l}{u}_N\in V_N\\ u_{2N}\in W_{2N}\end{array}}{\inf }\left|\Vert u-u_N-u_{2N}\Vert \right|+\underset{q_N\in M_N}{\inf }\Vert p-q_N\Vert$

下面我们根据Babuška不等式性质以及速度误差的相关结论有

(1)式减去(3)式

${\Vert e\Vert }_1^2-\left(p-p_N-p_{2N},\nabla \cdot v_N\right)=0$

由Babuška不等式，找任意 $q_N\in M_N$ ， $q_{2N}\in M_{2N}$ 则

$\begin{array}{l}{\beta }_N\Vert q-p_N-p_{2N}\Vert \le \underset{v_N\in H_0^1\left(\Omega \right)}{\sup }\frac{\left(q_N-p_N-p_{2N},\nabla \cdot v_N\right)}{\left|\Vert v_N\Vert \right|}\\ =\underset{v_N\in H_0^1\left(\Omega \right)}{\sup }\frac{\alpha \left(u-u_N-u_{2N},v_N\right)+\nu \left(\nabla \left(u-u_N-u_{2N}\right),\nabla v_N\right)-\left(p-q_N,\nabla \cdot v_N\right)}{\left|\Vert v_N\Vert \right|}\end{array}$

注意到 $$ ， $\forall v\in H_0^1\left(\Omega \right)$ 。

因此

$\begin{array}{c}{\beta }_N\Vert q-p_N-p_{2N}\Vert \le {\beta }_N\inf \left\{\Vert p-q_N\Vert +\Vert q_N-p_N-p_{2N}\Vert \right\}\\ \le \underset{q_N\in M_N}{\inf }\Vert p-q_N\Vert +\Vert u-u_N-u_{2N}\Vert \end{array}$

根据 ${\Vert e\Vert }_1^2\le N^{1-s}{\Vert u\Vert }_s$ ，考虑到 $P_N$ 与 $P_{N-2}$ 正交，文献 [5] 有以下结论

${\beta }_N\ge C\left(\alpha ,\nu \right)N^{-\frac{d-1}{2}}$ , $\left(d=2或3\right)$

因此

$N^{-\frac{d-1}{2}}\Vert p-p_N-p_{2N}\Vert \le N^{-m}{\Vert p\Vert }_m$

综上结论，我们有如下误差估计

${\left|u-{\tilde{u}}_N\right|}_1+{\Vert p-{\tilde{p}}_N\Vert }_{L^2}\le C{N}^{1-s}+M{N}^{\frac{d-1}{2}-m}+K\Delta t^2$

其中C是与 ${\Vert u\Vert }_s$ 有关的常数，K是与 $\Vert u_t\Vert$ 有关的常数，M是与 ${\Vert p\Vert }_m$ 有关的常数， $d\ge 2$ 。

4. 数值算例

我们考虑如下的N-S方程(飞机抖振数学模型)：

$\{\begin{array}{l}{u}_t-\nu \Delta u+u\cdot \nabla u+\nabla p=0\left(在\Omega 中\right)\\ divu=0\left(在\Omega 中\right)\\ u=0\left(在\partial \Omega 上\right)\end{array}$ (4)

其中 $\Omega =\left[-1,1\right]\times \left[-1,1\right]$ ， $0\le t\le T$ ，速度空间是一个狄利克雷边界，时间上使用差分法(显示欧拉格式)，压强空间上用次数不超过N的多项式空间近似逼近，得到如下格式：

$\left(\frac{u_N^{k+1}-u_N^k}{\Delta t},v_N\right)+\nu \left({\left(u_N^{k+1}\right)}_x,{\left(v_N\right)}_x\right)+\left(u_N^{k+1}{\left(u_N^k\right)}_x,v_N\right)+\left(\frac{p_N^{k+1}-p_N^k}{\Delta t},v_N\right)=0$ (5)

现对u做Legendre多项式插值，即

$u_N^k\approx {\displaystyle \underset{k=0}{\overset{N-2}{\sum }}{\tilde{u}}_k}{\varphi }_k\left(x\right)$ ,

对p做Legendre多项式插值，即

$p_N^k\approx {\displaystyle \underset{k=0}{\overset{N}{\sum }}{\tilde{p}}_k}{\varphi }_k$

${\tilde{u}}_k$ 、 ${\tilde{p}}_k$ 是插值系数， ${\varphi }_k$ 是插值基函数。

将这些近似多项式和代入上述(5)，则有

$\left[\begin{array}{cc}\frac{1}{\Delta t}M+K& D^{\text{T}}\\ D& 0\end{array}\right]\left[\begin{array}{c}{U}^{k+1}\\ p^{k+1}-p^k\end{array}\right]=\left[\begin{array}{c}\frac{1}{\Delta t}M{U}^k-D^{\text{T}}{p}^k\\ 0\end{array}\right]$

这是一个线性方程组形式，算法的步骤如下：

$\{\begin{array}{l}\text{Step}\text{1}:\left(\frac{1}{\Delta t}M+K\right){\tilde{U}}^{k+1}=\frac{1}{\Delta t}M{U}^k-D^{\text{T}}{p}^k\\ \text{Step}\text{2}:D{M}^{-1}{D}^{\text{T}}{\Phi }^{k+1}=\frac{1}{\Delta t}D{\tilde{U}}^{k+1}\\ \text{Step}\text{3}:U^{k+1}={\tilde{U}}^{k+1}-\Delta t{M}^{-1}{D}^{\text{T}}{\Phi }^{k+1}\\ \text{Step}\text{4}:p^{k+1}={\Phi }^{k+1}+p^k\end{array}$

已知方程(4) $\left(u,p\right)$ 的精确解

$u\left(x,y,t\right)=\text{π}\sin t\left(\sin 2\text{π}y{\sin }^2\text{π}x,-\sin 2\text{π}x{\sin }^2\text{π}y\right)$

$p\left(x,y,t\right)=\sin t\cos \text{π}x\sin \text{π}y$

处理不可压缩N-S方程的谱方法与传统的变分法以及有限元法处理N-S方程的区别很明显，从表1和表2可以看出相同的时间间隔，谱方法得出来的结果比经典格式的结果更加精确；表3中的误差表

Table 1. The comparison of error estimates for velocity in two different formats with N = 16 and T = 1

表1. N = 16和T = 1两种不同格式下的速度的误差估计的比较

Table 2. The comparison of error estimates for pressure in two different formats with N = 16 and T = 1

表2. N = 16和T = 1两种不同格式下的压强的误差估计的比较

Table 3. The error estimates and convergence order for different time intervals and N with T = 1 and T = 2

表3. T = 1和T = 2时不同时间间隔和不同N的误差估计以及收敛阶

示在T = 1、T = 2时刻下的速度误差和压力误差总和，从表3中可以看出不同时间间隔下，时间间隔对方程的误差估计影响更大，该方程的收敛阶可以用 $\Delta t^2$ 来刻画，并且速度收敛阶是1；在相同时间间隔和不同N下，N对方程的误差估计影响更大，该方程的收敛阶可以用N来刻画，且随着基函数的个数N的增加，误差越来越少，所得的精度越高，收敛结果也符合谱方法的收敛阶，以N的负指数次幂速度收敛精确解。

5. 结语

研究提出了一类求解定常纳维斯托克斯方程的谱方法，基于二维不可压缩N-S方程给出了LG格式，并且得出相关误差结果，比较了该方法与经典方法和有限元方法的区别，从数值结果我们可以发现LG谱方法下速度逼近的收敛率是最优的，而压强逼近的收敛率也很好。从表格上可以看出当解 $\left(u,p\right)$ 是光滑的，即使对较小的N也能得到很好的结果，而且计算量不大，这充分体现了LG谱方法的优越性。基于N-S方程的GL方法的数值分析，可以看出对于飞机抖振问题，可以通过减少与空气阻力，达到节约能源的目的。

基金项目

贵州省科学技术基金项目(黔科合J字[2013]2128)，贵州省大数据重点实验室开放课题 (2017BDKFJJ012)。

NOTES

^*通讯作者。