三維可壓縮Navier-Stokes方程的間斷Galerkin有限元方法研究

秦望龍，呂宏強(qiáng)，*，伍貽兆，陳正武

(1.南京航空航天大學(xué)航空宇航學(xué)院，江蘇南京210016;2.中國(guó)空氣動(dòng)力學(xué)研究與發(fā)展中心，四川綿陽(yáng)621000)

三維可壓縮Navier-Stokes方程的間斷Galerkin有限元方法研究

秦望龍1，呂宏強(qiáng)1，*，伍貽兆1，陳正武2

(1.南京航空航天大學(xué)航空宇航學(xué)院，江蘇南京210016;2.中國(guó)空氣動(dòng)力學(xué)研究與發(fā)展中心，四川綿陽(yáng)621000)

拓展了二維間斷Galerkin(DG)有限元方法研究，將該數(shù)值方法用于三維可壓縮歐拉方程和Navier-Stokes方程的求解。基于六面體網(wǎng)格單元，采用插值方法將物面的四邊形面網(wǎng)格單元構(gòu)造為彎曲面網(wǎng)格單元，更好地表述了真實(shí)物面特征;物面邊界相鄰體網(wǎng)格單元相應(yīng)構(gòu)造為高階體網(wǎng)格單元，其余體網(wǎng)格單元采用八節(jié)點(diǎn)六面體單元，以較小的計(jì)算代價(jià)使網(wǎng)格滿足DG方法計(jì)算需求。通過(guò)對(duì)三維帶bump管道內(nèi)流、圓球繞流以及旋轉(zhuǎn)流線體繞流進(jìn)行的數(shù)值求解，驗(yàn)證了邊界彎曲方法的可行性及DG方法的高精度特性。此外，由于采用了隱式計(jì)算方法，僅需較少的時(shí)間步就能迭代收斂。

間斷有限元;Navier-Stokes方程;高精度;隱式方法;邊界彎曲

0 引言

近年來(lái)，隨著科研工作者對(duì)計(jì)算精度的要求越來(lái)越高，高精度數(shù)值方法受到越來(lái)越多的關(guān)注。在眾多的高精度方法中，間斷有限元方法由于精度高且插值模板小，易于實(shí)現(xiàn)網(wǎng)格和階數(shù)自適應(yīng)，以及并行計(jì)算效率高等優(yōu)點(diǎn)，在計(jì)算流體力學(xué)、計(jì)算聲學(xué)、計(jì)算電磁學(xué)領(lǐng)域都得到了廣泛的關(guān)注和發(fā)展［1］。

間斷有限元方法發(fā)展于20世紀(jì)70年代，由Reed和Hill在其解決中子運(yùn)輸方程的論文中提出。20世紀(jì)80年代后期和90年代，Cockburn和Shu等人［2-5］發(fā)展了Runge-Kutta Discontinuous Galerkin (RKDG)方法，求解了非線性一維守恒律方程和方程組、高維守恒律方程和方程組，并給出了部分收斂性理論證明后，該方法才引起了人們注意，并開(kāi)始逐漸應(yīng)用于計(jì)算流體力學(xué)領(lǐng)域。目前DG方法不僅應(yīng)用于雙曲守恒律方程，還被廣泛應(yīng)用于求解橢圓型方程、對(duì)流擴(kuò)散方程、KdV方程、Maxwell方程、可壓縮Navier-Stokes(N-S)方程等［6-18］。Bassi等［6-8］采用DG方法結(jié)合k-ω兩方程湍流模型求解了雷諾平均N-S方程，Luo等［10-12］研究了基于加權(quán)本質(zhì)無(wú)振蕩(Weighted Essentially Non-Oscillatory，WENO)的重構(gòu)間斷有限元方法，邱建賢等［13］針對(duì)WENO限制器在DG方法中的應(yīng)用展開(kāi)大量研究;張來(lái)平等［14］發(fā)展了靜動(dòng)態(tài)混合重構(gòu)的DG/FV混合格式。

