Numerical Study of the b-Ξ Flame Wrinkling Combustion Model in Oracles Test Rig

AbstrAct: Numerical modeling of premixed combustion is important for a wide range of machines and systems focusing on compliance with the increasing pollutants reduction requirements. However good industrial numerical combustion models need a practical requiring, in this way, a balance between speed and accuracy. The flamelet models are suitable for this purpose providing a decoupling of the reactive and fluid dynamic problems, and an important model of this family is the b-Ξ flame surface wrinkling model. A specially challenging experiment to test this combustion model is the ORACLES test rig whose two independent parallel inlet channels consistently influence the turbulent combustion, injecting fuel and oxydizer at different equivalence ratios. The b-Ξ flamelet combustion model is known by the sensibility to numerical schemes and boundary conditions and, based on this, the present study proposes to investigate the coupling with the important SST k-ω turbulence model and achieve good balance among accuracy, boundedness, stability and efficiency using the ORACLES experiment.


IntroductIon
A truly instigating and challenging phenomenon of numerical modeling is the study of partially premixed turbulent combustion cases in closed environment.There are several factors responsible for that, for example, a bigger coupling between chemistry and turbulence on premixed turbulent combustion in relation to the non-premixed problem (Bilger et al. 2005).
In industrial applications, specially gas turbines and internal combustion engines, the study of premixed turbulent combustion becomes a central point.In this combustion type the reactive wave propagates towards the reactant molecular mixture, initially in a deflagration speed (Lipatnikov 2012).The large eddies wrinkle the flame and the deformations consequently increase the speed, and, as soon as the eddies scale declines, reaching the flame thickness, they may penetrate and modify the flame structure (Gulder 1984), which can even quench the flame.An interesting feature of premixed flames is that its flammability and extinction limits depend on the flame front wrinkling, buoyancy (Qiao et al. 2008), equivalence ratio, heat losses, among other factors.
It must be noticed that there is a great production of thermal energy inside the flame (Lipatnikov 2012), and, despite this, behind the flame, the speed grows about 20 times the front flame value due the great density difference between the reactants and products in the flame nearby, as says Moss (1980).
Turbulent reactive internal flows are usually more difficult to simulate than the external ones due many reasons as wave reflection on the walls, flow separation and existence of boundary layers, requiring a better mesh treatment and refinemen (Hirsch 2007).Consequently it is very important to choose adequate numerical schemes and turbulent models to reduce numerical errors propagation and also the computational costs.To build consistent experimental basis for numerical calculations on internal flows, many studies have been developed, as those by Moriyoshi et al. (1996).Pitz (1981), focusing on gas turbines applications, shows that a sudden expansion in the reactant flow can lead to a premixed flame stabilization and a better control for the resultant pollutant emissions.With the Large Eddy Simulations (LES) popularisation many validation needs arose, including the requirement to validate unsteady problems with flame anchoring in flow separation.Focusing on building an experiment for this validation type, a pioneering and interesting study was developed by Besson (2002) consisting of a turbulent fully or partially premixed combustion chamber with opposite sudden expansions on a perfectly controled conditions.The combustion chamber, with a very large aspect ratio and a rectangular profile, was built with two parallel injection ducts for fuel-oxydizer mixture injection with independent mixture ratio from each other.Ahead of these inlet rectangular ducts, two opposite backward steps were conceived for the flame stabilization and a long exhaust duct was instrumented in some stations for LES data accquisiton as the frequency energy spectrum.The author has showed as results that, with the same inlet parameters, a visible asimmetry exists in the inert flow as well as the existence of the strong deterministic component in some configurations.In this test bench, named One Rig for Accurate Comparisons with Large Eddy Simulations (ORACLES), Nguyen (2003) focuses on the construction of a database for many equivalence ratios and extintion characteristics, continuing Besson's study.
A complete mathematical model to describe a reactive turbulent flow is constituted by conservation equations for mass, energy and mass fractions.The most tradicional models are based on the idea of simplifying the complete chemistry kinectics aiming to attenuate the enormous computational cost for direct resolution of these equations (Weller 1993).Some reduced chemistry methods were developed recently with the aim of increase this chemistry simplification, as the Intrinsec Low Dimensional Manifolds (ILDM), which intends to eliminate slow time scales, reducing the dimensions of the system (Maas and Pope 1992).In this method all the data is classified in an in situ table, even the chemical reaction source term and the mass fractions (Pope 1997; Kröger et al. 2010).The Eddy Break-Up models, or EBU models, consider an infinitely fast chemistry for combustion processes entirely controled by turbulent mixture (Spalding 1971), where some of them incorporate curvature and local deformation effects.After the EBU models advent, Probability Density Function (PDF) combustion models came to avoid some gradient transport assumptions and insert more information in correlations and means (Law 2006).Nevertheless, these models are resolved in a multidimensional way by expensive methods as the Monte Carlo ones.By this reason, presumed PDF models were developed combining turbulent mixture effects and chemistry kinectics, where a PDF representing the reaction rate and kinectic effects is included by simplified mechanisms (Leoni 2010).
Flamelet models use the same concept as the presumed PDF models and many of these have an energy equation very similar to the one used in EBU models.The idea is that all reactive field is represented by many local laminar flames, not perturbed by the turbulent field.A problem in some EBU models is associated to a near wall poor performance when simulating complex reactive flows.For this reason Marble and Broadwell (1977) developed a flamelet model for tubulent flow with laminar flames that grow according to the flame front distortions and the wall distance, with transport equations for the progress variable and flame area.This area, the flame area density Σ, is given per unit of volume and the model formulation is a little bit different from the EBU ones.The flame wrinkling model developed by Weller (1993), in turn, is an alternative for Σ models using the flame wrinkle density factor Ξ, that is the laminar flame wrinkled area per unit of area, solved on the flow direction, which introduces more phisical realism in the formulation.However, one drawback is that the b-Ξ flamelet combustion model is very sensitive to boundary conditions as well as numerical schemes and Ξ models (CFD-Online 2011;Weller 1993;Kröger et al. 2010).

