李兆輝 胡 帆 吳 瓊,2 馮 毅,2 孫 群,2

（1. 中國(guó)船舶及海洋工程設(shè)計(jì)研究院 上海 200011； 2. 上海市船舶工程重點(diǎn)實(shí)驗(yàn)室 上海 200011）

0 引 言

目前，船舶領(lǐng)域的CFD計(jì)算主要集中在模型尺度。由于存在尺度效應(yīng)，計(jì)算結(jié)果需要通過(guò)經(jīng)驗(yàn)統(tǒng)計(jì)公式換算到實(shí)船，這些公式并不具有普適性，而且難以獲得實(shí)船雷諾數(shù)下的流場(chǎng)，對(duì)探索新的提升船舶能效的水動(dòng)力解決方案形成了制約。而直接開(kāi)展實(shí)尺度CFD計(jì)算無(wú)疑是一個(gè)很好的技術(shù)途徑，但所采用的實(shí)尺度CFD計(jì)算方法必須經(jīng)過(guò)仔細(xì)和獨(dú)立的驗(yàn)證，才能在業(yè)內(nèi)獲得普遍接受和認(rèn)可。

國(guó)外：2006年，HANNINEN等分析了EFFORT項(xiàng)目的算例，驗(yàn)證了RANS 求解器用于實(shí)尺度數(shù)值模擬的可行性；2009年，BHUSHAN等基于URANS 求解器，證明了在考慮壁面粗糙度和壓力梯度影響下，兩點(diǎn)多層壁面函數(shù)模型尺度與實(shí)尺度數(shù)值模擬中的通用性；2011年，CASTRO等采用重疊網(wǎng)格研究了實(shí)尺度下螺旋槳和船體的相互作用，并分析了船模尺度和實(shí)尺度對(duì)船的邊界層和螺旋槳推進(jìn)性能的影響。

國(guó)內(nèi)：2015年，陳天福等采用結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格相結(jié)合的計(jì)算方法，進(jìn)行了實(shí)尺度導(dǎo)管螺旋槳的敞水性能數(shù)值模擬；2016年，尹崇宏等應(yīng)用改進(jìn)的DES 模型（IDDES），采用CFD 方法同時(shí)對(duì)模型尺度和實(shí)尺度的32萬(wàn)t VLCC 進(jìn)行阻力預(yù)報(bào)，并分別與RANS 的計(jì)算結(jié)果和試驗(yàn)結(jié)果進(jìn)行比較分析；2016年，李亮等采用RANS方法和VOF模型，開(kāi)展了考慮自由液面的實(shí)船自航性能的數(shù)值模擬，分析發(fā)現(xiàn)自航點(diǎn)轉(zhuǎn)速和伴流分?jǐn)?shù)尺度效應(yīng)明顯；2018年，陳騫等進(jìn)行了實(shí)尺度阻力模擬，并與水池試驗(yàn)結(jié)果的實(shí)尺度換算數(shù)值進(jìn)行對(duì)比，指出實(shí)尺度下伴流更加均勻。

本文采用體積力和滑移網(wǎng)格這2種螺旋槳模擬方式對(duì)9 400 TEU集裝箱船進(jìn)行實(shí)尺度自航性能數(shù)值模擬，分別考察湍流模型、網(wǎng)格劃分、最大非線性迭代次數(shù)、時(shí)間步長(zhǎng)、激勵(lì)盤(pán)參數(shù)及動(dòng)區(qū)域大小等因素對(duì)計(jì)算結(jié)果的影響，并與模型試驗(yàn)換算及實(shí)船試航結(jié)果對(duì)比，分析2種螺旋槳模擬方式的優(yōu)劣并驗(yàn)證實(shí)尺度自航性能數(shù)值模擬的可靠性。

1 數(shù)學(xué)模型

本文使用FINE/Marine商用軟件，它是NUMECA公司推出的船舶與海洋工程專業(yè)CFD軟件包包括網(wǎng)格生成器HEXPRESS、黏流求解器ISISCFD和后處理器CFView。FINE/Marine軟件基于RANS方法求解船體周?chē)鲌?chǎng)，使用流體體積函數(shù)（volume of fluid，VOF）方法捕捉自由液面，采用動(dòng)網(wǎng)格技術(shù)處理船舶姿態(tài)變化。

1.1 控制方程

FINE/Marine的求解器ISIS-CFD流動(dòng)控制方程為不可壓縮非定常雷諾平均N-S方程（RANS），采用有限體積法對(duì)輸運(yùn)方程進(jìn)行空間離散。在等溫條件下，不可壓縮黏性流體的質(zhì)量和動(dòng)量守恒方程如下：

1.2 湍流模型

