规则波浪中舰船摇荡耦合切片计算方法

朱军，夏齐强，王智宇，朱韬，黄昆仑，葛义军

1海军工程大学舰船与海洋学院，湖北武汉430033

2海军研究院，上海200235

0 引 言

切片理论方法首次发表于20世纪50年代［1］，之后衍生出了各种不同的切片方法，并被应用到某些船舶上用于确定速度的效应。Salvesen等［2］提出了线性切片理论方法，该方法可用于描述满足交叉耦合项的对称定理。在线性和摇荡谐振假设下，线性切片势流理论方法满足物面和自由表面的边界条件，尽管该方法采用势流理论计算了水动力并提出了六自由度的计算模型，但由于物面边界条件是在平均位置上得到满足的，因此垂直面和水平面实际上是分离的。北约耐波性一般准则和共同程序专家组［3］的研究表明，采用各种计算机代码计算出的驱逐舰垂直平面的运动结果非常一致。Ogilvie和 Tuck［4］的研究表明，当波长在船宽量级的高频范围内时，切片方法是严格有效的。实践经验表明，对于船长与波长等量级的情况，切片方法能给出良好的结果。Timman等［5］采用摄动展开方法，将原本在瞬时物面上满足的边界条件变换成了在物面平均位置上满足的边界条件，这虽然给求解带来了极大的方便，但却使得垂荡、纵摇与横摇实际上在非线性条件下被解耦了。

为了计算甲板上浪和艏底砰击等非线性波浪载荷问题，需采用线性切片理论方法预报船体的摇荡运动。Kaplan等［6-8］］将线性切片方法预报的运动响应作为输入，在频域中有效解决了运动问题，而载荷预报则在时域中执行。Cong等［9］采用脉冲响应函数方法计算了辐射波浪力和衍射波浪力，同时修正了入射波浪压力分布的Smith效应。Gu等［10］采用线性切片理论方法预报了垂直面的摇荡（特指垂荡和纵摇）运动，并直接在时域中计算了运动、甲板上浪与砰击载荷。林超友和朱军等［11-12］为了预报波浪中的操纵运动，采用势流理论方法计算了波浪力。

线性切片理论方法本质上没有耦合船体的垂直面（垂荡和纵摇）和水平面（横摇）运动，这是因为在计算摇荡速度势和入射衍射波浪速度势时，船体切片只取了其静水的平均位置，导致垂荡、纵摇与横摇的相互影响被忽略。尽管速度势摄动方法为非线性方法，但因其船体切片只取了静水平均位置，故在本质上并没有改善垂荡、纵摇与横摇的耦合性。

为了能计入船体垂荡、纵摇与横摇的运动耦合效应，本文拟基于广义纵倾角和广义吃水增量参数［13］，在船体瞬时运动状态（ΔT，θ和φ）和瞬时波浪状态下提取船体切片，以满足瞬时波面压力为零的条件修正压力分布的计算公式，采用文献［13］所应用的图形面域技术计算切片的静水力与波浪扰动力（傅汝德—克雷洛夫力）之和，采用经验方法估算瞬时状态下的船体切片惯性和阻尼水动力，并由所建立的摇荡耦合运动方程预报规则波浪中的船体摇荡运动。

1 坐标系及波面方程

1.1 坐标系及广义参数

坐标系如图1所示，图中χ为浪向角，t为时间。惯性坐标系为E-ξηζ，坐标平面E-ξη与静水面重合；船体运动坐标系为O-xyz，坐标平面O-xy与静水面重合，沿惯性坐标纵向轴以速度U等速运动；船体固定坐标系O′-xO′yO′zO′与船体运动坐标保持等速前进运动。

将船体固定坐标系绕Ox轴转动φ角，然后沿Oz轴平移ΔT（吃水增量），再绕Oy轴（静水面轴线）旋转θ角，运动坐标系（x，y，z）和固定坐标系（xO′，yO′，zO′）的转换关系（图2）为：

