周小宇，李紅霞，黃 一

(大連理工大學(xué) 船舶工程學(xué)院，遼寧 大連 116024)

一直以來(lái)，船舶橫搖是船舶穩(wěn)性研究的關(guān)鍵問(wèn)題.近些年發(fā)生的一些航運(yùn)事故(如1993年，C11集裝箱船在北大西洋航行發(fā)生大幅度橫搖運(yùn)動(dòng))表明：船舶在滿足靜穩(wěn)性要求的情況下仍然可能發(fā)生大幅度橫搖.經(jīng)大量學(xué)者研究發(fā)現(xiàn)，參數(shù)橫搖是導(dǎo)致這種大幅橫搖的重要原因之一[1-5].

目前針對(duì)船舶參數(shù)橫搖的研究主要是對(duì)參數(shù)橫搖響應(yīng)進(jìn)行模擬和預(yù)測(cè).然而，工程上更加關(guān)心橫搖的極值響應(yīng).一方面，橫搖極值響應(yīng)對(duì)于計(jì)算船舶的傾覆概率至關(guān)重要；另一方面，橫搖極值響應(yīng)可以作為輸入來(lái)評(píng)估船體結(jié)構(gòu)的安全性和可靠性.而目前針對(duì)船舶參數(shù)橫搖的極值響應(yīng)問(wèn)題的研究較少.

本文提出了一種基于窄帶隨機(jī)過(guò)程理論和Hermite變換法的極值分析方法，并且利用該方法對(duì)C11集裝箱船的參數(shù)橫搖極值響應(yīng)進(jìn)行了研究.首先將船舶參數(shù)橫搖假設(shè)為窄帶隨機(jī)過(guò)程，再通過(guò)Hilbert變換得到其分量表達(dá)形式[6].根據(jù)窄帶隨機(jī)過(guò)程理論，窄帶隨機(jī)過(guò)程的分量過(guò)程和原過(guò)程的極大值具有緊密的關(guān)系，比原過(guò)程更加能夠反映極值的概率特性.因此，接下來(lái)可以利用Hermite變換得到分量過(guò)程和某個(gè)慢變高斯隨機(jī)過(guò)程之間的非線性關(guān)系.本文假設(shè)隨機(jī)參數(shù)橫搖過(guò)程也可近似利用這個(gè)非線性關(guān)系映射到一個(gè)底層高斯隨機(jī)過(guò)程.因此，當(dāng)根據(jù)高斯隨機(jī)過(guò)程的極值理論計(jì)算出底層過(guò)程的極值之后，便可利用上述的非線性關(guān)系計(jì)算船舶隨機(jī)參數(shù)橫搖的極值響應(yīng).本文利用窄帶分量過(guò)程是因?yàn)槠渑c極值聯(lián)系更為緊密，建立的非線性關(guān)系更適合極值的預(yù)測(cè).但是在計(jì)算極值時(shí)，由于橫搖過(guò)程并不嚴(yán)格滿足窄帶的條件，所以需考慮帶寬的影響.

1 C11船參數(shù)橫搖運(yùn)動(dòng)模型

根據(jù)Dostal等[7-8]的研究，船舶隨機(jī)橫搖運(yùn)動(dòng)可以采用單自由度方程來(lái)描述：

gΔlGZ(φ,η,φ)=M(t)

(1)

式中：Ixx為船舶橫搖轉(zhuǎn)動(dòng)慣量；δIxx為附加橫搖轉(zhuǎn)動(dòng)慣量；φ為橫搖角；c1、c3分別為一次和三次阻尼系數(shù)；g為重力加速度；η為波高；φ為波峰與船中相對(duì)位置；Δ為排水量；lGZ為復(fù)原力臂；M(t)(t為時(shí)間)為橫搖的強(qiáng)迫激勵(lì)力矩.由于本文主要研究的是C11船在縱浪情況下的參數(shù)橫搖，所以M(t)≡0.根據(jù)文獻(xiàn)[7],lGZ和φ、η及φ有關(guān).文獻(xiàn)[8]表明，lGZ可以近似表示為多項(xiàng)式與傅里葉級(jí)數(shù)混合的形式，因此船舶隨機(jī)橫搖運(yùn)動(dòng)方程可以簡(jiǎn)化為

