水面比降對船舶航行阻力影響的數(shù)值研究

冀 楠，楊 光，舒麟棹，錢志鵬，楊 春

(重慶交通大學(xué) 航運與船舶工程學(xué)院，重慶 400074)

0 引 言

船舶航行過程中經(jīng)常會遇到有水面比降的情況，例如：明渠河道、受洪水漲落或潮汐影響的河道(長江南京至鎮(zhèn)江河段[1])或者河底高程沿流向逐漸降低的河道水面均會產(chǎn)生比降；海洋中長波波面也可以認為是有比降水面；山區(qū)河流因河道高程和寬度多變, 水面比降現(xiàn)象更常見。船舶航行于帶有比降的水面時會產(chǎn)生不同于一般水面的水動力特性，導(dǎo)致船舶航行姿態(tài)改變，進而影響航速甚至操縱性能的惡化，容易造成事故。因此，為了提高船舶航行在有比降水面時的安全性，有必要對水面比降對船舶的航行阻力與浮態(tài)影響進行研究。

對水面比降影響船舶航行問題的研究主要有基于經(jīng)驗公式的理論分析方法和模型試驗法，研究范圍基本集中在內(nèi)河航道的上灘問題上?；诮?jīng)驗公式的理論分析方面，曹民雄等[2]利用急流灘模型的船模上灘試驗資料，對比分析了幾種內(nèi)河船舶上灘阻力計算方法；周曉杰等[3]根據(jù)航道規(guī)范和理論計算，對流速-比降組合值對應(yīng)的船舶上灘指標進行計算分析，得到了研究河段對應(yīng)的規(guī)范設(shè)計尺度；董曉韜[4]采用船舶推力-航行阻力平衡的方法，計算多種功率、載重條件下船舶自航上灘的臨界流速-比降組合，從而獲取多種工況下的急灘通航水力指標；張鵬等[5]結(jié)合船舶在洪水急流灘自航上灘的實測資料，確定了三峽庫區(qū)洪水急流灘河段船舶自航上灘的水力指標；李慶寧等[6]研究了三峽-葛洲壩航段航行的3種船型在泄洪期安全上水過灘的通航能力；童思陳等[7]以山區(qū)河流瀾滄江為代表，提出了適合于瀾滄江代表船舶的水流阻力計算式，并且結(jié)合坡降阻力計算, 建立了瀾滄江急流灘的通航水力指標。模型試驗方面，曹民雄等[8]進行了多種水力條件下的船模上灘試驗，發(fā)現(xiàn)水面比降、水流速度及船舶水深吃水比是船舶上灘的主要影響因素；劉曉菲等[9]開展了烏江小幺灘河段航道整治物理模型試驗，通過調(diào)整河段流速和比降，改善水流條件，以滿足船舶安全航行要求。

截至目前，國內(nèi)外學(xué)者對涉及船舶航行于比降水面的問題主要集中在評估水力特性和通航指標上，對航行阻力的計算采用的是經(jīng)驗公式或模型試驗。筆者采用計算流體力學(xué)(CFD)方法，基于商業(yè)軟件STAR-CCM+，以某型集裝箱船為研究對象，通過改變河道斷面縱向收縮角和入口流速獲得不同的水面比降工況，分析了水面比降對船舶航行阻力與浮態(tài)的影響，為船舶航行于帶有比降水面時的安全決策提供一定的參考。

1 船舶航行于比降水面阻力分析

船舶航行于比降水面時，由于水面坡度的影響，船舶的浮態(tài)及阻力都會區(qū)別于一般水面。如圖1，根據(jù)文獻[5]～文獻[7]的分析，船舶處于比降水面臨界狀態(tài)時，會產(chǎn)生縱傾，船舶重力沿船體縱向的分力被定義為比降阻力。此時船舶航行的總阻力和比降阻力的經(jīng)驗計算公式為：

R=Rv+Rj

(1)

式中：Rv為水流阻力，N；Rj為比降阻力，N。

比降阻力Rj的計算公式為：

Rj=βWJ

(2)

