Three-dimensional, rotational flamelet closure model with two-way coupling

Abstract

A new flamelet model is developed to be used for subgrid modelling and coupled with the resolved flow description for turbulent combustion. The model differs from current models in several critical ways. (i) Non-premixed flames, premixed flames or multi-branched flame structures are determined rather than prescribed. (ii) The effects of shear strain and vorticity on the flames are determined. (iii) The strain rates and vorticity applied at the subgrid level are directly determined from the resolved-scale strain rates and vorticity without the use of a contrived progress variable. (iv) The flamelet model is three-dimensional without need for assuming axisymmetry or planar geometry. (v) The effect of variable density is addressed in the flamelet model. Solutions to the Navier–Stokes equations and the associated scalar equations governing the flamelet model are obtained without the boundary-layer approximation. By appropriate coordinate transformation, a similar solution is found for the rotational flamelet model, reducing it to a system of ordinary differential equations. Vorticity is shown to create a centrifugal force on the subgrid counterflow that modifies the molecular transport rates and burning rate. Sample computations of the rotational flamelet model without coupling to the resolved flow are presented first to demonstrate the importance of the new features of the model. Scaling laws are presented for relating strain rates and vorticity at the subgrid level to those quantities at the resolved-flow level for coupling with large-eddy simulations or the time-averaged mean-flow level for Reynolds-averaged flows. The time-averaged behaviour of a simple turbulent flow is resolved with coupling to the rotational flamelet model. Specifically, a two-dimensional, multicomponent, time-averaged planar shear layer with variable density and energy release is employed using a mixing-length description for the eddy viscosity. Needs for future study are identified.

1. Introduction

Combustion in high-mass-flux chambers is the practical and major method for energy conversion for mechanical power and heating. Inherently, the high mass-flow rate leads to turbulent flow. Thereby, many length and time scales appear in the physics, making serious challenges for both computational and experimental analyses. For computations where the smallest scales typically cannot be resolved, the method of large-eddy simulation (LES) is employed wherein the smaller scales are filtered via integration over a window size commensurate with the computational mesh size that allows affordable computations. Consequently, the essential, rate-controlling, physical and chemical processes that occur on shorter scales than the filter size must be modelled. Those subgrid models must be properly coupled to the resolved LES flow field.

Current flamelet models that are used for LES or Reynolds-averaged Navier–Stokes (RANS) methods have some advantages. Typically, the flamelet equations are a system of ordinary differential equations (ODEs) that can be solved offline with solutions available in tabular form or through neural networks. The flamelet models can handle multi-species, multi-step oxidation kinetics without requiring small time steps during the solution of the resolved-scale fluid dynamics. Thus, for several reasons, savings of computational resources can be huge compared with direct numerical simulation. We aim here to retain these very attractive features while removing some less desirable features. Already, some progress has been made in extending the fundamental flamelet theory beyond its long-term limitation of a single-flame structure, two-dimensional (or axisymmetric) configuration and use of the uniform-density assumption. However, those advances still must be applied to LES or RANS. In addition, the flamelet theory must be advanced to consider shear strain and vorticity at the small scale of the flamelet; these are the vital forgotten physics in current flamelet modelling. Furthermore, the strain rates in the flamelet model are far from properly connected to the strain rates at the resolved scale. Attempts at corrections of these weaknesses are made here.

The focus here is on the types of turbulent flames found in the shear-driven flows of practical combustors. Thereby, we address flamelets as originally described by Williams (1975), namely ‘highly sheared small diffusion flamelets’ and ‘forming a turbulent flame brush which appears on the average to fill’ the flow domain. Our discussion will not cover some interesting work on corrugated or wrinkled flames. The goals in this paper are to improve the flamelet model by including several important physical effects that are commonly neglected in present models and to identify other issues, related to the coupling between the subgrid-scale physics and the resolved-scale (or time-averaged) physics, that require further study.

1.1. Existing flamelet theory

There is need to understand the laminar mixing and combustion that commonly occur within the smallest turbulent eddies. These laminar flamelet subdomains experience significant strain of all types: shear, tensile and compressive. Some important works exist here but typically for either counterflows with only normal strain or simple vortex structures in two dimensions or axisymmetry and often with a constant-density approximation; see Linan (1974), Marble (1985), Karagozian & Marble (1986), Cetegen & Sirignano (1988), Cetegen & Sirignano (1990), Peters (2000) and Pierce & Moin (2004). Linan and Peters focused on the counterflow configuration. Williams (1975) first established the concept of laminar flamelets in a turbulent diffusion flame structure. Karagozian and Marble examined a three-dimensional flow with radial inward velocity, axial jetting and a vortex centred on the axis. The flame sheet wrapped around the axis due to the vorticity. An interesting review of the early flamelet theory is given by Williams (2000). Generally, flamelet studies have focused on either premixed or non-premixed flames; a unifying approach to premixed, non-premixed and multi-branched flames has not been developed. That unifying approach is taken here.

Most flamelet studies have not directly considered vorticity interaction with the flamelet; see, for example, Linan (1974), Peters (2000), Williams (2000) and Pierce & Moin (2004). Mueller (2020) presented the flamelet model in a somewhat different mathematical framework but without the addition of new physical description. Williams (1975) first recognized the advantage of separating rotation (due to vorticity) and stretching by transformation to a rotating, non-Newtonian reference frame. He did not, however, examine the momentum consequences in the new reference frame which will be examined later here. The other works that have examined vortex–flame interaction have not separated the effects of stretching and rotation; see Marble (1985), Karagozian & Marble (1986), Cetegen & Sirignano (1988), Cetegen & Sirignano (1990) and Meneveau & Poinsot (1991).

Again, with the above-mentioned emphasis on turbulent flames in shear-driven flows as defined by Williams (1975), some interesting work on corrugated or wrinkled premixed laminar flames is not detailed here. The reader is guided to Swaminathan & Grout (2006), Steinberg & Driscoll (2009), Chakraborty & Swaminathan (2007) and Mura & Champion (2009) for examples of those studies.

The two-dimensional planar or axisymmetric counterflow configuration has become a foundation for the flamelet model. Local conversion to a coordinate system based on the principal strain-rate directions can provide the counterflow configuration in a general flow. Furthermore, the quasi-steady counterflow can be analysed by ODEs because the dependence on the transverse coordinate is either constant or linear, depending on the variable. Pierce and Moin modified the non-premixed flamelet counterflow configuration by fixing domain size and forcing flux to zero at the boundaries. Flamelet theory as a closure model for turbulent combustion is typically based on the tracking of two variables: a normalized conserved scalar and the strain rate. The latter is generally given indirectly through a progress variable. Mixture fraction is traditionally used for the conserved scalar.

The flamelet model has become a popular subgrid model for gas-fuelled combustors. Some development is also underway for the use of flamelets in spray combustion, but here we focus on the former type. The flamelet model for LES developed by Pierce & Moin (2004) was a substantial advancement through the introduction of the flamelet progress variable (FPV). Their approach has also been used by Ihme, Schmitt & Pitsch (2009), Nguyen, Popov & Sirignano (2018), Nguyen & Sirignano (2018, 2019) and others. Other works are based on the use of the original form developed by Peters (2000). Pierce and Moin extended that work in two ways. Firstly, the inclusion of both the upper and middle branches of the curve of flame temperature versus scalar dissipation rate allowed better representation of the unsteady details in turbulent combustion such as extinction and reignition. Secondly, they created a progress variable as a function (of the scalar dissipation rate) which is governed by a partial differential equation (added to the LES equations) with a chemical-rate source term determined through the flamelet model. Nguyen & Sirignano (2018) found that the inclusion of both branches for flame stability resulted in better agreement with experiment. Note that Nguyen & Sirignano (2018) addressed rocket combustion instability which places a greater demand on the flamelet model than most other applications. In addition to the velocity fluctuations due to turbulence, it becomes necessary to address very large fluctuations in velocity, pressure and temperature due to the nonlinear acoustics.

There are concerns with the existing model. While there is clear utility for the FPV approach, there are concurrently clear signs of incompleteness in the flamelet model and contradiction between the tenets of the model and the LES results produced with the model. A few examples of problems that require resolution are provided here. Firstly, the design of the flamelet model uses a counterflow configuration where only normal strain is imposed on the flame region by the ambient flow. Yet, the LES results show that the flame is embedded in a flow field with substantial shear strain rate, i.e. a flow with vorticity. Secondly, the current models assume that the inflowing streams in the counterflow are irrotational; yet, we know that the smaller scales in turbulent combustion are highly rotational. Thirdly, the classical flamelet model is two-dimensional (or axisymmetric) whereas it has been used in flows that are clearly three-dimensional. The model considers only a single-branched diffusion flame while its coupled use in LES commonly predicts multi-branched flames in qualitative agreement with experimental evidence. A fourth issue involves the quasi-steady assumption for the flamelet model that is used in a highly unsteady LES flow with a broadband noise.

These models are built around the postulate that the flamelets are always non-premixed (i.e. diffusion) flames or premixed flames. However, evidence of both non-premixed and premixed flames has been found in LES and experimental results. In fact, they can coexist in a multi-branched structure. Nguyen et al. (2018) and Nguyen & Sirignano (2018, 2019) employed the Pierce & Moin (2004) flamelet approach in the simulation of a single-injector rocket engine. They showed the importance of flamelets subject to high strain rates. However, contradictions occurred in that both premixed flames and non-premixed flames appeared in the predictions. In fact, they report multi-branched flames; in particular, the combination is often seen of a fuel-lean premixed flame branch with a branch consisting of a merged diffusion flame and fuel-rich premixed flame.

Note that the mixture fraction has been used widely as an independent variable to display non-premixed flamelet scalar variations; this cannot be useful for premixed flames. Sirignano (2021a) has shown that any conserved scalar can serve well as an independent variable to present scalar results for non-premixed and multi-branched flamelets.

Experiments and asymptotic analysis by Hamins, Thridandam & Seshadri (1985) showed that a partially premixed fuel-lean flame and a diffusion flame can coexist in a counterflow with opposing streams of heptane vapour and methane–oxygen–nitrogen mixture. Thus, a need exists for flamelet theory to address both premixed and non-premixed flames. Recently, Rajamanickam et al. (2019) provided an interesting three-dimensional triple-flame analysis, describing the effect of imposed normal strain on a multi-branched flame. While it did not consider shear strain, it was a helpful step followed by the work of Sirignano (2021c) where both shear strain and normal strain were considered.

The classical counterflow treatment by Linan (1974) and Peters (2000) has two opposing streams, fuel or fuel plus a chemically inert gas and oxidizer or oxidizer plus an inert gas. They considered uniform density. That critical assumption was relaxed by Sirignano (2019) for reacting flows and heated flows. Sirignano (2021a,c) with one-step kinetics and López-Cámara, Jordà Juanós & Sirignano (2020) with detailed kinetics address that single diffusion flame case. In addition, situations are addressed where the inflowing streams from $y_{\infty }$ and $y_{-\infty }$ may consist of a combustible mixture of fuel and oxidizer, thereby allowing another flame or two besides the simple diffusion flame to coexist. Sirignano (2021a) provides a counterflow analysis with three-dimensional strain and shows the possibility for a variety of flame configurations to exist depending on the compositions of the inflowing streams: (i) three flames including fuel-lean partially premixed, non-premixed (i.e. diffusion-controlled) and fuel-rich partially premixed; (ii) non-premixed and fuel-rich partially premixed; (iii) fuel-lean partially premixed and non-premixed; (iv) non-premixed; and (v) premixed. López-Cámara et al. (2020) extended the counterflow analysis to consider detailed kinetics for methane–oxygen detailed chemical kinetics and confirmed that combinations of premixed and non-premixed flames could exist in a multi-flame counterflow.

1.2. Relative orientations of principal strain axes, vorticity and scalar gradients

Both normal strain rate and shear strain rate are important. There is a strong need to study mixing and combustion in three-dimensional flows with both imposed normal strain and shear strain and therein imposed vorticity with global circulation. Shear strain can, in general, be decomposed into a normal strain and a rotation (whose rate is half of the vorticity magnitude). For example, a rectangular shape that is changed by shear strain can be viewed as a combination of deformation to a parallelogram caused by normal strain perpendicular to the diagonal and rotation of the diagonal caused by vorticity. The behaviour due to the strain and rotation becomes especially important on the smallest scales of turbulence where mixing and chemical reaction occur. The magnitudes of strain rate and vorticity increase as the eddy size (or wavelength) decreases in the turbulence energy cascade process. The Kolmogorov scale size is determined by the dissipation rate of turbulence kinetic energy and dynamic viscosity and is the smallest turbulence length scale. The final molecular mixing and chemical reaction in the combustion process occur on a still smaller scale, where there will be an axis (or direction) of principal compressive normal strain and an orthogonal axis for principal tensile strain; the third orthogonal axis could be either tensile or compressive. These axes would rotate under shear strain (or equivalently vorticity). Similarly, the direction of the scalar gradient rotates under shear. A useful flamelet model must have a statistically accurate representation of the relative orientations on this smallest scale of the vorticity vector, scalar gradients and the directions of the three principal axes for strain rate. Several studies exist that are helpful in understanding this important alignment issue.

