高超聲速飛行器邊界層外緣參數仿真分析*

孟竹喧，胡　凡,2，彭　科，張為華,2

(1.國防科技大學 航天科學與工程學院， 湖南 長沙　410073；2.高超聲速沖壓發動機技術重點實驗室， 湖南 長沙　410073)

?

高超聲速飛行器邊界層外緣參數仿真分析*

孟竹喧1，胡凡1,2，彭科1，張為華1,2

(1.國防科技大學 航天科學與工程學院， 湖南 長沙410073；2.高超聲速沖壓發動機技術重點實驗室， 湖南 長沙410073)

摘要：以高超聲速飛行器為研究對象，構建快速準確計算高超聲速飛行器無黏邊界層外緣參數的計算方法。擬合空氣比熱、比熱比隨溫度變化曲線，建立空氣屬性溫度劃分準則。基于不同空氣屬性建立高超聲速飛行器邊界層外緣參數工程與數值計算模型，采用鈍雙錐模型，對比分析工程估算、無黏數值及有黏數值計算方法的計算結果。結果表明，0°攻角狀態下，基于無黏流場的數值計算與工程估算和有黏數值計算的壓強最大差值分別為1.19%和2.39%；10°攻角狀態下，最大差值分別為5%和50%；從而證明所提出的無黏數值計算方法明顯優于工程計算方法，為進一步快速準確計算高超聲速飛行器氣動熱環境奠定了重要基礎。

關鍵詞：高超聲速飛行器；比熱/比熱比；空氣屬性；無黏流場仿真；邊界層外緣參數

高超聲速飛行器氣動加熱物理過程本質為邊界層內的復雜傳熱傳質過程，確定邊界層外緣參數是實現高超聲速飛行器氣動熱環境快速、準確計算的重要基礎。根據普朗特邊界層理論，大雷諾數條件下流場邊界層很薄，邊界層外的流場參數與飛行器無黏外流場參數基本一致，目前發展了兩類邊界層外緣參數計算方法，一是基于牛頓理論的工程計算方法，二是基于無黏假設的流場數值模擬方法，前者的計算效率優于后者，后者的計算精度則更高。

目前國內外研究主要集中在單純數值模擬或工程計算，以流場特性結合工程與數值方法的研究并不成熟，難以滿足兼顧計算準確性與效率的要求。

高超聲速飛行器外流場不同區域溫度分布差異較大，空氣屬性包含量熱完全氣體、完全氣體、平衡氣體、非平衡氣體四類，目前常采用的量熱完全氣體假設(即取比熱/比熱比為常數)在高溫環境下失效。根據溫度條件確定空氣屬性，準確構建比熱/比熱比計算模型，進而合理建立空氣物性參數計算模型，這對保證邊界層外緣參數計算精度至關重要。

本文研究高超聲速飛行器邊界層外緣參數計算問題。根據溫度條件確定空氣屬性，構建比熱/比熱比計算模型，建立無黏邊界層外緣參數工程計算與數值計算模型，完成算例分析，定量對比二者精度、效率等，研究快速準確確定高速飛行器邊界層外緣參數的計算方法，為氣動熱環境計算提供重要的參數支持。

1邊界層外緣參數計算模型

1.1空氣屬性分類與比熱/比熱比計算模型

Anderson對氣體模型劃分的定義為[1-3]：完全氣體包括量熱完全氣體、熱完全氣體和化學反應完全氣體混合物；化學反應完全氣體混合物中每一組分都遵從完全氣體狀態方程，達到化學反應平衡的完全氣體稱為平衡化學反應完全氣體。本文不考慮非平衡氣體。

大氣主要成分是N2，O2，0.1 MPa下空氣中反應發生溫度如圖1所示。恒壓狀態條件下，800 K分子振動激發，此時定比熱假設失效，應用熱完全氣體假設模型計算；2500 K時O2開始離解，這時應使用平衡化學反應完全氣體模型。氣體模型適用溫度范圍如表1所示。

圖1　0.1 MPa條件下空氣分子振動激發和O2，N2的離解電離Fig.1　Vibrational excitation of air molecularand the dissociative ionization of O2,N2 at 0.1 MPa

