Manual
User Manual:
Open the PDF directly: View PDF 
.
Page Count: 17
User Guide for Detailed Modeling of Laminar Flames
and High-Resolution LES in OpenFOAM R
Adhiraj Dasgupta
December 2, 2018
Contents
1 Numerical Model 1
1.1 The fully-compressible formulation ......................... 1
1.2 The low-Mach-number formulation ......................... 2
1.3 The transport model .................................. 2
1.4 Using reduced chemical mechanisms ......................... 6
2 Installation and Running 7
2.1 Installation ...................................... 7
2.2 Setting up and running ................................ 7
2.2.1 The transportProperties file .......................... 8
2.2.2 Modifications to the file fvSchemes ...................... 10
2.2.3 The keyword pCorr .............................. 10
2.2.4 Fitting the transport properties and running ................. 10
2.2.5 Using reduced chemical mechanisms .................... 11
2.2.6 Tips and tricks ................................ 11
i
Abstract
In this work a set of libraries and solvers have been developed within OpenFOAM R
’s code frame-
work that can be used to model laminar flames as well as perform high-resolution LES of reacting
flows. Transport properties are calculated using a detailed, mixture-averaged model, and both a
fully-compressible, and an ‘incompressible’, low-Mach-number solver have been developed. This
theory manual provides a brief summary of the governing equations and models used in the code;
it also serves as an installation guide, and provides a step-by-step guide for setting up a simulation
along with all the settings with their recommended values.
Additionally a new feature has been added to complement OpenFOAM R
’s chemistry model
that allows for the use of reduced chemical mechanisms that are available as Fortran subroutines.
Chapter 1
Numerical Model
The governing equations for reacting flows are essentially the statements of conservation of mass,
momentum, species, and energy. A detailed summary has been provided elsewhere [1,2], and only
a brief overview will be provided here.
1.1 The fully-compressible formulation
The fully-compressible version of the conservation equations may be summarized as
∂ρ
∂t +∂(ρui)
∂xi
= 0 ,(1.1)
∂ρui
∂t +∂(ρujui)
∂xj
=−∂p
∂xi
+∂τij
∂xj
+ρfi,(1.2)
∂ρYk
∂t +∂(ρujYk)
∂xj
=−∂Jk,j
∂xj
+ ˙ωk,(1.3)
∂ρhs
∂t +∂(ρujhs)
∂xj
+∂ρK
∂t +∂(ρujK)
∂xj
=−∂qj
∂xj
+∂p
∂t +ρuifi+qrad
−
N
X
l=1
˙ωl∆h◦
f,l +∂(τij ui)
∂xj
,(1.4)
p=ρRuT
N
X
l=1
Yl
Wl
.(1.5)
Here ρis the gas density, uiis the gas velocity, pis pressure, τij is the viscous stress tensor, Yk
is the mass fraction of species k,Jk,j is the molecular diffusion flux of the kth species, ˙ωkis the
production rate of species k,hsis the mixture sensible enthalpy, Kis the kinetic energy per unit
mass of the flow, qiis the heat flux, fiis the body force per unit volume (e.g., due to gravity),
qrad is the radiative heat source, ∆h◦
f,k is the standard state heat of formation of species k,Ruis
the universal gas constant, and Wkis the molecular weight of species k. Also, the term ∂(τij ui)
∂xj
contains the the viscous dissipation function Φ = τij ∂ui
∂xj.
1
1.2 The low-Mach-number formulation
In the low-Mach-number or the ‘incompressible’ formulation the pressure driving the flow is
decoupled from the thermodynamic state of the system; the pressure is split into two parts as
p=p0+pd[3,4]. Here p0is the thermodynamic pressure (assumed constant), and the hydrody-
namic pressure pddrives the flow field. The governing equations are also modified; the flow kinetic
energy Kis considered negligible, and the viscous dissipation term is also dropped. The resulting
equations are Eqs. (1.1), (1.3), and
∂ρui
∂t +∂(ρujui)
∂xj
=−∂pd
∂xi
+∂τij
∂xj
+ρfi,(1.6)
∂ρhs
∂t +∂(ρujhs)
∂xj
=−∂qj
∂xj
+qrad −
N
X
l=1
˙ωl∆h◦
f,l ,(1.7)
p0=ρRuT
N
X
l=1
Yl
Wl
.(1.8)
When running high-resolution, implicit LES all the variables in these equations may be inter-
preted as filtered quantities.
1.3 The transport model
In this work mass diffusion is modeled using the Hirschfelder and Curtiss approximation [5]. Mass
conservation is ensured by using a correction velocity approach [6,7]. Thus the diffusion flux for
species kis given as
Jk,i =ρ−Dk
∂Xk
∂xi
+YkVc
i,
where Xkis the mole fraction of species k,Dkis the mixture-averaged diffusivity for species k,
and Vicis a spatially-varying correction velocity, which is the same for all species. The heat flux
is given by
qi=−λ∂T
∂xi
+
N
X
l=1
hlJl,i,
where λis the thermal conductivity of the gas mixture, and hlis the enthalpy of species l. Lastly,
the viscous stress tensor τij is computed as
τij = 2µSij −2
3µSkk δij ,
where µis the dynamic viscosity of the gas mixture, and
Sij =1
2∂ui
∂xj
+∂uj
∂xi
2
is the strain rate tensor.
Evaluation of the transport properties follows the procedure outlined by Kee et al. [6]. The
viscosity of the gas-phase species kis given by [5,6]
µk=5
16
√πmkkBT
πσ2
kΩ(2,2)∗,(1.9)
where σkis the Lennard-Jones collision diameter, mkis the molecular mass for species k,kB
is Boltzmann’s constant, and Ω(2,2)∗is the reduced collision integral, which is a function of the
reduced temperature T∗
kgiven by
T∗
k=kBT
εk
,
where εkis the Lennard-Jones potential well depth for species k. The collision integrals are eval-
uated from the reduced temperature using curve-fits generated by Neufeld et al. [8].
The binary diffusion coefficient for gas-phase species jand kis given by [5,6]
Djk =3
16 p2πk3
BT3/mjk
pπσ2
jkΩ(1,1)∗,(1.10)
where mjk is the reduced molecular mass for species jand k, defined as
mjk =mjmk
mj+mk
,
σjk is the reduced colllision diameter, and the reduced collision integral Ω(1,1)∗is a function of the
reduced temperature T∗
jk
1. These parameters are calculated based on whether the two molecules
involved are both polar, both non-polar or one polar and the other non polar.
For the case of two non-polar molecules or two polar molecules interacting, the reduced quan-
titites are defined as
εjk
kB
=sεj
kBεk
kB,
σjk =1
2(σj+σk).
For a non-polar molecule interacting with a polar one, the reduced polarizability α∗
nand the reduced
dipole moment µ∗
pare defined as
α∗
n=αn
σ3
n
,
µ∗
p=µp
pεpσ3
p
,
1In the low-Mach-number formulation p0is used in place of p.
3
where αdenotes the polarizability and µdenotes the dipole moment; the subscript prefers to the
polar species and the subscript nto the non-polar species. The reduced temperature T∗
jk is defined
as
T∗
jk =kBT
εjk
.
The reduced quantities are calculated as
εnp
kB
=ζ2sεn
kBεp
kB,
σnp =1
2(σn+σp)ζ−1
6,
where
ζ= 1 + 1
4α∗
nµ∗
prεp
εn
.
The thermal conductivity is assumed to have contributions from translational, rotational, and
vibrational components. The thermal conductivity of species kis given by [6,9]
λk=µk
Wk
(ftrans,k Cvk,trans +frot,kCvk,rot +fvib,k Cvk,vib),(1.11)
where Cvk,trans,Cvk,rot, and Cvk,vib are respectively the translational, rotational, and vibrational
contributions to the constant-volume specific heat Cvk of species k,
ftrans,k =5
21−2
π
Cvk,rot
Cvk,trans
Ak
Bk,
frot,k =ρkDkk
µk1 + 2
π
Ak
Bk,
fvib,k =ρkDkk
µk
.
Here the self-diffusion coefficient Dkk is obtained by substituting j=kin Eq. 1.10, and the density
of gas-phase species kis obtained as
ρk=pWk
RuT.
Further,
Ak=5
2−ρkDkk
µk
,
Bk=Zrot,k +2
π5
3
Cvk,rot
Ru
+ρkDkk
µk.
4
Here Zrot,k is the rotational relaxation number for species kand is given by
Zrot,k (T) = Zrot,k(298)Fk(298)
Fk(T),
where
Fk(T) = 1 + π3
2
2εk/kB
T1
2
+π2
4+ 2εk/kB
T+π3
2εk/kB
T.
The translational and rotational components of the gas specific heat are dependent on the kind of
molecule under consideration. For a linear molecule,
Cvk,trans
Ru
=3
2,
Cvk,rot
Ru
= 1,
Cvk,vib =Cvk −5
2Ru.
For a non-linear molecule
Cvk,trans
Ru
=3
2,
Cvk,rot
Ru
=3
2,
Cvk,vib =Cvk −3Ru.
For a single atom (for example H, O, etc.) there are no internal contributions to the specific heat.
Hence
Cvk,trans
Ru
=3
2,
Cvk,rot
Ru
= 0,
Cvk,vib = 0.
Transport properties for the mixture are computed using a mixture-averaged approach. The
mixture-averaged diffusivities are obtained as [6,10]
Dk=1−Yk
PN
l=1
l6=k
Xl
Dkl
,(1.12)
where Dkl is the binary diffusivity for the species pair kand l. The mixture viscosity is obtained
as [6,11]
µ=
N
X
l=1
Xlµl
PN
j=1 XjΦlj
,(1.13)
5
where µlis the dynamic viscosity of species l, and
Φmn =1
√81 + Wm
Wn−1
2"1 + µm
µn1
2Wn
Wm1
4#2
.(1.14)
The mixture-averaged thermal conductivity is obtained as [6,12]
λ=1
2 N
X
l=1
Xlλl+1
PN
l=1 Xl/λl!,(1.15)
where λlis the thermal conductivity of species l.
To speed up calculations, the temperature-dependent parts of the transport properties for the
individual species are fitted using the method of least squares implemented in the GNU Scientific
Library [13]. The logarithms of the properties are fitted to polynomials in the logarithm of temper-
ature, as was done by Kee et al. [6]. It is to be noted that the viscosity and the thermal conductivity
of species kare a function of temperature only, and thus these curve-fits are straightforward and
unambiguous. However, the binary diffusivities are functions of temperature and pressure; the fits
are therefore performed at 1bar pressure. The values obtained from the curve fits are divided by
the pressure in bar in order to obtain the actual diffusivities. Details of the fitting process can be
found elsewhere [2].
1.4 Using reduced chemical mechanisms
In some cases the chemical mechanism to be used is available only as a Fortran subroutine that
evaluates the reaction rates for the species. The commonly used quasi-steady-state (QSS) assump-
tion provides such an example. In this case since fewer species are transported, some savings in
computational effort can be expected.
6
Chapter 2
Installation and Running
The present code framework is based on OpenFOAM R
version 2.3.x. It is assumed that the user
has a working installation of this version of OpenFOAM R
on their computer.
2.1 Installation
The transport model library depends on the mole fraction library for some calculations, and so the
latter must be built first. The easiest way to do that is to copy the moleFractions directory to the
user’s directory and issuing the command wmake libso in the terminal.
Once the mole fraction library is successfully built, the transport library may be built. The user
may need to modify the provided options file in order to specify the location of the user’s source
code and the user’s libraries. Then the command wmake libso should work.
If all goes well, as they should, at this stage the libraries have been built successfully; it remains
now to build the applications that use the libraries. The transport fitting program must be built and
executed before any solvers can be run. This utility is called fitTransport, and it has the GNU
Scientific Library (https://www.gnu.org/software/gsl) as a dependency. If GSL is installed in the
standard locations, as would be the case if it is built using the default options, the provided options
file should work. Issuing the command wmake at the command prompt should install the fitting
program.
The compressible solver is called laminarReactingFoam and the low-Mach-number solver
is called laminarReactingLMFoam which is influenced by the reactingLMFoam code [4]. For
each solver a sample options file is provided, and the user may need to change it to include the lo-
cations of the user’s code and libraries. Then the solvers may be built by issuing the command
wmake in the corresponding directory.
2.2 Setting up and running
Setting up a case for use with the solver is similar to the way the same case would be set up for
reactingFoam, with a few key differences:
7
1. A file called transportProperties must be placed in the constant directory of the case.
2. For the low-Mach-number solver the variable pReff (corresponding to p0) must be provided
in the file chemistryProperties in the constant directory.
3. Additional numerical schemes must be specified in the file fvSchemes in the system directory
of the case.
4. For running the low-Mach-number solver, solver options must be specified for the dynamic
pressure pdin the file fvSolutions in the system directory of the case.
5. Optionally, the extra term in the pressure equation can be turned on or off through a flag set
in the file thermoPhysicalProperties in the constant directory of the case.
2.2.1 The transportProperties file
The file transportProperties has the following entries:
transportModel
This option tells the code which transport model to use. The available options are
1. mixtureAverage,
2. LewisNumber.
For each of these options, a corresponding subdictionary must be specified in the file.
Subdictionary mixtureAverageProperties
This subdictionary contains the options specific to running with the mixture-averaged transport
model. The entries are:
thermophoresis This option indicates whether thermophoresis will be included for light
species. By default this is set to off.
CutOff This is the molecular weight below which thermophoresis will be computed.
The default value is 5.
8
gradX In the original formulation for Fickian diffusion velocities, the flux is taken to
be proportional to the gradient of the mole fraction, and the diffusion velocity
of species kis given as
Vk,j =−Dk
Xk
∂Xk
∂xj
=−Dk
Yk
∂Yk
∂xj−Dk
¯
M
∂¯
M
∂xj
,
where ¯
Mis the mixture molecular weight.
If this option is set to on, the above equation is used to calculate the diffusion
velocity. However, if it is set to be off, the flux is taken to be proportional to
the gradient of the mass fraction,
Vk,j =−Dkm
Yk
∂Yk
∂xj
.
By default this is set to be on.
Subdictionary LewisNumberProperties
This subdictionary contains the options specific to running with the constant Lewis number model.
The options are to be specified in the form of a subdictionary of LewisNumberProperties, named
Le.
Le This contains the names of the species with the corresponding Lewis numbers to be
used for each species. By default the Lewis numbers are taken to be unity, so a
simulation with an empty Le dictionary will enforce unity Lewis number for all species.
CHEMKINFile
This is the path to the kinetics file in Chemkin format [14].
CHEMKINThermoFile
This is the path to the thermo data in Chemkin format.
CHEMKINTransportFile
This is the path to the transport database in Chemkin format.
9
Tmin and Tmax
The transport properties are fitted between these two temperatures. By default they are taken to be
300 K and 3000 K.
viscousDissipation
For the compressible formulation, this option indicates whether or not the viscous dissipation term
τij ∂ui
∂xjshould be included in the energy equation. By default this is set to off.
2.2.2 Modifications to the file fvSchemes
The numerical schemes and solver options are very similar to the ones used by reactingFoam.
In addition, some new entries are added to the divSchemes subdictionary. The new keywords are
1. div(phi,Yi),
2. div((mu*dev2(T(grad(U))))),
3. div(JHs),
4. div(Hconduction) or div(Econduction), and
5. div(phi,he).
Additionally, if the low-Mach-number solver is used, the variable pd must be added to the fluxRe-
quired subdictionary.
2.2.3 The keyword pCorr
The implementation of the pressure equation in OpenFOAM R
contains an extra flux term, likely
for added stability. This term can cause problems in some cases. In this code the user has the
option of disabling or enabling this extra term through the keyword pCorr in the file thermoPhysi-
calProperties.
By default the term is dropped.
2.2.4 Fitting the transport properties and running
The utility fitTransport generates fits to the transport properties. The coefficients for the
fits are written to the following files in the constant directory of the case: conductivityProperties,
diffusivityProperties,thermophoreticProperties, and viscosityProperties. For each species, these
files contain the fit coeffients, and the fitting error.
In addition, the molecular parameters read from the transport database are summarized in the
file transportData in the constant directory. Also, the transport properties are written out as a
function of temperature in a new directory called transport.
Once the fits are generated, the code can be run simply by typing laminarReactingFoam,
or laminarReactingLMFoam.
10
2.2.5 Using reduced chemical mechanisms
As an added option, the present work implemented a way to incorporate reduced chemical mecha-
nisms within the OpenFOAM R
framework. The relevant code is placed under the provided chem-
istryModel directory that should be copied to the user’s own directory. For now this implemen-
tation only works for mechanisms that are available as Fortran routines called ckwyp, but may be
extended easily to include other forms. Compiling this piece of code produces the shared library
libCustomChemistryModel.so.
It should be noted that the library does not by itself contain the mechanism information; the
Fortran mechanism must be compiled separately into a shared library. Since the description of the
chemical mechanism is incomplete, directly linking the library libCustomChemistryModel.so
to the solvers will not work. Thus the two libraries—libCustomChemistryMode.so, and the
Fortran shared library must be added to the controlDict of the case.
To use the new model, the following steps should be taken:
1. The Fortran mechanism must be compiled to produce a shared library.
2. Both libraries should be added to the file controlDict in the system directory.
3. The file chemistryProperties should be modified:
•The keyword chemistrySolver should be set to QSS.
•A subdictionary called QSSCoeffs should be added. The content of this should be
similar to the odeCoeffs, since the stiff ODE system is integrated the same way; only
the reaction rates are computed in a new way.
4. For a mechanism using the QSS assumption there are not reactions specified in the Chemkin
mechanism file; however a dummy reaction must be specified to keep OpenFOAM R
’s Chemkin
parser happy. For instance, if modeling ethylene combustion, the reaction block can look like
this:
REACTIONS
C2H4+3O2 =>C2H4+3O2 0.0 0.0 0.0
END
2.2.6 Tips and tricks
Most of the other settings in use are analogous to the corresponding settings used by reactingFoam.
However, a few settings were found to work well, and as such may be used as rough guidelines.
They are mentioned here:
•Since a lot of terms in the equations are directly calculated from gradients, very sharp gradi-
ents can be problematic. The cellMDLimited scheme has been found to be stable.
11
•Whenever possible, running with the low-mach-number solver is to be preferred, for numeri-
cal stability even when using relatively coarse time steps. The savings in computational time
are likely to be at least an order of magnitude.
•The term pCorr may cause problems, and it may be better to keep it turned off unless there
is evidence to indicate otherwise.
12
Bibliography
[1] A. Dasgupta, E. Gonzalez-Juez, and D. C. Haworth, “Flame simulations with an open-source
code,” Computer Physics Communications, 2018.
[2] A. Dasgupta, Numerical Simulation of AxiSymmetric Laminar Diffusion Flames with Soot,
Ph.D. thesis, Pennsylvania State University, 2015.
[3] A. G. Tomboulides, J. C. Y. Lee, and S. A. Orszag, “Numerical Simulation of Low Mach
Number Reactive Flows,” Journal of Scientific Computing, vol. 12, no. 2, pp. 139–167, Jun
1997.
[4] K.-J. Nogenmyr, C. Chan, and C. Duwig, “Finite rate chemistry effects and combustor liner
heat transfer studies in a frame-work of LES of turbulent flames-investigation of pollutant
formation using OpenFOAM,” 5th OpenFOAM Workshop, Chalmers, Gothenburg, Sweden,
2010.
[5] J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird, Molecular Theory of Gases and Liquids,
John Wiley & Sons, New York, 1994.
[6] R. J. Kee, G. Dixon-Lewis, J. Warnatz, M. E. Coltrin, and J. A. Miller, “A Fortran Com-
puter Code Package for the Evaluation of Gas Phase Multicomponent Transport Properties,”
Sand86-8246 unlimited release, Sandia National Laboratory, 1986.
[7] T.P. Coffee and J.M. Heimerl, “Transport algorithms for premixed, laminar steady-state
flames,” Combustion and Flame, vol. 43, pp. 273 – 289, 1981.
[8] P. D. Neufeld, A. R. Janzen, and A. R. Aziz, “Empirical Equations to Calculate 16 of the
Transport Collision Integrals Ω(l,s)∗for the Lennard-Jones (12–6) Potential,” The Journal of
Chemical Physics, vol. 57, pp. 1100–1102, 1972.
[9] N. Peters and J. Warnatz, Eds., Numerical Methods in Laminar Flame Propagation, Vieweg
& Teubner Verlag, 1982.
[10] R. B. Bird, W. E. Stewart, and E. N. Lightfoot, “Transport phenomena,” 1960.
[11] C. R. Wilke, “A Viscosity Equation for Gas Mixtures,” The Journal of Chemical Physics,
vol. 18, no. 4, pp. 517–519, 1950.
13
[12] S. Mathur, P. K. Tondon, and S. C. Saxena, “Thermal conductivity of binary, ternary and
quaternary mixtures of rare gases,” Molecular Physics, vol. 12, pp. 569–579, 1967.
[13] Brian Gough, GNU Scientific Library Reference Manual - Third Edition, Network Theory
Ltd., 3rd edition, 2009.
[14] R. J. Kee, F. M. Rupley, and J. A. Miller, “CHEMKIN-II: A Fortran Chemical Kinetics
Package for the Analysis of Gas-Phase Chemical Kinetics,” Tech. Rep. SAND89-8009B,
Sandia National Laboratories, 1989.
14
