曲德才 于立偉 王樹青 馬 寧

（1.中國海洋大學(xué) 山東省海洋工程重點實驗室 青島 266101； 2.上海交通大學(xué) 國家重點海洋工程實驗室 上海 200030）

0 引 言

惡劣的尾隨浪中，在螺旋槳推力、阻力和波浪力的共同作用下會存在船舶縱蕩運動的平衡點。在一定的波浪和初始條件下，船舶可能處于波谷附近的穩(wěn)定平衡點，此時船速加速到波速，這被稱為“騎浪”。騎浪普遍被認(rèn)為是橫甩的先決條件，而橫甩會導(dǎo)致不可控的大幅艏搖和橫搖，甚至造成船舶傾覆。

騎浪橫甩運動作為一種強非線性現(xiàn)象，其在不規(guī)則波中的概率統(tǒng)計分析無法適用線性疊加原理。GRIM最早開展了不規(guī)則尾隨浪下騎浪運動的概率統(tǒng)計研究，并將騎浪的發(fā)生與航速超過臨界速度的增速現(xiàn)象聯(lián)系在一起。UMEDA基于規(guī)則波中確定性騎浪的臨界值和波高/周期聯(lián)合概率分布，提出了一種直接的理論方法，來估計不規(guī)則波中的騎浪概率。UMEDA等進一步將這種方法應(yīng)用于不規(guī)則波中橫甩的概率估計，并通過蒙特卡羅法驗證。該方法采用了規(guī)則波中橫甩的確定性閾值作為不規(guī)則波中的橫甩判定條件。然后，UMEDA等開展了不規(guī)則波中的自航模實驗，以進一步驗證所提出的橫甩概率估計方法，結(jié)果表明橫甩概率估計與實驗結(jié)果的差異處在95%置信區(qū)間內(nèi)。

近年來，受規(guī)則波浪中騎浪發(fā)生時船速達到波速的現(xiàn)象啟發(fā)，研究人員提出將不規(guī)則波浪中騎浪的發(fā)生與船速超過瞬時波速的現(xiàn)象聯(lián)系起來，這種船速超過瞬時波速的現(xiàn)象稱為增速（high run）。這就首先需要闡明不規(guī)則波中瞬時波速(也稱為局部波速)的物理含義。SPYROU等和BELENKY等提出通過波浪的過零和恒定斜率點的速度來表征局部瞬時波速。同時，SPYROU等和THEMELIS等提出通過希爾伯特變換獲得瞬時波速的方法，并對不同的瞬時波速計算方法進行了比較，證明了這些方法的相似性。BELENKY等進一步研究了不規(guī)則波中增速的概率統(tǒng)計結(jié)果。在此基礎(chǔ)上，SPYROU等和THEMILIS等研究了橫甩概率及其與增速的相關(guān)性。通過對大量數(shù)值模擬的直接計數(shù)，擬合了增速和顯著艏搖累積時間直方圖，研究了橫搖對增速橫甩相關(guān)性的影響，并對航向、舵控制增益以及運動閾值進行了靈敏度研究。由此可見，對增速現(xiàn)象的研究有助于騎浪和橫甩發(fā)生概率的研究，而隨著IMO二代完整穩(wěn)性衡準(zhǔn)的發(fā)布試用，不規(guī)則波中騎浪和橫甩發(fā)生概率的研究至關(guān)重要。

因此，本文針對不規(guī)則波中的騎浪/橫甩穩(wěn)性失效模式開展研究，以ITTC A2拖網(wǎng)漁船為例，采用了六自由度數(shù)值預(yù)報模型，選取具有不同航速、浪向、波浪譜截斷區(qū)域和值工況，進行大量密集重復(fù)模擬，再現(xiàn)了船舶發(fā)生騎浪/橫甩現(xiàn)象的過程。通過使用希爾伯特變換計算的瞬時波速來識別增速狀態(tài)，定量研究了船舶發(fā)生增速、橫甩、顯著艏搖和橫搖的概率，為IMO船舶二代完整穩(wěn)性衡準(zhǔn)的應(yīng)用提供技術(shù)支撐。

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

本文中建立了六自由度操縱性運動模型開展不規(guī)則波中騎浪/橫甩的數(shù)值模擬模型，模型中考慮了非線性的回復(fù)力和Froude-krylov（F-K）力以及舵力、螺旋槳動力的作用。

1.1 坐標(biāo)系統(tǒng)

船舶的運動坐標(biāo)系統(tǒng)中包含3個坐標(biāo)系：地球坐標(biāo)系O-xyz、參考坐標(biāo)系-和船體固定坐標(biāo)系O-xyz。各坐標(biāo)系及運動正方向如圖1所示。原點被選為船舶重心。

