高超声速飞行器翼面气动加热、辐射换热与瞬态热传导的耦合分析

陈 鑫，刘 莉，李昱霖，杨 武

（北京理工大学 宇航学院，北京100081）

随着飞行马赫数的提高，气动加热效应对高超声速飞行器的影响愈来愈不可忽略。准确地预测高超声速飞行器的热环境是飞行器设计及飞行试验成败的关键难点之一。根据国外在NASA研究开发计划中的经验，利用分层求解［1］的思路可以满足工程的需求，即进行气动加热计算时，只考虑与飞行器表面辐射换热达到热平衡，暂不考虑表面向结构内部的传导。

国外关注气动热、热传导及结构的耦合问题的工作开展较早，1989年文献［2］完成了流场、传热与结构的一体化耦合计算；文献［3］从数值模拟角度对气动加热与瞬态热传导耦合性的影响进行了分析，给出了部分模拟计算结果，指出了耦合性在分析持续气动加热问题中的必要性；文献［4］给出了一种把流场、热、结构耦合起来进行一体化数值模拟的方法，以简单的二维圆柱绕流为例，验证了该方法的可行性。文献［5］采用流场－热－结构耦合的计算方法，耦合气动加热，辐射换热及结构热传导，数值模拟了旋转体的外流场和壁面结构温度场的非定常过程。

上述参考文献指出了考虑气动加热、辐射换热与热传导耦合的必要性，但均采用了数值模拟气动加热和结构热传导过程，计算效率较低，计算模型多为二维简单模型。利用气动加热工程算法能极大地减少求解流场的计算时间，比传统的数值计算效率更高。本文利用以经验公式为基础的工程算法计算气动加热热流，考虑辐射换热热流，利用有限单元法计算瞬态热传导过程，并通过实时更新边界条件，实现了三维翼面的气动加热、辐射换热与结构热传导的耦合求解。

1 翼面气动加热热流

运用牛顿公式，气动加热热流为

式中：αh为焓值热交换系数，hr为恢复焓值，hw为飞行器壁面焓值。将飞行器升力面分为前缘与非前缘两部分分别进行气动加热热流计算。

1.1 飞行器翼面驻点热流

采用Kemp－Riddell公式计算翼面前缘驻点热流：

式中：vc＝7 900m／s，ρsl＝1.225kg／m3为大气海平面密度，RN为驻点曲率半径，ρ∞、v∞分别为无穷远来流密度及速度，h0为驻点总焓，h300为300K时空气的焓值。

1.2 飞行器翼面非驻点热流

采用片条理论，翼面弦向气动加热热流［6］为

式中：ρ＊、μ＊分别为参考焓下气体密度及粘性系数，ρe、μe分别为边界层外缘气体密度及粘性系数，Rex为参考焓下的气体雷诺数，带上标“＊”的参数为参考焓下的气流参数。

参考焓法的基本思想是利用不可压缩边界层理论求得的公式来估算可压缩边界层中的摩擦和传热。Eckert参考焓公式［7］最常用，在层流及湍流中均与试验吻合得较好。即：

式中：h＊为参考焓值，he为边界层外缘焓值。

1.3 边界层外缘及翼面壁面气流参数的确定

利用布泽曼理论可以考虑不同飞行参数及翼型参数对翼面表面温度分布的影响［8］。利用布泽曼理论得到迎风面边界层外缘气流参数：

式中：带下标e的参数为边界层外缘参数，带下标c的参数为翼面前缘参数，Cp为翼面表面压力系数，Cp0为翼面前缘压力系数，φ为迎风面某点处翼型切线与来流方向的夹角，λ为比热比，pe、p∞分别为翼面前缘和来流压强，Ma∞为来流马赫数，ρe、ρ∞为翼面前缘和来流密度，Te、T∞分别为翼面前缘和来流温度。

2 翼面辐射热流

飞行器壁面向外辐射的热流量为

式中：波尔兹曼常数σ＝57.6×10－9W／（m2·K4），ε为蒙皮表面辐射系数，Tw为壁面温度。

3 翼面热传导

热传导微分方程为

热流边界条件为

式中：k为导热系数，ρ为物体的密度，c为比热容，W为内部热源强度。

4 翼面气动加热、辐射换热与瞬态热传导耦合分析

