波浪中的螺旋槳水動力性能數值分析

張 耕, 姚建喜,b

(武漢理工大學 a. 高性能艦船技術教育部重點實驗室; b. 船海與能源動力工程學院, 武漢 430063)

螺旋槳是目前最常用的船舶推進器,其水動力性能與船舶的快速性和經濟性密切相關,對螺旋槳水動力性能進行預報分析是螺旋槳設計工作的重要內容之一[1].蔡榮泉等[2]采用雷諾平均NS(RANS)方法研究螺旋槳的敞水性能,發現在槳葉的葉背表面,最大負壓區域集中在內部的彎凸部位,葉梢輪廓兩側葉面與葉背間的壓力連續變化.劉志華等[3]采用RANS方法對不同葉數、不同側斜度和不同剖面形式的螺旋槳的敞水性能進行了數值預報,計算結果與試驗數據吻合較好.Yao[4]采用RANS方法對斜流中螺旋槳的水動力性能進行了數值預報,在多種進速系數和攻角組合的情況下,得到了螺旋槳的推力及轉矩,以及螺旋槳盤面在旋轉一周過程中的平均推力分布.王超等[5]采用RANS方法,對敞水中螺旋槳的槳葉表面壓力分布進行了考察與分析,發現槳葉壓力面上的壓力由葉根至葉梢先增大后減小.

以上研究僅考慮螺旋槳在無限水深條件下的水動力性能,而不考慮自由面、波浪等的影響.實際上,這些因素在一定條件下,對船后工作螺旋槳的性能有很大影響.近年來,國內外學者基于水池試驗結果并結合理論分析,針對存在自由面且無波浪條件下的螺旋槳水動力性能進行了研究.例如,Califano等[6]采用RANS方法研究螺旋槳在自由面附近的水動力性能,發現在螺旋槳旋轉過程中,與靜水中不存在自由面的情況相比,推力及轉矩出現周期性振蕩,且單片槳葉的振蕩幅度更大.Kozlowska等[7]研究了浸深對螺旋槳水動力性能的影響,發現當進速系數較大時,螺旋槳推力的數值計算結果與試驗數據較為吻合,而當進速系數較小時,數值計算結果與試驗數據相差較大.張磊等[8]采用RANS方法計算分析了浸深與進速系數對螺旋槳水動力性能的影響,發現當浸深與進速系數較小時,螺旋槳的推力及轉矩易發生突變,且槳葉受力脈動更加劇烈.

此外,船舶在海上航行時,不可避免地會受到波浪的作用,船體會產生搖蕩運動,此時船后工作的螺旋槳可能會刺穿水面,與空氣發生接觸,產生吸氣、飛車等現象,導致螺旋槳的推力及轉矩減小、效率下降.不僅如此,螺旋槳周期性出入水面,會導致槳葉產生負荷脈動變化,進而引起發動機轉速及功率的巨大變化,影響螺旋槳和船舶動力系統軸系的強度.在較大的波浪作用下,螺旋槳性能的惡化還可能導致船體失去控制,威脅船舶航行安全.曹梅亮[9]開展了波浪中不同浸深條件下的螺旋槳水動力性能試驗,發現隨著浸深減小,螺旋槳推力及轉矩的平均值也隨之減小,效率亦會下降,而當浸深較大時,螺旋槳在波浪中的推力及轉矩曲線與無波浪時相同浸深下的曲線很接近.董國祥等[10]開展了一系列波浪中的螺旋槳水動力性能試驗,研究了在浸深相同的條件下,波長對螺旋槳推力、轉矩和效率的影響,發現波浪中螺旋槳的推力增量、轉矩增量在短波段的變化比長波段大,且在波浪中選取不同的進速系數時,推力增量、轉矩增量的變化規律一致.Tokgoz等[11]采用RANS方法,對不同浸深與進速系數條件下波浪中螺旋槳的推力波動現象進行了研究,發現當浸深較大時,數值計算結果與試驗數據吻合良好,而當浸深較小時,誤差較大.盡管如此,與存在自由面且無波浪的靜水條件相比,目前針對波浪中螺旋槳水動力性能的研究還相對較少.

本文采用RANS方法,數值預報存在自由面且無波浪的靜水條件下和波浪中螺旋槳的水動力性能,分析浸深與進速系數對螺旋槳推力、轉矩、槳葉表面壓力分布和水面擾動變形的影響,在現有其他學者的研究基礎上,進一步分析波浪中螺旋槳水動力性能的特點.

1 螺旋槳幾何模型及計算工況

