柳占新 黃其柏

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

0 引言

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

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

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

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

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

為了便于計(jì)算,所有變量都采用以下基準(zhǔn)進(jìn)行量綱一處理：長(zhǎng)度基準(zhǔn)為管壁外徑D(即圖1中的D o);速度基準(zhǔn)為管道中的音速a0;密度基準(zhǔn)為管道中空氣密度 ρ0。時(shí)間為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)可以改寫(xiě)為

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

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

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

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

或者展開(kāi)寫(xiě)成

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

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

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

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

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

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

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

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

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

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

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

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

從而

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

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

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

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

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

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

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

2 求解βmn的方法

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

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

從而

式(23)的解為

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

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

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

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

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

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

3 管道中的聲功率

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

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

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

式中,p、u分別為聲擾動(dòng)的壓力和速度的實(shí)部。

對(duì)單個(gè)的聲模態(tài)有

所以管道中的聲功率為

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

式中,G*為G的共軛。

因此

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

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

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

4 算例

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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

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

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

[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.