武從海, 馬瑞軒,2, 王益民, 張樹海,*

(1.中國空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 綿陽 621000；2.中國空氣動(dòng)力研究與發(fā)展中心 氣動(dòng)噪聲控制重點(diǎn)實(shí)驗(yàn)室, 綿陽 621000)

0 引 言

氣動(dòng)聲學(xué)問題通常要求數(shù)值格式具有良好的多尺度分辨率。高精度緊致差分格式由于其良好的譜分辨率，在計(jì)算氣動(dòng)聲學(xué)中得到了廣泛應(yīng)用[2-3]。但是緊致格式在使用時(shí)一般需要求解全局方程組，這會導(dǎo)致計(jì)算效率的降低和并行處理上的困難。

優(yōu)化差分格式是另一類常用于氣動(dòng)聲學(xué)問題的數(shù)值格式[4]。優(yōu)化差分格式分為中心型和迎風(fēng)型兩種。中心型格式不含耗散，在實(shí)際應(yīng)用中需要增加濾波過程。而迎風(fēng)型優(yōu)化格式由于其內(nèi)在的耗散一般無需濾波。Tam認(rèn)為迎風(fēng)型優(yōu)化格式通常會存在一定程度的不穩(wěn)定性[5]。文獻(xiàn)[6]給出了一種穩(wěn)定性判據(jù)，發(fā)現(xiàn)了兩種迎風(fēng)型優(yōu)化格式存在不穩(wěn)定性，并采用這個(gè)判據(jù)作為約束條件分別構(gòu)造了相應(yīng)的穩(wěn)定優(yōu)化格式。但是，構(gòu)造迎風(fēng)優(yōu)化格式首先需要設(shè)計(jì)一個(gè)綜合考慮耗散誤差與色散誤差的目標(biāo)函數(shù)，而該目標(biāo)函數(shù)的設(shè)計(jì)存在較大的經(jīng)驗(yàn)性成分。

優(yōu)化格式通過犧牲收斂精度的階數(shù)而對某一特定范圍波數(shù)的波形(比如短波)的分辨率進(jìn)行優(yōu)化[4]。但是損失的精度階數(shù)會導(dǎo)致格式在波數(shù)0附近的波數(shù)誤差增大，因而格式在含大量低頻成分波形的遠(yuǎn)距離傳播問題中表現(xiàn)不佳。文獻(xiàn)[1]中，通過將優(yōu)化格式和最高階格式加權(quán)組合，得到的格式達(dá)到了最優(yōu)階收斂精度，并且對于無量綱波數(shù)不大于π/3的單色波的波數(shù)誤差幾乎為0。

由于該格式相比普通優(yōu)化格式具有更好的譜分辨率性質(zhì)，這里稱其為增強(qiáng)優(yōu)化差分格式(下文簡寫為增強(qiáng)優(yōu)化格式)。增強(qiáng)優(yōu)化格式由優(yōu)化格式和最高階格式非線性加權(quán)組合得到。本文將基于不同的線性優(yōu)化格式構(gòu)造增強(qiáng)優(yōu)化格式，并通過線性傳播方程的數(shù)值測試選取表現(xiàn)較優(yōu)的格式，最后針對聲散射問題進(jìn)行數(shù)值模擬與分析。

1 增強(qiáng)優(yōu)化格式

1.1 差分格式的波數(shù)誤差

本文討論的格式均為中心格式，采用的模板為7點(diǎn)對稱等距模板，網(wǎng)格間距為h。這種差分格式可寫為:

(1)

中心格式的系數(shù)滿足反對稱條件，即:

a0=0，a1=-a-1，a2=-a-2，a3=-a-3

對于波數(shù)為k的單色波,

f=eikx

其導(dǎo)函數(shù)

f′=ikf(x)

對于格式(1)，我們可以在節(jié)點(diǎn)處得到類似的公式:

(2)

考慮光滑函數(shù)f(x)，其傅立葉變換為:

F(f(x))=F(k)

則

F(f′(x))=ikF(k)

(3)

(4)

