杨立军 黄东骐 韩旺 李敬轩 富庆飞

(1. 北京航空航天大学 宇航学院, 北京 100083; 2. 北京航空航天大学宁波创新研究院, 宁波 315800)

推进剂雾化包括液体射流失稳、断裂、破碎、液滴生成等一系列过程。 推进剂的雾化性能会显著影响发动机性能;研究射流雾化机理有助于实现对雾化过程的精准调控,从而优化后续的蒸发燃烧过程,提升发动机的推力性能。 尽管前人对液体射流雾化开展了大量的研究,但对其雾化机理还存在认识不清的地方[1]。

射流雾化的研究方法主要分为2 类:实验测量和数值模拟。 实验测量方法的优势在于容易开展参数可调的大规模研究,但也面临着测量精度有限,测量物理量较少等挑战。 数值模拟方法虽然计算量较大,但能够提供更高的精度,同时也能提供更多实验中不易获得的数据,帮助开展更加细致的研究,因此受到广泛的关注。

在雾化的直接数值模拟(DNS)研究中,通常采用水平集[2]、流体体积[3]等方法来捕捉气液界面。 例如,Leboissetier 和Zaleski[4]最早对柴油燃料射流进行了DNS,该研究表明喷嘴内部的液体射流在层流条件下不会发生破裂雾化。 这一结果与实验对比显示出了较为良好的一致性。 对于喷嘴外部的射流,Klein[5]指出,当平面液体射流的韦伯数小于270 时,没有发现射流破裂及雾化的行为。 Sander 和Weigand[6]发现,除了喷嘴出口处的湍流特性会对射流液面的稳定性造成影响,速度型也会显著影响射流雾化的结果。 Ménard等[7]首次在两相环境下研究了柴油燃料射流DNS的网格分辨率问题;Desjardins 和Pitsch[8]则使用DNS 研究了中等韦伯数下(We=500 ～2 000)湍流射流液面的破裂行为。 Shinjo 和Umemura[9]在层流环境及静止空气中,在大韦伯数下(We=14 000)使用高分辨率网格对圆柱柴油燃料射流进行了数值模拟。

以上数值模拟研究展现了较为清晰的射流雾化过程,发现韦伯数和流场环境是影响射流雾化的重要因素,但是有关射流雾化的详细机制仍然不清。 基于此,本文尝试通过流动拓扑理论分析液体射流雾化机理。 流动拓扑理论由Chong等[10-11]提出。 流动拓扑理论主要是基于速度梯度张量的3 个不变量来定义的不同流场拓扑结构,这种方法被广泛运用于解决流体问题及燃烧问题。 Chong 等[12]运用该方法研究了壁面剪切流动,Elsinga 和Marusic[13]研究了均匀且各向同性的湍流流动。 Hasslberger 等[14]首次在气泡空气和水的两相流流动中研究了流动拓扑的分布,并发现代表旋转运动的拓扑结构主要存在于气泡内部。 Watanabe 等[15]研究了在湍流与非湍流界面附近的拓扑结构。 Josef 等[16]研究了柴油燃料射流雾化过程中的流动拓扑结构影响;Han 等[17]研究了湍流流动拓扑结构对湍流非预混火焰切向扩散的影响,从流动拓扑结构角度给出了经典火焰面模型假设不成立的条件。

基于上述讨论,本文采用水平集方法,对湍流液体平面射流直接进行三维数值模拟。 水平集方法有2 个优点:①相比于传统的欧拉法,水平集方法可以直接在直角网格中对曲线和曲面进行数值计算,而不用将曲线和曲面参数化;②水平集方法易于追踪物体结构的改变。 因此,使用水平集方法进行DNS 可以有效提高计算精度,并减少计算成本。 本文通过使用水平集方法得到高分辨率数值模拟数据,从流动拓扑角度揭示液体平面射流雾化的机理。

1 模型与算法

1.1 控制方程与算例设置

考虑在静止气体环境中的液体平面射流,其液体射流湍流初速度场在一个槽道中产生(见图1,右侧液体初始湍流场由左侧充分发展的槽道流速度场提供)。 射流的质量方程和动量方程由经典的不可压缩Navier-Stokes 方程给出:

图1 湍流液体平面射流算例设置Fig.1 Turbulent liquid plane jet arithmetic setup

式中:u为速度场;ρ为密度;p为压强;μ为黏度。

