孫金偉，邵 萌，邵珠曉，李選群，萬曉正，趙環(huán)宇，劉海豐

（1.中國(guó)海洋大學(xué)工程學(xué)院，青島266100；2.山東省海洋工程重點(diǎn)實(shí)驗(yàn)室，青島266100；3.齊魯工業(yè)大學(xué)（山東省科學(xué)院）海洋儀器儀表研究所，青島266001；4.山東省海洋環(huán)境監(jiān)測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，青島266001）

0 引 言

船舶橫搖阻尼的準(zhǔn)確估算是正確預(yù)報(bào)波浪中橫搖運(yùn)動(dòng)的前提。小角度橫搖時(shí)，船舶橫搖運(yùn)動(dòng)可以用線性方程描述，但當(dāng)橫搖幅值增大時(shí)，必須考慮非線性效應(yīng)。橫搖非線性主要表現(xiàn)在回復(fù)力矩的非線性和阻尼的非線性，很多學(xué)者對(duì)非線性的阻尼和回復(fù)力矩?cái)?shù)學(xué)模型開展了研究[1-4]。

船舶橫搖回復(fù)力矩可由不同階次的奇次多項(xiàng)式函數(shù)表示，回復(fù)力矩系數(shù)可通過與流體靜力學(xué)計(jì)算得到GZ曲線進(jìn)行擬合獲得[5]。而對(duì)于橫搖阻尼項(xiàng)，目前尚未有完備的理論計(jì)算方法，現(xiàn)有的船舶耐波性評(píng)估的勢(shì)流理論不能預(yù)報(bào)橫搖運(yùn)動(dòng)中由于摩擦、旋渦和流動(dòng)分離等因素產(chǎn)生的粘性阻尼[6]。基于模型試驗(yàn)獲得的經(jīng)驗(yàn)公式在實(shí)際使用中較為方便，但是僅對(duì)特定船型適用[7-9]。因此在工程實(shí)踐中常用模型試驗(yàn)的方法來確定船舶橫搖阻尼系數(shù)。

最早利用模型試驗(yàn)進(jìn)行橫搖阻尼系數(shù)估算的學(xué)者是Froude（1872）[10]。Wassermann（2016）等[11]對(duì)不同的模型試驗(yàn)技術(shù)，包括衰減橫搖試驗(yàn)和強(qiáng)迫橫搖試驗(yàn)方法等進(jìn)行了研究總結(jié)。Spouge（1988）[12]對(duì)于不同學(xué)者提出的各類不同阻尼系數(shù)識(shí)別方法和結(jié)果進(jìn)行了總結(jié)比較，并探討了各種方法的識(shí)別精度。Mathisen 和Price（1985）[13]提出了估算衰減橫搖和強(qiáng)迫橫搖試驗(yàn)阻尼系數(shù)的方法，研究發(fā)現(xiàn)線性加平方項(xiàng)阻尼模型比線性加立方項(xiàng)阻尼模型具有更好的優(yōu)越性。Roberts（1985）[5]提出了根據(jù)自由橫搖衰減曲線估算非線性橫搖阻尼項(xiàng)的計(jì)算方法，該方法適用于小阻尼情形，不限制回復(fù)剛度的非線性程度，但遺憾的是作者給出的例子中并沒有考慮強(qiáng)非線性的回復(fù)力矩。Chan 等（1995）[14]將徐兆（1985）[15]提出的一種新的漸進(jìn)法應(yīng)用于船舶非線性橫搖阻尼系數(shù)估算中，該方法適用于大角度橫搖以及強(qiáng)非線性回復(fù)力矩，但是其在不同初始橫搖角度下的阻尼系數(shù)識(shí)別結(jié)果存在一定計(jì)算誤差，識(shí)別精度仍有待提高。李紅霞等（2005）[16]根據(jù)衡量耗散的觀點(diǎn)利用橫搖試驗(yàn)衰減曲線，提出了一種非線性阻尼識(shí)別方法，通過算例分析認(rèn)為二次阻尼模型和三次阻尼模型對(duì)船舶橫搖的作用是等效的。馬新謀等（2013）[17]采用李紅霞等提出的能量法，依據(jù)兩棲車輛自由橫搖衰減曲線確定了兩棲車輛的非線性阻尼項(xiàng)。馬山等（2012）[6]基于模型試驗(yàn)的橫搖衰減測(cè)量數(shù)據(jù)，采用能量法進(jìn)行了船舶橫搖阻尼系數(shù)的研究，并通過實(shí)驗(yàn)證實(shí)了采用能量法獲取的橫搖阻尼系數(shù)結(jié)合勢(shì)流興波阻尼預(yù)報(bào)橫搖運(yùn)動(dòng)具有較高的精度。

