超/高超声速飞行器动态稳定性导数极快速预测方法

李正洲，高昌，肖天航，马自成，肖济良，朱建辉

1. 中国空气动力研究与发展中心 高超声速冲压发动机技术重点实验室，绵阳 621000

2. 南京航空航天大学 航空学院，南京 210016

目前获取动导数的主要方法有飞行试验、风洞试验[4-5]、数值计算[6-10]及工程近似方法等。飞行试验和风洞试验是飞行器动态特性判定的主要依据，但存在难度大、周期长、费用高的特点，需要与数值计算相互补充、互相验证；采用非定常CFD数值模拟方法，可以考虑到流场的非线性特性，便于开展复杂外形的气动力计算，但主要局限性在于计算量大，难以快速获得定性的结论和定量判据；以牛顿撞击理论[11]为代表的工程近似法考虑了线化空气动力学理论和经验关系，具有快捷高效的优势，但对复杂外形的大迎角非线性流动可能在数值上存在量级甚至符号的差别。

动导数的高效预测关键在于飞行器非定常气动力的快速、准确计算。在非定常气动力快速计算方面，近年来的热门方向是降阶模型(Reduced Order Modeling, ROM)[12]，这一方法在计算精度和计算效率上取得了很好的平衡。但由于降阶模型仍然部分依赖耦合非定常CFD数值计算，对非定常数值求解的鲁棒性有较高要求，因此一些学者转而采用定常CFD和工程方法相结合的思路构建非定常气动力模型，其中最典型的是基于CFD技术的当地流活塞理论(Local Piston Theory, LPT)[13]。陈劲松和曹军[14]、张伟伟[15-16]等基于定常欧拉(Euler)方程计算流场，利用获得的物面当地流动参数结合活塞理论来计算非定常气动力，这一方法解决了活塞理论只能计算一定马赫数范围内小迎角、尖头薄翼的缺点；刘溢浪[17]、秦之轩[18]等成功地将该方法应用于气动导数快速预测。

结合定常CFD与当地流活塞理论的方法来预测动导数，从结果来看，能够较好地揭示超声速/高超声速流动下飞行器动导数的变化规律，相比于完全非定常CFD时域推进方法可以节约大量计算时间。然而，该方法仍然依赖定常CFD流场数据，无法在飞行器设计早期阶段快速地给出定性的结果。本文针对上述问题，发展了一种结合当地表面斜度法和当地流活塞理论的高速飞行器非定常气动力快速计算方法，再对非定常气动力进行提取、辨识，可得到动导数。该方法不依赖CFD流场结果，在预测动导数方面具有极高的效率，同时具有较高的精度。

1 动导数高效预测策略及流程

本文动导数高效预测策略及流程如图1所示：① 所需输入为飞行器物面网格、来流参数及物面随时间的变形、振动历程；② 飞行器的非定常气动力根据当地流活塞理论分为为受自由来流引起的无附加扰动项以及飞行器物面变形或运动引起的附加扰动项；③ 采用基于牛顿撞击理论的当地表面斜度法求出气动力无附加扰动项，再根据激波后的等熵关系求解当地流动参数，结合物面变形或振动可求出当地物面下洗速度，继而求出气动力附加扰动项随时间变化；④ 对气动力无附加扰动项和扰动项沿物面积分，即可求出高速飞行器受到的非定常气动力。图中：Ma、H、α、β分别马赫数、高度、迎角、侧滑角；Vl、Vb为飞行器物面的当地速度矢量；α0、αm、k、t分别为强迫旋转运动的初始迎角、迎角振幅、缩减频率、时刻；h、hm分别为强迫沉浮运动的高度、高度振幅。

图1 高速飞行器动导数高效预测流程

本文方法的优势在于：采用当地表面斜度法及激波后的等熵关系式等求解飞行器表面流场参数，从而克服了传统方法对CFD流场参数的依赖和耦合(图1中虚线所示)。由于该方法不需要CFD方法求解定常流场，因此具有极高的效率；后文的算例也表明该方法具有较高的精度。因此本文方法可作为高速飞行器总体设计阶段布局选型的工具。