式中：W為船舶排水量，N；J為船長范圍內(nèi)平均水面比降；β為考慮船舶上灘時水面比降局部增大的修正系數(shù)，一般取β=1.1～1.2。

圖1 船舶航行于比降水面Fig. 1 The ship sailing on the water surface gradient

船舶航行時，由于水面存在一定比降，船首尾的水位差導(dǎo)致船舶傾斜。船傾斜后，浮心位置會發(fā)生偏移，致使浮心和重心不在同一鉛垂線上，于是作用方向相反的重力和浮力會促使船舶縱傾角增大，增大后的縱傾角設(shè)為θ′，θ′理論上要大于比降對應(yīng)的水面坡度角θ。因此，實際上船舶重力在船體縱向的分力應(yīng)該為Wsinθ′，這才是比降阻力。由此可知，采用式(2)計算得到的比降阻力會小于實際值，從而導(dǎo)致利用經(jīng)驗公式計算得到的船舶航行于比降水面時的總阻力會偏小。

2 數(shù)值方法

2.1 控制方程

在空間固定坐標系下，數(shù)值模擬采用的是不可壓縮流動RANS方程，控制方程為連續(xù)性方程和雷諾平均方程，如式(3)和式(4)：

(3)

(4)

2.2 湍流模型

考慮計算時間和計算精度，選用Realizablek-ε湍流模型，湍流動能k及湍流耗散率ε方程如式(5)和式(6)：

(5)

(6)

式中：k為湍流動能；ε為湍流耗散率；μt為湍流黏性系數(shù)；C1=1.44；C2=1.9;σk=1.0；σε=1.2。其中μt為：

(7)

式中：Cμ=0.09。

3 計算模型與驗證

3.1 計算模型

以國際研討會標準集裝箱船KCS為研究對象，如圖2。船模主要參數(shù)為：垂線間長Lpp=7.278 6 m，型寬B=1.019 m，型深D=0.341 8 m，設(shè)計吃水T=0.34 m，濕表面積S=9.438 m2，方形系數(shù)Cb=0.65，模型縮尺比λ=31.67。

KCS雖為海船船型，但其船型的長、寬、吃水與長江過閘船型集-18[10]的相應(yīng)參數(shù)基本呈2∶1的尺度比。雖然KCS船型的方形系數(shù)0.65比集-18的方向系數(shù)0.53大，但根據(jù)船舶原理[11]中關(guān)于方形系數(shù)對船舶阻力影響的敘述可知，在低于1.93 m/s的航速范圍內(nèi)，方形系數(shù)差異造成的船型阻力變化僅為10%左右。

圖2 KCS船幾何模型Fig. 2 KCS ship geometric model

3.2 船舶浮態(tài)與阻力之間關(guān)系的驗證

船舶在比降水面航行時會由于水面坡度而產(chǎn)生不同程度的尾傾，進而造成船舶阻力的變化。為了驗證文中數(shù)值模型的合理性，將船舶在不同初始縱傾角下的阻力數(shù)值計算結(jié)果與文獻[12]進行了對比。

表1 網(wǎng)格基本參數(shù)信息Table 1 Grid basic parameter information

表2 船模縱傾阻力數(shù)值計算與試驗對比Table 2 Comparison between numerical calculation and experiment of trim resistance of ship model

由表2可以發(fā)現(xiàn)，三套網(wǎng)格的計算結(jié)果與試驗的擬合度較好，最大誤差僅為3.21%。網(wǎng)格1在縱傾角小于0.5°后，與試驗結(jié)果存在一定偏離，特別是在縱傾角為0.25°和-0.5°時的誤差較大；網(wǎng)格3在縱傾角0.25°時誤差最大，大于0.25°時的計算結(jié)果較網(wǎng)格2偏差也較大。在綜合考慮計算精度與時間成本后，選擇網(wǎng)格2的屬性設(shè)置參數(shù)作為數(shù)值計算的網(wǎng)格劃分參數(shù)。

3.3 計算工況設(shè)置

