守恒律方程的導數人工壓縮法

董海濤 劉福軍 韓 沖

(北京航空航天大學 航空科學與工程學院，北京100191)

流體力學方程組解的間斷性在早期是一個主要難點，所有二階以上的經典格式 (如 Lax-Wendroff格式)在激波等間斷處會不可避免地出現偽振蕩.之后，由于TVD(Total Variation Diminishing)判據及限制器技術的出現，產生了一批無振蕩高精度格式，如TVD，ENO(Essentially Non-oscillatory)，NND(Non-Oscillatory and Non-free-parameter Dissipation)，無振蕩緊致格式、熵條件格式等，統(tǒng)稱為高分辨率格式，對激波的模擬取得了較好結果.許多復雜流動問題，流動參數變化劇烈，除可能存在激波間斷外，還有帶極值點的渦結構等，尤其對低速問題 (定常)沒有激波間斷，更多的 (高雷諾數)是極值點.而現有高分辨率格式若嚴格要求無振蕩，則都存在極值點精度退化問題，為了更好地模擬復雜流動問題，須考慮提高極值點處的精度.目前，對這個問題的研究還較少，本文提出的導數人工壓縮，就是為了解決20世紀80年代以來這些高階格式在極值點降階的問題，以提高對極值點的分辨率.本文先進行了導數人工壓縮方法的推導，并構造了人工壓縮和導數人工壓縮結合的混合格式，然后用單個標量方程的多極值間斷初值問題和Shu-Osher激波管算例進行了數值驗證，極值點處的精度得到較好的改善.

1 導數人工壓縮法的思想和推導

1977年文獻 ［1］提出人工壓縮方法并在1983年用于Harten-TVD［2］格式，以提高接觸間斷分辨率，本文受此啟發(fā)，獲得了導數人工壓縮方法構造思路.考慮守恒律方程初值問題

1.1 Harten的人工壓縮法

Harten的原始二階TVD格式為

其中絕對值號 Harten原格式為 Q(ν)［1］，因其增加數值粘性，并且已有嚴格熵條件方法［3］，本文忽略這種近似熵修正.

人工壓縮法既給原通量加一具有收縮特征場的小通量 (人工壓縮通量)，使特征線平行的線性 (退化)場變成非線性壓縮場，從而抑制數值耗散，達到提高 (接觸)間斷分辨率的目的(不同于提高精度).

Harten人工壓縮通量的構造方法為:通過引入人工壓縮通量增加間斷左端的特征速度，減小間斷右端的特征速度，從而把原特征場變成壓縮場［1］，將原格式的gj按如下公式替換，即可得到帶人工壓縮的Harten-TVD格式

1.2 導數人工壓縮法

值得注意的是這個格式是非守恒形式，只能用在原函數連續(xù)的地方，不能用于間斷處.即使線性方程也不能在間斷處用導數方程，因為間斷的導數為無窮，這種情況使用數值格式誤差非常大.最后綜合起來，數值方法用到兩種途徑推導的格式，需要研究一下如何匹配好.

1.3 導數人工壓縮法守恒修正

在計算中發(fā)現，對于非線性場有時甚至線性場，非守恒部分常常導致數值解的間斷有明顯相位誤差.如，用于激波管問題中時發(fā)現激波位置計算得不準確，非守恒通量權數θ越大，激波位置偏差越大.因此，本文又進行了守恒修正，用守恒格式的一階部分的主要項代替非守恒格式的一階部分的主要項，這樣導數方程的數值通量的非守恒部分是個小量，是近似守恒.將數值通量分為主要部分與小量部分

其中

式 (3)其余記號同式 (2)，守恒修正后格式為

1.4 方程組的導數人工壓縮

導數人工壓縮格式可按特征方向單個化的方法推廣于雙曲型守恒律方程組

方程組的導數人工壓縮格式為

2 算例驗證

單個方程算例取線性通量f(u)=u.方程組算例為一維完全氣體Euler方程

算例1 極值間斷算例

初值條件

計算時間為2，計算網格為100，CFL數為0.5，格式混合參數θ=0.12.分別采用Harten-TVD格式 (原格式)、帶人工壓縮的Harten-TVD格式 (人工壓縮)和本文新格式 (導數人工壓縮)進行計算，結果如圖1所示.可以看出，原格式在間斷處和極值點處抹光都很厲害，人工壓縮在極值點處和間斷處相對于原格式都有明顯改善.導數人工壓縮在間斷處同人工壓縮格式重合，在極值點處導數人工壓縮比人工壓縮格式更接近于精確解.

圖1 單個方程間斷問題的數值解

算例2 極值間斷相鄰算例［4］

初值條件

計算時間為0.5，計算網格為100，CFL數為0.3，格式混合參數θ=0.4.所采用的格式同上，在圖2中也可以看出人工壓縮和導數人工壓縮在間斷和極值點處的作用.

