连续液相碎化过程数值模拟研究

＊兰治科 李勇 苏光辉 昝元锋

（1.中国核动力研究设计院 四川 610041 2.西安交通大学核科学与技术学院 陕西 710049）

液体雾化过程包含了能量和质量转化、连续液体颗粒化等物理过程，在现代社会被广泛应用各个领域，如农业喷雾、喷涂、发动机、航空航天、核能等。根据用途不同，对雾化性能的要求也不尽相同。因此获得特定喷雾头的相关性能参数显得尤为重要，特别是雾化后的液滴粒径、流量密度分布和雾化锥角。

由于喷雾机理的复杂性，到目前为止还没有能够完整、准确描述喷雾头性能的好方法。目前的研究方法主要有三种，分别为实验研究、理论计算和计算流体力学（CFD）模拟。其中理论计算主要是建立在半经验公式的基础上，而关于实验研究，高昂的实验成本以及测量精度，仍是主要的制约因素。特别是应用于某些特殊领域的大型喷雾头，如稳压器、安全壳中的喷雾设备，其共同特点是：喷雾头尺寸和喷雾流量大，液滴粒径分布范围较广，普通粒径测量仪器不适用于该类型喷雾头测量。此外，对于喷雾头内部流场以及出口处液膜破碎过程等微观现象研究，目前尚未找到适合的实验方法。

随着计算机性能提高，CFD被广泛应用于雾化性能各项参数研究。采用数值模拟方法研究雾化现象有助于更加深入认识雾化机理和指导相关实验，同时可以弥补喷雾实验不足，提供喷雾实验无法得到的流场分析结果。

然而，对于液体雾化数值模拟，CFD技术仍然面临着较大挑战。主要是因为雾化过程不仅牵涉到从单相液体到气液共存两相流动的过渡，同时，对于液相来说，还存在从连续相的液体到离散相液滴的转变。目前几乎所有的雾化场模拟均使用Eulerain-Lagrangian两相流模型[1]。这种模型通常假定喷口附近的液滴粒径和速度分布作为已知的边界条件已经被确定，而获得这些边界条件的方法又往往依赖实验或者经验。因此，从严格意义上说，这种数值方法并不能实现喷雾过程的真正模拟。而Eulerain-Eulerain方法，相对于网格尺寸而言，仅局限于对两种连续相共存的流场进行计算，并不能获得完全离散化的液滴[2]。

基于以上原因，本文根据最大不稳定波增长率液膜破碎理论，在连续液相出口位置构建初级破碎模型-PBM，与喷雾头内部流场和雾化场的数值模拟进行耦合，提出了一种全流场的数值模拟方法，实现从喷雾头入口到雾化场的全过程数值模拟。最后，通过已有的实验数据对以上方法进行了验证，计算结果符合较好，证明了全流场数值模拟方法的有效性。

1.全流场模型构建

目前，在已有的两相流模拟中，Eulerain-Eulerain法将两相视为能够互相穿插，但不相溶的连续介质。该方法使用空间和时间上为连续函数的相体积分数来模拟两相在空间的分布。从喷雾头入口全液相流动到出口连续液膜与气相并存，符合Eulerain-Eulerain方法两相模拟的要求。

Eulerain-Lagrangian方法用于喷雾场模拟，将喷雾空间的气体介质作为连续介质，为其建立Eulerain坐标系；而将液滴颗粒当作离散介质，为其建立Lagrangian坐标系。该方法假定液相以离散液滴的形式进入连续气相中，通过质量交换、能量交换和动量交换进行两个相间的耦合。因此Eulerain-Lagrangian方法仅用于喷雾头出口连续液膜破碎后的两相流场模拟。

本文使用Eulerain-Eulerain法中VOF模型捕捉两相之间界面，模拟从单相液体进入喷雾头，直到出口位置形成液膜的内部流场；使用Eulerain-Lagrangian方法中DPM模型模拟离散液滴在气相中的碰撞、破碎以及再分布过程。本文采用最大表面扰动波增长率破碎理论，构建了适合于大型喷雾头的初级破碎模型（PBM）。通过在FLUENT求解器上使用UDF，编译相关程序与以上两个模型耦合，形成VOF+PBM+DPM的完整流场的数值模拟方法，模型结构见图1所示。

图1 连续液相碎化过程数值模拟模型结构图

2.全流场数值模拟过程及结果分析

(1)喷雾头内部流场数值模拟

VOF方法是在1981年由Hirt和Nichols[3]首先提出，它使用一个标量场函数γ体积份数来表征第二相流体在计算网格中的分布，其大小为第二相流体在网格中所占体积百分比。因此，VOF又被称为流体体积法。在本文雾化模拟中，环境介质空气作为第一相流体，雾化介质水作为第二相流体，所以体积分数γ代表网格控制体中雾化介质水的体积份额。

根据体积分数的标量场γ，全计算域流体物性密度ρ以及粘性μ的表达式如下：

其中ρ1、ρg分别为液相和气相所对应的密度；μ1、μg分别为液相和气相所对应的粘性应力。因此，在液相、气相以及两相混合区域均使用式（2）中流体物性，可以使在流体力学相关控制方程中按单流体模型仅使用一套方程组。

体积分数函数γ满足如下输运方程：