值得说明的是，虽然本文在预测飞行器气动导数过程中采用了基于牛顿撞击理论的当地表面斜度法，但与“采用牛顿撞击理论直接计算动导数”的工程近似法，两者的原理不同：后者是将压力系数对角速度(或迎角变化量)求导，再积分求出整个飞行器的动导数；本文方法是对飞行器施加强迫简谐运动求出周期非定常气动力，再对飞行器的非定常气动力进行提取、辨识气动导数。本文对飞行器施加强迫简谐运动为风洞试验、数值模拟计算动导数通用方法[19]。

2 非定常气动力快速预测

2.1 当地流活塞理论的推导

本文基于当地流活塞理论对超/高超声速飞行器非定常气动力进行快速预测，首先对当地流活塞理论进行推导。

采用动量定理推导当地流活塞理论[20]，如图2所示：考虑Ma≫1的气流经过飞行器表面，此时扰动可近似为沿物面法向传播，如同气缸中活塞队气流的扰动一样。图中：V、P、a分别为速度、压强、声速；下标∞、l、n分别表示自由来流、当地、法向的流动参数。

假设气流受到物面变形或振动引起的扰动后速度在dt时间内变化dW，扰动传播的距离为a·dt，则受到扰动气体质量为ρaSdt(ρ为气体密度，S为活塞面积)，受到扰动气体的总动量变化为ρaSdt·dW。另外，活塞对气体的冲量为dP·S·dt。根据动量定理可得

dP·S·dt=ρaSdt·dW

(1)

式(1)在当地流动参数下可简化为

dP=ρlaldW

(2)

对式(2)积分，可得当地流活塞理论的压力计算公式为

P=ρlalW+C

(3)

式中：常数C为附加扰动为0时的压力；ρlalW为物面扰动时的压力。

将附加扰动为0时的物面压力用Psteady表示，又物面扰动可分为变形或振动，则当地流活塞理论的压力计算公式可写为

(4)

式中：n0为物面变形前的外法线单位矢量；n为物面变形后的外法线单位矢量；Vl和Vb分别为当地流动速度和物面振动速度；物面扰动速度W由物面变形速度Vl·δn和振动速度Vb·n组成。

图2 当地流活塞理论示意图

与经典活塞理论相比，基于动量定理推导出的当地流活塞理论：① 推导过程中不需要级数展开，也就不存在忽略高阶项的近似，因此精度比忽略高阶项的活塞理论高；② 没有活塞理论所必须的扰动速度小于来流声速(W/a∞<1)的要求，因此当地流活塞理论只需满足“扰动可近似为沿物面法向传播”的假设即可，拓宽了马赫数适用范围的上限；③ 不需要等熵假设，这表明当地流活塞理论不仅可以用于大迎角、大厚度翼面问题，也可用于三维翼面、机身以及钝前缘或钝头体存在弓形激波的情况。上述特点是当地流活塞理论可推广用于机身和其他复杂外形的基础[21]。

式(4)表明超/高超声速飞行器非定常气动力可分为无扰动项和物面变形/振动后的扰动项。下面分别介绍这两项气动力的计算方法。

2.2 气动力无附加扰动项计算

对于无物面变形或振动扰动的高速飞行器，可看做定常流动，再用基于牛顿撞击理论的当地表面斜度法计算气动力无附加扰动项Psteady。具体方法为[22]：将飞行器表面离散、划分为三角形面元网格，对每个面元区分为迎风面或背风面，再分别采取不同的公式进行计算。迎风面、背风面的划分是根据自由来流条件及面元法矢共同确定的；此外，在面元判断时加入了对“遮挡面元”的判断，即采用光线投射算法(Ray-casting Algorithm)[23]判断当前面元是否被遮挡。若为遮挡面元，则将该面元作为背风面处理。

