A Lagrangian analysis approach is used to examine the effects of heat release on the dynamics of the enstrophy during highly turbulent premixed combustion. The analysis is performed using data from a direct numerical simulation of a statistically planar premixed methane–air flame at a Karlovitz number of 100. Through cumulative, conditional, and correlation analyses, we show, consistent with prior studies, that vortex stretching and baroclinic torque both increase enstrophy at these highly turbulent conditions, while viscous transport and dilatation both lead to enstrophy destruction. However, although vortex stretching and viscous transport are individually an order of magnitude greater than all other terms in the enstrophy budget, the cumulative and combined effect of these two terms along Lagrangian trajectories is roughly only twice as large as the combined cumulative effect of dilatation and baroclinic torque. Moreover, trajectories that exhibit an increase in enstrophy through the flame are found to frequently have cumulative contributions from budget terms outside a single standard deviation of the mean contribution, indicating that enstrophy production at such highly turbulent conditions is associated with relatively infrequent but large values of dynamical terms. Lagged correlations further reveal a small but measurable contribution of baroclinic torque in enstrophy production, but these increases are overwhelmed, on average, by concurrent decreases in enstrophy due to viscous transport and dilatation. Taken together, these results provide further understanding of enstrophy dynamics in highly turbulent premixed flames.

## I. INTRODUCTION

The effects of turbulent premixed flames on the characteristics and dynamics of the enstrophy, $\Omega =\omega i\omega i$ (where $\omega i=\epsilon ijk\u2202uk/\u2202xj$ is the vorticity vector, $\epsilon ijk$ is the cyclic permutation tensor, and *u _{k}* is the velocity vector) have received considerable attention over the past decade. Because it is a gradient-based quantity, the enstrophy provides a measure of turbulence properties at small scales, and the dynamical effects of the flame on the turbulence appear explicitly in the enstrophy transport equation given by

Here, $D/Dt=\u2202/\u2202t+uj\u2202/\u2202xj$ is the material (or Lagrangian) derivative, $Sij=(1/2)(\u2202ui/\u2202xj+\u2202uj/\u2202xi)$ is the strain rate tensor, *ρ* is the density, *p* is the pressure, and *τ _{kl}* is the viscous stress tensor given by

where *μ* and *μ _{B}* are the dynamic and bulk viscosities, respectively, and

*δ*is the Kronecker delta.

_{kl}The first two terms on the right-hand side of Eq. (1) are generally nonzero in both reacting and non-reacting three-dimensional (3D) turbulent flows and represent vortex stretching and viscous transport processes, respectively. The last two terms, representing dilatation and baroclinic torque, are compressible flow effects that can be significant in the overall dynamics of the enstrophy during premixed combustion due to heat release from chemical reactions.

A number of recent computational and experimental studies have documented the characteristics and dynamics of enstrophy across a range of premixed flame conditions (e.g., spanning different Karlovitz, Ka, and Damköhler, Da, numbers), configurations (e.g., statistically planar^{1–6} and jet^{7,8} flames, as well as bluff-body stabilized flames^{9} and swirl burners^{10,11}) and locations within the flame (e.g., in the preheat and reaction zones). These studies have consistently demonstrated the importance of flame-induced dilatation, baroclinic torque, and enhanced viscous transport (due to the increase in gas phase viscosity with increasing temperature) on the dynamics of the enstrophy during premixed combustion; the reviews by Lipatnikov and Chomiak^{12} and Steinberg *et al.*^{13} provide comprehensive summaries of these and other studies.

Although much has been learned about vorticity and enstrophy dynamics in turbulent premixed combustion over the past decade, certain questions remain, including: (*i*) What is the *cumulative* effect along fluid pathlines of different physical phenomena on the enstrophy?; (*ii*) When the creation or destruction of enstrophy is observed in a flame, what is the dominant physical cause?; and (*iii*) What is the temporal correspondence between different physical phenomena and changes in the enstrophy? These questions are of particular interest for the highly turbulent conditions (generally corresponding to high Ka) found in many next-generation aerospace propulsion systems.^{14,15} Recent results from Geikie and Ahmed^{9} and Kazbekov *et al.*^{10,11} indicate that, in realistic configurations where mean pressure gradients are present (specifically, bluff-body stabilized flames and high-swirl burners, respectively), baroclinic torque can be a substantial contributor to flame-induced turbulence at highly turbulent conditions. By contrast, baroclinic torque has been found to contribute only weakly on average to the enstrophy dynamics at high Ka in statistically planar flames,^{1,4,5} although Poludnenko^{16} showed that baroclinic torque can still lead to intermittently large enstrophy production, even in statistically planar configurations.

To address each of the questions noted above, the temporal history of the enstrophy evolution along fluid particle pathlines is required, thereby suggesting a Lagrangian, as opposed to purely Eulerian, analysis approach. There is a considerable historical precedent for the use of a Lagrangian framework in the study of turbulent premixed combustion, including the early application of the Lagrangian direct interaction approximation by O'Brien^{17} to model mixing of chemically reactive scalars. More recently, Lagrangian analyses have been used to study the dynamical history of transient combustion phenomena, such as extinction, auto-ignition, pocket formation, and deflagration to detonation transition. Such analyses have included the study of both flame features (described as flame “particle” or “surface” tracking by Uranakara *et al.*,^{18} Dave *et al.*,^{19,20} and You and Yang^{21}) as well as velocity and thermochemical properties within fluid particles, the latter being the more traditional definition of a “Lagrangian” analysis in a turbulent flow context.

With respect to fluid particle tracking, Day *et al.*^{22,23} found that thermochemistry along fluid pathlines was complex and non-monotonic (i.e., fuel consumption did not occur monotonically with increasing time, as is the case in laminar flames). Subsequently, Steinberg and collaborators^{24,25} performed an analysis of velocity gradient dynamics and fluid particle residence times in premixed jet flames, finding that residence times between reactant and product temperatures decrease in turbulent flames, as compared to laminar flames. Through an analysis of highly turbulent H_{2}–air premixed combustion in statistically planar configurations, Hamlington *et al.*^{26} showed that non-monotonicity and thermochemical complexity are caused by enhanced molecular diffusion due to large turbulence-induced gradients in thermochemical quantities. From the perspective of turbulence characteristics, Geikie and Ahmed^{9} and Kazbekov *et al.*^{10,11} used Lagrangian approaches to examine enstrophy dynamics in bluff-body stabilized flames and swirl burners, respectively. To date, however, a Lagrangian analysis of enstrophy dynamics in the physically simpler statistically planar configuration has yet to be performed.

