高速飛行器空氣舵前緣三維燒蝕/溫度耦合分析研究

楊凱威，梁 歡，趙小程，陳政偉，那 偉

（北京航天長征飛行器研究所，北京，100076）

0 引 言

空氣舵是高速滑翔和再入機動飛行器的關鍵部件，其可提供飛行器控制力矩和升力，以此實現高速飛行器機動變軌和長時間在大氣層內機動飛行。空氣舵通常突出飛行器表面且結構形式復雜，當以高速、大攻角、大側滑角在大氣層內長時間飛行時，空氣舵受到激波相互干擾，周圍流場非常復雜，力、熱環境惡劣，空氣舵出現明顯燒蝕，特別是舵前緣燒蝕更為嚴重[1]。空氣舵前緣熱防護問題突出，該問題一直是高速再入機動飛行器熱防護設計的關鍵點和難點[2]。而解決這一難點的關鍵之一是準確預測空氣舵前緣三維燒蝕/溫度場，只有準確預測了空氣舵前緣燒蝕外形、燒蝕量變化規律及舵前緣各組成結構溫度分布、溫度變化規律，才能精確設計空氣舵前緣防熱結構，合理評價空氣舵前緣防熱結構可靠性，以及外形變化對飛行器氣動特性影響等問題。美國Aerotherm公司在20世紀60年代，開發出了燒蝕熱響應計算程序（Charring Material Thermal Response and Ablation Program，CMA)，實現了燒蝕表面能量平衡方程與內部能量平衡及熱解方程之間的耦合。然而，CMA程序在熱解計算中采用的是顯式算法，在計算收斂性和穩定性上有所欠缺[3]。Chen[4]等又將該程序應用到簡單三維外形中，在穩定性和準確性方面有所提升，但該方法只能采用正交性良好的結構網格，不能很好地處理結構復雜的外形。2012 年，John A.Dec[5]利用有限元方法發展了系統的針對燒蝕材料內部多維熱響應問題的分析程序，此程序采用高階曲網格，大大增加了算法對多維曲邊外形的適應能力。2015年，Alexandre Martin[6]等嘗試在一維條件下高超聲速熱燒蝕非平衡流和材料熱響應之間的緊耦合，并在一些簡單的案例中對算法進行了驗證。2018年，Wang[7]等應用用戶子程序對一維條件下材料的熱燒蝕效應進行了預測。張濤[8]等研究了二維條件下熱解型碳化復合材料的燒蝕計算模型和計算方法，獲得了熱解氣體質量流量分布。楊德軍[9]等對碳/碳復合材料平頭柱體模型的瞬態溫度場進行了分析計算。2016年，劉驍[10]等對熱防護系統三維有限元計算方法進行了研究，其求解方法還需要通過試驗進行進一步驗證與確認。針對空氣舵前緣三維燒蝕溫度場耦合分析的研究尚未見報道。

本文基于高速氣動熱力學、傳熱傳質學、燒蝕理論等，采用自編燒蝕移動用戶子程序代碼方式，二次開發Abaqus有限元軟件，利用二次開發軟件仿真分析試驗條件下空氣舵前緣三維燒蝕/溫度場，給出空氣舵前緣不同時刻燒蝕外形、燒蝕量和三維溫度分布，理論計算結果和試驗結果進行對比的結果表明理論計算和實測結果吻合較好，驗證了方法的合理性和正確性。

1 控制方程和數值模型

空氣舵再入大氣層飛行時，舵前緣外表面受到強烈氣動加熱，氣體動能轉化為熱能，熱流通過空氣舵前緣表面進入結構內部，從而引起結構內溫度場變化，空氣舵前緣結構傳熱過程遵循三維熱傳導微分方程為

式中ρ為材料密度；Cp為材料比熱；T為溫度；t為時間坐標；xλ，yλ，zλ分別為x、y、z方向的導熱系數；Q為內熱源。