圖1 坐標(biāo)系與運動方向

位置、速度和力矢量定義如下：

式中：、、、、和是地球坐標(biāo)中的速度；和是變換矩陣。

1.2 6-DOF操縱性運動模型

在六自由度操縱性運動模型中，基于弱耦合理論將橫搖、垂蕩與縱搖MMG模型和耐波性模型在同一坐標(biāo)系下進行耦合，如式（3）所示。

式中：為船舶質(zhì)量，kg；m為附加質(zhì)量，kg；I為慣性矩，m；J為附加慣性矩；(X，Y，N)為舵力，kN；X為縱蕩波浪力，kN；X為阻力X為螺旋槳推力，kN；B是Ikeda半經(jīng)驗公式估計的橫搖粘性阻尼，N·s/m；(X，Y，N)是船體力，kN；z是船體力Y作用點的垂向坐標(biāo)；F()為非線性回復(fù)力，kN；F()為F-K力，kN；F()為輻射；F()為繞射力，kN。各力的計算方法如式（4）所示。

式中：、和是阻力曲線的擬合系數(shù)；Y為橫蕩力，kN；N為搖艏力矩一階水動力導(dǎo)數(shù)；X為縱蕩力；Y為橫蕩力，kN；N為搖艏力矩二階水動力導(dǎo)數(shù)。根據(jù)脈沖響應(yīng)函數(shù)（IRF）方法，將切片（STF）法得到的頻域結(jié)果通過時延函數(shù)R和Q轉(zhuǎn)移到時域，得到時域輻射和繞射力F()，F()。由于騎浪時遭遇頻率很小，因此在此計算中忽略了垂蕩和縱搖運動的繞射力。螺旋槳和舵力的計算方法可參見文獻[15]中所示。

1.3 非線性回復(fù)力和F-K力

本文通過瞬時濕表面上的壓力積分來計算船舶在波浪中所受到的非線性回復(fù)力和F-K力。在計算過程中，船體和上層甲板由幾個非均勻有理B樣條（NURBS）曲面組成，如圖2所示。

圖2 ITTC A2 漁船船體NURBS曲面

每個曲面的面積為A，在船體固定坐標(biāo)系下，其形心坐標(biāo)r=(x，y，z)，曲面的法向量n=(，，)。因此，非線性回復(fù)力和F-K力的計算公式如下：

式中：吃水d(x，)和波幅()的計算公式如下：

式中：A和A代表規(guī)則波和不規(guī)則波中的波幅；S()是不規(guī)則的波譜；，，，ω分別是波數(shù)、浪向角、波頻和波遭遇頻率；N是波浪譜的離散分量數(shù)；θ是隨機相位；上標(biāo)()和()表示船體固定坐標(biāo)-XYZ和地球坐標(biāo)O-XYZ中的 矢量。

2 騎浪橫甩數(shù)值模擬

2.1 船舶參數(shù)和計算工況

用于數(shù)值模擬的船模是 ITTC A2 漁船。船舶和舵的主尺度如表1所示。

表1 ITTC A2 漁船主尺度

在數(shù)值模擬中，采用了JONSWAP波浪譜，其有義波高為H= 2.68 m，譜峰周期T= 5 s，如下頁圖3所示。在仿真過程中，以Δ=0.004 rad/s的間隔對頻率范圍進行離散，這樣可以獲得2π/Δ=1 570 s的不重復(fù)仿真時間。不規(guī)則波產(chǎn)生的參數(shù)如下頁表2 所示。

表2 不規(guī)則波生成參數(shù)

圖3 JONSWAP波浪譜與頻率范圍

在上述不規(guī)則波下，選取了具有不同波浪譜截斷頻率范圍、值、標(biāo)稱弗勞德數(shù)（）和浪向角的工況。對于每個工況，均進行了225次具有不同隨機相位的獨立重復(fù)實現(xiàn)。基于這些密集的重復(fù)模擬結(jié)果，采用直接計數(shù)法來研究不同工況下騎浪、橫搖和顯著橫搖的概率。詳細工況設(shè)置如表3所示。

表3 計算工況

2.2 瞬時波速的計算

在數(shù)值模擬過程中，為了判定不規(guī)則波中船速是否達到波速發(fā)生騎浪，需要計算不規(guī)則波中瞬時波速。本文中基于希爾伯特變換計算瞬時波速，這一方法類似于SPYROU描述的方法，其是使用希爾伯特變換得到瞬時相位，然后計算瞬時相位相對于位置和時間的偏導(dǎo)數(shù)，來獲得瞬時波速。

