为简化计算,计算域只提取喷嘴以及燃烧室内的流域,不考虑身部冷却通道的流动。由于喷注器喷嘴分布具有周期性,因此选取整体的1/6模型进行计算,计算区域自喷嘴出口至喷管出口扩张段。由于喷管扩张段对燃烧室内的燃烧流场影响不大,因此扩张段只截取到与燃烧室直径相同处。计算域示意图如图2所示。为便于分析燃烧流场结构,定义了几处参考面(线),分别为对称面(A-A)、内圈喷嘴轴线(b)以及对称面上沿轴线不同位置(c1、c2、c3),c1、c2、c3位置分别为x=10d_O2、x=20d_O2和x=30d_O2,其中d_O2=6 mm为氧喷嘴出口直径。

图2 计算域示意图
Fig.2 Schematic diagram of computational domain

燃烧室内的流动为气液两相湍流燃烧流动过程,不考虑壁面与外界的换热,忽略热辐射、重力等因素的影响。气液掺混燃烧过程的控制方程采用湍流气粒两相反应流动的Navier-Stokes(N-S)方程组^[15]。

湍流模型采用SST k-ω模型,该模型在氢氧火箭发动机推力室的燃烧仿真中得到较好的应用^[16-17]。为分析氢氧燃烧流场的细节结构,燃烧模型采用涡耗散概念(EDC)模型,该模型可以将细致的化学反应机理应用到湍流反应流动中,其中的化学反应机理采用6组分/9步化学反应(6s9r)机理,组分分别为H₂、O₂、H₂O、H、O、OH,化学反应常数

参考文献[18]中的数值。

本文研究的室压范围跨越了液氧的临界压力,对于液氧喷嘴的喷雾过程,存在亚临界和超临界喷射两种状态。然而有研究表明,混合物的临界压力可能超过组分的临界压力,超临界喷射过程仍可能表现出亚临界的破裂特征^[19-20]。因此,本文忽略液氧液膜的破碎过程,利用文献[21]中同轴直流式喷嘴的经验公式,采用离散相模型直接在喷嘴出口截面定义氧液滴的速度及粒径分布,液滴尺寸服从Rosin-Rammler分布。液滴的轨迹采用Lagrange坐标下的随机轨道模型模拟。液氧的蒸发速率较大,液滴蒸发采用对流/扩散控制(convection/diffusion-controlled)模型计算,该模型可以很好地模拟高蒸发速率下液滴表面到气相的斯蒂芬流动。

对于物性参数,气相的黏度及导热系数均采用FLUENT软件中的分子动力学理论(kinetic-theory)计算各组分参数,再通过理想气体混合定律(ideal-gas-mixing-law)计算混合物的物性参数。液相物性参数的计算模型

参考文献[22],其中定压比热的计算公式为

(C_p)/R=(C⁰_p)/R+(C^r_p)/R(1)

C⁰_p/R=a₀+a₁T+a₂T²+a₃T³+a₄T⁴(2)

式中:C_p为定压比热; R为通用气体常数; C⁰_p为理想气体的热容,按照式(2)计算,其中系数a₀～a₄由文献[22]得到; C^r_p为剩余热容,由文献[22]得到。

蒸发潜热采用Watson方程计算,即

L₂=L₁(1-T_r2)^0.38/(1-T_r1)^0.38(3)

式中:L₂为所求温度T₂下的未知蒸发潜热; L₁为温度T₁下的已知蒸发焓; T_r2为T₂温度下的对比温度; T_r1为T₁温度下的对比温度。液相物性参数通过用户自定义函数UDF写入FLUENT软件。

对计算域进行六面体网格划分,在喷嘴内部、氢氧剪切层、边界层附近区域进行了加密处理,网格示意图如图3所示。

图3 网格示意图
Fig.3 Grid schematic diagram

采用FLUENT商业软件对模型进行了求解计算。采用有限体积法对控制方程进行离散和求解,压力和速度耦合关系采用SIMPLEC算法计算,压力离散格式采用二阶格式,其余项均采用二阶迎风格式。

气氢入口采用质量流量入口边界。喷管出口为超声速流动,出口参数采用上游参数插值得到。壁面采用绝热无滑移边界条件,对于离散相采用反射边界条件。

