淺水中浮體在波浪作用下運動的三維時域計算模型

張俊生，滕 斌，叢培文

(大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024)

淺水中浮體在波浪作用下運動的三維時域計算模型

張俊生，滕 斌，叢培文

(大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024)

港口中系泊船在波浪作用下運動問題的本質(zhì)是淺水波浪與浮體的相互作用。與深水情況不同，淺水問題應當考慮水底、水域邊界的影響及淺水波浪自身的特性，單一模型很難實現(xiàn)該模擬過程。為此，建立了Boussinesq方程計算入射波和Laplace方程計算散射波的全時域組合計算模型。有限元法求解的Boussinesq方程能使入射波充分考慮到水底、水域邊界的影響和淺水波浪的特性；散射波被線性化，采用邊界元法求解，并以浮體運動時的物面條件為入射波和散射波求解的匹配條件。該方法為完全的時域方法，計算網(wǎng)格不隨時間變動，計算過程較為方便。通過與實驗及其他數(shù)值方法的結果進行比較，驗證了本模型對非線性波面、浮體的運動都有比較理想的計算結果，顯示了本模型對非線性問題具有較好的計算能力。

淺水中浮體運動；非線性波浪作用；Boussinesq方程；Laplace方程；三維時域數(shù)值模型

Abstract:The calculation of the wave-induced motion of a moored ship in a harbor is essentially a problem of calculating the wave-induced motion of a floating body in shallow water.In this case,the influence of the sea bed and boundaries as well as waves must be considered,as it is different from the problem in deep water.It is difficult to solve the problem by using a single model.Therefore,a combination of an improved Boussinesq equation,to calculate incident waves,and the Laplace equation,to calculate scattered waves,is used to establish a numerical model.The Boussinesq equation solved by the finite element method can embody the characteristics in shallow water while simulating incident waves,and scattered waves are solved by the boundary element method under the linear approximation.The two parts are matched through the boundary condition on the surface of the floating body.The present model is a full time-domain one,but computational grids do not change with the time,which makes the calculation convenient.While calculating nonlinear wave elevations and the motion of the floating body induced by nonlinear waves,the comparison of present results with experimental data and results of other numerical methods shows a good agreement and proves the combined model feasible to deal with nonlinear problems.The model is worthy of being applied to calculate the motion of the moored ship in engineering.

Keywords:motion of the floating body in shallow water; action of nonlinear waves; Boussinesq equation; Laplace equation; 3D time-domain numerical model

計算港口內(nèi)系泊船在波浪作用下的運動是一個傳統(tǒng)而遠未解決的問題。準確計算需要考慮多方面問題：波浪的非線性、港口邊界和地形、浮體的運動響應、系泊系統(tǒng)等，整個計算是一個復雜的耦合系統(tǒng)。其中，核心問題就是淺水中波浪作用下浮體運動響應的計算。最早比較完整地給出該問題計算的是Sawaragi & Kubo[1]。其計算是一種頻域解析的方法，整個計算域被分為外域、船體與港池邊界間的內(nèi)域、船底下方與水底間計算域等三個部分，各部分之間通過速度勢和壓力連續(xù)相匹配，進行求解。計算只能假定水底為平地，港池也只能為規(guī)則的矩形，而計算條件完全采用線性假設。Weiler等[2]提出了一個計算岸壁前駐波對船體作用的頻域計算方法，計算中，駐波被分離為兩個反向的行進波進行考慮。鑒于頻域方法對入射波特性、計算域形狀等的限制，為使求解更能適用于真實、復雜的近岸港口水域，于是，一些非線性時域數(shù)值方法被不斷地提出。

最有代表性的是Bingham[3]提出的一種組合型的計算方法，其通過差分法求解Boussinesq方程(以下簡稱“B方程”)計算入射波，通過源匯法在頻域內(nèi)求解線性的散射波激振力、附加質(zhì)量和輻射阻尼，再通過傅里葉變換轉(zhuǎn)化為時域的散射波作用力、附加質(zhì)量和遲滯函數(shù)，最后把入射波和散射波的激振力相加，進而求解含有卷積積分的時域運動方程。B方程的特點是可以將三維問題簡化為二維計算，使計算量大為減少，從而能夠在大面積水底地形變化的淺水水域中模擬非線性波浪。因此，在系泊船運動的計算模型中應用B方程在此后的研究中被不斷地采用，以使計算模型能在一定程度上考慮水域地形、邊界的影響和波浪的非線性。Van Der Molen[4]也提出了類似的計算模型，但散射勢的求解采用了卷積積分求解時域線性解的方法。不過，歸根到底，該方法還是一種利用脈沖響應函數(shù)的算法，不是完全意義上的時域方法。王大國等[5]、Wang等[6]及王海龍等[7]給出了幾個耦合時域模型，盡管形式不同，但本質(zhì)上皆是內(nèi)外域耦合計算的方法，外域采用B方程求解，內(nèi)域采用Euler方程或Laplace方程(以下簡稱“L方程”)計算。該計算思路雖是非線性耦合，也是完全時域的，但最大的不足是這些模型只能適用于物體不動的情況，即模型不能考慮物體的運動，只能計算波面。故在其進行的實驗中，浮體也固定不動，從而僅驗證了波面計算的準確性，但該模型無法滿足實際工程的需求。

