农药液滴撞击移动液膜和非对称冠状水花形成机理分析

2021-06-09蒋彦龙

农业机械学报 2021年5期

赵可王瑜蒋彦龙

(1.南京航空航天大学飞行器环境控制与生命保障工业和信息化部重点实验室，南京 210016；2.南京工业大学城市建设学院，南京 210009)

0 引言

目前，喷雾防治是我国治理农业病虫害的主要手段。在喷洒过程中，通过喷头将农药输送到靶标植物上，该过程不仅要关注有害生物数量的变化，还应评估农药对周围环境的影响[1]。研究发现，喷洒的农药仅有20%～30%沉积在靶标植物表面[2]，大量农药流入土壤及水环境，从而引发邻近作物药害事故。研究农药在靶标植物上的沉积流失规律对获得农药最佳使用效益、降低农药负效应具有重要意义[3]。

农药在植物叶片上的最终沉积分布取决于药液雾化参数(粒径、运行速度)、药液的物理或化学特性(黏度、表面张力、密度)[4]、喷洒系统[5]、操作参数、气象环境、叶片表面结构、冠层结构等[6]。已有研究者对农药的雾化特性及液滴分布进行了研究，结果表明：农药的喷雾效率和有效性与喷嘴型号、喷量流速、喷雾压力密切相关，增大压力、风速容易引发农药脱靶飘移[7-12]。

农药的喷洒过程可归结为多个液滴冲击植物叶片的复杂物理场耦合，属于典型的多相流问题。因此，明晰单个液滴对植物叶片的冲击行为是研究农药沉积机理的基础。现有研究多集中于农药液滴的空间分布，缺乏对液滴冲击植物叶片的动力学行为研究。另外，受喷洒工况和农药种类的影响，药液在植物叶片上的沉积结构随之发生改变。研究发现，药液动态表面张力、粘度、密度对药滴反弹、铺展、飞溅扩散具有影响，静态表面张力决定药滴的铺展和最终形成的沉积结构[13]。在农药喷洒过程中，液滴碰触植物叶片并迅速铺展成液膜，随后液滴持续冲击液膜的过程属于流体力学中的液滴撞击液膜过程。

目前，撞击液膜研究主要借助高速摄像机和液滴发生器进行相关试验，或利用计算机技术对撞击过程进行数值分析。由于成像技术和微尺度测量技术的局限性，相关参数无法准确控制，如液滴撞击参数(直径、撞击速度)、液膜流动参数(液膜厚度、液膜流动速度、液膜表面微观结构)等，而且更难以获得撞击过程中涉及的流体内部微观数据。因此，数值方法已逐渐成为研究流体运动过程的重要手段。数值研究无需试验成本，不受试验过程中不确定性因素的影响，并且可以准确捕捉液膜内的压力场和速度场等信息。常用的数值模拟方法包括CLSVOF法(Coupled level set and volume of fluid)、格子-玻尔兹曼方法(Lattice-Boltzmann method, LBM)、光滑粒子动力学方法(Smoothed particle hydrodynamics method, SPH)、界面追踪法(The front tracking method，FTM)等。许多学者采用数值方法建立了液滴撞击液膜的二维物理模型，并取得令人满意的结果[14-21]。

在农药实际喷洒过程中，多采用飘移式喷雾。通常植物叶片呈一定倾角，这导致在植物叶片上形成的药液膜以不同速度沿径向向四周运动。目前，关于液滴撞击移动液膜的研究非常有限，由于液滴与液膜之间存在剪切应力，碰撞的动力学行为变得更加复杂，产生不同于冲击静态液膜的独特动态特征，因此针对液滴在静止液膜上的撞击研究具有更广泛的应用和物理意义。一些学者研究发现，液滴撞击运动液膜后会发生合并、形成喷射、喷射破碎、冠状飞溅和冠状喷射等现象，并认为这和液滴撞击速度与液膜移动速度之比、液膜厚度、液膜表面的结构(平面膜、毛细波和波峰)以及液体的物理性质有关，而且流动液膜的冲击还会引发壁面上的热流分布不均匀[22-26]。

可以看出，相关研究主要集中于冲击静态液膜的试验观测与理论分析上，针对流动液膜的研究非常有限。本文利用CLSVOF方法模拟单液滴撞击移动液膜的运动过程，结合冲击后的形貌演变和液体内部流场状态，揭示非对称水花的形成机理，并对二次液滴的产生和细射流的形成机制进行验证，深入研究农药在靶标植物上动态沉积过程，分析农药在植物叶片上的沉积机理，为探究农药的最终有效沉积和流失规律提供理论支撑。

