王超，張立新，鄭小龍，魏勝任

(1.哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江 哈爾濱150001;2.武漢第二船舶設(shè)計研究所，湖北 武漢430074)

螺旋槳噪聲是船舶三大噪聲之一，并在其中占有極其重要的位置，主要由旋轉(zhuǎn)噪聲、渦流噪聲和空泡噪聲等組成。螺旋槳的旋轉(zhuǎn)噪聲是一種周期性的噪聲，它與螺旋槳槳葉負荷和厚度有關(guān)。渦流噪聲不只是水流對螺旋槳的沖擊，還包括由于槳葉葉梢和螺旋槳轂渦流破裂而產(chǎn)生的噪聲。旋轉(zhuǎn)噪聲和渦流噪聲都與聲學(xué)理論中的單極子、偶極子和四極子聲源機理相聯(lián)系。對于均勻來流下的螺旋槳無空化噪聲，主要由單極子噪聲和偶極子噪聲所組成［1］。關(guān)于螺旋槳噪聲的理論預(yù)報大致分為兩種:一種是利用勢流理論程序預(yù)報，一種就是利用粘流軟件進行數(shù)值模擬。朱錫清等［2-4］結(jié)合非定常升力面理論和聲類比方法對船舶螺旋槳的低頻線譜噪聲和寬帶譜噪聲進行了研究;Seol［5-6］等采用面元法求解流場和FWH方法時域內(nèi)求解聲傳播的耦合方法，計算了非均勻進流條件下DTRC4119無側(cè)斜標(biāo)準(zhǔn)槳無空化和空化狀態(tài)下的總聲級中厚度噪聲和負載噪聲各自所占的比重以及聲指向性，對非均勻流螺旋槳噪聲做了全面的分析;楊瓊方［7］采用流場大渦模擬(large eddy simulation，LES)和聲場邊界元數(shù)值聲學(xué)的弱耦合方法在頻域內(nèi)進行了螺旋槳噪聲分析，并對不同頻率槳葉表面激勵對聲輻射的影響進行分析。這些計算結(jié)果大多沒有表現(xiàn)出葉頻處聲壓峰值等信息。對這一問題，楊瓊方在文獻［8］中做了很好的完善。本文采用FLUENT軟件進行了大渦模擬(LES)非定常運算，對比了敞水性能和槳葉表面壓力分布，并對螺旋槳周圍流場的脈動壓力幅值進行了分析，且利用CFD和聲學(xué)無限元方法相結(jié)合的耦合方法預(yù)報了螺旋槳噪聲特性。

1 流-聲耦合方法理論

1.1 大渦模擬SGS模型

SGS模型在LES方法中占有十分重要的地位，本文應(yīng)用Smagorinsky-Lilly模型來模擬亞格子應(yīng)力:［9-10］

式中:μt是亞格子尺度的湍動粘度，在文獻［11］中推薦用下式計算:

1.2 Lighthill聲學(xué)類比理論

Lighthill考慮的模型為:在無限大的均勻、靜態(tài)聲介質(zhì)中包含一個有限的湍流運動區(qū)域V，因此與流動有關(guān)的聲源都集中在該區(qū)域內(nèi)。在區(qū)域V外，遠離湍流區(qū)域的流體中密度的波動和聲波相似，因此整理連續(xù)方程和動量方程，并簡化得到遠離湍流區(qū)域流體中的勻質(zhì)聲學(xué)波動方程:

Lighthill聲類比理論是從流體力學(xué)基本方程納維-斯托克斯(Navier-Stokes)方程導(dǎo)出的:

式中:c0是等熵條件下的聲速值;ρ＇=ρ-ρ0為有聲擾動時的密度分量，ρ與ρ0分別是擾動與未擾動時的流體密度;Tij是 Lighthill應(yīng)力張量［12］:

式中，σij表示雷諾應(yīng)力張量的粘性部分。

