高超声速滑翔飞行器再入气动系数改进拟合模型

2021-11-14蔡光斌张胜修

宇航学报 2021年9期

徐慧，蔡光斌，张胜修

(火箭军工程大学导弹工程学院，西安 710025)

0 引言

高超声速滑翔飞行器凭借着全球快速到达、高空高速突防、大范围机动等优良特性，成为世界各军事强国争相研究的热门领域[1-4]。近年来，高超声速滑翔飞行器再入轨迹优化问题成为飞行控制领域内研究热点[5-7]。基于空气热力学外形的通用飞行器[8-9](Common aero vehicle，CAV)成为再入轨迹优化问题研究中常用的研究对象。基于此外形，美国波音公司和洛马公司分别设计了CAV-L和CAV-H飞行器，最大升阻比分别约为2.5和3.5。

再入机动飞行过程中，高超声速滑翔飞行器主要利用气动力在大气层进行远距离滑翔，并且受到非线性气动力和过载、动压、热流率等过程约束影响，使得再入轨迹优化成为一个极具挑战性的问题[10-14]。由于再入过程复杂，气动升力、阻力会影响到再入过程中飞行器速度、高度、航向角、航迹角的变化，进一步影响飞行器再入过程中的热流率、动压和过载变化。因此，建立精确的高超声速滑翔飞行器再入气动系数模型，对研究高超声速滑翔飞行器再入轨迹优化，具有较大的理论意义与应用价值[15]。

在较早的高超声速滑翔飞行器再入轨迹优化问题研究中，一般采用的是只考虑攻角影响得到的较为简化的气动系数模型。实际上，复杂的再入环境导致高超声速滑翔飞行器气动系数受到攻角、马赫数、大气密度等多种因素的影响。只考虑攻角这一单变量影响的气动系数模型会带来较大拟合误差。李惠峰[16]等对CAV气动外形进行了参数化建模，并分析了得到模型的气动特性，但得到的模型纵向静不稳定，俯仰通道气动控制效率较低。Gao等[17]考虑攻角和马赫数的共同影响，采用二元二次多项式的形式，进行高超声速滑翔飞行器气动系数的计算。孙勇[18]等通过对升力系数和阻力系数的分析，在拟合模型中使用马赫数的负指数幂形式，提高了模型精度，减小了拟合误差，但该模型考虑的攻角和马赫数影响仍不够充分，模型拟合精度有进一步提高的空间。

通过对上述多种拟合模型的对比分析，本文提出了一种改进的高超声速滑翔飞行器再入气动系数拟合模型。本文采用升阻比较大的CAV-H飞行器作为气动系数模型的研究对象，利用多元非线性最小二乘的方法对本文提出的改进模型进行了参数辨识。引入一种拟合优度指标，将本文模型与文献[17-18]中的模型进行对比研究。从升阻比的角度分析了本文气动系数模型对高超声速滑翔飞行器再入轨迹优化的影响。拟合仿真结果表明，本文提出的改进模型拟合精度较高。通过进一步的仿真实验讨论了改进气动系数模型的升阻比特性，并对三类模型进行再入轨迹优化的对比实验，仿真结果表明，基于本文提出的改进气动系数模型，得到的优化轨迹更为平稳。

1 CAV-H气动系数及对比模型分析

1.1 CAV-H气动数据

准确有效的真实气动数据和高精度的气动拟合模型，对于研究高超声速滑翔飞行器轨迹优化与制导具有十分重要的意义。由于世界范围内高超声速飞行器风洞实验数据保密性高，目前仅能查阅到极少部分的公开数据。因此，本文采用文献[17-19]中给出的CAV-H的升力系数和阻力系数数据开展此类研究。

CAV-H升力系数、阻力系数和升阻比分别在表1、表2和表3中给出[19-21]。

表1 CAV-H升力系数Table 1 CAV-H lift coefficient

表2 CAV-H阻力系数Table 2 CAV-H drag coefficient

表3 CAV-H升阻比Table 3 CAV-H Lift-to-drag ratio

图1和2是分别是升力系数和阻力系数随攻角变化的情况，可以看出升力系数几乎与攻角成线性关系，而阻力系数与攻角的关系接近二次函数，但由于攻角为自变量的数据较少，并不能简单地确认升力系数、阻力系数与攻角的关系。

图1 升力系数随攻角变化Fig.1 Lift coefficient varies with the angle of attack

