Large eddy simulation (LES) has the potential to predict turbulent combustion phenomena in modern practical combustors. As errors from sub-grid models may be comparable to the numerical errors in the LES approach, mitigating the impact of the numerical errors is as important as constructing accurate sub-grid models. Therefore, a low-dissipative, structure-preserving ROUND (Reconstruction Operators on Unified Normalized-variable Diagram) scheme is tested for the LES of reacting flows in this study. The high efficiency of this scheme is demonstrated by evaluating its accuracy, central processing unit cost, and structure-preserving property by simulating the convection–diffusion process of a passive scalar. Simulations of two bluff body stabilized flames are studied using this scheme. For low turbulence intensity, the current scheme improves the numerical resolution of the instantaneous and averaged flow fields. The prediction of flow statistics is also improved by the ROUND schemes compared to the conventional schemes. Moreover, the ROUND schemes preserve the axisymmetry of the averaged flow better than the conventional schemes for the cases investigated here. For the high turbulence intensity case, the ROUND scheme avoids nonphysical numerical oscillations. The flow and flame statistics obtained using this scheme compare well with measurements. Therefore, this work demonstrates the advantages of using ROUND schemes for LES of reacting flows.

## I. INTRODUCTION

Computational Fluid Dynamics (CFDs) is a valuable tool to investigate and design next-generation combustors with improved stability and reduced emissions. However, turbulent flows, combustion phenomena, and their interactions, which involve a wide range of flow and chemical scales, pose challenges to numerical methods and combustion modeling. To address these challenges, there are a variety of options in terms of methodologies for simulating turbulent reacting flows. On one end of the spectrum, Reynolds-averaged Navier–Stokes (RANS) methods achieve relatively low-cost simulations by modeling all of the spatiotemporal scales of the reacting flow. However, the effectiveness of the RANS approach is highly dependent on the model and the case. On the other end, Direct Numerical Simulation (DNS) aims to resolve all of the spatiotemporal scales, which requires tremendous computational cost and is thus usually prohibitively expensive for simulating practical combustors. Recently, DNS has been applied to simulate practical combustion reactors of medium Reynolds number and has provided high-resolution results for turbulence–flame interactions.^{1,2}

Large Eddy Simulation (LES) can be considered a compromise between effectiveness, accuracy, and computational cost. In the LES approach, only the statistically significant large scales are directly solved while the smaller scales with homogeneous structures are modeled, which is referred to as the Sub-Grid Scale (SGS) range. This separation of scales is achieved by filtering the governing equations. From the perspective of simulation, scale separation is achieved implicitly by the grid element, where scales smaller than the grid spacing are modeled by the physical sub-grid model. Since in turbulent reacting flows, SGSs contain sub-grid turbulence and its interaction with the reaction process, sub-grid models for turbulence and combustion are usually required. However, although there have been a number of developments in the LES approach for turbulent combustion over the decades, the accuracy and predictability of LES still require substantial improvements for industrial applications. As pointed out in the work,^{3} the performance of the LES approach is affected by several factors such as computational setup, numerical methods, chemical kinetics, and SGS models. Therefore, in order to improve the predictive capability of the LES approach, at least the following perspectives are essential:

Consistent and case-independent physical SGS models, which can represent the sub-grid flow and combustion physics in a consistent and accurate manner. The assumptions behind the SGS models should not limit their wide use, including in realistic conditions.

High-fidelity numerical schemes. Ideal high-fidelity schemes should have low numerical dissipation and structure-preserving properties. Low-dissipative schemes should also ensure that most of the resolvable flow scales larger than the grid size are captured. The structure-preserving property is necessary when solving for the transport of scalars such as species mass fraction.

Due to the fact that the physical SGS models are generally linked to the mesh size, the errors from SGS become comparable to the errors from the numerical schemes. As pointed out in Refs. 4 and 5, these two types of errors can interact with one another leading to convergence problems. It was observed in Ref. 6 that the impact of the numerical schemes for LES of reacting flows is significant. Mitigating this impact is essential to improve predictive capabilities. The benefits of low-dissipative schemes in solving combustion problems have also been shown in Refs. 7–11. Especially, the works^{8,9} have shown that high-fidelity schemes are essential to produce correct flame structures in high-speed flows. However, as indicated in Ref. 12, a generic and affordable approach to mitigate numerical errors has yet to be developed.

