基于POD方法的復雜外形飛行器熱環(huán)境快速預測方法

聶春生,黃建棟,王 迅,李 宇

(中國運載火箭技術(shù)研究院 空間物理重點實驗室,北京 100076)

基于POD方法的復雜外形飛行器熱環(huán)境快速預測方法

聶春生*,黃建棟,王 迅,李 宇

(中國運載火箭技術(shù)研究院 空間物理重點實驗室,北京 100076)

使用CFD數(shù)值方法獲得三維熱環(huán)境數(shù)據(jù)庫，利用本征正交分解對CFD數(shù)據(jù)庫進行降階處理，結(jié)合相應的基系數(shù)插值方法，快速預測出未知狀態(tài)下滿足精度要求的表面熱環(huán)境參數(shù)。以Hermes外形為研究對象，建立熱環(huán)境數(shù)據(jù)庫，使用本文提出的方法預測未知來流狀態(tài)下模型表面熱流，并和CFD結(jié)果進行對比。結(jié)果表明，該方法可大幅提高計算效率，不損失預測精度。分析了數(shù)據(jù)庫包絡(luò)性對重構(gòu)精度的影響。實現(xiàn)了沿給定彈道的三維熱環(huán)境快速預測。該方法能夠反映真實的熱流空間分布特征，快速獲得精確的激波干擾區(qū)熱流，有力地彌補了工程算法的不足。

高超聲速流動；氣動熱；本征正交分解；代理模型

0 引 言

目前高超聲速飛行器在國際和國內(nèi)處于蓬勃發(fā)展時期，具備翼身融合外形驗證機的研制被投入了大量人力物力，美國的HTV及X-37B就是典型代表[1]。這種氣動外形表面熱流分布差異大，氣流干擾造成高熱流區(qū)域較多[2]，這些區(qū)域往往是防熱設(shè)計的重點，在設(shè)計階段，需要能夠快速提供氣動熱數(shù)據(jù)，以確定相應的熱防護系統(tǒng)方案及總體結(jié)構(gòu)方案。

目前氣動熱環(huán)境計算常用的方法有兩類——工程計算方法[3-5]和CFD數(shù)值計算[6-7]。工程算法計算速度快，但在處理復雜外形飛行器，尤其是在處理復雜環(huán)境、復雜區(qū)域的氣動熱問題時無法滿足計算精度要求。目前CFD發(fā)展迅速，已經(jīng)可以對高超聲速復雜外形飛行器進行全機流場數(shù)值模擬，但是沿彈道開展大規(guī)模數(shù)值計算的代價相當昂貴，也無法滿足工程設(shè)計階段對氣動熱環(huán)境快速分析的要求。發(fā)展一種可以達到與CFD相當?shù)挠嬎憔鹊膹碗s外形高超聲速飛行器熱環(huán)境快速預測方法，是目前高超聲速飛行器防/隔熱設(shè)計亟待解決的問題。

國內(nèi)外目前采用代理模型的方式實現(xiàn)了對流場的快速重建。其思想是通過分析已有數(shù)值計算結(jié)果，通過某種重構(gòu)策略得到未知狀態(tài)的參數(shù)。該方法在優(yōu)化設(shè)計、實時仿真、氣動彈性設(shè)計等領(lǐng)域獲得了重要的發(fā)展。對于如表面壓力、熱流分布這類物理量，由于數(shù)據(jù)量大且內(nèi)部存在相關(guān)性，采用簡單的一維或多維插值無法獲得合理的結(jié)果，需要采用諸如響應面法[8]、Kriging[9]、POD[10]等代理模型處理。

