三維孔腔流激噪聲的大渦模擬與聲學類比預報與驗證研究

張 楠，沈泓萃，朱錫清，姚惠之，謝 華

（中國船舶科學研究中心，江蘇 無錫 214082）

1 前 言

流聲耦合是流體力學與聲學的交叉學科，是國際矚目的重點與難點。流激噪聲或稱流動發(fā)聲的預報是流聲耦合領(lǐng)域的重要組成部分，具有相當?shù)募夹g(shù)難度。這是因為水是可壓縮性極小的流體介質(zhì)，要研究的問題一般是低馬赫數(shù)的流動，而水中運載體的雷諾數(shù)都很大，聲轉(zhuǎn)換效率相當?shù)停虼怂性肼曨A報是十分復雜和困難的。研究水中流激噪聲的預報方法，探討渦旋與聲輻射之間的相互關(guān)系，具有重要意義。

近年來，大渦模擬（LES）方法的逐漸成熟以及計算機的快速發(fā)展，使得計算一些試驗狀態(tài)下的流動發(fā)聲成為可能。動態(tài)Smagorinsky模型（Germano 1991，Lilly 1992）[1-2]對于大渦模擬的可靠性與準確性做出了重大的改善。 Seror等人 （2000，2001）[3-4]，He等人 （2002，2004）[5-6]的工作都研究了動態(tài)Smagorinsky模型在預報流動輻射噪聲方面的可靠性。Wang，F(xiàn)reund與Lele（2006）[7]指出對于實際應用的雷諾數(shù)范圍，在可預見的將來，LES是一種最有希望取得重大突破的用于流激噪聲預報的CFD模擬方法。Wang與Moin（2000）[8]利用LES結(jié)合Lighthill理論模擬了翼型支桿的隨邊噪聲，弦長雷諾數(shù)為Re=2.1×106，其計算結(jié)果與 Blake（1975）的試驗結(jié)果進行了對比，達到了合理的吻合程度。Wang（2004）[9]通過計算弦長雷諾數(shù)為Re=1.5×105的風扇葉片模型非定常壓力的時空特性，進一步驗證了LES方法。這些工作都為預報復雜流動的輻射噪聲奠定了基礎(chǔ)。

孔腔流動從屬于與自持振蕩密切相關(guān)的一類基本流動范疇?？浊涣鲃幽軌蛞鸾Y(jié)構(gòu)振動與疲勞，激發(fā)很強的噪聲，其中包含了剪切層不穩(wěn)定性，湍流與結(jié)構(gòu)及噪聲之間的相互作用等流固耦合、流聲耦合復雜現(xiàn)象，給理論分析、數(shù)值模擬和試驗測量帶來了很大的挑戰(zhàn)。在過去的五十多年里，這一類型的流動被廣泛的研究過，最初的工作是Roshko在1955年所做的試驗測量與分析[10]。Rossiter（1964）[11]給出了辨識流動模態(tài)振蕩共振頻率的公式。近年來，人們又做了大量的試驗和計算研究。

Ask等人（2003）[12]采用不可壓縮CFD代碼CALC-BFC計算了長深比L/D=4的二維長方形孔腔的繞流，采用Curle方程計算了遠場9個接收點處的輻射噪聲。來流馬赫數(shù)M∞=0.15。計算充分說明了壁面脈動壓力是重要的偶極子聲源。計算得到的聲壓譜趨勢與參考結(jié)果完全一致，幅值最大誤差為6dB。

Ching（2004）[13]采用MILES方法，利用非結(jié)構(gòu)化網(wǎng)格研究了低馬赫數(shù)下的二維“L”型孔腔發(fā)聲問題（M∞=0.0788～0.147）。計算域中含85300非結(jié)構(gòu)化網(wǎng)格。在M∞=0.147時，對于流體動力振蕩引起的聲諧振峰，頻率的計算誤差為1%，幅值的計算誤差為1dB，對于Hemholtz共振引起的聲諧振，頻率的計算誤差為2.2%，幅值的計算誤差為17dB。對于其他諧振峰，計算得到的聲壓譜幅值都明顯大于試驗結(jié)果。