在雷諾平均N-S方程中，由于雷諾應(yīng)力是未知的，因此方程的封閉需要定義湍流雷諾應(yīng)力。根據(jù)對(duì)雷諾應(yīng)力的?；绞讲煌?，湍流模型分為渦黏模型和雷諾應(yīng)力模型。渦黏模型基于BOUSSINESQ渦黏假設(shè)，引入湍流渦黏系數(shù)，建立雷諾應(yīng)力與平均應(yīng)變率的關(guān)系；雷諾應(yīng)力模型則直接利用雷諾應(yīng)力輸運(yùn)方程，對(duì)方程中相關(guān)項(xiàng)進(jìn)行模化使方程封閉。由于RANS方法采用時(shí)間平均忽略了湍流脈動(dòng)，湍流模型也不具備普適性；而大渦模擬（large eddy simulation，LES）對(duì)大尺度渦結(jié)構(gòu)直接計(jì)算，對(duì)小尺度渦構(gòu)造亞格子應(yīng)力模型進(jìn)行模化，從而獲得流動(dòng)的動(dòng)態(tài)特性。但是，LES對(duì)網(wǎng)格尺寸的要求限制了工程應(yīng)用，1997年，SPALART提出了一種混合RANS/LES的分離渦模型（detached eddy simulation，DES）。

本文分別使用渦黏模型、雷諾應(yīng)力模型和分離渦模型這三類(lèi)湍流模型中的（SST-Menter）模型、EASM模型和DES-SST模型。

1.3 VOF方法

ISIS-CFD求解器采用VOF方法捕捉自由液面，VOF方法把空氣和水作為兩相流體進(jìn)行計(jì)算，多相流的體積分?jǐn)?shù)守恒方程如下：

式中：c是第相體積分?jǐn)?shù)。當(dāng)空間中全部是第相流體時(shí)，c=1；而沒(méi)有第相流體時(shí)，c=0。顯然體積分?jǐn)?shù)介于0～1時(shí)，表示是多相混合流體。由于船舶CFD計(jì)算中通常只包含空氣和水兩相，因此將c=0.5處定義為交界面。

多相混合流體的物性參數(shù)（黏性和密度）是各組成相的物性參數(shù)（黏性μ和密度ρ），由本構(gòu)關(guān)系得到：

1.4 體積力方法

式中各參數(shù)表達(dá)式為：

式中：K和K分別是推力系數(shù)和轉(zhuǎn)矩系數(shù)；是進(jìn)速系數(shù)；是轉(zhuǎn)速，r/s；Ω是旋轉(zhuǎn)速度，rad/s；R是螺旋槳半徑，m；R是槳轂半徑，m；是參考速度，m/s；b是投影到-平面的平均弦長(zhǎng)，m（或激勵(lì)盤(pán)厚度）；Y和Z是螺旋槳中心坐標(biāo)。

2 計(jì)算模型

2.1 計(jì)算對(duì)象及工況

以9 400 TEU集裝箱船為計(jì)算對(duì)象，船體主尺度和螺旋槳幾何參數(shù)分別如表1和表2所示。

表1 船體主尺度

表2 螺旋槳幾何參數(shù)

該螺旋槳是與設(shè)計(jì)槳敞水性征相近的備用槳，本文數(shù)值計(jì)算和模型試驗(yàn)采用的是備用槳，而實(shí)船試航采用的是設(shè)計(jì)槳。計(jì)算時(shí)除船體外，包含的附體有艏側(cè)推孔及舵（含假舵及舵球）。該集裝箱船及其附體的幾何外形如圖1所示，笛卡爾坐標(biāo)系下，以尾垂線和基平面交點(diǎn)為原點(diǎn)。軸沿船長(zhǎng)方向，正向指向船首；軸沿船寬方向，正向指向船的左舷（從后往前看）；軸沿船深方向，正向垂直向上。

圖1 幾何模型及坐標(biāo)系

對(duì)幾何模型進(jìn)行實(shí)船尺度的數(shù)值計(jì)算，計(jì)算工況如表3所示，雷諾數(shù)=2.693 6×10，弗勞德數(shù)=0.215。

表3 計(jì)算工況

2.2 計(jì)算域及邊界條件

由于螺旋槳的轉(zhuǎn)動(dòng)帶動(dòng)周?chē)黧w產(chǎn)生旋轉(zhuǎn)，導(dǎo)致船體周?chē)鲃?dòng)不對(duì)稱，因此實(shí)尺度自航性能計(jì)算選取全船進(jìn)行計(jì)算。計(jì)算域如圖2所示。

圖2 計(jì)算域

計(jì)算域各邊界面距坐標(biāo)原點(diǎn)的距離如表4所示，L為垂線間長(zhǎng)。計(jì)算域的長(zhǎng)度為5L，寬度為4L，高度為2L。

計(jì)算域邊界條件設(shè)定見(jiàn)表4。船體表面除甲板外為固壁非滑移邊界，使用標(biāo)準(zhǔn)壁面函數(shù)，甲板為滑移邊界。

表4 計(jì)算域大小與邊界條件

2.3 網(wǎng)格劃分

采用非結(jié)構(gòu)化六面體網(wǎng)格生成工具HEXPRESS進(jìn)行網(wǎng)格劃分。圖3為計(jì)算域全局網(wǎng)格、船體表面網(wǎng)格及邊界層網(wǎng)格分布。

圖3 計(jì)算域、船體表面及邊界層網(wǎng)格