對于低速等熵流動，粘性應(yīng)力張量σij對Lighthill應(yīng)力張量Tij的貢獻遠小于雷諾應(yīng)力項ρuiuj，可以忽略不計，同時可以得到

因此，可得到Lighthill應(yīng)力張量的近似式Tij=ρuiuj。

1.3 ACTRAN 聲學(xué)原理

ACTRAN/LA基于 Lighthill方法，并結(jié)合 Curle’s理論:

1)Curle’s方程的體積分作為有限元區(qū)域的體源;

2)Curle’s方程的面積分作為邊界條件;

3)自由場的格林函數(shù)作為其他的邊界條件。

對式(6)在Ω邊界上積分，乘以測試函數(shù)δρ，并應(yīng)用分布積分產(chǎn)生弱變分形式，在面積分上應(yīng)用應(yīng)力張量，可得到下式:

1.4 螺旋槳噪聲模擬基本流程

ACTRAN軟件處理流致噪聲問題時，CFD計算與聲學(xué)計算是耦合的，即首先進行CFD仿真，提取出湍流信息，然后再利用Lighthill聲類比方法分析聲場。對于聲學(xué)分析中，只要滿足每波長6網(wǎng)格的規(guī)則即可，使用積分法將流場信息加載到聲學(xué)網(wǎng)格上，因此不需要對聲源區(qū)的網(wǎng)格做特別的優(yōu)化。圖1為流-聲耦合方法的計算流程圖。

圖1 流－聲耦合計算流程圖Fig.1 Flowchart of fluid-acoustics interaction calculation

2 計算前處理

2.1 幾何模型建立和計算域的劃分

文章研究螺旋槳為DTRC4119槳，其具有比較詳細的流場數(shù)據(jù)，被ITTC選為考證數(shù)值方法預(yù)報精度的標(biāo)準(zhǔn)槳，數(shù)據(jù)取自文獻［13］，幾何參數(shù)如表1所示。

表1 螺旋槳幾何參數(shù)Table 1 Geometry coefficients of propeller

圖2 螺旋槳幾何模型和計算域Fig.2 Geometry model of propeller and Computational domain

在建模過程中使用的是直角坐標(biāo)系O-XYZ，X軸方向代表來流方向，它沿著螺旋槳的旋轉(zhuǎn)軸指向下游，Y軸沿螺旋槳某一葉片的母線，Z軸服從右手定則。為了方便對螺旋槳進行結(jié)構(gòu)化網(wǎng)格劃分，只需要建立一個槳葉的幾何體和附近小域即可。創(chuàng)建完成的螺旋槳三維模型如圖2(a)所示。

計算域采用與螺旋槳同軸的圓柱流域，并被劃分為兩個區(qū)域，內(nèi)域包含螺旋槳。上游速度入口位置為槳前2倍直徑，下游自由出口設(shè)定在槳后5倍直徑，外邊界半徑大小為3倍直徑，如圖2(b)所示。

2.2 網(wǎng)格劃分及參數(shù)設(shè)定

網(wǎng)格劃分是CFD模擬過程中比較耗時的環(huán)節(jié)，也是直接影響模擬精度和效率的因素之一。本文在劃分網(wǎng)格時采用全結(jié)構(gòu)化網(wǎng)格。根據(jù)螺旋槳周期性特點，首先劃分單個槳葉流域內(nèi)的網(wǎng)格，然后進行旋轉(zhuǎn)，得到全流域內(nèi)的網(wǎng)格。劃分單槳葉流域網(wǎng)格時，采用H型網(wǎng)格，在槳葉位置處單獨布置C型網(wǎng)格與槳葉形狀匹配，此外在槳葉位置處向葉面外一側(cè)開O網(wǎng)(如圖3所示)，此種方式可以達到在槳葉表面進行局部加密，同時還可控制網(wǎng)格數(shù)量，本文邊界層第一層網(wǎng)格尺寸取 0.000 3，Y+控制在 10～100。