圖2 單個方程極值間斷相鄰的數值解

算例3 多極值算例

初值條件

計算時間為2，計算網格260，CFL數為0.3，仍然采用上述3種格式，格式混合參數θ取0.07.計算結果如圖3、圖4所示.此算例顯示原格式完全失去對脈動波的分辨能力，人工壓縮可以分辨出脈動但有較大誤差，而本文提出的導數人工壓縮則可以很小的誤差分辨脈動.

圖3 單個方程多極值問題的數值解

圖4 單個方程多極值問題的數值解局部放大圖

算例4 振蕩激波管算例

文獻 ［5］中構造了初始值一端有密度脈動的激波管算例，初值條件為

此算例沒有精確解，本文采用10000網格上的TVD格式的解作為精確解，用400網格上的計算結果進行比較.CFL數為0.8，計算時間為1.8，格式混合參數θ取0.04，計算結果如圖5所示.

圖5 導數人工壓縮計算結果與精確解的比較

圖6 Harten人工壓縮與精確解的比較

圖7 未做守恒修正時的計算結果與精確解的比較

為顯示新方法的優(yōu)點，以密度圖為例，對不同格式進行對比，同時顯示計算過程發(fā)現的問題.圖6顯示的是 Harten人工壓縮過壓縮造成的階梯現象，本文按照文獻［6］的方法，只在線性場中進行壓縮明顯改善了這個現象.圖7是格式 (3)計算的Shu-Osher算例密度圖，因格式的非守恒性造成的激波位置的偏移很明顯，可以看出，隨著非守恒數值通量 f～#的增大，激波位置的偏移越大來越大，這種偏移在使用守恒修正(式 (4)～式 (6))后消除.從圖8、圖9可以看出，采用人工壓縮格式的解在極值處比原始格式要好得多，導數人工壓縮格式在極值處的解又比人工壓縮的格式有明顯改善.導數人工壓縮格式在極值點不再受TVD或保單調限制，但仍受保凸性限制，故仍然嚴格滿足無振蕩條件.

圖8 三種格式的密度圖比較

2.4 格式混合參數θ的選取

由于本文構造的新格式壓縮性較大，作為導數人工壓縮權重的混合參數θ都很小，綜合幾個算例來看都在0.1附近，算例2中取值為0.4，注意到精確解中有兩處同時存在間斷和極值，為了抵消由此帶來的雙重數值耗散，所以取值有所偏大.而算例1和算例4中盡管也存在間斷和極值，但不在同一點出現，因此數值耗散沒有那么大，所以θ取值較小.關于混合格式的研究近年來已有不少學者進行了研究，如文獻［7－8］等.本文對混合參數的詳細研究會在后續(xù)的工作中進行.

圖9 三種格式密度分布局部放大圖比較

3 結論

比較上述數值實驗，可以看出人工壓縮方法比原TVD格式在間斷處和極值點處都提高了分辨率，本文發(fā)展的導數人工壓縮格式又比人工壓縮方法在極值點處提高了分辨率.當脈動多而網格不夠密時，一般格式可能會完全或部分喪失對脈動的分辨力，本文提出的導數人工壓縮法則能在相同情形較準確地分辯脈動，這種特性對于網格數不能太大的三維流體計算很有價值，特別是對脈動比較多的流場 (如湍流)更有價值.

References)

［1］Harten A.The artificial compression method for computation of shocks and contact discontinuities［J］.Communications on Pure and Applied Mathematics，1977，30(5):611－638

［2］Harten A.High resolution schemes for hypersonic conservation laws［J］.Journal of Computational Physics，1983，49(2):357－393

［3］Dong Haitao，Zhang Lidong，Lee Chun-Hian.High order discontinuities decomposition entropy condition schemes for Euler equations［J］.Computational Fluid Dynamics Journal，2002，10(4):448－457

［4］Harten A，Osher S.Uniformly high-order accurate non-oscillatory schemes［J］.SIAM J Numer Anal，1987，24(2):279 －309

［5］Shu C W，Osher S.Efficient implementation of essentially nonoscillatory shock-capturing schemesⅡ ［J］.Journal of Computational Physics，1989，83(1):32 －78

［6］Lie K A，Noelie S.On the artificial compression method for second-order non-oscillatory central difference schemes for systems of conservation laws［J］.SIAM J Sci Comput，2003，24(4):1157－1174

［7］Ren Y X，Liu M E，Zhang H X.A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws［J］.Journal of Computational Physics，2003，192(2):365－386

［8］Ziegler J L，Deiterding R，Shepherd J E，et al.An adaptive high-order hybrid scheme for compressive，viscous flows with detailed chemistry ［J］.Journal of Computational Physics，2011，230(20):7598－7630