Mixing and combustion at low heat release in large eddy simulations of a reacting shear layer

Abstract

In this paper, the mixing and combustion at low-heat release in a turbulent mixing layer are studied numerically using large eddy simulation. The primary aim of this paper is to successfully replicate the flow physics observed in experiments of low-heat release reacting mixing layers, where a duty cycle of hot structures and cool braid regions was observed. The nature of the imposed inflow condition shows a dramatic influence on the mechanisms governing entrainment, and mixing, in the shear layer. An inflow condition perturbed by Gaussian white noise produces a shear layer which entrains fluid through a nibbling mechanism, which has a marching scalar probability density function where the most probable scalar value varies across the layer, and where the mean-temperature rise is substantially over-predicted. A more sophisticated inflow condition produced by a recycling and rescaling method results in a shear layer which entrains fluid through an engulfment mechanism, which has a non-marching scalar probability density function where a preferred scalar concentration is present across the thickness of the layer, and where the mean-temperature rise is predicted to a good degree of accuracy. The latter simulation type replicates all of the flow physics observed in the experiment. Extensive testing of subgrid-scale models, and simple combustion models, shows that the WALE model coupled with the Steady Laminar Flamelet model produces reliable predictions of mixing layer diffusion flames undergoing with fast chemistry.

1 Introduction

The accurate prediction of turbulent combusting flows through numerical simulation is of significant interest for combustion system designers. Combustion of fuels normally used in practical systems involves dozens of species undergoing hundreds of reaction steps; hence, numerical simulation of such systems is an extreme computational challenge. The turbulent nature of the flow in which the combustion occurs adds an extra level of complexity to the system. In order to understand the fundamental physics of mixing in combustion systems, it is common to use both simplified reaction mechanisms, and relatively simple flow geometries, such as mixing layers.

Mixing layers form when two parallel fluid streams of differing velocities (and/or densities) merge downstream of a splitter plate. For non-premixed combustion, the reactants are seeded into opposing freestreams and are brought together within the mixing layer. The combustion occurs once the reactants are mixed at the molecular scale. Two-stream shear layers have been studied extensively for almost 80 years, and it has been shown that the geometrical simplicity of the flow configuration belies its complexity. The flow displays a hypersensitivity to its initial conditions, with the state of the upstream boundary layers influencing the mixing layer for many thousands of boundary layer momentum thicknesses downstream of the origin [8]. Large-scale, spanwise-orientated coherent structures were discovered in the turbulent flow at high Reynolds number [13], and subsequently a statistically stationary streamwise-orientated structure was also observed [5, 6, 30]. The ubiquity of these structures for all initial conditions remains an open research question; the spanwise-orientated structure develops far downstream in mixing layers originating from turbulent upstream conditions [53], and the streamwise structure is absent in mixing layers developing from turbulent initial conditions [4].

For mixing layers originating from laminar upstream conditions, it has been shown that the probability density function (p.d.f.) of high-speed fluid concentration in the mixing layer can have a non-marching distribution, where a preferred concentration occupies the visual thickness of the layer [28, 30, 31, 50]. This p.d.f. shape is attributed to an engulfment entrainment mechanism [17], where pure unmixed fluid is entrained into the layer in the braid regions and mixed to an almost uniform composition within the two-dimensional coherent structures. When the mixing layer originates from tripped, or turbulent upstream conditions, the p.d.f. takes a marching form, where the most probable concentration varies linearly across the thickness of the layer [3, 28, 50]. This p.d.f. shape is attributed to a nibbling entrainment mechanism, where the freestream fluid is rapidly mixed at the outer edges of the mixing layer. The mixing layer entrains fluid asymmetrically, with a bias towards the high-speed stream [30, 31]. Models which describe the growth, entrainment, and mixing in the turbulent mixing layer have been developed, based on the assumption that the large-scale coherent structures are present in the layer [10, 17].

Over the last 35 years, a series of experiments studying the effects of heat release on the mixing layer have been performed. The experiments used dilute concentrations of H$_2$ and F$_2$ carried in nitrogen. The hypergolic reaction that occurs between the reactants permits the study of the mixing and entrainment processes that occur in the mixing layer. Low-speed mixing layers with both low-heat release [45], and high-heat release [23] have been studied. Measurements of the temperature rise within the layer showed the presence of hot regions within the flow that corresponded to the large-scale spanwise-orientated structures, the presence of tongues of cool fluid which penetrated far into the layer, and an almost uniform distribution of the temperature rise within the structures. These spanwise-orientated structures were considered to be of the Brown–Roshko form, being quasi-two dimensional in nature [13]. Small temperature gradients along the streamwise extent of the structures were inferred as evidence for the presence of a stationary streamwise structure in the mixing layer, but this structure was not measured directly. The growth rate of the mixing layer reduced with increasing adiabatic flame temperature rise [23], which was attributed to the effect of heat release on shear stress production. For all values of heat release studied, the mean-temperature rise in the layer did not attain the adiabatic flame temperature. Further experiments showed the effects of Reynolds number on the flow [46] and extended the work into the compressible regime [53].

Numerical simulation methods are now sufficiently mature that they can offer new insights into the physics of turbulent flows, that would otherwise be impossible to obtain experimentally. Early studies into the mixing layer were performed on the temporally evolving flow, as the doubly periodic domain substantially reduced the computational cost. The effect of heat release was shown to reduce the entrainment into the temporal layer, reducing the rate of product formation [41, 42]. Simulating the laboratory-frame mixing layer is computationally expensive; hence, early studies of the spatially developing, reacting mixing layer were confined to two-dimensional boxes [21, 22, 26, 55]. These simulations produced qualitative agreement with experiments, but they fail to capture the streaky streamwise structure known to exist in the flow. Simulations of the fully three-dimensional, isothermal mixing layer, have shown that the imposed inflow condition plays a key role in the development of the layer. For initially laminar mixing layers, it is common to impose an inflow condition based on a mean streamwise velocity profile which is perturbed by Gaussian white noise [1, 7, 39]—these simulations produced mean flow statistics that agree well with experiment, but a spatially stationary streamwise structure is absent [37, 40]. An inflow condition produced by an inflow generation method results in a mixing layer which produced excellent statistical agreement with experiment and captured both the quasi-two-dimensional coherent structures, and the spatially stationary streamwise structure, present in the laboratory mixing layer [40]. A recent study of an exothermically reacting, high-speed mixing layer has shown that good mean-temperature rise predictions can be obtained using LES, and that the state of the upstream flow heavily influenced the p.d.f. of the fluid in the layer [29]. To date, no numerical study has successfully replicated the flow structure features reported in the Mungal and Dimotakis experiments [45].

In this research, we simulate the low-speed, high Reynolds number mixing layer flame experiments of Mungal and Dimotakis [45]. The research aims to accurately capture the salient flow physics of the experiment, and to assess the effect of inflow conditions on the predicted flow. As the experiments originated from laminar upstream boundary layers, two commonly used methods to model this type of inflow are employed; the first is a mean streamwise velocity profile onto which Gaussian white noise is superimposed, and the second is an inflow condition obtained from an inflow generator. The inflow conditions obtained from both methods are statistically identical, in terms of the mean and fluctuating velocity profiles. It is expected that the combusting mixing layer which originates from laminar upstream conditions will contain spanwise-orientated turbulent structures and an organised streamwise vortex structure.

This paper is organised as follows. Details of the reference experiment are provided in Sect. 2. Numerical methods utilised in the research are provided in Sect. 3. The set-up of the simulations is described in Sect. 4. Main simulation results are given in Sect. 5. The effect of subgrid-scale modelling is described in Sect. 6, and the influence of the combustion model on the computed flow statistics is discussed in Sect. 7. Conclusions are drawn in Sect. 8.

2 Reference experiment

