基于三維時(shí)域混合源法的頂浪不規(guī)則波參數(shù)橫搖研究

卜淑霞，魯 江，顧 民，吳乘勝，儲(chǔ)紀(jì)龍

（中國(guó)船舶科學(xué)研究中心，江蘇省綠色船舶技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 無(wú)錫 214082）

0 引 言

船舶參數(shù)橫搖是波浪中復(fù)原力矩周期性變化引起的非線性大幅運(yùn)動(dòng)，是船舶在波浪中典型的傾覆現(xiàn)象（參數(shù)橫搖、純穩(wěn)性喪失和橫甩）之一，該現(xiàn)象也可能發(fā)生在比較溫和的海況下，因此對(duì)船舶的安全航行構(gòu)成了嚴(yán)重的威脅，也是國(guó)際海事組織（IMO）正在制定的船舶第二代完整穩(wěn)性衡準(zhǔn)中最重要的穩(wěn)性失效模式之一[1]。1998年，巴拿馬型C11集裝箱船APL CHINA號(hào)在北太平洋海域頂浪航行時(shí)遭遇到了嚴(yán)重的參數(shù)橫搖[2]，橫搖角高達(dá)40°，此次事故引起了船舶界對(duì)頂浪參數(shù)橫搖的關(guān)注。其后一艘汽車(chē)運(yùn)輸船在亞速爾群島水域頂浪航行時(shí)也遭遇到嚴(yán)重的參數(shù)橫搖[3]，這些船舶都具有較大的外飄船首和較寬的船尾，使得波浪中穩(wěn)性的變化較大。這兩起事故使參數(shù)橫搖的研究成為熱點(diǎn)，并將參數(shù)橫搖的研究重點(diǎn)轉(zhuǎn)移到頂浪航行狀態(tài)。對(duì)于規(guī)則波中參數(shù)橫搖的數(shù)值預(yù)報(bào)，目前常用的預(yù)報(bào)方法有二維切片法和三維面元法[2，4-5]，本文作者也分別基于二維切片法和三維時(shí)域面元法進(jìn)行了頂浪規(guī)則波中參數(shù)橫搖的數(shù)值預(yù)報(bào)，預(yù)報(bào)結(jié)果與模型試驗(yàn)吻合良好[6-7]。

參數(shù)橫搖的發(fā)生需要滿足特定的參數(shù)條件，不規(guī)則波中參數(shù)橫搖動(dòng)態(tài)系統(tǒng)的非線性使得參數(shù)橫搖的概率分布特征與常規(guī)運(yùn)動(dòng)有所不同[8-9]，這些特征顯示出各態(tài)歷經(jīng)性的假設(shè)不適用于該現(xiàn)象，導(dǎo)致不規(guī)則波中參數(shù)橫搖的定量預(yù)報(bào)十分困難。由于隨機(jī)波浪符合正態(tài)分布和各態(tài)歷經(jīng)，根據(jù)Weiner-Khinchin的理論，則船體的運(yùn)動(dòng)響應(yīng)也符合正態(tài)分布和各態(tài)歷經(jīng)的特征，這也是大部分耐波性和穩(wěn)性計(jì)算方法的理論基礎(chǔ)，但這些方法不能有效地衡量不規(guī)則波中的參數(shù)橫搖。

對(duì)于不規(guī)則波中參數(shù)橫搖的數(shù)值模擬方法，歐盟SAFEDOR項(xiàng)目[10]曾對(duì)14種勢(shì)流方法進(jìn)行了規(guī)則波和不規(guī)則波中參數(shù)橫搖數(shù)值預(yù)報(bào)的試驗(yàn)基準(zhǔn)研究，包括切片法和面元法，預(yù)報(bào)結(jié)果并不理想。Bulian等[11]針對(duì)巴拿馬集裝箱船開(kāi)展了系列不規(guī)則波參數(shù)橫搖模型試驗(yàn)，試驗(yàn)表明參數(shù)橫搖時(shí)歷曲線統(tǒng)計(jì)的不確定性比較大，表現(xiàn)出“非各態(tài)歷經(jīng)性”的特點(diǎn)。Umeda等[12]考慮了波浪增阻對(duì)頂浪不規(guī)則波參數(shù)橫搖的影響，基于切片理論對(duì)不規(guī)則波中的參數(shù)橫搖進(jìn)行了預(yù)報(bào)。Hashimoto等[13]針對(duì)C11集裝箱船進(jìn)行了不規(guī)則波中參數(shù)橫搖的模型試驗(yàn)和數(shù)值預(yù)報(bào)，對(duì)垂蕩-縱搖-橫搖相互耦合的3DOF數(shù)學(xué)模型進(jìn)行了研究，研究表明3DOF數(shù)學(xué)模型可以較好地預(yù)報(bào)頂浪不規(guī)則波中的參數(shù)橫搖。魯江等[14]基于二維切片法，考慮了波浪增阻影響，對(duì)頂浪不規(guī)則波參數(shù)橫搖開(kāi)展了數(shù)值計(jì)算研究。