In this study, we calculate enstrophy dynamics along Lagrangian fluid particle pathlines using data from a direct numerical simulation (DNS) of a statistically planar turbulent premixed methane–air flame. We first briefly describe the numerical simulation and the Lagrangian analysis approach used to calculate the enstrophy dynamics. We then outline results and conclusions from the Lagrangian analysis, with a particular focus on answering the three questions outlined above. Finally, conclusions and important directions for future research are provided at the end.

## II. METHODS

### A. Description of numerical simulations

The simulation performed here is identical in its configuration to the simulation described and examined previously by Darragh *et al.*^{27} for the analysis of fluid particle dispersion in a turbulent premixed flame. The simulation models a 3D statistically planar stoichiometric methane–air premixed flame in a prismatic unconfined domain, with an aspect ratio of $1\xd71\xd716$. Additional details of the simulation are provided by Darragh *et al.*,^{27} and here, we give a relatively brief overview of the simulation setup.

The simulation was performed using Athena-RFX,^{16,28,29} which is a second-order-in-time, third-order-in-space code that solves the compressible, reactive-flow Navier–Stokes equations^{30} using an unsplit corner transport upwind scheme^{31,32} on a fixed, equispaced, 3D mesh with a nonlinear HLLC Riemann solver and piecewise parabolic spatial reconstruction. The bulk viscosity, $\mu B$, was set equal to the dynamic viscosity, *μ*, and an ideal gas equation of state was used to relate thermodynamic variables. Large-scale periodic forcing was used throughout the domain to sustain turbulence during the simulation, as described and validated by Poludnenko and Oran,^{29} resulting in statistically stationary homogeneous isotropic turbulence that interacts with the premixed flame. A 19-step reduced chemical mechanism^{33} was used to model the stoichiometric methane–air reaction at atmospheric conditions, and mixture averaging was used to compute molecular transport coefficients. Despite the use of mixture averaging in this study, it should be noted that individual chemical species will each have different Lewis numbers. The effects of these different Lewis numbers are of greater significance in studies of scalar and scalar gradient dynamics, as well as in cases where the analysis is conditioned on specific chemical species compositions. The present analysis of enstrophy dynamics conditioned on temperature is likely to be less dependent on this variability in Lewis number between species, although this topic does warrant further attention in future research, particularly since the Lewis number has been shown by Dopazo *et al.*^{34} to affect enstrophy dynamics in premixed combustion, resulting in an increased contribution of baroclinic torque. Further details on the computational setup of Athena-RFX for simulations of multi-species chemically reacting systems are provided by Hamlington *et al.*^{26} and Xu *et al.*^{35}

A uniform $256\xd7256\xd74096$ grid with cell size $2.7\xd710\u22123$ cm was used in the simulation, resulting in roughly half a grid cell per Kolmogorov scale in the reactants and six grid cells per Kolmogorov scale in the products. There are 16 grid cells per laminar flame thermal width, $\delta L=(Tp\u2212Tr)/|\u2207T|L,max$, where $Tr$ and $Tp$ denote the reactant and product temperatures, and $|\u2207T|L,max$ is the maximum gradient of temperature, *T*, in the laminar case. Given the resolution of the turbulent Kolmogorov and flame scales in this simulation, it can reasonably be described as a DNS, and physical characteristics of the present DNS are summarized in Table I.

L | 0.700 37 cm | Domain width |

$Tr$ | 300 K | Unburnt temperature |

$pr$ | $1.013\u200925\xd7106$ erg/cm^{3} | Unburnt pressure |

$\rho r$ | $1.12\xd710\u22123$ g/cm^{3} | Unburnt density |

$\delta L$ | $4.38\xd710\u22122$ cm | Laminar thermal width |

$SL$ | 37.2 cm/s | Laminar flame speed |

$\tau L$ | $1.18\xd710\u22123$ s | Laminar flame timescale |

$\u2113$ | $2.05\xd710\u22121$ cm | Unburnt integral scale |

$U\u2113$ | 785 cm/s | Unburnt integral velocity |

τ _{L} | $5.922\xd710\u22124$ s | Eddy turnover time |

Da | 0.50 | Damköhler number |

Ka | 100 | Karlovitz number |

$\eta r$ | $1.382\xd710\u22123$ cm | Unburnt Kolmogorov length |

$\eta p$ | $1.734\xd710\u22122$ cm | Burnt Kolmogorov length |

$\tau \eta r$ | $1.174\xd710\u22125$ s | Unburnt Kolmogorov time |

$\tau \eta p$ | $6.340\xd710\u22125$ s | Burnt Kolmogorov time |

$Rer$ | 87 | Unburnt Reynolds number |

$Rep$ | 28 | Burnt Reynolds number |

L | 0.700 37 cm | Domain width |

$Tr$ | 300 K | Unburnt temperature |

$pr$ | $1.013\u200925\xd7106$ erg/cm^{3} | Unburnt pressure |

$\rho r$ | $1.12\xd710\u22123$ g/cm^{3} | Unburnt density |

$\delta L$ | $4.38\xd710\u22122$ cm | Laminar thermal width |

$SL$ | 37.2 cm/s | Laminar flame speed |

$\tau L$ | $1.18\xd710\u22123$ s | Laminar flame timescale |

$\u2113$ | $2.05\xd710\u22121$ cm | Unburnt integral scale |

$U\u2113$ | 785 cm/s | Unburnt integral velocity |

τ _{L} | $5.922\xd710\u22124$ s | Eddy turnover time |

Da | 0.50 | Damköhler number |

Ka | 100 | Karlovitz number |

$\eta r$ | $1.382\xd710\u22123$ cm | Unburnt Kolmogorov length |

$\eta p$ | $1.734\xd710\u22122$ cm | Burnt Kolmogorov length |

$\tau \eta r$ | $1.174\xd710\u22125$ s | Unburnt Kolmogorov time |

