季衛(wèi)棟，王江峰，唐國慶

（南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016）

高超聲速多體分離過程氣動加熱計算技術(shù)

季衛(wèi)棟，王江峰，唐國慶

（南京航空航天大學(xué) 航空宇航學(xué)院，江蘇 南京 210016）

發(fā)展了無粘外流解與氣動熱工程方法相結(jié)合的用于全機外形復(fù)雜流動氣動加熱計算技術(shù)，并將其與動態(tài)多體分離過程流場解相結(jié)合，提出了一種高超聲速飛行器多體分離過程氣動加熱特性計算技術(shù)，然后耦合結(jié)構(gòu)傳熱，實現(xiàn)了飛行器熱防護系統(tǒng)中結(jié)構(gòu)材料溫度分布特性的數(shù)值模擬。采用本文方法，針對假定高超聲速多體分離布局及飛行條件開展了分離過程中的結(jié)構(gòu)溫度分布特性與熱流密度分布特性計算，分析并給出了組合體分離過程中部件干擾對氣動熱的影響特性。結(jié)果表明本文發(fā)展的計算技術(shù)可為高超聲速飛行器的防熱方案設(shè)計、氣動熱特性定性分析及熱防護系統(tǒng)設(shè)計等提供技術(shù)支持。

高超聲速；多體分離；氣動加熱；工程方法；結(jié)構(gòu)傳熱

0 引 言

飛行器氣動加熱特性的計算與分析是高超聲速飛行器設(shè)計中的關(guān)鍵問題之一，隨著高超聲速技術(shù)的迅速發(fā)展，高馬赫中低空多體分離方案越來越受到重視，此設(shè)計方案的關(guān)鍵問題之一是分離過程中的氣動熱特性的快速準(zhǔn)確預(yù)測。

目前在高超聲速氣動加熱計算方法研究方面，針對單體外形布局或飛行器某一部件的工程方法［1-5］、數(shù)值方法［6-8］和工程與數(shù)值相結(jié)合的方法［9-10］的發(fā)展已相對較成熟。隨著高超聲速飛行器技術(shù)的發(fā)展，多體分離問題涉及到的高度越來越低、分離時的飛行速度越來越快，多體分離過程中的氣動加熱問題逐漸變得重要起來，目前國內(nèi)外對這一問題的研究多為實驗方面［11］，因此，發(fā)展一種多體分離動態(tài)氣動加熱計算技術(shù)就變得非常有必要。

本文將單體飛行器布局氣動加熱計算技術(shù)與多體分離過程的氣動熱計算相結(jié)合，提出了一種多體分離動態(tài)氣動加熱計算技術(shù)，對多體分離過程中的氣動加熱問題進行計算與分析。單體飛行器布局氣動加熱計算技術(shù)是基于普朗特邊界層理論，將整個流場分為邊界層外的無粘流動和邊界層內(nèi)的粘性流動兩部分，使用數(shù)值方法計算無粘解，氣動熱工程方法計算邊界層傳熱，并耦合結(jié)構(gòu)傳熱［12-13］，發(fā)展了一種可用于全機外形復(fù)雜流場氣動加熱的計算方法；針對多體分離問題，提出了不依賴于計算網(wǎng)格類型和分離體個數(shù)的分離體劃分策略與熱防護參數(shù)設(shè)置方法，結(jié)合本文發(fā)展的氣動熱計算技術(shù)，形成了針對高超聲速多體分離過程中各分離體表面熱流密度分布以及防熱層的結(jié)構(gòu)溫度分布的計算方法。這種方法結(jié)合了數(shù)值模擬和純工程方法各自的優(yōu)點，提高了數(shù)值模擬的效率，擴充了純工程方法的應(yīng)用范圍，能方便的應(yīng)用于高超聲速多體分離的氣動力、氣動熱的快速分析。

采用本文提出的數(shù)值計算技術(shù)，針對假定高超聲速多體分離布局，進行了多體分離過程中的結(jié)構(gòu)溫度分布特性與熱流密度分布特性的數(shù)值模擬，給出了動態(tài)分離過程中的計算與分析結(jié)果。

