Close Nozzle Spray Characteristics of a Preﬁlming Airblast Atomizer

: The formation of pollutant emissions in jet engines is closely related to the fuel distribution inside the combustor. Hence, the characteristics of the spray formed during primary breakup are of major importance for an accurate prediction of the pollutant emissions. Currently, an Euler–Lagrangian approach for droplet transport in combination with combustion and pollutant formation models is used to predict the pollutant emissions. The missing element for predicting these emissions more accurately is well deﬁned starting conditions for the liquid fuel droplets as they emerge from the fuel nozzle. Recently, it was demonstrated that the primary breakup can be predicted from ﬁrst principles by the Lagrangian, mesh-free, Smoothed Particle Hydrodynamics (SPH) method. In the present work, 2D Direct Numerical Simulations (DNS) of a planar preﬁlming airblast atomizer using the SPH method are presented, which capture most of the breakup phenomena known from experiments. Strong links between the ligament breakup and the resulting spray in terms of droplet size, trajectory and velocity are demonstrated. The SPH predictions at elevated pressure conditions resemble quite well the effects observed in experiments. Signiﬁcant interdependencies between droplet diameter, position and velocity are observed. This encourages to employ such multidimensional interdependence relations as a base for the development of primary atomization models.


Introduction
The design of aircraft propulsion systems, mainly gas turbines, is facing major challenges regarding the reduction of noise, pollutants and greenhouse gas emissions.The combustion chamber of a jet engine is the key component to meet these requirements.The design process of combustion chambers strongly relies on numerical methods, mostly Euler-Lagrangian simulations of the reacting spray.Combined with models for fuel injection, turbulent dispersion, secondary breakup and evaporation of droplets as well as models for combustion and pollutant formation, this approach enables the numerical prediction of pollutant emissions.Besides the air distribution, pressure, temperature and other parameters, the fuel atomization and the resulting spray have a significant influence on the pollutant formation processes [1].Hence, the reliable prediction of pollutant emissions strongly depends on the precise prediction of the droplet trajectory through the combustion chamber, which is influenced by its starting conditions and the airflow through the atomizer [1].The interaction of air and fuel droplets can be captured by the momentum equation [2].In recent Euler-Lagrangian simulations of aero-engine combustion chambers, different methodologies have been used to set the droplet starting conditions.The initial droplet diameters are usually computed using drop size distributions [3][4][5][6][7][8][9], which are derived by primary atomization models.Those correlations are mostly derived from measurements downstream from the injection location.Due to the annular shape of the injector nozzle, the droplets are assumed to be injected at multiple discrete injection points equally distributed on a circle [6,8].The position of the circle downstream from the injector nozzle is selected to fit experimental data.Usually, the initial droplet velocity is purely axial, neglecting the radial and tangential components.Some authors set a fixed initial axial velocity [4,5,7,8], while others use models to predict the initial axial velocity depending on the droplet diameter [3,6].
These common definitions of the droplet starting position and the initial velocity vector are illustrated in Figure 1.In most publications the initial diameter, position and velocity are treated independently, and sometimes a relation between droplet size and axial velocity is considered.Consequently, there is a lack of knowledge if and how droplet diameter, position and velocity are correlated.
Collecting such data for elevated pressure conditions and for realistic injector geometries may be quite challenging, due to complex experimental facilities and restricted optical access.Hence, some attempts have been made to predict the atomization process numerically from first principles.In the context of prefilming airblast atomization, Euler-Euler [11][12][13] and Lagrange-Lagrange [14] approaches have been developed to resolve the gas-liquid interface.Two-phase simulations at elevated pressure have been conducted using both fully Eulerian methods [12] as well as fully Lagrangian methods [15].
In the present work, the fully Lagrangian Smoothed Particle Hydrodynamics (SPH) method is applied for simulating the breakup process of a prefilming airblast atomizer with different trailing edge thicknesses and gas velocities.The main objectives of this study are (i) to provide spray characteristics close to the fuel nozzle and at elevated pressure, (ii) to analyse breakup and spray quantities using descriptive statistics, (iii) to explain relations between ligament breakup and spray formation, (iv) to highlight the multivariate interdependencies between the spray quantities, (v) to show the effects of trailing edge thickness and gas velocity on these multivariate interdependencies within the sprays, (vi) to give detailed statistical insights into the ligament and spray formation, which support the selection of appropriate modelling methods.The paper is organized as follows.First, the numerical method is introduced.Second, the investigated atomizer geometry and operating conditions as well as the details of the numerical setup are presented.Third, the characteristic quantities and the sampling procedure are introduced.Fourth, the airflow field around the prefilmer is described.Then the breakup phenomena as well as the resulting ligament and spray quantities including the multivariate interdependencies between each other are discussed for one baseline case.Subsequently, the effects of trailing edge thickness and gas velocity on these multivariate interdependencies are outlined.Finally, the impact of the present work on modelling primary atomization is discussed.

