2.1 Governing equations

The Favre-filtered transport equations for mass, momentum, thermochemical enthalpy (sum of sensible and chemical enthalpies) and the reaction progress variable are written respectively as

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}}\right)=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D70C}}\frac{\mathit{D}\widetilde{\boldsymbol{U}}}{\mathit{D}t}=-\unicode[STIX]{x1D735}\overline{p}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[\overline{\unicode[STIX]{x1D749}}-\left(\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}\boldsymbol{U}}-\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}}\widetilde{\boldsymbol{U}}\right)\right], & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D70C}}\frac{\mathit{D}\widetilde{h}}{\mathit{D}t}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[\overline{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D735}\widetilde{h}-\left(\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}h}-\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}}\widetilde{h}\right)\right], & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D70C}}\frac{\mathit{D}\widetilde{c}}{\mathit{D}t}=\overline{\dot{\unicode[STIX]{x1D714}}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[\overline{\unicode[STIX]{x1D70C}{\mathcal{D}}}\unicode[STIX]{x1D735}\widetilde{c}-\left(\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}c}-\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}}\widetilde{c}\right)\right]. & \displaystyle\end{eqnarray}$$

The overline and tilde operators respectively signify filtered and Favre-filtered variables. The term $\mathit{D}/\mathit{D}t=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ is the material derivative and the symbols $\unicode[STIX]{x1D70C}$ and $\boldsymbol{U}$ respectively denote the fluid density and velocity vector. The molecular shear stress term in (2.2) is $\overline{\unicode[STIX]{x1D749}}=\unicode[STIX]{x1D707}[\unicode[STIX]{x1D735}\widetilde{\boldsymbol{U}}+(\unicode[STIX]{x1D735}\widetilde{\boldsymbol{U}})^{\text{T}}-2/3\,(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\widetilde{\boldsymbol{U}})\,\boldsymbol{I}]$ , where $\unicode[STIX]{x1D707}$ is the molecular dynamic viscosity, which is related to temperature via Sutherland’s law, and $\boldsymbol{I}$ is the identity matrix. The pressure and thermochemical enthalpy are denoted using $p$ and $h$ , respectively. The symbols $\unicode[STIX]{x1D6FC}$ and ${\mathcal{D}}$ denote the molecular thermal diffusivity of the mixture and the mass diffusivity for $c$ in (2.3) and (2.4), respectively.

The reaction progress variable $c$ can be defined using either temperature or appropriate species mass fractions, so that $c$ takes a value of zero and unity in the unburnt and fully burnt mixtures, respectively. It is defined using $\text{CO}$ and $\text{CO}_{2}$ mass fractions following previous studies (Fiorina et al. Reference Fiorina, Baron, Gicquel, Thevenin, Carpentier and Darabiha2003; Ruan et al. Reference Ruan, Swaminathan, Isono, Saitoh and Saitoh2015) and is given here as

(2.5) $$\begin{eqnarray}c=\frac{Y_{\text{CO}}+Y_{\text{CO}_{2}}}{(Y_{\text{CO}}+Y_{\text{CO}_{2}})_{b}},\end{eqnarray}$$

where the subscript $b$ denotes the burnt mixture value. The last term in (2.2) $\overline{\unicode[STIX]{x1D749}}^{R}=(\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}\boldsymbol{U}}-\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}}\widetilde{\boldsymbol{U}})$ is the residual stress tensor; this is decomposed into its isotropic and anisotropic parts respectively as $\overline{\unicode[STIX]{x1D749}}^{R}=2/3\,(\overline{\unicode[STIX]{x1D70C}}\widetilde{k}_{sgs})\,\boldsymbol{I}+\overline{\unicode[STIX]{x1D749}}^{r}$ , where s.g.s. denotes the sub-grid scale. The isotropic part is typically absorbed into the filtered pressure $\overline{p}$ (Pope Reference Pope2000). The anisotropic part is modelled using the expression (Poinsot & Veynante Reference Poinsot and Veynante2012)

(2.6) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D749}}^{r}=-2\,\overline{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D708}_{T}\left[\widetilde{\unicode[STIX]{x1D64E}}-{\textstyle \frac{1}{3}}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\widetilde{\boldsymbol{U}})\,\boldsymbol{I}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D708}_{T}$ is the sub-grid eddy viscosity and $\widetilde{\unicode[STIX]{x1D64E}}=1/2\,[\unicode[STIX]{x1D735}\widetilde{\boldsymbol{U}}+(\unicode[STIX]{x1D735}\widetilde{\boldsymbol{U}})^{\text{T}}]$ is the symmetric strain rate tensor. The sub-grid eddy viscosity is modelled using a localised dynamic Smagorinsky model as (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Lilly Reference Lilly1992)

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{T}=(C_{s}\unicode[STIX]{x1D6E5})^{2}\Vert \widetilde{\unicode[STIX]{x1D64E}}\Vert ,\end{eqnarray}$$

where $\Vert \widetilde{\unicode[STIX]{x1D64E}}\Vert =(2\,\widetilde{\unicode[STIX]{x1D64E}}:\widetilde{\unicode[STIX]{x1D64E}})^{1/2}$ . The Smagorinsky constant $C_{s}$ is assigned following the study by Lilly (Reference Lilly1992). The sub-grid scalar fluxes in (2.3) and (2.4) are modelled using gradient hypotheses with a dynamic evaluation of the turbulent Schmidt number $Sc_{T}$ (Lilly Reference Lilly1992). The remaining term that requires closure in (2.1)–(2.4) is the filtered reaction rate $\overline{\dot{\unicode[STIX]{x1D714}}}$ in the progress variable equation (2.4) and this is described next.

2.2 Combustion closure

The combustion modelling used here is based on unstrained premixed flamelets, which has been described in many books (Libby & Williams Reference Libby and Williams1980, Reference Libby and Williams1994; Echekki & Mastorakos Reference Echekki and Mastorakos2011; Swaminathan & Bray Reference Swaminathan and Bray2011; Poinsot & Veynante Reference Poinsot and Veynante2012). A revised closure based on this concept has been tested for premixed combustion using Reynolds-averaged Navier–Stokes (RANS) and unsteady RANS (URANS) methodologies for laboratory flames (Kolla & Swaminathan Reference Kolla and Swaminathan2010; Swaminathan et al. Reference Swaminathan, Xu, Dowling and Balachandran2011; Darbyshire & Swaminathan Reference Darbyshire and Swaminathan2012; Ahmed & Swaminathan Reference Ahmed and Swaminathan2013, Reference Ahmed and Swaminathan2014) and for practical burners (Ruan et al. Reference Ruan, Swaminathan, Isono, Saitoh and Saitoh2015). This model for the LES of premixed combustion has been developed by Langella & Swaminathan (Reference Langella and Swaminathan2016) and tested for laboratory scale flames (Langella et al. Reference Langella, Swaminathan and Pitz2016a ,Reference Langella, Swaminathan, Williams and Furukawa b ). The use of this model for partially premixed combustion has also been developed by Chen, Ruan & Swaminathan (Reference Chen, Ruan and Swaminathan2017) and applied for gas turbine combustors (Chen et al. Reference Chen, Langella, Swaminathan, Stöhr, Meier and Kolla2019a ,Reference Chen, Swaminathan, Stöhr and Meier b ; Langella et al. Reference Langella, Chen, Swaminathan and Sadasivuni2018a ). The full description of this combustion model is available in these references, for example, see the study by Langella et al. (Reference Langella, Swaminathan and Pitz2016a ). Hence, only a brief description is given here for the sake of completeness.

The filtered reaction rate is modelled as

(2.8) $$\begin{eqnarray}\overline{\dot{\unicode[STIX]{x1D714}}}=\int _{0}^{1}\dot{\unicode[STIX]{x1D714}}(\unicode[STIX]{x1D701})\,P(\unicode[STIX]{x1D701};\widetilde{c},\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2})\,\text{d}\unicode[STIX]{x1D701}=\overline{\unicode[STIX]{x1D70C}}\int _{0}^{1}\frac{\dot{\unicode[STIX]{x1D714}}(\unicode[STIX]{x1D701})}{\unicode[STIX]{x1D70C}}\,\widetilde{P}(\unicode[STIX]{x1D701};\widetilde{c},\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2})\,\text{d}\unicode[STIX]{x1D701},\end{eqnarray}$$

where $\widetilde{P}(\unicode[STIX]{x1D701};\widetilde{c},\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2})$ is the density-weighted sub-grid probability density function (p.d.f.) of the reaction progress variable. The sample space variable for $c$ is $\unicode[STIX]{x1D701}$ and the flamelet reaction rate $\dot{\unicode[STIX]{x1D714}}(\unicode[STIX]{x1D701})$ is obtained from one-dimensional unstrained premixed laminar flame calculations with complex chemistry. The sub-grid p.d.f. is specified using the beta function for given values of $\widetilde{c}$ and its s.g.s. variance $\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}$ . Both of these variables are transported in the LES.

