充氣式再入減速器動態氣動載荷與結構特性研究

吳 杰，張 章，侯安平，王立武，王 潔，曹 旭

(1. 北京航空航天大學能源與動力工程學院，北京 100191；2. 北京空間機電研究所中國空間技術研究院航天器無損著陸技術核心專業實驗室，北京 100094)

0 引 言

隨著航天技術的不斷發展，空天往返運輸活動也更加頻繁，然而傳統剛性返回器的裝載能力容易受到結構尺寸的限制，導致有效載荷的局限性越來越突出[1-3]。針對這一缺陷，充氣式返回技術應運而生，其在發射過程中折疊以占有盡可能小的體積，返回時通過充氣使氣囊展開并形成氣動外形[4-5]。然而區別于剛性結構，柔性結構型面的氣動變形問題更加突出：高超聲速氣動力會導致薄膜結構的大變形，甚至產生結構斷裂與解體[6-7]；氣動熱載荷會使內充壓氣體劇烈膨脹并超過氣囊縫合強度，從而造成柔性結構破壞。因此研究再入過程中氣動載荷與結構特性的變化規律，校核危險點并提出優化方案，在充氣式再入減速器的設計過程中愈加重要。

針對充氣式再入減速器在氣動力與氣動熱作用下的結構特性，國內外學者已開展了一定的研究工作。現有研究中，在氣動力的求解上主要通過CFD計算個別高度工況下的氣動壓力，并采取單向壓力加載的方式求解結構應力[8-11]。然而，由于其通過插值來獲取其余高度的氣動結果，無法獲得更為連續的氣動載荷與結構特性變化歷程，并且單向耦合無法計入結構變形對于下一時刻流場分布的影響。另一方面，現有研究在求解氣動熱上主要依靠工程算法或CFD得到表面熱流分布，再通過熱傳導方程得到各功能層溫度分布[12-15]，然而其未能考慮內充壓氣體的熱膨脹這一重要因素對于熱應力的影響。

針對現有研究方法的不足，本文提出了將動態飛行參數加載至CFD邊界條件的技術手段，并且在CFD中通過數次迭代求解更為精確的阻力系數，實現飛行動力學與空氣動力學之間的雙向耦合。同時基于ANSYS建立了考慮內充壓氣體熱效應的流固雙向耦合模型，既考慮了結構變形對下一時刻流場的影響，又考慮了內充壓氣體溫度變化對熱應力的影響。利用此模型研究了氣動力與氣動熱分別作用下的結構變化規律，并通過改變半錐角和氣囊數目，研究結構參數對氣動載荷與結構特性的影響。

1 數值方法及驗證

1.1 分析流程

為對比氣動力與氣動熱對結構特性的影響，本文將兩者進行分層求解，如圖1所示。氣動力的分析中，將動態飛行參數作為CFD的邊界條件，并結合流固雙向耦合模型進行求解。其中，通過CFD求出的阻力系數能夠對面元積分法得到的初始阻力系數進行反饋修正。

在氣動熱分析中，將CFD計算得到的氣動熱作為邊界條件依次加載至熱力學及動力學模型中，并且將薄膜熱變形以及內充壓氣體的熱膨脹都計入熱應力的影響因素中。

1.2 幾何模型

充氣式再入減速器的幾何模型選取美國IRVE飛行器[10]，結構如圖2所示。七層氣囊獨立充氣以提供剛性，氣囊之間通過隔層連接以增加結構強度，蒙皮包裹在氣囊外部以形成穩定的氣動外形。柔性薄膜材料選擇凱夫拉(Kevlar)纖維膜，該材料密度為1440 kg/m3，楊氏模量131 GPa，泊松比0.35。

圖2 IRVE結構示意圖Fig.2 Schematic diagram of IRVE system

圖3展示了充氣式再入減速器的結構有限元模型，采用shell 181四節點面單元，網格數10萬。圖4是流體域網格模型的剖面圖，采用混合網格進行劃分，總網格800萬。

圖3 有限元模型Fig.3 Finite element model

圖4 CFD模型Fig.4 CFD model

