局部開孔深腔流噪聲發(fā)聲機理研究

李榮華，樓京俊，朱石堅

(1. 海軍工程大學(xué) 動力工程學(xué)院，湖北 武漢 430033；

2. 海軍工程大學(xué) 船舶振動噪聲重點實驗室，湖北 武漢 430033)

0 引 言

通風(fēng)孔等功能孔是水下航行器不可少的結(jié)構(gòu)，隨著水下航行器速度的提升，由流水孔孔腔振蕩產(chǎn)生的流噪聲成為不可忽視的因素。

在國內(nèi)外已經(jīng)有多位學(xué)者對孔腔流噪聲展開了研究，首先是在氣動聲學(xué)領(lǐng)域，Lighthill[1]推導(dǎo)了聲類比方程，將流體噪聲源等效為緊致極子聲源，開創(chuàng)了聲類比方法；Powell[2]提出渦聲理論，指出有渦結(jié)構(gòu)的地方才有聲音；Howe[3]指出在極低馬赫數(shù)下，孔腔氣動流噪聲的主要聲源是由非定常阻力產(chǎn)生的偶極子聲源；楊黨國[4]采用CFD計算和氣動聲學(xué)理論在來流馬赫數(shù)為0.64條件下對長深比為2:1的空腔進行自激振蕩發(fā)聲機理研究，提出了腔內(nèi)渦脫落-一次聲波-新的渦脫落-二次聲波的反饋機制發(fā)聲機理。隨著潛艇等水中航行器的發(fā)展，水動力聲學(xué)也逐漸得到了重視和研究，對于水動力聲學(xué)的研究在數(shù)值仿真方面也采用了與氣動聲學(xué)仿真類似的方法。Charib[5]采用實驗手段研究開口腔的水流激發(fā)振蕩問題，指出了邊界層動量損失厚度與流體開口腔振蕩模式的關(guān)系；張楠[6]中指出水中孔腔渦旋流場誘導(dǎo)輻射噪聲的主要聲源是偶極子聲源；王玉[7]提出陷落腔激發(fā)的渦流噪聲包括低頻純音和寬帶連續(xù)噪聲，其中低頻純音主要由大尺度渦運動引起，聲功率較高。高頻輻射聲呈現(xiàn)寬帶特性，是由小尺度渦運動變化的結(jié)果；徐俊[8]通過數(shù)值模擬得出空腔噪聲的4個聲源分別為：腔口剪切層小渦系列、后緣尾渦流系列、腔內(nèi)大渦以及空腔后緣的渦流撞擊與破碎，并且針對其分析結(jié)果給出了在孔腔內(nèi)添加隔板的降噪方案；高巖等[9]針對全開孔式孔腔模型對其聲模態(tài)進行研究，對比了剛性壁以及彈性壁條件下的聲耦合狀況；陳燦[10]通過大渦模擬-聲類比混合方法進行數(shù)值模擬，得出來流速度的增大會加劇空腔的振蕩，甚至?xí)箍浊挥刹徽袷幾優(yōu)檎袷?，但是隨著孔腔長深比的減小，振蕩會逐漸減弱，說明潛艇流水孔的深腔結(jié)構(gòu)設(shè)計在流場方面較為合理，因為深腔的振蕩會比淺腔振蕩更弱；Yuan Guo-qing[11]采用大渦模擬-聲類比混合方法對局部開孔的腔體模型進行仿真，對淺腔局部開孔腔體的流噪聲激勵機理進行研究，得到自持振蕩反饋環(huán)是造成孔腔流噪聲的原因。

楊黨國、王玉、徐俊、陳燦等采用的計算模型均為圖1（a）所示的全開孔孔腔模型，Yuan Guo-qing采用的模型則為淺腔局部開孔模型（長深比6:1），與水下航行器潛浮系統(tǒng)的實際模型中采用的局部開孔深腔模型均有較大差異。本文采用大渦模擬-聲類比混合方法對圖1（b）所示的局部開孔深腔模型的孔腔流噪聲發(fā)聲機理進行研究，根據(jù)流場渦量、壓強變化以及孔后壁檢測點脈動壓力的周期、相位異同，結(jié)合空腔內(nèi)聲學(xué)模態(tài)，對局部開孔深腔模型（長深比＜1）的孔腔流噪聲的發(fā)生機理進行研究。

圖1 孔腔模型Fig. 1 The model of open cavity

1 數(shù)值求解方法

1.1 大渦模擬

湍流N-S方程由于方程的不封閉，因此不能進行解析求解，所以對于N-S方程的求解都采用了數(shù)值計算方法。直接數(shù)值模擬方法（DNS）、雷諾平均NS方程方法（RNS）、大渦模擬方法（LES）是目前研究湍流和工程應(yīng)用的3種方法，其中LES綜合考慮了計算機資源和計算精度，對湍流進行過濾，分離湍流中的大渦和小渦，其中大尺度渦在湍流中起主要作用，而小尺度渦則采用特殊方法進行封閉。

