滑翔制导炮弹鸭舵的气动外形快速优化研究

0 引 言

滑翔制导炮弹是通过在常规尾翼弹的基础上增加制导控制系统而形成的远程精确打击武器[1]。在弹道升弧段，鸭舵折插在弹体内，经过弹道顶点某位置后，鸭舵张开进行滑翔控制，从而保证炮弹在正攻角下始终以最大升阻比状态飞行。可见，鸭舵是滑翔制导炮弹在滑翔过程中极为重要的部件，其气动特性影响了全弹的升阻比、滑翔效率及稳定性等。因此，有必要对滑翔制导炮弹的鸭舵外形开展优化设计，以期提高滑翔制导炮弹的性能。

目前，国内外关于制导炮弹气动外形优化开展了一些工作。文献[2]利用DATCOM软件和粒子群算法来优化远程制导炮弹控制舵的外形参数，并通过Cart3D软件分析评估了其飞行性能；文献[3]提出了一种基于试验设计和响应面方法的制导炮弹外形设计优化框架，得到了一组射程最远的帕累托最优设计集；文献[4]研究了一种多学科优化设计方法，用于设计并评估无尾式远程制导炮弹的性能。国内方面，文献[5]对某滑翔增程制导炮弹进行了气动数值计算，并通过风洞实验进行验证；文献[6]建立了优化设计模型，并采用复形调优法对制导炮弹的鸭舵外形进行优化；文献[7]针对常规气动布局的制导炮弹提出了适用于概念设计阶段的气动外形优化方法，采用参数化建模方法和非线性面元法快速准确地预估炮弹气动特性。

为提高优化效率，国内外目前对滑翔制导炮弹的外形优化设计常采用工程算法计算气动参数，但获得气动参数的方法还包括数值模拟、风洞试验和靶道自由飞行试验等。在不具备开展风洞试验或自由飞行试验条件的初步设计阶段，工程算法和计算流体力学(Computational Fluid Dynamics, CFD)数值模拟是重要的气动参数预测手段。一般而言，CFD数值模拟的计算精度高于工程算法，如文献[8]所述，采用工程算法数据计算某型炮弹的轨迹时，射程误差可达10%，而采用CFD数据所得的误差为2%左右。但是，CFD数值计算耗时很长，难以实现快速优化。为此，研究者们开发了一些代理模型来代替耗时的CFD仿真。其本质是采用机器学习技术，从设计空间内的少量样本信息中学习目标函数随设计变量变化的规律。目前比较常用的代理模型有响应面模型、径向基函数模型、人工神经网络模型和Kriging模型等[9]。文献[10]采用人工神经网络模型实现了对导弹的气动参数预测；文献[11]则发展了一套基于Kriging模型的机翼气动优化设计方法。构造一个精确的气动力代理模型需要大量的样本点，若每个样本点都调用CFD数值模拟，其计算量同样不容小觑。考虑引入计算快速的工程算法数据辅助建立代理模型，将CFD数值模拟和工程算法结合起来，用较少的CFD数据实现高精度的模型，从而提高计算效率。

本文以滑翔制导炮弹的鸭舵为对象开展外形参数优化，充分考虑工程算法的快速性和CFD数值模拟的高精度性，构建代理模型近似代替对不同外形方案滑翔制导炮弹的气动力计算，并以此为基础建立了鸭舵外形优化的数学模型，通过遗传算法求解该模型得到最优外形方案，验证所提方法的可行性和有效性。

1 气动力快速计算方法

对于复杂的多设计变量气动外形优化问题，气动力计算结果的精准度影响着炮弹外形优化后的气动特性，若采用CFD仿真，那么气动力的反复计算在迭代优化的过程中就会耗时很长；若采用工程算法，计算精度则较低。如何在保证计算精度的前提下提高计算速度，就显得尤为重要。

