郭海鵬, 鄒早建, 李廣年

(1. 寧波大學(xué) 海運學(xué)院, 浙江 寧波 315211; 2. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院, 上海 200240)

船舶的緊急制動能力與船舶的航行安全性直接相關(guān)，而緊急制動能力的優(yōu)劣很大程度上取決于緊急倒車工況下螺旋槳的水動力特性.相對于目前研究較多的螺旋槳正車前進(jìn)工況，螺旋槳緊急倒車工況下的繞流場更為復(fù)雜，具有明顯的非定常及非線性特性，這導(dǎo)致作用在螺旋槳上的水動力載荷同樣呈現(xiàn)出復(fù)雜的變化規(guī)律.因此，從提高船舶航行安全性的角度出發(fā)，對螺旋槳緊急倒車工況下的水動力特性及繞流機(jī)理的深入研究能夠為全面評估船舶的緊急制動能力提供依據(jù).

針對螺旋槳緊急倒車工況，早期開展的相關(guān)研究主要基于實驗流體動力學(xué)(EFD)技術(shù)[1-2].近年來，發(fā)展迅速的計算流體動力學(xué)(CFD)技術(shù)成為開展螺旋槳水動力性能研究的有力工具.Chen等[3]基于自主開發(fā)求解器對DTMB4381螺旋槳的四象限工況，即緊急正車、正車前進(jìn)、緊急倒車及倒車后退進(jìn)行了數(shù)值模擬，分析了不同工況下的葉片壓力、流線分布及尾流場特征；李理等[4]基于Fluent對AU5-65型螺旋槳的四象限水動力性能進(jìn)行了數(shù)值預(yù)報，分析了推力系數(shù)的變化規(guī)律及相應(yīng)的流場特征；Verma等[5]和Kumar等[6]基于自主開發(fā)求解器對DTMB4381螺旋槳的緊急倒車工況進(jìn)行了數(shù)值模擬，探討了上游船體對螺旋槳的影響、螺旋槳側(cè)向力與葉片附近流動分離間的聯(lián)系；王貴彪等[7]基于Fluent對某導(dǎo)管槳的四象限工況進(jìn)行了數(shù)值模擬，分析了不同工況下的推力、轉(zhuǎn)矩及周圍流場特征；張文璨等[8]基于Fluent對MARIN B系列螺旋槳的四象限工況水動力性能進(jìn)行了數(shù)值預(yù)報，并將數(shù)值結(jié)果與試驗圖譜進(jìn)行了對比驗證；Pergande等[9]基于ANSYS CFX對某調(diào)距槳的緊急倒車工況進(jìn)行了數(shù)值模擬，探討了該工況下的流動分離現(xiàn)象；Pontarelli等[10]基于ANSYS REX對DTMB4381和MARIN 7371R螺旋槳的緊急倒車工況進(jìn)行了數(shù)值模擬，對不同進(jìn)速條件下的流場特征進(jìn)行了分析；Chen等[11]基于Fluent對DTMB4381螺旋槳緊急正車、緊急倒車工況進(jìn)行了數(shù)值模擬，分析了不同工況下的槳葉壓力及流線分布等流場特征；王超等[12-13]基于STAR-CCM+對E1619螺旋槳和SUBOFF潛艇模型的緊急倒車工況進(jìn)行了數(shù)值模擬，分析了槳葉誘導(dǎo)環(huán)狀渦的產(chǎn)生機(jī)理、發(fā)展過程以及上游艇體產(chǎn)生的影響；樊翔等[14]基于STAR-CCM+對KCS船模的倒車工況進(jìn)行了數(shù)值模擬，分析了該工況下的螺旋槳載荷變化特性；Wang等[15]基于自主開發(fā)的求解器naoe-FOAM-SJTU對KCS船模的緊急停船運動進(jìn)行了數(shù)值模擬，分析了該過程中船體與螺旋槳之間的相互干擾特性.

