0 引　言

随着航速的提高，船舶的航行姿态因流体垂向支持力的比重与作用位置不同而发生改变。根据作用于船体的支持力和船舶首尾吃水，将船舶航态分为3种典型的状态：排水航行状态、过渡状态与滑行状态。排水型船舶以静浮力为垂向的主要支撑力，而高速滑行艇重量主要依靠水动升力支撑^[1]。高速滑行艇在航行时一般有排水状态和滑行状态2种状态。高速滑行艇在高速航行过程中受到水动升力的影响，船舶重量与静水浮力不再稳定平衡，仅有部分船体与水接触，船体湿表面相比于静浮时的结果有较大差别。处于滑行状态的滑行艇在波浪中运动的计算一直是滑行艇水动力特性研究的重点和难点。常规排水型船舶水动力的预报方法未能充分考虑流体动升力的影响，使其不再适用于高速滑行艇。基于势流理论针对滑行艇在波浪中时域运动的研究较为匮乏，尚未形成完整的研究理论与计算方法，因此需要建立高速滑行艇的时域运动计算模型。

滑行艇在高速滑行过程中受到水动升力的作用，艇体抬出水面，船体的吃水与排水体积发生显著变化，使得运动具有较强的非线性特征。此外，高速滑行艇的运动响应与波幅的非线性关系也限制了频谱分析方法预报不规则波中船舶运动的应用。针对高速滑行艇在波浪中的水动力计算分析方面，迄今为止尚无完善的研究理论与计算方法。围绕高速滑行艇耐波性研究的数值方法主要有2种：其一是基于理想流体假设的势流方法，包括排水型船舶运动计算的切片法结合湿表面积变化修正方法^[2]、基于2D+T理论结合二维剖面入水砰击模型建立非线性运动模型^[3]等；其二是基于粘流理论的CFD方法。

针对排水型船舶运动计算的切片法结合湿表面积变化修正方法，别所正利等^[2]基于切片法计算滑行艇的运动，数值结果表明切片法在低速排水状态时（Fn<0.5）预报结果与试验结果一致。MARTIN^[4]应用切片法的思想，结合波倾角引起的流体动力增量、与波浪轨圆运动相关的力、波浪引起的静水压力增量以及浸湿长度变化引起的力增量求解波浪干扰力，建立了求解滑行艇在规则波中线性运动的模型。KATAYAMA等^[5]基于模型试验，采用线性和非线性切片法对滑行艇静水直航与强迫振荡运动进行了计算研究，数值结果与试验吻合较好。董文才等^[6]在滑行艇纵向运动研究中，引入查洁法以考虑滑行动升力影响，基于切片法提出了同时考虑浮态变化和滑行升力影响的滑航法，结果表明滑航法可较为准确地预报棱柱型滑行艇在中高速阶段的纵向运动。

囿于切片法的局限性，有学者采用2D+T理论结合二维剖面入水砰击模型建立非线性运动模型，ZARNICK^[3]基于细长体假定，同时结合楔形体附加质量理论，推导了滑行面在波浪中所受到的流体动力，并建立了滑行艇在规则波中纵荡、垂荡和纵摇模态的非线性耦合运动模型。基于此方法预报Fridsma滑行艇高航速时的运动结果与试验值吻合良好。戴仰山等^[7]基于Zarnick非线性运动模型，分别采用时域方法和谱分析方法预报了滑行艇在不规则波中的运动响应。任慧龙等^[8]基于非线性运动模型结合辐射能量法探讨了滑行艇在规则波中的波浪增阻问题，并推导了滑行艇在不规则波中的波阻增加的计算方法。KEUNING^[9]基于模型试验，引入了经验公式考虑滑行艇升沉与纵倾的影响，同时拓宽了Zarnick模型对变斜升角滑行艇的应用。HICKS等^[10]、GARME^[11]应用Zarnick模型对滑行艇在波浪中的运动预报做了大量研究与改进。SUN等^[12]基于非线性边界元法求解二维物体入水砰击问题，计算中考虑了重力对喷溅与兴波的影响，结合2D+T理论研究了滑行艇的静水直航水动力性能与波浪中的运动问题。GHADIMI等^[13]在Zarnick模型基础上，结合非对称楔形体入水砰击模型建立了滑行艇在规则波中的六自由度时域运动计算方法。GHADIMI等^[14]进一步计算分析了滑行艇在恒定横倾角下的迎浪运动。PENNINO^[15]在无小纵倾角的假设下考虑二阶项的影响，详细推导了滑行艇非线性运动方程，并计算分析了滑行艇的时域运动，相比原始Zarnick方法的结果更为准确。TAVAKOLI等^[16]使用非对称楔形体入水砰击的速度势，结合2D+T理论建立了垂荡、纵摇与横摇耦合的时域运动模型，并探讨了有无横摇对附加质量与运动的影响。HOSSEINZADEH等^[17]采用强迫振荡试验给出了滑行艇附加质量系数的拟合表达式，以此改进了Zarnick方法在滑行艇非线性运动方面的预报精度。TAVAKOLI等^[18]基于模型试验、CFD方法与2D+T方法对比研究了滑行艇在波浪中的运动响应，并给出了加速度、波浪增阻与艇底压力分布的结果，证实了2D+T方法的准确可靠。BILANDI等^[19]将该方法拓展至双断阶滑行艇运动的预报中。