1 计算模型

CLSVOF方法包括Level Set函数和VOF函数初始化、流动控制方程求解、Level Set函数和VOF相函数(φ函数和F函数)的对流输运方程求解、相界面重构以及函数再次初始化等。在VOF方法中，F函数表示计算单元中液相占据单元空间的体积分数，Level Set方法中用φ函数代数值来区分计算域中的各相，CLSVOF方法可有效解决VOF和Level Set两种气液两相流动相界面追踪方法存在的计算收敛性、稳定性及准确性不高等问题[27-28]。

CLSVOF方法通过φ函数和F函数共同构造相界面，采用分段线性界面重构思想，在法向方向上移动界面使得单元液态区域面积比率和F函数值相匹配，单元中心到相界面的垂直距离由割线法迭代求解得到。

界面法向量n计算公式为

(1)

界面曲率κ(φ)为

(2)

对流输运方程为

(3)

(4)

式中u——流动速度，m/s

t——时间，s

φ函数经过对流输运方程求解后，将不再保持距离函数的性质，因此必须对 φ函数重新初始化，再次重新初始化主要包括φ函数符号以及φ函数值的确定。φ函数值是单元中心到相界面最小距离，而φ函数的符号由F函数值决定[27]

Sφ=sign(0.5-F)

(5)

当F<0.5时，单元中心在气态区域中，φ函数为正；F>0.5时，单元中心在液态区域中，φ函数为负；当F=0.5时，表示相界面，φ=0；sign()是符号函数。

流动控制方程为

(6)

(7)

其中

τ=2μ(φ)S

式中ρ——密度，kg/m3

p——压力，Pa

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

τ——黏性应力张量

Fs——表面张力源项

表面张力采用CSF(Continuum surface force)模型

(8)

式中σ——表面张力系数，N/m

H(φ)——Heaviside函数

引入Heaviside函数来光顺界面处的密度和黏度[28]

H(φ)=

(9)

式中h——网格尺寸，mm

μ——黏度，Pa·s

应变率张量S表达式为

(10)

则不同区域密度ρ和黏度μ可由Heaviside函数计算[28]，即

ρ(φ)=ρl(1-H(φ))+ρgH(φ)

(11)

μ(φ)=μl(1-H(φ))+μgH(φ)

(12)

式中ρl——液相密度，kg/m3

ρg——气相密度，kg/m3

μl——液相黏度，Pa·s

μg——气相黏度，Pa·s

初始时刻如图1所示，本文在xoz平面建立二维计算模型，初始时刻圆形截面底部与液膜截面相切，初始直径D0、初始速度U0的液滴垂直向下撞击厚度为hf、速度为uf的液膜，流动方向向右，液膜和液滴均以水作为流动介质，x表示铺展前端距离中心轴的长度，“-”、“+”分别代表运动方向，以对称中心为原点，左侧铺展长度记为负，右侧记为正，则铺展直径Dt=x+-x-。考虑重力的影响，顶部及两侧为压力出口边界，底部为无滑移壁面，将气液两相均视为不可压缩流体。

计算相关参数见表1，其中液滴撞击参数(直径、撞击速度)和液膜流动参数(液膜厚度、液膜流速)参照文献[29]试验数据，本文着重考察液膜流动的影响，涉及的液滴初始直径和速度均为恒定，为了便于对比分析，将相关计算参数无量纲化，具体公式为

表1 相关参数

式中φ——无量纲铺展因子

X——无量纲长度

T——无量纲时间

U——无量纲速度

H——无量纲厚度

We——液滴韦伯数

利用Fluent 17.0软件求解数值模型，控制方程采用有限体积法进行离散，压力与速度耦合采用Coupled方法，压力求解采用PRESTO方法，采用QUICK算法对Level Set方程进行求解，采用Geometric-Reconstruction Scheme 求解VOF方程，动量和能量方程采用二阶迎风格式，时间步长Δt=10-6s，Δt内迭代次数为20。

2 单液滴撞击壁面运动特性分析

2.1 计算模型验证