綜上，針對螺旋槳緊急倒車工況的相關(guān)數(shù)值研究大多借助商業(yè)CFD軟件實現(xiàn)，這類軟件存在價格昂貴、代碼不可見等局限性.相較而言，開源CFD平臺OpenFOAM以其開源免費、支持高度自定義等優(yōu)勢，得到了越來越多研究人員的關(guān)注.本文基于OpenFOAM對5葉螺旋槳DTMB4381的正車前進(jìn)和緊急倒車工況進(jìn)行數(shù)值模擬，通過與國際公開的模型試驗結(jié)果進(jìn)行對比，考察所采用數(shù)值方法的有效性，并基于數(shù)值模擬結(jié)果探討了螺旋槳緊急倒車工況下局部繞流場特征隨進(jìn)速的變化規(guī)律及其與螺旋槳整體水動力特性之間的聯(lián)系.

1 求解方法

1.1 控制方程

所考慮的螺旋槳繞流為黏性、不可壓縮流體流動，其控制方程為雷諾平均的連續(xù)性方程和動量方程，即雷諾平均納維斯托克斯(RANS)方程：

(1)

(2)

式中：ρ為流體密度；U為平均速度；p為平均壓力；μeff=ρ(ν+νt)為有效動力黏度，其中ν和νt為運動黏度和湍流渦黏度，后者通過湍流模型確定.

1.2 湍流模型

湍流渦黏度采用Menter等[16]提出的兩方程SSTk-ω模型進(jìn)行求解，其中k和ω為湍流動能和比耗散率.該模型通過剪切應(yīng)力輸運公式將k-ω模型與標(biāo)準(zhǔn)k-ε模型結(jié)合，在近壁面采用k-ω模型，而在遠(yuǎn)場采用標(biāo)準(zhǔn)k-ε模型，充分發(fā)揮了兩者的優(yōu)勢，是目前船舶與海洋工程領(lǐng)域應(yīng)用較為廣泛的湍流模型.在OpenFOAM中，SSTk-ω模型的輸運方程[17]如下：

(3)

(4)

(5)

(6)

式中：S為應(yīng)變率的不變量；β*、αk、γ、β、αω及αω2為方程常量；F1為結(jié)合k-ω模型和k-ε模型的混合函數(shù).變量的具體定義及常數(shù)的取值見文獻(xiàn)[16-17].

1.3 速度壓力耦合求解

控制方程求解過程中的速度與壓力的耦合通過OpenFOAM中的PIMPLE算法實現(xiàn).該方法將Patankar等[18]在1972年提出的壓力耦合方程組的半隱式方法(SIMPLE算法)與Issa[19]在1986年提出的PISO算法相結(jié)合，可以實現(xiàn)瞬態(tài)問題的大尺度時間步求解以提高計算效率.

1.4 物體運動處理

螺旋槳旋轉(zhuǎn)運動模擬采用OpenFOAM中的多重參考系法(MRF)及任意網(wǎng)格界面法(AMI)來實現(xiàn)[20]：MRF方法通過在不同區(qū)域內(nèi)建立不同的參考系與控制方程將螺旋槳動態(tài)旋轉(zhuǎn)轉(zhuǎn)化為靜態(tài)問題進(jìn)行求解；AMI方法通過在不同區(qū)域交界面插值進(jìn)行通量傳遞，從而實現(xiàn)螺旋槳動態(tài)旋轉(zhuǎn)運動的直接模擬.

2 數(shù)值計算

2.1 研究對象及工況

以5葉螺旋槳DTMB4381模型為研究對象，其直徑為0.305 m，螺距比為1.2，側(cè)斜角為0°，盤面比為0.725，轂徑比為0.2，設(shè)計進(jìn)速系數(shù)為0.889，螺旋槳為右旋，其幾何模型如圖1所示.首先對螺旋槳正車前進(jìn)工況進(jìn)行數(shù)值模擬，然后進(jìn)行螺旋槳緊急倒車工況數(shù)值模擬，倒車工況的進(jìn)速方向與正車前進(jìn)工況相同，但螺旋槳反向旋轉(zhuǎn).詳細(xì)模擬工況如表1所示.其中，J為進(jìn)速系數(shù)，UA為進(jìn)速，n為螺旋槳轉(zhuǎn)速.

