Abstract

A lifted turbulent jet flame issuing into a vitiated coflow is investigated using the conditional moment closure. The conditional velocity (CV) and the conditional scalar dissipation rate (CSDR) submodels are chosen such that they are fully consistent with the moments of the presumed probability density function (PDF). The CV is modelled using the PDF-gradient diffusion model. Two CSDR submodels based on the double integration of the homogeneous and inhomogeneous mixture fraction PDF transport equations are implemented. The effect of CSDR modelling is investigated over a range of coflow temperatures ( ) and the stabilisation mechanism is determined from the analysis of the transport budgets and the history of radical build-up ahead of the stabilisation height. For all , the balance between chemistry, axial convection, and micromixing, and the absence of axial diffusion upstream of the stabilisation height indicate that the flame is stabilized by autoignition. This conclusion is confirmed from the rapid build-up of ahead of , , and . The inhomogeneous CSDR modelling yields higher dissipation levels at the most reactive mixture fraction, which results in longer ignition delays and larger liftoff heights. The effect of the spurious sources arising from homogeneous modelling is found to be small but nonnegligible, mostly notably within the flame zone.

1. Introduction

The investigation of stabilisation mechanisms in lifted flames is an active area in combustion research. A substantial number of experimental and numerical studies have been dedicated to the understanding and analysis of the process by which lifted flames are stabilised. Several theories have been proposed in the literature. Vanquickenborne and van Tiggelen [1] suggest the premixed flame propagation stabilisation theory. Their experimental findings indicate that the fuel and the oxidiser are premixed ahead of the base of lifted diffusion flames and that stabilisation takes place at stoichiometric locations where the local mean velocity and the turbulent burning velocity of a stoichiometric premixed flame are equal. Peters and Williams [2] suggest that stabilisation is achieved by the local quenching of diffusion flamelets due to excessive straining. Cabra et al. [3] conduct experimental measurements for a lifted turbulent jet flame issuing into a coaxial vitiated coflow consisting of the hot combustion products of a lean premixed /air flame. The coflow conditions allow for the possibility of stabilisation by means of autoignition. It is postulated that autoignition takes place as inert fuel parcels are convected downstream and become well mixed with the hot coflow. However, premixed flame propagation remains another possible stabilisation mechanism.

The Cabra flame is attractive for the validation of combustion models due to its simple configuration and well-defined boundary conditions. It is extensively investigated using PDF methods [36]. Masri et al. [4] use the composition PDF approach coupled with the turbulence model and conclude that flame is largely controlled by chemical kinetics. They identify key reactions and show that the predicted liftoff height and the composition at the base of the flame are highly sensitive to the rate parameters of these reactions. Gordon et al. [6] use the same computational methodology in order to determine the underlying stabilisation mechanism. They develop two numerical indicators for the distinction between stabilisation via autoignition and premixed flame propagation. These indicators are based on the analysis of the transport budgets of convection, diffusion, and reaction and the history of radical build-up ahead of the reaction zone. A convection-reaction balance without significant contribution from axial diffusion, along with an early and rapid build-up of ahead of the remaining radicals, indicate that stabilisation occurs via autoignition, whereas an axial diffusion-convection balance with negligible reaction accompanied by a simultaneous build-up of all radicals points to stabilisation by means of premixed flame propagation. With the aid of these indicators, Gordon et al. [6] conclude that autoignition is the controlling mechanism for a wide range of coflow temperatures. Cao et al. [5] employ the joint velocity-turbulence frequency-composition PDF method and perform a detailed parametric study using several mixing models and chemical kinetic mechanisms, over a wide range of jet and coflow conditions. They show that although the modelling of mixing affects the prediction of liftoff height, the flame is mainly controlled by chemical kinetics. They also demonstrate that the liftoff height increases with increasing jet and coflow velocities and decreases with increasing jet and coflow temperatures. In all of the PDF studies described above, autoignition was found to be the controlling stabilisation mechanism. In the context of CMC, Stankovic and Merci [7] and Navarro-Martinez and Kronenburg [8] investigate this flame using Large Eddy Simulation (LES) with a first-order closure for the conditional chemical source. The main difference between the two LES-CMC studies lies in the treatment of the CSDR. The first employs the Amplitude Mapping Closure (AMC) [9], which produces symmetric profiles in mixture fraction space. The second employs a different conditioning technique described in [10], which allows for skewed profiles with the possibility of local minima and maxima in mixture fraction space. In both sets of calculations the coflow temperature ranges between 1020 and 1080 K. The first authors conclude that stabilisation is controlled by autoignition for all coflow temperature. The second authors draw the same conclusions for high coflow temperatures; however, they report a transition to flame stabilization by premixed flame propagation as the coflow temperature is reduced. The different conclusions are most likely attributed to the treatment of the CSDR and, to a lesser extent, to the employed numerical methods, discretization schemes, and grid resolutions in physical and mixture fraction spaces. Patwardhan et al. [11] also use the first-order CMC to study this flame. They employ the standard model for the computation of the turbulence and mixing fields. To a certain degree, their conclusions in regard to the controlling stabilisation mechanism are in line with those of Navarro-Martinez and Kronenburg [8]. Premixed flame propagations are found to be the controlling stabilisation mechanism of the lifted flame for low and intermediate coflow temperatures (1025 and 1045 K), whereas the flame is stabilised by autoignition for high coflow temperatures (1080 K). To close the CMC equations, they employ the -distribution for the PDF, the linear velocity model [12] for the CV, and the AMC for the CSDR. The assumption of homogeneous Gaussian turbulence is inherent to the adopted CV and CSDR closures. The linear velocity model is exact only if the velocity and the mixture fraction are jointly Gaussian [13, 14], which is in general not the case for inhomogeneous flows. Though being consistent with the first moment of the PDF, this CV model is inconsistent with the second moment [15]. As for the AMC, this closure is a particular solution of the homogeneous PDF transport equation, wherein the employed mapping function is Gaussian. Patwardhan et al. [11] use a -distribution instead. Further, the fact that the AMC does not account for inhomogeneous effects leads to an inconsistent CMC implementation. When such closure is employed, the integration of the PDF-weighted CMC equations over the mixture fractions space yields the unconditional equations in addition to extra spurious source terms [14, 16]. This term may lead to inaccurate results and misleading conclusions, depending on the magnitudes of these sources. As such, the CV and CSDR closures employed in [11] do not provide a fully consistent CMC implementation.