对于迎风面，采用活塞理论与修正牛顿理论结合的方法进行计算[24]：

(5)

式中：B=(2/π)arctan[(Ma2+6)/10]；Cp,max为驻点压力系数。

对于背风面，则采用普朗特-迈耶膨胀波方法：

(6)

式中：θ和γ分别为物面偏折角、气体比热比。

2.3 气动力附加扰动项计算

求解气动力附加扰动项时，需要用到当地流动参数。为求解方便，本文采用求解边界层外缘的扰动压力系数而不直接求解物面扰动压力系数。又根据边界层内的压力梯度较小，可近似认为边界层外缘的扰动压力等于飞行器表面的压力，因此所求出有效外形的非定常气动力即为飞行器的非定常气动力。

当地流动参数求解原理如下：对于完全气体或平衡气体，所有状态参数(压力、密度、温度、焓、声速、熵、流速、黏性系数)中只有两个是相互独立的，其他状态参数均可根据相应的热力学关系由这两个参数求出。

2.3.1 正激波完全气体边界层外缘参数

对于大钝头体飞行器形成的正激波完全气体，可认为进入边界层外缘的流线基本是通过弓形激波的正激波部分。由于正激波后的熵在来流条件给定后是唯一确定的，因此，利用所得出的物面定常压强和正激波后的熵这两个独立的状态参数，可以确定其他边界层外缘参数。

激波后的气体压力P2和密度ρ2可利用正激波前后关系式求得，即

(7)

(8)

式中：P∞为来流静压；ρ∞为来流密度。

利用等熵关系，求边界层外缘密度ρe：

(9)

再用边界层外缘的压力和密度求出边界层外缘的声速：

(10)

求出边界层外缘的当地密度和声速以后，结合物面的变形或振动的下洗速度，即可根据式(4)计算出气动力的附加扰动项。

2.3.2 斜激波完全气体边界层外缘参数

对于尖锐前缘、细长体飞行器外形，通过其头部的激波是一道或数道斜激波，斜激波完全气体后的压力和密度可分别为

(11)

(12)

式中：β为激波角，可由θ-β-Ma图表查出[25]，也可通过代数变换后的激波角显示表达式求出[26]。

求出激波后的气体压力P2和密度ρ2后，再采用正激波完全气体相同的方法求出边界外缘的密度ρe和声速ae。

对高温平衡气体而言，边界层参数的确定要比完全气体复杂。波后压力P2、焓H2和密度ρ2相互依赖，不能用波前参数简单求出，需要用迭代方法求解[27]。

3 气动导数提取及辨识

以求解高速飞行器俯仰方向组合动导数为例，求解组合动导数，通常对飞行器施加小幅强迫简谐定轴振动的运动形式。强迫简谐运动方程为

α=α0+αmsin(ωt)

(13)

根据线化气动力理论，俯仰力矩系数的泰勒展开表达式为

(14)

(15)

(16)

式中：Lref为参考长度。

由于式(15)中包含间接量Cm0、Cmα，因此无法直接求出组合动导数。文献[28]给出了一种求解方法，将俯仰力矩系数表达式改写为

Cm=A+Bsin(ωt)+Ccos(ωt)

(17)

式中：A、B、C为待定系数。在俯仰振荡非定常周期计算后，可得到俯仰力矩系数随时间的变化曲线，在曲线上任取3组不同时刻的数据即可求出A、B、C这3个系数。结合式(15)与式(17)，从而，俯仰组合动导数可表示为

(18)

式(18)即为本文结合当地流活塞理论与当地表面斜度法，采用强迫振动的运动方式求解俯仰组合动导数的公式。

为比较本文方法与“采用牛顿撞击理论直接计算动导数”的工程近似法之间的差异，下面给出工程近似法求动导数的方法[24,29]。

飞行器在体轴坐标系下的来流速度为

V∞=V∞[cosαcosβ,sinβ,cosβsinα]

(19)

物面的无量纲外法向速度投影可通过向量积求出，即