對于任意數(shù)據(jù)序列()，其希爾伯特變換()可以導(dǎo)出為：

式中：..即principal value的縮寫，是柯西主值積分；希爾伯特變換被定義為()和1/的卷積積分；()與()的相位差為π/2；此外，希爾伯特變換將真實數(shù)據(jù)序列傳輸?shù)綇?fù)數(shù)信號。復(fù)數(shù)信號定義為：

式中：()是瞬時振幅，即包絡(luò)；()是瞬時相位。

因此，瞬時頻率()可以定義為瞬時相位()的導(dǎo)數(shù)，如下所示：

將替換為式（10）、（11）和（12）中的，瞬時波數(shù)()可以計算如下：

最后，瞬時波速計算如下：

在本研究中，希爾伯特變換是使用快速傅里葉變換（FFT）實現(xiàn)的。

2.3 數(shù)值模擬結(jié)果

基于第2章中的數(shù)值模型，模擬了不同工況下船舶在不規(guī)則波中的運動。在每種工況下均進行了225次獨立的重復(fù)實現(xiàn)。部分運動響應(yīng)結(jié)果如圖4～9所示。每個圖包含有4個子圖，從上到下分別為縱搖與垂蕩、橫搖、艏搖，以及船速與波速的時歷曲線。

圖4中，在低情況下，船速達不到波速，幾乎未發(fā)生騎浪。如圖5所示的高下，出現(xiàn)大量航速高于波速的情況。在不規(guī)則波下，航速高于波速的狀態(tài)稱為增速，增速現(xiàn)象與騎浪的發(fā)生密切相關(guān)；同時，在0°浪向角下，橫搖和艏搖運動幾乎 為0。

圖4 Fr=0.3、浪向角為0°工況的模擬結(jié)果

圖5 Fr=0.45、浪向角為0°工況的模擬結(jié)果

圖6～8顯示了浪向角30°下的數(shù)值模擬結(jié)果。當(dāng)浪向角較大時，可以觀察到橫搖和艏搖運動的不穩(wěn)定性。對于某些工況，船舶在經(jīng)歷相對較長的增速后發(fā)生傾覆（如圖8中的=100 s左右所示），這表明船舶傾覆和增速之間可能存在一定的相關(guān)性。在=0.45、浪向角為10°的圖9中，橫搖和艏搖運動的不穩(wěn)定性也很明顯，在經(jīng)歷較長時間增速后，傾覆發(fā)生在=364.74 s左右。

圖6 Fr=0.40、浪向角為30°工況的模擬結(jié)果

圖7 Fr=0.43、浪向角為30°工況的模擬結(jié)果

圖8 Fr=0.43、浪向角為30°工況的模擬結(jié)果，船舶傾覆t =100 s

圖9 Fr=0.45、浪向角為10°工況的模擬結(jié)果，船舶傾覆t =364 s

3 增速和橫甩運動的統(tǒng)計分析

基于第2章中數(shù)值模型開展了大量重復(fù)數(shù)值模擬，本節(jié)中采用直接計數(shù)方法對數(shù)值模擬結(jié)果進行分析，開展不規(guī)則波下增速和橫甩運動的統(tǒng)計評估。首先，根據(jù)THEMELIS等的定義，給出了增速、增速后橫甩、顯著艏搖和顯著橫搖的閾值和概率的估計方法：

（1）增速現(xiàn)象：當(dāng)船速超過瞬時波速時，識別出增速發(fā)生，船速低于標(biāo)稱時，增速結(jié)束。其概率估計為增速的總持續(xù)時間與總模擬時間的比值。

（2）增速后橫甩：當(dāng)船舶的航向角超過指令航向±5°并且船舶正在經(jīng)歷增速時，認(rèn)為其發(fā)生增速后橫甩。其概率估計值定義為增速后橫甩的總持續(xù)時間與增速的總持續(xù)時間的比值。在分析過程中還計算了總模擬時間內(nèi)增速后橫甩的概率估計值，其定義為增速后橫甩的總持續(xù)時間與總持續(xù)時間的比值。

（3）顯著艏搖：當(dāng)船舶艏搖偏航的絕對值超過5°時，發(fā)生顯著艏搖。其概率估計值定義為偏航的總持續(xù)時間與總模擬時間的比值。

（4）顯著橫搖：當(dāng)船舶橫搖角的絕對值超過5°和10°時，發(fā)生顯著橫搖。其概率估計值定義為橫搖超限的總持續(xù)時間與總模擬時間的比值。

3.1 增速現(xiàn)象的概率估計

