章 東，王文濤，卜淑霞，劉 偉

（1.中國船舶科學(xué)研究中心，江蘇無錫 214082；2.水動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，江蘇無錫 214082）

0 引 言

橫搖是船舶在波浪中6個(gè)自由度運(yùn)動(dòng)中最難理解和顯示的復(fù)雜現(xiàn)象[1]，其運(yùn)動(dòng)阻尼是船舶耐波性與波浪穩(wěn)性預(yù)測精度的決定性因素。一般來說，船舶運(yùn)動(dòng)阻尼主要由船體輻射出的重力波所致，然而相比于其他自由度，船舶橫搖受流體粘性影響大。橫搖阻尼是由波浪、旋渦、升力、摩擦等因素共同產(chǎn)生，這導(dǎo)致基于勢流理論的邊界元方法（BEM）預(yù)測橫搖運(yùn)動(dòng)精度過低。針對此難題，ITTC（國際拖曳水池會(huì)議）的波浪穩(wěn)性委員會(huì)與耐波性委員會(huì)開展了大量工作[2-5]。目前橫搖阻尼估算方法主要包括三類：模型試驗(yàn)、CFD計(jì)算、經(jīng)驗(yàn)公式。中國船舶科學(xué)研究中心顧民等[6-8]將模型試驗(yàn)與CFD計(jì)算結(jié)合，對一系列標(biāo)準(zhǔn)模型做了深入研究。模型試驗(yàn)與CFD 計(jì)算結(jié)果都需要通過數(shù)據(jù)處理以獲得阻尼值，其中自由衰減試驗(yàn)具有試驗(yàn)難度小、成本低的特點(diǎn)，該方法被廣泛采用。

從橫搖衰減曲線估算橫搖阻尼的方法眾多，自1872 年Froude[9]的工作開始，人們在該領(lǐng)域做了大量工作。目前除經(jīng)典的線性方法（消滅曲線法）、Froude能量法[1]外，意大利的Bulian 等[10]、德國的Wassermann 等[11]在對數(shù)衰減法方面做了深入研究，英國的Roberts[12]提出的能量方法、國際海事組織推薦的最小二乘法[13]、上海交通大學(xué)曾智華[14]的粒子群方法、香港城市大學(xué)Chan 等[15]的漸近方法、韓國的Kim 等[16]的Hilbert 變換方法、中國船舶及海洋工程設(shè)計(jì)研究院封培元等[17]的改進(jìn)型消滅曲線法、中國海洋大學(xué)孫金偉等[18]的Prony-SS 方法、隨機(jī)減量法[19]等方法均成功用于橫搖阻尼估算。其中隨機(jī)減量法是從波浪中船舶橫搖運(yùn)動(dòng)數(shù)據(jù)中提取橫搖衰減運(yùn)動(dòng)的方法，該方法能從非規(guī)則波和規(guī)則波船舶橫搖運(yùn)動(dòng)中提取橫搖衰減曲線。

船舶橫搖運(yùn)動(dòng)特別是大角度橫搖運(yùn)動(dòng)非線性強(qiáng)，基于此特點(diǎn)，從橫搖衰減曲線估算橫搖阻尼的方法主要分為兩類：一是基于分段線性假設(shè)，將橫搖衰減曲線分成多段，分別處理各小段擬線性運(yùn)動(dòng)的阻尼關(guān)系；二是從橫搖運(yùn)動(dòng)微分方程出發(fā)，先假設(shè)阻尼模型及復(fù)原力臂模型（未知時(shí)），后通過漸近方法、數(shù)值方法等手段經(jīng)計(jì)算機(jī)迭代識(shí)別出阻尼模型參數(shù)，該類方法稱為參數(shù)識(shí)別技術(shù)（parameter identification technique，PIT）。本文的目的是針對這兩類方法進(jìn)行分析，以獲取高精度阻尼預(yù)報(bào)方法。

1 橫搖運(yùn)動(dòng)數(shù)學(xué)模型

橫搖運(yùn)動(dòng)時(shí)，粘性項(xiàng)與波浪項(xiàng)產(chǎn)生的阻尼量級相當(dāng)，直接求解船體周圍流體速度場以預(yù)報(bào)船舶橫搖運(yùn)動(dòng)難度大且精度低，現(xiàn)階段一般通過彈簧阻尼滑塊模型來模化橫搖運(yùn)動(dòng)。解耦并歸一化單自由度橫搖運(yùn)動(dòng)方程[3]如下：