這里也采用“組合”計算的方法，以B方程模擬入射波浪，以L方程計算散射波，通過運動浮體的物面條件連接兩部分的計算并求解運動方程。與Bingham[3]和Van Der Molen[4]的模型相比，在浮體運動響應和散射波的計算上，采用沿時間步遞進求解的方法，從而構建一個嚴格意義上的完全時域計算模型。在數(shù)值方法上，則采用有限元法求解B方程，以滿足曲邊界的計算，避免差分法處理邊界不靈活的弊端；采用邊界元法計算浮體的運動響應和散射波，從而既能滿足浮體形狀多樣性的需求，又能簡化三維計算。本文對散射波也做線性化處理。正如Van Der Molen[4]指出的：入射波非線性的影響是重要的，但船在港口內(nèi)受波浪作用的運動幅度比較小，輻射波和繞射波比入射波小很多，因此采用線性方法計算波浪對船體的作用是有效的。也正因此，計算散射波時，地形的影響可以忽略，視水底為平地。

通過與Wang等[6]實驗結果的對比，驗證了本模型對非線性波面的計算能力；通過與Bai等[8]的浮體運動時域計算模型的線性波作用結果和Teng等[9]的二階Stokes波作用下浮體運動頻域計算模型的結果相比，驗證了本模型的組合計算方法和銜接條件的正確性，考察了本模型對非線性波浪作用下的浮體運動的計算能力。

1 模型的建立

1.1入射波的計算

入射波浪的計算無需考慮浮體的存在。鑒于等參單元有限元法無法計算四階及其以上導數(shù)，故采用Beji等[10]的改進型B方程計算入射波浪，方程形式為：

空間離散可以采用四邊形或三角形單元；入射邊界條件則為Fenton[12]的高精度非線性波解析解；固邊壁上的邊界條件按照Engelman等[13]方法處理，可對曲邊界準確求解；根據(jù)計算的具體情況和需要，阻尼層和松弛層可分別布置在計算域出口和入射邊界前，以吸收計算域內(nèi)傳出的波浪。時間積分方面，前四步，采用Runge-Kutta法求解，此后，采用四階Adams-Bashforth-Moulton預報-校正方法進行計算。

1.2散射波的計算

采用L方程求解散射波，其形式為

其中，ηs為散射波波高。物面取浮體平均濕表面，物面條件為

式中：U為浮體上某一點運動速度矢量，n=(n1,n2,n3)為該點處指向物體內(nèi)的物面單位法向向量，φI為入射勢。浮體被看作剛體，其位移和轉(zhuǎn)角矢量分別為ξ=(ξ1,ξ2,ξ3)、χ=(χ1,χ2,χ3)，在小轉(zhuǎn)角假定下有

其中，x(x,y,z)為所求點坐標，xr(xr,yr,zr)為轉(zhuǎn)動中心坐標。將式(7)代入式(6)，并進一步轉(zhuǎn)化，可得

令(ξ4,ξ5,ξ6)=χ，(n4,n5,n6)=(x-xr)n，式(8)可簡寫為

其中，N=(n1,n2,n3,n4,n5,n6)為指向物體內(nèi)的廣義物面單位法向向量，ξ也擴為6維向量(ξ1,ξ2,ξ3,ξ4,ξ5,ξ6)。φI/n=φIn=uIn是散射波和入射波相銜接的條件。uI為入射波的質(zhì)點速度，可由求得，其水平分量uI、vI近似到二階的沿水深分布表達式為

對應的垂向分量wI的表達式為

wI=-

以邊界元法求解散射波時，水底被看作不透水的平地。采用Rankine源及其關于海底的像作為格林函數(shù)

