Approximate Closed-Form Solution of the Differential Equation Describing Droplet Flight during Sprinkler Irrigation

: Sprinkler irrigation is widely used in agriculture because it allows for rational use of water. However, it can induce negative effects of soil erosion and of surface waterproofing. The scholars of these phenomena use the numerical integration of the equation of motion, but if there was an analytical solution, the study would be facilitated, and this solution could be used as software for regulating sprinklers. Therefore, in this study, the solution of the differential equation of the flight of droplets produced by sprinklers in the absence of wind was developed. The impossibility of an exact analytical solution to the ballistic problem due to the variability of the drag coefficient of the droplets is known; therefore, to find the integrals in closed form, the following were adopted: a new formula for the drag coefficient; a projection of the dynamic’s equation onto two local axes, one tangent and one normal to the trajectory and some linearization. To reduce the errors caused by the latter, the linearization coefficients and their calculation formulas were introduced through multiple non-linear regressions with respect to the jet angle and the initial droplet speed. The analytical modeling obtained, valid for jet angles from 10 ◦ to 40 ◦ , was compared to the exact numerical solution, showing, for the total travel distance, a high accuracy with a mean relative error MRE of 1.8% ± 1.4%. Even the comparison with the experimental data showed high accuracy with an MRE of 2.5% ± 1.1%. These results make the analytical modeling capable of reliably calculating the travel distance, the flight time, the maximum trajectory height, the final fall angle and the ground impact speed. Since the proposed analytical modeling uses only elementary functions, it can be implemented in PLC programmable logic controllers, which could be useful for controlling water waste and erosive effects on the soil during sprinkler irrigation.


Introduction
Sprinkler irrigation is widely used in agriculture because it allows a rational use of water.However, the high specific power (W m −2 ), a product of kinetic energy (J m −3 ) with an application rate (m 3 s −1 m −2 ) of water droplets impacting soil, is a factor to be controlled because it can induce the negative effects of soil erosion and surface sealing [1].The landing angle of the droplets on the soil surface is also a factor to be controlled because low values induce the negative effects of soil erosion and surface sealing [1].
In recent years, scientific research in this area has been oriented towards study, mathematical modeling and experimentation to reduce the phenomena of soil sealing and erosion and to increase distribution uniformity.Hui et al. [1] studied the influence of low sprinkler pressure on soil erosion.Chen et al. [2] investigated erosion due to droplet kinetic energy for different soil types and indicated the importance of pressure, nozzle diameter and sprinkler spacing.Hui et al. [3], after confirming that the phenomenon of soil erosion depends more on the impact angle than on the kinetic energy of the droplets, identify low-pressure sprinklers as a possible solution.Zhu et al. [4] quantified the angle of impact on crops such as corn and the associated splash production resulting in soil erosion.Chen et al. [5] Inventions 2024, 9, 73 2 of 26 studied the influence of the kinetic energy of droplets produced by low-pressure sprinklers on soil erosion.Ge et al. [6] set up a method to evaluate the kinetic energy, the diameter of the droplets and the distribution of their speeds.Félix-Félix et al. [7] characterized the morphology and velocity of the droplets in the laboratory with the aim of improving sprinkler simulation models.Hua et al. [8] proposed three-dimensional mathematical modeling of sprinkler rainfall, based on numerically solved ballistics equations, to estimate the influence of rotation speed and performed an experimental verification of it.
In many of these works, it was necessary to use the differential equations of droplet motion, solved through numerical methods [9][10][11][12].In fact, the numerical solution allows for the simulation of the characteristics of the droplet trajectory and, therefore, the total travel distance, the height of the trajectory, the flight time, the landing angle and the terminal velocity.The differential equations of droplet motion during sprinkler irrigation are derived from Newton's second law of dynamics, in which the system of external forces applied to the droplet is equal to the product of mass and acceleration [13][14][15][16][17][18][19].The forces applied are gravity, buoyancy and drag.While the first two are constant during the motion, the drag force depends on the droplet velocity with a complex and non-unique law: this is the reason why the exact closed-form solution of the droplet flight dynamic equations is impossible.This is why the numerical solution is necessary; however, in the recent past, two proposals for approximate analytical solutions have been made.In the first proposal [20], the closed-form solution was possible because two simplifying hypotheses were adopted.The first was on the simplified formulation of the drag coefficient, while the second concerned the way to project the vector equation of dynamics onto the x and y axes, multiplying the drag force by cos 2 θ and by sin 2 θ instead of by cos θ and by sin θ, respectively.These simplifications, on the one hand, had allowed an easy closed-form solution, but on the other, had led to a certain error then calculated by Yan et al. [21].In the second proposal [22], the solution was obtained by imposing that the drag coefficient followed Stokes' law, i.e., the drag force linear function of the droplet velocity.This hypothesis, on the one hand, allowed for easy analytical integration but, on the other hand, introduced an easily reliable error if one considers that Stokes' law is valid for 0.1 ≤ Re ≤ 2 while during the motion of the droplets, the Reynolds number Re is many times higher than this range, often by two or three orders of magnitude.
Comparing the two methods of solving the equations of motion, the analytical one and the numerical one, the advantage of the closed-form solution is that this is a set of equations based on elementary functions (trigonometric, logarithmic and exponential) that can be solved in sequence and instantaneously, using a double dozen cells in a spreadsheet.A pocket scientific calculator can also be used.Furthermore, it must be considered that equations based on elementary functions can be implemented in PLC (programmable logic controller).This possibility is precluded by numerical integration methods and, furthermore, these methods are not instantaneous: for each condition to be analyzed, modern software requires, on average, tens of seconds of processing.Finally, software for numerical integration can have high investment costs.
Therefore, a closed-form solution of the droplet flight dynamic equation, approximated by minimizing errors, will be carried out through a series of measures: the development of a new single formula for the drag coefficient valid for the entire droplet flight; a projection of the dynamic equation onto two local axes, one tangent and one normal to the trajectory; some appropriate linearization [23] introducing linearization coefficients as a function of the jet angle and the initial speed, with formulas obtained through multiple non-linear regressions.