大渦模擬采用過濾方法來消除湍流中小尺度脈動，在物理空間中，過濾過程采用積分運算來實現(xiàn)。積分過程可用如下公式表示：

其中： ?為積分立方體的邊長，其大小稱為過濾長度； ui(ξ,t)為樣本瞬時速度； uˉi(x,t)為過濾后的大尺度速度； G(x-ξ)為過濾函數(shù)，其表達式如下：

經(jīng)過過濾函數(shù)之后，湍流速度場或其他湍流流量過濾后仍然為不規(guī)則量，只是這些不規(guī)則量中小尺度脈動已經(jīng)過過濾，只剩下尺度大于的脈動。經(jīng)過濾后的N-S方程（大渦模擬控制方程）如下：

其中： i,j=1,2,3； ui為與 xi關(guān)聯(lián)的速度分量；為過濾后的平均速度分量； ρ為流體密度； ν為流體運動粘性系數(shù)； τij為亞格子應(yīng)力

1.2 FW-H聲類比

FW-H方程本質(zhì)上是非齊次方程，其可以從連續(xù)方程以及N-S方程中演化得來，F(xiàn)W-H方程如下：

其中： ui為在 xi方向上的流體速度分量； un為垂直于平面的流體速度分量；為在方向上的物體速度分量；為垂直于物面方向的物體速度分量；為狄拉克函數(shù)；為階躍函數(shù)；為遠場聲壓。為物體表面，物體外部無界空間；為遠場聲速，為遠場密度，為Lighthill應(yīng)力張量：

其中：

對于本文中實例而言，有vi=0。

1.3 計算模型與CFD計算參數(shù)設(shè)置

結(jié)合水下航行器實際情況，將計算模型等效為帶有局部開孔的深腔模型，如圖2所示。長寬深分別為L=400 mm，W=300 mm，H=500 mm，其中開孔長寬Lhole=300 mm，Whole=200 mm，孔頸高度為Hhole=20 mm，按照Sarohia[12]所提出的劃分標(biāo)準(zhǔn)，該模型為深腔模型（長深比＜1）。圖3所示為深腔模型計算區(qū)域，其中計算區(qū)域模型尺寸為5 000 mm×1 800 mm×1 500 mm，開口中心距計算域前端及兩側(cè)的距離分別為2 000 mm和900 mm。采用SolidWorks進行建模，導(dǎo)入ICEM CFD中進行網(wǎng)格劃分。整個計算模型的網(wǎng)格總數(shù)為1 676 229，第1層網(wǎng)格距離為0.03，全域采用結(jié)構(gòu)化網(wǎng)格離散，為了在保證仿真準(zhǔn)確性同時減少計算量，如圖4所示，對邊界層網(wǎng)格、腔體開口處進行加密處理。

邊界條件設(shè)置如圖4所示。

速度入口：設(shè)定來流的大小為10 m/s，方向與X軸方向平行；

壓力出口：設(shè)定相對于參考壓力為0；

壁面：腔體壁面設(shè)置為無滑移光滑壁面；除腔體壁面之外的其余壁面為自由滑移壁面。

圖2 局部開口深腔示意圖Fig. 2 The model of deep cavity with local hole

圖3 計算域示意圖Fig. 3 Figure of calculation zone

圖4 網(wǎng)格及邊界條件示意圖Fig. 4 The flow field meshes and boundary conditions

2 計算分析

在CFD軟件Fluent中，先選用大渦模擬模塊進行計算，待計算穩(wěn)定后，打開FW-H聲計算模塊，設(shè)置局部開口腔體為聲源，進行聲場計算。

為方便結(jié)果分析，在計算過程中布置多個監(jiān)測點對流場中的孔壁、腔壁壓強的變化情況進行監(jiān)測，在開口中心正上方處布置一個接收點receiver1，坐標(biāo)為[0，1 020，0]，對聲場中聲壓進行監(jiān)測，并做圖2所示切面xOy，對仿真過程中切面xOy上的流場渦量、壓強變化情況進行分析。

2.1 流場分析

如圖5左側(cè)所示，1，2，3點分別為圖2中所示孔后壁的垂向中線處的點，監(jiān)測點1，2，3隨時間變化的壓強脈動信息變化如圖5右側(cè)所示，各監(jiān)測點壓力變化周期相同，但是其相位不同，其中1點和3點相位差值接近π。

圖5 監(jiān)測點壓強時間序列Fig. 5 Pressure-time diagram of measurement points