本文采用Co-Kriging模型来代替在优化过程中对滑翔制导炮弹气动参数的计算，Co-Kriging模型是在Kriging模型的基础上发展起来的地质统计学方法[12]，通过该模型可以用少量的CFD仿真数据修正大量易获取的工程算法数据，从而建立高精度代理模型。相较于直接使用CFD和工程算法计算，该方法不但可以保证较高的精度而且计算时间大大降低，从而实现对鸭舵气动外形的快速优化。

1.1 Co-Kriging代理模型

Co-Kriging模型是利用少量高可信度数据来修正大量的低可信度数据，从而建立模型来近似代替炮弹外形参数与气动参数之间的映射关系y=f(x)。那么如何获得该映射关系最佳的近似函数就是代理模型方法的关键所在。

假定在样本集Xe处抽样获得精度较高的高可信度数据ne个，其函数值记为Ye，在样本集Xc处抽样获得精度较低的低可信度数据nc个，其函数值记为Yc。将样本集X处的函数值视为一个稳态随机过程，于是，对于Co-Kriging模型，可得到如下的随机变量场：

(1)

构造的代理模型的近似估计可以表示为

Ze(x)=ρZc(x)+Zd(x)

(2)

式中: Zc(x)为对低可信度样本数据构造的Kriging模型；Zd(x)为对高、低可信度样本数据差值构造的Kriging模型。这里实质上是将高可信度数据近似看作低可信度数据乘以一个比例因子常量ρ后，再加上一个表示ρZc(x)和Ze(x)之间差别的高斯随机过程Zd(x)。考虑两组独立的高、低可信度数据，可以构建Co-Kriging完整的协方差矩阵：

(3)

式中:为Zc(x)和Zd(x)的方差；ψc(Xe, Xc)为Xe和Xc之间的相关性矩阵，这些相关性依赖于样本点间的绝对距离以及参数θ和p。与Kriging模型不同的是，这里有两种不同的空间相关函数ψc和ψd，就需要训练更多的模型参数：θc，θd，pc，pd以及比例因子ρ。为了估计参数θd, pd和ρ的值，首先定义：

d=Ye-ρYc(Xe)

(4)

其中，Yc(Xe)是在点Xe处低可信度模型的值Yc。这样就可以采用与标准Kriging类似的优化过程来求得对θd，pd及ρ的估计[13]。Co-Kriging模型预估值表达式如下：

(5)

考察一下采用该预估值表达式来预测某一高可信度样本点处的函数值。假设有 且c是C的第nc+i列，那么cTC-1是一个第nc+i个元素为1、其余元素为0的单位向量，则 可以认为，该式是对高可信度数据的一个插值模型。从另一个角度，一定程度上也可认为是对低可信度数据的回归，并使得在高可信度样本点处的预测值与已知的函数值Ye一致。

1.2 数值方法与验证

高可信度样本点的响应值由CFD数值计算得到。为保证计算精度，采用结构化网格划分流场并在弹箭头部和翼面等气流变化剧烈处加密网格，以提高网格质量和数量。求解器中，以N-S方程作为流场计算的主控方程，采用realisable K-epsilon湍流模型、二阶迎风差分格式以及密度基隐格式，对滑翔制导炮弹的绕流场进行数值模拟。

为了验证拟采用的计算方法，选取120 mm迫击炮弹为对象进行气动力计算验证，如图1所示。

图1 气动参数随马赫数变化曲线

Fig.1 Variation curve of aerodynamic parameters with Mach number

所得计算结果与文献[14]中提供的气动数据规律基本一致，且升力系数导数的误差小于5%、俯仰力矩系数导数的误差小于9%。这表明拟采用的CFD数值计算方法可行有效。

1.3 代理模型构建过程

作为代替真实气动力计算方法的工具，代理模型的建立必然离不开先验数据的获取与计算，其构建过程如图2所示。

图2 代理模型构建过程

Fig.2 Establishing process of surrogate model

(1)确定鸭舵各个外形参数的取值范围，形成一个多维设计空间。

(2)在设计空间中采用全因子试验选取大量的低可信度样本点构成样本集Xc，并均匀选取少量的高可信度样本点构成样本集Xe。其中，高可信度样本点是低可信度样本点的子集(Xe⊂Xc)。