間斷Galerkin(DG)方法通過(guò)提高單元階數(shù)來(lái)提高精度，可以在相對(duì)稀疏的網(wǎng)格單元上得到較高精度的數(shù)值解。但較稀疏的物面單元使得邊界表述不夠準(zhǔn)確，從而引入額外的數(shù)值熵增［6］，使得計(jì)算結(jié)果不準(zhǔn)確。為了對(duì)邊界進(jìn)行高階表述，文獻(xiàn)［6，15］采用高階等參單元對(duì)二維和三維Euler方程進(jìn)行了數(shù)值求解，文獻(xiàn)［16-17］則借助于CAD工具對(duì)全場(chǎng)計(jì)算單元進(jìn)行了高階表述。為了減少計(jì)算代價(jià)，文獻(xiàn)［18-21］采用一種邊界彎曲方法在較稀疏的二維網(wǎng)格上得到了高精度的數(shù)值解。

本文作為二維DG方法研究［22-23］的拓展，將DG方法用于求解三維可壓縮氣動(dòng)方程組。與以往工作不同點(diǎn)在于:將邊界網(wǎng)格彎曲方法拓展到三維六面體網(wǎng)格單元，根據(jù)四邊形面網(wǎng)格單元信息構(gòu)造出相應(yīng)的彎曲邊界信息，更準(zhǔn)確地表述了真實(shí)壁面特征，從而在相對(duì)稀疏的網(wǎng)格上求解氣動(dòng)方程組。為了提高DG方法的計(jì)算效率，文中采用了隱式計(jì)算方法。此外，文中發(fā)展的間斷有限元方法理論上可以拓展到任意高階精度。

1 控制方程

守恒形式的Navier-Stokes方程可以表示為

其中守恒變量U、對(duì)流通量F、黏性通量G可用下面的張量形式來(lái)表示:

其中ρ、p、e、h分別是流體密度、壓強(qiáng)、單位總能和單位總焓。ui=(u，v，w)是三維笛卡爾坐標(biāo)系下的速度分量，σij是黏性應(yīng)力分量，qj是熱通量分量。若不考慮黏性通量G，則方程(1)退化為守恒形式的歐拉方程［12］。

2 數(shù)值方法

2.1 空間離散

令Ω表示對(duì)計(jì)算區(qū)域的一個(gè)剖分，e為剖分單元體，文中為六面體單元，?Ωe表示單元e的邊界，ne表示單元邊界?Ωe的外法矢。在方程(1)兩邊同乘測(cè)試函數(shù)ω，在計(jì)算域積分后可得:

其中Pk(e)表示單元e上定義的至多k次多項(xiàng)式。求解方程(1)的半離散間斷有限元方法即為:求解Uh∈Vh，使得對(duì)于任意單元和任意測(cè)試函數(shù)ωh∈Vh，有:

方程(5)中存在二階偏導(dǎo)數(shù)項(xiàng)，這里采用混合方法進(jìn)行求解，將變量梯度看作額外變量，引入輔助變量Θ =U，得到如下方程組:

從式(8)可以看出，輔助變量可以看成單元梯度解和修正量R的和:

面的梯度修正量re可以表述為:

方程右邊表示在單元某個(gè)面上進(jìn)行面積分，可以看出;

采用以上方法求出方程(6)中的輔助變量，方程(7)中未知量?jī)H剩下單元變量U為未知量。數(shù)值通量II中通量函數(shù)包含無(wú)黏數(shù)值通量Hc和黏性數(shù)值通量Hv。文中無(wú)黏數(shù)值通量Hc采用LLF數(shù)值通量函數(shù)［24］，黏性數(shù)值通量Hv則取左右單元的通量平均值。

2.2 隱式時(shí)間離散

