发动机热喷流红外辐射计算与仿真

发动机热喷流红外辐射计算与仿真李建勋，童中翔*，王超哲，童奇，李贺，张志波【摘要】摘要利用求每个小视场视线方向辐射亮度的方法计算喷流红外辐射的光谱分布。以辐射传递方程数值和形式为基础，采用Malkmus统计窄谱带模型和Curtis-Godson(CG)近似求视线方向的辐射强度。采用CFD分析软件FLUENT模拟流场和组分摩尔分数分布。建立喷流红外成像仿真模型，仿真生成了液体火箭发动机热喷流红外图像。结果表明，该方法可以很好地分辨出流场的细微结构。该模型也适用于航空发动机喷流红外辐射计算与仿真。【期刊名称】光谱学与光谱分析【年(卷)，期】2013(033)001【总页数】7【关键词】关键词热学；发动机；红外辐射；热喷流；红外图像引言热喷流是飞机后向主要红外辐射来源之一，研究其红外辐射强度和光谱辐射强度，对飞机红外抑制技术、空空导弹及弹道导弹预警系统的红外探测系统的优化设计有重要意义。鉴于发动机热喷流的红外辐射特性研究的重要军事意义，国内外卜许多研究单位都开展了相关研究。国外已经建立了比较成熟的红外辐射特性计算软件［1］。国内的学者也在这方面进行了大量研究：邓洪伟等［2］进行了发动机喷管红外辐射数值模拟和实验研究，潘丞雄等［3,4］计算了发动机尾喷管的红外辐射特性，宗靖国等［5,6］研究了发动机排气系统红外特征。本工作以发动机热喷流为对象，研究热喷流红外辐射强度的光谱分布，建立热喷流红外辐射特性的通用计算模型，开发一种发动机热喷流红外辐射计算程序，实现了对红外热像仪的红外图像的仿真。1建模计算1.1流场计算喷流辐射取决于燃气的温度、压力和CO2，H2O等成分的组分浓度参数。因此，流场和组分浓度的计算是必要的。对于亚音速圆喷管热喷流，可以采用CFD和理想湍流射流模型两种方法。FLUENT商用软件［7］,是目前功能最全面、适用性最强、使用最广泛的CFD软件之一。其求解器基于有限体积法（finitevolumemethod,FVM）,包含丰富的物理模型，可以精确模拟无粘流、层流、湍流、传热和传质以及多相流等复杂的流动现象，为用户提供了二次开发接口，可开展从低速、跨音速到超声速和高超声速流动模拟。理想湍流射流模型［8］认为亚音速圆喷管热喷流的截面包括主体段和初始段部分［9］。1.2视线方向气体辐射亮度计算1.2.1吸收性介质辐射传递方程设定向光谱辐射亮度为Lv的辐射能在吸收/发射性介质内传递，按照布格尔定律，该能量通过厚度为ds的气体后，在原射线方向上辐射能的衰减量dLv正比于投射量Lv及厚度ds,即如果沿途介质处于局部热力学平衡之中，根据基尔霍夫定律，某方向厚度为ds的介质沿着该方向发射的光谱辐射亮度为考虑到吸收和发射，由式（1）和（2）可知，该方向光谱辐射亮度的变化量是式（3）即为吸收/发射性介质辐射传递微分方程。该式等价于X（S）=为光学厚度。式（4）两边乘以exp（X）从0积分到X从而有式（8）表明：定向光谱辐射亮度Lv（0）经过光学厚度为X（S）的吸收和发射性介质后，最终辐射亮度是两项辐射亮度的叠加，第一项是Lv（0）经过光学厚度X（L）,被衰减后的辐射亮度；第二项是X*处气体介质发射的辐射亮度Lvb（X*）经X-X*光学厚度的气体衰减后的光谱辐射亮度。式（8）是一个积分方程，为便于计算，还应推导它的数值和形式。1.2.2辐射传递方程数值和形式对式（1）从Lv（0）到Lv（S）进行积分因为X（S）,经计算，化简得气体透过率的定义是设气体介质厚度为L,以观测端为起点，辐射传递的起始位置为m=N,将传递路径分为N小段，如图1所示。N的值足够大，因此在每段内气体压力，温度和组分浓度可认为是相等的。则（8）式可以写成式中X（m）表示从1到m这段气体的光学厚度。将式（11）代入式（12）,得到吸收/发射性介质辐射传递的数值和形式在该式中，辐射传递方向从1到NoTq（m）表示从点m到点N的透过率。数值和形式利用透过率来求解气体辐射传递方程，便于计算机编程。对光谱辐射亮度在波长域或频率域进行积分可以得到波段积分辐射强度。1.3谱带模型利用辐射传递方程式（13）计算视线方向光谱辐射亮度，要用到光谱吸收系数，光谱吸收系数是各种分子所有谱线在该波长的吸收贡献之和，与温度和压力（浓度）等有关，需要采用逐线法来计算，计算量非常大。在谱带间隔内，吸收系数的变化相对普朗克函数变化很快，将普朗克函数看作常数，利用谱带平均透过率求谱带平均辐射亮度，计算量将大幅减小，也能得到比较准确计算结果。谱带模型按波数间隔的大小分为窄谱带模型和宽谱带模型。窄谱带模型［10］的波数间隔一般为5~50cm-1,而宽谱带模型则可包括整个振动-旋转谱带。根据谱线分布规律的假设，谱带模型分为正规谱带模型和统计谱带模型。正规谱带模型认为气体的吸收带是由一系列等谱线强度、等谱线半宽、等谱线间隔排列的洛伦兹谱线组成。正规谱带模型适用于谱线分布规则的气体分子谱带辐射特性计算。统计谱带模型认为：气体分子在波数间隔内分布着N条谱线，这些谱线的位置、强度的分布均是随机的，即认为在△n内任意位置上谱线出现的概率是相同的。按照谱线强度分布规律假设的类型不同，统计谱带模型可分为等线强度分布，指数线强度分布（Goody），指数尾倒数线强度分布（Malkmus）的统计模型。Malkmus模型在大气科学研究中应用广泛，后来也被气体高温辐射计算采用。对多原子分子气体Malkmus模型被公认为最佳谱带模型，当谱带宽度为25cm-1时和逐线计算的吸收率间的最大误差为10%，如果用最小二乘法从分子谱线库推导模型参数，在带宽为10cm-1、气体压力大于0.1个大气压时，Malkmus窄带模型的误差小于1%。在很多工程问题中，涉及非均匀气体辐射传递。因此，由均匀气体推导的谱带模型处理非均匀气体，要采用近似方法，即按照非均匀路径的参数分布计算出〃等效均匀路径”的参数，再采用谱带模型来计算。以辐射传递方程数值和形式为基础，本文采用Malkmus统计窄谱带模型［11］和非均匀气体Curtis-Godson近似［12］来求视线方向的辐射亮度。2喷流红外成像仿真2.1喷流红外图像仿真方法红外图像反映探测阵列中每个像素的辐射照度分布，辐射照度与瞬时视场内的辐射亮度分布是一样的，因此只要计算出所有像素对应瞬时视场的辐射亮度，并将其归一化，按照伪彩色编码关系，就能模拟红外热像仪拍摄喷流形成的热像。红外图像仿真流程如图2所示。2.2瞬时视场视线辐射传递路径瞬时视场气体辐射传递路径指的是瞬时视场中心的视线穿过喷流包络体的路径，由于瞬时视场很小，该传递路径的辐射亮度看作是瞬时视场内的平均辐射亮度。将CFD的计算域作为辐射源（完全包裹热喷流的圆柱形包络体），瞬时视场辐射传递路径起点与终点坐标由焦平面阵列探测元视线与圆柱包络体的表面交点确定。设红外热像仪镜头中心和喷流水平中心在同一面上，镜头外法线平行于喷流正侧向。红外热像仪视场（FOV）为ax&,其空间分辨率是mxn（即探测阵列有m行n列的探测元）。探测器位于焦平面上，设探测器阵列中心的坐标是（0,0,f）,喷流与热像仪的相对位置如图3所示。第i行，第j列的像素中心的坐标是系数a',b',c‘的表达式是设从像素中心经透镜中心引出视线的方向向量为（a,b,c）该视线方程也可写成如下形式其中t是形参。而圆柱包络面的方程是联立方程得若该方程有两个实数根t1，2,将两个根分别代入式（17）,得到气体辐射传递路径起点和终点坐标（xA，yA，zA）和（xB，yB，zB）。如果无根，就说明该探测面视线与包络面无交点，接收的是背景辐射。2.3辐射量伪彩色编码设红外热像仪的动态响应范围是E0min到E0max，探测器焦平面接收到的辐照度变化范围使探测元输出在响应范围之内，那么以下关系式应成立E0min<Emin<Emax<E0maxo计算探测器阵列各元素的灰色值其中，color_number是灰色级别数，表示从最强到最弱分成的等份数。由于计算机能够显示的是8bit的视频格式，所以范围一般从0到255，可以取color_number为255。下面采用三原色（RGB）模式来研究灰度图与伪彩色图生成的方法。可见光包含了在一定范围内的不同波长的电磁波，但并没有必要对每一种波长的可见光都用一种颜色来独立表示。因为几乎所有的颜色都可以用红、绿、蓝（RGB）这三种颜色的不同强度来组合而得。用I这个无量纲量来表示红外热像仪探测阵列各探测元接收的能量强弱的不同，再将I按一定函数关系映射成相应的彩色，不同的灰度I对应不同的彩色，因此图像的伪彩色生成首先要设定红、绿、蓝三个变色函数，对应每一个I都有相应的RGB输出，将这种变换关系称为颜色编码［13L3算例3.1液体火箭发动机喷流3.1.1流场参数采用文献［14］中的算例检验本算法准确性，文献中使用aerodyneradiationcode（ARC）计算了一台推力为1.125x106N的液体火箭发动机在2-5pm的红外辐射特性，发动机氧化剂为N2O4，还原剂为单甲基肼［14］,考虑喷流中H2O和CO2等组分的吸收和发射。将圆柱面包住喷焰，圆柱的半径是11m,长度是112m。整个喷流的静压保持为1个大气压不变，温度和组分摩尔分数随轴向和径向位置变化，如表1所示。3.1.2计算方法液体火箭发动机喷流不含固体颗粒，因此认为没有散射，并且认为它处于化学和热力学平衡状态。根据表中的流场参数计算正侧面（方位角为90°）喷流红外卜光谱辐射强度，所用的方法是：首先，建立喷流的圆柱包络面，划分其表面网格，计算出各面元面积Ai和法线方向。其中，cosei是面元发现与观测方向夹角的余弦。将光谱辐射强度对波长积分就能计算出某波段的辐射强度（W-Sr-1）。3.1.3计算结果图5是文献［14］计算的液体火箭喷焰沿正侧向的光谱辐射强度，计算过程可参见文献［14］,图6是本方法计算结果。比较两图可以看出，本计算光谱辐射强度及随波长的变化规按照观测方向，判断在观测方向可见的面元。对于每个可见面元，以面元中心为起点，按照观测方向弓I出射线与喷流相交得到气体辐射传递路径，计算出该路径的辐射亮度日，该辐射亮度的计算方法与1.2节相同。通过下式计算出该观测方向的光谱辐射强度［W-(Sr-pm)-1］律与文献计算的结果符合较好。但是在2.5-3pm波段和4.6pm处两种计算结果存在较大差异，分析原因可能是由于谱带模型参数不同所导致的。3.2热喷流伪红外图像仿真3.2.1流场计算设置计算域如图7所示。喷管内部从涡轮后开始建模，包括内外涵、波瓣混合器和中心锥等部件。其尺寸来自某型发动机的测量结果。以某型发动机最大状态作为喷管计算工况，其进口参数如表2所示。内外涵进口燃气为亚音速，设定为压力入口，给定进口总压、静压、总温。出口截面条件由流动状态决定，若为亚音速流，则给定反压，其余变量通过外插给出；若为超音速流，则外插给出所有流动参数。壁面采用绝热、无滑移物面条件。对称边界采用无流量通过条件，法向速度为0，其他参数外推。对于外流的入口，采用黎曼边界条件。空气和燃气的组分参数如表3所示。最大状态时喷管喉道半径R8=29.3cm，出口半径R9=30.6cm。其中喷流部分长度和半径为：1217cmx400cm，网格节点的数目是151x111(轴向，径向)。流场参数预处理：喷流辐射模拟需要的流场参数包括静压、静温、水蒸汽和CO2摩尔分数。热喷流网格如图8所示，划分网格时做了局部加密处理以适应实际流场参数的变化。这种非均匀网格的流场在使用时很不方便，先对该流场进行预处理，插值计算出新的均匀网格的流场作为辐射计算的输入条件。这里新的流场节点数目是151x120。3.2.2喷流流场图9所示为喷流静压等高线分布，内外涵气体掺混后经喷管加速膨胀后喷管的出口压力为104608Pa，由于大气压力为95930Pa，因此该喷管处于不完全膨胀状态，从图中可以看到膨胀波和压缩波组。图10所示为喷管出口局部静压等高线，可见由于喷管出口压力大于外界大气压，在出口边缘处产生两道膨胀波，气流经过后气压下降并与外界气流压力相等。这两道膨胀波相交后在交点处继续产生膨胀波与自由边界相交，气流经过后压力继续下降，以致压力稍低于外界大气，这时外界气体将压缩气流，迫使气流向内转折。由于压力相差较少，超音速气流向内转折时，在转折点处将产生弱压缩波。这两道压缩波相交后，会在交点处再产生两道压缩波。此后气流压力又稍大于外界大气压力，气流将继续膨胀，再度交替产生膨胀波和压缩波组，但强度逐渐衰减，直至消失。由图可知，在该工况下喷流在出现第二道膨胀波时，强度已经大为衰减。喷流的温度分布以及H2O,CO2的摩尔分数分布如图11—图13所示。3.2.3红外仿真图像红外热像仪的总视场(FOV):20°x15°;空间分辨率：320x240;波段：2000~5000cm-1。喷流横截面的尺寸是800cmx1217cm。与喷流中心在同一水平面，考虑了两种热像仪位置：使喷流计算域轴向恰好充满视场，距喷流中轴Z0=(1217/2)/tan(10/180xpi)=3451cm，热像仪距喷管出口Y0=608.5cm。(2)距喷流中轴Z0=1300cm，距喷管出口Y0=Z0xtan(10/180xpi)=229.23cm。如图15和图16所示，热像仪离得越近，喷流高温区所占的热像仪像素点越多，越能清楚的描述高温区的辐亮度变化。从图中可以看到，喷管出口(左侧)及下游不远处分别出现了辐亮度明显大于周边的球形区域，而连接两球形区域的地方，高辐射亮度区域向内收缩形成颈部。由于该部分喷流的温度和红外活性气体的摩尔分数几乎相等，产生该现象的原因是该处存在交替出现的膨胀波和压缩波，因此导致压力高低交替变化，在高压部位与附近低压区域红外活性气体的摩尔分数几乎相同，但是分子密度更大，因此同样体积分子数越多，红外辐射越大。4结论以发动机热喷流为对象，研究了发动机热喷流红外辐射特性和热喷流伪红外图像的生成方法。分析研究得出以下结论(1)由均匀气体推导的谱带模型处理非均匀气体，要采用近似方法，即按照非均匀路径的参数分布计算出〃等效均匀路径”的参数，再采用谱带模型来计算。提出的视线方向光谱辐射强度计算结果与文献计算结果基本一致。但在2.5-3pm波段和4.6pm处本计算结果与文献计算结果存在较大差异，可能是谱带模型参数不同导致的。发动机喷管出口正后向气流压力与外界大气压力不相等，引起膨胀波和压缩波组交替出现，但强度逐渐衰减，直至消失。热喷流伪红外图像中，喷管出口及下游不远处分别出现辐亮度明显大于周边的球形区域，而连接两球形区域的地方，高辐射亮度区域向内收缩形成颈部。红外仿真图中可以清晰地分辨出流场的细微结构。References[1]Murtagh,DonalP.35thInternationalConferenceonInfrared,Millimeter,andTerahertzWaves,IRMMW-THz2010,Rome,Italy.DocumentNo.20104813432516,2010,9:445.[2]DENGHong-wei,SHAOWan-ren,ZHOUSheng-tian,etal（邓洪伟，邵万仁，周胜田，等）.Aeroengine（航空发动机），2010,36（1）:44.[3]PANCheng-xiong,ZHANGJing-zhou,SHAOWan-ren,etal（潘丞雄，张靖周，邵万仁，等）.JournalofAerospacePower（航空动力学报），2010,25（7）:1518.[4]LIUChang-chun,JIHong-hu,LINa,etal（刘常春，吉洪湖，李娜，等）.JournalofEngineeringThermophysics（工程热物理学报），2010,31（9）:1567.[5]ZONGJing-guo,ZHANGJian-qi,LIUDe-lian（宗靖国，张建奇，刘德连）.ActaPhotonicaSinica（光子学报）,2011,40（2）:289.[6]WEIJun-guang,YANGQing-zhen,LIYue-feng（未军光，杨青真，李岳锋）.ComputerSimulation（计算机仿真）,2011,28（4）:66.[7]ZHOUZhang-gen,MADe-yi（周章根，马德毅）.MachineBuildingandAutomation（机械制造与研究），2010,39（1）:61.[8]WUHui-xia