基于NFFD 和高斯近似的螺旋桨多目标优化

2022-11-26周喜宁朱耀龙

舰船科学技术 2022年19期

刘旭，周喜宁，朱耀龙

(1.渤海船舶职业学院，辽宁葫芦岛 125101；2.招商局重工(江苏)有限公司，江苏南通 226000；3.华中科技大学，湖北武汉 423201)

0 引言

螺旋桨设计及优化一直是研究船舶快速化、经济性的紧迫任务之一。但由于螺旋桨参数众多、建模繁琐、水动力计算周期长等问题[1]，长期限制船舶设计初期对螺旋桨的高效优化设计。建立一套包括螺旋桨几何变形和重构、水动力性能预测近似模型以及多目标优化的螺旋桨构型优化设计方法[2]，对促进螺旋桨高效优化设计具有重要意义。

螺旋桨的快速建模、高效优化始终是国内外学者关注的热点研究。Ramin 等[3]采用Non-Gradinet-based Algorithm 和梯度法分别对弦长、厚度以及螺距沿半径的分布进行数值优化；Florian 等[4]基于B-spline 对螺距、拱度、厚度、弦长、曲率和倾角沿半径分布进行拟合，并以最大效率和最小脉动压力为目标进行了螺旋桨优化。叶礼裕[5]结合定常面元法和克里格近似理论建立螺旋桨参数与性能特征的预测模型，并利用粒子群算法对螺旋桨参数进行了多目标优化设计。王超[6]用试验设计方法为神经网络近似模型的建立提供足量且高质量的训练样本，利用遗传算法可以获得更广解集范围更广，而且耗时间大幅度缩减。马丹萍[2]借助于ISIGHT 优化平台，集成螺旋桨的FFD 方法和CFD水动力数值评估，提出了一套从构型-数值求解-优化的螺旋桨设计优化方法。

针对目前螺旋桨研究中存在的螺旋桨模型构造周期长、水动力性能评估耗时多的问题，本文结合基于NURBS 基函数的自由变形算法（NURBS-based freefrom deformation，NFFD）、多输出高斯近似模型以及NSGA-Ⅱ多目标优化算法，构建一套包括螺旋桨变形重构—水动力性能快速预测—多目标优化的螺旋桨高效自动优化方法，即以NFFD 技术实现螺旋桨三维模型的变形与重构，采用有限元数值仿真得出样本螺旋桨的水动力性能，基于多输出高斯近似理论建立螺旋桨水动力性能预测近似模型，结合近似模型和NSGA-Ⅱ算法对KP505 桨进行多目标优化设计，验证该方法的可行性。

1 NFFD 方法及螺旋桨变形控制

基于Sederberg 等[7]提出的自由曲面变形（free form deformation，FFD）方法，研究人员将FFD 技术扩展为很多不同形式[8–9]，基本包括：构建参数控制体；将变形物嵌入参数控制体；参数控制体变形；计算变形物形变后的几何坐标[10]。基于FFD 方法，Lamousin 等[8]提出NFFD 算法，将顶点加权与非均匀框架结合，引入权因子和节点作为变量，增加了变形操作的自由度，从而达到更好的变形。马丹萍[2]、张人会[10]、李冬琴[11]、WU[12]、马明生[13]、唐静[14]已将NFFD 技术应用于船型与螺旋桨变形控制与重构，并在船型与螺旋桨优化设计方面取得了丰硕成果。

螺旋桨上任意一点绝对坐标可以用局部坐标系(T,U,V)表示为：

式中：Q0为局部部坐标原点；(t,u,v)为点Q的局部坐标；(T,U,V)为局部坐标系标准方向向量；即控制体基本坐标系。螺旋桨曲面参数化控制节点Pi,j,k可表示为：

式中：l，m，n分别表示局部坐标在(T,U,V)3 个方向的等分份额，i=0,1,···l；j=0,1,···m；k=0,1,···n。

更进一步，螺旋桨上任一点绝对坐标可以用局部坐标表示为：

式中：（v，u，t）为局部坐标；Pi,j,k为参数化控制节点；Bil(t)，Bjm(u)，Bkn(v)为Bernstein 基函数，表示为：