Indeed, it is challenging to design affordable, low-dissipative, structure-preserving numerical schemes that can be used straightforwardly in the LES of realistic combustors. In the finite volume method (FVM), several high-order schemes have been proposed on arbitrary unstructured meshes to achieve low numerical dissipation, for example, the well-known *k*-exact least-square scheme developed in Refs. 13 and 14. This *k*-exact least-square scheme constructs a polynomial of degree *k* in a stencil consisting of the target cell and its neighbors. To construct high-degree polynomials, this stencil is extended to include more neighboring cells, which results in a wide stencil increasing the complexity of the numerical algorithm and parallelism. As a result, high-order FVM on arbitrary unstructured grids is not widely used for problems of practical interest.

There have also been large numbers of studies to impose the structure-preserving property on numerical schemes. Special techniques such as limiting projection or artificial viscosity have to be designed to overcome the Gibbs limitation, which states that no linear convection scheme of second-order accuracy or higher can be monotonic. In the FVM framework, TVD (Total Variation Diminishing) schemes have been developed in Refs. 15 and 16 and have been widely employed in commercial CFD codes. As pointed out in Refs. 17 and 18, conventional TVD schemes are likely to cause smearing, clipping, and squaring effects, which undermine the structure-preserving property. The WENO (Weighted Essentially Non-Oscillatory) schemes proposed in Refs. 19 and 20 are considered to be more accurate than the TVD schemes. The WENO scheme is able to achieve essentially oscillation-free results by constructing a weighted average of the polynomial interpolations over the candidate stencils. The capability of WENO schemes to resolve critical points has also been improved by recent studies.^{21,22} However, it is not trivial to construct efficient WENO schemes on 3D unstructured grids, as the selection process of candidate stencils increases the algorithmic complexity. The artificial viscosity method^{23} is relatively easy to implement. However, the artificial viscosity method is usually case-dependent and is likely to change the flame-related quantities computed in the simulation. As a result, most 3D unstructured-grid-based CFD codes used for LES employ conventional schemes.

Recently, a new family of schemes named ROUND (Reconstruction Operators on Unified Normalized-variable Diagram) was proposed by Deng.^{18} In the ROUND scheme family, the low-dissipation and structure-preserving properties can be easily imposed on a scheme by rigorously adjusting numerical errors in the unified normalized-variable diagram (UND) that unifies existing second-order schemes. This work aims to test the capabilities and accuracy of this scheme in comparison to the conventional schemes used for LES of reacting flows. First, a pure convection problem and a convection–diffusion process are investigated to verify the accuracy and the structure-preserving property of the ROUND schemes. Then, bluff body stabilized premixed flames experiencing two different turbulent intensities are simulated in order to address the objectives of this study. Comparisons with conventional schemes and experimental measurements are considered to evaluate the efficacy of the ROUND schemes.

## II. LARGE EDDY SIMULATION FOR TURBULENT COMBUSTION

### A. LES governing equations

The Favre-filtered conservation equations for mass, momentum, and total enthalpy are solved using a solver employing the pressure correction method. The sub-grid stress tensor $\tau sgs=\rho \u0304(uu\u0303\u2212u\u0303u\u0303)$ is modeled using the Smagorinsky model.^{24}

*Z*and the reaction progress variable

*c*. The mixture fraction

*Z*is defined using Bilger’s definition.

^{27}For methane-air combustion, the reaction progress variable is calculated as $c=\psi c/\psi cEq$, where $\psi c=(YCO+YCO2)$, and $\psi cEq$ is its equilibrium value for the local

*Z*. Then the Favre-filtered transport equations for these variables are

^{28},

*ν*

_{T}is the SGS eddy viscosity, which is closed using the dynamic Smagorinsky model. The sub-grid Scalar Dissipation Rate (SDR) terms $\chi \u0303Z,sgs$ and $\chi \u0303c,sgs$ require modeling. A linear relaxation model proposed in Ref. 29 is used to model $\chi \u0303Z,sgs$, which gives

*C*

_{Z}= 2 is commonly used. The SDR term $\chi \u0303c,sgs$ is calculated by the algebraic model proposed in Ref. 30 as

_{Δ}is the SGS Damköhler number, and

*C*

_{3},

*C*

_{4}, and

*β*

_{c}are model parameters as detailed in Ref. 30. The other parameters $Kc*$,

*τ*, $SL0$, and $\delta L0$ are obtained from unstrained planar laminar flame calculations.

### B. LES sub-grid combustion model

