An Eulerian-Lagrangian Coupled Model for Droplets Dispersion from Nozzle Spray

In this chapter, an Euler-Lagrangian double-way coupled model is presented for simulating the liquid particle dispersion ejected from a high-pressure nozzle. The Eulerian code is advanced regional prediction system (ARPS), developed by Center of Analysis and Prediction of Storm (CAPS) and Oklahoma University, USA, which is specialized in weather simulation. This code is the double way coupled with a Lagrangian one-particle model. The theoretical remarks of the double-way coupling, the simulation of the liquid droplet trajectory, and, finally, the droplet collision in the spray cloud using a binary collision model are descripts. The results of droplet velocities and diameters are compared with experimental laboratory measurements. Finally, agrochemical spraying over a cultivated field in weak wind and high air temperature conditions is showed.


Introduction
Numerous engineering applications are focused on solving the problems of dispersion of a sprayed droplet jet from a high-pressure nozzle into a gaseous medium. Spraying is used in internal combustion engines, application of agrochemicals over cultivated fields and greenhouses, irrigation systems, among others. Some questions raised by this topic of engineering can be studied using the computational simulation of multiphase flows. There are currently different ways to implement these tools. In the Eulerian approach, the physical domain is subdivided into cells of a grid space. Each cell has a portion of its volume filled by the liquid phase and the other part by the air. Continuity, momentum, energy and species, for a single-fluid mixture conservation equations, are solved in all pass time of the simulation [1,2]. On the other hand, the Lagrangian approach at one particle proposes the velocities and positions of particles simulation (solids, liquid, vapor, or scalar species) by solving a stochastic equation following the Markov chains. The deterministic term is obtained from the average air velocity values, while the random term is like a white noise following a Brownian motion. The coupled Eulerian large-eddy simulation (LES) with Lagrangian one-particle stochastic method (STO) has been proposed in order to obtain more details on the turbulent properties of the fluid carrying the particles. Several studies using this coupling methodology (LES-STO) can be found [3][4][5][6][7][8][9][10][11]. In this chapter, we focus on the ejection of droplets in air environment from a spray nozzle. The sprayed liquid is a water at 20°C temperature, and the ejection pressure reaches 3 bar. The atmosphere temperature is like the water ejected, but the air pressure is 1.013 bar. These conditions are like as Nuyttens experience [12]. The author carries out paired measurements of droplet diameters and velocities at 25 cm below the spray nozzle using phase Doppler particle analyzer (PDPA) instrument. The particle's Euler-Lagrangian double-way coupling code LES-STO is proposed for to simulate the trajectory of these particles in their liquid phase. The original finite-difference Eulerian LES code named advanced regional prediction systems (ARPS) developed by the University of Oklahoma's Center for Analysis and Forecasting of Storms (CAPS) [13] has been adapted by Aguirre [14] for the simulation of fluid particles in order to validate it with measurements of concentration of a passive gas made in a wind tunnel over flat ground [15] and in the presence of a gentle sloping hill [16]. First time, we present a random ejection algorithm of droplet diameters whose probability density function replies to the two-parameter Weibull distribution. These parameters are previously obtained using laboratory experimental data. Second time, we present the theoretical approach for obtaining the results of collision droplets into the spray. The binary collision droplet model [17][18][19] has been performed in the LES-STO code. This model uses the concept of symmetric weber number [20] to consider the relationship between the kinetic and surface energy of the two colliding droplets. Finally, an agrochemical spraying over a cultivated field in low wind velocity and high air temperature conditions is showed.

