劉亞柳 劉利琴 劉依倫

（天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室 天津 300072）

0 引 言

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

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

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

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

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

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

式中：為橫搖角，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；為時(shí)間，s。

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

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

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

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

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

采用JONSWAP波浪譜對隨機(jī)海況進(jìn)行描述：

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

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

圖2 波浪散布圖

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

表1 C11型集裝箱船主尺度

2 改進(jìn)能量包線隨機(jī)平均法

式中：

并作如下尺度變化：

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

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

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

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

圖3 不同海況下的KL散度

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

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

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

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

由公式（9）可知：

則：

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

式中：

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

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

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

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

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

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

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

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

3 計(jì)算結(jié)果分析

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

表2 計(jì)算工況

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

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

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

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

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

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

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

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

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

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

4 結(jié) 語

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

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