$\tau \eta p$ | $6.340\xd710\u22125$ s | Burnt Kolmogorov time |

$Rer$ | 87 | Unburnt Reynolds number |

$Rep$ | 28 | Burnt Reynolds number |

Periodic boundary conditions were used in the *x* and *y* directions for the duration of the DNS. Prior to initialization of the flame, periodic boundaries were used in the *z* direction, with open zeroth-order extrapolation boundaries used after ignition. The domain was initialized with cold reactants at a temperature of $Tr=300$ K, and turbulence was allowed to develop over one large-scale eddy turnover time, $\tau L=L/UL$, where *L *=* *0.700 37 cm is the width of the domain and *U _{L}* is the turbulent velocity at the scale of the domain width

*L*. After the turbulence was fully developed, a planar laminar flame was initialized near the center of the domain with the flame normal in the

*z*direction. The simulation was then allowed to develop for an additional

*τ*before 3D Eulerian fields of velocity, density, pressure, and temperature were output at a high rate for the Lagrangian analysis.

_{L}Given the relatively high turbulence intensity of $U\u2113/SL=21$ in these simulations, the Karlovitz number, $Ka=\tau L/\tau \eta r=100$, was high and the corresponding Damköhler number, $Da\u2261\tau L/\tau L=0.50$, was low, where $\tau L=\delta L/SL$ is the laminar flame characteristic timescale given in Table I. As a result, the small-scale turbulent fluctuations occur on short timescales relative to the characteristic timescale of the flame. These flame conditions place the present simulation in the broken reaction zone regime.^{36} Although the turbulence is intense, both the bulk flow and turbulent Mach numbers were quite small, and there were subsequently no shocks or shocklets in the domain (recent research^{37} has, however, shown that Lagrangian analyses can be accurately performed in flows with such compressible features).

### B. Description of Lagrangian analysis and trajectory initialization

Fluid particle trajectories were calculated using the Lagrangian algorithm described by Darragh *et al.*^{38} and Hamlington *et al.*,^{26} with velocities obtained from the 3D Eulerian velocity field output from the DNS. Fluid particles were independently evolved using fifth-order Runge–Kutta time integration to solve the ordinary differential equation $dxn/dt=u(xn)$, where $xn(t)$ denotes the trajectory of the *n*th particle and **u** is the 3D Eulerian velocity field from the DNS. Akima splines were used for spatial and temporal interpolation during the integration procedure, giving at least second-order interpolation accuracy for all variables. Verification and validation of this algorithm are provided by Darragh *et al.*,^{38} and the first demonstration of this algorithm for fluid particle trajectory calculations using Eulerian data from Athena-RFX is outlined by Hamlington *et al.*^{26}

In the present study, fluid particles were initialized in 256 × 256 grids within *x*–*y* planes at 15 different *z* locations near the reactants, roughly at the start of the preheat zone. Each of these 983 040 particles were then evolved forward in time for one *τ _{L}*. Because this study is primarily focused on the effects of heat release by the flame on enstrophy dynamics, only particles that underwent complete combustion were considered in the analysis. This was accomplished by restricting the analysis to particles that reach a temperature of 500 K (indicating that reactions have begun to occur), followed at a later time by a temperature of 1800 K. This resulted in a total of 239 068 particles considered in this study.

Through this conditioning procedure, we ensured that all particles examined in the analysis were of a similar type, regardless of their initial spatial location. In particular, each particle analyzed underwent combustion with a corresponding temperature rise from 500 to 1800 K, accomplished via heat release resulting from chemical reactions. There are other particles that do not undergo combustion during the analysis period, but these particles are of lesser interest here since the corresponding enstrophy evolution would be similar to that obtained for particles in homogeneous isotropic turbulence (the type of turbulence in the reactants) in the absence of combustion. Moreover, the requirement that the particles reach a temperature of 500 K in order to be included in the analysis is equivalent, for this statistically stationary flame, to examining only particles initialized at locations where the temperature is at 500 K, that then subsequently undergo combustion.

With the Lagrangian trajectories computed, terms in the enstrophy governing equation in Eq. (1) were calculated from the Eulerian fields output by the DNS and then interpolated to the trajectories. This allowed us to recover the dynamics of the enstrophy along fluid particle pathlines. Gradients were computed by using Akima splines to interpolate quantities of interest to cell edges, followed by the calculation of flux differences to obtain gradients. The time derivative of Ω required to compute $D\Omega /Dt$ was calculated by first smoothing the enstrophy and then using a first-order backward finite difference. Using this approach, we computed each of the four terms on the right-hand side of Eq. (1) along each of the Lagrangian trajectories, and the resulting normalized time-dependent quantities are denoted here by

where $A=\u27e8(\Omega SmnSmn)1/2\u27e9$ is a normalizing factor that accounts for the intensity of the turbulence prior to its interaction with the flame. A similar normalization was used by Hamlington *et al.*^{1} for an Eulerian analysis of enstrophy dynamics in premixed flames, and this factor is computed as an average over the portion of trajectories in the reactants (i.e., for temperatures less than 310 K). Although it is traditional in studies of non-reacting flows to separate the viscous transport term in Eq. (4) into separate diffusive and dissipative parts, an unambiguous separation of this type is difficult in premixed flames where both the density and viscosity are spatially and temporally varying, and the dilatation is non-zero. This results in a number of different terms when expanding $D$ from Eq. (4), including terms involving gradients of the density and viscosity. Consequently, for simplicity in the present study we examine the viscous transport term in Eq. (4) as a whole.

It should be noted that, as a result of the numerical forcing used to maintain statistical stationarity of the DNS, as well as numerical diffusion introduced by the computational algorithm, the transport equation for Ω in Eq. (1) should also, from a practical standpoint, include an additional residual term. This term is obtained in the present analysis from the Navier–Stokes momentum conservation equation, which is given by

where *N _{i}* represents the residual. This term is calculated in the present study by using the 3D Eulerian field output by the DNS to subtract the pressure gradient and viscous transport terms on the right-hand side of Eq. (7) from the material derivative of velocity; this approach has also previously been suggested by Schranner

*et al.*