充氣式再入減速器通過熱防護系統(Thermal protection system, TPS)抵御巨大的熱載荷，從外到內分為防熱層、絕熱層和氣密層，分別用以承受熱流、隔絕熱流和防止氣體泄漏，如圖5所示。本文采用solid 45單元建立各功能層的熱力學模型。

圖5 TPS結構示意圖Fig.5 Diagram of the thermal protection system

1.3 氣動力與氣動熱計算方法

流體域和固體域通過耦合交界面將位移及壓力依據協調條件進行轉換。采用基于有限體積法的雷諾平均N-S方程求解氣動力，并采用SST湍流模型。在求解氣動熱時需考慮邊界層內高焓氣流與結構之間的對流換熱，以及壁面產生的輻射。設σ為黑體輻射系數，給定結構發射率ε為0.89，其熱輻射方程為：

(1)

利用文獻[16]中的CFD結果進行了驗證，如圖6所示。駐點動壓和熱流密度的最大偏差均發生在各自峰值處，分別為5.8%和5.1%，均較為合理。

圖6 氣動力及氣動熱方法驗證Fig.6 Verification of aerodynamic and aeroheating methods

1.4 結構動力學方程

設薄膜厚度為h，其上作用有均勻壓力p，此時薄膜將產生橫向位移z以及由變形引起的內部張力Tx和Ty。給定彈性模量為E，則靜力學控制方程為：

(2)

考慮振動效應，設薄膜單位面積的質量為ms，承受橫向壓力載荷為f(x,y,t)，則受迫振動方程為：

(3)

文獻[17]給出了氣囊和隔層上最大應力的理論公式，將本文有限元計算結果和其進行對比，如圖7所示。兩種方法得到的最大應力都隨內壓線性增加，但有限元得到的斜率稍大，其偏差均在4%左右。

圖7 靜力學方法驗證Fig.7 Verification of static stress method

本文依托北京空間機電研究所回收著陸實驗室進行了IRVE試驗件的模態振動試驗。試驗件根據IRVE的結構進行了簡化，采取三層氣囊的形式，如圖8所示。通過錘擊法測量試驗件在不同充氣壓力下的模態振型及固有頻率。

圖8 模態試驗件結構Fig.8 Structure of modal test

將試驗結果和本文所用方法進行對比，如圖9所示，兩種方法得到的結果趨勢一致。搖擺和橫向伸縮頻率隨內壓變化很小，偏差均在9%左右；縱向伸縮頻率隨內壓呈近似線性變化，偏差在7%左右。由于試驗件存在銷釘、充氣閥等連接裝置，增加了附加質量，在一定程度上降低了固有頻率，因此該偏差在合理范圍內。

圖9 模態方法驗證Fig.9 Verification of modal method

1.5 熱力學方程

TPS熱傳導分析的微分方程如下：

(4)

其中,T是溫度，c是比熱容，λ是熱傳導率，φ是內熱源。對于內充壓氣體熱膨脹，設摩爾質量M，氣體常數R，則可通過理想氣體狀態方程進行描述：

pM=ρRT

(5)

本文根據文獻[18]中的參數建立了相應的模型，并和文獻結果進行了對比，如圖10所示。兩者溫度曲線基本一致，驗證了本文熱力學方法的可靠性。

圖10 熱傳導方法驗證Fig.10 Verification of heat conduction method

1.6 飛行動力學方程

假設充氣式再入減速器以穩定的姿態飛行，無攻角無側滑，則描述其飛行狀態需要六個參數：速度大小v，地心距r，緯度φ，經度θ，彈道傾角γ和飛行方向角Ψ。考慮地球自轉，設角速度為w，其三自由度質心運動方程如下所示：

(6)

編寫MATLAB程序并和文獻[19]給出的算例進行了對比。圖11給出了對比結果，由于選用的大氣模型存在差異，在再入后半段存在一定偏差，但總體基本吻合，驗證了本程序的準確性。

圖11 飛行動力學程序驗證Fig.11 Verification of flight dynamics program

2 飛行動力學與流場分析

2.1 飛行動力學分析

表1展示了再入過程的初始參數，編寫飛行動力學程序并設定時間步長為0.1 s進行迭代求解。

表1 飛行任務相關參數Table 1 Flight parameters

圖12展示了高度、速度隨時間的變化關系。再入總時間為712 s，其中前180 s為高超聲速區，在244 s時進入亞聲速區，此時高度降低至32 km。

