易文彬，王永生，楊瓊方，李劍

(海軍工程大學(xué) 動(dòng)力工程學(xué)院， 湖北 武漢 430033)

船阻力預(yù)報(bào)大部分是基于模型試驗(yàn)進(jìn)行的，但是模型試驗(yàn)一般是在傅汝德數(shù)相等的條件下進(jìn)行的，與實(shí)船雷諾數(shù)并不相等。在實(shí)際工程中，可以通過(guò)模型試驗(yàn)測(cè)得模型的總阻力，然后通過(guò)外推經(jīng)驗(yàn)公式(二因次法或三因次法)和一些船體粗糙度等有關(guān)的經(jīng)驗(yàn)修正項(xiàng)得到實(shí)船阻力。這些方法具有很強(qiáng)的工程實(shí)用性，但是包含了大量的經(jīng)驗(yàn)修正項(xiàng)，船模與實(shí)船相關(guān)換算缺乏嚴(yán)格的理論支撐，也不能給出實(shí)船雷諾數(shù)下的流場(chǎng)信息。

目前基于CFD的船舶數(shù)值計(jì)算主要集中于模型尺度的雷諾數(shù)(107)，對(duì)實(shí)尺雷諾數(shù)(109)下的研究還比較少。但是隨著船舶計(jì)算流體力學(xué)的發(fā)展和工程實(shí)踐的需要，實(shí)尺度雷諾數(shù)下阻力和流場(chǎng)的精確計(jì)算獲得了更加越來(lái)越多的關(guān)注[1]。Raven等[2]通過(guò)CFD的方法計(jì)算了模型和實(shí)船的阻力并且分析了船模實(shí)船阻力換算過(guò)程中各阻力成分的尺度效應(yīng)。Bhushan等[3]通過(guò)Athena船模和實(shí)船的數(shù)值計(jì)算與校驗(yàn)，驗(yàn)證了一種新型的壁面函數(shù)模型。其中實(shí)尺計(jì)算值考慮了粗糙度的影響，與船模試驗(yàn)外推值吻合良好。劉志華等[4]提出了一種基于設(shè)置流體運(yùn)動(dòng)粘度系數(shù)實(shí)現(xiàn)雷諾數(shù)相似的計(jì)算理論，能夠較為快速準(zhǔn)確地求取實(shí)船雷諾數(shù)下的自航因子。倪崇本等[5]提出將實(shí)尺度下勢(shì)流理論與模型尺度下湍流理論相結(jié)合的方法來(lái)求取實(shí)船的阻力，雖然阻力預(yù)報(bào)有較高的精度，但是不能給出實(shí)船雷諾數(shù)下的流場(chǎng)信息。

本文采用基于數(shù)值計(jì)算的3種方法預(yù)報(bào)了實(shí)船阻力、波形及流場(chǎng)信息，對(duì)計(jì)算結(jié)果進(jìn)行了比較和分析。在虛流體粘度方法中分析了粗糙度對(duì)各阻力成分的影響，比較了阻力與流場(chǎng)的計(jì)算結(jié)果與其他方法計(jì)算結(jié)果的差異。

1 基于CFD的船舶總阻力預(yù)報(bào)方法

1.1 數(shù)值計(jì)算模型及控制方程

采用有限體積法離散控制方程，不可壓縮粘性流體的控制方程如下：

(1)

(2)

式中：ui分別為流體質(zhì)點(diǎn)在i方向的速度分量，fi是質(zhì)量力，p為流體的壓力，μ是相體積分?jǐn)?shù)平均的動(dòng)力粘度系數(shù)，μt為湍流動(dòng)力粘性系數(shù)。

上述方程需要結(jié)合湍流模型封閉方程組，常見(jiàn)的兩方程湍流模型有k-ε、RNGk-ε、k-ω、k-ωSST、RSM等。本文選取的是k-ωSST湍流模型。

1.2 船舶阻力試驗(yàn)換算方法