表1 橫搖阻尼估算方法分類Tab.1 Classification of roll damping estimation methods

2 基于分段線性假設(shè)的橫搖阻尼估算方法

該類方法的基本思想是分段研究每個(gè)周期內(nèi)橫搖，假設(shè)每單個(gè)周期內(nèi)的運(yùn)動(dòng)均為正弦（或余弦）運(yùn)動(dòng)，再按如圖1所示流程估算橫搖阻尼。除消滅曲線法與改進(jìn)型消滅曲線法外，基于分段線性假設(shè)的其他方法，都要引入一個(gè)重要概念：等效線性阻尼系數(shù)Beq及其歸一化形式μeq，某一時(shí)間段內(nèi)式（1）中的阻尼項(xiàng)d(φ? )表示如下[13]：

圖1 基于分段線性假設(shè)的橫搖阻尼估算流程圖Fig.1 Procedure of roll damping estimation based on piecewise linear assumption

式中，μeq為該時(shí)間段內(nèi)等效阻尼系數(shù)，φa為該周期內(nèi)橫搖角幅值，ωE為該周期內(nèi)橫搖圓頻率。獲得等效阻尼μeq后，由能量關(guān)系可推導(dǎo)出非線性阻尼模型（如式（3））系數(shù)與等效阻尼的關(guān)系[10]：

2.1 消滅曲線法

假定橫搖復(fù)原力臂——GZ( )φ為線性的情況下，橫搖衰減運(yùn)動(dòng)由下式控制：

需要在時(shí)歷曲線的三個(gè)不同位置分別取值，獲得三個(gè)三元一次方程以求解式（15）中的三個(gè)非線性阻尼模型系數(shù)b1、b2、b3。

2.2 改進(jìn)型消滅曲線法

封培元等[17]提出了復(fù)原力臂——GZ( )φ為3階模型、阻尼模型為2項(xiàng)情況下的改進(jìn)型消滅曲線法：

式（19）即擴(kuò)展的改進(jìn)型消滅曲線關(guān)系，其系數(shù)ε1、ε2、b1、b2、b3均可為0，與文獻(xiàn)[17]相比能適應(yīng)更多復(fù)原力臂模型與阻尼模型。

2.3 對數(shù)衰減法

對數(shù)衰減法（Bulian）[10]將相鄰峰值時(shí)間段[ti+1-ti]內(nèi)橫搖衰減運(yùn)動(dòng)類比為二階線性系統(tǒng)振動(dòng)過程：

根據(jù)線性系統(tǒng)振動(dòng)理論易得

若式（25）～（26）中無阻尼自然橫搖頻率ωx未知，可利用式（27）代入橫穩(wěn)性高-- --GM與船寬BWL估算，

2.4 Hilbert變換法

Hilbert變換相當(dāng)于90°移相[16]，記橫搖衰減信號φ(t)的Hilbert變換為φ?(t)，則可以構(gòu)造復(fù)平面內(nèi)解析信號:

式中，μeq(t)為隨時(shí)間變化的等效阻尼項(xiàng)，ω2x(t)為隨時(shí)間變化的橫搖自然頻率項(xiàng)。將式（28）代入式（30）可得以下關(guān)系：

根據(jù)式（31）與式（32），只要從橫搖衰減曲線中擬合出A(t)與ωE( )t，即可獲得等效阻尼隨時(shí)間的變化關(guān)系，進(jìn)而得到每個(gè)橫搖周期的等效阻尼μeq。

2.5 Froude能量法

根據(jù)系統(tǒng)做功能量關(guān)系，相鄰橫搖幅值間阻尼做功應(yīng)等于系統(tǒng)勢能的變化量。當(dāng)從橫搖幅值φi運(yùn)動(dòng)到φi+2時(shí),系統(tǒng)勢能變化量與阻尼耗散的能量分別為

令式（33）與式（34）相等，可得

歸一化的等效阻尼系數(shù)μeq的估算同式（26）。

2.6 Roberts能量法