采用基于OpenFOAM的RANS求解器,數值預報螺旋槳在波浪中的水動力性能.研究對象為國際研究標準船模KVLCC2(The Second KRISO Very Large Crude Carrier)使用的螺旋槳,編號為KP458.數值計算時,螺旋槳的縮尺比取100,對應的KVLCC2垂線間長Lpp為3.2 m.螺旋槳的幾何模型如圖1所示.螺旋槳模型的主要參數為:直徑 0.098 6 m;螺距比 0.721 2;盤面比0.431;轂徑比0.155;槳葉數4.

圖1 螺旋槳幾何模型Fig.1 Geometry of propeller

為了將計算結果與文獻[11]中的試驗數據進行比較,驗證計算精度,計算考慮的波浪參數與該文獻中的參數一致,即規則波的波長取1倍船長(3.2 m),波幅取0.019 m,波陡為 0.011 9.計算時,螺旋槳轉速n取30 r/s,但考慮了不同的浸深比與進速系數,具體工況為:J=0.20,0.35,0.50;I/R=1.20,1.60,2.00,3.39,4.00.螺旋槳浸深比(I/R)定義為螺旋槳旋轉中心至無擾動水面的距離I與螺旋槳半徑R之比,如圖2所示.進速系數定義為

圖2 浸深比的定義Fig.2 Definition of immersion depth ratio

(1)

式中:u0為來流速度;D為螺旋槳直徑.

此外,為了對比分析波浪對螺旋槳水動力性能的影響,本文還計算了不考慮波浪的靜水工況,但考慮了自由面的影響,具體工況為:J=0.15,0.20,0.25,0.30,0.35,0.40,0.50,0.60;I/R=1.20,1.60,2.00,4.00.

2 數值方法

2.1 流動控制方程

在空間直角坐標系O-xyz下,描述螺旋槳周圍流體的流動,坐標原點O位于螺旋槳旋轉中心,x軸指向來流方向,z軸垂直向下,如圖3所示.

圖3 坐標系示意圖Fig.3 Definition of coordinate system

基于不可壓縮的牛頓流體假設,螺旋槳周圍的流體流動滿足質量守恒方程(連續性方程)和動量守恒方程(RANS方程),其張量形式如下:

(2)

根據Bousinesq假設,雷諾應力可表示為

(3)

式中:νt為渦黏系數;k為湍流動能;δij為Kronecker系數.

采用SSTk-ω湍流模型[12]近似式(3)中的渦黏系數νt.

采用流體體積分數法(Volume of Fluid, VOF)捕捉自由面,流體體積分數F滿足以下傳輸方程:

(4)

對于單個網格單元,若F=0,則該網格單元充滿空氣;若F=1,則該網格單元充滿水;若0

(5)

式中:下標w和a分別表示水和空氣.

2.2 計算域與邊界條件

數值計算的計算域取一長方體,其具體范圍為-2.5Lpp≤x≤0.5Lpp,-10.14D≤y≤10.14D,-Lpp≤z≤Lpp,如圖4所示.

圖4 計算域示意圖Fig.4 Computational domain

采用在遠場邊界上給定無擾動流場和波面的方法在計算域內數值造波,基于線性波浪理論,入口邊界處(即x=0.5Lpp處的邊界,見圖4)的流速為來流速度與造波流速的線性疊加[14],即

(6)

對應的波面表達式為

ζ=-Asin(ωet+kwx)

(7)

式中:ζ為波面升高.

為了保證計算精度,需要消除下游邊界(即x=-2.5Lpp處的邊界,見圖4)反射波的影響.本文采用阻尼消波法在下游區域消波,在消波區域內RANS方程中的fi≠0,其表達式為

fi=(0,0,d(x)w)

(8)

(9)

式中:w為Ui的第3個分量,即Ui=(u,v,w);α為消波系數;xs為消波區域的起始位置.

在計算域的4個側邊界上,即y=±10.14D、z=±Lpp處,設置滑移邊界條件;而在螺旋槳(槳葉、槳轂)表面上,設置壁面邊界條件.具體邊界條件如表1所示.表中:ω為比耗散率;β1=0.075,為湍流模型系數;yw為到壁面距離;kin和ωin由下式估算:

表1 邊界條件Tab.1 Boundary conditions

(10)

式中:It≈0.05,為湍流度;Cμ=0.09,為湍流模型系數;ls為湍流長度尺度,與網格尺度接近;uin為入流流速.

2.3 離散格式

