船舶運(yùn)動與液艙晃蕩耦合頻域計(jì)算方法

孫雷，羅賢成，劉昌鳳，姜勝超

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

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

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

0 引 言

船舶運(yùn)動與液艙內(nèi)流體的晃蕩會產(chǎn)生相互作用，特別是對于液化天然氣船（LNG）和浮式生產(chǎn)儲油卸油設(shè)備船（FPSO）等含有大型液艙的船舶，其具有較大的液艙，并且液艙的固有頻率接近于船舶運(yùn)動頻率，會導(dǎo)致液艙內(nèi)流體產(chǎn)生的晃蕩作用更加明顯。因此，在載液船舶的設(shè)計(jì)初期，就應(yīng)該考慮船舶運(yùn)動與液艙內(nèi)流體晃蕩的耦合作用影響。近年來，有關(guān)船舶與液艙晃蕩耦合作用的研究已有較大進(jìn)展。Francescutto和 Contento［1］通過對載液船舶在橫浪工況下進(jìn)行規(guī)則波作用下的模型實(shí)驗(yàn)，證明液艙內(nèi)液體晃蕩自由面的存在會對船舶運(yùn)動產(chǎn)生影響。Rognebakke等［2］通過對液艙進(jìn)行橫蕩實(shí)驗(yàn)，獲得了穩(wěn)定狀態(tài)下液艙橫蕩運(yùn)動與線性入射波頻率之間的線性關(guān)系。Molin等［3］通過計(jì)算載有液艙的LNG船，得到了船舶運(yùn)動與艙內(nèi)液體自由水面的變化情況。Newman［4］開發(fā)的WAMIT程序可用于計(jì)算線性液艙晃蕩與船舶運(yùn)動耦合情況。KIM等［5］采用脈沖響應(yīng)函數(shù)（IRF）求解線性船舶運(yùn)動，采用有限差分方法模擬非線性液艙的晃蕩問題，并與載有方形液艙的S175船舶橫搖運(yùn)動實(shí)驗(yàn)結(jié)果進(jìn)行了對比，結(jié)果顯示吻合較好。操戈等［6］考慮了液艙內(nèi)流體的粘性，發(fā)現(xiàn)所得結(jié)果與試驗(yàn)值更接近。

本文將采用頻域計(jì)算方法研究船舶運(yùn)動與液艙晃蕩耦合作用問題。在頻域方法中，對于船舶運(yùn)動方程，采用自由面格林函數(shù)方法進(jìn)行求解；對于液艙內(nèi)流體晃蕩過程，采用Rankine源方法進(jìn)行求解。為提高計(jì)算精度，本文將采用高階邊界元方法進(jìn)行離散，應(yīng)用8節(jié)點(diǎn)等參數(shù)單元進(jìn)行頻域求解，并以文獻(xiàn)［7］中裝載2個液艙的船舶模型為例，將本文計(jì)算結(jié)果與采用水動力軟件HYDROSTAR所得結(jié)果、文獻(xiàn)［8］的數(shù)值計(jì)算結(jié)果及文獻(xiàn)［7］的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，以驗(yàn)證本文計(jì)算方法的準(zhǔn)確性。

1 運(yùn)動預(yù)報(bào)理論

1.1 船舶運(yùn)動

船舶在波浪上的運(yùn)動預(yù)報(bào)是研究船舶運(yùn)動與液艙流體晃蕩的基礎(chǔ)，且波浪幅值與船舶特征尺度相比是小量，流體的粘性影響很小，波浪繞射力對船舶影響較大。因此，在研究船舶在波浪中的運(yùn)動姿態(tài)時(shí)，基于勢流理論，流體通?？梢约僭O(shè)為均勻、不可壓縮、無粘、無旋的理想流體。

波浪中的參考坐標(biāo)系如圖1所示。其中，固定坐標(biāo)系為大地坐標(biāo)系o0-x0y0z0，平面x0o0y0與船舶靜水面重合，o0z0軸垂直向上；船舶位置參考系為o′-x′y′z′，o′x′軸位于中縱剖面，指向船艏為正，o′y′軸指向右舷為正，o′z′軸垂直向上；船舶動坐標(biāo)系為o-xyz，oz軸垂直向上，當(dāng)船有航速時(shí)，該坐標(biāo)系隨船以平均速度移動，代表著船舶運(yùn)動形態(tài)。

