跨介质飞行器头部外形优化及入水性能分析

2022-12-01马文朝孟繁敏马诺孟军辉邹汝平

兵工学报 2022年10期

马文朝，孟繁敏，马诺，孟军辉,2，邹汝平

(1.北京理工大学宇航学院，北京 100081；2.飞行器动力学与控制教育部重点试验室，北京 100081；3.中国兵器工业集团公司西安现代控制技术研究所，陕西西安 710003)

0 引言

跨介质飞行器入水性能的研究对于解决鱼雷、导弹入水，航行器水上起降等实际问题的运动特性与结构分析具有重要的工程意义[1]。跨介质飞行器高速入水时，头部外形影响其与水面的接触角度，使空泡的形成、发展和溃灭过程发生改变，进而影响其运动学特性[2]；头部外形同样影响冲击瞬间与水面的接触面积，进而影响其所受冲击应力。如何在满足入水速度等约束条件下，实现跨介质飞行器最佳的头部外形设计对于减小跨介质飞行器入水阻力和冲击应力至关重要。

近年来，国内外学者对不同跨介质飞行器入水问题进行仿真和试验研究。May[3]对垂直入水问题进行试验研究，讨论了钢球的速度、直径大小与旋成体头部外形对入水过程中空泡产生、发展和闭合的影响。Bodily等[4]通过试验探究了圆锥、椭圆、扁平头部外形细长轴弹体的受力、入水轨迹和速度，结果表明圆锥头部外形的弹体所受阻力最小。张伟等[5]通过试验对比分析速度在35～160 m/s的平头部外形、卵型和截卵型弹体对入水过程中弹体弹道稳定性的影响；罗驭川等[6]通过试验研究截锥体头部外形的旋成体在低速入水时，头部直径大小对入水过程中空泡、旋成体的速度和俯仰角的影响；施红辉等[7]试验研究了4种不同头部外形的钝体在复杂工况下的入水过程，分析了头部外形和入水速度对弹体运动特性及空泡发展规律的影响。Shi等[8]通过试验与仿真，对比分析了入水参数、飞行器壳体厚度与头部形状等参数对跨介质飞行器入水过程中头部变形、应力分布情况的影响。国内外研究现状表明，跨介质飞行器入水过程的性能研究牵涉到气/液多相流体和跨介质飞行器结构的复杂耦合关系，特别是在触水阶段，头部外形对跨介质飞行器的入水速度与所受最大冲击应力有重要影响。由于对复杂多相流耦合问题仿真研究精度不足，且计算成本较高，学者更倾向于通过试验定性分析入水参数对入水性能的影响，同时对基于介质跨越阶段的头部外形多学科设计优化研究较少。

多学科设计优化的理论和方法的研究国内起步较晚[9]，近年来随着多学科优化设计方法的发展，学者针对回转体优化设计逐渐展开研究。宋保维等[10]利用现代流体力学计算技术，优化得出回转体最小阻力外形；李福新等[11]对双参数椭圆型曲线头部外形鱼雷进行优化设计，优化结果的流噪声明显小于初始结果；余德海等[12]建立鱼雷外形一体化优化设计模型，进行多学科多目标优化设计，综合考虑流体动力与声学两个学科的耦合关系，优化后鱼雷性能有较大提高；赵加鹏等[13]提出了一种多学科优化算法，适用于综合优化设计鱼雷外形，通过实例计算，获得了一种综合性能较优的方案；权晓波等[14]采用代理模型技术对水下航行体头部外形进行优化设计，得到具有良好水动特性的航行体头部外形；Luo等[15]综合考虑水下航行器的航行阻力与能耗，通过仿真与动态代理模型对水下航行器外形进行优化设计，阻力和能耗的优化结果证明了优化方案的有效性；现有文献多为在单一介质条件下对回转体进行多学科设计优化。以跨介质飞行器入水性能为目标的头部外形优化，需要在跨介质飞行器入水性能仿真分析的基础上，进行反复多次的迭代计算，有限元仿真与迭代计算需要消耗大量的计算成本，应在保证计算精度的前提下尽量减小计算成本，采用代理模型技术可以有效地节约时间，提升计算效率[15-16]。