實(shí)尺度下雷諾數(shù)較大，邊界層內(nèi)湍流劇烈，與模型尺度相比需要布置更多的網(wǎng)格層數(shù)。采用壁面函數(shù)處理近壁面一般要求第1層網(wǎng)格的無(wú)量綱化壁面距離滿足30＜＜200，但是在實(shí)尺度下滿足該要求將導(dǎo)致網(wǎng)格數(shù)激增。2020年，蘇玉民等分析了國(guó)內(nèi)外實(shí)尺度船舶快速性數(shù)值計(jì)算的，建議采用的范圍為30～500。第1層邊界層網(wǎng)格的可以按式（6）估算：

式中：為特征長(zhǎng)度，m，一般取=L；Δy為第1層邊界層網(wǎng)格厚度，m。

2.4 激勵(lì)盤(pán)設(shè)置

采用激勵(lì)盤(pán)即體積力的方法進(jìn)行實(shí)船尺度自航性能數(shù)值計(jì)算，施加的體積力由實(shí)槳的敞水性征值（如圖4所示）確定，本文中實(shí)槳的敞水性征值由敞水試驗(yàn)結(jié)果利用1978年ITTC推薦的修正公式換算得到。對(duì)于試驗(yàn)數(shù)據(jù)點(diǎn)之間的敞水性征值，ISIS-CFD求解器會(huì)自動(dòng)利用插值獲取。對(duì)于實(shí)船尺度數(shù)值計(jì)算，不需要考慮摩擦阻力修正。激勵(lì)盤(pán)設(shè)置在螺旋槳中心點(diǎn)，激勵(lì)盤(pán)厚度參考螺旋槳槳葉側(cè)投影輪廓確定，內(nèi)半徑為槳轂半徑，外半徑為螺旋槳半徑。

圖4 槳模經(jīng)ITTC修正到實(shí)槳的敞水性征曲線

FINE/Marine對(duì)于采用激勵(lì)盤(pán)來(lái)確定自航點(diǎn)有2種計(jì)算方式：一種是給定航速求解螺旋槳需要的轉(zhuǎn)速，另一種是給定螺旋槳的轉(zhuǎn)速求解可以達(dá)到的航速。后者航速的迭代需要不斷模擬船體加速的過(guò)程，因此計(jì)算相對(duì)耗時(shí)，而前者轉(zhuǎn)速的迭代只需要對(duì)敞水?dāng)?shù)據(jù)進(jìn)行插值。本文選取第1種計(jì)算方式，通過(guò)不斷改變轉(zhuǎn)速直至船體阻力與激勵(lì)盤(pán)推力平衡來(lái)確定實(shí)船自航點(diǎn)。最終設(shè)置的激勵(lì)盤(pán)如圖5所示，與螺旋槳旋轉(zhuǎn)方向保持一致即為右旋。

圖5 激勵(lì)盤(pán)示意圖

2.5 滑移網(wǎng)格設(shè)置

采用滑移網(wǎng)格進(jìn)行自航性能數(shù)值計(jì)算時(shí)，需在螺旋槳周?chē)O(shè)置動(dòng)區(qū)域，動(dòng)區(qū)域設(shè)置為圓柱體，如圖6所示。動(dòng)區(qū)域需包括整個(gè)螺旋槳幾何，其直徑厚度不宜過(guò)小，以避免動(dòng)區(qū)域外邊界到螺旋槳距離過(guò)小導(dǎo)致其附近網(wǎng)格質(zhì)量較差的問(wèn)題。動(dòng)區(qū)域與螺旋槳一同轉(zhuǎn)動(dòng)，動(dòng)區(qū)域外邊界與靜區(qū)域內(nèi)邊界產(chǎn)生相對(duì)旋轉(zhuǎn)，為實(shí)現(xiàn)計(jì)算過(guò)程中靜區(qū)域與動(dòng)區(qū)域在交界面進(jìn)行數(shù)據(jù)交換，將靜區(qū)域與動(dòng)區(qū)域交界面設(shè)置為完全非匹配邊界條件（full non-matching boundary）。螺旋槳及尾軸表面設(shè)置為固壁非滑移邊界，使用標(biāo)準(zhǔn)壁面函數(shù)。

螺旋槳網(wǎng)格劃分時(shí)，對(duì)槳葉導(dǎo)邊、隨邊和葉根處進(jìn)行加密，如圖6所示。

圖6 螺旋槳表面網(wǎng)格及動(dòng)區(qū)域

采用真實(shí)槳的自航數(shù)值模擬確定實(shí)船自航點(diǎn)，通常是在給定航速后計(jì)算多個(gè)螺旋槳轉(zhuǎn)速，根據(jù)船體阻力和螺旋槳推力隨螺旋槳轉(zhuǎn)速變化曲線確定2個(gè)曲線交點(diǎn)即為該航速的自航點(diǎn)。該方法由于需要計(jì)算多個(gè)轉(zhuǎn)速，相對(duì)耗時(shí)并且存在插值誤差。FINE/Marine開(kāi)發(fā)了一種動(dòng)態(tài)控制器，在計(jì)算過(guò)程中自動(dòng)調(diào)整轉(zhuǎn)速，使作用在船體上的合力為0，從而在一次計(jì)算中確定自航點(diǎn)，但仍然比固定航速和轉(zhuǎn)速時(shí)耗時(shí)更多。