系統(tǒng)機(jī)械能EC的耗散率應(yīng)等于阻尼做功的功率。

橫搖運(yùn)動(dòng)到幅值φi時(shí)，角速度為0，式（36）右側(cè)第一項(xiàng)為0。

根據(jù)式（37）可求得幅值φi處的EC值，擬合EC與時(shí)間關(guān)系曲線，再求導(dǎo)并結(jié)合式（22）可得

歸一化的等效阻尼系數(shù)μeq的估算同式（26）。

3 基于參數(shù)識(shí)別的橫搖阻尼估算方法

該類方法需預(yù)設(shè)未知項(xiàng)的模型，一般根據(jù)經(jīng)驗(yàn)選擇合適的阻尼和復(fù)原力臂模型，得到有待定參數(shù)的二階常微分方程，再通過數(shù)學(xué)手段尋找最接近實(shí)驗(yàn)數(shù)據(jù)的解的模型參數(shù)。

3.1 最小二乘法與粒子群方法

以采用式（9）所示阻尼模型為例。

步驟一：先猜測一組初始參數(shù)如b1、b2、b3。

步驟二：利用數(shù)值方法(如Runge-Kutta法)或解析方法求解橫搖運(yùn)動(dòng)方程，生成解曲線。

步驟三：定義優(yōu)化目標(biāo)函數(shù)，如式（39），利用最優(yōu)化方法（非線性最小二乘優(yōu)化中的Levenberg-Marquardt優(yōu)化、粒子群優(yōu)化[14]等方法均被成功應(yīng)用于橫搖阻尼識(shí)別）尋找使得該式最小的一組參數(shù)：

式中，φi為計(jì)算所得橫搖角數(shù)據(jù)，φexp,i為試驗(yàn)橫搖角數(shù)據(jù)，既可采用整條衰減曲線數(shù)據(jù)也可只選用幅值數(shù)據(jù)。為減少計(jì)算量，一般采用幅值數(shù)據(jù)。

3.2 Prony-SS方法與漸近方法

Prony-SS方法最先應(yīng)用于信號處理領(lǐng)域，對于給定橫搖衰減曲線，利用Prony-SS方法重構(gòu)出衰減曲線的近似表達(dá)式[18]，后通過矩陣運(yùn)算獲得阻尼參數(shù)值。Prony-SS 方法能提取出橫搖信號中的主要成分，具有內(nèi)在的噪聲抑制機(jī)制。

漸近方法是求解非線性微分方程的有力手段，對于利用衰減曲線估算橫搖阻尼，該方法[15]先利用攝動(dòng)方法求解出橫搖衰減運(yùn)動(dòng)方程的漸近解表達(dá)式，再通過曲線擬合獲得表達(dá)式系數(shù)，利用消滅曲線反復(fù)迭代優(yōu)化阻尼系數(shù)值。

4 橫搖阻尼估算-研究結(jié)果與討論

4.1 研究對象及數(shù)據(jù)簡介

為驗(yàn)證并比較上述方法的特點(diǎn)，本文采用標(biāo)模DTC[20]的橫搖衰減試驗(yàn)數(shù)據(jù)估算橫搖阻尼，DTC 標(biāo)模是德國杜伊斯堡埃森大學(xué)按照超巴拿馬型集裝箱船體以1∶59.407 的比例設(shè)計(jì)建造的，船體剖面如圖2 所示。船模數(shù)據(jù)如表2 所示，采用壓載（ballast）、滿載（full condition）兩種工況下的無航速自由橫搖衰減數(shù)據(jù)（如圖3 所示）作為樣本，其復(fù)原力臂曲線如圖4 所示，從圖中可以看出，在試驗(yàn)橫搖數(shù)據(jù)范圍內(nèi)，復(fù)原力臂曲線擬合效果良好。采用Wassermann[11]文中推薦的Froude 能量法復(fù)現(xiàn)其阻尼估算結(jié)果后，利用本文中所述方法處理該數(shù)據(jù)并進(jìn)行結(jié)果對比。

圖2 DTC標(biāo)模剖面圖[20]Fig.2 Hull sections of DTC[20]

圖3 橫搖衰減曲線Fig.3 Roll decay curve[20]