1 復(fù)雜外形氣動加熱計算技術(shù)

根據(jù)普朗特邊界層理論，在高超聲速條件下，流體的粘性和導(dǎo)熱可以看成集中作用在物體表面的薄層即邊界層內(nèi)。因此，可將整個流場分為邊界層外的無粘流動和邊界層內(nèi)的粘性流動兩部分，而在使用工程方法對邊界層內(nèi)傳熱問題進行分析時，只需要邊界層外緣參數(shù)的分布，由此，對氣動熱的求解也可以化為高超聲速無粘流場的求解和邊界層內(nèi)粘性主導(dǎo)區(qū)域的氣動熱的求解。

1.1 無粘外流場的求解

在表面熱流計算中，所關(guān)心的僅僅是邊界層外緣氣動特性的分布，并且在高超聲速下的邊界層很薄，因此本文在計算邊界層外緣參數(shù)時，采用了直接求解Euler方程的方法，將Euler方程數(shù)值解的物面參數(shù)分布作為粘性流動邊界層的外緣參數(shù)。

本文采用課題組已有重疊網(wǎng)格并行計算程序完成分離過程無粘外流解計算，輸出物面網(wǎng)格點上的流動參數(shù)，如速度、壓強、密度等，再根據(jù)完全氣體狀態(tài)方程，計算出其他所需流動參數(shù)，如溫度、馬赫數(shù)等。

1.2 表面熱流密度計算

表面熱流密度的計算分為駐點區(qū)和非駐點區(qū)兩部分分別計算。

駐點區(qū)熱流密度的計算，已有多種適用于不同條件的成熟計算方法，如：SCALA方法、LEES方法、Eay-Riddell方法［14］、Kemp-Riddl方法等。本文采用Eay-Riddell方法。Eay-Riddell公式主要用于軸對稱的平衡邊界層駐點熱流密度的計算，公式如下：

或

其中Pr＝0．71，Le＝1．0。

非駐點區(qū)的表面熱流計算，在工程上也有局部相似解法（LEES鈍頭體層流熱流密度分布公式）和參考焓方法等方法，本文采用平板傳熱模型計算傳熱系數(shù)及熱流密度，對較復(fù)雜外形繞流的計算進行簡化處理。

1.3 結(jié)構(gòu)傳熱計算

為了達(dá)到快速預(yù)測的目的，防熱層結(jié)構(gòu)溫度分布根據(jù)一維熱傳導(dǎo)方程構(gòu)建差分方程求解。根據(jù)防熱層畢奧數(shù)的大小，防熱層分為熱薄壁和熱厚壁分別使用不同的傳熱計算公式［15］。熱薄壁傳熱公式如下：

式中ρδ為材料密度；cδ為材料比熱容；δ為材料厚度；α為傳熱系數(shù)；ε為材料輻射系數(shù)。

2 氣動加熱與結(jié)構(gòu)傳熱耦合計算

高超聲速流動中，飛行器受到的氣動加熱過程是一個持續(xù)非穩(wěn)態(tài)過程，壁面熱邊界條件在達(dá)到平衡前處于時刻變化中。邊界層外的氣流向飛行器表面?zhèn)魅氲臒崃髅芏萹w除了與飛行器氣動參數(shù)有關(guān)，還與物面的溫度Tw有關(guān)，而Tw的計算又與熱流密度q w相關(guān)，因此，要準(zhǔn)確模擬飛行器受熱情況，氣動加熱和結(jié)構(gòu)傳熱必須進行耦合計算［16］。

其中Taw為絕熱壁溫。

氣動加熱和結(jié)構(gòu)傳熱耦合計算步驟如圖1所示。

圖1 氣動加熱和結(jié)構(gòu)傳熱耦合計算步驟Fig.1 Calculation steps of aerodynamic heating coupled with structure heating