界面Γ将液相和气相分隔开。 值得注意的是,在气液多相流中,物理量在两相界面处存在突变。 令液相密度ρ=ρl,气相密度ρ=ρg,液相动力黏度μ=μl,气相动力黏度μ=μg。 因此,在界面Γ处满足:[ρ]Γ=ρl-ρg,[μ]Γ=μl-μg。 而由于速度连续,在界面处[u]Γ=0。 影响气液界面传输计算精度的主要是如下的压力跳跃:

[P]Γ=σκ+ 2[μ]ΓnT·(Δu)·n(3)

式中:σ为表面张力系数;κ为界面曲率;n为界面法向量。

本文使用低马赫多相湍流燃烧求解器NGA[18]对式(1) ～式(3)方程进行数值求解。 在进行数值模拟时,为解决界面处的压力跳跃和密度跳跃,本文使用了精确水平集守恒方法ACLS,在处理表面张力时应用了虚拟流动方法GFM[19]。

图1 展示了流动计算域。 射流初始高度为H,入口处平均速度为U0,整个计算域在流向和展向上(即x,z方向)是周期性边界条件。 在流动高度方向(即y方向)是出口边界条件。 整个计算域在x-y-z方向上尺寸是4H×6H×4H,对应每个方向的网格数目是256 ×384 ×256,网格大小为H/64[8]。 而提供湍流初速度流场的槽道流计算域如图1 所示,其计算尺寸是4H×H×4H,对应每个方向的网格数目是256 ×64 ×256,槽道流的高度为H,与后续射流初始高度对应。 射流初始高度处于-H/2 ＜y＜H/2 的范围内,液体射流的雷诺数Re=5 000,韦伯数We=1 000,其中雷诺数Re=ρU0H/μ,韦伯数We=ρU20H/σ。 气液密度比被设置为ρl/ρg=40,为使气液运动黏度相等,设置μl/μg=40。

图2 展示了不同时刻下气液界面的演化。 本文气液界面通过液体体积百分数为0.5 的等值面表征,图2(a) ～(d)分别对应的时间是0. 01,0.03,0.06,0.10 ms 时刻,对应界面失稳,表面波增长,界面出现破碎,射流初次雾化的时间点。

图2 射流气液界面的演化Fig.2 Evolution of jet gas-liquid interface

1.2 流动拓扑理论

本节介绍流动拓扑分析方法。 根据Perry 和Chong 等[10-12]的研究,速度梯度张量Aij由式(4)给出:

同时,可将速度梯度张量分解为Sij和Wij两部分之和,其中Sij=0.5(Aij+Aji),Wij=0.5(Aij-Aji),分别称为对称张量和反对称张量。 速度梯度张量Aij的3 个不变量由式(5)给出:

显然,第1 不变量P代表流体体积的变化,由质量守恒知P=0,因此本文将主要在Q-R平面内分析流动拓扑结构。 式(5)是一元三次方程,其判别式为

判别式将Q-R平面分为2 个区域,当Δ＞0,方程有1 个实数根和2 个复数根,当Δ＜0,方程有3 个实数根,当Δ=0 时,对应2 条直线,这2 条直线分离不同的拓扑结构区域(见图3):

当2 个复数根为纯虚数时,对应r2=0。 这3条直线r1a、r1b、r2将整个Q-R平面分为4 个区域,分别命名为S1 ～S4,如图3 所示。 利用速度梯度张量的第2 不变量Q可以解释每个部分所表示的物理意义,Q可以被分成2 个部分:

图3 不同拓扑区域分布示意图Fig.3 Different topology area distribution diagram

式中:QS和QW分别为Sij和Wij的第2 张量 不 变量。QS代表每单位质量的动能耗散为热能,代表变形;QW与涡量相关,代表旋转,因此Q＞0,代表系统以旋转为主,Q＜0,则代表系统以变形为主。同时R＞0,代表压缩行为,R＜0,代表拉伸行为。通过流动拓扑分析,可以看出在不同的拓扑区域内,流动受到不同的机制控制。 总而言之,S1 拓扑区域代表不稳定焦点结构(UFC),主要受旋转压缩机制影响。 S2 拓扑区域代表不稳定节点/鞍点结构(UN/S/S),主要受到压缩变形机制影响。S3 拓扑区域代表稳定节点/鞍点结构(SN/S/S),主要受到拉伸变形机制影响。 S4 拓扑区域代表稳定焦点结构(SFS),主要受到旋转拉伸机制影响。

1.3 气液界面曲率

为了研究流动拓扑与气液界面的耦合关系,还需研究气液界面曲率。 下面对界面曲率给出定义。 本文中所考虑的界面是液体体积分数等值面,以α表示。α=1 时表示全为液相,α=0 时表示全为气相。 该等值面的法向量n定义为

