董海濤 劉福軍 韓 沖

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

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

1.1 Harten的人工壓縮法

Harten的原始二階TVD格式為

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

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

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

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

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

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

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

其中

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

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

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

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

2 算例驗(yàn)證

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

算例1 極值間斷算例

初值條件

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

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

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

初值條件

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

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

算例3 多極值算例

初值條件

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

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

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

算例4 振蕩激波管算例

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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