喷雾雾滴撞击茶树叶面流场特性的数值模拟与验证*

刘冬梅，陈中哲，周宏平，陈青，杨杭旭

(1. 金华职业技术学院机电工程学院，浙江金华，321017；2. 浙江省农作物收获装备技术重点实验室，浙江金华，321017；3. 南京林业大学机械电子工程学院，南京市，210037)

0 引言

雾滴撞击叶片表面后，会产生多种撞击行为如雾滴粘附、破碎、反弹和飞溅等，而反弹和飞溅存在沉积到非叶面的风险，同时叶面上的粘附小雾滴聚合后也存在滑落及溪流现象，这些都会造成雾滴脱离植物叶面，不利于药液的沉积持留[1]。对于雾滴与叶片间的相互作用，需分为两个部分进行研究，一是雾滴撞击叶面的初始碰撞行为，二是针对撞击后粘附在叶片上的雾滴或雾滴碎片，需跟踪它们的运动以了解它们是否保留以及如何保留，才能更完整地探索雾滴沉积规律[2]。

目前，对雾滴在叶面上撞击行为的研究方法有理论分析、试验研究和数值模拟，因理论分析做了简化和假设，无法反映雾滴运动的细节[3]；而喷雾试验是利用高速摄影仪多次拍摄微米级雾滴，试验成本高且仅能从表观上观察雾滴撞击叶面过程的现象[1]，不能提供全部流场范围的详细信息，对多雾滴撞击的复杂情况，试验及测量更为困难。为直观瞬时显示雾滴在运动过程中全部流场特性，利用数值模拟对不同工况下的雾滴运动进行仿真分析，对完善雾滴撞击机理及运动规律具有重要意义。Roman等[4]利用格子玻尔兹曼方法评估了表面润湿不均匀性对雾滴运动的影响，研究了连续雾滴对固体表面的冲击动力学。Dorr等[5]基于植物模拟软件L-studio预测了雾滴被植物叶截留的情况。Ahmad等[6]利用Lattice Boltzmann研究了两个连续雾滴在平坦表面上的斜碰撞。Fujimoto等[7]在低雾滴冲击惯性下模拟雾滴与基体的正常和倾斜碰撞，研究了冲击角对雾滴变形行为的影响。Mayo等[2]模拟了雾滴在虚拟棉花叶面的运动，以了解农业喷雾剂在植物上的持留。Delele等[8]利用CFD模型研究了雾滴在植物表面的动态撞击行为，预测了撞击、粘附、弹跳或飞溅过程后雾滴的运动。Zhu等[9]利用CLSVOF界面跟踪法研究不同药剂对Shu-ChaZao叶片表面的雾滴冲击行为。近年来，叶片表面特性对喷雾沉积的影响受到日益关注，因不同叶片结构如叶脉、蜡质、毛发、边缘和弯曲梯度等都是重要的撞击和持留变量[10]，雾滴在叶片表面碰撞扩散时，叶片表面会影响雾滴运动，并导致各种雾滴撞击行为[11-13]。

目前模拟研究雾滴撞击行为的文献多以固体表面或以棉花、水稻等农作物叶片作为载体[1, 8, 14-17]。因茶树叶片表面具有自身的表面特征，为完善雾滴在茶树叶面的撞击机理及运动规律，需拓宽茶树叶片雾滴撞击的研究维度将叶片表面特性纳入雾滴沉积特性研究中。本文通过对茶树叶片进行三维扫描，点云重构具有叶片表面形态特征的虚拟叶片，在虚拟叶面上进行雾滴撞击模拟，以探究撞击角、叶片曲率及重力对雾滴流场特性的影响规律。

1 点云重构茶树三维虚拟叶片

茶树叶片三维模型对模拟结果有较大影响，为准确反映茶树叶片表面细节，本文依托逆向工程技术[18]，利用Reeyee蓝光固定式三维扫描仪配套天远三维数字扫描系统获取安吉白茶成熟叶片的点云数据并进行拼接锐化。然后用Geomagic Design X 64软件对点云数据进行精简、去除噪声点等后处理[19]，通过赋厚将二维STL转化为三维STL文件。因Ansys划分网格时会把曲面打断成三角面片造成叶片曲面有一定失真，为避免这个问题，在Ansys SpaceClaim软件中采用不拼接面的转化方法进行STL立体模型导入成实体模型的转化，实现茶树叶面的逆向重构[20]，如图1所示。

