不可壓Navier—Stokes方程的一階有限元解法

卞正寧　羅建輝

摘 要：不可壓NavierStokes方程求解的困難之一在于如何確定壓力場(chǎng)并且同時(shí)要滿(mǎn)足不可壓條件.壓力項(xiàng)在連續(xù)性方程中并不出現(xiàn)，但是卻對(duì)速度起約束作用.為了解決這一問(wèn)題，對(duì)于粘性不可壓流動(dòng)，提出了以速度和應(yīng)力為基本變量，不含壓力項(xiàng)的一階流體動(dòng)力學(xué)方程系統(tǒng)及對(duì)應(yīng)的積分形式.采用有限元方法，對(duì)于速度和應(yīng)力進(jìn)行同階插值，對(duì)于非線(xiàn)性對(duì)流項(xiàng)，采用牛頓迭代法進(jìn)行處理，對(duì)于時(shí)間項(xiàng)采用后向歐拉方法.基于FreeFem++平臺(tái)，對(duì)兩平行平板間的穩(wěn)態(tài)粘性流動(dòng)及二維非定常圓柱繞流進(jìn)行了數(shù)值計(jì)算.分別通過(guò)和精確解及標(biāo)準(zhǔn)算例的對(duì)比，驗(yàn)證了方法的可行性和有效性.采用不含壓力項(xiàng)的一階系統(tǒng)，避免了連續(xù)性方程中不含壓力項(xiàng)給不可壓縮NavierStokes方程求解帶來(lái)的困難.

關(guān)鍵詞：一階系統(tǒng)；不可壓縮流體；壓力；應(yīng)力精度

中圖分類(lèi)號(hào)：O351.2 文獻(xiàn)標(biāo)識(shí)碼：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

一直以來(lái)，基于原始變量（速度和壓力）的NavierStokes方程是求解流體力學(xué)問(wèn)題的主流提法.因?yàn)樵诜匠讨星髮?dǎo)的最高階次為二階，NavierStokes方程也可以稱(chēng)作是二階系統(tǒng).對(duì)于不可壓NavierStokes方程的數(shù)值求解，其困難之一在于確定壓力場(chǎng)并且同時(shí)要滿(mǎn)足不可壓條件.壓力項(xiàng)在連續(xù)性方程中并不出現(xiàn)，但是卻對(duì)速度起約束作用.不可壓縮流動(dòng)求解的這一困難，已經(jīng)有一些學(xué)者提出了好幾種不同的解決方案，如Harlow和Welch^{[1]給出的求解壓力泊松方程得到壓力項(xiàng)的方法，Patankar和Spalding[2]提出的SIMPLE （SemiImplicit Method for PressureLinked Equation）方法.Chorin[3-4]提出的虛擬壓縮法和投影法，以及采用渦量和流函數(shù)[5]（二維）或渦量速度[6]（三維）作依賴(lài)變量，消去控制方程中的壓力項(xiàng)的方法等.}

將有限元方法用于不可壓縮流體動(dòng)力學(xué)問(wèn)題時(shí)，如果采用以速度和壓力為基本變量的二階系統(tǒng)混合有限元方法求解，速度和壓力的插值函數(shù)不能任取，要滿(mǎn)足LadyzhenskayaBabuskaBrezzi（LBB）條件，又稱(chēng)infsup條件^{[7].常用的處理不可壓流體問(wèn)題的有限元方法有Streamline Upwind PetrovGalerkin法（SUPG）[8]，Galerkin leastsquares （GLS）[9]，Pressure Stabilized PetrovGalerkin （PSPG）[10]， Characteristicbased split（CBS）[11-12] 等等.Schwarz[13]等提出了以速度、應(yīng)力和壓力為基本變量的一階系統(tǒng)最小二乘有限元方法.因?yàn)樵诜匠讨兄缓形粗康囊浑A導(dǎo)數(shù)，最小二乘有限元方法這一類(lèi)方法可以稱(chēng)作一階系統(tǒng).目前，無(wú)論是一階系統(tǒng)還是二階系統(tǒng)，現(xiàn)有的解法都將壓力項(xiàng)作為顯式的變量進(jìn)行求解.}