图2 阻力系数随攻角变化Fig.2 Drag coefficient varies with angle of attack

图3和4分别是升力系数和阻力系数随马赫数的变化情况。显然，升力系数和阻力系数随马赫数的增加而逐渐减小，在马赫数大于Ma10后，两者的下降趋势较为平缓。经分析可知，可以采用二次项或者负指数幂的形式描述升力系数及阻力系数与马赫数之间的关系。

图3 升力系数随马赫数变化Fig.3 Lift coefficient varies with Mach number

图4 阻力系数随马赫数变化Fig.4 Drag coefficient varies with Mach number

1.2 对比典型模型

常见的再入轨迹优化使用传统的气动系数模型，仅考虑攻角带来的影响，忽视了马赫数的影响，带来的拟合精度问题已经在文献[18]中进行了一定讨论。这里我们不再对单变量模型做深入研究，本文着重探讨综合考虑攻角和马赫数影响的拟合模型。文献[17]采用二元二次多项式形式的模型来拟合CAV-H气动系数，表现形式如下：

(1)

式中:CL,CD分别为升力系数、阻力系数；α为攻角,M为马赫数;ai,bi分别为升力系数、阻力系数拟合模型的待辨识参数，其中i=0,1,2,3,4。文献[18]提出使用马赫数负指数幂形式，来进行气动系数拟合，并且通过讨论分析表明，在阻力系数拟合模型中可以忽略攻角一次项。文献[18]中的气动系数拟合模型如下：

(2)

式中:CLi,CDi分别为升力系数、阻力系数拟合模型的待辨识参数，其中i=0,1,2,3。

高超声速滑翔飞行器一般具有高升阻比的气动外形，且升阻比越大，机动性能越好。在研究高超声速滑翔飞行器再入轨迹优化问题时，升阻比η计算如下：

(3)

上述式(1)和式(2)两种模型的系数如下：

[a0,a1,a2,a3,a4]T=[0.11139, -0.019871;4.161×10-4, 2.2991, 1.2292]T; [b0,b1,b2,b3,b4]T=[0.23462, -0.02421, 7.089×10-4, -0.17481, 2.7251]T; [CL0,CL1,CL2,CL3]T=[-0.2335, 2.9451, 0.2949, -3.3943×10-4]T; [CD0,CD1,CD2,CD3]T=[0.0234, 2.3795, 0.3983, -1.0794×10-4]T

表4～6中的数据，展现了两种模型的拟合结果。对表4～5分析来看，文献[18]中提出的模型升力系数和阻力系数拟合方面误差较大，但是对表6拟合升阻比结果分析时，文献[17]中模型的拟合误差远远超过文献[18]中的模型，可以看出，使用马赫数负指数幂形式进行拟合具有一定优势。

2 改进气动系数拟合模型及分析

2.1 改进气动系数拟合模型

通过对CAV公开数据的分析，以及对文献[17-18]中拟合模型结果的对比，可以发现文献[17]中通过使用攻角和马赫数的二次项进行拟合，得到的升力系数、阻力系数的拟合精度较高，但升阻比的拟合精度远小于文献[18]中采用马赫数负指数幂形式的模型。因此，综合考虑攻角和马赫数对升力系数、阻力系数的影响，本文提出了一种新的改进拟合模型如下：

(4)

式中:li,di分别为升力系数、阻力系数拟合模型的待辨识参数，其中i=0,1,2,3,4,5。

本文模型得到升阻比η计算公式为：

(5)

2.2 多元非线性最小二乘参数辨识

本文采用多元非线性最小二乘的方法辨识模型参数。非线性最小二乘法是以各数据点误差平方和最小为目标函数，来估计非线性模型中的参数。

本文需要辨识参数的模型，自变量为攻角α和马赫数M，因变量为气动系数C。则可以将气动系数C表示如下：

C(α,M)≈f(α,M,p)=p0+p1α+p2α2+p3M+p4ep5M

(6)

式中:p=(p0,p1,p2,p3,p4,p5),f(·)为拟合模型。现已知数据(αi,Mi,C(αi,Mi))，可以将多元非线性最小二乘模型参数辨识问题转化为寻找参数p，使得下式最小的优化问题：

(7)

式中:p为待优化参数，F(·)为拟合模型误差。

根据表1～3中的CAV-H的升力系数阻力系数数据，采用多元非线性最小二乘的方法对本文提出的模型进行参数辨识。得到CAV-H的新的气动系数模型的参数。

