王英第，陳彥臻，周偉健，張盛龍

(1. 渤海船舶職業(yè)學(xué)院 船舶工程學(xué)院，遼寧 葫蘆島 125001；2. 上海海事大學(xué) 商船學(xué)院，上海 201306；3. 哈爾濱工程大學(xué) 動(dòng)力與能源工程學(xué)院，黑龍江 哈爾濱 150001；4. 常熟理工學(xué)院 汽車學(xué)院，江蘇 常熟 215500)

0 引 言

自從能效設(shè)計(jì)指數(shù)（EEDI）標(biāo)準(zhǔn)實(shí)施以來，學(xué)者們紛紛開始關(guān)注船舶的能源利用和環(huán)境保護(hù)問題。在船舶運(yùn)營(yíng)中，燃料的消耗不僅和能源利用、環(huán)保等問題息息相關(guān)，燃料成本也是公司的巨大開銷之一。因此，設(shè)計(jì)出阻力更小的船舶以減少船舶運(yùn)營(yíng)成本、提高燃料利用效率，成為各國(guó)船舶研究人員的共同目標(biāo)。能夠在設(shè)計(jì)階段準(zhǔn)確預(yù)報(bào)出船體阻力變得非常重要。

隨著計(jì)算機(jī)技術(shù)的高速發(fā)展以及數(shù)值理論的日益成熟，昂貴的船模實(shí)驗(yàn)開始漸漸被建立在仿真基礎(chǔ)上的設(shè)計(jì)方法所替代。近年來計(jì)算流體力學(xué)（CFD）方法得到了廣泛的認(rèn)可，隨著許多優(yōu)秀的CFD軟件如Fluent、CFX、SHIPFLOW等在工程問題中被大量應(yīng)用，CFD方法體現(xiàn)出了其巨大價(jià)值，也使CFD方法應(yīng)用于船舶阻力預(yù)報(bào)成為可能[1-3]。CFD方法能夠獲得與試驗(yàn)數(shù)據(jù)相符的計(jì)算結(jié)果，同時(shí)計(jì)算周期短，成本低廉[4-6]。近年來，CFD方法已經(jīng)成為計(jì)算船舶阻力的一種有效途徑，有相當(dāng)多的CFD方法被應(yīng)用在船體阻力預(yù)報(bào)中，基于CFD的數(shù)值計(jì)算成為一種新的經(jīng)濟(jì)、高效的研究船舶水動(dòng)力性能的工具。倪崇本[7]采用基于CFD理論的動(dòng)網(wǎng)格技術(shù)對(duì)某多體船進(jìn)行了數(shù)值模擬，探討了多體船船體姿態(tài)對(duì)船體表面壓力及船體附近波形等方面的影響。張艷等[8]利用CFD方法算得船身周圍水壓分布情況，并預(yù)測(cè)了水流對(duì)船舶產(chǎn)生的航行阻力?？捉鹌降萚9]用CFD軟件對(duì)雙螺旋槳系統(tǒng)進(jìn)行了數(shù)值計(jì)算，并采用湍流模型預(yù)報(bào)螺旋槳的水動(dòng)力性能。

為了準(zhǔn)確預(yù)報(bào)船體阻力，首先以自由模為例，討論了邊界層劃分方式及湍流模型中系數(shù)的選擇對(duì)阻力預(yù)報(bào)的影響。然后，在合適的網(wǎng)格方式上，采用NLPQL方法建立基于湍流系數(shù)的阻力優(yōu)化系統(tǒng)來獲得最佳湍流系數(shù)。其次，根據(jù)最合適的網(wǎng)格劃分方式及數(shù)值計(jì)算方法模擬不同航速下KCS船在海面上的水動(dòng)力性能，討論約束模和自由模在流場(chǎng)捕捉及阻力預(yù)報(bào)方面的差異。

1 CFD數(shù)學(xué)模型建立

