农药雾滴空间运行中的变形特征分析

马学虎，薛士东，孙桐，奚溪，宋小沫，赵峻逸，兰忠，郝婷婷

（辽宁省化工资源清洁利用重点实验室，大连理工大学化学工程研究所，辽宁大连116024）

化学农药在重大病虫草害防控、保证国家粮食安全方面发挥着不可替代的作用，合理使用农药有助于提高农业生产力和保障农作物质量，使用农药可减少全球高达45%的粮食损失[1]。据国家统计局数据显示，2014 年中国农药的使用量为180.69×104t，但是由于喷雾技术和喷雾管理的限制，农药的利用率仅为20%～30%[2]。

农药的使用过程一般分为以下几个步骤[3-4]：①药液通过喷头完成雾化过程；②雾滴在空气中运动扩散，沉积到靶标冠层或飘失到靶标外；③雾滴在冠层表面撞击、铺展；④雾滴被冠层表面吸收，传导到为害部位并发挥药效。根据ISO 22856 标准所定义，雾滴飘移是指“在施药过程中，农药雾滴被气流携带脱离靶标冠层、流失到地面或者悬浮在大气中的行为”，被认为是导致农药污染的“罪魁祸首”，对周围河流中的水生动植物、大气环境、人类身体健康造成不可扭转的危害[5-10]。雾滴飘移是诸多因素共同作用的结果，包括药液的理化性质、雾滴的粒径及分布、作物结构（排与排的间距、冠层结构、冠层穿透性）、操作条件（喷雾流量、操作压力、农机行驶速度）、气候条件（风速、风向、温度、湿度）等[11-18]。

研究雾滴飘移的手段主要包括风洞试验、田间试验、数值模拟等。在风速可调、温湿度可控的低速风洞中，研究喷头类型、操作压力、助剂种类等对雾滴群飘移比例的影响[19-21]。在田间试验中，通过检测靶标冠层内示踪剂的含量来评价雾滴群的飘移沉积比例[22-24]。在数值模拟中，采用欧拉-拉格朗日离散相模型追踪雾滴群的运动轨迹，从而获得特定条件下雾滴群的飘移率[25-27]。

以往的数值模拟研究中多将农药雾滴视为不会发生形变的刚性颗粒，忽略了雾滴的形变过程。在雾滴撞击靶标和界面蒸发的研究中也多以完美球状作为雾滴的初始状态[28-30]。然而大量研究表明，雾滴形状对雾滴的聚并与破碎、靶标撞击、所受的曳力大小等都有着重要作用，会直接影响雾滴的飘移损失与界面流失等。QUAN 等[31]利用耦合了移动网格相界面追踪的三维有限体积交错网格法模拟了液滴在加速气流中的运动情况，分析了液滴的形变状态及附近气流的速度场变化。结果表明，变形过程中液滴的非稳态曳力总是比同等雷诺数下液滴的稳态曳力大，同时降低表面张力和黏度会增大液滴的形变量，导致曳力增加，而密度比对液滴动力学影响不显著。QU 等[32]用ANSYS Fluent 软件对非对称流中的减速液滴进行模拟，结果表明随着We 的增大，液滴的形变量及曳力会急剧增大，振荡周期延长；随着Oh 的增大，形变量与曳力略有增大而振荡周期几乎不变。

目前针对单个农药雾滴在风场中的运动研究鲜有报道，尤其是对微米或者亚毫米级别农药雾滴的空间变形特征更是缺乏研究。单雾滴是构成雾滴群的基础元素，明晰单雾滴的运行飘移规律，对于深入研究整个雾滴群的运动规律具有重要的指导意义。因此本文以农药喷雾中的单雾滴为研究对象，采用数值模拟的手段研究了雾滴在横向风场中的空间形变特性，揭示了不同因素作用下单雾滴的形态演化规律，并建立了多因素协同下雾滴形变量与形变周期的预测模型。

1 数值模拟计算模型与方法

1.1 几何模型

基于COMSOL Multiphysics 5.2建立数值模拟过程，选择带相初始化的层流-水平集物理场，对气液相界面进行追踪，获得雾滴的空间形变特性。如图1 所示，针对直径为100µm 的雾滴，选取1200µm×300µm二维矩形计算域，雾滴的初始位置坐标为(100µm，1100µm)。计算域的尺寸随着雾滴的大小进行相应调整，以保证计算域的合理利用并尽可能降低计算成本。