In this work, the lifted hydrogen jet flame of Cabra et al. [3] is revisited and investigated using the two-dimensional first-order CMC and a modified version of the model. The conditional submodels required for the closure of the CMC equations are selected such that they are fully consistent with the moments of the presumed PDF. Additionally, a transport equation is solved for the Favre-averaged scalar dissipation rate rather employing the traditional “equal time scales” algebraic approach. The simulations are set up such that the settings of Patwardhan et al. [11] are reproduced as closely as possible. This study is aimed at determining (1) the effect of micromixing on liftoff via the modelling of the CSDR, (2) the response of the flame to changes in the coflow temperature, (3) the controlling stabilisation mechanism, (4) the influence of the spurious sources, and (5) the impact of chemical kinetics.

2. Experimental Configuration

The flame under investigation is the lifted turbulent jet flame of Cabra et al. [3]. A schematic of the experimental setup is shown in Figure 1 and the flow conditions are summarized in Table 1. The burner consists of a central jet surrounded with a coaxial flow composed of the hot combustion products of a lean premixed /Air flame stabilized on a perforated disk. The perforated disk is surrounded by an exit collar for the purpose of delaying the entrainment of ambient air into the coflow. The nozzle exit is placed above the surface of the disk in order to allow the fuel to exit into the coflow with a uniform composition. The experimental criterion for the determination of the liftoff height ( ) is taken as the first location where the mass fraction of reaches 600 ppm. The measured height is equal to 10 nozzle diameters.

3. Mathematical Model

3.1. Conditional Moment Closure

The modelling of the turbulence-chemistry interactions is carried out using the first-order CMC [14]. In CMC, reactive scalars such as the species mass fractions and the temperature are conditionally averaged with respect to the mixture fraction, , and their conditional transport equations are solved. The conditional average of the mass fraction of a species is defined as where is the mass fraction of the species, denotes the conditional average of the quantity to the left of the vertical bar subject to the quantity to its right, and is a sample variable of , with . Following the decomposition approach [14], is written as the sum of its conditional average and a fluctuation : such that . Substitution of (2) into the transport equation of species , followed by the conditional averaging of the resulting expression, leads to [14]: where is the scalar dissipation rate. The species and mixture fraction diffusivities in (3)–(6) are assumed to be equal to the molecular diffusivity, , and the Lewis number is set to unity. In (3), is the CV, is the CSDR, and is the conditional chemical source. The expressions of (see (4)) and (see (5)) are unclosed. They are modelled using the primary closure hypothesis. Klimenko and Bilger [14] show that, given finite Schmidt numbers, all the terms of scale as the inverse of the Reynolds number (Re). Therefore, this term is neglected in high Re applications. By neglecting the conditional density and diffusivity fluctuations, they further show that and reduces to where is the Favre PDF of mixture fraction and is the conditional turbulent flux. Using these approximations, the conditional species transport equation takes the form: The left-hand side term of (8) is the rate of change of the conditional species mass fraction. The first term on the right-hand side (r.h.s.) represents transport by means of convection, the second corresponds to transport by conditional turbulent fluxes (diffusion in physical space), the third accounts for micromixing (diffusion in mixture fraction space), and the fourth is the chemical source. The conditional temperature equation is derived in a similar fashion and is given by    where is the chemical source term, is the total number of species in the mixture, and is the conditional radiative source term. The unsteady pressure term in (9) is neglected due to the fact that the investigated flame is open to the atmosphere, and hence temporal pressure changes are expected very small.

The unconditional (Favre) averages of the reactive scalars are calculated by integrating their conditional values weighted by over the mixture fraction space:

3.2. Conditional Submodels

The quantities , , , , , , and in (8) and (9) are unclosed and require further modelling. The submodels employed in this study are discussed in this section.

3.2.1. Probability Density Function

The PDF is presumed using the -distribution. The Favre -PDF is given by where the parameters and are related to the mean and variance of the mixture fraction by and with , and is the beta function.

3.2.2. Conditional Turbulent Fluxes

The conditional turbulent fluxes are modelled using the gradient diffusion assumption: where is the turbulent diffusivity. The constant is equal to 0.09 and the turbulent Schmidt number is set to 0.7 following Jones and Whitelaw [17]. Richardson and Mastorakos [18] suggest adding a correction to (12) in order to account for counter-gradient transport effects. When applied to a lifted flame, they conclude that the inclusion of counter-gradient effects leads to a slight increase in liftoff height and to a decrease in flame thickness. This extension is not employed in the current study. However, it is worth investigating in future studies.

3.2.3. Conditional Velocity

The PDF-gradient model proposed by Pope [19] is employed to model the CV fluctuations. The resulting expression for the CV takes the form: As tends to zero, diverges to , depending on the sing of the gradient of . This behavior is of minor importance, because low-probability events have a negligible effect on mixing [20]. One important feature of the PDF-gradient model is its consistency with both the first and second moments of the mixture fraction [15, 20]. This does not hold for the commonly used linear model [12], which is only consistent with the first moment [15].

3.2.4. Conditional Scalar Dissipation Rate

Girimaji’s Model. Girimaji [21] derives a model for by doubly integrating the homogeneous mixture fraction PDF transport equation. His formulation is based on the observation that a presumed -PDF is capable of accurately characterizing the evolution of the scalar PDF over all stages of two-scalar, constant-density mixing in statistically stationary, isotropic turbulence. The model is given by: where is the Favre-averaged scalar dissipation rate (discussed in more detail in Section 3.3) and is given by the integral: As mentioned above, Girimaji’s model is valid for homogeneous turbulence. Therefore, its usage in the CMC of inhomogeneous flows is inconsistent and results in spurious source terms.

