考慮邊界層轉(zhuǎn)捩的復(fù)雜外形火箭彈氣動熱計算

章 輝，張向洪

(南京理工大學(xué) 能源與動力工程學(xué)院, 南京 210094)

邊界層轉(zhuǎn)捩是指邊界層流動由層流發(fā)展為湍流的過程，是一個多因素耦合的強非線性復(fù)雜流動的物理現(xiàn)象。高超聲速飛行器表面邊界層轉(zhuǎn)捩長期以來一直是空氣動力學(xué)領(lǐng)域研究的重點和難點，考慮邊界層轉(zhuǎn)捩和氣動加熱相結(jié)合的研究結(jié)果則更難以預(yù)測。對于轉(zhuǎn)捩點的預(yù)測目前有半經(jīng)驗的eN方法[1]、轉(zhuǎn)捩模式[2]、轉(zhuǎn)捩準(zhǔn)則[3]、直接數(shù)值模擬[4]和大渦模擬[5]，而對轉(zhuǎn)捩和湍流氣動加熱的研究主要是通過風(fēng)洞實驗和轉(zhuǎn)捩準(zhǔn)則。陳鵬[6]通過Batt準(zhǔn)則預(yù)測了高速彈頭的轉(zhuǎn)捩起始點和轉(zhuǎn)捩長度，計算了彈頭表面的氣動熱。楊愷[7]通過給定臨界雷諾數(shù)，將飛行器劃分為層流和湍流區(qū)域，計算分析了鈍錐和航天飛機的氣動加熱情況。高瑩瑩等[8]針對Reentry+F飛行試驗測量結(jié)果，應(yīng)用變熵膨脹法和變熵流線法計算并與試驗結(jié)果對比，獲得了適用Reentry+F飛行器的轉(zhuǎn)捩判定準(zhǔn)則。李睿劬等[9]通過對高超平板邊界層轉(zhuǎn)捩進行氣動光學(xué)診斷，得出平板模型子午中心線的熱流密度。李素循等[10]通過對高速繞流平板邊界層特性的研究，得出不同馬赫數(shù)下平板中心線上的熱流密度分布。然而相比較風(fēng)洞試驗的高成本，轉(zhuǎn)捩準(zhǔn)則可以簡單有效地預(yù)測轉(zhuǎn)捩點的位置，結(jié)合彈體表面的流動狀態(tài)，可以不失真實性地預(yù)測彈體表面的氣動加熱情況。

1 基于無粘流的氣動熱計算方法

求解有限體積[11]的歐拉方程得出彈體表面參數(shù)，對時間采用隱式LU-SGS離散，對流通量項采用二階迎風(fēng)Roe格式離散。

1.1 氣動熱工程公式

火箭彈表面的熱環(huán)境與邊界層外緣的氣動參數(shù)密切相關(guān)。把數(shù)值求解無粘流場獲得的彈體表面參數(shù)作為邊界層外緣參數(shù)，代入工程公式計算火箭彈表面氣動熱。

將火箭彈分為駐點和非駐點區(qū)域進行氣動熱計算，在彈體駐點處采用Fay-Riddell計算熱流密度，在非駐點區(qū)域的層流段和湍流段采用平板熱流密度計算公式，根據(jù)間斷因子Γ[12]，轉(zhuǎn)捩過渡段的表達式為：

(1)

其中，Res為當(dāng)?shù)氐睦字Z數(shù)，Retrl為轉(zhuǎn)捩起始雷諾數(shù)，RetrT為轉(zhuǎn)捩終止雷諾數(shù)。

計算可得當(dāng)?shù)責(zé)崃髅芏龋?/p>

qwc=qwl(1-Γ)+qwtΓ

(2)

其中，qwc轉(zhuǎn)捩過渡段的熱流密度，qwl為層流狀態(tài)下的當(dāng)?shù)責(zé)崃髅芏?，qwt為湍流狀態(tài)下的當(dāng)?shù)責(zé)崃髅芏取?/p>