图1 安吉白茶成熟叶片逆向建模

2 计算模型及相关参数

2.1 物理模型

模拟仿真设计了单雾滴撞击和多雾滴撞击两种不同的物理模型。对于单雾滴，主要研究其在不同撞击角下的流场变化，因叶片表面具有不同曲率梯度，雾滴铺展情况跟叶片曲率有关，为脱离叶片曲率变量，单雾滴以虚拟叶片微元为研究对象，用具有叶片表面特性(接触角和粗糙度)的平面来表征微元表面。考虑计算区域既要观察雾滴运动，又要减少三维模拟计算量，单雾滴模拟计算区域为3 mm×3 mm×1.5 mm，对于多雾滴流动，采用三维虚拟叶片作为基底进行模拟，其模拟计算区域为7 mm×7 mm×3 mm，如图2所示。初始化利用Patch方法将单个或多个雾滴置入到流场环境中。

(a) 单雾滴试验物理模型

(b) 多雾滴试验物理模型

2.2 数学模型

雾滴在茶树叶面的运动过程，实质上是液相雾滴和气相空气的两相流问题[1, 21]，本文采用多相流VOF模型进行模拟。在雾滴撞击茶树叶面及后续运动过程的模拟中，不考虑相间传质和传热，其控制方程主要包括质量、动量守恒方程，雾滴在叶面上运动时，视为不可压缩流体。质量守恒方程在直角坐标系下可表示为[22]

(1)

式中：ρ——雾滴密度，kg/m3；

t——时间，s；

u,v,w——流体速度矢量在坐标x、y、z方向速度分量，m/s。

动量方程[22]

(2)

式中：p——压力，Pa；

U——流体速度矢量，m/s；

η——流体动力粘度，Pa·s；

Su、Sv、Sw——广义源项，N/m3。

2.3 参数设置及计算方法

模拟试验中采用的喷雾介质为水，室温下气、液两相的主要参数如表1所示，其中气液两相间的表面张力为0.073 5 N/m。水滴在安吉白茶成熟叶片上的静态接触角为77.5°(4 μL液滴)。考虑茶树叶面特征，其表面粗糙度利用粗糙常数Cs和粗糙表面厚度Ks进行表征，其中Cs=0.5，Ks=0[21]。

表1 数值模拟时气液两相主要参数Tab. 1 Main parameters of gas-liquidtwo-phase in numerical simulation

模拟计算区域的上边界设置为压力进口边界条件，其进口压力为0，下边界设置为无滑移壁面边界条件，其余边界设置为易于收敛的出流边界条件。湍流模型选择Realizablek-ε模型，采用有限体积法对控制方程进行离散化，模拟时速度与压力耦合算法选择SIMPLE算法。运算时为减少计算时间，同时又能精确观察雾滴运动，计算中改变运算时间步长，把雾滴在空中的计算步长设为1×10-5s，在雾滴接触界面后，时间步长设为1×10-6s，每个Time step最大迭代步为25。

3 仿真结果分析

3.1 撞击角对雾滴流场特性的影响

雾滴初始参数设置为雾滴粒径D=300 μm，撞击速度V=3 m/s，为分析雾滴入射速度与叶片表面间的夹角即撞击角ψ对流场的影响，将上述雾滴撞击角分别设置为45°、75°、90°。雾滴在下落过程中因受到空气阻力影响，雾滴形状会发生一定变形，雾滴运动方向的空气由于受到推力而具有一定速度。为尽量贴合雾滴运动实际情况，模拟时将雾滴初始位置定义为距离叶面2 mm。不同撞击角时的单雾滴撞击界面的液相分布图如图3～图5所示，当雾滴撞击角为45°、75°时，其在界面撞击行为均为粘附，当雾滴撞击角为90°时，其撞击界面有微小雾滴破碎现象，因破碎过程中雾滴只是由大雾滴变为小雾滴，且破碎小雾滴的铺展过程与正常粘附的雾滴铺展过程一致，故将其视为广义粘附，即三种撞击角下，雾滴都是粘附在界面表面，没有飞溅。