由于以上各種方法都將壓力項(xiàng)作為顯式的變量進(jìn)行求解，無(wú)法繞開(kāi)連續(xù)性方程中不含壓力項(xiàng)給求解帶來(lái)的難題.二階系統(tǒng)還有一個(gè)明顯的缺陷，即應(yīng)力是由速度求導(dǎo)得到的.相對(duì)于結(jié)點(diǎn)變量本身的精度，有限元方法通過(guò)求導(dǎo)運(yùn)算得到導(dǎo)數(shù)精度會(huì)比結(jié)點(diǎn)變量本身的精度降低一階^{[14].由于NavierStokes方程是非線(xiàn)性偏微分方程，在迭代求解過(guò)程中要反復(fù)使用到對(duì)速度的偏導(dǎo)數(shù).反復(fù)迭代會(huì)導(dǎo)致這種精度下降的多次積累.}

湖南大學(xué)學(xué)報(bào)（自然科學(xué)版）2013年

第7期卞正寧等：不可壓NavierStokes方程的一階有限元解法

對(duì)于粘性不可壓流動(dòng)，本文提出了以速度和應(yīng)力為基本變量的一階解法.基于有限元方法，利用FreeFem++^{[15]平臺(tái)，對(duì)于應(yīng)力和速度均采用線(xiàn)性元的同階插值，對(duì)于非線(xiàn)性對(duì)流項(xiàng)，采用牛頓迭代法進(jìn)行處理，對(duì)于時(shí)間項(xiàng)采用后向歐拉方法.對(duì)兩平行平板間的穩(wěn)態(tài)粘性流動(dòng)及二維非定常圓柱繞流進(jìn)行了數(shù)值計(jì)算.分別通過(guò)和精確解及標(biāo)準(zhǔn)算例[16]的對(duì)比，驗(yàn)證了方法的可行性和有效性.}

1 不可壓縮流體動(dòng)力學(xué)基本方程

1） 動(dòng)量方程：

ρ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）連續(xù)性方程：

ujxj=0（4）

4） 在邊界S=Su+Sσ上，速度邊界條件和應(yīng)力邊界條件分別為：

uiSu=i，σjiljSσ=i（5a，b）

其中Su為速度邊界，Sσ為應(yīng)力邊界.

5）初始條件為：

ui（0）=u0i（5c）

將式（3）代入式（1）得：

ρDuiDt=ρFi+μ2uixjxj-pxi

（6）

2 不含壓力的一階方程

由式（3），對(duì)正應(yīng)力求和得：

σii=2μuixi-pδii（7）

將連續(xù)方程（4）代入式（7）得壓力

p=-σiiδii（8）

將式（8）回代到式（3）得

σij=μuixj+ujxi+ σkkδkkδij

（9）

由式（1）和式（9）組成了不含壓力解法的基本方程，邊界條件和初始條件仍采用式（5）表示.在解出應(yīng)力后，壓力p可由式（8）表示.

在一階方程中連續(xù)性方程沒(méi)有以顯式的形式出現(xiàn).下面將證明，對(duì)于不可壓粘性流動(dòng)，一階方程隱含了滿(mǎn)足連續(xù)性方程.

利用式（9），對(duì)正應(yīng)力求和得：

σii=2μuixi+σkk（10）

即μujxj=0.若μ=0，則式（10）為恒等式.若μ≠0，則：

ujxj=0（11）

即連續(xù)性方程（4）得到滿(mǎn)足.以上分析表明，對(duì)于不可壓無(wú)粘流動(dòng)，滿(mǎn)足了式（9），不能導(dǎo)出隱含了滿(mǎn)足連續(xù)性方程（4）.對(duì)于不可壓粘性流動(dòng)，如果滿(mǎn)足了式（9），就隱含了滿(mǎn)足連續(xù)性方程（4）.所以，不含壓力的一階方程，只是對(duì)于不可壓粘性流動(dòng)成立.

3 不含壓力一階方程的積分形式