氣體模型完全氣體量熱完全氣體熱完全氣體平衡氣體溫度范圍800K以下800～2500K2500K以上

定容比熱Cv、定壓比熱Cp和比熱比γ是熱力參數計算的重要基礎，研究中普遍采用熱完全氣體假設計算邊界層外緣參數，此假設下比熱比均為常數，而化學反應平衡氣體中比熱比不僅是溫度的函數，還與每種組分的分壓有關[4]，此時直接求解比熱比較為困難，文獻[4]的附錄I中給出的空氣熱力參數插值表，對定容比熱Cv、定壓比熱Cp和比熱比γ作多項式擬合，各擬合系數如表2所示，Cv，Cp單位均為J/(kg·K)。

圖2～4給出了三個熱力參數隨溫度的變化。由圖可以看出，擬合得到的Cv，Cp和γ與文獻中給出的數值吻合較好，且變化趨勢與空氣隨溫度變化特性十分一致，能夠很好地滿足熱力參數計算需要。

表2　比熱/比熱比擬合多項式系數

圖2　Cv值比較Fig.2　Comparison of Cv

圖3　Cp值比較Fig.3　Comparison of Cp

圖4　Gamma值比較Fig.4　Comparison of gamma

Cv擬合值與文獻[4]給定值最大相差為0.085 6%，Cp最大差值為0.068 3%，γ最大差值為0.060 7%。通過分析可知，本文對氣體模型適用溫度范圍劃分與實際吻合良好，各溫度范圍適于作為氣體模型選擇標準。

1.2無黏邊界層外緣參數計算模型

1.2.1工程計算模型

文獻[5]將高超聲速流動分為低高超聲速區(3.0≤Ma∞≤6.5)和高高超聲速區(Ma∞≤8.5)，給出了高超聲速飛行器壓強系數計算方法選擇原則，如表3所示。

溫度高于800 K時，氣體特性不再單純表征為量熱完全氣體，計算中普遍采用的正激波前后關系式和等熵關系式不再適用，本文給出不同溫度狀態下邊界層外緣參數工程計算方法，計算流程如圖5所示。

對高超聲速流場正激波后參數，800 K以下的量熱完全氣體采用正激波前后關系式可以求解；而對非量熱完全氣體，Pv=RT仍然成立，但是R在此時是壓強P與溫度T的函數，通過簡單的求解不能求出正激波后參數，需要迭代方法求解。從三大守恒定律和等熵關系出發，將控制方程重新組合后，應用文獻[6]中的方法通過一系列迭代求解，得到邊界層外緣各參數。

表3高超聲速壓強系數計算方法

Tab.3Calculation method of hypersonic pressure coefficient

區 域低高超聲速高高超聲速迎風面背風面迎風面背風面頭部圓Dahlem-buckDahlem-buck修正牛頓Prandtl-Meyer平切楔Dahlem-buck修正牛頓Prandtl-Meyer身部圓有攻角錐ACM修正牛頓Prandtl-Meyer平切錐ACM修正牛頓Prandtl-Meyer升力面邊條切錐Dahlem-buck修正牛頓Prandtl-Meyer普通切楔Dahlem-buck修正牛頓Prandtl-Meyer

圖5　邊界層外緣參數工程計算方法流程圖Fig.5　Engineering calculation flow-process diagram of outer edge boundary parameters

對于完全氣體(包括800 K以下的量熱完全氣體與800～2500 K的熱完全氣體兩類)，焓與內能都只是溫度的函數，同時有狀態方程Pv=RT成立，且R是常數，這時正激波后壓強P2和密度ρ2，利用文獻[7]中給出的正激波前后等熵關系式得到邊界層外緣密度ρe，由完全氣體狀態方程求解邊界層外緣焓He，由能量方程得到邊界層外緣速度Ve，邊界層外緣黏性系數μe由薩瑟蘭(Sutherland)公式求得。

對實際高超聲速流動環境中高于2500 K的平衡氣體，由于其復雜的特征，單純的工程計算很難確切地得到理想結果，所以為了提高熱流密度的計算精度，由蘇聯高溫空氣動力學函數表擬合的誤差不大于10%的流場特性參數計算公式[8]求解得到邊界層外緣參數。