二维网格划分包括三角形单元网格和四边形单元网格，四边形均匀网格在离散域逼近，计算精度以及网格质量等方面均优于三角形单元网格，故本文采用四边形网格对计算域进行划分，为了对网格进行无关性验证，模型网格数分别划分为80 000、180 000、320 000、500 000，选择铺展因子作为比较参数，图2给出了不同网格数的无量纲铺展因子随无量纲时间的变化曲线，对比分析4种网格的无量纲铺展因子随无量纲时间的变化曲线可见，当网格数为320 000时计算结果与网格数不再相关，因此在后续计算中网格数均采用320 000。

另外，本文以文献[29]试验结果为参照，在相同试验条件下(无量纲厚度H=0.05，无量纲速度U=0.31)对模拟结果进行定量对比，对比结果见图2，试验采用高速摄像机记录撞击过程(帧率为5 000 f/s,快门速度30 μs)，通过注射泵产生小液滴，液滴直径D0由注射泵下端毛细管孔径控制，调节生成液滴高度来控制液滴撞击速度U0，利用加压罐垂直向下供应液体射流并在基底表面形成径向流动液膜，液膜的流动参数(液膜厚度hf和流速uf)根据WATSON理论[30]计算(该理论计算预测结果已被AZUMA等[31-32]采用多普勒相位仪验证)，计算式为

hf=

(13)

(14)

式中χ——撞击点到射流中心距离，mm

a——射流喷嘴的孔半径，mm

V——射流体积流率，cm3/s

由图2可以看出铺展因子曲线的模拟结果与试验结果线较为接近，说明计算结果与试验结果一致性良好。

图3a为相同试验条件下单液滴撞击静止液膜的模拟结果与试验结果[29]的定性对比(无量纲厚度H=0.26，无量纲速度U=0)。t=0.4 ms时，液滴撞击静止液膜后在液滴周围形成小水花，随着撞击过程持续，水花不断外延，并产生竖直向上运动的薄膜射流，t=6.0 ms时，形成完美对称冠状水花，在冠状水花顶部边缘形成指状结构，并在顶部边缘发生飞溅，产生二次小液滴，模拟结果与试验结果吻合较好。图3b为相同试验条件下单液滴撞击水平右移液膜模拟结果与试验结果[29]的定性比较，其中无量纲厚度H=0.058、无量纲速度U=0.479。不同于液滴在静止液膜上的冲击，在移动液膜上观察到非对称的冠状水花的独特特征，这种非对称行为在xoz平面上尤为显著；t=1.0 ms时，在液滴周围观察到不对称的细射流并逐渐演变成非对称的冠状水花，且伴有二次小液滴；在5.0 ms时，撞击上游细射流仍斜向上移动，而在撞击下游射流在流动液膜中逐渐淹没，并在铺展边缘作波状运动，与试验结果吻合。事实上，试验结果和模拟结果存在一定差距，原因可能是试验中液膜流动参数和液滴撞击参数的控制难度较大，流动液膜表面不可能为绝对平面且生成的液滴达不到完美球形，这与模型中液膜和液滴的几何参数存在一定差距。尽管如此，由图2和图3的定量及定性对比可以看出，本文建立的二维模型仍可以捕捉到相关现象的基本特征，如在xoz平面上，撞击上下游射流的不对称、二次液滴的产生等主要特征，说明本文建立的液滴撞击移动液膜的二维数值模型是可行的。

2.2 液膜流动速度的影响

农药喷洒工况的调整、大气环境的改变以及靶标植物叶片的倾斜都会引起液膜流动，为了深入研究液膜的流动对撞击行为的影响，图4给出了液滴撞击不同流速液膜的形貌演变(红色虚线为z轴)，其中无量纲速度U为0、0.4、0.8，无量纲厚度H=0.2，液滴韦伯数We=333，液膜流动方向水平向右。可以看出，液膜的流动对撞击行为影响显著，液膜流速越大，水花的非对称特征越明显，在撞击初期(0.2 ms)，液滴惯性力驱使液滴迅速冲击液膜，在液滴与液膜交界处产生小的细射流，形成冠状水花的冠壁，并伴随飞溅的二次液滴，液滴底部不断挤压周围液膜形成了冠基；随着撞击过程的持续(0.5～1.0 ms)，冠壁高度和冠基直径逐渐增大；提高液膜流速后，两股撞击流体在冠基中心的左侧聚集，撞击上游的冠壁呈现水平运动的趋势，而撞击下游的冠壁下部与液膜几乎垂直，非对称行为愈加显著；而且上下游射流端部的飞溅现象也不一致，液膜流速越大，撞击上游的飞溅现象越明显，而撞击下游则与之相反。

