基于聲擾動(dòng)方程的氣動(dòng)噪聲傳播積分計(jì)算方法

鄭雯斯，王芳

1.西北工業(yè)大學(xué)，陜西 西安 710072

2.北方民族大學(xué)，寧夏 銀川 750030

以FW-H 方程[1]為代表的經(jīng)典聲比擬理論假設(shè)介質(zhì)整場均勻，建立了聲壓與不同類型聲源之間的響應(yīng)關(guān)系，成為工程應(yīng)用中最受歡迎的氣動(dòng)噪聲預(yù)測方法[2]。但流動(dòng)在聲源區(qū)域通常存在梯度，對近場聲傳播起決定性作用的是聲源相對于緊臨周圍（約一個(gè)波長左右范圍）的非均勻平均流運(yùn)動(dòng)，而不是相對于無窮遠(yuǎn)處的均勻介質(zhì)運(yùn)動(dòng)。為了將非均勻平均流的作用包括在聲比擬理論中，需建立一個(gè)恰當(dāng)?shù)牟▌?dòng)方程來描述聲波產(chǎn)生與輻射過程。

文獻(xiàn)[3]～[6]對非均勻平均流中的聲傳播進(jìn)行了研究，提出了各自的波動(dòng)方程。Goldstein[7]提出廣義聲比擬理論，將N-S方程重組為對流形式的線化N-S方程組，以此分析了射流中“真實(shí)”的氣動(dòng)聲源。這些對流波動(dòng)方程深化了人們對氣動(dòng)噪聲產(chǎn)生機(jī)理的認(rèn)識，但面對復(fù)雜流動(dòng)問題，用于求解這些方程的格林函數(shù)不能獲得理論解。在實(shí)際應(yīng)用中，如何準(zhǔn)確求解格林函數(shù)成為噪聲預(yù)測中的一個(gè)關(guān)鍵。Tam[8]考慮了非均勻平均流的聲折射作用，發(fā)展了關(guān)聯(lián)近場聲源和遠(yuǎn)場聲壓的伴隨格林函數(shù)方法，計(jì)算量小且不存在積分奇點(diǎn)，目前僅應(yīng)用于固體邊界不起主要作用的氣動(dòng)噪聲問題，如Karabasov 等[9-10]采用伴隨格林函數(shù)對射流噪聲在非均勻介質(zhì)中的傳播進(jìn)行過研究。針對非均勻介質(zhì)中聲傳播，Goldstein[11]提出了線化勢函數(shù)模型，Pierce[12]發(fā)展了高頻近似模型，Bras 等[13]則應(yīng)用兩個(gè)模型分析了運(yùn)動(dòng)介質(zhì)中聲傳播的直接問題和伴隨問題。伴隨方法有益于解決非均勻介質(zhì)中聲傳播問題，Spieser 等[14]采用伴隨方法研究了Ma0.9 的湍流噪聲問題，趙雯[15]采用伴隨方法分析了矩形噴口射流噪聲。

流動(dòng)發(fā)聲與渦之間存在密切聯(lián)系，Powell[16]和Howe[17]提出的渦聲理論能從物理上清晰地解釋流體運(yùn)動(dòng)誘發(fā)氣動(dòng)噪聲的機(jī)理，認(rèn)為渦是流動(dòng)發(fā)聲的根本原因。目前，渦聲理論多應(yīng)用于低馬赫數(shù)流動(dòng)產(chǎn)生的氣動(dòng)噪聲問題，并認(rèn)為聲波是在均勻介質(zhì)中傳播，忽略了平均流與聲波間的相互作用。Ewert等[18]提出了考慮介質(zhì)非均勻效應(yīng)的聲擾動(dòng)方程，但需要利用高精度時(shí)空離散格式[19-20]和無反射邊界條件[21-22]在時(shí)域內(nèi)計(jì)算聲傳播，因而效率低。

高效、準(zhǔn)確預(yù)測聲傳播是氣動(dòng)聲學(xué)研究的一個(gè)重要方向。為了將渦聲思想應(yīng)用于非低馬赫數(shù)流動(dòng)，本文采用聲擾動(dòng)方程描述聲波運(yùn)動(dòng)，依據(jù)伴隨方法構(gòu)建關(guān)聯(lián)近場聲源與遠(yuǎn)場聲壓的頻域格林函數(shù)，從而發(fā)展一種高效、準(zhǔn)確的聲傳播積分方法，可為復(fù)雜流動(dòng)氣動(dòng)噪聲產(chǎn)生機(jī)理和復(fù)雜介質(zhì)氣動(dòng)噪聲傳播機(jī)制分析提供理論與方法支持。

1 方法描述

1.1 聲擾動(dòng)方程

忽略黏性和熱傳導(dǎo)對聲傳播的影響，N-S方程可簡寫為

