再入飛行器俯仰動態失穩的分叉理論與計算分析

袁先旭，陳 琦，＊，何 琨，謝昱飛，張涵信

（1.中國空氣動力研究與發展中心計算空氣動力研究所，四川綿陽 621000；2.國家CFD實驗室，北京 100191）

（1），易得：

再入飛行器俯仰動態失穩的分叉理論與計算分析

袁先旭1，陳 琦1，＊，何 琨1，謝昱飛1，張涵信2

（1.中國空氣動力研究與發展中心計算空氣動力研究所，四川綿陽 621000；2.國家CFD實驗室，北京 100191）

采用非線性自治動力系統分叉理論，耦合求解非定常Navier-Stokes方程和俯仰運動方程，研究了再入飛行器單自由度俯仰運動失穩問題。研究表明，航天飛行器再入時，如果僅有一個配平攻角，隨馬赫數降低，其配平攻角處的俯仰動態失穩一般對應于Hopf分叉，并存在亞臨界Hopf分叉和超臨界Hopf分叉兩種失穩形態。作為驗證實例，數值模擬了飛船返回艙外形和平頭有翼雙錐外形的俯仰動態失穩現象。結果表明，返回艙再入時，隨馬赫數降低將發生超臨界Hopf分叉，俯仰運動由點吸引子演化為周期吸引子，臨界Hopf分叉點發生在馬赫數2.2處；而平頭再入體隨馬赫數降低，發生亞臨界Hopf分叉，俯仰運動則是由周期吸引子演化為點吸引子，馬赫數6.8為臨界Hopf分叉點。

Hopf分叉；動態失穩；再入飛行器；非定常數值模擬

0 引 言

航天再入飛行器出于防熱的需求，通常設計成鈍頭體構型。這類飛行器再入時，由于要經歷高超聲速、超聲速至亞跨聲速的全流域飛行，特別是在中-低超聲速階段，俯仰運動將會出現動態不穩定現象［1］，表現為俯仰振動的振幅逐漸增加，最終形成極限環振蕩。對再入飛行器動態失穩現象的研究始于20世紀50年代［2-3］，如今六十多年過去了，動態失穩的物理機制以及動態特性的準確預測仍然困擾著研究人員［4］。美國“阿波羅號”載人飛船在研制過程中，雖然采用9套模型分別在14座風洞中進行了亞、跨、超和高超聲速風洞試驗，試驗時間累計高達700多小時，但其動態特性問題仍沒有完全解決。最近基于“阿波羅號”飛船研制的“獵戶座”CEV，在投放試驗中由于多次出現動態失穩現象而失效，經過Langley的深入研究才得以解決，并于2014年12月5日飛行試驗成功。

已有的研究表明，小升阻比的載人飛船返回艙外形以球冠倒錐形為最優［5-7］。該類外形的飛船返回艙以較大的配平攻角再入，后體有大尺度的分離和再附，流場中存在復雜的波系結構（圖1）。鈍體氣流分離效應、后體氣流再附效應、船尾近尾渦流效應和動態遲滯效應等對飛船返回艙的靜、動穩定性都有影響。

圖1 返回艙再入時典型波系Fig.1 Flow characteristics of reentry capsule

因此，有研究者試圖通過有限度的改變返回艙后體外形，以改變后體分離拓撲形態，從而改善返回艙的穩定性，如“聯盟號”和“神舟號”載人飛船返回艙在光體倒錐上增加穩定翼［8］。當前已有的大量試驗與計算研究，一般只能給出動態失穩現象的唯象描述，缺乏理論分析［9-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），易得：

這里，下標0表示在平衡攻角處的定態俯仰力矩，故有（Cm）0＝0。非線性項是關于θ、的高階項，假定當時，非線性項更高階地趨于0，即非線性項滿足Perron定理［12］。將式（6）代入式（5），即得：

1.2 Hopf分叉分析

假設上述非線性系統式（8）的非線性項滿足Perron定理，根據動力學方法［8，12］，可以用非線性系統的一次近似系統式（9）來分析其平衡解的穩定性。

則一次近似系統式（9）的Jacobian矩陣的特征值為：。根據非線性動力學系統的有關定理，可知在p-q平面上各類平衡點（也稱奇點或臨界點）的區域如圖3。

圖3p-q平面上各類平衡點的區域圖Fig.3 Distribution of critical point inp-qplane

（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時，滿足：

于是，在λ＝λ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分叉失穩，出現極限環振蕩的臨界條件，即：

飛行器自由飛行時，Cμ＝0，一般情況下，二階動導數的量值遠小于1，于是有：

2 數值模擬結果與分析

針對飛船返回艙外形和平頭雙錐外形，對上述Hopf分叉理論進行驗證。對于每個外形，需要求解三類流場：首先，計算定態繞流流場，以獲取平衡攻角，并為后續的非定常計算提供初場；然后，數值模擬強迫俯仰振蕩過程，通過參數辨識，獲取平衡攻角處的靜、動導數［13-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。

基于上述遲滯環，采用最小二乘法即可辨識出平衡攻角處的靜導數、動導數和二階動導數，辨識結果隨馬赫數的變化曲線見圖9。在所有馬赫數下，靜導數都小于0，即靜穩定。隨馬赫數降低，動導數由小于0演化大于0，發生變號的臨界馬赫數約為2.2。二階動導數量值很小，基本可忽略不計。

