守恒律方程的導(dǎo)數(shù)人工壓縮法

董海濤 劉福軍 韓 沖

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

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

1 導(dǎo)數(shù)人工壓縮法的思想和推導(dǎo)

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

1.1 Harten的人工壓縮法

Harten的原始二階TVD格式為

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

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

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

1.2 導(dǎo)數(shù)人工壓縮法

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

1.3 導(dǎo)數(shù)人工壓縮法守恒修正

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

其中

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

1.4 方程組的導(dǎo)數(shù)人工壓縮

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

方程組的導(dǎo)數(shù)人工壓縮格式為

2 算例驗(yàn)證

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

算例1 極值間斷算例

初值條件

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

圖1 單個方程間斷問題的數(shù)值解

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

初值條件

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

圖2 單個方程極值間斷相鄰的數(shù)值解

算例3 多極值算例

初值條件

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

圖3 單個方程多極值問題的數(shù)值解

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

算例4 振蕩激波管算例

文獻(xiàn) ［5］中構(gòu)造了初始值一端有密度脈動的激波管算例，初值條件為

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

圖5 導(dǎo)數(shù)人工壓縮計算結(jié)果與精確解的比較

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

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

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

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

2.4 格式混合參數(shù)θ的選取

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

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

3 結(jié)論

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

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