首先，為了驗證獨立模擬的數(shù)量是否足夠，使用式（16）和（17）計算增速概率的標(biāo)準(zhǔn)差隨重復(fù)實現(xiàn)樣本數(shù)增加的變化情況。

式中：是同一工況下重復(fù)實現(xiàn)樣本的數(shù)量；x是第個實現(xiàn)的概率估計值；μ為平均值；σ為標(biāo)準(zhǔn)差。0.40、0.45和浪向角=20的結(jié)果如圖10和下頁圖11所示。從這兩個圖中可得，200次模擬后概率的標(biāo)準(zhǔn)差是穩(wěn)定的。因此，同一工況下獨立重復(fù)實現(xiàn)的數(shù)量足以獲得具有統(tǒng)計意義的結(jié)果。

圖10 Fr=0.40，浪向角為20°，騎浪集成概率的標(biāo)準(zhǔn)差

圖11 Fr=0.45，浪向角為20°，騎浪集成概率的標(biāo)準(zhǔn)差

圖13～18顯示了不同工況下增速的概率估計結(jié)果，在這些圖中，均采用了箱型圖來表示增速的概率估計結(jié)果，見圖12。

圖12 箱型圖

根據(jù)THEMELIS等，箱型圖主要基于4分位數(shù)繪制。4分位數(shù)定義為將頻率分布劃分為25%、50%、75%和100%4個部分，如圖12所示，概率的中位數(shù)和平均值以粗實心和虛線水平線的形式繪于圖中?；疑珔^(qū)域下限表示第 1個4分位數(shù)的邊界，而灰色區(qū)域的上限表示第 3個4分位數(shù)3的邊界。此外，最高和最低的細實心水平線表示與中位數(shù)有最大偏差的值。

圖13顯示了不同下浪向角0°的增速概率估計結(jié)果。

圖13 浪向角0°時，不同F(xiàn)r的增速概率統(tǒng)計結(jié)果

可以發(fā)現(xiàn)：與規(guī)則波中的騎浪相似，隨著船舶標(biāo)稱的增加，增速的概率也顯著增加。因此，在不規(guī)則波浪中增速現(xiàn)象的概率與騎浪的概率高度相關(guān)。

圖14～16顯示了每個標(biāo)稱在不同浪向角下增速的概率估計結(jié)果。隨著標(biāo)稱從=0.30的圖14增加到=0.45的圖16，所有浪向角的增速概率顯著增加。對于每個標(biāo)稱下，不同浪向角的騎浪概率沒有顯示出很大的差異。在較大的浪向角下，騎浪的可能性往往會降低。

圖14 不同浪向角時，F(xiàn)r=0.3和Fr=0.35時的增速概率統(tǒng)計結(jié)果

圖16 不同浪向角時，F(xiàn)r=0.43和Fr=0.45時的增速概率統(tǒng)計結(jié)果

下頁圖17為不同波浪譜截斷頻率范圍下的增速概率統(tǒng)計結(jié)果，從圖中可以看出，當(dāng)波浪譜截斷頻率范圍增加到0.9時，增速的概率估計上升。

圖17 不同波浪譜截斷頻率范圍下，F(xiàn)r=0.40和Fr=0.45時的增速概率統(tǒng)計結(jié)果

圖15 不同浪向角時，F(xiàn)r=0.38和Fr=0.40時的增速概率統(tǒng)計結(jié)果

3.2 增速后橫甩的概率估計

圖18～20顯示了增速后橫甩的概率估計。從圖中可知，隨著浪向角的增加，增速后橫甩的概率估計明顯增加。因此，可能橫甩的概率與浪向角的大小密切相關(guān)。當(dāng)浪向角為10°時，增速中可能橫甩的總持續(xù)時間隨著的增加而緩慢增加。當(dāng)浪向角在20°時，其值大部分都在10%～25%的范圍內(nèi)。當(dāng)浪向角達到30°時，增速中可能橫甩的概率估計值先增加后降低。對于≥0.38的工況，可能橫甩的總持續(xù)時間可占增速總時間持續(xù)時間的30%～50%。

圖18 不同浪向角時，F(xiàn)r=0.35和Fr=0.38時的增速后橫甩統(tǒng)計結(jié)果

圖19 不同浪向角時，F(xiàn)r=0.40和Fr=0.43時的增速后橫甩統(tǒng)計結(jié)果

圖20 不同浪向角時，F(xiàn)r=0.45的增速后橫甩統(tǒng)計結(jié)果

3.3 顯著艏搖的概率評估