本文基于非線性力學(xué)中的漸進(jìn)法，發(fā)展了一種基于橫搖自由衰減數(shù)據(jù)的阻尼系數(shù)識(shí)別方法，通過數(shù)值仿真和對(duì)比研究分析，驗(yàn)證了該方法的高精度和適用性，最后將該方法應(yīng)用到船舶橫搖衰減實(shí)測(cè)數(shù)據(jù)，識(shí)別其非線性阻尼系數(shù)，數(shù)值模擬結(jié)果與試驗(yàn)數(shù)據(jù)表現(xiàn)出了較好的一致性。

1 漸進(jìn)法介紹

對(duì)強(qiáng)非線性系統(tǒng)，其運(yùn)動(dòng)可由下列微分方程描述

式中：ε為小參數(shù)，函數(shù)g( x )和f(x,x)滿足如下關(guān)系：

根據(jù)攝動(dòng)理論[18]，（1）式具有如下形式的廣義漸進(jìn)解

式中：a 是振幅，φ 是相位角，均為時(shí)間的緩變函數(shù)，x1,x2,…,xm-1是相位角φ 的周期函數(shù)，其周期為2π。a和φ由如下微分方程決定：

為簡(jiǎn)化計(jì)算，（1）式的解通常由其一階近似解表示[15]

式中：

因此要求解（6）式，則需要求解（7）式中參數(shù)A1( a )，Φ0( a,φ )和Φ1( a,φ )。具體步驟如下：

令ε = 0，則由（7）式可知：a= 0，Φ0( a,φ )= φ。當(dāng)ε = 0時(shí)，（1）式變?yōu)?/p>

（8）式兩邊同乘以x并積分可得

式中，C是積分常數(shù)，且有

假定初始條件為x(0)= a,x( 0 )= 0，則由（9）式可確定C = V(a)。將（6）式代入（9）式中可得

由（6）-（7）式可得

把函數(shù)g( x )展開成ε的冪級(jí)數(shù)

將（12）和（13）式代入（1）式，使等式兩邊ε 項(xiàng)的系數(shù)相等，得到

令φ = 2π，可以求出A1表達(dá)式

令φ = π，可以求出Φ1表達(dá)式。至此，（7）式中的未知項(xiàng)A1( a )，Φ0( a,φ )和Φ1( a,φ )全部求出，于是方程一階近似解的表達(dá)式（6）可求出。

2 阻尼系數(shù)識(shí)別

船舶橫搖非線性阻尼選用如（17）式所示的線性加平方阻尼模型，回復(fù)力矩用如（18）式所示的奇次多項(xiàng)式函數(shù)表示：

（18）式中的冪次數(shù)項(xiàng)j和回復(fù)力系數(shù)kj可由流體靜力學(xué)計(jì)算的GZ 曲線進(jìn)行擬合求得。將（17）式和（18）式代入（1）式，得到船舶橫搖非線性運(yùn)動(dòng)方程如下

由（10）式和（18）式可得

將（17）式代入（16）式，可得

式中，

P( a )和Q( a )是幅值a的函數(shù)，可以用多項(xiàng)式函數(shù)( a )和( a )近似表示：

式中，系數(shù)p11，p12，q11和q12可通過最小二乘法進(jìn)行數(shù)據(jù)擬合得到。于是（21）式可寫為

橫搖固有周期T0等于衰減數(shù)據(jù)的平均周期[13]，可由（26）式確定

式中，x0是初始橫搖角。根據(jù)橫搖衰減數(shù)據(jù)，采用間隔半個(gè)周期的橫搖角幅值差與相鄰兩次的平均橫搖角幅值繪制橫搖消滅曲線，對(duì)數(shù)據(jù)進(jìn)行擬合可得到曲線近似表達(dá)式為

由（25）和（27）式可求得阻尼系數(shù)c1和c2的估算值：

（28）式中阻尼系數(shù)的識(shí)別精度與初始橫搖角有關(guān)。為提高阻尼系數(shù)識(shí)別精度，本文進(jìn)行如下的迭代誤差修正：將（28）式中的阻尼系數(shù)初次估算值和代入船舶橫搖運(yùn)動(dòng)方程（19）式中，初始條件不變，仍為x( 0 )= x0和x(0)= 0，采用龍格庫(kù)塔方法進(jìn)行數(shù)值求解，可仿真生成自由橫搖衰減曲線，求得新的橫搖消滅曲線

