基于格子玻尔兹曼方法的双液滴撞击壁面液膜的热特性演化过程

宋翔, 何小泷

(1.重庆交通大学河海学院, 重庆 400074; 2.重庆交通大学西南水运工程科学研究所, 重庆 400074)

液滴撞击覆盖有液膜壁面的现象广泛存在于工业生产和工程建设中,如农业农药喷洒、喷雾冷却、内燃机汽缸内部喷射燃烧等[1-3]。其中喷淋降温被广泛应用于工业过程中,由于其重要应用价值以及复杂的热流耦合特性受到国内外研究学者广泛关注。这一现象通常以液滴群形式出现,其撞击液膜过程中存在液滴与壁面、液滴与液滴之间的相互作用,目前仍难以阐述其在液滴尺度上的热流耦合机理。而开展液滴尺度上的研究,有助于解释其热动力和水动力的演化机理,是研究喷淋降温的基础。目前在液滴尺度的研究中,主要依托高速摄影技术和数值模拟手段,着眼于液滴撞击液膜过程中的水动力过程中的研究,而对壁面与液滴、液膜之间的传热过程的研究则十分有限。

实验研究可以直观的展示出复杂界面的形态变化,但难以提供水动力和热动力演化的连续过程[4-7]。数值模拟作为实验研究的重要补充,可以提供流场与温度场的连续变化过程,有利于对这一现象中的热流耦合内在机理开展系统研究。研究者们对液滴撞击液膜过程中的传热现象开展了大量研究。Liang等[8-10]利用耦合Level-set和VOF的方法(CLSVOF)对单液滴和多液滴撞击不同状态液膜的热动力过程开展了系统研究,这些研究着眼于液滴直径、液滴排列、碰撞速度等相关几何参数对传热过程的影响,并未涉及液相物理性质(黏滞系数、表面张力等)对其热动力演化的影响机制,其研究结果表明:单液滴撞击覆盖液膜壁面的传热过程中,壁面热流密度与撞击速度呈正相关,液膜厚度和液滴直径对壁面传热过程的影响较小。Wang等[11]和郭亚丽等[12]在此基础上,分别对双液滴撞击移动液膜的不对称性和液滴连续撞击恒温壁面上的传热特性进行数值研究。研究表明液滴撞击速度、液膜厚度、双液滴间距是影响壁面热流密度的关键因素。

格子玻尔兹曼方法(lattice Boltzmann method,LBM)作为一种有效的介观数值模拟方法,通过求解一系列线性方程,避免了对传统宏观方程中非线性对流项的处理,兼顾其在并行运算和复杂边界处理上的优势被广泛应用于多种复杂多相流现象的数值模拟[13-17]。其中格子玻尔兹曼伪势模型因其独特的简洁性和计算高效性优势脱颖而出,通过粒子间伪势作用,自动形成两相界面,无需进行界面跟踪,被广泛应用于带有传热传质过程的复杂多相流现象模拟中,包括空化、液滴撞击液膜、结晶、蒸发和沸腾等复杂多相流现象中[18-20]。现采用双分布方程热流耦合伪势LBM,对并排双液滴撞击液膜过程中的水动力和热动力过程开展了系统研究,系统探讨液相黏滞系数、双液滴间距、撞击速度等因素在撞击过程中对壁面瞬时热流密度分布的影响。

1 计算模型与方法

1.1 格子玻尔兹曼伪势模型

由于原始的伪势模型(Shan-Chen模型)不满足热力学一致性且数值稳定性较差,故采用多松弛时间(multi-relaxation time, MRT)伪势模型对流体运动过程模拟,可实现大密度比两相流模拟,兼具良好数值稳定性[21-22]。带有MRT碰撞算子的粒子分布方程可以描述为

(1)

(2)

式(2)中:c=1为格子常数。根据粒子分布函数可求得各节点上宏观密度ρ和实际速度u为

(3)

式中:F为作用在粒子上的合力。Li等[21]提出的附加源项S被引入本研究中保证大密度比和大黏滞系数比条件下的数值模拟稳定性和热力学一致性,其表达式为

(4)

式(4)中:ux和uy为速度在x和y方向上的分量;Fx和Fy分别为合力F在x和y方向上的分量;ε为用于调整热力学一致性的参数;τe为与系统总能相关的松弛系数;Fm[22]为粒子间作用力,其计算式为

(5)

式(5)中:G=-1,为粒子间作用强度;ωi为权重;ψ为粒子间伪势,可表示为

(6)

(7)

式(7)中:T为温度;在本研究中,各参数分别选用a=1,b=4,R=1。

Λ为由松弛系数组成的对角矩阵,即

(8)

为了实现可调表面张力,引入源项C,即

(9)

式(9)中:Qxx、Qyy、Qxy作为源项C的相关离散项可以被表达为

(10)