Generally and always for incompressible flow, one principal strain rate $\gamma$ locally will be compressive (corresponding to inflow in a counterflow configuration), another principal strain rate $\alpha$ will be tensile (also named extensional and corresponding to outflow) and the third can be either extensional or compressive and will have an intermediate strain rate $\beta$ of lower magnitude than the other like strain rate. Specifically, $\alpha > \beta > \gamma, \;\alpha > 0 , \; \gamma <0$ and, for incompressible flow, $\alpha + \beta + \gamma = 0$. If the intermediate strain rate $\beta < 0$, there is inflow from two directions with outflow in one direction; a contracting jet flow occurs locally. Conversely, with $\beta > 0$, there is outflow in two directions and inflow in one direction; a counterflow or, in other words, the head-on collision of two opposed jets occurs. Betchov (1956) has shown that, for homogeneous, isotropic turbulence in an incompressible flow, the situation with $\beta > 0$ and resulting counterflow is the most important for production of vorticity and the turbulence energy cascade to smaller scales.

Several interesting findings result from direct numerical simulations for incompressible flows. Both Ashurst et al. (1987) and Nomura & Elghobashi (1992) compared a case of homogeneous sheared turbulence with a case of isotropic turbulence. They report that the vorticity alignment with the intermediate strain direction is most probable in both cases but especially in the case with shear. Furthermore, the intermediate strain rate is most likely to be extensive (positive) implying a counterflow configuration. Dresselhaus & Tabor (1991) use a kinematic approach to study the stretching of material and vorticity in a fluid flow and predict the tendency towards alignment of the intermediate strain direction with the vorticity. If the vorticity had strong alignment with the major compressive or major tensile strain direction, the magnitude of helicity, the dot product of velocity $\boldsymbol {u}$ and vorticity $\boldsymbol {\omega }$, would be large. Kerr (1987) reports that large values of helicity are not found in the turbulence cascade process. Ashurst et al. (1987) note further that the positive intermediate strain rate has a significantly smaller magnitude than either of the other two principal strain rates. Furthermore, the time for alignment of the vorticity with that intermediate direction is short compared to the eddy-turnover time.

Nomura & Elghobashi (1993) studied reacting flow and show that in regions of exothermic reaction and variable density, alignment of the vorticity with the most tensile strain direction can occur. Still though as the strain rates increase, the intermediate direction becomes more favoured for alignment with vorticity; that direction is also preferred in regions where mixing occurs without substantial divergence of the velocity due to chemical reaction.

We may also expect that a material interface most probably aligns to be normal to the direction of the compressive normal strain. That is, the scalar gradient and the direction of compressive strain are aligned; see Ashurst et al. (1987), Nomura & Elghobashi (1992), Nomura & Elghobashi (1993) and Boratav, Elghobashi & Zhong (1996, 1998). Authors agree that the most common intermittent vortex structures in regions of high strain rate are sheets or ribbons rather than tubes.

An important issue for flamelet modelling is the relative magnitudes of the vorticity and the rates of principal normal strain. For homogeneous, incompressible turbulence, Betchov (1956) showed that, for the average across all length scales, these quantities are of the same order of magnitude. Of course, the smallest scales contribute more to the average since the velocity derivatives are larger on those scales. Also, for shear flows, we expect that the turbulence at the smaller scales will be isotropic and behave more like the homogeneous flow. In our analysis for variable-density, reacting shear flows, we assume the same order-of-magnitude similarity between vorticity and the rates of principal normal strain applies for the smallest scales.

Based on those understandings concerning vector orientations, Sirignano (2021c) extended flamelet theory in a second significant aspect beyond the inclusion of both premixed and non-premixed flame structures; namely, a model was created of a three-dimensional field with both shear and normal strains. The three-dimensional problem is reduced to a two-dimensional form and then, for the counterflow or mixing-layer flow, to a one-dimensional similar form. The system of ODEs is presented for the thermochemical variables and the velocity components. Conserved scalars are determined and can become the independent variable if they behave in a monotonic fashion. Sirignano (2021c) also was able to use a velocity component as the independent variable for the flamelet model with shear strain and vorticity. The validity of the similar solution form for mixing layers with certain thin reaction zones was discussed using concepts from singular perturbation theory. The chemical-kinetic model appears as a source term for an ODE. These new findings are very helpful in improving the foundations for flamelet theory and its use in subgrid modelling for turbulent combustion. Still, however, there has not been a flamelet model that connects well the scales on the resolved LES level to the subgrid scales at the level of the flame structure.

Based on the observations of the needed improvements, the aim here is to develop a flamelet model that (i) determines rather than prescribes the existence of non-premixed flames, premixed flames or multi-branched flame structures; (ii) determines directly the effect of shear strain and vorticity on the flames; (iii) applies directly the resolved-scale strain rates and vorticity to the subgrid level without the use of a contrived progress variable; (iv) employs a three-dimensional flamelet model; and (v) considers the effect of variable density. Furthermore, some discussion and primitive analysis of the in situ function of the flamelet model in a turbulent shear flow is presented. The analysis uses one-step kinetics to avoid complications in this initial study; however, a clear template will exist for the employment of multi-step kinetics. The choice for the in situ study will be a mixing layer with the use of mixing-length theory. Here, the goal is not to advance the portion of the analysis for the resolved scale; rather, the coupling with the subgrid closure model will be made clear.

In § 2, the scaling and connections between the resolved shear flow and the flamelet behaviour on the subgrid scale are discussed. Section 3 has the description of a new subgrid flame model that better handles connection with strain and vorticity on the resolved scale, three-dimensional character and multi-branched flame structure. The application of the subgrid flamelet with a resolved shear flow is addressed in § 4. Concluding comments are made in § 5.

2. Scaling between flamelet scale and resolved scale

Current flamelet theory for use in LES or RANS makes no substantial attempt to scale properly between the small scale of the flamelet and the larger flow scales which pertain to the computational fluid dynamics analysis. A theoretical basis is needed to prescribe how to determine vorticity, strain rates and scalar gradients on the flamelet scale given those properties on the resolved scale. Certainly, as mentioned in § 1.2, a body of helpful literature exists on this subject; it is considered here. Specifically, we avoid the creation of arbitrary variables such as the FPV used in many publications. The FPV is actually not a good measure of progress. In order to relate increasing temperature during combustion to the flamelet theory, scalar dissipation rate is obliged to increase as the FPV increases. Instead, scalar dissipation rate on the flamelet scale should be related to temporally varying and spatially varying scalar and velocity gradients and not to temperature magnitude. The existing theory will not allow a hot gas to experience reductions in dissipation rates related to reductions in the values of velocity and scalar gradients which surely is not consistent with general flow patterns in turbulent combustors. Here, a new method is developed for determining burning rate from the flamelet theory and applying it to the resolved scale for LES or to the averaged flow field for RANS. During the burning process, temperature is able to increase even if strain rate and the associated scalar dissipation rate might be decreasing.

2.1. Scaling of velocity, strain rate and vorticity

Here, at first, we use approximate concepts which perhaps are reasonably well suited for a use of mixing-length theory to describe the time-averaged turbulent flow in a shear layer. The aim is to provide a simple framework for the first application and test of a new flamelet theory. More sophisticated and modern statistical approaches can be found (e.g. Pope 2000), and can be used in the future for examination of turbulent flows using RANS or LES.

The shear-driven flow on the larger scale can be characterized by a length $\delta$ and a time-averaged velocity difference $\Delta U$ across that particular length that ultimately relate to the magnitudes of the largest eddy size and the turbulence kinetic energy. Then, the velocity, length and time scales for the resolved scales are $\Delta U$, $\delta$ and $\delta /\Delta U$, respectively. The root-mean-square velocity fluctuation $u'$, the turbulence kinetic energy $k$ and the rate of dissipation of turbulence kinetic energy $\epsilon$ will have magnitudes of the order of $\Delta U$, $\Delta U^2$ and $\Delta U^3/\delta$, respectively. The Kolmogorov scale is the smallest scale in the turbulence energy cascade where the inertial and viscous effects balance each other (Pope 2000; White 2005). On that smallest scale, the characteristic velocity, length and time scales become $u_{\kappa } = (\nu \epsilon )^{1/4}$, $\kappa = (\nu ^3/\epsilon )^{1/4}$ and $t_{\kappa } = (\nu / \epsilon )^{1/2},$ respectively, where $\nu$ is the kinematic viscosity of the fluid. One can use these Kolmogorov scales to be the scales for flamelet analysis; a more refined approach might be to use the rate of dissipation for scalar quantities $\epsilon _c$ discussed by Elghobashi & Launder (1983) instead of $\epsilon$.

We estimate a Reynolds number for the resolved flow using $Re \equiv \Delta U \delta /\nu$. Vorticity, rate of normal strain and rate of shear strain will be $O(\Delta U/\delta )$. We estimate the magnitudes of rate of strain and vorticity on the flamelet scale to be given as

(2.1a,b)\begin{equation} S^* \equiv \frac{\partial u_{\kappa}}{\partial \kappa} = O\left(\frac{u_{\kappa}}{\kappa} \right); \quad \frac{u_{\kappa}}{\kappa} =\left(\frac{\epsilon}{\nu}\right)^{1/2} = \frac{\Delta U}{\delta} Re^{1/2}. \end{equation}

Clearly, for high $Re$ values, we may expect vorticity and rate of strain on the flamelet scale to be orders higher than found on the resolved scale. This scaling and the connection of the strain rates and vorticity on different scales have not been addressed in prior flamelet modelling.

Note that different estimates of a relevant resolved-scale or averaged-flow Reynolds number are used at later points in this discussion. The intention, however, is to maintain the same order of magnitude.

If we examine a time-averaged shear flow, the quantity $\Delta U/\delta$ can be replaced by the magnitude of a velocity gradient for the averaged flow treated through RANS simulations following a mixing-length concept. For LES, that quantity can be related to a velocity gradient on the smallest resolved scale, following an approach similar to the Smagorinsky model for Reynolds stress. Specifically,

(2.2a–c)$$\begin{gather} \frac{\Delta U}{\delta} \equiv \left|\frac{\partial u}{\partial x}\right| ; \quad Re \equiv \left|\frac{\partial u}{\partial x}\right| \frac{\delta^2}{\nu} ;\nonumber\\ \left|\frac{\partial u_{\kappa}}{\partial \kappa}\right| = \left|\frac{\partial u}{\partial x}\right|Re^{1/2} = \left|\frac{\partial u}{\partial x}\right|^{3/2}\frac{\delta}{\nu^{1/2}} = \frac{S_{rs}^{*3/2} \delta}{\nu^{1/2}}. \end{gather}$$

Here, one can make a choice about the interpretation of the resolved scale strain rate $S^{*}_{rs} \equiv |\partial u / \partial x|$. The largest component of shear strain rate on the resolved scale or for the time-averaged flow is recommended for use. Strain rate $S^{*}_{rs}$ can vary with location in the flow and, for unsteady RANS and LES, can vary with time as well. The quantity $\partial u_{\kappa }/\partial \kappa$ will be imposed as the compressive normal strain rate in the flamelet model; it will be a negative number $-(S^*_1 + S^*_2)$. Thus, (2.2a–c) can relate the normal strain rate on the flamelet scale to the strain rate on the resolved scale or for the time-averaged flow. In the case with the two-equation RANS model using $k, \epsilon$ theory, (2.1a,b) will yield $S^* = (\epsilon / \delta ^2)^{1/3}$.

A relation must be created between the dimensional vorticity $\omega ^*$ on the resolved scale and the dimensional vorticity $\omega ^*_{\kappa }$ on the flamelet subgrid scale. (Note that, in the development of the flamelet model, the subgrid $\omega _{\kappa }$ is dimensionless.) A reasonable relationship, mimicking the strain rate relation, is given as $\omega ^*_{\kappa } = \omega ^{*3/2}\delta / \nu ^{1/2} = O(\omega ^* Re^{1/2}) \gg \omega ^*.$

2.2. Scaling of the scalar properties

The inflow boundary conditions for the scalar quantities in the flamelet calculation must be determined from the resolved scale behaviour. The scalar gradients are much larger on the flamelet scale due to the dynamics of the turbulent flow. As a first approximation, consider that the scalar gradients scale in proportion to the strain rates. For example, using the mass fraction of species $m$, we state the rough approximation:

(2.3)\begin{equation} \frac{\dfrac{ \partial Y_{m, \kappa}} {\partial \kappa}}{\dfrac{ \partial Y_m} {\partial x} }= \frac{\dfrac{ \partial u_{\kappa}} {\partial \kappa} }{\dfrac{ \partial u} {\partial x} }= Re^{1/2}, \end{equation}

where the subscript $\kappa$ designates the flamelet scale. The domain sizes between the Kolmogorov scale and the resolved scale as $\kappa = \delta Re^{-3/4}$. Setting the change in scalar value across the given domain as the product of domain size and its gradient, the result is $\Delta Y_{m, \kappa } = \Delta Y_m Re^{-1/4}$. So, the variation in scalar properties across the smallest eddy is smaller than the variation across the larger eddies. Although the gradient of the scalar property is much greater on the flamelet scale due to turbulent mixing, the variation across the flamelet domain is smaller due to the greatly reduced domain size. Changes in enthalpy $h$ and density $\rho$ can be determined following the same pattern. The important implication is that, in general, the partial premixing on the smaller scales should be greater than experienced on the largest scales.

