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

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

（北京航天長(zhǎng)征飛行器研究所，北京，100076）

0 引 言

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

本文基于高速氣動(dòng)熱力學(xué)、傳熱傳質(zhì)學(xué)、燒蝕理論等，采用自編燒蝕移動(dòng)用戶子程序代碼方式，二次開(kāi)發(fā)Abaqus有限元軟件，利用二次開(kāi)發(fā)軟件仿真分析試驗(yàn)條件下空氣舵前緣三維燒蝕/溫度場(chǎng)，給出空氣舵前緣不同時(shí)刻燒蝕外形、燒蝕量和三維溫度分布，理論計(jì)算結(jié)果和試驗(yàn)結(jié)果進(jìn)行對(duì)比的結(jié)果表明理論計(jì)算和實(shí)測(cè)結(jié)果吻合較好，驗(yàn)證了方法的合理性和正確性。

1 控制方程和數(shù)值模型

空氣舵再入大氣層飛行時(shí)，舵前緣外表面受到強(qiáng)烈氣動(dòng)加熱，氣體動(dòng)能轉(zhuǎn)化為熱能，熱流通過(guò)空氣舵前緣表面進(jìn)入結(jié)構(gòu)內(nèi)部，從而引起結(jié)構(gòu)內(nèi)溫度場(chǎng)變化，空氣舵前緣結(jié)構(gòu)傳熱過(guò)程遵循三維熱傳導(dǎo)微分方程為

式中ρ為材料密度；Cp為材料比熱；T為溫度；t為時(shí)間坐標(biāo)；xλ，yλ，zλ分別為x、y、z方向的導(dǎo)熱系數(shù)；Q為內(nèi)熱源。

溫度場(chǎng)控制方程可從泛函變分求得，也可從導(dǎo)熱微分方程出發(fā)用加權(quán)余量法求得。加權(quán)余量法中更多采用的是伽遼金法（Calerkin）。假定空氣舵內(nèi)部無(wú)熱源，物性參數(shù)各向同性，初始溫度均勻。將空氣舵前緣溫度場(chǎng)計(jì)算域V離散為若干個(gè)單元，n個(gè)節(jié)點(diǎn)，對(duì) 式（1）加權(quán)積分得：

式中Wl為權(quán)函數(shù)，

對(duì)式（2）應(yīng)用高斯公式整理后可得：

式中D為計(jì)算面積區(qū)域；S為邊界面積；為溫度在n方向上的導(dǎo)數(shù)。

燒蝕計(jì)算一般取六面體單元，對(duì)于單元E內(nèi)任意點(diǎn)的溫度T（x,y,z）用插值函數(shù)表示：

式中N1～N8為形函數(shù)；Te～Tm為單元E節(jié)點(diǎn)溫度，對(duì)式（4）求導(dǎo)可知：

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

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

式（6）為三重積分，在笛卡爾坐標(biāo)系中進(jìn)行三重積分，通常采用坐標(biāo)變換的方法進(jìn)行，即將笛卡爾坐標(biāo)系轉(zhuǎn)化為局部坐標(biāo)系（ξ,,ηζ），即可將笛卡爾坐標(biāo)內(nèi)形狀不規(guī)則單元轉(zhuǎn)化為局部坐標(biāo)系內(nèi)形狀規(guī)則單元，從而使得積分運(yùn)算變得簡(jiǎn)化。首先引入坐標(biāo)變換：

在自然坐標(biāo)系下，式中的型函數(shù)須滿足：

式中ξi、ηi、ζi為局部坐標(biāo)系的節(jié)點(diǎn)坐標(biāo)，其中，ξ0=ξiξ，η0=ηηi，ζ0=ζiζ。

對(duì)式（4）、式（7）求偏導(dǎo)數(shù)，將式（8）代入可得式（9）、式（10）。

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

對(duì)于體積微元、面積微元應(yīng)用向量關(guān)系可得：

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

式（15）可簡(jiǎn)寫為

式中K為溫度剛度矩陣；T未知溫度列向量；N為非穩(wěn)態(tài)變溫矩陣；P為邊界載荷列向量；下標(biāo)t表示這些列向量都取同一個(gè)t時(shí)刻的值。

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

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

2 三維燒蝕移動(dòng)邊界

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

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

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

對(duì)于式（18）中的q0可由式（19）獲得：

式中qs為駐點(diǎn)熱流密度；Rn為端頭半徑；為流態(tài)因子；ρ∞為來(lái)流密度；V∞為來(lái)流速度；rc為舵前緣半徑；μ∞為氣體粘度；Λ為舵前緣后掠角；m，a，b為常數(shù)，由試驗(yàn)獲得。

對(duì)于hr有：

式中he為靜焓；r（ 0）為恢復(fù)系數(shù)，為普朗特?cái)?shù)。

ψ值可由下式獲得：

式中β為系數(shù)，層流下為0.62，湍流下為0.2；f為質(zhì)量分?jǐn)?shù)（包含樹(shù)脂質(zhì)量分?jǐn)?shù)、SiO2質(zhì)量分?jǐn)?shù)）。