Mortensen’s Model. In a similar fashion, Mortensen [22] derives an expression for by doubly integrating the mixture fraction PDF transport equation. However, in his derivation, the inhomogeneous terms are retained. The CV fluctuations appearing in the transport equation are modelled using the PDF gradient model of Pope [19] and the PDF is presumed using a functional form described by the vector of mixture fraction moments. When the -PDF is employed, this model reads as follows [23]: where is the incomplete -function. By applying some identities and integral properties of the incomplete -function, Mortensen [23] shows that simplifies to The term in (17) represents the source term of the mixture fraction variance transport equation (discussed in more detail in Section 3.3). It is given by Mortensen’s model ensures a fully consistent CMC implementation as it accounts for the inhomogeneous terms of the PDF transport equation and employs a consistent closure for the CV fluctuations. As such, this model does not yield any spurious sources. It is important to note that Girimaji’s model is the exact equivalent of the homogeneous portion of Mortensen’s, .

3.2.5. Conditional Chemical Source

The conditional chemical source term, , is modelled using a first-order closure. In this closure, it is assumed that the conditional fluctuations about the conditional averages of the reactive scalars are small. Accordingly, is modelled as a function of the conditional density, mass fractions, and temperature: where and . This closure has been successfully applied in the CMC of lifted flames [7, 8, 11, 24].

3.2.6. Conditional Radiative Source

The conditional radiative source term is modelled using the optically thin assumption with being the predominantly participating species. To a first order approximation, where  W/m2K4 is the Stefan-Boltzmann constant, and are the partial pressure, and the Planck mean absorption coefficient of , respectively, and is the background temperature (300 K). A Curve fit for is obtained from [25].

3.3. Turbulent Flow Field Calculations

The Favre averages of the velocity, mixture fraction mean and variance, turbulence kinetic energy, and turbulence eddy dissipation are obtained from an inert flow calculation using the commercial Computational Fluid Dynamics (CFD) software FLUENT [26]. This approach was successfully applied in the CMC calculations of lifted flames performed by Kim and Mastorakos [24]. They justify the validity of this simplification by the fact that the flow remains frozen until the stabilisation height. Nevertheless, they note that density changes caused by the flame are non-negligible and may affect the flow field.

A two-dimensional axisymmetric computational domain is employed to perform the inert flow calculations. The domain extends 55 nozzle diameters above the nozzle exit in the axial direction and 15 diameters in the radial direction with the origin of the domain placed at the centre of the nozzle. The thickness of the wall of the nozzle is neglected [3, 27]. A quadratic unstructured mesh is generated using ANSYS ICEM CFD [28]. It consists of nodes (axial radial). The model is employed to perform the calculations. Default model constants are employed, except for which is modified from the standard value of 1.44 to 1.6 in order to improve the prediction of the spreading rate of the jet [29]. The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm is used for pressure-velocity coupling. The PREssure STaggering Option (PRESTO) pressure interpolation scheme is employed to compute the pressure at cell faces. Spatial discretization is performed using the second-order upwind scheme. The inlet velocity boundary condition is specified using the 1/7 power law with the centreline velocity set to 5/4 times the experimentally estimated mean value of 107 m/s as in [27]. The coflow velocity is assumed to be uniform and is set to 3.5 m/s. The turbulence intensity is taken to be 5% at both the inlet and the coflow. The boundary conditions of the temperature and species mass fractions are specified following Table 1. Transport equations for the mean and variance of the mixture fraction are added to the solver. These are given by The gradient diffusion assumption is employed to close the turbulent fluxes as and . The last term on the r.h.s. of (25), , is the Favre-averaged scalar dissipation rate. It is a standard practice to model this term as where is a constant set to 2. In this algebraic model, it is assumed that the flow time ( ) is proportional to the time scale of scalar turbulence ( ), with being the proportionality constant (the time scale ratio). This modelling approach was adopted in the early stages of this study. However, it was abandoned later because the calculation yielded very short liftoff heights for moderate and low and nearly attached flames for sufficiently high . This behavior was not observed in the lifted flames simulations of Patwardhan et al. [11] and Kim and Mastorakos [24]. Alternatively, the transport equation derived for by Jones and Musonge [30] is solved. This equation reads where , , , , and are model constants. The turbulent flux is modelled using the gradient diffusion assumption. For consistency with the Boussinesq hypothesis, the production of turbulence kinetic energy appearing in the last term on the r.h.s of (27) is modelled as where is the turbulent viscosity and is the modulus of the mean strain rate tensor. The tensor appearing in the first term on the r.h.s. of (27) represents the anisotropic diffusivity. Initial trials showed numerical instabilities and suffered from severe convergence issues when accounting for anisotropic effects. Therefore, for simplicity, the diffusivity was assumed to be isotropic and standard modelling was employed ( used instead of ). The inclusion of the transport equation of (see (27)) in the context of CMC was first employed in the counterflow flames calculations of Kim and Mastorakos [31]. The authors solve an additional transport equation for rather than resorting to gradient diffusion modelling.

3.4. The CMC Implementation

