Evolution of Surface Density Function in an Open Turbulent Jet Spray Flame

A three-dimensional Direct Numerical Simulation of an open turbulent jet spray flame representing a laboratory-scale burner configuration has been used to analyse the statistical behaviours of the magnitude of reaction progress variable gradient ∇c\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left| {\nabla c} \right|$$\end{document} [alternatively known as the Surface Density Function (SDF)] and the strain rates, which affect its evolution. The flame has been found to exhibit fuel-lean combustion close to the jet exit, but fuel-rich conditions have been obtained further downstream due to the evaporation of fuel droplets, which leads to the reduction in the mean value of the SDF in the downstream direction. This change in mixture composition in the axial direction has implications on the statistical behaviours of the SDF and the strain rates affecting its evolution. The mean value of dilatation rate remains positive, whereas the mean normal strain rate assumes positive values where the effects of heat release are strong but becomes negative towards both unburned and burned gas sides. The mean values of dilatation rate, normal strain rate and tangential strain rate decrease downstream of the jet exit. However, the mean behaviours of displacement speed and its components do not change significantly away from the jet exit. The mean values of normal strain rate arising from flame propagation remain positive and thus act to thicken the flame. The mean tangential strain rate due to flame propagation (alternatively the curvature stretch rate) remains negative throughout the flame at all axial locations investigated. The mean effective normal strain rate assumes positive values throughout the flame and it increases in the downstream direction for the present case, which is consistent with the reduction in the peak mean value of the SDF in the axial direction. The mean effective tangential strain rate (alternatively stretch rate) assumes negative values throughout the flame at all axial locations.


List of symbols
Angles between ∇c and the eigenvectors associated with e , e and e ̇c Reaction rate term in the c transport equation 1 3 (a) to identify the statistical behaviour of |∇c| and its evolution for an open turbulent jet spray flame and provide physical explanations for these behaviours. (b) to demonstrate and explain the roles played by flame propagation on the evolution of |∇c| in turbulent spray flames (c) to provide modelling implications of the above findings.
The rest of the paper is organised as follows. The mathematical background and numerical implementation pertaining to this analysis are presented in the next section. This will be followed by the presentation of results and their discussion. The main findings are summarised and conclusions are drawn in the final section of this paper.

