水下對轉槳非空化線譜噪聲分析與數值研究

曾賽，杜選民，范威

(上海船舶電子設備研究所 水聲對抗技術國防科技重點實驗室，上海201108)

0 引言

在航空領域廣泛應用的對轉螺旋槳具有推進效率高、扭矩平衡等優點，也被應用于水下高速航行體、高速快艇等［1］。近年來，水下或者水上航行體的聲隱身指標越來越重要，作為重要噪聲源的螺旋槳始終都是被關注的重點。對轉槳由于前后槳反向旋轉，使得流場的相互作用干擾非常明顯，而且梢渦分離特征復雜，對其進行數值噪聲預報需要耗費大量計算資源，加之其軍事意義顯著，使得有關水下對轉槳噪聲預報分析的文獻非常少，大量對轉槳噪聲預報分析的研究集中在航空領域。

最突出的工作來自于Hanson［2］和Parry［3］，前者基于螺旋面理論和升力面理論，在頻域的范圍內建立了航空對轉槳的噪聲預報模型，并且實驗驗證了模型的準確性;Parry［3］則建立了適于工程計算的航空對轉槳噪聲計算公式。隨著計算機硬件技術的提高，計算氣動聲學(CAA)與計算流體力學(CFD)的耦合計算方法成為為計算航空螺旋槳噪聲的首選，Peters［4］正是基于該方法分析了氣動對轉槳的噪聲譜并闡述了線譜噪聲的產生機理。借助航空領域中的成功經驗，國內學者針對水下對轉槳的噪聲預報有過一些嘗試。20 世紀90年代，朱錫清等［5］首先討論了對轉槳線譜噪聲的理論計算方法，其主要基于升力面理論和聲學方法來實現。王順杰等［6］則對水下對轉槳空化線譜頻率進行了預報和數值模擬，其基本出發點是認為空泡呈現單極子輻射特征，利用聲類比方程推導了遠場的聲壓表達式。楊瓊方等［7］則研究了伴流場中對轉槳空化初生的判斷準則，提出了判斷空化初生的4 個充分條件，并且用大規模的數值計算預報了對轉槳空化初生時的線譜和寬帶譜噪聲。

不過，隨著螺旋槳設計水平的提高，對轉槳設計的趨勢是多槳葉、大盤面比、大側斜等，其最終目的是避免空化和降噪。因此水下對轉槳大部分工況下是非空化的。本文主要研究非空化線譜輻射噪聲特性，進行數值預報，并對結果進行對比分析。

1 對轉槳非空化線譜噪聲的理論與分析

1.1 廣義聲類比方程的變形

由廣義聲類比方程［8］可知，螺旋槳噪聲場可以用單極子、偶極子和四極子來描述。

式中:G 表示格林函數;vn、fi均取運動物體表面的外法線方向為正;Tij為ligthhill 應力張量;fi表示流體作用在槳葉表面單位面積上的力的大小;τ= -T 表示在時間平穩過程下，遙遠過去對現在的影響;τ =T 表示外伸波條件。右邊第1 項表示四極子作用，其本質是lighthill 應力源，在低馬赫數時輻射效率很低可以忽略;右邊第2 項為偶極子性質的力源;右邊第3 項為單極子聲源，只有當槳葉出現空化時，才是主要的發聲源。本文研究非空化對轉槳的非空化線譜噪聲問題，偶極子的力源是討論的重點。

槳葉的實際運動軌跡如圖1所示，為了將槳葉表面的力源與遠場聲壓聯系起來，建立4 套坐標系，即大地坐標系(x，r，φ)、隨動坐標系(x1，r，φ)、葉片固定坐標系(γ，r，ξ)和聲源螺旋坐標系(γ0，r0，，ξ0)，圖2為聲源螺旋坐標系示意圖，表示螺旋槳徑向距離為r0時葉切面的坐標系，γ0方向與聲源點運動方向相反，表示在螺旋方向上離開t=0 時刻螺旋槳母線位置的距離。其中，MCA 為槳葉側斜距離，FA 為槳葉縱傾距離。

圖1 槳葉螺旋面軌跡示意圖Fig.1 Schematic diagram of blade＇s helical track