采用有限體積法(Finite Volume Method, FVM)求解流體動力控制方程,其中RANS方程中的擴散項和對流項分別采用中心差分格式和二階迎風格式進行離散處理,時間項采用一階隱式歐拉格式進行離散處理,湍流方程采用二階迎風格式進行離散處理.壓力速度耦合算法采用PIMPLE算法,即將SIMPLE(Semi-Implicit Method for Pressure-Linked Equations)算法和PISO(Pressure Implicit with Splitting of Operators)算法相結合,利用“預報-校正”的策略,通過迭代的方式將速度場和壓力場的耦合方程組解耦.對于由離散方程得到的線性方程組,采用高斯-賽德爾迭代法迭代求解流速、湍流動能和比耗散率,采用GAMG(Generalized Geometric Multi-Grid)法迭代求解壓力[15].

2.4 滑移網格

采用OpenFOAM中的滑移網格方法模擬螺旋槳的旋轉運動.運用該方法時,計算區域被劃分為包含螺旋槳幾何形狀的旋轉區域和剩余的靜止區域.在計算過程中,隨著時間的步進,旋轉區域圍繞槳軸旋轉,靜止區域和旋轉區域之間存在一對幾何尺寸相同的滑移交界面(見圖5),兩個滑移交界面之間的流動信息(流速、壓力等)采用面積分數法進行插值交換.

圖5 滑移網格(x=0 m)Fig.5 Sliding grid (x=0 m)

3 計算結果分析

3.1 網格生成與依賴性分析

本文使用商業軟件Hexpress生成計算網格.圓柱體旋轉域的直徑為1.3D,高為0.7D;靜止區域為長方體,范圍如前所述:-2.5Lpp≤x≤0.5Lpp、-10.14D≤y≤10.14D、-Lpp≤z≤Lpp,其中Lpp=3.2 m.旋轉域的中心位于靜止區域內部x=0 m、y=0 m、z=0 m的位置,兩者之間存在一對幾何尺寸相同的滑移交界面.本文不采用壁面函數近似邊界層內的流動,因此在生成網格時,將靠近壁面第1層網格與壁面間的無量綱垂直距離,即y+,控制在y+<1.

為了確保計算精度并盡量減少計算時間,首先開展了網格依賴性分析,計算工況為J=0.20、I/R=2.00.使用商業軟件Hexpress生成了3套疏密程度不一的多面體非結構化網格,分別記為coarse、medium和fine.網格生成過程中,在x、y、z方向上使網格單元數量等比例增加.為了保證網格質量,在葉尖及葉邊處,進行了加密處理.3套網格的網格單元總數分別為 357 494、692 314 和 1 350 755,其中旋轉域的網格單元數分別為 269 196、491 250 和 932 883;靜止域的網格單元數分別為 88 298、201 064 和 417 872.圖6為fine網格的旋轉域網格和靜止域網格.

圖6 Fine網格Fig.6 Fine grid

網格依賴性分析時,定義螺旋槳旋轉一周的時間步數量為Ns.為了節省計算資源,在計算域內開始造波時,計算時間步長取0.005 s,直至流場穩定后再減小時間步長,取Ns=200.使用3套網格計算得到的螺旋槳推力系數(KT)及轉矩系數(10KQ)時歷曲線,如圖7所示.可以看出,隨著網格加密,時歷曲線變得接近,其中采用medium網格和fine網格得到的計算結果已經十分接近.

注意到,圖7中給出的是無因次推力及轉矩時歷曲線,其表達形式如下:

(11)

式中:T為推力;Q為轉矩.

采用3套網格計算得到推力系數KT的平均值分別為0.224 7、0.228 6 和 0.229 9,轉矩系數10KQ的平均值分別為 0.249 5、0.251 3 和 0.251 2.由此可見,采用medium網格和fine網格的計算結果已經十分接近,繼續加密網格,計算結果已不會出現顯著的變化.但是,本文為獲得更為細致的流場信息,仍采用fine網格對后續算例進行計算.

采用上述fine網格,對時間步長依賴性進行分析,Ns取100、200和400,對應的時間步長分別為 0.000 333 s、0.000 167 s和0.000 083 s.計算得到的時歷曲線如圖8所示.可以看到,隨著時間步長減小,螺旋槳推力系數及轉矩系數的時歷曲線逐漸接近.采用3個時間步長得到的推力系數計算平均值分別為0.224 5、0.229 9 和 0.231 0,轉矩系數計算平均值分別為 0.247 9、0.251 2 和 0.252 9,Ns取200和400得到的平均推力系數和轉矩系數已經十分接近.若進一步減小時間步長,計算得到的推力系數與轉矩系數的峰值估計會更加接近.

圖8 時間步長依賴性分析Fig.8 Time-step dependency study

