A Numerical Study of Aero Engine Sub-idle Operation: From a Realistic Representation of Spray Injection to Detailed Chemistry LES-CMC

High altitude relight is a matter of increasing importance for aero engine manufacturers, in which combustion plays literally a vital role. In this paper we want to evaluate the predictive capability of a combined Smoothed Particle Hydrodynamics (SPH) and Large Eddy Simulation with Conditional Moment Closure (LES-CMC) approach for a spray combustion process at these extreme conditions. The focus is on the SPH modelling of the kerosene primary atomization, the extraction of realistic spray boundary conditions for LES-CMC and the effect of the spray on combustion. Interestingly, it will be demonstrated that the fragment size distributions resulting from the airblast atomization are characterized by bimodal behaviour during the relight process and that small and large fragments differ significantly in their dynamical behavior. This is shown to affect the combustion in the Central Recirculation Zone (CRZ). Very large fragments are even able to supersede the flame from the CRZ, such that endothermic pyrolysis becomes dominant, but simultaneously essential to sustain and stabilize the remaining flame with reactive pyrolysis species. The study proves the ability of our methodology for extreme operating conditions, in which experimental insights are hardly possible.


Introduction
Although an in-flight flameout in civil aero engines occurs extremely rarely, it is nevertheless a safety topic of significant importance. As the volume of a combustor needs to be defined early on the design process, it is crucial to study and ensure the high-altitude re-ignition capability of the engine later on. This confirmation becomes especially relevant for modern aero engines, when higher efficiencies are pursued by higher bypass ratios implying increased fan diameters and, hence, rotational inertia. In this case, successful reignition and continuous, reliable heat release inside the combustor during the spool-up of the engine become a vital issue. As set forth by safety regulations, engine manufacturers Extended author information available on the last page of the article 1 3 have to ensure that a relight within a certain operating range of calibrated airspeed and flight altitude is possible ( Fig. 1 Left) (Read 2008;Mosbach et al. 2010). These ranges are termed as relight envelopes and for the engine certification three different relight procedures are distinguished (Zachos 2010;Klinger et al. 2011;Denton et al. 2018). Certainly, all are relevant from an engineering point of view, but in this study we will focus on the windmill relight envelope and in particular the circled upper left corner in Fig. 1 (Left).
These operating conditions are one of the most challenging for the spray combustion process inside the combustor. The thermodynamic conditions at such flight altitudes, of up to 10 km (Read 2008;Mosbach et al. 2010;Denton et al. 2018), are rather extreme. After passing the compressor, the air is still characterized by sub-atmospheric pressure and density, as well as temperatures below 0 • C (Mosbach et al. 2010) in the case of an extinguished flame inside the combustor. Hence, the liquid fuel atomization, evaporation, mixing of air and kerosene vapor, and the ignition chemistry are strongly compromised (Mosbach et al. 2010). Even if the first ignition is successful, a stable reliable heat release inside the combustor during the spool-up of the engine back to idle, termed as pull-away ( Fig. 1 Right), is far beyond a matter of course. Considering the fact that the pull-away can take up to 90 s (Klinger et al. 2011), a fundamental understanding of this complex spray combustion process seems indispensable.
The goal of this paper is to evaluate the predictive capabilities of state-of-the-art CFD methods for this transient spray combustion process, simultaneously pursuing a deeper understanding of the process itself. Therefore, three operating conditions representative of the pull-away process (OP1, OP2 and OP3 as shown in Fig. 1 Lind et al. 2020) and Large Eddy Simulation with Conditional Moment Closure (Klimenko and Bilger 1999;Navarro-Martinez et al. 2005;Mortensen and Bilger 2009;Giusti et al. 2016;Giusti and Mastorakos 2017;Mesquita et al. 2022) approach will be explored. First, the SPH modelling of the kerosene primary atomization and the extraction of realistic spray boundary conditions is presented. Then, the spray combustion process at high altitude conditions will be predicted by means of a LES-CMC approach, in which the spray results will serve as crucial data input.
The structure of the paper is as follows: In Sect. 2 the methodology will be presented, describing the geometry, operating conditions studied, the SPH and CFD methods. In Sect. 3, the results from the flow aerodynamics are presented, as they are an input for the  Zachos (2010). The circled upper left corner highlights the region of interest for this study. Right: Pull-away transient for the circled upper left corner of the windmill relight envelope. With increasing engine load, namely from the first ignition up to idle, most process relevant quantities are characterized by a monotonous behavior in time 1 3 SPH simulations. Then, two dimensional SPH simulations will be presented and critically discussed, including one comparative three dimensional validation run. Finally, the LES-CMC reacting simulations are presented, discussing the influence of the spray and operating conditions on the flame structure. The work will be concluded in Sect. 4.

Geometry
In this work we study a single sector of a developmental Rolls-Royce plc. Rich-Quench-Lean (RQL) industrial combustor. The full geometry includes the pre-diffuser, combustion chamber, and inner and outer annuli to correctly capture the air split between injector, dilution and cooling. We also explicitly model the air film cooling of the combustion chamber walls, to capture its influence on the flame and flow near the wall. Concerning the air flow split, air-to-fuel ratio and residence time, they are all representative of the operation of a RQL combustor in high-altitude relight conditions, when not otherwise specified. The fuel injector relies on the airblast strategy, using multiple concentric swirling flows (with swirl number S n > 1 ) to atomize the fuel and to generate a recirculation zone in the combustion chamber.

Windmill Relight Operating Conditions
This work studies three quasi-stationary operating points (OP1, OP2 and OP3), each representative of different instants of the pull-away transient as depicted in Fig. 1 (Right). The most important operating parameters, namely the air mass flow rate ṁ air relative to OP3, the thermodynamic pressure of air p air and its temperature T air are summarized in Table 1. The consecutive operating points are on a high-level characterized by increased engine load. In more detail, they represent the following conditions: • OP1 is representative of sub-idle high altitude relight immediately before the initiation of the flame. Therefore, the compressor is still in windmilling operation, with the air mass flow rate being significantly reduced compared to idle, which is the nominal condition of OP3. It amounts only 11 % of ṁ air,OP3 . The inlet air is still cold and of subatmospheric pressure. • OP2 corresponds roughly to the half-way through the pull-away transient. Conditions have improved, as the compressor has started working, but nominal idle is not yet reached. The air mass flow rate is 29 % of ṁ air,OP3 and the pressure is still sub-atmospheric. • OP3 represents the end of the pull-away, which, in practice, means that the combustor is operating close to idle nominal conditions.
The fuel mass flow rates ṁ fuel at each operating point are implicitly given by the air-fuel ratio AFR, which is increasing during pull-away. It is interesting to note that ṁ fuel is identical for OP1 and OP2 as sketched in Fig. 1. Moreover, we specify a kerosene temperature of T fuel = 353 K for all operating points, assuming the engine's fuel system thermal inertia being very large. Consequently, the surface tension, the density and the dynamic viscosity of the liquid kerosene are constant in this study at = 0.0186 N∕m , fuel = 757.3 kg∕m 3 and fuel = 0.658 mPas according to Rachner's work (1998).