对比雾滴在45°、75°、90°撞击角时的液相分布，发现雾滴从接触界面0 ms到1.25 ms时间段内，都经历铺展、收缩、二次铺展、二次收缩等形态变化[23-24]。当雾滴撞击角为45°、75°时，雾滴在铺展收缩过程中同时向雾滴运动方向滑移，且撞击角为45°时的雾滴滑移距离大于撞击角为75°时的滑移距离，而撞击角为90°时的雾滴无滑移现象，这主要是因为撞击角越小，雾滴接触界面后沿运动方向的切向速度分量越大于法向速度分量，即两者差值越大。当90°撞击角时的雾滴垂直撞击界面时，雾滴向四周同步铺展、收缩，没有滑移产生。因滑移易造成雾滴从叶片滑落，建议喷雾时采用90°撞击角。

(a) 0 ms

(e) 0.65 ms

(a) 0 ms

(e) 0.65 ms

(a) 0 ms

(e) 0.57 ms

45°、75°、90°撞击角时的雾滴内部及周围气流速度分布如图6～图8所示。

下落时，雾滴在空气反冲力作用下其形状发生一定变形，雾滴运动方向的空气被推向两侧，雾滴后方的空气被周围空气迅速补充过来，在雾滴两侧形成涡流[25-26]。

(a) 0 ms

(d) 0.30 ms

(g) 0.95 ms

雾滴在撞击界面前，其内部速度分布均匀，当雾滴撞击界面后，其动能开始转化为表面能，雾滴在界面上向四周铺展，随着雾滴铺展面积增大，雾滴运动速度逐渐减少，当铺展雾滴内部速度为0时，此时在表面张力作用下，雾滴开始反方向收缩。收缩时雾滴表面能转化为动能，雾滴有回到初始状态的趋势，此时雾滴速度方向由向内逐渐转变为向上。在惯性力作用下，雾滴从回缩的最小状态开始二次铺展及回缩，在铺展回缩过程中由于液膜内部粘滞力和表面摩擦力的能量消耗，雾滴铺展回缩的幅度逐渐变小。

(a) 0 ms

(d) 0.30 ms

(g) 1.15 ms

(a) 0 ms

(d) 0.42 ms

(g) 1.12 ms

从45°、75°、90°撞击角时的气流速度分布发现，90°撞击角时雾滴左右两侧铺展速度几乎相同，而45°、75°撞击角时雾滴左右两侧铺展速度存在差异性，且撞击角越小两侧铺展速度差异性越大。从速度分布图中可看出，45°、75°撞击角时雾滴是从右侧射入向界面撞击，铺展过程中雾滴左侧铺展速度因与运动方向一致其速度要大于右侧铺展速度，而在回缩过程中，右侧回缩速度方向因与雾滴运动方向一致其速度要大于左侧回缩速度，造成雾滴在回缩到最高处时，雾滴顶部有向左倾斜的趋势。在二次铺展过程中，雾滴左右两侧速度变化趋势跟首次铺展回缩的趋势相同。撞击角越小，雾滴左右两侧的铺展速度差异性越大，雾滴铺展得越细长。因撞击角对雾滴铺展的影响主要体现在雾滴沿运动方向的扩散，即对椭圆长轴影响较大，对椭圆短轴影响较小[27]，故文中只研究雾滴铺展长轴所在平面即ZX面的速度分布，对于短轴的速度分布不作分析。

3.2 叶片曲率及重力对雾滴流场特性的影响

理论分析及试验研究发现，雾滴粒径D<200 μm的小雾滴几乎都是粘附在叶面上，雾滴铺展面积较小[28]，小雾滴碰撞聚合后的体积仍很小，一般仍是粘附在叶面，很难体现叶片曲率及重力对多雾滴撞击叶面后流场特性的影响。文中多雾滴模拟选取粒径范围为200～400 μm的中等雾滴、粒径D>400 μm的大雾滴两组进行对比分析，每组设置16个雾滴，雾滴间隔1 mm，其中大雾滴初始粒径设置为500 μm，中等雾滴初始粒径设置为300 μm，两组雾滴的撞击速度均为3 m/s，且雾滴沿z轴负方向下落，在x、y方向上的速度分量为0。因撞击界面为茶树虚拟叶片，叶片表面不同位置具有不同曲率，故组内不同雾滴在自身撞击点处的叶片曲率均不相同。为尽量贴合雾滴运动实际情况，模拟时雾滴需距离叶片表面一定距离，文中雾滴所在平面距离xy平面为2.5 mm。