計算水深設(shè)定為8 m，水深吃水比為h/T=23.5。通過改變河道斷面縱向收縮角和入口流速V來獲得不同的水面比降。河道斷面縱向收縮角是邊坡與x方向的夾角α，如圖3。筆者設(shè)置夾角α分別為1°、2°、3°、5°，入口流速V分別為0.5、1.25、2 m/s，以α-V來表示對應(yīng)工況。數(shù)值計算一共設(shè)置3組工況，第一組工況：1-0.5、2-0.5、3-0.5、5-0.5；第二組工況：1-1.25、2-1.25、3-1.25、5-1.25；第三組工況：1-2、2-2、3-2、5-2。

船體坐標系如圖4，原點O位于船模重心，x軸指向船首，y軸指向左側(cè)岸壁，z軸為垂直水面向上。船模首傾時縱傾角為正值，尾傾時縱傾角為負值。

圖3 計算域平面Fig. 3 Plan of calculation domain

圖4 船體坐標系Fig. 4 Ship hull coordinate system

3.4 邊界條件及網(wǎng)格劃分

假定船舶沿航道中線航行，且航道關(guān)于航道中線對稱，因此，航道和船舶均可采用半模進行設(shè)置和計算。邊界條件的設(shè)置如圖5，計算域的進口、頂部及底部均設(shè)置為速度進口，出口設(shè)置為壓力出口，側(cè)邊及船體表面設(shè)置為壁面，計算域中縱剖面設(shè)置為對稱平面。計算域中，進口距船首2Lpp，出口距船尾3Lpp，出口的寬度固定為0.7Lpp，邊坡長度固定為3.3Lpp。設(shè)定船模距岸壁的距離超過一倍船長，水深吃水比h/T=23.5，認為船舶不受岸壁和淺水效應(yīng)的影響，只在帶有比降水面的敞水水域航行。數(shù)值計算放開了船體升沉和縱傾兩個自由度；求解時間離散采用一階；擴散項采用二階中心差分格式離散；對流項采用二階迎風(fēng)格式，以提高計算精度；壓力-速度耦合方程采用SIMPLEC算法進行求解。

網(wǎng)格采用切割體網(wǎng)格，對自由液面及船體周圍進行體加密，如圖6。船體首尾型線曲度較大，需進行體加密，如圖7。經(jīng)過多次試算，最終確定無量綱化的壁面距離y+值的范圍在30～60之間，說明船體邊界層的網(wǎng)格設(shè)置是合理的。

圖5 計算域邊界條件Fig. 5 Boundary conditions of calculation domain

圖6 計算域整體網(wǎng)格Fig. 6 Global grid of calculation domain

圖7 船首尾網(wǎng)格Fig. 7 Bow and stern grid

4 結(jié)果分析

4.1 阻力分析

表3匯總了3組工況計算的阻力值，其中水流阻力作為對比項，是在提取到對應(yīng)比降工況下船體附近平均流速數(shù)值后，另外計算的船舶在敞水中的阻力值。比降阻力的數(shù)值模擬值由總阻力與水流阻力的差值得到，比降阻力的經(jīng)驗公式值由式(2)求得，并且也給出了經(jīng)驗公式值與數(shù)值模擬值的差值百分比。從表3可知：每組工況中，隨著比降的增加，船舶的總阻力逐漸增大，意味著航行難度增加；當(dāng)入口流速一定，每組工況的水流阻力相差并不大，水流阻力主要由船速決定，或者說是固定船模數(shù)值計算中的入口水流流速決定。在每組工況中，船舶周圍的流速與入口流速相比變化也不大，這說明式(1)將船舶在有比降水面上航行時的總阻力分解為水流阻力和比降阻力兩個部分的思路是合理的。

表3給出了3組工況的比降阻力數(shù)值模擬值和經(jīng)驗公式值隨水面比降變化的對比。由表3可以發(fā)現(xiàn)，采用經(jīng)驗公式計算的比降阻力明顯偏小，這與第1節(jié)理論分析得出的結(jié)論一致。而且，在每組工況中，隨著比降的增大，相對誤差又在縮小，說明式(2)在大比降時計算準確度高，而并不適用于小比降時的情況。究其原因，應(yīng)該是因為比降較小時，船舶重力浮力力矩造成的船舶縱傾角調(diào)整相對影響較大，而在大比降時，由于比降數(shù)值大，因此船舶重力浮力力矩造成的船舶縱傾角調(diào)整相對較小，比降阻力的數(shù)值趨向于直接用比降坡度角計算的結(jié)果。