式(10)中:κ为用来调节表面张力的参数;G为作用强度,取值-1;ψ(x)和eα分别为对应位置的伪势和α方向速度分量;ωi表示对应于i方向上的权重因子。

1.2 LBM热动力模型

本文中粒子分布方程所对应的宏观温度控制方程为

∂tT+∇·(uT)=∇·(α∇T)+φ

(11)

式(11)中:∂tT表示温度在i方向上的偏导;α=k/ρcv为热扩散系数;k为热导率;cv为比热容;φ为源项,表达式为[24]如下

(12)

其对应的LBM热粒子分布方程为

(13)

O=(O0,0,0,0,0,0,0,0,0)T

(14)

利用矩阵M,可将式(13)投影到矩空间

(15)

neq=T(1,-2,2,ux,-ux,uy,-uy,0,0)

(16)

此外,n*中两个修正项[25]改写为

(17)

(18)

2 模型验证

2.1 Laplace定律和热力学一致性

Laplace定律定义当液滴达到平衡状态时,液滴内外压差和半径成反比例关系,即ΔP=δ/r0,ΔP为液滴内外压差,δ为表面张力,不同κ值对应直线的斜率表示不同的表面张力。本文中选取5种液滴初始半径r0=35,40,45,50,55在5个κ值条件下对模型进行验证,模拟结果如图1所示,本模型模拟获得的液滴内外压差与其平衡时刻的半径倒数具有良好的线性关系。此外,本模型热力学一致性与Maxwell理论解对比,如图2所示,本文模型模拟得到的共存密度与麦克斯韦理论解在液体和气体分组中一致。

图1 不同κ值条件下液滴内外压差与其平衡态半径倒数的关系Fig.1 Relationship of pressure difference outside and inside the droplet to the inverse of its equilibrium radius under different κvalues

Tc表示临界温度图2 LBM模型理论解与Maxwell理论解的热力学一致性验证Fig.2 Verification of thermodynamic consistency between the theoretical solution of the LBM model and Maxwell’s theoretical solution

2.2 D2定律

本研究模型的热动力过程由D2定律验证。当液滴在静止无限域中蒸发时,其无量纲直径平方[D(t)/D0]2会随时间t线性减小。将直径D0=80 lu的液滴放置在400 lu×400 lu正方形计算域中心,密度场和温度场的四周边界均采用周期边界,不考虑浮力和黏性散热,导热系数定为常数k=2/3,液滴温度设为Tl=0.86Tc。周围气相温度设为Tg=Tc以触发蒸发条件。如图3所示,无量纲直径平方[D(t)/D0]2随时间t线性减小,其中D(t)为t时刻液滴直径。数值模拟结果表明液滴半径减小过程符合D2定律。

图3 无量纲液滴直径[D(t)/D0]2与时间的关系Fig.3 The relationship between dimensionless droplet diameter [D(t)/D0]2 and time

2.3 单液滴撞击液膜形态变化

进一步利用液滴撞击液膜过程中的水动力演化对本模型进行验证,矩形计算域大小为1 000 lu×300 lu,计算域底部为无滑移边界,左右侧为周期边界,顶部设置为非平衡外推边界。其中计算域温度、壁面温度和边界温度均设置为恒温边界T=0.50Tc,对应密度比为ρl/ρg=720。气相黏滞系数为τg=1.062 5,液相黏滞系数为τl=0.537 5,其对应黏滞系数比为υl/υg=15。

液滴撞击液膜过程中各无量纲参数的定义为

(19)

式(19)中:vl和σ分别为液相黏度和表面张力;ρl为液相密度;本验证算例中U0=0.125 lu/ts为液滴碰撞速度,H=25 lu为液膜厚度,r0=50 lu为液滴半径;h*为无量纲液膜厚度,t*为无量纲时间,本研究中对应的雷诺数Re=500,韦伯数We=87.8,研究中无量纲液冠半径和无量纲液冠高度则定义为

(20)

式(20)中:rc表示液冠半径;h为液冠高度。根据Gao等[26]提出的半经验模型,液冠无量纲半径的演化过程是液膜厚度和无量纲时间的函数,即

(21)

式(21)中:β=(2λ2/3h*)1/4,由能量损失系数λ=0.26/Re-0.05We0.07h*0.34确定。

模拟结果与理论计算值对比如图4(a)所示,在无量纲时间t*=0.4～1.4期间,两者吻合良好,而在t*=1.4～2.0期间,受液冠破碎的影响,两者存在较小的差异。进一步对其液冠在垂直方向的延展率进行定量对比。Gao等[26]研究表明,液冠垂向延展率为

St=(t*+ti+2τ0)-1

(22)

○为液冠破碎前演化过程 ●为液冠破碎后演化过程图4 理论结果与模拟结果之间的对比Fig.4 Comparison between the numerical results and theoretical analysis

