渔船球鼻艏参数优化研究

李广年，于 群，杜 林，秦海瑞，袁文鑫

（1.宁波大学海运学院，浙江宁波 315211；2.上海市东方海事工程技术有限公司，上海 200011）

0 引 言

近年来，由于过度捕捞及环境污染导致近海渔业资源匮乏。为缓解我国渔业压力、建设新型现代化渔业体系、推进海洋强国战略部署，发展远洋渔业势在必行。2020年“中国远洋渔业履约白皮书”提出支持渔船更新改造，鼓励采用安全高效、节能减排的技术和船用装备。船舶参数化水动力构型为绿色节能船型的研发提供了极大的便利，是船舶行业实现“碳达峰”、“碳中和”目标的技术途径之一。球鼻艏可在一定程度上减小船舶阻力，但其泛化性较弱，即同一球鼻艏应用到其他船时并不一定能够获得满意的减阻效果。因此，球鼻艏形状、参数等与船体有恰当的配合是球鼻艏方案设计、优化的重点方向。

国内外研究人员针对球鼻艏做过一系列的研究。Chrismianto 等[1]使用4 个设计参数对球鼻艏进行参数化设计，利用多目标遗传算法对KCS（KRISO Container Ship）进行优化。Lu等[2]通过参数化方法对球鼻艏进行变形，并利用改进非支配排序遗传算法（Non-dominated Sorting Genetic Algorithms，NSGA-II）对一艘大型集装箱船进行总阻力优化，结果表明优化后的球鼻艏使船体总阻力降低2.845%。Huang 等[3]利用径向基函数（Radial Basis Function，RBF）和多目标人工蜂群优化算法对一艘货船球鼻艏进行优化，并通过计算和试验验证得到了最优船型。Hong等[4]结合CFD（Computational Fluid Dynamics）技术和DoE（Design of Experiment）技术，以最小阻力为目标对某渔船球鼻艏进行了优化设计。Yu等[5]使用粒子群算法对一艘散货船船艏进行优化研究，结果表明最优船型的兴波阻力降低13.2%。随着多目标与多学科优化的兴起以及设计变量维度的增加，优化设计往往需要海量的计算，基于近似模型的优化技术是解决这种复杂优化的方法之一。Lin 等[6]以Kriging 模型作为替代模型，使用NSGA-II对某双艉渔船阻力性能进行优化。Guerrero 等[7]以两个球鼻艏参数为设计变量，使用Kriging 模型作为替代模型并进行优化，最终使阻力降低约7%。Feng 等[8]采用支持向量回归（Support Vector Regression,SVR）替代模型和NSGA-II方法对某近海养殖船进行优化，优化船型的阻力比原始船型减小了1.59%。Coppedè 等[9]基于几何变形方法、CFD 计算、普通Kriging 模型以及高斯回归的响应面法搭建的优化框架对KCS进行优化，使总阻力减小4%。综上所述，尽管专家学者们针对球鼻艏优化做了大量的工作，但探究多个球鼻艏特征参数对渔船阻力的影响并将参数化方法运用到渔船球鼻艏优化的并不多见。

本文以一艘远洋渔船为研究对象，以设计航速（Fr=0.28）时总阻力最小为目标，介绍了一种渔船球鼻艏参数化优化方法。首先，选取渔船球鼻艏特征参数和特征曲线对其进行参数化建模，以球鼻艏的几何参数作为设计变量，利用敏感性分析方法结合CFD技术对各变量进行灵敏度分析，甄别出对阻力性能变化最为敏感的若干变量；然后，针对敏感变量使用穷举搜索法对样本空间进行采样，并通过数值仿真计算样本总阻力；最后，利用多项式构建替代模型得到设计变量与总阻力系数之间的函数关系并使用遗传算法进行优化，以获得阻力最优的球鼻艏参数。

1 基本理论与方法

本文将数值预报、灵敏度分析、最优化技术等方法应用于渔船球鼻艏优化问题。其中，数值预报是建立渔船球鼻艏优化问题数学模型的基础，数值计算精度直接影响了优化结果的质量；灵敏度分析可以去除球鼻艏模型中对阻力变化影响较小的参数，降低整个优化的计算成本及时间成本；最优化技术是解决渔船球鼻艏优化问题的必要手段，通过优化算法能够在设计空间内快速、准确地搜寻最优解，即最佳渔船球鼻艏。

采用Numeca公司的FINE/Marine 模块进行数值模拟，在等温条件下，针对所求解的不可压缩黏性流，质量和动量守恒方程分别为