船模阻力試驗(yàn)換算方法最早由傅汝德提出，該方法將船體總阻力分成摩擦阻力Rf和剩余阻力Rr，并認(rèn)為摩擦阻力僅與雷諾數(shù)Re有關(guān)，剩余阻力僅與傅汝德數(shù)Fr有關(guān)，且兩者互不干擾[6]。因此總阻力系數(shù)可以表示為

Ct(Fr,Re)=Cf(Re)+Cr(Fr)

(3)

因此實(shí)船尺度下的阻力系數(shù)可由船模阻力試驗(yàn)中測(cè)得的阻力系數(shù)表示為

(4)

式中：Cf0代表相當(dāng)平板摩擦阻力系數(shù)，Cr為剩余阻力系數(shù)，下標(biāo)s和m分別代表變量對(duì)應(yīng)于實(shí)船和模型。在此基礎(chǔ)上，休斯提出了三因次換算方法，引入形狀因子1+k來(lái)表示船體粘性阻力和相當(dāng)平板摩擦阻力之間的關(guān)系：

(5)

(6)

引入形狀因子的三因次法較二因次法合理，其中形狀因子可由低速拖曳或疊模試驗(yàn)確定。考慮到實(shí)船船體表面粗糙度的作用以及船模實(shí)船阻力換算過(guò)程中由于雷諾數(shù)不相等導(dǎo)致的尺度效應(yīng)，需要補(bǔ)貼一定的經(jīng)驗(yàn)系數(shù)，即船模實(shí)船換算補(bǔ)貼系數(shù)Ca，Ca可由下式確定：

(7)

式中：Lpp為垂線(xiàn)間長(zhǎng)，ks為粗糙度表觀高度，一般可取ks=0.15 mm。

1.3 船模實(shí)船雷諾數(shù)和傅汝德數(shù)全相似條件

為保證船模和實(shí)船總阻力系數(shù)相等，必須滿(mǎn)足實(shí)船和船模的雷諾數(shù)Re和傅汝德數(shù)Fr相等：

(8)

(9)

式中：L為船長(zhǎng)，υ為流體運(yùn)動(dòng)粘度，V為來(lái)流速度，g為重力加速度。在船模試驗(yàn)時(shí)不能滿(mǎn)足全相似定律，因此工程實(shí)踐中船模試驗(yàn)都是在傅汝德數(shù)相等的條件下進(jìn)行的。為實(shí)現(xiàn)實(shí)船和船模的雷諾數(shù)和傅汝德數(shù)全相似，在數(shù)值計(jì)算中可令

(10)

此時(shí)，式(8)、(9)同時(shí)成立。引入虛擬流體運(yùn)動(dòng)粘度系數(shù)之后，流體仍為牛頓流體，流場(chǎng)的基本運(yùn)動(dòng)規(guī)律不變，因此湍流模型及控制方程不需要調(diào)整。

2 計(jì)算結(jié)果及分析

2.1 基于模型尺度數(shù)值模擬外推方法

本文的主要研究對(duì)象是DTMB 5415，該船型是ITTC推薦的軍艦類(lèi)型中唯一有大量公開(kāi)試驗(yàn)數(shù)據(jù)的平臺(tái)。美國(guó)的DTMB，IIHR以及意大利的INSEAN三家研究機(jī)構(gòu)對(duì)該船型進(jìn)行了全面的試驗(yàn)。本文對(duì)DTMB 5415進(jìn)行了流場(chǎng)模擬和阻力計(jì)算，并以試驗(yàn)數(shù)據(jù)[7-8]為校驗(yàn)。因DTMB 5415沒(méi)有相關(guān)的實(shí)船數(shù)據(jù)，本文暫將模型放大20倍作為實(shí)船。