對于流場的關(guān)鍵區(qū)域(如槳葉隨邊、導(dǎo)邊、葉根與槳轂連接處、葉梢等)進行加密，以便捕捉到重要的流場信息;而對于距離螺旋槳較遠的區(qū)域網(wǎng)格將其密度適當(dāng)降低，便于控制總網(wǎng)格數(shù)，這樣，在網(wǎng)格模型總節(jié)點數(shù)一定的情況下可以提高計算精度，還可以避免流場變化平緩區(qū)域的計算資源浪費。圖4(a)是隨邊葉根處O型網(wǎng)格局部放大圖。槳葉網(wǎng)格如圖4(b)所示，整體計算域網(wǎng)格約為200萬。

圖3 槳葉外O網(wǎng)Fig.3 O-block out of propeller

螺旋槳壁面定義為不可滑移壁面條件，外壁面不考慮粘性作用，采用滑移壁面。計算時采用單一旋轉(zhuǎn)參考坐標(biāo)系模型，將旋轉(zhuǎn)參考坐標(biāo)系固定于螺旋槳中心，并以600 r/min的角速度進行旋轉(zhuǎn)，外域采用絕對靜止坐標(biāo)系，兩個域之間利用INTERFACE邊界進行連接，流場信息通過插值進行傳遞。

圖4 螺旋槳網(wǎng)格拓撲結(jié)構(gòu)Fig.4 Mesh topology of propeller

3 大渦模擬計算結(jié)果分析

計算螺旋槳流場，首先采用k-ε湍流模型進行定常運算，待獲得穩(wěn)定流場后，改用LES湍流模型進行非定常運算。計算時采用有限體積法進行離散，擴散項采用中心差分格式，壓力速度耦合采用SIMPLEC算法，連續(xù)性曲線小于0.000 01時認為計算收斂，時間步長設(shè)置為 0.002 5。

3.1 敞水性能計算結(jié)果驗證

為便于計算結(jié)果表述，定義如下參數(shù):

式中:J表示進速系數(shù)，KT表示推力系數(shù)，KQ表示扭矩系數(shù)。

對進速系數(shù)在0.5～0.9進行大渦模擬，對于非定常運算，流場并不是恒定不變，存在微小波動，待計算穩(wěn)定后，取整周期讀取其平均值作為計算結(jié)果，對比其敞水性能。

圖5為計算結(jié)果與試驗值［14］的比較，從中可以看出螺旋槳水動力性能計算值與試驗結(jié)果吻合較好，KT曲線最大誤差為1.5%，KQ曲線最大誤差1.417%，這說明本方法計算螺旋槳水動力性能具有很好的準(zhǔn)確性。

圖5 螺旋槳敞水性能曲線Fig.5 Open water performance of propeller

3.2 槳葉表面壓力系數(shù)對比分析

在研究螺旋槳的誘導(dǎo)脈動壓力、空泡和噪聲等問題時，求解出流場中的槳葉壓力分布是十分重要的。對比螺旋槳0.3R、0.7R葉剖面處的壓力系數(shù)CP，以進一步驗證計算方案的可靠性，如圖6所示。CP=，其中 (P-P0)為相對壓力，相對進流速度

圖6DTRC4119螺旋槳表面壓力系數(shù)分布(J=0.833)Fig.6 Pressure coefficient distribution of propeller DTRC4119(J=0.833)

由圖6可以看出，通過大渦模擬計算出的結(jié)果與試驗值之間誤差較小，0.7R處壓力分布與試驗值吻合較好，0.3R處相對吻合稍差，這可能是因槳轂簡化為圓柱體而省略了轂帽部分，而0.3R較貼近槳轂，使得該處表面壓力分布存在些許偏差。導(dǎo)邊(x/c=0.0)與隨邊(x/c=1.0)處偏差比葉面中部大。這是由于導(dǎo)邊與隨邊處流動梯度變化大，壓力變化明顯，造成模擬計算結(jié)果偏差較大，文獻［16］也揭示了相同的情況。