Theoretical framework and techniques of numerical simulation 2.1 Conditions and simulation of liquid particle ejection
In this section, we present a random ejection algorithm for simulating different diameters of droplets whose probability density function matches a Weibull distribution. The scale and shape parameters of Weibull distribution are previously obtained from laboratory experimental data using a phase Doppler particle analyzer (PDPA) performed by Nuyttens [12] from an HARDI™ spray nozzle. The sprayed liquid in laboratory experience has been water at 20°C temperature, and the ejection pressure reaches 3 bar. The atmosphere temperature is like the water ejected, but the air pressure is 1.013 bar and a calm wind.
• The minor semiaxis of A ellipse (transverse direction y) is dy 0 = 0.0046 m.
• The major semiaxis of A ellipse (longitudinal direction x) is • The initial vertical velocity of liquid particles is adopted from the confined fluid simulation [21] W 0 = 16.356 m s À1 .
• The initial horizontal component velocity of each liquid particle will depend on the initial position within the A ellipse ( Figure 1): where (x c , y c ) are the relative horizontal positions of the liquid particle at center of the A ellipse. So, the horizontal velocity of droplets at initial time of the simulation depends directly on these relative distances. To determine the initial horizontal position (x c , y c ), an algorithm of random variable is used: where χ is a continuous uniform random variable in the [0, 1] interval whose average value is μ χ ¼ 0:5 and standard deviation σ χ ¼ ffiffi ffi 3 p =6. All particles are located at the z = h À L c height within the A ellipse at the initial time of the simulation.

Initial distribution function of the liquid particles' diameters
The Rosin-Rammler (R-R) distribution function is a cumulative function of continuous random variable whose probability density function (p.d.f.) is a twoparameter Weibull. This distribution function is used [22,23] to adjust experimental data of droplet diameter measurements as a function of liquid-sprayed fraction volume in order to obtain the shape m and scale k Weibull parameters. The experimental data require very precise measurements of the diameters of liquid droplets. The Doppler phase particle analyzer (PDPA) meets the necessary requirements and has the advantage of obtaining paired velocity and diameter data of very small droplets, which are ejected from the nozzle. Nuyttens [12] presents laboratory Spray cone HARDI ™ ISO F 110-O3 nozzle follows Nuyttens laboratory experience [12]. measurements with this device for different working pressures with an HARDI ™ ISO F 110-O3 nozzle, among others, in calm air conditions. The values obtained in the Nuyttens [12] measurement experiences are as follows: average droplet diameters ϕ 0 ¼ 267:6 μm, standard deviation of droplet diameters σ ϕ ¼ 110:3 μm, and a modal value Mo ¼ 250:2 μm. These values allow finding the parameters of the Weibull p.d.f. for the initial conditions of ejected droplet diameter. In this work, we obtained the m and k parameters following the methodology presented [24] for the droplet diameters that are ejected with an internal pressure of 3 bar. The properties of the Weibull distribution and the fit with the experimental measurement data are described in Appendix A.

Initial randomization of the droplet diameters
Once the shape and scale parameters of the Weibull p.d.f. are obtained, which characterize the diameters of drops ejected from the spray nozzle, it is necessary to carry out a temporal sequence for the simulation of these diameters. An algorithm based on the function of the random variable χ, already used, is proposed for the initial position of droplets in the A ellipse ( Figure 1). The randomization algorithm should allow the diameters to be assigned to each ejected droplets such that the mean and standard deviation of p.d.f. droplets simulated over a long period time are close to those corresponding to the Weibull p.d.f. Using the central limit theorem for a set of values corresponding to the random variable χ, a normalized random variable Z of mean value Z ¼ 0 and standard deviation σ Z ¼ 1 can be obtained where n is the sample size. Michelot [25] performed several tests with different n values using 1 million particles to find an acceptable number from computational cost time. The author concludes that n = 50 is a good value to obtain the standard normal random variable for generating random numbers with χ whose μ χ ¼ 0:5 and σ χ ¼ ffiffi ffi 3 p =6. With this method, it is possible to simulate diameters of liquid particles that follow the normal distribution from random number generation where ϕ 0 and σ ϕ are the mean and standard deviation values of the Weibull p.d.f., respectively, whose expressions are given as: However, in Eq. (4), we do not consider the asymmetry of the Weibull p.d.f. For this, it is necessary to incorporate the mode (Mo) of Weibull p.d.f. whose expression is given as: Incorporating this value into Eq. (4), the random variable ϕ can be written as: Eqs. (1), (2), and (7) provide the initial conditions of velocities, positions, and diameters of the ejected liquid particles from the HARDI™ ISO F110-O3 nozzle.

