Impact of molecular mixing and scalar dissipation rate closures on turbulent bluff-body flames with increasing local extinction

Combustion and Flame 206 (2019) 51–67

Contents lists available at ScienceDirect

Combustion and Flame journal homepage: www.elsevier.com/locate/combustflame

Impact of molecular mixing and scalar dissipation rate closures on turbulent bluff-body flames with increasing local extinction Lu Tian, R. Peter Lindstedt∗ Department of Mechanical Engineering, Imperial College, Exhibition Road, London SW7 2AZ, UK

a r t i c l e

i n f o

Article history: Received 31 October 2018 Revised 3 January 2019 Accepted 17 April 2019

Keywords: Hybrid transported PDF methods Mixing models Local extinction Scalar dissipation rate Generalised Langevin Model

a b s t r a c t Bluff-body turbulent CH4 : H2 (1:1) flames at 50% (HM1), 75% (HM2) and 91% (HM3) of the blow-off velocity (235 m s−1 ) were studied experimentally by Masri and co-workers and found to exhibit gradually increasing periodic and shear layer instabilities. The latter are coupled with increasing levels of local extinction with subsequent re-ignition further downstream. This study provides a systematic evaluation of the sensitivity of predictions to molecular mixing and scalar dissipation rate closures. The latter include extended forms of the Euclidean Minimum Spanning Tree (EMST) and modified Curl’s (MC) models, applicable to premixed turbulent flames via a closure that accounts for local Damköhler number effects (EEMST and EMC), and a conceptually related blended scalar time-scale approach (BEMST and BMC). Computations are performed using a hybrid finite volume (FV) – transported Joint Probability Density Function (JPDF) algorithm featuring stochastic Lagrangian particles, a comprehensive 48-scalar systematically reduced C/H/N/O mechanism, and a second moment method based on the Generalised Langevin Model that provides a partial resolution of the unsteady fluid motion. The sensitivity to solution parameters affecting the temporal resolution is quantified using Fourier transforms of the time histories of velocity and scalar traces. Radial profiles, conditional means and scatter plots are compared to the experimental data along with burning indices based on the conditional mean temperature. Vortex related instabilities ∼ 1 kHz in the outer shear layer appear for all closures with EMC showing periodic local extinction and re-ignition in the neck region for HM3 and flame turbules (i.e., discrete pockets of hot gas) separating periodically at frequencies ∼ 85 Hz. Results are similar to well–resolved JPDF/LES simulations for HM1. It is shown that the EMC and (E)EMST models essentially enclose the experimental data for HM2 and HM3. For HM3, emissions of NO are controlled by local extinction events that become increasingly sensitive to the molecular mixing closure as blow-off is approached. © 2019 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Multiple modes of turbulent combustion may occur simultaneously in practical combustion devices. Examples include piloted premixed gas turbines, diesel engines and industrial burners [1]. A combustion regime independent model is hence of great practical importance. The transported Joint Probability Density Function (JPDF) method [2] has the advantage that the chemical source term appears in closed form. It has seen wide application to non-premixed combustion [3] and is increasingly used for premixed [4] and partially premixed flames [5]. The closure of the molecular mixing term remains a primary challenge that includes the modelling of the scalar dissipation rate. The current work presents a systematic investigation of such closures using

∗

Corresponding author. E-mail address: [email protected] (R.P. Lindstedt).

bluff-body stabilised flames with increasing levels of local extinction. The flames feature hydrodynamically complex flow fields that contain a recirculation zone attached to the bluff-body, a neck zone where the interactions between turbulence and chemistry are intense, and a jet-flame-like core structure [6]. The flames are relevant to the aforementioned industrial applications and also appropriate for evaluating modelling approaches. Accurate treatment of turbulence–chemistry interactions is required for LargeEddy Simulations (LES) [7] and Reynolds-Averaged Navier–Stokes (RANS) methods. Chemical reactions typically occur on the unresolved scales and, for example, correlations with pressure fluctuations can result in significant cross-scale (i.e., backscatter) energy transfer [8] that need to be modelled, e.g., via moment based methods [9]. The latter can also offer an efficient partial resolution of the unsteady fluid motion, particularly if closed at the second moment level, and have been used to study non–linear hydrodynamic and thermoacoustic oscillations in bluff-body stabilised turbulent premixed flames [10].

https://doi.org/10.1016/j.combustflame.2019.04.039 0010-2180/© 2019 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

52

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

Nomenclature BI Cφ , Cφ∗

CSe , C 1 , C 2 Db Fh Fl Jα NTA Rejet Sα Stsl Ster T Ucof Ujet Y c f uK k uL

 c η

φ κ ωθ ρ θ δ τT

burning index [−] coefficients for scalar time–scale [−] coefficients in the dissipation rate equation [−] diameter of the burner bluff–body [m] high-frequency [Hz] low-frequency [Hz] molecular flux for species α [kg m−2 s−1 ] under–relaxation factor [-] Reynolds number based on the fuel jet diameter and velocity [−] reaction rate for species α [s−1 ] Strouhal number based on the shear layer instability [−] Strouhal number based on local extinction and re–ignition events [−] temperature [K] co–flow air velocity [m s−1 ] bulk velocity at the fuel jet exit [m s−1 ] mass fraction [−] reaction progress variable [−] mixture fraction [−] Kolmogorov velocity [m s−1 ] turbulent kinetic energy [m2 s−2 ] Laminar burning velocity [m s−1 ] turbulent kinetic energy dissipation rate [m2 s−3 ] scalar dissipation rate of reaction progress variable c [s−1 ] segregation factor in the blended scalar time-scale closure [−] scalar dissipation rate for φ [s−1 ] von Karman constant 0.41 [−] vorticity pointing towards the circumference direction [s−1 ] density [kg m−3 ] boundary layer momentum thickness [m] boundary layer thickness [m] turbulence integral time–scale [s]

Comprehensive experimental databases for the HM1–HM3 bluff-body flame series have been made available by Masri and co-workers [6,11–18]. These flames feature increasing levels of local extinction along with the existence of premixed gas pockets. Periodic and shear layer instabilities also become more severe as blow-off is approached. The data is of particular relevance given the current interest in hydrogen enrichment [19,20]. Flow structures and thermochemical characteristics were reported for flames at 50% (HM1) and 70% (HM2) and 91% (HM3) of the experimental blow-off velocity of 235 m s−1 . The HM3 flame features particularly strong turbulence–chemistry interactions and is the current focus. The HM flame series has been studied computationally using different levels of accuracy for the turbulence–chemistry interactions [3,21–28]. HM1 has the lowest probability of local extinction and can be reasonably reproduced by simple thermochemistry (e.g., steady-state flamelets and a first-order Conditional Moment Closure (CMC)) and different flow field closures [21–23,27]. CMC/LES has been shown to provide good agreement with experimental data for flames with moderate amount of local extinction [27,29]. Flames HM2 and HM3 have proved more challenging due to the enhanced probability of local extinction

τφ

Fφ



scalar time–scale of scalar φ [s] mass density function of scalar φ [−] sample space

Abbrevations BEMST EMST model coupled with blended scalar time–scale closure BMC MC model coupled with blended scalar time–scale closure CMC Conditional Moment Closure EEMST EMST model coupled with extended scalar time– scale closure EMC MC model coupled with extended scalar time–scale closure EMST Euclidean Minimum Spanning Tree GGDH Generalised Gradient Diffusion formulation of Daly and Harlow GLM Generalised Langevin Model IEM The Interaction by Exchange with the Mean model JPDF Joint Probability Density Function LES Large Eddy Simulation MC modified Curl’s model MMC Multiple Mapping Conditioning MPI Message Passing Interface MT-PRNG Mersenne Twister Pseudorandom Number Generator RANS Reynolds-Averaged Navier–Stokes RMS Root Mean Square URANS Unsteady Reynolds-Averaged Navier–Stokes Subscripts ad adiabatic flame properties lean parameters at the lean flammability limit max maximum rich parameters at the rich flammability limit s stoichiometric value u unburnt mixture properties Superscripts  Favre - averaged fluctuation – Reynolds - averaged mean ∼ Favre - averaged mean

and the co-existence of diffusion flames and premixed gas pockets. Multiple Mapping Conditioning (MMC) [30] has also been used for diffusion flames [31] with local extinction [32] and has recently extended to premixed combustion [33]. Transported JPDF method [24–26,28] can in principle predict flames with high levels of local extinction as well as re-ignition. The Interaction by Exchange with the Mean (IEM) [34], modified Curl’s (MC) [35] and Euclidean Minimum Spanning Tree (EMST) [36,37] mixing models have been widely used in transported JPDF methods and applied to the modelling of HM2 and HM3. Liu et al. [24] applied the JPDF of velocity–turbulence frequency–scalars to HM1–HM3 combined with EMST and with the turbulence closed at the second moment level. The predictions for HM1 showed essentially good agreement. However, the experimentally recorded local extinction was not observed for HM2 and HM3. Merci et al. [25] compared the performance of MC and EMST for the same set of flames using a JPDF and Reynolds-stress approach. Only results with EMST were presented for HM3 with MC leading to blow-off. Improvements for scalar statistics were obtained using a JPDF/LES approach [28] combined with the IEM model [34]. Species such as NO remained overestimated downstream of the recirculation zone.

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

