基于改進能量包線隨機平均法求解窄帶譜激勵下船舶隨機參激橫搖響應(yīng)

劉亞柳 劉利琴 劉依倫

（天津大學(xué) 水利工程仿真與安全國家重點實驗室 天津 300072）

0 引 言

船舶在遭遇惡劣海況時，為了避免橫風(fēng)橫浪的影響，通常會調(diào)整航向，選擇縱浪或斜浪航行。但即使船舶在靜水中有足夠的穩(wěn)性，在縱浪和斜浪航行時，仍有發(fā)生大角度搖擺乃至傾覆的可能。隨著船-波相對位置的周期性變化，船舶橫搖復(fù)原力矩也呈現(xiàn)出周期性變化，從而引發(fā)顯著的參數(shù)橫搖運動。國際海事組織（IMO）制定的第二代完整穩(wěn)性衡準(zhǔn)包含了參數(shù)激勵橫搖等5種失效模式，更加深入地研究了船舶在波浪中的動穩(wěn)性問題。

隨機波浪中的船舶大幅橫搖運動非常復(fù)雜。首先，它是一個強非線性問題，包括了以下非線性項：非線性恢復(fù)力、非線性阻尼力和非線性波浪載荷。其次，由于海浪的隨機特性，船舶運動響應(yīng)為隨機過程。因此，本文研究的隨機波浪中船舶橫搖運動是一個強非線性隨機動力學(xué)問題，不能用確定性方法描述，也不能采用常規(guī)的譜分析方法。研究隨機波浪中船舶橫搖運動響應(yīng)和穩(wěn)性最直接的方法是數(shù)值模擬，許多學(xué)者基于數(shù)值模擬對參數(shù)橫搖問題進行了大量研究。ESPARZA等采用安全池理論衡量了船舶在參數(shù)激勵下的抗傾覆能力，并分析了安全池形狀面積隨環(huán)境參數(shù)的變化以及破損現(xiàn)象。SILVA等基于切片理論提出了研究迎浪船舶動穩(wěn)性的新方法，并將該方法拓展到船舶六自由度非線性運動模型，通過與試驗結(jié)果對比，驗證了該方法的可行性。MATSUDA等考慮橫搖恢復(fù)力矩的變化研究了規(guī)則和隨機縱浪和斜浪中船舶的參數(shù)激勵橫搖運動，其中橫搖恢復(fù)力矩系數(shù)根據(jù)模型實驗測出。BELENKY等對頂浪中C11型集裝箱的參數(shù)橫搖運動進行了數(shù)值模擬，并采用統(tǒng)計方法和傅里葉變化得到了響應(yīng)譜密度和自相關(guān)函數(shù)。

船舶傾覆是十分罕見的事故，因此解析方法對于評估船舶傾覆概率，衡量船舶橫搖穩(wěn)定性仍然十分必要。在應(yīng)用隨機平均法求解隨機響應(yīng)時，一般要求隨機激勵的功率譜為寬帶譜。朱位秋指出，隨機平均法的基本思想是將隨機運動響應(yīng)近似為擴散的馬爾可夫過程。隨機平均法的準(zhǔn)確性在很大程度上取決于激勵譜帶寬與響應(yīng)譜帶寬的比值。但在實際工程應(yīng)用中，系統(tǒng)的實際響應(yīng)帶寬并不是很大，因此也可以采用隨機平均法求解非白噪聲激勵下的響應(yīng)。ROBERTS采用隨機平均法求解了有色噪聲激勵下的橫搖運動方程，基于Fokker-Planck方程推導(dǎo)出了橫搖穩(wěn)態(tài)響應(yīng)的概率分布。然而，這種隨機平均方法僅能反映系統(tǒng)的非線性阻尼特性。當(dāng)船舶具有強非線性恢復(fù)力時，該方法在大橫搖角處表現(xiàn)出明顯的誤差。另一方面，MARUYAMA等采用能量包線隨機平均法借助Hamiltonia量則可以同時考慮非線性阻尼和非線性恢復(fù)力。DOSTAL等采用能量包線隨機平均法求解了以RORO船的橫搖穩(wěn)態(tài)概率密度函數(shù)。柴威等采用能量包線隨機平均法和路徑積分法研究了隨機海浪中船舶非線性橫搖運動響應(yīng)。研究表明：隨機平均法的基本思想是對原始系統(tǒng)進行降維處理，從而降低隨機響應(yīng)的計算難度；相比之下，路徑積分法則是基于系統(tǒng)的馬爾可夫性質(zhì)，在時間域內(nèi)求解相應(yīng)的Fokker-Planck方程，得到隨時間演變的船舶橫搖運動的轉(zhuǎn)移概率密度。

