張晨曦 ， 鄒早建 ，b， 楊 勇

（上海交通大學a.船舶海洋與建筑工程學院；b.海洋工程國家重點實驗室，上海 200240）

1 引 言

隨著航運業(yè)的快速發(fā)展，港口、航道內(nèi)的船舶密集程度越來越高，航行船舶與其他船舶的會遇、超越變得不可避免，過往船只對停泊的船舶產(chǎn)生的水動力干擾也不容忽略；另一方面，隨著人類開發(fā)海洋步伐的加快，海上補給以及兩船之間原油、天然氣等貨物的過駁輸運等水上作業(yè)日益增加，也使船船之間或船舶與其它海洋結(jié)構(gòu)物之間的水動力干擾問題呈現(xiàn)多樣化。

鄰近兩船在靠近航行時，由水動力相互干擾產(chǎn)生的偏航力和力矩會很大程度上影響船舶的正常航行，嚴重的甚至會使船舶失去控制導致碰撞，造成海難事故。隨著人們對于海上安全的關(guān)注程度不斷提高，和船舶航行安全密切相關(guān)的船舶操縱性也日益得到重視。而船舶操縱性會受到上述各種情況下船舶與其他船舶或海洋結(jié)構(gòu)物的水動力相互作用的不利影響，因此這些情況下的水動力干擾問題成為了當前的研究熱點。不僅如此，隨著現(xiàn)代船舶向大型化發(fā)展，相應發(fā)展較慢的港口、航道相對變窄、變淺，船舶經(jīng)常需要靠近其他船舶航行，船舶之間的水動力干擾作用顯著，帶來了嚴峻的航行安全問題，也使研究船舶之間的水動力相互作用變得更加重要。

研究船舶之間的水動力相互作用通?？梢圆捎媚Ｐ驮囼灧椒ê屠碚撚嬎惴椒?。國外學者很早就對兩船水動力干擾問題進行了試驗研究。采用模型試驗方法，Newton[1]首先研究了深水中兩船相互超越的問題，Müller[2]研究了兩船在狹窄航道中的會遇與超越問題，Remery[3]給出了航行船舶對一系泊船舶的水動力影響，Dand[4]研究了沿兩平行航線會遇與超越的兩船之間的相互作用。近年來，Vantorre[5]通過約束模試驗全面研究了各種情況下的兩船水動力相互作用。

在理論計算方面，兩船非定常水動力干擾問題的研究大多采用勢流方法。Tuck和Newman[6]、Yeung[7]、Kijima[8]、Tim[9]等人利用基于細長體理論的匹配漸進展開法分別對開闊水域淺水和深水中兩船之間水動力相互作用、淺窄航道中兩船之間水動力相互作用以及兩船水動力相互干擾下的下沉與縱傾進行了研究；鄭才土[10]利用奇點疊加法對淺窄航道中的兩船水動力相互干擾進行了研究；張謝東等人[11-12]利用邊界元方法研究了淺水中兩船超越和會遇的情況；陳波和吳建康[13]利用Green-Naghdi方程，考慮波浪影響研究了兩船交錯運行的波浪干涉作用。近年來，隨著CFD技術(shù)的不斷發(fā)展，對于船舶相互作用的研究也逐漸轉(zhuǎn)向粘性流方向。Chen等人[14]利用浸入邊界法，通過求解RANS方程，考慮風和岸壁的影響研究了淺水中兩船和多船水動力相互干擾問題。

本文應用FLUENT軟件，采用動網(wǎng)格及滑動交界面技術(shù)，對兩船在淺水中會遇的三維非定常粘性流場進行了數(shù)值模擬，計算了船-船水動力相互作用。通過將計算得到的水動力與試驗數(shù)據(jù)對比，驗證了本文方法的有效性。進而，通過對比整個會遇過程中船-船水動力相互作用隨船舶橫向間距、水深、船速以及船長的變化規(guī)律，分析了這些因素對兩船之間水動力相互作用的影響。本文研究結(jié)果對于指導淺水航道中的船舶安全操縱具有一定的參考價值。