CAV-H气动系数模型参数辨识结果如下：

(8)

2.3 拟合模型的评价指标

常用的模型拟合评价指标，包括和方差(Sum square error，SSE)、均方根误差、均方差和平均绝对误差。但这四种常见的性能指标，仅体现模型输出与真值之间的误差大小，不能衡量拟合模型输出与真实数据变化的趋势是否一致。拟合优度能够较好地刻画模型输出与真实数据的动态变化趋势[22]。本文在和方差拟合性能指标之外，引入拟合优度指标，进一步评价模型拟合好坏。拟合优度记作τ,τ值越接近1，表明提出的模型拟合效果越好，越能刻画真实的数据；反之，越接近0，表明提出的模型效果拟合越差，用该模型拟合出来的气动系数结果越不准确。拟合优度的计算如下：

(9)

注1.拟合优度是回归问题中常用到的表述模型对数据解释程度的指标，相对于仅使用SSE类型的指标，能够更准确评价模型的好坏。本文将其引入到CAV-H气动系数拟合模型的评价问题中，后续的仿真实验证明本文改进的气动系数拟合模型拟合优度更大，即能更准确的拟合和解释CAV-H气动系数数据。

2.4 模型升阻比特性分析

高升阻比飞行是高超声速滑翔飞行器的重要特点。对高超声速滑翔飞行器机动飞行的升阻比特性的充分讨论，有利于后续轨迹优化问题研究。因此，针对本文提出的改进气动系数拟合模型的升阻比特性进行分析。

1)改进气动模型的最大升阻比分析

首先讨论最大升阻比与攻角、马赫数的关系，可将式(5)改写为：

(10)

式中:lM=l0+l3M+l4el5M,dM=d0+d3M+d4ed5M。

令式(10)对攻角α求导且等于0，可得:

(11)

将α*代入式(5)，即可得最大升阻比与马赫数的关系为：

(12)

通过式(12)，可以得到本文气动系数模型的最大升阻比对应的攻角随速度变化的关系，在后续的仿真实验中会展示相应的仿真结果。

注2.高升阻比气动特性是高超声速再入滑翔飞行器一个重要特点，对本文改进气动系数模型最大升阻比的讨论，有利于掌握该模型最大升阻比条件下攻角与速度之间的关系，对后续应用该模型的再入轨迹优化与制导问题研究具有重要意义。

2)不同升阻比对再入飞行轨迹的影响分析

在2.4节的基础上，讨论本文改进气动模型不同升阻比对高超声速飞行器再入飞行的影响。

对于不同升阻比，可以得到相应的攻角-速度剖面。定升阻比情况下，式(10)可改写为：

(13)

注3.高升阻比的气动特性直接关系到高超声速再入飞行器的机动能力，对改进模型的不同升阻比再入飞行讨论，得到不同升阻比对应的攻角-速度剖面解析式形式，为了保持一定的机动能力，可以由式(13)得到攻角-速度剖面内的飞行走廊，能够对考虑禁飞区约束的高超声速滑翔飞行器再入轨迹优化与机动制导问题研究中攻角-速度剖面的设计，带来一定的指导意义。

2.5 高超声速滑翔飞行器再入轨迹优化模型

本文采用的高超声速再入轨迹优化模型与文献[19-20]中模型一致。高超声速滑翔飞行器再入轨迹优化问题的动力学方程为：

(14)

式中:r是飞行器的地心距;V是飞行器的地球相对速度;ψ与γ分别是飞行器的航向角与航迹角;飞行器所处的经度θ纬度φ是判断飞行器路径约束的主要因素;m和g是飞行器的质量和当前地心距的重力加速度;飞行器的侧滑角σ与攻角α分别控制飞行器制导策略中的横向、纵向制导剖面。D=ρV2SrefCD/2与L=ρV2SrefCL/2是飞行器在飞行过程中的气动阻力与升力，其中ρ是飞行器当前高度的空气密度，S是飞行器的参考横截面积，CL与CD分别是与飞行器攻角有关的空气动力学参数。该模型中，侧滑角σ与攻角α是控制量。

在高超声速飞行器再入轨迹优化问题中，飞行器飞行过程还需满足路径约束、终端约束、控制量约束和状态约束。

1) 路径约束

在飞行过程中需要满足热流率约束、过载约束、动压约束等常规路径约束:

(15)

2) 控制量约束

(16)

式中:σmin、σmax分别为倾侧角的最小值和最大值；αmin、αmax分别为攻角的最小值和最大值。

3) 状态约束

(17)

式中:下标“min”和“max”分别表示各状态量最小取值和最大取值。

4) 终端约束

(18)

式中:tf指终端时刻；下标“f”指期望的终端状态。

5) 性能指标

本文进行讨论的仿真实验中，3.3.2节中的选用最大航程作为性能指标，对本文模型的三种不同升阻比下飞行性能进行考察，即本节仿真实验的性能指标为：

J1=min(-θ)

(19)

3.3节中的选用最小时间作为性能指标，对三种不同模型，在完成任务情况下对飞行性能进行考察，即本节仿真实验的性能指标为：

J2=mintf

(20)

3 仿真实验与结果分析

采用多元非线性最小二乘法对本文提出的模型进行参数辨识，得到辨识结果后，将本文模型同文献[17-18]中模型进行对比。仿真实验包括：1)本文改进模型对CAV-H气动系数数据的拟合，并与文献[17-18]中的模型拟合结果进行对比；2)从升阻比的角度对本文气动系数模型的飞行特性进行讨论；3)应用本文改进气动系数模型的再入轨迹优化仿真实验，并与其它两种模型进行对比。

本文仿真实验环境CPU为Intel i7 9700处理器，运行速度3GHz。

3.1 模型对CAV-H气动系数数据的拟合

利用多元非线性最小二乘方法对本文所提模型进行参数辨识。并将本文模型与文献[17-18]模型的拟合结果进行对比，如图5(a)、(b)、(c)分别是升力系数、阻力系数和升阻比的拟合结果。

图5 三种模型对CAV-H数据拟合Fig.5 Three models fit CAV-H data

表4～6中数据分别展示了三种模型的升力系数、阻力系数与升阻比相对于真值的误差大小，表中粗体的数据为误差较小的数据。比较仿真实验结果可得，本文模型升力系数拟合平均误差相对于文献[18]降低10%，相对于文献[17]降低5%；本文模型阻力系数拟合误差相对于文献[18]降低11.4%，相对于文献[17]同样降低11.4%；特别地，本文模型升阻比拟合误差相对于文献[18]降低了26.6%，相对于文献[17]降低52.6%，拟合精度明显提升。可以看出本文提出的新的CAV-H气动系数拟合模型，能够准确地得到气动系数，明显地降低拟合误差。

表4 升力系数拟合误差Table 4 Lift coefficient fitting error

表5 阻力系数拟合误差Table 5 Resistance coefficient fitting error

表6 升阻比拟合误差Table 6 Lift-to-drag ratio fitting error

表7和图6展示了三种模型拟合优度的对比，表中粗体的数据为拟合优度较大的数据。可以看出本文提出的模型，与文献[17-18]中的模型相比，对CAV-H升力系数、阻力系数和升阻比的数据拟合优度更接近1。特别地，本文模型在升阻比拟合方面，拟合优度相对于文献[17]提高了52.4%。由此可知，本文模型对样本数据随自变量的变化趋势刻画更为精确，对样本数据拟合更为准确。

表7 拟合优度对比Table 7 Goodness of fit comparison

综上所述，本文提出的改进气动系数模型，能够更好地描述CAV-H高超声速飞行器的气动系数变化，得到的气动数据比传统模型结果更为精确。

3.2 本文改进模型的升阻比特性分析

3.1节的仿真结果表明，本文提出的模型对数据的解释性更好，拟合精度更高。利用式(10)～(12)对此模型的升阻比特性进行分析，有利于此模型更好的应用于高超声速飞行器再入轨迹优化问题研究。

3.2.1改进气动系数模型最大升阻比分析

对于本文提出新的气动系数拟合模型，进行升阻比特性的分析。图6和图7是利用式(5)进行的升阻比特性分析得到的结果。图6是升阻比随攻角变化的曲线，由图可知，当马赫数一定时，随着攻角的增大，飞行器升阻比先增后减。由图6可知当马赫数约为Ma14时，升阻比能取得较大的值。图7是升阻比随马赫数变化的曲线，由图可知，当攻角一定时，随着马赫数的增大，飞行器升阻比先增后减。由图7可知当攻角约10°时，升阻比能一直保持较大的值。图8～9是由式(11)、(12)，进行仿真实验得到的结果。分析图8和图9曲线可知，当马赫数为Ma13.6，攻角为9.583°时，本文提出的气动系数拟合模型能够达到最大升阻比约为3.679。这与图7和图8的结果分析相一致。需要指出的是，CAV公开的资料中最大升阻比约为3.5，本文提出的模型最大升阻比基本与此数据一致。