Equation (2.4) is used for determining $\widetilde{c}$ and $\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}$ is obtained using the transport equation

(2.9) $$\begin{eqnarray}\displaystyle \overline{\unicode[STIX]{x1D70C}}\frac{\mathit{D}\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}}{\mathit{D}t} & \simeq & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[\left(\overline{\unicode[STIX]{x1D70C}{\mathcal{D}}}+\overline{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x1D708}_{T}}{Sc_{T}}\right)\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}\right]\nonumber\\ \displaystyle & & \displaystyle +\,2(\overline{\dot{\unicode[STIX]{x1D714}}c}-\overline{\dot{\unicode[STIX]{x1D714}}}\widetilde{c})-2\,\overline{\unicode[STIX]{x1D70C}}\widetilde{\unicode[STIX]{x1D712}}_{c,\mathit{sgs}}+2\,\overline{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x1D708}_{T}}{Sc_{T}}(\unicode[STIX]{x1D735}\widetilde{c}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\widetilde{c}).\end{eqnarray}$$

The quantities $\unicode[STIX]{x1D708}_{T}$ and $Sc_{T}$ are obtained dynamically, as outlined earlier. The reaction-related term is closed using

(2.10) $$\begin{eqnarray}\overline{\dot{\unicode[STIX]{x1D714}}c}=\overline{\unicode[STIX]{x1D70C}}\int _{0}^{1}\left(\frac{\dot{\unicode[STIX]{x1D714}}\,\unicode[STIX]{x1D701}}{\unicode[STIX]{x1D70C}}\right)\widetilde{P}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701},\end{eqnarray}$$

which is similar to (2.8). The filtered reaction rate and the term $\overline{\dot{\unicode[STIX]{x1D714}}c}$ are precomputed using unstrained laminar flame results and these are tabulated as functions of $\widetilde{c}$ and $\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}$ in a look-up table, which is required for the LES.

The s.g.s. scalar dissipation rate for $c$ , denoted as $\widetilde{\unicode[STIX]{x1D712}}_{c,\mathit{sgs}}$ , is influenced by both combustion and turbulence in premixed flames and thus, its modelling should include those influences. The model proposed by Dunstan et al. (Reference Dunstan, Minamoto, Chakraborty and Swaminathan2013), which has been tested thoroughly in previous studies (Gao, Chakraborty & Swaminathan Reference Gao, Chakraborty and Swaminathan2014, Reference Gao, Chakraborty and Swaminathan2015; Langella et al. Reference Langella, Swaminathan, Gao and Chakraborty2015; Langella & Swaminathan Reference Langella and Swaminathan2016; Langella et al. Reference Langella, Swaminathan and Pitz2016a , Reference Langella, Swaminathan, Gao and Chakraborty2017), is used here and is written as

(2.11) $$\begin{eqnarray}\widetilde{\unicode[STIX]{x1D712}}_{c,\mathit{sgs}}={\mathcal{F}}\left[2\,K_{c}\frac{s_{L}^{0}}{\unicode[STIX]{x1D6FF}_{L}^{0}}+(C_{3}-\unicode[STIX]{x1D70F}C_{4}Da_{\unicode[STIX]{x1D6E5}})\left(\frac{2\,u_{\unicode[STIX]{x1D6E5}}^{\prime }}{3\,\unicode[STIX]{x1D6E5}}\right)\right]\frac{\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}}{\unicode[STIX]{x1D6FD}_{c}}.\end{eqnarray}$$

The function ${\mathcal{F}}=1-\exp (-0.75\,\unicode[STIX]{x1D6E5}^{+})$ ensures that the s.g.s. dissipation rate approaches zero when the filter width $\unicode[STIX]{x1D6E5}$ approaches zero and the normalised filter width is $\unicode[STIX]{x1D6E5}^{+}=\unicode[STIX]{x1D6E5}/\unicode[STIX]{x1D6FF}_{L}^{0}$ . The laminar flame speed and its thermal thickness, denoted as $s_{L}^{0}$ and $\unicode[STIX]{x1D6FF}_{L}^{0}$ , respectively, are obtained from the unstrained laminar flame calculations. The symbol $K_{c}=0.79\,\unicode[STIX]{x1D70F}$ is a thermochemical parameter for turbulence–combustion interactions. The other parameters are defined as $C_{3}=1.5\,\sqrt{Ka_{\unicode[STIX]{x1D6E5}}}/(1+\sqrt{Ka_{\unicode[STIX]{x1D6E5}}})$ and $C_{4}=1.1/(1+Ka_{\unicode[STIX]{x1D6E5}})^{0.4}$ , following earlier studies (Dunstan et al. Reference Dunstan, Minamoto, Chakraborty and Swaminathan2013; Langella et al. Reference Langella, Swaminathan, Gao and Chakraborty2015). The symbols $Da_{\unicode[STIX]{x1D6E5}}$ and $Ka_{\unicode[STIX]{x1D6E5}}$ are defined as $Da_{\unicode[STIX]{x1D6E5}}=s_{L}^{0}\unicode[STIX]{x1D6E5}/(u_{\unicode[STIX]{x1D6E5}}^{\prime }\unicode[STIX]{x1D6FF}_{L}^{0})$ and $Ka_{\unicode[STIX]{x1D6E5}}=(u_{\unicode[STIX]{x1D6E5}}^{\prime }/s_{L}^{0})^{3/2}(\unicode[STIX]{x1D6FF}_{L}^{0}/\unicode[STIX]{x1D6E5})^{1/2}$ , respectively. The s.g.s. velocity scale $u_{\unicode[STIX]{x1D6E5}}^{\prime }$ is modelled using a scale-similarity approach as $u_{\unicode[STIX]{x1D6E5}}^{\prime }=C_{q}\sum |\widetilde{\boldsymbol{U}}-\widehat{\widetilde{\boldsymbol{U}}}|$ , where $\widehat{\widetilde{\boldsymbol{U}}}$ is the velocity field obtained using a Gaussian test filter of width $\widehat{\unicode[STIX]{x1D6E5}}\simeq 2\unicode[STIX]{x1D6E5}$ . This model for $u_{\unicode[STIX]{x1D6E5}}^{\prime }$ is used to be consistent with the previous study of the bluff body burner considered here (Langella et al. Reference Langella, Swaminathan and Pitz2016a ), but different models are available, as outlined by Langella, Doan & Swaminathan (Reference Langella, Doan and Swaminathan2018b ).

It was established in past studies that the parameters in (2.11) and their values are closely connected to certain physical aspects of turbulence–combustion interactions and their influences on the scalar dissipation rate. Therefore, these parameters cannot be tuned and further detail on this is given by Kolla et al. (Reference Kolla, Rogerson, Chakraborty and Swaminathan2009), Dunstan et al. (Reference Dunstan, Minamoto, Chakraborty and Swaminathan2013) and Gao et al. (Reference Gao, Chakraborty and Swaminathan2014). The term $\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}/\unicode[STIX]{x1D6FD}_{c}$ in (2.11) comes from influences of flame curvature induced by wrinkling and hence, the scale-dependent parameter $\unicode[STIX]{x1D6FD}_{c}$ for this study is obtained dynamically (Gao et al. Reference Gao, Chakraborty and Swaminathan2015; Langella et al. Reference Langella, Swaminathan, Gao and Chakraborty2015).

