液液针栓多喷注单元喷雾场数值模拟

2023-05-05雷凡培杨岸龙周立新

火箭推进 2023年2期

王凯,唐亮,雷凡培,杨岸龙,周立新

(1.西安航天动力研究所液体火箭发动机技术重点实验室,陕西西安 710100;2.中国船舶工业集团有限公司,北京 100044)

0 引言

自1970年首次被提出以来,针栓式喷注技术就成为变推力技术研究的热点,并被成功应用于月球软着陆探测和火星软着陆探测。SpaceX的Tom Mueller首次将针栓式喷注器技术应用于猎鹰9号火箭的Merlin系列发动机上。不到十年间发展了3个改进版本,使得猎鹰9成为商业航天市场上最具竞争力的火箭[1]。正是由于梅林发动机喷注器结构简单,零组件可快速更换实现改进,固有燃烧稳定性好,变推力实现容易,为猎鹰火箭的可重复使用提供了坚定的动力技术支持。

以往的喷注器推力变比最大为2.5∶1,变比再增大时由于低工况喷注压降降低,会出现低频燃烧不稳定和燃烧效率显著降低的情况。与以往的直流撞击式及离心式喷注器不同,针栓式喷注器通过调节装置改变氧化剂和燃料的喷注面积,保持各个工况的喷注压降基本不变,进而保持高的燃烧效率和良好的燃烧稳定性,因而该喷注器技术成为深度变推力发动机的关键技术。针栓式喷注器具有独特的几何特性和喷注特性,它通过伸入推力室内部的针栓喷注器结构,使轴向的液膜与径向的液束(或液膜)垂直撞击形成喷雾扇雾化混合,其工作原理结构如图1所示[2]。

图1 针栓式喷注器原理图Fig.1 Schematic diagram of pintle injector

目前针栓式喷注器虽已成功应用于工程,然而相关的研究多见于工程研制,基础理论和数值研究很少。近年来,文献[3-6]针对轴对称针栓喷注器开展了很低速度或很低动量比下的喷雾水试试验研究,然而仅能获得描述喷雾场宏观结构的雾化角信息,难以获得喷雾场内部的详细喷雾场结构,雾化角不能够反映喷雾场内部的详细信息。这是因为针栓式喷注器的喷雾场非常独特,液雾极为稠密,且内部呈现复杂的分区结构特征。这给用现有光学测量设备进行试验观测带来了极大挑战,能获得的雾场特性参数非常有限。Sakaki等率先提出了平面针栓喷注器单元的构想,并通过试验清楚观测到了喷雾场和燃烧火焰结构[7]。后续Sakaki等又开展了同参数的轴对称喷注器对比研究,结果表明两者的特性一致,这为使用平面针栓式喷注器代替轴对称针栓式喷注器开展更深层次的研究提供了强有力的依据,表明平面针栓的研究方法非常有效[8]。Kim等借助平面针栓喷注单元的思路,针对径向圆孔的平面针栓多喷注单元开展了自燃推进剂的点火特性研究,结果表明阻塞率0.65的结构在中等动量比下的点火延迟时间最短[9]。成鹏借助同样的思路主要针对径向圆孔的单喷注单元开展了较多的雾化特性试验研究及气氧酒精的平面针栓燃烧过程光学观测研究[10]。王凯等也借助平面针栓喷注单元的思路开展了单喷注单元雾化过程高精度的数值模拟,首次分析了单喷注单元膜束撞击变形特征、流场涡结构特征及雾化混合分区结构特征等,为揭示针栓式喷注器喷雾燃烧精细过程奠定了基础[11]。

然而以上研究未关注实际针栓式喷注器中相邻喷注单元间的相互影响。实际针栓式喷注器相邻喷注单元之间的间距较小,存在较强的相互影响,与单喷注单元形成的喷雾场结构不同,这使得单喷注单元得到的结果不能直接应用于实际针栓式喷注器。因此,需要对多喷注单元与单喷注单元喷雾场之间的差异性进行深入的研究。

综上所述,本文基于AMR技术和分相识别的PLIC VOF新方法,以液液平面针栓多喷注单元为研究对象,通过对径向和轴向两路液体分别追踪,开展液液针栓多喷注单元雾化过程的高精度数值模拟。首先通过数值仿真结果获得了多喷注单元喷雾场特有的结构特征;然后深入分析了相邻喷注单元喷雾扇的撞击变过程,揭示了单元间的相互影响机制;最后对比分析了多喷注单元与单喷注单元喷雾场结构及特性的差异。这对于将单喷注单元研究得到的结果应用于实际针栓式喷注器具有重要的意义。

1 数值仿真方法

1.1 控制方程

