文章快速检索  
  高级检索
超/高超声速飞行器动态稳定性导数极快速预测方法
李正洲1,2, 高昌1, 肖天航2, 马自成1, 肖济良2, 朱建辉2     
1. 中国空气动力研究与发展中心 高超声速冲压发动机技术重点实验室, 绵阳 621000;
2. 南京航空航天大学 航空学院, 南京 210016
摘要: 飞行器设计早期阶段需要预测大量工况下的动导数。本文发展了一种面向超/高超声速飞行器的动导数极快速预测方法:首先基于当地流活塞理论,将飞行器进行小幅非定常运动所受到的气动力分为受自由来流引起的无附加扰动项以及受物面变形或运动引起的附加扰动项;通过当地表面斜度法、激波后等熵关系求解物面当地流动参数,进而结合非定常运动规律求出飞行器所受非定常气动力;再采用待定系数法对非定常气动力进行提取、辨识,最终得到超/高超声速飞行器动导数。该方法克服了传统方法对CFD流场参数的依赖和耦合,具有极高的计算效率;同时典型算例验证表明,该方法在超声速、高超声速工况下都能够很好预测动导数变化趋势。将该方法应用于复杂外形飞行器动导数预测,并讨论了与CFD方法的误差来源。本文方法可作为高速飞行器总体设计阶段布局选型的工具。
关键词: 动导数    强迫振动    动稳定性    当地流活塞理论    高超声速    超声速    
Extremely efficient prediction technique of dynamic derivatives for super/hypersonic flight vehicles
LI Zhengzhou1,2, GAO Chang1, XIAO Tianhang2, MA Zicheng1, XIAO Jiliang2, ZHU Jianhui2     
1. Science and Technology on Scramjet Laboratory, China Aerodynamics Research and Development Center, Mianyang 621000, China;
2. Aeronautics Engineering College, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
Abstract: In the early stage of aircraft design, it is necessary to predict the dynamic stability derivatives under a large number of operating conditions. An extremely efficient prediction technique of dynamic stability derivatives for high-speed flight vehicles is developed in this paper. Firstly, based on the local piston theory, the unsteady aerodynamics of high-speed flight vehicles can be divided into the non-disturbance aerodynamic component and the disturbance aerodynamic component. Secondly, the local surface inclination method and the isentropic flow relation after a shock wave are employed to solve the surface flow variables, and then the unsteady aerodynamic forces are calculated by combining the unsteady motion. Finally, the dynamic derivatives are extracted and identified from unsteady aerodynamic forces using the undetermined coefficient method. Our method overcomes the time-consuming of the traditional CFD method that is dependent and coupled with numerical simulation flow field variables. Typical test cases indicate that our method can predict the trend of dynamic derivatives well under both supersonic and hypersonic conditions. At last, our method is successfully applied to dynamic derivatives prediction of flight vehicles with complex shape, and the error sources with CFD method are discussed. The efficient dynamic derivatives prediction technique in this paper can be used to screen layout of high-speed flight in conceptual design stage.
Keywords: dynamic derivatives    forced vibration motion    dynamic stability    local piston theory    hypersonic    supersonic    

飞行器动态稳定导数是指飞行器绕重心转动时,其非定常气动力以及力矩系数相对于瞬时运动状态参数$ (p, q, r, \dot \alpha , \dot \beta ) $的导数,工程上简称“动导数”。动导数是飞行器姿态控制系统设计、弹道设计及动态稳定性分析中的关键参数[1]。随着现代飞行器飞行包线的扩展,设计人员对复杂飞行环境中飞行器的动态稳定性研究更加重视,对飞行器在各种飞行条件下的动导数都需要预先评估。对于超/高超声速飞行器,随着马赫数、迎角的增加,飞行器的动态稳定性问题变得更为复杂[2],气动布局设计早期阶段就必须考虑到飞行器的操稳特性[3]。因此,在超声速/高超声速飞行器的概念研究和初步设计阶段,高效、准确地评估飞行器动导数变得尤为重要。

