基于自融合方法的船体外形优化研究

洪智超，宗 智，赵玫佳，杜金刚

（1.江苏科技大学，江苏镇江 212100；2.大连理工大学，辽宁大连 116024；3.中国造船工程学会，北京 100861）

0 引 言

船体几何外形优化主要包含三个方面[1-2]：船体几何重构技术、船舶性能评估预报和优化算法。设计变量受到船体几何重构技术的直接影响，而且设计变量是船舶性能预报和优化的输入参数，因此船体几何重构技术是整个优化过程中的基础。良好的船体几何重构技术应该具备以下几个优点：设计变量少、几何外形变形幅度大、生成的曲面光顺以及生成的船体外形具有实用性[3]。近年来，CFD因其高精度得到了广泛应用，但其对计算资源的消耗也较大，因此参数数量少的船体几何重构技术受到了广泛关注。Tahara[4-5]使用融合方法对一艘船的几何外形进行了优化，通过对两艘船进行融合实现船体几何重构。该方法的优点是使用的参数数量少，但船体几何外形的变化限制其用于融合的母型船之间。特征参数方法[6]使用的参数数量也较少，与传统的基于CAD 的几何重构技术不同；特征参数法并不使用分布在船体表面的点，而是使用形状参数以及一组纵剖线来描述船体外形。该方法的缺点是形状参数的选择对设计员的经验依赖性较大。

变形幅度大的几何重构技术由船体外形局部优化拓展到全局优化时，设计变量会急剧增加，从而导致优化时间大幅增加，因此此类方法一般用于船体外形的局部优化。Peri等[7]使用CFD 方法对一艘油轮进行了优化，使用Perturbation Surface 方法对球鼻艏进行几何重构，该方法的特征是变形幅度大，但因设计变量数量的原因不适合用于全局优化；冯梅[8]提出了一种基于径向基函数插值（radial basis function，RBF）的曲面变形方法，并将改进的船体曲面变形模块嵌入团队开发的船舶水动力性能多学科综合优化平台（SHIPMDO-WUT V6.0）；Chrismianto 和Kim[9]使用三阶贝塞尔曲线、曲线平面交叉法对球鼻艏进行了优化；Perez等[10]使用B样条曲面方法对球鼻艏进行了优化，优化中使用CFD方法对目标函数进行数值预报。FFD（free form deformation）[11-12]方法也是适用于船体外形局部优化的方法，该方法的优点是几何曲面可以任意变形。倪其军[13]使用SBD 技术对“探索一号”科考船的艏部线型进行优化，其船体几何重构的实现使用FFD 方法，并利用粒子群优化算法对线型优化设计问题进行求解，其总阻力减小了7.7%；冯佰威[14]使用基于径向基插值技术的船体曲面变形方法对DTMB 5415 展开阻力优化，兴波阻力最多减小了49.66%；周广利[15]使用NURBS 技术对DTMB 5415 船体几何重构展开研究，实现了船体高精度几何重构，且具有曲面变形、曲面拼接和曲面分割等功能。

本文使用自融合方法对一艘高速船展开局部及全局优化。自融合方法由融合方法[4-5]发展而来，二者的区别在于自融合方法使用船体横剖面，而不是将整个船体作为融合的基本元素。融合方法的基本元素为整个船体，因此需要寻找两条以上相似船体作为融合的基本元素。而自融合方法的基本元素来源于原始船体的横剖面，省去了寻找相似船体的麻烦。使用横剖面作为融合基本元素会增加设计变量数量，但由此造成的数值计算量的增加有限，仍然在可接受的范围内。使用横剖面作为融合基本元素能大幅增加船体外形的变化范围，且能同时对船体的局部及全局进行变形操作。此外，自融合方法生成的船体表面较为光滑，可直接用于后续的数值计算，无需进行人工修正。