图4 使用平均曲率和高斯曲率对等值面结构进行分类[11]Fig.4 Classification of structure of equivalence surface using mean curvature and Gaussian curvature[11]

2 结果与讨论

本节讨论液体射流雾化过程中流动拓扑结构演化以及流动拓扑结构对气液界面曲率的影响。首先,图5 展示了第2 不变量Q和第3 不变量R的联合概率密度分布图(JPDF)。 为方便对比,Q使用QW进行无量纲化(Q/ ＜QW＞),同理R使用QW的3/2 次方进行无量纲化(R/ ＜QW＞3/2)。选取t=0.01 ms 和t=0.10 ms 的结果进行展示。如图5所示,2 个时刻最大的概率密度分布均出现在原点附近,因此将图5 中的结果原点处在图6中进行放大展示。 可以看出,在液体射流初始时刻,S1 拓扑结构出现的概率略微不同,S1 概率稍大,而随着时间推移,到了初次雾化的时刻,S1 拓扑结构出现的概率较大,表明流动主要受UFC 拓扑结构的影响,即在雾化过程中,流体微团主要受旋转压缩作用机制的影响;而且随着雾化过程的进行,这种机制的重要性越发明显。 值得注意的是,不同时刻的图像与之前研究中观察到的水滴形状比较相似,且在雾化后,区域面积变化较小,呈现出一定的自相似性质。

图5 第2、第3 张量不变量的联合概率密度分布Fig.5 Joint probability density function of the second and the third tensor invariance

图6 图5 在原点附近的放大图Fig.6 Fig.5 enlarged view near origin point

为了进一步研究流动拓扑对不同液体体积分数等值面包括气液界面的影响,本文引入条件应变率来对结果进行分析,定义如下:

式中:a为等值面切向应变率,该值的正与负分别对应于拉伸应变率和压缩应变率。

图7 展示了射流雾化过程中应变率在Q-R空间的条件平均分布。 如图7 所示,在UFC(S1)拓扑分区,主要引起拉伸应变率(应变大于0)。而在SFS(S4)拓扑分区,主要引起压缩应变率(应变小于0)。 而在UN/S/S(S2)、SN/S/S(S3)区域,则由应变主导雾化。 另外,图7 表明,初始时刻,应变能主要是拉伸应变率主导。 而随着雾化进行,压缩应变率与拉伸应变率效应共存。 其中,UFC(S1)流动拓扑区域产生拉伸应变率的概率最高。

图7 不同时刻下应变率在Q-R 平面内的条件平均分布Fig.7 Conditional average distribution of strain rate in Q-R plane in different moments

图8 展示了同一时刻t=0.06 ms 下,不同液体体积分数(α=0.2,0.5,0.8)等值面上应变率的差异。 如图8 所示,对于UFS 拓扑结构,当气相占比较大时(α= 0. 2),具有较高的拉伸应变率,而液相占比较大(α=0.8),拉伸应变率较小。对于其他流动拓扑结构,均为产生压缩应变率的概率较高。 其中SFS 结构,压缩应变率较小,而在变形主导的区域(S2,S3),压缩应变率较大。

图8 不同等值面应变率在Q-R 平面内的条件平均分布Fig.8 Conditional average distribution of strain rate in Q-R plane on different iso-surfaces

选取同一时刻t=0.06 ms 下不同的液体体积分数等值面(α=0.2,0.5,0.8) 曲率进行分析,通过κg-κm联合概率密度分布图像探讨不同界面的几何结构。 如图9 所示,等值面呈现出不同的结构。 在气液界面处(α=0.5 处),概率密度分布较为对称,最高概率密度出现在原点附近,说明主要呈现出片状的结构。 而在α=0. 2 时,气相占比较大,平均曲率大于0 概率较大,这种情况下主要呈现出凸形表面结构。 在α=0. 8 时,液相占比较大,平均曲率小于0 概率较大,主要呈现出凹形表面结构,这些表明,局部界面几何结构会受到局部拓扑结构影响很大。

图9 不同等值面上的κg-κm 联合概率密度分布Fig.9 Joint probability density distribution of κg-κm on different equivalence sufales

