与传统的飞行器不同，高超声速飞行器表现飞行出多任务、多工作模式、大范围高速机动的特点^[1-3]。此外高超声速飞行器飞行航程较长并且发动机燃油消耗量大，进而增加起飞规模，所以合理设计高超声速飞行器上升段飞行方案，可以保证飞行器爬升阶段的可控性以及合理的油耗规模。

目前，对于高超声速飞行器飞行轨迹设计的研究主要针对飞行剖面中的爬升段。因为这一阶段是保证飞行器顺利达到预设巡航点的基础，并且飞行时间较长、涉及到的约束条件较多。所以，合理设计高超声速飞行器爬升轨迹可以降低飞行成本和对飞行器各个子系统的要求^[4]。针对高超声速飞行器爬升过程中的不同约束以及性能指标条件，许多学者对爬升轨迹设计与优化进行了研究。采用定动压来模拟飞行器加速爬升过程，以平衡推力和质量变化之间的关系。将航迹角变化率作为约束条件之一，并基于自身结构与燃料质量之间的等式关系进行轨迹优化以获得最省燃油轨迹^[5-9]。

飞行器的轨迹优化主要包括解析法和数值法两大类，目前数值法在高超声速飞行器的轨迹优化问题中有较多的应用^[10-12]。高超声速飞行器通常采用数值算法^[13]进行轨迹优化求解，因为在高超声速飞行器的方程维度较高且模型复杂。能量状态法通过引入能量，对飞行器的运动学方程进行简化，减低了对高超声速飞行器进行轨迹优化时的模型维度和复杂性。David等^[14-15]采用能量状态法对高超声速飞行器爬升段轨迹进行优化设计。

但是传统能量状态法进行轨迹优化设计时仅考虑采用高超声速飞行器的质点模型，而高超声速飞行器的机身/发动机耦合设计，使得姿态对轨迹的影响不可忽视。所以，本文考虑了高超声速飞行器的姿态运动，将速度、高度、轨迹角和质量作为状态变量，将迎角和舵面作为控制变量，对能量状态法进行改进。

1 飞行器/发动机一体化模型 1.1 发动机数学模型

由于高超声速飞行器的水平起降能力，使得其具有高可用率，在军事和民用上均具有较高实用价值。其常采用涡轮基组合动力(turbine based combined cycle, TBCC)发动机，如图 1所示。TBCC发动机从起飞就处于涡轮模式，在马赫数为2.5~3内发动机转换为冲压发动机模式，然后由亚燃/超燃冲压发动机驱动加速度到高超声速。

根据动量定理，发动机的未安装推力可以写为：

$ T = {\dot m_9}{V_9} - {\dot m_0}{V_0} + ({P_9} - {P_0}){A_e} $

(1)

式中发动机位置0为自由流，位置9为发动机出口，$\dot m{\;_*}$为质量流量(下标“*”为0或9)。安装推力的入口气流受飞行器前体压缩角影响。飞行马赫数和高度的组合可以由特定值θ表示：

$ \theta = {T_t}/{T_{{\rm{std}}}} = T/{T_{{\rm{std}}}}\left( {1 + \frac{{\gamma - 1}}{2}M_0^2} \right) $

(2)

安装推力为：

$ T = \phi \alpha (Ma,h){T_{{\rm{SL}}}} = \phi \alpha ({\theta _0}){T_{{\rm{SL}}}}(TR) $

(3)

式中：T_SL是可用输出推力与最大推力的无量纲比值; TR是油门开度；α表示发动机工作状况，它由马赫数(Ma)、高度(h)决定，可以定义为无量纲温度θ的函数。

当飞行条件为θ≤TR时，发动机受Π_c(隔离室允许动压)限制；当飞行条件为θ>TR时，发动机受T_t4(燃烧室出口允许温度)限制。

考虑到上述控制逻辑，可用推力可写为关于油门开度的分段多项式形式：

