水滴撞击球面形态与压力的数值模拟

叶 琛，王关晴，郑恩华，徐江荣

(杭州电子科技大学理学院，浙江 杭州 310018)

0 引 言

液滴撞击固体表面的动态行为在工业生产中有着广泛应用，例如喷雾冷却[1]、喷墨打印[2]、燃油喷雾燃烧[3]以及石油化工行业的催化裂化装置[4-6]等，因此，液滴撞击一直是研究热点。Rein[7]和Yarin[8]对液滴冲击固体壁面的理论以及实验研究进行了详细综述，将液滴运动过程分为3个过程，即铺展、回缩、附着。Harlow等[9]采用标记网格法(Marker-and-Cell，MAC)的有限差分技术求解流动方程，对液滴撞击疏水平面进行研究，该方法主要适用于惯性效应占主导地位的冲击铺展阶段。Bussmann等[10]采用欧拉固定网格法对液滴撞击倾斜壁面进行仿真，随着倾斜度增大，液滴的铺展过程持续的时间越长。Roisman等[11]建立了二维轴对称模型，运用动量方程的守恒对液滴撞击平面进行研究，得出预测液滴铺展过程薄膜厚度的经验公式。秦缘等[12]通过改变基底尺寸，研究底面直径对液滴撞击的影响，当底面直径小于液滴直径时，液滴被刺穿并产生飞溅现象。Li等[13]采用基于水平集方法和改进的浸入边界方法数值模型，着重分析了We介于30～120的水滴撞击固体球面后的铺展过程。Monica等[14]采用两相流界面追踪法(Coupled Level-Set and Volume of Fluid，CLSVOF)对液滴撞击加热的固体曲面的流体动力学进行研究，分析了温度对铺展过程的影响。唐鹏博等[15]对液滴撞击球面展开实验研究，结果表明，液滴初始速度以及球面曲率对液滴的铺展过程无明显影响。

上述研究取得了较为丰硕的成果，但仍有不足。首先，在宏观角度能观察到液滴的运动过程，但无法获取运动过程中内部的压力变化，对液滴撞击后运动机理的探究不够深；其次，对于铺展过程的研究偏多，很少涉及回缩过程，特别是中高韦伯数液滴撞击球面及液滴回缩过程的研究甚少。针对上述不足，本文建立三维数学模型，对中高韦伯数下液滴撞击球面进行数值模拟，探究液滴的回缩过程和附着行为，分析液滴运动过程中的压力变化。

1 数理模型与CLSVOF方法

1.1 物理模型

本文数值模拟的物理模型以文献[15]的实验装置系统为基础，仿真中高韦伯数下的水滴撞击干燥球面后的形态变化及附着行为，物理模型如图1所示。假定喷雾水滴在下落过程中保持球形，撞击到球形基底上，水滴降落物理过程的实验图像如图1(a)所示，从中可以看出，水滴的初速度与直径并不会改变水滴运动的整体趋势，因此本文选取与实验工况相同的蒸馏水滴，水滴下降高度H0=20.0 cm，直径d0=3.0 mm，撞击速度u1=2.5 m/s；球形基底的直径D0=20.0 mm，材质为不锈钢，球面的润湿性通过静态接触角给定，静态接触角为65°，如图1(b)所示。

图1 水滴撞击球面的物理模型

1.2 数学模型

计算区域内的不可压缩流动Navier-Stokes方程为：

u=0

(1)

(2)

式中，t表示液滴运动时间，u表示液滴速度，g表示重力加速度，κ表示界面平均曲率,ρ表示计算单元中流体平均密度，ρ=ρg+α(ρl-ρg)，下标l,g分别表示液相与气相；p表示压力，α表示体积分数；μ表示流体的动力粘度，μ=μg+α(μl-μg)；σ表示液体与空气之间的表面张力系数，S表示由连续表面张力模型(Continuum Surface Force，CSF)引起的动力源项，

(3)

采用有限体积法(Volume of Fluid，VOF)对计算单元内气液两相分布进行计算。当计算单元中全部为液相时，α=1；计算单元中全部为气相时，α=0；同时存在液相与气相时，0<α<1，α的控制方程为：

