基于CFD和氣動(dòng)聲學(xué)理論的空腔自激振蕩發(fā)聲機(jī)理

楊黨國(guó),李建強(qiáng),梁錦敏

(中國(guó)空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,四川 綿陽(yáng) 621000)

0 引 言

空腔繞流廣泛存在于航空航天飛行器中如飛機(jī)起落架艙、燃燒室、飛機(jī)部件接縫、武器艙等。高速氣流流經(jīng)空腔,當(dāng)滿(mǎn)足一定的空氣動(dòng)力學(xué)條件和幾何形狀條件時(shí),由于腔口剪切流與腔內(nèi)流動(dòng)的相互作用,腔內(nèi)流動(dòng)可能出現(xiàn)強(qiáng)烈的自持振蕩,腔內(nèi)外存在復(fù)雜的非定常流動(dòng)。流場(chǎng)不僅包含渦生成、脫落與破裂,還包含流動(dòng)分離、膨脹波與激波及聲與流動(dòng)相互作用等。腔內(nèi)噪聲使空腔結(jié)構(gòu)承受較大的非定常載荷,嚴(yán)重時(shí)會(huì)危及腔內(nèi)的設(shè)備和電子器件,甚至?xí)鹂涨蛔陨斫Y(jié)構(gòu)疲勞損壞。

國(guó)外Rossiter于1964年提出了空腔流聲共振反饋模型,并給出預(yù)估振蕩頻率的半經(jīng)驗(yàn)公式[1],后來(lái)Heller提出空腔后緣處的反饋聲波速度應(yīng)為當(dāng)?shù)芈曀?對(duì)Rossiter公式進(jìn)行了修正[2]。Bilanin[3]和Tam[4]等提出了預(yù)測(cè)空腔自激振蕩的解析模型,這些模型都能在一定的精度范圍內(nèi)預(yù)測(cè)空腔振蕩的離散頻率,但都沒(méi)能解決振蕩幅值問(wèn)題。隨著CFD的發(fā)展,空腔聲學(xué)特性數(shù)值研究被廣泛開(kāi)展[5-6],國(guó)內(nèi)對(duì)空腔聲學(xué)特性的數(shù)值模擬研究起步較晚。羅柏華[7]和侯中喜[8]采用求解URANS方程對(duì)空腔聲學(xué)特性進(jìn)行了數(shù)值研究,空間采用WENO格式,時(shí)間采用Runge-Kutta法,湍流模型為S-A或B-L模型。李曉東等應(yīng)用高階低耗散低頻散格式求解二維URANS方程及標(biāo)準(zhǔn)k-ε湍流模型,并配合適當(dāng)?shù)臒o(wú)反射邊界條件,對(duì)亞聲速空腔聲學(xué)特性進(jìn)行了數(shù)值研究[9]。大渦模擬(LES)是近年來(lái)蓬勃發(fā)展的湍流數(shù)值模擬方法。與RANS采用時(shí)間平均不同,LES采用空間平均或?yàn)V波,因而可以保留相當(dāng)精度的湍流瞬態(tài)信息[10],在湍流的數(shù)值研究中有重要的發(fā)展前景,被認(rèn)為是今后可應(yīng)用于工程湍流問(wèn)題研究較有希望的數(shù)值方法[11]。

本文采用 LES與FW-H方程結(jié)合的方法著重分析了空腔自激振蕩發(fā)聲機(jī)理和腔內(nèi)聲學(xué)特性,為下一步研究空腔噪聲抑制方法提供參考與借鑒。

1 聲學(xué)模擬與計(jì)算公式

設(shè)在流體中的運(yùn)動(dòng)物體表面可用 f(x,t)=0描述,且物面內(nèi)的一個(gè)區(qū)域表示為 f(x,t)＜0,物面外的一個(gè)區(qū)域表示為 f(x,t)＞0,將FW-H方程可表示為[12]:

式中右邊三項(xiàng)分別代表厚度噪聲項(xiàng)、載荷噪聲項(xiàng)和四極子噪聲項(xiàng)。其中Δ2代表波的運(yùn)行控制,即波動(dòng)算子[13-14]:

設(shè)定模型無(wú)運(yùn)動(dòng)、無(wú)穿透,使用固定壁面作為聲源積分面,故聲源的發(fā)聲主要來(lái)自載荷噪聲p′L(xi,t),應(yīng)用的FW-H方程可簡(jiǎn)化為[15]:

2 延遲時(shí)間求解方法

對(duì)于在固定聲源積分表面上的單元,聲音輻射時(shí)間與聲壓接收時(shí)間的關(guān)系式可以表示如下[15-16]:

其中:τk為聲源積分表面上第k個(gè)聲源單元的聲音輻射時(shí)間,t為空間觀測(cè)點(diǎn)處的聲壓接收時(shí)間。在輻射時(shí)間,每一個(gè)聲源單元的聲音特性能夠計(jì)算得到,對(duì)所有聲源單元(N個(gè))進(jìn)行求和獲得聲源在輻射時(shí)間的聲音信息,并考慮輻射時(shí)間與接收時(shí)間的不同步,參考滯后時(shí)間運(yùn)算,可以獲得空間觀測(cè)點(diǎn)xi的聲壓如下:

采用迭代法求解包含延遲時(shí)間的聲波方程,獲得聲源積分表面的聲音輻射信息,從而求解空間觀測(cè)點(diǎn)聲壓。求解二維問(wèn)題時(shí),需指定一個(gè)聲源相關(guān)長(zhǎng)度(Lc)來(lái)完成聲學(xué)計(jì)算[17]。

3 流動(dòng)控制方程求解

3.1 大渦模擬

采用盒式濾波函數(shù)代入流動(dòng)控制方程,可得瞬時(shí)下大渦模擬主導(dǎo)方程[18-19],其中混合尺度的亞格子模型(SGS)定義見(jiàn)文獻(xiàn)[20]。本文采用格心形式的有限體積法求解大渦模擬主導(dǎo)方程,對(duì)流項(xiàng)采用三階迎風(fēng)偏置格式離散[21],擴(kuò)散項(xiàng)、亞格子項(xiàng)采用中心差分格式離散,空間離散為二階精度[22];采用雙時(shí)間推進(jìn)方法進(jìn)行時(shí)間項(xiàng)求解,離散格式為隱式三步二階精度格式[22]。

3.2 方法驗(yàn)證

驗(yàn)證算例是二維圓柱模型選取Revell等人試驗(yàn)研究過(guò)的圓柱模型[23],基于圓柱直徑的雷諾數(shù)90000,數(shù)值計(jì)算時(shí)間步長(zhǎng)為2×10-6s。圓柱壁面為聲源表面,聲源相關(guān)長(zhǎng)度取5D。兩個(gè)聲壓接收位置Receive-1和Receive-2,從圓柱最前端點(diǎn)順時(shí)針?lè)较?0°,離圓柱中心距離分別為35D和128D,計(jì)算區(qū)域總長(zhǎng)度為31D,寬度為15D,圓柱中心位于計(jì)算區(qū)域總寬度中點(diǎn)處(7.5D),離計(jì)算區(qū)域最前端7.5D,離計(jì)算區(qū)域計(jì)算最后端23.5D。網(wǎng)格總結(jié)點(diǎn)數(shù)為72800個(gè),圓柱壁面邊界層內(nèi)網(wǎng)格加密。邊界條件的相關(guān)參數(shù)設(shè)置均同參考文獻(xiàn)的試驗(yàn)條件一致。

表1給出了本文計(jì)算結(jié)果與參考文獻(xiàn)的平均阻力系數(shù)和描述圓柱尾部渦脫落頻率特性的斯托羅哈數(shù)St(St=f D/U,其中 f為渦脫落頻率)的對(duì)比關(guān)系。可見(jiàn),本文獲得的平均阻力系數(shù)比文獻(xiàn)[3]的計(jì)算結(jié)果更接近試驗(yàn)值;另外計(jì)算獲得的St數(shù)也更接近試驗(yàn)值;表明本文采用的大渦模擬方法正確、可靠。

表1 平均阻力系數(shù)和斯托羅哈數(shù)Table1 Averaged drag coefficient and Strouhal number

圖1給出了Receive-1和Receive-2的聲學(xué)特性,接收點(diǎn)1比接收點(diǎn)2更靠近聲源,從圖1可知接收點(diǎn)1的聲壓級(jí)較接收點(diǎn)2的大,符合聲音從近場(chǎng)向遠(yuǎn)場(chǎng)傳播耗散的規(guī)律。表2給出了聲學(xué)特性與參考文獻(xiàn)的對(duì)比,可知本文數(shù)值方法獲得的總聲壓級(jí)與聲壓峰值頻率與試驗(yàn)結(jié)果較接近,說(shuō)明此方法基本可行。

表2 兩個(gè)接收點(diǎn)處的聲學(xué)特性Table2 Aero-acoustic characteristics at receiver-1 and receiver-2

4 空腔發(fā)聲機(jī)理

4.1 計(jì)算模型

