船舶參數(shù)橫搖數(shù)值計算與力學機理分析

朱軍，朱韜，王智宇，夏齊強，黃昆侖，葛義軍

1 海軍工程大學艦船與海洋學院，湖北武漢430033

2 海軍研究院特種勤務(wù)研究所，上海200235

0 引 言

半個多世紀前，研究人員就注意到，船舶在沒有橫向作用力的縱向波浪中會發(fā)生大幅度的橫搖運動，即船舶參數(shù)橫搖現(xiàn)象。Kerwin［1］采用一階諧波形式表達了回復力矩的變化量，運用馬蒂厄方程構(gòu)建了橫搖運動模型；Paulling 等［2］認為，垂蕩和縱搖的非線性耦合至橫向會引發(fā)橫搖運動，其在垂蕩或縱搖為諧振運動模式的假定下，導出了一階諧波回復力矩形式的馬蒂厄方程，經(jīng)模型試驗，測定了靜水中由垂蕩引發(fā)的參數(shù)橫搖運動。

1998 年，巴拿馬型C11 大型集裝箱船因受極端氣象影響而發(fā)生了參數(shù)橫搖運動，導致甲板集裝箱遭受巨大損失，再次引起了人們對船舶參數(shù)橫搖運動的關(guān)注。有關(guān)船舶參數(shù)橫搖運動計算，可以歸納為2 類：一是采用馬蒂厄方程構(gòu)建橫搖運動模型，其中穩(wěn)性參數(shù)采用遭遇頻率的一階諧波模式；二是采用垂蕩、縱搖和橫搖的耦合數(shù)值計算模型，即直接的搖蕩耦合運動數(shù)值模擬。

France 等［3］對C11 大型集裝箱船遭遇的海況進行了調(diào)查，確認了參數(shù)橫搖發(fā)生的海況條件。Shin 等［4-6］針對船舶發(fā)生參數(shù)橫搖運動的現(xiàn)象給出了相同的物理解釋，即波浪通過船體時水線面面積會發(fā)生周期性的改變，并用遭遇頻率的一階諧波模式描述變化的穩(wěn)性參數(shù)，以垂蕩和縱搖的靜態(tài)或準靜態(tài)平衡方式估計穩(wěn)性參數(shù)的改變量，以馬蒂厄標準方程的數(shù)值圖解作為穩(wěn)定性的判據(jù)。Belenky 等［5］采用Spyrou［7］給出的參數(shù)橫搖判別的衡準取決于穩(wěn)性參數(shù)改變量的大小。

在數(shù)值計算方面，F(xiàn)rance 等［3］采用多個耐波性專業(yè)軟件對參數(shù)橫搖運動進行了計算，這些軟件盡管不是完全針對參數(shù)橫搖運動開發(fā)的，但其計算結(jié)果與試驗結(jié)果在趨勢上相同。傅汝德-克雷洛夫波浪力被認為是船舶發(fā)生參數(shù)橫搖運動的關(guān)鍵。Ma 等［8-10］采用瞬時三維壓力分布的濕表面積積分計算波浪力，計算研究了船舶的參數(shù)橫搖運動。魯江和卜淑霞等［11-12］運用單自由度（類似馬蒂厄方程）和三自由度數(shù)學模型模擬了參數(shù)橫搖運動，并在三自由度模型中采用STF 方法處理了波浪力，結(jié)果顯示，與單自由度模型相比，三自由度模型的效果并非更好。

上述研究均未指明參數(shù)橫搖運動的能量來源與機制（力學機理），Shin 和Belenky 等［4-5］只是泛泛地認為穩(wěn)性參數(shù)的變化是橫搖運動的能量來源，這個機理的認識是數(shù)學性的。另外，在國際海事組織（IMO）正在制定的船舶第2 代完整穩(wěn)性標準中，參數(shù)橫搖的薄弱性衡準是采用穩(wěn)性變化量的大小來衡量［6］，然而，不同類型的船舶其計算結(jié)果離散較大，這表明以穩(wěn)性變化量大小作為參數(shù)橫搖衡準只是權(quán)宜辦法［5］。

本文擬采用所提出的搖蕩耦合切片計算方法［13］，克服大幅運動時參數(shù)含義的模糊性，以及會誘發(fā)更多自由度運動耦合的缺陷，數(shù)值計算船舶的參數(shù)橫搖運動，由此分析發(fā)生參數(shù)橫搖運動的力學機理，并基于能量原理提出船舶參數(shù)橫搖的衡準。