式中，t为时间，ρ为密度，V为控制体，S为围成控制体的面积，Ud为S上n方向的速度，U和p分别为速度和压力，Ui为在xi坐标轴方向上的平均速度分量，τij和gi分别为黏性应力张量和重力矢量，Ii和Ij为方向向量，ci为i流体的体积份额。

使用k-ω(SST-Menter)湍流模型，其湍流能k和湍流耗散频率ω满足：

式中：xj为坐标轴；Uj为在xj坐标轴方向上的平均速度分量；μ为分子黏度；μt为湍流涡黏度；τtij为湍流雷诺应力张量；Sij为平均应变率张量；F1为辅助混合函数；Pω为ω的导出项；β*、σk和σω2分别为湍流模型常数。

在设计初期，可以通过灵敏度分析快速有效地确定灵敏度较大的参数，为后续球鼻艏优化工作节约时间及计算成本。一种典型的灵敏度分析方法是OAT(one-at-a-time)，它可以直观地观察输入因素对输出的影响。本文以球鼻艏的特征参数为设计变量，使用OAT 进行灵敏度分析，具体内容为：改变一个设计变量，其它设计变量保持不变，记录设计变量的变化对总阻力系数的影响；然后将该设计变量返回初始值，对其它设计变量进行重复操作；最后通过分析每个设计变量对总阻力系数Ct的影响，将灵敏度最高的设计变量作为优化的设计变量。

选择有限个点分析设计变量与目标方程的响应，即确定替代模型。可近似替代设计变量与目标方程的关系，从而减少计算量。建立替代模型的前提是对样本空间的选择，穷举搜索法（Exhaustive Search，ES）是一种具有规则结构，能够系统地列举多状态网格的所有可能情况，并检查每个情况是否符合问题描述的方法。本文利用ES方法将设计变量的约束区间细分为给定数量的等距间隔，逐一匹配确定样本点，并通过CFD计算对应的Ct。

采用二元多项式构建替代模型，其多项式拟合模型定义如下：

式中：f(x,y)即为总阻力系数Ct；x、y为自变量；a1、a2、…、ai、…、aj为拟合系数，a0为截距项；m、n为拟合次数。替代模型生成之后，需要对其进行精度验证，确保其可以精确地反应出设计变量对Ct的响应。拟合优度（goodness of fit）是指回归曲线对真实值的拟合程度，度量拟合优度的统计量是确定系数Rsquare。

式中，yi为样本的真实值，y^i为替代模型预测值，yˉi为样本均值，SSR为预测值与真实值均值差的平方和，SSE为真实值与预测值差的平方和，SST为真实值与真实值均值差的平方和。Rsquare给出了目标变量的可变性程度，它介于0 到1 之间。Rsquare越高，说明模型拟合得越好。在拟合模型过程中，当模型Rsquare≥0.95时，目标变量对总阻力系数的响应满足精度要求。

遗传算法（Genetic Algorithm，GA）是借鉴生物进化过程而提出的一种启发式搜索算法。其主要特点是直接对结构对象进行操作，且不存在求导和对函数连续性的限定，采用概率化的寻优方法，不需要确定的规则就能自动获取和指导优化的搜索空间，自适应地调整搜索方向。其通用表达式为

式中：f(xi)为目标函数；x1,x2,x3…xn为自变量；s.t.为约束条件。本文在使用GA进行优化时，首先，以总阻力最小为目标，随机生成目标的初代解作为初代种群；其次，对种群中的个体，即设计变量进行编码，以替代模型作为个体的适应度，计算种群中每个个体的适应度；接着，根据适应度的高低选择参与繁衍的父体与母体，并对被选出的父体与母体执行遗传操作，即复制父体与母体的基因，并采用交叉、变异等算子产生出子代；最后，找出所有子代中适应度最高个体，即总阻力最小的设计变量作为最优解。

2 问题描述

在进行优化前，首先要进行CFD 精度计算和参数化建模。选择KCS 船型进行CFD 精度验证，数值模拟计算域尺寸及边界条件如图1所示。采用全六面体非结构网格，网格的生成方式是由体到面，并将物体表面附近网格适当细化并投影到物面上形成贴体网格。为了准确模拟流动和船体阻力，在自由液面处和壁面处进行网格加密，整个计算域的网格数量约为134万，船体网格划分如图2所示。

图1 边界条件Fig.1 Simulation domain and boundary conditions

图2 网格划分Fig.2 Mesh division

通过对比KCS在不同傅汝德数下总阻力系数计算值与MOERI拖曳水池试验数据，如图3所示，可以看出二者相对吻合，平均误差在1.13%左右。说明本文所采用的数值分析方法适用于船舶阻力计算，其计算精度可达到预期效果。