通過以上分析,可以看出,對后續工況的計算,采用fine網格且Ns取200,能夠保證足夠的計算精度.

此外,對數值造波、消波的效果進行了考察.圖9給出了工況J=0.20、I/R=2.00在縱剖面y=0.99 m處的理論波形與數值波形.可以看出,在-0.5 m≤x≤1.6 m范圍內(螺旋槳旋轉中心位于x=0 m),計算波形與理論波形符合較好,波幅沒有出現明顯衰減.從x=-0.5 m起(消波起始位置xs=-0.5 m),波幅迅速衰減,在下游靠近出口邊界處,波幅幾乎衰減至0.

圖9 造波消波驗證Fig.9 Verification of wave-making and wave-elimination

3.2 靜水條件下的計算結果與分析

首先對存在自由面且無波浪的靜水條件下螺旋槳的水動力性能進行預報分析.為了驗證計算精度,將推力系數及轉矩系數的計算平均值與 SIMMAN 會議公開的KP458螺旋槳敞水試驗數據進行比較,如表2、表3和圖10所示.可以觀察到:當I/R>1.60時,螺旋槳推力系數及轉矩系數的計算平均值與試驗數據相比較為接近,說明此時自由面對螺旋槳水動力性能沒有顯著的影響;對于I/R=1.20、進速系數較小的情況,推力系數及轉矩系數的計算平均值顯著減小,說明此時自由面對螺旋槳水動力性能影響較大;計算時,沒有考慮螺旋槳空化的影響, 然而I/R=4.00的數值計算結果與試驗數據吻合良好,間接說明所考慮的計算工況沒有出現空化現象.

表2 靜水條件下的推力系數Tab.2 Thrust coefficient under the condition of calm water

表3 靜水條件下的轉矩系數Tab.3 Torque coefficient under the condition of calm water

圖10 靜水條件下的計算結果與試驗數據Fig.10 Computational results and experimental data under calm water conditions

3.3 波浪中的計算結果與分析

3.3.1浸水深度 為了分析浸深對波浪中螺旋槳水動力性能的影響,對一系列工況進行了計算:I/R=1.20,1.60,2.00,3.39,4.00,J=0.50,并將計算流體力學(CFD)計算結果與工程流體力學(EFD)試驗數據進行比較.計算得到的推力系數及轉矩系數時歷曲線,如圖11所示.圖12給出了不同浸深比條件下的水面變形情況,水面變形越大,說明波浪對螺旋槳的水動力性能影響越大.圖13給出了對應工況下的槳葉表面壓力分布.可以看出,浸深比越小,螺旋槳對水面的擾動越顯著;當浸深較小時,距離自由面最近的槳葉葉梢處壓力最小,且葉面表面存在受到負壓的區域,隨著螺旋槳的浸深逐漸增大,葉面表面受到的正壓也逐漸增大,受到負壓的區域逐漸減小.

圖11 推力系數及轉矩系數時歷曲線(J=0.50)Fig.11 Time history curves of thrust coefficient and torque coefficient (J=0.50)

圖12 波谷處水面形態(J=0.50)Fig.12 Disturbance of free surface at wave trough (J=0.50)

圖13 波谷處槳葉表面壓力分布(J=0.50)Fig.13 Distribution of blade surface pressure at wave trough (J=0.50)

注意到,采用無量綱化的壓力系數Cp來表示圖13中螺旋槳槳葉表面的壓力分布情況,其表達式為

(12)

式中:p0為靜壓.

從圖中可以看出,I/R=1.20的推力系數及轉矩系數時歷曲線在部分時間段出現劇烈振蕩,與其他浸深比的情況相比,呈相反的變化趨勢.此時,波谷經過螺旋槳附近,導致槳葉刺穿水面,直接暴露在空氣中,如圖12(a)所示.暴露在空氣中的槳葉無法撥水產生推力,槳葉表面的壓力分布如圖13(a)和圖13(b) 所示,因此推力及轉矩出現了劇烈振蕩且減小.

表4給出了推力系數及轉矩系數的計算平均值,將工況I/R=2.00,3.39的推力系數計算結果與大阪大學拖曳水池中測得的試驗數據進行比較,可以看出,兩種浸深比下的相對誤差分別為0.94%和0.09%,計算精度較高.圖14給出了推力系數平均值隨浸深比變化的情況,可以觀察到:螺旋槳推力隨浸深比的減小而減小,且浸深比越小,螺旋槳推力減小得越快.與靜水條件相比,波浪中螺旋槳推力系數及轉矩系數的平均值均有減小.

表4 推力系數及轉矩系數計算值(J=0.50)Tab.4 Computational results of thrust coefficient and torque coefficient (J=0.50)

