某发动机喷管内流动与换热的瞬态模拟

2020-09-27龚红兰

上海理工大学学报 2020年4期

龚红兰，李凌

（上海理工大学能源与动力工程学院，上海 200093）

喷管属于发动机的重要组成部分之一[1]，作为燃烧室燃烧产物的出口，它具有能量转换的作用，通过喷管可以使推进剂燃烧产生的高温高压气体膨胀加速，将燃烧产物的热能转换为气体动能，从而使发动机获得推力。在发动机工作条件一定的情况下，喷管的结构及性能会对发动机推力大小和推进效率高低产生影响。因此，对发动机喷管的特性进行深入研究，揭示其规律，对于喷管的设计及优化十分重要。

鉴于喷管内流场数值计算的重要性，国内外学者进行了大量研究与分析。文献[2-5]对喷管的分离流动进行了数值计算；文献[6-8]研究了喷管的瞬态流动特性；文献[9-11]对发动机喷管进行了分离流动实验和数值模拟研究。但是，大部分文献只模拟分析了喷管内流动特性，而对于发动机喷管内的传热规律及特性却涉及较少。故本文以某发动机喷管为研究对象，借助CFD[12]技术对喷管内部流动及换热情况进行了瞬态数值模拟，获得了不同时刻喷管内流场、压强、温度分布情况。

1 数值模型

1.1 计算模型及网格划分

本文模拟的喷管喉部直径为12 mm，收敛半角为30° ，扩张半角为10° ，收敛段半径为18 mm。考虑到喷管构型的轴对称性，建立二维轴对称模型，如图1 所示。对计算区域进行结构化网格划分，由于壁面存在边界层，对靠近喷管壁面以及喉部附近进行加密，以提高计算结果的精度，且网格均匀过渡，第一个内节点布置在对数律成立的范围内，网格划分情况如图2 所示。首先进行了网格无关性验证，不同网格下喷管喉部中心点的温度和速度大小如图3 所示。可以看出当网格数从7 788 变至15 356 时，随着网格数的增加，喉部中心点的温度和速度的变化很小，满足无关性要求。因此，选取网格数为7 788 的网格作为本文计算网格。

图 1 喷管计算区域简图Fig. 1 Schematic diagram of computational domain of the nozzle

图 2 计算区域结构网格划分Fig. 2 Division of the computational domain with structured grids

1.2 数学模型

利用计算流体力学软件Fluent 对喷管内部进行瞬态模拟，控制方程由守恒型Navier-Stokes 方程、湍流输运方程及状态方程组成。

a. 守恒型方程的通用形式。

式中：V为速度矢量； ϕ为通用变量； Γϕ为广义扩散系数；Sϕ为相应的源相；ρ为燃气的密度；t为时间。

图 3 网格无关性验证结果图Fig.3 Grid independence verification results

b. 状态方程。

式中，p，T分别为燃气的压强和温度。

c. 湍动能方程。

d. 湍流耗散方程。

式中：k为湍流动能； ε为湍流耗散率；xi，xj分别为i，j方向上的坐标；ui为流体时均速度在i方向上的分量； µ为分子黏性系数； µt为湍流黏性系数，取值为0.09； Gk为由于平均速度梯度产生的湍流动能；C1ε，C2ε，σk， σε均为经验常数，其取值分别为1.44，1.92，1.0，1.3。

1.3 初始和边界条件及燃气物性参数

以发动机燃烧室压强达到发动机工作压强的时刻为计算起点，初始压力为7 MPa，初始温度为300 K，计算区域边界条件设置情况如表1 所示。表中：p0为喷管入口总压；T0为喷管入口总温；pb为喷管出口压强；Tb为喷管出口温度；∂T/∂n为壁面法线方向上的燃气温度梯度；u，v分别为燃气速度在轴向和径向上的分量；∂u/∂r，∂p/∂r，∂T/∂r分别为速度、压强、温度在径向上的梯度。

为便于分析，引入假设：a. 喷管内燃气为理想气体；b. 不考虑燃气的化学反应；c. 不考虑辐射传热及内壁面烧蚀。

燃气物性参数：燃气密度按理想气体处理；燃气气体常数R为320 J/（kg·K），燃气比热比γ为1.2216；导热率由kinetic-theory 给定[13]；黏性按Sutherland 定律处理（三系数形式）[14]。

式中：T0为参考温度，T0=273.11 K；S为等效温度，S=110.56 K； µ0为T0时的参考黏性系数，µ0=1.716 ×10−5kg/（m·s）。

表 1 计算区域边界条件设置Tab.1 Boundary condition settings for the computational domain

1.4 数值方法