3.2.2定升阻比对再入飞行轨迹的影响分析

为分析本文模型在不同升阻比情况下对飞行轨迹产生的影响，选取下面三种定升阻比情况进行仿真：1)η=3.2，升阻比稍大；2)η=3.0，升阻比中等大小；3)η=2.7，升阻比稍小。

本文提出的气动系数模型理论最大升阻比为3.679，但是从图6～9中可以发现，当定升阻比大于3.5时，对应的马赫数、攻角的变化范围非常小，不利于对本文气动系数模型的飞行特点进行分析，因此本文并未采用最大升阻比进行仿真实验。

图6 本文模型升阻比与攻角的关系Fig.6 Relationship between lift-drag ratio and angle of attack

图7 本文模型升阻比与马赫数的关系Fig.7 Relationship between lift-drag ratio andMach number

图8 最大升阻比与马赫数的关系Fig.8 Relationship between maximum lift-to-drag ratio and Mach number

图9 最大升阻比对应攻角与马赫数的关系Fig.9 Relationship between the angle of attack of maximum lift-to-drag ratio and Mach number

图10～图15是在本文模型基础上进行不同升阻比的再入飞行仿真结果。可以看出，采用本文气动系数模型进行定升阻比飞行时，升阻比越大，高超声速滑翔飞行器再入时高度下降越慢；在相同飞行仿真时间内，升阻比越大，航程越远，速度下降越慢，动能越大。从图13可知，本文改进模型高升阻比飞行时，航迹角波动相对较小。另外，如图14 所示大升阻比飞行时，对应攻角较小，这与图9中展现的结果相一致。图15结果可以发现，升阻比η分别取3.2、3.0、2.7时，路径约束都能满足条件，但大升阻比对应的这些路径约束值相对较大。

图10 不同升阻比飞行时高度变化Fig.10 Altitude changes with different lift-to-drag ratio

图11 不同升阻比飞行时经度变化Fig.11 Longitude changes with different lift-to-drag ratios

图12 不同升阻比飞行时速度变化Fig.12 Speed changes different lift-to-drag ratios

图13 不同升阻比飞行时航迹角变化Fig.13 Track angle changes with different lift-to-drag ratios

图14 不同升阻比飞行时攻角变化Fig.14 Angle of attack changes with different lift-to-drag ratios

图15 不同升阻比飞行时路径约束变化Fig.15 Path constraints changs with different lift-to-drag ratios

从上述对定升阻比仿真结果分析来看，大升阻比有利于大航程飞行，但对应的攻角较小，会导致动压、过载和热流率较大，因此，在再入飞行的开始阶段，采用小升阻比，可以在满足约束的条件下，使得攻角的控制范围变的更大，有利于高超声速滑翔飞行器再入轨迹的优化。

3.3 本文气动系数拟合模型再入轨迹优化

采用本文提出的气动系数拟合模型，针对CAV-H飞行器，进行高超声速滑翔飞行器再入轨迹优化仿真实验，并与文献[17-18]气动系数模型进行对比。CAV-H飞行器质量为907.2 kg，参考横截面积为0.484 m2，再入滑翔飞行实验的任务参数见表8，路径约束中最大热流率为4000 kW/m2，最大负载为2.5g(g为重力加速度)，最大动压为60 kPa。实验在Matlab2018b编程环境下进行，采用GPOPS2软件进行轨迹求解。

表8 CAV-H飞行器再入飞行任务数据Table 8 CAV-H aircraft reentry mission data

图16～19为三种模型在CAV-H再入滑翔飞行任务中得到的结果。图16～18可以看出，三种模型都可以按预定飞行状态到达终点，且到达终点时间几乎一致。但仔细观察图16和图17，可以发现，在各飞行阶段，采用本文气动系数模型的飞行器跳跃飞行的高度变化幅值，相对其他两种模型较小，速度的下降也比文献[17-18]中的模型更为平稳，从而带来动压、过载、热流率的变化整体上幅值波动相对较小，对飞行器结构强度、热防护措施等设计具有参考价值。