高超声速飞行器飞行过程中，翼面壁面气动加热热流qw、辐射热流qr和壁面温度Tw相互耦合，通过边界条件相互作用，如图1所示。气动加热热流qw求解需要翼面壁面温度Tw，Tw求解需要以热流qw为边界条件；辐射换热热流qr求解需要翼面温度Tw，Tw求解需要以热流qr为边界条件求解。数学上表现为各部分的求解对象是其它求解对象的定解条件，即三部分求解对象之间是互为定解条件的。

图1 翼面传热各部分输入输出关系

传统处理气动加热、辐射换热、热传导需要联立三部分方程耦合求解，这样处理求解复杂，计算繁琐，实际操作困难。本文尝试采用一种准定常求解方法，其计算思路如下：

①给定高超声速飞行器典型弹道下随时间变化的飞行条件、初始时刻的翼面温度；

②由初始时刻的飞行条件、翼面温度求解翼面气动加热与辐射换热热流，并以此为热流边界条件求解结构热传导方程，得到单位时间步长的新的翼面温度分布；

③若新的温度分布与原有翼面温度分布满足收敛条件，则更新到下一步时间步长，回到①继续进行下一时间步迭代，直到整个弹道计算完毕，否则更新翼面温度返回②直到计算收敛。

给定典型弹道，考虑传热耦合性，气动加热、辐射换热与热传导耦合算法计算流程图如图2所示。

图2 翼面传热耦合计算流程图

5 算例分析

5.1 气动加热工程算法验证

对某一高超声速二维翼面进行气动加热数值及工程算法计算，飞行高度H＝30km，马赫数Ma＝5.0，来流温度T＝226.51K，攻角为0，翼面弦长为0.5m。CFD－FASTRAN是目前广泛应用于高超声速流的专用CFD商业软件［9］。本文高超声速气动热计算采用CFD－FASTRAN求解器计算热流分布。数值计算边界条件为压力远场、压力出口和壁面边界。翼面边界表面考虑辐射，辐射系数为0.85，时间离散采用隐式Jacobi点迭代方法进行推进，空间离散采用一阶Roe’s FDS格式。图3为利用数值计算的CFD－FASTRAN计算结果网格图。

图3 CFD二维翼面网格图

图4 为热流数值计算与工程计算的对比图，由图4可知，气动加热工程算法与数值计算热流在驻点及弦向上均较为接近，在驻点热流求解上，工程算法与数值求解结果吻合较好，约为60.7kW／m2。在x＝0.25m处，数值求解考虑翼型时出现热流较小的跳变，在弦向0.25m以后热流均略微小于工程算法热流。

图4 热流数值计算与工程计算对比

5.2 气动加热与热传导耦合分析

某一高超声速三维翼面如图5所示，翼面根弦长为0.5m，根梢比为0.45，展长为0.287 5m翼面表面蒙皮采用钛合金，厚度2.5mm。

图5 高超声速飞行器三维翼面结构图

飞行器助推飞行时间t＝0～50s，马赫数Ma＝1.5～5.5，飞行助推段H＝15～30km。

5.2.1 气动加热、辐射换热、热传导耦合分析与非耦合分析对比

不考虑气动加热与热传导耦合的输入输出关系，如图6所示。以0.001s为时间步长，在典型弹道上进行100步长时间后，考虑气动加热、辐射换热与热传导耦合的方法得到翼面温度，如图7所示，图中x为翼面弦长坐标，y为翼面展长坐标。非耦合翼面温度分布如图8所示，耦合与非耦合前缘最大相对温度差约为40%，且耦合方法求得的翼面分布更均匀，更能反映翼面真实的温度分布。非耦合求解时气动加热热流计算不能利用实时变化的翼面温度，随着飞行时间增加，气动加热热流将明显增大，不能反映真实的翼面热流。随着飞行时间的增加，非耦合求解的气动加热热流在翼面的积累效应将越来越明显，非耦合翼面分布温度越来越高于耦合翼面温度。故预测翼面热环境时必须考虑气动加热与瞬态热传导的耦合效应。

图6 非耦合各部分输入输出关系

图7 耦合100时间步长翼面温度分布

图8 非耦合100时间步长翼面温度分布

5.2.2 气动加热、辐射换热、热传导耦合分析与分层求解方法对比

分层求解思想指求解翼面温度时仅考虑气动加热热流与辐射换热热流瞬时平衡，即qw＝qr，从而得到表面辐射热平衡温度Tw。由美国NASA研究经验［1］可以看出，分层求解方法可以满足一定的工程需求。