Materials and Methods
Section 2.1 is dedicated to the setting of the two flight differential equations, one for the ascending portion of the trajectory and one for the descending portion, and to their partially approximate analytical solution.The integration procedure of these two flight equations is long and divided into many steps and is reported to allow interested readers to follow the procedure and the phases dedicated to the approximations.Instead, readers interested only in the final equations to determine the geometric and dynamic values of the trajectories useful for the study and mitigation of the soil erosion phenomenon can find them grouped in Section 3.4, which follows.They are to be used on the spreadsheet or inserted into the PLC software, i.e., into the microprocessors that can be installed in the sprinklers for their adjustment aimed at mitigating soil erosion during irrigation.Section 2.2 is dedicated to the description of the exact numerical integration of the two flight equations.The exact results of the numerical integration constitute the benchmark to evaluate (Section 3.1) the accuracy of the results of the analytical modeling equations developed in Section 2.1 and summarized in Section 3.4.

Closed-Form Solution and Mathematical Modeling
When the water jet comes out of the sprinkler nozzle, it is compact but then breaks into ligaments and then into droplets, thanks to the interaction with the ambient air.This is a complex phenomenon which also sees a subsequent interaction between the droplets themselves, already formed.As always happens in engineering, to study such a high complexity, some simplifying hypotheses [9] are posed to carry out mathematical modeling.If necessary, after comparing the solution and the experimental results, correction coefficients can be introduced.
Therefore, to set the differential equation of the flight of the droplets emitted from the nozzle, it was necessary to define some simplifying conditions: 1.
Flow of water completely disintegrates into droplets of various diameters in the outlet section of the sprinkler nozzle; 2.
All droplets have the same initial velocity, which is equal to the average velocity of the water flow in the exit section of the nozzle; 3.
Each droplet of water is considered alone along the trajectory, i.e., without collisions with the other droplets; 4.
Water droplets are considered rigid spheres upon exiting the nozzle and remain spherical during flight; 5.
Water droplets are not subject to evaporation, so their diameters remain unchanged throughout their trajectory; 6.
No wind.
Figure 1 shows the trajectory of a droplet emitted from the nozzle of a sprinkler.

Materials and Methods
Section 2.1 is dedicated to the setting of the two flight differential equations, one for the ascending portion of the trajectory and one for the descending portion, and to their partially approximate analytical solution.The integration procedure of these two flight equations is long and divided into many steps and is reported to allow interested readers to follow the procedure and the phases dedicated to the approximations.Instead, readers interested only in the final equations to determine the geometric and dynamic values of the trajectories useful for the study and mitigation of the soil erosion phenomenon can find them grouped in Section 3.4, which follows.They are to be used on the spreadsheet or inserted into the PLC software, i.e., into the microprocessors that can be installed in the sprinklers for their adjustment aimed at mitigating soil erosion during irrigation.Section 2.2 is dedicated to the description of the exact numerical integration of the two flight equations.The exact results of the numerical integration constitute the benchmark to evaluate (Section 3.1) the accuracy of the results of the analytical modeling equations developed in Section 2.1 and summarized in Section 3.4.

Closed-Form Solution and Mathematical Modeling
When the water jet comes out of the sprinkler nozzle, it is compact but then breaks into ligaments and then into droplets, thanks to the interaction with the ambient air.This is a complex phenomenon which also sees a subsequent interaction between the droplets themselves, already formed.As always happens in engineering, to study such a high complexity, some simplifying hypotheses [9] are posed to carry out mathematical modeling.If necessary, after comparing the solution and the experimental results, correction coefficients can be introduced.
Therefore, to set the differential equation of the flight of the droplets emitted from the nozzle, it was necessary to define some simplifying conditions: 1. Flow of water completely disintegrates into droplets of various diameters in the outlet section of the sprinkler nozzle; 2. All droplets have the same initial velocity, which is equal to the average velocity of the water flow in the exit section of the nozzle; 3.Each droplet of water is considered alone along the trajectory, i.e., without collisions with the other droplets; 4. Water droplets are considered rigid spheres upon exiting the nozzle and remain spherical during flight; 5. Water droplets are not subject to evaporation, so their diameters remain unchanged throughout their trajectory; 6.No wind.
Figure 1 shows the trajectory of a droplet emitted from the nozzle of a sprinkler.It is possible to divide the trajectory into a first part from the emission point of the nozzle 0, which is at a height h (m) from the ground up to the top, and into a second part from the top up to the landing.The jet angle θ 0 ( • ) is the initial angle of the droplet jet, equal to the nozzle angle.The most important expected results of the mathematical modeling are the droplet travel distance x tot (m), the flight time t tot (s), the terminal velocity v end (m s −1 ) and the ground impact angle θ end ( • ), which depend on droplet diameter D (m); initial velocity of the droplet v 0 (m s −1 ); jet angle θ 0 ; nozzle height from the ground h (m); characteristics of air, such as viscosity µ (Pa s) and density ρ air (kg m −3 ).

Ascending Portion of the Trajectory
Considering the ascending portion (Figure 1) from the sprinkler (0, h) to the top (x top , y top ), the differential equation of dynamics is where a = dv dt is the acceleration vector (m s −2 ), m D •g is the gravitational force vector (N), m D is the mass of the droplet (kg), m air •g is the buoyancy force (N), m air is the mass (kg) of the air moved by the droplet, R = c d π 4 D 2 ρ air v 2 v is the drag force vector (N), which is proportional to the drag coefficient c d , to the equatorial area of the droplet of diameter D and to the kinetic energy (J).The air density is ρ air (kg m −3 ).The buoyancy force is three orders of magnitude less than the gravitational force, so it can be neglected: (m D − m air ∼ = m D ).
Generally, to proceed with numerical integration and as was performed in a previous simplified analytical modeling [20], the vector Equation ( 1) is projected along the horizontal x and vertical y axes.Instead, in the present study, the projection of dynamics equation onto a local reference system having a tangent axis s and a normal axis n to the trajectory was chosen (Figure 2).
It is possible to divide the trajectory into a first part from the emission point of the nozzle 0, which is at a height h (m) from the ground up to the top, and into a second part from the top up to the landing.The jet angle θ0 (°) is the initial angle of the droplet jet, equal to the nozzle angle.The most important expected results of the mathematical modeling are the droplet travel distance xtot (m), the flight time ttot (s), the terminal velocity vend (m s −1 ) and the ground impact angle θend (°), which depend on droplet diameter D (m); initial velocity of the droplet v0 (m s −1 ); jet angle θ0; nozzle height from the ground h (m); characteristics of air, such as viscosity µ (Pa s) and density ρair (kg m −3 ).

