肖汶斌，董文才

(海軍工程大學(xué)艦船工程系，湖北武漢 430033)

近年來，三維移動脈動源格林函數(shù)已逐步應(yīng)用于數(shù)值求解單船和多體船的耐波性問題［1-3］。該類格林函數(shù)嚴(yán)格滿足自由面邊界條件，其作為描述船體一邊搖蕩運動且一邊向前航行狀態(tài)的擾動速度勢是較為合理的。三維移動脈動源格林函數(shù)的原型為一雙重積分形式，數(shù)值積分時存在奇異性規(guī)避和高頻振蕩等困難，之后出現(xiàn)了單重積分和混合積分的改進(jìn)形式［4］。針對格林函數(shù)的改進(jìn)形式，眾多學(xué)者開展了快速穩(wěn)定的數(shù)值積分方法研究［5-6］，有效推動了移動脈動源格林函數(shù)在船舶耐波性預(yù)報中的應(yīng)用。在實際工程問題中，諸如球鼻艏、折角線等船體外形對離散物面網(wǎng)格的要求十分精細(xì)，此時船型的復(fù)雜化必將導(dǎo)致網(wǎng)格數(shù)目的增加和系統(tǒng)計算量的增大，而移動脈動源格林函數(shù)及其偏導(dǎo)數(shù)的數(shù)值計算仍是最為主要的耗時部分。

在浮體水動力研究方面，Wu［7］通過引入復(fù)合格林函數(shù)的概念來表達(dá)線性水動力計算中格林函數(shù)的對稱性，以對速度勢的求解矩陣作簡化處理。劉日明［8］基于任一函數(shù)可寫為對稱函數(shù)和反對稱函數(shù)的思想，將格林函數(shù)的對稱性拓展至浮體多階速度勢的求解問題。針對具有多個對稱面的海洋浮式結(jié)構(gòu)物，王興剛［9］在系泊索纜系統(tǒng)的耦合動力分析中應(yīng)用了格林函數(shù)的對稱性。由于波浪中有航速浮體輻射與繞射問題的復(fù)雜性，目前還尚未見有關(guān)三維移動脈動源對稱性分析的相關(guān)報道。

對互為場點和源點的兩個空間坐標(biāo)點，一般的格林函數(shù)(如零航速格林函數(shù))具有對稱特性，而移動脈動源由于直接計及航速效應(yīng)對船體搖蕩運動的影響，場點和源點互換后格林函數(shù)的數(shù)值并不滿足對稱關(guān)系，這隱含著開展移動脈動源關(guān)于點對點的對稱性分析時不應(yīng)局限于空間中的兩點;另外，應(yīng)用移動脈動源格林函數(shù)法時，邊界元積分方程中增添了線積分項，此時求解速度勢的系數(shù)矩陣形式更為復(fù)雜。

為進(jìn)一步提高有航速船舶耐波性預(yù)報的計算效率，本文以Bessho型移動脈動源為積分內(nèi)核，詳細(xì)分析了格林函數(shù)點對點、線積分及面積分的對稱特性，并將其拓展至源強密度和速度勢的數(shù)值求解之中。針對不同船型在波浪中的水動力性能，開展了基于三維移動脈動源格林函數(shù)的船舶耐波性計算與數(shù)值驗證工作。

1 基本數(shù)學(xué)問題

設(shè)船舶以航速U在無限水深的規(guī)則波中航行，并建立隨船體平動的直角坐標(biāo)系o-xyz，其中ox軸與靜水面重合，且與航速U同向，oz豎直向上。假設(shè)流場中的流體為不可壓縮的理想流體，且流動為無旋流動。在船體搖蕩運動和定常興波均為小量的雙重假定下，并考慮流場的穩(wěn)態(tài)解，入射波作用下場點p(x，y，z)處的非定?？臻g速度勢Φ可寫為如下形式

式中:ζa為入射波的波幅值;ωe為船體遭遇頻率;φj(j=1～6)表示船體第j個模態(tài)單位振蕩運動的輻射勢，ηj為其相應(yīng)模態(tài)下的復(fù)振幅;φ7和φ0分別為空間繞射勢和入射波勢。

流場中的擾動速度勢φj(j=1～7)不僅要滿足場內(nèi)的拉普拉斯方程，而且還需要滿足線性自由面條件、物面不可穿透條件、水底固壁條件和遠(yuǎn)方輻射條件，上述方程和邊界條件即構(gòu)成了波浪中三維頻域船舶水動力的定解問題。