同理，再次根據(jù)（28）式，可以求得阻尼系數(shù)的第二次估算值：

前兩次阻尼系數(shù)估算值的相對(duì)誤差為：

3 有效性驗(yàn)證

為驗(yàn)證本文提出的阻尼系數(shù)識(shí)別方法的有效性，選用（33）式的橫搖動(dòng)力學(xué)非線性方程作為算例[14]，初始條件x( 0 )=0.3 rad( 0 )= 0。通過數(shù)值方法模擬生成自由橫搖衰減曲線，然后用本文方法識(shí)別其阻尼系數(shù)，并與精確值相比較以驗(yàn)證方法的有效性。橫搖運(yùn)動(dòng)響應(yīng)采用四階龍格庫(kù)塔法求解，時(shí)間間隔為0.01s，計(jì)算時(shí)長(zhǎng)為50s。圖1 和圖2 是模擬得到的橫搖衰減曲線及相應(yīng)的消滅曲線。回復(fù)力矩表達(dá)式為g( x )= 2.25x - x3，具有強(qiáng)非線性特征，其最大回復(fù)力矩對(duì)應(yīng)的橫搖角為0.866 rad，回復(fù)力矩曲線如圖3所示。

圖1 模擬自由衰減曲線 Fig.1 Simulated free-decay curve

圖2 橫搖消滅曲線Fig.2 Roll extinction curve

根據(jù)消滅曲線數(shù)據(jù)擬合結(jié)果，（27）式可寫為

根據(jù)（20）式和（11）式計(jì)算可得

根據(jù)（23）式可得

圖3 非線性回復(fù)力矩Fig.3 Nonlinear restoring moment

由（26）式可得橫搖衰減運(yùn)動(dòng)固有周期T0=4.253 16 s，根據(jù)（28）式計(jì)算可得到阻尼系數(shù)的第一次估算值為= 0.147 83 和= 0.200 70；將和作為新的阻尼系數(shù)代入到（33）式中，重復(fù)上述計(jì)算過程，得到阻尼系數(shù)第二次估算值=0.14567和=0.20144，前兩次阻尼系數(shù)估算值之間的相對(duì)誤差為η1=-1.458 74%和η2=0.367 16%；阻尼系數(shù)的第三次估算值=0.143 53和= 0.202 21，后兩次阻尼系數(shù)估算值之間的相對(duì)誤差為=-1.469 20%和= 0.383 33%。于是，根據(jù)（32）式可求得一階和二階阻尼系數(shù)的最終估算值c1 =0.150 00和c2= 0.200 00，與（33）式中阻尼系數(shù)準(zhǔn)確值相等，識(shí)別誤差為零。

為考察本文方法的精度，開展不同橫搖角幅值下的參數(shù)研究，將船舶橫搖方程的阻尼系數(shù)識(shí)別結(jié)果與采用其他方法得到的估算結(jié)果進(jìn)行對(duì)比分析。

采用本文方法對(duì)（33）式在不同初始橫搖角時(shí)的非線性阻尼系數(shù)進(jìn)行了識(shí)別，計(jì)算結(jié)果見表1，同時(shí)將本文中阻尼系數(shù)的識(shí)別結(jié)果和相對(duì)誤差等與文獻(xiàn)[14]的計(jì)算結(jié)果進(jìn)行了比較，如表2所示。由表1 可知：隨著初始橫搖角的增大，阻尼系數(shù)識(shí)別質(zhì)量逐漸變差，相比于一階阻尼項(xiàng)系數(shù)而言，二階阻尼項(xiàng)系數(shù)識(shí)別受初始橫搖角的影響更為敏感，表現(xiàn)在其相對(duì)誤差出現(xiàn)時(shí)對(duì)應(yīng)的橫搖角更小，且相對(duì)誤差值更大。同時(shí)由表2 的對(duì)比可以看出：無論是對(duì)小初始橫搖角還是大初始橫搖角，文獻(xiàn)[14]的阻尼系數(shù)識(shí)別值均存在相對(duì)誤差。而本文方法的識(shí)別結(jié)果在初始橫搖角不大于0.6 rad 范圍內(nèi)，均能夠較為精確地識(shí)別出一階和二階阻尼系數(shù)，識(shí)別誤差不大于0.03%；在大初始橫搖角度時(shí)，本文阻尼系數(shù)識(shí)別精度明顯高于文獻(xiàn)[14]，如在最大橫搖角0.866 rad 的初始條件下，本文的一階和二階阻尼系數(shù)識(shí)別誤差分別為0.26%和0.415%，而文獻(xiàn)[14]給出的阻尼系數(shù)識(shí)別誤差分別為1.22%和1%，本文方法的一階和二階阻尼系數(shù)識(shí)別誤差僅為文獻(xiàn)[14]的0.21倍和0.42倍。

