限制水域船舶水動力干擾快速計算方法

黃金， 鄧強， 許辰， 諸葛凌波， 任慧龍， 周學謙

(1.哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱 15001； 2.教育部船舶與海洋工程技術國際聯(lián)合合作實驗室，黑龍江 哈爾濱 150001； 3.中國船舶及海洋工程設計研究院，上海 200001)

船舶在港口、航道等擁擠水域內(nèi)的自主航行，作為船舶智能航行中的難點問題，吸引了國內(nèi)外學者的廣泛關注。大船型船舶由于其尺寸與載重量極大，其受到的水動力在船舶操縱過程中往往是不可忽略的。當船舶在港口、碼頭以及內(nèi)河中航行時，由于空間限制，船舶常常會以較小的距離實現(xiàn)趕超與會遇，而在此操縱過程中船舶受到的水動力干擾極大，并嚴重威脅船舶航行安全。因此，研究限制水船舶水動力干擾的實時預報方法可為船舶航行安全與自動航行提供一定的技術支持。

船-船水動力干擾作為船舶水動力干擾問題的主要內(nèi)容之一，其主要的研究方法主要分為2種：模型試驗與數(shù)值模擬，后者根據(jù)計算原理，可分為粘性流方法與勢流方法。基于模型試驗的方法更能準確可靠地預報實船的水動力性能，因此國外學者很早就對淺水中的船-船水動力干擾進行了模型試驗研究。Dand[1]采用模型試驗的方法模擬淺水中兩船相遇與趕超中的船舶水動力干擾，并研究了航速、船-船橫向距離與水深對干擾力和船體浮態(tài)(升沉、縱傾)的影響。Vantorre等[2]對4個不同的船模展開了全面的模型試驗，并使用回歸分析的方法建立了能夠預測淺水中水動力干擾力極值的數(shù)學模型。

數(shù)值模擬方面，根據(jù)考慮因素的不同，淺水中的船舶水動力干擾的研究方法主要可以分為3類：

1)模擬真實流體流動的粘性流方法；

2)忽略粘性效應的勢流方法；

3)忽略線性效應與自由液面效應的勢流方法。

粘性流方法主要是基于雷諾平均N-S方程(RANSE)求解船體周圍的流場的特征，該方法能夠考慮自由液面的變化與流體的粘性效應，故其計算精度高。劉曉艷[3]使用CFD軟件Fluent，計算了限制水域內(nèi)兩船會遇與追越過程中的船舶間水動力干擾，研究了湍流模型與時間步長對計算結果的影響。高智勇等[4]基于RANS方程研究近距并行兩船的相互干擾效應對船舶操縱性的影響。Pawar等[5]采用三維粘流方法探究了不同航道情況下航行船經(jīng)過系泊船的影響。Fonfach等[6]利用不同的流動模型，研究了粘性與自由液面效應對船舶間水動力干擾的影響，對比結果顯示：粘性對干擾力的影響較小，但當兩船橫向距離非常小時，自由液面的影響是較大的。粘性流方法存在低計算效率的特點，不能用于水動力干擾的實時預報，但相較于其他數(shù)值方法，粘性流模擬最能夠反映真實的流場特性，可用于驗證其他計算方法的可行性。

由于粘性效應在船-船水動力干擾中的影響較小，忽略粘性效應的勢流方法也被廣泛應用于船舶水動力干擾的研究工作中。S?ding等[7]使用Rankine源法分析了發(fā)生在歐洲Elbe River內(nèi)一起真實的船舶碰撞事故。王隸加[8]應用泰勒展開邊界元法和頻域格林函數(shù)，系統(tǒng)研究了淺水中兩船并行航行時相遇和超越情況下干擾力和干擾力矩的問題。Yuan[9]使用計及自由液面變化的三維邊界元法，分別計算了限制水域內(nèi)船舶沿岸壁行駛、船舶進入船閘及兩船相遇過程中的水動力干擾力，對比粘性流結果與勢流結果顯示：即使沒有考慮粘性作用，勢流方法仍能獲得較高精度的數(shù)值結果。但該方法的計算效率較低，求解一個時間步大概需要幾分鐘的時間，遠遠不滿足實時計算的要求。

由于港口、航道等限制水域內(nèi)的波浪較小，又因為限制水域內(nèi)的航速限制，船舶運動引起的興波較小，不計自由液面效應的勢流理論方法也可用于預報限制水域的船舶水動力干擾問題。Sutulo等[10]使用Hess & Smith面元法[11]求解深水中的船舶水動力干擾問題，由于基于低傅汝德數(shù)假設，使用了合模方法[12]處理自由液面條件，該算法的計算效率極快。徐華福[13]基于同樣的假設，使用高階面元法研究了水深、船船橫向距離與船速對淺水中的船舶水動力干擾的影響。相較于常值面元法，高階面元法具有更高的魯棒性，但其計算效率仍不滿足實時計算的要求，采用類似高階面元的研究還有王隸加[8]。Xu 等[14]研究了基于勢流理論的水動力干擾預報方法中不對稱網(wǎng)格的誤差特性。Huang 等[15]基于勢流理論提出了一個船岸干擾水動力預報方法，與實驗結果和基于RANSE的粘性數(shù)值模擬結果相比，除水深吃水比極小的工況外，基于勢流理論的數(shù)值計算精度均在可接受范圍內(nèi)。