本文在现有跨介质飞行器入水研究的基础上，用拉丁超立方采样法构建跨介质飞行器样本点，利用任意格拉朗日-欧拉(ALE)方法对入水过程进行多相流固耦合仿真，分析入水所受冲击应力与末速度；为降低优化过程中反复迭代有限元仿真所消耗的计算成本，构建径向基函数(RBF)代理模型；同时针对跨介质飞行器头部外形优化这一典型的多目标优化问题，采用多岛遗传算法(MIGA)与自适应优化策略，以避免寻优过程中收敛到局部最优解。以触水阶段跨介质飞行器最大末速度与所受冲击应力峰值最小为优化目标，获取跨介质飞行器头部外形曲线参数最优解，对最佳头部外形进行入水性能分析，相对于常见飞行器头部外形，有效地减少入水过程跨介质飞行器的能量损失，为跨介质飞行器的头部外形设计提供参考。

1 跨介质飞行器入水数学模型

本文采用常见的双参数立方多项式曲线方程对跨介质飞行器头部外形进行参数化建模，跨介质飞行器以60 m/s的速度垂直入水，采用ALE方法分析跨介质飞行器入水后5 ms时速度vend与入水过程中所受最大冲击应力Fσ，为后续代理模型的建立提供样本数据。

1.1 参数化建模

为简化计算，对跨介质飞行器外形参数进行无量纲化处理。跨介质飞行器参数化建模示意图如图1所示，按照回转半径R为1、总长L为20进行设计，跨介质飞行器头部外形曲线表达式采用常见的双参数立方多项式形式[14]，如(1)式所示进行参数化建模：

(1)

式中：X0为头部外所占长度，X0∈[0.5，4]；k1为头部外形前端的曲率变化率，k1∈[0，1]；k2为头部外形后端的曲率变化率，k2∈[0，15]，k1和k2为曲线的控制参数。对跨介质飞行器头部外形影响规律如图2所示。由图2可知，随着k1的减小，头部外形前端曲率增大，而随着k2的增大，头部外形后端曲率增大。

1.2 有限元模型

为获取跨介质飞行器入水性能，对其入水过程进行有限元仿真。本文假设空气和水是不可压缩流体，选用ALE耦合算法求解Navier-Stokes方程[17]，水和空气采用Euler公式求解，跨介质飞行器与自由界面采用Lagrange公式求解，ALE耦合算法的控制方程由下列3个守恒方程给出：

1)质量守恒方程

(2)

2)动量守恒方程

(3)

3)能量守恒方程

(4)

应力张量σij为

(5)

式中：ρ为流体密度；t为时间；wi为相对速度，wi=v-u，v为物质速度，u为网格速度；bi为单位体积力；E为比内能；vi、vj为i轴方向与j轴方向上的速度；xi、xj为ALE坐标；p为压力；δij表示Kronecker函数；μ为动黏性系数。

方程与下列边界条件联立求解：

(6)

式中：

(7)

求解ALE方程时对Lagrange结构进行约束，将结构的相关参量传递给流体单元，本文采用罚函数实现流固耦合算法的约束，通过罚函数耦合系数追踪跨介质飞行器的结构网格(Lagrange)节点和流体物质网格(Euler)位置间的相对位移d，每个计算时间步首先检查结构网格节点是否穿透流体物质表面，若穿透，则在该节点与流体物质表面间、流体节点与结构表面间添加较大的界面接触力，限制穿透，若无穿透则不进行处理。其中界面接触力F为

F=-k·d

(8)

式中：k为基于节点质量模型特性的刚度系数。

由于跨介质飞行器自身结构具有轴对称特性，在垂直入水过程中，所受阻力分布、结构应力分布、入水空泡等同样具有轴对称特性，为进一步减少计算量，基于(1)式建立四分之一垂直入水仿真模型，在对称面上进行节点运动方向约束，边界条件设置为无反射边界条件，模拟无限水域，使结果更为精确。有限元模型如图3所示。

ALE方法结合了Lagrange方法与Euler方法的特点，处理物质边界运动时，ALE方法采用Lagrange方法，能精确跟踪物质边界的运动状态；处理流体等大变形网格时，ALE方法借鉴Euler方法的优势，使网格与材料相互独立，能够处理大变形问题[18]。本文计算域内有限元模型均采用六面体单元进行分析。空气域大小为250 mm×250 mm×500 mm，水域大小为250 mm×250 mm×500 mm。水与空气采用Euler网格，水的材料参数参考汪振等[19]的设置，具体材料参数如表1所示。

表1 空气与水的材料参数

