0.   引　言

随着全球气候变暖，北极航线的开辟成为可能^[1]。北极航线具有天然的地缘优势，其航行距离短，且不必经过繁忙的苏伊士运河^[1-3]。但北极航线中存在大量的浮冰，这些浮冰的存在会显著增加船舶阻力^{[1, 4]}，进而增加船舶航行的时间和能耗。

船型优化是降低浮冰对船舶航行性能影响的一种有效解决方法^[5-11]，为此，有学者开展了一系列有关冰阻力船型优化的研究。刘经京等^[5]基于《芬兰−瑞典冰级规范》的碎冰阻力计算公式，并结合Shipflow仿真软件，对极地邮轮的静水阻力和碎冰阻力进行了优化，并采用CFD & DEM方法考虑了碎冰对螺旋桨性能的影响；刘瀛昊等^[6]采用Holtorp阻力计算公式，并结合船级社对破冰船的破冰功率和推进功率的要求，对某破冰船的型线进行了优化设计；段菲等^[7]基于CFD方法和冰阻力规范公式，采用NSGA-II算法对极地油轮进行了优化；王志鹏等^[8]基于CFD和代理模型（SM）对极地油船的静水阻力进行了优化，随后采用离散元黏接模型对船舶的破冰能力予以了验证；毛方云^[9]基于LS-Dyna商用软件，分析讨论了不同船首形式破冰船的破冰能力；陈超等^[10]基于Lindquist经验公式和势流方法，对极地科考破冰船进行了优化，随后又采用目标加权和NSGA-II算法进行了优化求解；Lu等^[11]基于Lindquist经验公式对极地船舶进行了优化，并采用DEM方法对优化后船型的破冰性能予以了验证。然而，目前有关冰阻力船型优化的研究均是基于经验和规范公式来优化船舶和层冰（level ice）碰撞时的冰阻力性能，针对碎冰环境（pancake ice）下的船型优化研究相对匮乏，但有关碎冰的优化研究又具有非常重要的价值^{[1-2, 4]}。碎冰是由大块浮冰破碎后形成的小块冰，在波浪冲刷和浮冰碰撞的作用下，通常呈现圆形^[2]。根据学者们的预测，这种碎冰环境将成为未来北极航线最普遍的航行环境，并且在南极也将呈现主导地位^[1]。

现有的船型优化方法应用于碎冰环境下时存在一定的局限性。首先，基于经验公式的优化方法无法考虑碎冰尺寸和分布的随机性，而已有研究表明，这些因素对冰阻力有着显著影响^[12]。其次，现有方法无法考虑流场对于碎冰的影响，而流场会改变碎冰的运动和速度，从而影响船舶的受力，尤其是在碎冰比较小的情况下^[1]。

针对上述不足，本文拟将CFD & DEM方法引入船型优化体系中，并结合CFD方法提出一种能够精确优化碎冰阻力和静水阻力的多目标船型优化方法。此外，针对CFD & DEM方法无法自主生成随机冰场^[1]以及仿真成本过高^[1-2]的问题，本文拟提出一种高效生成高浓度冰场的算法，以模拟计算真实冰场下的碎冰阻力；并通过集成学习中的XGBoost模型构建代理模型，以降低对数据量的需求，进而降低仿真成本。最后，以KCS集装箱船作为优化示例，验证所提方法的有效性。

1.   方法论

1.1   CFD & DEM方法

数值仿真是船型优化中的关键环节之一，其精度决定了优化结果的质量^[13-14]。CFD & DEM方法相比经验公式可以考虑流场、随机冰场和复杂船体形状对碎冰运动的影响，从而更加精确地计算冰阻力。CFD & DEM方法是通过在CFD方法的基础上引入DEM模型来实现在流场中模拟固体的碰撞^[1]。下面分别予以介绍。

在CFD方法中，船体周围的黏性流场采用雷诺平均纳维−斯托克斯（Reynolds-averaged Navier-Stokes，RANS）方法求解，并使用剪切应力传输（shear stress transport，SST）湍流模型来闭合RANS方程，该模型非常适合模拟高雷诺数下受壁面限制的层流和湍流。然后，采用欧拉多相流（Eulerian multiphase model，EMP）模型模拟多相流体，并选择体积法（volume of fluid，VOF）计算多相流体交界面的运动和变形^[3]。