由上述文獻可知，使用合模方法模擬自由液面是一種十分高效的方法。并且在文獻[10]可以看到，基于常值面元法與高階面元法獲得的橫向力與轉首力矩差異很小。因此，本文采用常值面元法實時預報限制水域內(nèi)的船舶水動力干擾。使用鏡像法處理水平水底處的邊界條件，基于低傅汝德數(shù)假設，使用合模方法模擬自由液面。分別模擬船舶經(jīng)過無航速船舶與兩船相遇過程中的船舶水動力干擾，通過與試驗結果對比，驗證數(shù)值算法的可行性與適用性，確定所用假設的合理性。

1 船舶水動力干擾實時計算方法

1.1 坐標系定義

如圖1所示，對于限制水域內(nèi)的兩船之間的水動力干擾問題，建立O-ξηζ，o1-x1y1z1和o2-x2y2z23個坐標系。其中，大地坐標系O-ξηζ的ξOη平面與靜水面重合，ζ軸豎直向下。o1-x1y1z1和o2-x2y2z2分別與Model 1和Model 2固連，x軸指向船艏，y軸指向船體右舷，z軸方向與ζ軸相同。

圖1 坐標系定義Fig.1 Definition of coordinate systems

1.2 控制方程與邊界條件

本文基于流體不可壓縮，無粘性，流動無旋的理想流體假設，采用勢流理論實時求解淺水中兩船趕超與相遇過程中的水動力干擾問題。流域內(nèi)合速度勢Φ(ξ,η,ζ,t)可表示為：

Φ=Vcurξξ+Vcurηη+φ

(1)

式中：Vcurξ和Vcurη分別是水中水平均勻流Vcur的縱向與橫向分量。φ(ξ,η,ζ,t)為誘導速度勢，誘導速度勢的梯度為誘導速度VI=φ。在任意時刻，流域中任意一點的誘導速度勢均滿足拉普拉斯方程：

Δφ=0

(2)

誘導速度勢φ在船體濕表面上滿足不可穿透條件：

(3)

式中：n是船體濕表面上單位外法向量;Vr是局部相對速度：Vr=Vi-Vcur;Vi是第i條船的速度。在水平水底上，速度勢φ同樣滿足不可穿透邊界條件：

(4)

基于低傅汝德數(shù)假設，忽略船體興波對船舶水動力干擾的貢獻，使用合模法簡化自由液面處的邊界條件：

(5)

勢流理論中，流域內(nèi)任意一點某物理量的值可由邊界上的值表示，邊界元積分方程為：

(6)

式中:S是2條船的濕表面;σ為源強;M(x,y,z)和P(x′,y′,z′)分別為場點和源點;f(M)=Vr(M)·n(M);G為格林函數(shù)。由于使用鏡像法和合模模型分別處理水平水底與自由液面處的邊界條件，如圖2所示，本文使用鏡像格林函數(shù)：

(7)

圖2 鏡像法示意Fig.2 The diagram of mirror image technique

1.3 數(shù)值求解

用Hess & Smith面元法將船體表面分割成N個四邊形網(wǎng)格，上述積分方程(7)轉化分線性方程組：

(8)

式中：Aij為影響系數(shù)矩陣，未知數(shù)σj為每個面元控制點處的源強。本文使用高斯-賽德爾迭代法求解該線性方程組的近似值。場點M處的誘導速度勢φ與誘導速度VI可分別表示為：

(9)

(10)

船體表面壓力分布由非定常伯努利方程可得：

(11)

式中：ρ是流體密度，VP=VI-Vr。通過沿著船體濕表面對壓力p進行積分，可獲得作用在第k條船的總干擾力：

(12)

1.4 船舶操縱與船間干擾力求解

根據(jù)船舶操縱性方程得到慣性項水動力：

(13)

(14)

對于本文所研究的問題，船體各方向的加速度都為零，式(13)獲得的干擾力即為由船-船水動力干擾力。

2 計算參數(shù)與結果分析

分別模擬淺水中船舶平行趕超另一無航速船與兩船平行會遇2種典型情況下的船-船水動力干擾問題。在模擬過程中，忽略水動力干擾對船舶姿態(tài)、軌跡與速度的影響。通過對比實驗結果，驗證實時計算方法的可行性與計算精度。