The reference experiments were performed in a blowdown rig specifically designed to measure two-dimensional, exothermically reacting mixing layers. The test section measured 0.76 $\times {}$ 0.15 $\times {}$ 0.2 (m) in the streamwise (x), vertical (y), and spanwise (z) directions, respectively. The lower guidewall of the test section could be angled in order to minimize the streamwise pressure gradient, such that the low-speed stream speed varied by less than 2%. Experiments were performed with the high-speed velocity $U_1 =$ 22 ms$^{-1}$, and the low-speed velocity $U_2 =$ 8.8 ms$^{-1}$, yielding a velocity ratio parameter value of $R = (U_1 - U_2) / (U_1 + U_2) =$ 0.428. Both streams contained laminar boundary layers at the trailing edge of the splitter plate, with the high-speed stream boundary layer momentum thickness estimated at $\theta _1 =$ 1.63 $\times {}$ 10$^{-4}$ m using Thwaites’s method. Temperature and velocity measurements were recorded at a streamwise distance of $x = 0.457$ m downstream of the splitter plate trailing edge. This equates to a non-dimensional distance of $x / \theta _1 \approx 2800$, which is far into the fully developed region of the self-similar mixing layer [8]. At the measurement station, the local Reynolds number of the mixing layer based on the velocity difference across the layer and its visual thickness was 6.5$\times {}$10$^{4}$.

The heat release in the experiments was provided by the irreversible reaction between concentrations of hydrogen and fluorine, carried in inert nitrogen gas.

(1)

(2)

where k is given in $\mathrm {cm}^3 \mathrm {mol}^{-1}$. The first reaction given in Eq. 1 is an order of magnitude faster than the second one, given in Eq. 2.

(3)

Additional reactions involving scavenging of atomic fluorine, and chain-termination reactions involving third-body collisions, complete the set of reactions [45]. It has been shown that the reactions listed in Eqs. 1–3 dominate the behaviour of the system for the adiabatic flame temperature rise considered here [20]. In this paper, the chemical system is therefore modelled on these three reactions. The Damköhler number is defined here as the ratio of the mixing time to the chemical time, $\text {Da} = \delta _1/\tau \Delta U$, where $\delta _1$ is the mixing layer visual thickness at the measurement station, $\tau$ is the chemical time, and $\Delta U$ is the velocity difference across the layer. For the current conditions, $\text {Da} \approx 10$, indicating that the chemistry is therefore sufficiently fast that the reactions are mixing-limited [44].

The equivalence ratio was defined the ratio of low-speed reactant concentration, $\xi _2$, to the high-speed molar concentration, $\xi _1$, such that

$$\begin{aligned} \phi = \frac{\xi _2 / \xi _1}{(\xi _2 / \xi _1)_s} \ \ = \frac{\xi _2}{\xi _1}, \end{aligned}$$

(4)

as the molar stoichiometric ratio for the H$_2$/F$_2$ reaction is unity. By seeding 1% F$_2$ in the low-speed stream, the concentration of hydrogen seeded in the high-speed stream was varied to obtain equivalence ratios of $\phi = 1$,$\frac{1}{2}$,$\frac{1}{4}$,$\frac{1}{8}$. Flipping the streams in which the reactants were seeded produced equivalence ratios of $\phi =$ 1, 2, 4, 8. In the stream carrying the hydrogen, a 0.03% concentration of nitric oxide was seeded to ensure rapid ignition of the mixture. Adiabatic flame temperature rises in the range of 93K $\le$ $T_{\text {flm}}$ $\le$ 165K were achieved. The maximum Richardson number of the experiments was $\sim 10^{-4}$, some two orders of magnitude lower than the value required for the effects of buoyancy to become important [32].

Temperature data were recorded by eight cold wires placed across the thickness of the layer. The mean-temperature profiles obtained from the experiments were repeatable to within 2%. The largest source of error in the measurement arose from an uncertainty of the absolute reactant concentration in the freestream of 3–5%, resulting in an uncertainty of the absolute value of heat release of the same magnitude. Velocity measurements were obtained using a rake of fifteen Pitot tubes placed across the flow, and no quantification of the measurement error was quoted [45].

3 Numerical methods

The LES code is based on the low-Mach number approximation of the spatially filtered governing equations. The filtered low-Mach number equations of continuity, momentum, and species transport are written as

$$\begin{aligned}&\displaystyle \frac{\partial \bar{\rho }}{\partial t} + \frac{\partial {\bar{\rho } {\tilde{u}}_{i}}}{\partial {x_{i}}} = 0, \end{aligned}$$

(5)

$$\begin{aligned}&\displaystyle \frac{\partial {\bar{\rho } {\tilde{u}}_{i}}}{\partial {t}} + \frac{\partial \bar{\rho }{\tilde{u}}_i{\tilde{u}}_j}{\partial x_j} = - \frac{\partial {\bar{p}}}{\partial {x_{i}}} + \frac{\partial {\tau _{ij}}}{\partial {x_{j}}} - \frac{\partial {T_{ij}}}{\partial {x_{j}}}, \end{aligned}$$

(6)

$$\begin{aligned}&\displaystyle \frac{ \partial \bar{\rho }{\tilde{Y}}_\alpha }{\partial t} + \frac{\partial \bar{\rho } {\tilde{u}}_i {\tilde{Y}}_\alpha }{\partial x_i} = \frac{ \partial {\tilde{J}}_i^\alpha }{\partial x_i} - \frac{\partial M_i^\alpha }{\partial x_i} + {\tilde{\omega }}_\alpha , \end{aligned}$$

(7)

where a quantity $\overline{\phi }$ denotes a spatially filtered variable, ${\tilde{\phi }}$ denotes a Favre-averaged variable. The velocity field is $u_i$, with $\bar{p}$, ${\tilde{Y}}_\alpha$, and ${\tilde{\omega }}_\alpha$ representing pressure, species mass fraction, and reaction rate terms for the species $\alpha$ ($\alpha = 1, 2, \ldots , n$), respectively. Under the assumption of a Newtonian fluid obeying Fick’s Law of diffusion, the filtered species mass flux is given by ${\tilde{J}}_i = \mu / Sc (\partial \tilde{Y_\alpha }/\partial x_i)$, where $\mu$ is the dynamic viscosity, and Sc is the Schmidt number. The filtering operation introduces extra terms into the governing equations which must be modelled to close the system—the subgrid stress and mass flux tensors are given by $T_{ij} = \overline{\rho }( \widetilde{u_i u_j} - \tilde{u_i} \tilde{u_j})$ and $M_{i} = \overline{\rho }( \widetilde{u_i \xi } - \tilde{u_i} {\tilde{\xi }})$.

A subgrid-scale model must be used to close the momentum equation. Applying the Boussinesq turbulent viscosity assumption leads to

$$\begin{aligned} T_{ij} = - 2 \mu _{sgs} {\tilde{S}}_{ij} + \frac{1}{3} \delta _{ij}T_{kk}. \end{aligned}$$

(8)

Two distinct models are employed in this study to calculate the subgrid viscosity, $\mu _{sgs}$. The first is the WALE model [47]—in this model, the subgrid viscosity is computed from

$$\begin{aligned} \mu _{sgs} = \bar{\rho }(C_w{\Delta })^2 \frac{(S^d_{ij}S^d_{ij})^{3/2}}{(\overline{S}_{ij}\overline{S}_{ij})^{5/2} + (S^d_{ij}S^d_{ij})^{5/4}}, \end{aligned}$$

(9)

where $S^d_{ij} = \frac{1}{2}(\overline{g}^2_{ij} + \overline{g}^2_{ji} - \frac{1}{3} \delta _{ij}\overline{g}^2_{kk})$, $\overline{g}^2_{ij} = \overline{g}_{ik}\overline{g}_{kj}$, and $\overline{g}_{ij} = \partial \overline{u}_i /\partial x_j$. This model is attractive for shear layer simulations, as the model predicts zero eddy viscosity in regions of pure shear. The second model is the standard Smagorinsky model [54], given by