^{26,31,32}where $\omega \u0307\u0304fp$ signifies the contribution of premixed mode combustion and $\omega \u0307\u0304np$ signifies the contributions from non-premixed mode combustion. $\omega \u0307\u0304cdr$ denotes the cross-dissipation contribution, which can be neglected since this term is an order of magnitude smaller than the contributions from the other two terms. $\omega \u0307\u0304fp$ is modeled as

^{26,31}

*ρ*(

*η*,

*ζ*) are the flamelet reaction rate and density, respectively, which are obtained through one-dimensional unstrained planar laminar premixed flame calculations. The flamelet reaction rate is $\omega \u0307\u2009fp\u2009=\omega \u0307\Psi /\Psi Eq$, where $\omega \u0307\Psi =\omega \u0307CO+\omega \u0307CO2$. Using an unstrained planar laminar premixed flame to obtain the flamelet reaction rate has been validated in many past works.

^{5,25,28,32–34}

*η*and

*ζ*are the sample space variables for

*Z*and

*c*, respectively. The joint probability density function (PDF) is approximated as $P\u0303(\eta ,\zeta )\u2248P\u0303(\eta ;Z\u0303,Z\u2032\u20322\u0303)P\u0303(\zeta ;c\u0303,c\u2032\u20322\u0303)$, where two PDFs are assigned using beta functions. The non-premixed contribution in Eq. (7) is modeled as

## III. NUMERICAL SCHEMES

### A. Discretization of convection terms

*ϕ*is given as

**u**is the velocity vector. To solve the above equation with the finite volume method, the computational domain is decomposed into cell elements. For a given cell

*C*, the spatial convection term of

*ϕ*in Eq. (10) is discretized by integration over the cell volume using the Gauss theorem, following as

*V*

_{c}is the volume of the cell

*C*and

*f*are the mesh faces bounding the cell

*C*.

*ϕ*

_{f}and

**u**

_{f}are the fluid variable and velocity vector at the face

*f*, respectively.

**n**

_{f}is the outward normal-vector of the face

*f*, and

*A*

_{f}is the face area of

*f*. The fluid variable

*ϕ*

_{f}at face

*f*is approximated by the convection interpolation schemes, which will be further discussed in the following sections.

### B. Conventional schemes

^{35}The fluid variable

*ϕ*

_{f}at the cell face

*f*is approximated by the linear central scheme using the adjacent cell values

*ϕ*

_{C}and

*ϕ*

_{D}as

*θ*is defined as the ratio of distances |

**d**

_{Cf}| and |

**d**

_{CD}|,

**d**

_{Cf}| and |

**d**

_{CD}| on unstructured grids are illustrated in Fig. 1. The central scheme is usually considered a good choice for LES. However, as shown in our numerical tests, the sub-grid viscosity of LES may not be high enough to suppress the numerical instability caused by under-resolved energy and numerical interpolations.

*ϕ*

_{f}is given as

*ψ*(

*r*

_{f}) is the flux limiter function and

*r*

_{f}is the gradient ratio of

*ϕ*at the face

*f*, which is defined as

*ϕ*

_{U}denotes the virtual cell value at the upwind cell regarding the face

*f*. As the upwind cell is usually not a real cell element for unstructured grids, the virtual cell value

*ϕ*

_{U}at the upwind cell

*U*is obtained by the following interpolation formulation:

*ϕ*

_{C}is the gradient at cell

*C*. One of the representative flux limiter functions is the vanLeer limiter $\psi (rf)vanLeer$,

^{36}which is written as

### C. Low-dissipative structure-preserving ROUND schemes

*ϕ*

_{f}in physical-variable space, the ROUND scheme proposed in Ref. 18 projects the reconstruction process into normalized-variable space. The normalized fluid variable $\varphi \u0302f$ at the face

*f*, the so-called normalized reconstruction value (NRV), is defined as

*ϕ*

_{C},

*ϕ*

_{D}, and

*ϕ*

_{U}are the same as those in Sec. III B. Once $\varphi \u0302f$ is calculated with the reconstruction operator $\varphi \u0302f=R\u0302(\varphi \u0302C)$,

*ϕ*

_{f}can be obtained with Eq. (18). In ROUND schemes, $R\u0302(\varphi \u0302)$ is devised in the framework that unifies existing second-order polynomial-based and non-polynomial-based schemes, which include TVD, NVD (Normalized Variable Diagram),

^{37,38}WENO,