(V⊥/V∞)=nxcosαcosβ+nysinβ+

nzcosβsinα

(20)

式中：nx、ny、nz为物面法向矢量的三分量；下标⊥表示速度在垂直物面方向的分量。

在俯仰方向，由角速度引起物面与气流相对运动的导数为

(V⊥/V∞)q=(xnz-znx)/l

(21)

式中：l表示物面与飞行器参考点之间的距离。

则根据牛顿撞击理论，物面压力系数的角速度导数为

Cpq=Cp,max[1/Ma∞+1.2(V⊥/V∞)]·

(V⊥/V∞)q

(22)

面元ds上力矩系数的角速度导数为

dMyq=dFaq·z-dFnq·x

(23)

式中：

dFaq=Cpq·nxds,dFnq=Cpq·nzds

(24)

(25)

式(25)即为“采用牛顿撞击理论直接计算动导数”的工程近似法。

4 典型算例验证

分别在超声速和高超声速来流条件下选取典型算例对本文方法进行验证，以考核本文方法在不同速域范围内的适应性。

4.1 超声速工况算例

在超声速工况下，对有翼导弹(Basic Finner Missile，BFM)的俯仰组合动导数进行分析，目的是验证本文方法超声速工况预测动导数的精度以及效率。BFM是国际上验证动导数计算的标准模型[30-31]，有比较成熟的实验结果和工程估算结果。如图3所示，BFM为尖锥形头部、圆柱形弹身并带有4个矩形尾翼的外形，尾翼为十字布局，图中d为头部直径。

图4为计算所用的有翼导弹面元网格，总网格量约为3.5×104。该算例的计算条件为：Ma=1.96，以头部直径为参考长度的雷诺数Red=0.187×106，质心位置取外形纵向总长的61%处。

图5为不同方法对有翼导弹的动导数预测结果对比，其中风洞试验数据来源于文献[30]；CFD数据为文献[32]采用非定常CFD差分法的结果；准定常方法为对飞行器施加准定常定轴旋转运动的工程近似结果。从图中可以看出：在超声速工况下本文方法能够较好地反映有翼导弹动导数随迎角变化规律，精度也比工程近似方法有较大的提高。

图3 有翼导弹几何外形示意图

图4 有翼导弹面元网格

图5 有翼导弹动导数预测结果对比

表1为不同方法动导数计算效率对比。其中：非定常周期计算、非定常差分法预测动导数的耗时来源于文献[32]；本文方法在Intel Core i7-8700 3.20 GHz单核上运行，仅需约5 s就能获得所有工况的结果，表明本文方法具有极高的效率。特别的是，本文方法仅需要提供面元网格、不需要生成体网格，这节约了复杂外形体网格生成工作量与时间。

表1 不同方法动导数计算效率对比

Table 1 Efficiency comparison of different methods for dynamic derivatives computation

方法耗时非定常CFD周期计算[32]61 h非定常CFD差分法[32]12 h本文方法～5 s

4.2 高超声速工况算例

在高超声速工况下，选择弹道外形(HyperBallistic Shape, HBS)进行俯仰组合动导数预测，目的是验证本文方法高超声速工况下预测动导数的精度。

HBS标模是AGARD用来测定高超声速飞行器稳定性的典型外形。如图6所示，HBS由三段长均为1.5d′(d′为球体直径)的球柱、5°半锥角的圆锥段和15°半锥角的圆锥段组成，其俯仰动导数具有较为精确的试验结果[33]，通常被用来衡量高超声速工况下动导数计算结果的准确程度。

图7为弹道外形面元网格，总网格量约为2×104。 该算例的计算条件为：Ma=6.85，以头部直径为参考长度的雷诺数为Red′=0.72×106，质心位置取外形纵向总长的72%处。

图6 弹道外形示意图

图7 弹道外形面元网格

