﻿ 二维不可压缩Navier-Stokes方程的并行谱有限元法求解
 计算机应用   2017, Vol. 37 Issue (1): 42-47  DOI: 10.11772/j.issn.1001-9081.2017.01.0042 0

### 引用本文

HU Yuanyuan, XIE Jiang, ZHANG Wu. Solution of two dimensional incompressible Navier-Stokes equation by parallel spectral finite element method[J]. JOURNAL OF COMPUTER APPLICATIONS, 2017, 37(1): 42-47. DOI: 10.11772/j.issn.1001-9081.2017.01.0042.

### 文章历史

1. 上海大学 计算机工程与科学学院, 上海 200444 ;
2. 上海大学 高性能计算中心, 上海 200444

Solution of two dimensional incompressible Navier-Stokes equation by parallel spectral finite element method
HU Yuanyuan1, XIE Jiang1, ZHANG Wu2
1. School of Computer Engineering and Science, Shanghai University, Shanghai 200444, China ;
2. High Performance Computing Center, Shanghai University, Shanghai 200444, China
Abstract: Due to a large number of computational grids and slow convergence existed in the numerical simulation of Navier-Stokes (N-S) equation, Triangular mesh Spectral Finite Element Method based on area coordinate (TSFEM) was proposed. And further, TSFEM was paralleled with OpenMP. Spectral method was combined with finite element method, and the exponential function with infinite smoothness was selected as the basis function to replace the polynomial function in the traditional finite element method, which can efficiently reduce the amount of computational grids as well as improve the convergence and accuracy of the proposed algorithm. Because area coordinates can facilitate the calculation of triangular units, which were selected as the computing units to enhance the applicability of the algorithm. The lid-driven cavity flow was used to verify the TSFEM. The experimental results show that, compared with the traditional Finite Element Method (FEM), the TSFEM greatly improves the convergence rate and the calculation efficiency.
Key words: incompressible Navier-Stokes (N-S) equation    OpenMP    lid-driven cavity flow    high-precision    infinite convergence
0 引言

Navier-Stokes (N-S)方程是流体力学中最重要的方程之一。数学家和物理学家深信，无论是微风还是湍流，都可以通过求解N-S方程来进行解释和预言，因此研究N-S方程具有广泛的应用价值。对N-S方程的研究距今已有200多年的历史，其弱解又称为Leray-Hopf弱解。关于N-S方程强解的局部适定性、存在性与光滑性被列为21世纪7个价值100万美元的数学难题之一。数学家断言，如果没有新的分析工具和数学思想，这个难题将很难得到解决。但是,到目前为止,证明弱解的唯一性和正则性,即强解的整体存在性,仍是一个极具挑战性的问题。只有极少数非常简单的流动问题才能求得其精确解，大多数还是要用离散的方法求得数值解。在利用数值方法对其进行求解时，计算格式的时间步长受制于CFL(Courant-Friedrichs-Lewy)数，CFL=λ·Δτ/Δη，λ是与当地声速成正比的特征函数。当流体近似不可压缩时，λ趋向于无穷大，而CFL数是一个定值，因此最大时间步长Δτ趋向于0。如果用可压缩N-S方程求解近似不可压缩流场，计算效率将会变得极低，因此只能通过不可压缩N-S方程来计算近似不可压流场。不可压缩N-S方程的主要优点是求解未知数的数量较少，将连续性方程和动量方程联立，求解得到速度和压力后，代入能量方程就可以得到温度，计算效率较高。