在空氣動力學領(lǐng)域，M 219孔腔是一種被廣泛研究的經(jīng)典模型。M219孔腔長（L）、寬（W）、深（D）之比為，L:W:D=5:1:1，來流馬赫數(shù)M∞=0.85，是接近于跨音速狀態(tài)下的一種高雷諾數(shù)流動問題。研究者做了大量的脈動壓力與噪聲試驗，并采用這些詳細的試驗結(jié)果來校驗各種數(shù)值計算方法的可信度。

瑞典國防研究署（SDRA）的Peng（2005）[14]對于M 219孔腔的非定常流動進行了數(shù)值模擬，計算工作被SDRA當作一種基準式的計算（baseline-type computation），為隨后所有該類型的計算提供可資參考的技術(shù)細節(jié)。研究者首先在跨音速風洞中測量了脈動壓力，給出了脈動壓力的頻譜，試驗清晰反映了不同頻率下的諧振峰。作者使用結(jié)構(gòu)化的NS求解器EURANUS，采用非定常RANS方法結(jié)合SA模型，數(shù)值模擬了二維孔腔的脈動壓力。計算中采用了4套網(wǎng)格，網(wǎng)格數(shù)分別為：19884（C）、20930（NM）、33036（M）、51410（F），分別進行了計算。研究表明，只有正確求解了來流邊界層和孔口剪切層以及孔中的流譜，才能獲得對壓力譜和聲諧振的準確預報。M網(wǎng)格和F網(wǎng)格計算得到的幅值與試驗值接近程度超過前兩套網(wǎng)格，特別是對于F網(wǎng)格，其預報的第一階模態(tài)的幅值與頻率都與試驗十分吻合。

Larcheveque（2003）[15]采用大渦模擬方法，研究了M219孔腔的流動與聲輻射，網(wǎng)格數(shù)為3.2×106。與試驗結(jié)果相比，其預報的聲諧振峰幅值的誤差分別為-3dB、-14dB、2dB、-4dB，只有第四階模態(tài)頻率的計算誤差在5%之內(nèi)。數(shù)值模擬計算表明最重要的諧振峰為第三階模態(tài)，而試驗測量表明第二階模態(tài)是最主要的。

Mendonca（2003，2004）[16-17]使用離散渦模擬（DES）方法研究了 M 219 孔腔流動與發(fā)聲，計算中采用了兩套網(wǎng)格1.1×106，2.8×106。從脈動壓力均方根的計算結(jié)果來看，精細的網(wǎng)格并沒有明顯提高計算精度。Ashworth（2005）[18]比較了非定常RANS（URANS）方法與DES方法在計算M219孔腔繞流時的差異，網(wǎng)格數(shù)為1.68×106。計算結(jié)果表明，URANS丟失了第二階模態(tài)，但是DES方法捕捉到了這一模態(tài)。

Lai等人（2006）[19]針對M 219孔腔，采用大渦模擬和FW-H聲類比技術(shù)研究了抑制流激噪聲的滲流壁面方法。來流馬赫數(shù)M＝0.85，雷諾數(shù)Re＝1.36×106。網(wǎng)格數(shù)為441萬。Lai等人在孔腔底部和后壁上采用了滲流介質(zhì)（porous media），流入孔腔的流體可以通過滲流墻流入或流出，從而孔腔之內(nèi)高度非定常的流動得到調(diào)節(jié)，渦強顯著減小。壁面滲流減小了流動與壁面的撞擊以及波反射，改變了流動振蕩的相位，有效降低了近場模態(tài)振蕩的幅值，因而可以明顯降低孔腔內(nèi)的各種脈動，從而使輻射噪聲減小。

