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.
- multiphase flow
- large-eddy simulation
- Lagrangian stochastic model
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 . 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)  has been adapted by Aguirre  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  and in the presence of a gentle sloping hill . 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  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.
2. 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  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.
2.1.1 Initial conditions of the droplet positions and velocities
The initial conditions of ejection droplets are as follows:
The nozzle height was located at h = 0.75 m over the ground.
The elliptical shape A (Figure 1) for exit droplets. It was located at Lc = 0.023 m below the nozzle because, according to the measurements of Nuyttens , it is the region of detachment of droplets from the sheet of the liquid film.
The angle of the spray in the transverse direction is α = 110° according to the technical specifications of the HARDI™ nozzle.
The minor semiaxis of A ellipse (transverse direction y) is dy0 = 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  W0 = 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 (xc, yc) 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 (xc, yc), an algorithm of random variable is used:
where χ is a continuous uniform random variable in the [0, 1] interval whose average value is and standard deviation .
All particles are located at the z = h − Lc height within the A ellipse at the initial time of the simulation.
2.1.2 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 two-parameter 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  presents laboratory 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  measurement experiences are as follows: average droplet diameters , standard deviation of droplet diameters and a modal value . 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  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.
2.1.3 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 and standard deviation can be obtained
where n is the sample size. Michelot  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 and .
With this method, it is possible to simulate diameters of liquid particles that follow the normal distribution from random number generation
where 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.
2.1.4 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 Fli is the force actuating over l liquid particle in i direction i = 1, 2, 3 (x, y, z) on Cartesian coordinate system (Figure 1), ml is the mass of liquid particle, Ui is the air velocity, Vi 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  as:
where CD 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 CD 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 :
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 . From Eqs. (8) and (9):
Note that is the drag coefficient at sedimentation regime. It depends on (Eq. (11)) and therefore on the sedimentation velocity itself (Eq. (10)). So, Vs 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.
2.2 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  is implemented in ARPS by Aguirre . 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.
2.2.1 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 Ui has a deterministic term and a random term:
The tensors and 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), denotes the random characteristic variable of zero mean and covariance:
2.2.2 The Eulerian flow model
It is necessary to simulate the air velocity at liquid particle position Ui. The LES method decomposes in a resolved component
The LES code advanced regional prediction system (ARPS) developed by Center of Analysis and Prediction of Storm (CAPS) and Oklahoma University  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 is a buoyancy force, 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 and are modeled using the dynamic Smagorinsky formulation .
2.2.3 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
2.2.4 Eulerian to Lagrangian coupling
Aguirre and Brizuela  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. :
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 C0 = 2.1. Therefore, we only need to evaluate the αij tensor. Aguirre and Brizuela  propose the expression of αij for inhomogeneous and anisotropic turbulence:
The subgrid turbulent kinetic energy is solved by 1.5 order transport equation  and is the Reynolds SGS stress tensor.
The unresolved velocity component 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 :
With Eqs. (16) and (22–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 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 time step has elapsed. For this reason, the simulation time step must be chosen less than the relaxation time of the smallest possible liquid particle. According to the experimental measurements by Nuyttens , 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).
2.3 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.
2.3.1 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)  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, and 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 as shown in Figure 2.
In Eq. (30), 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 on the plane perpendicular to the relative velocity . The impact factor will be equal to the cosine of γ angle:
So, inserting Eq. (34) into Eq. (33) results in:
It is necessary to obtain Dp. It is the distance from large droplet center to the perpendicular plane at velocity passing through the small droplet center . This plane P-P is shown in green color in Figure 2. The plane equation passing through the small drop center is where are the components of . In this expression, D is a constant of plane equation. This constant is obtained as: . So, the Dp distance can be obtained as:
2.3.2 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  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”.
2.3.3 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 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 (Impc-s) is according to Rabe  as follows:
The transition impact factor between coalescence and reflexive separation (Impc-r) is:
The transition impact factor between reflexive and stretching separation appears when the Wes > 2.5 and can be considered a constant value Impr–s = 0.28.
It should be noted that the boundary curve between coalescence and reflexive separation Impc-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 Impc-s frontier curve, which decreases the coalescence area as the Imp increases.
For bounce, the model proposed by Estrade  calculates the number of transition Weber Web according to the Imp, Δ and a shape parameter, φ:
where ξ is computed as:
and λ = (1 − Imp)(1 + ∆). The shape parameter φ can be computed as Zhang :
The transition bounces into Wes-Imp map collision droplets, and the Weber symmetric bounce frontier Wesb is used. So, it is obtained from Web (Eq. (40)) as:
The Wes-Imp map collision droplets define areas
2.3.4 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  (bounce outcome).
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 (Eq. (31)). For droplets of the same density, its diameter is .
Munnannur and Reitz  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 (Esurten), stretching (Estrtch), and viscous dissipation energies (Edissip) by using a separation coefficient (CVS):
Case CVS ≤ 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  proposes:
where is the coalescence efficiency, , and is the Weber number, which follow O’Rourke model .
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 CVS > 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 ro 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 η = ro. In this model, a time scale of temporal evolution ligament is proposed: . If T ≤ 2, the ligament contracts in a single satellite whose radius is ro. k2 = 0.45 and . Otherwise, it is necessary to compute the evolution time of ligament radius equation for obtaining the final value rbu:
k1 = 11.5 following Kim . Eq. (48) can be solved by iteration with an initial value and ∆t = 1 × 10−2.
The diameter of satellite droplets can be determined by following Georjon :
The number of satellite droplets is calculated from the mass conservation by assuming uniform satellites size . 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:
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  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 , Kim  affirms that no uniform droplet sizes are obtained for the end-droplets 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 . If Nsat ≤ 0, it is assumed that the ligament breaks up without satellite droplet and the two end-droplets have their own radius. If 0 < Nsat ≤ 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 Nsat > 1, the ligament breaks up into uniform droplets with identical diameters . The velocities of end-droplets and satellite droplets are computed with Eqs. (46) and (50).
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 en,p = 0.97 by following Almohammed . 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 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 is the normal impulse of a droplet-droplet collision: .
3. 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 R2 = 0.9965.
3.2 Effects of the Eulerian-Lagrangian double-way coupling
The trajectories of liquid particles are simulated with an Euler-Lagrangian double-way 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 . 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.
3.3 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, 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.
3.4 Simulation of droplet dispersion from a nozzle in a cultivated field
Figure 13 shows the drift simulation of droplets spraying over cultivate field with a nozzle at 0.75 m above the ground.
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  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.
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 R2 = 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.
The present work is funded through the PIO CONICET-UNER 2015–2016 project and the FCyT-UADER research project “Development of a simulation model for the study of the drift of agricultural sprayings, using a flat fan nozzle, from trailing 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.
A. Appendix A
The Weibull probability density function (Weibull p.d.f.) can be used for describing a lot of technical applications for which the distribution of ground material, particles dispersion, or droplet diameters in spray jet normally in the μm-band have behaviors with a random characteristic. In this case, the diameters of droplets ejected from a spray nozzle are simulated using a Weibull p.d.f.
Let us name the random variable that represents the initial diameter of the liquid particles, and the expression of Weibull p.d.f. for this random variable is:
In this p.d.f., m > 0 is the shape parameter and k > 0 is the scale parameter. The form of the density function of the Weibull distribution changes drastically with the value of m parameters. The k parameter does not change the shape of the distribution, but it extends along the random variable . In this way, if the parameters m and k are chosen correctly, it is possible to obtain the shape and stretch of the Weibull p.d.f. that fits the experimentally measured diameter data.
To take these data into account, the cumulative function of the Weibull distribution, named Rosin-Rammler (R-R), is used.
The R-R distribution function is expressed as:
Eq. (A.2) can be written as:
If we associate Eq. (A.3) with a linear equation f(X) = mX + a where the dependent variable is and the independent variable is , m is the slope and a = −mLn(k) is intercept, in which the line cuts on the Y = f(X) axis. The method to obtain the Weibull p.d.f. parameters consists of plotting over logarithmic scale; the results of the experimental measurements vs. and approximate the point cloud to a linear regression by the least squares to obtain the slope (shape parameter) m and the constant term a. The k scale parameter is obtained by .
With these two parameters (m, k), it is possible to obtain the Weibull p.d.f. that describes the droplet diameter’s distribution corresponding to the experiment measurements.