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

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

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

2 海軍研究院特種勤務研究所，上海200235

0 引 言

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

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

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

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

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

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

1 坐標系及廣義參數

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 廣義參數

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

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

圖2 坐標系轉動與平移Fig.2 Coordinate systems transformation by rotation and translation

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

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

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

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

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

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

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

3 水動力計算式

3.1 瞬時波面下壓力分布

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

水下壓力p的分布式為e 指數形式

式中：k為波數；χ為浪向角；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，這就意味著瞬時波面下的壓力用等效靜水壓力近似，等效靜水的密度ρ′在波峰區域（cosΘ＞0）較實際水的密度ρ小，在波谷區域（cosΘ＜0）較實際水的密度ρ大。

3.2 入射波浪的位移力

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

式中，F和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 部分構成，也即式（2）中的Mζ。由此，由式（7）和式（8）不難確定出垂向力的橫向坐標位置yζ=K/Fz。推導中，利用了頂浪航行條件sinχ=0。

3.3 船體附加質量和阻尼系數估算

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

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

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

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

2）切片垂蕩附加質量和阻尼。

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

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

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

3）船體附加質量和阻尼系數。

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

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

4 參數橫搖數值計算及力學機理分析

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

4.1 參數橫搖運動計算

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

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

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

表1 參數橫搖運動的預報結果Table 1 Prediction results of parametric rolling motion

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

4.2 參數橫搖力學機理分析

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

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

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

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

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

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

4.3 參數橫搖衡準

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

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

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

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

這就是參數橫搖的判別衡準。

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

5 結 論

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

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

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

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