船舶危險(xiǎn)對(duì)遇時(shí)水動(dòng)力數(shù)值研究

邱云明，陸冬青，鄧 銳

(1. 陸軍軍事交通學(xué)院鎮(zhèn)江校區(qū)，江蘇 鎮(zhèn)江 212003；2. 中山大學(xué)海洋工程與技術(shù)學(xué)院，廣東 珠海 519000)

0 引 言

近年來隨著航運(yùn)業(yè)的快速發(fā)展，船舶交通密度增大，兩船對(duì)遇態(tài)勢也增多，船舶間因水動(dòng)力相互干擾產(chǎn)生船吸現(xiàn)象，特別在兩船右舷對(duì)右舷相遇時(shí)，加上船舶避讓措施不當(dāng)，經(jīng)常發(fā)生碰撞事故。因此，對(duì)兩船對(duì)遇問題的研究引起了廣大學(xué)者的關(guān)注。目前，對(duì)于兩船會(huì)遇時(shí)水動(dòng)力性能研究大都采用數(shù)值計(jì)算和模型試驗(yàn)方法。在基于模型試驗(yàn)研究中，國外學(xué)者所做的工作較多。例如，Remery等[1]開展了航行船舶誘導(dǎo)的系泊力的相關(guān)試驗(yàn)，討論了兩船之間的橫向間距以及航行船舶的速度等有關(guān)因素對(duì)系泊船舶的水動(dòng)力性能的影響。Dand等[2]開展了航行船與停泊船的模型試驗(yàn)，并在此基礎(chǔ)上，推出船舶相互作用引起的力和力矩的經(jīng)驗(yàn)公式。Simmons等[3]研究了在岸壁效應(yīng)下兩船船間相互作用力干擾試驗(yàn)。Duffy等[4-5]基于模型試驗(yàn)研究了航行船經(jīng)過系泊船時(shí)，航道的不同因素對(duì)兩船之間相互作用力的干擾。在基于數(shù)值計(jì)算研究中，關(guān)于兩船會(huì)遇時(shí)非定常三維粘性流場的分析，更多學(xué)者采用勢流理論進(jìn)行了深入探究。Yeung等[6]基于細(xì)長體理論的匹配漸進(jìn)展開法，對(duì)開闊水域以及深淺水中兩船會(huì)遇時(shí)水動(dòng)力以及耐波性的變化等做了有關(guān)研究。鄭才土等[7]借助奇點(diǎn)疊加法相關(guān)理論，討論了在淺窄航道中兩船會(huì)遇時(shí)兩船性能彼此受到干擾。Gourla等[8]采用勢流理論，計(jì)算了兩船會(huì)遇和追越態(tài)勢下的兩船發(fā)生橫傾、縱傾等變化。

目前，隨著CFD技術(shù)迅速發(fā)展與應(yīng)用，研究兩船會(huì)遇時(shí)水動(dòng)力性能已經(jīng)考慮粘性流的存在。李學(xué)東等[9]討論了在狹窄水道內(nèi)兩船追越和對(duì)遇態(tài)勢下船間水動(dòng)力作用的問題。張晨曦等[10]基于數(shù)值模擬方法，深入討論兩船在淺水域中，追越船平行追越和會(huì)遇被追越船時(shí)的水動(dòng)力性能變化規(guī)律。鄒早建等[11]基于CFD流體軟件，通過求解RANS方程，重點(diǎn)分析了淺水中兩船會(huì)遇時(shí)水動(dòng)力的相互作用問題。

綜上所述，對(duì)于兩船會(huì)遇問題研究，已引起眾多學(xué)者的關(guān)注，并且開展了模型試驗(yàn)和數(shù)值模擬方法研究，取得了一些成果。本文為了精準(zhǔn)研究船舶危險(xiǎn)對(duì)遇時(shí)兩船間水動(dòng)力性能情況，采用動(dòng)網(wǎng)格及滑動(dòng)交界面等技術(shù)，并借助Star ccm+流體軟件研究右舷對(duì)右舷危險(xiǎn)對(duì)遇時(shí)的水動(dòng)力情況。

1 控制方程和湍流模型

1.1 控制方程

船舶危險(xiǎn)對(duì)遇時(shí)的水動(dòng)力的數(shù)值計(jì)算是通過求解三維粘性不可壓RANS方程，具體的方程如下：

連續(xù)方程為

式中： ui， uj為速度分量；μ為動(dòng)力粘度，ρ為粘度系數(shù)；為Reynolds應(yīng)力項(xiàng)。由于其存在使得方程組中的未知數(shù)大于已有方程數(shù)，需要對(duì)Reynolds應(yīng)力項(xiàng)做出某些假設(shè)，并建立相應(yīng)的應(yīng)力表達(dá)式封閉方程組。