圖4 采用5階級數(shù)擬合的——GZ曲線[11]Fig.4 Fitting curve of——GZ based on 5th oder series[11]

表2 DTC標(biāo)模參數(shù)[20]Tab.2 Calculation parameters of the DTC model[20]

4.2 基于分段線性假設(shè)的橫搖阻尼估算結(jié)果

圖5（a）為Wassermann[11]對圖3（a）所示的橫搖衰減曲線先進(jìn)行濾波處理后再基于Froude能量法獲得的無量綱阻尼系數(shù)結(jié)果，圖5（b）為作者對圖3（a）所示的橫搖衰減曲線直接使用Froude 能量法后獲得的無量綱阻尼系數(shù)結(jié)果，二者的大振幅結(jié)果基本一致，小振幅結(jié)果振蕩大，這與Wassermann 對數(shù)據(jù)進(jìn)行了濾波[11]操作有關(guān)。圖中無量綱阻尼系數(shù)由下式得到：

圖5 不同橫搖幅值下的無量綱阻尼（壓載（D=0.2018 m），無航速）Fig.5 Non-dimensional roll damping with different amplitudes,ballast condition(D=0.2018 m)，without forward speed

按照圖1所示流程分別處理圖3所示的橫搖衰減曲線，阻尼模型采用式（9）。表3及圖6～9展示了壓載工況下各方法（為方便比較對方法進(jìn)行了編號）的結(jié)果，表4 展示了滿載工況下非線性阻尼系數(shù)的擬合結(jié)果。圖6～8中數(shù)據(jù)點(diǎn)振蕩均較大，這與圖3（a）的數(shù)據(jù)尤其是相鄰正負(fù)幅值的振蕩密切相關(guān)，說明實(shí)驗(yàn)數(shù)據(jù)質(zhì)量對基于該類方法的阻尼系數(shù)估算有影響。

圖6 采用式（15）擬合的消滅曲線與采用式（19）擬合的改進(jìn)型消滅曲線Fig.6 Fitting curve of extinction curve method based on Eq.(15)and enhanced extinction curve method based on Eq.(19),ballast condition(D=0.2010 m),without forward speed

圖7 基于式（7）的等效阻尼擬合，壓載（D=0.2018 m），無航速Fig.7 Curve fitting of equivalent linear damping coefficients based on Eq.(7),ballast condition(D=0.2018 m),without forward speed

圖8 基于式（7）的等效阻尼擬合，壓載（D=0.2018 m），無航速Fig.8 Fitting curve of equivalent linear damping coefficients based on Eq.(7),ballast condition(D=0.2018 m),without forward speed

圖9 幅值與頻率隨時(shí)間的關(guān)系曲線擬合-Hilbert變換方法，壓載（D=0.2018 m），無航速Fig.9 Fitting curve of amplitudes and frequency with time-Hilbert transform method,ballast codition(D=0.2018 m),without forward speed

表3 非線性阻尼系數(shù)擬合結(jié)果，壓載（D=0.2018 m），無航速Tab.3 Results of nonlinear damping coefficients fitting,ballast condition(D=0.2018 m)

表4 非線性阻尼系數(shù)擬合結(jié)果，滿載（D=0.2354 m），無航速Tab.4 Results of nonlinear damping coefficients fitting,full loading condition(D=0.2354 m)

4.3 基于參數(shù)識(shí)別的橫搖阻尼估算方法結(jié)果

利用參數(shù)識(shí)別方法（PIT）估算橫搖阻尼需要先預(yù)設(shè)阻尼模型和復(fù)原力模型（未知時(shí)），為便于對比分析，本文中復(fù)原力采用圖4 所示的擬合模型，阻尼模型選用式（9）形式。所得歸一化的自由橫搖運(yùn)動(dòng)微分方程形式如下：

根據(jù)圖4 擬合結(jié)果，式（41）中ε1、ε2已知，ωx使用式（27）估算。根據(jù)文獻(xiàn)[11]可知，式（27）估算結(jié)果偏差不超過1%,則上式中未知參數(shù)有b1、b2、b3，本文選取最小二乘法、粒子群算法和Prony-SS 方法（為方便比較對方法進(jìn)行了編號）識(shí)別上述參數(shù)。圖10和表5展示了部分識(shí)別過程及非線性阻尼系數(shù)識(shí)別結(jié)果。

