有均勻平均流的圓環(huán)管道聲模態(tài)及其應(yīng)用

柳占新 黃其柏

1.華中科技大學(xué)數(shù)字制造裝備與技術(shù)國家重點實驗室,武漢,430074 2.中國船舶重工集團(tuán)710研究所,宜昌,443003

0 引言

氣動聲學(xué)是一門研究氣體運(yùn)動產(chǎn)生聲音以及聲音在運(yùn)動氣體中的傳播過程的一門科學(xué)。由于氣動聲學(xué)方程復(fù)雜,要解析地描述復(fù)雜幾何形體的氣動聲學(xué)場(包括聲的輻射與傳播)幾乎是一件不可能的事。計算氣動聲學(xué)(computational aeroacoustics,CAA)的出現(xiàn)為打開氣動聲學(xué)的神秘大門提供了神奇的鑰匙。一些簡單幾何形體的氣動聲學(xué)問題,相對較容易求出解析解或者進(jìn)行實驗驗證,這類問題對驗證CAA方法的準(zhǔn)確性提供了素材,如CAA的一系列基準(zhǔn)問題[1]。有平均流的圓環(huán)管道中的聲傳播問題就屬于這類相對簡單的問題,引起了大量學(xué)者的興趣[2-5]。有平均流的圓環(huán)管道中的聲傳播的研究意義不僅在于給CAA提供論據(jù),更重要的是它為研究更復(fù)雜的問題提供了思路。如人們用圓(環(huán))管道近似代替渦扇發(fā)動機(jī)進(jìn)氣道,然后研究吸聲材料添加方法[4,6-8]。本文將由線性歐拉方程出發(fā),推導(dǎo)有均勻平均流的情況下圓環(huán)管道內(nèi)的聲傳播、聲功率的解析表達(dá)式,利用單個聲模態(tài)來模擬渦扇發(fā)動機(jī)進(jìn)氣道的聲傳播現(xiàn)象。

1 有平均流的管道聲模態(tài)

圓環(huán)管道的模型如圖1所示,假定其軸向無限長,內(nèi)外徑分別為D i和D o,管壁都是剛性硬壁,管道內(nèi)的平均流場均勻,其馬赫數(shù)為Ma。對于圓環(huán)管道,自然坐標(biāo)應(yīng)該選用柱坐標(biāo),本文以(x,r,φ)代表柱坐標(biāo)在軸向、徑向和旋轉(zhuǎn)方向上的分量。管道內(nèi)聲波的控制方程為線性歐拉方程,表述為

式中,ρ、u、v、w、p分別為流體(氣體)的密度、x方向速度、r方向速度、φ方向速度、壓力;γ為氣體絕熱指數(shù);撇號“′”表示聲學(xué)擾動量,上橫線表示平均流。

圖1 固壁圓環(huán)管道模型示意圖

為了便于計算,所有變量都采用以下基準(zhǔn)進(jìn)行量綱一處理：長度基準(zhǔn)為管壁外徑D(即圖1中的D o);速度基準(zhǔn)為管道中的音速a0;密度基準(zhǔn)為管道中空氣密度 ρ0。時間為D/a0,故壓力為 ρ0a20。

量綱一處理后,管道的外徑為1。用σ代表圓環(huán)管道內(nèi)外徑之比,即σ=D i/D o,從而管道的內(nèi)徑Di=σ。所以管道外半徑Ro=1/2,內(nèi)半徑Ri=σ/2。

由等熵關(guān)系知

式(2)可以改寫為

將p看成是ρ的函數(shù),然后式(3)兩邊同時對ρ求導(dǎo)可得

由于d p和dρ表示密度變化的小量——聲學(xué)擾動量,所以 p′=d p,ρ′=dρ。同時 ,

將式(5)代入式(4),可得

其中V=(u,v,w)為速度矢量。消去式(8)中的速度矢量,可得

或者展開寫成

式(9)或式(10)就是有平均流的情況下圓環(huán)管道的聲傳播控制方程。為了得到式(10)的解,還需添加適當(dāng)?shù)倪吔鐥l件。對于剛性壁,應(yīng)該添加無穿透邊界條件,即

對于圓環(huán)管道里傳播的聲音,可以認(rèn)為其由很多個Fourier模態(tài)組成,而單個Fourier模態(tài)的形式為

式中,p(r)為實數(shù);m為周向模數(shù)(為整數(shù));k為軸向波數(shù);ω為角速度。