A smaller physical domain is employed in the CMC calculations for computational efficiency. The chosen domain is 30 nozzle diameters in length and 5 diameters in width. The mesh consists of (axial radial) nonuniformly distributed nodes, with the mesh density being highest near the experimentally measured stabilization height ( ) [3]. The mixture fraction grid consists of 80 points. These resolutions ensure grid independence. The CMC equations are discretized using the finite difference method. The first-order derivative appearing in the convective term is discretized using the second-order upwind difference scheme with the kappa flux limiter of Koren [32], following Patwardhan et al. [11]. The second-order derivatives in physical and mixture fraction spaces are discretized using the second-order central difference scheme. The steady-state solution is obtained from the transient CMC transport equations, (8) and (9), by relaxation of the time step, for which a value of 10−5 s is used. A three-step fractional method (Strang operator splitting [33]) is implemented in order to treat the stiff chemical source term separately. The first fractional step accounts for convection and diffusion in physical space over the first half of the time step. The second step handles the chemical source and micromixing in mixture fractions space over the whole time step. The third and last step accounts for convection and diffusion in physical space over the second half of the time step. Each step uses the solution of the previous one as initial conditions. Splitting errors are negligible due to the usage of a small time step. The Ordinary Differential Equation (ODE) solver VODPK [3436] is used to solve the system of equations. The ODEs in the first and third steps are nonstiff and hence solved with the implicit Adams method. Conversely, the ODEs in the second step are stiff and therefore treated with the Backward Differentiation Formulas (BDF). The boundary conditions for and at the inlet and the coflow are specified using the inert mixing solution. The composition of the fuel and coflow streams is obtained from the experiments [3]. Zero-gradient boundary conditions are specified at the remaining boundaries. The solution is initialized by igniting some nodes above the experimentally measured stabilization height with a steady laminar flamelet computed with a strain rate of 1000 s−1. The Inert mixing solution is set elsewhere. The Favre-averaged velocity, mixture fraction mean and variance, turbulence kinetic energy, and turbulence eddy dissipation are transferred from FLUENT to the CMC solver by applying bilinear interpolation. This interpolation scheme is accurate enough given the high resolution of the CFD grid. The integrals appearing in Girimaji’s model are computed using the quadrature integration package QUADPACK [37]. In Mortensen’s model, the incomplete beta-function (see (19)) is calculated using the method of continued fractions [38] and the partial derivatives of (see (20)) with respect to the mean and variance of the mixture fraction are computed using Ridders’ method of polynomial extrapolations [38, 39]. Due to the division by the PDF in the PDF-gradient and Mortensen’s models, numerical instabilities may arise at low probabilities. To resolve this issue, is set equal to its unconditional value, , and is set to zero when falls outside . Additionally, in the event where Moretnsen’s model yields unphysical negative values, is set to zero.

The chemical kinetics mechanisms developed by Mueller et al. [40], Li et al. [41], Burke et al. [42], and Ó Conaire et al. [43] are considered in this study. The corresponding numbers of species and reactions are summarised in Table 2. All four mechanisms have a comparable level of detail, with those of Mueller and Burke being the least and most detailed, respectively. The Mueller mechanism is the most extensively used throughout this work. It is provided in Table 3.

4. Results and Discussion

4.1. Flow Field Results

The radial profiles of the Favre-averaged axial velocity are compared to the experimental measurements of Kent [44] in Figure 2(a). The predictions are in very good agreement with the experimental data up to . The measurements are slightly underpredicted at for . No velocity measurements are available at . The profile at this location is provided for completeness. The radial profiles of the Favre-averaged mixture fraction mean are displayed in Figure 2(b) along with the experimental data of Cabra et al. [3]. The predictions are in very good agreement with the experiments at = 1, 8, 9, 10, and 11. Slight underpredictions are obtained at and 26 for . Figure 2(c) shows a comparison of the radial Favre-averaged mixture fraction variance profiles with the experimental measurements. A very good agreement is obtained at all axial locations except at . The discrepancies at this location are attributed to the underprediction of the mixture fraction mean (bottom pane in Figure 2(b)). The results obtained with the algebraic modelling of (see (26)) tend to be underpredictive near the centreline and overpredictive around the peaks of the experimental data (not shown here). Therefore, it is concluded that the alternative approach of incorporating an additional transport equation for (see (27)) provides more accurate mixture fraction variance predictions. Overall, the results presented in Figure 2 show that the inert flow calculations employed in this study provide reliable flow and mixing fields for the CMC calculations.

The centreline evolution of the time scale ratio, , is shown in Figure 3. Here, is computed using (26) with obtained from the solution of (27). In the vicinity of the nozzle, is well above the commonly used constant value of 2. Away from the nozzle, decays gradually and tends asymptotically to 2.76. This suggests that the usage of (26) with would generally underpredict the centreline scalar dissipation levels, especially near the nozzle area. The trend obtained in Figure 3 is in line with the experimental observations of Markides [45] who performs scalar dissipation measurements in autoignitive nitrogen-diluted fuel jets injected into coflowing heated air, a configuration similar to the one under investigation. His results show that the centreline values of decay from high levels near the injector to a constant value around 2 at downstream locations. Accordingly, in a subsequent CMC study in the same context, Markides et al. [46] tune in (26) in the vicinity of the injector in order to emulate the decay of the time scale ratio and set it to 2 downstream. In the absence of scalar dissipation measurements, as in the present case, such treatment is ad hoc and tedious: must be tuned by trial and error until a reasonable agreement between the predicted and measured profiles is achieved, and most likely, a nonconstant functional form for may be necessary. This process is easily avoided by means of solving (27).

4.2. Effect of the CSDR Modelling for Various Coflow Temperatures

The effect of the modelling of the CSDR is investigated using the homogeneous model of Girimaji [21] and the inhomogeneous model of Mortensen [22]. The reader is reminded that the latter model degenerates to the former when the inhomogeneous terms in (17) are discarded. Calculations are performed for three coflow temperatures. In addition to the experimentally reported value  K [3], two other temperatures are investigated: and 1060 K. Given the experimental uncertainty of 3% in the temperature measurements [3], these two values lie within the experimental error. The chemical kinetics mechanism of Mueller et al. [40] is employed throughout this section.

4.2.1. Results in Physical Space

