孫雷，陸婷婷，鄧瀟瀟，劉昌鳳

1 大連理工大學(xué) 船舶工程學(xué)院，遼寧 大連 116024

2 高新技術(shù)船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240

3 大連海洋大學(xué) 海洋與土木工程學(xué)院，遼寧 大連 116023

0 引 言

隨著海洋開發(fā)逐漸向深海領(lǐng)域發(fā)展，海洋環(huán)境更加惡劣，在大波幅波浪作用下，海洋結(jié)構(gòu)物載荷與運(yùn)動的非線性特性顯著增強(qiáng)。目前，常用的商用軟件（例如AQWA，Seasam，Hydrostar，Wamit等）都是基于線性或弱非線性的假定進(jìn)行數(shù)值模擬，這必然會導(dǎo)致計算精度不足。因此，針對波浪和海洋結(jié)構(gòu)物相互作用的問題，開展完全非線性數(shù)值模擬研究，對于提高海洋工程裝備的安全性具有重要意義。

1976 年，Longuet-Higgins 等[1]基于勢流理論提出了一種完全非線性方法，采用混合歐拉—拉格朗日方法追蹤自由水面質(zhì)點，自由水面和物面邊界條件得到實時滿足，每一時刻都對網(wǎng)格進(jìn)行重構(gòu)。近年來，隨著計算機(jī)運(yùn)算能力的不斷提高，完全非線性方法得到了快速發(fā)展。周斌珍[2]采用完全非線性方法建立了開敞水域模型，并對漂浮圓柱的受迫振動、繞射問題以及Ringing現(xiàn)象進(jìn)行了深入研究，揭示了非線性作用的機(jī)理。Bai 等[3-6]基于高階邊界元法（HOBEM），運(yùn)用區(qū)域分解技術(shù)，建立了波浪水池的完全非線性三維數(shù)值模型，并對圓柱的受迫振動以及非線性波浪作用下圓柱的運(yùn)動特性及波浪爬高特性等進(jìn)行了分析研究。寧德志[7]采用快速多極子去奇異邊界元法（Fast Multi-pole Boundary Element Method）分別建立了非線性三維開敞水域模型及非線性數(shù)值波浪水槽模型，并對圓球的強(qiáng)迫縱蕩、垂蕩和轉(zhuǎn)動問題以及波浪對圓柱的繞射等問題進(jìn)行了探討。Feng 等[8-9]利用完全非線性數(shù)值模型對并列雙箱在波浪作用下的水動力共振問題進(jìn)行了研究。然而，上述研究都僅局限于簡單幾何形狀物體，鮮有考慮船舶這種帶有復(fù)雜曲面形狀的海洋結(jié)構(gòu)物。

本文將針對某浮式生產(chǎn)儲卸油輪（Floating Production Storage and Offloading，F(xiàn)PSO）模型進(jìn)行數(shù)值模擬研究。首先，綜合考慮自由水面、物面以及入射波浪的非線性影響，基于勢流理論，結(jié)合高階邊界元法，建立波浪與海洋結(jié)構(gòu)物相互作用的三維完全非線性開敞水域模型。然后，采用混合歐拉—拉格朗日（MEL）方法追蹤瞬時自由水面，利用四階Runge-Kutta 法預(yù)報更新下一時間步的水面、物體位移和速度，每一時刻都對自由水面和物體濕表面網(wǎng)格進(jìn)行更新重構(gòu)。最后，通過對自由漂浮的圓柱、無航速的Wigley 船模型進(jìn)行模擬，驗證所提模型的準(zhǔn)確性，并模擬某FPSO模型的水動力特性。

1 數(shù)學(xué)模型

1.1 初邊值問題

圖1 所示為建立的完全非線性開敞水域模型。圖中：SF為瞬時自由水面邊界；SB為瞬時物體濕表面；SD為水底邊界；d為水深；n為物面單元單位法向矢量，以指出流體方向為正方向。定義2 組坐標(biāo)系，分別描述流體域和物體的運(yùn)動。一個是固定坐標(biāo)系O-xyz，其中z= 0 位于靜水面上，垂直向上為正；另一個是船體坐標(biāo)系O′-x′y′z′，其中原點O′在船體質(zhì)心上。當(dāng)物體靜止處于平衡位置時，2 組坐標(biāo)系相互重合。