Ascending Portion of the Trajectory
Considering the ascending portion (Figure 1) from the sprinkler (0, h) to the top (xtop, ytop), the differential equation of dynamics is where  = is the acceleration vector (m s −2 ),  • ̅ is the gravitational force vector (N), mD is the mass of the droplet (kg),  • ̅ is the buoyancy force (N), mair is the mass (kg) of the air moved by the droplet,  =    ̅ is the drag force vector (N), which is proportional to the drag coefficient cd, to the equatorial area of the droplet of diameter D and to the kinetic energy (J).The air density is ρair (kg m −3 ).The buoyancy force is three orders of magnitude less than the gravitational force, so it can be neglected: Generally, to proceed with numerical integration and as was performed in a previous simplified analytical modeling [20], the vector Equation ( 1) is projected along the horizontal x and vertical y axes.Instead, in the present study, the projection of dynamics equation onto a local reference system having a tangent axis s and a normal axis n to the trajectory was chosen (Figure 2). .The local reference system adopted, with the s and n axes that are, respectively, tangent and normal to the trajectory.The red vector is the drag force  , the black vector is the gravitational force ̅ projected in the two components along s ( ̅ sin ) and n ( ̅ cos ).The blue vector is the velocity ̅ .
Figure 3 shows that the scalar velocity variation in the s direction is dv, and therefore, the acceleration is , while in the n direction, the scalar velocity variation is v•dθ; therefore, the acceleration is  .Therefore, Equation (1) projected onto s and n, neglecting the buoyancy force, which is three orders of magnitude less than the gravitational force ( −  ≅  ), becomes .The local reference system adopted, with the s and n axes that are, respectively, tangent and normal to the trajectory.The red vector is the drag force R, the black vector is the gravitational force mg projected in the two components along s (m D gsin θ) and n (m D gcos θ).The blue vector is the velocity v.
Figure 3 shows that the scalar velocity variation in the s direction is dv, and therefore, the acceleration is dv dt , while in the n direction, the scalar velocity variation is v•dθ; therefore, the acceleration is v dθ dt .Therefore, Equation (1) projected onto s and n, neglecting the buoyancy force, which is three orders of magnitude less than the gravitational force The convention that the angle θ starts from the jet angle θ 0 value and takes on a positive value during the ascending portion of the trajectory is adopted.
In Equation ( 3), the component of the gravitational force has a positive sign because its direction is concordant with the n-axis.In Equation ( 2), the drag force appears in which the drag coefficient c d is not constant but depends on the velocity v through the Reynolds number.Bird et al. (1960) [24] indicated that, with laminar motion (0.1 ≤ Re ≤ 2) in the boundary layer, the Stokes formula c d = 24  Re applies; with transition motion (2 < Re ≤ 500), the formula c d = 18.5 Re 0.6 applies; with turbulent motion (Re > 500) , the formula c d = 0.44 applies.These formulas are valid for rigid spheres.The convention that the angle θ starts from the jet angle θ0 value and takes on a positive value during the ascending portion of the trajectory is adopted.
In Equation ( 3), the component of the gravitational force has a positive sign because its direction is concordant with the n-axis.In Equation ( 2), the drag force appears in which the drag coefficient cd is not constant but depends on the velocity v through the Reynolds number.Bird et al. (1960) [24] indicated that, with laminar motion (0.1 ≤  ≤ 2) in the boundary layer, the Stokes formula  = applies; with transition motion (2  ≤ 500), the formula  = .
. applies; with turbulent motion ( 500) , the formula  = 0.44 applies.These formulas are valid for rigid spheres.
Since, during the motion of the droplet, the Reynolds number varies, it may be that the choice of a relationship between the three indicated by Bird et al. [24] in order to perform the integration of Equation ( 2) produces important errors in the prediction of the droplet travel distance and the flight time.
To overcome this difficulty, in the present analytical model, a unique relationship has been applied.Recently [25], for 0.1 ≤  ≤ 200,000, the following formula was developed: = 26.31•  − 0.000037447• 4 −0.00066989• 3 +0.0016779• 2 −0.033243•+0.86961, where  = ln().However, since in the motion of the droplets emitted by the sprinklers, the Reynolds number does not exceed the value of 10,000, the previous relationship has been simplified, always using the Standard Drag Curve data (SDC) [25,26]: It is valid for 0.3 ≤  ≤ 10,000 and, compared to the SDC data, it has a relative mean error-MRE = 3.8% and a standard deviation-SD = 3.0%.In the same Reynolds range, the three relationships proposed by Bird et el.[24] present MRE = 3.1% and SD = 3.0%, i.e., values similar to those obtained with (4).
Equation ( 2) with (4), after dividing by the mass of the droplet mD, becomes Equation ( 3) gives = • cos , which is inserted into (5).Furthermore, the velocity variable v is replaced with the Reynolds number variable Re,  =