圖1 船舶運(yùn)動坐標(biāo)系Fig.1 The coordinate system for ship motions

根據(jù)三維線性勢流理論，采用一階近似，分布在船舶周圍的速度勢函數(shù)滿足線性邊界條件與拉普拉斯方程。當(dāng)船舶靜止在自由水面上搖蕩時(shí)（無航速），流場中的一階不定常速度勢?可以分解為

式中：?0為入射勢；?7為繞射勢；?j為j自由度上的輻射勢；ω為波浪圓頻率；ξj為在j自由度上的位移；i為虛數(shù)單位。

入射勢可取AIRY線性波，繞射勢與輻射勢的規(guī)范化定解條件滿足拉普拉斯方程（式（2））與邊界條件（式（3）和式（4））。

在船體表面SB上，

在自由水面SF上，

式中：g為重力加速度；n為網(wǎng)格單元垂直方向；j為各個自由度。

應(yīng)用格林第二定理，將源點(diǎn)布置在船體表面，并在船體表面建立邊界積分方程：

式中：α為自由項(xiàng)系數(shù)（固角系數(shù)），對于常單元，α=1/2，對于高階邊界元，α隨體表面變化；G為自由面格林函數(shù)，也即Havelock源格林函數(shù)。

本文采用邊界元法將頻域內(nèi)的積分簡化為僅在船體表面積分，大大減小了計(jì)算時(shí)間和數(shù)據(jù)存儲時(shí)間，并采用高階邊界元法離散上述積分方程，使表面離散成等參數(shù)單元，然后代入積分方程進(jìn)行求解，便可得到速度勢大小。得到速度勢值后，在船體瞬時(shí)濕表面上進(jìn)行伯努利積分，便可得到船舶濕表面瞬時(shí)波浪載荷

式中：fex為波浪激振力矩陣；fhydro為輻射力矩陣；frestor為回復(fù)力矩陣；μ為附加質(zhì)量矩陣；λ為輻射阻尼系數(shù)矩陣；C為靜回復(fù)力系數(shù)矩陣；ξ為船舶位移矩陣。

獲得船舶水動力系數(shù)后，可建立船舶運(yùn)動方程（式（7）），得到船舶在各個自由度上的運(yùn)動位移

式中：m為廣義質(zhì)量矩陣；B為船舶系統(tǒng)阻尼矩陣；K為船舶系統(tǒng)剛度矩陣。

1.2 液艙流體晃蕩

研究艙內(nèi)流體的流動時(shí)，選取如圖2所示的坐標(biāo)系。液艙計(jì)算用動坐標(biāo)系o-xyz，原點(diǎn)位于艙內(nèi)自由水面中心；參考坐標(biāo)系o′-x′y′z′不隨物體而搖晃，當(dāng)物體處于平衡位置時(shí)，其與動坐標(biāo)系重合；固定坐標(biāo)系o0-x0y0z0固定于空間某點(diǎn)，不隨模型運(yùn)動而變化?；诰€性勢流理論，艙內(nèi)流體運(yùn)動滿足拉普拉斯方程與以下邊界條件。

在液艙壁面ST上，

在自由水面SF上，

在船體表面SB上，

式中：Ξ =（ξ1，ξ2，ξ3），為液艙平動位移（動坐標(biāo)系相對于參考坐標(biāo)系）；A=（α1，α2，α3），為液艙轉(zhuǎn)動位移（動坐標(biāo)系相對于參考坐標(biāo)系）；rT為位置矢量（相對于參考坐標(biāo)系原點(diǎn)）；nT為方向矢量；η為艙內(nèi)液體波面升高；v=ωk，為考慮自由水面的阻尼系數(shù)［9］，其中k為經(jīng)驗(yàn)系數(shù)。