The Favre-filtered temperature is calculated as $\widetilde{T}=T_{0}+(\widetilde{h}-\widetilde{\unicode[STIX]{x0394}h_{f}^{0}})/\widetilde{c}_{p}$ , where $\widetilde{c}_{p}$ and $\widetilde{\unicode[STIX]{x0394}h_{f}^{0}}$ respectively represent the effective mixture specific heat capacity at constant pressure, as defined by Ruan, Swaminathan & Darbyshire (Reference Ruan, Swaminathan and Darbyshire2014), and the formation enthalpy of the gas mixture, and the reference temperature is $T_{0}=298.15\,\text{K}$ . For this study, the influence of temperature fluctuations on the specific heat capacity is taken to be small and this assumption is justified by the good comparisons shown in many previous studies (Langella et al. Reference Langella, Swaminathan and Pitz2016a ,Reference Langella, Swaminathan, Williams and Furukawa b ; Chen et al. Reference Chen, Langella, Swaminathan, Stöhr, Meier and Kolla2019a ,Reference Chen, Swaminathan, Stöhr and Meier b ; Langella et al. Reference Langella, Chen, Swaminathan and Sadasivuni2018a ) employing this approach. The filtered enthalpy $\widetilde{h}$ is transported in the LES using (2.3). The mixture density is computed using the state equation $\overline{\unicode[STIX]{x1D70C}}=\overline{p}\widetilde{M}/\Re ^{0}\widetilde{T}$ , where $\widetilde{M}$ represents the Favre-filtered molecular mass of the mixture and $\Re ^{0}$ is the universal gas constant. The three thermochemical quantities of the mixture, $\widetilde{\unicode[STIX]{x0394}h_{f}^{0}}$ , $\widetilde{c}_{p}$ and $\widetilde{M}$ are calculated in a manner similar to (2.8), as described in detail by Ruan et al. (Reference Ruan, Swaminathan and Darbyshire2014), and are included in the look-up table. These thermochemical properties for the flamelet are obtained from the unstrained premixed laminar flame calculations for methane–air combustion. The laminar flames used to build this table are computed using the PREMIX code of CHEMKIN–II (Kee et al. Reference Kee, Grcar, Smooke, Miller and Meeks1985) using the GRI–Mech 3.0 chemical mechanism for methane–air combustion.

The s.g.s. reaction rate closure described in this section uses an unstrained premixed flamelet formulation, which typically assumes that the chemical time scale is shorter than the relevant turbulent time scales. In the context of RANS modelling, it is often questioned whether this closure includes (or responds to) fluid dynamic strain effects. This has led to the proposition of strained flamelet closures in past studies, for example, see the study by Kolla & Swaminathan (Reference Kolla and Swaminathan2010), since the flame stretch effects can reduce the reaction rate which can eventually lead to flame extinction. However, the situation is different for LES modelling, since many of the fluid dynamic time scales, their interactions and mutual influences on the scalar fields ( $\widetilde{c}$ and $\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}$ for premixed combustion and mixture fraction related quantities must be included for partially premixed combustion) are resolved explicitly and captured by the LES equations. This means that the stretching caused by the resolved scales are captured, but the s.g.s. stretching needs to be modelled. Langella & Swaminathan (Reference Langella and Swaminathan2016) showed that including the s.g.s. stretching by using a strained flamelet closure for the filtered reaction rate led to substantially longer flames compared to those in the experiments. However, improved comparisons were obtained when using an unstrained flamelet closure. This is because the strained flamelet closure substantially underestimates the fuel consumption rate. The reasons for these behaviours became apparent through the multi-scale analysis studies by Doan, Swaminathan & Chakraborty (Reference Doan, Swaminathan and Chakraborty2017) and Ahmed et al. (Reference Ahmed, Doan, Lai, Klein, Chakraborty and Swaminathan2018), where it was demonstrated that turbulent eddies smaller than $2\,\unicode[STIX]{x1D6FF}_{L}^{0}$ to $3\,\unicode[STIX]{x1D6FF}_{L}^{0}$ contribute weakly to the overall straining of the flame. Thus, including a model for the s.g.s. stretching may overestimate its effects. If the LES grid resolves scales up to the above eddy sizes, then the unstrained flamelet closure works well. This is demonstrated in many past studies of premixed and partially premixed combustion in laboratory burners (Langella et al. Reference Langella, Swaminathan and Pitz2016a ,Reference Langella, Swaminathan, Williams and Furukawa b ) and model gas turbine combustors (Langella et al. Reference Langella, Chen, Swaminathan and Sadasivuni2018a ; Chen et al. Reference Chen, Langella, Swaminathan, Stöhr, Meier and Kolla2019a ,Reference Chen, Swaminathan, Stöhr and Meier b ) with and without thermo-acoustic instabilities. Furthermore, this s.g.s. combustion closure performed well in capturing the flame root dynamics (Chen et al. Reference Chen, Ruan and Swaminathan2017) including local extinction (Massey, Chen & Swaminathan Reference Massey, Chen and Swaminathan2019).

3.1 Experimental configuration

Figure 1(a) illustrates the schematic of the bluff body burner investigated here, which was studied experimentally by Pan et al. (Reference Pan, Vangsness, Heneghan and Ballal1991a ,Reference Pan, Vangsness, Heneghan and Ballal b , Reference Pan, Vangsness and Ballal1992), Nandula et al. (Reference Nandula, Pitz, Barlow and Fiechtner1996) and Nandula (Reference Nandula2003). Laser Doppler velocimetry (LDV), coherent anti-Stokes Raman spectroscopy (CARS), spontaneous Raman scattering and Rayleigh techniques were employed to measure velocity, temperature and species variations inside the combustor. These measurements were used as validation for the previous numerical study of this burner (Langella et al. Reference Langella, Swaminathan and Pitz2016a ).

A premixed methane–air mixture of equivalence ratio $\unicode[STIX]{x1D719}$ at a temperature $T_{u}=294\,\text{K}$ entered the combustion chamber, which consisted of a square duct with a constant area of $W^{2}=79\times 79\,\text{ mm}^{2}$ and had a length of $L=284\,\text{mm}$ . The conical bluff body at the base of the combustor had a cylindrical stem of diameter $D_{st}=12.7\,\text{mm}$ , a base diameter of $D=44.45\,\text{mm}$ and an apex angle of $\unicode[STIX]{x1D703}=45^{\circ }$ . The premixed methane–air mixture entered the combustor section with a bulk-mean velocity of $U_{b}=15\,\text{m}~\text{s}^{-1}$ at the bluff body base, as shown in figure 1(a). A turbulence generator grid with holes of diameter 3.46 mm was positioned 58 mm upstream of the bluff body base in the experiment. Different turbulence grids were used in the experiment to produce turbulence intensities of $\text{TI}=u^{\prime }/U_{b}=2\,\%$ , $17\,\%$ and $22\,\%$ . These values were measured at a radial location of $r/D=0.55$ at the bluff body base and were used as reference values for the respective experiments. The cases with $2\,\%$ and $22\,\%$ TI are considered for the LES cases, since an extensive set of experimental data is available for these two cases for a thorough model validation.

Figure 1. Schematic of the (a) experimental burner (not to scale) and (b) its computational model (Pan et al. Reference Pan, Vangsness and Ballal1992; Nandula Reference Nandula2003).

3.2 Computational model

The computational domain for the experimental burner is shown in figure 1(b), which is discretised using a structured multi-block mesh with a total of approximately $2.2\,\text{ M}$ cells. The adequacy of this grid for the flow and thermochemical conditions of interest here is tested thoroughly by using two different grids in an earlier study (Langella et al. Reference Langella, Swaminathan and Pitz2016a ), which demonstrated that the flow statistics, including the recirculation zone length, for isothermal flows obtained using this $2.2\,\text{M}$ grid agreed well with measurements. Good comparisons between the LES statistics and measurements for reacting conditions were reported in that study. It is estimated that this grid resolved at least $80\,\%$ of the turbulent kinetic energy in the isothermal flow simulations (Pope Reference Pope2000). Hence, this study employs the same grid, which has refined cells near the bluff body and in the regions where shear layers and filtered flames are expected to be present. The dimensionless wall distance $y^{+}$ can be used to assess the wall refinement of a grid for the viscous sub-layer $(y^{+}\leqslant 5)$ . The experimental data obtained by Nandula et al. (Reference Nandula, Pitz, Barlow and Fiechtner1996) and Nandula (Reference Nandula2003) were used to estimate the wall distance to guide the numerical grid refinement. Approximately two cells are placed within the viscous sub-layer in the grid used for this study. The lengths of the turbulent flames in the experiments exceeded the length of the combustion chamber and hence, difficulties arise with specifying meaningful and numerically stable boundary conditions for the combustor exit. This is overcome by including an additional domain of size $4.5D\times 4.5D\times 17.5D$ downstream of the combustor exit, as shown in figure 1(b). This additional domain allows for the effect of entrained air on the flow exiting the combustor to be captured and also helps to specify clear exit boundary conditions for the computational domain.