式(7)的半離散系統(tǒng)可以寫(xiě)成常微分方程組:

M是塊對(duì)角矩陣，R(u)是總殘值，u是待求未知變量。采用牛頓方法進(jìn)行線化，得到:

采用塊高斯賽德?tīng)柗椒?BGS)方法對(duì)方程(17)進(jìn)行求解:

其中，［Ae］表示單元e的對(duì)角塊矩陣，［Be］表示非對(duì)角塊矩陣，表示與單元e相鄰的每個(gè)單元f在本時(shí)間步的值。式(19)中，C∈［0，1］為穩(wěn)定因子，文中統(tǒng)一取0.5。

3 物面處理方法

3.1 彎曲物面構(gòu)造方法

與二維物面彎曲構(gòu)造的方法相似［22］，將三維物面邊界的四邊形單元從計(jì)算空間(x，y，z)轉(zhuǎn)換到參考空間(ξ，η，ζ)，在參考坐標(biāo)平面構(gòu)造曲面:

系數(shù)根據(jù)點(diǎn)的坐標(biāo)和加權(quán)法矢信息(式21)得到。

其中pi為計(jì)算空間四邊形單元的頂點(diǎn)，n1(pi)、n2(pi)為頂點(diǎn)的加權(quán)偏導(dǎo)數(shù)項(xiàng)轉(zhuǎn)換到參考平面后的數(shù)值。將參考空間構(gòu)造得到的曲面映射到原(x，y，z)計(jì)算空間即得到重構(gòu)的彎曲物面邊界。區(qū)別于二維曲面重構(gòu)后的連續(xù)和一階導(dǎo)數(shù)連續(xù)的特性［18，22］，三維重構(gòu)曲面僅有連續(xù)特性。具體曲面重構(gòu)步驟如下:

1)在(x，y，z)計(jì)算空間根據(jù)加權(quán)方法得到四邊形頂點(diǎn)的偏導(dǎo)數(shù)項(xiàng)Ni，i=1，2，3。

2)假設(shè)物面四邊形的四個(gè)頂點(diǎn)處于同一平面，取其中三點(diǎn)pi(i=1，2，3)從計(jì)算空間(x，y，z)映射到(ξ，η，ζ)空間使?jié)M足式(21)，可以得到兩個(gè)坐標(biāo)空間的映射矩陣T。根據(jù)映射矩陣T，將各頂點(diǎn)的法矢Ni轉(zhuǎn)換到參考平面得到ni，i=1，2，3。根據(jù)ni計(jì)算得到n1(pj)和n2(pj):

3)根據(jù)方程(21)～(23)，求出系數(shù)a1～a12。

4)將(ξ，η，ζ)空間下的曲面信息映射到原(x，y，z)計(jì)算空間，得到高階表述后的面單元(圖1)。

圖1 圓球的邊界彎曲示意圖Fig.1 Illustration of curved boundary representation method on a sphere

實(shí)際計(jì)算中，為了盡可能滿足步驟2中“四個(gè)頂點(diǎn)處于同一平面”的假設(shè)，在彎曲弧度較大的地方需要適當(dāng)加密網(wǎng)格。

圖2給出了曲面構(gòu)造后圓球Z=0截面與真實(shí)圓的誤差，圖3給出了誤差隨角度的分布情況。該算例中Z=0截面的誤差最大不足1%，實(shí)際計(jì)算中可以通過(guò)增加物面網(wǎng)格單元個(gè)數(shù)進(jìn)一步減小曲面構(gòu)造帶來(lái)的誤差。

圖2 曲面重構(gòu)后Z=0截面的誤差Fig.2 Illustration of the deviation on Z=0 plane

圖3 曲面重構(gòu)后Z=0截面的誤差隨角度的分布Fig.3 Illustration of the deviation according to the angle on Z=0 p lane

3.2 高階單元映射

對(duì)于非物面單元，根據(jù)六面體單元的八個(gè)頂點(diǎn)進(jìn)行線性映射:

由于物面邊界進(jìn)行了彎曲重構(gòu)，物面體網(wǎng)格單元采用32點(diǎn)進(jìn)行高階映射(圖4)。映射關(guān)系如下:

根據(jù)已有坐標(biāo)信息求出未知系數(shù)αxijk、αyijk和αzijk，即可得到單元映射雅克比矩陣及其他所需的映射關(guān)系。

圖4 六面體單元映射點(diǎn)示意圖Fig.4 Illustration of mapping points for 32-node hexahedral element

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

4.1 三維含Bump管道內(nèi)流

采用該內(nèi)流問(wèn)題驗(yàn)證邊界彎曲方法對(duì)間斷有限元格式的精度影響。該算例計(jì)算區(qū)域的長(zhǎng)度，寬度和高度分別為3、0.5、0.8。計(jì)算邊界示意圖見(jiàn)圖5，入口采用亞聲速來(lái)流，馬赫數(shù)M∞=0.5，出口為自由出流。兩側(cè)面和上表面為對(duì)稱(chēng)邊界條件，下表面為滑移邊界。下表面Bump的表達(dá)式為z=0.0625e－25x2，x∈(－1.5，1.5)。由于該流動(dòng)問(wèn)題是求解歐拉方程，且?guī)缀魏土鲃?dòng)均光滑，所以理論上流場(chǎng)等熵，文中采用熵增作為衡量精度的標(biāo)準(zhǔn)。

圖5 含Bum p管道內(nèi)流邊界示意圖Fig.5 Illustration of boundary conditions for subsonic flow through a channel w ith a bum p on the lower surface at M∞=0.5

熵增ε的定義如下:

采用三套連續(xù)剖分加密的網(wǎng)格，網(wǎng)格點(diǎn)分布為11×5×3、21×9×5、41×17×9，單元數(shù)分別為80、640、5120，如圖6。

圖6含Bum p管道內(nèi)流計(jì)算網(wǎng)格Fig.6 A sequence of three successively refined meshes used for com puting subsonic flow through a channel w ith a bump

圖7 為該算例的精度計(jì)算結(jié)果，橫坐標(biāo)表示網(wǎng)格的尺度。比較邊界彎曲修正前后的精度曲線可知，當(dāng)對(duì)彎曲幾何的表述不足時(shí)，計(jì)算得到的熵增較大，且階數(shù)較高時(shí)存在精度損失。邊界彎曲修正方法提高了該算例的幾何表述，從而空間離散達(dá)到了格式的設(shè)計(jì)精度。

圖7三維含Bump管道內(nèi)流精度測(cè)試收斂速率曲線Fig.7 Rates of convergence for subsonic flow through a channel w ith a bump

圖8 為邊界彎曲之后密網(wǎng)格上計(jì)算得到的三階精度(p=2)的壓力等值線圖，可以看出，流場(chǎng)的等值線較光滑且對(duì)稱(chēng)性較好。圖9為邊界彎曲后密網(wǎng)格上采用隱式計(jì)算方法(p=0～2)計(jì)算收斂所需的迭代步和計(jì)算時(shí)間，由于采用前一階的計(jì)算結(jié)果作為下一階數(shù)值計(jì)算的初值，計(jì)算時(shí)間和迭代步數(shù)大幅減少，密度殘值收斂到-10量級(jí)僅需20個(gè)牛頓步。隨著階數(shù)的提高，自由度呈倍增加，計(jì)算所需的時(shí)間相應(yīng)呈倍增加，在計(jì)算資源有限的情況下文中只計(jì)算到三階精度(p=2)。

圖8 三維含Bum p管道內(nèi)流壓力等值線圖Fig.8 Pressure contours for subsonic flow through a channel w ith a bump

圖9 三維含Bump管道內(nèi)流密度殘值收斂曲線Fig.9 Logarithm ic density residual versus time step and CPU time for subsonic flow through a channel w ith a bump

