空天飛行器高超聲速氣動熱特性計算方法

王 鵬, 房 帥, 金 鑫, 張衛(wèi)民

(中國航天空氣動力技術(shù)研究院, 北京 100074)

空天飛行器高超聲速氣動熱特性計算方法

王 鵬*, 房 帥, 金 鑫, 張衛(wèi)民

(中國航天空氣動力技術(shù)研究院, 北京 100074)

通過數(shù)值方法和基于無粘表面流線的工程快速方法，計算了空天飛行器基本型在Ma=8.0及Ma=10.2狀態(tài)下的氣動熱特性。數(shù)值方法計算格式選用Roe的FDS格式，工程快速方法中飛行器的表面流線是通過基于直角網(wǎng)格的無粘Euler方程計算得到的，采用參考焓理論沿流線積分即得到沿流線的表面熱流分布。結(jié)果表明，本文建立的氣動熱工程方法及數(shù)值方法得到的機身及機翼的熱流分布與試驗數(shù)據(jù)吻合較好，得到的駐點熱流值與試驗數(shù)據(jù)的誤差小于5%。

表面流線；氣動熱；空天飛行器；數(shù)值計算；工程方法

0 引 言

以往的飛行器設計中，飛行器的熱環(huán)境主要依靠試驗來確定，費用高昂，設計成本高，還存在地面數(shù)據(jù)往飛行狀態(tài)外推的問題，且難于分析流場細節(jié)[1-2]。隨著計算流體力學的發(fā)展，利用CFD技術(shù)結(jié)合風洞試驗來預測氣動加熱是飛行器設計的重要發(fā)展方向。但高超聲速氣動熱的數(shù)值模擬受數(shù)值格式、網(wǎng)格生成、收斂問題[3-6]等各方面因素的影響,多種因素之間的交錯影響決定了數(shù)值計算熱流的復雜性。為了適應高超聲速飛行器概念研究和設計需要,在發(fā)展高精度的CFD技術(shù)的同時，尋求一種快速、簡便、精度較高的氣動熱工程方法具有重要意義。

基于無粘表面流線的氣動熱計算方法是求解高超聲速鈍頭體飛行器表面熱流的有效工程方法。它的計算效率很高，對于簡單的鈍頭體外形，精確度也有保證，在高超聲速鈍頭體飛行器研制初期，發(fā)揮著較大作用。由DeJarnette和Hamiton[7]開發(fā)的基于軸對稱比擬的方法生成表面流線，沿流線積分得到熱流分布。流線求解方法有跟蹤流線法和幾何流線法，但兩種求解方法差別較大，且對于有翼面的復雜外形很難求解。隨后DeJarnette等[8]開發(fā)了UNLAUCH系列程序，逐步將該方法用于較復雜外形。但這類方法的缺陷是流線會間斷，從而影響熱流的連續(xù)求解(作者后續(xù)研究逐漸解決了該問題)。

國內(nèi)研究人員對基于表面流線的方法進行過研究[9]，對于簡單的鈍頭體進行過較為詳細的驗證，但并沒有與數(shù)值計算及風洞試驗進行過系統(tǒng)地比較；且對于具有翼面的較復雜的外形，由于存在流線間斷等原因，未有研究人員進行過系統(tǒng)驗證。

為此，本文建立了一種新的無粘表面流線求解方法，采用基于直角網(wǎng)格的無粘Euler方程計算得到飛行器表面流線，根據(jù)Eckert參考焓理論沿流線積分確定熱流，并與數(shù)值計算與風洞試驗結(jié)果進行了系統(tǒng)地比較。

1 計算模型

本文研究的空天飛行器基本型選自李素循主編的《典型外形高超聲速流動特性》一書[10]，是一種無立尾、機身與三角翼融合的氣動布局。

該空天飛行器外形及熱流測點分布如圖1所示。

位于機身上的熱流測點分布為：熱流測點1-12位于縱向?qū)ΨQ面的上表面，測點13為駐點熱流測點，14-40位于縱向?qū)ΨQ面的下表面。圖1(c)為位于機翼的兩個橫截面上的熱流測點分布。風洞試驗在JF48脈沖型風洞中進行，詳見文獻[10]。

2 工程快速計算理論

2.1基于直角網(wǎng)格的無粘歐拉方程計算理論

計算的流場區(qū)域取為大的長方體, 圖2為計算的直角網(wǎng)格。為了在物面附近加密網(wǎng)格，找出與物面相碰的單元，這里不僅將相碰的單元等分為8個單元， 而且將其鄰居以及鄰居的鄰居等等也劃分為8個,如此形成第2層加密網(wǎng)格。再對第2層網(wǎng)格做上述自適應加密, 形成第3層加密網(wǎng)格。繼續(xù)上述過程, 可以形成近壁面的多層加密網(wǎng)格。

