陳二云，趙改平，楊愛玲，卓文濤

(1.上海理工大學 能源與動力工程學院，上海 200093;2.上海理工大學 醫(yī)療器械與食品學院，上海 200093)

氣動聲學作為氣動力學和聲學之間的一門交叉性學科，主要研究流動及其與物體相互作用產(chǎn)生噪聲的機理［1］。氣動聲學所涉及的問題具有非定常、小量級和較大的空間尺度分布，而且其脈動頻率也非常寬。與傳統(tǒng)的計算流體力學(CFD)方法相比，氣動聲學的數(shù)值方法不僅要具有較高的計算精度，還要具有盡可能低的色散、耗散誤差。因此，發(fā)展一種適合氣動聲學數(shù)值模擬的計算方法，進而來理解、預測并最終能夠控制噪聲的研究工作具有非常重要的理論意義和工程應用價值。

針對氣動噪聲的傳播特性，近幾年一些研究者對計算氣動聲學(CAA)的離散方法進行了大量的研究。如Tam和Webb［2］提出了保色散中心有限差分DRP格式;LeLe［3］提出了高階精度緊致類格式;Kim 等［4－5］在此基礎上進一步提出了優(yōu)化緊致差分格式以及優(yōu)化邊界緊致差分格式。此外，我國的傅德薰先生［6］也進行了相關的研究，提出了高精度迎風緊致型格式，并在試驗中得到了驗證。但這些方法都是基于均勻分布的笛卡爾網(wǎng)格條件下構造的對稱形式的有限差分格式，不適于計算具有復雜外形的氣動聲學問題。

80 年代末，Cockburn 和舒其望等［7－9］提出一種間斷伽遼金(DG)方法。該方法吸收了有限元方法和有限體積方法的優(yōu)點，具有一致高階精度、網(wǎng)格適應性強、結構守恒性和容易向高維推廣等特點，在計算流體力學領域被廣泛采用，有些研究者已開始將該方法應用到氣動聲學的計算。因此，建立低耗散、低色散的DG方法對氣動聲場進行數(shù)值模擬有著廣闊的應用前景。

按照基函數(shù)的不同構造方式，可以將DG方法分為兩類［10－11］:Nodal-DG 方法和 Modal-DG 方法。在本文的研究中，僅考慮Nodal-DG方法。通過選取合適的基函數(shù)和插值節(jié)點，質(zhì)量矩陣是對角化的，有利于質(zhì)量矩陣求逆和并行計算。

本文工作主要是基于線性雙曲方程，運用本征值方法分析了高階DG半離散格式的色散、耗散特性。結果發(fā)現(xiàn)，對于任意給定的m階多項式基函數(shù)，數(shù)值波解有m+1個值，但僅有一個值能夠表示對應微分方程的物理波傳播方式，其余的都是寄生波，且兩種波型的傳播方向相反。此外，通過與DRP格式和六階緊致格式進行比較發(fā)現(xiàn)，在相同的計算精度下，DG方法的有效求解波數(shù)范圍介于DRP格式和六階緊致格式之間。通過對高斯波形傳播的計算比較發(fā)現(xiàn)，在較少的網(wǎng)格數(shù)下，DG方法的計算結果可以與緊致格式的計算結果相比，但優(yōu)于DRP格式的計算結果，非常適合于氣動聲學的數(shù)值模擬，值得進一步深入研究。

1 數(shù)學模型

考慮一維線性雙曲方程:)其中x∈R，f=au，a為常數(shù)。為便于討論，設方程的解具有周期性，初始條件為:u(x，0)=exp(ikx)，k是波數(shù)。

顯然，方程(1)具有如下形式的非平凡解:

2 數(shù)值離散

圖1 計算域剖分Fig.1 Computational domain

首先將一維計算區(qū)域剖分成互不重疊的等間距網(wǎng)格單元，其中網(wǎng)格步長，如圖1所示。下標r與l分別表示每個單元的左右邊界。在每個單元Dn上，將近似解表示成:

由于在DG方法中，并不要求解在單元邊界上是連續(xù)的，故等式右邊數(shù)值通量函數(shù)f*不是唯一定義的。本文計算中取數(shù)值通量的表達形式為Lax－Friedrichs通量:

為了建立數(shù)值色散關系式，設局部近似解滿足:

由于方程的解具有周期性，所以，

將上述各表達式進行等參變換后代入方程(5)，得:

其中，L=kh/(m+1)，Ω =ωh/［a(m+1)］。

顯然，間斷有限元空間半離散格式的數(shù)值色散關系式為:L=Ω。其中Ω為復數(shù)，Ω=Ωr+iΩt，實部代表數(shù)值格式的色散特性，虛部代表數(shù)值格式的耗散特性，該值通過在復數(shù)域空間求解關于方程(6)的廣義本征值問題而獲得。

3 色散特性分析

圖2表示DG方法的色散、耗散特性曲線分布圖。本文中以m=1為例進行說明。其中實線代表原方程的精確色散、耗散特性分布曲線，符號線代表間斷有限元方法的數(shù)值色散、耗散特性分布曲線。從圖中可以看出，對于任意給定的m階多項式基函數(shù)，數(shù)值波解有m+1個值，但僅有一個值能夠表示對應微分方程物理波的傳播方式(如圓圈線)，其余的都是寄生波(如四邊形線)，兩種波型的傳播方向相反，且寄生波在有效求解波數(shù)范圍之外具有較大的衰減率。

圖2 色散、耗散特性分布曲線Fig.2 Dispersion and dissipation characteristics