The radial profiles of the Favre-averaged reactive scalars are shown in Figures 4(a)4(e) at the axial locations , 9, 10, 11, 14, and 26. For  K (thick black lines), the results obtained using the two CSDR models are almost identical, except for some minor differences at and 11. When is decreased to 1030 K (thin black lines), the predictions remain in good agreement with the experiments up to and the two CSDR models yield very similar results. Farther downstream, the temperature and products mass fractions are underpredicted and the reactants mass fractions are overpredicted, most notably at and 14. The influence of the CSDR modelling becomes more apparent at these axial locations. Overall, the results remain in reasonable agreement with the experiments. When is increased to 1060 K (grey lines), the CSDR models yield distinct results for . Using both models, the temperature and products mass fractions are grossly overpredicted and the reactants mass fractions are underpredicted. All three coflow temperatures result in very similar profiles at , irrespective of the CSDR model. As a general observation in Figure 4, the usage of Girimaji’s model always yields lower temperatures and product mass fractions and higher reactant mass fractions when compared to Mortensen’s, which leads to a shorter liftoff height. This is attributed to the fact that Girimaji’s model results in earlier ignition delays and hence shorter ignition kernels for all coflow temperatures. This phenomenon will be discussed shortly in more detail. Before doing so, the liftoff heights obtained from the cases considered in Figure 4 are first presented. The contours of are shown in Figure 5(a) for  K, Figure 5(b) for  K, and Figure 5(c) for  K. In each subfigure, the left and right panes correspond to the results obtained with Mortensen’s and Girimaji’s models, respectively. As shown in Figure 5(a), stabilization takes place at slightly distinct locations due to the very small differentials (see the radial profiles in Figure 4(e) for  K). Mortensen’s model (left) results in (19.1% in excess of ) while Girimaji’s model (right) gives (14.8% in excess of ). The same scenario is observed in Figure 5(b) where  K. However, in this case, the location of stabilisation is better predicted. Mortensen’s model (left) yields (6.1% in excess of ) while Girimaji’s model (right) results in (5% in excess of ). The advantage of the inclusion of the effects of inhomogeneity Mortensen’s model is judged to be insignificant for these two coflow temperatures. Conversely, the differences in the stabilisation locations are more substantial for  K, as displayed in Figure 5(c). This is due to the much larger differentials (see the radial profiles in Figure 4(e) for  K). Mortensen’s model (left) yields (26.1% below ) while Girimaji’s model (right) gives (33.9% below ). As for the radial location of stabilisation (the radial distance from the stabilisation point to the centreline normalised , denoted here as ), the values obtained using the two CSDR models for  K and 1045 K are the same (1.58 and 1.47, resp.). As for  K, the two models yield slightly different values (Mortensen’s, 1.26; Girimaji’s, 1.16). It is obvious from the above results that the liftoff height becomes smaller and that flame base becomes narrower as is decreased. This is due to the fact that at lower , the mixture can autoignite in areas closer to the nozzle which are characterised by high scalar dissipation rate levels. This behaviour will be addressed in more detail in Section 4.2.3.

4.2.2. Stabilisation Mechanism

The numerical indicators developed by Gordon et al. [6] for the distinction between stabilisation by autoignition and premixed flame propagation are employed here to determine the underlying controlling mechanism. For this purpose, the transport budget of the temperature and the history of radical build-up ahead of the stabilisation height are analysed.

Budgets in Mixture Fraction Space. Figure 6 shows the transport budget of the steady-state conditional temperature equation, that is the individual contributions of the r.h.s. terms of (9). Only the results obtained using Girimaji’s model are presented here. Mortensen’s model produces similar trends. For each , the budgets are displayed at several axial locations around the corresponding liftoff height. The radial locations are held fixed and set to the value obtained with each . As such, the chosen coordinates cover the preflame, liftoff, and postflame regions. Before analysing the balance in these budgets, the axial evolution of the chemical source term, , is first examined. As shown in Figures 6(a), 6(b), and 6(c), common trends are observed in the evolution of this term for all three coflow temperatures. Below the stabilisation height (first two locations in the preflame regions, Pre-F1 and Pre-F2), as increases, shifts from very lean to less lean mixtures and its amplitude increases dramatically. At liftoff (middle panes), peaks at with a significantly higher amplitude. Within approximately two to three nozzle diameters downstream the stabilisation height (first locations in the postflame regions, Post-F1) the peak of increases further and occurs around the stoichiometric mixture fraction. Further downstream (second locations in the postflame regions, Post-F2), the magnitude of decreases substantially from the Post-F1 locations. The reaction zone becomes wider and presents two peaks. The first occurs around the stoichiometric mixture fraction and the second in rich mixtures. The second peak is most likely attributed to the propagation of a rich reaction zone towards stoichiometric mixtures [24]. All of the above trends are consistent with the direct numerical simulation Yoo et al. [47]. As in the evolution of , common trends are observed in the axial variation of the remaining terms contributing to the r.h.s. of (9). In the preflame region (locations Pre-F1 and Pre-F2), there is a clear balance in lean mixtures between , the axial convection term, , and micromixing, . The axial and radial diffusion terms, and , and the radial convection term, , are nonnegligible but have little contribution to the overall budget. This balance suggests that stabilisation by premixed flame propagation wherein balances in the preheat zone with negligible is not the case. Instead, stabilisation occurs via autoignition as can be seen from the - - balance. At liftoff, the balance shifts to . The role of in balancing diminishes but this term remains a major contributor to the overall balance, and prevails. As for the remaining terms, becomes much weaker while and remain important. In the postflame regions (locations Post-F1 and Post-F2), the role of diminishes further, becomes more dominant, and is virtually null. The contributions of and are as significant as that of at Post-F1 and supersede it at Post-F2. As increases in mixture fraction space, the contributions of , , and remain smaller than that of , which acts as the major heat sink. In contrast, as decays, increases and acts as a source (notice its positive contribution), and both of these terms are counterbalanced by , , and . Thus, beyond the stabilisation height, the flame budgets indicate the structure of a nonpremixed flame, which is largely characterised by a - balance. To be noted that the radiative source, , is negligible at all locations for all coflow temperatures due to the fact that hydrogen in the fuel stream is highly diluted with nitrogen. Therefore, this term can be safely neglected without loss of accuracy.

Budgets in Physical Space. The fact that autoignition is the controlling stabilisation mechanism can be further confirmed by analysing the transport budget of in physical space. To do so, the PDF-weighted integration of the r.h.s. terms of (9) is performed over the mixture fractions space. This yields , , , , , , and (“IT” stands for “Integrated Transport”). The axial profiles of these terms are plotted in Figure 7(a) for Girimaji’s model and in Figure 7(b) for Mortensen’s model at the respective locations obtained with each . The parallels between the observation made in Figure 6 and this figure are evident. is mostly balanced by and in the preflame regions, with the former being more dominant. The remaining terms are small but nonnegligible. The absence of eliminates the possibility of stabilisation by premixed flame propagation and the current balance indicates that autoignition is the controlling mechanism. Right ahead of liftoff, and emerge and becomes more important. These three terms start to dominate at the expense of , which remains the main heat sink. Downstream of the stabilisation height, , , and diminish gradually and becomes more important. Further downstream is mostly counterbalanced by , with smaller contributions , , and , which indicates the structure of a nonpremixed flame.

