翟錦果，辛榮斌，汪宗御，張繼鋒,2，張 海，紀(jì)玉龍

(1. 大連海事大學(xué) 輪機(jī)工程學(xué)院，遼寧 大連 116026；2. 浙江清華長三角研究院，浙江 嘉興 314006；3. 清華大學(xué) 能源與動力工程系，北京 100084)

0 引 言

在設(shè)計(jì)過程中，對船舶航行阻力的計(jì)算必不可少，如何快速準(zhǔn)確的估算出船舶航行阻力是船舶行業(yè)一直以來重點(diǎn)關(guān)注的問題[1]。目前，關(guān)于船舶阻力的估算方法有船模試驗(yàn)法、近似值估算法以及數(shù)值模擬法[2]。隨著現(xiàn)代船舶用途以及船型等的改變，很多近似值估算法已經(jīng)不能滿足實(shí)際工程的計(jì)算需要。而對船舶進(jìn)行水池試驗(yàn)雖然精確度較高，但是成本高且耗費(fèi)周期長。隨著計(jì)算流體力學(xué)以及計(jì)算機(jī)仿真技術(shù)的不斷發(fā)展，采用數(shù)值模擬方法來計(jì)算船舶阻力逐漸成為主流。

目前關(guān)于船舶數(shù)值模擬的軟件有兩類：一種是采用勢流理論、發(fā)展較為成熟的專業(yè)軟件，如Shipflow；另一類是采用粘性理論的商用數(shù)值計(jì)算軟件，如Ansys-Fluent、STAR-CCM+等。勢流理論對于以前比較簡單的船舶型線的阻力計(jì)算具有計(jì)算快、精度高的特點(diǎn)。但是，由于現(xiàn)在船舶的型線不斷發(fā)展改變，采用勢流理論計(jì)算船舶阻力的準(zhǔn)確性有待提高。因此，目前更多學(xué)者開始采用基于粘性理論的商用軟件STAR CCM+[3-6]。雖然商用軟件均具備給出誤差較小的計(jì)算結(jié)果的能力，但是在數(shù)值模擬過程中，由于計(jì)算的復(fù)雜度，精確度很大程度上取決于計(jì)算策略。例如，同樣是對DTMB45 船模計(jì)算，劉飛[7]的計(jì)算區(qū)域距船首前為1.5 倍船長，船尾后為2.5 倍船長，自由液面上下部分均為2 倍船長，對船舶的波形區(qū)域進(jìn)行加密并選取標(biāo)準(zhǔn)k-Epsilon 湍流模型，最后計(jì)算結(jié)果的誤差在3.31%以內(nèi)；而尹輝[5]選取的計(jì)算區(qū)域距船首前約1 倍船長，距船尾后約3 倍船長，整個(gè)計(jì)算域的上部為一倍船長，下部為2 倍船長，采用K-Omega 湍流模型，傅汝德數(shù)為0.2～0.35時(shí)，計(jì)算結(jié)果的誤差在2.65% 以下，但是在弗汝德數(shù)為0.05 時(shí)誤差則達(dá)到16%。羅良[8]采用了Star-CCM+對KCS 船舶進(jìn)行數(shù)值模擬，計(jì)算區(qū)域在船長方向?yàn)? 倍船長，SST KOmega 湍流模型，仿真誤差在2%以內(nèi)。倪崇本等[9]采用K-Espilon 湍流模型對KCS 進(jìn)行數(shù)值計(jì)算，計(jì)算誤差達(dá)到6.5%。

對于船舶數(shù)值模擬來說，目前常用的湍流模型有RANS 模型中的SST K-Omega 湍流模型、標(biāo)準(zhǔn)K-Epsilon 湍流模型、標(biāo)準(zhǔn)K-Omega 湍流模型、可實(shí)現(xiàn)的KEspilon 湍流模型4 種，總的來說，這些湍流模型預(yù)測船舶總阻力較為可靠[10-15]。