numerIcAl method
The used model for this investigation is the Weller b-Ξ flame surface wrinkling combustion model (Weller 1993), implemented in the premixed combustion solver XiFoam.This code is part of the OpenFOAM (ESI-OpenCFD 2015), that is a C++ toolbox for scientific simulation created originally to resolve Computational Fluid Dynamics (CFD) problems using finite volume discretization (Jasak 1996).The basic parameter to define the premixed flame using the present model is the progress variable, c, which indicates the reaction progress towards the unburned gas.The progress variable complement, the regress variable, models the flame front propagation and is calculated by Eq. 1: where: P is pressure; Φ is the equivalence ratio.For propane, the empirical coefficients are W = 0.446; η = 0.12; ξ = 4.95; α = 1.77; β = -0.2.The specific heat, entropy and enthalpy are calculated by polynomial approximation using the JANAF tables.
The flame wrinkle factor can be calculated by the transport equation described in Eq. 7: where: b = 1 for the fresh gas; b = 0 for the burned one; T is temperature; the subscripts b and u refer to burned and unburned gases, respectively.
The transport equation for this variable is given by Eq. 2: where: Sc = μ/ρD is the Schimidt number; Sb is a source term; ρ is the density; u is the velocity; D is the molecular diffusion rate.The index t refers to turbulence; μ is the viscosity, obtained by the Sutherland law, shown in Eq. 3: The empirical values are: A s = 1.67212 x 10 -6 and T s = 170.672.The source term for the regress equation on the right hand of Eq. 2 is modeled by: where: S is the laminar flame speed; Ξ is the flame wrinkle factor, which is represented by an algebraic equation (Eq.5): where: u' is the turbulence intensity; R η is the Reynolds number based on Kolmogorov length.
Still detailing Eq. 4, the laminar flame speed S u is calculated by the Gulder correlation and is given by Eq. 6: where: U is the average flame surface velocity; σ is the strain rate; the subscript s refers to surface.Thus in Eq. 7 G is given by: and R is given by: where: τ η is the Kolmogorov time scale.The turbulence model used is the SST k-ω, described by Menter (1994), as well the wall functions.The SST k-ω was chosen by the fact that it is a robust model that can provide good results with flow separation under adverse pressure gradients.
Among many numerical schemes for temporal term discretization available in OpenFOAM, two of them were chosen in this paper: the tradicional first-order Euler method and the second-order Backward Diferencing Scheme.The latter computes a temporal derivative using Taylor expansion starting in the current time going back two sucessive time steps (Maric et al. 2014).If δt is a variable time step, the code implementation is: where: ϕ c is a cell-centered field; the superscript o refers to the old time step and oo refers to the old-old time step.Thus, the other coefficients in Eq. 10 are: (1) To investigate the behavior of the convection terms in this simulation two schemes are chosen for divergence terms.Limited Linear, as named in OpenFOAM, is a second-order bounded Gauss method with a Sweby flux limiter function (Sweby 1984).In order to compare the results it is chosen the first-order/second-order Linear Upwinding Differencing Scheme (LUDS) that uses the flow direction to choose the interpolation points, merging the boundedness of the upwind scheme with the accuracy of the central differencing scheme.