Numerical Method
The Smoothed Particle Hydrodynamics (SPH) method is a mesh-free method that relies on a Lagrangian representation of the fluid by particles moving with the velocity of the fluid.The particles carry the physical properties of the fluid such as mass, momentum and energy.This method was originally developed for astrophysics [16,17] and later adapted to free surface flow applications [18].In recent years, the SPH method has been applied to air-assisted liquid atomization [14,19,20].
The main advantage of the SPH method in the context of non-miscible multiphase flows is the natural handling of the phase interface.In contrast to Eulerian methods, neither an interface reconstruction nor a threshold for the volume fraction [11] are required.Furthermore, the usage of discrete particles avoids interface diffusion and enables inherent mass conservation.
Since the state of a certain particle depends on its adjacent particles, a weighting function W [18] as depicted in Figure 2 (top) is used to evaluate a function f at the particle location r a by: where V b is the volume of the neighbour particles b, located inside the sphere of influence Ω of particle a, as illustrated in Figure 2 (bottom).The weighting function is defined on a compact support, the so-called 'sphere of influence', that depends on the smoothing length h, and must fulfil mathematical properties such as the unity integral ( W(r − r , h)dr = 1) and the convergence to the Dirac δ-function for h → 0. The Wendland C 2 kernel function [21] is used with h = 1.5∆xwhere ∆x is the mean particle spacing.
With this formalism the Navier-Stokes equations can be cast into an SPH formulation, as explained in [14].To close the system, the Tait equation of state is used.This approach is referred to as Weakly Compressible SPH [22].Surface tension is applied using the Continuum Surface Force (CSF) method [23,24].The liquid-wall interaction (contact angle) is modelled by a correction of the normal vector in the framework of the CSF model [25].
The numerical performance of the SPH code developed at the Institut für Thermische Strömungsmaschinen (ITS) has been demonstrated recently [26].A 3D SPH simulation with this code has proved the capability of the SPH method to predict the primary breakup process accurately compared to experimental data [14,26], in particular the drop size.Using 2D SPH simulations of the primary breakup at ambient pressure, the effect of the trailing edge thickness on the atomization process has been clarified [27].2D SPH simulations of a realistic fuel injector geometry at high pressure conditions [28] gave detailed insights into the film flow development, mixing and spray characteristics.Consequently, the SPH method seems to be an appropriate tool for capturing the atomization process at elevated pressure.

Numerical Setup
In this chapter, a brief introduction to the investigated geometry and operating conditions is given.Furthermore, the details of the numerical setup with respect to the SPH method are presented.

Investigated Fuel Injector
A planar prefilming airblast atomizer, Figure 3 (left), which has been extensively studied experimentally [11,29,30], is investigated numerically using a 2D segment at the trailing edge, as depicted in Figure 3 (right).The trailing edge thickness h a , which has a significant influence on the generated spray [31], is set to 230 and 640 µm.Any other dimension is kept constant in this study.

Investigated Operating Conditions
In order to consider real engine conditions, the air is assumed to be compressed isentropically from ambient conditions to a pressure p g of 5 bar before entering the atomizer.The air is at 464 K with a density ρ g of 8.22 kg/m 3 and a viscosity µ g of 2.575 × 10 −5 kg/m/s.The air velocity u g is chosen to be 70 and 90 m/s.As liquid, the fuel substitute Shellsol D70 is used [30].The liquid properties of D70 are almost identical to those of kerosene, namely density ρ l = 770 kg/m 3 , surface tension σ = 0.0275 kg/s 2 and dynamic viscosity µ l = 0.00156 kg/m/s at ambient conditions.The same liquid properties are used for all three simulations.Due to the short residence time of the liquid on the prefilmer, the liquid is assumed not to be heated up significantly.The liquid film loading is kept constant at V/b = 75 mm 2 /s.This is equivalent to an inlet velocity of u f ilm = 0.9375 m/s at a film height of h f ilm = 80 µm.
All three cases presented in this paper and the corresponding non-dimensional quantities are given in Table 1.The definitions can be found in Appendix A.