$$\begin{aligned} \mu _{sgs} = \bar{\rho }(C_s \Delta )^2 |{\tilde{S}}_{ij}|, \end{aligned}$$

(10)

where $\Delta$ is the filter width, and the magnitude of the strain-rate tensor is given by $|{\tilde{S}}_{ij}| = \sqrt{ 2 {\tilde{S}}_{ij} {\tilde{S}}_{ij}}$ with ${\tilde{S}}_{ij} = \frac{1}{2} (\partial {\tilde{u}}_i/\partial x_j + \partial {\tilde{u}}_j/\partial x_i)$. The limitations of this model are well-known, as the model coefficient $C_s$ is not universal, and the model predicts finite subgrid-scale viscosity in regions of laminar flow.

The species mass fraction equation is closed with the commonly used gradient-diffusion model. A constant turbulent Schmidt number is assumed in this model, such that

$$\begin{aligned} M_i^\alpha = -\frac{\mu _{sgs}}{Sc_t}\Big (\frac{\partial {\tilde{Y}}_\alpha }{\partial x_i}\Big ), \end{aligned}$$

(11)

and in this study the subgrid-scale Schmidt number is set to $Sc_t = 0.5$.

The chemical reactions of Eqs. 1–3 are sufficiently fast that they may be assumed to be significantly faster than the mixing rate of the reactants. Essentially this means that the combustion between the species occurs once the fluids mix. This permits the use of a mixture fraction approach [51] where the local mass fractions of the chemical species, the temperature, and the density, are uniquely related to a single scalar function called the mixture fraction, Z. The Favre-filtered mixture fraction transport equation is written as

$$\begin{aligned} \frac{ \partial \bar{\rho }{\tilde{Z}}}{\partial t} + \frac{\partial \bar{\rho } {\tilde{u}}_i {\tilde{Z}}}{\partial x_i} = \frac{ \partial {\tilde{J}}_{Z,i}}{\partial x_i} - \frac{\partial M_{Z,i}}{\partial x_i}, \end{aligned}$$

(12)

with ${\tilde{J}}_{Z,i}$ and $M_{Z,i}$ denoting the filtered, and subgrid-scale, mixture fraction mass fluxes, respectively. The subgrid term is closed in a similar manner to Eq. 11. In this research, there is one ‘fuel’ stream, and one ‘oxidiser’ stream; hence, the filtered mixture fraction can be expressed as

$$\begin{aligned} {\tilde{Z}} = \frac{{\tilde{Y}}_\alpha - {\tilde{Y}}_{\alpha ,O}}{{\tilde{Y}}_{\alpha ,F} - {\tilde{Y}}_{\alpha ,O}}, \end{aligned}$$

(13)

where ${\tilde{Y}}_{\alpha ,F}$ and ${\tilde{Y}}_{\alpha ,O}$ are the filtered species mass fractions of the fuel and oxidiser, respectively. For the hydrogen-fluorine considered in this research, the species which constitutes the fuel is not obvious [45]; hence, the stream carrying the hydrogen is assigned a mixture fraction value of unity, and the stream carrying the fluorine is assigned a mixture fraction value of zero.

The combustion is modelled using the well-known flamelet modelling approach of Cook et al. [15]. The chemical reactions of Eqs. 1–3 are computed a priori, and stored in a run-time look-up table. The species mass fractions, the fluid density, and the fluid temperature are functions of the mixture fraction, ${\tilde{Z}}$, the mixture fraction variance, ${\tilde{Z}}_v^2$, defined as

$$\begin{aligned} {\tilde{Z}}_v^2 = C_v \Delta ^2 | \nabla {\tilde{Z}}|^2, \end{aligned}$$

(14)

where $C_v$ = 0.1, and the scalar dissipation rate, $\chi$, defined as

$$\begin{aligned} \chi = C_\chi \frac{(\mu _{sgs} + \mu )}{\overline{\rho } \text {Sc}_t} |\nabla {\tilde{Z}}|^2, \end{aligned}$$

(15)

where $C_\chi =$ 2. The equation for $Z_v$ can be obtained from local equilibrium assumptions [9], and scalar dissipation rate can be interpreted as the inverse of a characteristic diffusion time [49].

Two variants of this model are utilised here; the first assumes a state of chemical equilibrium whereby the effects of flame stretching are ignored; therefore, only a single flamelet library is required for this model at $\chi =$ 0. This model is called the Chemical Equilibrium Model (ECM). The second model accounts for the scalar dissipation caused by flame stretching in the turbulent flow; hence, multiple flamelets are tabulated in this case. This model is called the Steady Laminar Flamelet Model (SLFM). Owing to the assumptions on which these models are based, it is known that the amount of chemical product formed may be over-predicted when compared to experiment [15]. The look-up tables are generated using the ANSYS CFD software suite.

The governing equations of motion are solved using a finite volume method on a staggered mesh. Second-order accurate central differencing schemes are used to evaluate terms in the momentum equation, and a third-order accurate Total Variation Diminishing (TVD) scheme is used for the convective and diffusive terms in the mixture fraction transport equation. A TVD scheme is used for the mixture fraction transport in order to minimise out of bounds errors in the solution of the equation [16]. The outflow condition applied in the simulations is a convective condition similar to that used in previous uniform density studies of the mixing layer flow [38]. The second-order accurate Adams–Bashforth method is used to temporally advance the governing equations. The pressure field is computed using a multi-grid method. The temperature, density, and species mass fractions are obtained from the look-up tables through bilinear interpolation from the ECM model, and trilinear interpolation for the SLFM model.

Two inflow condition generation methods are used in this study to assess the impact of the inflow conditions on the simulated flame. The first inflow condition is obtained by imposing a mean streamwise velocity profile at the inlet plane, which is perturbed at each time step by Gaussian white noise. This type of inflow condition is commonly used for the simulation of mixing layers which originate from initially laminar conditions [1, 7, 39]. The second inflow condition is produced by recycling and rescaling procedure of Xiao et al. [59]. This procedure is similar to the recycling method of Lund et al. [33], except that (a) virtual domains are used in which the inflow condition is generated, and (b) the flow within the entire virtual domain is rescaled, at intervals, to a target set of statistics. The method has been used to simulate mixing layers [37, 40, 59], liquid jet breakup [60], and accelerating nozzle flows [58], with excellent results obtained in all cases. A short summary of the procedure is included here. The inflow generation method requires small domains in which the inflow condition is generated for each stream. These are placed upstream of the main mixing layer domain and have a mesh resolution that matches the initial region of the mixing layer domain described in Sect. 4. The virtual domains extend 14$\delta _{0.99}$ upstream of the trailing edge, where $\delta _{0.99}$ is the high-speed side boundary-layer thickness. The flow at a sampling plane positioned $\delta _{0.99}$ upstream of the trailing edge is recycled onto the virtual domain inlet plane at each time step, and at a given interval the entire virtual domain is rescaled to a set of target flow statistics. The flow field that passes through the downstream end of the virtual domain is then fed into the main mixing layer simulation. Considering the streamwise velocity component only, the mean velocity in the virtual domain, $\bar{U}^{(n+1)}$, is computed using a weighted average such that

$$\begin{aligned} \bar{U}^{(n+1)}(y) = \frac{k\Delta t}{T_s} \langle U(x,y,z,t)\rangle _{x,z} + \Bigg (1 - \frac{k\Delta t}{T_s} \Bigg ) \bar{U}^{n}(y), \end{aligned}$$

(16)