Geometry And mesh
The internal fluid volume of the ORACLES combustion chamber is represented in the present study by the geometry shown in Fig. 1.A detail of the sudden expansion can be viewed in Fig. 2.
All the numerical studies are carried out with a pseudo bidimensional (one element in z axis) mesh with approximately 80,000 fully hexaedral elements.A refinement detail of the junction area between the two gas inlets is shown in Fig. 3.

test CAses
The numerical studies in this paper are compared with experimental results obtained from the ORACLES test rig (Besson 2002).The chosen experimental cases are shown in Table 1.The number after the dot refers to the upper inlet, named 1, and to the lower inlet, named 2. U is the velocity measured on the x axis and "axial" concerns to this velocity measured in the half height of the respective inlet channel, in a position 170 mm far from the origin, indicated in Fig. 2. The equivalence ratio Φ is present in the reactive case r2 and U deb refers to the inlet average speed.
The inert case i1 has Reynolds numbers approximately the same of the ones obtained with the reactive r2 case, as it can be seen in Table 1.This correspondece is chosen purposely to test the turbulence model in this Reynolds range for i1 before test the combustion model coupled with the turbulence model in the case r2.The numerical schemes and k -ω configurations are tested in higher Reynolds numbers by case i2.
Out of a large number of possible numerical schemes and combustion model configurations, five comparations are chosen in this paper, as shown in Table 2, with modifications in the temporal discretization scheme and divergence term discretization as well the used Ξ modelling.In the configurations 3, 4 and 5 it is used a Linear Upwind scheme for velocity (U), pressure (p), kinetic energy (k) and specific dissipation (ω) and Gauss Limited Linear for the other variables, whereas in the configurations 1 and 2 the predominant divergence scheme     for all variables is the Gauss Limited Linear.The comparison between configurations 4 and 5 leads to the analysis of the flame wrinkling model and, consequently, the flame shape and its characteristics.

BoundAry Conditions
The outlet pressure condition is atmospheric, the walls are isolated, the inlet temperatures are setted as 293 K and b and Ξ are considered unitary in inlet.The mixture fractions follow Table 1 in inlets.The ignition points are equally spaced, located in x = 0.005 mm and away from horizontal walls 0.005 mm.The ignition time takes 0.05 s.Wall functions are used for ω and k, detailed in Menter (1994).
Table 3 shows the detailed boundary conditions chosen for the numerical simulation; notice that it is a bidimensional simulation.Fully developed profiles are imposed as initial conditions for k, ω and U fields by mapping the results of simulations performed in the OpenFOAM standard solver BoundaryFOAM (ESI-Open CFD 2015).This is a one-dimensional steady-state OpenFOAM olver code to simulate fully developed turbulent flows.

