張來(lái)平，劉 偉，賀立新，鄧小剛，張涵信

(1.中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 綿陽(yáng) 621000;2.中國(guó)空氣動(dòng)力研究與發(fā)展中心，四川 綿陽(yáng) 621000)

0 引言

真實(shí)的流動(dòng)是很復(fù)雜的，流場(chǎng)中常常會(huì)有激波、旋渦、分離等現(xiàn)象。為了模擬這些復(fù)雜的流動(dòng)現(xiàn)象，高精度格式構(gòu)造方法越來(lái)越受到CFD學(xué)者們的關(guān)注。在眾多的高階精度計(jì)算方法中，間斷Galerkin有限元方法(DGM)備受關(guān)注。間斷有限元方法是最早由Reed和Hill為解決中子運(yùn)輸方程而提出的［1］。在此基礎(chǔ)上，特別是近年來(lái)，出現(xiàn)了種類豐富多樣的DGM方法。DGM方法保持了傳統(tǒng)有限元方法(FEM)和有限體積法(FVM)的優(yōu)點(diǎn)，融入了高分辨率有限差分方法(FDM)和有限體積方法中如數(shù)值通量、Riemann間斷分解、TVD和限制器等思想;同時(shí)它又是一種積分形式的計(jì)算方法，可以直接應(yīng)用于非結(jié)構(gòu)網(wǎng)格或混合網(wǎng)格，與其它數(shù)值方法相比，DGM具有突出的優(yōu)勢(shì)。

然而，高精度格式在激波附近會(huì)產(chǎn)生數(shù)值振蕩，引起非線性的不穩(wěn)定性，甚至出現(xiàn)非物理解，如負(fù)的壓力和密度。Harten在1983年首次提出了總變差減小(total variation diminishing，TVD)的概念［2］，為差分方法的高精度格式構(gòu)造開(kāi)拓了嶄新的方向，并通過(guò)限制器(limiter)技術(shù)，構(gòu)造出二階TVD格式。限制器的作用是引入非線性效應(yīng)。常見(jiàn)的限制過(guò)程是通過(guò)比較臨近單元解的性態(tài)，使得重構(gòu)的左右變量不至于“上沖”或“下沖”，以保證解的本質(zhì)無(wú)振蕩，我們把這種限制器稱為梯度限制器(slope limiter)，其本質(zhì)上是對(duì)一階導(dǎo)數(shù)進(jìn)行限制，TVD格式和NND［3］格式均屬此類，不同的是TVD格式是對(duì)原始變量或守恒變量進(jìn)行限制，而原始NND格式限制于通量。

這類梯度限制器的缺陷在于對(duì)光滑的極值點(diǎn)附近也會(huì)起作用，從而降低了解的一致高精度。間斷偵測(cè)器的作用就是偵測(cè)出激波所在的位置，從而只在“間斷區(qū)域”引入限制器，其它的“光滑區(qū)域”仍然采用線性格式，以改善數(shù)值解的模擬精度。Chi-Wang Shu等［4］通過(guò)TVB函數(shù)判斷出“問(wèn)題單元”(trouble cell);S.Adjerid等的研究［5］指出，DGM在單元的出流邊界具有很強(qiáng)的超收斂特性，并基于此特性構(gòu)造了一種間斷偵測(cè)器。

本文通過(guò)比較單元交界面的左右變量，構(gòu)造了一種新的間斷偵測(cè)器，并將該方法推廣應(yīng)用于一維和二維Euler方程的計(jì)算。首先利用間斷偵測(cè)器偵測(cè)出“間斷單元”(“問(wèn)題單元”)，然后對(duì)“問(wèn)題單元”進(jìn)一步引入梯度限制器，由于本文采用的是基于Taylor基的間斷Galerkin方法，梯度限制器能夠更直接地實(shí)現(xiàn)。針對(duì)一維問(wèn)題，我們采用 minmid限制器［3］，及其推廣的 Moment限制器［6］;二維情形采用魯棒性更好的Barth限制器［7］。數(shù)值計(jì)算結(jié)果表明本文的間斷偵測(cè)器能較好地捕捉“問(wèn)題單元”，從而盡可能減少限制器對(duì)計(jì)算精度的影響，同時(shí)保證在強(qiáng)間斷處沒(méi)有非物理振蕩。

1 間斷偵測(cè)器和DG計(jì)算方法

