薛絲丹，錢振昊，何嘉華，劉秋洪

（1.西北工業(yè)大學(xué)翼型葉柵空氣動力學(xué)重點實驗室，西安 710072；2.中國空氣動力研究與發(fā)展中心氣動噪聲控制重點實驗室，四川綿陽 621000）

在過去幾十年中，旋轉(zhuǎn)葉片產(chǎn)生的氣動噪聲問題一直是航空聲學(xué)的研究熱點，目前已經(jīng)在理論分析和數(shù)值計算方面取得一些重要的成果。早期理論研究，包括Ffowcs Williams 等[1]、Lowson[2]以及Morfey 等[3]的工作，都傾向于將旋轉(zhuǎn)葉片噪聲源簡化為理想的旋轉(zhuǎn)點源，以點源模型預(yù)測旋轉(zhuǎn)葉片噪聲。另一方面，越來越多的研究者從計算流體動力學(xué)（CFD）和聲比擬理論出發(fā)，將旋轉(zhuǎn)葉片噪聲源視為一系列旋轉(zhuǎn)運動點源分布，采用時域[4-5]或頻域[6-7]數(shù)值方法預(yù)測噪聲。

研究旋轉(zhuǎn)點源的聲場精確解有助于發(fā)展高精度的噪聲數(shù)值預(yù)測方法。吳九匯等[8]最早應(yīng)用球諧級數(shù)展開方法建立了旋轉(zhuǎn)點聲源在自由空間中的理論解析解。Poletti[9-10]提出了旋轉(zhuǎn)二維和三維單極子點源聲輻射精確時域球諧解。Carley[11]提出了一種用級數(shù)展開法預(yù)測旋轉(zhuǎn)單極子點源輻射噪聲的理論模型。Mao等[12]利用球坐標(biāo)系和柱坐標(biāo)系級數(shù)展開法，推導(dǎo)了旋轉(zhuǎn)單極子和偶極子輻射聲場的精確解。進(jìn)一步地，Mao 等[13]推導(dǎo)了旋轉(zhuǎn)四極子點源聲輻射精確解。這些理論模型都假設(shè)點源在自由場中旋轉(zhuǎn)。

實際應(yīng)用中，旋轉(zhuǎn)點源輻射的聲波會因周圍聲學(xué)邊界條件的存在而發(fā)生散射效應(yīng)，比如機(jī)翼和機(jī)身等結(jié)構(gòu)對螺旋槳噪聲的散射，使得實際聲場與自由空間輻射聲場存在明顯差異。在已發(fā)表的學(xué)術(shù)文獻(xiàn)中，使用理論解析方法對結(jié)構(gòu)的聲散射進(jìn)行了有限的關(guān)注。Glegg[14]將螺旋槳聲源建模為旋轉(zhuǎn)偶極子點聲源、中心體簡化為無限長的剛性圓柱體，得到一個預(yù)測螺旋槳遠(yuǎn)場聲壓的公式，其中包含了輪轂的散射。Kingan 等[15]擴(kuò)展了Glegg 的工作以評估中心體散射對開式轉(zhuǎn)子噪聲的影響。Hanson等[16]提出了一種計算飛行器機(jī)身散射開式轉(zhuǎn)子噪聲的方法，機(jī)身被簡化為無限長聲學(xué)剛性圓柱體。McAlpine等[17]應(yīng)用Hanson 等的分析進(jìn)一步評估機(jī)身散射對旋轉(zhuǎn)單極子或偶極子點源以及開式轉(zhuǎn)子噪聲的影響。Kingan等[18]也進(jìn)行了類似的研究。毛義軍等[19]對旋轉(zhuǎn)聲源輻射噪聲和邊界散射噪聲的預(yù)測方法進(jìn)行了綜述。