The scalar time-scale closure is also vital for the accurate modelling of turbulent flames. The constant mechanical-to-scalar time-scale ratio approach, common in isothermal flows, becomes questionable in chemically reacting flows. Extended forms have been proposed for premixed turbulent flames, e.g., via an algebraic scalar dissipation rate closure that accounts for local Damköhler number effects [4]. Gkagkas et al. [26] combined a hybrid JPDF/unsteady RANS (URANS) approach [3] with the extended scalar time-scale closure and MC models (EMC) to compute flame HM2. Better agreement was observed with predictions of NO substantially improved. The EMC model was also found to improve predictions of turbulent burning velocities in premixed flames [4,38]. A related blended scalar time-scale closure [39] that accounts for deviations from bimodality has been recently proposed and coupled with EMST (BEMST) for the simulations of turbulent premixed flames. The flow dynamics of bluff-body flames has been discussed in previous studies. Kim and Pitsch [27] used CMC/LES to analyse HM1 in detail and showed that the flame is governed by largescale structures including vortex shedding, self-sustained oscillations and interactions of different mixing layers. The structures were found to be associated with locally high scalar dissipation rates, typically associated with local extinction, with reduced CO and NO formation in the neck region and also near the tip of the expanding main jet. Vortex shedding was also observed experimentally in HM2 and HM3 by Dally et al. [16]. The observations are consistent with Roquemore et al. [40] who observed flame turbules (i.e., discrete pockets of hot gas) periodically separating from a bluff–body stabilised flame using direct high-speed imaging and records of CH emission. Brum et al. [41] also noted the appearance of flame turbules associated with vortex shedding. The current work systematically evaluates the performance of the MC and EMST mixing models when combined with the extended and blended scalar time-scale closures, also applicable to premixed flames, to form the EMC, Blended MC (BMC), Extended EMST (EEMST) and BEMST models. The assessment of their applicability is hence of relevance to the formulation of combustion regime independent calculation methods for flows with extensive local extinction (HM3) treated using a hybrid JPDF/URANS approach closed at a comprehensive second moment level. Velocity and scalar statistics are presented along with flow dynamics. The sensitivity to solution parameters affecting the temporal resolution is quantified using Fourier transforms of time histories of velocity/scalar traces and explores the applicability of well-resolved second moment closures to unsteady combusting flows. Such methods (e.g., Speziale [42] and Frölich and von Terzi [43]) have, in addition to their applicability to flames undergoing thermoacoustic oscillations [10], been shown by Palkin et al. [44] to reproduce complex shedding flows with good accuracy at high Reynolds numbers where the (intractable) cross-correlations between coherent and stochastic stress tensor components vanish. 2. Methodology The current study is focused on flame HM3, which shows high levels of local extinction associated with strong finite-rate chemistry effects. The current hybrid JPDF/URANS approach [3,26] consists of a Lagrangian stochastic particle simulation of the jointcomposition-enthalpy PDF equation and a second moment closure for the velocity field formulated for compressible flows. Past work suggests that the latter provides reasonable accuracy for the current geometry [3,24–26] and remains computationally tractable with sophisticated chemistry. The Generalised Langevin Model (GLM) of Haworth and Pope [45,46] was adopted for the redistribution terms, while the Generalised Gradient Diffusion formulation of Daly and Harlow (GGDH) [47] was used for triple-moments.

53

The turbulent kinetic energy dissipation rate followed the standard form [48] with the C 1 = 1.8 [3,26] round jet modification. The governing equations are listed in Appendix A. The density was obtained from the mass density function Fφ in Eq. (1) using a Monte Carlo method [2] and coupled with the elliptic finite-volume based CFD solver [49], which provides the velocity field and time-scale information [2].

∂ Fφ ∂  ∂ + [u F ] + [S ( )Fφ ] = ∂t ∂ xl l φ ∂ α α  ∂     ∂  1 ∂ Jlα   − [ ul Fφ ] + Fφ ∂ xl ∂ α ρ ∂ x l 

(1)

The flow field was solved by a time-dependent, compressible, twodimensional axi-symmetric elliptic solver based on the variable predictor-corrector variant of PISO [50] algorithm with splitting error control and a TVD scheme [51]. The terms on the righthand side of Eq. (1) were modelled. The first term represents transport in physical space due to turbulent fluctuations and was closed via a gradient diffusion approximation with the turbulent Prandtl/Schmidt number set to 0.7. The second term stands for transport in scalar space due to molecular mixing and is a focal point of the current work as discussed below. The chemical reaction term is naturally closed. The applied systematically reduced chemistry is identical to that used in previous studies [3,26,52– 54] and features 20 solved species (H, O, OH, HO2 , H2 O, H2 , O2 , CH4 , CH3 , CO, CO2 , C2 H2 , C2 H4 , C2 H6 , N2 , N2 O, NO, NO2 , HCN and NH3 ) and 28 steady-state species (C, CH, 1 CH2 , 3 CH2 , CHO, CH2 O, CH2 OH, CH3 O, C2 , C2 H, C2 H3 , C2 H5 , C2 HO, C2 H2 O, N, NH, NH2 , N2 H2 , N2 H, HNO, HNO2 , CN, NCO, HOCN, HNCO, HCNO, HCN and H2 CN). 2.1. Closures for the scalar dissipation rate The standard scalar time-scale (τ φ ) closure is based on with a linear relationship between the mechanical (τ T ) and scalar (τ φ ) turbulence time-scales,

 φ Cφ   Cφ −1 τφ−1 =  = = τ  2 2 2 T  k φ

(2)

∂ φ  ∂ φ  where 

φ = ρ Dc ∂ x ∂ x /ρ¯ is the mean scalar dissipation rate; j

j

 φ  2 the mean scalar variance;  k and   the mean turbulent ki-

netic energy and its dissipation rate. The value of Cφ depends on the structure of the turbulence [3,25,26,55,56] and the applied model. Xu and Pope [55] suggested 1.5 ≤ Cφ ≤ 2.0 for EMST. Merci et al. [25] showed a low sensitivity with Cφ = 2.0 used in the present study. For the MC model, Cφ = 2.3 was preferred in studies of CH4 jet flames [56] and for HM1 [3] and HM2 [26] and is retained here. Flame HM3 features strong turbulence–chemistry interactions. Calculations of HM2 [26] have shown that the standard MC model tends to overestimate local extinction and that the EMC variant with an extended scalar time-scale closure that accounts for local Damköhler number effects, see Section 2.1.1, significantly improves agreement for scalars. Kuron et al. [39] proposed a related form with the aim to model the transition of the scalar time-scale from turbulence-dominated mixing to the flamelet regime using EMST as a basis. The current implementation of the blended scalar timescale closure is introduced in Section 2.1.2. 2.1.1. Extended scalar time-scale closure The standard expression for the time–scale ratio given in Eq. (2) can be derived from the transport equation for the scalar dissipation rate in passive flows with the assumption that the evolution

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

is dominated by local source terms [57]. An extension is required in high-Damköhler-number (Da) flows [39,57,58],

ρ¯  φ

=β

∗

ρ¯  φ 2



ρ¯  φ φ  2

=

Ct∗

60

(3)

1000

50 800 ]

where Sφ represents the reaction rate and Dφ the molecular diffusion coefficient. Mantel and Borghi [58] included a reaction specific term in the transport equation of order of O(Da), as shown in Eq. (3), and an expression was suggested based on the scalar dissipation rate and scalar variance as shown in Eq. (4),

Sφ

1200

40

1

 φ Cφ   ∂φ  ∂ Sφ τφ−1 =  = + where Sφ = 2Dφ 

2 2 k ρ¯ φ ∂ xl ∂ xl φ 

b

70

