陳昆鵬，張青山，杜云龍

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

0 引 言

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

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

1 數(shù)值方法

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

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

(1)

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

(2)

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

(3)

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

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

(4)

(5)

(6)

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

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

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

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

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

y+的表達式為

y+=(yuτ)/ν

(7)

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

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

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

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

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

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

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

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

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

3 計算工況

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

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

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

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

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

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

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

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

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

表5 時間步長對數(shù)值模擬結(jié)果的影響

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

圖3 總阻力隨時間步長的變化趨勢

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

a) 時間步長為0.010 s

b) 時間步長為0.080 s

c) 時間步長為0.211 s

d) 時間步長為0.316 s

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

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

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

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

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

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

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

5 結(jié) 語

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