**Ahmad Alqallaf 1, Markus Klein <sup>2</sup> and Nilanjan Chakraborty 1,\***


Received: 30 November 2018; Accepted: 11 January 2019; Published: 16 January 2019

**Abstract:** The effects of Lewis number on the physical mechanisms pertinent to the curvature evolution have been investigated using three-dimensional Direct Numerical Simulation (DNS) of spherically expanding turbulent premixed flames with characteristic Lewis number of *Le* = 0.8, 1.0 and 1.2. It has been found that the overall burning rate and the extent of flame wrinkling increase with decreasing Lewis number *Le*, and this tendency is particularly prevalent for the sub-unity Lewis number (e.g., *Le* = 0.8) case due to the occurrence of the thermo-diffusive instability. Accordingly, the *Le* = 0.8 case has been found to exhibit higher probability of finding saddle topologies with large magnitude negative curvatures in comparison to the corresponding *Le* = 1.0 and 1.2 cases. It has been found that the terms in the curvature transport equation due to normal strain rate gradients and curl of vorticity arising from both fluid flow and flame normal propagation play pivotal roles in the curvature evolution in all cases considered here. The net contribution of the source/sink terms of the curvature transport equation tends to increase the concavity and convexity of the flame surface in the negatively and positively curved locations, respectively for the *Le* = 0.8 case. This along with the occurrence of high and low temperature (and burning rate) values at the positively and negatively curved zones, respectively acts to augment positive and negative curved wrinkles induced by turbulence in the *Le* = 0.8 case, which is indicative of thermo-diffusive instability. By contrast, flame propagation effects tend to weakly promote the concavity of the negatively curved cusps, and act to decrease the convexity of the highly positively curved bulges in the *Le* = 1.0 and 1.2 cases, which are eventually smoothed out due to high and low values of displacement speed *Sd* at negatively and positively curved locations, respectively. Thus, flame propagation tends to smoothen the flame surface in the *Le* = 1.0 and 1.2 cases.

**Keywords:** Lewis number; flame curvature; iso-scalar non-material surfaces; turbulent premixed spherical flame

#### **1. Introduction**

Spherically expanding turbulent premixed flames are of fundamental importance in Spark Ignition (SI) engines and for understanding accidental explosions. Hence, they are often used as a canonical configuration for laboratory-scale experiments [1–9] and numerical investigations [10–31]. The role of mean flame curvature [20–22,26] on Flame Surface Density (FSD) [20,25] and Scalar Dissipation Rate (SDR) [25,27], which are central to determine the fuel burning rate, has been demonstrated in several previous DNS studies on ignition kernels [11,22] and spherically expanding flames [15,17–27]. Moreover, LES of spherical flames using the flame wrinkling factor [14,16], FSD [27,29] and

combined FSD-probability density function (PDF) [28] sub-grid closures, and Reynolds Averaged Navier–Stokes (RANS) simulations using various combustion modelling approaches [10,13,30] showed good agreement with experimental measurements. However, several analyses [20–22,24,27,30,31] demonstrated that there are significant differences between statistically planar and spherical flames, specifically in terms of flame propagation and fuel burning rate. Although these past studies provided important physical insights on spherical flames, they seldom considered thermo-diffusive effects arising from differential diffusion of heat and mass, characterised by Lewis number *Le* (i.e., ratio of thermal and mass diffusivities). The presence of thermo-diffusive instabilities augments the burning rate as demonstrated in experiments [4,5] and thus spherical flames in lean hydrogen–air mixture grow quicker compared to those in hydrocarbon–air mixture under statistically similar turbulent flow conditions. The necessity to include these effects was demonstrated by computing stoichiometric and fuel-lean hydrogen–air and methane–air spherical flames [30,31].

The flame wrinkling is often characterised in terms of flame front curvature distribution, which plays a key role in determining the local flame propagation behaviour. This is reflected in the correlation between displacement speed *Sd* and curvature *κ<sup>m</sup>* [21,22,24,32–41]. Moreover, in non-unity Lewis number flames the consumption speed *Sc* also demonstrates correlation with local flame curvature *κ<sup>m</sup>* [42]. Furthermore, tangential strain rate and curvature have been found to be negatively correlated for small turbulence intensities and the strengths of the correlations of tangential strain rate and flame speed with curvature weaken with increasing turbulence intensity. The flame curvature *κ<sup>m</sup>* is also known to affect the curvature and propagation terms in the Surface Density Function (SDF = |∇*c*| with *c* being the reaction progress variable) transport equation [43–50], which in turn influence the evolutions of FSD and SDR [27,29–31,50]. It has been demonstrated in several previous analyses [27,39,44,47] that the curvature dependences of displacement speed, temperature and SDF are influenced by the characteristic Lewis number *Le*, and thus it is expected that the curvature evolution is also likely to be affected by *Le*.

Pope [51] derived a transport equation for a parameterised surface and recently, Dopazo et al. [52] derived a transport equation of flame front curvature *κ<sup>m</sup>* and analysed the statistical behaviours of the different terms of this transport equation for passive scalar mixing without the effects of heat release. This analysis has been extended by Cifuentes et al. [53] to analyse the statistical behaviours of the terms in the curvature transport equation in a bluff-body stabilised turbulent premixed flame burner configuration using a flame-resolved high-fidelity simulation. However, the effects of characteristic Lewis number on the flame curvature evolution in a configuration with a non-zero mean curvature, as in the case of spherically expanding turbulent premixed flames, are yet to be analysed in detail. For example, the reasons for the differences in flame topology in terms of curvature distribution in response to the variation of *Le*, and the predominance of positively curved bulges and the presence of intermediate sharply negatively curved cusps between positively curved bulges for the *Le* < 1 flames have not been sufficiently explained in the existing literature. This distribution has strong implications on the augmentation of flame wrinkling and overall burning rate with decreasing Lewis number [39,42,44,54–58]. Thus, it is important to gain a better understanding of the curvature evolution in the presence of thermo-diffusive instabilities and their interrelation with flame shape and propagation in order to be able to derive high fidelity models, which can predict the non-unity Lewis number effects on premixed turbulent combustion. The flame curvature transport equation in Cartesian coordinates is a relatively new concept and the effects of thermo-diffusive instabilities on curvature transport are yet to be understood in detail. In order to address the aforementioned deficits in the existing literature, the statistical behaviours of the curvature transport have been analysed in this paper for statistically spherical flames with characteristic Lewis numbers *Le* = 0.8, 1.0 and 1.2. In this respect, the main objectives of this paper are:

1. To demonstrate the effects of Lewis number on the terms of the transport equation of flame curvature *κ<sup>m</sup>* in spherically expanding turbulent premixed flames.

2. To identify the mechanisms, which lead to the effects of thermo-diffusive instability (i.e., augmentation of burning rate and flame wrinkling and a positive correlation between local burning rate and flame curvature) in the flames with *Le* < 1.

The mathematical and numerical background pertaining to this analysis are presented in the next two sections of this paper. Following that, the results will be presented and subsequently discussed. The main findings will be summarised and conclusions are drawn in the final section of this paper.

#### **2. Mathematical Background**

In the current analysis, the chemical mechanism has been simplified by a single-step irreversible reaction so that the effects of characteristic Lewis number *Le* can be investigated in isolation, as done in several previous analyses [39–42,44,54–58]. The present analysis considers three characteristic Lewis numbers *Le* = 0.8, 1.0 and 1.2 following previous analyses [39–42,44,54–58]. Furthermore, the current analysis only focuses on Lewis number effects (i.e., differential diffusion of heat and mass) and not the differential diffusion between species.

In the case of non-unity Lewis number flames, the scalar field can be characterised by the reaction progress variable *c* and non-dimensional temperature *θ*, which are defined as:

$$\mathcal{L} = \frac{\mathcal{Y}\_{R0} - \mathcal{Y}\_{R}}{\mathcal{Y}\_{R0} - \mathcal{Y}\_{R\infty}} \tag{1}$$

$$\theta = \frac{(T - T\_0)}{(T\_{ad} - T\_0)} \tag{2}$$

where *YR* is the mass fraction of a suitable reactant which is used for defining the reaction progress variable. The subscripts 0 and ∞ are used to refer to values in the unburned gas and fully burned products, respectively. In Equation (2), *T*, *T*<sup>0</sup> and *Tad* denote the dimensional temperature, unburned gas temperature and adiabatic flame temperature, respectively.

The transport equation of the reaction progress variable *c*(**x**, *t*) is given by:

$$\frac{\partial \mathbf{c}}{\partial t} + \boldsymbol{\mu}\_{\dot{\jmath}} \frac{\partial \mathbf{c}}{\partial \mathbf{x}\_{\dot{\jmath}}} = \frac{1}{\rho} \frac{\partial}{\partial \mathbf{x}\_{\dot{\jmath}}} \left(\rho D\_{\varepsilon} \frac{\partial \mathbf{c}}{\partial \mathbf{x}\_{\dot{\jmath}}}\right) + \dot{\boldsymbol{w}}\_{\varepsilon} = \frac{1}{\rho} \frac{\partial}{\partial \mathbf{x}\_{N}} \left(\rho D\_{\varepsilon} \frac{\partial \mathbf{c}}{\partial \mathbf{x}\_{N}}\right) + D\_{\varepsilon} \frac{\partial \mathbf{c}}{\partial \mathbf{x}\_{N}} \boldsymbol{n}\_{i,j} + \dot{\boldsymbol{w}}\_{\varepsilon} \tag{3}$$

