Complex Chemistry

In Simcenter STAR-CCM+, detailed Complex Chemistry is solved using a stiff CVODE ODE (Ordinary Differential Equation) solver to integrate the chemical source terms. To consider the turbulence effects on combustion, you can also select either the Laminar Flame Concept (LFC) model or the Eddy Dissipation Concept (EDC) model.

The Complex Chemistry model is suitable for introducing detailed chemistry information to the CFD simulation. This model can solve thousands of reactions among hundreds of species—hence the term complex chemistry. Since an ODE solver is used to integrate the chemical source terms, the Complex Chemistry model can handle stiff reaction systems (reaction systems with a wide range of reaction time scales).

The model applied to integrate the chemical source term in each computational cell is the Constant Pressure Reactor.

Reaction System

The Complex Chemistry model requires detailed reaction mechanism information about species, reactions, thermodynamics, and transport properties. These details are supplied by complex chemistry definition files that are imported in the Chemkin format using the Complex Chemistry model.

Source Term Definition

The general species transport equation is formulated as follows:

Figure 1. EQUATION_DISPLAY

\frac{\partial}{\partial t} ρ Y_{i} + \frac{\partial}{\partial x_{j}} (ρ {u_{j} Y}_{i} + F_{k, j}) = ω_{i}

(3410)

where $F_{k, j}$ is the diffusion flux component and source term $ω_{i}$ is the rate of production of species $i$ .

The operator splitting algorithm takes advantage of the different time scales involved for the chemical reactions and the flow field. The time integration of the chemical state (species mass fractions $Y_{i}$ and temperature $T$ ) is performed in two steps:

At the beginning of each time-step, the chemical state is integrated in each CFD cell from state ${(Y_{i}, T, p)}^{n}$ to ${(Y_{i}, T, p)}^{*}$ , accounting only for the chemical source terms:
Figure 2. EQUATION_DISPLAY

$Y_{i}^{*} = Y_{i} + \int_{0}^{τ} r_{k} (Y, T, p) d t$

(3411)
where $Y_{i}^{*}$ is the mass fraction at the end of a time integration $τ$ with a stiff ODE solver, as calculated in Eqn. (3411). $r_{k}$ is the reaction rate from Eqn. (3358), $Y$ is the mass fraction vector, and $T$ is temperature.
The system of Eqn. (3411) is solved using the stiff CVODE solver.
The species transport equation is solved with an explicit reaction source term $ω_{i}$ for the $i$ 'th species given as:
Figure 3. EQUATION_DISPLAY

$ω_{i} = ρ f (\frac{Y_{i}^{*} - Y_{i}}{τ})$

(3412)
where $ρ$ is the density, $f$ is the mean reaction rate multiplier, $Y_{i}$ denotes the current mass fraction in the cell. The time integration $τ$ is specified for unsteady simulations, and for steady simulations is taken as the residence time $τ_{r e s}$ as defined in Eqn. (3328).

For unsteady simulations, it is essential to maintain a low Courant number to ensure that any errors in the operator splitting scheme are small.

Both the Laminar Flame Concept (LFC) and Eddy Dissipation Concept (EDC) consider the turbulence effects on combustion implicitly through the increased turbulent diffusivity that is provided by the turbulence model.

In the case of premixed flames, this increased diffusivity results in a flame thickness that is larger than the laminar flame thickness, and a turbulent flame propagation speed that is faster than the laminar flame speed. The flame thickness and speed is controlled through the turbulent Schmidt and Prandtl numbers (set these equal to each other).

In addition to this implicit turbulence effect through the turbulent diffusivity, the EDC model also considers the "turbulence-chemistry interaction" in the species chemical sources. The mean species source term in the species transport equation for both LFC and EDC is modeled as Eqn. (3412).

For the LFC model, the mean reaction rate multiplier $f$ in Eqn. (3412) is 1, and the time scale, $τ$ in Eqn. (3412) is the residence time in the cell, which is calculated as Eqn. (3328).
For the EDC model, the mean reaction rate multiplier $f$ in Eqn. (3412) is calculated as:
Figure 4. EQUATION_DISPLAY

$f = {({[C_{l} {(\frac{υ τ_{t u r b}}{L^{2}})}^{0.25}]}^{- 3} - 1)}^{- 1}$

(3413)

where $C_{l}$ is the fine structure length constant with a default value of 2.1377, $υ$ is the kinematic viscosity, $τ_{t u r b}$ is the turbulent time-scale, and $L$ is the turbulent length scale.

For the EDC model, the time scale $τ$ in Eqn. (3412), is calculated as:
Figure 5. EQUATION_DISPLAY

$τ = τ_{t u r b}$