CFD & DEM方法是通过在CFD的基础上建立拉格朗日多相模型，将固体引入流体域中，并通过DEM模型生成大量的浮冰。由于双向耦合机制计算成本高且不会显著影响冰阻力，因此本文选择单向耦合机制来预测冰阻力。同时，通过线性弹簧接触模型建立碰撞力模型，其中弹簧负责响应的弹性部分，而阻尼器则负责碰撞过程中的能量耗散^[3]，其原理如图1所示。

图  1  线性弹簧接触模型

Figure  1.  Linear spring contact model

下载: 全尺寸图片 幻灯片

将数值仿真结果与Guo等^[4]的试验结果进行对比以验证数值方法的准确性。参照实验，采用同样缩比的船型，船体模型主要参数如表1所示，参考实验仅放开纵摇和横摇这2个自由度。

表  1  KCS船主要参数

Table  1.  Main parameters of KCS ship

参数	原船型	模型船
船长/m	230.000 0	4.367 1
船宽/m	32.200 0	0.611 4
吃水/m	10.800 0	0.205 1
排水体积/m3	52 030.000 0	0.356 2
缩尺比	1.000 0	52.667 0

下载: 导出CSV 

| 显示表格

数值仿真中，计算域和网格的设置如图2～图4所示。按照ITTC和文献[1]的建议进行网格加密，并对兴波区域、水表面区域和近船区域进行重点加密。船体外壳壁面的棱柱层层数设置为6层，过渡率为1.4，保证船体附近y⁺的取值为30～60，同时保证棱柱层最外层和附近网格过渡的合理性。在压力出口处和两侧设置阻尼消波。在计算冰阻之前，采用CFD方法初始化稳定的流场，因此，CFD仿真和CFD & DEM仿真均采用同等的网格数目和时间步长，网格数量为162×10⁴，网格无关性验证如表2所示，时间步参考文献[3]，设为0.02 s。碎冰参数沿则用文献[1]和文献[3]中的设置。数值仿真结果与实验的对照验证如图5所示，其中R_ice为冰阻力，R_t为静水阻力。

图  2  计算域设置

Figure  2.  Computational domain setting

下载: 全尺寸图片 幻灯片

图  3  整体网格加密展示

Figure  3.  Overall mesh refinement scheme presentation

下载: 全尺寸图片 幻灯片

图  4  船侧和船首网格加密展示

Figure  4.  Presentation of mesh refinement at the ship's side and bow

下载: 全尺寸图片 幻灯片

表  2  网格无关性验证结果

Table  2.  Grid independence verification results

网格密度	网格数	碎冰阻力误差/%	静水阻力误差/%
低	87.20×104	76.67	196.93
中	167.47×104	8.46	2.74
高	340.14×104	2.35	2.79

下载: 导出CSV 

| 显示表格

图  5  数值仿真结果验证

Figure  5.  Validation of numerical simulation results

下载: 全尺寸图片 幻灯片

1.2   HMIGA算法

在计算冰阻力之前，还需预先生成满足要求的碎冰场^[1]。在实际中，浮冰的位置和尺寸都是随机的，并据文献[1, 4]所述，浮冰的尺寸满足对数正态分布，且不存在重叠。浮冰的位置和尺寸决定了船舶阻力，为了更加精确地模拟船舶的冰阻力，必须在数值仿真之前生成符合实际的冰场。

现有的生成浮冰场的空间分割算法和不断循环试错的方法有难以生成指定分布的浮冰和高浓度冰场效率低下等问题，因此，本文采用Huang等^[1]的研究方法，即基于元启发式优化算法，提出一种混合循环试错局部搜索策略的、专用于生成浮冰场的混合多岛遗传算法（hybrid multi-island genetic algorithm，HMIGA）。