圖 1 坐標(biāo)系和計算域定義Fig. 1 Definition of coordinate system and computational domain

假定流體為無黏性、不可壓縮、無旋運(yùn)動的理想流體，則整個流域可用速度勢函數(shù)描述，總的流場速度勢?滿足Laplace 方程：

由物質(zhì)導(dǎo)數(shù)和伯努利方程推導(dǎo)可得瞬時自由水面的完全非線性運(yùn)動學(xué)和動力學(xué)邊界條件。本文采用混合歐拉—拉格朗日方法追蹤瞬時自由水面，其邊界條件在拉格朗日系統(tǒng)下可以寫為

式中： d/dt=?/?t+??·?，為物質(zhì)導(dǎo)數(shù)；X(x,y,z)為瞬時自由水面流體任意質(zhì)點的位置矢量；η為自由水面高程；g為重力加速度。

在瞬時物體濕表面SB上，滿足如下邊界條件：

式中，V為物面單元速度矢量。針對固定邊界，滿足V= 0；針對運(yùn)動的物體，V可以表示為

式中：U為結(jié)構(gòu)物的平動速度；Ω為結(jié)構(gòu)物的轉(zhuǎn)動速度；rb為位置矢量。

1.2 速度勢分離

本文采用Ferrant 等[10-11]提出的入射勢和散射勢 分 離 技 術(shù)，將 總 波 浪 (?,η)分 解 為 入 射 波 浪(?I,ηI)和 散射波浪 (?S,ηS)，即

式中的入射波浪通過高階Stokes 理論解給定。將式(6)代入式(1)～式(4)，整理可得關(guān)于散射波浪(?S,ηS)的初邊值問題。

式中：r為流體質(zhì)點與空間固定坐標(biāo)系原點O的水平距離；x0和y0為初始時刻流體質(zhì)點的水平坐標(biāo)；μ(r)為阻尼系數(shù)函數(shù)，其作用是對流體出口處的散射波浪進(jìn)行消波，以防止波浪反射對結(jié)構(gòu)物造成二次作用，具體表示為

式中：ω為波浪圓頻率；λ為入射波波長；r0為阻尼層起點半徑；α0和β0為關(guān)系數(shù)，本文中，α0=β0=1。

為了避免初始效應(yīng)導(dǎo)致數(shù)值不穩(wěn)定，干擾模擬進(jìn)程，本文在與(?I,ηI)有關(guān)的項上作用一個緩沖函數(shù)：

式中，Tm為緩沖時間，一般取入射波周期的整數(shù)倍，這里取為2 倍。

時域分析理論是一種完全模擬實際作用過程的理論，其控制方程中不僅含有空間項，還含有時間項，因此還需要給出初始條件：

1.3 邊界積分方程的建立

在整個計算域內(nèi)應(yīng)用格林第二定理，將上述關(guān)于散射波浪(?S,ηS)的邊界值問題轉(zhuǎn)化成如下的邊界積分方程進(jìn)行求解：

式中：p= (x0,y0,z0)，為源點坐標(biāo)；q= (x,y,z)，為場點坐標(biāo)；α(p)為固角系數(shù)；S1為整個計算域表面；G(p,q)為格林函數(shù)。

由于整個計算域是關(guān)于O-xz平面對稱，且海底是水平的，所以選取簡單Rankine 源及其關(guān)于對稱面（y= 0）和海底（z= ?d）的鏡像為格林函數(shù)。鑒此，在計算時只需考慮一半計算域，而且海底被排除在外，式（17）則為G(p,q)格林函數(shù)。

其中，

這里，采用HOBEM 法，將計算域離散為8 節(jié)點的四邊形單元，并將其表示為參數(shù)坐標(biāo)下的等參單元。在每個單元內(nèi)，則采用二次形函數(shù)來描述任意點的幾何坐標(biāo)、速度勢和速度勢的法向?qū)?shù)，即

式中： (ξ,?)為 參數(shù)坐標(biāo)；Xk，?k， (??/?n)

k，hk分別為單元上第k個節(jié)點的坐標(biāo)、速度勢、速度勢法向?qū)?shù)和形函數(shù)；K為單元內(nèi)節(jié)點的個數(shù)。將式(19)代入式(16)中，積分方程可以離散成如下形式：