圖12 高度及速度曲線圖Fig.12 Height and velocity curves

2.2 流場分析

由于亞聲速時過載及熱流密度都較低，因此本文只針對超聲速及高超聲速區進行流固耦合分析，結果如圖13所示。駐點壓力峰值位于60 km高度，大小為5033 Pa；駐點溫度峰值位置有所區別，位于66 km高度處，大小為1942 K。

圖13 駐點壓力與溫度曲線圖Fig.13 Pressure and temperature curves on stagnation point

選取有代表性的Ma5和Ma1.2進行流場對比，圖14是子午面溫度云圖。Ma5飛行時結構前緣形成一道強烈的弓形激波，到達駐點時溫度從240 K升高至730 K，除前緣處氣動溫度較高外，在飛行器尾部凹陷區亦會產生嚴重的氣動加熱。相比而言，Ma1.2飛行時激波曲率更小，脫體距離更大，到達駐點時溫度僅升高至326 K。但氣流經過迎風面兩側時產生更明顯的膨脹波，溫度迅速降低至200 K。

圖14 流場溫度云圖對比Fig.14 Temperature comparison of flow field

圖15是流場子午面速度云圖,Ma5時前緣來流速度從1602 m/s起迅速降低，動能轉化為氣動熱，飛行器尾部凹陷處產生了明顯的低速區。相比而言，Ma1.2時由于迎風面兩側產生了更明顯的膨脹波，速度升高明顯并迅速超過來流速度，達到500 m/s。

圖15 流場速度云圖對比Fig.15 Velocity comparison of flow field

圖16分別展示了表面壓力和溫度在各自峰值處(60 km與66 km)的分布，圖17是兩者沿子午線的變化。發現從駐點向周圍變化的過程中，溫度和壓力的下降速度逐漸增加。到達迎風面邊緣時，溫度從1942 K下降至1169 K，壓力從5033 Pa下降至2518 Pa。從背風面來看，溫度仍然保持在較高的水平，最大達到1390 K，而壓力卻僅有230 Pa左右，因此背風面上的氣動壓力作用可忽略不計。

圖16 迎風面壓力與溫度分布Fig.16 Pressure and temperature distribution on upward

圖17 子午線壓力及溫度分布Fig.17 Distribution of meridian pressure and temperature

3 考慮氣動力作用的結構動力學分析

圖18展示了氣動力作用下的最大應力變化，此時不考慮氣動熱的影響。其規律和駐點壓力相似，在60 km處存在峰值，大小為39.6 MPa。和再入初期相比，最大應力值提升了1.84倍。下文著重選取應力峰值點A研究氣動力對應力的影響。

圖19對比了無氣動力和氣動力最大狀態(A點)的應力云圖。圖19(a)中應力完全由20 kPa內充壓產生，此時氣囊應力普遍高于蒙皮，最大應力位于最內層氣囊與隔層連接處，大小為13.9 MPa，超過內壓三個數量級。圖19(b)中在氣動力的作用下，最大應力位于迎風面靠近內層的一圈蒙皮上，大小為39.6 MPa。

圖18 氣動力作用下結構最大應力變化Fig.18 Maximal stress variation under aerodynamic

圖19 結構應力分布對比Fig.19 Comparison of stress distribution

圖20為結構子午面變形云圖的對比。圖20(a)中結構最大變形發生在各氣囊處，其幅值為4.9 mm。

圖20 結構子午面變形Fig.20 Meridian deformation of structure

氣囊在蒙皮束縛位置會將部分力傳遞給蒙皮，從而改善自身受力效果。圖20(b)中蒙皮受氣動力向內部凹陷，最大變形為3.8 mm，在束縛位置會將部分壓力傳遞給氣囊，從而改善蒙皮受力。

為比較結構中各位置應力的大小，將蒙皮、氣囊和隔層劃分成30個部分，如圖21所示，通過提取各個面上最大應力，研究不同工況下應力的分布。

圖21 結構標注示意圖Fig.21 Structural labeling sketch