圖1 DTMB4381螺旋槳的后視圖及側(cè)視圖Fig.1 Back view and side view of DTMB4381 propeller

表1 數(shù)值模擬工況Tab.1 Numerical simulation condition

計算結(jié)果包括不同進(jìn)速系數(shù)下的推力系數(shù)和轉(zhuǎn)矩系數(shù)，定義如下：

(7)

(8)

(9)

式中：D為螺旋槳直徑；T為螺旋槳推力；Q為螺旋槳轉(zhuǎn)矩.

2.2 計算域及網(wǎng)格劃分

計算域及邊界條件設(shè)置如圖2所示.根據(jù)螺旋槳流場特性選取圓柱形計算域，其直徑為6D，長度為8D，其中槳盤面前、后分別為3D、5D，同時選取包含螺旋槳的圓柱體作為旋轉(zhuǎn)域，其直徑為1.5D，長度為2D，其中槳盤面前、后分別為1D.入口邊界設(shè)為fixed value邊界條件；出口邊界設(shè)為inlet outlet邊界條件，根據(jù)速度方向在fixed value與zeroGradient之間進(jìn)行切換；外側(cè)邊界設(shè)為no slip邊界條件；螺旋槳設(shè)為moving wall velocity邊界條件；旋轉(zhuǎn)域與靜止域的交界面設(shè)為cyclic AMI邊界條件，可以實現(xiàn)不同區(qū)域間的流場信息傳遞.

圖2 計算域及邊界條件設(shè)置Fig.2 Setup of computational domain and boundary condition

計算域采用OpenFOAM自帶的前處理工具進(jìn)行空間離散，網(wǎng)格劃分如圖3所示.首先，采用blockMesh工具將計算域離散為均勻的六面體背景網(wǎng)格，然后采用snappyHexMesh工具進(jìn)行螺旋槳附近區(qū)域的網(wǎng)格細(xì)化、表面對齊等操作.考慮到螺旋槳繞流場的特性，在螺旋槳附近及下游區(qū)域進(jìn)行體加密，以便捕捉螺旋槳尾流場特征.為了更好地求解螺旋槳近壁面流動，在螺旋槳壁面附近生成4層棱柱形網(wǎng)格，無因次壁面距離y+值控制在30左右.

圖3 計算域網(wǎng)格劃分Fig.3 Grid generation of computational domain

2.3 迭代計算設(shè)置

為了提高數(shù)值計算效率及穩(wěn)定性，在計算過程中首先應(yīng)用simpleFoam求解器開展定常計算，螺旋槳旋轉(zhuǎn)采用MRF方法模擬，以此獲得穩(wěn)定的初始流場.然后，以獲得的初始流場結(jié)果為初始值，應(yīng)用pimpleFoam求解器開展非定常計算，螺旋槳旋轉(zhuǎn)采用AMI方法模擬.在非定常計算過程中，時間步長設(shè)為1.0×10-4s，模擬物理時間為5.0 s，對應(yīng)螺旋槳旋轉(zhuǎn)50圈所需時間.

2.4 網(wǎng)格收斂性分析

圖4 槳葉網(wǎng)格劃分Fig.4 Grid generation on blade

網(wǎng)格收斂性分析結(jié)果見表2，其中KT和KQ的計算值通過對螺旋槳旋轉(zhuǎn)最后5圈的數(shù)據(jù)進(jìn)行平均獲得.從表中可以看出，KT和KQ的網(wǎng)格收斂率RG為負(fù)值且絕對值小于1，因而表現(xiàn)為振蕩收斂.同時，網(wǎng)格不確定性無量綱值UG%S1<1，表明當(dāng)前計算結(jié)果受網(wǎng)格尺度的影響較小.此外，不同尺度網(wǎng)格計算對比結(jié)果表明，細(xì)網(wǎng)格計算結(jié)果與中網(wǎng)格計算結(jié)果相差很小.因此，權(quán)衡數(shù)值計算的精度和效率，在后續(xù)數(shù)值研究中均采用中網(wǎng)格進(jìn)行計算.