二維空腔模型長(zhǎng)深比為2(見(jiàn)圖2),來(lái)流Ma為0.64,Re為1.16×105/m;采用約30萬(wàn)的結(jié)構(gòu)網(wǎng)格,近場(chǎng)稠密遠(yuǎn)場(chǎng)逐漸稀疏;邊界條件包括黎曼不變量無(wú)反射遠(yuǎn)場(chǎng)邊界和無(wú)滑移絕熱壁條件;聲源積分運(yùn)算時(shí)聲源相關(guān)長(zhǎng)度取5倍空腔長(zhǎng)度。

圖1 兩個(gè)接收點(diǎn)的聲學(xué)特性Fig.1 Aero-acoustic characteristics at receiver-1 and receiver-2

圖2 計(jì)算區(qū)域及邊界條件Fig.2 Computational domain and boundary conditions

4.2 自激振蕩

關(guān)于空腔自激振蕩發(fā)聲的機(jī)理,研究者們有不同的解釋,然而,目前并沒(méi)有徹底全面的公認(rèn)解釋,需要進(jìn)一步分析。Rossiter和Heller等人提出的針對(duì)開(kāi)式空腔(長(zhǎng)深比小于10)剪切層激發(fā)空腔自激振蕩誘導(dǎo)空腔噪聲得到了廣泛關(guān)注。圖3給出了一個(gè)周期內(nèi)不同時(shí)刻的渦量等值線(xiàn)圖,在四分之一周期,腔前緣的剪切層內(nèi)和腔后壁的下游各形成了一個(gè)渦(圖3a),在二分之一周期剪切層內(nèi)的渦脫落,經(jīng)過(guò)空腔中部后繼續(xù)向下游發(fā)展(圖3b),直到腔后壁處流向空腔下游(圖3c),在一個(gè)周期結(jié)束時(shí)在腔前緣的剪切層內(nèi)又形成了新的渦(圖3d),可見(jiàn),來(lái)流剪切層在腔前緣處分離,在空腔上方形成了具有一定脫落頻率的渦,并與空腔流場(chǎng)結(jié)構(gòu)相互作用,產(chǎn)生了復(fù)雜的非定常特性。

圖3 空腔瞬時(shí)渦量等值線(xiàn)Fig.3 The instantaneous vorticity contours in the cavity

圖4給出了計(jì)算10個(gè)周期后一個(gè)周期內(nèi)空腔流場(chǎng)的渦量等值線(xiàn)圖,可以看出空腔流場(chǎng)結(jié)構(gòu)在不同時(shí)刻的變化情況。四分之一周期時(shí),在腔前緣剪切層內(nèi)形成了渦,并脫落經(jīng)過(guò)空腔中部向后繼續(xù)向下游發(fā)展到腔后部,此時(shí)在腔前緣處又形成了新的脫落渦(如圖4a所示);在二分之一周期時(shí),剪切層內(nèi)的第一個(gè)脫落渦與腔后壁相撞,并產(chǎn)生了噪聲,一次聲波經(jīng)過(guò)腔內(nèi)反饋回路向上游傳播并到達(dá)腔前壁后,與前壁相撞,并激發(fā)新的渦生成(圖4b);在剪切層內(nèi)新的渦又脫落,經(jīng)腔中部后又與腔后壁相撞并導(dǎo)致二次聲波的產(chǎn)生,二次聲波經(jīng)腔內(nèi)的反饋回路又傳播到腔前壁,再次誘導(dǎo)新的渦生成(圖4c);故而在腔內(nèi)因剪切層內(nèi)的脫落渦與后壁相撞產(chǎn)生的一次聲波經(jīng)腔內(nèi)的反饋回路到達(dá)腔前壁,擾動(dòng)剪切層并誘導(dǎo)新的脫落渦生成,激發(fā)二次聲波,在腔內(nèi)形成了脫落渦—一次聲波—新的脫落渦—二次聲波的一個(gè)反饋機(jī)制,流動(dòng)形成自激振蕩(圖4d)。

圖4 瞬時(shí)渦量等值線(xiàn)發(fā)展歷程Fig.4 The developing process of vorticity contours

4.3 聲學(xué)特性

圖5分別給出了 Rowley等人[5]關(guān)于空腔繞流空間聲場(chǎng)聲壓級(jí)等值線(xiàn)的計(jì)算結(jié)果和本文關(guān)于空間六個(gè)不同觀測(cè)點(diǎn)的聲壓級(jí),從結(jié)果可看出,腔內(nèi)聲壓級(jí)較高,隨著離腔壁面距離的增大,聲音信號(hào)衰減,聲壓級(jí)降低;聲音主要向左上方輻射,這主要因?yàn)榧羟袑觾?nèi)的脫落渦與后壁相撞產(chǎn)生的聲波向上游傳播,到達(dá)腔前壁后又穿透剪切層向腔外空間場(chǎng)輻射?？梢?jiàn),剪切層內(nèi)的流動(dòng)特性是影響空腔聲學(xué)特性關(guān)鍵所在。