图8为不同方法对弹道外形的动导数预测结果对比，其中风洞试验数据来源于文献[33]；CFD数据为文献[34]采用非定常CFD方法的结果；准定常方法为对飞行器施加准定常定轴旋转运动的工程近似结果。从图中可以看出：在高超声速工况下本文方法能够较好地反映动导数随迎角变化规律，而工程近似方法只能反映动导数的量级。

为了考察网格密度对本文方法的影响，将本算例模型分别划分为粗糙、中等、加密3种不同密度的网格，以比较不同网格量下的动导数预测结果。不同网格在0°迎角下的动导数对比如表2所示。

图8 弹道外形动导数预测结果对比

表2 不同网格密度对计算结果的影响

从表2可以看出，在面元网格能够描述基本飞行器气动外形的前提下，网格量变化对本文方法影响十分微小(不超过0.8%)。

在计算效率方面：使用非定常的双时间方法计算本算例单个工况的俯仰动导数耗时为16 h，结合CFD和当地流活塞理论的方法耗时约为40 min[18]，本文方法计算单个工况的动导数平均不超过1 s，表明了本文方法具有极高的效率。

5 复杂外形飞行器气动导数预测

为了验证本文方法对复杂外形飞行器动导数预测的适用程度，以类X-37B飞行器为研究对象，分别采用本文方法和基于定常CFD的当地流活塞理论方法预测类X-37B飞行器的动导数，并与施加强迫振动的非定常CFD方法作对比。

X-37B是美国为了验证可重复使用空间技术和在轨空间飞行任务而启动的项目[35]，并计划在X-37B飞行器技术基础上继续进行能够投送6名宇航员进入太空的X-37C计划[36]。因此，以类X-37B飞行器为对象研究高速飞行器气动导数具有典型的意义。图9为类X-37B飞行器示意图。

图10为类X-37B飞行器动导数预测结果对比，其中CFD方法为文献[37]对飞行器施加小幅强迫振荡的非定常CFD结果；“无黏CFD+当地流活塞流理论”为基于定常CFD流场数据，采用当地流活塞理论计算非定常气动力，再对动导数进行提取、辨识方法得到的结果。通过图10可以看出：本文方法很好地预测出了类X-37B飞行器的动导数随迎角变化规律，精度与“无黏CFD+当地流活塞流理论”方法相当，但本文方法效率远高于后者。

图9 类X-37B飞行器示意图

图10 类X-37B飞行器动导数预测结果对比

由于“无黏CFD+当地流活塞流理论”提取、辨识动导数方法与本文方法相同，因此两者结果的误差来源可能是无黏CFD流场与本文“当地表面斜度法+激波后等熵关系”计算出的流场不一致造成。如式(4)所示：在使用当地流活塞理论计算非定常气动力时，用到了两个重要的流场参数：边界层外缘密度及边界层外缘声速。因此，本文对边界层外缘密度及边界层外缘声速进行对比分析。

图11、图12分别为上述两种方法得出的边界层外缘密度与边界层外缘声速对比(α=30°)。从图中可以看出，两种方法计算出的流场大面积分布较为一致，但在驻点、翼前缘附近有较大差别。这应是两种方法预测动导数误差的来源。

图11 边界层外缘密度对比

图12 边界层外缘声速对比

6 结 论

1) 本文动导数预测方法是基于“当地流活塞理论”及“牛顿撞击理论”等相关理论推导而出，因此“当地流活塞理论”及“牛顿撞击理论”的成立条件决定了本文方法适用的速域范围，即扰动可近似为沿物面法向传播的超声速/高超声速流动。

2) 当地流活塞理论推导过程中不需要等熵假设，拓宽了经典活塞理论对飞行器外形和迎角的适用范围，因此本文方法可用于复杂外形飞行器的动导数预测。

3) 本文方法避免了传统动导数预测方法对CFD流场的依赖、耦合，从而大幅提高了计算效率；同时在面对复杂外形时也能较好地得出动导数随迎角的变化规律，精度能够满足飞行器总体设计阶段的要求，可作为飞行器布局选型阶段的工具。