k1φ+k3φ3+k5φ5+q1φζ(t)=0

(2)

式中：k1、k3、k5及q1為模型系數(shù)，需根據(jù)lGZ的擬合而定；ζ(t)為隨機(jī)波面函數(shù).在擬合lGZ時(shí)，假設(shè)船舶垂蕩和縱搖運(yùn)動(dòng)為準(zhǔn)平衡態(tài)，作用在船體濕表面上的水壓力可表示為

(3)

式中：k為波浪波數(shù)；d為水深；x、z為波浪橫向和縱向坐標(biāo)；ke、ωe為遭遇波數(shù)和遭遇頻率.

可以根據(jù)式(3)對(duì)壓力在瞬時(shí)濕表面積分求出C11集裝箱船在不同φ，η和φ情況下的lGZ，接著根據(jù)式(2)的關(guān)系擬合其中的系數(shù).需要注意的是，為了獲得更準(zhǔn)確的系數(shù)，達(dá)到較好的擬合精度，應(yīng)當(dāng)適當(dāng)增加高次擬合項(xiàng).增加的項(xiàng)對(duì)應(yīng)的都是高頻激勵(lì)項(xiàng)，在動(dòng)力分析中可以忽略.圖1和2給出了不同φ，η及φ情況下橫穩(wěn)性臂的計(jì)算結(jié)果以及擬合結(jié)果.式(2)系數(shù)的擬合計(jì)算結(jié)果k1=0.060 9 s-2,k3=0.043 8 s-2,k5=-0.070 4 s-2,c1=0.008 4 s-1,c3=5.299 s,q1=0.021 3 m-1·s-1.

圖1 固定η下lGZ的計(jì)算和擬合結(jié)果Fig.1 Calculation and fitting results of lGZ at fixed η

圖2 固定φ下lGZ的計(jì)算和擬合結(jié)果Fig.2 Calculation and fitting results of lGZ at fixed φ

2 船舶參數(shù)橫搖極值響應(yīng)分析方法

由上文可知，C11集裝箱船的隨機(jī)參數(shù)橫搖運(yùn)動(dòng)可由式(2)來(lái)刻畫.假設(shè)海浪滿足PM譜：

(4)

式中:ζ為隨機(jī)波面函數(shù)；ω為圓頻率；Hs為義波高；Tp和ωp分別為譜峰周期和譜峰頻率.為了對(duì)比已有的試驗(yàn)結(jié)果，本文將有義波高10.43 m，譜峰周期9.99 s設(shè)定為所研究的海況，并且假設(shè)船舶迎浪，船速為0，在此海況下停留100 min.利用Monte Carlo模擬方法，可以得到C11集裝箱船的隨機(jī)參數(shù)橫搖的時(shí)間歷程樣本.所采用的Monte Carlo模擬方法的基本步驟如下：① 首先利用偽隨機(jī)波法生成一系列功率譜滿足所研究海況的波浪時(shí)間歷程；② 將這些波浪時(shí)程代入式(2)替換其中的ζ(t)；③ 利用4階Runge-Kutta法計(jì)算出在不同波浪時(shí)程下的響應(yīng)時(shí)程，這些響應(yīng)時(shí)程便是Monte Carlo模擬結(jié)果.忽略初始很短時(shí)間的不穩(wěn)定狀態(tài)后，C11船隨機(jī)參數(shù)橫搖是一個(gè)穩(wěn)態(tài)的隨機(jī)過(guò)程，可以按照定義計(jì)算出它的自相關(guān)函數(shù).利用譜分析理論，可以計(jì)算出C11船隨機(jī)參數(shù)橫搖運(yùn)動(dòng)的功率譜密度.圖3所示為C11船在所研究海況下橫搖的時(shí)間歷程曲線，圖4所示為其功率譜密度(PSD)，圖中f′為頻率.譜寬系數(shù)ε可以按下式計(jì)算：