Radical History ahead of the Stabilisation Height. The transport budgets of the temperature in mixture fraction and physical spaces are both indicative of stabilisation via autoignition. In order to ascertain this conclusion, the history of radical build-up ahead of the reaction zone is now investigated. Figure 8 displays the axial profiles of the normalised Favre-averaged temperature ( ) and mass fractions ( ) of , , , and . The subscripts “min” and “max” denote the minimum and maximum values of the reactive scalars at . The results are reported for Girimaji’s and Mortensen’s models in Figures 8(a) and 8(b), respectively. Common trends are observed ahead of the stabilisation height for the different combinations of values and CSDR models.(1)The production of the intermediate initiates upstream of the axial locations where the radicals , , and emerge.(2) , , and start building up before the consumption of begins. As shown, the mass fractions of these three species start to increase before reaches its peak.(3) builds up rapidly prior to the runaway of , , and , as can be seen from the axial locations of the maximum slopes.

Therefore, radicals do not build up simultaneously as in premixed flame propagation, but rather acts as a precursor to the production of , , and as in autoignition scenarios [6]. This confirms that the flame is stabilised by autoignition.

4.2.3. Impact of CSDR Modelling on Autoignition

To better understand the autoignition process, standalone zero-dimensional CMC (0DCMC) calculations are performed for each . In these calculations, the convective and (physical space) diffusive terms in (8) and (9) are turned off and the CSDR is modelled using the AMC [9]. This CSDR model is parametrised by its maximum value and it is given by where is the inverse error function. The resulting homogeneous CMC equations are identical to those of the unsteady laminar flamelet model (with unity Lewis number) [48]. For a given , the solution evolves from the inert mixing solution at up to the moment when autoignition occurs, . Autoignition is declared at the moment when the maximum OH mass fraction reaches at any point in mixture fraction space, following Stanković and Merci [49]. Figure 9 displays the ignition delay ( ) as a function of . These results are obtained using the mechanism of Mueller et al. [40] for  K, 1045 K, and 1060 K. It is evident from this figure that, for a given coflow temperature, reaches an asymptotic limit as is gradually increased. This means that the occurrence of autoignition becomes impossible beyond a critical value of , denoted here by . The approximate values of are 227 s−1, 482 s−1, and 991 s−1 for  K, 1045 K, and 1060 K, respectively. Therefore, increases with increasing . In order to link these findings to the flow field of the flame under investigation, is first obtained by integrating (29) weighted by over the mixture fraction space: which requires the usage of the readily available quantities , , and . Figure 10 shows the Contours of computed from (30) for  K. The contours for  K and 1060 K are virtually the same and are not included here for brevity. The contour of Figure 10 will be used in all of the subsequent analysis irrespective of the value of . The thick contours represent the levels 227 s−1, 482 s−1, and 991 s−1 which correspond to the values determined in Figure 9. Also shown are the ignition kernel locations obtained from the different combinations of CSDR models and coflow temperatures. As expected, for all three coflow temperatures, ignition occurs at locations where is lower that (the circles, squares, and triangles lie outside the 227 s−1, 482 s−1, and 991 s−1 contour levels, resp.). Since decays away from the nozzle, autoignition does not take place until (which depends on via (30)) falls below . This explains why increasing results in earlier ignition and shorter liftoff from a mixing point of view.

In another set of 0DCMC calculations, the “most reactive” mixture fraction, , is computed for each by setting in (29) to 10−20 s−1 (homogeneous reactor). The results are summarised in Table 4. Consistent with the findings of Stanković and Merci [49], decreases with increasing . The objective behind determining is to compare the local CSDR values obtained using Mortensen’s and Girimaji’s models at , and , respectively. As shown in Figure 11, for all given coflow temperatures and radial locations, the axial profiles of and reveal that Mortensen’s models generally result in higher scalar dissipation rate levels and therefore yield longer ignition delays and higher liftoff heights. Also displayed in Figure 11 are the critical CSDR values computed using the AMC model parametrised by at ; that is, (horizontal grey lines). These values are also provided in Table 4 for all coflow temperatures. As expected, increases with increasing . Figure 11(a) (  K) shows that ignition cannot occur as long as or is below . Indeed, when Mortensen’s model is used, ignition happens at the first axial and radial locations where falls below . This takes place at and (first pane from the top in Figure 11(a)). A similar scenario is observed when Girimaji’s model is used. Ignition occurs at the first axial and radial locations where falls below , which happens at and (second pane from the top in Figure 11(a)). Similar observation can be made in Figure 11(b) (  K) and Figure 11(c) (  K). However, in these two cases, although and fall quickly below at locations smaller than the ones where ignition is indicated, ignition occurs at the first axial location with lowest . As shown in Figure 11(b) (  K), ignition takes place at and using Mortensen’s model and at and using Girimaji’s model. Similarly, in Figure 11(c) (  K), ignition happens at and using Mortensen’s model and at and using Girimaji’s model.

4.2.4. Effect of Spurious Sources

Ideally, the integration of the PDF-weighted CMC equations over -space should yield the unconditional set of equations without any additional source terms [14]. This outcome is guaranteed when fully consistent CSDR models (e.g., the inhomogeneous model of Mortensen) are employed. To be noted that in order for an inhomogeneous CSDR model to be fully consistent, the CV model employed in the derivation should be consistent with the moments of the presumed PDF. This condition is satisfied in Moretnsen’s model due to the usage of the PDF-gradient model. An example of an inhomogeneous yet inconsistent CSDR closure is the one proposed in [50]. The discrepancy arises from the modelling of the CV using the linear velocity model, which is not consistent with the second moment of the presumed -PDF.

When inconsistent CSDR models are employed, the integration of the CMC equations results in spurious (or false) source terms. In this study, this occurs when the CSDR model of Girimaji is used. The spurious source associated with a species is calculated through [14, 16]: where is the conditional mass fraction of species obtained in the inconsistent realisation (using Girimaji’s model).