1 坐標系及廣義參數(shù)

1.1 坐標系

3 個坐標系分別為慣性坐標系、船體運動坐標系和船體固定坐標系，如圖1 所示。

圖1 坐標系Fig.1 Coordinate system

慣性坐標系E-ξηζ的原點E固定于靜水面，Eξ軸指向船舶前進方向，Eζ軸垂直向上，Eη軸指向船體右舷。

船體運動坐標系M-xyz在初始時刻與慣性坐標系重合，坐標系沿慣性坐標系的Eξ軸方向以速度U做直線運動。

船體固定坐標系B-xByBzB固定于船體，原點位于船體龍骨基線左右對稱面的船舯處，BzB軸線垂直向上，ByB軸線指向船體右舷，BxB軸線指向船艏。

1.2 廣義參數(shù)

繞船體軸線轉(zhuǎn)動定義的常規(guī)橫傾角和縱傾角只適合小角度情況，它存在轉(zhuǎn)動順序問題，即轉(zhuǎn)動順序不同，船體姿態(tài)亦不同。在極端情況，例如當橫傾角和縱傾角均為90°時，先后不同地轉(zhuǎn)動橫傾角和縱傾角，船體姿態(tài)會有很大的差異。為此，本文定義橫傾角繞船體縱向坐標軸轉(zhuǎn)動，縱傾角繞靜水面橫軸線（My 軸）轉(zhuǎn)動，具體定義如下。

假定船舶靜水漂浮時吃水為T，船體繞Mx軸線轉(zhuǎn)動橫搖角φ，再沿Mz軸平移吃水增量ΔT，最后再繞My軸轉(zhuǎn)動縱傾角θ（圖2）。這里，稱ΔT，θ分別為廣義吃水增量和廣義縱傾角。根據(jù)轉(zhuǎn)動關(guān)系，不難得到船體運動坐標系和船體固定坐標系的轉(zhuǎn)換關(guān)系為：

圖2 坐標系轉(zhuǎn)動與平移Fig.2 Coordinate systems transformation by rotation and translation

2 船體搖蕩耦合運動方程組

習慣上，將運動方程構(gòu)建在船體固定坐標系中，其優(yōu)點是繞坐標軸的轉(zhuǎn)動慣量不隨運動變化，其為常數(shù)，缺點是大幅度運動時易誘發(fā)更多自由度的運動耦合。例如，頂浪航行時，垂向波浪力在轉(zhuǎn)換到船體固定坐標系后會有沿船體橫向的分量，這個力會導致船體固定坐標下的橫蕩和搖艏運動。這就意味著原本在靜水面觀測到的3 個自由度運動，在船體固定坐標下會演變?yōu)? 個自由度的運動，這將導致交叉耦合項的水動力出現(xiàn)，使耦合運動的計算變得更加復雜。

因此，本文將搖蕩耦合運動方程構(gòu)建在船體運動坐標的慣性系中。在該慣性坐標系下，隨時間不斷變化的是船體水下形狀，在垂直面內(nèi)，船體重力和入射波浪力作用會產(chǎn)生垂蕩運動，重力和入射波浪力對船舯的力矩會產(chǎn)生縱搖運動。采用普通切片法，將切片的位移、速度和加速度的作用力沿船體縱向積分得到作用于船體的力/力矩，從而不難得到垂蕩和縱搖運動方程式為：

式中：ZΔT?，ZΔT?，Zθ?，Zθ?，ZUθ?，ZUθ為船體垂蕩的水動力系數(shù)；Mθ?，Mθ?，MΔT?，MΔT?，MUθ?，MUθ為船體縱搖的水動力系數(shù)；Fζ和Mζ，F(xiàn)ζ?和Mζ?，F(xiàn)ζ?和Mζ?分別為入射波浪相對船體位移、速度和加速度的垂向力及縱搖力矩；m為船體質(zhì)量；Jθ為船體縱向轉(zhuǎn)動慣量；xG為船體重心縱向坐標位置；ρg?為船體在靜水中的排水量，其中ρ為水的密度，g為重力加速度，?為靜水中的排水體積。

橫搖運動方程將忽略入射波浪對船體速度和加速度產(chǎn)生的偏心作用，只考慮入射波浪力的位移作用。因此，繞船體重心的橫搖運動方程式可寫為

