基于引力搜索算法的两级串联吻切锥乘波体设计

蔡致鹏,王新阁,张宴嘉

(空军航空大学， 长春 130000)

1 引言

高超声速武器是指飞行马赫数达到5马赫以上，能够在大气层或临近空间飞行的武器[1]。因其具有高空高速、不易拦截等优势而深受追捧。未来基于先进的超音速轰炸机平台，高超声速巡航导弹势必会得到更加广泛的应用。高超声速巡航导弹相对于亚音速巡航导弹有反应速度快、突防能力强、毁伤效果大等优点[1,2]，是目前各国都在研究的主流武器之一。美在高超声速武器研究中起步早，美国开发了多款高超声速武器。X-51A验证机进行了多次实验，验证机由B52轰炸机携带到1.5万m高空投放，之后火箭助推器把验证机加速到4.8Ma后脱落分离，机身内的超燃冲压发动机启动，速度保持在5.1Ma，飞行时长约为4 min。先进高超声速武器(advanced hypersonic weapons,AHW)计划的目标是发展一种射程6 000 km，飞行速度6马赫，精度在10 m以内的导弹，到2017年底前AWH进行过3次试飞，2次取得成功[3]。自美国退出《反弹道导弹条约》以来，俄罗斯为突破美国导弹防御系统，一直致力于研发高超声速武器。匕首高超音速导弹由陆基伊斯坎德尔-M近程弹道导弹改装而来，由米格-31K截击机或图-22M轰炸机携带投放，飞行速度10Ma，圆概率误差7 m。俄罗斯的“锆石”高超声速巡航导弹经过多次成功试射后已经服役[3-4]。

常规构型飞行器高超声速飞行时，下表面高压气体会上溢至上表面，上表面压力增大，导致升阻比降低。Kuchemann[5]根据风洞试验和试飞结果，提出了高超声速“升阻比屏障”，传统气动构型能达到的高超声速升阻比极限为4。乘波体构型是能突破“升阻比屏障”的飞行器外形之一，其设计思想是采用已知基准流场反设计飞行器外形。自1959年提出以来，发展出楔形乘波体、锥导乘波体、吻切锥乘波体和楔-锥乘波体4种外形[6]。面对日益复杂的周边环境和世界军事大国的领先地位，加紧对高超声速巡航导弹的研究具有重要意义。本文中基于低-高马赫数两级串联吻切锥乘波体对高超声速巡航导弹气动外形进行研究。

2 吻切锥乘波体模型建立

当超声速气流绕流无限长零迎角圆锥时，且气流马赫数和圆锥半顶角在一定范围内，圆锥前会产生锥面激波。乘波体能够将激波后的高压气体限制在乘波体的下表面，阻止下表面气体绕过前缘曲线与上表面形成连通，因而会得到更大的升力。

如图1所示，将激波出口型线划分为多个小段激波曲线，每小段曲线均可视为部分圆锥激波，产生该段激波的圆锥为吻切锥。预先选定的吻切锥乘波体流动捕捉曲线FCC(flow capture curve)通过自由流线法获得流动捕捉管FCT(flow capture ture)。FCT与每个吻切锥所产生的交线即为吻切锥乘波体的前缘曲线。前缘曲线与FCT合围部分为乘波体上表面。将前缘曲线按照激波出口型线划分成对应的几段，每段前缘曲线都进行离散化，分别进行流线追踪，最后将流线经过曲面拟合生成乘波体下表面。两级串联吻切锥乘波体由一级乘波体、衔接段、二级乘波体组成。一、二级长度h相同，采用桥接曲面连接，衔接段长度l为全长的十分之一。

图1 两级串联吻切锥乘波体截面示意图

当来流马赫数M，内圆锥半锥角δ，和流场长度L0固定时，乘波体的FCC曲线就成了改变乘波体外形的唯一变量。如图2所示，本文中设定吻切锥乘波体FCC曲线为：

(1)