Considered Case and Computational Configuration
In this study, the DNS configuration corresponds to the experimental investigation of the Ethanol spray EtF3 flame of Gounder et al. (2012). The spray and carrier gas are injected from a central jet nozzle ( D j = 10.5 mm ) with a bulk velocity U j = 24 m s −1 surrounded by a coaxial pilot annulus ( D p = 25 mm ; U p = 11.6 m s −1 ; T p = 2493 K ) and an air coflow ( U c = 4.5 m s −1 ). The pilot is composed of the fully burned products of a stoichiometric mixture of 5.08% Acetylene (C 2 H 2 ), 10.17% Hydrogen (H 2 ) and 84.75% air by volume and provides the heat necessary for the evaporation of the fuel droplets. The flame is stabilised in the shear layer formed between the pilot and inner jet streams. The mass flow rate of liquid Ethanol in the jet is 45 g/min, however, amongst the polydisperse droplet formed by the nebulizer, some droplets evaporate prior to reaching the nozzle exit, giving a partially gaseous fuel at the exit. The Ethanol mass flow rates at the nozzle exit are 14.3 g/min for the gaseous phase and 30.7 g/min for the droplets, resulting in a gaseous equivalence ratio of 0.85.
The DNS data used in this analysis has been reported earlier by Kurose (2018, 2019). The simulation has been conducted using the FK 3 code (Kitano et al. 2014a, b;Turquand d'Auzay et al. 2019;Ahmed et al. 2020) employing a two-step chemical mechanism. This chemical mechanism accurately captures the equivalence ratio dependence of unstrained laminar burning velocity (Westbrook and Dryer 1981). Although simple chemistry is used in this simulation, the configuration of this study is relatively complex in terms of DNS and took more than 1,100,000 CPU-hours. It is worthwhile to note that the SDF statistics obtained based on DNS with single step chemistry (Chakraborty and Klein 2008;Chakraborty et al. 2013) are found to be in good agreement with the corresponding statistics obtained from detailed chemistry DNS data (Chakraborty et al. , 2013. Therefore, it is expected that the simplification in terms of chemical mechanism should not have any major influence on the conclusions drawn from the current analysis. Moreover, it has been demonstrated elsewhere Kurose 2018, 2019) that the results obtained from DNS are consistent with experimental findings of Gounder et al. (2012). Details of the numerical implementation have been presented in detail elsewhere Kurose 2018, 2019;Turquand d'Auzay et al. 2019) and interested readers are directed towards these studies. Whilst these details will not be repeated here but a brief overview will be provided. The standard conservation equations for mass, momentum, energy and species for reacting flows are solved for the Eulerian gaseous phase. The standard transport equations of position, diameter, momentum, and energy are solved for individual liquid fuel droplets in a Lagrangian framework Kurose 2018, 2019;Kitano et al. 2014a, b;Turquand d'Auzay et al. 2019). Coupling between the Eulerian and Lagrangian phases is achieved by the Particle-Source-In-Cell approach Kurose 2018, 2019;Kitano et al. 2014a, b;Turquand d'Auzay et al. 2019) where particles are modelled as spherical point masses. A polydisperse spray with a diameter distribution matching that of the experiment is injected based on a log-normal distribution with diameters ranging from 1 μm to 80 μm (the most likely diameter is ∼ 20 μm ). Details of the particle diameter distribution considered have been discussed in detail elsewhere and interested readers are directed to Pillai and Kurose (2018) for further information. As this is a dilute flame with an inflow droplet volume fraction of ∼ 5 × 10 −4 , the collisions and break-up have been neglected.
A domain of 94D j × 49D j × 49D j is used and is discretised by a non-uniform 1160 × 400 × 400 Cartesian grid. A large stretching is applied in all directions towards the boundaries to form absorbing zones that minimize reflection and contamination of the acoustic field near the jet Kurose 2018, 2019). Due to the coupling strategy between the Eulerian and Lagrangian phases, the minimum cell size must be larger than the droplet size to accurately capture the evaporation process Kurose 2018, 2019). To ensure an appropriate resolution of both the turbulence and the flame structure, the smallest cell size at the nozzle exit is taken to be Δx = 150 μm . Further details on the boundary conditions and mesh can be found in Kurose 2018, 2019;Kitano et al. 2014a, b;Turquand d'Auzay et al. 2019). The turbulence intensity decreases and the integral length scale increases as turbulence decays in the downstream direction. The Kolmogorov length-scale k increases continuously in the downstream direction from k ≈ 170 μm (i.e. Δx∕ k ≈ 1.35 ) at the nozzle lip, and the ratio of grid spacing to the Kolmogoroov length scale throughout the domain remains consistent with the DNS grid requirement recommended by Pope (2000). The turbulent Reynolds number Re t is in the range of 50-100 in the shear layer of the axial locations considered as discussed in detail in previous studies of this simulation (Pillai and Kurose 2019;Turquand d'Auzay et al. 2019). Additionally, a two-step global chemistry is used as the reaction model, which requires much less spatial resolution in comparison to more detailed chemical mechanisms employed elsewhere. The thermal flame thickness of the corresponding laminar premixed stoichiometric flame is 536 μm. Therefore, the grid resolution is sufficient inside the jet flame region.