Dynamic parameters of liquid particles
Several simplifications are imposed at ejection and trajectory simulation of liquid particle phenomena: • Particles are considered to have a constant spherical shape in their trajectory.
• The rotating motion of the particle is not considered.
• The ratio between droplet and air densities is very large.
Assuming these simplifications, the force per unit mass to which the liquid particles are submitted is based on a balance between gravity and drag forces per unit mass: where F li is the force actuating over l liquid particle in i direction i = 1, 2, 3 (x, y, z) on Cartesian coordinate system (Figure 1), m l is the mass of liquid particle, U i is the air velocity, V i is the liquid particle velocity, and τ is the characteristic time response or relaxation time of liquid particle, which represents the time required for the liquid adapt to sudden changes in air velocity. This last parameter can be estimated [26] as: where C D is the dynamic drag coefficient due to the air viscosity, ρ l is the liquid particle density, and ρ is the air density. It is important to note that, if the spray injection pressure is increased, the relative velocity between the droplets and the air will be higher. This implies that the relaxation time will be decreased. In addition, if the diameter of the liquid particle decreases, the relaxation time will also be shorter. As the droplet is considered spherical, the drag coefficient C D depends on both, the diameter of droplet and the air viscosity, which is used for the calculation of the Reynolds number referred to the droplet ℜe l : Several drag coefficient expressions have been analyzed [27]. The authors showed that Turton expression [28] has given better results in this simulation case: When the drag and gravity forces are balanced, the liquid particles reach the sedimentation regime. In this case of free fall, the droplets have only vertical velocity component. This velocity is named sedimentation velocity (8) and (9): Note that C D, s is the drag coefficient at sedimentation regime. It depends on ℜe l (Eq. (11)) and therefore on the sedimentation velocity itself (Eq. (10)). So, V s can only be calculated iteratively.
The time elapsed until the particle reaches the sedimentation velocity can be written as: This is an important parameter of liquid particles because if the time elapsed until the liquid particle reaches the ground is longer than the sedimentation time, it will be exposed to drift.

Euler-Lagrangian double-way coupled model
The double-way coupled model presents a bidirectional coupling between the Eulerian and Lagrangian equation systems. Based on the Eulerian approach, the large-eddy simulation (LES) technique is proposed to obtain a detailed turbulent flow. The turbulent intensity of the fluid that transports the liquid particles is taken into account in the simulation of its trajectories. In this approach, it is not possible to obtain a full description of all eddies, so the LES technique is applied for resolving the larges scales of turbulence. The small scales are modeled by subgrid eddy viscosity model (SGS). A dynamic SGS model proposed by Germano [29] is implemented in ARPS by Aguirre [14]. On the other hand, the Lagrangian form is proposed to simulate the trajectories of the liquid particles. In the double-way-coupled LES-STO model, it is considered that the intensity of turbulent flow is taken into account in Lagrangian stochastic equation, and the presence of the liquid particles is taken into account in the momentum equation for LES.

The Lagrangian stochastic model
The governing equations of the liquid particles trajectories are based on a Lagrangian stochastic model at a one-particle and one-time scale following the classical equation of Langevin. The air velocity model at liquid particle position U i has a deterministic term and a random term: The tensors h ij U i ; t ð Þand q ij U i ; t ð Þare determined dynamically according to the characteristics of the turbulence at each position of the simulation and at each instant time. This requires that the LES equations are coupling with Lagrangian stochastic model. In Eq. (14), η t ð Þ denotes the random characteristic variable of zero mean and covariance:

The Eulerian flow model
It is necessary to simulate the air velocity at liquid particle position U i . The LES method decomposes in a resolved component and a fluctuation : The LES code advanced regional prediction system (ARPS) developed by Center of Analysis and Prediction of Storm (CAPS) and Oklahoma University [13] numerically integrates the time-dependent equations of mass balance, forces and energy of the largest turbulent scales. Filtered continuity as in Eq. (17), filtered momentum of fluid velocity as in Eq. (18), and filtered momentum of scalars as in Eq. (19) are described as follows: where B ⊕ is a buoyancy force,S a ij is the anisotropic deformation tensor, ν is the molecular viscosity, and Φ ψ represents the sink and sources of the scalars variables ψ. The variables with tilde indicate that they have been weighted by the density of air , which is only dependent of z (vertical) height. The pressure equation is obtained using the material derivative of the state equation for moist air and replacing the time derivative of density by velocity divergence using the continuity equation. The correlation terms containing unsolved scalesτ ij andτ iψ are modeled using the dynamic Smagorinsky formulation [29].

Lagrangian to the Eulerian coupling
From Lagrangian stochastic equation to the Eulerian LES model taken into account, the number of liquid particles is very large near the nozzle. The coupling has been computed by adding term at filtered momentum of fluid velocity Eq. (18), which expresses the additional momentum due to the presence of liquid particles per volume of carrier fluid: where is the force of the liquid particle l in i direction, ΔV = ΔxΔyΔz is the grid cell volume, and n Δ is the number of liquid particles within the grid cell. Using Eq. (8), the additional momentum is given by:

Eulerian to Lagrangian coupling
Aguirre and Brizuela [11] show that the coupling LES-STO model allows to find the expressions of the deterministic and random terms of Eq. (14) using the velocity-filtered density function (VFDF) proposed by Gicquel et al. [30]: The material derivatives of the velocity-filtered air flow, subgrid turbulent kinetic energy K ¬ , and energy molecular dissipation ε are calculated using the ARPS code at each position and time step of the simulation. The Kolmogorov constant value is C 0 = 2.1. Therefore, we only need to evaluate the α ij tensor. Aguirre and Brizuela [11] propose the expression of α ij for inhomogeneous and anisotropic turbulence: The subgrid turbulent kinetic energy is solved by 1.5 order transport equation [31] and R ij ¼ u À i u À j ⊕ is the Reynolds SGS stress tensor.
The unresolved velocity component u À j in Eq. (16) is obtained in discrete form using the Markov chains: where the subscript (n) denotes the value at present time of the simulation, while (nÀ1) is the value in the previous time step. The first pass time, subscript (0), is considered as isotropic homogeneous turbulence [32]: With Eqs. (16) and (22)(23)(24)(25), it is possible to calculate the air velocity at the liquid particle position. The equations describing the motion of the liquid particle in its discrete form are: > : It is necessary to note that in the first Eq. (26), the rate Δt=τ must be greater than one for convergence of the numerical solution. Cases in which the liquid particle diameter is very small, the second Eq. (26) must be used considering that the sedimentation velocity Vs has been reached before Δt time step has elapsed. For this reason, the simulation time step Δt must be chosen less than the relaxation time of the smallest possible liquid particle. According to the experimental measurements by Nuyttens [12], the liquid particles whose diameters are smaller to 50 μm are exposed to drift before reaching the ground. These particles have relaxation times less than 7.6 ms. So, a simulation time step Δt = 0.2 ms has been chosen due to little size of cell grid ΔV, whereby liquid particles whose diameters are smaller than 7 μm are being considered for calculation with the second Eq. (26).

Binary collision droplet model
Once the droplets ejected from the spray nozzle and having simulated their positions, velocities and diameters along their trajectory, it is necessary to consider the collision. Binary droplet collision models are a widely used theoretical approximation to obtain the outcome of the interaction droplets [17][18][19][20][33][34][35][36][37][38][39]. This model consists of estimating the positions, velocities, and diameters of droplets after the collision. In addition, satellite droplets can be created from the ligament breakup as a consequence of it.

