Mathematical Modelling of the Impact of Liquid Properties on Droplet Size from Flat Fan Nozzles Problem presented

Flat fan nozzles atomize crop protection products, breaking them into droplets. Droplet size matters - smaller droplets give better perfor- mance, but very small droplets drift. We want to use mathematical models to better understand how liquid properties affect droplet size. 
There are three types of breakup: wavy sheet, perforation, and rim. In wavy sheet breakup, increasing viscosity or surface tension increases droplet size. To investigate further, we carry out direct numerical simulations of jet breakup, which show that suface tension has little effect, but increasing viscosity leads to fewer droplets. Decreasing the jet velocity also results in fewer droplets, with a wider size distribution. 
Each type of breakup involves primary breakup into cylinders of fluid, then secondary breakup into droplets. We thus consider the breakup of a cylinder of fluid. Direct numerical simulations suggest that within the tested parameter range viscosity has little impact on droplet size, however it does influence the timescale on which the instability evolves considerably. Linear stability analysis suggests that increasing viscosity increases the wavelength of the most unstable mode, which we expect leads to larger droplets, and that it reduces the rate of breakup. 
Perforations - holes in the sheet - also lead to breakup. We find how the length fraction of the sheet that is void changes with time. 
After breakup, the droplets continue to evolve. We develop a model, based on a transport equation, for this process. A key parameter is the breakup rate constant - larger values lead to more breakup, fewer large droplets, and a narrower size distribution. 
Together, these mathematical approaches improve our understanding of how droplets form, and can be used to guide experimental work.

1 Introduction (1.1) Flat fan nozzles are used to atomise crop protection products, and their performance is affected by droplet size.Knoche [14] claims that small droplets are good for herbicide performance.If droplets are too small, however, they drift.We want to understand how the liquid's physical properties (viscosity, surface tension, etc.) affect droplet size.
(1.2) Droplets form when liquid exits the nozzle as a sheet and then breaks up via three different breakup mechanisms: wavy sheet, perforation, and rim.Rim breakup produces the largest droplets, followed by perforation, then wavy sheet instability.The objective of the Study Group was to investigate what causes these breakup mechanisms, and how they are affected by the liquid's physical properties (viscosity, surface tension, composition, etc.) (1.3)In this report, we present relevant results from the literature in section 2.
Then, section 3 is on direct numerical simulations.Section 4 is on the breakup of liquid ligaments after the sheet has ruptured, section 5 is on perforation, and section 6 is on the late stage droplet distribution.Finally, section 7 contains our conclusions and suggestions for future work.
2 Literature review 2. 1 The effects of liquid properties on droplet formation (2.1.1)As expected, results in the literature show that increasing viscosity increases droplet size.Gratkowski and Stewart [12] claim that "[thickening agents] increase the viscosity of the water phase of spray solutions, thus increasing the size of droplets during spray application".Different viscosities lead to different flows, as shown in Fig. 1, and to different types of breakup.Li and Ashgriz [17] show empirically that the type of breakup of a liquid jet is dependent on a dimensionless number called the Reynolds number, where ρ is the liquid's density, U 0 is the characteristic velocity scale, L is the characteristic length scale, and µ is the liquid's dynamic viscosity (notation is also defined in Appendix F).Ahmed et al [2] show similar results to [17], but with viscosity changing independently, rather than as part of Re.While the large increase in viscosity in [2] -a factor of 80may be infeasible for Syngenta's application, we expect viscosity to have an effect.
(2.1.2)Miller and Ellis show that an increase in surface tension increases droplet size [20], but they note that it is unclear whether this increase in droplet size is due solely to surface tension.They also consider the composition of the liquid, finding that a more homogeneous liquid gives more controlled Figure 3: Different Reynolds numbers give different breakup regimes, from [17] results for the distribution of the liquid, as shown in Fig. 2.

Rim breakup
(2.2.1) Changing Re results in different types of breakup.Specifically, Fig. 3 shows that, for Re <∼ 1000, a sheet has yet to form (instead the jet breaks up due to a capillary instability).For Re ∼2000-2400, a smooth (closed) sheet forms which exhibits a Plateau-Rayleigh type rim instability, but without waves.For higher Re, ∼2500-3000, an open sheet forms, which breaks up at the bottom, showing wavy sheet instability.
For Re > 3000 the flow becomes turbulent and unstable.
(2.2.2) We are ideally want an open sheet with minimal rim breakup (since rim breakup results in the largest droplet sizes).We thus want to adjust the liquid properties such that the corresponding value of Re is achievedfor example, we could increase the viscosity to a degree where the rim is stable, but the droplet size is small enough.

