鄭巢生 張志榮

中國船舶科學(xué)研究中心，江蘇無錫214082

基于OpenFOAM的螺旋槳敞水性能預(yù)報方法

鄭巢生 張志榮

中國船舶科學(xué)研究中心，江蘇無錫214082

為深入研究螺旋槳周圍的粘流問題，基于面向?qū)ο蟮拈_源CFD計算平臺OpenFOAM，選擇DTMB P4119槳作為對象，利用RANS方程計算了槳的敞水性能，并分別考察了網(wǎng)格依賴性和離散格式的影響。同時分析了不同半徑處葉剖面的壓力分布及槳前后的流場速度分布。為保證對流項離散的穩(wěn)定性和高精度性，采用了二階NVD格式Gamma混合差分格式，選取k－wSST模型作為湍流模型，壓力速度耦合采用SIMPLE松弛算法。計算結(jié)果與試驗數(shù)據(jù)吻合較好。研究建立了基于OpenFOAM的螺旋槳敞水性能預(yù)報方法。

粘流；螺旋槳；OpenFOAM；RANS；敞水性能；湍流模型

0 引 言

隨著現(xiàn)代船舶性能的快速發(fā)展，能否準(zhǔn)確預(yù)報螺旋槳性能至關(guān)重要。20世紀(jì)80年以來，基于RANS方程求解的粘流方法開始被廣泛應(yīng)用于模擬螺旋槳周圍的粘流［1-2］。目前，國內(nèi)進行螺旋槳方面的計算主要通過CFD商業(yè)軟件Fluent和自編程序來實現(xiàn)，而OpenFOAM與前兩者相比具有完全開源性、面向?qū)ο笮院土己玫睦^承性等優(yōu)點。作為一款開源CFD計算平臺，憑其靈活高效的C++模板自定義功能，OpenFOAM越來越受到CFD研究人員的關(guān)注。

本文基于OpenFOAM開源平臺，將選用多參考系MRFSimpleFoam穩(wěn)態(tài)求解器進行螺旋槳敞水性能預(yù)報方法的研究。

1 數(shù)值模擬

1.1 控制方程

在慣性參考系中，螺旋槳以設(shè)定的角速度繞固定軸旋轉(zhuǎn)，流場是非定常的。因本文中只存在單一螺旋槳旋轉(zhuǎn)，不存在定子的情況，所以既可以采用整個流場坐標(biāo)建立在槳葉上的單一旋轉(zhuǎn)坐標(biāo)系SRF方法，也可以采用多參考系MRF方法建立N－S方程，即在槳葉附近區(qū)域采用建立在槳葉上的旋轉(zhuǎn)參考系，在槳葉遠場區(qū)域采用慣性參考系，再在兩者交界面上通過將速度換算成絕對速度的形式進行流場數(shù)據(jù)交換。為便于今后考慮螺旋槳和定子同時存在的情況，本文采用MRF方法。

其慣性參考系下的N－S方程為：

其旋轉(zhuǎn)參考系下絕對速度形式的N-S方程［3］為：

式中，UI和UR分別為絕對速度和相對速度；p為壓力；ρ為密度；Ω為旋轉(zhuǎn)角速度；Ω×UI為包含向心力的廣義柯氏力項。

OpenFOAM中多參考系MRFSimpleFoam穩(wěn)態(tài)求解器即根據(jù)式（2）求解旋轉(zhuǎn)參考系下的絕對速度。

1.2 湍流模型

為了求解湍流應(yīng)力項，引入了k－wSST湍流模型［4］，其中，湍動能k和渦量脈動強度w的輸運方程為：

令 ?1表示常數(shù) (α1，β1，σk1，σw1)，?2表示常數(shù)(α2，β2，σk2，σw2)，? 表示常數(shù) (α，β，σk，σw)，它們之間關(guān)系為：

式中，α1=5/9；α2=0.44；β1=3/40；β2=0.082 8；σk1=0.85；σk2=1；σw1=0.5；σw2=0.856；β*=9/100。

1.3 離散方法

采用有限體積法求解任意多面體非結(jié)構(gòu)網(wǎng)格離散下的RANS方程。選取的控制體單元如圖1所示，圖中P為控制體單元中心，N為鄰近單元中心，f為單元界面。

圖1 有限體積離散Fig.1 The finite volume discretization

對流項在控制體上的積分和線性化如下：

式中，面場?f采用Jasak提出的非結(jié)構(gòu)網(wǎng)格下的Gamma混合差分格式：