^{39}as a method to estimate numerical diffusion. The time derivative of velocity is computed using a second-order centered finite difference, while the gradients required in the material derivative are again computed by interpolating quantities to cell interfaces using Akima splines and then computing flux differences.

Given the resulting Eulerian fields of *N _{i}*, we then compute the curl of these fields and construct the corresponding normalized residual term in the enstrophy transport equation, given by

This time-dependent quantity is combined in the following analysis of enstrophy dynamics along fluid particle pathlines with the other terms in Eqs. (3)–(6).

## III. RESULTS

### A. Eulerian fields

Instantaneous two-dimensional (2D) Eulerian fields of temperature and enstrophy near the flame are shown in Fig. 1. Consistent with prior studies of statistically planar premixed flames,^{1,4,5} there is more spatial variability in the temperature field on the unburnt reactant side of the flame, as compared to the burnt product side. This is also indicated by the greater “wrinkling” of the *T *=* *400 K isocontour, corresponding to higher values of the enstrophy, as compared to the *T *=* *2000 K isocontour. In general, enstrophy is larger in the reactants and decreases through the flame into the products, again consistent with prior studies of enstrophy magnitude in turbulent premixed flames.^{1,4,5} The enstrophy field also shows greater small-scale structure in the reactants, which is destroyed as the temperature increases due to the flame; this is indicative of the suppression of small-scale motions that has been quantified in prior studies of the multi-scale structure of turbulence during premixed combustion.^{30,40–42}

Figure 2 shows corresponding 2D Eulerian fields for each of the normalized budget terms from the right-hand side of Eq. (1), where the definitions of each term are given in Eqs. (3)–(6) and (8). The primary balance indicated by the fields in Fig. 2 is between vortex stretching, $S$, which is generally positive throughout the flow, and viscous transport, $D$, which is generally negative. This dominant balance is consistent with previous studies of highly turbulent statistically planar premixed flames.^{1,4,5}

Dilatation, $L$, and baroclinic torque, $B$, shown in Figs. 2(c) and 2(d), are generally smaller than the vortex stretching and viscous transport terms and have their largest magnitudes in the vicinity of the flame. Although both of these fields exhibit both positive and negative contributions to the enstrophy dynamics, the dilatation has more prevalent negative values, indicating enstrophy destruction, while the baroclinic torque is more generally positive, indicating enstrophy generation. Additionally, the baroclinic torque and dilatation appear to have opposite signs at each location; this correlation will be examined in more detail in Sec. III D.

The magnitudes of each of the terms in Fig. 2 are smaller in the products than in the reactants, and the characteristic sizes of the flow features are larger in the products than in the reactants. Both of these results are likely due to the associated changes in enstrophy through the flame shown in Fig. 1(b), particularly given the dependence of each of the budget terms on the vorticity magnitude. The connection between the enstrophy and budget terms will be made more quantitative in Sec. III D, where temporal cross-correlations are used to determine the relationship between changes in enstrophy and the budget terms.

The residual field, $N$, shown in Fig. 2(e) qualitatively resembles fields for the other budget terms, particularly with respect to the change in magnitude and feature size across the flame. However, the residual field also exhibits pixelated noise, which is most apparent in the reactants where the spatial resolution of the DNS is lower relative to the Kolmogorov scale and numerical diffusion is more pronounced. The coherent large-scale structures in the products in Fig. 2(e) are likely associated with the numerical forcing.

### B. Lagrangian analysis of enstrophy dynamics

As noted in Sec. II B, in this study we only consider particles that undergo combustion, determined to be particles that experience a temperature of 500 K followed by a temperature of 1800 K. The particles that remain in the reactants or otherwise do not reach a temperature of 1800 K will not be considered in the following analyses. Additionally, all mentions of a particle reaching a specific temperature are assumed to be the first time the particle reaches the given temperature, despite the fact that it is possible for the temperature of the particle to decrease and then increase at a later time. As will be seen in the results, the temperature is frequently non-monotonic along fluid particle trajectories, and so there can be several time instances at which the particle reaches a temperature of *T *=* *500 K, for instance. We thus choose the first time at which the particle reaches a specified temperature in the analysis. The duration between 400 and 2000 K will also include trajectories that reach 400 and 1800 K, but may not reach 2000 K before the end of data collection.

Figure 3 shows the temperature, enstrophy, and normalized enstrophy budget terms along a single representative trajectory. At lower temperatures, the particle is still in the reactants where the enstrophy magnitude is larger than in the products, and thus the budget terms are larger in magnitude than at later times and higher temperatures. The rapid variations observed in the budget terms correspond to the small-scale flow features evident in Figs. 1 and 2, which are more prevalent in the reactants. These variations are also characterized by timescales much smaller than the chemical timescale, due primarily to the high Ka and low Da nature of this flame. In particular, Fig. 3(a) shows that the fluid particle undergoes combustion (corresponding to a rapid increase in temperature) over a timescale of roughly 0.1 ms, but the fluctuations in the total budget terms shown in Fig. 3(c) take place over timescales that are over an order of magnitude smaller. Such a wide disparity in chemical and turbulence timescales is characteristic of highly turbulent flames where Ka is high. In the region of intermediate temperatures (i.e., between roughly *T *=* *400 and 1800 K), Fig. 3(b) shows that the baroclinic torque and dilatation become larger in magnitude, with the baroclinic torque tending to have positive values and the dilatation negative values. Additionally, all terms appear to change magnitude together, indicating the presence of correlations between budget terms, which are again examined further in Sec. III D.

The balance between the left- and right-hand sides of Eq. (1) are shown over the entire trajectory in Fig. 3(c), where the right-hand side also includes the residual. There is good agreement between the two resulting time histories, indicating that the enstrophy budgets are accurately calculated using the procedure outlined in Sec. II B.

Temperature, enstrophy, and the normalized enstrophy budget terms are shown for 25 randomly selected trajectories in Fig. 4, where each trajectory is centered on the time at which the fluid particle first reaches a temperature of 500 K in Figs. 4(a)–4(g) and a temperature of 1400 K in Figs. 4(h)–4(n). The centering on 500 K in Figs. 4(a)–4(g) provides an indication of the different particle evolutions from the start of the analysis period, which begins for each trajectory when the corresponding particle first reaches a temperature of 500 K. These figures also indicate the nearly isothermal nature of the particles as they travel through the reactants prior to the onset of combustion. The centering at 1400 K in Figs. 4(h)–4(n) more clearly shows the particle evolutions near the region of peak heat release in the flame, as determined from the laminar flame profile. Although it is difficult from Fig. 4 to discern the specific details of any individual trajectory, this figure does give a sense of the range of values attained along the different trajectories, as well as the abrupt change in the magnitudes of the different normalized budget terms across the flame.

