基于NSGA Ⅱ的七叶侧斜螺旋桨多目标优化设计

孙鸣远，刘昕瑶，张 皓，严 俊

(1.南方海洋科学与工程广东省实验室(湛江), 广东 湛江 524000； 2.中国核动力研究设计院，成都 610213)

1 引言

螺旋桨处于船舶后的非均匀流场中，螺旋桨旋转过程中，不同半径处桨叶剖面的来流攻角会由于流场的不均匀性而随时间变化，从而导致螺旋桨受到的水动力呈现一定的脉动特征，这是螺旋桨产生激振力的主要原因。螺旋桨噪声与激振力密切相关，而在大量的试验数据和理论研究表明，通过改善螺旋桨的侧斜分布可以很大程度上降低螺旋桨产生的激振力。这是由于具有良好侧斜分布的螺旋桨可以使不同半径处的桨叶剖面不同时进入船后的高伴流区域，从而有效的降低螺旋桨在非均匀流场中产生的激振力，同时也可以降低螺旋桨由于振动等产生的机械噪声[1-4]，从而减小船舶航行过程中产生噪声对海洋牧场、深海养殖、海洋生物生态等渔业产业的影响。

而螺旋桨的桨叶螺距分布直接影响着各半径桨叶剖面的迎流攻角，对螺旋桨的水动力性能有着重要影响。螺旋桨设计中要求振动和噪声尽可能低同时螺旋桨效率较高，然而在螺旋桨实际设计过程中，降低激振力和提高螺旋桨效率往往是矛盾的，这需要在效率指标和振动指标上进行权衡优化设计。同时随着船舶节能减排方向的研究发展，国际拖曳水池会议(ITTC)[5-6]对于螺旋桨的优化方向也逐渐受到关注。

螺旋桨优化中对螺旋桨桨叶参数拟合采取的方法包括样条曲线法[7-8]、多项式拟合[9]等方法，基本可以较为精确的对桨叶参数进行拟合。在螺旋桨数值计算方法上大部分采用面元法[10-12]或CFD方法[13-14]，仿真理论已经较为成熟，相对面元法而言，CFD方法由于可以较为真实的模拟螺旋桨空泡形成、脱落、溃灭等现象，已越来越多的应用在螺旋桨优化方向。在优化理论的选择上较为常用的包括模拟退火算法[15]、粒子群算法[16-17]、遗传算法[18]等，其中NSGA II算法由于其良好的全局最优性、种群多样性、计算复杂度低等优点成为多目标优化中的经典算法之一。

为探究七叶侧斜螺旋桨非定常空泡流中降低激振力和提高螺旋桨效率的优化平衡问题，采用CFD数值方法结合NSGA Ⅱ多目标优化算法对某七叶侧斜螺旋桨的敞水效率和激振力进行了优化，螺旋桨侧斜及螺距沿桨叶半径分布通过B样条曲线进行拟合，分析结果表明所采用优化策略对提高七叶侧斜螺旋桨敞水效率和降低激振力是有效的。

2 螺旋桨优化方法

2.1 基于精英策略的非支配排序遗传算法

遗传算法在求解非确定性多项式问题上相比传统的搜索算法有明显的优势。非支配排序遗传算法[19](non-dominated sorting genetic algorithm，NSGA)通过样本中个体之间的支配关系分层实现，求解得到的非劣最优解的分布较为均匀，但是计算过程比较复杂且无精英策略。改进的基于精英策略的非支配排序遗传算法NSGA Ⅱ[20](non-dominated sorting genetic algorithm Ⅱ，NSGA Ⅱ)对NSGA进行了改进，基于精英策略的算法大幅改进了程序计算复杂、共享参数依赖性较大等问题，大幅提高了计算效率。同时NSGA Ⅱ算法在选择交叉的个体时，省去了共享函数、共享半径等的求解过程，实现快速非支配排序。同时NSGA Ⅱ算法引入了拥挤距离的概念用来保证最优解集的多样性和均匀性，结合其精英策略保证种群进化过程中不会丢失最优解，NSGA Ⅱ算法相比于NSGA算法有很大的提升。

2.2 螺旋桨非定常水动力性能预报方法

螺旋桨非定常水动力性能预报采用CFD数值模拟方法，对三维不可压缩流动，控制方程如下：

(1)

对流场中流体微元的动量守恒方程可以写为：

(2)

式中:ρ为流体密度；u、v、w为流体微元的速度分量；p为流体微元压力；Fx、Fy、Fz分别为流体微元三个方向的质量力。

采用SSTk-ω两方程湍流模型，该模型结合了标准k-ω方程和标准k-ε方程的优点，在模拟螺旋桨旋转等复杂流动时可以在近场和远场模拟中都取得较好的效果，相关研究在文献[21-22]中有详细论述，在此不详细说明。

2.3 螺旋桨几何参数化表达及优化理论