3.3 脈動壓力計算分析

為進一步獲取流場信息，對槳葉周圍流場中特征點的脈動壓力進行監(jiān)控分析，特征點的分布情況如圖7所示。

圖7 特征點位置Fig.7 Positions of the feature points

圖8是經(jīng)過FFT變換后的各特征點的脈動壓力幅值，由于螺旋槳旋轉(zhuǎn)作用，能夠明顯看出槳葉的葉頻(blade passing frequency，BPF)及其各倍葉頻。螺旋槳轉(zhuǎn)速為10 r/s，葉數(shù)為3，故BPF=30，同時可以看出一階葉頻脈動壓力幅值遠高于其他倍葉頻處，這也是研究脈動壓力時一般只對一階脈動壓力進行研究的原因。通過對比槳葉徑向P1、P2、P3三點的脈動特性，得知P1和P2之間壓力幅值衰減速度明顯大于P2和P3之間的衰減速度，這是由于螺旋槳徑向作用范圍有限，距離槳葉較近范圍，螺旋槳影響程度隨距離增加迅速減弱，向外延伸，流場逐漸進入穩(wěn)定狀態(tài)。P4和P5兩點BPF幅值達到了300 Pa以上，明顯高于徑向上的P1點，說明該處受螺旋槳旋轉(zhuǎn)作用影響很大。P6和P7兩點較P4和P5兩點衰減相當(dāng)顯著，尤其是P7的BPF壓力幅值已降低到了90 Pa左右，說明沿槳葉軸向向后移動，槳葉尾流影響作用逐漸減小。對比P4～P7，可以看出，靠近槳轂的P5、P7兩點壓力幅值要小于各自同軸向位置的P4、P6兩點，說明P5、P7兩點受螺旋槳旋轉(zhuǎn)作用影響相對也略小。

圖8 特征點的脈動壓力幅值Fig.8 Pulsating pressure amplitude of feature points

4 螺旋槳噪聲數(shù)值計算

通過大渦模擬獲得穩(wěn)定的非定常流場，創(chuàng)建聲學(xué)網(wǎng)格，包括聲源區(qū)和聲傳播區(qū)，其中聲源區(qū)取自CFD計算區(qū)域，略小于CFD計算區(qū)域。聲場的外邊界是一層無限元邊界面，無限元可以傳播聲而不反射聲，因此較好地解決了聲音反射為計算帶來誤差的問題。同時為了有利于傳播聲波，聲學(xué)網(wǎng)格尺寸滿足每波長至少6個網(wǎng)格節(jié)點。CFD節(jié)點和聲學(xué)網(wǎng)格節(jié)點之間通過保守整合法進行信息傳遞以保證計算精度。

取3個距離槳中心5D距離的特征點分別位于軸向槳前和槳后、徑向，3點的頻譜曲線如圖9所示，由于時間步長設(shè)置為0.002 5，對應(yīng)有效頻率上限為200 Hz，參考聲壓為 1 μPa。

圖9 特征點的聲壓頻譜曲線Fig.9 Sound pressure spectrum at feature points

由圖9可以看出，螺旋槳各特征點的頻譜曲線能夠反映出葉頻BPF信息，且在各階葉頻處都存在聲壓峰值。徑向上槳前槳后兩點的聲壓頻譜曲線幾乎重合。徑向與軸向相比，各階葉頻處聲壓峰值相差不大，但在一階葉頻內(nèi)的較低頻段，徑向位置聲壓明顯小于軸向位置。

