類Starship 飛行器大迎角動態特性數值研究

呂凡熹，趙 飛，劉 瑜，杜若凡，石 泳，張宇佳

（中國空間技術研究院錢學森空間技術實驗室，北京 100094）

人類載人航天事業經歷了60 多年的發展歷程，形成了載人飛船和航天飛機為代表的兩類天地往返載人航天運輸器，此外便少有跨越性的突破。近年來，民營商用航天公司為該領域注入了新的活力，其中以SpaceX 公司研發的Starship 最為突出。有別于傳統傘降和機場水平降落的回收方式，Starship 利用操縱面的氣動力和火箭發動機矢量推力聯合控制完成再入減速和著陸過程，從而實現了可重復使用，顯著降低發射成本［1-2］。需要指出的是，再入過程的穩定飛行涉及很多關鍵技術，其中動態穩定性的水平將決定飛行任務的成?。?］。由于Starship 僅在著陸末端利用矢量推力進行飛行控制，要保證飛行器姿態受控，并且為乘員提供良好的生存環境，Starship 新穎氣動外形的動態穩定性是否良好至關重要。

為滿足快速再入減速需求，Starship 全程采用大迎角飛行，這必然帶來復雜的非定常流動問題。而非定常分離、旋渦和非定常運動激波及其與附面層的干擾使其流場呈現強烈的非線性效應，致使動導數辨識也十分困難［4］。同時，Starship 再入過程中要經歷寬速域、大空域的飛行條件，若通過風洞試驗或者在飛行試驗階段對其動態穩定性進行研究，將導致飛行器設計周期成倍增加，設計成本大幅增長。隨著近年來計算流體力學（Computational fluid dynamics，CFD）及大規模高性能計算技術的迅猛發展，國內外飛行器動態特性的研究廣泛應用了非定常CFD 方法。文獻［4-7］采用強迫振動法，數值模擬了彈道外形、尖錐、鈍錐、導彈以及返回艙等典型飛行器的動導數，結果表明飛行器在小迎角下動導數具備一定的線性特性。文獻［8］研究了彈道外形以及有翼導彈的大迎角下的動態特性，結果表明飛行器在大迎角下動導數具備非線性特性。上述主要還是針對簡單外形、單一速域下的研究工作，而復雜外形飛行器大迎角、寬速域非線性動態特性研究則相對較少。因此，應該充分利用非定常CFD 方法對Starship 這類創新構型飛行器的動態特性進行詳細分析，從而充分降低設計成本，盡可能減少設計風險。

非定常CFD 方法求解動導數充分評估流場的非線性特性，且能辨識更全面的動導數，主要方法有：強迫振動法［9］、自由振動法［10］、準定常法［11］和諧波平衡法［12］。強迫振動法可辨識的動導數種類較其他方法更為全面，對交叉導數、交叉耦合導數、加速度導數均能較好地開展計算評估；自由振動法的模擬過程與真實飛行情況狀態吻合度更高，但對辨識方法要求高，計算工作量大，僅適用于平衡狀態下的動導數計算，且動導數計算種類有限；準定常法計算精度高且速度快，但只能用于旋轉導數的評估；諧波平衡法計算效率高，但是對非定常流動的計算精度遠低于上述時域計算方法。

因此，本文采用基于小振幅強迫振動法的非定常CFD 方法對類Starship 飛行器寬速域、大空域、大迎角再入減速過程的動態特性進行系統評估。本文采用一種基于分區彈簧近似法動網格技術的三維非定常N-S 方程有限體積法進行了飛行器動態特性的仿真評估，經HBS 標準模型驗證，計算精度滿足需求。應用該方法深入評估了類Starship飛行器再入過程中的動態特性，包含不同馬赫數、重心位置、操縱面偏轉、迎角以及減縮頻率對飛行器俯仰組合導數的影響。

1 計算方法

1.1 控制方程

本文采用有限體積法求解非定?？蓧嚎sN-S方程，其三維積分形式為

本文有限體積法在空間方上采用Piecewise Linear 方法進行變量重構，對流通量的計算選擇二階精度Roe 格式以及Venkatakrishnan 限制器，湍流模型選取SST 兩方程模型。時間方向上采用基于GMRES 算法的全隱式迭代求解。

