基于有限体积法的近空间高超声速飞行器辐射特性仿真研究

基于有限体积法的近空间高超声速飞行器辐射特性仿真研究

1热辐射源红外特性的研究近距离飞机是一种新型的战略遏制和战术应用武器平台。这是世界上所有主要军事能力的集合和研究的热点。随着红外探测、跟踪和制导技术的发展,空中军事目标面临日益严峻的红外威胁,具备红外隐身能力已成为现代先进飞行器提高自身生存和突防能力的必然要求。开展近空间高超声速飞行器红外辐射特性研究对于其红外隐身设计具有十分重要的意义。获取目标的红外特性数据可通过实验测量和数值仿真来实现。由于实验测量耗资巨大,实验结果容易受诸多环境因素的影响,建立红外特性理论模型,通过数值仿真来获取目标的红外辐射特征被广泛采用。国外这方面的研究开展得较早,并已开发出相关计算软件,如北约部分国家联合开发的空中目标的红外辐射模型NIRATAMNATO;近年来,国内开始重视飞行器隐身及相关方面的研究工作,罗明东等研究了非加力涡轮发动机排气系统红外辐射强度的数值计算方法;张小英等研究了轴对称矢量喷管在光谱2100～5260cm-1范围内红外特性的数值计算方法;郝金波等在给定浓度场、温度场、粒子和气体组分辐射物性参数的条件下,将贴体坐标系下有限体积法用于主动段尾喷焰红外光谱特性计算。但研究重点主要集中在各类常规飞行器的发动机热部件及尾喷焰。对于近空间高超声速飞行器来说,气动加热引起的热辐射是其主要的红外辐射源。本文从红外隐身设计的角度,针对飞行马赫数Ma=4～7、飞行高度H=15km～30km、具有长时间热载荷和非烧蚀构型等特点的近空间高超声速飞行器气动热红外特性进行了理论建模和数值仿真研究。希望以此为基础能够进一步建立一套完整的数值预测体系,可以为近空间高超声速飞行器红外隐身设计提供参考依据。2辐射场-流场当飞行器在近空间大气层内以高超声速飞行时,强烈的气动热会引起周围介质(空气)各点的温度、压力等热力学状态在空间上的剧烈变化,从而使介质表现为强烈的参与性和不均匀性。但由于介质热力学状态只在飞行器近域空间变化比较剧烈,在无穷远处则趋向于外界环境。为此,在远场处引入一“假想面”将飞行器周围空间划分为两部分,如图1所示。“假想面”外围介质认为是常物性的或透明的,只要知道“假想界”处辐射参量,就可得到任意远处的辐射参量值,这样问题便转化为求解封闭系统内参与性介质的辐射换热问题。通常,辐射场和流场是耦合在一起的,但这里所研究的对象,速度不大于3km/s,流场最高温度一般不超过3000K。根据目前文献,在此情况下可以不考虑它们之间的干扰效应,认为辐射场与流场是可以解耦的。采用解耦计算的具体步骤为:首先基于气体动力学方程计算流场参量;然后以已知流场参量为基础,计算介质辐射物性;最后通过求解介质内辐射传输方程获得相应辐射参量。3流场数值法3.1navier-st对流通量方程的数学模型假设空气是由多种组分气体均匀混合而成的完全气体模型。为了能够获得较为精确的流场参量,考虑可压缩粘性流动的Navier-Stokes方程组。在直角坐标系下,其守恒形式为:∂U∂t+∂(FU+FV)∂x+∂(GU+GV)∂y+∂(ΗU+ΗV)∂z=0(1)∂U∂t+∂(FU+FV)∂x+∂(GU+GV)∂y+∂(HU+HV)∂z=0(1)式中,U为守恒变量;FU、GU、HU为对流通量;FV、GV、HV为粘性通量。3.2类tvd模型本文采用Jameson中心差分格式求解上述控制方程。人工粘性项采用二阶与四阶差分的组合形式;针对高超声速流动问题,引入类TVD变量修正;为了加速求解过程,采用当地时间步长和隐式残值光顺等加速收敛措施。湍流模型采用Baldwin-Lomax代数模型。4辐射场数值法4.1热流密度qkr的基本原理气体分子遵循Rayleigh散射,一般工程中,其散射能量很小,故忽略不计。基于波段Δλ内的发射、吸收性非灰介质内辐射传输方程表达式为:dΙk(→r,ˆs)ds=-κk(→r)Ιk(→r,ˆs)+κk(→r)Ιbk(→r)(2)dIk(r⃗,sˆ)ds=−κk(r⃗)Ik(r⃗,sˆ)+κk(r⃗)Ibk(r⃗)(2)式中,下标k表示波段Δλ区域(如3～5μm);ˆssˆ为单位方向矢量;→rr⃗为空间位置矢量;Ik(→rr⃗,ˆs)sˆ)为位置→rr⃗处、ˆssˆ方向上Δλ内的辐射强度;Ibk(→rr⃗)为位置→rr⃗处黑体辐射强度;κk(→r)为位置→r处介质Δλ内吸收系数平均当量。在已知位置→r处各方向辐射强度值后,位置→r处的辐射热流密度→qK(→r)的表达式为:→qk(→r)=∫Ω=4πΙk(→r,ˆs)ˆsdΩ(3)其中,Ω为空间立体角。认为飞行器壁面为不透明、漫发射、漫反射面,其边界条件为:Ιk,w(→rw,ˆs)=εk,wΙbk,w(→rw)+1-εk,wπ∫nw⋅si＜0Ιk‚w(→rw‚ˆsi)|ˆnw⋅ˆsi|dΩi(4)式中,εk,w为壁面Δλ波段内发射率平均当量;ˆnw为壁面单位法向矢量。4.2辅助网格单元在求解RTE的方法中,有限体积法是一种非常强有力的方法,它可确保辐射能量整体守恒,对不规则边界适应性强,可以使用与流场相同的计算网格,方便地与流场求解的程序相结合。因此本文采用FVM求解辐射传输方程。采用FVM求解辐射传输方程时,需要进行空间离散和角度离散。空间离散与流场计算相同;角度离散就是将整个4π空间划分为互不重叠的若干微元立体角,本文采用文献提出的角度离散方法,如图2所示。在控制体Vp和微元立体角Ωm(如图3所示)内对辐射传输方程积分得到其在→sm方向的离散形式为:∑j=e,w,s,n,t,b[Ιmk,c,jSjDmj]≈[-κk,ΡΙmk,Ρ+κk,ΡΙbk,Ρ]ΩmVΡ(5)式中,Ωm为第m个微元立体角;Imk,c,j为控制体j面上Ωm方向的k区域内辐射强度;Imk,Ρ为控制体中心P处Ωm方向的k区域内辐射强度;Ibk,P为中心P处k区域内黑体辐射强度;Sj为控制体j面的面积;Dmj为控制体j面的方向权值,其表达式为:Dmj=∫Ωm(ˆsm⋅ˆnj)dΩm(6)式中,→sm表示微元立体角Ωm中心的单位方向矢量;→nj表示表面Sj的单位外法矢量。为了将控制体表面上的辐射强度和体内节点处的辐射强度联系起来,本文选用阶梯格式(thestepscheme),即:Imk,c,j=Imk,p(7)则方程(5)可写成如下形式:amΡΙmk,Ρ+∑J=E,W,S,Ν,Τ,BamJΙmk,J=bmk,Ρ(8)式中各项系数分别为:amΡ=∑j=e,w,s,n,t,bmax[SjDmj,0]+κk,ΡVΡΩmamJ=min[SjDmj,0]bmk,p=κk,pBk,ΤpσΤ4pπVpΩm式中,下标J=E,W,S,N,T,B分别表示与控制体VP相邻的六个控制体中心节点;Bk,TP表示温度TP下Δλ波段内的黑体辐射能量占其总辐射能量的份额。壁面边界条件(4)的离散形式为:Ιmk,w=εk,wBk,ΤwσΤ4wπ+1-εk,wπ∑Dlw＜0Ιlk,w|Dlw|(9)对于“假想面”边界的处理,不考虑环境辐射源影响,则就没有外部辐射能量进入封闭系统内部,假设系统通过“假想面”辐射到外界去的辐射能量,将会被为外界大气完全吸收掉,不再反射回系统内部,也即系统内辐射能“只出不进”。为了实现这一点,引入一层辅助网格单元作为最外层边界,如图4所示。设DmΡ为系统内部P处控制体“单元交界面”上的方向权值,Imk,Ρ为P处的辐射强度,Imk,Μ为M处的辐射强度,则当DmΡ>0时,Imk,Μ=Imk,Ρ;当DmΡ≤0时,Imk,Μ保持原值不变。对于得到的离散方程组,本文采用CGSTAB算法求解。4.3气体辐射吸收系数计算在求解辐射传输方程之前需要先知道气体的辐射吸收系数。把每个控制体内的气体看作是等物性的,在已知其热力学参数的基础上,利用逐线计算法计算气体辐射吸收系数,逐个计算所有控制体内气体的吸收系数即得封闭系统内气体吸收系数的空间分布。在计算中考虑光谱线的压力增宽和多普勒增宽混合效应,所需谱线参数由HITRAN数据库获得。对于波段Δλ区域的吸收系数平均当量,本文采用普朗克平均法计算,即:κk=∫λ2λ1κλΙbλdλ∫λ2λ1Ιbλdλ(10)此外,假设组成空气的各组分气体之间互不影响,空气介质的吸收系数等于各组分气体吸收系数之和。5计算结果和分析5.1计算程序的初步验证1壁温ma模型为钝锥体:球头半径为27.94mm,半锥角为15°,底部圆半径为152.4mm。计算条件:Ma∞=10.6,P∞=131.97pa,T∞=47.34K,给定壁温为Tw=294.44K。图5为零攻角时钝锥体表面热流密度分布与文献提供的试验数据的比较,计算值与试验值符合得很好(图中q0为驻点热流密度值)。2辐射热流密度q为了在一定程度上验证RTE计算方法和程序的正确性和可靠性,设计了一个简单球形黑体模型:球体半径为1m,表面为黑体,温度为500K;周围为透明介质。设R为到球体中心的距离,图6表示辐射热流密度q随R的变化情况。从图中可以看出,计算结果比理论值稍低些,最大相对误差约为30%,这主要是由空间离散和角度离散引起的,但计算结果与理论值的变化趋势是一致的,完全能反映出球形黑体周围空间热辐射的分布规律,这表明方法和程序是正确可靠的。5.2空气流场模拟结果乘波构型被认为是高超声速飞行器的最佳构型。以乘波构型为实例,利用上述建立的方法和计算程序数值仿真了其某一飞行状态下的3～5μm波段内红外辐射特性的空间分布。基于高超声速细长圆锥绕流的近似解,根据文献中的设计方法构建了简单的锥导乘波构型。基本参数为:总长度为L=5m,设计马赫数为Ma=6,压缩角为θc=β-ωK=4°,外形如图7所示。利用Gridgen软件进行网格划分。空间网格的数量为121×51×61,在物面和前缘附近进行了加密,计算网格如图8所示。流场计算条件:攻角为0°,来流马赫数Ma=6,飞行高度H=20km。数值边界条件:给定壁面温度Tw=300K;在入口处取自由来流条件;在出口处采用一阶外推;远场条件为无反射边界条件。图9为乘波构型流场压强和温度的分布图,图中锥形激波清晰明显;压缩角θc值大约为5°左右,与设计值非常接近。以已知的流场数据为基础,计算介质的吸收系数。在这里只考虑了H2O和CO2两种组分气体,设H2O的体积比为2%,CO2的体积比为0.033%。图10为乘波构型流场空气介质的收系数分布。在已知流场数据和介质吸收系数数据的基础上,计算乘波构型周围的辐射强度场。在计算中,假设飞行器壁面在3～5μm波段内平均当量发射率为εk,w=0.8;空间角度离散数量为N(N+2)=8(8+2)=80。下面给出本次计算的部分计算结果。乘波构型在3～5μm波段内的辐射热流密度q(W/m2)的空间分布如图11所示;图12和图13分别为截取的辐射热流密度空间分布的侧视图、后视图。由图11、12、13可以看出:①乘波构型红外辐射场的空间分布基本上为一“锥形”,这与锥导乘波构型流场特征相一致;②在乘波构型上、下表面附近的空间区域,红外辐射分布强烈,这是由于粘性摩擦使物面附近空气介质具有较高的温度所致;同时由于下表面还存在强激波,空气温度相对更高,最高可达约700K,因而下方的红外辐射更为强烈,传输得也更远;③乘波构型前向迎头方向上的红外辐射比较微弱,很快就与环境基本上没有差别。这是符合物理实际的,因为在高超声速状态下,飞行器对波前空间的流场是没有影响的。综合以上分析可知,仿真结果完全能够反映出乘波构型红外辐射特性的分布规律。6辐射模型的建立本文以红