(3414)

where $τ_{t u r b}$ is the turbulent time-scale which is calculated as a constant with the default 0.4082, times the Kolmogorov turbulent time-scale, defined as $\sqrt{\frac{ν}{ε}}$ .

The equations Eqn. (3411) and Eqn. (3412) represent the average rate of change of the species in the cell over the time-step $τ$ . The LFC time-step is approximately the residence time in the cell Eqn. (3328). The EDC time-step in Eqn. (3414) is close to the Kolmogorov time-scale $τ_{k}$ (as described in the section, Turbulence Transfer between Phases), which is a model for the time-scale of the smallest turbulent eddy. Since the Kolmogorov time is typically less than the residence time, and the EDC mean reaction rate multiplier $f$ in Eqn. (3413) is less than 1, EDC typically models the effect of turbulence as a reduction in the mean reaction rate compared to LFC.
The Relax to Chemical Equilibrium method assumes that the chemical composition relaxes to the local equilibrium composition at a time scale that is determined by flow and chemistry time scales. The equilibrium composition can optionally be computed using ISAT.
The chemistry source term is computed based on the specified fuel species:
Figure 6. EQUATION_DISPLAY

$ω_{i} = ρ \frac{Y_{i} - Y_{i}^{e q}}{τ_{c h a r}}$

(3415)
where $Y_{i}$ is the mass fraction for species $i$ , $Y_{i}^{e q}$ is the local instantaneous thermodynamic equilibrium mass fraction for species $i$ , and $τ_{c h a r}$ is a characteristic time scale:
- When using models other than the Turbulent Flame speed Closure (TFC) model:
  Figure 7. EQUATION_DISPLAY
  
  $τ_{c h a r} = τ_{f l o w} + τ_{c h e m}$
  
  (3416)
  
  where:
  - $τ_{f l o w} = k / (ε c_{f l o w})$ , where $ε$ is the turbulent dissipation rate, $k$ is the turbulent kinetic energy, and $c_{f l o w}$ is a turbulent rate constant with a default value of 4.
    For laminar flow, $τ_{f l o w} = Δ / ({U c}_{f l o w})$ , where $Δ$ is the grid size and $U$ is the local velocity.
  - $τ_{c h e m} = c_{c h e m} \min_{f} (\frac{Y_{f}}{r_{f}})$ , where $c_{c h e m}$ is a model constant with a default value of 2, $f$ is an index for a user-specified list of fuel species, $\min_{f}$ is taking the minimum of all the user-specified fuel species, $Y_{f}$ is the fuel species mass fraction in a CFD cell, and $r_{f}$ is the net reaction rate of species $f$ .
- When using the TFC model:
  Figure 8. EQUATION_DISPLAY
  
  $τ_{c h a r} = c_{t} τ_{k}$
  
  (3417)
  where $c_{t}$ is the timescale constant and $τ_{t}$ is the Kolmogorov turbulent time-scale, defined as $\sqrt{\frac{ν}{ε}}$ .

Solution Procedure

In Simcenter STAR-CCM+, the reacting species transport equations are solved using the CVODE solver with the operator splitting algorithm to find an average reaction rate to remove stiffness. When running with the Steady State model, an artificial chemistry time step is introduced which is based on convection and diffusion fluxes in that cell.

The time step for integration of the chemistry in a steady state simulation $τ$ for each cell $p$ is determined by:

Figure 9. EQUATION_DISPLAY

τ = \frac{V_{p} ρ_{p}}{A_{p}}

(3418)

Figure 10. EQUATION_DISPLAY

= \frac{m}{\dot{m}}

(3419)

Figure 11. EQUATION_DISPLAY

= τ_{r e s}

(3420)

which is a function of the cell volume $V_{p}$ , density $ρ_{p}$ , and central coefficient $A_{p}$ which results from discretization and linearization. Eqn. (3418) is approximately the residence time Eqn. (3328) in the cell.

The Complex Chemistry model uses either the CVODE solver to integrate the stiff chemistry over a time-step.

The equations that are integrated in the CVODE solver correspond to a constant pressure reactor:

Figure 12. EQUATION_DISPLAY

\frac{d Y_{k}}{d t} = \frac{{\dot{ω}}_{k}}{ρ}

(3421)

Figure 13. EQUATION_DISPLAY

\frac{d T}{d t} = - \frac{Σ_{k = 1}^{n} {\dot{ω}}_{k} h_{k}}{ρ C_{p}}

(3422)

where

t

is time,

Y_{k}

is the mass fraction,

T

is temperature,

{\dot{ω}}_{k}

is the net reaction rate of species

k

, given by Eqn. (3362),

n

is the number of species,

ρ