Wavy sheet instability droplet size
(2.3.1)In the wavy sheet instability, waves form on the liquid sheet, then liga-ments of liquid separate off, breaking down further into droplets.(We present a simple stability analysis for liquid sheet breakup, exploring the stability of waves on the sheet, in appendix A.) By looking for the wave with the maximum growth rate, assuming that ligaments of half this wavelength break off, and that these ligaments break down by the Rayleigh-Plateau-Savart instability (see section 4), Fraser et al. [8] (neglecting viscosity) and Dombrowski and Johns [6] (including viscosity) find expressions for droplet sizes.The latter is where d D is the droplet diameter, µ is the liquid's dynamic viscosity, ρ is the liquid's density, σ is the surface tension coefficient, and d L is the ligament diameter, where K is the liquid sheet thickness x distance from the source, ρ a is the density of the surrounding air, and U is the sheet velocity.
(2.3.2) In equation ( 2) low surface tension leads to small droplets, and high viscosity leads to large droplets, although as noted in [6], "the effect of viscosity on drop size is dependent on the other operating conditions, being greater for liquids of low density and surface tension".Here, viscosity appears as part of the Ohnesorge number µ/ √ ρσL, a dimensionless number that characterises the relative importance of viscous and surface tension forces [32], where L is a characteristic length scale -here it is the ligament diameter.The relative importance of inertial and surface tension forces is characterised by the Weber number [33], ρU 2 0 L/σ, where U 0 is the characteristic velocity scaling.Varying these dimensionless numbers independently, as was discussed for Reynolds number in section 2.2, provides insight into the different breakup mechanisms.

Introduction
(3.1.1)The dynamics of liquid sheet breakup emerging in agricultural sprays as a result of using flat fan nozzles is a highly complex and nonlinear process.
The typical velocities of at least 10 m/s, as well as the large range of parameters to be varied make this problem particularly challenging.
(3.1.2)Experimental, analytical, and numerical approaches are all useful in providing insight into drop formation.Experimental methods give a direct line of investigation into the target flow properties, but are often restricted by laboratory time, equipment, video technology limitations and the sheer size of the parameter space.Analytical techniques offer a fundamental understanding, producing useful formulae, which can then drive innovation further in the design pipeline.They are sometimes, however, restricted to the early stages of the flow, when the processes are linear or weakly nonlinear and will make simplifying assumptions to make the problem tractable.Direct numerical simulations offer a bridge between the two.They are "direct" because the full Navier-Stokes equations are solved, taking into account both air and liquid properties, as well as the full set of interfacial conditions.In the early stages, numerical results can cross-validate analytical progress.In the late stages, accurate numerical experiments delve into nonlinear regimes beyond the reach of other techniques, offering detailed information at any timestep and from angles difficult to construct experimentally.Although they need significant computational resources (processing power, memory, data handling algorithms), they serve to reduce the search space, guiding experimental investigations in the right direction.
(3.1.3)We use direct numerical simulations (the numerical methodology is explained in Appendix B.1) to investigate two related problems: 1.The primary breakup of the liquid sheet (wavy sheet instability) into liquid cylinders.
2. The secondary breakup of the liquid cylinders into droplets through the Rayleigh instability.