由此可知，不同的計(jì)算策略對于船舶阻力性能的精確度有著不同程度的影響[16]。但是目前對于數(shù)值模擬中計(jì)算區(qū)域大小、網(wǎng)格密度以及湍流模型對于數(shù)值模擬的總阻力、摩擦阻力和耗費(fèi)時(shí)長的研究還存在不足，因此本文對船舶數(shù)值模擬過程中的這些影響因素進(jìn)行對比分析，提出一種滿足實(shí)際工程需求的數(shù)值模擬方法。

1 數(shù)值計(jì)算

1.1 數(shù)值計(jì)算模型

本文以某30 000 t 的散貨船為例進(jìn)行數(shù)值建模，船舶的設(shè)計(jì)參數(shù)如表1 所示，為了減少計(jì)算量，數(shù)值計(jì)算模型與實(shí)際船體比例尺采用29.3∶1。由于船舶在靜水中的直線航行是對稱的，為了減少計(jì)算量又不影響計(jì)算精度，僅對一半的船舶進(jìn)行仿真計(jì)算。將船尾軸中心線和船舶設(shè)計(jì)吃水線的交點(diǎn)設(shè)置為坐標(biāo)原點(diǎn)，船尾指向船首方向?yàn)閄軸正方向，船左舷方向?yàn)閅軸正方向，垂直向上為Z軸正方向。仿真模擬在處理器型號為Inter Xeon CPU E5-2687W v4@3.00 GHz、內(nèi)存128 GB的惠普服務(wù)器上進(jìn)行，數(shù)值模擬采用40 進(jìn)程，約占CPU 利用率的80%。

表1 某30 000 t 船舶參數(shù)Tab. 1 Parameters of 30 000 t ship

1.2 控制方程

數(shù)值計(jì)算模型基于RANS 湍流模型，使用商業(yè)CFD 軟件包STAR-CCM+開發(fā)，流場中采用的控制方程為不可壓縮流體的連續(xù)性方程和動量方程。

連續(xù)性方程：

動量方程：

1.3 邊界條件的選取

本文對計(jì)算流場的邊界條件設(shè)置如圖1 所示。將對稱面設(shè)置為對稱邊界條件，頂部設(shè)置為壓力出口邊界，其他面均設(shè)為速度入口邊界，船體表面設(shè)置為無滑移的壁面。對計(jì)算區(qū)域進(jìn)口、出口、以及側(cè)面設(shè)置進(jìn)行Wave Forcing 消波，Wave Forcing 可以在指定的距離內(nèi)，迫使離散化N-S 方程的解趨向于另一種解，相比于阻尼消波，這種方法既可以消除區(qū)域邊界上的波浪反射，又可以設(shè)置較小的計(jì)算區(qū)域從而減少網(wǎng)格數(shù)量，節(jié)省計(jì)算時(shí)長。

圖1 數(shù)值計(jì)算區(qū)域劃分Fig. 1 Division of numerical calculation area

1.4 網(wǎng)格劃分

通過創(chuàng)建自動網(wǎng)格，選擇面網(wǎng)格、切割體網(wǎng)格以及棱柱層網(wǎng)格。為了提高計(jì)算精確度、捕捉船舶航行中波浪的變化，采用切割體網(wǎng)格單元生成器，對自由液面處x、y和z方向的網(wǎng)格進(jìn)行劃分，其中x、y、z方向的網(wǎng)格比例為8∶8∶1。

1.5 物理模型的選擇

采用VOF（Volume of Fluid）法來模擬船舶航行中自由液面，VOF 方法在網(wǎng)格單元內(nèi)設(shè)置流體體積分?jǐn)?shù)的值 αi來定義自由液面，定義為：

式中：Vi為一個(gè)網(wǎng)格單元體中流體的體積分?jǐn)?shù)；V為一個(gè)網(wǎng)格單元的體積。