基函數(shù)的選取對(duì)于間斷Galerkin方法尤為重要，傳統(tǒng)的DGM對(duì)于三角形/四面體采用面積坐標(biāo)/體積坐標(biāo)構(gòu)造基函數(shù)，對(duì)于矩形單元?jiǎng)t采用雙線性插值構(gòu)造基函數(shù)［8］。其它形狀的單元?jiǎng)t需要坐標(biāo)變換至標(biāo)準(zhǔn)單元，格式中出現(xiàn)多種基函數(shù)，這使得限制器的構(gòu)造，隱式時(shí)間推進(jìn)，邊界條件處理等都存在諸多不便。Luo Hong基于Taylor展開(kāi)構(gòu)造了一組Taylor基［9］，其優(yōu)勢(shì)在于針對(duì)不同網(wǎng)格單元基函數(shù)的形式是一致的，能夠簡(jiǎn)單方便地應(yīng)用于帶有“懸空”節(jié)點(diǎn)的混合網(wǎng)格，而這種情況在多層次的Cartesian網(wǎng)格和自適應(yīng)網(wǎng)格中經(jīng)常出現(xiàn)。并且Taylor基具有先天性的層次性，更有利于p-multigrid策略的實(shí)施。本文中DGM的基函數(shù)采用的是Taylor基，限于篇幅，這里不再贅述，具體內(nèi)容可參考文獻(xiàn)［9］。這里重點(diǎn)介紹在DGM計(jì)算過(guò)程中，為了有效捕捉激波等強(qiáng)間斷的間斷偵測(cè)器，以及在一維和二維情況下的限制器。

1.1 間斷偵測(cè)器

DGM假設(shè)解是分段p次多項(xiàng)式，即使在光滑區(qū)域，通過(guò)左右單元插值得到的單元交界面的左右變量通常是不相等的，但它們滿足:

這里假設(shè)q(x)為精確解，xj為單元i的交界面上的第j個(gè)Gauss積分點(diǎn)，h為單元的特征長(zhǎng)度。則

也就是說(shuō)，在光滑區(qū)域，雖然假設(shè)單元交界面是間斷的，但左右變量的差別是p+1階小量。對(duì)于間斷區(qū)域或大梯度附近，左右變量的差別是很大的，通常為|qr－ql|∝O(h)。

為了方便書(shū)寫，我們記:變量q在xj的跳躍［q］j和平均分別為:

定義單元i的偵測(cè)因子

通過(guò)判斷Di的大小，以確定“間斷單元”，或稱“問(wèn)題單元(trouble cell)”:

θ為偵測(cè)閥值。(4)式中的Lp范數(shù)，在本文的算例中，均取L∞范數(shù)，即只要存在一個(gè)Gauss積分點(diǎn)的間斷偵測(cè)因子大于偵測(cè)閥值，則這個(gè)單元為“間斷單元”。本文沿用文獻(xiàn)［5］中的方法，分別選用密度和熵作為偵測(cè)的對(duì)象。

1.2 DGM格式中的限制器

本文的算例中，一維問(wèn)題采用Lilia Krivodonova構(gòu)造的Moment限制器［6］，二維問(wèn)題采用Barth和Jesperson的梯度限制器［7］。

Moment limiter的基本思想是判斷單元i的p階導(dǎo)數(shù)與若干相鄰單元的(p－1)階導(dǎo)數(shù)，從高階到低階依次對(duì)p階導(dǎo)數(shù)進(jìn)行限制，避免直接對(duì)梯度限制而帶來(lái)的過(guò)多耗散。在光滑區(qū)域，單元中心點(diǎn)的二階導(dǎo)數(shù)滿足(一維問(wèn)題):

其中ξ為局部坐標(biāo)。因此對(duì)二階導(dǎo)數(shù)采用如下限制:

同理對(duì)一階導(dǎo)數(shù)

如果二階導(dǎo)數(shù)被限制了，則進(jìn)一步限制一階導(dǎo)數(shù)，否則不對(duì)一階導(dǎo)數(shù)限制。這里的α1、α2為控制耗散的參數(shù)，它們的值越小，限制的強(qiáng)度越大，耗散也就越大，計(jì)算則更穩(wěn)定。數(shù)值實(shí)驗(yàn)表明取 α2∈［1，10］，α1∈［1，2］是合適的。

對(duì)于二維問(wèn)題，我們采用Barth限制器。Barth限制器對(duì)梯度進(jìn)行如下限制:

其中rc為單元中心點(diǎn)坐標(biāo)，rj為邊界上的Gauss積分點(diǎn)坐標(biāo)。限制系數(shù)α定義為:

限制器可以針對(duì)守恒變量、原始變量或特征變量進(jìn)行限制，一般來(lái)說(shuō)，限制特征變量魯棒性更好，而對(duì)守恒變量限制則更直接，不需要矩陣變換。本文中，限制的是守恒變量。對(duì)于高超聲速問(wèn)題，守恒變量型的限制器在數(shù)值通量的計(jì)算中，單元交界面左右的壓力可能出現(xiàn)負(fù)值，導(dǎo)致無(wú)法計(jì)算當(dāng)?shù)芈曀?。此時(shí)，我們用左右單元中心點(diǎn)的壓強(qiáng)值代替交界面的壓力。