^{19,20}and THINC (Tangent of Hyperbola for INterface Capturing)

^{39–42}schemes. Through rigorously adjusting numerical anti-dissipation errors in this framework, the overall numerical dissipation of ROUND schemes can be reduced while the structure-preserving property can be maintained. The details of deriving ROUND formulations can be found in the work of Deng.

^{18}Here, the ROUND schemes will be extended to unstructured grids and applied to enhance the LES approach for turbulence flame simulations. A low-dissipative structure-preserving ROUND scheme denoted as ROUND_L is employed to capture the resolvable flow scales of the conservation equations of mass, momentum, and energy. The formulation of ROUND_L is given as

*ω*

_{0}and

*ω*

_{1}are given as $\omega 0=1+\gamma 0(\varphi \u0302C)2\u22124$ and $\omega 1=1+\gamma 1(\varphi \u0302C\u22121)2\u22128$, where

*γ*

_{0}= 12.0 and

*γ*

_{1}= 5.0.

*Z*,

*c*and their variances, a ROUND scheme with rigorously adjusted numerical dissipation is applied. This scheme is denoted as ROUND_A+. Through rigorously adjusting numerical dissipation, ROUND_A+ can simultaneously obtain the TVD condition, low-dissipation, and structure-preserving properties. The additional TVD condition is favorable for solving bounded scalars on unstructured grids. The formulation of ROUND_A+ is written as

*γ*

_{0}= 1100.0,

*γ*

_{1}= 800.0, and

*λ*

_{1}= 0.15. As these parameters have been rigorously adjusted through the framework described in Ref. 18, the parameter values are general to all applications and case-independent. In Sec. IV, results obtained using LES with the above ROUND schemes are discussed.

## IV. NUMERICAL RESULTS AND DISCUSSIONS

In this section, the pure convection problem and the convection–diffusion equation are first solved. Then, the bluff body stabilized turbulent flames subject to two different turbulence intensities are simulated using the LES approach with the ROUND and conventional schemes.

For all the tests, the second-order backward scheme is used for time integration. The diffusion term is discretized using the second-order central difference scheme. The ROUND schemes are implemented within the OpenFOAM^{43} open-source code for LES, which will be used for all the simulations conducted in this study. (Notes: The source code of the ROUND scheme implemented in OpenFOAM is available to interested readers, a copy of which can be obtained by contacting the corresponding author.) For the simulations of the turbulent flames, the Favre-filtered Navier–Stokes equation system is solved using the pressure-based solver with a modified PIMPLE algorithm with density coupling to handle compressibility effects (rhoPimpleFoam solver) available in the OpenFOAM-v7 source code.^{44}

### A. Accuracy tests

The accuracy of conventional schemes and the ROUND schemes is evaluated first by solving a pure convection problem. The computational domain is given by [0, 1] × [0, 1] with periodic boundary conditions applied on all sides. The initial distribution of the passive scalar *ϕ* is given as *ϕ*(*x*, *y*) = sin(2*π*(*x* + *y*)). The convection velocity field is set to **u** = (1, 1). The simulation is carried out on a hexahedral mesh for a period of *t* = 10 s. In order to assess the accuracy and efficiency of the conventional and ROUND schemes, the numerical errors *L*_{1} and *L*_{∞} are calculated by gradually refining the mesh sizes *h* = 1/20, *h* = 1/40, *h* = 1/80, and *h* = 1/160. The corresponding central processing unit (CPU) costs on a single core of the Intel Xeon(R) Gold 6226R CPU at 2.90 GHz are also monitored. The variations of *L*_{1} errors and *L*_{∞} errors with CPU time are presented in Figs. 2 and 3 for different convection schemes. Figure 2 compares the vanLeer and the ROUND_A+ schemes, which are used for updating the bounded scalar. It can be seen that ROUND_A+ significantly reduces the numerical errors with a minor increase in CPU cost. Figure 3 compares two schemes that will be used for discretizing the convection terms in the Favre-filtered Navier–Stokes equation. Similarly, ROUND_L can achieve much smaller numerical errors. Therefore, it is expected that the ROUND schemes will produce higher accuracy than conventional schemes for capturing resolved flow fields.

### B. Convection–diffusion problems of a passive scalar

*ϕ*is

**u**= (1, 1). The Peclet number is Pe =

*UL*/

*D*

_{m}= 1.4 × 10

^{4}. The computation runs up to

