秦望龍，呂宏強(qiáng)，伍貽兆

（南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016）

彎曲網(wǎng)格上的間斷有限元湍流數(shù)值解法研究

秦望龍，呂宏強(qiáng)，伍貽兆

（南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016）

采用間斷有限元方法對(duì)雷諾平均Navier-Stokes（RANS）方程進(jìn)行了數(shù)值求解，對(duì)Spalart-Allmaras單方程湍流模型進(jìn)行了部分修正，使得求解器更加魯棒。構(gòu)造了分段高次多項(xiàng)式來逼近真實(shí)物面，同時(shí)物面附近采用多層彎曲網(wǎng)格來避免網(wǎng)格交叉，此外提出了一種快速計(jì)算積分點(diǎn)的曲面物面距的方法。采用混合網(wǎng)格對(duì)NACA0012翼型以及RAE翼型進(jìn)行了數(shù)值計(jì)算，并與實(shí)驗(yàn)數(shù)據(jù)以及前人數(shù)據(jù)進(jìn)行了對(duì)比。計(jì)算結(jié)果表明，采用物面彎曲網(wǎng)格結(jié)合修正的湍流模型方法在相對(duì)稀疏的網(wǎng)格上就能得到比較好的數(shù)值解。

間斷有限元；RANS；Spalart-Allmaras湍流模型；物面距；混合網(wǎng)格

0 引 言

近年來，隨著眾多學(xué)者的深入研究，間斷Galerkin有限元（Discontinous Galerkin，DG）方法求解雷諾平均Navier-Stokes（Reynolds-averaged Navier-Stokes，RANS）方程取得了迅速發(fā)展［1-11］。目前采用DG方法求解RANS方程的湍流模型主要有一方程Spalart-Allmaras模型［4-10］、兩方程Wilcoxk-ω模型［1-3］以及兩方程剪切力運(yùn)輸（SST）模型［11］。

由于RANS方程本身的非光順性質(zhì)，間斷有限元方法對(duì)其求解比較困難。尤其當(dāng)采用稀疏網(wǎng)格，單元內(nèi)會(huì)出現(xiàn)吉布斯（Gibbs）現(xiàn)象導(dǎo)致求解失敗。針對(duì)SA一方程模型的穩(wěn)定求解問題，一些學(xué)者提出了相應(yīng)的方法［4-7］。文中采用一種修正的S-A一方程模型較好的抑制了Gibbs現(xiàn)象，順利求解了RANS方程。

間斷有限元方法對(duì)物面形狀的表達(dá)精度非常敏感，如果物面采用分段線性網(wǎng)格進(jìn)行計(jì)算會(huì)導(dǎo)致求解不收斂或數(shù)值結(jié)果不光順［12-13］。為了解決以上問題，一般采用二次及以上多項(xiàng)式進(jìn)行高階近似以逼近真實(shí)物面，文獻(xiàn)［14-18］給出了不同的高階近似方法。然而在求解粘性流動(dòng)，尤其是雷諾數(shù)較大時(shí)，由于第一層網(wǎng)格厚度很小，物面高階近似會(huì)帶來網(wǎng)格的交叉纏繞，從而出現(xiàn)負(fù)體積。為了使求解順利進(jìn)行，文中采用多層彎曲網(wǎng)格遞推的方法，以很小的計(jì)算代價(jià)解決了以上問題。

Spalart-Allmaras單方程模型及SST兩方程模型的求解需要計(jì)算積分點(diǎn)的物面距。對(duì)于彎曲物面，積分點(diǎn)的物面距的求解比較困難，Liu等人采用有限元方法對(duì)Eikonal方程進(jìn)行數(shù)值離散求解了內(nèi)流問題中的物面距［19］。然而，該方法在求解外流問題時(shí)會(huì)不穩(wěn)定，Schoenawa采用流線-擴(kuò)散方法對(duì)Eikonal方程進(jìn)行離散求解，并引入人工粘性方法使其穩(wěn)定魯棒［11］。為了不額外引入計(jì)算方程，本文采用一種簡(jiǎn)單的數(shù)值方法對(duì)物面距進(jìn)行計(jì)算。

1 控制方程

可壓縮的雷諾平均Navier-Stokes方程耦合修正的一方程Spalart-Allmaras湍流模型可以寫成如下的守恒形式：

其中，Ω是有界域。守恒變量u、無粘數(shù)值通量Fc、粘性數(shù)值通量Fv可以寫成如下矢量形式：

