網格因素及湍流模型對三體船阻力計算影響分析

張明霞 韓兵兵 盧鵬程 趙正彬

(大連理工大學船舶工程學院 大連 116024)

0 引 言

近年來，三體船優良的水動力性能成為研究的熱點[1-2].利用計算流體力學研究船舶阻力預測、耐波性分析方面也取得成功應用.在進行三體船CFD數值模擬過程中，作為數值模擬和分析的載體，網格劃分的好壞直接影響到計算的可行性、收斂性、計算精度及計算效率，甚至關乎計算的成敗.

國內外學者深入開展了不同船型的CFD不確定度分析.文獻[3-4]研究了肥大型船伴流場數值模擬的網格劃分方法，并對影響該船型尾部伴流場的關鍵網格劃分進行了比較與分析，得出適用于肥大型船舶的伴流場網格劃分方法.孫偉華等[5]探究了不同網格因素對高速單體滑行艇阻力計算精度、收斂速度及結果穩定性的影響，得到適合滑行艇阻力計算的網格參數.對于三體船型的網格不確定度研究，同樣成為CFD領域一個重要課題.Caponnetto等[6-7]利用CFD軟件模擬了高速三體滑行艇周圍粘性流場，尋求獲得阻力計算結果誤差較小的網格劃分手段.陳康等[8-9]研究了網格因素對wigley三體數學船型阻力計算結果的影響并提出了該船型阻力改進方法.文獻[4]基于排水型三體船中體的網格劃分及其不確定度研究，得到一套可行的網格劃分方法.

在目前的船型CFD不確定度研究中，針對排水型三體船數值模擬過程中網格劃分方法研究的可參考文獻資料很少.不同于其他船型，三體船型由于片體間興波干擾導致周圍流場變化更為復雜，作為數值計算的重要部分，網格劃分不僅是對于計算域本身的離散化，而且對流場細節變化的捕捉和數值計算精度及計算成本會產生重要影響.基于STAR-CCM+平臺針對排水型三體船周圍流體流動特點，探討了不同網格因素及湍流模型選擇對阻力計算結果的影響，提出適用于排水航行的三體船型的網格劃分方法，并通過模型試驗數據進行驗證.

1 三體船型建模

文中利用CATIA平臺，按大連理工大學船池三體船實物模型1∶1船體建模作為數值計算對象，該模型主尺度參數見表1，型線圖、實物圖、CATIA三維視圖見圖1.

表1 三體船模主尺度參數[10]

圖1 三體船圖

2 控制方程與流域設置

2.1 控制方程

對于不可壓縮的黏性流體流動連續性方程為

(1)

張量形式的時均RANS方程為

(2)

三體船作為排水型船舶，需要考慮自由表面問題，文中采用VOF方法捕捉自由表面.VOF法通過定義一個流域體積函數F，來定義劃分的每個網格單元的狀態，F的值等于一個單元內流體體積與該單元體積之比.若F值等于1，則說明該單元全部為指定相流體所占據；若F值等于0，則該單元為無指定相流體單元，若F值介于0～1，說明該單元內含有自由表面[11].流域體積函數F的運算方程為

(3)

2.2 坐標系及計算流域設置

船體坐標系見圖2，X軸沿主體船長方向并指向主體艏部方向為正，Y軸沿船寬方向并指向主體左舷為正.側體中心線與主體中心線間的橫向跨距用a表示，a始終為正值；側體船舯與主體船舯的縱向距離用b表示，當側體在主體船舯之前時，b為正值，當側體在主體船舯之后時，b為負值.

圖2 船體坐標系

由于船體左右對稱，故采用單側模型進行計算，所得計算結果最終處理為整船體的阻力數值.文中流域設置為一長方體，由于試驗中三體船周圍的流場變化相對復雜，故在船體周圍和水線面附近設置加密區域.利用STAR-CCM+平臺的自動網格劃分功能，采用切割六面體網格，網格節點分布系數取r*≈1.2，網格層數為6層，實現網格自動劃分.流域設置見圖3，其中，L表示主體長度，下文流域設置與此相同.

圖3 流域設置

3 網格因素的影響性分析

3.1 船體表面網格尺寸影響

研究表明，船體表面網格劃分尺寸范圍為0.07%L～1%L之間，故分別采用5,15,25,35,45,55 mm(0.125%L,0.375%L,0.625%L,0.875%L,1.125%L,1.375%L)的網格對船體表面進行劃分并采用STAR-CCM+平臺進行靜水阻力數值計算.表2為試驗及數值模擬各參數設定.表3為不同船體表面網格尺寸下的計算結果與試驗結果[12]對比.