If, on the resolved scale at a particular point $\boldsymbol {x}, t$, the scalar property is $Y_m$, the bounding values for the inflow of the flamelet counterflow will be taken as $Y_{m, \kappa, \infty } = Y_{m}(\boldsymbol {x},t) + \Delta Y_{m, \kappa }$ and $Y_{m, \kappa, - \infty } = Y_{m}(\boldsymbol {x},t) - \Delta Y_{m, \kappa }$. Again, the boundary values for other scalars can be handled identically. A consequence here will be that the incoming streams of the flamelet counterflow are more likely to be fuel-rich or fuel-lean than pure fuel or pure oxidizer. Multi-branched flame structures of the type found by Hamins et al. (1985), Rajamanickam et al. (2019), Sirignano (2019, 2021c) and López-Cámara et al. (2020) can be expected.

2.3. Scaling of energy release rate and species consumption and production rates

The resolved scale will require input from the flamelet model for the quantities giving consumption (or production) rates per unit volume for the chemical species and energy release rate per unit volume due to chemical reaction and perhaps also viscous dissipation rate. The production and consumption rates within the flame are substantially higher than the average values over the counterflow volume. It is the average over the counterflow volume that should be used in the resolved-scale calculations. The same approach should be used for the viscous dissipation rate in high-speed flows where that is considered to be important. These quantities to be used on the resolved scale are given as integrals over the flamelet scale by (3.24).

The dimensional form of the rates indicates that the integrated chemical rates are proportional to the magnitude of the compressive strain rate $S^*\equiv S^*_1 + S^*_2$ for the counterflow. The portion of the dimensional energy release rate due to viscous dissipation is proportional to $S^{*2}$. That subgrid strain rate is much larger than the strain rate on the resolved scale as indicated by (2.3). The scalar gradients have similar scaling relation between the subgrid and the resolved scale. They impact diffusion and therefore, in diffusion-controlled combustion, they determine burning rates.

3. Subgrid flamelet analysis

We formulate the problem in a quasi-steady three-dimensional form. The following alignments are assumed. The direction of major compressive principal strain is orthogonal to the vorticity vector direction. Specifically, the intermediate principal strain direction is aligned with the vorticity while the scalar gradient aligns with the principal compressive strain direction. These assumptions are consistent with the statistical findings of Nomura & Elghobashi (1993).

3.1. Coordinate transformation

A transformation displayed in figure 1 is made from the Newtonian frame with rotating material (due to vorticity) to a rotating, non-Newtonian frame. Let the vorticity direction be the $z'$ direction in an orthogonal framework. Any $x', y'$ plane contains the directions of scalar gradients, major principal axis for compressive strain and major principal axis for tensile strain. Note that the $x',y', z'$ directions are not correlated with coordinates on the resolved scale. Here, $\omega _{\kappa }$ is the vorticity magnitude on this subgrid (Kolmogorov) scale. Directions $x', y', z'$ are transformed to $\xi, \chi, z'$ wherein the material rotation is removed from the $\xi, \chi$ plane by having it rotate at angular velocity ${\rm d}\theta /{\rm d}t = \omega _{\kappa }/2$ relative to $x', y'$. Here, $\theta$ is the angle between the $x'$ and $\xi$ axes and simultaneously the angle between the $y'$ and $\chi$ axes. Clearly, we take the subgrid domain to be sufficiently small to consider a uniform value of $\omega$ across it.

Figure 1. Transformation to $\xi, \chi, z'$ rotating coordinate system from $x', y', z'$ Newtonian system. Angle $\theta$ increases in the counterclockwise direction.

The following relations apply:

(3.1)\begin{equation} \left. \begin{aligned} \xi & = x' \cos \theta + y' \sin\theta ; \quad \chi = y' \cos \theta - x' \sin\theta;\\ \frac{\partial \xi}{\partial x'} & = \cos \theta ; \quad\frac{\partial \xi}{\partial y'}= \sin\theta ; \quad \frac{\partial \chi}{\partial x'} ={-} \sin \theta ; \quad \frac{\partial \chi}{\partial y'} = \cos \theta ;\\ u_{\xi} & = u \cos \theta + v \sin \theta + \chi\frac{\omega_{\kappa}}{2} ; \quad u_{\chi} = v\cos \theta -u \sin \theta - \xi\frac{\omega_{\kappa}}{2} ;\\ \frac{\partial u}{\partial x'} & = \frac{\partial u}{\partial \xi}\cos \theta - \frac{\partial u}{\partial \chi}\sin \theta ; \quad \frac{\partial u}{\partial y'} = \frac{\partial u}{\partial \xi}\sin \theta + \frac{\partial u}{\partial \chi}\cos \theta ;\\ \frac{\partial v}{\partial x'} & = \frac{\partial v}{\partial \xi}\cos \theta - \frac{\partial v}{\partial \chi}\sin \theta ; \quad \frac{\partial v}{\partial y'} = \frac{\partial v}{\partial \xi}\sin \theta + \frac{\partial v}{\partial \chi}\cos \theta . \end{aligned} \right\} \end{equation}

Since

(3.2)\begin{equation} \frac{\partial v}{\partial x'} - \frac{\partial u}{\partial y'} = \omega_{\kappa} , \end{equation}

it follows that

(3.3)\begin{equation} \frac{\partial u_{\chi}}{\partial \xi} - \frac{\partial u_{\xi}}{\partial \chi} = 0. \end{equation}

Namely, the flow in the rotating frame of reference does not have the vorticity imposed on it. However, two points must be understood. Firstly, the frame is not Newtonian and a reversed (centrifugal) force is imposed. Secondly, the expansions due to combustion and energy release can produce new vorticity but it will integrate to zero globally. The inflowing free-stream vorticity in the transformed coordinates is zero. Thus, in similar fashion to classical counterflow, the vorticity develops as an odd function so that the circulation on a contour surrounding the flow domain remains with zero value. This creation of vorticity is related to gas expansion with density variation; however, that expansion has a symmetry that yields an antisymmetry in the vorticity.

3.2. Governing equations

The governing equations for unsteady three-dimensional flow in the non-Newtonian frame can be written with $u_i = u_{\xi }, u_{\chi }, w ; \; x_i = \xi, \chi, z$. The centrifugal acceleration $a_i = \xi \omega _{\kappa }^2/4 , \chi \omega _{\kappa }^2 /4, 0$. The quantities $p, \rho, h, h_m, Y_m, \dot {\omega }, \mu, \lambda, D$ and $c_p$ are pressure, density, specific enthalpy, heat of formation of species $m$, mass fraction of species $m$, chemical reaction rate of species $m$, dynamic viscosity, thermal conductivity, mass diffusivity and specific heat, respectively. Furthermore, $\tau _{ij}$ is the viscous stress tensor and the Lewis number $Le \equiv \lambda /(\rho D c_p )$.

(3.4)\begin{gather} \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u_j)}{\partial x_j} =0, \end{gather}

(3.5)\begin{gather} \rho \frac{\partial u_i}{\partial t} + \rho u_j\frac{\partial u_i}{\partial x_j} +\frac{\partial p}{\partial x_i} = \frac{\partial \tau_{ij}}{\partial x_j} + \rho a_i, \end{gather}

where, following the Stokes hypothesis for a Newtonian fluid,

(3.6)\begin{gather} \tau_{ij} = \mu \left[ \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} -\frac{2}{3}\delta_{ij} \frac{\partial u_k}{\partial x_k}\right], \end{gather}

(3.7)\begin{gather} \rho \frac{\partial h}{\partial t} + \rho u_j\frac{\partial h}{\partial x_j} - \frac{\partial p}{\partial t} - u_j\frac{\partial p}{\partial x_j} = \frac{\partial}{\partial x_j} \left( \frac{\lambda}{c_p} \frac{\partial h}{\partial x_j} \right) \nonumber\\ + \frac{\partial}{\partial x_j} \left( \rho D (1 - Le) \varSigma^N_{m=1}h_m \frac{\partial Y_m}{\partial x_j} \right) - \rho \varSigma^N_{m=1}h_{f,m} \dot{\omega}_m + \tau_{ij}\frac{\partial u_i}{\partial x_j}, \end{gather}

(3.8)\begin{gather} \rho \frac{\partial Y_m}{\partial t} + \rho u_j\frac{\partial Y_m}{\partial x_j} = \frac{\partial }{\partial x_j}\left(\rho D \frac{\partial Y_m}{\partial x_j}\right) + \rho \dot{\omega}_m ;\quad m=1, 2,\ldots, N. \end{gather}

Equations (3.4)–(3.8) together with the equation of state and the relations describing fluid physicochemical properties give a complete description of behaviour in the non-Newtonian reference frame. These equations are used in the remainder of this section.

Equation (3.7) is a thermodynamic statement wherein, following a material element, we are stating the differential relation $\rho \,\textrm {d}h - \textrm {d}p = \rho T \,\textrm {d}s$. The terms on the right-hand side of (3.7) are entropy-producing terms. An alternative form of the energy equation can be developed to govern the total $H$ of the specific enthalpy, specific chemical energy and kinetic energy per unit mass. That is, $H \equiv h + \varSigma _{m=1}^NY_m h_{f,m} +u_k u_k/2$. Specifically, the vector dot product of $u_i$ with (3.5) is used to substitute for $u_j \partial p/ \partial x_j$ in (3.7) and (3.8) is used to substitute for $\dot {\omega }_m$ there. The Lewis number $Le =1$ is considered. It follows that

(3.9)\begin{equation} \rho \frac{\partial H}{\partial t} + \rho u_j\frac{\partial H}{\partial x_j} - \frac{\partial p}{\partial t} =\frac{\partial }{\partial x_j} \left( \frac{\lambda}{c_p} \frac{\partial (h+ \varSigma_{m=1}^NY_m h_{f,m}) }{\partial x_j} \right) +\frac{\partial (u_i \tau_{ij})}{\partial x_j}+\rho u_j a_j . \end{equation}

The energy source term $\rho u_j a_j = \rho (\omega _{\kappa }/2)^2(\xi u_{\xi } + \chi u_{\chi })$. If we neglect terms of the order of the kinetic energy per mass, this effect disappears. We might still retain the viscous dissipation rate $\tau _{ij} \partial u_i/\partial x_j$ for special cases where large strain rates are expected on the subgrid scale, as suggested by Drozda, Quinlan & Drummond (2020). The viscous dissipation rate will have exactly the same value whether it is calculated in the Newtonian frame or the rotating frame; this result is expected because it relates to the thermodynamics where the laws are independent of the reference frame. The resulting equation becomes

(3.10)$$\begin{gather} \rho \frac{\partial h}{\partial t} + \rho u_j\frac{\partial h}{\partial x_j} - \frac{\partial p}{\partial t} \approx \frac{\partial}{\partial x_j} \left( \frac{\lambda}{c_p} \frac{\partial h}{\partial x_j} \right) + \frac{\partial}{\partial x_j} \left( \rho D (1 - Le) \varSigma^N_{m=1}h_m \frac{\partial Y_m}{\partial x_j} \right) \nonumber\\ -\rho \varSigma^N_{m=1}h_{f,m} \dot{\omega}_m + \tau_{ij}\frac{\partial u_i}{\partial x_j}. \end{gather}$$

Here, we define the non-dimensional Prandtl, Schmidt and Lewis numbers: $Pr \equiv c_p \mu / \lambda$; $Sc \equiv \mu / (\rho D)$; and $Le \equiv Sc/ Pr$.

The non-dimensional forms of the above equations remain identical to the above forms if we choose certain reference values for normalization. In the remainder of this article, the non-dimensional forms of the above equations are considered. The superscript $^*$ is used to designate a dimensional property. The variables $u_i^*, t^*, x_i^*, \rho ^*, h^*, p^*$ and $\dot {\omega }_m^* ,$ and properties $\mu ^*, \lambda ^*/ c_p^*$ and $D^*$ are normalized respectively by $[(S_1^* +S_2^*) \mu _{\infty }^*/ \rho _{\infty }^*]^{1/2}, (S_1^* + S_2^*)^{-1}, [ \mu _{\infty }^*/ (\rho _{\infty }^*(S_1^* +S_2^*))]^{1/2}, \rho _{\infty }^*, (S_1^* +S_2^*) \mu _{\infty }^*/ \rho _{\infty }^*,$ $(S_1^* + S_2^*)\mu _{\infty }^*, (S_1^* + S_2^*), \mu _{\infty }^*, \mu _{\infty }^*$ and $\mu _{\infty }^*/ \rho _{\infty }^*$. The dimensional strain rates $S^*_1$ and $S^*_2$ and vorticity $\omega _{\kappa }^*$ are normalized by $S^*_1 + S^*_2$. It is understood that, for unsteady flow, the reference values for strain rates and far-stream variables and properties used for normalization are constants; for example, averages might be taken for fluctuating conditions. Note that the reference length $[ \mu _{\infty }^*/ (\rho _{\infty }^*(S_1^* +S_2^*))]^{1/2}$ is the estimate for the magnitude of the viscous-layer thickness. In the following flamelet analysis, the vorticity $\omega _{\kappa }$ and the velocity derivatives $\partial u_i/ \partial x_j$ are non-dimensional quantities; their dimensional values can be obtained through multiplication by $S^*_1 + S^*_2$. In § 2, the algorithms are given that relate dimensional vorticity and velocity derivatives on the resolved scale to dimensional vorticity and velocity derivatives on the subgrid scale.

3.3. Similar form for the velocity and pressure

