Influence of Thickening Factor Treatment on Predictions of Spray Flame Properties Using the ATF Model and Tabulated Chemistry

Different strategies to account for the heat and mass transfer between liquid droplets and their carrier phase within the Artificially Thickened Flame (ATF) approach are analyzed and compared. Herein, two approaches are introduced to take into account the droplet movement relative to the thickened flame front orientation. While the first approach achieves this behavior through scalar modifications in the droplet temperature and mass evolution equations, the second one introduces a trajectory modification within the thickened flame front. Both approaches, referred to as projection and refraction correction, are first compared to state of the art methods in a simplified two-dimensional configuration, and then in a complex turbulent spray flame. The investigated spray flame corresponds to the operating condition EtF6 of the Sydney Spray Burner. The analysis showed that: (1) A consideration of a simplified configuration is insufficient to fully uncover the performance of the different approaches. (2) While the proposed approaches performed best in the two-dimensional configuration, only the projection method outperforms the remaining ones in the turbulent spray flame. (3) The formulation to consider the flame thickening has a strong effect on global flame properties, combustion regime distribution as well as carrier and liquid phase statistics.


Introduction
One outstanding advantage of liquid fuels when compared with their gaseous counterpart is their high energy density at ambient conditions, making them generally the number one choice for transport and storage issues. Therefore, a lot of attention has been given to combustion devices fueled with liquids in the last decades. However, applications involving such fuels include processes like atomization, evaporation, and mixture formation coupled in a strong manner within a turbulent reacting flow. The complexities associated with these strongly coupled phenomena turn the characterization of the turbulent spray combustion to a quite challenging task. In view of the harmful emissions and pollutants resulting from the combustion of fossil fuels, and considering the increasingly restrictive environmental regulations, there is a strong need for more understanding of the turbulent spray combustion.
With respect to numerical simulations, Direct Numerical Simulation (DNS) could aid researchers and engineers to better characterize and understand complex flow processes. However, DNS generally involves prohibitive computational costs, which makes the development of computationally affordable models highly required. In this regard, Large Eddy Simulation (LES) has proven its capability to accurately handle complex unsteady flow phenomena (Janicka and Sadiki 2005;Pitsch 2006).
Regarding turbulent combustion simulations, it is essential to connect the realistic fuel chemistry needed and the computational costs related to the LES. For that purpose, various chemical mechanism reductions (Marinov 1999;Smooke 1991) and chemistry tabulation or storage/retrieval approaches have been suggested and widely utilized, for instance in Kuenne et al. (2012, Chrigui et al. (2013), Dressler et al. (2019). In this respect, a wide range of approaches, including the intrinsic low dimensional manifold (ILDM) (Maas and Pope 1992), the reaction-diffusion manifold (REDIM) (Maas 2007), the flamelet progress variable approach (FPV) (Pierce and Moin 2004), the flame prolongation of ILDM (FPI) (Gicquel et al. 2000;Fiorina et al. 2005) or the flamelet generated manifold (FGM) (Van Oijen and De Goey 2000), among others, have been proposed. These approaches make use of a so-called mapping of the thermochemical state onto a reduced set of n lookup scalars, usually reducing the number of transport equations to be solved to the Navier-Stokes-plus n additional scalar equations. In theory, the n-dimensional table can store an arbitrary amount of information, which is often in practice, limited by the available memory at runtime.
Another issue in the context of LES arises when considering the small flame thickness compared to the LES grid spacing, making the closure of the equations needed to model turbulent combustion a subgrid scale matter. Thereby, the accurate modeling of the turbulence-flame interaction plays a central role in model development. The fundamental approaches to assess this problem can generally be divided into three categories (see Poinsot and Veynante 2005), namely geometrical methods (Pitsch and de Duchamp Lageneste 2002;Hawkes and Cant 2000), statistical approaches, as the presumed (Ge and Gutheil 2008) or transported (Jones and Kakhi 1998;Hu et al. 2017) PDF method (Pope 1985), and algebraic approaches, as the thickened flame or artificially thickened flame (ATF) model. The latter, first applied to single phase, laminar combustion (Butler and O'Rourke 1977) has been extended to turbulent combustion (Charlette et al. 2002;Colin et al. 2000), to stratified combustion using a grid adaptive dynamic thickening strategy (Kuenne et al. 2012) and later to multiphase combustion Cheneau et al. 2019).

3
The present work aims at applying the ATF-model-based Large Eddy Simulation coupled to tabulated chemistry to investigate turbulent spray combustion under consideration of issues related to the thickening factor. Herein, simulations are performed for the Sydney Diluted Spray Burner configuration from the University of Sydney (Gounder et al. 2012;Gounder 2011, 2010). It is worth mentioning that this configuration constitutes a test case within the Workshop on Turbulent Combustion of Sprays (TCS, www.tcs-works hop.org).
The Sydney Diluted Spray Burner has been numerically investigated by many different groups (Chrigui et al. 2013;Heye et al. 2013Heye et al. , 2014Sacomano Filho et al. 2018;Rittler et al. 2015;De et al. 2011) as it provides an extensive experimental database for model validation as well as a reduced modeling effort due to the dilute nature of the spray. (Chrigui et al. 2013) performed LES of the ethanol Flames EtF3 and EtF8 using a Lagrangian approach coupled to tabulated chemistry and a presumed PDF based turbulence combustion interaction model. (Rittler et al. 2015) validated an approach combining a presumed PDF method, artificial thickening and tabulated chemistry for the cases EtF3, EtF6 and EtF8 and evaluated the influence of the chosen sub-filter distributions. (Heye et al. 2013) investigated the spray flame EtF1 numerically using LES with a PDF based Lagrangian Monte-Carlo approach for the flame-turbulence interaction. A valuable contribution is another work by Heye et al. (2014), where the authors compare simulations performed by four different research groups for the flames EtF2 and EtF6. The simulations feature -besides diverse numerical treatments and codes -different spray injection strategies, turbulence modeling, evaporation modeling, chemistry treatment and combustion modeling, including FGM with presumed and transported PDF methods. The authors point out that the determining factors are: (1) The inflow conditions chosen for the spray, (2) the evaporation modeling as well as (3)  In regard of the application of the ATF model in LES of turbulent spray combustion, the contributions are scarce. To mention are the work of Cheneau (2019) and Cheneau et al. (2019), where this method is applied in a swirled combustor. These works are grounded on the preparative work of Boileau (2007) and the validation by Paulhiac (2015). Also of importance is the contribution of Sacomano , where a lean partially prevaporized spray flame is numerically investigated. Further applications of this approach can be found in Rittler et al. (2015), Sacomano Filho (2017) (2020). However, all these contributions generally differ in their strategy to correct the interaction of the thickened flame with droplets.
The objectives of this work are on one hand to spread awareness towards the complexities that arise by combining flame thickening approaches with an additional disperse phase. In this regard, this paper proposes and evaluates two alternative formulations for the correction of droplet heat and mass transfer which takes into account alignment of droplet movement and flame front orientation. The methodologies are first investigated and compared with state of the art methods in a simplified two-dimensional setup. The approaches are then applied to predict the combustion properties of the dilute spray combustion configuration EtF6 of the Sydney Spray Burner (Masri and Gounder 2010;Gounder 2009). Thereby, the differences arising from the various particle-thickened flame interaction treatments are thoroughly discussed.
The remaining of this manuscript is structured as follows. First, the modeling approaches (Sect. 2), the investigated configuration as well as the numerical setup (Sect. 3) are presented. In Sect. 4, the applied methods are then compared in configurations with increasing complexities. Thereby, we start from a simplified configuration before going over to the complex turbulent spray flame EtF6. The achievements are finally discussed and summarized in Sect. 5.

Numerical Methods
All investigations are performed using the finite volume open-source C++ code Open-FOAM (Weller et al. 1998) including a low-Mach modification as described in Ries et al. (2017). Dealing with reacting turbulent two-phase flows, the carrier phase motion is represented by the transport equation for mass and momentum: Herein, ( ⋅) corresponds to Favre-filtered and ( ⋅) to filtered quantities. denotes the density, u i the Cartesian velocity component, the dynamic viscosity, g i the component of gravitational acceleration, p the hydrodynamic pressure. The terms S m,l and S u i ,l represent the liquid phase contribution source terms resulting from mass and momentum exchange with the Lagrangian parcels (two-way coupling). The unresolved turbulent fluxes, which are accounted for through the subgrid scale stress tensor sgs ij are modeled using the Sigma eddy-viscosity model (Toda et al. 2010). In this work, combustion modeling is realized by coupling a chemistry tabulation approach based on the FGM method (Van Oijen and De Goey 2000) to the Artificially Thickened Flame model (Butler and O'Rourke 1977). The later is employed in order to resolve the flame on the discretized domain, thus keeping the computational effort moderate. In regard of chemistry tabulation, a two-dimensional table is used which was generated from freely propagating premixed flamelets assuming a Lewis number of unity. The chemical mechanism proposed by Marinov (1999) was employed. The thermochemical state is tabulated as function of the mixture fraction Z (as formulated by Bilger (1989)) and a reaction progress variable (as adopted in Sacomano Filho (2017)): ( where Y i and W i denote the species mass fraction and molecular weight, respectively. In view of this, the previously introduced set of equations is extended by additional transport equations for mixture fraction Z and reaction progress variable PV: Thereby, ̇P V denotes the reaction progress variable source term and the terms S Z,l and S PV,l the respective source terms resulting from the presence of the liquid phase. Both equations follow the definition of a thickened scalar in the ATF context (Kuenne et al. 2012) as extended to multiphase reacting flows in Sacomano Filho (2017) where f stands for the laminar flame thickness at a specific mixture fraction, V 1 3 i the representative cell dimension and n CV the minimum number of control volumes to reproduce the laminar flame speed s f at a given mixture fraction (here, n CV = 4 ). denotes the flame sensor, which is zero in regions where no reaction occurs and evaluates to 1 near the flame. In this work, the formulation by Dressler et al. (2019) is employed for the flame sensor (see Eq. 8). The efficiency function E in Eqs. 4 and 5 corrects the distorted flame dynamics caused by the thickening of the flame. In this regards, the non-dynamic formulation by Charlette et al. (2002) is adopted with an exponential factor = 0.5 . A turbulent Schmidt number Sc t = Sc = 0.7 is assumed for both transported scalars.
The liquid phase is described in a Lagrangian manner, where the evolution of computational parcels, each representing a multitude of real droplets having the same properties, is determined by the following set of coupled ordinary differential equations for position, velocity, mass and temperature:

3
Here, ∑ i denotes the forces acting on the parcels (only drag (Schlichting 1960;Yuen and Chen 1976) and gravity are considered here), Sc, Sh, Pr and Nu the Schmidt, Sherwood, Prandtl and Nusselt number, respectively. Nusselt and Sherwood number are calculated according to Ranz and Marshall (1952). p = p d 2 p ∕18 stands for the particle relaxation time. The equations for parcel mass and temperature follow the form introduced by Miller et al. (1998). Herein, the parameter f 2 corresponds to a correction factor of heat transfer due to evaporation and H M to the mass transfer driving potential. In this work, the infinite conductivity model M7 from Miller et al. (1998) is employed which has been widely used in spray combustion studies and delivered decent results (Noh et al. 2018;. Within the adopted ATF-based modeling, an additional parameter F eff is especially introduced in Eqs. 11 and 12, corresponding to an effective thickening factor. The different approaches to calculate this parameter are presented later in this section. In order to derive the introduced equations for droplet heat and mass transfer, constant thermodynamic properties across the droplet boundary layer are assumed (Bellan and Summerfield 1978) and it is therefore crucial to estimate these quantities as accurately as possible. In view of this, the well-established "1/3"-rule ( r = s − ( s − ∞ )∕3 , = (T, Y i ) ) (Knudsen et al. 2015;Miller et al. 1998;Abramzon and Sirignano 1989) proposed by Yuen and Chen (1976) is applied here. Once the composition and temperature at the "reference" (subscript r) are known, the Wilke rule (Wilke 1950) is used to obtain the dynamic viscosity and thermal conductivity at reference conditions while the heat capacity is mass averaged. In the context of tabulated chemistry, the procedure can be summarized as follows: First, the composition and temperature are retrieved from the thermochemical lookup table. Secondly, reference conditions are determined according to Yuen and Chen (1976). Subsequently, all properties of interest are evaluated through NASA polynomials (McBride et al. 1993). Finally, the mixing rules are applied and mass and heat exchange computed.
The validation of the evaporation model is performed by comparing the diameter evolution predicted by the employed model with experimental data from Saharin et al. (2012). The experimental setup consists of ethanol droplets evaporating at different ambient temperatures in a nitrogen atmosphere. A comparison of simulated and experimentally measured data is presented in Fig. 1. It can be observed that the model reproduces the trends correctly and matches the experimentally determined droplet surface evolution quite well. The discrepancies between simulations and experiments which can be observed at lower diameters are caused by the presence of a second species in the liquid, as discussed in detail in Saharin et al. (2012) and Saharin (2013).
Differently from previous works Sacomano Filho 2017;, where the carrier gas is simplified as mixture of fuel and air, we followed the strategy proposed in Sacomano . Hereby, the mixture surrounding droplets is assumed as compound of 8 species (see Table 1). In Sacomano Filho et al. (2019), it is shown that the combination of the 6 major species ( CO 2 , H 2 O , H 2 , CO , C 2 H 5 OH , and air (as a single species)), would be enough to represent the mixture.
In that sense, the present approach applies the findings of Sacomano Filho et al. (2019) obtained from laminar flames propagating in droplet mists to turbulent combustion. However, using 8 instead of 6 in the present work results in no significant increase in computational costs. This leads to 8 additional fields that have to be stored in the thermochemical lookup table. Please note that other authors did already perform detailed mixture description in a turbulent context (e.g. Franzelli et al. 2016;Ma et al. 2015).
The main concern of this work refers to the interaction of the liquid phase with the flame. This is relevant as the interaction of a droplet with a thickened flame will definitely differ from the one where no thickening is applied. Differences are expected to be especially high for heat and mass exchange. In literature, this impact is either neglected (e.g. Rittler et al. 2015) or the correction procedure, initially proposed by Sacomano Filho et al. (2017) is employed. Usually, the correction procedure consists in defining a correction factor for heat and mass transfer to the droplet in order to retrieve a correct flame exit diameter for droplets crossing the thickened flame. In the present work, the correction is expressed in terms of an effective factor F eff . Hence, the uncorrected and corrected approaches result in a different effective thickening factor F eff = 1 and F eff = F in Eqs. 11 and 12. The validity of the second method is shown in Sacomano , where flame front and droplet are moving orthogonally to each other. The definition is consistent with the initial idea of the thickened flame approach as both, movement of droplet relative to the flame front and the direction in which thickening is applied (namely flame orthogonal) are equivalent. However, this procedure is expected to introduce distorted heat and mass transfer for droplets moving with a substantial velocity component parallel to the flame front, as it is often the case in spray jet flames. This is outlined in Fig. 2, which shows a droplet traveling with a velocity p across a non thickened flame with thickness f or a thickened flame with thickness F f . The nonthickened and thickened flame front are represented in Fig. 2 by vertical black or thickened red lines, respectively. The path of a particle (blue line in Fig. 2) traveling with the velocity p across the flame front for an unthickened flame can be described by the states (a) (before the flame), (b) (flame front entry) and (c) (flame front exit). Taking advantage of the monotonic increase of the reaction progress variable, it is possible to assess the flame front orientation using the reaction progress variable gradient ∇PV , which is used to define the flame front orthogonal direction (i.e. ∇PV |∇PV| ) and the flame front parallel direction. While the length of the droplet path in flame front orthogonal direction is defined by the actual flame thickness, its parallel component is determined by its flame parallel velocity. The relation of both components is specified through the particle-flame interaction angle . As the droplet enters the flame, its temperature rises quickly, yielding strong evaporation rates which are distributed along the droplet movement path and thus influencing mixture evolution in both flame front normal and parallel direction. This is an important characteristic that must be kept in mind before going over to the analogous scenario with a thickened flame. In this case, the droplet follows an elongated path and the flame front exit position is shifted to (c � ) in Fig. 2. As stated in Butler and O'Rourke (1977), the idea of thickening is to expand the coordinate normal to a discontinuity in order to resolve it, which is a one-dimensional transformation. However, a particle interacting with a flame in a specific angle can no longer be perceived as such problem. If the flame is thickened in one direction, the path of the droplet within the flame is not only elongated in flame normal, but also in flame parallel direction, which will alter the mixture evolution in flame parallel direction. It is likely that such a constellation will influence flame dynamics.
In view of the issues related to the droplet displacement along the parallel direction of a thickened flame, we introduce two approaches which aim at reducing the stretching of the vapor release in flame front parallel direction. The first of these is obtained by limiting the correction of heat and mass transfer through the parameter F eff to the orthogonal direction of droplet displacement to the flame surface. The strategy is derived by considering the limiting conditions of a droplet moving (1) orthogonally or (2) parallel to the flame front. While for the first scenario, the original formulation found in Sacomano  or Cheneau (2019) should be recovered, the second scenario would yield no correction at all. Such properties are achieved by weighting the thickening correction with the projection of droplet movement onto the flame propagation direction. It is important to mention that this correction assumes a "frozen" or stationary flame (which is only true in a statistical way for the investigated turbulent spray flame). Accordingly, the correction factor F eff is defined as follows: As previously pointed out, the quantity ∇PV refers to the gradient of the progress variable, which is interpolated onto the parcel position and relates to the flame front orientation perceived by the parcel. The quantity p | p | denotes the normal unity vector defining the parcel movement direction. Various limiting conditions shall be considered from Eq. 13: (1) Far away from the flame, i.e. where no thickening is applied, the term (F − 1) is zero and no correction of heat and mass transfer is applied; (2) supposing ∇PV and p are orthogonal to each other, corresponding to a parcel moving parallel to the flame front (while thickening is applied, in other words F ≠ 1 ), the last term of the right hand side of Eq. 13 also vanishes, yielding F eff = 1 ; (3) in the case of a droplet moving orthogonal to the flame front while thickening is performed, the factor initially proposed by Sacomano Filho et al. (2017) is retrieved. Subsequently, this method will be referred to as projection correction.
A different manner to approach such a problem represents a second strategy, which will be referenced in the remainder of this work as refraction correction due to its geometrical resemblance with the phenomenon of light refraction. Here, the same formulation for F eff as used in Sacomano Filho et al. (2017), Cheneau et al. (2019), i.e. F eff = F , is employed. However, the trajectory of the droplet in flame parallel direction is limited to the unthickened flame front parallel displacement. Such a tracking correction yields the green path and the corrected exit position and state (c) ′′ in Fig. 2. In fact, such a correction is equivalent to reducing the particle flame interaction angle to ′′ . From a practical point of view, such path correction requires a passive transformation of particle movement in a flame front orthogonal coordinate system, of which the axes can be defined by evaluating ∇PV |∇PV| and its orthogonal vectors at the droplet position. Once the transformed velocity is obtained, the flame parallel components are scaled by the thickening factor F, while the flame orthogonal component is left unchanged. Subsequently, the corrected velocity is retransformed into the original reference coordinate system. This corrected velocity is then employed to track the particle along the thickened flame. It should be noticed that this velocity correction does not change the particle velocity p . This correction is only applied to update the droplet position during the droplet-flame interaction time. Thus, no change on momentum exchange or convective heat and mass transfer is explicitly applied here.
The difference between refraction and projection correction can be summarized as follows: The former modifies the distribution of vapor release in flame front parallel direction and at the same time keeps the droplet-flame interaction time equal to the formulation found in Sacomano . In contrary, the projection approach avoids any trajectory correction but changes the vapor release along the droplet path, which causes a (13) reduction of its lifetime compared with the refraction correction and the approach found in Sacomano Filho et al. (2017).

Configurations and Numerical Setup
The investigations presented in this work are performed for the operating condition EtF6 of the well known Sydney Spray Burner (Gounder et al. 2012). The burner consists of three concentrically arranged streams (see Fig. 3(a)). The spray is generated using an ultrasonic nebulizer placed almost 21 jet diameters upstream of the burner exit. The droplets are then transported downstream in the central pipe with a diameter D = 10.5 mm. The central jet is surrounded by an annular pilot consisting of burnt products resulting from a stoichiometric acetylene/hydrogen/air reaction. The burner is centered in a primary coflow with a diameter of 103 mm. The coflow is composed of pure air with an exit bulk velocity of 4.5 m s −1 .
The burner and coflow are placed in a wind tunnel with a cross section of 290 × 290 mm with same composition and velocity as the primary coflow. For the operating condition investigated, two experimental data sets exists, namely A and B (see Table 2). While both sets delivered similar results, some deviations are observed, thus providing an initial   Masri and Gounder (2011). The cylindrical computational domain consists of 4.3 million hexahedral control volumes as shown in Fig. 4. It ranges 75 jet diameters downstream in axial and 10.5 diameters in radial direction. Therefore, the secondary coflow, namely the wind tunnel coflow, is only partially included in the computational domain. A section of 6 jet diameters upstream of the burner exit is included for the central jet to generate valuable conditions at the burner exit plane. Here, a fully developed turbulent flow, attained with a recycling method, is assumed for the central jet which is in agreement with Gounder et al. (2012), Masri and Gounder (2011). Figure 3b presents an a priori analysis of the velocity profile at the burner exit plane for EtF6. The fully developed turbulent pipe flow assumed in our simulations is approximated by a power law profile ( n = 7 ). One can observe an underestimation of the velocity profile when using the bulk velocity specified in Table 2. A least square fit reveals that a higher bulk velocity is needed to match the experimental profiles. The bulk velocities were therefore adjusted accordingly (the value used in our simulations can be taken from Fig. 3b).
The pilot and wind tunnel inflows are considered laminar while random turbulence is applied for the primary coflow. In difference to Gounder et al. (2012) and similarly to Chrigui et al. (2013), De et al. (2011), the present work simplifies the pilot flame as stoichiometric ethanol-air flame. In order to keep the computational cost of the LES reasonable, the near wall regions are modeled through the near wall formulation by Spalding (1961). All further boundary conditions were set to total pressure boundaries allowing for entrainment of the surrounding flow. The disperse phase is injected at the axial position x = 0.3D corresponding to the first axial plane measured in the experiments. In this context, a size distribution conditioned on the injection position, obtained from experimental set A (Clean Combustion Research Group Database 2019), is applied. The disperse phase velocity at injection was simplified to the carrier velocity interpolated at the droplet injection position. It should be noticed that the employed injection strategy assumes an equal mass for all injected parcels. This results in a number of real droplets represented by computational parcels varying from parcel to parcel. The parcel mass at injection was estimated a priori in order to obtain a computationally affordable number of simultaneously tracked parcels. This number varied from around 2.2 to 2.9 million, depending on the modeling approach employed. By using this strategy, it is also possible to inject parcels representing less than 1 real droplet. This is the case for droplets with larger diameters. The impact of the numerical grid has been investigated and is presented in the Appendix, where a systematic grid sensitivity study is shown. Thereby, it is demonstrated that the grid dependency has no dominant effect for the employed grid.
The simulations were first initialized from a non-reacting flow and run for around 0.1 s, which corresponds to approximately 60 jet pipe flow through times, to allow the turbulence in the central jet to develop. Subsequently, the pilot was ignited and simulations were run for another 0.05 s before starting with the liquid injection. Thereafter, the simulations continued running for another 0.1 s before averaging was started. All simulations are averaged over a time period of at least 0.15 s. Additionally, to reduce averaging time, all carrier phase properties shown here are averaged in circumferential direction.
Before going over to the presentation and subsequent discussion of the results, a critical assessment of the chosen modeling strategy has to be performed. Some simplifications are made, some more important than others and these shall be discussed. First, since we are using a two-dimensional lookup table, we are only able to take the effects of chemical reaction and mixing into account. Thus, heat losses due to evaporative cooling of the liquid phase are not considered. Besides extending the thermochemical lookup table to a third dimension, a possible alternative could be to incorporate the effect of heat losses through linearization of the transported energy, as done by Franzelli et al. (2016). Another possibility is to consider the evaporative cooling effect as done by Knudsen et al. (2015), which is based on the assumption that all fuel in the fresh mixture originates from droplet evaporation (and based on a Lewis number of unity). By setting the enthalpy of the fresh mixture in the generation of the various flamelets accordingly, it is possible to construct a thermochemical lookup table that takes liquid evaporation into account. Lately, Sacomano Filho et al. (2019) demonstrated that the consideration of evaporative cooling delivers better agreements with experiments than when it is neglected. The impact was investigated for case EtF5 of the Sydney Spray Burner. However, it was shown that the effect was not much pronounced. Nevertheless, in order to focus on the investigated modeling aspects and to avoid the introduction of additional uncertainties, an adiabatic strategy as employed in previous works is used Sacomano Filho 2017;Sacomano Filho et al. 2017, 2020. At that point, it should also be noted that since the employed table consists of purely premixed flamelets, deviations can be expected in diffusion flame mode dominated regions (Franzelli et al. 2013;Vreman et al. 2008). For combustion of sprays, this was discussed in Hu and Kurose (2018) with acetone as fuel. It was shown that both regimes, namely premixed and diffusion flames, may coexist for various operating conditions of the investigated Sydney Spray Burner. Additionally, the employed ATF method is a strategy dedicated to the prediction of premixed flames. Its application to address diffusion flame is therefore not rigorous. However, for the EtF6 flame, experimental investigations presented by Masri and Gounder (2010) and Gounder (2009) suggest that the flame burns predominantly in a premixed mode. Therefore, the modeling inconsistencies in the presence of diffusion flames are expected to be small.
Another uncertainty is related to the employed efficiency function formulation. In this work the common definition by (Charlette et al. 2002) is adopted due to its broad applicability. However, an adjustment of the exponential factor may influence simulation results substantially, as demonstrated recently in the sensitivity study by Sacomano Filho et al. Additionally, no interaction of liquid droplets are considered here due to the dilute nature of the spray (i.e. four-way-coupling). Finally, to be consistent with Eqs. 4 and 5, a unity Lewis number was assumed in the computation of the iso-mixture fraction flamelets. This results in a slower flame propagation speed as forward diffusion of small species like Hydrogen is neglected.
In view of these simplifications, it should be noted that deviations with experimental data are expected. The comparison with experiments shall, however, guide the present analysis and assist the following interpretation and discussion of the results.

Results
This section is divided into several parts. The first of these is dedicated to the evaluation and comparison of the different approaches in a simplified two-dimensional configuration (Sect. 4.1). The different methods listed in Table 3 are: The approach without correction (i.e. F eff = 1 , e.g. Rittler et al. (2015)); the more common method, which is denoted as standard correction, where F eff = F , Sacomano Filho et al. (2017), Cheneau (2019); and the two procedures introduced here, namely the projection and refraction correction. Thereafter, Sect. 4.2 briefly discusses the influence of the flame sensor formulation used for the present spray flame by means of instantaneous contour plots. Finally, the LES results for the EtF6 spray flame are presented and discussed in Sect. 4.3. Especially, we clarify emerging differences caused by the various strategies introduced to treat droplet heat and mass transfer in the presence of a thickened flame. This is performed in terms of instantaneous contour plots as well as radial profiles of quantities of interest for both, the carrier and disperse liquid phase.

Approach Assessment in a Simplified Two-Dimensional Domain
The goal of this section is to demonstrate and clarify the impact of the various modeling approaches previously introduced. The test setup is rather simple and similar to the onedimensional test configuration presented in Sacomano Filho et al. (2017) and Cheneau (2019). The main difference is that the computational domain is extended to a second spatial dimension to allow for two-dimensional droplet trajectories. The main objective is to evaluate the flame diameter after a flame crossing when compared to a reference case and-no less important-to assess the droplet diameter evolution through a flame front parallel direction.

3
The rectangular computational domain employed is sketched in Fig. 5. It extends 0.02 m in x-direction and 0.01 m in y-direction and consists of 1000 x 500 control volumes, thus ensuring that the flame is well resolved in the reference case, i.e. where no artificial thickening is applied. The premixed flame, sketched as red line in Fig. 5, is considered stationary, which is also different from the validation cases employed in Sacomano , Cheneau (2019). For all investigated cases, a droplet, represented as blue circle in Fig. 5, travels through the domain with a constant velocity | p | = 0.3 m s −1 . Thereby, the angle describes the relation between flame orthogonal and flame parallel velocity component and can also be used to quantify the angle of interaction between droplet and flame in this case. As the droplet travels along its path, it crosses the flame which, depending on the approach used (see Table 3), produces different heat and mass transfer rates along the droplet trajectory. To make all approaches more comparable and following , heat and mass transfer is only activated in regions were thickening is performed, which is achieved reasonably well by the flame sensor criterion > 0.01 (evaluated at the droplet position). Please note that, similarly to Sacomano Filho et al. (2017), no coupling between carrier and liquid phase is performed. Besides, the impact of different flame thickening factors F, namely 1 (the reference scenario), 2, 5 and 10, is also investigated. Four different case setups are investigated here, each referring to a specific flame interaction angle . Alongside to the two edge cases = 0 • and = 90 • , portraying the cases where a droplet moves orthogonal or parallel to the flame front, the setups with the angles = 30 • and 60 • are analyzed. It should be noted that, for the case where a droplet moves parallel to a flame front, the droplet is inserted directly into the flame.
The results for the first two cases are shown in Fig. 6, which shows the diameter evolution of the droplet over its distance within the flame. These two edge cases can be reduced to one dimension, either flame orthogonal or flame parallel direction. First considering the flame orthogonal case = 0 • , the diameter evolution along the trajectory is shown in Fig. 6a. As previously stated, the heat and mass transfer is turned off outside of the flame, explaining the regions prior to (and succeeding) the flame where the droplet diameter is assumed to remain constant. The reference solution, that is the unthickened flame case is exhibited by the solid black line. Here, as expected, all correction approaches reproduce the reference solution for F = 1 (all solid lines are superimposed). However, it becomes obvious that the droplet evaporation is overpredicted when thickening is applied but the correction of droplet heat and mass transfer omitted (dashed and dotted black lines in Fig. 6a). The approach without correction is, as a result, not able to reproduce the correct flame front exit diameter. For this case, the other approaches, namely the standard, refraction and projection methods are all equally able to reproduce the reference flame exit diameter. Going over to droplet parallel movement case ( = 90 • ) exhibited in Fig. 6b, the behavior is different. Here, all approaches, independent of the thickening factor are reproducing the same diameter evolution as the reference solution, except the standard approach (blue lines). Thereby, the deficiency of the standard thickening correction becomes evident as one can observe that the diameter evolution, in other words the droplet vapor source term, is elongated by the thickening factor F, despite the flame being not thickened in that direction. In difference, the projection and refraction approaches are both able to reproduce the reference solution. Nonetheless, as previously outlined, while these approaches are reproducing equal results, they are grounded on two different procedures and at a closer look, they are in fact, very different. While the projection method attains the correct behavior by reducing the effective thickening correction factor F eff to unity, the refraction method does so by performing the previously explained trajectory correction. As outlined in Sect. 2, for the refraction approach, F eff is the same as for the standard one (see Table 3). Following this argumentation, the approaches must at least reproduce equal droplet-flame interaction times. This fact demonstrates that, while the diameter evolves along the same path for the refraction and projection correction, they do it at different times. In the case of droplet movement parallel to the flame front, the refraction approach is, differently from the projection correction and the method without correction, elongating the droplet lifetime by the factor F.
For the two-dimensional cases = 30 • and 60 • , the differences between the approaches become more distinct. For these cases, it is no more enough to take a look at one direction. Both, flame parallel and orthogonal direction, must be considered. The results are displayed in Fig. 7, showing the flame orthogonal coordinate (left) and the flame parallel one (right). First, considering the case where = 30 • in the top of Fig. 7, one can already observe differences for the flame orthogonal direction. Here, the standard and refraction approaches produce equal results and are able to reproduce the reference flame exit diameter. This is a remarkable point since, under the condition that a droplet exits the flame, i.e. that it has a flame front exit diameter, the standard approach already fulfills the exit-diameter criterion. This is independent of the droplet flame front parallel movement. However, it should be kept in mind that comparing the flame exit diameter of a droplet crossing a thickened flame front with the exit diameter of its non thickened pendant does not reveal the potential impact of the evaporation through the flame parallel direction. In view of complex spray combustion simulations, this is a crucial aspect since, when compared to experiments, it is not only important to observe similar effects, but also at similar positions. If flame front orthogonal effects dominate, the displacement resulting from the thickening procedure is expected to not surpass several flame thicknesses, which is in most cases insignificant compared to the geometric dimensions of the configuration. However, as flame front parallel effects become relevant, such a dislocation is expected to grow quickly. This is due to the elongation of the trajectory parallel to the flame front, which is portrayed best in Fig. 6b and also visible in the top right of Fig. 7. For this direction, the refraction correction yields a similar distribution of vapor release as the reference solution, which is caused by the trajectory correction. This is very different from the standard approach. As for the orthogonal scenario, the approach without correction yields much faster evaporation rates and the droplet quickly disappears. In between the other methods lies the projection approach, for which slightly smaller exit diameters than for the standard or refraction approach can be observed. It is worth noticing that for this approach, the diameter at the flame front exit is depending on the thickening factor. Both previously mentioned effect are not unexpected since the projection can be reformulated as combination of the standard approach ( i.e. F eff = F std eff ≠ ) and the one without correction ( F eff = F none eff = 1 ), which becomes more obvious by reformulating Eq. 13 to: with the flame orthogonal weighting factor w ⊥ = | | | (∇PV) The visible but small differences between the projection approach exit diameter and the one resulting from the standard or refraction approach are attributed to the minor flame parallel movement compared to its flame orthogonal component, yielding values of w ⊥ close to unity. Evidently, this is also For the configuration = 60 • in the bottom of Fig. 7, it turns out that the small orthogonal velocity yields a complete evaporation of the droplet before the flame front exit. As expected, the diameter evolution for the projection method is shifted towards the solution obtained without correction due to the decreasing ratio between the flame orthogonal and flame parallel velocity. This results in a stronger concentration of the vapor release in both, flame parallel and orthogonal direction. In the flame normal direction, as for the previous case, the vapor release is distributed equally along the thickened flame for the standard and refraction correction while similar deviations are obtained in the flame parallel direction.
To summarize, the diameter evolution of a droplet crossing a thickened flame along its flame normal and flame parallel trajectory is strongly depending on the manner heat and mass transfer is treated. Following the presented analysis, the refraction correction clearly shows the most promising results as it conserves the flame front exit droplet diameter of a given non-thickened flame and at the same time keeps the droplet evolution across the flame parallel direction unaltered. It now remains to investigate how the approaches perform in a complex spray jet flame. The following sections aim at clarifying the impact of the various approaches in such configuration.

Influence of the Flame Sensor Formulation
The objective of this section is to briefly justify the choice of the flame sensor formulation. From Fig. 8, which exhibits temperature contours with isolines of the flame sensor superimposed with computational parcels, one can observe some differences between the formulation used in this work (i.e. Eq. 8) and a standard formulation commonly used for purely gaseous combustion initially proposed by Durand and Polifke (2007), where = 16(c(c − 1)) 2 . Hereby c denotes a normalized reaction progress variable ranging from 0 (unburnt) to 1 (burnt). While a different formulation for the flame sensor may impact combustion characteristics for single phase flows, it certainly influences mixing and reaction processes in the presence of a liquid phase if a thickening correction is applied . One can clearly see that on the left of Fig. 8, the latter approach detects the flame much earlier than the temperature rise, thus creating a region prior to the flame front where carrier temperatures are still low, but where thickening is already performed. Taking a look at Eqs. 11 and 12, it becomes clear that this circumstance yields a correction of droplet heat and mass transport in low temperature regions. Additionally, since the principal movement direction of the parcels is parallel to the flame front, droplets tend to dwell longer in that region. Both factors are acting towards a reduction of heat and mass transfer between liquid phase and carrier in direct proximity of the flame front. Evidently, this also yields an altered mixture development of the carrier in flame front parallel direction. Following this line of thought, one may come to the conclusion that the reduction of vapor release in this region may also influence the flame propagation due to this altered mixture evolution found in the carrier. In contrast to the formulation by Durand and Polifke (2007) and as exhibited on the right of Fig. 8, the onset of the thickening is much closer to the temperature rise for the definition of used in this work (Eq. 8). Thus, by using this formulation, it is likely to substantially reduce the previously described effect.

Approach Assessment in a Turbulent Spray Flame
In this section, the impact of the modeling approaches listed in Table 3 is now assessed for a turbulent case, namely the EtF6 spray flame. Each method refers to a different way of addressing the droplet heat and mass transfer in the presence of a thickened flame. Firstly, instantaneous contour plots of various scalar fields are shown and differences between the approaches discussed. Subsequently, radial profiles of temporal averages and rms values for velocity and scalar fields are addressed for carrier and liquid phases. It should be indicated that for all simulations performed, the values of the thickening factor F do not exceed 10.
The instantaneous contour plots of temperature, mixture fraction as well as OH and CO mass fraction fields for the standard, projection and refraction correction approaches are shown in Figs. 9, 10 and 11, respectively. Here, it should be mentioned that the species mass fractions are directly extracted from the thermochemical lookup table, which may introduce errors for species with higher chemical timescales (Ganter et al. 2017). However, the same chemistry treatment for all approaches makes at least a qualitative comparison possible. In direct proximity of the burner exit, no differences between the methods can be observed. In this region, the flame is predominantly found concentric to the main jet and does not propagate towards the centerline as the mixture of the central jet is below the lower flammability limit. As axial distance increases, the mixture evolves towards a flammable one since more fuel vapor is released from the droplets. This aspect summed up with the higher interaction of flame and turbulence allows the propagation of the reaction zone towards the jet centerline for all approaches. At this point, first differences caused by the various methods become evident. While the standard and refraction approaches extend the central non-reacting core up to almost 30 jet diameters, the projection procedure seems to produce a notably shorter cold core ( ≈ 20D ). This strong disparity also results in a seemingly shorter flame for the projection approach. A second observation is that all three approaches induce a different distribution of evaporated liquid shown through the instantaneous contour plots of mixture fraction. First taking a look at Fig. 10b for the projection correction, one can see some pockets of rich mixtures occurring at lower axial positions. While these richer regions are also present for the standard approach in Fig. 9, they only appear at higher distances from the burner exit and seem less intense. Also to mention is the large region with noticeably higher mixture fraction around x∕D = 25 for the projection approach which is on one hand evidently narrower and appears at higher axial positions on the other, for the standard approach. A completely different scenario is exhibited in Fig. 11b for the refraction correction, where distinct layers of fuel rich mixtures can be observed at lower axial positions. These layers, which are spread over almost 15 jet diameters, clearly indicate strong entrance of fuel through droplet evaporation in these regions. The corresponding contour plots of CO and OH illustrate that in these richer regions a high concentration of partially oxidized products are present. The rather thin OH zones suggest not only a double flame structure but seem to promote the development of partially-premixed or non-premixed combustion regions at these lower axial positions. This is confirmed by analyzing the volumetric ratio of premixed to non-premixed combustion at lower axial positions ( x < 15D ) which are presented in Table 4. This ratio is obtained by evaluating the flame index conditioned on regions of reactions, that is where the progress variable source term ̇P V exceeds 10% of its maximal possible value. The flame index follows the definitions used in Yamashita et al. (1996): Thereby, positive values of indicate a premixed and negative values non-premixed combustion. A total of 10 independent snapshots were used to allow meaningful statistics. As can be deduced from Table 4, while the ratios for the standard and projection approach are almost similar, the refraction correction shifts this ratio in direction of non-premixed combustion. The differences in combustion regime distribution are also present at higher axial positions, where an inverse behavior can be observed. In this region, while the premixed Fig. 9 Instantaneous contour plots of (a) temperature, (b) mixture fraction, (c) OH and (d) CO using the standard thickening correction regime is still dominant for all correction types, the refraction approach yields higher volumetric ratios of premixed to non-premixed combustion than the other methods.
Going back to the standard and projection approaches shown in Figs. 9 and 10, the OH and CO profiles illustrate that few pockets of partially oxidized products, suggesting combustion above stoichiometry, evolve behind the flame. This is more evident in Fig. 10d for the projection approach. A qualitative comparison with OH measurements (not shown here) found in Clean Combustion Research Group Database (2019) and Gounder (2009) (see Figs. 7.11-7.13 for this reference) revealed that, aside from the reduced flame wrinkling due to the artificial flame thickening, the standard and projection method are evenly Table 4 Volumetric ratio of premixed ( > 0 ) to diffusion flame ( < 0 ) regime up-and downstream of the axial positions x = 15D conditioned on regions of reaction, i.e. where ̇P V exceeds 10% of its maximal possible value Projection Refraction x < 15D 9.76 9.31 6.37 x > 15D 2.60 2.07 3.53

Fig. 10
Instantaneous contour plots of (a) temperature, (b) mixture fraction, (c) OH and (d) CO using the projection thickening correction able to predict the flame structure at lower axial positions. However, at higher distances, the projection correction results better match the experiments due to the earlier flame propagation towards the jet centerline. Differently, the OH structures observed in Fig. 11c for the refraction approach hardly resemble the experimental results. It becomes clear that the manner in which droplets interact with a thickened flame has a strong influence on overall flame structure and combustion regime distribution. The disparities between the approaches become comprehensible when considering the different correction schemes themselves. In the projection approach, the relative movement directions of droplet and thickened flame (i.e. their orientation) is taken into account. Since the main movement direction of the particles is parallel to the flame front, only little or no correction of heat and mass transfer is applied.
This effect is illustrated in Fig. 12, which shows on the left hand side the instantaneous distribution of computational parcels at the axial plane x = 5D superimposed with isolines of the flame sensor , on the other hand a probability histogram of the effective correction factor calculated for the parcels at this axial plane (collected over at least 0.15 s). The parcels are colored with the angle = tan(u p,ax ∕u p,r ) , which represents the angle between radial and axial momentum of the parcels. As the droplets travel mainly in axial direction, the high angles observed in Fig. 12a are not unexpected. More impressive is however the impact on the computation of the effective thickening factor  Fig. 12b. The high axial velocity of the parcels combined with a flame front that is predominantly radially oriented leads to a significant reduction of the effective thickening factor for the projection correction.
This promotes faster evaporation and a more intense evolution of the mixture towards flammability at lower axial positions for the projection approach. This is also valid for the case without thickening correction, for which details are presented in the next part of this section. In difference, the standard method does not account for any orientation between parcels/droplets and the thickened flame, which leads to correction factors higher than for the projection approach, hence reducing heat and mass transfer for droplets in the thickened flame region. The consequence is, as previously discussed, a mixture changing far less through the longitudinal direction. This is elongating the cold spray core to greater distances from the burner exit, which also affects the flame length. At the same time, it leads to different mixture compositions and resulting distribution of partially oxidized products, mainly due to evaporation taking place at higher distances from the burner exit. For the refraction approach, the trajectory correction of the flame front parallel motion causes a considerably richer mixture at lower axial positions compared with the other approaches, which goes hand in hand with the stronger concentration of the vapor release in flame front parallel direction. However, it seems that this approach induces an elongation of the cold central jet core compared with the projection correction. This at first sight unexpected effect will be further discussed in the remainder of this section.
Next, radial profiles of temporal mean and rms of different carrier phase quantities are presented in Figs. 13 and 14. First the carrier phase axial and radial mean and rms velocity shown in Fig. 13 are considered. The simulation results are compared with results from experimental set A and B. Hereby, it should be mentioned that the comparison is not fully consistent since experimental carrier data were produced from statistics obtained from droplets below 10 μm . In addition radial velocity statistics are only available for experimental set B. All presented fields are normalized with the bulk jet velocity u bulk displayed in Table 2. At the lowest measurement plane, namely x = 0.3D , the distinct peaks of central jet and pilot can be observed for the averaged axial velocity. As distance from the burner exit increases, these peaks blur quickly into each other leading to the smooth profile observed at the axial position x = 10D . Up to that position, the velocity averages are evenly able to reproduce the experimental profiles. However, downstream of that position, first differences between the approaches arise. At x = 20D , the refraction approach profile indicates a broader region of higher velocities than the other approaches. However, the more pronounced differences can be observed at the last measurement plane, where the centerline values for the standard and refraction approach are noticeably higher. Considering the standard approach, the reasons for the deviations are twofold: First, since this approach tends to reduce evaporation, liquid penetration will increase and thereby increase momentum exchange with the carrier phase up to higher axial positions. This is especially pronounced at the centerline. Secondly, as previously observed in Fig. 9, the non reactive core may exist up to axial distances slightly below x = 30D . A flame front in that region will clearly lead to thermal expansion resulting in a higher axial velocity (see Fig. 9). The combustion reaction taking place at higher axial positions is also responsible for the higher velocities observed for the refraction approach. However, the causal sequence leading to the observed behavior is very different. Here, the trajectory correction coming with this method shifts the residence probability of parcels -at a fixed axial position -towards higher radial positions compared to other approaches. This, combined with the absence of oxidizer, allows the development of the previously mentioned fuel rich layers at lower axial positions (Fig. 11). The mixing of these oxidizer-lacking-zones with coflow air at higher axial distances from the burner exit promotes combustion reaction yielding thermal expansion and an acceleration of the flow field. In contrast, the projection approach and the one without correction are better reproducing the experimental data. In regard of the radial velocity averages, all approaches are in the same order of magnitude as the experiment. However, some differences are also noticeable, which are most pronounced at x = 10D , where the peak values of the radial velocity average differ strongly. These higher peak values for the projection approach and the one without correction clearly indicate strong combustion reactions in that region. Looking at rms profiles of axial and radial velocity, one can perceive that, besides centerline values being underestimated at the lowest measurement plane and an overprediction of radial velocity rms at x = 30D for all but the refraction approach, the simulation results reproduce the experimental trends satisfactorily. In fact, all approaches are performing similarly with exception of the refraction correction which generates considerably smaller rms values around the centerline repeatedly. A possible explanation for that is the reduced interaction of droplets with the shear layer that surrounds the main jet stream, and which is promoted by the droplet trajectory modification performed for this approach. As a result, droplets interact for a shorter time period with this high turbulence production zone. Going over to the scalars presented in Fig. 14, where only the excess temperature T − T 0 was measured for experiment A, the impact of the various models becomes more apparent. As previously revealed in the instantaneous contour plots, the mixture evolution shows to be strongly depending on the approach to correct heat and mass transfer. This hypothesis is now investigated quantitatively. At the lowest plane where x = 0.3D the distinct levels of temperature and mixture fraction corresponding to central jet, stoichiometric pilot flame and primary coflow are clearly visible. This together with the fluctuations of the turbulent flow give rise to the two small peaks in the rms profiles of mixture fraction shown in  Fig. 14c. As distance from the burner exit increases, the evaporation of the liquid phase leads to an increase of mixture fraction which is expected to be dominant in proximity of the flame front. In this context, it is also important to mention that the general overestimation of temperatures is influenced by the adiabatic assumptions made in this work (see Sect. 2). At x = 10D , some differences between the approaches become noticeable. Starting with the average temperature profile, one can observe that all simulations tend to underestimate the temperature at the centerline, indicating a longer non-reacting sprayjet-core than in the experiments. This attests that flame propagation towards the jet center is underestimated in the simulations at lower axial distances. In our view, this is caused by the combined effects of an under prediction of turbulence-flame interaction and by the negligence of differential diffusion transport. As discussed at the end of Sect. 3, both could increase the resulting turbulent flame speed, thus favoring a higher entrance of the flame into the jet core. In regard of the different approaches, it can not only be observed that the averaged temperature rise takes place at different radial positions, but also that the flame structure varies from approach to approach. This is most evident for the refraction correction, where the double peak present in the average temperature profiles at x = 10D and x = 20D indicates a double flame structure. This confirms the raised hypothesis from the previous analysis of instantaneous contours (see Fig. 11 and subsequent discussions). In this regard, the other approaches better agree with the experimental profiles. At x = 10D , the temperature rise, which is closer to the centerline for the projection approach and the one without correction, indicates a different flame position. The lowest radial positions can be observed for the case without thickening correction as F eff has the smallest value for this approach. This is also confirmed by the averaged and rms mixture fraction profiles. Taking a look at the mean profiles, a similar behavior as for the averaged temperature can be observed: the projection method and the one without correction yield a mixture fraction rising closer to the centerline due to the different flame position. This is closely connected to the lower values of F eff compared with the other approaches. This is even more apparent in the variance charts, where undulated profiles with two distinct peaks are visible for all approaches. While the most outer peak, corresponding to the pilot-coflow mixing layer is similar for all modeling strategies, the most inner peak, representing a blurring of jet-pilot mixing layer and evaporation zone, differs in its strength and position. For the standard and refraction approaches, a small increase of the rms is perceived at lower radii. As the radial distance increases, a peak can be observed around r = 0.75D for the standard approach. This peak is much more distinct as well as shifted to higher radii for the refraction approach. Such radial shift is also apparent in the mixture fraction average profile. This is not surprising when considering that the trajectory correction results in higher radial positions of droplets (at a given axial position) compared to the other approaches. This shift clearly affects the strongly connected vapor release and the carrier phase mixture evolution. For the two remaining approaches, namely no-correction and projection correction, the mixture fraction variance peak is much closer to the centerline, suggesting a strong interaction of turbulence and evaporated fuel. The disparities between the approaches reach their maximum at x∕D = 20 . The findings from the instantaneous contour plots analysis are confirmed by the temperature and mixture fraction profiles, where one can observe that: (1) For the projection approach and the one without correction, the temperature profile indicates combustion (or hot products) at the centerline, which better matches the experimental measurements. This is not the case for the standard and refraction correction approach indicating a longer non-reacting cold core.
(2) Centerline mixture fractions are shifted towards richer mixtures for the projection method compared to the standard and refraction approach and are highest when no correction is applied. Here, the highest deviations between the proposed projection method and the one without correction emerge, which is not unexpected as a considerable part of the droplets move orthogonally to the flame front leading to higher correction factors F eff in this region. Finally, at the last measurement plane, a considerably higher mixture fraction average can be observed for the standard approach. This is certainly caused by the stronger evaporation rates resulting from a greater spray penetration for this case. Also interesting is the clearly lower mean mixture fraction for the refraction correction which is coupled to the lower centerline temperatures for this method. However, this observation must also be attributed to the trajectory correction of the droplet, since as a droplet within the central jet travels along its trajectory, the probability to cross a thickened flame region increases with the axially traveled distance. This in turn increases the probability of the droplet to experience a radial dispersion. The velocity statistics as well some characteristic scalar properties for the liquid phase are presented in Figs. 15 and 16. The statistics shown are based on number averages of the real droplets (as each computational parcels represents a specific number of real droplets with the same properties). In regard of the velocity profiles shown in Fig. 15, a first observation is that all modeling approaches are able to represent the main trends adequately. However, considering the last measurement plane, it appears that the standard and refraction approach overestimate axial averages, which is similar to the behavior observed in Fig. 13a for the carrier phase. The other approaches deliver better results in that regions. The similarity of the droplets averaged radial velocity profiles in Fig. 15c suggests that the radial velocity is not strongly connected to the correction approach. With respect to the variance profiles, a similar picture as for the carrier phase velocity profiles emerges, namely that the rms of the axial velocity component (Fig. 15b) are underestimated at the lowest measurement plane. This is resulting from the droplet injection strategy, i.e. that computational parcels are injected with the carrier phase velocity interpolated onto the injection position. However, while all approaches differ from the experimental results, the projection method and the one without correction seem to deliver smaller offsets than the other approaches. The radial velocity rms profiles also indicate a better agreement with the experiments for the two former approaches, with exception of the last measurement plane. Radial profiles at different axial positions for liquid phase scalar properties, namely Sauter Mean diameter D 32 (SMD), mean diameter D 10 and liquid volumetric flux are shown in Fig. 16. Excellent results are obtained for all properties at the lowest measurement plane, which is expected since this is the liquid injection plane in the simulations (see Sect. 2). Taking a look at the characteristic diameters displayed in the first two columns, the simulations have a tendency to underpredict their profiles at all positions and for all approaches with exception of the injection plane. It should be noticed that similar deviations have been observed in the comparative study of Heye et al. (2014) for all modeling frameworks compared. The offset goes hand in hand with the simulations predicting longer cold cores (see Fig. 14) than in the experiment. This is because smaller droplet can subsist up to greater axial distances without the presence of a flame or hot products at the centerline. Another possible reason is a mismatch between simulations and experiments for the injection velocity of the liquid phase. Since droplets are injected with the carrier velocity, the velocities of large droplets are likely to be overestimated, as these droplets are in general slower than the carrier phase. This higher injection velocity for large droplets leads to a preferential residence of smaller droplets around the centerline because larger droplets are able to exit 1 3 the central jet due to their higher inertia. The hypothesis is confirmed by the higher values perceived in the profiles at greater radial positions. This effect and the previously discussed longer non reacting jet core are both acting towards shifting the characteristic diameters around the centerline to smaller values. Nevertheless, one can observe that the projection correction as well as the approach without correction reduce the gap between simulations and experiments. This is most evident at x = 20D for the SMD and at x = 30D for the mean diameter. The differences between the approaches are even more pronounced when considering the liquid volumetric flux shown in Fig. 16c. At x = 10D , the profiles are not far apart. However, the slight peak around r = 1D for the refraction approach is remarkable as it does not appear for the other approaches at this radial position. This peak combined with the mixture fraction profile shown in Fig. 14b suggests that the deviation of the particles performed by this approach predominates at lower axial regions. Indeed, the aforementioned peak disappears when passing to higher axial distances. At these positions other differences between the correction approaches emerge, which are strongest around x = 20D . At this location, the standard approach yields centerline volumetric fluxes almost twice as high than in the experiments or when no thickening correction is performed. The refraction approach is also overestimating the centerline liquid volumetric flux. In between lies the projection approach which seems to slightly overestimate the volumetric flux in that region. This is however consistent with a reaction taking place at the centerline at higher axial distances than in the experiments. The differences between the projection approach and the one without correction also suggests that a considerable amount of droplets are moving with a significant velocity component in flame front orthogonal direction yielding higher values for F eff . This in turn affects the evaporation and leads to the differences between the two methods. At the last measurement plane, where x = 30D , the standard approach fits the experiments most accurately. However, the previous results as well as the overall order of magnitude for the volumetric flux at this axial positions clearly reduces its significance for interpretation. To summarize, the previous discussion clearly showed that the way in which droplet interacts with a thickened flame has a meaningful impact not only on global flame characteristics as flame length or spatial distribution of flame modes, but also on flow dynamics and scalar distribution of the carrier and liquid phase. The differences between the standard method and the introduced correction approaches, namely projection and refraction correction, are already evident at lower axial positions, which can be observed in instantaneous contour plots of carrier scalar fields as well as in radial profiles of temporal averages for various fields. The refraction correction, which performs best in a simplified one-way coupled configuration, can by far not be considered adequate when applied to a complex spray flame. Differently, the approach without thickening correction and the projection correction appear to be most promising. However, the contrast between both approaches is not pronounced, especially when looking at the radial profiles upstream of x = 20D . This is expected, since droplet movement is predominantly parallel to the flame front in that region. Differently, around x = 20D , a considerable part of the droplets move orthogonally to the flame front, resulting in higher correction factors F eff in that region. This yields the observed deviating results at higher axial positions which are most evident for mixture fraction averages (Fig. 14b) and for the liquid volumetric flux shown in Fig. 16c.

Conclusions
Four different strategies to describe heat and mass transfer of droplets in the presence of an artificially thickened reaction zone have been investigated. Thereby, the modeling combined the FGM method with the ATF approach, while the multiphase flow was treated within an two-way coupled Euler-Lagrange framework.
In a first stage, it has been shown by means of simplified two-dimensional tests that the approaches found in literature are not sufficient to describe the various scenarios of droplet-flame interaction that a particle/droplet may experience. The weaknesses of the two methods currently employed in literature, that is the approach where heat and mass transfer are corrected directly by the flame thickening factor (standard approach) and the strategy omitting any correction, were unveiled when considering the two edge cases of a droplet traveling either parallel or normal to a flame front. These shortcomings motivated the introduction of two additional approaches, referred to as projection and refraction correction. Both methods aim at including the orientation of a thickened flame front with respect to a droplet crossing it in the correction procedure of heat and mass transfer. While the projection approach attains the desired behavior by scalar modification of the effective thickening factor F eff in the droplet temperature and mass evolution equations, the refraction strategy accounts for a droplet trajectory correction. Both approaches are able to correctly predict the evolution of the droplet diameter for the two limiting conditions where at least one of the two approaches found in literature failed. For scenarios in between the two edge cases, the refraction approach performed best followed by the projection methods.
In a next step, the impact of the various approaches and their prediction ability have been evaluated in a complex spray flame, the configuration EtF6 of the Sydney spray burner. Here, a different behavior than for the two-dimensional test configuration was observed. Interestingly, the refraction approach not only fails to reproduce the experimentally determined flame structure and regime, but delivers also worse carrier and liquid statistics. In retrospective, it becomes clear that these discrepancies are caused by the trajectory correction which yields a radial deflection of the parcels. This is closely connected to the greater mixture fraction values, which occur at higher radial positions for the axial positions x = 10D and x = 20D . The projection approach appeared to reproduce the experimental trends and profiles more accurately than the standard approach and the refraction method. Thereby, the projection approach leads to a better prediction of the flame propagation towards the center of the spray jet, better estimations for the evolution of the disperse liquid phase, shown by means of radial profiles of liquid volumetric fluxes and characteristic diameters at several axial positions, compared with the standard approach used in Sacomano Filho et al. (2017), Cheneau (2019), Cheneau et al. (2019). The impact on carrier and liquid phase velocity statistics is especially pronounced at higher axial distances from the burner exit. The disparities between the proposed approach and the one without correction (i.e. F eff = 1 ) are, as expected, most pronounced in regions where flame-normal droplet movement is significant.
It can finally be concluded that: (1) The manner in which the interaction of droplet with a thickened flame is treated does not only strongly affect global flame quantities, for instance the length of the cold central spray core or the flame length, but also the prediction of the flame regime as well as carrier and liquid phase statistics. (2) By taking into account the relative orientation of flame front and droplet movement by means of the novel projection method, the overall consistency of the modeling framework has been improved.
(3) The analysis showed that a consideration of a simplified configuration is insufficient to fully uncover the performance of the different approaches. (4) In particular, the projection approach appears most suitable for treating the ATF technique in a complex spray flame. How the various approaches perform in a different flame configuration (e.g a swirling flame (Shum-Kivan et al. 2017)) is left for future research work.
From Fig. 17, some differences are perceivable but the grid dependency does not seem dominant when going to higher spatial resolutions. This can be deduced from the radial profiles of mean and rms carrier axial velocity (Fig. 17a) as well as averaged temperature values (Fig. 17b) obtained for both grids, which are very close to each other. Regarding the impact of the resolution of liquid phase statistics, similar conclusions can be drawn, as can be seen from the liquid volumetric flux profiles. Equivalent results are also obtained for other physical quantities and axial positions. This confirms that both grids are appropriate to describe the physics of the configuration. Consequently, in order to allow accurate longterm statistics within a reasonable computational effort, the grid with 4.3 million control volumes is selected for further analysis.