整個(gè)CFD數(shù)學(xué)模型的流場(chǎng)控制方程由連續(xù)方程和運(yùn)動(dòng)方程構(gòu)成。采用Realizable k-ε模型作為整個(gè)流場(chǎng)的湍流模型。采用VOF[10]方法精確捕捉自由水平面，VOF法是一種在固定的歐拉網(wǎng)格下的表面跟蹤方法，通過求解動(dòng)量方程和一處或多處流體的體積分?jǐn)?shù)來模擬多項(xiàng)流模型。采用SIMPLE法對(duì)壓力場(chǎng)和速度場(chǎng)進(jìn)行耦合。

為了精確模擬KCS船在海面上的運(yùn)動(dòng)情況，將船體縱搖和垂蕩運(yùn)動(dòng)添加到阻力預(yù)報(bào)之中。在建立船舶運(yùn)動(dòng)方程時(shí)，建立2個(gè)參考坐標(biāo)系：一個(gè)是固定于大地的固定坐標(biāo)系，另一個(gè)是固定于船體的隨船動(dòng)坐標(biāo)系，其中動(dòng)坐標(biāo)系的原點(diǎn)在船體的重心處。

2 幾何模型及網(wǎng)格劃分

本文以KCS船型作為研究對(duì)象，該船的主尺度如表1所示。

整個(gè)計(jì)算域尺度設(shè)置為5Lpp×3Lpp×4Lpp。為了減小網(wǎng)格數(shù)量，取半船作為研究對(duì)象。在水池入口邊界給定船體航行速度，水池兩側(cè)設(shè)為對(duì)稱邊界，水池頂部和底部設(shè)為速度入口，出口處設(shè)為壓力出口，邊界條件設(shè)置如圖1所示。

切割體網(wǎng)格技術(shù)同時(shí)包含了非結(jié)構(gòu)網(wǎng)格和結(jié)構(gòu)網(wǎng)格的優(yōu)點(diǎn)，對(duì)模擬大幅度物體運(yùn)動(dòng)問題具有較高的計(jì)算精度，是一種新發(fā)展起來的網(wǎng)格生成技術(shù)。本文采用切割體網(wǎng)格技術(shù)劃分計(jì)算域，網(wǎng)格劃分如圖2～圖4所示。圖2為整個(gè)流域網(wǎng)格圖，自由液面附近網(wǎng)格加密。圖3為船體及舵網(wǎng)格劃分圖，船體首尾及舵網(wǎng)格加密。圖4為自由液面網(wǎng)格劃分圖，船體附近網(wǎng)格加密，以準(zhǔn)確捕捉船體周圍的流場(chǎng)。

表1 KCS船型主尺度Tab. 1 Principal dimensions of KCS ship types

圖1 邊界條件設(shè)置Fig. 1 Boundary conditions

圖2 整個(gè)計(jì)算域網(wǎng)格劃分Fig. 2 Mesh on the computational domain

圖3 船體及舵表面網(wǎng)格劃分Fig. 3 Mesh on the hull and rudder surface

圖4 自由液面網(wǎng)格劃分Fig. 4 Mesh on the free surface

3 阻力預(yù)報(bào)流程

以KCS為例，采用RANS法模擬不同航速下船體水動(dòng)力性能，其中船體分別設(shè)置成約束模和自由模，具體計(jì)算流程如圖5所示。

圖5 計(jì)算流程圖Fig. 5 Computational flow chart

4 CFD模型參數(shù)影響分析

4.1 邊界層對(duì)阻力的影響

船體在水中直航時(shí)，船體壁面附近會(huì)受到流體黏性的影響。對(duì)于黏性很小的空氣和水來說，黏性對(duì)流動(dòng)的影響僅限于貼近船體表面的一個(gè)薄層內(nèi)，這一薄層以外的黏性影響可以完全忽略。因此，為了準(zhǔn)確預(yù)報(bào)KCS船型阻力，這一層網(wǎng)格的劃分顯得非常重要。以無(wú)因次參數(shù)y+表示船體表面第1層網(wǎng)格節(jié)點(diǎn)的高度，計(jì)算公式為：

