轨控系统热喷干扰效应数值模拟

2022-10-14傅杨奥骁丁明松刘庆宗江涛石润董维中高铁锁

航空学报 2022年9期

傅杨奥骁,丁明松,2,刘庆宗,2,*,江涛,石润,董维中,高铁锁

1. 中国空气动力研究与发展中心计算空气动力研究所,绵阳 621000 2. 中国空气动力研究与发展中心空气动力学国家重点实验室,绵阳 621000

近年来,高超声速飞行器工作区间呈现出跨流域、宽速域的趋势,常用的气动舵面控制技术常常不能满足跨流域高速飞行器的控制和机动需求,喷流反作用控制系统与气动舵面的复合控制逐渐成为高机动操控的必要手段之一。

喷流复合控制方法有姿控(力矩操纵)、轨控(力操纵)2种模式,其中轨控喷流一般具有推力大、流量大的特征,其与来流的相互干扰更加强烈、影响范围更大。轨控喷流反作用控制系统一般采用小型火箭发动机,喷管出口气体为高温燃气,其温度较高,一般称为“热喷”。而在风洞试验中,常采用空气作为喷射气体进行喷流试验模拟,由于其温度相对较低,一般称为“冷喷”。由于冷喷试验重复性好、无毒害、成本低,因此在风洞试验中被大量采用。这就存在热喷和冷喷差异性及其影响问题：一方面,热喷高温燃气与空气冷喷的气体性质存在较大差别,影响喷流与外流场相互干扰的流动结构,即异质气体喷射效应；另一方面,热喷燃气中常含有未反应完全的易燃组分,进入高速主流后会发生再次燃烧现象,即二次燃烧效应。这二者可统称为热喷干扰效应。

与喷流干扰相关的空气动力学问题研究在过去几十年引起了广泛的兴趣,但长期以来的研究多集中于姿控喷流干扰以及冷喷的试验和数值计算,而热喷干扰效应则是冷喷试验/数值计算难以模拟的。目前,对热喷干扰问题,已有学者开展了一些研究,但在这些研究中,大多对燃气喷流采用等效比热比异质流模型或者化学反应冻结流模型进行计算,未考虑喷流燃气进入外流后的二次燃烧效应。例如：Votta等通过数值计算研究了火星着陆器的喷流干扰问题,在计算中将喷流气体取为不反应的H、N和NH混合气体；Despirito通过数值计算研究了冷、热喷干扰问题中湍流模型的影响,在计算中对热喷流气体同样取为不反应的燃气混合气体；Gnemmi和Seiler通过风洞试验和数值计算研究了尖锥外形的喷流干扰问题,在计算中采用了等效比热比异质流模型；孙得川等通过数值计算研究了喷流气体性质对侧喷干扰流场的影响,在计算中同样对热喷燃气取为化学反应冻结的混合气体。

部分学者开展了热喷二次燃烧效应的探讨,但这些研究仍不够充分。例如,Ebrahimi研究了导弹轨控系统热喷干扰效应的影响,结果显示燃气的二次燃烧效应对整体气动力特性有较大影响,但这项研究只计算了单一飞行状态,且仅限于对气动力特性的影响研究；Dong等对冷、热喷干扰风洞试验中的表面压力分布进行了对比验证,结果显示,冻结化学反应时,热喷干扰分离区长度明显低于试验结果,考虑二次燃烧效应时,分离区长度与试验接近,但这项研究仅针对风洞试验和数值模拟的对比结果给出了定性结论。

可见,目前国内外对于热喷干扰问题的研究还相对不足：当前的研究大多采用化学反应冻结流或等效比热比异质流模型,未考虑喷流燃气进入外流后的二次燃烧效应；或者仅给出了少量状态的气动力方面的定性结论,缺乏足够的研究结果来总结规律。由于高温燃气热喷的二次燃烧效应和异质气体喷射效应会受到飞行速度、高度、姿态等因素影响,其作用大小和影响规律尚不十分明确,同时高温燃气喷流对喷口附近局部热环境的影响不容忽视,因此有必要进一步开展多种条件下热喷干扰效应作用机制及其对气动力热特性影响规律的研究。