表2 網(wǎng)格收斂性分析結(jié)果Tab.2 Results of grid convergence study

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

3.1 水動力載荷

螺旋槳正車前進(jìn)工況下的推力、轉(zhuǎn)矩以及效率的計算結(jié)果如圖5所示.通過與國際上公開的模型試驗結(jié)果[22]進(jìn)行對比發(fā)現(xiàn)，本文計算結(jié)果與模型試驗值吻合良好，推力誤差值基本在3%以內(nèi)，在高進(jìn)速(J=1.1)工況下誤差值最大達(dá)到14%；轉(zhuǎn)矩誤差整體相對略高，在低進(jìn)速(J=0.1)工況下誤差值最大約為7%；隨著進(jìn)速增大，誤差值基本在5%以內(nèi).總體而言，數(shù)值計算精度與常規(guī)的商業(yè)CFD計算軟件相當(dāng)[23]，驗證了所采用的計算平臺及數(shù)值方法在預(yù)報螺旋槳正車前進(jìn)工況下的水動力性能方面的有效性.

圖5 螺旋槳正車前進(jìn)工況模擬結(jié)果Fig.5 Simulation results of propeller ahead condition

螺旋槳緊急倒車工況下推力、轉(zhuǎn)矩的計算結(jié)果如圖6所示.通過與國際上公開的模型試驗結(jié)果[22]進(jìn)行對比發(fā)現(xiàn)，數(shù)值計算結(jié)果與模型試驗結(jié)果相比存在一定偏差，且隨著進(jìn)速系數(shù)的改變發(fā)生顯著變化.在低進(jìn)速(J=-0.5)工況下，推力及轉(zhuǎn)矩的計算值偏高，誤差值均在15%左右；在中進(jìn)速(J=-0.7)工況下，推力及轉(zhuǎn)矩的計算值略微偏低，誤差值分別為8%和2%；在高進(jìn)速(J=-0.9)工況下，推力及轉(zhuǎn)矩的計算值明顯偏低，誤差值分別達(dá)到20%和15%.根據(jù)前文網(wǎng)格收斂性分析結(jié)果可知，當(dāng)前螺旋槳緊急倒車工況下的水動力性能預(yù)報精度受網(wǎng)格尺度影響較小，而應(yīng)該與所采用的湍流建模方法難以準(zhǔn)確捕捉該工況下存在的大規(guī)模流動分離現(xiàn)象有關(guān).總體而言，盡管本文所采用的數(shù)值方法對螺旋槳緊急倒車工況下的水動力預(yù)報精度相對有限，但數(shù)值計算結(jié)果仍呈現(xiàn)出與模型試驗結(jié)果相一致的總體趨勢.

圖6 螺旋槳緊急倒車工況模擬結(jié)果Fig.6 Simulation results of propeller crashback condition

3.2 繞流場特征

3.2.1泄出渦系結(jié)構(gòu) 不同進(jìn)速系數(shù)下數(shù)值模擬得到的緊急倒車工況下螺旋槳繞流場中的渦系結(jié)構(gòu)如圖7所示.渦系結(jié)構(gòu)采用渦量等值面Q=1 000 表示，同時采用無因次軸向速度Ux/U0著色，其中U0為遠(yuǎn)場入流速度.結(jié)果表明，緊急倒車工況下螺旋槳梢部產(chǎn)生的泄出渦十分明顯，并呈現(xiàn)出與正車前進(jìn)工況截然不同的結(jié)構(gòu).由于緊急倒車工況下螺旋槳反向旋轉(zhuǎn)，螺旋槳產(chǎn)生的排出流與來流方向相反.在來流影響下，槳葉梢部泄出渦有沿螺旋槳徑向發(fā)展的趨勢，同時，螺旋槳葉背根部也產(chǎn)生了明顯的泄出渦.由渦量等值面上的速度分布可知，在靠近葉梢部位，螺旋槳誘導(dǎo)速度的軸向分量與來流方向相反，在葉梢外部區(qū)域，來流的影響逐漸占據(jù)主導(dǎo)地位，呈現(xiàn)出與來流相一致的流向.

