在计算中假定气液流动过程是等温的且不考虑蒸发过程,因此无需求解能量方程。求解的不可压、两相流、含有表面张力的N-S方程组为:

{ρ((u^)/(t)+u^·SymbolQC@u^)=-SymbolQC@p+SymbolQC@·(2μD)+σκδ_sn^

(ρ)/(t)+SymbolQC@·(ρu^)=0

SymbolQC@·u^=0

式中:ρ≡ρ(x^,t)为流体密度; u^=(u,v,w)为网格内流体速度; μ≡μ(x^,t)为动力粘度; D_为应变张量; D_ij≡1/2((u_i)/(x_j)+(u_j)/(x_i)); δ_s为集中在界面上的表面张力; κ和n^分别为界面的曲率和法向方向。

VOF方法最早由Hirt和Nichols^[16]提出,对于两相流动,通过定义第一相体积分数α(x,t)来描述气液界面。将密度连续方程用体积分数α的连续方程替换为:

(α)/(t)+SymbolQC@·(αu^)=0

密度和粘度也可由体积分数表示为:

{ρ(α)≡αρ₁+(1-α)ρ₂

μ(α)≡αμ₁+(1-α)μ₂

式中:ρ₁,ρ₂和μ₁,μ₂分别为第一相和第二相的密度和粘度。

Gerris计算三维算例的网格均为正方体单元,如图1所示每个Box由前、后、左、右、上、下共六个面组成,计算域可以根据计算模型由多个Box块连接而成。计算区域可以不用三维建模软件建立,直接在参数脚本文件中给定Box数目、Box连接方向和计算模型在计算域的位置。为了提高计算精度,文中采用八个边长为11.25 mm的Box,计算域如图2所示,上下两层,每层4个Box,计算域长、宽、高分别为45 mm,11.25 mm和22.5 mm,横向气流以一定的速度从两层最左侧Box的左侧面流入,液体射流以一定角度从下表面上的圆孔射入。计算时最高网格等级设为10级,最小网格为10.99 μm,计算域和边界条件的设置均通过编写大量代码实现。

图1 基本块Box离散方向的定义
Fig.1 Definition for discrete direction of basic block Box

图2 计算域
Fig.2 Calculation domain

射流破碎过程仿真边界条件如表1所示,算例1和算例2属袋式破碎,在算例1与文献[1]中同工况实验结果对比验证计算方法准确的基础之上,算例2捕捉射流柱表面细节特征。算例3～算例5属剪切破碎,在与文献[2]中同工况实验对比验证后,分析射流角度对破碎效果的影响。

横向气流中液体射流破碎过程的数值模拟属于多相流模拟的范畴,在计算过程中需要给定的物性参数有液相的密度、粘度系数、表面张力系数、气相的密度及粘性系数等参数。数值仿真中液相为水,气相为空气。物性参数如表2所示。

表1 边界条件
Tab.1 Boundary conditions

表2 物性参数表
Tab.2 Physical parameters

图 3给出了算例1工况下数值仿真结果和文献[1]中实验结果对比图,计算所得到的雾场宏观结构与实验结果非常相近。

图3 算例1的计算结果与试验结果对比
Fig.3 Comparison of numerical calculation and experiment results of Case 1

图 4为算例3～算例5工况下的数值仿真结果与文献[2]中实验结果的对比图,图4中颜色代表速度大小,从图4中可以看出,随着射流角度的增大,射流柱纵向穿透深度增大,横向破碎长度减小,由于气流速度与射流柱横向分速度差变大,射流柱上的波动和惯性力更加明显,雾化效果较好。当角度较小时,射流柱上的气动力相对作用将减小,射流柱变直,接近于自由射流。

图4 算例3～算例5的计算结果与试验结果对比
Fig.4 Comparison of numerical calculation and experiment results of Case 3, Case 4 and Case 5

为了进一步验证数值计算的准确性,本文取射流柱外围为射流轨迹,将数值计算得到的射流轨迹结果与实验结果对比,图5(a)和图5(b)分别为袋式破碎和剪切破碎射流轨迹图,数值计算结果与实验结果相对误差在5%左右。

图 6为射流柱在横向气流作用下变形、弯曲到进一步破碎成大液滴的破碎全过程,由于采用VOF相界面捕捉方法,显示结果是体积分数为α=0.5的等值面。从图6中可以看出,射流柱在初始发展阶段没有明显的液滴剥离,只在射流柱表面有小褶皱,随着射流柱的发展,射流柱表面的小褶皱在气动力作用下剥离形成小液滴,接着迎风面向内凹陷形成波谷,随着射流柱的变形,迎风面和背风面间的波谷厚度不断变薄,有袋式破碎的特征,最终在表面张力和气动力共同作用下断开,形成小液块和液滴。

图5 液体射流轨迹
Fig.5 Trajectories of liquid jets

图6 袋式破碎过程数值计算结果
Fig.6 Numerical calculation results of bag breakup process of jet column

为了更清晰的观察射流柱表面特征和截面特征,对算例2中射流柱弯曲前截面进行分析。将三维雾场沿着射流柱所在的xy平面和yz平面剖切,得到了前视图和左视图的射流柱二维体积分数分布云图,如图7所示。图7(a)显示了射流柱在横向气流的作用下沿射流方向越来越窄,图7(b)显示了射流柱沿射流方向逐渐变宽,迎风面稍有变大。这是因为横向气流遇到射流柱发生滞止,在迎风面上形成一个滞留区,滞留区内速度为零,使得迎风面压力增大,同时射流柱背风面压力较低。在压差作用下,射流柱变得宽而薄,随着射流柱变形,射流柱受到气流作用力越来越大,导致射流柱在流动方向表现为向下游弯曲。

为了进一步认识射流柱破碎机理,图8是不同视角观察到的射流柱弯曲前迎风面放大图,可以清晰看到射流柱表面形成的表面波,表面波是导致射流柱不稳定的根本原因。由于射流柱左侧迎风面受到气流垂直作用,气液密度差和交界面上的垂直加速度导致RT波的产生,RT波将射流柱两侧液体向前推移,以致射流柱变得薄且向下游弯曲。横向气、液相对速度由1逐渐变小,气液横向速度差导致KH波的产生。横向气、液相对速度差也是影响射流柱破碎的重要参数。箭头指出了射流柱发展过程中在射流柱表面产生的表面波,1为射流柱在射流发展方向与气流存在速度差通过相互剪切而生成的KH表面波,2为射流柱在横向与气流存在速度差产生的KH表面波,随着KH波的发展射流柱侧面有液滴剥离,从图中也可以清晰的看出表面波处颜色较深,即表面波处速度较大。

图7 体积分数云图
Fig.7 Contours of volume fraction α

图8 射流柱表面速度分布图
Fig.8 Contours of velocity distribution on jet column surfaces

以30°为例分析射流柱剪切破碎过程,图9是射流角度为30°工况下的射流破碎过程图,在射流柱刚喷入横向气流中的初始阶段,即有小液滴从射流柱表面剥离,这一现象始终贯穿于射流柱发展的全过程。射流柱迎风面和背风面均出现小尺度波动,但迎风面表面波尺度大于背风面表面波尺度,这是由于迎风面受到的气动力更大。小尺度波动在气动力作用下发展为明显的大尺度波动,最终由于气动力大于射流柱表面张力和粘性力而导致射流破碎。这类似于射流柱袋式破碎在气动力作用下迎风面向内凹陷形成波谷并被拉长后断裂的特征。同时,随着射流柱的发展,射流柱逐渐破碎成许多大液块,发生一次雾化,而不像射流柱袋式破碎先断裂为单一大液块后在气动力的作用下发生二次雾化破碎为众多小液滴。在整个射流柱破碎过程中,射流柱剪切破碎和袋式破碎一样存在R-T和K-H两种不稳定波,但剪切破碎的不稳定波更加明显。

图 10～图 12分别表征的是射流中间截面的体积分数云图、速度云图和涡量云图,从图中可以看出射流在出口处即有小液滴剥离,射流柱迎风面在气动力的作用下出现小尺度波,同时可以看到在射流柱迎风面很薄的区域内速度很小,这是因为横向气流遇到射流柱发生滞止,速度减小,压强增大,背风面比迎风面速度大,压强小,迎风面的压强大于背风面的压强使射流柱弯曲。随着射流柱的发展,非稳态波动变大,并形成大量涡结构,加剧雾化过程,缩短射流柱破碎长度。

Gerris可以对全场所有的液滴数目、体积进行统计后并输出。算例3～算例5工况下的全场液滴SMD分别为60.23 μm,59.96 μm和59.12 μm,射流角度对粒径影响很小。图 13是射流角度为30°时计算的全场液滴粒径分布图,从图 13中可以看出,液滴粒径分布范围广,但呈现两端少、中间多的分布。

图9 射流柱剪切过程数值计算结果
Fig.9 Numerical calculation results of shear breakup process of jet column

图 10 射流柱中心截面体积分数云图
Fig.10 Volume fraction contour on center section of jet column

图 11 射流柱中心截面速度云图
Fig.11 Velocity contour on center section of jet column

图 12 射流柱中心截面涡量云图
Fig.12 Vorticity contour on center section of jet column

图 13 全场液滴粒径概率密度分布
Fig.13 Probability density distribution of droplet diameter in whole flow field