式中：y為第1層網(wǎng)格節(jié)點(diǎn)距離船體表面的高度（首層邊界層厚度）；Lpp為船體垂線間長(zhǎng)；Re為相對(duì)船體長(zhǎng)度所定義的雷諾數(shù)；y+取值范圍在30～200。

在邊界層網(wǎng)格劃分中，首層邊界層厚度、網(wǎng)格增長(zhǎng)因子和邊界層數(shù)決定了網(wǎng)格劃分方式。因此，按照這3個(gè)因子依次討論邊界層劃分方式與阻力之間的關(guān)系，計(jì)算結(jié)果如表2～表4所示。

表2 首層邊界層厚度與阻力之間的關(guān)系Tab. 2 The relation between the thickness of the first boundary layer and the resistance

表3 網(wǎng)格增長(zhǎng)因子與阻力之間的關(guān)系Tab. 3 The relationship between grid growth factor and resistance

表4 邊界層數(shù)與阻力之間的關(guān)系Tab. 4 The relationship between the number ofboundary layers and resistance

從計(jì)算結(jié)果可以看出，誤差變化幅值均低于10%。其中，總阻力隨著首層邊界層厚度的增加變化明顯，而隨著網(wǎng)格增長(zhǎng)因子和邊界層數(shù)的增加變化較小。隨著首層邊界層厚度的增加，船體總阻力增加；隨著邊界層數(shù)的增加，船體總阻力減?。欢W(wǎng)格增長(zhǎng)因子與船體總阻力之間沒有明顯的線性關(guān)系。從表3可以看出，當(dāng)網(wǎng)格增長(zhǎng)因子為1.5時(shí)，計(jì)算結(jié)果與實(shí)驗(yàn)值更接近。

4.2 湍流系數(shù)對(duì)阻力的影響

根據(jù)經(jīng)驗(yàn)，網(wǎng)格數(shù)量的增加會(huì)提高阻力預(yù)報(bào)的準(zhǔn)確性。但是隨著網(wǎng)格數(shù)量的不斷增加，計(jì)算效率將大大降低，并且網(wǎng)格數(shù)量增加到一定程度反而導(dǎo)致計(jì)算精度下降。所以在合適的網(wǎng)格數(shù)量上探討高精度的求解方式成為CFD數(shù)值模擬的關(guān)鍵。

雖然許多湍流模型已經(jīng)取得了一定計(jì)算效果，但迄今為止仍舊沒有得到有效而統(tǒng)一的湍流模型，并且在分析湍流模型方程的一系列參數(shù)對(duì)計(jì)算結(jié)果及不確定度的影響這方面，很少有學(xué)者進(jìn)行深入研究。因此，以Realizable k-ε模型作為研究對(duì)象，其中該模型的推薦計(jì)算值如表5所示。在湍流系數(shù)取值范圍附近改變?nèi)≈挡⒎謩e用來計(jì)算船舶阻力，其計(jì)算工況及阻力結(jié)果如表6～表10所示。

從計(jì)算結(jié)果可以看出，隨著湍流系數(shù)的變化，船體總阻力變化范圍在2.5%以內(nèi)。隨著Cμ和αk的增加，船體總阻力減小；隨著C2ε和αε的增加，船體總阻力增加；而隨著C1ε的增加，阻力計(jì)算結(jié)果沒有明顯的線性關(guān)系。當(dāng)C1ε＞1.42時(shí)，阻力隨著C1ε的增加而減小。

表5 Realizable k-ε湍流模型系數(shù)Tab. 5 Realizable k-ε turbulence model coefficients

表6 Cμ與阻力之間關(guān)系Tab. 6 The relationship between Cμ and resistance

表7 C1ε與阻力之間關(guān)系Tab. 7 The relationship between C1ε and resistance

表8 C2ε與阻力之間關(guān)系Tab. 8 The relationship between C2ε and resistance