16个300 μm雾滴撞击叶面后的液相分布如图9所示，撞击时每个雾滴均粘附在叶面上且都呈现多次铺展和回缩的形态变化。因茶树叶面有主脉、侧脉及次生脉，这些结构相对于叶面上稀疏绒毛，属于较明显障碍物，其可能会导致雾滴分离，由图9(c)和图9(d)可看出，右侧第一列上部第二个雾滴撞击在主脉上，其在铺展回缩过程中破碎为2个雾滴。从不同时刻液相分布可看出，在不考虑叶脉因素的情况下，初始粒径300 μm的雾滴在叶面上不发生聚合时，其初始铺展收缩主要受惯性力及表面张力的影响，随着能量消耗，雾滴速度减小，叶片曲率对雾滴运动的影响变大[2]，其与雾滴惯性力及表面张力一起决定雾滴的稳定形态，当粒径300 μm的多个雾滴在铺展、回缩过程中接触时则会发生聚结。因粒径300 μm水滴的体积仅为0.014 μL，由此质量产生的重力对雾滴影响比较小[2]，由液相分布发现，没有发生聚结的雾滴主要由叶片曲率和表面张力决定其在叶面上沉积的最终位置及形态，没有因重力而下降到叶片最低处。

(a) 0 ms

(d) 0.6 ms

(g) 1.8 ms

16个500 μm大雾滴撞击叶面后的液相分布如图10所示，雾滴初始铺展为液膜后相互粘连，此时雾滴主要受惯性力和表面张力的影响并没有快速聚结在一起，而是先发生回缩及二次铺展，随着雾滴回缩变大，雾滴之间的聚结变快。粒径500 μm的多雾滴铺展后的轮廓线及运动趋势如图11所示，壁面剪切力的大小和方向反映了雾滴铺展和回缩的运动趋势，其中箭头长短代表壁面剪切力大小。雾滴运动初期其速度比较大，随着雾滴铺展回缩过程中能量的消耗，雾滴壁面剪切力及运动速度逐渐变小，直至雾滴处于相对稳定状态。

(a) 0 ms

(d) 1.2 ms

(g) 3.6 ms

由图10、图11可看出，粒径500 μm的多雾滴在扩散初期其雾滴质量相对较小，雾滴运动主要由惯性力、表面张力和叶片曲率驱动，随雾滴逐渐凝聚变大其质量也同时变大，重力对雾滴运动的影响逐渐体现，在雾滴扩散后期，雾滴质量进一步变大，此时雾滴运动主要由雾滴重力驱动[2]。

由图10(i)看到，聚结的大雾滴在重力作用下沉积到叶片最低处，而右侧小雾滴由于质量轻，主要由叶片曲率和表面张力决定其在叶面上的最终沉积位置及形态。

(a) t=0.4 ms

(c) t=2.0 ms

4 理论及试验验证

4.1 理论验证撞击角影响规律及仿真模型准确性

将水滴在茶树叶面的静态接触角77.5°代入文献[28]中的反弹模型公式，可知水滴在白茶叶面不发生反弹[28]。因静态接触角77.5°<90°，茶树叶面具有亲水性，根据亲水叶面雾滴飞溅临界曲线Kcrit=-0.584(CA20%acetone)+147[29]得茶树叶面的飞溅临界曲线Kcrit=108.4，如图12所示，其中曲线右侧为飞溅区间。根据雾滴在茶树叶面的反弹分析及飞溅临界曲线可知初始粒径300 μm，撞击速度3 m/s的雾滴其撞击结果为粘附。数值模拟时，撞击速度3 m/s的300 μm雾滴在三种撞击角时的撞击行为均为粘附，模拟结果与理论分析具有一致性，仿真模型准确。