在數(shù)值仿真過程中，一個(gè)網(wǎng)格單元體中不同的流體的體積分?jǐn)?shù)值的和為1。在船舶航行的過程中，自由液面處的流體包含水和空氣，因此，在數(shù)值計(jì)算過程中，在自由液面處設(shè)置值為0.5 來捕捉自由液面，如圖2 所示，水的體積分?jǐn)?shù)用于反映計(jì)算單元的狀態(tài)。值為0.5 表示計(jì)算單元中填充了50% 的水和50%的空氣，代表自由表面。值0 和1 表示計(jì)算單元分別充滿空氣和水。

圖2 船體水的體積分?jǐn)?shù)Fig. 2 Volume fraction of water on the hull

2 數(shù)值計(jì)算的影響因素分析

2.1 計(jì)算區(qū)域大小對計(jì)算結(jié)果的分析

首先以航速19 kn，試驗(yàn)阻力值924 600 N 為例，對30 000 t 散貨船模型尺度下的靜水阻力進(jìn)行數(shù)值模擬。采用了計(jì)算區(qū)域尺寸從4 倍船長到8 倍船長的5 種尺寸進(jìn)行數(shù)值模擬，計(jì)算區(qū)域內(nèi)寬度約為2 倍船長，自由液面以下約為2 倍船長，自由液面以上約1 倍船長。其中，網(wǎng)格基礎(chǔ)尺寸均為0.14 m，采用SST K-Omega 湍流模型。計(jì)算時(shí)間步長為0.03，計(jì)算時(shí)間為90 s，每個(gè)時(shí)間步長迭代次數(shù)為5，在整個(gè)數(shù)值模擬過程迭代1.5 萬次。

采用三因次換算法[2]將數(shù)值模擬總阻力換算為實(shí)船阻力，并與實(shí)船試驗(yàn)值進(jìn)行對比。為了減少換算誤差，摩擦阻力采用的模型阻力值，計(jì)算結(jié)果如表2 所示。其他條件不變時(shí)，隨著計(jì)算區(qū)域的增加，網(wǎng)格數(shù)量也增加，計(jì)算時(shí)長也不斷增加，誤差先減小后增大，在6 倍船長時(shí)誤差達(dá)到最小值1.39%。但是計(jì)算區(qū)域大小對于船舶摩擦阻力的影響基本不變。因此，在對船舶進(jìn)行數(shù)值計(jì)算過程中，推薦采用6 倍計(jì)算區(qū)域。

表2 計(jì)算區(qū)域?qū)?shù)值計(jì)算結(jié)果的影響Tab. 2 Influence of calculation area on numerical results

2.2 網(wǎng)格密度對計(jì)算結(jié)果的影響

在數(shù)值模擬過程中，分別取誤差最大的4 倍船長和誤差最小的6 倍船長計(jì)算區(qū)域，船速為19 kn，其他條件不變，通過改變網(wǎng)格的基礎(chǔ)尺寸的基礎(chǔ)設(shè)置，以此來改變網(wǎng)格密度。探究網(wǎng)格密度對于數(shù)值模擬的精確度的影響。

由表3 可知，在基礎(chǔ)尺寸是0.12 和0.14 時(shí)，6 倍船長的誤差均小于4 倍船長?；A(chǔ)尺寸為0.12 的4 倍船長誤差和基礎(chǔ)尺寸為0.14 的6 倍船長相差不大，因此減少基礎(chǔ)尺寸，適當(dāng)?shù)脑黾佑?jì)算區(qū)域，可以在保證計(jì)算精度的情況下，減少計(jì)算時(shí)長。計(jì)算區(qū)域6 倍船長，基礎(chǔ)尺寸為0.12 時(shí)，誤差最小，達(dá)到-0.62%，但是基礎(chǔ)尺寸為0.16 時(shí)的誤差達(dá)到1.39%，也可以滿足實(shí)際工程，且整體耗時(shí)也較短。

表3 網(wǎng)格密度對數(shù)值計(jì)算結(jié)果的影響Tab. 3 Influence of grid density on numerical results

2.3 湍流模型對計(jì)算結(jié)果的影響