where $\Delta t$ is the simulation time step, k is the integer number of time steps between rescaling operations, $T_s$ is a characteristic time scale in the flow, U(x, y, z, t) is the instantaneous flow field in the virtual domain, and $\langle \rangle _{x,z}$ denotes averaging in both the streamwise and spanwise directions. Here, $k=10$, and $T_s=10\delta _{0.99}/U$, where U is the freestream velocity of the flow. The root mean squared (r.m.s.) velocity fluctuation is similarly computed through

$$\begin{aligned} u'^{(n+1)}(y) = \sqrt{\frac{k\Delta t}{T_s}\langle [U(x,y,z,t) - \bar{U}^{(n+1)}(y)]^2\rangle _{x,z} + \Bigg (1 - \frac{k\Delta t}{T_s} \Bigg ) [u'^n(y)]^2}. \end{aligned}$$

(17)

At intervals of k time steps, the entire flow field in the virtual domain is rescaled to a target set of statistics through

$$\begin{aligned} {\tilde{u}}^\mathrm{new} (x,y,z,t) = \frac{u'_\mathrm{target}(y)}{u'^{(n+1)}(y)} [U(x,y,z,t) - \bar{U}^{(n+1)}(y)] + \bar{U}_\mathrm{target}(y), \end{aligned}$$

(18)

where ${\tilde{u}}^\mathrm{new}$ is the new velocity component, $\bar{U}_\mathrm{target}$ is the target mean velocity, and $u'_\mathrm{target}$ is the target rms fluctuation. The same operations are performed on the other two velocity components. Validation of this method for initially laminar mixing layers can be found elsewhere [35].

4 Simulation setup

The mixing layer is simulated in a computational domain that extends 0.8 $\times {}$ 0.61 $\times {}$ 0.12 (m) in the streamwise (x), vertical (y), and spanwise (z) directions, respectively. The domain is discretised into 896 $\times {}$ 256 $\times {}$ 256 cells. The grid is refined in the region of the splitter plate, such that the minimum grid spacings are $\Delta x_\mathrm{min} = 2\times {} 10^{-4}$m , $\Delta y_\mathrm{min} =4\times {} 10^{-5}$m in the streamwise and vertical directions, respectively. Geometric expansions are used to stretch the grids in both of these directions, in order to reduce the cell count in regions of low flow variability. The spanwise grid spacing is constant, with a value of $\Delta z = 4.7\times {}10^{-3}$ m. Grids of comparable resolution have been shown to produce reliable mixing layer simulation results [39, 40], and the spanwise domain extent is sufficient to avoid artificial confinement of the mixing layer [34, 36]. The grid is well-refined, when compared to other studies of the high Reynolds number mixing layer [29]. Where the recycling and rescaling method generates the inflow conditions, extra computational domains are required to produce the time-dependent inflow data, as described in Sect. 3. These virtual domains have a mesh resolution that matches the initial region of the mixing layer domain. The computational overhead of the virtual domains adds an extra 30% to the computational cost of the RRM-type simulations.

The upper and lower guidewalls of the computational domain are fixed horizontal. The vertical extent of the computational domain is substantially larger than the test section of the experiment (0.61 m in the simulation compared to 0.15 m in the experiment), in order to minimise the streamwise pressure gradient that arises owing to the growth of the mixing layer. In the simulations, the low-speed stream velocity changes by less than 1% along the streamwise extent of the domain, a value which is well within the experimental estimate of 2% [45]. A standard convective outflow condition is imposed at the outflow boundary, which has been shown to be passive [56]. The inflow boundary conditions are obtained from the methods described in Section 3, with the mean streamwise velocity profile, $\overline{U}$, and root mean squared fluctuation profiles, $u'_i$, imposed in both streams shown in Fig. 1. The spanwise boundaries are periodic.

Fig. 1

Inflow conditions for the simulations

Table 1 Chemical compositions, and adiabatic flame temperature rises, of the mixing layer simulations

Equivalence ratios in the range $\frac{1}{8} \le \phi \le 8$ are simulated, matching those reported experimentally. The runs are outlined in Table 1. In all simulations, the stream carrying hydrogen is seeded with a 0.03% concentration of nitric oxide. The time step for the simulations is $\Delta t = 6 \times {} 10^{-7}$s. The simulations are run for two convective flow through times in order to produce a statistically stationary flow-field. Statistical sampling and flow visualisation images are then collected over a period of twenty convective flow through times. A convective flow through time is defined as the time required for a fluid parcel to traverse the streamwise extent of the computational domain at the convection velocity, $U_c = 0.5(U_1 + U_2)$.

Simulations with the Gaussian white noise perturbations are hereinafter referred to as ‘WN’, and simulations with the recycling and rescaling method inflow are hereinafter referred to as ‘RRM’.

5 Results

Simulations are performed for conditions listed in Table 1 using both the WN-type, and the RRM-type, inflow conditions. Statistical information is recorded at $x = 0.457$ m, and at the mid-span of the computational domain, unless otherwise stated. The effects of subgrid-scale modelling, and combustion modelling, on the simulation results are discussed in Sects. 6 and 7, respectively.

5.1 Flow structure

Density gradient magnitude contour maps of the flow field are shown in Fig. 2 for the $\phi = 1$ case. The instantaneous visualisations are spanwise-averaged, for consistency with experimental visualisations [46]. Both visualizations show large-scale, spanwise-orientated structures present in the simulated mixing layer. The structures in the WN-type inflow simulation (Fig. 2a) are quite irregular, with thick interconnecting braid regions separating the structures. The structures present in the RRM-type calculation (Fig. 2b) have a well-defined central core, with thin interconnecting braid regions. There are deep incursions of freestream fluid into the mixing layer, suggesting that the intermittency at the edges of the mixing layer is lower in the RRM-type inflow case, when compared to the WN-type simulation. Qualitatively at least the structures present in the RRM-type inflow simulation bear remarkable similarity to those observed in a low heat release experimental visualisation [24].

Fig. 2

Spanwise-averaged density gradient magnitude flow visualisations of the $\phi = 1$ simulations

Fig. 3

Spanwise-averaged temperature fields of the $\phi = 1$ simulations, captured at the same time instant as the images in Fig. 2

The corresponding spanwise-averaged temperature rise field, in both $\phi =$ 1 simulations is shown in Fig. 3. The temperature rise, T, is measured as the temperature above ambient, such that the ambient temperature in the images is equivalent to $T/T_\mathrm{flm} = 0$. In the WN-type inflow simulation, the regions of elevated temperature in the cores of the structures are visible, and it is noticeable that the braid regions are also quite hot with respect to the ambient. In the RRM-type inflow simulation, the regions of elevated temperature in the structure cores are very well defined, and the thin interconnecting braid regions between the structures are quite cool. As noted in the schlieren image of Fig. 2b, there are incursions of cold freestream fluid deep into the layer from both freestreams. The images in Figs. 2 and 3 suggest that the pattern of mixing, and therefore reaction, is different in both simulations.

Fig. 4

Sections of time traces of the normalised temperature rise in the mixing layer, at $\phi =$ 1. Measurements recorded at $x = 0.457$ m. The samples shown here corresponds to the passage of six large-scale structures through the measurement plane

More quantitative information on the flow structure can be gleaned from time traces of the temperature field. At the $x =0.457$ m measurement station, the streamwise ($y-z$) cross plane is sampled at a frequency of 3.33 kHz. A short section of the time trace taken at mid-span ($z =$ 0.06 m) is shown in Fig. 4 for both simulation types, with $\phi =$ 1. The time axis in the images is reversed, so that the structures have the same orientation as those present in the flow visualisation images of Fig. 3. The traces show the passage of six large-scale spanwise structures through the measurement plane. In the WN-type simulation, the time trace shows an irregular structure in the mixing layer. There are no obvious well-defined cores in the trace, and there is no uniformity in the temperature rise within each structure. In addition, cold freestream fluid does not penetrate far into the mixing layer. In the RRM-type simulation, a fundamentally different flow structure is apparent. Three key features are evident in the trace shown in Fig. 4b; (1) There are clearly defined hot regions within the flow, which correspond to the structure cores; (2) The temperature rise within the structures is almost uniform across the width of the layer; (3) There are tongues of cold freestream fluid that penetrate deeply into the mixing layer. These features were also prominent in the experimental time trace data at $\phi =$ 1 [45], and they are consistent with other observations of the laboratory mixing layer [18, 30, 57].