$ \left\{ {\begin{array}{*{20}{l}} {\alpha = {\delta _0},\quad \theta \le TR}\\ {\alpha = {\delta _0}[1 - 3.5(\theta - TR)/\theta ],\quad \theta > TR} \end{array}} \right. $

(4)

式中δ₀是无量纲压力比：

$ {\delta _0} = {P_t}/{P_{{\rm{std}}}} = P/{P_{{\rm{std}}}}{\left( {1 + \frac{{\gamma - 1}}{2}M_0^2} \right)^{\frac{\gamma }{{\gamma - 1}}}} $

(5)

由于高超声速飞行器机身推进一体化的设计，发动机偏置推力也会对飞行姿态产生影响。假设C_Tx和C_Tz为飞行器x轴向和z轴向的推力系数。图 2为在30 km处对应的不同飞行马赫数下的x轴向和z轴向的推力系数随迎角的变化趋势。图 3为在Ma=3下不同飞行高度的比冲随迎角的变化趋势。

根据图 2和图 3可知，高超声速飞行器的发动机推进系统不仅与飞行速度、高度相关，也受飞行迎角的影响。所以在使用能量状态法对高超声速飞行器进行分析时，不能仅考虑质点模型。

1.2 气动数学模型

由参数化建模方法^[4]可以快速获取高超声速飞行器的气动数据库，然后建立相应的代理模型。其中α∈[-15°, 15°]，δ_e∈[-15°, 15°]，Ma∈[1.2, 7]。所有气动系数可以分离为2个部分：1)洁净体C_i(clean)；2)由控制舵面偏转而产生的增量ΔC_i^δ_eδ_e。

$ C_{i(\text { total })}=C_{i(\text { clean })}+\Delta C_{i}^{\delta_{e}} \delta_{e} $

(6)

洁净体纵向气动系数的表达式为：

$ C_{\text {Lclean }} =\left(a_{0}+a_{1} M a\right)+\left(a_{2}+a_{3} M a\right) \alpha $

(7)

$ C_{D_{\text {clean }}}=\left(b_{0}+b_{1} M a\right)+\left(b_{2}+b_{3} M a\right) C_{\text {Lclean }}^{2} $

(8)

$ C_{m_{\text {clean }}} =\left(c_{0}+c_{1} M a\right)+\left(c_{2}+c_{3} M a\right) \alpha $

(9)

增量纵向气动系数的表达式为：

$ \begin{array}{c} C_{L}^{\delta_{e}}=a_{4}+a_{5} M a+a_{6} \alpha \\ C_{D}^{\delta_{e}}=b_{5}+b_{6} \delta_{e}+b_{7} M a+b_{8} \alpha+b_{9} M a^{2} \delta_{e} \\ C_{m}^{\delta_{e}}=c_{4}+b_{5} \alpha^{2}+b_{6} M a+c_{7} M a^{2} \end{array} $

(10)

选择拟合优度(goodness of fit, GOF)和均方根误差(root mean square error, RMSE)作为模型精确性的评价指标，如表 1所示，表明代理模型具有较高的精度。

2 能量状态法 2.1 方法简介

对于采用能量状态法对飞行器进行性能预测，垂直平面上运动的质点模型通常是足够的。该模型的运动学方程为：

$ \begin{array}{c} m\dot V = T\cos (\alpha + \varepsilon ) - D - mg\sin \gamma \\ \begin{array}{*{20}{c}} {mV\dot \gamma = T\cos (\alpha + \varepsilon ) + L - mg\cos \gamma }\\ {\dot h = V\sin \gamma }\\ {\begin{array}{*{20}{c}} {\dot x = V\cos \gamma }\\ {\dot m = - T/{g_0}{I_{sp}}} \end{array}} \end{array} \end{array} $

(11)

式中：V为飞行器速度；α为飞行器迎角；γ为轨迹角；ε为推力轴线与零升力轴线之间的夹角；g₀为标准重力加速度；I_sp为发动机比冲；x、h为地面坐标系中飞机质心的前向和纵向坐标；m为飞机质量。

因为α和ε均较小，所以进行如下简化：

$ \sin (\alpha + \varepsilon ) \cong \alpha + \varepsilon ,\cos (\alpha + \varepsilon ) \cong 1 $

(12)

对于在亚音速飞行的飞行器，忽略其飞行姿态进行轨迹设计是足够的。然而高超声速飞行器机身推进一体化的设计，发动机的工作情况会受飞行迎角的影响。因此将速度、高度、轨迹角和质量作为状态变量，将迎角和舵面作为控制变量，采用改进能量状态法对高超声速飞行器进行分析。

总能量可以看作是系统连续的状态变量。定义单位质量飞行器的能量为：

$ E = h + \frac{{{v^2}}}{{2g}} $

(13)

进而得飞行器质点运动的能量状态模型：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{{\rm{d}}E}}{{{\rm{d}}t}} = \frac{{v(T\cos \alpha - D)}}{{mg}}}\\ {\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = v\cos \alpha + {v_w}} \end{array}} \right. $

(14)

2.2 基于能量状态法的飞行包线分析

能量状态法假定飞行器纵向剖面由爬升、巡航、下降3种飞行任务阶段组成，其能量变化率分别为dE/dt>0、dE/dt=0和dE/dt < 0。所以对于高超声速的飞行包线对应的是飞行器能量状态变化率为0的飞行状态。