按照上述方法，只要給定初始時刻壁溫，按時間步長推進求解，就可進行高超聲速全機外形復(fù)雜流動下的氣動力、氣動熱及結(jié)構(gòu)溫度分布特性的快速計算。在此方面，課題組已經(jīng)發(fā)展了相應(yīng)的軟件系統(tǒng)“高速飛行器氣動加熱計算與分析”V4．0。

3 多體分離過程氣動加熱計算技術(shù)

隨著高超聲速飛行器技術(shù)的發(fā)展，可以預(yù)測高馬赫數(shù)中低空多體分離方案必將成為發(fā)展方向之一，必然涉及多體分離過程中氣動加熱的快速預(yù)測技術(shù)，目前國內(nèi)外對此問題的數(shù)值模擬方法研究還較少。

本文在國內(nèi)首次對多體分離過程中的氣動加熱計算技術(shù)開展了研究，發(fā)展了分離體劃分策略、防熱材料設(shè)置等方法，并首次完成了多體分離過程氣動加熱計算軟件。

3.1 分離體劃分策略

本文提出的方法是：將飛行器的每一個分離的整體作為一個分離體，首先確定分離體的個數(shù)，然后確定每一個分離體所包含的部件序號，這些參數(shù)可在程序輸入卡中進行設(shè)置，對分離體個數(shù)、計算網(wǎng)格類型等沒有限制。

計算時，運行程序中的物面提取模塊，可提取出無粘外流解中的物面部件，根據(jù)物面網(wǎng)格信息，完成分離體的劃分。

3.2 分離體熱防護參數(shù)設(shè)置方法

完成分離體的劃分后，整個飛行器按部件分為了不同的分離體，對每一分離體單獨進行熱防護參數(shù)設(shè)置，設(shè)置方法與單體熱防護設(shè)置方法相同，用選擇計算部件加上坐標(biāo)范圍劃分熱防護區(qū)域，進行熱防護參數(shù)設(shè)置，如圖2。

圖2 熱防護系統(tǒng)設(shè)置Fig.2 Thermal protection system setting

4 算例與分析

4.1 RAMC-II典型外形氣動加熱驗證

飛行器外形采用RAMC-II試驗?zāi)Ｐ?，總體外形為球錐體，頭部半徑為0．1524m，半錐角為9°，總長度為1．295m。分別進行巡航狀態(tài)瞬時傳熱、長時加熱計算。本文圖中熱流密度單位W／m2，溫度單位K。

4.1.1 RAMC-II巡航狀態(tài)瞬時傳熱

巡航狀態(tài)瞬時傳熱，設(shè)定巡航狀態(tài)為：迎角0°、飛行高度60km，Ma＝4．0、5．0、6．0，初始壁溫247 K。表1給出了本文駐點熱流密度計算結(jié)果與工程算法Kemp-Riddell駐點加熱公式計算結(jié)果的對比，誤差在10%以內(nèi)。

表1 駐點熱流密度計算結(jié)果對比Table 1 Comparison of stagnation point’s heat flux

4.1.2 RAMC-II巡航狀態(tài)長時加熱

巡航狀態(tài)長時加熱，設(shè)定巡航狀態(tài)為：迎角6°、飛行高度60km，Ma＝6．0，初始壁溫247K，巡航時間1000s，飛行器防熱層為2mm的Ti，表面發(fā)射率為0．8。巡航第1000s時駐點溫度為945K，外表面溫度分布如圖3（a），子午面溫度分布如圖3（b），其中T＿RAD為使用同樣外流解直接計算輻射平衡溫度結(jié)果，可以看到，子午面上溫度分布與輻射平衡溫度符合很好。

圖3 巡航第1000s氣動熱參數(shù)分布Fig.3 Aerodynamic heating variables distribution（t＝1000s）

4.2 多體分離氣動加熱計算

由于缺乏類似算例對比，本文假定了一個機體／掛彈組合體方案。飛行方案為：組合體巡航飛行150s時釋放掛彈；計算要求為：從0s開始至安全分離時刻整個過程中的機體（JT）、掛彈（DD）氣動熱特性。設(shè)計計算狀態(tài)與初始參數(shù)設(shè)定如下：

