船舶在波浪中航行时，波面会导致船体水线面较静水中发生明显变化，进而改变船体的排水体积、水线面面积及其绕纵轴的转动惯性矩。当船随浪(或顺浪)航行时，由于船波遭遇周期变长，波面形状对船体复原特性的影响更为显著，可能明显地降低船舶的稳性，这种现象称为纯稳性丧失。目前，国际海事组织(IMO)正在制定的“第二代完整稳性横准”，将纯稳性损失作为5种波浪中船舶稳性失效的模式之一^[1-2]。

对波浪中船舶稳性损失的研究由来已久。早在1861年，Froude就已注意到了波浪中船舶的实际复原力矩可能与理想状态下的结果有着显著的不同，并试图从定量上估算波浪的影响。

Paulling等基于在旧金山湾进行的船模实验，首次提出了波峰线与船舶中线垂直状态下船舶稳性的实用计算方法。同时从理论上指出，波浪中稳性与静水中的稳性有很大的不同，主要是因为波浪中船型的淹湿形状改变和波浪引起的作用在船体上压力改变产生的，当船舶处于波峰时稳性有可能急剧减小，以致有可能导致倾覆^[3]。

Hamamoto发展了艉斜浪中的稳性计算方法，只是存在大横倾角时其计算结果有时将发散^[4]。Chiu等基于切片理论，采用瞬时湿表面的时域分析法研究了波浪对船舶稳性的影响^[5]。Fang等在讨论斜波中船舶稳性时，在波浪上船型淹湿表面积的计算中，不再采用“Froude-Krylov”假定，在计算中将波浪绕射和散射的影响计入，使问题大为完善^[6]。Kreuzer等在横摇、垂荡和纵摇3个自由度运动状态下计算淹湿表面积上的压力，以研究船舶稳性和船舶倾覆问题^[7]。

国内相关研究主要有：孔祥金等改善了Hamamoto艉斜浪中的稳性计算方法，但是这些计算并不包括绕射波的影响^[8]。董艳秋等应用时域分析法，加入了瞬时湿表面积的影响，给出了船舶在波浪中横摇稳性的计算方法^[9]。李子富等考虑纵倾对船舶横剖面左右浸深的影响，提出了一种新的计算船舶发生大倾角横摇时的复原力臂的方法^[10]。鲁江等基于波浪中船舶受力及力矩的瞬时静平衡方程，提出了波浪中船体横剖面与波面交点的求法及各浸水剖面面积的通用计算法，并基于“Froude-Krylov”假定研究了规则波、斜浪长峰规则波以及随机波中船舶的复原力变化，并对参数横摇进行数值模拟^[11-14]。卜淑霞和李红霞等采用垂荡-横摇-纵摇三自由度耦合数学模型，基于“Froude-Krylov”假定，考虑瞬时横摇角度对瞬时湿表面的影响计算辐射力和绕射力，对波浪中的船舶稳性进行计算，并对顶浪及斜浪规则波中的参数横摇进行数值模拟^[15-16]。胡开业等采用NURBS方法对船舶曲面进行拟合，计算了任意倾斜姿态下船体的几何要素及船舶的随浪稳性，分析了波长、波高、浪向角和波浪与船舶的相对位置等参数对船舶随浪稳性的影响^[17]。

上述研究多基于“Froude-Krylov”假定，认为船舶周围的水质点运动不会受到船舶影响，虽然有的研究针对辐射力和绕射力做了修正，但与实际情况仍有差距，尤其是在航速较高时，船行波对船体的瞬时湿表面积及稳性的影响，无法忽略。

随着计算机技术及粘性流理论的发展，计算流体力学方法(computational fluid dynamics，CFD)已经在船舶水动力学领域得到广泛运用，在船舶阻力等问题研究中，已经成为实验方法的有利补充。CFD方法基于粘性流理论，可以精确模拟船体在波浪航行时的流固耦合运动，同时可以处理甲板上浪、波浪破碎等强非线性现象，这使得研究船舶在波浪中的综合性能成为可能。目前，国内外学者基于CFD方法对船舶在波浪中的性能做了大量研究，但这些研究主要集中于顶浪、斜顶浪及横浪状态，对于船舶随浪航行性能的研究相对较少。