雖然有一些學(xué)者成功采用能量包線隨機平均法求解了某些特定船舶的橫搖運動響應(yīng)，但本文作者以C11型集裝箱船為例，研究高海況下的船舶橫搖運動響應(yīng)時，采用能量包線隨機平均法得到的結(jié)果與數(shù)值模擬的結(jié)果有明顯差異。MARUYAMA等采用該方法求解ITTC A1型集裝箱船的橫搖運動響應(yīng)時也出現(xiàn)了類似現(xiàn)象。本課題組研究了外部激勵強度、激勵帶寬和系統(tǒng)阻尼對能量包線隨機平均法計算精度的影響，給出了該方法求解窄帶海浪激勵下船舶橫搖運動的適用性條件。當(dāng)采用能量包線隨機平均法求解非白噪聲激勵的系統(tǒng)響應(yīng)時，要求系統(tǒng)的阻尼和外部激勵強度足夠小且激勵帶寬足夠大。如果不滿足這些條件，則需要采用數(shù)值方法來計算橫搖響應(yīng)，以獲得更準(zhǔn)確的結(jié)果。隨機平均法分為標(biāo)準(zhǔn)隨機平均法和能量包線隨機平均法。標(biāo)準(zhǔn)隨機平均法雖然應(yīng)用很廣泛但在平均過程中會損失掉非線性恢復(fù)力項從而不能計及其對響應(yīng)的影響。而能量包線隨機平均法雖然能彌補標(biāo)準(zhǔn)隨機平均法的這個缺陷，但一般要求外界激勵為理想白噪聲。船舶與海洋工程結(jié)構(gòu)物遭遇的實際海浪一般具有有限的帶寬以及給定的譜密度函數(shù)。當(dāng)外界激勵為窄帶譜，并采用非線性阻尼和非線性恢復(fù)力描述船舶橫搖運動方程時，2種隨機平均法的適用性范圍受到限制。

本文基于這個背景，通過改進能量包線隨機平均法求解非線性橫搖系統(tǒng)在窄帶隨機波浪激勵下的響應(yīng)。以C11型集裝箱船為例，考慮非線性阻尼力和非線性恢復(fù)力建立橫搖單自由度運動方程，分析其在隨機波浪下的運動響應(yīng)特性。采用BULIAN提出的改進Grim有效波理論來描述隨機波面升高以降低解析方法的求解難度，即用固定波長的規(guī)則波在一個船長范圍內(nèi)對不規(guī)則波波面進行最小二乘替換，將波幅和船波相對位置作為隨機過程。假設(shè)船舶在波浪中的運動為準(zhǔn)靜態(tài)過程，基于切片理論數(shù)值計算橫搖復(fù)原力矩，并采用解析表達式進行擬合。采用改進能量包線隨機平均法求解隨機激勵下船舶的橫搖運動響應(yīng)，得到橫搖角的概率分布，更直觀地展示船舶的穩(wěn)性特征。采用蒙特卡洛法數(shù)值計算橫搖響應(yīng)分布，驗證改進能量包線隨機平均法的準(zhǔn)確性。

1 參數(shù)激勵橫搖運動方程

式中：為橫搖角，rad；I為橫搖慣性矩，kg·m；A(ω)為橫搖附加慣性矩，kg·m；ω為橫搖固有頻率，rad/s；為線性阻尼系數(shù)，N·m·s；為立方非線性阻尼系數(shù)，N·m·s；Δ為排水量，kg；為重力加速度，m/s；(，，)為船舶復(fù)原力臂，m，由()、()和()決定，其中()為波幅，m，()∈［0，2π］為船和波的相對位置，rad；為時間，s。

