自適應(yīng)的高精度中心-WENO混合格式

田 奇郭 元梁 賢李新亮

(1.北方民族大學(xué) 數(shù)學(xué)與信息科學(xué)學(xué)院,銀川 750021；2.中國科學(xué)院大學(xué) 工程科學(xué)學(xué)院,北京 100049；3.中國科學(xué)院力學(xué)研究所 高溫氣體動力學(xué)國家重點實驗室,北京 100190)

0 引 言

在數(shù)值模擬復(fù)雜流體流動問題中,特別是在研究非定常多尺度流動的細(xì)微結(jié)構(gòu)及其機(jī)理時,高精度格式比低精度格式更能準(zhǔn)確地捕捉各種小尺度量,特別像湍流、旋渦、氣動聲學(xué)等流動問題。因而,近年來高精度低耗散數(shù)值方法的研究越來越受到重視。對于可壓縮復(fù)雜流動的精細(xì)模擬,除了需要高精度、低耗散地分辨小尺度流動細(xì)節(jié),還需要對激波進(jìn)行捕捉,并盡量避免數(shù)值振蕩。1994年,Liu等對ENO格式進(jìn)行了改進(jìn)[1],提出了 WENO(Weighted Essential Non-Oscillatory)格式,WENO格式成為目前應(yīng)用最廣泛的激波捕捉格式。隨即,Jiang和Shu對該方法進(jìn)一步改進(jìn)[2],提出了新的光滑度量因子,使得 WENO格式的精度進(jìn)一步得到提高。數(shù)值試驗表明,Jiang和Shu的WENO格式具有較強(qiáng)的魯棒性,因而被廣泛地使用。此后,人們對WENO格式進(jìn)行著一系列測試、評估和改進(jìn)[3-7]。在國內(nèi),高精度激波捕捉格式的構(gòu)造也得到了廣泛的重視。鄧小剛等[8-9]構(gòu)造了高精度加權(quán)緊致格式WCNS和高精度耗散加權(quán)緊致格式DWCNS。張涵信等[10]基于三階ENN格式,引入加權(quán)函數(shù),構(gòu)造了一種具有四階或五階精度的WENN高階格式。傅德薰、馬延文先后提出了耗散比擬法[11]和群速度控制方法[12],以及三階和五階精度的迎風(fēng)緊致格式以及超緊致格式[13]。任玉新等[14]提出了色散最優(yōu)耗散可調(diào)的高分辨率線性格式。李新亮、何志偉等[15-16]提出非線性譜分析優(yōu)化的保單調(diào)格式和群速度控制方法。

一般而言,低耗散特性及強(qiáng)激波捕捉能力不容易兼顧。因而近年來發(fā)展出了混合格式(Hybrid schemes),即將兩種格式混合使用,在間斷或大梯度區(qū)采用強(qiáng)魯棒的激波捕捉格式,而在相對光滑區(qū)采取無耗散的中心格式或低耗散的線性格式。Pirozzoli[17]運用光滑識別器作為判據(jù),在光滑區(qū)使用緊致格式而在非光滑區(qū)使用WENO格式,從而較好實現(xiàn)了分辨小尺度波與捕捉激波的效果。Ren等使用了加權(quán)思想[5],將WENO格式與緊致格式混合使用,避免了兩種格式之間的“硬切換”,從而起到了更好的效果。

本文的主要工作就是基于WENO格式中的光滑因子,構(gòu)建一種新型的加權(quán)函數(shù),并在此基礎(chǔ)上將經(jīng)典的7階WENO格式和8階中心格式組合成混合格式,實現(xiàn)數(shù)值耗散的自適應(yīng)控制,構(gòu)造出一種自適應(yīng)低耗散的數(shù)值格式。通過對格式進(jìn)行非線性色散、耗散特性分析,對比研究格式在不同波數(shù)下的色散和耗散特性,并經(jīng)若干一維和二維算例數(shù)值驗證新格式對流場間斷及小尺度波的捕捉能力,以期為可壓縮湍流的數(shù)值模擬提供一種新的方法。

1 數(shù)值方法

1.1 數(shù)值格式構(gòu)造方法

使用Jiang等提出的WENO7格式[2]為通量插值形式,本文使用原始變量重構(gòu)形式。WENO7格式[18]為:

其中:為模板各給出的通量,ωk為各個模板對應(yīng)的非線性權(quán)值,相應(yīng)的表達(dá)式為:

根據(jù)Taylor級數(shù)展開式容易給出光滑區(qū)域的理想加權(quán)因子:

光滑度量因子βk的表達(dá)式為:

八階中心格式的數(shù)值通量表達(dá)式為:

半點格式的表達(dá)式為:

由于WENO格式的耗散相對較大,而中心格式無耗散,因此本文的目的是將兩種格式進(jìn)行混合,從而構(gòu)造出一種低耗散的數(shù)值格式?；谶@一思想,在此引入加權(quán)函數(shù)σ,格式變?yōu)?

為實現(xiàn)對該混合格式數(shù)值耗散的自適應(yīng)控制,本文建立起加權(quán)因子σ與光滑度量因子βk的關(guān)系。根據(jù)βk在光滑區(qū)時很小,而在間斷區(qū)相對較大這一特點,本文在此構(gòu)造σ與βk的關(guān)系式為:

在流場光滑區(qū)域βk很小,當(dāng)βk趨于0時,σ也趨于0,此時格式為近似八階中心格式,將具有很低的耗散；在流場非光滑區(qū)域βk值相對較大,此時σ趨于1,因此格式將會保留足夠的耗散來抑制非物理振蕩。這樣,該格式將會具有中心格式與WENO格式共同的特點。

1.2 數(shù)值格式誤差分析

采用Fourier精度分析法[19]分析該空間離散格式,所得到的修正波數(shù)的虛部Ki和實部Kr分別表示數(shù)值格式的色散和耗散效應(yīng),α為波數(shù)。從圖1中可以看出,本文H-WENO7-CD8格式的色散誤差要比經(jīng)典的WENO7更小,而WENO格式在高波數(shù)區(qū)存在較大的耗散誤差,其耗散也在高波數(shù)區(qū)也比WENO7要小,原因是H-WENO7-CD8格式是通過經(jīng)典的WENO7格式與中心格式加權(quán)組合而成,而中心格式無耗散。這就在一定程度上降低了新格式的耗散。

2 數(shù)值算例

2.1 激波-密度波干擾問題

該問題為一個馬赫數(shù)為3的運動激波和正弦密度波相互干擾,使用各空間離散格式模擬一維激波-密度脈動干涉算例[20],以測試數(shù)值方法對小尺度結(jié)構(gòu)的分辨能力。該控制方程為一維的Euler方程組,初始條件為:

圖1 兩種不同格式的色散(a)和耗散(b)特性Fig.1 Spectral properties of two different schemes

以x＝1處為左右側(cè)分界點,最終結(jié)果在區(qū)間x＝[0,10]上進(jìn)行計算,結(jié)束時間為t＝1.8。網(wǎng)格點數(shù)為N,由于該問題不存在解析解,因此選用N＝4001時的結(jié)果作為參考的精確解(Exact),采用Runge-Kutta方法進(jìn)行時間推進(jìn),通量分裂使用Steger-Warming分裂,使用局部特征分解。

圖2給出了WENO7格式和H-WENO7-CD8格式的計算結(jié)果?？梢钥闯?WENO格式和新格式均能較好地捕捉x＝7.4處的間斷以及x＝5.5左側(cè)的密度波動,但在該網(wǎng)格數(shù)下 WENO7得到的幅值明顯偏??；而H-WENO7-CD8格式不但能得到與精確解符合更好的波形,而且幅值也更接近精確解,因此新格式比WENO7對流場間斷及脈動具有更強(qiáng)的分辨率。

圖2 Shu-Osher問題,t＝1.8時刻的密度分布,N＝201Fig.2 Density distribution of the 1D shock/turbulence interaction at t＝1.8,N＝201

2.2 瑞利-泰勒不穩(wěn)定問題

Rayleigh-Taylor(RT)不穩(wěn)定性問題[19]是當(dāng)兩種密度不同的流體界面有微小的擾動,且由于某種原因從重流體到輕流體的方向產(chǎn)生加速度時,在這兩種流體的界面上就會出現(xiàn)不穩(wěn)定性。RT不穩(wěn)定性包含很多細(xì)致結(jié)構(gòu),以此用來測試數(shù)值格式分辨率。其計算條件設(shè)置如下[19]:計算域為 [0,0.25]×[0,1]；初始時刻,輕重流體的交界面位于y＝0.5處,密度ρ＝2的流體位于界面之下,密度ρ＝1的流體位于界面之上。初始場為:

源項的表達(dá)式為S＝(0,ρ,0,ρv)T。計算最終時間t＝1.95。