4.2 圓球黏性流動(dòng)

選取低雷諾數(shù)圓球繞流算例驗(yàn)證邊界彎曲方法在黏性流動(dòng)中的適用性。自由來(lái)流條件為:馬赫數(shù)M∞=0.5，雷諾數(shù)Re∞=118。取半模進(jìn)行數(shù)值計(jì)算，網(wǎng)格單元總數(shù)為9300，物面單元總數(shù)為192(圖10)。

圖10圓球黏性流動(dòng)計(jì)算網(wǎng)格Fig.10 M esh used for computing lam inar flow past a sphere

圖11 為壁面彎曲前后三階精度(p=2)下計(jì)算得到的球表面及Z=0對(duì)稱(chēng)面的壓力等值線圖?？梢钥闯觯瑴?zhǔn)確的壁面表述對(duì)DG方法較為重要，壁面彎曲后計(jì)算得到的等值線較為光滑和對(duì)稱(chēng)。圖12為計(jì)算得到的Z=0對(duì)稱(chēng)面的壓力系數(shù)曲線，經(jīng)過(guò)壁面彎曲修正后得到的壓力系數(shù)曲線較為光滑連續(xù)。

圖13為壁面彎曲前后三階精度(p=2)的密度殘值下降曲線。可以看出，物面彎曲對(duì)隱式方法的收斂效率同樣有較大影響。

圖11 圓球黏性流動(dòng)壓力等值線Fig.11 Com puted pressure contours in the flow field

圖12 圓球黏性流動(dòng)Z=0截面計(jì)算壓力系數(shù)曲線比較Fig.12 Comparison of computed pressure coefficient on Z=0 plane

圖13 三維圓球黏性流動(dòng)密度殘值收斂曲線Fig.13 Logarithm ic density residual versus time step for flow past a sphere

4.3 旋轉(zhuǎn)流線體繞流

選取High-order CFD Workshop［25］的三維旋轉(zhuǎn)流線體層流算例驗(yàn)證文中的邊界彎曲方法對(duì)于較復(fù)雜曲面的適用性。計(jì)算自由來(lái)流條件為:馬赫數(shù)M∞= 0.5，雷諾數(shù)Re∞=5000，迎角α=1°。計(jì)算網(wǎng)格單元總數(shù)為18294，物面單元總數(shù)為582(圖14)，在物面前緣和后緣進(jìn)行網(wǎng)格局部加密。圖15和圖16分別給出了三階精度(p=2)的計(jì)算等密度線圖和等馬赫線圖，即使物面網(wǎng)格相對(duì)稀疏，采用邊界彎曲方法后得到的結(jié)果依然較為光滑。圖17為本算例的密度殘值收斂曲線，由于采用了高效的隱式計(jì)算方法，殘值在30個(gè)牛頓步內(nèi)均收斂到-6量級(jí)以下。

圖14 旋轉(zhuǎn)流線體黏性流動(dòng)計(jì)算網(wǎng)格Fig.14 M esh used for computing lam inar flow past a stream lined body

圖15 旋轉(zhuǎn)流線體黏性流動(dòng)密度等值線圖Fig.15 Density contours for lam inar flow past a stream lined body

圖16 旋轉(zhuǎn)流線體黏性流動(dòng)馬赫數(shù)等值線圖Fig.16 M ach number contours for lam inar flow past a stream lined body

圖17 三維旋轉(zhuǎn)流線體黏性流動(dòng)密度殘值收斂曲線Fig.17 Logarithm ic density residual versus time step for lam inar flow past astream lined body

5 結(jié)論