滿足上述邊界條件的速度勢可用點源形式的格林函數(shù)在物面上的合理分布來表達(dá)，從而將定解問題中的相應(yīng)微分方程轉(zhuǎn)換為物面上的積分方程來進(jìn)行數(shù)值求解。若以三維頻域移動脈動源格林函數(shù)為積分內(nèi)核，物面上的速度勢由分布在船體平均浸濕表面S及水線L上的點源誘導(dǎo)而成［10］，此時的邊界元積分方程為

其中:σj為第j個模態(tài)下的源強密度;線積分項的積分對象為各源點，此時n1表示水線L處物面單位外法線向量在oξ軸方向上的分量;G(p，q)為源點 q(ξ，η，ζ)對場點 p(x，y，z)的三維頻域移動脈動源格林函數(shù)，其Bessho型表達(dá)式［11］為

式中:

上述Bessho型格林函數(shù)表達(dá)式中，1/r和1/r'項及其關(guān)于場點p(x，y，z)的偏導(dǎo)數(shù)易于數(shù)值求解;而調(diào)和函數(shù)G*是本文探討的重點，其偏導(dǎo)數(shù)可寫為

將邊界元積分方程式(2)代入物面不可穿透條件之中，物面上第j(j=1～7)個模態(tài)下的源強密度可通過式(5)求解獲得

式中:

2 移動脈動源格林函數(shù)的對稱性分析

對于一般船舶而言，船體型線關(guān)于其中縱剖面具有對稱的幾何特性。在船舶的波浪載荷、運動響應(yīng)及自由面波形特征計算中，若能充分利用格林函數(shù)及其偏導(dǎo)數(shù)的對稱性，僅在船體的左舷或右舷部分布置點源，則此時的計算效率將明顯優(yōu)于全結(jié)構(gòu)的船體水動力模型。這里將隨船坐標(biāo)系的原點建立在船體重心處，在船舶有航速航行的工況下分析格林函數(shù)的對稱特性，進(jìn)而探討B(tài)essho型移動脈動源在快速求解有航速船舶水動力問題時的改進(jìn)措施。

2．1 點對點的對稱性

分析調(diào)和函數(shù)的偏導(dǎo)數(shù)可知，G*y除了由場點和源點的空間坐標(biāo)決定，還與符號函數(shù)的取值特性密切相關(guān)，由于，則G*的偏導(dǎo)數(shù)滿足如下的對稱關(guān)系式

式(7)～(10)的物理意義為，若單位強度的兩源點關(guān)于o-xz平面對稱，且這兩源點分別對應(yīng)著關(guān)于o-xz平面對稱的兩場點，則兩源點在空間場點處產(chǎn)生的擾動速度勢及其ox、oz方向分量也關(guān)于o-xz平面對稱，而oy方向的擾動速度勢分量則是反對稱的。

當(dāng)U=2m/s和ωe=1．4rad/s時，圖1和圖2分別給出了單位強度的兩點源和q2在自由面上的Re(G*)和 Re(G*y)圖像。由圖可見，在2．5和 z=0所確定的矩形區(qū)域內(nèi)，Re(G*)和Re(G*y)分別關(guān)于ox軸對稱和反對稱。對不同水平面下的調(diào)和函數(shù)及其偏導(dǎo)數(shù)開展了計算，相應(yīng)的數(shù)值結(jié)果有效驗證了式(7)～(10)所表達(dá)的格林函數(shù)點對點的對稱特性。

圖1 單位強度的兩點源在自由面上Re(G*)圖像Fig．1 Re(G*)at z=0 plane by two sourceswith unit strength

圖2 單位強度的兩點源在自由面上Re(G*y)圖像Fig．2 Re(G*y)at z=0 plane by two sources with unit strength

2．2 線積分的對稱性

假設(shè)線元上各控制點處的源強密度σj為一常數(shù)，下面推導(dǎo)證明線積分的對稱關(guān)系式。設(shè)qm和qc+m分別為線元ΔLm和ΔLc+m上的控制點，且ΔLm和 ΔLc+m關(guān)于o-xz平面對稱，則 qm和 qc+m也關(guān)于o-xz平面對稱，由點對點的對稱性可知