Relevant Mathematical Background
The evaporation of droplets leads to mixture inhomogeneities that can be characterised by the mixture fraction , which can be defined as Y and W are mass fraction and the molar mass of element , respectively. A reaction progress variable c based on the oxidiser mass fraction can be defined, following several previous analyses Ozel Erol et al. 2018, 2019aFujita et al. 2013;Wandel 2013Wandel , 2014, as where Y O 2 is the oxygen mass fraction, Y O 2 ∞ is the oxygen mass fraction in the pure oxidiser stream and Y Eq O 2 is the equilibrium oxidiser mass fraction (i.e. Y ). Furthermore, the definition of c is consistent with a number of previous studies of turbulent stratified flames where variations in mixture fraction exist (Hélie and Trouvé 1998;Bray et al. 2005; 1 3 Chakraborty 2010b, c, c, d, f). The transport equation for c takes the following form Ozel Erol et al. 2019a): where D is the progress variable diffusivity. The first term on the right-hand-side of Eq. (1) arises due to molecular diffusion, the second term represents the reaction rate, the third term is the source/sink term due to droplet evaporation, and the final term is the crossscalar dissipation term due to reactant inhomogeneity (Malkeson and Chakraborty 2010a;Ozel Erol et al. 2019a). The definitions of ̇c , Ṡ ev and Ṡ c depend on the local value of . In Eq. (1), ̇c is defined in the following manner Ozel Erol et al. 2019a): is the stoichiometric mixture fraction. In Eq. (1), Ṡ ev is defined as Ozel Erol et al. 2019a): are the droplet source/sink terms in the fuel and oxygen transport equations, respectively, and Γ m is the source term in the mass conservation equation due to evaporation. In Eq. (1), Ṡ c is defined as Ozel Erol et al. 2019a): N ⋅ ∇c ) and tangential i.e. − 2 D m |∇c| components Ozel Erol et al. 2019a;Echekki and Chen 1999) where ⃗ N = −∇c∕|∇c| is the flame normal vector, m = 0.5 ∇ ⋅ ⃗ N is the curvature of a given c = c * iso-surface. The transport equation of c (see Eq. 1) can be rewritten in kinematic form as (Dopazo et al. 2015;Dopazo and Cifuentes 2016): where S d is the displacement speed, which is given by (Ozel Erol et al. 2019a;Wang et al. 2017): The SDF transport equation takes the following form (Dopazo et al. 2015;Dopazo and Cifuentes 2016;Ozel Erol et al. 2019a;Wang et al. 2017): are the flame normal and tangential strain rates, respectively, V j = u j + S d N j is the jth component of the flame propagation velocity, ∇.⃗ u is the dilatation rate, N j S d ∕ x j is the normal strain rate induced by flame propagation and 2S d m is the tangential strain rate due to flame propagation (or curvature stretch rate) (Dopazo et al. 2015;Dopazo and Cifuentes 2016;Candel and Poinsot 1990). Furthermore, it can be shown that (Dopazo et al. 2015;Dopazo and Cifuentes 2016): where Δx n is the distance between two neighbouring c iso-surfaces and a eff N is the effective flame normal strain rate (Dopazo and Cifuentes 2016;Ozel Erol et al. 2019a;Wang et al. 2017;Sandeep et al. 2018). For (a N + N j ( S d ∕ x j )) > 0 , the separation between c iso-surfaces increases but |∇c| decreases. Moreover, the evolution of an elemental flame surface area A can be written as (Dopazo and Cifuentes 2016;Candel and Poinsot 1990): where a eff T is the effective tangential strain rate (or flame stretch rate) (Dopazo et al. 2015;Dopazo and Cifuentes 2016;Candel and Poinsot 1990). Equations (7)-(9) suggest that the statistical behaviours of ∇.⃗ u , a N , a T , N j S d ∕ x j and 2S d m determine |∇c| and flame area generation in turbulent spray flames, which will be discussed in the next section.