通過對比不同進(jìn)速下的渦量等值面圖可以發(fā)現(xiàn)，渦系結(jié)構(gòu)的尺度和位置隨入流速度的改變而發(fā)生變化.在低進(jìn)速(J=-0.5)工況下，螺旋槳梢部泄出渦沿周向均勻分布且整體向螺旋槳上游發(fā)展.隨著進(jìn)速的增加，螺旋槳葉梢部泄出渦逐步增強(qiáng)并呈現(xiàn)出周向非均勻特性，而且逐漸呈現(xiàn)出向螺旋槳下游發(fā)展的趨勢.與此同時，葉背根部的泄出渦也愈發(fā)顯著，葉面根部也有新的泄出渦生成.此外，由于來流速度影響的增強(qiáng)，螺旋槳梢渦上的軸向速度分布也發(fā)生了明顯變化，來流的影響逐漸在螺旋槳繞流場中占據(jù)主導(dǎo)地位.

3.2.2剖面流線分布 不同進(jìn)速系數(shù)下數(shù)值模擬得到的緊急倒車工況下螺旋槳附近流場剖面的流線分布如圖8所示.此外，圖中還給出了剖面上的無因次軸向速度Ux/U0分布.剖面流線分布情況清晰地反映了螺旋槳緊急倒車工況下渦系結(jié)構(gòu)的產(chǎn)生機(jī)理及特性：在緊急倒車工況下，螺旋槳反轉(zhuǎn)所誘導(dǎo)產(chǎn)生的流動與來流方向相反，在螺旋槳上游流場產(chǎn)生擠壓、剪切，形成沿螺旋槳徑向發(fā)展的匯聚流，結(jié)果與圖7中的槳葉梢部泄出渦相對應(yīng)；同時，由于螺旋槳反轉(zhuǎn)所產(chǎn)生的流動影響，來流在螺旋槳下游發(fā)生分離而形成分離流.在匯聚流與分離流的共同作用下，槳盤面后方形成了大尺度環(huán)狀渦結(jié)構(gòu).

圖7 緊急倒車工況下Q=1 000等值面Fig.7 Iso-surface of Q=1 000 under crashback condition

通過對比不同進(jìn)速下的流線圖可知，由匯聚流和分離流形成的環(huán)狀渦結(jié)構(gòu)的尺度和位置隨入流速度的改變而發(fā)生變化.低進(jìn)速(J=-0.5)工況下螺旋槳誘導(dǎo)速度的影響占據(jù)主導(dǎo)地位，匯聚流和分離流的作用范圍較廣，形成的環(huán)狀渦結(jié)構(gòu)尺度也較大.隨著進(jìn)速的增加，來流的影響逐漸增大，匯聚流和分離流的作用區(qū)域逐漸減小，所形成的大尺度環(huán)狀渦逐漸收縮，并呈現(xiàn)出向螺旋槳下游方向移動的趨勢.

3.2.3葉片壓力分布 數(shù)值模擬得到的緊急倒車工況下的螺旋槳葉背及葉面壓力分布如圖9和圖10所示.Cp表示無因次化的壓力系數(shù)，計算公式為

(10)

式中：Δp為相對壓力.緊急倒車工況下由于螺旋槳反向旋轉(zhuǎn)，導(dǎo)邊與隨邊也相應(yīng)改變.螺旋槳葉背以正壓為主，表現(xiàn)為壓力面，最大正壓出現(xiàn)在導(dǎo)邊附近區(qū)域；相反地，螺旋槳葉面以負(fù)壓為主，表現(xiàn)為吸力面，最大負(fù)壓出現(xiàn)在導(dǎo)邊及葉梢附近區(qū)域.由于壓力面與吸力面的壓差作用，螺旋槳緊急倒車工況下產(chǎn)生反向的推力.