根據(jù)泰勒展開，差分格式的截?cái)嗾`差為:

=Mf(q+1)(x)hq+O(f(q+2)(x)hq+1)

(5)

其中M是由差分格式?jīng)Q定的常數(shù)。對上式做傅立葉變換得:

M(ik)q+1F(k)hq+O(kq+2hq+1)

(6)

(7)

由該式可知，高階精度的格式在波數(shù)0附近的波數(shù)誤差更小。

1.2 線性優(yōu)化差分格式

本小節(jié)討論7點(diǎn)對稱模板上的最高階格式及優(yōu)化格式。在該模板上的最高階顯式差分格式為六階。根據(jù)泰勒展開，六階格式的系數(shù)(用上標(biāo)S表示)需滿足:

(8)

求解得到的結(jié)果為:

(9)

這里考慮的優(yōu)化格式為四階精度的中心格式，且滿足在某無量綱波數(shù)為κ0∈[0,π]時(shí)波數(shù)誤差為0。根據(jù)泰勒展開和有效波數(shù)的計(jì)算公式可得到系數(shù)(用上標(biāo)O表示)滿足的方程組:

(10)

文獻(xiàn)[5]中的七點(diǎn)色散關(guān)系保持(Dispersion Relation Preserving,DRP)格式可由κ0取1.42151084298得到，下文用DRP4表示。文獻(xiàn)[1]中的優(yōu)化格式為無量綱波數(shù)κ0取π/3時(shí)的情況，對應(yīng)單位波長所含網(wǎng)格點(diǎn)數(shù)(Points Per Wavelength, PPW)[8]PPW=6。六階中心格式Cen6也可視為κ0趨于0的極限情況。

(11)

Opt2的系數(shù)為:

(12)

Opt3的系數(shù)為:

(13)

(a)Modified wavenumbers

1.3 增強(qiáng)優(yōu)化差分格式

根據(jù)下文中該格式的設(shè)計(jì)，其波數(shù)誤差在[0,κ0]中非常小，而對于更高波數(shù)的單色波，其表現(xiàn)與普通的優(yōu)化格式無異。而對于更一般的波形，如果將波形的每一小段(如四個(gè)節(jié)點(diǎn))當(dāng)作單色波，由于相位誤差小，各段的波形應(yīng)可以得到較好的保持，那么整體波形理應(yīng)可以得到保持。而實(shí)際的波形中，由于一些非常高波數(shù)成分的存在，這使得該格式的表現(xiàn)不如在單色波中突出。

構(gòu)造增強(qiáng)優(yōu)化格式時(shí)，首先仿照文獻(xiàn)[9]給出一種局部波數(shù)指示子：

IWj=[(fj-2-4fj-1+6fj-4fj+1+fj+2)2+

(-fj-2+2fj-1-2fj+1+fj+2)2]/

[(fj-1-2fj+fj+1)2+(-fj-1+fj+1)2+ε]

(14)

ε是一個(gè)小的正數(shù)，用來防止分母為0，這里取1×10-40。那么對于單色波f=sin(kx+φ)+C(C為常數(shù)項(xiàng))，

(15)

其中κ=kh為無量綱波數(shù)。

對于如何由IWj構(gòu)造單色波相位誤差較小的權(quán)值，最直接的方法是先求解無量綱波數(shù)κ，再構(gòu)造權(quán)系數(shù)。令:

(16)

從IWj得到λ0的過程需要計(jì)算反三角函數(shù)和幾個(gè)三角函數(shù)，計(jì)算量相對較大。因此采用一個(gè)簡單的公式來逼近：

λj=min(1,γj)

(17)

這里min表示兩者的較小值，其中

(18)

則λj在0到1之間，這使得增強(qiáng)優(yōu)化格式為兩個(gè)子格式的凸組合。

并且,

λj=aκ2+O(κ4)

其中系數(shù)C6和C4可由泰勒展開計(jì)算得到。

根據(jù)截?cái)嗾`差與波數(shù)誤差之間的關(guān)系式(7)，兩個(gè)子格式的波數(shù)誤差為：

那么增強(qiáng)優(yōu)化格式的波數(shù)誤差為：

=-C6(1-λj)κ7+C4λjκ5+O(κ9)

=(aC4-C6)κ7+O(κ9)

為了減小波數(shù)誤差，令首項(xiàng)系數(shù)為0，則:

(19)

對于參數(shù)b，我們要求當(dāng)無量綱波數(shù)為κ0時(shí)γj為1。那么:

則:

(20)

本文構(gòu)造的增強(qiáng)優(yōu)化格式為：

(21)

表1 幾種優(yōu)化格式對應(yīng)的參數(shù)a和b的值

