基于映射函數(shù)的中心型三階格式

劉朋欣, 李 沁, 張涵信

(1. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室, 四川 綿陽 621000;2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 621000)

基于映射函數(shù)的中心型三階格式

劉朋欣1,2,*, 李 沁1,2, 張涵信1,2

(1. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室, 四川 綿陽 621000;2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 621000)

高階有限差分格式具有更小的數(shù)值耗散和色散誤差，可以更好地捕捉流場的精細(xì)結(jié)構(gòu)。在工程計算中，追求計算格式高分辨率的同時，也要兼顧其穩(wěn)定性。Henrick首先提出映射函數(shù)的方法，作用于WENO5格式的非線性權(quán)，使其滿足收斂階數(shù)要求并提高了格式的性能。不同于Henrick的構(gòu)造思路，李沁等提出了一種分段多項式形式的映射函數(shù)，改進(jìn)后的WENO5格式在保證穩(wěn)定性的前提下，擁有更好的分辨率。本文采用映射函數(shù)的思想，將一種五次分段多項式映射函數(shù)(Piecewise Polynomial Mapping Method, PPM)運用到中心型三階格式上，構(gòu)造出SWENO3-PPM5格式；將Henrick等提出的映射函數(shù)運用到中心型三階格式上，構(gòu)造出SWENO3-M格式，并與NND、WENO3格式的性能進(jìn)行比較。近似色散關(guān)系分析表明，SWENO3-PPM5和SWENO3-M的色散和耗散關(guān)系相比于WENO3得到明顯的提高。典型一維算例(Shu-osher問題)、二維算例(激波/層流邊界層干擾問題、激波/旋渦干擾問題、繞25°/55°尖雙錐流動問題)的計算結(jié)果表明，新格式具有比NND、WENO3更高的分辨能力，可以更好地捕捉到流場中的激波、接觸間斷面等復(fù)雜的流動現(xiàn)象，且由于分段多項式映射函數(shù)相比于Henrick提出的映射函數(shù)，可以使最終的非線性權(quán)值更加接近于格式的線性權(quán)，因此SWENO3-PPM5的分辨率要稍微高于SWENO3-M。

有限差分；映射函數(shù)；中心型三階格式；分辨率；穩(wěn)定性

0 引 言

隨著所研究的流動問題越來越復(fù)雜，研究和應(yīng)用對差分格式的計算精度也提出了更高的要求。相比于傳統(tǒng)的二階精度格式(如TVD格式)，高階精度格式具有更小的數(shù)值耗散和色散誤差，可以更好地刻畫流場的精細(xì)結(jié)構(gòu)。

高階有限差分格式是高階精度格式的重要組成部分。國外自從Harten提出ENO[1]和Liu、Osher和Chan等提出WENO[2]格式以來，眾多學(xué)者從不同的思路出發(fā)，提出各種形式的高階格式。其中，Jiang&Shu[3]提出了一種有效地實現(xiàn)WENO的方式，并將其用于有限差分，其五階形式(WENO5)成為最流行的高階格式之一，其非線性權(quán)求解過程及光滑度量因子(Indicator of Smoothness, IS)形式也逐漸成為后續(xù)研究的參考標(biāo)準(zhǔn)。盡管擁有很多的優(yōu)勢，但WENO5還存在一定的缺陷。Henrick等[4]首先指出，WENO5在極值點不能保持五階精度，他們分析了保持精度的充分必要條件，并提出了一種精心設(shè)計的映射函數(shù)對格式的非線性權(quán)進(jìn)行修正，相應(yīng)的格式稱為WENO-M。WENO-M可以保持精度一致，且分辨率有所提高。對非線性權(quán)進(jìn)行改進(jìn)并提出新格式的還有Borges[5]等提出的WENO-Z格式，Ha[6]等提出的WENO-NS格式，Ping Fan[7]等提出的WENO-η格式，Hui Feng[8]等提出的WENO-PM格式等。最近，Li等[9]從另一種構(gòu)造映射函數(shù)的思路出發(fā)，獨立地提出了一種分段多項式映射函數(shù)，對WENO5格式的非線性權(quán)進(jìn)行修正，提出了WENO-PPM格式。新格式在保持計算穩(wěn)定性的同時，進(jìn)一步提高了格式的分辨率。國內(nèi)張涵信[10]較早地提出了構(gòu)造差分格式應(yīng)該滿足的四個原則，提出了一種無振蕩、無自由參數(shù)耗散算法，即二階精度的NND格式。由于NND形式簡單、物理概念清楚、易于編程和計算量小，在國內(nèi)工程領(lǐng)域得到了廣泛的應(yīng)用。隨后，張涵信及其團(tuán)隊進(jìn)一步發(fā)展了一系列高精度計算格式，包括ENN[11]、WNND[12]、WENN[13]、中心型三階格式[14]等。