對應(yīng)頻譜曲線中特殊頻率，給出各頻率時通過螺旋槳軸線的截面聲壓云圖，如圖10所示。圖10(a)以5 Hz聲壓云圖為例，展示了在一階葉頻內(nèi)頻段的聲壓云圖，該頻段偶極子噪聲尤為顯著，聲壓分布呈現(xiàn)出橫8字形特征，螺旋槳軸向聲壓明顯高于徑向聲壓。圖10(b)～(d)為螺旋槳各階葉頻 BPF=30、60、90 Hz處的聲壓分布云圖，觀察螺旋槳周圍，可以發(fā)現(xiàn)聲壓分布同樣呈8字形，這也體現(xiàn)出螺旋槳周圍流場的脈動壓力幅值大小，軸向位置受槳葉旋轉(zhuǎn)影響很小，壓力幅值明顯小于槳葉周圍其他位置，此規(guī)律在一階葉頻BPF處表現(xiàn)尤其明顯。同時可以看出在二階葉頻和三階葉頻處8字形明顯減小，槳葉后方聲壓高于槳葉前方，說明后方受槳葉旋轉(zhuǎn)作用影響要強于前方。總體來說，可以看出螺旋槳噪聲主要由單極子噪聲組成，低頻段內(nèi)比重較大的偶極子噪聲在各階葉頻處相對于單極子噪聲要小很多，在云圖中已經(jīng)無法體現(xiàn)，這與圖9頻譜曲線表述是一致的。

圖10 螺旋槳噪聲的聲壓云圖Fig.10 Sound pressure contours of propeller noise

在通過螺旋槳軸線平面上距離槳中心5D距離取32個測點平均分布于螺旋槳周圍，根據(jù)總聲壓級計算:

式中:SLi是第i個1/3oct中心頻率點處聲壓級。

計算得到各點處總聲壓級，如圖11所示。可見，螺旋槳周圍5D處各點聲壓級大小相近，符合均勻來流單極子聲源的特性，各點最大聲壓級為112.3 dB，最小聲壓級為111.3 dB，誤差可能是流體網(wǎng)格和聲學(xué)網(wǎng)格插值時產(chǎn)生的，其差值僅為1 dB。

圖11 螺旋槳聲指向性(距離中心5D)Fig.11 Noise directivity of propeller at distance 5D

5 結(jié)論

本文針對NREL5MW風(fēng)力機設(shè)計了Spar式風(fēng)力機平臺，并采用頻域水動力分析方法對Spar式風(fēng)力機平臺的水動力性能進行數(shù)值模擬，分析平臺有無縱蕩板、平臺的重本文采用流-聲耦合方法對螺旋槳均勻流噪聲進行了分析，進而得出以下結(jié)論:

1)通過分析敞水性能、槳葉表面壓力分布，驗證LES湍流模型結(jié)合全結(jié)構(gòu)化網(wǎng)格技術(shù)的計算方案具有較好的準(zhǔn)確性;

2)通過對槳葉周圍流場的脈動壓力計算，可以顯示出各階葉頻處的峰值，均勻來流時的一階葉頻(BPF)脈動壓力幅值明顯高于其他葉頻;

3)通過計算得到的槳葉周圍特征點的聲壓頻譜曲線，可以看出，在各階葉頻處存在聲壓峰值，且在一階葉頻處達到最大值，總聲壓級大小和該值大小相近;

4)通過計算得到的聲壓云圖分析，驗證了均勻流螺旋槳噪聲主要由單極子噪聲和偶極子噪聲組成，其中單極子噪聲占主要部分。

通過本文大渦模擬計算螺旋槳水動力性能，結(jié)合Actran聲學(xué)軟件對螺旋槳均勻流噪聲進行數(shù)值模擬，為進一步非均勻流噪聲的數(shù)值計算打下了基礎(chǔ)。

［1］ROSS D.Mechanics of underwater noise［M］.New York:Pergamon Press，1976:113-119.