results inert CAses
In Fig. 4, it is shown the velocity field for case i2 using configuration 1.A non-physical behavior is observed due the Euler scheme use in the temporal term discretization, which creates an irreal wake in the inlet separating plate.The error is analyzed in Fig. 5, where it is shown the horizontal velocity profile at the origin.H is the backward step heigth and U deb is found in Table 2.
The wake is caused by the amplification of the pressure difference between the two inlets and the profile for Euler scheme shown in Fig. 5 does not represents the reality.The difference between configurations 2 to 3 is just the numerical scheme for the divergence term, changed to LUDS.Both results do not show visible difference in the total results; however configuration 3 converges in 0.7 s while configuration 2 converges in 1.5 s.Considering these results configuration 3 is the chosen for inert simulations.
Figure 6 shows the flow velocity field for i2 with flow streamlines.It can be seen in this same figure black vertical lines where the variables are analyzed and compared to experimental data, which are named as stations.By observing the flowfield it can be noted an asymmetry between the upper and lower recirculation zone, as shown by Besson (2002).In Fig. 7 the numerical and experimental horizontal velocity profiles are compared for the stations shown in Fig. 6 and in Fig. 8 the vertical velocity profiles in the stations 3, 4 and 5 are shown.The weak agreement of the vertical velocity profiles in this comparison shows that the turbulence model is not able to capture the recirculation zones using low velocity fluctuation levels but it provides very good agreement for horizontal velocity profiles.
Figures 9 and 10 show the time history for, respectively, the vertical and horizontal velocities in three points in the geometry.All the three measuring points are located 0.005 m far from the superior horizontal wall, in the stations 3, 5 and 6.The total simulation time for all simulations conducted in this study is 2 s and it is shown in these figures that it is sufficient for the stablishment of a fluctuation pattern.The low frequency waves, with wavelength of approximately 0.5 s, are due numerical errors because the experimental coherent movements from 50 to 100 Hz (Besson 2002).In Fig. 11 it can be noted the average vertical velocity evolution along two horizontal lines, located in Y/H = -2 and Y/H = 2.The recirculation zones lenghts are defined where the velocity reverts the direction in these lines.This results agree with Abbott and Kline (1962), which measures the upper and lower recirculation zone lengths in a square profile duct with sudden expansion.For an expansion ratio of 1.8, approximately the value for ORACLES test rig, the author obtained the length of 10 H and 4 H.In this study these values are 12 H and 4 H, representing a good agreement.By analysing case i1 now, in a lower Reynolds range, the behavior for velocity fields is similar to i2, as presented in Fig. 12 for horizontal velocity profiles in stations 3, 4, 5 and 6.In Fig. 13 it is seen the vertical velocity profiles in stations 4 and 6 showing the turbulence model error in the capture of the recirculation zone in this low velocity level.In Fig. 14 k and ω fields show the recirculation zones, the boundary    conditions used and visible asymmetry for case i1 in the same way as case i2.
reACtive CAse With the same boundary conditions of the i1 case, except the mass fraction boundary condition, case r2 is the reactive case i1.Two configurations are tested, one of them with the algebraic flame wrinkling model and another using a transported model.To define the best configuration, velocity profiles comparison is conducted in the burned gas acceleration zone, where the stations 4, 5 and 6 are located.
In Fig. 15 it is shown a comparison between the algebraic model, transport model and experimental results, giving large advantage for transported model when defining the flame length and correct velocity fields.By this finding, configuration 5 is chosen for obtaining all the results for reactive simulation.
Figure 16 shows the temperature distribution for r2 case.In Figs. 17 and 18 the temperature field and a front flame photography taken during Besson (2002) i2 experiment are shown, respectively, pointing out the similitude with the numerical simulation.
In Table 4 the temp eratures, numerically and experimentally obtained, are compared.Higher values are obtained in both cases, probably due the use of pure propane in this study, in contrast with the experimental study, which uses a mixture with approximately 85% of butane.Other ideal supposition is the null heat transfer through the walls.Figures 19 and 20 show, respectively, the horizontal and vertical profiles in the same stations 3, 4 and 5 adopted for i2 simulations.In Fig. 19 it is shown a noteworthy concordance between the numerical and experimental results using configuration 5, detailed in Table 2.The best results for vertical velocity profiles in Fig. 20 are due to high burned gases velocity relative to the combustion phenomena.

dIscussIon
In this study it is described the numerical simulation of some experiments carried out in the ORACLES test rig using b-Ξ flame wrinkling combustion model and SST k-ω turbulence model focusing on the coupling adjustment between them.This experiment is particularly chalenging for numerical schemes testing due the pressure flutuations among the parallel inlet channels requiring special atention to the temporal discretization schemes.It should be noted that backward scheme for temporal term has reached the best results and maintains plausible behavior whereas the Euler scheme amplifies the pressure oscilations and causes an irreal wake that starts in inlets separator blade.The results obtained from the reactive case show that the combination of the backwad Euler for temporal term and Linear Upwind for divergence terms is valid for reactive and inert flows in the Reynolds range tested in this study, and presents a better behavior than the inert case due the higher velocities.The transported Ξ model gives a large advantage for the reactive simulation where the algebraic model fails in determining correct velocity fields.With these results LES turbulence models associated with b-Ξ flame wrinkling combustion model can be used in the ORACLES test rig since there is a large database of temporal statistics available for this experiment.
Numerical Study of the b-Ξ Flame Wrinkling Combustion Model in Oracles Test Rig

Figure 4 .
Figure 4. Numeric error in the central wake using Euler scheme for temporal term.

Figure 8 .
Figure 8. Vertical velocity profiles in stations in case i 2.

Figure 12 .
Figure 12.Horizontal velocity profiles in case i 1.

Figure 13 .
Figure 13.Vertical velocity profiles in case i 1.

Figure 14 .
Figure 14.k and ω fields for case i 1.

table 1 .
Experimental cases for numeric comparison.