連續擺動吊艙推進器水動力性能數值模擬

胡健， 趙旺， 王子斌， 王雅楠， 張維鵬

(哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱 150001)

吊艙推進器是一種新型的電力推進形式，它集螺旋槳和操舵裝置于一體，使得船舶具有更好的操縱性能，吊艙推進器應用于船舶推進器領域，吊艙推進器在斜流中，所受水動力載荷會發生顯著變化，尤其在瞬時回轉工況下，水動力載荷變化更加劇烈，需要密切關注。

關于吊艙推進器的水動力性能研究，理論方法包括基于勢流理論的升力面法和面元法，Bal[1]分析了吊艙單元周圍的流動情況，并對螺旋槳在吊艙上的性能特性進行了研究。楊晨俊等[2]采用單槳吊艙推進器的定常計算方法，計算了吊艙與槳葉之間的相互影響。胡健等[3]用迭代方法求解了吊艙和螺旋槳之間的相互影響，并且還研究了在船后伴流中推進器的水動力性能。Liu等[4]采用面元法預測了不同方位角下吊艙推進器的非定常力、扭矩和彎矩。

國內外很多學者用CFD法做了很多工作，CFD 方法的基本思想是：用一系列連續的多面體網格分割原來連續的計算域，用網格點的變量值表示原來位置的物理量，再運用物理場中的控制方程，通過不斷迭代求解得到數值解，作為物理變量的近似解[5]。Shamsi等[6]基于雷諾平均Navier-Stokes(reynolds averaged Navier-Stokes,RANS)的求解器研究了不同角度下的吊艙推進器的水動力特性的變化。王智展等[7]結合RANS法計算了吊艙槳的瞬時推力系數和扭矩系數，并計算了吊艙單元水動力系數隨回轉角的變化。羅曉園[8]采用面元法計算了3種不同縱斜螺旋槳的敞水特性，并結合CFD計算了不同靜水系數下吊艙推進器的水動力性能。常欣等[16]研究了斜流中螺旋槳的非定常水動力性能，結果表明，斜流角度越大，單槳葉受力的脈動幅度越大。徐嘉啟[17]等研究了HCRSP推進器在操舵工況的空泡性能，通過吊艙后槳與前槳的組合用以展示空泡特征。張維鵬等[18]等對斜流工況中槳舵的干擾過程進行了分析，結果表明舵表面的脈動壓力受斜流和尾渦的雙重影響。孫聰等[19]基于分離渦模擬方法研究了導管槳在斜流中的水動力性能，結果表明斜向來流使得槳葉水動力載荷非定常性增強。Dubbioso等[20]基于雷諾時均方程研究了10°～30°斜流角下螺旋槳的水動力性能，研究結果表明槳的推力隨斜流角的增大而增大。Li等[21]對泵噴推進器的尾渦進行了比較研究，結果表明，混合RANS/LES方法能夠捕捉到豐富的湍流特征。

很多學者從試驗角度研究了吊艙推進器的水動力性能。Maciej[9]采用試驗方法對吊艙推進器在-45°～45°的斜流角下的水動力性能開展了研究，并測量了不同進速系數下的水動力特性。Islam[10]將數值預報方法與實驗評價相結合，用于探索吊艙推進器的水動力性能。趙大剛等[11]通過試驗研究了L型吊艙推進器在動態操舵時的螺旋槳載荷，研究表明：右旋槳L型吊艙推進器在向右操舵時螺旋槳整體上的推力要較靜態操舵狀態的大。沈興榮等[12]通過模型試驗研究了拖式吊艙推進器±30°舵角時的水動力性能,分析了舵角和進速系數對吊艙推進器水動力性能的影響。張志榮等[13]應用數值方法計算了直航條件下推式和拉式吊艙推進器的水動力性能，還給出了斜航條件下吊艙推進器的流場模擬，并將其推力、扭矩等與試驗值作對比。謝清程等[14]測量了0°～180°方位角下吊艙推進器的水動力特性，探索了吊艙推進裝置研制中可能遇到的技術難點。

上述研究的重點主要是在穩定斜流工況下吊艙推進器的水動力性能。吊艙推進器的主要貢獻之一是提高了船舶的操縱性。因此，有必要對全方位吊艙推進器的瞬態水動力載荷進行更詳細的分析。本文使用數值模擬軟體STAR-CCM+，進行數值模擬。在數值模擬廣泛收斂性研究的基礎上，分析了吊艙推進器產生的力和力矩。同時，還研究了螺旋槳載荷與水流入射角的關系。

1 吊艙推進器模型及計算域

1.1 控制方程

在選擇物理模型時，假設流體恒密度不可壓縮，忽略質量力，采用時間平均法建立了雷諾湍流平均方程。

連續性方程為：

(1)

動量方程為：

(2)

本文計算選用的模型為SSTk-ω湍流模型，它考慮了湍流剪應力的輸運特性，能夠準確預報由于逆壓梯度導致的流動分離點和分離區域。SSTk-ω湍流模型方程：