理論模型對實際物理問題進(jìn)行了太多簡化，需要發(fā)展可靠的數(shù)值方法以提高旋轉(zhuǎn)點源的聲場預(yù)測精度。本文基于點源模型，利用邊界積分方法研究二維聲學(xué)邊界對旋轉(zhuǎn)點源的聲散射效應(yīng)。為驗證數(shù)值方法的可靠性，對二維圓柱的聲散射進(jìn)行了理論建模分析。

1 聲散射計算的邊界積分方法

觀察點與源點的位置和時間分別用x和y以及t和τ表示。如圖1所示，剛性固體邊界S附近有一簡諧振動角頻率為ωs的單極子點源q（y,τ），以角頻率ωr繞o點逆時針旋轉(zhuǎn)，其初始相位和振幅分別為0和1，即點源可表示為：

圖1 旋轉(zhuǎn)點源聲傳播示意圖

那么，時域下的聲學(xué)控制方程為：

式中：c0為聲傳播速度。在剛性固體邊界上，聲場應(yīng)滿足聲學(xué)硬邊界條件：

方程式（2）的基本解為自由空間格林函數(shù)G0（x,y,t-τ ），滿足方程（4）：

利用方程式（2）～方程式（4）可得到：

根據(jù)Fourier變換定義：

頻域自由空間格林函數(shù)為：

因此方程式（5）可在頻域下寫成：

式中右端兩項分別為入射聲和散射聲。

以點源旋轉(zhuǎn)中心為坐標(biāo)原點建立極坐標(biāo)系，觀察點x和源點y的位置分別為（rx,θx）和（ry,θy），假設(shè)點源的初始位置為θ0，則有：

假設(shè)聲傳播介質(zhì)靜止、均勻，對二維聲傳播問題，頻域自由空間格林函數(shù)為：

其級數(shù)展開形式為：

其中：Jm（ · ）表示m階Bessel 函數(shù)，Hm（2）（ · ）表示m階第二類Hankel函數(shù)，r＜=min（rx,ry），r＞=max（rx,ry），k=ω c0為聲學(xué)波數(shù)。將式（1）、式（10）和式（12）代入方程式（9）得到：

利用Dirac函數(shù)性質(zhì)：

得到聲入射的表達(dá)式為：

將散射邊界離散為N個網(wǎng)格單元，采用等參邊界元方法求解方程（8）以獲得散射邊界上的聲場。當(dāng)觀察點置于散射邊界時，積分方程式（8）可離散為一系列線性方程組：

寫成矩陣形式為：

其中：

2 圓柱聲散射理論解析解

二維圓柱幾何形狀簡單，易獲取散射場的理論解析表達(dá)式。將觀察點的聲學(xué)壓力場分解為自由入射聲場和圓柱散射聲場兩部分，前者可由式（15）決定，還需得到后者的解析表達(dá)。類似于式（9），假設(shè)圓柱聲散射的表達(dá)為：

2.1 點源與圓柱同軸旋轉(zhuǎn)

點源q（y,τ）繞圓柱同軸旋轉(zhuǎn)如圖2 所示，假設(shè)散射格林函數(shù)可表示為：

圖2 點源繞圓柱同軸旋轉(zhuǎn)示意圖

當(dāng)觀察點位于圓柱表面時，將式（12）和式（20）代入聲學(xué)邊界條件（3），可解得：

式中：

即有：

將式（1）、式（10）和式（22）代入式（19），利用式（14）得到：

2.2 點源與圓柱異軸旋轉(zhuǎn)

如圖3 所示，二維圓柱附近有一旋轉(zhuǎn)單極子點源q（y,τ），其旋轉(zhuǎn)中心與圓柱中心的距離為d，以點源旋轉(zhuǎn)中心點為坐標(biāo)原點建立極坐標(biāo)系，對應(yīng)的觀察點x和源點y的位置分別為（rx,θx）和（ry,θy）；同時以圓柱中心點為原點建立另一極坐標(biāo)系，對應(yīng)觀察點x位置為當(dāng)觀察點位于圓柱表面時，顯然有rx＞ry和＜d。根據(jù)Graf加法定理，