式中：J(ξ,?)為雅各比行列式；Ne1和Ne2分別為自由水面和物面邊界上離散的曲面單元數(shù)。上述積分采用標(biāo)準(zhǔn)的Gauss-Legendre 積分，將8 節(jié)點四邊形單元離散為4×4 個樣本點進(jìn)行求解。最終積分方程轉(zhuǎn)換為如下線性方程組：

式中：A11，A12，A21，A22均為系數(shù)矩陣A的元素；B1和B2均為向量B的元素。上述積分邊界每一時刻都在變化，因此每一時間步的系數(shù)矩陣A和向量B都需重新建立再求解。計算過程中，假定當(dāng)前時刻物面上的速度勢導(dǎo)數(shù)和自由水面上的速度勢均已知，由積分方程(20)計算當(dāng)前時刻的物面速度勢和自由水面上的速度法向?qū)?shù)后，采用四階Runge-Kutta 法，基于自由水面邊界條件式(8)～式(11)，計算下一時刻水質(zhì)點位置的和速度勢。自由水面流體質(zhì)點的速度分量可由式(22)計算：

式中：nx，ny和nz為法向向量n沿x，y和z方向的分量；? ?S/?ξ 和? ?S/??可直接通過式(19)求解。需要注意的是，一個節(jié)點可能由幾個單元共用，因此需對節(jié)點周圍單元速度取平均值，得到該節(jié)點的速度分量。

2 數(shù)值實現(xiàn)

因為完全非線性方法的自由水面和物面邊界條件都均實時滿足，所以每個時間步都要重新捕捉及更新自由水面和物體濕表面網(wǎng)格。自由水面和物面的網(wǎng)格重構(gòu)是完全非線性模型工作的重點及難點。

2.1 自由水面網(wǎng)格更新

本文采用彈簧近似法[12-13]對自由水面動網(wǎng)格進(jìn)行處理。彈簧近似法的本質(zhì)是把網(wǎng)格節(jié)點看作由彈簧連接起來的，整個網(wǎng)格像彈簧系統(tǒng)一樣變形。因此，只需要在初始時刻形成一次網(wǎng)格即可。網(wǎng)格節(jié)點位移可由式(23)確定：

式中：δi和δj分別為節(jié)點i和j的水平位移矢量；Ni為網(wǎng)格節(jié)點總數(shù)；kij為節(jié)點i和j之間的彈簧剛度。kij的計算式為

式中，lij為節(jié)點i與j之間的距離。更新后，網(wǎng)格節(jié)點i的水平位置矢量可表示為

至此，新的均勻網(wǎng)格節(jié)點位置已經(jīng)完成更新。下面將通過二次形函數(shù)插值求解該網(wǎng)格節(jié)點上的速度勢以及波面高程。為方便描述，將通過上述方法求得的網(wǎng)格記為網(wǎng)格1，將通過混合歐拉—拉格朗日方法（MEL）方法求得的網(wǎng)格記為網(wǎng)格2。為進(jìn)行插值求解，需首先確定網(wǎng)格1 中的節(jié)點在網(wǎng)格2 中所屬的單元，可以通過式(26)來確定：

式中：Sm為網(wǎng)格2 中第m個單元的面積；Sj為三角形子元素的面積；該三角形由網(wǎng)格1 中所求節(jié)點和第m個單元的節(jié)點構(gòu)成；M1為三角形子元素的個數(shù)。當(dāng)式(26)成立時，網(wǎng)格1 所求節(jié)點一定在第m個單元的內(nèi)部，此時可以通過式(27)來求解節(jié)點的參數(shù)坐標(biāo) (ξ,?)：

式中：(xi,yi)為網(wǎng)格1 上所求節(jié)點的坐標(biāo)；(xkm,ykm)為網(wǎng)格2 上第m個單元節(jié)點的坐標(biāo)。

通過上述處理得到的自由水面的計算域網(wǎng)格相對均勻，可避免因網(wǎng)格變形、扭曲而導(dǎo)致的數(shù)值不穩(wěn)定性。

2.2 物面網(wǎng)格更新