溫度場控制方程可從泛函變分求得，也可從導熱微分方程出發用加權余量法求得。加權余量法中更多采用的是伽遼金法（Calerkin）。假定空氣舵內部無熱源，物性參數各向同性，初始溫度均勻。將空氣舵前緣溫度場計算域V離散為若干個單元，n個節點，對 式（1）加權積分得：

式中Wl為權函數，

對式（2）應用高斯公式整理后可得：

式中D為計算面積區域；S為邊界面積；為溫度在n方向上的導數。

燒蝕計算一般取六面體單元，對于單元E內任意點的溫度T（x,y,z）用插值函數表示：

式中N1～N8為形函數；Te～Tm為單元E節點溫度，對式（4）求導可知：

將式（5）代入式（3）式并考慮熱流和輻射邊界整理可得：

式中qn為壁面凈熱流；ε為輻射系數；σ為斯蒂芬-波爾茲曼常數；TW為表面溫度。

式（6）為三重積分，在笛卡爾坐標系中進行三重積分，通常采用坐標變換的方法進行，即將笛卡爾坐標系轉化為局部坐標系（ξ,,ηζ），即可將笛卡爾坐標內形狀不規則單元轉化為局部坐標系內形狀規則單元，從而使得積分運算變得簡化。首先引入坐標變換：

在自然坐標系下，式中的型函數須滿足：

式中ξi、ηi、ζi為局部坐標系的節點坐標，其中，ξ0=ξiξ，η0=ηηi，ζ0=ζiζ。

對式（4）、式（7）求偏導數，將式（8）代入可得式（9）、式（10）。

式中J為雅克比（Jacobi）矩陣，其表達式為

對于體積微元、面積微元應用向量關系可得：

將式（9）～（10）、式（12）～（14）代入式（6），并將單元矩陣總體合成即可得到整個求解區域n個代數方程，將代數方程寫為矩陣形式為

式（15）可簡寫為

式中K為溫度剛度矩陣；T未知溫度列向量；N為非穩態變溫矩陣；P為邊界載荷列向量；下標t表示這些列向量都取同一個t時刻的值。

對式（16）中｛?T?t｝t采用向后有限差分離散，整理后可得空氣舵前緣三維溫度場有限元方程：

式中 ｛T｝t-Δt為初始溫度場或前一時刻的溫度場，為已知量，采用牛頓-拉普森（Newton-Raphson）迭代方法即可求得空氣舵前緣各節點溫度值。

2 三維燒蝕移動邊界

2.1 燒蝕表面凈熱流邊界條件

上述固體導熱微分方程對于無燒蝕邊界易于求解，但對于有燒蝕移動邊界問題，求解就變得異常復雜。目前高速空氣舵熱防護通常采用燒蝕式熱防護，防熱材料大都為硅基類材料。該類材料通常以二氧化硅為主要成分按不同工藝方法與樹脂復合而成，其燒蝕防熱主要機理是依靠材料的質量損失和物理化學變化來消耗氣動加熱以達到防熱目的。依據能量守恒條件可寫出燒蝕條件下空氣舵前緣表面凈熱流密度：

式中qN為凈熱流密度；為表面溫度梯度；ψ為質 量引射系數；q0為不考慮壁面焓值影響的表面熱流密度；hw為熱壁焓；hr為恢復焓；ε為材料輻射系數；Tw為表面溫度；v-∞為線燒蝕速率；ΔH為損失單位質量帶走的熱量（包括液態層流失、樹脂分解熱、SiO2蒸發吸熱），其值可由試驗獲得。

對于式（18）中的q0可由式（19）獲得：

式中qs為駐點熱流密度；Rn為端頭半徑；為流態因子；ρ∞為來流密度；V∞為來流速度；rc為舵前緣半徑；μ∞為氣體粘度；Λ為舵前緣后掠角；m，a，b為常數，由試驗獲得。

對于hr有：

式中he為靜焓；r（ 0）為恢復系數，為普朗特數。

ψ值可由下式獲得：

式中β為系數，層流下為0.62，湍流下為0.2；f為質量分數（包含樹脂質量分數、SiO2質量分數）。

2.2 移動邊界的燒蝕速率