令(n1m，n2m，n3m)和(n1c+m，n2c+m，n3c+m)分別表示線元ΔLm和ΔLc+m的單位外法線向量，由于一般船體的浸濕表面具有對稱性，則對稱面元處的法向量也是對稱的，即

針對式(2)中的線積分項，調(diào)和函數(shù)的法向偏導(dǎo)數(shù)在線元上的積分可寫為

將式(11)～(14)及式(15)代入式(16)～(17)，可得

式(18)表明，若兩線元和兩場點均關(guān)于o-xz平面對稱，則線元上相同強度點源的調(diào)和函數(shù)法向偏導(dǎo)數(shù)對相應(yīng)場點的線積分滿足對稱關(guān)系。

下面以一對關(guān)于o-xz平面對稱的線元和場點為例來驗證上述結(jié)論的正確性。若直線型線元ΔLm的端點坐標(biāo)為(1，2，0)和(0，1，0)，則 ΔLc+m由(1，-2，0)和(0，-1，0)唯一確定，且假設(shè) ΔLm和ΔLc+m均位于直壁型船體面元上，而兩場點為和假設(shè)線元上均布置單位強度的點源格林函數(shù)，圖3給出了ΔLm在p1點和ΔLc+m在p2點處誘導(dǎo)的速度勢法向分量隨遭遇頻率的變化曲線，其中U=2m/s，ωe=0．5～2．5rad/s。從實部和虛部的角度可知，兩線元在不同頻率下的積分結(jié)果曲線完全吻合，從而驗證了式(18)的正確性。

圖3 線元上調(diào)和函數(shù)的法向偏導(dǎo)數(shù)積分結(jié)果Fig．3 Line integral result of the partial derivative at normal direction for the harmonic function

2．3 面積分的對稱性

依據(jù)Hess-Smith方法，船體平均浸濕表面S可離散為若干個平面四邊形面元［12］，則格林函數(shù)及其偏導(dǎo)數(shù)在面元上的積分可采用高斯積分方法求解，這里假設(shè)面元上各高斯節(jié)點處的源強密度相等［13］。

設(shè)面元 ΔSm和 ΔSc+m關(guān)于o-xz平面對稱，且qm和 qc+m分別為 ΔSm和 ΔSc+m的形心，則 qm和qc+m也關(guān)于o-xz平面對稱。由點對點的對稱性和對稱面元處法向量的對稱性，并參考線積分對稱性的推導(dǎo)過程，可得

式(19)的物理含義為，若兩場點和分布有相同源強密度的兩面元關(guān)于o-xz平面對稱，則面元上調(diào)和函數(shù)的法向偏導(dǎo)數(shù)對相應(yīng)場點的面積分也滿足對稱關(guān)系。

2．4 對稱性在邊界元積分中的應(yīng)用

假設(shè)船舶平均浸濕表面由N塊平面四邊形面元組成，且船體左側(cè)的面元編號為1，2，…，c(c=N/2)。由式(2)可知，第j模態(tài)下的邊界元積分方程可離散化為如式(20)所示的線性代數(shù)方程組。

上述方程中的系數(shù)矩陣Gkm(1≦k，m≦N)表示位于面元ΔSm或線元ΔLm上的源點對場點pk的積分值，其表達(dá)式為

分析可知，系數(shù)矩陣Gkm由四部分構(gòu)成，即船體左側(cè)、右側(cè)面元分別對其自身的影響系數(shù)和船體左側(cè)與右側(cè)面元的相互影響系數(shù)，可記為GLL、GRR、GLR和GRL。依據(jù)上述推導(dǎo)的點對點、線積分及面積分的對稱性，并結(jié)合三維頻域移動脈動源格林函數(shù)的數(shù)學(xué)表達(dá)式可知，子系數(shù)矩陣滿足如下的對稱關(guān)系式