2 一維Euler方程計(jì)算結(jié)果

2.1 算例1:Sod問(wèn)題

Sod問(wèn)題的初始條件由(11)式給出，在x=0.5處有初始的間斷。在計(jì)算區(qū)域［0，1］上，共用101個(gè)網(wǎng)格點(diǎn)，計(jì)算0.20個(gè)時(shí)間單位。圖1給出了p=1和p=2時(shí)利用DGM方法得到的計(jì)算結(jié)果，其中的“點(diǎn)標(biāo)+和*”分別代表密度和熵偵測(cè)器偵測(cè)到的“問(wèn)題單元”，可以看到，熵偵測(cè)器偵測(cè)到的“問(wèn)題單元”較少。圖2給出了密度和熵偵測(cè)器的偵測(cè)因子分布，可以看出，間斷偵測(cè)器能夠準(zhǔn)確地捕捉到激波所在位置。比起熵偵測(cè)器，密度偵測(cè)器的偵測(cè)敏感度更強(qiáng)，由于熵值在接觸間斷處仍然是連續(xù)的，熵偵測(cè)器對(duì)接觸間斷沒(méi)有響應(yīng)，而密度偵測(cè)器在接觸間斷處也有反應(yīng)(x≈0.23的位置)。

圖1 DGM方法求解Sod問(wèn)題的密度分布比較(T=0.20s，dx=1/100):(上)p=1;(下)p=2Fig.1 Density distribution for Sod problem with DGM，(Up)p=1;(Down)p=2

圖2 Sod問(wèn)題的密度和密度偵測(cè)因子分布(T=0.20s，dx=1/100)(上)p=1;(下)p=2Fig.2 The distribution of density and entropy indicators for Sod problem.(Up):p=1;(Down):p=2

2.2 算例2:Shu問(wèn)題

Shu問(wèn)題描述的是一個(gè)運(yùn)動(dòng)Mach數(shù)為3的激波和正弦波的相互干擾，從而誘發(fā)出高頻擾動(dòng)，被看作是激波與湍流相互作用的簡(jiǎn)化模型。該問(wèn)題流場(chǎng)中含有明顯的細(xì)致結(jié)構(gòu)，對(duì)于考察計(jì)算格式對(duì)流場(chǎng)細(xì)節(jié)和廣譜范圍內(nèi)的分辨率十分合適，一直是格式驗(yàn)證的標(biāo)準(zhǔn)算例。初始條件由式(12)給出，計(jì)算區(qū)域［－2，8］剖分為600個(gè)單元，時(shí)間演化至1.6個(gè)時(shí)間單位。圖3給出了利用三階DG方法(p=2)在不同的間斷偵測(cè)器作用下的計(jì)算結(jié)果，從圖3右邊的放大圖可以看出，不用間斷偵測(cè)器時(shí)限制器將“抹平”峰值，而熵偵測(cè)器能給出與精確解更為符合的計(jì)算結(jié)果，因?yàn)槠鋫蓽y(cè)到的“問(wèn)題單元”較少。

圖3 Shu問(wèn)題的間斷偵測(cè)器計(jì)算結(jié)果比較Fig.3 Shock detection for Shu problem.Left:Density indicator;Middle:Entropy indicator;

3 二維Euler方程計(jì)算結(jié)果

3.1 算例1:繞NACA0012翼型的跨聲速流動(dòng)

對(duì)于該外形，計(jì)算域采用非結(jié)構(gòu)的混合網(wǎng)格進(jìn)行剖分，即在壁面附近采用靈活性較好的三角形網(wǎng)格離散，而在外場(chǎng)采用非結(jié)構(gòu)的Cartesian網(wǎng)格，如圖4所示。網(wǎng)格單元總計(jì)5004個(gè)，其中三角形1100個(gè)，矩形單元2904個(gè)。計(jì)算方法采用了基于Taylor基的二階精度DG格式，并利用Barth限制器抑制激波附近的波動(dòng)。初始條件取來(lái)流馬赫數(shù)為0.85，攻角為1.25°。在此計(jì)算條件下，翼型的背部和腹部均會(huì)出現(xiàn)激波。圖5(a)給出了翼型壁面的壓力分布，與Shu等人的計(jì)算結(jié)果吻合較好。從圖5(b－c)給出的間斷偵測(cè)器判斷的“問(wèn)題單元”分布情況，可以看出密度偵測(cè)器和熵偵測(cè)器都能夠準(zhǔn)確地判斷激波所在的位置，對(duì)比熵偵測(cè)器，密度偵測(cè)器判斷的“問(wèn)題單元”更多，這和一維情形的結(jié)論一致。圖6顯示了兩種偵測(cè)器得到的流場(chǎng)Mach數(shù)等值線，可以看到，計(jì)算給出了清晰的激波結(jié)構(gòu)。