目前获取动导数的主要方法有飞行试验、风洞试验[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所示:①所需输入为飞行器物面网格、来流参数及物面随时间的变形、振动历程;②飞行器的非定常气动力根据当地流活塞理论分为为受自由来流引起的无附加扰动项以及飞行器物面变形或运动引起的附加扰动项;③采用基于牛顿撞击理论的当地表面斜度法求出气动力无附加扰动项,再根据激波后的等熵关系求解当地流动参数,结合物面变形或振动可求出当地物面下洗速度,继而求出气动力附加扰动项随时间变化;④对气动力无附加扰动项和扰动项沿物面积分,即可求出高速飞行器受到的非定常气动力。图中:MaHαβ分别马赫数、高度、迎角、侧滑角;VlVb为飞行器物面的当地速度矢量;α0αmkt分别为强迫旋转运动的初始迎角、迎角振幅、缩减频率、时刻;hhm分别为强迫沉浮运动的高度、高度振幅。

图 1 高速飞行器动导数高效预测流程 Fig. 1 Efficient prediction process of dynamic derivatives for high-speed flight vehicles

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

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

2 非定常气动力快速预测 2.1 当地流活塞理论的推导

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

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

图 2 当地流活塞理论示意图 Fig. 2 Sketch of local piston theory

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

$ {\rm{d}}P \cdot S \cdot {\rm{d}}t = \rho aS{\rm{d}}t \cdot {\rm{d}}W $ (1)
 

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

$ {\rm{d}}P = {\rho _1}{a_1}{\rm{d}}W $ (2)
 

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

$ P = {\rho _1}{a_1}W + C $ (3)
 

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

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

$ \left\{ {\begin{array}{*{20}{l}} {P = {P_{{\rm{steady}}}} + {\rho _1}{a_1}W}\\ {W = {\mathit{\boldsymbol{V}}_1} \cdot {\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{n}} + {\mathit{\boldsymbol{V}}_{\rm{b}}} \cdot \mathit{\boldsymbol{n}}}\\ {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{n}} = {\mathit{\boldsymbol{n}}_0} - \mathit{\boldsymbol{n}}} \end{array}} \right. $ (4)
 

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

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

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

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

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

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

$ {C_p} = {C_{p,\max}}\left( {\frac{{1 - B}}{{0.32Ma}} + B \sin \theta } \right)\sin\theta $ (5)
 

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

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

$ {C_p} = \frac{{ - \left( {\gamma - 1} \right){\theta ^2}}}{2}\left[ {\sqrt {1 + {{\left( {\frac{4}{{\left( {\gamma + 1} \right)Ma\theta }}} \right)}^2}} - 1} \right] $ (6)
 

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

2.3 气动力附加扰动项计算

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

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

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

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

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

$ {P_2} = {P_\infty }\left( {\frac{{2\gamma Ma_\infty ^2}}{{1 + \gamma }} - \frac{{\gamma - 1}}{{\gamma + 1}}} \right) $ (7)
 
$ {\rho _2} = {\rho _\infty }\left[ {\frac{{\left( {1 + \gamma } \right)Ma_\infty ^2}}{{2 + \left( {\gamma - 1} \right)Ma_\infty ^2}}} \right] $ (8)
 

式中:P为来流静压;ρ为来流密度。

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

$ {\rho _{\rm{e}}} = {\left( {\frac{{{P_{\rm{e}}}}}{{{P_{\rm{2}}}}}} \right)^{\frac{1}{\gamma }}}{\rho _2} $ (9)
 

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

$ {a_{\rm{e}}} = \sqrt {\gamma \frac{{{P_{\rm{e}}}}}{{{\rho _{\rm{e}}}}}} $ (10)
 

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

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

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

$ {P_2} = {P_\infty } + {P_\infty }\frac{{2\gamma }}{{\gamma + 1}}\left( {Ma_\infty ^2{{\sin }^2}\beta - 1} \right) $ (11)
 
$ {\rho _2} = {\rho _\infty }\frac{{\left( {\gamma + 1} \right)Ma_\infty ^2{{\sin }^2}\beta }}{{\left( {\gamma - 1} \right)Ma_\infty ^2{{\sin }^2}\beta + 2}} $ (12)
 

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

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

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

3 气动导数提取及辨识

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

$ \alpha = {\alpha _0} + {\alpha _{\rm{m}}}\sin \left( {\omega t} \right) $ (13)
 

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

$ {C_m} = {C_{m\alpha }}\Delta \alpha + {C_{m\dot \alpha }}\Delta \bar {\dot \alpha} + {Q_{mq}}\Delta \bar q + {C_{m\ddot \alpha }}\Delta \bar {\ddot \alpha} + \cdots $ (14)
 

式中:Cm为俯仰力矩系数;$ \bar {\dot \alpha} $q分别为归一化的迎角变化率和俯仰角速度;C$ {C_{m\dot \alpha }} $Cmq分别为俯仰力矩系数对迎角、迎角变化率、俯仰角速度的导数。

显然,飞行器俯仰振荡过程中迎角变化率等于角速度,即$ \Delta \overline{\dot{\alpha}}=\Delta \bar{q} $,将其代入式(14), 并忽略高阶项,则有

$ {C_m} = {C_{m0}} + {C_{m\alpha }}\Delta \alpha + \left( {{C_{m\dot \alpha }} + {C_{mq}}} \right)\Delta \bar {\dot \alpha} $ (15)
 

式中:$\Delta \alpha=\alpha_{\mathrm{m}} \sin (\omega t) ; \Delta \overline{\dot{\alpha}}=k \alpha_{\mathrm{m}} \cos (\omega t)$;缩减频率k的表达式为

$ k = \frac{{\omega {L_{{\rm{ref}}}}}}{{2{V_\infty }}} $ (16)
 

式中:Lref为参考长度。

式(15)中(${C_{m\dot \alpha }}$+Cmq)即为所求的俯仰组合动导数。以括号形式表示组合动导数的原因是:该方法能够直接体现俯仰组合动导数,而非求出动导数的分量(洗流时差动导数${C_{m\dot \alpha }}$、直接阻尼动导数Cmq)再将其相加求和。

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

$ {C_m} = A + B \sin\left( {\omega t} \right) + C \cos\left( {\omega t} \right) $ (17)
 

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

$ \left( {{C_{m\dot \alpha }} + {C_{mq}}} \right) = C/\left( {k{\alpha _{\rm{m}}}} \right) $ (18)
 

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

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

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

$ {\mathit{\boldsymbol{V}}_\infty } = {V_\infty }\left[ {\cos\alpha \cos\beta ,\sin \beta ,\cos \beta \sin \alpha } \right] $ (19)
 

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

$ \begin{array}{l} \left( {{V_ \bot }/{V_\infty }} \right) = {n_x} \cos \alpha \cos\beta + {n_y}\sin\beta + \\ \;\;\;\;\;{n_z}\cos \beta \sin \alpha \end{array} $ (20)
 

式中:nxnynz为物面法向矢量的三分量;下标⊥表示速度在垂直物面方向的分量。

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

$ {\left( {{V_ \bot }/{V_\infty }} \right)_q} = \left( {x{n_z} - z{n_x}} \right)/l $ (21)
 

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

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

$ \begin{array}{l} {C_{{p_q}}} = {C_{p,\max }}\left[ {1/M{a_\infty } + 1.2\left( {{V_ \bot }/{V_\infty }} \right)} \right] \cdot \\ \;\;\;\;\;\;\;{\left( {{V_ \bot }/{V_\infty }} \right)_q} \end{array} $ (22)
 

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

$ {\rm{d}}{M_{yq}} = {\rm{d}}{F_{aq}} \cdot z - {\rm{d}}{F_{nq}} \cdot x $ (23)
 

式中:

$ {\rm{d}}{F_{aq}} = {C_{{p_q}}} \cdot {n_x}{\rm{d}}s,{\rm{d}}{F_{nq}} = {C_{{p_q}}} \cdot {n_z}{\rm{d}}s $ (24)
 

同样,可推导出面元ds上力矩系数的加速度导数$ {\rm{d}}{M_{y\dot \alpha }} $。沿物面对加速度导数和角速度导数同时积分,可得飞行器的俯仰组合动导数为

$ \left( {{C_{m\dot \alpha }} + {C_{mq}}} \right) = \frac{{\sum {\left( {{\rm{d}}{M_{y\dot \alpha }}} \right)} + \sum {\left( {{\rm{d}}{M_{yq}}} \right)} }}{{l \cdot S}} $ (25)
 

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

4 典型算例验证

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

4.1 超声速工况算例

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

图 3 有翼导弹几何外形示意图 Fig. 3 Sketch of basic finner missile geometry

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

图 4 有翼导弹面元网格 Fig. 4 Surface mesh of BFM

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

图 5 有翼导弹动导数预测结果对比 Fig. 5 Comparison of dynamic derivatives prediction results for BFM

表 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],通常被用来衡量高超声速工况下动导数计算结果的准确程度。