為了進(jìn)行有限元分析，必須首先建立不含壓力一階方程的積分形式^[17].

設(shè)權(quán)函數(shù)為δui，δσij.假設(shè)速度邊界條件（5a）已經(jīng)事先滿(mǎn)足，則

δ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）亦可寫(xiě)為

∫Ωσjixj+ρFi-ρDuiDtδui dΩ+

∫Ωσi-σjiljδuidΩ+∫Ωσij2μδσij-

12μσkkδkkδσkk-12uixj+ujxiδσijdΩ=0（14）

4 基于一階方程的有限元分析

采用積分形式（13）建立有限元的離散格式.對(duì)于式（13）中非線(xiàn)性對(duì)流項(xiàng)的線(xiàn)性化，采用牛頓迭代方法，即：

ukjukixj=ukjuk-1ixj+uk-1jukixj-uk-1juk-1ixj（15）

其中k表示當(dāng)前迭代步，k-1表示前一迭代步.對(duì)于式（13）中時(shí)間項(xiàng)的處理采用后向歐拉方法，即：

ut=un-un-1Δt（16）

其中n表示當(dāng)前時(shí)間步，n-1表示上一時(shí)間步。經(jīng)過(guò)時(shí)間離散和線(xiàn)性化處理的動(dòng)量方程（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 基于一階方程的平板間粘性流動(dòng)精確解

如圖1所示，兩個(gè)相互平行的無(wú)限大平板間充滿(mǎn)不可壓流體.上板以常速度U滑動(dòng)，下板靜止不動(dòng).已知壓力梯度為：

6 基于FreeFem++的數(shù)值計(jì)算

6.1 無(wú)限大平板間定常粘性流動(dòng)數(shù)值解

對(duì)圖1所示問(wèn)題，取10 m×1 m的計(jì)算域，粘性系數(shù)μ和密度ρ均取作1， 常數(shù)P取為-8， 上板拖動(dòng)速度U=1 m/s，出口邊界條件為σx x=10=0. 網(wǎng)格在x方向平均分為100段，在y方向平均分為10段，單元采用P2單元（三角形六節(jié)點(diǎn)）.

取x=5 m處截面上u的數(shù)值結(jié)果與精確解進(jìn)行比較，見(jiàn)表1.由表1可看出，數(shù)值解與精確解相同.

6.2 二維非定常圓柱繞流數(shù)值計(jì)算

采用圓柱繞流標(biāo)準(zhǔn)算例^{[16]作為數(shù)值算例.計(jì)算域和邊界條件如圖2所示：}

其中，標(biāo)準(zhǔn)算例的結(jié)果^{[16]通過(guò)多家研究機(jī)構(gòu)分別采用有限差分法、有限體積法、有限元法、格子玻爾茲曼法等而得到.通過(guò)與標(biāo)準(zhǔn)算例的比較，計(jì)算結(jié)果吻合較好.}

7 結(jié) 論

對(duì)于粘性不可壓流動(dòng)，本文提出了不含壓力項(xiàng)的一階流體動(dòng)力學(xué)方程系統(tǒng).基于有限元方法，驗(yàn)證了一階系統(tǒng)的可行性和有效性.

本文的研究具有以下特色：

1） 在本文的方法中，壓力項(xiàng)沒(méi)有作為顯式的變量出現(xiàn)，避免了不可壓縮NS方程的連續(xù)性方程中不含壓力項(xiàng)給方程的求解帶來(lái)的難題.

2） 在不可壓縮NS方程常規(guī)的有限元求解中，一般采用二階系統(tǒng)，應(yīng)力精度比速度精度低一個(gè)階次.采用本文的一階系統(tǒng)，應(yīng)力精度和速度精度屬于同一階.避免了求導(dǎo)數(shù)運(yùn)算帶來(lái)的偏導(dǎo)數(shù)精度下降的問(wèn)題.

參考文獻(xiàn)

[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]袁駟， 肖嘉， 葉康生. 線(xiàn)法二階常微分方程組有限元分析的EEP超收斂計(jì)算 [J]. 工程力學(xué)，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.