陳昆鵬，張青山，杜云龍

(上海船舶運(yùn)輸科學(xué)研究所有限公司 航運(yùn)技術(shù)與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上海 200135)

0 引 言

目前，最常用的船舶有效功率預(yù)測(cè)方法有模型尺度的拖曳水池試驗(yàn)和模型尺度的數(shù)值模擬2種。然而，在模型尺度下通常只能滿足模型與實(shí)船的弗勞德數(shù)相同,不能滿足雷諾數(shù)相同，因此實(shí)尺度船舶數(shù)值模擬方法是目前國(guó)內(nèi)外學(xué)者研究的熱點(diǎn)[1]。CASTRO等[2]采用壁面函數(shù)和沙粒粗糙度對(duì)KCS集裝箱船進(jìn)行了實(shí)尺度自航數(shù)值模擬。DEMIREL等[3]提出了一種基于非穩(wěn)態(tài)雷諾平均納維-斯托克斯(Reynolds-Averaged Navier-Stokes,RANS)的計(jì)算流體動(dòng)力學(xué)(Computational Fluid Dynamics,CFD)模擬方法，該數(shù)值模擬方法能預(yù)測(cè)船用涂料和生物污底對(duì)船舶阻力的影響，并考慮粗糙度對(duì)實(shí)尺度KCS船數(shù)值模擬結(jié)果的影響，從而得出船舶的阻力和有效功率。張立等[4]對(duì)標(biāo)準(zhǔn)模型尺度船舶開展了CFD多維度仿真，得到了不同維度下的計(jì)算結(jié)果，復(fù)合航次試驗(yàn)結(jié)果表明，網(wǎng)格密度、無量綱壁面距離y+和湍流模型對(duì)計(jì)算結(jié)果有很大影響。邵文勃等[5]基于CFD技術(shù)對(duì)滑行艇的靜水阻力進(jìn)行了計(jì)算，分析了時(shí)間步長(zhǎng)和近壁面網(wǎng)格劃分對(duì)計(jì)算精度的影響，數(shù)值模擬結(jié)果表明，時(shí)間步長(zhǎng)和y+發(fā)生變化主要影響的是滑行艇的摩擦阻力和艇底表面的氣液分布，數(shù)值計(jì)算結(jié)果與模型試驗(yàn)結(jié)果能較好地吻合。胡海洲[6]開展了轂帽鰭和舵附推力鰭的實(shí)尺度節(jié)能效果研究，其y+設(shè)置在30～500范圍內(nèi)。NIKLAS等[7]將y+設(shè)置在30～200范圍內(nèi),對(duì)實(shí)尺度船舶開展了耐波性研究。

RANS方法通過壁面函數(shù)控制近壁面處的流體流動(dòng)，近壁面處第一層網(wǎng)格的高度用y+表示[4]。在隱式非定常模擬中，時(shí)間步長(zhǎng)往往取決于流動(dòng)特性，而與庫(kù)朗數(shù)無關(guān)。為得到合理的計(jì)算結(jié)果，應(yīng)根據(jù)模擬特征采用不同的時(shí)間步長(zhǎng)。文獻(xiàn)[1]提出,在實(shí)尺度船舶阻力性能模擬中，目前普遍采用的時(shí)間步長(zhǎng)范圍為0.005Lpp/U～0.010Lpp/U(Lpp為船舶垂線間長(zhǎng);U為船舶航速)。本文以某23萬載重噸散貨船為研究對(duì)象,根據(jù)其特點(diǎn)，結(jié)合船長(zhǎng)和航速選取時(shí)間步長(zhǎng)，建立不同的近壁面網(wǎng)格，以舯部為參考點(diǎn)讀取y+值，開展實(shí)尺度船舶阻力數(shù)值模擬研究，分析時(shí)間步長(zhǎng)和近壁面網(wǎng)格劃分對(duì)計(jì)算精度的影響。

1 數(shù)值方法

本文的控制方程為RANS方程和連續(xù)方程，滿足質(zhì)量守恒和動(dòng)量守恒。

1) 所有流動(dòng)問題都要滿足質(zhì)量守恒定律，該定律可表述為：?jiǎn)挝粫r(shí)間內(nèi)流體微元體中質(zhì)量的增加等于同一時(shí)間內(nèi)流入該微元體的凈質(zhì)量。根據(jù)該定律,可得出質(zhì)量守恒方程為

(1)

在笛卡爾坐標(biāo)系中，式(1)可表示為

(2)

對(duì)于不可壓縮流體，其密度為常數(shù),有

(3)

式(1)～式(3)中:ρ為流體密度;t為時(shí)間;v為速度矢量;u、v和w為速度矢量在x、y和z方向的分量。質(zhì)量守恒方程亦稱作連續(xù)方程。

