张富毅， 吴钦， 赵晓阳， 刘影， 王国玉

(北京理工大学 机械与车辆学院, 北京 100081)

0 引言

喷水推进器具有推进效率高、抗空泡能力强和附体阻力小等特点，广泛应用于高速舰艇和两栖车辆等水面航行体上[1-2]。为了实现水面航行体的机动性、高速性和操纵性，研究作为核心动力来源的喷水推进器显得尤为重要[3-4]。喷水推进器进水流道是将水流从船底引流到喷水推进泵进口的部件，是喷水推进系统的重要组成部分之一。在喷水推进器工作过程中，有7%～9%的输入功率在进水流道内损失，而喷水推进系统和船体结构之间的相互作用对推进效率的影响可以达到20%，其中进水流道对泵船相互作用有着重要影响[5]。因此进水流道的流动性能直接影响到喷水推进器整体的推进性能，优化进水流道是提高喷水推进器整体性能的重要途径之一。

早期对喷水推进器进水流道的研究主要采用模型试验方法，Brandner等[6]对平口式进水流道展开试验研究，基于压力探针观测了进水流道流动的显著特征，包括流动分离、空泡产生等。Robert等[7]采用平口式喷水推进装置研究进水流道的吸水效应，试验结果表明，针对这种形式的喷水推进器，惯用的设计方法可能会导致对推力的明显低估。吴民权等[8]采用风洞试验研究船体- 喷水推进器的相互作用，试验结果表明进口处的阻力随着喷速比的增大而减小，并且采用边界层抽吸方法可以减小进口处大尺度旋涡区，有利于改善泵进口截面的速度分布。

随着高性能计算机和计算流体力学(CFD)的快速发展，数值模拟方法被广泛地应用到喷水推进器流场分析评估及性能预报等研究中，并与模型试验方法结合来证明方法的可靠性。许慧丽等[9]对不同船舶航行姿态下的喷水推进器进行了模型试验和数值模拟，结果表明计算结果和试验数据吻合，与直航和斜航相比，倒航状态下流场更不均匀，压力更低，更易发生空化和流动分离。靳栓宝等[10]对某轴流式喷水推进器进行改型设计且进行实船试航试验，结果表明新设计的喷水推进器推进航速超设计航速9.4%，数值预报航速与试航结果误差为1.5%.

在此基础上，许多学者针对喷水推进器进水流道影响因素展开进一步研究。魏应三等[11]针对某平口式进水流道，在纵向总长度、宽度、高度给定的条件下，利用CFD技术建立了6种不同倾斜角的进水流道模型，研究倾斜角与流动性能之间的关系，结果表明倾斜角为40°时流动性能最优。丁江明等[12]基于进水流道11个动态关联的几何参数，提出了进水流道参数化设计方法，结果表明随着流道倾斜角的增加，流道背部产生明显的流动分离。吴娜等[13]通过分析进水流道效率、出流口加权平均角、不均匀系数和唇部驻点位置等参数，建立优选方案进行数值模拟，结果表明减小流道倾斜角可以改善流道出流均匀性，圆弧唇角的流动性能优于抛物线唇角，且唇角形状应适当尖锐。

本文针对喷水推进器进水流道的内部流场进行数值计算，并采用试验设计、响应面拟合以及多目标遗传算法对进水流道设计变量进行优化研究，分析不同进水流道设计变量对流道性能的影响，为喷水推进器进水流道快速优化设计提供一种高效的方法。

1 数值计算方法

1.1 基本控制方程

采用雷诺时均方法求解连续性方程和动量方程：

(1)

(2)

式中：i，j=1，2，3；xi、xj为坐标轴方向；ui、uj为速度分量；ρ为流体密度；p为压力；μ、μt分别为层流和紊流黏性系数。

1.2 湍流模型

采用剪切应力传输(SST)k-ω(k为湍动能，ω为湍流频率)湍流模型，该模型集合了k-ε(ε为湍动能耗散)模型和k-ω模型的优点：在近壁区域采用k-ω模型，湍流耗散率小，收敛性好；在湍流充分发展区域采用k-ε模型，计算效率高，对复杂流场的适应性更好[14]。湍动能k方程和湍流频率ω方程分别为

(3)

(4)

(5)