•
, where µ is the air viscosity (Pa s).Finally, instead of the mass of the spherical droplet, the expression is in- , where ρw is the density of water (kg m −3 ).Equation ( 5) becomes Since, during the motion of the droplet, the Reynolds number varies, it may be that the choice of a relationship between the three indicated by Bird et al. [24] in order to perform the integration of Equation ( 2) produces important errors in the prediction of the droplet travel distance and the flight time.
To overcome this difficulty, in the present analytical model, a unique relationship has been applied.Recently [25], for 0.1 ≤ Re ≤ 200, 000, the following formula was developed: 86961) , where G = ln(Re).However, since in the motion of the droplets emitted by the sprinklers, the Reynolds number does not exceed the value of 10,000, the previous relationship has been simplified, always using the Standard Drag Curve data (SDC) [25,26]: It is valid for 0.3 ≤ Re ≤ 10, 000 and, compared to the SDC data, it has a relative mean error-MRE = 3.8% and a standard deviation-SD = 3.0%.In the same Reynolds range, the three relationships proposed by Bird et el.[24] present MRE = 3.1% and SD = 3.0%, i.e., values similar to those obtained with (4).
Equation ( 2) with (4), after dividing by the mass of the droplet m D , becomes Equation ( 3) gives dθ dt = g v •cosθ, which is inserted into (5).Furthermore, the velocity variable v is replaced with the Reynolds number variable Re, v = Re µ ρ air •D , where µ is the air viscosity (Pa s).Finally, instead of the mass of the spherical droplet, the expression is inserted, m D = ρ w πD 3   6 , where ρ w is the density of water (kg m −3 ).Equation ( 5) becomes where the quantity Ar is the Archimedes number [25]: Equation ( 6) is a non-integrable differential equation due to the exponent (0.04925•lnRe − 0.90977 + 3).However, through some parallel numerical integration tests of Equation (6), it was seen that the maximum decrease in Re occurs with the minimum droplet diameter of 0.25 mm and is approximately 30-fold during the ascending part of the trajectory.Consequently, the maximum variation of the exponent (0.04925•lnRe − 0.90977 + 3) is 7%, while it reduces to 2% with large droplet diameter (D = 5 mm).
Consequently, the simplifying hypothesis of considering the exponent calculated with the Reynolds number equal to the initial value Re 0 = ρ air v 0 D µ was adopted.Therefore, Equation (6) becomes where β = 0.04925•lnRe 0 − 0.90977 + 3.By introducing the new variable u a = 1 Re β−1 , so du a dθ = −(β−1) Re β dRe dθ , Equation ( 7) can be transformed: Equation ( 8) is a first-order linear differential equation that can be solved by the Bernoulli method, that is, by introducing the variable u a as the product of two dummy variables w a and z a : u a = w a •z a .Therefore, Equation ( 8) becomes A function z a = z a (θ) must be found such that its derivative satisfies this equality: It is sufficient to integrate, by separation of the variables, Equation ( 10) to obtain the desired function: where C za is the constant of integration.The combination of Equations ( 9)- (11) gives Since β is between approximately 2.3 and 2.6 depending on the droplet diameter D and the initial velocity v 0 , (β − 1) is between 1.3 and 1.6, i.e., non-integer number, so the differential equation (12) has no solution in closed form.However, the expression (cos θ) (β−1) can be replaced by α (cos θ) 2 , where the coefficient α depends on the initial values of the Reynolds number Re 0 , the jet angle θ 0 and the water pressure in the nozzle p.The relationship α = α(Re 0 , θ 0 , p) will be obtained later, in Appendix A, using a nonlinear multiple regression.It will be performed in such a way as to also consider the errors due to the approximation introduced in Equation ( 7) with the exponent β = constant.Thus, Equation (12) becomes The solution of which, containing the quantity C wa , which is the value assumed by the variable w a when θ = θ 0 , is Combining Equations ( 14) and (11), the variable u a becomes where 1) , Equation (15) becomes (16) By imposing the initial conditions, θ = θ 0 , Re = Re 0 , the constant C ua is obtained, as follows: Setting θ = 0, Equation ( 16) gives the value of the Reynolds number at the top of the trajectory: From Equation ( 16), it is easy to obtain the velocity of the droplet with diameter D: The integral of Equation ( 19) provides the length s a (m) of the trajectory in the ascending portion as a function of the time and angle θ: The combination of Equation (19) with the expression dθ dt = g v •cosθ obtained from Equation (3) gives the elemental time interval dt as a function of angle θ: Finally, the combination of Equations ( 20) and (21) gives s a as a function of only the angle θ: Performing the integral, Equation ( 22) becomes where C s is the constant of integration, which is determined by imposing the initial conditions θ = θ 0 , s a = 0. Furthermore, by imposing θ = 0, Equation (23) gives the length of the ascending trajectory up to the top: If the function inside the integral of Equation ( 22) were multiplied by cos θ, the integration would provide the droplet travel distance x top or, if it were multiplied by sin θ, the integration would provide the height y top (Figure 1).Unfortunately, the two new functions cannot be integrated; therefore, simplified calculation methods of x a and y a are proposed.It is valid because it involves mean relative errors of approximately 1%.The variables s a and, therefore, x a and y a of the ascending portion originate at the sprinkler nozzle (Figure 1).
The method consists of calculating the arithmetic mean value, γ, of cos θ 0 and cos 0, which is γ = cos θ 0 +1 2 , and multiplying it by s a , obtaining the droplet travel distance x top .To obtain y top , first, the average angle θ ave is calculated, then its trigonometric tangent is calculated, which is the ratio between y top and x top .The formulas are To obtain the flight time to the top, Equation ( 21) could be used, which, however, cannot be integrated.However, with the precaution of multiplying the function inside the integral of Equation ( 21) by the cos θ, a solution is obtained, which must be corrected with a coefficient φ to be obtained (Appendix A) via the regression method, similar to the coefficient α.Therefore, the following equation gives the flight time t top :