水的状态方程采用Gruneisen状态方程：

(γ0+αμV)E

(9)

式中：V为变化的相对体积；α为对常数γ0的1阶体积修正；μV为体积变化率；C为vsw-vp曲线的截距，vsw为冲击波速度，vp为质点速度；S1、S2和S3为曲线斜率系数。质点速度vp与冲击波速度vsw的关系方程为

(10)

空气的状态方程采用Linear_Polynomial状态方程：

(11)

式中：μρ为密度变化率；c0～c6为材料拟合系数，假设空气为理想气体，c4=c5=0.4，其余为0。

跨介质飞行器采用Lagrange方法，材料模型选择弹塑性材料，材料参数如表2所示，跨介质飞行器直径为50 mm，径向长度为500 mm。

表2 跨介质飞行器的材料参数

1.3 载荷与应力理论分析

跨介质飞行器高速入水瞬间属于碰撞现象，会产生高能量冲击，水初期可忽略所受气体摩擦力与水的浮力[20]，假定跨介质飞行器仅受重力G与头部所受水对其的作用力Fd：

(12)

式中：m为跨介质飞行器的质量；g为重力加速度；Cd为阻力系数；ρw为流体密度；S为跨介质飞行器入水部分在速度方向上的横截面积；vs为跨介质飞行器的速度。

由牛顿第二定律，跨介质飞行器入水初期所受加速度a为

(13)

应力用于测量跨介质飞行器在承受冲击载荷时的内力大小，von Mises应力是根据第四强度理论获得的等效应力，在考虑主应力和剪应力效应同时，可以清楚地描述整个模型中应力的变化，应力大的部位更容易发生结构的损坏失效，设计时应考虑降低跨介质飞行器实际入水过程中所受应力，von Mises应力的表达式如下：

σ=

(14)

式中：σx、σy、σz分别为x、y、z轴方向上的主应力；τxy、τyz、τzx分别为各个面上的切应力。

1.4 模型正确性验证

徐思博[21]采用高压气体推进装置对回转体进行入水试验，试验装置如图4所示，试验方形水箱尺寸为800 mm×800 mm×600 mm，其中水深500 mm，试验模型采用铝质平头实心回转体，直径8 mm、长40 mm，入水过程由高速摄像机拍摄记录，经计算回转体入水速度为78.5 m/s，记录入水空泡并计算入水后跨介质飞行器的速度衰减结果。

参照以上试验方案建立回转体垂直入水仿真速度验证模型，回转体采用弹塑性材料模型，水域大小为100 mm×100 mm×500 mm，空气域大小为100 mm×100 mm×100 mm。回转体高速入水初期空泡形态对比结果如图5所示。数值仿真计算的回转体高速入水速度衰减结果与文献[21]试验速度衰减结果对比如图6所示。结果表明，文献[21]中试验可观察到的自由液面的雾状喷溅主要由微小雾状液滴组成，仿真所得空泡最大直径、长度与试验结果误差均在4%左右。同时，速度仿真结果与试验结果绝对误差仅相差2.81%，下降趋势一致。误差可能原因是网格尺寸不够精细、未考虑初始静压且未考虑在高速入水时水的汽化。空泡在液面位置与试验存在误差，入水初期仿真的速度衰减更快，与试验结果误差最大为6.41%。但由于头部外形对跨介质飞行器在触水阶段的入水速度与所受最大冲击应力有重要影响。本文对跨介质飞行器头部外形在触水阶段进行仿真优化设计，可忽略初始静压与液体汽化的影响。在触水阶段，本文仿真数据与试验数据的误差在可接受范围内。

为进一步验证本文冲击应力仿真分析的准确性，参照Shi等[8]所做仿真数据完成回转体入水仿真冲击应力验证。跨介质飞行器入水速度为30 m/s，入水角度为60°，跨介质飞行器为铝质壳体，直径200 mm、厚度8 mm、长1 722 mm，水域大小为3 000 mm×500 mm×2 400 mm，空气域大小为3 000 mm×1 500 mm×2 400 mm。最大冲击应力数值仿真计算结果与文献[8]仿真数据进行对比，对比结果如表3所示。结果表明，本文冲击应力仿真结果与文献数据相差仅4.71%，吻合较好。

表3 冲击应力对比

2 优化过程