A flat velocity profile of $U_{b,\mathit{in}}=11.5\,\text{m}~\text{s}^{-1}$ is prescribed at the inlet, which is based on the measured mass flow rate to give the required reference bulk-mean velocity of $U_{b}=15\,\text{m}~\text{s}^{-1}$ at the bluff body base, as marked in figure 1(a). A small velocity of $U_{\mathit{air}}=0.2\,\text{m}~\text{s}^{-1}$ is specified at the boundary in line with the combustor exit, as shown in figure 1, in order to mimic the ambient air entrainment. Langella et al. (Reference Langella, Swaminathan and Pitz2016a ) showed that the velocity close to the bluff body base could be affected by the heat losses from the recirculation zone to the bluff body base and this loss was reported to be approximately $5$ – $8\,\%$ by Pan et al. (Reference Pan, Vangsness, Heneghan and Ballal1991b ). The previous study by Langella et al. (Reference Langella, Swaminathan and Pitz2016a ) demonstrated that the computed recirculation zone length agreed well with measurements, although this heat loss was excluded in the LES. Hence, this study follows that approach by imposing adiabatic no-slip conditions for the bluff body and walls of the combustion chamber. The lateral walls of the additional domain are specified with slip conditions, while the outlet is specified to have zero streamwise gradients for all variables. The TI at the inlet is specified using the synthetic turbulence obtained using the digital filter technique (Klein, Sadiki & Janicka Reference Klein, Sadiki and Janicka2003) instead of including the turbulence generator used in the experiment.

The numerical simulations are conducted using PRECISE–MB (Anand et al. Reference Anand, Zhu, Connor and Razdan1999), which solves the five transport equations, along with the combustion modelling described in § 2, using the finite volume methodology. The spatial gradients are discretised using second-order accurate difference schemes and the discretised equations are time marched using a second-order scheme (Ferziger & Perić Reference Ferziger and Perić2002). A blended second-order central difference scheme is used to control small overshoots ( ${<}0.6\,\%$ ) and undershoots (not below $-0.002$ ) for the filtered progress variable and its s.g.s. variance. These are the worst values observed over the entire period of the 20 simulations listed in table 1. The velocity and pressure coupling is maintained using the SIMPLEC algorithm (Van Doormaal & Raithby Reference Van Doormaal and Raithby1984). The variables $\widetilde{c}$ and $\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}$ are set to be zero for both the inlet and co-flow boundaries. Following earlier studies (Langella & Swaminathan Reference Langella and Swaminathan2016; Langella et al. Reference Langella, Swaminathan, Williams and Furukawa2016b ), a thermochemical property of the mixture is determined using a mixing rule $\widetilde{\unicode[STIX]{x1D711}}_{mix}=\widetilde{Z}\widetilde{\unicode[STIX]{x1D711}}_{reac}+(1-\widetilde{Z})\widetilde{\unicode[STIX]{x1D711}}_{air}$ , where the subscripts ‘reac’ and ‘air’ denote the values of $\widetilde{\unicode[STIX]{x1D711}}$ for the flamelet and air, respectively. This is used only for the extended domain, as it includes the entrained air. The flamelet thermochemical properties and the reaction-related source terms for the LES are obtained from the look-up table, which has 101 and 51 evenly distributed points for $\widetilde{c}$ and $\unicode[STIX]{x1D70E}_{c,\mathit{sgs}}^{2}$ , respectively, for this study. The procedure allows for the inclusion of mixing or dilution of the burnt mixture with the entrained air. This is achieved by using a filtered transport equation for a passive fluid marker $\widetilde{Z}$ , which is set to be unity and zero, respectively, for the methane–air (combustor exit) and ambient air boundaries. This equation is given as

(3.1) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D70C}}\frac{\mathit{D}\widetilde{Z}}{\mathit{D}t}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[\overline{\unicode[STIX]{x1D70C}{\mathcal{D}}}\unicode[STIX]{x1D735}\widetilde{Z}-\left(\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}Z}-\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}}\widetilde{Z}\right)\right].\end{eqnarray}$$

The last term of (3.1) is modelled using a gradient hypothesis, as outlined previously for (2.3) and (2.4).

All of the simulations reported in this study were run on the Darwin cluster (University of Cambridge High Performance Computing Service). Each node had two 2.6 GHz eight-core Sandy Bridge E5–2670 processors. A constant time step of $7.5~\unicode[STIX]{x03BC}\text{s}$ is chosen to ensure that the Courant–Friedrichs–Lewy (CFL) number does not exceed 0.3. The simulations for $16\,t_{f}$ required 24 hrs of wall clock time on $80$ cores. The flow-through time is defined as $t_{f}=L/U_{b}$ , where $L$ is length of the combustion chamber. The time-averaged statistics are obtained for $8\,t_{f}$ after allowing the transients to escape the computational domain over an initial period of $8\,t_{f}$ .

4.1 Validation

Table 1. Database of simulations and their attributes used for analysis in this study. The bulk-mean turbulence level $u_{b}^{\prime }$ is the surface-averaged value at the bluff body base. The laminar flame speed $s_{L}^{0}$ is in $\text{m}~\text{s}^{-1}$ .

Before addressing the objectives of this study, the modelling framework and the combustion models are validated first using the time-averaged statistics. Experimental measurements are available only for 6 cases of the 20 simulations used for this study, which are listed in table 1. These 20 simulations include both isothermal and reacting flows. The isothermal cases and flames with $\unicode[STIX]{x1D719}=0.59$ are labelled, respectively, using ‘I’ and ‘L’, which are taken from the study by Langella (Reference Langella2016). The relatively richer flames with $\unicode[STIX]{x1D719}=0.8$ and $0.9$ are labelled using ‘R’ and ‘S’, respectively. Cases labelled with ‘1’ and ‘2’ after the letter correspond to turbulence intensities of $2\,\%$ and $22\,\%$ , respectively, which were used in the experiment. The methane–air flames have a Lewis number close to unity.

The values for $u_{b}^{\prime }/U_{b}$ at the bluff body base for these cases are also listed in table 1, which scale approximately as $u^{\prime }/U_{b}\sim 2.3\,u_{b}^{\prime }/U_{b}$ . These values are obtained by varying the axial root mean square (r.m.s.) velocity value at the computational inlet boundary and assigning lateral and longitudinal length scales, as required for the digital filter technique used in the LES (Klein et al. Reference Klein, Sadiki and Janicka2003). The Karlovitz number, defined as the ratio of the chemical time scale to the Kolmogorov time scale, is given as (Swaminathan & Bray Reference Swaminathan and Bray2005)

(4.1) $$\begin{eqnarray}Ka=\frac{\unicode[STIX]{x1D70F}_{c}}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}}=\left\{\left[2(1+\unicode[STIX]{x1D70F})^{0.7}\right]^{-1}\left(\frac{u_{b}^{\prime }}{s_{L}^{0}}\right)^{3}\left(\frac{\unicode[STIX]{x1D6FF}_{L}^{0}}{\unicode[STIX]{x1D6EC}}\right)\right\}^{0.5},\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}$ is the integral length scale at the bluff body base, which is estimated using $70\,\%$ of the flow passage width of $(W-D)/2$ and is equal to approximately $12.1\,\text{mm}$ . The values of $Ka$ listed in table 1, which are based on the turbulence characteristics at the bluff body base, suggest that the combustion conditions in the various cases range from the wrinkled flamelets regime to the thin reaction zones regimes in the turbulent premixed combustion regime diagram shown by Peters (Reference Peters2000). However, the local turbulent time scale, rather than one based on the incoming turbulence level, is more appropriate for understanding the local combustion condition, which can change with downstream distance from the bluff body base. Indeed, Langella et al. (Reference Langella, Swaminathan and Pitz2016a ) showed that there is multi-regime combustion behind the bluff body in cases L1 and L2, where the regimes ranged from the flamelets regime to the distributed flamelets regime. Similar behaviour is observed for the other reacting flow cases listed in table 1. Cases I1, I2, L1 and L2 were validated thoroughly by Langella et al. (Reference Langella, Swaminathan and Pitz2016a ) using velocity measurements (Pan et al. Reference Pan, Vangsness and Ballal1992; Nandula Reference Nandula2003) and scalar fields (Nandula et al. Reference Nandula, Pitz, Barlow and Fiechtner1996; Nandula Reference Nandula2003), where it was demonstrated that the modelling framework detailed in § 2 successfully captured the flame behaviour in the multi-regime combustion. For these reasons and for further validation, cases I1, I2, L1, L2, S1 and S2 are used for analysis in this section, in order to highlight how the equivalence ratio and the TI at the inlet affect the values for $L_{R}$ . These validation cases are marked in table 1 and the additional cases will be used to gain further insights that are required to identify a semi-empirical scaling relation for $L_{R}$ , which will be presented in § 4.3.

Figure 2. Time-averaged centreline axial velocity comparison between the LES results (lines) and experimental data (symbols) for (a) isothermal flow (I1 & I2), and reacting flows with (b) $\unicode[STIX]{x1D719}=0.59$ (L1 & L2) and (c) $0.90$ (S1). There are no experimental data for case S2 ( $22\,\%$ ) shown in (c).

Figure 3. Computed streamlines of the time-averaged velocity in the low and high TI cases of (a) isothermal flows (I1 & I2), and reacting flows with (b) $\unicode[STIX]{x1D719}=0.59$ (L1 & L2) and (c) $\unicode[STIX]{x1D719}=0.90$ (S1 & S2).