式中：Pk、Pω为湍流生成项；Dk为湍流耗散项；σk和σω2分别为湍动能k和湍流频率ω的普朗特数；F1、F2为混合函数；S为剪切应变率；Cω、βω和a1为模型常数，取值见文献[15]。

1.3 数值计算模型

计算采用平口式进水流道，进水口采用矩形+椭圆形，进水流道的轴面结构如图1所示，由入口过渡段、倾斜直管段、圆弧过渡段和水平直管段组成，各部分几何参数如表1所示。

图1 进水流道轴面结构参数示意图Fig.1 Structural parameters of axial plane of inlet duct

表1 进水流道几何参数Tab.1 Geometric parameters of inlet duct

由于进水流道的性能与船体结构以及航行条件密切相关，因此在对进水流道进行数值模拟时应充分考虑船底尾部进水口周围的流场区域。参考文献[16]，船底流场控制体的长宽高分别为20D、10D和8D，如图2所示。流场来流面设置为速度进口边界条件，出流面和流道出口设置为压力出口边界条件，流道壁面、泵轴和船底设置为无滑移壁面，控制域的两个侧面和底面设置为开放边界条件。

图2 计算区域及边界条件Fig.2 Computational domain and boundary conditions

进水流道以及船底控制域均采用非结构化网格，在进水流道以及进水口附近的船底区域进行局部网格加密，在流道和泵轴近壁区采用边界层网格。整个计算域网格数量为140万，如图3所示。动量方程、湍流动能方程和耗散率方程均采用2阶迎风格式，收敛精度为10-5.

图3 网格划分Fig.3 Mesh generation of water-jet inlet duct

2 喷水推进器进水流道的优化分析

本文采用响应面方法对喷水推进器进水流道进行参数设计及优化。响应面法是一种基于试验设计理论对指定的设计点集合进行试验，得到目标函数和约束函数的响应面模型，来预测非试验点响应值的方法[17-18]。响应面方法优化流程如图4所示，包括问题设置、试验设计、响应面拟合和Pareto寻优等过程。

图4 响应面法优化流程图Fig.4 Flow chart of response surface optimization

由图1可知，进水流道出口直径D由推进泵直径决定，流道高度H、进水口长度LIN和流道总长度L是由船体整体尺寸确定，结合文献[10]以及初步理论分析，倾斜角α和过渡段半径RT对流动性能影响较大，因此本文以倾斜角α和过渡段半径RT为设计变量，分析二者对流动性能的影响。设计变量取值范围如表2所示。

表2 设计变量取值范围Tab.2 Value ranges of design variables

为了定量评估进水流道流动性能，选取流道效率η、出口速度加权平均角θ和出口速度不均匀系数ξ为目标变量[13,19]，分别定义为

(6)

(7)

(8)

对于以上目标变量，流道效率η越高、出口速度加权平均角θ越接近90°、出口速度不均匀系数ξ越小，进水流道流动性能越好。

2.1 试验设计

试验设计是在设计变量取值范围内生成试验的样本空间，取样时应选取具有代表性的设计点，以尽可能少的点最大限度地反映出样本空间数值特点。本文采用中心复合设计(CCD)方法[20]，如图5所示，试验点由以下3部分组成：

1) 一个样本空间中心点(α，RT)=(40°，450 mm)；

2) 2×2个样本空间轴向点(α，RT)=(30°，450 mm)、(50°，450 mm)、(40°，300 mm)、(40°，600 mm)；

3) 22个样本空间析因点(α，RT)=(30°，300 mm)、(30°，600 mm)、(50°，300 mm)、(50°，600 mm)。

图5 样本空间示意图Fig.5 Sample space

为了提高响应面模型的精度，再结合拉丁超立(LHSD)方法随机生成均匀覆盖整个样本空间的26个样本点，如图5中的圆点，共计35个样本点。

2.2 响应面拟合

响应面拟合是以试验设计中的设计点为基础，通过一定的函数关系拟合出目标变量与设计变量之间的关系。本文采用2阶多项式拟合，具体表达[21]如下：

(9)

为了对拟合响应面精度进行定量评估，研究基于决定系数R2和均方根差σRMSE来评价响应面模型的预测能力[18]：

(10)

(11)