对螺旋桨进行多目标优化首先需要对螺旋桨的几何模型进行参数化表达，常用的参数化表达方法包括B样条曲线拟合[7-8,23-24]、贝塞尔函数曲线拟合[25]、多项式拟合[26]等。B样条曲线拟合方法由于有着局部可控性、曲线光顺性等优势从而有着广泛的应用。本文采用B样条函数对螺旋桨侧斜及螺距沿桨叶半径的分布进行曲线拟合，并将B样条曲线的控制点作为设计变量。图1所示为七叶侧斜螺旋桨曲线拟合分布与原始分布。图中螺距分布由N1～N4四个点进行控制，侧斜分布由N5～N9五个点进行控制。从图中可以看出，所选取B样条曲线控制点个数拟合螺旋桨几何参数分布的精度是可以接受的。

图.1 螺旋桨参数拟合与原参数分布拟合曲线Fig.1 Comparison of propeller parameter fitting and original parameter distribution

NSGA优化算法最初由Srinivas提出，其理论核心包括基于非支配排序原理的种群分级和基于共享小生境原理的虚拟拥挤度计算2个方面，在多目标优化领域有着广泛应用，但随着复杂度的增加，NAGA算法存在计算复杂度过高、需要人为制定共享参数等缺陷。因此，在前人的研究基础上，Deb等学者提出了改进型的NAGA-Ⅱ型算法，提出了精英保留、非支配排序等概念，同时提出了拥挤度的概念，保证了种群的多样性，大大提高了算法的适用度，同时也扩展了Pareto解集的分布范围。

拥挤度是指种群中给定个体的周围个体密度，直观上可表示为个体i周围仅仅包含个体i本身时的最大长方形的长，如图2所示，优化目标函数为f1和f2，个体i的拥挤度为d。拥挤度的引入实际上是为了保证在交叉运算中不至于出现基因过度集中的现象。

图2 拥挤度概念示意图Fig.2 Schematic diagram of congestion

两目标优化问题的拥挤距离的计算方法如下：

(3)

式中：L[i]d为第i个个体的第m个目标函数值；L为个体之间的距离。

2.4 螺旋桨侧斜及螺距优化方法

所进行的螺旋桨侧斜及螺距优化，是基于2.2节中螺旋桨非定常水动力数值计算方法及2.1节中多目标优化算法，将2.3节中的拟合参数(N1～N9)在合理范围内取值，从而得到一系列满足设计目标及约束条件的螺旋桨设计方案。该优化方法的具体做法如下：

优化变量：

N1～N9

优化目标：

min(KFX(N1,N2,…,N9))

(4)

min(-η0(N1,N2,…,N9))

(5)

约束条件：

(6)

因此螺旋桨优化的流程为：通过CFD数值方法预报得到螺旋桨敞水推进效率以及轴向一倍频非定常力系数等优化目标，由NSGA Ⅱ算法进行优化计算，根据优化目标及约束条件满足情况，重新设计优化变量，得到新的螺旋桨设计方案再次开始CFD数值计算，得到优化目标后再次进行优化结果评判、重新设计优化变量，直至满足优化迭代要求为止。

3 数值方法及计算模型

3.1 数值方法验证

螺旋桨优化选取原型桨为某七叶侧斜螺旋桨E1619，螺旋桨几何参数如表1所示，几何模型如图3所示。

表1 E1619桨几何参数

图3 E1619桨几何模型示意图Fig.3 E1619 propeller geometric model

螺旋桨非定常水动力数值计算域为圆柱体，计算中将计算域划分为静止域和旋转域2部分，旋转域同样为圆柱体，计算域尺寸及边界条件划分如图4所示。静止域与旋转域之间边界条件为Interface，旋转域与静止域之间的数值传递通过Interface进行，因此为保证数据传递的准确性，需保证Interface的网格大小基本一致，同时在Interface两侧各设置一层网格大小相同的边界层，网格划分结果如图5所示。

图4 E1619桨计算域划分示意图Fig.4 E1619 propeller computational domain division

图5 网格划分结果示意图Fig.5 Mesh division results

图6 数值计算结果与试验结果曲线Fig.6 Comparison of numerical calculation results and experimental results

3.2 网格收敛性分析

网格划分情况对计算结果的影响较大，因此数值计算需要对网格进行收敛性分析，只有当网格数量的增加对计算结果影响很小时，此时的计算结果才具有可信度。参照文献[26-27]的方法进行网格收敛性分析，不同网格方案满足一定细化率要求：

k=(Nf/Nc)1/d

(7)

其中，k为细化率，取为1.2；Nf为加密后网格数量；Nc为加密前网格数量；d为网格空间维度，取为3。通过上述网格细化策略得到3组不同的网格方案Mesh a、Mesh b和Mesh c，三组网格数量分别为6.53×106,1.16×107,2.10×107。网格划分基本尺寸均为螺旋桨直径D，桨叶表面网格边界层均设置为6层，网格划分参数如表2所示。

表2 网格划分参数

按照上述网格划分方案，计算螺旋桨设计工况(J=0.7)时的水动力特性，表3给出了3种网格方案计算得到的螺旋桨水动力参数。可以看出，随着网格的加密，数值计算结果与试验值的差异逐渐减小。采用Mesh b网格时，推力系数与试验值误差为2.16%，敞水效率与试验值误差为0.7%，已经满足计算精度要求，因此后续计算中均采用Mesh b网格进行，桨叶表面网格划分如图7所示。