計(jì)算區(qū)域入口取船艏向前延伸1倍船長(zhǎng)處，出口取船艉向后延伸2倍船長(zhǎng)處，側(cè)邊界及下方邊界均取1倍船長(zhǎng)，上方邊界取水線(xiàn)上方0.5倍船長(zhǎng)。采用VOF的方法來(lái)追蹤自由液面，湍流模型采用k-ωSST模型，對(duì)流項(xiàng)離散選用二階離散格式。本文采用分塊結(jié)構(gòu)化網(wǎng)格，船體周?chē)捎肙型網(wǎng)格以便捕捉船體附近的邊界層，同時(shí)對(duì)船艏、船艉及靜水面進(jìn)行網(wǎng)格加密，網(wǎng)格節(jié)點(diǎn)數(shù)320萬(wàn)。滿(mǎn)足計(jì)算要求。

Fr為0.15、0.28、0.41 共3個(gè)航速下阻力計(jì)算值與試驗(yàn)值如表1所示。其中Fr=0.28及Fr=0.41航速下考慮了船體姿態(tài)的變化。從表1可以看出：通過(guò)數(shù)值模擬能夠較好地預(yù)報(bào)船模的阻力，阻力計(jì)算值誤差在3%以?xún)?nèi)。若采用三因次換算方法求取實(shí)船阻力，還必須求得形狀因子1+k。ITTC建議用傅汝德數(shù)低于0.15的船模阻力試驗(yàn)數(shù)據(jù)來(lái)求取形狀因子。本文通過(guò)疊模計(jì)算求取傅汝德數(shù)在0.1～0.15時(shí)的形狀因子，計(jì)算結(jié)果如表2。

表1 DTMB 5415船?？傋枇ο禂?shù)

表2 DTMB 5415疊模阻力計(jì)算

取多個(gè)速度下形狀因子的平均值，1+k=1.083。由低速時(shí)拖曳船模試驗(yàn)數(shù)據(jù)通過(guò)普魯哈斯卡方法換算得到1+k=1.15。疊模計(jì)算得到的形狀因子1+k偏小，因?yàn)榇５退偻弦吩囼?yàn)時(shí)興波雖然不明顯，但是仍然存在興波阻力成分，所以得到的形狀因子比疊模計(jì)算值大?？紤]船模實(shí)船換算補(bǔ)貼Ca(按照式(7)本船取為0.000 509)，將模型試驗(yàn)與數(shù)值計(jì)算得到的船模阻力系數(shù)Ctm按照三因次法(1+k分別取1.15和1.083)外推得到的實(shí)船阻力系數(shù)Cts。其值如表3所示。

表3 1+k法預(yù)報(bào)的阻力

其中相對(duì)偏差定義為

(11)

基于模型尺度數(shù)值模擬外推方法得到的實(shí)船阻力與由模型試驗(yàn)外推得到的阻力差別在7%以?xún)?nèi)。因疊模計(jì)算得到的形狀因子1+k與通過(guò)低速拖曳試驗(yàn)得到的1+k相比偏小，如果采用相同的船模實(shí)船換算補(bǔ)貼，通過(guò)計(jì)算預(yù)報(bào)的實(shí)船阻力與通過(guò)試驗(yàn)預(yù)報(bào)的實(shí)船阻力存在較大的偏差。

2.2 “虛流體粘度”方法

為實(shí)現(xiàn)船模、實(shí)船雷諾數(shù)Re和傅汝德數(shù)Fr全相似條件，在對(duì)船模進(jìn)行數(shù)值模擬時(shí)將水的運(yùn)動(dòng)粘度按照全相似的條件設(shè)置為一個(gè)虛擬的粘度系數(shù)(本文將水的運(yùn)動(dòng)粘度系數(shù)取為1.067×10-8m2/s)。文獻(xiàn)[9]表明，在船模雷諾數(shù)下，如果要準(zhǔn)確模擬船模的摩擦阻力，船體表面第一層網(wǎng)格的厚度要滿(mǎn)足y+在30～ 300之間，但是在實(shí)船雷諾數(shù)下，y+最大值可能遠(yuǎn)遠(yuǎn)大于1 000。于是本文基于船體摩擦阻力變化對(duì)船體表面第一層網(wǎng)格的厚度做了相關(guān)的研究：本文針對(duì)Fr=0.28航速下，劃分了5套網(wǎng)格，船體表面第一層網(wǎng)格厚度分別為0.2、1、2、3、5 mm，計(jì)算得到的摩擦阻力系數(shù)如表4。

