Enstrophy transport conditional on local flow topologies in different regimes of premixed turbulent combustion

Enstrophy is an intrinsic feature of turbulent flows, and its transport properties are essential for the understanding of premixed flame-turbulence interaction. The interrelation between the enstrophy transport and flow topologies, which can be assigned to eight categories based on the three invariants of the velocity-gradient tensor, has been analysed here. The enstrophy transport conditional on flow topologies in turbulent premixed flames has been analysed using a Direct Numerical Simulation database representing the corrugated flamelets (CF), thin reaction zones (TRZ) and broken reaction zones (BRZ) combustion regimes. The flame in the CF regime exhibits considerable flame-generated enstrophy, and the dilatation rate and baroclinic torque contributions to the enstrophy transport act as leading order sink and source terms, respectively. Consequently, flow topologies associated with positive dilatation rate values, contribute significantly to the enstrophy transport in the CF regime. By contrast, enstrophy decreases from the unburned to the burned gas side for the cases representing the TRZ and BRZ regimes, with diminishing influences of dilatation rate and baroclinic torque. The enstrophy transport in the TRZ and BRZ regimes is governed by the vortex-stretching and viscous dissipation contributions, similar to non-reacting flows, and topologies existing for all values of dilatation rate remain significant contributors.

Turbulent flows are inherently rotational in nature and the extent of this rotation is often quantified in terms of vorticity ω → (i.e. curl of velocity vector → u ). An important and characteristic feature of turbulent flows is the vortex-stretching mechanism, which is responsible for spreading turbulent velocity fluctuations over different length scales 1 . Thus, fundamental understanding of the enstrophy (i.e. ω ω Ω = → ⋅ → /2) transport is of key importance for the purpose of characterisation of the energy cascade in turbulent flows. In order to better understand premixed flame-turbulence interaction, the seemingly chaotic turbulence may be split into different generic topologies. Understanding which topologies make significant contribution to the enstrophy transport in different regimes of combustion might therefore help to gain deeper insight into strongly coupled phenomena of chemical reaction and turbulent flow and mixing.The current analysis aims to identify the canonical flow configurations which make dominant contributions to the enstrophy transport in different regions of the flame brush for different combustion regimes. This information can be utilised to design simplified experimental configurations representing dominant flow topologies to gain further insight into the flame-turbulence interaction.
Turbulent flows exhibit eight generic canonical local flow configurations underneath an apparently random fluid motion 2, 3 , and these distinct flow topologies are categorised depending on the values of first, second and third invariants (i.e. P, Q and R, respectively) of the velocity gradient tensor, A ij = ∂u i /∂x j = S ij + W ij , where the symmetric strain-rate tensor is S ij = 0.5(A ij + A ji ) and the anti-symmetric rotation rate tensor is W ij = 0.5(A ij − A ji ).
probability of finding focal topologies decreases from the unburned to the burned gas side of the flame front, which was also observed by Wacks et al. 20 for droplet combustion. A recent analysis by Wacks et al. 21 focused on local flow topology distributions for a direct numerical simulations (DNS) database consisting of three nominally thermo-diffusively neutral H 2 -air turbulent premixed flames representative of the corrugated flamelets, thin reaction zones and broken reaction zones regimes 22 . This analysis revealed significant qualitative differences in flow topology distribution within the flame front for flames representing different regimes of premixed turbulent combustion.
In turbulent flow, the enstrophy (i.e. Ω = ω i ω i /2 where ω i is the i th component of vorticity) and vorticity fields affect the statistical behaviours of the second and third invariants, Q and R, which in turn affect the distributions of local flow topologies S1-S8. The instantaneous transport equation of Ω = ω i ω i /2 is given by [23][24][25][26][27][28][29][30] : The term T I is the vortex-stretching contribution to the enstrophy transport, whereas term T II arises due to misalignment of gradients of density and viscous stresses. The term T III is responsible for molecular diffusion and dissipation of enstrophy due to viscous action, whereas the term T IV signifies the dissipation of enstrophy due to dilatation rate. The term T V is the baroclinic torque term which arises due to misalignment between pressure and density gradients.
A number of recent studies [23][24][25][26][27][28][29] investigated the Reynolds-averaged enstrophy transport equation in turbulent premixed flames by Reynolds-averaging eq. 4: Here q indicates the Reynolds-averaged value of a general quantity q. In the current analysis all the Reynolds averaged quantities are evaluated by ensemble averaging the quantities in the statistically homogeneous directions (which are the directions normal to the direction of the mean flame propagation) when the flame attains a quasi-steady state. A similar approach was adopted in several previous analyses (e.g. refs 26 and 27 and references therein).
Hamlington et al. 23 concentrated on different mechanisms of vorticity generation and the alignment of vorticity with local principal strain rates in the flames representing the thin reaction zones regime of premixed turbulent combustion. Chakraborty 25 showed that the alignment of vorticity with local principal strain rates is significantly affected by the regime of combustion and the characteristic Lewis number (i.e. ratio of thermal diffusivity to mass diffusivity). Both studies 23,25 reported a predominant alignment of vorticity with the intermediate principal strain rate in turbulent premixed flames, which is similar to non-reacting turbulent flows 31 . However, the relative alignment of vorticity with the most extensive and most compressive principal strain rates is significantly affected by the strength of dilatation rate 25 , which is influenced by the regime of combustion and the characteristic Lewis number. Hamlington et al. 23 showed that enstrophy decays significantly in the burned gas across the flame brush, whereas Treurniet et al. 24 demonstrated an opposite trend for the flames with high density ratio between the unburned and burned gases. This behaviour has been explained by Lipatnikov et al. 26 by analysing the terms of enstrophy and vorticity transport equation for weakly turbulent premixed flames representing the corrugated flamelets regime.
Recent analyses by Chakraborty et al. 27 and Dopazo et al. 30 demonstrated that the characteristic Lewis number significantly affects the vorticity generation within the flame, and a combination of augmented flame normal acceleration and high extent of flame wrinkling for small values of Lewis number may give rise to a significant amount of vorticity generation within the flame brush due to baroclinic torque. Bobbit and coworkers 28,29 demonstrated that the enstrophy transport in statistically planar flames propagating in homogeneous isotropic turbulence for large values of Karlovitz number is governed by the relative balance between the vortex-stretching and viscous dissipation similar to its non-reacting counterpart, and also revealed that the choice of chemical reaction mechanism does not affect the qualitative nature of the enstrophy transport. The enstrophy field in turbulent premixed flames using cinema-stereoscopic particle image velocimetry (PIV) measurements of rim-stabilised turbulent premixed flames has been investigated 32-34 and confirmed some of the observations based on DNS data.
The aforementioned analyses 23-29, 31-34 provided invaluable insight into the flame-turbulence interaction, but the nature of the enstrophy transport conditional on generic local flow topologies S1-S8 in different premixed combustion regimes is yet to be analysed in detail. Such an analysis is expected to reveal the canonical flow configurations which make dominant contributions to the enstrophy transport for different combustion regimes.
In order to address the aforementioned gap in the existing literature and to meet the above objectives an existing detailed chemistry DNS database 21,35 has been considered consisting of three H 2 -air flames with an equivalence ratio of φ = 0.7, representative of the combustion processes in the corrugated flamelets, thin reaction zones, and broken reaction zones regimes of combustion. A detailed chemical mechanism 36 involving 9 steps and 19 chemical reactions is considered for this analysis. The unburned gas temperature T 0 is taken to be 300 K, which gives rise to an unstrained laminar burning velocity S L = 135.62 cm/s under atmospheric pressure. The inlet values of normalised root-mean-square turbulent velocity fluctuation u′/S L , turbulent length scale to flame thickness ratio l T /δ th , Damköhler number Da = l T S L /u′δ th , Karlovitz number Ka = (ρ 0 S L δ th /μ 0 ) 0.5 (u′/S L ) 1.5 (l T /δ th ) −0.5 and turbulent Reynolds number Re t = ρ 0 u′l T /μ 0 for all cases are presented in Table 1, where ρ 0 is the unburned gas density, μ 0 is the unburned gas viscosity, l T is the most energetic length scale, δ th = (T ad − T 0 )/max|∇T| L is the thermal flame thickness with T, T ad and T 0 being the instantaneous dimensional temperature, adiabatic flame temperature and unburned gas temperature, respectively, and the subscript 'L' is used to refer to unstrained laminar flame quantities.
The cases investigated in this study are representative of three regimes of combustion: case A: corrugated flamelets (Ka < 1), case B: thin reaction zones (1 < Ka < 100) and case C: broken reaction zones regime (Ka > 100) 22 . The Karlovitz number can be scaled as δ η ∼ Ka / th 2 2 (where η is the Kolmogorov length scale), which suggests that the turbulent eddies do not disturb the inner flame structure in case A, whereas turbulent eddies penetrate into the preheat zone in the thin reaction zones regime, but the reaction zone (with typical thickness of δ th /10) retains a quasi-laminar structure. In the broken reaction zones regime, the energetic turbulent eddies may penetrate into the reaction zone and disturb the chemical processes.