农药液滴在植物叶片上的沉积实质上是液滴在液膜上的铺展，图5为不同移动速度下无量纲铺展长度随无量纲时间的变化曲线，无量纲厚度H=0.2，液滴韦伯数We=333。由图5可知，液膜流动方向与撞击上游液滴的铺展方向相反，阻碍液滴的铺展，但与撞击下游的铺展方向相同，故液膜流速固定时撞击下游的无量纲铺展长度(X+)总大于撞击上游(X-)，液膜流速与上游的无量纲铺展长度成反比，与下游成正比；进一步分析发现，冠基中心随着液膜的流动逐渐向右推移，冠基中心与z轴不再重合，且冠基中心右移幅度随着流动速度增大而增大，当无量纲速度U=0.8，无量纲时间T=3.47时，非对称冠状水花完全迁移至z轴右侧，这是因为本文模拟工况中0

2.3 非对称冠状水花形成机理

由上述模拟结果可以看出，液滴撞击移动液膜产生的非对称水花主要体现在左右两侧射流发展行为不一致、冠基液体厚度不均匀以及射流端部的飞溅现象不对称。对于农药喷雾防治而言，这种非对称特征容易引发靶标植物叶片上面农药分布不均匀，尤其是飞溅行为将直接导致植物叶片上农药的流失。目前，关于液滴撞击移动液膜的试验或理论研究相对较少，更缺乏非对称水花演变机理研究。为了对液滴撞击移动液膜产生水花演变过程进行深入分析，图6(图中压差单位为Pa)分别给出了液滴撞击静止和移动液膜(无量纲速度U=0.8)的压差分布图，其中红线为相界面，黑线为等压线，无量纲厚度H=0.2。梁刚涛等[27]研究认为液滴撞击后的颈部位置较大的压差导致细射流的发生。观察图6a发现，在0.05 ms时，液滴冲击静止液膜后在两侧颈部出现压力较大的滞止点，越靠近滞止点压力越大，而在两相界面处压力接近大气压，在颈部与相界面产生13 700 Pa的压差，出现颈部凸起；较大的压差足以克服表面张力的束缚，当撞击持续至0.2 ms时，颈部凸起演变为高速运动的细射流，此时压差快速下降至4 250 Pa。

对于移动液膜，液体内的压力分布与静止液膜截然不同，由图6b可以看出，0.05 ms时流体内部高压区集中在撞击上游颈部，压差为25 000 Pa，远高于撞击下游，且约为静止液膜的1.8倍。如图6b所示，在较大压差的驱动下，聚集在高压区的大量液体向外部气流的低压区流动，在撞击上游颈部位置形成斜向上微小凸起。相反的是，在撞击下游颈部位置附近的压力梯度较小，液体凸起不明显。随着撞击过程的持续，撞击颈部位置的压力迅速下降，凸起不断增大演变成射流，撞击上下游射流的演化不再对称，说明射流产生的压差机制同样适应于流动液膜。

由图4和图6发现，液滴撞击静止液膜时，两侧射流顶部末端射流呈现水平运动的趋势，梁刚涛等[27]将其归因于液膜内部流体的径向运动作用。而当液滴撞击移动液膜时，随着液膜流速的增大，撞击上游的射流水平运动更显著，而撞击下游射流更倾向于垂直向上延伸，似乎径向运动不再主导下游射流的延伸方向。为了进一步探究上述现象，图7给出了0.2 ms时在不同液膜流速下的速度分布和涡量云图，其中红线为相界面，无量纲厚度H=0.2。

YARIN等[33]提出液滴撞击静止液膜产生的冠状水花是由于流体内部的运动间断(Kinematic discontinuity)导致的；由图7a可知，静止液膜流体内速度分布呈现较好的对称性，对于不可压缩的液体，液滴撞击液膜后导致流体向下运动，流体碰撞到固壁后出现流动滞止区[34](右侧蓝色虚线圈)，随后液体的流动方向发生改变，形成径向运动，径向运动的液滴流冲击周围静止液膜时导致运动间断，出现速度不连续现象(左侧蓝色虚线圈)，并沿着速度不连续前沿产生流出的细射流，射流速度达到8.6 m/s。