對于船舶這種具有復(fù)雜曲面的海洋結(jié)構(gòu)物來說，物面網(wǎng)格的更新較為困難。為了降低網(wǎng)格更新難度及工作量，本文將物面網(wǎng)格分為靠近水線區(qū)域和遠(yuǎn)離水線區(qū)域，如圖2 所示。

因為遠(yuǎn)離水線區(qū)域始終處于瞬時水線以下，所以這部分網(wǎng)格可以不用更新。針對靠近水線區(qū)域的網(wǎng)格，通過背景網(wǎng)格插值技術(shù)予以更新。初始時刻通過讀入船體各站的型線數(shù)據(jù)進(jìn)行網(wǎng)格劃分，包括水線以上的干舷部分和水線以下的部分，其中水線以上干舷部分的網(wǎng)格只用來更新物面網(wǎng)格，不參與計算。

圖 2 物面區(qū)域劃分示意圖Fig. 2 Schematic diagram of body surface domain division

物面網(wǎng)格均在船體坐標(biāo)系下更新。假設(shè)船體x坐標(biāo)不變，根據(jù)更新后自由水面的水線位置對船體z坐標(biāo)進(jìn)行重新劃分，y坐標(biāo)通過背景網(wǎng)格插值來求解，即：將船體網(wǎng)格投影到O-xz平面，通過式(26)確定當(dāng)前所求節(jié)點在背景網(wǎng)格中所在的單元m。當(dāng)式(26)成立時，所求節(jié)點一定在第m個單元的內(nèi)部。此時，首先通過式(27)求得所求節(jié)點的參數(shù)坐標(biāo) (ξ,?)，然后利用二次形函數(shù)插值求得節(jié)點的y坐標(biāo)。當(dāng)所有節(jié)點的坐標(biāo)均完成更新后，需要將船體坐標(biāo)系下的物面節(jié)點坐標(biāo)轉(zhuǎn)換到固定坐標(biāo)系下，2 組坐標(biāo)系下的坐標(biāo)滿足如下轉(zhuǎn)換關(guān)系：

其中，

式中：(x,y,z)為固定坐標(biāo)系下的坐標(biāo)；(x′,y′,z′)為船體坐標(biāo)系下的坐標(biāo)；(xc0,yc0,zc0)為物體的質(zhì)心位置；( ζ1, ζ2, ζ3)為物體沿坐標(biāo)軸的平動位移；(θ,β,γ)為物體繞坐標(biāo)軸的轉(zhuǎn)動角度。

3 數(shù)值計算及討論

3.1 模型驗證

3.1.1 漂浮圓柱模型

為驗證本文模型的準(zhǔn)確性，對文獻(xiàn)[10]中給出的漂浮圓柱模型進(jìn)行模擬。參數(shù)設(shè)置如下：水深d=1 m，圓柱半徑R/d=1/3，吃水深度D/d=1/6，波幅A/d=0.05，阻尼陣為0，波浪周期。入射波浪采用Stokes 5 階波。然后，將本文所建完全非線性模型的模擬結(jié)果與文獻(xiàn)[10]的模擬結(jié)果進(jìn)行對比，結(jié)果如圖3 所示。包括圖3 及以下各圖結(jié)果經(jīng)過無量綱化處理，以便對比。

由圖3 可見，兩者模擬結(jié)果吻合良好，驗證了本文模型的準(zhǔn)確性。另外，由垂向速度時間歷程曲線可以看出，本文結(jié)果中的二階波峰高度并不總是相等，存在能量傳遞的過程。這可能是因為本研究采用了混合歐拉—拉格朗日方法來追蹤瞬時流體質(zhì)點，相較于文獻(xiàn)[10]使用的半混合歐拉—拉格朗日方法，流體質(zhì)點不僅有垂向運(yùn)動，還有水平運(yùn)動，其水平運(yùn)動對垂向運(yùn)動產(chǎn)生了一定影響，進(jìn)而對物體的運(yùn)動特性也有所影響。