为了真实地反映高温燃气喷流干扰特性,针对轨控系统热喷干扰问题,开展典型外形热喷干扰化学非平衡流场的数值模拟,研究二次燃烧效应和异质气体喷射效应对飞行器气动力热特性影响的机理,分析不同马赫数、飞行高度和攻角条件下热喷干扰效应对飞行器气动力特性和表面气动热环境的影响规律。

1 数值计算方法

高温燃气进入主流后会发生二次燃烧现象,这些化学反应是在流动中发生的,这就存在化学反应的非平衡现象。为了分析高温燃气热喷干扰效应,控制方程采用完全气体和化学非平衡流的统一形式,其无量纲化为

(1)

式中：为守恒变量,完全气体模型为=[,,,,],化学非平衡气体模型为=[,,,,,];为混合气体总密度；为组分的密度；、、分别为、、方向的速度；为气体内能；为雷诺数,、、分别为、、方向的对流项；、、分别为、、方向的黏性项；为化学非平衡源项。

采用结构网格的有限差分方法离散控制方程式(1),对流项采用AUSMPW+(Advection Upstream Splitting Method by Pressure-based Weight Functions)格式离散,黏性项采用中心差分格式离散,时间离散采用LU-SGS(Lower-Upper Symmetric Gauss Seidel)隐式方法,湍流模型采用Menter SST(Shear Stress Transport)两方程模型,具体处理方法详见文献[12,14-15]。

2 气体模型

2.1 二次燃烧效应模拟方法

由于喷流燃气中的颗粒相杂质占比很小,对干扰流场的影响很小,因此这里忽略燃气喷流中少量固体杂质的影响,主要考虑气体介质的化学反应。

发动机产生的高温燃气进入空气后,需要考虑空气受热离解、复合、置换等化学反应,如表1所示,M为碰撞体。由于常用的火箭发动机推进剂一般包含C、H、O、N等元素,高温燃气中包括CO、C等易燃组分。对于含C、O等元素气体,需要考虑与CO和CO相关的7个反应,如表2所示。高温燃气中还可能包括OH、H、H等易燃组分,因此还需要考虑了HO、H相关的8个反应,如表3所示。

表1 高温空气化学反应模型Table 1 Chemical reaction model of high temperature air

表2 含C元素组分化学反应模型

表3 含H元素组分化学反应模型

根据燃气成分的不同,将表1～表3的化学反应进行组合,形成高温燃气与空气混合的燃烧化学反应模型,反应速率通过Arrhenius公式计算

(2)

式中：为化学反应速率；为温度；、、为系数常数。表1～表3中各个反应的系数常数可参考文献[17-22]。

燃烧过程中,化学非平衡效应的模拟通过流动方程与化学反应源项的强耦合自动实现,同时考虑组分的束缚电子激发效应和分子组分振动能量激发效应；湍流燃烧效应的模拟,采用耦合计算Menter SST湍流模型实现。相关处理方法详见文献[12,14-15]。

为了分析二次燃烧效应的影响,本文将采用无反应喷流模型进与高温反应喷流模型行对比计算,无反应喷流模型冻结了燃气组分相关的化学反应,即将喷流燃气视为不反应的多组分化学惰性气体,其余处理方法与高温反应喷流一致。

2.2 异质气体喷射效应模拟方法

在喷流干扰的风洞试验和数值模拟中,为了简化问题、降低成本,常将空气作为喷流气体,但实际喷管出口一般为高温燃气,其气体性质与空气较大差别。为了研究这种异质气体喷射效应的影响,本文还计算空气冷喷流状态,计算中采用完全气体模型。通过燃气热喷流(不考虑二次燃烧)与空气冷喷流的对比,即可单独分析异质气体喷射效应的影响。

为了保证空气冷喷流对燃气热喷流的正确模拟,需要通过相似准则对出口参数进行匹配：首先,冷、热喷流的喷管出口形状与面积保持一致；从以往的研究中可知,喷流干扰流场结构主要受压力比和动量比的影响,即式(3)～式(4)；但是仅凭这样仍无法确定空气冷喷流出口的密度和温度,采用与文献[9,23]中相同的做法,即保证冷、热喷流的质量流量一致,在已满足式(3)～式(4) 的条件下即是保证冷、热喷流的密度一致。

(3)

(4)

(5)