采用Jameson 的有限體積法, 將積分型Euler方程應用到任一單元, 得到

式中，Wk為守恒變量在單元中的平均值, 或格心處的值；Fi是單元表面的通量張量，可取相鄰單元處通量的平均值；Si為表面有效面積；Ωk為單元體積。求和是對單元的m個表面進行。在遠場邊界采用無反射邊界條件。Dk為附加的耗散項， 取二階與四階耗散的融合。

2.2表面流線生成方法

基于上述計算理論，考慮某一時刻t的速度場，在飛行器表面取任一點P0(x,y,z,t)，首先找到與該點最近的物面法向第一層的網(wǎng)格編號C0，將網(wǎng)格C0格心P1存儲的速度V(x,y,z)作為下游流線上P1點的切線方向，并取該網(wǎng)格步長dh作為流線的積分步長。在距離P1(x,y,z,t)點V(x,y,z)*dh范圍內(nèi)尋找一與C0鄰近網(wǎng)格格心點P2(x,y,z,t)，若P2(x,y,z,t)與P1(x,y,z,t)坐標向量與P1點的速度矢量形成的夾角最小且不大于60°(若無滿足條件的點則作為流線的終點)，則將P2點作為下游流線的第二點，以此類推求出下游流線上的各點，得到一條折線，記為經(jīng)過飛行器表面P0點的流線；P0點上游流線求法類似。圖3為利用本文建立的方法得到的空天飛行器的表面無粘流線示意圖。

2.3熱流計算方法

2.3.1 駐點熱流計算方法

工程設計中駐點熱流計算理論較多，主要有Fay-Riddell、Kemp-Riddell及Lees等公式[11-12]。Fay-Riddell理論利用相似性假定和多羅得尼津-曼格勒變換，將高溫氣體邊界層偏微分方程轉(zhuǎn)化為常微分方程，將氣體熱力學特性等有關(guān)的無因次參數(shù)假設為一系列常數(shù)，利用薩瑟蘭定律，得出的駐點熱流經(jīng)驗公式。本文駐點熱流計算采用的就是Fay-Riddell公式。平衡邊界層條件下的駐點熱流計算公式為：

(2)

其中，hD為空氣平均離解焓，Pr為普朗特數(shù)，Le為路易斯數(shù)，h∞為來流焓值，hs為駐點焓值，μs為駐點粘性系數(shù)，ρs為駐點密度，ζ為壁面焓值與駐點焓值之比，(due/dx)s為駐點的速度梯度。

2.3.2 沿流線的熱流計算方法

本文采用的表面熱流計算理論為Zoby[13]發(fā)展的基于參考焓理論的層流及湍流計算理論。層流氣動熱計算公式為：

qwl= 0.22(Reθe)-1(ρ*/ρe)(μ*/μe)·

其中，qwl為所求流線上任一點的熱流，ρ*、μ*為參考密度及粘度，ρe、μe為邊界層邊緣的密度及粘度，Haw為絕熱壁面焓，Hw為壁面焓，Reθe為根據(jù)邊界層外緣參數(shù)計算得到的層流動量厚度雷諾數(shù)，其中動量厚度根據(jù)下式確定：

(4)

其中，ue為所求流線上任一點對應的邊界層邊緣的速度，V∞為自由來流速度，r為尺度因子(拉梅系數(shù))，s為沿流線的積分長度。

湍流氣動熱的計算公式為

qw t=c1(Reθe)-m(ρ*/ρe)(μ*/μe)m·

其中，湍流動量厚度雷諾數(shù)的湍流厚度根據(jù)下式得到：

(6)

其中參數(shù)m及c1、c2、c3、c4的定義參見文獻[13]。

對于尺度因子的求解,假設邊界層內(nèi)的流動方向與物面繞流方向一致，且物面上垂直于流線的流動速度相對于流線方向的速度可以略去，將三維邊界層方程簡化為軸對稱的邊界層方程，方程中流線長度作為子午線的長度，垂直于流線的尺度因子即為式中的r，具體求解參見文獻[13]。

3 數(shù)值計算方法

3.1計算理論

不計質(zhì)量力情況下，三維直角坐標系中非定?？蓧嚎sN-S方程的守恒形式可以寫為下列向量式：

其中，U為守恒變量矢量，F(xiàn)、G、H為對流項，F(xiàn)v、Gv、Hv為粘性項，各項定義具體參見文獻[14]。