HMIGA的核心原理是在多岛遗传算法（MIGA）全局寻优的框架下加入不断循环试错的局部寻优方法，全局搜索能力和局部搜索能力的结合能有效提升算法在面对高维优化问题时的收敛质量。该算法采用实数编码将浮冰位置转化成个体，然后结合更适合实数编码的SBX（simulated binary crossover）交叉和多项式变异策略产生新个体，并采用精英保留策略保证算法的稳定收敛。在局部搜索中，HMIGA算法对循环试错方法进行了改进。首先，对大型碎冰进行位置调整，然后，寻找中、小型碎冰的合适位置，这是因为移动的大型碎冰会导致更多的浮冰重叠，因此在优化过程中不要轻易变动大小碎冰。更详细的算法实施流程如图6所示。图中：C为浮冰的浓度；R为满足对数正态分布的浮冰半径列向量；P_t为第t代的种群；F_t为第t代的适应度；R_i为单个种群中第i个浮冰的半径大小；d_圆心为浮冰i，j的圆心距离；t_max为最大迭代次数。

图  6  HMIGA算法流程图

Figure  6.  Flow chart of the HMIGA algorithm

下载: 全尺寸图片 幻灯片

采用HMIGA算法生成的浮冰场如图7所示。从图7中可以看出，浮冰场的质量更好，浮冰之间没有出现重叠，且尺寸满足对数正态分布。

图  7  基于HMIGA算法生成的浮冰场

Figure  7.  Ice floe generated by HMIGA algorithm

下载: 全尺寸图片 幻灯片

1.3   XGBoost 算法

集成学习方法的优势在于将多个弱学习器集成为一个强学习器，从而获得更好的性能^[15]。XGBoost是集成学习中具有代表性的一种算法，被广泛运用于Kaggle竞赛和工业应用。

XGBoost算法的核心是目标函数的推导和优化。目标函数的推导可以分为3个步骤：首先，针对损失函数进行二阶泰勒展开，并去除常数项；其次，将正则化项予以展开，去除常数项；最后，合并一次项次数和二次项系数，并将所有训练样本按照叶子节点进行分组。最终的目标函数如式（1）所示。

式中：$y$代表自定义控制叶子节点参数；$\lambda$为L₂正则项的参数；${{\boldsymbol{w}}_j}$为第j棵树的权值向量；$G_{j}=\sum_{i \in I_j} g_{i}$，为叶子节点$j$所包含样本的一阶偏导数的累加和，其中${{I}_{j}}$为所有样本的集合，${{g}_{i}}$为一阶梯度统计量；$H_{j}=\sum_{i \in I_j} h_{i}$，为叶子节点$j$所包含样本的二阶偏导数的累加和，其中${h_i}$为二阶梯度统计量。以上两者均为常数项。因此，此时目标函数中的变量只有第t棵树的权值向量w，优化式（1），即可获得最佳树模型。

XGBoost被认为是在分类和回归上都拥有超高性能的先进评估器^[16]，凭借其优秀的性能，可以实现仅在有限的数据集即可以达到良好的拟合能力和泛化性能，从而降低工程实际运用中对大量数据集的需求。加州房价数据集是被广泛用于机器学习回归分析验证测试中的一个经典数据库，将XGBoost模型与工程中运用广泛、比较具有代表性的ANN模型在加州房价数据集上进行训练对比，以验证其性能。为保证对比公平，两者均采用网格搜索法调整超参数，训练效果如表3所示。可以看出，XGBoost模型在同等数据量的情况下其性能显著优于ANN模型，其拟合指标$R^{2}$在训练集上为0.987 0，在测试集上为0.842 0，证明其泛化性能良好。

表  3  XGBoost模型和ANN模型的性能对比

Table  3.  Performance comparison between XGBoost model and ANN model

模型	$R^{2}$
	训练集	测试集
ANN	0.698 4	0.637 1
XGBoost	0.987 0	0.842 0

下载: 导出CSV 

| 显示表格

2.   优化流程及关键方法

2.1   整体优化流程

