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

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

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

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

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

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

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

1 吊艙推進(jìn)器模型及計(jì)算域

1.1 控制方程

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

連續(xù)性方程為：

(1)

動(dòng)量方程為：

(2)

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

k的運(yùn)輸方程：

(3)

ω的運(yùn)輸方程：

(4)

式中：u為速度；x為坐標(biāo)軸(i，j=1，2，3分別表示x、y、z3個(gè)空間坐標(biāo))；k為湍流動(dòng)能；ω為比耗散率。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)克算子。生成項(xiàng)Pω為：

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

(6)

1.2 幾何模型

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

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

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

1.3 計(jì)算域劃分

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

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

圖2 計(jì)算域Fig.2 Computational domain

1.4 網(wǎng)格劃分

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

(7)

式中：L以螺旋槳直徑D作為特征長(zhǎng)度(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è)置棱柱層外，為了精細(xì)網(wǎng)格，還特別對(duì)吊艙特征線及其包裹的特征面進(jìn)行了加密，生成的網(wǎng)格如圖3所示。

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

2 網(wǎng)格和時(shí)間步驗(yàn)證

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

圖4 吊艙推進(jìn)器受力示意Fig.4 Force sketch of podded propulsor

螺旋槳的推力和轉(zhuǎn)矩可用無因次化系數(shù)來表示，對(duì)于螺旋槳的吊艙單元的推力系數(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部分。

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

(9)

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

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

2.1 網(wǎng)格驗(yàn)證

針對(duì)吊艙單元設(shè)置3組不同的網(wǎng)格數(shù)量，計(jì)算域內(nèi)的其他網(wǎng)格尺寸按照同等百分比去變化，在整體網(wǎng)格數(shù)量發(fā)生變化時(shí)，保證y+值和時(shí)間步不變。通常在不同的網(wǎng)格數(shù)量下，計(jì)算得到的值之間的誤差應(yīng)在5%以下，特別的，對(duì)于3種及以上的情況，求解的不確定度應(yīng)滿足Gi+2-Gi+1≤Gi+1-Gi，即隨著網(wǎng)格數(shù)量增加，所求的解應(yīng)有逐漸穩(wěn)定的趨勢(shì)。表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 時(shí)間步收斂性分析

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

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

圖5 不同時(shí)間步長(zhǎng)下螺旋槳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個(gè)時(shí)間步所繪曲線雖不完全重合，但誤差在5%以下，考慮到是在±50°這種極端工況下，這樣的誤差是可以接受的。兼顧運(yùn)算速度以及計(jì)算精確度，最終選取的時(shí)間步長(zhǎng)為Δ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 操縱工況載荷計(jì)算

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

圖7展示了螺旋槳載荷隨著進(jìn)速系數(shù)的變化情況。在圖7中可以看到KTPx自中間向兩邊逐漸增大，這是因?yàn)殡S著斜流角的增大，經(jīng)過螺旋槳軸向的水流速度越小，進(jìn)而使得螺旋槳推力增加；還可以發(fā)現(xiàn)在同一斜流角度下，隨著進(jì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個(gè)周期內(nèi)的載荷變化。從圖8中可以看到，KTx曲線有1個(gè)波峰和1個(gè)波谷，并且隨著斜流角度的增加，KTx曲線變化更為劇烈。圖8中的進(jìn)速系數(shù)選取為J=0.2，斜流角選取5個(gè)角度，分別為10°、20°、30°、40°和50°。

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

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

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

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

圖10 在z=0截面上吊艙單元在瞬時(shí)回轉(zhuǎn)工況和穩(wěn)定斜流工況的速度場(chǎng)對(duì)比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)在瞬時(shí)回轉(zhuǎn)工況下，吊艙單元在向左和向右回轉(zhuǎn)2種狀態(tài)下的螺旋槳載荷不同，以y向推力系數(shù)為例，在0°斜流角時(shí)，向左轉(zhuǎn)時(shí)其值為-0.01，向右轉(zhuǎn)時(shí)其值為-0.04。

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