1.1.1 流線長度的確定

為了準(zhǔn)確地預(yù)測彈體表面每個點上的當(dāng)?shù)乩字Z數(shù)Res，就要準(zhǔn)確計算彈體表面每個格點的流線長度?；谂ｎD的最速下降法[13]，假設(shè)火箭彈表面的流線都是沿著與火箭彈表面形狀相切的方向流動，利用反向追蹤流線直至駐點，可以得出火箭彈表面的流線和流線長度S。

駐點單元是氣流的滯止點，氣動加熱最為嚴(yán)重。根據(jù)數(shù)值求解的計算結(jié)果，若網(wǎng)格單元上四條邊中點的速度方向均指向單元體外側(cè)，那么這個網(wǎng)格就是駐點單元，如圖1所示。

(3)

1.1.2 流線追蹤程序驗證

本文采用NASA-TN5450鈍錐[14]，onera M6機翼[15]和某翼身融合體[15]為驗證飛行器，用流線追蹤程序作出的流線(如圖3～圖5之左圖)與Tecplot作出的表面流線(如圖3～圖5之右圖)進行比較，從而證明流線追蹤程序的正確性。

1.1.3 轉(zhuǎn)捩過渡段的確定

高超聲速飛行器轉(zhuǎn)捩[16]與馬赫數(shù)、壁面溫度、鈍度和表面粗糙度以及防熱涂層的燒蝕氣體引射等多種因素有關(guān)，準(zhǔn)確地預(yù)測高超聲速飛行器邊界層內(nèi)轉(zhuǎn)捩點的位置有相當(dāng)大的難度。

大量實驗數(shù)據(jù)表明，在壁面粗糙度小的情況下，光壁轉(zhuǎn)捩準(zhǔn)則對轉(zhuǎn)捩點的預(yù)測能滿足工程需要。本文按先后順序采用光壁轉(zhuǎn)捩準(zhǔn)則1和轉(zhuǎn)捩準(zhǔn)則2進行經(jīng)驗分析驗證，如式(4)～式(6)所示。其中轉(zhuǎn)捩準(zhǔn)則1是根據(jù)光壁實驗歸納所得，轉(zhuǎn)捩準(zhǔn)則2是根據(jù)有燒蝕的小粗糙度壁面實驗歸納所得。

RetrL～RetrT為轉(zhuǎn)捩過渡段的雷諾數(shù)區(qū)間，將RetrL、RetrT代入式(1)、式(2)可得出轉(zhuǎn)捩過渡段的熱流密度。

RetrL=1.8×105e0.394Mae

(4)

(5)

(6)

1.1.4 轉(zhuǎn)捩準(zhǔn)則確定

本文所選火箭彈可分為彈頭部鈍錐段、圓柱彈身段和舵片組。舵片組平面和圓柱段沿來流方向逆壓梯度小，可近似認(rèn)為是空間平面。為了選取對不同物面形狀綜合吻合能力更好的轉(zhuǎn)捩準(zhǔn)則，在不確定火箭彈轉(zhuǎn)捩點具體位置的情況下，選取來流條件與本文火箭彈來流條件相近的二維平板[17]和三維鈍錐[18]兩個算例進行轉(zhuǎn)捩驗證，相關(guān)模型尺寸見參考文獻。二維平板轉(zhuǎn)捩分析選取了參考文獻中Re/cm為0.27×106，Tw/Tt為0.2時的來流條件。三維鈍錐轉(zhuǎn)捩分析選取了參考文獻中Run20的來流條件。運用轉(zhuǎn)捩準(zhǔn)則所得的熱流密度預(yù)測值與實驗值如圖6、圖7。

從圖6可知關(guān)于平板轉(zhuǎn)捩起止點和表面熱流密度，采用轉(zhuǎn)捩準(zhǔn)則1計算得到的預(yù)測值與實驗值吻合更好。