船型优化整体流程如图8所示。首先，基于循环最优拉丁超立方采样（circulation optimal Latin hypercube sampling，COLHS）方法，针对选定的设计变量进行采样；接着，采用参数化建模技术，根据采样结果自动建立目标船型，随之将目标船型放入商业求解器STAR-CCM+中分别求解其阻力性能和冰阻力性能；然后，以仿真结果建立数据库并用于训练XGBoost代理模型；最后，将代理模型作为优化算法中的适应度函数，以NSGA-III算法作为优化算法进行寻优得到Pareto解集，采用赋权加和方式处理Pareto解集得到最终优化后的船型。

图  8  船型优化整体流程

Figure  8.  The overall process of ship hull optimization method

下载: 全尺寸图片 幻灯片

优化流程中的代理模型和数值仿真方法已在1.1和1.2中详细介绍，本节将介绍优化流程中的其他关键技术：采样方法、参数化建模技术和多目标优化算法。

2.2   COLHS采样方法

样本点的采样数目和均匀性对代理模型的拟合精度及泛化能力至关重要^[13]。传统的最优拉丁超立方采样（optimal Latin hypercube sampling，OLHS）方法采用强制分层的策略和CL₂差异准则保证了采样点的均匀分布性^[13]。然而，采用拉丁超立方采样（LHS）方法必须预先确定采样点数目，因采样过多会造成资源浪费，过少则会导致代理模型训练不充分、精度不足，而追加新的采样点又会导致新、旧采样点间的重叠和冲突，完全重新采样虽然可以解决这些问题，但会造成已有数据集被全部浪费。

为解决上述问题，提出一种简单、实用的采样方法，即COLHS。该方法是基于OLHS分层思想^[13]的一种改进方法，其按照某一维度将相邻的采样点归为一组。通过这种方式，将整个数据集划分为多个部分，并在每个部分内随机抽取一个点，从而将整体数据集划分为多个子集，其分层抽取的原理如图9（a）所示。由于子集是基于均匀分布的初始采样结果划分而成，故任意子集的集合仍能保持均匀分布的特性。这使得在追加数据集时，新、旧数据集不会存在冲突和重叠问题，且在划分训练集和测试集时更具代表性，从而提升训练结果的泛化能力。基于COLHS方法对设计变量进行采样，随机抽取两维的采样结果如图9（b）所示。

图  9  设计变量采样结果

Figure  9.  Sampling results of design variables

下载: 全尺寸图片 幻灯片

在整体优化流程中，COLHS方法的实施步骤如下：

1） 采用LHS进行初始采样，并根据CL₂差异准则优化采样结果，提升其均匀性。

2） 按照COLHS的原理，将优化后的采样结果划分为多个子集，并存储为“.csv”格式文件。

3） 先对一定初始数目子集进行数值仿真，放入代理模型中进行训练，若代理模型满足精度要求，则停止仿真和训练，若不满足，则继续对新的子集进行仿真和训练，直至满足精度要求。

从以上流程可以看出，采用COLHS方法时不需要提前确定采样数量，可以根据拟合情况逐步对子集进行仿真与训练，从而避免不必要的资源浪费或重新训练，有助于提升整体优化效率。

2.3   参数化建模方法

参数化建模方法是船型优化方法中必不可少的一环，可避免繁琐、低效的人工建模过程。NFFD（non-uniform free-form deformation）方法^[17]具有高灵活性、连续性和良好的局部控制能力，能够以少量的控制点实现复杂形状的变形，并且在形状优化过程中不会出现突变与断裂，因此被广泛应用于船舶、航天和汽车等领域。

NFFD方法的核心是以式（2）建立控制体与实体之间的关系，然后再通过控制体变形带动实体变形。

式中：$D(s,t,u)$为物体上变形后的点的坐标；$P_{i, j, k}$为控制顶点；$W_{i, j, k}$为权值因子； $B_{i, l}(s)$,$B_{j, m}(t)$和$B_{k, n}(u)$分别为$l, m, n$次规范B样条基函数，其中$s, t, u$为局部坐标下的坐标；$p, q, r$为节点矢量的个数。$B_{i, k}(u)$是由节点向量U和$de \;Boor-Cox$递推公式计算得出，其中节点向量U采用准均匀的节点向量。