計算中，網(wǎng)格固定不變。通過式(4)、(5)可求得每一時間步自由水面上的ηs、φs；通過式(9)，可求得物面上的φs/n。通過這些邊界條件，利用邊界元法，可求得自由水面上的φs/n(線性情況下，為z=0面上的φs/z)和物面上的φs，前者結合ηs通過式(4)、(5)可進一步求得下一時間步自由水面上的ηs、φs，后者可用于計算浮體受力，利用運動方程求得下一時間步的浮體運動位移，再通過式(9)計算物面上的φs/n，以此往復，在時間域上持續(xù)求解。對于計算散射波時的時間積分，這里采用四階Runge-Kutta法。

1.3浮體運動的計算

對于波浪與浮體作用的問題，浮體運動以及其產(chǎn)生的輻射波必須予以考慮。不同于Bingham和Van Der Molen在計算中使用含有遲滯函數(shù)和附加質(zhì)量的運動方程，本模型在時域內(nèi)直接計算輻射波及其作用力，浮體的運動方程采用下式求解：

其中，M、B、C分別為浮體質(zhì)量陣、黏性阻尼系數(shù)陣、恢復力矩陣，均為66的方陣，M和C需要根據(jù)物體質(zhì)心和質(zhì)量分布計算得出，具體求解過程參見李玉成和滕斌[14]給出的方法；F(t)為六個分量的廣義波浪激振力，包括入射波和散射波的共同作用；G(t)為系泊系統(tǒng)等其它外部作用力和力矩。該方程包含了6個廣義方向的運動：縱蕩(ξ1沿x軸)、橫蕩(ξ2沿y軸)、垂蕩(ξ3沿z軸)、橫搖(ξ4繞x軸)、縱搖(ξ5繞y軸)、回轉(zhuǎn)(ξ6繞z軸)。系泊系統(tǒng)等外部作用力G(t)，如纜繩、護舷等，在具體工程估算中，可以由計算中求出的纜繩、護舷變形結合產(chǎn)品力學性能曲線獲得；理論計算中，可由纜繩、護舷理論模型求得。本文不對G(t)進行具體分析求解，只集中考慮波浪激振力F(t)對浮體的作用，計算后文算例時，僅在式(13)中加入剛度陣K(式左側加入Kξ(t)項)，代替系泊系統(tǒng)，以保證浮體能圍繞某一平衡位置運動。

通過流體動壓強p在物面的平均濕表面Ωb上積分求解波浪激振力：

同樣，p可分為入射波和散射波產(chǎn)生的動壓強pI、ps。根據(jù)B方程理論，pI沿水深分布的二階表達式為

ps的線性結果為

運動方程的時間積分也采用Runge-Kutta法，與散射波時間積分同步求解。如果把式(13)抽象表述為

則浮體運動位移和速度的時間積分表達式為

2 模型驗證

2.1波面計算的驗證

Wang等[6]開展了一個波浪對箱型船繞射的實驗，實驗中箱型船被固定住，并通過浪高儀測定了一些位置的波高?，F(xiàn)利用該實驗結果檢驗本模型對非線性問題的處理能力。

實驗水池長31 m、寬14 m，實驗水深為0.3 m，實驗所用箱型船長l=2.0 m、寬B=0.6 m、高h=0.45 m、吃水dr=0.24 m。具體船與浪高儀的布置如圖1所示。在靜水面上，以箱形船的幾何中心為原點，建立方向如圖1所示的坐標系，表1給出了24個測點的坐標。實驗中，入射波浪周期為T=3.0 s，波高為H0=0.03 m。

由實驗布置可以看出，水池的上下兩個邊壁與船模型的距離并不遠，對散射波浪場亦會產(chǎn)生影響，故計算散射波時需要對兩側邊壁作些處理。如1.2節(jié)所述，在格林函數(shù)中，加入源點關于一側邊壁的像，從而去掉該側邊壁的面積積分，簡化計算。另一側邊壁，則需進行網(wǎng)格劃分，參與計算，其上的邊界條件為φs/n=0。計算中為保證船不動，需取K足夠大，文中取其對角線元素為1020。

圖1 Wang等實驗布置圖 Fig.1 Layout of Wang et al.’s experiment

表1 測點位置Tab.1 Positions of the test points

Wang等給出了波浪場穩(wěn)定后，兩個周期內(nèi)各測點上波高變化曲線。圖2展示了本文模型計算結果與實驗結果的對比。盡管本模型對散射波做了線性化處理，致使在某些測點上存在次波峰的差異，但從波高變化曲線的整體上觀察，本計算與實驗結果相當一致。可見，在散射波線性化處理的情況下，本模型對波浪非線性還是有較強的模擬能力，這從一方面體現(xiàn)了本模型的有效性，也證明了Van Der Molen的說法：散射波相比于入射波很小，可線性化處理。

圖2 各測點上計算波高與實驗結果比較Fig.2 Comparison of wave elevations at the test points