本征正交分解方法(POD)[10]可以提供一組經(jīng)驗基，這些基包含了所描述系統(tǒng)的絕大部分信息，用這些基的線性組合去近似初始問題，從而將一個高維問題轉(zhuǎn)化成近似的低維問題。Duan等[11-12]基于該方法求解二維翼型的反設(shè)計問題，在對基準翼型添加擾動并進行CFD計算獲得采樣解基礎(chǔ)上，基于數(shù)據(jù)擬合POD系數(shù)的方法就可以進行翼型反設(shè)計，可得到與目標壓力分布對應的翼型形狀。Bui Thanh Tan[10]提出在已知來流狀態(tài)(如馬赫數(shù)、攻角)對應壓力系數(shù)分布的情況下，對POD系數(shù)進行插值，可快速獲得未知流場結(jié)點的近似壓力系數(shù)分布，極大地提高了計算效率。此類文獻多針對壓力系數(shù)、黏性力系數(shù)等進行預測，并且大多進行單個狀態(tài)預測，針對復雜外形高超聲速飛行器表面熱流及激波干擾區(qū)熱流重構(gòu)的研究較少。

本文建立了基于POD方法的復雜外形飛行器表面熱流計算代理模型。使用該方法對Hermes航天飛機進行氣動熱環(huán)境的快速預測，分析了大面積、前緣等飛行器關(guān)鍵部件及激波干擾區(qū)的熱環(huán)境POD重構(gòu)結(jié)果和CFD結(jié)果的匹配程度和耗時特性，分析了數(shù)據(jù)庫分布對計算精度的影響規(guī)律。針對給定彈道，分析了該外形典型截面特征點的熱流規(guī)律，并與工程方法結(jié)果進行比較，驗證了該方法在飛行器預先設(shè)計階段對熱環(huán)境預測的適應性。

1 氣動熱計算代理模型

1.1 本征正交分解方法

數(shù)據(jù)庫中的任何樣本都可以投影到彼此正交的基函數(shù){φi}上，并用式(3)表示：

1.2 插值POD方法

(4) 通過插值方法，獲得自變量空間內(nèi)φi的近似連續(xù)函數(shù)；

(5) 在自變量空間中所關(guān)心點q處的所有POD基系數(shù)bi，并依據(jù)1.1節(jié)式(3)求得q處的近似流場U(q)。

1.3 基系數(shù)插值方法

如1.2節(jié)所述，假設(shè)流場信息是關(guān)于自變量空間來流參數(shù)的連續(xù)函數(shù)，則插值POD方法不再依賴原流動方程，只需要通過合適的插值方法構(gòu)造出POD基系數(shù)與自變量空間來流參數(shù)的函數(shù)關(guān)系，就可以獲得自變量空間取任意來流參數(shù)值時的近似流場信息。

通常來說樣條插值的精度和魯棒性最好，但要求在來流參數(shù)的各個維度上構(gòu)造出網(wǎng)格型分布的采樣，使得采樣數(shù)據(jù)庫異常龐大，所需計算量與自變量的維數(shù)成幾何量級式增長。本文采用的徑向基函數(shù)(Radial Basic Function，RBF)插值，是一種較為適合解決該問題的插值方法，其精度僅次于三次樣條，對于多變量的插值問題，其采樣數(shù)目的增加明顯少于三次樣條，其采樣量級的增加通常只與自變量的維數(shù)成正相關(guān)[13]。

2 基于插值POD方法的表面熱流計算分析

2.1 幾何模型與計算狀態(tài)

對Hermes基本型航天飛機建立幾何模型。如圖1所示，模型總長290 mm，半展長92.4 mm，后掠角68°[15]。熱環(huán)境數(shù)值計算采用結(jié)構(gòu)化網(wǎng)格，總網(wǎng)格量控制在400萬左右。

基于POD方法表面熱流重構(gòu)應用中，一個關(guān)鍵的問題是如何采用經(jīng)過驗證后的數(shù)值計算方法建立最佳的CFD數(shù)據(jù)庫。因為POD分析的優(yōu)劣不僅取決于數(shù)據(jù)庫的精度，而且很大程度上取決于作為數(shù)據(jù)庫中數(shù)據(jù)集合的分布特性。