树状结构的自然分级特性非常适合多尺度网格的应用,而且四叉树/八叉树网格的灵活性与简单性在处理界面流动时有着巨大的优势。基于这种考虑,Popinet提供了一种处理不可压、变密度、带有表面张力的N-S方程的数值方法,结合四叉树/八叉树离散、投影法与多级泊松求解器[12]。通过使用面向对象的C语言编程实现后即为Gerris软件。Gerris所采用的数值方法可表示为[13]

(1)

本文主要为了研究撞击变形及流动过程,采用分三相VOF方法对两种液体推进剂和环境气体的相界面分别进行追踪。为了分别识别两路液体的变形运动过程,认为两种液体互不相溶,且两者之间存在液相界面。于是分别定义液相1的体积分数c1(x,t)和液相2的体积分数c2(x,t),对应得到的流体密度和动力黏性系数为

工作倦怠发生的原因有很多种，然而我认为员工的分配公平感是现阶段最值得关注的问题之一。分配公平感，指的是职工自身对于资源的分配是否合理的个人性判断，是以个人的公平观为基础的。不公平感是比较出来的，所以其具有相对性，但并无绝对标准。公平感是具有主观性的，完全因个人的价值取向，这样就可能造成因人而异的后果。不公平感是具有扩散性的，职工在有不公平感时，往往会牵扯到整个营销团队，负面情绪波及，这样就会对公司的业绩造成严重的影响，从而影响公司的正常发展。

(2)

式中:ρl1,ρl2和μl1,μl2分别是液相1和液相2的密度和黏度;ρg和μg分别为气相的密度和黏度[16-18]。

Gerris中各项采用相应的格式进行离散。时间采用经典的分裂投影法,达到二阶精度。空间采用四叉树/八叉树进行离散,达到二阶精度,这使得自适应加密算法实现更简易灵活,在不损失计算精度的情况下显著降低了计算量,非常适合处理存在多尺度流动的瞬态流场计算问题[18]。PLIC VOF几何重构方法非常适合用于包含破碎、聚合现象的多相流计算过程,图2为雾化过程仿真自适应生成的网格。Francois在采用连续表面张力模型(continuum surface force,CSF)计算液相表面张力的基础上,通过将表面张力转化为相界面附近的连续体积力源项,结合高度函数曲率估计实现了表面张力的精确求解,并对压力进行修正,确保表面张力和压力梯度之间的平衡[19]。采用单调集成大涡模拟MILES(又称隐式大涡模拟ILES)[20]近似模拟亚格子SGS的能量传递。

图2 基于PLIC VOF方法的AMR技术Fig.2 AMR based on PLIC VOF method

图3 PLIC VOF方法Fig.3 PLIC VOF method

1.2 气液相界面捕捉及表面张力

