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

（中國(guó)空間技術(shù)研究院錢學(xué)森空間技術(shù)實(shí)驗(yàn)室，北京 100094）

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

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

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

因此，本文采用基于小振幅強(qiáng)迫振動(dòng)法的非定常CFD 方法對(duì)類Starship 飛行器寬速域、大空域、大迎角再入減速過程的動(dòng)態(tài)特性進(jìn)行系統(tǒng)評(píng)估。本文采用一種基于分區(qū)彈簧近似法動(dòng)網(wǎng)格技術(shù)的三維非定常N-S 方程有限體積法進(jìn)行了飛行器動(dòng)態(tài)特性的仿真評(píng)估，經(jīng)HBS 標(biāo)準(zhǔn)模型驗(yàn)證，計(jì)算精度滿足需求。應(yīng)用該方法深入評(píng)估了類Starship飛行器再入過程中的動(dòng)態(tài)特性，包含不同馬赫數(shù)、重心位置、操縱面偏轉(zhuǎn)、迎角以及減縮頻率對(duì)飛行器俯仰組合導(dǎo)數(shù)的影響。

1 計(jì)算方法

1.1 控制方程

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

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

1.2 動(dòng)網(wǎng)格方法

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

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

圖1 基于分區(qū)彈簧近似法的計(jì)算網(wǎng)格Fig.1 Computational grid based on partitioned spring analogy method

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

圖2 網(wǎng)格變形過程Fig.2 Mesh deformation process

1.3 動(dòng)導(dǎo)數(shù)辨識(shí)

本文采用小振幅強(qiáng)迫振動(dòng)法辨識(shí)飛行器動(dòng)導(dǎo)數(shù)。以俯仰動(dòng)導(dǎo)數(shù)為例，在體軸坐標(biāo)系下，給定繞質(zhì)心的強(qiáng)迫運(yùn)動(dòng)形式（無量綱形式）為

2 算例驗(yàn)證

為了驗(yàn)證數(shù)值方法的準(zhǔn)確性，選取高超聲速彈道外形（Hyper ballistic shape，HBS）來作為驗(yàn)證算例。HBS 由1 個(gè)半球柱和2 段擴(kuò)張尾裙組成，該外形有較為精確的試驗(yàn)結(jié)果［19］，常被用作驗(yàn)證標(biāo)模。

驗(yàn)證算例的流動(dòng)條件為：Ma=6.85，基于底部直徑的Red=0.72×106，縱向質(zhì)心位置Xcg/L=0.72。數(shù)值計(jì)算時(shí)取流動(dòng)條件為：壓力461.51 Pa，密度0.005 545 kg/m3，溫度290 K，計(jì)算時(shí)參考長(zhǎng)度取頭部直徑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 給出了本文計(jì)算結(jié)果和試驗(yàn)數(shù)據(jù)以及文獻(xiàn)［15］中的CFD 計(jì)算結(jié)果的對(duì)比。結(jié)果可見，所有計(jì)算工況動(dòng)導(dǎo)數(shù)均為負(fù)值。當(dāng)迎角小于13°俯仰動(dòng)導(dǎo)數(shù)變化相對(duì)較小，當(dāng)迎角大于13°俯仰動(dòng)導(dǎo)數(shù)迅速減小，俯仰動(dòng)穩(wěn)定性增加。本文結(jié)果與試驗(yàn)數(shù)據(jù)吻合度較好，誤差小于9%，同時(shí)對(duì)于缺少試驗(yàn)數(shù)據(jù)的20°迎角工況，本文結(jié)果與文獻(xiàn)CFD 結(jié)果也誤差也較小，證明本文方法對(duì)于求解飛行器動(dòng)導(dǎo)數(shù)的有效性。

圖4 HBS 動(dòng)導(dǎo)數(shù)計(jì)算結(jié)果對(duì)比Fig.4 Comparison of HBS dynamic derivative calculation results

3 類Starship 動(dòng)態(tài)特性分析

3.1 幾何外形

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

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

3.2 計(jì)算網(wǎng)格

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

圖6 表面網(wǎng)格及邊界層網(wǎng)格Fig.6 Surface and boundary layer meshes

3.3 馬赫數(shù)

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

表1 計(jì)算工況Table 1 Computation cases

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

圖7 1 個(gè)周期初始時(shí)刻表面壓力和流場(chǎng)馬赫數(shù)云圖Fig.7 Surface pressure and flow field Mach number contour at initial moment in one period

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

圖9 給出了計(jì)算得到的飛行器俯仰組合導(dǎo)數(shù)隨飛行馬赫數(shù)的變化，可見再入過程中飛行器從高超聲速逐漸減速至Ma=2.0 時(shí)，其俯仰組合導(dǎo)數(shù)均為負(fù)值且變化較小，非定常氣動(dòng)力對(duì)飛行器俯仰振動(dòng)運(yùn)動(dòng)起阻尼作用。當(dāng)Ma＜2.0 飛行器減速進(jìn)入跨聲速區(qū)，俯仰組合導(dǎo)數(shù)開始增大，在Ma=0.8 附近由負(fù)值變?yōu)檎?，Ma=0.6 附近達(dá)到峰值。飛行器進(jìn)一步減速進(jìn)入亞聲速區(qū)域，俯仰組合導(dǎo)數(shù)由峰值開始減小但仍為正值。

圖9 大迎角俯仰動(dòng)導(dǎo)數(shù)隨不同馬赫數(shù)的變化Fig.9 Variation of pitch dynamic derivative with different Mach number at high angle of attack

3.4 重心位置

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