圖3分別給出了WENO7和H-WENO7-CD8兩種不同格式在網(wǎng)格60×240和240×960上計算得到的數(shù)值結(jié)果。從圖3中可以看出,在網(wǎng)格60×240上不同的格式捕捉到的流動結(jié)構(gòu)和相關(guān)文獻(xiàn)基本保持一致,密度大的流體沿中軸線流入密度小的流體,中軸線兩側(cè)形成了若干小尺度的渦結(jié)構(gòu),兩種格式均能保證質(zhì)量守恒,而新格式在此粗網(wǎng)格上能分辨出更多的小尺度結(jié)構(gòu)。網(wǎng)格加密到240×960之后,新格式捕捉到的小尺度渦結(jié)構(gòu)仍然要比WENO7要多,這就表明了在復(fù)雜結(jié)構(gòu)的數(shù)值計算中,H-WENO7-CD8格式擁有更高的數(shù)值結(jié)構(gòu)分辨率。

圖3 RT不穩(wěn)定性密度分布Fig.3 Rayleigh-Taylor instability：density profile

這是一個入射角為60°的馬赫數(shù)為10的激波的雙馬赫反射問題[21]。具體求解設(shè)置為:計算區(qū)域[0,0]×[4,1],下邊界y＝0,平板前緣位于x＝0.1667處,一直延伸到x＝4前,從x＝0到x＝0.1667給定波后條件；平板處采用波后壁面條件；上邊界各點的值由馬赫數(shù)為10的激波的精確運動來確定。計算最終時間為t＝0.2。

由于中心格式無耗散,而在雙馬赫反射的計算中需要數(shù)值格式有足夠的耗散,因此需要設(shè)定一個開關(guān)函數(shù),在間斷區(qū)時,要給予WENO格式足夠的權(quán)重,在此961×241的密網(wǎng)格的計算中,限制式(9)中的σ≥0.94,這樣就可以有效的防止發(fā)散。

圖4給出了在網(wǎng)格數(shù)相同情況下,采用不同格式所計算得到的密度等值線及其在接觸間斷和馬赫桿附近的放大圖(a:WENO7；b:H-WENO7-CD8),可以看出,兩種格式均能很好的捕捉到該問題基本的流動特征,如馬赫桿、激波和滑移線。另外通過比較兩幅局部放大圖可以清楚地看到,H-WENO7-CD8格式能分辨出更多的近壁面結(jié)構(gòu),更清晰的捕捉到了滑移線的卷曲現(xiàn)象和小尺度的渦,表明新格式對于小尺度結(jié)構(gòu)具有良好的分辨能力,具有較高的分辨率。

2.3 雙馬赫反射問題

圖4 密度分布,網(wǎng)格點數(shù)961×241,30條等值線,從ρ＝1.731到ρ＝20.92Fig.4 Density contours of the Double Mach Reflection problem,961×241 grid

3 結(jié) 論

本文通過新型加權(quán)函數(shù)將經(jīng)典的WENO格式和中心格式結(jié)合起來,利用光滑度量因子設(shè)定中心格式和WENO格式上的混合權(quán)重,實現(xiàn)了對格式數(shù)值耗散的自適應(yīng)控制,構(gòu)造出一種自適應(yīng)低耗散的數(shù)值格式。通過對一維激波密度脈動干涉問題和二維的瑞利-泰勒不穩(wěn)定性問題和雙馬赫問題進(jìn)行數(shù)值模擬,結(jié)果表明新的格式在繼承了原WENO格式良好的激波捕捉特性的同時,具有更低的耗散,對小尺度結(jié)構(gòu)更高的分辨率,并且有較好的穩(wěn)定性,這也為可壓縮湍流的高分辨率數(shù)值模擬提供了一種新的備選方法。此外,本文構(gòu)造格式是在均勻網(wǎng)格上推導(dǎo)的,這也是目前大多數(shù)高精度差分格式的做法。在實際應(yīng)用中,對于曲面非光滑網(wǎng)格,格式會產(chǎn)生幾何誘導(dǎo)誤差,幾何守恒問題對計算格式的精度和穩(wěn)定性都有較大的影響,將在今后的工作中考慮格式的幾何守恒律問題。

致謝：本項目得到國家自然科學(xué)基金(Nos.11472010,91441103,11372330,11472278)、國家重點研發(fā)計劃(2016YFA0401200)、科學(xué)挑戰(zhàn)專題項目(JCKY2016212A501)以及民用飛機(jī)專項科研(MJ-2015-F-028)的資助。感謝中國科學(xué)院力學(xué)研究所的傅德薰、馬延文研究員對本研究的幫助。感謝國家超級計算天津中心以及國家超級計算廣州中心提供計算機(jī)時。