图6给出了基于样本空间设计点的2阶多项式响应面拟合优度图，其中横坐标是设计点目标变量数值计算结果，纵坐标是设计点目标变量响应面预测值，可以看出，不同目标变量的响应面预测值随设计点观测值的变化基本呈现出斜率为1的线性变化，说明响应面拟合精度较高。流道效率η、出口速度加权平均角θ和出口速度不均匀系数ξ的决定系数均大于0.96，均方根差均小于0.1，如表3所示，因此，响应面拟合精度符合要求。

图6 拟合优度图Fig.6 Goodness of fit

表3 拟合度评价结果Tab.3 Evaluated results of fit degree

图7通过比较自变量的偏方差和总方差分别给出了各自变量对目标变量的局部灵敏度和全局灵敏度，即自变量α和RT对目标变量η、θ和ξ的相对重要程度。根据文献[22]，灵敏度指标定义为

(12)

(13)

图7 自变量对目标变量局部、全局灵敏度图Fig.7 Distribution of local and global sensitivities

对3个目标变量η、θ和ξ的灵敏度分析中发现，倾斜角α对目标变量流道效率η、出口速度加权平均角θ和出口速度不均匀系数ξ的全局灵敏度系数分别为63%、78%和65%，均大于过渡段半径RT对目标变量的全局灵敏度系数，说明目标变量受倾斜角α的影响较大，即倾斜角α对进水流道的流动性能具有较大影响。

图8分别给出了流道效率η、出口速度加权平均角θ和出口速度不均匀系数ξ随倾斜角α和过渡段半径RT的变化趋势图。由图8(a)图可以看出，流道效率η随着流道倾斜角α的增加而降低，而流道效率η随着过度段半径RT的变化基本不变。同理，由图8(b)可以看出，出口速度加权平均角θ随着倾斜角α的增大有减小的趋势。由图8(c)可以看出，当α为30°～40°时，出口速度不均匀系数ξ随着倾斜角α的增加显著增大，当α为40°～50°时，出口速度不均匀系数ξ随着倾斜角α的增大无明显变化，维持在0.15左右。

图8 目标变量随自变量变化趋势图Fig.8 Variation of object variable with design variable

2.3 多目标优化

为了提高流道效率η、出口速度加权平均角θ最接近90°和减小出口速度不均匀系数ξ，建立进水流道的多目标优化数学模型如下：

(14)

式中：Xi为设计变量；Xmin为设计变量下限值；Xmax为设计变量上限值。图9给出了目标变量相互之间的变化趋势：如图9(a)所示，当流道效率η从0.91增加到0.92时，出口速度加权平均角θ由85°减小到83°；同理，由图9(b)可知，当出口速度加权平均角θ从82°增加到85°时，出口速度不均匀系数ξ由0.04增加到0.06；由图9(c)可知，当出口速度不均匀系数ξ从0.08减小到0.055时，流道效率η由0.93减小到0.92. 因此，在多目标优化中，无法实现流道效率η、出口速度加权平均角θ和出口速度不均匀系数ξ同时达到最优值，所有可能解构成的解集为Pareto最优解集，它是由那些任一个目标变量的提高都必须以牺牲其他目标变量为代价的解组成的集合[23]。

图9 目标变量之间相互变化趋势图Fig.9 Variations of object variables