1.2 湍流模型

在數(shù)值計(jì)算中，湍流模型選擇Realizable k -ε模型，此模型可以避免標(biāo)準(zhǔn) k-ε模型在處理時(shí)均應(yīng)變率因特別大而出現(xiàn)正應(yīng)力為負(fù)的問題，其中Realizable k-ε模型中湍動(dòng)能k輸運(yùn)方程如式（3），湍流耗散率ε輸運(yùn)方程如式（4）。

式中： ?k= 1.0； ?ε= 1.2；C2=1.9。

Realizable k -ε 模型在計(jì)算時(shí)引入的 ωk（旋轉(zhuǎn)角速度）和(流體時(shí)均轉(zhuǎn)動(dòng)速率張量)考慮到了流體旋轉(zhuǎn)的相關(guān)情況，包含流耗散率 ε輸運(yùn)方程中的項(xiàng)，在k為0時(shí)也不具奇異性，這是與標(biāo)準(zhǔn) k-ε 型和RNG k -ε模型之間存在的差異并結(jié)合以前學(xué)者的研究經(jīng)驗(yàn)，本文采用Realizable k -ε模型模擬兩船相遇時(shí)的水動(dòng)力性能。

2 數(shù)值模型

2.1 坐標(biāo)系統(tǒng)

研究中為了簡化計(jì)算模型，采用相對(duì)運(yùn)動(dòng)方法，設(shè)定一船相對(duì)靜止，進(jìn)行研究另一船右舷對(duì)右舷通過相對(duì)靜止船時(shí)的水動(dòng)力性能影響。如圖1所示，SPth表示兩船的橫向距離，ST表示兩船的縱向間距,船與兩岸間距分別表示為Sp1和Sp2，航道寬度W=Sp1+Sp+Sp2，水深h。在船舶三維坐標(biāo)中，采用空間固定右手坐標(biāo)系統(tǒng)o-xyz，x軸指向航行船運(yùn)動(dòng)方向，y軸指向航行船左舷，z軸指向水面，空間固定坐標(biāo)系原點(diǎn)位于靜止船中心處。

圖 1 計(jì)算坐標(biāo)系統(tǒng)Fig. 1 Calculates the coordinate system

2.2 計(jì)算對(duì)象

目前，在兩船相遇問題研究中，Duffy等[4]為了研究在航道寬度、水深等多種因素存在時(shí)，靜止船受到航行船經(jīng)過時(shí)的受力影響情況，做了船模試驗(yàn)，分析了停泊船所受力以及力矩等無量綱值。本文研究中，以其中的一種試驗(yàn)工況為參考，停泊船相當(dāng)于靜止船，進(jìn)行相關(guān)的數(shù)值模擬，Duffy等研究的試驗(yàn)中航道幾何尺寸如圖2所示。

圖 2 航道的橫切面幾何Fig. 2 Channel cross-section geometry

試驗(yàn)中航行船與停泊船采用相同船型，皆為MarAd L系列散貨船，船模參數(shù)見表1。航行船模速度為0.518 m/s，對(duì)應(yīng)的Fr為0.074。本次計(jì)算中，由于缺乏試驗(yàn)船型資料，因此，選取了與試驗(yàn)船型相近的KVLCC2作為數(shù)值計(jì)算船型，船型參數(shù)見表2。

表 1 MarAd L系列散貨船模型主尺度參數(shù)Tab. 1 Main scale parameters of MarAd L series bulk carrier model

表 2 KVLCC2實(shí)船與船模的主尺度參數(shù)Tab. 2 Main scale parameters of KVLCC2 real ship and ship model

2.3 計(jì)算前處理

計(jì)算域應(yīng)足夠大，以避免岸壁對(duì)流場的影響，但是考慮到計(jì)算實(shí)際情況，本文將計(jì)算域設(shè)為長方體，船舶前后流域選取為8倍船長，岸壁與水底面方向以水底和無擾自由面為界，計(jì)算域如圖3所示。網(wǎng)格單元為正六面體，在物面邊界同時(shí)設(shè)置棱柱層網(wǎng)格，船艏和船尾處由于曲率變化較大，采用網(wǎng)格局部加密。船體網(wǎng)格如圖4所示。

圖 3 計(jì)算域示意圖Fig. 3 Schematic diagram of the calculation domain

邊界條件的示意圖如5所示。對(duì)于流域入口與出口的邊界均設(shè)定為壓力邊界條件，船體的表面、水底面以及岸壁設(shè)置統(tǒng)一定為無滑移壁面邊界條件。在此數(shù)值模擬中，由于船舶的航行速度較低，可以忽略自由液面帶來的影響，從而可使用剛性自由面假定，自由表面就設(shè)定為對(duì)稱邊界條件。