Details of the SPH Setup
The 2D domain is discretized by three different types of particles representing the air, the fuel and the walls.In previous simulations (page 4 in [27]), the influence of the mean particle distance, i.e., the spatial resolution on the airflow field and the spray, was studied at ambient conditions and at a medium air velocity (u g = 50 m/s).From these simulations it was concluded that a mean particle distance of 5 µm is a good compromise between accuracy and computational effort.For the simulations presented in this paper the mean particle distance has been reduced to 2.5 µm, since elevated air pressure (p g = 5 bar) and high air velocity (u g = 90 m/s) are covered.This fine spatial resolution is required for capturing the atomization process properly.It results into 11.3Mio particles within the domain.Due to an adaptive time stepping which takes into account the Courant-Friedrichs-Lewy (CFL) condition and the maximum acceleration as well as viscous effects and surface tension, the mean step time varies between 64 and 68 ns.
No turbulence model has been applied.The smallest turbulence length and time scales can be estimated using the channel height of 8.11 mm [30] as large scale, a turbulence intensity of 10 % [30] and taking into account the variation of air velocity.This results into Kolmogorov length scales from 6.7 to 8 µm and Kolmogorov time scales of 6.4 to 9.4 µs (Table 1 in [32]).Hence, the smallest turbulence scales in time and space are adequately resolved by the present simulations, meaning Direct Numerical Simulations (DNS) are conducted.
In the Tait equation of state, the polytropic exponent γ is taken as 7 and the fictive speed of sound are c g = 250 m/s for the air and c l = 25 m/s for the liquid fuel.The background pressure p back is set to 0.42 bar.The contact angle of the wall-liquid-gas interface is kept constant at an angle of 60 • .
Using 1000 CPUs for 72 h for each simulation, 511 to 633 through flow times were simulated.

Characteristic Quantities and Sampling Procedure
The ligament breakup can be characterized using five quantities as illustrated in Figure 4.The SPH particles which form the liquid structure attached to the trailing edge (red and green) are identified using a Connected Component Labeling technique [33].The liquid on the prefilmer (red) is neglected in the following analysis.A length scale d Lig characterizing the size of the ligament (green) is derived from the equivalent area of a circle.The size of the ligament in stream-wise direction is equal to its projected length l Lig from the trailing edge to its maximum x-coordinate.The height of the ligament h Lig is defined by the difference between the minimal and maximal y-coordinates of the ligament.The axial and vertical velocity of the ligament u Lig and v Lig are obtained at the crest of the ligament.The crest is defined by the particle of the ligament with the maximum x-coordinate.To investigate the spray quantitatively, clusters of liquid particles are identified using a Connected Component Labeling technique [33].Each particle cluster with 5 or more particles is considered to represent a droplet.The droplet diameter d Drop is determined by the equivalent area of a circle.These two assumptions lead to a minimum detectable diameter of 6.3 µm.The centroid of the clustered particles is assumed to represents the droplet positions x Drop and y Drop .The droplet velocities u Drop and v Drop are equivalent to the mean velocities of all particles representing the cluster.Furthermore, only droplets located downstream from the trailing edge and which do not touch the top and bottom walls are taken into account.In the present work, the droplet set of each operating point comprises the droplets of all time steps and consists of 137,000 to 243,000 droplets.Due to the high sampling frequency, most droplets are taken into account multiple times during their flight through the domain.Since no recirculation zone is established in the airflow field (Figure 5), the droplets will move towards the outlet to the domain.Droplets hitting the top and bottom wall boundaries of the computational domain are removed from the further analysis by filtering (|y Drop | > 4 mm).

Airflow Field
The prefilming airblast atomization is driven by the momentum transfer from the airflow to the liquid fuel.The airflow around and downstream from the trailing edge is presented in Figure 5 using both the instantaneous velocity magnitude as well as the time averaged axial velocity.
Two instantaneous flow fields are recorded 2.4 ms after the start of the simulation, which is 0.2 ms after the first droplets are formed downstream from the prefilmer (Figure 5a,d).Two more instantaneous flow fields are recorded at the end of the simulations (Figure 5b,e).The time averaged flow fields are obtained from 7.8 ms to 60.6 ms and 58.5 ms, respectively (Figure 5c,f).Hence, for all plots, the prefilmer is wetted and liquid undergoes breakup.It is clearly visible that vortex shedding occurs downstream from the trailing edge (Figure 5a,b,d,e).As expected, the size of the wake downstream from the prefilmer is proportional to the height of the trailing edge (Figure 5c,f).Due to the liquid breakup, no recirculation zone is established.