本文CFD計算格式采用Roe的FDS格式，因其具有較高的粘性分辨率，且收斂性較好。物理流動模型分別為無湍流模型的薄層假設與B-L湍流模型。時間離散采用LU-SGS方法。

3.2計算網(wǎng)格及邊界條件

采用的計算網(wǎng)格為多塊結(jié)構(gòu)網(wǎng)格，如圖4所示，并保證第一層網(wǎng)格的壁面網(wǎng)格雷諾數(shù)小于10。網(wǎng)格數(shù)量約為350萬。

采用的邊界條件包括：遠場采用外插邊界條件；物面采用無滑移的等溫壁面條件(Tw=300K,與試驗條件一致)；對稱面采用鏡面反射邊界條件。

4 結(jié)果分析

利用建立的數(shù)值及工程方法計算了空天飛行器在Ma=8.0、10.2，迎角為0°條件下的熱流分布，并與試驗數(shù)據(jù)對比。具體計算條件如表1所示。

表1 計算條件Table 1 Computational conditions

熱流采用無量綱值，參考熱流采用書中給定的熱流試驗值：Ma=8.0時，qref=476.30kW/m2；Ma=10.2時，qref=585.50kW/m2。Ma=8.0計算得到的縱向?qū)ΨQ面上、下表面的熱流如圖5及圖6所示。

圖5為通過工程方法、數(shù)值方法得到的機身縱向?qū)ΨQ面下表面熱流分布與試驗數(shù)據(jù)的對比(圖中l(wèi)aminar-，表示層流氣動熱的工程計算；turbulent-，表示湍流氣動熱的工程計算；NS-laminar，表示層流模型數(shù)值計算；NS-bl，表示bl湍流模型數(shù)值計算)。迎角0°時，下表面是背風面，采用工程方法及數(shù)值方法，都可以揭示熱流的變化趨勢，在兩種方法中采用湍流模型的計算結(jié)果明顯高于無湍流模型的結(jié)果。 因為實驗數(shù)據(jù)是在層流狀態(tài)下的氣動熱，層流的計算結(jié)果與實驗數(shù)據(jù)比較吻合，且數(shù)值計算的結(jié)果比工程方法更接近于實驗值。圖6為縱向?qū)ΨQ面上表面熱流分布，數(shù)值與工程方法及實驗數(shù)據(jù)在后半段吻合較好，前半段差別較大。層流的數(shù)值計算結(jié)果明顯優(yōu)于工程方法，且數(shù)值方法與實驗數(shù)據(jù)吻合較好。Ma=10.2得到的縱向?qū)ΨQ面熱流分布如圖7所示。

工程方法或數(shù)值方法都可以揭示熱流的變化趨勢，數(shù)值方法中采用層流氣動熱的計算方法與風洞試驗試驗數(shù)據(jù)吻合的較好。工程方法與數(shù)值方法及實驗數(shù)據(jù)存在少許差異。圖7為縱向?qū)ΨQ面上表面熱流分布的對比，數(shù)值方法及工程方法與實驗數(shù)據(jù)吻合較好，在拐角處差異較大。Ma=8.0得到的機翼II-II截面的熱流如圖8所示。

通過兩種方法得到的計算曲線規(guī)律一致，且基本重合，但在機翼前緣附近數(shù)值計算的結(jié)果高于工程計算得到數(shù)值，在存在試驗數(shù)據(jù)的一段內(nèi)與試驗數(shù)據(jù)吻合較好。表2為通過數(shù)值方法與工程方法得到的駐點熱流值與試驗值的比較，兩種方法得到的駐點熱流與試驗值之間的誤差均小于5%。

表2 駐點熱流對比Table 2 Aerodynamic heating comparisons for stagnation point

5 結(jié) 論

針對空天飛行器基本形，建立了高超聲速氣動熱的數(shù)值方法及基于表面流線的氣動熱快速工程方法，表面流線通過基于直角網(wǎng)格的無粘Euler方程得到，得到的結(jié)論為：

建立的數(shù)值方法及快速工程方法都較好的反映了空天飛行器在高超聲速時的氣動熱分布規(guī)律。通過與試驗數(shù)據(jù)的對比，兩種方法在層流狀態(tài)下計算得到的機身與機翼的熱流分布與試驗數(shù)據(jù)吻合較好；通過數(shù)值計算及工程計算得到的駐點熱流與試驗值的誤差在5%以內(nèi)。

本文建立的基于表面流線的氣動熱工程方法計算速度快，計算結(jié)果較好，雖然與數(shù)值計算方法及實驗數(shù)據(jù)相比，在個別點存在較大的差異，在下一步工作中力求對其進行較為系統(tǒng)的驗證，使其成為氣動熱快速計算的工具。

