孫 強(qiáng)，呂宏強(qiáng)，伍貽兆

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

基于高階物面近似的自適應(yīng)間斷有限元法歐拉方程數(shù)值模擬

孫 強(qiáng)，呂宏強(qiáng)*，伍貽兆

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

將高階間斷有限元與網(wǎng)格自適應(yīng)相結(jié)合，于非結(jié)構(gòu)網(wǎng)格上求解二維Euler方程。將數(shù)值解多項(xiàng)式的高階項(xiàng)貢獻(xiàn)用人工粘性項(xiàng)系數(shù)的形式進(jìn)行量化，網(wǎng)格自適應(yīng)過程中以人工粘性項(xiàng)系數(shù)作為網(wǎng)格自適應(yīng)的指示器。在系數(shù)達(dá)到設(shè)定的上限閥值的區(qū)域進(jìn)行網(wǎng)格加密，在系數(shù)達(dá)到設(shè)定的下限閥值的區(qū)域?qū)⒌^程中加密過的網(wǎng)格稀疏以減少網(wǎng)格量。所有自適應(yīng)均在高階曲線逼近真實(shí)物面的基礎(chǔ)上進(jìn)行，以保證數(shù)值結(jié)果的精度。典型數(shù)值算例結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行了對比，表明采用該自適應(yīng)間斷有限元法可以保證以盡可能小的計(jì)算量得到高精度結(jié)果。

高階間斷有限元;自適應(yīng)方法;Euler方程;人工粘性;物面高階近似

0 引言

間斷Galerkin有限元法(DGM)發(fā)展于20世紀(jì)70年代［1］。由于該方法能夠構(gòu)造高精度的格式，并且易于實(shí)現(xiàn)自適應(yīng)計(jì)算和并行計(jì)算，已成為計(jì)算流體力學(xué)領(lǐng)域研究熱點(diǎn)之一。Bassi等在1997年用DGM成功求解了二維Euler方程和N-S方程［2-3］;Cockburn等提出Runge-Kutta間斷Galerkin方法(RKDG)［4-5］; Nguyen等將S-A湍流模型成功應(yīng)用于DGM［6］。國內(nèi)呂宏強(qiáng)等在稀疏的非結(jié)構(gòu)網(wǎng)格上求解了二維Euler方程和線化Euler方程［7-9］;張來平等發(fā)展了靜動態(tài)混合重構(gòu)的DG/FV混合格式［10］;于劍等對人工粘性和NS方程進(jìn)行了研究［11-12］。

自適應(yīng)方法根據(jù)研究對象的不同和解的性質(zhì)自動調(diào)整網(wǎng)格分布或者數(shù)值解的逼近精度。在流體力學(xué)計(jì)算中，解通常在局部變化較大，其他大部分區(qū)域是光滑的［12］，此時使用自適應(yīng)方法可以有效降低計(jì)算量。自適應(yīng)方法通常被分為三類:對局部網(wǎng)格加密或稀疏的h型自適應(yīng);調(diào)整局部單元內(nèi)逼近多項(xiàng)式次數(shù)的p型自適應(yīng);調(diào)整網(wǎng)格單元位置保持單元總數(shù)不變的r型自適應(yīng)。

由于間斷有限元方法通常將未知數(shù)表達(dá)為高階形式，在相同的網(wǎng)格上計(jì)算量與有限體積法和有限差分法相比成倍增加，故而自適應(yīng)方法與間斷有限元方法相結(jié)合勢在必行。Yang等提出了一種于非結(jié)構(gòu)網(wǎng)格上模擬界面的自適應(yīng)方法［13］;Flaherty和Remacle研究小組則首次將網(wǎng)格度量場引入自適應(yīng)間斷有限元方法中［14］。徐云、吳迪、蔚喜軍［15-16］等則根據(jù)后誤差估計(jì)提出了不同的h型網(wǎng)格自適應(yīng)方法。