*t*= 10.0 s with 200 × 200 uniform hexahedra cells. Figure 4(a) presents the reference solution, which is obtained by only taking the diffusion process into account. The numerical solutions calculated by the TVD vanLeer limiter, the TVD SuperBee limiter, and ROUND_A+ are presented in Figs. 4(b)–4(d).

The results have shown that both the TVD vanLeer limiter and SuperBee limiter distort the shape of the distribution of the passive scalar compared with the reference solution. The vanLeer limiter produces a more diffusive solution due to excessive numerical diffusion. On the contrary, nonphysical anti-diffusion is observed from the result calculated by SuperBee, which indicates that the SuperBee limiter introduces anti-diffusion errors. The numerical solution obtained by ROUND_A+ is close to the reference solution regarding the distribution shape and physical dissipation. Therefore, the ROUND_A+ can better reproduce the convection–diffusion process of the scalar than the TVD schemes. Since convection–diffusion is one of the essential processes in combustion, it can be expected that ROUND_A+ can achieve high-fidelity results compared with conventional TVD schemes.

### C. Numerical simulation of bluff body stabilized premixed flames

#### 1. Case description

The bluff body stabilized methane-air flame is considered here to validate the performance of the LES approach using the ROUND schemes over the conventional scheme. This configuration has been studied experimentally^{45–48} and numerically^{32,33,49} in the past, and it is shown schematically in Fig. 5. A methane-air mixture at 294 K with an equivalence ratio *φ* = 0.59 enters the combustion chamber, which has a dimension of 79 × 79 × 284 mm^{3}. The bluff body burner has a base diameter of *D* = 44.45 mm, a stem diameter of *D*_{st} = 12.7 mm, and an apex angle of *θ*_{a} = 45°. The bluff body is placed at the center of the combustion chamber. A turbulence generator with 3.46 mm diameter holes is placed 58 mm upstream of the bluff body base. A flat velocity profile of *U*_{b,in} = 11.5 m/s is prescribed at the inlet, which is based on the measured mass flow rate to give the required reference bulk-mean velocity *U*_{b} = 15 m/s at the bluff body base, as shown in Fig. 5(a). To overcome the difficulties arising with the boundary conditions for the combustor exit, an additional domain of size 4.5*D* × 4.5*D* × 17.5*D* is included downstream of the combustor exit, as shown in Fig. 5(b). Simulations are undertaken for two different turbulence intensities, i.e., 2% and 22%, for which experimental data are available. These intensities are measured at the bluff body base and in the middle of the jet opening.

The simulation is conducted using the LES approach described in Sec. II, and both the conventional and ROUND schemes are employed. When using the conventional schemes, the central linear scheme is applied to discretize the convection terms in the Favre-filtered Navier–Stokes equation, while the TVD vanLeer limiter is used for solving the strictly bounded variable in Eqs. (1)–(4). When the computation is set up by the ROUND schemes, the ROUND_L scheme is applied for the Favre-filtered Navier–Stokes equation, and the ROUND_A+ scheme is employed to update the strictly bounded scalar variables.

The computational model consists of a hexahedral grid with 3.2 M cells. The finest grid resolution is downstream of the bluff body. For the 2% turbulence intensity case, no inflow turbulence is specified. For the high turbulence intensity case, the velocity boundary conditions use an inflow turbulence generator based on the synthetic eddy method.^{50} The LEMOS inflow generator^{51} is used to prescribe synthetic turbulence by assigning the required rms velocities and turbulence length scale to give a 22% turbulence intensity at the plane in line with the bluff body base. The second-order backward time scheme is used for time marching with an adaptive time step technique and a maximum CFL (Courant–Friedrichs–Lewy) number of 0.4. Pope’s criterion^{52} is used to evaluate how much turbulent fluctuation is resolved by the LES grid. Figure 6 shows the Pope’s criterion for resolved turbulent kinetic energy at the mid-plane for the high turbulence intensity case. It can be seen that more than 90% of turbulent fluctuations are resolved by the LES grid near the downstream of the bluff body.

#### 2. Low turbulence intensity case

The bluff body burner with 2% turbulence intensity is simulated. In this case, it is observed that the flame is thin near the base, and the flame brush becomes thicker moving downwind.^{46} The instantaneous temperature fields in the mid-plane obtained by the two schemes are presented in Fig. 7. It can be seen that the ROUND scheme is able to resolve flame front wrinkles, which are caused by the small-scale vortical structures along the shear layer. However, the conventional scheme yields smoother variation, as in Fig. 7(a). The same behavior can also be observed at other times. As reported in Refs. 7 and 10, resolving small-scale turbulent flame structures requires low-dissipative schemes or very fine mesh grids. Since the same grid is used for both of these schemes, the low-dissipative nature of the ROUND scheme is apparent.