is the density,

h_{k}

is the enthalpy, and

C_{p}

is the specific heat.

The evolution of mean species mass fractions and temperature in a given CFD cell is calculated by integrating at the cell conditions using Jacobian matrices to integrate the ordinary differential equation:

Figure 14. EQUATION_DISPLAY

Y_{i} = Y_{i} (t_{0}) + \int_{0}^{τ_{r e s}} {\dot{ω}}_{i (t)} d t

(3423)

The Jacobian

J_{ω}

is defined as:

Figure 16. EQUATION_DISPLAY

J_{ω} = \frac{\partial ω}{\partial Y}

(3425)

For

n

species, the Jacobian matrix is written as:

Figure 17. EQUATION_DISPLAY

J_{ω} = [\begin{matrix} \frac{\begin{matrix} \partial ω_{1} \end{matrix}}{\partial Y_{1}} & \frac{\begin{matrix} \partial ω_{1} \end{matrix}}{\partial Y_{2}} & \dots & \frac{\begin{matrix} \partial ω_{1} \end{matrix}}{\partial Y_{n}} \\ \frac{\begin{matrix} \partial ω_{2} \end{matrix}}{\partial Y_{1}} & \frac{\begin{matrix} \partial ω_{2} \end{matrix}}{\partial Y_{2}} & \dots & \frac{\begin{matrix} \partial ω_{2} \end{matrix}}{\partial Y_{n}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{\begin{matrix} \partial ω_{n} \end{matrix}}{\partial Y_{1}} & \frac{\begin{matrix} \partial ω_{n} \end{matrix}}{\partial Y_{2}} & \dots & \frac{\begin{matrix} \partial ω_{n} \end{matrix}}{\partial Y_{n}} \end{matrix}]

(3426)

A numerical Jacobian describes a Jacobian that is derived using the finite differences method for each species in the system. Since the chemical source term, $ω$ , is calculated repeatedly as the CVODE solver iterates, the cost of solving the ODE is expensive.

An analytical Jacobian is calculated with the chain rule to obtain an analytical expression for $J_{ω}$ .

Since this process is computationally expensive, ISAT, Clustering, or Dynamic Mechanism Reduction are used to reduce the expense and reach a converged solution faster.

ISAT

In-Situ Adaptive Tabulation (ISAT) dynamically tabulates functions that are computationally expensive to evaluate, and interpolates many times after tabulation [766]. The ISAT method approximates the solution of the ODE (Eqn. (3411)) for a given initial condition through:

Thermodynamic variable (temperature or total enthalpy).
Composition variable (mass fraction or molar number).
Possibly pressure and time step (depending on the type of the simulation).

These values constitute an independent variable or a query point (

ϕ

). The solution

f

of the ODE equation (Eqn. (3411)) is called reaction mapping.

The ISAT method starts with an empty table and gradually adds points to it. To solve the first ODE equation for

ϕ_{0}

(a first query), ISAT does the following:

It calls the ODE solver to calculate $f_{0}$ .
It creates an entry and adds it to the table.

Every table entry consists of:

An independent variable $ϕ$
A reaction mapping variable $f$
A mapping gradient matrix $A = \frac{\partial f}{\partial ϕ}$
An Ellipsoid Of Accuracy (EOA), derived by a singular value decomposition of the mapping gradient.

For all subsequent queries, ISAT does the following for query point

ϕ_{n}

It traverses the table and finds a suitable point to check for the possible approximation.
For the suitable point $ϕ_{t}$ it checks the value: $f_{q} = f_{t} + A_{t} (ϕ_{n} - ϕ_{t})$

The following two scenarios can happen:

Either $f_{q}$ lies in the EOA, which means that the value $f_{q}$ approximates, within the user-specified tolerance, the true value of $f (ϕ_{n}))$ . In other words, query $ϕ_{n}$ is retrieved.
Or $f_{q}$ lies outside of the EOA. In this case, the true value of $f (ϕ_{n}))$ is calculated and the difference between the two values is deduced. If this difference is smaller than the user-defined tolerance, this means that the EOA at point $ϕ_{n}$ is too strict. In this case, the ellipsoid is grown to accept point $f_{q}$ . If the difference is higher than the specified tolerance, the point is added to the table.

Points in the table are stored in the leaves of the binary space partition tree. This tree divides the region in the

n

-dimensional space (

n

being the dimension of the query point) currently occupied by the points into smaller subregions for faster search. When a new query point is added to the tree, an original leaf becomes now two leafs, split by a hyperplane. One new leaf contains the old point, which lies on one side of the hyperplane. The second new leaf, which contains the new query point, lies on the other side of the hyperplane.