Jet atomisation in two dimensions
(3.2.1)In order to better understand the target regime and because the liquid sheet is thin, we consider the breakup problem in two dimensions: a cross-section through the main body of the sheet, excluding any rim effects.Liquid fragments break off from the main jet as columns, which subsequently breakup into droplets, as discussed in section 3.3.Here the primary focus is the wavy sheet instability and the dynamics in the main body of the jet. (3.2. 2) The domain and an example snapshot are shown in Fig. 4, which also outlines the size of the inlet region (around 1 mm) and the velocity of the jet, taken to be 16 m/s, in line with typical operational conditions.Outflow conditions are prescribed on either side of the domain, as well as at the bottom.We turn our attention towards gathering and analysing information about the formation of droplets, their number and their size distribution during the breakup, as in Fig. 5.It illustrates the tracking of O(1000) droplets, which, at a given time step far enough into the flow evolution, have a size characterised by a log-normal distribution.Similarly, we target variations in the parameters of interest: Figure 5: Example quantitative data obtained from a typical direct numerical simulation.Top: number of droplets in the computational domain in time, with a specific timestep selected in order to study the drop size distribution (bottom).
• Jet velocity 4 m/s < U < 16 m/s • Surface tension coefficient 0.025 N/m < σ < 0.072 N/m.The upper bound represents water, while the lower bound is given by the declared threshold given any alterations of the liquid formula (for example through surfactants).
The velocity is fixed at 16 m/s in all cases.
• Viscosity 0.001 kg/m•s < µ < 0.016 kg/m•s, with the same rationale as above and using the values declared as being of interest.
Fig. 6-Fig.8 summarise the results, indicating the number of droplets formed as a function of time in each of the three studies, with three subcases each.Distributions of the drop sizes at selected timesteps are provided in order to complement the raw numerical data, with bar colours on the right of each figure corresponding to the counterpart line plots on the left.
Lowering the velocity of the jet (Fig. 6) reduces the number of droplets formed, and impacts their distribution.In the largest velocity case, a relatively clean spread of the droplets between 100 µm and 500 µm is observed, which becomes progressively flatter as the jet velocity is lowered, with three times fewer droplets when the jet speed is set to 4 m/s.In this case droplets larger than 500 µm in diameter are also clearly represented, which is less pronounced at higher velocities.Interestingly, the variation in surface tension (Fig. 7) has little impact on the number of droplets (with the current grid resolution at least), but it does make the size distributions cleaner and shift slightly towards the smaller size.This may in part be due to the fact that the window of variation is quite small and the regime even for the highest surface tension coefficient chosen is already quite violent for the instability targeted by these computations.Surface tension may of course affect the rest of the dynamics (perforations, rim instabilities, Rayleigh breakup) in a different manner.
As viscosity is increased (Fig. 8) we reach a configuration with one large liquid volume (the body of the jet), out of which is the occasional pinch-off, with fewer droplets as this parameter is increased.We thus conclude that viscosity is one of the primary factors affecting the initial wavy sheet instability.(3.3.

The Rayleigh instability
2) The stability of an axisymmetric column of liquid is one of the most stringent tests for computational fluid dynamics packages, with the multi-scale nature of the flow and the interfacial rupture posing numerical challenges that are often tested to ensure robustness of implementation.For the package we use, Gerris, a discussion on this topic and extensive validation with the results of Rayleigh and Weber have been described by Popinet [22].The key parameter in the flow is the so-called Laplace number with σ denoting the surface tension coefficient, L the reference lengthscale, ρ the liquid density and ν = µ/ρ its kinematic viscosity.We consider four different cases of µ ranging from 10 −3 Pa•s to 8 • 10 −3 Pa•s.We note that much higher viscosities are possible in the target application. (3.3. 3) The remaining parameters are for the classical water-air case, with σ = 0.072 N/m, ρ = 998 kg/m 3 and L = 50 µm chosen as a reference lengthscale.Two snapshots of the flow evolution are presented in Fig. 9, with the plot at the later time step suggesting the formation of two droplets of different sizes, with radii that are tracked over time.To ensure a fast runtime for the full three-dimensional implementation, we have used all available symmetries and are solving over one eighth of the domain.Each direct numerical simulation ran for 8 hours on high performance computing clusters; the adaptive mesh refinement (carefully crafted around curvature criteria for the pinch off) allows us to run these on a reasonable timescale.
(a) Evolution in time of the normalised radius of the largest drop.
(b) Comparison to experimental results [16] in the inviscid case.(3.3.4)Fig. 10 provides an outline of the results, starting with a comparison of the least viscous case with the experimental work of Lafrance [16], as well as previous analytical results on the right hand side.We find excellent agreement between both the large and the small drop radii for the most unstable wavenumber in this setting, giving us confidence moving forward with the parameter variation.First, we note that the range of viscosities chosen is well within the allowed operational envelope, with the largest one limited by the timescale of the instability.A much larger liquid viscosity results in a very slow development of the flow.There are two noticeable features of the results in Fig. 10a: • Viscosity has virtually no impact on the final drop size, with the radius changing by less than 1% from the least to the most viscous case, albeit in the right direction (larger viscosity leads to larger radius).
• Increasing the viscosity damps the oscillation considerably.An idealised inviscid case (zero viscosity) would oscillate indefinitely without any damping, however the µ = 10 −3 Pa•s case is a low viscosity, allowing it to be compared to theory.
Since the results of the previous subsection indicated significant variation of the ligament formation properties with respect to changes in viscosity, it is interesting to note that the secondary breakup appears to be much less sensitive, at least in terms of final drop sizes.We do underline however that this is in an ideal setting in which the air flow surrounding the liquid column is stationary.

Summary
(3.4.5)The main outcomes of this section are: • Development of a versatile implementation capable of qualitative and quantitative study of jet breakup and drop dynamics.
• Variation of velocity, density, viscosity, surface tension within the desired operational envelope is possible, and several numerical experiments on the wavy sheet instability have been conducted in two dimensions.
• Lowering the velocity of the jet decreases the number of droplets from the primary breakup and significantly widens their size distribution.
• Increasing viscosity strongly reduces primary breakup in that the liquid sheet remains almost intact and very few liquid columns form.
• Varying the surface tension coefficient within the relatively small given range produced a negligible effect on the primary wavy sheet instability.
• The secondary instability (breakup of liquid columns into droplets) has been studied in a full three-dimensional context in the absence of an external flow.
• No significant effect of viscosity on the drop radii has been found, however the evolution to a steady state is strongly affected, with viscosity strongly damping the oscillations following the interfacial rupture and relaxation of the drop shapes.
Many of the results and their subsequent post-processing and analysis have been adjusted to account for the time window of the study group, in which tangible progress was obtained in the order of a few days.On longer timescales, several exciting research branches present themselves on both fundamental and practical sides of the project, and they are outlined in section 7.1.
4 Linear stability of an axisymmetric column of liquid  The three breakup mechanisms (wavy sheet, perforation, rim) all result in a cylinder of liquid breaking up into droplets.In wavy sheet breakup, ligaments of the sheet break off and recoil into cylinders before disintegrating into droplets [6].In rim breakup, surface tension forces at the edge of the sheet lead to a cylindrical rim, which then breaks into droplets.For perforations, the growing holes slowly merge, forming a web of filaments that breaks up into individual ligaments.Thus, we revisit the classical studies of the breakup of a cylinder of fluid that undergoes small perturbations to its steady state.
(4.1.2)Rayleigh [24] demonstrated that a cylinder of inviscid fluid breaks up due to surface tension if the cylinder is slightly perturbed (the Rayleigh-Plateau-Savart instability).We begin by considering the more general case of two viscous fluids: a cylinder of fluid 1 in an infinite pool of fluid 2 (think water-in-air), and derive the linear stability problem following Tomotika [27].We then neglect the air at leading-order, since the viscosity and density ratios are so small, leading to Weber's findings [31].
(4.1.3)We aim, in this section, to investigate how changing the viscosity changes the wavelength of the most unstable mode.Droplet sizes could then be estimated for a given initial cylinder radius.This analysis is also useful for comparison with the numerical simulations of section 3. short timescales mean its effect is small.The analysis for this general case is given in Appendix C.1.

