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

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

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

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

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

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

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

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

1.1 控制方程

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

連續(xù)性方程為：

(1)

動量方程為：

(2)

本文計算選用的模型為SSTk-ω湍流模型，它考慮了湍流剪應(yīng)力的輸運特性，能夠準(zhǔn)確預(yù)報由于逆壓梯度導(dǎo)致的流動分離點和分離區(qū)域。SSTk-ω湍流模型方程：

k的運輸方程：

(3)

ω的運輸方程：

(4)

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

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

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

(5)

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

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

(6)

1.2 幾何模型

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

表1 螺旋槳主要參數(shù)Table 1 Main parameters of propeller

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

1.3 計算域劃分

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

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

圖2 計算域Fig.2 Computational domain

1.4 網(wǎng)格劃分

在湍流模型下，棱柱層的設(shè)置至關(guān)重要，對于螺旋槳而言，y+值應(yīng)取小于1的值，y+值過大會導(dǎo)致棱柱層過厚，計算結(jié)果不準(zhǔn)確，通過經(jīng)驗公式可以求得第1層棱柱層的厚度[15]：

(7)

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

圖3 吊艙推進器表面的網(wǎng)格Fig.3 Mesh of podded propulsor surface

2 網(wǎng)格和時間步驗證

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

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

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

(8)

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

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

(9)

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

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

2.1 網(wǎng)格驗證

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

2.2 時間步收斂性分析

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

表2 網(wǎng)格獨立性分析Table 2 Mesh independency

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

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

在轉(zhuǎn)向角為-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向力矩系數(shù)Fig.6 y direction moment coefficient of propeller under the condition of two oblique flow angles

3 操縱工況載荷計算

3.1 進速系數(shù)對吊艙性能的影響

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

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

3.2 斜流工況

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

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

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

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

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

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

4 結(jié)論

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

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

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