2.1 計算模型

以Dand[1]公開發(fā)表的實驗報告中的貨船與油船作為研究對象，兩船的主要參數(shù)見表1，兩船型線見圖3、圖4。

表1 船模參數(shù)表Table 1 Parameters of the models

圖3 貨船型線Fig.3 The body profile of cargo ship

圖4 油船型線Fig.4 The body profiles of tanker

2.2 計算工況

本文從Dand[1]的報告中選取了5個實驗工況，用于驗證船舶水動力干擾實時計算方法。其中2個工況屬于船舶趕超干擾，另外3個屬于會遇問題。這5個工況的主要參數(shù)見表2和表3。

表2 趕超工況參數(shù)表Table 2 Configuration of the overtaking case

表3 會遇工況參數(shù)表Table 3 Configuration of the encounter case

2.3 計算參數(shù)的確定

面元數(shù)量和時間步長是影響計算結果的重要因素，為保證快速計算要求，并保證一定的精度要求，需要對計算參數(shù)進行確定。

面元數(shù)量和時間步長是影響船舶水動力干擾快速計算的重要因素。實時預報要求程序計算時間小于時間步長，而程序計算耗時主要在求解線性方程組，面元數(shù)量是最主要的因素。因此找出滿足面元數(shù)量和計算時間的關系，就能確定不同面元數(shù)量對應的最小時間步長。通過多組數(shù)據(jù)得到如圖5關系曲線。

圖5 面元數(shù)量與計算時間關系Fig.5 Number of panels versus calculation time

由圖可得，本文船型計算時間與面元數(shù)量成二次曲線關系，經(jīng)過最小二乘擬合得到如下關系式：t=0.349 07-4.74×10-4x+3.88×10-7x2。其中t為計算時間，x為面元數(shù)量。

由以上結果可得，為保證一定的精度要求，當面元數(shù)量較小時須采用較小時間步長，反之使用較大步長。因此，為了找到合適的面元數(shù)量和時間步長，對不同面元數(shù)量計算得到的干擾力進行計算并得到相對誤差結果，滿足誤差最小的即為最佳面元數(shù)量，時間步長由上圖確定。

選取面元數(shù)量為2 000，時間步長為0.5 s作為基準，由圖5確定一系列不面元數(shù)量和時間步長，取不同工況計算干擾力結果并計算相對誤差。相關參數(shù)見表4和表5，結果見圖6和圖7。

表4 面元數(shù)與時間步長表Table 4 Number of panels and time step

由圖可得，誤差最小出現(xiàn)在面元數(shù)量1 600附近，根據(jù)5%誤差標準，面元數(shù)取1 520～1 680時，誤差最小。取總面元數(shù)1 586進行計算，其中貨船濕表面離散成798個面元，油船表面離散成788個面元，此時時間步長為0.5 s。兩船的濕表面網(wǎng)格見圖8、圖9。

鏡像次數(shù)n影響水平水底的建模精度，從而影響水動力干擾的預報精度。由圖10可見，當鏡像數(shù)n=7時，水動力干擾的預報結果趨于穩(wěn)定，也驗證了Suluto等[12]的結論。

表5 工況參數(shù)表Table 5 Configuration of the caculation case

圖6 橫向力相對誤差結果Fig.6 Relative error of sway force

圖7 搖艏力矩相對誤差結果Fig.7 Relative error of yaw moment

圖8 貨船濕表面網(wǎng)格劃分(面元數(shù)量818)Fig.8 Discretization of the wetted surface of the cargo ship

圖9 油船濕表面網(wǎng)格劃分(面元數(shù)量768)Fig.9 Discretization of the wetted surface of the tanker

2.4 數(shù)值結果與實驗結果對比

由于Dand的實驗報告中只有作用于油船干擾力的試驗結果，因此本文只對作用于油船的干擾力進行對比。數(shù)值結果與試驗結果都作無因次化處理：

(15)

式中：Cy和Cn分別是橫向力與轉艏力矩系數(shù);B0和T0分別是貨船的型寬與吃水。兩船的相對位置也進行無因次化處理：

(16)

式中：ξp和ξo分別為貨船與油船的縱向位置，L0為油船的垂線間長。當ξ<0時，兩船的相對縱向間距逐漸變小。ξ=0時，兩船船中對齊。ξ>0時，兩船縱向間距逐漸變大。

圖10 不同鏡像次數(shù)下油船受到的水動力干擾力Fig.10 The hydrodynamic interaction forces acting on the tanker obtained with different numbers of images

2.4.1 趕超工況

圖11和圖12是貨船平行趕超無航速油船時，油船受到的船舶水動力干擾力，其中圖11是工況1的模擬結果，圖12是工況2的模擬結果。

圖11 工況1中油船受的水動力干擾力Fig.11 The hydrodynamic interaction forces acting on the tanker in Case 1

