基于橫搖衰減數據的橫搖阻尼估算方法對比研究

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

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

0 引 言

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

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

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

1 橫搖運動數學模型

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

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

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

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

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

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

2.1 消滅曲線法

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

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

2.2 改進型消滅曲線法

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

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

2.3 對數衰減法

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

根據線性系統振動理論易得

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

2.4 Hilbert變換法

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

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

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

2.5 Froude能量法

根據系統做功能量關系，相鄰橫搖幅值間阻尼做功應等于系統勢能的變化量。當從橫搖幅值φi運動到φi+2時,系統勢能變化量與阻尼耗散的能量分別為

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

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

2.6 Roberts能量法

系統機械能EC的耗散率應等于阻尼做功的功率。

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

根據式（37）可求得幅值φi處的EC值，擬合EC與時間關系曲線，再求導并結合式（22）可得

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

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

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

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

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

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

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

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

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

3.2 Prony-SS方法與漸近方法

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

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

4 橫搖阻尼估算-研究結果與討論

4.1 研究對象及數據簡介

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

圖2 DTC標模剖面圖[20]Fig.2 Hull sections of DTC[20]

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

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

表2 DTC標模參數[20]Tab.2 Calculation parameters of the DTC model[20]

4.2 基于分段線性假設的橫搖阻尼估算結果

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

圖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展示了壓載工況下各方法（為方便比較對方法進行了編號）的結果，表4 展示了滿載工況下非線性阻尼系數的擬合結果。圖6～8中數據點振蕩均較大，這與圖3（a）的數據尤其是相鄰正負幅值的振蕩密切相關，說明實驗數據質量對基于該類方法的阻尼系數估算有影響。

圖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 幅值與頻率隨時間的關系曲線擬合-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 非線性阻尼系數擬合結果，壓載（D=0.2018 m），無航速Tab.3 Results of nonlinear damping coefficients fitting,ballast condition(D=0.2018 m)

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

4.3 基于參數識別的橫搖阻尼估算方法結果

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

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

圖10 基于Prony-SS方法的橫搖衰減曲線重構，壓載（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 非線性阻尼系數參數識別結果Tab.5 Results of nonlinear damping coefficients based on PIT

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

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

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

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

5.1 橫搖阻尼評估效果

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

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

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

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

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

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

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

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

6 結 論

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

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

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

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

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