本文使用具有6个变量的自融合方法对高速船展开局部及全局优化。使用CFD 方法预报船体总阻力系数（优化的目标函数）。CFD 方法对计算资源的消耗较大，为了降低计算量，节约计算资源，对该问题展开灵敏度分析并建立响应面（RSM）。灵敏度分析的结果表明，全局变量（控制全船曲面变形的变量）对目标函数的影响大幅高于局部变量（控制船体局部变形的变量）的影响。使用RSM 进行优化可大幅降低数值计算的工况数量，由56下降到100 以内。本文使用MIGA（multi-island genetic algorithm）对该问题进行优化，与传统的GA 相比，该方法能避免最优解陷入局部陷阱，可靠性得到大幅提高。优化后，该船的总阻力系数下降了9.76%。

1 船体自融合方法

船体自融合方法是船体几何外形优化的重要组成部分，该方法选取母型船的横剖面为基本元素并使用融合因子给其赋予权重，生成新的横剖面，最后使用新的横剖面重构船体三维曲面，获得新的船型。文中使用6个自融合因子和12个由原船中选取的基础横剖面进行船体横剖面的融合及船体三维曲面的几何重构。

1.1 自融合方法基本流程

自融合方法能快速高效地修改船体几何外形，其关键在于融合因子和从原始船型表面提取的横剖面，自融合方法的计算式见式（1），基本流程见图1。首先，根据原船基本特性及优化目标选取第一类横剖面，第一类横剖面用于后续的三维船体几何重构，所选取的第一类横剖面应尽可能完整地包含原船的外形及性能信息；然后，针对每个第一类横剖面，从原船选取两个第二类横剖面，用于融合产生第一类横剖面，并作为第一类横剖面的变形空间的边界；最后调整融合因子，产生新的第一类横剖面，并进行船体重构获得新船的几何外形。自融合方法的三个重要因素有：（1）横剖面的表达方式；（2）横剖面的选取；（3）约束条件（本文中的约束条件为船体总排水量保持不变）。

图1 自融合方法工作流程Fig.1 Flow chart of self-blending method

式中，CS表示横剖面，ω表示融合因子，i、j表示横剖面的编号。

1.2 横剖面表达方法

横剖面的表达应同时具有涉及参数少和几何外形描述精确的特性。本文通过定义横剖面轮廓线上的节点，并使用样条曲线连接各节点实现船体横剖面的表达，其示意图见图2。横剖面轮廓线上的节点分布方式是轮廓线构建的关键，使用水线或纵剖线与轮廓线的交点作为节点是一种较为简便的节点选取方式，见图3（a）～（b），但使用这两种方法选取节点会导致不同横剖面上相同序号节点的横坐标或纵坐标均为同一数值，降低了后续融合计算中节点移动的自由度，即降低了船体变形的幅度。本文使用轮廓线与辐射线的交点作为横剖面的节点，见图3（c），其融合计算式为

图2 船体横剖面示意图Fig.2 Sketch of hull section

图3 横剖面轮廓线节点分布方式示意图Fig.3 Sketch of section line nodes

式中，yN表示融合后所得到的新横剖面的轮廓线节点的y坐标，zN表示融合后所得到的新横剖面轮廓线节点的z坐标，ω表示融合因子，y表示用于融合的基础横剖面的轮廓线节点的y坐标，z表示用于融合的基础横剖面的轮廓线节点的z坐标，下标i、j分别为新横剖面、基础横剖面的序号，下标k为横剖面轮廓线节点的编号。

1.3 横剖面的选取

基础横剖面的选取对船体外形优化具有重要意义，直接关系到船体几何变形的空间。从原船选取的基础横剖面应包含船体外形的关键信息，尤其是与船舶性能高度相关的几何特征。本文的研究对象为高速方尾船，其主尺度见表1，横剖面图见图4，船体三维轮廓图见图5。该船后半段的形状变化较小，横剖面轮廓由两条直线组成。因此，在0～10站之间的横剖面中选取第0站和第9站两个横剖面作为第一类横剖面。对于船体前半部分，尤其是球鼻艏部分，曲率变化较大，选取4 个横剖面作为第一类横剖面，全船共计6个第一类横剖面（19.692、19.4、19、16、9、0），第二类横剖面选取与第一类横剖面相邻的横剖面。