Motivation
A priori, it is very likely that the spray quality at the considered relight conditions will be affected. The liquid fuel atomization in aero engines is typically realized by airblast nozzles, which accomplish the surface enlargement of the liquid kerosene by means of the air flow momentum (Lefebvre 1980;Lefebvre and McDonell 2017). Although this concept produces excellent sprays at high load, strong deviations from the required spray statistics can be anticipated for relight conditions due to the hampered compressor action, leading to a reduced air momentum (Beck et al. 1991). Spray quantification mostly relies on statistical methods, in light of the fact that a complete deterministic description of myriads of liquid fragments is, up to the present day, elusive. Still, most of this statistical knowledge is purely empirical and relies on complex experimental setups, which date back to the works of Lefebvre and coworkers (Lefebvre 1980;Rizkalla and Lefebvre 1975;Lorenzetto and Lefebvre 1977;El-Shanawany and Lefebvre 1980;Rizk and Lefebvre 1982;Beck et al. 1991). The few fundamental and general frameworks developed in the past to describe spray statistics are asymptotic theories (Cheng and Redner 1990;Gorokhovski and Saveliev 2003;Villermaux 2007;Gorokhovski and Saveliev 2008;Villermaux 2020), which only hold for very simplistic conditions. Therefore, their general application to aero engine combustors, representing strong non-equilibrium reactors, is often too restrictive. This scientific deficit motivated researchers in the last decades to conduct high-fidelity, scale-resolving simulations of the first fragmentation period, called primary atomization, in order to develop a fundamental theory for the spray generation process (Gorokhovski and Herrmann 2008). Especially for airblast atomization, the works presented in Refs. (Sauer et al. 2014(Sauer et al. , 2016Warncke et al. 2017;Bilger and Cant 2018;Braun et al. 2019;Warncke et al. 2020;Wetherell et al. 2020;Carmona et al. 2021;Dauch et al. 2021;Mukundan et al. 2022;Palanti et al. 2022) had been a major breakthrough, with different authors focusing on different methods to unravel the physical foundations of the process. While most authors focus on purely Eulerian methods with interface tracking, such as the Volume of Fluid (VoF) method (Sauer et al. 2014(Sauer et al. , 2016Warncke et al. 2017Warncke et al. , 2020Carmona et al. 2021;Palanti et al. 2022), the Robust Conservative Level-Set (RCLS) method (Bilger and Cant 2018) and Coupled Level Set Volume of Fluid (CLSVoF) simulations (Wetherell et al. 2020;Mukundan et al. 2022), some authors employ a fully Lagrangian Smoothed Particle Hydrodynamics (SPH) approach Dauch et al. 2021). All these methods have their own advantages and drawbacks, but generally they are capable to provide the expected detailed insights into the spray formation process. However, detailed predictions share the conclusion that statistically converged results are nowadays not in reach, due to the computational cost of such simulations. Only a few highly resolved disintegration events can be captured and in most cases for a restricted domain size. The only exception is the full atomizer simulation presented by Warncke et al. (2020).
Due to these computational limitations, we will present a sophisticated two dimensional modelling approach in this work, in order to study the relight atomization problem. Such two dimensional modelling approaches were already employed by Dauch et al. (2017), by Holz et al. (2019) and Mingalev et al. (2020), concluding that the main spray characteristics can still be captured despite neglecting the flow instabilities in the third spatial dimension. Our approach can be interpreted as an extension of the first ideas presented in the work of Dauch et al. (2017) and accordingly solves the resulting multiphase flow problem by means of SPH. Finally, we will extract qualitative representations of the spray statistics, which will serve as realistic spray boundary conditions for the LES-CMC simulations. In most cases, these spray boundary conditions are just guessed. Hence, the statistical estimates from two dimensional simulations constitute a significant improvement.

Definition of Primary Atomization Domain
As explained in the last section, resolved, statistically converged three dimensional simulations of the primary atomization process are nowadays not in reach. Hence, we propose a two dimensional modelling of the primary atomization process, which, however, requires a careful domain definition in order to capture the local flow characteristics relevant for the disintegration process. The aim of this section is to describe the workflow for the two dimensional domain definition as an extension of the work of Dauch et al. (2017) with its features and limitations. It is schematically illustrated in Fig. 2, which also introduces the cylindrical coordinates used in this work, with r = 0 corresponding to the central axis of the atomizer.
The key elements of the workflow are 1. To focus only on the region where primary atomization takes place. 2. To neglect circumferential inhomogeneities instead aiming at the circumferential average. 3. To simultaneously incorporate the effect of the Central Recirculation Zone (CRZ) created by the swirling flow of the airblast atomizer despite neglecting the circumferential flow component.
All these criteria can only be met if at least the temporally averaged aerodynamic flow field inside the airblast atomizer is known a priori, described by its velocity v ∈ ℝ 3 , density ∈ ℝ and pressure p ∈ ℝ . This flow field can be provided either by an experiment or a CFD simulation. In general it is inhomogeneous in circumferential direction. Having this information as input, the whole atomizer can be reduced to a significantly smaller SPH modelling domain as highlighted in red in Fig. 2. However, this requires to introduce a proper circumferential averaging. We consider a density weighted average on small annuli F(x, r) ∶= [r − Δr 2 , r + Δr 2 ) × [0, 2 ) with width Δr as appropriate, since it represents the mass flow rate exactly. The latter follows directly from the x-component of the velocity field v in Eq. (1). Eventually, the circumferential average of the velocity field is defined by 1 3 with dA denoting a differential surface element. Using this definition, it is straightforward to extract inlet velocity profiles in the inner and outer air duct of the atomizer for the SPH modelling domain. The final position of these inlets is arbitrary, but should be sufficiently upstream of the prefilmer to ensure that turbulent fluctuations in the air flow can develop from the static inlet boundary conditions imposed. These fluctuations are likely to be important for the disintegration process.
In order to incorporate the averaged effect of the CRZ in the two dimensional model, we use the axisymmetrical concept of the Stokes streamfunction Ψ ∈ ℝ (Lamb 2009). Though this concept cannot be applied to the inhomogeneous velocity field v , the circumferentially averaged velocity field v is axisymmetrical by nature. From the exact differential of Ψ , one can derive the expression with Ψ(x, r) = const representing streamlines. If one chooses two streamlines such that the primary atomization region is embedded between the inner and outer air duct, the radial expansion of the CRZ is inherently captured in the SPH modelling domain (see Fig. 2).
For the remaining outlet boundary we choose a r = const line sufficiently far downstream the primary atomization region. This is motivated by ensuring that the prescribed boundary conditions at the outlet do not affect the first disintegration process. Fig. 2 Schematic of the workflow for the extraction of the SPH modelling domain. Based on three dimensional aerodynamic flow field inside an airblast atomizer, a proper circumferential average and the concept of the Stokes Streamline, a two dimensional primary atomization domain can be constructed. It also includes the flow deflection caused by the CRZ Nevertheless, it should be stressed that the chosen streamlines are converging, finally causing a nonphysical acceleration of the flow towards the outlet. In contrast, in a three dimensional model, the flow between the corresponding stream surfaces would decelerate for increasing r as the lateral surface would increase. Hence, the post-processing must be focused on the region close to the prefilmer.
All in all, the procedure results in the two dimensional SPH modelling domain as depicted in Fig. 2 with the kerosene film being initialized as a block. This setup will be mostly considered in this study. However, we want to emphasize that it can be easily extended to a three dimensional model with periodic boundary conditions in circumferential direction. We will use such a three dimensional 20 • sector model (Fig. 3) to demonstrate that circumferential instabilities do not alter the main characteristics of the fragment size distribution for OP1 in Table 1.

SPH Modelling
As this work focuses on the prediction of the primary atomization by means of the SPH method [Refs. Monaghan (2005) Contrary to standard Eulerian methods, the idea of SPH is to decompose the fluid domain into N Lagrangian constant mass M i particles with i ∈ {1, ..., N} and size Δl . Hence, the particles move with the flow. Then, based on the momentum balance for each individual particle i with its N ngb neighbors, namely and the kinematic condition  Table 1 the trajectories of the particles can be computed. Subsequently, the approximation of the individual terms in Eq.
(3) will be shortly explained. The pressure gradient term in Eq.
(3) will be modelled by an approach, which for example is presented in the work of Colagrossi and Landrini (2003). It reads in which the individual particle densities are estimated by usage of a quintic B-spline kernel W h,ij (Dehnen and Aly 2012) with smoothing length h = Δl . The latter implies a kernel radius of R K = 3Δl used throughout this work. The expression for the individual densities is given by and appropriate for multiphase flows (Español and Revenga 2003;. Moreover, the term in Eq. (5) requires the individual pressures on the particles. Stating that the fluid flow behaves barotropic, the following equation of state (EOS) will be utilized, which is often termed as Cole equation (1948), reading The numerical parameters , ref , p ref and c s describe a compressibility parameter, the reference density of the fluid, a reference pressure and an artificial speed of sound, respectively. Despite the parameter , which for all cases is air = 1 and fuel = 7 , the remaining parameters depend on the considered case and the fluid phase. The reference densities are defined by the thermodynamic conditions in Sect. 2.2 for each phase. For the reference pressure in the two dimensional cases p ref = is specified, which represents a good compromise between numerical dissipation introduced by the pressure term and stability (Price 2012;Colagrossi et al. 2012 (3) we will follow the work of Szewc et al. (2012), which utilizes the concepts introduced by Brookshaw (1985), Cleary and Monaghan (1999) and . The term is defined as with representing the kinematic viscosity, n ∈ {2, 3} the dimension of the problem and = 0.01h in order to prevent singularities. The surface tension based acceleration a ,i in Eq.
(3) is modelled as a continuum surface force (CSF) according to Brackbill et al. (1992) and was adapted for SPH by Adami et al. (2010) such that density ratios in the order of 1000 can be captured. For details, we refer to the original works.
Lastly, the acceleration term a ,i in Eq. (3) is exclusively considered for the two dimensional SPH simulations to mimic the centrifugal forces of the swirling flow. This model was first introduced in the work of Dauch et al. (2016) and validated later by the same author (Dauch et al. 2017). The key idea is to conserve angular momentum for each particle along its trajectory starting from the inlet. Knowing the circumferentially averaged swirl component of the particle i at the inlet as function of the radius, i.e. v ,i,in (r i,in ) , as explained in Sect. 2.3.2, the centrifugal force can be easily computed as a function of r i only. Since angular momentum conservation implies v (r)r = const , one finds the following expression from the definition of the centrifugal acceleration with the unit vector in radial direction e r .
The above set of equations is solved with the proprietary in-house SPH code turboSPH, which was developed for the prediction of primary atomization. Details regarding the modified Velocity-Verlet time integration scheme, as well as for the boundary conditions are available in the work of Chaussonnet et al. (2020). However, we want to note that the boundary conditions prescribed on the wall defined by the streamlines (see Sect. 2.3.2) must be slip boundary conditions contrary to the geometrical walls on which the usual noslip boundary conditions hold.

Resolution, Statistical Methods and Convergence
Before presenting the statistical insights gained in our analysis, it is of paramount importance to discuss the formal aspects of spatial resolution, temporal convergence as well as how the statistical data was acquired. This is inevitable to confirm the reliability of the statistical results.
Concerning the spatial resolution, all simulations presented in the following were performed with a mean particle distance of Δl = 10 μm and a kernel diameter of D K = 6Δl = 60 μm . Effectively, this mean particle distance results into approx. 1.3 Mio. particles in the two dimensional setups, whereas approx. 230 Mio. particles are required for the corresponding three dimensional setup. The chosen resolution was first heuristically determined for OP1 performing comparative two dimensional simulations with three different mean particle distances of Δl ∈ {10, 20, 40} μm . From the results it could be concluded that only the simulation run with Δl = 10 μm was sufficient to capture the sheet breakup mode. The latter is evident at all operating points and schematically illustrated in Fig. 4 for OP1. In case of coarser particles, the development of these characteristic sheets was strongly inhibited. However, one can also analytically underpin that this resolution is required to resolve the crucial surface tension effects. Therefore, the size of the smallest, spherical structures that emerge at the end of the atomization process must be ascertained. Assuming a mechanical equilibrium, in which the created droplets are balanced by the dynamical pressure Δp dyn and the capillary pressure Δp , one finds an estimate for their equilibrium diameter D eq . It reads with the characteristic velocities of the air flow v air and the droplet v D . It is straightforward to show that this equilibrium corresponds to a subcritical Weber number of representing a stable but probably oscillating droplet for Ohnesorge number Oh < 0.1 (Guildenbecher et al. 2009). In order to obtain a conservative estimate for the size of the smallest structures according to Eq. (10), the relative velocity v air − v D should be as large as possible. The remaining values for and air follow from Sect. 2.2. Since in airblast atomizers only a very small part of the kinetic energy of the air flow is converted into surface energy, i.e. the atomization efficiency is ≤ 1 % (Jedelsky and Jicha 2014), it is reasonable to assume that the characteristic velocity of the air flow is unchanged during the atomization process. Defining v air as the maximum velocity magnitude from the velocity profiles prescribed at the inlets (Sect. 2.3.2) and setting v D = 0 , as the most extreme case, one finds the ratio which was verified to hold for all operating conditions in Table 1. Hence, even the surface tension effects of the smallest structures will be resolved by several SPH particles and is in agreement with Fig. 4.
Nevertheless, the analytical procedure solely justifies that surface tension effects can be properly resolved. Certainly, it must be also addressed whether the resolution is sufficient to account for velocity fluctuations which initiate the primary atomization process. In other words, it has to be clarified if and how well subsonic turbulence can be resolved with the SPH method and the chosen spatial resolution. There exists eligible criticism in the literature that SPH struggles to capture turbulence statistics well, namely the work of Bauer and (10)  Springel (2012). However, we could recently demonstrate that SPH can be interpreted as a discrete Lagrangian LES method using a coarse-graining perspective (Okraschevski et al. 2021a(Okraschevski et al. , 2021b(Okraschevski et al. , 2022. This implies that second-order statistics of turbulent velocity fluctuations can be captured in the best case up to the kernel diameter D K = 2R K , which would correspond to D K = 6Δl = 60 μm for the parameters chosen above. Our experience shows that this is sufficient to capture the primary atomization characteristics Dauch et al. 2021;Chaussonnet et al. 2020). Moreover, it is important to discuss how the statistical data was acquired. Therefore, it must be specified how the liquid fragments are counted and which metrics will be utilized to describe the statistical behavior. Initially, we have to define the term fragment itself. Subsequently, a fragment is defined as an agglomerate which comprises at least four SPH particles in the two dimensional cases or eight particles in the three dimensional case. Smaller fragments are interpreted as numerical artefacts and discarded. In principle, the liquid fragments produced during the process can be either counted inside the whole SPH modelling domain or passing a properly placed sampling plane for each time step. This is identical to experimental investigations and comes with similar problems (Tropea 2011). One of the main issues of the domain based method is that fragments are counted multiple times causing a systematic bias of the statistics. This is caused by the fact that smaller fragments move with different characteristic velocities compared to large fragments, which can be empirically accounted for through velocity weighting, however, it lacks rigorous justification. Hence, we will evaluate the statistics using a sampling plane (Fig. 5).
Nevertheless, care must be taken in terms of the proper placement of this sampling plane. On the one hand, as explained in Sect. 2.3.2, our two dimensional SPH domain definition procedure becomes inaccurate in the vicinity of the outlet, which requires that the position is as close as possible to the prefilmer. On the other hand it must also be ensured that even the longest lamellae and their primary atomization products will be counted in the statistics. Therefore, the temporal evolution of the axial lamella length L Lam for OP1 was evaluated. This is depicted in Fig. 6a. Except for some rare single events, the strongly fluctuating lamella length is mostly below L Lam < 3 mm . Consequently, we placed the measurement plane 3 mm downstream the atomization edge, as highlighted in red in Fig. 5. This placement will be used subsequently also for OP2 and OP3 as the lamella length decreases with increasing engine load during the pullaway (not shown herein).
Having defined the method of fragment sampling, we can proceed with the description of the statistical metrics. Therefore, random variables z characterizing the liquid fragments have to be chosen first. In this work, z will either be 1. The equivalent circular area diameter D ∶= √ 4A D ∕ in two dimensions, with A D as area of the liquid fragment, or the equivalent spherical diameter D ∶= 3 √ 6V D ∕ in three dimensions, with V D as volume of the liquid fragment, 2. The radial position of the center of mass r, 3. The axial velocity of the center of mass v x or 4. The radial velocity of the center of mass v r .
Then, one metric which will be analyzed is the univariate volume based probability density function f(z) (pdf) of the variable z. Restricting ourselves to N bin = 10 equidistant bins with width Δz = (z max − z min )∕N bin , in which z min , z max describe values close to the minimum and maximum found in all liquid fragment collectives considered herein, and bin centers z k = z min + (k − 1∕2)Δz, k ∈ {1, ..., N bin } , the volume based pdf reads In Eq. (13), the upper index N D is the overall count of liquid fragments in the collective and the quantity ΔV k represents the overall volume collected in bin z k . Whereas the fragment volume V D,j of fragment j is directly accessible in the three dimensional case, in the two dimensional cases the assumption of sphericity will be used to estimate V D,j based on D .
(13) We will discuss this choice in Sect. 3.2.3. The second metric that will be considered is the bivariate scatter plot of different pairs of random variables as listed above. These scatter plots qualitatively represent the bivariate statistics of the fragment collective. A final important aspect concerns the temporal convergence of the statistics. Therefore, we have listed relevant metrics as the simulated time ranges T, the number of liquid fragments N D counted and the dominant frequencies f max of the smoothed axial lamella length spectra resulting from a discrete Fourier transform (Fig. 6b) in Table 2 for all simulations. Only for the three dimensional simulation the frequency could not be determined as only one breakup event was captured. The dominant frequency is roughly indicative for the minimal number of breakup events N B = Tf max captured in the simulations, although faster, less important modes (see Fig. 6b) will likely cause additional breakup events. Evidently, the two dimensional simulations contain at least an order of magnitude more liquid fragments than the three dimensional run and also significantly more breakup events. Consequently, their statistical reliability exceeds the one from the three dimensional simulation. However as detailed in Sect. 3, even for the two dimensional runs, only a few of the large and rare liquid fragments are captured, which might indicate some lack of convergence for these structures.

CFD Methods
As discussed in Sect. 2.3.2, the SPH simulations require the aerodynamic flow field in the combustor as an input to perform the primary atomization simulations of the fuel. Therefore, the first step is to predict the combustor aerodynamics. This is done in this study by conducting RANS simulations for each of the three operating points, as only the averaged flow is necessary. Afterwards, SPH spray statistics are acquired, which are used in subsequent LES-CMC simulations as spray starting conditions.
The Rolls-Royce proprietary finite-volume code PRECISE-UNS (Anand et al. 2013) is used for all these CFD simulations (RANS and LES), solving the low Mach number approximation for the gas-phase. A second-order scheme is used for spatial discretization, while, for time derivatives in LES, a second-order implicit backward scheme is used. Moreover, for the LES, the Vreman model (2004) was used to close the sub-grid scale stress tensor, as it is able to provide a vanishing turbulent viscosity near to the walls. The utilized mesh is hexa-dominant and unstructured, composed of approximately 32 million cells, and refined inside the injector's passage and the swirler. The time step used is Δt = 0.1 μs , ensuring CFL number under 1 in all of the domain.
One essential element to assure the accuracy of each of these CFD simulations, and the subsequent SPH simulations, is the knowledge of the air flow at the combustor inlet for each operating point. A major challenge of studying the phenomena occurring in the combustor during pull-away is the limited availability of data. To address this limitation  (2021) at Cranfield University (CU). These simulations provided realistic velocity profiles to be imposed at the inlet for each case. The imposed air inlet boundary conditions for the RANS correspond to the time averaged data from the compressor simulations. These boundary conditions were compared against air mass flow rate rig data and validated for each case. For the outlet boundary conditions at the combustor's exit, typical conditions are used.
Once the spray statistics are acquired from the SPH runs, reacting LES-CMC simulations (Giusti et al. 2016;Giusti and Mastorakos 2017;Zhang et al. 2019;Foale et al. 2019) using an Euler-Lagrange framework, similar to Thari et al. (2021), are performed in order to show the importance of realistic spray starting conditions. In the reacting LES performed for the full combustor, fragments are directly injected as spherical parcels with the spray statistics imposed from the SPH simulations of each case, namely diameter, velocity, injection angle. Evaporation is modelled with the Abramzon-Sirignano model (1989). Secondary breakup is considered through the Schmehl model (2002). The interaction of fragments with the wall is modelled as elastic rebound.
To model combustion, we employ the Conditional Moment Closure (CMC) combustion model (Zhang et al. 2019;Giusti and Mastorakos 2017), which is adapted to tackle the difficult non-premixed conditions of these cases. The CMC model solves the conditionally filtered species mass fraction and the conditionally filtered temperature or enthalpy. The CMC equation for a given reacting scalar Q is defined as where Q =Ỹ | represents the -th species filtered conditional mass fraction and represents the variable sample in mixture fraction space. The quantity represents the sub-grid scales contribution, which is modelled with the gradient assumption and neglects the subgrid conditional joint fluctuations of the droplet evaporation rate and species (Borghesi et al. 2011;Tyliszczak et al. 2014). Further, ũ j | represents the conditional velocity, assumed � u j | =ũ j . The conditional source terms S(Π| ) are related to the spray evaporation, modelled as in Tyliszczak et al. (2014). The quantity ̃ | represents the conditionally filtered chemical source term and is modelled with the first order CMC model that considers ̃ | ≈ (Q 1 , ..., Q n , Q T ) , with n the number of species and Q T ≡ T| as conditionally filtered temperature. The scalar dissipation rate Ñ | in mixture fraction space is closed with the Amplitude Mapping Closure (AMC) (O'Brien and Jiang 1991) In the physical space, the filtered scalar dissipation rate Ñ is calculated from the LES solution, with contributions from resolved and sub-grid fields where D = ∕(̄Sc) is the molecular diffusivity with an assumed Schmidt number of 0.7, t is the sub-grid kinematic viscosity and C N is a model constant, taken as 42.0 (Triantafyllidis et al. 2009).
As the CMC model relies on a mixture fraction approach, the transport of the filtered mixture fraction and its variation are solved with additional equations (Pera et al. 2006;Triantafyllidis and Mastorakos 2010) where Π represents the evaporation source term and the sub-grid term J SGS is modelled by a gradient approach J SGS =̄D t̃∕ x j with turbulent diffusivity D t . The latter is computed from the turbulent viscosity with a turbulent Schmidt number of 0.4 (Kempf et al. 2008;Triantafyllidis et al. 2009;Tyliszczak et al. 2014). The scalar dissipation rate Ñ is modelled as aforementioned and the ̃ Π terms are modelled as � Π =̄sΠ , following Tyliszczak et al. (2014), where ̄s is the saturation mixture fraction. The CMC model is solved on a second, coarser mesh composed of approximately 46,000 cells, discretizing only the flame tube. This results in cells of around 5 mm near the injector, which are progressively coarsened. The coupling between the two solvers is done through the density and the temperature, following the implementation discussed in Garmory and Mastorakos (2015). The operator splitting technique was employed to solve the CMC equations, treating the chemical source terms with the VODPK implicit solver (Brown and Hindmarsh 1989). A first-order upwind scheme is used for convective terms, a second-order one for diffusion terms and a first-order one for time discretization. Filtered Probability Density Functions (FDF) are then used to compute the unconditional value of a given filtered quantity f in physical space from the respective conditional value f | with where P ( ) is the FDF. The modelling used in this work presumes that the FDF has a -function shape calculated from the resolved mixture fraction ̃ and its variance ′′2 . LES-CMC has been already validated for flame ignition and extinction in a range of configurations, which should be consulted for more details on the model and the code [e.g. Zhang et al. (2019); Giusti and Mastorakos (2017); Borghesi et al. (2011)].
The Jet-A Kerosene "A2" from the NJFCP (Colket et al. 2017) is employed as fuel. We apply the detailed high-temperature HyChem mechanism Xu et al. 2018) to model it. The HyChem mechanism considers Jet-A as a single component fuel with the average fuel properties of C 11 H 22 , and comprises 119 species and 843 reactions. A particularly interesting feature of this mechanism for this work is the explicitly modelling of fuel pyrolysis. This is one of the main reasons this mechanism is used in this work, as this capability allows for the evaluation of the influence of fuel pyrolysis on the flame stabilisation under the extreme operating conditions we investigate. The HyChem mechanism uses seven semi-global reactions steps to lump the fuel pyrolysis, with the following decomposition products: H , CH 4 , C 2 H 4 , C 3 H 6 , 1-C 4 H 8 , i-C 4 H 8 , C 6 H 6 , C 6 H 5 CH 3 and the methyl radical CH 3 . The USC Mech-II detailed mechanism (with reactions updated for i-C 4 H 8 ) is then used to oxidise these decomposition products. This numerical setup (LES-CMC approach with HyChem mechanism) has been successfully used in similar applications (Foale et al. 2019;Mesquita et al. 2022), ensuring its use for this case.

3
The reacting simulations of this work are continued after the realisation of the full ignition sequence (Mesquita et al. 2022). After the ignition, which are done under OP1 conditions, the simulation is run for more 15 ms to ensure stabilisation. After this period, the operating conditions are changed again to match OP2, first, and then OP3. In each case 15 ms are run to reach stabilisation before performing the analysis.

Flow Aerodynamics
We present in this section the results of the RANS simulations for the three operating points defined in Sect. 2.2. First, we look at the flow at the inlet and in the pre-diffuser to study the impact of the initial and boundary conditions. Looking at Fig. 7, one can see on the left how the axial velocity profile of the flow coming from the compressor strongly changes as a function of the operating point. At the OP1 sub-idle conditions, the velocity field at the inlet (Fig. 7a Left) is much more irregular even having a low velocity region at its center. This irregular profile affects the flow inside the prediffuser creating flow separation at its top walls (blue region at the centre of Fig. 7a Right). As conditions improve going from OP1 to OP2 (Fig. 7b) and then to OP3 (Fig. 7c), one can see that the inlet profiles become more uniform and the separation region disappears.
Nevertheless, even if the aerodynamics inside the pre-diffuser is strongly affected by the inlet conditions, inside the combustion chamber, the average flow aerodynamics are less affected (Fig. 8). Indeed, the change of the air mass flow rate does not affect the flow topology, all cases presenting a Central Recirculation Zone (CRZ) with the same shape and a Conic Vortex Breakdown (CVB). Furthermore, the change in the incoming flow at the pre-diffusor does not significantly modify the flow split between the primary zone (i.e. the air flow that goes inside the swirler) and the dilution jets (i.e. the flow that goes through the Fig. 7 Axial velocity field inside the pre-diffuser for a OP1, b OP2 and c OP3. The left fields of each image represent the axial velocity at the perpendicular inlet plane indicating the flow entering pre-diffusor, while the right fields represent the longitudinal cut of the pre-diffuser, the flow passing from left to right inner and outer annuli surrounding the combustion chamber). In all cases, approximately 2% of the inlet air mass flow rate goes through the central stage of the swirler, where the air blast atomizer is localized (black arrows in Fig. 8). Therefore, the major impact of the operating conditions on the atomization quality is a consequence of the change in the air mass flow rate arriving at the combustion chamber from the compressor ( ṁ air,OP1 = 0.11ṁ air,OP3 and ṁ air,OP2 = 0.29ṁ air,OP3 ). Since the air split and the percentage of air that reaches the atomizer is constant for all cases, it seems that the velocity profile at the inlet and the flow separation inside the pre-diffuser does not affect the fuel atomization nor the average combustion chamber flow. It is the change in the air momentum that will then affect the quality of the atomization accordingly.
Applying the circumferential average in Eq.
(1) to the RANS solution inside the atomizer in order to extract the Stokes streamlines for the SPH domain definition (Sect. 2.3.2), confirms the observations just discussed from another perspective. In Fig. 9 the streamlines extracted for each OP, starting from the same initial point in the inner and outer air duct, are illustrated. The coordinates are normalized. The dots represent the data points resulting from the post-processing routine (blue for the inner and red for the outer boundary), whereas the black solid lines are least-square fits with the ansatz functions r inner (x) = A ⋅ exp(cx 3 + bx + a) and r outer (x) = cx 2 + bx + a and parameters A, a, b, c ∈ ℝ . Evidently, all resulting streamlines represent the radial deflection caused by the CRZ and they are very similar. With increasing engine load from OP1 to OP2 (cp. Fig. 9a, b) as well as from OP2 to OP3 (cp. Fig. 9b, c) only slight shifts of the streamlines towards larger radii are noticeable. However, the characteristics of this averaged deflection caused by the CRZ is robust during the whole pull-away process despite the very different inlet conditions at the pre-diffusor. Hence, for the definition of the SPH modelling domain according to Sect. 2.3.2, only the streamlines extracted for OP1 (Fig. 9a) will be used.

Two Dimensional Statistical Analysis for OP1
In this section we will analyze the spray statistics for OP1 as explained in Sect. 2.3.4. The univariate volume based pdfs are depicted in Fig. 10 (blue curves) as well as the number of fragments ΔN k contained in each bin (red curves). Since merely for the fragment size distribution a significant, qualitative difference between the volume based pdf and the ΔN k distribution exists, only these differences will be explained. For the remaining volume based pdfs the ΔN k distributions are added for completeness. Complementary, the scatter plots of selected pairs of random variables are illustrated in Fig. 11.
The most interesting observation in the univariate pdfs concerns the fragment size distribution in Fig. 10a, which reveals a bimodal shape with local maxima around D max,1 ≈ 200 m and D max,2 ≈ 400 m . This is physically reasonable despite the fact that most of the fragments are smaller than D max,1 (red curve in Fig. 10a). Since the pdf scales with f k ∼ ΔV k ∼ ΔN k D 3 k according to Eq. (13), even rare liquid fragments with D > D max,1 are of significant importance for the considered fragment statistics. Taking asymptotic theories as reference, it is obvious that the measured fragmentation collective is far from an equilibrium distribution and that the atomizer is not operated in its actual design point at OP1. Usually, one would expect to observe monomodal behavior close to equilibrium Cheng and Redner (1990), Gorokhovski and Saveliev (2003;2008), Villermaux (2007;2020). This is a valuable insight for the spray initialization in the Euler-Lagrange based LES-CMC simulations in Sect. 3.3.
The remaining univariate volume based pdfs are either bell shaped in case of r (Fig. 10b) or sharp peaked for v x and v r (Fig. 10c and d) for the considered binning. It is interesting to note that the peak of the radial position distribution in Fig. 10b is around r max ≈ 7 mm and, hence, above the radius of the prefilmer, i.e. r max > r prefilmer . This reflects the mean effect of the CRZ on the atomization process, which pushes the liquid fragments towards larger radii. The axial velocity distribution in Fig. 10c is mostly characterized by positive values. However, some fragments seem to move against the main flow direction, which is likely caused by aerodynamic fluctuations present in the gaseous phase. The radial velocity distribution behaves very similar being centered around v r,max ≈ 0 and slightly skewed towards positive velocities, which again is due to the CRZ.
The scatter plots in Fig. 11, representing the bivariate statistics between pairs of random variables, support the previous observations, additionally providing further insights into the spray statistics. For visualization purposes the diagrams in Fig. 11 only include D < 250 m and a least-square fit with the ansatz function g(D) = a ⋅ exp(−bD) + c and parameters a, b, c ∈ ℝ (red curve). The latter is indicative for the qualitative trend of statistical correlations. In Fig. 11a the resulting point cloud for the pair {r, D} is shown. Evidently, most fragments are smaller than ∼ 100 m and characterized by strong dispersion regarding the radial position. This dispersion decreases for larger D , which is likely related to the fact that the inertia of the fragment increases, as well as the relaxation time D ∼ D 2 (Balachandar and Eaton 2010). Consequently, smaller fragments do respond faster to aerodynamic fluctuations causing this dispersion. As the Stokes number is proportional to the relaxation time, i.e. St ∼ D (Balachandar and Eaton 2010), we will call this converging trend in the points clouds a Stokes number effect. Moreover, the exponentially decaying red trend line in Fig. 11a reveals that small fragments are on average pushed stronger towards larger radii than large fragments, due to the lower inertia and the effect of the CRZ. Large fragments converge towards a constant value of around r ≈ 7 mm . The latter is in accordance with the bell shaped behavior in the univariate volume based pdf in Fig. 10b with its peak at r max ≈ 7 mm.
Overall, the situation is very similar for the pairs {v x , D} and {v r , D} shown in Fig. 11b and c. Again, small fragments are characterized by more dispersion of the axial and radial velocity than the large fragments, reflecting the Stokes number effect. For the axial velocity in Fig. 11b the red trend line is exponentially decaying, being representative for the averaged correlation. Apparently, small fragments tend to move faster in axial direction than large fragments. The trend line converges to a value slightly below v x ≈ 5 m∕s , consistent with the peak of the univariate pdf in Fig. 10c. Focusing on the pair {v r , D} in Fig. 11c, the only difference is that the red trend line is constant. Its value, which is marginally above zero, is again in accordance with the univariate pdf in Fig. 10d. It is interesting to note that the shape of the point clouds for {v x , D} and {v r , D} , including the Stokes number effect, is a general behavior of all sprays also known from experiments, e.g. Refs. Holz et al. (2019), Soni and Kolhe (2021), Vallon et al. (2021). We interpret this as one of the biggest strengths of our two dimensional SPH modelling procedure, which can reproduce this characteristic behavior at low computational cost.
So far, the topic of validation has not been addressed. Despite the fact that the numerical observations seem to be reasonable up to this point, a validation with experimental results would significantly substantiate our whole modelling approach. However, it is hardly possible to experimentally extract primary atomization statistics from airblast atomizers similar to the one considered, due to the limited optical accessibility (Menon and Ranjan 2016). Nevertheless, more arguments can be gathered to underpin the results described above. Hence, we will explain in Sects. 3.2.2 and 3.2.3 why the bimodality in the fragment size distribution, as one of the most important quantities in atomization, seems reasonable notwithstanding the lacking statistical significance in Fig. 10a for the large fragments. After that the effect of the pull-away process on the spray quality will be analyzed in Sect. 3.2.4.

Lamella-Aero-Interaction for OP1
In this part we will discuss why the observed bimodality in the fragment size distribution in Fig. 10a seems physically sound based on a qualitative analysis of the interaction between the liquid lamella and the aerodynamics. We will term it as lamella-aero-interation. For the analysis, the backward finite-time Lyapunov exponent (FTLE) (Sun et al. 2016;Dauch et al. 2018) is utilized with a integration time of Δt = 0.0005 s in order to qualitatively visualize the vortex dynamics of the air flow. A similar FTLE study for airblast atomization was presented by Dauch et al. (2019). The FTLE fields of three consecutive snapshots, superimposed with liquid kerosene (blue), are depicted in Fig. 12. These snapshots represent a typical disintegration sequence observable in the atomizer studied.
As soon as a lamella or sheet has developed, the FTLE fields reveal that a vortex develops close to it, which is a consequence of the film waves created on the prefilmer (Fig. 12a). These waves act as a fluid dynamic analogy of the backward facing step, eventually leading to a local backflow once the air leaves the wave crest behind. Consecutively, these vortices move further downstream causing the kerosene lamella to be bent towards the outer air duct (Fig. 12b). The upwards bent but still growing lamella effectively results in an increasing blockage of the outer air duct such that the air there is pilled up. As a reaction, the influence of the outer air flow on the lamella grows in time, first pushing the lamella into the main flow direction again but finally also causing an inflation of the lamella up to the point when disintegration takes place. From Fig. 12c it becomes clear that the disintegration is related to chaotic vortex dynamics, but more importantly to a characteristic spray signature. Most fragments created in this primary atomization process are small, numerous and almost spherical, as highlighted by the circles, whereas only a few large, aspherical fragments are created, as indicated by the arrows. These two size classes of fragments are likely to be the reason why the volume based pdf in Fig. 12a is characterized by a bimodal shape. As stated in Sect. 3.2.1, since the pdf cubically scales with the characteristic length scale of the fragment, i.e. f k ∼ ΔV k ∼ ΔN k D 3 k , large fragments become important for the statistics even if they are rarely created.
Finally, we can state that the analysis of the lamella-aero-interaction possibly explains the statistical observations of the fragment size distribution, indicating that the two dimensional modelling pursued herein is in itself reasonable. However, it should be borne in mind that this analysis is not a prove for the observed bimodality in the fragment size distribution, but only a possible explanation in case this bimodality truly exists. Hence, it does not resolve the issue of lacking convergence for the large liquid fragments.

Three Dimensional Analysis for OP1
To ensure that the two dimensional simulation results are qualitatively correct despite neglecting the circumferential flow, a comparative three dimensional simulation for OP1 was conducted. As explained in Sect. 2.3.3, a circumferentially periodic 20 • sector was simulated, which comprises approx. 230 Mio. particles of size Δl = 10 μm . The chosen sector angle is arbitrary, but the focus is on including the swirling flow, causing circumferential instabilities. However, it should be mentioned that the chosen sector will not be sufficient to capture the transient effect of the precessing vortex core (PVC). The simulation was run on the HoreKa supercomputer of KIT on 988 cores for around 1345 h ≈ 2 months with the same domain decomposition strategy as presented by Dauch et al. (2021). This makes sure that the simulation was conducted as efficient as possible with the turboSPH code. Overall, the result represents exactly one disintegration event as listed in Table 2 in Sect. 2.3.4, which simultaneously highlights the biggest drawback of these detailed three dimensional simulations. The corresponding two dimensional model covers a 62 times larger time range and the simulation is a matter of only a few days. Hence, the three dimensional results must be understood as a qualitative plausibility for the actual atomization process. Fig. 12 FTLE fields superimposed with liquid kerosene (blue) in order to visualize the lamella-aero-interaction for a typical desintegration sequence. a Vortex development behind the prefilmer close to the lamella ( t = 0.0844 s ). b Lamella bending towards the outer air duct as consequence of the vortex transport ( t = 0.0848 s ). c Chaotic vortex dynamics after desintegration ( t = 0.0856 s ). The generated spray is characterized by many, almost spherical droplets (circles) and single, large aspherical fragments (arrows) Nevertheless, the data set is sufficient to give an idea which effects can be captured by the two dimensional model and which cannot. Most evidently, in the three dimensional disintegration process a thin liquid lamella is created, which is also characteristic for the two dimensional simulation of OP1. However, as depicted in Fig. 13, the lamella is additionally exposed to circumferential instabilities triggered by the swirling flow as expected. Consequently, it becomes clear that the actual physical transition of the lamella to the final spray will likely be different from the aero-lamella-interaction described in Sect. 3.2.2. However, the fragmentation collective resulting from the primary atomization process (Fig. 14) is still characterized by two size classes, as it was the case in the two dimensional simulation. Mostly, very small, spherical droplets are created, albeit not exclusively. There are also a few larger fragments, which can be observed.
Evaluating the volume based fragment size distribution inside the whole domain for the time step illustrated in Fig. 14, which departs from the two dimensional procedure owing to the lack of disintegration events captured, one finds the result shown in Fig. 15b. A comparison with the two dimensional pdf in Fig. 15a shows that the actual shape of the pdfs is quite different. Please also note the difference in the scales of the y-axes. Still, the pdf in Fig. 15b reflects the two fragment size classes by a bimodal shape. Interestingly, the local maxima are also around D max,1 ≈ 200 μm and D max,2 ≈ 400 μm in accordance with the two dimensional simulation. It is very likely that this is just a statistical coincidence, however, a small chance remains that this is attributable to the volume estimation we introduced for the two dimensional simulation. Therefore, we assumed sphericity in order to define V D,j based on the equivalent circular area diameter D (see Sect. 2.3.4). Figure 14 qualitatively indicates that for the given atomizer in OP1 high level of sphericity prevails independent of the fragment size, which might justify this choice.
We like to conclude this section with the observation that the two dimensional model seems to be able to qualitatively capture a sheet breakup process and the bimodality in the fragment size distribution. Indeed, the actual transition from the lamella to the spray evolves differently due to neglecting the circumferential instabilities. However, we are convinced that the two dimensional model is sufficient to study how the spray statistics is qualitatively affected during pull-away. Hence, in the next section we will proceed with this objective. Fig. 13 Three dimensional kerosene lamella disintegrating into the primary atomization products. Evidently, the process is influenced by circumferential instabilities, which are neglected in the corresponding two dimensional model

Spray Statistics During Pull-Away
In this section, the spray statistics resulting from the two dimensional simulations for the three pull-away conditions described in Sect. 2.2 will be compared. Information about the temporal metrics are listed in Table 2 in Sect. 2.3.4. It is interesting to note that from OP1 to OP3 the fragment count increases as a consequence of the increased aerodynamic load and larger air momentum during pull-away (Sect. 3.1), although the simulated time ranges decrease. The univariate pdfs for D, r and v x , with the same representation as in Sect. 3.2.1, are given in Fig. 16. The first row illustrates the results for OP1, the second for OP2 and the third for OP3. In Fig. 17, the scatter plots for {r, D} , {r, v x } and {r, v r } are juxtaposed. Again, the rows represent the individual operating conditions. The biggest impact of the increased aerodynamic load during the pull-away process is evident in the univariate fragment size distribution in Fig. 16a, d, g. Starting from the bimodal pdf in Fig. 16a at OP1 with large fragments at D max,2 ≈ 400 μm being the most import ones, an increase of the air momentum to OP2 first causes a merging of both local maxima accompanied by a reduction of the biggest fragments observed. The resulting pdf Is is interesting to note that most fragments are qualitatively characterized by a high level of sphericity independent of their size. Moreover, two fragment classes can be distinguished. One class of small, numerous droplets and the other class comprising single, larger fragments  Fig. 16d is a broad and skewed monomodal distribution, still comprising fragments with length scales around D ≈ 400 μm . Further intensification of the aerodynamics to OP3 demonstrates that the focus of the pdf is shifted increasingly to smaller scales, finally creating a distribution, which approaches the ideal shape as reported by asymptotic theories Fig. 16g. This is expected as the atomizer at OP3 is closer to its design point compared to OP1 and OP2, which demonstrates that the two dimensional SPH predictions correctly capture this trend. Nonetheless, it must be acknowledged that a weak bimodality with a second peak at D ≈ 400 μm persists, even if the engine is operated close to idle conditions. This is one of the most important observations of our study and an useful insight for the LES-CMC simulations of the spray combustion process during pull-away in Sect. 3.3. Although, it cannot be completely ruled out that this behaviour might be a consequence of the lacking convergence for the large and rare liquid fragments, the observation seems reasonable in respect of Sects. 3.2.2 and 3.2.3.
The pdfs of the radial position presented in Fig. 16b, e, h are quite robust in their shape during the pull-away process. For all operating conditions a bell shaped pdf with a maximum around r max ≈ 7 mm is present. However, the pdf becomes narrower as the aerodynamic load is increased. This is likely related to the fact that the lamella length decreases on average from OP1 to OP3 (not shown herein), which leads to less radial fragment dispersion considering the lamella-aero-interaction presented in Sect. 3.2.2. For the axial velocity pdf in Fig. 16c, f, i, the influence of the pull-away is also not as pronounced as Fig. 16 Estimated volume based pdfs during pull-away (blue). The first row represents OP1, the second OP2 and the third OP3. a, d, g Fragment size distributions including the number of fragments in each bin (red). b, e, hRadial position distributions. c, f, i Axial velocity distributions for the fragment size pdf. The sharp peaked distribution overall keeps its shape, though it is shifted on average towards larger axial velocities and is characterized by a larger variance. These observations are reasonable, since the averaged air momentum increases from OP1 to OP3 and accordingly the aerodynamic fluctuations in the gaseous phase must grow, eventually transporting momentum fluctuation to the liquid fragments as well. We will not show the results for the radial velocity pdfs as they are similar to the axial velocity component.
The global shape of the scatter plots in Fig. 17, being representative for the bivariate statistics during pull-away, is extremely robust to the increased aerodynamic load. All pairs include the characteristic Stokes number effect as explained in Sect. 3.2.1. Mainly, the results of the univariate pdfs described above are confirmed. Starting with the pair {r, D} in Fig. 17a, d, g, it is obvious that the point clouds approach the red trend lines more closely with increased engine load. Hence, the point clouds become more compact, which is likely related to the mitigation of radial dispersion caused by shorter lamellae at larger aerodynamic load. At all operating conditions, the point clouds approach a value of r ≈ 7 mm for Fig. 17 Scatter plots representing bivariate statistics during pull-away. The red curve is a least-square fit, which indicates the averaged correlation between two random variables. The first row represents OP1, the second OP2 and the third OP3.  Fig. 17b, e, h, the Stokes number effect during pull-away is enhanced, which is also reflected by the stronger curvature of the red trend line. Especially for OP3 in Fig. 17h the axial velocity dispersion of the small fragments is strongly increased compared to Fig. 17b, e. This maybe rooted back to the increase of mean air momentum that causes an intensification of aerodynamics fluctuations, dominantly affecting the small liquid structures. Moreover, the points clouds are shifted towards larger axial velocity values as direct consequence of the higher aerodynamic load, as the asymptotics of the trend lines reveal. Again, these observations confirm the previously discussed coherences for the univariate pdfs in Fig. 16b, e, h. The point clouds for the pair {v r , D} in Fig. 17c, f, i behave almost identical except for the fact that the absolute levels of the red trend lines are smaller in general.

Reacting Flow Analysis
Following the acquisition of spray statistics from the SPH simulations, the results were used to perform reacting LES-CMC simulations as explained in Sect. 2.4. This allows for a characterization of the flame at the three different stages of pull-away investigated herein. The SPH spray statistics are used as starting conditions for the spray in the full combustor CFD, where the liquid-phase is modelled with a Lagrangian discrete particle approach, thus all fragments are modelled as spherical particles. Hence, this section aims to demonstrate the importance of incorporating realistic spray starting conditions to obtain physically sound results. First, OP1 is investigated, where the high-altitude relight of the combustor is simulated. However, the ignition sequence is not discussed in this work but instead in detail in Mesquita et al. (2022). We focus here on the stable flame obtained after the ignition sequence. For the subsequent two operating points of interest (i.e. OP2 and OP3), the operating conditions of the simulation are progressively changed to correspond to OP2 and OP3. All flames are stabilized for 20 ms before being analyzed. Figure 18 shows the flame for all the three conditions. One can see that for all cases the flame has a similar shape with two parts: in the primary zone, the flame is rooted near the injector and mostly stabilized at the exterior limits of the Central Recirculation Zone (CRZ) at the Inner Shear Layer (ISL); in the dilution zone (after the dilution jets), the flame spreads shapeless through the rest of the flame tube. The similarity of flame shape between the three cases is also retrieved on the mixture fraction fields (Fig. 18). Indeed, all three points present the behavior expected for a RQL burner such as this, with the primary zone presenting a rich mixture with the flame burning around the stoichiometric mixture region, while in the subsequent dilution zone the mixture and combustion are mostly lean. However, this similarity in flame shape and mixture fraction field hides a strong difference in flame structure. This difference in flame structure can be visualized looking at the temperature versus mixture fraction plots of Fig. 19 for all the three points. For OP1, various levels of reaction are present between "low" temperature pyrolysis (roughly occurring for temperatures under 1200 K , region under the dotted black line path in Fig. 19a) and a fully-burning flame (red regions in Fig. 19a). As the operating conditions change to OP2 (Fig. 19b), these "low" temperature reacting states reduce and for OP3 the fully-burning state is mostly exclusive found, with only some restricted regions presenting isolated pyrolysis (blue path shown in Fig. 19c). We can see, then, that OP1 is a very spatially heterogeneous cases, chemically-wise, while OP3 is mostly homogeneous. This results confirms the hypothesis of this work that detailed chemistry is relevant for relight conditions (OP1), while it also shows that for OP3 a simpler chemical modelling is very likely enough. This difference in 1 3 flame structure between operating points, marking a very distinct general behavior of the system in each case, is a consequence of the very different operating conditions and spray atomization levels. We further analyze these aspects in the following paragraphs. The Sauter Mean Diameter fields (Fig. 18) show how the different atomization quality in the three operating points (summarized in Fig. 16) translate to different spatial distributions of spray. For the OP1 case, very large fragments are found at the immediate vicinity of the injector. These fragments are mostly injected straight inside the CRZ, where the smallest accumulate and the largest keep travelling to reach the end of the flame tube. We see a drastic difference in terms of fragments size when conditions are changed to OP2, with the diameter of fragments inside the CRZ being roughly halved in comparison to OP1. The fragments' diameter are further reduced moving to OP3, where the fragments have all similar sizes and spread rather evenly inside the primary zone. While the Liquid Volume Fraction plots (Fig. 18) point out to a rather similarly concentrated liquid-phase inside the CRZ, one can see how it also shows the spray angle and penetration increasing progressively, leading to a much better spread of the spray at OP3 conditions. The concentration of cold liquid fuel in the form of large fragments inside the CRZ for the OP1 case results in the strong concentration of cold fuel vapor therein (Fig. 20a). The cold temperature and the high concentration of both liquid and gaseous fuel prevents the mixture from reacting inside the CRZ, requiring the breakup of the fuel by pyrolysis first to allow for combustion. Indeed, we can see from Fig. 20a that fuel concentrates in the cold regions inside the CRZ (or that the regions where fuel concentrate are colder) and is only completely consumed in regions with higher temperatures (roughly over 1200 K , which is an approximate threshold for the initiation of the fuel pyrolysis). This can be better seen on the mixture fraction space, analyzing fuel mass fraction over all the combustion chamber domain (Fig. 21). The analysis in mixture fraction space shows that it is in the high-temperature zone near the flame (Fig. 21a) that the most of the pyrolysis is occurring, as expected, as we see the production of both C 2 H 4 (shown here as a tracer of the pyrolysis products) and CO (shown as a tracer of the high heat release rate reactions, as its oxidation is one of most exothermic ones). However, when looking at the regions with temperature under 1200 K (Fig. 21b), one can see that in the colder and fuel-rich region where very little reaction is found, the fuel also pyrolyses. This corresponds in the physical space to the CRZ, where the fuel pyrolyses in the absence of the flame. Indeed, as conditions are too cold and too rich to sustain a flame, the recirculated heat enables pyrolysis nevertheless. This capability of breaking the fuel into smaller species by pyrolysis in the CRZ seem to Fig. 19 Volume-averaged probability of encountered states on mixture fraction space of temperature for a OP1, b OP2 and c OP3. The values on the temperature axis were hidden due to the confidentiality of this data be an important mechanism to help maintain a flame in such harsh conditions. Indeed, the pyrolysis of the fuel in the CRZ provides an initial source of reactants for the high-release rate reactions that helps feeding the flame around the CRZ, which does not need to evaporate and pyrolyse the inject fuel on the spot. Therefore, it seems that this observed partial Fig. 20 Mid-plane cross-sectional cut of the combustion chamber: Instantaneous fields of Temperature (normalized, blue regions T < 1200 K , yellow to pink, T > 1200 K ), Mass Fraction of fuel, of C 2 H 4 and of CO . For the species fields, shades of blue represent values under the fuel stoichiometric proportions, and yellow to pink, values over the stoichiometric proportions spatial decoupling of pyrolisis-only in the CRZ (where conditions are too cold and too fuel-rich to burn, but pyrolysis is possible) and fully-reacting flame combustion in the ISL and dilution zone (where the products of the CRZ pyrolysis are essential to allow the stabilization of a flame) are two parts of the stabilizing mechanism for sustained combustion in high-altitude windmilling conditions for this combustor, highlighting another benefit of the swirling flow topology. Looking at the spatial distribution of C 2 H 4 and CO , illustrated in Fig. 20a, helps visualising this behavior. First, one can see how C 2 H 4 is mostly present in the regions of temperature inferior to 1200 K , but is almost absent in high temperature regions, as it is one of the source for the high-release rate reactions. Finally, the Fig. 21 OP1 volume-averaged probability of encountered states on mixture fraction space for fuel mass fraction (top row), C 2 H 4 mass fraction and CO mass fraction. a represents all CMC cells, while b considers only regions with T < 1200 K , to highlight the conditions at which pyrolysis is happening 1 3 CO presence mostly at the vicinity of the flame and, therefore, at the end of the CRZ and downstream of the fuel-rich zone, further highlights the spatial decoupling between the pyrolysis-only region and high heat-release rate reactions (fully burning regions).
Changing to OP2 conditions (Fig. 20b), representing an operating point a few seconds after the restart of a flame and of the compressor operation (therefore, with a considerable increase in air mass flow rate with a slight increase in temperature- Table 1), shows an improvement of this partial spacial decoupling of fuel consumption. In these conditions, the pyrolysis-only zone moves upstream and is restricted to the region immediately downstream of the injection. This is the only region of the CRZ to present low temperatures (inferior to 1200 K ), fuel-rich mixture, and very low HRR and CO concentration. This shows how the quality of the spray atomization and quantity of fuel being injected are central to the stabilization of the flame and its optimal performance. The plots in mixture fraction space (Figs. 19b and 22a) confirm this analysis, showing the reminiscence of "low" temperature pyrolysis (blue paths in the three graphs of the left column) Finally, looking at OP3 (Fig. 20c), one confirms what was previously discussed concerning the similarity in flame shapes between cases, but distinct flame structures. Indeed, one can see how despite the three cases having similar fields of mixture fraction, in each of them the mixture fraction represents different things inside the CRZ. For OP1, the CRZ is mostly filled with cold fuel vapor, while for OP3, the CRZ is filled by final combustion products, with OP2 somewhere in the middle. The same is seen on the mixture fraction space (Figs. 19c and 22b), where the shows an almost exclusive fully-burning behavior with the pyrolysis happening very intensely close to the flame front. It is also worth noting how the high temperature and low AFR conditions of OP3 allow for the intense decomposition of fuel in richer conditions when compared to the two previous cases, as in OP3 the CRZ now is uniformly filled with hot gases. These results highlight very well the impact of operating conditions and spray characteristics on the combustion performance, as they portray how the flame evolves from a non-ideal sub-idle situation to close-to-optimal conditions, much improving the overall system performance.

Conclusion
In this work, the spray combustion process during aero engine relight was studied at one of the most challenging conditions. The focus was on the prediction of the primary atomization statistics by means of SPH, accompanied by corresponding combustor LES-CMS simulations, which utilize the statistics as spray starting conditions. Three operating points during the pull-away process were investigated (OP1, OP2 and OP3), where for each case primary atomization was simulated with the SPH method, providing spray statistics subsequently used as spray starting conditions for the following reacting LES-CMC simulations. This study enabled a deep understanding of operating conditions on the spray atomization and the impact, of both spray and operating conditions, on the flame structure of an industrial RQL burner.
In order to enable computationally affordable SPH simulations of the primary atomization process, the two dimensional simulation framework of Dauch et al. (2017) was extended. It captures the most important characteristics of the inherently three dimensional primary atomization process, as we demonstrate in our work. The statistical spray analysis revealed that the fragment size distribution, as one of the most important spray metrics, is characterized by a bimodal shape reflecting two distinct fragmentation classes produced in the disintegration process. This could be vividly explained using a FTLE analysis to study the aerodynamic interaction with the kerosene lamella. One of these fragmentation classes is represented by small, almost spherical droplets with low inertia, whereas the other class comprises a small number of large fragments, which are of statistical importance due to their volume. This behaviour was also evident in a three dimensional, circumferentially periodic simulation of OP1 with a 20 • sector model. It must be stressed, however, that the bimodality could also be a consequence of lacking convergence for the large and rare liquid fragments, even for the two dimensional simulations. A significant extension of the simulated time ranges in future work could resolve this issue. Moreover, the three dimensional simulation reveals that the transition of the kerosene lamella to the spray, resulting from the primary atomization, departs significantly from the picture developed from the two dimensional FTLE analysis. This is rooted in the circumferential instabilities developing in three dimensions, which are unaccounted for in the two dimensional setup. Another important insight is that the two fragment classes are characterized by different dynamics, as reflected by the bivariate statistics. The small droplets follow the air flow and its fluctuations well, leading to significant velocity dispersion in the corresponding scatter plots, whereas the large and rare fragments are characterized by high inertia, and hence are likely to follow the movement of the former kerosene lamella. This behaviour is related to a Stokes number effect, eventually leading to characteristic scatter plots, which is also observable in experiments Soni and Kolhe 2021;Vallon et al. 2021). We understand the reproduction of this characteristic bivariate statistics as one of the biggest strengths of our two dimensional modelling approach.
From an engineering point of view, the statistical insights gained, and especially the bimodal fragments size distribution with its statistical provisos, strongly indicate that the considered airblast atomizer is operated far from its design point during pull-away. Taking asymptotic theories as reference, one would expect to observe a monomodal shape of the fragment size distribution (Cheng and Redner 1990;Gorokhovski and Saveliev 2003;Villermaux 2007;Gorokhovski and Saveliev 2008;Villermaux 2020). This is an interesting observation, especially seeing how the combustion process at these conditions are affected by the spray. The Euler-Lagrange LES-CMC simulations in Sect. 3.3 show that for these extreme operating conditions spray initialization do matter. We observe that directly after the first relight during pull-away the very large fragments fill the CRZ with liquid and fuel vapor. These adverse conditions for the flame prevent any burning inside the CRZ, which becomes a pyrolysis-only region. This partial spatial decoupling of combustion and pyrolysis is, however, useful to help maintain the flame in such hard conditions, feeding the flame with smaller, easier-to-burn species. It is also important to mention that combustion efficiency is very low and well predicted by the LES-CMC simulations, with only a 4% error when compared to Rolls-Royce rig data. Consequently, it is likely that the engine at these conditions will be sensitive to external perturbations, which is also in accordance with the engine manufacturer's experience. However, with increasing engine load the spray quality improves significantly as idle operation is approached. As a consequence, the pyrolysis reactions increasingly subordinate to the usual combustion reactions in the CRZ, with the flame reaching the usual fully burning structure at OP3, finally proving that the standard RQL combustor mode can be restored. The combustion efficiency for both OP2 and OP3 are also well predicted, showing the effectiveness of the detailed coupled numerical strategy employed in this work. This methodology presents itself then as an interesting tool when exceptional circumstances prevent detailed experimental measurements of flow boundary conditions and primary atomization, allowing for good predictions and insights on these challenging conditions for the engine.