因此，本文在開(kāi)展基于滑移網(wǎng)格的實(shí)尺度自航性能數(shù)值計(jì)算方法研究時(shí)，首先固定航速和轉(zhuǎn)速，以收到功率P為指標(biāo)評(píng)判計(jì)算精度；然后根據(jù)數(shù)值計(jì)算方法研究結(jié)果，選取合適網(wǎng)格劃分方案和數(shù)值參數(shù)，開(kāi)啟控制器，在給定航速的情況下預(yù)報(bào)實(shí)船自航點(diǎn)；最后，與體積力方法預(yù)報(bào)的實(shí)船自航點(diǎn)進(jìn)行分析對(duì)比。

3 基于體積力的實(shí)尺度自航性能數(shù)值計(jì)算方法研究

3.1 網(wǎng)格依賴性研究

針對(duì)邊界層網(wǎng)格的不同劃分方案進(jìn)行計(jì)算，分析邊界層網(wǎng)格劃分對(duì)基于體積力自航性能計(jì)算結(jié)果的影響。

選取第1層邊界層網(wǎng)格厚度分別為0.000 02L、0.000 01L和0.000 005L，根據(jù)公式（6），其計(jì)算前預(yù)估值分別為1 058、528和264。相應(yīng)邊界層網(wǎng)格計(jì)算后得到的船體表面值實(shí)際分布如圖7所示。

圖7 船體表面y+分布

船體濕表面（水下部分）平均值依次為714、359和180，表明采用公式（6）估算的值偏大。船體表面高區(qū)集中在舵迎流面，船身大部分區(qū)域在相應(yīng)平均值左右。

如表5所示，第1層邊界層網(wǎng)格厚度0.000 02L時(shí)預(yù)報(bào)的螺旋槳轉(zhuǎn)速偏差較大，而第1層邊界層網(wǎng)格厚度為0.000 01L和0.000 005L時(shí)，兩者計(jì)算精度相近。因此，根據(jù)公式（6）在約為500計(jì)算得到的第1層邊界層網(wǎng)格厚度用于實(shí)尺度數(shù)值計(jì)算時(shí)，可達(dá)到較高計(jì)算精度。綜合計(jì)算時(shí)間和精度，其他研究中，將船體表面第1層邊界層網(wǎng)格厚度設(shè)置為0.000 01L。

表5 不同邊界層網(wǎng)格方案計(jì)算結(jié)果

3.2 湍流模型依賴性研究

選取渦黏模型-（SST-Menter）、雷諾應(yīng)力模型EASM和分離渦模型DES-SST，分析湍流模型對(duì)自航性能計(jì)算結(jié)果的影響。

如表6所示，-（SST-Menter）模型和EASM模型的螺旋槳轉(zhuǎn)速預(yù)報(bào)精度相差不大，DES-SST模型螺旋槳轉(zhuǎn)速預(yù)報(bào)精度更高，偏差為0.52%。每次迭代所消耗計(jì)算資源從小到大依次為-（SSTMenter）模型、EASM、DES-SST模型，在相同設(shè)置下，-（SST-Menter）模型計(jì)算耗時(shí)最短。綜合計(jì)算時(shí)間和精度，在其他研究中，湍流模型設(shè)置為-（SST-Menter）模型。

表6 不同湍流模型計(jì)算結(jié)果

3.3 激勵(lì)盤(pán)幾何參數(shù)影響研究

激勵(lì)盤(pán)幾何參數(shù)設(shè)置包括中心點(diǎn)坐標(biāo)、厚度、內(nèi)半徑和外半徑。激勵(lì)盤(pán)中心點(diǎn)坐標(biāo)為螺旋槳中心點(diǎn)坐標(biāo)，外半徑為螺旋槳半徑，這兩者取值是確定的；激勵(lì)盤(pán)內(nèi)半徑為槳轂半徑，而槳轂半徑通常沿槳軸方向是變化的；激勵(lì)盤(pán)厚度需要參考螺旋槳槳葉形狀，而各種類(lèi)型螺旋槳水動(dòng)力性能各異，其厚度值選取難以統(tǒng)一。因此，本文專門(mén)考察激勵(lì)盤(pán)內(nèi)半徑和厚度的取值對(duì)自航性能計(jì)算結(jié)果的影響。

槳轂半徑沿槳軸方向由0.800 m變化到0.850 m。另外，如果槳轂存在時(shí)，需保證激勵(lì)盤(pán)不與物面相交，因此要增加考慮激勵(lì)盤(pán)內(nèi)半徑大于槳轂最大半徑的情形，于是選取3個(gè)內(nèi)半徑：0.80 m、0.85 m和0.90 m。

如表7所示，在槳轂半徑變化范圍內(nèi)的內(nèi)半徑取值對(duì)計(jì)算結(jié)果影響很小，內(nèi)半徑越大預(yù)報(bào)轉(zhuǎn)速越高。內(nèi)半徑為0.80 m時(shí)，螺旋槳轉(zhuǎn)速預(yù)報(bào)偏差最小。由于本文中用于體積力自航計(jì)算的船體幾何中不包含槳轂，內(nèi)半徑按槳轂最小半徑取值時(shí)不存在激勵(lì)盤(pán)與槳轂相交的情況，所以在其他研究中，將激勵(lì)盤(pán)內(nèi)半徑設(shè)置為該值。