1.2.2數值計算模型

邊界層外緣參數計算基于無黏歐拉方程如式(1)～(3)所示，采用有限體積法求解；空間離散采用二階迎風總變差非增(Total Variation Diminishing, TVD)格式，時間推進采用二階隱式格式[9]；為提高計算效率，采用非結構網格劃分流場；計算采用平衡氣體模型，對于非熱完全氣體狀態下流場同樣適用[10]。

連續性方程：

(1)

動量方程：

(2)

能量方程：

(3)

其中：ρ是流場當地密度；u,v,ω分別是當地速度在x,y,z方向的分量；P是當地壓強；μ是黏性系數；h是當地焓值。

2邊界層外緣參數計算與分析

對于處于平衡狀態下的完全氣體來說，所有狀態參數(例如ρe,ue,Te,SL,Se,Pe)中只有兩個是獨立的，其他參數均可由熱力學關系通過這兩個獨立變量確定。首先求解邊界層外緣壓強，以此作為獨立變量進一步求解邊界層外緣其他五個參數，本節僅給出壓強計算結果分布云圖。

2.1鈍雙錐模型邊界層外緣參數計算

鈍雙錐模型是國內外氣動熱計算相關文獻中普遍使用的計算模型，且這一模型在風洞實驗中較為常見，擁有詳盡的氣動熱實驗數據，本文使用這一模型為后續氣動熱計算模型驗證奠定了基礎。鈍雙錐模型尺寸為：頭錐曲率半徑為3.835 mm，前錐半錐角為12.84°，后錐半錐角為7°，前錐距頭部距離為69.55 mm，后錐距頭部距離為122.24 mm[11]。

本文計算了0°，10°攻角狀態下，來流P∞=59.92 Pa，T∞=48.88 K，Ma∞=9.86的鈍雙錐模型邊界層外緣無黏流場壓強、速度、溫度、熵與密度以及流線長度。來流方向為X軸負方向，攻角是在XZ平面內與X軸正方向的夾角。

應用邊界層外緣無黏流場工程與數值計算模型對鈍雙錐模型進行仿真計算，與有黏數值方法計算的迎風面母線壓強計算結果進行對比分析。

如圖6所示，0°攻角條件下，無黏數值計算得到壓強最大值位于頭部頂點處，為7.65 kPa；如圖7所示，10°攻角條件下，壓強最大值位置較0°攻角條件下稍向后移，為7.657 kPa。同樣地，由云圖可以明顯看出，10°攻角條件下，迎風面壓強較0°攻角條件下略有增大。

圖6　0°攻角無黏數值壓強分布云圖Fig.6　0° inviscid numerical pressure

圖7　10°攻角無黏數值壓強分布云圖Fig.7　10° inviscid numerical pressure

圖8　0°攻角有黏數值壓強分布云圖Fig.8　0° viscid numerical pressure

圖9　10°攻角有黏數值壓強分布云圖Fig.9　10° viscid numerical pressure

如圖8所示，0°攻角條件下，有黏數值計算壓強最大值位于頭部頂點處，為7.734 kPa；如圖9所示，10°攻角條件下，壓強最大值位于頭部頂點稍向后移，為7.620 kPa。相較之前工程計算方法與無黏數值計算方法一樣，10°攻角條件下，迎風面壓強較0°攻角條件下略有增大。

圖10和圖11分別給出模型處于0°與10°攻角條件下的邊界層外緣無黏流場采用工程計算方法得到的壓強分布云圖。0°攻角條件下，壓強最大值位于頭部頂點處，為8.218 kPa；10°攻角條件下，壓強最大值仍位于頭部頂點處，為8.218 kPa。由云圖可以明顯看出，10°攻角條件下，迎風面壓強較0°攻角條件下略有增大。這與文獻中迎風面攻角隨攻角增大而增大且迎風面與背風面壓強差異隨攻角增大而增大的描述吻合。

圖10　0°攻角無黏工程壓強分布云圖Fig.10　0°inviscid engineering pressure