圖21 ～25顯示了不同工況下增速后橫甩和顯著艏搖的概率估計值。在每個圖中，（a）是在總模擬時間內(nèi)增速后橫甩的概率估計，而（b）是顯著艏搖的概率估計。首先，造成顯著艏搖有2個原因：一是波浪力引起的正常偏航，另一個是騎浪引起的艏搖不穩(wěn)定。后者被稱為橫甩。因此，在總模擬時間內(nèi)，所有增速后橫甩事件都包含在艏搖角偏航事件中。也就是說，在總模擬時間（a）中，在增速后橫甩的概率估計總是小于顯著艏搖（b）的概率估計，這可以從圖21～25中觀察到。此外，（b）-（a）產(chǎn)生由波浪力引起的正常偏航的概率估計值。

從圖21～25中發(fā)現(xiàn)，隨著的增加，總時間（a）中增速后橫甩的概率估計顯著增加，而波浪力（b）-（a）誘導(dǎo)的正常偏航運動的概率估計值降低。對于≥0.43，增速后橫甩事件可以占顯著艏搖角事件的50%～80%，如圖24和圖25所示。此外，從這些圖中可以看出，隨浪向角的增加，增速后橫甩和顯著艏搖的概率估計均增加。

圖21 不同浪向角時，F(xiàn)r=0.35的總時間下，增速后橫甩（a）與顯著艏搖（b）統(tǒng)計結(jié)果

圖22 不同浪向角時，F(xiàn)r=0.38的總時間下，增速后橫甩（a）與顯著艏搖（b）統(tǒng)計結(jié)果

圖24 不同浪向角時，F(xiàn)r=0.43的總時間下，增速后橫甩（a）與顯著艏搖（b）統(tǒng)計結(jié)果

圖25 不同浪向角時，F(xiàn)r=0.45的總時間下，增速后橫甩（a）與顯著艏搖（b）統(tǒng)計結(jié)果

圖23 不同浪向角時，F(xiàn)r=0.40的總時間下，增速后橫甩（a）與顯著艏搖（b）統(tǒng)計結(jié)果

3.4 顯著橫搖概率評估

顯著橫搖閾值為10°時不同標(biāo)稱的顯著橫搖概率估計值如下頁圖26和圖27所示。從中可以觀察到，顯著橫搖概率降低至0.03以下，顯著橫搖的概率估計值均隨浪向角增加而增加，這與增速后橫甩的概率變化規(guī)律相同。不同值下的顯著橫搖概率估計值請參見下頁圖28和下頁圖29?？梢?，小值下，顯著橫搖概率估計值會增加。

圖26 不同浪向角時，F(xiàn)r=0.38和Fr=0.40時的顯著橫搖（10°）統(tǒng)計結(jié)果

圖27 不同浪向角時，F(xiàn)r=0.43和Fr=0.45時的顯著橫搖（10°）統(tǒng)計結(jié)果

圖28 不同GM值時，F(xiàn)r=0.38和Fr=0.40時的顯著橫搖（10°）統(tǒng)計結(jié)果

圖29 不同GM值時，F(xiàn)r=0.43和Fr=0.45時的顯著橫搖（10°）統(tǒng)計結(jié)果

4 結(jié) 論

本研究中，基于不規(guī)則波中騎浪橫甩六自由度數(shù)值預(yù)報模型，對不同浪向角、航速、波浪譜截斷范圍和GM值工況下的船舶騎浪橫甩運動進行了密集的數(shù)值模擬。然后，在數(shù)值模擬結(jié)果基礎(chǔ)上應(yīng)用直接計數(shù)法，估計了增速、增速后橫甩、顯著艏搖和顯著橫搖的概率，并得出以下結(jié)論：

（1）使用希爾伯特變換計算的瞬時波速用于識別增速。在不同浪向角下，船舶經(jīng)歷長時間的增速后橫甩傾覆，這表明傾覆和增速之間存在相關(guān)性。

（2）類似于規(guī)則波中的騎浪，增速的概率隨著標(biāo)稱的增加而顯著增加，可見不規(guī)則波中增速與騎浪的發(fā)生高度相關(guān)。同時，波浪譜頻率截斷范圍對增速的概率估計有較大的影響。

（3）增速后橫甩的概率隨著浪向角的增加而顯著增加。當(dāng)浪向角達到30°時，橫甩的總持續(xù)時間可占增速總持續(xù)時間的30%～50%，特別是對于高。

（4）顯著艏搖包括增速后橫甩和波浪力誘導(dǎo)的艏搖2部分；而在高下，在總模擬時間中，增速后橫甩事件可以占顯著艏搖角事件的50%～80%。

（5）顯著橫搖的概率估計值均隨浪向角增加而增加，這與增速后橫甩的概率變化規(guī)律相同。