Stability analysis of jet breakup
(4.2.2) In the limit in which the effects of air are small (Weber's limit), we obtain, as shown in Appendix C.1, a dispersion relation (62).We set where s is the growth rate, k is wavenumber, and an overline indicates scaled variables.We find, after some further manipulation of the dispersion relation, that where the Oh is the Ohnesorge number with characteristic lengthscale a, the radius of the cylinder, and I 0 , I 1 are modified Bessel functions of the first kind.
(4.2.3) Varying the Ohnesorge number allows us to examine the impact of viscosity on the maximum wavenumber, giving a dominant wavelength for the breakup of the cylinders into droplets.Firstly, Oh = 0 gives Rayleigh's solution [24]: We compare this inviscid result with various small values of Oh -note, for water and a cylinder of radius of 50 micrometers, Oh ≈ 5.3 × 10 −3 -in Fig. 12.By increasing the viscosity, we decrease the wavenumber for which s is largest: hence, the wavelength of the most unstable mode increases.However, it is also apparent that increasing Oh also decreases the growth rate, so viscosity appears to increase the wavelength of the dominant mode, but to reduce the rate of breakup.We plot the maximum wavenumber and growth rate as a function of the Ohnesorge number in Fig. 13, which confirms our findings.Denoting the void radius as r v , the rate at which the void expands radially is given by

Summary
where γ is the surface energy per unit area, h is the sheet thickness (varying with distance z from the aperture, assumed constant across the sheet width) and ρ is the liquid density (assumed constant).
(5.1.2)If the initial sheet width is w 0 and the sheet fans at angle θ f from either side of the aperture, then the sheet width w at a distance z downstream from the aperture is, By conservation of mass, the cross-sectional area of the sheet must remain constant, so, denoting the initial thickness as h 0 , Far from the aperture z >> w 0 , we approximate the radial velocity of void expansion by where t is time.With a constant liquid sheet velocity U , if the perforation starts at z0 , at time t its centre is at z(t) = z0 + U t.We then integrate (9) to get From this expression we obtain the length fraction of the sheet that is void at time t, where again we have used z >> w 0 .

Webbing Thickness and Volume
(5.2.1) Next, we consider what happens to the liquid that leaves the voids and moves to the sheet.Here, we suppose that this displaced liquid forms a torus encircling the void that it left.Denoting the volume of liquid that would fill the perforation as V p , the volume of an annular ring around the void existing prior to the torus as V r and the torus volume as V t , conservation of volume gives Here R o is the outer radius of the annulus, and is the inner toroidal radius.Thus we obtain the outer radius by solving the cubic equation 6 Late stage droplet size distribution (6.0.1)The droplet breakup has three stages.The early stage breakup -the primary breakup into cylinders, and the secondary breakup into droplets -occurs near the nozzle, and the late stage breakup occurs further downstream, where a droplet distribution already exists.We model this late stage breakup as the evolution of the droplet distribution, coupling it with the early stage breakup through the initial distribution.We aim to predict the size distribution of the droplets that actually reach the ground or crops.
6.1 Model formulation (6.1.1)We model the late stage droplets with a transport equation: where V is droplet volume, x is distance across the flow, t is time, f (V, x, t) is the droplet size distribution function, u D is the droplet velocity, B(V, x, t) is the breakup rate, C(V, x, t) is the coalescence rate and S(V, x, t) is a source term.
(6.1.2)We take slices across the flow to get the domain shown in Fig. 14.We assume a steady velocity field that fans away from the nozzle and is symmetric around the midpoint.The periodicity of the nozzle array suggests u D (0) = u D (1) = 0 and symmetry about x = 1/2 suggests u D (1/2) = 0.The imposed form, as shown in Fig. 15, is therefore 3) The literature surveys many phenomenological models -they have parameters that must be fit by experimental data -describing the dynamics for both the coalesence and breakup of droplets in a mist (see [19,30] and references within).In this work, we neglect coalescence (Appendix E.1 explains why) and, since we have no experimental data, we implement a simple breakup kernel, following [13], to demonstrate evolution of the distribution.The derivation of the breakup term is given in Appendix E.3.  the ratio of the typical time for a droplet to travel from the nozzle to the crops, to the characteristic timescale of the evolution of the droplet distribution, (10).When this ratio is small, the droplet size distribution evolves only a small amount during secondary breakup.When this ratio is large, however, the initial distribution evolves substantially.We find, in Appendix E.2, that a steady state cannot exist for a nontrivial velocity, but we are unconcerned since the evolution occurs only for finite time.