本文將間斷有限元方法拓展到三維可壓縮氣動(dòng)方程組的求解中。與以往一些數(shù)值方法的區(qū)別在于:在三維情況下對(duì)邊界四邊形網(wǎng)格單元進(jìn)行了彎曲重構(gòu)，更準(zhǔn)確地表述了物面特征，從而使得DG方法在相對(duì)稀疏的網(wǎng)格上就能得到高精度的數(shù)值解。同時(shí)采用了魯棒可靠的隱式方法，縮短了計(jì)算時(shí)間，提高了計(jì)算效率。對(duì)Euler方程和Navier-Stokes方程數(shù)值求解的結(jié)果表明，文中的壁面彎曲方法能較好地應(yīng)用于間斷有限元方法且有很好的魯棒性。后續(xù)將設(shè)計(jì)高效的并行算法進(jìn)一步提高計(jì)算效率，同時(shí)考慮加入湍流模型研究更復(fù)雜的流動(dòng)問(wèn)題。

［1］Wang Z J.High-order methods for the Euler and Navier-Stokes equations on unstructured grids［J］.Progress in Aerospace Sciences，2007，43(1):1-41.

［2］Cockburn B，Shu C W.The Runge-Kutta discontinuous Galerkin method for conservation laws V:multidimensional systems［J］.Journal of Computational Physics，1998，141(2):199-224.

［3］Cockburn B，Shu C W.The local discontinuous Galerkin method for time-dependent convection-diffusion systems［J］.SIAM Journal on Numerical Analysis，1998，35(6):2440-2463.

［4］Cockburn B，Shu C W.Runge-Kutta discontinuous Galerkin methods for convection-dominated problems［J］.Journal of scientific computing，2001，16(3):173-261.

［5］Cockburn B，Li F，Shu C W.Locally divergence-free discontinuous Galerkin methods for the Maxwell equations［J］.Journal of Computational Physics，2004，194(2):588-610.

［6］Bassi F，Rebay S.High-order accurate discontinuous discontinuous finite element solution of the 2D Euler equations［J］.Journal of Computational Physics，1997，138(2):251-285.

［7］Bassi F，Rebay S.A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations［J］.Journal of Computational Physics，1997，131(2):267-279.

［8］Bassi F，Crivellini A，Rebay S，et al.Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes and k-ω turbulence model equations［J］.Computers＆Fluids，2005，34(2):507-540.

［9］Lyu Hongqiang，Sun qiang，Qin Wanglong.3D numerical solution of aero-noise with high-order discontinuous Galerkin method［J］.Journal of Nanjing University of Aeronautics＆Astronautics，2013，30(3):227-231.

［10］Xia Y，Luo H，Nourgaliev R.An implicit Hermite WENO reconstruction-based discontinuous Galerkin method on tetrahedral grids［J］.Computers＆Fluids，2014，96:406-421.

［11］Luo H，Xia Y，Li S，et al.A Hermite WENO reconstruction-based discontinuous Galerkin method for the Euler equations on tetrahedral grids［J］.Journal of Computational Physics，2012，231(16): 5489-5503.

［12］Luo H，Luo Luqing，Nourgaliev R，et al.A reconstructed discontinuous Galerkin method for the compressible Navier-Stokes equations on arbitrary grids［J］.Journal of Computational Physics，2010，229(2):6961-6978

［13］Zhu J，Qiu J.WENO Schemes and their application as limiters for RKDG methods based on trigonometric approximation spaces［J］.Journal of Scientific Computing，2012:1-39.

［14］Zhang Laiping，Li Ming，Liu Wei，et al.Recent development of high order DG/FV hybrid methods［J］.Acta Aerodynamica Sinica，2014，32(6):717-726.(in Chinese)張來(lái)平，李明，劉偉，等.基于非結(jié)構(gòu)/混合網(wǎng)格的高階精度DG/FV混合方法研究進(jìn)展［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2014，32 (6):717-726.

［15］Li S.A parallel discontinuous Galerkin method with physical orthogonal basis on curved elements［J］.Procedia Engineering，2013，61:144-151.