2 數(shù)學模型

2.1 坐標系

圖1 坐標系Fig.1 Coordinate systems

考慮兩船在淺水中分別以定常航速U1和U2沿平行直線相向而行的會遇情況，如圖1所示。兩船（Ship1和Ship2）的橫向間距為Sp，縱向間距為ST，水深為h。作用在Ship1上的橫向力和首搖力矩分別為F1和M1，作用在Ship2上的橫向力和首搖力矩分別為F2和M2。采用三個右手直角坐標系：空間固定坐標系o-xyz、固定在Ship1上并以U1運動的運動坐標系o1-x1y1z1和固定在Ship2上并以U2運動的運動坐標系o2-x2y2z2，其中xoy平面、x1o1y1平面和x2o2y2平面均位于無擾自由面上。在運動初始時刻，空間固定坐標系的原點位于（ST0/2,Sp/2,0）,ST0為初始縱向間距。

2.2 控制方程

由于兩船在會遇狀態(tài)下的船速相對于正常航行狀態(tài)下的船速較低，因此本文忽略興波的影響，對自由面使用剛性自由面假定，整個流場屬于單相、三維、粘性和非定常的不可壓縮流體流動。對湍流的模擬采用雷諾平均方法，在空間固定坐標系下，控制方程為如下連續(xù)性方程和動量方程：

式中，xi為坐標分量，ρ為流體質(zhì)量密度，ui為平均速度分量，p為流體平均壓力，為雷諾應力項。

綜合考慮計算精度和計算時間，本文選擇RNG k-ε兩方程湍流模型進行計算。RNG k-ε湍流模型的湍流動能k及湍流耗散率ε方程如下：

式中Gk為平均速度梯度引起的湍流動能產(chǎn)生項，表達式為：

其中μt為湍流粘度系數(shù)，其他項的表達式以及常數(shù)為：

3 數(shù)值方法

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

本文研究的是淺水中兩船會遇問題，前后左右流域范圍應設置得足夠大以消除邊界對流場的影響。計算中設定船體表面和水底面為無滑移壁面邊界。由于本文采用的是固定坐標系且考慮的是靜水中的兩船會遇問題，所以，對于流域左右兩側(cè)的邊界，取其距離船體足夠遠，可以認為對流場不產(chǎn)生干擾，因此也設置為無滑移壁面邊界；對于流域前后的邊界，將其取在距離船體前后足夠遠的地方，流動充分發(fā)展，因而在其上設定壓力出口邊界條件。同時，考慮到剛性自由面假定，將自由面（即計算域上表面）設置為對稱邊界。經(jīng)數(shù)值比較驗證，本文計算流域取左右寬度為2.6倍船長，前后長度為8倍船長的范圍；兩船初始位置以區(qū)域中心點等距離放置。邊界條件的設定如圖2所示。

圖2 邊界條件Fig.2 Boundary condition

船體周圍采用四面體網(wǎng)格與六面體網(wǎng)格的混合網(wǎng)格：在船體首部和中部采用六面體網(wǎng)格；在船體尾部幾何比較復雜的部分采用適應性強的四面體網(wǎng)格，同時也便于控制動網(wǎng)格區(qū)域中網(wǎng)格的數(shù)量；在其余區(qū)域采用六面體網(wǎng)格。

3.2 動網(wǎng)格方法

動網(wǎng)格方法適合用來模擬流場形狀由于邊界運動而隨時間改變的問題，包括預定義的運動和事先未知的運動。對于動網(wǎng)格，F(xiàn)LUENT提供了三種網(wǎng)格更新方法：彈簧光順模型（Spring Smoothing），局部重劃模型（Local Remeshing）以及動態(tài)層模型（Layering）。結(jié)合本文所考慮的問題，船體沿固定直線單自由度運動，在動網(wǎng)格區(qū)域中船體前后使用的是六面體網(wǎng)格，適合采用動態(tài)層模型來進行網(wǎng)格更新。動態(tài)層模型的基本原理為根據(jù)運動網(wǎng)格在相鄰邊界處的高度變化來對邊界處網(wǎng)格層進行分割或合并，如圖3所示。