Figure 2 compares the measured (Pan et al. Reference Pan, Vangsness and Ballal1992) and computed time-averaged axial velocity along the centreline for five cases; the computational results for case S2 are shown for comparison purposes. For the following analysis, the time-averaged statistics, denoted using the angle brackets, are obtained using the samples collected over $8\,t_{f}$ and the averaging is done both in time and in the azimuthal direction. The negative values for the axial velocity imply reverse flow within the recirculation zone and hence, the values for $L_{R}$ are given by the axial distance of the zero crossing of the normalised axial velocity; this location corresponds to the rear stagnation point of the recirculation zone. The results in figure 2 show that the measured axial variation of the streamwise velocity is captured well in the computations and hence, the values for $L_{R}$ are close to the values seen in the experiment. There are no scalar or temperature measurements available for cases L1 and S1 shown respectively in figures 2(b) and 2(c).

The streamlines of the time-averaged velocity for the six cases are presented in figure 3. The two isothermal cases are shown in figure 3(a) and the streamline patterns are similar for both low and high TI cases, although there is a decrease in $L_{R}$ from $1.45D$ to $1.27D$ when the TI is increased from $2\,\%$ to $22\,\%$ . There is a non-monotonic variation of $L_{R}$ with $\unicode[STIX]{x1D719}$ , or $\unicode[STIX]{x1D70F}$ , for low TI, as shown in figure 3. The value for $L_{R}$ increases to $1.98D$ in case L1 from its isothermal value of $1.45D$ . When the equivalence ratio is increased to $\unicode[STIX]{x1D719}=0.9$ , the value for $L_{R}$ decreases to $1.3D$ for case S1 with $\text{ TI }=2\,\%$ , as shown in figure 3(c). However, this behaviour is different when the TI is increased to $22\,\%$ . The value for $L_{R}$ for case L2 is $1.25D$ , which is very similar to the length of $1.27D$ for the isothermal case I2. The general flow patterns for these cases are also very similar, as shown in figures 3(a) and 3(b). This is consistent with the experimental observation reported by Pan et al. (Reference Pan, Vangsness and Ballal1992), where it was suggested that $L_{R}$ in reacting flows approached its isothermal value when the TI at the bluff body base was increased. However for case S2, the value for $L_{R}$ is $0.7D$ , which is significantly lower than for the isothermal case I2 for $\text{ TI }=22\,\%$ , as seen when comparing figure 3(a) and 3(c).

These variations presented in figure 3 show that $L_{R}$ generally decreases when the TI increases, but it is highly sensitive to changes in the TI for a given $\unicode[STIX]{x1D719}$ in reacting flow. In addition, the study by Bill & Tarabanis (Reference Bill and Tarabanis1986) concluded that the effect of combustion is to increase $L_{R}$ , which does not seem to hold according to the measurements obtained in the study by Pan et al. (Reference Pan, Vangsness and Ballal1992) and the current LES results. The measured values of $L_{R}$ in reacting flows behind a backward facing step were observed to decrease with an increase in $\unicode[STIX]{x1D719}$  (Hong et al. Reference Hong, Shanbhogue and Ghoniem2015), which also contradicts the observation of Bill & Tarabanis (Reference Bill and Tarabanis1986). However, the value of $L_{R}$ is larger than the isothermal counterpart for the low TI cases with moderate heat release. When the heat release is stronger, the value for $L_{R}$ is smaller than the corresponding isothermal value. This intriguing behaviour is shown in figure 4(a), which is to be described in § 4.2. The physical reasoning for this non-monotonic behaviour is unclear and this study attempts to provide that. As noted in § 1, the recirculation zone is in the near-field wake behind the bluff body, which is governed by the momentum transfer. This translates into forces acting on the recirculation zone. Hence, the value of $L_{R}$ is dictated by the various forces acting on the recirculation zone and these forces are not only influenced by turbulence but also by combustion and its nonlinear interaction with turbulence. Therefore, it is imperative to investigate the balance of these forces acting on the recirculation zone in reacting flows.

4.2 Force balance

Figure 4. Variation of (a) the recirculation zone length, and the time-averaged normalised source ${\mathcal{S}}=-\langle (\widetilde{\boldsymbol{U}}-\langle \widetilde{\boldsymbol{U}}\rangle )\rangle \boldsymbol{\cdot }\langle \unicode[STIX]{x1D735}\overline{p}\rangle /(\unicode[STIX]{x1D70C}_{u}U_{b}^{3}/D)$ at (b) $x=L_{R}$ and (c) the maximum width of the recirculation zone, for the six LES cases shown in figure 3 with $\unicode[STIX]{x1D70F}$ .

The previous analysis suggests that $L_{R}$ is influenced by the TI at the bluff body base and the mixture equivalence ratio, or $\unicode[STIX]{x1D70F}$ , but the physical mechanisms for this are unclear. However, Pan et al. (Reference Pan, Vangsness and Ballal1992) speculated that the turbulence production through the interaction of the velocity fluctuations with the pressure gradients could play an important role. This production is given by $-\langle \boldsymbol{u}^{\prime \prime }\rangle \boldsymbol{\cdot }\langle \unicode[STIX]{x1D735}\overline{p}\rangle$ , where $\boldsymbol{u}^{\prime \prime }$ denotes the Favre fluctuation of velocity. It was also suggested that this quantity could be of high importance at the maximum width and at the rear stagnation point of the recirculation zone, since $\boldsymbol{u}^{\prime \prime }$ and $\unicode[STIX]{x1D735}\overline{p}$ are large at those locations. The maximum width of the recirculation zone is taken as twice the distance between the centreline to the farthest radial position on the zero axial velocity contour.

This source can be extracted from the LES results for the two locations in the six cases shown in figure 3 and its variation with $\unicode[STIX]{x1D70F}$ is shown in figure 4. The source term is normalised using $\unicode[STIX]{x1D70C}_{u}$ , $U_{b}$ and $D$ , and is written as ${\mathcal{S}}=-\langle (\widetilde{\boldsymbol{U}}-\langle \widetilde{\boldsymbol{U}}\rangle )\rangle \boldsymbol{\cdot }\langle \unicode[STIX]{x1D735}\overline{p}\rangle /(\unicode[STIX]{x1D70C}_{u}U_{b}^{3}/D)$ . The variation of $L_{R}/D$ with $\unicode[STIX]{x1D70F}$ is also shown in figure 4(a). If this source term is responsible for the variation of $L_{R}/D$ with $\unicode[STIX]{x1D70F}$ , a similar behaviour of ${\mathcal{S}}$ with $\unicode[STIX]{x1D70F}$ must be seen, but it is apparent that no such trend is seen in figure 4. This source term directly contributes to the turbulent stress and hence the turbulent shear force, which is one among many forces acting on the recirculation zone. Thus, it is worthwhile to conduct a force balance on the stationary (time-averaged) recirculation zone. The control volume considered for this analysis is shown in figure 5, which spans the entire length of the recirculation zone and covers the bluff body base in the radial direction. Since the mean structure is axisymmetric, only one half is considered.

Figure 5. Control volume for the force balance of the recirculation zone.

Using (2.2), the stationary form of the mean momentum equation is

(4.2) $$\begin{eqnarray}\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}}\widetilde{\boldsymbol{U}}\rangle +\langle \unicode[STIX]{x1D735}\overline{p}\rangle \approx \langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D749}}_{eff}\rangle ,\end{eqnarray}$$

where $\overline{\unicode[STIX]{x1D749}}_{eff}$ is an effective stress tensor, which includes the molecular and anisotropic turbulent stresses. The averaged molecular stresses are smaller than the turbulent stresses, except on the surface $S_{1}$ marked in figure 5, and this is verified using the LES results.

By applying Gauss’s theorem to (4.2) over the control volume CV1, the force balance is written as

(4.3) $$\begin{eqnarray}\underbrace{\int _{S}\langle \overline{\unicode[STIX]{x1D70C}}\widetilde{\boldsymbol{U}}\widetilde{\boldsymbol{U}}\rangle \boldsymbol{\cdot }\widehat{\boldsymbol{n}}\,\text{d}S}_{\boldsymbol{F}_{i}}+\underbrace{\int _{S}\langle \overline{p}\rangle \,\widehat{\boldsymbol{n}}\,\text{d}S}_{\boldsymbol{F}_{p}}=\underbrace{\int _{S}\langle \overline{\unicode[STIX]{x1D749}}_{eff}\rangle \boldsymbol{\cdot }\widehat{\boldsymbol{n}}\,\text{d}S}_{\boldsymbol{F}_{T}},\end{eqnarray}$$