where *uj* is the *<sup>j</sup>*th component of the flow velocity, *<sup>ρ</sup>* is the fluid density, *Dc* is the diffusivity of *<sup>c</sup>* and . *wc* is chemical reaction rate. The reaction rate . *wc* variation with *c* for the present thermo-chemistry is shown in Figure 1 of ref. [59] and thus will not be repeated in this paper. The flame normal vector, *<sup>n</sup>*, of a *<sup>c</sup>* iso-surface is defined as:*ni* <sup>=</sup> <sup>−</sup>(*∂c*/*∂xi*)/|∇*c*<sup>|</sup> <sup>=</sup> <sup>−</sup>*c*,*i*/|∇*c*<sup>|</sup> <sup>=</sup> *<sup>c</sup>*,*i*/(*∂c*/*∂xN*), where *xN* is the coordinate in the normal direction to the iso-surface. The quantity, *ni*,*i*/2 = 0.5*∂ni*/*∂xi* = *κ<sup>m</sup>* is the mean value of two principal curvatures of the iso-surface and will henceforth be referred to as the flame curvature in this paper. According to the convention used here, the flame normal points towards the reactants and the flame surface has a positive (negative) curvature where it is convex (concave) to the reactants. The transport equation of *c*(*x*, *t*) can alternatively be presented in the kinematic form as [21,22,24,32–41]:

$$\frac{\partial c}{\partial t} + \mu\_{\dot{\jmath}} \frac{\partial c}{\partial x\_{\dot{\jmath}}} = S\_d |\nabla c| \tag{4}$$

where *Sd* is the displacement speed, which is expressed as [21,22,24,32–41]:

$$\mathcal{S}\_d(\mathbf{x}, t) = \frac{1}{\rho |\nabla c|} \frac{\partial}{\partial \mathbf{x}\_N} \left(\rho D\_\varepsilon \frac{\partial \mathbf{c}}{\partial \mathbf{x}\_N}\right) - 2D\_\varepsilon \mathbf{x}\_m + \frac{\dot{w}\_c}{|\nabla c|} \tag{5}$$

**Figure 1.** Instantaneous isosurfaces of reaction progress variable *c* = 0.8 colored by non-dimensional temperature *θ* (**a**), local values of normalised curvature *κ<sup>m</sup>* × *δth* (**b**) and normalised Gauss curvature *<sup>κ</sup><sup>g</sup>* <sup>×</sup> *<sup>δ</sup>*<sup>2</sup> *th* (**c**) for cases with *Le* = 0.8, 1.0 and 1.2 (1st–3rd row).

Differentiating Equation (3) with respect to *xi* yields [52,53]:

$$\frac{1}{|\nabla c|} \left( \frac{\partial |\nabla c|}{\partial t} + v^{\epsilon}\_{\not{j}} \frac{\partial |\nabla c|}{\partial x\_{\not{j}}} \right) = -n\_{i} v^{\epsilon}\_{\not{j},i} n\_{j} \tag{6}$$

$$\frac{\partial n\_i}{\partial t} + v^c\_{\dot{j}} \frac{\partial n\_i}{\partial x\_{\dot{j}}} = - \left(\delta\_{i\dot{j}} - n\_i n\_{\dot{j}}\right) S^c\_{\dot{j}k} n\_k + \mathcal{W}^c\_{\dot{i}\dot{j}} n\_{\dot{j}};\tag{7}$$

where *v<sup>c</sup> <sup>j</sup>* = *uj* + *Sdnj* is the *j*th component of local propagation velocity of a given *c* isosurface. The quantity *niv<sup>c</sup> j*,*i nj* = *niS<sup>c</sup> ijnj* = *<sup>a</sup><sup>c</sup> <sup>N</sup>* = *aN* + *∂Sd*/*∂xN* is often referred to as the total or effective normal strain rate [46,48,49,53,56,57], whereas *aN* = *niSijnj* is the fluid-dynamic normal strain rate and *∂Sd*/*∂xN* is an added normal strain rate induced by flame propagation [46,48,49,53,60,61]. The gradient of the flame propagation velocity *v<sup>c</sup> <sup>j</sup>* takes the following form [52,53]:

$$\begin{split} w^{\epsilon}\_{j,i} = \frac{\partial v^{\epsilon}\_{j}}{\partial x\_{i}} &= \frac{1}{2} \left( \frac{\partial v^{\epsilon}\_{j}}{\partial \overline{x\_{i}}} + \frac{\partial v^{\epsilon}\_{i}}{\partial \overline{x\_{j}}} \right) + \frac{1}{2} \left( \frac{\partial v^{\epsilon}\_{j}}{\partial \overline{x\_{i}}} - \frac{\partial v^{\epsilon}\_{i}}{\partial \overline{x\_{j}}} \right) \\ &= \frac{1}{2} \left( \frac{\partial u\_{j}}{\partial \overline{x\_{i}}} + \frac{\partial u\_{i}}{\partial \overline{x\_{j}}} \right) + \frac{1}{2} \left( \frac{\partial u\_{j}}{\partial \overline{x\_{i}}} - \frac{\partial u\_{i}}{\partial \overline{x\_{j}}} \right) + \frac{1}{2} \left( \frac{\partial S\_{d}}{\partial \overline{x\_{i}}} n\_{j} + \frac{\partial S\_{d}}{\partial \overline{x\_{j}}} n\_{i} \right) \\ &+ \frac{1}{2} \left( \frac{\partial S\_{d}}{\partial \overline{x\_{i}}} n\_{j} - \frac{\partial S\_{d}}{\partial \overline{x\_{j}}} n\_{i} \right) + S\_{d} \frac{1}{2} \left( \frac{\partial n\_{i}}{\partial \overline{x\_{i}}} + \frac{\partial n\_{i}}{\partial \overline{x\_{j}}} \right) + S\_{d} \frac{1}{2} \left( \frac{\partial n\_{i}}{\partial \overline{x\_{i}}} - \frac{\partial n\_{i}}{\partial \overline{x\_{j}}} \right) \\ &= S\_{ij}^{\epsilon} - \mathcal{W}\_{ij}^{\epsilon} \end{split} \tag{8}$$

Table 1 presents the nomenclature associated with Equation (8), where the total strain rate and rotation rate tensors are expressed as *S<sup>c</sup> ij* = *Sij* + *<sup>S</sup><sup>a</sup> ij* and *<sup>W</sup><sup>c</sup> ij* = *Wij* + *<sup>W</sup><sup>a</sup> ij* respectively [52,53]. The strain rate *Sij* and rotation rate *Wij* tensors originate due to the fluid motion, and the additional strain rate *S<sup>a</sup> ij* and rotation *W<sup>a</sup> ij* rate tensors originate due to the gradients of spatially dependent *Sd*.


**Table 1.** Nomenclature associated with the velocity gradient tensor *v<sup>c</sup> <sup>i</sup>*,*<sup>j</sup>* of an iso-surface element.

Taking the derivative with respect to *xi* on both sides of Equation (7) yields a transport equation of *κ<sup>m</sup>* = 0.5(*κ*<sup>1</sup> + *κ*2) = 0.5*ni*,*<sup>i</sup>* (with *κ*<sup>1</sup> and *κ*<sup>2</sup> being the two principal curvatures) [52,53]:

$$\begin{split} \underbrace{\frac{\partial \mathbf{x}\_{n}}{\partial t} + \boldsymbol{\upsilon}\_{j}^{c}}\_{}^{\text{dir}} \underbrace{\frac{\partial \mathbf{x}\_{n}}{\partial \mathbf{x}\_{j}}}\_{} &= \underbrace{\underbrace{\frac{\partial \mathbf{x}\_{N}(n\_{i,i})}{\mathbf{T}\_{1}}}\_{T\_{1}} + \underbrace{\frac{1}{2} \frac{\partial \mathbf{a}\_{N}}{\partial \mathbf{x}\_{N}} \underbrace{-S\_{ij}n\_{j,i}}\_{T\_{3}} - \underbrace{\frac{1}{2} \frac{\partial S\_{ij}}{\partial \mathbf{x}\_{i}}n\_{j}}\_{T\_{4}} + \underbrace{\frac{1}{2} \frac{\partial \mathbf{W}\_{ij}}{\partial \mathbf{x}\_{i}}n\_{j}}\_{T\_{5}} \\ &+ \underbrace{\frac{1}{2} \frac{\partial S\_{d}}{\partial \mathbf{x}\_{N}}n\_{i,i}}\_{T\_{6}} + \underbrace{\frac{1}{2} \underbrace{\frac{\partial^{2} S\_{d}}{\partial \mathbf{x}\_{N}} - S\_{ij}^{d}n\_{i,i}}\_{T\_{7}} - \underbrace{\frac{1}{2} \frac{\partial S\_{ij}^{d}}{\partial \mathbf{x}\_{i}}n\_{j}}\_{T\_{9}} + \underbrace{\frac{1}{2} \frac{\partial \mathbf{W}\_{ij}^{d}}{\partial \mathbf{x}\_{i}}n\_{j}}\_{T\_{10}} \end{split} \tag{9}$$

The terms *T*1−<sup>5</sup> arise due to fluid motion, whereas *T*6−<sup>10</sup> originate due to flame propagation. The positive (negative) contributions of these terms tend to increase the convexity (concavity) of iso-surfaces. The physical significances of the terms on the right side of Equation (9) are summarised in Table 2.

**Table 2.** Description of the various terms in the curvature transport Equation (9).


It can be appreciated from Equation (9) that curvature transport depends mainly on the statistics of fluid velocity/vorticity, scalar gradient and displacement speed. It has been demonstrated in the past that displacement speed statistics from simple chemistry [38–41,45] and detailed chemistry [32–37] DNS are qualitatively similar. The same is true for the statistics of the reactive scalar gradient obtained from simple chemistry [43–45,62] and detailed chemistry [47,62] DNS studies. Moreover, the vorticity and sub-grid flux statistics obtained from simple chemistry [59,63,64] DNS are found to be qualitatively consistent with those obtained from detailed chemistry [65,66] DNS. Furthermore, several models developed based on simple chemistry data [64,67,68] have been found to perform equally well in the context of detailed chemistry and transport [66,69,70].

#### **3. Numerical Implementation**

In the present analysis, DNS simulations of spherically expanding turbulent premixed flames with *Le* = 0.8, 1.0 and 1.2 have been carried out using a well-known compressible code SENGA [71] where the conservation equations of mass, momentum, energy and reaction progress variable have been solved in non-dimensional form. The *Le* = 0.8 case is representative of hydrogen-blended methane-air mixtures (e.g., 10% by volume hydrogen blended methane-air flames with overall equivalence ratio of 0.6) and the Lewis number 1.2 case is representative of a hydrocarbon-air mixture involving a hydrocarbon fuel which is heavier than methane (e.g., ethylene-air mixture with equivalence ratio of 0.7) [72,73]. The unity Lewis number flames are analogous to the stoichiometric methane-air flame [72,73]. The spatial discretisation for internal grid points has been carried out using a 10th order central difference scheme and the order of differentiation decreases gradually to a one-sided 2nd order scheme at the non-periodic boundaries [71]. The time-advancement has been carried out using a 3rd order explicit Runge–Kutta scheme [74]. The computational domain is taken to be a cube of 58.10*δth* × 58.10*δth* × 58.10*δth* (where *δth* = (*Tad* − *T*0)/max|∇*T*|*<sup>L</sup>* is the thermal flame thickness and the subscript 'L' is used to refer to conditions in the unstrained laminar premixed flame), which is discretised by a uniform Cartesian grid of 650 × 650 × 650. In all the cases considered here, the boundaries of the computational domain are taken to be partially non-reflecting and are specified using the Navier–Stokes Characteristic Boundary conditions (NSCBC) technique [75]. The reacting scalar fields obtained from a steady state unstrained laminar flame have been utilised to create a burned gas sphere with its centre initially at the centre of the domain. This reacting scalar field is allowed to evolve in a quiescent environment at least for one chemical time scale (i.e., *t* = *δth*/*SL*). For all simulations, standard values are taken for Prandtl number *Pr* and Zel'dovich number *<sup>β</sup>* = *Tac*(*Tad* − *<sup>T</sup>*0)/*T*<sup>2</sup> *ad* (i.e., *Pr* = 0.7 and *β* = 6.0 with *Tac* being the activation temperature) and the heat release parameter *τ* = (*Tad* − *T*0)/*T*<sup>0</sup> is taken to be 4.5 for all cases considered here. The spherical laminar flame kernels for different Lewis numbers with a normalised radius of *rSL*/*αT*<sup>0</sup> = 10.6 (where *αT*<sup>0</sup> is the thermal diffusivity of the unburned gas) based on the region corresponding to *c* ≥ 0.85 have been considered as the initial condition for turbulent simulations. A standard pseudo-spectral method [76] has been adopted to initialise homogeneous isotropic turbulent velocity fluctuations following the model spectrum by Pope [77]. For all cases the initial values of the normalised root-mean-square velocity fluctuation and longitudinal integral length scale are given by: *u* /*SL* = 7.5 and *l*/*δth* = 4.58, respectively. These values of *u* /*SL* and *l*/*δth* are representative of the thin reaction zones regime of combustion according to the regime diagram by Peters [78]. All the turbulent simulations have been continued for 2 initial eddy turnover times (i.e., *t* = *l*/ <sup>√</sup>*<sup>k</sup>* where *<sup>k</sup>* is the turbulent kinetic energy based on the whole domain), which is equivalent to 1.0 chemical timescale (i.e., *tchem* = *δth*/*SL*). By this time, the turbulent kinetic energy was not varying rapidly with time and *u* decayed by 40% in comparison to its initial value.

#### **4. Results and Discussion**

#### *4.1. Curvature Characterisation*

The flame topology is often characterised by the mean of principal curvatures (i.e., *κ<sup>m</sup>* = 0.5(*κ*<sup>1</sup> + *κ*2)) and Gauss curvature *κ<sup>g</sup>* = *κ*1*κ*<sup>2</sup> [79] where *κ*<sup>1</sup> and *κ*<sup>2</sup> are the principal curvatures. The *c* = 0.8 isosurface coloured by local values of non-dimensional temperature *θ*, normalised mean curvature *<sup>κ</sup><sup>m</sup>* × *<sup>δ</sup>th* and normalised Gauss curvature *<sup>κ</sup><sup>g</sup>* × *<sup>δ</sup>*<sup>2</sup> *th* for the *Le* = 0.8, 1.0 and 1.2 flames are shown in Figure 1. It is worth noting that the reaction rate . *wc* assumes maximum value close to *c* = 0.8 for the present thermochemistry and this isosurface can be taken to be the flame surface for the purpose of this analysis.

It is evident from Figure 1 that the extent of wrinkling and burned gas volume increase with decreasing *Le*, and especially the occurrences of saddle point topologies (i.e., *κ<sup>g</sup>* < 0) and sharply negatively curved cusps (i.e., large magnitudes of negative *κm*) increase with decreasing Lewis number. The augmentation of flame wrinkling with decreasing *Le* has implications on the extent of burning and flame surface area, which can be substantiated from Table 3 where the normalised volume-integrated burning rate Ω*T*/Ω*<sup>L</sup>* and flame surface area *AT*/*AL* are listed. Here, Ω and *A* are evaluated by volume-integrals Ω = *V ρ* . *wcdV* and A = *<sup>V</sup>*|∇*c*|*dV*, respectively and the subscripts T and L refer to the values in turbulent flame and initial laminar flame kernels, respectively. The fresh reactants diffuse faster into the reaction zone than the rate at which heat diffuses out in the *Le* = 0.8 flame. This gives rise to higher burning rate in the *Le* = 0.8 flame than in the corresponding unity Lewis number case. The diffusion of reactants is slower into the reaction zone than the rate at which heat diffuses out in the *Le* = 1.2 flame, which in turn reduces the burning rate in the *Le* = 1.2 case than in the corresponding unity Lewis number case. This is qualitatively consistent with several previous findings [1,39,41,42,44,54–58].

**Table 3.** Lewis number dependence of the normalised volume-integrated burning rate Ω*T*/Ω*<sup>L</sup>* and flame surface area *AT*/*AL* when the statistics were extracted.


It can be seen from Figure 1 that high temperature zones are associated with positive *κ<sup>m</sup>* values in the *Le* = 0.8 case, whereas the high temperature zones are associated with negatively curved zones (i.e., *κ<sup>m</sup>* < 0) in the *Le* = 1.2 case. The non-dimensional temperature *θ* remains uniform and equal to 0.8 (i.e., *θ* = *c* = 0.8) on the *c* = 0.8 isosurface in the unity Lewis number case. The statistically spherical flames have predominantly positive curvatures where the combination of strong focussing of reactants and weak defocussing of heat gives rise to high reaction rate in the positively curved regions in the *Le* = 0.8 flame. Just the opposite mechanisms lead to low reaction rate at the negatively curved regions in the *Le* = 0.8 flame. This also gives rise to *θ* > *c* at the positively curved locations in the *Le* = 0.8 flame. As a result of this, the *Le* = 0.8 case exhibits stable positively curved bulges with large radius of curvature separated by sharply negatively curved cusps. In the *Le* = 1.2 case, the focussing of heat is faster than the defocussing of reactants at the negatively curved locations, which increases the reaction rate magnitude in these locations and may give rise to local occurrences of *θ* > *c*. Thus, the sharply negatively curved pockets burn faster in the *Le* = 1.2 case and as a result sharp negatively curved cusps with small radius of curvature are unlikely to survive in this case. The unity Lewis number case is thermo-diffusively neutral and thus, for low-Mach number adiabatic conditions, the non-dimensional temperature *θ* remains equal to the reaction progress variable.

The high temperature values at the positively curved zones significantly increase the fuel consumption rate per unit area in the turbulent *Le* = 0.8 case in comparison to the corresponding unstrained laminar flame value and this gives rise to a Ω*T*/Ω*<sup>L</sup>* value, which is greater than *AT*/*AL*. The consumption rate per unit area in the turbulent *Le* = 1.0 case remains comparable to the corresponding unstrained laminar flame value and thus Ω*T*/Ω*<sup>L</sup>* and *AT*/*AL* remain close to each other. The combination of strong thermal diffusion and weak diffusion of reactants into the reaction zone reduces the consumption rate per unit area in the positively stretched zones in turbulent *Le* = 1.2 case in comparison to the corresponding unstrained laminar flame value, which results in Ω*T*/Ω*<sup>L</sup>* < *AT*/*AL* in this case.

The probability density functions (PDFs) of normalised curvature *κ<sup>m</sup>* × *δth* for different *c*-isosurfaces across the flame front are shown in Figure 2. It is evident from Figure 2 that the width of the curvature PDFs tends to increase with decreasing Lewis number due to increased flame wrinkling with decreasing Lewis number *Le*. The most probable value of *κ<sup>m</sup>* has been found to be close to zero for all flame kernels considered here but the mean value of *κ<sup>m</sup>* remains positive due to the statistically spherical configuration. In the *Le* ≥ 1 flames, the PDFs of normalised curvature *κ<sup>m</sup>* × *δth* remain

almost symmetric around *κ<sup>m</sup>* × *δth* = 0 towards the unburned gas side of the flame front (e.g., *c* = 0.1) due to the strong deformation of the flame surface as a result of the interaction between energetic eddies with the preheat zone. However, the PDFs of normalised curvature *κ<sup>m</sup>* × *δth* exhibit higher propensity to obtain positive values for the major part of flame front. It can be seen from Figure 2 that the high magnitudes of negative *κ<sup>m</sup>* are more frequent for the *c* = 0.5, 0.7 and 0.9 isosurfaces in the *Le* = 0.8 case than in the corresponding *Le* = 1.0 and 1.2 cases. This tendency for the *Le* = 0.8 flame could be a consequence of thermo-diffusive instability associated with *Le* < 1, where sharply negative curved cusps with small radii of curvature are found in between positively curved bulges with large or moderate radii of curvature.

**Figure 2.** Probability density functions (PDFs) of *κ<sup>m</sup>* × *δth* on (**a**–**e**) *c* = 0.1, 0.3, 0.5, 0.7 and 0.9 isosurfaces for *Le* = 0.8, 1.0 and 1.2.

The Joint Probability Density Functions (JPDFs) of normalised mean curvature *κ<sup>m</sup>* × *δth* and normalised Gauss curvature *<sup>κ</sup><sup>g</sup>* × *<sup>δ</sup>*<sup>2</sup> *th* for different *c* iso-surfaces across the flame are shown in Figure 3. The region corresponding to *κ<sup>g</sup>* > *κ*<sup>2</sup> *<sup>m</sup>* on the *κ<sup>m</sup>* − *κ<sup>g</sup>* plane corresponds to complex principal curvature values and thus cannot be realised in practice. The combination of *κ<sup>m</sup>* > 0 (*κ<sup>m</sup>* < 0) and *κ<sup>g</sup>* > 0 corresponds to cup convex (concave) shapes. By contrast, the combination of *κ<sup>m</sup>* > 0 (*κ<sup>m</sup>* < 0) and *κ<sup>g</sup>* < 0 corresponds to saddle convex (concave) shapes. Finally, the combination of *κ<sup>m</sup>* > 0 (*κ<sup>m</sup>* < 0) and *κ<sup>g</sup>* = 0 corresponds to tile convex (concave) flame topologies. For a schematic diagram of these topologies interested readers are referred to Figure 1 of ref. [53]. The spread of *κ<sup>m</sup>* × *δth* values on both positive and negative sides for the *Le* = 0.8 case is greater than that in the corresponding *Le* = 1.0 and 1.2 cases, which is indicative of larger extent of flame wrinkling in the *Le* = 0.8 case (see Figure 1 and Table 3). The highest values of joint PDFs are obtained for a positive value of *κ<sup>m</sup>* close to *κ<sup>m</sup>* ≈ 0, which is consistent with the observations made from Figure 2. Figure 3 demonstrates that the joint PDFs in the *Le* = 0.8 case show skewness for the negative values of *κ<sup>m</sup>* × *δth* on the unburned gas side of the flame, but the joint PDFs eventually become skewed towards the positive values of *κ<sup>m</sup>* as the burned gas side is approached. The distribution depicted by the joint PDF also shows a propensity to obtain predominantly positive *κ<sup>m</sup>* towards the burned gas side of the flame front for both *Le* = 1.0 and 1.2 cases but this tendency is more prevalent for the *Le* = 1.2 case. The high probability of finding negatively curved cusps (*κ<sup>m</sup>* < 0) for the *Le* = 0.8 case is a reflection of the thermo-diffusive instability in this flame. Both cup and saddle topologies appear considerably for convex and concave curvatures throughout the flame front in all cases. However, the cup like topology, whether for concave or convex structure, is more probable than the saddle and tile topologies for all cases considered here. The probability of finding saddle topologies decreases from the unburned to the burned gas side for all flames considered here. Moreover, the probability of finding saddle topologies decreases with increasing *Le*, which is consistent with the observations made from Figure 1.

**Figure 3.** Joint PDF (coloured from blue to dark red) between normalised curvature *κ<sup>m</sup>* × *δth* and normalised Gauss curvature *<sup>κ</sup><sup>g</sup>* <sup>×</sup> *<sup>δ</sup>*<sup>2</sup> *th* for *c* = 0.3, 0.5, 0.7 and 0.9 (1st–4th column) isosurfaces for cases *Le* = 0.8,1.0 and 1.2 (1st–3rd row).

The budgets of the various terms in the curvature transport equation conditional on *c* will be discussed in the next section. Based on this analysis it will be possible to judge the magnitude of the different terms and which terms possibly balance each other. Beside this, another focus of this paper is to establish a relation between the terms of the curvature transport equation and flame instabilities that might occur under certain conditions. It will therefore be instructive to analyse the terms *T*1, ... , *T*<sup>10</sup> conditional on mean curvature in the next but one section. Only by this additional representation, it will be possible to understand if positive or negative flame curvature elements are damped or possibly amplified in the presence of flame instabilities.

#### *4.2. Mean Profiles of Source/Sink Terms of the Curvature Transport Equation*

The profiles of mean values of *Ti* × *<sup>δ</sup>*<sup>2</sup> *th*/*SL* conditional on *c* for *Le* = 0.8, 1.0 and 1.2 cases are shown in Figure 4a–c respectively. The objective of Figure 4 is to provide insights into the mean behaviour of the terms on the right-hand side of the curvature transport equation (i.e., Equation (9)) across the flame. It can be seen from Figure <sup>4</sup> that the term *<sup>T</sup>*<sup>1</sup> × (*δ*<sup>2</sup> *th*/*SL*), originating from the correlation between curvature and normal strain rate, assumes weakly negative values for the major part of the flame front for the for *Le* = 1.0 and 1.2 cases. In these cases, the mean value of *<sup>T</sup>*<sup>1</sup> × (*δ*<sup>2</sup> *th*/*SL*) remains weakly positive towards the burned gas side of the flame front. However, the mean values of *<sup>T</sup>*<sup>1</sup> × (*δ*<sup>2</sup> *th*/*SL*) remain positive throughout the flame front for the *Le* = 0.8 case. Therefore, the term *T*<sup>1</sup> tends to increase the convexity of iso-surface in the *Le* = 0.8 case, while the opposite behaviour is obtained in the *Le* = 1.0 and 1.2 cases. In the *Le* = 1.0 and 1.2 flames the reaction progress variable gradient ∇*c* aligns predominantly with the most compressive principal strain rate due to stronger turbulent straining than the strain rate arising from flame normal acceleration induced by chemical heat release. This leads to predominantly negative normal strain rate *aN* = (*e<sup>α</sup>* cos2 *θα* + *e<sup>β</sup>* cos2 *θβ* + *e<sup>γ</sup>* cos2 *θγ*) (where *eα*,*e<sup>β</sup>* and *e<sup>γ</sup>* are the most extensive, intermediate and the most compressive eigenvalues of the strain rate tensor *Sij* = 0.5(*∂ui*/*∂xj* + *∂uj*/*∂xi*) and *θα*, *θβ* and *θγ* are the angles of their eigenvectors with ∇*c*, respectively) in the *Le* = 1.0 and 1.2 flames considered here. In these cases, ∇*c* exhibits some tendency for local preferential collinear alignment with *e<sup>α</sup>* only in the heat release zone. However, the alignment of ∇*c* with the most extensive principal strain rate weakens especially in the positively curved locations due to defocussing of heat [27,39,40,43,80–82], which gives rise to a situation where the highly negative *aN* values are associated with highly positively curved locations. This gives rise to negative mean value of *<sup>T</sup>*<sup>1</sup> × (*δ*<sup>2</sup> *th*/*SL*) in the *Le* = 1.0 and 1.2 cases considered here. In the case of *Le* = 0.8, the alignment of ∇*c* with the eigenvector associated with *e<sup>α</sup>* (*eγ*) is stronger (weaker) than the *Le* = 1.0 and 1.2 cases owing to stronger heat release effects. Furthermore, the alignment of ∇*c* with the eigenvector associated with *eα* increases in the positively curved regions due to high reaction rate and heat release associated with these locations. As a result, high positive values of *aN* are obtained at the positively curved locations in the *Le* = 0.8 case to yield positive mean values of *<sup>T</sup>*<sup>1</sup> × (*δ*<sup>2</sup> *th*/*SL*).

The mean contribution to the curvature transport due to normal gradients of the flow normal strain rate, *<sup>T</sup>*<sup>2</sup> × (*δ*<sup>2</sup> *th*/*SL*) = [(*∂aN*/*∂xN*)/2]/(*δ*<sup>2</sup> *th*/*SL*), assumes negative values on the reactant side but becomes positive towards the burned gas side of the flame. The normal strain rate *aN* increases gradually from the reactant side toward the product side due to the heat release and reaches its maximum close to the reaction zone, and hence its normal gradient *∂aN*/*∂xN* becomes zero, then it decreases again near the product side [61]. Accordingly, the tendency of (*∂aN*/*∂xN*) < 0 is high on the reactant side, while the opposite occurs on the product side of the flame front. Consequently, *T*<sup>2</sup> tends to act as a sink (source) term toward the reactant (product) side of the flame front. However, the profile and the magnitude of *aN* change due to differential diffusion effects induced by non-unity Lewis number [61]. This also alters the location at which the highest value of *aN* is obtained [61]. As a result of this, the magnitude of *T*<sup>2</sup> increases with decreasing *Le* and in all cases the magnitude of the positive mean value of *<sup>T</sup>*<sup>2</sup> × (*δ*<sup>2</sup> *th*/*SL*) remains greater than the negative mean contribution towards the unburned gas side.

**Figure 4.** Profiles of the normalised mean values of the fluid flow induced terms (i.e., *<sup>T</sup>*1−<sup>5</sup> <sup>×</sup> *<sup>δ</sup>*<sup>2</sup> *th*/*SL*) (left column) and flame propagation induced terms (i.e., *<sup>T</sup>*6−<sup>10</sup> <sup>×</sup> *<sup>δ</sup>*<sup>2</sup> *th*/*SL* ) (right column) of the curvature transport equation conditional upon *c* for cases (**a**–**c**) *Le* = 0.8, 1.0 and 1.2.

The mean value of the curvature flow stretching term *<sup>T</sup>*<sup>3</sup> × (*δ*<sup>2</sup> *th*/*SL*) = −(*Sijnj*,*i*)(*δ*<sup>2</sup> *th*/*SL*) is found to be positive throughout the flame brush in all cases considered here. The magnitude of *T*<sup>3</sup> becomes comparable to *T*<sup>2</sup> for the *Le* = 1.2 case but its magnitude in comparison to *T*<sup>2</sup> remains small in the *Le* = 0.8 case. Using (*xT*1, *xT*2) as the local principal axes along the tangential directions on a given *c* isosurface, the curvature flow stretching term can be expressed as *T*<sup>3</sup> = −(*S*11*κ*<sup>1</sup> + *S*22*κ*2), where *S*<sup>11</sup> and *S*<sup>22</sup> are the tangential strain rates along axes *xT*<sup>1</sup> and *xT*2, respectively. The mean tangential strain rate *aT* = *S*<sup>11</sup> + *S*<sup>22</sup> assumes positive values throughout the flame front for all flames considered here [43–49]. A positive magnitude of *T*<sup>3</sup> is obtained for a combination of cup and saddle concave iso-scalar topologies. For *S*<sup>11</sup> > 0 and *S*<sup>22</sup> > 0, one can obtain *T*<sup>3</sup> > 0 for a cup concave structure where *κ*<sup>1</sup> < 0 and *κ*<sup>2</sup> < 0 or for a saddle topology where either *κ*<sup>1</sup> < 0 and *κ*<sup>2</sup> > 0 when *S*<sup>11</sup> > *S*<sup>22</sup> or *κ*<sup>1</sup> > 0 and *κ*<sup>2</sup> < 0 when *S*<sup>11</sup> < *S*22.

The statistical behaviour of the mean value of the term *<sup>T</sup>*<sup>4</sup> × (*δ*<sup>2</sup> *th*/*SL*) = −{[(*∂Sij*/*∂xi*)*nj*]/2}(*δ*<sup>2</sup> *th*/*SL*) conditional on *c* shows an opposite trend in comparison to that of *T*2, in other words it exhibits positive contribution toward the reactant side but becomes negative on the burned gas side. This means that the two vectors with components *∂Sij*/*∂xi* and *nj*, respectively are in the same direction on the reactant side, while they point in the opposite direction on the burned gas side. The magnitudes of *T*<sup>2</sup> and *T*<sup>4</sup> are mostly comparable and remain in balance for all cases considered here.

The normalised contribution to the curvature transport due to vorticity gradients, *<sup>T</sup>*<sup>5</sup> × (*δ*<sup>2</sup> *th*/*SL*) = {[(*∂Wij*/*∂xi*)*nj*]/2}(*δ*<sup>2</sup> *th*/*SL*), can alternatively be written as: *<sup>T</sup>*<sup>5</sup> = {[*niεijk*(*∂ωk*/*∂xj*)]/4}(*δ*<sup>2</sup> *th*/*SL*). The mean value of *T*<sup>5</sup> can be expressed as *<sup>T</sup>*<sup>5</sup> = (*∂ω*2/*∂x*<sup>1</sup> <sup>−</sup> *∂ω*1/*∂x*2)/4, using the local principal axes (*xT*1, *xT*2, *xN*) and *<sup>n</sup>* = (0, 0, 1), where *ω*<sup>1</sup> and *ω*<sup>2</sup> are the components of the flow vorticity tangent to the *c* iso-surface. Thus, co-rotating parallel vortices of different intensity and counter-rotating parallel vortices of the same intensity can potentially curve a planar local iso-surface structure, leading to positive or negative curvatures. For the cases considered here, the mean value of *T*<sup>5</sup> shows negligible contribution throughout the flame front.

In the *Le* = 1.2 case, the mean values of the term *<sup>T</sup>*<sup>6</sup> × (*δ*<sup>2</sup> *th*/*SL*) = {(*∂Sd*/*∂xN*)*ni*,*i*/2}(*δ*<sup>2</sup> *th*/*SL*) which arises due to the correlation between (*∂Sd*/*∂xN*) and *κm*, exhibit positive contributions throughout the flame front, while a weakly positive trend is seen in the *Le* = 1.0 case for the major part of the flame. However, this term shows negative contribution throughout the flame in the *Le* = 0.8 case. It has been demonstrated elsewhere that the contributions of (*∂Sd*/*∂xN*) and *aN* are of same order of magnitude but opposite in behaviour [61], and hence an opposite behaviour for the mean values of *T*<sup>1</sup> and *T*<sup>6</sup> is observed. It is also worth noting that (*∂Sd*/*∂xN*) is predominantly negative and its magnitude decreases with increasing *Le* [61]. Furthermore, small values of the flame thickness are associated with the positively curved locations in the *Le* = 0.8 case [60], which tends to increase the magnitude of negative (*∂Sd*/*∂xN*) at positive *κ<sup>m</sup>* values. This leads to predominantly negative mean contribution of *T*<sup>6</sup> in the *Le* = 0.8 flame. By contrast, small values of the flame thickness are associated with the negatively curved locations in the *Le* = 1.2 case [61], and this increases the magnitude of negative (*∂Sd*/*∂xN*) at negative *κ<sup>m</sup>* values, leading to positive mean values of *T*6. The quantity (*∂Sd*/*∂xN*) remains weakly correlated with *κ<sup>m</sup>* in the *Le* = 1.0 case [61], which leads to weak positive mean value of this term throughout the flame front.

The normalised mean value of the term *<sup>T</sup>*<sup>7</sup> × (*δ*<sup>2</sup> *th*/*SL*) = [(*∂*<sup>2</sup>*Sd*/*∂xN*<sup>2</sup>)/2](*δ*<sup>2</sup> *th*/*SL*), due to the normal gradient of the added normal strain rate, assumes a positive mean value towards the unburned gas side but becomes negative toward the burned gas side of the flame. The magnitude of *T*<sup>7</sup> decreases with increasing *Le*. As it has been shown elsewhere [61] that the contributions of (*∂Sd*/*∂xN*) and *aN* are of the same order of magnitude but opposite in sign, their normal gradients (*∂*<sup>2</sup>*Sd*/*∂xN*<sup>2</sup> and *∂aN*/*∂xN*) are expected to exhibit the opposite behaviour, and thus, *T*<sup>7</sup> exhibits an opposite trend to *T*2.

In the thin reaction zones regime flames considered here, the normalised added stretching term *<sup>T</sup>*<sup>8</sup> × (*δ*<sup>2</sup> *th*/*SL*) = −(*S<sup>a</sup> ijnj*,*i*)(*δ*<sup>2</sup> *th*/*SL*) acts as a leading order contributor in all cases where it remains negative throughout the flame. Using the local principal axes (*xT*1, *xT*2, *xN*) and for *n* = (0, 0, 1), one can express the term as *<sup>T</sup>*<sup>8</sup> = −*Sd*(*κ*<sup>2</sup> <sup>1</sup> + *<sup>κ</sup>*<sup>2</sup> <sup>2</sup>), which explains the predominant negative contribution of *T*<sup>8</sup> throughout the flame front. The magnitude of *T*<sup>8</sup> increases with decreasing *Le* due to a larger spread of displacement speed *Sd* and principal curvatures *κ*<sup>1</sup> and *κ*<sup>2</sup> as a result of greater extent of flame wrinkling.

The second derivatives of *Sd* and *ni* determine the normalised mean contributions due to the added strain rate gradients *<sup>T</sup>*<sup>9</sup> × (*δ*<sup>2</sup> *th*/*SL*) = −{[(*∂S<sup>a</sup> ij*/*∂xi*)*nj*]/2}(*δ*<sup>2</sup> *th*/*SL*), which exhibit negative (positive) contribution towards reactant (product) side in all cases. The mean contribution of the added vorticity contribution term *<sup>T</sup>*<sup>10</sup> = [−0.25(*δij* − *ninj*)*∂*<sup>2</sup>*Sd*/*∂xi∂xj* + 0.25(*∂Sd*/*∂xN*)*ni*,*i*] remains negligible in all cases considered here.

It can be inferred from Figure 4 that the terms *T*2, *T*4, *T*7, *T*<sup>8</sup> and *T*<sup>9</sup> are the leading order terms of the mean curvature transport for all cases. In addition, the term *T*<sup>3</sup> plays a leading order role for the *Le* = 1.2 case. However, in all cases considered here, the mean contributions of *T*<sup>2</sup> and *T*<sup>4</sup> assume comparable magnitudes with opposite signs and remain mostly in balance.

#### *4.3. Mean Profiles of the Terms of the Curvature Transport Equation Conditioned Upon Curvature*

The normalised mean values of flow-induced terms (i.e., *T*1−5) and added flame propagation terms (i.e., *T*6−10) in the mean curvature transport equation conditional upon the normalised curvature *κmδth* for the *c* = 0.8 isosurface are shown in Figure 5. It can be seen from Figure 5 that the magnitudes of the mean values of *T*1−<sup>10</sup> increase with decreasing *Le*. In Figure 5, the positive contributions of the terms in the positively curved locations tend to increase the convexity of the flame surface, whereas the negative contributions of the terms act to reduce the convexity at *κ<sup>m</sup>* > 0. By the same token, the positive contributions of the terms in the negatively curved locations tend to decrease the concavity of the flame surface, whereas the negative contributions of the terms act to increase the concavity at *κ<sup>m</sup>* < 0.

It is evident from Figure 5 that the mean contribution of *T*<sup>1</sup> = [(*aNni*,*i*)/2] assumes negative values for *κ<sup>m</sup>* < 0 and positive values for *κ<sup>m</sup>* > 0 for all cases on the *c* = 0.8 isosurface. In the reactive region *aN* assumes predominantly positive values due to predominant preferential alignment between ∇*c* with the eigenvector associated with *eα*, and this leads to negative (positive) values of *T*<sup>1</sup> = [(*aNni*,*i*)/2] at the negatively (positively) curved zones on the flame surface.

In all cases, the term *T*<sup>2</sup> exhibits negative mean values associated with *κ<sup>m</sup>* > 0 and positive mean values are obtained for *κ<sup>m</sup>* < 0. The flow divergence ahead of the positively curved zones and flow convergence ahead of the negative curved regions lead to positive mean values of *T*<sup>2</sup> for *κ<sup>m</sup>* < 0 and negative values of *T*<sup>2</sup> for *κ<sup>m</sup>* > 0. The mean contribution of the flow stretching term *T*<sup>3</sup> = −(*Sijnj*,*i*) shows large positive values for convex topologies and weakly positive mean values for topologies with concave curvatures in the *Le* = 1.0 and 1.2 cases, whereas large positive mean values of *T*<sup>3</sup> are obtained at both highly positive and negative curved locations for the *Le* = 0.8 case. This is in agreement with Figure 4, which shows that the mean value of *T*<sup>3</sup> conditional on *c* remains positive throughout the flame. The probability of finding both cup concave and saddle concave topologies is greater in the *Le* = 0.8 case than in the *Le* = 1.0 and 1.2 cases, where *T*<sup>3</sup> can assume positive values, and thus the mean value of *T*<sup>3</sup> conditional on *κ<sup>m</sup>* assumes large positive values in the negative curved zones. The mean value of *T*<sup>4</sup> conditional on *κ<sup>m</sup>* shows an opposite behaviour in comparison to *T*2, and thus it exhibits positive (negative) mean values for *κ<sup>m</sup>* > 0 (*κ<sup>m</sup>* < 0). The mean value of *T*<sup>5</sup> conditioned on curvature shows almost similar behavior to *T*<sup>4</sup> but the conditional mean value of *T*<sup>5</sup> remains negligible for *κ<sup>m</sup>* < 0.

**Figure 5.** Profiles of the normalised mean values of the fluid flow induced terms (i.e., *<sup>T</sup>*1−<sup>5</sup> <sup>×</sup> *<sup>δ</sup>*<sup>2</sup> *th*/*SL*) (left column) and flame propagation induced terms (i.e., *<sup>T</sup>*6−<sup>10</sup> <sup>×</sup> *<sup>δ</sup>*<sup>2</sup> *th*/*SL* ) (right column) of the curvature transport equation conditional upon normalised curvature *κ<sup>m</sup>* × *δth* on the *c* = 0.8 isosurface for cases (**a**–**c**) *Le* = 0.8, 1.0 and 1.2.

The mean contributions of the additional normal strain term *T*<sup>6</sup> = {(*∂Sd*/*∂xN*)*ni*,*i*/2} and the term due to the normal gradient of the added normal strain rate *T*<sup>7</sup> conditional upon curvature remain negligible in comparison to the conditional mean values of *T*8, *T*<sup>9</sup> and *T*10. The mean contribution to the added stretching term, *<sup>T</sup>*<sup>8</sup> = −(*S<sup>a</sup> ijnj*,*i*) conditioned on curvature plays a leading order role. It has already been mentioned that in the coordinate aligned with principal axes of curvature, *T*<sup>8</sup> can be expressed as: *<sup>T</sup>*<sup>8</sup> = −*Sd*(*κ*<sup>2</sup> <sup>1</sup> + *<sup>κ</sup>*<sup>2</sup> 2), which implies a non-linear (e.g., cubic) curvature dependence

of *T*<sup>8</sup> due to the curvature dependence of displacement speed *Sd* [32–41]. Consequently in all cases considered here, *T*<sup>8</sup> exhibits an asymmetric trend with respect to *κ<sup>m</sup>* = 0, where its negative (positive) mean values correlating to *κ<sup>m</sup>* < 0 (*κ<sup>m</sup>* > 0), with *κ<sup>m</sup>* = 0 being the inflection point for the *Le* = 0.8, 1.0 and 1.2 cases. However, in the *Le* = 0.8 case, the negative mean contribution of *T*<sup>8</sup> conditioned on *κ<sup>m</sup>* for negatively curved regions remains much greater than the positive mean contribution in the positively curved zones. By contrast, the positive mean contribution of *T*<sup>8</sup> conditioned on *κ<sup>m</sup>* for positively curved regions remains much greater than the negative mean contribution in the negatively curved zones in the *Le* = 1.0 and 1.2 cases. The mean contribution of the term *T*<sup>9</sup> due to the added strain rate gradients conditional upon curvature *κ<sup>m</sup>* exhibits large positive values for topologies with highly negative curvatures but its mean value assumes small negative values in the positively curved zones for the *Le* = 1.0 and 1.2 cases. In the *Le* = 0.8 case, the negative mean contribution of *T*<sup>9</sup> remains negligible at highly positively curved locations. The mean values of the added vorticity curl contribution *T*<sup>10</sup> conditional upon *κ<sup>m</sup>* exhibit similar qualitative trend as that of *T*9.

#### *4.4. Overall Behaviour of the Terms in the Curvature Transport Equation*

Figure 6 (left column) shows the normalised mean values of net fluid flow contributions (*T*<sup>1</sup> + ··· + *<sup>T</sup>*5) × (*δ*<sup>2</sup> *th*/*SL*), flame propagation induced added contributions (*T*<sup>6</sup> + ··· + *<sup>T</sup>*10) × (*δ*<sup>2</sup> *th*/*SL*) and the total contribution (*T*<sup>1</sup> + ··· + *<sup>T</sup>*10) × (*δ*<sup>2</sup> *th*/*SL*) conditional upon the reaction progress variable *c*. As *T*<sup>1</sup> and *T*<sup>5</sup> have negligible contributions and *T*<sup>4</sup> nullifies *T*2, the net mean values of fluid motion terms (*T*<sup>1</sup> + ··· + *<sup>T</sup>*5) × (*δ*<sup>2</sup> *th*/*SL*) follow the statistical trend of the curvature flow stretching term *<sup>T</sup>*<sup>3</sup> × (*δ*<sup>2</sup> *th*/*SL*) in all cases considered here. It is evident from Figure 6 (left column) that the net contribution of the added terms due to flame propagation (*T*<sup>6</sup> + ··· + *<sup>T</sup>*10) × (*δ*<sup>2</sup> *th*/*SL*) dominates over the net contribution of the terms (*T*<sup>1</sup> + ··· + *<sup>T</sup>*5) × (*δ*<sup>2</sup> *th*/*SL*) arising from fluid flow. The mean contributions of (*T*<sup>6</sup> + ··· + *<sup>T</sup>*10) × (*δ*<sup>2</sup> *th*/*SL*) and (*T*<sup>1</sup> + ··· + *<sup>T</sup>*10) × (*δ*<sup>2</sup> *th*/*SL*) remain negative throughout the flame front for all cases. This behaviour originates because the negative contributions of *T*<sup>8</sup> and *T*<sup>9</sup> dominate over the positive contribution arising from *T*<sup>7</sup> on the unburned gas side, whereas the combined negative contributions of *T*<sup>7</sup> and *T*<sup>8</sup> dominate over the positive contribution of *T*<sup>9</sup> on the product side of the flame front.

The statistical behaviours of the normalised mean values of flow contributions (*T*<sup>1</sup> + ··· + *T*5) × (*δ*<sup>2</sup> *th*/*SL*), flame propagation induced contributions (*T*<sup>6</sup> + ··· + *<sup>T</sup>*10) × (*δ*<sup>2</sup> *th*/*SL*) and the total contribution (*T*<sup>1</sup> + ··· + *<sup>T</sup>*10) × (*δ*<sup>2</sup> *th*/*SL*) to the curvature transport, conditioned on curvature *κ<sup>m</sup>* are shown in Figure 6 (right column) for the *c* = 0.8 isosurface. It is clear from Figure 6 (right column) that the mean values of the flow contributions (*T*<sup>1</sup> + ··· + *T*5) are negligible in comparison to the mean values of the overall added flame propagation terms (*T*<sup>6</sup> + ··· + *T*10) in the negative curved locations, whereas the flow terms (*T*<sup>1</sup> + ··· + *T*5) dominate over (*T*<sup>6</sup> + ··· + *T*10) for positive curvatures in the *Le* = 0.8 case. However, in the *Le* = 1.0 and 1.2 cases, the net mean flame propagation contribution (*T*<sup>6</sup> + ··· + *T*10) dominates over the mean contribution of the flow terms (*T*<sup>1</sup> + ··· + *T*5) especially for the highly positive and negative curvatures. It can clearly be seen that the net contribution of the added terms induced by flame propagation conditional upon mean curvature *κ<sup>m</sup>* shows the same trend of the added strain rate and vorticity gradients terms *T*<sup>9</sup> and *T*<sup>10</sup> in the positively curved locations, whereas the added curvature stretching term *T*<sup>8</sup> determines the net mean behaviour of the added flame propagation contribution in the negatively curved zones. Thus, the net mean contribution of (*T*<sup>6</sup> + ··· + *T*10) remains negative for both highly negative and positive curvatures. However, the magnitudes of mean negative (*T*<sup>6</sup> + ··· + *T*10) in the positively curved zones are much greater than that in the negative curvatures in the *Le* = 1.0 and 1.2 cases. Just the opposite behaviour can be seen for the *Le* = 0.8 case, and the magnitude of the negative mean contribution of (*T*<sup>6</sup> + ··· + *T*10) at the negatively curved locations is especially large in this case due to large negative values of *T*8. It is recalled that positivity (negativity) of (*T*<sup>1</sup> + ··· + *T*10) leads to increasing positivity (negativity) of mean curvature. As a result of such an amplification of positive (negative) curvature values, the position of flame elements in the diagrams shown in the right column of Figure 6 moves towards the right

(left) and the opposite happens for a damping of curvature values. In other words, Figure 6 implies a characteristic motion of the location of flame elements in this diagram until they either behave in a neutral manner or small-scale wrinkles with characteristic length scales considerably smaller than the inner cut-off scale are smoothed out by molecular diffusion effects.

**Figure 6.** Profiles of normalised mean values of the flow contributions (*T*<sup>1</sup> <sup>+</sup> ··· <sup>+</sup> *<sup>T</sup>*5) <sup>×</sup> (*δ*<sup>2</sup> *th*/*SL*), flame propagation induced contributions (*T*<sup>6</sup> <sup>+</sup> ··· <sup>+</sup> *<sup>T</sup>*10) <sup>×</sup> (*δ*<sup>2</sup> *th*/*SL*) and the total contribution (*T*<sup>1</sup> <sup>+</sup> ··· <sup>+</sup> *<sup>T</sup>*10) <sup>×</sup> (*δ*<sup>2</sup> *th*/*SL*) to the curvature transport, conditioned on the reaction progress variable *c* (left column) and normalised curvature *κ<sup>m</sup>* × *δth* on the *c* = 0.8 isosurface (right column) for cases (**a**–**c**) *Le* = 0.8, 1.0 and 1.2.

Figure 6 (right column) shows that the mean contribution of (*T*<sup>1</sup> + ··· + *T*10) remains strongly negative at negative values of *κ<sup>m</sup>* and mildly positive at positive values of *κ<sup>m</sup>* in the *Le* = 0.8 case. However, in the *Le* = 1.0 and 1.2 cases, the mean contribution of (*T*<sup>1</sup> + ··· + *T*10) assumes negative values at highly positive and negative curved locations and positive values are obtained at moderately positive curved locations. The above findings indicate that the concavity of the negatively curved cusps, and the convexity of the positively curved bulges are promoted by the curvature transport in the *Le* = 0.8 case, which is further aided by the presence of high and low temperature (and thus burning rate) regions at positively and negatively curved regions, respectively. As a result, positively curved bulges with large radii of curvature and intermediate negatively curved cusps with small radii of curvature are likely to be observed in the *Le* = 0.8 case, which is consistent with the observations made from Figure 1 and the expected picture of the flame surface associated with *Le* < 1. In the *Le* = 1.0 and 1.2 cases, the negative mean contribution of (*T*<sup>1</sup> + ··· + *T*10) acts to reduce the convexity of the positively curved bulges. By contrast, the negative mean (*T*<sup>1</sup> + ··· + *T*10) values at negatively curved locations in the *Le* = 1.0 and 1.2 cases act to increase the concavity of the negatively curved cusps. However, these structures are unstable and are eventually smoothed out due to increased *Sd* at negatively curved locations. This effect is particularly strong in the *Le* = 1.2 case because of large temperature (and thus also burning rate) at the negatively curved regions.

#### *4.5. Relations of the Terms of the Curvature Transport Equation with Local Curvature*

The mean values of *<sup>T</sup>*1−<sup>5</sup> <sup>×</sup> (*δ*<sup>2</sup> *th*/*SL*) conditional on *<sup>κ</sup><sup>m</sup>* × *<sup>δ</sup>th* and *<sup>κ</sup><sup>g</sup>* × *<sup>δ</sup>*<sup>2</sup> *th* for the *Le* = 0.8, 1.0 and 1.2 cases are shown in Figure 7 for the *c* = 0.8 isosurface. The corresponding variations of *<sup>T</sup>*6−<sup>10</sup> <sup>×</sup> (*δ*<sup>2</sup> *th*/*SL*) conditional on *<sup>κ</sup><sup>m</sup>* × *<sup>δ</sup>th* and *<sup>κ</sup><sup>g</sup>* × *<sup>δ</sup>*<sup>2</sup> *th* on the *c* = 0.8 isosurface are shown in Figure 8 for these cases. Figure 7 indicates high magnitudes of *T*<sup>1</sup> and *T*<sup>3</sup> are mostly obtained for concave topologies on the flame surface for the *Le* = 0.8 case. The same trend is seen for *T*1, *T*<sup>2</sup> and *T*<sup>4</sup> in the *Le* = 1.0 and 1.2 flames. The large magnitudes of *T*<sup>5</sup> are associated with the cup convex flame topologies in all cases, and the magnitudes of these terms increase with decreasing Lewis number.

It is evident from Figure 8 that the terms *T*6−<sup>10</sup> in the *Le* = 0.8 case show high magnitudes for cup concave (i.e., *κ<sup>m</sup>* < 0 and *κ<sup>g</sup>* > 0) and saddle concave (i.e., *κ<sup>m</sup>* < 0 and *κ<sup>g</sup>* < 0) topologies. However, large magnitudes of these terms are mostly obtained for cup convex and saddle convex flame topologies in the *Le* = 1.0 and 1.2 cases. Moreover, large negative magnitudes of *T*<sup>8</sup> and large positive magnitudes of *T*<sup>9</sup> and *T*<sup>10</sup> can be discerned for small (both positive and negative) values of *κ<sup>g</sup>* for weakly negatively curved zones (i.e., *κ<sup>m</sup>* < 0) in the *Le* = 1.0 case.

#### *4.6. Modelling Implications*

It can be seen from Figures 4–8 that the net flame propagation contribution to the curvature *κ<sup>m</sup>* transport plays significant roles for all flames considered here. This suggests that displacement speed *Sd* is not only interlinked with curvature (i.e., a negative correlation exists between *Sd* and *κ<sup>m</sup>* [32–41]) but it also affects flame wrikling through the curvature evolution. Moreover, it was demonstrated elsewhere [67,69,83–86] that the curvature *κ<sup>m</sup>* and its interrelation with displacement speed *Sd* plays a key role in the Flame Surface Density (FSD) and scalar dissipation rate (SDR) transports. Furthermore, the analysis of curvature evolution also reveals that the displacement speed induced terms are principally responsible for the generation of negatively curved cusps in the flames with *Le* < 1.

**Figure 7.** Variations of the mean values of *<sup>T</sup>*1−<sup>5</sup> <sup>×</sup> (*δ*<sup>2</sup> *th*/*SL*) (1st–5th row) conditional on *κ<sup>m</sup>* × *δth* and *<sup>κ</sup><sup>g</sup>* <sup>×</sup> *<sup>δ</sup>*<sup>2</sup> *th* for *Le* = 0.8, 1.0 and 1.2 flames (1st–3rd column) for *c* = 0.8 isosurface.

**Figure 8.** Variations of the mean values of *<sup>T</sup>*6−<sup>10</sup> <sup>×</sup> (*δ*<sup>2</sup> *th*/*SL*) (1st–5th row) conditional on *κ<sup>m</sup>* × *δth* and *<sup>κ</sup><sup>g</sup>* <sup>×</sup> *<sup>δ</sup>*<sup>2</sup> *th* for *Le* = 0.8, 1.0 and 1.2 flames (1st–3rd column) for *c* = 0.8 isosurface.

#### **5. Conclusions**

The effects of characteristic Lewis number on the statistical behaviours of the different terms of the curvature transport equation have been analysed based on three-dimensional compressible DNS data of spherically expanding turbulent premixed flames with *Le* = 0.8, 1.0 and 1.2. The statistically spherical flames with *Le* = 0.8, 1.0 and 1.2 had the same initial radius before they were allowed to interact with initially homogeneous isotropic decaying turbulence. It has been found that the flame surface area and volume-integrated burning rate increase with decreasing *Le*, which is consistent with several previous findings. The greater extent of flame wrinkling in the *Le* = 0.8 case is reflected in the wider range of both positive and negative curvatures than in the corresponding *Le* = 1.0 and 1.2

cases where the joint PDF remains almost symmetrical about *κ<sup>m</sup>* = 0 on the reactant side, but skews gradually toward positive values of *κ<sup>m</sup>* in the reaction and hot product zones. The PDFs of curvature for the *Le* = 0.8 case show higher probabilities of finding sharply negatively curved cusps than in the corresponding *Le* = 1.0 and 1.2 cases. Moreover, the saddle topologies have been found to be more frequent in the *Le* = 0.8 case than in the other cases considered in this analysis.

It has been found that the mean contributions of flame normal gradient of normal strain rate and the added strain rate due to flame displacement speed (i.e., *T*<sup>2</sup> and *T*9) assume negative values (i.e., promote concavity of the flame surface) towards the unburned gas and positive values (i.e., promote convexity of the flame surface) on the burned gas side of the flame. By contrast, the mean contributions arising from flow strain rate gradient and the flame normal gradient of the added strain rate (i.e., *T*<sup>4</sup> and *T*7) tend to promote positive and negative curvatures on reactant and product sides of the flame. The mean added curvature stretch term *T*<sup>8</sup> exhibits negative mean values throughout the flame, while the curvature flow stretching term *T*<sup>3</sup> remains positive throughout the flame front.

The mean added curvature stretch term *T*<sup>8</sup> conditional on curvature shows different behaviour in response to the changes in *Le*. The added curvature stretch term *T*<sup>8</sup> shows negative mean values for *κ<sup>m</sup>* < 0 and negligible mean value of *T*<sup>8</sup> is obtained for *κ<sup>m</sup>* > 0 in the *Le* = 0.8 case but positive (negative) mean values of *T*<sup>8</sup> are obtained for positive (negative) *κ<sup>m</sup>* values in the *Le* = 1.0 and 1.2 cases. The mean value of curvature flow stretching term *T*<sup>3</sup> remains positive for both *κ<sup>m</sup>* > 0 and *κ<sup>m</sup>* < 0 but the magnitude of *T*<sup>3</sup> associated with *κ<sup>m</sup>* < 0 increases significantly with decreasing *Le*. The terms due to normal gradient of added strain rate *T*<sup>9</sup> and the curl of added vorticity *T*<sup>10</sup> assume positive (negative) values for negative (positive) *κ<sup>m</sup>* values. The net mean contribution of the terms arising from flame propagation (i.e., (*T*<sup>6</sup> + ··· + *T*10)) dominates over the net mean contribution of the terms due to background fluid motion (i.e., (*T*<sup>1</sup> + ··· + *T*5)) for the negatively curved locations but the opposite behaviour has been observed for the positively curved zones in the *Le* = 0.8 case. However, in the *Le* = 1.0 and 1.2 cases the net mean flame propagation contribution (*T*<sup>6</sup> + ··· + *T*10) dominates over the mean contribution of the flow terms (*T*<sup>1</sup> + ··· + *T*5) especially for high magnitudes of *κm*. For the *Le* = 1.0 and 1.2 cases, the net mean contribution of (*T*<sup>1</sup> + ··· + *T*10) remains negative for high positive curvatures but this contribution assumes weak negative values for negatively curved regions. By contrast, in the *Le* = 0.8 case, the net mean contribution of (*T*<sup>1</sup> + ··· + *T*10) assumes negative values for *κ<sup>m</sup>* < 0 and mild positive values are obtained for *κ<sup>m</sup>* > 0, which is indicative of promoting sharply negatively curved cusps and positively curved bulges with large radii of curvature. This tendency is further augmented by high temperatures (also burning rates) in the positively curved and low temperatures (also burning rates) in the negatively curved regions in the *Le* = 0.8 case, which is expected in the presence of thermo-diffusive instability. By contrast, highly positively curved bulges are not promoted in the *Le* = 1.0 and 1.2 cases, and weak negative mean contributions of (*T*<sup>1</sup> + ··· + *T*10) at the negatively curved zones tend to produce negatively curved cusps which are eventually smoothed by large values of *Sd* in these regions. Thus, flame propagation tends to smoothen the flame wrinkles induced by turbulence in the *Le* = 1.0 and 1.2 cases and this effect is stronger in the *Le* = 1.2 case due to high temperature (and thus high reaction rate) at the negatively curved regions. Furthermore, it has been found that flame propagation plays a pivotal role in the curvature evolution irrespective of the characteristic Lewis number and thus the interrelation between displacement speed and curvature needs to be explicitly accounted for in the context of FSD and SDR closures, especially in order to predict thermo-diffusive instability effects for *Le* < 1 flames. This model development along with further analysis of curvature evolution using detailed chemistry and transport-based DNS data will form the basis of further investigations.

**Author Contributions:** Conceptualization, N.C. and M.K.; simulation, M.K.; Postprocessing code development, A.A., N.C.; formal analysis, A.A., N.C.; writing—original draft preparation, A.A., N.C.; writing—review and editing, N.C., M.K.; visualization, A.A., N.C. M.K.; supervision, N.C.; funding acquisition, N.C., M.K.

**Funding:** This research was funded by Engineering and Physical Sciences Research Council (EPSRC) (EP/K025163/1) and the German Research Foundation (DFG, KL1456/5-1).

**Acknowledgments:** We received computational support from ARCHER, Rocket High Performance Computing and Gauss Centre for Supercomputing (grant: pn69ga).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
