郭永恒，楊 永，劉逸鷗

（西北工業(yè)大學(xué) 翼型葉柵空氣動(dòng)力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室，陜西西安 710072）

其中，ˉH 代表單元交界面處的數(shù)值通量，它反映了當(dāng)前單元和相鄰單元之間的信息傳遞，這里采用 Roe格式對它進(jìn)行處理。這樣，整個(gè)離散系統(tǒng)就可以寫成常微分方程的形式

一種基于間斷Galerkin格式的新型限制器設(shè)計(jì)

郭永恒，楊 永，劉逸鷗

（西北工業(yè)大學(xué) 翼型葉柵空氣動(dòng)力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室，陜西西安 710072）

使用高階間斷 Galerkin格式求解守恒律方程組時(shí)，激波附近的 Gibbs效應(yīng)容易導(dǎo)致非物理解的產(chǎn)生。為抑制這一現(xiàn)象，必須構(gòu)造合理的限制器對數(shù)值解進(jìn)行處理。目前間斷 Galerkin格式中的限制器多源于有限體積法，在非結(jié)構(gòu)網(wǎng)格上只對低階導(dǎo)數(shù)項(xiàng)進(jìn)行限制，對高階導(dǎo)數(shù)項(xiàng)則很難給出普適判據(jù)。文章對間斷 Galerkin解進(jìn)行廣義Fourier展開，實(shí)現(xiàn)不同頻域范圍內(nèi)的譜分解；在新的模板坐標(biāo)系下描述各階方向?qū)?shù)的變化規(guī)律；結(jié)合當(dāng)前單元和相鄰單元的信息，分層限制各階方向?qū)?shù)，實(shí)現(xiàn)對非物理解的抑制。通過求解Euler方程，對二維 Riemann問題、翼型周圍的亞、跨聲速流動(dòng)問題、前臺(tái)階問題以及超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流場激波反射問題進(jìn)行數(shù)值模擬，檢驗(yàn)了新型限制器的可靠性以及向高階格式推廣的可行性。

間斷 Galerkin格式；限制器；廣義Fourier展開；模板坐標(biāo)系；方向?qū)?shù)

0 引 言

航空航天事業(yè)的迅速發(fā)展，對計(jì)算流體力學(xué)提出了越來越高的要求，特別是在研究復(fù)雜流動(dòng)問題方面，傳統(tǒng)的有限體積法已經(jīng)很難適應(yīng)新形勢的需要。于是，許多新型的高精度算法應(yīng)運(yùn)而生，其中，間斷Galerkin格式［1］的研究與應(yīng)用在最近幾年逐漸成為一個(gè)熱點(diǎn)問題，并在計(jì)算流體力學(xué)領(lǐng)域取得了一系列令人矚目的結(jié)果。然而，當(dāng)處理間斷解問題時(shí)，Gibbs效應(yīng)很容易導(dǎo)致非物理解的產(chǎn)生，使間斷 Galerkin格式的計(jì)算過程出現(xiàn)中斷。為了避免這種現(xiàn)象的發(fā)生，人們在間斷Galerkin格式中應(yīng)用大量的限制器，如Shu的斜率限制器［2］、Hoteit改進(jìn)的 Van Leer限制器［3］、Tu和Aliabadi改進(jìn)的Van Albada-typed限制器［4］、Qiu改進(jìn)的WENO／HWENO 限制器［5］以及Krivodonova根據(jù)多元Taylor公式推導(dǎo)的限制器［6］。它們大多來源于傳統(tǒng)的有限體積法，一般只對有限元解的低階導(dǎo)數(shù)項(xiàng)進(jìn)行處理，而對高階導(dǎo)數(shù)項(xiàng)產(chǎn)生的非物理振蕩抑制效果不佳甚至無法抑制，這使它們很難推廣到高階的空間離散格式中。從某種意義上說，限制器的設(shè)計(jì)已經(jīng)成為制約高階間斷Galerkin格式應(yīng)用的一個(gè)“瓶頸”。在這一背景下，本文結(jié)合間斷Galerkin格式的具體形式，把握“局部解的方向?qū)?shù)”這一信息，并在新的模板坐標(biāo)系下對其進(jìn)行分層限制，有效避免了非物理解的產(chǎn)生，在保持高階精度不變的同時(shí)，使計(jì)算過程得以順利進(jìn)行。

