不可压Navier—Stokes方程的一阶有限元解法
2013-04-29卞正宁罗建辉
卞正宁 罗建辉
摘 要：不可压NavierStokes方程求解的困难之一在于如何确定压力场并且同时要满足不可压条件.压力项在连续性方程中并不出现，但是却对速度起约束作用.为了解决这一问题，对于粘性不可压流动，提出了以速度和应力为基本变量，不含压力项的一阶流体动力学方程系统及对应的积分形式.采用有限元方法，对于速度和应力进行同阶插值，对于非线性对流项，采用牛顿迭代法进行处理，对于时间项采用后向欧拉方法.基于FreeFem++平台，对两平行平板间的稳态粘性流动及二维非定常圆柱绕流进行了数值计算.分别通过和精确解及标准算例的对比，验证了方法的可行性和有效性.采用不含压力项的一阶系统，避免了连续性方程中不含压力项给不可压缩NavierStokes方程求解带来的困难.
关键词：一阶系统；不可压缩流体；压力；应力精度
中图分类号：O351.2 文献标识码：A
Oneorder Algorithm of Incompressible NavierStokes
Equations in Finite Element Method
BIAN Zhengning， LUO Jianhui
（College of Civil Engineering， Hunan Univ， Changsha， Hunan 410082，China ）
Abstract： One of the difficulties of the numerical solution of incompressible NavierStokes equations is the determination of the pressure field and the fulfillment of the incompressibility condition. In fact， the pressure variable is not present in continuity equation， but a constraint for the velocity field is present. In this paper， the basic variables of velocity and stress were proposed for incompressible viscous fluid， a oneorder fluid dynamics equation system without pressure term was proposed and its integral form was given to handle this problem. The stress and the velocity were interpolated by equal order finite element. The Newton iterative method was used to handle the nonlinear convective term. The backward Euler method was used to discretize the time term. A steady flow of incompressible viscous fluid between two infinite parallel plates and a Benchmark problem of incompressible viscous fluid flow around a cylinder were computed on the basis of FreeFem++. The feasibility and the effectivity of the method were verified by comparing with the analytic solution and the Benchmark results respectively. The difficulty of pressure term which is not present in continuity equation is circumvented by using oneorder system without pressure term.
Key words：oneorder system； incompressible fluid； pressure； stress accuracy
一直以来，基于原始变量（速度和压力）的NavierStokes方程是求解流体力学问题的主流提法.因为在方程中求导的最高阶次为二阶，NavierStokes方程也可以称作是二阶系统.对于不可压NavierStokes方程的数值求解，其困难之一在于确定压力场并且同时要满足不可压条件.压力项在连续性方程中并不出现，但是却对速度起约束作用.不可压缩流动求解的这一困难，已经有一些学者提出了好几种不同的解决方案，如Harlow和Welch[1]给出的求解压力泊松方程得到压力项的方法，Patankar和Spalding[2]提出的SIMPLE （SemiImplicit Method for PressureLinked Equation）方法.Chorin[3-4]提出的虚拟压缩法和投影法，以及采用涡量和流函数[5]（二维）或涡量速度[6]（三维）作依赖变量，消去控制方程中的压力项的方法等.
将有限元方法用于不可压缩流体动力学问题时，如果采用以速度和压力为基本变量的二阶系统混合有限元方法求解，速度和压力的插值函数不能任取，要满足LadyzhenskayaBabuskaBrezzi（LBB）条件，又称infsup条件[7].常用的处理不可压流体问题的有限元方法有Streamline Upwind PetrovGalerkin法（SUPG）[8]，Galerkin leastsquares （GLS）[9]，Pressure Stabilized PetrovGalerkin （PSPG）[10]， Characteristicbased split（CBS）[11-12] 等等.Schwarz[13]等提出了以速度、应力和压力为基本变量的一阶系统最小二乘有限元方法.因为在方程中只含有未知量的一阶导数，最小二乘有限元方法这一类方法可以称作一阶系统.目前，无论是一阶系统还是二阶系统，现有的解法都将压力项作为显式的变量进行求解.
由于以上各种方法都将压力项作为显式的变量进行求解，无法绕开连续性方程中不含压力项给求解带来的难题.二阶系统还有一个明显的缺陷，即应力是由速度求导得到的.相对于结点变量本身的精度，有限元方法通过求导运算得到导数精度会比结点变量本身的精度降低一阶[14].由于NavierStokes方程是非线性偏微分方程，在迭代求解过程中要反复使用到对速度的偏导数.反复迭代会导致这种精度下降的多次积累.
湖南大学学报（自然科学版）2013年
第7期卞正宁等：不可压NavierStokes方程的一阶有限元解法
对于粘性不可压流动，本文提出了以速度和应力为基本变量的一阶解法.基于有限元方法，利用FreeFem++[15]平台，对于应力和速度均采用线性元的同阶插值，对于非线性对流项，采用牛顿迭代法进行处理，对于时间项采用后向欧拉方法.对两平行平板间的稳态粘性流动及二维非定常圆柱绕流进行了数值计算.分别通过和精确解及标准算例[16]的对比，验证了方法的可行性和有效性.
1 不可压缩流体动力学基本方程
1） 动量方程：
ρDuiDt=ρFi+σjixj，DuiDt=uit+ujuixj（1）
2）几何方程和物理方程：
sij=12（uixj+ujxi）， σij=2μsij-pδij（2a，b）
几何方程和物理方程可以合并为：
σij=μuixj+ujxi-pδij（3）
3）连续性方程：
ujxj=0（4）
4） 在边界S=Su+Sσ上，速度边界条件和应力边界条件分别为：
uiSu=i，σjiljSσ=i（5a，b）
其中Su为速度边界，Sσ为应力边界.
5）初始条件为：
ui（0）=u0i（5c）
将式（3）代入式（1）得：
ρDuiDt=ρFi+μ2uixjxj-pxi
（6）
2 不含压力的一阶方程
由式（3），对正应力求和得：
σii=2μuixi-pδii（7）
将连续方程（4）代入式（7）得压力
p=-σiiδii（8）
将式（8）回代到式（3）得
σij=μuixj+ujxi+ σkkδkkδij
（9）
由式（1）和式（9）组成了不含压力解法的基本方程，边界条件和初始条件仍采用式（5）表示.在解出应力后，压力p可由式（8）表示.
在一阶方程中连续性方程没有以显式的形式出现.下面将证明，对于不可压粘性流动，一阶方程隐含了满足连续性方程.
利用式（9），对正应力求和得：
σii=2μuixi+σkk（10）
即μujxj=0.若μ=0，则式（10）为恒等式.若μ≠0，则：
ujxj=0（11）
即连续性方程（4）得到满足.以上分析表明，对于不可压无粘流动，满足了式（9），不能导出隐含了满足连续性方程（4）.对于不可压粘性流动，如果满足了式（9），就隐含了满足连续性方程（4）.所以，不含压力的一阶方程，只是对于不可压粘性流动成立.
3 不含压力一阶方程的积分形式
为了进行有限元分析，必须首先建立不含压力一阶方程的积分形式[17].
设权函数为δui，δσij.假设速度边界条件（5a）已经事先满足，则
δuiSu=0（12）
由式（1），（9） 和（5b）得不含压力一阶方程的积分形式：
∫Ωσjixj+ρFi-ρDuiDtδui dΩ+
∫Ωσi-σjiljδuidΩ+∫Ωσij2μ-
12μσkkδkkδij-12uixj+ujxiδσijdΩ=0（13）
式（13）亦可写为
∫Ωσjixj+ρFi-ρDuiDtδui dΩ+
∫Ωσi-σjiljδuidΩ+∫Ωσij2μδσij-
12μσkkδkkδσkk-12uixj+ujxiδσijdΩ=0（14）
4 基于一阶方程的有限元分析
采用积分形式（13）建立有限元的离散格式.对于式（13）中非线性对流项的线性化，采用牛顿迭代方法，即：
ukjukixj=ukjuk-1ixj+uk-1jukixj-uk-1juk-1ixj（15）
其中k表示当前迭代步，k-1表示前一迭代步.对于式（13）中时间项的处理采用后向欧拉方法，即：
ut=un-un-1Δt（16）
其中n表示当前时间步，n-1表示上一时间步。经过时间离散和线性化处理的动量方程（1）为：
ρuni-un-1iΔt+un，kjun，k-1ixj+un，k-1jun，kixj-
un，k-1jun，k-1ixj=ρFni+σnjixj（17）
将式（17）代入式（13），得一阶解法的最终积分形式为：
∫Ωσnjixj-ρuni-un-1iΔt+un，kjun，k-1ixj+un，k-1jun，kixj-
un，k-1jun，k-1ixj+ρFniδuidΩ+
∫Ωσi-σnjiljδuidΩ+
∫Ωσnij2μ-12μσnkk3δij-12unixj+unjxiδσijdΩ=0（18）
利用式（18），即可建立有限元的离散格式.
5 基于一階方程的平板间粘性流动精确解
如图1所示，两个相互平行的无限大平板间充满不可压流体.上板以常速度U滑动，下板静止不动.已知压力梯度为：
6 基于FreeFem++的数值计算
6.1 无限大平板间定常粘性流动数值解
对图1所示问题，取10 m×1 m的计算域，粘性系数μ和密度ρ均取作1， 常数P取为-8， 上板拖动速度U=1 m/s，出口边界条件为σx x=10=0. 网格在x方向平均分为100段，在y方向平均分为10段，单元采用P2单元（三角形六节点）.
取x=5 m处截面上u的数值结果与精确解进行比较，见表1.由表1可看出，数值解与精确解相同.
6.2 二维非定常圆柱绕流数值计算
采用圆柱绕流标准算例[16]作为数值算例.计算域和边界条件如图2所示：
其中，标准算例的结果[16]通过多家研究机构分别采用有限差分法、有限体积法、有限元法、格子玻尔兹曼法等而得到.通过与标准算例的比较，计算结果吻合较好.
7 结 论
对于粘性不可压流动，本文提出了不含压力项的一阶流体动力学方程系统.基于有限元方法，验证了一阶系统的可行性和有效性.
本文的研究具有以下特色：
1） 在本文的方法中，压力项没有作为显式的变量出现，避免了不可压缩NS方程的连续性方程中不含压力项给方程的求解带来的难题.
2） 在不可压缩NS方程常规的有限元求解中，一般采用二阶系统，应力精度比速度精度低一个阶次.采用本文的一阶系统，应力精度和速度精度属于同一阶.避免了求导数运算带来的偏导数精度下降的问题.
参考文献
[1] HARLOW F H， WELCH J E. Numerical calculation of timedependent viscous incompressible flow with free surface[J]. Physics of Fluids， 1965， 8（12）： 2182-2189.
[2] PATANKAR S V， SPALDING D B. A calculation procedure for heat， mass and momentum transfer in threedimensional parabolic flows[J]. International Journal of Heat and Mass Transfer， 1972， 15： 1787-1806.
[3] CHORIN A J. A numerical method for solving incompressible viscous flow problems[J]. Journal of Computational Physics， 1967， 2（1）： 12-26.
[4] CHORIN A J. Numerical solution of the NavierStokes equations[J]. Mathematics of Computation， 1968， 22（104）： 745-762.
[5] HOU T Y. Stable fourthorder streamfunction methods for incompressible flows with boundaries[J]. Journal of Computational Mathematics， 2009， 27（4）： 441-458.
[6] WU X H， WU J Z， WU J M. Effective vorticityvelocity formulations for threedimensional incompressible viscous flows[J]. Journal of Computational Physics， 1995， 122： 68-82.
[7] GELHARD T， LUBER G， OLSHANSKII M A，et al. Stabilized finite element schemes with LBBstable elements for incompressible flows[J]. Journal of Computational and Applied Mathematics，2005，177（2）：243-267.
[8] BROOKS A， HUGHES T J R. Streamline upwind/PetrovGalerkin formulations for convection dominated flows with particular emphasis on the incompressible NavierStokes equations[J]. Computer Methods in Applied Mechanics and Engineering， 1982， 32（1/3）： 199-259.
[9] HUGHES T J R， FRANCA L P， HULBERT G M. A new finite element formulation for computational fluid dynamics： VIII. The galerkin/leastsquares method for advectivediffusive equations[J]. Computer Methods in Applied Mechanics and Engineering， 1989， 73（2）： 173-189.
[10]TEZDUYAR Te. Stabilized finite element formulations for incompressible flow computations[J]. Advances in Applied Mechanics， 1992， 28（1）： 1-44.
[11]ZIENKIEWICZ O C， CODINA R. A general algorithm for compressible and incompressible flow—Part I. the split， characteristicbased scheme[J]. International Journal for Numerical Methods in Fluids， 1995， 20（8/9）： 869-885.
[12]ZIENKIEWIC O C， MORGAN K B. A general algorithm for compressible and incompressible flow—Part II. tests on the explicit form[J]. International Journal for Numerical Methods in Fluids，1995， 20（8/9）： 887-913.
[13]SCHWARZ Alexander， SCHRDER Jrg. A mixed leastsquares formulation of the NavierStokes equations for incompressible Newtonian fluid flow[J]. Proceedings in Applied Mathematics and Mechanics， 2011，11（1）：589-590.
[14]袁驷， 肖嘉， 叶康生. 线法二阶常微分方程组有限元分析的EEP超收敛计算 [J]. 工程力学，2009， 21（11）：1-9.
YUAN Si ， XIAO Jia ， YE Kangsheng. EEP superconvergent computation in FEM analysis of FEMOL second order ODES[J]. Engineering Mechanics，2009，21（11）：1-9.（In Chinese）
[15]HECHT F， PIRONNEAU O. FreeFem++ 3.13 user manual[M/OL]. http：//www.freefem.org， 2011.
[16]SCHFERA F M， TUREK S， DURST F，et al. Benchmark computations of laminar flow around a cylinder[J]. Notes on Numerical Fluid Mechanics，1996， 52： 547-566.
[17]LUO Jianhui， LI Qiusheng， LIU Guangdong. A biorthogonality relationship for threedimensional couple stress problem[J]. Science in China Series GPhysics， Mechanics & Astronomy， 2009， 52（2）： 270-276.