式中，標(biāo)量ρ，p及矢量u，分別表示當(dāng)?shù)厮矐B(tài)流動(dòng)的密度、壓力和速度。將當(dāng)?shù)亓鲃?dòng)分解為平均流和非穩(wěn)態(tài)擾動(dòng)兩部分，即

式中，ρˉ，pˉ和uˉ表示無窮遠(yuǎn)均勻來流分量；ρ′，p′和u′表示非穩(wěn)態(tài)擾動(dòng)量。利用式（2）線化式（1），得到

1.2 伴隨格林函數(shù)法

對圖1（a）所示聲波直接輻射問題，觀察點(diǎn)和聲源位置分別用矢量x和y表示，關(guān)聯(lián)兩者的格林函數(shù)為G(x,y)。為獲取滿足式（6）的格林函數(shù)，需將圖1（a）所示的直接問題轉(zhuǎn)化為圖1（b）所示的伴隨問題，伴隨格林函數(shù)用G(y,x)表示。對線性聲學(xué)問題，根據(jù)聲學(xué)互易定理有G(x,y) =G(y,x)。用G4和Gi（二維問題i=1, 2，三維問題i=1, 2, 3）分別表示壓力型和速度型伴隨格林函數(shù)。根據(jù)式（6），G4和Gi在頻域內(nèi)滿足

圖1 直接問題與伴隨問題示意圖Fig.1 Schematic of direct and adjoint problems

式中，i和ω分別為虛數(shù)單位和圓頻率。

采用有限差分方法求解式（8）以獲得壓力型和速度型伴隨格林函數(shù)的數(shù)值解?？臻g離散采用五點(diǎn)四階中心差分格式。將聲傳播區(qū)域分為常系數(shù)和變系數(shù)波動(dòng)區(qū)域兩部分，在常系數(shù)區(qū)域使用緩沖區(qū)（Buffer Zone）邊界條件結(jié)合拉伸網(wǎng)格的方法來吸收和衰減聲波，在固體邊界施加聲學(xué)硬邊界條件。

遠(yuǎn)場聲傳播可采用下述積分方程描述

2 算例驗(yàn)證

利用二維圓柱和NACA0012翼型繞流算例來驗(yàn)證上述氣動(dòng)噪聲預(yù)測方法的可靠性。采用Fluent軟件計(jì)算非定常流動(dòng)以提取氣動(dòng)聲源，湍流模型選擇基于S—A模型的IDDES方法。

2.1 二維圓柱氣動(dòng)噪聲預(yù)測

二維圓柱直徑D=0.04m，無窮遠(yuǎn)均勻來流馬赫數(shù)為0.3，基于直徑的雷諾數(shù)Re為100。流動(dòng)遠(yuǎn)場距離圓柱中心200D，對圓柱尾渦脫落區(qū)域進(jìn)行網(wǎng)格加密處理，結(jié)構(gòu)化網(wǎng)格總數(shù)約為20 萬個(gè)，聲學(xué)積分區(qū)域位于紅色線條內(nèi)部，如圖2所示。

圖2 二維圓柱網(wǎng)格Fig.2 Configuration of mesh gird of 2D cylinder

圖3 為二維圓柱繞流展向渦量的空間分布云圖，低雷諾數(shù)層流運(yùn)動(dòng)在圓柱尾跡區(qū)域出現(xiàn)了周期性的渦脫落。流動(dòng)計(jì)算得到的尾渦脫落頻率為416Hz，相應(yīng)的無量綱頻率St=0.161。

圖3 二維圓柱瞬時(shí)展向渦量云圖Fig.3 Snapshots of the instantaneous spanwise-vorticity structure of 2D cylinder

將聲學(xué)觀察點(diǎn)置于以圓柱中心為原點(diǎn)、半徑為20D的圓上。對圓柱尾渦脫落頻率噪聲，選取120°觀察點(diǎn)，非均勻平均流和均勻平均流條件下壓力型伴隨格林函數(shù)虛部的空間分布如圖4 所示，兩者結(jié)果幾乎一致。尾渦脫落頻率波長與圓柱直徑比約為20，意味著近場非均勻流動(dòng)區(qū)域的特征尺寸遠(yuǎn)小于一個(gè)波長，這使得近場流動(dòng)非均勻性對聲傳播的影響可以忽略不計(jì)。

圖4 二維圓柱壓力型伴隨格林函數(shù)虛部云圖Fig.4 Snapshots of the imaginary part of the pressure pattern adjoint Green's function for 2D cylinder

數(shù)值預(yù)測的90°觀察點(diǎn)聲場頻譜如圖5所示，與計(jì)算流體力學(xué)（CFD）直接計(jì)算結(jié)果相比，尾渦脫落頻率St=0.161及其一階諧波St=0.322 處，兩者的噪聲大小具有很好的一致性。圖6進(jìn)一步顯示了尾渦脫落頻率的聲壓級指向性分布，各觀察點(diǎn)處的數(shù)值預(yù)測結(jié)果與CFD直接計(jì)算值非常吻合，最大誤差約為1dB。