分別對垂向位移、速度和加速度的時間歷程曲線進(jìn)行傅里葉分頻處理，可以得到垂向位移、速度以及加速度的幅值譜，如圖4 所示。從圖中可以看出：物體的運(yùn)動除了主頻之外還有很大的二階成分，加速度甚至還有其他更高階的成分，因此物體在非線性波浪中的運(yùn)動具有很強(qiáng)的非線性效應(yīng)。而實際海浪大多為非線性波浪，因此，采用完全非線性模型對海洋結(jié)構(gòu)物的運(yùn)動特性進(jìn)行模擬，有利于提高海洋結(jié)構(gòu)物的安全性。

圖 3 垂向位移、速度和加速度的時歷曲線Fig. 3 Time history curves of vertical displacement, velocity and acceleration

圖 4 垂向位移、速度和加速度幅值譜Fig. 4 Amplitude spectrum of vertical displacement, velocity and acceleration

3.1.2 零航速Wigley 船

為驗證本文模型對不同形狀物體的模擬效果，對文獻(xiàn)[14]中采用的無航速Wigley-Ⅲ船模型進(jìn)行模擬，并將模擬結(jié)果與該文獻(xiàn)的實驗結(jié)果以及文獻(xiàn)[15]的數(shù)值計算結(jié)果進(jìn)行對比。

該模型的船長L=3 m，船寬B=0.3 m，吃水深度D=0.187 5 m，重心距船底0.17 m，縱向慣性半徑Ryy=0.75 m，排水量0.078 m3。船型函數(shù)由式(30)確定：

本文只考慮Wigley 船的升沉和縱搖2 個自由度運(yùn)動，其他方向都受到限制。

圖5 所示為Wigley 船的網(wǎng)格劃分示意圖，基于船體對稱性，將Wigley 船的一半劃分為14×6個單元。

圖 5 Wigley 船網(wǎng)格劃分示意圖Fig. 5 Schematic diagram of mesh generation on Wigley hull

圖6 是Wigley 船的升沉運(yùn)動和縱搖運(yùn)動的振幅響應(yīng)（RAO）結(jié)果，該結(jié)果經(jīng)過無量綱化處理。從圖中可以看出：1）本文模型模擬結(jié)果與文獻(xiàn)[14]實驗結(jié)果以及文獻(xiàn)[15]數(shù)值計算結(jié)果基本吻合；2）在入射波波長較小時，本文的數(shù)值計算結(jié)果沒有出現(xiàn)實驗結(jié)果中的波動現(xiàn)象，這可能由模型實驗的精度誤差或者是流體黏性導(dǎo)致；3）當(dāng) λ/L＞2.0時無實驗結(jié)果，這可能是因為當(dāng)入射波波長較長時，物理模型實驗中的波浪反射現(xiàn)象較為嚴(yán)重，對實驗結(jié)果有較大影響。

圖 6 Wigley 船運(yùn)動響應(yīng)對比Fig. 6 Comparison of motion responses for Wigley hull

模型驗證結(jié)果表明，本文的模型對不同形狀的物體均能得到較好的模擬結(jié)果。

3.2 零航速FPSO 水動力特性研究

利用本文所建完全非線性模型對某FPSO模型的水動力特性進(jìn)行模擬，將本文結(jié)果進(jìn)行無量綱化處理，并與AQWA 計算結(jié)果以及文獻(xiàn)[16]的實驗結(jié)果進(jìn)行對比。其中，本文實驗是在水池中完成的，水池的有效使用尺寸為40 m，最大工作水深為0.7 m，配有實驗室自制的可移動式多向不規(guī)則波造波機(jī)以及微機(jī)控制的數(shù)據(jù)采集系統(tǒng)。

表1 所示為FPSO 的主尺度參數(shù)，為與實驗結(jié)果[16]進(jìn)行對比，選取縮尺比為1∶45，本文對7.5 和9 m 這2 種吃水深度的模型尺度進(jìn)行模擬計算，網(wǎng)格劃分示意圖如圖7 所示。一般情況下，船舶在頂浪或者順浪時的縱搖和升沉運(yùn)動幅度最為劇烈，因此，本文主要研究FPSO 在頂浪和順浪時的升沉及縱搖運(yùn)動響應(yīng)特性。根據(jù)對稱性，在一半的FPSO上劃分了175 個單元，來浪方向為0°（順浪）和180°（頂浪），所模擬的波浪頻率為0.35～1.09 Hz。

表 1 FPSO 主尺度參數(shù)Table 1 Particulars of the FPSO

