楊黨國,李建強,梁錦敏

(中國空氣動力研究與發(fā)展中心空氣動力學國家重點實驗室,四川 綿陽 621000)

0 引 言

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

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

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

1 聲學模擬與計算公式

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

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

設定模型無運動、無穿透,使用固定壁面作為聲源積分面,故聲源的發(fā)聲主要來自載荷噪聲p′L(xi,t),應用的FW-H方程可簡化為[15]:

2 延遲時間求解方法

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

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

采用迭代法求解包含延遲時間的聲波方程,獲得聲源積分表面的聲音輻射信息,從而求解空間觀測點聲壓。求解二維問題時,需指定一個聲源相關(guān)長度(Lc)來完成聲學計算[17]。

3 流動控制方程求解

3.1 大渦模擬

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

3.2 方法驗證

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

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

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

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

表2 兩個接收點處的聲學特性Table2 Aero-acoustic characteristics at receiver-1 and receiver-2

4 空腔發(fā)聲機理

4.1 計算模型

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

圖1 兩個接收點的聲學特性Fig.1 Aero-acoustic characteristics at receiver-1 and receiver-2

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

4.2 自激振蕩

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

圖3 空腔瞬時渦量等值線Fig.3 The instantaneous vorticity contours in the cavity

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

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

4.3 聲學特性

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

圖5 不同測點的聲壓級Fig.5 SPL at different measurement positions

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

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

圖6 兩個接收點的聲學特性Fig.6 Aeroacoustic characteristics at two receivers

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

5 結(jié) 論

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

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

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

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

[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].空氣動力學學報,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]侯中喜,夏剛,秦子增.三維超聲速開式空腔振蕩特性研究[J].國防科學技術(shù)大學學報,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].力學學報,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.