Fig. 5

Perspective views of the large-scale structure passing through the $x = 0.457$ m measurement station, in the $\phi = 1$ simulations

The geometry of the large-scale structures is quantified using pseudo-three-dimensional perspective views of the flow passing through the measuring station. Iso-surfaces of the mixture fraction value ${\tilde{Z}} =$ 0.01 define the upper edge of the mixing layer, and the iso-surface is coloured by the normalised streamwise velocity ${\tilde{u}} / U_1$. Local streamwise velocities with a value higher than the freestream value have lighter colours, and darker colours denote a local streamwise velocity lower than the freestream value. An overshoot in the local streamwise velocity indicates the presence of a coherent motion of spanwise vorticity (i.e. a large scale turbulent vortex structure) passing through the measurement station [11, 12]. The temperature contour map within the structure is shown to aid interpretation of the position of the structures.

Short sections of the perspective views at $x = 0.457$ m from the $\phi =$ 1 calculations are shown in Fig. 5. For the WN-type inflow simulation (Fig. 5a), the overshoots in the streamwise velocity have no clear spanwise coherence, which implies that the spanwise vorticity is irregularly distributed within the structures. The overshoots in the local streamwise velocity are also present in the interconnecting braid regions between structures. This allows spanwise vorticity to be exchanged between neighbouring structures as the structures convect downstream, facilitating the continuous linear structure growth observed in simulations that originate from WN-type inflow conditions [39, 40]. There is no obvious secondary streamwise structure that convects along with the primary structure, and it is clear that the turbulent vortex structures present in the WN-type simulations have an inherently three-dimensional geometry. The perspective view of Fig. 5a is similar to that obtained from a non-reacting mixing layer simulation originating from WN-type inflow conditions [40].

The perspective view from the RRM-type simulation is shown in Fig. 5b. In this simulation, the local streamwise velocity overshoots run parallel to the span, demonstrating that the turbulent vortex structures are spanwise coherent. There are clear undulations across the span of the mixture fraction iso-surface, which indicate the presence of an organised streamwise structure in the flow. These undulations can easily be traced from one structure to the next, providing evidence that the streamwise structure is, in a mean sense, spatially stationary. These features have been observed in experiments of the aqueous mixing layer [6, 27], and in other simulations of spatially developing mixing layers originating from a RRM-type inflow condition [40]. The structure present in the RRM-type simulations is of the Brown–Roshko form, where large quasi-two-dimensional turbulent vortex structures have a secondary streamwise structure superimposed upon them [14].

Fig. 6

Time traces from the WN-type inflow simulations, at the extremes of equivalence ratio. Traces recorded at mid span of the $x =0.457$ m measurement station

Fig. 7

Time traces from the RRM-type inflow simulations, at the extremes of equivalence ratio. Traces recorded at mid span of the $x =0.457$ m measurement station

Time traces recorded in the WN-type inflow simulations at $\phi = \frac{1}{8}$, 8 are shown in Fig. 6. The effect of the change in reactant concentrations is to shift the peak temperature rise in the layer towards the stream carrying the lean reactant. Both traces in Fig. 6 show that the mixing layer has an irregular structure, with elevated temperatures present in both the structure cores, and the interconnecting braids. There is a well-defined gradient in the temperature rise from the lean reactant side, to the rich reactant side. In both simulations, the temperature towards the stream carrying the rich reactant is quite cool with respect to the adiabatic flame temperature. The traces recorded in the RRM-type inflow simulations for the same equivalence ratios are shown in Fig. 7. The peak temperature rise is also shifted towards the lean reactant stream in these simulations. The three clear features observed at $\phi =$1 (Fig. 4b) are also present in these traces; well-defined structure at elevated temperature, a near-uniform temperature rise within the structures, and penetration of cold fluid deep into the layer in the braids. The so-called ramps in the temperature along the streamwise extent of each structure are visible in the traces shown in Fig. 7; for $\phi = \frac{1}{8}$, the large-scale spanwise structures are rich in high-speed reactant and the structure trailing edge (on the left of each structure in Fig. 7a) is hotter than the leading edge (on the right of each structure in the image). These ‘ramps’ in the streamwise direction of the time trace were also observed by Mungal and Dimotakis [45]. The ramps have the opposite orientation for the $\phi =$ 8 case shown in Fig. 7b—the structures are rich in low-speed reactant, and the leading edge of the structure is hotter than its trailing edge. This behaviour has been associated with the presence of an organised streamwise vortex structure in the mixing layer [45].

The above results show significant differences in the distribution of the instantaneous temperature field within the mixing layer, as a function of the imposed inflow condition.

Fig. 8

Vorticity thickness profiles obtained from the simulations

5.2 Mean flow statistics

The effects of the inherent differences in the large-scale structure on the computed flow statistics are now analysed.

The vorticity thickness of the mixing layer is obtained from

$$\begin{aligned} \delta _\omega = \frac{U_1 - U_2}{\partial \overline{U}/\partial y |_\mathrm{max}}, \end{aligned}$$

(19)

and the profiles from all simulations are shown in Fig. 8. From the WN-type inflow calculations, presented in Fig. 8a, b, an increasing heat release results in a decrease in the vorticity thickness of the flow. This trend holds regardless of the stream in which the lean reactant is placed. In the RRM-type inflow condition simulations, the vorticity thickness evolution of the flow is only very weakly a function of the temperature rise. The trends observed in the RRM-type simulations agree favourably with the experiments, where no particular change in vorticity thickness for flame temperature rises in the range 93K $< T_{\text {flm}}<$ 165 K was observed [23].

Profiles of the mean streamwise velocity in the $\phi =$ 1, and $\phi = \frac{1}{8}$, are shown in Fig. 9a. The vertical co-ordinate is normalised by the distance of the virtual origin, $y / (x - x_0)$. Generally good agreement between the predicted mean streamwise velocity profile, and the experimental data, is observed. A slight thinning of the layer with increasing heat release is observed in the RRM-type simulations, with this trend being more pronounced in the WN-type inflow cases.

The Reynolds stress tensor of the reacting flow is given by $R_{ij} = \overline{\rho u'_i u'_j} / (\rho _1 U_0^2$), where $\rho _1$is the density of the high-speed stream, and $U_0$ is the velocity difference across the layer. Figure 9b shows the predictions of the primary Reynolds normal stress along with comparable experimental [4], and numerical simulation [48, 52], data. The vertical co-ordinate here is normalised by the vorticity thickness to allow comparison with the data from the literature. The simulation predictions of are within the range of those reported elsewhere, and it is noticeable that the magnitude of the stress decreases with increasing heat release in the WN-type simulations. This effect is less marked in the RRM-type calculations. Similar behaviour is observed for the primary Reynolds shear stress, as shown in Fig. 9c.

Profiles of the predicted mean-temperature rise, $\overline{T}$, normalised by the adiabatic flame temperature rise, $T_{\text {flm}}$, are shown in Fig. 10a for $\phi =$ 1. The mean-temperature profiles are corrected to the position of the dividing streamline, $y^*$, for an equivalence ratio of unity, obtained from