式中：ΔIφ和Kφ?分別為橫搖附加慣性矩和阻尼系數(shù)；Iφ為船體橫向轉(zhuǎn)動慣量；yζ為作用力的橫向坐標位置；yG船體重心橫向坐標位置。

式（2）和式（3）即構(gòu)成了船體垂蕩、縱搖和橫搖的耦合運動方程組。值得注意的是，F(xiàn)ζ是將船體姿態(tài)（ΔT，θ，φ）和瞬時波浪位置ζ作為一個整體來計算的，方程組通過力Fζ=f(ΔT，θ，φ，ζ)耦合了船體的垂蕩、縱搖和橫搖運動，具體計算公式將在后文給出。

3 水動力計算式

3.1 瞬時波面下壓力分布

由勢流理論，可知船體運動坐標系下小振幅規(guī)則波升ζ可表達為

水下壓力p的分布式為e 指數(shù)形式

式中：k為波數(shù)；χ為浪向角；t 為時間；A 為波幅；ωe為遭遇頻率；Θ為計入了時間和位置的波浪相位。為方便起見，令Θ=ωet+k(xcosχ+ysinχ) 。式（4）并不適合直接做壓力積分計算，為此，通常采用靜水壓力和拉伸修正的近似方法［14］，或不做修正。本文采用線性近似，為此，做一個變量代換：

z=z1+AcosΘ

該代換的含義是，z坐標是以靜水面為零點，而z1坐標則是以瞬時波面為零點。將其代入式（3），忽略高階項后，可得壓力的近似表達式為

p=-ρ′gz1（5）

式中，ρ′=ρ(1-α0cosΘ) ，為等效水的密度，其中α0=Ak，為最大波傾角。式（4）滿足瞬時波面（z1=0）處壓力p=0，這就意味著瞬時波面下的壓力用等效靜水壓力近似，等效靜水的密度ρ′在波峰區(qū)域（cosΘ＞0）較實際水的密度ρ小，在波谷區(qū)域（cosΘ＜0）較實際水的密度ρ大。

3.2 入射波浪的位移力

利用高斯定理，將瞬時波面下船體濕表面的壓力積分轉(zhuǎn)換成體積分：

式中，F(xiàn)和M分別為作用于船體的力和繞坐標原點的力矩；s為船體濕表面的面積；V為波浪下船體濕表面的體積；n為船體表面的單位外法向矢量；r為微分濕表面ds到坐標原點的位置矢量；i，j，k分別為船體運動坐標系的3 個單位坐標矢量。將壓力的表達式（4）代入到上述力的計算公式中，不難導出瞬時波面下船體入射波浪3 個方向的作用力為：

式中：L 為船長；Sx為瞬時波面下船體切片的面積（圖3 中陰影部分面積）。式（7）表明，瞬時波面下船體入射波浪的作用力只有垂向力和縱向力，橫向力為0。顯然，縱向力Fx與波傾角α0成正比，因波傾角α0為小量，因此縱向力Fx是α0級小量，波浪的主要作用力是垂向力Fz。

圖3 瞬時波面下船體切片示意圖Fig.3 Diagram of hull strip under transient wave surface

同樣，可以導出繞坐標軸的橫搖力矩K 和縱搖力矩M：

式中，xs，ys和zs分別為瞬時波面下船體切片面積中心在船體運動坐標系的坐標值。其中，縱搖力矩M由縱向力和垂向力2 部分構(gòu)成，也即式（2）中的Mζ。由此，由式（7）和式（8）不難確定出垂向力的橫向坐標位置yζ=K/Fz。推導中，利用了頂浪航行條件sinχ=0。

3.3 船體附加質(zhì)量和阻尼系數(shù)估算

式（2）需要計算船體的附加質(zhì)量和阻尼系數(shù)等，本文則將按經(jīng)驗公式來估算。本節(jié)中的經(jīng)驗公式和近似圖譜來自于文獻［15］。

1）船體橫搖附加慣性矩和阻尼系數(shù)。

船體橫搖的慣性矩、附加慣性矩和阻尼系數(shù)按照常規(guī)單自由度運動方式，以船體平均位置（靜水平衡狀態(tài)）近似估算。橫搖附加慣性矩和阻尼系數(shù)由經(jīng)驗公式估算：