圖2給出了幾種格式在0≤κ≤3π/4之間波數(shù)誤差的絕對值，橫坐標(biāo)為無量綱波數(shù)，間隔π/60，縱坐標(biāo)采用對數(shù)坐標(biāo)。其中Com6表示6階中心緊致格式[2]。從圖2中可以看出，在波數(shù)0附近DRP格式由于只有四階精度因而誤差最大，其次是Cen6，再次是Com6，而幾種增強(qiáng)優(yōu)化格式的誤差較小。在優(yōu)化格式及增強(qiáng)優(yōu)化格式各自的κ0處，根據(jù)其設(shè)計(jì)原理，波數(shù)誤差為機(jī)器0。而DRP和EOT由于相應(yīng)的κ0不在采樣點(diǎn)上，圖中該處僅出現(xiàn)小幅下探。對于幾種增強(qiáng)優(yōu)化格式，隨著相應(yīng)κ0的變大，格式在無量綱波數(shù)較小時(shí)的誤差隨之變大，而在無量綱波數(shù)較大時(shí)誤差則會變小。

注意上述波數(shù)誤差的結(jié)果是基于單色波得到的。對于更一般的波形，其各單色波成分的波數(shù)誤差與圖2所示并不一致。

圖2 幾個(gè)格式波數(shù)誤差絕對值的比較

2 收斂精度驗(yàn)證

文獻(xiàn)[1]中已經(jīng)對增強(qiáng)優(yōu)化格式的收斂精度進(jìn)行了分析，其結(jié)論對于本文中幾種增強(qiáng)優(yōu)化格式依然成立。增強(qiáng)優(yōu)化格式在流場光滑區(qū)域基本為6階精度，但是由于格式的非線性性質(zhì)，其收斂精度可能會在極值點(diǎn)附近降低。具體為：若求解問題是光滑的，且極值點(diǎn)處二階導(dǎo)數(shù)不為0，則增強(qiáng)優(yōu)化格式除極值點(diǎn)附近為5階外均為6階，則L1收斂精度為6階，L∞收斂精度為5階；若極值點(diǎn)處二階導(dǎo)數(shù)為0，則增強(qiáng)優(yōu)化格式L1收斂精度為5階，L∞收斂精度為4階。

下面通過數(shù)值實(shí)驗(yàn)驗(yàn)證收斂精度?？紤]如下線性對流方程的初值問題：

(22)

表2 增強(qiáng)優(yōu)化格式收斂精度測試

3 一維算例對比測試

本節(jié)通過一維算例來測試對比幾種增強(qiáng)優(yōu)化格式，同時(shí)參與計(jì)算的還有六階中心格式Cen6、優(yōu)化格式DRP4及六階中心緊致格式Com6。時(shí)間格式采用三階TVD龍格庫塔方法[10]。在實(shí)際計(jì)算中常常采用低耗散低色散龍格庫塔方法[11-12]。

3.1 正弦波傳播問題

采用對流方程進(jìn)行測試，對流速度設(shè)為1，其方程為：

ut+ux=0

(23)

初值為u0(x)=sin(4πx)，求解區(qū)域?yàn)閇0,1]，計(jì)算的CFL數(shù)均取0.1，PPW分別取6和12，推進(jìn)的時(shí)間分別為T=20和1000。圖3中給出了[0,0.5]間的計(jì)算結(jié)果。兩種情況下，六階中心格式Cen6和優(yōu)化格式DRP4的誤差都較大。當(dāng)PPW為6時(shí)，其它格式除EO3之外均得到了較好的結(jié)果，相比而言，EO2更為準(zhǔn)確。這是由于EO3僅對無量綱波數(shù)在[0,π/4]的單波誤差很小，PPW=6的正弦波的無量綱波數(shù)為π/3，超出了此范圍，從圖2中也可以看出此時(shí)EO3的波數(shù)誤差較大。當(dāng)正弦波的PPW為12時(shí)，四種增強(qiáng)優(yōu)化格式及緊致格式均取得了較好的結(jié)果。

(a)PPW=6, T=20

兩種情況下，EOT的結(jié)果大幅均優(yōu)于其組成成分格式DRP4。實(shí)際上，從圖2的結(jié)果可以看出，二者的波數(shù)誤差相差至少一個(gè)量級。這也說明了增強(qiáng)優(yōu)化格式在單色波問題中有著上佳的表現(xiàn)。

3.2 高斯波傳播問題

這里計(jì)算一個(gè)高斯波的傳播問題[13]，其控制方程與前一問題一致，而初始波形為：

(24)