式（2）中的NURBS基函数确保了NFFD方法的局部控制能力。以图10所示的曲线来说明该问题。由图可见，在移动单一控制点后，基于NURBS基函数的NURBS曲线相比基于Bernstein基函数的Bezier曲线对其他区域变形影响更小，变形仅发生在变化控制点附近，对整体曲线形状变化影响很小。在船型优化过程中，希望能够对目标区域进行精确控制而不影响整体形状，而基于NURBS基函数的NFFD方法能够针对特定区域进行细致的调整来实现更精细化的建模，因此，本文采用NFFD方法作为参数化建模方法。

图  10  NURBS基函数的局部控制能力

Figure  10.  Local control capability of NURBS basis functions

下载: 全尺寸图片 幻灯片

NURBS基函数的引入虽然增强了FFD方法的性能，但却使得初始的局部坐标计算方法^[17]不再适用。而采用MIGA则可以解决该问题，即以最小化全局坐标(x, y, z)和局部坐标 (s, t, u)代入式（2）反算后的全局坐标间差距为优化目标，以局部坐标(s, t, u)为优化变量，最终可获得控制体坐标系下的精确局部坐标。获取局部坐标后，通过变化控制点，再根据局部坐标以及实体和控制体之间的映射关系进行计算，即可获取变化后的船型。

基于NFFD方法原理，采用MATLAB编程语言开发了参数化建模平台。该平台能够读取“.csv”文件中指定控制点的位移，然后生成大批量变形船型。所有船型相关数据以“.stl”格式存储，便于直接导入数值仿真软件中进行后续的分析，提升了不同模块间数据交互的便捷程度。采用该平台进行参数化建模的过程和效果如图11所示。

图  11  基于NFFD方法的参数化建模过程与效果

Figure  11.  Parametric modeling process and effect based on NFFD method

下载: 全尺寸图片 幻灯片

2.4   NSGA-III 优化算法

考虑碎冰阻力和静水阻力的船型优化问题是一个典型的多目标优化问题，因此，需要采用多目标优化算法进行求解。而在众多多目标优化算法中，NSGA-III算法因其Pareto排序效率高、种群多样性大^[18]而被选做本文的优化算法。

NSGA-III算法的关键在于ENS（efficient nondominated sort）排序操作和基于小生境的选择操作。相较于传统的全部个体两两对比然后逐级排序的操作，NSGA-III中的ENS排序策略仅需和已比较过的个体进行比较即可，因此计算效率更高，尤其是在面对高维多目标优化问题时。其次是在选择个体时，基于小生境和参考点的选择策略相对于基于拥挤度的选择策略，选择的个体更加均匀，进而保证了求解质量。

3.   优化示例

3.1   优化问题定义

本文中，固定航速（v = 1 m/s）和浮冰浓度（C = 70%）下的静水阻力R_c与碎冰阻力R_i被选为优化目标，并将船首变形区域的体积变化幅度$\Delta V_{\mathrm{amp}}$作为限制条件。船型优化问题的定义如式（3）所示：

在选取优化变量时，NFFD方法的一个便捷之处在于可以直接将控制点位移作为设计变量，被选作设计变量的控制点位置如图12所示，其移动方向和移动范围如表4所示。其中， x₁，x₂和x₈控制球鼻艏形状， x₃控制船首水线区域的流线型形状。这些设计变量中，仅 x₃和 x₄为单个控制点，其余均为多个控制点的组合。

图  12  设计变量说明

Figure  12.  Definition of design variables

下载: 全尺寸图片 幻灯片

表  4  设计变量及其变形范围

Table  4.  Design variables and their range of deformation

变形区域	设计变量	移动方向	变形范围/m
球鼻艏区域	$x_{1}$	z	[−0.09, 0.09]
	$x_{2}$	x	[−0.15, 0.15]
	$x_{8}$	$y$	[−0.10, 0.10]
水线截面区域	$x_{3} ～ x_{6}$	y	[−0.10, 0.10]
艏柱角	$x_{7}$	x	[−0.15, 0.15]