(4)

采用界面追踪法(Level-Set)计算单元边界。引入距离函数φ表示计算单元的中心到界面的距离，计算单元的边界定义为φ=0，计算单元内部的距离函数φ为：

(5)

1.3 网格、边界条件及离散方法

三维模拟的计算域如图2所示。计算域采用底面半径为15 mm、高度为15 mm的圆柱体计算区域。为了更好地反映计算域内部的网格划分，对网格沿YZ轴进行剖面处理，球面采用墙边界，固体球之外的其他区域皆为空气，将其设定为压力输出边界，网格剖视图以及边界条件的设定如图3所示。计算域内部采用六面体结构网格，由于水滴的运动区域主要集中在球面的上方，因此对球面上方进行加密处理，整个计算域中的网格数量达到74万。

图3 网格剖视图

采用有限体积法对控制方程进行离散。针对动量方程的离散采用二阶迎风(Second Order Upwind)格式，VOF方程的离散采用几何重构(Geo-Reconstruct Scheme)格式，同时动量方程中的压力速度的耦合采用Coupled格式；控制方程中对于梯度以及导数的计算，采用最小二乘单元的梯度评价(Least Squares Cell Based)格式；压力的计算采用适用于大体积或者强压力的Presto方法。对于能量方程以及Level-Set方程的离散，均采用收敛性较好的二阶迎风(Second-order-upwind Scheme)格式以保证计算的收敛。为了保证计算的收敛性计算过程中时间步长为10-7～10-6s，每个时间步长内计算45步。

2 模拟仿真与分析

采用CLSVOF方法模拟水滴撞击固体球面，在数值模拟过程中，所有物性与实验测量得到的数据相同：水滴密度ρl=998.2 kg/m3，动力粘度μl=0.001 m2/s；空气密度ρg=1.225 kg/m3，动力粘度μg=1.789 4×10-5m2/s，水滴与空气之间的表面张力系数σ=0.072。球面为干燥的常温球面，接触角为本文1.1节中提到的65°，并且壁面与水滴之间不涉及热量交换，因此，水滴温度、球面温度以及室温均为25 ℃。

2.1 模型验证与形态分析

水滴撞击球面的实验与仿真结果如图4所示，上方为实验图像，下方为仿真图像。

图4 水滴撞击球面的实验与仿真结果

从图4可以看出，当t=0.0 ms时，水滴撞击球面后进入铺展过程；t=2.0 ms 时，水滴逐渐在球面形成一层液膜，水滴铺展的形态以及在同一时间下对应的铺展直径均吻合良好；随着水滴的持续铺展，t=5.2 ms时，达到最大铺展状态，并且在液膜外缘出现明显的“冠指”结构；随后，液膜开始向球面顶部撞击点附近回缩，t=30.0 ms时，液膜重新回缩到撞击点附近，产生向上的中心凸起，随后回缩液膜与球面接触底面基本不变，上部开始震荡；t=55.0 ms时，水滴稳定在撞击点附近。在同一时间下，仿真获得的水滴撞击形态与实验结果基本吻合，铺展半径、润湿面积、稳定附着高度等均吻合。

水滴撞击球面的三维形态运动变化的仿真结果如图5所示。

图5 水滴撞击球面三维运动变化过程

从图5可以看出，水滴撞击初期，在惯性与重力作用下，水滴以撞击点为中心，沿着球面曲率方向向外铺展，并在球体表面形成液膜；当t=5.2 ms时，水滴达到最大铺展状态，在球面摩擦力与水滴自身粘度的影响下，位于液膜边缘部分凸起，形成与主体液膜的高度差，呈“冠指”结构；当t>5.2 ms时，水滴开始回缩；在10.0 ms

2.2 铺展过程的压力变化

水滴撞击球面后，铺展过程的内部压力变化过程如图6所示。

