孔為平 王金寶 馮學梅

（中國船舶及海洋工程設計研究院 上海200011）

0 引 言

螺旋槳敞水性能是船舶快速性能的一個基本組成部分，與船舶阻力性能預報一樣，其準確的定量預報是迄今為止船舶計算流體力學（CFD）發(fā)展比較受關注的一個方面。

螺旋槳敞水性能、尾流場的數(shù)值計算方法的開發(fā)與驗證，構成了目前螺旋槳水動力性能研究內(nèi)容的主要部分，受到歷屆ITTC、歐盟框架計劃項目（FP5-LEADING EDGE項目、FP6-VIRTUE項目）以及許多計算流體力學工作組（CFD Workshop）的持續(xù)關注。

本部分計算對象為E779A型螺旋槳（E779A槳）和30.8萬噸超級油輪螺旋槳（30.8萬噸VLCC槳）,進行了敞水性能計算,并考察了網(wǎng)格、計算域、壓力離散格式等對敞水性能預報結果的影響。其中，E779A槳為意大利船模水池中心INSEAN的標槳，其敞水性能測量數(shù)據(jù)也從該水池獲得。INSEAN圍繞E779A槳敞水性能和空化特性進行了一系列完整細致的試驗研究，歐盟第六框架計劃（EU-FP6，即6th Framework Program of the European Union）的項目 VIRTUE（The Virtual Tank Utility in Europe）下的第四工作包 “數(shù)值空泡水池”（WP4，The Numerical Cavitation Tank）也選擇E779A槳作為研究對象，組織多個單位，采用多種求解器，對處于均勻來流中和船后尾流場中螺旋槳的空泡現(xiàn)象進行計算比較，并分別在2007年和2008年針對該槳組織了兩次相關的CFD Workshop［1］。 30.8萬噸VLCC螺旋槳的試驗數(shù)據(jù)則由中國船舶及海洋工程設計研究院水池提供。

馮學梅等曾針對E779A及PPTC槳，成功地進行了槳的敞水、尾流、空泡計算［2-3］。計算結果與試驗結果都相當一致，成為與會19份計算結果中最為突出者。在此基礎上，本文對其中一些數(shù)值計算作了分析和調(diào)整，以期使相關的算法更為合理，并形成更趨成熟、快捷的流程。

1 E779A螺旋槳敞水性能計算

1.1 計算對象和工況

歐盟框架計劃項目將E779A槳的測量數(shù)據(jù)資料視為機密，但隨時日推延，E779A槳水動力性能數(shù)值和水聲測量數(shù)據(jù)已有部分公開，并被廣泛用作基準校驗的對象，進行數(shù)值計算的比較研究。該槳帶有縱斜，為固定螺距的右旋四葉槳。槳模及其主參數(shù)見圖1和表1。

圖1 E779A槳模

表1 E779A槳模主要參數(shù)

表2 E779A槳敞水性能計算工況

1.2 計算方案

鑒于來流均勻和螺旋槳幾何上的周期性，為進行敞水計算，只需取單個槳葉所在的單通道（對于4葉槳而言，周向跨90°）作為計算域即可進行計算分析，由此可節(jié)省計算時間,提高計算效率。計算域的進口面、出口面均為90°的扇面，與槳軸垂直，彼此沿周向有所錯位；以槳轂和軸為內(nèi)邊界面；外邊界面取在圓柱體表面上，其直徑為螺旋槳直徑的數(shù)倍；周向側(cè)面由進、出口面附近的子午面和連貫其間的螺紋面組成。本部分敞水計算中，取進口在上游2.11D處，出口在下游3.08D處，外邊界所在圓柱面半徑為1.54D，槳軸為從進口一直延伸至出口的圓柱，其直徑為0.123 5D。來流方向與z軸方向相反，具體參照圖2。

圖2 計算區(qū)域及坐標

這里的螺旋槳敞水計算中，整個計算區(qū)域處于一個以轉(zhuǎn)速n繞槳軸作旋轉(zhuǎn)運動的參考坐標系中，槳葉相對此參考坐標系靜止不動。因此可選用FLUENT軟件提供的運動參考坐標系模型（即MRF模型）［4］。

采用基于壓力耦合的粘性求解器。壓力與速度的耦合引用SIMPLE方法，湍流模型為SSTκ-ω模型。對壓力項，考察采用了標準（Standard）和二階差分（Second Order）格式離散。動量方程其余項、湍流模型方程湍流動能（Turbulence Kinetic Energy）和湍流耗散率（Specific Dissipation Rate）項均采用二階迎風格式（Second Order Upwind）作離散。

在進口邊界處設置為速度進口條件，給定均勻來流的各速度分量；出口邊界給定表壓為0，即與參考點靜壓相等；外邊界同樣設為速度進口邊界；而兩個周向側(cè)面置為旋轉(zhuǎn)周期性邊界。

