0 引　言

船舶在波浪中航行会受到比在静水中航行更大的阻力，为抵抗这部分波浪增阻，一般会使用储备功率来维持波浪中航速的不变。但是为了减少温室气体的排放，在2011年MEPC的第62届会议上，IMO^[1]（International Maritime Organization）正式通过了船舶能效设计指数（Energy Efficiency Design Index，EEDI），EEDI的含义为船舶单位载重量单位航速下的二氧化碳排放量。所以，提高海上航行船舶的燃油效率至关重要。造船工程师需要不断开发高性能的船型来满足需求，包括减少波浪增阻。其中，重要的是要了解波浪增阻与船体形状之间的关系。

为了有效减小波浪增阻，许多计算方法和船首形状被提出。1980年，Faltinsen^[2] 基于直首假定提出了适用于肥大型船的短波增阻渐进修正公式。2003年Yamasaki等^[3]的研究中指出普通船首会向前反射入射波使波浪中阻力增大，以Ax-bow为代表的水线面形状为向前尖锐的船首形状，左右分割反射入射波，可以减少波浪中的阻力增加。Kuroda等^[4]用三维的壁面模型进行了迎浪规则波下的试验，从而探究速度、钝性系数与速度，以及吃水和频率对于波浪增阻的影响，最后提出修正遭遇频率和前进速度这2项公式。2014年，Kim等^[5]对KVLCC2分别布置原始船首、AX-bow和leadge bow等3种不同船首的情况进行了水池试验，AX-bow和leadge bow对波浪的反射比原始船首低，在短波中表现更为明显，对比其波浪增阻，研究发现leadge bow对波浪增阻的减小比AX-bow更明显。2019年，Lee等^[6]研究了66000 t散货船在静水和波浪条件下，船首形状对于船舶阻力的影响，通过数值计算和试验研究在迎浪条件下，尖首型船体在前肩附近的浪高比钝首船低。尖首的停滞压力区比钝首要小，使得前肩附近的压力比钝首小，在波浪中尖首型船比钝首型船总阻力减小了4.5%。2020年，Yasukawa等^[7]为研究船舶方形系数为0.87的船首形状对波浪增阻的影响，并分别在规则波和不规则波下进行了阻力试验，得到不突出弓形的直首可以减少波浪增阻。2021年，Trung-Kien Le等^[8]应用商业计算流体力学（CFD）程序，研究了规则波迎浪下船首形状对船体阻力的影响。对钝首和球首这2种船型的水动力性能和阻力进行了模拟，提出了一种新的球根状船首形状，大大减少了波浪增阻。最后给出了船舶水动力性能的CFD计算结果。顶浪中的波浪增阻研究较为丰富，但全浪向下的波浪增阻研究相对较少。Valanto等^[9]在2015年的论文中对HSVA船舶进行了全浪向的试验，发现最大的波浪增阻发生不是在顶浪方向上，而是在首斜浪方向上。本研究将分析比较5种典型的船首在首斜浪向下对于波浪增阻的影响。

在Kim等^[5]的研究中发现切片法和经验公式法总体上不如三维面源方法和CFD方法。在三维波浪增阻的计算中，Ferreira^[10]根据动量守恒原理建立了一个可以计算有航速下波浪增阻的中场方法。Pan等^[11]根据Ferreira的方法，采用时域Rankine源方法进行了波浪增阻的计算。本研究采用此三维时域Rankine源方法对船舶波浪增阻进行计算分析。

1 船型参数

本研究针对某超大型集装箱船进行研究，在保证船中、船尾、船长和船宽基本不变的情况下，对船首形状进行改变，本计算仅考虑平均吃水线下的船体形状，形成了5种不同船首，分别是前倾艏、球鼻艏、直艏、飞剪艏和斧型艏，如图1所示。

本研究以装有上述船首的超大型集装箱船为研究对象，各船模主尺度参数见表1。