式中：为静压；为速度；为气体常数；下标j表示喷管出口,air表示空气冷喷流,gas表示燃气热喷流。对于量热完全气体,内能只是温度的线性关系式,而输运系数采用Sutherland公式。

2.3 喷流干扰气动力特性计算方法

工程上通常采用喷流干扰附加力、附加力矩及放大因子来衡量喷流干扰效应对飞行器气动力特性的影响。定义无喷时气动力、力矩,喷流开启时气动力、力矩,则喷流干扰附加气动力、力矩可表示为

(6)

轨控喷管位于质心,喷流沿法向喷出,因此定义喷流推力放大因子为

(7)

式中：为喷管推力。

在本文的计算中,力矩参考点均取质心位置,俯仰力矩均以抬头为正。

3 计算方法的验证

为了验证计算方法的正确性,这里针对锥柱裙外形侧向喷流试验进行数值模拟,并与试验结果进行对比。风洞试验模型外形如图1所示,其中导弹身部直径=90 mm,喷管出口直径4.6 mm。试验中风洞来流为空气,来流静压92.3 kPa,来流静温105 K,来流马赫数3.0,攻角0°；喷流总压12.0 MPa,喷流压比130,冷喷时喷流气体为空气,喷流总温293 K,喷流出口温度244 K；热喷时喷流气体为燃气,喷流总温2 300 K,喷流出口温度2 058 K,喷流出口比热比1.233,喷口处组分及质量分数为：CO(38.21%),HO(10.47%),CO(35.47%),N(14.12%),H(1.73%)。

壁面处取绝热壁、无滑移、零压力梯度、非催化条件,冷喷计算时采用用完全气体模型(比热比1.4),热喷计算时,采用前述的考虑二次燃烧效应的化学非平衡流模型。

图1 试验模型外形Fig.1 Configuration of test model

为了保证计算结果的网格无关性,这里对疏、密2套网格进行对比计算,其中密网格Grid1壁面法向第1层网格高度为10m,疏网格Grid2壁面法向第1层网格高度5×10m,计算状态为冷喷试验状态。图2(a)给出了计算网格Grid1的示意图,图2(b)给出了Grid1和Grid2计算得到的子午线上压力系数分布对比,可见,二者计算结果差异很小,满足网格无关要求。

图2 网格示意图及不同网格得到的压力系数分布对比Fig.2 Schematic of grid and comparison of pressure coefficient distribution

图3给出了冷喷时子午线上表面压差系数(Δ=,jet-,nojet,其中,jet为喷流开启时压力系数,,nojet为无喷时压力系数)分布与文献中试验结果的对比,可以看出,本文冷喷计算结果与试验数据基本吻合。图4给出了热喷流场的CO质量分数分布。图5给出了热喷试验中子午线上表面压差系数分布与文献中试验结果的对比,热喷计算结果与试验数据基本吻合。从对比结果来看,本文的数值计算方法适用于冷、热喷干扰流场的计算。

图3 冷喷模型表面压差系数分布结果对比Fig.3 Comparison of pressure coefficient distribution along model surface for cold jet test

图4 热喷流场CO2质量分数分布云图Fig.4 Contour of CO2 mass fraction in flow field for hot jet test

图5 热喷模型表面压差系数分布结果对比Fig.5 Comparison of pressure coefficient distribution along model surface for hot jet test

4 轨控系统热喷干扰效应影响分析

4.1 计算状态、网格及典型状态流场结构

针对一种锥柱裙拦截弹外形开展控系统热喷干扰效应的数值模拟,如图6所示,弹体总长2 250 mm,身部直径为386 mm,底部最大直径为502 mm,轨控发动机位于质心,距头部1 177 mm。轨控发动机为工程应用中的典型小型火箭发动机,其推进剂为NO/MMH(一甲基肼),燃烧室总温3 200 K,总压4 MPa,总推力为6 800 N,喷管出口直径78 mm。喷管出口燃气组分及质量分数为CO(8.67%)、HO(29.76%)、CO(17.44%)、N(41.91%)、H(1.52%)、其他(0.7%)；喷管出口马赫数2.796,出口静压0.133 3 MPa,静温1 596 K。

计算状态为高度=20,30,40,50 km,马赫数=2～10,攻角=-10°～10°。在计算中,采用Menter SST湍流模型,对壁面取等温壁300 K、无滑移、零压力梯度条件,不考虑壁面催化效应,喷口处取为喷管出口条件。