式中：ρφ為橫搖慣性半徑，ρφ=CB，其中B為船寬，經(jīng)驗系數(shù)C的取值范圍為0.33～0.45；2μφ為無量綱橫搖衰減系數(shù)，其取值范圍為0.11～0.14；GM0為初穩(wěn)性高；Fn 為傅汝德數(shù)。

2）切片垂蕩附加質(zhì)量和阻尼。

用二維浮體的寬度、寬吃水比和面積系數(shù)這3 個參數(shù)來估算船體切片的垂向附加質(zhì)量和阻尼系數(shù)，沿船體縱向積分得出船體的垂蕩附加質(zhì)量、縱搖附加慣性矩和對應(yīng)的阻尼系數(shù)。

二維浮體單位長度的附加質(zhì)量Δm和阻尼系數(shù)ΔNμ估算公式為

式中：Bn為二維浮體水面處的寬度；C1為附加質(zhì)量系數(shù)；Aˉ為衰減系數(shù)，是輻射波幅值與垂蕩幅值之比。附加質(zhì)量系數(shù)C1和衰減系數(shù)Aˉ的大小與遭遇頻率ωe、寬吃水比和水下面積系數(shù)有關(guān)，本文采用文獻［15］的圖譜插值近似。

3）船體附加質(zhì)量和阻尼系數(shù)。

對式（2）等號左邊的船體垂蕩水動力系數(shù)和縱搖水動力系數(shù)，按照切片在瞬時波面下的船體姿態(tài)（ΔT，θ，φ），由式（10）來估算切片附加質(zhì)量Δm和阻尼系數(shù)ΔNμ，然后再沿船體縱向積分得到船體垂蕩和縱搖的附加質(zhì)量及附加慣性矩等。

入射波浪相對于船體速度和加速度的垂向力及縱搖力矩也可按照同樣的方法積分計算。

4 參數(shù)橫搖數(shù)值計算及力學機理分析

本文采用上述動力學模型和基于圖形面域技術(shù)開發(fā)的軟件計算了一艘船舶在頂浪規(guī)則波中的搖蕩運動，以此分析船舶發(fā)生參數(shù)橫搖運動的力學機理。

4.1 參數(shù)橫搖運動計算

計算船舶的船長L=142 m，航速U=14 kn，波陡范圍Hw/λ=0.015～0.050。

圖4 所示為波陡Hw/λ=0.025 時的垂蕩、縱搖和橫搖計算時歷曲線。計算結(jié)果顯示：在該狀態(tài)下呈現(xiàn)出了參數(shù)橫搖運動，橫搖幅值不斷增加；同時還表明，橫搖頻率為縱搖頻率的一半，且位于橫搖共振頻率范圍內(nèi)，這與發(fā)生參數(shù)橫搖的頻率條件是吻合的。

圖4 頂浪中的搖蕩運動計算時歷曲線（Hw/λ=0.025）Fig.4 The calculated time history curves of oscillating motions in head wave（Hw/λ=0.025）

表1 參數(shù)橫搖運動的預(yù)報結(jié)果Table 1 Prediction results of parametric rolling motion

表1 所示為不同波陡和波長條件下的參數(shù)橫搖運動預(yù)報結(jié)果。表中：符號“●”表示橫搖運動振蕩增加，呈現(xiàn)參數(shù)橫搖運動；符號“×”表示橫搖運動衰減，無參數(shù)橫搖運動。預(yù)報結(jié)果顯示，除較小的波陡不發(fā)生參數(shù)橫搖運動外，隨著波陡的增加，發(fā)生參數(shù)橫搖的波段向高頻區(qū)偏移，整體上波長/船長（λ/L）在0.9～1.4 范圍。

4.2 參數(shù)橫搖力學機理分析

將式（3）第3 項的回復力矩寫成通常的形式：

式中，參數(shù)Kφ為回復力矩系數(shù)，相當于船舶橫搖系統(tǒng)剛度，其作用是在每半個橫搖周期內(nèi)，在橫搖角增大的過程中吸收能量，在橫搖角減小的過程釋放所吸收的能量。

圖5 回復力矩系數(shù)和橫搖角計算曲線（λ/L=1.2）Fig.5 calculated curves of recovery moment coefficient and rolling angle（λ/L=1.2）

圖6 回復力矩系數(shù)和橫搖角計算曲線（λ/L=1.5）Fig.6 calculated curve of recovery moment coefficient and rolling angle（λ/L=1.5）