2) 所有流動(dòng)問題都要滿足動(dòng)量守恒定律，該定律可表述為：微元體中流體的動(dòng)量對(duì)時(shí)間的變化率等于外界作用在該微元體上的各種力之和。動(dòng)量守恒方程(又稱納維-斯托克斯方程)為

(4)

(5)

(6)

式(4)～式(6):Su=Fx+sx;Sv=Fy+sy;Sw=Fz+sz;μ為流體的動(dòng)力黏度。Fx、Fy和Fz為微元體上的體力;對(duì)于黏性為常數(shù)的不可壓縮流體,sx、sy和sz為0。

雷諾平均方法又稱RANS方法,其核心是不直接求解瞬時(shí)的納維-斯托克斯方程，而是在時(shí)均化的方程中通過某種模型將瞬態(tài)的脈動(dòng)量體現(xiàn)出來。

關(guān)于湍流脈動(dòng)值的雷諾應(yīng)力項(xiàng)屬于新的變量，為使RANS方程封閉，本文采用兩方程渦黏模型可實(shí)現(xiàn)的k-ε(k為紊流脈動(dòng)動(dòng)能;ε為紊流脈動(dòng)動(dòng)量的耗散率)湍流模型。

采用VOF(Volume Of Fluid)方法捕捉自由液面，該方法通過研究網(wǎng)格單元中流體和網(wǎng)格體積分?jǐn)?shù)確定和追蹤自由液面。若水的體積分?jǐn)?shù)為1，則該網(wǎng)格單元全部充滿水；若水的體積分?jǐn)?shù)為0，則該網(wǎng)格單元為空氣；若水的體積分?jǐn)?shù)在(0,1)范圍內(nèi)，則該網(wǎng)格單元為交界面。

由于船體壁面附近存在較大的速度梯度，只有劃分高質(zhì)量的邊界層網(wǎng)格才能求解近壁面的流場(chǎng)信息。與模型尺度數(shù)值模擬相比，實(shí)尺度船舶數(shù)值模擬中的邊界層相對(duì)較薄，第一層邊界層對(duì)計(jì)算精度有重要影響，其厚度通常用y+表示。

y+的表達(dá)式為

y+=(yuτ)/ν

(7)

式(7)中:y為該點(diǎn)到壁面的距離;ν為流體的運(yùn)動(dòng)黏性系數(shù);uτ為壁面摩擦速度。[8]

2 計(jì)算域與網(wǎng)格劃分

數(shù)值模擬中的計(jì)算域?yàn)榱⒎襟w，與模型尺度不同，實(shí)尺度船舶的計(jì)算域相對(duì)較大，長(zhǎng)、寬和高分別約為3 330 m、860 m和1 210 m，船首距離邊界超過4倍船長(zhǎng)，船尾距離邊界超過5倍船長(zhǎng)，船體到側(cè)邊界的距離和水深均超過2.5倍船長(zhǎng)。由于船體是左右對(duì)稱的，可計(jì)算半船，船體中縱剖面所在面為對(duì)稱面。在計(jì)算域邊界中，頂部、底部和遠(yuǎn)離船體的側(cè)面為滑移壁面，入口設(shè)置為速度入口，出口設(shè)置為壓力出口。圖1為數(shù)值模擬中的計(jì)算域。

圖1 數(shù)值模擬中的計(jì)算域

網(wǎng)格采用非結(jié)構(gòu)化網(wǎng)格;船體附近的網(wǎng)格設(shè)置得較為密集;對(duì)船首、船尾和水線面的網(wǎng)格作加密處理；遠(yuǎn)離船體位置的網(wǎng)格設(shè)置得較為稀疏。通過建立矩形等加密區(qū)域控制體網(wǎng)格局部加密，對(duì)船體附近、艏艉和水線面網(wǎng)格進(jìn)行局部加密，可平衡計(jì)算資源與網(wǎng)格精細(xì)度。圖2為船體附近網(wǎng)格劃分情況。

a) 船體附近網(wǎng)格

b) 自由液面區(qū)域網(wǎng)格加密

表1 網(wǎng)格無關(guān)性驗(yàn)證結(jié)果

由表1可知,最密集的第4套網(wǎng)格與第3套網(wǎng)格的差距為-2.1%，因此選取第3套網(wǎng)格進(jìn)行數(shù)值模擬。

3 計(jì)算工況

