實船阻力及流場數(shù)值預報方法

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

(海軍工程大學 動力工程學院， 湖北 武漢 430033)

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

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

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

1 基于CFD的船舶總阻力預報方法

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

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

(1)

(2)

式中：ui分別為流體質點在i方向的速度分量，fi是質量力，p為流體的壓力，μ是相體積分數(shù)平均的動力粘度系數(shù)，μt為湍流動力粘性系數(shù)。

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

1.2 船舶阻力試驗換算方法

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

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

(3)

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

(4)

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

(5)

(6)

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

(7)

式中：Lpp為垂線間長，ks為粗糙度表觀高度，一般可取ks=0.15 mm。

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

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

(8)

(9)

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

(10)

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

2 計算結果及分析

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

本文的主要研究對象是DTMB 5415，該船型是ITTC推薦的軍艦類型中唯一有大量公開試驗數(shù)據(jù)的平臺。美國的DTMB，IIHR以及意大利的INSEAN三家研究機構對該船型進行了全面的試驗。本文對DTMB 5415進行了流場模擬和阻力計算，并以試驗數(shù)據(jù)[7-8]為校驗。因DTMB 5415沒有相關的實船數(shù)據(jù)，本文暫將模型放大20倍作為實船。

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

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

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

表2 DTMB 5415疊模阻力計算

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

表3 1+k法預報的阻力

其中相對偏差定義為

(11)

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

2.2 “虛流體粘度”方法

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

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

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

hs/Ls=hm/Lm

(12)

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

(13)

表5 虛流體粘度方法阻力預報

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

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

圖1 粗糙度對摩擦阻力的影響

圖2 粗糙度對壓阻力影響

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

2.3 實尺度船舶RANS計算

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

表6 實尺度計算阻力預報

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

2.4 波形和流場比較和分析

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

圖3 船體附近興波

波高的比較結果見圖4。定義原點為首垂線與水線的交點，x軸指向船艉，y軸指向船體右側，z軸指向船體上方。

圖4 Y/Lpp=0.172處的波高

通過y/Lpp=0.172處波高的對比，可以看出：

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

2)虛流體粘度給出的結果與實尺度計算得到的波高吻合較好，表明虛流體粘度的方法可以較好地預報實船的波高。

3)波高在0～Lpp處都吻合較好，但是在船艉處實船雷諾數(shù)下波幅要大于模型雷諾數(shù)下的波幅，這與之前波形的分析結果也是一致的。因為在模型雷諾數(shù)下，船艉處的邊界層厚度相對較大，粘性對興波的影響更加明顯，導致波幅較小。

圖5為槳盤面處伴流場，可以看出：

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

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

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

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

3 結論

本文基于CFD對DTMB 5415實船尺度下的阻力和流場進行了計算和預報。通過對計算結果的進一步分析，得到以下結論：

1)通過合理地設置船體表面第一層網(wǎng)格厚度和粗糙度，可以通過虛流體粘度方法或實尺度計算準確地求取實船的阻力。

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

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

隨著尺度效應和粗糙度相關研究和試驗的開展，本文的研究工作有待進一步深入。

參考文獻：

[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]劉志華, 熊鷹, 韓寶玉. 雷諾相似船模預報實船推進因子的數(shù)值方法[J]. 哈爾濱工程大學學報,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]倪崇本, 朱仁傳, 繆國平, 等.一種基于CFD的船舶總阻力預報方法[J].水動力學研究與進展,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]盛振邦, 劉應中. 船舶原理[M]. 上海：上海交通大學出版社,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.