2. 国家CFD实验室, 北京 100191
2. National Laboratory for CFD, Beijing 100191, China
0 引 言
航天再入飞行器出于防热的需求，通常设计成钝头体构型。这类飞行器再入时，由于要经历高超声速、超声速至亚跨声速的全流域飞行，特别是在中-低超声速阶段，俯仰运动将会出现动态不稳定现象[1]，表现为俯仰振动的振幅逐渐增加，最终形成极限环振荡。对再入飞行器动态失稳现象的研究始于20世纪50年代[2, 3]，如今六十多年过去了，动态失稳的物理机制以及动态特性的准确预测仍然困扰着研究人员[4]。美国“阿波罗号”载人飞船在研制过程中，虽然采用9套模型分别在14座风洞中进行了亚、跨、超和高超声速风洞试验，试验时间累计高达700多小时，但其动态特性问题仍没有完全解决。最近基于“阿波罗号”飞船研 制的“猎户座”CEV，在投放试验中由于多次出现动态失稳现象而失效，经过Langley的深入研究才得以解决，并于2014年12月5日飞行试验成功。
已有的研究表明，小升阻比的载人飞船返回舱外形以球 冠倒锥形为最优[5, 6, 7]。该类外形的飞船返回舱以较大的配平攻角再入，后体有大尺度的分离和再附，流场中存在复杂的波系结构(图 1)。钝体气流分离效应、后体气流再附效应、船尾近尾涡流效应 和动态迟滞效应等对飞船返回舱的静、动稳定性都有影响。
 |
| 图 1 返回舱再入时典型波系Fig. 1 Flow characteristics of reentry capsule |
因此，有研究者试图通过有限度的改变返回舱后体外形，以改变后体分离拓扑形态，从而改善返回舱的稳定性，如“联盟号”和“神舟号”载人飞船返回舱在光体倒锥上增加稳定翼[8]。当前已有的大量试验与计算研究，一般只能给出动态失稳现象的唯象描述，缺乏理论分析[9, 10, 11]。本文采用非线性自治动力系统分叉理论，耦合求解非定常Navier-Stokes方程和俯仰运动方程，研究了球冠倒锥飞船返回舱外形和平头有翼双锥外形的单自由度俯仰动态失稳问题，给出了Hopf分叉理论描述和数值计算结果，分析了两种外形发生不同Hopf分叉的机制[8]。 1 数学模型 1.1 平面自治动力系统
这里只研究轴对称飞行器绕重心作单自由度俯仰振荡的运动情况，其坐标系见图 2。
 |
| 图 2 单自由度俯仰运动的坐标系统Fig. 2 Coordinate system of pitching |
设ξ-η-ζ是与飞行器相固连的正交坐标系，ξ轴与飞行器的体轴重合，η、ζ两轴构成俯仰对称平面。 x-y-z是惯性静止坐标系，x轴与来流方向一致。设 αt为平衡攻角，θ(t)是由平衡攻角起算的俯仰振荡角，则瞬态攻角为α(t)=αt+θ(t)。此时，描述飞船返回舱俯仰振荡运动和非定常动态绕流以及所受俯仰力矩的无量纲化的耦合方程组为[3]：
式(1)为描述飞行器单自由度俯仰运动方程，I表示无量纲的转动惯量，右端项包括气动俯仰力矩系数Cm和机械阻尼力矩系数Cμ。在自由飞行中Cμ为0，在风洞实验中应当计及它的贡献。
分别表示俯仰角θ对时间的一阶和二阶导数。式(2)为非定常Navier-Stokes方程，空间离散采用原始变量NND格式，时间离散采用双时间步方法。式(3)为俯仰力矩系数的积分关系式。非定常计算方法的详细描述及方程中各符号的含义见文献[8]。
对任意非定常运动，Etkin认为[9]非定常气动力和力矩是状态变量的泛函。据此可给出动态俯仰力矩系数表达式：
在多数情况下，可只保留到一阶导数项，但在某些情况下，为了保证足够的精确度，必须保留到二阶导数项。已经证明，这一假设在实际应用中的适用范围极为广泛。将式(4)保留至二阶导数，并将无量纲的转动惯量I吸收到Cm和Cμ中，且吸收后的力矩系数和机械阻尼系数仍用Cm和Cμ表示，代入式(1)，易得：
将俯仰力矩系数在平衡攻角αt处展开，可得：
这里，下标0表示在平衡攻角处的定态俯仰力矩，故有(Cm)0=0。非线性项G(θ,)是关于θ、的高阶项，假定当(θ2+
)1/2→0时，非线性项G(θ,)比(θ2+)1/2更高阶地趋于0，即非线性项满足Perron定理[12]。将式(6)代入式(5)，即得：
令x=
,y=θ，代入上式就得到：
分别称为静导数、动导数和二阶动导数。飞行器一般是静稳定的，有
<0；二阶动导数的量值远小于1，也即
≤1。至此，就建立了描述单自由度俯仰运动的非线性平面自治动力系统式(8)。该动力系统的平衡点对应于飞行器在平衡攻角处的定态绕流，此时有(θ,)0=(0,0,0)。
1.2 Hopf分叉分析假设上述非线性系统式(8)的非线性项满足Perron定理，根据动力学方法[8, 12]，可以用非线性系统的一次近似系统式(9)来分析其平衡解的稳定性。
设：
则一次近似系统式(9)的Jacobian矩阵的特征值为：
。根据非线性动力学系统的有关定理，可知在p-q平面上各类平衡点(也称奇点或临界点)的区域如图 3。
 |