圖11　10°攻角無黏工程壓強分布云Fig.11　10°inviscid engineering pressure

由上文建模方法計算，可由工程與數值方法分別得出0°和10°攻角狀態下邊界層外緣無黏流場溫度、速度、密度、熵和流線長度的計算結果，在此不一一列出。

2.2計算結果分析

由于本文計算模型缺少實際邊界層外緣參數試驗數據，本文將考慮黏性的流場數值計算結果作為對比標準，將無黏數值、有黏數值和工程計算結果進行對比分析。

取鈍雙錐母線上的點，取其距原點長度為X軸變量，以坐標點對應的壓強值為Y軸變量作曲線如圖12和圖13所示。0°攻角狀態下工程計算結果和無黏數值計算結果都與有黏數值計算結果吻合較好，工程方法與有黏數值方法結果最大差值為2.39%，無黏數值方法與有黏數值方法最大差值為1.19%。10°攻角狀態下，無黏數值計算結果與有黏數值計算結果最大誤差不超過5%，而工程計算結果最大誤差達到50%。無黏數值方法計算精度優于工程方法。

對比曲線圖12和圖13可知，無黏數值計算結果整體稍小于有黏數值計算結果，且在飛行器后部差異大于頭部。分析誤差產生的原因如下：

1)在近壁面區域，考慮黏性的流場相較于無黏流場邊界層外緣向外流場稍有擴張，使得空氣壓縮效應更加劇烈，故無黏數值方法求解的邊界層外緣壓強較小；

2)有黏流場中流體黏性作用使得飛行器后部流體流速較無黏流場低，由此壓強會較無黏流場稍小，離頭部位置越遠，流體減慢程度越大，則有黏流場與無黏流場計算結果差異越大。

無黏數值計算時間周期以小時為單位，略大于工程計算周期，而有黏流場數值計算周期以天為單位，無黏數值計算效率明顯優于有黏數值計算。故在求解飛行器氣動熱環境計算中，采用無黏數值計算方法求解邊界層外緣參數，可以兼顧計算精度與效率。

圖12　0°攻角迎風面母線壓強對比Fig.12　Comparison of pressure on generatrix of windward at 0°

圖13　10°攻角迎風面母線壓強對比Fig.13　Comparison of pressure on generatrix of windward at 10°

綜上分析，所提出的高超聲速邊界層外緣參數無黏數值求解方法適用于高超聲速條件下大雷諾數流場，但在飛行器后部區域稍有誤差，計算模型還需通過飛行試驗進一步修正。

3結論

1)建立了比熱/比熱比計算模型，提出了空氣屬性隨溫度變化劃分準則，為不同溫度條件下選擇適當無黏邊界層外緣參數計算模型奠定了基礎；

2)基于溫度空氣屬性劃分準則，建立了邊界層外緣參數的工程與無黏數值計算模型，以鈍雙錐模型為對象，求解了其邊界層外緣壓強、溫度、速度、度、熵和流線長度；

3)將鈍雙錐模型壓強無黏數值計算和工程計算結果與有黏計算結果進行對比，無黏數值計算結果的準確性明顯優于工程計算結果，計算效率明顯優于有黏數值計算，采用無黏數值方法計算邊界層外緣參數能夠兼顧計算精度與效率，為氣動熱環境分析提供重要借鑒。

參考文獻(References)

[1]Anderson J D. Hypersonic and high temperature gas dynamics[M]. New York: McGraw-Hill Book Company, 1989.

[2]張志成. 高超聲速氣動熱和熱防護[M]. 北京: 國防工業出版社, 2003: 1-5.

ZHANG Zhicheng. Hypersonic aerodynamic heat and thermal protection[M]. Beijing: National Defence Industry Press, 2003: 1-5. (in Chinese)

[3]張魯民. 再入飛船返回艙空氣動力學[M]. 北京: 國防工業出版社, 2002.

ZHANG Lumin. Aerodynamics of reentry capsule[M]. Beijing: National Defence Industry Press, 2002. (in Chinese)

[4]徐華舫. 空氣動力學基礎[M]. 北京: 北京航空學院出版社, 1986.