圖3 研究海況下C11船橫搖運(yùn)動(dòng)時(shí)間歷程Fig.3 Time history of C11 roll in investigated sea state conditions

(5)

式中：m0、m2及m4分別為0階、2階及4階譜矩.計(jì)算得到其譜寬系數(shù)為0.61，說(shuō)明橫搖過(guò)程是個(gè)介于窄帶和寬帶之間的隨機(jī)過(guò)程.根據(jù)上述Monte Carlo模擬方法，生成104條時(shí)間歷程樣本，得到C11船橫搖極值響應(yīng)的均值為35.7°.文獻(xiàn)[9]提供了C11船在本文選用的海況和航行狀態(tài)下10組隨機(jī)參數(shù)橫搖模型試驗(yàn)的橫搖角極值(見(jiàn)表 1).該試驗(yàn)于2010年在日本大阪大學(xué)拖曳水池中完成，模型縮尺比為1∶100，每組試驗(yàn)時(shí)長(zhǎng)約為10 min.這10個(gè)極值的均值為36.7°，與本文Monte Carlo結(jié)果相差2.72%，驗(yàn)證了本文計(jì)算模型的準(zhǔn)確性.

圖4 研究海況下C11船橫搖運(yùn)動(dòng)功率譜密度Fig.4 PSD of C11 roll in investigated sea state conditions

表1 研究海況下C11船橫搖模型試驗(yàn)結(jié)果Tab.1 Model test results of C11 roll in investigated sea state conditions

觀察C11船隨機(jī)參數(shù)橫搖運(yùn)動(dòng)的功率譜密度可以發(fā)現(xiàn)，該運(yùn)動(dòng)非常接近窄帶隨機(jī)過(guò)程.因此，利用窄帶隨機(jī)過(guò)程的相關(guān)理論對(duì)其進(jìn)行研究.根據(jù)窄帶隨機(jī)過(guò)程的理論，C11船隨機(jī)橫搖過(guò)程可以寫成如下的分量形式：

φ(t)=φc(t)cosω0t-φs(t)sinω0t

(6)

式中：φc(t)為同相分量；φs(t)為正交分量；ω0為中心頻率.根據(jù)窄帶隨機(jī)過(guò)程理論，φc(t)和φ(t)的1階矩和2階矩是相同的.本文假設(shè)如果一個(gè)映射f能將φc(t)映射為一個(gè)慢變標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程Uc，那么這個(gè)映射f也能將φ(t)映射為一個(gè)標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程U(t).

在C11船隨機(jī)橫搖過(guò)程φ(t)已知的情況下，兩個(gè)分量可以通過(guò)下式計(jì)算得到：

(7)

(8)

圖5 C11船橫搖運(yùn)動(dòng)同相分量Fig.5 In-phase component of C11 roll

許多學(xué)者已證明參數(shù)橫搖響應(yīng)是非高斯的，因此本文在研究C11船橫搖極值響應(yīng)時(shí)利用了Hermite變換法.Hermite變換法是由Winterstein首先提出的一種將非高斯隨機(jī)變量轉(zhuǎn)化為高斯隨機(jī)變量的多項(xiàng)式形式的變換方法[10].該變換方法最終可以將一個(gè)非高斯隨機(jī)變量轉(zhuǎn)化為標(biāo)準(zhǔn)高斯隨機(jī)變量的Hermite多項(xiàng)式級(jí)數(shù)的形式.隨后，一些學(xué)者在Winterstein研究基礎(chǔ)上對(duì)Hermite變換法進(jìn)行了更深入研究[11].但是目前，最為成熟的Hermite變換法仍是3階Hermite變換方法[12].按照上文所述，同相分量φc(或正交分量φs)是將要進(jìn)行Hermite變換的非高斯隨機(jī)過(guò)程.通過(guò)試算可知，φc(或φs)的峰度小于3，因此采用硬化過(guò)程的3階Hermite變換法.將φc進(jìn)行歸一化處理得到Z=(φc-μφc)/σφc(μφc為φc的均值，σφc為φc的標(biāo)準(zhǔn)差).這樣，硬化過(guò)程3階Hermite變換表達(dá)式可表示為[13]