相比于基礎(chǔ)和機(jī)理研究所采用的高階格式，工程應(yīng)用的需求通常具有以下特點：1) 實際問題中網(wǎng)格結(jié)構(gòu)復(fù)雜，希望計算格式使用較小的網(wǎng)格模板；2) 由于實際問題流動的復(fù)雜性，要求格式具有類TVD格式的強(qiáng)魯棒性；3) 格式具有較好的耗散特性，能夠提升對旋渦等流動結(jié)構(gòu)的描述能力。已有計算顯示(如直升機(jī)旋翼的翼梢渦、渦流發(fā)生器產(chǎn)生流向渦的數(shù)值模擬等)，當(dāng)采用TVD格式對這些結(jié)構(gòu)進(jìn)行數(shù)值模擬時，結(jié)構(gòu)往往很快被耗散掉[14]。因此有必要進(jìn)一步發(fā)展符合工程需求的三階格式，在保證計算穩(wěn)定性的前提下，改善計算的精度和分辨能力，且能夠方便地推廣到有限體積程序。

本文采用映射函數(shù)的思想，將五次分段多項式映射函數(shù)應(yīng)用到中心型三階格式上，構(gòu)造了新型三階格式SWENO3-PPM5；并與Henrick等提出的映射函數(shù)作對比(SWENO3-M)。通過對典型算例的計算，驗證比較了新格式的性能。

1 數(shù)值計算方法

1.1 中心型三階格式

李沁等[14]通過分析格式的耗散關(guān)系，提出了一類耗散可調(diào)的中心型三階格式，提升了三階格式的性能。本文基于此格式，做進(jìn)一步的改進(jìn)。首先簡單介紹中心型三階格式。

以一維雙曲守恒律方程為例：

考慮守恒型差分格式，式(1)的數(shù)值求解可以寫成如下半離散化形式：

數(shù)值通量表示為：

本文計算中取p=2，ε=10-16。其中，ISk代表第k個模板的光滑度量因子。參考WENO系列格式非線性加權(quán)的基本框架ISk具體形式為：

由于中心型三階格式基于NND格式增加了一個下風(fēng)模板，格式的線性權(quán)系數(shù)中就引入了一個自由參數(shù)，記為n3rd。線性權(quán)系數(shù)見式(8)。在約化波數(shù)為π處，中心型三階線性格式的耗散為WENO3線性格式的1/n3rd。通過改變n3rd的值就可以調(diào)整中心型三階格式的耗散。具體分析可參考文獻(xiàn)[14]。本文的計算中取n3rd=10。

為了增強(qiáng)格式的穩(wěn)定性，在使用(7)式計算得到各個模板的光滑因子IS后，要對最下風(fēng)模板進(jìn)行一次修正：

1.2 分段多項式映射函數(shù)

通過使用映射函數(shù)來修正原來的非線性權(quán)ωk，可以使WENO5格式的精度保持五階。這個函數(shù)定義域為[0,1]，且滿足：

gk(0)=0,gk(1)=1

此函數(shù)將使用式(6)計算得到的非線性權(quán)ωk進(jìn)行一次映射，再將得到的值進(jìn)行歸一化作為最終的非線性權(quán)。圖1給出了Henrick型映射權(quán)函數(shù)與分段多項式型映射函數(shù)的分布。可以看出，相比于Henrick型映射權(quán)函數(shù)，采用五次分段多項式映射函數(shù)方法使非線性權(quán)在ω=Ck處的分布更加平坦，即更加接近于格式的理想線性權(quán)值，這將有利于提高格式分辨精細(xì)結(jié)構(gòu)的能力。

本文將分段多項式映射函數(shù)(式(14))運用于中心型三階格式，相應(yīng)的格式記為SWENO3-PPM5。并與采用Henrick等提出的映射函數(shù)的三階格式作比較，記為SWENO3-M格式。