表10 αε與阻力之間關(guān)系Tab. 10 The relationship between and resistance

Caseαε 相對(duì)偏差/%CtError/%Case330.96-200.003 64-1.91 Case341.08-100.003 665-1.24 Case351.200.003 685-0.7 Case361.32100.003 695-0.43 Case371.44200.003 7230.32

5 模型參數(shù)優(yōu)化及多航速阻力預(yù)報(bào)

5.1 模型參數(shù)優(yōu)化

為了尋找到更適合KCS船體阻力預(yù)報(bào)的湍流系數(shù)，構(gòu)建基于NLPQL算法[11]的優(yōu)化系統(tǒng)。由之前的討論可以看出，C2ε和αε對(duì)阻力的影響較大，而C1ε和αk的變化對(duì)阻力影響較小。為了減小設(shè)計(jì)變量提高優(yōu)化效率，固定C1ε=1.44，αk=1.0，選取其余3個(gè)參數(shù)作為優(yōu)化設(shè)計(jì)變量。目標(biāo)函數(shù)||RtCFD-RtEFD||為最小值。優(yōu)化流程圖如圖6所示，優(yōu)化過程如圖7所示。最優(yōu)解集如表11所示。可以看出，經(jīng)過一系列的優(yōu)化得到了最優(yōu)湍流系數(shù)。按照最優(yōu)解計(jì)算的總阻力與實(shí)驗(yàn)值最接近，誤差僅有0.054%。

圖6 優(yōu)化流程Fig. 6 Optimization process

圖7 序列二次規(guī)劃算法優(yōu)化過程Fig. 7 The optimization process for the NLPQL algorithm

表11 模型系數(shù)選取Tab. 11 The selection of Realizable k-ε turbulence model coefficients

5.2 多航速阻力預(yù)報(bào)

通過以上研究與討論，明確了船體阻力預(yù)報(bào)的邊界層劃分方式和最優(yōu)湍流系數(shù)取值，以下采用該方法預(yù)報(bào)KCS船在靜水中直航的總阻力。

分別根據(jù)推薦系數(shù)和最優(yōu)系數(shù)對(duì)船體總阻力進(jìn)行預(yù)報(bào)，表12給出了CFD計(jì)算結(jié)果與實(shí)驗(yàn)值的對(duì)比?？梢钥闯觯w設(shè)置成自由模型比約束模型計(jì)算精度高，平均誤差為0.95%。且在高航速時(shí)，約束模型誤差較大而自由模型阻力預(yù)報(bào)值與實(shí)驗(yàn)值吻合較好，誤差只有1.38%。

表12 總阻力計(jì)算結(jié)果與實(shí)驗(yàn)值對(duì)比Tab. 12 The comparison of the Ct calculated by CFD and EFD methods

圖8和圖9為船體表面壓力分布圖。從圖中可以發(fā)現(xiàn)，在每個(gè)航速下，船體表面壓力等值線分布變化不大，但在球鼻首部分壓力等值線有明顯的不同，此處產(chǎn)生的壓阻是興波阻力的主要組成部分。由此可知，只有考慮了船體自身的運(yùn)動(dòng)姿態(tài)才能更全面的反應(yīng)船體實(shí)際航行特性。在自由模計(jì)算中，最優(yōu)系數(shù)在船體總阻力預(yù)報(bào)上平均誤差只有0.43%以內(nèi)，而推薦系數(shù)的阻力計(jì)算結(jié)果誤差在0.95%?？梢钥闯觯顑?yōu)系數(shù)在阻力預(yù)報(bào)精度明顯高于推薦系數(shù)，只有在Fr=0.282時(shí)計(jì)算結(jié)果較差，以此驗(yàn)證了該湍流系數(shù)在多航速船體阻力預(yù)報(bào)上的實(shí)用性與可靠性。

圖8 船首表面壓力分布（約束模—推薦系數(shù)）Fig. 8 Static pressure distribution on the bow section（restraint mode-recommendation coefficients）