Figure 12 shows the radial profiles of the mean chemical sources ( ) and spurious sources of , , , and for  K. In comparison to the mean chemical sources, the magnitudes of the spurious sources are small but nonnegligible, in particular within the flame zone ( and 26). The effect of the spurious sources ahead of the predicted stabilisation heights ( and 10) is not sufficiently large to change the nature of the stabilisation mechanism, as the flame is stabilised by autoignition using both Mortensen’s and Girimaji’s models. However, the fact that Girimaji’s model yields a relatively earlier ignition and a smaller liftoff height indicates that the spurious sources are influential. Although this is not the case in the calculations of the current flame, inconsistent CMC implementations may lead erroneous results and conclusions.

4.3. Effect of Chemical Kinetics

The mechanism of Mueller et al. [40] has been used in all of the results reported so far. The sensitivity of the flame to chemical kinetics is now assessed by testing the performance of the mechanisms of Li et al. [41], Burke et al. [42], and Ó Conaire et al. [43]. The coflow temperature is set to the experimentally reported value of 1045 K.

4.3.1. Comparison of the Chemical Kinetic Mechanisms

0DCMC A Priori Analysis. 0DCMC calculations are first performed using the three mechanisms in order to determine the corresponding values. The results are shown in Figure 13, along with those obtained with the Mueller mechanism in Section 4.2.3 for K. By inspecting the calculated values and referring to Figure 10, it can be postulated that the Li mechanism would yield the shortest ignition kernel and smallest liftoff height followed by the Mueller mechanism, then the Ó Conaire mechanism, and finally the Burke mechanism. This analysis provides a preliminary qualitative description of the sensitivity of liftoff height to the considered chemical kinetics.

Results in Physical Space. To ascertain these findings, the flame is simulated using the three newly considered mechanisms. The CSDR is closed using Girimaji’s model in all of the following calculations. The results are displayed in Figure 14 at the axial locations = 9, 10, 11, 14, and 26, along with those obtained using Mueller mechanism in Section 4.2.1. Unlike the Burke mechanism, the other three kinetic schemes are capable of reproducing the experimental trends. The results of the Mueller and Ó Conaire mechanisms show best agreement with the experimental data. The Mueller mechanism yields better predictions at and 10, while Ó Conaire’s results in better agreement with the measurements at . The differences between the two are small at and 26. When the Li mechanism is employed, the temperature rises at early axial locations owing to the ability of the mixture to ignite at high scalar dissipation rate levels (  s−1). This effect propagates downstream and results in the deviation of the numerical predictions from the experiments. The performance of the Burke mechanisms is unsatisfactory, as the mixture fails to achieve early ignition ( s−1) and remains inert up to approximately . The contours of ( ) are displayed in Figure 15. As predicted from the 0DCMC calculations, the Li and Burke mechanisms yield the smallest ( ) and largest ( ) liftoff heights, as shown in Figures 15(a) and 15(b), respectively. Also being consistent with the 0DCMC findings, the liftoff heights obtained using the Mueller and Ó Conaire mechanisms fall between those of the Li and Burke mechanisms, as demonstrated in Figures 15(c) and 15(d). However, the Mueller mechanism results in a slightly larger liftoff height ( ) in comparison to the Ó Conaire mechanism ( ), which contradicts with the 0DCMC predictions. Apart from the Burke mechanism, the liftoff heights computed from the remaining mechanisms are similar and are in close agreement with the experimental value . As for the radial location of liftoff, the Mueller, Li, and Ó Conaire mechanisms result is , whereas Burke’s yields a wider flame base with . As in the Mueller mechanism calculations, the flame is also stabilised by autoignition when the Ó Conaire, Li, and Burke mechanisms are employed. This can be clearly seen in Figure 16 where , , , and exhibit the same trends observed in Figure 8(a).

Results in Mixture Fraction Space. For completeness, the conditional profiles are shown in Figure 17. The calculation of the conditional data from the experimental scatter is described in the Appendix. The numerical results are reported at = 9, 10, 11, 14, and 26, and ( obtained with all mechanisms except Burke’s). Overall, the results of the Mueller and Ó Conaire mechanisms are in reasonable agreement with the experiments, especially within mixture fraction intervals where the PDF is large (see Figure 17(f)). This explains the good agreement these two mechanisms yield with the measurements in physical space (Figure 14). The Li mechanism predicts the highest levels up to (Figure 17(e)) which leads to the occurrence of the earliest ignition and consequently results in the smallest liftoff height among the considered kinetic schemes. In turn, this early ignition produces higher than expected temperature levels, an effect which propagates downstream resulting in the deviation of the conditional predictions from the experiments. In the three mechanisms discussed up to this point, the peak of in Figure 17(a) shifts gradually from lean to stoichiometric mixtures as increases, which is consistent with the axial evolution of the conditional heat release term in Figure 6. The profiles of the Burke mechanism remain very close to the inert mixing solution up to . This behaviour results in the most delayed ignition and leads to the largest liftoff height. The conditional profiles at where liftoff occurs undergo very similar trends (not shown here).

4.3.2. Diagnosis of the Burke Mechanism