.2. Descending Portion of the Trajectory
The descending portion presents a differential equation of droplet motion like Equation ( 1), but with the difference that the sign of the gravitational force is positive because it becomes a driving force (Figure 1): Therefore, Equation ( 7) also changes the sign as follows: where the exponent δ depends on the Reynolds number at the top of the trajectory Re top , already determined by Equation ( 18): δ = 0.04925•lnRe top − 0.90977 + 3.
The convention is adopted that the angle θ starts from a zero value at the top of the trajectory and takes on a positive value during the descending portion.
Also, considering the dummy variable for the descending segment, 9)-( 12) become Equation ( 32), like Equation (12), cannot be integrated because δ is not integer.But, even in this case, an approximation such as (cos θ) δ ∼ = ε(cos θ) 2 , where ε is a coefficient obtainable (Appendix A) through regression method like α, makes Equation (32) integrable, as follows: The solution is where C wd is the value of the variable w d at the beginning of the descending trajectory, i.e., at the top when θ = θ top = 0.The multiplication of Equations ( 31) and (34) gives the dummy variable u d relating to the descending section: The new constant of the descending portion is 35) becomes for the descending portion: The constant C ud is determined with the initial condition of the descending portion: θ = 0, Re = Re top : From Equation (36), it is easy to obtain the elemental displacement ds d of the droplet of diameter D along the descending portion of the trajectory: The combination of Equation (36) with the expression dθ dt = g v •cosθ, obtained from Equation (3), gives the following: By inserting Equation (39) into (38) and proceeding on to integration, the path length s d (θ) of the droplet along the descending trajectory can be obtained: Unfortunately, Equation (40) has no solution due to the term (cos θ) 3 and the non- integer exponent 2  δ−1 .To overcome the first difficulty, (40) was projected onto the x-axis, given that the main objective is to have a formula for the droplet travel distance x.To overcome the second difficulty, a linearization was introduced: + η (41) The coefficients λ and η will be determined using the regression method (Appendix A), as expected for the previous coefficients α, ε and φ.Equation (40) transforms as follows: The variables x d , y d and s d originate at the top of the trajectory (Figure 1).The solution of Equation ( 42), including the constant of integration, is For Equation (43) to provide the droplet travel distance x d = x end up to the point of impact to the ground, the landing angle of the trajectory at the point of impact θ end must be known (Figure 1).An effective way to solve the problem is to start from the observation that the tan θ coincides with the derivative dy d /dx d .Therefore, it is sufficient to highlight the tan θ and proceed to the integration to obtain the coordinate y d = y d (x d ).
To facilitate the operation, some groups of parameters will be collected under two constant symbols, K 1 and K 2 : Thus, Equation (43) becomes From which it is easy to highlight tan θ: Integration of Equation (47) gives The combination of Equations ( 46) and (48) gives If this Equation ( 49) is applied to the landing point of the droplet, then the total height of the descending trajectory is known and equals y end = y top + h (Figure 1), where y top is obtained from Equation (25).Therefore, from (49), the only remaining unknown can be found: θ = θ end .
With the value of θ end entered into Equation (46), the complete travel distance of the descending portion x end can be obtained.The total travel distance is the sum of the travel distance of the ascending portion x top (25) and the descending portion x end (46): Concerning the flight time during the descending trajectory, Equation (39) needs to be integrated.However, the exponent 1 δ−1 is not an integer and, therefore, prevents a solution.As already performed to solve Equation (40), a linearization is necessary: + ψ (51) The coefficients σ and ψ will be determined using the regression method (Appendix A), as expected for the previous coefficients λ, η, α, ε and φ.Equation (39) transforms as follows: After integrating (52) and entering the value of the landing angle θ = θ end , the result obtained is the flight time from the top to the point of landing.
Finally, the total flight time t tot results from the sum of t top (26) with t end (53):

Numerical Solution
The approximate closed-form solution, which led to the analytical modeling of the trajectory and, therefore, the travel distance, represented by Equations (18) for the calculation of Re top , Equations ( 24) and ( 25) for the calculation of x top and y top , Equation (49) for the calculation of θ end and Equation (46) for the calculation of x end , was evaluated by comparing its results with the exact ones obtained with a numerical solution of the differential equation (ODE) of motion: This Equation ( 55) is that of the ascending portion of the trajectory and is different from Equation ( 7) because the variable exponent with ln Re has been kept, and therefore, it is the exact ODE.A similar differential equation was used for the descending portion of the trajectory.It also differs from Equation ( 28) because the variable exponent with ln Re has been maintained: Both ODEs were integrated numerically for the first time to obtain Re as a function of θ.Finally, a subsequent numerical integration led to travel distance x and flight time t.
Numerical Adams' method, especially suitable for this type of differential equation and easily implementable in Excel, was used for solving the two ODEs (55) and (56).