在高階情況下(p≥2)，間斷有限元對于物面形狀的表達(dá)十分敏感［3］，若物面邊界以直線表達(dá)，則無法得到收斂解。本文設(shè)計(jì)了一種h型自適應(yīng)高階間斷有限元法，物面邊均采用高階曲線逼近真實(shí)物面以保證數(shù)值解的精度。將數(shù)值解多項(xiàng)式的高階項(xiàng)貢獻(xiàn)以人工粘性［17］系數(shù)ε的形式量化，以系數(shù)ε作為網(wǎng)格自適應(yīng)的指示器。數(shù)值模擬過程中產(chǎn)生的非線性方程組采用牛頓法進(jìn)行迭代，每個迭代步產(chǎn)生的大型線性系統(tǒng)采用塊高斯-賽德爾方法求解［18］。

1 間斷有限元數(shù)值離散

二維守恒形式的Euler方程為

其中，U為守恒變量，對流通量F(U)=(f(U)，g (U))，表達(dá)式分別為:

式中，ρ為密度;v1和v2為x、y方向的速度分量;E為單位總能;P為壓強(qiáng);h為單位總焓。

式(1)兩邊乘檢驗(yàn)函數(shù)V，在計(jì)算域運(yùn)用分部積分后，弱解形式為:

將計(jì)算域劃分為網(wǎng)格{e}，其中e表示網(wǎng)格單元。在每個單元e內(nèi):

式(4)中，φi(x，y)為基函數(shù)。將式(4)代入式(3)中，得

式(6)稱為p階間斷有限元離散［7］，p為式(4)中基函數(shù)的階數(shù)。

由于間斷有限元法允許解在單元的邊界處不連續(xù)，因而數(shù)值通量的定義不是唯一的。此處的數(shù)值通量直接借用傳統(tǒng)有限體積法的處理方法，即將式(6)中的通量函數(shù)F(Uh)·n用數(shù)值通量代替，其中和分別表示該單元和其相鄰單元在該邊界處的值。本文所有計(jì)算都采用LLF格式數(shù)值通量:

在物面和遠(yuǎn)場分別采用無穿透和無反射邊界條件。式(6)最終的離散形式可寫為:

其中，M是質(zhì)量矩陣，如果間斷有限元采用正交基函數(shù)，則是對角矩陣;R(u)是殘值;u包含Uh中所有的未知系數(shù)。對式(8)的求解采用牛頓-塊高斯賽德爾方法求解［18］。

2 人工粘性

人工粘性法［18］的思想就是在原Euler方程中加入一個耗散修正項(xiàng)，從而抑制激波附近的非物理振蕩。方程(1)變?yōu)槿缦滦问?

式(9)中，參數(shù)ε控制粘性項(xiàng)的大小。

式中，(·，·)e表示L2(Ωe)空間的標(biāo)準(zhǔn)內(nèi)積。這里預(yù)估Se的值能達(dá)到1/p4數(shù)量級。

單元的光順指示器確定后，粘性項(xiàng)系數(shù)通過以下函數(shù)定義:

式(11)中，se= log10Se，參數(shù) ε0= De/p，s0= log10(1/p4)，κ是人為定義常數(shù)。其中，De是計(jì)算元單元e的特征尺度。

由式(10)和式(11)可以看出，se是由U中高階項(xiàng)決定，表征了U的高階項(xiàng)的貢獻(xiàn)大小。解的光順程度越差se越大。εe只是作為se的分段函數(shù)，在數(shù)值表征上意義相同，它抹平光順區(qū)域的較小的se，突出非光順區(qū)域，使計(jì)算和表達(dá)更簡便。

3 網(wǎng)格自適應(yīng)

3.1 網(wǎng)格加密

在計(jì)算過程中單元的ε達(dá)到設(shè)定好的上限閥值則認(rèn)為該單元需要進(jìn)行加密，采用三角形單元四分法對網(wǎng)格單元進(jìn)行加密，如圖1所示。

圖1 單元加密Fig.1 Element refinement

圖中實(shí)線單元為非物面單元，可以直接取三條邊的中點(diǎn)連線剖分加密。