表1 不同初始橫搖角下（33）式的阻尼系數(shù)識(shí)別結(jié)果Tab.1 Estimation of damping coefficients from Eq.(33)under different initial roll angles

同時(shí)，選取如下船舶橫搖方程進(jìn)行研究[5]

表2 本文識(shí)別的阻尼系數(shù)與文獻(xiàn)[14]結(jié)果比較Tab.2 Comparison of the present work and the results from Ref.[14]

4 實(shí)例分析

將本文阻尼系數(shù)識(shí)別方法應(yīng)用到實(shí)際船舶橫搖衰減數(shù)據(jù)中。船模試驗(yàn)得到的橫搖角峰值如表4所示[5]，船模的自搖橫搖周期為2 s，對(duì)應(yīng)實(shí)船的船模回復(fù)力矩為g( x )= 9.87x + 0.001 974x3。阻尼模型選用線性加二次阻尼模型，于是船舶模型的橫搖運(yùn)動(dòng)方程如下：

根據(jù)回復(fù)力矩表達(dá)式，可求得Φ0( a,φ )的表達(dá)式

計(jì)算得到（23）式中擬合系數(shù)p11=0.500 0,p12=0,q11=0,q12=1.333 4；由橫搖衰減峰值數(shù)據(jù)可得到船模的橫搖消滅曲線，進(jìn)行數(shù)據(jù)擬合計(jì)算得到（27）式中擬合系數(shù)λ1= 0.678 4，λ2= 0.073 9。

圖4 模擬的自由橫搖衰減曲線和試驗(yàn)峰值數(shù)據(jù)比較Fig.4 Comparison of experimental data on free-roll decay motion with the numerically simulated results

根據(jù)（28）式，可得阻尼系數(shù)c1和c2的第一次估算值= 0.147 9 和= 0.508 8。重復(fù)前文所述計(jì)算流程，得到阻尼系數(shù)最終估算值c1= 0.138 7，c2= 0.528 2。將該阻尼系數(shù)代入到（39）式中，通過數(shù)值模擬得到自由衰減曲線，并與衰減試驗(yàn)峰值數(shù)據(jù)進(jìn)行比較，對(duì)比結(jié)果見圖4。由圖4 可知，本文估算得到的橫搖衰減曲線峰值與船模試驗(yàn)得到的橫搖峰值吻合較好，進(jìn)一步驗(yàn)證了本文方法的有效性。

5 結(jié) 論

本文基于漸進(jìn)理論，發(fā)展了一種利用試驗(yàn)衰減曲線進(jìn)行非線性橫搖阻尼系數(shù)識(shí)別的方法。通過數(shù)值仿真模擬及與已有研究成果的對(duì)比分析研究，驗(yàn)證了本文提出的阻尼系數(shù)識(shí)別方法的有效性和高精度，最后將該方法應(yīng)用于實(shí)際船舶算例，模擬計(jì)算得到的自由衰減結(jié)果與試驗(yàn)數(shù)據(jù)吻合較好。通過本文的研究，可以得到以下結(jié)論：

（1）僅需要自由橫搖衰減曲線峰值數(shù)據(jù)，利用本文提出的阻尼系數(shù)識(shí)別方法可精確估算出非線性橫搖阻尼系數(shù)；

（2）本文提出的非線性阻尼系數(shù)識(shí)別方法適用于強(qiáng)非線性的回復(fù)力矩情形，回復(fù)力矩可以是任意形式奇次多項(xiàng)式函數(shù)；

（3）本文提出的非線性阻尼系數(shù)識(shí)別方法適用于大角度橫搖，且具有較高的識(shí)別精度。雖然隨著初始橫搖角幅值的增大，阻尼系數(shù)識(shí)別質(zhì)量下降，二階阻尼系數(shù)的識(shí)別誤差相對(duì)更高、出現(xiàn)誤差時(shí)對(duì)應(yīng)的初始橫搖角度也更小。但是在最大回復(fù)力矩對(duì)應(yīng)的橫搖角范圍內(nèi)，一階和二階阻尼系數(shù)識(shí)別的相對(duì)誤差仍非常小，在可接受范圍內(nèi)，特別是在小角度初始橫搖角幅值條件下，阻尼系數(shù)的識(shí)別誤差接近為零。