1 間斷Galerkin格式的建立

直角坐標(biāo)系下二維Euler方程的守恒形式為

首先，把整個(gè)求解區(qū)域Ω 劃分成三角形非結(jié)構(gòu)單元的集合，即其中 NE 為網(wǎng)格單元的個(gè)數(shù)；同時(shí)，在第k號(hào)單元上，設(shè)局部有限元解為

其中PHj是當(dāng)前單元上的基函數(shù)，NP 為基函數(shù)的個(gè)數(shù)。根據(jù)加權(quán)余量法［7］的思想，令測試函數(shù)空間等于解函數(shù)空間，作如下內(nèi)積，得

根據(jù)Gauss公式，可以得到間斷Galerkin方程的弱守恒形式

其中，ˉH 代表單元交界面處的數(shù)值通量，它反映了當(dāng)前單元和相鄰單元之間的信息傳遞，這里采用 Roe格式對它進(jìn)行處理。這樣，整個(gè)離散系統(tǒng)就可以寫成常微分方程的形式

沿時(shí)間方向使用Runge-Kutta方法推進(jìn)，就得到式（5）的數(shù)值解。

2 對問題單元的初篩

在式（2）中，依 Gram-Schmidt過程［8］構(gòu)造規(guī)范正交基，并通過線性變換，把每個(gè)物理單元映射成計(jì)算域

對應(yīng)的基函數(shù)的具體表達(dá)式為當(dāng)選取前6個(gè)基函數(shù)時(shí)，對應(yīng)的空間離散精度為三階。從廣義Fourier級(jí)數(shù)的角度來看，正交基函數(shù)連同前面的系數(shù)QHj構(gòu)成了有限元解的低頻分量和高頻分量，這是新型限制器設(shè)計(jì)的出發(fā)點(diǎn)。按照Gibbs效應(yīng)［9］的相關(guān)理論，非物理振蕩對應(yīng)的廣義 Fourier系數(shù)QHj不能快速的衰減。 對于三階格式，在此約定，把有限元系數(shù)分為三層，即平均層｛QH1｝、一次層｛QH2，QH3｝和二 次層｛QH4，QH5，QH6｝，它 們 對應(yīng) 于 守恒變量的積分平均值、一階和二階空間導(dǎo)數(shù)分量。接著，定義無量綱的振蕩參數(shù)

這里，把β1和β2稱為一次、二次相對振蕩強(qiáng)度。首先判斷一次層，當(dāng)β1大于一定的控制閥值μ1時(shí)，需要對一次層進(jìn)行限制；然后判斷二次層，當(dāng)β2大于一定的控制閥值μ2時(shí)，需要對二次層進(jìn)行限制。在實(shí)踐過程中，我們選取密度項(xiàng)對應(yīng)的廣義Fourier系數(shù)對β1和β2進(jìn)行計(jì)算。這樣，我們就完成了對問題單元的初篩。在這個(gè)階段，控制閥值越小，初篩越敏感，在接下來的限制過程中判斷是否需要進(jìn)行重構(gòu)的計(jì)算量越大。在實(shí)踐過程中，當(dāng)這些參數(shù)處于同一量級(jí)（μ1為10-2量級(jí)，μ2為10-4量級(jí)）時(shí)，計(jì)算結(jié)果幾乎是一致的。也就是說，計(jì)算精度主要不是由問題單元的初步篩查來單獨(dú)決定的，而是結(jié)合接下來的限制過程決定的。

3 模板坐標(biāo)系下的限制器設(shè)計(jì)

非結(jié)構(gòu)網(wǎng)格單元排列無序且缺乏正交性。如果單純采用計(jì)算坐標(biāo)系ξ-o2-η，就很難反映當(dāng)前單元和相鄰單元之間的空間聯(lián)系；如果單純采用物理坐標(biāo)系x-o1-y，就很難像結(jié)構(gòu)網(wǎng)格那樣對信息傳遞的方向進(jìn)行很好的判斷。這里，我們建立模板坐標(biāo)系，并結(jié)合方向?qū)?shù)的相關(guān)信息構(gòu)造新型限制器。如圖1所示，記當(dāng)前單元和相鄰單元的重心分別為G0、G1、G2、G3，相應(yīng)的三個(gè)模板依次為（0，1，2）、（0，2，3）、（0，3，1），特別的，對于邊界處的單元，只存在一個(gè)模板。