where $S$ denotes the boundaries of the control volume and $\widehat{\boldsymbol{n}}$ is the unit normal vector that points outwards from each of the boundaries. The subscripts $i$ , $p$ and $T$ for the force vector $\boldsymbol{F}$ denote the inertial, pressure and turbulent shear forces, respectively. These forces can be extracted from the LES results and since the inertial and pressure forces can influence the turbulent velocity fluctuations, and hence the turbulent shear forces, the values of inertial and pressure forces are investigated here. Furthermore, the turbulent shear force can be extracted from these two forces using (4.3). The net inertial and pressure forces acting in the axial and radial directions are listed in table 2 for the six cases analysed previously in figures 3 and 4. The signs for these forces are according to the coordinate system shown in figure 5. The values listed in table 2 show that the pressure forces are significantly larger than the inertial forces in both radial and axial directions. In addition, the radial pressure forces are significantly larger than the axial pressure forces. When the TI is increased for the isothermal cases, I1 and I2, the axial pressure force remains relatively the same, but the radial pressure force decreases by approximately 7 % and this leads to the decrease in the recirculation zone length. It is seen that there is substantial increase in $F_{p,r}$ in case L1 compared to the isothermal case I1, which corresponds to the increase in the value for $L_{R}/D$ . The radial pressure force for case S1 is approximately 35 % lower in comparison to the radial pressure force for case L1 and leads to an approximate 34 % decrease in the value for $L_{R}/D$ , as shown in table 2. All of these trends and the interplay between the radial pressure force and $L_{R}/D$ are shown in figure 6(c), along with the corresponding variation of the axial pressure force, shown in figure 6(b). The axial pressure force is marginally affected when increasing the TI for a given value of $\unicode[STIX]{x1D70F}$ and this force increases with $\unicode[STIX]{x1D70F}$ before decreasing slightly for the flames close to stoichiometric conditions (S1 and S2). There is a direct correspondence between the variations of $L_{R}/D$ and $F_{p,r}$ with $\unicode[STIX]{x1D70F}$ , which suggests that the changes in the recirculation length are dictated by the changes in the radial pressure force in confined bluff body stabilised flames.

Figure 6. Variations of (a) $L_{R}/D$ and the pressure force in (b) axial and (c) radial directions with $\unicode[STIX]{x1D70F}$ . These are deduced from LES results.

Table 2. Net inertial and pressure forces acting on the control volume CV1 in the axial and radial directions.

The radial pressure force is influenced by the heat release from the flame and the corresponding momentum transport is influenced by the turbulent stresses produced through mean shear. In addition, the predominant balance for the radial momentum equation is influenced by these two forces, as shown in (4.3). This is investigated by plotting the shear layers and the flame brush for these cases; these are displayed in figure 7. The edges of the shear layer are marked using 10 % of the maximum positive shear $\unicode[STIX]{x2202}U/\unicode[STIX]{x2202}r$ observed for every $x/D$ location and the non-smoothness seen along the inner edges is caused by the limited sample size available for azimuthal averaging at positions with small $r$ . The time-averaged progress variable values of $\langle \widetilde{c}\rangle =0.1$ and $0.9$ are used to mark the flame brush. The results are shown in the region of $r/D\leqslant 0.7$ for clarity and the recirculation zone length is also marked. The shear layers for isothermal flows are shown in figure 7(a), where it is seen that the shear layer for case I2 is thicker in comparison to case I1, particularly in the region close to the bluff body. This is related to the increased momentum transport that results from the higher TI in case I2, which leads to the decrease in the radial pressure force (see table 2). This causes the recirculation zone to be shorter in case I2, as marked in figure 7(a).

This situation is more complicated for reacting flow, due to the interplay between the induced effects of turbulence and combustion. This interplay depends on the TI and heat release, which is related to the equivalence ratio of the mixture. The averaged shear layer and flame brush isolines are shown in figure 7(b,c). Only the outer edge of the shear layer is influenced by the heat release when the TI is low because the major portion of the flame brush is located close to this edge, as seen in figure 7(b,c). Radial forces are exerted by the combustor wall, due to thermal expansion effects and these forces push the outer edge of the shear layer inward, which is clearly visible in figure 7(b,c) for cases L1 and S1 in comparison to figure 7(a). Increasing the equivalence ratio causes the heat release rate to increase, which leads to the outer edge moving further inward in the near-field region $(x\leqslant D)$ . The flame brush is also thinner in case S1 in comparison to case L1 in the near-field region but widens when moving downstream because of turbulent (diffusive) transport. The flame brush width increases further to accommodate the burning of excess fuel in case S1 compared to case L1. This causes a further inward shift of the shear layer edge, as seen in figure 7(c), which also leads to a drop in the radial pressure force acting in the control volume CV1, as described previously. These interactions cause substantial changes in the net radial pressure forces exerted on the recirculation zone, leading to significant changes in the recirculation zone lengths. The axial pressure force is also influenced by this interplay, but the magnitudes are relatively smaller compared to the radial pressure force, as seen by their values in table 2.

Figure 7. Edges of the time-averaged shear layers (black lines) in (a) isothermal flows and reacting flows. The flame brushes are marked along $\langle \widetilde{c}\rangle =0.1$ and $0.9$ (grey lines) for (b)  $\unicode[STIX]{x1D719}=0.59$ and (c) $\unicode[STIX]{x1D719}=0.90$ cases. These are deduced using LES results.

The interactions of the shear layer and flame brush become more complex when the turbulence intensity is increased to 22 % for cases L2 and S2, as seen in figures 7(b) and 7(c), respectively. The increased turbulent transport widens the flame brush substantially and the majority of the flame brush is inside the shear layer in the near field of case L2, whereas the flame brush is still located in the vicinity of the outer edge for case S2 (with larger heat release). The flame brush is relatively thinner in case S2 when compared to case L2, as the flame in case S2 is closer to stoichiometric conditions. For all of these reasons, the outer edge is shifted further inward in the near-field region for case S2. Moving further downstream, the heat release effects become stronger in case S2, where it is clearly seen that the flame brush is located outside the shear layer in figure 7(c). These strong nonlinear interactions between the flow and thermochemical effects yield further changes to the radial pressure force, which significantly changes the recirculation length.

It has become clear that the changes in the values for $L_{R}$ are predominantly caused by the radial forces exerted on the recirculation zone rather than the turbulence production through the source term $-\langle \boldsymbol{u}^{\prime \prime }\rangle \boldsymbol{\cdot }\langle \unicode[STIX]{x1D735}p\rangle$ , as speculated in the earlier study (Pan et al. Reference Pan, Vangsness and Ballal1992). Further understanding and support to the role of the radial pressure force can be obtained if a scaling relation for $L_{R}$ could be obtained using the insights gained in this analysis.

4.3 Scaling equation for the recirculation zone length

Figure 8. Control volume used to deduce a scaling expression for $L_{R}$ .

Based on the analysis presented in § 4.2 and for the reasons outlined in § 1, it is necessary to consider all of the forces, including those exerted by the confinement of the heat release. Hence, the appropriate control volume used to deduce a scaling relation for $L_{R}$ is shown in figure 8. This modified control volume now spans the whole width of the combustion chamber. It should be noted that the control volume CV2 is not cylindrical, since the outer walls of the combustion chamber form a square duct. The various surfaces of this control volume are marked in figure 8. The surface-averaged radial velocity on $S_{1}$ is observed to be significantly (more than an order of magnitude) smaller than the corresponding axial velocity and hence, the inertial forces acting in the radial direction at $S_{1}$ are neglected. In addition, the turbulent shear force $\int _{S}\langle \overline{\unicode[STIX]{x1D707}}_{T}\rangle \,(\unicode[STIX]{x2202}\langle \widetilde{U}\rangle /\unicode[STIX]{x2202}r)\,\text{d}S$ on surface $S_{1}$ is assumed to be smaller compared to its value on surface $S_{2}$ . This is verified using the LES results, which showed that the magnitude of this force on $S_{2}$ is nearly two orders of magnitude larger. Therefore, the force balance includes only the radial forces acting on surfaces $S_{2}$ and $S_{3,2}$ of the control volume CV2. Applying conservation of mass in the axial direction across CV2 gives

(4.4) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{u}U_{b}\left(W^{2}-\frac{\unicode[STIX]{x03C0}D^{2}}{4}\right)=\unicode[STIX]{x1D70C}_{2}U_{2}W^{2},\end{eqnarray}$$

where the subscript 2 denotes the condition on surface $S_{2}$ and $\unicode[STIX]{x1D70C}_{2}U_{2}=\int _{S_{2}}\!\unicode[STIX]{x1D70C}U\,\text{d}S/W^{2}$ . The angle brackets used to denote the time-averaged quantities are omitted in this and the following expressions for simplicity. Rearranging the mass conservation equation gives

