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

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

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

0 引 言

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

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

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

1 三體船型建模

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

表1 三體船模主尺度參數(shù)[10]

圖1 三體船圖

2 控制方程與流域設(shè)置

2.1 控制方程

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

(1)

張量形式的時均RANS方程為

(2)

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

(3)

2.2 坐標(biāo)系及計算流域設(shè)置

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

圖2 船體坐標(biāo)系

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

圖3 流域設(shè)置

3 網(wǎng)格因素的影響性分析

3.1 船體表面網(wǎng)格尺寸影響

研究表明，船體表面網(wǎng)格劃分尺寸范圍為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)的網(wǎng)格對船體表面進行劃分并采用STAR-CCM+平臺進行靜水阻力數(shù)值計算.表2為試驗及數(shù)值模擬各參數(shù)設(shè)定.表3為不同船體表面網(wǎng)格尺寸下的計算結(jié)果與試驗結(jié)果[12]對比.

表2 三體船模型各參數(shù)設(shè)定

表3 船體表面不同網(wǎng)格尺寸結(jié)果對比

由表3船體表面不同網(wǎng)格尺寸變化后的結(jié)果對比可知：

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

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

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

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

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

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

3.2 第一層網(wǎng)格厚度影響

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

圖5 第一層網(wǎng)格節(jié)點分布示意圖

通常，定義船體表面以外的第一層網(wǎng)格節(jié)點到船體表面的無因次距離為y+，則無因次系數(shù)y+為

(4)

式中：Δy為第一層網(wǎng)格節(jié)點到船體表面的距離；L為主體長度；Re為相對主體長度的雷諾數(shù).

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

表4 不同y+取值時結(jié)果對比

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

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

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

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

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

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

3.3 加密區(qū)域網(wǎng)格尺寸影響

三體船型周圍流場變化復(fù)雜，將船體周圍區(qū)域加密處理，能夠更加準(zhǔn)確地捕捉流場細(xì)節(jié).文獻[9]表明船體周圍加密區(qū)域網(wǎng)格尺寸取值為1.2%L～2%L時更為合理，因此，三體船周圍加密區(qū)域的網(wǎng)格分別劃分為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).計算時，船體表面網(wǎng)格尺寸取為35 mm，第一層網(wǎng)格節(jié)點到船體表面的無因次距離y+=170，計算結(jié)果見表5.

表5 加密區(qū)域不同網(wǎng)格尺寸結(jié)果

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

3.4 網(wǎng)格因素敏感性對比

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

圖7 阻力計算結(jié)果對網(wǎng)格因素尺度變化敏感性對比

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

4 湍流模型影響

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

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

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

表6 不同湍流模型結(jié)果對比

5 組合方案驗證

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

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

表7 阻力結(jié)果對比

圖8 總阻力結(jié)果對比

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

6 結(jié) 論

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

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

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