Figures 4(a) and 4(h) show that the temperature time histories along the particle trajectories are frequently non-monotonic, which has also been noted in previous studies of highly turbulent premixed flames.^{22,23,26} In general, Figs. 4(b) and 4(i) show that enstrophy tends to decrease through the flame; however, the decrease is once again not monotonic, and some periods of enstrophy increase are observed as different terms compete to increase and decrease the enstrophy through the flame. A conditional analysis of only trajectories that experience an increase in enstrophy while in the flame is presented in Sec. III C.

Figures 4(c)–4(g) and 4(j)–4(n) show the normalized enstrophy budget terms for each of the 25 trajectories. From the trajectories centered at 1400 K, shown in Figs. 4(j)–4(n), the magnitudes of all budget terms tend to be greater prior to peak heat release (which occurs close to 1400 K), with all terms reducing in magnitude at later times, corresponding to higher temperatures. While all terms are capable of both increasing and decreasing the enstrophy, these trajectories more often have enstrophy increased by stretching and baroclinic torque, and decreased by dilatation and viscous dissipation. The residual term tends to be more evenly distributed in positive and negative contributions.

To aggregate these Lagrangian results and develop an overall understanding of the enstrophy characteristics and dynamics, each trajectory is centered at the time when 1400 K is first reached and averages over the temperature, enstrophy, and normalized budget terms are taken over all burnt trajectories, as shown in Fig. 5. This method results in fewer trajectories included in the averages at early or late times because not all trajectories reach 1400 K at the exact same time, although at *t *=* *0, all trajectories are included in the average.

In Fig. 5(a), the average enstrophy is largest at early times when particles are in the reactants and decreases as the particles move through the flame, with almost all of the decrease occurring prior to *t *=* *0 (corresponding to *T *=* *1400 K). At these early times, the average vortex stretching is positive and the average viscous transport is negative, with both having similar magnitudes that decrease as *t *=* *0 is approached. The average baroclinic torque and average dilatation are also mirrors of each other, with opposite values of similar magnitude both peaking immediately prior to *t *=* *0. Away from *t *=* *0, both of these terms make much smaller contributions to the overall enstrophy dynamics. However, the contribution near *t *=* *0 does not result in a large average effect on the enstrophy, suggesting that these terms may be more relevant for individual trajectories, rather than the average evolution.

The range of values attained by each of the normalized enstrophy budget terms is indicated by the probability density functions (pdfs) shown in Fig. 6(a). The pdfs, which are calculated over all trajectories during the period of time that each particle spends between 400 and 2000 K, show that all terms tend to peak near a value of zero, although the pdfs of stretching and baroclinic torque are skewed more toward positive values, while the pdfs of viscous transport and dilatation are skewed more toward negative values. It should be noted that the pdf of the residual term is also skewed to negative values, reflecting its association with numerical diffusion and the destruction of enstrophy.

The enstrophy budget along Lagrangian trajectories can also be examined by considering the cumulative effect of each budget term on the evolution of the enstrophy. This cumulative analysis is motivated by the first of the three questions posed in Sec. I; namely, what is the *cumulative* effect along fluid pathlines of different physical phenomena on the enstrophy? This cumulative effect is determined by integrating each budget term between 400 and 2000 K along a trajectory, and Fig. 6(b) shows pdfs of the resulting cumulative effects aggregated over all trajectories. As compared to the pdfs of the instantaneous budget terms in Fig. 6(a), the cumulative effect pdfs in Fig. 6(b) more clearly show the net positive contributions of vortex stretching and baroclinic torque to enstrophy production, and the net negative contributions of viscous transport and dilatation to enstrophy destruction. The cumulative dilatation effect is the weakest of all effects, and the residual (associated with numerical diffusion) provides the second largest contribution to enstrophy destruction, after the viscous transport effect.

From a more quantitative perspective, the mean of each of the cumulative effects averaged over all trajectories between 400 and 2000 K is shown in Table II. Both dilatation and baroclinic torque are similar in magnitude and opposite in sign, and the same is true for vortex stretching and viscous transport. The magnitude of both baroclinic torque and dilatation is roughly an order of magnitude below the effects of stretching and viscous transport. In a similar type of comparison, Bobbitt *et al.*^{5} and Hamlington *et al.*^{1} both found around an order of magnitude difference when comparing baroclinic torque and dilatation with stretching and viscous dissipation, consistent with the present observations.

Mean integrated quantity . | Effect . |
---|---|

$\u27e8\u222bT=400T=2000S\u27e9$ | $2.62\xd7105$ |

$\u27e8\u222bT=400T=2000D\u27e9$ | $\u22122.41\xd7105$ |

$\u27e8\u222bT=400T=2000L\u27e9$ | $\u22123.67\xd7104$ |

$\u27e8\u222bT=400T=2000B\u27e9$ | $4.92\xd7104$ |

$\u2212\u27e8\u222bT=400T=2000N\u27e9$ | $\u22125.63\xd7104$ |

$\u27e8\u222bT=400T=2000(S+D)\u27e9$ | $2.07\xd7104$ |

$\u27e8\u222bT=400T=2000(L+B)\u27e9$ | $1.26\xd7104$ |

$\u27e8\u222bT=400T=2000(S+L+B+D\u2212N)\u27e9$ | $\u22122.30\xd7104$ |

Mean integrated quantity . | Effect . |
---|---|

$\u27e8\u222bT=400T=2000S\u27e9$ | $2.62\xd7105$ |

$\u27e8\u222bT=400T=2000D\u27e9$ | $\u22122.41\xd7105$ |

$\u27e8\u222bT=400T=2000L\u27e9$ | $\u22123.67\xd7104$ |

$\u27e8\u222bT=400T=2000B\u27e9$ | $4.92\xd7104$ |