General Flame Behaviour and Structure
The instantaneous iso-surface of reaction progress variable c = 0.75 coloured with normalised mixture fraction ∕ st at t sim = 0.024s is presented in Fig. 1a, which shows the variations in ∕ st due to droplet evaporation as well as due to the interaction between the flame, droplets and turbulent velocity field. The instantaneous field of normalised mixture fraction ∕ st on the central x-y plane is presented in Fig. 1b. It can be noticed in Fig. 1b that the evaporation process occurs in the mixing layer, which increases continuously from the nozzle lip and shows large values of up to ∕ st = 2.0 at x∕D j = 15 and ∕ st = 2.5 at x∕D j = 20 before decreasing due to mixing. Around the pockets of high fuel content created by the droplet evaporation (as shown in Fig. 1b), the burning occurs predominantly in a non-premixed mode as the hot fuel does not have the time to fully mix with the surrounding air, and leads to partial-premixing, which is characteristic of spray flames. This trend strengthens further in the downstream direction. The nature of the combustion (e.g. premixed, non-premixed) can be characterised by considering a Flame Index, FI . A modified Flame Index, FI , has previously been proposed by Briones et al. (2006) building upon the work of Yamashita et al. (1996) in the following manner: where the modified Flame Index, FI , assumes values between −1 and +1. A value of FI = −1.0 indicates a lean premixed mode of combustion, whereas a value of FI = 1.0 indicates a rich premixed mode of combustion and a value of FI = 0 indicates a diffusion flame. Probability density functions of the modified Flame Index, FI , at different reaction progress variable values (i.e. c = 0.1 , 0.3 , 0.5 , 0.7 and 0.9 ) at x∕D j = 2 , 6 and 10 are shown in Fig. 2a-c, respectively. It can be seen from Fig. 2 that close to the nozzle exit ( x∕D j = 2 ) the lean premixed mode of combustion remains dominant but rich premixed combustion becomes more dominant moving towards the burned gas side of the flame as droplets evaporate. However, moving further downstream (i.e. x∕D j = 6 and 10 ), greater contributions of non-premixed mode of combustion can be seen towards the unburned gas side of the flame (i.e. c = 0.1 ) due to a greater number of droplets beginning to evaporate far downstream of the jet exit. The non-premixed mode of combustion decreases (i.e. the PDF peak at FI = 0 decreases) with increasing c , as mixing progressively takes place within the flame.
The increased probability of finding fuel rich mixtures in the downstream direction, as observed in Fig. 1b and Fig. 2, is reflected in the variations of the mean values of conditional upon c , which are shown in Fig. 3a at x∕D j = 2 , 6 and 10 . The locations at x∕D j = 2 , 6 and 10 will be considered to analyse the behaviour of |∇c| and its evolution in the current study. Figure 3a indicates that likelihood of obtaining fuel-rich mixture increases towards the burned gas due to droplet evaporation at each of the locations considered, which is also consistent with Fig. 2. It should be noted that there are no occurrences of spontaneous ignition in the upstream shear layer. It was neither found in the simulation (Pillai and Kurose 2018, 2019) or nor was it reported in the experiments by Gounder et al. (2012). The absence of sufficient flammable mixture in the shear layer before a propagating  Fig. 3b at x∕D j = 2 , 6 and 10 as well as for the corresponding stoichiometric laminar unstretched premixed flame, showing that |∇c| × th,st decreases moving downstream. The inverse of the peak mean value of |∇c| provides a measure of flame thickness and thus Fig. 3b suggests that the flame thickens in the downstream direction. As the characteristic flame thickness increases with increasing ∕ st for the fuel-rich mixture, the flame thickening in the downstream direction is consistent with the variation of mixture composition in the gaseous phase with increasing x∕D j . This behaviour arises due to the smaller droplets being evaporated closer to the jet exit but larger droplets evaporating downstream leading to more fuel-rich conditions in these locations (see Figs. 1b and 3a). Figure 3b shows that as x∕D j increases, the peak mean value of |∇c| shifts towards the burned gas side of the flame. Equations (6)-(9) indicate that the statistical behaviour of |∇c| and its evolution are affected by the flame normal and tangential strain rates induced by fluid motion (i.e. a N and a T ) and flame propagation (i.e. N j S d ∕ x j and 2S d m ). Note that the additional strain rates induced by flame propagation (i.e. N j S d ∕ x j and 2S d m ) are dependent upon the variations of S d and its components (i.e. S r , S n , S t , S ev and S c ) which will be discussed next in this paper.