Parameters of binary collision
The binary droplet collision is simulated using three important parameters. The ratio of the droplet diameters Δ (Eq. (27)), the dimensionless symmetric Weber number (Wes) [20] relating kinetic energy vs. surface energy (Eq. (28)), and the dimensionless impact parameter (Imp) takes into account the way in which the two droplets impact (Eq. (29)): where the subscripts S and L indicate the smaller and larger droplet, respectively, σ denote the surface tension, V ! mS and V ! mL are the relative velocities to the mass center of the incoming droplets, and X is the projection of the distance between the droplet centers in the normal direction to the relative velocity Figure 2.
In Eq. (30), V ! mR is the velocity of mass center. If the droplets have the same density: The relative velocity droplets of the mass center can be resumed using Eq. (27): For the impact parameter in Eq. (29), it is necessary to compute X variable. This variable is the projection of segment b ¼ 0:5 ϕ S þ ϕ L ð Þon the plane perpendicular to the relative velocity V ! R . The impact factor will be equal to the cosine of γ angle: So, inserting Eq. (34) into Eq. (33) results in: Figure 2.

Scheme of small droplet S and large droplet L before collision (dashed line) and at contact instant (solid line).
It is necessary to obtain D p . It is the distance from large droplet center x L ; y L ; z L À Á to the perpendicular plane at V ! R velocity passing through the small droplet center x S ; y S ; z S À Á . This plane P-P is shown in green color in Figure 2.

The plane equation passing through the small drop center is
In this expression, D is a constant of plane equation. This constant is obtained as: D ¼ Àu R x S À v R y S À w R z R . So, the D p distance can be obtained as:

Numerical resolution of the impact coefficient
In each time of numerical simulation, it checks whether the collision between two droplets occurs. For obtaining a more optimize algorithm, collision boxes are placed around and inside the liquid particle ejection spray. The sizes of grid boxes vary dynamically, adjusting to the boundaries of the particle domain as shown in Figure 3a. The size of the boxes is the same as the Eulerian calculation grid in horizontal direction Δb = Δx = Δy = 0.1 m. In this way, every drop inside this box will be questioned about whether it collided with the other drops that are in the same box. When a binary collision is successfully found (e.g., 5-6; 3-7 in Figure 3b), the pairs are marked and removed from the next iteration of detection. This technique was proposed by Michelot [25] and used by Aguirre [4,14] to consider the diffusion of chemical species that are carried by fluid particles. It is evident that due to the temporal discretization of the numerical solution of droplet motion equations, it is almost impossible that for an instant of discretized time t (n) = t (n-1) + Δt, the contact between two drops can be concurrent. Most likely, by that instant time t (n-1), the contact is about to occur and at instant time t (n), it has already occurred. When the distance between the centers of the two drops inside a collision box is less than the sum of their radii, then the drops collided. In this case, the time Δt" elapsed from the collision to the computation instant time t (n) to be calculated. The particles are repositioned at the moment of collision t" = t (n) À Δt", and the impact factor is calculated according to the positions of their centers as shown in Figure 4. The lapse time Δt" for repositioning the droplets at instant of collision is computed in resolving: The outcomes of collision droplets are computed using the map collision theory. Once the droplets collided and the effects of collision are into account on the droplets, they are repositioned by advancing the same pass time Δt".