表1 原船主尺度（模型尺度）Tab.1 Principal dimensions of original vessel(model)

图4 船体横剖面曲线图Fig.4 Body plan of the vessel

图5 船体的三维曲面图Fig.5 3D surface of original ship hull

2 数值计算方法验证

本文使用CFD 方法预报船体阻力，求解器为基于有限体积法的STAR-CCM+。使用VOF 方法捕捉水和空气的交界面，湍流模型为k-ω模型，网格划分使用STAR-CCM+的切割体网格。此外，本文使用重叠网格技术模拟船体的纵摇和垂荡运动。数值计算的不确定度分析使用Stern 提出的收敛性分析方法[13]。阻力预报试验由大连理工大学完成。

本文对设计航速下（Fn=0.367）船舶阻力预报数值模拟的迭代收敛性和网格收敛性展开分析。迭代收敛性通过总阻力系数的时历曲线进行分析，网格收敛性分析以三套不同尺度网格的数值计算结果为基础。总阻力系数的时历曲线见图6。迭代不确定度为时历曲线最后几个周期的振幅，即为时历曲线收敛后周期性振荡的幅值，为0.28%。不同尺度网格的数量及阻力预报结果见表2，网格尺寸的变化率为 2，该值的选取参照Stern[16]的论文。网格收敛性验证的结果见表3。网格收敛因子RG为0.123，小于1，即网格收敛性符合单调收敛。误差E为1.45%S，网格不确定度UG为0.37%，迭代不确定度为0.28%，模拟不确定度为0.22%，均为较低水平，说明使用该方法进行阻力预报是可行的。

图6 原始船型的CT时历曲线Fig.6 Time history of CT of original vessel

表2 不同网格尺度下原始船型的总阻力系数（Fn=0.367）Tab.2 Results of CTunder different grid sizes（Fn=0.367）

表3 网格收敛性分析（S为网格1尺度下总阻力系数CT的数值计算结果）Tab.3 Grid convergence analysis(S refers to the simulation result of CT under Grid 1)

3 CFD计算结果

本文开展了一系列数值计算以构建关联设计变量和CT的响应面关系式。在讨论响应面之前首先介绍一些基本概念。数值模拟空间是一个由设计变量的最小及最大值定义的范围，变量的水平是指变量的不同取值，响应是指总阻力系数CT的值。使用响应面进行优化的流程见图7。首先，使用灵敏度分析研究各设计变量对响应的影响，并以此为依据定义响应面的数值模拟空间。其次，使用中心复合设计（CCD）方法确定数值模拟计算工况并展开计算，根据结果构建响应面，并评估其可靠性。最后，以响应面为基础展开优化，获得最优解。

图7 使用响应面的流程图Fig.7 Flow chart of the application of RS

3.1 灵敏度分析结果

本文使用L25（56）正交试验表及直观分析方法对设计变量展开灵敏度分析，共包含25 个计算工况、6 个因子（设计变量），每个因子包含5个水平（均布在设计空间）。6个因子及灵敏度分析的数值模拟空间见表4，其中的4个设计变量（ω3，ω4，ω5，ω6）及其对应的横剖面19、16、9、0 用于控制船体的全局形状变化，3 个设计变量（ω1，ω2，ω3）及其对应的横剖面19.672、19.4、19用于船体的局部外形（球鼻艏）控制。

表4 设计变量空间Tab.4 Simulation domain of the design variables

正交试验表中所有工况的数值计算结果见图8，由图中的结果可知，25个工况中CT的波动幅度为22.98%。对正交试验的结果进行直观分析，其结果见图9，图中的结果表明因子ω6对总阻力系数的影响最大，说明方尾的形状对总阻力系数的影响最为显著。因子ω1、ω2、ω3对总阻力系数的影响远小于因子ω4、ω5、ω6的影响，说明局部形状的变化对总阻力系数的影响小于全局形状变化的影响，这与实际情况相符。根据图9的结果，最优解在融合因子的水平为[1，4，2，3，5，5]的邻域内。

图8 正交试验数值计算结果Fig.8 Calculation results of the orthogonal test

