Predicting the Ignition Sequences in a Separated Stratified Swirling Spray Flame with Stochastic Flame Particle Tracking

Stochastic flame particle tracking in conjunction with non-reacting combustor simulations can offer insights into the ignition processes and facilitate the combustor optimization. In this study, this approach is employed to simulate the ignition sequences in a separated dual-swirl spray flame, in which the newly proposed pairwise mixing-reaction model is used to account for the mass and energy transfer between the flame particle and the surrounding shell layer. Based on the flame particle temperature, the particle state can be classified in to burnt, hot gas, and extinguished. The additional state of hot gas is introduced to allow the flame particles with high temperature to survive from nonflammable region and then potentially to ignite the nearby favorable regions. The simulations of the separated stratified swirl spray flame reveal two different ignition pathways for flame stabilization. The first showed that some flame particles from the spark would directly enter the main recirculation zone resulting from the velocity randomness and then ignite both sides of the combustor simultaneously. The second showed that flame particles from the spark would ignite the traversed regions following the swirl motion inside the combustor. The predicted ignition sequences were compared with the evolution of flame morphology recorded by high-speed imaging from experiments, showing qualitative agreement.


INTRODUCTION
Low pollutant emissions and high combustion efficiency can be achieved by the lean-burn concepts such as Lean Premixed Prevaporised (LPP) and Lean Direct Injection (LDI) . The pollutant emissions and flame stabilization can be balanced by the centrally staged LPP technology, as demonstrated by the Twin Annular Premixing Swirler (TAPS) (Foust et al., 2012) and the TeLESS-II (Technology of the Low Emission Stirred Swirl) combustor (Wang et al., 2016. For these technologies, more air is put into the main stage, making it harder to ignite due to stronger convection effect that inhibits the transportation of fuel droplets to the igniter tip. Furhter investigation of the forced ignition process in the centrally LPP combustors are needed for practical importance. As observed in recent experiments (Yang et al., 2020), the forced ignition of the swirl spray flame can be divided into five phases: (1) kernel generation, (2) flame propagation, (3) flame residence, (4) flame growth, (5) stable combustion. The flame residence has been reported in some realistic combustors (Mosbach et al., 2010, Read et al. 2010, during which the flame kernel remains small or even invisible before the kernel grows obviously after a certain time This phenomenon has also been observed in bluff-body flames (Ahmed et al., 2007) and was named as "ignition delay" (Naegeli and Dodge, 1991). This residence phase is not included in the usual classification of ignition process (Eyssartier et al., 2013). Although studies have observed the flame residence phenomenon, there is lack of understanding of the controlling physiochemical mechanism of flame kernel evolution during this phase.
The flame evolutions in the early stages of the forced ignition process in nonpremixed flames have been suggested to be mostly influenced by convection and turbulent diffusion (Neophytou et al., 2010). And based on the observations from experiments , Heeger et al., 2009) and simulations (Chakraborty et al., 2006, Neophytou et al., 2010 for both gaseous and liquid fuels, the assumption that the displacement speed of flame in non-premixed flows is much smaller than the turbulent fluctuation speed (Neophytou et al., 2012) is made to develop a low-order ignition model to predict the ignition probability in different nonpremixed flames. For such model, only cold flow simulations are needed to provide the background environments, and the stochastic Lagrangian particles are used to mimic the motion of flame elements, and the local Karlovitz number is computed and used as the criterion to check if the flame particles can ignite the regions they pass through. This method has been extended to predict the lean ignition limit (Soworka et al., 2014) and to simulate the light-round process (Machover and Mastorakos, 2017). Similar efforts have been made to use local quenching criteria combined with Fokker-Planck equation to simulate the stochastic evolution of flame kernels (Esclapez et al., 2020). And this method has been tested in premixed, non-premixed, and spray flames and have shown good agreement of ignition probability with experiments.
The reliability of the distribution of ignition probability predicted by these low-order ignition models depends on the correct capture of the ignition sequences, i.e., the pathways that the flame kernels evolve along. However, the assumptions made by these low-order models do not recognize the importance of the chemical kinetics, which may cause large error of the predicted time of different phases in the whole ignition processes. And it is highly suspected that chemical kinetics may play an important role for capturing the "ignition delay" or the residence phase in the centrally staged LPP combustors mentioned above. Recent work (Xie et al., 2021) has proposed the pairwise mixing-reaction (PMR) model to include the influence of chemical kinetics in describing the evolutions of particle states. And the particle state is decided by its own temperature instead of the Karlovtiz number. This new method has been tested in the nonpremixed methane/air bluff-body flame and the predicted ignition sequences and ignition probability agree well with experiments. But it still needs further validations in the swirl spray flames.
So the objective of this work is to investigate whether the low-order ignition model can capture the residence phase of the ignition sequences in the centrally staged LPP combustor, TeLESS-II, with the help of the PMR model. This paper is organized as follows. First, the stochastic Lagrangian flame particle method is briefly overviewed. Next, the PMR model is elaborated with details. Then the flame configuration is briefly introduced and the cold flow simulation is compared with experiment. Finally, the predicted representative ignition sequences are discussed and compared with experiments, and potential further development of the low-order methods are provided in the conclusion.

Igniter
Cooling Air Dome Cooling Main Main Pilot number is 0.5. The swirl numbers of the two pilot stage radial swirlers are 0.8 and 0.6 respectively. The igniter is located 50mm downstream of the pilot swirler exit. The test fuel is Chinese aviation kerosene RP-3 (Zhang et al., 2019). Only the pilot stage fuel was supplied in the ignition test. Detailed characteristics of the ignition system can be found in previous work . The distribution of spray was obtained by a Planar Mie Scattering system. The ignition process is recorded by a high-speed camera. Figure 1b illustrates the evolution of flame area from one representative successful ignition test. As shown, the ignition process in the separated stratified spray flame can be divided into four phases: the flame propagation phase S1 from the moment 0ms to the inflection point P1; the flame residence phase S2 from P1 to the inflection point P2; the flame growth phase S3 from P2 to the inflection point P3; and the stable combustion phase S4 from P3 to the end. In phase S1, a high temperature region is generated after the spark, and it heats up the local flow to form a strong kernel. Then the hot kernel propagates toward the centreline of the combustor, and its size shrinks rapidly. Then the flame enters the residence phase S2, and the flame size remains at a small scale. After some time, the flame enter phase S3 and begins to grow large rapidly. Then the flame reaches the phase S4, where stable combustion is achieved and the flame size may change periodically.

Lagrangian Flame Particle Tracking
As pointed out (Neophytou et al., 2012), a flame element can be represented by a flame particle that travels following the flow motion if the flame motion relative to the flow is slow compared to convection and turbulent diffusion. The cold flow simulation and fuel/air mixing fields provide the flame particles with the necessary background information to evolve. The evolutions of flame particle location and velocities are governed by the Langevin model for turbulent flows (Haworth and Pope, 1986) , , , . , where , is the flame particle location, , is the particle velocity in direction , is the mean velocity of the cold flow at the particle position, , is a random variable of normal distribution with zero mean and unity variance, is turbulence frequency calculated by turbulent kinetic energy and dissipate rate at the particle position: / , and is the model constant usually set to be 2.0.

Pairwise Mixing-Reaction Model
The flame particle is assumed to be a gaseous mixture being surrounded by a shell layer that consists of the cold flow mixture and possible fuel droplets in the local CFD cell. Mass and energy will transfer between the flame particle and the shcell layer. The liquid fuel droplets are allowed to evaporate when the shell layer is heated up by the hot flame particles. The composition of flame particle is described by the species mass fraction and the sensible enthalpy ℎ of the gas mixture, , , ⋯ , , ℎ , where is the number of species. The gaseous composition of the shell layer consists of the same state variables. The structure of the Pairwise Mixing-Reaction (PMR) model is depicted by the threelayer of the flame particle, surrounding shell and the cold flow. The chemical reaction, turbulent mixing and evaporation of the droplets are thus modelled to account for the mass and energy interaction between the flame particle and its surrounding shell. The states of local cold flow is used to initialize the shell layer, while mass and energy transfer between the shell and local cold flow in the cell is not considered. Over a characteristic coupling time , the interactions between the flame particle and the surrounding shell mixture are described by where and are the chemical reaction rates for the flame particle and the gaseous shell mixture , respectively, is the scalar mixing timescale and it is modelled as / with being the mixing constant. The turbulent timescale is defined as / with and being the turbulent dissipation rate and kinetic energy of the gaseous mixture, respectively. The evaporation source is only non-zero for the single fuel species. The evaporation models developed by Abramzon and Sirignano (Abramzon and Sirignano, 1989) are used to evaluated by the Lagrangian droplets information together with the mean gas phase properties interpolated at the flame particle position. The droplet information can be extracted from the Discrete Phase Model (DPM) in the non-reacting flow simulation.
Described by the above pairwise mixing-reaction model, during its evolution, a flame particle could have three different states of burnt, hot gas and extinguishe according to the particle temperature. The extinguished particles are no longer tracked, while only the burnt particles are allowed to ignite more cells.

Predicting the Ignition Sequence
At the beginning of the ignition simulation, the flame particles are uniformly distributed in the specified ignition region. The particle mixture fractions are generated using the Beta Probability Density Function (Beta-PDF) given the mean and variance of the mixture fraction that are extracted from the local cold flow. The particle velocities are initialized using the Gaussian Probability Density Function (Gaussian-PDF) given the mean and fluctuation velocities that are also extracted from the local cold flow. Once released, the flame particles are tracked through the following steps: Step.1 The cells are allowed to have only two states: unburnt and burnt. All cells are initialized with unburnt state. Each cell that overlaps the specified ignition region will release a "flame particle". The initial of flame particles are obtained through chemical equilibrium calculations using the generated random mixture fractions and user-specified ignition energy, . The flame particle is labeled as burnt if , or hot gas if , where is particle temperature, is the ignition temperature (the minimum temperature that allow local air and fuel mixture to ignite), and is the extinguish temperature. The Jet-A mechanism with 16 species and 23 reactions (Kundu et al., 1998) is used to simulate the combustion processes of the flame particles and their shells and the in situ adaptive tabulation (ISAT) (Lu and Pope, 2009) is used to accelerate the chemistry integration.
Step.2 The Langevin model Eq. (1) and Eq. (2) is used to simulate the turbulent motion of flame particles. The flame particle composition is updated via Eq. (3) and Eq. (4). With the interactions between flame particles and the shells, the flame particle could evolve from burnt state to hot gas and may eventually extinguish and no longer be tracked.
Step.3 Once a flame particle visits an unburnt cell, and the particle temperature stays beyond , and the local cold flow mixture fraction is within the flammability limits, the cell switches to the burnt state. Then, a new particle will be released at the cell center, and its initial states are generated by the Beta-PDF and Gaussian-PDF mentioned earlier. And each ignited cell are allowed to release only one new particle.
Step.4 For each simulated ignition process, the burnt and unburnt status of cells can be counted. The Ignition Progress Factor (IPF), defined as the fraction of cells in burnt state for a prescribed domain of interest, can be obtained to indicate the possibility of successful ignition.
Step.5 For each spark location, 60 repetitions of flame particle tracking with distinct realizations are performed to obtain the probability of the successful ignitions.

The Cold Flow Characteristics of the TeLESS-II Combustor
The cold flow Reynolds-Averaged Navier-Stokes (RANS) simulation of the TeLESS-II combustor is performed using the CFD software ANSYS Fluent. The standard -turbulent model is used to properly capture the turbulent swirl and recirculation flow. And the partially premixed combustion model is used to obtain the gas phase mixture fraction and variance fields, deactivating the reactions. The liquid phase is simulated using Discrete Phase Model (DPM). The Rosin-Rammler drop-size distribution is used with an averaged droplet diameter of 35μm and the spread parameter of 3.0. The half spray cone angle is 40deg. The discrete random walk model is used for stochastic tracking. The temperature of injected fuel is 300K and the fuel vaporization temperature is 288K. The simulated total fuel to air ratio (FAR) is 0.025 and the total pressure drop is 3% under the conditions of 101kPa and 300K.   Figure 2b. For the experiment, the averaged position of stagnation point is located around Y=50mm on the axis (Figure 2a). Yet for the RANS simulation, the stagnation point is located much further downstream around Y=78mm on the axis (Figure 2b). The overall characteristics of the non-reacting flow is well captured by RANS simulation. Figure 3 shows the spatial distribution of fuel spray droplets. From PMie, a large layer of fuel droplets with relatively low intensity extend downstream along the inner swirl jet. And the boundary of droplet mists reaches stagnation point on the axis at around Y=50mm (Figure 3a). Yet for the RANS simulation, the DPM particles majorly concentrate along the swirl jets and the distribution of droplets on the axis is largely underestimated (Figure 3b). This deviation of the predicted distribution of fuel droplets on the axis could cause significant difference in predicting the ignition sequences, which would be further discussed in the next section.   The computational domain for stochastic flame particle tracking is sketched in Figure 4a. The spatial ranges along X, Y and Z axises are (-55,55)mm, (13,120)mm, and (-70,70)mm, respectively. Uniform structured mesh is generated for particle tracking with the spatial resolution of 1mm. The simulation time is 15ms with the step size of 0.04ms. The spark kernel is located at (0,50,50)mm with the initial radius of 1cm. According to the experiment (Yang et al., 2020), the ignition energy is set to be 2J unless otherwise specified, which will make the initial particles of the ignition kernel reach equilibrium states. The ignition sequences being displayed by the particles' positions and temperature will be shown in the front, side and top views.
The ignition sequences are simulated with two sets of model parameters. The first set has a high ignition energy of 2J and a low mixing constant of 2.0. The second set has a lower ignition energy of 1J and a larger mixing constant of 8.0, i.e., faster turbulent mixing. With the two different sets of parameters, the simulated evolutions of the number of active (burnt and hot) particles, , from 60 stochastic ignition processes, are displayed in Figure 4. As shown, will eventually decrease to zero, since each cell is allowed to release only one stochastic flame particle when it is ignited. For the first set of parameters (Figure 4b), shows two different forms: one has only one peak, the other has two peaks, indicating two different ignition pathways. For the second set of parameters (Figure 4c), shows majorly two peaks. And with faster turbulent mixing, is overall much smaller than that from the first sets of parameters, since a number of burnt particles will be more likely to cool down rapidly and extinguish before they ignite the next favourable cell. The value of indicates how much area of the combustor is ignited. And for the first set of parameters, the increases immediately after the ignition energy release. While for the second set of parameters, the will first drop to a very low value and then increase, which agrees qualitatively with the experimental observation of the flame propagation and residence phases as shown in Figure 1b. However, the simulated time duration of these two phases and the whole time of the ignition process are overall much shorter than the experimental observations (Yang et al., 2020). Two typical ignition sequences being displayed by the particles' positions and temperature are shown in Figure 5 and Figure 6, respectively. The first is from the 48th spark event being marked by dash-dot line as shown in Figure 4b. The second is from the 59th spark event as shown in Figure 4c. As can be seen from the side views the first ignition sequence, the ignition kernel first moves to the left of the combustor following the swirl motion (Figure 5b, 1.8ms). Then, some of the flame particles move directly across the middle plane thus igniting both sides of the combustor simultaneously ( Figure  5b, 2.8ms). The particles on the right move against the swirl and are much slower than those from the left (Figure 5b (Figure 5c, 5.6ms). As mentioned above, the of the first ignition sequence shows only one peak. This is because, under the influence of stochastic turbulent fluctuation of velocity field, some of the ignited particles move across the middle plane against the swirl in the early stage of ignition process, and both sides of the combustor are ignited simultaneously. As can be seen from the second ignition sequence, the ignition kernel first shrinks and moves along the swirl (Figure 6, 0.6ms). And the number of ignited particles remains very low (Figure 4c) until the particles meet with suitable region to ignite more cells (Figure 6, 1.4ms). Due to more rapid turbulent diffusion, the temperature of ignited particles drops faster, so that they can only ignite much smaller region of the combustor (Figure 6, 4.6~9.2ms) compared with the first ignition sequence shown in Figure 5. The of the second ignition sequence shows two peaks. And the first peak corresponds to the moment that the bottom left region is ignited (Figure 6b, 4.6ms). The second peak corresponds to the moment that the upper right region is ignited (Figure 6b, 9.2ms). And the second ignition sequence totally follows the swirl motion and ignite the combustor counterclockwise as seen from the side views (Figure 6b). In the two typical simulated ignition sequences shown above, the ignited regions are displayed by the stochastic flame particles being coloured by high temperature. However, the residence phase, when a small hot region in the middle of the combustor lasts several milliseconds observed in experiments (Yang et al., 2020), cannot be seen in the simulations. One possible reason is that the simulated non-reacting flow underestimated the concentration of fuel droplets in the middle of the combustor as mentioned in the last section, so that the simulated particles would rapidly extinguish when they reach the axis. Another possible reason is that the stagnation point of the RANS simulation is located much further downstream with a greater distance from the fuel droplets, so that the stochastic particles would stay much shorter time in the middle of the combustor. Furthermore, for the consideration of computational efficiency, each ignited cell can only be allowed to release one particle, so the stochastic particles cannot be generated from the same ignited region continuously.

CONCLUSIONS
The ignition sequences in a separated stratified swirling spray combustor are simulated with the stochastic flame particle tracking method and the pairwise mixing-reaction model. Two different ignition sequences are identified. For both ignition sequences, the flame particles will generally follow the swirl motion. For the first ignition sequences, due to velocity randomness, some of the particles will move across the middle plane in the early stage of ignition process and ignite both sides the combustor. And the number of active particles shows only one peak. For the second ignition sequences, the particles strictly follow the swirl and ignite different sides of the combustor in sequence. And the number of active particles show two peaks. And with stronger turbulent diffusion effect, the number of active particles in the second ignition sequences will first drop to a low value and then increase, which shows qualitative similarity with the flame propagation and residence phases observed in experiments. The flame particle tracking with pairwise mixing-reaction model shows great potential in analysing the possible ignition pathways for realistic aero engine combustors. Yet a lot more validations are still needed for the turbulent spray combustors. Further attention should be paid in the future studies, including improving the accuracy of non-reacting flow simulation, analysing the detailed evolution of individual particle states, and testing different particle release strategies.