Gerris使用PLIC VOF方法对界面进行重构。PLIC VOF方法定义函数c作为第一相体积分数来描述相界面(包含相界面的网格为0

Gerris采用的PLIC VOF方法对气液相界面捕捉分两步进行:①界面重构,即如图3所示;②几何通量计算和界面对流的处理[13]。界面法向n可以通过考虑相邻单元的体积分数,在八叉树空间离散产生的3×3×3模板中使用Aulisa等提出的Mixed-Youngs-Centred(MYC)格式进行计算[24]。一旦实现界面重构,在规则笛卡尔网格中采用交替方向方法很容易计算得到几何通量,体积分数对流项采用文献[25]实现的方向交替格式计算,进而对式(2)中体积分数输运方程进行离散计算。这种对流格式可以捕捉到锐利的界面,并且已在实际应用中证明达到二阶精度。

1.3 计算模型

对于平面针栓多喷注单元,建立如图4所示计算域模型,至少选取3个喷注单元作为计算对象,径向孔间距根据阻塞率计算确定。其中h为轴向液膜厚度,u1为轴向速度,d为径向液束直径,u2为径向速度。计算域由10～46个L×L×L基本结构Box构成,根据计算需求选用不同的Box数量构成不同大小的计算域。计算域左部和底部分别为轴向与径向的速度入口边界。壁面均为无滑移壁面条件,两侧设置为对称面边界,其余为出口,采用Outflow边界条件,背压为大气环境。设置第一相和第二相均为水,对应两路液体,第三相为空气。计算域Box的长度为L=10 mm,最高网格等级为9,对应的最小网格约为19.5 μm,最低等级为6,同时绝大部分区域网格不超过100 μm。将两种液相的体积分数及其梯度设置为网格自适应函数,即网格会自动随着体积分数及其梯度大小而加密/粗化。具体结构参数如表1所示,工况参数为u1=20 m/s,u2=25 m/s。

图4 计算域模型Fig.4 Model of simulation zone

表1 径向圆形孔的结构参数

1.4 数据处理方法

图5 测点设置示意图Fig.5 Schematic diagram of measuring point arranged

2 仿真验证与结果分析

2.1 多喷注单元喷雾场结构特征分析

为了降低多喷注单元数值仿真的计算量,同时又能说明多喷注单元特有的喷雾场结构,本节一开始先选用两个径向圆孔的喷注单元开展分析,重点关注喷注单元间相互影响造成的特征差异,后续将针对特征差异形成的原因开展分析。

从图6所示的多喷注单元的喷雾场数值仿真结果可以看到,与单喷注单元喷雾场仿真结果(如图7所示)相比,由于相邻喷注单元之间存在相互影响,多喷注单元的喷雾场主要存在以下几个特有结构:相邻两雾扇相撞背部呈脊状结构,该结构使得雾化区域大于雾化角;两个雾扇相撞会在中间截面汇聚成一道薄液膜,使得液滴在整个雾化角范围内均有分布;相邻两孔之间会有一定宽度的液膜下漏,下漏液膜宽度和下漏率是重要的参数;液滴空间分布及混合比分布也发生变化。

图6 多喷注单元喷雾场结构特征Fig.6 Spray field structures of multi-injector elements

图7 单喷注单元喷雾场结构Fig.7 Spray field structures of injector element

图8 受相邻喷注单元影响的喷雾扇变形过程Fig.8 Deformation process of spray fan affected by adjacent injector elements

图9 相邻两个喷注单元的喷雾扇结构Fig.9 Spray fan structure for adjacent two injector elements

图10 多喷注单元的喷雾扇结构Fig.10 Spray fan structure for multi-pintle injector

图11和图12分别为表1中多喷注单元工况参数下的数值仿真和高速摄影试验拍摄的结果,可以看到三者外形非常相似。

笔者在之前的研究结果中还对试验和数值仿真测得的雾化角及单喷注单元的液滴平均粒径进行了定量对比验证[26],相同工况下雾化角的试验测量值和仿真值分别为72.88°和70.28°,相对误差为3.57%;单喷注单元液滴平均粒径的试验测量值和仿真值分别为175 μm和158 μm,相对误差为9.71%。上述结果表明该仿真方法的准确性较高,也表明上述分析的变形过程及揭示的相邻喷注单元间相互影响的机制是准确合理的。

图11 多喷注单元的喷雾扇数值仿真结果Fig.11 Numerical result of spray fan for multi-pintle injector

图12 多喷注单元的喷雾扇试验拍摄结果Fig.12 Experimental result of spray fan for multi-pintle injector

另外,从多喷注单元的结果可以看出,相邻喷雾扇由于受到相互影响被挤压,最终形成的喷雾场呈扁平的多凹腔状,每个喷注单元将整个喷雾场内部分隔成一个个相对独立的凹腔空间,而外围拼接成一个大雾锥,与单喷注单元的喷雾场不同;同时,液膜破碎后的液雾在空间也呈现相应的不规则曲面分布,从而使得喷雾场的液滴空间分布及混合比分布也发生相应变化。

2.2 单喷注单元与多喷注单元的对比

针对表1相同结构参数(d=0.8 mm,h=0.45 mm)和工作参数(u1=20 m/s,u2=25 m/s)的单喷注单元和多喷注单元,通过数值仿真获得了喷雾场液滴粒径信息。获得的单喷注单元与多喷注单元的喷雾场计算结果对比如图13所示,对应的全场液滴粒径概率密度分布如图14所示,液滴粒径沿径向的空间分布如图15所示,相对流强(即流强概率密度分布)及混合比沿径向的空间分布分别如图16和图17所示。通过对比可以看到,正如2.1节所述,单喷注单元与多喷注单元的喷雾场存在显著不同。

图13 喷雾场计算结果Fig.13 Simulation results of spray field

图14 全场液滴粒径概率密度分布Fig.14 Droplet diameter probability density distribution in the whole spray field

从图14可以看出,与单喷注单元相比,多喷注单元液膜路和液束路的液滴分布峰值均偏向更大粒径方向,大液滴数目占比明显增加,液滴平均粒径S均明显增大,分别从119.38 μm和127.44 μm增大至170.29 μm和171.73 μm,分别增大了42.6%和34.8%。

图15 液滴粒径沿径向的空间分布Fig.15 Radial distribution of droplet diameter

如图16所示,实际的多喷注单元喷雾扇仅有较小半径处的部分发生重叠,液膜厚度加倍。在较小半径位置采用虚线液膜厚度加倍的理论结果,在大半径位置采用实线单喷注单元原来液膜厚度的理论结果,中间半径位置介于两种结果之间。由于中间过渡区雾场结构十分复杂,目前还没更好的理论模型进行解析描述,还不能形成一个统一的解析表达式。目前只能用实线和虚线组合起来描述多喷注单元。这样的模型虽存在不足,但相对简单,可以满足一定精度要求下的理论预估使用需求。