In order to diagnose the unsatisfactory performance of the Burke mechanism, the rate parameters of two of the key reactions are investigated. The recombination reaction and the chain-branching reaction are two of the most important reactions in oxidation chemistry. and (labelled here as in [4043]) compete for atoms and therefore control the overall branching ratio and determine the second explosion limit in homogeneous systems [42]. The treatment of and in the Li and Burke mechanisms differs significantly from the Ó Conaire and Mueller mechanisms. The rate parameters of are provided in Table 5. The high-pressure limit rate parameters are identical in the Ó Conaire, Mueller, and Li mechanisms [51]. Burke’s mechanism uses an updated set of parameters [52]. In the low-pressure limit, the Ó Conaire and Mueller mechanisms use the same parameters [53] except for the third-body efficiency factors of and . The Li mechanism employs a completely different set of parameters [54]. The Burke parameters differ from Li’s by the broadening factor (decreased by a factor of 1.6 and equal to the ones employed in the Ó Conaire and Mueller mechanisms) and the third-body efficiency factor for (increased by a factor of 1.3) [42]. As for , the corresponding rate parameters are listed in Table 6. The parameters employed in the Ó Conaire and Mueller mechanisms are identical [55]. Two distinct sets of updated parameters are employed in the mechanisms of Li [56] and Burke [57]. In light of the above description, the variability of the predictions in Figure 14 and the close agreement between the results of the Mueller and Ó Conaire mechanisms may be explained in part by the treatment of and . This can be demonstrated, for instance, by adjusting the rate parameters of these two reactions in the Burke mechanism, in an attempt to improve its predictions. Since this mechanism is an updated version of Li’s, two separate experiments are carried out as follows. (1)The parameters of   are substituted by those of the Li mechanism, with the parameters of   left unchanged.(2)The parameters of both and are substituted by those of the Li mechanism.

The preliminary 0DCMC calculations shown in Figure 18 reveal that increases significantly when these modifications are applied. Therefore, the mixture is now anticipated to ignite earlier and the liftoff height is expected to decrease. Two additional CMC runs are performed using these modifications while keeping the same settings for the CSDR and the coflow temperature. The results are shown in Figure 19 at = 10, 11, and 14. It is evident that the modification of alone brings a substantial improvement over the rate parameters of the original mechanism. The temperature (Figure 19(a)) and the species mass fractions (Figures 19(b)19(e)) are much better predicted, especially at and 14. This further proves the importance of as a precursor to autoignition reactions. Better agreement with the experimental data is obtained when both and are modified. As shown in Figure 19(f), the liftoff height decreases by 26.69% when is modified and by 29.24% when both and are modified. In both cases the computed liftoff heights are substantially closer to the experimentally measured value. Figure 19(f) also shows a remarkable decrease in as liftoff occurs closer to the centreline. With the modified settings, both the axial and radial positions of liftoff become in line with those obtained using the Ó Conaire, Mueller, and Li mechanisms. There should be further room for improvement by applying additional modifications to the rate parameters of other key reactions. However, extreme care must be taken while doing so because chemical kinetic mechanisms are usually optimised as a whole, rather than on an individual reaction basis. A sensitivity analysis is beyond the scope of this work and was not performed here. The modifications applied above are ad hoc.

5. Conclusions

A lifted jet flame issuing into a vitiated coflow was investigated using the first-order CMC. The flow and mixing fields were obtained using a modified version of the model. The calculations included an additional transport equation for the Favre-averaged scalar dissipation rate. Two formulations were implemented for the CSDR, the homogeneous model of Girimaji, and the inhomogeneous model of Mortensen. The CV was modelled using the PDF-gradient model. Calculations were performed for three coflow temperatures, 1030, 1045, and 1060 K, and four chemical kinetic mechanisms were assessed. In light of the previous results, the following conclusions are drawn.(1)The modification of the constant in the model from the standard value of 1.42 to 1.6 improves the predictions of the spreading rate of the jet.(2)The solution of a transport equation for the scalar dissipation rate in lieu of using the traditional algebraic modelling approach provides a more reliable mixing field.(3)The flame is very sensitive to small changes in the coflow temperature. A variation of roughly in the experimentally reported value of 1045 K results in substantial changes in the predictions. K yields best agreement with the experiments. The calculated liftoff height at the 1030 K level remains close to the experimentally measured value. A drastic decrease occurs when the coflow temperature is increased to 1060 K.(4)For all coflow temperatures and using both CSDR models, the transport budgets in mixture fraction and physical spaces and the history of radical build-up ahead of the reaction zone indicate that the flame is stabilised by autoignition rather than premixed flame propagation. These findings fully agree with previous PDF calculations [46], with the CMC results of Stankovic and Merci [7], and in part with those of Navarro-Martinez and Kronenburg [8] and Patwardhan et al. [11].(5)Standalone zero-dimensional CMC calculation indicates that the mixture is capable of igniting at higher scalar dissipation levels as the coflow temperature is increased. This provides an explanation for the occurrence of ignition at locations closer to the nozzle exit and the decrease in liftoff height with increasing coflow temperature.(6)In comparison to Mortensen’s model, Girimaji’s model always results in lower CSDR at the “most reactive” mixture fraction and therefore results in earlier ignition and smaller liftoff heights for all coflow temperatures.(7)The spurious source terms resulting from the modelling of the CSDR are in general small but nonnegligible, mostly notably within the flame zone.(8)The flame shows high sensitivity to chemical kinetics. The results obtained with the Mueller, Li, and Ó Conaire mechanisms reproduce the experimental trends with varying degrees of accuracy. The predictions are in general in good agreement with the experimental measurements, particularly those of the Mueller and Ó Conaire mechanisms. The three kinetic schemes predict very similar liftoff heights. Conversely, the performance of the Burke mechanism is unsatisfactory due to the inability of the mixture to achieve early ignition. The temperature and species mass fractions are not well predicted and the liftoff height is grossly overestimated. The modification of the rate parameters of some of the key reactions improves the predictions substantially.

Further improvements are possible by modelling the conditional chemical sources with a second-order closure and by inclusion of counter-gradient and differential diffusion effects.

Appendix

Calculation of the Conditional Experimental Data

For each of the axial locations = 1, 8, 9, 10, 11, 14, and 26, measurements are reported at several radial locations. For a given axial location, the data is conditionally averaged at each radius as described in [14].(1)The range of the mixture fraction is subdivided into 40 bins.(2)The data is sorted into the bins.(3)The average of the data in each bin is taken to be the conditional average at the centre of the bin.

As an example, the radial profiles of the experimental conditional temperature at are shown in Figure 20(a). The weak radial dependence of the conditional data is obvious. Therefore, Steps 1–3 are repeated by using all of the experimental scatter at , irrespective of the radial location. The results are displayed in Figure 20(b).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This project is funded by the Ministry of Training, Colleges and Universities of Ontario through the OGS program, and the Natural Sciences and Engineering Research Council of Canada. The authors are grateful for their financial support.