通過坐標系的變換，可將力源點坐標(γo，ro，ξo)用直角坐標(y1，y2，y3)表示，觀測點所在的柱坐標(x，r，φ)所對應的直角坐標為(x，rcos φ，rsin φ)，可以計算得到觀測點與力源點之間的距離為

圖2 槳葉切面螺旋坐標系示意圖Fig.2 Helical coordinate system of blade＇s tangent plane

式中:Ω 為槳葉的旋轉角速度。

物體運動的法向速度可以表示為:vn= u0·為中間面在運動方向的斜率。由曲面積分可得dγ0dr0，式中:J 為雅可比算子。另外，自由空間中的格林函數為式中:R 表示源點到觀察者的距離;c0表示水中聲速。由δ 函數的性質，可以將h'表示為對于槳葉上作用力亦可以寫成:fj(y，τ)=等式中的ω0表示槳葉上脈動力的圓頻率，下標0 表示作用力為定常力。另外，由于聲源以u0速度沿著-γ0方向運動，實際的源坐標應該用γ0+u0τ 代替，綜合這些考慮，可將(1)式變為

式中:g 為源函數。引入傅里葉變換［9］，

將(4)式代入(3)式對τ 積分可以得到:

若考慮遠場情形，將(2)式作泰勒展開，可以得到:

故可以將(5)式中1/R 用1/r 代替，如圖3所示。

圖3 觀察點示意圖Fig.3 Schematic diagram of viewport in field

1.2 干涉作用以及周向諧波流場作用產生的遠場聲壓

對轉槳非空化線譜噪聲主要由以下7 種作用產生:非均勻進流與前槳導邊相互作用發生、槳葉面附近大尺度渦對槳葉非定常誘導力作用、前槳葉面二次非定常流對后槳作用、后槳抽吸作用形成的非定常壓力場與前槳相互作用、前槳梢渦脫落與后槳導邊相互作用、前槳隨邊尾渦脫落與后槳導邊相互作用、前槳槳轂尾渦及邊界層流與后槳相互作用。這7 種作用又可以分為兩類:前后槳相互作用和流場對槳的作用，如圖4所示。

干涉作用可以細分為3 個部分:粘性流場、隨邊勢流場和邊界層效應。可以用非定常力來衡量這3 個部分作用。設前槳葉片數為B1，旋轉角速度為Ω1，后槳葉片數為B2，旋轉角速度為Ω2，前槳葉在通過后槳靜止葉片時會產生一個干涉流場，表現形式為干涉脈沖，這樣在后槳上會產生一個前槳葉片通過頻率的單頻聲源，這樣可以將前槳對后槳作用產生的干涉擾動頻率表示為

圖4 前后槳葉干涉作用(左)與周向諧波流場作用(右)示意圖Fig.4 Schematic diagram of interference effect and circumferential harmonic field effect

式中:kg1表示前槳對后槳干擾頻率的諧波階數。

取時間因子eiωt，后槳第1 個槳葉上的作用力可以表示為

那么第N 個后槳葉上的作用力應該等于第1 個槳葉在位置X、時間T 時的作用力，其中:

將(9)式代入到(8)式既可得到第N 槳葉的作用力。另一方面，由時間引起的相位變化為

由坐標位置引起的相位變化為

由(10)和式(11)式可知，后槳所有槳葉的相位變化疊加滿足

再由特殊函數［10］變形，

將(6)式～(14)式變形代入到(5)式中，即可以得到后槳在干涉作用下所產生的遠場聲壓:

考慮流場對槳葉輻射聲壓的影響。后槳工作在前槳的非均勻尾流之中，由于后槳的切水作用，會帶動周圍流場一起旋轉使其具有明顯的周向調諧性質，設對后槳的調諧分量為W2，對前槳的調諧分量為W1，以后槳為研究對象，每個調諧分量會對后槳葉產生一個激勵，會有對應的通過頻率，可以將該周向諧波流場的作用頻率記為:ω0=k2wW2Ω2(k2w、W2=0，1，2，…)，k2w為后槳調諧諧波分量的階數。與干涉作用的推導類似，經過一系列的特殊變形可以得到后槳由于流場非均勻流場而產生的遠場聲壓輻射為(15)式和(16)式中積分下、上限分別表示槳葉葉根和葉梢。