针对上述高超声速翼面，以0.001s为时间步长进行耦合算法迭代，为了与分层求解对比研究，耦合求解初始温度为零时刻分层求解翼面的稳态温度，整个典型弹道求解需时约为55h。

耦合与分层求解根部弦向温度对比如图9所示，耦合与分层求解根部前缘温度－时间曲线如图10所示。如图9、图10所示，考虑耦合时整个翼面弦向温度分布与分层求解温度分布相当，前缘温度分布和弦向温度分布相差不大，但考虑耦合求解的温度分布低于分层求解的温度分布，对比发现耦合求解翼面弦向温度比分层求解温度降低约5%。同时，随着飞行时间的增加，耦合求解的温度越来越低于分层求解温度。通过与分层求解结果对比，可以看出，本文的耦合求解方法能够反映翼面的温度分布，是一种可行的计算气动加热、辐射换热、热传导耦合分析的方法，同时比分层求解所得的温度分布更低。考虑气动加热、辐射换热与瞬态热传导耦合后，仅考虑气动加热热流与辐射换热流平衡的分层求解方法求得的翼面温度将高于耦合求解翼面的温度。

图9 耦合与分层求解根部弦向温度对比

图10 耦合与分层求解根部前缘温度－时间图

6 结束语

本文给出了一种求解气动加热、辐射换热与瞬态热传导耦合求解的方法，首先保证了气动加热的有效性，其次实现了高超声速飞行器在典型弹道飞行条件下三维机翼的气动加热、辐射换热、瞬态热传导的耦合求解，与非耦合的求解方法对比，指出了考虑气动加热、辐射换热与瞬态热传导耦合的必要性。与工程上的分层求解方法相比，该耦合求解方法能获得相近的翼面温度分布，且耦合求解的翼面温度更低。

［1］SPAIN C，SOISTMANN D，PARK E.An overview of selected NASP aeroelastic studies at the NASA Langley Research Center，AIAA－90－5218［R］.1990.

［2］DECHAUMPHAI P，THORNTON E A，WIETING A R.Flow－thermal－structural study of aerodynamically heated leading edges［J］.J Spacecraft，1989，26（4）：201－209.

［3］程克明，吕英伟.飞行器持续气动加热的耦合性分析［J］.南京航空航天大学学报，2000，32（3）：150－155.CHENG Ke－ming，LV Ying－wei.An analysis of coupling feature in continuing gasdynamic heating over flight vehicle［J］.Journal of Nanjing University of Aeronautics and Astronautics，2000，32（2）：150－155.（in Chinese）

［4］黄唐，毛国良，姜贵庆，等.二维流场、热、结构一体化数值模拟［J］.空气动力学学报，2000，18（1）：115－119.HUANG Tang，MAO Guo－liang，JIANG Gui－qing，et al.Two dimensional coupled flow－thermal－structural numerical simulation［J］.Acta Aerodynamica Sinica，2000，18（1）：115－119.（in Chinese）

［5］杨荣，王强.高超声速旋转体气动加热、辐射换热与结构热传导的耦合数值分析［J］.上海航天，2009（4）：25－29.YANG Rong，WANG Qiang.Coupled numerical study on aroheating，radiative heat transfer and structure heat conduction for hypersonic bodies of revolution［J］.Aerospace Shanghai，2009（4）：25－29.（in Chinese）

［6］张志成.高超声速气动热和热防护［M］.北京：国防工业出版社，2003.ZHANG Zhi－cheng.Hypersonic heating and thermal protection［M ］.Beijing：National Defense Industry Press，2003.（in Chinese）

［7］ECKERT R G.Engineering relations for friction and heat transfer to surfaces in high velocity flow［J］.AIAA Journal，1955，22（8）：585－587.

［8］陈鑫，刘莉，李昱霖，等.高超声速飞行器翼面气动加热的工程计算方法［J］.弹箭与制导学报，2013，33（3）：133－137.CHEN Xin，LIU Li，LI Yu－lin，et al.Engineering calculation of aerodynamic heating for airfoils of hypersonic vehicles［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2013，33（3）：133－137.（in Chinese）

［9］LEI Jian，ZHAO Xiang，ZHANG S J.Hypersonic non－equilibrium computations for ionizing air，AIAA 2009－1591［R］.2009.