物流企業(yè)不同于其他的制造性企業(yè)，其發(fā)生的成本費(fèi)用沒有專門的會(huì)計(jì)制度來規(guī)范。在進(jìn)行成本核算過程中，各企業(yè)往往按照現(xiàn)行的《企業(yè)會(huì)計(jì)制度》來核算。但是物流企業(yè)提供的服務(wù)時(shí)無形的，其成本核算對象不易選取，在核算方面必然存在缺陷，管理層也無法得到最精確的成本數(shù)據(jù)，無法達(dá)到管理的目的。在賬務(wù)核算時(shí)，沒有專門的會(huì)計(jì)科目來進(jìn)行歸集和分類，不易對某一環(huán)節(jié)發(fā)生的費(fèi)用進(jìn)行控制。一般來說，企業(yè)根據(jù)財(cái)務(wù)的會(huì)計(jì)科目來做成本預(yù)算，這樣不利于物流企業(yè)對成本進(jìn)行控制。大部分物流企業(yè)采用傳統(tǒng)成本法下的科目來核算，如將不同環(huán)節(jié)發(fā)生的物流成本統(tǒng)一歸集到“管理費(fèi)用”等期間費(fèi)用中，無法核算各個(gè)作業(yè)下的成本信息。

圖1 模板坐標(biāo)系示意圖Fig.1 A diagram for the mode coordinate system

以第1號(hào)模板為例，令當(dāng)前單元的重心 G0為原點(diǎn)o3，分別連接該點(diǎn)與1、2號(hào)單元的重心 G1、G2，構(gòu)成兩個(gè)線性無關(guān)的方向r和s，于是，新的模板坐標(biāo)系即定義為r-O3-S。對于一次層，我們分析其一階方向?qū)?shù)在當(dāng)前單元上可以通過復(fù)合函數(shù)的鏈導(dǎo)法則 進(jìn) 行計(jì) 算；另 一方 面的數(shù)值也可以通過有限差分法直接進(jìn) 行 計(jì)算，其 表 達(dá) 式為和 Q（2）H 1分 別 為 第 1、2 號(hào) 單 元的第一項(xiàng)廣義Fourier系數(shù)，r01和s02分別為當(dāng)前單元重心G0和第1、2號(hào)單元重心G1、G2之間的距離。設(shè)修正后的一次層系數(shù)為那么

這樣，聯(lián)立式（9）、（10）就可以求出修正后的一次層系數(shù)。對于模板1，我們定義一個(gè)加權(quán)因子w1=β1（3），即第3號(hào)單元的一次相對振蕩強(qiáng)度?？梢钥闯?，當(dāng)?shù)?號(hào)單元上的振蕩相對較強(qiáng)時(shí)，這樣的定義方式模板1中重構(gòu)的信息所占的加權(quán)比重相對更大一些。同理，可定 義于是最后的加權(quán)修正結(jié)果

為

類似的，我們可以進(jìn)一步比較高階方向?qū)?shù)的解析結(jié)果和差分結(jié)果，建立如下三個(gè)限制方程

聯(lián)立方程（13）、（14）、（15）就可以求出修正后的二次層系數(shù)，對于模板1，我們定義一個(gè)加權(quán)因子w1為w1=β（3）

2，即第3號(hào)單元的二次相對振蕩強(qiáng)度。

就這樣，我們通過對不同階次的方向?qū)?shù)的逐層分析，最終完成了對廣義Fourier系數(shù)的限制?？梢宰C明，一方面，當(dāng)略去高階項(xiàng)的限制部分后，本文的限制器形式非常接近于Shu的斜率限制器［2］，所以該限制器可以視為Shu限制器向更高階格式的推廣；另一方面，當(dāng)問題退化到一維情況或者結(jié)構(gòu)化網(wǎng)格情況時(shí)，方向?qū)?shù)可以用偏導(dǎo)數(shù)直接代替，這又和Krivodonova設(shè)計(jì)的限制器［6］非常相似，只是后者從高階部分向低階部分依次進(jìn)行限制，本文的限制器的處理次序與之相反，所以它又可以視為Krivodonova限制器向非結(jié)構(gòu)網(wǎng)格的推廣。很顯然，這種分層限制的處理方式可以方便的推廣到更高階的格式中。