下载: 导出CSV 

| 显示表格

3.2   数据库与代理模型

确定了优化问题后，就需根据设计变量的取值范围依据COLHS方法进行采样，并将采样结果放入NFFD建模平台中生成目标船型，然后通过Java宏方式自动读取“.stl”格式的建模结果，最终计算100个目标船型的碎冰阻力和静水阻力，存储为“.csv”格式并建立数据库用于训练代理模型。

基于代理模型取代成本高昂的数值仿真过程是提升优化效率的一个重要环节^[13]。以设计变量作为输入，碎冰阻力和静水阻力作为输出，由此，代理模型的输入为8维，输出为2维。采用XGBoost作为代理模型并进行训练。将均方根误差（root mean square error，RMSE）作为XGBoost的损失函数，按照8:2的比例划分训练集和测试集，并添加L₁正则化项和L₂正则化项限制树模型的复杂度，然后采用网格搜索方法寻找最佳的超参数配置方案。XGBoost模型的拟合效果如图13（a）所示，其中蓝色为测试集数据，橘色表示训练集数据，测试集数据的拟合指标$R^{2}$ = 0.917 6。采用同样的数据集对ANN（artificial neural network），CNN（convolutional neural network）和RNN（ recurrent neural network）模型进行训练以对比其拟合能力，为保证对比的公平性，其他模型也采用同样的损失函数、数据划分方式和超参寻优方式。ANN，CNN和RNN的测试集的拟合指标$R^{2}$分别为0.861 9，0.721 7和0.664 1，图13展示了不同模型拟合效果的对比结果。从中可看出，在同等训练数据下，仅XGBoost模型达到了精度要求，能够支持下一步的优化，而其他模型均未达到相应的精度要求，需进一步训练更多的数据量，这也说明基于XGBoost的船型优化方法所需的数据量更少，优化效率更高。

图  13  不同代理模型拟合效果对比

Figure  13.  Comparison of regression performance of different surrogate models

下载: 全尺寸图片 幻灯片

3.3   船型优化模块

采用结合XGBoost模型的NSGA-III算法执行优化求解。具体地，NGGA-III算法中的适应度函数被替换为训练好的XGBoost模型以取代成本高昂的数值仿真过程，算法中的种群数目设定为200，变异概率为0.15，交叉概率为0.70。执行100次迭代后，优化结果如图14所示，图中蓝色点集表示多目标优化问题的Pareto前沿解。从图中可以看出，种群中的所有个体都收敛于前沿解，表明算法的收敛性能良好。

图  14  船型优化Pareto解集

Figure  14.  Pareto solution set for hull form optimization

下载: 全尺寸图片 幻灯片

Pareto前沿解意味着在解集中不存在所有目标上均有优于他们的解，从而为设计者提供了可供权衡选择的船型方案。获得前沿解后，可通过归一化后加权赋值的方法^[14]确定最终方案。本文给出了一个优化示例，其中碎冰阻力R_i和静水阻力R_c的权重相等，通过该方式选择的相对最优船型如图14中的橙色点所示。

为验证优化结果的可靠性，将优化后的船型代入数值仿真软件中进行测试。由优化算法得到的船型碎冰阻力R_i = 17.95 N，静水阻力R_c = 6.71 N，而采用数值仿真得到的结果则为R_i = 18.15 N，R_c = 6.71 N。结果表明，优化算法结果与数值仿真结果之间的差异极小，进一步验证了优化结果的稳定性与可靠性。

3.4   优化效果与讨论

经优化方法求解后，船型设计变量的变化如表5所示。图15展示了优化前、后船型的几何形状对比。优化后的船型变化主要体现在3个方面：一是球鼻艏变得更加尖锐，并且上端有所下移；二是水线处的艏倾角减小；三是水线形状扩张，前端和后端略有增加，中间部分显著增宽。具体的船型尺寸变化详见表5。

表  5  优化后的船型参数

Table  5.  Optimized hull form parametes

参数	数值		参数	数值
x1/mm	−82.51		x5/mm	98.99
x2/mm	−26.60		x6/mm	32.02
x3/mm	−82.51		x7/mm	−148.50
x4/mm	33.56		x8/mm	−99.00