圖9 船首表面壓力分布（自由?！扑]系數(shù)）Fig. 9 Static pressure distribution on the bow section（free moderecommended coefficients）

表13給出了船體穩(wěn)定時(shí)刻的縱傾角度和垂蕩位移隨Fr變化的數(shù)值大小?？梢钥闯?，本文計(jì)算結(jié)果與文獻(xiàn)[12]中的計(jì)算結(jié)果較吻合。隨著航速的增加，船體升沉明顯增大，而縱傾角增加較為緩慢。圖10為計(jì)算穩(wěn)定時(shí)刻船體航行姿態(tài)。

表13 縱傾值和垂蕩變化值與文獻(xiàn)值對(duì)比Tab. 13 The comparison of pitch and heave values obtained by CFD method and previous literature

圖10 穩(wěn)定時(shí)刻船體姿態(tài)（自由模—最優(yōu)系數(shù)）Fig. 10 Ship attitude at stable time（free mode-optimal coefficients）

圖11 推薦系數(shù)波形等高線對(duì)比圖Fig. 11 The comparison of the waveform contour obtained by recommended coefficients

圖11為采用推薦系數(shù)計(jì)算的自由模與約束模波形等高線對(duì)比圖。可以看出，2種方法均能較清晰地捕捉到船體周圍的波浪形狀，說明VOF方法在自由液面捕捉上的優(yōu)勢(shì)。雖然2種波形圖相似，但是在船首附近的波形相差較大。這主要是由于穩(wěn)定后自由模在流場(chǎng)中的航行姿態(tài)與約束模不同，使得自由模的船體首尾壓力差大于約束模。圖12為采用自由模計(jì)算的KCS船波形等高線對(duì)比圖。可以看出，在Fr=0.282工況下波形等高線較為相似，而在Fr=0.260時(shí)波形等高線差別較大。這與表12中阻力計(jì)算結(jié)果相符。可見波形等高線的捕捉同樣可以作為船體阻力預(yù)報(bào)正確與否的評(píng)價(jià)標(biāo)準(zhǔn)之一。

圖12 自由模波形等高線對(duì)比圖Fig. 12 The comparison of the waveform contour obtained by free-mode

6 結(jié) 語(yǔ)

本文采用SIMPLE算法求解不可壓縮RANS方程，對(duì)KCS船在靜水中的阻力預(yù)報(bào)精度進(jìn)行研究，模擬了不同F(xiàn)r時(shí)KCS船受到的總阻力、船體表面及周圍流場(chǎng)分布情況。結(jié)果表明：

1）船體附近邊界層劃分方式對(duì)阻力預(yù)報(bào)有較大的影響。隨著首層邊界層厚度的增加，阻力也隨之增加；隨著邊界層數(shù)的增加，阻力減?。欢W(wǎng)格增長(zhǎng)因子與阻力之間沒有明顯的線性關(guān)系。

2）C2ε和αε對(duì)阻力的影響較大，而C1ε和αk的變化對(duì)阻力影響較小。隨著Cμ和αk的增加，阻力減??；隨著C2ε和αε的增加，阻力增加；而C1ε與阻力之間沒有明顯的線性關(guān)系，但當(dāng)C1ε＞1.42時(shí)，阻力隨著C1ε的增加而減小。其次，構(gòu)建了基于NLPQL法的優(yōu)化系統(tǒng)，得到了適合KCS船體阻力預(yù)報(bào)的湍流系數(shù)：Cμ=0.09，C1ε=1.44，C2ε=1.950 5，αk=1.0，αε=1.187 5。

3）自由模在船體阻力預(yù)報(bào)精確度上優(yōu)于約束模，誤差只有0.95%。采用最優(yōu)湍流系數(shù)對(duì)船體阻力進(jìn)行預(yù)報(bào)時(shí)，除了Fr=0.282以外，其他航速阻力計(jì)算精確度均高于推薦湍流系數(shù)。