1 不可压缩N-S方程的谱解法

 $\left\{ \begin{matrix} -v{{\nabla }^{2}}u+(u\cdot \nabla )u+\nabla p+f=\frac{\partial u}{\partial t}, & x\in \Omega ,t\in (0,T) \\ divu=0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ & x\in \Omega ,t\in (0,T) \\ u(x,0)=\bar{u}(x)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ & x\in \Omega \\ \end{matrix} \right.$ (1)

 $\int_{\Omega }{p(x,0)}dx=0;t\in [0,T]$ (2)

 ${{u}^{N}}(x,t)=\sum\limits_{k=-N}^{N-1}{{{{\hat{u}}}_{k}}{{\Phi }_{k}}(x)}$ (3)

 \begin{align} & \int_{0}^{2\pi }{{{\Phi }_{k}}(x){{{\bar{\Phi }}}_{l}}}(x)dx=2\pi {{\delta }_{kl}}= \\ & \left\{ \begin{matrix} 0, & k=l, \\ 2\pi , & k\ne l \\ \end{matrix}k,l=0,1,2,... \right. \\ \end{align} (4)

 $\int_{0}^{2\pi }{(\frac{\partial N}{\partial t}-L{{u}^{N}}){{{\bar{\Phi }}}_{k}}(x)dx=0,k\in [-N,N-1]}$ (5)

 \begin{align} & \sum\limits_{s=-N}^{N-1}{\int\limits_{0}^{2\pi }{{{\Phi }_{s}}(x){{{\bar{\Phi }}}_{k}}(x)dx}\frac{d{{{\hat{u}}}_{k}}}{dt}}-\int\limits_{0}^{2\pi }{L}{{u}^{N}}\cdot {{{\bar{\Phi }}}_{k}}(x)dx=0, \\ & k\in [-N,N-1] \\ \end{align}

 $2\pi \frac{d{{{\hat{u}}}_{k}}}{dt}-\int\limits_{0}^{2\pi }{L{{u}^{N}}\cdot {{{\bar{\Phi }}}_{k}}(x)dx=0,k\in [-N,N-1]}$

 $2\pi \frac{d{{{\hat{u}}}_{k}}}{dt}-(kv{{\hat{u}}_{k}}+ik{{\hat{u}}_{k}}+\nabla p+f)=0$

2 TFSEM

2.1 三角形单元谱有限元法

 图 1 面积坐标图 Figure 1 Area coordinate system

P点的坐标p(li,lj,lm)可由三个比值来唯一确定:

lii

ljj

lmm

Δijm

u=(u1,u2)T和p满足：

 $\left\{ \begin{matrix} \text{-}v\Delta u+\sum\limits_{i=1}^{2}{{{u}_{i}}\frac{\partial {{u}_{i}}}{\partial x}+\nabla p=f,} & 在\Omega 内, \\ div\ u=0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ & 在\Omega 内, \\ u=0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ & 在\partial \Omega 内. \\ \end{matrix} \right.$ (6)

2.2 混合变分

 $M=\left\{ q:q\in {{L}^{2}}(\Omega )and\int_{\Omega }{q}dx=0 \right\}$

V≡(H01(Ω))2

 \begin{align} & a(u,v)=v\int_{\Omega }{\nabla u\cdot \nabla vd\Omega } \\ & {{a}_{1}}(w;u,v)=\frac{1}{2}\sum\limits_{i,j=1}^{n}{\int_{\Omega }{{{w}_{j}}(\frac{\partial {{u}_{i}}}{\partial {{x}_{j}}}{{v}_{i}}-}}\frac{\partial {{u}_{i}}}{\partial {{x}_{j}}}{{u}_{i}})d\Omega \\ & b(v,q)=-\int_{\Omega }{qdiv\ vd\Omega } \\ & (f,v)\int_{\Omega }{f\cdot vd\Omega } \\ \end{align}

 $\left\{ \begin{matrix} a(u,v)+{{a}_{1}}(w;u,v)+b(v,p)=(f,v), & \forall v\in V, \\ b(u,q)=0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ & \forall q\in M \\ \end{matrix} \right.$

 \left\{ \begin{align} & {{\xi }_{j}}=-\cos \left( {{\pi }_{j}}/N_{x}^{i} \right) \\ & \ \eta =-\cos \left( {{\pi }_{k}}/N_{y}^{i} \right) \\ \end{align} \right.
2.3 选取插值基函数

 \begin{align} & {{u}^{i}}(\xi ,\eta )=\sum\limits_{j=0}^{N_{x}^{i}}{\sum\limits_{k=0}^{N_{y}^{i}}{u_{jk}^{i}\Phi _{j}^{i}(\xi )}}\Phi _{k}^{i}(\eta ) \\ & {{v}^{i}}(\xi ,\eta )=\sum\limits_{p=0}^{N_{x}^{i}}{\sum\limits_{q=0}^{N_{y}^{i}}{v_{pq}^{i}\Phi _{p}^{i}(\xi )}}\Phi _{q}^{i}(\eta ) \\ \end{align}

 $a({{u}^{i}},{{v}^{i}})+{{a}_{1}}(w;{{u}^{i}},{{v}^{i}})+b(v,p)=(f,{{v}^{i}})$

 $\sum\limits_{j=0}^{N_{x}^{i}}{\sum\limits_{k=0}^{N_{y}^{i}}{u_{jk}^{i}(A_{jkpq}^{i}+a\cdot B_{jkpq}^{i})=}}\sum\limits_{j=0}^{N_{x}^{i}}{\sum\limits_{k=0}^{N_{y}^{i}}{f_{jk}^{i}B_{jkpq}^{i}}}$

3 数值算例

3.1 串行数值实验

 $Re={{U}_{up}}L/v$

 \begin{align} & \Delta u_{\max }^{n}= \\ & \max [\sqrt{{{(u_{i,j}^{n+1}-u_{i,j}^{n})}^{2}}+{{(v_{i,j}^{n+1}-v_{i,j}^{n})}^{2}}}]＜{{10}^{-7}} \\ \end{align}

 图 2 不同雷诺数情况下的流线 Figure 2 Streamline in different Reynolds number
 图 3 不同雷诺数情况下的收敛 Figure 3 Convergence graph in different Reynolds number

3.2 并行实验设计

OpenMP采用fork-join并行执行模式，如图 4所示：OpenMP程序首先由主线程执行。需要并行计算时，主线程派生出子线程来执行并行任务。并行代码执行结束，子程序退出或者阻塞，不再工作，控制流程回到单独的主线程中。定义在成对的fork和join之间的区域，称为并行域。

 图 4 fork-join并行执行模式 Figure 4 fork-join parallel execution model

 图 5 整体刚度矩阵并行计算流程 Figure 5 Parallel computing process of global stiffness matrix
3.3 并行实验性能分析

 图 6 不同雷诺数下的加速比 Figure 6 Speed-up ratio in different Reynolds number
4 结语

1) 插值精度较高，对于具有波动性的方程，只需要很少的点就可也模拟出原函数。

2) 基函数是无限光滑的指数函数，具有无限可微性，适用于大雷诺数下湍流的计算。

3) 基函数形式简单，与待插函数形式无关。

4) 运用并行计算可以在一定程度上缩短计算时间，提高计算效率。

 [1] 樊枫, 徐国华, 史勇杰. 基于N-S方程的剪刀式尾桨前飞状态气动力计算研究[J]. 空气动力学学报, 2014, 32 (4) : 527-533. ( FAN F, XU G H, SHI Y J. Computation research on aerodynamic force of scissors tail-rotor in forward flight based on the N-S equations[J]. Acta Aerodynamica Sinica, 2014, 32 (4) : 527-533. ) [2] REKATSINAS C S, NASTOS C V, THEODOSIOU T C, et al. A time-domain high-order spectral finite element for the simulation of symmetric and anti-symmetric guided waves in laminated composite strips[J]. Wave Motion, 2015, 53 : 1-19. doi: 10.1016/j.wavemoti.2014.11.001 [3] SAMARATUNGA D, JHA R, GOPALAKRISHNAN S. Wavelet spectral finite element for modeling guided wave propagation and damage detection in stiffened composite panels[J]. Structural Health Monitoring, 2016, 15 (3) : 317-334. doi: 10.1177/1475921716640468 [4] AXELSSON O, FAROUQ S, NEYTCHEVA M. A preconditioner for optimal control problems, constrained by Stokes equation with a time-harmonic control[J]. Journal of Computational & Applied Mathematics, 2016, 310 : 5-18. [5] 王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2003 : 5 -9. ( WANG M C. Finite Element Method[M]. Beijing: Tsinghua University Press, 2003 : 5 -9. ) [6] RANNACHER R. Finite Element Methods for the Incompressible Navier-Stokes Equations[M]//Fundamental Directions in Mathematical Fluid Mechanics. Berlin:Springer, 2000:191-293. [7] 杨建宏. 定常Navier-Stokes方程的三种两层稳定有限元算法计算效率分析[J]. 数值计算与计算机应用, 2011, 32 (2) : 117-124. ( YANG J H. Computation efficiency on three kinds of two-level stabilized finite element methods for stationary Navier-Stokes equations[J]. Journal of Numerical Methods and Computer Applications, 2011, 32 (2) : 117-124. ) [8] 于佳平, 史峰, 李剑, 等. 应用四边形单元稳定化有限元方法求解定常不可压缩流动问题[J]. 工程数学学报, 2012, 29 (2) : 309-316. ( YU J P, SHI F, LI J, et al. Numerical applications of the new stabilized quadrilateral finite element method for stationary incompressible flows[J]. Chinese Journal of Engineering Mathematics, 2012, 29 (2) : 309-316. ) [9] 黄鹏展, 何银年, 冯新龙. 解Stokes特征值问题的一种两水平稳定化有限元方法[J]. 应用数学和力学, 2012, 33 (5) : 588-597. ( HUANG P Z, HE Y N, FENG X L. A two-level stabilized finite element method for the Stokes eigenvalue problem[J]. Applied Mathematics and Mechanics, 2012, 33 (5) : 588-597. ) [10] 许进超. 数值格式的稳定性,相容性和收敛性[J]. 中国科学(数学), 2015, 45 (8) : 1205-1216. ( XU J C. On the stability, consistency and convergence of numerical schemes[J]. Science in China(Series A), 2015, 45 (8) : 1205-1216. ) [11] CHACON REBOLLO T, GIRAULT V, MURAT F, et al. Analysis of a coupled fluid-structure model with applications to hemodynamics[J]. SIAM Journal on Numerical Analysis, 2016, 54 (2) : 994-1019. doi: 10.1137/140991509 [12] BIRGERSSON F, FERGUSON N S, FINNVEDEN S. Application of the spectral finite element method to turbulent boundary layer induced vibration of plates[J]. Journal of Sound & Vibration, 2003, 259 (4) : 873-891. [13] GOPALAKRISHNAN S, CHAKRABORTY A, MAHAPATRA D R. Spectral Finite Element Method[M]. London: Springer, 2008 : 177 -217. [14] 刘宏, 傅德薰, 马延文. 迎风紧致格式与驱动方腔流动问题的直接数值模拟[J]. 中国科学(A辑), 1993, 23 (6) : 657-665. ( LIU H, FU D X, MA Y W. Direct numerical simulation on the problem of upwind compact scheme and driven cavity flow[J]. Science in China (Series A), 1993, 23 (6) : 657-665. ) [15] 向新民. 谱方法的数值分析[M]. 北京: 科学出版社, 2000 : 48 -90. ( XIANG X M. Numerical Analysis of Spectral Method[M]. Beijing: Science Press, 2000 : 48 -90. ) [16] 黄冬冬.用Fourier谱方法求解二维椭圆方程Dirichlet边值问题[D].长春:吉林大学,2010. ( HUANG D D. The Fourier spectral method for solving two-dimensional elliptic partial differential equations with Dirichlet boundary condition[D]. Changchun:Jilin University, 2010. ) [17] 王健平. 有限谱有限元法求解二维不可压缩Navier-Stokes方程[J]. 力学研究, 2014, 3 (3) : 33-42. ( WANG J P. A finite spectral finite element method for incompressible Navier-Stokes equations[J]. International Journal of Mechanics Research, 2014, 3 (3) : 33-42. doi: 10.12677/IJM.2014.33004 ) [18] GRILLET A M, YANG B, KHOMAMI B, et al. Modeling of viscoelastic lid driven cavity flow using finite element simulations[J]. Journal of Non-Newtonian Fluid Mechanics, 1999, 88 (1/2) : 99-131. [19] 陈雪江, 秦国良, 徐忠. 谱元法和高阶时间分裂法求解方腔顶盖驱动流[J]. 计算力学学报, 2002, 19 (3) : 281-285. ( CHEN X J, QIN G L, XU Z. Spectral element method and higher-order fractional step method for lid driven flow in closed square cavity[J]. Chinese Journal of Computational Mechanics, 2002, 19 (3) : 281-285. ) [20] GAROOSI F, BAGHERI G, RASHIDI M M. Two phase simulation of natural convection and mixed convection of the nanofluid in a square cavity[J]. Powder Technology, 2015, 275 : 239-256. doi: 10.1016/j.powtec.2015.02.013 [21] OSSWALD G A, GHIA K N, GHIA U. Simulation of dynamic stall phenomenon using unsteady Navier-Stokes equations[J]. Computer Physics Communications, 1991, 65 .