将直艏船模的船首下沉深度进行变化，如图2所示，保持半进流角基本不变，直艏底部距离基线分别为13、11、9、7 m。4种直艏船模的主尺度参数：船长均为365 m，船宽均为51 m，半进流角保持不变为0.202，重心位置变化较小。其中，浪向角定义为船舶航行方向与波浪方向的夹角，所以180°为顶浪航行，0°为随浪航行。

以装有上述船首的超大型集装箱船为研究对象，开展船模吃水14.5 m，航速23 kn工况下顶浪和首斜浪下波浪增阻的计算。

2 计算方法 2.1 三维时域Rankine源法

船舶在大地参考系下的前进速度为$ {\boldsymbol{U}} $，滑移速度为${\boldsymbol{V}}$，旋转速度为${\boldsymbol{\varOmega}}$。在随船坐标系下，可以描述船舶的非定常运动。非定常位移$ {\boldsymbol{X}}\left( {{{x}},t} \right) = {\boldsymbol{T}}\left( t \right) + {\boldsymbol{R}}\left( t \right) \times {\boldsymbol{x}} $，其中，$ {\boldsymbol{x}} = \left( {x,y,z} \right) $为位置矢量，刚体平移为$ {\boldsymbol{T}} = \left( {{\xi _1},{\xi _2},{\xi _3}} \right) $，刚体旋转为$ {\boldsymbol{R}} = \left( {{\xi _4},{\xi _5},{\xi _6}} \right) $。在流体体积内无粘性、无旋转、不可压缩假设下，在此假设下，平均速度场${\boldsymbol{W}} = \left( {{\boldsymbol{U}} - {\boldsymbol{\varOmega}}x} \right)\hat i + \left( {{\boldsymbol{V}} - {\boldsymbol{\varOmega}}y} \right)\hat j$，流体速度${\boldsymbol{V}} = \nabla {\varPsi} \left( {{\boldsymbol{x}},t} \right)$。总扰动速度势$ \varPsi \left( {{\boldsymbol{x}},t} \right) $在流域中满足拉普拉斯方程和边界条件：

$ {\nabla ^2}\varPsi \left( {{\boldsymbol{x}},t} \right) = 0，$

(1)