1.2 動網格方法

強迫振動方法需要使用動網格技術，而大量相關研究中強迫振動方法的數值模擬常采用剛性網格法［13］、滑移網格法［14］以及重疊網格法［15］，這3 種方法能有效保持網格品質，但剛性網格法在遠場邊界、滑移網格法在滑移面和重疊網格法在重疊區，均會因插值誤差而導致計算誤差，并將隨時間累積不斷放大。而強迫振動法振動幅度較小，彈簧近似法既能滿足動網格的需求又可避免插值誤差的產生。文獻［16］中提出了線性彈簧近似法并將其應用于機翼顫振問題，但是線性彈簧模型易導致負體積現象的出現。為此，文獻［17］中增加了扭轉彈簧模型進行完善。此后彈簧近似法不斷改進并開展了大量的應用。

但是，傳統的彈簧近似法無論改進與否，對復雜問題的處理效率和穩定性仍然有限［18］。因此，本文采用一種基于分區思想的彈簧近似法動網格技術，相比傳統的彈簧近似法，該方法計算量小且可保持較高網格質量，進而在避免插值誤差的基礎上保證數值計算的穩定性。如圖1 所示，計算網格被分成剛性網格區和彈性網格區。剛性網格區是包含整個飛行器及其周邊邊界層網格的球形區域，區域幾何中心為飛行器質心，也是強迫振動中網格的轉動中心。該區域要求足夠大確保該剛性網格區外邊界附近的網格質量較高。剛性網格區之外其他計算區域均為彈性網格區。在強迫振動過程中，剛性網格區與飛行器固連，其網格整體隨飛行器做強迫振動，從而保證壁面附近的網格質量不因邊界運動而變差。彈性網格區則以剛性網格區外邊界作為運動邊界，該區域網格中每條邊為彈簧模型，各個網格點在周圍網格點的作用力下達到受力平衡，每個網格點的靜態平衡方程為

圖1 基于分區彈簧近似法的計算網格Fig.1 Computational grid based on partitioned spring analogy method

圖2 給出了采用分區彈簧近似法從平衡位置振動到最大振幅過程中的網格變形過程，為方便展示，圖2（b）最大振幅達到了45°。由圖2 可見，在網格變形過程中剛性網格區網格沒有形變，有效保證了區內網格質量。彈性網格區由于本身網格質量較好，雖然在變形過程中會有所下降，但能保證非定常動態計算的魯棒性。以圖2 中展示的NACA0012 翼型的強迫振動為例，振幅達45°時不同工況的計算仍能保持時間迭代正常推進，雙時間步內迭代能有效收斂，流場結構合理可信。但振幅達60°時，由于剛性區大幅度地扭轉，彈性區網格質量已不能確保非定常計算的正常進行。但是，對于動導數計算所采用的小振幅強迫運動法，求解過程中振幅較小，彈性網格區的網格變形后質量不會出現顯著下降，更不會出現負體積等極端情況，本文方法能充分滿足強迫振動法的網格變形需求。

圖2 網格變形過程Fig.2 Mesh deformation process

1.3 動導數辨識

本文采用小振幅強迫振動法辨識飛行器動導數。以俯仰動導數為例，在體軸坐標系下，給定繞質心的強迫運動形式（無量綱形式）為

2 算例驗證

為了驗證數值方法的準確性，選取高超聲速彈道外形（Hyper ballistic shape，HBS）來作為驗證算例。HBS 由1 個半球柱和2 段擴張尾裙組成，該外形有較為精確的試驗結果［19］，常被用作驗證標模。

驗證算例的流動條件為：Ma=6.85，基于底部直徑的Red=0.72×106，縱向質心位置Xcg/L=0.72。數值計算時取流動條件為：壓力461.51 Pa，密度0.005 545 kg/m3，溫度290 K，計算時參考長度取頭部直徑d=1 m，參考面積為0.785 498 m2，幾何外形如圖3 所示。迎角分別取α0=0°、3.068 9°、6.443 7°、9.054 5°、13.388 1°、15°、20°，振幅取αm=1°，減縮頻率k=0.01。

圖3 HBS 幾何外形Fig.3 HBS geometry