本文采用Grim有效波理論描述隨機波面升高，因此整個船長范圍內(nèi)隨機波面，可以用波長等于船長的規(guī)則波進行擬合。在波長等于船長的規(guī)則波中，基于切片理論和Froude-Krylov假設(shè)數(shù)值計算不同、和時C11集裝箱船的復(fù)原力臂函數(shù)(，，)。采用高斯積分求解作用在船體上的流體力，考慮瞬時濕表面計算靜水恢復(fù)力和F-K力。復(fù)原力求解示意圖如圖1所示，船舶橫搖復(fù)原力矩為：

圖1 復(fù)原力臂計算示意圖

式中：為海水密度，kg/m；為船長，m；y()與z()為各浸水橫剖面形心在參考坐標(biāo)系下的坐標(biāo)，m，參考坐標(biāo)系以船舶重心為原點，不隨船的運動而旋轉(zhuǎn)。()為各浸水橫剖面面積，m；為航向角，°；()為各橫剖面的壓力梯度系數(shù)。

為了使用解析方法求解運動方程（1），將(，，)展開為關(guān)于橫搖角的多項式：

式中：q (1，2，3)為展開項的系數(shù)，使用最小二乘法確定。

采用JONSWAP波浪譜對隨機海況進行描述：

式中：為能量尺度度量；為譜峰升高因子；ω為譜峰頻率，rad/s，ω=2π/T，T為譜峰周期，s；H為有義波高，m；為峰形參數(shù)。當(dāng)波浪頻率ω＜ω時，0.07；當(dāng)ω＞ω時，0.09。

圖2是根據(jù)有義波高和譜峰周期的聯(lián)合概率密度繪制的波浪散布圖。它表示100 000次海浪觀測中不同有義波高和譜峰周期組合的數(shù)量。

圖2 波浪散布圖

本文以C11型集裝箱船為例，采用改進能量包線隨機平均法求解其在隨機波浪中的橫搖運動響應(yīng)。C11型集裝箱船的主尺度參數(shù)如表1所示。

表1 C11型集裝箱船主尺度

2 改進能量包線隨機平均法

式中：

并作如下尺度變化：

為了書寫簡便，以下用替代，則運動微分方程（5）變?yōu)椋?/p>

系統(tǒng)總能量的具體表達式為：

基于原始能量包線隨機平均法計算不同海況下的橫搖角概率分布，并與蒙特卡洛數(shù)值模擬得到的結(jié)果進行對比。采用Kullback-Leibler（KL）散度描述能量包線隨機平均法的橫搖角度概率分布于蒙特卡羅法計算結(jié)果的差異。如果2個分布相同，KL散度為0。當(dāng)2個分布之間的差異增加時，KL散度也會增加。

從圖3可以直觀地看出，隨著有義波高和峰值因子的增大，原始能量包線隨機平均法計算的準(zhǔn)確性降低。因此，采用原始能量包線隨機平均法計算高海況下船舶橫搖運動響應(yīng)時，存在很大誤差，嚴重影響橫搖穩(wěn)性評估。

圖3 不同海況下的KL散度

基于以上背景，對原始能量包線隨機平均法進行改進。借鑒標(biāo)準(zhǔn)隨機平均法引入相角，進行坐標(biāo)變換：

公式（8）兩邊對時間求導(dǎo)可得：

公式（7）下式的左右兩邊同時乘以可得：

聯(lián)立公式（10）和公式（11）可得：

由公式（9）可知：

則：

根據(jù)公式（9）可以將、表示為和的函數(shù)，代入公式（12）和公式（16）中可以得到標(biāo)準(zhǔn)形式的隨機微分方程為：

式中：

式中：()為ξ的自相關(guān)函數(shù)；

式中：()和()為相互獨立的維納過程。

對于上述的平均伊藤隨機微分方程（20），其響應(yīng)為擴散過程，系統(tǒng)能量對應(yīng)的轉(zhuǎn)移概率密度(，，)滿足如下的FPK方程：

在船舶不發(fā)生傾覆的前提下，式（21）對應(yīng)的平穩(wěn)概率密度函數(shù)為：

式中：參數(shù)由邊界條件確定。

根據(jù)能量的平穩(wěn)概率密度函數(shù)，可以很容易地得到橫搖角和橫搖角速度的聯(lián)合概率密度函數(shù)。其對應(yīng)的解析表達式為：