(4.5) $$\begin{eqnarray}U_{2}=U_{b}\frac{\unicode[STIX]{x1D70C}_{u}}{\unicode[STIX]{x1D70C}_{2}}\left(1-\frac{\unicode[STIX]{x03C0}D^{2}}{4\,W^{2}}\right)=U_{b}G\frac{\unicode[STIX]{x1D70C}_{u}}{\unicode[STIX]{x1D70C}_{2}},\end{eqnarray}$$

where $G$ is a geometrical parameter. The radial force balance is now written as

(4.6) $$\begin{eqnarray}\int _{S_{2}}\unicode[STIX]{x1D70C}UV\,\text{d}S+\int _{S_{3,2}}p\,\text{d}S\approx \int _{S_{2}}\unicode[STIX]{x1D707}_{T}\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}r}\,\text{d}S,\end{eqnarray}$$

where $V$ represents the radial velocity. The velocity gradient $\unicode[STIX]{x2202}U/\unicode[STIX]{x2202}r$ is approximated as $\unicode[STIX]{x0394}U/\unicode[STIX]{x0394}r\approx {\mathcal{A}}\,U_{2}/W$ , where ${\mathcal{A}}\approx 8$ because $U$ varies from zero to a maximum over a length of approximately $W/4$ and the maximum velocity is roughly twice the value of $U_{2}$ . Furthermore, the gradient $\unicode[STIX]{x2202}V/\unicode[STIX]{x2202}x\ll \unicode[STIX]{x2202}U/\unicode[STIX]{x2202}r$ on surface $S_{2}$ for this flame configuration. Integrating equation (4.6) and rearranging gives

(4.7) $$\begin{eqnarray}4\,p_{w}WL_{R}\approx 8\,\unicode[STIX]{x1D707}_{T,2}U_{2}W-\unicode[STIX]{x1D70C}_{2}U_{2}V_{2}W^{2}.\end{eqnarray}$$

The pressure on surface $S_{3,2}$ is integrated along the walls of the combustion chamber to give a surface-averaged value of $p_{w}=\int _{S_{3,2}}p\,\text{d}S\,/(4\,L_{R}W)$ . The dynamic viscosity $\unicode[STIX]{x1D707}_{T,2}$ is the average over surface $S_{2}$ . Substituting the expression for $U_{2}$ from (4.5) into the radial force balance in (4.7) and rearranging the resulting expression gives

(4.8) $$\begin{eqnarray}\frac{L_{R}}{W}\simeq \frac{U_{b}G}{4\,p_{w}}\left(\frac{8\,\unicode[STIX]{x1D707}_{T,2}}{W}\frac{\unicode[STIX]{x1D70C}_{u}}{\unicode[STIX]{x1D70C}_{2}}-\unicode[STIX]{x1D70C}_{u}V_{2}\right).\end{eqnarray}$$

The mixture on surface $S_{2}$ will consist of unburnt, partially burnt and fully burnt mixtures, as suggested by the flame brush contours in figure 7 and therefore, the mixture density $\unicode[STIX]{x1D70C}_{2}$ must be in the range $\unicode[STIX]{x1D70C}_{b}\leqslant \unicode[STIX]{x1D70C}_{2}\leqslant \unicode[STIX]{x1D70C}_{u}$ . Hence, the density $\unicode[STIX]{x1D70C}_{2}$ can be expressed as a fraction of burnt mixture density $\unicode[STIX]{x1D70C}_{b}$ using $\unicode[STIX]{x1D70C}_{2}=f\unicode[STIX]{x1D70C}_{b}$ , with $f$ bounded as $1\leqslant f\leqslant \unicode[STIX]{x1D70F}+1$ ; this yields $\unicode[STIX]{x1D70C}_{u}/\unicode[STIX]{x1D70C}_{2}=(\unicode[STIX]{x1D70F}+1)/f$ . Using this expression and noting the flow symmetry on the surface $S_{2}$ , equation (4.8) is written as

(4.9) $$\begin{eqnarray}\frac{L_{R}}{D}\simeq \left(\frac{2\,GU_{b}\unicode[STIX]{x1D707}_{T,2}}{Dp_{w}}\right)\left(\frac{\unicode[STIX]{x1D70F}+1}{f}\right).\end{eqnarray}$$

For isothermal flow, where $\unicode[STIX]{x1D70F}=0$ and $f=1$ , it is observed that $L_{R}/D$ is influenced by the bulk strain $U_{b}/D$ , the eddy viscosity and the wall pressure. Although the TI does not appear explicitly in (4.9), the latter two terms, $\unicode[STIX]{x1D707}_{T,2}$ and $p_{w}$ , will vary with the TI at the combustor entry for a given geometry and $U_{b}$ . Hence, the drop in $L_{R}$ for increasing TI suggests that the rise in the wall pressure is larger than the increase in $\unicode[STIX]{x1D707}_{T,2}$ . Furthermore, it is known that $\unicode[STIX]{x1D707}_{T,2}\sim \unicode[STIX]{x1D70C}_{2}k_{2}^{2}/\unicode[STIX]{x1D700}_{2}$ , where $k_{2}$ and $\unicode[STIX]{x1D700}_{2}$ are the surface-averaged turbulent kinetic energy and its dissipation rate on surface $S_{2}$ , respectively. The values of $k_{2}$ and $\unicode[STIX]{x1D700}_{2}$ would depend on the incoming turbulence and this dependence may be represented using a power law. The pressure distribution on $S_{3,2}$ will also be influenced by the incoming turbulence level and this distribution determines $p_{w}$ .

For reacting flows, there is an additional factor $\unicode[STIX]{x1D70F}$ (related to heat release) that influences the values for $L_{R}$ . The expression (4.9) shows that the recirculation zone length increases linearly with $\unicode[STIX]{x1D70F}$ for a given turbulence level, which is superficial since the heat release will also influence $p_{w}$ and $\unicode[STIX]{x1D707}_{T,2}$ . The experimental results presented in § 4.2 show that $L_{R}$ increases with $\unicode[STIX]{x1D70F}$ for low turbulence levels and moderate values of $\unicode[STIX]{x1D70F}$ . This is because the thermochemical effects are stronger compared to turbulence effects at low TI values. This behaviour also suggests that the increase in $p_{w}$ due to thermal expansion is smaller in comparison to the influences of heat release on $L_{R}$ through $\unicode[STIX]{x1D707}_{T,2}$ and $\unicode[STIX]{x1D70F}$ . The increase in $\unicode[STIX]{x1D707}_{T,2}$ is consistent with the known behaviour of the increase in the turbulence level and its integral length scale across premixed flames with low turbulence intensity (high Damköhler number flames). However, this increase in $L_{R}$ is observed only for a moderate value of $\unicode[STIX]{x1D70F}$ and the recirculation length becomes smaller for higher values of $\unicode[STIX]{x1D70F}$ , as seen in figure 6. This is because the combustion effects coming through the eddy viscosity and $\unicode[STIX]{x1D70F}$ are offset by the rise in $p_{w}$ , as suggested by (4.9), due to larger thermal expansion effects. For higher turbulence levels, the observed variation of the recirculation zone length with $\unicode[STIX]{x1D70F}$ is different. The influences of combustion, on average, are overwhelmed by the effects of turbulence at higher TI and moderate $\unicode[STIX]{x1D70F}$ values, which yield almost no variation of $L_{R}$ with $\unicode[STIX]{x1D70F}$ . It is apparent that these behaviours are contained implicitly in (4.9). The wall pressure and the eddy viscosity are influenced not only by turbulence, but also by dilatation and turbulence–chemistry interactions, as described above. The influences of the latter two thermochemical effects on the spatial evolution and distribution of $k$ and $\unicode[STIX]{x1D700}$ are nonlinear. If the change in the wall pressure is relatively large compared to the change in $\unicode[STIX]{x1D707}_{T,2}$ for a given value of $\unicode[STIX]{x1D70F}$ and TI, then $L_{R}$ will decrease. Hence, a relatively larger thermal expansion coming from near-stoichiometric flames will lead to shorter recirculation zones, irrespective of the TI values, which is also observed in the values listed in tables 1 and 2.

Figure 9. Variation of normalised $L_{R}$ with normalised TI at the bluff body base for the 20 cases listed in table 1. The result of the scaling relation in (4.10) is shown for $a=2.5$ and $b=-0.25$ . The open flame data points shown are from the study by Kariuki, Dawson & Mastorakos (Reference Kariuki, Dawson and Mastorakos2012). The recirculation zone lengths for a backward facing step configuration measured by Hong et al. (Reference Hong, Shanbhogue and Ghoniem2015) for reacting flows of propane–air mixtures with 0 % (black dots) and 50 % hydrogen (grey dots) by volume are also shown.