使用基于密度的耦合隐式求解器求解Navier-Stokes 方程组。耦合隐式方法虽然比分离算法需要更多的内存，但更易收敛，对于高速可压缩流来说，耦合算法更加适合。对流项采用二阶迎风格式进行离散，湍流模型选用RNGk−ε两方程模型，该模型对应变率较大、存在分离现象及流线弯曲程度较大的流动问题比较适合，采用标准壁面函数（standard wall function）来求解近壁区域物理量。

2 计算结果分析

2.1 模型验证

为验证本文模型，采用本文方法模拟了文献[15]的双圆弧喷管，将其计算结果与文献中的结果进行了比较，如图4 所示。从图中可以看出采用本文方法得到的结果与参考文献的结果比较一致，验证了文中数值方法的正确性。

2.2 流场分析

本文模拟分析了在 Pa 的燃烧室压强条件下喷管内流场及温度场随时间的变化情况，图5为喷管内燃气在不同时刻的速度与马赫数沿轴向变化曲线，图6 为喷管内不同时刻的压力沿轴向变化曲线。从图5 可以看出，在发动机燃烧初期，喷管内燃气的流速随时间增大，但仍为亚音速，马赫数小于1，燃气以亚音速流动，其压力沿喷管轴线均呈单调下降趋势。这个过程一直持续到燃气在喉部达到音速之前，如图6 所示。0.5 ms时，喷管喉部处的气流流速达到临界速度，此时马赫数等于1，但在喉部之后速度又降低，燃气压力在喷管收敛段呈下降趋势，但在扩张段压力略微上升并趋向于一个相对稳定的平台压强，如图6所示。0.6 ms 时，燃气在喷管收敛段以亚音速加速流动，在临界位置达到声速，然后在喷管扩张段经过一段距离的超音速加速后速度突然减小，而后燃气流速又增加，如图5（a）所示。由于扩张段中出现了激波，导致通道内燃气突然减速流动，而且由于激波的产生使得气体压力在激波之前一直处于下降状态，遇到激波后压力突升。速度又增加是因为管口膨胀波传出管外，管口压力下降，进而速度上升，如图6 所示。图7 是不同时刻通道内燃气的速度分布图。从图5 和图7可以看出，在0.6～1.3 ms 时间段内，随着通道内燃气的膨胀，其在喷管扩张段超音速加速流动的距离逐渐增加。同时，激波在慢慢向喷管出口方向移动，燃气压力仍呈现出先下降，然后突然升高，再下降至外界大气压的规律，如图6 所示。从1.4 ms 开始，燃气流动速度继续增大，喷管内激波已经消失不见，燃气压力再次呈单调下降趋势。但越是靠近喷管出口，压力变化越平缓，喷管内气体流速和压力按照这种变化规律一直持续到稳定状态为止。2.2 ms 开始，速度的增幅减小，到3.0 ms 时流动达到稳态，此时出口处燃气的流速达到最大，面平均马赫数为3.35。

图 4 模拟结果与参考文献的比较Fig.4 Comparison of simulation results with references

图 5 不同时刻速度及马赫数沿轴向变化曲线Fig.5 Axial variation curve of velocity and Mach number at different time

图 6 不同时刻轴向压力变化曲线Fig. 6 Axial pressure curve at different time

图 7 不同时刻喷管内流场速度云图Fig.7 Flow field velocity map in the nozzle at different time

图 8 不同时刻轴向温度变化曲线Fig.8 Axial temperature curve at different time

2.3 温度场分析

图8 是不同时刻喷管内燃气温度随时间的变化曲线。从图8（a）可以看出：在0.6 ms 之前，喷管收敛段前段温度高于300 K，且随着时间的增加，高于300 K 的部分增多；收敛段后段温度等于初始温度300 K，而扩张段部分温度低于300 K，这是由于高温燃气热量还未传递至扩张段时燃气在扩张段已经开始膨胀造成的，且随着时间的推移，扩张段温度开始低于300 K 的点逐渐向喷管喉部靠近，且扩张段温度逐渐减小；0.6 ms 时，正如前面所说，喷管扩张段中出现了激波，这对喷管内流场温度分布也会产生影响，喷管内燃气温度在激波之前一直处于下降趋势，遇到激波后温度突然上升。同样，从图8（b）可以看出：在0.6～1.3 ms 时间段内，随着通道内气体的膨胀，激波逐渐往喷管出口方向移动；到1.4 ms 时，喷管内已无激波存在，燃气温度呈下降趋势，收敛段前端温度已有了明显的提升，扩张段后段燃气温度仍然低于300 K，如图9（a）所示；直到1.8 ms 时，高温燃气热量传至扩张段中部使得整个扩张段燃气温度都高于300 K，如图9（b）所示。此后，随着燃气的膨胀，喷管内燃气温度不断升高，直到3.0 ms 时达到稳定状态。