本文采用三維時(shí)域混合源法進(jìn)行頂浪不規(guī)則波中參數(shù)橫搖的研究。三維時(shí)域混合源法，即內(nèi)域采用Rankine源，在外域采用瞬態(tài)的時(shí)域Green函數(shù)，這樣可以保留Rankine源易于計(jì)算、可以得到近場(chǎng)定常速度勢(shì)，以及Green函數(shù)法僅需在物體表面進(jìn)行離散，函數(shù)自動(dòng)滿足線性自由表面和遠(yuǎn)場(chǎng)輻射條件的優(yōu)點(diǎn)，消除了兩者的缺點(diǎn)，在計(jì)算非線性大幅運(yùn)動(dòng)時(shí)具有明顯的優(yōu)勢(shì)[15]。

文中首先基于混合源法求解頂浪不規(guī)則波中的輻射勢(shì)和繞射勢(shì)，對(duì)船體平均濕表面積上的源強(qiáng)積分得到船體濕表面上的擾動(dòng)速度勢(shì)，通過(guò)伯努利方程得到水動(dòng)壓力；其次沿不規(guī)則波中瞬時(shí)濕表面對(duì)壓力積分求解出FK力和靜水力；然后通過(guò)垂蕩-縱搖-橫搖三自由度耦合方程求解參數(shù)橫搖。最后以國(guó)際標(biāo)模C11集裝箱船為研究對(duì)象，進(jìn)行了頂浪不規(guī)則波中的參數(shù)橫搖研究，探索了三維時(shí)域混合源法計(jì)算不規(guī)則波參數(shù)橫搖的有效性。

1 數(shù)學(xué)模型

1.1 三維時(shí)域混合源法

混合源法在數(shù)值求解中引入了控制面SC，將流場(chǎng)分為內(nèi)域I和外域II，內(nèi)域I是由船體濕表面Sb、部分自由液面Sf1和控制面SC包圍的閉合區(qū)域，流場(chǎng)分布和船體表面網(wǎng)格分布如圖1所示。

圖1 流場(chǎng)區(qū)域劃分和船體表面網(wǎng)格示意圖Fig.1 Domain definitions and meshes schematic

假設(shè)流體無(wú)粘、無(wú)旋和不可壓縮，水深為無(wú)限水深，則流場(chǎng)非定常的速度勢(shì)可表示為：

其中：Φw（p,t）為入射波速度勢(shì)，是由于船體存在引起的總擾動(dòng)勢(shì)，由于入射波是已知的，則只要求解由繞射勢(shì)和輻射勢(shì)組成的總擾動(dòng)勢(shì)即可。

在內(nèi)域，記內(nèi)域總擾動(dòng)勢(shì) Φ（ p,t）為 ΦI（p,t），那么 ΦI（p,t）應(yīng)該滿足以下條件：

則內(nèi)域I中Rankine源的邊界積分方程如下：

其中：ΦI是內(nèi)域I總擾動(dòng)速度勢(shì)，G=1/rpq為簡(jiǎn)單格林函數(shù)，p（ x,y, z ）為場(chǎng)點(diǎn)，q （ ξ,η, ζ）為源點(diǎn)，rpq=

在外域，記外域總擾動(dòng)勢(shì) Φ（ p,t）為 ΦII（p,t），那么 ΦII（p,t）應(yīng)滿足以下條件：

為了求解該定解問(wèn)題，我們引入如下的時(shí)域格林函數(shù)：

外域II中使用時(shí)域格林函數(shù)，面元分布在控制面SC上，邊界積分方程如下：

其中：ΦII為外域II總擾動(dòng)速度勢(shì)，w（τ）是控制面的水線面，VN是w（τ）的法向速度。

控制面隨船體一起運(yùn)動(dòng)，因此在控制面上內(nèi)外域連續(xù)，采用面元法對(duì)邊界積分方程（3）和（6）進(jìn)行數(shù)值離散，可以獲得當(dāng)前時(shí)刻船體濕表面積Sb上的ΦI，自由表面Sf1上的以及控制面SC上的ΦI和。然后就可以利用物面上的ΦI，通過(guò)伯努利方程計(jì)算船體表面的壓力以及相應(yīng)的水動(dòng)力，利用內(nèi)域I中線性自由面獲得下一時(shí)刻整個(gè)流場(chǎng)的擾動(dòng)勢(shì)和下一時(shí)刻內(nèi)域的速度勢(shì)ΦI。