图 6 弹道外形示意图 Fig. 6 Sketch of hyperballistic shape

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

图 7 弹道外形面元网格 Fig. 7 Surface mesh of hyperballistic shape

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

图 8 弹道外形动导数预测结果对比 Fig. 8 Comparison of dynamic derivatives prediction results for HBS

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

表 2 不同网格密度对计算结果的影响 Table 2 Influence of different mesh density on calculation results
网格类型 网格量 动导数计算结果
粗糙 ~4 000 -31.9341
中等 ~20 000 -31.6621
加密 ~100 000 -31.1773

表 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飞行器示意图。

图 9 类X-37B飞行器示意图 Fig. 9 Sketch of X-37B analog

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

图 10 类X-37B飞行器动导数预测结果对比 Fig. 10 Comparison of dynamic derivatives prediction results for X-37B analog

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

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

图 11 边界层外缘密度对比 Fig. 11 Comparison of boundary layer edge densities
图 12 边界层外缘声速对比 Fig. 12 Comparison of boundary layer edge sound speeds
6 结论

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

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

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

参考文献
[1] 刘绪, 刘伟, 柴振霞, 等. 飞行器动态稳定性参数计算方法研究进展[J]. 航空学报, 2016, 37(8): 2348-2369.
LIU X, LIU W, CHAI Z X, et al. Research progress on the numerical method of dynamic stability derivatives[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2348-2369. (in Chinese)
Cited By in Cnki (16) | Click to display the text
[2] 祝立国, 赵俊波, 叶友达. 高速飞行器耦合失稳分析及应用[M]. 北京: 国防工业出版社, 2015: 71-76.
ZHU L G, ZHAO J B, YE Y D. Coupling departure analysis of high speed aircrafts[M]. Beijing: National Defence Industry Press, 2015: 71-76. (in Chinese)
[3] 闵昌万, 王颖. 高超气动外形设计的控制稳定性准则研究[M]. 北京: 科学出版社, 2019: 4-5.
MIN C W, WANG Y. Research on control stability criteria of aerodynamic shape design for hypersonic vehicles[M]. Beijing: Science Press, 2019: 4-5. (in Chinese)
[4] 郭雷涛. Φ1米高超声速风洞动导数试验技术研究[D].绵阳: 中国空气动力研究与发展中心, 2013: 3-8.
GUO L T. Investigation on dynamic derivative test technique in Φ1 m hypersonic wind tunnel[D]. Mianyang: China Aerodynamics Research and Development Center, 2013: 3-8(in Chinese).
Cited By in Cnki (6) | Click to display the text
[5] TOMEK D M, SEWALL W G, MASON S E, et al. The next generation of high-speed dynamic stability wind tunnel testing[C]//25th AIAA Aerodynamic Measurement Technology and Ground Testing Conference. Reston, VA: AIAA, 2006.
[6] GUO C, REN Y. The computation of the pitch damping stability derivatives of supersonic blunt cones using unsteady sensitivity equations[J]. Advances in Aerodynamics, 2019, 1(1): 1-16.
Click to display the text
[7] MI B G, ZHAN H, CHEN B. Calculating dynamic derivatives of flight vehicle with new engineering strategies[J]. International Journal of Aeronautical and Space Sciences, 2017, 18(2): 175-185.
Click to display the text
[8] 陈琦, 陈坚强, 袁先旭, 等. 谐波平衡法在动导数快速预测中的应用研究[J]. 力学学报, 2014, 46(2): 183-190.
CHEN Q, CHEN J Q, YUAN X X, et al. Application pf a harmonic balance method in rapid predictions of dynamic stability derivatives[J]. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(2): 183-190. (in Chinese)
Cited By in Cnki (25) | Click to display the text
[9] STERN E, GIDZAK V, CANDLER G V. Estimation of dynamic stability coefficients for aerodynamic decelerators using CFD[C]//AIAA Applied Aerodynamics Conference. Reston, VA: AIAA, 2012.
[10] HASHIMOTO A, MURAKAMI K, AOYAMA T, et al. Dynamic stability analysis of a reentry lifting capsule with detached eddy simulation[C]//54th AIAA Aerospace Sciences Meeting. Reston, VA: AIAA, 2016.
[11] TOBAK M, WEHREND W R. Stability derivatives of cones at supersonic speeds: NACA-TN-3788[R]. Washington, D.C.: NASA, 1956.
[12] LUCIA D J, BERAN P S, SILVA W A. Reduced-order modeling:new approaches for computational physics[J]. Progress in Aerospace Sciences, 2004, 40(1-2): 51-117.
Click to display the text
[13] MORGAN H G, RUNYAN H GL, HUCKEL V. Theoretical considerations of flutter at high Mach numbers[J]. Journal of the Aeronautical Sciences, 1958, 25(6): 371-381.
Click to display the text
[14] 陈劲松, 曹军. 超声速和高超声速翼型非定常气动力的一种近似计算方法[J]. 空气动力学学报, 1990, 8(3): 339-343.
CHEN J S, CAO J. An approximate calculating method of supersonic/hypersonic unsteady aerodynamic forces of airfoils[J]. Acta Aerodynamica Sinica, 1990, 8(3): 339-343. (in Chinese)
Cited By in Cnki (62) | Click to display the text
[15] 张伟伟, 史爱民, 王刚, 等. 结合定常CFD的当地流活塞理论[J]. 西北工业大学学报, 2004, 22(5): 546-549.
ZHANG W W, SHI A M, WANG G. On determining unsteady aerodynamic loads accurately and efficiently[J]. Journal of Northwestern Polytechnical University, 2004, 22(5): 546-549. (in Chinese)
Cited By in Cnki | Click to display the text
[16] ZHANG W W, YE Z Y, ZHANG C A, et al. Supersonic flutter analysis based on local piston theory[J]. AIAA Journal, 2009, 47(10): 2321-2328.
Click to display the text
[17] 刘溢浪, 张伟伟, 田八林, 等. 一种超音速高超音速动导数的高效计算方法[J]. 西北工业大学学报, 2013, 31(5): 824-828.
LIU Y L, ZHANG W W, TIAN B L, et al. Effectively calculating supersonic and hypersonic dynamic derivatives[J]. Journal of Northwestern Polytechnical University, 2013, 31(5): 824-828. (in Chinese)
Cited By in Cnki (5) | Click to display the text
[18] 秦之轩, 史爱明, 安富强, 等. 典型超声速/高超声速动导数计算方法研究[J]. 飞行力学, 2016, 34(3): 17-20.
QIN Z X, SHI A M, AN F Q, et al. Computing dynamic derivatives for supersonic and hypersonic models based on local piston theory[J]. Flight Dynamics, 2016, 34(3): 17-20. (in Chinese)
Cited By in Cnki | Click to display the text
[19] 刘绪, 刘伟, 周云龙, 等. 吸气式内外流一体化飞行器动导数数值模拟[J]. 空气动力学学报, 2015, 33(2): 147-155.
LIU X, LIU W, ZHOU Y L, et al. Numerical simulation of dynamic derivatives for air-breathing hypersonic vehicle[J]. Acta Aerodynamica Sinica, 2015, 33(2): 147-155. (in Chinese)
Cited By in Cnki (5) | Click to display the text
[20] 史晓鸣.基于当地流活塞理论的全机组合体颤振及气动伺服弹性分析[D].上海: 复旦大学, 2011: 13-14.
SHI X M. Flutter analysis of wing-fuselage complete vehicle based on local piston theory[D]. Shanghai: Fudan University, 2011: 13-14(in Chinese).
Cited By in Cnki (4) | Click to display the text
[21] LIU W, ZHANG C A, HAN H Q, et al. Local piston theory with viscous correction and its application[J]. AIAA Journal, 2016, 55(3): 942-954.
Click to display the text
[22] LI Z Z, XIAO T H, LV F X, et al. A rapid analysis tool for aerodynamics/aerothermodynamics of hypersonic vehicles[J]. Transactions of Nanjing University of Aeronautics&Astronautics, 2017, 34(4): 398-404.
Click to display the text
[23] BONET J, PERAIRE J. An alternating digital tree (ADT) algorithm for 3D geometric searching and intersection problems[J]. International Journal for Numerical Methods in Engineering, 1991, 31(1): 1-17.
Click to display the text
[24] 张鲁民, 叶友达, 纪楚群. 航天飞机空气动力学分析[M]. 北京: 国防工业出版社, 2009: 100-113.
ZHANG L M, YE Y D, JI C Q. Space shuttle aerodynamic analysis[M]. Beijing: National Defense Industry Press, 2009: 100-113. (in Chinese)
[25] ANDERSON J D. Hypersonic and high temperature gas dynamics[M]. 2nd ed. Reston, VA: AIAA, 2006: 40-41.
[26] 张堃元. 高超声速曲面压缩进气道及其反设计[M]. 北京: 国防工业出版社, 2018: 27-28.
ZHANG K Y. Hypersonic curved compression inlet and its inverse design[M]. Beijing: National Defense Industry Press, 2018: 27-28. (in Chinese)
[27] 彭科.飞行器气动力/热高精度快速计算方法及应用研究[D].长沙: 国防科技大学, 2016: 45-49.
PENG K. Accurate and rapid method of aircraft aero-force&aero-heating calculation and research on applications[D]. Changsha: National University of Defense Technology, 2016: 45-49(in Chinese).
Cited By in Cnki | Click to display the text
[28] 叶川, 马东立. 利用CFD技术计算飞行器动导数[J]. 北京航空航天大学学报, 2013, 39(2): 196-200.
YE C, MA D L. Aircraft dynamic derivatives calculation using CFD Techniques[J]. Journal of Beijing University of Aeronautics and Astronautics, 2013, 39(2): 196-200. (in Chinese)
Cited By in Cnki (21) | Click to display the text
[29] 刘伟, 沈清, 张鲁民, 等. 钝倒锥体动导数数值工程模拟[J]. 国防科技大学学报, 1998, 20(1): 5-8.
LIU W, SHEN Q, ZHANG L M, et al. Numerical and analytic simulation of the dynamic stability derivative of blunt cone[J]. Journal of National University of Defense Technology, 1998, 20(1): 5-8. (in Chinese)
Cited By in Cnki (6) | Click to display the text
[30] USELTON B L, JENKE L M. Experimental missile pitch and roll-damping characteristics at large angles of attack[J]. Journal of Spacecraft and Rockets, 1977, 14(4): 241-247.
Click to display the text
[31] OKTAY E, AKAY H U.CFD predictions of dynamic derivatives for missiles: AIAA-2002-0276[R]. Reston, VA: AIAA, 2002.
[32] 米百刚, 詹浩. 先进飞行器动导数数值模拟新方法[J]. 上海交通大学学报, 2016, 50(4): 619-624.
MI B G, ZHAN H. New calculation methods for dynamic derivatives of advanced flight vehicles[J]. Journal of Shanghai Jiaotong University, 2016, 50(4): 619-624. (in Chinese)
Cited By in Cnki (2) | Click to display the text
[33] EAST R A, HUTT G R. Comparison of predictions and experimental data for hypersonic pitching motion stability[J]. Journal of Spacecraft and Rockets, 1988, 25(3): 225-233.
Click to display the text
[34] 席柯, 阎超, 黄宇, 等. 俯仰阻尼导数分量的CFD数值模拟[J]. 北京航空航天大学学报, 2014, 41(2): 222-227.
XI K, YAN C, HUANG Y, et al. Numerical simulation of individual components of pitch-damping coefficient sum[J]. Journal of Beijing University of Aeronautics and Astronautics, 2014, 41(2): 222-227. (in Chinese)
Cited By in Cnki (9) | Click to display the text
[35] PAEZ C. The development of the X-37 re-entry vehicle[C]//40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. Reston, VA: AIAA, 2004.
[36] GRANTZ A C. X-37B orbital test vehicle and derivatives[C]//AIAA SPACE 2011 Conference&Exposition. Reston, VA: AIAA, 2011: 27-29.
[37] 张庆, 叶正寅. 基于气动导数的类X-37B飞行器纵向稳定性分析[J]. 北京航空航天大学学报, 2020, 46(1): 77-85.
ZHANG Q, YE Z Y. Dynamic characteristics analysis for X-37B like trans-atmospheric orbital test vehicle based on dynamic derivatives[J]. Journal of Beijing University of Aeronautics and Astronautics, 2020, 46(1): 77-85. (in Chinese)
Cited By in Cnki | Click to display the text
http://dx.doi.org/10.7527/S1000-6893.2019.23545
中国航空学会和北京航空航天大学主办。
0

文章信息

李正洲, 高昌, 肖天航, 马自成, 肖济良, 朱建辉
LI Zhengzhou, GAO Chang, XIAO Tianhang, MA Zicheng, XIAO Jiliang, ZHU Jianhui
超/高超声速飞行器动态稳定性导数极快速预测方法
Extremely efficient prediction technique of dynamic derivatives for super/hypersonic flight vehicles
航空学报, 2020, 41(4): 123545.
Acta Aeronautica et Astronautica Sinica, 2020, 41(4): 123545.
http://dx.doi.org/10.7527/S1000-6893.2019.23545

文章历史

收稿日期: 2019-10-08
退修日期: 2019-11-08
录用日期: 2019-11-19
网络出版时间: 2019-11-29 10:54

相关文章

工作空间