由于反复迭代有限元仿真需要消耗大量计算成本，同时头部外形优化可能存在多组局部最优解，采用基于RBF代理模型和MIGA的跨介质飞行器头部外形自适应优化流程，如图7所示。图7中，vmin、vmax分别为样本集中速度vend的最小值和最大值，Fmin、Fmax分别为样本集中冲击应力Fσ的最小值和最大值。具体流程包括：

步骤1进行计算试验设计，采用拉丁超立方在跨介质飞行器头部外形曲线参数X0、k1、k2设计区域内采样，并进行三维建模。

步骤2对模型进行前处理并进行有限元入水仿真，获取以在60 m/s垂直入水条件下，5 ms时跨介质飞行器的速度vend与入水过程中跨介质飞行器所受最大冲击应力Fσ。

步骤3以现有样本集的跨介质飞行器头部外形曲线控制参数X0、k1、k2为输入，对应速度vend与冲击应力Fσ归一化后的联合目标函数QM为输出,构建RBF代理模型。

步骤4归一化跨介质飞行器头部外形曲线系数X0、k1、k2为设计变量，速度vend与冲击应力Fσ归一化后的联合目标函数QM为优化目标，确定约束条件，采用MIGA对构建的代理模型寻优。

步骤5求解头部外形曲线控制参数优化结果X0op、k1op、k2op，以及最优结果对应的优化目标函数QMop，并计算与有限元模型对应归一化联合目标函数QM的绝对误差ERR，判断其是否小于收敛精度ε，是则停止，否则将寻优结果X0op、k1op、k2op及对应的联合目标函数QM作为新的样本点，加入样本集以构建新的RBF代理模型,返回步骤3继续迭代，直至满足收敛精度。

2.1 数据点采样

针对跨介质飞行器头部外形曲线的控制参数X0、k1、k2，在设计范围内采用拉丁超立方进行采样，能够从变量的设计区间进行均匀采样，在控制参数X0、k1、k23个变量的规定区间中取出18个样本，每个变量的累计分布被分成相同的18个小区间，从每一个区间随机的选择一个值，每一个变量的18个值和其他变量的值进行随机组合。随机采样数据较为集中，相比随机采样，拉丁超立方采样法通过小区间分别采样，保证采样的均匀性。如图8所示，调用模型计算，获取18组跨介质飞行器头部外形曲线系数X0、k1、k2，对应跨介质飞行器的速度vend与入水过程中所受最大冲击应力Fσ，构建联合目标函数QM，作为构建代理模型基础样本集。

2.2 径向基函数神经网络

以归一化跨介质飞行器头部外形曲线控制参数X0、k1、k2为输入，将入水后5 ms时跨介质飞行器的速度vend与入水过程中跨介质飞行器所受最大冲击应力Fσ归一化处理，构建联合目标函数QM为输出，构建RBF代理模型，其中QM如(15)式所示：

(15)

通过函数逼近理论模型，RBF代理模型是一种具有隐层的三层前馈型神经网络。RBF代理模型通过隐层将非线性问题输出为线性问题加权和，其中权重比可调，因此在处理系统内难以解析的规律性时，具有很好的通用性和很快的学习收敛速度，在不同领域内均得到广泛应用，能够逼近跨介质飞行器头部外形优化的非线性函数问题。以初始采样获取18组数据为基础样本集，训练获取初代RBF代理模型，其复相关系数R2远小于0.9，无法满足应用需求，选取添加近似最优解的自适应优化策略[22]。

在自适应优化策略中，全局近似精度不再是关注的重点，而是仅提高代理模型在兴趣区域(ROI)的局部近似精度，ROI为当前RBF代理模型存在全局最优解的区域，迭代过程中将上一迭代步的近似最优解添加到样本集中，将不断提高RBF代理模型在ROI附近的近似精度，可以降低优化过程的采样数量与次数，快速收敛，通过多岛遗传优化算法解出最优解。在本文的自适应近似优化策略中，只将每次迭代过程中的近似最优解添加至样本集中，不再进行单独采样，即使用MIGA对每次迭代新构建的RBF代理模型进行寻优，获得近似最优解(X0op，k1op，k2op)，(X0op，k1op，k2op)可视为在现有样本集下求解的最接近真实最优解的结果，若未达到收敛条件，将(X0op，k1op，k2op)作为新增样本点，通过仿真求得对应的速度vend、最大冲击应力Fσ与联合目标函数值QM构建新的样本集更新代理模型，并重复上述过程直至收敛，获取全局最优解。