图6 锥柱裙拦截弹外形Fig.6 Sketch of configuration of interceptor missile

为了对网格收敛性进行研究,这里分别取4套网格计算,如表4所示,主要对网格规模、第1层网格高度进行对比。图7(a)给出了计算网格示意图,图7(b)给出了喷口前方高热流区域的热流分布对比,其中为表面热流值。表5给出了气动力系数对比,其中为轴向力系数,为法向力系数,为俯仰力矩系数。可以看出：① Grid2 与Grid4第1层网格高度相同,但Grid4网格总量更大,从结果来看,二者的气动力热特性结果差别较小,说明在当前状态下,网格规模对计算结果的影响不大；② Grid1与Grid4的网格总量一致,但Grid1的第1层网格高度更小,从结果来看,二者的热流峰值存在差别,说明在当前状态和网格规模下,第1层网格高度的影响较大；③ Grid3 较Grid1网格规模更大、第1层网格高度更小,从结果来看,二者差异较小,满足网格无关性要求,为了兼顾计算效率和精度,本文采用Grid1进行计算。

表4 不同网格对比Table 4 Comparison of different grids

图7 网格示意图及不同网格得到的热流分布对比Fig.7 Schematic of grid and comparison of heat flux

图8给出了=20 km、=5、攻角=0°条件下轨控系统热喷干扰流场的结构示意图,图9 给出了喷口附近的表面极限流线分布。从图中可以看到,由轨控发动机喷管产生的喷流与高速来流相互作用,形成了弓形激波、桶型激波、马赫盘、尾迹涡等典型的波系结构和涡结构。在喷口前方可以看见一条明显的主分离线和再附线,在喷口下游,可以看见由喷流尾迹形成的另一组分离线和再附线。

表5 不同网格计算得到的气动力系数对比

图8 干扰流场结构示意图Fig.8 Schematic of jet interaction flow flied

图9 喷口附近极限流线分布Fig.9 Distribution of surface streamline near nozzle exit

4.2 二次燃烧效应的影响分析

图10为高温反应喷流(Reaction)与无反应喷流(No Reaction)得到的CO和HO质量分数分布对比,其计算状态为=20 km、=5、=0°。由于CO和HO是二次燃烧反应的主要产物,因此通过其质量分数的对比可以反映出二次燃烧反应的强度,从图8中可见,高温反应喷流得到的喷口上游回流区略大于无反应喷流,流场中CO和HO的质量分数也明显更高,尤其是在弓形激波附近、喷口周围以及喷口下游二次分离/再附线附近,这表明高温燃气进入主流后,在这些位置发生了二次燃烧反应,导致CO和HO的质量分数增大。

图10 反应喷流与无反应喷流的组分质量分数分布对比Fig.10 Comparison of species mass fraction distribution between reacting and non-reacting flow

图11给出了流场中化学反应导致的化学能密度变化,其计算方法为

(8)

图11 流场中化学反应导致的化学能密度变化Fig.11 Chemical energy density variation caused by chemical reaction in flow field

表6为高温反应喷流与无反应喷流计算得到的气动力特性对比。可见,考虑流场中的二次燃烧效应后,轴向力系数变化不大,但喷流附加推力增大,附加低头力矩减小,从而导致法向力系数绝对值增大、俯仰力矩系数绝对值减小、推力放大因子增大。

表6 二次燃烧效应引起的气动力特性变化

图12给出了流场对称面及弹体表面的压力分布对比,图13进一步给出了弹体表面不同位置的压力分布对比,图中表示沿周向角度。可以看出,考虑二次燃烧效应后,喷流干扰区面积增大,且干扰区内压力明显升高,高温反应喷流较无反应喷流的干扰区流场结构整体呈现出更加膨胀的特点,结合图10可知,考虑化学反应后在干扰区内CO和HO的质量分数明显增大,而生成这2种组分的反应均为放热反应,释放的化学能使周围流场温度升高,在同等条件下温度升高会导致气体膨胀,从而使干扰区范围变大、干扰区内压力值升高。由于干扰区大多集中于弹体上半部分,整体的压力积分效果呈现为一个向下的附加法向力,与喷管推力同向,因此高温燃气二次燃烧效应的影响表现为增加了附加推力。同时由于喷口下游二次分离/再附线附近的二次燃烧反应增强(见图10(b)、图10(c))导致这一区域压力明显升高(见图13(b)),这一部分的压力积分呈现为1个抬头力矩,因此高温燃气二次燃烧效应的影响还表现为增加了1个抬头的附加力矩。