圖2 液艙運(yùn)動坐標(biāo)系Fig.2 The coordinate system of the tank

根據(jù)線性關(guān)系，可將速度勢函數(shù)進(jìn)行分解：

液艙的邊界條件可以寫為：

在液艙壁面ST上，

在自由水面SF上，

應(yīng)用格林定理，建立邊界積分方程：

式中，G為格林函數(shù)，為簡便計(jì)算，采用簡單Rankine源格林函數(shù)。

同樣，采用高階邊界元法對上述積分方程進(jìn)行離散求解，可以得到速度勢。將速度勢代入伯努利方程，沿著艙內(nèi)瞬時(shí)濕表面進(jìn)行積分，便可得到艙內(nèi)波浪激振力

1.3 船舶運(yùn)動與液艙流體晃蕩耦合作用

通過上述研究，綜合考慮船舶運(yùn)動與液艙晃蕩影響，建立頻域運(yùn)動下船舶與液艙晃蕩耦合運(yùn)動方程為

式中：μext為外部波浪激勵的附加質(zhì)量矩陣；μslosh為液艙流體晃蕩的附加質(zhì)量；Bext為外部波浪激勵的阻尼矩陣；Bslosh為液艙流體晃蕩的阻尼矩陣；Cext為外部波浪激勵的靜回復(fù)力系數(shù)矩陣；Cslosh為液艙流體晃蕩的靜回復(fù)力系數(shù)矩陣。

本文考慮兩者耦合作用計(jì)算的流程圖如圖3所示。

圖3 計(jì)算流程示意圖Fig.3 Flow chart for numerical simulation

2 驗(yàn)證計(jì)算

2.1 波浪與三維結(jié)構(gòu)物作用

本文以某艘FPSO船型為例，計(jì)算在規(guī)則波浪中，船舶在各自由度上的運(yùn)動幅值響應(yīng)算子（RAO），并將本程序計(jì)算結(jié)果與采用HYDROSTAR，AQWA軟件等計(jì)算的結(jié)果進(jìn)行了對比，以驗(yàn)算本程序計(jì)算的準(zhǔn)確性。

2.1.1 計(jì)算模型

該計(jì)算船型的主尺度參數(shù)如表1所示，網(wǎng)格模型如圖4所示。

表1 FPSO船型主尺度參數(shù)Table 1 The main parameters of FPSO ship

圖4 船體網(wǎng)格計(jì)算模型Fig.4 The computational model of ship

2.1.2 計(jì)算結(jié)果與分析

選取橫浪（波浪與船艏方向夾角為90°）與迎浪（波浪與船艏方向夾角為180°）這2種工況，用本程序計(jì)算船舶在規(guī)則波作用下的運(yùn)動響應(yīng)，并分別與HYDROSTAR和AQWA軟件計(jì)算的水動力計(jì)算結(jié)果進(jìn)行對比，結(jié)果如圖5～圖8所示。圖中，PRESENT指本文計(jì)算結(jié)果，HYDROSTAR指采用HYDROSTAR軟件計(jì)算的結(jié)果，AQWA指采用AQWA軟件計(jì)算的結(jié)果，圖中縱坐標(biāo)分別為3個自由度上的運(yùn)動幅值響應(yīng)算子［10］。

圖5 90°浪向下橫搖運(yùn)動響應(yīng)Fig.5 Rolling RAO in 90°wave direction angle

圖6 90°浪向下垂蕩運(yùn)動響應(yīng)Fig.6 Heaving RAO in 90°wave direction angle

圖7 180°浪向下縱搖運(yùn)動響應(yīng)Fig.7 Pitching RAO in 180°wave direction angle

圖8 180°浪向下垂蕩運(yùn)動響應(yīng)Fig.8 Heaving RAO in 180°wave direction angle

圖9～圖14所示為橫浪作用下船舶6個自由度上的波浪激振力（力矩）對比。

圖9 x軸的波浪激振力Fig.9 Wave exciting force in the x axis

圖10 y軸的波浪激振力Fig.10 Wave exciting force in the y axis

