卜淑霞，魯 江，顧 民，吳乘勝，儲紀龍

（中國船舶科學研究中心，江蘇省綠色船舶技術重點實驗室，江蘇 無錫 214082）

0 引 言

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

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

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

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

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

1 數(shù)學模型

1.1 三維時域混合源法

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

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

假設流體無粘、無旋和不可壓縮，水深為無限水深，則流場非定常的速度勢可表示為：

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

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

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

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

在外域，記外域總擾動勢 Φ（ p,t）為 ΦII（p,t），那么 ΦII（p,t）應滿足以下條件：

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

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

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

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

對船體平均濕表面積Sb上的源強積分即可得到船體濕表面上的擾動速度勢ΦI，且已知入射波速度勢Φw，參見公式（9），最后通過伯努利方程可得到相應的壓力項：

求得每個面元控制點的壓力后，對每個面元積分即可求得作用于該面元上的流體作用力F和力矩M。

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

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

將波浪頻譜ω等分成N份，此時N個波浪譜對應的波幅可表示為：

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

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

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

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

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

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

在統(tǒng)計分析中，隨機序列的平均值和標準差采用如下公式計算：

2 計算模型

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

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

該船型是參數(shù)橫搖研究的國際標模，大阪大學的Umeda和Hashimoto等人基于該船型開展了不規(guī)則波中的參數(shù)橫搖模型試驗[13]，文中采用該模型試驗數(shù)據(jù)對數(shù)值模擬結(jié)果進行了對比分析。在頂浪狀態(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進行研究。針對表1中所示的狀態(tài)，作者曾開展過規(guī)則波中的模型試驗以及靜水、有航速下的自由橫搖衰減試驗[6]，因此數(shù)值模擬中的線性和立方阻尼采用模型試驗數(shù)據(jù)結(jié)果。

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

3 結(jié)果與分析

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

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

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

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

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

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

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

圖7 垂蕩幅值的標準偏差（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 有義橫搖幅值的標準偏差（H1/3=7.82 m，F(xiàn)n=0.1）Fig.9 Standard deviation of roll angles

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

3.2 統(tǒng)計分析

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

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

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

圖13 參數(shù)橫搖隨時間的變化（左：時歷曲線；右：最大橫搖幅值和標準差的平均值）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é) 論

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

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

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

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

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