下载: 导出CSV 

| 显示表格

图  15  原始船型和优化船型对比

Figure  15.  Comparison of original and optimized hull forms

下载: 全尺寸图片 幻灯片

图14所示的多目标优化结果验证了所提方法在改善船舶静水阻力和碎冰阻力方面所具有的显著效果，表6给出了原始船型与优化船型的静水阻力和碎冰阻力数值仿真结果对比。由表可看出，经所提方法优化后的船型能够能有效改善静水阻力和碎冰阻力，其中碎冰阻力的优化效果尤为明显，降阻效果可达10%以上。碎冰阻力的优化效果明显而静水阻力的优化效果较弱的原因在于，静水阻力与船首形状、船体流线型设计以及船尾形状等因素均相关，故仅考虑船首形状对静水阻力的影响效果有限；而碎冰阻力主要受船首形状的影响，与船中和船尾形状无关^[3]，故其优化效果非常显著。

表  6  原始船型和优化船型阻力计算结果对比

Table  6.  Comparison of resistance calculation results between original and optimized hull forms

船型	静水阻力 Rc /N	静水阻力 变化幅值/%	碎冰阻力 Ri /N	碎冰阻力 变化幅值/%
原始	6.87		20.32
优化	6.71	−2.32	18.17	−10.58

下载: 导出CSV 

| 显示表格

图16展示了原始船型与优化船型的时历碎冰阻力计算结果。从中可以看出，尽管两者的碎冰阻力随时间变化的总体趋势相近，但优化船型的阻力峰值更低，峰值出现的次数也更少；在原始船型出现碎冰阻力峰值的时刻，优化船型并未出现相应的峰值。并且在大多数情况下，优化船型的碎冰阻力要略低于原始船型。这些因素使得优化船型的碎冰阻力得到了改善。尽管优化后船型的碎冰阻力性能得到了改善，但是由于受流体影响的碎冰运动具有随机性，因此依然不可避免地出现了少量的碎冰阻力峰值，不过其出现峰值阻力的数目频率很低，不会对整体阻力产生很大影响。

图  16  原始和优化船型碎冰阻力仿真结果对比

Figure  16.  Comparison of ice resistance simulation results between original and optimized hull forms

下载: 全尺寸图片 幻灯片

图17展示了原始船型与优化船型的兴波对比。在数值仿真过程中，采用STAR-CCM+的标量幅值记录功能记录最大和最小兴波高度。结果显示，优化船型在航行过程中所产生兴波的高度幅值显著低于原始船型，说明优化船型在航行中因波浪引起的能量损失更少，从而解释了优化船型静水阻力更低的原因。兴波高度幅值降低的原因在于，优化后的船首球鼻艏更加尖锐、细长，这种设计可以更加有效地分割和引导水流，减少流体分析并抑制波浪的生成，从而降低兴波高度。

图  17  船舶兴波对比图

Figure  17.  Comparison of ship wave making

下载: 全尺寸图片 幻灯片

由原始船型与优化船型的峰值阻力对比可以看出，优化后船型的峰值阻力数目明显更低，这可以从CFD & DEM的标量场来细致分析其原因。标量场细致地展示了每一时刻碎冰受流体的影响，以及碎冰和船舶的相互作用。选取3个峰值时刻进行分析。在前2个时刻，原始船型出现了显著的阻力峰值，而优化船型却未出现类似的峰值，而到第3个时刻时情况相反。细节放大图显示，这是因为优化船型的艏倾角更小，在面对浮冰时能够更有效地引导浮冰下压并排开；相比之下，原始船型则是直接与大块浮冰发生正面碰撞（图18），从而导致显著的阻力峰值。这种减阻策略与船舶在与层冰碰撞时所采用的策略类似。在第3个时刻，优化船型的阻力显著增加，这是因为部分碎冰大幅翻转并贴附于船体表面（图19），从而导致了异常峰值阻力的出现。然而，这种阻力峰值并非因浮冰直接撞击船体而产生，所以对船舶的安全性能影响较小。此外，这种情况发生的频率有限，对整体阻力的均值影响不大。综合以上3种情况，说明艏倾角减小的优化措施可以降低船舶所承受的峰值阻力，能避免船舶在一瞬间因承受过大的阻力而导致结构发生损坏，保证了船舶结构安全。