共選飛行高度H、飛行馬赫數(shù)Ma、攻角α三個變量構(gòu)建數(shù)據(jù)庫，馬赫數(shù)取值14～22，攻角取值0°～15°，高度選擇(45～70) km，最終的CFD計算狀態(tài)點為34個。 圖2給出了數(shù)據(jù)庫采樣點對應的高度、馬赫數(shù)和攻角分布情況。Object point1和Object point2是未知狀態(tài)，下文對這兩個狀態(tài)的熱流密度進行POD預測，并與CFD結(jié)果進行比較。Object point1在數(shù)據(jù)庫的包絡(luò)內(nèi)，Object point2在數(shù)據(jù)庫的包絡(luò)外。圖2中Trajectory為2.3節(jié)算例所用彈道在高度-馬赫數(shù)和高度-攻角平面的分布，可見彈道完全被數(shù)據(jù)庫所包絡(luò)。

2.2 重構(gòu)結(jié)果影響因素分析

在基于2.1節(jié)所述獲得CFD數(shù)據(jù)庫的基礎(chǔ)上，對該外形Ma=19、H=54 km、α=14°來流條件下(即圖2中Object point1)進行POD重構(gòu)，將結(jié)果和CFD結(jié)果進行對比。圖3～圖5分別是迎風面、背風面、前緣熱流密度等值面?？梢奝OD重構(gòu)能夠清晰地反映飛行器迎風/背風大面積、翼前緣氣動加熱特征，在飛行器頭部、大面積均可獲得和CFD一致的結(jié)果。對于存在激波干擾的翼前緣區(qū)域，熱流密度分布特征及熱流量值與CFD結(jié)果也基本相同。圖6給出了沿x=65 mm、220 mm和z=5 mm、75 mm四個截面熱流密度的對比圖。在迎風/背風大面積等不受激波干擾影響的區(qū)域兩者幾乎完全一致；z=75 mm截面在機翼前緣激波干擾區(qū)域熱流較高，POD重構(gòu)和精確CFD結(jié)果有一定偏差，但偏差不超過15%，而翼面處POD結(jié)果與CFD結(jié)果幾乎一致。

下面對Ma=22、H=73 km、α=15°來流狀態(tài)下(即圖2中Object point2)進行POD重構(gòu)，此時待求狀態(tài)點在數(shù)據(jù)庫的包絡(luò)之外。沿4個給定截面的熱流密度分布見圖7，與精確CFD計算結(jié)果的對比表明，POD重構(gòu)結(jié)果誤差較大，特別是頭部、翼前緣激波干擾區(qū)POD結(jié)果達到CFD結(jié)果的兩倍以上?？梢姰敶鬆顟B(tài)點在POD采樣解集的包絡(luò)之外時，POD重構(gòu)的結(jié)果是不可信的。

2.3 沿彈道熱流預測

采用POD重構(gòu)的方法進行飛行器沿彈道熱流預測，并與工程方法預測結(jié)果進行對比，典型機動彈道參數(shù)如圖8所示。

圖9給出了6個熱流觀測點，包括迎風面軸向x=220 mm截面的展向位置z=5 mm、20 mm、35 mm共3個特征點和翼前緣激波干擾區(qū)域(Point2、Point3)、非干擾區(qū)(Point1)3個觀測點。

POD重構(gòu)結(jié)果表明飛行器迎風面熱流呈現(xiàn)從前緣高熱流到中心線低熱流逐漸變化的過程，這一結(jié)論與飛行器真實的氣動加熱規(guī)律相符，但目前常用的工程方法[4-5,15]無法很好地預示出迎風面熱環(huán)境這一空間變化特性。圖10(a)給出了不同方法得到的飛行器迎風面x=220mm截面三個觀察點沿彈道的熱流密度變化對比曲線，其中工程方法采用了基于跟蹤流線的軸對稱比擬工程算法[16]，可以看出工程計算結(jié)果與POD結(jié)果變化規(guī)律基本一致，但是工程結(jié)果無法反映該截面展向的熱流差異，而POD方法能夠獲得典型截面上各點的熱環(huán)境沿彈道變化特征，分辨率更高，反映了迎風面各點更加真實的氣動加熱水平。