图12 二次燃烧效应引起的对称面及表面压力分布变化Fig.12 Variation of pressure distribution on symmetry plane and surface caused by aferburning effect

图13 二次燃烧效应引起的表面压力分布变化(H=20 km, Ma=5, α=0°)Fig.13 Variation of surface pressure distribution on surface caused by afterburning effect (H=20 km, Ma=5, α=0°)

图14为高温反应喷流与无反应喷流的热流分布云图对比,图15为上表面子午线上的热流分布对比。可见,壁面热流极值位置与压力极值基本重合,均为喷口上游的再附点附近,二次燃烧效应使热流峰值升高了21%左右。图16给出了表面化学能密度变化,可见,化学能释放剧烈的区域与高热流区基本重合,结合图8可知,这些区域均是生成CO和HO反应较剧烈的区域,说明二次燃烧反应释放的能量造成了干扰区表面热流值升高。

图14 热流分布云图对比Fig.14 Comparison of heat flux contour

图15 二次燃烧效应引起的上表面子午线上热流分布变化Fig.15 Variation of heat flux distribution along upper side meridian line caused by afterburning effect

图16 化学反应导致的化学能密度变化Fig.16 Chemical energy density variation caused by chemical reaction

由于二次燃烧效应会受到飞行速度、高度、飞行姿态等因素的影响,因此,接下来将开展不同飞行条件下二次燃烧效应的影响分析。

图17给出了不同条件下推力放大因子的对比。可以看出,高度越低、马赫数越高,二次燃烧效应对气动力特性的影响越大,这符合一般的认知：随着高度增加,来流密度迅速降低,因此来流氧气质量流量迅速降低(如图18所示,其中,为来流密度；O为来流中氧气的质量流量；为喷流燃气的质量流量),当来流氧气质量流量过低时,不足以使二次燃烧反应充分,会导致流场中二次燃烧效应明显减弱；马赫数越高,喷流与来流相互干扰作用越强,类似于激波诱导燃烧,在强干扰的影响下会诱导发生更强的化学反应,因此二次燃烧效应的影响增强。在负攻角条件下,二次燃烧效应的影响略大于正攻角,这是由于负攻角时,喷口位于迎风面,由于受到来流的压迫增强,激波干扰强度增加,导致二次燃烧效应增强,正攻角时情况则相反,喷口位于背风面,导致二次燃烧效应有所减弱。最后值得指出的是,高度30 km以上时,二次燃烧效应对气动力特性的影响已经较小,说明在本文的高空状态下,可以不考虑二次燃烧效应的影响。

图17 不同条件下二次燃烧效应对推力放大因子的影响Fig.17 Influence of afterburning effect on thrust amplification factor under different conditions

图18 来流密度和氧气质量流量随高度变化规律Fig.18 Inflow density and O2 mass flux for different altitude

图19给出了不同条件下干扰区热流峰值的对比。可见,二次燃烧效应对气动热环境的影响规律与气动力特性是类似的：高度越低、马赫数越高、喷口越迎风,二次燃烧效应的影响越大。值得注意的是,攻角变化时,相对于气动力特性,二次燃烧效应对表面气动热环境的影响更加显著。

总的来看,高温燃气喷流的二次燃烧效应对飞行器气动力热特性有较大影响,表现为增加附加推力、增加抬头的附加力矩(以本文设置的力矩参考点),以及使干扰区热流峰值上升。同时,高度越低、马赫数越高、喷口越迎风,二次燃烧效应的影响越大。

4.3 异质气体喷射效应的影响分析

喷流燃气的气体性质与空气差别较大,如表6 所示,喷流燃气的气体常数为400 J/(kg·K),比热比为1.27,这些物理参数的差别会对燃气喷出后的膨胀过程产生影响。在本节的计算中,燃气喷流采用真实的多组分燃气混合气体,但人为冻结了燃气进入外流后的化学反应,这样可以将真实燃气喷流中喷流介质热物性参数的影响提取出来。下面将通过燃气喷流(Gas Jet)与空气喷流(Air Jet)的对比,来分析喷流燃气的异质气体喷射效应。

