王 鑫， 余培汛， 楊 海， 潘 凱

(中國飛機(jī)強(qiáng)度研究所 航空聲學(xué)與動強(qiáng)度航空科技重點(diǎn)實(shí)驗(yàn)室， 西安 710065)

抑制短波的有效手段主要有兩種：① 迎風(fēng)格式，通過固有的耗散性來抑制短波，該方法的不足之處在于無法識別短波長波范圍，將兩者一起耗散，這樣影響到計(jì)算結(jié)果的準(zhǔn)確度；② 在對稱模板的差分方程中添加人工黏性項(xiàng)，所選擇的人工黏性項(xiàng)必須在有效抑制短波的同時(shí)，保證對長波的影響很小。Tam在1993年提出了人工選擇耗散法[3]，通過對人工黏性耗散項(xiàng)的系數(shù)選擇，使得其耗散項(xiàng)只對短波起作用，而對長波幾乎沒有影響，用于氣動聲學(xué)計(jì)算時(shí)取得了很好的結(jié)果。但是該方法只限于均勻?qū)ΨQ網(wǎng)格，尚未在非均勻網(wǎng)格中進(jìn)行拓展。Engquist[4]提出用非線性濾波法來捕獲激波的方法，它在每一時(shí)間步用一個(gè)非線性過濾因子對數(shù)值解進(jìn)行處理,只作用在部分網(wǎng)格點(diǎn)上。王澤暉等[5]在Bjrn提出的非線性濾波法基礎(chǔ)上，根據(jù)一個(gè)檢驗(yàn)函數(shù)，分析了非線性濾波效率與波長的關(guān)系，并引入了具有TVD性質(zhì)的非線性濾波法，通過數(shù)值算例驗(yàn)證了DRP格式與非線性濾波法結(jié)合的可行性，表明該方法對于在網(wǎng)格交界面、初始值間斷、非線性聲波傳播等計(jì)算氣動聲學(xué)常見情況下產(chǎn)生的短波有很好的抑制作用。該方法對短波具有很好的抑制作用，但同時(shí)對長波的耗散也相對較大，影響數(shù)值的計(jì)算精度。

本文在直角坐標(biāo)人工選擇耗散法的基礎(chǔ)上，將其轉(zhuǎn)換成曲線坐標(biāo)系下的表達(dá)形式，并在多塊網(wǎng)格交界面處、交界點(diǎn)、出口邊界等處按照當(dāng)?shù)鼐W(wǎng)格雷諾數(shù)的倒數(shù)進(jìn)行正態(tài)分布，并通過NASA標(biāo)模算例進(jìn)行了該方法的驗(yàn)證分析。

1 數(shù)值模擬方法

1.1 線化歐拉方程

本文的主動控制方程采用的是二維線化歐拉方程作為聲傳播方程，其表達(dá)形式如下

(1)

式中:Q0=(ρ′,u′,υ′,p′)T為瞬時(shí)變化量(密度，速度分量和壓力)。Q0=(ρ0,u0,υ0,p0)T為時(shí)均變量。式(1)的矩陣系數(shù)表示為

(2)

式中:γ=1.4為比熱比，S為聲源項(xiàng)。

經(jīng)過曲線坐標(biāo)轉(zhuǎn)換后，式(1)可寫為下式

(3)

1.2 空間離散格式

計(jì)算氣動聲學(xué)所要模擬的是物理波的傳播過程，為此所選擇的差分格式要能準(zhǔn)確反映其傳播特征(色散性，耗散性，方向性，相速度，群速度等)。高階泰勒級數(shù)截?cái)喾ㄋ@得的差分模板關(guān)注的是聲波的耗散性。它能準(zhǔn)確的逼近物理波的幅值，保證波幅計(jì)算值與實(shí)際值偏差較小。但對于聲傳播問題，差分格式不僅需要關(guān)注聲波的幅值，而且需要關(guān)注如何準(zhǔn)確描述聲波的傳播，降低振蕩，避免發(fā)散。

本文在對于式(1)中通量項(xiàng)的離散采用了Tam提出的7點(diǎn)DRP差分格式，其x方向的通量項(xiàng)可表示為式(4)，其余兩個(gè)方向的通量采用同樣的差分格式。

(4)

其中,插值系數(shù)為

(5)

圖1 有效波數(shù)隨波數(shù)變化曲線(Tam方法)

1.3 無反射邊界條件