而對于其他流體變量,也可以作相似的假設(shè)。將式(12)代入式(10)可得

式(13)是一個典型的貝塞爾方程,其通解為

式中,A、B為常數(shù);Jm、Ym分別為m階的第一類和第二類貝塞爾函數(shù);下標(biāo)n為徑向模數(shù)。

從式(13)可以看出,如果βmr為其解,那么-βmn也是其解,因此可以只考慮非負(fù)的βmn。從而在式(14)中,m和βmn有以下幾種組合方式：

(1)m=βmn=0。由于第二類貝塞爾函數(shù)在0處趨于-∞,為了保證解有界必有B=0,此時解退化為一個常數(shù),也就是p=A,這個解就是平面波。

(2)m ＞0且βmn=0。此時p為零,我們只關(guān)注非零解,所以可以排除此類情況。

(3)m ≥0且βmn＞0。此時根據(jù)邊界條件式(11)可得

式(15)中,為了使A和B有解,必須讓左邊矩陣的行列式等于零,即

式(16)就是圓環(huán)管道內(nèi)聲波的色散關(guān)系式,求βmn的根的方法將在下一節(jié)討論。找出了 βmn的值,也就確立了波數(shù)k和角速度ω的關(guān)系。將β2mn=(ω-kMa)2-k2改寫為

從而

式中,k+、k-分別為聲波向下(游)和上(游)傳播。

聲波在傳播過程中應(yīng)該既不被衰減也不被放大,因此 k必須為實數(shù),并且必須有ω2≥(1-Ma2)β2mn。ω2=(1-Ma2)β2mn對應(yīng)的頻率稱為截止頻率,ω2＞(1-Ma2)β2mn對應(yīng)的聲模態(tài)稱為傳播聲模態(tài),ω2＜(1-Ma2)β2mn對應(yīng)的聲模態(tài)稱為截止聲模態(tài)。

在式(18)中,聲波的角速度ω一般為已知數(shù),也可以通過聲波的頻率求得,即 ω=2πf。但需要注意的是,由于這里所有物理量都是量綱一單位,頻率也應(yīng)該為量綱一單位。頻率為 f 0(Hz)的聲波對應(yīng)的量綱一頻率為

通過式(16)求得βmn,ω也已知,則可通過式(18)求得波數(shù)k,從而式(14)中聲壓的解變?yōu)?/p>

式中,Amn為聲壓幅值,由聲功率決定。

將式(20)代入式(12)可得

將式(21)代入式(7)中,可以得到聲波各個變量的表達(dá)式：

2 求解βmn的方法

為了求解圓環(huán)管道中的聲模態(tài),必須首先通過式(16)才能得到βmm。然而,此方程為超越方程,不可能得到簡單的解析解,所以需要利用數(shù)值解法,如牛頓迭代法求解。在利用牛頓迭代法求解的過程中,需要給定初值。為了得到這些初值,可用薄管近似方法。假定非常小,然后用R c代替式(13)中的r,可得

令 p=eλr,則式(23)變?yōu)?/p>

從而

式(23)的解為

將邊界條件式(11)代入式(26),可得

為了得到非零解,式(27)中系數(shù)矩陣的行列式必須為零,即

將式(25)代入式(28)可得

當(dāng)管道內(nèi)外徑之比σ比較大的時候,這種薄管近似給出的初始值比較精確。然而當(dāng)σ比較小的時候,圓環(huán)管道和薄管有很大的差異,如果用式(29)估計初值,則會產(chǎn)生根的排列混亂和丟根的情況。此時可以采用多次薄管近似逐漸逼近的方法或者二分法來估計初值。設(shè)函數(shù)

則式 (16)的根就是由式(30)確定的函數(shù)f(β)的零點。m=22且σ=0.4時,f隨β變化的曲線如圖2所示。根據(jù)函數(shù)圖像的形狀,完全可以用二分法求根的初值或者直接求根。

圖2 m=22且σ=0.4時 f(β)的函數(shù)曲線圖

3 管道中的聲功率

式(22)中與聲學(xué)變量幅值相關(guān)的常數(shù)Amn由聲功率決定。如果聲功率為P,則聲功率級定義為

其中P0為基準(zhǔn)聲功率,P0=10-12W。令P為量綱一單位的聲功率,則

另一方面,平均流的聲強(qiáng)為[9]