綜合式(7)、式(18)、式(19)和式(22)可知，格林函數(shù)的對稱性有利于數(shù)值計算的快速實施主要體現(xiàn)在以下幾個方面:1)在船體面元網(wǎng)格劃分時，僅需對船體的左舷或右舷部分進(jìn)行離散，這樣可節(jié)省網(wǎng)格預(yù)處理的時間并減小網(wǎng)格的存儲規(guī)模，即格林函數(shù)的對稱性可指導(dǎo)船體面元的網(wǎng)格生成;2)輻射勢和繞射勢求解的耗時大部分集中于移動脈動源格林函數(shù)及其偏導(dǎo)數(shù)的計算，利用格林函數(shù)的對稱性可有效減少系統(tǒng)的內(nèi)存量和大幅提高計算的效率;3)格林函數(shù)的對稱性涉及船體輻射問題和繞射問題的各個方面，如各模態(tài)下源強密度及輻射勢、繞射勢的相應(yīng)系數(shù)矩陣的生成、附加質(zhì)量及阻尼系數(shù)的求解、波浪主干擾力中繞射力的計算等，隨著計算流程的深入，格林函數(shù)的對稱性可持續(xù)傳遞下去。

3 算例分析

3．1 波浪干擾力

船體受到的波浪干擾力為入射波作用力(即傅汝德-克雷洛夫力)和繞射波作用力之和。由于波浪干擾力不涉及輻射問題，而僅計及入射波和繞射波與船體的相互作用，因此常以波浪干擾力作為耐波性數(shù)值計算的首要驗證對象。針對雙體船的船型特點，并依據(jù)移動脈動源格林函數(shù)的對稱特性，僅對雙體船一側(cè)片體的物面上進(jìn)行網(wǎng)格離散，而另一側(cè)的物面網(wǎng)格則通過幾何對稱關(guān)系獲得。分析可知，雙體船左側(cè)和右側(cè)片體上分布點源的相互關(guān)系與單體船的左舷和右舷類似，此時左右片體上對稱物面處的點對點、線積分和面積分也應(yīng)滿足對稱性。

圖4 Lewis型雙體船的物面網(wǎng)格劃分Fig．4 Surfacemesh of Lewis-type catamaran

以規(guī)則波中頂浪航行的Lewis型雙體船為對象，分別計算了船體受到的垂蕩干擾力和縱搖干擾力矩，圖4為該船型的物面網(wǎng)格劃分情況，而圖5和圖6則給出了數(shù)值計算和試驗測試［14］的作用力幅值結(jié)果。這里，Lewis型雙體船模型的主尺度為1．5m ×0．75m ×0．125m，每個片體的寬度為0．25m，速度傅汝德數(shù)Fr=0．3。對比分析圖中的數(shù)據(jù)可知，垂蕩干擾力與縱搖干擾力矩的數(shù)值計算結(jié)果與試驗值十分接近，進(jìn)而驗證了應(yīng)用移動脈動源格林函數(shù)的對稱性在求解分析有航速船體波浪干擾力的可行性。

圖5 頂浪下Lewis型雙體船的垂蕩干擾力幅值Fig．5 Amplitude of heave force on Lewis-type catamaran in head waves

圖6 頂浪下Lewis型雙體船的縱搖干擾力矩幅值Fig．6 Amplitude of pitch moment on Lewis-type catamaran in head waves

3．2 運動響應(yīng)

斜浪航行工況下，船體左右兩側(cè)對稱面元處的源強密度和速度勢必定不相等，這與迎浪和背浪面的邊界條件不一致有關(guān)，然而格林函數(shù)的對稱性仍然適用于邊界元積分方程中系數(shù)矩陣的計算，因為在相同參數(shù)τ下單位強度的點源對固定位置場點的格林函數(shù)及其偏導(dǎo)數(shù)必定相等，即三維頻域移動脈動源格林函數(shù)的對稱性與船型、航行速度和遭遇頻率聯(lián)系密切。

下面以系列60船(CB=0．7)為例來驗證格林函數(shù)的對稱性在求解船舶波浪運動響應(yīng)問題中的可靠性，其中船模主尺度為3．048m×0．435m× 0．174m，速度傅汝德數(shù) Fr=0．2，浪向角為135°(頂浪為180°)。圖7～圖9分別為系列60船在規(guī)則波中斜浪航行時垂蕩、橫搖及縱搖運動響應(yīng)幅值的計算結(jié)果，并給出了文獻(xiàn)［15］的模型試驗值。由圖可知:在無因次遭遇頻率等于3．2附近，各模態(tài)運動響應(yīng)的計算幅值略大于試驗結(jié)果，這與諧振區(qū)內(nèi)船體的大幅運動有關(guān);另外橫搖運動時流體的黏性效應(yīng)也不應(yīng)忽略;然而，在整個試驗頻段內(nèi)數(shù)值計算與模型試驗結(jié)果還是基本吻合的，從而表明格林函數(shù)的對稱性可適用于求解斜浪航行船舶的運動響應(yīng)。