[s

Sφ

a

[c m s]

54

400 20 200

10

Sφ   2

600

30

(4)

ρ¯ φ

0

0 0.05

β∗

where ρ¯ stands for the mean density and with recently shown to be close to unity [9]. The second equality in Eq. (4) results from the linear relationship between the scalar dissipation rate and the mean reaction rate in the flamelet regime of combustion [59]. The constant Ct∗ is expected to be in the range 0.216 ≤ Ct∗ = β ∗ (2Cm − 1 )/2 ≤ 0.416 for a Cm variation from 0.716 to 0.916 consistent with BML theory [59]. Kuan et al. [57] closed Eq. (4) by using the fractalbased algebraic reaction-rate model [49] shown in Eq. (5).

  1  ) Sφ = CR ρu uL φ (1 − φ  k uK

(5)

where ρ u donotes the unburnt density; uL the unstrained laminar burning velocity; uK the Kolmogorov velocity. Bray et al. [60] evaluated a number of algebraic and transport equation closures for the rate of reaction and suggested that the majority failed to perform satisfactorily irrespective of the choice of modelling constants with Eq. (5) an exception. The closure has been successfully used in a number of premixed flame studies with CR ࣃ 3.0 [9,49]. The expression leads to,

Sφ  φ Cφ  Cφ   ρu u L   τφ−1 =  = + = 1.0 + Cφ∗  

2 2 ρ ¯  2 ρ ¯ u K  φ k k φ

(6)

Consistency with BML theory for the flamelet regime of combustion implies that Cφ∗ = 2Ct∗CR /Cφ has values in the range 0.648 ≤

Cφ∗ ≤ 1.248 with Cφ = 2.0. Cφ∗ = 1.2 [57] calibrated from the use of the MC model for molecular mixing is close to the upper limit. The implied higher value of Cm = 0.916 is in agreement with models for dilatation, scrambling and pressure transport in turbulent premixed flames [9]. Here, Cφ = 2.3 (compared to the standard value of Cφ = 2.0) is retained for consistency with the application of standard MC model and the previous studies of turbulent diffusion flames (e.g., [56]). The linear variation of the time-scale ratio of mechanical to scalar (τ T /τ φ ) turbulence with uL /uK is further consistent with the experimental findings of O’Young and Bilger [61]. The laminar burning velocity was obtained from a premixed laminar flame calculation with an unburnt mixture of 50% CH4 and 50% H2 . The value of uL varies with the mixture fraction and hence a β -PDF fit [53] was used to interpolate the laminar burning velocity in the flammable region as shown in Fig. 1(a). The extended scalar time-scale closure in Eq. (6) has been coupled with MC model in studies of premixed [38], non-premixed [53] and partially premixed [62] flames, and also provides a suitable scaling behaviour of turbulent burning velocities for flames in the flamelet regime of combustion, while being compliant with the standard approach for the limiting case of passive scalars. 2.1.2. Blended scalar time-scale closure An alternative form that also combines flamelet and turbulence mixing effects has been proposed on the basis of a detailed analysis of DNS data [39]. Rather than using the consumption rate in

0.1 [ ]

0.15

0

0.2

0.4 [

0.6 ]

0.8

1

Fig. 1. Profiles of (a) laminar burning velocities against mixture fraction ( f s = 0.05) and (b) conditional scalar dissipation rate against reaction progress variable obtained for a premixed stoichiometric CH4 :H2 (1:1) flame. Symbols: (◦) laminar flame calculations; lines: (—–) (a) a fit using a β -PDF (b) a fit using the 8th order polynomial.

laminar flamelets to limit the influence of the second term, an ad hoc segregation factor η, normalised by the bimodal variance, was used with the EMST mixing model and found to provide encouraging agreement. The blended scalar time-scale closure can be written as Eq. (7),

 φ

  2

φ

=

Cφ   (1 − η ) + η 2 k˜



1 0

c (ζ )dζ /c  2 < c |ζ > P

(7)

 where Cφ = 2.0; η = c 2 /[ c (1 −  c )];  c the scalar dissipate rate in laminar flamelets. Following the approach of Kuron et al. [39], c = Dc (dc/dx )(dc/dx ) was obtained by reconstructing the conditional scalar dissipation rate using a 1-D freely propagating premixed stoichiometric CH4 :H2 (1:1) flame [63]. As shown in Fig. 1(b), the data was fitted via a polynomial expression against the reaction progress variable. A consistent definition of the reaction progress variable [5] was used in the laminar and JPDF calculations,

c=

f −(YCH4 +YH2 ) fs

c ∈ (0, cmax )

(8)

where f and fs represent the mixture fraction and its stoichiometric value; YCH4 and YH2 are the mass fractions of CH4 and H2 ; cmax has a form of,

cmax =

f fs

cmax =

1− f 1− f s

f ≤ fs

f > fs

(9)

(10)

The mixture fraction f ≡ fs in the laminar case, but varies in the JPDF simulations. The extended closure in Eq. (6) is derived from the scalar dissipation rate for a reactive scalar and returns to the constant mechanical-to-scalar time-scale ratio for mixtures outside the flammability range (e.g., uL ( f ) = 0). The blended closure in Eq. (7) is a bimodal combination of passive and reactive scalar dissipation rates and depends on the computed shape of the blending coefficient (η). 3. Case configuration The bluff-body geometry consists of an axi-symmetric fuel jet with a 50:50 blend of CH4 :H2 , a diameter of D j = 3.6 mm and a bluff-body with an outer diameter of Db = 50 mm. Flame HM3 features a mean velocity of 214 m/s and a co-flow air with a mean

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

55

Table 1 Experimental configurations, where BO indicates the blow-off limit. Flame

Fuel

Ujet (m/s)

Ucof m/s

Rejet

% BO

HM3 HM3E HM1 HM2

H2 /CH4 (1:1) H2 /CNG (1:1) H2 /CH4 (1:1) H2 /CH4 (1:1)

214 195 118 178

40 35 40 40

28,700 26,300 15,800 23,900

91 90 50 75

velocity of 40 m/s. The Reynolds number based on the jet velocity and jet diameter is 28,700. Scalar fields were measured under these conditions with velocity fields obtained in a similar flame HM3E [6] at a slightly lower jet velocity of 195 m/s and a co-flow air velocity of 35 m/s. For flame HM3E, the methane content was substituted by CNG (91% of CH4 ). Conditions for flames HM1 and HM2 are also included for reference in Table 1. The axi-symmetric computational domain extends axially from the bluff-body surface to 160 mm downstream and radially from the centreline of the burner to 100 mm towards the co-flow air. The domain was discretised with 189 axial and 187 radial cells. The half width of the fuel jet was resolved with 12 cells and the area over the bluff-body surface by 97 cells. A coarser 189 × 138 grid was used in the simulations of flames HM1 and HM2 [3,26]. The refinement in the radial direction over the bluff-body (97 cells) for HM3 is significantly greater than that in previous JPDF/URANS simulations (60 cells) [24,25] and similar to the CMC/LES and JPDF/LES simulations (around 85 cells) [27,28]. Grid independence was verified through a finer grid featuring 282 (axial) × 267 (radial) cells. Boundary conditions are consistent with previous studies [3,26]. The turbulence is approximately fully developed in the fuel jet and co-flow at the burner exit [6] and the 1/7th power law was used for the inlet mean axial velocity profiles in the fuel jet and the co-flow boundary layer close to the bluff-body. The thickness of the latter was estimated as 3.75 mm from measurements [6]. The Reynolds stresses at the fuel jet inlet were set to follow a fully developed turbulent pipe flow [64] and the Reynolds stresses at the co-flow inlet were extrapolated from measurements of the HM3E flame [16]. A low sensitivity to boundary conditions in the co-flow was observed in the study of Kuan and Lindstedt [3]. The length scale at the burner exit was prescribed using the Prandtl’s mixing length hypothesis L = κ y where κ = 0.41 is the von Karman constant and y is the perpendicular distance from the bluff-body surface. The prescribed length scale affects the dominant frequencies [3] with the current configuration capturing the decay and spreading rates. Wall functions were applied at the solid boundaries and atmospheric pressure was specified along with transmissive (non-reflecting) outflow boundaries with a symmetry condition on the centreline. Different values for the under-relaxation factor NTA and time step were evaluated in order to assess the impact on the computed flow field dynamics. The base case features NTA = 400 and a time step t ≈ 10−6 s. Comparisons were made with either increasing NTA to 10 0 0 or decreasing t to 0.5 × 10−6 s. If the computed timedependent signals are evaluated according to the Nyquist criterion, the base case provides a resolved frequency of ࣃ 1250 Hz. The joint-composition-enthalpy PDF equation was solved using a Monte-Carlo method with 120 Lagrangian particles initialised in each finite volume cell [3,26] with results validated against 200 particles per cell. The stochastic particle simulations were performed in a parallel environment featuring 200 CPU cores and the Message Passing Interface (MPI) library [65]. The multistream Mersenne Twister Pseudorandom Number Generator (MT PRNG) [66] was used to generate multi-steam uncorrelated random numbers.

Fig. 2. Instantaneous contours of vorticity pointing to the circumferential direction for the low-frequency instability period obtained using the EMC model. The time inside the parentheses represents the actual computational time as shown in Fig. 5.

4. Results and discussions Results are presented in three parts: (i) flow field dynamics analysed via dominant frequencies, the impact of mixing models and extended/blended scalar time-scale closures on (ii) velocity statistics and (iii) scalar statistics (including scatter plots). For HM3, statistics were obtained over five flame instability cycles (corresponding to a time interval > 250 · Db /Ujet ) after the initial adjustment from the flamelet solution. 4.1. Flow field dynamics The flow dynamics for HM3 with U jet = 214 m/s (blow-off at 235 m/s) are strongly influenced, arguably dominated, by turbulence-chemistry interactions. The periodic behaviour obtained with the EMC model is illustrated in Figs. 2–4. Figure 2 shows that vortices are released from the outer shear layer and that the inner shear layer evolves to a jet-wake-like structure. The outer-shearlayer vortices move towards the jet centreline due to the low pressure in the recirculation zone. It can be seen that, at t = 5 ms (the time inside the parentheses represents the actual computational time), some of the vortices move back towards the recirculation zone. Such effects accumulate until the vortices are strong enough to interact with the main fuel jet wake at t = 5 ms. This interaction results in a vortex released from the jet wake convecting at nearly the same speed as the outer-shear layer vortices. The accumulated vortices at the end of the recirculation zone start to

56

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

Fig. 3. Instantaneous contours of temperature for the low-frequency instability period obtained using the EMC model. The time inside the parentheses represents the actual computational time as shown in Fig. 5.

fade and flow pattern reverts. The associated temperature variation in the same period can be seen in Fig. 3. The flame stabilised by the recirculation zone is pulled out towards the neck zone. A flame turbule [40] is separated and convected downstream simultaneously with the release of a vortex from the jet wake as shown in Fig. 2. In addition, it can be observed that local extinction occurs beyond the recirculation zone at t = 0 ms; however, the temperature downstream increases before the release of a flame turbule at t = 4 ms. The corresponding mixture fraction field shown in Fig. 4 suggests that values within the mixing layer are within the flammability limits (0.016 ≤ f ≤ 0.15) throughout the whole period. Hence, apart from flow dynamics, the local chemistry impacts the re-ignition behaviour. A flame turbule separates at the neck zone (around x = 90 mm) and the cycle is then repeated. The ability of EMC and EMST to reproduce the observed local extinction and reignition behaviour at the neck zone is discussed in Section 4.3.3. The instabilities in the outer shear layer and the frequency of shear-layer vortices/jet wake interactions were quantified using velocity and temperature traces obtained using a sampling rate of 5 kHz at different axial and radial locations. The three mixing layers: (i) the inner shear layer, (ii) the outer jet wake and (iii) the interaction between these two flow structures [21,27,67] are shown with sampling at three corresponding radial locations (r = 2 mm, 22 mm and 12 mm). The obtained with the EMC model are shown in Fig. 5. A high-frequency perturbation can be observed at r = 22 mm and at x = 9 mm in the bottom two panels of Fig. 5, while a low-frequency periodic motion can be detected at all ra-

Fig. 4. Instantaneous contours of mixture fraction for the low-frequency instability period obtained using the EMC model. The time inside the parentheses represents the actual computational time as shown in Fig. 5.

dial locations beyond x = 45 mm. Hence, the dominant frequencies were investigated using Fourier transforms of the velocity and temperature traces at r = 22 mm. As shown in the top panel of Fig. 6(a), a high-frequency (Fh ) of 1049 Hz and a low-frequency (Fl ) of 85 Hz are detected at x = 9 mm. The higher frequency is not as distinct due to complexities in the flow field. However, the value is broadly consistent with the frequency of vortex shedding in the outer shear layer, while the low frequency relates to the release of flame turbules. The bimodal feature in the top figure of Fig. 6(a) is lost along the axial direction. The magnitude of the high-frequency instability decreases significantly along the flow as the low-frequency feature becomes more prominent. The same trend is present in the temperature traces in Fig. 6(b) with multiple frequencies at x = 45 and 65 mm due to flame interactions with the recirculation zone. The HM3 flame is at 91% of the blow-off velocity and sufficiently close to relate local extinction and re-ignition to studies where detailed measurements have not been performed. The instability in the outer shear layer corresponding to the highfrequency observed in Fig. 6 has a shear layer Strouhal number Stsl = Fh θ /Uco f  0.010 based on the momentum thickness (θ ) and co-flow velocity. Consistent with the inlet boundary condition, the momentum thickness was estimated to be θ = 0.365 mm based on the boundary layer thickness of δ = 3.75 mm and the 1/7th power law (θ /δ = 7/72) [64]. The momentum thickness estimation [68,69] is not sensitive to small differences in velocity profiles [69] with 1/5 or 1/9 profiles resulting in differences below

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

a 180

b 2500

x=9mm x=45mm x=65mm x=90mm x=120mm

150 2000

120

1500

60

1000

30

500

0 90

0 2500

60

2000 T [K ]

U [m s]

90

30 0

1500 1000

−30

500

−60 60

0 2500

40

2000 1500

20 1000 0 −20 390

500 410

430 t [ms]

0 390

450

410

430 t [ms]

450

Fig. 5. Time histories of (a) mean axial velocity and (b) mean temperature for flame HM3 based on the EMC model. Top: r = 2 mm; Middle: r = 12 mm; Bottom: r = 22 mm.

a

b 20

x = 9 mm

0.4 0.3

15

0.2

10

0.1

5

0 3

x = 9 mm

0 50

x = 45 mm

x = 45 mm

Magnitude

Magnitude

40 2

1

30 20 10

0 3

x = 65 mm

0 80

x = 65 mm

60

2

40 1 20 0 1000

1500

10%. Shear-layer behaviours in ramjet combustors with backwardfacing steps have been investigated [70–72] using linear instability theory providing Stsl ࣃ 0.017. A Strouhal number for the local extinction and re-ignition cycle based on the bluff-body diameter and main jet velocity (Ster = Fl Db /U jet ), provides a value ࣃ 0.020. There are few studies involving such frequencies as the flame dynamics is related to acoustic-combustion interactions [73,74] and ignition time scales [75,76]. However, the local extinction and reignition process of a premixed methane flame near the lean blow out limit in a swirl-stabilised combustor has been analysed using proper orthogonal decomposition. A mode related to the interaction between the shear layer and flame showed Ster ࣃ 0.019 [77]. The dynamic behaviour of a bluff-body propane–air diffusion flame generating flame turbules has also been studied [40,41] and the behaviour attributed to interactions between the outer-shear-layer vortices and the fuel jet wake. Roquemore et al. [40] suggested that the fuel jet has a dominant influence on the frequency and obtained 0.040 ≤ Ster ≤ 0.095 depending on the fuel flow rate. The convective velocity of turbules in the present study can be found by analysing the instantaneous contours or calculating the time between two peaks in temperature traces as shown in Fig. 5(b). Consistent with Roquemore et al. [40], the mean convective velocity of the released flame turbule is close to the co-flow air velocity (here 40 m/s). The same methodology has been applied to analyse the flow field frequencies in the simulations based on the EMST based models and the flamelet solution. Results show that although there is no obvious local extinction in either flow field, a high shear-layer frequency can be detected. The values obtained from the EMST based models (e.g., 978 Hz for EMST) and the flamelet calculation (1097 Hz) are similar to that based on the EMC model (1049 Hz). The flamelet solution does not feature under-relaxation via Eq. (6). This suggests that the shear-layer instability is not related to turbulence-chemistry effects and that the applied under-relaxation does not substantially alter the observed vortex-shedding frequency. Potential inaccuracies are more likely related to the axisymmetric assumption and an increased instability amplitude obtained with second moment based URANS compared to LES [44]. The sensitivity to solution parameters was further investigated by varying the under-relaxation factor (NTA ) and time step. Cases with the EMC model were accordingly also computed (i) with an increase in NTA by 250% and (ii) a decrease in t by a factor of two. Results show deviations in the high-frequency of less than 5%. The low-frequency is insensitive to t in the current range, but decreases by 27% when NTA is increased by 250%. The low-frequency instability is related to vortex interactions and corresponds to local extinction and re-ignition events leading to the release of flame turbules. The latter affect the evolution of the mixing field and hence the turbulence-chemistry interactions. Furthermore, an improved understanding of interactions with the high-frequency shear layer instabilities is important in the context of the blow-off point. Detailed experimental information is currently not available for HM3. However, methods such as electrical connectivity have been used as a measure of intermittency with a frequency of 100 Hz detected at 50% intermittency for bluff-body LPG flames [11]. The current work suggests that an extension to include such experimental data for comprehensive data sets, such as the HM1-HM3 series, would be invaluable in understanding the stability limits of bluff-body flames. 4.2. Velocity statistics

0 500

57

500

1000

1500

Fig. 6. Dominant frequencies based on (a) velocity and (b) temperature traces at r = 22 mm for flame HM3 based on the EMC model.

Velocity statistics obtained from all mixing models/scalar timescale closures investigated (MC, EMC, BMC, EMST, EEMST and BEMST) are compared in the present section. It is noted that in the previous section the MC based models predict large periodic

58

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

x = 4 mm

1

x = 10 mm

x = 20 mm

0.6

0.01

0.2

−0.01

x = 10 mm

x = 20 mm

x = 30 mm

x = 40 mm

x = 50 mm

x = 60 mm

x = 90 mm

x = 120 mm

−0.03

−0.2 x = 30 mm

x = 40 mm

x = 50 mm

0.03

[- ]

[- ]

1

0.01 V U

0.6 U U

x = 4 mm 0.03

0.2

−0.01 −0.03

−0.2 x = 60 mm

1

x = 90 mm

x = 120 mm

0.03

0.6

0.01

0.2

−0.01 −0.03

−0.2 0

10

20 30 r [mm]

40

10

20 30 r [mm]

40

10

20 30 r [mm]

Fig. 7. Radial profiles of the normalised mean axial velocity for flame HM3-HM3E. (—– )

Symbols: (◦) measurements of flame HM3E [6]; lines: ) BMC; ( ) EMST; ( ) EEMST; ( ) BEMST. (

MC; (

0

40

10

10

20 30 r [mm]

40

10

20 30 r [mm]

40

) EMC;

x = 4 mm

x = 10 mm

x = 20 mm

x = 30 mm

x = 40 mm

x = 50 mm

x = 60 mm

x = 90 mm

x = 120 mm

0.1 0.05 0 [- ]

0.15

U

0.1 u

motions of flow field. Hence, as shown in Eqs. (11)–(12), the coherent fluctuation component (uc and vc ) was included along with the stochastic counterpart (us and vs ) to constitute Root Mean Square (RMS) data [44]. The former represents the fluctuation generated from periodic motions, while the latter can be obtained from a direct average of RMS at each time step. Their crosscorrelations have been found to be negligible at high Reynolds number [44]. A comparison of coherent and stochastic contributions on the axial and radial RMS velocities based on the EMC model suggests that the coherent component of the axial RMS velocity cannot be neglected after x = 40 mm where periodic motions become significant. In particular, at x = 120 mm, the coherent component is around 20% of the stochastic counterpart. The coherent contribution on radial RMS velocity is, however, modest with only a slight impact around r = 20 mm where outer shear layer vortices are present. To be consistent, the coherent component has been included in the RMS data for all mixing models/scalar timescale closures though the effects are found to be negligible for the EMST based models.

0.05 0 0.15 0.1 0.05

(11) 0 0

v = vc + vs

40

Fig. 8. Radial profiles of the normalised mean radial velocity for flame HM3-HM3E. Symbols and lines as Fig. 7.

0.15

u = uc + us

20 30 r [mm]

(12)

Figures 7–10 show radial profiles of velocity statistics at different axial locations obtained from all mixing models/scalar timescale closures investigated (MC, EMC, BMC, EMST, EEMST and BEMST). Predictions of flame HM3 are compared to measurements of flame HM3E [11]. Flame HM3E has a lower bulk jet velocity ( ࣃ 10% reduction) with a similar jet to co-flow velocity ratio (5.6 compared to 5.4) as shown in Table 1. Accordingly, flow field parameters have been scaled with the respective bulk jet velocities. Mean axial velocities, as shown in Fig. 7, are in excellent agreement with measurements except slight discrepancies beyond the neck zone (x = 90 mm). The different mixing models/scalar time-scale closures cause modest differences due to the different levels of lo-

10

20 30 r [mm]

40

10

20 30 r [mm]

40

10

20 30 r [mm]

40

Fig. 9. Radial profiles of the normalised axial RMS velocity for flame HM3-HM3E. Symbols and lines as Fig. 7.

cal extinction affecting the mean density. The MC based models slightly underestimate the mean axial velocities close to the centreline at x = 90 mm and x = 120 mm. Mean radial velocities are also in reasonable agreement with the experiment. As shown in Fig. 8, differences between the mixing models are mainly present within the recirculation zone and the MC based models capture the experimental features up to x = 40 mm before a slight misalignment for the second peak at x = 50 mm and 10 ≤ r (mm) ≤ 20. As shown in Fig. 2, the interactions of the outer shear layer and

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

x = 4 mm

x = 10 mm

x = 20 mm

0.6

0.12

x = 13 mm

x = 30 mm

x = 45 mm

x = 65 mm

x = 90 mm

x = 120 mm

0.4

0.08

0.2

0.04

0

0 x = 30 mm

x = 40 mm

0.6

x = 50 mm ]

