0 引　言

随着计算机性能的不断发展，计算流体力学（CFD）在船舶阻力、操纵性、耐波性和自航等方面得到更广泛应用，其计算方案得到进一步完善，计算效率得到大幅提升，计算结果也更加可靠。相较于传统模型试验中高昂模型成本和试验成本，数值计算主要依靠高性能计算机，成本相对低廉，可完成多工况同步计算，效率较高。更重要的是数值计算能够展示模型试验无法展示细节，为科学研究和工程应用提供更全面的支撑。

目前，国内外学者对系列船型或者标准船型采用不同手段和方法展开数值计算研究，反馈结果均与试验结果吻合较好。沈兴荣等^[1]通过对4 000 TEU巴拿马型集装箱船的自航数值模拟结果进行总结分析，证明了对船模进行自航试验数值模拟可行。何惠明等^[2]利用CFD对某艘艇在波浪中的自航进行数值模拟，得到该艇在规则波下自航状态的自航因子变化曲线。程宣凯等^[3]使用Shipflow对螺旋桨敞水试验、船舶阻力试验和自航试验进行数值模拟，螺旋桨的推力系数、扭矩系数，船模的阻力系数、自航因子及推进效率与试验值均吻合较好。Carrica P M等^[4]使用CFD Ship-lowa v4分别对DTMB5613船、KVLCC1船和KCS船的自航试验进行数值模拟，并利用DES湍流模型和动网格技术对KCS船的自航试验进行数值模拟，模拟结果与试验结果吻合较好。Castro A M等^[5]对KCS船进行了敞水数值模拟，根据敞水、阻力和自航的模拟计算结果，得到了KCS船的自航因子。

基于减小计算工作量提高效率考虑，上述研究主要以模型尺度为研究对象，对于实尺度下船舶自航性能的数值模拟研究较少，同时模拟研究中涉及到的标模主要为KCS船型。本文基于CFD商业软件STAR-CCM+对REGAL货船进行实尺度下的自航数值模拟，并将船体阻力、螺旋桨推力与扭矩的模拟结果与劳氏船级社的海上试验测量结果进行比较，得到实尺度自航数值模拟结果的误差范围，验证了实尺度下自航数值模拟的精度，可为实尺度下船舶自航数值计算的工程应用提供参考。

1 数值方法

CFD的数值计算主要涉及到质量守恒、动量守恒两大物理守恒定律^[5]，实尺度船舶自航数值模拟的本质为对N-S方程进行求解^[6-7]。其中流体的连续性方程和N-S方程分别如下：

$ \nabla \cdot u=0 ，$

(1)

$ \frac{\partial \rho u}{\partial t}+\nabla\cdot(\rho(u-u_g)u)-\nabla \cdot(\mu\nabla u)=-\nabla P \text{。}$

(2)

式中： $ \dfrac{\partial \rho u}{\partial t}$ 表示压力对时间的微分； $ \nabla P$ 为压力梯度； $ \nabla\cdot(\rho(u-u_g)u)$ 为扩散项；u为流体速度； $ \ \rho$ 为密度；t为时间；P为压力； $ \ \mu$ 为动力粘度系数。

数值模拟过程中，开启六自由度监测对船体自身运动进行求解，船体的浮船平动惯量和转动惯量表示如下：

$ (F_{CG})=\frac{m({\rm{d}}^2(x_{CG})_i)}{{\rm{d}}t^2}=m \ddot{x}_{CG}，$

(3)

$ (I_{CG})=\frac{{\rm{d}}}{{\rm{d}}t}\left(I_{ij}\frac{{\rm{d}}(\theta_{CG})_i}{{\rm{d}}t}\right)=I_{ij}\ddot{\theta}_{CG} 。$

(4)

式中：m为船体的质量；I_ij为船体的转动惯量； $ \ddot{x}_{CG}$ 和 $ \ddot{\theta}_{CG}$ 分别为船体平动的加速度和转动的角加速度； $ (x_{CG})_i$ 和 $ (\theta_{CG})_i$ 分别为船体平动的位移和转动的角位移； $ (F_CG)_i$ 和 $ (I_{CG})_i$ 分别为船体转动和平动作用在重心上的力和力矩。当i = 1，2，3分别表示纵荡、横荡和垂荡，当i = 4，5，6分别表示横摇、纵摇和首摇。

2 网格及边界条件设置