圖11給出兩個攻角狀態(tài)下翼前緣熱流分布。15°攻角狀態(tài)下前緣干擾不明顯，而5°攻角狀態(tài)下前緣激波干擾嚴重，出現(xiàn)局部高熱流區(qū)。

對于飛行器機翼前緣激波干擾引起的局部熱流增大，干擾特征與彈道參數(shù)直接相關(guān)，工程快速預測方法無法考慮激波干擾。圖10(b)給出了不同方法得到的飛行器機翼前緣三個觀察點沿彈道的熱流密度變化對比，其中工程方法采用了基于后掠圓柱的前緣駐點線熱流預測方法[16]。可以看出在Point1處，工程計算結(jié)果與POD結(jié)果基本一致，說明兩種方法都能很好的預測前緣無干擾熱環(huán)境；在Point2、Point3處，由于工程算法不能直接考慮激波干擾的影響，工程結(jié)果與POD結(jié)果存在較大差異，小攻角狀態(tài)下Point3熱流顯著增高，POD結(jié)果反映出了圖11所示的翼前緣在小攻角狀態(tài)下激波干擾嚴重這一氣動加熱特征。

在熱環(huán)境數(shù)據(jù)庫完備的前提下，使用Intel Xeon主頻3.6 GHz服務器，針對所述彈道采用快速預測方法時間為33 s，用所述工程算法預測彈道熱環(huán)境需要28 s，而開展大規(guī)模并行數(shù)值計算僅僅單狀態(tài)最少需要48 h。由此可見使用提出的POD重構(gòu)策略，可以在與工程算法耗時相當?shù)幕A(chǔ)上，獲得更為精細的飛行器表面的熱流密度分布，能夠獲得較高精度的激波干擾區(qū)熱流。而使用CFD計算預測沿給定彈道熱流變化耗時巨大，不適合飛行器設(shè)計階段需求。

3 結(jié) 論

通過結(jié)合本征正交分解和RBF插值，建立了一種適宜復雜外形表面熱流預測的代理模型，給出了合理的氣動熱采樣點數(shù)據(jù)庫構(gòu)建方法。以Hermes飛行器為例，使用本文提出的基于POD方法的熱環(huán)境快速預測方法，對未知狀態(tài)以及某典型機動彈道下表面熱流分布進行預測。得到以下結(jié)論：

(1) 在未知來流狀態(tài)參數(shù)被數(shù)據(jù)庫包絡(luò)的情形下，使用POD方法可以獲得與CFD方法精確相當?shù)慕Y(jié)果；前緣激波干擾區(qū)熱流POD預測和CFD計算誤差不超過15%；在未知來流狀態(tài)參數(shù)不被數(shù)據(jù)庫包絡(luò)時，POD結(jié)果不可信。

(2) POD重構(gòu)方法可以獲得飛行器各個部位熱流沿給定彈道的變化規(guī)律，彌補工程算法不能考慮激波干擾的不足，能夠更加精細的反映飛行器表面熱流的空間分布特征。

(3) 在建立了熱流采樣數(shù)據(jù)庫的基礎(chǔ)上，POD重構(gòu)方法耗時和工程算法相近，可以滿足高超聲速飛行器初步設(shè)計階段對熱流快速預測的要求。

[1]李建林.臨近空間高超聲速飛行器發(fā)展研究[M].北京:中國宇航出版社,2012.

[2]國義軍,曾磊,張昊元,等.HTV2第二次飛行試驗氣動熱環(huán)境及失效模式分析[J].空氣動力學學報,2017,35(4):496-503.