0.12 0.4

U

[

[- ]

59

0.08

v

0.2

0.04 0

0

0.6

x = 60 mm

x = 90 mm

x = 120 mm

0.12

0.4

0.08

0.2

0.04

0 0

0 0

10

20 30 r [mm]

40

10

20 30 r [mm]

40

10

20 30 r [mm]

40

5

10

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

30

Fig. 11. The segregation factor (see Eq. (7)) obtained from the BEMST model.

Fig. 10. Radial profiles of the normalised radial RMS velocity for flame HM3-HM3E. Symbols and lines as Fig. 7.

the fuel-jet are intense in this region. The prediction of the decay of the second peak depends on the mixing rate and the MC based models reproduce this feature adequately at x = 60 mm. The predicted mean radial velocity peak at x = 90 mm is closer to the centreline than measured and the measured peak at x = 120 mm is higher than the predictions. Such discrepancies can be found in the simulation of HM1 using CMC/LES [27]. Hence, a slight asymmetry in the experimental data may play a role. Similar trends can be also observed in the Reynolds stresses components shown in Figs. 9 and 10, though with a clear tendency to over–predict fluctuations close to the centreline for 40 ≤ x (mm) ≤ 50. The RMS velocities predicted by BMC tend to be lower than other closures at x = 90 mm and x = 120 mm. The velocity statistics obtained with the different models are similar with variations in the flow field caused by different levels of local extinction impacting the mean density as discussed in detail in the Supplementary Material. 4.3. The impact of mixing models on scalar statistics Results for scalar statistics are presented below in physical and conditional space. The impact of different treatments of molecular-mixing is highlighted. The blended scalar time-scale closure [39] was proposed for premixed turbulent flames. However, the conceptual similarity with the extended scalar time-scale closure suggests that the approach may have greater applicability in the context of (increasingly) combustion regime independent methods. Figure 11 shows the evolution of the segregation factor (η) obtained from the BEMST model at different axial locations. The flamelet behaviour comes into effect away from the bluffbody surface as η increases after initially featuring a double-peak at the location of the two shear layers. The two peaks then begin to merge at x = 13 mm and the single peak aligns around 10 ≤ r (mm) ≤ 15 downstream. 4.3.1. Radial profiles The impact of the different scalar dissipation rate closures is evaluated in two groups based on the MC and EMST mixing mod-

els. The mixture fraction [78] is a passive scalar and predominantly affected by mixing processes driven by the flow field. Radial profiles of the predicted mean and fluctuating mixture fractions by the MC based models at different axial locations are compared to the measurement of flame HM3 in Fig. 12. In addition to the current simulations, previous results based on JPDF/LES [28] are also included for reference. As shown in Fig. 12, all closures provide initially good predictions of the mean and fluctuating mixture fractions with some discrepancies appearing downstream. The impact of the extended scalar time-scale closure is modest, while the blended closure tends to overestimate mixing far downstream. This relates to the magnitude of η (see Fig. 11) and an overestimation of the contribution from the reactive part. The MC based models tend to under-predict the fluctuations at x = 90 mm and x = 120 mm with BMC aggravating the issue. The latter model leads to global extinction of the flame beyond the recirculation zone. Similarly, EMST based models provide accurate predictions of mean mixture fractions as shown in Fig. 13. The EEMST variant improves mixture fraction fluctuations close to the bluff-body, but slightly underestimates the value on the fuel lean side far downstream. Similar to BMC, BEMST tends to underestimate the fluctuations beyond x = 65 mm. Overall, apart from the blended closures, the current predictions of mixture fraction statistics are similar to the JPDF/LES simulations. Based on the above observations, EMC, EMST and EEMST predictions of temperature and species mass fractions are compared below. EMC has been shown to outperform MC for flame HM2 [26]. The performance of EMST and EEMST is also compared to further investigate the impact of extended scalar time-scale closure. Results based on the other closure combinations are available in the Supplementary Material. Radial profiles of the mean and fluctuating temperatures are compared with experimental data at different axial locations in Fig. 14. The mean temperature obtained from a flamelet solution is also included as a reference for both cases. As shown in Fig. 14, the EMC model provides good predictions of mean temperature close to the bluff-body, but still underpredicts the mean temperature beyond the recirculation zone due to less resistance to local extinction. The EMC model captures the double-peak feature of the fluctuating temperature at x = 13 mm,

60

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

1

x = 13 mm

1

x = 30 mm

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0 1

x = 45 mm

0 1

x = 65 mm

0.8

x = 30 mm

x = 45 mm

x = 65 mm

x = 90 mm

x = 120 mm

]

]