Chen（2007）[20]采用CFD數(shù)值模擬結(jié)合氣動聲學預報方法研究了三維M219孔腔流動的近場與遠場流/聲特性。研究者旨在通過商用代碼SBLI和SotonCAA來發(fā)展一種有效的數(shù)值方法來完成從近場流動到遠場噪聲的完整的數(shù)值模擬。對于M219孔腔的流動和噪聲進行了大渦模擬（LES）和離散渦模擬（DES）并與試驗結(jié)果進行了對比。二維LES算例的網(wǎng)格數(shù)為7.2×104，三維LES算例的網(wǎng)格數(shù)為5.74×106。 二維DES算例的網(wǎng)格數(shù)為8.09×104與 8.17×104，三維DES算例的網(wǎng)格數(shù)為 1.05×106與4.08×106。與試驗結(jié)果相比，可以看出三個特點：（1）LES預報精度高于DES；（2）三維算例預報精度高于二維算例；（3）數(shù)值計算低估了頻率而高估了噪聲。

本文采用大渦模擬方法求解非定常流動，進而計算出流體動力聲源。在此基礎(chǔ)上，采用基于Kirchhoff積分（源于電磁場求解的Kirchhoff-Helmholtz定理）的FW-H聲學類比方法計算出遠場輻射噪聲?；贙irchhoff積分的FW-H方法可以克服傳統(tǒng)聲學類比方法計算體積分的困難，具有很高的靈活性和實用性。目前，在國際氣動聲學界此方法得到一定程度的運用，而在國內(nèi)，特別是水動力學界，此方法的應用還非常稀少。本文即采用上述方法求解了水中孔腔的流動輻射噪聲。

2 流激噪聲計算方法

2.1 大渦模擬（LES）方法

通過將非定常的NS方程進行濾波，得到大渦模擬的控制方程。濾波過程有效地過濾掉了那些尺度小于濾波寬度（或網(wǎng)格尺度）的小渦。

濾波函數(shù)G（x,x′）取為：

則濾波后的連續(xù)性方程和NS方程可以表示為：

其中，σij為分子粘性引起的應力張量，τij為亞格子應力

本文采用動態(tài)Smagorinsky模型，其推導過程簡述如下：

Germano（1991）在推導動力模型的過程中定義了兩種濾波算子。

另一個是測試濾波算子（test filter operator），假設(shè)測試濾波寬度比網(wǎng)格濾波大，也就是說測試濾波相當于較粗網(wǎng)格的濾波。Germano等人經(jīng)研究槽道湍流發(fā)現(xiàn)較優(yōu)的測試濾波寬度與網(wǎng)格濾波寬度的比為 2，亦即

考慮可解的湍流應力Lij。Lij表征介于網(wǎng)格濾波寬度與測試濾寬度之間的雷諾應力的貢獻，亦即小的可解尺度的貢獻。

（14）式就相當于在最小二乘意義下的最優(yōu)值。這個辦法提供了一個自洽（self－consistent）的動態(tài)模型，直接通過大渦模擬本身的計算結(jié)果，可以在每個空間網(wǎng)格點、在每一時間步上計算模型參數(shù)C。

2.2 Ffowcs Williams-Hawkings（FW-H）方程

Ffowcs Williams與Hawkings（1969）對于Lighthill理論進行了推廣，考慮了流體中任意移動壁面的發(fā)聲問題。他們使用了廣義函數(shù)的數(shù)學理論，依照Lighthill的推導步驟，將連續(xù)性方程和動量方程整理為含有兩個面源和一個體源的有源波動方程，即FW-H方程：

其中，ui：xi方向的流體速度分量；un：垂直于物體表面f＝（0）的流體速度分量；vi：xi方向的物體表面速度分量；vn：物體表面法向速度分量；δ（f）：Dirac delta函數(shù)；H（f）：Heaviside階躍函數(shù)。

p′為遠場聲壓。f=0表示物體表面，f>0表示外部無界的自由空間。ni為壁面外法線，指向流體內(nèi)部。c0為遠場聲速，ρ0為遠場密度，Tij為Lighthill應力張量。

利用自由空間格林函數(shù)δ（g）/4πr，再結(jié)合Kirchhoff可滲透表面積分法可以得到遠場解。遠場輻射聲壓中包含厚度噪聲pT′（x,t），負荷噪聲pL′（x,t）以及四極子噪聲pQ′（x,t）。