2.2浮體運動計算的驗證

Bai等[8]和Teng等[9]分別給出了開敞水域中浮體在二階Stokes波作用下運動的時域和頻域計算模型，并對模型的有效性進行了充分驗證。本文計算浮體運動的結果將分別與該兩種模型的結果進行比較，以驗證本文對浮體運動計算的能力。

圖3 被測試浮體的示意Fig.3 Sizes of the floating body tested

如圖3所示，取開敞的計算水域水深為d=0.5 m，漂浮于水面上的箱形船長l=4.2 m、寬B=0.6 m、高h=0.5 m、吃水dr=0.3 m，計算中取箱形船為均質(zhì)，質(zhì)心oc為其形心。以質(zhì)心垂直投影到靜水面所得之點為原點o，建立坐標系，x軸沿船長方向，y軸沿船寬方向，z軸沿垂直方向。計算時，取原點o為轉(zhuǎn)動中心。

測試波浪為4組，周期皆為T=1.8 s，波高由低到高分別取H0=0.012 m、0.036 m、0.072 m和0.108 m，相應的波陡分別為H0/L=0.34%、1%、2%、3%，均沿y方向入射。根據(jù)《The Shore Protection Manual (1984)》，可以大體判斷出四組波浪的非線性和適用理論，如圖4所示。很明顯，Case 1(H0=0.012 m)可以被看作線性波，適用于線性理論求解；Case 2(H0=0.036 m)、Case 3(H0=0.072 m)、Case 4(H0=0.108 m)則皆適用于二階Stokes波理論求解。

為便于計算結果的比較，需要保證浮體在波浪激振力作用下于某一平衡位置做簡諧振動，故在計算中需要給出剛度陣K和黏性阻尼系數(shù)陣B。剛度陣中，取k22=450，k44=24，其余元素取0；黏性阻尼系數(shù)陣中，取b22=470，b44=25，其余元素取0。

圖4 測試波浪的適用理論Fig.4 Appropriate theories for the tested waves

圖5給出了本模型和Bai等的浮體運動時域計算模型對Case 1的計算結果?？v蕩ξ1、縱搖ξ5、回轉(zhuǎn)ξ6始終為0，故只給出橫蕩ξ2、垂蕩ξ3、橫搖ξ4的歷時曲線?？梢悦黠@地看到，除因初始狀態(tài)不同而致使前面若干個周期有差異外，兩者的計算結果在浮體運動穩(wěn)定之后完全一致。這證明了本模型作為一種組合模型，計算入射波的B方程和計算散射波的L方程相互間的銜接條件是合理的，匹配過程是準確的。圖6給出了本模型和Teng等的二階Stokes波作用下浮體運動頻域計算模型對Case 2、3、4的計算結果，這里只計算浮體的周期運動，未考慮浮體因平均漂移力而產(chǎn)生的定常位移。同樣只給出ξ2、ξ3、ξ4的歷時曲線，其余方向的運動始終為0?？梢灾庇^地看到，整體上，本模型與二階Stokes理論的頻域計算結果吻合較好，顯示出：盡管對散射波做了線性化處理，但本模型仍能在一定程度上對浮體在非線性波浪作用下的運動進行模擬和計算。進一步觀察可以發(fā)現(xiàn)，在波陡從1%增加到3%的過程中，文中計算的運動幅值會逐漸小于二階Stokes理論的頻域結果。波浪非線性越強，偏小的程度就越略大一些，這是散射波線性化的必然結果。由圖4可知，Case 4的波浪非線性已經(jīng)比較可觀，但這里計算結果與二階Stokes理論結果仍比較接近；而另一方面，一般認為在港池中波浪的非線性不會太強，所以，本模型就系泊船的運動而言，對非線性波浪的計算能力是可以接受的。

圖5 線性波作用下浮體運動的計算結果比較(T=1.8 s,H0=0.012 m)Fig.5 Comparison of calculated results for the motion of the floating body induced by a linear wave (T=1.8 s,H0=0.012 m)

圖6 二階Stokes波作用下浮體運動的計算結果比較Fig.6 Comparison of calculated results for the motion of the floating body induced by 2nd order Stokes waves

3 結 語

港口中波浪作用下系泊船運動問題的核心是計算淺水中浮體在波浪作用下的運動。而相比于深水問題，海底地形、水域邊界的影響和淺水波浪特性成為應該考慮的必要因素，因此，建立含有B方程的組合模型成為一種有效的方法。相比于已有的模型，文中建立的B方程和L方程的組合模型是一個能有效計算浮體運動的完全時域模型。采用有限元法求解B方程模擬入射波，在體現(xiàn)淺水波浪特性和水底影響之外，與差分法相比又能較好地處理曲邊界的影響。散射波的線性化使計算得到簡化，浮體的物面積分也只需在平均濕表面上進行，在時域求解過程中網(wǎng)格也始終保持不變，從而使整個求解過程相對簡單、方便。