表3 不同网格划分下KT and η0与试验值参数

图7 螺旋桨表面网格划分示意图Fig.7 Mesh division of the propeller surface

4 计算结果及讨论

对N1～N9共计9个优化变量在螺旋桨设计工况(J=0.7)处对螺旋桨进行多目标优化，优化目标如式(4)、式(5)，约束条件ε=5%，即螺旋桨推力系数KT范围为KT≥0.232 8。9个优化变量在±10%范围内进行变化，设定种群大小为15，迭代次数为20次，共生成300个个体，其中满足设定约束条件的可行解有232个。图8为可行解的气泡状分布图，图中气泡直径表示螺旋桨扭矩系数，气泡颜色表示螺旋桨推力系数。图中红色虚线圈中颜色为红色及粉红色圆圈表示不满足推力系数约束条件的解，剩余部分为Pareto解，可以看出Pareto解集在Pareto前沿上分布较为均匀。

图8 可行解气泡状分布图Fig.8 Bubble diagram distribution of feasible solutions

参照图8计算结果在Pareto前沿上选取四型螺旋桨进行分析：P236、P177、P102和P42。其中P236为满足推力系数约束条件下一阶轴向脉动力指标最优的方案；P177为满足推力系数约束条件下敞水效率性能最优的方案；P102为满足推力系数约束条件下综合考虑一阶轴向脉动力指标及敞水效率性能的折中方案；P42为不考虑推力系数约束条件下的最优方案。表4为上述几型螺旋桨设计参数及设计目标参数。

表4 优化桨与原桨几何参数

从表4中可以看出：优化螺旋桨敞水效率均较原桨有一定程度提高，一阶轴向力脉动值均较原桨下降，同时都相对原桨有一定程度的推力损失。在考虑推力系数约束条件限制下，P236桨敞水效率提高0.8%，一阶轴向力脉动水平下降10%；P177桨敞水效率提高1.7%，一阶轴向力脉动水平下降2.8%；折中方案P102桨敞水效率提高1.4%，一阶轴向力脉动水平下降6.7%。P42桨敞水效率提高2.4%，一阶轴向力脉动水平下降9.6%，是所有优化方案中优化效果最好的，可惜不满足推力系数限制条件，因此该方案不可取。

图9和图10分别为3组优化桨和原桨的螺距比分布和侧斜分布曲线。

图9 优化桨与原桨螺距分布曲线Fig.9 Comparison of the pitch distribution between the optimized propeller and the original propeller

从图9中可以看出，除个别半径处优化桨螺距比原桨略大外，优化桨整体螺距比分布相较原桨总体呈减小的趋势。其中P236桨和P177桨0.7R处螺距相比原桨降低了9.9%和11.3%，折中方案P102桨较原桨降低了4.8%。

图10 优化桨与原桨侧斜分布曲线Fig.10 Comparison of the skew distribution between the optimized propeller and the original propeller

从图10可以看出，在r大于0.7R后满足推力系数约束条件的3组优化桨的侧斜均较原桨有所提高。P236桨和P177桨0.9R处的桨叶侧斜角度较原桨增加4.2%和4.8%，折中方案P102桨较原桨增加了2.2%。整体来看优化桨与原桨的侧斜分布变化幅度较小。上述数据在一定程度上说明，将原型桨降低各半径处螺距比、提高桨叶靠近叶梢位置处的侧斜角度是提高推进效率、降低一阶轴向脉动力水平有益的探索方向。

图11为原桨与3组优化桨推力系数时历曲线，推力系数在一个稳定值附近随时间脉动，推力系数时历曲线符合螺旋桨推力脉动特性。从图中已标出推力系数脉动幅值来看，3组优化桨中P236桨脉动幅值最小，P177桨脉动幅值最大，P102桨则介于上述二者之间，该结果与3组优化桨一阶轴向力脉动计算结果相符。

图11 原桨和优化桨的推力系数时历曲线Fig.11 Comparison of thrust coefficient time-history curves of the original propeller and the optimized propeller

5 结论

1) 采用基于CFD数值方法与多目标优化方法结合的方式，对某七叶侧斜螺旋桨进行螺距及侧斜分布等多目标全局优化， Pareto解集在Pareto前沿上分布均匀，表明所采用的NAGS-Ⅱ多目标优化方法可有效提高螺旋桨推进效率，降低一阶轴向力，所采取的优化策略可行。

2) 优化后的螺旋桨总体螺距相比原型桨较小，侧斜分布在r大于0.7R后有所增加，在考虑螺旋桨激振力更优情况下选择P236桨，一阶轴向力脉动值下降10%，考虑效率更优时选择P177桨，敞水效率可提高1.7%。

3) 七叶侧斜螺旋桨优化中，将原型桨降低各半径处螺距比(r<0.7R)、提高桨叶靠近叶梢位置处(r>0.7R)的侧斜角度是提高推进效率、降低一阶轴向脉动力的研究方向。