圖11 z軸的波浪激振力Fig.11 Wave exciting force in the z axis

圖12 繞x軸的波浪激振力矩Fig.12 Wave exciting moment around the x axis

圖13 繞y軸的波浪激振力矩Fig.13 Wave exciting moment around the y axis

圖14 繞z軸的波浪激振力矩Fig.14 Wave exciting moment around the z axis

由圖可知，本文程序計(jì)算所得船舶在各自由度上的運(yùn)動響應(yīng)和波浪激振力與HYDROSTAR和AQWA水動力軟件計(jì)算的結(jié)果吻合較好，證實(shí)了本文計(jì)算程序的準(zhǔn)確性，也為下文進(jìn)行船舶與液艙晃蕩耦合運(yùn)動提供了較好的基礎(chǔ)。

2.2 考慮液艙晃蕩的船舶運(yùn)動

本文以某艘?guī)?個液艙的FPSO船型為例，分別計(jì)算前、后液艙在不同裝載率下的運(yùn)動響應(yīng)。

2.2.1 計(jì)算模型

選取的計(jì)算船型與液艙的主尺度參數(shù)分別如表2和表3所示。該船型與液艙（左艙為前艙）的平面視圖如圖15所示（圖中數(shù)值單位為m），網(wǎng)格計(jì)算模型如圖16所示。

表2 FPSO船型主尺度參數(shù)（帶2個液艙）Table 2 The main parameters of FPSO ship with two tanks

表3 液艙主尺度參數(shù)Table 3 The main parameters of the tank

圖15 船型與液艙的平面視圖Fig.15 The plane view of the ship and tanks

圖16 船型與液艙耦合網(wǎng)格計(jì)算模型Fig.16 The computational model of the ship and tanks

2.2.2 計(jì)算結(jié)果與分析

本文選取裝載液體為水，計(jì)算無液艙裝載（裝載率為0）以及前、后液艙裝載率分別為（20%，20%），（30%，30%）和（57.5%，43.3%）時(shí)，船舶在浪向?yàn)?90°，120°，180°時(shí)各自由度下的運(yùn)動幅值響應(yīng)算子（基于船舶重心處），并與水動力軟件HYDROSTAR計(jì)算結(jié)果、文獻(xiàn)［8］的數(shù)值模擬結(jié)果以及文獻(xiàn)［7］的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，計(jì)算結(jié)果如下。

圖17～圖22所示為該該耦合船舶在無裝載情況下，浪向分別為 90°，120°和 180°時(shí)的運(yùn)動響應(yīng)幅值。

圖17 90°浪向下橫搖運(yùn)動響應(yīng)（無裝載）Fig.17 Rolling RAO in 90°wave direction angle（0 filling）

圖18 90°浪向下垂蕩運(yùn)動響應(yīng)（無裝載）Fig.18 Heaving RAO in 90°wave direction angle（0 filling）

圖19 120°浪向下橫搖運(yùn)動響應(yīng)（無裝載）Fig.19 Rolling RAO in 120°wave direction angle（0 filling）

圖20 120°浪向下縱搖運(yùn)動響應(yīng)（無裝載）Fig.20 Pitching RAO in 120°wave direction angle（0 filling）

圖21 180°浪向下垂蕩運(yùn)動響應(yīng)（無裝載）Fig.21 Heaving RAO in 180°wave direction angle（0 filling）

圖22 180°浪向下縱搖運(yùn)動響應(yīng)（無裝載）Fig.22 Pitching RAO in 180°wave direction angle（0 filling）

圖23～圖28所示為該耦合船舶在前、后艙裝載率為（20%，20%）情況下，浪向分別為90°，120°和180°時(shí)的運(yùn)動響應(yīng)幅值。

圖29～圖34所示為該耦合船舶在前、后艙裝載率為（30%，30%）情況下，浪向分別為90°，120°和180°時(shí)的運(yùn)動響應(yīng)幅值。

圖35～圖40所示為該耦合船舶在前、后艙裝載率為（57.5%，43.3%）情況下，浪向分別為90°，120°和180°時(shí)的運(yùn)動響應(yīng)幅值。