Therefore the following properties characterize the ISAT method:

Number of retrieves: the number of queries for which the reaction mapping was approximated
Number of grows: the number of queries for which the reaction mapping was calculated but the point was not added to the table because the EOA of the old point in the table was grown
Number of adds: the number of queries for which the reaction mapping was calculated and the point was added to the table

Growing of the EOAs can be deactivated through a property within the ISAT approximation method. Since EOA is only an approximation of the actual region of accuracy, growing of the EOA can introduce errors in the approximation table. Furthermore, it is possible that adding a point to the table will cover more future queries than growing the existing EOA.

Clustering

Clustering reduces the computational expense of complex chemistry calculations by averaging together cells with similar chemical compositions, integrating the reduced ODE set, and then interpolating the clusters back to the cells.

A low-dimensional, uniformly spaced clustering grid is formed from a few composition variables. By default, Simcenter STAR-CCM+ uses two clustering variables for unsteady simulations—equivalence ratio and temperature. For steady-state simulations, the log of the chemical time scale is also used as an additional clustering variable.

The clustering tolerance determines the distance of the clustering composition variable between the clustering grid points. For relative tolerances

ε_{r e l}

the number of grid points for the clustering composition variable

ϕ

is:

Figure 21. EQUATION_DISPLAY

\frac{1}{ε_{r e l}}

(3430)

For absolute tolerances

ε_{a b s}

, the number of grid points for the clustering composition variable

ϕ

is:

Figure 22. EQUATION_DISPLAY

\frac{ϕ_{\max} - ϕ_{\min}}{ε_{a b s}}

(3431)

where

ϕ_{\max}

and

ϕ_{\min}

represent the maximum and minimum values of

ϕ

in the computational domain.

Following Babajimopolous [750], the clustering equivalence ratio is defined as:

Figure 23. EQUATION_DISPLAY

ϕ = \frac{2 C^{*} + \frac{H^{*}}{2}}{O^{*}}

(3432)

where:

$C^{*}$ is the number of $C$ atoms in all species except ${C O}_{2}$ , multiplied by the species mass fraction and divided by the species molecular weight.
$H^{*}$ is the number of $H$ atoms in all species except $H_{2} O$ , multiplied by the species mass fraction and divided by the species molecular weight.
$O^{*}$ is the number of $O$ atoms in all species except ${C O}_{2}$ and $H_{2} O$ , multiplied by the species mass fraction and divided by the species molecular weight.

In addition to the default clustering variables, you can select additional clustering variables such as mixture fraction and species. The clustering mixture fraction is defined as the elemental mass fraction of Hydrogen and Oxygen atoms in the mixture.

Within each clustering grid point, species mass fractions, energy, pressure—and for steady-state, the chemistry time-step—are cell mass-weighted averaged to form the cluster mean. This composition (

Y_{i}^{*}

in Eqn. (3411)) is integrated with the CVODE solver, as in Eqn. (3411). The source term is then applied to all cells in the cluster.

Dynamic Mechanism Reduction

Dynamic Mechanism Reduction allows the CVODE solver in Simcenter STAR-CCM+ to solve for a reduced number of species in a full chemical mechanism. The Directed Relation Graph (DRG) algorithm [755] considers the reaction rates of all species in a cell at each time-step or iteration and removes species from the mechanism if they contribute insignificantly (based on the error tolerance) to production/destruction of any other species over the time-step.

When deciding if a species can be removed from the full chemical mechanism, the DRG algorithm considers how significant an influence that species has on the rate of production of other species.

For any species A, the effect of removing species B on the rate of production of species A is determined as follows:

Figure 24. EQUATION_DISPLAY

r_{A B} = \frac{\max_{i} | ν_{A, i} ω_{i} δ_{B, i} |}{\max_{i} | ν_{A, i} ω_{i} |}

(3433)

if the $i$ th elementary reaction involves species B, $δ_{B, i} = 1$
if the $i$ th elementary reaction does not involve species B, $δ_{B, i} = 0$

When the normalized contribution of species B to the production rate of species A (

r_{A B}

) is greater than the error tolerance

ε

that is specified, the DRG algorithm also retains species B. Species A and B are said to be strongly coupled and removal of species B would impact the production rate of species A significantly.

Each node (A to F in the diagram below) in a DRG represents a species in the full chemical mechanism.

All arrows in the diagram above represent directed edges—these are edges in the DRG that exist between species with a dependence. The arrows show the direction of dependence of one species on another—the width indicates the strength of the dependence. Since A is dependant on B and C, and C is dependant on F, the group of species A, B, C, and F must be maintained to correctly calculate A. The remaining species can be removed as they do not have a significant effect on the production rate of A.