图  18  船舶在碎冰中航行比较

Figure  18.  Comparison of ship navigation in floating ice

下载: 全尺寸图片 幻灯片

图  19  碎冰贴附导致的异常阻力峰值

Figure  19.  Abnormal resistance peak caused by ice adhesion

下载: 全尺寸图片 幻灯片

此外，由图18中的细节方法图和俯视对比图还可以看出，优化船型周围的浮冰运动更加平稳，并未出现剧烈的翻转和扰动现象，这使得船舶能够更加平稳地通过冰区，减小了浮冰与船体表面之间的接触和摩擦，从而降低了相应的阻力。这一现象也从侧面表明优化船型周围的流体环境更加稳定，湍流、旋涡以及局部高压区域的产生更少，有助于降低船体的静水阻力。

适当扩大水线形状有助于降低碎冰阻力，因为水线角度会影响船体与冰块的横向相互作用。与原始尖锐的水线相比，优化后的扩张水线能形成波浪轮廓，如图20所示，这有助于有效地将浮冰推离船体，从而减少冰块对船体的直接冲击和摩擦，这一结果与Huang等^[2]的研究结论一致。波浪轮廓的形成可以减少船体与冰块的接触，从而改善船舶在碎冰区域的航行性能。

图  20  船舶水线附近横向速度场

Figure  20.  Transverse velocity field near the waterline

下载: 全尺寸图片 幻灯片

4.   结　论

本文提出了一种基于CFD & DEM方法的多目标船型优化方案，综合优化了碎冰阻力和静水阻力，并在优化流程中考虑了流场及随机碎冰场对碎冰阻力的影响，从而获得了更精确的优化结果。为了模拟真实环境中的碎冰阻力，建立了一种HMIGA算法用来生成符合实际分布的碎冰场，并采用集成学习中的XGBoost模型来降低CFD & DEM方法成本过高的问题，提升了优化方法的应用效果。最后，以极区航线的运输主体集装箱船为例进行验证，主要得到如下结论：

1） 所提方法能够有效解决碎冰环境下的阻力优化问题。在优化过程中，考虑了流场和碎冰场分布的随机性对冰阻力的影响，克服了基于经验公式的优化方法在碎冰环境中的局限性。基于KCS船的优化实例验证了所提方法能够在保持船舶快速性的同时显著降低碎冰阻力。优化后的船舶在碎冰环境下阻力降低10.58%，在静水环境下阻力降低2.32%。

2） 优化后船型球鼻艏的上部位置下移且更加尖锐，艏倾角减小，船舶水线形状得到适当的扩张。

3） 尖锐且细长的船首设计降低了兴波高度，减少了能量损失，进而减少了静水阻力。艏倾角的降低使得优化船型在峰值阻力频率上明显较低，从而整体降低了阻力水平，提高了结构安全性。此外，适当的型线改变使得优化船型能够产生波浪轮廓而将浮冰推离船体，从而进一步降低了碎冰阻力。

4） 与其他代理模型相比，XGBoost模型在小数据量的情况下性能表现更好，降低了所需仿真数据量，能提供整体优化效率。

研究结果显示所提方法能够显著优化船舶的碎冰阻力和静水阻力，可为在碎冰航区航行的船舶设计提供重要的技术支持。这种方法能够在概念设计阶段快速生成可行的设计方案，有助于设计出更适合在碎冰航区航行的船舶。不过，本文的研究仅假设碎冰为不可破碎单元，并且也仅在无波浪条件下对船舶的碎冰阻力进行了优化，而北极航区的波浪虽然浪级较小，但对碎冰运动还是存在一定的影响，进而影响船舶的阻力性能。因此，未来将进一步探讨浮冰破碎以及波浪对极区船舶设计的影响，以完善本方法的适用性和实际运用价值。