Open access peer-reviewed chapter

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

Written By

Carlos G. Sedano, César Augusto Aguirre and Armando B. Brizuela

Submitted: 18 June 2018 Reviewed: 23 August 2018 Published: 18 December 2018

DOI: 10.5772/intechopen.81110

Chapter metrics overview

864 Chapter Downloads

View Full Metrics


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.


  • droplets
  • spray
  • multiphase flow
  • large-eddy simulation
  • Lagrangian stochastic model

1. 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.


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 [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.

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 [12], 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 dx0=tgα/2Lc=0.0328m.

  • The initial vertical velocity of liquid particles is adopted from the confined fluid simulation [21] 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 μχ=0.5 and standard deviation σχ=3/6.

Figure 1.

Spray cone HARDI ™ ISO F 110-O3 nozzle follows Nuyttens laboratory experience [12].

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 [12] 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 [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.

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 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 σχ=3/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.

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 [26] 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 el:


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 UiViδi3=Vs. From Eqs. (8) and (9):


Note that CD,s is the drag coefficient at sedimentation regime. It depends on el (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 [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.

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 hijUit and qijUit 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:


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

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˜ija 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 u˜i=ρ¯ui, 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 τ˜ are modeled using the dynamic Smagorinsky formulation [29].

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

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:


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:


2.2.4 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 C0 = 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 Rij=uiuj is the Reynolds SGS stress tensor.

The unresolved velocity component uj 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–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).

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) [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, VmS and VmL 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 VR=VSVL as shown in Figure 2.


Figure 2.

Scheme of small droplet S and large droplet L before collision (dashed line) and at contact instant (solid line).

In Eq. (30), VmR 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 VR. 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 xLyLzL to the perpendicular plane at VR velocity passing through the small droplet center xSySzS. This plane P-P is shown in green color in Figure 2. The plane equation passing through the small drop center is uRxS+vRyS+wRzS+D=0, where uRvRwR are the components of VR. In this expression, D is a constant of plane equation. This constant is obtained as: D=uRxSvRySwRzR. 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 [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:


Figure 3.

(a) Collision boxes around and inside the spray ejection of droplets and (b) droplets inside the collision box at t(n) instant time of the simulation.

Figure 4.

Positions of droplets before and after the collision for a time lapse Δt. Droplets before collision at t(n−1) instant time. Droplets after collision at t(n) = t(n−1) + Δt instant time. Repositioning of droplets at the instant of collision (t″ = t(n) – Δt″).

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:

  1. 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.

  2. 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.

  3. Stretching: the two drops collide tangentially, so they separate and can generate satellite droplets.

  4. 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.

Figure 5.

Time sequence diagram of the binary droplet collision and its outcomes.

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 [20] 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 [35] 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 [19]:


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

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
change their position as a function of Δ.

Figure 6.

Map collision droplets with areas of outcomes collision. Coalescence, reflexive separation, stretching separation, and bounce. Bounce frontiers: Δ = 1, Δ = 0.75, Δ = 0.5, and Δ = 0.1.

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 [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 Vnew=VmR (Eq. (31)). For droplets of the same density, its diameter is ϕnew=ϕL3+ϕS313. 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 (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 [18] proposes:


where ecoal=min1.02.4fWe1 is the coalescence efficiency, f=Δ32.4Δ2+2.7Δ1, and We=ρlϕSΔ3σ1VR2, 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 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: T=0.75k2We0. If T ≤ 2, the ligament contracts in a single satellite whose radius is ro. k2 = 0.45 and We0=2r0ρl/σVR2. Otherwise, it is necessary to compute the evolution time of ligament radius equation for obtaining the final value rbu:


Figure 7.

Collision model for the stretching outcome. (a) Formation instant time of ligament and (b) temporal evolution ligament.

k1 = 11.5 following Kim [18]. Eq. (48) can be solved by iteration with an initial value rbu/r0=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 Nsat=6r0/ϕsat3. 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:

ϕLnew=1ΨL16πϕL3ϕSnew=1ΨS16πϕS3.E51 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 r03=ϕL/23+ϕS/23. 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 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 Nsat=6r0/ϕsat32. 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 ϕLnew=ϕSnew=ϕsat. The velocities of end-droplets and satellite droplets are computed with Eqs. (46) and (50). 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 en,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 xiR=xiSxiL 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 fndd is the normal impulse of a droplet-droplet collision: fndd=1+en,pVnSVnL.


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.

Figure 8.

Regression line of R-R distribution function and data measurement of droplet diameters.

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.

Figure 9.

Influence of droplets to air velocity at different instant times of the simulation. (a) t = 1 s, (b) t = 20 s, (c) zoom of eddy formed at 2.5 m from center spray at 20 s, and W is a vertical component of air velocity around the spray.

Figure 10.

Position of droplets at t = 20 s of the simulation classified by their diameters.

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.

Figure 11.

Distribution of vertical droplet velocities in (m/s) as a function of the diameters (μm). Droplet simulation. Mean and extreme range values measured by Nuyttens [12].

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.

Figure 12.

Map outcomes from binary droplet collision model. Coalescence, reflexive, stretching, and bounce. The numbers next to the symbols indicate the number of satellite droplets formed.

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.

Figure 13.

Drift of spraying droplets from a nozzle at 0.75 m over 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 [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.


4. 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 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 ϕ0 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 ϕ0. 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 Fϕ0 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 fX=LnLn1Fϕ0 and the independent variable is X=Lnϕ0, 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 LnLn1Fϕ0 vs. Lnϕ0and 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 k=expa/m.

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.


  1. 1. Battistoni M, Xue Q, Som S. Large-eddy simulation (LES) of spray transients: Start and end of injection phenomena. Oil & Gas Science and Technology—Rev IFP Energies Nouvelles. 2015;2015:1-24. DOI: 10.2516/ogst/2015024
  2. 2. Sabelnikov V, Lipatnikov A, Chakraborty N, Nishiki S, Hasegawa T. A transport equation for reaction rate in turbulent flow. Physics of Fluids. 2016;28(8):081701. DOI: 10.1063/1.4960390
  3. 3. Wei G, Vinkovic I, Shao L, Simoëns S. Scalar dispersion by a large-eddy simulation and a Lagrangian stochastic subgrid model. Physics of Fluids. 2006;18:095101. DOI: 10.1063/1.2337329
  4. 4. Aguirre C, Brizuela A, Vinkovic I, Simoëns S. A subgrid Lagrangian stochastic model for turbulent passive and reactive scalar dispersion. International Journal of Heat and Fluid Flow. 2006;27(4):627-635. DOI: 10.1016/j.ijheatfluidflow.2006.02.011
  5. 5. Vinkovic I, Aguirre C, Simoëns S. Large-eddy simulation and Lagrangian stochastic modelling of passive scalar dispersion in a turbulent boundary layers. Journal of Turbulence. 2006;7(30):1468-5248
  6. 6. Vinkovic I, Aguirre C, Simoëns S, Gorokhovski M. Large-eddy simulation of droplet dispersion for inhomogeneous turbulent wall flow. International Journal of Multiphase Flow. 2006;32(3):344-364
  7. 7. Vinkovic I, Aguirre C, Ayrault M, Simoëns S. Large-eddy simulation of the dispersion of solid particles in a turbulent boundary layers. Journal of Boundary Layers-Meteorology. 2006:1472-1573. DOI: 10.1007/s10546-006-9072-6
  8. 8. Stephens DW, Keough S, Sideroff C. Euler-Lagrange large eddy simulation of a square cross-sectioned bubble column. In: APCChE 2015 Congress incorporating Chemeca; 27 Sept–01 Oct 2015; Melbourne, Victoria; Paper no. 3134665. 2015. pp. 1-12
  9. 9. Zamansky R, Vinkovic I, Gorokhovsky M. Acceleration in turbulent channel flow: Universalities in statistic, subgrid stochastic models and an application. Journal of Fluid Mechanics. 2013;721:627-668. DOI: 10.1017/jfm.2013.48
  10. 10. Orcellet E, Berri J, Aguirre C, Müller G. Atmospheric dispersion study of TRS compounds emitted from a pulp mill plant in coastal regions of the Uruguay river, South America. Aerosol and Air Quality Research. 2016;16(6):1473-1482. DOI: 10.4209/aaqr.2015.02.0112
  11. 11. Aguirre C, Brizuela A. Computational tools for the simulation of atmospheric pollution transport during a severe wind event in Argentina. In: Atmospheric Hazards—Case Studies in Modeling, Communication, and Societal Impacts. Chapter 6. InTech Open Science, Open Mind.; 2016. pp. 111-136. DOI: 10.5772/63552
  12. 12. Nuyttens D. Drift from field crops prayers: The influence of spray application technology determined using indirect and direct drift assessment means [thesis]. Germany: Faculteit Bio-ingenieurs wetenschappen, Katholieke Universiteit Leuven; 2007
  13. 13. Xue M, Droegemeier K, Wong V. The advanced regional prediction system (ARPS). A multi-scale nonhydrostatic atmospheric simulation and prediction model. Part I: Model dynamics and verification. Meteorology Atmospheric Physics. 2000;75:161-193
  14. 14. Aguirre C. Dispersión et Mélange Atmosphérique Euléro-Lagrangien de Particules Fluides Réactives. Application à des cas simples et complexes [thesis]. Lyon, France: Université Claude Bernard; 2005. 337 p
  15. 15. Fackrell J, Robins A. Concentration fluctuation and fluxes in plumes from point sources in a turbulent boundary layers. Journal of Fluid Mechanics. 1982;117:1-26
  16. 16. Gong W. A wind tunnel study of turbulent dispersion over two- and three-dimensional gentle Hills from upwind point sources in neutral flow. Boundary Layers Meteorology. 1991;54:211-230
  17. 17. Munnannur A, Reitz R. A new predictive model for fragmenting and non-fragmenting binary droplet collisions. International Journal of Multiphase Flow. 2007;33(8):873-896. DOI: 10.1016/j.ijmultiphaseflow.2007.03.003
  18. 18. Kim S, Lee D, Lee C. Modeling of binary droplet collisions for application to inter-impingement sprays. International Journal of Multiphase Flow. 2009;35:533-549. DOI: 10.1016/j.ijmultiphaseflow.2009.02.010
  19. 19. Zhang H, Li Y, Li J, Liu Q. Study on separation abilities of moisture separators based on droplet collision models. Nuclear Engineering and Design. 2017;325:135-148. DOI: 10.1016/j.nucengdes.2017.09.030
  20. 20. Rabe C, Malet J, Feuillebois F. Experimental investigation of water droplet binary collisions and description of outcomes with a symmetric Weber number. Physics of Fluids. 2010;22:047101. DOI: 10.1063/1.3392768
  21. 21. Sedano C, Aguirre C, Brizuela A. Simulación de la Eyección de Spray Líquido desde un Pico de Pulverizadora para Aplicación de Herbicidas. Asociación Argentina de Mecánica Computacional AMCA - Mecánica Computacional. 2017;VXXIII:1049-1068
  22. 22. Ayres D, Caldas M, Semião V, da Graça Carvalho M. Prediction of the droplet size and velocity joint distribution for sprays. Fuel. 2001;80:383-394
  23. 23. Baetens K. Development and application of drift prediction models in fields praying [thesis]. Leuven: Faculteit Bio-ingenieurs wetenschappen, Katholieke Universiteit; 2009
  24. 24. Macías-García A, Cuerda-Correa E, Díaz-Díez M. Application of the Rossin-Rammler and Gates-Gaudin-Schuhmann models to the particle size analysis of agglomerated cork. Material Characterization. 2004;52:159-164. DOI: 10.1016/j.matchar.2004.04.007
  25. 25. Michelot C. Développement d’un modèle stochastique lagrangien. Application à la dispersion et à la chimie de l’atmosphère [thesis]. Lyon, France: Université Claude Bernard; 1996. 180 p
  26. 26. Holterman H. Kinetic and evaporation of waterdrops in air, Wageningen: IMAG. Report 2003-12; Wageningen UR, InstituutvoorvMilieu en Agritechniek; 2003. 67 p
  27. 27. Sedano C, Aguirre A, Brizuela B. Numerical simulation of spray ejection from nozzle for herbicide application: Comparison of drag coefficient expressions. Computers and Electronics in Agriculture. Elsevier. Forthcoming
  28. 28. Turton R, Levenspiel O. A short note on the drag correlation for sphere. Powder Technology. 1986;47:83-86. DOI: 10.1016/0032-5910(86)80012-2
  29. 29. Germano M, Piomelli U, Moin P, Cabot WH. A dynamic subgrid-scale eddy viscosity model. Physics of Fluids. 1991;A(3):1760-1765
  30. 30. Gicquel LYM, Givi P, Jaberi FA, Pope SB. Velocity filtered density function for large-eddy simulation of turbulent flow. Physics of Fluids. 2002;14(3):1196-1213
  31. 31. Deardorff JW. Stratocumulus-capped mixed layer derived from a three dimensional model. Boundary Layer Meteorology. 1980;18:495-527
  32. 32. Pope SB. Lagrangian PDF methods for turbulent flows. Annual Review of Fluid Mechanics. 1994;26:23-63
  33. 33. Brazier-Smith PR, Jennings SG, Latham J. The interaction of falling water drops: Coalescence. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences. 1972;326(1566):393-408. DOI: 10.1098/rspa.1972.0016
  34. 34. Ashrgiz N, Poo JY. Coalescence and separation in binary collisions of liquid drops. Journal of Fluid Mechanics. 1990;221:183-204
  35. 35. Estrade JP, Carentz H, Lavergne G, Biscos Y. Experimental investigation of dynamic binary collision of ethanol droplets—A model for droplet coalescence and bouncing. International Journal of Heat and Fluid Flow. 1999;20(5):486-491. DOI: 10.1016/S0142-727X(99)00036-3
  36. 36. Brenn G, Valkovska D, Danov K. The formation of satellite droplets by unstable binary drop collisions. Physics of Fluids. 1994;13(9):2463, 2001-2477. DOI: 10.1063/1.1384892
  37. 37. Kollár LE, Farzaneh M, Karev AR. Modeling droplet collision and coalescence in an icing wind tunnel and the influence of these processes on droplet size distribution. International Journal of Multiphase Flow. 2005;31(1):69-92. DOI: 10.1016/j.ijmultiphaseflow.2004.08.007. ISSN 0301-9322
  38. 38. Ko GH, Ryou HS. Modeling of droplet collision-induced breakup process. International Journal of Multiphase Flow. 2005;31:723-738. DOI: 10.1016/j.ijmultiphaseflow.2005.02.004. ISSN 0301-9322
  39. 39. Hu C, Xia S, Li C, Wu G. Three-dimensional numerical investigation and modeling of binary alumina droplet collisions. International Journal of Heat and Mass Transfer. 2017;113:569-588. DOI: 10.1016/j.ijheatmasstransfer.2017.05.094
  40. 40. O’Rourke P. Collective drop effects in vaporizing liquid sprays [PhD dissertation]. Princeton, NJ: Dept. Mech. Aerospace Engg., Princeton University; 1981
  41. 41. Georjon TL, Reitz RD. A drop-shattering collision model for multidimensional spray computations. Atomization and Sprays. 1999;9(3)
  42. 42. Almohammed N. Modelling and simulation of particle agglomeration droplet coalescence and particle-wall adhesion in turbulent multiphase flow [thesis]. Hamburg, Germany: Helmet-Schmidt University; 2018; 402 p

Written By

Carlos G. Sedano, César Augusto Aguirre and Armando B. Brizuela

Submitted: 18 June 2018 Reviewed: 23 August 2018 Published: 18 December 2018