將二者相加既可以得到后槳的遠場聲壓輻射表達式，同理可以得到前槳的遠場聲壓輻射表達式(分析過程不再贅述)。再將前、后槳遠場聲壓輻射公式相疊加即可以得到總的遠場聲壓輻射表達式。上面等式中，干涉作用參數:

周向調諧流場對槳的作用參數:

式中:l=1，2，分別表示前、后槳;B1、B2分別表示前、后槳的槳葉數;Dl表示槳的直徑，前后槳取各自的直徑;rl表示微元源點與測點之間的距離;bl表示前后槳槳葉片對應的弦長;MCA、FA 分別為槳葉側斜和縱傾的距離;r 為源點到觀測點之間的距離;φ(l)表示計算的初始相位角;帶下標w 的量表示是流場非均勻性非定常性的影響;Max、Mar和Matip分別表示前進馬赫數、微元槳葉切片馬赫數和葉片葉梢馬赫數;q、m、k1g和k2w表示諧波階數;θ 表示軸向平面上，測點與軸向的角度，由于軸對稱性，這個角度和表示初始角的角度可以將該對轉槳整個輻射聲場的特性表示出來;CDk、CLk、CDw和CLw分別表示在兩種作用下的升阻力系數;ψ 表示形狀函數［2］。

本文利用螺旋面理論結合廣義聲類比方程和特殊函數理論得到對轉槳在兩種作用下各自的非空化線譜輻射噪聲表達式，與文獻［5］中利用升力面理論和聲學方法得到的表達式有4 個方面的不同。其一，本文同時考慮了升阻力作用，這顯然更符合真實情形，而文獻［5］只考慮了非定常升力作用。其二，本文變形得到的貝塞爾函數宗量與文獻［5］不同，本文中，干涉作用下貝塞爾函數的宗量中所體現的干擾頻率是加性的，加性頻率為kg1B1Matip1+mB2Matip2，這與(7)式的物理意義是相符合的，而文獻［5］中干涉作用下貝塞爾函數的宗量所體現的干擾頻率是減性的，有差頻的物理意義，文獻［5］中沒有給出物理解釋。其三，無量綱的波數不同，無量綱波數Ky反應了在干涉作用下前后槳輻射噪聲重要的相位信息，其對輻射指向性也有重要影響，在氣動對轉槳噪聲領域，Bradley［11］證明了貝塞爾函數階數的正負(mB2-kg1B1)對于決定聲場的源分布有重要作用，文獻［5］并沒有強調該無量綱波數的作用。其四，本文選用時間因子為eiωt，和文獻［5］不同，因而得到的最終輻射表達式也會與文獻［5］不同。

1.3 非空化線譜輻射噪聲的頻率分析

分析對轉槳輻射噪聲表達式，考慮等式中波的傳播因子可以給出不同模式成分波的頻率，不失一般性，考慮(15)式中后槳在前槳干涉作用下的傳播因子:

該式表征對轉槳的不同模式空間旋轉輻射波在φ方向、特征值為(mB2- kg1B1)時，其旋轉角速度為(mB2Ω2+kg1B1Ω1). 等式中的時間項是與特征非空化線譜對應的頻率。可以看出4 類非空化線譜噪聲頻率，分別為軸頻APF、前槳葉頻BPF1及其諧波、后槳葉頻BPF2及其諧波、干涉作用頻率pBPF1+qBPF2.線譜頻率可以看成這4 類非空化線譜頻率的組合，即f=sAPF+pBPF1+hBPF2.

1.4 輻射噪聲聲壓方向性討論

通過輻射聲壓表達式可以看出，對轉槳輻射聲壓幅值以及方向性均是由貝塞爾函數決定的。討論(15)式中后槳在干涉作用影響下輻射聲壓的方向性。定義聲壓方向函數

式中:N=mB2-k1gB1，表示貝塞爾函數的階數，這里定義為輻射模態，(m，k1g)表示非空化線譜的噪聲諧波數。將0 階和前6 階輻射模態的無量綱聲壓方向性曲線給出，如圖5所示。

圖5 不同模態聲壓方向性示意圖Fig.5 Schematic diagram of sound pressure directivity at different modes

從圖5可以看出，除0 階輻射模態無方向性，前6 階輻射模態聲壓方向均表現為“8”字形，呈現偶極子特性，隨著階數的增加，“8”字形越來越尖銳。