图12 茶树叶面的雾滴飞溅临界曲线Fig. 12 Droplet splashing critical curve on the surface of tea leaves

为验证撞击角ψ对雾滴运动特性的影响，将撞击角90°、75°、60°、45°、30°时水滴撞击茶树叶面的参数代入文献[28]中的雾滴铺展面积公式，得到不同撞击角时雾滴铺展面积曲线，如图13所示，由图13中铺展曲线可看出，在雾滴粒径及撞击速度相同的情况下，撞击角增大则雾滴铺展面积也增大，更有利于雾滴沉积，同时发现在雾滴粒径较小时，撞击角对铺展面积影响较小，随着雾滴粒径变大，撞击角对铺展面积的影响也逐步变大。模拟发现，雾滴撞击角越小其发生雾滴滑移的趋势越明显，因滑移存在脱离叶面的风险，模拟结果建议采用90°撞击角。因90°撞击角时雾滴铺展面积最大，从有利于雾滴沉积考虑，理论分析建议喷雾作业时采用90°撞击角。模拟与理论分析得到的撞击角选择结果具有一致性，仿真模型准确。

图13 茶树叶面的雾滴铺展面积曲线

4.2 试验验证撞击角影响规律及仿真模型准确性

利用文献[28]中粘附雾滴的试验数据进行主体效应和相关性分析。通过主体效应分析，可知雾滴粒径D≤0.3 mm时，撞击角对铺展面积显著性P值为0.919，大于0.05，即影响不显著；在雾滴粒径D>0.3 mm时，P值为0.045，小于0.05，说明撞击角对铺展面积影响显著[28]。通过相关性分析得P值为0.270，即整体上撞击角对铺展面积无显著性影响，相关系数为0.145，即撞击角对铺展面积为正相关。因对于D>0.3 mm的雾滴，撞击角对铺展面积有显著性影响，结合撞击角对铺展面积为正相关，即撞击角越大铺展面积越大，故喷雾试验时建议采用90°撞击角。模拟与试验分析得到的撞击角选择结果具有一致性，仿真模型准确。由图13的铺展曲线看出，D=0.10 mm时撞击角对雾滴铺展几乎无影响，当雾滴粒径增大时撞击角对雾滴铺展的影响逐渐明显，试验分析同步验证了理论模型的准确性。模拟时,多雾滴模型中初始粒径为300 μm、500 μm的雾滴在撞击速度3 m/s情况下，均粘附在叶片表面，没有发生飞溅。由飞溅临界曲线Kcrit=108.4可看出，300 μm和500 μm的雾滴在撞击速度3 m/s情况下均在粘附区间，撞击行为与数值模拟中的结果保持一致。

5 结论

1) 撞击角越小，雾滴在叶面上沿运动方向的切向速度分量越大于法向速度分量，即两者差值越大。雾滴斜撞击叶面时，存在滑移现象且撞击角越小雾滴滑移距离越大。因雾滴滑移存在脱离叶面的风险，结合90°撞击角时雾滴铺展面积最大，从利于雾滴沉积的角度考虑，建议喷雾时采用90°撞击角。

2) 阐明了叶片曲率及重力对雾滴流场特性的影响。当撞击雾滴为100～400 μm中小雾滴时，其初始铺展收缩主要受惯性力及表面张力影响，随能量消耗雾滴运动速度减小，叶片曲率对雾滴运动影响变大，其与雾滴惯性力及表面张力一起决定雾滴的稳定形态。当中小雾滴相互接触时会发生聚结，因其体积较小重力对雾滴的影响比较小。对于400 μm以上大雾滴，在铺展初期其雾滴质量相对较小，雾滴运动主要由惯性力、表面张力和叶片曲率驱动，随雾滴逐渐凝聚其质量同时变大，重力对雾滴运动的影响逐渐体现，在铺展后期雾滴质量进一步变大此时雾滴运动主要由重力驱动。

3) 雾滴在茶树虚拟叶面上的流场特性与雾滴撞击茶树叶面的理论分析和试验结果具有一致性，其为完善雾滴在茶树叶面上撞击机理及运动规律提供了新的途径。