詳細的公式推導以及公式中符號含義和公式表達的物理意義見文獻[21]。

2.3 數(shù)值計算方法

時間項采用二階隱式格式離散，動量方程采用有界中心差分格式離散，壓力速度耦合采用PISO算法。計算中時間步長Δt=1×10-5s。壁面y+≈1～5。采用代數(shù)多重網(wǎng)格法加速收斂。

3 計算結(jié)果分析

3.1 計算模型

本文計算了流經(jīng)孔腔的非定常流動以及孔腔的流動輻射噪聲?？浊婚_設(shè)在載體上，載體長4m，寬0.8m，厚0.3m，首尾及兩側(cè)光順過渡?？浊痪幪?、幾何特征以及在算例中所用的網(wǎng)格數(shù)見表1。五個孔腔在載體上的位置以及相互之間的大小關(guān)系見圖1。計算中水速取V=5m/s，與試驗保持一致，基于載體長度的雷諾數(shù)為Re=2×107。試驗中每次只開設(shè)一個孔腔。

表1 孔腔編號與尺寸Tab.1 Serial number and size of cavity

3.2 非定常流動計算結(jié)果分析

Colonius（2001）[22]指出三維孔腔繞流與發(fā)聲問題的復雜程度遠遠超過二維孔腔。本文所研究的正是三維孔腔，下面的計算結(jié)果足以說明這一問題的復雜程度。流動充分發(fā)展之后，某時刻下孔腔內(nèi)的流譜見圖2～6，圖中左側(cè)為孔腔中縱剖面內(nèi)的流線分布，右側(cè)為孔腔中的三維流線。載體開設(shè)4#孔腔后表面渦量分布見圖7。

從圖2～6可知，孔腔內(nèi)的流動是十分復雜的非定常三維渦旋流動，孔內(nèi)充滿了多種形式的旋渦，它們的翻卷和擾動影響十分顯著。在孔口處，有著明顯的內(nèi)外流動交換。這些旋渦產(chǎn)生、輸運并相互干擾，構(gòu)成了主要的發(fā)聲根源。從圖7可知，孔腔極大地影響了后方及其周圍的渦量分布，其影響范圍和增大與減小的趨勢十分復雜。我們下面分析一下引起遠場輻射噪聲的主要成因。

Rossiter（1964）將孔腔中的自持振蕩（self-sustaining oscillation、self-controlled oscillation）歸結(jié)為孔腔流動發(fā)聲的根源。孔穴自持振蕩可以分為三種類型：流體—動力（fluid-dynamic）振蕩、流體—共振（fluid-resonant）振蕩；流體—彈性（fluid-elastic）振蕩。下面我們就上述三種類型逐一分析，找出本研究中孔腔輻射噪聲的來源。

首先，計算中采用的是剛性壁面，孔腔沒有變形，也不存在結(jié)構(gòu)振動。試驗模型也進行了加固和剛性處理，通過加速度傳感器測量分析，可知也不存在結(jié)構(gòu)振動。所以，那種由于結(jié)構(gòu)變形和結(jié)構(gòu)振動而引起的流體—彈性振蕩在本研究中沒有出現(xiàn)，因而流體—彈性振蕩的影響可以不予考慮。

其次，我們知道，當流體動力振蕩與孔腔的某階聲駐波共振的固有頻率相吻合時，將發(fā)出極強的單頻音（tone），產(chǎn)生聲駐波耦合共振，這種流/聲共振將產(chǎn)生強烈的輻射噪聲。Rossiter（1964）給出了辨識流體動力模態(tài)振蕩頻率的半經(jīng)驗公式：

其中，U0為來流速度，L為孔腔沿流向的長度，M為馬赫數(shù)，n為模態(tài)階數(shù)，fn為模態(tài)頻率。

第一階孔腔聲駐波共振的固有頻率為：

其中，C為環(huán)境聲速，L為孔腔沿流向的長度，H為孔腔深度。