圖5 和圖6 所示為Hw/λ=0.025 時回復力矩系數(shù)Kφ與橫搖角φ隨時間變化的計算曲線（其中Kφ的單位為N·m，φ的單位為（°），圖5、圖6 的縱坐標為它們各自乘以不同系數(shù)后的尺度）。圖5呈現(xiàn)的是參數(shù)橫搖運動，其λ L=1.2；而圖6 呈現(xiàn)的則是橫搖運動振蕩衰減，其λ L=1.5。從圖5 所示的計算曲線可以看出，在橫搖角增加過程中，參數(shù)Kφ低于平均值，處于較低水平；而在橫搖角減小的過程中，參數(shù)Kφ高于平均值，處于較高水平。這就意味著在每半個橫搖周期內(nèi)，在橫搖角增大過程中所吸收的能量小于橫搖角減小過程中釋放的能量；圖6 的情況與之相反，因而不會出現(xiàn)參數(shù)橫搖運動。

因此，參數(shù)橫搖的力學機理是：回復力矩系數(shù)Kφ隨著波浪的作用呈現(xiàn)出圍繞均值的波動，在橫搖角減小過程中，其釋放的能量要大于橫搖角增大過程中回復力矩系數(shù)吸收的能量。該能量差值也即參數(shù)橫搖運動持續(xù)放大的能量來源。

4.3 參數(shù)橫搖衡準

圖5 和圖6 顯示，在不同波長/船長、相同波陡情況下，回復力矩系數(shù)Kφ的均值與振蕩幅值沒有明顯差異。如上述分析，是否發(fā)生參數(shù)橫搖運動取決于系統(tǒng)吸收和釋放能量的平衡與否，系數(shù)Kφ與橫搖角φ隨時間變化的相對位置，即相位關(guān)系決定了吸收與釋放能量的相對大小。現(xiàn)用遭遇頻率ωe的諧波函數(shù)擬合系數(shù)Kφ：

任意一個橫搖半周期內(nèi)，假定橫搖角為正弦函數(shù)，其頻率為用遭遇頻率ωe的一半，即

式中，φ0為半周期內(nèi)的橫搖角幅度，為常數(shù)。因此不難得到，橫搖角半周期內(nèi)回復力矩系數(shù)Kφ吸收和釋放的能量總和為

式（13）表明，只有系數(shù)Kφ的一階諧波分量的積分不為0，常數(shù)項和高階諧波項的積分均等于0。發(fā)生參數(shù)橫搖的條件是能量J＜0。因此，滿足此條件的相位角ε1可由式（12）得到

這就是參數(shù)橫搖的判別衡準。

如果將式（11）只保留常數(shù)項和一階諧波項，則與France 等［3-5］研究的參數(shù)橫搖所采用的馬蒂厄方程相類似。但是，判別參數(shù)橫搖的衡準不同，是否發(fā)生參數(shù)橫搖不取決于一階諧波項系數(shù)的大小，而是一階諧波項的相位角ε1。相位角ε1的變化取決于船波相對運動，其與波浪中船舶的垂蕩、縱搖和由入射波浪引起的水線面形狀變化以及水下排水體積變化相關(guān)。

5 結(jié) 論

本文建立了慣性坐標下的垂蕩和縱搖耦合運動方程，以及船體坐標下的橫搖運動方程，由此構(gòu)成垂蕩、縱搖和橫搖的混合動力學模型，并采用普通切片方法計算了入射波浪力、附加質(zhì)量和阻尼系數(shù)。數(shù)值計算了一艘船舶的垂蕩、縱搖和橫搖耦合運動，預(yù)報了參數(shù)橫搖運動，通過對數(shù)值計算結(jié)果的分析，得出了參數(shù)橫搖運動的力學機理和參數(shù)橫搖判別衡準，研究認為：

1）參數(shù)橫搖的力學機理是：由于回復力矩系數(shù)Kφ是隨時間變化的，導致在橫搖過程中回復力矩系數(shù)吸收和釋放的能量不等，釋放能量大于吸收能量是發(fā)生參數(shù)橫搖的力學機理，也是發(fā)生參數(shù)橫搖后運動持續(xù)放大的能量來源。

2）判別參數(shù)橫搖的衡準是，回復力矩系數(shù)Kφ的一階諧波分量具有超前的相位角ε1?[0，π]。

船舶參數(shù)橫搖力學機理的明確有助于對參數(shù)橫搖運動的認識，而如何建立相位角ε1與船體參數(shù)的關(guān)系尚需深入研究。