圖10 基于Prony-SS方法的橫搖衰減曲線重構(gòu)，壓載（D=0.2018 m），無航速Fig.10 Reconstruction of roll decay curves based on Prony-SS method,ballast condition(D=0.2018 m),without forward specd

表5 非線性阻尼系數(shù)參數(shù)識(shí)別結(jié)果Tab.5 Results of nonlinear damping coefficients based on PIT

最小二乘法參數(shù)初值采用[0；0；0]，粒子群算法的參數(shù)采用默認(rèn)設(shè)置。圖10（a）展示了使用Prony-SS 方法時(shí)對橫搖衰減數(shù)據(jù)構(gòu)成的Hankel 矩陣進(jìn)行奇異值分解（singular value decomposition）后獲得的奇異值量級分布。根據(jù)奇異值量級大小可判斷重構(gòu)的衰減曲線近似表達(dá)式的主要項(xiàng)，圖10（b）為取前11個(gè)奇異值對應(yīng)的項(xiàng)重構(gòu)出的近似表達(dá)式與實(shí)驗(yàn)數(shù)據(jù)的對比，可見重構(gòu)精度良好。

表5 中的非線性阻尼系數(shù)識(shí)別結(jié)果為未經(jīng)調(diào)試檢驗(yàn)的結(jié)果，均與基于分段線性假設(shè)的方法識(shí)別結(jié)果存在較大差異，說明基于參數(shù)識(shí)別的方法可能存在多解的現(xiàn)象。該類方法結(jié)果受識(shí)別過程中參數(shù)設(shè)置的影響，為排除多解和得到足夠精度的阻尼系數(shù)，既要經(jīng)驗(yàn)也需要根據(jù)每條實(shí)驗(yàn)數(shù)據(jù)反復(fù)比對調(diào)試。

5 橫搖阻尼估算方法特點(diǎn)對比

橫搖阻尼評估是船舶橫搖運(yùn)動(dòng)預(yù)報(bào)簡化的重要步驟，由于流動(dòng)的復(fù)雜性、船舶各自由度運(yùn)動(dòng)的耦合，式（1）形式的運(yùn)動(dòng)方程只是對真實(shí)橫搖運(yùn)動(dòng)的近似。好的評估結(jié)果應(yīng)能提取出橫搖阻尼的內(nèi)在特性，即真實(shí)的橫搖運(yùn)動(dòng)應(yīng)為該特性疊加上各種隨機(jī)因素的耦合結(jié)果。從表3～5 的阻尼參數(shù)擬合結(jié)果來看，基于分段線性假設(shè)的方法一致性較好，基于參數(shù)識(shí)別的方法結(jié)果差別較大，顯示橫搖阻尼估算結(jié)果準(zhǔn)確度與使用的方法強(qiáng)相關(guān)。

5.1 橫搖阻尼評估效果

為評估各方法的阻尼估算效果，參照式（39）定義誤差R，表達(dá)式如下：

式中，φi為評估得到的阻尼模型參數(shù)代入橫搖運(yùn)動(dòng)微分方程后，使用四階Runge-Kutta法求得的橫搖曲線衰減數(shù)據(jù)（本文取橫搖幅值數(shù)據(jù)），φexp,i為對應(yīng)實(shí)驗(yàn)數(shù)據(jù)。各方法的R值見圖11，可見R值均不為0。R不為0 有三方面原因：存在測量誤差及噪聲，阻尼模型不夠精確，方法自身引入誤差。

圖11 DTC標(biāo)模不同數(shù)據(jù)處理方法的誤差RFig.11 Comparison of roll prediction errors by roll damping of DTC model from different estimation methods

5.2 基于自由橫搖衰減曲線的橫搖阻尼評估方法誤差來源

對于從自由衰減曲線估算橫搖阻尼而言，有效橫搖周期越多越有利于獲取準(zhǔn)確的橫搖阻尼信息。然而，一般小振幅橫搖測量噪聲與誤差占比大，大幅度橫搖運(yùn)動(dòng)衰減太快，有效橫搖周期少，如圖5所示。為盡可能提高橫搖阻尼評估精度，需要選擇對測量噪聲與誤差不敏感且評估過程誤差引入少的方法。