圖7 斜浪下系列60船(CB=0．6)的垂蕩運動響應(yīng)幅值Fig．7 Amplitude of heavemotion of Ser．60 with CB=0．6 in oblique waves

圖8 斜浪下系列60船(CB=0．6)的橫搖運動響應(yīng)幅值Fig．8 Amplitude of rollmotion of Ser．60 with CB=0．6 in oblique waves

圖9 斜浪下系列60船(CB=0．6)的縱搖運動響應(yīng)幅值Fig．9 Amplitude of pitchmotion of Ser．60 with CB=0．6 in oblique waves

3．3 自由面波形特征

有航速航行船舶在自由面產(chǎn)生的波形包括非定常的輻射波和繞射波，采用移動脈動源格林函數(shù)法求解自由面波形特征時，依據(jù)輻射問題和繞射問題先分別求解船體浸濕表面處的源強分布密度，并將相應(yīng)模態(tài)下的速度勢代入自由面動力學(xué)邊界條件，進(jìn)而數(shù)值模擬船體周圍的非定常波形。在頂浪航行工況下，船體左右對稱物面上的源強分布密度和速度勢均滿足對稱關(guān)系，故右舷物面的源強和速度勢可由左舷物面通過對稱關(guān)系映射獲得，進(jìn)而模擬自由面處的波形特征，此時系統(tǒng)內(nèi)存量和計算耗時量都大大減少。

圖10 Wigley船模在自由面產(chǎn)生的繞射波形Fig．10 Diffraction waves produced by Wigley ship model on free surface

以一修正的Wigley船模所產(chǎn)生的繞射波為對象，開展了格林函數(shù)對稱性在模擬自由面波形特征中的適用性驗證。該船模的主尺度為2．5m×0．5m ×0．175m，規(guī)則波波長與船長比為 λ/L=0．5，速度傅汝德數(shù)為 Fr=0．2，圖10 給出了船模在規(guī)則波中頂浪航行情形下產(chǎn)生的繞射波波形的模型試驗［16］和數(shù)值計算結(jié)果，其中波形的余弦和正弦部分分別與瞬時非定常波 cos(ωet)和sin(ωet)的幅值分量相對應(yīng)。由圖可知:數(shù)值預(yù)報的繞射波形在測試區(qū)域內(nèi)與試驗結(jié)果基本吻合，這表明格林函數(shù)的對稱性適用于船體周圍非定常波的數(shù)值模擬;另外，繞射波的余弦和正弦部分均限制在關(guān)于船舶航行方向?qū)ΨQ的楔形區(qū)域內(nèi)，這與格林函數(shù)的傳播模態(tài)密切相關(guān)，因為移動脈動源的遠(yuǎn)場波形包含計及航速效應(yīng)的環(huán)形波系和計及搖蕩運動下的開爾文波系兩部分。

4 結(jié)論

基于Bessho型移動脈動源格林函數(shù)的對稱特性，本文重點探討分析了波浪中有航速船舶水動力的快速求解措施，并獲得了以下主要結(jié)論:

1)Bessho型移動脈動源格林函數(shù)中的調(diào)和函數(shù)及其ox、oz方向的偏導(dǎo)數(shù)關(guān)于o-xz平面對稱，調(diào)和函數(shù)在oy方向的偏導(dǎo)數(shù)關(guān)于o-xz平面反對稱;結(jié)合船體對稱物面處的法向量取值特性，在關(guān)于o-xz平面對稱的線元和面元處，推導(dǎo)獲得了調(diào)和函數(shù)的法向偏導(dǎo)數(shù)積分的對稱關(guān)系表達(dá)式。

2)將線積分和面積分的積分對稱性應(yīng)用于邊界元積分方程之中，可有效降低線性方程組中系數(shù)矩陣的求解規(guī)模;在波浪中的三維頻域有航速船舶水動力問題，應(yīng)用格林函數(shù)的對稱性可顯著降低數(shù)值計算的運算量和減小變量的存儲量，并可指導(dǎo)船體物面的網(wǎng)格劃分。

3)算例分析表明移動脈動源格林函數(shù)的對稱性適用于求解有航速船舶在波浪中的水動力問題，且驗證了所編制的相關(guān)計算程序的可靠性。