［16］Wang L，Anderson W K，Erwin J T，et al.Solutions of high-order methods for three-dimensional compressible viscous flows［C］// 42nd AIAA Fluid Dynamics Conference and Exhibit.New Orleans，2012.

［17］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

［18］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.

［19］Lübon C，Ke?ler M，Wagner S.A parallel CFD solver using the discontinuous Galerkin approach［M］//High Performance Computing in Science and Engineering，Garching/Munich 2007.Springer Berlin Heidelberg，2009:291-302.

［20］Yu Jian，Yan Chao.Discontinuous Galerkin method based on artificial viscosity［J］.Acta Aerodynamica Sinica，2013，31(3): 371-375.(in Chinese)于劍，閻超.基于人工粘性的間斷Galerkin有限元方法［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2013，31(3):371-375.

［21］Xia Yidong，Wu Yizhao，Lyu Hongqiang，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］Qin Wanglong，Lyu Hongqiang，Wu Yizhao.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.

［23］Qin Wanglong，Lyu Hongqiang，Wu Yizhao.Discontinuous Galerkin solution of RANS equations on curved mesh［J］.Acta Aerodynamica Sinica，2014，32(5):581-586.(in Chinese)秦望龍，呂宏強(qiáng)，伍貽兆.彎曲網(wǎng)格上的間斷有限元湍流數(shù)值解法研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2014，32(5):581-586.

［24］Toro E F，Spruce M，SpearesW.Restoration of the contact surface in the HLL-Riemann solver［J］.Shock Waves，1994，4(2):25-34.

［25］High-order CFD Workshop.Problem C 2.3.Analytical 3D Body of Revolution［OL/EB］.http://www.as.dlr.de/hiocfd/case_c2.3.html

Discontinuous Galerkin method for 3-D com pressible Navier-Stokes equations

Qin Wanglong1，Lyu Hongqiang1，*，Wu Yizhao1，Cheng Zhengwu2
(1.College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China;2.China Aerodynamics Research and Development Center，Mianyang 621000，China)

A curved-boundary based discontinuous Galerkin(DG)method is developed for solving three-dimensional compressible Euler and N-S equations on hexahedral grids.In this method，the quadrilateral face elements are reconstructed to be curved with polynomial interpolation approach，which is better to represent the real boundary.With high-order volume elements clustering only around the boundary surface，this method is easy to implement and requires a small amount of extra computations.Numerical experiments on a variety of flow problems demonstrate that DG method can obtain high-order accurate solutions on relatively coarse grids with the presented curved boundary representation approach.It is worth noting that with an implicit time integration method，converging solutions can be achieved within several time steps.

discontinuous Galerkin method;Navier-Stokes equations;high-order method;implicit method;curved boundary

V211.3

Adoi:10.7638/kqdlxxb-2015.0060

0258-1825(2016)05-0617-08

2015-05-08;

2015-09-09

國(guó)家自然科學(xué)基金(11272152);航空基金(20152752033);氣動(dòng)噪聲控制重點(diǎn)實(shí)驗(yàn)室開(kāi)放課題;江蘇高校優(yōu)勢(shì)學(xué)科建設(shè)工程資助項(xiàng)目

秦望龍(1988-)，男，江蘇南京，博士研究生，研究方向:計(jì)算流體力學(xué)，間斷有限元方法.E-mail:qinwanglong@126.com

呂宏強(qiáng)*，博士，教授.E-mail:hongqiang.lu@nuaa.edu.cn

秦望龍，呂宏強(qiáng)，伍貽兆.三維可壓縮Navier-Stokes方程的間斷Galerkin有限元方法研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2016，34(5): 617-624.

10.7638/kqdlxxb-2015.0060 Qin W L，Lyu H Q，Wu Y Z.Discontinuous Galerkin method for 3-D compressible Navier-Stokes equations［J］.Acta Aerodynamica Sinica，2016，34(5):617-624.