基于高效數值方法的高速飛行器氣動力熱特性快速預測研究

王 榮，張學軍，紀楚群

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

基于高效數值方法的高速飛行器氣動力熱特性快速預測研究

王 榮*，張學軍，紀楚群

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

結合空間推進數值模擬方法和流線追蹤法發展了氣動力、熱快速預測技術。針對高速飛行器的算例研究表明，相對常規時間推進方法，基于空間推進法的氣動快速預測方法計算效率提高了一個量級，而兩者氣動力計算精度相當，相對實驗熱流預測誤差在20%以內。所發展的技術為適應氣動外形快速選型和優化設計需求提供了有效的方法。

空間推進;流線追蹤;高速飛行器

0 引言

飛行器設計面向的氣動外形越來越復雜，對氣動力和氣動熱性能預測評估技術也提出了更高的要求，一方面要求預測精度更好，另一方面要求預測評估時間愈短愈好。尤其是在選型設計階段和氣動外形優化設計方面更迫切地需要建立氣動力熱特性的快速預測評估技術。

精度和效率往往是有沖突的，目前氣動力熱預測廣泛采用的方法有工程方法和基于時間相關法(時間推進法)的CFD數值方法。純工程計算方法計算量小，可以在短時間內完成大量的計算，但是計算精度不夠理想。CFD方法計算精度與工程方法相比具有非常明顯的優勢，但缺點是計算量非常大，計算成本較高。

對于主流方向無分離的超聲速流動，可采用空間推進方法沿主流方向高效推進求解。相對純工程方法和時間推進方法，該方法相對高效同時能夠保持精度。因此，近些年國內外針對拋物化NS方程(PNS)采用空間推進方法的計算工作很多［1-5］。由于PNS需要計算粘性項，對網格要求較高，邊界層網格密度大。根據邊界層理論［6］，粘性主要作用于薄邊界層，流場外層可看作是無粘流動。因此可先采用無粘歐拉方程通過空間推進法快速求解外層流場，以此得到的邊界層外緣參數為基礎，再采用流線追蹤工程方法快速求出熱流。這種將無粘高效數值方法和工程方法相結合預測氣動力熱性能的方法，相對PNS方法無需求解粘性項可減少一定的計算量，同時也降低了對網格的要求，相對于采用時間推進法［6-8］獲得邊界層外緣無粘流場參數的做法耗時少，更加快速高效。

對于超聲速流動，擾動僅向下游傳播，此時定常Euler方程為流向初值適定［1］的擬線性雙曲型方程，在給定初始橫流面流場解條件下，可通過空間推進求解，且推進一次就可得到流場解。初值面可通過時間相關法求出，由于求解初值面僅需很少的網格，耗時是比較少的。因此總的來講，基于定常Euler方程空間推進求解方法計算效率高，魯棒性好，在一定馬赫數和迎角范圍可滿足工程型號復雜外形氣動特性快速評估要求。

沿流線氣動加熱積分計算的方法是由Zoby［9］于上世紀80年代初提出的，Zoby基于對軸對稱邊界層動量方程的分析并結合Eckert參考焓概念給出了可壓縮條件下物面摩阻與邊界層動量厚度之間的積分關系式，而后運用“修正雷諾比擬”給出適用于變熵條件的層流、湍流氣動加熱率計算的積分關系式;這些關系式結合“軸對稱比擬法”技術可用于一般三維外形氣動加熱的預測。“軸對稱比擬法”最早由Cooke［10］提出并由DeJarnette［11］等人通過實踐而成為一種實用的技術。運用軸對稱比擬技術可將三維邊界層問題轉化為軸對稱邊界層問題，所有適用于軸對稱條件的氣動加熱率的計算公式都可用于計算一般的三維外形。該氣動加熱快速預測方法與純粹的數值模擬技術相比具有耗時少、精度輕微損失的特點，適于預測比較復雜的三維外形的熱流。

因此，本文針對工程上幾何形狀比較復雜的三維高速飛行器外形，基于空間推進求解無粘方程的高效數值方法和流線追蹤工程方法，將二者結合建立了有一定精度的可滿足工程應用要求的氣動力熱特性快速預測方法。氣動力預測基于高效空間推進數值方法，氣動熱預測采用流線追蹤法基于空間推進無粘數值解與軸對稱比擬的氣動加熱積分關系式。

1 氣動力熱性能高效預測方法

1.1 無粘流場空間推進法