本文基于CFD方法，对船舶随浪航行时，稳性损失最严重的情况——船舶骑浪航行，波峰位于船舯时的稳性进行研究。基于边界造波法构建适合船舶骑浪航行的规则波数值波浪水池；在此基础上，综合计入船行波及波浪影响，数值模拟船舶骑浪航行流场，对不同横倾角下的浮力及复原力矩进行计算，求得复原力臂，并与水池实验数据进行比对，以验证本文方法的精度。本文为研究船舶随浪航行时的稳性性能提供了新的思路，也丰富了基于CFD的船舶波浪性能研究方法。

1 数学模型 1.1 控制方程

将水视为不可压缩粘性流体，控制方程有连续性方程和动量方程(Navier-Stokes方程)，其张量表示为

$\frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \left( {\rho {u_i}} \right)}}{{\partial {x_i}}} = 0$

(1)

$\begin{array}{l} \frac{{\partial \left( {\rho {u_i}} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {u_i}{u_j}} \right)}}{{\partial {x_j}}} = \\ - \frac{{\partial P}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left[ {\mu \left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)} \right] + \rho {f_i} \end{array}$

(2)

式中：u_i、u_j为流体速度矢量u在x_i、x_j方向的分量，t为时间，P为压力，ρ为流体密度，f_i为质量力，μ为流体动力粘性系数。

1.2 波浪数学模型

本文采用船舶耐波性水池实验中常用的微幅波模型模拟波浪。当船舶骑浪航行时，航向与波浪传播方向相同，航速与波速相等，波面与船舶位置相对固定。假设船舶静止，建立以船舶重心为原点，X轴正向指向船尾，Y轴正向指向船舶右舷，Z轴正向垂直向上的船体坐标系，基于相对运动原理，可得波浪波高方程为^[3]

$\eta = {A_0}{\rm{cos}}\left[ {k\left( {x + {U_0}t} \right) - \omega t + {\varepsilon _0}} \right]$

(3)

速度方程为