圖 7 FPSO 網(wǎng)格劃分示意圖Fig. 7 Schematic diagram of the mesh generation on FPSO

圖8 和圖9 分別為7.5 m 吃水深度時FPSO 在順浪和頂浪時的運(yùn)動時歷曲線。從圖中可以看出，在短波情況下（0.90 Hz），其升沉和縱搖運(yùn)動響應(yīng)幅值都較??；在中等波長情況下（0.70 Hz），升沉運(yùn)動幅值并未達(dá)到最大值，但是縱搖固有周期約等于波浪對船體的擾動周期，并發(fā)生諧搖，而縱搖幅值達(dá)到最大；在最大波長情況下（0.35 Hz），其升沉運(yùn)動幅值比中等波長時大，縱搖運(yùn)動幅值比中等波長時小，此時已經(jīng)處于隨波逐流狀態(tài)。

圖 8 波浪0°（順浪）入射工況下7.5 m 吃水深度時FPSO 的運(yùn)動響應(yīng)時歷曲線Fig. 8 Time histories of motion responses of FPSO with 7.5 m draft and 0° incident angle wave

圖 9 波浪180°（頂浪）入射工況下7.5 m 吃水深度時FPSO的運(yùn)動響應(yīng)時歷曲線Fig. 9 Time histories of motion responses of FPSO with 7.5 m draft and 180° incident angle wave

由此可見，針對升沉和縱搖運(yùn)動，波浪對船體擾動的大小與波長和船長的比值有較大關(guān)系：波浪頻率越大，升沉和縱搖運(yùn)動越緩和，故大船在小波浪中的升沉和縱搖運(yùn)動幅度都不會很大；最大縱搖角度出現(xiàn)在中等波長下，此時波長略大于船長，最大波長的縱搖幅值介于小波長和中等波長之間。

圖10 和圖11 分別為9 m 吃水深度時FPSO在順浪和頂浪時的運(yùn)動時歷曲線。從圖中可以看出，其運(yùn)動規(guī)律與7.5 m 吃水深度時的運(yùn)動規(guī)律保持一致，表明吃水深度的改變對FPSO 的運(yùn)動特性影響較小。

圖12 和圖13 分別為順浪和頂浪情況下FPSO的升沉運(yùn)動響應(yīng)結(jié)果。由圖可以看出，本文模擬結(jié)果與AQWA 計算結(jié)果以及文獻(xiàn)[16]的實驗結(jié)果趨勢一致，升沉運(yùn)動響應(yīng)隨著入射波浪周期的增大而增大。在順浪時，不論是長周期波還是短周期波，本文模擬結(jié)果均與AQWA 計算結(jié)果吻合良好；但是在頂浪時，本文模擬結(jié)果與AQWA 計算結(jié)果存在一定差異，尤其是對高頻波浪。從圖中可以看到：在高頻波浪區(qū)段，本文模擬結(jié)果仍出現(xiàn)了局部小峰值，與文獻(xiàn)[16]的實驗結(jié)果吻合得更好，而AQWA 計算結(jié)果未能體現(xiàn)這一現(xiàn)象，表明本文完全非線性模型的模擬結(jié)果更加符合實際。