表2 三體船模型各參數設定

表3 船體表面不同網格尺寸結果對比

由表3船體表面不同網格尺寸變化后的結果對比可知：

1) 網格尺寸取5 mm(0.125%L)時，計算結果雖與試驗結果接近，但網格數量急劇增加，導致計算成本較高.

2) 船體表面網格尺寸大小為15，25，35 mm時，網格數量變化幅度平緩，計算結果誤差范圍在5%以內，計算結果與試驗結果較為接近.

3) 船體表面網格尺寸為45，55 mm時，計算結果相對誤差大于5%.分析原因：船體表面網格劃分粗糙，貼體性較差，對船體模型特征難以進行細致表達.

考慮不同船體表面網格尺寸下的阻力計算結果收斂速度及時間成本隨物理時間的變化，見圖4.

圖4 不同網格尺寸下阻力收斂性隨物理時間變化

由圖4可知，收斂后的阻力計算值明顯小于試驗值，但船體表面網格尺寸為5～35 mm時，試驗值與試驗值更為接近，其中，船體表面網格尺寸為35 mm時，計算精度最高.此外，不同網格尺寸下的結果收斂速度和穩定性無明顯差異.

3.2 第一層網格厚度影響

船舶CFD數值模擬通常采用壁面函數法[13]近似處理船體表面附近流動速度和湍流動能的梯度變化.通過在固壁附近物理梯度大的地方分布大量網格節點，來捕捉船體壁面附近流體流動的物理特性和流場細節.緊貼船體近固壁區域的流體可以分為黏性底層、過渡層和對數律層[14].當劃分網格的第一層網格節點位于過渡層內時(見圖5)，能夠有利于捕捉船體表面附近流場中的物理量及其梯度，所以第一層網格節點的高度是影響計算結果的重要因素.

圖5 第一層網格節點分布示意圖

通常，定義船體表面以外的第一層網格節點到船體表面的無因次距離為y+，則無因次系數y+為

(4)

式中：Δy為第一層網格節點到船體表面的距離；L為主體長度；Re為相對主體長度的雷諾數.

研究表明，y+值應該滿足11≤y+≤500.針對文中船型，在采用船體表面網格尺寸為35 mm的基礎上，劃分y+值分別為20，70，120，170，220，270六組網格進行探究,結果對比見表4.

表4 不同y+取值時結果對比

由表4可知，六組網格計算結果整體相對誤差范圍在1.94%～12.63%，具體分析如下.

1) 當20≤y+≤70時，計算結果與試驗結果誤差范圍為6.90%～12.63%，分析原因：第一層網格節點高度較小，節點落入粘性底層內部，分子粘性效應高于湍流效應導致計算誤差偏大；當y+=270時，計算結果與試驗結果的偏差為7.70%，是由于第一層網格節點接近對數律層，使得對船體表面周圍復雜流場中的物理量的捕捉度變低.

2) 當120≤y+≤220時，第一層網格節點位于過渡層附近，粘性效應與湍流效應相對均衡，節點分布有利于高效捕捉船體表面周圍流場中的物理量.此時，數值計算結果與試驗結果相對誤差在1.94%～2.74%.

對于不同取值的第一層網格節點到船體表面的無因次距離，考慮結果收斂速度、收斂時間及穩定性，四種y+值下阻力收斂性隨物理時間變化見圖6.

圖6 不同y+時阻力收斂性隨物理時間變化

由圖6可知，當y+取值為120～220時，即Δy取值為1.28～2.35 mm時，計算結果更為精確，且結果收斂時間、穩定性相近.結合表4結果對比分析，能夠確定當y+=170時，計算精度最高，且計算成本較低.

3.3 加密區域網格尺寸影響

三體船型周圍流場變化復雜，將船體周圍區域加密處理，能夠更加準確地捕捉流場細節.文獻[9]表明船體周圍加密區域網格尺寸取值為1.2%L～2%L時更為合理，因此，三體船周圍加密區域的網格分別劃分為50，55，60，65，70，75 mm(1.25%L，1.375%L，1.5%L，1.625%L，1.75%L，1.875%L).計算時，船體表面網格尺寸取為35 mm，第一層網格節點到船體表面的無因次距離y+=170，計算結果見表5.

表5 加密區域不同網格尺寸結果

