0 引　言

随着温室效应和全球变暖现象的日益突出，温室气体的排放控制受到越来越多关注。国际海事组织（IMO）推出“能效设计指数”（EEDI）对船舶节能减排进行强制约束，船舶的波浪增阻是EEDI评估公式中海洋气象因子fw计算的关键，因而在船舶型线设计中受到了广泛的重视。

学术界对船舶波浪增阻的研究已有60年的历史。预报手段上分为模型试验和数值预报等。根据控制方程的不同，数值预报方法又可以分为被分为不考虑粘性影响的势流方法和RANS方法。RANS方法需考虑流体的粘性效应和非线性效应，预报精度高；但计算量大，预报速度较慢。势流方法预报成本低、速度快，在工程设计中运用较多。势流方法经历了从二维到三维、从频域向时域发展的过程。Havelock^[1]最先通过势流方法进行波浪增阻预报，基于Froude-Krylov 假设，认为影响船体运动的主要作用力是惯性力。Maruo^[2]对该方法进行了修正，将速度势速度势分解为入射效应、绕射效应和辐射效应，通过边界元法计算获得波浪增阻。Gerritsma等^[3]从能量守恒的角度出发，认为船舶纵向运动所产生的辐射波是引起波浪增阻的原因，据此提出了波浪增阻计算的辐射能量法。Matulja等^[4]采用切片法获得船体运动响应，并对采用近场法和远场法获得波浪增阻进行了比较。段文洋等^[5]基于高速细长体理论进行绕辐射问题，获得了计算高速船阻力波浪增阻的方法。Liu等^[6]通过三维时域势流理论获得船体运动响应并通过远场法获得船舶波浪增阻。Song等^[7]采用三维时域Rankine源法计算船舶船舶在波浪中的运动响应，并分别采用近场和中场法获的了船舶的波浪增阻。

势流方法对于中长波中阻力增加的预报有一定精度，在船舶早期型线设计中获得了广泛地应用。随着船舶设计逐渐向大型化方向发展，短波中波浪增阻预报受到了越来越多的关注。航行船舶的耐波性研究表明，短波中船舶运动响应小，辐射波引起的波浪增阻较小。波浪增阻主要受到波浪绕射效应及非线性等因素的影响。高桥等^[8]在考虑船体运动引起阻力增加的同时，考虑了船首反射波引起的阻力。Faltinsen等^[9]针对短波中的绕射问题给出了一个渐近表达式，适用于短波，钝型外形和傅汝德数低于0.2的情况。Liu等^[10]基于130多种不同船型的3000多组波浪增阻数据库总结了波浪增阻预报的半经验公式。Liu等^[11]针对该公式在具有极限尺度比的船舶波浪增阻预报不准问题上进行了修正。赵博文等^[12]利用切片理论对全浪向下船舶波浪增阻进行了预报与分析。袁帅等^[13]采用RANS 方法对有限水深中船舶的波浪增阻开展了研究，并给出了池壁效应对波浪增阻预报的影响。刘树魁^[14]对各种波浪增阻预报方法的研究现状进行了总结，并对其发展历程和未来发展趋势进行了评述和展望。随着计算机技术的发展和人工智能的逐渐成熟，Mittendorf 等^[15]采用机器学习对船舶在各个浪向下的波浪增阻进行了预报。

本文针对航行船舶在波浪中的阻力增加进行研究，针对现有势流方法计算短波增阻精度的不足，提出了一种波浪增阻计算的混合方法。在中长波范围内，基于频域切片法计算获得航行船舶在波浪中的绕射势、辐射势及运动响应。进一步地，采用Salvensen法^[9,16−17]计算获得波浪增阻的各个成分，即漂移力、阻尼力和绕射力。在短波范围内，采用藤井和高桥^[8]的近似公式或者Faltinsen^[9]经验公式获得短波增阻。引入一个过渡参数实现长短波中波浪增阻预报结果平稳过渡。对Wigley3、KVLCC2船和一艘50000 t油船在不同航速、波长下的阻力增加进行数值预报，数值结果与试验值吻合良好，验证了本文方法在波浪增阻计算中的适用性。

