二維圓柱對(duì)旋轉(zhuǎn)單極子點(diǎn)源的聲散射

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

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

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

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

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

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

1 聲散射計(jì)算的邊界積分方法

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

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

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

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

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

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

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

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

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

式中右端兩項(xiàng)分別為入射聲和散射聲。

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

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

其級(jí)數(shù)展開(kāi)形式為：

其中：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個(gè)網(wǎng)格單元，采用等參邊界元方法求解方程（8）以獲得散射邊界上的聲場(chǎng)。當(dāng)觀察點(diǎn)置于散射邊界時(shí)，積分方程式（8）可離散為一系列線性方程組：

寫成矩陣形式為：

其中：

2 圓柱聲散射理論解析解

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

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

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

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

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

式中：

即有：

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

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

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

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

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

其中：?=π+α。

進(jìn)一步令：

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

式中：

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

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

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

3 二維圓柱算例驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)語(yǔ)

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

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

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