表7 不同激勵(lì)盤(pán)內(nèi)半徑計(jì)算結(jié)果

FINE/Marine推薦的激勵(lì)盤(pán)厚度為0.1，實(shí)際上本文螺旋槳槳葉側(cè)投影輪廓厚度約為0.15，于是選取3個(gè)厚度：0.92 m(0.1)、1.38 m(0.15)和1.84 m(0.2)。如表8所示，激勵(lì)盤(pán)厚度對(duì)計(jì)算結(jié)果影響較大，激勵(lì)盤(pán)厚度偏大或偏小均導(dǎo)致螺旋槳轉(zhuǎn)速預(yù)報(bào)偏差增大，而根據(jù)螺旋槳槳葉實(shí)際厚度確定的激勵(lì)盤(pán)厚度計(jì)算精度最高，所以在其他研究中，將激勵(lì)盤(pán)厚度設(shè)置為1.38 m(0.15)。

表8 不同激勵(lì)盤(pán)厚度計(jì)算結(jié)果

3.4 數(shù)值計(jì)算控制參數(shù)影響研究

針對(duì)不同的時(shí)間步長(zhǎng)和最大非線性迭代次數(shù)進(jìn)行計(jì)算，分析這些數(shù)值計(jì)算控制參數(shù)對(duì)自航性能計(jì)算結(jié)果的影響。

根據(jù)0.01L /、0.005L /和0.002 5L /選取3個(gè)時(shí)間步長(zhǎng)：0.25 s、0.125 s和0.062 5 s。如表9所示，在（0.002 5～0.01）L /范圍內(nèi)時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果影響不大，在其他研究中，將時(shí)間步長(zhǎng)設(shè)置為0.125 s（0.005L /）。

表9 不同時(shí)間步長(zhǎng)計(jì)算結(jié)果

最大非線性迭代次數(shù)指每個(gè)時(shí)間步長(zhǎng)內(nèi)在求解非線性控制方程過(guò)程中所設(shè)置的最大迭代次數(shù)。選取3個(gè)最大非線性迭代次數(shù)：5次、10次和15次。如下頁(yè)表10所示，最大非線性迭代次數(shù)對(duì)計(jì)算結(jié)果幾乎無(wú)影響，而計(jì)算時(shí)間與迭代次數(shù)正相關(guān)。綜合計(jì)算時(shí)間和精度，在其他研究中，最大非線性迭代次數(shù)設(shè)置為5次。

表10 不同最大非線性迭代步數(shù)計(jì)算結(jié)果

4 基于滑移網(wǎng)格的實(shí)尺度自航性能數(shù) 值計(jì)算方法研究

4.1 網(wǎng)格依賴性研究

采用滑移網(wǎng)格模擬真實(shí)槳轉(zhuǎn)動(dòng)時(shí)，由于螺旋槳相對(duì)船體旋轉(zhuǎn)，其螺旋槳表面流動(dòng)與船體表面流動(dòng)存在較大差異。因此，專門(mén)針對(duì)螺旋槳表面邊界層網(wǎng)格的不同劃分方案進(jìn)行計(jì)算，分析螺旋槳表面邊界層網(wǎng)格劃分對(duì)基于滑移網(wǎng)格自航性能計(jì)算結(jié)果的影響。船體表面第1層邊界層網(wǎng)格厚度設(shè)置為0.000 01L。

計(jì)算過(guò)程中，除航速固定為22 kn外，為節(jié)省計(jì)算時(shí)間，未啟用轉(zhuǎn)速控制器，而直接將螺旋槳轉(zhuǎn)速設(shè)置為與實(shí)船試航結(jié)果相同，即80.7 r/min。

螺旋槳表面第1層邊界層網(wǎng)格厚度依次設(shè)置為0.000 02L、0.000 01L和0.000 005L，相應(yīng)邊界層網(wǎng)格計(jì)算后得到的螺旋槳表面值實(shí)際分布如圖8所示，螺旋槳表面平均值依次為1 768、888和454。

圖8 螺旋槳表面y+分布

螺旋槳表面高區(qū)位于槳葉導(dǎo)邊靠近葉梢側(cè)，低區(qū)位于槳轂表面靠近舵球側(cè)，槳葉表面沿徑向由葉根到葉梢逐漸增加。

數(shù)值計(jì)算結(jié)果與實(shí)船試航結(jié)果的對(duì)比如表11所示，螺旋槳收到功率均用實(shí)船試航結(jié)果進(jìn)行了無(wú)量綱化。

表11 不同螺旋槳邊界層網(wǎng)格方案計(jì)算結(jié)果