Bassi指出高階間斷有限元對于物面形的表達(dá)精度十分敏感［3］，因而本文中的物面單元的物面邊均采用高階曲線逼近真實(shí)物面:在計(jì)算過程中設(shè)計(jì)物面邊

的高階曲線(本文中為6階)，即使網(wǎng)格十分稀疏，構(gòu)造好的物面曲邊依然可以十分精確的表達(dá)真實(shí)物面，如圖1中左圖虛線所示。

在物面單元需要進(jìn)行剖分時，利用已構(gòu)造好的高階曲線找到逼近真實(shí)物面的物面邊的中點(diǎn)，然后連接另外兩條邊的中點(diǎn)生成四個新的網(wǎng)格單元，并且新生成的兩個物面單元的物面邊也用高階曲線擬合，如圖1中右圖虛線所示。

網(wǎng)格加密需要遵守以下原則:

1)相鄰單元的被剖分次數(shù)相差最多為1。

2)當(dāng)單元的幾何尺寸細(xì)化到設(shè)定好的閥值之后就不可繼續(xù)剖分。

3.2 網(wǎng)格稀疏

在計(jì)算過程中，部分單元被加密后可能在進(jìn)一步的迭代過程中解變得十分光滑，需要對其進(jìn)行粗化以減少網(wǎng)格量。當(dāng)加密過的單元的四個子單元ε均達(dá)到設(shè)定好的下限閥值時，則這四個子單元進(jìn)行合并。如圖2所示。

圖2 單元的粗化Fig.2 Element coarsening

圖2中，實(shí)線表示非物面單元的合并過程。若粗化后的單元為物面單元，物面邊依舊采用高階曲線逼近真實(shí)物面，如圖中虛線所示。

網(wǎng)格粗化是網(wǎng)格加密的逆向操作，需要遵守以下原則:

1)初始單元不可被粗化。

2)只有被加密過的單元才可以被粗化。

3)當(dāng)鄰單元已經(jīng)比判斷需要粗化的四個單元剖分次數(shù)大1時，不可進(jìn)行粗化。

3.3 網(wǎng)格數(shù)據(jù)結(jié)構(gòu)

為使程序的可移植性更高，將網(wǎng)格自適應(yīng)作為模塊進(jìn)行操作。雖然二叉樹結(jié)構(gòu)存儲加密和粗化網(wǎng)格直觀高效，但自適應(yīng)過程中增加與減少的點(diǎn)和線編號復(fù)雜且對于程序的結(jié)構(gòu)不利。因而本文所有單元、點(diǎn)、線的存儲結(jié)構(gòu)均為鏈表結(jié)構(gòu)，加密過程中單元的增加操作如圖3所示。

E表示需要加密的單元，生成的四個子單元中間單元依舊為編號E，其余三個則排在總單元后三個，點(diǎn)與邊的操作類似。在這種操作下，網(wǎng)格模塊與程序流場計(jì)算模塊的接口不受自適應(yīng)操作的影響，從而提高程序的可移植性。

同樣的，網(wǎng)格稀疏的過程單元減少過程如圖4所示。

圖3 單元增加Fig.3 Elements increase

圖4 單元減少Fig.4 Elements decrease

其中E和三個fi表示由加密操作生成的四個子單元，合并后生成的單元依舊編號為E，而fi后的單元編號依次減3。

在網(wǎng)格加密過程中，每個子單元均設(shè)置兩個標(biāo)記，用于記錄所在的初始單元號和本身的加密次數(shù)。在后續(xù)的網(wǎng)格稀疏過程中，根據(jù)這兩個信息和單元相鄰關(guān)系進(jìn)行單元合并。

3.4 數(shù)據(jù)傳遞

網(wǎng)格加密和粗化的過程中，使各物理量在子單元內(nèi)的積分值等于合并單元內(nèi)的積分值，從而保證物理量的守恒。假設(shè)合并單元為E，各子單元為ei，uE(x，y)為合并單元的守恒變量的函數(shù)，uei(x，y)為子單元守恒變量的函數(shù)，則需滿足以下方程:

其中m為合并單元剖分后的子單元個數(shù)。

4 數(shù)值結(jié)果

4.1 圓柱繞流測試算例

為驗(yàn)證物面構(gòu)造與網(wǎng)格自適應(yīng)結(jié)合的效果，本文首先對圓柱繞流(Ma=0.38、α=0)這一經(jīng)典算例進(jìn)行了數(shù)值模擬。

前人［3，8］做圓柱繞流計(jì)算為獲得精度較高的結(jié)果，網(wǎng)格均人為設(shè)計(jì)為規(guī)則分布。本文為驗(yàn)證物面構(gòu)造和自適應(yīng)間斷有限元的優(yōu)點(diǎn)，初始網(wǎng)格使用了分布較隨意、單元尺寸較極端的計(jì)算網(wǎng)格。

圖5給出了圓柱繞流網(wǎng)格，該網(wǎng)格共有80個單元，物面僅有8個點(diǎn)。

圖5 圓柱繞流計(jì)算網(wǎng)格Fig.5 Mesh around a circle

圖6 ε分布和自適應(yīng)后的網(wǎng)格Fig.6 ε distribution and mesh after three times adaption

圖6(a)為在原始網(wǎng)格上p=4時ε分布，可以看出主要分布在流場數(shù)值可能出現(xiàn)數(shù)值解變化比較劇烈的區(qū)域。對這一區(qū)域進(jìn)行網(wǎng)格自適應(yīng)，2次自適應(yīng)后的網(wǎng)格如圖6(b)所示?？梢钥闯鼍W(wǎng)格主要在ε集中區(qū)域自適應(yīng)，從而僅需對這一區(qū)域的少量網(wǎng)格進(jìn)行自適應(yīng)操作。自適應(yīng)前后的網(wǎng)格對比，發(fā)現(xiàn)自適應(yīng)操作確實(shí)基于真實(shí)物面曲線，自適應(yīng)后的網(wǎng)格更加逼近真實(shí)物面。

圖7給出了p=4情況下的圓柱繞流等馬赫線圖和壓強(qiáng)系數(shù)分布，可以看出壓強(qiáng)系數(shù)分布光滑并且對稱。

圖7 圓柱繞流等馬赫線圖和Cp分布Fig.7 Mach isolines and Cpdistribution

4.2 跨聲速NACA0012翼型繞流

采用上述自適應(yīng)間斷有限元法給出了經(jīng)典算例跨聲速(Ma=0.8、α=1.25)下NACA0012翼型繞流的數(shù)值結(jié)果。

圖8給出的NACA0012翼型的網(wǎng)格共有448個單元，而物面單元只有32個。

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

圖9給出了引入人工粘性未加入網(wǎng)格自適應(yīng)的粘性項(xiàng)系數(shù)ε的分布與馬赫云圖。

圖9 未加自適應(yīng)的ε分布及等馬赫線圖Fig.9 ε distribution and Mach isolines without adaption

從圖中可以看出雖然網(wǎng)格十分稀疏，但依然可以在一個網(wǎng)格尺度內(nèi)捕捉到激波，并且ε的分布集中在激波區(qū)域。

圖10給出了1次、2次自適應(yīng)后的網(wǎng)格變化。圖11給出了網(wǎng)格自適應(yīng)過程中局部網(wǎng)格的細(xì)化和粗化的動態(tài)變化情況。

從圖中可以看出，自適應(yīng)過程中激波區(qū)域的網(wǎng)格加密，而遠(yuǎn)離激波區(qū)域在精度低的時候加密了的單元則合并降低網(wǎng)格量。

圖12給出了自適應(yīng)結(jié)束時的網(wǎng)格分布和ε分布。

圖13給出了添加自適應(yīng)后的五階精度(p=4)的等馬赫線圖和表面壓強(qiáng)系數(shù)分布。