圖3 點源與圓柱異軸旋轉(zhuǎn)示意圖

可將自由空間格林函數(shù)改寫為：

其中：?=π+α。

進(jìn)一步令：

利用式（1）和式（10），分別將式（25）代入式（9）、式（26）代入式（19），根據(jù)邊界條件式（3）和關(guān)系式（14），可解得：

式中：

從而散射格林函數(shù)可寫為：

將式（28）代入式（19），積分可得：

旋轉(zhuǎn)單極子源的解對m具有無限項，在實際應(yīng)用過程中應(yīng)以截斷方式采用有限項計算，截斷項數(shù)M與觀察點感興趣的聲波頻率ω有關(guān)，當(dāng)ω小于源諧振頻率ωs時,截斷項數(shù)應(yīng)保證ω＜Mωr-ωs，也就是；當(dāng)ω大于源諧振頻率ωs時，應(yīng)滿足。除此之外，根據(jù)Bessel 函數(shù)的指數(shù)衰減性質(zhì)，M不小于，[ · ]表示向上取整數(shù)?；谏鲜龅姆治觯財鄶?shù)M應(yīng)該滿足：

3 二維圓柱算例驗證

采用等參邊界元方法計算單極子點源與圓柱同軸和異軸旋轉(zhuǎn)的聲輻射特性，將結(jié)果與理論解析解對比，驗證數(shù)值方法的正確性。取圓柱半徑a=1.0 m，單極子點源的初始位置θ0=0°，其諧振頻率和旋轉(zhuǎn)頻率分別為

3.1 點源與圓柱同軸旋轉(zhuǎn)算例驗證

如圖2 所示，取ry=1.5 m。采用邊界元方法進(jìn)行聲傳播計算時，需要足夠的網(wǎng)格數(shù)量以保證數(shù)值結(jié)果與網(wǎng)格無關(guān)。在圓柱表面均勻分布180個網(wǎng)格單元，對500 Hz 以下的聲波，一個波長范圍內(nèi)不少于20個網(wǎng)格單元，可以滿足網(wǎng)格無關(guān)性要求。數(shù)值積分采用四點高斯積分方法，線性方程組的求解采用全選主元高斯消去法。

圖4為觀察點（rx=10 m,θx=0°）的聲壓幅值譜，邊界元數(shù)值解與理論解吻合一致。與固定點源產(chǎn)生的單頻聲場不同，旋轉(zhuǎn)點源因存在多普勒效應(yīng)，產(chǎn)生的噪聲分布在一系列離散頻率上，離散頻率滿足關(guān)系式ω=ωs-mωr，其中m為諧波階次；最高峰值頻率為f=160 Hz，對應(yīng)的諧波階次為m=4。在低頻（小于100 Hz）和高頻（大于400 Hz）處，聲壓幾乎為0。

圖4 點源繞圓柱同軸旋轉(zhuǎn)時觀察點聲壓幅值譜（rx =10 m,θx =0°）

圖5為圓柱表面上160 Hz聲壓幅值指向性分布的邊界元數(shù)值結(jié)果與理論解析結(jié)果的對比圖。對不同的網(wǎng)格單元，數(shù)值解與理論解的一致性非常好。對相同極半徑的觀察點，聲壓的幅值是相同的，呈現(xiàn)出典型的單極子聲場特性。

圖5 點源繞圓柱同軸旋轉(zhuǎn)圓柱表面160 Hz聲場指向性分布

圖6 為邊界元方法計算得到的140 Hz、200 Hz和260 Hz總聲壓實部的空間分布云圖，對應(yīng)的諧波階次分別為m=6、m=0 和m=-6。從圖中可以看到，m=0的聲波呈同心圓狀向外輻射；而m≠0的聲波則以螺旋狀向外輻射，螺旋瓣的數(shù)量與 |m|相同，且m＞0 時螺旋瓣的方向與點源旋轉(zhuǎn)方向相反，而m＜0 時則相同。

