由于国际海事组织(IMO)要求估算因风和波浪引起的船舶附加阻力相对于平静海面阻力的大小^[1]。改变船型和环境条件，波浪附加阻力具有明显的变化。因此，估算附加阻力是航运业的一个重要问题。在过去的几十年中，用实验和数值方法对波浪引起的附加阻力问题进行了广泛的研究。

通过水池试验的研究，附加阻力与船舶运动、船速、波长、波高、航向、船体形状、船艏形状有关。船艏形状通常提供更大的附加阻力^[2]。Kashiwagi等^[3]研究了波高以及球鼻艏对船舶运动的影响。Steen^[4]专注于考虑低海况的短波区域的船舶运动实验。计算流体力学(CFD)的优点是不使用附加阻力分析公式或粘性效应经验值的情况下进行非线性计算，但是比势流理论方法花费更多时间。Sadat-Hosseini等^[5]通过建模分析不同船型的附加阻力和船舶运动，并将初步结果与实验数据进行了比较。黄德波等^[6]基于计算流体动力理论建立了数值波浪水池对单船和双体船迎浪航行进行了模拟研究。

短波条件下附加阻力是国内外研究重点。由于球鼻艏饶射的水动力非线性效果加剧^[7]，用已有的计算方法精确计算附加阻力是十分困难的。为了解决这个问题，Fujii等^[8]通过采用一些互补系数得出了一个半经验公式。Faltinsen等^[9]通过假设船舶在水面上有一个垂直面，得出了一个简化的渐近公式。日本国家海事研究所(NMRI)^[10]提出了一种改进表达式，利用实验数据修正了互补系数。

对于油轮和散货船，最常见的2种操作条件是满载和压载条件。然而，压载条件下的附加阻力的研究并不普遍。朱仁传等^[11]对VLCC附加阻力的半经验公式进行改进，并考虑了吃水对附加阻力经验公式的影响。Orihara等^[12]将压载条件下的实验数据与数值结果进行了比较，但并未提供一种适用于压载条件的数值方法。

本研究为了研究不同吃水情况下的附加阻力，考虑了4个条件：满载、压载以及满载和压载条件之间的2个吃水条件。采用计算流体力学对船舶运动进行了数值计算，分析吃水对船运动以及附加阻力的影响并与实验结果分析正确性。对短波区域的附加阻力对比分析，得到计算流体力学相对于经验公式以及势流理论更符合实验结果，CFD数值方法在不同吃水情况下具有更好的适用性。

1 船型与数值方法

以KVLCC2船为研究船型，船模缩尺比为1:80，实船与模型的主要数据如表 1所示。

1.1 吃水条件

船舶吃水根据船舶类型和作业条件的变化而变化。对于油轮和散货船，压载条件占总运行时间的近50%。对于其他船舶，吃水取决于操作条件。

因此，研究4种吃水条件下的附加阻力和船舶运动：满载、满载和压载条件之间的2种情况以及压载条件，吃水分别设定为D₁、D₂、D₃、D₄。图 2显示了各吃水条件下的吃水线。表 2总结了4种情况下KVLCC2船的艏艉吃水、排水量、重心高度。对于船艉而言，水上平面随吃水深度而显著变化。然而对附加阻力影响较大的船艏水平面形状随吃水变化不大。

在4种吃水条件下，入射波的波长为0.4 L~2.0 L，振幅为0.005 L。在短波区，选择振幅不要太高的陡度波，设定波陡为1/60。对于设计速度，所有情况下的弗劳德数均为0.142。

1.2 CFD数值模拟设置

采用国际船模试验池会议(international towing tank conference，ITTC)建议，对于存在入射波的模拟，入口边界应位于离船体1~2 L_BP的位置，而出口应位于下游3~5L_BP的位置，以避免边界墙的任何波反射。本文选择具上部边界0.5L_BP，下部边界1L_BP，侧壁边界1L_BP，入口边界1L_BP，出口边界2L_BP，并在出口处增加消波阻尼，以减少艉部长度，从而减少计算时间。