4 數(shù)值算例分析

為了檢驗(yàn)新型限制器在光滑解附近對精度的影響，首先計(jì)算NACA0012翼型的亞聲速流場問題，設(shè)定計(jì)算狀態(tài)為：Ma∞=0.63，α=2°，對應(yīng)的物面壓強(qiáng)系數(shù)分布曲線和等壓線分布如圖2、3所示。在不使用限制器的條件下，這個(gè)算例依然能夠得到收斂的光滑解。從圖2中可以看出，加入分層限制器后，壓強(qiáng)系數(shù)分布基本沒有發(fā)生變化；從圖3中可以看出，使用限制器后，等壓線仍然是很光滑的。所以，可以認(rèn)為分層限制器有利于光滑解區(qū)域高精度的保持。

圖2 NACA0012翼型壓強(qiáng)系數(shù)分布（亞聲速）Fig.2 The distribution of pressure coefficient for NACA0012（subsonic）

圖3 NACA0012翼型等壓線分布（亞聲速）Fig.3 Lines of constant pressure for NACA0012（subsonic）

為了驗(yàn)證新型限制器在捕捉激波時(shí)的有效性，本文計(jì)算了兩個(gè)二維Riemann問題［10］，它的物理背景為：把［0，1］×［0，1］區(qū)域等分為四個(gè)彼此隔離的正方形子區(qū)域，每個(gè)區(qū)域上設(shè)定一個(gè)初始的流動(dòng)狀態(tài)，然后撤掉隔板，研究整個(gè)區(qū)域的流場變化規(guī)律。初始時(shí)刻算例1和算例2的密度ρ、速度u、v和壓強(qiáng)p的空間分布以及輸出流場數(shù)據(jù)時(shí)的計(jì)算時(shí)刻t見表1。相應(yīng)時(shí)刻密度ρ的分布情況如圖4、5所示。在新型限制器的輔助下，高階間斷Galerkin格式所計(jì)算出的激波位置以及形狀和文獻(xiàn)［10］的結(jié)果吻合得很好，激波干涉的細(xì)致結(jié)構(gòu)得以生動(dòng)的體現(xiàn)。

表1 二維Riemann問題初始狀態(tài)（ρ，u，v，p）及計(jì)算時(shí)刻tTable 1 Initial state（ρ，u，v，p）and final time for the 2D Riemann problem

接著，計(jì)算NACA0012和RAE2822翼型的跨聲速流場，無窮遠(yuǎn)處自由來流狀態(tài)分別為Ma∞=0.8，α=1.25°和Ma∞=0.729，α=2.31°，對應(yīng)的物面壓強(qiáng)系數(shù)分布曲線和等壓線分布如圖6～圖9所示，其中也給出了利用Shu限制器和HWENO限制器得到的對比曲線?？梢园l(fā)現(xiàn)，相對于Shu限制器，本文的限制器在激波附近更好的抑制了非物理振蕩；相對于HWENO限制器來說，本文的限制器對激波的識(shí)別更加敏感，對應(yīng)的曲線間斷部分更加陡峭。

圖4 二維Riemann問題密度等值線，t=0.23（算例1）Fig.4 Lines of constant density for 2D Riemann problem，t=0.23（case 1）

圖5 二維Riemann問題密度等值線，t=0.25（算例2）Fig.5 Lines of constant density for 2D Riemann problem，t=0.25（case 2）

圖6 NACA0012翼型壓強(qiáng)系數(shù)分布（跨聲速）Fig.6 The distribution of pressure coefficient for NACA0012（transonic）

圖7 NACA0012翼型等壓線分布（跨聲速）Fig.7 Lines of constant pressure for NACA0012（transonic）

圖8 RAE2822翼型壓強(qiáng)系數(shù)分布（跨聲速）Fig.8 The distribution of pressure coefficient for RAE2822（transonic）

圖9 RAE2822翼型等壓線分布（跨聲速）Fig.9 Lines of constant pressure for RAE2822（transonic）

