CFD Simulations AC2-09

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

SANDIA Flame D

Application Challenge AC2-09   © copyright ERCOFTAC 2024

Overview of CFD Simulations

All the calculations presented below were obtained within the MOLECULES FP5 project Contract N° G4RD-CT-2000-00402 by the team of the Institute of Thermal Machinery, Częstochowa University of Technology. The computations were performed with BOFFIN-LES code developed at Imperial College by the group of Prof. W.P. Jones. BOFFIN-LES computer code utilizes a boundary conforming general curvilinear coordinate system with a co-located storage arrangement. It incorporates a fully implicit formulation and is second order accurate in space and time. For the convection terms an energy conserving discretization scheme is used and matrix preconditioned conjugate gradient methods are used to solve the equations for pressure and velocity etc. The CFD simulations are all LES predictions with various subgrid scale models and turbulence/combustion interaction approaches and neither RANS nor URANS methods are studied in this document.

In the LES calculations two models of turbulence/combustion interaction were applied: steady flamelet model and simplified Conditional Moment Closure (CMC) neglecting the convection term in physical space (The CMC module was developed by Prof. E. Mastorakos from Cambridge University). In both cases the standard subgrid-scale (SGS) Smagorinsky model was used. Then in order to evaluate the importance of the subgrid-scale models the LES calculations were also performed using steady flamelet and dynamic (Germano) SGS model.

Table CFD – A  Summary description of all test cases

Name	GNDPs	PDPs (Problem Definition Parameters)		SPs (Simulated Parameters)
	Re	Fuel jet composition	Pilot flame composition	Detailed data	DOAPs
CFD1 (steady flamelet with subgrid-scale model) CFD2 (CMC with subgrid-scale model) CFD3 (CMC with dynamic Germano subgrid-scale model)	22400	25% of methane (CH4) and 75% of air	C2H2, H2, air, CO2 and N2	$\langle U\rangle ,\langle U_{rms}\rangle ,\langle V\rangle ,\langle V_{rms}\rangle ,\langle \eta \rangle ,$ $\langle T\rangle ,\langle Y_{H_{2}0}\rangle ,\langle Y_{O_{2}}\rangle ,\langle Y_{N_{2}}\rangle ,\langle Y_{H_{2}}\rangle ,$ $\langle Y_{CO}\rangle ,\langle Y_{CO_{2}}\rangle ,\langle Y_{CH_{4}}\rangle ,\langle \eta _{rms}\rangle ,$ $\langle T_{rms}\rangle ,\langle Y_{{H_{2}O}_{rms}}\rangle ,\langle Y_{{O_{2}}_{rms}}\rangle ,$ $\langle Y_{{N_{2}}_{rms}}\rangle ,\langle Y_{{H_{2}}_{rms}}\rangle ,\langle Y_{{CO}_{rms}}\rangle ,$ $\langle Y_{{CO_{2}}_{rms}}\rangle ,\langle Y_{{CH_{4}}_{rms}}\rangle$	Axial profiles Tmax , x/D (Tmax ) Lconst(η , YCH4 , YO2) Lconst(YH2O , YCO2) YH2, max , z/D (YH2, max ) YCO, max , z/D (YCO, max ) RMSmax z/D (RMSmax ) Radial profiles z/D = 15, 30, 45 Fmax , Umax r½(η) , r½(U )

Table CFD – B Summary description of all available data files and simulated parameters

	SP1 $\langle U\rangle ,\langle U_{rms}\rangle ,\langle V\rangle ,\langle V_{rms}\rangle$ (ms-1) axial profiles	SP2 radial profiles x/D = 15	SP3 radial profiles x/D = 30	SP4 radial profiles x/D = 45	SP5	DOAPs
CFD1	File:AC2-09 cfd11.dat	File:AC2-09 cfd12.dat	File:AC2-09 cfd13.dat	File:AC2-09 cfd14.dat	File:AC2-09 cfd15.dat
CFD2	File:AC2-09 cfd21.dat	File:AC2-09 cfd22.dat	File:AC2-09 cfd23.dat	File:AC2-09 cfd24.dat	File:AC2-09 cfd25.dat
CFD3	File:AC2-09 cfd31.dat	File:AC2-09 cfd32.dat	File:AC2-09 cfd33.dat	File:AC2-09 cfd34.dat	File:AC2-09 cfd35.dat