在自由表面上每个波长至少应使用80个单元。为了捕获清晰的自由表面流，如砰击现象，在船体自由表面附近，每个波长至少使用150个网格点。耐波性模拟网格与静水模拟使用的网格相比，自由表面的精细网格面积相对较小，选择0.058 5%L_BP的单元尺寸描述波浪特征。根据长短波要求，最终分别采用网格数量为2.76×10⁶和3.42×10⁶个网格。

在RAN求解中，应用分离流模型以非耦合方式(PRESTO)求解流动方程，采用二阶逆风格式离散Rans公式中的对流项。整个求解过程是根据SIMPLE算法得到。

对于波浪的模拟，利用Fluent软件中波浪模拟功能进行造波。通过设置入射波的基本条件如水深、波速、波高、波长、波头角、相位角。其中关键设置利用Fluent中的Open Channel wave进行波浪模拟。如图 3所示，其中理论值是利用正弦波子程序对速度入口进行波浪设定的结果，是波浪设定常规方法。波浪模拟功能的准确性如图 3所示，可以准确且简便地模拟入射波浪。

为了模拟真实的船舶行为，采用动态流固耦合模型，使船舶在纵摇和升沉方向自由移动。动态流固耦合模型使Rans解算器能够计算波浪作用在船体上的激振力和力矩，并求解刚体运动的控制方程，以重新定位刚体。采用一阶时间格式对N-S方程中的非定常项进行离散。

1.3 计算流体力学理论公式

使用商用CFD软件FLUENT计算规则波中的附加阻力和船舶运动。假设流体不可压缩，则连续性方程和动量方程为：

$\frac{\partial \bar{u}_{i}}{\partial x_{i}}=0 $

(1)

$\frac{\partial\left(\rho \bar{u}_{i}\right)}{\partial t}+\frac{\partial}{\partial x_{j}}\left(\rho \bar{u}_{i} \bar{u}_{j}\right)=-\frac{\partial \bar{p}}{\partial x_{i}}+\frac{\partial}{\partial x_{j}}\left(\mu \frac{\partial \bar{u}_{i}}{\partial x_{j}}-\rho \overline{u_{i}^{\prime} u_{j}^{\prime}}\right)+\rho \bar{g}_{i} $

(2)

式中:t为时间; ρ为流体密度; μ为流体的动力粘性系数; u_i和u_j为平均速度的分量; p为平均压力; x_i和x_j分别为i和j方向上的位置坐标; g_i为单位质量的重力; $\rho \overline{u_{i}^{\prime} u_{j}^{\prime}}$为雷诺应力方程中雷诺应力项属于未知量。

为使方程(1)、(2)封闭引进SSTk-ω湍流模型：

$\frac{\partial}{\partial t}(\rho k)+\frac{\partial}{\partial x_{i}}\left(\rho k \bar{u}_{i}\right)=\frac{\partial}{\partial x_{j}}\left(\mathit{\Gamma}_{k} \frac{\partial k}{\partial x_{j}}\right)+\widetilde{G}_{k}-Y_{k} $

(3)

$\frac{\partial}{\partial t}(\rho \omega)+\frac{\partial}{\partial x_{i}}\left(\rho \omega \bar{u}_{i}\right)=\frac{\partial}{\partial x_{j}}\left(\mathit{\Gamma}_{\omega} \frac{\partial k}{\partial x_{j}}\right)+\widetilde{G}_{\omega}-Y_{\omega}+D_{\omega} $

(4)

式中：Γ_k、Γ_ω为湍动能k和比耗散率ω的有效扩散系数；Y_k、Y_ω为k和ω的湍流耗散；$\widetilde{G}_{k}$为平均速度梯度引起的k的产生项；$\widetilde{G}_{\omega}$为ω的产生项；D_ω为交叉扩散项。

2 CFD模型结果分析 2.1 网格收敛性检验

研究了短波和长波情况下的网格收敛性问题，分别对船体和波浪区域选择3种不同网格数量进行计算分析, 并对结果无因次处理，其中取短波λ/L=0.5、长波λ/L=2.0作为典型情况计算。