（a）組合體以馬赫數(shù)5．0、迎角3°、飛行高度40km巡航飛行150s；計算從飛行時間0s～150s的全機氣動加熱特性；初始時刻組合體表面溫度為250K。

（b）第150s掛彈分離，在氣動力與重力作用下按六自由度剛體運動。分離時刻參數(shù)為：彈拋出速度0m／s、彈初始姿態(tài)角3°、彈初始姿態(tài)角速度為0°／s。分離時間段為150s～160s。計算整個分離過程中機體與掛彈的氣動加熱特性。分離初始時刻分離體結(jié)構(gòu)溫度分布及氣動熱環(huán)境參數(shù)為為巡航段尾端第150s時刻的值。

（c）分離體劃分策略：機體和掛彈分別作為一個分離體，根據(jù)網(wǎng)格信息，每個分離體包含一個網(wǎng)格部件。

（d）機體與掛彈采用相同的熱防護方案，防熱材料均為厚度為2mm的Ti，表面發(fā)射率為0．8。

圖4 飛行器外形Fig.4 Aircraft configuration

4.2.1 組合體分離過程氣動加熱

圖5為飛行器巡航第150s時全表面溫度分布和熱流密度分布。由計算結(jié)果可知機體在頭部和機翼前緣氣動加熱嚴(yán)重，溫度較高，最高頭部溫度為1004 K，掛彈在彈頭處達(dá)到最高溫度1088K，巡航第150s時機體頭部和掛彈頭部溫度基本達(dá)到平衡，熱流密度在頭部均接近為0W／m2。

圖5 巡航第150s全表面氣動熱參數(shù)分布Fig.5 Distribution of surface aerodynamic heating variables（t＝150s）

圖6為飛行器從150s開始分離，直至安全分離整個過程中飛行器外表面溫度分布和熱流密度分布，分離過程中機體在頭部達(dá)到最高溫度1004 K，掛彈在彈頭達(dá)到最高溫度1112K。

圖6 分離過程全表面氣動熱參數(shù)分布Fig.6 Distribution of surface aerodynamic heating variablesduring separation

圖7（a）和圖7（b）分別給出了機體頭部一點和掛彈頭部一點的表面溫度和熱流密度隨時間變化曲線（JT表示機體、DD表示外掛導(dǎo)彈）。

整個氣動加熱計算過程在普通臺式計算機上CPU耗時約25min（不記無粘外流解計算時間），可知本文發(fā)展的多體分離動態(tài)氣動加熱計算技術(shù)，可對復(fù)雜外形高速飛行器分離過程氣動加熱進行快速計算與分析。

圖7 機體頭部和掛彈頭部特征點氣動熱參數(shù)變化曲線Fig.7 Variation of aerodynamic heating variables of characteristic point on JT and DD

4.2.2 分離過程干擾區(qū)內(nèi)氣動加熱

分離過程中，各分離體對其他分離體的流場會產(chǎn)生干擾，對其氣動加熱產(chǎn)生影響。

圖8為分離過程中機體與掛彈相接附近不同時刻的熱流密度分布。圖9（a）為機體受掛彈干擾區(qū)內(nèi)兩個點上的熱流密度隨時間變化，Point1在左，Point2在右；圖9（b）為掛彈彈身上兩個點上的熱流密度隨時間變化，Point3在左，Point4在右。

圖8 分離過程機體與掛彈熱流密度分布（單位：W／m2）Fig.8 Heat fluxdistributionduring separation

圖9 機體和掛彈點上熱流密度隨時間變化（單位：W／m2）Fig.9 Variation of heat flux of points on JT and DD

從圖8和圖9（a）中可以看出，機身受干擾區(qū)域內(nèi)，分離開始，高熱流區(qū)向后移動，Point1最大熱流出現(xiàn)時間為151.89s，最大值7095W／m2，Point2最大熱流出現(xiàn)時間為153.69s，最大值7668 W／m2。從圖9（b）中可以看到，彈身上兩點的在分離過程中熱流密度值出現(xiàn)躍增，峰值熱流密度均超過初始時刻熱流密度，其中Point4達(dá)到37505 W／m2，遠(yuǎn)遠(yuǎn)大于初始時刻熱流密度22771 W／m2。

