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分别为流体微元三个方向的质量力。
采用SST k-ω两方程湍流模型，该模型结合了标准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)
式中: KFX为螺旋桨一倍频轴向脉动力系数，该系数可表示螺旋桨轴向激振力,定义螺旋桨的轴频为 APF(axial passing frequency)Hz，一倍频轴向脉动力系数即对轴向力T进行傅里叶分析得到APF频率下的轴向力系数(KFX = 
其中
为傅里叶分析后APF频率下的轴向力)；η0为螺旋桨敞水效率，-η0为优化目标即敞水效率尽可能大；约束条件中，KT为优化桨推力系数
为原始桨推力系数，ε为约束容差，表示优化后螺旋桨与原始桨推力系数的误差，该约束条件即保证优化桨与原桨的推力系数差值在可接受范围，否则因降低螺旋桨轴向激振力而使螺旋桨推进性能大幅下降也不符合螺旋桨的优化需要。
因此螺旋桨优化的流程为：通过CFD数值方法预报得到螺旋桨敞水推进效率以及轴向一倍频非定常力系数等优化目标，由NSGA Ⅱ算法进行优化计算，根据优化目标及约束条件满足情况，重新设计优化变量，得到新的螺旋桨设计方案再次开始CFD数值计算，得到优化目标后再次进行优化结果评判、重新设计优化变量，直至满足优化迭代要求为止。
3 数值方法及计算模型
3.1 数值方法验证
螺旋桨优化选取原型桨为某七叶侧斜螺旋桨E1619，螺旋桨几何参数如表1所示，几何模型如图3所示。
表1 E1619桨几何参数
Table 1 E1619 propeller geometry parameters
参数数值直径(D/mm)485桨叶数(N)7θ0.9R(°)16.55P/D at r=0.7R1.067 3
图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
选取Nathan[25]论文中的试验数据验证数值模型，图6给出了数值结果与试验值。图中，KT为推力系数
其中T为螺旋桨推力)，KQ为扭矩系数
螺旋桨扭矩)，J为螺旋桨进速系数
为进流速度，n为螺旋桨转速)。从图6中可以看出，数值计算结果与试验值相比偏小，这可能由于数值计算中流场设置为理想流场与真实流场有一定差异，数值结果与试验值误差最大值时J=0.8，此时KQ与试验值误差为3.5%，整体与试验值吻合良好，验证了所采用数值计算方法的准确性。
图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 网格划分参数
Table 2 Meshing parameter settings
螺旋桨直径桨叶表面网格尺寸首层边界层厚度网格数量Mesh a5.0%×D0.5%×D0.1%×D6.53E6Mesh b2.5%×D0.25%×D0.05%×D1.16E7Mesh c1.25%×D0.125%×D0.025%×D2.10E7
按照上述网格划分方案，计算螺旋桨设计工况(J=0.7)时的水动力特性，表3给出了3种网格方案计算得到的螺旋桨水动力参数。可以看出，随着网格的加密，数值计算结果与试验值的差异逐渐减小。采用Mesh b网格时，推力系数与试验值误差为2.16%，敞水效率与试验值误差为0.7%，已经满足计算精度要求，因此后续计算中均采用Mesh b网格进行，桨叶表面网格划分如图7所示。
表3 不同网格划分下KT and η0与试验值参数
Table 3 Comparison of KT and η0 with experimental results
KTΔKT/%η0Δη0/%试验结果0.264—0.628—Mesh a0.2534.20.6102.9Mesh b0.2583.00.6161.9Mesh c0.2601.50.6203.3
图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 优化桨与原桨几何参数
Table 4 Comparison of geometric parameters and performance between the optimized propeller and the original propeller
螺旋桨(P/D)0.7Rθ0.9R(°)KFX×104ΔKFX/%η0Δη0/%KTΔKT/%E16191.067 316.550.978—0.616 3—0.258 5—P2360.961 617.240.880-10.00.621 50.80.242 1-6.3P1770.946 217.350.951-2.80.626 51.70.233 0-9.9P1021.016 216.910.912-6.70.625 01.40.235 6-8.8P420.936 815.870.884-9.60.631 12.40.215 0-16.8
从表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)的侧斜角度是提高推进效率、降低一阶轴向脉动力的研究方向。
[1] 熊远辉,赵道荣,李爱平,等。船舶螺旋桨对白鳍豚等大型水生动物的伤害及防护初探[J].淡水渔业,2004,34(05):15-16.
Xiong Y H,Zhao D R,Li A P,et al.Preliminary study on damage and protection of ship propellers to large aquatic animals such as Baiji[J].Freshwater Fisheries,2004,34(05):15-16.
[2] 蔡昊鹏,马聘,陈科,等.确定螺旋桨侧斜分布的一种数值优化方法[J].船舶力学,2014(07):771-777.
Cai H P,Ma P,Chen K,et al.A numerical optimization method to determine the skew distribution of propellers[J].Ship Mechanics,2014(07):771-777.
[3] 王超，熊鹰，何苗，等.粒子群算法在大侧斜螺旋桨侧斜分布优化中的应用[J].上海交通大学学报,2012,46(03):398-409.
Wang C,Xiong Y,He M,et al.Application of particle swarm optimization in the optimization of skew distribution of large skew propellers[J].Journal of Shanghai Jiao Tong University,2012,46(03):398-409.
[4] ANDERSEN P,KAPPEL J J,SPANGENBERG E.Aspects of propeller developments for a submarine[C]//Norway：First International Symposium on Marine Propulsors,2009.
[5] The propulsion committee final report and recommen-dations to the 26th ITTC[R].Rio de Janeiro,Brazil,Propulsion Committee,ITTC,2011.
[6] Final report and recommendations to the 27th ITTC[R].Copenhagen,Denmark,Propulsion Committee,ITTC,2014.
[7] MISHIMA S,KINNAS S A.Application of a numerical optimization technique to the design of cavitating propellers in non-uniform flow[J].Journal of ship research,1997,41(02):93-107.
[8] BERTETTA D,BRIZZOLARA S,GAGGERO S,et al.CPP propeller cavitation and noise optimization at different pitches with panel code and validation by cavitation tunnel measurements[J].Ocean Engineering,2012,53(53):177-195.
[9] GRIFFIN P E,KINNAS S A.A design method for high speed propulsor blades[J].Journal of Fluids Engineering,Transactions of the ASME.1998,120(03):556-562.
[10] 曾志波,丁恩宝,唐登海.基于BP人工神经网络和遗传算法的船舶螺旋桨优化设计[J].船舶力学,2010,14(Z1):20-27.
Zeng Z B,Ding E B,Tang D H.Optimal design of ship propeller based on BP artificial neural network and genetic algorithm[J].Ship Mechanics,2010,14(Z1):20-27.
[11] 潘光,肖志勇,黄明明,等.对转螺旋桨设计参数优化[J].鱼雷技术,2010,18(01):6-8.
Pan G,Xiao Z Y,Huang M M,et al.Design parameter optimization of counter-rotating propeller[J].Torpedo Technology,2010,18(01):6-8.
[12] 叶礼裕,王超,黄胜,等.不同螺距拟合方式对螺旋桨优化效果的影响分析[J].中国造船,2015,56(03):96-107.
[13] Ye L Y,Wang C,Huang S,et al.Analysis of the influence of different pitch fitting methods on the propeller optimization effect[J].China Shipbuilding,2015,56(03):96-107.
[14] 齐江辉，郭健，郑亚雄，等.七叶侧斜螺旋桨设计参数对空泡性能的影响研究[J].推进技术，2019,40(10):2366-2372.
Qi J H,Guo J,Zheng Y X,et al.Research on the influence of design parameters of seven-blade skew propeller on cavitation performance[J].Propulsion Technology,2019,40(10):2366-2372.
[15] 常欣,孙帅,王超,等.考虑非设计工况的螺旋桨优化设计方法[J].哈尔滨工程大学学报,2015,36(12):1544-1548.
Chang X,Sun,Wang C,et al.Propeller optimization design method considering non-design conditions[J].Journal of Harbin Engineering University,2015,36(12):1544-1548.
[16] 郭春雨,王超,熊鹰,等.两种优化算法在螺旋桨优化设计中的应用与比较[J].船舶工程,2013,35(03):23-26.
Guo C Y,Wang C,Xiong Y,et al.Application and comparison of two optimization algorithms in propeller optimization design[J].Ship Engineering,2013,35(03):23-26.
[17] DAI C,HAMBRIC S,MULVIHILL L,et al.A prototype marine propulsor design tool using artificial intelligence and numerical optimization techniques[J].Transactions Society of Naval Architects and Marine Engineers,1994,102:57-69.
[18] 黄斌，熊鹰，王波.基于粒子群算法的螺旋桨侧斜分布优化[J].中国舰船研究,2016,11(06)：83-89.
Huang B,Xiong Y,Wang B.Optimization of propeller skew distribution based on particle swarm algorithm[J].Chinese Ship Research,2016,11(06):83-89.
[19] VESTING F,BENSOW R.Propeller optimisation considering sheet cavitation and hull interaction[C]//Hamburg:Second International Symposium on Marine Propulsors,2011.
[20] SRINIVAS N，DEB K.Multi-objective function optimization using non-dominated sorting genetic algorith[J].Evolutionary Computation,1995,2(03):221-248.
[21] DEB K，AGRAWAL S，PRATAP A，et al.A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization:NSGA Ⅱ[C]//Proceedings of the 6th International Conference on Parallel Problem Solving from Nature Vi.Paris,Springer,2000:849-858.
[22] 魏春阳，庞之洋，王永生，等.喷水推进器水动力性能分析与优化设计研究[J].船海工程，2014，43(06)：35-39.
Wei C Y,Pang Z Y,Wang Y S,et al.Hydrodynamic performance analysis and optimization design of water jet propulsion[J].Ship and Marine Engineering,2014,43(06):35-39.
[23] Ji B，Luo X W，Wu Y L，et al.Numerical investigation of unsteady cavitating turbulent flow around a full scale marine propeller[J].Journal of Hydrodynamics，2010,22(05):705-710.
[24] 王波,张以良,熊鹰.船用螺旋桨水动力性能优化设计[J].武汉理工大学学报(交通科学与工程版),2015(04):773-777.
Wang B,Zhang Y L,Xiong Y.Hydrodynamic performance optimization design of marine propeller[J].Journal of Wuhan University of Technology (Transportation Science and Engineering Edition),2015(04):773-777.
[25] 饶志强,杨晨俊.七叶侧斜螺旋桨水动力性能优化[J].中国造船,2016,57(04):1-13.
Rao Z Q,Yang C J.Optimization of hydrodynamic performance of seven-blade skew propeller[J].China Shipbuilding,2016,57(04):1-13.
[26] Nathan Chase.Simulations of the DARPA Sub-off submarine including self-propulsion with the E1619 propeller[D].University of Iowa,2012.
[27] 王超.螺旋桨水动力性能、空泡及噪声性能的数值预报研究[D].哈尔滨：哈尔滨工程大学，2010.
Wang C.Research on numerical prediction of propeller hydrodynamic performance,cavitation and noise performance[D].Harbin：Harbin Engineering University,2010.
[28] ROBERT W，JUN S，FRED S.Discussion：Criticisms of the“Correction Factor” Verificatio Method[J].Journal of Fluids Engineering，2004，126(04)：704-706.