Displacement Speed S d Behaviour
The  Figure 4a-c exhibit that the mean S d remains positive with its magnitude increasing with c for the major part of the flame before reducing and potentially exhibiting negative values towards the burned gas side (see Fig. 4a, b). Figure 4 shows that the mean behaviour of S d does not change significantly with axial distance from the nozzle, and the mean value of S r acts as a leading-order contributor, exhibiting positive values increasing in magnitude towards the burned gas side. Figure 4a-c show that the mean contribution of S n also plays a leading-order role for all axial locations, exhibiting positive values towards the unburned gas side of the flame before exhibiting large negative values towards the burned gas side with a transition close to c ≈ 0.55 . However, the mean value of S t has been found to be small in comparison to the mean contributions of S r and S n at all axial locations, though its magnitudes have been found to increase with increasing x∕D j . Moreover, the mean values of S ev and S c remain small in comparison to the mean contributions of S r and S n at all axial locations. Accordingly, Fig. 4a-c indicate the behaviour of S d is principally determined by the competition between S r and S n irrespective of the location, which is consistent with the behaviour of S d and its components in canonical configurations Ozel Erol et al. 2019a).

Behaviour of the Fluid Dynamic Strain Rates
The effects of fluid-dynamic straining on the SDF evolution can be understood by analysing the statistical behaviours of ∇ ⋅ ⃗ u, a N and a T = ∇ ⋅ ⃗ u − a N , and the profiles of the normalised mean values of these quantities conditional upon c are presented in Fig. 5a-c, respectively. Figure 5a shows that the mean value of ∇ ⋅ ⃗ u remains positive with decreasing magnitudes exhibited with increasing x∕D j . The evaporation of the larger droplets further downstream leads to richer mixtures for the larger x∕D j , which reduces the strength of the thermal expansion and this is reflected in the reduced magnitude of the mean values of ∇ ⋅ ⃗ u . Figure 5b shows that a N exhibits predominantly mean positive values but exhibits negative mean values towards the burned gas side of the flame for x∕D j = 2 . For x∕D j = 6 and 10 , a N shows negative mean values towards the unburned gas before exhibiting positive mean values towards the middle of the flame, followed by negative values again towards the burned gas. The normal strain rate can be expressed as: where e , e and e are the most extensive (or positive), intermediate and the most compressive (or negative) principal strain rates, and , and are the angles between ∇c and the eigenvectors associated with e , e and e , respectively. It has been discussed elsewhere Ozel Erol et al. 2019a;Chakraborty and Swaminathan 2007;Ahmed et al. 2014;Kim and Pitsch 2007) that ∇c exhibits collinear alignment between the eigenvector associated with e (i.e. | cos | ≈ 1.0 ) when turbulent straining plays a dominant role, whereas a preferential alignment between ∇c and the eigenvector associated with e (i.e. | cos | ≈ 1.0 ) is observed when the strain rate induced by flame normal acceleration arising from thermal expansion overcomes turbulent straining. A preferential alignment between ∇c and the eigenvector associated with e ( e ) yields positive (negative) values of a N . The preferential alignment of ∇c and the eigenvector associated with e in the regions where the effects of thermal expansion are weak leads to negative (11) a N = e cos 2 + e cos 2 + e cos 2 mean values of a N towards the burned gas side for all axial locations shown here and also on the unburned gas side for x∕D j = 6 and 10. The strain rate induced by thermal expansion overcomes turbulent straining giving rise to a preferential alignment of ∇c and the eigenvector associated with e , which leads to positive mean values of a N within the flame where heat release effects are strong. The extent of ∇c alignment with the eigenvector associated with e increases in the downstream direction, and thus the effects of thermal expansion weaken in the axial direction (see Fig. 5a). Accordingly, the mean value of a N decreases in the in the axial direction and becomes more negative for larger axial distances. The relative magnitudes of the mean values of ∇ ⋅ ⃗ u and a N determine the mean behaviour of a T = ∇ ⋅ ⃗ u − a N . Figure 5c shows that the mean variation of a T exhibits similar qualitative behaviour at all axial locations (i.e. positive across the flame) but with a decreasing magnitude with increasing x∕D j . It is worth noting that the mean value of a T can be scaled as a T ∼ u � ∕ (Meneveau and Poinsot 1991) where u ′ is the rms turbulent velocity fluctuation and is the integral length scale. Turquand d' Auzay et al. (2019) have demonstrated that for the current cases, u ′ ( ) decreases (increases) with increasing axial distance resulting in reduced magnitudes of a T with increasing x∕D j .