$\left\{ {\begin{array}{*{20}{l}} {u = {A_0}\omega {{\rm{e}}^{kz}}{\rm{cos}}\left[ {k\left( {x + {U_0}t} \right) - \omega t + {\varepsilon _0}} \right] + {U_0}}\\ {v = 0}\\ {w = {A_0}\omega {{\rm{e}}^{kz}}{\rm{sin}}\left[ {k\left( {x + {U_0}t} \right) - \omega t + {\varepsilon _0}} \right]} \end{array}} \right.$

(4)

式中：u、v、w分别为x、y、z三个方向的速度，A₀为波幅，k为波数，ω为波浪圆频率，U₀为航速，ε₀为初始相位。

2 数值模拟方案

船舶骑浪航行时，当波峰位于船舯时，稳性降低最为严重。针对此种工况，本文采用商用CFD软件Fluent，对船模在静水及3种不同波高波浪中骑浪航行时的周围流场进行了数值模拟。进行数值模拟时，分别将船模以不同横倾角度(从5°~80°，间隔5°)置于流场中，计算船模在不同横倾条件下的复原力矩，绘制船模稳性曲线。通过稳性曲线变化情况，分析船行波及波浪对船舶骑浪航行稳性的影响。

2.1 船模选择及实验参数

选择某瘦长型船模为研究对象，其主要参数为船长3.000 m，船宽0.329 m，吃水0.096 m，排水量0.049 t，方形系数0.494，缩尺比1:44。数值模型如图 1所示。

船舶在骑浪状态下，船长等于波长，航速等于波速。同时，为考察波高对稳性影响，分别进行不同波高条件下实验，由此确定波浪基本参数如表 1所示。

2.2 计算域及网格划分

计算域为长方体，长、宽、深为5L×2L×1.5L(L为船长)，其中入口距船艏1.5L，船艉距消波区1.5L，另有长1L的消波区，水深1L，自由面距上边界0.5L，示意图如图 2所示。

采用混合网格，在船体周围区域采用非结构网格以较好地表现船型，其余区域采用结构网格以减少网格数量并提高计算效率，在自由面附近进行网格加密以满足造波边界要求，总网格数180万，舰艏部分及自由面附近网格如图 3所示。

2.3 边界条件设置

采用边界造波法模拟波浪，计算域的边界条件具体设置如下：

入口边界——速度入口，按照式(3)、(4) 给定波高及速度值；

出口边界——压力出口，设置静水压力；

船体——壁面，有剪切力无滑移；

外边界(包含水池底部、顶部及侧壁)——壁面，剪切力为0。

采用VOF方法模拟自由面，RNG k-ε模型模拟湍流，采用壁面函数法模拟近壁面流动。

3 数值模拟结果及分析 3.1 波浪数值模拟结果

为验证数值造波方法的精度，采用不放船模的计算域进行数值造波，计算50 s流动时间后，瞬时波面的模拟值与理论值的对比如图 4所示。图中可以看出，工作区波浪模拟值与理论值吻合良好，消波区波浪衰减明显。经计算，波峰位置的波高平均误差为1.73%，波谷位置的波高平均误差为-3.02%，波长误差为-0.59%。

3.2 复原力臂计算

在数值模拟过程中记录船模所受浮力(F)及复原力矩(M)，典型的力及力矩变化曲线如图 5所示。从图中可以看出，力矩和力数值很快收敛，在小范围内波动，这主要因为船体周围流场不断变化所致。取力矩和力稳定后一段时间的均值作为该种工况下的浮力及复原力矩数值，两者比值即为复原力臂。经缩尺比换算后，可得不同工况下的复原力臂如表 2所示(表中A₀为实船对应波幅值)。

3.3 数值模拟结果分析

文献[18]给出了该型船在静水中及在波幅A₀=4.4 m波浪中骑浪航行时的稳性曲线，具体数值如表 3所示。其中静水稳性曲线通过水池实验得到，骑浪稳性曲线通过计算方法得到，在计算过程中没有考虑船行波影响。

3.3.1 数值模拟精度

图 6是静水稳性曲线水池实验值与数值模拟值的对比图。各角度复原力臂数值模拟值与水池实验值的相对误差如表 4所示。从表 4中可以看出，数值模拟值与水池实验值吻合良好，相对误差基本在2%~3%，最大值小于5%；靠近曲线两端误差有所增大，是因为复原力臂值本身较小，容易导致较大的相对误差。

3.3.2 船行波对稳性损失的影响

图 7是船模在静水以及A₀=4.4 m波浪中骑浪航行时的稳性曲线数值模拟值与计算值的对比图。从图中可以看出，船模骑浪航行时，复原力臂明显减小，稳性明显降低；但同计算值相比，数值模拟值的减小幅度要小于计算值。

图 8是船模骑浪航行时的瞬时波面图，从图中可以看出，船行波改变了骑浪航行波浪场，使得船模的湿表面积及水下排水体积发生变化，改变了船模的受力，使船模的稳性发生变化，在一定程度上削弱了波浪的影响。实际上，按照缩尺比计算，当该型船骑浪航行时，航速需达到27.9 kn，如此高航速下产生的船行波的影响是不应该忽略的。而计算值并没有考虑船行波的影响而做相应的修正，这是导致两者差别的主要原因。

3.3.3 波高对稳性损失的影响

图 9是不同波高骑浪航行时船模复原力臂对比图，不同角度的复原力臂减小值如表 5所示，可以看出：

1) 船模的稳性损失随波高增大而增加；

2) 当波高较小时(A₀=2.2 m)，由于船行波的影响，稳性损失并不是很大，基本上在10%以内；

3) 对同一波高而言，稳性损失基本上随横倾角增大而增大(个别横倾角较小时的复原力臂减小程度较大，可能与相对误差较大有关系)，尤其是当达到最大复原力臂之后，稳性损失尤为严重，这是在海上航行尤其需要注意的。

4 结论

1) 数值计算得到的静水复原力臂与水池实验值吻合良好，说明本文所用方法可以较为精确地计算船模复原力臂，这为CFD方法用于船舶的稳性研究提供了可靠的依据；

2) 大型船舶骑浪航行时，航速一般较高，船行波会削弱波浪对稳性损失的影响，在计算稳性损失的时候不应忽略船行波的影响；

3) 船舶稳性损失随波高增大而增大，对同一波高而言，横倾角越大，稳性损失越严重；

本文在进行数值计算时，并未考虑船舶垂荡与纵摇的影响，这与海上实际有一定差距，需要在以后的研究中加以改进，并进一步开展随浪航行横摇模拟的研究。