從圖7可知關(guān)于鈍錐轉(zhuǎn)捩起止點和表面熱流密度，采用兩種轉(zhuǎn)捩準(zhǔn)則計算的預(yù)測值基本相同，與實驗值的誤差相差不大。

由表1可見根據(jù)轉(zhuǎn)捩準(zhǔn)則計算預(yù)測的鈍錐表面轉(zhuǎn)捩起止點值與實驗值。

表1 鈍錐轉(zhuǎn)捩起止點對比

從表1可知：相比鈍錐實驗值，根據(jù)轉(zhuǎn)捩準(zhǔn)則1計算預(yù)測的轉(zhuǎn)捩起始點值優(yōu)于根據(jù)轉(zhuǎn)捩準(zhǔn)則2計算預(yù)測的轉(zhuǎn)捩起始點值；兩種轉(zhuǎn)捩準(zhǔn)則預(yù)測的轉(zhuǎn)捩終止點值一樣，都小于實驗值。

文中選取的兩個算例中，來流Ma都為6，考慮壁面燒蝕的轉(zhuǎn)捩準(zhǔn)則2預(yù)測的轉(zhuǎn)捩起始點都要大于實驗值，而由光壁實驗歸納的轉(zhuǎn)捩準(zhǔn)則1預(yù)測的轉(zhuǎn)捩起始點與實驗值吻合良好。本文火箭彈來流Ma為5，并且火箭彈外形為多形狀的組合，綜合考慮以上因素得：在不確定火箭彈表面轉(zhuǎn)捩起始點的情況下，轉(zhuǎn)捩準(zhǔn)則1能夠預(yù)測火箭彈表面的轉(zhuǎn)捩和氣動加熱情況。

下文取轉(zhuǎn)捩準(zhǔn)則1與完全層流狀態(tài)下預(yù)測的氣動熱進行對比分析。

1.2 高溫防護涂層耦合計算[19]

防熱設(shè)計的首要任務(wù)是針對不同類型的飛行器熱環(huán)境，選擇合理的飛行器防熱方案。本文采用輻射防熱結(jié)構(gòu)方案，表2給出了火箭彈表面3層低熱導(dǎo)率的陶瓷防護涂層的基本參數(shù)。防護涂層按從上到下依次是外層、中間層和內(nèi)層，厚度均為5 mm，a為熱擴散系數(shù)；λ為熱導(dǎo)率。

表2 防護涂層材料基本性能常數(shù)

在防護涂層的導(dǎo)熱過程中，材料的厚度與另外兩個方向的尺度相比小得多，所以將防護涂層內(nèi)的傳熱簡化為沿厚度方向的一維不穩(wěn)定導(dǎo)熱。假設(shè)材料的基本性能參數(shù)為常數(shù)，不計材料受熱線變形和最外層涂層的輻射散熱，對一維導(dǎo)熱方程進行離散，得：

(7)

對式(7)進行有限差分求解溫度分布，對平板的厚度方向和時間進行離散，得到式(8):

(8)

式中：i為空間層上的節(jié)點，n為時間層上的節(jié)點，Δτ為時間步長，取為0.01 s，Δx為空間步長取為5×10-4m。邊界條件采用第二類邊界條件：忽略防護涂層表面的輻射散熱，防護涂層最外層熱流密度q1為氣動加熱熱流密度q1=qw，防護層最內(nèi)層與彈體表面為絕熱，q2=0。

2 火箭彈氣動熱計算研究

2.1 火箭彈數(shù)值計算

本文火箭彈飛行高度為12 km，飛行馬赫數(shù)為5，攻角為0°，數(shù)值計算結(jié)果如圖8、圖9。由于氣體的強烈壓縮滯止，高超聲速彈體頭部形成了一道激波，彈頭駐點區(qū)域和舵片前緣的壓力值遠高于彈身段的壓力值，彈體尾部形成一個低壓死水區(qū)。