| 图 3 p-q平面上各类平衡点的区域图Fig. 3 Distribution of critical point in p-q plane |
现在研究平衡点(θ=
=0)处的动态稳定性及其相关问题。当Δ<0,参数λ由λ<0经过λ=0变化至λ>0时，一次近似系统的Jacobian矩阵的特征值为一对共轭复根，进一步可得到以下结论：
(1) 若λ<0，Δ<0，易得到：p>0，q>0，则在(x,y)相平面上，平衡点为稳定的螺旋点，也称焦点。在平衡点(0,0)附近，非线性动力系统式(8)的轨线为稳定的螺旋点形态(图 4(a))。俯仰振荡角的时间历程 曲线是收敛的，即随时间t增加，θ是减小的，最后趋于0(图 4(b))。因此，λ<0和Δ<0可作为平衡点的动稳定性的判据。
 |
| 图 4 平衡点动态稳定性示意图(λ<0)Fig. 4 Sketch of dynamic stabilization at balance angle of attack(λ<0) |
(2) 若λ>0，Δ<0，可得到：p<0，q>0，则在(x,y)相平面上，平衡点为不稳定的螺旋点，在平衡点(0,0)附近，非线性动力系统式(8)的轨线为不稳定的螺旋点形态(图 5(a))。俯仰振荡角的时间历程曲线是发散的，即随时间t增加，θ是增大的(图 5(b))。
 |
| 图 5 平衡点动态稳定性示意图(λ>0)Fig. 5 Sketch of dynamic stabilization at balance angle of attack(λ>0) |
(3) 若λ=0，Δ<0，可以得到：p=0，q>0，非线性动力系统式(8)的一次近似系统Jacobian矩阵的特征值 
，在λ=λcr=0时，满足：
① 特征值的实部Re λ1(λcr),λ2(λcr) =0；
② 特征值的虚部Im λ1(λcr),λ2(λcr) ≠0；
③ 因为Δ<0，则Re λ1(λ),λ2(λ) = 1/2 λ，
故有：
于是，在λ=λcr=0时，系统的特征值满足Hopf分叉的三个条件[12]，其中第三个条件称为Hopf分叉的横截条件(transversality condition)。即当λ由λ<0经过λ=0变化至λ>0时，动力系统式(8)将发生Hopf分叉，在(x,y)相平面上平衡点(0,0)失稳，出现稳定的极限环(图 6(a))。俯仰振荡角的时间历程曲线在达到稳定后既不收敛也不发散，即随时间t增加θ出现周期振荡(图 6(b))。
 |