图1 几何模型

1.2 控制方程及边界

对含有分界面的不可压缩、不混溶的气液两相流动系统进行模拟时，关键问题在于如何将界面离散化并且准确地追踪界面的演化过程。Level Set方法通过设置一个距离函数隐式追踪自由面，并且将界面条件融入控制方程中，可在整个域内用统一的差分格式进行求解。与移动网格方法和相场方法相比，可以避免显式地追踪界面，不需要进行复杂的界面重构，可方便地追踪界面的拓扑结构改变，两相控制方程组如式(1)、式(2)所示。

连续性方程

Navier-Stokes方程

式中，Fσ为单位体积流体所受的表面张力；Fd为单位体积流体所受的曳力。由于气相密度远小于液相密度，因此忽略雾滴所受浮力的作用。

雾滴在横向风场中所受的曳力大小如式(3)所示。

式中，A为雾滴在风向上的投影面积，在二维计算中降维为投影长度；ud为雾滴在空间某一位置处的速度矢量；ua为相应位置处空气的速度矢量；Cd为曳力系数，计算依据选取文献[33]中提供的农药喷洒过程中雾滴雷诺数在0.25＜Re＜104范围内的经验公式，如式(4)所示。所受曳力方向与雾滴的相对运动方向相反。

表面张力可以在整个计算域中连续求解，表达为一个光滑函数的形式，见式(5)。

其中，κ为相界面曲率，其值见式(6)。

n为气液相界面法向向量，其值见式(7)。

φ为经Heaviside函数平滑后的水平集函数，代表网格内空气相的体积分数，其在界面一侧的值恒为1，另一侧恒为0，在相界面很小的宽度内由1过渡到0。分别求解各相内的流体物性，在相界面处借助水平集函数对物性参数进行定义，使物性在界面处得以连续。计算过程见式(8)、式(9)。

δ表示规整后的delta函数，其值见式(10)。

δ为一维Dirac函数，通过δ可将边界上的积分改写为体积上的积分。δ只在相界面附近有值，从而使表面张力的作用限定在界面处。

1.3 边界条件

二维矩形计算域的左侧设置为速度入口，如式(11)。

右侧边界为压力出口，与大气连通，如式(12)。

上下边界定义为周期性流动条件，以减小计算域的面积，缩短计算时间，减少占用内存。域内空气初速度大小等于其入口速度ud,0，雾滴在初始化位置具有竖直向下的初速度，如式(13)。

本文中所有算例雾滴的初速度设定为10m/s[34]。

1.4 网格划分

为提高网格规整度，防止液滴初始化导致的网格不规则分布，使用网格控制边对计算域分区，在各子域内划分网格，如图1中虚线所示。雾滴初始化位置四周使用自由三角形网格，余下区域应用较规整的四边形单元map网格，以提高计算精度并减少所需计算资源。本模拟中通过设置网格最大单元尺寸和最小单元尺寸，在COMSOL 中自动进行网格的生成，计算前进行网格无关性检验，以优化网格的划分。

1.5 后处理

通过雾滴在空间运行中的轨迹变化来判断网格选择的合理性，需要得到不同时刻雾滴的位置。如图2所示，将计算得到的体积分数图像用MATLAB进行编程处理，将图像转化为二值图像，提取相界面的轮廓。通过计算轮廓各点坐标的均值可获得各个时刻雾滴的质心坐标，从而连成雾滴的运动轨迹。

雾滴在空间运行中会由球状变成椭球状，为了描述雾滴在横向风场中下落时的形态演化规律，定义雾滴的形变比E=b/a，其中a 代表椭圆长轴的长度，b代表椭圆短轴的长度。定义雾滴的形变量为e=1-E，当雾滴呈现球状时e=0，椭球状时0＜e＜1。

图2 后处理过程

2 结果与讨论

2.1 网格无关性验证

以直径为100µm 的雾滴为例，设定最小及最大单元尺寸自动划分网格。雾滴下落是一个瞬态过程，以雾滴在0.1ms内的运动轨迹及速度变化作为评价指标来验证网格无关性，结果如图3所示。

当网格数达到45569时，再提高网格精度对雾滴的轨迹及速度影响不大，可认为该网格精度足够高，具有较好的计算性价比。此时网格最大单元尺寸为3µm，最小为0.01µm。所有算例的域内最大速度基本一致，对不同尺寸的雾滴选取相同的网格划分精度，保证了所有雾滴都具有相同的界面厚度。