$$\begin{aligned} \int _{y^*}^{\infty } \frac{\overline{\rho U}}{\rho _1 U_1} \Bigg ( \frac{\overline{U}}{U_1} -1 \Bigg ) {\text {d}}y = \int _{-\infty }^{y^*} \frac{\overline{\rho U}}{\rho _1 U_1} \Bigg ( \frac{{U_2}}{U_1} - \frac{\overline{U}}{U_1} \Bigg ) {\text {d}}y, \end{aligned}$$

(20)

such that the normalised vertical co-ordinate is $\eta = (y - y^*) / (x - x_0)$.

In both simulations, the adiabatic flame temperature is not reached, in the mean, at any point across the layer. In the WN-type inflow simulation, the maximum temperature rise, $\overline{T}_{\text {max}} = 0.722 T_{\text {flm}}$, over-predicting the experimental value obtained from a curve-fit of the measured data by 10%. The RRM-type simulation produces excellent maximum temperature rise data, when compared with the experimental data. It is worth noting, however, that no experimental point resides where the maximum temperature rise occurs for this equivalence ratio.

Figures 10b-c show that mean-temperature rise profile, along with the maximum and minimum temperatures recorded across the layer at any time during the simulation, for the WN-type and RRM-type simulations, respectively. The maximum temperature data in the WN-type simulation show that the layer can be quite hot across its thickness, and also that the minimum temperature within the central part of the layer is significantly above ambient. The relatively high minimum temperature across the layer is caused by the elevated temperature in the braid regions, as seen in Fig. 4a. The prolonged temperature rise in the core of the layer leads to the over-prediction of the mean-temperature rise. In the RRM-type simulation, the maximum temperature approaches the adiabatic flame temperature rise across the entire thickness of the layer, and the minimum temperature is barely above ambient. The duty cycle of alternating periods of hot fluid within the quasi-two-dimensional structures, and of cold fluid in the interconnecting braid regions, observed in Fig. 4b, results in a mean-temperature rise profile which agrees extremely well with the experiment. The duty cycle described above was a key feature of the flame reported in the experiment [45], and its presence in the RRM-type simulation indicates that the flow physics is correctly captured in this case. The 10% change in the predicted temperature rise between the two simulation types is therefore linked to the pattern of mixing within the structures present in each simulation—in the WN-type simulation, mixed fluid appears to be more prevalent within the layer, which results in a higher mean-temperature rise. This is explored further in Sect. 5.3.

Fig. 9

Flow statistics obtained from simulations at equivalence ratios $\phi = 1, \frac{1}{8}$. Comparable experimental data obtained from archival literature

Fig. 10

Temperature rise data obtained from the $\phi =$ 1 simulations

The mean-temperature rise profiles from the $\phi =$ 8 simulations are shown in Fig. 11a. The WN-type inflow simulation over-predicts the maximum temperature rise in the layer by 22%, and the thickness of the layer is substantially smaller than the experimental observation. The RRM-type simulation predicts a temperature distribution that agrees well with the experiment, and the maximum temperature rise is over-predicted by 9%. The minimum recorded temperature in the WN-type simulation again shows the presence of elevated temperature within the central part of the layer throughout the simulation. The shift of the maximum temperature towards the lean reactant side, and the drop in the maximum temperature towards the low-speed side, are evident in the data. For the RRM-type simulation, the maximum temperature is near-uniform across the layer, and the minimum temperature barely exceeds ambient. These trends are repeated in the $\phi = \frac{1}{8}$ simulations shown in Fig. 12, with the bias shifted towards the opposing stream.

The mean-temperature rise profiles of the remaining simulations are shown in Fig. 13. It is noticeable that all WN-type simulations markedly over-predict the maximum temperature rise, and that the layer is substantially thinner than the experimental flow. The RRM-type simulations generally produce predictions which are in excellent agreement with experiment, and the over-prediction of the maximum temperature rise increases with increasing equivalence ratio.

Table 2 outlines the properties of the mixing layer at $x =0.457$ m. The area under the mean-temperature rise profile is computed from $\frac{1}{T_\mathrm{flm}}\int \overline{T} d\eta$ and gives a measure of the amount of product formed in the layer. The visual thickness, $\delta _1$, is obtained from the vertical points where the mean temperature rise is 1% of the maximum. The RRM-type simulations produce values for both the area, and the normalised visual thickness, that are in good agreement with the reference data [45]. The WN-type simulations predict an area that is much lower than the experimental data, and the visual thickness of the simulated layer is substantially under-predicted for all equivalence ratios.

The product thickness, $\delta _{P1}$, can be obtained from the mean-temperature rise profile, such that

$$\begin{aligned} \delta _{P1} = \int _{-\infty }^{\infty } \frac{C_p \overline{T}(y)}{\xi _{1} \Delta Q} {\text {d}}y, \end{aligned}$$

(21)

where $C_p$ is the molar heat capacity of the carrier gas, and $\Delta Q$ is the heat release per mole of reactant. The product thickness, normalised by the visual thickness, is shown in Fig. 14 for all simulations. When compared to the experimental data, the WN-type simulations over-predict the normalised product thickness for all values of the equivalence ratio. The reasons for this discrepancy are discussed in Sect. 5.3. The RRM-type simulations show good agreement with experiment for all values of the equivalence ratio, except for $\phi =$ 4, 8. The deviation for these equivalence ratios is caused by the over-prediction of the maximum temperature rise in the layer, and hence an over-estimation of the product formed.

Fig. 11

Temperature rise data obtained from the $\phi =$ 8 simulations

Fig. 12

Temperature rise data obtained from the $\phi = \frac{1}{8}$ simulations

Fig. 13

Mean-temperature rise profiles from the $\phi = \frac{1}{4}$,$\frac{1}{2}$, 2, 4 simulations

5.3 Mixing and entrainment

Time samples of the flow passing through the $x =0.457$ m measurement station are used to construct p.d.f.’s of the high-speed fluid concentration, $\xi$. The measured high-speed concentrations are binned into intervals of $\epsilon =$ 0.025, such that concentration in the range 1-$\epsilon \le \xi \le$ 1 is marked as pure high-speed fluid, and concentration in the range 0 $\le \xi \le \epsilon$ is denoted as pure low-speed fluid.

The p.d.f.’s recorded in the $\phi =$ 1 simulations are shown in Fig. 15a, b. The WN-type inflow simulation shows a peak at $\xi =$ 0 in the low-speed stream, a peak at $\xi =$ 1 in the high-speed stream, and the most probable value of the high-speed concentration varies in a marching fashion across the vertical extent of the mixing layer. The so-called marching p.d.f.’s have been observed in experiments where the upstream flow was highly disturbed [50], and in other numerical simulations of spatially developing mixing layers with inflow conditions perturbed by Gaussian white noise [1]. In the RRM-type inflow simulation, peaks are observed in the freestreams, and a peak is also observed within the mixing layer. The most probable concentration, $\xi = 0.593$, is invariant across the thickness of the mixing layer, yielding what is known as a ‘non-marching’ p.d.f. distribution. A non-marching p.d.f. was inferred from the temperature data in the reference experiment [45], and in several other mixing layer experiments which originated from laminar upstream conditions [28, 30, 31, 50].

Table 2 Effects of subgrid-scale modelling on the simulated mixing layer properties

The method of Koochesfahani and Dimotakis [31] is used to compute mixed fluid statistics from the probability density functions. The probability of finding mixed fluid at any concentration is defined as

$$\begin{aligned} P_m(y) = \int _\epsilon ^{1-\epsilon } P(\xi , y) {\text {d}}\xi . \end{aligned}$$

(22)

The probability of finding pure unmixed fluid from the high-speed stream, $P_1$, and the probability of finding pure unmixed fluids from the low-speed stream, $P_2$ are defined as

$$\begin{aligned} P_1(y) = \int _{1-\epsilon }^1 P(\xi , y) {\text {d}}\xi , \ \ \ P_2(y) = \int _0^{1-\epsilon } P(\xi , y) {\text {d}}\xi . \end{aligned}$$