对不同网格下的船舶附加阻力结果分析收敛性，不同网格下的船舶附加阻力系数结果如表 3所示。在短波和长波情况下，对于船舶附近网格尺寸细化比分别为1.608和1.644，并满足单调收敛条件，分别采用网格数量2.76×10⁶和3.42×10⁶个网格。根据表格结果可知网格收敛性符合计算要求。因此，本研究针对短波λ/L=0.5和长波λ/L=2.0情况下，选择了基本网格系统进行CFD模拟。但是在短波中，需要更多的单元来满足波浪模拟的收敛性，因此，与长波模拟相比，短波模拟的计算成本更高。

2.2 CFD正确性与优越性分析

不同波长条件在15.5 kn航速下对满载船舶进行CFD数值计算。分别取波长船长比为0.4~2.0取10个值，数值计算得到船舶在波浪运动中的升沉和纵摇运动的响应振幅算子，并与Oliveira等^[14]常规迎浪试验数据进行比较，如图 3所示。其中，λ是深水中的波长，L是船长。其中势流理论结果为运用势流理论，不考虑流体粘性下，势流理论软件AQWA计算不同吃水下船舶响应振幅算子的结果。

由图 4可知，对于升沉运动，CFD数值模拟结果与实验数据结果吻合度较高，仅在波长船长比为1~1.4内略微低估了共振周期周围的升沉运动。三维势流理论方法在共振周期(1.0 < λ/L < 1.5)附近放大了升沉运动响应。纵摇运动CFD数值模拟与实验数据结果在短波和中波范围内误差较小，仅在波长船长比为1.5~2.0内略微低估了共振周期周围的纵摇运动。三维势流理论方法在λ/L>1.5时高估了纵摇运动。CFD误差出现的原因可能由于CFD中数值计算采用了非惯性参考系，当出现振幅较大的运动时会导致自由表面捕捉不准确，从而出现误差。CFD计算得到附加阻力可以很好地反映实验结果，无明显误差存在。综上，CFD数值计算可以准确地计算船舶运动响应与附加阻力。

如图 5所示，三维势流理论与实验和CFD结果相比，低估了附加阻力。这是由于势流理论船艏区域网格要求网格尺寸需要小于波长的1/7，但由于船艏具有球鼻艏，网格划分困难，准确性较低。综合可知，与势流理论结果相比，CFD数值计算附加阻力更准确，但是CFD数值计算需要更多时间。

3 吃水对船舶阻力的影响 3.1 不同吃水下的船舶运动与附加阻力

图 6描述不同吃水下的船舶运动响应，对于升沉运动在短波区域运动响应随吃水减小而增大，在长波区域升沉运动响应随吃水减小而减小。对于纵摇运动4种吃水条件下的运动响应基本相同，说明纵摇运动响应对吃水条件的改变不敏感。在4种吃水条件下，观测到纵摇运动的最大运动响应几乎处于相同位置。

图 7表示不同吃水条件下CFD船舶附加阻力的结果。不同吃水条件下的附加阻力具有相似的趋势，短波区域附加阻力随波长增加而增加，当波长与船长相近时附加阻力达到最大值，长波区域附加阻力随波长增加而减少。在短波区，压载吃水D₄的附加阻力比满载吃水D₁的附加阻力大。由于在短波区，衍射是主要影响因素，而绕射成分与球鼻艏的形状有关。尽管船型未改变，但是由于吃水的变化，船艏有效的球鼻艏的形状也随之改变。与满载吃水D₁相比，D₂、D₃和压载吃水D₄球鼻艏水线面形状更为尖锐，故船艏形状改变对短波区域的附加阻力有较大的影响。

3.2 短波区域吃水对船舶附加阻力的影响

在短波区域内分析船舶附加阻力，分别使用CFD，势流理论以及采用国际船模试验池会议试航速度—功率分析规程中的(national maritime research institute, NMRI)半经验公式结果与满载与压载吃水条件下的Larsson实验数据^[14]作对比。由于Larsson实验并未进行吃水为D₂、D₃的实验，故无法进行实验结果对比。