References)

［1］ Maury C，Delhommeau G，Boin J P，et al．Comparison between numerical computations and experiment for seakeeping on ship models with forward speed［J］．Journal of Ship Research，2003，47(4):347-364．

［2］ Inoue Y， Kamruzzaman M． Analysis of hydrodynamic characteristics for arbitrary multihull ships advancing in waves［J］．Journal of Marine Science and Technology，2008，13(3):231-243．

［3］ Xu Y，Dong W C．A fast seakeeping prediction method for a high speed catamaran advancing in waves［J］．Journal of Ship Mechanics，2014，18(6):644-661．

［4］ Xu Y，Dong W C．Study on characteristics of 3-D translatingpulsating source Green function of deep-water Havelock form and its fast integration method［J］．China Ocean Engineering，2011，25(3):365-380．

［5］ 洪亮，朱仁傳，繆國平，等．三維頻域有航速格林函數(shù)的數(shù)值計算與分析［J］．水動力學(xué)研究與進(jìn)展A輯，2013，28(4):423-430．HONG Liang， ZHU Renchuan， MIAO Guoping， et al．Numerical calculation and analysis of3-DGreen's function with forward speed in frequency domain［J］． Journal of Hydrodynamics(Ser．A)，2013，28(4):423-430．(in Chinese)

［6］ Yao C B，DongW C．A fast integrationmethod for translatingpulsating source Green's function in Bessho form［J］．Journal of Zhejiang University(Science A:Applied Physics＆Engineering)，2014，15(2):108-119．

［7］ Wu Y S．Hydroelasticity of floating bodies［D］．Middlesex:Brunel University，1984．

［8］ 劉日明．基于B樣條面元法的浮體二階水動力計算［D］．哈爾濱:哈爾濱工程大學(xué)，2009．LIU Riming．The second order hydrodynamic calculation of floating bodies based on B-spline panelmethod［D］．Harbin:Harbin Engineering University，2009．(in Chinese)

［9］ 王興剛．深海浮式結(jié)構(gòu)物與其系泊纜索的耦合動力分析［D］．大連:大連理工大學(xué)，2011．WANG Xinggang．Coupled dynamic analysis of deep-water floating structure and themooring system［D］．Dalian:Dalian University of Technology，2011．(in Chinese)

［10］ 劉應(yīng)中，繆國平．船舶在波浪上的運動性能［M］．上海:上海交通大學(xué)出版社，1986．LIU Yingzhong，MIAOGuoping．Theoriesof ship'smotions in waves［M］．Shanghai:Publishing House of Shanghai Jiao Tong University，1986．(in Chinese)

［11］ Bessho M．On the fundamental singularity in the theory of ship motion in a seaway［C］//Proceedings of Memoirs of the Defense Academy of Japan，Tokyo，1977，17(8):95-105．

［12］ 張國華，袁乃昌，莊釗文．基于面元法的航母雷達(dá)散射截面計算［J］．國防科技大學(xué)學(xué)報，2001，23(5):79-83．ZHANG Guohua， YUAN Naichang， ZHUANG Zhaowen．Calculation the radar cross section of aaircraftcarrier based on the plate-elementmethod［J］．Journal of National University of Defense Technology，2001，23(5):79-83． (in Chinese)

［13］ 仝佩．基于水彈性力學(xué)理論的船體動態(tài)形變研究［D］．長沙:國防科學(xué)技術(shù)大學(xué)，2011．TONG Pei．Research of dynamic distortion on a ship based on the theory of hydroelasticity［D］． Changsha:National University of Defense Technology，2011．(in Chinese)

［14］ Molland A F，Wellicome JF，Temarel P，etal．Experimental investigation of the seakeeping characteristics of fast displacement catamarans in head and oblique waves［J］．Transactions of the Royal Institution of Naval Architects，2001，143:78-97．

［15］ Bailey P A，Hudson D A，Price W G，et al．Comparisons between theory and experiment in a seakeeping validation study［J］．Transactions of the Royal Institution of Naval Architects，2001，143:44-77．

［16］ Iwashita H，Kashiwagi M，Elangovan M，et al．On an unsteady wave pattern analysis of ships advancing in waves［J］．Journal of the Japan Society of Naval Architects and Ocean Engineers，2011，6(13)，95-106．