2.4 動(dòng)網(wǎng)格技術(shù)

基于計(jì)算流體軟件Star ccm+中的動(dòng)網(wǎng)格技術(shù)動(dòng)態(tài)鋪層法進(jìn)行危險(xiǎn)對(duì)遇時(shí)的兩船網(wǎng)格的實(shí)時(shí)更新。動(dòng)網(wǎng)格區(qū)域中的網(wǎng)格會(huì)隨著船舶的航行而隨之整體的移動(dòng)，在船體的前后計(jì)算域邊界處會(huì)進(jìn)行相應(yīng)的分割和合并，其他對(duì)應(yīng)的區(qū)域就會(huì)作為靜止網(wǎng)格區(qū)域，動(dòng)網(wǎng)格與靜止網(wǎng)格區(qū)域之間采用滑動(dòng)網(wǎng)格交界面技術(shù)來實(shí)現(xiàn)流場信息的交換。

由于模型中的一船航行，一船靜止，因此，只需設(shè)定航行船的動(dòng)網(wǎng)格區(qū)域，具體示意圖如圖6所示。圖6中虛線部分表示靜止區(qū)域與運(yùn)動(dòng)區(qū)域之間的交界面，虛線區(qū)域內(nèi)為航行船的動(dòng)網(wǎng)格區(qū)域，兩端橢圓部分為網(wǎng)格的合并或者分裂。需要定義此處邊界為靜止區(qū)域。

2.5 數(shù)值求解

基于有限體積法對(duì)流體運(yùn)動(dòng)的控制方程和湍流方程展開離散，其中壓力項(xiàng)離散采用標(biāo)準(zhǔn)格式，動(dòng)量、湍流動(dòng)能以及湍流耗散采用2階迎風(fēng)格式，速度、壓力的耦合迭代選擇Simplec算法，擴(kuò)散項(xiàng)選擇中心差分格式，湍流模型選擇Realizablek-ε模型。

圖 6 動(dòng)網(wǎng)格示意圖Fig. 6 Schematic diagram of moving grid

3 數(shù)值計(jì)算及分析

為了便于比較分析，將計(jì)算結(jié)果進(jìn)行無量綱化，以便與試驗(yàn)數(shù)據(jù)、理論曲線等進(jìn)行對(duì)比，其中靜止船的橫向力系數(shù)CF和偏航力矩系數(shù)CM定義如下：

式中：F為靜止船所受的橫向力；M表示靜止船偏航力矩。U為航行船速度，Lpp為靜止船舶垂線間長，B為船舶型寬，T為船舶吃水。

3.1 數(shù)值計(jì)算結(jié)果分析

通過求解RANS方程，對(duì)航行船經(jīng)過靜止船時(shí)的三維非定常粘性流場進(jìn)行數(shù)值模擬，主要計(jì)算靜止船所受力和力矩，再將數(shù)據(jù)進(jìn)行無量綱化，并轉(zhuǎn)化為與相對(duì)位置P的函數(shù)關(guān)系，其中：定義P=SP/LC，SP（見圖1）為兩船中心縱向距離，Lc為特征船長，計(jì)算中靜止船舶與航行船舶為相同船型條件下，因此，LC=LPP，LPP為船舶垂線間長。最后將計(jì)算結(jié)果與文獻(xiàn)[10]中的試驗(yàn)結(jié)果進(jìn)行對(duì)比，結(jié)果如圖7和圖8所示。

由圖7和圖8可知，總體上而言，計(jì)算值與試驗(yàn)值的曲線基本吻合，具有較好的一致性。圖中橫向力為正，代表排斥力，為負(fù)則表示吸引力；偏航力矩為正表示使船首向外，為負(fù)則表示使船首向內(nèi)。

對(duì)橫向力和偏航力矩的數(shù)值曲線進(jìn)行綜合分析可知，P=-1～1的范圍是兩船會(huì)遇過程中的危險(xiǎn)區(qū)間。首先，在此范圍前段P=-0.5處，船舶所受橫向力和偏航力矩都出現(xiàn)一個(gè)極值點(diǎn)，此階段最會(huì)出現(xiàn)船舶因穩(wěn)性不夠?qū)е码y以抵消干擾力而出現(xiàn)傾覆現(xiàn)象。之后在P=-0.5～0的過程中偏航力矩和橫向力迅速下降，在0.5附近，兩者均為負(fù)值，并且偏航力矩處于極值點(diǎn)，非常容易引起船舶的相互吸引而導(dǎo)致船首相碰撞。