如图 8可知在满载吃水条件下CFD，势流理论和NMRI半经验公式均与实验数据吻合，仅在短波λ/L为0.3时势流理论和半经验公式结果略小于实验结果，但CFD数值模拟结果仍然准确。压载吃水条件下CFD与实验结果存在较小误差，但势流理论和NMRI半经验公式结果在短波λ/L小于0.6时与实验结果趋势不同，存在较大的误差。

NMRI方法辐射增阻计算精度受三维势流速度势求解方法求解精度影响，NMRI方法虽然对船型水线面形状修正，但是并不完善，由于满载吃水D₁和压载吃水D₄球鼻艏水线面形状区别较大，形状修正算法本身求解精度无法满足。使用CFD进行数值模拟计算短波区域附加阻力具有更好的合理性。根据已有实验数据和CFD模拟数据进行数据分析，对已有的Takahashi^[8]绕射增阻经验公式进行有关吃水系数的修正，本文提出的公式拟合结果如图 8所示。增阻公式为辐射增阻与绕射增阻的和：

$R_{\mathrm{AW}}=R_{\mathrm{AWR}}+R_{\mathrm{AWM}} $

(5)

辐射增阻为基于Jinkine^[13]的半解析半经验公式:

$R_{\mathrm{AWM}}=4 \rho g \zeta^{2} B^{2} / L_{P P} \cdot r_{\mathrm{AW}} $

(6)

式中：r_AW=a₁a₂a₃ω^-b₁exp(b₁/d₁(1-ω^-d₁)); r_AW为辐射增阻系数；ω为无量纲化频率；a₁为主尺度效应系数；a₂为航速效应系数；a₃为新增吃水形状系数：

$a_{3}=0.7+0.31 \theta $

(7)

式中：θ为纵倾角；a₃吃水形状系数是在已有公式的基础上考虑了由于纵倾带来的艏部吃水变化，吃水形状系数会影响辐射增阻的峰值。

绕射增阻计算方法基于Takahashi^[8]的半解析半经验公式，绕射增阻为:

$\left\{\begin{array}{l} R_{\mathrm{AWR}}=0.5 \rho g \zeta^{2} B \alpha_{1}\left(k_{e} T_{M}\right) f_{1} \\ \alpha_{1}\left(k_{e} T_{M}\right)=\frac{\pi^{2} I_{1}^{2}\left(1.5 k_{e} T_{M}\right)}{\pi^{2} I_{1}^{2}\left(1.5 k_{e} T_{M}\right)+K_{1}^{2}\left(1.5 k_{e} T_{M}\right)} \end{array}\right. $

(8)

式中：α₁(k_eT_M)无量纲化绕射增阻系数，其中k_e为遭遇频率；f₁体现了航速、吃水船型特征对绕射增阻的影响：

$f_{1}=0.692\left(\frac{U_{0}}{\sqrt{T_{M} g}}\right)^{0.769}+1.81 C_{b}^{6.95} $

式中：U₀为船速；T_M为平均吃水；C_b为方形系数。

4 结论

1) CFD数值方法计算不同吃水下的升沉纵摇运动响应以及附加阻力与实验结果吻合较好。CFD数值模拟船舶运动与附加阻力具有较高的准确性。

2) 升沉运动在短波区域运动响应随吃水减小而增大，在长波区域升沉运动响应随吃水减小而减小。对于纵摇运动4种吃水条件下的运动响应基本相同，说明纵摇运动响应对吃水条件的改变不敏感。

3) 不同吃水条件下船舶附加阻力呈现相似的趋势。在短波区域压载条件下的附加阻力最大；在中等波长区域随吃水减小附加阻力最大值减小；吃水越小，附加阻力最大值对应波长越短。

4) 在船舶设计中不仅要考虑满载状态下的阻力和附加阻力，还要考虑其他吃水状态下的附加阻力。针对不同吃水条件，提出改进的波浪增阻半经验公式。