在實(shí)尺度船舶阻力性能模擬中，目前普遍采用的時(shí)間步長(zhǎng)范圍為0.005Lpp/U～0.010Lpp/U[1]，當(dāng)該船的垂線間長(zhǎng)為314.6 m，航速為14.5 kn時(shí)，其時(shí)間步長(zhǎng)范圍為0.211～0.316 s。為研究時(shí)間步長(zhǎng)對(duì)實(shí)尺度數(shù)值模擬結(jié)果的影響，在設(shè)計(jì)航速下從小到大設(shè)置若干個(gè)時(shí)間步長(zhǎng)，其范圍包含文獻(xiàn)[1]中的0.005Lpp/U～0.010Lpp/U，本文的時(shí)間步長(zhǎng)最小為0.000 24Lpp/U(0.010 s)，最大為0.015 00Lpp/U(0.633 s)，表2為數(shù)值模擬中不同工況下的時(shí)間步長(zhǎng)。

為研究y+對(duì)實(shí)尺度船舶數(shù)值模擬結(jié)果的影響，在設(shè)計(jì)航速下從小到大設(shè)置若干個(gè)邊界層厚度和層數(shù)。船體不同位置的y+值不同，因此統(tǒng)一取平行中體舯部作為參考點(diǎn)。根據(jù)時(shí)間步長(zhǎng)的數(shù)值模擬結(jié)果，以時(shí)間步長(zhǎng)0.080 s為例開展y+對(duì)比計(jì)算。

表2 數(shù)值模擬中不同工況下的時(shí)間步長(zhǎng)

在對(duì)時(shí)間步長(zhǎng)進(jìn)行對(duì)比研究時(shí)，近壁面網(wǎng)格、y+和其他參數(shù)不變;在對(duì)y+進(jìn)行對(duì)比研究時(shí)，時(shí)間步長(zhǎng)和其他參數(shù)不變。表3為數(shù)值模擬中不同工況下的y+。

表3 數(shù)值模擬中不同工況下的y+

4 數(shù)值模擬結(jié)果分析

表4為通過水池模型試驗(yàn)預(yù)報(bào)的14.5 kn航速下的實(shí)船阻力Rt和有效功率Pe。試驗(yàn)在上海船舶運(yùn)輸科學(xué)研究所有限公司的船模拖曳水池內(nèi)進(jìn)行，船模長(zhǎng)7.3 m，水池長(zhǎng)192 m，寬10 m，水深4.2 m。根據(jù)以往的預(yù)報(bào)結(jié)果，該船模拖曳水池試驗(yàn)預(yù)報(bào)結(jié)果與實(shí)船試航結(jié)果的吻合性較好，因此通過模型試驗(yàn)預(yù)報(bào)的實(shí)船阻力和有效功率可作為實(shí)尺度船舶數(shù)值模擬對(duì)比驗(yàn)證的依據(jù)。

表4 通過水池模型試驗(yàn)預(yù)報(bào)的14.5 kn航速下的實(shí)船阻力Rt和有效功率Pe

當(dāng)控制變量為時(shí)間步長(zhǎng)時(shí)，其對(duì)數(shù)值模擬結(jié)果的影響見表5。

表5 時(shí)間步長(zhǎng)對(duì)數(shù)值模擬結(jié)果的影響

圖3為總阻力隨時(shí)間步長(zhǎng)的變化趨勢(shì)。從圖3和表5中可看出:當(dāng)時(shí)間步長(zhǎng)小于0.211 s(即小于0.005 00Lpp/U)時(shí)，阻力和有效功率隨時(shí)間步長(zhǎng)的變化在1.5%以內(nèi);當(dāng)時(shí)間步長(zhǎng)大于0.211 s時(shí),阻力和有效功率隨時(shí)間步長(zhǎng)的增加大幅度上升;當(dāng)時(shí)間步長(zhǎng)為0.422 s(即0.010 00Lpp/U)時(shí),通過數(shù)值模擬得到的阻力相比作為基準(zhǔn)的時(shí)間步長(zhǎng)為0.010 s(即0.000 24Lpp/U)時(shí)的阻力數(shù)值模擬結(jié)果增加7.7%。

圖3 總阻力隨時(shí)間步長(zhǎng)的變化趨勢(shì)

當(dāng)時(shí)間步長(zhǎng)為0.010～0.211 s(即0.000 24Lpp/U～0.005 00Lpp/U)時(shí),宏觀總阻力變化幅度有限，但通過數(shù)值模擬得到的流場(chǎng)細(xì)節(jié)隨時(shí)間步長(zhǎng)的增加而有損失。圖4為自由液面隨時(shí)間步長(zhǎng)的變化趨勢(shì)，顏色表示自由液面高度在大地坐標(biāo)系下的變化。在網(wǎng)格和其他參數(shù)完全相同的情況下，當(dāng)時(shí)間步長(zhǎng)為0.211 s(即0.005 00Lpp/U)和0.316 s(即0.007 50Lpp/U)時(shí),很多波形細(xì)節(jié)完全消失，且自由液面出現(xiàn)不真實(shí)的異常鋸齒狀。