1.3 不同格式的近似色散關(guān)系(ADR)對比

圖2給出了三種非線性格式的近似色散關(guān)系[15](Approximate Dispersion Relations, ADR)對比，可以看到：采用映射函數(shù)方法的中心型三階格式的色散和耗散關(guān)系得到明顯的提高，且SWENO3-PPM5的表現(xiàn)略優(yōu)于SWENO3-M。

2 數(shù)值計算結(jié)果與分析

本節(jié)計算了典型的一維和二維算例來檢測新型三階格式的性能，對比了SWENO3-PPM5、SWENO3-M、WENO3、NND四種格式。所有的算例均使用Steger-Warming分裂方式，并在計算中使用特征投影。時間推進(jìn)方面，對于非定常問題使用具有TVD性質(zhì)的三階Runge-Kutta格式，對于定常問題則使用LU-SGS格式。

2.1 Shu-Osher問題

該算例描述了一維情形下馬赫數(shù)為3.0的右行激波和密度波發(fā)生干擾作用，其目的在于分辨主激波波后的系列小尺度結(jié)構(gòu)，通常用于驗證格式的分辨率。初始條件為：

(ρ,u,p)=

計算采用801個網(wǎng)格點的均勻網(wǎng)格，結(jié)果如圖3和圖4所示。圖中把使用WENO5在1601個網(wǎng)格點上的計算結(jié)果作為準(zhǔn)確解。從計算結(jié)果來看，由于中心型三階格式的耗散低于WENO3格式，且使用映射函數(shù)可以使格式的權(quán)值更加接近于線性權(quán)，因此SWENO3-PPM5和SWENO3-M的分辨率要高于WENO3格式和NND格式；且由于PPM5映射函數(shù)相比于Henrick提出的函數(shù)，可以使非線性權(quán)在ω=Ck處的分布更加平坦，即更加接近于線性權(quán)，因此SWENO3-PPM5的分辨率也要略高于SWENO3-M。

2.2 激波/層流邊界層干擾問題

激波/層流邊界層干擾流動包含豐富的流動結(jié)構(gòu)，常作為檢驗格式性能的算例之一。來流狀態(tài)為M∞=2.0，T∞=293K，Re∞=2.96×105，Pr=0.72，β=32.585°，β為激波入射角；壁面為黏性絕熱條件；進(jìn)出口皆為超聲速邊界。計算域為：0.0≤x≤2.02，0.0≤y≤1.30，網(wǎng)格為nx×ny=102×121，x方向均勻分布，y方向在壁面附近進(jìn)行網(wǎng)格加密。

為了研究這幾種格式分辨膨脹波系的能力，本文取y=0.3045處的壓力分布(即圖5中紅線位置)，并

與WENO5的計算結(jié)果進(jìn)行比較，如圖6所示。可以看出，在x=1.4附近，NND格式完全不能分辨該區(qū)域；WENO3格式亦不能有效分辨該膨脹波系；雖然與模板較寬的WENO5格式的結(jié)果有一定差距，但是SWENO3-PPM5和SWENO3-M格式也能夠較好分辨到此區(qū)域，分辨率明顯優(yōu)于NND格式和WENO3格式。

2.3 激波/旋渦干擾問題

圖7給出了四種格式計算得到的密度等值線。可以看到，SWENO3-PPM5和SWENO3-M格式分辨出的流場結(jié)構(gòu)明顯比WENO3和NND更加豐富。由于初始流場含有激波間斷，為了適應(yīng)Euler方程的解算器，在x=-3.6處會形成垂直方向的二次干擾波。NND和WENO3格式由于耗散較大，并沒有較好地捕捉到此干擾波；而經(jīng)過映射后的中心型三階格式則較為清晰地分辨出此處的干擾波。同樣還有圖中方框區(qū)域內(nèi)的流場結(jié)構(gòu)，可以看出，SWENO3-PPM5分辨出的結(jié)構(gòu)較SWENO3-M更加細(xì)致。

2.4 繞25°/55°尖雙錐流動問題