[1]范緒箕. 氣動熱與熱防護系統(tǒng)[M]. 科學出版社. 2004.

[2]Liu Xin, Deng Xiaogang. Application of high-order accurate algorithm to hypersonic viscous flows for calculating heat transfer distributions[R]. AIAA-2007-691, 2007.

[3]Klopfer G H, Yee H C. Viscous hypersonic shock-on-shock interaction on blunt cowl lips[R]. AIAA-88-0233, 1988.

[4]周偉江, 姜貴慶. 迎風TVD 格式在粘流計算中的應用研究與改進[J]. 計算物理, 1999, 16(4): 401-408.

[5]潘沙, 馮定華, 丁國昊, 等. 氣動熱數(shù)值模擬中的網(wǎng)格相關(guān)性及收斂[J]. 航空學報, 2010, 31(3): 493-499

[6]Cohen C B, Reshotko E. Similar solution for the compressible laminar boundary layer with heat transfer and pressure gradient[R]. NACA Report 1293, 1956.

[7]DeJarnette F R, Jones M H. Calculation of inviscid surface streamlines and heat transfer on shuttle type configurations. Part 1 Description of basic method[R]. NASA CR-111921, 1971.

[8]Dejarnette F R, Hamilton H H, Weilmuenster K J. New method for computing convective heating in stagnation region of hypersonic vehicles[R]. AIAA-2008-1261, 2008.

[9]賀國宏, 高曉斌, 龐勇. 高超聲速再入體表面熱流數(shù)值模擬研究[J]. 空氣動力學學報, 2001, 19(2): 177-185.

[10]李素循. 典型外形高超聲速流動特性[M]. 北京: 國防工業(yè)出版社, 2007.

[11]Fay J A, Riddell F R. Theory of stagnation point heat transfer in dissociated air[J]. Journal of the Aeronautical Science, 1958, 25(2): 112-120.

[12]Kemp N H, Rose R H, Detra R W. Laminar heat transfer around blunt bodies in dissociated air[J]. Journal of the Aerospace science, 1959, 26(7): 68-75.

[13]Zoby E V, Moss J N, Sutton K. Approximate convective heating equations for hypersonic flows[J]. Journal of Spacecraft and Rockets, 1981, 18(1): 64-70.

[14]傅德薰, 馬延文, 等. 計算流體力學[M]. 北京: 高等教育出版社, 2002.

Computationalmethodonhypersonicaerodynamicheatingforaerospacevehicle

WANG Peng*, FANG Shuai, JIN Xin, ZHANG Weimin

(ChinaAcademyofAerospaceAerodynamics,Beijing100074,China)

Hypersonic aerodynamic heating characteristics for an aerospace vehicle atMa=8.0 andMa=10.2 was calculated by computational fluid dynamics (CFD) method and engineering rapid computational method based on in-viscid surface streamlines. The Roe flux difference scheme(FDS) was employed in the CFD method. For the engineering rapid computational method, the in-viscid surface streamlines were determined by the Euler equations on the Cartesian grid, and the surface heating distribution can be computed by the integration along the in-viscid surface streamlines according to the Eckert reference enthalpy theory. Both laminar and turbulent flows were simulated in the present study. Comparisons were carried out between the engineering rapid computational results, the CFD results, and the experimental results. The results show that the heating distribution can be well predicated for the airframe and the wing of the aerospace vehicle by the CFD method and the rapid engineering method developed in this paper. The error of the aerodynamic heating distribution at the stagnation point is less than 5%.

surface streamline; aerodynamic heating; aerospace vehicle; numerical simulation; engineering method

V211.3

A

10.7638/kqdlxxb-2015.0067

0258-1825(2017)05-0640-05

2015-05-24;

2015-11-26

國家自然科學基金(61273153)

王鵬*(1984-)，男，山東濰坊人，工學碩士，工程師，研究方向：高超聲速氣動力/熱的數(shù)值計算和工程估算及嵌入式大氣數(shù)據(jù)傳感系統(tǒng)的研究. E-mail：pengwang0413@163.com

王鵬, 房帥, 金鑫, 等. 空天飛行器高超聲速氣動熱特性計算方法[J]. 空氣動力學學報, 2017, 35(5): 640-644.

10.7638/kqdlxxb-2015.0067 WANG P, FANG S, JIN X, et al. Computational method on hypersonic aerodynamic heating for aerospace vehicle[J]. Acta Aerodynamica Sinica, 2017, 35(5): 640-644.