(3)对低可信度样本点采用文献[6]中的工程算法计算其气动参数。高可信度样本点的气动参数则采用CFD数值模拟来求解。

(4)基于以上高、低可信度数据组成的样本库，采用第1.1节方法构建代理模型。

所建立的Co-Kriging代理模型，其输入为鸭舵的外形设计变量xi，输出为炮弹的气动参数。利用该模型，设计空间内的其他外形方案，就无需经过复杂的CFD仿真而能快速预估其气动参数并保证足够的精度。相较于传统的气动力计算方法，这样一个有两种数据来源的多可信度代理模型不仅有着更快的计算速度，而且可以过滤数值噪声并显著提高优化效率。

2 鸭舵优化设计模型

滑翔制导炮弹通常是在尾翼弹的基础上通过加装制导装置和控制机构而来，其气动外形设计可在弹体确定的前提下，在弹头部添加具有适当外形的舵面。由于鸭舵对制导炮弹滑翔性能具有重要影响，本文选取鸭舵作为优化对象，对其主要外形参数进行优化设计，从而使炮弹具有较高的滑翔控制能力。

2.1 设计变量及目标函数的选取

依据舵面参数确定的原则[15]，为了兼顾舵面有较好的气动特性、较小的舵面铰链力矩和舵面结构的强度和刚度等，选用后掠的梯形舵面。如图3所示，其外形参数主要包括舵面的展长Bc(本文以展长的一半作为设计变量，即 根弦长Cc、舵面前缘至弹顶的距离Xc和前缘后掠角χc，以上4个主要几何参数就是鸭舵外形优化的设计变量，即

图3 某滑翔制导炮弹鸭舵的外形设计变量

Fig.3 Shape design variables of canard for the gliding guided projectile

X=[x1, x2, x3, x4]T

(6)

对于滑翔制导炮弹，为提高射程，通常保持在最大升阻比状态下滑翔飞行。因此，本文选取制导炮弹的升阻比K作为优化目标函数。当4个设计变量确定后，就可以在一定的马赫数和雷诺数下，计算出制导炮弹的升阻比，其函数关系为

(7)

式中: K为全弹升阻比; Ry, Rx为飞行弹道上全弹的升力和阻力。

2.2 约束条件

约束条件是对优化设计问题提出的一些条件限制，以使设计方案在满足这些限制下达到较优。对于滑翔增程制导炮弹，所参加调优的设计方案应保证如下条件。

(1)稳定性约束条件

鸭舵张开不偏转时，全弹静稳定性储备量ε应满足：c1≤ε≤c2，其中c1，c2为预定常数。静稳定储备量可以表示为

(8)

式中: XCP为全弹压心至弹顶距离；XCG为全弹重心至弹顶距离；l为弹长。

(2)机动性约束条件

制导炮弹在飞行过程中，假设在任意时刻都处于平衡状态，即舵面偏转时作用的力矩在每个瞬时都处于平衡状态。由此可得

(9)

即

(10)

式中:为弹体的特征面积；为弹身和鸭舵的升力系数导数；为单位舵偏角所能引起的平衡攻角，该值越大，则稳定性越小，机动性增加。鸭舵张开进行控制时，应满足机动性要求： 其中c3，c4为预定常数。对于制导炮弹在舵面实施控制时，需要兼顾炮弹的稳定性和机动性。

(3)尺寸约束条件

制导炮弹在发射前，通常将尾翼和舵面折叠在弹体内，因此，希望舵面的尺寸、体积满足火炮发射条件，即

(11)

式中：为1对鸭舵展长的最大允许值；为1个舵面的面积最大值。

(4)几何约束条件

受鸭式气动布局限制及要求，鸭舵位置受一定几何关系约束，即

(12)

式中：为舵面前缘至弹顶的允许最近距离；XCG为全弹重心至弹顶距离。

2.3 优化流程

综上所述，滑翔制导炮弹鸭舵的外形优化数学模型可以表示为

(13)

式中：E4为四维欧氏空间；R为满足约束条件的实数空间；gi(X)为上述约束条件。

综合考虑，采用遗传算法进行目标函数寻优，该算法具有通用性强、稳定性好等特性，且有着良好的全局搜索能力。滑翔制导炮弹鸭舵外形的优化过程大致如下：

(1)确定鸭舵外形优化的设计变量、约束条件及目标函数等；

(2)利用遗传算法在鸭舵外形参数构成的设计空间中生成初始种群；

(3)通过训练好的Co-Kriging代理模型快速计算种群中对应不同鸭舵外形的升阻比K；

(4)在约束条件下对目标函数，即升阻比K进行适应度评估，若不满足收敛条件则使种群交叉变异生成新种群并返回步骤(3)，直至收敛输出最优结果；

(5)采用CFD数值模拟对遗传算法的收敛结果进行验证，若误差较大则重新构建代理模型，若误差较小则确定为最优外形方案。

3 算例仿真

3.1 算例分析

本文以某大口径滑翔制导炮弹为例，该炮弹尾翼3对、舵面2对，舵面与尾翼以“++”形式配置。该炮弹鸭舵的各个设计变量及其取值范围如表1所示。

表1 鸭舵外形参数样本取值范围

Table 1 Sample range of canard shape parameters

在鸭舵外形参数构成的4维设计空间中(表1)，采用4×5×3×3的全因子试验设计均匀采样。全因子设计可以均匀填充样本空间，当因子数不超过5时较为合适。共采集低可信度样本点180个，并在设计状态：来流马赫数Ma=0.7、攻角α=6°的飞行条件下选用工程算法计算。

为方便分析，利用控制变量法筛选部分样本点，得到关于各个外形参数对全弹升阻比K及静稳定储备量ε的影响趋势曲线，如图4所示。由图可知，全弹升阻比随着鸭舵展长的增大而不断增大，随着鸭舵前缘后掠角的增大而减小，并且随着鸭舵根弦长的增大而缓慢升高。鸭舵至弹顶的距离越短，全弹的压心不断前移，则炮弹的静稳定度就越小。通过对工程算法的计算结果分析，可以认为在4维设计空间内存在一个点，可以使得升阻比达到一个最大值。

图4 鸭舵外形参数对全弹升阻比及稳定性的影响

Fig.4 Influence of canard shape parameters on lift-to-drag ratio and stability

如表2所示，在设计空间中以样本点2为中心点均匀采集高可信度样本点10个。对表2中每一个高可信度样本点所对应的滑翔制导炮弹分别建立模型、划分流场网格并采用CFD计算气动参数，计算条件与上述低可信度样本点一致。以某外形为例，弹体表面网格如图5所示，结构化网格数量900万左右，网格质量0.7以上，以此保证足够的精度。

表2 高可信度样本点采样分布

Table 2 Sampling distribution of high reliability sample points

图5 弹体表面网格

Fig.5 Surface grids of projectile body

利用上文10个高可信度数据和180个低可信度数据构建Co-Kriging代理模型，实现在优化设计过程中对不同鸭舵外形的滑翔制导炮弹的气动力快速计算。该模型输入为鸭舵的外形参数，输出为气动特性参数(即全弹升阻比)。采用遗传算法作为优化方法在约束条件下搜索最大升阻比，经反复试算，设置种群规模为20、变异概率为0.2、超过进化代数30代时优化停止，可取得较好优化效果。

3.2 仿真结果

对构建的代理模型进行验证。取9个点(与构建代理模型的样本点不重复)分别计算其高精度CFD仿真结果和Co-Kriging代理模型预测值，结果如图6所示。观察图6可以发现，代理模型的预测结果与CFD仿真计算结果差距不大，样本点5处的误差最大，为3.4%，误差最小处(样本点2)只有0.08%，平均误差为1.94%。同时，代理模型对随外形参数变化的样本点的预测结果趋势与CFD计算趋势相同。这表明，本文构建的多可信度代理模型有较高的精度，可以代替CFD数值模拟用于对制导炮弹的气动参数预测。

图6 预测值和计算值对比

Fig.6 Comparison between the predicted values and the calculated values

在多可信度代理模型基础上，按照制导炮弹的鸭舵气动外形优化模型编制计算程序。遗传算法对目标函数的优化收敛历程如图7所示，当种群更新至第七代时开始逐渐收敛。

图7 遗传算法收敛曲线

Fig.7 Genetic algorithm convergence curve

表3为优化前后具体结果，优化后鸭舵外形如图8所示，实体部分为鸭舵外形优化后的炮弹头部，虚影部分为优化前鸭舵的初始外形方案。与初始方案相比，优化方案的半展长增加了24%，根弦长增加了10.76%，至弹顶距离缩小至199.85 mm，后掠角变为12.78°。

表3 优化结果

Table 3 Optimization results

图8 头部鸭舵外形优化前后对比

Fig.8 Comparison before and after optimization of head canard shape

据优化结果，其升阻比为2.34，采用CFD数值仿真验证其阻力系数0.41、升力系数0.95，升阻比为2.317，升阻比误差在0.99%，表明优化结果有较高可信度。优化方案相较于初始方案，升阻比提升了11.4%。该外形下的制导炮弹，单位舵偏角下产生的平衡攻角为α/δ=0.73，依据文献[15-16]判断其满足操纵性约束条件。鸭舵张开后其稳定储备量为5.76%，满足在2%～6%之间的稳定性设计要求[17]。

对较优外形采用CFD数值仿真计算，图9为全弹升阻比在不同攻角下随马赫数的变化曲线。由图可知，全弹升阻比在亚音速范围内较大，在跨声速时期随马赫数增加迅速减小，这是因为在跨声速时由于激波的出现阻力陡增。全弹的最大升阻比可以达到2.5左右，满足鸭式布局制导炮弹在亚声速段机动飞行的总体要求。

图9 不同攻角下升阻比变化曲线

Fig.9 Lift-to-drag ratio curves at different angles of attack

经以上计算获得制导炮弹最优外形下气动参数后，采用有控质点弹道模型[18]对该滑翔制导炮弹进行弹道仿真：假设炮弹初始质量为45 kg，火箭助推发动机平均推力为1 200 N，发动机装药质量耗散速率为0.535 7 kg/s，采用标准气象条件，初速930 m/s，射角61°，出炮口后，助推火箭点火工作14 s，到达弹道最高点附近时张开鸭舵保持全弹在6°正攻角下滑翔飞行。飞行弹道轨迹如图10，无控弹射程仅42 km，在鸭舵控制下初始方案射程为64 km，优化方案的射程可达到72 km。相比于初始方案，鸭舵外形优化后的射程提升了12.5%，在仅仅优化鸭舵外形参数的条件下，这一增程效果是较为可观的。

图10 滑翔制导炮弹弹道仿真

Fig.10 Trajectory simulation of gliding guided projectile

值得注意的是，对于上述算例，对10个高可信度样本点的仿真评估过程共花费30 h，对180个低可信度样本点的评估、训练代理模型以及程序优化时间一共只需25 s左右(采用Intel Core四核处理器、主频2.6 GHz、8 G内存)。若所有样本点均采用CFD数值模拟来评估气动特性，则计算量是巨大的，需要耗时数天乃至数月[8]。后续利用该快速优化方法，可以考虑更多的设计变量、约束条件及目标函数等，以满足复杂多样的作战需求。

4 结 论

为了兼顾气动参数的计算精度和效率，本文设计了一种基于多可信度代理模型的气动力高精度快速计算方法，能够快速预估不同外形参数对应的炮弹气动力系数。与基于CFD的气动外形优化相比，本文所建立的优化过程大幅缩短了优化时间，升阻比的计算误差约为1.94%。以某大口径滑翔制导炮弹为算例，对其鸭舵外形的主要几何参数进行优化，在满足机动性和稳定性的条件下，炮弹的滑翔升阻比提升了11.4%，射程增加了12.5%，验证了本文方法的可行性和有效性。