圖4 給出了本文計算結果和試驗數據以及文獻［15］中的CFD 計算結果的對比。結果可見，所有計算工況動導數均為負值。當迎角小于13°俯仰動導數變化相對較小，當迎角大于13°俯仰動導數迅速減小，俯仰動穩定性增加。本文結果與試驗數據吻合度較好，誤差小于9%，同時對于缺少試驗數據的20°迎角工況，本文結果與文獻CFD 結果也誤差也較小，證明本文方法對于求解飛行器動導數的有效性。

圖4 HBS 動導數計算結果對比Fig.4 Comparison of HBS dynamic derivative calculation results

3 類Starship 動態特性分析

3.1 幾何外形

本文研究的類Starship飛行器根據SpaceX 官網公布的Starship尺寸信息進行三維建模［20］，飛行器總長50 m，直徑d=9 m，如圖5 所示。前翼展長15 m，后翼展長18 m，前后翼舵偏范圍均為向上偏轉0°～60°?？紤]大迎角再入過程中，前后翼在避免承受過大的力、熱載荷的同時，均需保留足夠舵效。因此，除下文分析前后翼偏轉對動態特性影響的工況，其他計算工況中使用的模型前后翼均向上偏轉30°。

圖5 類Starship 飛行器三維模型Fig.5 Starship-like spacecraft 3-D model

3.2 計算網格

本文計算網格采用非結構混合網格，總網格量1 200 萬個，物面網格和邊界層網格如圖6 所示。如無特別說明，下文計算工況質心位于縱向質心位置Xcg/L=0.6，強迫振動的減縮頻率k=ωd/2u∞取0.05，振幅取αm=1°。飛行器采用高空高速大迎角的再入方式，故迎角α0取70°，側滑角β取0°。

圖6 表面網格及邊界層網格Fig.6 Surface and boundary layer meshes

3.3 馬赫數

類Starship 飛行器再入過程是一個高度和馬赫數逐漸減小的過程，參考有翼再入飛行器典型再入彈道，選取不同高度下，70°迎角、0°側滑、前后翼操縱面偏角均為30°，Ma范圍為0.3～10 的計算工況，如表1 所示。

表1 計算工況Table 1 Computation cases

圖7 給出了Ma=0.3、0.9、1.5、5.0 時，1 個周期初始時刻的計算流場，圖中展示了飛行器表面和流場對稱面的壓力分布及馬赫數等值線的分布情況，迎風面整體高壓、背風面低壓，背風面尾部存在大分離區，從高速結果中可以清晰識別流場中的波系，仿真結果的流動結構合理可信。圖8（a）給出了Ma=5.0 工況俯仰力矩系數的遲滯環曲線。遲滯環形狀飽滿展示了非定常運動時物體表面氣動力的作用和物體運動之間存在較強的時間延遲和滯后效應，并在計算開始后1 個周期內即進入遲滯環。從時間歷程上看非定常氣動力矩系數隨時間步的增長很快按周期規律達到諧振狀態，如圖8（b）所示。

圖7 1 個周期初始時刻表面壓力和流場馬赫數云圖Fig.7 Surface pressure and flow field Mach number contour at initial moment in one period

圖8 俯仰力矩系數的遲滯環和時間歷程(Ma=5.0)Fig.8 Hysteresis loop and time history of pitching moment coefficient(Ma=5.0)

圖9 給出了計算得到的飛行器俯仰組合導數隨飛行馬赫數的變化，可見再入過程中飛行器從高超聲速逐漸減速至Ma=2.0 時，其俯仰組合導數均為負值且變化較小，非定常氣動力對飛行器俯仰振動運動起阻尼作用。當Ma＜2.0 飛行器減速進入跨聲速區，俯仰組合導數開始增大，在Ma=0.8 附近由負值變為正值，Ma=0.6 附近達到峰值。飛行器進一步減速進入亞聲速區域，俯仰組合導數由峰值開始減小但仍為正值。

圖9 大迎角俯仰動導數隨不同馬赫數的變化Fig.9 Variation of pitch dynamic derivative with different Mach number at high angle of attack

3.4 重心位置