圖5 二維圓柱聲場的聲壓級頻譜Fig.5 Spectrum of sound pressure level of 2D cylinder

圖6 二維圓柱尾渦脫落頻率聲場空間指向性分布Fig.6 Directivity of sound pressure level at the frequency of vortex shedding of 2D cylinder

2.2 二維NACA0012翼型氣動(dòng)噪聲預(yù)測

NACA0012翼型弦長C=0.1m、迎角為7°的均勻來流馬赫數(shù)Ma=0.4，基于弦長的雷諾數(shù)Re=5×104。流動(dòng)遠(yuǎn)場距離翼型中心200C，網(wǎng)格總數(shù)約為32 萬個(gè)，在流動(dòng)分離區(qū)域進(jìn)行網(wǎng)格加密處理。圖7 為瞬時(shí)展向渦量空間分布，流動(dòng)計(jì)算所得尾渦脫落頻率為St=2.92，對應(yīng)的聲波波長與翼型弦長比為0.85，波長接近弦長預(yù)示非均勻平均流對聲波傳播的影響不能忽略。

圖7 二維NACA0012翼型瞬時(shí)展向渦量云圖Fig.7 Snapshots of the instantaneous spanwise-vorticity structure of 2D NACA0012 airfoil

觀察點(diǎn)布置在以翼型流向中心為原點(diǎn)、半徑為10C的圓上。對渦脫落頻率St=2.92，選取60°觀察點(diǎn)，非均勻平均流和均勻平均流條件下壓力型伴隨格林函數(shù)幅值的空間分布如圖8 所示。兩者的空間分布基本相似，但在渦脫落區(qū)域存在明顯差異。

圖8 二維NACA0012翼型壓力型伴隨格林函數(shù)云圖Fig.8 Snapshots of the pressure pattern adjoint Green's function for 2D NACA0012 airfoil

圖9 進(jìn)一步給出了非均勻平均流和均勻平均流條件下，60°觀察點(diǎn)對應(yīng)的速度型伴隨格林函數(shù)G1實(shí)部的空間分布。類似壓力型伴隨格林函數(shù)，渦脫落區(qū)域的非均勻平均流對聲傳播存在顯著影響。

圖9 二維NACA0012翼型速度型伴隨格林函數(shù)實(shí)部云圖Fig.9 Snapshots of the real part of thevelocity pattern adjoint Green's function for 2D NACA0012 airfoil

圖10 給出了90°觀察點(diǎn)聲場頻譜數(shù)值預(yù)測結(jié)果，與CFD 直接計(jì)算結(jié)果對比，兩者在頻率分布和峰值大小上基本一致。三個(gè)主要峰值頻率（St=1.96, 2.92, 3.87）的噪聲結(jié)果對比見表1，可以看出，St=2.92 處的噪聲差異最大，聲壓級的數(shù)值預(yù)測結(jié)果比CFD直接計(jì)算值高1.3dB。圖11進(jìn)一步展示了不同頻率下的聲壓級指向性分布，在各觀察點(diǎn)處，本文方法的數(shù)值預(yù)測結(jié)果與CFD直接計(jì)算值吻合較好。

表1 二維NACA0012翼型主要峰值噪聲Table 1 Main peak noise of 2D NACA0012 airfoil

圖10 二維NACA0012翼型聲場頻譜Fig.10 Spectrum of sound field of 2D NACA0012 airfoil

圖11 二維NACA0012翼型不同頻率聲場指向性分布圖Fig.11 Directivity patterns of 2D NACA0012 airfoil noise at different frequencies

3 結(jié)論

為了考慮非均勻介質(zhì)聲折射對氣動(dòng)噪聲傳播的影響，本文采用聲擾動(dòng)方程描述聲波在非均勻介質(zhì)中的傳播，基于聲學(xué)互易原理構(gòu)建關(guān)聯(lián)近場聲源與遠(yuǎn)場聲壓的頻域伴隨格林函數(shù)。進(jìn)而利用伴隨格林函數(shù)和聲擾動(dòng)方程源項(xiàng)的數(shù)值計(jì)算結(jié)果，通過頻域積分計(jì)算噪聲傳播。通過研究可以得到以下結(jié)論：

（1） 二維圓柱和NACA0012翼型繞流噪聲的頻譜和指向性預(yù)測結(jié)果與CFD直接計(jì)算值吻合。

（2） 對聲波波長遠(yuǎn)大于非均勻流動(dòng)區(qū)域范圍的聲輻射問題，可以直接忽略非均勻流動(dòng)對聲傳播的影響；而當(dāng)聲波波長與非均勻流動(dòng)區(qū)域特征尺寸相近時(shí)，非均勻流動(dòng)對聲波的折射作用明顯。