本文無粘流場計算通過空間推進法數值求解Euler方程完成。推進方向采用預測-校正兩步法，側面通量求解基于Godunov格式，根據Van Leer方法構造二階精度，對梯度用Minmod限制器進行單調性限制。方法細節見文獻［12-14］。

1.2 氣動加熱預測方法

所有氣動加熱快速預測方法的基礎是基于對邊界層方程的變換和求解。本文氣動加熱預測方法采用流線追蹤法，首先通過采用上述空間推進的方法數值求解Euler方程獲得無粘流場，然后利用無粘流場計算結果確定物面流線及尺度因子，沿流線運用Zoby［9］的方法積分得到物面氣動加熱率。在積分的每一步計算當地的邊界層厚度，插值得到無粘流場在該位置處的流動參數作為邊界層外緣參數。

2 網格劃分

本文在推進計算前，先生成整體三維網格，然后根據推進過程，下一個站位二維網格從三維網格中插值得到，這樣可以有效地解決推進過程中生成網格的問題。網格生成采用單體結構網格劃分技術，劃分步驟為，首先對每個軸向站位生成二維網格面，然后連接從頭部至底部各站位網格面最終建立三維空間網格。各站位的二維網格面通過數值求解拋物化方程從物面逐步推進到外流場邊界完成，由于拋物化推進方法較成熟，容易控制網格的稀密及正交等特性，因此對一般外形均可得到質量較好的網格分布。

3 高速飛行器外形氣動力熱算例研究

本文氣動力熱快速預測方法在高速飛行器外形氣動特性計算中得到成功應用，以下給出應用算例。

3.1 升力體外形

針對一種扁尖型升力體外形，運用本文所發展的方法完成了氣動力熱特性計算，外形如圖1所示。圖2給出了M=3.0和10.0，采用空間推進法(sm)計算的法向力系數CN，軸向力系數CA和對質心俯仰力矩系數Cmzg隨迎角變化曲線與時間相關法(tm)計算結果比較。兩種方法得到的氣動力系數符合很好，證明了空間推進法氣動特性預測的準確性。

圖1 升力體外形Fig.1 Lift-body configuration

圖2 空間推進法與時間相關法計算比較Fig.2 Aerodynamic coefficient compare between space and time marching methods

圖3為M=10.0、α=20°時采用空間推進法計算的升力體外形物面和流場典型截面壓力分布圖，圖4為物面熱流分布情況。

圖3 升力體外形物面和典型截面壓力分布Fig.3 Pressure distribution of the vehicle

圖4 升力體外形物面熱流分布(M=10.0、α=20°)Fig.4 Heat flux on the vehicle surface

3.2 空天飛機外形

根據文獻［15］給出的空天飛機外形(基本型)地面試驗結果，對本文所發展的氣動力熱快速預測方法進一步進行了驗證計算。該外形機身采用半圓和矩形組合橫截面，機翼為大后掠三角翼，翼端上彎。實驗模型全長0.29 m，翼前緣后掠角 68°，翼展0.1848 m，頭部半錐角10°，球頭半徑0.015 m，外形如圖5所示。

圖5 空天飛機試驗外形Fig.5 Space shuttle

圖6給出了M=8.0，α=20°時，使用空間推進法計算的機身中心迎風線(windward)和背分線(leeward)與實驗的比較，可以看出，計算(sm cal)和實驗(exp)規律一致，數據符合性也比較好;圖7和圖8分別為機身兩個橫截面和機翼弦向兩個截面計算與實驗比較結果，在背風面計算偏小，總體上兩者一致性較好。圖9為該狀態的計算網格，軸向站位網格點為30×58，網格總數約17萬。圖10為推進計算得到的五個站位壓力云圖，各站位壓力分布和激波形態是合理的。

圖6 M=8.0、α=20°空間推進法和實驗物面軸向壓力比較Fig.6 Pressure compare between experiment data and computational results of the symmetry plane

圖7 M=8.0、α=20°空間推進法和實驗物面周向壓力比較Fig.7 The circular pressure of two cross sections

圖8 M=8.0、α=20°空間推進法和實驗翼面弦向壓力比較Fig.8 The pressure of wing chord section

圖9 空間推進計算網格(M=8，α=20°)Fig.9 The corresponding grids

圖10 五個站位空間推進壓力云圖(M=8，α=20°)Fig.10 Pressure of five cross sections

以上計算工況在主頻為2.8GHz的CPU上計算時間見下表1，相對時間相關法，空間推進時間消耗下降接近一個量級。

表1 空間推進與時間迭代法計算效率對比Table 1 The efficiency of the two methods