通過式（18）可發現，要實現表面凈熱流密度加載，就需確定當前熱環境條件下的燒蝕速度。目前計算硅類防熱材料燒蝕速度方法主要是針對球-錐飛行器外形，而對于凸出飛行器表面的空氣舵外形特別是空氣舵前緣該方法已不適用。本文提出采用質量燒蝕速率與來流熱流密度和壓力關聯式來計算燒蝕速度的方法[11]。選取空氣舵前緣材料加工為駐點模型，試驗設備采用電弧加熱器，試驗共進行了30余次，試驗條件見表1。

表1 試驗條件 Tab.1 Experimental Conditions

圖1為對應狀態下模型燒蝕后典型駐點形貌，可以看出，模型在表1條件下均發生燒蝕，低狀態下模型表面的白色液態層少，較高狀態下模型表面的液態層有所增加，高狀態下除被氣動剪切力吹掉的液態層外，大部分的白色液態層仍附著在模型表面。

圖1 駐點燒蝕模型 Fig.1 Stagnation Ablation Model

根據表1測得的質量燒蝕速度與熱流密度和壓力的關系，采用多元回歸的方法，得到硅基材料質量燒蝕速率與熱流密度和壓力的關聯公式：

式中mt為質量燒蝕速度；k為常數；0P為壓力。根據質量燒蝕速度和線燒蝕速度的關系，可得節點坐標的移動速度：

對燒蝕速度積分可得節點燒蝕移動量：

2.3 燒蝕/溫度場耦合計算

當表面溫度達到空氣舵前緣材料燒蝕溫度時，空氣舵前緣發生燒蝕，其表面不斷向后退縮，熱量通過表面向內部傳導，從而引起內部結構溫升。燒蝕/溫度場耦合計算的總體思路是：某一時刻表面溫度達到材料燒蝕溫度時，通過式（22）計算出材料燒蝕量，進而獲得燒蝕帶走的熱量ρv-∞ΔH，通過式（24）可獲得各節點燒蝕移動量，保持網格拓撲關系不變，運用網格移動技術，重新形成新外形，并采用插值方式將各節點溫度信息插值到新外形節點上，再將新的熱流邊界加載至新外形上，重復上述步驟，即可實現各個時刻燒蝕/溫度耦合計算。對于熱流加載，本文開發了Abaqus有限元用戶子程序 DFLUX，對于節點燒蝕移動量開發了用戶子程序UMESHMOTION。由于燒蝕屬于網格大變形問題，通過采用ALE（任意拉格朗日歐拉）網格技術，實現了燒蝕區域網格平順變形。開發后的軟件計算耦合條件下燒蝕/溫度主要分為3步：a）讀入熱環境參數熱流q0、恢復焓hr、壁面壓力p；b）根據初始溫度，利用DFLUX加載凈熱流邊界計算出Tw以及空氣舵前緣內部溫度；c）用計算所得表面溫度與燒蝕判別溫度比較，通常硅基材料的燒蝕判別溫度是指材料表面出現熔融液態層時的表面溫度，不同的硅基材料燒蝕判別溫度有所不同，應通過試驗獲得燒蝕判別溫度。表面溫度若小于燒蝕判別溫度，溫度場按式（18）中ρv-∞ΔH=0的qN加載計算；若大于燒蝕判別溫度則完全按式（18）中的qN加載計算溫度場。

3 方法驗證

3.1 計算模型和條件

選取空氣舵前緣局部建立有限元模型，如圖2所示。舵前緣母線長度165 mm，高度60 mm，舵前緣后掠角Λ=68.7°，端頭半徑Rn=20 mm，舵前緣半徑rc= 13 mm。計算模型為硅基復合材料，密度ρ=1670 kg/m3，導熱系數λ=0.65 W/(m·K)，比燒蝕焓ΔH=45 kJ/kg，發射率ε=0.85，舵前緣熱流密度q=16 208 kW/m2，焓值hr=5898 kJ/kg，壓力p=0.7 MPa，計算時間為 8 s。

圖2 舵前緣局部有限元模型 Fig.2 A Portion of Rudder Leading Edge Finite Element Model

