一種基于間斷Galerkin格式的新型限制器設計

郭永恒，楊 永，劉逸鷗

（西北工業(yè)大學 翼型葉柵空氣動力學國防科技重點實驗室，陜西西安 710072）

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

一種基于間斷Galerkin格式的新型限制器設計

郭永恒，楊 永，劉逸鷗

（西北工業(yè)大學 翼型葉柵空氣動力學國防科技重點實驗室，陜西西安 710072）

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

間斷 Galerkin格式；限制器；廣義Fourier展開；模板坐標系；方向導數

0 引 言

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

1 間斷Galerkin格式的建立

直角坐標系下二維Euler方程的守恒形式為

首先，把整個求解區(qū)域Ω 劃分成三角形非結構單元的集合，即其中 NE 為網格單元的個數；同時，在第k號單元上，設局部有限元解為

其中PHj是當前單元上的基函數，NP 為基函數的個數。根據加權余量法［7］的思想，令測試函數空間等于解函數空間，作如下內積，得

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

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

沿時間方向使用Runge-Kutta方法推進，就得到式（5）的數值解。

2 對問題單元的初篩

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

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

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

3 模板坐標系下的限制器設計

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

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

圖1 模板坐標系示意圖Fig.1 A diagram for the mode coordinate system

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

這樣，聯立式（9）、（10）就可以求出修正后的一次層系數。對于模板1，我們定義一個加權因子w1=β1（3），即第3號單元的一次相對振蕩強度?？梢钥闯?，當第3號單元上的振蕩相對較強時，這樣的定義方式模板1中重構的信息所占的加權比重相對更大一些。同理，可定 義于是最后的加權修正結果

為

類似的，我們可以進一步比較高階方向導數的解析結果和差分結果，建立如下三個限制方程

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

2，即第3號單元的二次相對振蕩強度。

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

4 數值算例分析

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

圖2 NACA0012翼型壓強系數分布（亞聲速）Fig.2 The distribution of pressure coefficient for NACA0012（subsonic）

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

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

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

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

圖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翼型壓強系數分布（跨聲速）Fig.6 The distribution of pressure coefficient for NACA0012（transonic）

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

圖8 RAE2822翼型壓強系數分布（跨聲速）Fig.8 The distribution of pressure coefficient for RAE2822（transonic）

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

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

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

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

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

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

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

5 小結與展望

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

［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）王烈衡，許學軍.有限元方法的數學基礎［M］.北京：科學出版社，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）程其襄，張奠宙，魏國強.實變函數與泛函分析［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è)大學博士研究生，主要從事理論與計算流體力學的研究.

郭永恒，楊永，劉逸鷗.一種基于間斷Galerkin格式的新型限制器設計［J］.空氣動力學學報，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.