基于能量状态法的飞行包线计算步骤如下：

1) 在给定马赫数下，计算dE/dt=0，则飞行器需要满足下式的要求：

$ \left\{\begin{array}{l} T\left(V_{\text {given }}, h, \alpha\right) \cos \alpha-D\left(V_{\text {given }}, h, \alpha, \delta_{e}\right)-m g \sin \gamma=0 \\ T\left(V_{\text {given }}, h, \alpha\right) \sin \alpha+L\left(V_{\text {given }}, h, \alpha, \delta_{e}\right)-m g \cos \gamma=0 \end{array}\right. $

(15)

在迎角约束(0° < α < 8°)和升降舵面约束(|δ_e| ≤20°)的要求下，通过求解式(15)，获得2个高度值，dE/dt=0所对应的高度上界h_maxi和下界h_mini。当飞行器的飞行高度在区间[h_maxi, h_mini]内时，飞行器的能量变化率大于零，即可以进行爬升/加速运动；

2) 计算给定马赫数范围内，不同马赫数下dE/dt=0所对应的高度上下界，得出给定马赫数范围内的飞行包线。

根据上述步骤，计算马赫数在1.2~7的飞行包线范围，结果如图 4所示。

根据图 4所示的飞行包线结果可以看出发动机由涡轮模态转换为冲压模态时会出现一个“推力鸿沟”，会急剧缩小可飞行范围。高超声速飞行器在马赫数为2.5完成从涡轮模态向冲压模态的转换。此外，图 5所示推阻关系也与图 4的飞行包线范围相对应。

根据图 5可以得出在马赫数为1.2~7，飞行器在当前气动和推进系统状态下存在可飞行范围。根据图 5(b)可以看出，冲压发动机初始的净推力较小，随着马赫数增大冲压发动机逐渐进入较好的工作状态。并且高超声速飞行器在Ma为1.2的超声速初始飞行阶段的可行动压范围在5~25 kPa范围内，可飞行范围相对较窄，随着飞行速度增加，飞行包线范围也逐渐扩展。通过对比图 4和图 5所示的推阻关系，可以看出两者得出的结论一致。所以，采用能量状态计算高超声速飞行器的飞行包线范围具有一定的可行性。

3 基于能量状态法的节油轨迹设计

本文对传统能量状态法进行改进，以速度、迎角和升降舵偏转角为控制量，考虑飞行过程中的多约束等条件对高超声速飞行器的爬升轨迹进行寻优。

3.1 优化模型建立

能量状态法选取飞机的耗油量为优化性能指标，具体定义为：

$ J = \int\limits_{t0}^{{t_f}} {{{\dot m}_f}} {\rm{d}}t $

(16)

式中：${\dot m_f}$为燃油消耗率；t_f为终端时间。为了对计算过程进行进一步简化，将式(14)代入性能指标式(16)中，将状态变量E转化为性能指标泛函中的独立变量，得到：

$ \frac{\mathrm{d} x}{\mathrm{~d} E} =\frac{v \cos \alpha}{v(T \cos \alpha-D) / m g} $

(17)

$ J = \int_{{E_0}}^{{E_f}} {\frac{{{{\dot m}_f}}}{{v(T\cos \alpha - D)/mg}}} {\rm{d}}t $

(18)

对于整个飞行阶段，式(17)、(18)可以分解为：

$ \begin{array}{c} \frac{\mathrm{d} x}{\mathrm{~d} E}=\frac{\mathrm{d}\left(x_{\mathrm{up}}+x_{\mathrm{dn}}\right)}{\mathrm{d} E}= \\ \left.\frac{V_{\mathrm{up}} \cos \gamma+V_{\text {wup }}}{\dot{E}}\right|_{\dot{E}>0}+\left.\frac{V_{\mathrm{dn}} \cos \gamma+V_{\mathrm{wdn}}}{\dot{E}}\right|_{\dot{E}<0} \end{array} $

(19)

$ J=\left.\int_{E_{0}}^{E_{c}}\left(\frac{\dot{m}_{f}}{\dot{E}}\right) \mathrm{d} E\right|_{\dot{E}>0}+\lambda d_{c}+\left.\int_{E_{c}}^{E_{f}}\left(\frac{\dot{m}_{f}}{|\dot{E}|}\right) \mathrm{d} E\right|_{\dot{E}<0} $

(20)