3.2 計算與試驗結果對比

圖3為數值計算的空氣舵前緣三維燒蝕溫度云圖。由圖3可知，舵前緣最高溫度達到2363 K，超過了硅基材料熔融溫度，空氣舵出現明顯燒蝕，舵前緣網格出現后退，由于加載的熱流較為均勻，所以沿前緣母線方向的燒蝕量變化也較為均勻。

圖3 舵前緣三維燒蝕溫度云圖Fig.3 Three-dimensional Ablation Temperature Contour of Rudder Leading Edge

為更好觀察空氣舵前緣燒蝕溫度分布情況，沿舵前緣1/2位置處截取溫度分布圖，如圖4所示。由圖4可知，舵前緣和舵面溫度均達到了材料燒蝕溫度，發生燒蝕，前緣溫度高于舵面溫度，計算模型法線方向溫度梯度最大，高溫區均集中在燒蝕表面，主要是由于硅基材料熱導率低。

圖4 舵前緣截面溫度云圖Fig.4 Cross-sectional Temperature Contour of Rudder Leading Edge

圖5為舵前緣1/2位置處燒蝕量截面云圖。

圖5 舵前緣截面燒蝕云圖Fig.5 Cross-sectional Ablation Contour of Rudder Leading Edge

由圖5可知，舵前緣和舵面均出現燒蝕，燒蝕量沿前緣弧頂至舵面逐漸減小，前緣燒蝕量最大，為5.744 mm，位置在舵前緣圓弧頂點處，燒蝕后的舵前緣外形與初始外形相比，曲率半徑變小，外形變平緩。

為驗證方法的正確性，選取舵前緣局部模型（小三角），在電弧加熱器上進行了前緣定態模擬燒蝕試驗。模型材料與計算模型相同，試驗模擬了計算熱環境條件的100%。圖6a為電弧加熱器試驗模型，試驗采用矩形噴管，將試驗模型通過工裝固定在噴管前方；圖6b為燒蝕后的舵前緣模型，從燒蝕試驗后的模型可知，前緣圓弧與舵面接觸處出現白色熔融液態層，表明舵前緣發生了嚴重燒蝕，燒蝕量沿母線方向較為均勻。試驗過程利用紅外點溫儀測量了前緣表面溫度，其值為2255 K。

圖6 舵前緣局部燒蝕試驗Fig.6 Ablation Test of Rudder Leading Edge

試驗后測量了前緣法向燒蝕量，最大燒蝕量為5.39 mm（虛線為原始外形），試驗測量結果和仿真計算結果的對比見表2。

表2 表面溫度和燒蝕量仿真結果和試驗值對比Tab.2 Comparision of Simulation Results and Experimental Values of Surface Temperature and Ablation

由表2可知，理論計算表面溫度和燒蝕量均大于實測值，表面溫度高出實測值108 K，理論值和實測值相對偏差為 4.79%，燒蝕量理論值高出實測值0.354 mm，理論值和實測值相對偏差為 6.56%，滿足高速熱防護設計燒失量計算偏差要求，一般要求 10%以內。

4 結 論

本文通過對高速空氣舵前緣三維燒蝕/溫度場耦合分析研究，可得如下結論：

a）基于通用有限元軟件，采用二次開發自編熱流密度加載和燒蝕移動邊界用戶子程序方法，可仿真分析空氣舵前緣等復雜結構的三維燒蝕/溫度場，為分析空氣舵前緣熱應力和評價空氣舵前緣熱結構設計合理性和可靠性提供了較為精細的理論分析方法；

b）提出用熱流密度和壓力關聯式計算空氣舵前緣燒蝕后退量工程方法，解決了后掠空氣舵前緣燒蝕量計算和邊界節點移動控制量問題，同時與ALE方法有機結合，可合理地給出復雜結構三維燒蝕外形；

c）通過舵前緣局部燒蝕試驗可知，仿真結果和試驗實測值較為接近，且仿真結果高于實測值，偏差小于10%，表明方法正確合理。