Table CFD – B Summary description of all available data files and simulated parameters

	SP1	SP2	SP3	SP4	SP5
	$\langle U\rangle ,\langle U_{rms}\rangle ,\langle V\rangle ,\langle V_{rms}\rangle$ (ms-1)				$\langle F\rangle ,\langle T\rangle ,\langle Y_{H_{2}O}\rangle ,\langle Y_{O_{2}}\rangle ,$ $\langle Y_{N_{2}}\rangle ,\langle Y_{H_{2}}\rangle ,\langle Y_{CO}\rangle ,\langle Y_{CO_{2}}\rangle$ $\langle Y_{CH_{4}}\rangle ,\langle F_{rms}\rangle ,\langle T_{rms}\rangle ,\langle Y_{{H_{2}O}_{r}ms}\rangle$
	Axial profiles	Radial profiles x/D = 15	Radial profiles x/D = 30	Radial profiles x/D = 45	Axial profiles
CFD1	File:AC2-09 cfd11.dat	File:AC2-09 cfd12.dat	File:AC2-09 cfd13.dat	File:AC2-09 cfd14.dat	File:AC2-09 cfd15.dat
CFD2	File:AC2-09 cfd21.dat	File:AC2-09 cfd22.dat	File:AC2-09 cfd23.dat	File:AC2-09 cfd24.dat	File:AC2-09 cfd25.dat
CFD3	File:AC2-09 cfd31.dat	File:AC2-09 cfd32.dat	File:AC2-09 cfd33.dat	File:AC2-09 cfd34.dat	File:AC2-09 cfd35.dat

SIMULATION CASE CFD1

Solution Strategy

In the CFD1 the steady flamelet concept was applied with the standard Smagorinsky SGS model for turbulence.

Computational Domain

The CFD1 results were obtained with computational meshes 80×80×160 nodes. The computational domain at the inlet and outlet plane extended to 5.5D and 18.3D respectively in both horizontal directions. The length of the domain was equal to 50D. The mesh was stretched in axial direction by exponential function and in radial directions by hyperbolic tangent function. The grid refinement studies for the LES calculations showed that the grid resolution with 80×80×160 nodes in the proposed computational domain is sufficient and further grid refinement leads to minor changes of the statistically converged parameters. The computational domain is shown in Fig.4.

Fig. 4. Computational domain for Sandia flame D (left); mesh resolution in the inlet plane (right)

Boundary Conditions

The boundary conditions in the inlet plane were assumed to be as follows:

• the mean and RMS profiles of the axial velocity component were interpolated from experimental data - Fig.3 presents comparison of the experimental data with boundary profiles applied in computations; the random disturbances, introduced as a white noise, were scaled by RMS profile and next they were superimposed on the mean profile;
• the mean and RMS values of the radial velocity components were assumed equal to zero;
• the mixture fraction was assumed equal to 1.0 in the main jet; 0.27 in the pilot jet and zero in the co-flowing air.

At the lateral boundaries the axial velocity was assumed equal to the velocity of coflowing air (0.9 m/s) while the remaining components were equal to zero. At the outlet the convective type boundary conditions were assumed which do not require specification of any variables.

Application of Physical Models

In the most general case modeling of the combustion processes is very expensive computationally since together with the solution of the flow field it requires solution of additional transport equations for particular N species (e.g. CO, CO₂, H₂O, H₂, etc.) produced in chemical reactions. The transport equations for species have the following form:

${\frac {\partial \rho Y_{k}}{\partial t}}+{\frac {\partial \rho u_{i}Y_{k}}{\partial x_{i}}}={\frac {\partial }{\partial x_{i}}}\left(\rho D_{k}{\frac {\partial Y_{k}}{\partial x_{i}}}\right)+{\dot {\omega }}_{k}\quad k=1,2,\dots ,N$