圖3 動態(tài)層模型示意圖Fig.3 Sketch of layering

第j層網(wǎng)格作擴大變化時，如果其高度的最大值 δmax大于臨界值 （1+ αs）δ0，則根據(jù)設定將其劃分為兩層；第j層網(wǎng)格作縮小變化時，如果其高度的最小值δmin小于臨界值αcδ0，則將其與上一層網(wǎng)格合并。在此δ0為理想網(wǎng)格單元高度，αs，αc分別為分割因子與合并因子。同時，為了保證計算的穩(wěn)定和精度，本文船體周圍的網(wǎng)格保持不變，而在前后兩端采用頂部鋪層的動態(tài)層模型。

由于兩船只是沿固定直線作相對運動，因此只需在流域中分別以兩船航行的直線為中心，將寬度稍寬于船體的區(qū)域設定為動網(wǎng)格區(qū)域，如圖2所示。動網(wǎng)格區(qū)域內(nèi)的網(wǎng)格隨船體整體移動，在船體前后的流域邊界處進行分割與合并；其他區(qū)域為靜止區(qū)域，動網(wǎng)格區(qū)域與靜止區(qū)域之間采用滑動網(wǎng)格交界面技術(shù)進行數(shù)據(jù)交換。

3.3 UDF的使用

本文使用用戶自定義函數(shù)（UDF），通過調(diào)用FLUENT自帶的DEFINE宏實現(xiàn)以下目標：

（1）對兩船的速度進行設定；

（2）考慮到隨著船體的移動，力矩中心也隨之移動，以計算時間和船速為參數(shù)設定力矩中心的位置；

（3）調(diào)用FLUENT自帶函數(shù)對每一時間步中船體表面所受的力與力矩進行計算；

（4）將所計算的結(jié)果寫入文件以便進行分析處理。

3.4 數(shù)值求解方法

本文采用FLUENT軟件，以有限體積法離散流體運動的控制方程和湍流方程，其中擴散項采用中心差分格式，對流項、湍流動能方程以及湍流耗散率方程采用二階迎風差分格式，速度、壓力的耦合迭代采用SIMPLEC算法。

4 方法的驗證與計算結(jié)果影響因素的分析

4.1 方法驗證與計算結(jié)果分析

由于缺乏試驗船型數(shù)據(jù)，本文計算所用船型為和試驗船型相近的船型。計算及試驗所用的兩艘船的船型有關(guān)參數(shù)見表1（其中*為用于水動力相互作用比較驗證的船型）。表1中，L、B和D分別為船長、船寬和吃水；CB為方形系數(shù)，λ為縮尺比。

表1 計算及試驗所用船型主要參數(shù)Tab.1 Main parameters of ships used in calculation and experiment

本文在模型尺度下進行了計算，計算中兩船模航速相等，其傅汝德數(shù)分別為0.086 6和0.085 3；計算工況滿足實船和船模的傅汝德數(shù)相等。

為了便于比較，計算結(jié)果表達成無因次的形式。作用于船上的橫向力、首搖力矩以及縱向距離的無因次表達為：

式中，F(xiàn)，M分別為橫向力與首搖力矩，U為船速；船間縱向距離ST在會遇前定義為負值；下標1和2分別表示Ship1和Ship2。

圖4和圖5給出了計算得到的船舶所受橫向力和首搖力矩與試驗結(jié)果的對比。從圖中可以看出，計算結(jié)果與試驗結(jié)果吻合較好。

圖4 橫向力系數(shù)計算值與實驗值的比較Fig.4 Comparison of sway force coefficient between numerical and experimental results

