升力体气动布局，由于其良好的高超声速气动特性、有效的内部空间和有益的防热特性，越来越受到气动工作者的重视，该文应用张涵信院士的NND格式思想，用数值模拟的方法对类似于美国X-33的升力体外形进行了气动力、气动热研究。特别是在粘性项的离散和表面热流率的计算中应用了散度和梯度的积分定义，避免了数值计算的奇异性，保证了通量守恒，提高了计算效率。结果表明气动热计算结果与美国同类飞行器结果接近，驻点热流率与经典理论结果一致

刘宏

姚文秀

沈月阳

王发民

雷麦芳

English

Institute Of Mechanics,Chinese Academy of Sciences

第 19 卷 第 4 期2001 年 12 月　　　　　　　　　　空 　气 　动 　力 　学 　学 　报　ACTA AERODYNAMICA SINICA　　　　　　　　　　Vol . 19 ,No. 4Dec. ,2001收稿日期 :2000208207 ;修订日期 :2001205208.作者简介 :王发民 (19492) ,研究员 ,博士生导师 1　　文章编号 :025821825 (2001) 0420439207高超声速升力体气动力气动热数值计算王发民 , 　沈月阳 , 　姚文秀 , 　刘 　宏 , 　雷麦芳(中国科学院力学研究所 , 　北京 　100080)摘要 :升力体气动布局 ,由于其良好的高超声速气动特性、有效的内部空间和有益的防热特性 ,越来越受到气动工作者的重视 ,本文应用张涵信院士的 NND 格式思想 ,用数值模拟的方法对类似于美国 X233 的升力体外形进行了气动力、气动热研究。特别是在粘性项的离散和表面热流率的计算中应用了散度和梯度的积分定义 ,避免了数值计算的奇异性 ,保证了通量守恒 ,提高了计算效率。结果表明气动热计算结果与美国同类飞行器结果接近 ,驻点热流率与经典理论结果一致。关键词 : 升力体 ; 有限体积 TVD 格式 ; 气动力 (气动热)中图分类号 : V211. 3 　　文献标识码 : A 　　0 　引　言升力体构型航天器由于其高热荷载、低热流率再入物理特性 ,在大迎角下和高超声速时良好的气动力特性以及高效的内部体积利用率 ,是航天器气动构形的首选方案。美国国家航天局 (NASA)和空军 (USAF)自 1975 年以来 ,已经进行了有关升力体构型飞行器的 225 次飞行试验。包括 NASA 研究的 M2 - F 和 HL - 10 系列 ,美国空军的 X - 24 系列以及 X - 33、X - 34 系列。日本于 1996 年用 J - 1 火箭做了升力体构型的高超声速实验 ( M = 14) 。俄罗斯近年则做了用两级火箭推进到 200km 高空的亚轨道升力体实验 ,其再入大气层速度为 7. 25km/ s。本文以数值模拟的方法研究了升力体飞行器流场的气动热、气动力特性 ,其物理模型选择了类似于美国“冒险星”升力体气动布局方案 ,包括双垂尾、鸭式平尾和尾部的四个控制面。该模型飞行器周围流场比较复杂 ,物理内涵丰富 ,包括了流动分离、激波边界层干扰和激波相互作用等主要物理特征。在计算方法上 ,以张涵信院士的 NND 格式为依据 ,应用了一种对称性强的隐式格式 ,计算结果表明 ,该算法有较高的分辨率 ,且稳定性好 ,节省机时。特别是为了克服在复杂区域中因网格等因素造成的数值奇异性 ,在粘性项的离散中应用了流场的散度和梯度 ,保证了格式的通量守恒 ,使计算方法有比较好的网格适应性。1 　有限体积求解方法可压缩流动的三维雷诺平均方程的积分形式方程为　 　55 t∫Ωρ·dΩ + ϕΓρq ·d S = 0∫Ω55 t (ρq) dΩ + ϕΓρq ( q ·d S) = ϕΓτn1d S∫Ω5 e5 t dΩ + ϕΓ eq ·d S = ϕΓτn1 qd S + ϕΓ C 55 n ( pρ) d S[τ] = - p + 23μdiv( q) +μ[ε](1)其中 : q = ui + vj + wk ,Γ为控制体Ω的边界曲面 ,τn 为微元单位面积 d S 上的应力 , [ε]为应变张量。取来流参数 p ∞、a ∞、L 、μ∞,将变量无纲化 ,则雷诺数定义为 Re ∞=ρ∞a2∞/μ∞。对于无粘通量应用张涵信院士的 NND 格式 , 粘性应力没有采用差分格式替代粘性项中的微分形式 ,而是采用积分形式的粘性应力求法 ,因为在积分方程中使用差分格式会产生较大的误差 ,破坏扩散项的守恒性。让我们回到积分方程式 (1) ,其中应变张量表达式为[ε] =2 5 u5 x 5 v5 x + 5 u5 y 5 w5 x + 5 u5 z5 u5 y + 5 v5 x 2 5 v5 y 5 w5 y + 5 v5 z5 w5 x + 5 w5 y 5 w5 y + 5 v5 z 2 5 w5 z (2)设 n 为积分面外法向单位矢量 , n = nxi + ny j + nzk 则[ε] ·n = ( n ·grad u) i + ( n ·gradv) j + ( n ·gradw) k + nx ·grad u + ny ·gradv + nz ·gradw (3)　　对于任一标量 < ,由梯度定义grad < = li mV  0 ∮<d SV (4 a)可以得六面体小体积单元上梯度的表达式grad < = 1V 66l = 1 <lS (4 b)式中 V 为六面体体积 ,d S 是微元面积积分单元 , S l 为各面上的面积矢量 , <l 为 <在 l 上的值 ,则grad <≈ 1V 66l = 1 <lS lxi + 66l =1 <lS ly j + 66l = 1 <lS lzk (4 c)用 u , v , w 代替 <可得[ε] ·n = μ5 q5 n +μ·grad qn (5)方程 (1)中 p + 23μdiv ( q) 的离散 ,要用到散度的定义 ,div ( q) = li mv →0∮( n ·q) d SV(6)044 空 　气 　动 　力 　学 　学 　报 　　　　　　　　　　　　第 19 卷同样 ,应用有限体积近似 ,可以得 q 散度的近似表达式div ( q) ≈ 1V 66l = 1 ( q ·S) = 1V 66l = 1 [ ulS lx + vlS ly + wlS lz ] (7)粘性应力单位时间对单元所做的功为( [ε] ·n) ·q = μ 5 q5 n + grad qn ·q (8)进入该六面体单元的热量为C 55 n ( p/ρ) d s = Cn ·grad ( p/ρ) (9)可以证明该离散方程是守恒的 ,依 Nicolaides R a. 等的分析式 (2) ,有下面的误差估计 ,如果U ∈H′(Ω) ,则有‖Uh - U ‖w ≤ Cmax( S l) ( l = 1 ,6) | U| 1 ,Ω (10)这里 H′表示 Sobolev 空间 ,且| U| 21 ,Ω = ∮| ΔU| 2d S (11)说明在网格划分上 S l 不能太大。在时间域的离散采用一阶隐式格式。2 　湍流逼近在热流率计算中 ,湍流的壁面近似是很重要的。由于求解的流场相当复杂 ,作为升力体外形的先导性研究 ,需要计算 M = 27 ,α= 40°的比较苛刻的极限状况 ,因而选取了对流场适用面比较广的 B2L 模式。在三维湍流计算中 ,为了较好地计算壁面热流率 ,必须满足这样一个条件 :即离开壁面的第一个网格点的 y + 必须小于 5 ,且在 y + < 30 的区域内必须有相当多的网格点来模拟出近壁区迅速变化的速度分布和压力分布。在现有计算机条件下 ,对于这一复杂外形 ,几乎是不可能计算的。在近壁区 ,依据边界层内的量纲关系 ,定义速度和长度的特征量u3= τw/ρw (12)y 3 = μw/ ρwτw(13)则速度和长度的无量纲形式为 : u + = u/ u 3 , y + = y/ y 3 ,应用壁面函数近似 ,即u+= y + 　　　　　　　y + < 5 (14)u+= 5ln (0. 543 y +) 　　　　5 ≤ y + < 30 (15)而 u + 和 y +都是壁面摩擦力τw 的函数。可以看出 ,在数值计算中较好地求得τw 是正确反映壁面对流场影响的一个重要因素 ,本文应用了如下迭代。设离开壁面的第一个网格点的速度为 u1 ,距离为 y1。依据混合长度理论 u1 应满足即 　　μ1μ3 = 1k ln E y1y 3 (16)u1τw/ρw=1k ln Ey1τwρwμw(17)144第 4 期 　　　　　　　　　　王发民等 :高超声速升力体气动力气动热数值计算上式是一个非线性代数方程 ,仅有一个未知量τw ,在计算中 ,当 y + ≤1. 15 ,取τw =τk =μw| ¢ | 1为初值 ,应用迭代法求解。在湍流计算中 ,湍流尺度是用涡量的计算确定的。在壁面附近 ,本文应用积分法求ω¢ = rotq = ¢ xi + ¢ y j + ¢ zk (18)以求单元中心点的 ¢ x 为例 ,应用下式近似¢ x = 1V∫Ω¢ xdΩ = 1V∫Ω 5 w5 y - 5 v5 z dΩ = 1V∮s ( wd xd z - vd xd y) (19)采用有限体积离散得 ¢ x ≈ 1v 66l = 1 ( w lS ly - vlS ly) (20)同样表面热流率的计算中也应用了积分的方法 ,在靠近壁面的单元体积内对能量方程积分 (5 e5 t = 0) ,即∮S1eq ·d S = ∮S1τn ·qd S +∮S1qnd S (21)qn ·ΔS | w = - ∮S1eq ·d S +∮S1τn ·qd S + ∮S1 ≠wqnd S (22)离散后为qn =1ΔS | w ( 66i = 1 ( Fxu + Fyv + Fzw) - 66i = 1 eq ·ΔS + 65i =1 γγ - 1 μPr ·Δs (23)　　由于几何形状比较复杂 ,本文应用了网格分区处理技术 ,各区域间的信息交换是通过交界处的通量加权完成的。设区域 A 中 ( i , j , k) 单元的 P 面与区域 B 中 1～N 单元对应 ,用 Sp , S i分别表示所对应的面积 ,则 P 面的通量为 : Qp = 1Sp 6Ni = 1 ( Qi·S i) 。为了估计该类飞行器的热载荷 ,对 M ∞ = 27 ,α= 0°～40°的设计极限情况进行了计算 ,当M ∞= 27 ,α= 0°时 ,由 Rosing 近似公式 (6)求得的钝锥体的驻点热流率为 14. 3423 J / cm2 . s ,本文计算的热流率为 14. 3662 J / cm2 . s ,国外同类飞行器表面热流率在 M ∞ = 14. 3 时的计算结果为 (3) 16. 6 J / cm2 . s , 可以看出本文的计算在选型阶段是可信的。3 　计算结果及分析作为算例 ,我们对高马赫、大迎角的极限情况进行了计算 ,计算条件 : M ∞ = 27 　T ∞ =180. 65K　p ∞= 0. 4125 (N/ m2) 　ρ∞= 7. 955 e - 6 (kg/ m3) 　μ∞= 1. 216 e25 (kg/ m. . s) 采用等温壁面条件 : Tw = 1000K, (参考重心位于距机头 19. 2m) 。表 1 是该升力体在计算条件下三个迎角的气动力和最大热流率的数值计算结果244 空 　气 　动 　力 　学 　学 　报 　　　　　　　　　　　　第 19 卷表 1 　三个迎角的气动力和最大热流率Table 1 　Aerodynamic and maxima thermal flow rate for three angles of attackα= 0° α= 30° α= 40°轴向力系数 CA =FA0. 5ρ∞V2∞S ↑0. 08492 0. 11906 0. 12866法向力系数 CN =FA0. 5ρ∞V2∞S ↑ - 0. 005746 0. 503908 0. 776778俯仰力矩系数 Cm =M Z0. 5ρ∞V2∞SL ↑(相对于参考重心)0. 001483966 0. 0092166 0. 00848最大物面热流率 qmax(W/ cm2) 14. 3662 18. 4772 17. 4056　　从表中可以看出随着迎角α的增加轴向力系数 CA 、法向力系数 CN 都上升 ,在零迎角时法向力系数很小 ,这是因为上下物形总体差别不大 ,在迎角很大时压心变化不大 ;最大物面热流率在不同迎角情况下都是在同一量级 ,在零迎角时 ,最大热流率 qmax的位置位于机头的驻点 ;在大迎角时 ,最大热流率 qmax的位置位于机翼前缘。图 1 (a) 、(b)分别为α= 0°、α= 40°时对称面热流率分布 ,图中横坐标表示距机头距离 x 与图 1 　α= 0°、α= 40°时对称面热流率分布Fig. 1 　Thermal flow rate distribution for symmetric side whenα= 0°图 2 　α= 40°典型物型的热流率分布Fig. 2 　Thermal flow rate distribution for　 symmetric side whenα= 40°机身总长 L 的比 ,纵坐标表示该点的热流率 q与该迎角的最大热流率 qmax的比。在零迎角时上下表面热流率分布差别不大 ,在头部驻点处热流率约为 q/ qmax = 0. 77 ,其余区域热流率都比较小 ;而在大迎角时除头部的整个背风面热流率也都很小 ,迎风表面热流率都很高 ,在α=40°时驻点后移到 x/ L = 0. 1 处 , 热流率约为q/ qmax = 0. 74。图 2 为迎角α= 40°时 ,包括水平翼和垂尾的典型截面的物型和热流率分布。对迎风表面 ,在水平翼和机身结合处附近的热流率为迎344第 4 期 　　　　　　　　　　王发民等 :高超声速升力体气动力气动热数值计算风面的最小值 ,在水平翼的翼尖附近的热流率达到该截面的最大热流率 ;背风面热流率都很小 ,从此图还可以看出 ,在垂尾表面有一个小的峰值。图 3 为α= 40°时表面压力分布 ,从图可以总体看出升力体构型在高马赫数、大迎角时压力分布特征 ,整个背风面基本为低压区 ,下机身表面几乎都为高压区 ,从而产生很大的升力 ,这是升力体飞行器升力的主要来源 ,翼面产生的升力所占的比例很小。图 4 为α= 40°时对称面的压力分布。可以看出 , 在大迎角情况下的波系结构 ,在大迎角时迎风面有很强的贴体激波 ,产生高压。图 3 　α= 40°时表面压力分布Fig. 3 　Pressure distribution for surface when α= 40°　　　　图 4 　α= 40°时对称面的压力分布Fig. 4 　Pressure distribution for symmetric side when α= 40°图 5 和图 6 为α= 0°时 ,对称面和典型横截面的温度分布。显示出在α= 0°时包括尾迹的流场结构 ,清楚看出头激波压缩所产生的高温区 ,以及在尾迹膨胀后再压缩所产生的尾迹激波后的高温区。图 5 　α= 0°时表面温度分布Fig. 5 　Temperature distribution symmetric side whenα= 40°　　　图 6 　α= 0°时典型横截面温度分布Fig. 6 　Temperature distribution for typical transverse section whenα=0°本文所得的数值结果表明 :文中所采用的数值方法具有良好的稳定性 ,对激波有很高的分辨率 ,能够满足模拟大迎角、复杂外形的高超声速流场的要求 ,计算结果也反映了升力体飞行器的高超声速流动特点。所提出的湍流近壁近似处理 ,通量守恒的粘性项离散和表面热流积分计算方法在一定程度上保证了气动热计算的精度。所得结果说明了升力体构形飞行器在零迎角和大迎角高马赫数下的气动力和气动热特点。可作为气动选型的理论依据。参 　考 　文 　献 :[1 ] 　ANDERSON , JOHN D JR1Hypersonic and high temperature gas dynamics[M] . McGraw2Hill Book Company ,1989.[2 ] 　NICOLAIDES R A , PORSCHING T A , HALL C A. Covolume methods in computational fluid dynamics[ M] ,Computational fluiddynamics review 1995.444 空 　气 　动 　力 　学 　学 　报 　　　　　　　　　　　　第 19 卷[3 ] 　DAVID OLYNICK Trajectory2based thermal protection system sizing for an X233 winged vehicle concept[J ]1 Journal of spaceraft androckets 1998 ,35 ,May2June1[4 ] 　姚文秀、赵桂萍、王发民、于欣 1 整机高超声速流场数值模拟[ C]1 第九届高超声速气动力 (热) 学术交流会议论文集 ,广西北海 ,19971[5 ] 　陈让福 1 三维高超声速无粘定常流的数值模拟[D]1[博士学位论文 ] ,北京大学 1[6 ] 　卞荫贵 1 钟家康 1 高温边界层传热[M]1 北京 :科学出版社 ,19861Aerodynamic and aerothermal numerical simulationof hypersonic lifting body configurationWANG Fa2min , 　SHEN Yue2yang , 　YAo Wen2xiu , 　LIU Hong , 　LEI Mai2fang( Institute of Mechanics , Chinese Academy of Sciences , Beijing 　100080 , China)Abstract : The researchers are paying more attention to lifting body aerodynamic configurationbecause of its favorable aerodynamic character , effective inner space and beneficial thermal protectioncharacter. This paper uses Zhang Hanxin’s NND scheme to study the aerodynamic and aerothermalcharacter of a lifting body like X233 of America. Specially we utitize the integral definition of gradient anddivergence when dispersing the viscosity term and calculating the surface thermal flow rate , so we mayavoid the simulation singularity , assure the flux conservation and increase the calculation efficiency. Itshows that the aerothermal result approaches the result of similar vehicle of America and the thermal flowrate at the stagnation point is consistent with the classical result .Key words : lifting body ; finite volume TVD scheme ; aerodynamic ; aerothermodynamic544第 4 期 　　　　　　　　　　王发民等 :高超声速升力体气动力气动热数值计算

高超声速升力体气动力气动热数值计算

http://dspace.imech.ac.cn/bitstream/311007/17387/1/%e9%ab%98%e8%b6%85%e5%a3%b0%e9%80%9f%e5%8d%87%e5%8a%9b%e4%bd%93%e6%b0%94%e5%8a%a8%e5%8a%9b%e6%b0%94%e5%8a%a8%e7%83%ad%e6%95%b0%e5%80%bc%e8%ae%a1%e7%ae%97.pdf

高超声速升力体气动力气动热数值计算

Abstract

Similar works

Full text

Available Versions

Institute Of Mechanics,Chinese Academy of Sciences