圖5 不同測(cè)點(diǎn)的聲壓級(jí)Fig.5 SPL at different measurement positions

Rossiter關(guān)于開(kāi)式空腔聲學(xué)特性的振蕩頻率給出了半經(jīng)驗(yàn)預(yù)測(cè)公式,后來(lái) Heller等提出在腔后壁處向上游傳播的聲波速度應(yīng)為當(dāng)?shù)芈曀?不應(yīng)是遠(yuǎn)場(chǎng)聲速,獲得了修正的Rossiter公式,分別如式(16)和式(17),其中k為渦遷移速度與自由流速度之比,α為渦通過(guò)與產(chǎn)生聲壓之間的時(shí)間延遲因子,n為模態(tài)階數(shù)。

圖6給出了接收點(diǎn)1和2的聲壓頻譜特性與功率譜密度特性,空腔流動(dòng)自激振蕩一階模態(tài)和二階模態(tài)的St分別為0.29和0.73,同Rossiter和Heller的預(yù)測(cè)基本吻合,如表所示?？涨粌?nèi)的自激振蕩主要以一階模態(tài)為主,聲壓級(jí)幅值均出現(xiàn)在一階振蕩模態(tài);一階振蕩模態(tài)是因聲波從后壁接收點(diǎn)2傳播至接收點(diǎn)1而激發(fā)剪切層中渦脫落與后壁相撞形成,故接收點(diǎn)2的一階模態(tài)聲壓峰值大于接收點(diǎn)1,而二階模態(tài)是因接收點(diǎn)1處剪切層中的脫落渦與腔后壁相撞,形成反饋回路后再次傳播至腔前壁形成,故接收點(diǎn)1的二階模態(tài)聲壓峰值大于接收點(diǎn)2(如圖6a和6b)。腔內(nèi)噪聲的能量主要集中在峰值頻率303.07Hz(如圖6c和6d所示)。

圖6 兩個(gè)接收點(diǎn)的聲學(xué)特性Fig.6 Aeroacoustic characteristics at two receivers

表3 空腔自激振蕩模態(tài)分析Table3 Self-sustained oscillation modes in the cavity

5 結(jié) 論

通過(guò)本文研究,可得出以下結(jié)論:

(1)圓柱繞流噪聲計(jì)算結(jié)果與國(guó)外試驗(yàn)結(jié)果基本吻合表明本文采用的CFD技術(shù)和氣動(dòng)聲學(xué)時(shí)域理論方法基本正確、可靠。

(2)空腔自激振蕩形成和腔內(nèi)噪聲產(chǎn)生的主要原因是氣流在空腔前緣發(fā)生分離,在空腔上方形成剪切層,剪切層中存在渦的生成、脫落,脫落渦與空腔后壁相撞后產(chǎn)生一次聲波,一次聲波通過(guò)空腔向前傳播至空腔前緣氣流分離處,激發(fā)剪切層中新的渦脫落,新的脫落渦與空腔后壁再次相撞產(chǎn)生二次聲波,二次聲波又通過(guò)空腔向前傳播至空腔前緣再次激發(fā)新的渦脫落;一定的相位條件下,在空腔內(nèi)形成了具有一定頻率的脫落渦-聲波-新的脫落渦-新的聲波的流致噪聲反饋回路。

(3)研究的空腔內(nèi)聲壓幅值主要出現(xiàn)在一階和二階振蕩模態(tài),聲音能量主要集中在低頻300Hz左右,與Rossiter和Heller等人的試驗(yàn)結(jié)果類(lèi)似,為下一步開(kāi)展空腔降噪方法研究提供了一定的參考。

[1]ROSSITER J E.Wind tunnel experiments of the flow over rectangular cavities at subsonic and transonic speeds[R].A RCR and M.,1964:3458.

[2]HELLER H H,BLISS D B.Aerodynamically induced pressure oscillations in cavities[R].Physical Mechanisms and Suppression Concepts.AFFDL-T R-74-133,1975.[3]BILANIN A J,COVERT E E.Estimation of possible excitation frequencies for shallow rectangular cavities[J].AIAA Journal,1973,11(3):347-351.

