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

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

0 引 言

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

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

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

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

1 數(shù)值計算

1.1 數(shù)值計算模型

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

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

1.2 控制方程

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

連續(xù)性方程：

動量方程：

1.3 邊界條件的選取

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

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

1.4 網(wǎng)格劃分

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

1.5 物理模型的選擇

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

式中：Vi為一個網(wǎng)格單元體中流體的體積分數(shù)；V為一個網(wǎng)格單元的體積。

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

圖2 船體水的體積分數(shù)Fig. 2 Volume fraction of water on the hull

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

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

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值計算方法驗證

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

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

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

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

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

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

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

桑海公式：

柏蘭特-許立汀公式：

休斯公式：

ITTC1957 公式：

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

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

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

圖6 阻力占比Fig. 6 Resistance ratio

3.3 KCS 船舶航速驗證

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

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

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

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

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

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

4 結(jié) 語

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

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

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