Phenomena Observed During Breakup
Five snapshots of a breakup sequence are depicted in Figure 6 (bottom), including the recording time of the snapshots (top).Please note that a video of this breakup sequence is available online in the Supplementary Materials of this paper.In this sequence the film waves (i), the temporal evolution of one ligament (iii) as well as its breakup (iv) and finally the dispersion of the resulting droplets (v) can be observed.Due to the wavy film, most liquid is advected within the film waves (i) to the trailing edge (ii).The volume of the wave is either integrated into the liquid accumulated at the trailing edge or the wave has enough momentum to stretch out the accumulated liquid resulting in the formation of a ligament.In most cases, the tip of the ligament (Snapshot 1) reaches into the co-flowing high speed air.This causes a breaking up of the tip and a forcing of the middle part of the ligament to evade into the wake of the trailing edge.Consequently, the ligament is stretched in both the axial and vertical direction, whereas the tip mostly moves only in axial direction (Snapshots 2 and 3).The middle part dips into the lower high speed airflow (Snapshots 2 and 3).This results in the breakup of the tip and the middle part of the ligament (Snapshot 4).The resulting larger liquid structures will disintegrate further into droplets that are widely spread in the vertical direction (Snapshot 5).Meanwhile, the shortened ligament moves upwards, where it is further disintegrated by the upper airflow.This sequence will be repeated with the following film wave reaching the trailing edge (Snapshot 5).Shadowgraphy images of the liquid breakup at the tip of a prefilming airblast atomizer, as studied by Gepperth et al. [11,29,30], are depicted in Figure 7.The geometry is identical to that depicted in Figure 3 (left).The side view perspective on the tip of the prefilmer is presented in the first and third column and corresponds to the sequence shown in Figure 6.The green lines represent the tip of the prefilmer.The images showing the side view are blurry since the depth of field of the camera is very small compared to the width of the wetted prefilmer (|z| ≤ 25 mm).Consequently, the liquid on the prefilmer and the liquid structures attached to the trailing edge are superimposed by each other.
The corresponding top view perspective of the prefilmer is presented in the second and fourth column.From the top view perspective, a wavy film flow, the accumulation of liquid at the trailing edge, the elongation of ligaments and formation of liquid bags with a thick rim and a thin interior lamella as well as their breakup can be observed.These processes are classified as primary atomization.Furthermore, the subsequent breakup of non-spherical liquid blobs into spherical droplets is depicted, which corresponds to secondary atomization.In the side view perspective the stretching, flapping and breakup of both the ligaments as well as the bags can be observed.
Please note that the shadowgraphy images correspond to operating conditions of p g = 1 bar, u g = 70 m/s, h a = 230 µm and V/b = 25 mm 2 /s.Compared to the results of the simulation as presented in Figure 6 only two out of four boundary conditions are met by the shadowgraphy images.This can be explained as follows: First, at higher pressure (p g = 5 bar) and higher film loading ( V/b = 75 mm 2 /s) the number of concurrent breakup events and the amount of liquid present at the trailing edge increases rapidly.Hence, the side view perspective would not provide sufficient insights.Second, images at these conditions are not available.Since the underlying atomization mechanisms are the same for both ambient and elevated pressure, the same phenomena can be observed during breakup.Please note that experiments at low air velocity (u g = 20 m/s) provide many insights in the side view perspective (Figure 4 in [34]).The present numerical observations are in good agreement with the experimental observations.Despite the fact that the present SPH simulations are 2D, most of these effects can be reproduced.In the following we focus on the statistical analysis of the resulting ligaments and droplets and not on the formation process itself.The dynamics of the primary breakup process have been investigated in a previous study [35].