Numerical results
(6.2.1)To simulate the evolution, we implemented a high-resolution, central, conservative finite-volume scheme [15].The domain is discretised to N x × N v resolution, where N x and N v are the numbers of grid points for x and V respectively.The values for the parameters used in the simulations are given in Table 1.
(6.2.2) As discussed above, f (V, x, t = 0) should come from predictions of the early stage breakup.Here, we assume a simplified constant distribution across the domain, i.e. f (V, x, t = 0) = 1 for all (V, x).Results for the evolution of the full predicted droplet size distribution function are given in Fig. 16 for four instances in time. (6.2. 3) The evolution of the size distribution in Fig. 16 acts in the x-dimension to advect droplets from the centre point at x = 1/2, moving them out towards the edges.In the V -dimension we observe that large droplets, that is, large V , break up to produce smaller droplets, so the distribution thic-  kens above smaller values of V .This happens in a nonmonotonic manner: the critical Weber number prevents droplets below some threshold from breaking up, since surface tension forces dominate the dynamics.Very small droplets are, therefore, only produced in very asymmetric (and less likely) breakup events.(6.2.4) Fig. 17  distribution is narrower.The shapes of these distributions qualitatively match the distributions produced by the direct numerical simulations, and are produced in a fraction of the computing time.
(6.2.5)From these plots we see that the two key parameters driving the evolution of the distribution are the breakup rate constant b 0 and the critical Weber number.These two parameters could be tuned, using the droplet material properties, to obtain the most desirable distribution, once the functional dependence of b 0 on these parameters is experimentally determined.

Summary
We have presented a model governing the evolution of the droplet size distribution.The equation is driven by breakup alone, as an assumed large Weber number results in negligible coalescence.We derived an expression for the rate of breakup for the droplets.Note that the validity of the predictions of this model rely heavily on the initial distribution at t = 0 from the early stage breakup, and the value taken for b 0 .We expect b 0 to depend on the physical parameters in the system, and suggest an empirical correlation be used to prescribe a value.
(6.3.2) We imposed an assumed velocity field, but a more realistic approach is coupling the evolution of the droplet distribution with the governing equations for the velocity field.This, however, is computationally expensive.Although simplified, our approach provides fast predictions for the final drop size distribution.
7 Conclusions (7.0.1)We used several mathematical models to improve understanding of how liquid properties affect droplet size.Direct numerical simulations showed that, within the tested parameter range, decreasing velocity gives fewer droplets, but with a bigger size distribution, in the primary breakup.
Increasing viscosity reduces the primary breakup, but has little effect on the secondary breakup, and surface tension has a negligible effect.Linear stability analysis suggests that increasing viscosity leads to larger drops, but this analysis is only valid in the early stages of breakup.Perforations may form in the liquid sheet, and we found an equation for how much of the sheet (the length fraction) is void, due to perforations, at time t.We modelled the late stage drop evolution using a transport equation, and solved this numerically for a uniform initial droplet distribution, finding that the key parameters driving the evolution of the distribution are a breakup rate constant b 0 and the critical Weber number.
(7.0.2) These mathematical models can be used together -for example, the linear stability analysis can be compared with the early stage direct numerical simulations when looking at the secondary breakup, and the intiial droplet distribution for the transport equation comes from the results of modelling that secondary breakup.Using the mathematical models improves our understanding of the mechanisms by which breakup occurs, and can guide experimental work.
7.1 Future work Experiments to: • Isolate the influence of particular dimensionless numbers.
• Characterise the air flow -the air density could be reduced, for example, in experiments in a vacuum chamber.
• Detect the thickness of the sheet -a critical parameter that is missing from our work at present -using interferometry.
• Measure the flow in the sheet, using Particle Image Velocimetry [34].
• Gather more information about inhomogeneities in the flow -knowing where they are is important for understanding their effect.