圖22是氣動力作用前后各位置應力的對比，圖22(a)中隔層應力最高，氣囊次之，蒙皮應力最低。圖22(b)中蒙皮最大應力由3.6 MPa突增至39.6 MPa，氣囊最大應力也有小幅提升，由11.6 MPa提升至16.8 MPa，僅隔層最大應力有所下降，但幅度很小，從13.9 MPa下降至12.8 MPa。因此內充壓對隔層應力的影響最大，而氣動力對蒙皮上的應力影響最大。

圖22 結構各部分應力幅值對比Fig.22 Comparison of stress in different parts

研究內充壓作用下的固有振動特性，如圖23和表2所示。可以看出，前兩階固有頻率一致，為17.132 Hz，對應搖擺振型，表明此振型在X和Y兩個方向都有體現；第三階為伸縮振型，頻率為56.488 Hz。

圖23 前三階模態振型Fig.23 The first three modes

表2 前三階固有頻率Table 2 The first three natural frequencies

圖24展示了考慮氣動力時一階固頻的變化。飛行前半段氣動力不斷增大，一階固頻逐漸從17.5 Hz下降至14.1 Hz。分析該現象的原因，氣動力會削弱內充壓氣體的作用，從而減小了隔層的拉應力和結構的剛度，進而減小了固有頻率。

圖24 氣動力作用下一階固頻變化Fig.24 First order frequency variation under aerodynamic

4 考慮氣動熱作用的熱力學及結構動力學分析

圖25展示了各功能層的溫度變化，盡管駐點溫度最高達到1976 K，但TPS結構使熱流傳導減慢，最內部氣密層的溫度在有限時間內僅達到488 K，為駐點溫度的1/4，從而滿足了材料許用溫度值。

圖25 TPS各功能層溫度變化Fig.25 Temperature variation of TPS functional layers

圖26對比了氣動熱和氣動力各自對結構應力的影響。熱應力的變化存在兩個峰值，分別由于結構熱變形和內充壓氣體熱膨脹占主導所致，如圖中B點和C點所示，其應力幅值分別為33.5 MPa和31.1 MPa，分別發生在67 km 和45 km處。

圖26 氣動熱與氣動力對結構應力的影響對比Fig.26 Comparison of aeroheating and aerodynamic on stress

圖27著重選取B點和C點分析熱變形和熱膨脹對結構應力的影響規律。圖27(a)中結構最大應力發生在后蒙皮上靠近內端位置，而蒙皮越靠外，應力下降越明顯。圖27(b)中應力分布和只考慮內充壓時的應力相似，最大應力都位于最內層氣囊與隔層連接處。

圖27 結構熱應力分布對比Fig.27 Comparison of thermal stress distribution

圖28展示了B點和C點子午面變形圖。圖28(a)中最大變形位于最外層氣囊處，為4.3 mm。將變形放大10倍后可知此時結構相對形狀不發生變化，僅以內層為中心進行放大。圖28(b)中最大變形發生在氣囊處，大小為6.9 mm，而蒙皮上變形很小。

圖28 結構子午面熱變形Fig.28 Meridian thermal deformation of structure

圖29對比了B狀態和C狀態各部分應力。圖29(a)中蒙皮、氣囊和隔層最大應力分別為33.5 MPa、25.5 MPa和23.3 MPa，且各部分靠近內層的應力增加明顯，說明熱變形對最內層結構的影響最大。圖29(b)中蒙皮、氣囊和隔層最大應力分別為8.2 MPa、21.4 MPa和31.1 MPa，隔層應力最為顯著，而蒙皮應力很低。

圖29 結構各部分熱應力幅值對比Fig.29 Comparison of thermal stress in different parts

圖30展示了氣動熱作用下一階固頻的變化。其影響主要分為三部分：結構熱變形、內充壓氣體熱膨脹、材料彈性模量下降。在飛行前半段，溫度升高帶來的彈性模量下降占主導，固有頻率從16.1 Hz下降至9.6 Hz；在飛行后半段，氣體熱膨脹占主導，使固有頻率升高至24.8 Hz。

圖30 氣動熱作用下一階固頻變化Fig.30 First order frequency variation under aeroheating

5 結構參數對飛行、流場及結構特性的影響

結構參數依次影響飛行動力學、空氣動力學與結構動力學特性。本文先通過控制變量研究結構參數對各模塊的單獨影響，再利用流固耦合法研究各模塊綜合影響下的變化規律。