为了找出目标变量之间的最佳平衡点，本文采用多目标遗传算法(MOGA)对(14)式的数学模型进行优化[24]。图10为MOGA的基本优化流程，采用偏移哈默斯利抽样(SHS)技术[25选取10N=20个样本点为初始种群，其中N为自变量个数，计算流道效率η、出口速度加权平均角θ和出口速度不均匀系数ξ，若满足终止准则则产生终止种群，若不满足则进行适值分配、遗传操作和插入父代等过程得到下一代种群，直至满足终止准则。

图10 MOGA流程图Fig.10 Flow chart of MOGA

为了判断遗传算法是否收敛，引入最大允许Pareto百分比(MAPP)和收敛稳定性百分比(CSP)。Pareto百分比表示每次迭代样本数中Pareto点所占的百分比，当Pareto百分比达到MAPP(99%)时，算法收敛；稳定性百分比是基于种群均值和均方差表示全局稳定性的参数，当稳定性百分比小于CSP(0.01%)时，算法收敛，数学表达式为

(15)

(16)

式中：S为收敛稳定性百分比；Yi和Yi-1分别为第i和第i-1代种群平均值；σi和σi-1分别为第i和第i-1代种群均方差；Ymax和Ymin分别为初始种群的最大值和最小值。图11为Pareto百分比和稳定性百分比随迭代次数的变化趋势，可以看出，当迭代36次后，Pareto百分比为100%，大于99%，稳定性百分比为0.03%，趋近0.01%，满足收敛条件。

图11 收敛准则随迭代次数变化趋势图Fig.11 Convergence criterion versus iterations

图12给出了Pareto解集的分布图，不同颜色的点表示不同Pareto向前指数(PFI)，PFI越小表示Pareto解越符合要求，可以看出，设计空间内解的分布主要集中在Pareto前沿附近(图12中红色区域)，其余区域相对比较稀疏，说明优化过程中设计变量α和RT逐渐向最优目标变量η、θ和ξ靠拢。择优选出的3个Pareto最优解如表4所示，其中流道效率η均高于0.94，出口速度加权平均角θ均大于87°，出口速度不均匀系数ξ均小于0.067.

图12 Pareto解集分布图Fig.12 Distribution of Pareto solution sets

表4 Pareto最优解Tab.4 Pareto optimal solutions

2.4 优化结构分析

为了验证进水流道优化结果的准确性，对优化点1进行数值计算。表5为进水流道优化前后目标变量数值计算结果的对比，优化后的流道效率η、出口速度加权平均角θ和出口速度不均匀系数ξ都优于原有流道，表明优化后进水流道的流动性能优于原有流道。通过对比分析表4和表5优化点1目标变量的数值，表明基于多目标遗传算法的目标变量预测值与数值计算的目标变量结果误差较小，优化结果可靠。

表5 优化前后流动性能参数对比Tab.5 Comparison of flow performance parametersbefore and after optimization

图13 进水流道三维流线图Fig.13 3D streamline of inlet duct

图13为优化前后的进水流道三维流线图，可以看出：优化前进水流道三维流线存在明显扰动，尤其是泵轴上方，存在明显的周向速度；优化后，泵轴上方无明显扰动，周向速度基本消除，因此提高了出口速度加权平均角θ，进而提高了流道效率η.

图14对比了优化前后进水流道轴面压力图，可以看出：优化前在流道背部和唇部存在明显的低压区，流道的抗空化性能较差，泵轴上方压力分布不均匀，存在较大的压力梯度；优化后流道背部和唇部的低压区基本消除，轴面的压力分布更加均匀，从而提高了流道效率η.

图14 进水流道轴面压力分布图Fig.14 Pressure contour on axial plane of inlet duct

图15对比了优化前后进水流道出口截面速度分布图，可以看出:优化前出口截面在泵轴上方存在明显的低速区，局部区域速度仅1～3 m/s，其余区域速度保持较大值，速度分布不均匀，出口速度不均匀系数ξ较大；优化后出口截面的低速区明显减小，速度分布更加均匀，出口速度不均匀系数ξ较小。

图15 进水流道出口截面速度分布图Fig.15 Velocity contour of outlet section of inlet duct

3 结论

本文基于响应面方法对喷水推进器进水流道的内部流场进行优化设计，通过试验设计、响应面拟合以及多目标遗传算法对进水流道设计变量进行优化研究，结论如下：

1) 建立的关于流道效率η、出口速度加权平均角θ和出口速度不均匀系数ξ随倾斜角α和过度段半径RT变化的响应面模型具有较高的拟合精度，响应面模型预测值与仿真计算值的相对误差在合理范围之内。

2) 流道效率η、出口速度加权平均角θ和出口速度不均匀系数ξ受倾斜角α的影响均大于过渡段半径RT，优化设计过程中需重点对倾斜角α进行优化。流道效率η和出口速度加权平均角θ随着倾斜角α的增加而减小； 当α为30°～40°时，出口速度不均匀系数ξ随着倾斜角α的增加显著增大，当α为40°～50°时，出口速度不均匀系数ξ随着倾斜角α的增大无明显变化，维持在0.15左右。

3) 基于多目标遗传算法的优化方法可以实现进水流道的快速优化，优化后的进水流道较优化前的进水流道相比，流道效率、出口速度加权平均角以及出口速度不均匀系数都有明显的提高。