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 Reference Pope2000), 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 Reference Pope2000; White Reference White2005). 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 (Reference Elghobashi and Launder1983) 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. (Reference Hamins, Thridandam and Seshadri1985), Rajamanickam et al. (Reference Rajamanickam, Coenen, Sanchez and Williams2019), Sirignano (Reference Sirignano2019, Reference Sirignano2021c) and López-Cámara et al. (Reference López-Cámara, Jordà Juanós and Sirignano2020) 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.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 (Reference Drozda, Quinlan and Drummond2020). 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 (Reference Illingworth1949) 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. Reference Linan, Martinez-Ruiz, Vera and Sanchez2017; Weiss et al. Reference Weiss, Vera, Linan, Sanchez and Williams2018) 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 Reference Peters2000; Pierce & Moin Reference Pierce and Moin2004) together with an incompressible-flow assumption which gives an overly simplistic solution for the velocity field. Sirignano (Reference Sirignano2019, Reference Sirignano2021a,Reference Sirignanoc) 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 (Reference Sirignano2019, Reference Sirignano2021a). 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 (Reference Westbrook and Dryer1984) 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 (Reference Nomura and Elghobashi1993) 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 (Reference Sirignano2019, Reference Sirignano2021a).

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 (Reference Sirignano2021a), 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 (Reference Sirignano2021b), 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 (Reference Nomura and Elghobashi1993).

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 (Reference Sirignano2021b). Thus, they survive extinction at much lower $Da$ values than independent premixed flames can.

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 (Reference Pope2000). 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}|$.