The stagnation point in the steady counterflow is taken as the origin $\xi = \chi =z=0$. Along the line $\xi =z=0$ normal to the interface, we can expect the first derivatives of $u_{\chi }, \rho, h, T$ and $Y_m$ with respect to either $\xi$ or $z$ to be zero-valued. For unsteady cases, only symmetric situations are considered so that the stagnation point remains at the origin and the interface remains at $\chi =0$. The velocity components $u_{\xi }$ and $w$ are odd functions of $\xi$ and $z$, respectively, going through zero and changing sign at that line. Consequently, upon neglect of terms of $O(\xi ^2)$ and $O(z^2)$, the variables $u_{\chi }, \rho, h, T$ and $Y_m$ can be considered to be functions only of $t$ and $\chi$. For steady flow, the density-weighted Illingworth (1949) transformation of $\chi$ can be used to replace $\chi$ with $\eta \equiv \int ^{\chi }_0 \rho (\chi ') \,\textrm {d}\chi '$. Neglect of the same order of terms implies that $u_{\xi }= S_1 \xi (\textrm {d}f_1/\textrm {d}\eta )$ and $w= S_2 z(\textrm {d}f_2/\textrm {d}\eta )$. Note that $u_{\xi }$ is independent of $z$ and $w$ is independent of $\xi$ in this case where no shear strain is imposed on the incoming stream(s). At the edge of the viscous layer at large positive $\eta$, $\textrm {d}f_1/\textrm {d} \eta \rightarrow 1, \textrm {d}f_2/\textrm {d} \eta \rightarrow 1, f_1 \rightarrow \eta$ and $f_2 \rightarrow \eta$. Ordinary differential equations are created here through the variable $\eta$ and the convenient definition is made that $( )' \equiv \textrm {d}( )/\textrm {d}\eta$. Note that other transformations of the $\chi$ coordinate can be made, e.g. weighting by transport properties (Linan et al. 2017; Weiss et al. 2018) rather than density.

In the non-dimensional form given by (3.4)–(3.10), the dimensional strain rates $S_1^*$ and $S_2^*$ are each normalized by the dimensional sum $S_1^* + S_2^*$. Thus, the non-dimensional relation is $S_2 = 1 - S_1$ and only one independent non-dimensional strain-rate parameter is needed. Nevertheless, two strain rates are presented above and in the following analysis with the understanding that one depends on the other such that $S_1 + S_2 = 1$. Values of $S_1 + S_2$ are explicitly stated in our analysis without substitution of the unity value. This choice clarifies whether a particular term when converted to a dimensional form depends on $S_1^*, S_2^*$ or the sum of the two strain rates.

For steady state, the continuity equation (3.4) is readily integrated to give

(3.11)\begin{equation} \rho u_{\chi} ={-}S_1 f_1(\eta) - S_2 f_2(\eta). \end{equation}

For derivatives of variables that are functions of $\eta$ only, the prime is used to indicate differentiation with respect to $\eta$. Then,

(3.12)\begin{equation} u'_{\chi} = \frac{S_1 f_1(\eta) + S_2 f_2(\eta)}{\rho^2 } \rho' -\frac{S_1 f_1'(\eta)+ S_2 f_2'(\eta)}{\rho}. \end{equation}

Thus, the incoming inviscid flow outside the boundary layer is described by $u_{\chi }=-(S_1 + S_2)\eta$ for positive $\eta$ and $u_{\chi }=-(S_1 + S_2)\eta / \rho _{-\infty }$ for negative $\eta$. Note that the same result is found for the unsteady or steady incompressible state where there is no need to use $\eta$ in place of $\chi$ since $\rho =1$ everywhere. Then, $u_{\chi } = -(S_1 + S_2 ) \chi$ for the external incoming flow.

Equations (3.11) and (3.12) may be substituted into (3.5) to determine the pressure gradient:

(3.13)\begin{equation} \left. \begin{aligned} \frac{\partial p}{\partial \xi} & = \rho [ \rho \mu S_1 f_1''' + S_1 f_1'' (\rho \mu)' +(S_1 f_1 + S_2f_2)S_1 f_1'' + \left(\frac{\omega_{\kappa}}{2} \right)^2 - (S_1f_1')^2]\xi ,\\ \frac{\partial p}{\partial \eta} & = \frac{4}{3} ( \rho \mu u_{\chi}')' \!-\!\frac{2}{3}(S_1f_1' \!+\! S_2 f_2')\mu' \!+\!\frac{ \mu}{3}(S_1 f_1'' \!+\! S_2 f_2'') + \left(\frac{\omega_{\kappa}}{2} \right)^2\chi +(S_1 f_1 + S_2 f_2)u_{\chi}' ,\\ \frac{\partial p}{\partial z} & = \rho [ \rho \mu S_2 f_2''' + S_2 f_2'' (\rho \mu)' +(S_1 f_1 + S_2f_2)S_2 f_2''- (S_2f_2')^2]z . \end{aligned} \right\} \end{equation}

It follows from the $\eta$ pressure gradient in (3.13) that $\partial p/ \partial \eta$ is a function only of $\eta$. Therefore, $\partial ^2 p/ \partial \xi \partial \eta = 0$ and $\partial ^2 p/ \partial z \partial \eta = 0$. Now, the coefficient of $\xi$ on the right-hand side of the $\xi$ pressure gradient in (3.13) must be constant. The same conclusion is made for the coefficient of $z$ on the right-hand side of the $z$ pressure gradient in (3.13). At $\eta =\infty$, $f_1' =f_2'=1$ and $f_1'' = f_2'' = f_1''' = f_2''' =0$ which allows the two constants to be determined. Specifically, we obtain

(3.14)\begin{equation} \left. \begin{gathered} \rho \mu f_1''' + f_1'' (\rho \mu)'+(S_1 f_1 + S_2f_2) f_1'' + S_1\left(\frac{1}{\rho} -(f_1')^2\right) +\frac{\omega_{\kappa}^2}{4S_1}\left(1 - \frac{1}{\rho}\right)= 0 ,\\ \rho \mu f_2''' + f_2'' (\rho \mu)'+(S_1 f_1 + S_2f_2) f_2'' + S_2\left(\frac{1}{\rho}-(f_2')^2\right) = 0. \end{gathered} \right\} \end{equation}

The boundary conditions will use the assumption that two velocity components asymptote to the constant values $u_{\xi }(\infty ), u_{\xi }(-\infty ), w(\infty )$ and $w(-\infty )$ at large magnitudes of $\eta$. The stream function bounding the two incoming streams is arbitrarily given a zero value and placed at $\eta =0$.

(3.15)\begin{equation} \left. \begin{gathered} f'_1(\infty) = 1 ; \quad f'_1(-\infty)= \sqrt{\frac{1}{\rho(-\infty)}+ \left(\frac{\omega_{\kappa}}{2 S_1}\right)^2 \left( 1 - \frac{1}{\rho(-\infty)} \right) } ; \quad f_1(0) = 0 ;\\ f'_2(\infty) = 1 ; \quad f'_2(-\infty)= \frac{1}{\sqrt{\rho(-\infty)}}; \quad f_2(0) =0. \end{gathered} \right\} \end{equation}

In the incompressible case where density is uniform throughout the flow, i.e. $\rho =1$, the solutions become simply that $f'_1(\eta )=1$ and $f'_2(\eta ) =1$ everywhere. When density varies through the flow because of heating or variation of composition, $u_{\xi }$ and $w$ vary with $\chi$, thereby creating a shear stress and vorticity albeit that the frame transformation removed vorticity and shear from the incoming flow. The vorticity will be created in an antisymmetric manner since the two velocity components are odd functions of $\xi$ and $z$, respectively. Thereby, the circulation remains zero for the transformed counterflow.

For steady flows, $S_1 + S_2 =1$. The dependence of $u_{\chi }$ on $f\equiv S_1f_1 + S_2 f_2$ is shown by (3.11). Thus, the function $f$ will be important in determining both the field for $u_{\chi }$ and the scalar fields. From (3.14), an equation for $f$ can be formed:

(3.16)\begin{equation} \rho \mu f''' + f'' (\rho \mu)'+f f'' + \frac{1}{\rho} -(f')^2 +\frac{\omega_{\kappa}^2}{4}\left(1 - \frac{1}{\rho}\right)= 2S_1S_2\left(\frac{1}{\rho} -f_1'f_2'\right). \end{equation}

Consequently, $f$ as well as $f_1$ and $f_2$ depend on both $S_1$ and $S_2$, not merely on $S_1 + S_2$. That is, the particular distribution of the normal strain rate between the two transverse direction matters. Functions $f$ and $f_1$ also depend directly on $\omega _{\kappa }$ (unless $S_1 = 0$). Function $f_2$ depends on $\omega _{\kappa }$ through its coupling with $f_1$. In our calculations, we consider a planar case ($S_1 = 1.0 , S_2 = 0$) where the product $S_1S_2$ is minimized and the vorticity vector is normal to the plane with strain. The case where $S_1 = S_2 = 0.5$ has the maximum value for the product $S_1 S_2$ would be axisymmetric if $\omega = 0$. However, symmetry is lost if the rotation exists.

In the outer flow where variability of viscosity and density may be neglected and $u_{\chi } \rightarrow \eta / \rho _{\infty }$ as $\eta \rightarrow \infty$, the $\chi$-momentum equation from (3.13) becomes

(3.17)\begin{equation} \left. \begin{aligned} \frac{\partial }{\partial \eta} \left[p - \frac{4\rho \mu}{3}u_{\chi}' \right] & \rightarrow (S_1 f_1 + S_2 f_2)u_{\chi}' + \left( \frac{\omega_{\kappa}}{2}\right)^2 \chi \rightarrow (S_1 + S_2 )\eta u_{\chi}' + \left(\frac{\omega_{\kappa}}{2}\right)^2 \frac{\eta}{\rho_{\infty}} ,\\ \frac{\partial p}{\partial \eta} & \rightarrow \left[\left(\frac{\omega_{\kappa}}{2}\right)^2 -(S_1 + S_2)^2 \right] \frac{\eta}{\rho_{\infty}}. \end{aligned} \right\} \end{equation}

The negative pressure gradient will be reduced by the centrifugal effect. Essentially, the pressure gradient serves to decelerate the incoming stream in a counterflow; here, it is helped by the imposed acceleration.

The variable-density-and-viscosity case requires some couplings with (3.7) and (3.8) and with an equation of state and fluid-property laws which affect $\rho$ and $\mu$.

The pressure derivative can be determined by substituting from (3.11) and (3.12) into (3.13):

(3.18)\begin{align} \frac{\partial p}{\partial \eta} &=\frac{4}{3}\left( \rho \mu\left[ \frac{S_1 f_1(\eta) + S_2 f_2(\eta)}{\rho^2 } \rho'-\frac{S_1 f_1'(\eta)+ S_2 f_2'(\eta)}{\rho} \right]\right)' + \left(\frac{\omega_{\kappa}}{2} \right)^2\chi \nonumber\\ &\quad +\frac{ \mu}{3}(S_1 f_1'' + S_2 f_2'')-\frac{2}{3}(S_1f_1' + S_2 f_2')\mu' \nonumber\\ &\quad +(S_1 f_1 + S_2 f_2)\left[\frac{S_1 f_1(\eta) + S_2 f_2(\eta)}{\rho^2 } \rho'-\frac{S_1 f_1'(\eta)+ S_2 f_2'(\eta)}{\rho}\right] , \nonumber\\ p(x, \eta, z) &= p_{ref} + \frac{4\mu(\eta)}{3\rho(\eta)} \left[ (S_1 f_1(\eta) + S_2 f_2(\eta)) \rho' -\rho(\eta)(S_1 f_1'(\eta)+ S_2 f_2'(\eta))\right] \nonumber\\ &\quad + \frac{4 \mu(0)}{3}(S_1f_1'(0) + S_2f_2'(0)) \nonumber\\ &\quad + \int_{0}^{\eta}\left[\frac{\mu(\zeta)}{3}(S_1 f_1''(\zeta) + S_2 f_2''(\zeta)) -\frac{2}{3}(S_1f_1'(\zeta) + S_2 f_2'(\zeta))\mu' \right] \,{\rm d}\zeta \nonumber\\ &\quad +\int_{0}^{\eta}(S_1 f_1(\zeta) + S_2 f_2(\zeta))\left[\frac{S_1 f_1(\zeta) + S_2 f_2(\zeta)}{\rho(\zeta)^2 } \rho' -\frac{S_1 f_1'(\zeta)+ S_2 f_2'(\zeta)}{\rho(\zeta)}\right] \,{\rm d}\zeta \nonumber\\ &\quad- \frac{S_1^2 x^2}{2} - \frac{S_2^2 z^2}{2} + \left(\frac{\omega_{\kappa}}{2} \right)^2 \int_0^{\eta} \chi(\zeta) \,{\rm d}\zeta . \end{align}

The viscous dissipation $\varPhi$ can be determined as follows:

(3.19)\begin{align} \varPhi & \equiv \tau_{ij}\frac{\partial u_i}{\partial x_j} \!=\! \mu \left[\frac{4}{3} \left( \left(\frac{\partial u_{\xi}}{\partial \xi}\right)^2 + \left(\frac{\partial u_{\chi}}{\partial \chi}\right)^2 + \left(\frac{\partial w}{\partial z}\right)^2 - \frac{\partial u_{\xi}}{\partial \xi}\frac{\partial u_{\chi}}{\partial \chi} - \frac{\partial w}{\partial z}\frac{\partial u_{\chi}}{\partial \chi} - \frac{\partial w}{\partial z}\frac{\partial u_{\xi}}{\partial \xi} \right) \right. \nonumber\\ &\quad \left.+ \left(\frac{\partial u_{\xi}}{\partial \chi}\right)^2 + \left(\frac{\partial w}{\partial \chi}\right)^2 \right] \nonumber\\ &= \mu\left( 4(S_1f'_1)^2 + 4(S_2f'_2)^2 + \frac{8}{3}S_1 S_2 f'_1 f'_2 - 4\frac{(S_1f_1 + S_2f_2)(S_1f'_1 + S_2f'_2)}{\rho}\rho' \right.\nonumber\\ &\left.\quad+ \frac{4}{3}\left(\frac{S_1f_1 + S_2f_2}{\rho} \rho' \right)^2 \right) . \end{align}

Here, neglect is made of the quantities $\mu (\partial u_{\xi }/\partial \chi )^2$ and $\mu (\partial w/\partial \chi )^2$ because they are of $O(\xi ^2)$ and $O(z^2)$, respectively.

An exact solution of the variable-density Navier–Stokes equation has been obtained subject to determination of $\rho$ and $\mu$ through solutions of the energy and species equations as discussed below. There has been no need for use of a boundary-layer approximation. Thus, the solution here is the natural solution, subject to neglect of terms of $O(\xi ^2)$ and $O(z^2)$. Unlike the incompressible counterflow, a viscous layer exists with the three normal strains and normal viscous stresses varying through the layer due to varying density and viscosity. Shear strain does not appear explicitly in the incoming flow because it has been removed through the coordinate transformation to the rotating frame. Nevertheless, its effect appears through the rotational term.

3.4. Similar form for the scalar fields

Consider the reacting, steady case keeping viscous dissipation. Assume Prandtl number $Pr$ is constant. Substitution from (3.7) and (3.8) yields the following ODEs:

(3.20)$$\begin{gather} ( \rho \mu h' )' + Pr f h' + ( \rho^2 D (Pr - Sc) \varSigma^N_{m=1}h_m Y'_m )' = Pr \varSigma^N_{m=1}h_{f,m} \dot{\omega}_m - Pr \frac{\varPhi}{\rho} , \nonumber\\ (\rho^2 D Y'_m )' + f Y'_m ={-} \dot{\omega}_m ;\quad m=1, 2,\ldots, N . \end{gather}$$

If $Le =1$, i.e. $Pr = Sc$, the new scalar $\tilde {h} \equiv h + \varSigma _{m=1}^N h_{f,m} Y_m$ is governed by

(3.21)\begin{equation} (\rho \mu \tilde{h}')' + Pr f \tilde{h}'={-} Pr\frac{\varPhi}{\rho}. \end{equation}

When the viscous dissipation is negligible, $\tilde {h}$ is a conserved scalar indicating that the sum of thermal energy plus chemical energy is in an advective–diffusive balance. It remains to use thermodynamic relations to substitute for $\rho$ and $\mu$ in terms of $h$ and $p$.

Now, we address the special case where $\rho \mu =1$. The perfect gas law and the assumption of constant specific heat $c_p$ will give the relation that $1/ \rho = h$. Then, (3.14), (3.15), (3.20) and (3.21) yield

(3.22)\begin{equation} \left. \begin{gathered} f_1''' +f f_1'' + S_1[ h -(f_1')^2] + \frac{\omega_{\kappa}^2}{4S_1}(1 - h) =0, \\ f_2''' +f f_2'' + S_2[h-(f_2')^2] = 0,\\ Y''_m + Pr f Y'_m ={-} Pr\dot{\omega}_m ;\quad m=1, 2,\ldots, N, \\ h'' + Pr f h' + \left(\frac{1}{Le} -1\right) \varSigma^N_{m=1}h_m Y''_m = Pr \varSigma^N_{m=1}h_{f,m} \dot{\omega}_m - Pr \frac{\varPhi}{\rho}, \\ \tilde{h}'' +Pr f \tilde{h}' ={-}Pr \frac{\varPhi}{\rho}. \end{gathered} \right\} \end{equation}

The boundary conditions are

(3.23)\begin{equation} \left. \begin{gathered} f_1'(\infty)= 1 ;\quad f_1'(-\infty) = \sqrt{h_{-\infty}+ \left(\frac{\omega_{\kappa}}{2 S_1}\right)^2 ( 1 - h_{-\infty}) } ; \quad f_1(0) =0 ; \\ f_2' (\infty)=1 ;\quad f_2'(-\infty) = \sqrt{h_{-\infty}} ; \quad f_2(0) =0 ;\\ h(\infty) = 1 ; \quad h(-\infty)= \frac{1}{\rho_{-\infty}}; \\ Y_m (\infty) = Y_{m,\infty} ; \quad Y_m (-\infty) = Y_{m,-\infty} ;\\ \tilde{h}(\infty) = 1 + \varSigma_{m=1}^N h_{f,m}Y_{m, \infty} ; \quad \tilde{h}(-\infty) = h_{-\infty}+ \varSigma_{m=1}^N h_{f,m}Y_{m, -\infty}. \end{gathered} \right\} \end{equation}

Equations (3.14) indicate a dependence of the heat and mass transport on $f\equiv S_1f_1 + S_2f_2$. Manipulation of the first two equations of (3.22) leads to an ODE for $f$ with $S_1 S_2$ and $S_1S_2f_1'f_2'$ as parameters, clearly indicating that generally $f$ will have a dependence on $S_1S_2$. Thus, the behaviour for the counterflow can vary from the planar value of $S_1 =1, S_2 =0$ (or vice versa) or from the case $S_1 = S_2 =1/2$. This clearly shows that distinctions must be made amongst the various possibilities for three-dimensional strain fields as $S_1S_2$ varies between large negative numbers and $1/4$. An exception will be the incompressible case with constant properties where the $S_1S_2$ terms cancel in the equation for $f$.

The vorticity $\omega _{\kappa }$ impacts directly $f_1$ and $f$; thereby, it is affecting the velocity field. Then, through the advection of the scalar properties, there is impact on mass fractions and enthalpy. If the vorticity $\omega _{\kappa } =0$, a simple inspection of the governing ODEs leads to the conclusion that the values for $f_1, f'_1, f_2, f'_2,u/x$ and $w/z$ can be interchanged with the values for $f_2, f'_2, f_1, f'_1,w/z$ and $u/x$, respectively, when $S_1$ and $S_2$ are replaced by $1-S_1$ and $1-S_2$, respectively.

Note that, for $S_1 >1$ or $S_2 >1$ (which imply $S_2 < 0$ or $S_1 <0$, respectively), there would be incoming streams from two directions. One incoming stream would have a prescribed velocity profile in the viscous layer determined as a local exact solution to the Navier–Stokes conditional on its matching the profile determined by upstream conditions for the flow in that direction; this situation is too highly contrived and is not considered here. Thus, $S_1$ and $S_2$ are always each non-negative and bounded above by unity value in our considerations here. The figures show results for three strain rates: $S_1 = 0$ (planar case); $S_1 =0.25$ (three-dimensional strain); and $S_1 = 0.5$ (axisymmetric case).

If the viscous dissipation is negligible, $\tilde {h}$ becomes a conserved scalar. Other conserved scalars can be created. Mixture fraction is a popular choice. For the case of one-step kinetics with $Le =1$ and negligible viscous dissipation, the Shvab–Zel’dovich variables, $\alpha \equiv Y_F - \nu _S Y_O$ and $\beta \equiv h + \nu _S Y_O Q$, become conserved scalars. Parameters $Y_F, Y_O, \nu _S$ and $Q$ are fuel mass fraction, oxygen mass fraction, stoichiometric mass ratio and chemical energy per unit mass of fuel, respectively. For steady-state and time-averaged flows, these conserved scalars vary monotonically across a flow field and can be used to replace one of the spatial coordinates. This coordinate transformation in the counterflow configuration results in a new form of the scalar equations where the advective term is not present; a reactive–diffusive balance results. This result has classically been used (Peters 2000; Pierce & Moin 2004) together with an incompressible-flow assumption which gives an overly simplistic solution for the velocity field. Sirignano (2019, 2021a,c) has shown the transformation gains little when the variable density is considered; furthermore, the variable density significantly affects the reacting counterflow.

Consider the production or consumption rate of a particular species over the counterflow volume. We can integrate over a volume either using the original form in (3.8) or, more conveniently, using (3.22) to get exactly the same result. Consider the volume $-a < \xi < a, -b < y < b, -c < z < c$. The choices of lengths $a$ and $c$ will not matter on a per-unit-volume basis since mass fraction $Y_m$ and reaction rate $\dot {\omega }_m$ do not vary with $x$ or $z$. Length $c$ is chosen to be of the order of the Kolmogorov scale. Volume $V = 8abc$, $\tilde {\varPhi }$ is the volume-averaged viscous dissipation rate and $\widetilde {\rho \dot {\omega }_m}$ is the average mass production rate over the volume. It follows from integration of (3.20) after multiplication by density $\rho$ and division by $Pr V$ that

(3.24)\begin{equation} \left. \begin{gathered} \int^a_{{-}a}\int^b_{{-}b}\int^c_{{-}c} \frac{\rho}{Pr V} [Y''_m + Pr f Y'_m + Pr\dot{\omega}_m]\,{{\rm d} x}\,{{\rm d} y}\,{\rm d}z = 0;\quad m=1, 2,\ldots, N ,\\ \widetilde{\rho\dot{\omega}_m} \equiv \frac{1}{V}\int_V \rho \dot{\omega}_m \,{\rm d}V={-} \frac{1}{2b} \int^{\eta(b)}_{\eta({-}b)}fY'_m \,{\rm d}\eta ;\quad m=1, 2,\ldots, N ,\\ \int^a_{{-}a}\int^b_{{-}b}\int^c_{{-}c} \frac{\rho}{Pr V} \left[ h'' + Pr f h' - Pr \varSigma^N_{m=1}h_{f,m} \dot{\omega}_m + Pr\frac{\varPhi}{\rho}\right] \,{{\rm d}x} \,{{\rm d} y} \,{\rm d}z =0 ,\\ \varSigma^N_{m=1}h_{f,m} \widetilde{\rho \dot{\omega}_m }- \tilde{\varPhi} = \frac{1}{2b}\int^{\eta(b)}_{\eta({-}b)} fh' \,{\rm d}\eta. \end{gathered} \right\} \end{equation}

The quantity $\widetilde {\rho \dot {\omega }_m }$ is the only value that need be transferred from the subgrid flamelet analysis to the resolved-scale simulation. Note that the reaction rate is invariant under the change of reference frame. Here, $b$ has been considered large enough so that $Y'_m =0$ and $h'=0$ at those boundaries are good approximations. However, the value for $\widetilde {\rho \dot {\omega }_m }$ depends on the chosen domain size $2b$, which has a value of $O(10)$ typically in our analysis. Future studies should examine the optimization of the domain size $b$. Also, $Le =1$ has been considered. The volume-averaged viscous dissipation rate $\tilde {\varPhi }$ may be obtained by integration of (3.19).

Consider a species $m$ that is flowing inward away from $\eta = \infty$ towards $\eta = 0$. If it is being produced (consumed), the derivative $Y'_m$ in (3.24) will be negative (positive) for $\eta > 0$ where velocity $v<0$ and $f>0$. The signs are opposite for a species flowing inward away from $\eta = -\infty$ and towards $\eta = 0$. The equation provides two ways to evaluate the average production (consumption) rate for species $m$. The volume integral of the reaction rate has highly nonuniform integrand values over the space while the outflow integral over $\eta$ has a smoother variation of the integrand.

The flamelet model requires inputs that are scaled from the resolved flow. Specifically, rate of strain and vorticity, pressure and the inflowing scalar values for the counterflow are needed. The magnitude of the resolved-scale velocity is not relevant because the subgrid velocities are measured relative to a frame moving with a Galilean transformation. The flamelet can give back to the resolved flow the instantaneous value for energy release rate.

3.5. Chemical kinetics model

The above equations can be readily applied for diffusion flame counterflows and partially premixed flame counterflows, as explained in the following sections. They can also describe situations where multi-branched flames exist. For the case of zero vorticity, i.e. $\omega =0$, the generality has been shown by Sirignano (2019, 2021a). Although the analysis allows for the use of detailed chemical kinetics in the future, we focus here on propane–oxygen flows with one-step kinetics. However, results are expected to be qualitatively more general, applying to situations with more detailed kinetics and to other hydrocarbon–oxygen or hydrocarbon–air combinations. Westbrook & Dryer (1984) kinetics are used; they were developed for premixed flames but any error for non-premixed flames is viewed as tolerable often here because diffusion would generally be rate-controlling. Using asterisks to denote dimensional quantities, one may deduce that

(3.25)\begin{equation} \omega^*_F ={-} A^* {\rho^*}^{0.75} Y_F^{0.1} Y_O^{1.65} e^{{-}50.237/\tilde{h}}, \end{equation}

where the ambient temperature is set at 300 K and density $\rho ^*$ is given in units of kilograms per cubic metre. Here, $A^* =4.79\times 10^8\, (\textrm {kg}\,\textrm {m}^{-3})^{-0.75}\,\textrm {s}^{-1}$. The dimensional strain rate $S^*_1 + S^*_2$ (at the subgrid scale) is used to normalize time and reaction rate. In non-dimensional terms,

(3.26)\begin{equation} \left. \begin{aligned} \dot{\omega}_F & ={-} \frac{A^* {\rho_{\infty}^*}^{0.75}}{S_1^* + S^*_2}\tilde{h}^{{-}0.75} Y_F^{0.1} Y_O^{1.65} e^{{-}50.237/\tilde{h}},\\ \dot{\omega}_F & ={-} \frac{Da}{\tilde{h}^{0.75}} Y_F^{0.1} Y_O^{1.65} e^{{-}50.237/\tilde{h}}. \end{aligned} \right\} \end{equation}

The above equation defines the Damköhler number $Da$. Furthermore, we set $Da \equiv K Da_{ref}$, where

(3.27a,b)\begin{equation} Da_{ref} \equiv \frac{ \tilde{A}(10 \,{\rm kg}\,{\rm m}^{{-}3})^{0.75}}{(10^4\,{\rm s}^{{-}1})} =2.693 \times 10^5 ; \quad K \equiv \left[\frac{\rho_{\infty}^*}{10 \,{\rm kg}\,{\rm m}^{{-}3}}\right]^{0.75}\frac{10^4\,{\rm s}^{{-}1}}{S_1^* + S_2^*}. \end{equation}

Values of $10\, \textrm {kg}\,\textrm {m}^{-3}$ and $10\,000\,\textrm {s}^{-1}$ are conveniently chosen as reference values for density and strain rate, respectively. The reference value for density implies an elevated pressure. The strain-rate reference value is in the middle of an interesting range for this chemical reaction. Suppose the resolved scale has a velocity, a mixing length and kinematic viscosity with the following orders of magnitude, respectively: $10\,\textrm {m}\,\textrm {s}^{-1}$, $0.1\,\textrm {m}$ and $10^{-4}\,\textrm {m}^2\,\textrm {s}^{-1}$. Then, we may estimate $Re =10^4$ and resolved-scale strain rate $\partial u /\partial x = 100\,\textrm {s}^{-1}$. Thereby, (2.3) yields that $\partial u_{\kappa }/\partial \kappa = Re^{1/2}\partial u / \partial x = S^*_1 +S^*_2 = 10^4\,\textrm {s}^{-1}$.

Clearly, there is no need to set pressure (or its proxy, density) and the strain rate separately for a one-step reaction. For propane and oxygen, the mass stoichiometric ratio $\nu =0.275$.

The non-dimensional parameter $K$ increases (decreases) as the strain rate decreases (increases) and/or the pressure increases (decreases). Our reference case is $K =1$ and the range covered includes $O(10^{-1}) \leq K\leq O(10^2)$, allowing for the needed variation in strain rate and pressure to sustain premixed flamelets, diffusion flamelets and multi-branched flamelets.

The system of ODEs is solved numerically using a relaxation method and central differences. Solution over the range $-5 \leq \eta \leq 5$ provides adequate fittings to the asymptotic behaviours. The parameters that are varied are $K, Pr$ and $S_1$ (and thereby $S_2 = 1 -S_1$). Most calculations have $Pr = 1$ with emphasis on the effect of variation in $K$, i.e. pressure and strain rate.

3.6. Uncoupled diffusion flamelet calculations

Now, we treat a situation with a three-dimensional diffusion flame structure at the subgrid level. Figures 2 and 3 show the influence of vorticity on the flamelet stability near the extinction limit. The rotation of the flamelet due to vorticity causes a centrifugal effect on the counterflow velocity and thereby on the residence time in the vicinity of the reaction zone. The dashed red and solid red curves represent $K$ values of 0.195 and 0.196, respectively, both without vorticity. Despite the very modest difference in $K$, one survives as a strong flame and the other is basically extinguished. By applying $K=0.195$ and $\omega _{\kappa } =1.0$ and therefore rotational speed of value $\textrm {d}\theta / \textrm {d}t =0.5$, the solid blue curve is obtained, indicating a strong flame with regard to both reaction rate and peak temperature or enthalpy is induced by the rotation. Actually, the $K=0.196, \,\omega _{\kappa }= 1.0$ curve falls right on the solid blue curve as well, except in figure 2(d) where it appears as a dashed purple curve. The rotation causes a lower mass flux $f$ than was obtained without rotation. The inflow rate is diminished because the centrifugal effect creates a more adverse pressure gradient. However, the slower flow rate allows a longer residence time and modifies the extinction limit. As $K$ is lowered to values of 0.180 or 0.185, the flame is essentially extinguished in spite of the increased residence time due to the rotational effects.

Figure 2. Diffusion flame with varying vorticity and Damköhler number, for $S_1 =0.75 ; S_2 = 0.25$. Dashed red, $K= 0.195, \omega _{\kappa } = 0$; solid red, $K= 0.196, \omega _{\kappa } = 0$; solid blue, $K= 0.196, \omega _{\kappa } = 1.0$; solid purple, $K= 0.185, \omega _{\kappa } = 1.0$; dashed blue, $K= 0.180, \omega _{\kappa } = 1.0$. The $K=0.195, \omega _{\kappa }=1.0$ case is covered by the solid blue curve, except in (d) where it is given by the dashed purple curve. (a) Enthalpy, $h/h_{\infty }$. (b) Fuel mass fraction, $Y_F$. (c) Product of mass ratio and oxygen mass fraction, $\nu Y_O$. (d) Integral of reaction rate, $\int \dot {\omega }_F \,\textrm {d} \eta$. (e) Reaction rate, $\dot {\omega }_F$.

Figure 3. Diffusion flame with varying vorticity and Damköhler number, for $S_1 =0.75 ; S_2 = 0.25$. Dashed red, $K= 0.195, \omega _{\kappa } = 0$; solid red, $K= 0.196, \omega _{\kappa } = 0$; solid blue, $K= 0.196, \omega _{\kappa } = 1.0$; solid purple, $K= 0.185, \omega _{\kappa } = 1.0$; dashed blue, $K= 0.180, \omega _{\kappa } = 1.0$. The $K=0.195, \omega _{\kappa }= 1.0$ case is covered by the solid blue curve. (a) Mass flux per area, $f= - \rho u_{\chi }$. (b) Velocity component, $f_1' = u_{\xi }/(S_1 \xi )$. (c) Velocity component, $f_2' = w/(S_2z)$. (d) Velocity component, $u_{\chi }$.

The rotation also causes a decrease in $f_1'$ and therefore in the $\xi$ component of velocity as shown in figure 3. There is an associated increase in $f_2'$ and therefore an increase in the $z$ component of velocity. Basically, the low-density products of combustion more easily flow in the $z$ direction wherein no centrifugal effect is applied. This gives some physical explanation to the findings of Nomura & Elghobashi (1993) where, for reacting flows but not for non-reacting flows, the major extensional axis tends to align parallel with the vorticity vector. Specifically, when reaction and energy release occur, the outflow for the counterflow configuration will have the greater extensional strain rate in the principal direction aligned with the vorticity. The outflow velocity in the principal direction orthogonal to the vorticity might still experience a higher magnitude than the inflow velocity in the compressive-strain direction due to the density decrease in the flow. However, its increase will be less profound than for the velocity component aligned with the vorticity.

The asymptotic value as $\eta \rightarrow \infty$ in figure 2(d) gives the integrated burning rate in the subgrid volume. It should be realized that the factor $2b$ from (3.24) has not yet been factored into the integral result in the figures; it will reduce the values by an order of magnitude in order to give the average over the volume. It is seen that rotation affects burning and burning rate is negligible in several cases. The decrease in mass flux rate with increasing rotation rate (i.e. vorticity) results in a decreased burning rate through the subgrid volume.

The values of $S_1$ and $S_2$ do have some consequence on the behaviour. In figures 4 and 5, $S_1$ varies between 0.750 and 0.333; $S_2 = 1 -S_1$ and varies accordingly. As $S_1$ decreases and $S_2$ increases, the flame zone moves slightly, the integrated burning rate decreases and the normalized mass flux $f$ through the counterflow decreases. Very interestingly, as $S_1$ decreases, both the $u_{\xi }$ velocity component and $f_1'$ decrease while the $w$ velocity component and $f_2'$ increase. Velocity component $u_{\chi }$ also decreases. As the $S_1$ value moves from 0.500 to 0.333, some reversal of the $u_{\xi }$ velocity occurs in the region of highest temperature and lowest density. In that region, there is inflow (compressive strain) in two directions with outflow (tensile strain) only in the $z$ direction. This implies that, for $S_1 =0.333$, an inflowing particle of material enters at first with decreasing magnitude of $\chi$ and increasing values of $\xi$ and $z$ in Lagrangian time. Then, it changes to decreasing values of both $\chi$ and $\xi$ in that time but remaining with an increasing value for $z$. Note that a case with $K =0.195, \; \omega _{\kappa } = 1.0$ and $S_1 = 0.250$ is not plotted here. It resulted in extinction of the flame.

Figure 4. Diffusion flame with varying strain rate, for $K= 0.195, \omega _{\kappa } = 1.0$. Blue, $S_1 =0.750, S_2 = 0.250$; red, $S_1 =0.500, S_2 = 0.500$; purple, $S_1 =0.333, S_2 = 0.667$. (a) Enthalpy, $h/h_{\infty }$. (b) Fuel mass fraction, $Y_F$. (c) Product of mass ratio and oxygen mass fraction, $\nu Y_O$. (d) Integral of reaction rate, $\int \dot {\omega }_F \,\textrm {d} \eta$. (e) Reaction rate, $\dot {\omega }_F$.

Figure 5. Diffusion flame with varying strain rate, for $K= 0.195, \omega _{\kappa } = 1.0$. Blue, $S_1 =0.750, S_2 = 0.250$; red, $S_1 =0.500, S_2 = 0.500$; purple, $S_1 =0.333, S_2 = 0.667$. (a) Mass flux per area, $f= -\rho u_{\chi }$. (b) Velocity component, $f_1' = u_{\xi }/(S_1 \xi )$. (c) Velocity component, $f_2' = w/(S_2 z)$. (d) Velocity component, $u_{\chi }$.

As seen from figures 3 and 5, because of density gradients, velocity gradients exist in the new rotating coordinate system such that vorticity components result in the $z$ and $\xi$ direction from $f_1''$ and $f_2''$, respectively. However, those velocity components and thereby the associated vorticity components have an antisymmetry. Thus, the induced circulation around the flamelet due to density gradients is zero. Only the circulation due to the imposed vorticity $\omega _{\kappa }$ will exist. These findings are consistent with the results of Sirignano (2019, 2021a).

Clearly, the combination of fluid rotation, variable density and three-dimensional structure has major consequences for flamelet behaviour.

3.7. Uncoupled premixed flame calculations

The premixed flame is treated with consideration of a mixture that is 12 % propane by mass and 43.5 % oxygen by mass (stoichiometric proportion), with the remainder an inert gas. This flame structure requires a much larger Damköhler number to exist when compared with the diffusion flame. It extinguishes under high strain rate (i.e. low residence time in the counterflow). In figures 6 and 7, the rate of strain values are held constant at $S_1 =0.750$ and $S_2 = 0.250$ while $Da$ and/or $\omega _{\kappa }$ are varied. The value of $K$ varies between $15$ and $40$ while $\omega _{\kappa }$ varies between $0$ and $1$. At lower values of the imposed vorticity, the flame extinguishes. Increased vorticity strengthens the centrifugal effect, increasing the efflux of the counterflow in the $z$ direction aligned with the vorticity vector.

Figure 6. Premixed flame: influence of Damköhler number and vorticity, for $S_1 = 0.750 , S_2 = 0.250.$ Blue, $\omega _{\kappa } = 1.0, K= 40$; red, $\omega _{\kappa } = 1.0, K= 30$; dashed blue, $\omega _{\kappa } = 1.0, K= 15$; purple, $\omega _{\kappa } = 0, K=30$. (a) Enthalpy, $h/h_{\infty }$. (b) Fuel mass fraction, $Y_F$. (c) Product of mass ratio and oxygen mass fraction, $\nu Y_O$. (d) Integral of reaction rate, $\int \dot {\omega }_F \,\textrm {d} \eta$. (e) Reaction rate, $\dot {\omega }_F$.

Figure 7. Premixed flame: influence of Damköhler number and vorticity, for $S_1 = 0.750 , S_2 = 0.250.$ Blue, $\omega _{\kappa } = 1.0, K= 40$; red, $\omega _{\kappa } = 1.0, K= 30$; dashed blue, $\omega _{\kappa } = 1.0, K= 15$; purple, $\omega _{\kappa } = 0, K= 30$. (a) Mass flux per area, $f= - \rho u_{\chi }$. (b) Velocity component, $f_1' = u_{\xi }/(S_1 \xi )$. (c) Velocity component, $f_2' = w/(S_2 z)$. (d) Velocity component, $u_{\chi }$.

Figure 6 shows that modest changes for the Damköhler number $Da$ can result in differences in burning rate while figure 7 shows modification in the flame velocity. As $Da$ increases, the flame shifts to a position in the counterflow with higher incoming velocity and mass flux. All three components of the velocity generally increase as $Da$ increases. There is an increase in the chemical reaction rate leading to an increase in the premixed flame speed and a flame position farther from the counterflow interface where the incoming velocity is larger in magnitude. At the lowest value of $Da$ shown in the plots, the flame extinguishes as the reaction rate decreases because, in the spatially varying velocity field, it cannot obtain the needed residence time. The higher velocities result from the combination of an increase in mass burning rate and a decrease in density.

Vorticity produces the centrifugal effect that decreases outflow rate in the $\chi$ direction while increasing the outflow rate in the $z$ direction. It also causes the flame to stabilize at a larger distance from the interface at $\eta =0$. The inflow mass flux, reflected through $f$, decreases as $\omega _{\kappa }$ increases, implying an increase in the residence time.

It is noteworthy that the premixed flames require $Da$ values that are orders of magnitude larger than required for diffusion flames. We can expect that they therefore are less likely to appear in highly strained turbulent flows.

3.8. Uncoupled multi-branched flamelet calculations

In figures 8 and 9, results are shown for a multi-flame configuration with $K= 0.200$ where a fuel-rich mixture with $Y_F = 2/3$ and $Y_O = 1/3$ exists on one side of the counterflow flame and a fuel-lean mixture with $Y_F =1/12$ and $Y_O = 11/12$ exists on the other side. This allows for the possibility of multi-flame structures as found by Sirignano (2021a), in contrast to the FPV approach which disallows multi-branched or premixed flame behaviour by using a mixture fraction $Z$ domain that extends for $0 \leq Z \leq 1$. The figures show that three flame structures appear based on this new theory: a fuel-lean premixed flame on the left, a diffusion flame in the centre and a fuel-rich premixed flame on the right. As noted by Sirignano (2021b), the two premixed flames are not independent structures but rather depend substantially on heat flux from the diffusion flame. The figures also show that, as vorticity and rotation rate increase, the burning rate for each of the three flame branches and mass flux through the counterflow decrease. The location of each of the flames shifts towards the fuel-lean side as rotation rate increases. Also, the fraction of the outflow in the direction aligned with the vorticity increases substantially with increasing rate of rotation.

Figure 8. Multi-branched flame with $K \!=0.200$ and varying vorticity, for $S_1 =0.75 ; S_2 = 0.25$. Blue, $\omega _{\kappa } \!= 0$; red, $\omega _{\kappa } = 1.50$; purple, $\omega _{\kappa } = 1.0$. (a) Enthalpy, $h/h_{\infty }$. (b) Fuel mass fraction, $Y_F$. (c) Product of mass ratio and oxygen mass fraction, $\nu Y_O$. (d) Integral of reaction rate, $\int \dot {\omega }_F \,\textrm {d} \eta$. (e) Reaction rate, $\dot {\omega }_F$.

Figure 9. Multi-branched flame with $K \!=0.200$ and varying vorticity, for $S_1 =0.75 ; S_2 = 0.25$. Blue, $\omega _{\kappa } \!= 0$; red, $\omega _{\kappa } = 1.50$; purple, $\omega _{\kappa } = 1.0$. (a) Mass flux per area, $f= - \rho u_{\chi }$. (b) Velocity component, $f_1' = u_{\xi }/(S_1 \xi )$. (c) Velocity component, $f_2' = w/(S_2z)$. (d) Velocity component, $u_{\chi }$.

For the multi-branched flame as well as the simple diffusion flame, fluid rotation, variable density and three-dimensional structure combine to have major consequences for the behaviour.

For both the multi-branched flame and the single diffusion flame, values of $S_1$ and $S_2$ have a consequence on the behaviour. In figures 10 and 11, $S_1$ varies between 0.750 and 0.333. Also, $S_2 = 1 -S_1$ and varies accordingly. Again, as $S_1$ decreases and $S_2$ increases, the flame zone moves slightly, the integrated burning rate decreases and the normalized mass flux $f$ through the counterflow decreases. As $S_1$ decreases, the $u_{\chi }$ velocity component increases while the $u_{\xi }$ velocity component decreases and, below $S_1 = 0.500$, the $u_{\xi }$ velocity reverses direction in the region of highest temperature and lowest density. So, again in that low-density region, there is inflow (compressive strain) in two directions with outflow (tensile strain) only in the $z$ direction. The results show in figure 12 that for $S_1 =0.333$, an inflowing particle of material enters at first with decreasing magnitude of $\chi$ with Lagrangian time and increasing values of $\xi$ and $z$ in time. However, a reversal for the $\xi$ direction is seen with all the outflow going in the $z$ direction. As noted earlier, the behaviour is consistent with direct numerical simulation findings concerning alignment of vorticity and principal strain axes for reacting flows by Nomura & Elghobashi (1993).

Figure 10. Multi-branched flame with varying strain rate, for $K= 0.200, \omega _{\kappa } = 1.0$. Blue, $S_1 =0.750, S_2 = 0.250$; red, $S_1 =0.500, S_2 = 0.500$; purple, $S_1 =0.333, S_2 = 0.667$. (a) Enthalpy, $h/h_{\infty }$. (b) Fuel mass fraction, $Y_F$. (c) Product of mass ratio and oxygen mass fraction, $\nu Y_O$. (d) Integral of reaction rate, $\int \dot {\omega }_F \,\textrm {d} \eta$. (e) Reaction rate, $\dot {\omega }_F$.

Figure 11. Multi-branched flame with varying strain rate, for $K= 0.200, \omega _{\kappa } = 1.0$. Blue, $S_1 =0.750, S_2 = 0.250$; red, $S_1 =0.500, S_2 = 0.500$; purple $S_1 =0.333, S_2 = 0.667$. (a) Mass flux per area, $f= -\rho u_{\chi }$. (b) Velocity component, $f_1' = u_{\xi }/(S_1 \xi )$. (c) Velocity component, $f_2' = w/(S_2 z)$. (d) Velocity component, $u_{\chi }$.

Figure 12. Multi-branched flamelet with reversed flow, for $K= 0.200, \omega _{\kappa } = 1.0$; $S_1 =0.333, S_2 = 0.667$.

Interestingly, the partially premixed flames in the multi-branched structure are driven by heat flux from the diffusion flame as shown by Sirignano (2021b). Thus, they survive extinction at much lower $Da$ values than independent premixed flames can.

4. Shear flow on the resolved scale: planar jet

The ultimate objective is to connect this new flamelet model with RANS, unsteady RANS and LES computations. Thereby, a limit on dimensionality or time dependence of the large-scale resolved flow would not be sought. The longer-term goal is to couple the subgrid flamelet model to a multi-scale, fully unsteady LES. However, in this first application which is aimed towards gathering understanding of the method and not yet towards solving a new problem, a much simpler resolved scale is sought. So, we seek a resolved scale that is time-averaged with minimal dimensionality. This drives the consideration towards a simple mixing-length examination of two-dimensional mixing and reacting flows. This simple example is only a demonstration of the use of the new flamelet model; in the longer term, many of the simplifying assumptions in this section should be abandoned for the model to achieve its full potential. In the following discussion, the turbulent planar mixing layer is considered.

We desire a solution for the two-dimensional, time-averaged, turbulent shear layer with variable density. The density can vary due to variations in temperature and/or composition. Pressure gradients in the mixing layer are considered negligible and the boundary-layer approximation is employed. Pope (2000) and White (2005) present useful overviews on the topic. Görtler (1942) treated a non-reacting incompressible case, finding a one-dimensional similar solution with $y/x$ as the independent similarity variable. We have a long-standing Illingworth (1949) template with laminar flow to modify the relations for the similar compressible flow. However, the reacting flow does not allow for a similar solution in either a laminar or a turbulent case (except for infinite reaction rate with diffusion flames). Here, we take a pathway using a time-averaged two-dimensional turbulent planar mixing layer. We use a mixing-length concept for the eddy viscosity and diffusivities.

4.1. Reacting shear-layer analysis

We consider the averaged turbulent flow, e.g. steady-state two-dimensional flow. The density is variable. The pressure gradient is zero and the boundary-layer approximation is used. The governing equations for the time-averaged velocity components $u, v$ in $x, y$ space are

(4.1)\begin{equation} \left. \begin{gathered} \frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} = 0,\\ u\frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \frac{1}{\rho}\frac{\partial}{\partial y}\left( \rho \nu_t \frac{\partial u}{\partial y} \right),\\ u\frac{\partial h}{\partial x} + v \frac{\partial h}{\partial y} = \frac{1}{\rho Pr_t}\frac{\partial}{\partial y}\left( \rho \nu_t \frac{\partial h}{\partial y} \right) - \varSigma_{m=1}^N h_{f, m} \dot{\omega}_{m,rs} +\frac{1}{\rho}\varPhi, \\ u\frac{\partial Y_m}{\partial x} + v \frac{\partial Y_m}{\partial y} = \frac{1}{\rho Sc_t}\frac{\partial}{\partial y}\left( \rho \nu_t \frac{\partial Y_m}{\partial y} \right) + \dot{\omega}_{m,rs} ; \quad m=1, 2,\ldots, N .\end{gathered} \right\} \end{equation}