圖10 1次和2次自適應(yīng)后的網(wǎng)格Fig.10 Mesh after once and twice adaption

圖11 1次和2次自適應(yīng)后的局部網(wǎng)格Fig.11 Local mesh after once and twice adaption

圖12 自適應(yīng)結(jié)束的網(wǎng)格和ε分布Fig.12 Mesh and ε distribution after adaption

圖13 等馬赫線圖和Cp分布曲線Fig.13 Mach isolines and Cpdistribution

從圖中可以看出，自適應(yīng)后的數(shù)值結(jié)果較原始網(wǎng)格上的數(shù)值精度有顯著提高，并且激波捕捉更精確接近實(shí)驗(yàn)值。

4.3 超聲速NACA0012翼型繞流

為驗(yàn)證方法的適應(yīng)性，采用上述自適應(yīng)間斷有限元法給出了超聲速(Ma=5.0、α=0)下NACA0012翼型繞流的數(shù)值結(jié)果(p=4)。

初始網(wǎng)格如圖14，初始網(wǎng)格數(shù)為660。自適應(yīng)計(jì)算結(jié)束后的馬赫云圖及網(wǎng)格、ε分布如圖15所示。

從圖中可以看出，自適應(yīng)操作主要集中在頭部激波區(qū)域。

圖14 局部網(wǎng)格Fig.14 Local mesh

圖15 馬赫云圖及ε分布Fig.15 Mach contour and ε distribution

5 結(jié)論

(1)將數(shù)值解的高階項(xiàng)貢獻(xiàn)量化后作為自適應(yīng)指示器可行并直觀有效。

(2)即使在很稀疏的初始網(wǎng)格基礎(chǔ)上，使用基于真實(shí)物面逼近的網(wǎng)格自適應(yīng)方法也可以獲得高精度的數(shù)值解。

(3)本文設(shè)計(jì)的自適應(yīng)方法提高了數(shù)值解精度的同時有效地控制了網(wǎng)格量。

對于非定常問題和三維問題，自適應(yīng)方法在節(jié)省網(wǎng)格工作量和提高計(jì)算效率方面作用更為突出。同時具有加密和粗化的功能在非定常計(jì)算中會得到更有意義的應(yīng)用，該部分研究正在進(jìn)行中。在粘性流動的情況下可以在不引入人工粘性項(xiàng)的前提下計(jì)算人工粘性項(xiàng)系數(shù)，將其作為網(wǎng)格自適應(yīng)判據(jù)，這部分工作也正在進(jìn)行。

［1］ Reed W H，Hill T R.Triangular mesh methods for the neutron transport［R］.Los Alamos Scientific Laboratory Report LA-UR-73-479，1973.

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

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

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

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

［6］ Nguyen N C，Persson P O，Peraire J.RANS solutions using high order discontinuous Galerkin methods［R］.AIAA 2007-0914.

［7］ Lu Hongqiang，Wu Yizhao，Zhou Chunhua，et al.High resolution of subsonic flows on coarse grids［J］.Acta Aeronautica et Astronautica Sinica，2009，30(2):200-204.(in Chinese).

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

［8］ Lu Hongqiang，Zhu Guoxiang，Song Jiangyong，et al.High-order discontinuous Galerkin solution of linearized Euler equations［J］.Chinese Journal of Theoretical and Applied Mechanics，2011，43 (3):621-624.(in Chinese).

呂宏強(qiáng)，朱國祥，宋江勇，等.線化歐拉方程的高階間斷有限元數(shù)值解法研究［J］.力學(xué)學(xué)報(bào)，2011，43(3):621-624.

［9］ Lu H.High-order discontinuous Galerkin solution of low-Re viscous flows［J］.Modern Physics Letters B，2009，23(03):309-312.

［10］Zhang Laiping，Liu Wei，He Lixin，et al.A class of discontinuous Galerkin/finite volume hybrid schemes based on the“static re-construction”and“dynamic re-construction”［J］.Chinese Journal of Theoretical and Applied Mechanics，2010，42(6):1013-1022.(in Chinese).