圖6 點源繞圓柱同軸旋轉(zhuǎn)時聲場分布云圖

3.2 點源與圓柱異軸旋轉(zhuǎn)算例驗證

對如圖3 所示的點源與圓柱異軸旋轉(zhuǎn)算例，選取ry=0.5 m、0°和d=2.0 m。邊界元方法的數(shù)值計算采用與同軸旋轉(zhuǎn)算例一致的網(wǎng)格，以及相同的數(shù)值積分和線性方程組求解方法。

圖7 為觀察點（=10 m，=0°）的聲壓幅值譜，邊界元數(shù)值解與理論解再次吻合一致。異軸旋轉(zhuǎn)點源產(chǎn)生的噪聲分布在一系列離散頻率上，這與同軸旋轉(zhuǎn)的頻譜特性相似，最高峰值頻率為f=210 Hz，對應(yīng)的諧波階次為m=-1。

圖7 點源繞圓柱異軸旋轉(zhuǎn)時觀察點聲壓幅值譜（=10 m, =0°）

在=10 m 的圓周上取不同的觀察點，將數(shù)值計算的聲壓幅值與理論解析結(jié)果進(jìn)行對比。圖8為210 Hz聲場的聲壓幅值指向性分布，對不同的觀察點，數(shù)值解與理論解是一致的。由于圓柱對聲波的散射作用，聲場的空間分布趨于復(fù)雜化，不像同軸旋轉(zhuǎn)那樣具有顯著的規(guī)律特征。

圖8 點源繞圓柱同軸旋轉(zhuǎn)時聲場指向性圖（ =10 m, f =210 Hz）

圖9為邊界元方法計算得到的210 Hz聲場的自由入射聲壓、圓柱散射聲壓和總聲壓實部的空間分布云圖。入射聲波以明顯的螺旋狀向外輻射，這是由點源的旋轉(zhuǎn)運動特性決定的。與同軸旋轉(zhuǎn)算例不同的是，圓柱對異軸旋轉(zhuǎn)點源入射聲波的散射不再具有單極子的螺旋狀特性，而呈現(xiàn)顯著的偶極子特征，使得總聲場在多個方向存在強(qiáng)傳播。

圖9 點源繞圓柱異軸旋轉(zhuǎn)時210 Hz聲場分布云圖

4 結(jié)語

對任意形狀的剛性散射邊界，利用自由空間格林函數(shù)的級數(shù)展開式求解線性聲學(xué)波動方程，建立了旋轉(zhuǎn)單極子點源聲輻射的頻域邊界元數(shù)值計算方法。同時，根據(jù)貝塞爾函數(shù)的加法定理，推導(dǎo)了點源與二維剛性圓柱同軸、異軸旋轉(zhuǎn)狀態(tài)下的聲場理論解析表達(dá)式。

對點源與圓柱同軸、異軸旋轉(zhuǎn)兩個算例，邊界元方法的數(shù)值解均與理論解析解吻合一致。點源向外輻射聲波的頻率由點源諧振頻率、旋轉(zhuǎn)頻率和諧波階次共同決定。對同軸旋轉(zhuǎn)算例，諧波階次為0 的聲波以同心圓狀向外輻射，不為0 的聲波則以螺旋狀向外輻射，螺旋瓣的數(shù)量與諧波階次的絕對值相同。對異軸旋轉(zhuǎn)算例，圓柱的聲散射呈現(xiàn)顯著的偶極子特征，聲場的指向性復(fù)雜，不像同軸旋轉(zhuǎn)算例具有顯著的螺旋輻射特征。

在后續(xù)的研究中，將進(jìn)一步考慮旋轉(zhuǎn)偶極子、四極子點源的聲輻射特性和運動介質(zhì)的影響，以及將二維聲傳播的研究結(jié)果推廣到三維聲學(xué)問題的研究。