其中，ρ、p、e0、h分別是密度、壓強(qiáng)、單位總能和單位總焓，ui＝（vx，vy）是笛卡爾坐標(biāo)系下的速度分量。壓強(qiáng)p滿足理想氣體方程：

γ是定壓比熱容Cp與定容比熱容Cv的比值，文中γ取1.4。σij是粘性應(yīng)力張量，對(duì)于牛頓流體，

其中，μ是動(dòng)力粘性系數(shù)，可以通過Sutherland公式得出。μT是湍流粘性系數(shù)：

湍流生成項(xiàng)和消散項(xiàng)中的參數(shù)分別如下：

2 數(shù)值方法

2.1 數(shù)值離散

將計(jì)算域表示成單元集合Ω＝∪kΩk，文中Ωk是四邊形單元或三角形單元。定義Vh，p是單元h上采用多項(xiàng)式基函數(shù)表示的p階有限元空間，p≥0。單元變量uh可以在Vh，p下表示，即uh∈Vh，p。對(duì)方程（1）離散，有

式（10）中，F(xiàn)（uh，?uh）是對(duì)流通量和粘性通量的和。邊界通量H=F（uh，?uh）·n的求解需要引入合適的數(shù)值通量函數(shù)。文中采用LLF數(shù)值通量求解邊界無粘數(shù)值通量［21］。對(duì)于粘性項(xiàng)，采用BR2格式（the second scheme of Bassi and Rebay）［1］進(jìn)行數(shù)值離散。

2.2 時(shí)間積分

數(shù)值離散之后，方程最終可以寫成一個(gè)常微分方程系統(tǒng)：

式（11）中，M是質(zhì)量矩陣，R（Uh）是殘值，Uh是未知變量。文中采用隱式向后歐拉方法對(duì)式（11）進(jìn)行數(shù)值求解。每個(gè)迭代步產(chǎn)生一個(gè)線性方程組：

對(duì)于式（12）的求解文中采用ILU預(yù)處理的廣義最小余量方法［22］。

3 網(wǎng)格彎曲處理

高階間斷有限元方法對(duì)物面形狀的表達(dá)精度非常敏感，文中采用三次多項(xiàng)式對(duì)物面邊進(jìn)行分段高階近似，從而逼近真實(shí)物面。求解高雷諾數(shù)粘性流動(dòng)時(shí)，在物面進(jìn)行高階近似的同時(shí)需要對(duì)相鄰?fù)鈱泳W(wǎng)格進(jìn)行高階近似才能避免網(wǎng)格交叉纏繞。

區(qū)別于文獻(xiàn)［18］中的映射方法，文中采用參數(shù)方程實(shí)現(xiàn)物面高階近似和曲面網(wǎng)格陣推。如圖1，將曲線AB表示成參數(shù)方程形式：

根據(jù)邊界點(diǎn)A、B的坐標(biāo)及加權(quán)方法算出的法矢nA、nB，可以求出上述方程的八個(gè)參數(shù)，從而實(shí)現(xiàn)物面彎曲。在此基礎(chǔ)求出AB線段1／3、2／3點(diǎn)a、b處的彎曲距離α、β。根據(jù)每一層的陣推距離K1α、K2β以及邊界點(diǎn)A1、B1，我們可以求出所有彎曲邊的高階參數(shù)方程形式。這里，K1、K2是調(diào)節(jié)因子，用于調(diào)節(jié)該面彎曲的程度，文中均取1．0。采用以上方法可以求出任意彎曲層數(shù)下彎曲邊的表達(dá)形式，進(jìn)一步提高了網(wǎng)格彎曲的自動(dòng)化，并且為計(jì)算曲面物面距中搜索點(diǎn)的法矢提供了保證。需要指出的是，在實(shí)現(xiàn)高階近似過程中，對(duì)于物面一些特殊點(diǎn)，如翼型后緣尖點(diǎn)，由于過該點(diǎn)的相鄰邊不是圓滑過渡，所以該點(diǎn)的法矢不能簡(jiǎn)單采用加權(quán)平均的方法進(jìn)行計(jì)算，否則逼近的外形會(huì)偏離真實(shí)外形。

圖1 網(wǎng)格單元高階近似示意圖Fig.1 Illustration of high order approximation

4 曲面物面距

為了不額外引入方程而準(zhǔn)確計(jì)算積分點(diǎn)的曲面物面距，我們對(duì)物面邊進(jìn)行遍歷，求出積分點(diǎn)到物面的最小距離d。

對(duì)于固定曲面物面，根據(jù)積分點(diǎn)的位置可以將其分為三類（圖2），然后根據(jù)式（14）分別求解。