圖14 和圖15 分別為順浪和頂浪情況下FPSO的縱搖運(yùn)動響應(yīng)結(jié)果，從圖中可以看到，本文模擬結(jié)果與AQWA 計算結(jié)果以及文獻(xiàn)[16]實驗結(jié)果的趨勢一致，F(xiàn)PSO 縱搖運(yùn)動響應(yīng)隨著入射波浪周期的增大先增大，達(dá)到諧搖狀態(tài)，而后開始減小。從圖中還可以看到：本文完全非線性模型在共振區(qū)段的模擬結(jié)果要大于AQWA 計算結(jié)果，與文獻(xiàn)[16]的實驗結(jié)果吻合得更好，且共振區(qū)段的頻帶寬度相較于AQWA 計算結(jié)果變寬了，與文獻(xiàn)[16]的實驗結(jié)果所表現(xiàn)出來的現(xiàn)象一致。這是因為當(dāng)波浪的頻率與FPSO 的固有振蕩頻率接近時，F(xiàn)PSO 會產(chǎn)生共振，導(dǎo)致運(yùn)動的非線性效應(yīng)加強(qiáng)；另外由于考慮了非線性效應(yīng)，當(dāng)入射波浪與FPSO 相互作用時，可能會誘發(fā)產(chǎn)生其他頻率的波浪，從而導(dǎo)致了共振區(qū)的頻帶寬度變大。傳統(tǒng)的線性方法無法考慮這些非線性因素的影響，對海洋結(jié)構(gòu)物運(yùn)動的預(yù)報偏于危險，不利于海洋結(jié)構(gòu)物的安全性。此外，在低頻區(qū)段，本文模擬結(jié)果以及AQWA 計算結(jié)果和文獻(xiàn)[16]的實驗結(jié)果相差較大，這可能是因為在水池的有限范圍內(nèi)進(jìn)行模型實驗時，難以完全忽略低頻波長波浪反射等因素對實驗的影響（模型實驗水池長40 m，最大低頻波長約為7 m。雖然水池設(shè)有斜坡和消浪材料，仍難以完全消除反射波的影響），而數(shù)值模擬方法可以很好地排除這些因素對結(jié)果的影響。

圖 10 波浪0°（順浪）入射工況下9 m 吃水深度時FPSO 的運(yùn)動響應(yīng)時歷曲線Fig. 10 Time histories of motion responses of FPSO with 9 m draft and 0° incident angle wave

圖 11 波浪180°（頂浪）入射工況下9 m 吃水深度時FPSO 的運(yùn)動響應(yīng)時歷曲線Fig. 11 Time histories of motion responses of FPSO with 9 m draft and 180° incident angle wave

圖 12 波浪0°（順浪）入射工況下FPSO 升沉RAOFig. 12 The heave RAO of FPSO with 0° incident angle wave

圖 13 波浪180°（頂浪）入射工況下FPSO 升沉RAOFig. 13 The heave RAO of FPSO with 180° incident angle wave

圖 14 波浪0°（順浪）入射工況下FPSO 縱搖RAOFig. 14 The pitch RAO of FPSO with 0° incident angle wave

圖 15 波浪180°（頂浪）入射工況下FPSO 縱搖RAOFig. 15 The pitch RAO of FPSO with 180° incident angle wave

表2 給出了不同吃水深度和不同波浪入射角時FPSO 的縱搖最大響應(yīng)周期，可以看到本文數(shù)值結(jié)果與AQWA 計算結(jié)果以及文獻(xiàn)[16]的實驗結(jié)果均有所差異，但誤差都在可接受的范圍內(nèi)，表明本文模型能夠得到較為精確的結(jié)果。

圖16 和圖17 分別為順浪和頂浪工況下FPSO的橫搖運(yùn)動響應(yīng)結(jié)果，文獻(xiàn)[16]的實驗結(jié)果與其他2 種計算結(jié)果出現(xiàn)了較大偏差：模型實驗中出現(xiàn)了一定的橫搖角，而其他結(jié)果中基本沒有橫搖響應(yīng)。針對此問題，大多數(shù)學(xué)者[17]認(rèn)為是模型實驗產(chǎn)生的誤差，主要原因有：1）模型難以保證絕對的左右對稱；2）波浪由造波機(jī)產(chǎn)生，在較寬的水池中并排分布多塊造波板，造波板之間可能產(chǎn)生微小的相位偏差；3）其他可能的系統(tǒng)誤差或人為誤差。本文采用開敞水域模型，入射波浪直接由解析解給定，并采用對稱模型計算，避免了產(chǎn)生上述誤差。

表 2 不同入射波時縱搖最大響應(yīng)周期Table 2 Maximum response period of pitch with different incident angle wave

圖18 和圖19 為7.5 m 吃水深度、周期t/T=5時，不同波浪頻率下FPSO 附近的流場圖。從圖中可以看到，遠(yuǎn)場區(qū)域自由水面相對規(guī)則，表明本文使用的人工阻尼層消波方法效果良好。FPSO 對短波產(chǎn)生了較大的擾動，這是因為船長大于波長，F(xiàn)PSO 易對短波的傳播產(chǎn)生影響。由前面的分析可知：最大縱搖角度出現(xiàn)在中等波長下（0.70 Hz），此時波長略大于船長，從圖中可以更直觀地得到該結(jié)論。由圖還可以看出，在頂浪時FPSO 周圍的擾動相較于順浪時要更劇烈一些，會對船體的運(yùn)動產(chǎn)生一定影響，非線性效應(yīng)會更強(qiáng)一些，這與圖12～圖15 所呈現(xiàn)的現(xiàn)象是一致的。

