周小宇,李紅霞,黃 一
(大連理工大學 船舶工程學院,遼寧 大連 116024)
一直以來,船舶橫搖是船舶穩(wěn)性研究的關鍵問題.近些年發(fā)生的一些航運事故(如1993年,C11集裝箱船在北大西洋航行發(fā)生大幅度橫搖運動)表明:船舶在滿足靜穩(wěn)性要求的情況下仍然可能發(fā)生大幅度橫搖.經(jīng)大量學者研究發(fā)現(xiàn),參數(shù)橫搖是導致這種大幅橫搖的重要原因之一[1-5].
目前針對船舶參數(shù)橫搖的研究主要是對參數(shù)橫搖響應進行模擬和預測.然而,工程上更加關心橫搖的極值響應.一方面,橫搖極值響應對于計算船舶的傾覆概率至關重要;另一方面,橫搖極值響應可以作為輸入來評估船體結構的安全性和可靠性.而目前針對船舶參數(shù)橫搖的極值響應問題的研究較少.
本文提出了一種基于窄帶隨機過程理論和Hermite變換法的極值分析方法,并且利用該方法對C11集裝箱船的參數(shù)橫搖極值響應進行了研究.首先將船舶參數(shù)橫搖假設為窄帶隨機過程,再通過Hilbert變換得到其分量表達形式[6].根據(jù)窄帶隨機過程理論,窄帶隨機過程的分量過程和原過程的極大值具有緊密的關系,比原過程更加能夠反映極值的概率特性.因此,接下來可以利用Hermite變換得到分量過程和某個慢變高斯隨機過程之間的非線性關系.本文假設隨機參數(shù)橫搖過程也可近似利用這個非線性關系映射到一個底層高斯隨機過程.因此,當根據(jù)高斯隨機過程的極值理論計算出底層過程的極值之后,便可利用上述的非線性關系計算船舶隨機參數(shù)橫搖的極值響應.本文利用窄帶分量過程是因為其與極值聯(lián)系更為緊密,建立的非線性關系更適合極值的預測.但是在計算極值時,由于橫搖過程并不嚴格滿足窄帶的條件,所以需考慮帶寬的影響.
根據(jù)Dostal等[7-8]的研究,船舶隨機橫搖運動可以采用單自由度方程來描述:
gΔlGZ(φ,η,φ)=M(t)
(1)
式中:Ixx為船舶橫搖轉動慣量;δIxx為附加橫搖轉動慣量;φ為橫搖角;c1、c3分別為一次和三次阻尼系數(shù);g為重力加速度;η為波高;φ為波峰與船中相對位置;Δ為排水量;lGZ為復原力臂;M(t)(t為時間)為橫搖的強迫激勵力矩.由于本文主要研究的是C11船在縱浪情況下的參數(shù)橫搖,所以M(t)≡0.根據(jù)文獻[7],lGZ和φ、η及φ有關.文獻[8]表明,lGZ可以近似表示為多項式與傅里葉級數(shù)混合的形式,因此船舶隨機橫搖運動方程可以簡化為
k1φ+k3φ3+k5φ5+q1φζ(t)=0
(2)
式中:k1、k3、k5及q1為模型系數(shù),需根據(jù)lGZ的擬合而定;ζ(t)為隨機波面函數(shù).在擬合lGZ時,假設船舶垂蕩和縱搖運動為準平衡態(tài),作用在船體濕表面上的水壓力可表示為
(3)
式中:k為波浪波數(shù);d為水深;x、z為波浪橫向和縱向坐標;ke、ωe為遭遇波數(shù)和遭遇頻率.
可以根據(jù)式(3)對壓力在瞬時濕表面積分求出C11集裝箱船在不同φ,η和φ情況下的lGZ,接著根據(jù)式(2)的關系擬合其中的系數(shù).需要注意的是,為了獲得更準確的系數(shù),達到較好的擬合精度,應當適當增加高次擬合項.增加的項對應的都是高頻激勵項,在動力分析中可以忽略.圖1和2給出了不同φ,η及φ情況下橫穩(wěn)性臂的計算結果以及擬合結果.式(2)系數(shù)的擬合計算結果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的計算和擬合結果Fig.1 Calculation and fitting results of lGZ at fixed η
圖2 固定φ下lGZ的計算和擬合結果Fig.2 Calculation and fitting results of lGZ at fixed φ
由上文可知,C11集裝箱船的隨機參數(shù)橫搖運動可由式(2)來刻畫.假設海浪滿足PM譜:
(4)
式中:ζ為隨機波面函數(shù);ω為圓頻率;Hs為義波高;Tp和ωp分別為譜峰周期和譜峰頻率.為了對比已有的試驗結果,本文將有義波高10.43 m,譜峰周期9.99 s設定為所研究的海況,并且假設船舶迎浪,船速為0,在此海況下停留100 min.利用Monte Carlo模擬方法,可以得到C11集裝箱船的隨機參數(shù)橫搖的時間歷程樣本.所采用的Monte Carlo模擬方法的基本步驟如下:① 首先利用偽隨機波法生成一系列功率譜滿足所研究海況的波浪時間歷程;② 將這些波浪時程代入式(2)替換其中的ζ(t);③ 利用4階Runge-Kutta法計算出在不同波浪時程下的響應時程,這些響應時程便是Monte Carlo模擬結果.忽略初始很短時間的不穩(wěn)定狀態(tài)后,C11船隨機參數(shù)橫搖是一個穩(wěn)態(tài)的隨機過程,可以按照定義計算出它的自相關函數(shù).利用譜分析理論,可以計算出C11船隨機參數(shù)橫搖運動的功率譜密度.圖3所示為C11船在所研究海況下橫搖的時間歷程曲線,圖4所示為其功率譜密度(PSD),圖中f′為頻率.譜寬系數(shù)ε可以按下式計算:
圖3 研究海況下C11船橫搖運動時間歷程Fig.3 Time history of C11 roll in investigated sea state conditions
(5)
式中:m0、m2及m4分別為0階、2階及4階譜矩.計算得到其譜寬系數(shù)為0.61,說明橫搖過程是個介于窄帶和寬帶之間的隨機過程.根據(jù)上述Monte Carlo模擬方法,生成104條時間歷程樣本,得到C11船橫搖極值響應的均值為35.7°.文獻[9]提供了C11船在本文選用的海況和航行狀態(tài)下10組隨機參數(shù)橫搖模型試驗的橫搖角極值(見表 1).該試驗于2010年在日本大阪大學拖曳水池中完成,模型縮尺比為1∶100,每組試驗時長約為10 min.這10個極值的均值為36.7°,與本文Monte Carlo結果相差2.72%,驗證了本文計算模型的準確性.
圖4 研究海況下C11船橫搖運動功率譜密度Fig.4 PSD of C11 roll in investigated sea state conditions
表1 研究海況下C11船橫搖模型試驗結果Tab.1 Model test results of C11 roll in investigated sea state conditions
觀察C11船隨機參數(shù)橫搖運動的功率譜密度可以發(fā)現(xiàn),該運動非常接近窄帶隨機過程.因此,利用窄帶隨機過程的相關理論對其進行研究.根據(jù)窄帶隨機過程的理論,C11船隨機橫搖過程可以寫成如下的分量形式:
φ(t)=φc(t)cosω0t-φs(t)sinω0t
(6)
式中:φc(t)為同相分量;φs(t)為正交分量;ω0為中心頻率.根據(jù)窄帶隨機過程理論,φc(t)和φ(t)的1階矩和2階矩是相同的.本文假設如果一個映射f能將φc(t)映射為一個慢變標準高斯隨機過程Uc,那么這個映射f也能將φ(t)映射為一個標準高斯隨機過程U(t).
在C11船隨機橫搖過程φ(t)已知的情況下,兩個分量可以通過下式計算得到:
(7)
(8)
圖5 C11船橫搖運動同相分量Fig.5 In-phase component of C11 roll
許多學者已證明參數(shù)橫搖響應是非高斯的,因此本文在研究C11船橫搖極值響應時利用了Hermite變換法.Hermite變換法是由Winterstein首先提出的一種將非高斯隨機變量轉化為高斯隨機變量的多項式形式的變換方法[10].該變換方法最終可以將一個非高斯隨機變量轉化為標準高斯隨機變量的Hermite多項式級數(shù)的形式.隨后,一些學者在Winterstein研究基礎上對Hermite變換法進行了更深入研究[11].但是目前,最為成熟的Hermite變換法仍是3階Hermite變換方法[12].按照上文所述,同相分量φc(或正交分量φs)是將要進行Hermite變換的非高斯隨機過程.通過試算可知,φc(或φs)的峰度小于3,因此采用硬化過程的3階Hermite變換法.將φc進行歸一化處理得到Z=(φc-μφc)/σφc(μφc為φc的均值,σφc為φc的標準差).這樣,硬化過程3階Hermite變換表達式可表示為[13]
Uc=Z-c3(Z2-1)-c4(Z3-3Z)
(9)
式中:Uc為慢變的標準高斯隨機過程;c3、c4為模型系數(shù).文獻[13]提到,這2個系數(shù)與Z的3階中心矩α3Z和4階中心矩α4Z近似滿足下述關系:
(10)
c3和c4確定之后,可以反向求解出Z和Uc之間的關系,設為Z=g(Uc).那么φc和Uc便滿足如下關系:
φc=μφc+σφcg(Uc)=f(Uc)
(11)
同樣假設原過程φ和某個標準高斯隨機過程U(t)也滿足這種關系,即φ=f(U(t)).這里的標準高斯隨機過程U(t)不再是慢變的,而是和原過程一樣能量處于中心頻率附近.
本文提出的分析船舶參數(shù)橫搖極值響應的方法就是基于上述動力學模型,窄帶隨機過程理論以及Hermite變換法.其基本的計算流程如下.
(1)按照式(3)計算船舶在不同φ,不同η和不同φ時的lGZ.
(2)擬合計算得到的lGZ,得到式(2)中的系數(shù),即得到了描述船舶隨機橫搖運動的簡化單自由度動力學模型.
(3)根據(jù)式(2),利用Monte Carlo方法計算船舶隨機橫搖的幾條時間歷程樣本.
(4)根據(jù)窄帶隨機過程理論,計算出橫搖響應過程的φc(或φs)的樣本數(shù)據(jù).
(5)根據(jù)式(9)和(11),利用硬化過程3階Hermite變換法得到φc和Uc之間的非線性關系為:φc=f(Uc).
(6)按照上文的假設,原過程φ(t)和某個標準高斯隨機過程U(t)也近似滿足φ=f(U(t)).假設Ue表示U(t)的極值,那么對于C11船參數(shù)橫搖極值響應的均值,可近似利用Hermite變換得到,如下式:
E(φe)=f[E(Ue)]
(12)
式中:E為期望.由于橫搖過程不是嚴格的窄帶過程,所以與之對應的標準高斯隨機過程U(t)也不是嚴格窄帶的.在計算U(t)的極值時,需要用到零均值平穩(wěn)高斯隨機過程極大值的精確概率密度Pp(x)[14]:
(13)
圖6 標準高斯隨機過程U(t)極大值的精確概率密度Fig.6 Exact probability density of maximum value of U(t)
式中:x為概率密度函數(shù)的自變量;σ為高斯隨機過程的標準差.本文假設非線性變換對功率譜的改變不大,因此變換后得到的標準高斯隨機過程U(t)的ε也為0.61.圖6給出了式(13)在σ=1、ε=0.61時極大值的精確概率密度.從圖中可以看出,由于帶寬的影響,極大值分布已經(jīng)不再是瑞利分布.
在假設極大值之間相互獨立的條件下,極值分布Fm(x)可以按照下式計算:
(14)
式中:Fp(x)為極大值分布;n為所研究時間內(nèi)極大值的個數(shù).
通過表2還可以發(fā)現(xiàn),利用傳統(tǒng)Gumbel極值模型[15]對C11船參數(shù)橫搖極值響應的預測是偏大的,若以數(shù)值計算結果為標準,傳統(tǒng)Gumbel方法的預測誤差為21%,造成誤差的原因主要是傳統(tǒng)Gumbel方法將響應假設為高斯隨機過程.根據(jù)本文的計算結果可以再次證實C11船參數(shù)橫搖響應是非高斯隨機過程.同時,本文的計算結果還可證明傳統(tǒng)Gumbel方法對于參數(shù)橫搖極值的預測是不適用的.
表2 C11船橫搖極值響應均值預測結果Tab.2 Estimations of mean value of C11 roll extreme response
通過表2還可以看出,前3種極值預測方法的結果都與文獻[9]中的模型試驗結果有一定偏差.若以模型試驗結果作為標準,本文所提出的方法和數(shù)值計算的誤差約為10%,傳統(tǒng)Gumbel方法的預測誤差約為30%.由于式(14)的計算方法和數(shù)值計算方法均考慮了帶寬的影響,但未考慮極大值之間的相關性,所以推測造成這一誤差的原因主要是忽略了極大值之間的相關性.Naess[16]在其研究中討論了在高斯隨機過程中考慮相鄰極大值相關性對于預測結果的影響.他的結論是,考慮相關性后,極值的預測結果會變小.這個結論對于非高斯隨機過程也是適用的.但應該看到,忽略極大值之間相關性的極值模型預測結果偏于保守,符合工程的安全性原則.因此,對于工程應用而言,本文提出的方法具有參考價值.
本文基于窄帶隨機過程理論和Hermite模型,提出了一種預測船舶橫搖極值響應的方法.采用窄帶隨機過程理論主要是為了求出響應的窄帶分量,其原因是窄帶分量能夠更好地體現(xiàn)隨機響應極值的概率特性.接下來采用Hermite模型主要是為了解決響應過程的非高斯性問題.以C11集裝箱船為例,預測出某一特定海況下隨機參數(shù)橫搖極值響應的均值并與模型試驗結果進行了比較,本文提出的方法在極值預測上能夠給出偏于安全的預測值,可為實際工程提供參考.主要結論如下:
(1)傳統(tǒng)的Gumbel極值預測模型由于假設響應為高斯隨機過程,對C11船橫搖極值響應的預測誤差在18%左右,不適合預測參數(shù)橫搖響應極值.
(2)Monte Carlo預測的極值響應均值與模型試驗誤差很小,驗證了本文數(shù)學模型的準確性.
(3)本文預測方法在預測精度上相當于基于104樣本的數(shù)值計算,但是在計算成本上只有數(shù)值計算的1/500,驗證了本文方法的高效性.
(4)本文方法預測的極值響應比模型試驗結果大10%.推測其原因是忽略了響應極大值之間的相關性,這將會是作者后續(xù)研究的重點.