圖12 工況2中油船受的水動力干擾力Fig.12 The hydrodynamic interaction forces acting on the tanker in Case 2

由圖11(a)和圖12(a)可見，當兩船距離不斷變小時，作用于油船的橫向力先表現(xiàn)為排斥力，后表現(xiàn)為吸引力，并在ξ=0附近達到最大值。然后橫向力隨著兩船的距離不斷變大也逐漸降低。由圖11(b)和圖12(b)可見，轉艏力矩的峰值大約出現(xiàn)在ξ=±1時。以ξ=-1為例，貨船的船艏與油船的船尾位于同一縱向位置，由于兩船之間的流速較高，貨船的船艏與油船的船艉的壓力會降低。由于此刻低壓區(qū)對應的力臂最大，因此轉首力矩在此時達到極值。當ξ=0時，兩船相對距離最小，船體表面的低壓區(qū)轉移到船體舷側。此時低壓區(qū)的壓力最低，低壓區(qū)面積最大，因此船舶收到的橫向力最大，但由于此時低壓力臂最小，所以作用于船體的搖艏力矩很小。對比兩工況的干擾力，可發(fā)現(xiàn)：工況1的水深大約是工況2的一半，但干擾力的峰值卻增大了4～5倍，說明了水深對船舶水動力干擾的顯著影響，也印證了研究限制水域內(nèi)船舶水動力干擾的重要性。

對比試驗結果可發(fā)現(xiàn)，不計自由液面效應與粘性效應的勢流方法能夠較好地預報船舶趕超過程中的干擾力趨勢；但由于忽略了粘性和自由液面效應，在干擾力極值附近，數(shù)值結果數(shù)值結果均小于試驗結果。對比兩工況的預報極值可發(fā)現(xiàn)：本文使用的數(shù)值方法能夠更好地預測大水深情況下的船舶水動力干擾。當水深較小時，由于忽略了上述因素且未考慮船舶升沉姿態(tài)變化，該勢流方法在干擾力極值附近的計算精度較低。

2.4.2 會遇工況

圖13、圖14和圖15分別是兩船平行相遇過程中的油船受到的船舶水動力干擾，其中圖13是工況3的模擬結果，圖14對應工況4。圖15對應工況5的結果。與2.3.1節(jié)相比，在本節(jié)模擬的會遇工況中，兩船均有一定航速。

圖13 工況3中油船受到的水動力干擾力Fig.13 The hydrodynamic interaction forces acting on the tanker in Case 3

圖14 工況4中油船受到的水動力干擾力Fig.14 The hydrodynamic interaction forces acting on the tanker in Case 4

圖13～15分別是油船在工況3～5中受到的橫向力。與上節(jié)對比發(fā)現(xiàn)，兩船面對面會遇過程中的干擾力變化趨勢與趕超過程相同。對比3個工況發(fā)現(xiàn)：工況4中兩船橫向距離最大，兩船速度最小，所以工況4中的干擾力最小；工況5相較于其他2個工況，兩船的水深傅汝德數(shù)都較大，所以工況5中的轉艏力矩曲線的第1、4極值與2、3相近。

圖15 工況5中油船受到的水動力干擾力Fig.15 The hydrodynamic interaction forces acting on the tanker in Case 5

與試驗結果對比發(fā)現(xiàn)：數(shù)值方法能夠有效地預測兩船會遇過程中船舶水動力干擾的變化趨勢，但數(shù)值結果與試驗結果在干擾力極值區(qū)域仍存在一定差異，且均小于試驗結果。對于數(shù)值結果與試驗結果存在差異的原因，主要為沒有考慮粘性效應和沒有考慮自由液面變化，此外船舶姿態(tài)對船舶水動力干擾也有一定的影響。由于文獻可用于方法驗證的工況較少，本文并沒有進行系統(tǒng)地對比分析，研究不同水深、橫向距離與船舶航速下，以上4個因素對船舶水動力干擾的影響與作用。

3 結論

1)通過與實驗結果對比，驗證了不計自由液面變化與粘性效應的勢流理論方法能夠捕獲限制水域內(nèi)船舶水動力干擾的主要成分，并準確預測其變化趨勢，但該方法在預報干擾力極值時，計算精度仍有不足。從工程應用的角度講，該方法計算效率高，計算精度可接受，從一定程度上可為限制水域船舶避碰與自主航行提供支撐。

2)由于忽略了粘性效應、自由液面效應以及船舶姿態(tài)對船舶水動力干擾的影響，在某些極端工況下該方法精度有限。在中淺水 (h/T>1.5)，中低航速下本方法有著更好的精度。

3)本文僅對該方法進行了初步的驗證，只計算了典型工況，實際應用中需通過系列模型試驗與粘性流數(shù)值模擬等手段對本方法進行更深入和系統(tǒng)的驗證。