NFFD 方法利用NURBS 基函数替换Bernstein 函数，螺旋桨任意一点Q坐标可表示为：

式中，Wi,j,k为参数化控制节点Pi,j,k的权重因子。

在式（5）基础上采用deBoor-Cox 递推定义[15]和Deboor 算法[16]推导得出NURBS 体三次矩阵表示方法：

最后将螺旋桨嵌入以上参数控制体，并对采样点进行归一化处理，利用牛顿法逆向搜索非线性方程组的最优解，即可得到螺旋桨采样点局部坐标。

在建立螺旋桨与参数体的相互关系后，保持每个采样点在参数体中相应的局部坐标(t,u,v)不变，通过改变参数化控制节点Pi,j,k或其权重因子Wi,j,k即可得到新的带权控制点矩阵，代入式（6）中即可求解得出螺旋桨变形后的局部坐标。

以KP505 桨为例，利用Matlab 实现了螺旋桨几何形变与重构，初始输入螺旋桨表面离散点集的坐标值，并选取0.2R，0.4R，0.6R，0.7R和0.9R处的叶宽、厚度、螺距比、曲率、倾角、最大厚度位置共30 个参数作为螺旋桨形状控制设计变量，输出螺旋桨变形后的表面离散点集的坐标值。KP505 桨形状控制效果和网格模型如图1 所示。

图1 KP505 桨形状控制效果和网格模型Fig.1 KP505 propeller shape control effect and mesh model

2 性能预测高斯近似模型

2.1 多输出高斯近似模型

近似模型是基于不同模型函数对设计变量和设计目标组成的样本数据进行近似拟合，具有计算成本低、函数关系明确等特点[17]。常用的近似模型方法包括：响应面法（RSM）、克里格法（Kringing）和人工神经网络技术[18]。高斯近似是基于Kringing 和贝叶斯定理的一种机器学习方法，对少量样本的高阶非线性问题具有良好的适应性[19–22]。

给定训练样本：

其中，任意一组样本(xi,yi)是随机且相互独立的，通过近似逼近，输出向量yi∈Y⊂R与输入向量xi∈X⊂RP之间存在一些未知函数关系：

式中：Kn(X,X)为n×n阶对称半正定协方差矩阵，该矩阵的元素Kij=k(xp,Xq)表示xp与Xq之间的相关性；I为n×n阶单位矩阵。

多输出高斯近似模型协方差矩阵K表示为：

式中，Kf和Kx分别表示输出和输入向量的相关性。

训练样本的输出向量yi和测试点xtest对应的目标输出f(xtest)服从如下联合先验分布：

式中，k(xtest,X)和k(X,x)分别为n×1 和1×n 阶向量。在得到y和f(xtest)之间的联合先验概率分布之后，可以分析得到f(xtest)的后验概率分布如下式：

本研究中的输入向量X和输出向量Y分别表示为：

式中：a～f分别表示某一剖面处的叶宽、厚度、螺距比、曲率、倾角、最大厚度位置；i=0.2R，0.4R，0.6R，0.7R和0.9R，KT，10KQ，η0以及Pmin分别表示螺旋桨推力系数、扭矩系数、效率以及最小空化压力。

2.2 敞水性能数值仿真

以STAR-CCM+13.0 流体动力学软件对KP505桨的敞水性能进行仿真计算。利用NFFD 建立的原始螺旋桨曲面生成实体模型，构建非结构化六面体网格，采用k-ωSST 模型，并基于MRF 方法求解雷诺平均Navier-Stokes 方程计算不同进速下的螺旋桨性能。

以进速系数J=0.893 为例，在计算单元达到26 000 000时，螺旋桨的端涡和根涡现象如图2 所示，即在此时实现了网格收敛。

图2 高雷诺数下螺旋桨的端涡和根涡现象Fig.2 End vortex and root vortex phenomenon of propeller at high Reynolds number

图3 为进速系数J=0.2～1.1 时，螺旋桨推力系统、仿真计算结果与试验结果对比，仿真与试验所得推力系数KT、扭矩系数10KQ、效率η0的均方差分别为1.63%，2.41%，0.63%，说明以该软件和网格数量计算的螺旋桨敞水性能精度可信。

图3 螺旋桨仿真计算结果与试验结果对比Fig.3 Comparison between the simulation calculation results of the propeller and the test results