表4 不同網(wǎng)格下摩擦阻力計(jì)算結(jié)果

從表4可以看出，針對(duì)本算例情形，船體表面第一層網(wǎng)格厚度應(yīng)該使y+在5 000～ 10 000時(shí)能夠準(zhǔn)確模擬實(shí)船雷諾數(shù)下的摩擦阻力。進(jìn)行實(shí)船阻力計(jì)算，還必須要考慮粗糙度的影響。目前在數(shù)值計(jì)算中，一般是采用平均沙粒粗糙度來(lái)代替船舶實(shí)際粗糙度的影響，文獻(xiàn)[10-11]表明將實(shí)船平均沙粒粗糙度取為hs=0.03 mm與粗糙度表觀高度ks=0.15 mm是等效的。若船模實(shí)船相對(duì)粗糙度相等，則

hs/Ls=hm/Lm

(12)

將船模平均沙粒粗糙度hm取為0.001 5 mm。光船阻力系數(shù)和考慮粗糙度計(jì)算得到的實(shí)船阻力系數(shù)如表5。定義相對(duì)偏差為

(13)

表5 虛流體粘度方法阻力預(yù)報(bào)

從表5可以看出：考慮粗糙度的虛流體粘度系數(shù)方法阻力計(jì)算值與船模試驗(yàn)外推值在傅汝德數(shù)從0.15～0.41航速范圍內(nèi)相差不到4.5%，表明這種方法能夠較好的預(yù)報(bào)實(shí)船的阻力。在計(jì)算中同時(shí)考慮粗糙度和實(shí)船雷諾數(shù)的作用，克服了船模實(shí)船換算的尺度效應(yīng)，可以直接求得實(shí)船的阻力系數(shù)。此方法計(jì)算所需的網(wǎng)格與普通模型數(shù)值模擬的網(wǎng)格量接近，且不需要進(jìn)行疊模計(jì)算求取形狀因子，可以快速準(zhǔn)確地實(shí)現(xiàn)實(shí)船阻力預(yù)報(bào)。

為進(jìn)一步分析粗糙度的影響，將考慮粗糙度前后各阻力成分的變化用圖1，圖2表示如下。

圖1 粗糙度對(duì)摩擦阻力的影響

圖2 粗糙度對(duì)壓阻力影響

不考慮粗糙度時(shí)，由虛流體粘度方法計(jì)算得到的摩擦阻力系數(shù)與相當(dāng)平板摩擦阻力系數(shù)是很接近的，壓阻力系數(shù)與模型雷諾數(shù)下計(jì)算得到的壓阻力系數(shù)差別不大。考慮粗糙度之后，發(fā)現(xiàn)摩擦阻力系數(shù)增大了12%～31%，壓阻力系數(shù)基本不變，表明粗糙度對(duì)摩擦阻力影響較明顯，對(duì)壓阻力影響不大。

2.3 實(shí)尺度船舶RANS計(jì)算

針對(duì)實(shí)尺度計(jì)算，采用與模型尺度下相似的控制域，并且對(duì)船體和水線(xiàn)面附近的網(wǎng)格進(jìn)行加密，網(wǎng)格數(shù)目為1 500萬(wàn)。為縮短計(jì)算時(shí)間，本文以層流的方法得到的結(jié)果作為高精度湍流計(jì)算的初值，以光滑船體的計(jì)算結(jié)果作為粗糙船體計(jì)算的初值，實(shí)踐證明，按照這種方法可以極大地縮短計(jì)算時(shí)間。實(shí)船的平均沙粒粗糙度取0.03 mm。不考慮粗糙度和考慮粗糙度的實(shí)船阻力計(jì)算結(jié)果如表6。