繞尖雙錐流動是一類典型的高超聲速繞流。本文對繞25°/55°尖雙錐流動進(jìn)行模擬，來流條件為M∞=9.59，T∞=185.6K，Tw=293.3K，Re=139436m-1。網(wǎng)格為nstream×nnormal=262×148。計算結(jié)果與Gnoffo[17]給出的實驗和數(shù)值模擬結(jié)果進(jìn)行對比。Gnoffo數(shù)值模擬時使用的網(wǎng)格是512×256。

圖8給出了幾種格式計算得到的密度梯度等值線和流線分布。從密度梯度等值線可以看出，幾種激波捕捉格式基本上都捕捉到了流場中的激波、接觸間斷面、超聲速射流區(qū)、分離區(qū)等復(fù)雜的流動現(xiàn)象，但是SWENO3-PPM5和SWENO3-M分辨出的結(jié)構(gòu)要更加清晰一些，尤其是在分離區(qū)及超聲速射流區(qū)處；且在分離區(qū)，SWENO3-PPM5分辨出的旋渦數(shù)最多，NND格式最少。

圖9和圖10比較了幾種格式計算的壁面壓力和熱流分布。可以看出，SWENO3-PPM5的計算出的分離渦較大，且壁面壓力和熱流分布更加接近于Gnoffo的計算和實驗結(jié)果，顯示了SWENO3-PPM5較好的熱流預(yù)測能力。

3 結(jié) 論

本文采用映射函數(shù)的思想，將一種五次分段多項式映射函數(shù)運用到中心型三階格式上，構(gòu)造了SWENO3-PPM5格式，同時與Henrick等提出的映射函數(shù)作對比(SWENO3-M格式)。Shu-Osher算例、激波/層流邊界層干擾、激波/旋渦干擾以及繞25°/55°尖雙錐流動的計算表明，使用映射函數(shù)的新三階格式具有比NND、WENO3更高的分辨能力，可以有效地捕捉到流場中的激波、接觸間斷面等復(fù)雜的流動現(xiàn)象。且由于較Henrick的映射函數(shù)，分段多項式映射函數(shù)可以使最終的非線性權(quán)更加接近于格式的線性權(quán)，因此SWENO3-PPM5的分辨能力也要稍微高于SWENO3-M。

[1]Harten A, Enquist B, Osher S, et al. Uniformly high order essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1987, 71: 231-303.

[2]Liu X, Osher S, Chan T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994, 115: 200-212.

[3]Jiang G, Shu Qiwang. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126: 202-228.

[4]Henrick A K, Aslam T D, Powers J M. Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points[J]. Journal of Computational Physics, 2005, 207: 542-567.

[5]Borges R, Carmona M, Costa B, et al. An improved weighted essentially non-oscillatory scheme for hyperbolilc conservation laws[J]. Journal of Computational Physics, 2008, 227: 3191-3211.

[6]Ha Y, Kim C H, Lee Y J, et al. An improved weighted essentially non-oscillatory scheme with a new smoothness indicator[J]. Journal of Computational Physics, 2013, 232: 68-86.

[7]Ping Fan, Yiqing Shen, Naolin Tian, et al. A new smoothness indicator for improving the weighted essentially non-oscillatory scheme[J]. Journal of Computational Physics, 2014, 269: 329-354.

[8]Hui Feng, Fuxing Hu, Rong Wang. A new mapped weighted

essentially non-oscillatory scheme[J]. Journal of Scientific Computing, 2012, 51: 449-473.

[9]Li Q, Liu P X, Zhang H X. Piecewise polynomial mapping method and corresponding WENO scheme with improved resolution[J]. Communication in Computational Physics, 2015, 18(5): 1417-1444.

[10]Zhang H X. Non-oscillatory and non-free-parameter dissipation difference scheme[J]. Acta Aerodynamica Sinica, 1988, 6(2): 143-165.(in Chinese)張涵信. 無波動、無自由參數(shù)的耗散差分格式[J]. 空氣動力學(xué)學(xué)報, 1988, 6(2): 143-165.

[11]Zhang H X, He G H, Zhang L. Some important problems for high order accurate difference scheme solving gas dynamic equations[J]. Acta Aerodynamica Sinica, 1993, 11(4): 341-356.(in Chinese)張涵信, 賀國宏, 張雷, 等. 高精度求解氣動方程的幾個問題[J]. 空氣動力學(xué)學(xué)報, 1993, 11(4): 341-356.