1.5 非空化線譜的強弱分析

任意諧波處的聲壓幅值大小均是由貝塞爾函數決定的，根據貝塞爾函數的性質，對于第1 類貝塞爾函數Jn(x)，其階數與宗量的差別越大，貝塞爾函數的幅值越小，即當n?x 時，貝塞爾函數值趨于0. 根據貝塞爾函數的漸進公式:

可見，貝塞爾函數在理論上是有無窮多個零點。當干涉作用的貝塞爾函數和周向諧波流場作用的貝塞爾函數同時取0 時，一定沒有強線譜出現，可能出現很弱的線譜以至于混淆在連續譜之中。當二者同時取峰值時，一定有強線譜出現。其他的情形都介于二者之間，比如在強線譜附近會出現幅值較弱的一些線譜，但這些線譜又比連續譜強的多，可以認為是較強線譜。這可以作為對轉槳非空化強弱線譜出現的條件。

2 對轉槳數值模擬以及分析

2.1 數值模型以及網格劃分

分別采用Reliable k-ε 湍流模型和RNG k-ε 湍流模型［12］來計算模型槳的水動力性能，并與實驗值［13］進行比較，從而得出模擬流場更精確的湍流模型，然后與滑移網格模型相結合，利用FW-H 方程計算目標槳在流場中的聲場。本文所用的對轉槳基本幾何參數如表1所示，模型如圖6所示。

表1 對轉槳的幾何參數Tab.1 The geometry parameters of CR propeller

圖6 對轉槳示意圖Fig.6 Schematic diagram of counter-rotation propeller

將對轉槳置于與其同軸的圓柱體流場之中，上游進流面取為6D，下流面取為10D，徑向半徑取為5D，D 為前槳直徑。螺旋槳在旋轉時會帶動周圍的流體一起旋轉，因此在對轉槳周圍取一個與其同軸的圓柱體作為旋轉區域，徑向半徑取1.3D. 網格劃分采用結構-非結構網格。非結構網格1 471 892 個，結構網格2 064 998 個。

2.2 邊界條件的設置

整個計算域可以分為入口、出口、靜止區域、旋轉區域、靜止圓柱面、旋轉圓柱面和螺旋槳。邊界條件分別設置為速度入口、壓力出口、靜止區域、前槳旋轉區域、后槳旋轉區域、動靜區域交界面、靜止圓柱面無滑移墻邊界、對轉槳葉面設置為運動墻邊界和旋轉邊界，動靜域用滑移網格連接。流場的變化通過改變對轉槳的進速系數來確定。

2.3 敞水性能驗證

進速系數分別選取0.3、0.5、0.7、0.9、1.1，對轉槳前后槳轉速設為1 072 r/min，轉向相反。對轉槳的進速系數、推力系數kt以及扭矩系數kq按照文獻［14］處理，得到對轉槳的敞水性能曲線與試驗結果［13］的比較，如圖7所示。利用RNG k-ε 模型計算的推力系數與試驗值相對誤差為3.4%、3.7%、2.2%、1.0%、3.1%，扭矩系數與試驗值相對誤差為0.07%、1.5%、0.3%、3.9%、7.5%. 利用Realizable k-ε 模型計算的推力系數與試驗值相對誤差為9.1%、6.3%、1.1%、2.8%、4.6%，扭矩系數與試驗值相對誤差為0.9%、2.0%、3.0%、4.7%、8.3%.可見RNG k-ε 模型要優于Realizable k-ε，下面用RNG k-ε 模型結合FW-H 方程數值計算對轉槳在均勻進流條件下非空化線譜噪聲。

圖7 對轉槳敞水性能曲線Fig.7 The open water performance curves of CR propeller

2.4 對轉槳數值非空化離散線譜模擬及分析

