考慮邊界層轉捩的復雜外形火箭彈氣動熱計算

章 輝，張向洪

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

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

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

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

1.1 氣動熱工程公式

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

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

(1)

其中，Res為當地的雷諾數，Retrl為轉捩起始雷諾數，RetrT為轉捩終止雷諾數。

計算可得當地熱流密度：

qwc=qwl(1-Γ)+qwtΓ

(2)

其中，qwc轉捩過渡段的熱流密度，qwl為層流狀態下的當地熱流密度，qwt為湍流狀態下的當地熱流密度。

1.1.1 流線長度的確定

為了準確地預測彈體表面每個點上的當地雷諾數Res，就要準確計算彈體表面每個格點的流線長度。基于牛頓的最速下降法[13]，假設火箭彈表面的流線都是沿著與火箭彈表面形狀相切的方向流動，利用反向追蹤流線直至駐點，可以得出火箭彈表面的流線和流線長度S。

駐點單元是氣流的滯止點，氣動加熱最為嚴重。根據數值求解的計算結果，若網格單元上四條邊中點的速度方向均指向單元體外側，那么這個網格就是駐點單元，如圖1所示。

(3)

1.1.2 流線追蹤程序驗證

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

1.1.3 轉捩過渡段的確定

高超聲速飛行器轉捩[16]與馬赫數、壁面溫度、鈍度和表面粗糙度以及防熱涂層的燒蝕氣體引射等多種因素有關，準確地預測高超聲速飛行器邊界層內轉捩點的位置有相當大的難度。

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

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

RetrL=1.8×105e0.394Mae

(4)

(5)

(6)

1.1.4 轉捩準則確定

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

從圖6可知關于平板轉捩起止點和表面熱流密度，采用轉捩準則1計算得到的預測值與實驗值吻合更好。

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

由表1可見根據轉捩準則計算預測的鈍錐表面轉捩起止點值與實驗值。

表1 鈍錐轉捩起止點對比

從表1可知：相比鈍錐實驗值，根據轉捩準則1計算預測的轉捩起始點值優于根據轉捩準則2計算預測的轉捩起始點值；兩種轉捩準則預測的轉捩終止點值一樣，都小于實驗值。

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

下文取轉捩準則1與完全層流狀態下預測的氣動熱進行對比分析。

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

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

表2 防護涂層材料基本性能常數

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

(7)

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

(8)

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

2 火箭彈氣動熱計算研究

2.1 火箭彈數值計算

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

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

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

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

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

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

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

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

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

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

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

3 結論

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

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

3) 對彈體表面溫度的分析可知，添加防護涂層有效地延緩涂層表面氣動熱向彈體傳遞，本文設計的3層防護涂層布置方案合理。

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

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

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

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

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

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

[7] 楊愷.飛行器熱防護系統氣動熱預測與傳熱分析[D].南京:東南大學,2011.

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

[9] 李睿劬，官建，畢志獻，等.高超平板邊界層轉捩的氣動光學診斷技術[J].空氣動力學學報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] 黃志澄編著.航天空氣動力學[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] 紀兵兵，陳金瓶.ANSYS ICEM CFD網格劃分技術實例詳解[M].北京：中國水利水電出版社,2012.

[16] 王希季.[M].航天器進入與返回技術上.北京：中國宇航出版社,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] 周宇航，張向洪，張敏.高速火箭彈不同防熱方案對氣動加熱影響規律的研究[J].兵器裝備工程學報，2016,37(10):24-30.