取不同时刻研究所有液体体积百分数等值面上的κg-κm联合概率密度分布,以此探讨雾化过程的不同时刻整个射流流场的曲率分布,从而获得整个流场界面的大致形状。 除了选取初次雾化之前的t=0.01ms,t=0.03 ms 和t=0.06 ms,还选取了雾化发生后的时间t=0.10 ms,t=0.15 ms和t=0.20 ms。 如图10(a) ～(f)所示,雾化初始时刻主要存在的是凸面、凹面、片状结构,其中片状结构概率最高。 随着雾化进行,出现了管状的等值面结构,同时平均曲率及高斯曲率的绝对值不断增加,这意味着随着雾化进程,也会出现凹双曲面、凹抛物面、凸双曲面、凸抛物面结构,但在每个时刻,仍然是出现片状结构的概率最高。 而在初次雾化后,逐渐呈现出图像不关于κm=0 轴对称的现象,产生凹双曲面的概率略微大于产生凸双曲面的概率。 由图9 可知,液相占比较大时产生凹面的概率较大,可能是由于雾化后期液滴弥散至整个空间所致。

图10 不同时刻下整个流场的κg-κm 联合概率密度分布Fig.10 Joint probability density distribution of κg-κm for whole flow field at different moments

图11 展示了应变率在不同拓扑区域内的概率分布图(PDF)。 可以看出,所有流动拓扑结构都有助于在不同区域内形成压缩应变率和拉伸应变率。 但在每个区域内,产生零应变的概率均为最高。 而随着雾化进行,每个拓扑区域产生的压缩应变率都会大于其产生的拉伸应变率,即偏向于负的应变率。 由图11 可以看出,S1 和S3 的概率密度显著大于S2 和S4 的概率密度。 随着雾化进行,S1 拓扑会产生更大的压缩应变率。 另外,S1 和S3 对拉伸与压缩应变率的贡献最高,由于S1 拓扑比S3 拓扑概率更大,S1 概率增大,应变率变大,曲率变小。

图11 不同拓扑区域内的应变率概率密度分布Fig.11 Probability density distribution of strain rate in different topological regions

图12 展示了平均曲率在不同拓扑区域内的概率分布。 可以看出,所有流动拓扑结构都会产生正平均曲率和负平均曲率,每个区域内,产生零曲率的概率最大,表明产生片状结构的概率最高,这与之前的结论一致。 同时,曲率的概率分布大致关于κm=0 对称,表明在雾化过程中,产生凸界面结构与凹界面结构的概率大致相等。

图12 不同拓扑区域内的平均曲率概率密度分布Fig.12 Mean curvature probability density distribution in different topological regions

图13 展示了t=0.10 ms 时刻应变率和平均曲率的联合概率密度函数分布。 图13 表明,随着应变增加,等值面界面趋向于片状,管状结构,在压缩应变率区域(应变小于0),曲率下降的速度比拉伸区域(应变大于0)快,在压缩应变率区域,应变和曲率的负相关关系更为明显。 这里负相关关系指一个变量随着另一个变量的减少而增加的关系。基于上述结论,可以总结出如下因果链,其中,上下箭头代表该量绝对值的增加或者减少,水平箭头代表因果链:PDF(UFC)↑→应变率↑→曲率↓→雾化性能↓。

图13 整个流场中应变率与平均曲率的联合概率密度函数分布Fig.13 Joint probability density function distribution of strain and mean curvature over the entire flow field

3 结 论

本文对湍流液体平面射流进行了三维直接数值模拟。 通过流动拓扑理论,分析了射流雾化过程中流动拓扑的影响机制。 通过分析气液界面的曲率,揭示了流动拓扑与界面结构之间的关系。从流动拓扑的角度阐述了射流雾化气液界面不断演化的机制。 研究结果表明:

1) 所有流动拓扑结构都有助于产生压缩应变率和拉伸应变率。 射流雾化过程主要受到UFC 拓扑结构的影响,从而产生拉伸应变率,且UFC 的增多会导致更多的等值面界面压缩应变率,进而导致气液界面曲率下降,降低雾化性能。

2) 随着雾化进行,气相区域主要产生凸界面,液相区域主要产生凹界面,气液界面处界面主要为片状结构。 除UFC 外,其余拓扑结构主要产生压缩应变率,在SFS 区域产生的压缩应变率较小,而在UN/S/S 区域和SN/S/S 区域,产生的压缩应变率较大。

3) 随着雾化进行,产生压缩应变率的概率会略微大于拉伸应变率概率。

根据上述结论, 总结出如下的因果链:PDF(UFC)↑→应变率↑→曲率↓→雾化性能↓。该因果链说明,随着UFC 拓扑结构的增加,会使得拉伸应变率增大,界面受到拉伸,曲率则会下降,不利于雾化过程的进行,因此雾化性能下降率。 为了验证以上因果链的普适性,本文作者团队正在进行多参数变化的直接数值模拟研究,旨在全面地揭示不同雷诺数与韦伯数下流动拓扑对液体射流雾化的影响机制。