當(dāng)積分點(diǎn)處于b情況時(shí)，對(duì)該曲面物面進(jìn)行二分法遍歷，直到在曲面上找到點(diǎn)g滿足n（g）·l（bg）＝cosα＝0。n（g）表示曲面上g點(diǎn)的單位法矢，l（bg）表示b點(diǎn)到g點(diǎn)的單位向量，α是兩向量的夾角。為了減少遍歷次數(shù)，文中α取不小于80°。采用以上方法對(duì)RAE翼型的物面距的計(jì)算結(jié)果如圖3，等勢(shì)線顯示范圍為［0．001，0．5］。該算例計(jì)算了2695個(gè)網(wǎng)格單元中心點(diǎn)到物面的距離，其中1320個(gè)點(diǎn)需要計(jì)算case b情形。由于對(duì)α的設(shè)置，該算例case b情形的遍歷次數(shù)基本在2次及以下。

圖2 高斯積分點(diǎn)的分類Fig.2 Possible cases of Gaussian quadrature points

圖3 RAE翼型的物面距分布Fig.3 Wall distance solution of RAE airfoil

5 數(shù)值結(jié)果與分析

采用上述高階間斷有限元方法結(jié)合物面彎曲混合網(wǎng)格單元對(duì)RANS方程進(jìn)行數(shù)值求解。

5.1 NACA0012翼型算例

計(jì)算來流條件為Ma＝0.3，α＝0°，Re＝1．85×106。高雷諾數(shù)條件決定了第一層網(wǎng)格的厚度很薄，文中第一層網(wǎng)格取弦長(zhǎng)的2×10-5，對(duì)應(yīng)Y＋＝4。計(jì)算網(wǎng)格單元如圖4，其中非結(jié)構(gòu)網(wǎng)格單元由結(jié)構(gòu)網(wǎng)格單元剖分而來。由于物面采用大長(zhǎng)寬比結(jié)構(gòu)網(wǎng)格單元，該算例計(jì)算網(wǎng)格總數(shù)控制在2592。

圖4 NACA0012翼型計(jì)算網(wǎng)格Fig.4 Mesh around NACA0012 airfoil

圖5給出了計(jì)算空間點(diǎn)的物面距的等勢(shì)線，顯示范圍為［0．001，0．5］。

圖5 NACA0012翼型的物面距分布Fig.5 Wall distance solution of NACA0012 airfoil

圖6給出了本算例的馬赫云圖及湍流粘性系數(shù)云圖，圖7給出了本算例四階精度的表面壓力系數(shù)Cp分布及摩擦阻力系數(shù)Cf分布曲線?？梢钥闯?，Cp計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好。由于Cf沒有實(shí)驗(yàn)數(shù)據(jù)，這里我們采用Fluent軟件在14000個(gè)單元上計(jì)算得到的Cf數(shù)據(jù)作為比較，可以看出曲線走勢(shì)大致吻合。

圖6 NACA0012翼型馬赫云圖及湍流粘性系數(shù)云圖Fig.6 Mach contours and eddy viscosity contours of NACA0012 airfoil

圖7 Cp和Cf分布曲線Fig.7 Cpand Cfdistribution of NACA0012 airfoil

5.2 RAE翼型算例

計(jì)算來流條件為Ma＝0.4，α＝2．79°，Re＝6．5×106。采用混合網(wǎng)格單元，總數(shù)為2695（圖8）。計(jì)算精度從1階（p＝0）到4階（p＝3），為使p階初值盡可能的靠近最終收斂解，采用p－1階的收斂解作為其初值。圖9給出了本算例p＝3時(shí)的馬赫云圖及湍流粘性系數(shù)云圖。圖10給出了本算例的表面壓力系數(shù)Cp分布及摩擦阻力系數(shù)Cf分布曲線，計(jì)算結(jié)果與參考文獻(xiàn)［8］的計(jì)算結(jié)果吻合較好。由于文中曲線數(shù)據(jù)均未采用任何光順方法，所以在單元邊界處會(huì)有局部不光滑現(xiàn)象。

圖8 RAE翼型計(jì)算網(wǎng)格Fig.8 Mesh around RAE airfoil

圖9 RAE翼型馬赫云圖及湍流粘性系數(shù)云圖Fig.9 Mach contours and eddy viscosity contours of RAE airfoil

圖10 Cp和Cf分布曲線Fig.10 Cpand Cfdistribution of RAE airfoil

6 結(jié) 論