表6 實(shí)尺度計(jì)算阻力預(yù)報(bào)

考慮粗糙度的實(shí)船阻力CFD計(jì)算值與船模試驗(yàn)外推值偏差在3.5%以?xún)?nèi)，表明實(shí)船的數(shù)值模擬能夠提供較精確的實(shí)船阻力。實(shí)尺度計(jì)算不需要相關(guān)的經(jīng)驗(yàn)系數(shù)補(bǔ)貼，能夠克服船模實(shí)船換算的尺度效應(yīng)，提供實(shí)尺度下的流場(chǎng)信息。隨著計(jì)算機(jī)技術(shù)的發(fā)展，實(shí)尺度下的數(shù)值模擬表現(xiàn)了越來(lái)越強(qiáng)的工程應(yīng)用意義。

2.4 波形和流場(chǎng)比較和分析

圖3為Fr=0.28航速下船體的興波，船模的數(shù)值模擬與船模試驗(yàn)值的波形基本吻合。數(shù)值計(jì)算預(yù)報(bào)的波形比較接近，但是波幅有細(xì)微的差別，其中實(shí)船雷諾數(shù)下船艉的波幅更大。

圖3 船體附近興波

波高的比較結(jié)果見(jiàn)圖4。定義原點(diǎn)為首垂線(xiàn)與水線(xiàn)的交點(diǎn)，x軸指向船艉，y軸指向船體右側(cè)，z軸指向船體上方。

圖4 Y/Lpp=0.172處的波高

通過(guò)y/Lpp=0.172處波高的對(duì)比，可以看出：

1)模型試驗(yàn)與船模數(shù)值計(jì)算得到波峰和波谷的位置吻合較好，但是計(jì)算得到的首波峰較小，可能是數(shù)值耗散的原因。在1.5Lpp以后由于網(wǎng)格較疏，波高與試驗(yàn)值有一定的偏離。

2)虛流體粘度給出的結(jié)果與實(shí)尺度計(jì)算得到的波高吻合較好，表明虛流體粘度的方法可以較好地預(yù)報(bào)實(shí)船的波高。

3)波高在0～Lpp處都吻合較好，但是在船艉處實(shí)船雷諾數(shù)下波幅要大于模型雷諾數(shù)下的波幅，這與之前波形的分析結(jié)果也是一致的。因?yàn)樵谀Ｐ屠字Z數(shù)下，船艉處的邊界層厚度相對(duì)較大，粘性對(duì)興波的影響更加明顯，導(dǎo)致波幅較小。

圖5為槳盤(pán)面處伴流場(chǎng)，可以看出：

1)模型尺度數(shù)值模擬能夠和試驗(yàn)數(shù)據(jù)較好的吻合，但是在一些數(shù)據(jù)點(diǎn)處仍有偏離。

2)虛粘度方法與實(shí)船計(jì)算給出了較為一致的結(jié)果，且伴流作用較模型雷諾數(shù)下弱。主要原因是模型雷諾數(shù)下船尾的邊界層較厚，對(duì)伴流的影響更加顯著。

圖5 槳盤(pán)處縱向伴流(X/Lpp=0.935,Z/Lpp=-0.02)

通過(guò)對(duì)波形、波高和伴流場(chǎng)的比較可以看出：虛流體粘度方法和實(shí)尺度計(jì)算得出的結(jié)果較為一致且體現(xiàn)出了實(shí)船雷諾數(shù)下流場(chǎng)特征；粘性的作用相對(duì)較弱，導(dǎo)致伴流較弱，興波更強(qiáng)。

3 結(jié)論

本文基于CFD對(duì)DTMB 5415實(shí)船尺度下的阻力和流場(chǎng)進(jìn)行了計(jì)算和預(yù)報(bào)。通過(guò)對(duì)計(jì)算結(jié)果的進(jìn)一步分析，得到以下結(jié)論：

