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

劉亞柳 劉利琴 劉依倫

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

0 引 言

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

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

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

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

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

1 參數激勵橫搖運動方程

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

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

圖1 復原力臂計算示意圖

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

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

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

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

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

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

圖2 波浪散布圖

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

表1 C11型集裝箱船主尺度

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

式中：

并作如下尺度變化：

為了書寫簡便，以下用替代，則運動微分方程（5）變為：

系統總能量的具體表達式為：

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

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

圖3 不同海況下的KL散度

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

公式（8）兩邊對時間求導可得：

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

聯立公式（10）和公式（11）可得：

由公式（9）可知：

則：

根據公式（9）可以將、表示為和的函數，代入公式（12）和公式（16）中可以得到標準形式的隨機微分方程為：

式中：

式中：()為ξ的自相關函數；

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

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

在船舶不發生傾覆的前提下，式（21）對應的平穩概率密度函數為：

式中：參數由邊界條件確定。

根據能量的平穩概率密度函數，可以很容易地得到橫搖角和橫搖角速度的聯合概率密度函數。其對應的解析表達式為：

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

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

3 計算結果分析

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

表2 計算工況

隨船坐標系下遭遇頻率ω的表達式為：

式中：為波數；為航速，kn。

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

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

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

圖4 工況1計算結果

圖7 工況4計算結果

圖5 工況2計算結果

圖6 工況3計算結果

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

4 結 語

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

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