模拟获得延展率与理论求解的对比如图4(b)所示,模拟结果在t*=0.4～1.4与理论解吻合良好,而在t*>1.4后,受到液冠尖端液滴脱落的影响,与理论解之间出现差异。

3 结果与讨论

3.1 物理模型及参数设置

如图5所示,并排分布双液滴撞击液膜的计算域模型示意图中计算域大小为1 200 lu×500 lu,初始时刻将并排双液滴处于和液膜表面相切的位置,双液滴圆心距L0,液滴直径D0=2R=100 lu,双液滴关于横坐标中心轴(x=600)对称分布。液膜厚度h0=30 lu,以相同初始速度U0垂直向下撞击液膜。本研究中忽略重力作用,气液两相密度比为ρl/ρg≈750,与290 K条件下水与空气密度比接近。水动力计算边界条件中,上边界采用非平衡外推边界格式,下边界采用无滑移边界,左右边界采用周期边界。而热动力计算中四周均采用恒温边界,其中顶部边界和左右边界均为0.50Tc,而底部边界高温壁面恒定温度为Twall=0.70Tc,计算域内其余初始温度均设为0.50Tc。

图5 双液滴撞击液膜示意图Fig.5 Schematic diagram of double droplets impinging on liquid film

本研究中壁面各节点瞬时热流密度计算公式为

(23)

式(23)中:n为壁面法向方向;ΔT为第一层网格和壁面之间的温差;Δy为网格高度。

3.2 双液滴撞击液膜的动力学和热力学特征

如图6所示,并排分布双液滴在初速度U0=0.15 lu/ts,液滴间距L0=1.0D0,液滴直径D0=100 lu的条件下撞击覆有液膜的高温壁面时温度场和速度场的演变过程。内侧双液冠碰撞相融后,部分动能转化成热能,故在t*=0.75时观察到中心射流处存在局部高温,对应中心射流处的高流速区域。两内侧液冠撞击后,导致部分表面低温液体下潜,增大了壁面与流体之间的温度梯度,同时中心射流处的高流速促使流体与壁面之间的对流传热增强,进而在中心射流处对应壁面上出现高热流密度区,其峰值达到1.19 mu/ts3,如图7所示。同样,双液滴撞击形成的外侧液冠附近出现高流速区域导致对流传热增强,同样促使其对应位置的壁面对流传热增强,在其瞬时热流密度曲线出现对称分布的两个峰值,热流密度大小略小于中心峰值,达到1.14 mu/ts3。值得注意的是在t*=0.75时,壁面热流密度曲线无论是静态区和撞击区热流密度值均处于整个过程中的最大值,这是由于在初始撞击阶段,扩散传热占据主导地位。初始阶段无论撞击区还是静态区,液膜尚未被壁面加热,两者之间温度梯度较大,也是这一时刻热流密度较大的主要原因。

图6 双液滴撞击湿润高温壁面过程中温度场和速度场演化图Fig.6 Evolution during double droplets impinged on a high temperature wall covered with liquid film temperature field and velocity field

图7 双液滴撞击湿润高温壁面过程中壁面热流密度随时间的变化Fig.7 Variation of wall heat flux with time during double droplets impingement on a wetted high temperature wall

随着撞击液滴与液膜的持续性相互作用,在t*=1.50时,外侧液冠末端因Rayleigh-Taylor不稳定性发生破碎,形成二次液滴。中心射流高度持续上升,伴随着射流厚度变薄,同样产生二次液滴脱落。中轴线附近液膜内流速发生碰撞后迅速衰减,导致对流换热强度减弱,中轴线上热流密度峰值也随之减小。而撞击后低温液滴扩张和下潜导致壁面附近的温度梯度增大,撞击点附近的热流密度峰值随之减小,但其峰值范围则随着液冠半径的增大而增大。且由于液冠的不对称演化,在撞击点两侧其热通量呈现不对称特性。对比t*=1.50, 2.25, 3.75时热流密度瞬时分布曲线,在撞击区域,由于撞击的低温液滴与壁面附近高温流体的持续性热交换,导致其热流密度随着时间增加而增大。而在静态区,受到速度不连续性的影响,低温液滴与高温液膜之间的热交换仅靠扩散传热。随着液膜被壁面加热,两者之间的温度梯度逐渐减小,热流密度也随时间减小。

如图8所示,双液滴间距分别为1.5L0和2.0L0以撞击速度U0=0.15 lu/ts撞击液膜过程中在t*=0.75,2.25,3.75时壁面热流密度分布曲线。当双液滴间距L0=2.0D0时壁面热流密度和高热流密度区域面积更大。这是因为双液滴间距较大时,受内侧液冠碰撞相融带来的速度不连续性的影响,低温液滴撞击液膜后铺展范围随着两液滴之间的间距增大而增大,导致大温度梯度的区域随之增大。内侧液冠出现双液冠碰撞相融的时间延后,液膜吸收更多液滴带来的动能,液冠碰撞损耗动能更少,对壁面液膜产生更大扰动,增强对流传热。故壁面热流密度最大值随着间距增大而增大,同时也增大了高热流密度撞击区域面积。