值得注意的是，这里的吃水T和横倾角φ与通常的定义是一致的，但吃水增量ΔT为垂直于静水面的变化量，纵倾角θ为绕船舯剖面与静水面交线的转动角度，文献［13］分别称其为广义吃水增量ΔT和广义纵倾角θ。

1.2 波面方程

令波长为λ，波高为Hw，浪向角为χ，惯性坐标系下的波面方程为：

将式（2）代入到式（1）的第3项中，得到船体固定坐标系下的波面方程为

这里的x，y由式（1）的第1项和第2项确定。显然，这是一个隐式方程，方程右端前4项是静水面方程在船体坐标系下的表达式，第5项是波面相对静水面的升高量。

2 船体摇荡耦合动力学方程及水动力

2.1 船体摇荡耦合运动方程

令瞬时波面下船体静水力与入射波浪力之和的水动力为FFK+B，即将静水力与傅汝德—克雷洛夫波浪力作为一个整体计算，不难写出船体横摇的二阶运动方程式为

在船体纵向平面内，仅考虑纵摇与垂荡运动的相互耦合项，按类似方法可得到简化后的垂荡与纵摇耦合运动方程式为

式中，P为船体重力。水动力FFK+B作用点（xFK+B,yFK+B,zFK+B）和船体重心（xG,yG,zG）系指船体运动坐标，由式（1）确定与船体坐标系下的转换关系。

式（4）和式（5）构成了船体垂荡与纵摇、横摇耦合的摇荡运动方程。

2.2 瞬时波面下船体静水力与波浪扰动力之和

普通线性切片方法是将摇荡复原力与波浪扰动力分开处理，且船体切片只取平均水线位置，这导致波峰和波谷通过切片时无差异化。本文基于切片处的瞬时波面方程（3），根据切片瞬时波面下压力分布积分计算水动力，并将静水力与波浪扰动力之和作为整体处理。由势流理论可知，小振幅规则波的压力分布为：

式中，ρ为水的质量密度。该压力分布式忽略了速度平方的二阶小量，瞬时波表面压力等于零的条件并不满足，波峰处的压力估计偏大，波谷处的压力估计偏小［14］，即史密斯效应。为此，令

该转换就是将波面作为压力计算的起点。将瞬时波面下压力分布改写为

式中，A=Hw/2。式（9）满足瞬时波面（z1=0时）压力p=0，同时还满足z1→-∞时p→-ρgz。在t时刻x处的船体切片，对式（9）的指数项进行线性化处理并用y=0近似，则近似的压力分布为

于是，由式（10）对瞬时波面下的切片进行积分，然后再沿船体纵向积分，便可得到船体水动力FFk+B。

显然，式（10）的压力修正体现了波峰和波谷位于切片时的差异，同时，还包含了船体运动姿态（ΔT，θ和φ）的效应，即垂荡、纵摇与横摇的耦合效应。

2.3 船体水动力系数

式（4）和式（5）中的附加惯性水动力系数及阻尼水动力系数，将由所估算的切片水动力系数沿船体纵向积分计算得到。

1）船体横摇附加惯性矩和阻尼系数。

横摇附加惯性水动力和阻尼系数取切片静水平均位置，按通常的经验公式估算：

式中：ρφ=CB，为横摇惯性半径，其中B为船宽，经验系数C取0.40；h为初稳性高；Fn为傅汝德数。船体横摇运动的非线性效应明显，Himeno［15］针对横摇阻尼预测进行了综述，更加有针对性的估算公式可参看文献［15］。

2）船体垂荡、纵摇附加惯性矩和阻尼系数。

式（5）中的船体水动力系数由瞬时波面下切片的水动力系数积分得到，瞬时波面下切片的垂荡附加质量和阻尼系数Δb估算为：

式中，系数C1和由遭遇频率ωe、切片宽—吃水比Bn/dn及 切 片 面 积An的 系 数βn（βn=An/(Bndn)）按二维浮体试验图谱计算［16］。