可見，分離過程中各分離體對其他分離體在氣動熱方面會產(chǎn)生不容忽視的影響。

4.3 不同熱防護參數(shù)下分離過程氣動加熱

采用與算例4．2相同的計算外形與外流解，對機體和掛彈分別使用不同的熱防護方案，對分離過程氣動加熱進行分析。

4.3.1 分離參數(shù)

組合體巡航飛行1000s后，釋放掛彈，巡航參數(shù)和分離時刻參數(shù)與算例4．2相同。

4.3.2 分離體劃分及熱防護區(qū)域劃分

機體和掛彈分別作為一個分離體，每個分離體按坐標(biāo)劃分為兩個熱防護區(qū)域，如圖10所示，機體將頭部和機身分為兩個熱防護區(qū)域（Tps1和Tps2），掛彈將彈頭和彈身分為兩個熱防護區(qū)域（Tps3和Tps4）。

圖10 熱防護區(qū)域劃分Fig.10 Thermal protection areadivision

4.3.3 熱防護參數(shù)設(shè)置

對實際熱防護系統(tǒng)進行模擬，不同的熱防護區(qū)域使用不同的熱防護方案，如表2，初始物面溫度250 K。

表2 熱防護類型和材料設(shè)置Table 2 Thermal protection type and material settings

4.3.4 計算結(jié)果

圖11（a）、圖11（b）和圖11（c）分別為分離初始時刻1000s時飛行器防熱層外、內(nèi)表面溫度分布和熱流密度分布。圖12為分離過程熱流密度分布。

圖11 分離初始時刻氣動熱參數(shù)分布Fig.11 Aerodynamic heating variablesdistribution（t＝1000s）

圖12 分離過程熱流密度分布（單位：W／m2）Fig.12 Heat fluxdistributionduring separation

5 結(jié) 論

本文耦合課題組發(fā)展的高超聲速長時氣動加熱快速計算技術(shù)與多體分離數(shù)值模擬結(jié)果，完成了高超聲速多體分離過程飛行器全機、部件等的氣動加熱特性計算技術(shù)。論文結(jié)論如下：

（1）結(jié)合無粘外流解和氣動熱工程算法，耦合結(jié)構(gòu)傳熱計算，發(fā)展了一種針對全機復(fù)雜外形的長時氣動加熱計算技術(shù)，完成了“高速飛行器氣動加熱計算與分析”軟件，并在相關(guān)航空航天型號設(shè)計單位得到了成功應(yīng)用。

（2）提出了不依賴于計算網(wǎng)格類型和分離體個數(shù)的多體分離動態(tài)氣動熱加熱計算方法，提出了分離體劃分策略和按分離體進行熱防護參數(shù)設(shè)置的方法，將多體分離和氣動加熱計算技術(shù)相結(jié)合，實現(xiàn)了分離過程中氣動加熱的數(shù)值模擬。

［1］楊愷，高效偉．高超聲速氣動熱環(huán)境工程算法［J］．導(dǎo)彈與航天運載技術(shù)，2010，（4）：19-23．

［2］ HAMILTON H H，WEILMUENSTER K J，DEJARNETTE E R．Approximate method for computing laminar and turbulent convective heating on hypersonic vehicles using unstructured grids［R］．In：41st Annual AIAA Thermophysics Conference［C］．AIAA 2009-4310，2009．

［3］薛鵬飛，龔春林，谷良賢．基于軸對稱比擬的氣動熱計算方法研究［J］．計算機仿真，2012，29（6）：97-100．

［4］ZHAO J S，GU L X，MA H Z．A rapid approach to convective aeroheating prediction of hypersonic vehicles［J］．Science China Technological Sciences，2013：1-15．

［5］中國人民解放軍總裝備部軍事訓(xùn)練教材編輯工作委員會．高超聲速氣動熱和熱防護［M］．北京：國防工業(yè)出版社，2003：46-136．

