斜流中艇后螺旋槳水動力數(shù)值計算方法

張文照，肖昌潤

(1.海軍工程大學,湖北 武漢 430033； 2.海軍后勤技術裝備研究所,北京 100072)

0 引 言

斜流條件下，艇后螺旋槳水動力性能研究一直是進行潛艇操縱性研究的關鍵，目前主要通過模型試驗的方法來進行研究。但模型試驗存在尺度、伴流場的模擬、結(jié)果換算和經(jīng)費等眾多問題的困擾，受到極大限制。數(shù)值方法可以克服上述試驗方法存在的不足，勢流方法是螺旋槳性能分析中常用的數(shù)值方法。在非均勻流中螺旋槳性能分析也常采用勢流方法，比較突出的有Kinnas[1]和Hosino[2]。國內(nèi)有陳家棟[3]、熊鷹[4]、劉小龍[5]等用勢流方法分析給定不均勻流場中螺旋槳性能。

但采用勢流方法處理艇后螺旋槳這種多體結(jié)構(gòu)的水動力問題變得十分困難，并且勢流方法在處理非均勻進流中螺旋槳性能時，需要處理庫塔條件和假定尾渦形狀，有一定的人為因素。而粘流CFD方法則更方便處理復雜的幾何模型，實現(xiàn)對流場的真實模擬。Funeno[6]用粘流CFD方法分析了大側(cè)斜槳非定常性能，WATANABE[7]使用商業(yè)軟件Fluent研究了非均勻伴流中的螺旋槳性能。當前粘性CFD的迅速發(fā)展，使得用CFD方法研究艇后螺旋槳的水動力性能變得越來越容易。本文應用Fluent軟件計算了斜流條件下艇后螺旋槳的水動力，為開展?jié)撏Р倏v性研究奠定了基礎。

1 研究對象

本文所研究的潛艇模型為加拿大DRDC潛艇水動力研究計劃中的通用艇模[8]，模型長4.36 m，直徑0.49 m。螺旋槳模型為5葉槳，直徑0.22 m，盤面比0.5。圖1給出了艇體帶槳的幾何模型。

圖1 計算所用帶槳潛艇模型Fig.1 Submarine model

2 計算方法描述

本文采用基于粘流理論的數(shù)值方法模擬艇后螺旋槳周圍的流動，以商用軟件Fluent為求解器求解不可壓三維RANS方程計算艇后螺旋槳的推力和扭矩。

2.1 控制方程

螺旋槳水動力計算所選取的坐標系為固定在螺旋槳上的旋轉(zhuǎn)直角坐標系，定義坐標軸x為螺旋槳的旋轉(zhuǎn)軸，指向下游為正。控制方程采用Reynolds平均的NS(RANS)方程，其不可壓縮的無量綱形式為：

(1)

式中：ui和uj為絕對速度各分量；p為壓力；Re為雷諾數(shù)；

(4)

式中：νt為湍流運動引起的渦運動粘性系數(shù)；

采用SSTk-ω湍流模型來計算雷諾應力項，以實現(xiàn)控制方程組的封閉。SSTk-ω湍流模型混合了k-ω模型計算近壁面區(qū)域粘性流動的可靠性和k-ε模型計算遠場自由流動的精確性，并在螺旋槳性能計算中反復驗證其較高的精度。

控制方程采用有限體積法離散，有限體積法具有良好的守恒性，能夠采用各種形式的網(wǎng)格以適應復雜的幾何輪廓。控制方程中的對流項采用二階迎風格式離散，速度與壓力的耦合迭代采用SIMPLE 算法。離散后形成的代數(shù)方程用逐點Gauss-Seidel迭代法求解，并且采用代數(shù)多重網(wǎng)格方法加快求解的收斂速度。

2.2 滑移網(wǎng)格法

在進行艇后螺旋槳水動力計算時，螺旋槳旋轉(zhuǎn)，艇體(包含附體)靜止。對于這種問題，F(xiàn)luent軟件有混合面法、多參考系法、滑移網(wǎng)格法和動網(wǎng)格法4種處理方法，前2種方法都假設流動為定常，后2種方法都假定流動為非定常。動網(wǎng)格法需要配合諸如層生成吞并等技術，設置復雜，而滑移網(wǎng)格法在處理定轉(zhuǎn)子問題上更易于實現(xiàn)。本文對于艇后螺旋槳周圍流場的計算采用滑移網(wǎng)格法，包含艇體和附體的部分設置為靜止區(qū)域，螺旋槳部分為滑動區(qū)域，計算網(wǎng)格隨螺旋槳轉(zhuǎn)動。在2個區(qū)域的交界面上物理量的傳遞采用插值方法處理。