然后，我們計(jì)算了一個(gè)前臺(tái)階問題，這是一個(gè)比較難于計(jì)算的非穩(wěn)態(tài)測試模型，因?yàn)榕_(tái)階角點(diǎn)會(huì)產(chǎn)生壓力奇性。相關(guān)幾何模型的數(shù)據(jù)詳見文獻(xiàn)［11］，計(jì)算區(qū)域上三角形單元數(shù)為17000，節(jié)點(diǎn)數(shù)為8751，超聲速入口條件為Ma∞=3，α=0°，計(jì)算時(shí)刻為t=4，其密度等值線和壓強(qiáng)等值線如圖10、11所示。與文獻(xiàn)［11］的結(jié)果相比較，可以發(fā)現(xiàn)，結(jié)合新型限制器的間斷Galerkin格式能夠很精確的計(jì)算出激波的反射位置和反射角。

圖10 前臺(tái)階問題密度等值線分布（t=4）Fig.10 Lines of constant density for the flow in the wind-tunnel with a front step（t=4）

圖11 前臺(tái)階問題壓強(qiáng)等值線分布（t=4）Fig.11 Lines of constant pressure for the flow in the wind-tunnel with a front step（t=4）

最后，我們計(jì)算一個(gè)超燃沖壓發(fā)動(dòng)機(jī)的內(nèi)部流場，這是一個(gè)定常問題，相關(guān)幾何模型數(shù)據(jù)詳見文獻(xiàn)［12］，計(jì)算區(qū)域上三角形單元數(shù)為12828，節(jié)點(diǎn)數(shù)為6785，超聲速入口條件為Ma∞=3，α=0°，出口處采用超聲速邊界條件，穩(wěn)定狀態(tài)下的密度等值線和壓強(qiáng)等值線如圖12、13所示。與文獻(xiàn)［12］的計(jì)算結(jié)果相比，可以發(fā)現(xiàn)，在新型限制器的幫助下，間斷Galerkin格式處理復(fù)雜流動(dòng)問題的能力得到了進(jìn)一步的展示。

圖12 超燃沖壓發(fā)動(dòng)機(jī)密度等值線分布Fig.12 Lines of constant density for the flow in the scramjet engines

圖13 超燃沖壓發(fā)動(dòng)機(jī)壓強(qiáng)等值線分布Fig.13 Lines of constant pressure for the flow in the scramjet engines

5 小結(jié)與展望

本文建立了模板坐標(biāo)系，溝通了非結(jié)構(gòu)網(wǎng)格單元與相鄰單元的聯(lián)系。進(jìn)一步的，在分析Gibbs現(xiàn)象的數(shù)值特征的基礎(chǔ)上，對各階方向?qū)?shù)進(jìn)行了分層限制，從而對各階廣義Fourier系數(shù)進(jìn)行了重構(gòu)。大量的數(shù)值結(jié)果表明，本文設(shè)計(jì)的新型限制器能夠有效的抑制間斷解附近產(chǎn)生的非物理振蕩，從而使計(jì)算過程得以順利進(jìn)行，同時(shí)，它保持了間斷Galerkin格式高精度逼近的優(yōu)越性，為其在工程中的廣泛應(yīng)用打下了堅(jiān)實(shí)的基礎(chǔ)。今后，可以在本文的基礎(chǔ)上，對 Gibbs現(xiàn)象初篩部分做更加深入的研究并提出新的判據(jù)，在限制器構(gòu)造部分研究更高精度的差分格式，對加權(quán)平均因子進(jìn)行新的改進(jìn)，從而進(jìn)一步提高間斷Galerkin格式處理復(fù)雜流動(dòng)問題的能力。

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

［2］ COCKBURN B，SHU C W.The Runge-Kutta discontinuous Galerkin method for conservation laws v：multidimensional systems［J］.J Comput Phys，1998，141：199-224.

［3］ HOTEIT H，ACKERER P，MOSE R，et al.New two-dimensional slope limiters for discontinuous Galerkin methods on arbitrary meshes［R］.Institut National de Recherche en Informatique et en Automatique（France），Rapport de Recherche，n 4491，July 2002.

［4］ TU S，ALIABADI S.A slope limiting procedure in discontinuous Galerkin finite element method for gasdynamics applications ［J］.Int.J.Numer.Anal.Model.，2005，2（2）：163-178.