對橫搖角速度進行積分，進而得到橫搖角的概率密度分布。

打底層為了避免熔池在正面焊縫產(chǎn)生墜瘤，坡口兩側(cè)熔合良好，焊接過程中焊條向坡口根部用力送進，并在坡口兩側(cè)稍許停留，以保證背面的穿透和坡口兩側(cè)的熔合性。由于焊條電弧焊時焊條的電阻大、焊條易發(fā)紅，后半部分焊條的使用性能較差，所以在焊接過程中應(yīng)留有稍長的焊條頭，以5cm為宜，避免出現(xiàn)焊條沾粘、保護不良等焊接缺陷。

3 計算結(jié)果分析

本文重點研究高海況下采用改進能量包線隨機平均法預(yù)報船舶大幅非線性橫搖運動。因此，選取海況使船舶遭遇譜的譜峰頻率在橫搖固有頻率的2倍附近，以激發(fā)參數(shù)激勵橫搖運動。采用JONSWAP波浪譜對隨機海況進行描述，=1。船舶在隨機波浪中頂浪航行，具體計算工況如表2所示。

表2 計算工況

隨船坐標(biāo)系下遭遇頻率ω的表達式為：

式中：為波數(shù)；為航速，kn。

遭遇頻率ω下的海浪譜(ω)的表達式為：

采用四階龍格庫塔方法對方程（5）進行數(shù)值計算，并采用蒙特卡洛法統(tǒng)計橫搖角的概率分布，隨機參數(shù)激勵下的橫搖角具有普遍非各態(tài)歷經(jīng)性。為了克服這一現(xiàn)象，需要在同一海浪譜下進行多次不同的實現(xiàn)。因此，選取 500個隨機種子數(shù)，每次數(shù)值模擬的時長為50 000 s。將解析方法得到的橫搖角概率密度分布與蒙特卡洛法得到的結(jié)果進行比較以驗證改進能量包線隨機平均法的準(zhǔn)確性。

不同工況下原始和改進能量包線隨機平均法的計算結(jié)果以及和蒙特卡洛法的對比如圖4至圖7所示。

圖4 工況1計算結(jié)果

圖7 工況4計算結(jié)果

圖5 工況2計算結(jié)果

圖6 工況3計算結(jié)果

由圖4至圖7的計算結(jié)果可知，改進能量包線隨機平均法能夠大幅提升隨機平均法的準(zhǔn)確度，計算高海況下船舶參數(shù)激勵橫搖運動響應(yīng)分布，對船舶穩(wěn)性進行直接評估。但是從對數(shù)坐標(biāo)軸中能夠看出，橫搖角概率密度曲線在尾部區(qū)域仍存在一定的計算誤差。

4 結(jié) 語

本文提出改進能量包線隨機平均法，在原有的能量包線隨機平均法的基礎(chǔ)上借鑒標(biāo)準(zhǔn)隨機平均法引入相角。將非線性系統(tǒng)對隨機激勵的響應(yīng)近似看成擴散的馬爾柯夫過程，推導(dǎo)能量與相角滿足的伊藤隨機微分方程，將漂移系數(shù)和擴散系數(shù)在快變量的一個周期內(nèi)進行平均，得到支配系統(tǒng)能量的轉(zhuǎn)移概率密度所滿足的FPK方程，其中快變量對慢變量的影響由于其隨時間迅速變化的性質(zhì)被平均掉。改進的能量包線隨機平均法既可以考慮非線性恢復(fù)力，又同時放寬了外部激勵帶寬的限制，使其可以為窄帶隨機激勵。最后通過求解平均后的 FPK方程得到橫搖角穩(wěn)態(tài)響應(yīng)概率分布。

改進能量包線隨機平均法能大幅提升隨機平均法的計算精度，在高海況下計算得到的橫搖運動響應(yīng)分布與蒙特卡洛法對比良好。但是從對數(shù)坐標(biāo)軸中能夠看出，橫搖角概率密度曲線的尾部區(qū)域仍存在一定的計算誤差。由于大橫搖角區(qū)域?qū)τ跍?zhǔn)確預(yù)報船舶橫搖傾覆概率或者首次穿越時間非常重要，因此后續(xù)仍然需要研究如何提高尾部區(qū)域預(yù)報的準(zhǔn)確性。