NUMERICAL SIMULATION OF ATOMIZATION PROCESS OF SELF-IMPINGEMENT AND IMPINGEMENT BIPROPELLANTS
-
摘要: 肼类双组元推进剂在卫星和导弹等各类航天器动力系统液体火箭发动机中应用广泛. 为了深入理解自击撞击式肼类双组元推进剂的雾化过程, 基于OpenFOAM开源平台, 通过引入有效液体相分数, 发展了三维三相VOF-LPT转换算法, 并进一步开发了该转换算法的耦合求解器; 结合网格局部加密技术与自适应加密技术, 实现了对液体界面的动态捕捉; 通过与前人横向燃料喷射实验结果的对比, 验证了所开发耦合求解器的准确性和有效性. 利用所开发的耦合求解器数值研究了入射速度扰动、燃料偏靠角度和射流速度对自击撞击式肼类双组元推进剂雾化效果的影响. 结果表明, 入射速度扰动使液滴尺寸分布整体减小, 比表面积增加, 有助于提升雾化效果; 对于燃料液膜和氧化剂液膜的撞击段, 增大燃料偏靠角度或提高射流速度均会增加燃料液膜与氧化剂液膜之间的有效撞击动量, 从而提高了燃料液膜与氧化剂液膜之间的混合度; 在喷雾场完全发展后, 液滴尺寸分布随燃料偏靠角度的增大而整体增大, 随射流速度的增加而整体减小. 此外, 由于燃料偏靠角度限制了混合区域的范围, 导致液膜混合破碎程度随射流速度的增加逐渐趋于定值.Abstract: Hydrazine-based bipropellants are widely employed in liquid rocket engines of various aerospace propulsion systems for the space crafts such as satellites and missiles. To deepen the understanding of the atomization process of self-impingement and impingement hydrazine-based bipropellants, a three-dimensional three-phase volume of fluid (VOF) to Lagrangian particle tracking (LPT) conversion algorithm is developed based on the open-source platform OpenFOAM by introducing the effective liquid phase fraction. A coupled solver is further established for the VOF-LPT conversion algorithm. By combining the local grid refinement technology with the adaptive mesh refinement technology, the liquid interface is dynamically captured. The accuracy and validity of the developed coupled solver are well verified by comparing with the previous experimental results on a fuel jet in a cross flow. With the developed coupled solver, the effects of injecting velocity disturbance, fuel biased angle and jet velocity on the atomization progress of self-impingement and impingement hydrazine-based bipropellants are numerically investigated. The results indicate that the introduction of a moderate injecting velocity disturbance reduces the overall droplet size distribution and increases the specific surface area, thereby beneficial to improve the atomization efficiency. During the impingement stage of the fuel film and the oxidant film, the increase of the fuel biased angle and the jet velocity both promote the mixing between the fuel film and the oxidizer film due to the enhanced effective impingement momentum between the films. After the spray field is fully developed, the overall droplet size distribution increases with the increase of the fuel biased angle, and decreases with the increase of the jet velocity. In addition, since the mixing region of the fuel film and the oxidizer film is restricted by the fuel biased angle, the degree of mixing and breakup of the films is approaching a certain level with the increase of the jet velocity.
-
Keywords:
- bipropellant /
- VOF-LPT /
- droplet /
- spray
-
引 言
双组元液体推进剂以其高比冲、可控推力和多次点火等优点, 在航天器的飞行姿态和轨道控制中发挥着关键作用. 工程实践表明, 自击撞击式双组元推进系统中喷雾过程对发动机的燃烧和效率具有显著影响. 该过程主要包括一次破碎和二次破碎两个阶段. 一次破碎是指射流从喷嘴中喷出, 气液界面在空气扰动的作用下形成振动波, 振动波逐渐发展最终导致射流破裂成液片. 在二次破碎过程中, 这些液片进一步分解为由较小液滴组成的液雾[1], 这一过程一直持续到表面张力占主导地位, 形成稳定的球形液滴[2].
喷雾模拟一直是计算流体力学的研究热点. 典型的多相流模拟方法, 如水平集(level set, LS)方法[3]、流体体积(volume of fluid, VOF)方法[4]、耦合水平集和流体体积(coupled level set and volume of fluid, CLSVOF)方法[5-6]等界面捕捉方法, 通常用于模拟一次破碎. 只要液体结构, 如液体射流、液片或单个液滴, 能够被计算网格充分解析, 就可以使用这些方法精确描述雾化过程. 然而, 采用界面捕捉方法解析离散的小液滴云将导致巨大的计算成本. 与此相反, 拉格朗日粒子追踪(Lagrangian particle tracking, LPT)方法是专门为跟踪大量单个液滴或颗粒而开发的. 该方法要求计算网格尺寸必须明显大于颗粒的特征直径, 因此计算成本非常低. 结合液滴破裂[7-8]和碰撞模型[9-10], LPT方法特别适合于模拟二次破裂及随后的小液滴云的发展过程.
由于界面捕捉方法和拉格朗日方法各自涵盖了雾化过程的不同阶段, 因此需要将两者结合起来才能建立完整的喷雾模拟模型. 然而, 这两种方法对网格分辨率和控制方程求解的要求存在很大不同, 导致在耦合时需要引入特殊处理. 一种处理方法是使用二维耦合平面, 并在该平面上以映射的方法将上游欧拉框架(Euler system)中的液滴转换为下游的拉格朗日粒子, 如Grosshans等[11]和Spitzenberger等[12]的方法. 然而, 这种方法需要了解射流破碎的特征并提前确定耦合平面的位置. 另一种处理方法是采用三维耦合算法, 其优点是可以在指定的体积区域内跟踪和转换液滴, 因此无需事先了解射流破裂的情况. Herrmann[13]和刘昌波等[14]分别采用水平集方法和VOF方法处理大结构液片问题, 同时对较小的液片及液滴进行三维跟踪; 对于满足一定条件的液滴, 储存其位置、速度和尺寸信息, 并将其从欧拉框架中移除, 注入到拉格朗日框架下. 由于仅包含从欧拉框架到拉格朗日框架的转换过程, 通常称该方法为单向耦合算法. 除了从欧拉到拉格朗日的转换过程, 双向耦合算法也包括从拉格朗日粒子转换回欧拉框架的过程. Edelbauer等[15]在喷雾模拟中使用了双向耦合算法. 然而, Ling等[16]的测试结果表明, 单向耦合和双向耦合算法对喷雾模拟结果的影响差别很小.
以上提到的欧拉-拉格朗日三维耦合方法主要是针对单一组元的喷雾. 然而, 近来对双组元喷雾系统尤其是双组元自击撞击式喷雾系统的关注日益增加. 撞击式喷雾一方面可以加速液滴与环境流体之间的相互作用, 减少液滴冲击壁面现象[17]; 另一方面, 在双组元液体推进系统中, 撞击式喷雾还能实现不同的反应性和浓度分层, 有利于控制雾化和燃烧过程[18].
本文在OpenFOAM框架下发展了三维三相VOF-LPT耦合算法, 模拟研究了双组元自击撞击式喷雾过程, 揭示了入口速度扰动、自击偏靠角度和射流速度对液滴粒径分布与空间分布的影响机制, 以期为双组元推进剂系统的优化设计提供理论参考.
1. 耦合算法与实施
三相VOF-LPT耦合算法涉及三个部分, 即三相VOF方法、LPT方法以及三相VOF-LPT转换算法, 下面对它们一一进行介绍.
1.1 三相VOF方法
采用三相VOF方法在欧拉框架下捕捉燃料、氧化剂和空气3种流体的界面. 相分数${\alpha _1}$, ${\alpha _2}$和$ {\alpha _3} $分别代表燃料、氧化剂和空气在网格单元中的体积分数. 此外, 燃料、氧化剂和空气的密度和运动黏度分别用${\rho _1}$, ${\rho _2}$和${\rho _3}$及${\nu _1}$, ${\nu _2}$和${\nu _3}$来表示. 在VOF方法中整个计算域的流体被视为单一的连续体, 流体的密度及运动黏度为
$$\qquad\qquad\qquad \rho = \sum\limits_{i = 1}^3 {{\alpha _i}{\rho _i}} $$ (1) $$\qquad\qquad\qquad \nu = \frac{1}{\rho }\sum\limits_{i = 1}^3 {{\alpha _i}{\rho _i}{\nu _i}} $$ (2) 此外, 假设流体为不可压缩牛顿流体. 因此, 流体流动可由下面的连续性方程和动量守恒方程来控制
$$ \nabla \cdot {\boldsymbol{u}} = 0 $$ (3) $$ \frac{{\partial \rho {\boldsymbol{u}}}}{{\partial t}} + \nabla \cdot (\rho {\boldsymbol{uu}}) = - \nabla p + \rho{\boldsymbol{ g}} + \nabla \cdot {\boldsymbol{\tau}} + {{\boldsymbol{f}}_{st}} + {{\boldsymbol{f}}_s} $$ (4) 式中, $ {\boldsymbol{u}} $是流体速度, ${\boldsymbol{g}}$是重力加速度, ${\boldsymbol{\tau}} $是黏性应力张量, $ {{\boldsymbol{f}}_{st}} $是表面张力项. 最后一项$ {{\boldsymbol{f}}_s} $表示流场中拉格朗日粒子对流体的反作用力, 计算式为${{\boldsymbol{f}}_s} = - \dfrac{1}{{{V_i}}}\displaystyle\sum\limits_j {\left( {{m_{pj}}\frac{{{\mathrm{d}}{{\boldsymbol{u}}_{pj}}}}{{{\mathrm{d}}t}}} \right)} $, 其中${V_i}$为目标单元网格的体积, ${{\boldsymbol{u}}_{pj}}$和${m_{pj}}$分别为目标单元网格中第j个拉格朗日粒子的速度和质量. 此外, 引入3个相分数输运方程用以描述燃料、氧化剂和空气在欧拉框架下的空间分布以及相界面
$$ \left. \begin{split} & {\frac{{\partial {\alpha _1}}}{{\partial t}} + \nabla \cdot \left( {{\alpha _1}{\boldsymbol{u}}} \right) = 0} \\ & {\frac{{\partial {\alpha _2}}}{{\partial t}} + \nabla \cdot \left( {{\alpha _2}{\boldsymbol{u}}} \right) = 0} \\ & {\frac{{\partial {\alpha _3}}}{{\partial t}} + \nabla \cdot \left( {{\alpha _3}{\boldsymbol{u}}} \right) = 0} \end{split} \right\} $$ (5) 且这3种流体的相分数满足守恒性条件
$$ \sum\limits_{i = 1}^3 {{\alpha _i} = 1} $$ (6) 表面张力使用连续表面力(continuum surface force, CSF)模型[19-20]来建模
$$ {{\boldsymbol{f}}_{st}} = \sum\limits_i {\sum\limits_{j \ne i} {{\sigma _{ij}}{\kappa _{ij}}{{\boldsymbol{\delta}} _{ij}}} } $$ (7) 其中, ${\sigma _{ij}}$是相$i$和相$j$之间的表面张力系数; $ {\kappa _{ij}} $是界面$ij$的曲率, 可以通过相分数来计算, 即[20]
$$ {\kappa _{ij}} = - \nabla \cdot \frac{{{\alpha _j}\nabla {\alpha _i} - {\alpha _i}\nabla {\alpha _j}}}{{\left| {{\alpha _j}\nabla {\alpha _i} - {\alpha _i}\nabla {\alpha _j} + sf} \right|}} $$ (8) 式中, $sf$是考虑网格不均匀性的稳定因子, 定义式为$sf = \varepsilon {\left[ {\displaystyle\sum\limits_{i = 1}^{{n_{cv}}} {{V_i}{{({n_{cv}})}^{ - 1}}} } \right]^{1/3}}$, ${n_{cv}}$为计算单元的总数, ${V_i}$是单元$i$的体积, $\varepsilon $是一个很小的正数, 取值为$\varepsilon = 1.0 \times {10^{ - 8}}$. 最后一项${{\boldsymbol{\delta}} _{ij}}$定义为: ${{\boldsymbol{\delta}} _{ij}} = {\alpha _j}\nabla {\alpha _i} - {\alpha _i}\nabla {\alpha _j}$, 表示仅仅在相界面上施加表面张力.
1.2 LPT方法
LPT方法用于追踪拉格朗日包裹(Lagrangian parcel)在整个计算域中的运动. 每个包裹代表一定数量的液滴, 这些液滴具有相同的直径和速度. 由于液滴尺寸足够小, 拉格朗日包裹被视为质量点, 忽略其体积的影响. 根据牛顿第二定律和轨迹方程, 在每个时间步更新包裹的速度和位置
$$\qquad\qquad {m_p}\frac{{{\mathrm{d}}{{\boldsymbol{u}}_p}}}{{{\mathrm{d}}t}} = {{\boldsymbol{F}}_D} + {{\boldsymbol{F}}_G} $$ (9) $$\qquad\qquad \frac{{{\mathrm{d}}{{\boldsymbol{x}}_p}}}{{{\mathrm{d}}t}} = {{\boldsymbol{u}}_p} $$ (10) 其中, ${{\boldsymbol{u}}_p}$为包裹速度, ${m_p}$为质量, ${{\boldsymbol{x}}_p}$为包裹位置矢量; 对于本文等温系统, 包裹上的作用力仅考虑曳力和重力. 曳力的计算公式如下[21]
$$ {{\boldsymbol{F}}_D} = {C_D}\frac{{\text{π} D_p^2}}{8}{\rho _3}\left( {{{\boldsymbol{u}}_c} - {{\boldsymbol{u}}_p}} \right)\left| {{{\boldsymbol{u}}_c} - {{\boldsymbol{u}}_p}} \right| $$ (11) 其中, ${D_p}$代表包裹直径, ${\rho _3}$代表气相密度, ${{\boldsymbol{u}}_c}$代表包裹位置${{\boldsymbol{x}}_p}$处的速度(取自欧拉流场). 阻力系数${C_D}$根据Schiller-Naumann曳力模型计算[21]
$$ {C_D} = \left\{\begin{split} & {\frac{{24}}{{R{e_p}}}\left( {1 + 0.15Re_p^{0.687}} \right)},\quad {{\text{if }}\;R{e_p} \leqslant 1000} \\ & {0.44},\quad {{\text{if }}\;R{e_p} > 1000} \end{split} \right. $$ (12) 其中$R{e_p} = \left| {{{\boldsymbol{u}}_p}} \right|{D_p}/{\nu _p}$, ${\nu _p}$为包裹运动黏度. 重力和浮力作为一个合力计算
$$ {{\boldsymbol{F}}_G} = {m_p}{\boldsymbol{g}}\left( {1 - {\rho _3}/{\rho _p}} \right) $$ (13) 流体动力导致的二次破裂采用Reitz-Diwaker模型[22]进行模化处理, 包裹与包裹之间的碰撞以及随后的合并采用Nordin的碰撞模型[23]进行模化计算.
1.3 VOF-LPT转换算法
为了模拟整个雾化过程, 需要实现小液滴从VOF到LPT的转换, 对应的算法包括3个步骤: (1) 识别液体相分数场${\alpha _l}$(${\alpha _l} = {\alpha _1} + {\alpha _2}$)中的单个液滴; (2) 计算液滴属性; (3) 如果液滴符合转换要求, 则将液滴注入拉格朗日框架(Lagrangian system)进行转换, 并在欧拉框架下删除该液滴. 详细的实施过程如下.
如图1所示, 第1步是在全计算域中使用三维四连通空间连通分量标记(connected component labeling, CCL)算法[24]在全计算域中识别液体相分数场中的独立液滴. CCL算法具有能够较好地处理不同形状和尺寸的液滴及易于实现等特点, 广泛应用于液体射流雾化破碎模拟研究[8,24-25]. 其步骤如下: 首先定义一阈值${\alpha _t}$, 如果网格的液体相分数${\alpha _l}$大于${\alpha _t}$, 则被定义为液体, 标记为1; 否则, 它们将被定义为气体, 标记为0, 并在后续步骤中不予考虑. 不失一般性, 本文取${\alpha _t} = 0.01$, 它能较好地保证耦合过程中气液质量守恒. 依以上操作循环处理所有网格后, 为每个液体网格分配唯一的编号. 最终, 通过将每个液滴所包括的最小的编号赋值给液滴所包含的其他所有网格, 完成对液滴的识别和标记.
第2步, 计算每个液滴的体积、位置、速度与有效直径. 对于给定液滴, 计算方法如下
$$\qquad\qquad\qquad\quad {V_p} = \sum\limits_i {{\alpha _{li}}{V_i}} $$ (14) $$\qquad\qquad\qquad\quad {{\boldsymbol{x}}_p} = \frac{1}{{{V_p}}}{\alpha _{li}}{V_i}{{\boldsymbol{x}}_i} $$ (15) $$ \qquad\qquad\qquad\quad {{\boldsymbol{u}}_p} = \frac{1}{{{V_p}}}{\alpha _{li}}{V_i}{{\boldsymbol{u}}_i} $$ (16) $$\qquad\qquad\qquad\quad {d_p} = \sqrt[3]{{\frac{{6{V_p}}}{\text{π} }}} $$ (17) 其中, 索引i遍历属于液滴的所有网格. 在计算液滴体积${V_p}$、位置${{\boldsymbol{x}}_p}$和速度${{\boldsymbol{u}}_p}$时考虑了液体相分数${\alpha _l}$, 确保位于界面上的网格($0 < {\alpha _l} < 1$)仅计算其有效液相质量.
第3步, 对每个液滴进行检查, 看是否满足转换要求. 由于流动复杂且高度紊乱, 采用易于定义和计算的液滴几何特征为转换要求, 主要要求有两个[24]:
(1) 液滴的直径必须小于设定的某个临界直径${d_{\min }}$(${d_{\min }}$取值为网格特征长度的5倍), 即${d_p} < {d_{\min }}$;
(2) 液滴的球度${\phi _p}$必须小于2: ${\phi _p} < 2$.
球度的计算方法如下
$$ {\phi _p} = \frac{{2{r_{\max }}}}{{{d_p}}} $$ (18) 式中, ${r_{\max }}$为液滴质量中心到液滴表面的最大距离. ${\phi _p} = 1$表示液滴是一个完美的球体, 而${\phi _p} > 1$表示液滴为非球形.
如果满足以上两个转换要求, VOF液滴就会转换成拉格朗日包裹. 同时, 相应的空间位置上的速度以及液滴相分数被设置为0, 以保证流场的质量和动量守恒. 值得注意的是: 大变形液滴的运动无法用LPT模型准确预测, 它仍会保留在欧拉框架中, 不会进行转换.
图 1 在液体相分数场$ {\alpha }_{l} $中识别、标记液滴的过程. 不失一般性, 液体相分数阈值取为$ {\alpha }_{t} = 0.01 .$ 左侧图中蓝色线表示气液相界面Figure 1. Algorithm for identifying and labeling the droplets in the phase fraction field $ {\alpha }_{l} .$ Without loss of generality, the threshold value of the liquid phase fraction is taken as $ {\alpha }_{t} = 0.01 .$ The interface between gas and liquid is represented by the blue line in the left graph2. 模型验证
由于缺乏双组元自击撞击式喷雾的实验数据, 使用Gopala等[26]横流燃料喷射的实验结果来验证我们发展的三维三相VOF-LPT耦合求解器.
图2展示了与实验设置一致的计算域示意图. 该计算域由尺寸为80D × 80D × 160D的矩形空气通道组成, 其中蓝色部分为射流入口, D代表喷嘴直径, 大小为8 mm. 假设左边入口处气流充分发展, 喷嘴位于气体入口下游, 燃料入口处的速度是均匀的. 定义燃料和空气之间的动量比为动量通量比, 其值为10. 表1给出了相关的流动和物性参数.
表 1 流动和物性参数Table 1. Flow parameters and fluid propertiesProperty Fuel Air velocity/(m·s−1) 32 116.6 density/(kg·m−3) 780.6 5.88 kinematic viscosity/(m2·s−1) 8.96 × 10−7 3.15 × 10−6 surface tension/(N·m−1) 0.07 计算网格设置如下: 整个计算域采用正六面体背景网格, 网格大小为$\Delta x = D/10$. 在射流出口区域以及射流发展的主要区域, 进行网格局部加密, 加密等级为3级, 以确保在欧拉框架下捕捉喷雾发展的最小网格尺寸为$\Delta x = D/80$. 同时采用自适应加密技术, 在局部加密区域之外, 使用$\Delta x = D/80$的网格对液滴进行动态捕捉.
如图2所示, 在喷嘴出口下游60D处设置一个后处理平面, 所有的VOF液滴在到达该处之前已全部转换为拉格朗日粒子, 统计穿过该平面的每个拉格朗日粒子的信息, 并用于计算索特平均直径(Sauter mean diameter)${D_{32}}$, 其中
$$ {D_{32}} = \frac{{\displaystyle\int_{{D_{\min }}}^{{D_{\max }}} {{d^3}{\mathrm{d}}N} }}{{\displaystyle\int_{{D_{\min }}}^{{D_{\max }}} {{d^2}{\mathrm{d}}N} }} $$ (19) 式中, $N$为直径为$d$的液滴数目. 图3(a)给出了后处理平面上y方向不同位置处本文模拟获得的索特平均直径与实验结果的对比. 可以看出模拟所得到的索特平均直径完全落在实验结果的标准偏差范围($ \pm 5$μm)之内. 此外, 提取模拟所得到的燃料入口下游不同位置处的射流穿透距离${x_p}$, 并与实验结果进行了比较. 如图3(b)所示, 模拟与实验结果吻合良好. 以上结果表明本文所提出三维三相VOF-LPT耦合求解器模拟喷雾过程的准确性.
3. 双组元推进剂雾化过程研究
3.1 问题描述
图4展示了一种自击撞击式喷注单元示意图, 这里燃料和氧化剂先自击形成液膜, 后液膜之间相互撞击. 燃料喷注口直径为${D_1}$, 氧化剂喷注口直径为${D_2} = 1.4{D_1}$. 图中给出了自击结构, 其中燃料自击角${\alpha _1}$ = ${60^ \circ }$, 撞击高度${H_1} = 4.3{D_1}$, 氧化剂自击角${\alpha _2}$ = ${40^ \circ }$, 撞击高度${H_2} = 8.6{D_1}$. A-A剖面展示了燃料和氧化剂的偏靠结构, 燃料偏靠角度为${\beta _1}$, 氧化剂偏靠角度为${\beta _2}$, 撞击距离$L = 6.45{D_1}$为喷注盘面上燃料/燃料(氧化剂/氧化剂)喷注中心轴线的距离.
图5给出了双组元推进剂自击撞击雾化模拟所使用的网格, 计算域取为40 mm × 24 mm × 20 mm的长方体区域. 网格设置类似于验证算例: 背景网格大小为$\Delta x = {D_1}/10$, 网格局部加密应用于圆柱射流所在的区域以及液膜发展的主要区域, 同时在局部加密区域之外采用网格自适应加密技术, 二者加密等级均为3级, 即最小网格尺寸为$\Delta x = {D_1}/80$. 以偏二甲肼为燃料、四氧化二氮为氧化剂, 模拟所需的流速和物性参数(如密度、黏度系数以及表面张力系数等)详见表2.
表 2 流动和物性参数Table 2. Flow parameters and fluid propertiesProperty Fuel Oxidizer Air velocity/(m·s−1) 12.96 9 density/(kg·m−3) 791 1458 3.69 kinematic viscosity/(m2·s−1) 8.8 × 10−7 3.1 × 10−7 2.7 × 10−6 surface tension/(N·m−1) 0.033 0.027 3.2 扰动对雾化过程的影响
在实际应用中, 液体射流喷口可能受到各种形式的扰动影响, 使燃烧室中液态燃料的雾化破碎过程呈现出一定的非定常特性[27-28]. 因此, 本节设计了两个对比算例, 其中一算例有扰动, 一算例无扰动. 对于扰动算例, 扰动以正弦波的形式加在入射速度上, 即$u = {u_0}\left[ {1 + u' _n\sin (2\text{π} ft)} \right]$. 这里的$u$表示喷口实时速度, ${u_0}$表示未施加扰动时的入射速度, $u'_n$为无量纲扰动振幅, 取为0.05; $f$为扰动频率, 取为100 kHz. 不失一般性, 我们选择燃料/氧化剂偏靠角度${\beta _1} = {15^ \circ }$, ${\beta _2} = {10^ \circ }$.
图6展示了有、无扰动工况下的喷雾场. 其中, 红色代表燃料, 青色代表氧化剂. 从图6中可以观察到, 有、无扰动工况下喷雾场的分布基本是一致的. 此外, 液体比表面积的定义为: 指定区域内所有液体的面积(m2)与液体的体积(m3)之比, 其可在一定程度上代表液膜的破碎程度. 针对如图6红框区域所示的雾化特征区域, 我们定量分析了扰动对液滴比表面积的影响. 该特征区域x, y和z方向长度分别取为0.010, 0.016和0.010 m, 区域体心位于燃料、氧化剂喷注口之间中心位置下游0.017 m处. 统计发现有扰动工况下液体比表面积(13325 m−1)比无扰动工况下液体比表面积(13032 m−1)略高, 表明有扰动工况下液膜的破碎程度更剧烈, 有助于燃料和氧化剂之间的混合.
该喷雾整体过程如下, 由于撞击高度较小且速度较大, 燃料射流率先完成了自击, 并形成了液膜, 在燃料液膜发展一段时间之后, 与氧化剂射流撞击在一起, 两者相互挤压混合. 燃料液膜的头部, 由于速度较快, 没有受到氧化剂的影响, 仍然向下游运动了一段距离之后才发生破碎. 除此之外, 燃料与氧化剂完全混在一起. 当喷雾场达到稳定时, 统计整个计算域中液滴的粒径分布. 由于在喷雾场达到稳定时刻VOF液滴数目少且尺寸较大, 其对液滴尺寸分布的影响极小, 因此这里为简单起见仅统计拉格朗日液滴尺寸, 其中拉格朗日液滴尺寸通过式(17)计算. 如图7所示, 在有、无扰动两种工况下, 液滴尺寸的概率分布趋势相同. 大部分液滴的尺寸都位于50 ~ 160 μm之间, 有、无扰动工况占比分别为75%和76%. 在30 ~ 90 μm区间, 有扰动工况下在该区间的液滴数目占比为37.6%, 高于无扰动工况下在该区间的液滴数目占比(34.8%); 而在120 ~ 210 μm区间, 则是无扰动工况下液滴数目的占比更高, 为31.9%, 高于有扰动工况下液滴数目的占比27.4%. 这表明, 在自击撞击式双组元推进系统中, 适度的扰动有助于雾化效果的提升.
3.3 燃料偏靠角度对雾化过程的影响
在自击撞击雾化过程中, 影响雾化效果的关键因素是燃料和氧化剂的自击过程, 以及所形成的液膜之间相互撞击的过程. 现有的研究已广泛探讨了自击过程, 并明确了撞击角度对雾化效果的影响[29-30]. 具体来说, 随着自击角度的增加, 液膜的破碎长度减少, 喷雾场中液滴的平均粒径也减小. 当撞击角度过大时, 燃料和氧化剂可能会反向流向喷嘴所处壁面. 因此, 本节专注于探讨燃料偏靠角度对雾化效果的影响规律. 在3.2节有扰动算例的基础上, 设计了两个新算例. 这两个新算例中燃料偏靠角度分别设置为${\beta _1} = {20^\circ}$, ${10^\circ}$, 而其他参数与3.2节有扰动的算例完全相同.
图8展示了不同燃料偏靠角度下雾化过程的模拟结果. 从全局视角看, 喷雾场的发展轨迹和结构在各个偏靠角度下都相似. 但在红色框标注的区域, 明显可见随着偏靠角度的增加, 在低偏靠角度下并未完全破碎的氧化剂液膜逐渐破碎. 定量对比红框所示特征区域内液体比表面积发现, 燃料偏靠角度为10°, 15°和20°时, 液体比表面积分别为13026, 13325和13976 m−1. 以上结果均表明液膜破碎程度随燃料偏靠角度的增加而增强.
在两个液膜发展的平面上, 将液膜速度分解为x和y方向的分量. 实际上, 参与撞击的是y方向的速度分量. 随着燃料偏靠角度的增加, y方向的速度分量也相应增大, 从而增加了有效的撞击动量. 这导致整个雾化过程的不稳定性增加, 同时燃料和氧化剂的混合程度也随之增加. 从图9中液滴粒径的统计结果可见, 燃料偏靠角度的增加导致液滴粒径整体上升. 这是因为在小偏靠角度条件下更多的液膜是在发展中受到空气扰动作用而破碎, 而在大偏靠角度条件下更多的液膜是受液膜之间撞击作用而破碎. 液膜在发展过程中是逐渐变薄的, 在燃料液膜和氧化剂液膜相互撞击时, 二者液膜厚度仍较厚, 因此导致液膜之间撞击破碎更容易形成大尺寸的液滴.
3.4 射流速度对雾化过程的影响
射流速度在撞击式喷雾系统中扮演着至关重要的角色[31-32]. 在3.2节中有扰动的算例的基础上, 设计了两个新算例, 射流速度${U_0}$提高至原来的1.5倍与2倍, 同时保持燃料和氧化剂混合比与动量比均不变.
图10展示了不同射流速度下的喷雾场结果. 计算红框特征区域内的液体比表面积发现, 速度为${U_0}$, $1.5{U_0}$和$2{U_0}$时的液体比表面积分别为13325, 13497和13519. 速度为$1.5{U_0}$和$2{U_0}$时液体比表面积较接近, 均明显大于速度为${U_0}$时的液体比表面积. 该结果表明, 类似于3.3节的分析, 随着射流速度的增加, 有效的撞击动量也随之增大, 使得整个混合和雾化过程更加剧烈; 然而, 由于燃料偏靠角度限制了混合区域的范围, 导致混合破碎程度随着射流速度的增加而逐渐趋于一个阈值.
在液滴整体粒径随射流速度的变化方面, 如图11所示, 在30 ~ 90 μm的范围内, 随着射流速度的增加, 液滴尺寸的概率密度变大, 说明小尺寸液滴占比更大. 而在120 ~ 210 μm范围内, 液滴尺寸的概率密度随射流速度的增加而减小, 说明大尺寸液滴更少. 综合而言, 液滴尺寸随射流速度的增加整体减小. 然而, 在更大的液滴尺寸上, 3个工况并未显示出明显的差异.
4. 结 论
本文在OpenFOAM框架下成功开发了一种适用于双组元推进系统喷雾过程模拟的三相VOF-LPT耦合求解器, 并通过横流燃料喷射算例进行了测试, 证实了该算法能够定量模拟包含一次破碎到二次破碎的全雾化过程. 然后, 利用该耦合算法对肼类双组元推进系统的自击撞击雾化过程进行了高精度数值模拟, 系统研究了入口速度扰动、自击偏靠角度和射流速度对雾化效果的影响, 发现在肼类双组元推进系统的雾化过程中, 扰动并未显示出明显的影响, 但适度的扰动对提高雾化效果具有积极作用. 随着燃料偏靠角度的增大, 燃料与氧化剂之间的混合程度增加. 然而, 较大的偏靠角度导致生成更多尺寸较大的液滴. 在保持动量比和混合比不变的前提下, 提高射流速度将增强燃料和氧化剂的混合程度. 然而, 当射流速度达到一定阈值时, 混合程度不再随射流速度的增加而增加. 此外, 随着射流速度的增加, 生成液滴中小液滴的占比也随之增加. 本文研究不仅提供了适用于喷雾全过程精确模拟的三相流VOF-LPT耦合算法, 而且可为肼类双组元推进系统的设计和性能优化提供理论参考.
-
图 1 在液体相分数场$ {\alpha }_{l} $中识别、标记液滴的过程. 不失一般性, 液体相分数阈值取为$ {\alpha }_{t} = 0.01 .$ 左侧图中蓝色线表示气液相界面
Figure 1. Algorithm for identifying and labeling the droplets in the phase fraction field $ {\alpha }_{l} .$ Without loss of generality, the threshold value of the liquid phase fraction is taken as $ {\alpha }_{t} = 0.01 .$ The interface between gas and liquid is represented by the blue line in the left graph
表 1 流动和物性参数
Table 1 Flow parameters and fluid properties
Property Fuel Air velocity/(m·s−1) 32 116.6 density/(kg·m−3) 780.6 5.88 kinematic viscosity/(m2·s−1) 8.96 × 10−7 3.15 × 10−6 surface tension/(N·m−1) 0.07 表 2 流动和物性参数
Table 2 Flow parameters and fluid properties
Property Fuel Oxidizer Air velocity/(m·s−1) 12.96 9 density/(kg·m−3) 791 1458 3.69 kinematic viscosity/(m2·s−1) 8.8 × 10−7 3.1 × 10−7 2.7 × 10−6 surface tension/(N·m−1) 0.033 0.027 -
[1] 曹建明. 液体喷雾学. 北京: 北京大学出版社, 2013 (Cao Jianming. Liquid Sprays. Beijing: Peking University Press, 2013 (in Chinese) Cao Jianming. Liquid Sprays. Beijing: Peking University Press, 2013 (in Chinese)
[2] Gorokhovski M, Herrmann M. Modeling primary atomization. Annual Review of Fluid Mechanics, 2008, 40(1): 343-366 doi: 10.1146/annurev.fluid.40.111406.102200
[3] Poblador-Ibanez J, Nocivelli L, Magnotti GM, et al. A physics-driven Σ-Y atomization model for heavy-duty engine simulations. International Journal of Multiphase Flow, 2023, 167: 104523 doi: 10.1016/j.ijmultiphaseflow.2023.104523
[4] Li HF, Feng JH, Cao XY, et al. Simulation study of the swirl spray atomization of a bipropellant thruster under low temperature conditions. Energies, 2022, 15(23): 8852
[5] Chang JL, He LJ, Chen LH, et al. Atomization of liquid pulsed jet in subsonic crossflow. AIP Advances, 2023, 13(5): 055117
[6] Lyras P, Hubert A, Lyras KG. A conservative level set method for liquid-gas flows with application in liquid jet atomisation. Experimental and Computational Multiphase Flow, 2023, 5(1): 67-83 doi: 10.1007/s42757-021-0119-1
[7] Berni F, Sparacino S, Riccardi M, et al. A zonal secondary break-up model for 3D-CFD simulations of GDI sprays. Fuel, 2022, 309: 122064 doi: 10.1016/j.fuel.2021.122064
[8] Bhatia B, Johny T, De A. Understanding the liquid jet break-up in various regimes at elevated pressure using a compressible VOF-LPT coupled framework. International Journal of Multiphase Flow, 2023, 159: 104303 doi: 10.1016/j.ijmultiphaseflow.2022.104303
[9] Feng SQ, Zhang SL, Shi JD. Using radius of influence and adaptive collision mesh to capture droplet collision of a diesel engine. Fuel, 2023, 352: 129071 doi: 10.1016/j.fuel.2023.129071
[10] Zhang QK, Zhang P, Chi YC, et al. Eulerian-Lagrangian simulation and validating experiment of n-butanol/biodiesel dual-fuel impinging sprays. Fuel, 2023, 350: 128761 doi: 10.1016/j.fuel.2023.128761
[11] Grosshans H, Szász RZ, Fuchs L. Development of an efficient statistical volumes of fluid-Lagrangian particle tracking coupling method. International Journal for Numerical Methods in Fluids, 2014, 74(12): 898-918
[12] Spitzenberger A, Neumann S, Heinrich M, et al. Particle detection in VOF-simulations with OpenFOAM. SoftwareX, 2020, 11: 100382 doi: 10.1016/j.softx.2019.100382
[13] Herrmann M. On simulating primary atomization using the refined level set grid method. Atomization and Sprays, 2011, 21(4): 283-301 doi: 10.1615/AtomizSpr.2011002760
[14] 刘昌波, 雷凡培, 周立新. 两股湍流射流撞击雾化过程的数值研究. 推进技术, 2014, 35(12): 1669-1678 (Liu Changbo, Lei Fanpei, Zhou Lixin. Primary atomization simulation of two turbulent impinging jets. Journal of Propulsion Technology, 2014, 35(12): 1669-1678 (in Chinese) Liu Changbo, Lei Fanpei, Zhou Lixin. Primary atomization simulation of two turbulent impinging jets. Journal of Propulsion Technology, 2014, 35(12): 1669-1678 (in Chinese)
[15] Edelbauer W, Birkhold F, Rankel T, et al. Simulation of the liquid break-up at an AdBlue injector with the volume-of-fluid method followed by off-line coupled Lagrangian particle tracking. Computers & Fluids, 2017, 157: 294-311
[16] Ling Y, Zaleski S, Scardovelli R. Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model. International Journal of Multiphase Flow, 2015, 76: 122-143 doi: 10.1016/j.ijmultiphaseflow.2015.07.002
[17] Akagawa H, Miyamoto T, Harada A, et al. Approaches to solve problems of the premixed lean diesel combustion. SAE Transactions, 1999, 108: 120-132
[18] Lim JH, Reitz RD. High load (21 bar IMEP) dual fuel RCCI combustion using dual direct injection. Journal of Engineering for Gas Turbines and Power, 2014, 136(10): 101514 doi: 10.1115/1.4027361
[19] Brackbill JU, Kothe DB, Zemach C. A continuum method for modeling surface tension. Journal of Computational Physics, 1992, 100(2): 335-354 doi: 10.1016/0021-9991(92)90240-Y
[20] Dritselis C, Karapetsas G. Open-source finite volume solvers for multiphase (n-phase) flows involving either Newtonian or non-Newtonian complex fluids. Computers & Fluids, 2022, 245: 105590
[21] Schiller L. A drag coefficient correlation. Zeit Ver Deutsch Ing, 1933, 77: 318-320
[22] Reitz R. Modeling atomization processes in high-pressure vaporizing sprays. Atomisation and Spray Technology, 1987, 3(4): 309-337
[23] Nordin PAN. Complex Chemistry Modeling of Diesel Spray Combustion. Gothenburg, Sweden: Chalmers University of Technology Sweden, 2001
[24] Heinrich M, Schwarze R. 3D-coupling of volume-of-fluid and Lagrangian particle tracking for spray atomization simulation in OpenFOAM. SoftwareX, 2020, 11: 100483 doi: 10.1016/j.softx.2020.100483
[25] Zhou WY, Chen B, Zhu QB, et al. Numerical simulation of angled-injected liquid jet breakup in supersonic crossflow by a hybrid VOF-LPT method. International Journal of Multiphase Flow, 2023, 166: 104503 doi: 10.1016/j.ijmultiphaseflow.2023.104503
[26] Gopala YS, Zhang P, Bibik O, et al. Liquid fuel jet in crossflow-trajectory correlations based on the column breakup point//48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Orlando, Florida, 2010
[27] 周程林, 邹建锋, 张阳等. 喷口速度扰动下液体射流的频率响应. 航空动力学报, 2022, 37(4): 743-754 (Zhou Chenglin, Zou Jianfeng, Zhang Yang, et al. Frequency response of liquid jet under nozzle velocity disturbance. Journal of Aerospace Power, 2022, 37(4): 743-754 (in Chinese) Zhou Chenglin, Zou Jianfeng, Zhang Yang, et al. Frequency response of liquid jet under nozzle velocity disturbance. Journal of Aerospace Power, 2022, 37(4): 743-754 (in Chinese)
[28] 李珍妮. 纵向扰动控制下液体射流破碎机理的研究. [硕士论文]. 天津: 天津大学, 2016 (Li Zhenni. Investigation on the breakup mechanism of the liquid jet under longitudinal disturbance. [Master Thesis]. Tianjin: Tianjin University, 2016 (in Chinese) Li Zhenni. Investigation on the breakup mechanism of the liquid jet under longitudinal disturbance. [Master Thesis]. Tianjin: Tianjin University, 2016 (in Chinese)
[29] 张蒙正, 张泽平, 李鳌等. 两股互击式喷嘴雾化性能实验研究. 推进技术, 2000, 21(1): 57-59 (Zhang Mengzheng, Zhang Zeping, Li Ao, et al. Experimental research on spray properties of unlike impinging injectors. Journal of Propulsion Technology, 2000, 21(1): 57-59 (in Chinese) doi: 10.3321/j.issn:1001-4055.2000.01.017 Zhang Mengzheng, Zhang Zeping, Li Ao, et al. Experimental research on spray properties of unlike impinging injectors. Journal of Propulsion Technology, 2000, 21(1): 57-59 (in Chinese) doi: 10.3321/j.issn:1001-4055.2000.01.017
[30] 李佳楠, 费俊, 杨伟东等. 直流互击式喷注单元雾化特性准直接数值模拟. 推进技术, 2016, 37(4): 713-725 (Li Jianan, Fei Jun, Yang Weidong, et al. Quasi-Direct numerical simulation on atomization. Journal of Propulsion Technology, 2016, 37(4): 713-725 (in Chinese) Li Jianan, Fei Jun, Yang Weidong, et al. Quasi-Direct numerical simulation on atomization. Journal of Propulsion Technology, 2016, 37(4): 713-725 (in Chinese)
[31] Kuznetsov GV, Strizhak PA, Valiullin TR, et al. Atomization behavior of composite liquid fuels based on typical coal processing wastes. Fuel Processing Technology, 2022, 225: 107037
[32] Coratella C, Sahu A, Parry L, et al. Experimental investigation of the orifice-to-orifice variability in jet velocities in a diesel injector. Fuel, 2020, 280: 118353 doi: 10.1016/j.fuel.2020.118353