Behaviour of the Strain Rates Induced by Flame Propagation
The profiles of the mean values of N j S d ∕ x j × th,st ∕S L,st and its components (i.e. N j (S r + S n )∕ x j × th,st ∕S L,st , N j S t ∕ x j × th,st ∕S L,st , N j (S ev + S c )∕ x j × th,st ∕S L,st ) conditional upon c are shown in Fig. 6a-d, respectively. Figure 6a shows the mean value of N j S d ∕ x j remains mostly positive but the magnitude decreases in the downstream direction, which is consistent with previous findings for round jet in premixed flames (Wang et al. 2017). The mean behaviour of S d is principally determined by the competition between S r and S n which is reflected in the mean behaviours of N j S r + S n ∕ x j and N j S d ∕ x j (see Fig. 6b). The mean contribution of N j S t ∕ x j is shown to be smaller than that of N j S r + S n ∕ x j , and the mean contribution of N j S ev + S c ∕ x j is shown to be negligible. As the mean value of N j S d ∕ x j remains mostly positive across the flame, it acts to increase the distance between two neighbouring c iso-surfaces (see Eq. 8).
The profiles of the mean values of 2S d m × th,st ∕S L,st and its components (i.e.

Modelling Implications
The statistical behaviours identified in the previous sub-sections are of significant importance from the perspective of modelling turbulent spray flames. It should be noted that the transport equation for the SDF |∇c| can be rewritten in the following form by expanding Eq. (7i): Furthermore, the transport equation of the generalised FSD (i.e. Σ gen = |∇c| ) (Boger et al. 1998;Klein et al. 2018) can be obtained by Reynolds averaging/LES filtering Eqs. (12i) and (12ii). By multiplying Eq. (12i) by 2 |∇c| one obtains: Additionally, algebraic manipulation of Eq. (13) leads to the transport equation of the SDR ( N c = D|∇c| 2 ) which can be written as follows ): Upon inspection of Eqs. (12)-(14), it is evident that the statistics of the strain rates a N , a T , N j S d ∕ x j , 2S d m , a eff N and a eff T are likely to play a significant role in the transport behaviour of |∇c| . As such, the statistical behaviours of these strain rates across turbulent spray flames need to be considered for the satisfactory modelling of the FSD and SDR in the case of turbulent spray combustion. It should be noted that similar behaviours of these strain rates have been observed in the current analysis as those found in turbulent spray flames in canonical configurations (14i)  Ozel Erol et al. 2019a), which suggests that the geometrical configuration might not play a major role in the statistical behaviour of SDF and its evolution. Moreover, the SDF statistics show qualitative similarities to the statistics obtained for a bluff body stabilised premixed flame representing experimental burner configurations (Wang et al. 2017;Sandeep et al. 2018). As such, given the found similarities, it is anticipated that the modelling methodologies proposed for turbulent premixed/stratified flames could be extended for turbulent spray flames. It should further be noted that a priori analyses for the development of models for the FSD or SDR approaches in the context of turbulent spray flames will utilise the statistics of the SDF and the strain rates which determine its evolution but the development of new closures for spray combustion requires further consideration of the modelling of unclosed correlations. Whilst such analyses are beyond the scope of the current study, the insights observed from this present analysis are important not only to determine the fundamental insights potential for the differences to those observed in turbulent spray flames of canonical configurations but also as a prelude to developing novel models for turbulent spray flames. Previous studies have investigated a priori modelling of the generalised Flame Surface Density Σ gen = |∇c| in the context of Reynolds Averaged Navier-Stokes and Large Eddy Simulations for turbulent premixed gaseous flames (Chakraborty and Cant 2007;Chakraborty and Cant 2009). The generalised FSD transport equation, which can be obtained by Reynolds averaging/LES filtering the SDF transport equation can be written in the following manner (Candel and Poinsot 1990;Cant et al. 1990;Hawkes 2000;Hawkes and Cant 2001;Cant 2007, 2009): where (…) s = (…)|∇c|∕|∇c| represents a suitable surface averaging, overline suggests Reynolds averaging/LES filtering and (� …) = (…) ∕̄ represents Favre averaging. In Eq. (15), the first term on the left-hand-side is the transient term and the second is the resolved convection term. The first term on right-hand-side of Eq. (15) is the turbulent transport term, the second term is the strain rate term, the third is the propagation term and the fourth is the curvature term. All of the terms on the right-hand-side of Eq. (15) are unclosed and require modelling. It should be noted that the FSD modelling is not a standard practice for droplet-laden combustion and one needs to understand the SDF statistics first before venturing into the FSD modelling. Moreover, FSD modelling is one of the possibilities that originates from the SDF statistics and one can also obtain SDR modelling related information from the SDF statistics. That said, one can further consider the available statistics to gain some insight into the future modelling of the FSD transport equation. For example, the combined contribution of the curvature and propagation terms in the FSD transport equation can be decomposed into three terms in the following manner Cant 2007, 2009;Hawkes 2000;Hawkes and Cant 2001): where P mean and C mean are the resolved contributions of the propagation and curvature terms, respectively, and C sg is the subgrid component. Expressions for P mean and C mean previously considered for turbulent premixed flames include Cant 2007, 2009;Hawkes 2000;Hawkes and Cant 2001): It has previously been noted that non-linear dependence of curvature on the curvature and propagation terms observed in turbulent premixed flames (Chakraborty and Cant 2007) can only be accounted for if the curvature effect on S d is taken into account in the modelling. Table 1 exemplarily presents the S d − m correlations at c = 0.7 for the different axial locations considered. It is evident from Table 1 that S d − m are negatively correlated and qualitatively similar behaviour has been observed for other values of c . This is consistent with previous analyses of turbulent premixed flames (Chakraborty and Cant 2007) and suggests that the correlation between S d and m = 0.5 N j ∕ x j cannot be neglected while modelling the FSD curvature term 2S d m s Σ gen , which is the area-weighted value of curvature stretch rate (or alternatively the area-weighted value of tangential strain rate induced by flame propagation). Moreover, the FSD curvature term 2S d m s Σ gen cannot be adequately approximated by C mean which disregards the interrelation between S d and m . Therefore, the modelling of C sg will be of critical importance for the purpose of modelling the FSD curvature term. Given the qualitative similarities between the curvature dependences of the flame displacement speed for premixed and spray flames, it can be expected that the modelling methodologies for C sg for turbulent premixed combustion have the potential to be applied for turbulent spray flames.
The expression for P mean (i.e. Eq. 17i) implicitly assumes that the displacement speed S d and the flame normal vector components (i.e. N 1 , N 2 , N 3 ) are uncorrelated. The correlations of displacement speed S d and the flame normal vector components (i.e. N 1 , N 2 , N 3 ) are exemplarily presented in Table 1 at c = 0.7 for each of the axial locations considered but similar qualitative behaviour has been observed for other values of c . It is evident from Table 1 that displacement speed S d and the flame normal vector components (i.e. N 1 , N 2 , N 3 ) are weakly correlated at each of the axial locations considered. This implies that P mean can be used to close the FSD propagation term − S d N j s Σ gen ∕ x j without any additional modelling provided that (S d ) s and (N j ) s are suitably modelled.
It should be noted that the jth component of the flame normal vector, (N j ) s , is given as Cant 2007, 2009;Hawles 2000;Hawkes and Cant 2001): In the context of the generalised FSD, Eq. (18) is an exact expression and has previously been proposed in the context of RANS and LES simulations Cant 2007, 2009;Hawkes 2000;Hawkes and Cant 2001;Cant et al. 1990). From the perspective of modelling of S d s , in the context of FSD modelling, the following relationship is often invoked (Hawkes 2000;Hawkes and Cant 2001): where 0 is the unburned reactant density and S b( ) is the laminar burning speed as a function of the local equivalence ratio . The variations of the mean values of ( S d |∇c|) and ( 0 S b( ) |∇c|) across c at each of the axial locations considered are shown in Fig. 9. It is evident from Fig. 9 that the mean value of ( 0 S b( ) |∇c|) remains different from S d |∇c| across c for all of the axial locations considered. Thus, the approximation S d s Σ gen ≈ 0 S b( ) Σ gen (Hawkes 2000;Hawkes and Cant 2001) may not be valid for turbulent spray flames. Finally, the statistics of a T in Fig. 5 and the associated discussion suggest that the modelling of the tangential strain rate term (a T ) s Σ gen in Eq. (15) depends on the accurate incorporation of dilatation rate ∇ ⋅ ⃗ u and the alignment of ∇c with local principal directions, as discussed in the context of Fig. 5. Such a priori modelling of the FSD is indeed beyond the scope of the current study. However, these aspects will form the basis of future investigations. The above discussion with regards to the insights into the statistical aspects of the study of the SDF evolution and the a priori modelling of the FSD indicates that further analysis is still required. The unclosed terms of the FSD transport equation (i.e. turbulent transport term, strain rate term, propagation term and curvature term) all require modelling. Such a priori modelling has not been considered based on an open turbulent jet spray flame in the presence of a shear layer and will contribute important novel insights. However, such analyses are not trivial and are, therefore, beyond the scope of the current study but will form the basis of future investigations.