圖6 NACA0012翼型跨聲速繞流的Mach數(shù)等值線Fig.6 Mach number contours over NACA0012 airfoil

3.2 算例2:雙馬赫反射問(wèn)題

該問(wèn)題描述的是Mach數(shù)為10的強(qiáng)運(yùn)動(dòng)斜激波以與x軸方向呈60°角的方向入射，入射點(diǎn)在(1/6，0)，計(jì)算區(qū)域?。?，4］×［0，1］。激波波前靜止空氣的密度為1.4，壓力為1，波后按激波關(guān)系給定初始條件。下邊界在x＞1/6的區(qū)域給定壁面邊界條件，其它邊界區(qū)域按照激波運(yùn)動(dòng)所在的位置，分別給定波前或波后的值。時(shí)間推進(jìn)至t=0.2。利用Barth限制器和本文構(gòu)造的間斷偵測(cè)器，網(wǎng)格密度為480×120。采用基于Taylor基的二階DG方法(p=1)得到的計(jì)算結(jié)果如圖7所示，從圖中可以看出，引入間斷偵測(cè)器后既能較好地捕捉強(qiáng)激波結(jié)構(gòu)，又能給出比較精細(xì)的局部剪切結(jié)構(gòu)。圖8顯示了兩種間斷偵測(cè)器偵測(cè)到的“問(wèn)題單元”，同樣地，熵偵測(cè)器偵測(cè)到的“問(wèn)題單元”較少。圖9顯示了局部放大的流場(chǎng)結(jié)構(gòu)，可以看到，在沒(méi)有引入間斷偵測(cè)器時(shí)，接觸間斷處的細(xì)致結(jié)構(gòu)被限制器“抹平”，而采用間斷偵測(cè)器后，流場(chǎng)結(jié)構(gòu)的分辨率顯著提高。

4 結(jié)論

本文根據(jù)單元交界面左右變量的差別，提出了一種新的間斷偵測(cè)器構(gòu)造方法。該間斷偵測(cè)器的構(gòu)造原理簡(jiǎn)單，編程實(shí)現(xiàn)容易。針對(duì)一維和二維Euler方程，我們將此間斷偵測(cè)器用于間斷Galerkin格式的數(shù)值計(jì)算中，通過(guò)與適當(dāng)?shù)南拗破髋浜鲜褂?，可以有效抑制激波附近的?shù)值振蕩，并能保持在光滑區(qū)的計(jì)算精度。對(duì)一維的Sod問(wèn)題、Shu問(wèn)題，二維NACA0012翼型跨聲速繞流、雙馬赫反射問(wèn)題方程進(jìn)行了數(shù)值實(shí)驗(yàn)，計(jì)算結(jié)果證明了方法的有效性。

［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］HARTEN A.High resolution schemes for hyperbolic conservation laws［J］.Journal of Computational Physics，1983，49:357－393.

［3］張涵信.無(wú)波動(dòng)、無(wú)自由參數(shù)的耗散差分格式［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1986，6:143 －165.

［4］QIU J，SHU C W.Hermite WENO schemes and their application as limiters for Runge－Kutta discontinuous Galerkin method:one dimensional case［J］，Journal of Computational Physics，2003，193:115 －135.

［5］ADJERID S，DEVINE K，F(xiàn)LAHERTY J and KRIVODONOVA L.A posteriori error estimation for discontinuous Galerkin solutions of hyperbolic problems［J］.Computer methods in applied mechanics and engineering，2002，191:1097 －1112.

［6］KRIVODONOVA L.Limiter for high－order discontinuous Galerkin methods［J］.Journal of Computational Physics，2007，226:879－896.

［7］BARTH T，JESPERSON D.The design and application of upwind schemes on unstructured meshes［A］.In 27th Aerospace Sciences Meeting［C］.AIAA 89 －0036，Reno，Nevada，1989.

［8］賀立新，張來(lái)平，張涵信.間斷Galerkin有限元和有限體積混合計(jì)算方法研究［J］.力學(xué)學(xué)報(bào)，2007，39(1):15－21.

［9］LUO H，BAUM J D，LOHNER R.A discontinuous Galerkin method based on a Taylor basis for the compressible flows on arbitrary grids ［J］.Journal of Computational Physics，2008，227:8875－8893.