圖8是3組工況的摩擦阻力Rf和壓阻力Rp隨不同比降變化的曲線。

圖8 不同比降下的摩擦阻力和壓阻力曲線Fig. 8 The frictional resistance and pressure resistance curves under different gradient

從圖8中可以看出，阻力的主要成分是壓阻力Rp(粘壓和興波阻力之和)，且壓阻力Rp會隨著水面比降的增大而明顯增大；摩擦阻力Rf隨水面比降的變化較小，在高流速工況下，甚至有一定的減小。摩擦阻力Rf與船體附近流速以及船模的濕表面積Sm有關(guān)。表4給出了3組工況船模在不同比降下的濕表面積Sm。由表4可以發(fā)現(xiàn)：每組工況內(nèi)，濕表面積Sm變化很小；工況之間，由于入口流速變大，船模濕表面積Sm有一定的增加，但增加數(shù)值仍然很小。由此可知，船舶在有比降水面航行時，船舶摩擦阻力Rf的變化主要由船體附近流速決定，而與濕表面積Sm的變化關(guān)系不大。

表4 不同比降下KCS船模濕表面積Table 4 The wet surface area of KCS ship model under different gradient

圖9展示了第一組工況船體底部動壓分布情況，從圖9中可知：比降為0.1‰ 時，低壓區(qū)覆蓋船底首部到尾部的全部區(qū)域；比降為0.2‰ 時，船首部壓力增大，船底低壓區(qū)縮減到僅覆蓋尾部區(qū)域；比降為0.4‰ 時，船底低壓區(qū)進一步縮減到僅覆蓋尾肩部；比降為0.6‰ 時，雖然流速幾乎沒變，但比降的增大仍然使得船底的壓力數(shù)值整體增大。由此可見，比降的增大，不但會增大船體首尾的壓力差，也會增大船體所受水壓力的數(shù)值。

圖9 船體底部壓力隨比降變化分布Fig. 9 Distribution of hull bottom pressure changing with gradient

4.2 自由液面波形和船舶浮態(tài)分析

圖10為第三組工況的自由液面波形等值線。圖11為3組工況距船中縱剖面1倍的船寬B處波高H沿船長X方向的分布。由圖10、圖11可知：每組工況下，隨著比降的增大，波面的平均傾角會增大。3組工況間的對比規(guī)律為：隨著船速的提高，船體興波對水面的擾動作用逐漸顯著，船首壅水導(dǎo)致首波峰的幅值逐漸增大，船尾后部的波浪起伏也逐漸明顯。從能量的角度看，波浪起伏越大，需要消耗越多的能量，這一部分能量最終以阻力的形式體現(xiàn)，這也是船舶航行于比降水面時阻力增大的原因之一。

圖10 自由液面波形等值線Fig. 10 Free surface waveform contour

圖12是3組工況下船體縱傾角隨水面比降變化的曲線，并且給出了水面坡度角作為對比，因為尾傾時船體縱傾角為負值，為了便于比較，圖12中的水面坡度角取負值。由圖12可發(fā)現(xiàn)，船舶在上灘航行時船身會產(chǎn)生尾傾現(xiàn)象。隨著水面比降的增大，船身尾傾更嚴重，并且尾傾角度大于水面坡度角，這也進一步證明采用式(2)計算出的比降阻力值會偏小。

圖13是不同流速下，船尾自由面隨水面比降變化分布，數(shù)值仿真中采用VOF方法捕捉船體自由液面。從圖13中可以明顯看到，隨著水面比降增大，自由液面高度在逐漸抬升，說明船體尾傾角度在增大。

圖11 不同比降下的波高分布Fig. 11 Wave elevation distribution under different gradient

圖12 縱傾隨比降的變化曲線Fig. 12 Variation curves of trim changing with gradient