圖11給出了M=4.0時空間推進法計算的CN，CA和Xcp隨迎角變化曲線與時間推進法及實驗結果比較。可以看出，空間推進法與時間相關法得到的氣動力系數均符合很好，計算與實驗法向力系數一致性很好，軸向力系數在大迎角時偏差增大，總的來看計算與實驗有較好的一致性。

圖12給出馬赫8、迎角20°，試驗狀態單位雷諾數為Re=3.54×107/m時，航天飛機機身迎風中心線熱流的計算值與實驗值的對比情況，計算時認為在機身前部就已完成轉捩并達到完全湍流流態。物面無粘流線分布見圖13。對比表明，雖然實驗數據本身散布度比較大，但計算的曲線經過散布點中間區域，且兩者趨勢相符合。實驗值與計算值的相對偏差范圍在2.7%～27.4%之間，21個測點中有17個測點的相對偏差小于20%。

圖11 M=4.0空間推進法計算與時間相關法和實驗數據比較Fig.11 Compare results of the two methods and tests

圖12 機身中心迎風線熱流計算與實驗比較(湍流)Fig.12 Central line heat flux on the windward side(turbulence)

圖13 空天飛機外形上、下表面無粘流線分布圖(M=8.0、α=20°)Fig.13 Stream lines on the surface(M=8.0、α=20°)

馬赫10、迎角18°時，試驗狀態的單位雷諾數為Re=2.2×106/m，計算時認為機身迎風線為完全層流流態，翼面迎風側已完成轉捩達到了完全湍流流態。圖14和圖15分別給出機身和翼面迎風線熱流的計算值與實驗值的對比。對于機身，實驗數據本身散布度不大，計算與實驗符合較好。實驗值與計算值的相對偏差范圍在0.1%～33.9%之間，13個測點中有12個測點的相對偏差小于20%。翼面實驗值與計算值的相對偏差范圍在0.3%～16.3%之間。

圖14 機身中心迎風線熱流計算與實驗比較(層流)Fig.14 Central line heat flux(laminar)

圖15 迎風翼面熱流計算與實驗對比曲線(湍流)Fig.15 Heat flux of the wing windward side(turbulence)

4 結論

本文建立的氣動力熱特性快速預測方法有較好的預測精度和效率。通過算例研究表明氣動力特性空間推進快速預測方法計算與實驗有很好的一致性，并且計算效率較高，同樣情況下空間推進法比時間相關法耗時可降低一個量級。氣動加熱快速預測方法與純粹的數值模擬技術相比具有耗時少、預測精度輕微損失的特點，計算結果與地面試驗數據的對比表明，此方法可以比較準確的預測迎風區的氣動加熱率，二者的相對偏差不高于20%。

本文所發展的氣動力熱特性快速預測技術，為高速飛行器氣動外形快速選型設計和氣動外形優化設計奠定了基礎。

［1］ Lawrence S L，et al.Application of an upwind algorithm to thethree-dimensional parabolized Navier-Stokes equations［R］.AIAA-87-1112，1987.

［2］ Gao Zhi.Significance and use of basic equation system governing high Reynolds number flows and diffusion-parabolized Navier-Stokes equations［J］.Advances in Mechanics，2005，35(3):427-438.(in Chinese))

高智.高雷諾數流動的控制方程體系和擴散拋物化Navier-Stokes方程組的意義和用途［J］.力學進展，2005，35(3):427-438.

［3］ Birch T J，Prince S A，Simpson G M.IMPNS:A space-marching solver for predicting theaerodynamic characteristics of high-speed missiles［R］.AIAA 2002-4512.

［4］ Chen Bing，Xu Xu，Cai Guobiao.Numerical simulation of hypersonic Inlet flows using a parabolized Navier-Stokes algorithm［J］.Journal of Astronautics，2006，27(6):1319-1323.(in Chinese)

陳兵，徐旭，蔡國飆，等.用空間推進算法模擬高超聲速進氣道流場［J］.宇航學報，2006，27(6):1319-1323.

［5］ He Xuzhao，Le Jialing.Solving PNS equations with space marching method and it’s validation［J］.Acta Aerodynamic Sinica，2007，25(2):189-193.(in Chinese)

賀旭照，樂嘉陵.空間推進方法求解拋物化Navier-Stokes方程及其驗證［J］.空氣動力學學報，2007，25(2):189-193.

［6］ Han Dong，Fang Lei.Study of aerodynamic heating predictions for hypersonic aircrafts by tracing the surface stream trace［J］.Journal of Aerospace Power，2009，24(1):65-69.(in Chinese)