2.2 雾滴空间形变过程

具有向下初速度的雾滴在左侧来流横向风中运动的形态演化过程，如图4所示。雾滴在一个振荡周期Δtd内经历了“最初的球→椭球→球→椭球→球”的形变过程，椭球长轴方向与竖直方向有一定的夹角。雾滴由起始时刻的球形变为椭球状再回到球形的这一过程简称为“初次变形”，所经历的时间记为Δt1；继而再由球形变为椭球到振荡周期结束这一过程简称为“二次变形”，所经历的时间记为Δt2，Δtd=Δt1+Δt2。雾滴受到左侧来流横向风的作用，同时自身具有向下的初动能，在雾滴下缘产生一个低压区，从左侧来流的横向风更容易进入雾滴下侧，气流对雾滴产生剪切作用，使雾滴形变为长轴偏向上风向的椭球形。随着雾滴不断下降，雾滴的形状逐渐变为具有较大曲率的椭球形，导致雾滴的表面张力增大，逐渐增大的表面张力与曳力相平衡使得形变程度慢慢减弱，在一段时间后几乎又恢复到球状，但由于雾滴的形变仍在继续，二次形变中雾滴的拉伸与压缩方向与第一次相反，曳力的阻碍作用导致形变程度相对于第一次发生了衰减。随着雾滴继续下落，再一次恢复到近似球状，至此一个振荡周期结束。

图3 网格无关性验证(ua=3m/s;μd=1.0mPa·s;σ=72.7mN/m)

图4 150µm的雾滴在横向风中运动的形变过程

雾滴在横向风场中运动的形变过程是多种因素共同作用下的结果。如图5所示，为不同条件下雾滴达到最大形变量时的形态特征。下面依次讨论雾滴直径、横向风速、药液的理化性质等对雾滴形变特性的影响。

图5 不同条件下雾滴达到最大形变量的计算结果图

2.2.1 雾滴直径对雾滴形变的影响

根据美国农业工程师学会（ASAE）对农药雾滴尺寸的分类标准，取直径为100µm、150µm、200µm、250µm、300µm 的细雾滴、中等雾滴、粗雾滴，将它们置于速度为3m/s的横向风场中，液相黏度为1.0mPa·s，表面张力为72.7mN/m。考察雾滴直径对空间形变特性的影响，结果如图6所示。

图6 不同直径雾滴在横向风场中形变量e随时间的变化

结果表明，雾滴的振荡周期Δtd和形变量e 均随着直径的增大而增大。当直径增大时，雾滴的重力和曳力均会增加，在与黏性力及表面张力的竞争中取得优势，保持原有运动状态的趋势增强，雾滴形变的周期会延长，同时雾滴自身重力的提高会增大形变程度。在一个振荡周期Δtd内，形变量的第二个峰值小于第一个峰值emax1＞emax2，同时有Δt1＞Δt2，表明振荡形变发生了衰减。一方面是由于雾滴在下落过程中速度减慢，在一定程度上抑制了形变过程；另一方面是由于气流方向向右，空气对雾滴的拖拽作用总是沿着左上右下这一方向，使得雾滴在此方向上被拉长得更多。另外可以发现，雾滴从椭球变为球的时间略短于从球变为椭球的时间，这是因为对于椭球而言，由于它发生了形变，在风的投影方向长度拉长，空气给予椭球雾滴的曳力要大于球形雾滴，所以在雾滴呈现椭球形时速度衰减更快，形变也就更快。

2.2.2 横向风速对雾滴形变的影响

取直径为100µm、液相黏度为1.0mPa·s、表面张力为72.7mN/m 的雾滴，分别将其置于横向风速为1m/s、3m/s、5m/s、7m/s的稳定气流场中，考察横向风速对雾滴形变特性的影响。

如图7 所示，为雾滴形变量e 随时间的变化情况。随着横向风速的增加，雾滴的形变程度也呈现增大的趋势。雾滴受到空气的曳力作用而发生形变，曳力大小与气液相对速度呈正相关，风速越大意味着空气拖拽雾滴的能力越强，形变越强烈。尤其当风速达到7m/s 时，形变量有明显提升，第一次形变幅值emax1达到了风速为5m/s时的1.6倍左右，第二次形变幅值emax2达到了5m/s 时的1.2 倍左右。此外，分析雾滴的振荡周期可以发现，随着风速的增大，Δt1、Δt2和Δtd几乎不发生改变，即4 种情况下雾滴的振荡周期都为1.32×10-4s，达到第一个形变峰值的时间都在0.048ms左右。结果表明，雾滴的振荡周期受横向风速的影响不明显，振荡周期为一个本征量，由雾滴自身特性决定。