0.8

0.6

0.6 f [

f [

x = 13 mm

0.4 0.2

0.4 0.2

0 1

x = 90 mm

0 1

x = 120 mm

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0 0

5

10

0.2

15 20 r [mm]

25

30

5

10

x = 13 mm

15 20 r [mm]

25

0

30

0

x = 30 mm

5

10

0.2

0.15

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

x = 13 mm

x = 30 mm

x = 45 mm

x = 65 mm

x = 90 mm

x = 120 mm

30

0.15

0.1

0.1 0.05

0.05 x = 45 mm

0 0.2

x = 65 mm

0.15

[

]

]

0 0.2

0.15

f

f

[

0.1 0.05

0.1 0.05

0 0.2

x = 90 mm

0 0.2

x = 120 mm

0.15

0.15 0.1

0.1

0.05

0.05

0 0

5

10

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

30

Fig. 12. Radial profiles of the (top) mean and (bottom) fluctuating mixture fractions for flame HM3 based on MC with different scalar dissipation rate closures. Symbols: (◦) measurements of flame HM3 [16]; lines: (—–) MC; (- - -) EMC; (- · -) BMC; ( · ·  · · ) JPDF/LES by Popov and Pope [28].

which indicates the presence of outer and inner shear layers, and provides reasonable agreement at other axial locations. The discrepancies at x = 90 mm and x = 120 mm are associated with the (modest) under-prediction of mixture fraction statistics shown in Fig. 12. The periodic features predicted by the MC based models may be excessive and correspondingly suppress mean values. The mean temperatures predicted by EMST and EEMST are closer to the flamelet solution. The EMST model slightly underestimates the mean temperature at the outer shear layer close to the bluff-body, but overestimates the fluctuating temperature in the same region. Both the mean and fluctuating temperatures are over-predicted by EMST beyond x = 90 mm. The discrepancies beyond x = 90 mm are finite-rate chemistry related and the EMST model is arguably too resistant to local extinction. The EEMST model captures the

0 0

5

10

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

30

Fig. 13. Radial profiles of the (top) mean and (bottom) fluctuating mixture fractions for flame HM3 based on EMST with different scalar dissipation rate closures. Symbols: (◦) measurements of flame HM3 [16]; lines: (—–) EMST (- - -) EEMST; (- · -) BEMST; ( · ·  · · ) JPDF/LES by Popov and Pope [28].

two-peak feature of fluctuating temperature close to the bluffbody and notably improves the prediction of peak locations beyond x = 65 mm. Results based on EMC, EMST and EEMST are shown for major and minor species in Figs. 15–17. Radial profiles of the mean and fluctuating CO mass fractions, shown in Fig. 15, share similar trends with temperature. The mean CO mass fractions predicted by the EMC model agree well with the experiment close to the bluffbody, but become underestimated beyond x = 45 mm. The underprediction is associated with the reduced temperature at the corresponding locations as shown in Fig. 14, which is partially due to the contribution of the periodic motion to the statistics. The EMST model shows generally good agreement with experimental data for

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

2500

x = 13 mm

0.06

x = 30 mm

61

x = 13 mm

x = 30 mm

x = 45 mm

x = 65 mm

x = 90 mm

x = 120 mm

2000

0.04

1500 1000

0.02

500 0 2500

x = 45 mm

0 0.06

x = 65 mm

[

]

1500

0.04

1000

Y

T [K ]

2000

0.02

500 0 2500

x = 90 mm

0 0.06

x = 120 mm

2000

0.04

1500 1000

0.02

500 0 0

5

10

800

15 20 r [mm]

25

30

5

10

x = 13 mm

15 20 r [mm]

25

0

30

0

x = 30 mm

5

10

0.03

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

x = 13 mm

x = 30 mm

x = 45 mm

x = 65 mm

x = 90 mm

x = 120 mm

30

600

0.02 400

0.01

200

0 x = 45 mm

x = 65 mm

0.03 [

]

600 400

0.02

Y

T

[K ]

0 800

0.01

200 0 800

0 x = 90 mm

x = 120 mm

0.03

600

0.02

400

0.01

200 0 0

5

10

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

30

0 0

5

10

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

30

Fig. 14. Radial profiles of the (top) mean and (bottom) fluctuating temperatures for flame HM3. Symbols: (◦) measurements of flame HM3 [16]; lines: (—–) EMC; (- - -) EMST; (- · -) EEMST; ( · · · ) flamelet solution.

Fig. 15. Radial profiles of the (top) mean and (bottom) fluctuating CO mass fractions for flame HM3. Symbols: (◦) measurements of flame HM3 [16]; lines: (—–) EMC; (- - -) EMST; (- · -) EEMST.

both the mean and fluctuating CO mass fractions up to x = 65 mm, but tends to overestimate both beyond this point due to a lack of local extinction. The EEMST model improves the prediction of the fluctuating CO mass fraction on the fuel rich side close to the bluffbody and the peak locations far downstream. Predictions for the main products CO2 and H2 O share the same trend as the temperature and are shown in the Supplementary Material. Radial profiles of minor species, such as OH and NO, are presented in Figs. 16–17. All models tend to under-predict the mean OH mass fraction close to the outer shear layer at x = 30 mm. This is consistent with the under-estimation of the mean temperature in the same region as illustrated in Fig. 14. As shown in Fig. 16, the EMC model shows good agreement for the fluctuating OH mass

fraction at all axial locations, but underestimates the mean OH mass fraction at x = 45 mm and x = 65 mm due to the underprediction of temperature (excessive local extinction). The EMST model overestimates the mean and fluctuating OH mass fractions for x ≥ 90 mm due to a lack of local extinction. The EEMST model increases mean OH mass fractions, as compared to EMST, and provides an improved profile shape. Predictions of NO present a long-term challenge due to the strong finite-rate chemistry effects. The NO chemistry, which was successfully applied to flame HM2 [26], is here used for all simulations of HM3. It can be seen in Fig. 17 that the EMC model provides reasonable agreement for both the mean and fluctuating NO mass fractions. The NO levels predicted by the EMST/EEMST mod-

62

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

x = 13 mm

0.04

x = 30 mm

0.015

x = 13 mm

x = 30 mm

x = 45 mm

x = 65 mm

x = 90 mm

x = 120 mm

0.03 0.01

0.02 0.005

x = 65 mm

]

0.03

0

x = 45 mm

100 [

]

0 0.04

10[

0.01

0.015 0.01

0.01

Y

Y

0.02

0 0.04

x = 90 mm

0.005 0

x = 120 mm

0.015

0.03 0.01

0.02 0.01

0.005

0 0

5

10

0.03

15 20 r [mm]

25

30

5

10

x = 13 mm

15 20 r [mm]

25

0

30

0

x = 30 mm

5

10

0.06

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

x = 13 mm

x = 30 mm

x = 45 mm

x = 65 mm

x = 90 mm

x = 120 mm

30

0.05

0.02

0.04 0.03

0.01

0.01

x = 65 mm

]

x = 45 mm

1000 [

]

0.02

0 0.06 0.05 0.04 0.03

0.01

0.02 Y

Y

0 0.03

10[

0.02

0 0.03

x = 90 mm

0.01 0 0.06

x = 120 mm

0.05

0.02

0.04 0.03

0.01

0.02 0.01

0 0

5

10

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

30

Fig. 16. Radial profiles of the (top) mean and (bottom) fluctuating OH mass fractions for flame HM3. Symbols and lines as in Fig. 15.

els are substantially overestimated due to a lack of local extinction. However, the EEMST model improves the predictions of fluctuations close to the bluff-body and the peak locations far downstream. The wavy shape of the fluctuating NO mass fractions is consistent with the prediction of the fluctuating temperatures. It is relevant that recent JPDF/LES simulations [28] based on the IEM model also over-predict the means while using a GRI-based chemistry model. 4.3.2. Conditional statistics The performance of EMC and EMST for HM1 and HM2, which are at 50% and 75% of the blow-off velocity, as well as for HM3 is relevant due to the increasing levels of local extinction. Results using EMC for HM2 were presented in a previous study [26] and have been recomputed for consistency. Conditional mean temperatures

0 0

5

10

15 20 r [mm]

25

30

5

10

15 20 r [mm]

25

30

Fig. 17. Radial profiles of the (top) mean and (bottom) fluctuating NO mass fractions for flame HM3. Symbols and lines as in Fig. 15.

obtained with a bin width of  f = 0.005 are presented for flames HM1 and HM2 alongside previous CMC/LES or JPDF/LES simulations in Fig. 18. The performance of the two mixing models is similar for HM1 and provides good agreement with the experiment. The level of agreement remains for HM2, except that the two mixing models envelop the experimental data far downstream. Overall, the current JPDF/URANS results show similar agreement with experimental data as CMC/LES or JPDF/LES simulations. Comparisons of the EMC, EMST and EEMST models for conditional mean temperatures using the same bin width are shown for HM3 in Fig. 19. JPDF/LES results [28] at available locations are included for comparison. The laminar results of opposed jets (CH4 :H2 ; 1:1 against air) are included for reference. All models provide reasonable agreement with experimental data at x =

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

2500

x = 13 mm

x = 30 mm

2000

1500

1500

1000

1000

500

500

0 2500

x = 45 mm

x = 65 mm

2000

0 2500

x = 13 mm

x = 30 mm

x = 45 mm

x = 65 mm

x = 90 mm

x = 120 mm

2000 T [K ]

T [K ]

2500

2000

63

1500

1500

1000

1000

500

500

0 2500

x = 90 mm

x = 120 mm

0 2500

2000

2000

1500

1500

1000

1000

500

500

0 0

2500

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 [- ] [- ] x = 13 mm

x = 30 mm

2000

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 [- ] [- ]

Fig. 19. Conditional mean temperature for flame HM3. Symbols: (◦) measurements; lines: (—–) EMC; (- - -) EMST; (- · -) EEMST; ( · · · ) flamelet line.

1500 1000

0.08

500

0.06

0 2500

x = 13 mm

x = 30 mm

x = 45 mm

x = 65 mm

x = 90 mm

x = 120 mm

0.04 x = 45 mm

x = 65 mm

0.02

T [K ]

2000

0 ]

0.08

10 [

1500

0.06

500 0 2500

x = 90 mm

x = 120 mm

2000

Y OH

1000

0.04 0.02

1500

0

1000

0.08

500

0.06

0 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 [- ] [- ]

0.04 0.02

Fig. 18. Conditional mean temperature for flames (top) HM1 and (bottom) HM2. Symbols: (◦) measurements; lines: (—–) EMC; (- - -) EMST; (- · -) (top) CMC/LES by Kim and Pitsch [27] or (bottom) JPDF/LES by Popov and Pope [28].

13 mm. Discrepancies start at x = 30 mm where EMC and EMST again bracket the experimental mean temperature. A substantially larger underestimation by the EMC model can be seen at x = 45 mm, partially due to statistics averaged to include periodic motions (e.g., flame turbules). The EEMST model increases the conditional mean temperatures. The JPDF/LES method [28], which can provide a better prediction of scalar dissipation rate, is more accurate in the recirculation zone. The different behaviours of the two types of mixing models persist downstream (conditional experimental data is not available at x = 65 mm). Consistent trends can be found for conditional OH mass fractions in Fig. 20. The locations of the maxima are accurately reproduced. EMC and (E)EMST models envelop the experimental data. Figures 19–20 confirm that

0 0

0.1 [- ]

0.2

0.1 [- ]

0.2

Fig. 20. Conditional mean OH mass fraction for flame HM3 (dashed vertical line represents the stoichiometric value). Symbols: (◦) measurements; lines: (—–) EMC; (- - -) EMST; (- · -) EEMST; ( · · · ) flamelet line.

EMC tends to overestimate local extinction, while (E)EMST shows the opposite trend. The performance of EMC and EMST was further investigated by comparing burning indices (BI) [52,55] based on the conditional mean temperature defined by Eq. (13),

BI =

T | flean ≤ f ≤ frich − Tu Tad − Tu

(13)

where T|flean ≤ f ≤ frich represents the mean temperature conditioned on the flammable mixture fraction range with Tu (= 298 K) the reactant temperature and Tad (= 2272 K) the adiabatic flame

64

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

2400

1

T [K ]

0.4 x = 13 mm

x = 30 mm

1500 1200

600 300

0.8 BI[

]

0

0.1

0.6 0.4

0.2 [- ]

0.3

0.1

0.2 [- ]

0.3

0.1

0.2 [- ]

0.3

Fig. 22. Scatter plots of temperature for flame HM3 at x = 90 mm.

x = 45 mm

x = 65 mm 0.01

0 1

Exp.

EMC

EMST

0.008 Y OH [- ]

0.8 0.6 0.4 0.2

EMST

900

0 1

0.2

EMC

1800

0.6

0.2

Exp.

2100

0.8

x = 90 mm

x = 120 mm

0.006 0.004 0.002

0 100 120 140 160 180 200 220 240 120 140 160 180 200 220 240 U [m s] U [m s]

0 0

Fig. 21. Burning indices of conditional mean temperature against main jet velocities. Symbols: (◦) measurements; lines: (—–) EMC; (- - -) EMST; () MC; () EEMST; () BEMST;.

0.1

0.2 [- ]

0.3

0.1

0.2 [- ]

0.3

0.1

0.2 [- ]

0.3

Fig. 23. Scatter plots of OH mass fraction for flame HM3 at x = 90 mm.

4.3.3. Local extinction and re-ignition Burning indices shown in Fig. 21 suggest that the EMC and (E)EMST models exhibit different sensitivities to local extinction and tend to envelop the experimental data. In addition, as discussed in Section 4.1, local extinction and re-ignition occurs at the neck zone (around x = 90 mm). Hence, scatter plots based on EMC and EMST are also compared to experimental data in mixture fraction space. It can be seen from Fig. 22 that the EMC model tends to overestimate the percentage of extinguished samples. By contrast, the EMST model significantly underestimates the amount of local extinction. Similar trends can be found in OH and NO mass fractions as illustrated in Figs. 23 and 24. The shape of predicted OH mass fractions is the same for EMC and EMST, but the EMC model shows a higher percentage of extinguished samples. The current results suggest that local extinction plays a crucial role in determining NO levels for HM3. The experimental data for OH mass fractions may be affected by interference as samples appear outside flammability limits.

Exp.

0.01

EMC

EMST

0.008 0.006

Y NO

temperature. The latter was obtained via an opposed-flow laminar flame calculation with a strain rate a = 100 s−1 . The mixture fraction range (0.016 ≤ f ≤ 0.15) approximately corresponds to the flammability limits. A burning index of unity corresponds to a full burning and zero implies global extinction. Burning indices based on EMC and EMST are plotted against experimental data for HM1– HM3 in Fig. 21 with MC, EEMST and BEMST included for flame HM3. The experimental data is lower than the predictions close to the bluff-body. This may be attributed to heat losses to the bluff-body surface. The conditional experimental data for HM3 at x = 65 mm is not available. The experimental data shows a clear trend towards blow-off beyond x ≥ 45 mm with increasing jet velocity. Consistent with radial profiles, EMC and EEMST/BEMST affected by enhanced mixing tend to increase the means compared to the standard closures. Overall, EMC and (E)EMST tend to bracket the experimental data and further improvements can likely be obtained through optimised modelling constants.

1 0 0[ - ]

0.012

0.004 0.002 0 0

0.1

0.2 [- ]

0.3

0.1

0.2 [- ]

0.3

0.1

0.2 [- ]

0.3

Fig. 24. Scatter plots of NO mass fraction for flame HM3 at x = 90 mm.

5. Conclusions The emphasis of the current study is to clarify the impact of different mixing and scalar dissipation rate closures on velocity and scalar statistics for a comprehensively characterised CH4 :H2 (1:1) flames series (HM1,HM2,HM3) as global extinction is approached and multiple combustion modes co-exist [6,11–18]. Wellresolved computations featuring a compressible hybrid transported JPDF and second moment based URANS approach [3,26] that resolves the unsteady bulk motion of the flow were performed and the flow dynamics was analysed using Fourier transforms of velocity and temperature traces. Results show that all closures capture the mean and RMS of the velocity field reasonably well. The outershear-layer vortex shedding frequency around 1 kHz was found to be comparatively insensitive to the applied temporal resolution and the thermochemical closure. However, the flame dynamics differ close to global extinction with EMST based models less prone to incipient blow-off. By contrast, EMC showed a prominent periodic pattern of local extinction and re-ignition resulting in the release of flame turbules for HM3 (at 91% of the blow-off velocity) with a frequency of ∼ 85 Hz. The instability is induced by the interaction between outer-shear-layer vortices and the fuel jet wake, resulting in flame turbules convected downstream at a speed close to the co-flow air velocity. Similar behaviour has been observed in

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

related experimental studies e.g., [40,41]. The pattern is comparatively insensitive to temporal resolution parameters with a 250% increase in under-relaxation for the density field resulting a 27% change in the dominant frequency. The impact of mixing models/scalar time-scale closures on scalar statistics was further investigated in physical and mixture fraction space. Radial profiles suggest that EMC provides reasonable agreement close to the bluffbody surface, but tends to underestimate the downstream temperature and major species mass fractions for HM3 partly as a result of excessive periodic flame-turbule releasing cycles. By contrast, EMST underestimates local extinction beyond the neck zone and EEMST shows substantially increased mixing within flamelets leading to increased mean values and suppressed fluctuations close to the bluff-body and improved predictions of scalar peak locations beyond the neck zone. The MC based models are more sensitive to the scalar time-scale closure with BMC leading to global extinction. Overall, the blended approach tends to over-estimate the contributions from reactive scalars beyond x = 90 mm. Excellent agreement for conditional temperatures with experimental data and JPDF/LES or CMC/LES simulations is obtained for HM1 and HM2. Predictions of HM3 for the two types of mixing models (EMC and EMST/EEMST) essentially envelop the experimental data. It can be noted that some species more directly linked to the (faster) chemistry of the H2 fuel component appear closer to the EMST based solutions, while others such as NO show comparatively large deviations. Based on a comparison of scatter plots of temperature and NO mass fractions, emissions of NO are controlled by local extinction events. Overall, the results suggests that the EMC/EEMST approaches can form the basis for combustion regime independent closures for molecular mixing. EEMST, in particular, provides good predictions apart from reduced levels of local extinction. The blended closures can likely be enhanced with improved blending coefficient (η) models. It is also encouraging that the very different underlying molecular mixing closures essentially bracket the experimental data for the current complex flame operating close to global extinction. It is further concluded that experimental investigations providing detailed information on flame instabilities and flow dynamics are required to further advance the associated modelling.

65

Reynolds stresses

   u  u  u  ∂ ρ¯ u ∂ ρ¯ ul u ∂ u ∂ k i j i j i j   + = C ρ¯ uk ul ∂t ∂ xl ∂ xk s   ∂ xl   ∂ uj  ∂ ui      − ρ¯ ui ul + u j ul ∂ xl ∂ xl  u + ρ¯ G u    + ρ¯ Gik u ¯ C0 δi j   (A.4) jk k ui + ρ k j

where Cs = 0.22 and C0 = 2.1, and the form of Gij is a second– rank tensor derived on the basis of the Generalised Langevin Model (GLM) [2],

  ∂ u Gi j = (α1 δi j + α2 bi j ) + Hi jkl k  ∂ xl k Hi jkl =

(A.5)

β1 δi j δkl + β2 δik δ jl + β3 δil δ jk + γ1 δi j bkl + γ2 δik b jl +γ3 δil b jk + γ4 bi j δkl + γ5 bik δ jl + γ6 bil δ jk

(A.6)

and

1 3 1 α1 = − + C0 + b2ii (α2 + α3 ) + b3ii α3 2 4 3 1 ∗ +(β2 + β3 + γ )I1 + γ ∗ I2 , 3

α2 = 3.78, α3 = 0, β1 = −0.2, β2 = 0.8, β3 = −0.2, γ2 = 1.04, γ3 = −0.34, γ5 = 1.99, γ6 = −0.76, γ ∗ ≡ γ2 + γ3 + γ5 + γ6 , γ1 = γ4 = 0.0,  u   ∂ u u δ 1 k ∂u bi j =

i

j

 u u k k

−

ij

3

, Si j ≡

2 

i

∂xj

+

j

∂ xi

, I1 ≡ bi j Si j ,

I2 ≡ b2i j Si j , b2ii = bi j b ji

(A.7)

Dissipation

   ∂ ρ¯   ∂ ρ¯ ul   ∂ k ∂    + = C ρ¯ uk ul ∂t ∂ xl ∂ xk se   ∂ xl      ∂ uk  − C1 ρ¯ u u − C2 ρ¯   k l ∂x   l k k

(A.8)

where CSe = 0.18, C1 = 1.44, C2 = 1.80. Joint scalar pdf transport equation

Acknowledgments The authors wish to gratefully acknowledge discussions with H. Zhou (Tsinghua University) on the implementation of EMST as well as K. Gkagkas and T.S. Kuan. Lu Tian would like to acknowledge the financial support offered by the Imperial College President’s PhD Scholarship for exceptional students.

   ∂ Fφ ∂ ui Fφ ∂ Sα ( )Fφ ∂  1 ∂ Jiα   + + = F  φ ∂t ∂ xi ∂ α ∂ α ρ ∂ x i    ∂ ui | Fφ − ∂ xi

(A.9)

where Appendix A. Governing equations for the hybrid JPDF/URANS approach

∂ ρ¯ ∂ ρ¯ ul + =0 ∂t ∂ xl

τil = μ

(A.10)

(A.1)

Supplementary material

(A.2)

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.combustflame.2019.04. 039.

Momentum

where

  ∂ ui | Fφ ∂  μT ∂ (Fφ /ρ¯ ) = ∂ xi ∂ x i σT ∂ xi

and σ T is the turbulent Schmidt or Prandtl number, assumed to 0.7 in the present study.

Continuity

∂ ρ¯ ui + ∂t

−

∂ ρ¯ ul ui ∂ p¯ =− + ∂ xl ∂ xi

∂  u ) (τ − ρ¯ u i l ∂ xl il

References

∂ ui ∂ ul + ∂ xl ∂ xi

 −

2 ∂ uk μ δ 3 ∂ xk il

and where μ is the dynamic laminar viscosity.

(A.3)

[1] R.S. Barlow, S. Meares, G. Magnotti, H. Cutcher, A.R. Masri, Local extinction and near-field structure in piloted turbulent CH4 /air jet flames with inhomogeneous inlets, Combust. Flame 162 (10) (2015) 3516–3540. [2] S.B. Pope, PDF methods for turbulent reactive flows, Prog. Energy Combust. Sci. 11 (2) (1985) 119–192.

66

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67

[3] T.S. Kuan, R.P. Lindstedt, Transported probability density function modeling of a bluff body stabilized turbulent flame, Proc. Combust. Inst. 30 (1) (2005) 767–774. [4] R.P. Lindstedt, E.M. Váos, Transported PDF modeling of high-Reynolds-number premixed turbulent flames, Combust. Flame 145 (3) (2006) 495–511. [5] L. Tian, R.P. Lindstedt, Evaluation of reaction progress variable - mixture fraction statistics in partially premixed flames, Proc. Combust. Inst. 37 (2) (2019) 2241–2248. [6] A.R. Masri, Databases for bluff-body flows and flames (2006). [online] [7] H. Pitsch, Large-eddy simulation of turbulent combustion, Annu. Rev. Fluid Mech. 38 (2006) 453–482. [8] J. O’Brien, C.A.Z. Towery, P.E. Hamlington, M. Ihme, A.Y. Poludnenko, J. Urzay, The cross-scale physical-space transfer of kinetic energy in turbulent premixed flames, Proc. Combust. Inst. 36 (2016) 1967–1975. [9] L. Tian, R.P. Lindstedt, The impact of dilatation, scrambling, and pressure transport in turbulent premixed flames, Combust. Theor. Model. 21 (6) (2017) 1114–1147. [10] C.Y. Lee, L.K.B. Li, M.P. Juniper, R.S. Cant, Nonlinear hydrodynamic and thermoacoustic oscillations of a bluff-body stabilised turbulent premixed flame, Combust. Theor. and Model. 20 (1) (2016) 131–153. [11] A.R. Masri, R.W. Bilger, Turbulent diffusion flames of hydrocarbon fuels stabilised on a bluff body, Symp. (Int.) Combust. 20 (1) (1985) 319–326. [12] A.R. Masri, R.W. Dibble, R.S. Barlow, Raman–Rayleigh measurements in bluff– body stabilised flames of hydrocarbon fuels, Symp. (Int.) Combust. 24 (1992) 317–324. [13] A.R. Masri, B.B. Dally, R.S. Barlow, C.D. Carter, The structure of the recirculation zone of a bluff-body combustor, Symp. (Int.) Combust. 25 (1994) 1301–1308. [14] B.B. Dally, A.R. Masri, R.S. Barlow, G.J. Fiechtner, D.F. Fletcher, Measurements of NO in turbulent non-premixed flames stabilized on a bluff body, Symp. (Int.) Combust. 26 (1996) 2191–2197. [15] A.R. Masri, R.W. Dibble, R.S. Barlow, The structure of turbulent nonpremixed flames revealed by Raman–Rayleigh-LIF measurements, Prog. Energy Combust. Sci. 22 (4) (1996) 307–362. [16] B.B. Dally, A.R. Masri, R.S. Barlow, G.J. Fiechtner, Instantaneous and mean compositional structure of bluff-body stabilized nonpremixed flames, Combust. Flame 114 (1) (1998) 119–148. [17] A.R. Masri, J.B. Kelman, B.B. Dally, The instantaneous spatial structure of the recirculation zone in bluff-body stabilized flames, Symp. (Int.) Combust. 27 (1998) 1031–1038. [18] B.B. Dally, A.N. Karpetis, R.S. Barlow, Structure of turbulent non-premixed jet flames in a diluted hot coflow, Proc. Combust. Inst. 29 (2002) 1147–1154. [19] C. Fotache, T. Kreutz, C. Law, Ignition of hydrogen-enriched methane by heated air, Combust. Flame 110 (4) (1997) 429–440. [20] A.M. Briones, S.K. Aggarwal, V.R. Katta, Effects of H2 enrichment on the propagation characteristics of CH4 air triple flames, combust, Flame 153 (3) (2008) 367–383. [21] B.B. Dally, D.F. Fletcher, A.R. Masri, Flow and mixing fields of turbulent bluff– body jets and flames, Combust. Theor. Model. 2 (2) (1998) 193–219. [22] G.X. Li, B. Naud, D. Roekaerts, Numerical investigation of a bluff-body stabilised nonpremixed flame with differential Reynolds-stress models, Flow Turbul. Combust. 70 (1-4) (2003) 211–240. [23] A. Kempf, R.P. Lindstedt, J. Janicka, Large-eddy simulation of a bluff-body stabilized nonpremixed flame, Combust. Flame 144 (1-2) (2006) 170–189. [24] K. Liu, S.B. Pope, D.A. Caughey, Calculations of bluff-body stabilized flames using a joint probability density function model with detailed chemistry, Combust. Flame 141 (1) (2005) 89–117. [25] B. Merci, D. Roekaerts, B. Naud, S.B. Pope, Comparative study of micromixing models in transported scalar PDF simulations of turbulent nonpremixed bluff body flames, Combust. Flame 146 (1) (2006) 109–130. [26] K. Gkagkas, R.P. Lindstedt, T.S. Kuan, Transported PDF modelling of a high velocity bluff-body stabilised flame (HM2) using detailed chemistry, Flow Turbul. Combust. 82 (4) (2009) 493–509. [27] S.H. Kim, H. Pitsch, Mixing characteristics and structure of a turbulent jet diffusion flame stabilized on a bluff-body, Phys. Fluids 18 (7) (2006) 075103. [28] P.P. Popov, S.B. Pope, Large eddy simulation/probability density function simulations of bluff body stabilized flames, Combust. Flame 161 (12) (2014) 3100–3133. [29] H. Zhang, A. Giusti, E. Mastorakos, LES/CMC modelling of ignition and flame propagation in a non-premixed methane jet, Proc. Combust. Inst. 37 (2019) 2125–2132. [30] A.Y. Klimenko, S. Pope, The modeling of turbulent reactive flows based on multiple mapping conditioning, Phys. Fluids 15 (7) (2003) 1907–1925. [31] S.K. Ghai, S. De, A. Kronenburg, Numerical simulations of turbulent lifted jet diffusion flames in a vitiated coflow using the stochastic multiple mapping conditioning approach, Proc. Combust. Inst. 37 (2) (2019) 2199– 2206. [32] A.P. Wandel, R.P. Lindstedt, A mixture-fraction-based hybrid binomial Langevin-multiple mapping conditioning model, Proc. Combust. Inst. 37 (2019) 2151–2158. [33] B. Sundaram, A. Klimenko, A PDF approach to thin premixed flamelets using multiple mapping conditioning, Proc. Combust. Inst. 36 (2) (2017) 1937–1945. [34] C. Dopazo, Relaxation of initial probability density functions in the turbulent convection of scalar fields, Phys. Fluids 22 (1) (1979) 20–30. [35] J. Janicka, W. Kolbe, W. Kollmann, Closure of the transport equation for the probability density funcfion of turbulent scalar fields, J. Non-Equilib. Thermodyn. 4 (1) (1979) 47–66.

[36] S. Subramaniam, S.B. Pope, A mixing model for turbulent reactive flows based on Euclidean minimum spanning trees, Combust. Flame 115 (4) (1998) 487–514. [37] Z. Ren, S. Subramaniam, S.B. Pope, Implementation of the EMST mixing model [online] (2002). [38] R.P. Lindstedt, V.D. Milosavljevic, M. Persson, Turbulent burning velocity predictions using transported PDF methods, Proc. Combust. Inst. 33 (1) (2011) 1277–1284. [39] M. Kuron, Z. Ren, E.R. Hawkes, H. Zhou, H. Kolla, J.H. Chen, T. Lu, A mixing timescale model for TPDF simulations of turbulent premixed flames, Combust. Flame 177 (2017) 171–183. [40] W.M. Roquemore, R.L. Britton, S. Sandhu, Dynamic behavior of a bluff-body diffusion flame, AIAA J. 21 (10) (1983) 1410–1417. [41] R.D. Brum, L.M. Ikioka, G.S. Samuelson, Axial flowfield characteristics of reacting and non-reacting flows in a centerbody configuration, Proceedings of the Fall Meeting of the Western States Section of the Combustion Institute, Arizona State University (1981). WSS/CI 81–34. [42] C. Speziale, Turbulence modelling for time-dependent RANS and VLES: a review, AIAA J. 36 (2) (1998) 173–184. [43] J. Frölich, D. von Terzi, Hybrid LES/RANS methods for the simulation of turbulent flows, Prog. Aerosp. Sci. 44 (2008) 349–377. [44] E. Palkin, R. Mullyadzhanov, M. Hadžiabdic, K. Hanjalic, Scrutinizing URANS in shedding flows: the case of cylinder in cross-flow in the subcritical regime, Flow Turbul. Combust. 97 (2016) 1017–1046. [45] D.C. Haworth, S.B. Pope, A generalized Langevin model for turbulent flows, Phys. Fluids 29 (2) (1986) 387–405. [46] D.C. Haworth, S.B. Pope, A PDF modeling study of self-similar turbulent free shear flows, Phys. Fluids 30 (4) (1987) 1026–1044. [47] B.J. Daly, F.H. Harlow, Transport equations in turbulence, Phys. Fluids 13 (11) (1970) 2634–2649. [48] W. Jones, B. Launder, The prediction of laminarization with a two-equation model of turbulence, Int. J. Heat Mass Transf. 15 (2) (1972) 301–314. [49] R.P. Lindstedt, E.M. Váos, Modeling of premixed turbulent flames with second moment methods, Combust. Flame 116 (4) (1999) 461–485. [50] R.I. Issa, Solution of the implicitly discretised fluid flow equations by operator-splitting, J. Comput. Phys. 62 (1) (1986) 40–65. [51] P.K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal. 21 (5) (1984) 995–1011. [52] R.P. Lindstedt, S.A. Louloudi, Joint scalar transported probability density function modeling of turbulent methanol jet diffusion flames, Proc. Combust. Inst. 29 (2) (2002) 2147–2154. [53] R.P. Lindstedt, H.C. Ozarovsky, R.S. Barlow, A.N. Karpetis, Progression of localized extinction in high Reynolds number turbulent jet flames, Proc. Combust. Inst. 31 (1) (2007) 1551–1558. [54] R.P. Lindstedt, H.C. Ozarovsky, Joint scalar transported PDF modeling of nonpiloted turbulent diffusion flames, Combust. Flame 143 (4) (2005) 471–490. [55] J. Xu, S.B. Pope, PDF calculations of turbulent nonpremixed flames with local extinction, Combust. Flame 123 (3) (20 0 0) 281–307. [56] R.P. Lindstedt, S.A. Louloudi, E.M. Váos, Joint scalar probability density function modeling of pollutant formation in piloted turbulent jet diffusion flames with comprehensive chemistry, Proc. Combust. Inst. 28 (1) (20 0 0) 149–156. [57] T.S. Kuan, R.P. Lindstedt, E. Váos, Higher moment based modelling of turbulene enhanced explosion kernels in confined fuel-air mixtures, in: G.D. Roy, S.M. Frolob, R.J. Santoro, S.A. Tsyganov (Eds.), Confined Detonations and Pulse Detonation Engines, Torus Press, Moscow (2003), pp. 17–40. [58] T. Mantel, R. Borghi, A new model of premixed wrinkled flame propagation based on a scalar dissipation equation, Combust. Flame 96 (1994) 443– 457. [59] K.N.C. Bray, Laminar flamelets and the bray, moss, and Libby model, in: N. Swaminathan, K.N.C. Bray (Eds.), Turbulent Premixed Flames, Cambridge University Press, New York, 2011, pp. 4–60. [60] K.N.C. Bray, M. Champion, P.A. Libby, Premixed flames in stagnating turbulence: part V evaluation of models for the chemical source term, Combust. Flame 127 (1) (2001) 2023–2040. [61] F. O’Young, R.W. Bilger, Scalar gradient and related quantities in turbulent premixed flames, Combust. Flame 109 (4) (1997) 682–700. [62] R.S. Barlow, H.C. Ozarovsky, A.N. Karpetis, R.P. Lindstedt, Piloted jet flames of CH4 /H2 /air: experiments on localized extinction in the near field at high Reynolds numbers, Combust. Flame 156 (11) (2009) 2117–2128. [63] P. Kraus, R.P. Lindstedt, Microkinetic mechanisms for partial oxidation of methane over platinum and rhodium, J. Phys. Chem. C 121 (17) (2017) 9442–9453. [64] J. Laufer, The structure of turbulence in fully developed pipe flow, Technical Report, National Bureau of Standards, Washington, DC, USA, 1954. NACA-TR-1174 [65] M.P.I. Forum, MPI: A message-passing interface standard, Int. J. Supercomput. Appl. High Perform. Comput. 8 (3/4) (1994) 159–416. [66] M. Matsumoto, T. Nishimura, Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator, ACM Trans. Model. Comput. Simul. 8 (1) (1998) 3–30. [67] B.B. Dally, A.R. Masri, Modelling of bluff-body recirculation flows, Proceedings of the Twelfth Australasian Fluid Mechanics Conference (1995), pp. 529–532. [68] P. Bradshaw, D. Ferriss, N. Atwell, Calculation of boundary-layer development using the turbulent energy equation, J. Fluid Mech. 28 (3) (1967) 593–616. [69] B. Sakiadis, Boundary–layer behavior on continuous solid surfaces: II. the boundary layer on a continuous flat surface, AIChE J. 7 (2) (1961) 221–225.

L. Tian and R.P. Lindstedt / Combustion and Flame 206 (2019) 51–67 [70] K. Schadow, K. Wilson, J. Crump, J. Foster, Interaction between acoustics and subsonic ducted flow with dump, Proceedings of the 22nd Aerospace Sciences Meeting, Aerospace Sciences Meetings, American Institute of Aeronautics and Astronautics, Reno, NV, U.S.A. (1984). AIAA–84–0530. [71] K. Schadow, E. Gutmark, T. Parr, D. Parr, K. Wilson, J. Crump, Large-scale coherent structures as drivers of combustion instability, Combust. Sci. Technol. 64 (4-6) (1989) 167–186. [72] K. Schadow, E. Gutmark, Combustion instability related to vortex shedding in dump combustors and their passive control, Prog. Energy Combust. Sci. 18 (2) (1992) 117–132. [73] A.L. Birbaud, D. Durox, S. Ducruix, S. Candel, Dynamics of confined premixed flames submitted to upstream acoustic modulations, Proc. Combust Inst. 31 (1) (2007) 1257–1265.

67

[74] D. Durox, T. Schuller, S. Candel, Combustion dynamics of inverted conical flames, Proc. Combust. Inst. 30 (2) (2005) 1717–1724. [75] E.E. Zukoski, Flame stabilization on bluff bodies at low and intermediate Reynolds numbers, California Institute of Technology, Pasadena, CA, 1954 Phd thesis. [76] D. Spalding, Theoretical aspects of flame stabilization: an approximate graphical method for the flame speed of mixed gases, Aircr. Eng. Aerosp. Technol. 25 (9) (1953) 264–276. [77] S. Zhu, S. Acharya, Flame dynamics with hydrogen addition at lean blowout limits, J. Eng. Gas Turbines Power 136 (5) (2014) 051506. [78] R.W. Bilger, The structure of turbulent nonpremixed flames, Symp. (Int.) Combust. 22 (1) (1989) 475–488.