Analysis of the Ligament Evolution
The evolution of multiple ligaments can be analysed using bivariate histograms.In Figure 8a the interdependencies between the ligament length l Lig and height h Lig are depicted.
As the full temporal evolution from minimum to maximum elongation is covered by the data, most of the detected ligaments are smaller than 0.3 mm in height and 1 mm in length.In the following discussion, ligaments smaller than 0.2 mm in length are ignored, as they represent the liquid accumulated at the trailing edge instead of a ligament.Generally, the length to height ratio ranges from 1 to 10, which is reasonable due to the drag force imposed by the co-flowing air.Ligament heights h Lig ≥ 0.4 mm are very unlikely to occur.This value is larger than h a due to the boundary layer which is established on the prefilmer (Figure 5a).The interdependency of ligament length l Lig and ligament size d Lig is presented in Figure 8b.It can be divided into two sections using two lines.One with a steep slope (dash-dot line) and one with a flat slope (dashed line).In the first section the evolution of the ligament is dominated by the film waves.Consequently, d Lig , which is based on the area or "volume" represented by the liquid, strongly increases.In the latter section the stretching of the ligament caused by the drag force of the co-flowing air dominates.Hence, l Lig strongly increases, whereas d Lig changes only moderately.
The relation between the axial velocity of the crest of the ligament u Lig and its length l Lig is depicted in Figure 9a.Obviously, the crest velocity u Lig increases with increasing l Lig .This is reasonable since the probability of the ligament to flap and dive into the high speed airflow increases with increasing l Lig .Due to the drag forces the crest is accelerated.These findings are supported by the positive relation between u Lig and h Lig (Figure 9b).Furthermore, it can be observed that only a few events with u Lig ≥ 7 m/s occur.These result from events where the tip of the ligament is accelerated strongly by the drag forces and evades into the wake of the main part of the ligament.The negative axial crest velocities result from the contraction of the remaining ligament back to the trailing edge.The vertical crest velocity v Lig is symmetrically distributed around zero with a slight trend to negative velocities as depicted in Figure 10a.With increasing l Lig no significant change can be observed, indicating that the flapping mechanism of the ligaments is independent of their length.In a previous numerical study [36] it has been observed that the vortex shedding frequencies of the airflow are at least one order of magnitude higher than the flapping frequencies of the ligaments attached to the trailing edge.However, the dynamics of gas and liquid are not part of the present work.The relation between axial and vertical crest velocity is depicted in Figure 10b.While the remaining ligament contracts back to the trailing edge (u Lig ≤ 0 m/s), the ligaments flap vertically, resulting in both positive and negative v Lig .The ratio of axial to vertical crest velocity ranges from 2 to 20, which is reasonable due to the high speed airflow.In summary, the evolution of the ligaments is based on two mechanisms: The ligament growth driven by the film waves and the stretching of the ligaments by the drag imposed.For the investigated operating point, the acceleration of the crest of the ligament is mainly induced by drag forces.Furthermore, flapping of ligaments is observed for all ligament lengths.