對(duì)船體平均濕表面積Sb上的源強(qiáng)積分即可得到船體濕表面上的擾動(dòng)速度勢(shì)ΦI，且已知入射波速度勢(shì)Φw，參見(jiàn)公式（9），最后通過(guò)伯努利方程可得到相應(yīng)的壓力項(xiàng)：

求得每個(gè)面元控制點(diǎn)的壓力后，對(duì)每個(gè)面元積分即可求得作用于該面元上的流體作用力F和力矩M。

在不規(guī)則波中入射波的速度勢(shì)采用線性疊加的形式：

其中：Ai、ki、ωi、εi和 βi分別是長(zhǎng)峰不規(guī)則波中第 i個(gè)規(guī)則諧波的振幅、波數(shù)、頻率、相位和浪向角，εi=2πPi（Pi是 0～1 之間的隨機(jī)數(shù)）。

將波浪頻譜ω等分成N份，此時(shí)N個(gè)波浪譜對(duì)應(yīng)的波幅可表示為：

將波形疊加，得到不規(guī)則波的波面方程，表示如下：

不規(guī)則波浪譜采用ITTC雙參數(shù)譜：

其中：H1/3為有義波高；T01為波浪特征周期；ω為波浪圓頻率。

1.2 參數(shù)橫搖數(shù)學(xué)模型

根據(jù)IMO船舶第二代完整穩(wěn)性最新提案[1]的框架要求，本文選取了垂蕩—縱搖—橫搖相互耦合的3DOF數(shù)學(xué)模型。

其中：m為船舶質(zhì)量；Ixx為橫搖慣性矩；Iyy為縱搖慣性矩；Aij、Bij為附加質(zhì)量和阻尼系數(shù)；ζ為垂蕩位移；θ為縱搖；φ為橫搖；N1、N3為線性和立方的橫搖阻尼系數(shù)，文中采用模型試驗(yàn)結(jié)果；FFK+H為FK力和靜水力，通過(guò)對(duì)不規(guī)則波表面的瞬時(shí)濕表面壓力積分得到；FDF為繞射力，沿船體平均濕表面積分得到。船體運(yùn)動(dòng)的偏微分方程利用Runge-Kutta方法求解。

在統(tǒng)計(jì)分析中，隨機(jī)序列的平均值和標(biāo)準(zhǔn)差采用如下公式計(jì)算：

2 計(jì)算模型

選取C11集裝箱船作為研究對(duì)象，主尺度如表1所示，模型縮尺比為1：65.5，幾何外形和船體型線圖如圖2所示。

表1 C11集裝箱船主尺度（縮尺比1/65.5）Tab.1 Principal particulars of the C11 containership(Scale:1/65.5)

該船型是參數(shù)橫搖研究的國(guó)際標(biāo)模，大阪大學(xué)的Umeda和Hashimoto等人基于該船型開(kāi)展了不規(guī)則波中的參數(shù)橫搖模型試驗(yàn)[13]，文中采用該模型試驗(yàn)數(shù)據(jù)對(duì)數(shù)值模擬結(jié)果進(jìn)行了對(duì)比分析。在頂浪狀態(tài)下， 選取工況 T01=9.99 s、H1/3=7.82 m、Fn=0.0，T01=9.99 s、H1/3=10.43 m、Fn=0.0，T01=9.99 s、H1/3=7.82 m和Fn=0.1進(jìn)行研究。針對(duì)表1中所示的狀態(tài)，作者曾開(kāi)展過(guò)規(guī)則波中的模型試驗(yàn)以及靜水、有航速下的自由橫搖衰減試驗(yàn)[6]，因此數(shù)值模擬中的線性和立方阻尼采用模型試驗(yàn)數(shù)據(jù)結(jié)果。

圖2 C11集裝箱船（左：幾何外形；右：型線）Fig.2 C11 containership(Left：hull geometry;Right：body lines)

3 結(jié)果與分析

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

模型試驗(yàn)中零航速和有航速的工況均采用了不同的隨機(jī)種子序列，因此，數(shù)值模擬中也采用不同的隨機(jī)種子序列，以更好地匹配模型試驗(yàn)結(jié)果。數(shù)值模擬中部分不規(guī)則波波譜與ITTC目標(biāo)譜的對(duì)比如圖3所示，可以看出文中采用的方法和目標(biāo)譜吻合良好。