要注意的是水中聲速約為1500m/s，水中航行體的運動都是低馬赫數(shù)或極低馬赫數(shù)下的運動（一般M≈0.001～0.01），再考慮到水中航行體上常見孔腔的尺寸，可知上述流體動力振蕩的模態(tài)頻率和聲駐波的固有頻率之間相差很遠，一般而言，不太可能出現(xiàn)上述的耦合共振，以2#孔腔為例，則可以算得第一階fA=14.2Hz，fB=1421.8Hz，二者在數(shù)量級上相差甚遠，其他各階模態(tài)也有很大差異，不太可能出現(xiàn)耦合共振。無論是從計算得到的孔腔輻射噪聲頻譜上，還是從試驗得到的孔腔輻射噪聲頻譜上，都沒有出現(xiàn)聲諧振尖峰，說明耦合共振現(xiàn)象不存在，因而流體-共振振蕩的影響可以不予考慮。

綜上所述，可以肯定地講，本研究中孔腔發(fā)出的輻射噪聲完全是由于流體動力振蕩引起的?？浊婚_口處存在自由湍流剪切層，孔腔周圍為不連續(xù)的壁面湍流邊界層，孔腔內(nèi)為非定常的渦旋運動，它們?nèi)呦嗷プ饔茫餐M成輻射噪聲的聲源。歸根到底，還在于三者所包含的脈動剪應力、脈動壓力與脈動速度，它們構(gòu)成了最根本的偶極子聲源與四極子聲源。

3.3 流動輻射噪聲計算結(jié)果分析

本節(jié)分析了計算得到的孔腔流激輻射噪聲。在試驗中，單水聽器位于開孔部位的正下方。文中所使用的試驗結(jié)果已經(jīng)過自由場修正和距離修正，測量到的流動輻射噪聲聲級統(tǒng)一進行歸一化處理，即測量距離為1m，分析帶寬1Hz，參考聲壓Pref=1μPa的聲源級。詳細的處理過程和結(jié)果見文獻[23]，本文不再贅述。計算與試驗狀態(tài)保持一致。計算與試驗結(jié)果的比較見圖8～12。

從圖8～12可知，計算反映出的孔腔輻射噪聲頻譜趨勢與量級與試驗吻合很好，各孔腔發(fā)聲的大小關(guān)系與試驗略有差異。在2000～10kHz孔腔主要發(fā)聲頻段，1#孔聲壓譜級計算誤差為-4.4～1.7dB；2#孔聲壓譜級計算誤差為-0.6～3.3dB；3#孔聲壓譜級計算誤差為-0.9～4.4dB；4#孔聲壓譜級計算誤差為-0.8～4.6dB；5#孔聲壓譜級計算誤差為-3.5～0.8dB?？梢娫诖祟l率范圍內(nèi)計算與試驗結(jié)果十分吻合，表明建立的計算方法是可靠的，可以用于復雜流動結(jié)構(gòu)誘發(fā)輻射噪聲的預報研究。在較低頻段，計算誤差增大，但還在可接受的范圍內(nèi)。

4 結(jié) 語

流聲耦合是流體力學與聲學的交叉學科，該領(lǐng)域的研究需要深厚的數(shù)學基礎(chǔ)、繁瑣的數(shù)值計算和精細的試驗驗證，其難度是不言而喻的。孔腔流動與流激噪聲的計算與試驗研究是流聲耦合領(lǐng)域中的重要研究課題，其深刻的流/聲內(nèi)涵以及計算與試驗上的難度已經(jīng)得到國際上的普遍公認，具有相當?shù)募夹g(shù)難度。本文在前期工作的基礎(chǔ)上，首先利用大渦模擬方法預報了三維孔腔在水中繞流的流動特征，然后利用FW-H聲學類比方法預報由此產(chǎn)生的遠場輻射噪聲，分析了噪聲產(chǎn)生的流動根源，并與相應的試驗結(jié)果進行了對比，證明本文方法的計算流程切實可行，計算量級可靠，可用于水中流激噪聲的預報分析問題，并應在將來的工作中做進一步的分析和驗證。

[1]Germano M,Piomelli U,Moin P,Cabot W H.A dynamic subgrid-scale eddy viscosity model[J].Phys.Fluids,1991,A3(7):1760-1765.