圖13 船尾自由面隨比降變化分布Fig. 13 Distribution of free surface of stern changing with gradient

4.3 流場特性分析

槳盤處伴流場是計算螺旋槳推力及主機功率的重要參數(shù)，因此需研究水面比降對槳盤處伴流場的影響[14]。單槳船的槳盤面關(guān)于中縱剖面左右對稱，為方便對比，采用圖14進行分析，圖14中Wx為槳盤處標稱伴流分數(shù)。由圖14可知：圖14(a)中，槳盤面中心處的高伴流區(qū)隨著比降的增大呈現(xiàn)擴大的趨勢，槳盤面下部的伴流隨著比降的增大有向外擴展的趨勢，即隨著比降的增大，槳盤處的平均伴流增大；圖14(b)中，槳盤面下部的伴流隨著比降增大仍呈現(xiàn)向外擴展的趨勢，但上部水流卻隨著比降的增大呈收縮趨勢，說明高伴流區(qū)在向下移動，伴流的變化必將對螺旋槳的性能造成影響。

圖14 第三組工況下標稱伴流分數(shù)Wx隨比降變化的分布Fig. 14 Distribution of nominal wake fraction changing with gradient in the third group of working conditions

圖15～圖17分別展示了第三組工況船底、船首、船尾渦量細節(jié)分布，渦量顯示根據(jù)Q準則得到。從圖15中可以看到船首舭渦、側(cè)渦、肩渦和船尾舭渦、船尾渦等紊亂的渦流。對比4個工況下的船底渦量分布可以發(fā)現(xiàn)，隨著水面比降的增大，船底渦量強度會逐漸減弱。水面比降為1.7‰ 時，船底渦量強度最高，船首、尾舭渦附著在船底；水面比降增大至7.4‰ 時，舭渦渦量減少，強度也變?nèi)?。從圖16可以發(fā)現(xiàn)，水面比降為1.7‰ 時，船首側(cè)渦和肩渦清晰可見，并且呈現(xiàn)融合狀態(tài)；比降增大至3.3‰ 時，側(cè)渦、肩渦分離，且側(cè)渦強度已開始減弱；比降增大至4.6‰ 時，肩渦強度開始減弱；當(dāng)水面比降為7.4‰ 時，船首側(cè)渦分離成兩小段，兩側(cè)翼狀結(jié)構(gòu)消失。從圖17可以發(fā)現(xiàn)，水面比降增加時，船尾渦的強度也會減弱，水面比降增大到4.6‰ 時，船尾渦的尾翼狀結(jié)構(gòu)消失。

圖15 船底渦系分布Fig. 15 Diseribution of vontex system at bottem of hull

圖16 船首渦量Fig. 16 Bow vorticity

圖17 船尾渦量Fig. 17 Stern vorticity

5 結(jié) 論

筆者以實驗船模為研究對象，對船舶航行于有比降水面進行了數(shù)值模擬，分析了船舶的阻力、興波、浮態(tài)和黏性流場的變化。得出了以下主要結(jié)論：

1)計算流體力學(xué)(CFD)方法可以用來模擬船舶遭遇有比降水面航行時的阻力變化和航行特性，評估船舶航行能力。

2)船舶航行在帶有比降的水面時，水面比降越大，船舶所受到的總阻力也越大，增加了航行難度；水面比降對摩擦阻力的影響小，對壓阻力的影響較大。

3)船舶在有比降水面航行會產(chǎn)生一定量的尾傾。隨著水面比降的增大，船身尾傾會更嚴重，并且尾傾角度會大于比降對應(yīng)水面坡度角，這將導(dǎo)致比降阻力的真實結(jié)果比經(jīng)驗公式計算出的數(shù)值大。

4)船舶航行在帶有比降的水面時的流場特性表現(xiàn)為：船尾槳盤處的伴流分數(shù)會隨著水面比降的增大而增大，影響螺旋槳的運轉(zhuǎn)效率，增加航行難度；水面比降的變化會影響船體渦量強度和渦系的發(fā)展，渦量強度隨著水面比降的增大會減弱。