式中：a、b、c均为实数，保证曲线的指数段和水平段二阶导数连续。式中的参数a、b、c、L1并不相互独立。改变两级串联吻切锥乘波体外形仅有c、L1两个变量。为满足宽速域要求，兼顾投放加速起点，一级乘波体设计马赫数为4，二级乘波体设计马赫数为8。

图2 乘波体底面截面图

3 基于引力搜索算法的两级串联吻切锥乘波体设计

3.1 多目标优化

多目标优化问题(multi-objective optimization，MPO)一般由N个决策变量参数，M个目标函数和J个约束条件组成，目标函数、约束条件与决策变量之间呈函数关系。对于最大化多目标优化问题的描述有：

Maxf(x)=[fi(x),i=1,2,…,M]

gj(x)≤0,j=1,2,…,J

(2)

3.2 万有引力搜索算法

引力搜索算法(gravitational search algorithm,GSA)是Rasshedi等人受牛顿的万有引力定律启发，于2009年提出的一种算法，引力算法的个体位置代表着问题的解，个体质量是对结果优劣的评价[9]。

设在D维空间中有N个粒子，则第i个粒子的位置为：

(3)

(4)

(5)

其中：fiti(t)为t时刻粒子i的适应度；worst(t)表示最差适应度；best(t)表示最佳适应度。在求解最大化问题中：

best(t)=maxfitj(t)

(6)

worst(t)=minfitj(t)

(7)

t时刻粒子i受到粒子j的引力为：

(8)

(9)

其中：G(t)为引力常数；为了防止分母为零，添加一个大于零的小常数ε；Rij(t)为粒子之间的欧氏距离。粒子i受到的引力为粒子群内除本身外其他粒子吸引力的合力。通过采用各粒子施加引力的随机加权和代表粒子所受的合力，从而加强随机性。

(10)

(11)

(12)

但在寻找最优解的过程中，结果可能很快收敛在次优解附近，并反复跳跃，难以得到全局最优解。为此，任意粒子i上所受合外力不再是种群中除粒子i外全部粒子所提供的引力和，算法将选择部分粒子发挥作用，该部分粒子的集合记为Kbest。调整后粒子i所受合力为：

(13)

Kbest内的元素为种群中除粒子i外适应度值排在前Kbest个的全部粒子，且集合内的元素会随着迭代步数的增加逐渐减少。从初始化的种群大小为N，经过多次迭代，Kbest内部的元素可能会减少到1个。为了算法能够更加准确高效，GSA算法要求引力常数G随时间t变化，调整后的万有引力常数G(t)为：

(14)

其中：G0为万有引力常数的初始值；β为引力常数缩减系数。整体算法流程如图3所示[10]。

图3 万有引力搜索算法流程框图

3.3 升阻比和容积率估算

高超声速飞行器的设计应同时满足高升阻比、高升力系数和高容积率需求[11-12]。本文中选择升阻比和容积率作为优化目标，容积率η的计算公式为[13]：

(15)

其中:V为乘波体飞行器的体积；Sw为乘波体飞行器的湿润面积。

由于种群数量N较大，为精简计算量，提高效率。采用正交实验和半经验公式推导总结设计参数与性能之间关系，并建立二元非线性回归模型，使用对升阻比和容积率的估算结果来代替模拟结果。通过改变c和L1这2个设计参数可以改变两级串联吻切锥乘波体外形，每个参数选取4个水平，故采用L8(42)正交表，共需要8次计算。在Fluent中进行8次计算得到结果汇总为表1。根据结果建立多元非线性回归模型。参考升阻比估算和容积率的计算公式[14]，可以得到升阻比L/D和容积率η都与L1和c成二阶关系。最后按照正交表的试验结果，最终建立函数模型：

0.037 6c-0.039 9L1

3.4 升阻比和容积率优化

容积率η和升阻比L/D是一对存在矛盾关系的参数。取维度d=2，粒子群数N=200，Pareto解集的最大规模为100，c的边界条件为[1,2.1 838]，L1边界条件为[2,4.867 6]，容积率η大于0.128，升阻比L/D大于3.2，设定迭代次数为200。得到容积率η与升阻比L/D的Pareto解集分布如图4所示。选取点A(0.140 001，4.740 26)和点B(0.128 073，5.514 019)建立优化模型。点A对应的c和L1为(2.1，3.7)，点B对应的c和L1为(1.8，3.3)。建立相应的模型A和模型B，如图5、图6所示。