圖5 首搖力矩系數(shù)計算值與實驗值的比較Fig.5 Comparison of yaw moment coefficient between numerical and experimental results

如圖4所示，橫向力的受干擾范圍大致為ST′=-2～2，并且其變化經(jīng)歷了三次極值：在ST′=-2～-0.8處為正，并不斷增加；在船首剛剛交會后（ST′≈-0.83）達到極值并迅速減小，于ST′≈-0.43處過零后直至兩船船舯縱向位置重合附近（ST′=0）達到負的最大值，其幅值為前一次極值的1.6倍以上；后半段基本上是一個對稱的過程，但從相位上看稍有推遲，從負的極值上升到ST′≈0.6處改變方向，在ST′≈0.9附近達到正的第三個極值，但幅值只有第一個極值的1/4；之后隨著繼續(xù)航行，橫向力趨向于零。從圖5可以看出，首搖力矩的變化比較復雜，受干擾范圍大致也為ST′=-2～2，有四次極值，其中第一次極值的幅值明顯大于其他三次：在ST′=-2～-0.9的范圍內(nèi)，力矩為正，并不斷增加至ST′=-0.9處到達極值后迅速下降，在ST′≈-0.45處轉(zhuǎn)變方向并于ST′≈-0.23處達到負的第二個極值；隨著兩船進一步接近，首搖力矩增加至ST′≈0.33時達到正的第三個極值，之后又繼續(xù)下降，在ST′≈0.9處達到負的第四個極值；最后隨著兩船距離的增加趨于零。與橫向力的變化對比，首搖力矩的第一個極值發(fā)生在橫向力第一個極值附近，第二個極值與橫向力第一個過零點基本同步，之后第三個極值稍早于橫向力第二個過零點，最后一個極值與橫向力第三個極值基本同時出現(xiàn)。

結(jié)合橫向力和首搖力矩的變化曲線來看，ST′=-1～1范圍是兩船會遇過程中的危險區(qū)間：首先，在此范圍前段（ST′=-1～-0.5），船舶所受橫向力和首搖力矩都處在第一個極值點附近，極易出現(xiàn)船舶本身的控制力不足以抵消干擾力的情況而造成失控。之后在ST′=-0.5～0的過程中橫向力迅速下降，力矩也達到了負的第二個極值（ST′≈-0.23），此時橫向吸引力雖然沒有達到最大值，但也已經(jīng)超過極值的一半，因此非常容易引起船舶的相互吸引而導致船首碰撞；在后半段，在吸引力還沒有完全消失的時候，力矩又達到了正的第三個極值（ST′≈0.3），這時對應兩船船體后部相對的狀態(tài)，容易發(fā)生船尾的相互碰撞。

圖6 ST′=-1.673時的無擾自由面上壓力分布Fig.6 Pressure on the undisturbed free surface at ST′=-1.673

圖7 ST′=-0.827時的無擾自由面上壓力分布Fig.7 Pressure on the undisturbed free surface at ST′=-0.827

圖8 ST′=0.019時的無擾自由面上壓力分布Fig.8 Pressure on the undisturbed free surface at ST′=0.019

圖9 ST′=0.865時的無擾自由面上壓力分布Fig.9 Pressure on the undisturbed free surface at ST′=0.865

圖10 ST′=1.711時的無擾自由面上壓力分布Fig.10 Pressure on the undisturbed free surface at ST′=1.711

圖11 ST′=2.557時的無擾自由面上壓力分布Fig.11 Pressure on the undisturbed free surface at ST′=2.557