图8 不同液滴间距对壁面热流密度的影响Fig.8 Influence of different droplet spacing on wall heat flux

研究表明,双液滴撞击湿润高温壁面期间,壁面热流密度分布随着速度增大而增大[27-28]。并排分布的双液滴以不同撞击速度U0=0.10 lu/ts和U0=0.15 lu/ts撞击液膜过程中t*=0.75,2.25,3.75时壁面热流密度变化曲线如图9所示。双液滴撞击液膜过程中的温度分布如图10所示。结果表明,在保持液膜厚度、双液滴间距不变的情况下,撞击区壁面热流密度随双液滴撞击速度增大有明显增长的趋势。当撞击速度U0=0.10 lu/ts时,中心射流对应位置的壁面邻近网格节点上的液膜温度高于撞击速度在U0=0.15 lu/ts时的温度分布,这说明在t*=0.75时,撞击速度U0=0.10 lu/ts对应的中心射流处液冠撞击后形成的低温下潜流动未能影响到壁面附近的高温流体。因此,当撞击速度U0=0.10 lu/ts在t*=0.75时刻,与U0=0.15 lu/ts工况相比,此刻在中心射流位置未形成壁面高热流密度区域。随着撞击速度增大,液滴的初始撞击动能更大,这也导致液冠处不连续性增强,液冠动能增大,进而导致中心射流溅射高度更高。当U0=0.15 lu/ts时,外侧液冠末端和中心射流尖端均发生二次液滴脱落。但随着撞击速度减小到U0=0.10 lu/ts,外侧液冠和中心射流未能克服表面张力和黏性力发生断裂形成二次液滴。更大初始速度的低温液滴撞击后,其在液膜中铺展范围更大,同时也造成低温流体下潜更深。这增强了液膜中的对流传热强度,进而导致壁面附近流体温度降低,温度梯度随之增大,撞击区壁面热流密度峰值增大。但对于静态区,由于液滴撞击后在液冠附近液膜中存在速度不连续性,这导致液滴撞击后动能无法传递到静态区,此时静态区的传热仍以扩散传热为主。因此不同撞击速度条件下,液膜静态区的温度分布相似,同时其热流密度分布曲线基本重合。

图9 不同撞击速度对壁面热流密度分布的影响Fig.9 Influence of different impact velocities on wall heat flux distribution

图10 不同撞击速度撞击后液膜温度分布演化Fig.10 Temperature distribution evolution of liquid film after impact at different impact velocities

黏滞系数增大会直接影响到液滴撞击后的下潜深度和黏滞耗散。图11展示了不同黏滞系数条件下壁面热流密度分布随时间的变化。研究表明,液冠扩散过程中的液冠半径并不随着黏滞系数ν的变化而变化[29]。这导致不同黏滞系数条件下,其相同时刻的热流密度分布曲线形态相似。但在更小的黏滞系数条件下,其在液滴撞击后黏滞耗散掉的动能更小,液滴下潜深度略大。这将增加液膜内的对流传热强度和壁面附近的温度梯度,故黏滞系数较小时,壁面热流密度略大。

图11 不同黏滞系数对壁面热流密度的影响Fig.11 Effect of different viscosity coefficients on wall heat flux

4 结论

利用耦合水动力伪势模型和LBM热动力模型,考虑流场和温度场的耦合作用。研究表明,本文提出的模型可以有效模拟液滴撞击液膜这一形态演化过程和液膜与壁面的传热过程。进而对并排分布双液滴撞击湿润高温壁面过程进行模拟,探究了液滴间距、撞击速度和黏滞系数对不同时刻瞬时热流密度分布的影响,得到以下结论。

(1)水平分布双液滴撞击湿润高温壁面时,会引起液膜内部径向流动,改变其速度场,从而增强壁面与液膜的对流换热过程,造成液滴撞击区域壁面热流密度骤增。而在静态区,由于速度不连续性,液膜内的传热过程仍以扩散传热为主。

(2)撞击速度增大,引起的液膜扰动强,液滴下潜和扩展程度增大,对流传热效果增强,瞬时壁面热流密度峰值也随之增大。而当液滴间距适当增大时,内侧双液冠碰撞相融时间延后,低温液滴在液膜内扩展程度增大,导致瞬时高壁面热流密度面积增加。

(3)黏滞系数虽然不影响液冠半径发展,但其增大了液滴撞击液膜后的黏滞耗散,也影响到低温液滴的下潜程度和撞击区的对流强度,壁面热流密度峰值随着黏滞系数的增大而减小。