1)通過(guò)合理地設(shè)置船體表面第一層網(wǎng)格厚度和粗糙度，可以通過(guò)虛流體粘度方法或?qū)嵆叨扔?jì)算準(zhǔn)確地求取實(shí)船的阻力。

2)虛流體粘度方法和實(shí)尺度計(jì)算均能直接求取實(shí)尺度雷諾數(shù)下的阻力和流場(chǎng)，但虛流體粘度方法網(wǎng)格量少，計(jì)算量較小，表現(xiàn)了較強(qiáng)的工程意義。

3)基于船模尺度數(shù)值模擬外推值與船模試驗(yàn)外推值偏差較大，主要原因是形狀因子的計(jì)算值與試驗(yàn)采用的數(shù)值差別較大。有關(guān)形狀因子的計(jì)算有待進(jìn)一步探討。

隨著尺度效應(yīng)和粗糙度相關(guān)研究和試驗(yàn)的開(kāi)展，本文的研究工作有待進(jìn)一步深入。

參考文獻(xiàn)：

[1]STERN F, YANG Jianming, WANG Zhaoyuan, et al. Com-

putational ship hydrodynamics: nowadays and way forward[C]//29th Symposium on Nnaval Hhydrodynamics. Gothenburg,Sweden, 2012:26-31.

[2]RAVEN H C, VAN D, STARKE A R, et al. Towards a CFD-based prediction of ship performance—Progress in predicting full-scale resistance and scale effects[J]. International Journal of Maritime Engineering, 2009,150(5):14-18.

[3]BHUSHAN S, XING T, CARRICA P, et al. Model- and full-scale URANS simulation of Athena resistance, powering, seakeeping, and 5415 maneuvering[J]. Journal of Ship Research,2009,53(4): 179-198.

[4]劉志華, 熊鷹, 韓寶玉. 雷諾相似船模預(yù)報(bào)實(shí)船推進(jìn)因子的數(shù)值方法[J]. 哈爾濱工程大學(xué)學(xué)報(bào),2008, 29(7):658-662.

LIU Zhihua, XIONG Ying, HAN Baoyu. Accurately predicting the propulsion factors of a ship using a Reynolds-similarity ship model[J]. Journal of Harbin Engineering University, 2008,29(7):658-662.

[5]倪崇本, 朱仁傳, 繆國(guó)平, 等.一種基于CFD的船舶總阻力預(yù)報(bào)方法[J].水動(dòng)力學(xué)研究與進(jìn)展,A輯, 2010, 25(5):579-586.

NI Chongben, ZHU Renchuan, MIAO Guoping, et al. A method for ship resistance prediction based on CFD computation[J]. Chinese Journal of Hydrodynamics, 2010,25(5):579-586.

[6]盛振邦, 劉應(yīng)中. 船舶原理[M]. 上海：上海交通大學(xué)出版社,2003:154-157.

[7]OLIVIERI A, PISTANI F, AVANZINI A, et al. Towing tank experiments of resistance, sinkage and trim, boundary layer, wake, and free surface flow around a naval combatant INSEAN 2340 model[R].Iowa:The University of Iowa,2005.

[8]STERN F, LONGO J, PENNA R, et al. International collaboration on benchmark CFD validation data for surface combatant DTMB model 5415[C]//Proceedings of 23rd ONR symposium on naval hydrodynamics. Iowa,2000:17-22.

[9]ANSYS Inc. ANSYS CFX-solver theory guide[M]. (s.l.):ANSYS Inc, 2007.

[10]HE Chengbin, HUANG Jianguo, HAN Jin, et al. Cyclic shift keying spread spectrum underwater acoustic communication[J]. Acta Phys Sin, 2009, 58(12): 8379-8385.

[11]SCHULTZ M P. Effects of coating roughness and biofouling on ship resistance and powering[J]. Biofouling, 2007,23(13):331-341.

[12]ECA L, HOEKSTRA M, RAVEN H C. Quantifying roughness effects by ship viscous flow calculations[C]// 28th Symposium on Naval Hydrodynamics. California, 2010:12-17.