2.3 多岛遗传优化算法

传统遗传算法具有较高的鲁棒性，通过模拟生物遗传过程可以快速求解非线性的复杂组合优化问题[23-24]。MIGA是一种基于遗传算法的改进算法，通过并行遗传算法提高计算效率，进行样本交叉提高全局优化能力。MIGA将子代种群均匀分布到多个岛上，每个岛上的种群并行运算传统遗传算法，一定的代数后每个岛上随机选择一些种群进行相互交换，再并行传统遗传算法运算，提高种群的多样性，避免传统遗传算法收敛到局部最优解的问题。跨介质飞行器头部外形曲线涉及多控制参数，这使得跨介质飞行器头部外形优化问题可能存在局部最优解，MIGA相比传统的遗传算法更适用于求解本问题。本文的MIGA设置中，将代数设置为100代，分为10个岛，每个岛取20组子代种群参数来寻找最优解。以归一化跨介质飞行器头部外形系数X0、k1、k2为设计变量，速度vend与最大冲击应力Fσ归一化后的联合目标函数QM为优化目标函数，确定约束条件，如(16)式所示：

(16)

s.t. 1≤X0≤4

0≤k1≤1

0≤k2≤15

2.4 优化结果

基于2.3节的思路，采用基于添加近似解的自适应优化策略，采取RBF代理模型及MIGA，经过6次迭代解得跨介质飞行器头部外形参数全局最优解为X0=4.000、k1=0.314、k2=14.223，最优解头部外形对应代理模型预测值与仿真结果对比如表4所示，归一化后联合目标函数QM的收敛精度为0.001，可以认为代理模型与优化算法对最优值的预测较为准确。

表4 代理模型预测值与仿真结果对比

最优解参数对应跨介质飞行器头部外形与3种基本头部外形参数如表5所示，头部外形对比如图9所示，最优解头部外形处于椭圆型与尖头部外形之间，前端为类椭圆，头部后端变化斜率相对于尖头部外形更平缓，可见跨介质飞行器的最优解头部外形结合了椭圆型与尖头部外形不同区域的特点。

表5 基本头部外形参数

3 最优结果分析

跨介质飞行器的最优解头部外形与3种基本头部外形的对比仿真结果如表6所示，对比可知，入水时跨介质飞行器所受冲击应力与头部外形关系较大。入水瞬间平头部外形跨介质飞行器与水面接触面积最大，入水过程中所受冲击应力峰值最大，表明相同条件下更容易发生局部结构的损坏失效。最优解头部外形跨介质飞行器比平头部外形跨介质飞行器所受最大冲击应力减少60.16%，可大幅提升跨介质飞行器入水时的结构安全性。

表6 最优解结果与基本模型结果对比

以入水速度方向为正方向，跨介质飞行器在触水阶段为减速过程，最优解头部外形与三种基本头部外形的跨介质飞行器的速度数值大小对比结果如图10所示，加速度数值大小对比结果如图11所示。结果表明，在跨介质飞行器头部接触水面瞬间，会受到极大的冲击加速度，由(13)式可知，所受冲击加速度与跨介质飞行器入水部分在速度方向上的横截面积和速度的平方呈正比，平头部外形跨介质飞行器由于触水面积最大，所受加速度数值最大，最优解头部外形跨介质飞行器相对平头部外形跨介质飞行器的速度衰减量大幅减小，最优解头部外形跨介质飞行器比平头部外形跨介质飞行器受到的最大冲击加速度绝对降低81.81%。虽然尖头部外形跨介质飞行器速度衰减量初期更小，但最优解头部外形跨介质飞行器后期速度衰减趋势更为平缓，约2 ms后最优解头部外形的跨介质飞行器的速度大于尖头部外形跨介质飞行器，由于尖头部外形跨介质飞行器入水瞬间接触面积最小，且触水面积变化更为平缓，所以尖头部外形跨介质飞行器整个过程中没有较大的加速度峰值，但由于2 ms内尖头部外形跨介质飞行器的速度更大，当跨介质飞行器触水部分横截面积相同时，尖头部外形跨介质飞行器所受阻力更大，因此加速度平均值高于最优解头部外形的跨介质飞行器。对比可知，最优解所对应的头部外形能够有效的减少入水过程跨介质飞行器的能量损失与所受冲击应力。