综上所述，目前研究高速滑行艇在波浪中运动的势流方法主要有2种：一是排水型船舶运动计算的切片法结合湿表面积变化修正方法；二是应用2D+T理论结合二维剖面入水砰击模型建立非线性运动模型。修正的切片法在滑行艇高速航行阶段的运动预报结果与试验值存在较大偏差。基于2D+T模型的势流方法计算滑行艇运动效率较高，同时可以保证数值预报的精度。但在楔形体附加质量的求解过程中，相关参数的应用方面并未统一。

本文针对高速滑行艇，基于2D+T理论结合对称楔形体入水砰击模型，给出修正后的附加质量系数表达式，详细推导并建立滑行艇在迎浪中的非线性时域运动方程，自主开发数值计算程序，并以Fridsma滑行艇作为验证算例，开展高速滑行艇在迎浪规则波中的运动响应计算研究。

1 滑行艇运动方程

常规低速排水型船舶的重量与静浮力平衡，而高速滑行过程中的滑行艇的重量依靠水动升力和静浮力联合支撑。同时由于滑行艇高速航行时具有强非线性特点，流动在折角线和方尾底部发生分离现象，其水动力作用较为复杂，使得常规排水型船舶水动力的预报方法不再适用于滑行艇。本节考虑在固定大地坐标系下研究滑行艇在波浪中的运动，x轴在未扰动的自由面上，正方向选定为滑行艇前进的方向，z轴垂直向下为正方向。图1为滑行艇在波浪中运动的受力示意图。根据牛顿第二定律，运动方程可以表达为：

