Bachelorthesis

Karlsruhe Institute of Technology
Department of Chemical and Process Engineering

Engler-Bunte-Institute
Division of Combustion Technology
Prof. Dr.-Ing. Henning Bockhorn
Bachelor Thesis
Implementation and Validation of a Solver

for Direct Numerical Simulations of
Turbulent Reacting Flows in OpenFOAM
by
cand. chem. ing. Henning Bonart
Evaluator: Prof. Dr.-Ing. Henning Bockhorn

Supervisor: Dipl.-Ing. Feichi Zhang
October 2012
Karlsruher Institut fr Technologie

Fakultt fr Chemieingenieurwesen und Verfahrenstechnik
Engler-Bunte-Institut
Bereich Verbrennunsgtechnologie
Aufgabenstellung zur Bachelorarbeit

cand. chem. ing. Henning Bonart
Implementierung und Validierung eines Lsers fr Direkte

Numerische Simulationen von turbulenten, reagierenden
Strmungen in OpenFOAM
Ziel dieser Arbeit ist die Implementierung eines Lsers fr die dreidimensionale Direkte
Numerische Simulation (DNS) von turbulenten, reagierenden Strmungen in der Entwicklungsumgebung OpenFOAM. Dabei soll die Berechnung der molekularen Flsse auf der
kinetischen Gastheorie basieren. Darber hinaus sollen verschiedene komplexe chemische
Reaktionsmechanismen eingebunden werden knnen.
Mittels der Simulation einer eindimensionalen, laminaren Vormischflamme soll der
Lser mit Ergebnissen von CHEMKIN/PREMIX im Hinblick auf Diffusions- und Reaktionsvorgnge validiert werden. Danach soll gezeigt werden, das dreidimensionale DNS
von turbulenten, reagierenden Strmungen mit dem Lser mglich sind. Zur Durchfhrung der Simulationen soll der hauseigene Rechencluster sowie der CRAY XE6 vom
Hochleistungsrechenzentrum Stuttgart genutzt werden.
Die Arbeit gliedert sich in folgende Arbeitsschritte:
1. Einarbeitung in die Theorie der turbulenten Strmung und Verbrennung
2. Aussuchen eines geeigneten Lsers in OpenFOAM als Basisstrmungslser
3. Erweiterung der Transportgleichungen des Basislsers auf multikomponenten Reaktionsmischungen
4. Einbindung von Cantera zur Berechnung der Eigenschaften von multikomponenten Reaktionsmischungen wie Dichte, spezifische Wrmekapazitten, Wrmeleitkoeffizient, Diffusionskoeffizienten, chemische Quellterme etc.
5. Validierung des neuen Lsers anhand einer eindimensionalen laminaren Vormischflamme
6. Durchfhrung einer dreidimensionalen DNS einer turbulenten Vormischflamme
7. Beschreibung von Theorie und Implementierung sowie Darstellung und Diskussion
der Ergebnisse in einer schriftlichen Ausarbeitung
8. Vorstellung der Bachelorarbeit in einem ffentlichen Seminarvortrag
Betreuer: Dipl.-Ing. Feichi Zhang
Aufgabensteller: Prof. Dr.-Ing. Henning Bockhorn
I declare that I have developed and written the enclosed thesis completely by myself,
and have not used sources or means without declaration in the text.
Karlsruhe, October 2012 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

(Henning Bonart)
Contents
Figures, Tables and Listings
III
Nomenclature
1. Introduction
1.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2. Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3. Structure of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1
3
3
2. Basic Aspects of Turbulent Combustion

2.1. Structure of Laminar Premixed Flames . . . . . . . . . . .
2.2. Characteristics of Turbulent Flow . . . . . . . . . . . . . .
2.2.1. Energy Cascade and Kolmogorov Scales . . . . . .
2.3. Interaction between Premixed Flames and Turbulent Flow
2.3.1. Turbulent Flame Speed . . . . . . . . . . . . . . . .
2.3.2. Flame Stretch Rate . . . . . . . . . . . . . . . . . .
2.3.3. Combustion Regimes . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
7
7
8
9
10
11
3. Mathematical Description of Chemically Reacting Flows

3.1. Statistical Thermodynamics and the Rigorous Kinetic Theory
3.1.1. The Boltzmann Equation . . . . . . . . . . . . . . . .
3.1.2. Enskogs General Transport Equation . . . . . . . . . .
3.1.3. Chapman-Enskog Theory . . . . . . . . . . . . . . . .
3.2. Transport Equations for Chemically Reacting Flows . . . . . .
3.2.1. Equations of Mass . . . . . . . . . . . . . . . . . . . .
3.2.2. Equation of Momentum . . . . . . . . . . . . . . . . .
3.2.3. Equation of Energy . . . . . . . . . . . . . . . . . . . .
3.3. Multicomponent Molecular Transport . . . . . . . . . . . . . .
3.3.1. Species Mass Flux . . . . . . . . . . . . . . . . . . . .
3.3.2. Momentum Flux . . . . . . . . . . . . . . . . . . . . .
3.3.3. Heat Flux . . . . . . . . . . . . . . . . . . . . . . . . .
3.4. Reaction Kinetics . . . . . . . . . . . . . . . . . . . . . . . . .
3.5. Chemical Time Scales . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
13
13
14
15
15
16
16
17
17
18
18
20
21
21
22
.
.
.
.
.
.
.
4. Numerical Solution of Partial Differential Equations

23
4.1. The Finite Volume Method . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.1.1. Approximation of Integrals . . . . . . . . . . . . . . . . . . . . . . 24
4.1.2. Interpolation and Differentiation Procedures . . . . . . . . . . . .

4.2. Methods for Unsteady Problems . . . . . . . . . . . . . . . . . . . . . . .
4.3. Solution of Linear Equation Systems . . . . . . . . . . . . . . . . . . . .
25
26
27
5. Implemented Solver in OpenFOAM and the Cantera Interface

5.1. Software Packages . . . . . . . . . . . . . . . . . . . . . . . . .
5.2. Connection of OpenFOAM with Cantera . . . . . . . . . . . .
5.2.1. Structure of the Interface . . . . . . . . . . . . . . . . .
5.2.2. Integration of the Interface into OpenFOAM . . . . . .
5.2.3. Exchange of Data between OpenFOAM and Cantera .
5.3. Implementation of the Solver . . . . . . . . . . . . . . . . . . .
5.3.1. OpenFOAMs Approach for Transport Equations . . .
5.3.2. Implementation of Species Mass Equations . . . . . . .
5.3.3. Implementation of Sensible Enthalpy Equation . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
29
29
30
31
31
32
34
34
35
36
6. Validation of the Solver

6.1. Case Description . . . . . . . . . . . . . . . . .
6.2. Numerical conditions . . . . . . . . . . . . . . .
6.3. Comparison with CHEMKIN/PREMIX . . . . .
6.3.1. Temperature and Species Mass Fractions
6.4. Influence of the Grid Resolution . . . . . . . . .
6.5. Species Production Rates . . . . . . . . . . . . .
6.6. Computation of Chemical Time Scales . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
37
37
38
40
40
42
44
44
7. Direct Numerical Simulation of a Turbulent Premixed Flame in 3D

7.1. Numerical Setup . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2. Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3. Topological Results . . . . . . . . . . . . . . . . . . . . . . . . . .
7.4. Parallel Performance of the Cray XE6 (HERMIT) . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
46
46
48
49
52
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
8. Summary and Perspective
54
A. Code Fragments for OpenFOAM in C++3
55
B. Reaction Mechanism
57
II
Figures
1.1. World energy-related CO2 emissions by scenario . . . . . . . . . . . . . .
2.1. Species and temperature profiles for a laminar, premixed flat methaneoxygen flame. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2. Schematic representation of the turbulent kinetic energy spectrum E as
a function of the wavenumber k. . . . . . . . . . . . . . . . . . . . . . . .
2.3. Kinematic interaction between a turbulent eddy and a propagating flame
front. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4. Schematic representation of the turbulent flame velocity. . . . . . . . . .
2.5. Schematic classification of turbulent combustion regimes. . . . . . . . . .
9
10
11
3.1. Graphical expression of the Boltzmann equation. . . . . . . . . . . . . . .
14
4.1. Parameters used in the Finite Volume Method. . . .

4.2. Example of a 2D, structured, non-orthogonal grid to
through a duct. . . . . . . . . . . . . . . . . . . . . .
4.3. Approximation of the time integral over an interval. .
24
. . . . . . . . . . .
simulate the flow
. . . . . . . . . . .
. . . . . . . . . . .
5.1. Inheritance and dependency diagram for the coupling library. .

5.2. Initialization sequence of the coupling library during runtime.
5.3. Data exchange between OpenFOAM and Cantera through the
canteraFoamModel. . . . . . . . . . . . . . . . . . . . . . . . .
6
8
25
27
. . . . . .
. . . . . .
interface
. . . . . .
32
33
6.1. Physical domain of the laminar, premixed flame. . . . . . . . . . . . . . .

6.2. Numerical domain of the laminar, premixed flame. . . . . . . . . . . . . .
6.3. Temperature and species mass fraction profiles obtained from OpenFOAM
and from CHEMKIN/PREMIX. . . . . . . . . . . . . . . . . . . . . . . .
6.4. Profiles obtained from OpenFOAM with 600 Cells and 3000 Cells. . . . .
6.5. Reaction rates obtained from OpenFOAM. . . . . . . . . . . . . . . . . .
6.6. Calculated chemical time scales for T , O2 and OH. . . . . . . . . . . .
38
39
7.1.
7.2.
7.3.
7.4.
47
48
48
Numerical domain of the turbulent, premixed flame. . . . . . . . . . . . .

Initial velocity and temperature fields of the premixed turbulent flame. .
From the initial conditions calculated Kolmogorov scales. . . . . . . . . .
Isosurface of YCH4 . Bottom with eddy dissipation rate and backside with
vorticity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.5. Two dimensional slices of the temperature field and the vorticity field with
heat release rate of the turbulent premixed flame. . . . . . . . . . . . . .
III
34
41
43
44
45
49
50
7.6. Different mass fractions and reaction rates. . . . . . . . . . . . . . . . . .

7.7. Scale-up on Cray XE6 with a grid of 2 106 cells. . . . . . . . . . . . . .
51
53
Tables
6.1.
6.2.
6.3.
6.4.
Physical conditions of the simulation with OpenFOAM. . . . . . . . . . .

Numerical setups. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
In OpenFOAM used discretization schemes for the mathematical terms. .
Comparison of distinctive results obtained with OpenFOAM and CHEMKIN/PREMIX. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.5. Comparison of distinctive results obtained with 600 Cells and 3000 Cells.
37
38
39
40
42
7.1. Physical and numerical conditions of the turbulent, premixed flame in 3D. 46
Listings
A.1. Implementation of the species equation 5.8 in OpenFOAM. . . . . . . . .
A.2. Implementation of the energy equation 5.12 in OpenFOAM. . . . . . . .
A.3. Initialization during run time and performing a downcast to obtain access
to derived class functions in OpenFOAM. . . . . . . . . . . . . . . . . . .
55
56
B.1. The used reaction mechanism in CHEMKIN format . . . . . . . . . . . .
57
IV
56
Nomenclature
Latin symbols
h0f . . . . . . . . . . . .
D ..............
M .............
M ..............
A ..............
A? . . . . . . . . . . . . .
cp . . . . . . . . . . . . . .
Co . . . . . . . . . . . . .
D ..............
d ...............
d ...............
DkT . . . . . . . . . . . . .
Da . . . . . . . . . . . . .
E ..............
Ea . . . . . . . . . . . . .
F ..............
f ...............
f ...............
g ...............
H ..............
h ...............
hc . . . . . . . . . . . . . .
hs . . . . . . . . . . . . . .
I ...............
j ...............
K ..............
k ...............
k ...............
Ka . . . . . . . . . . . . .
L ...............
l ...............
Le . . . . . . . . . . . . .
m ..............
N ..............
n ...............
p ...............
Specific standard-state heat of formation

Binary diffusion coefficient
Name
Molecular weight
Surface
Pre-exponential constant
Specific heat at constant pressure
Courant number
Ordinary multicomponent diffusion coefficient
Distance
Molecular driving force
Thermal diffusion coefficient
Damkhler number
Energy
Activation energy
External force
Cell face
Velocity distribution function
Gravitational force
Terms of higher order
Specific enthalpy
Specific chemical enthaly
Specific sensible enthaly
Total number of chemical reactions
Diffusive mass flux
Total number of chemical species
Rate constant
Wavenumber
Karlovitz number
Matrix for the computation of transport quantities
Length
Lewis number
Mass
Number of molecules
Number of moles
Pressure
[J/kg]
[m2 /s]
[-]
[kg/kmol]
[m2 ]
[varies]
[J/(kgK)]
[-]
2
[m /s]
[m]
[1/m]
[kg/(ms)]
[-]
[J]
[J/mol]
[kgm/s2 ]
[-]
[s3 /m6 ]
[kgm/s2 ]
[-]
[J/kg]
[J/kg]
[J/kg]
[-]
3
[kg/(m s)]
[-]
[varies]
[1/m]
[-]
[-]
[m]
[-]
[kg]
[-]
[-]
[Pa]
Q ..............
q ...............
r ...............
Rs . . . . . . . . . . . . .
Re . . . . . . . . . . . . .
S ...............
s ...............
Sc . . . . . . . . . . . . . .
T ..............
t ...............
TO . . . . . . . . . . . . .
V ..............
V ..............
v ...............
w ..............
X ..............
x, y, z . . . . . . . . . .
Y ..............
Source
Heat flux
Chemical reaction rate
Specific gas constant
Reynolds number
Surface
Flame speed
Schmidt number
Temperature
Time
Inner reaction layer temperature
Mass diffusion velocity
Volume
Velocity
Convective or diffusive flux
Mole fraction
Direction
Mass fraction
Greek symbols
...............
...............
...............
...............
...............
...............
...............
0 . . . . . . . . . . . . . .
...............
...............
...............
...............
...............
...............
...............
...............
...............
...............
..............
...............
Temperature exponent
Molecular collision
Thickness of inner reaction layer
Rate of dissipation of turbulent kinetic energy
Perturbation function
Transport coefficient
Air equivalence number
Thermal conductivity
Dynamic viscosity
Kinematic viscosity
Stoichiometric coefficient
Ordering parameter
Mass density
Collision diameter
Stress tensor
Time scale
Physical quantity
Quantity
Reduced collision integral
Rate of progress
Indices and Accents