圖3 不同隨機(jī)序列下的波浪譜Fig.3 Wave spectra under different realization numbers

工況Fn=0.0、H1/3=7.82 m時(shí)最大橫搖幅值和橫搖偏差對(duì)比結(jié)果如圖4-5所示，從圖中可以看出，文中數(shù)值模擬結(jié)果較好地重現(xiàn)了模型試驗(yàn)結(jié)果。從圖5所示的該工況下的橫搖幅值偏差可以看出，橫搖運(yùn)動(dòng)的標(biāo)準(zhǔn)偏差比較分散，這是由于參數(shù)橫搖具有實(shí)際非各態(tài)歷經(jīng)特性。該工況下的垂蕩和縱搖的標(biāo)準(zhǔn)偏差如圖6所示，由于文獻(xiàn)中沒(méi)有給出垂蕩運(yùn)動(dòng)的試驗(yàn)值，因此，垂蕩運(yùn)動(dòng)僅有數(shù)值模擬結(jié)果。從圖中可以看出，垂蕩和縱搖運(yùn)動(dòng)的標(biāo)準(zhǔn)偏差比較集中，呈現(xiàn)各態(tài)歷經(jīng)的特征，與橫搖運(yùn)動(dòng)有所不同。

圖5 橫搖幅值的標(biāo)準(zhǔn)偏差（H1/3=7.82 m，F(xiàn)n=0.0）Fig.5 Standard deviation of roll angles

有航速時(shí)工況Fn=0.1、H1/3=7.82 m對(duì)應(yīng)的最大橫搖幅值對(duì)比結(jié)果如圖8所示，從圖中可以看出，文中采用的方法也可以較好地模擬有航速時(shí)頂浪不規(guī)則波中的參數(shù)橫搖。對(duì)于該工況文獻(xiàn)中未給出標(biāo)準(zhǔn)偏差的試驗(yàn)結(jié)果，但從圖9所示的標(biāo)準(zhǔn)偏差的計(jì)算結(jié)果也可以看出，有航速時(shí)橫搖運(yùn)動(dòng)的標(biāo)準(zhǔn)偏差也比較分散，呈現(xiàn)非各態(tài)歷經(jīng)的特征，與零航速的特征一致。

圖6 縱搖幅值的標(biāo)準(zhǔn)偏差（H1/3=7.82 m，F(xiàn)n=0.0）Fig.6 Standard deviation of pitch motions

圖7 垂蕩幅值的標(biāo)準(zhǔn)偏差（H1/3=7.82 m，F(xiàn)n=0.0）Fig.7 Standard deviation of heave motions

圖8 最大橫搖幅值（H1/3=7.82 m，F(xiàn)n=0.1）Fig.8 Maximum roll angles of parametric roll

圖9 有義橫搖幅值的標(biāo)準(zhǔn)偏差（H1/3=7.82 m，F(xiàn)n=0.1）Fig.9 Standard deviation of roll angles

零航速時(shí)另一工況Fn=0.0、H1/3=10.43 m對(duì)應(yīng)的最大橫搖幅值的平均值計(jì)算結(jié)果和模型試驗(yàn)對(duì)比如圖10所示，也可以證明本文計(jì)算方法的可靠性，從圖中也可以看出，H1/3=10.43 m時(shí)的橫搖幅值要比H1/3=7.82 m時(shí)大，表明該工況更容易發(fā)生參數(shù)橫搖。

3.2 統(tǒng)計(jì)分析

圖10 最大橫搖幅值的平均值（H1/3=10.43 m，F(xiàn)n=0.0）Fig.10 Average of maximum roll angles

圖11 概率密度分布（Fn=0.0，H1/3=10.43 m）Fig.11 Probability density function for Fn=0.0 H1/3=10.43 m

圖12 概率密度分布（Fn=0.1，H1/3=7.82 m）Fig.12 Probability density function for Fn=0.1 H1=7.82 m

統(tǒng)計(jì)分析重點(diǎn)研究不規(guī)則波中瞬時(shí)角度的分布概率，零航速下統(tǒng)計(jì)概率如圖11所示，有航速下的統(tǒng)計(jì)概率如圖12所示，從圖中可以看出，數(shù)值模擬的ITTC雙參數(shù)波浪譜符合高斯分布，這與實(shí)際情況相符。在數(shù)值模擬中，雖然采用了非線性的Froude-Krylov力和靜水力，但縱搖和垂蕩運(yùn)動(dòng)仍然符合高斯分布，此時(shí)可采用基于線性理論的方法計(jì)算。