［2］朱錫清，李亞，孫紅星.船舶螺旋槳葉片與艉部湍流場互作用噪聲的預(yù)報研究［J］.聲學(xué)技術(shù)，2006，25(4):361-364.ZHU Xiqing，LI Ya，SUN Hongxing.Ship propeller blades and stern turbulence field interaction noise prediction research［J］.Journal of Acoustic Technology，2006，25(4):361-364.

［3］朱錫清，唐登海，孫紅星，等.船舶螺旋槳低頻噪聲研究［J］.水動力學(xué)研究與進展(A輯)，2000(1):74-81.ZHU Xiqing，TANG Denghai，SUN Hongxing，et al.Ship low-frequency noise of propeller study［J］.Water Dynamics Research and Progress(A)，2000(1):74-81.

［4］熊紫英，朱錫清，劉小龍，等.船尾伴流場-導(dǎo)管-螺旋槳互作用噪聲預(yù)報研究［J］.聲學(xué)學(xué)報，2009(2):117-123.XIONG Ziying，ZHU Xiqing，LIU Xiaolong，et al.Stern wake field-catheter-interaction noise prediction research［J］.Acta Acustica，2009(11):117-123.

［5］SEOL H，JUNG B，SUH J C，et al.Prediction of non-cavitation underwater propeller noise［J］.Journal of Sound and Vibration，2002，257(1):131-156.

［6］SEOL H，SUH J C，LEE S.Development of hybrid method for the prediction of underwater propeller noise［J］.Journal of Sound and Vibration，2005，288:345-360.

［7］楊瓊方，王永生，曾文德，等.大側(cè)斜螺旋槳負載噪聲的邊界元數(shù)值聲學(xué)方法頻域內(nèi)計算分析［J］.兵工學(xué)報，2011(9):1118-1125.YANG Qiongfang，WANG Yongsheng，ZENG Wende，et al.Big side of the noise of the propeller load boundary element numerical acoustic method in the frequency domain analysis［J］.Acta Armamentarii，2011(9):1118-1125.

［8］楊瓊方，王永生，魏應(yīng)三，等.基于URANS模擬的伴流場中大側(cè)斜槳無空化噪聲的時域和頻域預(yù)報［C］//第十三屆船舶水下噪聲學(xué)術(shù)討論會論文集.北京，2011:325-338.

［9］KOBAYASHI T.Large eddy simulation for engineering applications［J］.Fluid Dynamics Research，2006，38:84-107.

［10］HOLM D D，GEURTS B.Commentator errors in large-eddy simulation［J］.Journal of Physics A:Mathematical and General，2006，39:2213-2229.

［11］FELTEN F，F(xiàn)AUTRELLE Y，TERRAIL Y D，et al.Numerical modelling of electromagnetically-driven turbulent flows using LES methods［J］.Applied Mathematical Modelling，2004，28(1):15-27.

［12］ESCOBAR M.Finite element simulation of flow-induced noise using Lighthill＇s acoustic analogy［D］.Bayen:University Erlangen Nurnberg，2007:36-57.

［13］KOYAMA K.Comparative calculations of propeller by surface panel method-workshop organized by 20th ITTC propulsion committee［C］//Papers of Ship Research Institute，(s.l.)，1993:179-186.

［14］譚廷壽.非均勻流場中螺旋槳性能預(yù)報和理論設(shè)計研究［D］.武漢:武漢理工大學(xué)，2003:53-66.TAN Tingshou.Theory of propeller operating in non-uniform flow field performance prediction and design research［D］.Wuhan:Wuhan University of Science and Engineering，2003:53-66.

［15］STAINER M J.The application of a RANS code to model propeller DTRC 4119［R］.Grenoble:ITTC report，1998.

［16］王超.螺旋槳水動力性能、空泡及噪聲性能的數(shù)值預(yù)報研究［D］.哈爾濱:哈爾濱工程大學(xué)，2010:78-79.WANG Chao.Propeller hydrodynamic performance and cavitation and noise performance of the numerical prediction research［D］.Harbin:Harbin Engineering University，2010:78-79.