| 图 6 平衡点动态稳定性示意图(λ=0)Fig. 6 Sketch of dynamic stabilization at balance angle of attack(λ=0) |
这样，理论上就获得了单自由度俯仰运动发生动态Hopf分叉失稳，出现极限环振荡的临界条件，即：
针对飞船返回舱外形和平头双锥外形，对上述Hopf分叉理论进行验证。对于每个外形，需要求解三类流场：首先，计算定态绕流流场，以获取平衡攻角，并为后续的非定常计算提供初场；然后，数值模拟强迫俯仰振荡过程，通过参数辨识，获取平衡攻角处的静、动导数[13, 14, 15, 16, 17, 18, 19, 20]；最后，计算自激俯仰振荡过程，验证Hopf分叉理论。 2.1 飞船返回舱数值模拟
日本的轨道再入实验飞船OREX是一个球冠加倒锥外形的旋成体，计算网格和外形见图 7。质心在体轴上，再入时只存在一个大头朝前的平衡攻角αt=0°。
 |
| 图 7 返回舱外形和计算网格Fig. 7 Configuration of the capsule and computational grid |
风洞试验和飞行试验均表明，该返回舱在再入经过跨、超声速阶段时，会发生俯仰方向的动态失稳现象，出现绕平衡攻角的准极限环振荡。这里选取的马赫数Ma变化范围为1.5～6.0，雷诺数均取Re=1.0×105，以单独考察俯仰动态特性随马赫数的演化规律。强迫俯仰运动给定为简谐振动，即：
这里，平衡攻角αt=0°，振幅Am=1°，无量纲的减缩频率k=0.2。图 8给出了不同马赫数时俯仰运动绕平衡攻角的迟滞环。可以看到，随着马赫数降低，迟滞环由顺时针旋转演化至逆时针旋转，包围的面积先减小后增大，转化的临界马赫数约为2.2。
 |
| 图 8 不同马赫数时的俯仰运动迟滞圈Fig. 8 Hysteresis loop of pitching at different Mach numbers |
基于上述迟滞环，采用最小二乘法即可辨识出平衡攻角处的静导数、动导数和二阶动导数，辨识结果随马赫数的变化曲线见图 9。在所有马赫数下，静导数都小于0，即静稳定。随马赫数降低，动导数由小于0演化大于0，发生变号的临界马赫数约为2.2。二阶动导数量值很小，基本可忽略不计。
 |
| 图 9 俯仰静、动导数随马赫数降低的演化情况Fig. 9 Static and dynamic derivatives vs. Mach number |
由静导数、动导数和二阶动导数可计算出相平面上的特征参数p、q和Δ，进而得到分叉参数λ(Ma)随马赫数的变化曲线，其结果见图 10。图 11还给出了 随马赫数降低，平衡攻角附近的俯仰动态特性在p-q平面上的性态变化示意图。
 |
| 图 10 分叉参数λ随马赫数变化曲线Fig. 10 Bifurcation parameter λ vs. Mach number |
 |
| 图 11 平衡攻角在p-q平面上的性态变化示意图Fig. 11 Dynamic stabilization at balance angle of attack in p-q plane |
根据第二节的Hopf分叉理论分析，可以预测， 随 着马赫数降低，该飞船返回舱的俯仰运动将发生临 界Hopf分叉失稳，平衡攻角将由稳定的点吸引子演 化为极限环周期吸引子，分叉的临界马赫数约为2.2。这个结论正确与否即可验证第二节的Hopf分叉理论，下面通过数值模拟自激俯仰振荡的过程来进行验证。
自激振荡的起始攻角为3°或5°，计算马赫数为2.5、2.2和2.0，通过耦合求解方程(1)~(3)，数值模拟自激俯仰振荡的时间历程。图 12给出了不同马赫数 时俯仰角的时间历程曲线。图 13给出了不同马赫数 时俯仰角速度-俯仰角相图。可以看到，马赫数2.5时，返回舱从3°攻角释放后，振幅逐渐减小；马赫数2.2时，从3°和5°攻角释放后，振幅均不增加也不减小；而马赫数2.0时，俯仰振动的振幅不断增大，最后发展为极限环振荡。也即，随着马赫数降低，平衡攻 角由稳定的点吸引子状态演化为极限环周期吸引子 状态，从而验证了本文动态失稳的Hopf分叉理论。
 |
| 图 12 自激俯仰振荡的时间历程曲线Fig. 12 Time history of pitch-angle of self-induced pitching at different Mach numbers |
 |