表1 正交实验计算结果

图4 多目标优化的Pareto解集曲线

图5 模型A示意图

图6 模型B示意图

使用Fluent对模型进行升阻比和容积率计算，得到如表2所示的结果。经过计算，仿真计算结果与估算结果差值较小，证明估算结果具有可靠性。

表2 模型A、B的容积率、升阻比Table 2 Volume ratio and lift-drag ratio of model A and B

由此能够提出一种反设计思路，可以在设计导弹之前优先依据任务的需要，明确携带燃料和弹药部质量，由此选择合适容积率，然后通过图4中的解集，找到最大的升阻比。最后再依据选择的升阻比与容积率约束两级串联吻切锥乘波体的设计参数。

3.5 优化结果

选取模型B进行气动分析。根据图7(a)可以看出，在速度为6Ma、迎角为4°的条件下，气流经过模型B能够生成激波且提供较大的升力。在考虑底阻的情况下，尾部离体后的区域会形成一段低压区。可以通过在尾部安装动力装置的方式，利用尾部喷出的高压气体来抵消乘波体前后压力差，从而能够进一步增大升阻比。图7(b)中，可以在多个截面上观察到，在乘波体下表面会形成一个高压封闭区，虽然边缘位置有少量与上表面形成连通区，但在整体上的压力分布上保持一个较好水平。

图7 模型B的压力云图

如图8(a)、(b)、(c)和(d)分别为在来流速度为4Ma时模型B的一级乘波体底面、二级乘波体底面的压力云图和来流速度为8Ma时模型B的衔接段、二级乘波体底面的压力云图。

图8 乘波体截面压力云图

从4张图都可以看出，压力最大区域集中在下表面边缘位置。因而在边缘位置的封闭性就直接影响了升力的大小。对比图8(a)图8(b)可以发现，一级乘波体的下表面封闭性优于二级乘波体。一级乘波体的设计马赫数为4，二级乘波体的设计马赫数为8，故在4Ma条件下，一级乘波体的压力封闭性能更好。分析图8(c)与图8(d)，发现衔接段的封闭性要好于二级乘波体。将图8(b)与图8(d)进行对比，在8马赫的条件下，压力外泄的情况得到了明显改善。上述结果都与预期设定吻合，证明达到了反设计的预期效果。

图9为两级串联吻切锥乘波体在不同速度下的升阻比曲线。从图9中看出，从4～8Ma部分升阻比随着速度的增大会逐步增大，但是在马赫数大于8部分，升阻比趋于稳定且略有下降。根据马赫数无关原理[15-16]，当马赫数增大到一定程度之后，飞行器表面的升力系数与阻力系数将不随马赫数发生变化，所以升阻比在8Ma后逐渐趋于稳定。但是底部阻力会略有增大，从而导致升阻比略有下降。

图9 模型B升阻比与马赫数之间的关系曲线

4 结论

1) 采用正交实验研究了设计参数与升阻比容积率的关系，得到升阻比L/D和容积率η都与L1和c成二阶关系，并建立多元非线性回归模型。分析发现容积率高的乘波体，其升阻比性能不佳，而升阻比高的乘波体，容积率又无法满足要求，2个参数是矛盾关系。

2) 采用多目标万有引力搜索算法对两级串联吻切锥乘波体进行优化，在满足容积率需求的前提下得到最大升阻比的优化模型，并通过CFD方法验证，容积率和升阻比的误差在1%以内，证明了优化方法的可靠性。

3) 对优化后的乘波体模型仿真分析，一、二级乘波体在各自的设计马赫数下都表现出较好的封闭性，在6Ma时整体的压力分布较好，有较大升阻比，在马赫数超过8后，升阻比略有下降后趋于稳定。