圖14 推力系數隨浸深比的變化(J=0.50)Fig.14 Thrust coefficient versus immersion depth ratio (J=0.50)

對進速系數J=0.35的工況進行了對比分析,I/R=1.60,2.00,3.39,4.00,推力系數及轉矩系數的時歷曲線和計算平均值分別如圖15和表5所示.推力系數隨浸深比的變化如圖16所示.將工況I/R=2.00,3.39 的推力系數計算結果與大阪大學拖曳水池中測得的試驗數據進行比較,兩者之間的相對誤差分別為3.19%和4.52%.與靜水條件相比,波浪中螺旋槳推力系數及轉矩系數的平均值均有減小.

表5 推力系數及轉矩系數計算值(J=0.35)Tab.5 Computational results of thrust coefficient and torque coefficient (J=0.35)

圖15 推力系數及轉矩系數時歷曲線(J=0.35)Fig.15 Time history curves of thrust coefficient and torque coefficient (J=0.35)

圖16 推力系數隨浸深比的變化(J=0.35)Fig.16 The change of thrust coefficient with immersion depth ratio (J=0.35)

圖17給出了不同浸深比條件下的水面變形情況,圖18給出了對應工況下的槳葉表面壓力分布.同樣可以看出,浸深比越小,螺旋槳對水面的擾動越顯著.

圖17 波谷處水面形態(J=0.35)Fig.17 Disturbance of free surface at wave trough (J=0.35)

圖18 波谷處槳葉表面壓力分布(J=0.35)Fig.18 Distribution of blade surface pressure at wave trough (J=0.35)

3.3.2進速系數 為了分析進速系數對波浪中螺旋槳水動力性能的影響,對一系列工況進行了計算:J=0.20,0.35,0.50,I/R=2.00.計算得到的推力系數及轉矩系數時歷曲線如圖19所示,可以看出:J=0.20的推力系數計算平均值與試驗數據相比偏大,但兩者的振蕩幅度一致;J=0.50的推力系數計算平均值與試驗數據更加接近,但計算所得時歷曲線的振蕩幅度偏大.推力系數及轉矩系數的計算平均值如表6所示,將推力系數計算結果與試驗數據進行比較,J=0.20,0.35,0.50的相對誤差分別為5.12%、3.13%和0.94%,精度良好.進速系數越大,相對誤差越小,如圖20所示.

表6 推力系數及轉矩系數計算值(I/R=2.00)Tab.6 Computational results of thrust coefficient and torque coefficient (I/R=2.00)

圖19 推力系數及轉矩系數時歷曲線(I/R=2.00)Fig.19 Time history curves of thrust coefficient and torque coefficient (I/R=2.00)

圖20 推力系數計算值與試驗值(I/R=2.00)Fig.20 Computational results and experimental data of thrust coefficient (I/R=2.00)

當I/R=2.00且波谷經過螺旋槳附近時,水面形態與槳葉表面壓力分布分別如圖21和圖22所示.可以觀察到:當浸深比相同時,進速系數越小,螺旋槳對水面的擾動越明顯,槳葉的葉面表面主要受正壓,葉背表面主要受負壓,且進速系數越小,葉面表面受到的正壓越大.

圖21 波谷處水面形態(I/R=2.00)Fig.21 Disturbance of free surface at wave trough (I/R=2.00)

圖22 波谷處槳葉表面壓力分布(I/R=2.00)Fig.22 Distribution of blade surface pressure at wave trough (I/R=2.00)

4 結論

采用基于OpenFOAM的RANS求解器,計算了螺旋槳在靜水和波浪中受到的推力及轉矩,分析了浸深與進速系數對螺旋槳受力特性、水面擾動變形和槳葉壓力分布的影響.通過對比分析,得出以下結論:

(1) 在靜水條件下,當I/R>1.60時,自由面對螺旋槳水動力性能沒有顯著的影響.

(2) 對于存在波浪的情況,隨著浸深逐漸減小,螺旋槳推力及轉矩的平均值也隨之減小,波浪對螺旋槳水動力性能的影響變大,當浸深較小時(I/R=1.20),槳葉會在波谷處刺穿水面,此時推力及轉矩發生突變且顯著減小.

(3) 當浸深較小時,距離自由面最近的槳葉葉梢處壓力較小,且葉面表面存在受到負壓的區域,隨著浸深逐漸增大,葉面表面受到的正壓也逐漸增大,受到負壓的區域逐漸減小.

(4) 與試驗數據相比,數值計算的結果精度較高,相對誤差大多在5%以內.