Comparison between Approximate Closed-Form Solution and Exact Numerical Solution
To evaluate the accuracy of the analytical modeling constituted by the approximate closed-form solution obtained in Section 2.1, the relevant equations to determine the travel distance were applied to droplets of seven different diameters (D = 0.25, 0.5, 1, 2, 3, 4 and With the same values of D, θ 0 and p, the numerical integrations described in Section 2.2 were carried out to obtain the exact values of the travel distance x tot .
The comparison results are shown in Figure 4a-d, respectively, for each of the four jet angles θ 0 but with the same pressure p = 2.5 bar.To verify the accuracy of the approximate closed-form solution compared to the exact numerical one, with p = 4 bar, the numerical experimentation was limited to the variation of the diameter D, maintaining a single jet angle of 30 • .The comparison results are shown in Figure 5.
Both ODEs were integrated numerically for the first time to obtain Re as a function of θ.Finally, a subsequent numerical integration led to travel distance x and flight time t.
Numerical Adams' method, especially suitable for this type of differential equation and easily implementable in Excel, was used for solving the two ODEs (55) and (56).

Comparison between Approximate Closed-Form Solution and Exact Numerical Solution
To evaluate the accuracy of the analytical modeling constituted by the approximate closed-form solution obtained in Section 2.1, the relevant equations to determine the travel distance were applied to droplets of seven different diameters (D = 0.25, 0.5, 1, 2, 3, 4 and 5 mm) produced by a sprinkler with different jet angles (θ0 = 40°, 30°, 20° and 10°) and with two different pressures (p = 2.5 and 4 bar).The reference environmental conditions were air temperature Tair = 20 °C and altitude 50 m, which corresponded to air density ρair of 1.2 kg m −3 and air viscosity µ of 0.000018 Pa s.
With the same values of D, θ0 and p, the numerical integrations described in Section 2.2 were carried out to obtain the exact values of the travel distance xtot.
The comparison results are shown in Figure 4a-d, respectively, for each of the four jet angles θ0 but with the same pressure p = 2.5 bar.To verify the accuracy of the approximate closed-form solution compared to the exact numerical one, with p = 4 bar, the numerical experimentation was limited to the variation of the diameter D, maintaining a single jet angle of 30°.The comparison results are shown in Figure 5.    • 100 (57) Considering the set of 35 error% values, the mean relative error MRE is 1.8%, and the standard deviation SD is 1.4%.These values show how analytical modeling, i.e., based on the approximate closed-form solution, is a valid calculation tool.
To have a further qualitative evaluation of the accuracy of the approximate analytical solution compared to the exact numerical one, in Figure 6a-d, the trajectories were plotted and calculated both analytically and numerically with the four different jet angles (40°, 30°, 20° and 10°).The curves are limited to the case of the 2 mm droplet because, for other droplet diameters, the differences between the two trajectories are similar.These differences between the two trajectories, one calculated analytically (black line) and one calculated numerically (red line), confirm the relative error% calculated on the total travel distance, which was shown in Figures 4 and 5.
For diameters less than approximately 0.45 mm, it may happen that the constant K2 calculated with Equation (45) is negative, just below zero, due to the approximation of the values of the linearization coefficients introduced with the regression (Appendix A) of λ and η, despite R 2 = 0.99.This makes the solution of Equation ( 49) impossible because the argument of the logarithm function becomes negative.
The meaning of this event is that K2 should tend to 0 due to the small size of the droplet and, therefore, its tendency to end its flight in a direction close to the vertical ( → 90°), i.e., mathematically, the trajectory tends toward an asymptote.Therefore, for very small variations of θend, the variations of xend are infinitesimal, but those of yend become very large.This extreme situation would require absolute accuracy in the calculation of  → 0, which the regression of λ and η cannot ensure.However, if (45) gives a negative K2 just below zero, the difficulty is easily overcome by replacing this negative value with a positive K2 value close to 0. Empirically, it has been found that the optimal value is K2 = 10 −4 .The histograms of all the Figures show, through horizontal bars, the values of the travel distance x tot obtained analytically (blue color) and numerically (red color) for each of the seven diameters of the droplet.In the histograms, for each diameter, there is a third bar (green color), which represents the relative percentage error, equal to Considering the set of 35 error% values, the mean relative error MRE is 1.8%, and the standard deviation SD is 1.4%.These values show how analytical modeling, i.e., based on the approximate closed-form solution, is a valid calculation tool.
To have a further qualitative evaluation of the accuracy of the approximate analytical solution compared to the exact numerical one, in Figure 6a-d, the trajectories were plotted and calculated both analytically and numerically with the four different jet angles (40 • , 30 • , 20 • and 10 • ).The curves are limited to the case of the 2 mm droplet because, for other droplet diameters, the differences between the two trajectories are similar.These differences between the two trajectories, one calculated analytically (black line) and one calculated numerically (red line), confirm the relative error% calculated on the total travel distance, which was shown in Figures 4 and 5.
For diameters less than approximately 0.45 mm, it may happen that the constant K 2 calculated with Equation (45) is negative, just below zero, due to the approximation of the values of the linearization coefficients introduced with the regression (Appendix A) of λ and η, despite R 2 = 0.99.This makes the solution of Equation (49) impossible because the argument of the logarithm function becomes negative.
The meaning of this event is that K 2 should tend to 0 due to the small size of the droplet and, therefore, its tendency to end its flight in a direction close to the vertical ( θ end → 90 • ), i.e., mathematically, the trajectory tends toward an asymptote.Therefore, for very small variations of θ end , the variations of x end are infinitesimal, but those of y end become very large.This extreme situation would require absolute accuracy in the calculation of K 2 → 0 , which the regression of λ and η cannot ensure.However, if (45) gives a negative K 2 just below zero, the difficulty is easily overcome by replacing this negative value with a positive K 2 value close to 0. Empirically, it has been found that the optimal value is K 2 = 10 −4 .

Comparison between Approximate Closed-Form Solution and Experimental Data
Analytical modeling based on the approximate closed-form solution was also applied to the experimentally measured droplet motion by Tompson et al. [27,28].These authors measured the travel distance xtot and the flight time ttot for droplets with diameters D equal to 0.3, 0.9, 1.8, 3.0 and 5.1 mm emitted by a sprinkler with a jet angle θ0 of 25° and with an initial velocity v0 of 30.9 m/s.Nozzle height above ground h was 4.5 m and environmental conditions were air temperature Tair = 38 °C, relative humidity RH = 20% and

Comparison between Approximate Closed-Form Solution and Experimental Data
Analytical modeling based on the approximate closed-form solution was also applied to the experimentally measured droplet motion by Tompson et al. [27,28].These authors measured the travel distance x tot and the flight time t tot for droplets with diameters D equal to 0.3, 0.9, 1.8, 3.0 and 5.1 mm emitted by a sprinkler with a jet angle θ 0 of 25 • and with an initial velocity v 0 of 30.9 m/s.Nozzle height above ground h was 4.5 m and environmental conditions were air temperature T air = 38 • C, relative humidity RH = 20% and no wind.Tompson's experimental results compared to those of this study are shown in Tables 1 and 2.About the travel distance (Table 1), the mean relative error MRE is 17%, mainly due to a prediction of excessive travel distance x tot for large droplets, from 3 to 5.1 mm.
For the flight time (Table 2), an MRE of 28% occurs due to the excessive time t tot predicted by the analytical modeling for the larger droplets and a reduced t tot for the 0.3 mm droplet.

Modification of the Drag Coefficient Formula to Consider the Deformation of the Droplets
Tables 1 and 2 showed, in the case of medium and large droplets, too high values of the total travel distance and total flight time predicted by the analytical modeling.This is attributable to the reduced value of the drag coefficient obtained from Equation (4).In fact, it was obtained with a regression procedure starting from experimental data of the Standard Drag Curve (SDC) [25,26], which are valid for rigid spheres.
Yan et al. [21] and Kincaid [29], however, suggest using the approach of Park et al. [30,31], who, analyzing the data of Laws [32] and Gunn and Kinzer [33], found that the drag coefficient c d for Re ≥ 1000 begins to increase compared to the constant value of 0.44 due to the deformation of the droplets during flight.
The equation that Park et al. [30,31] suggested to take this fact into account is the following: In this formula, there is a first constant term equal to 0.438, which is like the value of 0.44 predicted by Bird's model [24] and 0.45 predicted by Fukui's model [14].Then, there is a second term in parentheses, which is a factor greater than 1 for Re ≥ 1000.Due to this increasing factor, for example, with Re = 10, 000, the drag coefficient c d is 0.591.
In Table 3, the c d values from Equations ( 4) and (58) have been compared considering the Reynolds numbers Re 0 of the experimental tests of Tompson et al. [27,28].Table 3 shows that the values of Equation ( 4) are close to the experimental ones of the Standard Drag Curve-SDC [25,26], while those of Equation ( 58) are increasing with Re 0 .The solution to this problem is to introduce the increasing factor present in the Equation (58) of Park et al. [30,31], also in Equation ( 4).
But Equation (58) has the term inside the brackets, which for Re ≤ 1000 cannot be calculated because it contains a negative number with a non-integer exponent, and in any case, the term in the brackets should remain constant and equal to 1. Therefore, a new empirical formula has been identified, which has a polynomial in place of the Park Equation (58), so that, for Re ≤ 1000, Equation (4) remains valid, and for Re > 1000, it provides increasing values like those of Park.Ultimately, Formula (4), also considering the hypotheses written for Equation (7), becomes The new c d values calculated with Equation (59) are, in the order of increasing D from 0.9 to 5.1 mm, 0.450, 0.456, 0.502 and 0.560.They differ from those obtained with Equation (58) of Park et al. (Table 3), with an MRE of only 1.5%.It must be noted that (58) is valid only for 1000 ≤ Re ≤ 10, 000, while (59) is valid for 0.3 ≤ Re ≤ 10, 000.
Formula (59) keeps expression (4) within it; therefore, the only modification to be made to the analytical modeling of Section 2.1.is to replace, in all the equations in which it is present, the number 25.1 with the constant K equal to Because of this modification, the new values of total travel distance x tot and total flight time t tot are shown in Tables 4 and 5, respectively.Table 4 shows the new travel distance values x tot and the relative error% compared to Tompson's experimental data.They are significantly better than those in Table 1.Now, the mean relative error MRE is 2.5%, and the standard deviation SD is 1.1%.Therefore, the new values of c d calculated with Equation (59), based on the indications of Park et al. [30,31], and the consequent replacement in the analytical model of Section 2.1.of the number 25.1 with the constant K of Equation ( 60), made the model capable of correctly describing the trajectory and therefore the total travel distance x tot of the droplets experimentally analyzed by Tompson et al. [27,28].
Table 5 shows the new flight time values t top .Excluding the flight time of the 0.3 mm droplet, the others were improved compared to those in Table 2.In particular, the mean relative error-MRE was 9.8%, and the SD was 3.2%.
The error of −60.9%, relating to the smallest droplet, seems to have no explanation.Perhaps the flight time was longer during the experimentation because the droplet was also involved, in the last part of the almost vertical trajectory, in a convective upward motion of the air due to temperature differences with the ground; very small mass of the droplet (14 µg); high height of the nozzle from the ground (4.5 m); final angle calculated with Equation (49) very high (θ 0 = 89.8• ), which confirms a last part of the trajectory that is almost vertical.In this section, only the final equations of the closed-form integration carried out in Section 2.1 are collected, together with Equation (60) of Section 3.3 and the regression equations of Appendix A. They represent the mathematical model to be used on the spreadsheet or to be inserted into the PLC software (microprocessors dedicated to the control and regulation of sprinklers) to determine the geometric and dynamic values of the trajectories useful for the study and mitigation of the phenomena of soil sealing and erosion.

Synoptic Table of the
The travel distance of the ascending portion of the trajectory (0-top) is (Figure 1): where where the quantities K, Ar, β, α and C ua are where the initial value of Reynolds number (Figure 1) is Re 0 = ρ air v 0 D µ .The height from the nozzle to the top of the trajectory y top is (Figure 1) where The number of Reynolds at the top of trajectory Re top is The flight time from the nozzle to the top of the trajectory t top is where φ = 1.247•(cos θ 0 ) −0.5865 •Re −0.02691 0 (A2) The travel distance of the descending portion of the trajectory (top-end) is (Figure 1): where the quantities K 1 , K 2 , δ, ε, λ, η and C ud , are The tan θ end can be obtained with an iterative method from the following: The terminal velocity v end is The total travel distance x tot results from the sum of x top (25a) and x end (46): The flight time from top to end of the trajectory t end is where The total flight time t tot is the sum of t top (26) and t end (51): The nomenclature of the quantities present in this Section 3.4 is collected in Table 6.

Results of the Application of Definitive Analytical Modeling
The definitive approximate analytical modeling, presented in the previous Section 3.4, was applied with the same diameters D (0.25, 0.5, 1, 2, 3, 4 and 5 mm) and with the same jet angles θ 0 (40 • , 30 • , 20 • and 10 • ) and with the pressure p of 2.5 bar adopted during the comparison (Section 3.1) between the exact numerical solution and the approximate analytical one.For each diameter and for each jet angle, the result was the trajectories diagram (Figure 7a-d) and the total travel distance x tot , the maximum height y top , the impact angle on the ground θ end and the terminal velocity v end (Table 7).
Figure 7a-d, respectively, for the jet angles of 40 • , 30 • , 20 • and 10 • , each show the trajectories of the seven droplets differing in diameter, from 0.25 to 5 mm.Note that the total travel distances are smaller than those reported in Figure 4a-d due to the correction introduced with Equation (60), which is necessary to take into account the deformation of the larger droplets during the motion as predicted by Park et al. [30,31].
Finally, Table 7 shows the values of the total travel distance x tot , obtained with Equation (50); of the maximum height y top , obtained with (25d); of the landing angle θ end , obtained with (49) and of the terminal velocity v end , obtained with (36).
In the application of the definitive analytical model presented in Section 3.4, Figure 7 and, even more so, Table 7 show how the smallest jet angle (10 • ) has the smallest landing angles and quite high landing speeds.Scholars [1][2][3][4] point out that these values of the two parameters, the high dynamic one (landing velocity) and the low geometric one (landing angle), are the most dangerous for erosion and the formation of the superficial crust with relative reduction in infiltration.The travel distance achieved with a jet angle of 10 • is also the lowest, therefore showing a lower efficiency in the use of irrigated water.
The largest landing angles, at less risk of erosion and the formation of surface crust, are obtained with jet angles of 40 • .However, accepting a slight reduction in the landing angle but with an improvement in the travel distance and, above all, with a reduction in the landing velocity and, therefore, in the impact kinetic energy, the jet angle of 30 • seems the most appropriate.In any case, the application carried out here is limited to discrete jet angle values and, therefore, has only an example value.The reader will be able to carry out the application with different angles if they are within the 10-40 • range.
As regards the influence of the diameter, Table 7 shows how the smaller droplets have much less impact and, therefore, are less dangerous for erosion and the formation of surface crust, both in terms of kinetic energy and landing angle.

Analytical Modeling and Droplet Evaporation
The analytical modeling summarized in Section 3.4 has, among the hypotheses, zero evaporation of the droplets during the flight.How far is this hypothesis from reality?
The results shown in Table 4, where the total travel distance calculated with the analytical modeling without evaporation are compared to the experimental one of Tompson et al. [27,28], where the droplets evaporated, highlights that the zero-evaporation hypothesis of the modeling is acceptable.
Confirmation can also be found in the experimentation and mathematical modeling of Kincaid and Longley [34].They showed the influence of the droplet diameter D and environmental conditions (dry bulb temperature T air and relative humidity RH%) on the loss rate of water mass through evaporation [35] during the flight.
The loss rate increases with increasing temperature, decreasing humidity and droplet diameter.The situation with the greatest loss rate occurred with the smallest droplet The analytical modeling was built for a range of jet angle values between 10 • and 40 • ; therefore, outside of it, the use of the model could lead to important errors.
Ultimately, the proposed analytical modeling is reliable and can predict, in the absence of wind, the travel distance, the maximum height of the trajectory, the landing angle and the terminal velocity before impact on the ground.This is all useful information for controlling water waste and the erosive effects on the soil due to the kinetic energy and landing angle of the droplets.
Furthermore, the proposed analytical modeling contains only elementary functions such as trigonometric ones, the logarithm and the exponential, which allow it to easily become the software for PLC control systems to be applied to sprinklers.In the future, it will, therefore, be interesting to verify the application to sprinkler irrigation systems of PLC-based hardware implemented with software attributable to this analytical modeling, which should also be expanded to consider the effects of the wind.

Figure 1 .Figure 1 .
Figure 1.Flight trajectory of a droplet emitted by a sprinkler at a height h above the ground.The trajectory is composed of a first portion ascending to the top and a second portion descending from the top to the landing.

Figure 2
Figure2.The local reference system adopted, with the s and n axes that are, respectively, tangent and normal to the trajectory.The red vector is the drag force  , the black vector is the gravitational force ̅ projected in the two components along s ( ̅ sin ) and n ( ̅ cos ).The blue vector is the velocity ̅ .

Figure 2
Figure2.The local reference system adopted, with the s and n axes that are, respectively, tangent and normal to the trajectory.The red vector is the drag force R, the black vector is the gravitational force mg projected in the two components along s (m D gsin θ) and n (m D gcos θ).The blue vector is the velocity v.

Figure 3 .
Figure 3. (a) In an infinitesimal interval of time dt, the droplet travels along an elementary trajectory section ds, which subtends an angle dθ and the velocity vector ̅ becomes ̅ + ̅ .(b) Graphic construction to highlight the components along s and n of the velocity variations and, therefore, of the accelerations.

Figure 3 .
Figure 3. (a) In an infinitesimal interval of time dt, the droplet travels along an elementary trajectory section ds, which subtends an angle dθ and the velocity vector v becomes v + dv.(b) Graphic construction to highlight the components along s and n of the velocity variations and, therefore, of the accelerations.
5 mm) produced by a sprinkler with different jet angles (θ 0 = 40 • , 30 • , 20 • and 10 • ) and with two different pressures (p = 2.5 and 4 bar).The reference environmental conditions were air temperature T air = 20 • C and altitude 50 m, which corresponded to air density ρ air of 1.2 kg m −3 and air viscosity µ of 0.000018 Pa s.