圖3表示格式精度對色散、耗散特性影響分布曲線圖。其中實線代表原方程的精確色散、耗散特性分布曲線，符號線代表基函數(shù)取不同的m階多項式時對應的數(shù)值色散、耗散特性分布曲線。從圖中可以看出，當m=1時，即數(shù)值格式在空間上具有二階精度，間斷有限元方法的有效求解波數(shù)范圍約為0.62(以相對誤差小于0.5%為標準)，當m=3時，有效求解波數(shù)范圍約為1.09，當m=6時，有效求解波數(shù)范圍約為1.40。由此可以看出，隨著數(shù)值格式精度的提高，可以有效地降低DG方法的色散誤差和耗散誤差。

圖4表示相同計算精度下(六階)，DG方法與DRP格式和緊致格式的色散誤差對比圖，為了便于比較，對波數(shù)進行了歸一化處理。其中直線代表原方程精確的色散特性分布曲線，圓點、圓圈與三角形分別代表DG方法、DRP格式和緊致格式的數(shù)值色散特性曲線。從圖中可以看出，在相同的計算精度下，DG方法的有效求解波數(shù)范圍(約為1.25)介于DRP格式(約為1.15)和六階緊致格式(約為1.32)之間，非常適合于氣動聲學的數(shù)值模擬，值得深入研究。

4 計算結果

為了驗證本文方法求解波傳播問題的有效性，對文獻［12］中的經(jīng)典算例進行了求解，并與DRP格式和六階緊致格式計算結果進行了比較。

例1:控制方程是一維對流方程(1)，取a=1。

初始擾動為高斯波:u(x，0)=exp［－ln(2)(x2/4)］，計算區(qū)域為:x∈［－200，200］。

圖5表示m=1，無量綱時間t=40，網(wǎng)格步長不同時數(shù)值解與精確解的對比圖。其中實線表示精確解，圓圈表示數(shù)值解。從圖中可以看出，即使在低階精度情況下，隨著網(wǎng)格的逐步加密，數(shù)值解與精確解的吻合程度會迅速得以改善。

圖6表示DG方法、DRP格式和緊致格式在相同精度(六階)下的計算結果對比圖，計算時間t=40。其中實線表示精確解，圓圈表示數(shù)值解。從圖中可以看出，在相同的計算精度下，DG方法和六階緊致格式的計算結果與精確解的吻合程度是令人滿意的，但采用DRP格式獲得的結果具有明顯的色散特性，表明DG方法和緊致格式的有效求解波數(shù)超過了DRP格式的有效求解波數(shù)范圍，與上文的色散、耗散特性分析相一致。且值得需要指出的是在DG方法計算中，網(wǎng)格步長h=2，與其它兩種方法相比，僅采用了一半的網(wǎng)格數(shù)量。由此可以看出，DG方法具有較低的色散、耗散誤差，非常適合于波傳播問題的數(shù)值求解。

5 結論

運用本征值方法分析了高階Nodal-DG半離散格式的色散、耗散特性。結果發(fā)現(xiàn):

(1)對于任意給定的m階多項式基函數(shù)，數(shù)值波解在空間上有兩種傳播方式，其中之一能夠表示對應微分方程的物理波傳播方式，另一個是非物理振蕩，且兩種波型的傳播方向相反。

(2)通過與DRP格式和六階緊致格式進行比較發(fā)現(xiàn)，在相同的計算精度下，DG方法的有效求解波數(shù)范圍介于DRP格式和緊致格式之間，具有較低的色散、耗散誤差。

(3)通過對高斯波形傳播的計算比較發(fā)現(xiàn)，在相同的計算精度下，當用較少的網(wǎng)格數(shù)時，采用DG方法獲得的計算結果就可以與緊致格式的計算結果相比，且與精確解的吻合程度是令人滿意的，但DRP格式的計算結果有明顯的色散特性。

［1］ 王澤暉，羅柏華，劉宇陸.關于氣動聲學數(shù)值計算的方法與進展［J］.力學季刊，2003，24(2):219－226.

［2］ Tam K W，Webb J C.Dispersion-relation-preserving finite difference schemes for computational acoustics［J］.J.Comp.Phys.，1993，107:262－281.

［3］ Lele S K.Compact finite difference schemes with spectrallike resolution［J］.J.Comp.Phys.，1992，103:16－42.

［4］ Kim J W，Lee D J.Optimized compact finite difference schemes with maximum resolution［J］.AIAA Journal，1996，34:887－893.

［5］ Kim J W.Optimized Boundary Compact Finite Difference Schemes for Computational Aeroacoustics［J］.J.Comp.Phys.，2007，225:995－1019.

［6］ 馬延文，傅德薰.群速度直接控制四階迎風緊致格式［J］.中國科學，2001，36(6):554－561.

［7］ Cockburn B，Hou S C，Shu C W.TVB Runge-Kutta discontinuous Galerkin method for conservation lawsⅣ:the multidimensional case［J］.Math.Comp.，1990，54:545－581.

［8］ Cockburn B，SHu C W.TVB Runge-Kutta discontinuous Galerkin method for conservation laws V:Multidimensional systems［J］.J.Comp.Phys，1998，141:199－224.

［9］ Cockburn B，Shu C W.Runge-kutta discontinuous galerkin methods for convection-dominated problems［J］.J. Sci.Comp.，2001，16:173－261.

［10］ Hu F Q，Hussaini M Y，Rasetarinera P.An analysis of the discontinuous Galerkin method for wave propagation problems［J］.J.Comp.Phys，1999，151，921－946.

［11］ Hesthaven J S，Warburton T.Nodal discontinuous Galerkin methods:algorithms， analysis， and applications［M］.Springer，2007.

［12］ 楊愛玲，謝翠麗，陳康民.氣動聲學直接數(shù)值模擬的空間差分格式分析［J］.航空動力學報，2004，19(5):630－635.