0 . . . . . . . . . . . . . . . Unburnt
VI
[-]
[J/(m s)]
[kg/m3 ]
[J/(kgK)]
[-]
[m2 ]
[m/s]
[-]
[K]
[s]
[K]
[m/s]
[m3 ]
[m/s]
[-]
[-]
[m]
[-]
2
[-]
[-]
[m]
[m2 /s3 ]
[-]
[-]
[-]
[W/(mK)]
[kg/(ms)]
[m2 /s]
[-]
[-]
[kg/m3 ]
[m]
[kg/(ms2 )]
[s]
[-]
[-]
[-]
[kg/(m3 s)]
...............
...............
...............
..............
~ ...............
c ...............
i ...............
k ...............
L ...............
T ..............
Mean value
Transposed
Time derivation/rate
Surrounding
Vector
Characteristic
Species i
Species k
Laminar
Turbulent
Constants . . . . .
kB . . . . . . . . . . . . . Boltzmann constant (1.38064881023 )
NA . . . . . . . . . . . . . Avogadro constant (6.022141291023 )
R . . . . . . . . . . . . . . Universal gas constant (8.3144621)
Abbreviations
1D . . . . . . . . . . . . .
3D . . . . . . . . . . . . .
fvc . . . . . . . . . . . . .
fvm . . . . . . . . . . . .
CFD . . . . . . . . . . .
BLAS . . . . . . . . . .
CFL . . . . . . . . . . . .
CG . . . . . . . . . . . . .
CHEMKIN . . . . .
DNS . . . . . . . . . . .
FDM . . . . . . . . . . .
FVM . . . . . . . . . . .
GPL . . . . . . . . . . .
LAPACK . . . . . . .
LES . . . . . . . . . . . .
LU . . . . . . . . . . . . .
MATLAB . . . . . .
MPI . . . . . . . . . . . .
MPI . . . . . . . . . . . .
NRBC . . . . . . . . .
ODE . . . . . . . . . . .
OpenFOAM . . . .
[J/K]
[mol1 ]
[J/(molK)]
One dimension
Three dimensions
finite volume calculus
finite volume method
Computational Fluid Dynamics
Basic linear algebra subprograms
CourantFriedrichsLewy condition
Conjugate gradients
Software for chemical kinetics
Direct numerical simulation
Finite Difference Method
Finite Volume Method
General Public License
Linear algebra package
Large eddy simulation
Lower upper
Matrix laboratory
Message Passing Interface
Message passing interface, e.g. MPICH2/OpenMPI
Non-reflecting boundary conditions
Ordinary differential equation
Open Source Field Operation and Manipulation, registered trademark
of OpenCFD Limited
PBiCG . . . . . . . . . Preconditioned biconjugated gradient
PCG . . . . . . . . . . . Preconditioned conjugate gradient
PDE . . . . . . . . . . . Partial differential equation
VII
PISO . . . . . . . . . . .
RANS . . . . . . . . . .
SIMPLE . . . . . . .
SUNDIALS . . . .
Pressure implicit with split operator

Reynolds-averaged NavierStokes
Semi-implicit method for pressure-linked equations
Suite of nonlinear and differential/algebraic equation solvers
VIII
1. Introduction
1.1. Motivation
Chemically reacting flows impact many aspects of human life via, for example, medicinal
chemistry, chemical synthesis, material processing and combustion processes [1]. Combustion of hydrocarbon fossil fuels (or fuels stemming from renewable sources) is still an
important part of almost every energy conversion [2].
Technical applications based on combustion processes are important for transportation,
power generation and chemical engineering. At the same time pollution caused by
combustion processes leads to environmental problems. For instance, the contribution
of CO2 emissions to global warming is now considered to be proven and may be a
main threat to our standard of living. Figure 1.1 shows the world-energy related CO2
emissions by different scenarios. To reach the 450 Scenario1 an improvement of direct
and indirect energy conversion efficiency is without alternative [2].
Fig. 1.1.: World energy-related CO2 emissions by scenario. Figure taken from [2].
Other pollutants, such as unburnt hydrocarbons, soot, nitrogen oxides and sulfur oxides, can contribute to health hazards and cause smogs or acid rain [3]. However, due to
stricter government regulations for, e.g., the automotive industry and the energy sector,
pollutants emission has been reduced. In order to be conform to such legal regulations
1
The 450 Scenario describes pathways to a long-range CO2 concentration in the atmosphere of 450
parts per million, see [2].
1.1. MOTIVATION
and to meet the public expectation about cleaner and more efficiently industrial processes as well as to reach the 450 Scenario target of the Internal Energy Agency, it is
still of crucial importance to improve combustion technology.
Since most of the industrial applications of combustion involve turbulent flows, a deeper
understanding of the interplay between turbulence and combustion is needed. In many
cases turbulence increases combustion. On the other hand the heat release leads to gas
expansion and density variations and therefore influences the turbulent flow [4]. This
interaction of turbulence and chemical reaction is still not fully understood. Due to the
small length and time scales, details of the mixing process, as well as the structures in
turbulent flows, have not been fully investigated, yet. Therefore, a deeper understanding
of the interplay between turbulence and chemical reaction processes is of large interest
in order to improve efficiency in combustion.
The mechanisms involved in such turbulent combustion flows have been the object of
numerous theoretical, experimental and numerical scientific research in the last century [46]. Thereby, experimental studies have greatly advanced our knowledge on turbulent combustion. Recently, experiments that simultaneously measured the turbulent
velocity fields and the reacting scalar concentrations, together with the local temperature, have contributed to a more realistic understanding of such turbulent burning
processes. However, in order to gain detailed insight into the fundamental physics of the
turbulence-chemistry interaction on small space and time scales, numerical simulations
are to date without alternative.
In numerical simulations of turbulent combustion, additionally to the standard complexities of turbulent non-reacting simulation, other difficulties, such as the strong heat
release and the complex chemistry behind the combustion process, have to be taken
into account [7]. The computational grid has to be sufficient small in order to resolve
the smallest eddies in the turbulent flow. Due to numerical stability reasons very small
time steps are necessary. Hence, direct numerical simulations of turbulent reacting flows
demand a prohibitive amount of computational time. However, the huge progress in
computer technology in the past two decades now offers the opportunity for direct numerical simulations of three-dimensional (3D) turbulent reacting flows for certain cases.
This allows for studying the chemistry-turbulence interaction in more detail (see e.g. [8,
9]).
For a direct numerical simulation the underlying physical and chemical models which
describe chemical kinetics or molecular transport are required to be as precise as possible.
Thanks to extensive research in chemical kinetics a lot of chemical mechanisms can now
be modeled with great precision. Consequently, highly realistic chemical kinetic models
are implemented in many fluid dynamics codes. However, the use of multicomponent
transport equations based on the Maxwell-Stefan diffusion model or the rigorous kinetic
theory of ideal gases is still rare in computational fluid dynamics. The reason is probably
the additional complexity of the implementation and the higher computational cost.
Numerous studies have shown that, when omitting multicomponent diffusion equations,
significant errors occur in simulations of laminar flames [1, 1013]. As to turbulent flow,
1.2. OBJECTIVES
scientific publications comparing different diffusion models are very rare but it has been
argued by [9, 14] that accuracy should be improved when using the full multicomponent
transport equations.
1.2. Objectives
In the present thesis the main objective is to implement and validate a DNS solver which
is able to simulate turbulent, chemically reacting flows with a detailed diffusion modeling
and complex chemical kinetics. The equations are implemented with as little simplifications as possible and mostly derived from the rigorous kinetic theory of gases. The C++
toolbox OpenFOAM is used as the underlying structure for the DNS solver. The main
advantages of OpenFOAM are its high numerical accuracy and its free availability [15].
Moreover, multicomponent transport coefficients, which are obtained from the extended
Chapman-Enskog theory [16] and based on the molecular properties of the gases, are
used. These multicomponent diffusion coefficients are calculated with Cantera, an open
source chemical kinetics software. Cantera features a high calculation speed when calculating chemical source terms based on complex chemical mechanisms as well as efficient
calculation routines for multicomponent coefficients [17]. The interconnection between
OpenFOAM 2.0.1 and Cantera is newly programmed based on an library by [18] for
OpenFOAM 1.5.
To validate the solver laminar premixed flames are simulated and compared to calculations obtained with CHEMKIN/PREMIX with respect to numerical accuracy and
stability. Furthermore, the simulation of a turbulent premixed flame in 3D is carried
out. For that purpose, the solver, as well as the OpenFOAM framework and the Cantera
routines, are adapted for the hardware structure of the high performance cluster Cray
XE6 in Stuttgart, Germany, where the simulations have been carried out.
1.3. Structure of the Thesis

The thesis is organized in the following way: After the current introduction a general
discussion on the physics of the combustion and flow dynamics in chapter 2 is presented.
In the following chapter 3 the mathematical description of chemically reacting flows is
presented with an introduction to statistical thermodynamics and the kinetic theory of
gases, from which the transport equations, as well as the equations for the molecular
fluxes are derived. A brief discussion of chemical kinetics follows, to close the mathematical description of the partial differential equations. The organization of the thesis from
sections 3.1 to 3.4 has been designed in a way to gradually guide the reader from the
most basic physical model as a starting point towards the later implemented equations.
The numerical methods used in the simulation of reacting flows are presented in chapter
4. The finite volume method is discussed, as well as methods for unsteady problems
and solution strategies for linear equation systems. Chapter 5 gives at first an overview
over the used software. In the following, the coupling library between OpenFOAM and
1.3. STRUCTURE OF THE THESIS
Cantera is presented and the new solver for DNS of chemically reacting flows with multicomponent diffusion modeling is introduced. In chapter 6 the developed solver is applied
to numerical simulations of flat premixed flames for validation purpose. The results are
compared with calculations obtained from CHEMKIN/PREMIX. Finally the solver is
used for DNS of 3D turbulent premixed flames and the results are briefly discussed in
chapter 7. The thesis finishes with a conclusion and a perspective.
2. Basic Aspects of Turbulent

Combustion
Combustion is a highly exothermic, chemical reaction between a fuel and an oxidizer.
There are two general types of combustion processes which can be distinguished by the
nature of their mixing state. In a premixed combustion the composition of fuel and
oxidizer is spatially homogeneous while in a non-premixed combustion the unburnt fuel
composition varies in space such that mixing and chemical reaction occur simultaneously. Between these two extremes one often encounters so called partially-premixed
flames. Secondly, combustions may be classified with respect to the underlying flow
regimes which can be either laminar or turbulent. Latter is encountered in most practical combustion systems such as combustion engines, aircraft engines, industrial burners
and furnaces.
In the following sections a brief description of the different flame types with respect to
their mixing state and fluid motion is presented and some fundamental definitions are
introduced.
2.1. Structure of Laminar Premixed Flames

A flow is defined laminar if the fluid layers move smoothly in such a way that they remain
ordered and do not mix in the plane perpendicular to the direction of the flow [19].
Whether a flow is laminar or not depends on the ratio of the inertial and viscous forces
which is used to define the Reynolds number
Re =
lc v
,
(2.1)
where lc is a characteristic length, v is the flow velocity, is the mass density and
is the dynamic viscosity. If a critical value of Re is reached the laminar flow starts to
become unstable and changes to turbulent flow. For example in an internal channel flow,
this happens at Re 2000.
Laminar premixed flames propagate towards a premixed mixture of fuel and air. They
occur in gas ranges, heating appliances and Bunsen burners. The understanding of
laminar premixed flames is a prerequisite to study turbulent premixed flames [10]. Due
its simple configuration laminar premixed flat flames have been an important object
in numerical combustion science. Therefore a lot of specialized numerical codes, like
CHEMKIN/PREMIX [20] and Cantera [17], exist.
2.1. STRUCTURE OF LAMINAR PREMIXED FLAMES
Fig. 2.1.: Species and temperature profiles for a laminar, premixed flat methane-oxygen flame. Figure
taken from [4].
In order to distinguish different mixtures of fuel and oxidizer the air number is defined
as [10]
air =
Xair /Xf uel

,
Xair,stoich. /Xf uel,stoich.
(2.2)
where Xair (Xf uel ) are the mole fractions of the air (fuel) in the mixture, and Xair,stoich.
(Xf uel,stoich. ) are the mole fractions of the air (fuel) in the stoichiometric mixture. Based
on the value of air premixed combustions can be divided into different classes: rich
combustion with air < 1, stoichiometric combustion with air = 1 and lean combustion
with air > 1.
The schematic structure of a premixed, flat, methane-air flame is shown in Figure 2.1.
It consists of three main regions. In the chemically inert preheating zone, heat released
from the reaction is transported by conduction. The reaction layer is typically thin and
called fuel consumption layer or inner layer with thickness and temperature TO . This
layer is responsible for keeping the reaction process alive. Here the fuel is consumed and
the radicals are depleted by chain-breaking reactions. The inner layer temperature TO
corresponds to the crossover temperature between chain-branching and chain-breaking
reactions. In the oxidation layer the final oxidation to the products is accomplished
and the temperature reaches its maximum. The thickness of the reaction zone can be
calculated from the temperature profile with [21]
=
T2 T1
.

max T

x
(2.3)
2.2. CHARACTERISTICS OF TURBULENT FLOW
In a steady flow of premixed gas, the flame propagates upstream until the flow velocity
normal to the flame front is equal to the laminar flame speed sL . The laminar flame
speed can be calculated from the integral of the burning rate across the flame brush
with [21]
Z +
1
sL =
rf uel dx .
(2.4)
0 Yf uel,0
The laminar burning velocity dependence on the fuel type, the air number, the pressure, the unburnt flow temperature, etc. and is a very important characterization of
combustions [4].
2.2. Characteristics of Turbulent Flow

As mentioned in the section before the occurrence of turbulent flow essentially depends
on the ratio of inertial and viscous forces in the fluid. In a turbulent flow the destabilizing inertial forces exceed the stabilizing viscous forces such, that the generation rate
of vortices is higher then the viscous dissipation rate [19]. Therefore, turbulence arises
from instabilities associated with large Re numbers.
To further characterize the flow, a turbulent Reynolds number ReT is introduced [22]:
ReT =
v 0 lT
(2.5)
Here, v 0 is the turbulence intensity or velocity fluctuation, is the kinematic viscosity

and lT is the turbulent length scale. The turbulent Reynolds number represents the
ratio of turbulent transport to molecular transport of momentum.
Compared to laminar flows, turbulent flows have no well-defined values for velocities and
physical scalars. The flow regime is characterized by transient chaotic and stochastic
property changes. Variations of velocity lead to fluctuations in physical scalars such as
density, temperature or transport quantities. Since the physics of turbulence are not
fully understood the integration of chemical reactions is even more difficult.
2.2.1. Energy Cascade and Kolmogorov Scales

Turbulence can be considered to be composed of eddies of different sizes. Large eddies
are unstable and break up, transferring their energy E to somewhat smaller eddies which
break up again. This energy cascades continuous until the Reynolds number of the eddy
is sufficiently small that the eddy motion is stable and molecular viscosity is effective in
dissipating the kinetic energy [22].
Figure 2.2 shows the energy spectrum measured for all wavelengths k. For small
wavenumbers corresponding to large scale eddies the energy per unit wavenumber increases with a power law between k 2 and k 4 . This range is not universal and is determined by large scale instabilities, which depend on the boundary conditions of the flow.
The spectrum attains a maximum at a wavenumber that corresponds to the integral or
2.3. INTERACTION BETWEEN PREMIXED FLAMES AND TURBULENT FLOW
Fig. 2.2.: Schematic representation of the turbulent kinetic energy spectrum E as a function of the
wavenumber k.
turbulence macroscale l0 , since eddies of that scale contain most of the kinetic energy.
For larger wavenumbers corresponding to the inertial sub-range the energy spectrum decreases following the k 5/3 law. Finally, there is a cutoff at the Kolmogorov microscale
lK . In the viscous sub-range, the energy per unit wavenumber decreases exponentially
owing to viscous effects. Therefore, the Kolmogorov microscale represents the smallest
length scales associated with a turbulent flow [4].
To determine the necessary grid resolution and time step in numerical simulations of
turbulent flows one has to know about the value of the smallest scales. The Kolmogorov
microscale for length lK and time K are defined in [22] to
lK = ( 3 /)1/4 and K = (/)1/2 ,
(2.6)
in which is the rate of dissipation of turbulent kinetic energy defined to

(*
= 2 hsij sij i =
vj
vi
+
xj xi
+*
vj
vi
+
xj xi
+)
(2.7)
and is the kinematic viscosity. Here, sij is the fluctuating strain rate and hi is the
expectation.
2.3. Interaction between Premixed Flames and

Turbulent Flow
A highly exothermic chemically reaction and turbulence can have large influences on each
other. The flow is strongly accelerated by passing the flame front due to the thermal
expansion. In addition, turbulent eddies wrinkle the flame front (as illustrated in figure
2.3) and enhance the chemical reaction. In some cases on the other hand, the flow can
completely inhibit the chemical reaction and lead to flame quenching [4].
Fig. 2.3.: Kinematic interaction between a turbulent eddy and a propagating flame front.
2.3.1. Turbulent Flame Speed