Relations between Ligaments and Droplets
In the following, the interaction of ligaments and droplets is discussed using histograms of both ligament and droplet quantities (blue and orange bins).The number distributions of the ligament length l Lig and the axial position of the droplets x Drop are depicted in Figure 11a.The probability for l Lig ∈ [0 0.5] mm is quasi constant (solid line) and then decreases linearly with a steep and flat slope (dashed line).The first region (solid line) is related to the accumulation of liquid at the trailing edge and the ligament growth.The latter region (dashed line) is a consequence of the ligament stretching and breakup (Figure 6).A characteristic length of the ligaments at breakup cannot be provided since for the simulated time too few breakup events occur.However, the probability of the presence of a drop increases asymptotically from x Drop = 0 to 3.5 mm and then decreases slightly.The presence of droplets in the first region (solid line) can be related to stripping phenomena where some droplets are sheered off from smaller ligaments or film waves.Longer ligaments break up fully or partially, releasing a significant number of droplets.As most of these droplets are too large to withstand the harsh airflow, subsequent breakup into smaller droplets occurs, which results in the increase of the probability from x Drop = 2.0 to 3.5 mm.The slight decrease of the probability more downstream is related to the different axial velocities of the droplets.These observations regarding the distributions of l Lig and x Drop are in good agreement with experimental results (Figure 29 in [29]).The probability of the projected heights of the ligaments h Lig as well as the vertical droplet positions y Drop are depicted in Figure 11b.The ligament height follows a Log-logistic distribution within h Lig = 0 to 0.6 mm.This reflects the flapping mechanism of the prefilmer.Large ligaments reach into the co-flowing air and are disrupted (Figure 6).The resulting droplets are normally distributed in the vertical direction.They are more spread than the ligament crests due to the interaction with the airflow, as to be discussed later.
The amount of liquid contained in the ligaments and droplets can be analysed by using the number distribution of the ligament size d Lig and the droplet diameter d Drop , as depicted in Figure 12, since both quantities are derived from the area of the ligaments and droplets.The dimension of the ligaments is smaller than 400 µm with a peak at 235 µm following approximately a Gaussian distribution with a tail at smaller sizes.The tail results from the fact that the crest of the smaller ligaments is moving slower than the crest of large ligaments (Figures 8b and 9a).Consequently, small ligaments are recorded more often than large ligaments.The largest recorded droplet is about 200 µm.The droplet size distribution shows the typical peak at small diameters and an exponential decline to larger diameters as is expected for the atomization processes.In general, both distributions reflect the disintegration of very large ligaments into much smaller droplets.The sprays generated by the simulations are compared to the corresponding experiment using the Sauter Mean Diameter.As the full temporal evolution of the droplets is captured by the simulations, slowly moving droplets are taken into account more often than fast moving droplets.This bias is counterbalanced by a weight based on the droplet velocity (page 77 in [19]).It is defined as w CFL,i = u Drop,i • ∆t/l.The length scale l is equal to the length of the sampling region (l = 6 mm) and the time scale ∆t is the median between all saved snapshots of one simulation.Thus, the velocity weighted SMD can be defined as: ).The experimental D 32,corr are derived from an experimental correlation which is based on non-dimensional quantities (Appendix B).The resulting SMDs are presented in Table 2.For all three cases, a good agreement between the numerical results and the experimental correlation can be observed.Thus, the trends indicated by the correlation can be reproduced very well by the simulations.The velocity distributions of the crest of the ligament and the resulting droplets are presented in Figure 13.The axial crest velocity u Lig (Figure 13a) follows approximately a Gaussian distribution with a tail at high velocities.A peak is observed at u Lig = 1.5 m/s and negative velocities up to −3 m/s.Negative axial crest velocities can be triggered by the flapping of the ligaments around the trailing edge when its size is not changed and by the contraction of the ligaments after breakup.The maximum crest velocity is about 15 m/s and much smaller than the bulk velocity of the gas (u g = 70 m/s).In contrast, the maximum axial velocity of the droplets u Drop is 76 m/s and thus even higher than u g .This is reasonable, since by the flapping ligament the effective channel height is reduced when it reaches into the airflow (Figure 6).Consequently, the airflow is locally accelerated.The distribution of u Drop is bimodal and can be represented by two Gaussian distributions with peaks at u Drop = 5.5 m/s and 30 m/s.The Gaussian distribution at the slower velocities (dashed line) represents the droplets that are generated directly after the disintegration from the ligament crest.Hence, the range of u Drop ∈ [−2, 12] m/s is fully covered by u Lig .The Gaussian distribution at the faster velocities (dash-dot line) corresponds to the droplets that are accelerated by the fast co-flowing air.This bimodal character of the distribution was also observed in experiments (Figure 33 in [29]).
The vertical velocities of ligament crest v Lig and droplets v Drop are depicted in Figure 13b.Both follow a Gaussian distribution.However, the distribution of v Drop is much broader than that of v Lig .When the crests of the ligaments dive into the fast co-flowing air, their movement in vertical direction is slowed down due to the drag force imposed by the airflow on the body of the ligament (Figure 6).Immediately after the detachment from the ligament, the droplet is no more pulled down by the ligament.Hence, they can move faster than the ligament crest.
As to be expected, the ligaments and the droplets immediately after breakup are closely related, for example, in terms of axial velocity.However, the longer a droplet exist, the more its movement is influenced by the drag forces imposed by the co-flowing air.

