Euler–Lagrangian simulation of the fuel spray of a planar prefilming airblast atomizer

The pollutant emissions of aircraft engines are strongly affected by the fuel injection into the combustion chamber. Hence, the precise description of the fuel spray is required in order to predict these emissions more reliably. The characteristics of a spray is determined during the atomization process, especially during primary breakup in the vicinity of the atomizer nozzle. Currently, Euler-Lagrangian approaches are used to predict the droplet trajectories in combustor simulations along with reaction and pollutant formation models. To be able to reliably predict pollutant emissions in the future, well-defined starting conditions of the liquid fuel droplets close to the atomizer nozzle are necessary. In the present work, Euler-Lagrangian simulations of a generic airblast atomizer are presented. The starting conditions of the droplets are varied in the simulations by means of a primary breakup model, which takes into account the local gas velocity when predicting the droplet diameter. The objective of this work is to determine the optimal parameters of the probability density functions for the starting position and the starting velocity of the droplets. Spray properties observed in the simulations are used to qualitatively evaluate the major effects of the distribution parameters on the spray and the suitability of the primary breakup model being applied. Hence, the spatial distribution of an experimental spray can be reproduced using a statistical model for the droplet starting conditions.


Introduction
Increasingly stringent emission regulations are the driving forces for the development of civil aircraft jet engines.
To comply with environmental regulations and to prevent anthropogenic climate change, the substantial reduction of noise, pollutant and greenhouse gas emissions is indispensable [1]. In this regard, the combustion chamber is the key component featuring a great potential to satisfy these requirements. Fuel injection is a crucial element of the combustion process, since the formation of soot and NO x emissions is strongly related to the atomization and properties of the liquid fuel spray [2]. In practice, prefilming airblast atomizers are mostly used due to fine atomization and little change in performance over a wide range of operating conditions [3]. A typical swirl atomizer nozzle employed in modern aircraft engines is shown in Fig. 1.
The fuel supplied to the prefilmer surface is transported to the atomizing edge by the high-speed air flow via momentum transfer. At the prefilmer tip the fuel accumulates before bags and ligaments are formed and disintegrate [5,6] (see Fig. 2). For the specific improvement of combustion characteristics, a comprehensive understanding of the fuel atomization process and in particular of primary breakup is required. Understanding the fuel disintegration is one of the keys for reducing pollutant emissions to meet the challenges of future aircraft engine generations.
Experimental investigations aiming at the measurement of pollutant emissions at real engine conditions are complicated by high pressures and temperatures inside the combustor. Frequently, numerical simulations are performed in addition to experimental test series to predict emissions during the conceptual design stage of combustion chambers more reliably. Since the formation of pollutants is significantly influenced by the fuel spray features, accurate numerical predictions of initial droplet properties and trajectories are critical [8].
The mechanisms of primary breakup in prefilming airblast configurations leading to wide droplet size distributions are not understood in detail. In recent years, this lack of knowledge has been addressed by numerical simulations, which highly resolve the liquid disintegration process close to the trailing edge of prefilming airblast atomizers. To this end, different numerical approaches have been employed. Using the Smoothed Particle Hydrodynamics (SPH) method, characteristic breakup phenomena known from experimental findings of Gepperth et al. [6] like liquid accumulation, ligament flapping and bag breakup are captured well numerically [8,9]. In simulations based on the Volume of Fluid method [10,11] and the Level Set method [12] the same good agreement with experimental observations is achieved. Precise numerical investigations enable a comprehensive understanding of primary atomization and serve as a basis for the development of new primary breakup models. However, very fine resolution in time and space is required, since primary atomization is a multi-scale phenomenon. Furthermore, special interface tracking techniques, which are capable of handling high density ratios possibly leading to numerical instabilities, are essential. Consequently, such simulations are related to extremely high numerical costs [8] and, hence, they are not viable for full combustor flow predictions including atomization and spray propagation. To overcome the aforementioned deficiencies, mostly Euler-Lagrangian simulations are performed to study the distribution of the reacting spray by computing the trajectories of the liquid droplets. The trajectories through the combustor strongly depend on the air flow through the atomizer and the starting conditions of the droplets in the simulation [2].    [7] Therefore, well-defined starting conditions is the missing link for predicting spray characteristics properly.
Holz et al. [9] discuss different methodologies for the definition of droplet starting conditions in Euler-Lagrangian simulations of aero-engine combustors. In most cases, drop size distributions are used to determine initial diameters. These probability density functions are computed by means of primary atomization models [13][14][15], which rely on correlations derived from experimental measurements. The droplets are injected into the flow field at several discrete locations homogeneously distributed over a circular ring in order to mimic the annular shape of the injector nozzle. The axial position of this ring is shifted relative to the nozzle exit plane [16][17][18]. The initial axial velocity is either kept constant [17,19] or does correlate with the droplet diameter [18,20]. The imposition of a tangential velocity component is also investigated [16]. In most of the recent studies an assumption is made that initial diameters, velocities and positions of the droplets are independent from each other. Furthermore, the droplet starting conditions in state-of-theart methods are a priori limited to constant values and do not rely on statistical distributions.
In the present work, Euler-Lagrangian simulations of the fuel spray emerging from a planar prefilming airblast atomizer are performed. To consider effects of primary atomization influencing the spray characteristics, the primary breakup model PAMELA [15] is employed in the simulations. The authors propose an extended version of this model featuring parameters for droplet starting position and starting velocity. These initial droplet properties are based on normal distributions each determined by user input. Consequently, the stochastic nature of primary atomization is modeled in the present novel approach. Furthermore, a systematic variation of droplet starting parameters is presented. The main objectives of this study are (1) to assess the influence of the starting parameters on the evolution of the spray, (2) to explain relations between starting parameters and spray angle as well as Sauter mean diameter, (3) to find an optimal setup which is capable of reproducing experimental spray data of a generic prefilming airblast atomizer and (4) to highlight the potential of PAMELA extended to be implemented in CFD solvers, which are used for combustor emissions predictions. The present study clearly focuses on droplet starting conditions in the Euler-Lagrangian framework and on the capability of the introduced set of starting parameters as a calibration tool to accurately predict measured spray details. Highly resolving the aerodynamic features in the wake of the prefilmer trailing edge as well as varying operating conditions or geometric configurations is out of scope of this work. This paper is structured as follows. First, the numerical model used for the simulations is introduced. In that context, the strategy of droplet injection is explained. Second, the results of the gas phase simulation used for validation of the flow field are presented. Third, the variation of the droplet starting parameters is studied and the simulated spray properties are compared to experiments. Finally, the conclusions of the presented work are given.