XU Huafang. Fundamentals of aerodynamics[M]. Beijing: Beijing Aviation Institute Press, 1986. (in Chinese)

[5]Moore M,Williams J. Aerodynamic prediction rationale for analyses of hypersonic configuration[J]. AIAA, 1989.

[6]蔣友娣. 高超聲速飛行器氣動熱和表面瞬態溫度計算研究[D]. 上海: 上海交通大學, 2008: 1-4.

JIANG Youdi. Calculation of aerodynamic heating and transient surface temperature for hypersonic aircraft[D]. Shanghai: Shanghai Jiaotong University, 2008: 1-4. (in Chinese)

[7]瞿章華, 曾明, 劉偉, 等. 高超聲速空氣動力學[M]. 長沙: 國防科技大學出版社, 2001: 1-4.

QU Zhanghua, ZENG Ming, LIU Wei, et al. Hypersonic aerodynamics[M]. Changsha: National University of Defense Technology Press, 2001: 1-4. (in Chinese)

[8]耿艷棟， 肖建軍. 關于空天一體化的初步研究[J]. 裝備指揮技術學院學報, 2004, 15(6): 19-52.

GENG Yandong, XIAO Jianjun. Preliminary study on air and space integration[J] Journal of the Academy of Equipment Command & Technology, 2004, 15(6): 19-52. (in Chinese)

[9]董維中, 高鐵鎖. 帶座艙飛船高超聲速再入非平衡流場的數值研究[J]. 空氣動力學學報, 2002, 20(2):239-245.

DONG Weizhong, GAO Tiesuo. Numerical studies of non-equilibrium flow-fields over a capsule-type hypersonic reentry vehicle[J]. Acta Aerodynamica Sinica, 2002, 20(2): 239-245. (in Chinese)

[10]黃志澄. 高超聲速飛行器空氣動力學[M]. 北京: 國防工業出版社, 1995.

HUANG Zhicheng. Hypersonic aircraft aerodynamics[M]. Beijing: National Defence Industry Press, 1995. (in Chinese)

[11]方磊. 基于流線跟蹤法的氣動熱工程計算研究[D]. 南京: 南京航空航天大學, 2008.

FANG Lei. The study of aerodynamic heating predictions for hypersonic aircrafts based upon tracking the surface streamtrace[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2008. (in Chinese)

doi:10.11887/j.cn.201602006

*收稿日期：2015-02-17

基金項目：國家自然科學基金資助項目(51406230)

作者簡介：孟竹喧(1990—)，女，吉林白山人，博士研究生，E-mail：m13687361976@163.com；張為華(通信作者)，男，教授，博士，博士生導師，E-mail：zwh_kjs@163.com

中圖分類號:V411

文獻標志碼:A

文章編號：1001-2486(2016)02-031-06

Simulation analysis of outer edge boundary parameters for hypersonic-glide vehicle

MENG Zhuxuan1, HU Fan1,2, PENG Ke1, ZHANG Weihua1,2

(1.College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China；2.Science and Technology on Scramjet Laboratory, Changsha 410073, China)

Abstract：Method of inviscid flow field simulation of the outer edge boundary layer for hypersonic vehicle was created, which was faster than the numerical method and more accurate than the engineering method. And the curve of specific heat and the specific heat ratio changed with temperature were matched, which deduced the law of air attribute on different temperatures. Based on the law, the model of parameters of the outer edge boundary layer for hypersonic vehicle was built, and the comparisons of the results of the proposed method and the engineering method as well as the viscid method of blunt double-cone model were analyzed. Result shows that at the 0° angle of attack, compared with the viscid method, the maximum differences of the pressures between engineering method and viscid method are about 1.19% and 2.39% respectively, while at the 10°angle of attack, they are about 5% and 50% respectively. And the proposed inviscid numerical method obtained higher accuracy is superior to the engineering method, which lays the foundation for the thermal environment calculation of hypersonic-glide vehicle.

Key words：hypersonic vehicles;specific heat and specific heat ratio;air attribute;inviscid flow field simulation;outer edge boundary parameters

http://journal.nudt.edu.cn