基于橫搖衰減數(shù)據(jù)的橫搖阻尼估算方法對比研究

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

（1.中國船舶科學(xué)研究中心，江蘇無錫 214082；2.水動力學(xué)重點實驗室，江蘇無錫 214082）

0 引 言

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

從橫搖衰減曲線估算橫搖阻尼的方法眾多，自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è)計研究院封培元等[17]的改進型消滅曲線法、中國海洋大學(xué)孫金偉等[18]的Prony-SS 方法、隨機減量法[19]等方法均成功用于橫搖阻尼估算。其中隨機減量法是從波浪中船舶橫搖運動數(shù)據(jù)中提取橫搖衰減運動的方法，該方法能從非規(guī)則波和規(guī)則波船舶橫搖運動中提取橫搖衰減曲線。

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

1 橫搖運動數(shù)學(xué)模型

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

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

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

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

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

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

2.1 消滅曲線法

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

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

2.2 改進型消滅曲線法

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

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

2.3 對數(shù)衰減法

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

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

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

2.4 Hilbert變換法

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

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

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

2.5 Froude能量法

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

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

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

2.6 Roberts能量法

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

橫搖運動到幅值φi時，角速度為0，式（36）右側(cè)第一項為0。

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

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

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

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

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

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

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

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

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

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

3.2 Prony-SS方法與漸近方法

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

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

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

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

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

圖2 DTC標模剖面圖[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標模參數(shù)[20]Tab.2 Calculation parameters of the DTC model[20]

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

圖5（a）為Wassermann[11]對圖3（a）所示的橫搖衰減曲線先進行濾波處理后再基于Froude能量法獲得的無量綱阻尼系數(shù)結(jié)果，圖5（b）為作者對圖3（a）所示的橫搖衰減曲線直接使用Froude 能量法后獲得的無量綱阻尼系數(shù)結(jié)果，二者的大振幅結(jié)果基本一致，小振幅結(jié)果振蕩大，這與Wassermann 對數(shù)據(jù)進行了濾波[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展示了壓載工況下各方法（為方便比較對方法進行了編號）的結(jié)果，表4 展示了滿載工況下非線性阻尼系數(shù)的擬合結(jié)果。圖6～8中數(shù)據(jù)點振蕩均較大，這與圖3（a）的數(shù)據(jù)尤其是相鄰正負幅值的振蕩密切相關(guān)，說明實驗數(shù)據(jù)質(zhì)量對基于該類方法的阻尼系數(shù)估算有影響。

圖6 采用式（15）擬合的消滅曲線與采用式（19）擬合的改進型消滅曲線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 幅值與頻率隨時間的關(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ù)識別的橫搖阻尼估算方法結(jié)果

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

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

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

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

5 橫搖阻尼估算方法特點對比

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

5.1 橫搖阻尼評估效果

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

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

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

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

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

基于分段線性假設(shè)的估算方法涉及的曲線擬合、求導(dǎo)過程都會產(chǎn)生誤差，表6所示為各方法曲線擬合與求導(dǎo)狀況。對于消滅曲線法與改進型消滅曲線法，分段線性假設(shè)后，為了能實現(xiàn)消滅曲線擬合（如式（12）），一般假設(shè)所有周期頻率都為橫搖自然頻率，然而，如圖9（b）所示，真實的橫搖頻率會隨時間不斷變化，即式（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ù)識別的估算方法涉及的阻尼模型選擇同樣會引入誤差，根據(jù)最優(yōu)化理論，非線性最小二乘法、粒子群對測量誤差與噪聲敏感且尋優(yōu)結(jié)果可能為局部最小，即不一定是真實結(jié)果，且參數(shù)設(shè)置對結(jié)果影響大（如表5 所示）。Prony-SS 方法重構(gòu)曲線后，求解阻尼參數(shù)的線性方程組為超定方程組，解該類方程組用到線性最小二乘法，存在唯一解，不存在該問題，但只有阻尼模型足夠精確其參數(shù)識別結(jié)果才可信（作者使用自定義阻尼參數(shù)生成衰減曲線后，使用Prony-SS方法驗證后得出此結(jié)論）。

6 結(jié) 論

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

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

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

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

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