基于STAR-CCM+的30 000 t 散貨船水阻力計算

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

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

0 引 言

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

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

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

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

1 數值計算

1.1 數值計算模型

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

表1 某30 000 t 船舶參數Tab. 1 Parameters of 30 000 t ship

1.2 控制方程

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

連續性方程：

動量方程：

1.3 邊界條件的選取

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

圖1 數值計算區域劃分Fig. 1 Division of numerical calculation area

1.4 網格劃分

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

1.5 物理模型的選擇

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

式中：Vi為一個網格單元體中流體的體積分數；V為一個網格單元的體積。

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

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

2 數值計算的影響因素分析

2.1 計算區域大小對計算結果的分析

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

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

表2 計算區域對數值計算結果的影響Tab. 2 Influence of calculation area on numerical results

2.2 網格密度對計算結果的影響

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

由表3 可知，在基礎尺寸是0.12 和0.14 時，6 倍船長的誤差均小于4 倍船長。基礎尺寸為0.12 的4 倍船長誤差和基礎尺寸為0.14 的6 倍船長相差不大，因此減少基礎尺寸，適當的增加計算區域，可以在保證計算精度的情況下，減少計算時長。計算區域6 倍船長，基礎尺寸為0.12 時，誤差最小，達到-0.62%，但是基礎尺寸為0.16 時的誤差達到1.39%，也可以滿足實際工程，且整體耗時也較短。

表3 網格密度對數值計算結果的影響Tab. 3 Influence of grid density on numerical results

2.3 湍流模型對計算結果的影響

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

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

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

3 數值計算方法驗證

3.1 數值計算結果可視化分析

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

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

圖4 數值計算結果對比Fig. 4 Comparison of numerical calculation results

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

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

船舶靜水阻力可以分為摩擦阻力和粘壓阻力，根據相關的平板拖曳試驗結果，摩擦阻力系數計算有以下幾種常用的經驗公式。

桑海公式：

柏蘭特-許立汀公式：

休斯公式：

ITTC1957 公式：

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

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

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

圖6 阻力占比Fig. 6 Resistance ratio

3.3 KCS 船舶航速驗證

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

表5 KCS 船舶參數Tab. 5 KCS ship parameters

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

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

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

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

4 結 語

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

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

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