圖6至圖11給出了不同時刻無擾自由面上的壓力分布圖，其中顏色越深表示該區(qū)域壓力值越小。從圖中可以看出，隨著兩船的接近，船體兩側(cè)的壓力差也不斷變化。按照前進方向，船體在接近過程中右側(cè)壓力先是逐漸變大，此時橫向力為正并經(jīng)歷第一次極值；接下來隨著兩船縱向距離的進一步減小，右側(cè)壓力迅速減小并遠低于左側(cè)壓力，在兩船平行時達到吸引力的極值；之后，右側(cè)壓力逐漸增加并再次超過左側(cè)壓力，出現(xiàn)第三個極值，最后兩側(cè)壓力隨著兩船距離的增加而趨于平衡。

4.2 船間水動力干擾的影響因素分析

在此就兩船橫向間距、水深、兩船相對速度以及船長的變化對淺水中會遇船舶船間水動力干擾的影響進行分析。這里使用兩艘相同的KCS船進行計算，只給出Ship1的受力計算結(jié)果。表2給出了計算工況。

表2 不同參數(shù)的計算工況Tab.2 Conditions of parameter variation

圖12 不同橫向間距下的橫向力系數(shù)Fig.12 Sway force coefficient at different Sp

圖13 不同橫向間距下的首搖力矩系數(shù)Fig.13 Yaw moment coefficient at different Sp

圖12和圖13給出了不同橫向間距下的船舶所受橫向力和首搖力矩系數(shù)曲線。從圖中可以看出，橫向間距的減小對橫向力和首搖力矩的每個極值都起到了放大作用，而且距離越小放大作用越明顯。從相位上看，橫向力的第一個極值和第三個極值的出現(xiàn)時刻隨著橫向間距的減小分別推遲和提前，第二個橫向力的極值出現(xiàn)的時刻基本相同。首搖力矩的相位變化也有相同的規(guī)律，中間的兩個較小的極值的相位基本不變，而前后兩個極值的出現(xiàn)分別隨著橫向間距的減小而推遲和提前。圖14和圖15給出了不同水深情況下的船舶橫向力和首搖力矩系數(shù)曲線。從圖中可以看出，與橫向間距的變化類似，水深的減小對于兩船的水動力干擾同樣起到了放大作用，而且水越淺，放大作用也越明顯。由圖15可知，水深對首搖力矩的第一個極值和最后一個極值的幅值的影響比較明顯，對于中間兩個極值附近范圍的影響相對較小。從極值出現(xiàn)的時刻看，水深的影響與橫向間距的影響稍有不同：就橫向力而言，不同水深對于各個極值出現(xiàn)的時刻基本沒有影響；就首搖力矩而言，水深的變化只對第二、三次極值出現(xiàn)時刻有些影響，對于其他兩次極值的影響較小。

圖16和圖17給出了船速比的變化對兩船水動力干擾的影響。從圖中可以看出，隨著Ship2船速的增加，Ship1船的受力會被明顯放大，而且速度越快，放大作用越明顯。從極值出現(xiàn)的時刻看：隨著船速的增加，Ship1所受橫向力的第二次極值出現(xiàn)的時刻不變，前后兩個極值出現(xiàn)的時刻分別稍有提前和推遲；首搖力矩的前兩個極值和后兩個極值的出現(xiàn)時刻隨著速度的增加分別提前和推遲。

圖14 不同水深下的橫向力系數(shù)Fig.14 Sway force coefficient at different h/D

圖15 不同水深下的首搖力矩系數(shù)Fig.15 Yaw moment coefficient at different h/D

圖16 不同船速比下的橫向力系數(shù)Fig.16 Sway force coefficient at different U2/U1

圖17 不同船速比下的首搖力矩系數(shù)Fig.17 Yaw moment coefficient at different U2/U1

圖18和圖19給出了船長比的變化對船間干擾力的影響。這里的工況是保持Ship1的船長不變，按照比例縮放Ship2的大小。從橫向力的對比圖可見，隨著Ship2船長的增加，Ship1所受橫向力在各時刻的絕對值也變大，各曲線交點在過零點附近；由首搖力矩的對比圖可知，除了在第三個極值附近各曲線幾乎重合以外，Ship2船長的增加對于其他三個極值都起到了不同程度的放大作用，其中第一個極值放大最明顯，其次是第四個極值。船長的變化總體來說對橫向力和首搖力矩極值的出現(xiàn)時刻影響較小。