为给高斯近似模型提供初始样本，以0.2R，0.4R，0.6R，0.7R和0.9R处的叶宽、厚度、螺距、曲率、倾角、最大厚度位置共30 个参数作为螺旋桨形状控制设计变量，采用NFFD 技术生成100 个螺旋桨计算样本，并计算进速系数J=0.893 时KT，10KQ，η0以及Pmin。

2.3 近似模型的验证

对100 个初始样本采用交叉验证方法对训练后的多输出高斯代理模型进行评估。通过交叉验证得到的100 个测试样本预测误差如图4 所示。

图4 高斯近似模型的交叉验证预测误差Fig.4 Cross-validation prediction error of Gaussian approximation model

结果表明，各测试样本KT，10KQ，η0以及Pmin预测值与数值模拟结果的相对误差在3%以内，精度满足要求。训练后的多输出高斯代理模型可用于多目标优化中目标函数值的评估。

3 螺旋桨多目标优化

3.1 NSGA-Ⅱ多目标优化

NSGA-Ⅱ是由Kalyanmoy 等[23]针对NSGA 提出的一种改进进化算法，该算法具有多目标适度简化、最优解拓展、求解精度高、高鲁棒性等优点，尤其针对少量样本、多变量多目标优化问题就较快的求解速度。

螺旋桨敞水性能多目标优化问题可以描述为：

式中：叶宽0≤xa≤0.27，叶片厚度0.08 ≤xb≤0.095，螺距比1.2 ≤xc≤2.9，曲率−0.01 ≤xd≤0.02，倾角−0.5 ≤xe≤0.3，最大厚度位置0.45xb≤xf≤0.75xb，Ae螺旋桨展开面积，Ad面螺旋桨盘面积。

利用以100 个样本建立的多输出高斯近似模型计算螺旋桨敞水性能，并以NSGA-Ⅱ算法进行螺旋桨的优化设计，将种群交叉遗传概率设置为0.92，变异概率设置为0.08。经过1 000 代的演化，可得到Pareto 非支配最优解集。

初始种群和Pareto 解集的扭矩系数和效率如图5所示。可以看出最优解是在原始设计空间之外的，其中Pareto 前沿上A点为综合最优设计，表示Pareto 前沿中的最优综合性能点，同时兼顾了扭矩系数和效率。

图5 初始种群和Pareto 解集扭矩系数和效率Fig.5 Torque coefficient and efficiency of initial population and Pareto solution set

3.2 优化设计结果及分析

选取优化结果中A点进行根据优化结果分析，原模型的效率提高了2.6%、扭矩降低了3.3%。此外，盘比降低了2.5%，但螺旋桨模型的空化特性并未恶化，因为推进器旋转区域的最小压力增加了7%。优化方案4 桨叶r/R=0.7 处剖面轮廓如图6 所示。可以看出，经过优化，叶片的厚度和宽度发生了明显的变化，在相对半径0.7 处变小，但倾角有所增加。

图6 优化方案4 桨叶r/R=0.7 处剖面轮廓Fig.6 Section profile at r/R=0.7 of the optimized scheme 4 blade

桨叶宽度、剖面厚度、倾斜角沿径向的分布如图7所示。可以看出，经过优化，桨叶厚度、宽度、倾斜角在相对半径r/R=0.7 处明显变小，而曲率在此处所有增加。

图7 KP505 模型桨叶宽度、厚度和倾斜角沿半径分布Fig.7 KP505 model blade width,thickness and inclination angle distribution along the radius

优化前后螺距比和曲率沿径向的分布如图8 所示。优化前后螺旋桨的螺距比趋势相似，而且各半径处螺距比较对相近。此外，优化结果表明，螺旋桨螺距比在靠近桨叶根部较大。从推进设计的角度来看，螺旋桨螺距比在靠近桨叶根部的逐渐增大可以增加推力，从而提高效率。但是，通过增加螺距比，需要提高螺旋桨轴功率以及螺旋桨流动区域的最小压力。优化前后螺旋桨曲率沿半径分布特征相似，优化后的螺旋桨桨叶根部曲率变为负值，而且优化后的曲率变小，有利于降低螺旋桨噪声。