當(dāng)來流流經(jīng)開口前緣時，由于開口的存在導(dǎo)致來流在開口前緣處發(fā)生邊界層分離，不斷產(chǎn)生一系列分離渦，在來流帶動下不斷向孔后壁方向運動遷移，來流在運動過程中與壁面發(fā)生碰撞，從而引起壁面壓強的變化。取點3壓力達到最小值的tS=0.02時刻為周期起點t=0，每個周期時長用T表示，結(jié)合時域內(nèi)流體運動、聲場內(nèi)聲波傳播、聲壓變化等因素，對腔體開孔處的渦量、壓強周期性變化及不同監(jiān)測點之間的脈動周期相同、相位不同的原因進行分析。

t=0時刻xOy切面上流場渦量-壓強分布云圖如圖6所示。在t=0時刻，此時孔壁上緣3點附近渦量分布集中，下緣1點處渦量分布較少。由渦運動相關(guān)知識可知，在漩渦中，由于離心力被壓力所平衡，渦核處離心力最低，壓力也最小，所以此時在開口上緣區(qū)域壓強分布較小，下緣區(qū)域壓強分布較大。但是由于渦運動對開口上緣附近壁面產(chǎn)生較大沖擊，有較大脈動壓力，等效為偶極子聲源，此時在孔后壁上緣處所產(chǎn)生的聲壓也較大，隨著聲波向孔前壁方向傳播，對渦結(jié)構(gòu)運動產(chǎn)生排擠，促使后續(xù)向孔后壁遷移的渦結(jié)構(gòu)同時向開口下緣移動。

圖6 t=0時刻渦量-壓強云圖Fig. 6 Vorticity-pressure contours of t=0

t=0.25 T時刻xOy切面上流場渦量-壓強分布云圖如圖7所示。在上一時刻腔體后壁上緣聲源所產(chǎn)生聲壓的作用下，渦結(jié)構(gòu)向孔后壁移動的同時，逐漸向下緣移動，且隨著渦的不斷耗散和向后移動，在孔后壁上緣的渦量得不到補充，因此上緣渦量分布減少，而下緣的渦量分布增多。在渦的作用下，下緣1點處附近的壓強減小，上緣3點處附近的壓強增大，雖然此時3點附近下緣區(qū)域的渦量分布增大，對下緣作用產(chǎn)生的脈動壓力也增大，但是由此聲源產(chǎn)生的聲壓仍小于由上緣壁面產(chǎn)生的聲壓，因此在壁面上緣聲源產(chǎn)生的大聲壓作用下，向孔后壁遷移的渦結(jié)構(gòu)繼續(xù)向開口下緣移動。

圖7 t=0.25 T時刻渦量-壓強云圖Fig. 7 Vorticity-pressure contours of t=0.25 T

t=0.4 T時刻xOy切面上渦量-壓強分布云圖如圖8所示。隨著孔壁上緣渦的不斷耗散，同時由孔前壁向后壁移動的渦向下緣移動，在孔后壁上緣3點處的渦量持續(xù)得不到補充，因此上緣3點處的渦量分布達到周期內(nèi)的最小值，周邊壓強較大。而下緣1點處由于渦的持續(xù)補充，此時渦量達到較大值，周邊壓強較小。同理，由于此時在孔后壁下緣處的渦量較大，渦結(jié)構(gòu)運動過程中對開口下緣腔壁造成沖擊，產(chǎn)生較大脈動壓力，從而作為聲源產(chǎn)生較大聲壓，隨著聲波向開口前緣傳播，對向后運動的渦結(jié)構(gòu)產(chǎn)生排擠，導(dǎo)致向孔后壁遷移的渦結(jié)構(gòu)逐漸向開口上緣移動。

圖8 t=0.4 T時刻渦量-壓強云圖Fig. 8 Vorticity-pressure contours of t=0.5 T

圖9 t=0.75 T時刻渦量-壓強云圖Fig. 9 Vorticity-pressure contours of t=0.75 T

t=0.75 T時刻xOy切面上渦量-壓強分布云圖如圖9所示。由于孔后壁腔體下緣發(fā)聲產(chǎn)生的聲壓隨著聲波向孔前壁方向傳播對渦結(jié)構(gòu)產(chǎn)生排擠，因此從孔前壁向后壁遷移的渦逐漸向開口上緣移動，1點處的渦量得不到補充，分布減少，3點處的渦量分布增多，因此1點附近壓強增大，3點周邊壓強減小，此時渦在孔壁下緣聲源大聲壓作用下，繼續(xù)向開口上緣移動。

t=T時刻xOy切面上渦量-壓強分布云圖如圖10所示。此時在腔體下緣的渦量分布達到較小值，1點附近壓強達到較大值，而在開口上緣3點出的渦量分布達到最大值，附近壓強達到最小值，與t=0時刻一致。此時由于上緣渦量分布集中，渦結(jié)構(gòu)運動產(chǎn)生的壓力脈動在腔體后緣作為聲源產(chǎn)生較大聲壓，隨著聲波向前傳播，對渦結(jié)構(gòu)移動產(chǎn)生排擠，致使由開口前緣向孔后壁遷移的渦向后遷移的同時向開口下緣移動，從而進入下一個渦結(jié)構(gòu)運動周期。