图3 KCS阻力数据对比Fig.3 Resistance curve of KCS

本文研究对象为一艘金枪鱼船，如图4 所示，该船球鼻艏最大横剖面位于首柱附近，球鼻艏的长度Lb为2.43 m，球鼻艏的最大横剖面处的最大宽度Bb为1.07 m，球鼻艏纵剖面最前端点的高度Zb为1.88 m，球鼻艏中纵剖面面积Abl为5.63 m2，球鼻艏中最大横剖面面积Abt为5.14 m2。

图4 球鼻艏示意图Fig.4 Profile of bulbous bow

金枪鱼船的缩尺比λ=16，主要船型参数如表1所示。在球鼻艏建模过程中，针对球鼻艏主要参数及几何特征，共有7个球鼻艏特征参数，其初始值及取值范围如表2所示。表2中，X1和X2是与Abl相关的特征参数，X4与Abt相关。同时还设置了特征曲线，它们分别为最大宽度曲线以及最大宽度高度曲线，即Bb在船长方向上，从最大横剖面处至球鼻艏最远点宽度和高度的变化曲线。

表1 金枪鱼船的主要船型参数Tab.1 Principal hull form parameters of the tuna seiner

表2 球鼻艏特征参数Tab.2 The original seven parameters of bulbous bow

本文使用OAT方法，在不考虑多个变量耦合的情况下对每个变量进行分析，在每个参数所在的取值区间内均分取样，每个参数都得到5 个样本，其样本取值随着编号的增大而增大。如图5 所示，在变化区间内随着每个设计变量的改变，其对应Ct的变化趋势和变化范围都有所不同。每个设计变量对应的Ct曲线起点不同，但最终近似汇集于两处。第一处相关变量为X6和X7，X7对应的Ct不减反增，X6对应的Ct变化是波动的；X6和X7对应的Ct总体变化量分别为0.957%和-1.451%。第二处相关变量为X1～X5，各自对应的Ct总体变化量分别为-3.650%、-2.317%、-4.035%、-1.363%和-1.663%；X3对应的Ct在变化区间内单调下降，X4、X5对应的Ct变化趋势大致相同，都是先增加再减小，X1和X2对应的Ct变化也呈下降趋势，但在第三个样本处略有增加。

图5 球鼻艏特征参数灵敏度分析Fig.5 Sensitivity analysis of seven parameters

综上所述，在不考虑变量之间的耦合情况下，对每个设计变量进行灵敏度分析。就整体而言，X1和X3对应的Ct总体变化量比其它五个参数要大得多。进一步考虑X1和X3的耦合效果对总阻力性能的影响，其他变量均取初始值。以X1和X3为优化问题的设计变量，使用ES 取样，共生成100 个样本，经过ES抽样的空间点和平面投影点分布如图6所示。

图6 ES取样分布Fig.6 Uniform sampling on the design space

可以看出，在X1和X3耦合的情况下，其样本点的变化趋势与单一变量灵敏度分析大致相同，Ct随X1和X3增大而减小，尤其是在0.65～0.80 范围内，Ct下降的趋势较大。综合考虑模型精度与计算量的情况下，最终选择二元五次多项式作为本文的代理模型，即

其中拟合系数如表3所示。

表3 拟合系数Tab.3 Fitting coefficients

在保证精度满足要求情况下，剔除了样本集中个别数据奇异点，经过计算，最终替代模型的拟合优度Rsquare为0.95，拟合曲面如图7所示。

图7 拟合曲面Fig.7 Polynomial regression surface

从图7 可以看出，随着X1和X3增大，Ct呈减小的趋势较为明显，在f( 0.6,0.6 )附近出现极大值，最小值则在f( 0.9,0.9 )附近。

引入替代模型后，优化问题可表示为

式中，目标函数为在设计航速（Fr=0.28）时Ct最小，约束条件为|Δo-Δb|≤1%Δb（其中，Δb为母型船排水量，Δo为优化船排水量）及对应变量的取值范围。

设置种群数量为50，迭代次数为100，交叉率和突变率分别为0.8 和0.01。GA 的进化终止条件为达到最大进化代数或进化过程中连续10代种群之间的适应度平均值变化小于设定值1E-6，优化结果如图8所示。

图8 GA迭代过程Fig.8 Optimization convergent curve

从图8（a）可以看出，优化在40代以后开始收敛，且种群适应度平均值的变化小于1E-6，从图8（b）可以看出，第40代之前个体的平均距离呈振荡状态，第40代以后个体的平均距离逐渐收敛，个体间的平均距离为0。

3 结果与讨论