式中,p、u分別為聲擾動的壓力和速度的實部。

對單個的聲模態(tài)有

所以管道中的聲功率為

不難證明,當(dāng)T→∞時,對于任意復(fù)數(shù)C和G有

式中,G*為G的共軛。

因此

將式(37)和式(35)代入式(33)可得單個聲模態(tài)的聲功率：

將式(22)代入式(38),可得

至此,如果給定單個聲模態(tài)的聲功率、頻率、周向模態(tài)m以及平均流的馬赫數(shù)Ma,則可以完全確定圓環(huán)管道中聲波各個變量的解析表達(dá)式。

4 算例

下面利用有平均流的圓管聲模態(tài)來模擬渦扇發(fā)動機(jī)進(jìn)氣道的噪聲傳播過程。某發(fā)動機(jī)進(jìn)氣道的結(jié)構(gòu)如圖3所示,其中 x代表進(jìn)氣道的軸向,r代表進(jìn)氣道的半徑方向(計算將在柱坐標(biāo)中進(jìn)行)。發(fā)動機(jī)工作時風(fēng)扇將氣體由進(jìn)氣道吸入到內(nèi)外涵道。風(fēng)扇產(chǎn)生的噪聲是發(fā)動機(jī)的一個主要噪聲源,這些噪聲可以沿著氣流或者逆著氣流傳播。這里僅研究風(fēng)扇噪聲逆著氣流在進(jìn)氣道傳播的過程。為了便于計算,作如下簡化：①整個進(jìn)氣道關(guān)于x軸對稱;②將圓環(huán)管道中單個聲模態(tài)作為輸入的風(fēng)扇噪聲。從圖3可以看出,這個進(jìn)氣管道在

圖3 渦扇發(fā)動機(jī)進(jìn)氣道結(jié)構(gòu)圖

風(fēng)扇面附近近似為一個圓環(huán)形的管道,因此可以認(rèn)為,風(fēng)扇產(chǎn)生的噪聲是由無數(shù)個圓環(huán)管道聲模態(tài)疊加而成的。假定其中的單個聲模態(tài)在周向解析(m已知,為輸入?yún)?shù)),從而所有的聲學(xué)量都可以寫成如下形式：

其中帶有上尖括號的量均為復(fù)數(shù)。

在以上簡化下,控制方程變?yōu)?/p>

從而將一個三維問題轉(zhuǎn)化為二維問題。計算中先利用保形映射方法將物理空間(x-r)中轉(zhuǎn)子和機(jī)匣的壁面曲線映射到計算空間(ζ-η)中的兩條相互平行的直線上,如圖4所示。由于進(jìn)氣道滿足軸對稱條件,所有計算只需要在第二象限進(jìn)行,所以圖4中僅顯示轉(zhuǎn)子和機(jī)匣在第二象限中的像。映射后轉(zhuǎn)子和機(jī)匣的像之間的距離為h。計算時首先在計算平面中設(shè)計網(wǎng)格,然后再將網(wǎng)格逆變換到物理空間中。考慮到計算需求和效率,計算區(qū)域在計算平面中的大小選為6h×4h,這個區(qū)域?qū)?yīng)在物理坐標(biāo)中的范圍如圖5所示。

圖4 轉(zhuǎn)子和機(jī)匣壁面曲線在計算空間中的像

圖5 物理空間中的計算區(qū)域圖

計算中時間和空間離散采用多網(wǎng)格尺度 -多時間步色散關(guān)系保持(DRP)方法[10]進(jìn)行。計算中采用如下邊界條件：①A-B-C-D邊界采用輻射邊界條件[11],使聲波能無反射地穿過計算區(qū)間;②A-K采用軸對稱邊界;③機(jī)匣和轉(zhuǎn)子壁面上采用滑移邊界條件,利用外伸點(ghost point)法[12]實施;④風(fēng)扇表面(E-H)采用理想匹配層(PML)邊界[13]輸入聲波并吸收反射聲波,保持?jǐn)?shù)值穩(wěn)定。

作為模擬算例,這里僅計算輸入頻率為1614Hz,聲功率級為120dB,周向模數(shù)m=22時n=1和n=4兩種情況。在計算聲場之前,先計算平均流場,使得風(fēng)扇表面的馬赫數(shù)為0.3。計算后的瞬時聲壓云圖如圖6和圖7所示,可以看出,風(fēng)扇面輸入的不同聲模態(tài)在進(jìn)氣道以及進(jìn)氣道附近的區(qū)域傳播過程差別較大,這也正是眾多學(xué)者研究進(jìn)氣道噪聲傳播機(jī)理的原因。