Unlike the laminar flame speed, which depends mostly on the thermal and chemical
properties of the unburnt mixture, the turbulent flame speed sT also depends on the
transient interaction between turbulence and chemical reaction. As shown in figure 2.4
the instantaneous flame front represents a large area AT propagating with the locally
laminar flame speed sL . One can see, that the turbulent flame front is much greater than
For an observer traveling with the flame, the turbulent
the mean flame area AT A.
flame speed can be defined as the velocity at which reactants enter the flame zone in a
direction normal to the flame. Thus, the turbulent flame speed can be expressed as
T ,
m
= 0 AT sL = 0 As
(2.8)
with m
equal the mass flow rate of the unburnt gas mixture and 0 equal the unburnt
mixture mass density. Consequently the ratio of turbulent flame speed to local laminar
flame speed is equal to the ratio of wrinkled flame area to the mean flame area from
which sT can be calculated to
AT
sT
AT
= sT = sL .
sL
A
A
(2.9)
10
In addition, the turbulent burning velocity can be calculated with

Z +
1
sT =
rF l dV ,
A0 Yf uel,0
(2.10)
which is similar to equation (2.4) for the laminar burning velocity.
Fig. 2.4.: Schematic representation of the turbulent flame velocity. Figure similar to [4].
2.3.2. Flame Stretch Rate

The flame stretch rate is an important quantity in the understanding of flame phenomena
such as extinction and the local structure of turbulent flames. In general, the flame
stretch rate reads [21]
1 dAT
,
(2.11)
=
AT dt
which describes the fractional rate of change of a flame surface element AT . Note, that
this equation uses the substantial derivation of the surface to include the convective
change. Because of the highly unsteady flame surface, AT can not be easily calculated.
Therefore, different approaches for equation (2.11) have been developed (see e.g. [23]).
Here, an expression for the flame stretch rate of stoichiometric premixed flames is derived
based on chemical reaction rates. Assuming, that the total inlet mass of reactants is
consumed by chemical reactions with rate r in the flame volume VF l. and taking equation
(2.8) into account, one can write
m
f uel =
Z
VF l.
T .
rf uel (V ) dV = 0 Yf uel,0 AT sL = 0 Yf uel,0 As
(2.12)
With A = const., sL = const. and 0 = const. it must apply, that

sT AT
Z
VF l.
rf uel (V ) dV
rf uel,i Vi = VCV
rf uel,i ,
(2.13)
with VCV as the volume of a single cell which is equal for all cells in an equidistant mesh.
Consequently, it must apply locally in cell i, that
(sT )i (AT )i VCV (rf uel )i .
(2.14)
11
Finally, by inserting the expressions in equation (2.11) the flame stretch rate can be
expressed as
1 dAT
1 drf uel
=
=
.
(2.15)
AT dt
rf uel dt
This equation describes the fractional rate of change of the fuel consumption rate rf uel .
Note, that this equation is valid for every reactant.
2.3.3. Combustion Regimes

Diagrams defining regimes of premixed turbulent combustion in terms of velocity and
length scale ratios have been proposed by many authors (see e.g. [4] for details). Here, a
regime diagram, figure 2.5, is discussed following [21]. Therefore, two additional reduced
numbers are used, the turbulent Damkhler number DaT and the Karlowitz number Ka,
defined as:
lT /v 0
T
=
c
/sL
c
2
Ka =
= 2
K
lK
DaT =
(2.16)
(2.17)
DaT describes the ratio between the turbulent time scale T and the characteristic chemical time scale c . The interaction of the chemical reaction with the dissipative turbulent
structures of the flow field is described with Ka, which is the ratio of the characteristic
chemical time scale and the characteristic time scale of the smallest Kolmogorov eddy
K . Additionally, the turbulent Reynolds number introduced before is used in the diagram 2.5.
Fig. 2.5.: Schematic classification of turbulent combustion regimes.
12
By using the reduced numbers ReT , Ka and Da, turbulent premixed combustions can be
classified into three major categories. The thin reaction zone regime or flamelet regime
exists for Ka < 1, which means, that the flame thickness is smaller than the Kolmogorov
length scale. It assumes, that the flame structure is not effected by turbulence. The
flame sheet is wrinkled by vortices not small enough to enter it. At Da > 1 and Ka > 1,
the smallest eddies can enter into the flame front and enhance the diffusion inside of the
flame, which leads to a thickened flame. Here, the turbulent flow is intense enough to
generate eddies able to effect the structure of the reaction zone. Since it is expected,
that the lifetime of such eddies is very short, their impact on the reaction zone is thus
limited. If the reaction time is greater than the time needed for turbulent fluctuation,
a well stirred reactor regime occurs. It is characterized by Da < 1 and Ka > 1.
Consequently, most of the eddies can enter into the reactive-diffusive layer during their
lifetime and enhance the diffusivity within the reaction zone. It is even possible, that
nearly all of the eddies are embedded into the reaction zone. Hence, the flame front is
very thick and behaves similar to an ideally stirred reactor. Since in a turbulent flow a
wide range of length scales occur, a turbulent premixed flame is not represented by a
single point in figure 2.5, but by a zone that may cross the different regimes.
3. Mathematical Description of
Chemically Reacting Flows
In the following an introduction to statistical thermodynamics and the rigorous kinetic
theory of gases (see [1, 24, 25]) is presented. The focus in this section is on the ChapmanEnskog formulation and will not go through detailed derivations of all results. After,
the transport equations, as well as the equations for the molecular fluxes are derived
from the Chapman-Enskog formulation. A brief discussion of chemical kinetics follows,
to close the mathematical description of the partial differential equations.
3.1. Statistical Thermodynamics and the Rigorous

Kinetic Theory
Statistical thermodynamics apply probability theory to a large number of particles in order to add a molecular-level interpretation of the macroscopic thermodynamic quantities
described by classical thermodynamics. It therefore provides a mathematical description
of thermodynamic, chemical kinetic and transport quantities which are needed in the
numerical simulation of chemically reacting flows [24].
The rigorous kinetic theory of gases is used to relate the motion or kinetic energy of
a large number of molecules to macroscopic thermodynamic quantities [1]. The theory
makes the assumptions that
(i) only binary collisions occur and that therefore the pressure of the gas mixture is
moderate
(ii) quantum-mechanical effects are negligible and that therefore temperature has to
be moderate
(iii) the dimensions of the containing vessel and any obstacle therein are large compared
to the molecule mean free path which means that the average distance separating
the gas particles is large compared to their size.
Strictly speaking, the rigorous kinetic theory of gases applies only for elastic collisions
but it can be extended to take account for the effect of inelastic collisions.
Since the assumptions made by the rigorous kinetic theory include the properties for an
ideal gas the equation of state is given by
pV = N kB T = nNA kB T = nRT
(3.1)
where p is the pressure, V is the volume, N is the number of molecules, kB the Boltzmann
constant, T the temperature, n the number of moles, NA is the Avogadro constant and
R is the universal gas constant.
3.1. STATISTICAL THERMODYNAMICS AND THE RIGOROUS KINETIC THEORY
14
3.1.1. The Boltzmann Equation

The Boltzmann equation given by [26] or [25]1 is used as a starting point:
fk
~ k + 1 F~k fk
+ ~vk f
t
mk
~vk
X h (+)
j
()
kj kj
(3.2)
As shown in figure 3.1, equation (3.2) describes the time evolution of the velocity distribution function fk (~x, ~vk , t). The quantity fk (~x, ~vk , t)d~xd~vk is the mean number of
molecules of species k at time t, which are located in the volume element d~x around ~x,
and which have velocities within the range d~vk around ~vk . The population of molecules
propagating to a different position in phase space after a time dt is increased by some
(+)
(+)
collisions kj and decreased by others kj . The kj involve the intermolecular potential energy function, e.g. the Lennard-Jones-Potential, and all details of the collision
trajectories [19]. The quantity F~k is an external force acting on a molecule of species k.
Fig. 3.1.: Graphical expression of the Boltzmann equation.
The Boltzmann equation can be viewed as a transport equation in the six-dimensional

position-velocity space, where the right hand side of equation (3.2) serves as a source
term. In cases where the function f can be calculated all properties of a gas mixture,
like its local temperature of the chemical properties are completely known [19].
Different from [25] from now vector-tensor notation is used, e.g. [19]
3.1. STATISTICAL THERMODYNAMICS AND THE RIGOROUS KINETIC THEORY
15
In the following subsections, the Boltzmann equation is used to derive the transport
equations for mass, momentum and energy, as well as formulations for the molecular
fluxes and the transport properties.
3.1.2. Enskogs General Transport Equation

A general transport equation for a physical quantity can be derived from the Boltzmann
equation without actually determining the form of the distribution functions fk . By
multiplying equation (3.2) with the quantity k of species k and by integrating over ~vk
one obtains Enskogs general transport equation for the physical quantity k associated
with the i-th kind of molecule [25]

F~k k
(nk k ) ~
k
~
+ (nk k~vk )nk
+ ~vk k +
=
t
t
mk ~vk
Z
X h (+)
j
kj
()
kj
(3.3)
d~vk ,
where nk is the number of moles of species k and where the overbars indicate averaged
quantities. A summation of equation (3.3) over all species k gives the transport equation
for quantity of the entire gas mixture.
P
For conserved quantities, such as momentum ( = k mk~vk ), kinetic plus internal energy
P
P
(int)
( = k 21 mk vk2 + ek ) or the global mass ( = k mk ) it is explicitly shown in [25],
that the right hand side of equation (3.3) vanishes as it should. Note that, however, that
the mass of one species k is not conserved in case of a chemical reaction: By choosing
k = mk the sum of the right hand side of equation (3.3) represents the conversion rate
of molecules of species k which is in general not zero.
In section 3.2 Enskogs general transport equation (3.3) is used to derive explicit formulations for the transport equations of mass, momentum and energy.
3.1.3. Chapman-Enskog Theory

To obtain rigorous expressions for the molecular fluxes and the corresponding transport
coefficients a solution to the Boltzmann equation (3.2) has to be found. Here, a practical solution given by Enskog is presented [26]. By expanding the velocity distribution
[r]
function fk with different orders of approximation fk one obtains
[0]
[1]
[2]
[r]
fk = fk + fk + 2 fk + r fk ; ,
(3.4)
with as an ordering parameter. One can find an analytical solution for the 0-th
approximation of fk [25]. Now the first-order approximation to fk can be written in
terms of a perturbation function k as
[1]
[0]
fk = fk k ,
[1]
(3.5)
and substitute fk in equation (3.2) to obtain an integro-differential equation for . As

shown by [25], has a specific form with two unknown scalar functions [1]. By expanding
3.2. TRANSPORT EQUATIONS FOR CHEMICALLY REACTING FLOWS
16
the unknown scalar functions in finite series of Sonine polynomials, one can use the
results to obtain first-order approximations of the molecular fluxes and the transport
properties [25]. Approximations of higher order are not worthwhile since the calculation
is very time consuming and the additionally gained accuracy is small [16].
3.2. Transport Equations for Chemically Reacting