對(duì)比發(fā)現(xiàn)不同螺旋槳邊界層網(wǎng)格方案螺旋槳收到功率偏差預(yù)報(bào)變化幅度達(dá)到6.99%，表明螺旋槳邊界層網(wǎng)格對(duì)自航計(jì)算精度影響較大。螺旋槳第1層邊界層網(wǎng)格厚度為0.000 01L時(shí)螺旋槳收到功率預(yù)報(bào)精度最高，幾乎與實(shí)船試航結(jié)果完全吻合。綜合計(jì)算時(shí)間和精度，其他研究中，將螺旋槳表面第1層邊界層網(wǎng)格厚度設(shè)置為0.000 01L。

4.2 湍流模型依賴性研究

選取不同的3種湍流模型，分析湍流模型對(duì)自航性能計(jì)算結(jié)果的影響。如表12所示，-（SSTMenter）模型螺旋槳收到功率預(yù)報(bào)精度最高，在其他研究中，湍流模型設(shè)置為-（SST-Menter）模型。

表12 不同湍流模型計(jì)算結(jié)果

4.3 動(dòng)區(qū)域幾何參數(shù)影響研究

采用滑移網(wǎng)格進(jìn)行自航計(jì)算時(shí)，需要在螺旋槳周?chē)?dòng)區(qū)域。動(dòng)區(qū)域除了包圍螺旋槳之外，要避免船體相交，如果動(dòng)區(qū)域半徑過(guò)小會(huì)使螺旋槳葉梢距離動(dòng)區(qū)域圓柱面過(guò)近，還會(huì)使船尾表面距離動(dòng)區(qū)域圓柱面過(guò)近，這2種情況均會(huì)限制動(dòng)區(qū)域附近物面網(wǎng)格尺寸，進(jìn)而影響網(wǎng)格質(zhì)量，由于船尾布置后者相對(duì)前者更少出現(xiàn)。因此，本文專門(mén)考察動(dòng)區(qū)域半徑的取值對(duì)自航性能計(jì)算結(jié)果的影響。

根據(jù)螺旋槳半徑與船尾幾何形狀，選取3個(gè)動(dòng)區(qū)域半徑：4.8 m、5.0 m和5.2 m。如表13所示，動(dòng)區(qū)域半徑會(huì)影響其周?chē)W(wǎng)格劃分，對(duì)計(jì)算結(jié)果影響較大，動(dòng)區(qū)域半徑為4.8 m螺旋槳收到功率預(yù)報(bào)偏差最大，另外兩者偏差相近，所以動(dòng)區(qū)域半徑選取時(shí)不宜過(guò)小。動(dòng)區(qū)域半徑為5.0 m時(shí)，螺旋槳轉(zhuǎn)速預(yù)報(bào)偏差最小，在其他研究中，將動(dòng)區(qū)域半徑設(shè)置為該值。

表13 不同動(dòng)區(qū)域半徑計(jì)算結(jié)果

4.4 數(shù)值計(jì)算控制參數(shù)影響研究

同樣針對(duì)不同的時(shí)間步長(zhǎng)和最大非線性迭代次數(shù)進(jìn)行計(jì)算，分析這些數(shù)值計(jì)算控制參數(shù)對(duì)自航性能計(jì)算結(jié)果的影響。

采用滑移網(wǎng)格進(jìn)行自航計(jì)算時(shí)，為了加快收斂，先采用大的時(shí)間步長(zhǎng)獲得初解，然后轉(zhuǎn)為小的時(shí)間步長(zhǎng)獲得精確解。大的時(shí)間步長(zhǎng)參考L /，此處與體積力方法一致取為0.005L /，即0.125 s；小的時(shí)間步長(zhǎng)參考螺旋槳旋轉(zhuǎn)周期1/，由于螺旋槳轉(zhuǎn)速為80.7 r/min時(shí)，轉(zhuǎn)過(guò)1°所需時(shí)間約為0.002 s，于是選取3個(gè)時(shí)間步長(zhǎng)：0.004 s、0.002 s和0.001 s。如表14所示，在螺旋槳轉(zhuǎn)過(guò)0.5°～2°所需時(shí)間范圍內(nèi)，時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果影響很小，在其他研究中，將時(shí)間步長(zhǎng)設(shè)置為0.002 s。

表14 不同時(shí)間步長(zhǎng)計(jì)算結(jié)果

同樣選取3個(gè)最大非線性迭代次數(shù)：5次、10次和15次。如表15所示，最大非線性迭代次數(shù)較小時(shí)，螺旋槳收到功率預(yù)報(bào)偏差較大，另外兩者偏差相近。這是由于真槳旋轉(zhuǎn)使船尾流場(chǎng)更加復(fù)雜，導(dǎo)致收斂變慢，所以有必要增加最大非線性迭代次數(shù)。綜合計(jì)算時(shí)間和精度，在其他研究中，最大非線性迭代次數(shù)設(shè)置為10次。

表15 不同最大非線性迭代步數(shù)計(jì)算結(jié)果

5 實(shí)船自航性能預(yù)報(bào)結(jié)果分析

為了對(duì)比體積力方法和滑移網(wǎng)格方法的實(shí)尺度自航性能數(shù)值計(jì)算結(jié)果，針對(duì)基于滑移網(wǎng)格的實(shí)尺度自航性能數(shù)值計(jì)算，在給定航速22 kn下，開(kāi)啟轉(zhuǎn)速控制器，根據(jù)船體阻力與螺旋槳推力平衡求解實(shí)船自航點(diǎn)。