Theory on:
• A sheet in the generalised lubrication (thin film) theory, looking either at the 1D problem for a sheet or disk, or at the 2D problem to investigate edge effects, with the thickness of the sheet being a function of distance and angle from the source.1 • Each of the three breakup modes -wavy sheet, perforation, and rim -to gain more understanding of relevant parameters.
• Providing an accurate initial condition to the late stage droplet distribution model, examining its breakup rate constant b 0 , and looking further into how to incorporate the velocity field.
Numerics with2 : • An improved numerical methodology to target lower drop sizes (below 100 µm) and a more extended study on the variation in the parameters of interest, obtaining detailed information on the drop statistics.
• More specific (than the two-dimensional jet) configurations in either free or bounded environments.
• A non-zero air velocity in the background, which requires a full three-dimensional investigation.
• Simple non-Newtonian extensions of the developed flow models.
• A full two-stage model in which the droplets emerging from the primary instability are treated as cylinders, which are then subject to breakup under the Rayleigh instability.A A simple instability analysis of fluid sheet breakup (A.0.1)In this section, we find how the stability of a fluid sheet changes with surface tension and velocity.Consider, as shown in Fig. 18, an infinite fluid sheet of thickness h = 2 h moving at speed U in the positive x direction through ambient fluid that is moving at a speed U a .For simplicity, we ignore gravity and atmospheric pressure.We assume there is a symmetric configuration about the centre of the fluid sheet, meaning that we need to consider equations of motion for both the fluid in the sheet and the ambient fluid.We also assume potential flow, so we have the following equations for the velocity potential φ in the fluid sheet: where p is pressure in the fluid and c s is a constant. (A.0. 3) The equations for the velocity potential ψ in the ambient fluid are: lim where p a is the pressure in the ambient, and c a is a constant.
(A.0.4)The interfacial condition for the pressures is: The combined interfacial condition is: , y = h(t, x).
(21) The two constants are related by: The next problem is to understand the role of linear stability.We linearise about the following state: The linear equations of motions then become: The final interfacial condition becomes: To obtain a stability condition, we write: where k is a wavenumber and (the real part of) s is a growth rate.The solutions for φ and ψ are: Inserting these equations into the kinematic boundary conditions yields: On inserting this into the final equation, we obtain a quadratic equation for s, whose real part is the growth rate: (33) To further simplify the problem, we set the ambient fluid velocity U a = 0, and we denote the ratio of ambient and sheet densities by ρ = ρ a /ρ.We then find that The growth rate is then the real part of For large surface tension σ, there is no real part to s, i.e. the growth rate is zero, and so the flow is neutrally stable.For large velocity U , there is a real part to s, i.e. there is sa non-zero growth rate, and the flow becomes unstable.proved to be one of the historical advances in the computational study of interfacial flows.In this section we present relevant concepts of the technique.We also introduce technical aspects behind the implementation of the method in Gerris and highlight prominent numerical and algorithmical features of the package.

B Direct Numerical Simulations: Appendices
(B.1. 2) The Gerris Flow Solver is an open-source stand-alone software package with tremendous resources in terms of both numerical modelling and computational performance, with an architecture oriented towards solving fluid flow problems.The source code is based on the C language and is freely available under the Free Software GPL license.It was created by Popinet [21], with the support of the National Institute of Water and Atmospheric research (NIWA) in New Zealand and Institut Jean le Rond d'Alembert in Paris.
(B.1.3)We now introduce the mathematical formulation as implemented in Gerris.We briefly sketch the approach for our particular problem in the present section and refer the reader to Popinet et al. [22] for details.The equations of motion are where D is the rate of strain tensor D ij = (1/2)(∂ ũi /∂x j + ∂ ũj /∂x i ).All interfacial forces are transferred to the momentum equations in what is commonly called the "one-fluid" formulation -see Tryggvason [28].The physical properties describing each fluid (density, viscosity, permittivity etc.) are included by singular distributions and the same set of equations (36) accounts for the entire domain.The Dirac distribution δ s isolates the surface tension effects to the interface alone, and external forces (such as gravity) are included as appropriate via the F e term.
Under this treatment, a density equation (the others properties are treated in a similar manner) of the general form ρt + ∇ • (ρũ) = 0 (39) which is solved for c and the results substituted into (37)-( 38).The value of c is interpolated across the interface by introducing a small transition layer in its vicinity to smooth the variation of quantities from one region to the other.The sharp transitions in physical properties at the interface are relaxed numerically by applying a filter operator (smoothing) within a single grid cell, thus making grid refinement a crucial factor in obtaining accurate solutions.
(B.1.5)Gerris is one of the few available open-source packages incorporating dynamic mesh adaptivity.Once a base mesh is created, the level of refinement is described using quadtrees in 2D or octrees in 3D.This tree structure provides both a solid tracking possibility of how refinement occurs, and also allows for efficient parallel processing.The user is able to specify a maximum degree of refinement which will represent the depth of the graph in Fig. 19.The number of degrees of freedom is reduced by several orders of magnitude using this very beneficial feature.
(B.1.6)The quadtree (octree in three dimensions) structure was adopted from the onset with the aim to provide a suitable numerical environment for massively parallel computations.The MPI (Message Parsing Interface) library is employed for parallelisation purposes and solid scalability is implemented using a load-balancing algorithm.For further information regarding the technical aspects of the parallelisation in Gerris we refer the reader to Agbaglah et al. [1].
(B.1.7)The Gerris ( [21], [22]) and, more recently, Basilisk Flow Solvers are opensource freeware packages with tremendous resources in terms of both numerical modelling and computational performance.Their qualities, in Syngenta ESGI130 particular in terms of adaptive mesh refinement (AMR), become evident in the study of the problem of drop dynamics in high speed flows.Capturing the evolution of the fluid-fluid interfaces demands an appropriate resolution, accounting for possible topological transitions.Furthermore, vortical structures inside the droplets have been shown to drive the interfacial dynamics [26] and must therefore benefit from a specialised meshing procedure.Perhaps most importantly, we return to the essence of the drop motion and fragmentation, which entail the creation and subsequent tracking of a large number of secondary droplets, which possibly coalesce with other bodies of fluid.A highly refined uniform mesh spanning a large region of computational domain would make the problem intractable.Thus the combination between the volume-of-fluid technique and adaptive mesh refinement is perfectly suited for the problem at hand; in their absence highly accurate calculations at realistic parameter values would become unrealisable with present day computing architectures.Several notable high impact publications have emerged in recent years (see [10,26,29]), combining experimental examinations with insight from numerical solutions using Gerris on drop dynamics problems, thus further reflecting the suitability of the implementation for the proposed lines of investigation.For an extended review of the methodology, as well as applications ranging from electrohydrodynamic control in microfluidic devices to modelling high speed drop impact on aircraft surfaces, relevant introductory sections of the doctoral thesis of the contributer of this section of the report [4] and associated publications therein offer an accessible starting point into the area.