通過權(quán)重系數(shù)γ混合二階中心差分(?f)CD和一階迎風(fēng)差分(?f)CD，其表達式分別為：

由文獻［5］可知 γ∈（0 1），且由式（6）可知 γ越大，精度越高，但穩(wěn)定性越差；γ越小，精度越低，但穩(wěn)定性越好。本文中γ取值為0.2。

擴散項在控制體上的積分和線性化如下：

正交網(wǎng)格時，控制體單元中心P和N的位移向量d與平面 f正交，即平行于 Sf時，面梯度為隱式離散。

非正交網(wǎng)格時，先通過中心差分單元中心值得到單元中心梯度，再插值單元中心梯度引入顯式補償項。

在方程組求解過程中，壓力與速度耦合采用SIMPLE松弛算法，其中 p，U，k和w的松弛因子分別為0.3，0.5，0.5和0.5。

1.4 幾何模型

本文選取DTMB P4119槳作為計算對象。它是一個無傾斜、無縱傾的三葉槳，槳葉直徑D＝0.314 8 m（圖2）。

圖2 DTMB P4119槳Fig.2 DTMB P4119 propeller

由于槳繞軸旋轉(zhuǎn)具有周期性，因此文中采用的周期邊界條件只需計算單槳流道，計算區(qū)域為如圖3所示的1/3圓柱區(qū)域：上游入口位于槳前3.6D處，下游出口位于槳后5D處，側(cè)邊位于槳軸心3.8D處。

圖3 計算區(qū)域Fig.3 Computational region

1.5 網(wǎng)格生成

由于采用多參考系，利用GAMBIT分塊劃分網(wǎng)格時，將旋轉(zhuǎn)參考系下的區(qū)域劃分為直徑d＝1.08D，高h＝0.86D的1/3圓柱區(qū)域，采用非結(jié)構(gòu)網(wǎng)格，并在槳葉表面均勻加密。慣性參考系下的區(qū)域采用結(jié)構(gòu)網(wǎng)格，并向外逐漸稀疏。整個區(qū)域網(wǎng)格總數(shù)為1 086 986（圖4）。

圖4 分塊網(wǎng)格劃分Fig.4 Multi-block grid partition

為了方便OpenFOAM自動捕捉和處理滑移壁面及旋轉(zhuǎn)邊界條件，將旋轉(zhuǎn)部分和固定部分分別命名為rotor和stator，然后再通過終端輸入fluent?MeshToFoam即可。

1.6 邊界條件

在OpenFOAM中，變量的初始化需要在每個以變量名命名的文本中進行，其初始化分別如下：

速度U ：入口為固定值U＝（5.078 0 0），出口為零梯度，槳葉、槳轂和側(cè)邊為固定值U＝（0 0 0）；

壓力 p：入口為零梯度，出口為固定值 p＝0，槳葉、槳轂和側(cè)邊為零梯度；

湍動能 k：入口為固定值 k=3.87×10-5，出口為零梯度，槳葉、槳轂和側(cè)邊采用壁面函數(shù)；

渦量脈動強度w：入口為固定值w＝38.5，出口為零梯度，槳葉、槳轂和側(cè)邊采用壁面函數(shù)。

其中，湍動能k和渦量脈動強度w的初始值分別由湍流密度I和渦粘比 μt/μ得到。

式中，I＝0.1%；μt/μ=1。

2 計算結(jié)果與分析

2.1 敞水性能

分別計算進速系數(shù) J＝0.5，0.7，0.833，0.9和1.1工況下DTMB P4119槳的敞水性能，計算結(jié)果如圖 5和表 1 所示，并與 Jessup［6］的試驗結(jié)果相比較。其中：

式中，U為入口速度值；T和Q分別為槳推力值和扭矩值；ρ為密度；n為轉(zhuǎn)速；D為螺旋槳直徑。

圖5 DTMB P4119槳敞水性能曲線Fig.5 Open-water performance curves of DTMB P4119 propeller

表1 DTMB P4119槳敞水性能比較Tab.1 Comparisons of the open-water performance of DTMB P4119 propeller

由圖5可看出，OpenFOAM計算的槳的推力系數(shù)比實驗值略大，而扭矩系數(shù)比試驗結(jié)果略小一些。

從表1的具體數(shù)值看，推力系數(shù)與試驗值的誤差非常小，最大誤差僅為2.5%，而扭矩系數(shù)與試驗值的誤差相對較大一些。造成較大誤差的主要原因可能是混合差分格式中的權(quán)重系數(shù)γ取值較低，即應(yīng)取較大的值，但同時還要兼顧計算的穩(wěn)定性問題。