The time-averaged velocity magnitude field in the mid-plane of the bluff-body burner is presented in Fig. 8. A comparison is made between the conventional and ROUND schemes. From Fig. 8, it can be seen that the computed recirculation zone length using the conventional scheme is larger than that obtained using the ROUND scheme. This is because of the excessive numerical dissipation and diffusion in the conventional scheme. Furthermore, one observes that the averaged velocity magnitude field is more axis-symmetric for the ROUND scheme compared to the conventional scheme. This shows that the ROUND scheme can preserve the statistical axis-symmetry better, which is because of the structure-preserving property of the ROUND scheme.

Figure 9(a) compares the centerline variation of measured and computed averaged axial velocity, $\u27e8U\u0303\u27e9$, normalized using *U*_{b}. The centerline variation of measured and computed normalized mean temperature $\u27e8T\u0303+\u27e9$ is compared in Fig. 9(b). The normalized mean temperature is defined as $\u27e8T\u0303+\u27e9=(\u27e8T\u0303\u27e9\u2212Tu)/(Tad\u2212Tu)$, where *T*_{u} = 294 K is the unburnt temperature and *T*_{ad} is the adiabatic flame temperature for *φ* = 0.586. The conventional and ROUND schemes yield similar normalized temperature variations as in Fig. 9(b). However, for normalized axial velocity, the ROUND scheme gives a much less dissipative result. Therefore, the simulation results for the low turbulence intensity case demonstrate that the LES approach with the ROUND scheme improves the numerical resolution as well as the axis-symmetry of the averaged field. However, one observes that there are some differences between the computed and measured statistics. For example, the computed mean axial velocity along the centerline is under-predicted, and the reason for this is unclear at this time. The slight overestimate of the mean temperature near the bluff body base [see Fig. 9(b)] is due to the adiabatic boundary condition used for the bluff body base, while there could be some heat loss in the experiments. This heat loss, or the bluff body base temperature, is not quantified in the experimental studies.

#### 3. High turbulence intensity case

The performance of the ROUND scheme is further tested by considering a case with high turbulence intensity. An increase in the turbulence level will change the flow and combustion characteristics behind the bluff body. For example, the recirculation zone length is observed to decrease when the incoming turbulence intensity increases. Moreover, the high turbulence intensity usually poses a challenge to the stability of numerical schemes. Therefore, a numerically robust approach is required to handle the case of high turbulence intensity.

In this section, the bluff body stabilized flame with 22% turbulence intensity is considered. The computation is performed with different schemes but with the same CFL number and the same grid number. It is found that although the numerical solution can be obtained using the ROUND scheme, the simulation blew up at the initial stage for the conventional scheme. To investigate the causes of the computational failure of the conventional scheme, the instantaneous velocity magnitude and pressure fields at the instant just before the divergence are examined. The instantaneous velocity magnitude fields solved by the conventional and ROUND schemes are presented in Fig. 10, which shows that the conventional scheme substantially over-predicts the velocity magnitude around the base of the bluff body burner. Figure 11 shows the instantaneous pressure field. It is observed in Fig. 11(a) that the conventional scheme produces non-physical pressure oscillations. However, overshoot of velocity magnitude and non-physical pressure oscillations are not observed in the numerical solution obtained by the ROUND scheme. Therefore, it can be concluded that the conventional scheme fails to obtain physically meaningful results for the high turbulence intensity case. It also suggests that the SGS viscosity from the LES model fails to stabilize the numerical solution when the conventional scheme is used for the high turbulence intensity case.

Compared with the low turbulence intensity case, turbulence and flame interactions become stronger in the high turbulence intensity case, which leads to more under-resolved flow energy and a larger gradient flow field such as pressure waves. The second-order linear scheme tends to produce numerical oscillations across the large gradient field when the flow contains a large amount of under-resolved scales. For the low turbulence intensity case, SGS viscosity is able to suppress these numerical oscillations. However, as turbulence intensity increases, SGS viscosity fails to stabilize under-resolved flow energy.