Uc=Z-c3(Z2-1)-c4(Z3-3Z)

(9)

式中：Uc為慢變的標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程；c3、c4為模型系數(shù).文獻(xiàn)[13]提到，這2個(gè)系數(shù)與Z的3階中心矩α3Z和4階中心矩α4Z近似滿足下述關(guān)系：

(10)

c3和c4確定之后，可以反向求解出Z和Uc之間的關(guān)系，設(shè)為Z=g(Uc).那么φc和Uc便滿足如下關(guān)系：

φc=μφc+σφcg(Uc)=f(Uc)

(11)

同樣假設(shè)原過(guò)程φ和某個(gè)標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程U(t)也滿足這種關(guān)系，即φ=f(U(t)).這里的標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程U(t)不再是慢變的，而是和原過(guò)程一樣能量處于中心頻率附近.

本文提出的分析船舶參數(shù)橫搖極值響應(yīng)的方法就是基于上述動(dòng)力學(xué)模型，窄帶隨機(jī)過(guò)程理論以及Hermite變換法.其基本的計(jì)算流程如下.

(1)按照式(3)計(jì)算船舶在不同φ，不同η和不同φ時(shí)的lGZ.

(2)擬合計(jì)算得到的lGZ,得到式(2)中的系數(shù)，即得到了描述船舶隨機(jī)橫搖運(yùn)動(dòng)的簡(jiǎn)化單自由度動(dòng)力學(xué)模型.

(3)根據(jù)式(2)，利用Monte Carlo方法計(jì)算船舶隨機(jī)橫搖的幾條時(shí)間歷程樣本.

(4)根據(jù)窄帶隨機(jī)過(guò)程理論，計(jì)算出橫搖響應(yīng)過(guò)程的φc(或φs)的樣本數(shù)據(jù).

(5)根據(jù)式(9)和(11)，利用硬化過(guò)程3階Hermite變換法得到φc和Uc之間的非線性關(guān)系為：φc=f(Uc).

(6)按照上文的假設(shè)，原過(guò)程φ(t)和某個(gè)標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程U(t)也近似滿足φ=f(U(t)).假設(shè)Ue表示U(t)的極值，那么對(duì)于C11船參數(shù)橫搖極值響應(yīng)的均值，可近似利用Hermite變換得到，如下式：

E(φe)=f[E(Ue)]

(12)

式中：E為期望.由于橫搖過(guò)程不是嚴(yán)格的窄帶過(guò)程，所以與之對(duì)應(yīng)的標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程U(t)也不是嚴(yán)格窄帶的.在計(jì)算U(t)的極值時(shí)，需要用到零均值平穩(wěn)高斯隨機(jī)過(guò)程極大值的精確概率密度Pp(x)[14]：

(13)

圖6 標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程U(t)極大值的精確概率密度Fig.6 Exact probability density of maximum value of U(t)

式中：x為概率密度函數(shù)的自變量；σ為高斯隨機(jī)過(guò)程的標(biāo)準(zhǔn)差.本文假設(shè)非線性變換對(duì)功率譜的改變不大，因此變換后得到的標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程U(t)的ε也為0.61.圖6給出了式(13)在σ=1、ε=0.61時(shí)極大值的精確概率密度.從圖中可以看出，由于帶寬的影響，極大值分布已經(jīng)不再是瑞利分布.

在假設(shè)極大值之間相互獨(dú)立的條件下，極值分布Fm(x)可以按照下式計(jì)算：

(14)

式中：Fp(x)為極大值分布；n為所研究時(shí)間內(nèi)極大值的個(gè)數(shù).

3 C11船參數(shù)橫搖極值響應(yīng)分析