2.2 網(wǎng)格依賴性

為了研究網(wǎng)格的依賴性，選取更密的網(wǎng)格計算不同工況下的敞水性能（圖6），即在槳葉表面加密，得到網(wǎng)格總數(shù)為1 614 109。計算表明，兩種網(wǎng)格得到的結(jié)果幾乎相同，但為了節(jié)省計算時間，本文將選取較粗的網(wǎng)格完成所有計算。

2.3 離散格式影響

為了研究離散格式的影響，選用一階迎風(fēng)格式計算不同工況下的敞水性能，并與試驗值相比較。從圖7可以看出，與Gamma格式相比，一階迎風(fēng)格式下的扭矩系數(shù)相對較小，而推力系數(shù)的誤差則相對較大。

圖6 網(wǎng)格依賴性比較Fig.6 Comparisons of grid dependence

圖7 離散格式比較Fig.7 Comparisons of discretization scheme

2.4 表面壓力分布

Jessup利用LDV測量的速度場，通過伯努利方程換算得到不同半徑上槳葉表面的壓力分布。圖8～圖10分別給出了設(shè)計工況 J＝0.833時OpenFOAM計算得到的壓力系數(shù)分布，并與經(jīng)換算得到的結(jié)果進行了比較。圖中，橫坐標(biāo)為x/C（x為葉剖面上的點沿弦長方向到導(dǎo)邊的距離，C為葉剖面處弦長），縱坐標(biāo)Cp為壓力系數(shù)，定義如下：

圖8 r/R＝0.3處葉剖面壓力系數(shù)分布Fig.8 Pressure coefficient distributions on blade section atr/R＝0.3

圖9 r/R＝0.7處葉剖面壓力系數(shù)分布Fig.9 Pressure coefficient distributions on blade section atr/R＝0.7

圖10 r/R＝0.9處葉剖面壓力系數(shù)分布Fig.10 Pressure coefficient distributions on blade section atr/R＝0.9

式中，P為靜壓；P0為遠場參考壓力；U為入口速度值；ρ為密度；Ω為旋轉(zhuǎn)角速度值；r為葉剖面處半徑。

由圖8～圖10可看出，OpenFOAM計算得到的不同半徑葉剖面上的壓力系數(shù)分布與Jessup的試驗結(jié)果十分吻合，因此能準(zhǔn)確預(yù)報推力系數(shù)和扭矩系數(shù)。經(jīng)比較可以看出，導(dǎo)邊附近的壓力系數(shù)比試驗值略小，這主要可能是由于導(dǎo)邊附近的網(wǎng)格不夠密。在葉根附近，r/R＝0.3處的葉剖面壓力系數(shù)分布與Jessup的試驗結(jié)果差別相對較大，這可能是因為Jessup在用伯努利方程換算時假設(shè)流體為理想流體，而沒有考慮葉根處粘性及漩渦的影響。

2.5 流場分布

圖11～圖18分別給出了設(shè)計工況J＝0.833時槳盤面前方x＝-0.3R處和槳盤面后方x＝0.328 1R處軸向、徑向和切向速度的周向平均分布。圖中，橫坐標(biāo)為r/R，縱坐標(biāo)用入口速度無量綱化：

式中，U為入口速度值；Ux、Ur和Ut分別表示周向平均后的軸向、徑向和切向速度值。

圖11 槳盤面前x＝-0.3R處軸向速度的周向平均值Fig.11 Circumferentially-averaged value of axial velocity atx＝-0.3R

圖12 槳盤面前x＝-0.3R處徑向速度的周向平均值Fig.12 Circumferentially-averaged value of radial velocity atx＝-0.3R

圖13 槳盤面后x＝0.328 1R處軸向速度的周向平均值Fig.13 Circumferentially-averaged value of axial velocity atx＝0.328 1R

圖14 槳盤面后x＝0.328 1R處徑向速度的周向平均值Fig.14 Circumferentially-averaged value of radial velocity atx＝0.328 1R

圖15 槳盤面后x＝0.328 1R處切向速度的周向平均值Fig.15 Circumferentially-averaged value of tangential velocity at x＝0.328 1R

圖16 槳盤面后x＝0.328 1R，r＝0.7R處軸向速度的周向分布Fig.16 Circumferential distributions of axial velocity at x＝0.328 1R，r＝0.7R

圖17 槳盤面后x＝0.328 1R，r＝0.7R處徑向速度的周向分布Fig.17 Circumferential distributions of radial velocity at x＝0.328 1R，r＝0.7R