張來平，劉偉，賀立新，等.基于靜動態(tài)混合重構(gòu)的DG/FV混合格式［J］.力學(xué)學(xué)報(bào)，2010，42(6):1013-1022.

［11］Yu J，Yan C.An artificial diffusivity discontinuous Galerkin scheme for discontinuous flows［J］.Computers＆Fluids，2013，75(20): 56-71.

［12］Yu Jian，Yan Chao.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.

［13］Yang X，James A J，Lowengrub J，et al.An adaptive coupled levelset/volume-of-fluid interface capturing method for unstructured triangular grids［J］.Journal of Computational Physics，2006，217 (2):364-394.

［14］Remacle J F，Li X，Shephard M S，et al.Anisotropic adaptive simulation of transient flows using discontinuous Galerkin methods［J］.International Journal for Numerical Methods in Engineering，2005，62(7):899-923.

［15］Xu Yun，Yu Xijun.Adaptive discontinuous Galerkin methods for hyperbolic conservation laws［J］.Chinese Journal of Computational Physics，2009，26(2):159-168.(in Chinese).

徐云，蔚喜軍.自適應(yīng)間斷有限元方法求解雙曲守恒律方程［J］.計(jì)算物理，2009，26(2):159-168.

［16］Wu Di，Yu Xijun.Adaptive discontinuous Galerkin method for euler equations［J］.Chinese Journal of Computational Physics，2010，27 (004):492-500(in Chinese).

吳迪，蔚喜軍.自適應(yīng)間斷有限元方法求解三維歐拉方程［J］.計(jì)算物理，2010，27(004):492-500.

［17］VonNeumann J，Richtmyer R D.A method for the numerical calculation of hydrodynamic shocks［J］.Journal of Applied Physics，1950，21(3):232-237.

［18］Xia Y D，Wu Y Z，Lv H Q，et al.Parallel computation of a highorder discontinuous Galerkin method on unstructured grids［J］.Acta Aerodynamica Sinica，2011，29(5):537-541.(in Chinese).

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

［19］Lu H，Berzins M，Goodyer C E，et al.Adaptive high-order discontinuous Galerkin solution of elastohydrodynamic lubrication point contact problems［J］.Advances in Engineering Software，2012，45 (1):313-324.

Adaptive discontinuous Galerkin method to solve Euler equations based on high-order approximative boundary

Sun Qiang，Lu Hongqiang*，Wu Yizhao
(College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China)

A high-order discontinuous method(DGM)is integrated with adaptive method to solve Euler equations on unstructured mesh.Contribution of the polynomial’s highest-order terms is quantified in the form of artificial viscous coefficient.The coefficient is regarded as the indicator of h-adaptivity.Elements where the coefficients are greater than the upper limit are refined.Those where the coefficients are less than the lower limit are coarsened if they have been refined.A high-order geometric approximation of curved boundaries is adopted to ensure the convergence.Numerical results of test cases are consistent with corresponding experimental ones.High accurate numerical results can be obtained with the h-adaptive method at low expense.

high-order DGM;adaptive method;Euler equations;artificial viscosity;high-order boundary approximation

V211.3

Adoi:10.7638/kqdlxxb-2014.0027

0258-1825(2015)04-0446-08

2014-04-23;

2014-07-08

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

孫強(qiáng)(1987-)，男，山東萊陽市人，博士研究生，研究方向:計(jì)算流體力學(xué).E-mail:submarineultra@126.com

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

孫強(qiáng)，呂宏強(qiáng)，伍貽兆.基于高階物面近似的自適應(yīng)間斷有限元法歐拉方程數(shù)值模擬［J］.空氣動力學(xué)學(xué)報(bào)，2015，33(4):446-453.

10.7638/kqdlxxb-2014.0027 Sun Q，Lu H Q，Wu Y Z.Adaptive discontinuous Galerkin method to solve Euler equations based on high-order approximative boundary［J］.Acta Aerodynamica Sinica，2015，33(4):446-453.