$\left(1\right)$

where:

$\left.\rho \right.$	is the density
$\left.u_{i}\right.$	is the velocity component
$\left.Y_{k}\right.$	is the mass fraction of species $\left.k\right.$
${\dot {\omega }}_{k}$	is the reaction rate (speed of creation/destruction of a given species)
$\left.D_{k}\right.$	is the diffusion coefficient usually taken the same (denoted by $\left.D\right.$) for each species and defined as $\left.\rho D=\mu /Pr\right.$, where $\left.\mu \right.$ is the molecular viscosity and $\left.Pr\right.$ is the Prandtl number.

Reaction rate of a given species k is a sum of the reaction rates in all M reactions in which species k occurs. It is defined as:

${\dot {\omega }}_{k}=\sum _{j=1}^{M}{\dot {\omega }}_{kj}=\sum _{j=1}^{M}\left(\nu _{kj}''-\nu _{kj}'\right)\cdot W_{k}Q_{j}$

$\left(2\right)$

where:

$\left.{\dot {\omega }}_{kj}\right.$	is the reaction rate of species $\left.k\right.$ in reaction $\left.j\right.$
$\left.\nu _{kj}''\right.$ $\left.\nu _{kj}'\right.$	are the molar stochiometric coefficients after and before reaction $\left.j\right.$ respectively
$\left.W_{k}\right.$	is the atomic weight of species $\left.k\right.$
$\left.Q_{j}\right.$	is the rate of progress of reaction $\left.j\right.$; it is a function of temperature, density and species mass fraction and may be obtained from chemical kinetics or experiment; methods of determination of the rate of progress of reaction are beyond the scope of this report.

With the assumption of low Mach number flow the equation of energy may be expressed as the transport equation for the temperature which is given as :

${\frac {\partial \rho T}{\partial t}}+{\frac {\partial \rho u_{i}T}{\partial x_{i}}}={\frac {\partial }{\partial x_{i}}}\left({\frac {\lambda }{C_{p}}}{\frac {\partial T}{\partial x_{i}}}\right)+{\dot {\omega }}_{T}$

$\left(3\right)$

where:

$\left.T\right.$	is the temperature
$\left.\lambda \right.$	is the heat diffusion coefficient
$\left.C_{p}\right.$	is the specific heat
$\left.{\dot {\omega }}_{T}\right.$	is the heat release defined as a sum of product of reaction rates and formation enthalpies

The equations given above together with the Navier-Stokes equations, the equation of state and the continuity equation form a closed system which allows computing the flow field together with combustion process. However, their direct implementation in a computer code with regard for tens of species and tens (or even several hundred) of chemical reactions is still impossible from the point of view of capability of available computers, and for this reason significant simplifications have to be made.

The turbulence/combustion interaction steady flamelet concept implemented in the BOFFIN code was introduced by Peters , which stated that the flame can be seen as an ensemble of laminar flamelets. It used the equation for the conserved scalar referred to as the mixture fraction. The equation for the mixture fraction (denoted as $\left.\eta \right.$) has the simple form of convection-diffusion equation and is given as:

${\frac {\partial \rho \eta }{\partial t}}+{\frac {\partial \rho u_{i}\eta }{\partial x_{i}}}={\frac {\partial }{\partial x_{i}}}\left(\rho D{\frac {\partial \eta }{\partial x_{i}}}\right)$

$\left(4\right)$

The mixture fraction is a normalized quantity ( $\eta \in [0,1]$ ) and represents a local fuel to oxidizer ratio ( $\left.\eta =0\right.$ means pure oxidizer, $\left.\eta =1\right.$ means pure fuel ). The assumption that one conserved scalar is sufficient to describe thermochemical state of the flow decouples the modelling of reactive phenomena from that of flow modelling. Assuming that particular species and temperature are functions of the mixture fraction $\left(Y_{k}(\eta (x,t),t),T(\eta (x,t),t)\right)$ the equations $\left(1\right)$ and $\left(3\right)$ may be transformed into the mixture fraction space resulting in:

$\rho {\frac {\partial Y_{k}}{\partial t}}={\dot {\omega }}_{k}+\rho D\left({\frac {\partial \eta }{\partial x_{i}}}{\frac {\partial \eta }{\partial x_{i}}}\right){\frac {\partial ^{2}Y_{k}}{\partial \eta ^{2}}}={\dot {\omega }}_{k}+{\frac {1}{2}}\rho \chi {\frac {\partial ^{2}Y_{k}}{\partial \eta ^{2}}}$

$\left(5\right)$

$\rho {\frac {\partial T}{\partial t}}={\dot {\omega }}_{k}+\rho D\left({\frac {\partial \eta }{\partial x_{i}}}{\frac {\partial \eta }{\partial x_{i}}}\right){\frac {\partial ^{2}T}{\partial \eta ^{2}}}={\dot {\omega }}_{T}+{\frac {1}{2}}\rho \chi {\frac {\partial ^{2}T}{\partial \eta ^{2}}}$

$\left(6\right)$

In the above transformations spatial derivatives parallel to the iso-surface of mixture fraction have been neglected as they are small compared to the gradients in normal direction. In Equations $\left(5\right)$ and $\left(6\right)$ the only quantity depending on the flow field is the scalar dissipation rate $\left.\chi \right.$ reflecting the mixing process. Equations $\left(5\right)$ and $\left(6\right)$ constitute the unsteady flamelet approach in which the dependence of the temperature and species on time is retained. Assuming the structure of the flamelet to be steady, even though the mixture fraction itself depends on time, the functional dependence of the thermodynamic variables on the mixture fraction can be formulated in the form $Y_{k}(x)=Y_{k}\left(\eta (x,t)\right),\ \ T(x)=T\left(\eta (x,t)\right)$. These relations can be obtained from chemical equilibrium assumption or from laminar flamelet calculations. The latter approach is applied in the BOFFIN code in which the functional dependences $\left.Y_{k}=Y_{k}(\eta ),\ \ T=T(\eta )\right.$ and also $\rho =\rho \left(Y_{k}(\eta )\right)=\rho (\eta )$ are provided from the solution of the following system of equations:

${\dot {\omega }}_{k}=-{\frac {1}{2}}\rho \chi {\frac {\partial ^{2}Y_{k}}{\partial \eta ^{2}}}$

$\left(7\right)$

${\dot {\omega }}_{T}=-{\frac {1}{2}}\rho \chi {\frac {\partial ^{2}T}{\partial \eta ^{2}}}$

$\left(8\right)$

In the context of LES method of turbulence modelling, the mixture fraction equation has the form:

${\frac {\partial {\bar {\rho }}{\tilde {\eta }}}{\partial t}}+{\frac {\partial {\bar {\rho }}{\tilde {u_{i}}}{\tilde {\eta }}}{\partial x_{i}}}={\frac {\partial }{\partial x_{i}}}\left({\overline {\rho D{\frac {\partial \eta }{\partial x_{i}}}}}-\left({\bar {\rho }}u_{i}^{\pm }\eta -{\bar {\rho }}{\tilde {u_{i}}}{\tilde {\eta }}\right)\right)$

$\left(9\right)$

where bar and tilde represents the LES filtered and Favre-filtered variables according to the general definition:

${\bar {f}}(x)=\int f(x')G(x-x')\ dx'$

$\left(10\right)$

${\bar {\rho }}{\tilde {f}}(x)=\int \rho f(x')G(x-x')\ dx'$

$\left(11\right)$

The nonlinear interaction of the subgrid scales $(f'=f-{\bar {f}})$, in the diffusive term of Eq. $\left(9\right)$ is usually neglected as small comparing to interaction of the subgrid scales in the convective part represented by the term in brackets on the right hand side of Eq. $\left(9\right)$. In the BOFFIN code this term is modelled using the gradient hypothesis given as:

$\left({\bar {\rho }}u_{i}^{\pm }\eta -{\bar {\rho }}{\tilde {u_{i}}}{\tilde {\eta }}\right)=-{\frac {\mu _{t}}{Pr_{t}}}{\frac {\partial {\tilde {\eta }}}{\partial x_{i}}}$

$\left(12\right)$

Numerical Accuracy

The BOFFIN code is second order accurate in space and time.

References

3. Branley N., Jones W.P.: "Large eddy simulations of a turbulent non-premixed flame", Combustion and Flame, vol. 127, 2001

4. Perry R.H. and Green D.W.: Chemical Engineer's Handbook, McGraw-Hill, 1993

SIMULATION CASE CFD2

Solution Strategy

In the CFD2 case the simplified Conditional Moment Closure approach was applied. In this simplification the convective term in the physical space was neglected and the model was equivalent to the unsteady flamelet one. As a subgrid-scale model the standard Smagorinsky model was used.

Computational Domain

The computational domain for CFD2 was the same as for CFD1.

Boundary Conditions

Boundary conditions are the same as for CFD1.

Application of Physical Models

For the CMC model the simplified Smooke mechanism was applied considering 16 species: CH₄, CO, CO₂, H₂, H₂O, O₂, N₂, NO, OH, CH₃, CH₃O, CH₂O, HO₂, H₂O₂, H, HCO, O.

The CMC model independently proposed in the nineties by Bilger and Klimenko consists of solution of the conditional species concentration balance equations in the mixture fraction space. The equation proposed by Bilger and Klimenko is shown below:

${\bar {\rho }}{\frac {\partial \langle Y_{k}|\eta \rangle }{\partial t}}+\rho \langle {\tilde {\mathbf {u} }}|\eta \rangle \nabla \langle Y_{k}|\eta \rangle =\langle {\dot {\omega }}|\eta \rangle +\langle {\tilde {\chi }}|\eta \rangle {\frac {\partial ^{2}\langle Y_{k}|\rangle }{\partial \eta ^{2}}}$

$\left(1\right)$

In the present calculations the convective term with a conditional velocity was neglected and the turbulence/combustion interaction is equivalent to the unsteady flamelet model. Coupling between the balance equation of conditional species concentration corresponding to the chemical reaction and turbulent flow field is limited to the information contained in the scalar dissipation rate ${\tilde {\chi }}$. In the CMC module conditional scalar dissipation rate in the mixture fraction space is taken as the analytical solution for steady strained one-dimensional non premixed laminar flame [3] :

$\langle {\tilde {\chi }}|\eta \rangle ={\tilde {\chi }}_{0}F(\eta )={\tilde {\chi }}_{0}\exp \left\{-2\left[\operatorname {erf} ^{-1}(2\eta -1)\right]^{2}\right\}$

$\left(2\right)$

The maximum value of scalar dissipation rate is computed from the presumed pdf (β-function), which requires the information about the large scale mixture fraction (resolved scale from the LES code) and its SGS variance:

${\tilde {\chi }}_{0}={\frac {\tilde {\chi }}{\int _{0}^{1}F(\Psi )\ p(\Psi ;\mathbf {x} ,t)\ d\Psi }}$

$\left(3\right)$

The combustion model based on conditional average species concentration requires the following quantities from the turbulent flow field computed using LES:

• instantaneous values of mixture fraction
• SGS mixture fraction variance
• SGS scalar dissipation rate

In the RANS models the scalar variance is computed from its transport equation, while in LES the SGS mixture fraction is computed assuming the gradient-type approximation:

${\overline {\eta ''^{2}}}=C_{\eta }\Delta ^{2}\left({\frac {\partial {\tilde {\eta }}}{\partial x_{i}}}{\frac {\partial {\tilde {\eta }}}{\partial x_{i}}}\right)$

$\left(4\right)$

In order to close the CMC model the SGS dissipation rate must be estimated. The SGS scalar dissipation rate in BOFFIN is computed following dimensional arguments:

${\tilde {\chi }}\propto \lVert S_{ij}\rVert \ {\overline {\eta ''^{2}}}$

$\left(5\right)$

In order to solve the CMC equation in mixture fraction space the mass fractions for all chemical species must be computed in physical space, using the presumed pdf (β-function):

${\bar {\rho }}{\tilde {Y}}_{k}=\int _{0}^{1}\langle Y_{k}|\Psi \rangle \ p(\Psi ;{\tilde {\eta }},{\overline {\eta ''^{2}}})\ d\Psi$

$\left(6\right)$

Similarly, the following relation is valid for the mean temperature:

${\bar {\rho }}{\tilde {T}}=\int _{0}^{1}\langle T|\Psi \rangle \ p(\Psi ;{\tilde {\eta }},{\overline {\eta ''^{2}}})\ d\Psi$

$\left(7\right)$

Using the concentration of species at the new iteration step the density can be evaluated. The species concentrations at the new iteration step are transmitted later on to the BOFFIN code.

Solving the CMC equation in each LES cell would be prohibitively expensive. However, variation of scalar dissipation rate in space is much lower than variation of the velocity field, it was decided to solve the CMC equations only for the 16 subdomains shown in Fig. 5.

Fig. 5. Domain decomposition for parallel computations - simplified domain

Numerical Accuracy

The BOFFIN code is second order accurate in space and time.

CFD Results

CFD2 - CMC model with standard Smagorinsky SGS model for turbulence; all data files are in the Tecplot format

Simulation Case CFD3

Solution Strategy

In the CFD3 case the steady flamelet concept was applied with the dynamic SGS model for turbulence.

Computational Domain

The computational domain for CFD3 was the same as for CFD1 and CFD2.

Boundary Conditions

The boundary conditions were the same as for CFD1 and CFD2.

Application of Physical Models

The same chemical kinetics as for CFD1 and steady flamelet approach as for CFD1 were used. For turbulence, the standard Germano dynamic model was used.

Numerical Accuracy

The BOFFIN code is second order accurate in space and time.

CFD Results

CFD3 - steady flamelet model with dynamic Germano subgrid model for turbulence; all data files are in the Tecplot format.

File:AC2-09 cfd31.dat – axial and radial mean and fluctuating velocity profiles along the jet flame axis

ASCII file; columns: $x/D\ {\textrm {[-]}},\ \langle U\rangle \ {\textrm {[m/s]}},\ \langle U_{rms}\rangle \ {\textrm {[m/s]}},\ \langle V\rangle \ {\textrm {[m/s]}},\ \langle V_{rms}\rangle \ {\textrm {[m/s]}}$

File:AC2-09 cfd32.dat – axial and radial mean and fluctuating velocity profiles along radius at z/D = 15

ASCII file; columns: $R/D\ {\textrm {[-]}},\ \langle U\rangle \ {\textrm {[m/s]}},\ \langle U_{rms}\rangle \ {\textrm {[m/s]}},\ \langle V\rangle \ {\textrm {[m/s]}},\ \langle V_{rms}\rangle \ {\textrm {[m/s]}}$

File:AC2-09 cfd33.dat – axial and radial mean and fluctuating velocity profiles along radius at z/D = 30

ASCII file; columns: $R/D\ {\textrm {[-]}},\ \langle U\rangle \ {\textrm {[m/s]}},\ \langle U_{rms}\rangle \ {\textrm {[m/s]}},\ \langle V\rangle \ {\textrm {[m/s]}},\ \langle V_{rms}\rangle \ {\textrm {[m/s]}}$

File:AC2-09 cfd34.dat – axial and radial mean and fluctuating velocity profiles along radius at z/D = 45

ASCII file; columns: $R/D\ {\textrm {[-]}},\ \langle U\rangle \ {\textrm {[m/s]}},\ \langle U_{rms}\rangle \ {\textrm {[m/s]}},\ \langle V\rangle \ {\textrm {[m/s]}},\ \langle V_{rms}\rangle \ {\textrm {[m/s]}}$

---

Contributed by: Andrzej Boguslawski — Technical University of Częstochowa

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

© copyright ERCOFTAC 2024