图7 直径为100µm的雾滴在不同横向风场速度下形变量e随时间的变化

2.2.3 液相黏度对雾滴形变的影响

液体的黏度反映出流体内部的黏性效应，黏度越大的液体表明其内部的凝聚力越大，因此黏度的大小势必影响着雾滴的形变程度。考察雾滴初始直径为150µm，横向风速3m/s，表面张力为72.7mN/m，液相黏度分别为0.3mPa·s、0.5mPa·s、1.0mPa·s、2.0mPa·s时，雾滴的空间形变特性。

如图8所示，随着液相黏度增大，雾滴的形变量减小，但振荡周期Δtd没有明显的变化，说明雾滴的振荡周期与液相黏度几乎无关。当液相黏度为0.3mPa∙s 和0.5mPa∙s 时，Δt1与Δt2几乎相等，各占一个完整周期的一半，同时雾滴的振荡幅度也没有衰减，两次形变量最大值几近相同。然而当液相黏度为1.0mPa∙s 和2.0mPa∙s 时，emax1＞emax2，Δt1＞Δt2，形变程度发生了衰减。分析原因认为，当液相黏度较大时，雾滴的形变主要由两相间的剪切作用决定，故前半周期比后半周期长；而当液相黏度较小时，内部凝聚力小，为影响形变的主导作用，形变的衰减现象不再显著。此外，黏度较小的雾滴在第一个振荡周期结束时难以恢复到近球形，黏度越小这一趋势越明显，这是由于雾滴自身的黏度较小，内聚力较弱，在运动过程中难以与曳力相抗衡，很难恢复到最初的形状，只能维持一定的形变量在空气中继续运动。

图8 不同黏度的150µm雾滴在横向风场中形变量e随时间的变化

2.2.4 表面张力对雾滴形变的影响

根据文献[16,35-36]报道，按照一定的配比添加不同种类的助剂可以调控雾滴的表面张力在20～73mN/m 之间，在田间实际操作的范围内。考察雾滴直径为150µm，横向风速为3m/s，液相黏度为1.0mPa·s，表面张力分别为28.5mN/m、37.4mN/m、56.9mN/m、72.7mN/m时空间形变特性。

根据能量最低原理，在气液相界面处的液相表面分子会自发趋向于液相内部而非表面，这使得表面张力对塑造雾滴形态起着十分关键的作用[37]。如图9所示，不同表面张力下雾滴的形变规律与前述算例相同，形变幅度逐渐衰减（emax1＞emax2），形变时间逐渐缩短（Δt1＞Δt2）。当液相表面张力较小时，在与曳力、重力的竞争中处于弱势地位，雾滴恢复到表面能较低的球形状态更加困难，从而延长了雾滴振荡形变周期，增加了形变幅度。无论是前半个周期结束时还是整个振荡周期结束时，表面张力越小的雾滴越难以恢复到初始的球形，如表面张力为28.5mN/m的雾滴在0.392ms时完成整个振荡周期，此时形变量仍有0.111，呈现出长轴与水平方向的椭球状，这一点与低黏度的影响结果较为相似。

图9 不同表面张力的150µm雾滴在横向风场中形变量e随时间的变化

2.3 雾滴形变预测关联式

本文中影响雾滴空间变形特征的因素有雾滴直径、环境风速、雾滴黏度和表面张力，并将雾滴的空间形变归结于重力、曳力、黏性力和表面张力的共同作用。而Re 表征惯性力与黏性力之比，其中惯性力反映出雾滴速度对空间形变过程的影响，同时雾滴与空气的相对速度又影响着曳力，因此曳力对形变过程的影响可通过惯性力进行描述；We 表征惯性力与表面张力之比；Oh 表征黏性力与表面张力之比。结合本文中影响雾滴空间形变的物理参数以及上述量纲为1 参数的意义，选择了韦伯数We、奥内佐格数Oh和雷诺数Re的关联式对形变量和形变周期作出预测。计算公式分别见式(14)～式(16)。