图9 设计变量直观分析图Fig.9 Main effect of screening design

3.2 响应面结果

直观分析能给出各因子对目标函数的影响趋势，但无法考虑因子间的相互影响，因此本文以灵敏度分析为基础，在其得出的最优解的邻域内进行CCD数值试验设计，CCD 数值试验设计中各因子的取值范围见表4，每个因子取5 个水平，共计77 个工况，其数值计算结果见图10。根据CCD 的数值计算结果构建一个基于二次多项式的响应面，该响应面建立了设计变量与总阻力系数之间的映射关系，且该响应面包含了因子间的相互影响，响应面的计算式为

图10 CCD数值计算结果Fig.10 Results of CT by CCD

为检验响应面的有效性，检验其是否能真实反映设计变量与目标函数之间的关系，对响应面展开方差分析，其结果见表5（表中SS为平方和，DF为自由度，MS为均方）。方差分析主要用于评估响应面，检验其是否能真实反映设计变量与目标函数之间的关系。方差分析的中心思想是通过对比由设计变量变化引起的目标函数波动与由随机误差引起的目标函数波动来判断响应面的结果是否显著。文中响应面的F值为21.9，说明该响应面是显著的，即该响应面是有效的。

表5 响应面的方差分析Tab.5 Variance analysis of RS

4 船体总阻力系数优化

本文以控制船体曲面变形的融合因子为设计变量对船舶的总阻力系数CT进行优化，其数学描述见式（4）。总阻力系数由式（3）定义，该式为以CCD 的数值计算结果为基础构建的响应面，其设计变量为6个融合因子。响应面在水平2 到水平4 之间具有较高的精度，因此本文的优化空间即限定在水平2和水平4之间的设计空间内，各因子的具体范围见表4。优化的约束条件是保持船体排水体积不变。优化算法为多岛遗传算法，多岛遗传算法的各参数见表6。初始种群的总阻力系数均布在0.005 441到0.005 594之间。

表6 MIGA参数Tab.6 MIGA parameters

式中，gi为约束条件。

多岛遗传算法具有一定的偶然性，本文对该优化问题进行10 次优化计算，结果表明10 次优化的最优解极为接近，设计变量之间的误差小于0.05%，目标函数之间的误差小于0.000 1%。其中一个优化计算的结果见图11 和图12。结果表明，经过20 代的进化之后，总阻力系数最终收敛到0.005 427。最优解所对应的各设计变量值见表7。为验证该最优解，本文使用CFD 方法对该最优解进行数值计算，总阻力系数CT为0.005 422，与基于响应面的最优解吻合良好。优化后，总阻力系数下降9.76%。最优船型与原始船型的对比见图13，由图中的结果可知，优化后型宽、型深及球鼻艏变化较小，船体折边线垂向位置升高且船长略有增加。原船与优化后船型的波形对比见图14，由图可知，优化后船艏兴波小于原船。

图11 CT迭代过程Fig.11 Changes of CT in optimization

图12 CT迭代曲线Fig.12 Convergent results of CT

图13 最优船型和原始船型的对比图Fig.13 Comparison of geometies between original and optimized ship hull forms

图14 原船与最优船波形图对比（上半部为原船，下半部分为最优船）Fig.14 Wave patterns of original and optimized vessels(upper half:original hull;lower half:optimized hull)

表7 最优解Tab.7 Optimized solution

5 结 语

本文提出了一种船体变形方法——自融合方法，其核心思想是从原始船型中提取特征横剖面用于融合操作，生成新的横剖面并进行三维船体的几何重构，产生新的船型，具有使用变量少、生成的船体外形光顺等优点。使用6个融合因子和12个从原船提取的基准横剖面，结合CFD方法、响应面方法和多岛遗传算法，以总阻力系数最小为优化目标，排水量保持不变为约束条件，对一艘高速方尾船展开局部及全局优化，即对其全船外形和球鼻艏形状展开优化。优化后该船在保持排水量不变的前提下，设计航速状态下的总阻力系数下降了9.76%，证明了该方法的工程应用价值。