Figure 5 .
Figure 5.Total travel distance xtot obtained with closed-form solution (blue bar) and with numerical solution (red bar) and relative error% (green bar) vs. droplet diameter D with nozzle pressure p = 4 bar and jet angle θ0 = 30°.The histograms of all the Figures show, through horizontal bars, the values of the travel distance xtot obtained analytically (blue color) and numerically (red color) for each of the seven diameters of the droplet.In the histograms, for each diameter, there is a third bar (green color), which represents the relative percentage error, equal to Error% = ( ) − ( ) ( )• 100 (57)

Figure 5 .
Figure 5.Total travel distance x tot obtained with closed-form solution (blue bar) and with numerical solution (red bar) and relative error% (green bar) vs. droplet diameter D with nozzle pressure p = 4 bar and jet angle θ 0 = 30 • .

Figure 6 .
Figure 6.Trajectories of the droplet with a diameter D of 2 mm and with different jet angle values θ0 of 40° (a), 30° (b), 20° (c) and 10° (d): black line for the closed-form solution; red line for numerical solution.Pressure p of 2.5 bar.

Figure 6 .
Figure 6.Trajectories of the droplet with a diameter D of 2 mm and with different jet angle values θ 0 of 40 • (a), 30 • (b), 20 • (c) and 10 • (d): black line for the closed-form solution; red line for numerical solution.Pressure p of 2.5 bar.