圖41所示為在90°浪向下，該耦合船舶在液艙不同裝載率下的橫搖運(yùn)動響應(yīng)對比圖。

圖23 90°浪向下橫搖運(yùn)動響應(yīng)（前、后艙裝載率為20%，20%）Fig.23 Rolling RAO in 90°wave direction angle（filling ratio of fore and back tank is 20%，20%）

圖24 90°浪向下垂蕩運(yùn)動響應(yīng)（前、后艙裝載率為20%，20%）Fig.24 Heaving RAO in 90°wave direction angle（filling ratio of fore and back tank is 20%，20%）

圖25 120°浪向下橫搖運(yùn)動響應(yīng)（前、后艙裝載率為20%，20%）Fig.25 Rolling RAO in 120°wave direction angle（filling ratio of fore and back tank is 20%，20%）

圖26 120°浪向下縱搖運(yùn)動響應(yīng)（前、后艙裝載率為20%，20%）Fig.26 Pitching RAO in 120°wave direction angle（filling ratio of fore and back tank is 20%，20%）

圖27 180°浪向下垂蕩運(yùn)動響應(yīng)（前、后艙裝載率為20%，20%）Fig.27 Heaving RAO in 180°wave direction angle（filling ratio of fore and back tank is 20%，20%）

圖28 180°浪向下縱搖運(yùn)動響應(yīng)（前、后艙裝載率為20%，20%）Fig.28 Pitching RAO in 180°wave direction angle（filling ratio of fore and back tank is 20%，20%）

圖29 90°浪向下橫搖運(yùn)動響應(yīng)（前、后艙裝載率為30%，30%）Fig.29 Rolling RAO in 90°wave direction angle（filling ratio of fore and back tank is 30%，30%）

圖30 90°浪向下垂蕩運(yùn)動響應(yīng)（前、后艙裝載率為30%，30%）Fig.30 Heaving RAO in 90°wave direction angle（filling ratio of fore and back tank is 30%，30%）

圖31 120°浪向下橫搖運(yùn)動響應(yīng)（前、后艙裝載率為30%，30%）Fig.31 Rolling RAO in 120°wave direction angle（filling ratio of fore and back tank is 30%，30%）

圖32 120°浪向下縱搖運(yùn)動響應(yīng)（前、后艙裝載率為30%，30%）Fig.32 Pitching RAO in 120°wave direction angle（filling ratio of fore and back tank is 30%，30%）

圖33 180°浪向下垂蕩運(yùn)動響應(yīng)（前、后艙裝載率為30%，30%）Fig.33 Heaving RAO in 180°wave direction angle（filling ratio of fore and back tank is 30%，30%）

圖34 180°浪向下縱搖運(yùn)動響應(yīng)（前、后艙裝載率為30%，30%）Fig.34 Pitching RAO in 180°wave direction angle（filling ratio of fore and back tank is 30%，30%）

圖35 90°浪向下橫搖運(yùn)動響應(yīng)（前、后艙裝載率為57.5%，43.3%）Fig.35 Rolling RAO in 90°wave direction angle（filling ratio of fore and back tank is 57.5%，43.3%）

圖36 90°浪向下垂蕩運(yùn)動響應(yīng)（前、后艙裝載率為57.5%，43.3%）Fig.36 Heaving RAO in 90°wave direction angle（filling ratio of fore and back tank is 57.5%，43.3%）

圖37 120°浪向下橫搖運(yùn)動響應(yīng)（前、后艙裝載率為57.5%，43.3%）Fig.37 Rolling RAO in 120°wave direction angle（filling ratio of fore and back tank is 57.5%，43.3%）

圖38 120°浪向下縱搖運(yùn)動響應(yīng)（前、后艙裝載率為57.5%，43.3%）Fig.38 Pitching RAO in 120°wave direction angle（filling ratio of fore and back tank is 57.5%，43.3%）

圖39 180°浪向下垂蕩運(yùn)動響應(yīng)（前、后艙裝載率為57.5%，43.3%）Fig.39 Heaving RAO in 180°wave direction angle（filling ratio of fore and back tank is 57.5%，43.3%）