圖10 t=1 T時刻渦量-壓強云圖Fig. 10 Vorticity-pressure contours of t=T

從以上分析可得，流體在流經(jīng)開口時，由渦上（下）運動—碰撞發(fā)聲—渦下（上）運動—碰撞發(fā)聲—渦上（下）運動整個過程構(gòu)成的渦-聲反饋是造成腔體脈動壓力周期性變化及不同相位的主要原因。

2.2 受力分析

從圖5可以看出，各個監(jiān)測點之間的周期一致均為T=0.057 5 s，即期頻率為 f測=1/T=1/0.057 5 s=17.4 Hz。圖11為各監(jiān)測點的壓強功率譜曲線，其中點1，2，3位于圖2所示孔后壁，點4，5，6位于孔前壁，點7，8，9位于孔側(cè)壁，點10，11和12，13分別位腔后壁和腔前壁。從圖中可以得出各個監(jiān)測點擁有相同頻率的特征線譜，1階特征線譜頻率為Hz，2階特征線譜頻率為Hz，說明其所受來自流體的壓力頻率一致。但是孔后壁的1階、2階特征線譜峰值遠大于開孔前壁、孔側(cè)壁及開口下方腔體壁的壓強功率譜峰值，在孔后壁的上緣點3功率譜密度遠高于其他采集點的功率譜密度，表明孔壁后緣壓力脈動是聲場能量的主要來源，其中以上緣對其影響最大，即孔后壁上緣是開口自持振蕩發(fā)聲的主要聲源。

2.3 聲場分析

圖12（a）為開口正上方1 000 mm（2倍腔體深度）處接收點receiver1的聲壓頻譜特性曲線，截取其中的低頻成分進行分析，如圖12（b）所示，該聲壓頻譜的前2階頻率分別為Hz，34.96 Hz，與各采集點壓強功率譜的前2階特征線譜頻率一致，說明由流體產(chǎn)生的壓力脈動為其主要聲源。

圖11 各采集點壓力功率譜Fig. 11 Exciting force power spectral density of measurement points

由于該腔體為局部開孔腔體，且為深腔，滿足聲波模態(tài)出現(xiàn)的條件[13]，因此需要考慮某一特定頻率下出現(xiàn)孔腔聲模態(tài)共振現(xiàn)象所導(dǎo)致相應(yīng)頻率上強烈的波動對測量結(jié)果的影響。在孔腔內(nèi)聲場模態(tài)的簡正頻率滿足：

圖12 Receiver1聲壓頻譜特性Fig. 12 Sound pressure spectral of receiver1

式中 ： n1=0,1,2,3··； n2=0,1,2,3 ··;n3=0,1,2,3 ··；對其聲波模態(tài)進行求解，得到該腔體最低階聲波模態(tài)頻率為1 500 Hz，遠大于接收點所測得的各低階模態(tài)的頻率。對比圖12（a）可得，在1 500 Hz沒有出現(xiàn)特征線譜，因而可以排除腔體聲耦合共振對接收點receiver1所測各階模態(tài)的影響。

Rossiter JE[14]最早對矩形開口腔體繞流振蕩機理進行研究，并給出自持振蕩頻率預(yù)測半經(jīng)驗公式：

其中： α=0.25； kc=0.57； U∞=10 m/s； Lhole=0.3 m。根據(jù)本文實例計算可得 f1=14.29 Hz，仿真計算得到的fs1=17.41Hz，與經(jīng)驗公式的預(yù)測結(jié)果基本吻合，驗證了本文仿真結(jié)果的準(zhǔn)確性。

3 結(jié) 語

本文通過大渦模擬-聲類比混合方法，進行數(shù)值仿真，對局部開孔深腔體的流噪聲產(chǎn)生機理進行分析。

1）根據(jù)仿真得到一個周期內(nèi)空腔內(nèi)部各個時刻的渦量、壓強分布云圖，首次對孔后壁不同監(jiān)測點脈動壓力變化的周期、相位異同之處進行監(jiān)測分析，用渦上（下）運動—碰撞發(fā)聲—渦下（上）運動—碰撞發(fā)聲—渦上（下）運動的渦-聲反饋模型解釋孔腔流噪聲發(fā)聲機理。

2）對各個采集點的壓強功率譜密度進行對比分析，得出局部開孔深腔的各個部位對聲場能量的貢獻量不同，其中開孔后壁的貢獻量最大。

3）將計算結(jié)果與Rossiter經(jīng)驗公式的計算結(jié)果進行對比，驗證了仿真結(jié)果的準(zhǔn)確性。