Numerical Analysis of Spreading Process of Ellipsoidal Spraying Droplet Impacting on Superhydrophobic Surface

Agricultural spray deposition is especially important for pesticide application because low efficiency can lead to environmental pollution, poor biological efficiency and economic loss. The deposition of pesticide spray on the leave surfaces is related to the impact kinetic behavior of droplets. But after considering the deformation of the droplet, how impingement will affect the deposition is an interesting research. In this study, a superhydrophobic surface was used to replace the plant surface that the pesticide droplets may affect. An interface tracking method was proposed to characterize the impingement dynamics behaviors of different ellipsoid droplets impacting on the surface. The maximum spreading coefficient and time of ellipsoidal droplets increased with the raise of their size. A lower sized droplet has a faster spreading rate, while the center of a higher sized droplet is thinner. As the velocity of pesticide increases, maximum spreading coefficient of droplet increases with a decrease in the maximum spreading time of droplet. The simulation results can contribute to provide theoretical basis for improving spray efficiency.


Introduction
Spraying pesticide is one of the indispensable methods to prevent and control pests in agricultural production. The pesticide preparation is usually diluted with water and sprayed onto the leaves of the target plant after atomizing with plant protection machinery. But many spray programs widely used in the agricultural industries seem to provide poor pest control compared to experimental results (Zwertvaegher, et al., 2014). Statistics results show that hundreds of millions tons of liquid pesticide are applied to farmland every year in China. However, most pesticide droplets bounce or float away and enter the environment, which not only reduces the operation efficiency, but also deteriorate the ecological environment. Therefore, it is essential to study the potential mechanism of droplet impact behavior to improve pesticide utilization rate.
At present, a large number of scholars have carried out related researches on the impact of liquid droplets on the plant surface. From the view of fluid mechanics, the impact kinetic energy and surface energy of pesticide droplets are converted into viscous dissipation energy and surface energy during the collision. Subsequently, the droplet oscillates with gradually attenuated amplitude until it reaches a sessile state (Kim, et al., 2001;Subhasish, et al., 2013). There are a lot of interesting things happening when droplets hit the leaves of plants. For example, when the droplets hit a barley leaf, the bounce of the droplet is the main dynamic behavior of the impact, but in Wenzel's model, there are also spatter phenomena (Boukhalfa, et al., 2014). Dorr et al. (2014) studied that droplets hit the surface of different plant leaves, they would also break, bounce and stick according to different impact conditions. Andrade et al. (2012) experimentally studied the surface of cabbage and banana skin impacted by different droplets, and found that viscosity had a greater impact on the dynamic characteristics of the droplets, and droplet occurrence: no bounce, partial bounce and all bounce. Zhu et al. (2018) used CLSVOF method to numerically simulate the process of liquid droplets hitting tea leaves and proved that this method could be well applied in agricultural research.
The phenomenon of droplets impacting on hydrophobic or superhydrophobic surfaces is more complicated. Bange et al. (2015) numerically investigated bouncing and non-bouncing of droplets during impact on superhydrophobic surfaces. They found that the droplet bounces off the surface because the total energy of the droplet at maximum recoiling exceeds the initial surface and gravitational energy, while the non-bouncing droplet was due to the competition between kinetic energy and surface tension. Wang et al. (2020) found that the contact time of bouncing droplets on high solid fraction surfaces can be reduced by reducing the texture size to ~100 nm. They demonstrated that this phenomenon could be attributed to the advantage of line tension on nano-textures, and that compact arrangement of nano-textures is essential to withstand the impact pressure of raindrops.
When the droplets move towards the surface, it will deform under the influence of actual factors such as airflow and vibration, making the droplet non-standard spherical (Kékesi, et al., 2014). Yun et al. (2013Yun et al. ( , 2014 used an annular electrode to stretch the droplet along the direction of the parallel wall into a long ellipsoid and hit the wall in an asymmetric way. The effects of droplet morphology, impact velocity, wall wettability and other physical parameters on the rebound height were studied in combination with experiments and numerical simulation. Furthermore, they studied the process of egg-shaped droplet impacting (Yun, et al., 2018).
In the process of pesticide spraying, the influence of deformation on its impact on the surface cannot be ignored. At the same time, the spreading process of droplets on the leaves has a great influence on the retention. The study of the spreading characteristics of ellipsoidal droplets after impact on the wall can not only deepen the understanding of the droplet impact characteristics under complex conditions, but also promotes the application of a numerical simulation method in agricultural research to improve the efficiency of pesticide spray.

Numerical Method
In the CLSVOF method, the volume function α is used for interface reconstruction, and the level set function ϕ is used to accurately calculate the curvature and normal vectors on the interface (Sussman, et al., 2000). The continuity equation and momentum equation of two phases are respectively expressed as: where U is the velocity vector, g is the gravity acceleration, P is the pressure, ( ) φ ρ is the density, ( ) φ μ is the viscosity, F is the surface tension.
The normal vector n and the curvature ( ) φ κ can be defined as follows: During calculation, Heaviside function was used to smooth the density and viscosity of the interface, which could be expressed as: . 1 ε , and r Δ is the minimum size of the grid. The density and viscosity at the gas-liquid interface after smoothing are expressed as follows: where g and l represent the gas and liquid phases, respectively.
The continuous surface force (CSF) model proposed was used to solve the surface tension (Brackbill, et al., 1992). It can be described by the following equation: where σ is the surface tension coefficient.

Feasibility Verification of the Numerical Model
When the droplet was not deformed, its diameter was 2.13 mm, the impact velocity was 0.447 m/s, and the static contact angle was taken as 160°. The calculation area was 5 × 5 mm and a 2D axisymmetric model was adopted for numerical calculation, the grid was quadrangular uniform, the pressure velocity coupling was PISO method, the Level Set equation was discretized by the finite volume method, QUICK format was adopted, and the pressure term was PRESTO!, geometric reconstruction method was selected for gas-liquid interface interpolation. Figure 1a is an experimental morphology of the ellipsoid droplet which generated though vibration of the equipment impacting the wall (Paola, et al., 2015) and Figure 1b shows the numerical simulation results. It can be shown that the numerical simulation results are in good agreement with the experimental results. According to the test, the difference between forward and backward contact angles and the static contact angles of 2.25 mm droplet on the superhydrophobic wall surface was not more than 5° (Lin, et al., 2017), so the static contact angles were adopted.

Computational Domain
Under the action of forces, the spherical droplet will be compressed or stretched into ellipsoid. In the following text, the radius of the spherical droplet with the same volume will be taken as the equivalent radius (r e ) of the ellipsoidal droplet. The ellipsoidal droplet has an asymmetry along x and y axes, and the radius is a and b respectively, as shown in Figure 2. To facilitate the study, the ratio of different axial radii is defined as the axial diameter ratio AR = a/b. A 2D axisymmetric model is adopted for calculation, in which the x-axis is the one of symmetry and the size of domain is 10 × 20 mm.
We meshed the calculation area by ICEM and selected 6 mesh sizes to verify the independence of grids . The area of the AR = 0.65 droplet is 2469, 4380, 6835, 9835 and 13379 times of the unit grid area. Comparing the spreading radius at the same time shows that the refinement of the grid has a certain improvement on the simulation results, but when the mesh size reaches (1/9835) AR = 0.65 droplet area, the spreading radius no longer changes. Therefore, the mesh size of (1/9835) AR = 0.65 droplet area is chosen as the computing grid. The method used in the calculation is the same as that in section 2.2.
The initial contact between droplet and wall surface is the initial moment of impact. According to the equal volume, the following formula can be obtained:  Because the droplet size is small and many plants have smooth waxy layer or tomentum on their leaves, the impact surface is simplified to a superhydrophobic flat wall. The air and water are incompressible fluid, and the ambient pressure is maintained at 1 standard atmosphere (101325 Pa). Other relevant parameters are shown in Table 2.  Figure 3 shows the morphological characteristics of 5 different droplets during spreading. By comparing different AR droplets in Figure 3 at the same time, it can be found that each droplet has a similar form at the initial stage of as.ideasspread.org Vol. 2, No. 2; impacting. At the time of 0.2 -0.8 ms, the top of each droplet did not change significantly, because the droplet was mainly affected by surface tension and viscous force so remained intact in its initial shape. By observing the morphology of each droplet at 0.4 ms, it can be seen that liquid at the bottom forms a relatively thin layer with a certain thickness, which extends outwards with the development of time. Due to the small influence of wall wettability at the initial stage of impact, the size and shape of each liquid thin layer were similar. However, by comparing the morphology of droplet maximum spreading, there was less liquid in the center of the droplet with a high AR, and the liquid film formed was relatively thinner. The average thickness of the liquid film in the center area of AR = 1.40 droplet is 0.28 mm, which is about 43% of AR = 0.65 droplet. This means that droplets with high AR are distributed more evenly when acting on the target surface.  Figure 4 shows changes of spreading radius (R x ) of 5 AR droplets with time. The R x of the ellipsoidal droplet increases rapidly, and the spreading speed gradually decreases. After reaching the maximum spreading radius, the spreading stops under the action of surface tension and viscous force (Eggers, et al., 2010). With the increase of AR, the initial value and curve position gradually decrease, which indicates that droplet with high AR spreads slowly in a period of time after the collision. From calculation, the spreading radius of AR = 0.65 droplet at 0.1 ms is 57.4% higher than that of AR = 1.40 droplet, and the maximum difference between both droplets at 0.6 ms is 0.48 mm. By observing maximum spread radius data, it was shown that with increase in AR, the maximum value was larger. This indicates that lower sized droplets being rapidly spread at early stage do not have the maximum spreading radius and reach the maximum spreading time relatively early.

Agricultural Science
A total of 21 ellipsoidal droplets with the same liquid volume (AR in the range of 0.35 ~ 2.02, r e = 1) were further analyzed numerically and found that with increase in AR, the maximum spreading time (T x ) and the maximum spreading coefficient (β x ) of ellipsoid droplet as a whole showed an increasing trend, as displayed in Figure 5. The maximum spreading coefficient of the ellipsoid droplets defined as:

Ellipsoid Droplet Flow Field Analysis
The left side of Figure 6 shows the isoline of pressure inside the 1/2 droplet in unit Pa, while the right side is trace line and velocity isoline with unit of m/s, the black solid line is droplet free surface. By comparing the velocity distribution, there is an obvious velocity gradient inside each droplet. Flow velocity near the phase interface in the droplet is relatively larger, while that near the wall and center is smaller. This flow distribution makes liquid inside droplet constant on both sides of the impact point, forming a thin layer at the bottom. Velocity comparison of 5 AR droplets on the outside of thin layers, shows that the lower AR the higher liquid velocity at the outside of the thin layer.
The internal pressure of droplet forms a distinct gradient, and wall pressure is maximum. The difference is that at 0.4 ms, with increase in AR of ellipsoidal droplets, the area of high-pressure shows a trend of gradual increase. With the pressure distribution of each droplet at 0.4 ms, the internal pressure changes more rapidly at early stage of the impact of the low AR droplet. In conclusion, the higher spreading rate of the low AR droplet at early stage of impact is due to higher pressure on the central point and the rapid change in pressure.

Effect of Droplet Velocity (We) on Droplet Spreading
Weber number (We) is the ratio of inertial force and surface tension. The impact velocity is usually converted to the We when research the effect of the droplet's velocity, so we got five different We by changed the impact velocity of the droplet. Figure 7 shows variations of maximum spreading coefficient (β x ) and time (T x ) depending on the velocity of droplets after hitting the surface wall. The maximum spreading coefficient of different droplets increases gradually with raise in We as shown in Figure 8a. β x regarding droplet size ratio AR = 1.40 increases by 0.44, 0.48, 0.57, 0.63 and 0.77, respectively, compared with that of AR = 0.65. This indicates that droplet velocity has a more significant effect on higher-sized droplets. As droplet velocity increases, the droplet spreading process is significantly affected by inertia forces, with a larger value of maximum spreading coefficient.
From Figure 7b, the maximum spreading time regarding a droplet size ratio AR = 1.40 is always longer under different droplet velocity data, which correspond to larger spreading radius. However, with increase in droplet velocity, the time for AR = 1.40 droplet to reach maximum spreading always decreases. For example, when We = 53.70, its maximum spreading time T x decreases by 25.7% compared with that of We = 6.85. From We = 17.53, the maximum spreading time did not change significantly despite the increase in droplet velocity regarding droplet size ratios like AR = 1.00 and AR = 1.40. This is mainly due to the fact that, in a certain range of We, droplet size ratios less than unity (AR ≤ 1) are significantly affected by surface tension inherent to the spreading process.

Conclusion
To some extent, the maximum spreading coefficient and time of ellipsoid droplets increased with the raise of their size. Lower-sized droplet has a faster spreading rate due to the large pressure and its fast change in droplet central point. The influence of inertia force on liquid in the center of higher-sized droplet makes the maximum spreading coefficient larger. At maximum spreading time, central point of higher-sized droplet is thinner. Within a certain range, maximum spreading coefficient of different droplets increases significantly with increase in droplet velocity. The detailed research on the spreading of ellipsoidal droplets not only provides a foundation for the regulation of droplet size, but also promotes the use of interface tracking method in pesticide spraying research.