圖10 不同質(zhì)心位置Fig.10 Different centroid positions

圖11 給出了Ma=0.3 工況不同質(zhì)心位置時(shí)俯仰力矩系數(shù)的遲滯環(huán)曲線和時(shí)間歷程，力矩系數(shù)經(jīng)過了大約1 個(gè)周期后才呈現(xiàn)穩(wěn)定的周期振蕩。同時(shí)強(qiáng)迫振動(dòng)過程中初始迎角位置處的力矩系數(shù)大幅度偏離定常解，表明低速流動(dòng)下氣動(dòng)力的非定常效應(yīng)對(duì)氣動(dòng)特性影響極為顯著。

圖11 不同質(zhì)心位置俯仰力矩系數(shù)的遲滯環(huán)和時(shí)間歷程（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 工況不同質(zhì)心位置時(shí)俯仰力矩系數(shù)的結(jié)果，其力矩系數(shù)經(jīng)歷不到半個(gè)周期即發(fā)展成穩(wěn)定的周期振蕩，初始迎角位置處的力矩系數(shù)偏離定常解的幅度相比Ma=0.3 工況結(jié)果要小很多，高速流動(dòng)下非定常效應(yīng)的影響要小于低速流動(dòng)。

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

圖13 給出了不同質(zhì)心位置對(duì)俯仰組合導(dǎo)數(shù)的計(jì)算結(jié)果。對(duì)于Ma=0.3，隨著質(zhì)心位置后移，俯仰動(dòng)導(dǎo)數(shù)逐漸增大，當(dāng)質(zhì)心位于Xcg/L=0.53 之后，俯仰動(dòng)導(dǎo)數(shù)變?yōu)檎担蕜?dòng)態(tài)不穩(wěn)定狀態(tài)。對(duì)于Ma=5.0，隨著質(zhì)心位置后移，俯仰動(dòng)導(dǎo)數(shù)先增大后減小，最大值出現(xiàn)在質(zhì)心位于Xcg/L=0.6 左右，俯仰動(dòng)導(dǎo)數(shù)均為負(fù)值，呈動(dòng)態(tài)穩(wěn)定狀態(tài)。

圖13 大迎角下質(zhì)心位置對(duì)俯仰動(dòng)導(dǎo)數(shù)的影響Fig.13 Effect of centroid positions on pitch dynamic derivative at high angle of attack

3.5 操縱面偏轉(zhuǎn)角

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

圖14 大迎角下俯仰動(dòng)導(dǎo)數(shù)隨操縱面偏轉(zhuǎn)的變化Fig.14 Variations of pitch dynamic derivative with deflection angle of control surface at high angle of attack

3.6 迎 角

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

圖15 迎角對(duì)俯仰動(dòng)導(dǎo)數(shù)的影響Fig.15 Effect of angle of attack on pitch dynamic derivative

3.7 減縮頻率

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

圖16 減縮頻率對(duì)動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果的影響Fig.16 Effect of reduced frequency on dynamic derivative prediction result

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

圖17 減縮頻率對(duì)遲滯環(huán)的影響Fig.17 Effect of reduced frequency on hysteresis loop

4 結(jié) 論

本文采用基于分區(qū)彈簧近似法動(dòng)網(wǎng)格技術(shù)的強(qiáng)迫振動(dòng)方法計(jì)算研究了類Starship 飛行器在典型大迎角再入飛行狀態(tài)下（相對(duì)質(zhì)心位置0.6，前后翼向上偏轉(zhuǎn)30°，迎角70°，減縮頻率0.05）不同飛行馬赫數(shù)、重心位置、操縱面偏轉(zhuǎn)角、迎角、減縮頻率的俯仰動(dòng)態(tài)特性變化規(guī)律，結(jié)果表明：

（1）采用典型再入飛行狀態(tài)時(shí)，飛行器超聲速的俯仰動(dòng)導(dǎo)數(shù)均為負(fù)值，而Ma＜0.9 時(shí)的結(jié)果則轉(zhuǎn)變?yōu)檎怠?/p>

（2）在相對(duì)質(zhì)心位置0.3～0.75 的變化范圍內(nèi)，隨著質(zhì)心后移，Ma=0.3 時(shí)俯仰動(dòng)導(dǎo)數(shù)由負(fù)值變?yōu)檎?，而Ma=5.0 時(shí)俯仰動(dòng)導(dǎo)數(shù)則均為負(fù)值。

（3）對(duì)于Ma=5.0 工況，隨操縱面上偏，俯仰動(dòng)導(dǎo)數(shù)逐漸增大，阻尼效應(yīng)隨之減小，且后翼偏轉(zhuǎn)的影響相比于前翼更顯著。

（4）Ma=5.0、迎角0°～90°范圍內(nèi)俯仰動(dòng)導(dǎo)數(shù)均為負(fù)值，Ma=0.3 時(shí)俯仰動(dòng)導(dǎo)數(shù)為負(fù)值的迎角范圍分布在0°～65°和125°～170°之間。

（5）Ma=0.3 時(shí)俯仰動(dòng)導(dǎo)數(shù)的大小與減縮頻率密切相關(guān)，Ma=5.0 時(shí)俯仰動(dòng)導(dǎo)數(shù)則對(duì)減縮頻率不敏感。

本文研究系統(tǒng)掌握了類Starship 飛行器典型大迎角再入飛行狀態(tài)下的俯仰動(dòng)態(tài)特性，動(dòng)導(dǎo)數(shù)計(jì)算方法適合應(yīng)用于各類寬速域大空域復(fù)雜外形飛行器的動(dòng)態(tài)特性數(shù)值仿真。