Analysis of the Spray Evolution
The marginal distributions of the droplet quantities, as discussed in the previous section, revealed several interdependencies between these quantities.These can be further analysed using bivariate histograms which represent the joint probability.In Figure 14a the joint probability between droplet diameter d Drop and its axial position downstream from the trailing edge x Drop is depicted.This plot can be interpreted as the evolution of the droplet size distribution over time.A primary breakup zone (x Drop ∈ [0, 2] mm) and a secondary breakup zone (x Drop ∈ [2, 6] mm) can be identified.The primary breakup zone is characterized by large droplets (d Drop ≥ 100 µm) and is superimposed by the ligament breakup (Figure 11a).These large droplets will break up into smaller ones.Consequently, the number of small droplets will increase in the secondary breakup zone, whereas the number of large droplets will decrease.
The velocity of the droplets in the streamwise direction is illustrated by Figure 14b.Within the primary breakup zone most of the droplets move rather slow, u Drop ≤ 10 m/s, which is in the same range as the ligament crest velocities as presented in Figure 13a.Furthermore, the velocity distribution becomes significantly broader, because existing droplets are accelerated, while actually generated droplets have the same velocity as the crest of the ligament.In the secondary breakup zone most of the droplets are accelerated to up to 50 m/s.The maximum speed of a droplet depends on its size as depicted in Figure 14c, where the relation between the axial droplet velocity u Drop and the droplet diameter d Drop is presented.It can be observed that the maximum speed of a droplet is inverse proportional to its size.Large droplets (d Drop ≥ 100 µm) experience a significant amount of drag force, but as the inertia dominates, the relative velocity remains high, resulting in the breakup of these droplets.Consequently, these droplets cannot reach high velocities (u Drop ≥ 30 m/s).The resulting smaller droplets and the already existing ones can handle the acceleration better.Thus, these reach higher velocities.The maximum axial droplet velocity increases with u Drop ∝ d −0.5 drop , which is in good agreement with experimental and numerical investigations (Figure 26 in [11]).The structure of the airflow field (Figure 5) strongly influences the relation between axial droplet velocity u Drop and the vertical droplet position y Drop as presented in Figure 14d.Most of the droplets remain in the wake of the prefilmer (y ∈ ±200 µm), where they reach velocities of 30 to 60 m/s.Droplets entering the co-flowing air get significantly faster (u Drop ≥ 60 m/s).The minimum velocity of the droplets increases linearly with increasing distance from the prefilmer.As all droplets entering the high speed airflow move from the inner to the outer region of the computational domain (increasing |y Drop |), Figure 14d represents the acceleration of these droplets over time by the airflow.
The relations between the droplet diameter and both the vertical droplet position and velocity are depicted in Figure 15a,b.Both quantities follow a Gaussian distribution around zero and their width is inversely proportional to the droplet size.Large droplets (d Drop ≥ 100 µm) can only exist in and close to the wake downstream from the prefilmer (y Drop ∈ ±1 mm).Thus, their vertical velocity y Drop ∈ [−10, 10] m/s is relatively small.
It is actually the same range as the vertical velocities of the crest of the ligaments v Lig (Figure 13b).The smaller droplets are expelled into the fast co-flowing air (|y Drop | ≥ 1 mm) and reach higher vertical velocities (|v Drop | ≥ 10 m/s).The resulting droplet velocity vectors are depicted in Figure 15c.The ratio u Drop /v Drop ranges from 0.5 to 80.However, for most of the droplets, the axial component dominates.These vectors are related to the droplet trajectories as depicted in Figure 15d.The droplets are widely spread over the domain, forming a typical spray cone.Usually their trajectory follows a parabolic curve (dashed line).Some droplets cross the domain from top to bottom and vice versa, resulting in more chaotic trajectories.In summary, strong relations between all droplet quantities, namely diameter, position and velocity, are observed in the presented data.Hence, the correct representation of the spray characteristics near the injector is required, for example, for simulations focusing on soot production [37].Consequently, the multivariate interdependencies of the spray quantities as outlined previously have to be taken into account.

Effect of Gas Velocity and Prefilmer Geometry on the Spray Evolution
In this section the effect of gas velocity u g and the trailing edge thickness of the prefilmer h a on the interdependencies within the spray are outlined.Therefore, the kernel densities corresponding to Figures 14 and 15 are estimated and analysed using lines of constant probability (isolines).For the baseline case (p g = 5 bar, u g = 70 m/s, h a = 230 µm) up to three different isolines are plotted in red to illustrate the bivariate distribution.The value of the isolines are indicated by corresponding labels.For the cases to be compared (p g = 5 bar, u g = 70 m/s, h a = 640 µm and p g = 5 bar, u g = 90 m/s, h a = 230 µm) only one line corresponding to a low probability is plotted, since these are much more strongly affected by u g and h a than the lines of high probability.
In Figure 16a the effect of u g and h a on the droplet diameter as well as on the primary and secondary breakup zone is presented using isolines of P = 0.01.Obviously, the increase of u g results in smaller droplet sizes d Drop and the increase of h a results in larger d Drop compared to the baseline case.This is in good agreement with experimental investigations [29,30].For the P = 0.01 isolines, peaks of droplet size d Drop can be observed at x Drop = 1.2 and 1.8 mm.These peaks represent the primary breakup zone.The size of the primary breakup zone is enlarged with increasing h a .This is reasonable, since the wake region of the prefilmer is extended further downstream, resulting in longer ligaments.In contrast, u g seems to have no significant effect on the location of the primary breakup zone.Experimental results (Figure 24 (right) in [29]) indicate that for an increase of u g from 70 to 90 m/s the ligament length will not decrease significantly.Figure 16b reflects the acceleration of the droplets due to the imposed drag force.With increasing u g , the faster droplets reach higher axial velocities u Drop compared to the baseline case, whereas the slower droplets are not affected.This is reasonable, since the faster droplets are located in the high speed airflow.The slower droplets remain in the wake of the prefilmer, where they are just slightly accelerated.The increase of h a increases the size of the wake zone.Consequently, the slower droplets are less accelerated compared to the baseline case.Moreover, the faster droplets reach lower velocities compared to the baseline case.This can be explained by the larger size of the droplets and the fact that the maximum droplet velocity is inversely proportional to its size.In Figure 16c the correlation between u Drop and d Drop is illustrated in more detail.For the slow droplets the influence of u g and h a on d Drop is clearly visible.In addition, the maximum velocity of the small droplets is affected by u g and h a .The relation between the vertical droplet position y Drop and its axial velocity u Drop is depicted in Figure 16d.Increasing u g obviously increases the maximum velocity of the droplets, whereas an increase in h a results into smaller axial droplet velocities.Increasing both u g and h a decreases the vertical spread of the droplets (y Drop ).When increasing u g , the spray cone becomes slender, since the axial velocity increases, whereas the vertical velocity remains the same.The increase of h a leads to a larger wake zone resulting in a longer residence time of the droplets in this region, which will cause a narrowing of the spray cone.Hence, different mechanisms will lead to the same effects in terms of the vertical droplet position y Drop .
The change of the velocity vectors due to an increase of u g and h a is depicted in Figure 17a.Obviously, the former leads to higher axial velocities and the latter to lower ones compared to the baseline case.The lines of constant probability (P = 0.1) scale in size and their shape is identical since all three cases are subject to the same breakup mechanism.Consequently, the droplet trajectories and the spray cones as depicted in Figure 17b are only slightly effected by u g and h a .
In summary, the qualitative characteristics, i.e., the shape of the joint probabilities of the spray quantities do not change with u g and h a .This is reasonable, as all simulations belong to the same regime of prefilming airblast atomization.However, u g and h a have a major effect on the shape of the distributions of d Drop and u Drop .A change in h a also effects x Drop .Regarding y Drop and v Drop , only minor effects of u g and h a can be observed.