其中h=1為空間步長，取時(shí)間步長Δt=0.4。計(jì)算終止時(shí)間為T=400和1000, 計(jì)算結(jié)果顯示在圖4中。整體而言，除優(yōu)化格式DRP4外，幾種格式均獲得了較好的結(jié)果，其中六階緊致格式Com6的結(jié)果最接近精確解。優(yōu)化格式DRP4在波后發(fā)生了較明顯的數(shù)值波動(dòng)，這是由于這些波動(dòng)的成分對應(yīng)的相速度偏大所致。圖4(b)中局部的放大圖中可以看到六階中心格式Cen6存在相對明顯的數(shù)值波動(dòng)，與DRP4格式相反，其原因?yàn)橄鄳?yīng)成分的相速度偏小。四種增強(qiáng)優(yōu)化格式中，其結(jié)果最好的是EO3，然后依次是EO2、EOT和EO1，其中EOT和EO1偏離精確解較大。這是由于該高斯波在步長為1時(shí)主要由長波(波數(shù)較小)構(gòu)成，而此時(shí)這兩種格式有一定的過度優(yōu)化，導(dǎo)致對PPW較大的波數(shù)誤差相對較大。

(a)T=400, Δt=0.4

圖4(c)給出了時(shí)間步長取Δt=0.1時(shí)T=1000的計(jì)算結(jié)果。相比圖4(b)，計(jì)算結(jié)果在高斯波兩端的振蕩幅度明顯加大。這說明當(dāng)Δt=0.4時(shí)，時(shí)間格式的耗散較大，消除了部分振蕩。該算例中，EO2和EO3取得了相對較好的結(jié)果。由于中心格式?jīng)]有耗散，實(shí)際計(jì)算中通常采用濾波方法[2,14]來消除這種連續(xù)流場中的非物理振蕩。

4 聲散射問題數(shù)值研究

當(dāng)聲波穿過旋渦時(shí)，會發(fā)生強(qiáng)烈的非線性散射現(xiàn)象，聲波的幅值和相位都會發(fā)生顯著的變化[15]。下面采用增強(qiáng)優(yōu)化差分格式及其他常用的氣動(dòng)聲學(xué)計(jì)算格式直接數(shù)值求解二維歐拉方程，計(jì)算平面聲波穿過等熵渦的散射問題。

以無窮遠(yuǎn)處聲速和密度對問題進(jìn)行無量綱化。背景流動(dòng)中二維等熵泰勒渦的渦心位于計(jì)算區(qū)域的中心(0,0)，其流動(dòng)參數(shù)如下[16]：

uθ=Mvrexp[(1-r2)/2]

(25)

ur=0

(26)

(27)

(28)

其中，uθ為背景流動(dòng)切向速度，ur為徑向速度，ρ和p分別為密度和壓強(qiáng)。最大馬赫數(shù)Mv=0.125，γ為氣體絕熱指數(shù)(空氣中一般取為1.4)。

如圖5所示，計(jì)算區(qū)域?yàn)閇-20,20]×[-20,20]。左邊界入射平面波形式為：

圖5 平面聲波穿過泰勒渦示意圖

(29)

其中，ρ′、u′、v′、p′分別為入射聲波的密度、x方向速度、y方向速度、壓強(qiáng)。聲壓振幅ε=1×10-5，波長λ分別為0.8和1。在計(jì)算區(qū)域右側(cè)和上下兩邊設(shè)置吸收層，減小反射聲波對計(jì)算的影響，滿足無反射邊界條件，具體形式參見文獻(xiàn)[17-18]。

根據(jù)一維算例的結(jié)果，我們采用綜合表現(xiàn)最好的EO2格式計(jì)算此問題。為了進(jìn)行比較，同時(shí)還采用了六階中心格式Cen6、優(yōu)化格式DRP4以及六階緊致格式Com6進(jìn)行計(jì)算。為了便于比較，計(jì)算中統(tǒng)一采用了8階緊致濾波格式[2]進(jìn)行濾波。由于該問題沒有理論精確解，這里采用時(shí)間步長為0.004，空間步長為0.02的六階緊致格式計(jì)算結(jié)果作為該問題的參考精確解(圖中用Ref.表示)。λ取0.8和1時(shí)計(jì)算終止時(shí)間分別取t=56和60。圖6給出了λ=1、t=60時(shí)參考精確解的聲壓分布?？梢钥吹铰暡ù┻^旋渦后會發(fā)生明顯的畸變。在旋渦正后方有明顯的相移，聲波幅值快速衰減，在局部形成“寂靜區(qū)”。下游±15°的方向形成兩條一級強(qiáng)噪聲干擾條帶，在一級條帶外側(cè)，形成若干個(gè)次級噪聲干擾條帶，強(qiáng)度依次減弱。由于泰勒渦逆時(shí)針旋轉(zhuǎn)，干擾條帶關(guān)于y軸不對稱，且旋渦上方干擾強(qiáng)度大于下方。