a) 時(shí)間步長(zhǎng)為0.010 s

b) 時(shí)間步長(zhǎng)為0.080 s

c) 時(shí)間步長(zhǎng)為0.211 s

d) 時(shí)間步長(zhǎng)為0.316 s

在網(wǎng)格和計(jì)算機(jī)資源相同的條件下，時(shí)間步長(zhǎng)與數(shù)值模擬計(jì)算的時(shí)間成本呈線性關(guān)系，時(shí)間步長(zhǎng)越小,數(shù)值計(jì)算的時(shí)間成本越高。綜合考慮計(jì)算速度和精度，認(rèn)為類似的肥大型船的時(shí)間步長(zhǎng)應(yīng)取0.002Lpp/U～0.005Lpp/U,對(duì)應(yīng)本文所述船舶0.08～0.21 s，時(shí)間步長(zhǎng)0.005Lpp/U～0.010Lpp/U對(duì)于本船來說過大。

表6為y+對(duì)總阻力數(shù)值模擬結(jié)果的影響;圖5為總阻力隨y+的變化趨勢(shì)。由表6和圖5可知:當(dāng)y+的值在88～679范圍內(nèi)時(shí),通過數(shù)值模擬得到的阻力和有效功率隨y+的增加而增大;當(dāng)y+值增加到679時(shí),阻力值相比y+=175時(shí)的阻力值增加2.3%;當(dāng)y+從175增大到415時(shí),阻力值增大2%;當(dāng)y+從175減小到88時(shí),阻力值減小2.5%。

此次數(shù)值模擬的對(duì)象為實(shí)船，其雷諾數(shù)遠(yuǎn)大于模型，對(duì)數(shù)律區(qū)域較大，因此沒有必要將y+控制在50～60，對(duì)于本文所述船舶和類似船型而言，y+在175～679范圍內(nèi)取何值理論上都是合理的，當(dāng)舯部y+=679時(shí),艏部會(huì)出現(xiàn)y+>679的區(qū)域，因此認(rèn)為該船通過近壁面網(wǎng)格設(shè)置的舯部y+應(yīng)為約415。

表6 y+對(duì)總阻力數(shù)值模擬結(jié)果的影響

該船在y+=415，時(shí)間步長(zhǎng)為0.08 s時(shí)的數(shù)值計(jì)算結(jié)果與模型預(yù)報(bào)結(jié)果相差3%，有較好的流場(chǎng)細(xì)節(jié)。

圖6為舯部參考點(diǎn)y+=259時(shí)的全船y+分布，艏部和艉部因速度變化較大,產(chǎn)生了1個(gè)艏部y+大于參考點(diǎn)y+的區(qū)域和1個(gè)艉部y+小于參考點(diǎn)y+的區(qū)域，船體大部分區(qū)域的y+都與參考點(diǎn)y+接近。

當(dāng)時(shí)間步長(zhǎng)為0.010～0.211 s，y+在175～305范圍內(nèi)時(shí)，相比時(shí)間步長(zhǎng)，實(shí)船數(shù)值模擬中的總阻力和有效功率對(duì)y+更敏感，時(shí)間步長(zhǎng)變化21倍對(duì)總阻力的影響與y+變化1.74倍相近。當(dāng)時(shí)間步長(zhǎng)在0.211～0.633 s范圍內(nèi)時(shí)，實(shí)船數(shù)值模擬中的總阻力對(duì)時(shí)間步長(zhǎng)更敏感，總阻力隨時(shí)間步長(zhǎng)的最大變化幅度為9.8%。

5 結(jié) 語(yǔ)

本文以RANS方程和連續(xù)方程為控制方程，采用有限體積法進(jìn)行空間離散，在控制時(shí)間步長(zhǎng)和y+等參數(shù)的情況下建立了某23萬載重噸散貨船的實(shí)尺度船舶阻力數(shù)值模擬計(jì)算模型，重點(diǎn)研究了時(shí)間步長(zhǎng)和y+值對(duì)計(jì)算精度的影響。結(jié)果顯示：在較大范圍內(nèi),時(shí)間步長(zhǎng)和y+對(duì)阻力計(jì)算結(jié)果的影響顯著，最大達(dá)9.8%，在一定范圍內(nèi)影響不超過2%;當(dāng)舯部y+為200～420，時(shí)間步長(zhǎng)取0.002Lpp/U～0.005Lpp/U時(shí)，合理設(shè)置時(shí)間步長(zhǎng)和y+可相對(duì)準(zhǔn)確地開展實(shí)尺度數(shù)值模擬，預(yù)報(bào)船舶的靜水阻力性能和有效功率。研究成果可供類似船舶開展實(shí)尺度船舶快速性數(shù)值模擬研究參考。