k的運輸方程：

(3)

ω的運輸方程：

(4)

式中：u為速度；x為坐標軸(i，j=1，2，3分別表示x、y、z3個空間坐標)；k為湍流動能；ω為比耗散率。F1、β、γ、σk、σω均為模型參數;β*為模型常數，取0.09。

在式(3)和式(4)中，雷諾應力的渦粘性模型為：

τtij=2μt(Sij-SnnSij/3)-2ρkSij/3

(5)

式中：μt=ρk/ω為渦粘性;Sij為平均速度應變率張量;Snn為克羅內克算子。生成項Pω為：

Pω=2γρ(Sij-ωSnnSij/3)Sij

(6)

1.2 幾何模型

吊艙推進器包括螺旋槳、支架和吊艙，計算時所用的螺旋槳為四葉定距槳，旋向為右旋，相關參數如表1所示，圖1給出了吊艙推進器吊艙和支架的尺寸參數。

表1 螺旋槳主要參數Table 1 Main parameters of propeller

圖1 吊艙和支架的幾何參數Fig.1 Geometric parameters of pod and support

1.3 計算域劃分

所設置的計算域能夠模擬吊艙推進器所處的工作環境，計算域應盡量設置足夠大，可以避免邊界水波反射對尾波系造成影響，速度進口和壓力出口應與吊艙推進器保持一定的距離，這樣能夠得到均勻的入射流和槳后發育完全的流場。

設吊艙螺旋槳的直徑為D，以下用螺旋槳直徑作為度量長度的單位，用長方體代替吊艙推進器工作的水域，入口距吊艙單元的距離為3D，出口距吊艙單元的距離為6D，考慮到吊艙單元瞬時回轉工況回轉的工況，因此水域側向長度設為8D，設有2個圓柱形旋轉域，螺旋槳旋轉域直徑為1.12D，吊艙單元旋轉域直徑為2.4D。圖2展示的是大域和近場域的位置，定義螺旋槳初始位置的軸向與x軸重合，y軸和z軸的方向如圖2所示。

圖2 計算域Fig.2 Computational domain

1.4 網格劃分

在湍流模型下，棱柱層的設置至關重要，對于螺旋槳而言，y+值應取小于1的值，y+值過大會導致棱柱層過厚，計算結果不準確，通過經驗公式可以求得第1層棱柱層的厚度[15]：

(7)

式中：L以螺旋槳直徑D作為特征長度(L=0.25 m)；一般來說，雷諾數Rn越大，棱柱層厚度越小，此工況下雷諾數為Rn=UrefLref/ν=8.68×105，Lref=D/2=0.125 m，Uref為螺旋槳梢端的線速度，Uref=nπD≈7.85 m/s。除了在吊艙表面設置棱柱層外，為了精細網格，還特別對吊艙特征線及其包裹的特征面進行了加密，生成的網格如圖3所示。

圖3 吊艙推進器表面的網格Fig.3 Mesh of podded propulsor surface

2 網格和時間步驗證

從船艉向船艏看，吊艙推進器向左舷偏轉為正，向右舷偏轉為負。示意圖如圖4所示。

圖4 吊艙推進器受力示意Fig.4 Force sketch of podded propulsor

螺旋槳的推力和轉矩可用無因次化系數來表示，對于螺旋槳的吊艙單元的推力系數和轉矩系數，分別為：

(8)

式中：n為螺旋槳的轉速；i=x,y,z分別代表x向、y向和z向；Ti和Qpi分別為吊艙槳i向的受力和力矩；Fi和QU分別為吊艙單元i向的受力和轉舵力矩；KTPi和KQPi代表吊艙槳i向的推力系數和力矩系數；KTUi和KQU代表吊艙單元i向的推力系數和轉舵力矩系數。吊艙單元包括螺旋槳、吊艙和支架3部分。

在模擬中，通過改變來流速度來實現進速系數的變化，進速系數定義為J=VA/nD，VA代表來流速度。螺旋槳轉速為n=10 r/s，吊艙單元旋轉域運動的角速度函數為：

(9)

即吊艙單元的在-50°～+50°運動，運動周期為T=20/3 s，吊艙單元所處的初始位置為0°。

由于缺少高斜流角下的實驗數據，因此本研究通過將用不同網格數量和時間步條件下的計算結果作收斂性驗證。

2.1 網格驗證

針對吊艙單元設置3組不同的網格數量，計算域內的其他網格尺寸按照同等百分比去變化，在整體網格數量發生變化時，保證y+值和時間步不變。通常在不同的網格數量下，計算得到的值之間的誤差應在5%以下，特別的，對于3種及以上的情況，求解的不確定度應滿足Gi+2-Gi+1≤Gi+1-Gi，即隨著網格數量增加，所求的解應有逐漸穩定的趨勢。表2中展示了穩定斜流工況中3種網格數量下所求解(J=0.2,θ=0°)，可以看到，在Coarse Grid下，所求解與Medium Grid和Fine Grid下的結果相差較大，Medium Grid和Fine Grid下的結果相差不大，綜合求解精確度和求解速度來考慮，最后網格條件選用Medium Grid。