1.3 計算結果及分析

以下討論網(wǎng)格選擇、計算域的考察、壓力標準離散格式和壓力二階離散格式對計算收斂情況影響的比較以及收斂過程、時間的考察等工作，最后給出敞水計算的相對偏差和相關曲線。因E779A槳的敞水性能測量數(shù)據(jù)從意大利INSEAN水池獲得，并未授權公開，故以下計算所得敞水性能數(shù)據(jù)與試驗結果的比較都以相對偏差給出。

1.3.1 網(wǎng)格選擇

為考察網(wǎng)格對敞水計算結果的影響，對E779A槳生成3套網(wǎng)格進行計算比較。在CAD軟件Gambit中生成約84萬個四面體非結構化網(wǎng)格（G84），在Gambit和TGrid中生成約103萬個和125萬個帶邊界層網(wǎng)格的混合網(wǎng)格（G103和G125），其中第一層邊界層網(wǎng)格離壁距離為0.000 022D。網(wǎng)格特征信息見表3；網(wǎng)格G84、網(wǎng)格G103與網(wǎng)格G125之比較見圖3。

表3 網(wǎng)格特征

圖3 網(wǎng)格比較

針對上述3套網(wǎng)格，選擇進速系數(shù)J=0.397進行 計算考察（見表4），其中壓力項離散采用二階格式。

表4 J=0.397時，不同網(wǎng)格計算結果與試驗值相對偏差比較

表4中：

Δ表示相對偏差；

對葉面壓力分布進行了考察，可以看出：壓力量級達到1e4時，導邊及外緣壓力梯度非常大，參見圖4。所以，為更好捕捉真實的流動特性，槳葉邊緣應加密處理。相比較而言，從網(wǎng)格分布和邊界層布置來看，G125都更適應該敞水計算。

圖4 E7779A槳葉面壓力分布

由表4可以看出：G84沒用邊界層，Kq相差較大；G103計算出的Kt、Kq相對偏差都超過了 3%；G125網(wǎng)格的葉面輪廓進行了加密，計算結果都在2%以內(nèi)。下面將對G125在其他進速系數(shù)條件下進行詳細考察。

1.3.2 計算域影響

前面提到本次敞水計算區(qū)域。馮學梅等人針對該E779A槳的敞水計算在計算域上采用以下方法：取進口在上游2D處，出口在下游6D處，外邊界所在圓柱面半徑為3D，槳軸為從進口一直延伸至出口的圓柱，其直徑為0.123 5D。計算域相應位置明顯大于本計算域，為考察計算的可信度，在該計算域中提取出本計算域邊界相應位置上的壓力分布等相關量。

由圖5和圖6可以看出，在z=-0.7 m處壓力速度等值線都幾乎是周向均勻分布，并且沿徑向只在靠近壁面的小范圍內(nèi)變化，壓力量級基本為1e2，與槳葉上的壓力量級1e4相比，若取其為0，其影響應可忽略不計。速度分布也都在1.1 m/s左右，這與在J=0.397下所給定的速度入口條件基本吻合?？梢哉J為在z=-0.7 m處已經(jīng)達到了無窮遠條件。

圖5 z=-0.7 m時，槳葉面上的壓力分布（左：等值線；右：沿徑向分布）

圖6 z=-0.7 m時，槳葉面上的速度分布（左：等值線；右：沿徑向分布）

同樣，由圖7可以看出，在側(cè)邊界（r=0.35 m）上，無論是壓力分布還是速度分布，都和所給定的側(cè)面邊界條件基本吻合。因此，本文所取計算域能滿足計算要求。

圖7 r=0.35 m時，壓力和速度分布

1.3.3 壓力離散格式

對于壓力二階離散格式和標準離散格式的考察。標準離散格式是用計算出的單元體中心壓力和速度來插值到單元面上。二階壓力格式直接計算出單元面上壓力和速度值。二階格式一般適用于高速旋轉(zhuǎn)、高雷諾數(shù)、復雜幾何曲面的流動計算。而標準離散格式，則有較好的數(shù)值性能，比如收斂穩(wěn)定。

由圖8和圖9可以看出，當垂軸量度幅值相同時（一般在迭代500次左右），殘差、推力和扭矩就開始進入穩(wěn)定狀態(tài)；迭代1 000次時的推力和扭矩則非常穩(wěn)定。壓力標準離散格式收斂更好，殘差量級都降到1e-3以下，而且推力、扭矩都收斂于一條直線。相比之下，壓力二階離散格式收斂不太理想，連續(xù)方程和Omega殘差都在1e-3以上，推力和扭矩呈波動式收斂，數(shù)值波動約為±0.68%。從數(shù)值上來看，壓力標準格式推力、扭矩絕對值略低于壓力二階格式；與試驗數(shù)據(jù)相比較，壓力標準格式更接近試驗值。綜合看來，對于本敞水計算雷諾數(shù)1.9e6，轉(zhuǎn)速11.79 rad/s，使用標準壓力滿足計算精度要求，并且能收到更好的收斂效果。