通過(guò)表2還可以發(fā)現(xiàn)，利用傳統(tǒng)Gumbel極值模型[15]對(duì)C11船參數(shù)橫搖極值響應(yīng)的預(yù)測(cè)是偏大的，若以數(shù)值計(jì)算結(jié)果為標(biāo)準(zhǔn)，傳統(tǒng)Gumbel方法的預(yù)測(cè)誤差為21%，造成誤差的原因主要是傳統(tǒng)Gumbel方法將響應(yīng)假設(shè)為高斯隨機(jī)過(guò)程.根據(jù)本文的計(jì)算結(jié)果可以再次證實(shí)C11船參數(shù)橫搖響應(yīng)是非高斯隨機(jī)過(guò)程.同時(shí)，本文的計(jì)算結(jié)果還可證明傳統(tǒng)Gumbel方法對(duì)于參數(shù)橫搖極值的預(yù)測(cè)是不適用的.

表2 C11船橫搖極值響應(yīng)均值預(yù)測(cè)結(jié)果Tab.2 Estimations of mean value of C11 roll extreme response

通過(guò)表2還可以看出，前3種極值預(yù)測(cè)方法的結(jié)果都與文獻(xiàn)[9]中的模型試驗(yàn)結(jié)果有一定偏差.若以模型試驗(yàn)結(jié)果作為標(biāo)準(zhǔn)，本文所提出的方法和數(shù)值計(jì)算的誤差約為10%，傳統(tǒng)Gumbel方法的預(yù)測(cè)誤差約為30%.由于式(14)的計(jì)算方法和數(shù)值計(jì)算方法均考慮了帶寬的影響，但未考慮極大值之間的相關(guān)性，所以推測(cè)造成這一誤差的原因主要是忽略了極大值之間的相關(guān)性.Naess[16]在其研究中討論了在高斯隨機(jī)過(guò)程中考慮相鄰極大值相關(guān)性對(duì)于預(yù)測(cè)結(jié)果的影響.他的結(jié)論是，考慮相關(guān)性后，極值的預(yù)測(cè)結(jié)果會(huì)變小.這個(gè)結(jié)論對(duì)于非高斯隨機(jī)過(guò)程也是適用的.但應(yīng)該看到，忽略極大值之間相關(guān)性的極值模型預(yù)測(cè)結(jié)果偏于保守，符合工程的安全性原則.因此，對(duì)于工程應(yīng)用而言，本文提出的方法具有參考價(jià)值.

4 結(jié)論

本文基于窄帶隨機(jī)過(guò)程理論和Hermite模型，提出了一種預(yù)測(cè)船舶橫搖極值響應(yīng)的方法.采用窄帶隨機(jī)過(guò)程理論主要是為了求出響應(yīng)的窄帶分量，其原因是窄帶分量能夠更好地體現(xiàn)隨機(jī)響應(yīng)極值的概率特性.接下來(lái)采用Hermite模型主要是為了解決響應(yīng)過(guò)程的非高斯性問(wèn)題.以C11集裝箱船為例，預(yù)測(cè)出某一特定海況下隨機(jī)參數(shù)橫搖極值響應(yīng)的均值并與模型試驗(yàn)結(jié)果進(jìn)行了比較，本文提出的方法在極值預(yù)測(cè)上能夠給出偏于安全的預(yù)測(cè)值，可為實(shí)際工程提供參考.主要結(jié)論如下：

(1)傳統(tǒng)的Gumbel極值預(yù)測(cè)模型由于假設(shè)響應(yīng)為高斯隨機(jī)過(guò)程，對(duì)C11船橫搖極值響應(yīng)的預(yù)測(cè)誤差在18%左右,不適合預(yù)測(cè)參數(shù)橫搖響應(yīng)極值.

(2)Monte Carlo預(yù)測(cè)的極值響應(yīng)均值與模型試驗(yàn)誤差很小，驗(yàn)證了本文數(shù)學(xué)模型的準(zhǔn)確性.

(3)本文預(yù)測(cè)方法在預(yù)測(cè)精度上相當(dāng)于基于104樣本的數(shù)值計(jì)算，但是在計(jì)算成本上只有數(shù)值計(jì)算的1/500，驗(yàn)證了本文方法的高效性.

(4)本文方法預(yù)測(cè)的極值響應(yīng)比模型試驗(yàn)結(jié)果大10%.推測(cè)其原因是忽略了響應(yīng)極大值之間的相關(guān)性，這將會(huì)是作者后續(xù)研究的重點(diǎn).