为验证网格无关性,首先采用了3套网格进行了计算。分别在喷嘴内部、氢氧剪切层以及燃烧室区域进行了加密,网格单元数分别为205 786、744 872、1 668 856,分别记为工况1、工况2、工况3。为节省计算资源,采用氢氧总包反应对燃烧流场进行了计算,并提取了沿内圈喷嘴轴线上的O₂质量分数分布和无量纲温度分布,如图4所示。由图4可知,工况2和工况3的计算结果基本一致,采用网格数较少的工况2计算能达到网格无关性的要求。

图4 网格无关性验证结果
Fig.4 Verification results of grid independence

仿真得到的部分工况对称面上温度和O₂质量分数分布云图如图5所示。由图5可以看出,不同室压工况宏观上的燃烧流场结构基本相似。

不同室压工况下的颗粒轨迹如图6所示。由图6可以看出,室压5.1 MPa以上工况颗粒轨迹长度基本相同,与室压5.0 MPa以下工况存在明显差异。液氧的临界压力为5.043 MPa^[23],当室压大于液氧临界压力时,氧液滴的蒸发潜热为0。通过仿真可知,当室压超过液氧临界压力时,氧液滴自喷嘴喷出后,吸热升温至临界温度并完全蒸发的轨迹长度基本相同; 而室压低于临界压力时,氧液滴的蒸发长度随室压的降低而增大。

图5 不同室压工况对称面上温度和O2质量分数分布
Fig.5 Temperature and O2 mass fractional distribution on the symmetric surface under different combustor pressure

图6 不同室压工况下颗粒轨迹
Fig.6 Particle trajectories under different chamber pressure

为定量分析燃烧流场结构,分别提取了参考线b上的OH质量分数分布以及温度分布,分别如图7和图8所示。由图可知,室压分别处于液氧亚临界和超临界压力范围时,OH质量分数以及温度在距喷注面75～175 mm范围内,分别表现出两种不同的分布形式。这说明液氧在亚临界和超临界压力下的燃烧流场结构存在差异。而同处于超临界压力下的不同工况,OH质量分数以及温度分布基本相同,说明燃烧流场结构是基本相似的。

图7 沿内圈喷嘴轴线OH质量分数分布
Fig.7 Mass fraction distribution of OH along the axis of inner injector

图8 沿内圈喷嘴轴线无量纲温度分布
Fig.8 Distribution of dimensionless temperature along the axis of inner injector

对称面不同位置c1、c2和c3处的O₂质量分数分布如图9所示。由图9可知,在10d_O2(60 mm)处,不同工况下的O₂质量分数分布基本相同; 在20d_O2(120 mm)处,超临界压力下各工况的O₂质量分数明显大于亚临界压力各工况; 在30d_O2(180 mm)处,不同工况下的O₂质量分数分布又基本相同。在75～175 mm范围内,处于液氧超临界压力下各工况的氧液滴蒸发量明显大于亚临界压力下的各工况,这也造成了两者之间燃烧流场结构的差异。

图9 对称面上不同位置O2质量分数分布
Fig.9 Mass fraction distribution of O2 at different positions on symmetrical plane

由仿真结果可知,在推进剂介质、温度、混合比以及喷注速度相同的前提下,不同室压下影响气氢/液氧燃烧流场动力结构的主要因素为液氧的蒸发特性。当燃烧室压力高于液氧的临界压力时,交界面温度高于液氧的临界温度,液氧从喷嘴进入燃烧室后,液滴温度经历跨临界变化,这时液氧液滴表面张力和蒸发潜热减小,气液两相间的差别消失,导致混合扩散过程控制燃烧过程^[24-25]。当燃烧室压力由2.8 MPa提高到9.8 MPa时,液氧射流状态经历了跨临界变化,燃烧过程由雾化蒸发控制转变为由混合扩散过程控制,因此二者燃烧流场形态存在差异。而当室压由5.1 MPa提高到9.8 MPa时,燃烧过程由混合扩散过程控制的状态未发生变化,因此燃烧流场形态相似。当室压继续提高时,燃烧过程仍由混合扩散过程控制,由此推断,燃烧室内的流场结构变化不会太大,仍将保持相似。因此,对于处于超临界压力下的高压推力室的燃烧性能,在保证推进剂介质、混合比、温度以及喷注速度相同的前提下,可以用低压缩尺试验来评估实际工作压力下推力室的燃烧特性,但低压缩尺试验的室压不应低于液氧的临界压力。