在計(jì)算氣動聲學(xué)中，由于受到計(jì)算資源的約束，計(jì)算區(qū)域不可能無限大，因此需要人為地在計(jì)算區(qū)域中采用一種邊界條件來取代遠(yuǎn)場區(qū)域的信息。這種取代遠(yuǎn)場區(qū)域的邊界條件，其需要起到兩個(gè)作用：① 讓計(jì)算區(qū)域內(nèi)的各種物理波能無反射的通過邊界；② 讓邊界以外的任何擾動不能傳入計(jì)算區(qū)域，防止污染數(shù)值解。

論文中采用了Hu所提出的PML(Perfectly Matched Layers)邊界條件，對無分裂形式的PML方程[7-9]在曲線坐標(biāo)上進(jìn)行了轉(zhuǎn)化，并改善邊界條件的數(shù)值離散格式，增強(qiáng)其穩(wěn)定性。

遠(yuǎn)場邊界處的二維LEE(Linearized Euler Equations)方程組通過擬線性化處理，可寫為

(6)

根據(jù)Hu的方法，可知由于對流效應(yīng)，波的相速度和群速度可能沿著不同的方向發(fā)展。如果此時(shí)在PML方程中僅引入數(shù)值黏性，反射波將不能完全被消除，以致計(jì)算的穩(wěn)定性受到影響。為了解決這一難題，Hu采用了保持空間不變而對時(shí)間進(jìn)行變換的方法，即：

(7)

(8)

將其代入到式(6)中，可得以下方程

(9)

(10)

然后通過復(fù)數(shù)變易法，引入阻尼，即：

(11)

將式(11)代入式(10)可得

(12)

式(12)兩端同時(shí)乘以(ω+iσξ)(ω+iση)/ω2，可得式(13)，即：

(13)

傅里葉逆變換得

(14)

(15)

將上述偏導(dǎo)數(shù)關(guān)系式(15)代入式(14)，展開可得：

(σξ+ση)Q+σξσηq+βQ(σξξxA1+ηxA2ση)+

β(ξxA1+ηxA2)σξσηq=0

(16)

2 曲線坐標(biāo)系下人工黏性在間斷處的分部策略

對于聲學(xué)問題的數(shù)值模擬，在劃分網(wǎng)格的過程中不可避免地帶來網(wǎng)格間斷等問題。而這些間斷(物面邊界、遠(yuǎn)場邊界、網(wǎng)格交界面)易引起差分格式的色散，產(chǎn)生數(shù)值偽波。以物面邊界為例，對于空間離散采用7點(diǎn)DRP格式，假如物面不采用虛網(wǎng)格的形式，如圖2所示。垂直于物面的網(wǎng)格點(diǎn)A、B空間離散將無法采用對稱模板的7點(diǎn)DRP格式。為此在添加人工黏性項(xiàng)時(shí)，網(wǎng)格點(diǎn)A、B的人工黏性模板也是有所區(qū)別的，對于物面上方網(wǎng)格點(diǎn)A采用3點(diǎn)人工黏性模板，網(wǎng)格點(diǎn)B采用5點(diǎn)人工黏性模板，網(wǎng)格點(diǎn)C和D采用7點(diǎn)模板。網(wǎng)格交界面處可通過采用3層虛網(wǎng)格添加7點(diǎn)模板的人工黏性項(xiàng)，也可類似于物面邊界方式處理(注：程序在網(wǎng)格交界面采用了3層虛網(wǎng)格的形式)。

另外為了保證計(jì)算的穩(wěn)定性，在公共邊的奇異點(diǎn)(非常規(guī)的1-1交接網(wǎng)格頂點(diǎn)處)附近需要添加額外的人工黏性項(xiàng),如圖3所示。接下來將介紹程序中在網(wǎng)格內(nèi)部、網(wǎng)格交界面、物面邊界、公共邊奇異點(diǎn)處所采用的人工黏性項(xiàng)分布方式。

圖2 間斷處的人工黏性函數(shù)分布示意圖

圖3 奇異點(diǎn)示意圖

假設(shè)采用7點(diǎn)人工黏性模板，在一維控制方程中，添加人工黏性項(xiàng)后，可表示為

(17)

其中RΔx=a0Δx/υa，為單位網(wǎng)格雷諾數(shù)。當(dāng)網(wǎng)格存在拉伸和變形時(shí)，即計(jì)算坐標(biāo)與直角坐標(biāo)不一致，存在坐標(biāo)變換時(shí)，式(18)中人工黏性的添加需要轉(zhuǎn)換到計(jì)算坐標(biāo)系下

(18)

如果保證式(17)與式(18)中人工黏性添加的效果相同，需使得RΔξ=xξRΔx。Δξxξ為曲線坐標(biāo)中長度Δx在計(jì)算坐標(biāo)中相對應(yīng)的值。與一維方法類似，在三維空間中添加人工黏性可表示為