运用GA 对替代模型进行优化得到的最优解见表4，最优解所对应的球鼻艏形状如图9 所示。从球鼻艏的横剖面可以看出，在横向上其形状变化不大，球鼻艏的最大宽度不变；在纵向上球鼻艏上半部分的形状较母型更加饱满，但是球鼻艏的长度没有变化，最远端的高度却有所下降。

表4 优化结果Tab.4 Optimization results

图9 球鼻艏形状对比Fig.9 Comparison of profiles

根据表5 阻力成分（Fr= 0.28）对比可知，在球鼻艏优化过程中，与Re相关的粘性阻力没有下降，反而增加0.16%，但与Fr相关的压阻力下降了5.43%。

表5 阻力成分（Fr=0.28）Tab.5 Occupation changes of total resistance（Fr=0.28）

图10 为母型船与优化船在Fr为0.20～0.30 工况下的总阻力系数对比。图中显示，在Fr为0.22、0.24、0.30 时总阻力下降较多，分别减小了3.32%、3.27%和3.57%；而在Fr为0.20、0.26、0.28时，总阻力下降相对较小，分别下降1.35%、2.96%和2.68%。可以看出，本文优化后的球鼻艏在低航速时（Fr≤0.20）减阻效果相对较差，Fr=0.30 时优化效果最为明显，设计工况下（Fr=0.28）下降2.68%。

图10 不同Fr下阻力系数对比Fig.10 Ct changing with Fr

图11 给出了母型船和优化船在Fr为0.20～0.30 时的艏部波高对比图。可以看出，在Fr=0.20～0.24时优化后的球鼻艏处的兴波明显低于母型船，随着航速的增大，优化船船首处的波高与母型船的差距逐渐缩减；当Fr＞0.28 以后，优化船艏部的波形与母型船逐渐相近。优化船周围的波形较母型船得到改善，波高低于母型船，波峰波谷影响的区域也较母型船有所减小，见图12。

图11 艏部波高对比Fig.11 Comparison of bow wave heights

图12 波形图对比Fig.12 Comparison of wave elevation contours

图13给出了母型船与优化船在Fr为0.20～0.30时所受到压力的对比图。可以看出，随着Fr的增加，优化船和母型船在球鼻艏下部受力基本一致，而球鼻艏上部的差别较为明显，优化船所受压力相较于母型船有所降低。另外，球鼻艏首端受到的压力随着Fr增大而逐渐增大，优化船的球鼻艏首端压力相对大于母型船压力。

图13 压力云图Fig.13 Pressure force contours

4 结 论

本文提出了一种利用CFD 计算与替代模型相结合的优化设计方法并成功应用于渔船球鼻艏减阻。该方法以某金枪鱼船球鼻艏为研究对象，进行参数化建模，对设计变量进行灵敏度分析，寻找对总阻力系数最敏感的两个设计变量；接着使用穷举搜索法进行抽样生成样本空间，并利用CFD进行数值计算；然后通过多项式拟合生成以特征参数为自变量，总阻力系数为因变量的二元五次多项式替代模型；最后利用遗传算法求解替代模型的最优解，获得最优球鼻艏形状。优化结果考虑了多个工况下的阻力性能，得出以下结论：

（1）通过灵敏度分析确定了两个对总阻力系数最敏感的设计参数，优化后的球鼻艏形状在垂向上变化较大，纵向上最远端点高度较母型船降低，横向上的变化较小；

（2）阻力分析表明，在设计工况条件下，优化后的船舶受到的兴波阻力下降5.43%，粘性阻力增加了0.16%；优化后的球鼻艏在各个工况下的总阻力系数都较母型船低，总阻力平均降低2.86%，在Fr=0.22～0.30时，优化的效果较好，在设计航速（Fr=0.28）条件下，总阻力下降2.68%，Fr=0.30时优化效果最为明显，总阻力减小3.57%；

（3）优化后船体周围的波形较母型船得到改善，艏部波高低于母型船，波峰波谷影响的区域也较母型船有所减小。随着Fr的增加，优化船和母型船在球鼻艏下部受力基本一致，而球鼻艏上部的差别较为明显，优化船所受压力相较于母型船有所降低。球鼻艏首端受到的压力随着Fr增大而逐渐增大，优化船的球鼻艏首端压力相对大于母型船压力。

通过替代模型建立球鼻艏特征参数与总阻力系数之间的关系，考虑到计算量和模型精度的影响，只对两个对总阻力系数最敏感的设计变量进行优化。后期将会考虑更多变量之间的耦合效应，同时也会采用精度更高的替代模型，更加深入研究球鼻艏优化问题。