圖18 不同船長比下的橫向力系數(shù)Fig.18 Sway force coefficient at different L1/L2

圖19 不同船長比下的首搖力矩系數(shù)Fig.19 Yaw moment coefficient at different L1/L2

5 結(jié) 論

本文針對淺水中兩船會遇問題，采用CFD軟件FLUENT，選取RNG k-ε湍流模型，對非定常RANS方程進行求解，應用UDF和動網(wǎng)格技術(shù)，在剛性自由面的假定下進行了粘性流場數(shù)值模擬和船間水動力相互作用計算。主要結(jié)論如下：

（1）通過將計算結(jié)果與試驗結(jié)果進行對比，驗證了本文方法能夠較好地應用于計算兩船會遇狀態(tài)下的水動力相互作用。

（2）對會遇過程中橫向力和首搖力矩變化情況的分析表明，ST′=-1～1為兩船會遇過程中的容易失控和發(fā)生碰撞的階段。

（3）分別對比淺水中兩船會遇情況下的橫向力和首搖力矩計算結(jié)果隨船間橫向間距、水深、船速及船長的變化規(guī)律，表明這四種因素對于船舶會遇情況下的水動力相互作用都有著重要的影響。

本文方法及研究結(jié)果可以為研究復雜的船間水動力干擾問題提供參考，為淺水航道中船舶安全操縱提供理論依據(jù)。

[1]Newton R N.Some notes on interaction effects between ships close aboard in deep water[C]//First Symposium on Naval Maneuverability.Washington DC,1960:1-24.

[2]Müller E.Untersuchungen über die gegenseitige Kursbeeinflussung von Schiffen auf Binnenwasserstrassen[J].Schiff und Hafen,1967,19(6):393-406.

[3]Remery G F M.Mooring forces induced by passing ships[C]//Offshore Technology Conference.Dallas,1974,Paper No.2066,1974.

[4]Dand I W.Some measurements of interaction between ship models passing on parallel courses[R].National Maritime Institute,1981,Report R 108.

[5]Marc Vantorre,Erik Larforce,Ellada Verzhbitskaya.Model test based formulations of ship-ship interaction forces for simulation purposes[C]//IMSF-28th Annual General Meeting.Genova,Italy,2001.

[6]Tuck E O,Newman J N.Hydrodynamic interactions between ships[C]//10th ONR Symposium on Naval Hydrodynamics.Cambridge,Mass.USA,1974:35-70.

[7]Yeung R W.On the interactions of slender ships in shallow water[J].Journal of Fluid Mechanics,1978,85:143-159.

[8]Kijima K.Manoeuvrability of ships in confined water[C]//International Conference of Ship Manoeuvrability.London,1987,Paper No.20,1987.

[9]Gourlay T.Sinkage and trim of two ships passing each other on parallel courses[J].Ocean Engineering,2009,36:1119-1127.

[10]鄭才土.淺水中船舶相互作用的一種計算方法[J].船舶工程,1994(6):8-14.

[11]張謝東,劉祖源,吳秀恒.船舶超越時相互作用力理論計算[J].武漢交通科技大學學報,1997,21(3):236-241.

[12]Zhang Xiedong,Wu Xiuheng.Study of hydrodynamic forces of ships in narrow waterway[J].Journal of Hydrodynamics,Ser.B,2003,4:113-117.

[13]陳 波,吳建康.淺水域中兩船交錯運行時的非定常波浪干涉[J].水動力學研究與進展,A輯,2005,20(4):486-491.

[14]Chen H C,Liu T,Huang E T.Application of chimera RANS method for multiple-ship interactions in a navigation channel[C]//Proceedings of the Twelfth International Offshore and Polar Engineering Conference.Kitakyushu,Japan,2002:26-31.