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

Abstract

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.

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 L_c = 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 dy₀ = 0.0046 m.

• The major semiaxis of A ellipse (longitudinal direction x) is ${\mathit{dx}}_0=\mathit{tg}\left(\alpha /2\right)L_c=0.0328\mathrm{m}$.

• The initial vertical velocity of liquid particles is adopted from the confined fluid simulation [21] W₀ = 16.356 m s⁻¹.

• The initial horizontal component velocity of each liquid particle will depend on the initial position within the A ellipse (Figure 1):

$$\left\{\begin{array}{l}{U}_0=\frac{x_c}{L_c}{W}_0,\\ V_0=\frac{y_c}{L_c}{W}_0,\end{array}\right.$$E1

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:

$$\left\{\begin{array}{l}{x}_c=2{\mathit{dx}}_o\left(\chi -0.5\right),\\ y_c=2{\left(\left({\mathit{dx}}_o^2-x_c^2\right)\left(\frac{{\mathit{dx}}_o^2}{{\mathit{dy}}_o^2}\right)\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(\chi -0.5\right),\end{array}\right.$$E2

where χ is a continuous uniform random variable in the [0, 1] interval whose average value is ${\mu }_{\chi }=0.5$ and standard deviation ${\sigma }_{\chi }=\sqrt{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 − L_c 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 $\overline{{\varphi }_0}=267.6\mathrm{\mu m}$, standard deviation of droplet diameters ${\sigma }_{\varphi }=110.3\mathrm{\mu }\mathrm{m},$ and a modal value $\mathit{Mo}=250.2\mathrm{\mu 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 $\overline{Z}=0$ and standard deviation ${\sigma }_Z=1$ can be obtained

$$Z=\frac{\overline{\chi }-{\mu }_{\chi }}{\frac{{\sigma }_{\chi }}{\sqrt{n}}},$$E3

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 ${\mu }_{\chi }=0.5$ and ${\sigma }_{\chi }=\sqrt{3}/6$.

With this method, it is possible to simulate diameters of liquid particles that follow the normal distribution from random number generation

$${\varphi }_0=\overline{{\varphi }_0}+{\sigma }_{\varphi }Z,$$E4

where $\overline{{\varphi }_0}$ and ${\sigma }_{\varphi }$ are the mean and standard deviation values of the Weibull p.d.f., respectively, whose expressions are given as:

$$\left\{\begin{array}{l}\overline{{\varphi }_0}=k.\mathrm{\Gamma }\left(1+\frac{1}{m}\right),\\ {\sigma }_{\varphi }=k.{\left\{\mathrm{\Gamma }\left(1+\frac{2}{m}\right)-{\left[\mathrm{\Gamma }\left(1+\frac{1}{m}\right)\right]}^2\right\}}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}.\end{array}\right.$$E5

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:

$$\mathit{Mo}=k{\left(\frac{m-1}{m}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$m$}\right.}.$$E6

Incorporating this value into Eq. (4), the random variable ϕ can be written as:

$$\left\{\begin{array}{ll}{\varphi }_0=\mathit{Mo}+{\sigma }_{\varphi }.Z,& \overline{{\varphi }_0}=\mathit{Mo},\text{symmetric}\left(\text{normal}\mathrm{p}.\mathrm{d}.\mathrm{f}.\right),\\ {\varphi }_0=\mathit{Mo}+{\sigma }_{\varphi }.Z.exp\left(\frac{Z}{2\left(\overline{{\varphi }_0}-\mathit{Mo}\right)}\right),& \overline{{\varphi }_0}\ne \mathit{Mo},\text{asymmetric}\left(\text{Weibull}\mathrm{p}.\mathrm{d}.\mathrm{f}.\right).\end{array}\right.$$E7

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:

$$\frac{F_{li}}{m_l}\hspace{0.17em}=\underset{\begin{array}{l}\text{Drag}\\ \text{acceleration}\end{array}}{\underbrace{\frac{U_i-V_i}{\tau }}}-\underset{\begin{array}{l}\text{Gravity}\\ \text{acceleration}\end{array}}{\underbrace{g_i{\delta }_{i3}}},$$E8

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:

$$\tau =\frac{4}{3}\frac{\phi }{C_D}\frac{{\rho }_l}{\rho }\frac{1}{\left|U_i-V_i\right|},$$E9

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 $\mathrm{\Re }{e}_l$:

$$\mathrm{\Re }{e}_l=\varphi \frac{\left|U_i-V_i\right|}{\nu }.$$E10

Several drag coefficient expressions have been analyzed [27]. The authors showed that Turton expression [28] has given better results in this simulation case:

$$\left\{\begin{array}{l}{C}_D=\frac{24}{\mathrm{\Re }{e}_l}\left(1+0.173\mathrm{\Re }{e}_l^{0.657}\right)+\frac{0.413}{1+16300\mathrm{\Re }{e}_l^{-1.09}},\text{if}\mathrm{\Re }{e}_l<2\times {10}^5\\ C_D=0.465,\text{if}\mathrm{\Re }{e}_l\ge 2\times {10}^5\end{array}\right..$$E11

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 $\left|U_i-V_i\right|{\delta }_{i3}=V_s$. From Eqs. (8) and (9):

$$V_s={\left(\frac{4}{3}\frac{{\rho }_l\varphi g_i{\delta }_{i3}}{\rho C_{D,s}}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}.$$E12

Note that $C_{D,s}$ is the drag coefficient at sedimentation regime. It depends on $\mathrm{\Re }{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:

$${\tau }_s=\frac{V_s}{g}.$$E13

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 U_i has a deterministic term and a random term:

$$\frac{{\mathit{dU}}_i}{\mathit{dt}}=\underset{\text{Deterministic}}{\underbrace{h_{\mathit{ij}}\left(U_i,t\right)}}+\underset{\text{Random}}{\underbrace{q_{\mathit{ij}}\left(U_i,t\right){\eta }_j\left(t\right)}}.$$E14

The tensors $h_{\mathit{ij}}\left(U_i,t\right)$ and $q_{\mathit{ij}}\left(U_i,t\right)$ 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), $\eta \left(t\right)$ denotes the random characteristic variable of zero mean and covariance:

$$⟨{\eta }_i\left(t^{\prime }\right){\eta }_j\left(t^{″}\right)⟩={\delta }_{\mathit{ij}}\delta \left(t^{\prime }-t^{″}\right).$$E15

2.2.2 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

:

$$U_i=u_i^{\oplus }+u_i^{-}.$$E16

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:

$$\frac{\partial {\tilde{u}}_i^{\oplus }}{\partial x_i}=0,$$E17

$$\frac{\partial {\tilde{u}}_i^{\oplus }}{\partial t}+\frac{\partial \left({\tilde{u}}_i^{\oplus }{u}_j^{\oplus }\right)}{\partial x_j}=\overline{\rho }{g}_i{B}^{\oplus }-\frac{\partial p{\prime }^{\oplus }}{\partial x_i}-\frac{\partial {\tilde{\tau }}_{\mathit{ij}}}{\partial x_j}+2\nu \frac{\partial {\tilde{S}}_{\mathit{ij}}^{a\oplus }}{\partial x_j}-{\mathfrak{I}}_i,$$E18

$$\frac{\partial {\tilde{\psi }}^{\oplus }}{\partial t}+\frac{\partial \left({\tilde{u}}_j^{\oplus }{\psi }^{\oplus }\right)}{\partial x_j}={\mathrm{\Phi }}_{\psi }-\frac{\partial {\tilde{\tau }}_{\mathit{i\psi }}}{\partial x_j},$$E19

where $B^{\oplus }$ is a buoyancy force, ${\tilde{S}}_{\mathit{ij}}^a$ is the anisotropic deformation tensor, $\nu$ 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 $\left({\tilde{u}}_i^{\oplus }=\overline{\rho }{u}_i^{\oplus }\right)$, 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 ${\tilde{\tau }}_{\mathit{ij}}$ and ${\tilde{\tau }}_{\mathit{i\psi }}$ 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:

$${\mathfrak{I}}_i=\frac{1}{\mathrm{\Delta }V}\sum _{l=1}^{n_{\mathrm{\Delta }}}{F}_{li}$$E20

where

is the force of the liquid particle l in i direction, ΔV = ΔxΔyΔz is the grid cell volume, and $n_{\mathrm{\Delta }}$ is the number of liquid particles within the grid cell. Using Eq. (8), the additional momentum is given by:

$${\mathfrak{I}}_i=\frac{1}{\mathrm{\Delta }V}\sum _{l=1}^{n_{\mathrm{\Delta }}}\left[m_l\left(\frac{U_i-V_i}{\tau }-g_j{\delta }_{j3}\right)\right].$$E21

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]:

$$\left\{\begin{array}{l}{h}_{\mathit{ij}}\left(U_i,t\right)=\frac{{\mathit{du}}_j^{\oplus }}{\mathit{dt}}+{\alpha }_{\mathit{ij}}{u}_j^{-},\\ q_{\mathit{ij}}\left(U_i,t\right)=\sqrt{C_0\epsilon }{\delta }_{\mathit{ij}}.\end{array}\right.$$E22

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₀ = 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:

$${\alpha }_{\mathit{ij}}=\frac{1}{2K^{-}}\frac{{\mathit{dK}}^{-}}{\mathit{dt}}{\delta }_{\mathit{ij}}-\left(\frac{3}{4}{C}_0\right)\frac{\epsilon }{K^{-}}{\delta }_{\mathit{ij}}+\left(\frac{R_{\mathit{ij}}}{2K^{-}}-\frac{{\delta }_{\mathit{ij}}}{3}\right)\frac{\epsilon }{K^{-}}.$$E23

The subgrid turbulent kinetic energy is solved by 1.5 order transport equation [31] and $R_{\mathit{ij}}={\left(u_i^{-}{u}_j^{-}\right)}^{\oplus }$ is the Reynolds SGS stress tensor.

The unresolved velocity component $u_j^{-}$ in Eq. (16) is obtained in discrete form using the Markov chains:

$$u_{j\left(n\right)}^{-}=u_{j\left(n-1\right)}^{-}+{\alpha }_{\mathit{ij}\left(n\right)}{u}_{j\left(n-1\right)}^{-}\mathrm{\Delta }t+\sqrt{C_0\epsilon \mathrm{\Delta }t}{\chi }_{\left(n\right)},$$E24

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]:

$$u_{j\left(0\right)}^{-}=\sqrt{\frac{2}{3}{K}_{\left(0\right)}^{-}}{\chi }_{\left(0\right)}.$$E25

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:

$$\begin{array}{l}\left\{\begin{array}{ll}{V}_{i\left(n+1\right)}=V_{i\left(n\right)}+\frac{\mathrm{\Delta }t}{\tau }\left(U_{i\left(n\right)}-V_{i\left(n\right)}\right)-g_i\mathrm{\Delta }t{\delta }_{i3},& \text{if}\tau >\mathrm{\Delta }t,\\ V_{i\left(n+1\right)}=U_{i\left(n\right)}-V_s{\delta }_{i3},& \text{if}\tau \le \mathrm{\Delta }t,\end{array}\right.\\ \left\{X_{i\left(n+1\right)}=\frac{V_{i\left(n+1\right)}-V_{i\left(n\right)}}{2}\mathrm{\Delta }t.\right.\hfill \end{array}$$E26

It is necessary to note that in the first Eq. (26), the rate $\mathrm{\Delta }t/\tau$ 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 $\mathrm{\Delta }t$ time step has elapsed. For this reason, the simulation time step $\mathrm{\Delta }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)):

$$\mathrm{\Delta }=\frac{{\varphi }_S}{{\varphi }_L},$$E27

$$\mathrm{Wes}=\frac{{\rho }_l{\varphi }_S{\mathrm{\Delta }}^3{\left|{\overrightarrow{V}}_{\mathit{mS}}\right|}^2+{\left|{\overrightarrow{V}}_{\mathit{mL}}\right|}^2}{12\sigma \mathrm{\Delta }\left(1+{\mathrm{\Delta }}^2\right)},$$E28

$$\mathrm{Imp}=\frac{2X}{{\varphi }_L+{\varphi }_S},$$E29

where the subscripts S and L indicate the smaller and larger droplet, respectively, σ denote the surface tension, ${\overrightarrow{V}}_{\mathit{mS}}$ and ${\overrightarrow{V}}_{\mathit{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 ${\overrightarrow{V}}_R={\overrightarrow{V}}_S-{\overrightarrow{V}}_L$ as shown in Figure 2.

$$\left\{\begin{array}{l}{\overrightarrow{V}}_{\mathit{mL}}={\overrightarrow{V}}_L-{\overrightarrow{V}}_{\mathit{mR}},\\ {\overrightarrow{V}}_{\mathit{mS}}={\overrightarrow{V}}_S-{\overrightarrow{V}}_{\mathit{mR}},\end{array}\right.$$E30

Figure 2.

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

In Eq. (30), ${\overrightarrow{V}}_{\mathit{mR}}$ is the velocity of mass center. If the droplets have the same density:

$${\overrightarrow{V}}_{\mathit{mR}}=\frac{{\overrightarrow{V}}_S{\varphi }_S^3+{\overrightarrow{V}}_L{\varphi }_L^3}{{\varphi }_S^3+{\varphi }_L^3}$$E31

The relative velocity droplets of the mass center can be resumed using Eq. (27):

$$\left\{\begin{array}{l}{\overrightarrow{V}}_{\mathit{mL}}=+\frac{{\mathrm{\Delta }}^3{\overrightarrow{V}}_{\mathit{mR}}}{{\mathrm{\Delta }}^3+1},\\ {\overrightarrow{V}}_{\mathit{mS}}=-\frac{{\overrightarrow{V}}_{\mathit{mR}}}{{\mathrm{\Delta }}^3+1}.\end{array}\right.$$E32

For the impact parameter in Eq. (29), it is necessary to compute X variable. This variable is the projection of segment $b=0.5\left({\varphi }_S+{\varphi }_L\right)$ on the plane perpendicular to the relative velocity ${\overrightarrow{V}}_R$. The impact factor will be equal to the cosine of γ angle:

$$\mathrm{Imp}=\mathit{cos}\gamma =\frac{X}{b},$$E33

$$\mathit{sin}\gamma =\frac{D_p}{b},$$E34

So, inserting Eq. (34) into Eq. (33) results in:

$$\mathrm{Imp}=\mathit{cos}\gamma ={\left[1-{\left(\frac{2D_p}{{\varphi }_S+{\varphi }_L}\right)}^2\right]}^{\frac{1}{2}}.$$E35

It is necessary to obtain D_p. It is the distance from large droplet center $\left(x_L,y_L,z_L\right)$ to the perpendicular plane at ${\overrightarrow{V}}_R$ velocity passing through the small droplet center $\left(x_S,y_S,z_S\right)$. This plane P-P is shown in green color in Figure 2. The plane equation passing through the small drop center is $u_R{x}_S+v_R{y}_S+w_R{z}_S+D=0,$ where $\left(u_R,v_R,w_R\right)$ are the components of ${\overrightarrow{V}}_R$. 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:

$$D_p=\left|\frac{u_R{x}_L+v_R{y}_L+w_R{z}_L+D}{{\left({u_R}^2+{v_R}^2+{w_R}^2\right)}^{\frac{1}{2}}}\right|.$$E36

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:

$${\left(\frac{{\varphi }_S+{\varphi }_L}{2}\right)}^2={\left(x_S-x_L-u_R\mathrm{\Delta }{t}^{\prime \prime }\right)}^2+{\left(y_S-y_L-v_R\mathrm{\Delta }{t}^{\prime \prime }\right)}^2+{\left(z_S-z_L-w_R\mathrm{\Delta }{t}^{\prime \prime }\right)}^2.$$E37

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 (Imp_c-s) is according to Rabe [20] as follows:

$${\mathrm{Imp}}_{\mathrm{c}‐\mathrm{s}}=\frac{\sqrt{{0.53}^2+4.24\mathrm{Wes}}-0.53}{4\mathrm{Wes}}.$$E38

The transition impact factor between coalescence and reflexive separation (Imp_c-r) is:

$${\mathrm{Imp}}_{\mathrm{c}‐\mathrm{r}}=0.3059\sqrt{1-\frac{0.45}{\mathrm{Wes}}}.$$E39

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, φ:

$${\text{We}}_{\mathrm{b}}=\frac{\mathrm{\Delta }\left(1+{\mathrm{\Delta }}^2\right)\left(4\phi -12\right)}{\xi \left(1-{\mathrm{Imp}}^2\right)},$$E40

where ξ is computed as:

$$\xi =\left\{\begin{array}{ll}1-\frac{{\left(2-\lambda \right)}^2\left(1+\lambda \right)}{4}& \text{if}\lambda >1,\\ \frac{{\lambda }^2\left(3-\lambda \right)}{4}& \text{if}\lambda \le 1,\end{array}\right.$$E41

and λ = (1 − Imp)(1 + ∆). The shape parameter φ can be computed as Zhang [19]:

$$\phi =3.351{\left(\frac{{\rho }_l}{1.16}\right)}^{\frac{2}{3}}.$$E42

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:

$${\mathrm{Wes}}_{\mathrm{b}}={\text{We}}_{\mathrm{b}}\frac{{\mathrm{\Delta }}^2}{12\left(1+{\mathrm{\Delta }}^3\right)\left(1+{\mathrm{\Delta }}^2\right)},$$E43

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

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

2.3.4.1 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 ${\overrightarrow{V}}_{\left(\mathit{new}\right)}={\overrightarrow{V}}_{\mathit{mR}}$ (Eq. (31)). For droplets of the same density, its diameter is ${\varphi }_{\left(\mathit{new}\right)}={\left({\varphi }_L^3+{\varphi }_S^3\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$.

2.3.4.2 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):

$$C_{\mathit{VS}}=\frac{{\mathrm{E}}_{\text{strtch}}-{\mathrm{E}}_{\text{surten}}-{\mathrm{E}}_{\text{dissip}}}{{\mathrm{E}}_{\text{strtch}}+{\mathrm{E}}_{\text{surten}}+{\mathrm{E}}_{\text{dissip}}},$$E44

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:

$$Z=\frac{\mathrm{Imp}-\sqrt{e_{\mathit{\text{coal}}}}}{1-\sqrt{e_{\mathit{\text{coal}}}}},$$E45

where $e_{\mathit{\text{coal}}}=min\left[1.0,2.4f{\text{We}}^{-1}\right]$ is the coalescence efficiency, $f={\mathrm{\Delta }}^{-3}-2.4{\mathrm{\Delta }}^{-2}+2.7{\mathrm{\Delta }}^{-1}$, and $\text{We}={\rho }_l{\varphi }_S{\mathrm{\Delta }}^3{\sigma }^{-1}{\left|{\overrightarrow{V}}_R\right|}^2,$ 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:

$$\left\{\begin{array}{l}{\overrightarrow{V}}_{\mathit{mS}\left(\mathit{new}\right)}=Z{\overrightarrow{V}}_{\mathit{mS}}\\ {\overrightarrow{V}}_{\mathit{mL}\left(\mathit{new}\right)}=Z{\overrightarrow{V}}_{\mathit{mL}}\end{array}\right..$$E46

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:

$$\frac{1}{6}\pi \left({\mathrm{\Psi }}_S{\varphi }_S^3+{\mathrm{\Psi }}_L{\varphi }_L^3\right)=\pi r_0^2\eta ,$$E47

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\sqrt{{\text{We}}_0}$. If T ≤ 2, the ligament contracts in a single satellite whose radius is r_o. k₂ = 0.45 and ${\text{We}}_0=2r_0\left({\rho }_l/\sigma \right){\left|{\overrightarrow{V}}_R\right|}^2$. Otherwise, it is necessary to compute the evolution time of ligament radius equation for obtaining the final value r_bu:

$$\frac{3}{4\sqrt{2}}{k}_1{k}_2\sqrt{{\text{We}}_0}{\left(\frac{r_{\mathit{bu}}}{r_0}\right)}^{\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+{\left(\frac{r_{\mathit{bu}}}{r_0}\right)}^2-1=0,$$E48

Figure 7.

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

k₁ = 11.5 following Kim [18]. Eq. (48) can be solved by iteration with an initial value $\left(r_{\mathit{bu}}/r_0\right)=1$ and ∆t = 1 × 10⁻².

The diameter of satellite droplets can be determined by following Georjon [41]:

$${\varphi }_{\mathit{sat}}=3.78r_{\mathit{bu}}.$$E49

The number of satellite droplets is calculated from the mass conservation by assuming uniform satellites size $N_{\mathit{sat}}=6{\left(r_0/{\varphi }_{\mathit{sat}}\right)}^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).

$${\overrightarrow{V}}_{\mathit{sat}}=\frac{\left({\varphi }_L^3{\overrightarrow{V}}_L-{\varphi }_{L\left(\mathit{new}\right)}^3{\overrightarrow{V}}_{L\left(\mathit{new}\right)}\right)+\left({\varphi }_S^3{\overrightarrow{V}}_S-{\varphi }_{\mathrm{S}\left(\mathit{new}\right)}^3{\overrightarrow{V}}_{\mathrm{S}\left(\mathit{new}\right)}\right)}{N_{\mathit{sat}}{\varphi }_{\mathit{sat}}^3},$$E50

where the diameters of parent droplets after collision are:

$$\left\{\begin{array}{l}{\varphi }_{L\left(\mathit{new}\right)}=\left(1-{\mathrm{\Psi }}_L\right)\frac{1}{6}{\mathit{\pi \varphi }}_L^3\\ {\varphi }_{\mathrm{S}\left(\mathit{new}\right)}=\left(1-{\mathrm{\Psi }}_S\right)\frac{1}{6}{\mathit{\pi \varphi }}_S^3\end{array}\right..$$E51

2.3.4.3 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 $r_0^3={\left({\varphi }_L/2\right)}^3+{\left({\varphi }_S/2\right)}^3$. 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 $N_{\mathit{sat}}=6{\left(r_0/{\varphi }_{\mathit{sat}}\right)}^3-2$. If N_sat ≤ 0, it is assumed that the ligament breaks up without satellite droplet and the two end-droplets 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:

$${\varphi }_{L\left(\mathit{new}\right)}={\varphi }_{S\left(\mathit{new}\right)}={\left(\frac{8r_0^3-{\varphi }_{\mathit{sat}}^3}{2}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}.$$E52

When N_sat > 1, the ligament breaks up into uniform droplets with identical diameters ${\varphi }_{L\left(\mathit{new}\right)}={\varphi }_{S\left(\mathit{new}\right)}={\varphi }_{\mathit{sat}}$. The velocities of end-droplets and satellite droplets are computed with Eqs. (46) and (50).

2.3.4.4 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:

$$\left\{\begin{array}{l}{V}_{nL}=V_{iL}\frac{\left|x_{iR}\right|}{b}\\ V_{nS}=V_{iS}\frac{\left|x_{iR}\right|}{b}\end{array}\right.,$$E53

where $x_{iR}=x_{iS}-x_{iL}$ 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:

$$\left\{\begin{array}{l}{V}_{nL\left(\mathit{new}\right)}=V_{nL}-\frac{{\mathrm{\Delta }}^3}{1+{\mathrm{\Delta }}^3}{f}_n^{\mathit{dd}}\\ V_{nS\left(\mathit{new}\right)}=V_{nS}+\frac{1}{1+{\mathrm{\Delta }}^3}{f}_n^{\mathit{dd}}\end{array}\right.,$$E54

where $f_n^{\mathit{dd}}$ is the normal impulse of a droplet-droplet collision: $f_n^{\mathit{dd}}=-\left(1+e_{n,p}\right)\left(V_{nS}-V_{nL}\right)$.

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 R² = 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⁻¹ 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 R² = 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⁻¹ 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.

Acknowledgments

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 ${\varphi }_0$ that represents the initial diameter of the liquid particles, and the expression of Weibull p.d.f. for this random variable is:

$$\left\{\begin{array}{ll}f\left({\varphi }_0\right)=\frac{m}{k^m}{\varphi }_0^{m-1}exp\left[-{\left(\frac{{\varphi }_0}{k}\right)}^m\right],& \text{if}{\varphi }_0\ge 0,\\ f\left({\varphi }_0\right)=0,& \text{if}{\varphi }_0<0.\end{array}\right.$$E55

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 ${\varphi }_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\left({\varphi }_0\right)$ is expressed as:

$$F\left({\varphi }_0\right)=1-exp\left[-{\left(\frac{{\varphi }_0}{k}\right)}^m\right].$$E56

Eq. (A.2) can be written as:

$$\mathrm{Ln}\left\{-\mathrm{Ln}\left[1-F\left({\varphi }_0\right)\right]\right\}=m.\mathrm{Ln}\left({\varphi }_0\right)-m.\mathrm{Ln}\left(k\right).$$E57

If we associate Eq. (A.3) with a linear equation f(X)  = mX + a where the dependent variable is $f\left(X\right)=\mathrm{Ln}\left\{-\mathrm{Ln}\left[1-F\left({\varphi }_0\right)\right]\right\}$ and the independent variable is $X=\mathrm{Ln}\left({\varphi }_0\right)$, 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 $\mathrm{Ln}\left\{-\mathrm{Ln}\left[1-F\left({\varphi }_0\right)\right]\right\}$ vs. $\mathrm{Ln}\left({\varphi }_0\right)$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 $k=exp\left(-a/m\right)$.

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.