圖18 槳盤面后x＝0.328 1R，r＝0.7R處切向速度的周向分布Fig.18 Circumferential distributions of tangential velocity at x＝0.328 1R，r＝0.7R

由圖11～圖15可看出，OpenFOAM計算得到的槳盤面前方 x＝-0.3R處及槳盤面后方 x＝0.328 1R處軸向速度和徑向速度的分布與試驗結(jié)果基本吻合。

圖16～圖18分別給出了槳盤面后方 x＝0.328 1R，r/R=0.7處軸向、徑向和切向速度的周向分布。OpenFOAM的計算結(jié)果能較好地預(yù)報速度的變化趨勢，但沒有精確捕捉到速度峰值，這可能是由于此處網(wǎng)格不夠密而導(dǎo)致數(shù)值耗散過大［7］。要捕捉峰值變化，可以通過試驗測量到的伴流峰值分布，也可以根據(jù)勢流方法中假設(shè)的尾渦模型，在相應(yīng)的位置布置十分精細(xì)的網(wǎng)格。

3 結(jié) 語

本文基于開源CFD計算平臺OpenFOAM，計算了不同工況下DTMB P4119槳的敞水性能、槳葉壓力分布和流場速度分布，并考察了網(wǎng)格依賴性和離散格式的影響，計算結(jié)果與試驗數(shù)據(jù)能很好地吻合，建立了基于OpenFOAM的螺旋槳敞水性能預(yù)報方法，為今后研究更復(fù)雜的粘流問題奠定了基礎(chǔ)。

［1］TANG D H，CHEN J D，ZHOU W X.Comparative cal?culations of propeller performance by RANS/Panel Method［C］//22nd ITTC Propulsion Committee，Propel?ler RANS/Panel Method Workshop，Grenble，1998.

［2］CHEN J D，ZHOU W X，TANG D H，et al.Comparative calculations of unsteady propeller performance by pan?el method［C］//22nd ITTC Propulsion Committee，Pro?peller RANS/Panel Method Workshop，Grenble，1998.

［3］See the MRF development［EB/OL］.［2011-11-14］.http：//openfoamwiki.net.index.php.See_the_MRF_de?velopment.

［4］MENTER F R.Two-equation eddy-viscosity turbu?lence models for engineering applications［J］.AIAA Journal，1994，32（8）：1598-1605.

［5］JASAK H，WELLER H G，COSMAN A D.High resolu?tion NVD differencing schemes for arbitrarily unstruc?tured meshes［J］.International Journal for Numerical Methods in Fluids，1999，31：431-449.

［6］JESSUP S D.Experimental data for RANS calculations and comparisons［C］//22nd ITTC Propulsion Commit?tee，Propeller RANS/Panel Method Workshop，Gren?ble，1998.

［7］韋喜忠，黃振宇，洪方文.基于非結(jié)構(gòu)網(wǎng)格的螺旋槳周圍流場大渦模擬［J］.水動力學(xué)研究與進展，2008，23（4）：419-425.

WEI X Z，HUANG Z Y，HONG F W.Large eddy simu?lation of flow field about marine propeller on unstruc?tured meshes［J］.Chinese Journal of Hydrodynamics，2008，23（4）：419-425.

Prediction Method for Open-Water Performance of Propeller Based on OpenFOAM

ZHENG Chao-sheng ZHANG Zhi-rong

China Ship Scientific Research Center，Wuxi 214082，China

In order to further investigate the viscous flow around the propeller，open-water performance prediction of DTMB P4119 propeller was performed using RANS flow solver in OpenFOAM software.And the grid dependence and effects of discretization schemes were investigated.The pressure distributions on the blade sections at different radii and velocity distributions ahead and behind the propeller were also ana?lyzed.To ensure the stability and accuracy of convective term discretization，the approaches are as fol?lows：（a）adopting second order NVD Gamma blended differencing scheme；（b）using k－wSST model as the turbulence model；（c）applying the SIMPLE algorithm for pressure-velocity coupling.The predicted results agree well with the experimental data.

viscous flow；propeller；openFOAM；RANS；open water performance；turbulence model

U661.33

A

1673－3185（2012）03－30－06

10.3969/j.issn.1673－3185.2012.03.006

2011－11－14

鄭巢生（1987－），男，碩士研究生。研究方向：計算流體力學(xué)。E?mail：zcszcs2005＠163.com

張志榮（1966－），男，博士，研究員。研究方向：計算流體力學(xué)。E?mail：zhangzr8＠public.wx.js.cn

張志榮。

［責(zé)任編輯：李勇群］