Figure 7 .
Figure 7. Trajectories of the droplet with the droplet diameters D of 0.25, 0.5, 1, 2, 3, 4 and 5 mm and with different jet angle values θ0 of 40° (a), 30° (b), 20° (c) and 10° (d); for each figure, the trajectories increase the travel distance as the diameter increases from 0.25 up to 5 mm.Pressure p of 2.5 bar.

Figure
Figure7a-d, respectively, for the jet angles of 40°, 30°, 20° and 10°, each show the trajectories of the seven droplets differing in diameter, from 0.25 to 5 mm.Note that the total travel distances are smaller than those reported in Figure4a-d due to the correction introduced with Equation (60), which is necessary to take into account the deformation of the larger droplets during the motion as predicted by Park et al.[30,31].

Figure 7 .
Figure 7. Trajectories of the droplet with the droplet diameters D of 0.25, 0.5, 1, 2, 3, 4 and 5 mm and with different jet angle values θ 0 of 40 • (a), 30 • (b), 20 • (c) and 10 • (d); for each figure, the trajectories increase the travel distance as the diameter increases from 0.25 up to 5 mm.Pressure p of 2.5 bar.

Table 2 .
[27,28]light time t tot of sprinkler droplets: data of this study compared to those of Tompson et al.[27,28].

Table 4 .
[27,28]ravel distance x tot of the droplets: data from the analytical modeling completed with Equation (60) of this study, compared to the experimental ones of Tompson et al.[27,28].

Table 6 .
Nomenclature relating to the quantities of the equations in the synoptic table of definitive analytical modeling.

Table 7 .
Geometric and kinematic characteristics of droplet trajectories: total travel distance x tot , maximum height from ground y top + h, landing angle θ end and terminal velocity at the landing v end .Pressure p of 2.5 bar.