[4]TAM CKW,BLOCK PJW.On the tones and pressure oscillations induced by flow over rectangular cavities[J].J.Fluid Mech.,1978,89(2):373-399.

[5]ROWLEY C W,COLONIUS T,BASU A J.On selfsustained oscillations in two-dimensional compressible flow over rectangular cavities[J].J.Fluid Mech.,2002,455:315-346.

[6]LILLBERG E,FUREBY C.Large eddy simulations of supersonic cavity flow[R].AIAA 2000-2411.

[7]羅柏華.二維高亞聲速空腔流激振蕩的數(shù)值模擬研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),2002,20(1):84-88.(LUO B H.Numerical simulation of 2D high subsonic cavity flow oscillation[J].Acta Aerodynamic Sinica,2002,20(1):84-88.(in Chinese))

[8]侯中喜,夏剛,秦子增.三維超聲速開(kāi)式空腔振蕩特性研究[J].國(guó)防科學(xué)技術(shù)大學(xué)學(xué)報(bào),2004,26(6):1-4.(HOU Z X,XIA G,QIN Z Z.The numerical analysis of oscillatory characteristics in 3D supersonic open cavity[J].Journal ofNational University ofDefense Technology,2004,26(6):1-4.(in Chinese))

[9]李曉東,劉靖東,高軍輝.空腔流激振蕩發(fā)聲的數(shù)值模擬研究[J].力學(xué)學(xué)報(bào),2006,38(5):599-604.(LI X D,LIU J D,GAO J H.Numerical simulation of flow-induced oscillation and sound generation in a cavity[J].Chinese Journal of Theoretical and Applied Mechanics,2006,38(5):599-604.(in Chinese))

[10]GALPERIN A,ORSZAG S A.Large eddy simulation of complex engineering and geophysical flows[R].Cambridge University Press,1993.

[11]KYOUNG S C,SEUNG O P.Hybrid RANS/LES simulation of deep cavity flow[R].42nd AIAA Aerospace Sciences M eetingand Exhibit[C].Reno,Nevada,2004.

[12]FFOWCS WILLIAMS J E,HAWKINGS D L.Sound generation by turbulence and surfaces in arbitrary motion[R].London:Proc.Roy.Soc.,A264,1969:321-342.

[13]BRENTNER K S.Numerical algorithms for acoustic integrals-the devil is in the details[R].AIAA1996-1706.[14]FA RASSAT F,BREN TNER K S.The acoustic analogy and the prediction of rotating blades[J].Theoretical and Computational Fluid Dynamics,Springer-Verlag,1998,10:155-170.

[15]CASALINO D.An advanced time approach for acoustic analogy predictions[J].J.S.V.,2003,261:583-612.

[16]KIM S E,YI D,EVANGELOS K K.A versatile implementation of acoustic analogy based noise prediction method in a general-purpose CFD[R].AIAA 2003-3202.

[17]SARIGUL K N,DIETZ D,KARNOPP D,DUMMER J.A computational aeroacoustic method for near and far field vehicle noise predictions[R].AIAA 2001-0513.

[18]FELTEN F,FAU TRELLE Y,TERRAIL Y D,et al.Numerical modeling of electrognetically-riven turbulent flows using LES methods[J].Applied Mathematical Modelling,2004,28(1):15-27.

[19]DORIS L,CHRISTIAN T,PHUOC L T.LES of spatially developing 3D compressible mixing layer[J].Computational Fluid Mechanics,2000,328(7):567-573.

[20]WU H W,PERNG S W.LES analysis of turbulent flow and heat transfer in motored engines with various SGS models[J].International Journal ofHeat and Mass transfer,2002,45(11):2315-2328.

[21]ARADAG S,KNIGHT D D.Simulation of supersonic cavity flow using 3D RANS equations[R].AIAA 2004-4966.

[22]KIM S E,MATHUR S R,MURTHY J Y,CHOUDHURY D.A Reynolds-Averaged Navier Stokes solver using unstructured mesh based finite-Volume scheme[R].AIAA 1998-0231.

[23]REVELL J D,PRYDZ R A,HAYS P.Experimental study of airframe noise vs.drag relationship for circular cylinders[R].Lockheed Report 28074,1997.

[24]BRENTNER K S,COX J S,RUMSEY C L,et al.Computation of sound generated by flow over a circular cylinder:An acoustic analogy approach[R].Presented at the Second Computational Aeroacoustics Workshop on Benchmark Problems[C].Tallahassee,FL.,1996.

[25]GHARIB M,ROSHKO A.The effect of flow oscillations on cavity drag[J].J.Fluid Mech.,1987,177:501-530.