圖 8 偏航力矩計(jì)算值與試驗(yàn)值的比較Fig. 8 Comparison of calculated values of yaw moments with experimental values

圖 7 橫向力計(jì)算值與試驗(yàn)值的比較Fig. 7 Comparison of calculated values of transverse forces with experimental values

由圖6可知靜止船舶所受橫向力在P為-1，0，1時(shí)附近出現(xiàn)了3次峰值。在p=-2～-21區(qū)間，試驗(yàn)值與計(jì)算值表示船舶都受到正向橫向力并逐漸遞增至峰值。在p=-1～0區(qū)間，船舶所受橫向力從正的峰值遞減至負(fù)的峰值，在p=0.8附近，縱向力接近零，受力方向開始由負(fù)向變?yōu)檎颉５傮w看來，2條曲線的基本趨勢一致。

從圖7偏航力矩曲線變化來看，P=0～2區(qū)間偏航力矩曲線與P=-2～0區(qū)間相似，但方向相反。在P=-2～-1區(qū)間，船舶受到負(fù)的偏航力矩，且力矩值也隨之增大。在p=-1～-0.5，3種曲線均逐步上升，力矩由負(fù)極值增至正極值。在P=0時(shí)，偏航力矩出現(xiàn)零值，此種現(xiàn)象與文獻(xiàn)[11]的研究結(jié)論基本一致，均在此處出現(xiàn)零值。

3.2 流場分析

圖9～圖14給出了不同時(shí)刻自由面上的流場壓力分布圖。所有流場壓力云圖均在同一壓力范圍下進(jìn)行對(duì)比。從圖中可以看出，隨著兩船的接近，船體兩側(cè)的壓力差也不斷變化，干擾現(xiàn)象逐漸明顯。圖9中，航行船向靜止船靠近，但因?yàn)榭v間距較大，對(duì)靜止船的影響較小。圖10中兩船中心縱向相距約1倍特征船長，此刻航行船的船首高壓區(qū)對(duì)靜止船造成影響，靜止船受到排斥力。圖11中兩船相對(duì)位置p=-0.5，靜止船左舷區(qū)域受航行船船首高壓區(qū)和船中低壓區(qū)的共同影響，此刻橫向力很小。圖12中p=0，沿船長方向上靜止船舶壓力平衡，橫向處在航行船舶船中低壓區(qū)，受到一個(gè)極大的指向航行船舶的橫向力。圖13中，p=0.5，靜止船左舷側(cè)受航行船船中低壓區(qū)和船尾高壓區(qū)共同影響，橫向力接近零。圖14航行船的船尾壓力區(qū)對(duì)靜止船產(chǎn)生的影響隨著船舶間的縱向距離的加大逐漸變小。

圖 9 P=-1.5流場壓力分布Fig. 9 P=-1.5 flow field pressure distribution

圖 10 P=-1流場壓力分布Fig. 10 P=-1 flow field pressure distribution

圖 11 P=-0.5流場壓力分布Fig. 11 P=-0.5 flow field pressure distribution

4 結(jié) 語

本文研究采用CFD技術(shù)，基于Star ccm+計(jì)算流體軟件，通過求解RANS方程，選取Realizablek-ε湍流模型，應(yīng)用動(dòng)網(wǎng)格技術(shù)，在剛性自由面的假定下進(jìn)行了三維粘性流場的數(shù)值模擬和船舶危險(xiǎn)對(duì)遇態(tài)勢下的水動(dòng)力計(jì)算，得到如下結(jié)論：

1）本文計(jì)算值與已有試驗(yàn)值基本吻合，表明了本文數(shù)值方法能夠用于兩船危險(xiǎn)對(duì)遇狀態(tài)下水動(dòng)力的數(shù)值模擬。

圖 12 P=0流場壓力分布Fig. 12 P=0 flow field pressure distribution

圖 13 P=0.5流場壓力分布Fig. 13 P=0.5 flow field pressure distribution

圖 14 P=-1流場壓力分布Fig. 14 P=-1 flow field pressure distribution

2）對(duì)兩船對(duì)遇過程中橫向力和偏航力矩變化情況的分析，可以看出P=-1～1為兩船會(huì)遇過程中的難以操縱和發(fā)生碰撞的階段，在航行中一定要謹(jǐn)慎駕駛，加強(qiáng)避讓。

本文研究的成果能夠?yàn)橄乱徊綇?fù)雜情景下的船間水動(dòng)力干擾問題深入討論奠定一定的研究基礎(chǔ)，也為船舶相遇時(shí)的安全操縱提供一定的理論依據(jù)。