值得注意的是，式（12）所估算的是瞬时波面下的船体切片，此时，切片状态与船体纵倾角θ、横倾角φ、吃水增量ΔT以及波浪瞬时相位有关，故估算所得到的船体水动力系数是时变的，同时也包含了垂荡、纵摇与横摇的耦合影响。

由式（12）沿船体纵向积分，不难得到式（5）等号左端的惯性和阻尼水动力系数。

2.4 船体摇荡耦合运动方程特征

式（4）和式（5）构成的非线性时变耦合的船体摇荡耦合运动方程组具有时变和耦合的非线性特征。

1）水动力系数的时变性。

式（5）中的垂荡和纵摇的附加惯性矩、阻尼系数及附加惯性静矩共6个水动力系数均是针对瞬时波面下切片经积分而得到，瞬时波面下切片的大小和形状随船体运动状态（ΔT，θ和φ）和波浪相位的不同而变化。所以，式（5）具有系数时变的特征，而水动力系数的时变性则意味着运动方程包含了船体运动状态和波浪相位的效应。

2）垂荡、纵摇和横摇的运动耦合性。

水动力FFK+B为静水力与波浪扰动力（傅汝德—克雷洛夫波浪力）之和，其值取决于瞬时船体运动状态（ΔT，θ，φ）和波浪相位。水动力FFK+B既是构成船体横摇的激励力，也是垂荡、纵摇运动的激励力，可见水动力FFK+B耦合了船体垂荡、纵摇和横摇的运动。

3 数值计算及分析

3.1 船型主尺度及波浪参数

计算舰船的水线长L=142 m，静水吃水T=4.535 m，航速为0。计算的波浪参数为：

1）浪向角χ=0°，45°，90°，135°，180°，即顶浪、首斜浪、正横浪、尾斜浪和顺浪的典型浪向角。

2）波长船长之比λ/L=0.2～5.0，波浪频率从低频到高频，涵盖共振频率。

3）波陡Hw/λ=0.01，0.04，对应的最大波倾角α0=1.8°，7.2°。

3.2 数值计算及分析

数值计算初始运动参数：T=4.535 m，ΔT0=0.0 m，θ0=0.0°，φ0=-2.0°；其他运动参数的值取为0。其中T是静水平衡状态值，为一常数，φ0的初值不为0，相当于一个初始的外部扰动作用。

图3～图7所示为所计算的各个浪向角的摇荡幅频曲线，其计算结果的变化规律与理论分析一致。图3～图5的结果显示，垂荡和纵摇对横摇的耦合呈现出明显的非线性效应，图8和图9则表明，在横摇共振区内耦合还呈现出偏摇现象。

1）横摇幅频曲线计算分析。

由图3～图5中的横摇幅频曲线可以看出，波陡较高（Hw/λ=0.04）的共振区的幅频曲线明显低于波陡较低（Hw/λ=0.01）共振区的幅频曲线，这表明了横摇共振区幅频曲线对波陡的非线性，即波高对横摇幅频曲线具有非线性效应。由于计算中的阻尼系数按式（11）取常数值，因此可以认为该效应是因垂荡、纵摇耦合所导致，即垂荡、纵摇对横摇运动的耦合非线性效应。

图4和图5显示，高频区（ωe≈1.0～1.2）的横摇幅频曲线有一个较小的峰值，该频率范围大致对应垂荡、纵摇的共振频率，因此高频小峰值表明了横摇对垂荡和纵摇耦合的响应。

2）垂荡、纵摇幅频曲线计算分析。

图6和图7所示分别为顶浪和顺浪时的垂荡、纵摇幅频曲线。图4～图7所示的垂荡、纵摇幅频曲线的计算结果显示，不同波陡的垂荡、纵摇幅频曲线在所有频率下基本重合，这说明波高对垂荡、纵摇幅频曲线不具有非线性效应。但在正横浪时，在横摇共振区内，不同波陡情况下的垂荡、纵摇幅频曲线具有一定程度的差异，该差异被认为是横摇对垂荡、纵摇运动耦合的非线性效应。