［6］SINHA K K，REDDY D S．Effect of chemical reaction rates on aeroheating predictions of re-entry flows［J］．Journal of Thermophysics and Heat Transfer，2011，25（1）：21-33．

［7］鞏偉杰，唐碩，李世珍．高超聲速飛行器氣動加熱三維數(shù)值分析方法研究［J］．飛行力學(xué)，2011，29（2）：78-81．

［8］PRAKASH A，ZHONG X．Numerical simulation of planetary reentry aeroheating over blunt bodies with non-equilibrium reacting flow and surface reactions［R］．In：Proceeding of the 47th AIAA Aerospace Sciences Meeting Including The New Horizons Eorum and Aerospace Exposition［C］．Orlando，Elorida．AIAA 2009-1542．

［9］H AMILTON H H，GREENE E A，DEJARNETTE E R．Approximate method for calculating heating rates on three-dimensional vehicles［J］．Journal of Spacecraft and Rockets，1994，31（3）：345-354．

［10］王杰．高超聲速飛行器氣動加熱計算技術(shù)［D］．［碩士學(xué)位論文］．南京：南京航空航天大學(xué)，2011．

［11］LIECHTY D S．Aeroheating characteristics for a two-stage-toorbit conceptduring separation at Mach 6［R］．AIAA Paper 2005-5139，2005．

［12］BOVA S W，HOWARD M A．Coupling strategies for highspeed aeroheating problems［R］．Sandia National Laboratories，2011．

［13］任青梅，楊志斌，成竹，等．氣動加熱與熱響應(yīng)耦合分析及試驗研究［J］．航天器環(huán)境工程ISTIC，2012，29（3）：235-242．

［14］EAY J A．Theory of stagnation point heat transfer indissociated air［J］．Journal of the Aerospace Sciences，1958，25（2）：73-85．

［15］QJ1734-89．戰(zhàn)術(shù)導(dǎo)彈氣動加熱工程計算方法［S］．中華人民共和國航天工業(yè)部，1989．

［16］吳潔，閻超．氣動熱與熱響應(yīng)的耦合研究［J］．導(dǎo)彈與航天運載技術(shù)，2009，（4）：35-39．

Calculating method of aerodynamic heating for hypersonic multi-body separation process

JI Weidong，WANG Jiangfeng，TANG Guoqing
（Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China）

In this paper，a method of aerodynamic heating prediction for hypersonic multi-body separation process is presented．The inviscid outflow solutionduring multi-body separation and the engineering method for aerodynamic heating are combined，and then the structure heating coupling effect is considered to calculate the temperaturedistribution of thermal protection systems．Eor assumed layout of multi-body separation aircraft，this method accomplishes the calculation of temperaturedistribution characteristics and heat fluxdistribution characteristicsduring the separation，and analyses the influence of multi-body separation on aerodynamic heating．The research results show that the computing technologydeveloped in this paper can provide technical support for hypersonic aircraftdesign both on the thermal protection systems and aerodynamic heating characteristics analysis．

hypersonic；multi-body separation；aerodynamic heating；engineering method；structure heating

V211.3

Adoi：10.7638／kqdlxxb-2014.0080

0258-1825（2014）06-0761-06

2014-08-13；

2014-10-11

江蘇高校優(yōu)勢學(xué)科建設(shè)工程資助項目

季衛(wèi)棟（1987-），男，江蘇人，博士研究生，研究方向：高超聲速氣動力／熱／結(jié)構(gòu)多場耦合．E-mail:kongzhizi9595＠163.com

季衛(wèi)棟，王江峰，唐國慶．高超聲速多體分離過程氣動加熱計算技術(shù)［J］．空氣動力學(xué)學(xué)報，2014，32（6）：761-766．

10．7638／kqdlxxb-2014．0080 JI W D，WANG J E，TANG G Q．Calculating method of aerodynamic heating for hypersonic multi-body separation process［J］．ACTA Aerodynamica Sinica，2014，32（6）：761-766．