2.2 火箭彈最外層起始氣動加熱

火箭彈防護涂層表面的完全層流狀態(tài)下的云圖如圖10的左圖；轉(zhuǎn)捩分段狀態(tài)下的云圖如圖10的右圖。由于轉(zhuǎn)捩湍流段邊界層流體微團的強烈脈動和摻混，轉(zhuǎn)捩分段比完全層流的彈身段氣動加熱嚴(yán)重。

圖11為火箭彈防護涂層表面Z=0截面上，起始狀態(tài)下的完全層流和轉(zhuǎn)捩分段下的熱流密度曲線。

從圖11可以看出無論是完全層流還是轉(zhuǎn)捩分段，彈頭前部的熱流都是從駐點開始急劇下降。當(dāng)熱流下降至150 kw/m2時，在完全層流下彈身段熱流密度趨于平穩(wěn)，穩(wěn)定在100 kw/m2。而考慮轉(zhuǎn)捩湍流之后，由于轉(zhuǎn)捩階段流體微團開始呈現(xiàn)一定的脈動和不穩(wěn)定性，彈頭轉(zhuǎn)捩起始處熱流密度會有一個陡然上升，并且在轉(zhuǎn)捩終止點處熱流密度最大，熱流密度達到680 kw/m2。在湍流段，由于邊界層流體微團形成了有規(guī)律的脈動，沿著彈身方向熱流密度緩慢下降，最后趨于穩(wěn)定，熱流穩(wěn)定在200 kw/m2。

在兩組舵片處，由于氣流滯止，舵片前緣的氣動熱有一個突躍，并且后面一組舵片受到前一組舵片下洗的影響，氣動加熱較前一組舵片嚴(yán)重。并且湍流狀態(tài)下舵片前緣處的熱流密度是完全層流情況下的7～10倍，所以邊界層的流動狀態(tài)對氣動熱的預(yù)測和后續(xù)的氣動熱防護影響較大。

2.3 30 s火箭彈防熱涂層氣動加熱

由于防護涂層表面受到嚴(yán)重的氣動加熱，涂層表面溫度不斷上升。圖12為防熱涂層表面Z=0截面在30 s時完全層流和轉(zhuǎn)捩分段下的溫度曲線。

由圖12可知涂層表面的溫度與氣動熱分布的變化趨勢相同，在30 s時涂層表面駐點的溫度已經(jīng)達到了1 100 K；在完全層流下，舵片前緣的溫度也達到了400 K，在轉(zhuǎn)捩分段下，舵片前緣處于湍流氣動加熱區(qū)，溫度達900～1 000 K，所以對彈體頭部和舵片前緣進行熱防護十分必要。

由于防護涂層受到長時間的氣動加熱，防熱涂層內(nèi)部的溫度也會隨之上升，而合理的彈體表面的溫度對內(nèi)部電子元器件的正常運行起著關(guān)鍵作用，圖13所示為彈體表面Z=0截面上，30 s時完全層流和轉(zhuǎn)捩分段下的溫度曲線。

由圖13可知：彈體表面與涂層表面溫度變化規(guī)律相同，由于防護涂層的低導(dǎo)熱率，30 s時彈體表面溫度都被控制在307 K。表面溫度的分布也證明了3層防護涂層布置是合理的。

3 結(jié)論

1) 高超聲速火箭彈飛行時頭部駐點區(qū)域和舵片前緣氣動加熱最為嚴(yán)重，并且由于后一組舵片受到下洗氣流的影響，氣動加熱比前一組嚴(yán)重，所以需要對彈體頭部和舵片前緣進行相應(yīng)的熱防護。

2) 轉(zhuǎn)捩分段計算得出的氣動熱曲線在彈體頭部的轉(zhuǎn)捩起始點處熱流密度有一個陡升，在轉(zhuǎn)捩終止點后緩慢下降，這與完全層流計算的持續(xù)下降的結(jié)果有區(qū)別，而在實際飛行過程中的確存在轉(zhuǎn)捩過渡段，說明考慮轉(zhuǎn)捩過渡段的氣動熱可更加真實全面地計算彈體表面的氣動熱。