圖8 不同馬赫數時的俯仰運動遲滯圈Fig.8 Hysteresis loop of pitching at different Mach numbers

圖9 俯仰靜、動導數隨馬赫數降低的演化情況Fig.9 Static and dynamic derivatives vs.Mach number

圖10 分叉參數λ隨馬赫數變化曲線Fig.10 Bifurcation parameterλvs.Mach number

圖11 平衡攻角在p-q平面上的性態變化示意圖Fig.11 Dynamic stabilization at balance angle of attack inp-qplane

由靜導數、動導數和二階動導數可計算出相平面上的特征參數p、q和Δ，進而得到分叉參數λ（Ma）隨馬赫數的變化曲線，其結果見圖10。圖11還給出了隨馬赫數降低，平衡攻角附近的俯仰動態特性在p-q平面上的性態變化示意圖。

圖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

根據第二節的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分叉理論。

2.2 平頭雙錐外形的數值模擬

飛船返回艙是典型的大鈍體再入飛行器，另一類再入飛行器為細長體，這里選用一種平頭雙錐細長體外形，圖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 不同馬赫數時動導數隨攻角的變化曲線Fig.15 Dynamic derivatives with respect to angles of attack at different Mach numbers

圖15給出了動導數隨馬赫數和俯仰角的變化曲線。圖16給出了Hopf分叉參數隨馬赫數降低的變化規律，變號的臨界馬赫數約為6.8。圖17給出平衡攻角隨馬赫數降低在p-q平面上的性態變化示意圖。同理根據第2節的Hopf分叉理論分析，隨馬赫數降低，該平頭雙錐細長體外形的俯仰運動將發生臨界Hopf分叉失穩，平衡攻角將由極限環周期吸引子演化為穩定的點吸引子，分叉的臨界馬赫數為Ma＝6.8。

同樣對平頭雙錐外形，其俯仰運動隨馬赫數降低將發生臨界Hopf分叉失穩的理論預測，也可用數值模擬自激俯仰振蕩來驗證。初始姿態角為1°和4°，俯仰角速度為1°/s。圖18為馬赫數為6.0和7.0的時間歷程曲線，當Ma＞6.8時，俯仰姿態角是發散的；當M＜6.8時，俯仰姿態角是收斂的，進一步驗證了關于飛行器再入時動態失穩的臨界Hopf分叉理論。

圖16 分叉參數λ隨馬赫數變化曲線Fig.16λwith respect to Mach number

圖17 平衡攻角在p-q平面上的性態變化示意圖Fig.17 Dynamic stabilization at balance angle of attack inp-qplane

圖18 自激俯仰振蕩的時間歷程曲線Fig.18 Time history of pitch-angle at different Mach numbers

3 結 論

針對航天飛行器再入時的動態失穩問題，基于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.

Dynamic destabilization analysis of the reentry vehicles using bifurcation theory and unsteady numerical simulation

Yuan Xianxu1，Chen Qi1，He Kun1，Xie Yufei1，Zhang Hanxin2
（1.Computational Aerodynamics Institute of China Aerodynamics Research and Development Center，Mianyang 621000，China；2.National Laboratory for CFD，Beijing 100191，China）

To study the problem about dynamic destabilization of the reentry vehicles，the bifurcation theory of nonlinear autonomous dynamic system is used together with the unsteady Navier-Stokes equations simulating the pitching process.The investigations shown that，the astronautic vehicles reentering into the atmosphere would generally occur a dynamic destabilization of pitching motion as their flight Mach number decreased，and the phenomenon of dynamic destabilization befallens at the point of the Hopf bifurcation，furthermore，there are two types of Hopf bifurcation named subcritical Hopf bifurcation and supercritical Hopf bifurcation.To validate the theory of the Hopf bifurcation，two reentry configurations are numerically studied.One is the Japanese capsule OREX，composed of a spherical cap facing forward and a reversed cone facing backward.With Mach number decreases in the reentry process for the capsule，the pitching motion would evolve from a point attractor to a periodic attractor，and the critical Mach number is about 2.2when the subcritical Hopf bifurcation occurs，the theory analysis and simulation result is in good agreement with the experiment and flight test results.The other configuration is a flat-nose winged double-cone body，with the decrease of Mach number，its pitching motion would evolve from a periodic attractor to a point attractor，and the critical Mach number at which the supercritical Hopf bifurcation occurs is about 6.8.

hopf bifurcation；dynamic destabilization；reentry vehicle；unsteady numerical simulation

V211.3

：Adoi：10.7638/kqdlxxb-2015.0022

0258-1825（2015）02-0162-08

2015-01-10；

：2015-03-24

國家自然科學基金（11372341，912162001，11172315）

袁先旭（1974-），男，研究員，研究方向：非定常空氣動力學、飛行器動態特性.E-mail：yuanxianxu＠mail.cardc.cn

陳琦＊，男，助理研究員，研究方向：計算空氣動力學.E-mail：chenqi＠mail.ustc.edu.cn

袁先旭，陳琦，何琨，等.再入飛行器俯仰動態失穩的分叉理論與計算分析［J］.空氣動力學學報，2015，33（2）：162-169.

10.7638/kqdlxxb-2015.0022 Yuan X X，Chen Q，He K，et al.Dynamic destabilization analysis of the reentry vehicles using bifurcation theory and unsteady numerical simulation［J］.Acta Aerodynamica Sinica，2015，33（2）：162-169.