| 图 13 自激俯仰振荡的俯仰角速度-俯仰角相图Fig. 13 Pitch velocity vs. pitch angle of self-induced pitching at different Mach numbers |
飞船返回舱是典型的大钝体再入飞行器，另一类再入飞行器为细长体，这里选用一种平头双锥细长体外形，图 14为外形示意图。外形呈轴对称，重心在体 轴上，因此，在不同马赫数下，均有0°平衡攻角。计算马赫数分别为7.5、7.0、6.0、5.0、4.0，计算方法、过程与上节相同，这里不再重复，只给典型的计算结果。
 |
| 图 14 平头双锥外形示意图Fig. 14 Sketch of the flat-nose winged double-cone body |
图 15给出了动导数随马赫数和俯仰角的变化曲线。图 16给出了Hopf分叉参数随马赫数降低的变 化规律， 变号的临界马赫数约为6.8。图 17给出平衡攻角 随马赫数降低 在p-q平面上的性态变化示意图。同理根据第2节的Hopf分叉理论分析，随马赫数降低，该平头双锥细长体外形的俯仰运动将发生临界Hopf分叉失稳，平衡攻角将由极限环周期吸引子演化为稳定的点吸引子，分叉的临界马赫数为Ma=6.8。
 |
| 图 15 不同马赫数时动导数随攻角的变化曲线Fig. 15 Dynamic derivatives with respect to angles of attack at different Mach numbers |
 |
| 图 16 分叉参数λ随马赫数变化曲线Fig. 16 λ with respect to Mach number |
 |
| 图 17 平衡攻角在p-q平面上的性态变化示意图Fig. 17 Dynamic stabilization at balance angle of attack in p-q plane |
同样对平头双锥外形，其俯仰运动随马赫数降低将发生临界Hopf分叉失稳的理论预测，也可用数值模拟自激俯仰振荡来验证。初始姿态角为1°和4°，俯仰角速度为1°/s。图 18为马赫数为6.0和7.0的时间历程曲 线，当Ma>6.8时，俯仰姿态角是发散的；当M<6.8时，俯仰姿态角是收敛的，进一步验证了关于飞行器再入时动态失稳的临界Hopf分叉理论。
 |