Binary droplet collision map
The outcomes of the binary droplet collision model propose different scenarios: Coalescence: the two droplets that collide to form a single drop as a result of the collision. In this case, the surface energy is relatively greater than the kinetic energy.
Reflexive: the two colliding droplets almost head-on, so they join together as one, but the kinetic energy is large enough to separate again and can generate satellite droplets.
Stretching: the two drops collide tangentially, so they separate and can generate satellite droplets.
Bouncing: the two colliding drops remain separated after collision without exchanging mass between them. Figure 5 shows a time sequence of the binary droplet collision for each outcome described above. It is important to note that the result of the binary collision depends Droplets before collision at t (nÀ1) instant time.
not only on the velocity of both drops but also on their relative size and impact coefficient. Two of the four possible outcomes of the binary collision are susceptible to generating satellite droplets. These droplets are usually much smaller in size than the parent-drops and are, therefore, more prone to drift and evaporation. If these droplets are composed of a phosphonate-acid solution (such as glyphosate), then after evaporation, the solute will drift away from the airflow very quickly.
The outcomes of collision droplets are defined using a map collision. This map is the graphic representations between the Wes vs. Imp (Wes-Imp) frontier curves among the different outcomes of binary collision that are displayed on this map. Several researchers proposed equations for frontier curves. The transition impact factor between coalescence and stretching separation (Imp c-s ) is according to Rabe [20] as follows: Imp c-s ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0:53 2 þ 4:24Wes p À 0:53 4Wes : (38) The transition impact factor between coalescence and reflexive separation (Imp c-r ) is: The transition impact factor between reflexive and stretching separation appears when the Wes > 2.5 and can be considered a constant value Imp r-s = 0.28.
It should be noted that the boundary curve between coalescence and reflexive separation Imp c-r increases with the increase of Wes to the value of Imp = 0.28. This behavior indicates that for low Imp values (on-head collision) and relatively low droplet velocities before collision, surface energy is greater than kinetic energy and the result of the collision is stable coalescence. However, for the same Imp values but with higher velocities, the kinetic energy is predominant; the droplets have an unstable coalescence and then separate. This separation can generate satellite droplets. On the other hand, if the Imp is higher (tangential collision of the droplets), then coalescence as a result of the collision is more improbable since only a fraction of the volume of the drops interacts during the collision. The contact surface of both drops is smaller and therefore the surface energy as well. This reduces the likelihood of stable coalescence as a consequence of the collision. This behavior is evident in the Imp c-s frontier curve, which decreases the coalescence area as the Imp increases.
For bounce, the model proposed by Estrade [35] calculates the number of transition Weber We b according to the Imp, Δ and a shape parameter, φ: where ξ is computed as: if λ > 1, and λ = (1 À Imp)(1 + Δ). The shape parameter φ can be computed as Zhang [19]: The transition bounces into Wes-Imp map collision droplets, and the Weber symmetric bounce frontier Wes b is used. So, it is obtained from We b (Eq. (40)) as: The Wes-Imp map collision droplets define areas , , and where the outcomes of the binary droplet collision are represented. These areas are bounded by frontiers curves as proposed in Eqs. (38)-(43). The areas with frontier curves are shown in Figure 6. The frontier curves between and change their position as a function of Δ.

Numerical models of the binary droplet collision
The binary droplet collision model allows obtaining the diameters and velocities of the droplets after the collision. The values of these variables are obtained according to the proposed models [17][18][19]34] (coalescence, reflexive, and stretching separations) and [35] (bounce outcome).

Coalescence
For coalescence outcome, the two droplets coalesce into one. This occurs preferably at low Weber numbers as surface tensions exceed kinetic energy. The new droplet velocity is the velocity of mass center before the collision V (31)). For droplets of the same density, its diameter is ϕ new