在飛行器飛行過程中，隨著推進劑的消耗，重心會隨之變化。為了研究重心位置對飛行器俯仰動導數的影響，選取來流馬赫數Ma=0.3 和5.0，迎角α0=70°，縱 向 質 心 相 對 位 置Xcg/L=0.3、0.45、0.6、0.75 進行俯仰組合導數辨識，如圖10所示。

圖10 不同質心位置Fig.10 Different centroid positions

圖11 給出了Ma=0.3 工況不同質心位置時俯仰力矩系數的遲滯環曲線和時間歷程，力矩系數經過了大約1 個周期后才呈現穩定的周期振蕩。同時強迫振動過程中初始迎角位置處的力矩系數大幅度偏離定常解，表明低速流動下氣動力的非定常效應對氣動特性影響極為顯著。

圖11 不同質心位置俯仰力矩系數的遲滯環和時間歷程（Ma=0.3）Fig.11 Hysteresis loop and time history of pitching moment coefficient at different centroid positions（Ma=0.3）

圖12 給出了Ma=5.0 工況不同質心位置時俯仰力矩系數的結果，其力矩系數經歷不到半個周期即發展成穩定的周期振蕩，初始迎角位置處的力矩系數偏離定常解的幅度相比Ma=0.3 工況結果要小很多，高速流動下非定常效應的影響要小于低速流動。

圖12 不同質心位置俯仰力矩系數的遲滯環和時間歷程（Ma=5.0）Fig.12 Hysteresis loop and time history of pitching moment coefficient at different centroid positions（Ma=5.0）

圖13 給出了不同質心位置對俯仰組合導數的計算結果。對于Ma=0.3，隨著質心位置后移，俯仰動導數逐漸增大，當質心位于Xcg/L=0.53 之后，俯仰動導數變為正值，呈動態不穩定狀態。對于Ma=5.0，隨著質心位置后移，俯仰動導數先增大后減小，最大值出現在質心位于Xcg/L=0.6 左右，俯仰動導數均為負值，呈動態穩定狀態。

圖13 大迎角下質心位置對俯仰動導數的影響Fig.13 Effect of centroid positions on pitch dynamic derivative at high angle of attack

3.5 操縱面偏轉角

為評估飛行器高速再入過程中前后翼操縱面偏轉對俯仰動導數的影響，選取馬赫數Ma=5，迎角α0=70°，選取前翼、后翼、前后翼同步3 種模式，按0°、30°、60°的舵偏狀態進行評估，動導數計算結果如圖14 所示?？梢?，隨著前/后翼偏轉角逐漸增大，俯仰動導數呈單調上升的趨勢，前/后翼舵偏均將導致飛行器俯仰動穩定性減弱。對比前/后翼偏轉的兩種模式，前翼偏轉對俯仰動導數的影響相對較小，而后翼偏轉則影響更為顯著。同時，前后翼同步偏轉對俯仰動導數的影響與后翼偏轉相當，并沒有呈現出前后翼偏轉作用累加的現象。主要原因是前翼投影面積明顯小于后翼，因此前翼對動導數的影響較后翼會小很多。同時大迎角高超聲速流場非線性程度大，激波和分離同時存在，動導數不會隨前后翼大幅度偏轉呈現簡單的線性變化特征。

圖14 大迎角下俯仰動導數隨操縱面偏轉的變化Fig.14 Variations of pitch dynamic derivative with deflection angle of control surface at high angle of attack

3.6 迎 角

為了研究不同迎角對該類飛行器俯仰動導數的影響，根據不同速度下再入過程可能的迎角范圍，評估了Ma=0.3、迎角α0=0°～180°和Ma=5.0、迎角α0=0°～90°下的俯仰動導數，其中計算工況的縱向質心位置位于Xcg/L=0.6。圖15（a）為Ma=0.3 時，不同迎角俯仰動導數的計算結果。當迎角α0=0°時，俯仰動導數為負值，后隨迎角增大，俯仰動導數先后在65°、125°、170°迎角附近經歷了3 次變號，動導數為負值的迎角范圍 分 布 在0°～65°和125°～170°之 間。圖15（b）為Ma=5.0 時，俯仰動導數隨迎角的變化曲線。當迎角α0=0°時，俯仰動導數為負值，后隨迎角增大，俯仰動導數逐漸減小，在70°迎角附近達到最小值。對于Ma=5.0 所評估的不同迎角的工況，俯仰動導數均為負值，氣動力對飛行器俯仰動起阻尼效果。