但從圖中可以看出，橫搖運(yùn)動(dòng)分布與常規(guī)運(yùn)動(dòng)分布明顯不同：概率密度的峰值具有顯著的尖點(diǎn)，明顯不符合高斯分布，呈現(xiàn)非各態(tài)歷經(jīng)的特點(diǎn)，說(shuō)明參數(shù)橫搖具有顯著的非線性特征。橫搖角度為0°附近的概率密度較大，也就說(shuō)明在整個(gè)時(shí)歷過(guò)程中，不發(fā)生橫搖運(yùn)動(dòng)的概率較大。橫搖運(yùn)動(dòng)不符合高斯分布主要是由于在頂浪狀態(tài)下，船體遭遇波浪后會(huì)立即產(chǎn)生垂蕩和縱搖運(yùn)動(dòng)，但是僅當(dāng)遭遇的條件超出參數(shù)橫搖發(fā)生的閾值后，才會(huì)產(chǎn)生橫搖運(yùn)動(dòng)。該現(xiàn)象說(shuō)明了基于線性或者弱非線性假設(shè)的常規(guī)運(yùn)動(dòng)方程可以用于垂蕩和縱搖的計(jì)算，但是不能用于參數(shù)橫搖的計(jì)算，結(jié)論和文獻(xiàn)[9]的結(jié)論一致。

在數(shù)值模擬中也發(fā)現(xiàn)，橫搖運(yùn)動(dòng)的統(tǒng)計(jì)特性與模擬時(shí)間有關(guān)，選取工況Fn=0.0、H1/3=10.43 m進(jìn)行研究。從圖13（右）所示的某隨機(jī)序列下的時(shí)歷曲線可以看出，短時(shí)間內(nèi)可能不能展現(xiàn)所有的參數(shù)橫搖特征。對(duì)所有隨機(jī)序列進(jìn)行足夠時(shí)間的模擬（t＞8 000 s），然后分別統(tǒng)計(jì)t=2 000 s、4 000 s、6 000 s和8 000 s時(shí)最大橫搖幅值的平均值以及標(biāo)準(zhǔn)偏差的平均值，統(tǒng)計(jì)結(jié)果如圖13（左）所示，從結(jié)果可以看出，隨著時(shí)間的增加，橫搖運(yùn)動(dòng)的統(tǒng)計(jì)特性趨于穩(wěn)定，也就說(shuō)明參數(shù)橫搖的數(shù)值模擬需要足夠的隨機(jī)序列以及足夠的模擬時(shí)間。

圖13 參數(shù)橫搖隨時(shí)間的變化（左：時(shí)歷曲線；右：最大橫搖幅值和標(biāo)準(zhǔn)差的平均值）Fig.13 The parametric roll as a function of time(Left：time history of parametric roll under one realization number;Right：average of maximum roll angle and standard deviation)

4 結(jié) 論

文中提出了一種采用三維混合源法預(yù)報(bào)船舶頂浪不規(guī)則波中參數(shù)橫搖的方法，并選取國(guó)際標(biāo)模C11集裝箱船進(jìn)行了研究，得出如下結(jié)論：

（1）參數(shù)橫搖的發(fā)生具有非各態(tài)歷經(jīng)的特點(diǎn)，且單個(gè)隨機(jī)序列的統(tǒng)計(jì)結(jié)果與計(jì)算時(shí)間有關(guān)，時(shí)間足夠長(zhǎng)后，統(tǒng)計(jì)特性趨于穩(wěn)定，因此，不規(guī)則波中參數(shù)橫搖的數(shù)值模擬需要足夠的隨機(jī)序列以及足夠長(zhǎng)的模擬時(shí)間；

（2）參數(shù)橫搖具有強(qiáng)非線性的特征，基于線性或者弱非線性的理論不能用于參數(shù)橫搖的數(shù)值預(yù)報(bào)，因此，需要改進(jìn)目前常規(guī)耐波性和穩(wěn)性計(jì)算模型；

（3）三維時(shí)域面元法可以較好地預(yù)報(bào)頂浪不規(guī)則波中的參數(shù)橫搖，并且符合IMO參數(shù)橫搖穩(wěn)性直接評(píng)估框架的要求，因此，可用于IMO參數(shù)橫搖穩(wěn)性直接評(píng)估。

不規(guī)則波中的輻射力和繞射力是數(shù)值求解的難點(diǎn)，文中求解了不規(guī)則波中無(wú)橫傾時(shí)的輻射力和繞射力，后續(xù)需要進(jìn)一步研究不規(guī)則波中不同時(shí)刻不同橫傾時(shí)的輻射力和繞射力對(duì)參數(shù)橫搖的影響。