Numerical model
In this section, the numerical model of the planar prefilmer is presented. The model includes the geometry derived from an experiment and the setup of the simulations. The extended PAMELA model, which is used to set the droplet starting conditions, is introduced.

Experiment and computational domain
At the Institut für Thermische Strömungsmaschinen (ITS) experimental investigations have been performed focused on liquid atomization at the atomizing edge of a prefilming airblast injector [6,21,22]. The experimental setup is shown in Fig. 3. A two-dimensional generic atomizer was used to ensure optical access. This atomizer represents an abstraction of an annular airblast nozzle typically found in aircraft engines. Air passes by on both the top and the bottom of the symmetrical planar prefilmer in longitudinal direction. The liquid is ejected through holes and a thin film is formed on the top surface of the prefilmer. By the momentum transfer from the air flow to the liquid phase, the fuel is driven downstream and atomized at the trailing edge. Optical measurements of the air velocity field (LDA and PIV) as well as of the liquid spray (shadowgraphy and PLTV [23]) are used subsequently to compare the experiment with the numerical simulations. The shadowgraphy images cover the far field of the downstream region of the prefilmer up to 70 mm. From these images, probabilities of residence of the droplets and spray angles are obtained. The PLTV measurements provide Experimental test section of the planar prefilmer, adapted from [6] droplet diameters, axial positions and axial velocities at near field of the prefilmer up to 12 mm downstream the atomizing edge. The geometry of the planar prefilmer employed in the numerical investigations is derived from the experimental setup. A mid-plane cross section of the three-dimensional computational domain is shown in Fig. 4. The prefilmer is attached to an upstream inlet section and a downstream plenum. The outlet is located at the end downstream of the plenum. The dimensions of the plenum are larger than the prefilmer length in order to avoid back flow effects at the outlet. Beside the inlet and outlet patches, the other boundaries of the domain are considered as walls. The geometrical dimensions describing the computational domain are listed in Table 1.