1 理论基础 1.1 基于切片法的船舶水动力计算

船舶以定速U₀航行时，在波浪中的运动计算采用的坐标系o-xyz如图1所示。其中，坐标原点位于静止自由面上且以速度U₀沿船长方向移动，oxy平面与静水面重合，z轴竖直向上。船体在波浪激励下作六自由度运动，ξ_j (j=1,2,3)分别表示沿x轴、y轴、z轴的线位移，ξ_j (j=4,5,6)分别表示沿x轴、y轴、z轴的角位移。

假设海波为微幅规则波，则流场内非定常速度势Φ=ϕe^iωt，非定常速度势可以表示为：

$\phi=\phi_0+\phi_7+i\omega\sum\limits_{j=1}^6\xi_j\phi_j。$

(1)

式中：ω为遭遇频率；ϕ₀为入射势；ϕ₇为绕射势；ϕ_j为第j阶模态运动引起的辐射速度势，其运动幅值为ξ_j。

入射势ϕ₀可以被表达为如下形式：

$\phi_0=\frac{ig\zeta_a}{\sigma}e^{kz-ik\left(x\cos\beta+y\sin\beta\right)}。$

(2)

式中：β为波浪的浪向角；σ 为波浪的自然频率。与遭遇频率ω存在如下关系：

$\omega=\sigma-kU_0\cos\beta。$

(3)

式中：k为波数，k=σ²/g。

假设船体细长，在船速U₀为低速，遭遇频率ω较高时，速度势满足的控制方程可简化为二维形式，速度势满足的控制方程和边界条件如下：