| 图 18 自激俯仰振荡的时间历程曲线Fig. 18 Time history of pitch-angle at different Mach numbers |
针对航天飞行器再入时的动态失稳问题，基于Etkin假设，建立了描述单自由度俯仰运动的平面自治动力系统数学模型，引入Hopf分叉理论，分析俯仰运动的动态特性随马赫数变化的演化规律，并进行了非定常数值模拟验证。主要结论如下：
(1) 发展了再入飞行器动态失稳的Hopf分叉理论分析方法，结合静、动导数计算，即可预测动态失稳现象、临界参数和失稳后的极限环运动形态。
(2) 航天飞行器再入俯仰动态Hopf分叉失稳现象，存在亚临界和超临界Hopf分叉两种形态。对于亚临界Hopf分叉，再入时随Ma降低，由点吸引子演化至周期吸引子，对飞行安全有影响；对于超临界Hopf分叉，再入时随Ma降低由周期吸引子演化至点吸引子，一般不影响飞行安全。
(3) 针对钝体和细长体两类典型再入飞行器，通过数值模拟自激俯仰振荡的演化历程，分别对应于亚临界和超临界Hopf分叉失稳。验证了所发展的动态失稳Hopf分叉理论分析方法。
| [1] | Cole D K, Robert D B, Mark S. Dynamic stability analysis of blunt body entry vehicles through the use of a time-lagged aftbody pitching moment[R]. AIAA 2013-0226. |
| [2] | Allen J H. Motion of a ballistic missile angular misaligned with the flight path upon entering the atmosphere and its effect upon aerodymic heating, aerodynamic loads, and miss distance[R]. NACA TN-4048, 1957. |
| [3] | Wehrend W R Jr. An experimental evaluation of aerodynamic damping moments of cones with different centers of rotation[R]. NASA TND-1768, 1963. |
| [4] | Cole D K, Robert D B, Ian G C. Survey of blunt body dynamic stability in supersonic flow[R]. AIAA 2012-4509. |
| [5] | James S G, Benjamin S K, Randolph P L, et al. Crew Exploration Vehicle (CEV) crew module shape selection analysis and CEV aeroscience project overview[R]. AIAA 2007-603. |
| [6] | Orlik-Rueckemann K J. Dynamic stability parameters[R]. AGARD CP-235, London:1978. |
| [7] | 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). |
| [8] | Yuan Xianxu. Numerical simulation for unsteady flows and research on dynamic characteristics of vehicle[D].[PhD.Thesis]. Mianyang:China Aerodynamics Research and Development Center, 2002. (in Chinese)袁先旭. 非定常流动数值模拟及飞行器动态特性分析研究[D].[博士学位论文]. 绵阳:中国空气动力研究与发展中心, 2002. |
| [9] | Etkin B. Dynamics of atmospheric flight[M]. Science Press, 1979: 147-158. |
| [10] | Zhao M X. Dynamic stability of manned reentry capsules[J]. Aerodynamic Experiments:Measurements and Control, 1995, 9(2):1-8. |
| [11] | Takashi Yoshinage, et al. Orbital re-entry experiment vehicle ground and flight dynamic test results comparison[J]. Journal of Spacecraft and Rockets, 1996, 33(5):635-642. doi:10.2514/3.26813 |
| [12] | Zhang Jinyan. Geometric theories and bifurcation problems of the ordinary differential equations[M]. Beijing:Peking University Press, 1981. (in Chinese)张锦炎. 常微分方程几何理论与分叉问题[M]. 北京:北京大学出版社, 1981. |
| [13] | Yuan Xianxu, Zhang Hanxin, Xie Yufei. Calculation of static and dynamic derivatives of pitching motions using CFD methods[J]. Acta Aerodynamica Sinica, 2005, 23(4): 458-463. (in Chinese)袁先旭, 张涵信, 谢昱飞. 基于CFD 方法的俯仰静、动导数数值计算[J]. 空气动力学学报, 2005, 23(4): 458-463. |
| [14] | Chen Qi, Chen Jianqiang, Yuan Xianxu, et al. Application of 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)陈琦, 陈坚强, 袁先旭, 等. 谐波平衡法在动导数快速预测中的应用研究[J]. 力学学报, 2014, 46(2):183-190. |
| [15] | Sun Tao, Gao Zhenghong, Huang Jiangtao. Identify of aircraft dynamic derivatives based on CFD technology and analysis of reduce frequency[J]. Flight Dynamics, 2011, 29(4): 15-18. (in Chinese)孙涛, 高正红, 黄江涛. 基于CFD 的动导数计算与减缩频率影响分析[J]. 飞行力学, 2011, 29(4):15-18. |
| [16] | Liu Wei, Shen Qing, Zhang Luming, 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)刘伟, 沈清, 张鲁民, 等. 钝倒锥体动导数数值工程模拟[J]. 国防科技大学学报, 1998, 20(1):5-8. |
| [17] | Thomas J P, Dowell E H, Hall K C. Nolinear inviscid aerodynamic effiects on transonic divergence, flutter and limit cycle oscillations[J]. AIAA Journal, 2002, 40(4):638-646. |
| [18] | Moor F G. Dynamic derivatives for missile configurations to Mach number three[J]. Journal of Spacecraft and Rockets, 1978, 15(2):65-66. |
| [19] | Fan Jingjing, Yan Chan, Li Yuejun. Computation of vehicle pitching static and dynamic derivatives at high angles of attack[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(10):1846-1850. (in Chinese)范晶晶, 阎超, 李跃军. 飞行器大迎角下俯仰静、动导数的数值计算[J]. 航空学报, 2009, 30(10):1846-1850. |
| [20] | Liu Wei, Zhao Haiyang, Yang Xiaoliang. Analysis of dynamic stability and research of passive control method for capsule[J]. Scientia Sinica(Phys., Mech. & Astron.), 2010, 40(9):1156-1164.(in Chinese)刘伟, 赵海洋, 杨小亮. 返回舱动态稳定性分析及被动控制方法研究[J]. 中国科学:物理学、力学、天文学, 2010, 40(9):1156-1164. |