LES和無限元耦合方法預報螺旋槳均勻流噪聲

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

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

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

1 流-聲耦合方法理論

1.1 大渦模擬SGS模型

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

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

1.2 Lighthill聲學類比理論

Lighthill考慮的模型為:在無限大的均勻、靜態聲介質中包含一個有限的湍流運動區域V，因此與流動有關的聲源都集中在該區域內。在區域V外，遠離湍流區域的流體中密度的波動和聲波相似，因此整理連續方程和動量方程，并簡化得到遠離湍流區域流體中的勻質聲學波動方程:

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

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

式中，σij表示雷諾應力張量的粘性部分。

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

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

1.3 ACTRAN 聲學原理

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

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

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

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

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

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

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

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

2 計算前處理

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

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

表1 螺旋槳幾何參數Table 1 Geometry coefficients of propeller

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

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

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

2.2 網格劃分及參數設定

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

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

圖3 槳葉外O網Fig.3 O-block out of propeller

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

圖4 螺旋槳網格拓撲結構Fig.4 Mesh topology of propeller

3 大渦模擬計算結果分析

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

3.1 敞水性能計算結果驗證

為便于計算結果表述，定義如下參數:

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

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

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

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

3.2 槳葉表面壓力系數對比分析

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

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

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

3.3 脈動壓力計算分析

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

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

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

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

4 螺旋槳噪聲數值計算

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

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

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

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

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

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

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

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

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

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

5 結論

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

1)通過分析敞水性能、槳葉表面壓力分布，驗證LES湍流模型結合全結構化網格技術的計算方案具有較好的準確性;

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

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

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

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

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

［2］朱錫清，李亞，孫紅星.船舶螺旋槳葉片與艉部湍流場互作用噪聲的預報研究［J］.聲學技術，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］.水動力學研究與進展(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］熊紫英，朱錫清，劉小龍，等.船尾伴流場-導管-螺旋槳互作用噪聲預報研究［J］.聲學學報，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］楊瓊方，王永生，曾文德，等.大側斜螺旋槳負載噪聲的邊界元數值聲學方法頻域內計算分析［J］.兵工學報，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］楊瓊方，王永生，魏應三，等.基于URANS模擬的伴流場中大側斜槳無空化噪聲的時域和頻域預報［C］//第十三屆船舶水下噪聲學術討論會論文集.北京，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，FAUTRELLE 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］譚廷壽.非均勻流場中螺旋槳性能預報和理論設計研究［D］.武漢:武漢理工大學，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］王超.螺旋槳水動力性能、空泡及噪聲性能的數值預報研究［D］.哈爾濱:哈爾濱工程大學，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.