$\u2212\u27e8\u222bT=400T=2000N\u27e9$ | $\u22125.63\xd7104$ |

$\u27e8\u222bT=400T=2000(S+D)\u27e9$ | $2.07\xd7104$ |

$\u27e8\u222bT=400T=2000(L+B)\u27e9$ | $1.26\xd7104$ |

$\u27e8\u222bT=400T=2000(S+L+B+D\u2212N)\u27e9$ | $\u22122.30\xd7104$ |

Several combinations of normalized budget terms are also shown in Table II. These combinations (i.e., stretching and viscous dissipation, dilatation and baroclinic torque, and all five budget terms) are added together before the integration and averaging are performed. By combining the terms in this way, the combined cumulative effect of stretching and viscous dissipation are found to be similar to the combined cumulative effect of dilatation and baroclinic torque over this range of temperatures, with a difference of slightly less than a factor of two, instead of an order of magnitude when considering each of these terms alone. The combination of baroclinic torque and dilatation represents both terms that are non-zero due to the presence of the flame, and thus their net positive contribution can be considered a manifestation of flame-generated vorticity. That is, the non-negligible, positive, average, cumulative contribution to enstrophy confirms that enstrophy can be enhanced by the flame.

### C. Conditional analysis

In order to address the second question posed in Sec. I and determine the dominant physical cause when the creation or destruction of enstrophy is observed in a flame, we perform a conditional analysis of the trajectories based on whether the enstrophy increases or decreases during combustion. While most trajectories decrease in enstrophy through the flame due to the increasing temperature and corresponding increase in viscosity, it is also possible for the enstrophy to increase along a fluid particle pathline that passes through the flame. In particular, 28% of the trajectories that underwent combustion and had a temperature range of 400 to 2000 K experienced an enstrophy increase, while 72% decreased in enstrophy. Consequently, despite the overall average decrease in enstrophy indicated in Fig. 5, over a quarter of the trajectories undergo an increase in enstrophy during the combustion process, potentially suggesting the presence of flame generated enstrophy. Figure 7 shows the distribution of ratios between the enstrophy at the end of the analysis period to the enstrophy at the start. Trajectories that decrease in enstrophy are more likely than those that increase, and the decrease is more likely to be large. Conversely, increases in enstrophy are more likely to be small (i.e., the ratio is more likely to be close to unity, with large increases being less probable).

Table III separates the full set of trajectories into two cases: trajectories that increase in enstrophy between 400 and 2000 K, and trajectories that decrease in enstrophy over the same temperature range. The fraction of trajectories with a positive cumulative effect from each budget term for both cases is then given, along with the fraction of trajectories outside a single standard deviation (std) of the mean of the given budget term. These fractions are based on the total number of trajectories that increase or decrease in enstrophy, not the total number of trajectories in the study; however, the mean and standard deviation are computed using the set of all trajectories. Additionally, Table III also shows the fraction of trajectories with at least one budget term outside one standard deviation.

. | Enstrophy increase . | Enstrophy decrease . |
---|---|---|

Fraction of trajectories | 0.28 | 0.72 |

Positive cumulative $S$ | 0.99 | 0.98 |

Positive cumulative $L$ | 0.17 | 0.05 |

Positive cumulative $B$ | 0.82 | 0.79 |

Positive cumulative $D$ | 0.01 | 0.01 |

Positive cumulative $N$ | 0.44 | 0.27 |

Outside one std of $S$ mean | 0.17 | 0.05 |

Outside one std of $L$ mean | 0.16 | 0.10 |

Outside one std of $B$ mean | 0.19 | 0.13 |

Outside one std of $D$ mean | 0.17 | 0.06 |

Outside one std of $N$ mean | 0.11 | 0.05 |

At least one term outside one std of mean | 0.32 | 0.20 |

. | Enstrophy increase . | Enstrophy decrease . |
---|---|---|

Fraction of trajectories | 0.28 | 0.72 |

Positive cumulative $S$ | 0.99 | 0.98 |

Positive cumulative $L$ | 0.17 | 0.05 |

Positive cumulative $B$ | 0.82 | 0.79 |

Positive cumulative $D$ | 0.01 | 0.01 |

Positive cumulative $N$ | 0.44 | 0.27 |

Outside one std of $S$ mean | 0.17 | 0.05 |

Outside one std of $L$ mean | 0.16 | 0.10 |

Outside one std of $B$ mean | 0.19 | 0.13 |

Outside one std of $D$ mean | 0.17 | 0.06 |

Outside one std of $N$ mean | 0.11 | 0.05 |

At least one term outside one std of mean | 0.32 | 0.20 |

Table III indicates that, among the trajectories that increase in enstrophy, all show a greater fraction of positive cumulative budget effects except for the viscous transport term, which is roughly equivalent to the set of trajectories that decreased in enstrophy. Additionally, for all budget terms, a greater fraction of cumulative effects lie outside one standard deviation of the mean for trajectories that increase in enstrophy when compared to those that decrease. Furthermore, the fraction of trajectories with at least one term outside a single standard deviation of the mean is roughly 1.6 times greater than the fraction amongst trajectories with decreasing enstrophy. This suggests that these less common but large contributions are associated with increases in enstrophy along a trajectory during the combustion process. The importance of these extreme values suggests a connection between enstrophy production during combustion and turbulence intermittency, whereby large values of velocity gradient quantities are observed with increasing probability (as compared to Gaussian probabilities) for increasing Reynolds number.^{43} Thus, although enstrophy production is not observed on average, there are locally large instances of enstrophy production during the combustion process, leading to an overall increase in enstrophy in more than a quarter of the trajectories.

With respect to the probabilistic distributions of the different normalized budget terms, Fig. 8 shows pdfs of both the raw value and cumulative effect of each term for trajectories between 400 and 2000 K. Trajectories that increased in enstrophy are separated from those that decreased in enstrophy. Figure 8 shows that pdfs for both the raw values and the cumulative effects have similar peaks, both in height and location, for both sets of trajectories. The largest difference between the two sets of trajectories is found in the tails of the pdfs, again consistent with a connection between enstrophy production and turbulence intermittency. In particular, pdfs for all budget terms show that larger-in-magnitude values are more probable for trajectories that increase in enstrophy than those that decrease in enstrophy. This, along with the results shown in Table III, suggests that, in the highly turbulent statistically planar premixed flame examined here, the greatest increases in enstrophy are associated with relatively infrequent, large magnitude enstrophy budget effects.