2.2 移動(dòng)邊界的燒蝕速率

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

表1 試驗(yàn)條件 Tab.1 Experimental Conditions

圖1為對(duì)應(yīng)狀態(tài)下模型燒蝕后典型駐點(diǎn)形貌，可以看出，模型在表1條件下均發(fā)生燒蝕，低狀態(tài)下模型表面的白色液態(tài)層少，較高狀態(tài)下模型表面的液態(tài)層有所增加，高狀態(tài)下除被氣動(dòng)剪切力吹掉的液態(tài)層外，大部分的白色液態(tài)層仍附著在模型表面。

圖1 駐點(diǎn)燒蝕模型 Fig.1 Stagnation Ablation Model

根據(jù)表1測(cè)得的質(zhì)量燒蝕速度與熱流密度和壓力的關(guān)系，采用多元回歸的方法，得到硅基材料質(zhì)量燒蝕速率與熱流密度和壓力的關(guān)聯(lián)公式：

式中mt為質(zhì)量燒蝕速度；k為常數(shù)；0P為壓力。根據(jù)質(zhì)量燒蝕速度和線燒蝕速度的關(guān)系，可得節(jié)點(diǎn)坐標(biāo)的移動(dòng)速度：

對(duì)燒蝕速度積分可得節(jié)點(diǎn)燒蝕移動(dòng)量：

2.3 燒蝕/溫度場(chǎng)耦合計(jì)算

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

3 方法驗(yàn)證

3.1 計(jì)算模型和條件

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

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

3.2 計(jì)算與試驗(yàn)結(jié)果對(duì)比

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

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

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

圖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可知，舵前緣和舵面均出現(xiàn)燒蝕，燒蝕量沿前緣弧頂至舵面逐漸減小，前緣燒蝕量最大，為5.744 mm，位置在舵前緣圓弧頂點(diǎn)處，燒蝕后的舵前緣外形與初始外形相比，曲率半徑變小，外形變平緩。

為驗(yàn)證方法的正確性，選取舵前緣局部模型（小三角），在電弧加熱器上進(jìn)行了前緣定態(tài)模擬燒蝕試驗(yàn)。模型材料與計(jì)算模型相同，試驗(yàn)?zāi)M了計(jì)算熱環(huán)境條件的100%。圖6a為電弧加熱器試驗(yàn)?zāi)Ｐ?，試?yàn)采用矩形噴管，將試驗(yàn)?zāi)Ｐ屯ㄟ^(guò)工裝固定在噴管前方；圖6b為燒蝕后的舵前緣模型，從燒蝕試驗(yàn)后的模型可知，前緣圓弧與舵面接觸處出現(xiàn)白色熔融液態(tài)層，表明舵前緣發(fā)生了嚴(yán)重?zé)g，燒蝕量沿母線方向較為均勻。試驗(yàn)過(guò)程利用紅外點(diǎn)溫儀測(cè)量了前緣表面溫度，其值為2255 K。

圖6 舵前緣局部燒蝕試驗(yàn)Fig.6 Ablation Test of Rudder Leading Edge

試驗(yàn)后測(cè)量了前緣法向燒蝕量，最大燒蝕量為5.39 mm（虛線為原始外形），試驗(yàn)測(cè)量結(jié)果和仿真計(jì)算結(jié)果的對(duì)比見(jiàn)表2。

表2 表面溫度和燒蝕量仿真結(jié)果和試驗(yàn)值對(duì)比Tab.2 Comparision of Simulation Results and Experimental Values of Surface Temperature and Ablation

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

4 結(jié) 論

本文通過(guò)對(duì)高速空氣舵前緣三維燒蝕/溫度場(chǎng)耦合分析研究，可得如下結(jié)論：

a）基于通用有限元軟件，采用二次開(kāi)發(fā)自編熱流密度加載和燒蝕移動(dòng)邊界用戶子程序方法，可仿真分析空氣舵前緣等復(fù)雜結(jié)構(gòu)的三維燒蝕/溫度場(chǎng)，為分析空氣舵前緣熱應(yīng)力和評(píng)價(jià)空氣舵前緣熱結(jié)構(gòu)設(shè)計(jì)合理性和可靠性提供了較為精細(xì)的理論分析方法；

b）提出用熱流密度和壓力關(guān)聯(lián)式計(jì)算空氣舵前緣燒蝕后退量工程方法，解決了后掠空氣舵前緣燒蝕量計(jì)算和邊界節(jié)點(diǎn)移動(dòng)控制量問(wèn)題，同時(shí)與ALE方法有機(jī)結(jié)合，可合理地給出復(fù)雜結(jié)構(gòu)三維燒蝕外形；

c）通過(guò)舵前緣局部燒蝕試驗(yàn)可知，仿真結(jié)果和試驗(yàn)實(shí)測(cè)值較為接近，且仿真結(jié)果高于實(shí)測(cè)值，偏差小于10%，表明方法正確合理。