目前高超声速飞行器在国际和国内处于蓬勃发展时期,具备翼身融合外形验证机的研制被投入了大量人力物力,美国的HTV及X-37B就是典型代表[1]。这种气动外形表面热流分布差异大,气流干扰造成高热流区域较多[2],这些区域往往是防热设计的重点,在设计阶段,需要能够快速提供气动热数据,以确定相应的热防护系统方案及总体结构方案。
目前气动热环境计算常用的方法有两类——工程计算方法[3-5]和CFD数值计算[6-7]。工程算法计算速度快,但在处理复杂外形飞行器,尤其是在处理复杂环境、复杂区域的气动热问题时无法满足计算精度要求。目前CFD发展迅速,已经可以对高超声速复杂外形飞行器进行全机流场数值模拟,但是沿弹道开展大规模数值计算的代价相当昂贵,也无法满足工程设计阶段对气动热环境快速分析的要求。发展一种可以达到与CFD相当的计算精度的复杂外形高超声速飞行器热环境快速预测方法,是目前高超声速飞行器防/隔热设计亟待解决的问题。
国内外目前采用代理模型的方式实现了对流场的快速重建。其思想是通过分析已有数值计算结果,通过某种重构策略得到未知状态的参数。该方法在优化设计、实时仿真、气动弹性设计等领域获得了重要的发展。对于如表面压力、热流分布这类物理量,由于数据量大且内部存在相关性,采用简单的一维或多维插值无法获得合理的结果,需要采用诸如响应面法[8]、Kriging[9]、POD[10]等代理模型处理。
本征正交分解方法(POD)[10]可以提供一组经验基,这些基包含了所描述系统的绝大部分信息,用这些基的线性组合去近似初始问题,从而将一个高维问题转化成近似的低维问题。Duan等[11-12]基于该方法求解二维翼型的反设计问题,在对基准翼型添加扰动并进行CFD计算获得采样解基础上,基于数据拟合POD系数的方法就可以进行翼型反设计,可得到与目标压力分布对应的翼型形状。Bui Thanh Tan[10]提出在已知来流状态(如马赫数、攻角)对应压力系数分布的情况下,对POD系数进行插值,可快速获得未知流场结点的近似压力系数分布,极大地提高了计算效率。此类文献多针对压力系数、黏性力系数等进行预测,并且大多进行单个状态预测,针对复杂外形高超声速飞行器表面热流及激波干扰区热流重构的研究较少。
本文建立了基于POD方法的复杂外形飞行器表面热流计算代理模型。使用该方法对Hermes航天飞机进行气动热环境的快速预测,分析了大面积、前缘等飞行器关键部件及激波干扰区的热环境POD重构结果和CFD结果的匹配程度和耗时特性,分析了数据库分布对计算精度的影响规律。针对给定弹道,分析了该外形典型截面特征点的热流规律,并与工程方法结果进行比较,验证了该方法在飞行器预先设计阶段对热环境预测的适应性。
1 气动热计算代理模型 1.1 本征正交分解方法流场变量U(如表面压力、热流密度)可由来流参数(q1, q2, …, qn)唯一确定。对于M个采样流场集合构成的数据库{Uj}j=1M,通过POD方法可以获得数据库中采样流场的本征正交基{φi}i=1L。
基的求解中假设φi为数据库样本点的线性组合,如式(1)。POD方法的实质是寻找一组最优的基函数φi来使
$ {\varphi _i} = \sum\limits_{j = 1}^M {a_j^i} {U^j} $ | (1) |
$ \frac{1}{M}\sum\limits_{k = 1}^M {\left\langle {{U^j}, {U^k}} \right\rangle } a_k^i = {\lambda _i}a_j^i, \left( {j = 1, \cdots, M} \right) $ | (2) |
数据库中的任何样本都可以投影到彼此正交的基函数{φi}上,并用式(3)表示:
$ {U^j} = \sum\limits_{i = 1}^L {b_i^j} {\varphi _i} $ | (3) |
通过1.1节所述本征正交分解方法可以获得数据库{Uj}j=1M的一组最优标准正交基{φi}i=1L及数据库中采样流场在正交基上的投影系数bi。在此基础上通过某种插值方法,估计出bi在自变量空间(q1, q2, …, qn)内的近似连续函数,并依据1.1节中的式(3),组合出自变量空间内任意一点的流场。具体计算流程:
(1) 通过合理的采样策略,获得N维来流参数的M种组合方式{δkj}j=1M, k=(1, …, N),并依靠CFD获得相应的M个采样流场{Uj}j=1M;
(2) 由采样流场{Uj}j=1M,通过POD方法可以获得L(L≤M)组线性无关的正交基{φi}i=1L,简称POD基;
(3) 将采样流场{Uj}j=1M向POD基{φi}i=1L投影,获得投影系数bij。自变量空间内,第i个POD基的投影系数在M个离散点处的函数值的计算公式为bij=〈φi, Uj〉, j=(1, …, M);
(4) 通过插值方法,获得自变量空间内φi的近似连续函数;
(5) 在自变量空间中所关心点q处的所有POD基系数bi,并依据1.1节式(3)求得q处的近似流场U(q)。
1.3 基系数插值方法如1.2节所述,假设流场信息是关于自变量空间来流参数的连续函数,则插值POD方法不再依赖原流动方程,只需要通过合适的插值方法构造出POD基系数与自变量空间来流参数的函数关系,就可以获得自变量空间取任意来流参数值时的近似流场信息。
通常来说样条插值的精度和鲁棒性最好,但要求在来流参数的各个维度上构造出网格型分布的采样,使得采样数据库异常庞大,所需计算量与自变量的维数成几何量级式增长。本文采用的径向基函数(Radial Basic Function,RBF)插值,是一种较为适合解决该问题的插值方法,其精度仅次于三次样条,对于多变量的插值问题,其采样数目的增加明显少于三次样条,其采样量级的增加通常只与自变量的维数成正相关[13]。
2 基于插值POD方法的表面热流计算分析 2.1 几何模型与计算状态对Hermes基本型航天飞机建立几何模型。如图 1所示,模型总长290mm,半展长92.4mm,后掠角68°[15]。热环境数值计算采用结构化网格,总网格量控制在400万左右。
![]() |
图 1 Hermes基本型航天飞机构型 Figure 1 Basic configuration of shuttle Hermes |
基于POD方法表面热流重构应用中,一个关键的问题是如何采用经过验证后的数值计算方法建立最佳的CFD数据库。因为POD分析的优劣不仅取决于数据库的精度,而且很大程度上取决于作为数据库中数据集合的分布特性。
共选飞行高度H、飞行马赫数Ma、攻角α三个变量构建数据库,马赫数取值14~22,攻角取值0°~15°,高度选择(45~70)km,最终的CFD计算状态点为34个。图 2给出了数据库采样点对应的高度、马赫数和攻角分布情况。Object point1和Object point2是未知状态,下文对这两个状态的热流密度进行POD预测,并与CFD结果进行比较。Object point1在数据库的包络内,Object point2在数据库的包络外。图 2中Trajectory为2.3节算例所用弹道在高度-马赫数和高度-攻角平面的分布,可见弹道完全被数据库所包络。
![]() |
图 2 热环境数据库包络图 Figure 2 Database envelope and distribution |
在基于2.1节所述获得CFD数据库的基础上,对该外形Ma=19、H=54km、α=14°来流条件下(即图 2中Object point1)进行POD重构,将结果和CFD结果进行对比。图 3~图 5分别是迎风面、背风面、前缘热流密度等值面。可见POD重构能够清晰地反映飞行器迎风/背风大面积、翼前缘气动加热特征,在飞行器头部、大面积均可获得和CFD一致的结果。对于存在激波干扰的翼前缘区域,热流密度分布特征及热流量值与CFD结果也基本相同。图 6给出了沿x=65mm、220mm和z=5mm、75mm四个截面热流密度的对比图。在迎风/背风大面积等不受激波干扰影响的区域两者几乎完全一致;z=75mm截面在机翼前缘激波干扰区域热流较高,POD重构和精确CFD结果有一定偏差,但偏差不超过15%,而翼面处POD结果与CFD结果几乎一致。
![]() |
图 3 迎风面热流密度对比 Figure 3 Comparison of heat flux of windward surface |
![]() |
图 4 背风面热流密度对比 Figure 4 Comparison of heat flux of leeward surface |
![]() |
图 5 前缘热流密度对比 Figure 5 Comparison of heat flux density of leading edge |
![]() |
图 6 Object point1典型截面热流密度对比 Figure 6 Comparison of the typical cross-section heat flux of Object point1 |
下面对Ma=22、H=73km、α=15°来流状态下(即图 2中Object point2)进行POD重构,此时待求状态点在数据库的包络之外。沿4个给定截面的热流密度分布见图 7,与精确CFD计算结果的对比表明,POD重构结果误差较大,特别是头部、翼前缘激波干扰区POD结果达到CFD结果的两倍以上。可见当待求状态点在POD采样解集的包络之外时,POD重构的结果是不可信的。
![]() |
图 7 Object point2典型截面热流密度对比 Figure 7 Comparison of the typical cross-section heat flux of Object point2 |
采用POD重构的方法进行飞行器沿弹道热流预测,并与工程方法预测结果进行对比,典型机动弹道参数如图 8所示。
![]() |
图 8 弹道参数 Figure 8 Trajectory parameters |
图 9给出了6个热流观测点,包括迎风面轴向x=220mm截面的展向位置z=5mm、20mm、35mm共3个特征点和翼前缘激波干扰区域(Point2、Point3)、非干扰区(Point1)3个观测点。
![]() |
图 9 热流观测点 Figure 9 Heat flux observation points |
POD重构结果表明飞行器迎风面热流呈现从前缘高热流到中心线低热流逐渐变化的过程,这一结论与飞行器真实的气动加热规律相符,但目前常用的工程方法[4-5, 15]无法很好地预示出迎风面热环境这一空间变化特性。图 10(a)给出了不同方法得到的飞行器迎风面x=220mm截面三个观察点沿弹道的热流密度变化对比曲线,其中工程方法采用了基于跟踪流线的轴对称比拟工程算法[16],可以看出工程计算结果与POD结果变化规律基本一致,但是工程结果无法反映该截面展向的热流差异,而POD方法能够获得典型截面上各点的热环境沿弹道变化特征,分辨率更高,反映了迎风面各点更加真实的气动加热水平。
![]() |
图 10 沿弹道特征点热流密度对比 Figure 10 Comparison of the heat flux along the trajectory |
图 11给出两个攻角状态下翼前缘热流分布。15°攻角状态下前缘干扰不明显,而5°攻角状态下前缘激波干扰严重,出现局部高热流区。
![]() |
图 11 不同攻角状态前缘干扰热流(CFD结果) Figure 11 Heat flux of the leading edge interference area under different angle of attack (CFD results) |
对于飞行器机翼前缘激波干扰引起的局部热流增大,干扰特征与弹道参数直接相关,工程快速预测方法无法考虑激波干扰。图 10(b)给出了不同方法得到的飞行器机翼前缘三个观察点沿弹道的热流密度变化对比,其中工程方法采用了基于后掠圆柱的前缘驻点线热流预测方法[16]。可以看出在Point1处,工程计算结果与POD结果基本一致,说明两种方法都能很好的预测前缘无干扰热环境;在Point2、Point3处,由于工程算法不能直接考虑激波干扰的影响,工程结果与POD结果存在较大差异,小攻角状态下Point3热流显著增高,POD结果反映出了图 11所示的翼前缘在小攻角状态下激波干扰严重这一气动加热特征。
在热环境数据库完备的前提下,使用Intel Xeon主频3.6GHz服务器,针对所述弹道采用快速预测方法时间为33s,用所述工程算法预测弹道热环境需要28s,而开展大规模并行数值计算仅仅单状态最少需要48h。由此可见使用提出的POD重构策略,可以在与工程算法耗时相当的基础上,获得更为精细的飞行器表面的热流密度分布,能够获得较高精度的激波干扰区热流。而使用CFD计算预测沿给定弹道热流变化耗时巨大,不适合飞行器设计阶段需求。
3 结论通过结合本征正交分解和RBF插值,建立了一种适宜复杂外形表面热流预测的代理模型,给出了合理的气动热采样点数据库构建方法。以Hermes飞行器为例,使用本文提出的基于POD方法的热环境快速预测方法,对未知状态以及某典型机动弹道下表面热流分布进行预测。得到以下结论:
(1) 在未知来流状态参数被数据库包络的情形下,使用POD方法可以获得与CFD方法精确相当的结果;前缘激波干扰区热流POD预测和CFD计算误差不超过15%;在未知来流状态参数不被数据库包络时,POD结果不可信。
(2) POD重构方法可以获得飞行器各个部位热流沿给定弹道的变化规律,弥补工程算法不能考虑激波干扰的不足,能够更加精细的反映飞行器表面热流的空间分布特征。
(3) 在建立了热流采样数据库的基础上,POD重构方法耗时和工程算法相近,可以满足高超声速飞行器初步设计阶段对热流快速预测的要求。
[1] |
李建林. 临近空间高超声速飞行器发展研究[M]. 北京: 中国宇航出版社, 2012.
( ![]() |
[2] |
国义军, 曾磊, 张昊元, 等. HTV2第二次飞行试验气动热环境及失效模式分析[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. DOI:10.2514/3.26446 ( ![]() |
[6] |
Gnoffo P A. An upwind-biased, point-implicit relaxation algorithm for viscous, compressible perfect-gas flows[R]. NASA TP-2953, 1990.
( ![]() |
[7] |
王发民, 沈月阳. 高超声速升力体气动力气动热数值计算[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. DOI:10.2514/1.2159 ( ![]() |
[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. DOI:10.2514/1.J050997 ( ![]() |
[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] |
王仁宏. 数值逼近[M]. 北京: 高等教育出版社, 2005.
( ![]() |
[14] |
李海燕. 高超声速高温气体流场的数值模拟[D]. 绵阳: 中国空气动力学研究与发展中心, 2007
( ![]() |
[15] |
杨恺, 高效伟. 高超声速气动热环境工程算法[J]. 导弹与航天运载技术, 2010(4): 19-23. ( ![]() |
[16] |
王国雄. 弹头技术(上)[M]. 北京: 中国宇航出版社, 2009.
( ![]() |