2.2 時間步收斂性分析

時間步的選擇對計算結果的收斂與否至關重要，驗證過程中保證網格條件不變(Medium Grid)，根據螺旋槳的轉速選取了3個時間步長，分別對應每一時間步螺旋槳旋轉了1°、2.5°和5°，分別對應Δt=2.778×10-4，Δt=6.944×10-4，Δt=1.388 9×10-3，以在瞬時回轉工況下螺旋槳的y向力矩系數曲線為例，如圖5所示。

表2 網格獨立性分析Table 2 Mesh independency

圖5 不同時間步長下螺旋槳y向力矩系數Fig.5 y direction moment coefficient of propeller at different time steps

選取圖5中低斜流角(-2.5°～2.5°)和高斜流角(40°～45°)2種情況觀察。

在轉向角為-2.5°～2.5°的過程，在圖6中可以看到3條曲線幾乎重合；再選取40°和45°過程，可以看到Δt=1.388 9×10-3的曲線與其他2條曲線相差較大，另外2個時間步所繪曲線雖不完全重合，但誤差在5%以下，考慮到是在±50°這種極端工況下，這樣的誤差是可以接受的。兼顧運算速度以及計算精確度，最終選取的時間步長為Δt=6.944×10-4。

圖6 2種斜流角條件下的螺旋槳y向力矩系數Fig.6 y direction moment coefficient of propeller under the condition of two oblique flow angles

3 操縱工況載荷計算

3.1 進速系數對吊艙性能的影響

圖7展示了螺旋槳載荷隨著進速系數的變化情況。在圖7中可以看到KTPx自中間向兩邊逐漸增大，這是因為隨著斜流角的增大，經過螺旋槳軸向的水流速度越小，進而使得螺旋槳推力增加；還可以發現在同一斜流角度下，隨著進速系數增加，KTPx在減小，這是由于水流速度增加，使得經過螺旋槳軸向的水流速度變大，因此導致KTPx降低。

圖7 螺旋槳x向推力系數Fig.7 x direction thrust coefficient of propeller

3.2 斜流工況

圖8為斜流工況中的單槳葉受力，選取的是單槳葉0°～360° 1個周期內的載荷變化。從圖8中可以看到，KTx曲線有1個波峰和1個波谷，并且隨著斜流角度的增加，KTx曲線變化更為劇烈。圖8中的進速系數選取為J=0.2，斜流角選取5個角度，分別為10°、20°、30°、40°和50°。

圖8 單槳葉x向推力系數Fig.8 x direction thrust coefficient of single propeller blade

圖9展示了吊艙單元在瞬時回轉工況和斜流工況兩種情況下的載荷比較，進速系數選取為J=0.2，選取了吊艙單元1個完整的運動周期從圖9(a)可以看到，在斜流角為-40°～-50°時的曲線規律性不明顯，這是因為吊艙單元在-50°斜流角處轉向導致流場紊亂，從而使得此處力的變化更加劇烈。從圖9(b)、(c)中可以看到，吊艙單元在向左轉和向右轉的過程中載荷并不完全重合，這說明2個過程還是有區別的。對于斜流工況的吊艙單元載荷而言，可以發現，載荷曲線與瞬時回轉工況下的載荷曲線變化趨勢一致，且位于其中間。

圖9 瞬時回轉工況和斜流工況中吊艙單元載荷對比Fig.9 Comparisons of pod unit loads under maneuvering and oblique flow conditions

圖10展示的是在z=0截面處的槳后速度場，分別是在斜流角為-30°和30°下瞬時回轉工況和穩定斜流工況的對比，瞬時回轉工況選取的是從左舷轉向右舷的過程，可以看到2種工況下的速度場存在明顯差別，這是由于瞬時回轉工況下的吊艙單元一直在作回轉運動，導致了不定常入射流的產生。

圖10 在z=0截面上吊艙單元在瞬時回轉工況和穩定斜流工況的速度場對比Fig.10 Comparison of velocity field of pod unit in instantaneous rotary condition and steady oblique flow conditions at z=0

4 結論

1)在穩定斜流工況下，螺旋槳x向載荷隨斜流角增加而增加，y向載荷隨斜流角增加而增加。

2)在瞬時回轉工況下，吊艙單元在向左和向右回轉2種狀態下的螺旋槳載荷不同，以y向推力系數為例，在0°斜流角時，向左轉時其值為-0.01，向右轉時其值為-0.04。

3)在穩定斜流工況和瞬時回轉工況下并不相同，以z向的推力系數為例，在0°斜流角時，穩態載荷為0.035，而吊艙單元在向左轉時載荷為0.03，向右轉時為0.044。