其中，当量直径De选取为雾滴的初始直径，ud和ua分别为雾滴的初始速度和横向空气速度。经过计算后，拟合所涉及的量纲为1 数范围依次为：138.874＜We＜573.466， 66.730＜Re＜207.950，2.873×10-3＜Oh＜1.916×10-2。

采用多元非线性回归方法对计算结果进行拟合，结果如式(17)所示。

其中确定系数R2= 0.9497，如图10 所示对比分析了回归模型的结果与实际模拟的结果，雾滴形变量预测模型的平均相对偏差为8.40%。雾滴的最大形变量与液相黏度和表面张力呈负相关，说明液相黏性力和表面张力能够阻碍雾滴形变，使其维持原有形态，而最大形变比与相对速度和初始直径呈正相关，反映出曳力与重力对形变的增强作用。拟合分析结果与上文讨论结果基本一致。

图10 雾滴形变量模拟结果与回归模型的偏差比较

同时，对初次形变时间占一次振荡形变周期的比例（Δt1/Δtd） 进行拟合，得到如下结果，见式(18)。

其中确定系数R2=0.9150，如图11 所示，对比分析了回归模型的结果与实际模拟的结果，该预测模型的平均相对偏差为5.83%。当雾滴特征参数及环境风速已知时，可以通过式(17)和式(18)分别预测出雾滴的最大形变量和形变周期，以此明晰雾滴的空间形态演化规律，为提出雾滴群的减飘手段奠定理论基础。

图11 雾滴初次形变时间占比模拟结果与回归模型偏差比较

3 结论

通过数值模拟手段对多因素耦合下单雾滴的空间形态演化过程进行研究，结论如下。

（1）雾滴在下落过程中的形变存在振荡周期，且形变会发生衰减。雾滴的尺寸严重影响着振荡周期和幅度，直径越大，雾滴的重力和曳力越大，形变周期越长，形变幅度越大；横向风速几乎不改变形变周期，但风速增大会加剧雾滴形变程度；雾滴黏度主要影响形变程度，黏度越小，形变程度越大，同时在低黏度下雾滴完成一个形变周期后难以恢复到初始的球形，只能保持椭球状继续运动；表面张力值增大导致形变周期缩短，形变程度减弱，当表面张力很小时雾滴也难以恢复球形。

（2）雾滴的形变量受雾滴直径、横向风速、药液黏度、表面张力等因素的共同影响，即形变量是一个由重力、曳力、黏性力及表面张力共同作用的结果；曳力和重力促使雾滴发生形变，而黏性力和表面张力阻碍雾滴发生形变。

（3）通过回归分析得到的雾滴形变量与形变周期的量纲为1预测模型，可以计算出雾滴的空间形变量及形变周期，可明晰雾滴的空间形态演化规律，为提出农药的减飘手段奠定理论基础。

符号说明

A—— 雾滴在风向的投影长度，m

a—— 椭圆长轴，m

b—— 椭圆短轴，m

Cd—— 曳力系数

De—— 雾滴的当量直径，m

E—— 形变比

e—— 形变量

emax1—— 第一个振荡周期内，初次变形过程的最大形变量

Emax2—— 第一个振荡周期内，二次变形过程的最大形变量

Fd—— 单位体积雾滴所受曳力，N/m3

Fσ—— 单位体积流体的表面张力，N/m3

g—— 重力加速度，m/s2

n—— 相界面法向向量

Oh—— 奥内佐格数

p—— 单位体积压力，N/m3

patm—— 大气压，Pa

pout—— 出口压力，Pa

R2—— 拟合相关系数

Re—— 雷诺数

u—— 流体的速度，m/s

ua—— 空气的速度，m/s

ua,0—— 空气的初始速度，m/s

ud—— 雾滴的速度，m/s

ud,0—— 雾滴的初始速度，m/s

uin—— 入口速度，m/s

We—— 韦伯数

x—— 到相界面距离，m

Δt1—— 第一个振荡周期内，初次变形时间，s

Δt2—— 第一个振荡周期内，二次变形时间，s

Δtd—— 第一个振荡周期，s

δ—— Delta函数

κ—— 相界面曲率

μ—— 流体黏度，Pa∙s

μa—— 空气的黏度，Pa∙s

μd—— 雾滴的黏度，Pa∙s

ρ—— 流体的密度，kg/m3

ρa—— 空气的密度，kg/m3

ρd—— 雾滴的密度，kg/m3

σ—— 表面张力，N/m

φ—— 水平集函数

下角标

a—— 空气相

d—— 雾滴相

max—— 最大值