圖40 180°浪向下縱搖運(yùn)動響應(yīng)（前、后艙裝載率為57.5%，43.3%）Fig.40 Pitching RAO in 180°wave direction angle（filling ratio of fore and back tank is 57.5%，43.3%）

圖41 90°浪向下液艙不同裝載率下的橫搖運(yùn)動響應(yīng)Fig.41 Rolling RAO with different filling ratio in 90°wave direction angle

由以上圖中對比可以看出：

1）在船舶縱搖和垂蕩運(yùn)動響應(yīng)方面，以上4種計(jì)算結(jié)果均吻合較好，而在計(jì)算船舶橫搖運(yùn)動響應(yīng)時(shí)，幅值結(jié)果和共振頻率結(jié)果誤差較大。這主要是因?yàn)榇翱v、橫向尺寸比較大，因而在橫浪作用下橫搖運(yùn)動較為劇烈。

2）當(dāng)計(jì)算船舶橫搖運(yùn)動響應(yīng)時(shí)，在船舶和液艙共振頻率處，采用HYDROSTAR軟件計(jì)算的結(jié)果其幅值與其他計(jì)算結(jié)果相差較大，這可能是因?yàn)椴捎肏YDROSTAR軟件計(jì)算的沒有考慮自由水面阻尼因素；而本文通過增加自由水面阻尼值，船舶橫搖運(yùn)動響應(yīng)和縱搖運(yùn)動響應(yīng)等計(jì)算結(jié)果均與實(shí)驗(yàn)數(shù)據(jù)以及文獻(xiàn)［7］的數(shù)值計(jì)算結(jié)果吻合較好，這也證實(shí)了本文計(jì)算方法的正確性。

3）由上分析可知，船舶在橫浪（90°浪向）下的橫搖運(yùn)動響應(yīng)較劇烈。因此，通過對比橫浪作用下船舶在不同液艙裝載率情況下的橫搖運(yùn)動響應(yīng)可以發(fā)現(xiàn)，船舶的共振頻率出現(xiàn)了偏移，且當(dāng)裝載率較低時(shí)，船舶的共振頻率偏于低頻，并隨裝載率的上升而向高頻方向移動。當(dāng)液艙裝載液體時(shí)，船舶的橫搖運(yùn)動出現(xiàn)了多個峰值，且當(dāng)前、后液艙裝載率為（20%，20%）時(shí)，第2個峰值最大，說明此時(shí)液艙內(nèi)流體晃蕩得最為劇烈；當(dāng)裝載率高于或者低于此時(shí)的裝載率時(shí)，峰值又降低，說明液艙對船舶運(yùn)動的影響存在一個“最佳”裝載率，在該裝載率下，液艙內(nèi)的流體晃蕩對于船舶運(yùn)動響應(yīng)最大，在船舶設(shè)計(jì)時(shí)需著重予以注意。

3 結(jié) 語

本文通過頻域計(jì)算方法，對船舶在不同液艙液體裝載率情況下的運(yùn)動響應(yīng)進(jìn)行了計(jì)算分析。首先，將本文計(jì)算結(jié)果與其他水動力軟件計(jì)算結(jié)果進(jìn)行對比，驗(yàn)證了本文計(jì)算船舶水動力系數(shù)程序的準(zhǔn)確性。然后，計(jì)算了船舶耦合液艙晃蕩作用時(shí)的運(yùn)動響應(yīng)，將本文計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)及已發(fā)表文獻(xiàn)的數(shù)值結(jié)果進(jìn)行了對比，并采用HYDROSTAR軟件進(jìn)行了計(jì)算。通過對比以上計(jì)算結(jié)果，證實(shí)了本文計(jì)算方法的正確性，并對結(jié)果數(shù)據(jù)進(jìn)行了分析，指出了液艙晃蕩對船舶所造成的影響。本文所做研究對在船舶設(shè)計(jì)時(shí)需考慮液艙晃蕩的影響具有一定的參考意義。