按照RNG k-ε 模型、滑移網格模型與FW-H 方程相結合的方法來預報目標對轉槳的非空化離散噪聲。首先計算穩態條件下對轉槳所在流場，計算收斂后保存計算結果并將其作為非定常流動初始值計算非穩態條件下的對轉槳流場［6］，當監視的無量綱系數穩定后進入聲學模塊的計算。如圖8中前、后槳升力系數保持動態穩定，圖9表示半個計算周期內前、后槳推力系數ktf、kta均保持動態穩定，表明能夠進行聲學計算。由于螺旋槳噪聲頻譜通常在0 ～10 kHz，故設置計算時間步長為5 ×10-5s，計算時長大于0.5 s，所得頻率分辨率小于2 Hz，滿足分析要求。待計算存儲足夠多時間步長上測點的脈動聲壓值后，進入聲場分析。按照文獻［15］的結論，對轉槳線譜大多集中在前13 階葉頻范圍內，按照本文設計的對轉槳轉速估算，目標槳線譜會在0 ～2 kHz 的范圍內。

圖8 計算無量綱升力系數變化曲線Fig.8 The calculated lift coefficient curves of CR propeller

圖9 半個周期內前、后槳推力系數變化Fig.9 The thrust coefficient curves of front and rear propellers in half period

圖10 為測點布置示意圖，其中測點P1 ～P35距槳軸中心為6R1，R1為前槳半徑。測點沿著圓周等距分布且沿著軸向面布置，P37、P39 測點分別布置在前、后槳的徑向方向上，P38 布置在前、后槳中間平面上。分別觀察5 個典型測點的聲壓頻譜特性，如圖10 所示，f0、f1、f2分別為軸頻和前、后槳葉頻。在軸向測點P1 的聲壓頻譜上，線譜主要出現在葉頻上，峰值線譜出現在組合葉頻5f1+4f2上，此組合對應于貝塞爾函數的階數為0 階，這說明前、后槳之間的干涉作用在軸向表現得較為明顯;在徑向測點P38 的聲壓頻譜上，線譜幾乎都出現在葉頻及其諧波上，考慮到葉頻與軸頻的倍數關系，可以認為槳葉徑向處的線譜都出現在軸頻及其諧波處;在軸向和徑向之間的測點的聲壓頻譜上，如測點P7、P14，聲壓譜峰值幾乎都出現在前、后槳相互干擾的頻率上，這說明在非徑向和非軸向處，前、后槳之間的干涉作用非常顯著。如圖11 所示，在5f1+4f2頻率對應的功率譜密度最大，與聲壓頻譜相吻合。如圖12 所示，在1 階葉頻以及f1+ f2處對應的功率譜密度最大，隨著諧波階數增加、功率譜密度逐漸減小。如圖13 所示，在f2和f1+f2頻率對應的功率譜密度最大。結合聲壓譜可以看出，在0 ～2 000 Hz 的范圍內存在多條與前、后槳葉頻相關的線譜。如圖14 所示為各測點所得聲壓級連接曲線，可以看出與1.4 節所得的聲壓方向性大致吻合，總聲壓級呈“8”形狀。

圖10 典型測點的聲壓譜示意圖Fig.10 Schematic diagram of typical sound pressure spectra at measuring point

圖11 P1 點的功率譜Fig.11 Power spectrum at P1

圖12 P38 點的功率譜Fig.12 Power spectrum at P38

圖13 P39 點的功率譜Fig.13 Power spectrum at P39

圖14 測點聲壓級變化曲線Fig.14 SPL curve at measuring point

3 結論

通過分析對轉槳非空化噪聲的特性，以期為水下目標識別和分類提供判斷依據。通過理論和數值計算分析，得到如下結論:

1)水下對轉槳非空化線譜噪聲由7 種作用共同作用形成。從機制上來看可以分為兩類，即前、后槳干涉作用和非均勻流場作用。這些因素共同作用使得對轉槳噪聲線譜非常豐富。

2)水下對轉槳非空化線譜噪聲預報頻率為f=sAPF+pBPF1+hBPF2，系數s、p、h 為整數。

3)水下對轉槳非空化噪聲輻射聲壓具有方向性。徑向方向的總聲壓要強于軸向方向。

4)計算對轉槳水動力性能方面，RANS 方法能夠取得符合工程應用的精度，且RNG k-ε 模型要優于Reliable k-ε 模型。

5)通過uRANS 方法結合FW-H 方程的混合數值算法能夠合理地模擬對轉槳非空化噪聲問題，模擬的非空化線譜與理論吻合度較好。

References)