对于流动液膜而言，液滴撞击移动液膜时流体内速度分布会出现非对称特征。由图7b可知，当无量纲速度U=0.4时，撞击上游液膜流动方向与液滴的铺展方向相反，液滴的铺展受阻，运动间断现象非常明显，两股流体在运动间断处汇聚，汇聚的流体因壁面的法向应力(向上蓝色箭头)垂直向上运动[34-35]，由于流动液膜和液滴的铺展相互挤压，导致汇聚的流体与液滴接触区域增大，液滴内部径向运动对射流延伸方向的影响被加剧，因而表现出上游射流沿水平方向运动较明显；但在撞击下游液滴铺展方向与液膜流动方向相同，在液膜流动剪切力的作用下，运动间断被明显削弱，汇聚的流体沿着液膜流动方向迅速迁移，此时液滴与汇聚的流体与液滴接触区域减小，径向运动效应对下游的射流影响随之减弱，故表现在图4中撞击下游的射流更倾向于向上运动。另外，运动间断处汇聚的流体既是质量又是动量的汇聚[23]，显然上游汇聚流体的质量和动量高于下游，这意味着有更多的液体进入上游射流，且上游射流末端速度达到10.1 m/s，高于下游射流末端速度8.1 m/s。分析图7c可知，提高液膜流速至无量纲速度U=0.8时，液膜水平流动惯性显著增大，运动间断区逐渐沿着液膜流动方向迁移，最终在撞击上游运动间断区和滞止区发生合并(蓝色虚线圈)，大量液体在此处聚集，这与图4中冠基左侧液膜隆起，而右侧液膜凹陷是一致的；液膜流速的增加导致汇聚的液体质量和动量均迅速增加，因而更多的流体进入射流，两侧射流速度分别增长至10.6 m/s和8.5 m/s；图7中速度分布总体上表现出撞击上游射流末端速度、静止液膜射流末端速度、撞击下游射流末端速度依次增大的趋势。

从涡量的角度分析，如图7所示，当液滴撞击静止液膜时，在射流端部上下两侧出现一对顺时针(上方)和逆时针(下方)的反向旋涡。可以看出，随着液膜的流动和射流的发展，在撞击上游射流的下方出现逆时针旋涡比上方的顺时针旋涡更剧烈，并且随着液膜流速的增大，在下游射流下方的旋涡逐渐消失，而上方的旋涡逐渐增强，LIU等[36]研究认为旋涡方向和大小对射流发展有较大的影响；在撞击上游较强逆时针旋涡驱使射流向液膜方向弯曲[37]，即有利于射流沿着水平方向运动，而撞击下游仅在射流上方出现逆时针旋涡有利于射流端部向上运动；总体看来上游射流向水平方向运动显著，下游射流更倾向于垂直向上发展，但在下游射流末端仍有水平运动的趋势；可见在移动液膜上，射流发展受径向运动和涡量共同作用，但径向运动对射流发展趋势影响较大。结合分析涡量云图和速度分布可知，增大液膜流速加剧射流向液膜流动相反方向延伸的趋势，这对药液在植物叶片上有效沉积及沉积均匀性是不利的。因此在实际农药喷洒过程中，需尽可能保持农作物叶片表面液膜处于静止或低速流动状态，减少农药的射流损失。

飞溅行为直接影响药液在靶标植物叶片上的沉积量，本文研究发现飞溅行为在撞击上、下游表现不一致，戴剑锋等[19]在研究液滴撞击倾斜液膜时也发现两侧射流端部飞溅现象不同，并将其归因于被撞击液膜流动方向的改变，而SHETABIVASH等[38]提出层间流体的剪切作用导致流体中发生界面不稳定是射流端部发生飞溅的关键，THORODDSEN等[39]发现剪切失稳源于铺展的速度与液滴径向流速度的速度差，二者速度差越大越容易发生飞溅；这说明射流端部发生飞溅与液膜的流速密切相关；对液膜施加水平运动速度后，撞击上游的液滴径向流动方向与液膜流动方向相反，二者速度差较大，沿着射流端部分离出较多的二次小液滴，而在撞击下游二者同向运动导致速度差减小，撞击下游的飞溅迟缓，故水花两侧的飞溅现象不一致。对于农药喷洒而言，撞击下游的飞溅迟缓意味着药液脱靶率低，这对于提高农药喷洒效率是有利的。