2.3 計算域與網(wǎng)格劃分

計算域的進口邊界設置為艇首上游0.8L處(L為艇模長度)的半球面，出口邊界設置在螺旋槳下游2L處，柱面外邊界也為進口，柱面半徑取1L。整個計算域分割成兩大部分，如圖2所示，包含螺旋槳的那一部分區(qū)域為旋轉(zhuǎn)區(qū)域，其他區(qū)域為靜止區(qū)域。

圖2 計算區(qū)域示意圖Fig.2 Computing field

包含螺旋槳的旋轉(zhuǎn)部分由于螺旋槳復雜的幾何外形，劃分為六面體網(wǎng)格較難實現(xiàn)，所以對旋轉(zhuǎn)區(qū)域仍采用四面體網(wǎng)格劃分(見圖3)。

靜止和旋轉(zhuǎn)區(qū)域交界面上的網(wǎng)格劃分見圖4。

圖3 旋轉(zhuǎn)區(qū)域的網(wǎng)格劃分Fig.3 Mesh grid of rotation field

圖4 靜止和旋轉(zhuǎn)區(qū)域交界面處的網(wǎng)格劃分Fig.4 Mesh grid of the interface of static field with rotation field

包含吊艇體、圍殼和十字舵的靜止區(qū)域在進行合理的分塊后較易實現(xiàn)六面體網(wǎng)格劃分。靜止區(qū)域的分塊如圖5所示，分塊會給網(wǎng)格劃分帶來一個問題，即指揮臺圍殼和尾翼表面一層的網(wǎng)格尺寸會延續(xù)到整個計算域的外邊界，然而在計算域的外邊界，其他方向上的網(wǎng)格劃分已經(jīng)很疏，這樣就會使得網(wǎng)格長寬比過大。于是，在網(wǎng)格劃分用如圖6所示的處理方法，即采用“斜向擴展”的方式來防止網(wǎng)格在計算域外邊界上長寬比過大。

圖5 靜止區(qū)域網(wǎng)格劃分Fig.5 Mesh grid of static field

圖6 指揮臺圍殼網(wǎng)格斜向擴展法Fig.6 Oblique extension method

圖7 圍殼和尾翼端面網(wǎng)格劃分Fig.7 Mesh grid of shell and end face

此外，對于圍殼和尾翼端面網(wǎng)格的劃分，將上端面先分成如圖7所示的2個端部面和1個中間面，當支柱上端面邊的節(jié)點固定后，這樣的分塊能有效控制網(wǎng)格層數(shù)。

2.4 邊界條件和初始條件的設置

計算域的進口為圖2中右側(cè)的圓面和徑向的柱面，進口設置為速度進口，給定來流速度的大小和方向。計算域的出口為左側(cè)的圓面，出口設置為壓力出口，給定壓力大小。螺旋槳、艇體、指揮臺圍殼和尾附體表面設為物面無滑移條件，靜止區(qū)域與旋轉(zhuǎn)區(qū)域相重疊的面設置為交界面，完成2個區(qū)域的數(shù)據(jù)傳遞。用進口速度初始化整個流場。

3 計算結(jié)果的比較與分析3.1 0°攻角不同進速非定常計算結(jié)果

首先分別計算了艇體0°攻角下，螺旋槳艇后進速Jb為0.8，1.0，1.2，1.4，1.6工況下，螺旋槳轉(zhuǎn)動1周過程中推力和扭矩的變化值。選取設計點Jb=1.0，畫出該工況下螺旋槳轉(zhuǎn)動1周所產(chǎn)生的推力和扭矩的變化曲線(見圖8和圖9)。通過曲線可以看出螺旋槳轉(zhuǎn)動1周，推力和扭矩呈現(xiàn)出規(guī)律的周期性，周期數(shù)等于螺旋槳的頁數(shù)。螺旋槳轉(zhuǎn)速通常很高，因此推力和扭矩的這種周期性的脈動相對于潛艇運動而言十分快速，因此在研究潛艇操縱運動時，給定艇體運動狀態(tài)的情況下(如潛艇航速、縱傾角等)，可以忽略螺旋槳推力和扭矩的這種周期脈動，而用一個周期的平均值作為艇體該運動狀態(tài)下螺旋槳產(chǎn)生的推力和力矩計算值。