3) 對彈體表面溫度的分析可知，添加防護涂層有效地延緩?fù)繉颖砻鏆鈩訜嵯驈楏w傳遞，本文設(shè)計的3層防護涂層布置方案合理。

[1] 曹偉.高超聲速邊界層的轉(zhuǎn)捩問題[J].空氣動力報,2009,27(5):516-523.

[2] 戚瓊,韓慶.基于Spalart-Allmaras-γ-Reθt轉(zhuǎn)捩模型的橫流不穩(wěn)定性轉(zhuǎn)捩預(yù)測方法[J].氣體物理,2016,1(3):19-24.

[3] 孔維萱,張輝,閻超.適用于高超聲速邊界層的轉(zhuǎn)捩準(zhǔn)則預(yù)測方法[J].導(dǎo)彈與航天運載技術(shù),2013,328(5):54-58.

[4] 朱海濤,單鵬.超聲速平板邊界層旁路轉(zhuǎn)捩直接數(shù)值模擬[J].空氣動力學(xué)學(xué)報,2015,33(3):345-352.

[5] 劉海旭,鄭晨煒,陳林,等.可壓縮平板邊界層轉(zhuǎn)捩的大渦模擬[J].南京航空航天大學(xué)學(xué)報,2014,46(2):246-251.

[6] 陳鵬.高速彈頭氣動熱工程算法與數(shù)值計算研究[D].南京:南京理工大學(xué),2014.

[7] 楊愷.飛行器熱防護系統(tǒng)氣動熱預(yù)測與傳熱分析[D].南京:東南大學(xué),2011.

[8] 高瑩瑩,史增民，李旭東，等.基于Reentry+F的氣動熱工程方法及轉(zhuǎn)捩研究[J].導(dǎo)彈與航天運載技術(shù),2015,339(3):48-52.

[9] 李睿劬，官建，畢志獻，等.高超平板邊界層轉(zhuǎn)捩的氣動光學(xué)診斷技術(shù)[J].空氣動力學(xué)學(xué)報2017，35(1):137-140.

[10] 李素循，師軍，郭孝國.高速繞流平板邊界層特性研究[J].氣體物理2016,1(6):1-4.

[11] BLAZEK J.Computational Fluid Dynamics:Principles and Applications (Second Edition)[M].Elsevier,2005.

[12] 黃志澄編著.航天空氣動力學(xué)[M].北京:中國宇航出版社,1994.

[13] DETERS K.preliminary design estimates of high-speed streamlines on arbitary shaped vehicles Defined by quadrilateral elements[C].11th Applied Aerodynamics conference,1993.

[14] CLEARY J W.Effects of angle of attack and bluntness on laminar heating-rate distributions of a 15 deg cone at a Mach number of 10.6[R].Ames Research Center,1969.

[15] 紀(jì)兵兵，陳金瓶.ANSYS ICEM CFD網(wǎng)格劃分技術(shù)實例詳解[M].北京：中國水利水電出版社,2012.

[16] 王希季.[M].航天器進入與返回技術(shù)上.北京：中國宇航出版社,1991.

[17] CARY JR A M.Turbulent-boundary-layer heat-transfer and transition measurements with surface cooling at Mach 6[R].Langley Research Center,1970.

[18] WIDHOPF G F.Laminar,Transitional,and Turbulent Heat Transfer Measurements on a Yawed Blunt Conical Nosetip[R].Aerospace Corp San Bernardino CA San Bernardino Operations,1972.

[19] 周宇航，張向洪，張敏.高速火箭彈不同防熱方案對氣動加熱影響規(guī)律的研究[J].兵器裝備工程學(xué)報，2016,37(10):24-30.