C Linear
where |η| 1, k is the real wavenumber of the disturbance and ω ∈ C is the complex frequency.If Im(ω) > 0, we get unbounded growth.
(C.1.2)We seek axisymmetric solutions for the disturbance velocities in the r− and z-directions, u i (r, z, t) and w i (r, z, t), and pressures, p i (r, z, t), of the form u i (r, z, t) = u i (r)e i(kz−ωt) , w i (r, z, t) = w i (r)e i(kz−ωt) , p i (r, z, t) = p i (r)e i(kz−ωt) .
(42) Upon substituting these into the Navier-Stokes equations and linearising 3) The conservation of mass equation (45) implies the existence of a disturbance streamfunction, ψ i (r), such that Therefore, after eliminating the pressure terms from ( 43)-( 44), we deduce that ψ i (r) satisfies These operators commute and are independent for nonzero ω, so we find solutions of the form which have solutions where I 1 (s), K 1 (s) are first-order modified Bessel functions of the first and second kind, k i = (k 2 − iρ i ω/µ i ) 1/2 and α, β, γ, δ are constants to be Syngenta ESGI130 determined.Since we require the solutions to be bounded at r = 0 in fluid 1 and as r → ∞ in fluid 2, we find that the general solution is (C.1.4)In order to determine the eigenvalue ω, we must apply conditions of continuity of velocity and stress at the interface between the two fluids.Continuity of velocity requires that so that, using (51), we find that After linearising, continuity of tangential stress on r = a gives (C.1.5)Since the perturbed interface is at r = a + ηe i(kz−ωt) , the first two terms of the curvature are The first term gives the curvature of the steady cylindrical profile and is accounted for in the hydrostatic pressure jump.Thus, continuity of normal stress requires where (σ rr ) i = −p i + 2µ i u i (r) is the rr-component of the stress tensor.After some algebraic manipulation and noting that the linearised form of the kinematic boundary condition on the interface gives the relation we derive the final relation between A 1 , B 1 , A 2 , B 2 and ω: To find ω, we therefore most apply the Fredholm alternative to the system (53)-( 54), ( 56), (60), finding that (61) where In general, (61) must be solved numerically for ω(k).
C.1.1 Neglecting the air: Weber's limit (C.1.1)In the limit in which the effects of air are small, we can greatly simplify (61).In this case, the relevant boundary conditions are simply continuity of stress with the air terms neglected.Thus the matrix is simply twodimensional and the resulting dispersion relation for ω(k) is given by I 0 (ka) I 1 (ka) .(62) Note that for k real, it is clear that ω = is, for s real.following equations: This leads to an approximate dispersion relation for the perturbations in the sheet.Instability occurs for sufficiently large fluctuations in surfactant concentration on the surface.It is unclear, however, whether this instability would be a valid mechanism of breakup for fan nozzle jets, since the model assumes that all surfactant is at the surface of the sheet -this may be unrealistic.Even if this assumption is reasonable there is some complex phase behaviour between the surface and surfactant concentration that make the mechanics of any such instability unclear.
E Late Stage Droplet Size Distribution: Appendices E.1 Neglecting Coalescence (E.1.1)Various scaling laws characterise when coalescence matters, based on the capillary and Weber numbers.The capillary number is the ratio between viscous forces and surface tension: where µ is the liquid's viscosity, U 0 in the characteristic velocity scaling and σ is surface tension.The Weber number We is the ratio between inertial forces and surface tension, given by where ρ a is the surrounding air density and D is the typical length scaling for the droplet diameter.Frostad et al. [9] consider the characteristic time it takes for the thin film of air separating two colliding droplets to drain, allowing these to touch and coalesce, arriving at the scaling τ ∼ Ca −1/2 , where τ is the typical drainage time.Other studies [3,23,25] consider models in which the angle of incidence between two droplets is taken into account by a so-called impact parameter, and then for any such parameter value, coalescence occurs within some range of Weber numbers (see, for example, Figs. 1 and 2 in [23]).Crucially, these ranges of Weber numbers are bounded by a critical Weber number beyond which no coalescence will occur.
(E.1.2)Direct numerical simulations suggest that the velocity of the mist of droplets is large.For the range of parameter values of interest here, the Capillary number is sufficiently small to result in typical drainage times much larger than the time it takes most droplets to reach the crops.Similarly, the Weber number is sufficiently large so as to suggest that coalescence will be negligible.Droplets that remain airborne in the wake of the vehicle, however, lie within a slower ambient velocity field, and in this region of the mist, which we neglect in this study, coalescence may play a more significant role.