圖20 和圖21 為9 m 吃水深度、周期t/T=5時，不同波浪頻率下FPSO 附近的流場圖。從圖中可以看到，其流場特性與7.5 m 吃水深度時的流場特性相似，再次表明吃水深度對FPSO 的運(yùn)動特性及流場特性影響較小。

圖 16 波浪0°（順浪）入射工況下FPSO 橫搖RAOFig. 16 The roll RAO of FPSO with 0° incident angle wave

圖 17 波浪180°（頂浪）入射工況下FPSO 橫搖RAOFig. 17 The roll RAO of FPSO with 180° incident angle wave

圖 18 波浪0°（順浪）入射工況下7.5 m 吃水深度時FPSO 附近流場圖(t/T=5)Fig. 18 Flow field near the FPSO with 7.5 m draft and 0° incident angle wave (t/T=5)

圖 19 波浪180°（頂浪）入射工況下7.5 m 吃水深度時FPSO 附近流場圖(t/T=5)Fig. 19 Flow field near the FPSO with 7.5 m draft and 180° incident angle wave (t/T=5)

圖 20 波浪0°（順浪）入射工況下9 m 吃水深度時FPSO 附近流場圖(t/T=5)Fig. 20 Flow field near the FPSO with 9 m draft and 0° incident angle wave (t/T=5)

圖 21 波浪180°（頂浪）入射工況下9 m 吃水深度時FPSO 附近流場圖( t/T =5)Fig. 21 Flow field near the FPSO with 9 m draft and 180° incident angle wave (t/T=5)

4 結(jié) 論

本文利用高階邊界元（HOBEM）方法建立了波浪與海洋結(jié)構(gòu)物相互作用的時域三維完全非線性開敞水域模型，采用混合歐拉—拉格朗日方法追蹤瞬時自由水面，利用四階Runge-Kutta 法更新下一時間步的自由水面、物體位移和速度，且每一時刻都對自由水面和物體濕表面網(wǎng)格進(jìn)行更新?；谕耆蔷€性模型的數(shù)值模擬結(jié)果，得到如下結(jié)論。

1) 通過對非線性波浪作用下自由漂浮的圓柱的模擬，驗證了本文模型的準(zhǔn)確性，證明了物體在非線性波浪中的運(yùn)動具有很強(qiáng)的非線性。因此，在對物體運(yùn)動進(jìn)行預(yù)報時，采用線性波代替非線性波不利于物體的安全性。

2) 通過對Wigley 船和某FPSO 模型的模擬，表明了本文模型對不同形狀的復(fù)雜幾何物體均有較好的模擬結(jié)果，表明了本文模型在實際工程中應(yīng)用的可行性。

3) 通過對某FPSO 模型在順浪和頂浪中升沉及縱搖運(yùn)動響應(yīng)的模擬，結(jié)果表明，本文模擬結(jié)果與AQWA 軟件計算結(jié)果以及文獻(xiàn)實驗所得到的運(yùn)動特性基本吻合，這進(jìn)一步驗證了本文模型的準(zhǔn)確性。此外，在考慮非線性效應(yīng)的情況下，本文模型的升沉運(yùn)動受到的影響較小，基本上仍與線性理論所得結(jié)果一致?？v搖運(yùn)動表現(xiàn)出明顯的非線性特性：共振區(qū)段的幅值明顯增大且頻帶寬度變大，這與實驗結(jié)果的現(xiàn)象一致，說明傳統(tǒng)的線性方法尤其是在共振區(qū)段，低估了FPSO 的運(yùn)動響應(yīng)。由于本文采用對稱模型計算，避免了橫搖運(yùn)動模擬產(chǎn)生的誤差。因此，在實際工程應(yīng)用中有必要考慮非線性效應(yīng)，不僅要關(guān)注共振點，還要對其附近區(qū)域予以足夠的重視。