5 結(jié)論

(1)聲波在平均流中的傳播的控制方程為線性歐拉方程,在均勻流的圓環(huán)管道中,平均流只沿著軸向方向流動,從而可以使方程大大簡化。

(2)傳播中的聲波可以看成是多個Fourier模態(tài)的組合。對于圓環(huán)管道中的單個Fourier模態(tài),其軸向波數(shù)k、周向模數(shù)m、角速度ω以及平均流的馬赫數(shù)Ma等必須滿足色散關(guān)系式。

(3)在其他條件給定時,圓環(huán)管道中的聲波存在一個截止頻率,大于截止頻率的聲模態(tài)可以在管中傳播,而小于截止頻率的聲模態(tài)則不能在管中傳播。

(4)給定了聲功率級,平均流的馬赫數(shù),入射聲波的頻率、旋轉(zhuǎn)模數(shù)、徑向模數(shù)和管的幾何尺寸以后,便唯一確定了單個Fourier模態(tài)的聲波的解析表達(dá)式。

(5)可利用圓環(huán)管道的聲模態(tài)作為風(fēng)扇產(chǎn)生的噪聲來研究渦扇發(fā)動機(jī)進(jìn)氣到噪聲傳播規(guī)律,非均勻流場對不同聲波模態(tài)的散射效果不同。

圖6 輸入聲波m=22,n=1,PWL=120d B時的瞬時聲壓圖

圖7 輸入聲波m=22,n=4,PWL=120d B時的瞬時聲壓圖

[1] Dahl M D,Envia D,Huff D,et al.Fourth Computational Aeroacoustics(CAA)Workshop on Benchmark Problems[C]//Brook Park,Ohio：NASA Conference Publication,2004：212954.

[2] Rienstra S W.A Classification of Duct M odes Based on Surface Waves[J].Wave Motion,2003,37(2)：119-135.

[3] Ayub M,Tiwana M H,Mann A B.Propagation of Sound in a Duct with Mean Flow[J].Communications in Nonlinear Science and Numerical Simulation,2009,14(9/10)：3578-3590.

[4] Tam C K W,Ju H,Chien EW.Scattering of Acoustic Duct Modes by Axial Liner Splices[J].Journal of Sound and Vibration,2008,310(4/5)：1014-1035.

[5] Gabard G,Astley R J.A Computational Modematching Approach for Sound Propagation in Three-dimensional Ducts with f low[J].Journal of Sound and Vibration,2008,315(4/5)：1103-1124.

[6] Yang B,Wang T Q.Investigation of the Influenceof Liner Hard-splices on Duct Radiation/Propagation and Mode Scattering[J].Journal of Sound and Vibration,2008,315(4/5)：1016-1034.

[7] McAlpine A,Astley R J,Hii V J T,et al.Acoustic Scattering by an Axially-segmented Turbofan Inlet Duct Liner at Supersonic Fan Speeds[J].Journal of Sound and Vibration,2006,294(4)：780-806.

[8] Brun C,Goossens T.3D Coherent Vortices in the Turbulent Near Wake of a Square Cylinder[J].Comptes Rendus Mcanique,2008,336(4)：363-369.[9] Atassi O V.Computing the Sound Power in Nonuniform Flow[J].Journal of Sound and Vibration,2003,266(1)：75-92.

[10] Tam C K W,Kurbatskii K A.Multi-size-mesh Multi-time-step Dispersion-relation-preserving Scheme for Multiple-scales Aeroacoustics Problems[J].International Journal of Computational Fluid Dynamics,2003,17：119-132.

[11] Tam C K W,Webb JC.Dispersion-relation-preserving Finite Difference Schemes for Computational Acoustics[J].Journal of Computational Physics,1993,107(2)：262-281.

[12] Tam C K W,Dong Z.Wall Boundary Conditions for High-order Finite-difference Schemes in Computational Aeroacoustics[J].Theoretical and Computational Fluid Dynamics,1994,6(5/6)：303-322.

[13] Hu F Q.Development of PML Absorbing Boundary Conditions for Computational Aeroacoustics：A Progress Review[J].Computers and Fluids,2008,37(4)：336-348.