计算对象为劳氏船级社2016年国际研讨会所采用的REGAL货船，通过激光扫描船舶得到的三维模型数据，通过实船海上测试获得阻力、螺旋桨推力和扭矩等数据，为模拟计算结果的误差分析提供参考^[8]。表1和表2分别为REGAL货船和螺旋桨的主要参数，图1为螺旋桨的三维模型图。

为便于与已有试验结果比较，自航数值模拟的3种计算工况如表3所示。

2.1 计算域网格划分

数值模拟前，对计算域网格进行划分，实尺度下船舶自航的计算域大小与模型尺度下的不同，参考Castro A M等^[5]对KCS船自航数值计算设置，计算域大小取4L × 4L × 3L，其中L为船长，计算域如图2和图3所示。同时以桨轴为圆心建立一个完全包裹住螺旋桨的圆柱体，底面直径为1.3D，其中D为螺旋桨直径，此区域称为流体旋转域。将流体旋转域与螺旋桨区域进行布尔减运算后得到旋转计算域，以模拟螺旋桨的旋转运动，如图4所示。

采用CFD软件STAR-CCM+进行数值模拟，模型网格的基础尺寸取1/20L ~ 1/50L，本研究中网格基础尺寸取4 m，网格模型选择表面重构对网格曲率和兼容性进行优化，采用切割体网格单元生成器控制不同方向的网格增长率，棱柱层网格生成器控制边界层的网格y⁺大小控制在60～200之间，以增加自航模拟精度^[9]。同时运用自定义控制方式对流场变化较为剧烈的区域如水线面、船体周围的兴波等进行局部网格加密，以便于对流场细节的捕捉，网格的总数量为500万。图5为旋转计算域网格图，图6为计算域网格。

2.2 边界条件设置

计算入口设置为速度入口并给定初始速度；两侧计算域定义为对称平面；船体表面条件定义为无滑移的壁面，出口选择为压力出口并定义阻尼消波条件。在物理模型的选择上，选用K-Epsilon湍流、VOF波、三维、两层全y⁺壁面处理、欧拉多相、流体域体积（VOF）、重力和隐式不定常设置。在计算域内定义1个以船舶重心为原点，3个轴的方向分别与x，y，z轴平行的局部坐标系作为六自由度运动的参考坐标系，在六自由度DFBI设置中定义惯性矩和重心位置，在求解器里面设置时间步长、迭代求解次数及最大求解物理时间等，初始化算例后开始计算。

自航试验数值计算可以分为2个过程。首先保持螺旋桨固定，计算船体阻力，等待船体阻力收敛之后启动螺旋桨旋转运动，同时为了观测更为真实的桨后流场，自航计算中旋转运动选择Sliding Mesh方法，时间步长设置保证在每个时间步下螺旋桨的旋转角度为1°。

3 结果分析

3种计算工况下REGAL货船的实尺度自航数值模拟的阻力F，螺旋桨推力T和螺旋桨扭矩Q与劳氏船级社的海上试验测量结果对比如表4 ~ 表6所示。取实船航速为12 kn，螺旋桨转速91.1 r/min（工况2）下的阻力F，螺旋桨推力T和螺旋桨扭矩Q的对比曲线分别如图7 ~ 图9所示。

从上述对比结果可以得出：采用数值模拟方法对REGAL货船进行实尺度自航模拟整体效果良好，阻力最大误差为4.5%，推力和扭矩最大误差分别为8.6%与2.7%。

计算工况2下的船体阻力F，螺旋桨推力T和螺旋桨扭矩Q的收敛曲线分别如图10～图12所示。从图中可知，曲线均收敛良好。

螺旋桨附近的压力云图速度云图如图13和图14所示。

由于模型边界层的存在，螺旋桨模型与船体在装配面之间存在装配间隙，其长度为网格的最小长度，导致数值模拟流场与实际流场之前存在偏差，因此敞水效率会产生误差。同时螺旋桨的旋转运动导致了船体阻力的变化，而受到敞水效率误差的影响，螺旋桨扭矩结果误差较大，达到8.6%；若减小网格最小长度可使由此造成的误差减小，但会造成计算网格增加较多，计算效率降低。

4 结　语

采用STAR-CCM+软件对劳氏船级社2016年国际研讨会所采用的REGAL货船进行实尺度数值模拟，将3种不同工况船体阻力、螺旋桨推力及扭矩的数值模拟结果与海上实测结果进行对比分析，阻力最大误差为4.5%，推力和扭矩最大误差分别为8.6%与2.7%。研究结果表明实尺度自航数值模拟在一定范围内具备可行性，计算过程可为相关工程应用提供参考。