Stretching
Munnannur and Reitz [17] calculate the interaction volume between the droplets. This volume is released from both drops creating a ligament that gives rise (or not) to satellite droplets. This volume is computed and taken into account the magnitude of the opposing surface (E surten ), stretching (E strtch ), and viscous dissipation energies (E dissip ) by using a separation coefficient (C VS ): Case C VS ≤ 0: if this coefficient is carried at negative value, it is assumed that fragmentation of droplets does not occur. The drops only lose kinetic energy. The center mass of droplet velocities is affected by Z coefficient indicating the fraction of energy that is dissipated during collision. For this case, Kim [18] proposes: where e coal ¼ min 1:0; 2:4 f We À1 Â Ã is the coalescence efficiency, , is the Weber number, which follow O'Rourke model [40]. The relative velocities of mass center after collision can be written by using momentum conservation equation: The velocities after collision can be obtained by using Eqs. (30) and (31). The diameters of droplets after collision are unaltered.
Case C VS > 0: the separation volumes from droplets determine the evolution of the temporary fluid ligament that would form between them. In this model, it is assumed that the ligament has a uniform cylindrical shape, and the radius r o of ligament at initial instant time of the temporal evolution with a momentum balance equation can be obtained: where Ψ S and Ψ L are the fraction of volumes lost from the smaller and large droplets to form the ligament [17,18], and η is its initial time instantaneous length ( Figure 7). Another assumption is η = r o . In this model, a time scale of temporal evolution ligament is proposed: T ¼ 0:75k 2 ffiffiffiffiffiffiffiffiffiffi We 0 p . If T ≤ 2, the ligament contracts in a single satellite whose radius is r o . k 2 = 0.45 and We 0 ¼ 2r 0 it is necessary to compute the evolution time of ligament radius equation for obtaining the final value r bu : k 1 = 11.5 following Kim [18]. Eq. (48) can be solved by iteration with an initial value r bu =r 0 ð Þ¼1 and Δt = 1 Â 10 À2 . The diameter of satellite droplets can be determined by following Georjon [41]: The number of satellite droplets is calculated from the mass conservation by assuming uniform satellites size N sat ¼ 6 r 0 =ϕ sat ð Þ 3 . The velocities of satellite droplets can be obtained from momentum equation where the velocities of parent droplets after collision are computed by Eq. (46). where the diameters of parent droplets after collision are:

Reflexive
The volume of ligament is the entire temporarily merged mass of two droplets. The model of satellite droplet formation is similar at stretching outcome, but the initial radius of ligaments is The model proposed by Munnannur [17] for reflexive outcome uses the time scale of temporal evolution ligament. When T ≤ 3, a single satellite droplet is formed and the three droplets (considering the ligament breaks up into two end-droplets and one-satellite droplet) have same size. However, according to the experimental studies of Ashrgiz [34], Kim [18] affirms that no uniform droplet sizes are obtained for the enddroplets and a single satellite droplet after reflexive collision. We have adopted the last criteria, so the satellite droplet diameters are computed with Eqs. (48) and (49), but the number of satellite droplets is N sat ¼ 6 r 0 =ϕ sat ð Þ 3 À 2. If N sat ≤ 0, it is assumed that the ligament breaks up without satellite droplet and the two enddroplets have their own radius. If 0 < N sat ≤ 1, it is assumed that a single satellite droplet is formed that is smaller than the two end-droplets after collision. The diameter of single satellite droplets is ϕ sat and the end-droplets after collision have identical diameters: When N sat > 1, the ligament breaks up into uniform droplets with identical diameters ϕ L new

Bounce
In this case, the droplets bounce maintaining their diameters after the impact. In the general case, oblique collision between droplets is considered. The droplet velocities after collision must be decomposed into a normal component and a tangential component to the plane of impact. The tangential component after impact remains unchanged, but the normal component is affected by a soft inelastic rebound assuming a restitution coefficient e n,p = 0.97 by following Almohammed [42]. This restitution coefficient takes into account the dissipation of kinetic energy during the impact. The normal velocities of droplets at instant of collision are as follows: where x i R ¼ x i S À x i L is the relative position between the droplets in i = 1, 2, 3 (x, y, z) and b is the center droplet distance (Figure 2). The normal component velocities after collision are given as: where f dd n is the normal impulse of a droplet-droplet collision:

Results and discussion
3.1 Ejection of liquid particle simulation In order to obtain the Weibull p.d.f. corresponding to the droplet diameters as described in Section 2.1, the scale k = 301.228 and shape m = 2.606 parameters were obtained with R-R method. Figure 8 shows the minimum square adjusted of regression line with a correlation coefficient R 2 = 0.9965.