$ \left\{ \begin{aligned} &m{{\ddot x}_{_{{{CG}}}}} = {T_x} - N\sin \theta - {D_F}\cos \theta，\\ &m{{\ddot z}_{_{{{CG}}}}} = {T_z} - N\cos \theta + {D_F}\sin \theta + mg，\\ &{I_{55}}\ddot \theta = N{x_c} - {D_F}{x_d} + T{x_p} 。\end{aligned} \right. $

(1)

式中：m为滑行艇质量；I₅₅为滑行艇纵摇惯性矩；N为法向力；D为摩擦阻力；$ {T_x} $为滑行艇x方向受到的推力；$ {T_z} $为滑行艇z方向受到的升力；$ {x_c} $为重心与法向压力中心的距离；$ {x_d} $为重心到摩擦阻力作用线的距离；$ {x_p} $为推力关于重心位置的转动力臂。

滑行艇以恒定速度u在波浪中航行，假定作用于滑行艇的推力与阻力平衡。本文基于Zarnick模型，采用2D+T理论结合楔形体入水砰击模型，同时基于长波近似假设，船体处于滑行状态的绕射波相对于船体较小，在计算波浪激励力时忽略绕射力的影响，仅考虑入射波浪力的作用，建立滑行艇的纵荡、垂荡与纵摇模态耦合的非线性运动方程如下：

$ \left\{ \begin{aligned} &m{{\ddot x}_{_{{{CG}}}}} = - {F_{{{HD}}}}\sin \theta - {F_{{{CD}}}}\sin \theta，\\ &m{{\ddot z}_{_{{{CG}}}}} = - {F_{{{HD}}}}\cos \theta - {F_{{{FK}}}} - {F_{{{CD}}}}\cos \theta - {F_{{{HS}}}} + mg，\\ &{I_{55}}\ddot \theta = {M_{{{HD}}}} + {M_{{{FK}}}} + {M_{{{CD}}}} + {M_{{{HS}}}}。\end{aligned} \right. $

(2)

式中：$ x{_{CG}} $、$ z{_{CG}} $与$ \theta $分别为滑行艇重心位置处的纵荡、垂荡与纵摇运动的位移；F_HD为水动力；F_FK为入射波浪力；F_CD为粘性升力；F_HS为静回复力；M_HD、M_FK、M_CD与M_HS均为对应的纵摇力矩。

2 非线性时域运动模型 2.1 2D+T理论

图2为滑行艇在波浪中运动的2D+T理论示意图。在固定于大地的横平面内，将艇体对此平面内流体的扰动转化为对称楔形体入水问题的研究，确定作用在二维截面处的力，沿船长方向积分计算滑行艇受到的作用力。

滑行艇的各个横剖面随时间推移依次穿过该横平面：假设初始时刻t=t₀，滑行艇的横剖面尚未入水，此时正穿过横平面，相当于剖面恰好位于自由面上方；随时间推移，t=t₁时刻的横剖面接触自由面，折角线未浸湿，称为干折角线情况（dry chine）；t=t₂时刻的横剖面折角线浸湿，流动在折角线处发生分离，称为湿折角线情况（wetted chine）。因此，研究高速滑行艇在波浪中的运动问题转化为在固定于大地的横平面视角下的剖面入水流动问题，该方法的关键在于剖面入水过程中水动力的计算。

本文在计算入射波浪力时，同时考虑滑行艇瞬时湿表面与入射波垂向速度$ {w_{{z}}} $的影响。由于入射波水平速度相对于滑行艇前进速度较小，忽略入射波水平速度对船体运动的影响。平行于龙骨的切向速度$ U\left( {\xi ,t} \right) $与垂直于龙骨的法向速度$ V\left( {\xi ,t} \right) $可以表示为前进速度u、纵向运动$ \left( {{\eta _{\text{3}}},{\eta _{\text{5}}}} \right) $与入射波垂向速度$ {w_{\text{z}}} $的函数：

$ \left\{ \begin{aligned} &U\left( {\xi ,t} \right) = u\cos {\eta _5} - \left( {{{\dot \eta }_3} - {w_{{z}}}} \right)\sin {\eta _5} ，\\ &V\left( {\xi ,t} \right) = u\sin {\eta _5} + \left( {{{\dot \eta }_3} - {w_{{z}}}} \right)\cos {\eta _5} - {{\dot \eta }_5}\xi 。\end{aligned} \right. $

(3)

2.2 对称楔形体附加质量

本文采用2D+T理论计算滑行艇迎浪航行运动时，二维横剖面垂荡的附加质量则是一个重要参数，其计算基于对称楔形体入水砰击模型。本文采用Zarnick给出的对称楔形体附加质量的表达式：

$ \left\{ \begin{aligned} &{m_{{a}}} = {k_{{a}}}\frac{{\text{π}} }{2}\rho {b^2}，\\ &\dfrac{{{\rm{d}}{m_{{a}}}}}{{{\rm{d}}t}} = {{\dot m}_{{a}}} = {k_{{a}}}{\text{π}} \rho b\dot b 。\end{aligned} \right. $

(4)

式中：$ {k_{{a}}} $为附加质量系数；$ b $为剖面的瞬时半宽值。

Zarnick基于Wagner理论采用附加质量系数$ {k_{{a}}} $=1.0计算剖面附加质量。本文基于Hosseinzadeh等^[17]提出的附加质量系数拟合表达式，结合Zhao和Faltinsen的数值结果，给出修正后的附加质量系数表达式。图3为不同方法获得的楔形体的附加质量系数的比较。可以看出，本文计算的附加质量系数处于Zhao和Faltinsen与Vorus获得的数值结果之间，吻合较好。计算式为：

$ {k_{{a}}} = 0.35 + \left[ {\frac{{0.73}}{{\cos {\beta _{{\text{rad}}}}}} \times {{\left( {1 - \frac{{2{\beta _{{\text{rad}}}}}}{{\text{π}} }} \right)}^{{9 \mathord{\left/ {\vphantom {9 2}} \right. } 2}}}} \right]。$

(5)

式中：$ {\beta _{{\text{rad}}}} $为弧度值的底部斜升角。

本文分为2种情况计算剖面的附加质量：针对湿折角线情况，即浸没深度超过折角线位置时，附加质量表达式为：

$ \left\{ \begin{gathered} {m_{{a}}} = {k_{{a}}}\frac{{\text{π}} }{2}\rho b_{\max }^2 = {\text{constant}}，\\ {{\dot m}_{{a}}} = 0。\\ \end{gathered} \right. $

(6)

式中：$ {b_{\max }} $为折角线处的半宽值。

针对干折角线情况，考虑自由面爬升的影响，则剖面浸湿深度d与瞬时半宽值b存在如下关系：

$ b = {C_{{{pu}}}}d\cot \beta。$

(7)

式中：$ \beta $为剖面的底部斜升角；$ {C_{{{pu}}}} $为飞溅系数。Zarnick基于Wagner理论采用飞溅系数$ {C_{{{pu}}}} = {{\text{π}} \mathord{\left/ {\vphantom {{\text{π}} 2}} \right. } 2} $。本文采用Payne给出的飞溅系数表达式：

$ C_{pu}=\left\{\begin{aligned} & {{\text{π}}}/2，\beta=0，\\ & {\text{π}}/2-\beta\left(1-2/{\text{π}}\right)，0 < \beta < {\text{π}}/2，\\ & 1.0，\beta={\text{π}}/2。\end{aligned}\right. $

(8)

船体剖面中心线处$ P\left( {\xi ,0,\zeta } \right) $的坐标位置表示为：

$ \left\{ \begin{gathered} x(t) = {x_{{{cg}}}}(t) + \xi \cos {\eta _5}(t) + \zeta \sin {\eta _5}(t)，\\ z(t) = {z_{{{cg}}}}(t) - \xi \sin {\eta _5}(t) + \zeta \cos {\eta _5}(t)。\\ \end{gathered} \right. $

(9)

入射波的波长相对于高速滑行艇的吃水较大，针对小波陡工况，船体剖面垂直于基线的浸湿深度可表示为：

$ \begin{split} {\rm{d}}(t) =& \frac{{z(t) - {\zeta _0}(\xi ,t)}}{{\cos {\eta _5}(t) - \nu \sin {\eta _5}(t)}} = \\ &\frac{{{z_{{{cg}}}}(t) - \xi \sin {\eta _5}(t) + \zeta \cos {\eta _5}(t) - {\zeta _0}(\xi ,t)}}{{\cos {\eta _5}(t) - \nu \sin {\eta _5}(t)}} 。\end{split} $

(10)

式中：$ \nu = {{2{\zeta _a}} \mathord{\left/ {\vphantom {{2{\zeta _a}} \lambda }} \right. } \lambda } $，则剖面的浸湿深度的导数可表示为：

$ \begin{split} &\dot {\rm{d}}(t) = \frac{{\dot z(t) - {{\dot \zeta }_0}(\xi ,t)}}{{\cos {\eta _5}(t) - \nu \sin {\eta _5}(t)}} + \\ & \frac{{z(t) - {\zeta _0}(\xi ,t)}}{{\left[ {\cos {\eta _5}(t) - \nu \sin {\eta _5}(t)} \right]{}^{\text{2}}}}\frac{{\partial \left[ {\cos {\eta _5}(t) - \nu \sin {\eta _5}(t)} \right]}}{{\partial t}}。\end{split} $

(11)

由于$ z(t) - {\zeta _0}(\xi ,t) $是小量，上式可以简化为：

$ \dot {\rm{d}}(t) \approx \frac{{\dot z(t) - {{\dot \zeta }_0}(\xi ,t)}}{{\cos {\eta _5}(t) - \nu \sin {\eta _5}(t)}}。$

(12)

代入附加质量偏导数的表达式，则有：

$ {\dot m_{{a}}} = {k_{{a}}}{\text{π}} \rho b\left( {{C_{{{pu}}}}\cot \beta } \right)\frac{{\dot z(t) - {{\dot \zeta }_0}(\xi ,t)}}{{\cos {\eta _5}(t) - \nu \sin {\eta _5}(t)}}。$

(13)

2.3 滑行艇波浪力与运动方程推导

根据动量定理，作用于单位长度的船体剖面，与流体动量变化相关的水动升力可以表示为：

$ \begin{split} &{f_{{{HD}}}}\left( {\xi ,t} \right) = \frac{{\text{d}}}{{{\text{d}}t}}\left[ {{m_{{a}}}\left( {\xi ,t} \right)V\left( {\xi ,t} \right)} \right]= \\ &{m_{{a}}}\left( {\xi ,t} \right)\dot V + V{{\dot m}_{{a}}}\left( {\xi ,t} \right)- U\left( {\xi ,t} \right)\frac{{\partial \left[ {{m_{{a}}}\left( {\xi ,t} \right)V\left( {\xi ,t} \right)} \right]}}{{\partial \xi }} 。\end{split} $

(14)

与横向流动相关的垂向粘性升力可以表示为：

$ {f_{{{CD}}}}\left( {\xi ,t} \right) = {C_{{D}}}\rho b{V^2}\left( {\xi ,t} \right)。$

(15)

式中：$ {C_{{D}}} $为横向流动阻力系数，其数值与横剖面的底部斜升角的大小相关，本文采用Keuning给出的关系式$ {C_{{D}}} = 1.33\cos \beta $确定。

将上述作用力沿船长方向积分，则纵荡模态的水动力可表达为：

$ {\begin{split} &{F_{{{X - HD}}}} = - \displaystyle\int\limits_{L(t)} {{f_{{{HD}}}}\left( {\xi ,t} \right){\rm{d}}\xi \sin {\eta _5}(t)} - \displaystyle\int\limits_{L(t)} {{f_{{{CD}}}}\left( {\xi ,t} \right){\rm{d}}\xi \sin {\eta _5}(t)}= \\ &\left\{ {\displaystyle\int\limits_{L(t)} { - {m_{{a}}}\sin {\eta _5}{{\ddot \eta }_1}{\rm{d}}\xi - } \displaystyle\int\limits_{L(t)} {{m_{{a}}}\cos {\eta _5}{{\ddot \eta }_3}{\rm{d}}\xi } + \displaystyle\int\limits_{L(t)} {{m_{{a}}}\xi {{\ddot \eta }_5}{\rm{d}}\xi } } + \right. \\ &\displaystyle\int\limits_{L(t)} {{m_{{a}}}{{\dot \eta }_5}\left( {{{\dot \eta }_3}\sin {\eta _5} - {{\dot \eta }_1}\cos {\eta _5}} \right){\rm{d}}\xi } + \displaystyle\int\limits_{L(t)} {{m_{{a}}}\dfrac{{{\text{d}}{w_{{z}}}}}{{{\text{d}}t}}\cos {\eta _5}{\rm{d}}\xi }- \\ &\displaystyle\int\limits_{L(t)} {{m_{{a}}}{w_{{z}}}\sin {\eta _5}{{\dot \eta }_5}{\rm{d}}\xi } - \displaystyle\int\limits_{L(t)} {{m_{{a}}}V\dfrac{{\partial {w_{{z}}}}}{{\partial \xi }}\sin {\eta _5}{\rm{d}}\xi } + \displaystyle\int\limits_{L(t)} {{m_{{a}}}U\dfrac{{\partial {w_{{z}}}}}{{\partial \xi }}\cos {\eta _5}{\rm{d}}\xi } - \\ &\left. { UV{m_{{a}}}\left| {_{{\text{stern}}}} \right. - \displaystyle\int\limits_{L(t)} {{{\dot m}_{{a}}}V{\rm{d}}\xi } } \right\}\sin {\eta _5} - \int\limits_{L(t)} {\left\{ {{C_{{D}}}\rho b{V^2}} \right\}{\rm{d}}\xi } \sin {\eta _5}。\\[-13pt] \end{split} }$

(16)

垂荡模态的水动力可表达为：

$ {\begin{split} &{F_{{{Z - HD}}}} = - \displaystyle\int\limits_{L(t)} {{f_{{{HD}}}}\left( {\xi ,t} \right){\rm{d}}\xi \cos {\eta _5}(t)} - \displaystyle\int\limits_{L(t)} {{f_{{{CD}}}}\left( {\xi ,t} \right){\rm{d}}\xi \cos {\eta _5}(t)} = \\ &- \displaystyle\int\limits_{L(t)} {\left\{ {{m_{{a}}}\dot V + V{{\dot m}_{{a}}} - U\dfrac{{\partial {m_{{a}}}V}}{{\partial \xi }}} \right\}{\rm{d}}\xi } \cos {\eta _5} - \displaystyle\int\limits_{L(t)} {\left\{ {{C_{{D}}}\rho b{V^2}} \right\}{\rm{d}}\xi } \cos {\eta _5} = \\ &\left\{ {\int\limits_{L(t)} { - {m_{{a}}}\sin {\eta _5}{{\ddot \eta }_1}{\rm{d}}\xi - } \displaystyle\int\limits_{L(t)} {{m_{{a}}}\cos {\eta _5}{{\ddot \eta }_3}{\rm{d}}\xi } + \displaystyle\int\limits_{L(t)} {{m_{{a}}}\xi {{\ddot \eta }_5}{\rm{d}}\xi } } + \right. \\ &\displaystyle\int\limits_{L(t)} {{m_{{a}}}{{\dot \eta }_5}\left( {{{\dot \eta }_3}\sin {\eta _5} - {{\dot \eta }_1}\cos {\eta _5}} \right){\rm{d}}\xi } + \displaystyle\int\limits_{L(t)} {{m_{{a}}}\dfrac{{{\text{d}}{w_{{z}}}}}{{{\text{d}}t}}\cos {\eta _5}{\rm{d}}\xi } - \\ &\displaystyle\int\limits_{L(t)} {{m_{{a}}}{w_{{z}}}\sin {\eta _5}{{\dot \eta }_5}{\rm{d}}\xi } - \displaystyle\int\limits_{L(t)} {{m_{{a}}}V\dfrac{{\partial {w_{{z}}}}}{{\partial \xi }}\sin {\eta _5}{\rm{d}}\xi } + \displaystyle\int\limits_{L(t)} {{m_{{a}}}U\dfrac{{\partial {w_{{z}}}}}{{\partial \xi }}\cos {\eta _5}{\rm{d}}\xi } -\\ &\left. { UV{m_{{a}}}\left| {_{{\text{stern}}}} \right. - \displaystyle\int\limits_{L(t)} {{{\dot m}_{{a}}}V{\rm{d}}\xi } } \right\}\cos {\eta _5} - \displaystyle\int\limits_{L(t)} {\left\{ {{C_{{D}}}\rho b{V^2}} \right\}{\rm{d}}\xi } \cos {\eta _5} {。}\\[-3pt] \end{split} } $

(17)

纵摇模态的水动力可表达为：

$ \begin{split} &{M_{{{Y - HD}}}} = \int\limits_{L(t)} {{f_{{{HD}}}}\left( {\xi ,t} \right)\xi {\rm{d}}\xi } + \int\limits_{L(t)} {{f_{{{CD}}}}\left( {\xi ,t} \right)\xi {\rm{d}}\xi } =\\ & - \int\limits_{L(t)} {{m_{{a}}}{\xi ^2}{{\ddot \eta }_5}{\rm{d}}\xi } + \int\limits_{L(t)} {{m_{{a}}}\xi \cos {\eta _5}{{\ddot \eta }_3}{\rm{d}}\xi }+\int\limits_{L(t)} {{m_{{a}}}\xi \sin {\eta _5}{{\ddot \eta }_1}{\rm{d}}\xi } + \\ & \int\limits_{L(t)} {{m_{{a}}}\xi {{\dot \eta }_5}\left( {{{\dot \eta }_1}\cos {\eta _5} - {{\dot \eta }_3}\sin {\eta _5}} \right){\rm{d}}\xi }-\int\limits_{L(t)} {{m_{{a}}}\frac{{{\text{d}}{w_{{z}}}}}{{{\text{d}}t}}\cos {\eta _5}\xi {\rm{d}}\xi } + \\ & \int\limits_{L(t)} {{m_{{a}}}{w_{{z}}}\sin {\eta _5}{{\dot \eta }_5}\xi {\rm{d}}\xi } +\int\limits_{L(t)} {{m_{{a}}}V\frac{{\partial {w_{{z}}}}}{{\partial \xi }}\sin {\eta _5}\xi {\rm{d}}\xi } -\\ & \int\limits_{L(t)} {{m_{{a}}}U\frac{{\partial {w_{{z}}}}}{{\partial \xi }}\cos {\eta _5}\xi {\rm{d}}\xi } +UV{m_{{a}}}\xi \left| {_{{\text{stern}}}} \right. +\\ & \int\limits_{L(t)} {{m_{{a}}}UV{\rm{d}}\xi } + \int\limits_{L(t)} {{{\dot m}_{{a}}}V\xi {\rm{d}}\xi } + \int\limits_{L(t)} {\left\{ {{C_{{D}}}\rho b{V^2}\xi } \right\}{\rm{d}}\xi } 。\\[-3pt] \end{split} $

(18)

分离加速度相关项，写成与滑行艇附加质量相关的形式：

$ {\left\{ \begin{gathered} {A_{11}} = \int\limits_{L(t)} {{m_{{a}}}{{\sin }^2}{\eta _5}{\rm{d}}\xi } \;{A_{31}} = \int\limits_{L(t)} {{m_{{a}}}\sin {\eta _5}\cos {\eta _5}{\rm{d}}\xi } \;{A_{51}} = \int\limits_{L(t)} {{m_{{a}}}\xi \sin {\eta _5}{\rm{d}}\xi }{，} \\ {A_{13}} = \int\limits_{L(t)} {{m_{{a}}}\sin {\eta _5}\cos {\eta _5}{\rm{d}}\xi } \;{A_{33}} = \int\limits_{L(t)} {{m_{{a}}}{{\cos }^2}{\eta _5}{\rm{d}}\xi }\; {A_{53}} = \int\limits_{L(t)} {{m_{{a}}}\xi \cos {\eta _5}{\rm{d}}\xi } {，} \\ {A_{15}} = \int\limits_{L(t)} {{m_{{a}}}\sin {\eta _5}\xi {\rm{d}}\xi } \;{A_{35}} = \int\limits_{L(t)} {{m_{{a}}}\cos {\eta _5}\xi {\rm{d}}\xi } \;{A_{55}} = \int\limits_{L(t)} {{m_{{a}}}{\xi ^2}{\rm{d}}\xi }\,{。} \\ \end{gathered} \right.} $

(19)

则水动力可写成如下形式：

$ \left\{ \begin{gathered} {F_{{{X - HD}}}} = - {A_{11}}{{\ddot \eta }_1} - {A_{13}}{{\ddot \eta }_3} + {A_{15}}{{\ddot \eta }_5} + {{F'}_{{{X - HD}}}}，\\ {F_{{{Z - HD}}}} = - {A_{31}}{{\ddot \eta }_1} - {A_{33}}{{\ddot \eta }_3} + {A_{35}}{{\ddot \eta }_5} + {{F'}_{{{Z - HD}}}}，\\ {M_{{{Y - HD}}}} = {A_{51}}{{\ddot \eta }_1} + {A_{53}}{{\ddot \eta }_3} - {A_{55}}{{\ddot \eta }_5} + {{M'}_{{{Y - HD}}}}。\\ \end{gathered} \right. $

(20)

同时考虑静水力作用，采用波浪截切方式获取滑行艇的瞬时湿表面以求解非线性回复力。将上式代入运动方程，可写出矩阵形式的滑行艇运动方程如下：

$ {\left[ \begin{gathered} m + {A_{11}}{\text{ }}{A_{13}}{\text{ }} - {A_{15}} \\ {\text{ }}{A_{31}}{\text{ }}m + {A_{33}}{\text{ }} - {A_{35}} \\ {\text{ }} - {A_{51}}{\text{ }} - {A_{53}}{\text{ }}{I_{55}} + {A_{55}} \\ \end{gathered} \right]\left[ \begin{gathered} {{\ddot \eta }_1} \\ {{\ddot \eta }_3} \\ {{\ddot \eta }_5} \\ \end{gathered} \right] = \left[ \begin{gathered} {\text{ }}{{F'}_{{\text{X - HD}}}} \\ {{F'}_{{\text{Z - HD}}}} - {F_{{\text{HS}}}} + mg \\ {\text{ }}{{M'}_{{\text{Y - HD}}}} + {M_{{\text{HS}}}} \\ \end{gathered} \right]。} $

(21)

本文采用四阶Runge-Kutta法求解上述滑行艇运动方程。基于上述理论，本文开发了一套高速滑行艇在迎浪中非线性运动计算程序，根据流程图4进行滑行艇时域运动计算求解，获得船体运动响应结果。

3 非线性运动数值结果与分析

本文以Fridsma滑行艇为研究对象^[20]，Fridsma滑行艇主尺度如表1所示，计算工况为无限水深迎浪规则波。图5为Fridsma滑行艇三维面元网格与纵向划分的二维剖面示意图，可以看出Fridsma滑行艇具有恒定的底部斜升角。选取航速Fn=1.2，入射波$ \lambda /L $=2.0工况为例，时间步长取为∆t=Te/100，研究本文2D+T方法中纵向划分剖面数量对时域运动计算结果的影响。

图6为4种不同剖面数量对应的Fridsma滑行艇垂荡与纵摇运动时历，可以看到运动结果随剖面数量的增加而逐渐收敛。此外，滑行艇在高速航行时的平衡位置在静浮吃水位置以上，这与中高速排水型船舶的运动具有显著差异。图7为不同剖面数量对应的Fridsma滑行艇的运动响应结果，计算结果随剖面数量的增加而趋于稳定一致。因此在本文后续计算中，三维滑行艇艇体纵向离散的剖面数量选取为60。

图8与图9为本文采用2D+T方法计算Fridsma滑行艇在航速Fn=1.2时迎浪工况下的垂荡与纵摇运动结果，同时与试验数据^[20]和BAHRAMI等^[21]基于CFD模拟获得的结果的对比。可以看出，本文数值方法计算的纵摇运动与试验结果吻合良好；对于垂荡运动，数值预报结果的变化趋势与试验值相一致，但是大波高工况下的垂荡峰值预报结果偏大。

图10为Fridsma滑行艇在航速Fn=0.6时的垂荡与纵摇运动结果，数值结果与试验数据^[20]趋势吻合较好，进一步验证了本文数值方法的准确性。

4 结　语

本文在大地坐标系下，基于长波近似假设忽略绕射力的影响，建立滑行艇的纵荡、垂荡与纵摇模态耦合的非线性运动方程，给出了修正的楔形体附加质量系数表达式用于确定二维横剖面垂荡的附加质量；应用波切法获取滑行艇的瞬时湿表面以求解非线性回复力；在2D+T势流理论的框架下结合对称楔形体入水砰击模型，详细推导了高速滑行艇在迎浪中的非线性时域运动方程，并自主开发了数值计算程序。

本文选取Fridsma滑行艇作为研究对象，探讨了滑行艇纵向离散剖面数量对运动结果的影响，确定了合适的纵向离散剖面数量。选取航速为Fn为0.6和1.2工况，计算了滑行艇在迎浪规则波中的运动响应结果并与试验数据进行了对比，验证了所开发的时域运动程序的准确性。滑行艇高速航行时的平衡位置在静浮位置以上，这与中高速排水型船舶的运动具有显著差异，体现了滑行艇高速滑行过程中垂向水动升力的影响。本文的研究工作对高速滑行艇运动分析具有重要的理论意义和工程实用价值。