图4～图7所示的纵摇幅频曲线具有相同的响应规律，即在共振频率（ωe≈1.0～1.2）范围内有一个小的峰值响应。斜浪（首、尾斜浪）的纵摇幅频曲线与纵向波浪（顶浪和顺浪）的纵摇幅频曲线具有相似的响应规律，表明横摇对纵摇不存在耦合效应。

垂荡幅频曲线有明显的差异说明横摇对垂荡的耦合效应显著，该耦合效应出现在频率约大于横摇共振频率之后。纵向波浪中的稳态横摇趋近于0，也即图6和图7所示的垂荡幅频曲线不包含横摇的效应，对比图4和图5斜浪中响应曲线所呈现出的差异，说明横摇对垂荡具有抑制的耦合作用。

3）横摇偏摇现象。

图8所示为波陡Hw/λ=0.04时的正横浪摇荡运动时历曲线。横摇角时历曲线存在着明显的左、右摇荡幅值间的偏差，即在波浪传播方向，摇荡幅值大于其反向幅值，本文称其为横摇偏摇现象。横摇偏摇现象也是垂荡、纵摇耦合的一种效应，计算表明，当波高较大且在横摇共振区内时，会呈现出显著的横摇偏摇。图9所示为Hw/λ=0.04时的横摇偏摇幅值φa0幅频曲线，其中正横浪的横摇偏摇幅值最大。

4）摇荡耦合的非线性效应机理初步分析。

上述摇荡幅频曲线计算分析表明，垂荡、纵摇与横摇的耦合会导致运动响应的非线性效应，其中，最为明显的是垂荡、纵摇对横摇幅频曲线的非线性效应，以及横摇偏摇现象。由于船体线型具有上、下和前、后不对称的特点，所以摇荡耦合的非线性效应即来自于该船体线型的不对称性。

就切片而言，波陡越大，平均水线以上和水线以下线型的不对称性也越大，船体平均水线上、下线型的不对称性是共振区横摇幅频曲线非线性的主要根源。

船体前、后线型的不对称性是产生横摇偏摇现象的原因。图8所示为垂荡（ΔT）、纵摇（θ）和横摇（φ）的时历计算曲线。由纵摇（θ）时历曲线图可以看出，纵摇运动呈现不均匀的摇荡，埋艏（θ＜0）滞留时间较长，抬艏（θ＞0）滞留时间较短，此现象源自船体艏、艉线型的不对称性，符合预期。对比横摇与纵摇时历曲线表面发现，二者基本是反向同步，纵摇的不均匀性传递给横摇就会导致其左、右横摇不对称运动的出现，从而引起偏摇现象——左、右横摇幅值不等。可以预见，前、后线型对称的几何船体将不会出现横摇的偏摇现象。

4 结 论

本文基于广义纵倾角θ和广义吃水增量ΔT，以及船体坐标系下的波面方程，采用切片理论方法建立了垂荡、纵摇与横摇耦合的动力学方程组；基于瞬时波面方程，为满足波面压力处为0的条件，合理修正了瞬时波面下压力分布的计算公式，将静水力与波浪扰动力之和FFK+B作为整体积分予以了计算。规则波浪中，舰船耦合摇荡运动响应数值计算结果表明：

1）采用耦合摇荡切片计算方法可预测出垂荡、纵摇对横摇的耦合非线性效应；

2）采用耦合摇荡切片计算方法可预测出共振区内的横摇偏摇现象；

3）船体线型的上、下不对称性是引起共振区横摇非线性的主要原因；

4）船体线型的前、后不对称性是引起共振区横摇偏摇现象的主要原因。

本文提出了实用性的耦合摇荡运动计算方法，开发了相应的计算软件，较好地计入了垂荡、纵摇与横摇之间的相互耦合效应，对高海况下耐波性预报有积极的作用。后续，还将深入研究并完善计算方法。