為探究不同湍流模型對船舶數(shù)值計(jì)算精度及耗時(shí)的影響，在航速19 kn，計(jì)算區(qū)域?yàn)? 倍船長，網(wǎng)格大小為157.6 萬的情況下，分別采用SST K-Omega、可實(shí)現(xiàn)的K-Espilon 湍流和標(biāo)準(zhǔn)K-Omega 等3 種模型，進(jìn)行數(shù)值計(jì)算和時(shí)間估算，結(jié)果如表4所示。

表4 湍流模型對數(shù)值計(jì)算結(jié)果的影響Tab. 4 Influence of turbulence model on numerical results

可知，SST K-Omega 在船舶數(shù)值模擬過程中，對船舶的數(shù)值計(jì)算的精確度較高，誤差僅有1.39%，而且耗時(shí)最短，這是由于SST K-Omega 是標(biāo)準(zhǔn)的KOmega 和標(biāo)準(zhǔn)K-Epsilon 湍流模型結(jié)合優(yōu)化后的結(jié)果。不同的湍流模型對于摩擦阻力的影響較大，可實(shí)現(xiàn)的K-Espilon 湍流和SST K-Omega 對于摩擦阻力的數(shù)值計(jì)算較為接近，這是由于這2 個(gè)模型相對于標(biāo)準(zhǔn)的標(biāo)準(zhǔn)K-Omega 模型增加了額外的計(jì)算項(xiàng)。因此，本文推薦采用SST K-Omega 湍流模型來進(jìn)行船舶總阻力的數(shù)值計(jì)算。

3 數(shù)值計(jì)算方法驗(yàn)證

3.1 數(shù)值計(jì)算結(jié)果可視化分析

以表2 中6 倍船長計(jì)算結(jié)果為例，圖3 為船舶數(shù)值計(jì)算過程中的波形圖?？芍?，隨著時(shí)間的增加，船舶航行過程中的開爾文波形逐漸擴(kuò)大成型并且趨于穩(wěn)定狀態(tài)，波形圖符合船舶實(shí)際航行過程中的開爾文波形。從圖4(d)中船首、尾放大圖知，船首波高大于船尾波高，這也是粘壓阻力產(chǎn)生的重要原因。

圖3 船舶航行波形圖Fig. 3 Ship motion waveform

圖4 數(shù)值計(jì)算結(jié)果對比Fig. 4 Comparison of numerical calculation results

3.2 不同航速下船舶阻力變化

為了驗(yàn)證上述數(shù)值模擬的準(zhǔn)確性以及適用性，對船舶在14～21 kn，即船模速度為1.3312～1.9968 m/s 海況下的阻力進(jìn)行數(shù)值計(jì)算，實(shí)船與船模的船速采用傅汝德?lián)Q算法[2]，船模到實(shí)船阻力值采用三因次換算法，計(jì)算結(jié)果如圖4 所示。隨著船速的增加，船舶總阻力逐漸增大。數(shù)值模擬得到的船舶航行阻力變化趨勢和試驗(yàn)值一致，誤差均在5%以內(nèi)。導(dǎo)致誤差的原因可能是實(shí)船到船模尺寸阻力換算過程中產(chǎn)生的尺度效應(yīng)。由于誤差較小，表明該方法對船舶靜水阻力的計(jì)算精確度較高，滿足工程要求。

船舶靜水阻力可以分為摩擦阻力和粘壓阻力，根據(jù)相關(guān)的平板拖曳試驗(yàn)結(jié)果，摩擦阻力系數(shù)計(jì)算有以下幾種常用的經(jīng)驗(yàn)公式。

桑海公式：

柏蘭特-許立汀公式：

休斯公式：

ITTC1957 公式：

式中：Re=VsLwl/v為雷諾數(shù)；Vs為船速，m/s；Lwl為船舶設(shè)計(jì)水線長；v為海水運(yùn)動粘度。