韓冬，方磊.基于流線跟蹤法的氣動熱工程計算研究［J］.航空動力學報，2009，24(1):65-69.

［7］ Lv Lili，Zhang Weiwei，Ye Zhengyin.Predicting heating distributions for hypersonic reentry bodies［J］.Chinese Journal of Applied Mechanics，2006，23(2):259-262.(in Chinese)

呂麗麗，張偉偉，葉正寅.高超聲速再入體表面熱流計算［J］.應用力學學報，2006，23(2):259-262.

［8］ Kang Honglin，Yan Chao，Li Tinghe.Numerical study of aeroheating predictions for hypersonic reentry bodies［J］.Journal of Beijing University of Aeronautics and Astronautics，2006，32(12):1395-1398.(in Chinese)

康宏琳，閻超，李亭鶴.高超聲速再入鈍頭體表面熱流計算［J］.北京航空航天大學學報，2006，32(12):1395-1398.

［9］ Zoby E V.Approximate convective-heating equations for hypersonic flows［R］.AIAA 79-1078R.

［10］Cook J C.An axially symmetric analogue for general three-dimensional boundary layers［J］.R＆M No.3200，British A.R.C.1961.

［11］DeJarnette，Fred R.Aerodynamic heating on complex configurations［C］//Technical papers conference on Advanced Technology for Future Space Systems，May 1979:179-188(AIAA paper 79-0891).

［12］Ji Chuqun.Missile aerodynamics［M］.Beijing:China Astronautic Publishing House，1996.(in Chinese))

紀楚群主編.導彈空氣動力學［M］.北京:宇航出版社，1996.

［13］Ji Chuqun，Li Jun.Application of the Godunov method in the numerical simulation of complicated flow fields［J］.Acta Aerodynamica Sinica，2000，18(2):132-137.(in Chinese)

紀楚群，李駿.Godunov方法在復雜流場數值模擬中的應用［J］.空氣動力學學報，2000，18(2):132-137.

［14］Mendenhall M R.Tactical missile aerodynamics:prediction methodology［M］.Beijing:China Astronautic Publishing House，1999.(in Chinese) Mendenhall M R(美)主編，洪金森等譯.戰術導彈空氣動力學［M］.北京:宇航出版社，1999.

［15］Li Suxun.Hypersonic flow characteristics for typical shapes［M］.Beijing:National Defense Industry Press，2007.(in Chinese)

李素循.典型外形高超聲速流動特性［M］.北京:國防工業出版社，2007.

Fast prediction based on effective numerical method for aerodynamic force and heat of high speed vehicles

Wang Rong*，Zhang Xuejun，Ji Chuqun
(The First Institute，China Academy of Aerospace Aerodynamics，Beijing 100074，China)

A fast aerodynamic characteristics prediction technique is developed to predict aerodynamic force and heat，combining effective space marching numerical method with engineering method based on surface stream tracking technique.Hypersonic aerodynamic forces are obtained effectively by solving the Euler equations numerically with fast space marching method.In order to calculate the heat flux，an engineering method，called tracking the surface stream trace based on the axisymmetric analogue technique and boundary layer theory，is adopted and developed to be applicable for complex geometry.The inviscid flow parameters and surface stream lines required in the engineering method are obtained from above mentioned inviscous numerical flow fields.The proposed method for aerodynamic force and heat prediction is applied to high speed vehicles，the results show that the space marching method saves the computational cost enormously，approximately one order less comparing to the time marching method，while the two methods have similar precision in aerodynamic force evaluation.The predicted error of heat flux is within 20%comparing to that of wind tunnel test.The technique promoted is an effective approach suitable for the needs of fast aerodynamic configuration design and optimization.

space marching;stream trace;high speed vehicle

V211.3

A

10.7638/kqdlxxb-2014.0003

0258-1825(2015)04-0530-06

2014-02-11;

2014-05-09

王榮*(1981-)，男，高級工程師，研究方向:氣動計算與布局優化設計.E-mail:dilect@126.com

王榮，張學軍，紀楚群.基于高效數值方法的高速飛行器氣動力熱特性快速預測研究［J］.空氣動力學學報，2015，33(4):530-535.

10.7638/kqdlxxb-2014.0003 Wang R，Zhang X J，Ji C Q.Fast prediction based on effective numerical method for aerodynamic force and heat of high speed vehicles［J］.Acta Aerodynamica Sinica，2015，33(4):530-535.