### D. Cross-correlation analysis

Finally, to determine the temporal correspondence between different physical phenomena and changes in the enstrophy (thereby addressing the third question posed in Sec. I), we examine here the temporal cross-correlations between temperature, enstrophy, and the normalized budget terms. For two generic time-dependent quantities $a(tn)$ and $b(tn)$ defined along a Lagrangian trajectory, where *t _{n}* denotes the $nth$ discretized time along the trajectory, the cross correlation for a temporal separation of $k\Delta t$, where

*k*is the discrete number of lagged time steps and $\Delta t=1\xd710\u22126$ s is the time step between successive values along the Lagrangian trajectories, is defined as

^{44,45}

where *N* is the total number of time steps in each series and $||\xb7||$ denotes the L2 norm. Using Eq. (9), we compute a temporal cross correlation for each trajectory individually, spanning over all possible values of the number of lagged time steps *k* along the trajectory. Figure 9 shows the resulting average correlation, denoted $\u27e8cab(k\Delta t)\u27e9$, which is determined by averaging the cross-correlations $cab(k\Delta t)$ from Eq. (9) over all trajectories, giving an average that is itself a function of the number of lagged time steps, *k*. This is done for several selected correlations between pairs of normalized budget terms, temperature, and enstrophy. The notation is such that a negative time shift implies that the first subscript (e.g., *a* in *c _{ab}*) is large before the second subscript (e.g.,

*b*), while a positive time shift implies the opposite. Correlations between raw values are shown in Fig. 9(a), while Fig. 9(b) shows correlations between the absolute value of quantities that are not strictly positive.

The trajectory-averaged correlation between temperature and dilatation, $\u27e8cTL\u27e9$, is found to be large with a positive time lag, indicating that large temperatures follow large dilatation (i.e., dilatation is larger on the reactant side of the flame). This matches the results for the time evolutions of the temperature and enstrophy budget terms along trajectories shown in Figs. 3–5, where dilatation is large in magnitude before the temperature becomes large. A similar result is found for $\u27e8cT|L|\u27e9$ in Fig. 9(b).

For the correlation between the enstrophy and dilatation, $\u27e8c\Omega L\u27e9$ in Fig. 9(a) is larger in magnitude and more negative at negative time shifts, indicating that enstrophy tends to be larger in the reactants, as expected, and that when it is destroyed, dilatation tends to be large. A stronger correlation with zero time shift is found between enstrophy and the absolute value of dilatation, $\u27e8c\Omega |L|\u27e9$, in Fig. 9(b).

The correlation between dilatation and baroclinic torque, $\u27e8cLB\u27e9$ in Fig. 9(a), is negative with little to no time shift, thus supporting the findings in Fig. 5, where both terms peak at a similar time and work in opposition, with baroclinic torque enhancing enstrophy while the dilatation diminishes it. With respect to the magnitudes of dilatation and baroclinic torque, Fig. 9(b) further shows that $\u27e8c|L||B|\u27e9$ is even greater than $\u27e8cLB\u27e9$. These results indicate that enstrophy destruction by dilatation and production by baroclinic torque are typically concurrently large along trajectories, such that even when these two effects are large individually, they oppose each other and result in little overall change in the enstrophy. This corresponds to an overall reduction in the effect of the flame on enstrophy, despite the fact that the dilatation and baroclinic torque may each be individually large in magnitude.

Figure 9(a) shows that there is a small negative peak with no time shift in the correlation between the residual term and baroclinic torque, $\u27e8cNB\u27e9$, indicating little to no correlation between the two terms. However, when $\u27e8c|N||B|\u27e9$ is considered in Fig. 9(b), the correlation strengthens. This is explained by the fact that baroclinic torque is typically large where the pressure and density gradients are large. The residual term also tends to be largest where large gradients are located, since part of the residual term is due to numerical diffusion acting to smooth out gradients. Density is one of the solution variables in Athena-RFX, and thus gradients of density are attenuated such that $N$ should be large wherever gradients of density are large (i.e., where baroclinic torque is also large).

The correlation between temperature and baroclinic torque, $\u27e8cTB\u27e9$ in Fig. 9(a), is similar to the correlation between temperature and dilatation, with a stronger correlation at positive time lags, indicating that baroclinic torque is large before increases in temperature occur. Consequently, we find that baroclinic torque is most prominent immediately prior to heat release from chemical reactions, corresponding to flame locations roughly at the transition between the preheat and reaction zones (close to *T *=* *1400 K in the present case). The same conclusion is obtained from $\u27e8cT|B|\u27e9$ in Fig. 9(b), and is also consistent with the observation in Hamlington *et al.*^{1} that the conditional average of the baroclinic torque (where the conditioning is based on the reactant mass fraction) is largest near the transition between the preheat and reaction zones.

Finally, the correlation between enstrophy and baroclinic torque, $\u27e8c\Omega B\u27e9$ in Fig. 9(a), is relatively weak, with a small peak at zero time lag. When $\u27e8c\Omega |B|\u27e9$ in Fig. 9(b) is instead considered, a stronger correlation is found between enstrophy and the magnitude of baroclinic torque at negative time lags, indicating that enstrophy is smaller after baroclinic torque acts. However, despite the fact Fig. 9 shows that enstrophy decreases after baroclinic torque acts, we cannot conclude that baroclinic torque is responsible for the decrease in enstrophy. Instead, the observed decrease in enstrophy is more likely due to concurrently large dilatation and enhanced viscous transport, both of which have large negative magnitudes when the baroclinic torque has a large positive magnitude; this conclusion is supported, for example, by the strongly negative correlation between $L$ and $B$ shown in Fig. 9(a). Once more, despite the fact that baroclinic torque does act to enhance enstrophy through the flame, it is typically large and positive where other effects (e.g., dilatation and viscous transport) are also large but negative, thereby mitigating the overall enhancement of enstrophy by combustion processes.

## IV. CONCLUSIONS