圖5 為船模尺寸下數(shù)值模擬值摩擦阻力與各個(gè)經(jīng)驗(yàn)公式值的對比?？芍?，桑海公式與仿真計(jì)算值更接近，這是由于經(jīng)驗(yàn)公式都是根據(jù)船舶試驗(yàn)值分析得到，桑海公式更適合與本文相似船型的船舶。因此對同本文相似類型的船舶阻力進(jìn)行經(jīng)驗(yàn)公式計(jì)算時(shí)，推薦采用桑海公式進(jìn)行船舶摩擦阻力的計(jì)算。圖6 為數(shù)值計(jì)算中不同航速情況下船舶摩擦阻力和粘壓阻力的阻力值以及占比情況?？芍S著船速的增加，船舶摩擦阻力和粘壓阻力逐漸增大，但是隨著船速增加摩擦阻力占總阻力比值逐漸減小，粘壓阻力逐漸增大。這是因?yàn)槟Σ磷枇χ饕屠字Z數(shù)有關(guān)，航速越大，雷諾數(shù)越大，摩擦阻力系數(shù)越小，而粘壓阻力是因?yàn)檎承砸鸬拇w前后壓力不平衡而產(chǎn)生，船速越高，壓差越大，粘壓阻力越大，但是摩擦阻力還是處于主體地位。

圖5 摩擦阻力對比Fig. 5 Friction resistance comparison

圖6 阻力占比Fig. 6 Resistance ratio

3.3 KCS 船舶航速驗(yàn)證

為了驗(yàn)證本文的建模方法對其他船舶阻力計(jì)算同樣適用，因此本文對KSC 船舶進(jìn)行數(shù)值模擬，采用了6 倍船長計(jì)算區(qū)域，由于該船模長度略長于上述船模，故等比例采用0.15 m 的基礎(chǔ)尺寸，SST K-Omega湍流模型。模型尺度和實(shí)船尺度比例為31.6∶1，表5為KCS 船舶實(shí)船尺度相關(guān)參數(shù)。

表5 KCS 船舶參數(shù)Tab. 5 KCS ship parameters

KCS 船舶實(shí)驗(yàn)數(shù)據(jù)來自文獻(xiàn)[17]，由于文獻(xiàn)中將實(shí)驗(yàn)數(shù)據(jù)轉(zhuǎn)化為阻力系數(shù)，因此本文也將計(jì)算結(jié)果按照相同方法轉(zhuǎn)換為阻力系數(shù)進(jìn)行對比。

式中：Cd為總阻力系數(shù)；Fd為模型總阻力； ρ為水的密度；v為船模船速；A為船模表面積。

數(shù)值模擬誤差如表6 所示。船舶數(shù)值計(jì)算誤差均在3.5%以內(nèi)，滿足實(shí)際工程需求，說明本文計(jì)算策略適合和本文船型相似的船舶。

表6 KCS 船舶模型計(jì)算對比Tab. 6 Comparison of KCS ship model calculation

4 結(jié) 語

1）在船舶數(shù)值模擬過程中，數(shù)值計(jì)算精確度并不是隨著網(wǎng)格密度的增大而增大，但是隨著計(jì)算區(qū)域的增加，網(wǎng)格密度對計(jì)算結(jié)果的影響減小。可以適當(dāng)?shù)夭捎幂^大的計(jì)算區(qū)域和較小的網(wǎng)格密度，以節(jié)省計(jì)算時(shí)長。對于和本文相似的散貨船，推薦網(wǎng)格基礎(chǔ)尺寸為0.14，計(jì)算區(qū)域?yàn)? 倍船長，并采用SST K-Omega湍流模型。

2）采用經(jīng)驗(yàn)公式法計(jì)算和本文船型相似的船舶靜水阻力時(shí)，建議采用桑海公式可取得更高精度。

3）采用30 000 t 散貨船以及KCS 船，對本文提出的數(shù)值計(jì)算策略進(jìn)行驗(yàn)證分析，實(shí)驗(yàn)誤差均在5%以內(nèi)，滿足實(shí)際的工程需求。