(19)

(20)

式中:ξ0為網(wǎng)格邊界處坐標(biāo);nΔξ在文中取為當(dāng)前網(wǎng)格塊在ξ方向的半寬。同理在其余兩個(gè)方向網(wǎng)格雷諾數(shù)的倒數(shù)分布表示為

(21)

(22)

在需要添加人工黏性的奇異點(diǎn)(ξ0,η0,ζ0)上網(wǎng)格雷諾數(shù)的倒數(shù)分布為

(23)

3 算例驗(yàn)證

為了驗(yàn)證曲線坐標(biāo)系下人工黏性格式的合理性，文中選用了NASA第二屆氣動噪聲會議[11]中初始擾動繞圓柱傳播算例。同時(shí)，該算例同時(shí)也驗(yàn)證了遠(yuǎn)場無反射邊界條件處理曲面邊界的能力。

如圖4所示，直徑r=1的圓柱的坐標(biāo)系原點(diǎn)位于圓柱的中心，點(diǎn)A、B、C分別是三個(gè)聲壓觀測點(diǎn)，坐標(biāo)分別為(r=5,θ=90°)、(r=5,θ=135°)、(r=5,θ=180°)。在距離原點(diǎn)r=4處分布了一個(gè)人工點(diǎn)聲源S，其初始賦值如下

(24)

圖4 圓柱計(jì)算模型示意圖

圓柱壁面采用的是滑移邊界條件，遠(yuǎn)場采用的是無反射邊界條件。對于控制方程的離散，采用7點(diǎn)模板的DRP格式，人工黏性的分布策略如圖5～圖7所示，從網(wǎng)格塊邊界向內(nèi)部網(wǎng)格雷諾數(shù)的倒數(shù)逐漸降低。文中計(jì)算網(wǎng)格的最小尺度為0.01，聲場計(jì)算采用統(tǒng)一的無量綱時(shí)間步長Δt=0.002。圖8(a)、圖8(b)、圖8(c)、圖8(d)分別為t=50Δt、t=500Δt、t=2 000Δt、t=5 000Δt四個(gè)時(shí)刻的聲傳播布云圖。從圖8可知，聲波隨著時(shí)間推移，向空間逐步散開。當(dāng)聲波到達(dá)圓柱時(shí)，聲波會發(fā)生部分反射現(xiàn)象，形成二次聲源向空間傳播。

圖5 ξ方向無量綱分布

圖6 η方向無量綱分布

圖7 奇異點(diǎn)處無量綱分布

(a) t=50Δt

(b) t=500Δt

(c) t=2 000Δt

(d) t=5 000Δt

圖8 不同時(shí)刻的聲傳播云圖

Fig.8 Schematic diagram of sound propagation at different times

為了定量上對比分析所采用的人工黏性項(xiàng)分布方式的合理性，圖9(a)、圖9(b)、圖9(c)分別表示觀測點(diǎn) A、B、C處的聲壓隨時(shí)間的變化曲線。圖9同時(shí)給出了計(jì)算解和解析解，其中圓圈表示解析解，實(shí)線為計(jì)算結(jié)果。由圖9可知，文中計(jì)算得到的結(jié)果與解析解吻合得很好，對比結(jié)果說明了聲波傳播程序添加人工黏性項(xiàng)的合理性及處理曲邊邊界能力。

(a) A處

(b) B處

(c) C處

圖9 觀測點(diǎn)處聲壓隨時(shí)間變化曲線

Fig.9 Curve of sound pressure at observation point

4 結(jié) 論

通過文中理論的闡述及算例的驗(yàn)證分析，可得出以下結(jié)論：

(1) 本文基于笛卡爾坐標(biāo)系的人工黏性添加策略推導(dǎo)了曲線坐標(biāo)系下的表示形式，并通過在網(wǎng)格間斷等處的正態(tài)分布策略，提出了一種適用于網(wǎng)格間斷處的人工黏性添加策略，通過標(biāo)模算例的計(jì)算分析驗(yàn)證了該方法可成功應(yīng)用于曲線坐標(biāo)系下復(fù)雜網(wǎng)格的計(jì)算聲學(xué)短波發(fā)散問題中。

(2) 針對初始擾動繞圓柱傳播算例中的計(jì)算結(jié)果，可知聲波首先向空間傳播。當(dāng)聲波到達(dá)圓柱時(shí)，形成二次聲源向空間傳播。這一聲傳播現(xiàn)象以及監(jiān)測點(diǎn)的對比分析，在驗(yàn)證網(wǎng)格添加策略合理性的同時(shí)也驗(yàn)證了遠(yuǎn)場PML邊界條件的無反射特性。