喷雾头内部结构和流道对雾化特性起决定作用，因此在计算过程中需要对网格质量进行独立性验证。本文定义了4种不同的网格方案，对网格逐渐加密对流量特性及雾化锥角的影响进行了比较。当网格精度进一步提高，以上特性参数不再明显变化时即认为网格质量满足独立性条件。

本文为螺旋槽式旋流喷头提供连续液相入口，在大多数情况下均表现为强烈的旋流特征，为了避免标准κ-ε模型的缺陷，湍流模型考虑κ-ε模型的修正方案：RNGκ-ε模型。FLUENT求解器三维旋流计算中通过旋流因子as对RNG κ-ε模型中的μt进行修正。计算前对as=0.070，0.075和0.080三种不同旋流因子时流动特性计算值与实验值进行了比较，认为as=0.080作为湍流模型旋流因子比较适合该类型喷雾头的旋流强度。

连续液相破碎数值模拟，使用VOF模型从喷雾头入口为时间和迭代开始起始点。图2为喷雾头内部流场瞬态图。喷雾头在流动扰动过程中，在喷雾头混合腔接近出口位置出现较为明显的旋流、直流分层流动。其中直流部分位于正中心竖直向下，而旋流则沿四周壁面绕流。而在旋流作用下，混合腔内部易形成负压区，容易发生大气压力下环境介质中空气倒流现象。而直流与旋流交界处则成为相对较为“脆弱”的位置。

图2 连续液相碎化为液带过程图

(2)喷雾头出口初级破碎模拟

在喷雾头出口建立液膜射流的几何模型，假定液膜以锥形射进入环境空气。液膜与空气相互作用，在液膜表面上形成扰动波。

通常认为，扰动波满足如下关系式：

假设扰动波是导致液膜破碎[4]的主要原因，当扰动波增长率达到最大值时，液膜发生破裂，形成初始液滴。通过色散方程对式（4）进行求解，最终得到最大增长率对应下液膜破碎破碎模型PBM。

使用UDF中的DEFINE_EXECUTE_AT_END宏，在瞬态计算的每个时间步长统计喷雾头孔口所在面每个网格单元的流动信息。通过采集两相VOF模型在喷雾出口的计算结果作为边界条件输入雾化初级破碎模型。对旋流喷雾头，假定孔口液膜以环状形式存在，通过液相截面积相等的方法获得PBM初始液膜厚度。

通过PBM模型获得破碎距离后，在FLUENT求解器中通过UDF功能统计某一时间段内在该距离横截面上瞬态流量分布，获得该截面上雾化场流量密度分布。UDF由编辑宏DEFINE_EXECUTE_AT_END程序以及定义F_UDMI存储瞬态数据完成，在横截面每个face类型的网格单元中累积该时间步下经过该网格的质量流量，最终得到某一时间段内截面上的流量分布如图3所示。

图3 液膜破碎截面流量密度分布

(3)喷雾场数值模拟

在DPM模型中，通过积分拉式坐标下的颗粒作用力微分方程来求解离散相颗粒的轨道。颗粒运动方程：

其中u、up分别为流体和颗粒的速度；ρ、ρp分别为流体和颗粒的密度；FD(u-up)为颗粒的单位质量曳力。此外，液滴碰撞模型基于O'Rourke理论；液滴破碎模型适用TAB模型。

为了降低计算成本，本文构建2-D几何模型对喷雾场竖直截面内离散液滴流动特性进行数值模拟。该几何模型在原有VOF模型的正下方，选取过喷雾头轴心的竖直截面作为计算区域，其大小为4000mm×2200mm。

通过DPM模型可得到喷雾场温度、速度等各项参数分布。图4为雾化场液滴粒径分布计算结果。从图中可以看出，喷雾场中较大粒径液滴主要位于雾化场外侧。原因在于喷雾头雾化过程以旋流特征为主，即大部分雾化液体以旋流方式沿雾化锥角外侧方向进入喷雾场。其次，雾化场粒径分布图中左半侧液滴粒径均值明显大于右半侧，且颗粒分布较密。主要是因为在离散相边界条件上，喷雾头雾化液体并非对称分布，大部分液体分布在边界条件的左侧区域。

图4 喷雾场液滴粒径分布

3.全流场数值模拟与实验结果对比

针对本文研究对象，在数值模拟工况下开展了相应实验研究，具体实验研究内容和方法详见文献[5]，此处不再赘述。

图5给出了通过全流场数值模拟方法计算得到的喷雾场液滴粒径概率密度与实验数据的比较。图中横坐标为液滴粒径尺寸，曲线为计算得到粒径分布，离散点为实验测试得到的实际粒径概率分布，两者分布趋势符合较好，表明了全流场数值模拟的有效性。

图5 液滴粒径概率密度计算值与实验数据比较

4.结论

本文针对连续液相，提出一种全流场液相破碎数值模拟方法。通过构建初级破碎模型PBM，与用于内部流场计算的VOF模型和雾化场计算的DPM模型进行耦合，实现了从连续液相到雾化颗粒的喷雾全过程数值模拟。最后与实验结果对比，证明了VOF+PBM+DPM全流场模型的有效性。解决了目前雾化数值模拟中，Eulerain-Eulerain和Eulerain-Lagrangian方法仅能模拟局部过程的问题，为连续液相破碎及雾化场研究提供了一套完善的数值计算方法。