表1和表2分別列出了螺旋槳艇后進速Jb為0.8，1.0，1.2，1.4，1.6工況下計算得到平均推力系數(shù)Kt和平均扭矩系數(shù)10Kq，并和試驗值做比較，計算誤差在3%以內(nèi)。

圖8 Jb=1.0推力系數(shù)1周內(nèi)的變化Fig.8 Trust coefficient in a period when Jb=1.0

圖9 Jb=1.0扭矩系數(shù)1周內(nèi)的變化Fig.9 Torque coefficient in a period when Jb=1.0

進速JbKt計算值Kt試驗值相對誤差/%080275028019100242024408120211020816140175017123160136013325

表2 不同進速螺旋槳平均扭矩系數(shù)

3.2 艇體攻角對螺旋槳水動力性能影響計算

為了研究艇體攻角對艇后螺旋槳推力和扭矩的影響，計算設計進速Jb=1.0， 艇速u=0.7 m/s，艇體處于0°，5°，10°，15°，30°攻角工況下艇后螺旋槳平均推力、扭矩系數(shù)(見圖10)。由于整個計算對象幾何對稱特性，負攻角與正攻角下的槳模水動力對稱。

與試驗值相對比，計算值的相對誤差隨著攻角的增大而增大，總體控制在7%以內(nèi)，該誤差是由于在大攻角時，艇體流體分離比較嚴重，而計算網(wǎng)格不能有效地捕捉到該現(xiàn)象導致的。

圖10 Jb=1.0時艇后螺旋槳推力、扭矩系數(shù)隨艇體攻角的變化Fig.10 Behind-boat trust of the propeller with different pitch angle when Jb=1.0

從圖10可看出，螺旋槳推力和扭矩是隨著艇體攻角呈現(xiàn)出先減小后增大的趨勢。究其原因，是因為當艇體產(chǎn)生攻角時，螺旋槳進流速度會減小，伴流分數(shù)會增大，當攻角比較小的時候伴流分數(shù)的增大導致的螺旋槳實效進速的提高占主要作用，因此在小攻角時螺旋槳推力和扭矩均隨著攻角的增加而降低。但是，隨著攻角的不斷變大，攻角導致的螺旋槳進流速度的降低越來越嚴重，進而超過伴流分數(shù)減小的影響，此后，螺旋槳的推力和扭矩隨著攻角的增大而有所增加。

4 結(jié) 語

本文采用CFD方法對艇后螺旋槳的水動力性能進行計算，并對計算結(jié)果進行比較和分析后得出結(jié)論如下：

1)本文采用的數(shù)值計算方法對帶攻角的艇后螺旋槳推力和扭矩是有效的；

2)隨著攻角的增大，螺旋槳推力和扭矩的計算誤差有所增加，總體控制在7%以內(nèi)；

3)在同一進速下，艇后螺旋槳推力和扭矩隨著艇體攻角的增大呈現(xiàn)出先減小后增大的趨勢。

[1] KINNAS S A,HSIN C J, KEENAN D.A potential based panel method for the unsteady flow around open and ducted propellers[C].Proceedings of the 8thSymposium on Naval Hydrodynamics,1990.

[2] HOSHINO T.Hydrodynamic analysis of propeller in unsteady flow using a surface panel method[C].Journal of Society of Naval Architectures of Japan,1993.

[3] 陳家棟,董世湯.非定常螺旋槳表面壓力面元法預報計算[J].中國造船，1998(1)：21-26.

[4] 熊鷹.非均勻流中螺旋槳空泡和脈動壓力地數(shù)值和試驗研究[D]，武漢：武漢理工大學，2002.

[5] 劉小龍,王國強.螺旋槳非定常性能預估的面元法[J].船舶力學,2006,10(2).

[6] FUNENO I. Analysis of unsteady viscous flows around a highly skewed propeller(in Japanese)[J]. Kansai Society of Naval Architects,2002,237：39-45.

[7] WATANABE T.Simulation of steady and unsteady cavitation on a marine propeller using a RANS code，F(xiàn)ifth International[J].Symposium on Cavitation (CAV2003)，Osaka, Japan,2003：1-4.

[8] FOURNIER E.Wind tunnel investigation of the DREA po-wered MkI submarine model[R].CANADA:DREA,1993.