Conclusions
A three-dimensional Direct Numerical Simulation database of an open turbulent jet spray flame representing a laboratory-scale burner configuration (Gounder et al. 2012) has been analysed to investigate the statistical behaviours of the SDF |∇c| and the quantities affecting its evolution (i.e. S d , ∇ ⋅ ⃗ u , a N , a T , N j S d ∕ x j , 2S d m , a eff N = a N + N j S d ∕ x j and a eff T = a T + 2S d m ). The flame in question has been found to exhibit predominantly fuelrich conditions and this tendency strengthens further in the downstream direction due to the evaporation of fuel droplets. This behaviour significantly affects the statistical behaviours of S d , ∇ ⋅ ⃗ u , a N , a T , N j S d ∕ x j , 2S d m , a eff N and a eff T . The mean SDF decreases downstream of the jet exit indicating progressive flame thickening in the downstream direction. For all axial locations considered, the behaviour of S d is principally determined by the competition between S r and S n with the mean contributions S t , S ev and S c all being small in comparison to leading order mean contribution of (S r + S n ) , which in turn affect the mean behaviours of N j S d ∕ x j and 2S d m . The relative magnitudes of ∇ ⋅ ⃗ u and a N result in positive mean values of a T with reducing magnitudes moving downstream of the jet exit. The normal strain rate induced by flame propagation N j S d ∕ x j has been found to exhibit positive mean values across the flame for all axial locations but the magnitude reduces in the downstream direction. The curvature stretch rate 2S d m has been found to exhibit negative mean values across the flame for all axial locations. The mean effective normal strain rate a eff N has been found to be positive across the flame but reduces in magnitude moving away from the jet exit. These positive mean values of a eff N indicate flame thickening behaviour in a mean sense. The effective tangential strain rate (or stretch rate) a eff T has been found to exhibit mean positive values across the flame close to the jet exit (indicating an increase in the flame surface area) but negative mean values across the flame in the downstream (indicating a decrease in the flame surface area). Similar behaviours of the strain rates affecting the SDF evolution have been observed in the current analysis as those found in turbulent spray flames in canonical configurations Ozel Erol et al. 2019a), which suggests that the geometrical configuration might not play a major role in the statistical behaviour of SDF and its evolution in the case of turbulent spray flames. However, further examination of different flame configurations will be required to test the veracity of this hypothesis. Furthermore, the SDF statistics discussed in this paper show qualitative similarities to the corresponding statistics obtained for a premixed jet flame and a bluff body stabilised premixed flame where both flames represent experimental burner configurations (Wang et al. 2017;Sandeep et al. 2018). Given these similarities, it is expected that the modelling methodologies proposed for turbulent premixed/stratified flames might be extended for turbulent spray flames, which will form the basis of future investigations.