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

秦望龍，呂宏強，伍貽兆

（南京航空航天大學 航空宇航學院，江蘇 南京 210016）

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

秦望龍，呂宏強，伍貽兆

（南京航空航天大學 航空宇航學院，江蘇 南京 210016）

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

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

0 引 言

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

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

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

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

1 控制方程

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

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

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

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

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

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

2 數(shù)值方法

2.1 數(shù)值離散

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

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

2.2 時間積分

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

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

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

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

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

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

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

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

4 曲面物面距

為了不額外引入方程而準確計算積分點的曲面物面距，我們對物面邊進行遍歷，求出積分點到物面的最小距離d。

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

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

圖2 高斯積分點的分類Fig.2 Possible cases of Gaussian quadrature points

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

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

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

5.1 NACA0012翼型算例

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

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

圖5給出了計算空間點的物面距的等勢線，顯示范圍為［0．001，0．5］。

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

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

圖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翼型算例

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

圖8 RAE翼型計算網(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）采用高階間斷有限元方法對雷諾平均Navier-Stokes方程進行了數(shù)值離散和求解，湍流模型采用Spalart-Allmaras單方程湍流模型，并對其進行了部分修正，保證了數(shù)值求解的穩(wěn)定性和魯棒性。

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

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

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

［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）．

呂宏強，伍貽兆，周春華等．稀疏非結(jié)構(gòu)網(wǎng)格上的亞音速流高精度數(shù)值模擬［J］．航空學報，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］．力學學報，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）

秦望龍，呂宏強，伍貽兆．基于混合網(wǎng)格的高階間斷有限元黏流數(shù)值解法［J］．力學學報，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）

夏軼棟，伍貽兆，呂宏強，等．高階間斷有限元法的并行計算研究［J］．空氣動力學學報，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

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

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

呂宏強，博士，副教授．E-mail：hongqiang．lu＠126．com

秦望龍，呂宏強，伍貽兆．彎曲網(wǎng)格上的間斷有限元湍流數(shù)值解法研究［J］．空氣動力學學報，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．