Effects of the Eulerian-Lagrangian double-way coupling
The trajectories of liquid particles are simulated with an Euler-Lagrangian doubleway coupled model descript in Section 2.2. The influence of droplets to air velocity is shown in Figure 9 (a) for t = 1 s instant time and (b) for t = 20 s instant time of the simulation. The vertical velocities of air W are shown in color scale. Eddies around the spray plume are formed and extend up to 3 m from the center of the spray (c). It is observed that eddies are formed by influence of jet droplets. This effect should be taken into account as droplets of very small diameters are captured by eddies and are prone to contribute to drift. This effect can be seen in Figure 10 where the small droplets follow the streamlines of eddies on both sides of the sprayer.
The results of the vertical droplet velocities distribution as a function of the droplet diameters obtained at 0.35 m below the nozzle are shown in Figure 11. These are compared with the laboratory measurements of Nuyttens [12]. It is  observed that the dispersion of velocity values for each diameter class is greater in laboratory measurements than in simulation. In addition, for diameters less than 200 μm, the model slightly underestimates the vertical velocity values relative to the laboratory results.

Binary collision droplet map
The collision map for binary droplet model descripted in Section 2.3 is shown in Figure 12. The map allows showing the events of coalescence, bounce, reflexive,  Mean and extreme range values measured by Nuyttens [12]. and stretching separation. When considering the total number of droplet binary collision events, 21.1% corresponds to coalescence, 0.6% to reflexive separation, 8.8% to stretching separation, and 69.5% to bounce. The amount of satellite droplets arising from the separation by reflexive and stretching is displayed with numbers. It is noted that the number of satellite drops increases with the number of symmetrical Weber for both separately. This behavior indicates that the greater velocity the droplets are ejected from the spray nozzle, the more likely it is that satellite droplets will appear as a result of reflexive and stretching separation. As mentioned above,  this can cause an increase in the proportion of sprayed product not reaching its destination, leaving it adrift. Figure 13 shows the drift simulation of droplets spraying over cultivate field with a nozzle at 0.75 m above the ground.

Simulation of droplet dispersion from a nozzle in a cultivated field
The meteorological conditions of air temperature at nozzle level are 30°C with 2 m s À1 velocity wind. The simulation time shown in Figure 13 is 20 s after the start of spraying. The drift of small droplets (less than 50 μm in diameter) exceeds 8 m in the area of application. Of the total liquid sprayed, 0.43% corresponds to droplets smaller than 50 μm measured by Nuyttens [12] in wind tunnel at 50 cm below the spray nozzle. In the simulation shown in Figure 13, this percentage does not change because the satellite droplets generated are greater than 80 μm. The number of satellite droplets generated by stretching and reflective separation in these conditions was obtained. Of the 120 satellite droplets analyzed, 35.3% have diameters less than 150 μm, 61.3% have diameters between 150 and 250 μm, and 3.4% have diameters between 250 and 350 μm. There were no satellite droplets with diameters larger than 350 μm.

Conclusions
In the present work, it was possible to simulate and validate the ejection velocity of the liquid particles from an HARDI™ ISO F110 03 nozzle placed at 0.75 m over ground. The diameters of the drops were randomized to the volume applied following a procedure of Rosin-Rammler distribution function for obtaining the parameters of Weibull probability density function with a correlation coefficient R 2 = 0.997. The double-way coupled Euler-Lagrangian model has been used for obtaining the trajectory of droplet spraying. Eddies at both sides of spraying have been captured by the model. These extend up to 3 m from the center of the spray. The vertical component droplet velocity was simulated and validated with laboratory measurements. The velocity of droplets smaller than 200 μm slightly underestimates with respect to laboratory data. The binary collision models have been implemented into the code to consider particle collision events. A collision detection algorithm using collision boxes was presented and used to optimize computation times. The drift of droplets with air temperature of 30°C and wind speed value 2 m s À1 has been simulated in cultivated fields. The drift of droplets smaller than 50 μm diameter exceeds 8 m of the application area. No satellite droplets smaller than 80 μm are generated under field simulation conditions. The largest proportion of satellite droplets generated as a result of the droplet collision has a diameter between 80 and 250 μm. equipment." We are also grateful to the Laboratory of Prototyping of Electronics and 3D Printing of the School of Engineering, UNER, for the work of printing the nozzle.