从图6可以看出，t=0.1 ms时，水滴动能瞬间转化为驱动静压能，最大驱动压力位于水滴的外侧；t=0.3 ms时，驱动压力最大值由水滴两侧逐渐转移到液膜中心；t=1.0 ms时，随着水滴的持续铺展，液膜面积增大，但最大压力仍位于液膜中心；t=2.0 ms时，液膜外缘开始凸起，此时为“冠指”结构的雏形，最大驱动压力位置突然阶跃到压力液膜最外端“冠指”内部，驱动液膜继续铺展，“冠指”结构不断增大；t=5.2 ms时，水滴达到最大铺展状态，液膜内部最大压力仍位于“冠指”内部。纵观水滴的铺展过程，随着时间的推移，液膜内部的驱动压力逐渐减小。t=0.1 ms可视为水滴撞击初始时刻，水滴内部驱动压力高达10 000 Pa；当t=2.0 ms时，最大驱动压力仅为280 Pa，驱动压力损失显著，并在随后的铺展过程中持续降低；t=5.2 ms时，最大驱动压力仅为170 Pa。造成驱动压力损失的首要原因是水滴自身的粘性耗散以及水滴与球面之间的摩擦耗散；其次是随着液膜面积的持续增大，水滴内部的一部分静压能转化为液膜的表面能。

图6 铺展过程的水滴内部压力变化

对图6中t为2.0 ms，3.0 ms，5.2 ms中的“冠指”进行局部放大，结果如图7所示。

图7 铺展过程中，“冠指”附近压力局部放大

从图7可以看出，当t=2.0 ms时，基体液膜与“冠指”交界处出现负压，这是由于水滴自身的粘度与表面张力产生向内收缩的趋势。此时，水滴初始动能与重力作用仍主导铺展过程，粘度与表面张力抑制作用并未改变水滴运动趋势，继续向外铺展。当t=3.0 ms时，水滴内部的压力分布趋势和最大压力与t=2.0 ms时的差别不大，无明显改变。当t= 5.2 ms时，在重力和惯性力作用下，铺展液膜与液膜表面张力作用达到平衡，液膜铺展的流动动能降为0，球面液膜表面能达到最大值。在表面张力作用下，下一时刻，液膜开始向撞击中心回缩。

2.3 回缩过程的压力变化

回缩过程中，水滴液膜内部压力分布情况如图8所示。

图8 回缩过程的水滴内部压力变化

从图8可以看出，液膜内部的压力分布中，最大值仍是外缘液膜的“冠指”内压力，并在主体液膜与“冠指”结构的交界处存在压力负值。在这个过程中，随着回缩过程的能量耗散，液膜内部压力不断减小。当t=30.0 ms时，两侧的液膜重新聚集到撞击点附近，使得液膜内部的最大压力较之前时刻略微提升。

回缩过程中，“冠指”的局部放大如图9所示。

图9 回缩过程中，“冠指”附近压力局部放大

从图9可以看出，t=10.0 ms和t=20.0 ms时，“冠指”内部压力明显大于周围液膜处的压力，且最大压力处位于液膜外缘的“冠指”内部，压力由高向低传导，驱使液膜回缩。此过程中，球面的摩擦力和重力对回缩过程起到抑制作用，使得水滴内部的压力持续降低，最大压力由120 Pa降低至110 Pa。水滴自身的粘度与表面张力对回缩的驱动作用表现为“冠指”与液膜边缘所形成的负压力，此时回缩的驱动力转变为水滴自身的粘度与表面张力。且当t=20.0 ms时，由于“冠指”的逐渐聚集，使得撞击点附近的负压达到-120 Pa。当t=30.0 ms时，水滴经过回缩再次聚集在撞击点附近，聚集作用使得水滴内部的压力小幅增加，最大压力由液膜外缘转移到水滴中心。

3 结束语

本文采用CLSVOF方法研究中高韦伯数下水滴撞击固体球面的动态行为问题。通过Ansys Fluent建立了水滴三维运动形态图以及不同时间水滴内部的压力变化，分析了水滴撞击球面后的形态变化趋势以及运动机理。研究发现，受水滴内部静压能的驱动，水滴运动的压力由高向低传导，从而推动水滴的运动。由于本文仿真实验涉及的工况较少，后续将开展不同工况下水滴撞击球面的研究。