Here, $\nu _t = \mu _t / \rho$ and $\mu _t$ are the kinematic turbulent viscosity and dynamic turbulent viscosity, respectively. It is assumed that the turbulent Prandtl number $Pr_t$ and turbulent Schmidt number $Sc_t$ are uniform through the flow. The above equations are in a non-dimensional form that uses the resolved scales as reference length and velocity. Thus, quantities such as $\dot {\omega }_{m,rs}$ and $\varPhi$ must be properly scaled before calculations (see § 2). Note that we can replace $\nu _t$ in the equations by $\nu +\nu _t$, where $\nu$ is the molecular kinematic viscosity. In most portions of the flow, $\nu _t > > \nu$. For low mean velocities, the dissipation $\varPhi$ can be neglected.

We transform in standard fashion from $x, y$ to $\bar {x} \equiv x, \bar {y}\equiv \int _0^y (\rho /\rho _{\infty })\,{\textrm {d}y}'$. The transverse component of velocity, $v$, is transformed to the variable $w$ to mimic incompressible flow:

(4.2)\begin{equation} \left. \begin{gathered} w \equiv \frac{\rho}{\rho_{\infty}} v +\frac{u}{\rho_{\infty}}\int_{y_{ref}}^y \frac{\partial \rho}{\partial x} {{\rm d}\kern0.05em y}',\\ \frac{\partial u}{\partial \bar{x}} + \frac{\partial w}{\partial \bar{y}} = 0, \\ u\frac{\partial u}{\partial \bar{x}} + w\frac{\partial u}{\partial \bar{y}} = \frac{\partial}{\partial \bar{y}}\left( \frac{ \rho^2}{\rho_{\infty}^2 } \nu_t \frac{\partial u}{\partial \bar{y}} \right), \\ u\frac{\partial h}{\partial \bar{x}} +w \frac{\partial h}{\partial \bar{y}} = \frac{1}{ Pr_t}\frac{\partial}{\partial \bar{y}}\left( \frac{ \rho^2}{\rho_{\infty}^2 } \nu_t \frac{\partial h}{\partial \bar{y}} \right) - \varSigma_{m=1}^N h_{f, m} \dot{\omega}_{m,rs} +\frac{1}{\rho}\varPhi, \\ u\frac{\partial Y_m}{\partial \bar{x}} + w \frac{\partial Y_m}{\partial \bar{y}} = \frac{1}{ Sc_t}\frac{\partial}{\partial \bar{y}}\left( \frac{ \rho^2}{\rho_{\infty}^2 }\nu_t \frac{\partial Y_m}{\partial \bar{y}} \right) + \dot{\omega}_{m,rs}; \quad m=1, 2,\ldots, N. \end{gathered} \right\} \end{equation}