不同進(jìn)速下的葉片壓力云圖對比結(jié)果表明，葉片壓力分布隨入流速度的改變而發(fā)生變化.在低進(jìn)速(J=-0.5)工況下，葉片上的高壓區(qū)及低壓區(qū)主要位于導(dǎo)邊附近區(qū)域，且葉片壓力呈現(xiàn)出較好的周向均勻性.隨著進(jìn)速的增大，壓力面的高壓區(qū)及吸力面的低壓區(qū)均有所增強(qiáng)，特別是壓力面的高壓區(qū)逐漸發(fā)展到整個葉背大部分區(qū)域.值得注意的是，在壓力面靠近導(dǎo)邊的葉根部分區(qū)域呈現(xiàn)出相對較低的壓力，這與高進(jìn)速工況下的葉背根部位置的泄出渦有關(guān).類似地，吸力面靠近導(dǎo)邊的葉根部位形成了新的低壓區(qū)，同樣與葉面根部的泄出渦有關(guān).此外，隨著進(jìn)速的增加，葉片上的壓力分布還呈現(xiàn)出一定的周向非均勻性，這與槳葉泄出渦的結(jié)構(gòu)特征相對應(yīng).

圖9 緊急倒車工況下螺旋槳葉背壓力分布Fig.9 Pressure distribution on the back of propeller under crashback condition

圖10 緊急倒車工況下螺旋槳葉面壓力分布Fig.10 Pressure distribution on the face of propeller under crashback condition

4 結(jié)論

基于開源CFD平臺OpenFOAM對螺旋槳DTMB4381模型的正車前進(jìn)及緊急倒車工況進(jìn)行數(shù)值模擬，通過與模型試驗數(shù)據(jù)的對比以及對數(shù)值計算結(jié)果的分析，可以得到以下結(jié)論：

(1) 對于螺旋槳正車前進(jìn)工況，本文數(shù)值模擬結(jié)果與模型試驗數(shù)據(jù)吻合良好，除高進(jìn)速(J=1.1)工況外，推力和轉(zhuǎn)矩的誤差均在5%以內(nèi)；對于螺旋槳緊急倒車工況，數(shù)值模擬結(jié)果與模型試驗數(shù)據(jù)相差較大，推力和轉(zhuǎn)矩在低進(jìn)速(J=-0.5)和中進(jìn)速(J=-0.7)工況下誤差在10%左右，在高進(jìn)速(J=-0.9)工況下誤差達(dá)到20%.

(2) 緊急倒車工況下，螺旋槳梢部泄出渦在來流影響下呈現(xiàn)沿螺旋槳徑向發(fā)展的趨勢，同時葉片根部也有泄出渦產(chǎn)生.隨著來流速度的增大，槳葉泄出渦有所增強(qiáng)并向螺旋槳下游發(fā)展，此外還呈現(xiàn)出一定的周向非均勻性.

(3) 緊急倒車工況下，由于螺旋槳的誘導(dǎo)速度與來流相互作用，形成了螺旋槳上游的匯聚流及下游的分離流，兩者共同作用在槳盤面后方發(fā)展成為大尺度環(huán)狀渦.該環(huán)狀渦隨著來流速度的增加而逐漸收縮，并呈現(xiàn)出向螺旋槳下游方向移動的趨勢.

(4) 緊急倒車工況下，螺旋槳葉背為壓力面而葉面為吸力面，壓力面的高壓區(qū)和吸力面的低壓區(qū)主要分布在導(dǎo)邊附近.隨著來流速度的增加，槳葉吸力面上的低壓區(qū)及壓力面上的高壓區(qū)均有所增強(qiáng)，同時壓力分布呈現(xiàn)出與槳葉泄出渦相類似的周向非均勻性.

由以上研究結(jié)果可知，螺旋槳緊急倒車工況下的流場特征較為復(fù)雜，涉及大規(guī)模流動分離及回流現(xiàn)象，而本文所采用的RANS方法在該方面求解精度相對有限.因此，需要針對這一問題開展更為深入的研究，將采用具有更高精度的數(shù)值方法如大渦模擬(LES)方法、分離渦模擬(DES)方法等進(jìn)行湍流模擬，從而更準(zhǔn)確地揭示螺旋槳緊急倒車工況下的水動力特性及相應(yīng)的流動機(jī)理.