圖8 J=0.397時的壓力標準離散格式收斂曲線

圖9 J=0.397時的壓力二階離散格式收斂曲線

1.3.4 收斂時間

本文對該敞水計算的收斂時間進行了考察分析，選擇了幾個進速系數(shù)，給出了相關計算收斂歷史曲線如圖 10～14。

圖10 J=0.397時的收斂曲線

圖11 J=0.498時的收斂曲線

圖12 J=0.695時的收斂曲線

圖13 J=0.845時的收斂曲線

圖14 J=0.946時的收斂曲線

表5 收斂結果整體對比（500～1 200步或1 500步）

表6 相關迭代步和試驗值結果對比

綜合上述圖表可知，本敞水計算經(jīng)500步后，無論推力T還是轉(zhuǎn)矩Q都已漸趨穩(wěn)定，其后的變化均達到千分之一左右甚至更小的相對偏差。使用6核，500步計算時間大約為45 min，即最少45 min可以完成一個進速系數(shù)點計算。

1.3.5 結果比較和分析

由表7和圖15可見（其中cfd為計算值，efd為試驗值），除去進速系數(shù)最大的J=0.995和J=1.094以外，G125的推力系數(shù)、扭矩系數(shù)以及敞水效率的相對偏差基本保持在2.5%以內(nèi)。其中，邊界層網(wǎng)格G125的敞水性能曲線幾乎與試驗曲線重疊。

表7 敞水性能計算結果相對偏差

圖15 敞水曲線比較

由圖15可知，敞水效率在J=0.95時最大；而計算結果相對試驗結果偏差過大的情況，則出現(xiàn)在J=0.995和J=1.094。此處，因推力、扭矩本身很小，故無論計算還是試驗的難度均較大。

2 30.8萬噸VLCC槳敞水性能計算

30.8 萬噸VLCC槳為備用槳，該槳沒有縱斜，為可變螺距的右旋四葉槳（見圖16）。

圖16 30.8萬噸VLCC槳模

其主要參數(shù)為：

槳模直徑D=0.2 m，螺距比P0.7/D=0.667，盤面比Ae/Ao=0.54，轂徑比 d/D=0.146，縱斜角=0°，葉數(shù)=4，右旋向。

其敞水性能計算工況為：

進速系數(shù) J=0.0～0.8，轉(zhuǎn)速 n=25.947 r/s，水溫為10.5℃，水運動粘性系數(shù)v=1.295 3×10-6m2/s，水密度 ρ=1 000.6 kg/m3。

根據(jù)E779A槳的計算經(jīng)歷，對30.8萬噸VLCC槳的計算不再考察網(wǎng)格的影響，而是直接采用帶邊界層網(wǎng)格的網(wǎng)格布置。設置葉片上最小網(wǎng)格尺寸為0.001D，第一層邊界層網(wǎng)格離壁大小為0.000 02D，在Gambit和TGrid軟件中生成約80萬個帶邊界層的四面體非結構化網(wǎng)格。計算方案與邊界條件的設置與E779A槳相同。

由表8和圖17可以看出計算所得敞水性能與試驗結果的比較?？梢姡ㄋ阅芮€與相應的試驗曲線相當吻合，當進速系數(shù)J=0.3～0.6時，無論推力、扭矩還是效率偏差都在1%以內(nèi)。再次驗證了這套敞水性能數(shù)值模擬方案的可行性和可靠性。

表8 30.8萬噸VLCC槳敞水性能結果比較

圖17 30.8萬噸VLCC槳敞水性能曲線

3 結 論

在敞水計算中，為得到精確結果，除了邊界條件和計算域合理外，網(wǎng)格布置也很關鍵，應該根據(jù)流動特性，合理密布關鍵區(qū)域尤其是槳葉上的網(wǎng)格，以及采用邊界層網(wǎng)格技術，而標準壓力插值格式的運用，也使計算結果精確且效率更高。

4 致 謝

本文在寫作過程中，得到了蔡榮泉研究員的悉心指導，也得益于于海、吳瓊等同事的熱心幫助，在此表示衷心感謝。

［1］MARNET-CFD. ［2012-03-10］.Final report and review of the state-of the-art in the application of CFD in the Maritime and Offshore Industries［EB/OL］.https：//pronet.wsatkins.co.uk/marnet/.

［2］馮學梅.水面艦船水動力快速性能數(shù)值研究［D］.上海：上海交通大學，2001.

［3］蔡榮泉，馮學梅.淺談艦船計算流體力學（CFD）實用化［J］.船舶，2012，23（2）：75-84.

［4］Fluent 6.1 Tutorial Guide［Z］.2003.