On the contrary, the LES approach using the low-dissipative structure-preserving ROUND scheme is able to produce results without obvious numerical oscillations. The averaged velocity magnitude field and the averaged pressure field in the mid-plane of the bluff body burner are presented in Fig. 12, which is calculated using the ROUND scheme. The statistically axis-symmetrical solutions are obtained. Compared with the low turbulence intensity case, the recirculation zone length of the high intensity case decreases, as shown in Fig. 12(a). In Fig. 12(b), one sees that the thin shear layer near the edge of the base, where the mean pressure decreases, is able to be resolved by the low-dissipative ROUND scheme.

The solution quality using ROUND is further examined by comparing numerical solutions with measurements. Figure 13 compares the centerline variation of $\u27e8U\u0303\u27e9/Ub$. A good agreement between the experimental and numerical results is observed. The radial variations of $\u27e8U\u0303\u27e9/Ub$ are shown in Fig. 14 for nine axial locations from *x*/*D* = 0.1 to *x*/*D* = 2.0. In addition, the radial variations of $\u27e8Ur\u0303\u27e9/Ub$, where *U*_{r} is the radial velocity, are presented in Fig. 15. The measured and computed results agree well, although there is an acceptable small difference. The corresponding variations of normalized mean temperature are shown in Fig. 16. A good agreement is observed when *x*/*D* < 0.8. However, there are some differences between the measurements and the simulation results for *x*/*D* > 1. Reducing these differences may require more accurate LES sub-grid combustion models. To summarize, the LES using the low-dissipative structure-preserving ROUND scheme can provide numerically stable results for the high turbulence intensity bluff body stabilized flame, which is challenging for the LES using conventional or other high-resolution numerical schemes.

## V. CONCLUDING REMARKS

This work proposed and demonstrated an enhanced LES approach using the ROUND scheme, which improves numerical resolution, structure-preserving properties, and numerical robustness for simulating the bluff body stabilized flame. The high efficiency of ROUND is demonstrated by evaluating its accuracy and CPU costs. The structure-preserving property of the ROUND scheme is verified by simulating the convection–diffusion process of a strictly bounded scalar. Then, this work applied LES using the ROUND scheme to simulate premixed combustion in a bluff body stabilized flame with low turbulence intensity and high turbulence intensity, respectively. For the low turbulence intensity case, the current approach greatly improves the numerical resolution of the instantaneous and time-averaged flow fields. Moreover, the axis-symmetry of the averaged field can be better preserved compared with LES using the conventional scheme. For the high turbulence intensity case, LES using the conventional approach fails to obtain physically meaningful results due to non-physical numerical oscillations. However, these oscillations are not observed in the simulations using the low-dissipative structure-preserving ROUND schemes. The simulation results produced using the ROUND schemes are also compared with the measurements, which show good agreement. The results presented here show the efficacy of the ROUND scheme for the premixed flames.

It should be noted that the errors from SGS are as important as the errors from the numerical schemes since two types of errors can interact by adding or partially canceling. Therefore, in some cases, the LES approach with its overall low-order accuracy can also produce acceptable simulation results. A consistent and accurate LES approach, however, should take these two errors into account and mitigate them in a similar order, which will be addressed in future work.

## ACKNOWLEDGMENTS

X.D. and N.S. acknowledge the support of EPSRC through Grant No. EP/S025650/1. J.C.M. and N.S. acknowledge the financial support from Mitsubishi Heavy Industries, Ltd., Takasago, Japan. This work used the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk). The authors are grateful to the EPSRC (Grant No. EP/R029369/1) and ARCHER2 for financial and computational support as a part of their funding to the UK Consortium on Turbulent Reacting Flows (https://www.ukctrf.com).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Xi Deng**: Conceptualization (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **James C. Massey**: Conceptualization (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Nedunchezhian Swaminathan**: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*Towards predictive reacting flow LES*

*Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction*

*28th Aerospace Sciences Meeting*

*Quadratic reconstruction finite volume schemes on 3D arbitrary unstructured polyhedral grids*

*14th Computational Fluid Dynamics Conference*

*Shock capturing with higher-order, PDE-based artificial viscosity*

*18th AIAA Computational Fluid Dynamics Conference*

*A priori*assessment and a posteriori validation

*Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics*

*Spring Technical Meeting*

*Rayleigh/Raman/LIF measurements in a turbulent lean premixed combustor*

*34th Aerospace Sciences Meeting and Exhibit*

52nd Aerospace Sciences Meeting