圖15 迎角對俯仰動導數的影響Fig.15 Effect of angle of attack on pitch dynamic derivative

3.7 減縮頻率

本節評估了不同減縮頻率K對大迎角再入過程中俯仰動導數的影響。針對Ma=0.3、5.0，迎角70°，縱向質心位置X/L=0.6，減縮頻率分別選擇0.01、0.05、0.25、0.50、1.00 開展動導數計算。如圖16 所示，對Ma=0.3 的計算工況，動導數在高頻不敏感，但在低頻時，動導數隨減縮頻率的降低迅速增大，飛行器從負值轉變為正值。對Ma=5.0 的計算工況，動導數對減縮頻率的變化不敏感，結果符合對高速飛行器動態特性的經驗認識。需要指出的是，減縮頻率是影響動導數計算結果的重要因素，適當的減縮頻率才能得到正確的計算結果，而減縮頻率應由典型飛行狀態的氣動特性和質量特性估算給出。本文目標是研究類Starship 外形飛行器的動態特性，而通常飛行器超聲速下的動態特性動導數隨頻率變化不大，圖16 中，Ma=5.0時的仿真結果也驗證了該經驗認識。因此在質量特性未知的前提下，上文在減縮頻率的常用取值范圍內采用0.05 開展了動導數評估，這也是大量研究任務和文獻資料中所采取的相對適中的減縮頻率。

圖16 減縮頻率對動導數預測結果的影響Fig.16 Effect of reduced frequency on dynamic derivative prediction result

圖17 為不同減縮頻率下俯仰運動對應力矩系數的遲滯環曲線的計算結果。當強迫振動過程中力矩系數達到諧振解時，振動頻率越高，其遲滯環越飽滿，其力矩系數的峰值也越大。同時在相同迎角下，由于運動方向的不同，對應有2 個不同的力矩系數。減縮頻率的越小，遲滯環越逼近定常解。而減縮頻率的越大，2 個力矩系數的差值也增大，呈現了流場遲滯效應對氣動力的影響。對比Ma=0.3 和Ma=5.0 的遲滯環結果，Ma=0.3 遲滯環的長軸會隨減縮頻率的增大而偏轉。減縮頻率的物理含義可以理解為在自由來流運動參考長度的時間內物體做簡諧振動位置的相位角變化率，因此相同減縮頻率下低速結果遲滯環長軸的更強偏轉也體現了低速擾動波傳播速度更慢，偏離穩態解更顯著的流動特征。

圖17 減縮頻率對遲滯環的影響Fig.17 Effect of reduced frequency on hysteresis loop

4 結 論

本文采用基于分區彈簧近似法動網格技術的強迫振動方法計算研究了類Starship 飛行器在典型大迎角再入飛行狀態下（相對質心位置0.6，前后翼向上偏轉30°，迎角70°，減縮頻率0.05）不同飛行馬赫數、重心位置、操縱面偏轉角、迎角、減縮頻率的俯仰動態特性變化規律，結果表明：

（1）采用典型再入飛行狀態時，飛行器超聲速的俯仰動導數均為負值，而Ma＜0.9 時的結果則轉變為正值。

（2）在相對質心位置0.3～0.75 的變化范圍內，隨著質心后移，Ma=0.3 時俯仰動導數由負值變為正值，而Ma=5.0 時俯仰動導數則均為負值。

（3）對于Ma=5.0 工況，隨操縱面上偏，俯仰動導數逐漸增大，阻尼效應隨之減小，且后翼偏轉的影響相比于前翼更顯著。

（4）Ma=5.0、迎角0°～90°范圍內俯仰動導數均為負值，Ma=0.3 時俯仰動導數為負值的迎角范圍分布在0°～65°和125°～170°之間。

（5）Ma=0.3 時俯仰動導數的大小與減縮頻率密切相關，Ma=5.0 時俯仰動導數則對減縮頻率不敏感。

本文研究系統掌握了類Starship 飛行器典型大迎角再入飛行狀態下的俯仰動態特性，動導數計算方法適合應用于各類寬速域大空域復雜外形飛行器的動態特性數值仿真。