2種方法所預(yù)報(bào)實(shí)船自航點(diǎn)螺旋槳轉(zhuǎn)速、收到功率如表16所示。船模自航試驗(yàn)結(jié)果根據(jù)1978年ITTC方法換算到實(shí)船。與實(shí)船試航結(jié)果相比，2種方法計(jì)算得到的螺旋槳轉(zhuǎn)速、收到功率與試驗(yàn)換算結(jié)果偏差更小，其中實(shí)槳轉(zhuǎn)速偏差均在3%以內(nèi)，實(shí)槳收到功率偏差均在7%以內(nèi)。

表16 2種方法螺旋槳轉(zhuǎn)速與收到功率的計(jì)算結(jié)果

由于計(jì)算和試驗(yàn)采用的是與設(shè)計(jì)槳敞水性征相似的備用槳，而實(shí)船試航采用的是設(shè)計(jì)槳。船舶自航時(shí)，螺旋槳在船處于非均勻來(lái)流條件下工作，設(shè)計(jì)槳與備用槳雖然在均勻來(lái)流中的敞水性能相近，但是在非均勻來(lái)流中它們性能可能存在差異，從而導(dǎo)致計(jì)算結(jié)果與實(shí)船試航結(jié)果的偏差比其與試驗(yàn)換算結(jié)果的偏差大。體積力方法得到的螺旋槳轉(zhuǎn)速、收到功率相比于滑移網(wǎng)格方法更接近試驗(yàn)換算結(jié)果和實(shí)船試航結(jié)果，表明實(shí)尺度下采用體積力模型替代真實(shí)螺旋槳進(jìn)行實(shí)船性能預(yù)估仍能具有較高精度。

結(jié)合圖4的實(shí)槳敞水性征曲線與實(shí)尺度阻力性能數(shù)值計(jì)算結(jié)果，2種方法根據(jù)各自實(shí)尺度自航性能數(shù)值計(jì)算結(jié)果，采用等推力法求解的自航推進(jìn)因子如表17所示。表中：是根據(jù)等推力法從實(shí)槳敞水性征曲線插值所得進(jìn)速系數(shù)、K是推力系數(shù)、K是轉(zhuǎn)矩系數(shù)、t是推力減額、w是伴流分?jǐn)?shù)、η是推進(jìn)效率、η是船身效率、η是相對(duì)旋轉(zhuǎn)效率、是敞水效率。

表17 2種方法推進(jìn)因子的計(jì)算結(jié)果

除推力減額和伴流分?jǐn)?shù)外，2種方法計(jì)算得到的其余各項(xiàng)推進(jìn)因子與試驗(yàn)換算偏差值均在5%以內(nèi)，體積力方法過(guò)于低估了螺旋槳推力，而滑移網(wǎng)格方法過(guò)于低估了螺旋槳轉(zhuǎn)矩，兩者計(jì)算推進(jìn)效率相近。體積力方法不能得到相對(duì)旋轉(zhuǎn)效率，因?yàn)闊o(wú)法確定螺旋槳在船后時(shí)的轉(zhuǎn)矩，因此體積力方法的相對(duì)旋轉(zhuǎn)效率η=1。

2種方法計(jì)算得到的尾流場(chǎng)渦結(jié)構(gòu)如下頁(yè)圖9所示，取速度梯度二階不變量=5的等值面。體積力方法捕捉到了轂渦，在槳盤(pán)面處的渦環(huán)表明體積力模擬出了螺旋槳的螺旋性?；凭W(wǎng)格方法不僅捕捉到了轂渦，還捕捉到了槳葉葉梢處的梢渦。

圖9 體積力方法（左）與滑移網(wǎng)格方法（右）Q=5等值面的渦結(jié)構(gòu)

為了進(jìn)一步比較體積力與真實(shí)螺旋槳對(duì)尾流場(chǎng)的影響，在槳盤(pán)面前后0.15處選取截面（如圖10所示），截面形狀為與螺旋槳同軸、內(nèi)直徑為0.30、外直徑為1.00的圓盤(pán)。

圖10 尾流場(chǎng)截面位置

2種方法計(jì)算得到的槳前0.15處截面無(wú)量綱化軸向、切向和徑向速度場(chǎng)如下頁(yè)圖11所示。從船尾向前看，軸向速度V以指向船尾為正，切向速度V以順時(shí)針?lè)较驗(yàn)檎?，徑向速?span id="j5i0abt0b" class="emphasis_italic">V以向外為正。體積力方法得到的軸向、切向和徑向速度場(chǎng)基本呈對(duì)稱分布，未體現(xiàn)切向體積力對(duì)槳前近槳盤(pán)面處水流的作用。滑移網(wǎng)格方法模擬真實(shí)螺旋槳的旋轉(zhuǎn)使槳前近槳盤(pán)面處水流產(chǎn)生了扭轉(zhuǎn)，導(dǎo)致槳前軸向、切向和徑向速度場(chǎng)非對(duì)稱分布?；凭W(wǎng)格方法槳前徑向速度場(chǎng)分布結(jié)構(gòu)與體積力方法類(lèi)似。