由表5可知，船體周圍加密區域網格尺寸取值為50～60 mm(1.25%L～1.500%L)時，計算值與試驗值更為接近，誤差范圍在5%以內.但實際計算時，船體周圍加密區域網格尺寸為50 mm，生成網格數量急劇增加，導致計算成本增加.因此，對于文中三體船模的船體周圍加密區網格尺寸取值為1.375%L～1.500%L時計算成本較低，結果也更為精確.

3.4 網格因素敏感性對比

考慮到阻力計算結果對不同網格因素尺度變化的敏感程度，繪制了阻力計算結果對網格因素尺度變化敏感性對比曲線，見圖7.

圖7 阻力計算結果對網格因素尺度變化敏感性對比

由圖7中曲線緩急程度可知，不同網格因素對阻力計算結果的影響程度不同.其中，第一層網格節點到船體表面的無因次距離y+對結果的影響最為顯著，結果對船體周圍加密區域網格尺度變化及船體表面網格尺寸變化的敏感性依次降低.因此，在進行三體船CFD數值計算時應重點確定第一層網格節點的落點位置.

4 湍流模型影響

對于湍流問題的求解，不同的湍流模型的求解適用范圍有所不同.RANS湍流模型中，Spalart-Allmaras模型作為單方程模型，只需求解湍流粘性的輸運方程，無需求解當地剪切層厚度的長度尺度，對于流動尺度變化明顯的問題適用性較低.同時，Spalart-Allmaras模型中的輸運變量在近壁處的梯度要比k-ε中的小，因此，對網格粗糙導致的數值誤差敏感性低.k-ε模型是目前工程中應用最為廣泛的湍流模型，并在此基礎上發展衍生出一系列的湍流模型(RNG，Realizable)，該模型更加符合湍流的物理特性，適用于多尺度湍流擴散、旋轉流動、流動分離問題，但在逆壓梯度強度超限時精度降低.k-ω模型相比于k-ε模型，最明顯的變化就是關于耗散率ε的輸運方程被關于比耗散率ω的輸運方程代替，改進后的SSTk-ω模型能夠更好地捕捉粘性底層的流動，但對入口湍流參數過于敏感[15-16].

對于三體船周圍粘性流場模擬，基于上述網格因素結論，即船體表面網格尺寸取為35 mm；第一層網格節點到船體表面的無因次距離y+取值為170；船體周圍加密區域網格尺寸取值為60 mm.基于STAR-CCM+平臺中，分別采用k-ε、k-ω、Spalart-Allmaras三種湍流模型對其進行計算對比分析，計算結果見表6.

通過三種不同湍流模型下，阻力計算值與試驗值的對比發現，在網格劃分較為細致時，二方程湍流模型(k-ε，k-ω)計算結果顯著優于單方程湍流模型(Spalart-Allmaras)的計算結果.其中，k-ε湍流模型的適用度更高，其計算誤差為2.32%.

表6 不同湍流模型結果對比

5 組合方案驗證

通過上述分析，針對文中排水型三體船，可以得出網格因素及湍流模型選擇的優化組合，即船體表面網格尺寸為15～35 mm、第一層網格節點到船體表面的無因次距離y+=120～220、船體周圍加密區域網格尺寸為55～60 mm、湍流模型選取k-ε模型.

為進一步驗證該組合方案的準確性及更好地與試驗結果對比，設置與文獻[15]相同的系列航速及側體布局方案(側體橫向跨距a=700 mm，縱向距離b=-700 mm).計算靜水中，系列航速下的阻力值，并與試驗值進行對比，見表7和圖8.

表7 阻力結果對比

圖8 總阻力結果對比

由表7和圖8可知，該組合方案下阻力的相對誤差范圍在5.86%～8.62%，滿足計算精度要求.從二者曲線的對比來看，隨著航速的增加，阻力試驗值要略高于計算值，分析原因：試驗過程中隨著航速增加，噴濺阻力的影響增加導致總阻力增大，而數值計算中，對噴濺阻力影響尚未考慮.

6 結 論

1) 對于排水航行的三體船型，建議船體表面網格尺寸劃分取值為0.375%L～0.875%L；第一層網格節點到船體表面的無因次距離y+為120～220；船體周圍加密區域的網格尺寸取值為1.375%L～1.500%L.

2) 不同網格因素對阻力計算結果影響的重要程度不同，第一層網格節點厚度因素對結果影響最為顯著，船體周圍加密區域網格尺度因素和船體表面網格尺度因素對結果影響的敏感性逐漸降低.

3) 對于三體船周圍復雜的流場變化，采用k-ε湍流模型作為求解模型，結果誤差更小.但由于不完善的湍流模型理論導致的模擬模型誤差存在，湍流封閉模型理論缺陷仍須進一步研究.