With the knowledge from (2.1a,b), (2.3) and (3.24) about the impact of the Kolmogorov strain rate on the reaction rate and the relation between strain rates at different levels, we may write

(4.3)$$\begin{align} \rho \dot{\omega}_{m,rs} &=\rho_{dimensional} \dot{\omega}_{m,dimensional} = \rho_{\infty, dimensional}S^*\widetilde{\rho \dot{\omega}}_{m} \nonumber\\ &= \rho_{\infty, dimensional} \frac{S^{*3/2}_{rs}\delta}{\nu^{1/2}}\widetilde{\rho \dot{\omega}}_{m}. \end{align}$$

Note that the resolved-scale reaction rate $\omega _{m,rs}$ is dimensional while the subgrid averaged term $\widetilde {\rho \dot {\omega }}_{m}$ is non-dimensional in (4.3).

Here, the layer thickness is described by a constant value of the slope $\textrm {d}\delta / {\textrm {d} x}$ since a linear growth rate is expected. The constant value can be estimated and adjusted based on the results but is expected to be close to the value for the incompressible, non-reacting mixing layer. Experiments have indicated the range of $0.06$ to $0.11$ according to Pope (2000). Here, a constant $\sigma$ is considered as the reciprocal of the slope so that $\delta (x) = x/\sigma.$

(4.4)\begin{equation} \left. \begin{aligned} S^*_{rs} & \equiv \left|\frac{{\rm d}u}{{\rm d} y}\right| = \frac{\rho}{\rho_{\infty}}\left|\frac{{\rm d}u}{{\rm d}\bar{y}}\right| ,\\ \nu_t & = (\delta(x))^2 \left| \frac{{\rm d} u}{{\rm d} y} \right| . \end{aligned} \right\} \end{equation}

Accordingly, the value of $\nu _t$ will vary with both $x$ and $y$.

From (4.3) and (4.4), the dimensional resolved-scale reaction rate $\dot {\omega }_{m, rs}$ can be related to the non-dimensional subgrid-scale reaction rate $\widetilde {\rho \dot {\omega }_m}$. Equation (3.24) indicates that the non-dimensional reaction rate is proportional to an integral with $f$ in the integrand. Thus, the dimensional reaction rate is proportional to the subgrid-scale strain rate $S^*_1 + S^*_2$. That subgrid strain rate is larger than the resolved-scale strain rate by a factor $Re^{1/2}$ as shown by (2.1a,b).

The boundary conditions on each of the second-order partial differential equations are given by prescribing the $u$ component of velocity and the scalar properties in the two free streams. In addition, there are the upstream inflow conditions.