［1］常煜，洪方文，張志榮，等.對轉槳水動力性能的數值分析［C］∥2008年船舶水動力學學術會議暨中國船舶學術界進入ITTC30 周年紀念會論文集.杭州:中國造船工程學會，2008.CHANG Yu，HONG Fang-wen，ZHANG Zhi-rong，et al. Numerical study of contra-rotating propellers hydrodynamic performance［C］∥The Ship Hydrodynamics Conference and the Shipbuilding Academic of China Step into the ITTC30 Anniversary Conference.Hangzhou:China Shipbuliding Engineering Society，2008. (in Chinese)

［2］Hanson D B. Noise of counter-rotation propellers［J］. Journal of Aircraft，1985，22(7):609 -617.

［3］Parry A B. Theoretical prediction of counter-rotation propeller noise［D］. UK:University of Leeds，1988.

［4］Peters A. Assessment of propfan propulsion systems for reduced environmental impact［D］. US:Massachusetts Institute of Technology，2010.

［5］朱錫清，吳武生. 水下高速航行體對轉螺旋槳線譜噪聲預報研究［J］.聲學學報，1998，23(2):123 -134.ZHU Xi-qing，WU Wu-sheng. Prediction of line-spectrum noise induced by high speed vehicle counter-rotation propellers in water［J］. Acta Acustica，1998，23(2):123 -134.(in Chinese)

［6］王順杰，程玉勝，高鑫.水下對轉螺旋槳空化線譜頻率預報與數值模擬［J］.兵工學報，2013，34(3):311 -319.WANG Shun-jie，CHENG Yu-sheng，GAO Xin. Prediction and numerical simulation of cavitation noise line-spectrum frequency induced by underwater counter-rotation propeller［J］. Acta Armamentarii，2013，34(3):311 -319.(in Chinese)

［7］楊瓊方，王永生，張志宏，等.伴流場中對轉槳空化初生的判斷與輻射噪聲預報和校驗［J］. 聲學學報，2014，39(5):590 -606.YANG Qiong-fang，WANG Yong-sheng，ZHANG Zhi-hong，et al.Numerical prediction of cavitation inception radiated noise of contra-rotating propeller with non-uniform inflow［J］. Acta Acustica，2014，39(5):590 -606.(in Chinese)

［8］孫小峰，周盛.氣動聲學［M］.北京:國防工業出版社，1993.SUN Xiao-feng，ZHOU Sheng. Aeroacoustic［M］. Beijing:National Defense Industry Press，1993. (in Chinese)

［9］喬渭陽.航空發動機氣動聲學［M］. 北京:北京航空航天大學出版社，2010.QIAO Wei-yang. The aeroacoustic of aircraft engine［M］. Beijing:Beijing University of Aeronautics and Astronautics Press，2010. (in Chinese)

［10］梁昆淼.數學物理方法［M］.北京:高等教育出版社，2011.LIANG Kun-miao. The methods of mathematical physics［M］.Beijing:Higher Education Press，2011. (in Chinese)

［11］Bradley A J. A study of the rotor/rotor interaction tones from a contra-rotating propeller driven aircraft［R］. Reston，VA:AIAA，1986:1886 -1894.

［12］王福軍.計算流體動力學分析-CFD 軟件原理與應用［M］.北京:清華大學出版社，2010.WANG Fu-jun. Analysis of computational fluid dynamics-CFD software principle and application［M］. Beijing:Tsinghua University Press，2010. (in Chinese)

［13］Miller M L. Experimental determination of unsteady forces on counter-rotating propellers in uniform flow，SPD-659-01［R］.Maryland:David Naval Ship Research and Development Center，1976.

［14］張濤，陳彥勇，楊晨俊.RANS 方法預報對轉槳定常敞水性能研究［J］.船舶工程，2011，33(5):23 -26.ZHANG Tao，CHEN Yan-yong，YANG Chen-jun. Study on the prediction of contra-rotating propellers open-water steady performance using RANS formula［J］. Ship Engineering，2011，33(5):23 -26.(in Chinese)

［15］馬徐琨.淺析水下高速航行體對轉螺旋槳輻射噪聲線譜建模［J］.聲學學報，2002，27(6):503 -506.MA Xu-kun. Preliminary modeling of line spectrum for radiated noise induced by high speed counter rotation of underwater vehicle propellers［J］. Acta Acustica，2002，27(6):503 -506.(in Chinese)