圖6 λ=1、t=60時(shí)聲壓云圖

該問題計(jì)算時(shí)間步長和空間步長分別取0.02和0.1。圖7給出了r=10的圓周上散射聲壓的指向性分布。首先我們?nèi)∪肷洳úㄩL為0.8。圖7(a)為幾種格式在各個(gè)方向的散射聲壓均方根。從圖7(a)中可以看出，在這種相對較粗的網(wǎng)格下，數(shù)值結(jié)果相比參考解有較明顯的誤差。其中增強(qiáng)優(yōu)化格式EO2與緊致格式Com6的結(jié)果相差很小，它們與參考解之間的差距較小。在聲波散射的主要區(qū)域，其幅值與參考解基本符合，在25°到70°之間，其幅值隨角度變化有一定幅度的振蕩。聲波散射影響較小的區(qū)域，其結(jié)果的幅值誤差略微偏大。六階中心格式的誤差較大，而優(yōu)化格式DRP4表現(xiàn)最差。然后我們將入射波波長設(shè)為1，圖7(b)給出了計(jì)算結(jié)果散射聲壓均方根的指向性分布。結(jié)果與圖7(a)相似，但是誤差變得更小了。說明對于該問題，隨著入射波的波長的增加，數(shù)值計(jì)算更容易獲得好的結(jié)果。

(a)λ=0.8

5 時(shí)間效率討論

本文數(shù)值模擬均采用Fortran程序在CPU為Intel Core i7 6700K@4GHz的臺式機(jī)上單線程計(jì)算。編譯優(yōu)化選項(xiàng)為“-O2”。

對于3.2節(jié)中高斯波傳播問題Δt=0.1、T=1000的情況，即圖4(c)對應(yīng)的算例，我們記錄了幾種格式該算例的計(jì)算時(shí)間。然后將空間間距兩次減半，分別記錄了其計(jì)算時(shí)間，三次計(jì)算的時(shí)間均列入表3中。從表3中可以看出增強(qiáng)優(yōu)化格式EO2僅需緊致格式Com6不到一半的時(shí)間。

表3 幾種格式計(jì)算高斯波傳播問題所需時(shí)間(Δt=0.1, T=1000, 單位: s)

表4中給出了聲散射問題幾種格式的計(jì)算時(shí)間。從表4中可以看出增強(qiáng)優(yōu)化格式EO2所需時(shí)間少于緊致格式Com6，但是該算例中EO2取得了與緊致格式相近的結(jié)果。并且，EO2是顯式格式，相比緊致格式更容易處理邊界條件并且更適于并行計(jì)算。

表4 幾種格式計(jì)算聲波散射問題所需時(shí)間(單位：s)

6 結(jié) 論

增強(qiáng)優(yōu)化格式為最高階格式與線性優(yōu)化格式的加權(quán)平均。加權(quán)系數(shù)通過極小化單色波的相位誤差得到，加權(quán)的過程保證了格式的最優(yōu)階收斂精度，并且提高了分辨率。該類格式對于單色波問題的誤差很小，而對于其它連續(xù)波形的傳播問題，也有較好表現(xiàn)。

本文給出了七點(diǎn)增強(qiáng)優(yōu)化格式的統(tǒng)一構(gòu)造方法，并比較了幾種增強(qiáng)優(yōu)化格式。線性傳播問題的數(shù)值實(shí)驗(yàn)表明增強(qiáng)優(yōu)化格式EO2綜合表現(xiàn)相對較好。針對平面聲波穿過泰勒渦的散射問題，采用了增強(qiáng)優(yōu)化格式EO2進(jìn)行數(shù)值模擬，并比較了幾種不同類型的格式。結(jié)果表明增強(qiáng)優(yōu)化格式EO2的結(jié)果與六階中心緊致格式相近，但是需要的計(jì)算時(shí)間更少。另外，由于該格式為顯式格式，因此相比緊致格式更容易處理邊界條件，并且更適于并行計(jì)算。

下一步將采用更多算例對該格式進(jìn)行進(jìn)一步測試，并將相應(yīng)算法加入國家數(shù)值風(fēng)洞(NNW)工程擬開發(fā)的氣動(dòng)噪聲數(shù)值模擬軟件。NNW工程是我國在2018年底啟動(dòng)的項(xiàng)目，旨在自主研發(fā)功能先進(jìn)、種類齊全的CFD軟件系統(tǒng)，相關(guān)介紹及詳細(xì)進(jìn)展參見文獻(xiàn)[19]。