[12]Liu W. Nonlinear dynamics analysis for mechanism of slender wing rock and study of numerical simulation method[D]. National University of Defense Technology, 2004.(in Chinese)劉偉. 細(xì)長機(jī)翼搖滾機(jī)理的非線性動力學(xué)分析及數(shù)值模擬方法研究[D]. 湖南: 國防科學(xué)技術(shù)大學(xué), 2004.

[13]Li Q, Zhang H X, Gao S C. Numerical simulations on supersonic shear layer flow[J]. Acta Aerodynamica Sinica, 2000, 18(S): 67-77. (in Chinese)李沁, 張涵信, 高樹椿. 關(guān)于超聲速剪切流動的數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報, 2000, 18(增刊): 67-77.

[14]Li Q, Sun D, Zheng Y K, et al. On a class of center-type third order difference scheme orienting to engineering utilizations[J]. Acta Aerodynamica Sinica, 2013, 31(4): 466-472. (in Chinese)李沁, 孫東, 鄭永康, 等. 一類中心型三階格式及應(yīng)用[J]. 空氣動力學(xué)學(xué)報, 2013, 31(4):466-472.

[15]Sergio Pirozzoli. On the spectral properties of shock-capturing schemes[J]. Journal of Computational Physics, 2006, 219: 489-497.

[16]Inoue O, Hattori Y. Sound generation by shock vortex interactions[J]. Journal of Fluid Mechanics, 1990, 380: 81-116.

[17]Gnoffo P. CFD validation studies for hypersonic flow prediction[R]. AIAA 2001-1025.

A kind of third order central scheme based on mapping functions

Liu Pengxin1,2,*, Li Qin1,2, Zhang Hanxin1,2

(1.StateKeyLaboratoryofAerodynamicsofChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China;2.ComputationalAerodynamicsInstituteofChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China)

High-order finite difference schemes have low numerical dissipation and dispersion error. They can capture flow structures in detail. In engineering applications, robustness should be considered, while high resolution is pursued. The method of mapping function was first proposed by Henrick to adjust nonlinear weights in [0,1] for the fifth-order WENO scheme. With the help of this method, the requirement of convergence order is satisfied, and the performance of the scheme is improved. Different from Henrick’s method, a concept of piecewise polynomial function is proposed by authors in this paper, and the corresponding improved WENO schemes have better resolution with the preservation of robustness. In this study, a fifth-order piecewise polynomial function (PPM5) is applied to a third order central scheme, which is referred as the SWENO3-PPM5. Henrick's mapping function is also applied to the third order central scheme, which is referred as the SWENO3-M. The performance of the SWENO3-PPM5 and the SWENO3-M are compared with that of the NND and the WENO3. The analyses of approximate dispersion relations (ADR) indicate that the dissipation and dispersion relations of SWENO3-PPM5 and SWENO3-M have an improvement compared with those of the WENO3. The performance of the new schemes is evaluated by testing typical problems, i.e. one-dimensional Shu-osher problem, two-dimensional ones including shock-boundary layer interaction problem, shock/vortex interaction problem, and hypersonic flows of the sharp double cone. The results indicate that the new schemes have better resolution than the WENO3 and the NND and can present a better description on complex flow structures, such as shock waves and contact discontinuity. Because the PPM5 can make the final weight closer to the idea weight of the scheme, the SWENO3-PPM5 possesses a little better resolution than the SWENO3-M.

finite difference scheme; mapping functions; third order central scheme; resolution; robustness

0258-1825(2017)01-0071-07

2015-09-01;

2016-04-13

國家重點基礎(chǔ)研究發(fā)展計劃(2014CB744100)，國家自然科學(xué)基金重大研究計劃(91541105)

劉朋欣*(1990-)，男，河南商丘人，碩士，主要從事計算流體力學(xué)及爆震燃燒研究. E-mail：liupengxin2008@163.com

劉朋欣, 李沁, 張涵信. 基于映射函數(shù)的中心型三階格式[J]. 空氣動力學(xué)學(xué)報, 2017, 35(1): 71-77.

10.7638/kqdlxxb-2015.0172 Liu P X, Li Q, Zhang H X. A kind of third order central scheme based on mapping functions[J]. Acta Aerodynamica Sinica, 2017, 35(1): 71-77.

V211.3

A doi: 10.7638/kqdlxxb-2015.0172