E.2 Steady state
(E.2.3)If a steady state exists, it satisfies The solution of (72) with u D as given in (11) is given by which is undefined at the zeros of u, and thus a steady state cannot exist for a nontrivial velocity, since droplets overfill some regions of the distribution.

E.3 Breakup term
(E.3.4)We now consider the form of the breakup term in (10)

Figure 4 :
Figure 4: Left: Sketch of the two-dimensional computational domain and key parameters.Right: Snapshot of the flow, from left to right: adaptive computational grid, refined with respect to vorticity and interfacial position; fluid phases (liquid and air); vertical velocity; horizontal velocity.

Figure 6 :
Figure 6: Varying velocity U from 4 m/s -16 m/s: its effect on the number of droplets (left) and their distributions at a selected time (right).

Figure 7 :
Figure 7: Varying surface tension σ from 0.025N/m -0.072N/m: its effect on the number of droplets (left) and their distributions at a selected time (right).

( 3 . 3 . 1 )
Once liquid columns are formed as a result of the primary/wavy sheet instability, they undergo a secondary breakup in the form of a Rayleigh instability, resulting in interfacial rupture and drop formation with ti-Syngenta ESGI130 (a) t = 0.1.(b) t = 0.8.

Figure 9 :
Figure 9: Snapshots of the Rayleigh instability developing on a liquid column as part of a full three-dimensional simulation at two selected timesteps.

Figure 10 :
Figure 10: Variation of liquid viscosity and its effect on the Rayleigh instability.

( 4 . 2 . 1 )Figure 11 :
Figure11: A cylinder of radius of a of fluid 1 lies within an infinite bath of fluid 2 (dashed line).This steady state undergoes sinusoidal perturbations in the zdirection and we look for axisymmetric solutions of the resulting motion subject to suitable boundary conditions.

= 1 Figure 12 :
Figure 12: The scaled growth rate of perturbations s as a function of the scaled wavenumber k for various values of the Ohnesorge number.The solid black line indicates the inviscid Rayleigh solution.

Table 1 :
Number of grid points for x 100 -N v Number of grid points for V 70 -Parameter values used in numerical simulations.

Figure 16 :
Figure 16: Droplet size distribution function at four times using b 0 = 1.

Figure 17 :
Figure 17: Droplet size distribution function as a function of V with varying values for b 0 .

( 7 . 1 . 1 )
The following future work would further understanding of the effects of liquid properties on droplet sizes.

(B. 1 . 4 )
In volume-of-fluid methods, relevant properties such as density, viscosity or permittivity are represented in terms of a volume fraction c(x, t), where c is a generic colour function that takes the value 0 in one fluid and 1 in the other.More specifically we write ρ(c) Stability of an Axisymmetric Column of Liquid: Appendices C.1 A fully two-fluid stability analysis of jet breakup (C.1.1)Taking the setup described in section 4.2, at time t = 0, small sinusoidal perturbations displace the interface between the fluids to r = ηe i(kz−ωt) , More physics could be included in the model in (61), although investigations of that determinant are likely to require numerical treatment.Figure13:The scaled maximum wavenumber kmax and the corresponding growth rate for a range of Ohnesorge numbers.Increasing viscosity decreases both kmax and smax .
4.3.1)Inthissection, we have shown that viscosity increases the wavelength of the most unstable mode in Rayleigh jet breakup, suggesting that the resulting droplets are larger.This linear analysis, however, is only valid in the early stages of breakup, and nonlinear effects will become important as the perturbations grow, but it does provide a useful benchmark for simulations.(5.1.1)Consider a liquid sheet from a nozzle's aperture.At some distance from the aperture, the sheet perforates, causing a void, which grows in size.

Table 4 :
. The number of droplets of volume V is decreased by breakup, and also increased by Notation used in this report: Greek letters

Table 5 :
Notation used in this report: dimensionless numbers