Flows
For numerical simulation of chemically reacting flows the transport equations constitute
the mathematical and physical core of the simulation. Here, the transport equations
are presented in a different way from the traditional continuum mechanics approach,
e.g. [19]. Following [25], it is shown, that, based on the Boltzmann equation respectively
Enskogs general transport equation presented before, the equations of mass, momentum
and energy are direct consequences of the conservation laws for mass, momentum and
energy.
Generally, a transport equation can be written as
()
|t {z }
Transient term
~ (~v )
~ ( )
~ +
=
|
{z
Convection term
{z
Diffusion term
(3.6)
|{z}
Source term
with the conserved quantity and transport coefficient . The transient term accounts
for the accumulation of and the source term represents any sources or sinks in the
concerned control volume. The fluxes over the control volume faces are described by
either the convective term due to the flow velocity field ~v or the diffusion term due to
its gradients.
3.2.1. Equations of Mass

By choosing k = mk in equation (3.3) one obtains, after some intermediate transformations, the transport equation of the mass fraction for each species k in a mixture of
K species [25]
~ (Yk~v )
~ ~jk + rk
(Yk ) =
t
k = 1, 2, 3, , K ,
(3.7)
where Yk = mk /m is the mass fraction of k-th species, is the mixture density, ~v is

the velocity, rk is the reaction rate of species k (see section 3.4) and ~jk is the molecular
diffusion flux (see section 3.3).
Summing equation (3.7) over all Yk and using the relations
N
X
k=1
Yk = 1,
N
X
k=1
~j k = 0,
N
X
k=1
rk = 0
(3.8)
3.2. TRANSPORT EQUATIONS FOR CHEMICALLY REACTING FLOWS
17
one obtains the transport equation for the total mass of a multicomponent reacting
mixture [19, 25]
~ (~v ) .
=
(3.9)
t
3.2.2. Equation of Momentum

The momentum of a system of two colliding molecules is always conserved in any collision, even if a chemical reaction occurs. By choosing k = mk~vk in equation (3.3) one
obtains the equation of momentum with gravitation ~g as the only external force as [25]
~ (~v~v )
~ + ~g ,
(~v ) =
t
(3.10)
~ = p
~ +
~ ~
(3.11)
where
is the rate of momentum addition by molecular transport with ~ as the molecular momentum flux vector (compare section 3.3).
3.2.3. Equation of Energy

(int)
into equation (3.3) and summing over all species k one

Applying k = 12 mk vk2 + uk
obtains the transport equation for energy. By using the thermodynamic definition of
enthalpy the transport equation for specific enthalpy is given as [1, 25]
~ (h~v )
~ ~q : ~
~ v + Dp + Q source ,
(h) =
t
Dt
(3.12)
where ~q is the molecular heat flux (compare section 3.3), Q source is a combination of
is the reversible rate of enthalpy due
all heat source terms like the radiative flux, Dp
Dt
~
to compression and ( : ~v ) is the irreversible rate of enthalpy due to viscous dissipation.
Further with the sensible enthalpy hs and the chemical enthalpy hc the total enthalpy
h can be expressed as [21]
h = hs + hc =
Z T
T0
N
X
cp (T ) dT +
{z
h0f,k Yk ,
(3.13)
k=1
sensible enthalpy
{z
chemical enthalpy
with the heat capacity for an ideal gas at constant pressure [24]
cp (T ) =
h
T
(3.14)
Neglecting viscous dissipation and heat source terms like the radiative flux, one obtains
the equation of energy in terms of sensible enthalpy to
~ (hs~v )
~ ~q + Dp + ~qreaction .
(hs ) =
t
Dt
(3.15)
3.3. MULTICOMPONENT MOLECULAR TRANSPORT
18
The separation of sensible enthalpy from the heat of formations gives this equation the
advantage over equation (3.12), that the heat of formations become the only source term
~qreaction =
N
X
h0f,k rk .
(3.16)
k=1
Additionally, this equation can be used as a definition for the heat release rate.
To obtain the temperature from the sensible enthalpy, the definition of the sensible
enthalpy (3.13) is solved for T . In general, the mixture averaged specific heat cP (T ) =
PK
k=1 cp,k (T )Yk is a function of T and therefore, the integral in (3.13) has to be solved
iteratively [21]. The temperature dependencies of the pure species specific heat capacities
cp,k are fitted by 4-th order NASA polynomials:
k
cp,k (T )M
= a1,k + a2,k T + a3,k T 2 + a4,k T 3 + a5,k T 4
R
(3.17)
Here, an,k are coefficients of the k-th species which has to be given as input.
3.3. Multicomponent Molecular Transport

Gradients exist in a gas under non-equilibrium conditions in one or more physical quantities: composition, mass averaged velocity and temperature. These gradients cause the
molecular or diffusive transport of mass, momentum and energy through the gas.
In this section, expressions for the diffusive flux vectors in the transport equations (3.7),
(3.10) and (3.15) will be derived based on the Chapman-Enskog formulation in subsection
3.1.3. In addition, a formulation for multicomponent transport coefficients laid out
by [16] for computational purpose will be given.
3.3.1. Species Mass Flux

Inserting the first approximation of the distribution function in terms of the perturbation
function mentioned in subsection 3.1.3 and making use of the polynomial expansions and
orthogonality relations [25] one obtains the mass diffusion velocity [1, 16]
V~k =
N
1 X
DT 1 ~
M j Dkj d~j k T
,
Yk T
Xk M j=1
(3.18)
where Xk = nk /n is the molar fraction of species k, M is the mean molecular weight, M j

is the molecular weight of species j, Dkj is the ordinary diffusion coefficient of species
k through species j, DkT is the thermal diffusion coefficient of species k and the driving
force d~j is defined as
~ .
~ j + (Xj Yj ) 1 p
d~j = X
(3.19)
p
19
Finally, the multicomponent mass flux in equation (3.7) is given by

N
X
1~
~jk = Yk V~k = Yk
M j Dkj d~j DkT T
.
T
Xk M j=1
(3.20)
Equation (3.20) shows that the mass flux can be caused by three different phenomena.
The ordinary diffusion flux due to a gradient in concentration, the flux due to a pressure
gradient and the thermal diffusion flux due to a temperature gradient. The flux due to
a pressure gradient is very small compared to other effects and can be neglected [25].
Additionally, a flux due to external forces can occur for example in charged mixtures.
Mass Diffusion Coefficients
Based on the theory provided by [25] and laid out for computational purpose by [16],
the multicomponent diffusion coefficients Dkj and the multicomponent thermal diffusion
coefficients DkT in equation (3.20) are computed from a system of equations defined as
(L)-matrix which consists of nine sub-matrices:
0
L00,00 L00,10
0
a1

00
L10,00 L10,10 L10,01 a1 = X

10
X
0
L01,10 L01,01
a101
(3.21)
with the right hand side vector composed by the mole fraction vectors Xk .
Every component of the L matrix is a K K, with total number of species K, matrix.
For example the elements of the matrix L00,00 are calculated with
L00,00
=
jk
N
Xl
16T X
[M k Xk (1 jl M j Xj (jk kl )] ,
25p l=1 M j Djl
(3.22)
where is a small value due to prohibit species concentration values of exactly zero and
Djk is the binary diffusion coefficients calculated from Chapman-Enskogs theory with
3 3
T /mjk
3 2kB
=
.
(1,1)?
2
16 pjk
jk
Djk
Here kB is the Boltzmann constant, mjk =
mk mj
mk +mj
(3.23)
with mk the mass of a molecule k,

(1,1)?
is the reduced molecular mass, jk is the collision diameter and jk

is the collision
integral. The collision diameter and the collision integral can be found in transport data
files, e.g GRI-Mech 1.2 [27].
With the inverse (P ) of the L00,00 -block the first-order approximation of the multicomponent diffusion coefficients is given by
Djk = Xj
16T
M (Pjk Pjj ) ,
25pM k
(3.24)
20
and the thermal diffusion coefficients are given by

DkT =
8M k Xk 1
ak00 .
5R
(3.25)
Here, R is the universal gas constant and a1k00 is calculated from the (L)-matrix. Equation (3.24) shows two characteristics of the multicomponent diffusion matrix predicted
by [25],
Djj = 0 and Djk 6= Dkj .
(3.26)
For further details on the (L)-matrix one should consult [1, 16].
3.3.2. Momentum Flux

The laminar momentum flux ~ for a Newtonian fluid is given by [19]:
~ ~v )
~ v + (~
~ v ) ] + 2 (
~ = [~
3
(3.27)
Here, is the mixture viscosity and is the identity matrix. Note, that the momentum
flux is not derived from expressions obtained from the rigorous kinetic theory of gases,
but from hydrodynamic equations [25].
The viscosity of a mixture is calculated from Wilkes mixture rule [28] modified by [19]
=
K
X
k=1
Xk k
,
j=1 Xj kj
with
kj
1
Mk
= 1+
Mj
8
(3.28)
PK
!1/2
k
1 +
j
!1/2
Mj
Mk
!1/4 2
and the pure species viscosity derived from the rigorous kinetic theory of gases
5 mk kB T
k =
.
16 k2 (2,2)?
kk
(3.29)
(3.30)
kB is the Boltzmann constant, k is the collision diameter for the k-k interaction poten(2,2)?
tial, mk is the mass of molecule k and kk the reduced collision integral.
It has been shown by [25] that Wilkes mixture rule is a very good approximation to
the equation derived from the rigorous kinetic theory. One should note, that the pure
species viscosities calculated with equation (3.30) are still based on the rigorous kinetic
theory.
3.4. REACTION KINETICS
21
3.3.3. Heat Flux

By using the rigorous kinetic theory of gases for the energy flux ~q in a multicomponent
gas mixture, the heat flux vector is given by [16] to
~ +
~q = 0 T
N
X
k=1
~jk hk
N
X
RT
DkT d~k ,
M
X
k
k
k=1
(3.31)
where 0 is the special mixture thermal conductivity, ~jk is the diffusive mass flux discussed before, hk is the specific enthalpy of species k and dk is the diffusive driving force
see equation (3.19)
As one can see from equation (3.31) heat can be transported through a multicomponent
gas caused by a temperature gradient known as Fouriers law, a diffusive flux by species,
and the reciprocal process to thermal diffusion called Dufour effect [25]. Like the ordinary
and thermal diffusion coefficients, the thermal conductivity 0 is calculated with a1i10 and
a1i01 from the L-matrix (3.21) with
0 = 0,trans. + 0,inter. = 4
Xk (a1i10 + a1i01 )
(3.32)
It consists of a translational and a internal part due to the molecular vibration.
3.4. Reaction Kinetics

Dependent on physical conditions and used species a chemical reaction is proceed by a
series of elementary reactions, called the reaction mechanism. Generally, each elementary
reaction can be described as
K
X
k=1
K
X
ki Mk *
)
k=1
ki 00 Mk
(i = 1, , I) ,
(3.33)
where K is the total number of chemical species, I is the total number of chemical reactions, Mk is the name of species k and ki is the stoichiometric coefficient of species
k in the forward direction respectively reverse direction in the i-th reaction.
To obtain the rate of progress i of a two-body reaction i the reaction rates of for- and
backward reaction have to be combined to [1]
i = kf,i
K
Y
[Mk ]ki kr,i
k=1
K
Y
00
[Mk ]ki ,
(3.34)
k=1
where kf,i and kr,i are the rate constants for the forward and reverse direction of reaction
i and [Mk ] is the concentration of species k.
The rate constant k is computed by a modified three-parameter Arrhenius form
k = A? T exp(Ea /RT ) ,
(3.35)
3.5. CHEMICAL TIME SCALES
22
where A? is the pre-exponential constant, is the temperature exponent and Ea is the

activation energy. Obviously, the rate constant depends strongly on the temperature.
Equation (3.34) together with the reaction equations (3.33) and the rate constant equations (3.35) build the chemical kinetic model. This system of ordinary differential equations have to be solved numerically [29]. The parameters needed for chemical kinetic
models are listed in data files like in the appendix B.
Finally, the chemical source for species rk is given by
rk =
I
X
i
(ki00 ki0 )i .
(3.36)
3.5. Chemical Time Scales

The needed time step in a numerical simulation is determined by the smallest time
scale on which the regarded physical and chemical phenomena occur [3, 21]. Therefore,
knowledge about the magnitudes of chemical time scales is important and offers the
opportunity to adjust the time step.
Because of the many approaches a definition by [30] will be used to calculate chemical
time scales:
!1
(x)
.
(3.37)
(x) = (x)
t
Here, the time scale (x) is determined with the variation of a characteristic quantity
on a specific location with time t.
For highly unsteady chemical reactions like a moving flame front and negative mass
fraction gradients, this equation can be further converted with
=
=
v,
t
x t
x
(3.38)
into a more suitable equation,

(x) = (x)

!1
(x)

v
x
(3.39)
This equation describes the variation of a quantity with location x multiplied by the
velocity v. One should note, that equation (3.39) is mathematically not clearly defined
for mass fraction gradients equal zero.
For turbulent premix flames it is assumed, that the chemical time scale is not affected
by turbulence as long as the smallest eddies can not penetrate into the thin reaction
zone [4].
4. Numerical Solution of Partial

Differential Equations
The transport equations in section 3.2, together with the expressions for the fluxes and
transport properties in section 3.3 and the expressions for the chemical sources in section 3.4, discussed in the chapter before, form a complicated system of partial differential
equations. Due to the demand on accuracy one can not include simplifications to obtain
an analytic solution but has to approximate the solution numerical.
The components of a numerical solution method include a physio-mathematical model,
which was introduced in the sections before, a discretization method, a numerical grid
and a solution method [7].
A spatial discretization method has to be used to approximate the position dependent
parts of the partial differential equations by a system of algebraic equations at discrete
locations in space. There are many approaches like finite difference, finite element or
finite volume methods. Due to its importance in computational fluid dynamics and its
R
implementation in OpenFOAM
only the finite volume method (FVM) is presented following [7, 31]. Additionally, in unsteady problems a discretization of the time dependent
parts of the partial differential equations has to be used.
4.1. The Finite Volume Method

In the FVM the solution domain is first subdivided into a finite number of contiguous control volumes. Then the conservation equations are applied in integral form to
each control volume and volume integrals are transformed into surface integrals by using Gausss Theorem. After, the integrals are approximated using suitable quadrature
formulas. Since the variable values are calculated only at the centroid of each control
volume, interpolation is used to express variable values at the control volume surfaces.
As a result, one obtains an algebraic equation for each control volume. Finally, the
resulting matrix is solved directly or iteratively.
As mentioned above, by integrating the stationary part of the generic transport equation
(3.6) over a volume and using Gausss Theorem to transform volume integrals containing
a divergence term to surface integrals, a stationary generic conservation equation for a
quantity in integral form is obtained to
Z
S
~=
~v dS
Z
S
~ dS
~+
()
Z
V
q dV ,
(4.1)
4.1. THE FINITE VOLUME METHOD
24
where S is a surface, V is a volume, is a transport coefficient, q is a source term and

is some physical quantity. The used parameters in FVM are described in figure 4.1.
Fig. 4.1.: Two control volumes with centers P and N are connected through face f with face normal
~f . Figure taken from [31].
S
Here, P and N are the control volume centers, f is the connecting cell face, d~ is the
~f is the face area vector.
distance between the cell centers and S
Figure 4.2 shows how a structured grid divides the solution domain of a duct into a finite
number of small control volumes. Equation (4.1) applies to each single control volume.
A sum over all control volumes results in the global conservation equation, since surface
integrals over inner control volumes cancel out. This provides a principal advantage of
FVM thus global conservation is build in.
4.1.1. Approximation of Integrals

To obtain an algebraic equation for a control volume from equation (4.1), the integrals
need to be approximated.
Surface Integrals
By adding the integrals over the six control volume faces together one obtains the net
flux through the control volume boundary:
Z
S
~=
w
~ dS
XZ
k
Sk
w
~ dS~k ,
(4.2)
where w is the convective or diffusive flux vector from (4.1). Usually, values for velocity
and physical quantities are taken from the last time step. To preserve conservation,
control volumes are not allowed to overlap so that each control volume face is unique to
the two control volumes which lie on either site. Each cell face has only one owner and
4.1. THE FINITE VOLUME METHOD
25
Wall
Inlet
Wall
Outlet
Wall
Wall
Fig. 4.2.: Example of a 2D, structured, non-orthogonal grid to simulate the flow through a duct. Figure
similar to [7].
one neighbor [31].

R
The cell face integral Sk in equation (4.2) can not be calculated exactly, since only discrete values of f at the control volume center exist. Therefore w has to be approximated,
e.g. by the midpoint rule:
Z
Sf
wf dSf = wf Sf wf Sf
(4.3)
Here, the integral is approximated as a product of the mean value over the surface
wf , which has to be approximated itself, and the cell-face area Sf . To obtain higher
approximations, the flux has to be calculated at more than two locations.
Volume Integrals
The last term in equation (4.1) requires integration over the volume. The simplest
second-order approximation can be obtained by replacing the volume integral by the
product of the mean value of the integrand, which is approximated as the value at the
center of the control volume:
Z
V
q dV = qV qP V ,
(4.4)
where qP is the value of q at the center of the control volume. To obtain approximation
of higher order, one has to interpolate between values of q besides the value at the center.
4.1.2. Interpolation and Differentiation Procedures

In order to obtain the approximations to the integrals, values of variables at locations
other than control volume centers are needed. Since numerous possibilities are available,
only a few while be further discussed.
4.2. METHODS FOR UNSTEADY PROBLEMS
26
Linear Interpolation Scheme (CDS)

The simplest second-order scheme to obtain e at the control volume face center is the
linear interpolation between the two nearest nodes,

f = N
xf xP
xN xP
xf xP
1
xN xP
+ P
(4.5)
A Taylor series expansion of N about the point xP gives

xf xP
xf xP
f = N
+ P 1
xN xP
xN xP

2
x2
(xf xP )(xN xf )
+H ,
(4.6)
where H denotes higher-order terms. One can see, that the leading truncation error is
proportional to the square of the gird spacing. Therefore the scheme is second order
accurate.
Upwind Interpolation Scheme (UDS)
In upwind interpolation, f is approximated as
f =
if(~v S~f )f > 0

if(~v S~f )f < 0
(4.7)
A Taylor series expansion about P gives for a Cartesian grid and (~v S~f )f > 0
f = P + (xf xP )
x
!
P
(xf xP )2
+
2
2
x2
+H ,
(4.8)
One can see, that the UDS is of first order. Its leading truncation error term is diffusive
and will therefore never yield oscillatory solutions but at the expense of accuracy. Peaks
or rapid variations in the variables will be smeared out.
4.2. Methods for Unsteady Problems

Since many physical phenomenal, like turbulence, are highly unsteady, the time derivations in the transport equations are not zero. Therefore, just as the spatial derivations,
the time must be discretized. The main difference between spatial and time coordinates
is the direction of influence. Unsteady flows are elliptic in space but parabolic-like in
time. That means, that a force at any space may influence the flow anywhere else but
a force at a given instant will only affect the flow in the future.
As an example, an first order ordinary differential equation with an initial condition is
considered:
d(t)
= f (t, (t)); (t0 ) = 0
(4.9)
dt
If one finds a solution a short time t after the initial point, the solution at t1 = t0 +t
can be regarded as a new initial value for the next time step. Therefore the solution
4.3. SOLUTION OF LINEAR EQUATION SYSTEMS
27
Fig. 4.3.: Approximation of the time integral f (t) over an interval t. From left to right: explicit Euler,
implicit Euler, trapezoidal rule, midpoint rule [7].
methods advance in a step-by-step or marching manner.

Integrating equation (4.9) from tn to tn+1 = tn + t one obtains
Z tn+1
tn
Z tn+1
d
n+1
n
f (t, (t)) dt ,
dt =
=
dt
tn
(4.10)
with n+1 = (tn+1 ).

To calculate at discrete values in time, an approximation of the integral on the right
hand side has to be found. Figure 4.3 shows four relatively simple methods, the gray
area is the approximation of the integral.
Using the value of integrand at the initial point one obtains the explicit or forward Euler
method. Instead, by using the final point in the estimation one obtains the implicit or
backward Euler method. By using the midpoint one obtains the midpoint rule which
is the basis to the leapfrog method. Finally, by using a straight line interpolation one
obtains the trapezoid rule, which is the basis to the Crank-Nicolson method. Implicit
methods, which means including future times steps into the approximation, are more
stable than explicit.
A necessary condition for convergence is the CourantFriedrichsLewy (CFL) condition.
It is defined with the Courant number Co
Co =
~v t
Comax ,
x
(4.11)
which represents the fraction of the cell that the flow advances during a time step.
Normally, the maximum Courant number is taken to Comax = 1 [7].
4.3. Solution of Linear Equation Systems

Summing all discretized equations obtained from equation (4.1) for each control volume
leads to system of algebraic equations
~=Q
~ ,
A
(4.12)
4.3. SOLUTION OF LINEAR EQUATION SYSTEMS
28
~ is the solution vector and Q

~ contains the
where A is the matrix with the coefficients,
source terms.
To solve these equations, direct methods and iterative methods can be used. With direct
methods, the solution is determined in one single step. The basic method is the Gaussian elimination, which divides large systems of equations systematically into smaller
ones. For full matrices, it is one of the fastest and most accurate methods, but usually
matrices in fluid mechanics are sparse. In addition, the rounding errors can grow rapidly
and therefore the accuracy may be insufficient.
~ is corrected until a convergence criterion is
In iterative methods the solution vector
~ n is obtained which does not
satisfied. After n iterations an approximate solution
satisfy the equations (4.12) exactly. Instead, there is a non-zero residual n :
~n = Q
~ ~n
A
(4.13)
To iteratively drive the residual to zero one can write

~ n+1 = N
~ n + B.
~
M
(4.14)
~ n+1 =
~n =
~ it must be
Since, at convergence,
P A=M N
and
~ =PQ
~ ,
B
(4.15)
where P is a non-singular pre-conditioning matrix. Suitable for iterative methods are

the Jacobi-Iteration, the Gauss-Seidel-Iteration and relaxation methods. There are also
methods which work partly directly and partly iteratively, for example, the block iteration, LU (Lower-Upper) iteration and the CG (conjugate gradient) method.
5. Implemented Solver in OpenFOAM

and the Cantera Interface
This chapter starts with a short discussion of the used software packages OpenFOAM
2.0.1 and Cantera 1.8. Thereafter, the developed interconnection library between OpenFOAM and Cantera is presented. Finally, the implementation of the species and sensible
enthalpy equations as well as the integration of the coupling library into existent solvers
is showed.
5.1. Software Packages

OpenFOAM Version 2.0.1
R 1
OpenFOAM
is an open source object-oriented software package written in C++03.
It provides a framework to develop numerical solvers in continuum mechanics. Main
advantages of the design are the overloaded operators, allowing expressive and versatile
syntax for implementations of complex physical models, the extensive pre- and postprocessing tools including complex geometry handling and data in-/output and a wide range
of implemented discretization schemes and boundary conditions. Additionally, OpenFOAM provides a good parallelization, a high accuracy and is available free of charge.
It has been validated many times and is widely used. Fundamental developments are
done by the OpenFOAM Foundation with contributions from a large community. The
complete source code is licensed under the GNU General Public License (GPL) for free
use and customization [15].
In OpenFOAM the discretization of the solution domain is accomplished by the finite

volume method (FVM) presented before. For the discretization of mathematical equations OpenFOAM provides two basic classes. Implicit discretization is handled by the
class fvm. Using the different overloaded operators, e.g. laplacian or grad, of the class
fvm, adds directly to the solving matrix. Explicit discretization is accomplished by using
the class fvc. Operators from fvc add a source term to the solving matrix. The user
has the possibility to chose during runtime from a wide range of discretization schemes
which differ in accuracy, stability and error order (e.g. linear Gauss schemes in second or
fourth order for gradients, upwind or cubic divergence schemes, Crank-Nicholson for the
time and many more). OpenFOAM offers many solvers for specific cases like compressible, multiphase or reacting flows. Additionally, one can chose during runtime which
turbulence, combustion or thermodynamic model the solver should use [32].
1
Open Source Field Operation and Manipulation, registered trademark of OpenCFD Limited
5.2. CONNECTION OF OPENFOAM WITH CANTERA
30
OpenFOAM provides a framework for parallelization via domain decomposition using

different schemes or third-party applications, e.g. Scotch, and the possibility to include
different implementations of the Message Passing Interface (MPI), e.g. Open MPI or
MPICH2. The solving matrices can be solved using different strategies, e.g. the linear
preconditioned conjugate gradient (PCG) solver for symmetric matrices or the linear
preconditioned biconjugated gradient (PBiCG) solver for asymmetric matrices where
the pressure-velocity coupling can be solved using a semi-implicit method for pressurelinked equations (SIMPLE) algorithm [7, 32].
Cantera Version 1.8
Cantera from the California Institute of Technology is an object oriented software toolkit
for problems involving chemical kinetics, thermodynamics and transport processes [17].
It consists of a kernel written in C++03 which is accessible through a wide range of environments like MATLAB and C++. Cantera provides fast, efficient algorithms and is
highly customizable [33]. The complete source code is licensed under the Berkeley Software Distribution 3-Clause License (BSD-3) which has been verified as a GPL-compatible
free software license [34].
In simulation codes for turbulent, reacting flows, like OpenFOAM, Cantera can be used
to evaluate thermodynamic properties, chemical sources or transport coefficients that
appear in the transport equations. Cantera includes highly accurate equations derived
from the rigorous kinetic theory of gases presented in chapter 3. The implementation
is accomplished following [16]. Due to the mathematically expensive operations on very
stiff ordinary differential equation systems and matrices inversion, Cantera makes use
of highly specialized software packages like SUNDIALS2 , BLAS3 and LAPACK4 . To
simulate complex systems, a number of ideal reactors is available [33]. Cantera can read
thermophysical, transport and chemical data from data files in the CHEMKIN format,
see [35]. Additionally, it is possible to perform sensitivity analysis.
5.2. Connection of OpenFOAM with Cantera

Since the numerical computation of multicomponent transport coefficients is not trivial,
see [1, 16], a new implementation into OpenFOAM is not worthwhile. As mentioned in
the section before, Cantera is already capable of calculating accurate multicomponent
transport coefficients needed in the transport equations, see chapter 3. Additionally,
Cantera can read transport data files in the CHEMKIN format and the calculation of
chemical reaction rates is faster as in OpenFOAM [18], which is an important factor
while performing time-consuming direct numerical simulations. Both OpenFOAM and
Cantera use C++ as the main language. Hence, it is possible to include source code
from Cantera directly into OpenFOAM-classes respectively use Cantera functions in
OpenFOAM. Based on the Alternate Chemistry Library for OpenFOAM 1.5 (see [18])
2
Suite of nonlinear and differential/algebraic equation Solvers (SUNDIALS)

Basic Linear Algebra Subprograms (BLAS)
4
Linear algebra package (LAPACK)
3
31
the coupling interface between OpenFOAM and Cantera was rewritten to meet the requirements of OpenFOAM in Version 2.0.1.
The interface makes use of templates to work with existent solvers through the run-time
selection mechanism. The extensive use of OpenFOAM respectively C++ programming
principles enables easy development of further extensions. Additionally, the coupling
offers access to sensible enthalpies, multicomponent or mixture-averaged transport properties or chemical reaction rates. The input files containing molecular properties and
chemical mechanism data are in the classical CHEMKIN compatible format.
5.2.1. Structure of the Interface

As shown in figure 5.1, the new coupling interface consists of several classes. The structure is similar to OpenFOAMs with numerous classes derived from original source files.
The classes canteraChemistryModel and canteraHsPsiMixtureThermo offer access to cell
dependent values like the chemical source terms, the enthalpy source term or the transport properties. Therefore they are derived from the class canteraMixture, which sets
up the physical conditions in Cantera, like temperature, pressure and concentration. It
therefore contains an object of the class canteraThermo, which calculates cell dependent
thermophysical properties through Cantera while making use of an object of the class
canteraGasMixWrapper. The class canteraGasMixWrapper deploys the coupling with
original classes from Cantera. It is derived from Canteras class IdealGasMix and contains an object of Canteras class Multitransport.
To guarantee full compatibility to existent reacting flow solvers, the class canteraHsPsiMixtureThermo is derived from the original class hsCombustionThermo and the class
canteraMixture is derived from the original class basicMulticomponentMixture. Consequently one can perform downcasts5 in existent solvers to obtain access to functions
which are newly introduced through classes of the interface, e.g. mass diffusion coefficients.
5.2.2. Integration of the Interface into OpenFOAM

Since the interface offers access to nearly all functions known from the original OpenFOAM classes, one can select during runtime which library calculates physical properties
like sensible enthalpy or viscosity. This is realized with dictionary files in the case directory through OpenFOAMs runtime-selection mechanism. Note, that in general no
changes on solvers source code is needed. However, like already mentioned before, OpenFOAMs original classes do not offer mass diffusion coefficients. Hence, to implement
transport equations like in the section before, the source code of the solver has to be
changed more deeply.
Type casting is converting an expression of a given type into another type. Downcasting is converting
a bass-class reference or pointer to a derived-class. The opposite process is called upcasting [36, 37].
32
Fig. 5.1.: Inheritance and dependency diagram for the coupling library. Dotted arrow contains an
object of, normal arrow derived from.
Figure 5.2 shows the initialization sequence of the coupling library in any solver during
runtime. The order in which the class constructors are called indicates the high level of
interconnection between original classes and classes from the interface. It is important
that the start of the initialization is performed by the original function psiChemistryModelNew. This function is already implemented in existent solvers and guarantees
the connection to OpenFOAMs runtime selection mechanism. It reads input files and
chooses the chemistry class, here the canteraChemistryModel. The function hsCombustionThermoNewType chooses the thermo class from input files. Here, the class canteraHsPsiMixtureThermo is used, which calls the subclass constructors until finally the class
canteraGasMixWrapper deploys the connection to Canteras original classes. Listing A.3
in the appendix shows the implementation of the initialization.
5.2.3. Exchange of Data between OpenFOAM and Cantera

To obtain access to derived class properties, like the mass diffusion coefficients, a downcast shown in listing A.3 in the appendix, has to be performed. Figure 5.3 shows the
exchange of data between OpenFOAM and Cantera through the interface canteraFoamModel. OpenFOAM solves the implemented transport equations for energy, species,
velocity and pressure, and passes the thermodynamic data in vectors of dimension number of control volumes N to the interface. The interface then constructs for each single
control volume an independent reactor. For every reactor Cantera calculates cell depen-
33
Fig. 5.2.: Initialization sequence of the coupling library in a solver during runtime (order: top left first,
bottom right last).
dent data like the chemical reaction rates rk or transport properties like , , Dkj , DkT .
Cantera then returns per reactor vectors of dimension number of species K for rk or DkT ,
simple scalars for or , and matrices of dimension (K K) for Dkj . The interface
collects the values of every cell and constructs vectors or matrices with cell values for
every control volume and returns matrices of dimension (N K) for the reaction rates
or of dimension (N K K) for the ordinary diffusion coefficients Dkj and vectors of
dimension (N ) for or .
5.3. IMPLEMENTATION OF THE SOLVER
34
Fig. 5.3.: Data exchange between OpenFOAM and Cantera through the interface canteraFoamModel.
5.3. Implementation of the Solver

The new solver is based on the OpenFOAM solver reactingFoam. Transport equations
for species mass and energy based on equation (3.7) and (3.15) replace the approach
in the original solver discussed in the next subsection. To use functions from Cantera
like multicomponent mass diffusion coefficients or reaction rates, the interface presented
above has to be loaded dynamically during runtime or linked while compiling. Note,
that accurate expressions for the continuity equation 3.9 as well as for the momentum
equation 3.10 are already implemented in the original solver and therefore not altered.
5.3.1. OpenFOAMs Approach for Transport Equations

The original transport equations for species and sensible enthalpy are
~ (Yk~v ) +
~ (Y
~ k ) + rk
(Yk ) =
t
~ (hs~v ) +
~ ( h
~ s ) + Dp + ~qreaction
(hs ) =
t
cp
Dt
(5.1)
(5.2)
One can see, that equation (5.1) is based on the assumption that the ratio of momentum
diffusivity and mass diffusivity is unity. Thus, the dimensionless Schmidt number Sc
yields
Sc =
= 1 D = .
(5.3)
D
Therefore, the mass diffusivity D has to be the same for all species and the species
diffusion process is entirely controlled by the mixture viscosity, which is calculated from
35
Sutherlands law. Additionally, cross diffusion process based on pressure and temperature gradients are neglected.
In equation (5.2) a unity Lewis number Le as well as the disregards from above are
assumed:
Le =
=1
= D
(5.4)
Dcp
cp
It has been shown by [12], that equations (5.1) and (5.2) yield large errors while simulating chemically reacting flows.
5.3.2. Implementation of Species Mass Equations

In nearly isobar combustion the flux caused by the pressure gradient in equation (3.20)
can be neglected. External forces occur for example in the diffusion of electrically charged
particles and will be neglected, too. With these simplifications and the relation
N
X
M ~
~ j Yj M
~ j = M Y
Yi ,
X
Mj
M j i=1 M i
(5.5)
equation (3.20) can be rearranged to

~jk = M k
M
N
X
~ j Mk
Dkj Y
M
j=1
N
X
Dkj Yj
j=1
N
X
M ~
1~
Yi DkT T
,
T
i=1 M i
(5.6)
to obtain the mass diffusive flux in terms of mass fractions.

For numerical stability, including a discretized second derivation or diffusion term into
the solution matrix is essential [7]. A direct implementation of the mass diffusive flux
equation (5.6) is therefore not worthwhile, because of the missing diffusion terms which
can be discretized implicit by OpenFOAMs class fvm. Through the use of the, in
section 3.3 already mentioned, characteristics of the ordinary diffusion coefficient matrix,
equation (5.6) is rearranged to
~jk =
N
X
N X
N
H
X
X
Mk
Mk
~ k DT 1 T
~
~
~ ,
Dkj Yj
Dkj Yj Yi
Dkj Yj Y
k
T
M
M
i
j6=k i6=k
j6=k
j6=k
(5.7)
so that the third term on the right hand side can be discretized implicit by OpenFOAM.
Finally, by inserting the expression for the mass flux above into equation (3.7) the
transport equation for species k reads
N
X
~ (Yk~v )
~ (Dkj Yj Y
~ k) =
(Yk ) +
t
j6=k
N X
N
X
Mk
~
~
~
~ M k Dkj Yj Y
~ i +
~ DT 1 T
Dkj Yj +
+ rk ,

k
T
M
M
i
j6=k
j6=k i6=k
(5.8)
N
X
36
in which all terms discretized implicit by OpenFOAMs class fvm stand on the left hand
side and all terms included and calculated as explicit source terms by OpenFOAMs class
fvc stand on the right hand side. Listing A.1 in the appendix shows the implementation
of equation (5.8) in OpenFOAM.
5.3.3. Implementation of Sensible Enthalpy Equation

As mentioned above a diffusive flux term has to be included into the solution matrix.
Therefore, it is necessary to express Fouriers Law in terms of sensible enthalpy. Starting
with an expression of the enthalpy gradient for a multicomponent gas mixture given
by [21],
~ = cp T
~ +
h
N
X
k=1
~ k,
hk Y
(5.9)
one can, by express the specific enthalpy with the sum of sensible and chemical enthalpy
(3.13), modify Fouriers law to
~
T
=
N
X
~
~
h
hk Yk
cp
k=1 cp
N
N
N
X
X
X
~
~
~ k
~
hs +
hc,k Yk
hs,k Y
hc,k Yk
=
cp
cp
c
c
p
p
k=1
k=1
k=1
N
X
~
~ k.
=
hs
hs,k Y
cp
c
p
k=1
(5.10)
Now, the first term on the right hand side can be discretized implicit by OpenFOAM.
The diffusive species flux from equation (3.31) in terms of sensible Enthalpy gives
N
X
k=1
~jk hk =
N
X
~jk (hc,k + hs,k ) .
(5.11)
k=1
Combining equations (3.15), (3.31), (5.10), (5.11) and neglecting viscous dissipation
~ v and the Soret effect one obtains
: ~
!
~ (hs~v )
~ h
~ s =
(hs ) +
t
cp
N
N
N
X
X
Dp X
~
~
~
~
hs,k Yk
(hs,k jk )
hc,k rk ,
Dt k=1
cp
k=1
k=1
(5.12)
in which all terms discretized implicitly by OpenFOAMs class fvm stand on the left hand
side and all terms included and calculated as explicit source terms by OpenFOAMs
class fvc stand on the right hand side. The diffusive species flux ~jk is calculated by
interpolation from the species equation noted in listing A.1 as J[k]. Listing A.2 in the
appendix shows the implementation of equation 5.12 in OpenFOAM.
6. Validation of the Solver

In this chapter, the implemented solver introduced in chapter 5 is validated with focus
on the transport properties and reaction rates. The reference solution is obtained from
CHEMKIN/PREMIX with equal physical conditions. In addition, the effect of the grid
resolution onto the solution and the chemical time scales are investigated.
6.1. Case Description

Figure 6.1 shows the physical domain of the considered case. Here, the fresh gas mixture
flows from the left into the domain. Due to thermal conduction, the mixture is preheated
(dark blue) until it ignites (light blue). In the following, the outburn proceeds (red zone)
and the adiabatic flame temperature is reached.
1.0
YCH4 ,in
YO2 ,in
YN2 ,in
Tin [K]
p [bar]
vx,in [m/s]
0.055
0.2202
0.7248
300
1.0
0.4
Tab. 6.1.: Physical conditions of the simulation with OpenFOAM.
The mixture is stoichiometric with = 1 and the mass fractions at the inlet are YO2 ,in =
0.2202, YCH4 ,in = 0.055 and YN2 ,in = 0.7248. The fluids temperature at the inlet is
Tin = 300 K with an inflow velocity vx,in = 0.4 m/s. The pressure in the whole domain
as well as outside is constantly p = 1 bar. The physical conditions (shown in table 6.1)
are chosen to be the same as in the reference calculation.
6.2. NUMERICAL CONDITIONS
38
Fig. 6.1.: Physical domain of the laminar, premixed flame.
6.2. Numerical conditions

Table 6.2 shows the numerical setup of the various simulations performed with the new
solver in OpenFOAM. The different grids consist of 600 and 3000 equidistant rectangular
cells built with the blockMesh tool from OpenFOAM. The grid resolutions are x =
5 105 m and x = 1 105 m for each cell. The time steps are set to t = 5
107 s on the coarse grid and t = 1 107 s on the fine mesh. For all simulation,
including the reference calculation, the mechanism by [38], which contains 17 species
and 58 reactions (see Appendix B), is used to calculate the chemical kinetics. For
thermodynamic and molecular properties the same input files are used in OpenFOAM
and CHEMKIN/PREMIX.
Coarse
Fine
Domain [m]
Cells
t [s]
x [m]
0.03
0.03
600
3000
5 107
1 107
5 105
1 105
Tab. 6.2.: Numerical setups.
Figure 6.2 shows the numerical domain. At the inlet, fixed values are chosen as boundary condition for all physical properties. The boundary conditions at the outlet differ
because of different requirements. For mass fraction and temperature a simple zero gradient boundary condition is adequate. To control possible backflow at the outlet a mixed
boundary condition, the inletOutlet condition, is needed for the velocity. It switches the
boundary condition between a fixed value and zero gradient depending on direction of
the velocity. Since the numerical domain is just a clipping of the infinity large physical
domain, the front, back, top and bottom sides are defined empty. Additionally, the numerical domain is only divided into cells, marked by the dashed lines, in the x-direction.
Therefore, only a solution in this direction is approached leading to a one dimensional
case.
6.2. NUMERICAL CONDITIONS
39
Fig. 6.2.: Numerical domain of the laminar, premixed flame.
In the calculations the implicit Euler method is used to discretize the time integral.
Like mentioned in chapter 4 the implicit Euler method is of order O(1)t . The spatial
discretization is accomplished via Gaussian integration and using linear interpolating
schemes of order O(2)x . The discretization schemes used in OpenFOAM for the different
terms in the transport equations are listed in table 6.3.
Mathematical Term
/t
()
()
Expression in OpenFOAM
Discretization
Order
ddt
grad
div
laplacian
interpolation
Euler implicit
Gauss linear
Gauss linear/limitedLinear 1
Gauss linear corrected
linear
O(1)t
O(2)x
O(2)x
O(2)x
O(2)x
Tab. 6.3.: In OpenFOAM used discretization schemes for the mathematical terms.
To obtain similar conditions like in CHEMKIN/PREMIX, no transport equation for velocity and pressure is solved. The remaining transport equations are solved using different iterative solver and preconditioners, namely for density the Preconditioned conjugate
gradient (PCG) solver with Diagonal incomplete-Cholesky (DIC) preconditioning and
for sensible enthalpy and species mass fractions the Preconditioned biconjugate gradient
(PBiCG) solver with Diagonal incomplete-LU (DILU) preconditioning.
For the initialization of the simulation a homogeneous flow field with conditions like at
the inlet is shortly ignited. Therefore, an spatial limited source term is included in the
transport equation of enthalpy with an ignition strength of 5 109 J/m3 s for 0.0002 s.
In the ignited cells two flame fronts develop. The downstream traveling flame front
leaves the domain and is not of further interest.
6.3. COMPARISON WITH CHEMKIN/PREMIX
40
6.3. Comparison with CHEMKIN/PREMIX

To validate the OpenFOAM solver a simulation with 3000 cells is performed. The focus
of the validation are the newly implemented equations for species and sensible enthalpy
as well as the transport coefficients and reaction rates obtained from Cantera. The
calculations are compared to results obtained from CHEMKIN/PREMIX with similar
conditions. CHEMKIN/PREMIX is a highly specialized computer code for premixed
laminar flames in 1D [20]. The code presumes isobar conditions and a steady state flame.
R
To get similar conditions in OpenFOAM
no pressure correction is performed and the
velocity field is calculated from
0 v0
p
m
= 0 v0 = v = const. v =
and =
.
(6.1)
A
T Rs
in which Rs = R/M is the specific gas constant of the mixture. For this reason, the
pressure is forced to remain constant at all times in the complete domain. Furthermore,
in CHEMKIN/PREMIX the inlet mass flux m
is changed during the simulation in such
a way to hold the flame front steady. Due to missing moving mesh features implemented
into the new solver it is therefore not possible to completely prevent a motion of the
flame front. CHEMKIN/PREMIX is capable of dynamic mesh refinement and applies
for the current case a total number of only 185 grid points with the minimal distance of
xmin 3 106 m.
6.3.1. Temperature and Species Mass Fractions

In figure 6.3 the temperature profile from 3000 Cells (solid) is compared to the reference calculation (circles). The solutions are nearly the same, a part from two notable
disagreements at the beginning and the end of the inner reaction zone with maximum
variations of about 17 K. One can see, that the grid used by CHEMKIN/PREMIX at
the beginning of the reaction zone is too low against the fine mesh from OpenFOAM.
Additionally, the discretization scheme used by CHEMKIN/PREMIX is the backwarddifference scheme. It is of order one and introduces, similar to the upwind interpolation
(UDS) mentioned before, artificial diffusion. Therefore, CHEMKIN/PREMIX spreads
out the solution on a coarse mesh [7, 20] and can not resolve the large gradient as accurately as OpenFOAM.
Case
YCOmax
YCH3 max
YOHmax
YH2 max
sL [m/s]
CHEMKIN/PREMIX
OpenFOAM
0.051967
0.05240
0.00277
0.00274
0.00510
0.00497
0.00163
0.00163
0.373
0.376
0.83 %
1.08 %
2.54 %
0.00 %
0.8 %
Deviationa
a
Deviation based on CHEMKIN/PREMIX calculations

Tab. 6.4.: Comparison of distinctive results obtained with OpenFOAM and CHEMKIN/PREMIX.
In figure 6.3 the species mass fractions YCH4 , YO2 , YCO2 , YCO and YH2 O (solid) are compared to the reference calculation (circles). As one can see, the solution obtained from
6.3. COMPARISON WITH CHEMKIN/PREMIX
41
OpenFOAM is very similar to the reference solution. The profiles of intermediate species
mass fractions YOH , YO , YH2 , YCH3 and YCH2 O as well as the minimal species mass fractions YCH2 , YH , YH2 O2 , YHCO and YHO2 confirm that conclusion, although there are tiny
discrepancies (see table 6.4) between the reference calculation and OpenFOAM which
have a maximum variation at YOH of about 1.3 104 . Again, the different grid resolution between OpenFOAM and CHEMKIN/PREMIX in these regions may be the reason. By using equation (2.4) the laminar burning velocity calculates to sL = 0.376 m/s.
As one can see in table 6.4 this is very similar to the value obtained from CHEMKIN/PREMIX. Additionally, a comparison with experimental results, e.g. [39], confirms
a very good agreement of OpenFOAMs value, too.
Axial distance x [cm]
Temperature T [K]
1.05
1.1
1.15
1.2
OpenFOAM
CHEMKIN
1.05
1.1
1.15
1.2
0.2
O2
1500
CO2
1000
CH4
0.1
H2 O
CO
Mass fraction Y [-]
1
2000
500
0
104
103
4
OH
CH3
HO2
O
2
HCO
H2 O2
H2
CH2
CH2 O
0
1
1.05
1.1
1.15
Mass fraction Y [-]
Mass fraction Y [-]
0
1.2
1.05
1.1
1.15
1.2
Fig. 6.3.: Temperature and species mass fraction profiles obtained from OpenFOAM (lines) and from
CHEMKIN/PREMIX (circles).
6.4. INFLUENCE OF THE GRID RESOLUTION
42
6.4. Influence of the Grid Resolution

The calculation time in a numerical simulation highly scales with the number of cells [7].
Therefore, knowledge about the influence of the grid resolution and the time step onto
the solution is essential. In this section the validated grid from the previous section with
3000 cells is compared to the coarse grid with 600 cells. The CFL number was kept the
same leading to time steps up to t = 5 107 s.
In figure 6.4 the temperature profile obtained with 600 Cells (dashed) is compared to the
one obtained with 3000 Cells (cross). The two noticeable disagreements at the beginning
and the end of the inner reaction zone show maximum variations of about 15 K.
The heat release calculated from equation 3.16 is shown at the bottom of figure 6.4. One
can see, that the reaction zone is slightly displaced and thicker in the solution obtained
with the coarse grid compared to the one obtained with the finer grid. This is consistent
with the mass fraction profiles mentioned further down. Moreover, one can identify, that
most of the heat is released at higher temperatures. Upstream from the flame front, no
chemical reactions occur and the temperature rises due to thermal conduction.
Case
YCOmax
YCH3 max
YOHmax
YH2 max
sL [m/s]
600 Cells
3000 Cells
0.05080
0.05240
0.00248
0.00274
0.00497
0.00497
0.00162
0.00163
0.385
0.376
Deviationa
3.14 %
10.0 %
0.0 %
0.61 %
2.33 %
Deviation based on 600 Cells calculations

Tab. 6.5.: Comparison of distinctive results obtained with 600 Cells and 3000 Cells.
Figure 6.4 shows the major species mass fractions YCH4 , YO2 , YCO2 , YCO and YH2 O for 600
Cells (dashed) compared to 3000 Cells (line). One can see, that the solution obtained
from the coarser grid is very similar to the one obtained with the finer grid. The
profiles of the minor species mass fractions YOH , YO , YH2 , YCH3 and YCH2 O as well as the
minimal species mass fractions YCH2 , YH , YH2 O2 , YHCO and YHO2 in figure 6.4 confirm
this conclusion. From the differences between the mass fractions profiles of 600 Cells
and 3000 Cells one can see, that the peaks of the intermediate species differ only slightly
in magnitude and spatial position. The laminar burning velocity sL is nearly the same
for all grids.
6.4. INFLUENCE OF THE GRID RESOLUTION
43
Temperature T [K]
2000
1.05
1.1
1.15
1.2
1.05
Coarse
Fine
1.1
1.15
1.2
0.2
O2
1500
CO2
1000
0.1
H2 O
CH4
CO
Mass fraction Y [-]
500
Mass fraction Y [-]
104
OH
CH3
HO2
O
HCO
CH2 O
H2
H2 O2
CH2
Mass fraction Y [-]
103
0
1
1.05
1.1
1.15
1.2
1.05
1.1
1.15
1.2
Heat release rate q [J/m3 s]
109
4
3
2
1
0
1.05
1.1
1.15
1.2

Fig. 6.4.: Profiles obtained from OpenFOAM with 600 Cells (dashed line) and 3000 Cells (line).
6.5. SPECIES PRODUCTION RATES
44
6.5. Species Production Rates

In figure 6.5 the computed reaction rates for the species are shown. By comparing the
reaction rates with the mass fraction profiles from the section before, one can see, that
species with high diffusion coefficients (e.g. H radicals) diffuse upstream against the flow
into the preheat region. Here, HO2 begins to form which is later, at higher temperatures,
transformed into OH radicals. The OH radicals become relatively abundant versus the
O or H radicals to latter stages of the reaction zone. Throughout the reaction zone, CO
molecules build up until they are rapidly consumed by the abundant OH radicals. This
is responsible for the major heat release of the combustion process and also marks the
end of the reaction zone.
Reaction rate r [kg/m3 s]
1
200
1.05
1.1
1.15
1.2
1.05
1.1
H2 O
OH
CH3
CO2
1.15
0
CO
200
H2
CH4
CH2 O
O2
400
1.2
40
20
0
20
40
5
H
CH2
HO2
H2O2
HCO
1.05
1.1
1.15
1.2

Fig. 6.5.: Reaction rates obtained from OpenFOAM.
6.6. Computation of Chemical Time Scales

To determine the necessary time step for the numerical simulation of the turbulent
premixed flame the chemical time scales were calculated from equation (3.39). Figure
6.6. COMPUTATION OF CHEMICAL TIME SCALES
45
6.6 shows the time scales for T , O2 and OH against the temperature.
Chemical time scale [s]
100
101
OH
102
T
O2
103
104
105
300
600
900
1200
1500
1800
2100
Temperature T [K]
Fig. 6.6.: Calculated chemical time scales for T , O2 and OH.
One can see, that in an region, which corresponds to the main reaction zone, the time
scales are almost constant. For intermediate species the time scales show characteristic
peaks which correspond to formation and oxidation reactions. The minimum of the
calculated time scales determines the chemical time scale of the simulation. Here, the
time scale ranges between 105 s and 104 s. The time scale of temperature, which is
the most important for combustion modeling, is 2 104 s in this case.
7. Direct Numerical Simulation of a

Turbulent Premixed Flame in 3D
In this chapter, the simulation of a turbulent premixed flame in 3D is presented. The
main purpose is to show that a direct numerical simulation of turbulent reacting flow
using the validated solver is in principle possible. Therefore, the results are only briefly
discussed with main focus on the flame topology. In addition, since the simulations are
carried out on a Cray XE6 supercomputer [40], the parallel performance of the solver in
conjunction with the supercomputer is investigated.
7.1. Numerical Setup

Based on the results from the last chapter, the physical and numerical conditions are
chosen, see table 7.1. It has been shown before, that a numerical grid with x = y =
z = 5 105 m and a time step t = 5 107 s is a good compromise between the
desired accuracy and calculation time. This conclusion is later approved by the calculated Kolmogorov scales. Figure 7.1 shows the computational domain for the three
dimensional case. The dimension is (1.5 0.75 0.75) cm3 , which is thus half the length
to the one dimensional case, with 300 150 150 = 6.75 106 equidistant cells of size
x = y = z = 5 105 m. The physical conditions listed in table 7.1 are similar
to the one dimensional case in chapter 6 except for the inlet velocity field which has
0
= 2 m/s. The Reynolds number
the mean value ~vin = 4.0 m/s and fluctuates with ~vin
based on integral dimensions at the inlet is Re = 2000. The main differences, besides
the turbulence, are the three dimensional resolution of the shorter physical domain and
the needed solution of three additional transport equations for velocity and one equation
for pressure-correction. For chemical kinetics and thermodynamic properties as well as
molecular quantities the same input files are used as before. The total simulation time
is approximately 5 ms, which leads to a wall-clock time of approximately 48 hours on
the Cray XE6 with 2048 CPUs.
YCH4 ,in
YO2 ,in
YN2 ,in
Tin [K]
p [bar]
~vin [m/s]
0
~vin
[m/s]
Re
1.0
0.055
0.2202
0.7248
300
1.0
4.0
2.0
2000
Domain [m]
Cells
t [s]
x [m]
Discret. Order
(0.015 0.0075 0.0075)
6.75 106
5 107 s
5 105
O(2)x O(1)t
Tab. 7.1.: Physical and numerical conditions of the turbulent, premixed flame in 3D.
7.1. NUMERICAL SETUP
47
Since the flow is three dimensional, solutions in directions outside the domain (in addition to the in- and outlet) are needed. To ensure, that no mass is lost, the front, back,
top and bottom faces of the domain are defined as symmetry planes. The boundary
conditions at the in- and outlet for temperature T and species mass fractions Yi are
the same as for the one dimensional case. To avoid reflections of pressure waves at
the in- and outlet, non-reflective boundary conditions (NRBC) are used for pressure.
To generate a turbulent flow at the inlet a turbulent inflow generator by [41] is used
for the velocity field. The boundary condition is based on digital filtering of random
data to provide temporally and spatially correlated velocities. The mean flow velocity
is ~vin = 4.0 m/s with a turbulence level of 50 %. The turbulent length scale lT and the
turbulent time scale T are set to lT = 1 mm and T = 0.5 ms.
Fig. 7.1.: Numerical domain of the turbulent, premixed flame.
7.2. INITIAL CONDITIONS
48
7.2. Initial conditions

The initial conditions have been produced by a DNS on the same computational grid,
however, without chemical reactions. Then, the flat flame front with suitable temperature and mass fraction fields from the one dimensional case has been superimposed into
the flow. Figure 7.2 shows the initial velocity and the temperature fields right before
the simulation starts. One can see, that the flame front (marked with red colors) is still
planar since the passed time is short.
Fig. 7.2.: Initial velocity and temperature fields of the premixed turbulent flame.
Figure 7.3 shows the Kolmogorov scales. They are calculated from the produced initial
conditions by using equations (2.6). One can see, that the Kolmogorov length scale lL
varies between 2 105 m and 2.5 104 m. Since the small values occur mostly near
the inlet the used grid spacing can be expected to resolve the smallest eddies near the
flame front sufficiently. The Kolmogorov time scale K varies between 2 105 s and
2.83 103 s. Since K is relatively great, the used time step is limited by the Courant
number due to numerical stability reasons.
Fig. 7.3.: From the initial conditions calculated Kolmogorov scales.
7.3. TOPOLOGICAL RESULTS
49
7.3. Topological Results

The main purpose is to show, that a direct numerical simulation of turbulent reacting
flow using the validated solver is possible. Therefore, the results are given only briefly
with main focus on the flame topology. Figure 7.4 shows an isosurface of the heat
release rate, which can be taken as the position of the flame front. The heat release
rate is calculated by using equation (3.16). A corrugated surface, which is caused by
the turbulent flow shown in the back can be identified. On the bottom of figure 7.4 the
eddy dissipation rate is shown. The eddy dissipation rate is very high at the flame front
and near zero downstream, since no eddies are present anymore.
Fig. 7.4.: Isosurface of YCH4 . Bottom with eddy dissipation rate and backside with vorticity.
The ratio of mean flame surface and turbulent flame surface calculates to AAT = 2.33.
Additionally, by using equation 2.10 the turbulent burning velocity calculates to sT =
0.704 m/s, which is less than two times the laminar burning velocity calculated for the
laminar flame in the chapter before. The difference in the laminar flame velocities is
caused by the curved and stretched turbulent flame.
Figure 7.5 shows contour plots of the temperature and velocity field. The flame front
is marked by iso-contours of the heat release rate. The structure of the turbulence is
visible at the inlet and the flow is accelerated by passing the flame front. The temperature increases over the flame front from 300 K up to 2200 K. The velocity gradients
in the burnt mixture are clearly smaller caused by high viscosity. A distinct feature
of the flame topology is the presence of an instantaneous wrinkled flame front. In areas in which the flame front wrinkles with a convex geometry into the burnt gas, the
50
Fig. 7.5.: Two dimensional slices of the temperature field and the vorticity field with heat release rate
of the turbulent premixed flame.
heat release is higher. Additionally, due to the turbulent flow, the flame front varies in
thickness. In this case, small eddies penetrate the flame front and enhance the diffusion
process, leading to a thickened flame front. In areas where the cold flow is slow the flame
front even propagates into the unburnt gas. It is evident that the vorticity is greatly
damped in the burned gas side because the gas viscosity increases with temperature.
Additionally, in figure 7.6 the mass fractions of CH4 , O2 , H2 and OH as well as the
corresponding reaction rates are shown. One can see, how the reaction zone varies in
thickness and the strain of the flame front influences the reaction rates.
Fig. 7.6.: Different mass fractions and reaction rates.
51
7.4. PARALLEL PERFORMANCE OF THE CRAY XE6 (HERMIT)
52
7.4. Parallel Performance of the Cray XE6 (HERMIT)

To perform direct numerical simulations of chemically reacting flows, at least three velocity equations, one mass equation, one enthalpy equation and K 1 species equations
have to be solved iteratively. The number even increases if pressure-velocity coupling
schemes are used [7]. In addition, calculations for the chemical reaction rates or transport properties have to be performed at least once per time step and cell. In particular,
the calculation time of multicomponent mass diffusion coefficients, shown in section 3.3,
requires the inversion of a K K matrix [16]. Further more, if turbulent flow is considered, it can be shown, that the number of grid points in each direction must be at least
l/lK [7].
Due to the small time steps and the high grid resolution, direct numerical simulations of
chemically reacting flows have to be performed in parallel on high performance computing clusters. Since the architecture of each cluster is very unique and an efficient usage
of existing resources is essential one has to know about the speed-up particularly of high
performance clusters.
The numerical simulation of the turbulent premixed flame in 3D presented before was
performed on the Cray XE6 (HERMIT) operated by the High Performance Computing
Center Stuttgart (HLRS). It is a massively parallel supercomputer with 3552 compute
nodes. Each compute node has two sockets with each 16 AMD Interlagos processors.
Each XE6 node has either 32 or 64 GB of DDR3 SDRAM memory. The node-node
interconnection is based on CRAY Gemini. The XE6 is operated with the Cray Linux
Environment in version 3, including SUSE Linux Enterprise Server and Crays Compute
Node Linux [40, 42].
For the current investigation, several calculations have been carried out to obtain information about parallel efficiency and scalability. The new OpenFOAM solver (see
chapter 5) has been used and the conditions for all calculations where the same except
for the number of cores. The mesh consists of about 2 106 cells.
7.4. PARALLEL PERFORMANCE OF THE CRAY XE6 (HERMIT)
53
Figure 7.7 shows the intranode scale-up (left) and internode scale-up (right) calculated
with Sn = ts /tn [7]. Here ts is the execution time with the reference processor number
and tn is the execution time using n processors. The reference processor number for the
intranode scale-up is 1 CPU and for the internode scale-up equal the number of CPUs
per node (32 CPUs). One can see, that the intranode speedup is far from the linear
speed up marked as a solid line. Thus, it is not worthwhile using the Cray XE6 with
only a few CPUs on one node. On the contrary, the internode performance scales very
well up to 2048 processors. Between 64 processors and 1024 processors the scale-up
is even superlinear1 [7]. The scale-up indicates a very good scalability of Cray XE6 in
conjunction with the solver even for 6.75 106 Cells. Therefore, in the DNS of the
turbulent flame in the section before 2048 CPUs have been used.
ideal
ideal
24
Internode scale-up
Intranode scale-up
25
23
22
1
24
22
20
20
20
21
22
23
24
No. of CPUs
25 = 32
25
27
29
211 = 2048
No. of CPUs
Fig. 7.7.: Scale-up on Cray XE6 with a grid of 2 106 cells.
Superlinear is a speedup of more than p if p processors are used. It occurs rarely in parallel computing.
[43] showed with a simple example how a superlinear speedup is achieved.
8. Summary and Perspective

A solver for direct numerical simulations of chemically reacting flows has been implemented and validated for the free CFD framework OpenFOAM 2.0.1. The solver takes
detailed multicomponent diffusive fluxes into account. The molecular transport coefficients as well as the chemical kinetics and thermodynamic quantities were obtained via
a coupling interface from Cantera 1.8. The solver as well as the coupling library were
successfully compiled with an MPI implementation on the Cray XE6 for the parallelized
computation of DNS of turbulent reacting flows.
To validate the solver, numerical calculations of a laminar premixed flat flame have been
carried out. The results were subsequently compared to reference data obtained from
the CHEMKIN package. The agreement with respect to the computed flame profiles for
T and Y was very good. The calculated flame speed showed reasonable results with experimental data. Additional calculations were performed on a coarse grid to investigate
the influence of the gird resolution on the computed flame structure. The results were
shown to depend very little on whether the coarse or fine mesh had been used. Subsequently, direct numerical simulation of a turbulent premixed flame in three dimensions
have been performed on the Cray XE6 supercomputer with 2048 CPUs. The interaction
between turbulence and the flame front has been identified. Thereby, the capability of
the solver for complex 3D turbulent flows with chemical reactions as well as the very
good scale-up on the Cray XE6 was shown.
For future work, the validated DNS solver will be used for further investigations of the
interplay between turbulence and chemical reactions as well as characteristic values like
strain rate, the curvature or the turbulent burning velocity. The influence of selective
diffusion mechanisms (like e.g. the Soret effect) on the numerical solution of turbulent
premixed flames can be studied. Furthermore, essential statistical information, which is
needed to develop and validate models used in future RANS and LES simulation, can
be provided by the DNS solver.
Thanks to the usage of C++ programming principles it is possible to extend the solver by
including new numerical or computational utilities in order to study additional physical
phenomena. The investigation of details like the radiative heat transport or the Dufour
effect can further enrich such a study. In order to reduce the CPU time in numerical
simulations a mixture average approach for the diffusive fluxes can be included as well as
an adaptive mesh refinement based on characteristic gradients. Equations for the stretch
rate or flame velocity and models for soot or N Ox formation can also be incorporated
and validated for further examination of turbulent premixed flames.
A. Code Fragments for OpenFOAM in

C++3
The following listings show code fragments from the implemented OpenFOAM solver.
Listing A.1 shows the implementation of the transport equation (5.8) for species and
listing A.2 shows the implementation of the transport equation (5.12) for sensible enthalpy. In addition, listing A.3 shows the implementation of the initialization of the
coupling interface and a downcast to obtain access to derived class functions. Note, that
for the sake of clarity, the codes are simplified and will not compile without modifications.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
f v S c a l a r M a t r i x YEqn
(
fvm : : ddt ( rho , Y[ k ] )
+ mvConvection>fvmDiv ( phi , Y[ k ] )
==
c h e m i s t r y .RR( k )
);
f o r A l l (Y, j )
{
i f ( j !=k )
{
f o r A l l (Y, i )
{
i f ( i !=k )
{
YEqn = f v c : : l a p l a c i a n (W[ k ] /W[ i ] rho c o m p o s i t i o n .D( k , j ) Y[ j ] ,
Y[ i ] , " l a p l a c i a n (D,Y[ k ] ) " ) ;
J [ k ] = f v c : : i n t e r p o l a t e (W[ k ] /W[ i ] rho c o m p o s i t i o n .D( k , j ) Y[ j ] )
( f v c : : i n t e r p o l a t e ( f v c : : grad (Y[ i ] ) ) & mesh . S f ( ) ) ;
}
}
21
YEqn = fvm : : l a p l a c i a n ( rho c o m p o s i t i o n .D( k , j ) Y[ j ] , Y[ k ] ,

" l a p l a c i a n (D, Yi ) " ) ;
J [ k ] = f v c : : i n t e r p o l a t e ( rho c o m p o s i t i o n .D( k , j ) Y[ j ] )
( f v c : : i n t e r p o l a t e ( f v c : : grad (Y[ k ] ) ) & mesh . S f ( ) ) ;
22
23
24
25
26
27
28
29
30
YEqn += f v c : : l a p l a c i a n ( rho W[ k ] UnitMole Wm c o m p o s i t i o n .D( k , j ) , Y[ j ] ,

" l a p l a c i a n (D,Y[ k ] ) " ) ;
J [ k ] += f v c : : i n t e r p o l a t e ( rho W[ k ] UnitMole Wm c o m p o s i t i o n .D( k , j ) )
( f v c : : i n t e r p o l a t e ( f v c : : grad (Y[ j ] ) ) & mesh . S f ( ) ) ;
YEqn = f v c : : l a p l a c i a n ( c o m p o s i t i o n .DT( k ) /T, T, " l a p l a c i a n (D,Y[ k ] ) " ) ;
56
31
32
33
34
J [ k ] = f v c : : i n t e r p o l a t e ( c o m p o s i t i o n .DT( k ) /T)
( f v c : : i n t e r p o l a t e ( f v c : : grad (T) ) & mesh . S f ( ) ) ;
YEqn . r e l a x ( ) ;
s o l v e (YEqn , mesh . s o l v e r ( "Y[ k ] " ) ) ;
Listing A.1: Implementation of the species equation 5.8 in OpenFOAM.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
f v S c a l a r M a t r i x hsEqn
(
fvm : : ddt ( rho , hs )
+ mvConvection>fvmDiv ( phi , hs )
fvm : : l a p l a c i a n ( t u r b u l e n c e >a l p h a E f f ( ) , hs )
==
DpDt
+ chemistrySh
);
f o r A l l (Y, k )
{
hsEqn = f v c : : l a p l a c i a n ( thermo . a l p h a ( ) h s i [ k ] , Y[ k ] ) ;
hsEqn = f v c : : d i v ( J [ k ] , h s i [ k ] , " d i v ( J i , h s i ) " ) ;
}
hsEqn . r e l a x ( ) ;
hsEqn . s o l v e ( ) ;
Listing A.2: Implementation of the energy equation 5.12 in OpenFOAM.
1
2
3
4
5
autoPtr<psiChemistryModel> pChemistry ( psiChemistryModel : : New( mesh ) ) ;

psiChemistryModel& c h e m i s t r y = pChemistry ( ) ;
hsCombustionThermo& thermo = c h e m i s t r y . thermo ( ) ;
c a n t e r a M i x t u r e <canteraThermo>& c o m p o s i t i o n =
dynamic_cast<c a n t e r a M i x t u r e <canteraThermo>&>(thermo . c o m p o s i t i o n ( ) ) ;
Listing A.3: Initialization during run time and performing a downcast to obtain access to derived class
functions in OpenFOAM.
B. Reaction Mechanism
Listing B.1 shows the chemical reaction mechanism by Kee et al [38], which was used in
the numerical simulations of the 1D flame as well as in the direct numerical simulations.
It consists of 17 species and 58 reactions.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
! Kee , Gracr , Smooke , M i l l e r : SAND Report 85 8240 , 1985 Livermore ,

C a l i f o r n i a , USA
! 5 8 R e a c t i o n s , 17 S p e c i e s
ELEMENTS
H O C N
END
SPECIES
H2 O2 H2O H O OH HO2 CO
CH4 CH3 CH2 CH CH2O HCO
END
REACTIONS
! U n i t s i n cm / s / K /
CH3+H+M=CH4+M
CH4+O2=CH3+HO2
CH4+H=CH3+H2
CH4+O=CH3+OH
CH4+OH=CH3+H2O
CH3+O=CH2O+H
CH3+OH=CH2O+H2
CH3+OH=CH2+H2O
CH3+H=CH2+H2
CH2+H=CH+H2
CH2+OH=CH2O+H
CH2+OH=CH+H2O
CH+O2=HCO+O
CH+O=CO+H
CH+OH=HCO+H
CH+CO2=HCO+CO
CH2+CO2=CH2O+CO
CH2+O=CO+H+H
CH2+O=CO+H2
CH2+O2=CO2+H+H
CH2+O2=CH2O+O
CH2+O2=CO2+H2
CH2+O2=CO+H2O
CH2+O2=CO+OH+H
CH2+O2=HCO+OH
CH2O+OH=HCO+H2O
CH2O+H=HCO+H2
CH2O+M=HCO+H+M
CH2O+O=HCO+OH
CO2 H2O2
N2
c a l / mol
8 . 0 0 E26
3.000
7 . 9 0 E13
0.000
2 . 2 0 E04
3.000
1 . 6 0 E06
2.360
1 . 6 0 E06
2.100
6 . 8 0 E13
0.000
1 . 0 0 E12
0.000
1 . 5 0 E13
0.000
9 . 0 0 E13
0.000
1 . 4 0 E19
2.000
2 . 5 0 E13
0.000
4 . 5 0 E13
0.000
3 . 3 0 E13
0.000
5 . 7 0 E13
0.000
3 . 0 0 E13
0.000
3 . 4 0 E12
0.000
1 . 1 0 E11
0.000
3 . 0 0 E13
0.000
5 . 0 0 E13
0.000
1 . 6 0 E12
0.000
5 . 0 0 E13
0.000
6 . 9 0 E11
0.000
1 . 9 0 E10
0.000
8 . 6 0 E10
0.000
4 . 3 0 E10
0.000
3 . 4 3 E09
1.180
2 . 1 9 E08
1.770
3 . 3 1 E16
0.000
1 . 8 1 E13
0.000
0.
56000.
8750.
7400.
2460.
0.
0.
5000.
15100.
0.
0.
3000.
0.
0.
0.
690.
1000.
0.
0.
1000.
9000.
500.
1000.
500.
500.
447.
3000.
81000.
3082.
58
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
HCO+OH=CO+H2O
5 . 0 0 E12
0.000
0.
HCO+M=H+CO+M
1 . 6 0 E14
0.000
14700.
HCO+H=CO+H2
4 . 0 0 E13
0.000
0.
HCO+O=CO2+H
1 . 0 0 E13
0.000
0.
HCO+O2=HO2+CO
3 . 3 0 E13
0.400
0.
CO+O+M=CO2+M
3 . 2 0 E13
0.000
4200.
CO+OH=CO2+H
1 . 5 1 E07
1.300
758.
CO+O2=CO2+O
1 . 6 0 E13
0.000
41000.
HO2+CO=CO2+OH
5 . 8 0 E13
0.000
22934.
H2+O2=2OH
1 . 7 0 E13
0.000
47780.
OH+H2=H2O+H
1 . 1 7 E09
1.300
3626.
H+O2=OH+O
5 . 1 3 E16
0.816
16507.
O+H2=OH+H
1 . 8 0 E10
1.000
8826.
H+O2+M=HO2+M
3 . 6 1 E17
0.720
0.
H2O/ 1 8 . 6 / CO2/ 4 . 2 / H2 / 2 . 8 6 / CO/ 2 . 1 1 / N2 / 1 . 2 6 /
OH+HO2=H2O+O2
7 . 5 0 E12
0.000
0.
H+HO2=2OH
1 . 4 0 E14
0.000
1073.
O+HO2=O2+OH
1 . 4 0 E13
0.000
1073.
2OH=O+H2O
6 . 0 0 E08
1.300
0.
H+H+M=H2+M
1 . 0 0 E18
1.000
0.
H+H+H2=H2+H2
9 . 2 0 E16
0.600
0.
H+H+H2O=H2+H2O
6 . 0 0 E19
1.250
0.
H+H+CO2=H2+CO2
5 . 4 9 E20
2.000
0.
H+OH+M=H2O+M
1 . 6 0 E22
2.000
0.
H2O/ 5 . 0 /
H+O+M=OH+M
6 . 2 0 E16
0.600
0.
H2O/ 5 . 0 /
H+HO2=H2+O2
1 . 2 5 E13
0.000
0.
HO2+HO2=H2O2+O2
2 . 0 0 E12
0.000
0.
H2O2+M=OH+OH+M
1 . 3 0 E17
0.000
45500.
H2O2+H=HO2+H2
1 . 6 0 E12
0.000
3800.
H2O2+OH=H2O+HO2
1 . 0 0 E13
0.000
1800.
!
END
Listing B.1: The used reaction mechanism in CHEMKIN format
Bibliography
[1] R. J. Kee, M. E. Coltrin, and P. Glarborg. Chemically reacting flow : theory and
practice. Hoboken, NJ: Wiley Interscience, 2003.
[2] World energy outlook. Paris, 2010.
[3] S. R. Turns. An introduction to combustion: concepts and applications. McGrawHill series in mechanical engineering. Boston: McGraw-Hill, 2000.
[4] N. Peters. Turbulent combustion. Cambridge monographs on mechanics. Cambridge: Cambridge University Press, 2000.
[5] E. O. Lawrence. Mixing and chemical reactions in turbulent flow reactors. PhD
thesis. University of California, 1964.
[6] T. Echekki, ed. Turbulent combustion modeling : advances, new trends and perspectives. Fluid mechanics and its applications ; 95. Springer, 2011.
[7] M. Ferziger J. H. ; Peri. Numerische Strmungsmechanik. Berlin: Springer, 2008.
[8] J. H. Chen. Petascale direct numerical simulation of turbulent combustion. In:
Proceedings of the Combustion Institute 33 (2011), pp. 99123.
[9] C. Cabrit and F. Nicoud. Direct simulations for wall modeling of multicomponent
reacting compressible turbulent flows. In: Physics of Fluids 21.055108 (2009).
[10] J. Warnatz, U. Maas, and R. W. Dibble. Combustion : physical and chemical
fundamentals, modeling and simulation, experiments, pollutant formation. Berlin:
Springer, 2006.
[11] M. Dsilets, P. Proulx, and G. Soucy. Modeling of multicomponent diffusion in
high temperature flows. In: International Journal of Heat and Mass Transfer
40.18 (1997), pp. 4273 4278.
[12] J. Keller. Influence of Molecular Diffusion Flux Modeling on the Structure of
Non-Premixed and Partially Premixed Flames. MA thesis. Technische Universitt
Bergakademie Freiberg, 2011.
[13] A. Spille-Kohoff, E. Preu, and K. Bttcher. Numerical solution of multi-component
species transport in gases at any total number of components. In: International
Journal of Heat and Mass Transfer 55.1920 (2012), pp. 5373 5377.
[14] A. W. Cook. Enthalpy diffusion in multicomponent flows. In: Physics of Fluids
21.055109 (2009).
[15] H. Jasak, A. Jemcov, and Z. Tukovic. OpenFOAM: A C++ Library for Complex
Physics Simulations. In: International Workshop on Coupled Methods in Numerical Dynamics. 2007.
Bibliography
60
[16] G. Dixon-Lewis. Flame Structure and Flame Reaction Kinetics. II. Transport
Phenomena in Multicomponent Systems. In: Proceedings of the Royal Society of
London. Series A, Mathematical and Physical Sciences 307 (1986), pp. 111135.
[17] D. Goodwin. Cantera: Object-oriented software for reacting flows. September 2012.
url: \\http://code.google.com/p/cantera/.
[18] M. Rehm. Numerische Strmungssimulation der Hochdruckvergasung unter Bercksichtigung detaillierter Reaktionsmechanismen. PhD thesis. Technische Universitt Bergakademie Freiberg, 2010.
[19] R. B. Bird, W. E. Stewart, and E. N. Lightfoot. Transport phenomena. New York:
Wiley, 2002.
[20] F.M. Rupley, R.J. Kee, and J.A. Miller. PREMIX: A Fortran Program for Modeling Steady Laminar One-dimensional Premixed Flames. In: Sandia National
Laboratories Tech. Report SAND85-8240 (1995).
[21] T. Poinsot and D. Veynante. Theoretical and numerical combustion. Philadelphia,
PA: Edwards, 2001.
[22] S. B. Pope. Turbulent flows. Cambridge University Press, 2000.
[23] L. P. H. de Goey and J. H. M. ten Thije Boonkkamp. A Mass-Based Definition
of Flame Stretch for Flames with Finite Thickness. In: Combustion Science and
Technology 122.1-6 (1997), pp. 399405.
[24] R. Brdika. Grundlagen der physikalischen Chemie. Berlin: Dt. Verl. d. Wiss.,
1977.
[25] J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird. Molecular theory of gases and
liquids. Structure of matter series. New York [u.a.]: Wiley [u.a.], 1954.
[26] S. Chapman and T. G. Cowling. The Mathematical theory of non-uniform gases :
an account of the kinetic theory of viscosity, thermal conduction and diffusion in
gases. Cambridge: Univ.Pr., 1970.
[27] M. Frenklach et al. GRI-Mech - An Optimized Detailed Chemical Reaction Mechanism for Methane Combustion. In: Gas Research Institute Topical Report GRI95/0058 (1995).
[28] C. R. Wilke. A Viscosity Equation for Gas Mixtures. In: The Journal of Chemical
Physics 18 (1950), pp. 517519.
[29] J. I. Steinfeld, J. S. Francisco, and W. L. Hase. Chemical kinetics and dynamics.
Englewood Cliffs, NJ: Prentice Hall, 1989.
[30] H.-P. Schmid. Ein Verbrennungsmodell zur Beschreibung der Wrmefreisetzung
von vorgemischten turbulenten Flammen. PhD thesis. Universitt Karlsruhe (TH),
1995.
[31] OpenFOAM - The Open Source CFD Toolbox, Programmers Guide. 2.1.0. OpenFOAM Foundation. 2011.
[32] OpenFOAM - The Open Source CFD Toolbox, User Guide. 2.1.1. OpenFOAM
Foundation. 2012.
Bibliography
61
[33] D. G. Goodwin. Cantera C++ Users Guide. California Institute of Technology.

2002.
[34] Regents of the University of California. BSD 3-Clause License. September 2012.
url: http://opensource.org/licenses/BSD-3-Clause.
[35] R.J. Kee et al. CHEMKIN-III: A FORTRAN chemical kinetics package for the
analysis of gas-phase chemical and plasma kinetics. Tech. rep. SAND96-8216. Sandia National Laboratories Report, 1996.
[36] B. Stroustrup. The C++ programming language. Boston [u.a.]: Addison-Wesley,
2009.
[37] Scott Meyers. Effective C++ : 55 specific ways to improve your programs and
designs. Addison-Wesley professional computing series. Upper Saddle River, NJ:
Addison-Wesley, 2007.
[38] R. J. Kee et al. Tech. rep. SAND85-8240. Sandia National Laboratories Report,
1985.
[39] K.J. Bosschaart and L.P.H. de Goey. The laminar burning velocity of flames
propagating in mixtures of hydrocarbons and air measured with the heat flux
method. In: Combustion and Flame 136 (2004), pp. 261269.
[40] Cray XE6. Cray Inc. 2010.
[41] M. Klein, A. Sadiki, and J. Janicka. A digital filter based generation of inflow data
for spatially developing direct numerical or large eddy simulations. In: Journal of
Computational Physics 186 (2003), pp. 652665.
[42] B. Krischok. Cray XE6 (HERMIT). 2012. url: http://www.hlrs.de/systems/
platforms/cray-xe6-hermit/.
[43] D. Parkinson. Parallel efficiency can be greater than unity. In: Parallel Computing 3 (1986), pp. 261262.

Bachelorthesis

Hochgeladen von

DokumentinformationenKlicken, um Dokumentinformationen aufzuklappen

Copyright:

Verfügbare Formate

Bachelorthesis

Hochgeladen von

Dokumentinformationen

Originalbeschreibung:

Originaltitel

Copyright

Verfügbare Formate

Dieses Dokument teilen

Dokument teilen oder einbetten

Freigabeoptionen

Stufen Sie dieses Dokument als nützlich ein?

Sind diese Inhalte unangemessen?

Copyright:

Verfügbare Formate

Bachelorthesis

Hochgeladen von

Copyright:

Verfügbare Formate

Karlsruhe Institute of Technology

Department of Chemical and Process Engineering

Implementation and Validation of a Solver

Evaluator: Prof. Dr.-Ing. Henning Bockhorn

Karlsruher Institut fr Technologie

Aufgabenstellung zur Bachelorarbeit

Implementierung und Validierung eines Lsers fr Direkte

Karlsruhe, October 2012 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2. Basic Aspects of Turbulent Combustion

3. Mathematical Description of Chemically Reacting Flows

4. Numerical Solution of Partial Differential Equations

4.1.2. Interpolation and Differentiation Procedures . . . . . . . . . . . .

5. Implemented Solver in OpenFOAM and the Cantera Interface

6. Validation of the Solver

7. Direct Numerical Simulation of a Turbulent Premixed Flame in 3D

8. Summary and Perspective

A. Code Fragments for OpenFOAM in C++3

3.1. Graphical expression of the Boltzmann equation. . . . . . . . . . . . . . .

4.1. Parameters used in the Finite Volume Method. . . .

5.1. Inheritance and dependency diagram for the coupling library. .

6.1. Physical domain of the laminar, premixed flame. . . . . . . . . . . . . . .

Numerical domain of the turbulent, premixed flame. . . . . . . . . . . . .

7.6. Different mass fractions and reaction rates. . . . . . . . . . . . . . . . . .

Physical conditions of the simulation with OpenFOAM. . . . . . . . . . .

B.1. The used reaction mechanism in CHEMKIN format . . . . . . . . . . . .

Specific standard-state heat of formation

Indices and Accents

Pressure implicit with split operator

1.3. Structure of the Thesis

1.3. STRUCTURE OF THE THESIS

2. Basic Aspects of Turbulent

2.1. Structure of Laminar Premixed Flames

2.1. STRUCTURE OF LAMINAR PREMIXED FLAMES

Xair /Xf uel

2.2. CHARACTERISTICS OF TURBULENT FLOW

2.2. Characteristics of Turbulent Flow

Here, v 0 is the turbulence intensity or velocity fluctuation, is the kinematic viscosity

2.2.1. Energy Cascade and Kolmogorov Scales

2.3. INTERACTION BETWEEN PREMIXED FLAMES AND TURBULENT FLOW

in which  is the rate of dissipation of turbulent kinetic energy defined to

2.3. Interaction between Premixed Flames and

2.3. INTERACTION BETWEEN PREMIXED FLAMES AND TURBULENT FLOW

2.3.1. Turbulent Flame Speed

2.3. INTERACTION BETWEEN PREMIXED FLAMES AND TURBULENT FLOW

In addition, the turbulent burning velocity can be calculated with

which is similar to equation (2.4) for the laminar burning velocity.

2.3.2. Flame Stretch Rate

With A = const., sL = const. and 0 = const. it must apply, that

2.3. INTERACTION BETWEEN PREMIXED FLAMES AND TURBULENT FLOW

2.3.3. Combustion Regimes

Fig. 2.5.: Schematic classification of turbulent combustion regimes.

2.3. INTERACTION BETWEEN PREMIXED FLAMES AND TURBULENT FLOW

3.1. Statistical Thermodynamics and the Rigorous

3.1. STATISTICAL THERMODYNAMICS AND THE RIGOROUS KINETIC THEORY

3.1.1. The Boltzmann Equation

Fig. 3.1.: Graphical expression of the Boltzmann equation.

The Boltzmann equation can be viewed as a transport equation in the six-dimensional

in which is the rate of dissipation of turbulent kinetic energy defined to