The system of equations given as (4.2) can be made non-dimensional by using $u_{\infty }$ to normalize velocity components, a downstream length $x_0$ to normalize the $\bar {x}$ and $\bar {y}$ coordinates and $\rho _{\infty }$ and $h_{\infty }$ to normalize the corresponding scalar quantities. The turbulent viscosity $\nu _t$ is non-dimensionalized using the kinematic viscosity $\nu$. The perfect-gas assumption and constant specific heat are assumed so the non-dimensional relation $\rho = 1/h$ holds through this shear layer where the pressure is approximately uniform. We also define here a Reynolds number $Re \equiv u_{\infty } \delta _0 / \nu = u_{\infty } x_0 / (\sigma \nu )$. We consider cases where the fractional difference between $\Delta U = u_{\infty } - u_{-\infty }$ and $u_{\infty }$ is not major. The viscous dissipation is considered negligible because the Mach number is low.

(4.5)\begin{equation} \left. \begin{aligned} \frac{\partial u}{\partial \bar{x}} + \frac{\partial w}{\partial \bar{y}} & = 0, \\ u\frac{\partial u}{\partial \bar{x}} + w\frac{\partial u}{\partial \bar{y}} & = \frac{\partial}{\partial \bar{y}}\left( \frac{ 1}{h^2 } \left(\frac{x}{x_0}\right)^2\frac{1}{\sigma^2}\left|\frac{\partial u}{\partial \bar{y}}\right| \frac{\partial u}{\partial \bar{y}} \right),\\ u\frac{\partial h}{\partial \bar{x}} +w \frac{\partial h}{\partial \bar{y}} & = \frac{1}{ Pr_t}\frac{\partial}{\partial \bar{y}}\left( \frac{ 1}{h^2 } \left(\frac{x}{x_0}\right)^2\frac{1}{\sigma^2}\left|\frac{\partial u}{\partial \bar{y}}\right| \frac{\partial h}{\partial \bar{y}} \right) \\ & \quad - \varSigma_{m=1}^N h_{f, m} \frac{Re^{1/2}}{(\sigma h)^{1/2}}\frac{x}{x_0} \left|\frac{\partial u}{\partial \bar{y}}\right|^{3/2}\widetilde{\rho \dot{\omega}}_{m} , \\ u\frac{\partial Y_m}{\partial \bar{x}} + w \frac{\partial Y_m}{\partial \bar{y}} & = \frac{1}{ Sc_t}\frac{\partial}{\partial \bar{y}}\left( \frac{ 1}{h^2 } \left(\frac{x}{x_0}\right)^2\frac{1}{\sigma^2}\left|\frac{\partial u}{\partial \bar{y}}\right| \frac{\partial Y_m}{\partial \bar{y}} \right) \\ & \quad + \frac{Re^{1/2}}{(\sigma h)^{1/2}}\frac{x}{x_0} \left|\frac{\partial u}{\partial \bar{y}}\right|^{3/2} \widetilde{\rho \dot{\omega}}_{m} ; \quad m=1, 2,\ldots, N. \end{aligned} \right\} \end{equation}

Here, in the non-dimensional processing, we have taken

(4.6)\begin{align} (x_0/u_{\infty}) \dot{\omega}_{m,rs} & = (x_0/u_{\infty})(1/\rho) \frac{S^{*3/2}_{rs}\delta}{\nu^{1/2}}\widetilde{\rho \dot{\omega}}_{m} \nonumber\\ &= [ (x_0/u_{\infty})\frac{u_{\infty}^{3/2}}{x_0^{3/2}}\frac{\delta}{\nu^{1/2}}]_{dimensional} (\rho^{3/2}/\rho) \left|\frac{{\rm d}u}{{\rm d}\bar{y}}\right|^{3/2}\widetilde{\rho \dot{\omega}}_{m} \nonumber\\ &= \left(\frac{u_{\infty}x_o}{\nu}\right)^{1/2}\frac{x}{x_0}\frac{1}{\sigma h^{1/2}}\big|\frac{{\rm d}u}{{\rm d}\bar{y}}\big|^{3/2}\widetilde{\rho \dot{\omega}}_{m} \nonumber\\ &= Re^{1/2}\frac{x}{x_0}\frac{1}{(\sigma h)^{1/2}}\left|\frac{{\rm d}u}{{\rm d}\bar{y}}\right|^{3/2}\widetilde{\rho \dot{\omega}}_{m} . \end{align}

The resolved-scale calculation will provide certain information as inputs to the subgrid flamelet computation. Specifically, some information is needed to determine through reasonable scaling principles (i) constraints on the scalar properties (mass fractions of the reactants and the temperature) and (ii) strain rates and vorticity. Then, the subgrid calculation can give the burning rate and associated energy release rate to the resolved scale. Since the subgrid calculations are quasi-steady, they may be performed in advance with the input–output relation organized in the form of a look-up table or a neural network. Then, the table or neural network may be coupled with the resolved-scale computation.

4.2. Shear-layer results and discussion

Computations are performed for a turbulent diffusion flame in the shear layer. One free stream is pure propane gas while the other is pure oxygen. We assume that the flamelet scale also displays a diffusion flame character. (Clearly, a need exists for further direct numerical simulation studies to inform us about the degree of partial mixing in the turbulence cascade.) At the flamelet scale, we take $\int \dot {\omega }_F \,\textrm {d}\eta =0.4$ which is consistent with results shown in figures 2 and 4. Furthermore, considering $2b = O(10)$ in (3.24), we assign $\widetilde {\rho \dot {\omega }}_F = 0.04$. The burning rate is determined by scaling the flamelet-scale burning rate based upon the large-scale strain rate $S^*_{rs}$; thus, the averaged, large-scale burning rate will vary with $x$ and $y$. Implicitly, the dimensional subgrid-scale burning rate is thereby varying with $x$ and $y$. The burning rate is forced to zero value where $Y_F = 0$ and/or $Y_O = 0$. In the regions where both mass fractions have positive values, the burning rate is prescribed by (4.6). Thereby, discontinuities of the resolved-scale reaction rate can result. This can be rectified in the future by relating boundary conditions for the subgrid flamelet calculation to resolved-scale mass fractions.

The eddy diffusivity is estimated based on the assumption that the shear-layer density-weighted width grows linearly as 0.11$x$. Solutions start at $x_0 = 10$ and are marched to $x = 20.$ Hyperbolic tangents are used for the $\bar {y}$ profiles at $x_0$, except for the enthalpy where a Gaussian profile is superimposed to serve as an igniter. Reynolds number $Re =1000$ based on layer width at $x_0$. One free stream, at positive $Y$ values, is five times faster than the other free stream at negative $y$ values.

Figures 13 and 14 give results for a case where the fuel stream is the faster free stream while the oxygen stream is the slower stream. The profiles widen with increasing downstream distance for both velocity and scalar properties. The peak enthalpy and temperature values remain approximately unchanged but the locations of the peak value and of the centre of the reaction zone shift towards the oxygen-stream side. The sharp cut-offs in burning rate at certain transverse positions occur because the rate is non-zero only where both reactants exist on the large scale but that rate does not depend on the precise mass-fraction values on that scale.

Figure 13. Flow-scale diffusion flame: blue, $x = 20$; red, $x = 15$; purple, $x = 12.5$. Fuel in higher-speed stream, oxygen in lower-speed stream. (a) Enthalpy, $h/h_{\infty }$. (b) Fuel mass fraction, $Y_F$. (c) Product of mass ratio and oxygen mass fraction, $\nu Y_O$. (d) Burning rate, $(x_0/u_{\infty })\dot {\omega }_{m,rs}$.

Figure 14. Flow-scale diffusion flame: blue, $x = 20$; red, $x = 15$; purple, $x = 12.5$. Fuel in higher-speed stream, oxygen in lower-speed stream. (a) Velocity, $u/u_{\infty }$. (b) Velocity gradient, $|\textrm {d}u/{\textrm {d} y}|$.

Figures 15 and 16 give results for a case where the oxygen stream is the faster free stream while the fuel stream is the slower stream. Again, the profiles widen with increasing downstream distance for both velocity and scalar properties. The peak enthalpy and temperature values increase slightly as the locations of the peak value and of the centre of the reaction zone shift towards the oxygen-stream side. Note that the burning zone moves towards the oxygen-rich side with an increase of the downstream distance whether that stream is the faster or slower stream. This presumably occurs because nearly four times the mass of oxygen (compared to propane mass) is consumed in the reaction.

Figure 15. Flow-scale diffusion flame: blue, $x = 20$; red, $x = 15$; purple, $x = 12.5$. Oxygen in higher-speed stream, fuel in lower-speed stream. (a) Enthalpy, $h/h_{\infty }$. (b) Fuel mass fraction, $Y_F$. (c) Product of mass ratio and oxygen mass fraction, $\nu Y_O$. (d) Burning rate, $(x_0/u_{\infty })\dot {\omega }_{m,rs}$.

Figure 16. Flow-scale diffusion flame: blue, $x = 20$; red, $x = 15$; purple, $x = 12.5$. Oxygen in higher-speed stream, fuel in lower-speed stream. (a) Velocity, $u/u_{\infty }$. (b) Velocity gradient, $|\textrm {d}u/{\textrm {d} y}|$.

5. Concluding remarks

A new flamelet model is developed for use in subgrid modelling for analysis of turbulent combustion by RANS and LES. This flamelet model presents certain key advances: (i) non-premixed flames, premixed flames or multi-branched flame structures are allowed to appear naturally without prescription; (ii) the impacts of shear strain and vorticity (and associated centrifugal effects) on the flames are determined; (iii) the applied subgrid strain rates and vorticity are directly related to the resolved-scale strain rates and vorticity without the use of a contrived progress variable; (iv) the flamelet model is three-dimensional without need for assuming axisymmetry or planar geometry, allowing the physically correct counterflow under the vorticity constraint; and (v) variable density is addressed in the flamelet model. The results indicate that each of these five features introduces consequential, vital physics that is missed by current two-dimensional, irrotational, constant-density flamelet models that assume a priori a non-premixed or premixed flame structure and make no direct connection to shear strain or vorticity on the larger turbulence scales.

Information from direct numerical simulations concerning the relative alignments of the vorticity vector, scalar gradients and principal strain axes provides a basis for a set of assumptions. The analytical framework allows for multi-step, detailed kinetics although the calculations here are limited to one-step propane–oxygen kinetics. The quasi-steady assumption used in previous flamelet models is maintained here.

A similar solution is found for the flamelet model. Sample computations of the flamelet model without coupling to the resolved flow are presented first to demonstrate the importance of the new features of the model. The rotation due to vorticity creates a centrifugal force that generally decreases the mass flux through the flamelet counterflow. Thereby, an increase in residence time and a decrease in burning rate occur. Rotation can thereby allow flames, which would otherwise extinguish, to survive. Premixed flames can only exist with orders-of-magnitude larger values of $Da$ compared to diffusion flamelets or multi-branched flames. Thereby, we can expect that premixed flames will be less likely to survive extinction in a turbulent situation with high strain rates. The partially premixed flames in a multi-branched flame structure are driven by heat flux from the diffusion flame; thereby, they survive at lower values of $Da$.

Scaling laws relate subgrid strain rates and vorticity to resolved-scale quantities for coupling with LES or Reynolds-averaged flows. Given these relations, the burning rate is determined by the rotational flamelet physics and the energy release rate and rate of change of species mass fraction are given to the resolved scale by the rotational flamelet algorithms. The theory does not introduce any new contrived variable such as a FPV. Connection between the two scales is made using long-established variables.

For this initial study, a highly challenging turbulent flow is deliberately avoided. Rather, a simple turbulent flow is resolved with coupling to the rotational flamelet model in order to explore the interaction across the different scales. A two-dimensional, multicomponent, time-averaged planar shear layer with variable density and energy release uses a mixing-length description for the eddy viscosity and is coupled to the new rotational flamelet model. The eddy diffusivity is proportional to the local magnitude of the velocity gradient and grows with downstream distance. These simplifying assumptions are only used in this first demonstration of the application of the subgrid flamelet model. They should not be generalized to be an inherent component of the approach.

The profiles for velocity and scalar properties are seen to broaden with downstream distance in the shear layer, indicating the growing width of that layer. With increasing downstream distance, the zone where burning occurs moves towards the oxygen-rich side (for the chosen propane–oxygen kinetics) whether or not it is the faster stream. That burning zone becomes more narrow with increasing downstream position.

In the future, multi-step kinetics should be utilized with the rotational flamelet model. The model can be used to produce look-up tables or neural networks that can be employed with LES or RANS calculations.

There are several important issues that should be addressed in future studies to allow better matching between the closure model for the subgrid mixing and combustion with the resolved-scale (or time-averaged) flow: (i) the lags, due to the turbulence cascade, in spatial position and time for the coupling between the large-scale strain and the energy release that it affects; (ii) the degree of partial mixing of reactants during the cascade; (iii) the proper choice of the $b$ parameter to represent the optimal subgrid volume size for averaging; and (iv) relative magnitudes of the vorticity and normal principal strain rates on the subgrid scale. These matters can be examined through carefully designed direct numerical simulations and perhaps clever experiments. Attention is needed to improve our knowledge of the statistical relations between resolved-scale quantities for vorticity, strain rates and scalar gradients and those quantities on the flamelet scale. Of course, flamelets can in principle exist across a range of the smaller scales, not just the smallest scales. Perhaps more attention is needed for the details of the turbulence cascade; for example, strain-rate self-amplification as well as vortex stretching could be relevant (Johnson 2021). Also, better determination is needed of the range of scales where mixing and reaction are prominent.