� 2006年 12月第 32卷 第 12期
北 京 航 空 航 天 大 学 学 报
Journal of Beijing U niversity of A eronau tics and Astronau tics
Decem ber� 2006
Vo.l 32� N o�12
� 收稿日期: 2005�12�12
� 作者简介: 康宏琳 ( 1979- ),女,江西新余人,博士生, kang lin@ ase. buaa. edu. cn.
高超声速再入钝头体表面热流计算
康宏琳 � � 阎 � 超 � � 李亭鹤
(北京航空航天大学 航空科学与工程学院, 北京 100083 ) � �
郭迪龙
(中国科学院 力学研究所, 北京 100080)
� � 摘 � � � 要: 基于跟踪流线的轴对称比拟法,采用纯工程算法、Euler数值计算与边界层
内工程算法相结合的
方法
快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载
,对高超声速再入钝头体的表面热流进行了计算, 并将计算结果与风
洞实验数据进行了对比,两者吻合较好,验证了两种工程算法在计算高超声速飞行器热环境方
面的正确性.将两种工程算法与数值求解 N�S方程进行对比,表明工程算法在迎风面的热流计
算方面有较高的精度,节约了计算时间,很好地满足高超声速飞行器概念研究和设计的需要.
关 � 键 � 词: 热流; 高超声速; 数值计算; 工程算法
中图分类号: V 211. 3
文献标识码: A � � � � 文 章 编 号: 1001�5965( 2006) 12�1395�04
Numerica l study of aeroheating predictions for hypersonic reentry bodies
K ang Honglin� Yan Chao� LiT inghe
( S chool of Aeronau tic Science and T echnology, Be ijing Un iversity ofA eronau tics andA stronaut ics, Beijing 100083, Ch ina)
Guo D ilong
( Inst itu te ofM echan ics, Ch inese A cadem y of S cien ces, B eijing 100080, Ch ina)
Abstract: In o rder to ca lculate the heating rates on hypersonic b lunt bodies qu ick ly and effective ly, two
eng ineering methods w ere deve loped based upon the ax isymmetric analogue wh ich required tracking the invis�
cid surface stream lines. Eng ineering pred ictions in the boundary layer w ere used in both methods. How ever,
method No. 1 used eng ineering predict ions out of the boundary layer, and in methodNo. 2 the boundary�layer
edge conditions w ere determ ined by the CFD method so lved by Eu ler equations of three d imensions. Resu lts o f
the tw o methods are in good agreem entw ith the ex ist ing experim ental data, wh ich va lidates the correctness o f
them ethods in calcu lating the hea ting rates on hyperson ic vehic les. Results o f these two eng ineering methods
w ere also compared w ith those of CFD method solved by N�S equations, wh ich indicates that the eng ineering
methods save much ca lculating t ime and also have prec ision on predicting the heating rates on w indw ard.
Key words: heat ing rate; hyperson ic; num erical simulat ion; eng ineering code
� � 对于导弹、飞船、航天飞机、空天飞机以及高
超声速飞机等再入飞行器而言, 气动热及其防护
是整个飞行器成败的关键难题之一,因此关于再
入气动热防护计算的研究是十分必要的. 气动热
的计算方法有两类:数值求解 N �S方程以及工程
计算方法.近年来国内外发展了通过数值求解 N�
S方程及其近似形式来预测高超声速飞行器热环
境的 CFD程序 [ 1 - 4] , 但影响热流计算值的因素颇
多,且数值计算耗费机时巨大. 而为了充分模拟热
环境, 往往要计算飞行器在不同马赫数、攻角下的
多个状态,直接采用 CFD方法求解所有状态下的
表面热流是不太实用的. 国外在表面热流工程算
法上做了系统、深入的研究,涌现出 AEROHEAT,
INCHES, M IN IVER, CBAERO等多套工程应用软
件 [ 5- 7] ,以纯工程算法 CBAERO软件为例, 计算
X�33这类复杂外形飞行器的热环境, 在 2. 6GH z
Pentium 4台式 PC机上计算一个状态耗时仅
152 s,且与风洞实验吻合较好 [ 5 ] . 将这类软件应
用到航天器的概念研究和初步设计中,能大大节
约研制时间.国内在这方面的研究却鲜有报道.
Administrator
高亮
Administrator
高亮
本文基于跟踪流线的轴对称比拟法, 用纯工
程算法、Eu ler数值计算与边界层内工程算法相结
合的方法,对高超声速飞行器的热环境进行了计
算.这两种方法根据普朗特的边界层理论,将流场
分为边界层和边界层外缘的无粘流两部分, 无粘
流场采用工程计算或数值求解 Euler方程获得,
所得结果为边界层工程计算提供外缘条件.
1� 边界层外缘参数的计算
再入飞行器表面热环境, 与表面边界层外缘
的气动参数 (如速度、压力 )密切相关, 纯工程算
法和 Eu ler数值计算与边界层内工程算法相结合
的方法之间的差别,就在于两者在计算边界层外
缘参数的方法不同.
1. 1� 纯工程算法
在纯工程算法中, 为了与边界层内工程计算
方法匹配,采用近似的工程算法对边界层外缘参
数进行计算.以高超声速再入体普遍采用的球钝
锥和钝双锥外形为例, 本文将表面压力的计算分
为 3个区域:在球头部分声速点前,采用修正牛顿
压力公式; 从声速点到球锥相切点前, 采用普朗
特 �迈耶膨胀压力计算公式;对于球钝锥和钝双锥
的锥面,采用修正牛顿压力公式.速度则采用 Lees
的边界层边缘速度公式进行计算.
1. 2� Euler数值计算
在 Euler数值计算与边界层内工程算法相结
合的方法中,采用 CFD方法数值求解边界层外缘
参数. 由于高超声速下的边界层很薄,因此在计算
边界层外缘参数时忽略了粘性项的影响, 直接求
解 Euler方程, 这种做法同时降低了对网格密度
的要求,大大加快了计算速度. 这种方法较纯工程
算法而言,计算时间增加, 但更适用于求解复杂外
形物面的边界层外缘参数.
控制方程为三维可压非定常 Euler方程, 空
间离散采用 AUSM +
格式
pdf格式笔记格式下载页码格式下载公文格式下载简报格式下载
计算, 它兼有 Roe格式
的间断高分辨率和 V an Leer格式的计算效率, 而
且克服了二者的缺点,无需熵修正.时间离散采用
LU�SGS方法, 这种方法在三维情况下无条件稳
定,且三维求解只需 L及 U两次扫描, 标量求逆,
计算效率高.
边界条件:物面满足无穿透条件,即物面的法
向速度分量为 0. 采用 R iemann不变量关系来处
理远场无反射边界条件.
2� 表面热流计算
对于轴对称体外形, 无攻角情况下表面热流
的计算方法比较成熟.对于有攻角情况下的表面
热流, 常用的有轴对称比拟法、等价锥法、实验数
据关联法 3种处理方法. 本文采用基于跟踪流线
的轴对称比拟法,其基本思想是:基于小横向流近
似,认为物体表面边界层内流体流动方向基本上
与绕流物体的三维无粘流表面流线方向一致. 这
样,采用 Mager�s变换, 在物面以无粘表面流线作
为正交坐标系的一个坐标轴, 三维边界层方程就
可以化简为轴对称比拟边界层方程,可沿用轴对
称零攻角物体边界层内热流密度近似计算公式来
计算有攻角轴对称物体沿流线的热流密度分布.
这样, 跟踪多条流线, 即可得到整个物面上的热流
密度分布.相比等价锥等方法, 基于跟踪流线的轴
对称比拟法,理论基础更为坚实,在国外气动热工
程软件中得到广泛应用 [ 5- 7] .
2. 1� 无粘表面流线的计算
由于热流工程计算必须沿流线进行, 因此必
须先计算出绕流物体壁面上的流线.
在纯工程算法中,为表示物面无粘流线形状,
采用沿流线的笛卡尔坐标系, 如图 1. 令物面沿流
线的微分弧长为 ds= h�d�, 沿物面且垂直流线的
另一微分弧长为 d�= h d ( h�及 h 分别是坐标 �
及 的尺度因子 ), 垂直物面的法向微分距离为
dn,子午角为 !, 把流线方向与物面子午线方向的
夹角定义为流线方向角 ∀r,轴向距离记作 x.
图 1� 流线坐标系
在物面上垂直流线方向的横向流动速度与流
线方向的主流速度相比, 认为足够小,可忽略. 这
时将无粘运动欧拉方程应用于物面上 ( n = 0), 在
定常流动条件下,可得到物面无粘流线方程 [ 8] .
在 Euler数值计算与工程算法相结合的方法
中,流线计算直接采用流场计算的贴体三维坐标
系 ( �, �, #), 从无粘 Eu ler方程解中获得流场在物
面网格点上的速度分量 u, v, w 以及热力学变量
(压力 p和焓 h ),物面其余地方的流场变量值可
通过网格点上进行插值获得. 积分求解方程组
( 1),可以得到物面上的流线.
d�/d∃ = h�( u�x + v�y + w �z )
d�/d∃ = h�( u�x + v�y + w �z ) ( 1)
其中, 积分变量 ∃通过公式 ( 2)和沿流线的位移
1396 北 京 航 空 航 天 大 学 学 报 � � � � � � � � � � � � � � 2006年 �
Administrator
高亮
Administrator
高亮
Administrator
高亮
Administrator
高亮
Administrator
高亮
Administrator
高亮
Administrator
高亮
Administrator
高亮
相联系:
ds /d∃ = Vh� ( 2)
其中 V是物面上当地的无粘速度值.
V = ( u
2
+ v
2
+ w
2
)
1 /2
( 3)
h�由公式 ( 4)确定:
h� = ( x2� + y2� + z2�) 1 /2 ( 4)
2. 2� 热流估算公式
无粘表面流线求得后, 将轴对称边界层方法
应用在物面无粘流线上,从而获得三维边界层解.
2. 2. 1� 驻点热流密度计算公式
高速飞行器的驻点加热特别重要, 通过对多
个经验公式的对比, 最后选择了 Kemp�R iddell修
正公式:
qw s =
110 311. 7
RN
%
%0
1
2 V
V c
3. 15
h s - hw
h s - h300K
( 5)
式中, %0 = 1. 225 kg /m 3; V c= 7 900m /s; RN为驻点
曲率半径, m; qw s为驻点热流密度, kW /m 2; h s为滞
止焓值; hw 为壁面焓值; h300K为温度 300K时空气
的焓值.
为了验证该驻点热流公式的精度, 本文对钝
双锥在不同攻角下的驻点热流的计算值和实验
值 [ 9]进行了比较, 如表 1.可见 K emp�R idde ll修正
公式能够准确地计算出钝头体的驻点热流.
表 1� 钝双锥驻点热流值计算精度
攻角 /
( ! )
实验值 /
(W∀ m - 2 )
计算值 /
(W∀ m - 2 )
相对误差 /
%
0 443 200 433649. 5 2. 15
10 488 000 460678. 8 5. 59
20 421 300 443123. 7 5. 18
2. 2. 2� 非驻点区热流密度计算公式
非驻点区热流密度的计算采用参考焓方法.
该方法假设高速边界层与低速边界层结构相同,
直接将不可压缩流的摩擦系数公式应用于高速可
压缩流中,不同之处是将其中的气体性能参数取
为参考焓下的数值.参考焓是马赫数、壁温及来流
温度的函数, 本文选用在工程中被广泛应用的
Eckert参考焓方程:
h
*
= 0. 22h r + 0. 28he + 0. 5hw
h r = he + 5 # 10- 4 u2eP r1 /n ( 6)
式中, h* 为参考焓; h r为恢复焓; he为边界层外缘
焓值; ue为边界层外缘速度; P r为普朗特数.
对于非驻点区绕流, 由于在重要的峰值加热
区,边界层为层流状态,因此本文主要研究层流边
界层的加热问题.根据参考焓方法, 选取 Lees层
流热流密度分布公式:
qw
qw s
=
%* &* ue r
2 %*s &*s due
dx s
∃x0%* &* ue r2 ds - 12 � ( 7)
式中, %*为参考密度; &* 为参考动力粘性系数.
3� 方法验证及分析
为了验证上述方法的正确性, 本文计算了 2
个有详尽实验数据的三维模型的高超声速绕流情
况 [ 9- 10] .为叙述方便, 将纯工程算法称为方法 1,
Euler数值计算与边界层内工程算法相结合的方
法称为方法 2. 对于外形 1(钝锥 ), 其外形及方法
2计算时所采用的对称面网格如图 2.模型长度为
447mm, 头部曲率半径为 27. 94mm, 半锥角为
15!; 来流条件: 马赫数 Ma = 10. 6, 压力 P =
132. 1 Pa, 温度 T = 47. 3K, 壁温 Tw = 294. 44K,
攻角 ∋= 20!.对于外形 2 (钝双锥 ), 其外形及方
法 2计算时所采用的对称面网格如图 3. 模型长
度为 122. 24mm,头部曲率半径为 3. 835mm,前锥
半锥角为 12. 84!, 后锥半锥角为 7!; 来流条件:
M a = 9. 86, P = 59. 92 Pa, T = 48. 88K, Tw =
300. 0K, ∋= 0!, 10!, 20!.
图 2� 钝锥对称面网格 图 3� 钝双锥对称面网格
图 4� 钝锥 20!攻角下表面热流
图 4为钝锥在攻角为 20!时, 采用上述 2种方
法计算的沿迎风子午线 ( != 0!), 侧面子午线
(!= 90!)及背风子午线 (!= 180!)的热流分布与
相应实验结果 [ 9]的比较, 图中 !为子午角, 横坐
标表示轴向距离与头部曲率半径的比值, 纵坐标
为归一化热流,即当地热流与驻点热流的比值.计
算值与实验值吻合较好.可见, 采用基于跟踪流线
的轴对称比拟法, 能够较为准确地模拟高超声速
1397� 第 12期 � � � � � � � � � � � 康宏琳等:高超声速再入钝头体表面热流计算
Administrator
高亮
Administrator
高亮
Administrator
高亮
Administrator
高亮
Administrator
高亮
Administrator
矩形
Administrator
高亮
Administrator
高亮
Administrator
矩形
Administrator
高亮
钝头体在有攻角情况下的表面热流.
图 5~ 图 7为钝双锥在 0!, 10!, 20!攻角下,
采用方法 1、方法 2计算的迎风子午线表面热流
值,并与通过数值求解 N �S方程的计算结果 [ 3]以
及实验结果 [ 8]进行了对比. 可见, 在迎风面的热
流计算方面,本文的 2种工程算法都有较高的精
度,刻画了双锥交接处气流膨胀引起的热流变化,
可以为高超声速再入钝头体的设计提供准确的热
环境. 在计算速度方面, 在 1. 8 GH z AMD A thlon
( tm ) 64 Processo r的 PC机上,用方法 1计算耗时
约 45 s,方法 2计算耗时约 1 h, 而采用数值求解
N�S方程计算耗时达 4 h. 可见,本文的 2种方法
大大节约了计算时间, 适应了高超声速飞行器概
念研究和初步设计的需要.
图 5� 钝双锥 0!攻角下迎风子午线热流
图 6� 钝双锥 10!攻角下迎风子午线热流
图 7� 钝双锥 20!攻角下迎风子午线热流
4� 结 束 语
由于采用数值求解 N�S方程来预测高超声速
飞行器热环境耗费机时巨大, 难以应用到航天器
的概念研究和初步设计中, 因此有必要发展快速
有效的工程算法.本文基于跟踪流线的轴对称比
拟法, 采用纯工程算法、Eu ler数值计算与边界层
内工程算法相结合 2种方法, 对高超声速再入钝
头体在有攻角情况下的表面热流进行了计算, 与
数值求解 N�S方程相比, 2种方法大大节约了计
算时间;将计算结果与实验数据进行了比较,结果
表明, 2种方法均可较为准确地模拟迎风面的热
流分布.另外, 2种方法也各有利弊, 方法 1更迅
速快捷,方法 2则更易于推广至复杂外形飞行器
的计算.必须指出的是,本文 2种方法均采用小横
向流近似,但随着攻角的增大, 背风面流动出现分
离,横向流动增大, 小横向流近似不再成立,故背
风面热流误差较大.但由于在较大攻角情况下,背
风面热流远小于迎风面热流, 故在模拟飞行器热
环境时,关键还是在于求解迎风面热流.因此本文
2种方法能够很好地满足工程设计的需要.
参考文献 (References)
[ 1] Peter A G. C om pu tat ion al flu id dynam ics technology for hyp er�
son ic app lications[ R] . A IAA�2003�2829, 2003
[ 2] D ietik er J F, H offm ann K A. Challenges in th e h eat trans fer
com putations ofh igh sp eed f low s [ R] . A IAA�2004�0993, 2004
[ 3 ] 李君哲, 阎超. 气动热 CFD计算的格式效应研究 [ J] . 北京
航空航天大学学报, 2003, 29 ( 11) : 1022- 1025
L i Jun zhe, Yan Chao. Research on sch em e effect of com puta�
tional flu id dynam ics in aerotherm odynam ics[ J]. Journal ofB ei�
jing U n ivers ity ofAeronaut ics andA stronau t ics, 2003, 29( 11 ):
1022- 1025 ( in C h inese)
[ 4] 王发民, 沈月阳, 姚文秀. 高超声速升力体气动力气动热数
值计算 [ J] . 空气动力学学报, 2001, 19( 4 ): 439- 445
Wang Fam in, Shen Yu eyang, Yao W enx iu. A erodynam ic and
aerotherm al num erical s im ulation of hyperson ic lifting body con�
figurat ion [ J]. ACTA Aerodyn am ica S in ica, 2001, 19 ( 4 ) : 439
- 445( in Ch in ese)
[ 5] K inney D J. Aero�th erm odynam ics for con cep tual design [ R ].
AIAA�2004�0041, 2004
[ 6] Wu rster K E, Zoby E V, Th omp son R A. F low field and veh icle
param eter inf luence on resu lts of engineering aerotherm alm eth�
ods[ J] . Jou rnal of Spacecraft, 1991, 28 ( 1) : 16- 22
[ 7] H am ilton H H, Fran ceisA G. App rox imate m ethod for calculat�
ing heat ing rates on three�d im ens ional veh icles [ J ] . Journ al of
S pacecraft and Rockets, 1994, 31 ( 3) : 345- 354
[ 8] Dejarnette F R, Tai T C. A m ethod for calcu lating lam inar and
turbu len t convective heat transfer over bod ies at an angle of at�
tack [ R] . NASA�CR�101678, 1969
[ 9] M il ler C G. Experim ental and pred icted heating d istribut ion s for
bicon ics at in cidence in air atm ach 10 [ R] . NASA�TP�2334,
1984
[ 10] C leary JW. E ffects of angle of attack and blun tn ess on lam inar
h eating rate d istribut ion of a 15!con e at a mach number of10. 6
[R ] . NASA�TND�5450, 1969
(责任编辑:文丽芳 )
1398 北 京 航 空 航 天 大 学 学 报 � � � � � � � � � � � � � � 2006年 �
Administrator
高亮