［5］ QIU J，SHU CW.A Comparison of trouble cell indicators for Runge-Kutta discontinuous Galerkin methods using WENO limiters［J］.SIAM Journal on Scientific Computing，2005，27 （3）：995-1013.

［6］ LILIA KRIVODONOVA.Limiters for high-order discontinuous Galerkin methods［J］.Journal of Computational Physics，2007，226：879-896.

［7］ WANG L H，Xu X J.The mathematical foundations of the finite element method［M］.Beijing：Science Press，2004.（in Chinese）王烈衡，許學(xué)軍.有限元方法的數(shù)學(xué)基礎(chǔ)［M］.北京：科學(xué)出版社，2004.

［8］ CHENG Q X，ZHANG D Z，WEI G Q.Real variable function and functional analysis［M］.Beijing：Higher Education Press，2003.（in Chinese）程其襄，張奠宙，魏國強(qiáng).實(shí)變函數(shù)與泛函分析［M］.北京：高等教育出版社，2003.

［9］ JAN S HESTHAVEN，TIM WARBURTON.Nodal discontinuous Galerkin methods：algorithms，analysis，and applications ［M］.Springer，2008：135-145

［10］KURGANOV A，TADMOR E.Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers［J］.Numerical Methods Partial Differential Equations，2002，18：584-608

［11］WOODWARDP，COLELLA P.The numerical simulation of two-dimensional fluid flows with strong shocks［J］.J.Comput. Phys.，1984，54：115-173

［12］SHUANGZHANG Tu，SHAHROUZ ALIABADI.A slope limiting procedure in discontinuous Galerkin finite element method for gas dynamics applications［J］.International Journal of Numerical Analysis and Modeling，2005，2（2）：163-178.

An original designed limiter for discontinuous Galerkin method

GUO Yongheng，YANG Yong，LIU Yiou

（National Key Laboratory of Science and Technology on Aerodynamic Design and Research，Northwestern Polytechnical University，Xi′an 710072，China）

When the high-order discretization scheme is applied for solving the hyperbolic conservation laws，the Gibbs oscillations emerging around the discontinuities may cause the fields to take unphysical values，which could result in the undesired interruption of calculation process.This problem is especially severe for the discontinuous Galerkin scheme.Therefore，suitable limiting calculators are supposed to be designed in order to damp the appearance of unphysical oscillations.However，the existed limiters for the discontinuous Galerkin scheme originate mostly from the finite volume methods，thus for the unstructured grids they can only limit the spatial solutions of the solutions，namely the lower-order derivatives，rather than offer an general criterion for the higher order derivatives.Due to this background，the discontinuous finite element solutions are firstly processed into general Fourier expansion series for the realization of spectral decomposition in different frequency domains.Afterwards，a modal coordinate system is originally proposed，describing the changing patterns of the directional derivatives in different orders in a new form.Finally，after the combination of the information of the current mode and its adjacent mode，and the limitation of the directional derivatives in different orders for different levels，the unphysical solutions are successfully damped.Based on this method，the Euler equations are solved for the numerical simulations of the 2-D Riemann problems，the flow field of the wind-tunnel with a front step，the reflection of shockwaves in the inner field of scramjet engines as well as the subsonic and transonic flow surrounding airfoils.The obtained results demonstrate that the new limiter is relatively reliable and it can be generalized to high-order schemes.

discontinuous Galerkin scheme；limiting calculators；general Fourier expansion；mode coordinate system；directional derivatives.

V211.3

Adoi：10.7638／kqdlxxb-2012.0157

0258-1825（2014）04-0462-06

2012-10-16；

2013-03-05

郭永恒（1982-），男，西北工業(yè)大學(xué)博士研究生，主要從事理論與計(jì)算流體力學(xué)的研究.

郭永恒，楊永，劉逸鷗.一種基于間斷Galerkin格式的新型限制器設(shè)計(jì)［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2014，32（4）：462-467.

10.7638／kqdlxxb-2012.0157. GUO Y H，YANG Y，LIU Y O.An original designed limiter for discontinuous Galerkin method ［J］.ACTA Aerodynamica Sinica，2014，32（4）：462-467.