[2]Lilly D K.A proposed modification of the Germano Subgrid-scale closure method[J].Phys.Fluids,1992,A4(3):633-635.

[3]Seror C,Sagaut P,Bailly C,Juv’e D.Subgrid-scale contribution to noise production in decaying isotropic turbulence[J].AIAA J,2000，38(10):1795-1803.

[4]Seror C,Sagaut P,Bailly C,Juv’e D.On the radiated noise computed by large eddy simulation[J].Phys.Fluids,2001,13(2):476-487.

[5]He G W,Rubinstein R,Wang L P.Effects of subgrid scale modeling on time correlations in large eddy simulation[J].Phys.Fluids,2002,14(7):2186-2193.

[6]He G W,Wang M,Lele S K.On the computation of space-time correlations by large-eddy simulation[J].Phys.Fluids,2004,16(11):3859-3867.

[7]Wang Meng,Freund J B,Lele S K.Computational prediction of flow-generated sound[J].Annu.Rev.Fluid Mech.,2006,38:483-512.

[8]Wang M,Moin P.Computation of trailing-edge flow and noise using large eddy simulation[J].AIAA J,2000,38(12):2201-2209.

[9]Wang M,Moreau S,Iaccarino G,Roger M.LES prediction of pressure fluctuations on a low-speed airfoil.Annu[Z].Res.Briefs,Cent.Turbul.Res.,Stanford Univ./NASA Ames,Stanford,CA,2004:183-193.

[10]Roshko A.Some measurements of flows in a rectangular cutout[K].NACA TN 3488,1955.

[11]Rossiter J E.Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[R].Rep.Mem.3438,Aeronautical Research Council,1964.

[12]Jonas Ask,Lars Davidson,Hans Enwald,Johan Larsson.An acoustic analogy applied to incompressible flow fields[C]//Computational Aeroacoustics:From Acoustic Sources Modeling to Far-Field Radiated Noise Prediction Colloquium EUROMECH 449,December 9-12,2003.Chamonix,France,2003.

[13]Loh Ching Y.Computation of low speed cavity noise[C]//42nd Aerospace Sciences Meeting and Exhibit sponsored by the American Institute of Aeronautics and Astronautics Reno.Nevada,January 5-8,2004.(AIAA-2004-0680)

[14]Peng Shiahui.Unsteady RANS simulation of turbulent cavity flow:Summary of 2D baseline computations[R].Swedish Defence Research Agency.FOI-R-1915-SE,Technical report,December 2005.

[15]Larcheveque L,Sagaut P,Le T H.Large-eddy simulations of flows in weapon bays[G].AIAA Paper 2003-0778,2003.

[16]Mendonca F,Allen R,Charentenay J,Kirkham D.CFD prediction of narrowband cavity acoustics at M=0.85[G].AIAA Paper 2003-3303,2003.

[17]Allen R,Mendonca F.DES validations of cavity acoustics over the subsonic to supersonic range[G].AIAA Paper 2004-2862,2004.

[18]Ashworth R M.Prediction of acoustic resonance phenomena for weapon bays using detached eddy simulation[J].Aeronautical Journal,2005:631-638.

[19]Lai H,Luo K H.Large-eddy simulation and control of cavity aeroacoustics[C]//Conference on Turbulence and Interactions TI2006,May 29-June 2,2006.Porquerolles,France,2006.

[20]Chen Xiaoxian,Sandham Neil D,Zhang Xin.Cavity Flow Noise Predictions[R].University of Southampton,Report No.AFM-07/05,2007.

[21]Brentner K S.An efficient and robust method for prediction helicopter high-speed impulsive noise[J].JSV,1997,203(1):87-100.

[22]Tim Colonius.An overview of simulation,modeling,and active control of flow/Acoustic resonance in open cavities[G].AIAA 2001-0076,2001.

[23]熊紫英,孫紅星.孔腔模型流動噪聲測量試驗報告[R].無錫:中國船舶科學研究中心科技報告，2009.