式中：下标up、dn表示爬升和下降；V_W为风速；λ为巡航飞行中单位距离上的油耗量；d_c为巡航里程。此时，轨迹优化问题就转换为了求解满足约束(19)，使性能指标式(20)最小的优化问题。

3.2 爬升段寻优

总能量E在爬升段是单调增加的，从初始能量E₀开始，至末端能量E_f结束，增量为ΔE。采用能量状态法进行寻优的步骤如下。

1) 对每一能量水平E_j，计算最优爬升下的速度、迎角和升降舵偏转角如下：

$ \lambda_{E}=\min \limits_{V, \alpha, \delta_{e}}\left\{\frac{\dot{m}_{f} m g \Delta E}{v(T \cos \alpha-D)}\right\} $

(21)

计算相应的飞行高度公式为：

$ h_{j}=E_{j}-V_{j}^{2} / 2 g $

(22)

2) 能量水平E_j-1、E_j之间的飞行时间Δt_j，轨迹角γ_j，飞行距离Δd_j，以及燃油消耗量ΔF_j公式为：

$ \left\{\begin{array}{l} \Delta t_{j}=\Delta E / \dot{E}_{j} \\ \gamma_{j}=\arcsin \left[\left(\left(h_{j}-h_{j-1}\right) / \Delta t_{j}\right) /\left(V_{j}+V_{w j}\right)\right] \\ \Delta d_{j}=\left(V_{j}+V_{w j}\right) \Delta t_{j} \cos \gamma_{j} \\ \Delta F_{j}=f_{j} \Delta t_{j} \end{array}\right. $

(23)

式中${\dot E_j} = {\rm{d}}{E_j}/{\rm{d}}t$，为第j段的能量梯度，由式(17)的第1个方程给出。当设计轨迹达到末端状态时，结束寻优；否则，重复以上过程。

由于爬升段不允许飞机降低高度，故h_j≥h_j-1，即在能量水平E_j上，高度的最小允许值为h_{min, j}=h_j-1，根据能量方程，这个限制可以转换为相应的速度上限为：

$ {V_{{\rm{max}},h,j}} = \sqrt {2g({E_j} - {h_{j - 1}})} $

(24)

为了防止爬升高度超过末端高度h_f，高度的最大允许值为h_f，其相应的速度下限为：

$ {V_{{\rm{max}},h,j}} = \sqrt {2g({E_j} - {h_f})} $

(25)

除了飞行器的速度约束范围外，使用能量状态法寻优的过程中还需满足如下条件：动压约束30~100 kPa；高度约束，不允许飞行高度降低；推阻约束(T)>D；末端约束。

3.3 爬升轨迹寻优仿真 3.3.1 迎角影响分析

首先采用能量状态法分析飞行器迎角对最优轨迹和飞行燃油消耗产生的影响。基于改进能量状态法给出不同迎角下的最省燃油轨迹图(图 6)以及飞行器质量变化图(图 7)。

根据图 7和表 2中不同迎角下的燃油量可以得知，迎角对发动机的工况和燃油消耗情况影响较大。即使相差较小的迎角，也会产生较大的燃油消耗差异。根据图 6和图 7可以得出，飞行迎角对最优飞行轨迹影响不大，但会对高超声速飞行器飞行过程中的燃油消耗产生影响很大。

3.3.2 最优轨迹设计

在此采用改进能量法，考虑飞行过程中的动压约束、高度约束、推阻约束和末端约束等条件对高超声速飞行器的爬升轨迹进行优化设计。

采用改进能量状态法，获取最优的飞行高度和迎角，进而得到最省燃油轨迹，并且该轨迹满足飞行过程中的动压和飞行器舵面物理限制约束。寻优过程设计飞行器从Ma在1.2 ~6.6的最省燃油飞行轨迹。Ma在1.2~6.6，根据上升段寻优过程得到的最省燃油轨迹，如图 8所示。

由能量状态法获得的最省燃油爬升轨迹与沿100 kPa等动压爬升轨迹相比，飞行器沿该轨迹飞行可节省燃油1.87 t。此外，图 8所示的初始时刻沿30 kPa等动压爬升，与图 5中此时的最大净推力所对应动压相吻合。

4 结论

1) 采用基于能量状态法的飞行包线分析，为燃油最省爬升轨迹设计划定范围，以避免设计的轨迹在切换至冲压模态下处于可飞行包线以外。

2) 分析迎角对最优轨迹和燃油消耗的影响，以速度、迎角和升降舵偏转角为控制量，并且考虑飞行过程中的动压约束、高度约束、推阻约束和末端约束等条件的改进能量状态法，比常规能量状态法更适合高超声速飞行器的爬升轨迹优化设计。