（1）采用高階間斷有限元方法對(duì)雷諾平均Navier-Stokes方程進(jìn)行了數(shù)值離散和求解，湍流模型采用Spalart-Allmaras單方程湍流模型，并對(duì)其進(jìn)行了部分修正，保證了數(shù)值求解的穩(wěn)定性和魯棒性。

（2）采用混合網(wǎng)格單元，物面采用參數(shù)方程實(shí)現(xiàn)高階物面近似及曲面網(wǎng)格陣推，由于物面采用多層大長(zhǎng)寬比彎曲結(jié)構(gòu)網(wǎng)格單元，外部流場(chǎng)計(jì)算使用非結(jié)構(gòu)網(wǎng)格單元，在較少的網(wǎng)格上就可以得到相當(dāng)精度的數(shù)值解。

（3）給出了一種計(jì)算曲面物面距的方法，以很少的計(jì)算時(shí)間得到了滿足計(jì)算需要的曲面物面距。

（4）采用以上方法對(duì)NACA0012翼型和RAE翼型高雷諾數(shù)算例進(jìn)行了數(shù)值計(jì)算，驗(yàn)證了文中的物面距求解方法和湍流模型方法。

［1］BASSI F，CRIVELLINI A，REBAY S，et al．Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes andkωturbulence model equations［J］．Computers&Fluids，2005，34（2）：507-540．

［2］BASSI F，CRIVELLINI A，GHIDONI A，et al．High-order discontinuous Galerkin discretization of transonic turbulent flows［C］．In Proceedings of the 47th Aerospace Sciences Meeting，Oriando；AIAA-2009-180．

［3］HARTMANN R，HELD J，LEICHT T．Adjoint-based error estimation and adaptive mesh refinement for the RANS andk-ω turbulence model equations［J］．Journal of Computational Physics，2011，230（11）：4268-4284．

［4］NGUYEN N C，PERSSON P O，PERAIRE J．RANS solutions using high order discontinuous Galerkin methods［C］．The 45th Aerospace Sciences Meeting，Reno，2007；AIAA-2007-914．

［5］OLIVER T，DARMOFAL D．An unsteady adaptation algorithm for discontinuous Galerkin discretizations of the RANS equations［C］．The 18th AIAA Computational Fluid Dynamics Conference，Miami，2007；AIAA-2007-3940．

［6］LANDMANN B，KESSLER M，WAGNER S，et al．A parallel，high-order discontinuous Galerkin codes for laminar and turbulent flows［J］．Computers&Fluids，2008，37（2）：427-438．

［7］MORO D，NGUYEN N C，PERAIRE J．Navier-Stokes solution using hybridizable discontinuous Galerkin methods［C］．The 20th Computational Fluid Dynamics Conference，Honolulu，2011；AIAA-2011-3407．

［8］BURGESS N K，MAVRIPLIS D J．Robust computation of turbulent flows using a discontinuous Galerkin method［C］．The 50th Aerospace Sciences Meeting，Nashville，2012；AIAA-2012-0457．

［9］BURGESS N K，MAVRIPLIS D J．High-Order Discontinuous Galerkin Methods for Turbulent High-lift Flows［C］．The 7th International Conference on Computational Fluid Dynamics（ICCFD7），Big Island，2012，ICCFD7-4202．

［10］WANG L，ANDERSON W K，ERWIN J T，et al．Solutions of high-order methods for three-dimensional compressible viscous flows［C］．The 42nd AIAA Fluid Dynamics Conference and Exhibit，New Orleans，2012．

［11］SCHOENAWA S，HARTMANN R．Discontinuous Galerkin discretization of the Reynolds-averaged Navier-Stokes equations with the shear-stress transport model［J］．Journal of Computational Physics，2014．

［12］LüH Q，WU Y Z，ZHOU C H，et al．High resolution of subsonic flows on coarse grids［J］．Acta Aeronautica et Astronautica Sinca，2009，30（2）：200-204（in Chinese）．

呂宏強(qiáng)，伍貽兆，周春華等．稀疏非結(jié)構(gòu)網(wǎng)格上的亞音速流高精度數(shù)值模擬［J］．航空學(xué)報(bào)，2009，30（2）：200-204．

［13］LU H．High-order discontinuous Galerkin solution of low-Re viscous flows［J］．Modern Physics Letters B，2009，23（03）：309-312．

［14］LüBON C，et al．High-order boundary discretization for discontinuous Galerkin codes［C］．In Proceedings of the 24th Applied Aerodynamics Conference，F(xiàn)rancisco，2006；AIAA-2006-2822．