Extended PAMELA model
The primary atomization model PAMELA (Primary Atomization Model for prEfilming airbLAst injectors) was developed for prefilming airblast atomization by Chaussonnet et al. [15] and calibrated with experimental data obtained by Gepperth et al. [6,21]. The model is based on phenomenological observation of the liquid breakup mechanism and provides the complete drop size distribution of the emerging droplets. A Rosin-Rammler distribution function is used for describing the droplet number probability density function: where m is the scale parameter indicating the characteristic diameter of the droplets and q is the shape parameter controlling the width of the density function. The selection of the two parameters is explained in detail in [15]. The diameter of a droplet is designated by D. In Fig. 5 the procedure for defining droplet starting parameters in the extended PAMELA framework is shown schematically. The part of the model as proposed by Chaussonnet et al. [15] is highlighted in blue. The input parameters of the model are the fluid properties of the gaseous ( g , g ) and liquid ( l , ) phase, the gas velocity u g and the geometry of the atomizer ( h a , L film ). The two Rosin-Rammler parameters m and q are determined by means of calibrated model constants ( C 1 to C 5 ). The probability density function (PDF) describing the number of droplets is then obtained from Eq. (1). Integration yields the corresponding cumulative distribution function (CDF): Using the inverse transform sampling method, actual droplet diameters can be determined. To this end, a random variable X is used, which is equally distributed over the interval of [0, 1]. The inverse function of the cumulative distribution F −1 0 ∶ X ↦ D then provides a diameter D that is used for droplet injection in Euler-Lagrangian simulations.
In the extended PAMELA model, starting conditions for the initial position and the initial velocity of a droplet are additionally specified (indicated in yellow in Fig. 5). A mean value and a standard deviation are defined for the vertical and axial components of position and velocity via user input. The corresponding starting parameters y, z, v and w are determined for each droplet to be injected based on the respective Gaussian distribution by inverse transform sampling. Hence, together with the diameter D, a total of five independent starting parameters of each droplet is employed in the sense of univariate statistical modeling. The position of a droplet in lateral direction is identical with the x-coordinate of the center of the respective injection face at the atomizing edge (see Fig. 6). The lateral starting velocity is zero. Consequently, for the initial position with i and i , i ∈ {y, z}, and for the initial velocity with j and j , j ∈ {v, Fig. 4 Cross section of the three-dimensional computational domain of the planar prefilmer (not true to scale)

Setup of the simulations
The CFD code used to perform the Euler-Lagrangian simulations is a Rolls-Royce in-house code [24,25]. For the prediction in this work, a low Mach number approach with SIMPLE pressure correction method is applied [26]. Subgrid scale modeling as part of the large eddy simulations (LES) is done by the dynamic Smagorinsky model proposed by Germano et al. [27]. This model provides a correct asymptotic wall behavior and a good behavior in shear regions, which are present in the wake of the prefilmer. In this study, the low-pass filter length is implicitly related to the mesh cell size h. When performing LES calculations, it is imperative to take account of quality assessment measures ensuring a good resolution. Therefore, an index of quality proposed by Celik et al. [28] is investigated in Appendix A.
The transport equations are discretized and solved on an unstructured computational mesh with 5 million tetrahedral elements. In contrast to the plenum where h = 20 mm , the area of the prefilmer is represented by a finer mesh with a cell size of h = 1 mm . The refinement close to the atomizing edge is important to ensure adequate flow physics in that area (for mesh topology cf. Appendix A). To prove a grid independent solution in terms of an agreement with experimentally measured velocity profiles, further gas simulations are performed on computational meshes with each 10 and 20 million elements. In comparison with the coarse mesh, the latter two meshes are mainly refined in the downstream area of the prefilmer trailing edge. Results obtained on the advanced meshes are not presented in the following of this paper, since only insignificant differences in flow profiles between all of the meshes exist. For the discretization of the temporal as well as the spatial derivatives only 2nd order schemes are employed. The simulated physical time for the gas phase is 30 ms with a time step size of 1 × 10 −6 s . The temporal averaging of the flow quantities starts after a system run-in period of 15 ms.
At the inlet a constant velocity of u in = 19.18 m∕s is imposed. This value results in an air bulk velocity of u b = 60 m∕s averaged over the height of the flow channels above and below the prefilmer [29]. This operating point is selected for the simulations in the present work and corresponding experimental results are used for comparison. Furthermore, a turbulence intensity of 10% is imposed at the inlet. At the outlet the normal velocity gradient is set to zero. All solid walls are represented by no-slip conditions.
In this work, the Euler-Lagrangian approach is used to predict the droplet dynamics. The gas flow as the continuous carrying fluid phase is described by transport equations  within the classical Eulerian framework. In contrast, the droplets represent the disperse phase, for which a Lagrangian system of reference is used. The position of each individual droplet is calculated by means of a particle tracking code based on the underlying flow field. In this regard, the Lagrangian particles are treated as point-sources. The equation of motion of the droplets is used to predict the positions at each time step and, hence, the trajectories. Mainly assuming that the density of the droplets is much higher than that of the carrier fluid and the droplet size is small compared to turbulence integral length scale, drag and gravitational forces are determining for droplet dynamics [30]. Hence, the system of equations of particle motion can be written: where x d , v d and m d are the droplet position, velocity and mass. The term t represents time, F D and F g the drag and gravitational force.
When employing the Euler-Lagrangian approach, special emphasis has to be put on the preparation of the computational mesh. Farzaneh et al. [31] state that on the one hand, the mesh must be fine enough to accurately solve the governing equations and predict distinct features of the flow. On the other hand, the size of the mesh elements used to resolve the continuous phase should be larger than the size of the particles approximated as point-sources. Otherwise the drag force of the particulate phase is not calculated correctly. Askarishahi et al. [32] indicate that when flow details around individual particles are not resolved, the mesh cell size varies between 2 and 15 times the particle diameter.
The interaction between the gas phase and the disperse phase is accounted for by a momentum exchange using a source term in the balance equation. In the present work, a two-way coupling approach is employed. Hence, the mutual influence of the two phases is taken into account. Models accounting for secondary breakup and evaporation of the droplets are not used in the present study.
The spray simulations are started after the run-in period of the pure gas phase predictions, i.e. after 30 ms. The simulation time of the liquid spray covers at least 20 ms for all scenarios of droplet starting conditions in order to obtain a sufficiently large data base for the statistical analysis. The simulations were performed each with 448 CPU cores and for 48 h computational time on a HPC cluster.
The primary atomization model PAMELA is implemented in the CFD code with extensions regarding the initial position and velocity of the droplets. In Fig. 6 the principle of droplet injection, which is carried out at each time step during a spray simulation, is shown schematically. Beside the prefilmer, the patch of the atomizing edge is illustrated. The face center of each surface cell constitutes a virtual injection point for the droplets. From the face center, the droplets are shifted axially and vertically within one time step towards the actual starting position of the trajectories. In order to keep the numerical effort at a minimum, parcels are used which represent a monodisperse accumulation of individual droplets [33]. These parcel droplets do not move independently from each other, but are combined as cluster. Hence, each parcel represents the diameter and number of multiple droplets. The droplet diameter, which is statistically assigned to each parcel, is determined by the PAMELA model. Consequently, the droplet count of each parcel n d ∼ D −3 can be determined using the total liquid volume to be injected and the number of injection faces.
In the present study, the liquid film flow on top of the planar prefilmer (see Fig. 3) is not predicted and the film thickness is not considered in the geometry of the atomizer profile. Jones et al. [17] choose a similar approach in numerical investigations of spray combustion as liquid film breakup and droplet formation are not modeled directly. By tuning spray inlet conditions for the Lagrangian droplets, measured profiles downstream the swirl injector can be reproduced. Sanjosé et al. [16] take account of the liquid film thickness in the initial droplet positions by adjusting the radial extent of an annular injection ring. Furthermore, it is noted that the liquid film thickness is not an input parameter of the PAMELA breakup model (see Fig. 5). Sattelmayer and Wittig [5] concluded from experimental investigations that the film thickness is not to be considered a relevant parameter for the atomization process. In more recent work, experimental evidence was found by Gepperth et al. [6] that the film flow is decoupled from the breakup process due to the liquid accumulation at the prefilmer trailing edge.
The operating point of the experiments considered in the simulations is characterized by a liquid mass flow of ṁ fuel = 2 g∕s . To avoid injecting too many extremely small droplets, a minimum permissible diameter of D min = 5 μm is defined. Furthermore, for reasons of numerical stability, droplets larger than D max = 600 μm are not generated. The selected diameter range represents a typical characteristic drop size distribution for the planar prefilmer at the operating point considered [22,34]. Based on the underlying drop size distribution calculated by the PAMELA model, only few parcels ( 7% ) assigned with a diameter larger than 100 μm are injected. Hence, the Lagrangian point-source approach is valid, since the smallest cell with h = 1 mm is in general one order of magnitude greater than the droplet diameter. The reference plane (cf. Fig. 6) is used to evaluate the air velocity U 0 of the free flow at each time step. It is an arithmetic mean value measured over the reference plane. With the velocity U 0 , the gas velocity u g is calculated according to u g = 0.7 U 0 [15]. This is the input velocity parameter to the PAMELA model. Thereby, the interaction between the film flow and the high-speed air flow is modeled, since the gas velocity seen by the liquid is estimated by a fraction of the bulk velocity [29]. Hence, the influence of a reduced air flow velocity on the drop size distribution is factored in without predicting the film flow itself.
Different scenarios of droplet starting conditions are studied by varying the eight initial droplet parameters. The effects of the different parameters on the spray characteristics are to be elaborated. The scenarios for the mean values and standard deviations of the droplet starting parameters are listed in Table 2. In the following, the meaning of the parametric values is explained. The mean vertical position of y = 0.395 mm is composed of half the thickness of the atomizing edge of 320 μm and the film thickness of 75 μm measured experimentally [35]. The mean axial position of z = 4 mm is based on the mean breakup length of the ligaments as found in the experiment [35]. In recent SPH predictions focusing on primary breakup characteristics, the interdependencies between droplet initial properties are presented. In this context, the vertical flapping mechanism of the ligaments at the atomizing edge of the planar prefilmer due to interaction with the highspeed air flow is captured. Statistical analysis reveal that the vertical ligament velocity is symmetrically distributed around zero resulting in both positive and negative values [9]. Therefore, the standard deviation of the vertical velocity in this study is set to v = 4 m∕s with a mean value of v = 0 to mimic the flapping ligament behavior related to vortices detaching at the trailing edge. The mean axial velocity of w = 11 m∕s is derived from an empirical correlation for the velocity of primary droplets (cf. Appendix B) determined in the experiment [35].
In Table 3 the material properties of the fluids used in the simulations are presented. The properties of the fuel Jet-A are similar to the substitute fuel Shellsol D70 as used in the experiments.

Results of gas simulation
A contour plot of the time averaged velocity magnitude U of the gas flow is shown in Fig. 7. The mid-plane ( x = 0 ) of the prefilmer section is examined. This section contains the relevant downstream region of the prefilmer which is of interest for the subsequent analysis. The flow enters the domain with the velocity imposed at the inlet. After passing the inlet zone, the flow is accelerated through the convergent nozzle. The accelerated and homogenized flow hits the leading edge of the prefilmer and is split up between the upper and lower channel. Due to the wing-shape of the profile, the flows are accelerated again in the front area of the prefilmer. Further downstream, turbulence is developed and the characteristic bulk velocity of the investigated operating point of 60 m/s is reached. At the atomizing edge the flow expands and due to the high momentum, the air jet penetrates deep into the plenum, while it opens up slightly. In the plenum, backflow occurs around the jet in the form of large recirculation zones.
A comparison between LES data and the experiment is carried out to validate the flow field of the gas phase. Again it is pointed out, that it is not the present objective to precisely resolve features like vortex shedding and recirculation in the prefilmer wake considering the computational effort    Fig. 8 and compared to LDA measurements. The axial distance to the atomizing edge is ẑ = 0.3 mm . In the axial velocity profile, the wake of the prefilmer is clearly visible in the range of −1 mm <ŷ < 1 mm . The periodic boundary layer separation at the atomizing edge due to the geometrical discontinuity creates a small recirculation zone. The backflow in this area at ŷ = 0 is reflected in the velocity profile as a negative ditch which is captured quite well by the simulation. With increasing distance from the symmetry line ŷ = 0 , the bulk velocity in the air jet of u b = 60 m∕s is almost reached in the simulation. The reason for the slight underestimation of the velocity in this area is the prefilmer channel height d, which is about 2% higher than in the experimental test rig ( d = 8.11 mm ). The sharp transition from the wake to the jet zone in the LES profiles is probably due to a too coarse mesh resolution close to the prefilmer walls and the resulting inaccurate representation of the boundary layer flow.
Characteristic for the profile of the mean vertical velocity near the atomizing edge is a shape symmetrical to the origin. This can be explained by the non-parallel orientation and the inclination of the flow channels towards the center (see Fig. 4). The two local peaks at ŷ = ± 0.5 mm can be attributed to the increasing curvature and deflection of the streamlines in vertical direction when approaching the atomizing edge. In this region an overestimation of the velocity is observed, which can be explained by an insufficient spatial resolution. The apparent discrepancies between experiment and simulation below the prefilmer result from the fact that the experimental profile is not ideally symmetrical. This might be caused by uncertainties in the laser positioning during LDA measurements.
In Fig. 9 the profiles of the mean axial velocity w are presented together with PIV measurements for the positions ẑ = 5 mm and ẑ = 20 mm . It can be concluded from both profiles that the experimental data is captured well by the simulation. The two local velocity maxima are located within the air jet because of the merging of the two individual channel flows. At the center ( ŷ = 0 ) a minimum can be observed, which is due to the wake downstream the prefilmer. A difference in the geometry is the reason for the deviation from the experiment in the outer jet regions. The cavities above and below the atomizer module in the experiment (see Fig. 3) trigger large recirculation vortices in the plenum. In contrast, the plenum in the computational domain is bounded by vertical walls (see Fig. 4).
The data presented in this chapter demonstrate that there is a good agreement between the velocity profiles of the experiment and the simulation. Hence, the flow field downstream the prefilmer can be considered as validated. The present numerical setup based on large eddy simulation provides a robust and accurate approach for the following Euler-Lagrangian simulation of the liquid spray. A finer mesh with a higher number of cells at the prefilmer wake is obviously not required, since small improvements of the velocity profiles do not justify the much higher computational effort during spray simulations. Furthermore, reducing the smallest cell size below the largest droplet diameter will result into the violation of the Lagrangian point-force model employed in this study. Hence, the numerical mesh with 5 million elements is used for subsequent spray simulations.

Results of spray simulation
In this chapter, the results of the Euler-Lagrangian simulations of the fuel spray are presented and discussed. The main effects of the droplet starting conditions on the evolution of the spray characteristics are qualitatively evaluated by means of a parametric study. A setup is derived which best matches the experiment.
For illustration of the spatial spray topology in the far field of the prefilmer, probabilities of residence of the droplets up to an axial distance of 70 mm are used. These maps are generated from 51 drop samples over a simulation observation period of 5 ms. More than 90,000 parcels representing 1.4 million droplets are taken into account at each time step, which is considered to be a representative sample size. A grid with a resolution of 25 × 25 μm is superimposed over the region of interest. If there is at least one parcel in a grid cell at any time step, this cell is flagged. Consequently, a distinction can be made between liquid and gas. The relative frequency of the flags yields the corresponding probability map.
For determining the spray angle , marker on the outer boundaries on both sides of the spray are determined in vertical direction ( ẑ = const. ) using the 10% and 90% quantile of a statistical regression function. The angle of intersection of two best-fit straight lines of the boundary points is considered as the spray angle.
The calculation of the Sauter mean diameter D 32 of the spray in the prefilmer near field is carried out for discrete droplet sizes in a spray with N droplets as follows: The Sauter mean diameter represents a characteristic diameter. It describes the diameter of a droplet whose ratio of volume to surface area is identical to that of the entire spray.
Spray angle and Sauter mean diameter for the simulations S1 to S5 and for the experiment are compiled in Fig. 10. In the subsequent analysis the effect of the different parameters used to characterize the droplet starting conditions will be discussed.

Effect of vertical starting position
To assess the effect of the vertical starting position y on the spray properties, simulations S1 and S2 are examined in this section (cf. Table 2). In both simulations the droplets are injected with starting velocities of zero behind the atomizing edge. Thus, all droplets are accelerated from the state of rest by the gas flow.
In Fig. 11 the probabilities of residence of the droplets for S1 and S2 are shown. Due to the symmetrical starting parameters, both cases result in almost symmetrical spray distributions. Given the axial distance of the injection from the prefilmer of z = 4 mm , droplets are detected only downstream of this position. The injection on the center line ( ŷ = 0 ) leads to a negligible vertical deflection of the droplets in simulation S1. Hence, the spray cone opens up just slightly and is very narrow. The dispersion of the vertical position of y = 1 mm in simulation S2 results in a much wider spray topology. The width remains approximately constant in axial direction in the far field. With regard to the spray angles, overall a strong underestimation of the experimental spray angle can be observed for both cases S1 ( −88% ) and S2 ( −95% ). The Sauter mean diameters are in good agreement with the experiment showing a relative deviation of less than 4% . The selection of a small starting velocity for both cases leads to a high percentage of large droplets in the near field. Their relaxation time d ∼ D 2 is much greater than that of small droplets. Consequently, they will only weakly be accelerated by the air flow and remain at their current state of rest for a longer time after injection. In contrast, small droplets are accelerated faster due to their lower inertia.
In conclusion, the dispersion of the vertical position y has the effect of increasing the spray width. The higher variation of the vertical position slightly reduces the spray angle. However, this influence can be regarded as small. The comparison of Sauter mean diameters shows that the selection of initial velocities of zero for S1 and S2 leads to a droplet spectrum in the near field of the prefilmer which matches well the experiment. Fig. 10 Spray properties for simulations S1 to S5 and experiment Fig. 11 Probability of residence of droplets for simulations S1 (left) and S2 (right) 1 3

Effect of initial velocities
The simulations S3 and S4 are compared with each other in order to investigate the effect of the initial vertical velocity v on the spray properties (cf. Table 2). In both cases, the droplets are inserted at a fixed position in the wake of the prefilmer with an axial velocity of 11 m/s. The comparison of S3 with S1 also enables conclusions to be drawn on the effect of the axial velocity w.
To illustrate the spray topology, Fig. 12 shows the probabilities of residence of the droplets for the simulations S3 and S4. The dispersion of the vertical velocity v to mimic the flapping characteristic causes a significant increase in the spray width in the case of simulation S4. Hence, droplets are present at a greater distance from the center line than in simulation S3. The vertical component of the velocity vector leads to curved trajectories. The droplets gradually drift away from the axis in the perpendicular to the flow direction. Furthermore, the initial vertical momentum results in a rapid opening of the spray in the near nozzle region. Moreover, it can be noted that simulation S4 shows a much larger spray angle than S3 ( −88% compared to experiment). The spray angle for S4 is 3.5 • and comes close to the experimental value of 4.1 • (underestimation of 15% ). Regarding the Sauter mean diameter, no difference can be observed between simulations S3 and S4. The relative deviation from experimental data is −23% . However, the value of D 32 for both cases is lower than in simulation S1, in which the droplets are injected at the state of rest. As a result, the axial starting velocity causes a shift of the drop size distribution towards smaller diameters in the near field. This effect can be explained by the significant influence of the initial velocity on the large and inert droplets. Due to the initial momentum, the percentage of large droplets in the near field decreases.
In summary, it can be stated that the vertical starting velocity is an effective tool for increasing the spray angle. It also affects the spray width. An initial velocity in axial direction shifts the weight of the drop size distribution towards smaller diameters. The Sauter mean diameter of the droplet spectrum in the near field is consequently reduced.

Optimal setup of initial droplet parameters
In this section, the setup based on the previous parametric study is presented which will best match the experiment. It is intended to tune simulation S5 to most accurately reproduce the experimental results (cf. Table 2).
In order to qualitatively evaluate the spatial composition of the spray of the S5 simulation, Fig. 13 shows the probabilities of residence and snapshots of the droplet spray for S5 and the experiment. The following similarities can be observed. Regarding the vertical velocity v, the same parameter settings are used for S5 as in S4. Hence, the global width of the spray as well as the spray angle ( −12% ) are in good agreement with the experiment. The slight increase of the spray angle in simulation S5 compared to S4 can be explained by the combination of the vertical starting velocity and the injection position of y = 0.395 mm . Additionally, the slight curvature of the probability distribution of the experiment and the opening of the spray close to the prefilmer is captured by the simulation. Discrepancies between the experiment and the simulation S5 are discussed in the following.
The formation of a bulge close to the nozzle, which is observed in the experiment, is to be attributed to the flapping character of the ligament detachment at the atomizing edge in the turbulent co-flowing air. As the ligaments move up and down, a varying vertical momentum is imposed on the detached droplets. This phenomenon can only be reproduced to a limited extent by the simulation via the dispersion of the vertical starting velocity. Hence, the simulation shows a continuously expanding spray close to the prefilmer trailing edge. A further essential difference between simulation and experiment is the decrease of the probability of residence of droplets in downstream direction in the experiment. This decrease and, consequently, the dilution of the spray in the far downstream region can only be partially reproduced in the simulation. Once a droplet is injected, it exists continuously and unaltered. Effects such as secondary break up or evaporation, which can reduce the size of a droplet, are currently not considered in the simulation. In the experiment, on the other hand, the decrease of droplet diameter leads to an increasing sensitivity to velocity fluctuations. As a result, the effect of the turbulent dispersion on the droplet trajectories is increasingly underestimated in the simulation. Furthermore, the computational mesh becomes coarser in downstream direction, which leads to a damping of the effect of turbulence.
In Fig. 14 the probability density functions for the axial droplet position and the axial velocity in streamwise direction as well as for the droplet volume are shown. In each case, the distributions of the simulation S5 and the experiment are compared with each other. To obtain these distributions, all droplets which are located in the range up to 12 mm downstream the atomizing edge are recorded. The axial position of droplet injection is defined by means of the parameters z = 4.5 mm and z = 1.5 mm . These two parameters enable the adjustment of the spray position in streamwise direction. The range of the droplet injection location is estimated via z ± 2 z to be 1.5 mm <ẑ < 7.5 mm . The investigation of the position PDF indicates that the distribution in this range increases linearly and matches the experimental curve almost perfectly. Further downstream, an agreement is not possible, since the S5 function declines. Small droplets are accelerated faster in the gas phase due to their lower relaxation time than large droplets leading to an axial dispersion of the droplet accumulation. Secondary breakup in the experiments is probably another cause of the discrepancy.
Regarding the velocity PDF, a distinction is made between considering all droplets or just those droplets with a diameter greater than 25 μm . In the experiment, only droplets with D > 25 μm were detectable by PLTV diagnostic. As the graph shows, not taking into account the very small droplets ( D < 25 μm ) gives a good agreement with the experiment. Additionally, the initial Gaussian distribution of the axial velocity w in streamwise direction, as defined by the parameters w = 6 m∕s and w = 5 m∕s , of the S5 simulation is displayed. It is evident that in the near field the prescribed initial distribution results into a larger mean value and more pronounced spread.
For the volume PDF, the volume of each droplet is calculated through its diameter. The sum of the volume of all droplets for a certain diameter range together with the total volume of the near field spray yields the probability distribution. The PDF in the simulation shows a steeper increase and shifted decrease to smaller diameters. Consequently, the Sauter mean diameter of the simulation of 145 μm is below the experimental value of 175 μm ( −17% ). This deviation is the result of the imposed droplet starting velocity (see simulations S3 and S4). The initial velocity parameters are set for simulation S5 with the objective to match the PDF of experimental velocity. Hence, the good agreement of the Sauter mean diameter with the experiment in the simulations S1 and S2 cannot be achieved. At this point it is obvious that not all spray properties can be optimized at the same time using the present set of starting parameters due to contradicting dependencies.
The interdependencies between the droplet diameter D and axial velocity w d in the prefilmer near field for simulation S5 and the experiment are depicted in Fig. 15. The two distributions show that the correlations are similar. The maximum velocity of the droplets of a certain size increases with a smaller diameter because of lower inertia. The highest probabilities in the histograms are found for diameters smaller than 100 μm . The lower percentage of large droplets in the simulation (cf. Fig. 14, bottom) is reflected by a narrow histogram for large diameters. Moreover, the deviating velocities for small diameters between simulation and experiment can be attributed to uncertainties of PLTV measurements.
In the case of interdependencies between the axial droplet velocity w d and axial position z d as plotted in Fig. 16, qualitative similarities between S5 and the experiment can be observed. The velocity of the droplets increases in axial direction due to the momentum exchange with the gas phase. As a result, the smaller droplets are accelerated faster and reach higher velocities at the end of the near field. On average, the velocity as a function of the axial coordinate is reproduced by the simulation.

Conclusions
In this work, Euler-Lagrangian simulations of the fuel spray evolution using the primary atomization model PAMELA [15] have been presented. The model initializes the actual droplet diameters based on predefined drop size distributions and the local gas velocity. As sample geometry a planar airblast prefilmer is investigated. Experimental data [6,21,22] of the gas phase and the liquid spray is used for comparison to the predictions.
For the pure gas phase, the measured and simulated velocity profiles are in good agreement. The profiles of the dominant axial component are reproduced almost perfectly. This velocity component affects the droplet transport most significantly and is decisive for the acceleration of the droplets in flow direction.
Starting from the predicted gas flow, the droplets are injected into the flow field at discrete locations in the region close to the atomizing edge. In an extended version of the PAMELA model, additional droplet starting conditions accounting for the variation of the initial position and velocity are defined. The respective starting parameters for the axial and vertical component are statistically independent and normally distributed. By imposing such statistics to the droplet starting conditions, the essential effects on the spray are identified:  • Reduction of the Sauter mean diameter of the droplet collective in the near field by increasing the axial starting velocity. • Increase of the spray width by imposing a stronger dispersion of the vertical starting position. • Increase of the spray width and the spray angle by imposing an increased vertical starting velocity.
With the parameters of the axial starting position and the axial starting velocity, the density functions of the droplets in the near field can also be influenced and controlled. However, due to contradicting dependencies, it is not possible to optimize all spray properties at the same time with the initial parameters. The present study demonstrates that an optimal definition of droplet starting conditions can be derived which matches the experimental measurements of the atomization at the planar prefilmer. This highlights the potential of the PAMELA model for predicting the primary breakup characteristics in Euler-Lagrangian simulations. Discrepancies to experimental data are probably due to neglecting secondary breakup and evaporation in the simulations.
With the proposed extension of the PAMELA model, it is possible to vary and calibrate initial droplet properties by means of normal distributions to take stochastic breakup mechanisms into account. Furthermore, the requirements for precise gas flow prediction close to the atomizing edge are reduced since the droplet parameters are calibrated in order to match experimental findings. Currently, the starting conditions of the particles are treated independently. In future work, interdependencies between droplet diameter, position and velocity will be modeled by means of multivariate statistical methods.

Appendix A. LES Index of Resolution Quality of Celik
Celik et al. [28] proposed a quality assessment measure as an indicator to evaluate the resolution of LES calculation. This index of quality is based on the mesh cell size h relative to the Kolmogorov length scale k representing the smallest turbulent scales: Assuming that an index of quality greater than 80% is considered as a good LES and h ≅ 25 k , this yields = 0.05 and m ≅ 0.5.
In Fig. 17 the contour plot of the index of quality IQ is shown for the gas simulation. Furthermore, the topology of the tetrahedral mesh with 5 million elements used in the present work is pictured. Since values above 0.8 indicate good quality and sufficient mesh resolution, it can be stated that in the prefilmer far field downstream the atomizing edge the mesh meets the resolution requirements. In this area, the statistical analysis of the spray simulations is performed.

Appendix B. Correlation for primary droplet velocity of Gepperth
Gepperth [35] derived a correlation for the velocity of droplets after primary breakup in a prefilming airblast configuration using shadowgraphy technique: During the experimental campaign, the mean air velocity was varied by a factor of 4 from 20 m/s up to 90 m/s and the air pressure was increased up to 7 bar. The following non-dimensional quantities are used in the correlation. Reynolds number based on the boundary layer thickness at the end of the prefilmer: Weber number based on the boundary layer thickness at the end of the prefilmer: As characteristic length scale, the thickness of the turbulent boundary layer that develops along the prefilmer is applied as expressed by White [36]: .