$\left\{ \begin{array}{l}\begin{array}{cc}\left[L\right] & \dfrac{{\partial }^{2}{\phi }_{\text{j}}}{\partial {y}^{2}}+\dfrac{{\partial }^{2}{\phi }_{\text{j}}}{\partial {z}^{2}}=0,\qquad 流体域内,\end{array} \\ \begin{array}{cc}\left[F\right] & \dfrac{\partial {\phi }_{\text{j}}}{\partial z}-\dfrac{{\omega }^{2}}{g}{\phi }_{\text{j}}=0,\qquad 在z=\text{0}上,\end{array} \\ \begin{array}{cc}\left[S\right] & \dfrac{\partial {\phi }_{j}}{\partial n} = i\omega {n}_{j} + {U}_{0}{m}_{j},j = 1,\dots 6; \dfrac{\partial {\phi }_{7}}{\partial n} = - \dfrac{\partial {\phi }_{0}}{\partial n}\text{ } ,\end{array}\text{ } \\ \begin{array}{cc}\left[B\right] & \underset{y\to -\infty }{lim}\nabla {\phi }_{\text{j}}=0\text{ },\end{array} \\ \begin{array}{cc}\left[R\right] & \underset{x\to \pm \infty }{\mathrm{lim}}\left(\dfrac{\partial {\phi }_{\text{j}}}{\partial x}\mp i{k}_{0}{\phi }_{\text{j}}\right)=0\text{ }。\end{array}\text{ } \end{array}\right.$

(4)

式中：n_j（j=1，...6）为第 j 个模态的法向量，有：

$\left\{ \begin{gathered} \left( {{n_1},{n_2},{n_3}} \right) = \vec n，\\ \left( {{n_4},{n_5},{n_6}} \right) = \vec r \times \vec n。\\ \end{gathered} \right.$

(5)

m_j为航速对物面条件的影响，具体有：

$\left\{ \begin{gathered} \left( {{m_1},{m_2},{m_3}} \right) = \left( {0,0,0} \right)，\\ \left( {{m_4},{m_5},{m_6}} \right) = \left( {0,{U_0}{n_3}, - {U_0}{n_2}} \right) 。\\ \end{gathered} \right.$

(6)

绕射势和辐射势采用如下边界积分方程进行求解:

$\phi \left( P \right) = \frac{1}{{2{\text{π}} }}\int_C {\left[ {\frac{{\partial \phi \left( Q \right)}}{{\partial {{\boldsymbol{n}}_Q}}}G\left( {P,Q} \right) - \partial \phi \left( Q \right)\frac{{\partial G\left( {P,Q} \right)}}{{\partial {{\boldsymbol{n}}_Q}}}} \right]} {\mathrm{d}}l。$

(7)

式中：积分下标C表示二维切片流体的边界；G(P,Q)为二维格林函数，表示源点Q在场点P处产生的速度势；格林函数所满足的控制方程为$\nabla ^2$ G(P,Q)=−2πδ(P,Q)（δ为狄拉克函数）；n_Q为Q点处流体边界的法向向量。

求解上述边值问题后，表征辐射力的水动力系数可表示为：

${\mu _{jk}} - \frac{i}{\omega }{\lambda _{jk}} = \rho \iint_{{S_0}} {{\phi _k}{n_j}{\mathrm{d}}S}。$

(8)

根据上述速度势的线性分解，波浪扰动力可表示为入射波浪力和绕射波浪力的和：

${F_{Wj}} = F_{Wj}^I + F_{Wj}^D。$

(9)

由伯努利方程可知：

$F_{Wj}^I = - \rho i\sigma \iint_{{S_0}} {{\phi _0}{n_j}{\mathrm{d}}S}，$

(10)

$F_{Wj}^D = - \rho i\omega \iint_{{S_0}} {{\phi _7}{n_j}{\mathrm{d}}S} + \rho {U_0}\iint_{{S_0}} {\frac{{\partial {\phi _7}}}{{\partial x}}{n_j}{\mathrm{d}}S} 。$

(11)

求解上述方程获得各个运动模态的附加质量、阻尼系数以及波浪力，根据动量定理与动量矩定理，建立船舶六自由度的运动方程：

$\sum\limits_{k = 1}^6 {\left[ { - {\omega ^2}\left( {{{\boldsymbol{M}}_{jk}} + {\mu _{jk}}} \right) + i\omega {\lambda _{jk}} + {C_{jk}}} \right]{\xi _k}} = {{\boldsymbol{F}}_{Wj}}，j = 1,2 \ldots 6 。$

(12)

式中：M_jk 为广义质量矩阵；μ_jk 和 λ_jk 为附加质量系数和阻尼系数；C_jk 为回复力系数；F_j为 j 方向的波浪力。

1.2 计算波浪增阻的Salvensen法

Salvensen^[17]将航行中船舶的阻力增加表达为一阶势的形式，阻力增加由3部分组成，如下：

$R_{aw}=R^I+R^D+R_7。$

(13)

式中：${R_{aw}}$为总波浪增阻；${R^I}$为漂移力；${R^D}$为阻尼力；R₇为绕射力。

漂移力的计算方法如下：

${R^I} = \sum\limits_{j = 3,5} {\frac{1}{2}\rho k{\xi _j}\sigma \cos \beta \iint_{{S_B}} {{n_j}\phi _0^ * {\mathrm{d}}s}} 。$

(14)

式中：$\phi _0^ *$为入射速度势ϕ₀的共轭复数；S_B为船体湿表面。

阻尼力的计算方法如下：

${R^D} = \sum\limits_{j = 3,5} {\frac{i}{2}\rho k{\xi _j}\cos \beta } \iint_{{S_B}} {{\varphi _j}\frac{{\partial \varphi _0^ * }}{{\partial n}}}。$

(15)

绕射力的计算方法如下：

$R_7=\frac{i}{2}\rho k\cos\beta\iint_{S_B}^{ }\left(\varphi_7\frac{\partial}{\partial n}-\frac{\partial\varphi_7}{\partial n}\right)\varphi_0^*\mathrm{d}s。$

(16)

1.3 短波中波浪增阻的近似计算

Faltinsen^[5]短波阻力计算经验公式如下：

$\begin{split} {R_{{{aw}}}} =& \frac{1}{2}\rho g\zeta _a^2\int_L \Biggr\{ {{\sin }^2}\left( {\theta - \beta } \right) +\\ &\frac{{2\sigma {U_0}}}{g}\left[ {1 + \cos \theta \cos \left( {\theta - \beta } \right)} \right] \Biggr\} \sin \theta {\mathrm{d}}l。\end{split}$

(17)

各夹角的示意图如图2所示。

高桥等^[4]的短波增阻计算近似公式如下：

${R_{aw}} = {\alpha ^2}{R^2}\frac{1}{2}\rho g\zeta _a^2\int_{ - B/2}^{B/2} {{{\sin }^2}\beta {\mathrm{d}}y}。$

(18)

$其中，\alpha^2=\left(1+U_0\frac{c}{2}\right)^2或取\alpha^2=1+5F_n^{1\mathord{\left/\vphantom{12}\right.}2}，$

(19)

$R={\rm{\pi}}I_1\sqrt{{\rm{\pi}}^2I_1^2\left(kd\right)+K_1^2\left(kd\right)}。$

(20)

式中：B为船宽；d为船舶吃水；c²=g/k；I₁为第一类修正后贝塞尔函数；K₂为第二类修正后贝塞尔函数。

1.4 阻力增加计算的混合方法

1.2节长波阻力增加计算方法和1.3节短波阻力增加计算方法计算所得阻力增加存在差异，引入一个过渡参数实现长短波中波浪增阻预报结果平稳过渡，具体如下：

${R_{{{aw}}}} = \left( {1 - {\alpha _d}} \right){R_{{\text{aw1}}}} + {\alpha _d}{R_{{\text{aw2}}}}。$

(21)

式中：R_aw1为采用Salvensen法计算所得波浪增阻；R_aw2为短波增阻。

$\alpha_d=\frac{\mathrm{\pi}^2I_1^2\left(k_ed\right)}{\mathit{\mathrm{\pi}}^2I_1^2\left(k_ed\right)+K_1^2\left(k_ed\right)}。$

(22)

式中：I₁和K₂的定义同式(16)；k_e为遭遇频率对应的波数，k_e= ω ²/g。

2 数值计算结果

为检验上述方法的正确性，本文选取了细长型船Wigley3及2艘肥大型船（KVLCC2船和1艘万吨油船）进行阻力增加数值计算，并与试验结果和其他数值计算结果进行了比较。

2.1 Wigley3

Wigley3船的主尺度如表1所示，对其在迎浪海况、F_n=0.2、0.3航速下在波浪中的运动响应及波浪增阻进行计算，验算本文提出的波浪增阻计算混合方法的准确性。Wigley3的航速在F_n=0.2、0.3时波浪增阻计算结果如图3所示，图中横轴为波长与船长的比值λ/L，纵轴为无因次化的阻力增加，即$\dfrac{{{R_{aw}}}}{{\rho g\zeta _a^2{{{B^2}}/ L}}}$。其中${\zeta _a}$为波幅，B为型宽，L为船长。

图中“Exp”试验数据来自文献[18]，Wigley3船为细长体船，本计算为仅采用Salvensen法计算所得波浪增阻，总波浪增阻为Salvensen法中3种成分阻力之和。可知，波浪增阻预报结果与试验值吻合良好。这说明对于Wilgey3这种细长体船，仅采用Salvensen法即可较为准确地预报船舶在波浪中的阻力增加。从波浪增阻的3种成分分布可以看出，漂移力占波浪增阻的绝大部分成分。对于Wilgey3这种细长体船，绕射效应—绕射力的大小对波浪增阻的结果影响较小。影响波浪增阻的主要成分是辐射运动引起的阻力增加。

2.2 KVLCC2

KVLCC2船的主尺度如表2所示。对其在F_n=0.142航速下进行频域数值计算，浪向角为180°迎浪。波浪增阻计算结果如见图4所示，横轴为波长与船长的比值λ/L，纵轴为无因次化的阻力增加，即${{{R_{aw}}} / {( {\rho g\zeta _a^2{{{B^2}} / L}} )}}$。其中${\zeta _a}$为波幅，B为型宽，L为船长。

图中“Exp”试验数据来自文献[19]，“CFD”数值结果来自文献[20]，“Salvensen”数值结果表示仅采用Salvensen法计算波浪增阻，“Faltinsen修正”和“藤井修正”表示分别在短波增阻计算中采用Faltinsen经验公式及藤井和高桥经验公式的波浪增阻混合计算方法。从图中可以看出，仅采用Salvensen法计算波浪增阻在长波区域与试验数据及CFD方法计算结果吻合良好，但在短波区域数值计算结果偏小。这是因为在长波中，引起波浪增阻的主要成分是船体运动，即辐射效应引起的。本文采用切片法可较为准确地预报船体的运动，因而采用Salvensen法可较为准确地波浪增阻。而在短波中，船体运动较小，引起波浪增阻的主要成分是波浪反射，即绕射效应引起的波浪增阻。而切片法基于线性势流理论，无法准确描述波浪的非线性效应和粘性效应，因而采用Salvensen法预报的波浪增阻误差较大。基于此，本文提出混合波浪增阻方法，在短波中采用经验公式对波浪增阻进行预报，弥补Salvensen法从原理上预报波浪增阻的缺陷。无论是使用Faltinsen经验公式及藤井和高桥经验公式进行修正，都可以在保留长波区域波阻计算精度的同时提高短波区域的预报精度。因而能够在全波长范围内获得较为准确的波浪增阻预报结果

2.3 5万吨油船

5万吨油船主尺度见表3。计算航速分别为F_n=0.146、0.159、0.171、0.183，浪向角为180°迎浪。

图5为5万吨油船在迎浪工况下以4种不同航速航行时的阻力增加的数值结果，图中纵轴为无因次化的阻力增加，试验数据来自上海船舶运输科学研究所的试验报告。可以看出，在中长波范围内（λ/L>0.5）仅采用Salvensen法能够较为准确的预报出船舶在波浪中的阻力增加；而在短波范围内（λ/L<0.5），仅采用Salvensen法预报出船舶在波浪中的阻力增加值偏小，而采用本文提出的混合波浪增阻计算方法，无论是使用Faltinsen经验公式及藤井和高桥经验公式进行修正，获得的计算结果与试验数据吻合良好。这是因为，在中长波范围内波浪增阻的主要成分是辐射增阻（即辐射运动产生的阻力增加），而势流理论能够较为准确地预报船舶在波浪中的辐射运动。而在短波范围内，波浪增阻的主要成分是绕射增阻，受波浪非线性及粘性效应影响较大，因而需要在短波范围内对Salvensen法预报所得波浪增阻进行修正。对50000 t油船在不同航速下的阻力增加预报验证了本文提出的混合计算方法的准确性和适用性。

3 结　语

本文提出波浪增阻计算的混合方法，在中长波区域采用Salvensen法计算波浪增阻，在短波范围内采用Faltinsen经验公式或者藤井和高桥近似公式计算波浪增阻，通过一个过渡参数实现长短波中波浪增阻预报结果平稳过渡。计算分析了Wigley3、KVLCC2和50000 t油船迎浪状态下不同航速中航行的波浪阻力增加，有如下结论：

1） 对于肥大性船舶，采用Salvensen法的数值计算结果在长波区域与试验值吻合良好，但在短波区域预报结果偏小。

2） 本文提出的混合波浪增阻计算方法吸取了Salvensen法在长波区域及近似公式在短波区域波浪增阻预报的优点并规避的缺陷，实现了全波长范围内波浪增阻的高效准确预报。