基于分段線性假設(shè)的估算方法涉及的曲線擬合、求導(dǎo)過程都會(huì)產(chǎn)生誤差，表6所示為各方法曲線擬合與求導(dǎo)狀況。對于消滅曲線法與改進(jìn)型消滅曲線法，分段線性假設(shè)后，為了能實(shí)現(xiàn)消滅曲線擬合（如式（12）），一般假設(shè)所有周期頻率都為橫搖自然頻率，然而，如圖9（b）所示，真實(shí)的橫搖頻率會(huì)隨時(shí)間不斷變化，即式（15）、式（16）與式（19）都是近似等式，引入了額外誤差。

表6 基于分段線性假設(shè)的估算方法涉及的曲線擬合及求導(dǎo)次數(shù)Tab.6 Times of curve fitting and derivation of estimation methods based on piecewise linear assumption

基于參數(shù)識(shí)別的估算方法涉及的阻尼模型選擇同樣會(huì)引入誤差，根據(jù)最優(yōu)化理論，非線性最小二乘法、粒子群對測量誤差與噪聲敏感且尋優(yōu)結(jié)果可能為局部最小，即不一定是真實(shí)結(jié)果，且參數(shù)設(shè)置對結(jié)果影響大（如表5 所示）。Prony-SS 方法重構(gòu)曲線后，求解阻尼參數(shù)的線性方程組為超定方程組，解該類方程組用到線性最小二乘法，存在唯一解，不存在該問題，但只有阻尼模型足夠精確其參數(shù)識(shí)別結(jié)果才可信（作者使用自定義阻尼參數(shù)生成衰減曲線后，使用Prony-SS方法驗(yàn)證后得出此結(jié)論）。

6 結(jié) 論

本文以獲取高精度阻尼預(yù)報(bào)方法為目標(biāo)，對9種（漸近方法未考慮，其阻尼估算過程復(fù)雜，自身誤差項(xiàng)多）橫搖阻尼評估方法的原理進(jìn)行分析并利用DTC 標(biāo)模橫搖衰減數(shù)據(jù)進(jìn)行驗(yàn)證后，得出以下結(jié)論：

（1）基于分段線性假設(shè)的橫搖阻尼估算方法物理意義明確，不需預(yù)設(shè)阻尼模型，但涉及的曲線擬合與求導(dǎo)過程易引入額外誤差；除消滅曲線法外，其他方法均能將非線性橫搖力臂納入考慮；真實(shí)的橫搖頻率隨時(shí)間不斷變化，消滅曲線關(guān)系式（15）、改進(jìn)型消滅曲線關(guān)系式（16）和式（19）只是近似成立；Froude能量法對測量噪聲與誤差不敏感，自身引入誤差最少，具有良好性能。

（2）基于參數(shù)識(shí)別的橫搖阻尼估算方法直接從微分方程與實(shí)現(xiàn)數(shù)據(jù)出發(fā)，原理清晰，但需預(yù)設(shè)阻尼模型且識(shí)別過程為黑箱子；非線性最小二乘優(yōu)化、粒子群優(yōu)化迭代參數(shù)對尋優(yōu)結(jié)果影響大即可能為局部最??；Prony-SS方法存在唯一解，但只有阻尼模型足夠精確其參數(shù)識(shí)別結(jié)果才可信。

（3）推薦使用對測量噪聲與誤差不敏感[11]，只有一次曲線擬合、無求導(dǎo)、物理意義明確、識(shí)別過程簡潔的Froude能量法擬合出等效阻尼曲線，若能找出足夠精確的阻尼模型，再利用Prony-SS方法獲得更精確的阻尼模型參數(shù)。

橫搖阻尼估算是一項(xiàng)復(fù)雜的工作，本文主要研究了基于衰減曲線估算橫搖阻尼的方法，并未對不同影響因素下（航速、船型參數(shù)、尺度效應(yīng)、流動(dòng)記憶效應(yīng)等）獲得的自由衰減曲線是否能提供足夠全的信息、是否反映了橫搖的本質(zhì)進(jìn)行討論，因此使用以上方法估算橫搖阻尼時(shí)應(yīng)予以考慮。