(23)

The profiles of $P_m$, $P_1$, and $P_2$ obtained from the p.d.f.’s of Fig. 15a, b are shown in Fig. 15c. In these plots, the vertical co-ordinate is normalised by the mixing layer visual thickness, $\delta _1$, for consistency with previous experiments [31, 50]. It is significant that there is a substantially greater probability of finding mixed fluid within the layer in the WN-type inflow simulation, when compared to its RRM-type inflow counterpart. In the WN-type inflow case, $P_m$ attains a value of unity across the inner 40% of the mixing layer, meaning that the pure freestream fluid does not penetrate the core of the mixing layer at any time. The probability of finding pure unmixed fluid within the layer rapidly drops from unity in the freestream, to zero within the layer. This rapid mixing of pure freestream fluid is evident in the time trace images of Fig. 4a, and leads to the nonzero minimum temperature rise within the layer observed in Fig. 10b.

In the RRM-type inflow simulation, a value of $P_m =$ 1 is not attained anywhere in the layer (the maximum value attained is ($P_m =$ 0.992), signifying that pure freestream fluid can penetrate far into the mixing layer. The penetration of pure freestream fluid occurs in the braid regions, where tongues of unmixed fluid enter the mixing layer—features than can be observed in Fig. 4b, and which result in a near-zero minimum temperature rise within the layer as shown in Fig. 10c. The $P_1$ and $P_2$ curves show that the probability of finding freestream fluid in the layer is much higher than in the counterpart WN-type calculation. This higher probability is caused by the duty cycle of alternating periods of hot fluid (structure cores), and cool fluid (braid regions), owing to the dynamics of the coherent structures present in this simulation.

The average mixed-fluid concentration within the layer can be computed from

$$\begin{aligned} \xi _m(y) = \frac{\int _\epsilon ^{1-\epsilon } \xi P(\xi ,y) {\text {d}}\xi }{\int _\epsilon ^{1-\epsilon } P(\xi ,y) {\text {d}}\xi } \end{aligned}$$

(24)

and is shown in Fig. 15d for the $\phi =$ 1 cases, along with the mean concentration profile, $\overline{\xi }$. The ‘S-shaped’ curve present in the RRM-type simulation is caused by the presence of pure unmixed fluid within the layer, and the shape of this curve is in good agreement with previous experiments [30, 45, 50]. The $\xi _m$ curve of the WN-type simulation more closely follows the mean concentration, as the pure freestream fluid is more rapidly mixed at the outer edges of the layer in this case.

Fig. 14

Dependence of the product thickness on the equivalence ratio

Fig. 15

Mixed fluid statistical information for simulations with $\phi =$ 1

Fig. 16

Mixed fluid statistical information for simulations with $\phi = \frac{1}{8}$

Mixed fluid statistics extracted from the $\phi = \frac{1}{8}$ are shown in Fig. 16. The trends noted above for $\phi =$1 are present in these data; the p.d.f. of the WN-type simulation is marching, whilst the p.d.f. of the RRM-type simulation is non-marching with a preferred concentration at $\xi =$ 0.583. There is a higher probability of finding mixed fluid in the WN-type simulation, as shown in the $P_m$ and $\xi _m$ curves. For all equivalence ratios studied here, the p.d.f. obtained from the WN-type simulations is marching in nature. Similarly a non-marching p.d.f. is obtained from the RRM-type simulations at all equivalence ratios.

The r.m.s. concentration fluctuation is given by

$$\begin{aligned} \xi '(y) = \sqrt{\int _0^1 (\xi - \overline{\xi })^2 P(\xi ,y){\text {d}}\xi }, \end{aligned}$$

(25)

and is shown in Fig. 17 for the $\phi =$ 1, $\frac{1}{8}$ cases. The unmixedness parameter is also shown on these plots, defined as $\xi '^2 / \overline{\xi }(1-\overline{\xi })$ [19]. An unmixedness value of zero signifies perfect mixing, and values approaching unity signify no mixing, but stirring of the fluid. For both equivalence ratios shown in Fig. 17, the RRM-type inflow simulations show a bimodal r.m.s. concentration profile, with peaks towards each freestream. The WN-type inflow simulations have a much flatter r.m.s. concentration profile across the vertical direction. The unmixedness profiles for $\phi =$ 1 show that the layer in the WN-type inflow simulations is well-mixed, as the value of this parameter is quite low, and fairly uniform across the vertical extent of the layer. In the RRM-type inflow simulations, high values of the unmixedness parameter are reported towards the outer edges of the mixing layer at $\phi =$ 1, owing to the large tongues of cool fluid that penetrate deep into the layer in the braid regions. At $\phi = \frac{1}{8}$, the unmixedness of the RRM-type simulation is lower than at $\phi =$ 1 , although peaks are still observable towards the other edges of the layer.

The entrainment ratio, $E_v$ is defined as the relative proportion of fluid that is entrained into the layer from each freestream. The entrainment ratio can be computed directly from the p.d.f. [31]. The total mixed fluid concentration for the layer is given by

$$\begin{aligned} {\xi }_M = \frac{\int _{-\infty }^{+\infty } P_m(y) {\xi }_m(y){\text {d}}y}{ \int _{-\infty }^{+\infty } P_m(y) {\text {d}}y}. \end{aligned}$$

(26)

The entrainment ratio, $E_v$, is then linked to the total mixed concentration of the layer by $\xi _M = E_v / (1 + E_v)$. Values of $E_v > 1$ mean that more fluid is entrained from the high-speed stream—a feature that is present in mixing layer experiments [30, 31, 43].

For the p.d.f.’s shown in Fig. 15a, b, the computed entrainment ratios at $x =0.457$ m are $E_v =$ 1.56 for the WN-type inflow simulation, and $E_v =$ 1.44 for the RRM-type inflow simulation. At $\phi =$8, the entrainment ratios are $E_v =$ 1.36, 1.28 for the WN-type, and RRM-type inflow simulations, respectively. The computed values are within the range of values reported experimentally for this velocity ratio parameter [23, 30, 43]. The entrainment model proposed by Dimotakis predicts a value of $E_v =$ 1.29, for $R = 0.428$ [17].

Given that the entrainment ratio between the simulations types is similar for a given value of the equivalence ratio, the wide discrepancies in the statistical data relating to the mixing statistics in the predicted shear layers warrants further investigation.

Key to the discrepancies in the mixing statistics is the evolution of the turbulent coherent structures. In the WN-type inflow simulations, the large-scale turbulent vortex structures have an inherently three-dimensional geometry. The irregular distribution of spanwise vorticity within the mixing layer, inferred from perspective view images, permits the exchange of spanwise vorticity between neighbouring structures. This exchange process facilitates the continuous and linear growth of the large-scale structures [40]. The irregular structure geometry permits spanwise vorticity to exist in the thick braid regions between structures. The presence of spanwise vorticity in the braids leads to the rapid mixing of the freestream fluid entrained into the layer. This is clearly shown in the mixed fluid probability plots of Figs. 15c and 16c, where no pure freestream fluid is found within the central part of the layer. The rapid mixing leads to low values of the unmixedness parameter within the layer (Fig. 17), the over-prediction of the normalised product thickness, and the over-prediction of the maximum mean-temperature rise in the layer. This type of entrainment process is termed ‘nibbling’, as pure freestream fluid is drawn into the layer and rapidly mixed. A schematic of this entrainment process is shown in Fig. 18a. The nibbling entrainment mechanism has been reported in mixing layer simulations originating from WN-type inflow conditions [2], but the computed concentration p.d.f.’s contradict experimental observations of mixing layers originating from laminar upstream conditions [28, 30, 50].

Fig. 17

Root mean squared concentration fluctuation, and unmixedness parameter, profiles for selected equivalence ratios

Fig. 18

Schematic of the entrainment mechanisms present in the simulations

The results presented above show that the pattern of mixing and entrainment in the RRM-type simulations is markedly different to that which occurs in the WN-type cases. The quasi-two-dimensional structures present in the RRM-type simulations have been shown to grow with a square-root-of-time dependency, and pairing-type interactions must take place between them in order to produce the overall mean self-similar growth of the layer [40]. The concentration of spanwise vorticity within the structures, and the relative absence of spanwise vorticity in the braid regions (inferred from the perspective view of Fig. 5b), means that fresh freestream fluid is brought into the layer in the raid regions (as shown schematically in Fig. 18b), before it is mixed to a roughly uniform composition within the structures. This mechanism of entrainment and mixing in the shear layer is termed ‘engulfment’, or ‘gulping’ [17] and gives rise to the non-marching p.d.f. observed in the RRM-type simulations. This p.d.f. shape has been reported in experiments which originate from ‘clean’ laminar upstream conditions [28, 30, 43, 45, 50].

Given that both simulation types originate from upstream conditions with identical mean and r.m.s. flow statistics, the RRM-type calculations correctly capture the flow physics present in the real flow, because the structures present in these simulations have the correct internal geometry [37, 40]. In order to obtain accurate predictions of the mixing and heat release in turbulent mixing layers, it is clear that the inflow conditions for the simulation must adequately mimic the upstream flow conditions present in the experiment.

6 Influence of subgrid-scale model

Simulations are performed at $\phi = 1$ to assess the influence of the subgrid-scale model on the simulated flow. Six simulations are performed for each inflow condition type; three with the Smagorinsky model at a coefficient value of $C_s =$ 0.1, 0.15, 0.18, and a further three with the WALE model at a coefficient value of $C_w =$ 0.3, 0.48, 0.56. All of these simulations use the Equilibrium Chemistry Model for the chemical modelling. The mean flow statistics from these simulations are shown in Fig. 19, along with the experimental data. Regardless of the inflow condition used, the mean streamwise velocity profiles agree extremely well with the experimental data. The choice of subgrid-scale model has little influence on the velocity field. The r.m.s. streamwise velocity fluctuation profiles, shown in Fig. 19c, d, also show that the choice of subgrid scale has only a weak influence on the turbulence statistics. The mean-temperature rise profiles at $x = 0.457$ m are shown in Fig. 19e, f. The profiles are normalised by the adiabatic temperature rise. For the WN-type inflow simulations, the general shape of the mean-temperature profiles is in reasonable agreement with experiment. The profiles are thinner than the reference data, and the maximum temperature rise is over-predicted. These trends are present for all subgrid-scale model values tested. In the RRM-type inflow simulations, the predictions are substantially improved; the thickness of the layer is closer to the experimental data, and the maximum temperature rise is in close agreement with the reference data.

Fig. 19

Mean flow statistics obtained from the subgrid-scale model testing, at $\phi =$ 1. $T_\mathrm{flm}$ is the adiabatic flame temperature. Experimental data comes from Mungal and Dimotakis [45]

Table 3 outlines the flow properties computed from the mean-temperature rise profiles. In the WN-type inflow simulations, the maximum mean-temperature rise is over-predicted by up to 12% when compared to the reference data, regardless of the choice of subgrid scale model. The width of the layer is under-predicted by up to 23% when compared to the experiment.

For the RRM-type inflow simulations, the maximum mean-temperature rise is predicted to within 3–8% for the ECM model. The prediction of the width of the layer also falls within 8% of the experimental data. The WN-type simulations all show a smaller value for the area when compared to their RRM-type counterparts, owing to the smaller width of the mixing layer in the WN-type cases.

Table 3 Effects of subgrid-scale modelling on the simulated mixing layer properties

Fig. 20

Mean flow statistics obtained from combustion model validation simulations at $\phi =$ 1

The results presented here show that the grid is sufficiently well resolved that the choice of subgrid-scale model has little influence on the flow predictions, for a given inflow condition type. On poorly resolved grids, however, the WALE model is recommended for use in the simulation of initially laminar mixing layers as its ability to predict zero subgrid viscosity in regions of laminar flow is advantageous over the Smagorinsky model [25, 39].

7 Influence of chemistry model

The results from the SLFM simulations are now compared to those obtained using the ECM. For both inflow condition types, the WALE model is used, with the model coefficient set to $C_w =$ 0.56. Fig. 20 shows the mean flow statistics obtained the ECM and SLFM runs. The mean streamwise velocity profiles are once again in excellent agreement with the experimental data. The turbulence statistics for the SLFM cases do not vary substantially from their ECM counterparts in Fig. 19c, d and are not shown here. The mean-temperature rise profiles show that the maximum mean-temperature rise in the SLFM cases is reduced when compared to the ECM cases, and that they are in much better agreement with the experimental data. The improvement in the SLFM case predictions is caused by the fact that the influence of the scalar dissipation rate on the temperature field is accounted for in the SLFM cases. The mean-temperature profiles in the SLFM cases show that the nature of the imposed inflow condition continues to have an influence on the predicted temperature field. The maximum mean-temperature rise in the WN-type simulation is 11% higher than that obtained in the RRM-type simulation. The profile obtained in the RRM-type simulation profile is in excellent agreement with the experimental data. Table 3 shows that $\overline{T}_\mathrm{max}/T_\mathrm{flm}$ is within 1.5% of the experimental data, and that the width of the layer is within 2.5% of the reference value.

On the well-refined grid used in this study, the improvement of the SLFM model over the ECM model in terms of predicting the mean-temperature rise, is marginal. On coarse grids, however, we expect the SLFM choice to be the superior model for the reaction mechanism used here as the effect of the strain rate on the flame is incorporated into the model.

8 Conclusions

In this paper, we have studied the low speed, high Reynolds number, initially laminar mixing layer undergoing exothermic reaction with low heat release. The low heat release permits the study of mixing within the layer, without the fluid dynamics being affected by substantial changes in density. The high Damköhler reactions are modelled using tabulated chemistry, based on a mixture-fraction approach. The effect of inflow conditions on the simulated flow is studied through using two common methods for prescribing laminar upstream conditions—a mean profile perturbed by Gaussian white noise, and an inflow condition generated using a recycling and rescaling method.

The choice of inflow condition has a significant effect on the mixing and heat release in the simulated layer. Simple modelling of the inflow condition using white noise perturbations results in a mixing layer that substantially over-predicts the mean-temperature rise in the layer, and substantially under-predicts the mixing layer thickness. In addition, the marching nature of the high-speed concentration probability density function is entirely at odds with experimental observation of mixing in the high Reynolds number mixing layer originating from clean laminar flow conditions. These problems are shown to reside in the fact that the internal geometry of the turbulent vortex structures in such simulations promotes a nibbling entrainment mechanism, which rapidly mixes the fluid drawn into the layer, resulting in the over-prediction of product formation. Such simulations, therefore, do not correctly capture the flow physics present in the laboratory flow.

When a more sophisticated inflow modelling technique is employed, the structures present in the flow are of a quasi-two-dimensional form, on which a secondary streamwise structure rides passively. These structures produce a duty cycle of alternating regions of hot and cold fluid in the layer, which correspond to the structure core, and braid regions, respectively. The dynamics of the large-scale structures is such that the engulfment mechanism draws pure freestream fluid deep into the layer in the braid regions, which is then mixed to a roughly uniform composition in the structure core. Simulations utilising the inflow generation technique successfully replicate all of the flow features observed in the reference experiment and capture the mean-temperature rise, and produce formation, in the mixing layer to a good degree of accuracy.

The research in this paper shows that accurate modelling of inflow conditions is essential to obtain reliable information on the mixing of reactants in shear flow, particularly for flows originating from laminar upstream conditions. Further work will study the mixing layer at high-heat release.