［15］HARTMANN R，HELD J，LEICHT T，et al．Discontinuous Galerkin methods for computational aerodynamics－3D adaptive flow simulation with the DLR PADGE code［J］．Aerospace Science and Technology，2010，14（7）：512-519．

［16］PERSSON P，PERAIRE J．Curved mesh generation and mesh refinement using Lagrangian solid mechanics［C］．In Proceedings of the 47th Aerospace Sciences Meeting，Oriando，2009；AIAA-2009-949．

［17］YU J，YAN C．Study on discontinuous Galerkin method for Navier-Stokes equations［J］．Chinese Journal of Theoretical and Applied Mechanics，2010，42（5）：962-969（in Chinese）

于劍，閻超．Navier-Stokes方程間斷Galerkin有限元方法研究［J］．力學(xué)學(xué)報(bào)，2010，42（5）：962-969．

［18］QIN W L，LüH Q，WU Y Z．High-order discontinuous Galerkin solution of N-S equations on hybrid mesh［J］．Chinese Journal of Theoretical and Applied Mechani.，2013，45（6）：987-991．（in Chinese）

秦望龍，呂宏強(qiáng)，伍貽兆．基于混合網(wǎng)格的高階間斷有限元黏流數(shù)值解法［J］．力學(xué)學(xué)報(bào)，2013，45（6）：987-991．

［19］LIU C B，NITHIARASU P，TUCKER P G．Wall distance calculation using the Eikonal／Hamilton-Jacobi equations on unstructured meshes：a finite element approach［J］．Engineering Computations，2010，27（5）：645-657．

［20］SPALART P，ALLMARASS．A one-equation turbulence model for aerodynamic flows［C］．In Proceedings of the 30th Aerospace Sciences Meeting and Exhibit，Reno，1992；AIAA-92-0439．

［21］XIA Y D，WU Y Z，LüH Q，et al．Parallel computation of a high-order discontinuous galerkin method on unstructured grids［J］．ACTA Aerodynamica Sinica，2011，29（5）：537-541．（in Chinese）

夏軼棟，伍貽兆，呂宏強(qiáng)，等．高階間斷有限元法的并行計(jì)算研究［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2011，29（5）：537-541．

［22］SAAD Y，SCHULTZ M H．GMRES：A generalized minimal residual algorithm for solving nonsymmetric linear systems［J］．SIAM Journal on Scientific and Statistical Computing，1986，7（3）：856-869．

Discontinuous Galerkin solution of RAN Sequations on curved mesh

QIN Wanglong，Lü Hongqiang，WU Yizhao
（College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China）

Discontinuous Galerkin（DG）finite element method was adopted for the numerical approximation of the Reynolds-averaged Navier-Stokes（RANS）equations with the Spalart-Allmaras turbulence model．In order to make the solver robust，the original turbulence model equation were modified accordingly．Furthermore，high order approximation of the real solid boundary was used and several layers of curved meshes were constructed to avoid inconsistent mesh cross-overs．For the computation of the distance of each quadrature point to the nearest curved wall boundaries，a fast straightforward numerical method was proposed．The DG discretization of the RANS equations were demonstrated for turbulent flows past a NACA0012 airfoil and RAE airfoil based on hybrid mesh．Numerical results indicate that highly accurate solutions can be obtained with the modified turbulent equation on coarse curved hybrid mesh．

Discontinuous Galerkin（DG）；RANS；Spalart-Allmaras turbulent model；high order approximation；hybrid mesh

V211.3

Adoi:10.7638／kqdlxxb-2014.0058

0258-1825（2014）05-0581-06

2014-06-24；

2014-08-10

國(guó)家自然科學(xué)基金（11272152）；航空科學(xué)基金（20101552018）；江蘇高校優(yōu)勢(shì)學(xué)科建設(shè)工程資助項(xiàng)目

秦望龍（1988-），男，江蘇南通人，博士研究生，主要研究方向：計(jì)算流體力學(xué)，高階間斷有限元．E-mail：qinwanglong＠126．com

呂宏強(qiáng)，博士，副教授．E-mail：hongqiang．lu＠126．com

秦望龍，呂宏強(qiáng)，伍貽兆．彎曲網(wǎng)格上的間斷有限元湍流數(shù)值解法研究［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2014，32（5）：581-586．

10．7638／kqdlxxb-2014．0058． QIN W L，LU H Q，WU Y Z．Discontinuous Galerkin solution of RAN Sequations on curved mesh［J］．ACTA Aerodynamica Sinica，2014，32（5）：581-586．