Deducing an explicit expression for $L_{R}$ as a function of $u^{\prime }$ and $\unicode[STIX]{x1D70F}$ is not quite straightforward, since the effects of combustion on turbulence and related quantities are nonlinear. These nonlinear effects may be represented using a power law for the reasoning presented previously and hence, the scaling expression may be generalised by writing a functional form of the expression as $f(p_{w},\unicode[STIX]{x1D707}_{T,2},\unicode[STIX]{x1D70F})\approx f(u^{\prime }/U_{b},\unicode[STIX]{x1D70F})=(u_{b}^{\prime }/U_{b})^{a+1}(\unicode[STIX]{x1D70F}+1)^{b}$ , where the exponents $a$ and $b$ are constants, and $u_{b}^{\prime }$ is the surface-averaged turbulence level at the bluff body base; these values are listed in table 1. Therefore, the generalised form of the scaling expression is written as

(4.10) $$\begin{eqnarray}\frac{L_{R}}{D}=\widehat{L}_{R}\simeq {\mathcal{C}}\left(\widehat{u}_{b}^{\prime }\right)^{a+1}(\unicode[STIX]{x1D70F}+1)^{b},\end{eqnarray}$$

where ${\mathcal{C}}$ is a constant accounting for a given bulk-mean velocity and the combustion chamber geometry, and $\widehat{u}_{b}^{\prime }=u_{b}^{\prime }/U_{b}$ . Figure 9 shows the results for all of the simulations listed in table 1 that are collapsed using the scaling equation in (4.10). It should be noted that $\widehat{L}_{R}/\widehat{L}_{R,\mathit{ref}}$ is plotted instead of $\widehat{L}_{R}$ , so that the constant ${\mathcal{C}}$ can be eliminated. The subscript ‘ref’ represents the isothermal case I1, which is used as the reference case. The exponents $a=2.5$ and $b=-0.25$ are used for the scaling equation in (4.10) and a line of best fit is shown in figure 9, where it is seen that the proposed expression collapses the data well. The results show that the same behaviour may be seen for higher equivalence ratios, since these results collapse together well, as seen in figure 9. The variation of unscaled $L_{R}/D$ with $u_{b}^{\prime }/U_{b}$ is shown in the inset of this figure.

The values of $\widehat{L}_{R}/\widehat{L}_{R,\mathit{ref}}$ deduced using the results in figure 8 in the study by Hong et al. (Reference Hong, Shanbhogue and Ghoniem2015) for reacting flows over a backward facing step are also scaled as per equation (4.10) and the results are shown in figure 9. Hong et al. (Reference Hong, Shanbhogue and Ghoniem2015) considered propane mixed with 0 %, 30 % and 50 % hydrogen by volume and air mixtures for their experiments; all cases had a TI of approximately 6 %. The scaled recirculation zone lengths shown here are for the two extreme cases with 0 % and 50 % hydrogen. The values for $\unicode[STIX]{x1D719}$ considered in the experiments with the pure propane–air mixture are 0.65, 0.72, 0.79 and 0.88, while the values of $\unicode[STIX]{x1D719}$ for the 50 % hydrogen case are 0.63, 0.67 and 0.72. The scaled values shown in figure 9 range from 1.2 to 0.57 for the pure propane and from 0.88 to 0.53 for the 50 % hydrogen cases. These mixtures have Lewis numbers of approximately 1.9 (pure propane) and 1.6 (mixtures with 50 % hydrogen). The turbulence–chemistry interactions in these non-unity Lewis number flames are substantially different from those in unity Lewis number flames. Furthermore, the recirculation zone behind a backward facing step is constrained by the bottom and side walls and thus the corresponding shear forces cannot be ignored. These additional effects in the flames considered by Hong et al. (Reference Hong, Shanbhogue and Ghoniem2015) may need a different exponent and, perhaps, a Lewis number scaling also. However, the scaled values for these flames shown in figure 9 are of the same order as that given by the scaling relation, which is quite remarkable given the complexity involved in this problem.

4.4 Scaling equation for the recirculation zone length for open flame

An illustration for a bluff body stabilised flame without confinement is shown in figure 10, where it is seen that an outer shear layer will form at the top edges of the walls. The appropriate cylindrical control volume is extended beyond this outer shear layer into the air entrainment region with zero streamwise velocity. If the observations from § 4.3 regarding the role of radial forces acting on the recirculation zone are extended, then the radial force balance in (4.6) for this control volume is now written as

(4.11) $$\begin{eqnarray}4\,p_{\infty }W_{\infty }L_{R}\approx 8\,\unicode[STIX]{x1D707}_{T,2}U_{2}W_{\infty }-\unicode[STIX]{x1D70C}_{2}U_{2}V_{2}W_{\infty }^{2}.\end{eqnarray}$$

Figure 10. Control volume for an open bluff body burner configuration.

Applying the mass balance across the cylindrical control volume gives

(4.12) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{u}U_{b}\left(\frac{W^{2}-D^{2}}{4}\right)+\unicode[STIX]{x1D70C}_{\mathit{air}}\hat{\unicode[STIX]{x1D6FC}}U_{b}W_{\infty }L_{R}=\unicode[STIX]{x1D70C}_{2}U_{2}\frac{W_{\infty }^{2}}{4},\end{eqnarray}$$

where the entrained mass flow rate is taken to be $\unicode[STIX]{x03C0}W_{\infty }L_{R}\unicode[STIX]{x1D70C}_{\mathit{air}}\hat{\unicode[STIX]{x1D6FC}}U_{b}$ , with $\hat{\unicode[STIX]{x1D6FC}}$ representing an entrainment coefficient for air that accounts for scaling the velocity to $U_{b}$ and the area change to include the entrainment at the surface of the control volume in line with the bluff body base. After taking $\unicode[STIX]{x1D70C}_{air}\approx \unicode[STIX]{x1D70C}_{u}$ and then using $\unicode[STIX]{x1D70C}_{u}/\unicode[STIX]{x1D70C}_{2}=(\unicode[STIX]{x1D70F}+1)/f$ , equation (4.12) is rearranged for $U_{2}$ and is written as

(4.13) $$\begin{eqnarray}U_{2}=U_{b}\frac{\unicode[STIX]{x1D70F}+1}{f}\left(G+\frac{4L_{R}}{W_{\infty }}\hat{\unicode[STIX]{x1D6FC}}\right),\end{eqnarray}$$

where the geometrical parameter is expressed as $G=(W^{2}-D^{2})/W_{\infty }^{2}$ . Substituting the expression for $U_{2}$ into (4.11) and rearranging for $\widehat{L}_{R}=L_{R}/D$ gives

(4.14) $$\begin{eqnarray}\widehat{L}_{R}\simeq \frac{U_{b}^{2}G^{\ast }}{4\,p_{\infty }}\left(\frac{8}{W_{\infty }}\frac{\unicode[STIX]{x1D707}_{T,2}}{U_{b}}\frac{\unicode[STIX]{x1D70F}+1}{f}+\unicode[STIX]{x1D70C}_{u}\frac{|V_{2}|}{U_{b}}\right){\mathcal{E}},\end{eqnarray}$$

where ${\mathcal{E}}$ is a correction term accounting for the air entrainment, which can lead to some cross-stream velocities and hence, the velocity $V_{2}$ is retained in (4.14).

Equation (4.14) is similar to (4.8) besides the change in variables for the geometrical parameter, defined for the open flame as $G^{\ast }=G(W_{\infty }/D)$ , the pressure $p_{\infty }$ and the width $W_{\infty }$ . The influence of heat release and the TI on the entrained air is signified by the factors $(\unicode[STIX]{x1D70F}+1)/f$ and $\unicode[STIX]{x1D707}_{T,2}$ appearing in the first part. Thus, the justifications given earlier can be used to introduce a functional dependence on the TI and $\unicode[STIX]{x1D70F}$ and hence, equation (4.14) can be written in a form similar to (4.10) but perhaps with different values for the exponents $a$ and $b$ . This suggests that the values for $L_{R}/D$ for open flames are also likely to follow the scaling shown in figure 9. This is demonstrated in figure 9 by including the values of $L_{R}/D$ for four open flames, which are also normalised using the corresponding isothermal case in that experimental study. The values for the recirculation zone lengths are taken using the centreline particle image velocimetry (PIV) measurements made in the study by Kariuki et al. (Reference Kariuki, Dawson and Mastorakos2012). These results suggest that the scaling given by (4.10) works well for both confined and open flames.