圖11 體積力方法（上）與滑移網(wǎng)格方法（下）槳前0.15D處截面無(wú)量綱化軸向、切向和徑向速度分布

2種方法計(jì)算得到的槳后0.15處截面無(wú)量綱化軸向、切向和徑向速度場(chǎng)如下頁(yè)圖12所示。體積力方法得到的槳后軸向、切向和徑向速度場(chǎng)不再呈對(duì)稱分布，在切向體積力作用下速度場(chǎng)分布呈現(xiàn)出一定偏轉(zhuǎn)。滑移網(wǎng)格方法得到的槳后軸向、切向和徑向速度場(chǎng)在螺旋槳旋轉(zhuǎn)作用下也呈非對(duì)稱分布。兩者槳后軸向速度分布結(jié)構(gòu)類(lèi)似，盤(pán)面內(nèi)徑和外徑附近軸向速度小，內(nèi)徑與外徑之間軸向速度較大，但是體積力方法槳后軸向速度沿周向分布不均勻，而滑移網(wǎng)格方法槳后軸向速度沿周向均勻分布。兩者槳后切向速度分布結(jié)構(gòu)也比較相似，盤(pán)面左上側(cè)切向速度較大，且兩者切向速度較大區(qū)域的形狀相差不大，而盤(pán)面右下側(cè)切向速度由內(nèi)到外逐漸減小，到外徑附近切向速度與螺旋槳旋轉(zhuǎn)方向相反。兩者槳后徑向速度分布相差較大，主要是因?yàn)轶w積力方法只考慮了軸向力和切向力而未考慮徑向力。

圖12 體積力方法（上）與滑移網(wǎng)格方法（下）槳后0.15D處截面無(wú)量綱化軸向、切向和徑向速度分布

2種方法計(jì)算得到的槳前和槳后0.15處截面無(wú)量綱化橫向(切向、徑向)速度矢量分布分別如圖13和圖14所示。兩者槳前橫向速度矢量分布差異明顯，體積力方法槳前橫向速度矢量分布幾乎對(duì)稱，而滑移網(wǎng)格方法的橫向速度矢量周向分布不均勻，在槳葉附近橫向速度較大，且其盤(pán)面上方兩側(cè)的小渦比前者更靠近槳轂。兩者槳后橫向速度矢量分布結(jié)構(gòu)較為相似，由于自航點(diǎn)轉(zhuǎn)速差異，兩者橫向速度大小存在差異，此外，滑移網(wǎng)格方法的徑向流動(dòng)更為明顯。

圖13 體積力方法（左）與滑移網(wǎng)格方法（右）槳前0.15D處截面無(wú)量綱化切向、徑向速度矢量分布

圖14 體積力方法（左）與滑移網(wǎng)格方法（右）槳后0.15D處截面無(wú)量綱化切向、徑向速度矢量分布

本文采用體積力方法和滑移網(wǎng)格方法對(duì)9 400 TEU集裝箱船進(jìn)行實(shí)尺度自航性能數(shù)值模擬。通過(guò)數(shù)值計(jì)算方法研究發(fā)現(xiàn)，網(wǎng)格劃分對(duì)計(jì)算精度有較大影響，其中船身表面實(shí)際平均值在300左右、螺旋槳表面實(shí)際平均值在500左右時(shí)，計(jì)算精度較高；采用激勵(lì)盤(pán)方法時(shí)，激勵(lì)盤(pán)厚度選取對(duì)計(jì)算精度影響較大，建議參考螺旋槳槳葉側(cè)投影輪廓厚度確定；采用滑移網(wǎng)格方法時(shí)，動(dòng)區(qū)域大小選取對(duì)計(jì)算精度影響較大，需避免動(dòng)區(qū)域距離槳葉過(guò)近。

6 結(jié) 語(yǔ)

體積力方法和滑移網(wǎng)格方法的實(shí)船自航性能預(yù)報(bào)結(jié)果與模型試驗(yàn)換算及實(shí)船試航結(jié)果吻合較好，表明2種方法用于實(shí)尺度數(shù)值計(jì)算預(yù)估實(shí)船性能時(shí)具有較高的可靠性。相對(duì)于滑移網(wǎng)格方法，體積力方法由于網(wǎng)格數(shù)量更少計(jì)算效率更高，但是無(wú)法得到螺旋槳相對(duì)旋轉(zhuǎn)效率并且額外需要螺旋槳敞水?dāng)?shù)據(jù)。兩者對(duì)螺旋槳模擬方式不同造成尾流場(chǎng)的差異，對(duì)于槳前流動(dòng)，兩者主要差異在于軸向和切向流動(dòng)，體積力方法未能體現(xiàn)螺旋槳旋轉(zhuǎn)效應(yīng)；而對(duì)于槳后流動(dòng)，兩者主要差異在于徑向流動(dòng)，體積力方法未考慮徑向力的作用。

本文進(jìn)行實(shí)尺度自航性能數(shù)值模擬時(shí)，未考慮船體表面粗糙度和空氣阻力的影響，后續(xù)需進(jìn)一步研究。