表7给出了燃气喷流与空气喷流计算得到的气动力特性对比,计算状态为=20 km、=5、=0°。可以看出,喷流燃气的异质气体喷射效应对气动力特性的影响与二次燃烧效应类似：对轴向力影响不大,增加向下的附加法向力,增加抬头的附加力矩,但对比表6可以看出,其影响幅度明显小于二次燃烧效应。

表7 异质气体喷射效应引起的气动力特性变化Table 7 Variation of aerodynamic characteristics caused by heterogeneous gas injection effect

下面分析燃烧喷流与空气喷流气动力特性差异产生的原因：图20给出了在对称面及弹体表面的压力分布对比,图21给出了燃气喷流的比热比分布云图。可以看出,由于燃气喷流的比热比低于空气,可压缩性更强,因此其喷流干扰区面积略小于空气喷流。图22进一步给出了弹体表面的压力分布对比,可见,除再附点附近的压力峰值区域外,燃气喷流的整体压力分布与空气喷流较接近,高压区更多集中于弹体上表面,平均压力略有升高,因此使附加法向推力略有增大、附加低头力矩略有减小。但通过图20、图12,以及图22、图13 来对比2种效应引起的压力分布变化,可知燃气喷流与空气喷流二者压力分布的差别明显小于高温反应喷流与无反应喷流的差别,因此异质气体喷射效应对气动力特性的影响整体较小。

图20 异质气体喷射效应引起的对称面及表面压力分布变化Fig.20 Variation of pressure distribution on symmetry plane and surface caused by heterogeneous gas injection effect

图21 无反应燃气喷流比热比分布云图Fig.21 Specific heat ratio contour for non-reacting flow

图23为弹体上表面子午线上热流分布对比,与图15不同的是,二者的热流分布非常接近,热流峰值差别不大,表明异质气体喷射效应对表面热流分布的影响较小。

图23 异质气体喷射效应引起的上表面子午线上热流分布变化Fig.23 Variation of heat flux distribution along upper side meridian line caused by heterogeneous gas injection effect

综合来看,当前计算状态(=20 km、=5、=0°)下,喷流燃气的热喷干扰效应以二次燃烧效应为主,异质气体喷射效应的影响较小。由于异质气体喷射效应会受到飞行速度、高度、飞行姿态等因素的影响,因此,接下来将开展不同飞行条件下异质气体喷射效应的影响分析。

图24给出了不同条件下的推力放大因子对比,图25 给出了干扰区热流峰值对比。可见：① 在各个高度,异质气体喷射效应对气动力热特性的影响均是一个小量;② 当马赫数较低时,异质气体喷射效应对气动力热特性的影响较小,当马赫数增大时,影响程度有所增大。下面分析其原因：随着马赫数增加,激波干扰强度增加,造成激波层及干扰区内的温度升高,如图26所示,由于燃气喷流考虑了内能不同模态的激发以及热物性随温度的变化,从图21和图27可以看出,虽然喷口处燃气比热比是一致的,但相比于低马赫数状态,高马赫数下整个干扰区内气体比热比出现了明显变化,而这将导致气体膨胀特性发生变化,从而导致燃气喷流与空气喷流的差异增大。③ 当攻角为正时,异质气体喷射效应对气动力热特性的影响较小,当攻角为负时,其影响逐渐增大,尤其是对气动热环境的影响。究其原因：正攻角时,喷口位于背风面,喷流干扰强度减小,由于燃气和空气喷流保持了压比、动量比、推力和流量一致,因此二者对气动力热特性的影响差别不大；负攻角时,喷口位于迎风面而激波干扰结构增强,由于燃气、空气喷流的比热比不同、膨胀特性存在差异,在强压缩条件下,这种差异的影响会更加明显。从图28可见,正攻角时,燃气、空气喷流的分离区范围和压力分布差别较小,负攻角时,燃气喷流喷口上游的分离区明显减小,压力峰值位置更接近喷口,表明燃气喷流干扰结构受压缩更明显、激波干扰作用更强,因此干扰区热流峰值更高。