Results and Discussion
The distributions of natural logarithm (to cover the wide dynamic range) of normalised enstrophy field  Fig. 2. It is clearly seen that there are significant qualitative differences in the enstrophy and topology distributions between these cases. For example, cases B and C show a much larger range of structures, whereas mostly large scale structures are evident in case A. As seen in Table 1, Re t in case B is greater than in case A, whereas the value of l T remains the same. Thus, case B shows smaller structures than case A. For case C, the length scale remains smaller than case A, which along with higher Re t leads to smaller structures than in case A (note the different scaling of Case C). Figure 2 (top row) further shows that in case A the enstrophy value increases locally from the unburned to the burned gas side of the flame due to flame-induced vorticity generation. In contrast, the magnitude of enstrophy decreases from unburned to burned gas side of the flame for cases B and C. Furthermore, Fig. 2 indicates that S1, S2, S3 and S4 topologies are predominantly present in the burned gas in all cases, but the topology distribution within the flame is significantly different between cases A-C. This can be substantiated from Fig. 3 which shows the probability of finding each flow topology for different values of Favre-averaged reaction progress variable  c T , where the Favre-average of a general quantity q is given by ρ ρ =  q q/ . In actual RANS simulations, one obtains Favre-averaged reaction progress variable  c T and not the Reynolds averaged reaction progress variable c T . Furthermore, in statistically planar flames both Reynolds averaged and Favre-averaged values of reaction progress variable remain unique functions of the mean direction of flame propagation and thus all the terms in this configuration can be presented as a function of either c T or  c T . The findings and conclusions would have been qualitatively similar if the results were shown as a function of c T . As  c T is readily obtained from RANS simulations, the results are shown as a function of Favre-averaged reaction progress variable in this paper. Figure 3 indicates that S1-S4 topologies remain major contributors within the flame front for all three cases. As shown in Fig. 1, S1-S4 appear for all values of P, and thus they remain dominant contributors, which is consistent with previous findings by Cifuentes et al. 19 . The dilatation rate ∇ ⋅ → = − u P is predominantly positive within the flame and thus the topologies S5 and S6, which are typical of positive P (i.e. = −∇ ⋅ → P u), rarely occur within the flame, and a similar observation was recently made by Cifuentes et al. 19 . The probability of obtaining S7 topology, which is typical of positive dilatation rate (i.e. ∇ ⋅ → = − > u P 0), increases from the unburned to the burned gas side of the flame brush (i.e. with increasing  c T ) in all cases considered here. The strength of dilatation rate and likelihood of obtaining high positive values of ∇ ⋅ → u are significantly smaller in case C than that in cases A and B because energetic turbulent eddies penetrate into the reaction zone and disturb the chemical reaction, which in turn reduces the heat release rate and the magnitude of ∇ ⋅ → u 21 . The weakening of dilatation rate effects in case C is reflected in the fact that the S8 topology completely disappears in case C.
The variation of the normalised Reynolds-averaged enstrophy where i = 1, 2, ..., 8 for S1-S8 respectively) conditional on  c T is shown in Fig. 4. It is seen from Fig. 4 that δ Ω × S ( / ) th L 0 2 increases for the major part of the flame brush from the unburned gas side towards the burned gas side, before decreasing towards the far edge of the burned gas side ( ≈ . th L 0 2 decreases across the flame brush for cases B and C, which is consistent with the observation made from Fig. 2. The decay of enstrophy from the unburned to the burned gas side of the flame brush for thermo-diffusively neutral flames within the thin reaction zones regime is consistent with previous findings 23,27,30 . Figure 4 also shows the contributions of each individual topology (i.e. Ω i where i = 1, 2, ..., 8 for S1-S8) to the mean enstrophy as a function of  c T for cases A-C, indicating that the topologies S1-S4 (with the focal topologies Scientific REPORts | 7: 11545 | DOI:10.1038/s41598-017-11650-x S1 and S4 as the leading contributors), which are obtained for all possible values of P (see Fig. 1), contribute significantly to the mean enstrophy in all cases. Furthermore, the focal topology S7, which occurs only for positive dilatation rates, contributes significantly to the mean enstrophy Ω 0 in case A, but this trend weakens from case A to case C and the contribution of S7 remains negligible in case C.
In order to identify the physical mechanisms responsible for the difference in the enstrophy statistics between cases A-C, which are representative of combustion situations in three different combustion regimes, the variation of T I , T II , T III , T IV and T V with  c T for cases A-C are shown in Fig. 5. Significant differences are seen in the behaviour of the enstrophy transport between cases A, B and C. The term T III due to the viscous diffusion and dissipation acts as a leading order sink, and the magnitude of the viscous torque contribution T II remains small in comparison to T III for all cases. In case A, the baroclinic torque contribution T V acts as the major source term, and its magnitude remains comparable to that of T III . The contribution of the vortex-stretching term T I remains weak in case A and its magnitude remains small in comparison to that for the terms of T III and T V , whereas the dilatation term T IV   acts a major sink term in this case and its magnitude remains comparable to that of the combined viscous diffusion and dissipation term T III . The relative contributions of the dilatation term T IV and the baroclinic torque term T V are somewhat diminished in the flame representing the thin reaction zones regime (e.g. case B) and become negligible in the flame representing the broken reaction zones regimes of combustion (e.g. case C). In case B, the vortex-stretching term T I and the combined viscous diffusion and dissipation term T III remain the major source and sink terms, respectively. Although the baroclinic term T V and the dilatation term T IV continue to be source and sink terms in case B, their magnitudes are generally small in comparison to T I and T III , respectively. In case C, the transport of enstrophy is determined by the source and sink contributions of the vortex-stretching term T I and the combined viscous diffusion and dissipation term T III , respectively and the magnitudes of T II , T IV and T V remain negligible in comparison to T I and T III . The behaviour of the enstrophy transport in case C is qualitatively similar to a non-reacting flow in this configuration. The effects of dilatation rate and density change on the background fluid motion weaken with increasing Karlovitz number Ka, and thus the influences of the dilatation and baroclinic terms (i.e. T IV and T V ) weaken from the corrugated flamelets regime to the thin reaction zones regime to the broken reaction zones regime. The weaker contributions of T IV and T V in the thin reaction zones and broken reaction zones regime flames are consistent with previous findings [27][28][29][30] .
The vortex-stretching term T I can be expressed as: T e cos e cos e cos 2( ) where e α , e β and e γ are the most extensive, intermediate and the most compressive principal strain rates respectively, and the angles between ω → and the eigenvectors associated with e α , e β and e γ are given by α, β and γ, respectively. It was shown by Chakraborty 25 that ω → predominantly aligns with e β (i.e. high probability of |cosβ| ≈ 1), but the extent of alignment with e α increases in regions of high chemical heat release. The extent of alignment of ω → with e α (e γ ) decreases (increases) with decreasing Karlovitz number especially in corrugated flamelets regime flames 25 and a similar qualitative behaviour has been observed here. Interested readers are referred to Ref. 25 for the relevant explanations. The probability of increased alignment of ω → with the most compressive (i.e. most negative) and intermediate principal strain rates is responsible for the weak contribution of T I in case A. In premixed flames dilatation rate ∇ ⋅ → u remains predominantly positive as a result of thermal expansion due to heat release, which gives rise to a predominantly negative value of dilatation contribution It was demonstrated by Wacks et al. 21 that the magnitude and influence of dilatation rate weakens from case A to C and accordingly T IV plays a progressively less important role from case A to case C. It has been shown earlier that the cosine of the angle between ω → and ∇p × ρ ρ ∇ / 2 remains mostly positive within premixed flames 26,27 , and a similar behaviour has been found here (not shown for brevity), which gives rise to a positive contribution of the baroclinic torque in the mean sense (i.e. T V ) for all cases. A comparison between the enstrophy transport behaviours in cases A-C reveals that the baroclinic torque is responsible for the enstrophy generation in case A, which gives rise to local augmentation of enstrophy across the flame in this case (see Fig. 2). In cases B and C, only the vortex-stretching term T I acts as the major source and its magnitude decreases from the unburned to the burned gas side of the flame brush.
The term-by-term contributions of each individual topology to the enstrophy transport (i.e. T j i where i = 1, 2, …, 8 for S1-S8 and j = I, II, III, IV, V for T I , T II , T III , T IV and T V , respectively, such that ) as a function of  c T for cases A-C are shown in Fig. 6. It is seen that the focal topology S7 (which is typical of P < 0) associated with stretching is a leading positive contributor to the vortex-stretching term T I in the corrugated flamelets regime case A. Moreover, the focal topology S4 associated with stretching remains a significant positive contributor to the vortex-stretching term T I in case A. In this case the magnitudes of the contributions of the topologies S1, S2 (partially obscured by S3) and S3 to the vortex-stretching term T I remain significant, but these topologies act to induce negative values. The topologies S1, S2 and S4, which can be obtained irrespective of the value of P, play leading order roles in the vortex-stretching term T I in the thin reaction zones regime case B and also in the broken reaction zones regime case C.
The topologies S2, S7, S8 and S4 contribute significantly to the viscous torque term T II in the decreasing order of importance in the corrugated flamelets regime case A. The contribution of T II remains positive in case A and mostly positive except towards the unburned gas side of the flame brush in the thin reaction zones regime case B. The topologies S1, S2, S8, S7 and S4 contribute significantly to T II in case B and the positive contribution of T II towards the burned gas side of the flame brush originates principally due to the S7 topology. The topologies S4, S1 and S2, which can be obtained irrespective of the value of P, contribute significantly to the viscous torque term T II (in decreasing order of significance) in the broken reaction zones regime case C. Figure 6 suggests that the S4, S1 and S2 topologies are the leading order contributors to the combined viscous diffusion and dissipation term T III for all three cases, but S7 and, to a lesser extent, S3 topologies also contribute significantly to this term for the case representing the corrugated flamelets regime (i.e. case A). In the thin reaction zones and broken reaction zones regime cases (i.e. in cases B and C), S4, S1 and S2 topologies are responsible for most of the viscous diffusion and dissipation.
The S2, S7, S1, S4 and S8 (in decreasing order of significance) topologies contribute significantly to the dilatation rate T IV and baroclinic torque T V terms in the corrugated flamelets regime case A, whereas S7, S2, S1, S4 and S8 (in decreasing order of significance) topologies remain leading order contributors to T IV in the thin reaction zones regime case B. The topologies S2, S7, S4, S1 and S8 (in decreasing order of significance) are the major contributors to the baroclinic torque T V term in the thin reaction zones regime case B. In the broken reaction zones regime case C, the dilatation rate T IV and the baroclinic torque T V terms do not have much significance, but S4, S1 and S7 remain major contributors to these terms in this case.
The topologies which are significant contributors to T I , T II , T III , T IV and T V in cases A-C are summarised in Table 2. From the foregoing discussion and Table 2 it is evident that S1-S4, S7 and S8 topologies play key roles in the enstrophy transport in the corrugated flamelets regime flame, but the contributions of S7 and S8 topologies, which are obtained only for positive dilatation rates, weaken with an increase in Karlovitz number and for the broken reaction zones regime flame the enstrophy transport is governed by S1, S2 and S4 topologies which are obtained for conditions irrespective of the value of dilatation rate. Thus, simple experimental configurations satisfying canonical features of S7, S8 and S1-S4 are likely to be necessary and sufficient in order to extract fundamental information regarding the enstrophy transport in the corrugated flamelets regime, whereas a simple configuration with predominant features of S1, S2 and S4 will be sufficient for physical understanding of the enstrophy transport in the broken reaction zones regime.

Summary
An existing three-dimensional DNS database 21,35 containing three freely propagating statistically planar H 2 -air flames representing typical flame-turbulence interaction in the corrugated flamelets, thin reaction zones and broken reaction zones regimes of turbulent premixed combustion has been used to analyse the enstrophy transport conditional on flow topologies within the flame. The flow topologies have been characterized in terms of three invariants of velocity gradient tensor (P, Q and R), where the first invariant is the negative of dilatation rate and second and third invariants can be linked to strain rate and enstrophy, and their generation mechanisms 2,3,11 . In this analysis the flow topologies have been categorized in 8 types (i.e. S1-S8) depending on the location of velocity gradient tensor in P − Q − R space 2, 3 . It has been found that the weakening of dilatation rate with increasing Karlovitz number plays a key role in the enstrophy transport in turbulent premixed flames. The contributions to the enstrophy transport conditional on topology have been analysed in detail and it has been found that the enstrophy generation due to baroclinic torque weakens from the corrugated flamelets regime case to the broken reaction zones regime case. Further, the flow topologies S1-S4, which can be obtained for all values of dilatation rate, contribute significantly to the enstrophy and its transport in the thin reaction zones and broken reaction zones regimes of premixed turbulent combustion. However, the topologies (i.e. S7 and S8), which are obtained only for positive values of dilatation rate, also contribute significantly to the enstrophy transport in the corrugated flamelets regime.

Methodology
A three-dimensional complex chemistry (9 steps and 19 chemical reactions according to a detailed chemical mechanism 36 ) compressible flow DNS database of freely-propagating statistically planar turbulent H 2 -air premixed flames with φ = 0.7 has been considered for this analysis. An equivalence ratio of 0.7 is chosen because an H 2 -air mixture for this equivalence ratio is known to be thermo-diffusively neutral 36 , such that the additional effects of the preferential diffusion are eliminated. Interested readers are referred to Im et al. 35 and Wacks et al. 21 for detailed information on the DNS database used for the current analysis and here a brief description of the numerical methodology is provided. The spatial discretisation is carried out using an 8 th order central difference scheme for internal grid points and the order of differentiation gradually decreases to a one-sided 4 th order scheme at the non-periodic boundaries. A fourth order Runge-Kutta scheme is used for explicit time advancement. The flame is initialised by a steady 1D planar laminar flame profile, and a pre-computed auxiliary divergence free, homogeneous, isotropic turbulence field generated using a pseudo-spectral method 37 following Passot-Pouquet spectrum 38 is injected through the inlet. The mean inlet velocity has been gradually modified to match turbulent flame speed as the simulation progresses. The temporal evolution of flame area has been monitored and the flame is taken to be statistically stationary when the flame area no longer varies with time. Turbulent inflow and outflow boundaries are specified in the direction of mean flame propagation and the transverse boundaries are taken to be periodic. The non-periodic boundaries are specified using an improved Navier Stokes characteristic boundary conditions (NSCBC) technique 39 .
The domain size is taken to be 20 mm × 10 mm × 10 mm (8 mm × 2 mm × 2 mm) in cases A and B (case C) and the domain has been discretised by a uniform Cartesian grid of 512 × 256 × 256 (1280 × 320 × 320). The grid spacing was determined by the flame resolution, ensuring about 10 grid points across δ th , and in all cases the Kolmogorov length scale remains bigger than the grid spacing (i.e. η ≥ 1.5Δx where η and Δx are the Kolmogorov length scale and DNS grid spacing, respectively). Simulations have been carried out for 1.0t e , 6.8t e and 6.7t e (i.e. t e = l T /u ′ ) for cases A-C respectively, and this simulation time remains comparable to several previous analyses [40][41][42] .