Enstrophy dynamics in a highly turbulent statistically planar premixed methane–air flame have been examined using a Lagrangian analysis approach. In addition to vortex stretching, viscous transport, dilatation, and baroclinic torque terms in the enstrophy governing equation, a residual term from the DNS was included in the analysis, allowing for improved agreement between the left- and right-hand sides of the governing equation and ensuring the accuracy of the results. Through the resulting analysis, progress has been made in addressing each of the three questions identified in Sec. I.

For the first question, related to the cumulative effects of different dynamical terms on the enstrophy along trajectories, we found, as in previous studies of highly turbulent statistically planar premixed flames, that the dominant terms in the enstrophy dynamics were vortex stretching and viscous transport. In particular, the cumulative effects of each of these terms along Lagrangian trajectories were roughly an order of magnitude greater than dilatation and baroclinic torque. Consistent with previous studies, the average cumulative effect of stretching and baroclinic torque was found to increase enstrophy, whereas viscous transport and dilatation tended to decrease enstrophy along fluid pathlines. However, although vortex stretching and viscous transport were found to be the dominant individual terms in the enstrophy dynamics, when considered together, the combined cumulative effects of these terms were slightly less than a factor of two larger than the combined effects of dilatation and baroclinic torque, which are both nonzero primarily due to heat release from the flame. These results are consistent with those found in previous studies of statistically planar premixed flames,^{1,4,5} although they differ from results for statistically planar flames with smaller Lewis numbers^{34} and for more practical configurations with persistent mean pressure gradients,^{9–11} both of which reveal a greater importance of baroclinic torque in the overall enstrophy dynamics. There is, however, no inconsistency in any of these results, and the primary conclusion gained is that care should be taken when generalizing results for statistically planar flames with near-unity Lewis numbers to other conditions and configurations.

To address the second question in Sec. I, we performed a conditional analysis of the trajectories, where the conditioning was based on whether each trajectory exhibited an increase or decrease in enstrophy through the flame. This analysis showed that over a quarter of all trajectories had an increase in enstrophy through the flame and that these increases were due primarily to relatively infrequent, but quite large, values of the enstrophy budget terms. Thus, while it is true that, on average, enstrophy decreases through the present highly turbulent, statistically planar premixed flame, we have also shown that individual trajectories are capable of increasing in enstrophy and that these increases are not especially rare events, even in highly turbulent statistically planar premixed flames. Moreover, the importance of extreme events in the overall enstrophy budget suggests a connection to turbulence intermittency, which is more pronounced at the high Reynolds numbers characteristic of many highly turbulent combustion problems. By contrast, decreases in enstrophy were found to be less likely the result of extreme values of enstrophy budget terms.

The observation that the baroclinic torque is large immediately prior to increases in temperature suggest that this term is largest near the preheat zone immediately prior to the transition to the reaction zone. This observation is most likely due to the fact that large values of the baroclinic torque require both large gradients in thermodynamic properties, as well as substantial misalignment between these gradients. These two conditions are approximately satisfied at the transition to the reaction zone where the turbulence, which leads to the misaligned gradients, has not yet weakened appreciably due to viscous transport and dilatation, but the thermodynamic gradients are beginning to increase due to heat release by reactions.

Finally, to address the third question in Sec. I related to the temporal correspondence between different physical phenomena, we performed a temporal cross correlation analysis of temperature, enstrophy, and normalized enstrophy budget terms along particle trajectories. This analysis showed that baroclinic torque and dilatation are strongly negatively correlated, and that smaller vorticity magnitudes are associated with large magnitudes of both the baroclinic torque and dilatation. Correlations between baroclinic torque and enstrophy show that baroclinic torque has a small but measurable effect in increasing enstrophy, although the simultaneously occurring reductions in enstrophy due to viscous transport and dilatation tend to overwhelm the increases due to baroclinic torque, resulting in an average net decrease in enstrophy through the flame.

This study suggests a number of directions for future research. The importance of pressure and density gradients in setting the baroclinic torque, and the observation that the residual tends to be large when the baroclinic torque is large, suggests the need for further increases in spatial resolution to further reduce the residual and obtain an even more accurate analysis of different dynamical effects on the enstrophy. It is not expected, however, that the calculations of the baroclinic torque have been substantially impacted in the present analysis, but rather that further increases in resolution will reduce the contribution of numerical dissipation to the residual.

With respect to other future directions, although we have included viscous transport from Eq. (4) in the analysis, this term can be further separated into contributions from viscous dissipation and viscous diffusion, with and without viscosity and density gradients, for both solenoidal and dilatational modes. This analysis would serve to further elucidate the precise viscous effects on enstrophy through the flame. The present Lagrangian enstrophy analysis should also be extended to other statistically planar premixed flames occurring at lower turbulence intensities and Karlovitz numbers, where effects due to the flame are expected to be more pronounced.

The effects of different Lewis numbers, including flames where species Lewis numbers vary considerably, as well as different fuels, should also be considered, particularly given the dependence of vorticity dynamics on Lewis number^{34} and of the turbulence kinetic energy on fuel type^{46} noted in previous studies. It would also be beneficial to repeat the cumulative, conditional, and correlation analyses outlined here in more realistic configurations where mean pressure-gradients are present and baroclinic torque is expected to be much stronger, particularly at high turbulence intensities. Recent research by Geikie and Ahmed^{9} and Kazbekov and Steinberg,^{10,11} in particular, has shown that baroclinic torque becomes an increasingly dominant contributor to enstrophy production at highly turbulent conditions, potentially resulting in a net increase in enstrophy through the flame, by contrast to the reduction in enstrophy observed in statistically planar configurations. The extension of the present analysis to these and other flows would also enable more general statements to be made regarding the importance of baroclinic torque and flame generated vorticity at highly turbulent conditions, as well as provide insights for the development of subgrid-scale models for large eddy simulations of practical combustion problems involving high turbulence intensities.

## ACKNOWLEDGMENTS

R.D. was supported by a National Defense Science and Engineering Graduate Fellowship. M.A.M. was supported by a NSF Graduate Research Fellowship. C.A.Z.T. and P.E.H. were supported, in part, by AFOSR Award No. FA9550–17-1–0144 and NSF Award No. 1847111. Helpful discussions with Dr. Alexei Poludnenko are gratefully acknowledged.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

_{2}O PLIF