$ {\left\{ \begin{array}{l} \nabla p - {p}_{a} = -\rho \left(\displaystyle\frac{{\mathrm{d}}}{{\mathrm{d}}t}\varPsi + \frac{1}{2}\nabla \varPsi \cdot \nabla \varPsi + gz\right) \quad 伯努利方程，\\[-8pt] \\ \dfrac{\partial \varPsi }{\partial n}=\boldsymbol{V}_{B}\cdot \boldsymbol{n} \quad 物面条件，\\ \left(\displaystyle\frac{\partial }{\partial t} + \nabla \varPsi \cdot \nabla \right)\left[z - \zeta \left(x,y,t\right)\right] = 0 \; 自由表面运动学条件，\\ \displaystyle\frac{{\mathrm{d}}\Psi }{{\mathrm{d}}t}=-g\zeta -\frac{1}{2}\nabla \varPsi \cdot \nabla \varPsi \quad 自由表面动力学条件。\end{array} \right.}$

(2)

$p$为流域中某一点的压力；${p_a}$为大气压强；$\rho $为海水密度；$g$为重力加速度，$z$为该点处的相对高度。$ {{\boldsymbol{V}}_B} = {{\boldsymbol{V}}_{Bs}} + {{\boldsymbol{V}}_{Bu}} $为船体表面流体的总速度。$ {{z}} = \zeta \left( {x,y,t} \right) $是自由表面波形位置。此外，完整的定解问题还应包括合适的水底条件，辐射条件和初始条件。三维时域Rankine源方法中将总速度势$\varPsi \left( {{{x}},t} \right)$分为定常基本流速度势${\phi _{{B}}}$、局部速度势$\phi $、记忆速度势$\psi $。

$ \varPsi \left( {{\boldsymbol{x}},t} \right) = {\phi _B}\left( {\boldsymbol{x}} \right) + \phi \left( {{\boldsymbol{x}},t} \right) + \psi \left( {{\boldsymbol{x}},t} \right)。$

(3)

对于Double-body线性化，$ {\phi _B} $满足的边界条件为：

$ \left\{ {\begin{array}{*{20}{l}} & \displaystyle {\frac{{\partial {\phi _B}}}{{\partial n}} = 0}，{P \subset {S_F}}，\\ & \displaystyle {\frac{{\partial {\phi _B}}}{{\partial n}} = {{U}} \cdot {{n}}}，{P \subset {S_B}} 。\end{array} } \right. $

(4)

式中，$ {{n}} $为湿船体表面的单位法向量，指向湿船体表面。将式(3)代入式(2)，略去高阶小量并在平均自由表面展开得到自由表面运动学条件和自由表面动力学条件：

$ \left\{ \begin{gathered} \left( {\frac{\partial}{{\partial t}} - \left( {{\boldsymbol{W}} - \nabla {\phi _B}} \right) \cdot \nabla } \right)\zeta = \frac{{{\partial^2}{\phi _B}}}{{\partial{z^2}}}\zeta + \frac{{\partial\left( {\phi + \psi } \right)}}{{\partial z}}，\\ \left( {\frac{\partial}{{\partial t}} - \left( {{\boldsymbol{W}} - \nabla {\phi _B}} \right) \cdot \nabla } \right)\left( {\phi {\text{ + }}\psi } \right) = - g\zeta + \left( {{\boldsymbol{W}} \cdot \nabla {\phi _B} - \frac{1}{2}\nabla {\phi _B}\nabla {\phi _B}} \right) 。\\ \end{gathered} \right. $

(5)

物面上的非定常部分是由$\phi $导致，所以在物面上S_B满足：

$ \frac{{\partial \phi }}{{\partial n}} = \frac{{\partial {\boldsymbol{X}}}}{{\partial t}} \cdot {\boldsymbol{n}} + \left( {{\boldsymbol{W}} - \nabla \Phi - \nabla \psi } \right) \cdot {\boldsymbol{n}} 。$

(6)

将式(6)在平均船体湿表面展开，可以得到局部速度势$\phi $满足的物面条件为：$ \left( {{n_1},{n_2},{n_3}} \right) = {\boldsymbol{n}} $；$ \left( {{n_4},{n_5},{n_6}} \right) = {\boldsymbol{x}} \times {\boldsymbol{n}} $；$ \left( {{m_1},{m_2},{m_3}} \right) = \left( {{\boldsymbol{n}} \cdot \nabla } \right)\left( {{\boldsymbol{W}} - \nabla {\phi _B}} \right) $；$ \left( {{m_4},{m_5}, {m_6}} \right) = \left( {{\boldsymbol{n}} \cdot \nabla } \right)( {{{x}} \times \left( {{\boldsymbol{W}} - \nabla {\phi _B}} \right)} 。$

$ \frac{{\partial \phi }}{{\partial n}} = \sum\limits_{j = 1}^6 {\left( {\frac{{\partial {\xi _j}}}{{\partial t}} + {\xi _j}{m_j}} \right)} 。$

(7)

线性自由表面条件为：在$z = 0$平面上，$ \phi = 0 $。将局部速度势$\phi $进一步分解: 对于每一个模态，规范化的速度势$ {N_k} $和$ {M_k} $都要满足边界条件。

$ \phi = \sum\limits_{k = 1}^6 {\left( {{N_k}\left( {\boldsymbol{x}} \right){{\dot \xi }_k}\left( t \right) + {M_k}\left( {\boldsymbol{x}} \right){\xi _k}\left( t \right)} \right)}。$

(8)

速度势的求解用Rankine源方法。P、Q为布置的源汇点，根据第二格林公式，将前述的边值问题转化为边界积分方程，可得：

$ \begin{aligned}[b] 2{\text π} \varPhi \left( P \right) - &\iint_{{S_F}\cup{S_B}} \frac{\partial \varPhi \left( Q \right)}{\partial n\left( Q \right)}\left( {P,Q} \right) + \\ & \varPhi \left( Q \right)\frac{\partial G\left( {P,Q} \right)}{\partial n{\left( Q \right)}}{\mathrm{d}}S\left( Q \right) = 0 。\end{aligned} $

(9)

式中，$ G\left( {P,Q} \right) $为Rankine源，有

$ G\left( {P,Q} \right) = \frac{1}{{\left| {PQ} \right|}} 。$

(10)

总速度势分解得到的${\phi _{{B}}}$、$\phi $、$\psi $均满足式(9)。基本速度势和局部速度势的求解是边值问题，仅需在空间上离散。求解记忆速度势需要在时间和空间上进行离散。先通过初始方程和初始条件求出初始解，将初始解代入运动学条件求出下一个时刻的波面升高，然后通过此波面升高代入动力学条件可以求得该时刻自由表面的$\psi $。在空间离散化之后，可以利用该时刻的自由表面上的$\psi $和船体表面${{\partial \psi } \mathord{\left/ {\vphantom {{\partial \psi } {\partial n}}} \right. } {\partial n}}$代入离散的积分方程中，得到该时刻自由表面上的${{\partial \psi } \mathord{\left/ {\vphantom {{\partial \psi } {\partial z}}} \right. } {\partial z}}$和船体表面$\psi $^[12]。

2.2 中场方法

中场方法指在距离船体较近的控制面上积分得到波浪增阻。基于三维时域Rankine源方法，根据由控制面、船体表面以及控制面和船体表面之间的自由表面所包围的流体动量关系，可以得到如下波浪增阻的计算公式：

$ \begin{split} &{{{F}}^{\left( 2 \right)}} = - \rho \iint_{{S_0}} {\left[ { - \frac{{\left( {\nabla \varphi \nabla \varphi } \right){\boldsymbol{n}}}}{2} + \nabla \varphi \frac{{\partial \varphi }}{{\partial n}}} \right]dS} - \frac{{\rho g}}{2}\oint_{{\Gamma _0}} {{\zeta ^2}{\boldsymbol{n}}{\mathrm{d}}l} - \\& \rho \oint_{{\varGamma _0}} {\left[ {\nabla \varphi \left( {\nabla {\phi _B} \cdot {\boldsymbol{n}} - {\boldsymbol{W}} \cdot {\boldsymbol{n}}} \right) + \nabla {\phi _B}\left( {\nabla \varphi \cdot {\boldsymbol{n}}} \right)} \right]\zeta {\mathrm{d}}l}。\\[-3pt] \end{split} $

(11)

式中：$ \varphi $为入射势、局部速度势和记忆速度势一阶的总和；${S_0}$为控制面；${\varGamma _0}$为控制面和平均自由表面的交线。该方法的详细介绍见文献[11]。一阶的波面升高$ \zeta $如下式:

$ \zeta = - \frac{1}{g}\left( {\frac{{\partial \varphi }}{{\partial t}} - \left( {W - \nabla {\phi _B}} \right) \cdot \nabla \varphi } \right) 。$

(12)

入射波波幅为$ \zeta _a^{} $，将式(11)计算得到的规则波中波浪增阻传递函数值无因次化得到：

$ {K_{aw}} = \frac{{{F^{\left( 2 \right)}}}}{{\rho g\zeta _a^2\left( {{{{B^2}} \mathord{\left/ {\vphantom {{{B^2}} L}} \right. } L}} \right)}}。$

(13)

3 计算结果 3.1 5种船首对波浪增阻的影响

利用三维时域Rankine源方法数值计算所得到的5种不同典型船首的船模在规则波顶浪航行时，自由面波形如图3左侧图所示。船模首部绕射产生自由液面波形图如图3右侧图所示。

对比图3左侧船舶周围自由液面波形图：直艏和斧型艏船身周围的波浪起伏较其他船模较小，船首处波浪更容易与船首分离，自由液面幅值变化更小。斧型艏水线面较小，船舶对于流体的绕射显著减小，所以类似的船首形状更容易穿过波浪，不对波浪的幅值带来较大影响。从图3(b)可以看到，球鼻艏周围波浪的幅值变化区域较其他船首更大，说明由于球鼻艏前伸的球体直接将船体周围的流体绕射至更远的范围，球鼻艏兴波对自由液面直接产生了影响。

从图3(a)右侧首部绕射波浪场对比图可以得到，前倾艏对波浪的绕射影响较大，其较丰满的首部将波浪在端部堆积抬升，而肩部区域流体减少，说明前倾艏能够产生明显的兴波和波浪叠加垂荡幅值更大的波浪。对比图3(b)，球鼻艏对波浪也有一定的堆积作用，但由于球鼻艏体积较小且在水面以下，球鼻艏半进流角较前倾首要小的多，对波浪的反射作用要相对较小，在肩部处球鼻艏兴起的波浪抵消一部分原有波浪。直艏绕射波浪场云图如图3(c)，由于直艏的前端下沉深度不够，没有飞剪型和斧型艏类似的尖锐型形状，所以在船首端部对波浪依然存在堆积和反射现象，导致艏部端部波浪升高，也同时升高了船舶肩部的波浪。对比图3(d)和(e)右侧飞剪艏和斧型艏绕射波浪场，斧型艏的半进流角小于飞剪首，斧型艏对周围流场的影响更小，破开波浪的效果要比飞剪艏更好。

图4为相对于图3自由液面变化时刻船模所受的流体动压力图。结合图3(a)和图4(a)可以得到，前倾艏将波浪堆积在肩部，使得船舶肩部小部分区域压力值明显增大，而后逐步减小。同理，在球鼻艏的上方也由于流体堆积而负压幅值明显增大的，但球鼻艏兴波与波浪进行抵消，其肩部的压力没有明显增大。由图4(b)可见，球鼻艏带来的兴波作用使得船舶肩部及首部所受到的压力变化较为平缓，兴波作用范围较其他船首更长，但总体较直艏带来的压力更大，且图4(c)中，直艏所受到的表面压力较为平缓，没有明显的压力下降或升高，使得船身受到的压力保持在较低的范围内波动。图4(d)飞剪型船首相较于前倾艏和球鼻艏而言，在减小船舶肩部压力的同时，压力变化更加平缓。由图4(e)，斧型艏受到的压力较其他4种船模首部顶端受到的负压幅值明显偏大，原因是斧型艏的下沉深度更深，导致其受到的水压更大，但是斧型艏肩部压力变小速率快，说明斧型艏的半进流角较小，使得来流对肩部的冲击压力小。

图5为5种船模在3种浪向下的无因次波浪增阻系数对比图。可知，5种船模的无因次波浪增阻系数变化趋势一致，在短波区（$\lambda /L$为0～0.5），$ {K_{aw}} $随$\lambda /L$的增大先增大后减小，呈现一个小峰值；在长波区（$\lambda /L$为0.5～2.0），$ {K_{aw}} $在$\lambda /L$为1.0左右达到峰值，然后在$\lambda /L$为2.0左右逐步减小至0。其中，前倾艏的峰值较其他4种船模更大。球鼻艏和直艏计算所得的波浪增阻系数峰值较前倾艏和飞剪艏向高频短波区移动。而斧型艏计算所得的波浪增阻系数峰值较前倾艏和飞剪艏向低频长波区移动。在$\lambda /L$为1.0附近，船模波浪增阻由高到低：前倾艏>飞剪艏>球鼻艏>直艏>斧型艏，该顺序与船模半进流角减小顺序相同。

图5(b)为首斜浪（150°）下5种船模的$ {K_{aw}} $对比图。对比迎浪情况下，首斜浪下$ {K_{aw}} $的峰值向高频短波区偏移。在短波区，斜浪下由于船舶迎浪面积的增大使得绕射增阻大于迎浪下的增阻。斧型艏在斜浪下的波浪增阻系数高于迎浪下波浪增阻系数，斧型艏的半进流角过小导致船首舷侧曲率较小，其等效平板的面积更大，对斜浪的反射作用更强。在长波区峰值处，5种船首对波浪增阻的影响与迎浪下规律一致。该球鼻艏在180°和150°浪向角下，相对于前倾艏分别能起到16.7%和12.7%左右的减小波浪增阻效率。图5(c)指浪向角为120°下5种船模的$ {K_{aw}} $对比图。在长波区，5种船首船模的无因次波浪增阻系数趋于一致，说明在120°浪向角时，船首变化在长波区对于船舶波浪增阻的影响较小。在短波区，球鼻艏船模的$ {K_{aw}} $明显高于其他船模，说明该球鼻艏无法在短波区内良好减阻，甚至增大了阻力。

3.2 直首下沉深度对波浪增阻的影响

根据上述结论，直艏和斧型艏由于其半进流角比其余3种船首更小而具有更为优良的减阻性能，但斧型艏并没有因为其半进流角比直艏更小而展示出比直艏更强的减阻效果，且斧型艏对船舶利用空间影响较大，所以在不改变进流角下，改变直艏下沉深度，探究直艏到斧型艏的参数变化对于船舶波浪增阻系数的影响。

直艏船模根据底部距离基线分别为13、11、9、7 m对应4种船模。4种船模在180°、150°和120°浪向角下的无因次波浪增阻系数对比如图6所示。

浪向角为180°，在峰值处，$ {K_{aw}} $随着船首距离基线的长度的减小而减小，即直艏的长度越长（距离基线的距离越短），无因次波浪增阻的峰值越小，浪向角为150°时也基本满足此规律。距离基线7 m的船首比距离基线13 m的船首在峰值处$ {K_{aw}} $减小了10%。结合图3，可以说明在半进流角不变的情况下，船首下沉深度越深对波浪的堆积作用越小，从而对船舶运动的影响也比较小。

波浪增阻系数在短波区中与峰值处呈现相反的变化趋势，直艏的长度越长，$ {K_{aw}} $的峰值越大，结合图4可以说明，下沉深度越深，船首的湿表面积越大，船首底部受到的水压力也越大。与图5中各个浪向下斧型艏在短波区的无因次波浪增阻系数明显大于直艏波浪增阻的现象一致。

4 结　语

本文基于三维时域Rankine源理论对装有5种不同船首的船模在迎浪和首斜浪下的波浪增阻进行计算。探讨了船首形状对于波浪增阻影响机理，同时，本文研究了直艏长度对于船舶波浪增阻的影响，并对于本研究所采用的船型得出以下结论：

1） 对比首斜浪和迎浪下的无因次波浪增阻系数，斜浪下由船舶迎浪面积的增大引起的绕射增阻明显要大于迎浪下的对应增阻，所以首斜浪下的波浪增阻可能会超过顶浪下的波浪增阻。船首的变化对波浪增阻的影响随着浪向角的减小（从迎浪趋向于横浪）不断减小，120°浪向时，船首形状对于波浪增阻的影响已不到8%。

2） 球鼻艏通过兴波可以在峰值处减小船舶的波浪增阻，相比前倾艏可以减小20%，但是在短波区球鼻艏的兴波可能会增大船舶波浪增阻。

3） 船舶波浪增阻基本随着半进流角的增大而增大。在波浪增阻峰值对应频率附近，船模波浪增阻由高到低为前倾艏>飞剪艏>球鼻艏>直艏>斧型艏，该顺序与船模半进流角减小顺序相同。

4） 在半进流角不变的情况下，直艏的下沉深度越深，船舶在峰值处的波浪增阻越小，本文中距离基线7 m的船首比距离基线13 m的船首在峰值处波浪增阻系数减小了10%。但在短波区，由于等效迎浪面积增大，船舶波浪增阻系数会随着直首下沉深度的增加而增加。