通過與Wang等實驗結果的對比，體現(xiàn)了本模型對非線性波浪的計算能力；在計算線性波作用下浮體的運動時，本文結果與Bai等的時域模型結果完全一致，說明了本模型作為組合模型，各計算部分的匹配和計算過程是準確的。在計算二階Stokes波作用下浮體的運動時，本文結果與Teng等的二階Stokes理論頻域模型的結果也基本一致，表現(xiàn)出本模型對非線性問題有較理想的計算能力。

綜上可見，本模型是一個計算過程較簡單、能在較大程度上計算出淺水中浮體在非線性波浪作用下運動響應的組合模型?？梢灶A期，該模型在不斷完善之后，能夠成為工程中估算系泊船運動的有效工具。

[1] SAWARAGI T,KUBO M.The motion of a ship in a harbor basin[C]//Proc.18th ICCE.1982,3：2743-2762.

[2] WEILER O M,DEKKER J.Mooring container ships exposed to long waves[C]//Proc.13th Int.Harbor Congress.2003.

[3] BINGHAM H B.A hybrid Boussinesq-panel method for predicting the motion of a moored ship[J].Coastal Engineering,2000,40：21-38.

[4] VAN DER MOLEN W.Behavior of moored ships in harbors[D].Delft：The Water Research Center of The Delft University of Technology,2006：55-68.

[5] 王大國，鄒志利，唐春安.波浪對箱形船作用的三維耦合計算模型[J].船舶力學，2007，11(4)：533-544.(WANG Daguo,ZOU Zhili,TANG Chunan.Time stepping solutions of nonlinear wave forces on a three-dimensional box-shaped ship in a harbor[J].Journal of Ship Mechanics,2007,11(4)：533-544.(in Chinese))

[6] WANG D G,ZOU Z L,TANG C A.A three-dimensional coupled numerical model of nonlinear waves in a harbor[J].Science in China Series E,Technological Sciences,2008,51(12)：2195-2206.

[7] 王海龍，鄒志利，周亞龍，等.波浪對圓柱作用的三維耦合計算模型[J].中國海洋平臺，2010，25(5)：38-48.(WANG Hailong,ZOU Zhili,ZHOU Yalong,et al.The three-dimensional coupled numerical model of non-linear wave acting on a circular cylinder[J].China Offshore Platform,2010,25(5)：38-48.(in Chinese))

[8] BAI W,TENG B.Second-order wave diffraction around 3-D bodies by a time-domain method[J].China Ocean Engineering,2001,15(1)：73-85.

[9] TENG B,EATOCK R.A BEM method for the second-order wave action on a floating body in monochromatic waves[C]//Proc.1st Int.Conf.on Hydrodynamics.1994：274-281.

[10] BEJI S,NADAOKA K.A formal derivation and numerical modeling of the improved Boussinesq equations for varying depth[J].Ocean Engineering,1996,23(8)：691-704.

[11] LI Y S,LIU S X,YU Y X,et al.Numerical modeling of Boussinesq equations by finite element method[J].Coastal Engineering,1999,37：97-122.

[12] FENTON J D.The numerical solution of steady water wave problems[J].Computers & Geosciences,1988,14：357-368.

[13] ENGELMAN M S,SANI R L.The implementation of normal and/or tangential boundary conditions in finite element codes for incompressible fluid flow[J].International Journal for Numerical Methods in Fluids,1982,2：225-238.

[14] 李玉成，滕斌.波浪對海上建筑物的作用(第二版)[M].北京：海洋出版社，2002：70-73.(LI Y C,TENG B.Wave action on maritime structures (2nd edition)[M].Beijing:China Ocean Press,2002：70-73.(in Chinese))

A 3D time-domain numerical model for floating body motion
induced by shallow water waves

ZHANG Junsheng,TENG Bin,CONG Peiwen

(State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116024,China)

TV139.2

A

10.16483/j.issn.1005-9865.2016.01.001

1005-9865(2016)01-0001-09

2014-11-23

國家自然科學基金面上項目(51379032)；國家自然科學基金重大項目(51490672)

張俊生(1980-)，男，遼寧大連人，助理研究員，博士生，從事波浪對海上結構物作用的研究。E-mail：yzjz009@126.com

滕 斌。E-mail:bteng@dlut.edu.cn