[3]Caitlin H.USAF assesses X-37B OTV-1 data following successful test[R].Defense Weekly,U.S,2010.

[4]DeJarnette F R.A method for calculation laminar and turbulent convective heat transfer over bodies at an angle of attack[R].NASA CR-101678,1969,8-19.

[5]Hamilton H H,Greenet F A,DeJarnette F R.Approximate method for calculating heating rates on three-dimensional vehicles[J].Journal of Spacecraft and Rockets,1994,31(3):345-354.

[6]Gnoffo P A.An upwind-biased,point-implicit relaxation algorithm for viscous,compressible perfect-gas flows[R].NASA TP-2953,1990.

[7]王發(fā)民,沈月陽.高超聲速升力體氣動力氣動熱數(shù)值計算[J].空氣動力學學報,2001,19(4):439-445.

[8]Simpson T W,Booker A J,Ghosh D,et al.Approximation methods in multidisciplinary analysis and optimazation:a panel discussion[J].Structural and Multidisplinary Optimization,2004,27(5):303-313.

[9]Han Z H,Stefan G.Hierarchical kriging model for variable-fidelity surrogate modeling[R].AIAA 2012-5019,2012.

[10]Tan B T,Damodaran M,Willcox K E.Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition[J].AIAA Journal,2004,42(8):1505-1516.

[11]Duan Y H,Cai J S,Li Y Z.Gappy proper orthogonal decomposition-based two-step optimization for airfoil design[J].AIAA Journal,2012 ,50(4):968-971.

[12]Legresley P,Alonso J.Airfoil design optimization using reduced order models based on proper orthogonal decomposition[J].Proceedings of AIAA Space,2000,40(10):1954-1960.

[13]王仁宏.數(shù)值逼近[M].北京:高等教育出版社,2005.

[14]李海燕.高超聲速高溫氣體流場的數(shù)值模擬[D].綿陽:中國空氣動力學研究與發(fā)展中心,2007

[15]楊愷,高效偉.高超聲速氣動熱環(huán)境工程算法[J].導彈與航天運載技術(shù),2010,(4):19-23.

[16]王國雄.彈頭技術(shù)(上)[M].北京:中國宇航出版社,2009.

Fastaeroheatingpredictionmethodforcomplexshapevehiclesbasedonproperorthogonaldecomposition

NIE Chunsheng*,HUANG Jiandong,WANG Xun,LI Yu

(ScienceandTechnologyonSpacePhysicsLaboratory,CALT,Beijing100076,China)

Based on three dimensional aeroheating database calculated by Computational Fluid Dynamics (CFD) and degenerated by Proper Orthogonal Decomposition (POD),the aeroheating parameters in unspecified flow conditions were rapidly predicted by a proper interpolation method.The heat flux was solved by the current method for the Hermes configuration.The result proves that the proposed method,compared with CFD methods,can greatly improve efficiency without losing accuracy.The aeroheating along a given trajectory was predicted.The result shows that the proposed method can be applied to provide a relatively realistic aeroheating distribution on the surface of a whole vehicle,including the interfering regions of shocks.The proposed method overcomes the disadvantages in current empirical methods.

hypersonic flow; aerodynamic heating; POD; surrogate model

0258-1825(2017)06-0760-06

V211.22

A

10.7638/kqdlxxb-2015.0157

2015-08-25;

2015-09-30

聶春生*(1987-)，男，工程師，主要研究方向為高超聲速氣動熱防護研究.E-mail:niechunsheng6@163.com

聶春生,黃建棟,王迅,等.基于POD方法的復雜外形飛行器熱環(huán)境快速預測方法[J].空氣動力學學報,2017,35(6):760-765.doi:10.7638/kqdlxxb-2015.0157 NIE C S,HUANG J D,WANG X,et al.Fast aeroheating prediction method for complex shape vehicles based on proper orthogonal decomposition[J].Acta Aerodynamica Sinica,2017,35(6):760-765.