Read 高超声速再入钝头体表面热流计算

� 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数值计算与边界层 内工程算法相结合的,对高超声速再入钝头体的表面热流进行了计算, 并将计算结果与风 洞实验数据进行了对比,两者吻合较好,验证了两种工程算法在计算高超声速飞行器热环境方 面的正确性.将两种工程算法与数值求解 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 +格式计算, 它兼有 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 高亮