5.1 結構參數對各模塊的單獨影響

1)對飛行動力學的影響

通過表3可知，隨著半錐角或氣囊個數的增加，過載峰值有下降趨勢。原因在于隨著半錐角或氣囊的增加，迎風面積增大導致阻力相應增加，速度衰減程度加大，過載峰值處的速度及氣動載荷變小。

表3 結構參數對過載峰值的影響Table 3 Effect of structure on overload peak value

2)對空氣動力學的影響

令來流條件相同，由表4可知,當半錐角從50°增大至65°時，駐點壓力從4710 Pa增加至5234 Pa；平均壓力增幅更大，共增加1410 Pa。其原因在于半錐角增大時，激波增壓能力提高，強度增加，如圖31所示。此外物面法向角的增大也使蒙皮上的平均壓力有明顯的升高。氣囊個數增加時，由于前緣結構不發生變化，駐點壓力和平均壓力變化幅度很小。

表4 結構參數對氣動壓力的影響Table 4 Effect of structure on aerodynamic pressure

圖31 不同半錐角流場壓力對比Fig.31 Comparison of pressure at different half taper angles

3)對結構動力學的影響

令蒙皮上氣動力與氣動熱相同，表5說明當半錐角增大時，兩種應力均呈現上升趨勢，但氣動力的作用更為顯著，使應力從16 MPa增加至60 MPa，這是由于半錐角變大，應力在縱軸上的分量減小，為了平衡外壓，需要產生更大的應力。而增加氣囊個數不改變前緣結構，對兩種應力的影響較小。

表5 結構參數對最大應力的影響Table 5 Effect of structure on maximal stress

5.2 結構參數對各模塊的綜合影響研究

綜合上述各模塊的影響，圖32展示了不同結構對再入過程中兩種應力的影響。兩種應力的變化趨勢相同，但氣動力產生的應力變化更為明顯：半錐角每增大5°，最大應力平均增加8.3 MPa，但在60°半錐角以上時變化有所減緩；氣囊每增加一層，最大應力平均減小12.1 MPa。

圖32 再入過程下結構參數對最大應力的影響Fig.32 Effect of structure on maximal stress in reentry

綜合之前的分析，半錐角增大時，激波強度與物面法向角增大導致的氣動載荷上升占據主導，氣囊個數增加時，由于其對流場和結構的影響改變較小，此時阻力增大帶來的減速效果起主要決定作用。

6 結 論

1) 本文提出的基于飛行參數的CFD動態邊界加載方法在模擬再入過程中的氣動力時具有良好的能力。突破了現有研究中只通過對部分高度的計算結果進行插值來獲得其余高度處氣動力的思路，可以得到更為連續的氣動力變化曲線。同時阻力系數的反饋修正進一步加強了飛行動力學計算的準確性。

2) 建立的流固耦合模型既考慮了結構變形對流場的影響，還涵蓋了結構熱傳導與內充壓氣體熱膨脹的作用，能有效模擬氣動力與氣動熱對結構特性的影響。較已有方法而言該模型更全面地考慮了溫度導致的內充壓氣體狀態參數的改變。

3) 利用上述模型對典型軌跡下的結構特性進行了計算，并從氣動力與氣動熱的角度分別進行了分析。發現氣動力作用下結構最大應力從13.9 MPa升高至39.6 MPa，而由于其減弱了充壓氣體形成的剛度，使一階固頻下降3.4Hz。氣動熱作用下結構熱變形和內充壓氣體熱膨脹使應力曲線產生兩處峰值，分別為33.5 MPa和31.1 MPa，而一階固頻受熱膨脹的作用更明顯，最大升高15.2 Hz。

4) 對不同結構尺寸下的結構特性進行了對比研究，發現氣動力與氣動熱作用下的結構應力變化趨勢相同，但氣動力作用下其應力變化更明顯：半錐角每增加5°，氣動力作用下的最大應力平均增加8.3 MPa，此時激波強度與物面法向角增大導致的氣動載荷上升占據主導；氣囊數目每增加1層，氣動力作用下的最大應力平均減小12.1 MPa，此時阻力增大帶來的減速效果起主要決定作用。