Conclusions
Detailed simulations of a prefilming airblast atomizer using the Smoothed Particle Hydrodynamics (SPH) method have been presented.While the simulations are 2D, most of the breakup characteristics observed in experiments can be reproduced.Different effects of the film waves and the drag forces on the ligaments are discussed.Strong links between the ligament breakup and the resulting spray in terms of size, position and velocity of the droplets are demonstrated.The SPH simulations resemble quite well the effects observed by experiments.Hence, reliable spray characteristics close to the fuel nozzle and in particular at elevated pressure can be generated by the SPH method.This enables the usa of SPH data for the development of primary atomization models.Regarding the spray, significant interdependencies between droplet diameter, position and velocity are observed.The effects of gas velocity and trailing edge thickness on these interdependencies are outlined.In future work, these interdependencies will be further analysed and modelled using multivariate statistical methods.

Figure 1 .
Figure 1.Sketch of annular prefilming airblast nozzle illustrating common droplet starting positions and initial velocity vector in magenta, adapted from [10].

Figure 2 .
Figure 2. Top: Illustration of a 2D kernel.Bottom: Particle distribution superimposed with the kernel weight and illustration of the sphere of influence.

Figure 5 .
Figure 5. Instantaneous and time averaged airflow field for p g = 5 bar, u g = 70 m/s, h a = 230 and 640 µm.

Figure 6 .
Figure 6.Breakup sequence at the prefilmer (bottom) and recording times of snapshots (top).

Figure 7 .
Figure 7. Shadowgraphy images of the liquid breakup at the prefilmer tip in side view perspective (first and third column) and top view perspective (second and fourth column) for p g = 1 bar, u g = 70 m/s, h a = 230 µm, V/b = 25 mm 2 /s

Figure 8 .
Figure 8. Interdependencies between ligament length scales for p g = 5 bar, u g = 70 m/s, h a = 230 µm (a) Axial crest velocity and length (b) Height and axial crest velocity

Figure 9 .
Figure 9. Interdependencies between axial crest velocity and ligament length scales for p g = 5 bar, u g = 70 m/s, h a = 230 µm.
(a) Vertical crest velocity and length (b) Vertical and axial crest velocity

Figure 10 .
Figure 10.Interdependencies of the vertical crest velocity for p g = 5 bar, u g = 70 m/s, h a = 230 µm.
(a) Axial crest velocity and axial droplet velocity (b) Vertical crest velocity and vertical droplet velocity

Figure 16 .Figure 17 .
Figure 16.Effect of gas velocity and prefilmer geometry on the interdependencies of the spray quantities for p g = 5 bar.

Table 2 .
Sauter Mean Diameters of Smoothed Particle Hydrodynamics (SPH) simulations and experimental correlation.