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

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

0 引 言

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

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

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

1 三體船型建模

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

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

圖1 三體船圖

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

2.1 控制方程

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

(1)

張量形式的時(shí)均RANS方程為

(2)

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

(3)

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

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

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

由于船體左右對(duì)稱，故采用單側(cè)模型進(jìn)行計(jì)算，所得計(jì)算結(jié)果最終處理為整船體的阻力數(shù)值.文中流域設(shè)置為一長(zhǎng)方體，由于試驗(yàn)中三體船周圍的流場(chǎng)變化相對(duì)復(fù)雜，故在船體周圍和水線面附近設(shè)置加密區(qū)域.利用STAR-CCM+平臺(tái)的自動(dòng)網(wǎng)格劃分功能，采用切割六面體網(wǎng)格，網(wǎng)格節(jié)點(diǎn)分布系數(shù)取r*≈1.2，網(wǎng)格層數(shù)為6層，實(shí)現(xiàn)網(wǎng)格自動(dòng)劃分.流域設(shè)置見圖3，其中，L表示主體長(zhǎng)度，下文流域設(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)格對(duì)船體表面進(jìn)行劃分并采用STAR-CCM+平臺(tái)進(jìn)行靜水阻力數(shù)值計(jì)算.表2為試驗(yàn)及數(shù)值模擬各參數(shù)設(shè)定.表3為不同船體表面網(wǎng)格尺寸下的計(jì)算結(jié)果與試驗(yàn)結(jié)果[12]對(duì)比.

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

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

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

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

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

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

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

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

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

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

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

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

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

(4)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 湍流模型影響

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

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

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

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

5 組合方案驗(yàn)證

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

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

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

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

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

6 結(jié) 論

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

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

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