1 Introduction

Discrete kinetic theory, and notably the Lattice Boltzmann method has met
with significant success in the recent years for the numerical simulation of a
arXiv:comp-gas/9609003v1 20 Sep 1996

wide host of fluid phenomena, ranging from multiphase flows in grossly ir-
regular geometries, up to highly turbulent homogeneous incompressible flows
[1, 2, 3]. Main drivers behind the method are ideal amenability to parallel
computing, easy handling of grossly irregular geometries and physical sound-
ness.
The Lattice Boltzmann Equation (LBE) is more than a plain Navier-
Stokes solver; actually, it is a minimal hyperbolic superset of the Navier-Stokes
equations. In fact, in the process taking from the kinetic to the hydrodynamic
level, along with standard primitive fluid fields such as density, speed and
temperature, there is enough information left to describe the dynamics of
a genuinely hydrodynamic stress tensor as well as that of a passive scalar.
This property permits to extend the applicability of LBE techniques beyond
’plain’ hydrodynamics, with only minor modification to the basic discrete
equation. To date, this feature has been exploited in various ways, to describe
thermal convection, magnetohydrodynamics ([4] and references therein) and
also subgrid modeling [5].
All of the above refers to the 24-speed FCHC (Face Centered HyperCube)
lattice [6], the ancestor of modern Lattice Boltzmann schemes, introduced
by D’Humieres, Frisch and Lallemand back in 1987. Recently, FCHC has
been superseded by swifter Lattice BGK schemes [7, 8] which can boast a
simpler (diagonal) collision operator as well as a significantly smaller number
of discrete speeds (14 instead of 24).
In spite of these undeniable advantages, we shall argue that the LBE-
FCHC scheme still holds a certain interest for it provides the unique oppor-
tunity to integrate passive and active scalars dynamics within the very same
scheme tracking the flow variables, with no need to extend the corresponding
set of discrete speeds.
The scope of this paper is twodfold; first, to revisit the symmetry proper-
ties of LBE-FCHC schemes, second to present a novel potential application
of these symmetries to the field of reactive flow dynamics.

1
2 Reactive flow dynamics
The set of hydro-chemical equations governing the dynamics of a S-species
reactive flow read as follows (k, l run over the spatial dimension):

∂t ρs + ∂l (ρs ul ) = ∂l Ds ρ∂l (ρs /ρ) + Ss , s = 1, S (1)

∂t ρuk + ∂l (ρuk ul ) = ∂l νρ∂l uk − ∂k p (2)

∂t ρT + ∂k (ρuk T ) = ∂l χρ∂l T + Q̇ (3)

where ρs is the density of the s − th species, ρ is the fluid density (ρ =

P
s s ), T the fluid temperature, p the fluid pressure, uk the flow field, D, ν
ρ
and χ the mass, momentum and energy kinetic diffusivity of the fluid. Here
Ss is the mass density input to species s from chemical reactions, and Q̇ is the
P
heat input rate as obtained by summing all over the reactions: Q̇ = r Qr ω̇r
, where ω̇r is the progress rate (moles/s) of reaction r.
Combustion is governed by the competition between hydrodynamics (dif-
fusion and convection) and chemistry; the relative strength of this competion
is measured by a dimensionless parameter known as the Damkohler number,
defined as Da = ττhc where τh and τc are typical hydrodynamic and chemical
timescales respectively,
Low Damkohler’s characterize quasi-inert flows while high Damkholer’s
are associated with fast (τc << τh ) chemical reactions.
In this paper we shall be concerned with the ’fast’ combustion limit, i.e.
reactive flows in which chemistry proceeds much faster than hydrodynamics
so that as soon as mixing is completed the fuel instantly burns (”mixed-is-
burned”).
In this limit, combustion tends to localize on thin regions (”reaction
sheets”) whose thickness becomes vanishingly small as the inverse Damkohler
number goes to zero [9]. This brings about a drastic simplification in the
mathematical model because the focus can be shifted onto the topology of
the reaction sheet. It can be shown that for the case of non-premixed flames,

2
i.e. flames where fuel and oxidizer are fed into the system in separate streams,
and under the equidiffusional assumption (D = ν = χ), the topology of the
reaction sheet is properly described by a passive scalar, say Z, evolving ac-
cording to the usual advection-diffusion equation:

ρDt Z = ∂k ρD∂k Z (4)

where Dt = ∂t + u∂x + v∂y + w∂z is the fluid substantial derivative.
The underlying idea is that while species transform one into another un-
der reactive events, there are quantities (typically elemental mass fractions)
which are unaffected by chemistry and serve therefore as slow modes en-
slaving the system dynamics. As a result, in the limit where chemistry
proceeds much faster than hydrodynamics, all species concentrations are
’frozen-in’ byproducts of the conserved scalars. With reference to a diffu-
sive non-premixed flame, it is readily checked that a convenient choice for Z
is provided by the fuel mixture fraction, defined as the mass ratio of fuel in
the unburned mixture versus fuel mass ratio in the stream,

Z = Yf u /Yf s (5)
where Y stands for mass fraction and the subscripts u, s refer to the
unburned and stream fuel respectively.
For a single-step irreversible reaction of the form
F uel + Oxidizer → P roducts
combustion takes place around the surface where reactants and products
come in stoichiometric proportion. The stoichiometric surface is defined by
the relation:
Z(x, y, z) = Zst = 1/(1 + AYf,s /Yo,s )
where A is the ratio of the fuel to oxidizer stoichiometric coefficients
times the corresponding molecular weights [10]. In the limit of infinitely fast
chemistry, the stoichiometric surface represents a front of discontinuity for
the spatial derivatives of the hydrodynamic variables. Within this limit, and
with the additional assumption of steady state and adiabatic conditions, all
quantities become a simple piecewise linear map of Z. With reference to lean
(0 < Z < Zst , superscript ”L”) and rich (Zst < Z < 1), superscript ”R”)
portions of the mixture, these mappings read as follows:

YfLb (Z) = 0 (6)

3
 YfRb (Z) = Yf s (Z − Zst )/(1 − Zst ) (7)

YobL (Z) = Yos (1 − Z/Zst) (8)

YobR (Z) = 0 (9)

TbL (Z) = Tin + (Z/Zst)(Tst − Tin ) (10)

1−Z
TbR (Z) = Tin + ( )(Tst − Tin ) (11)
1 − Zst
where subscript ”b” means ”burned”.
Here Tin is the inlet stream temperature (we assume stream fuel and
oxidizer at the same temperature) and the actual value of Tst depends on the
specific chemical reaction. Note that in the unburned portion of the mixture
the dependence on Z disappears (Tu = Tin ).
Within this approximation the heat source reduces to Q̇ = Qr ρ∆Yf /Wf δ(Z−
Zst ) where Qr is the reaction heat (Joule/mole) and ∆Yf represents the
change of the fuel mass fraction between the unburned (Yf = ZYf s )) and
the burned (Yf = Yf b ) state. The presence of the Dirac’s delta reflects the
instantaneous nature of the chemical reactions.

3 Extended LBE for reactive flows

The LBE is a minimal discrete Boltzmann equation reproducing Navier-
Stokes hydrodynamics in the limit of small Knudsen numbers, i.e. particle
mean free path much smaller than typical macroscopic variation scales [2].
We shall refer to the 24-speed, single-energy, four-dimensional FCHC (Face
Centered HyperCube) lattice defined by the condition [6]
P 2
k=1,4 cik = 2 cik = (0, ±1)
The FCHC-LBE takes the form of a set of first-order finite difference
equations:

X
24
fi (xk + cik , t + 1) − fi (xk , t) = Aij (fj − fje ) (12)
j=1

4
 where fi , i = 1, 24 is a set of 24 populations moving along a corresponding
set of discrete speeds cik and k = 1, 4 runs over the spatial dimensions.
The eq.(12) is naturally interpreted as a multi-relaxation scheme in which
non-equilibrium gradients are brought back to the local equilibrium fje by
scattering events mediated by the matrix Aij . The local equilibrium popula-
tions are chosen in the form of a discretized Maxwellian expanded to second
order in the Mach number in order to retain convective effects.
ρ 1
fie = (1 + 2cik uk + 2(cik cil − δkl )uk ul ) (13)
24 2
Under the constraint of point-wise mass and momentum conservations
P P
( i Aij = i cik Aij = 0) the equation (12) describes the motion of a four-
dimensional fluid. The reason for working in four-dimensions is that no
single-energy discrete lattice exists in three-dimensional space which ensures
P
isotropy of the fourth order tensor Tklmn = i cik cil cim cin .
Three-dimensional motion is then recovered by quenching the propagation
along the discrete speeds projecting upon the fourth dimension x4 . Quench-
ing the fourth dimension releases an excess of populations versus discrete
speeds (24 against 18) introducing ’de facto’ a degeneracy of the ’momentum
eigenstates’ (cik ) propagating along the fourth dimension, since each of these
links carries two populations with opposite speeds along x4 .
Generally, to the purpose of describing sheer hydrodynamics, there’s no
point of keeping two distinct populations (”doublet”) along the same discrete
speed. A 18 speed-18 populations (18s18p for brevity) scheme is perfectly ad-
equate for three-dimensional Navier-Stokes equations. If, on the other hand,
doublets are retained, then the current J4 is readily shown to be passively
trasported by the three-dimensional fluid; this means that the 24s18p model
delivers the 3D Navier-Stokes equations plus a passive scalar. The same
reasoning can be down-iterated to lower dimensions leading to 9s24p, 9s18p
(two-dimensional fluid plus one/two passive scalars respectively) and 3s24p,
3s18p 3s9p(one dimensional fluid and one/two/three passive scalars) in 1D.
A simple sketch illustrates the tree structure associated with dimensional
quenching of the FCHC scheme:

5
 24 4D: 4D fluid, 0 passive scalar (24s24p)
* -
18 6 3D: 3D fluid, 0 passive scalar (18s18p)
* - -
9 9 6 2D: 2D fluid, 0 passive scalar (9s9p)
* - - - 2 (9s24p)
3 6 9 6 1D: 1D fluid, 0 passive scalar (3s3p)
1 (3s9p)
2 (3s18p)
3 (3s24p)
Since passive/active scalars are in heavy use in numerical combustion to
describe fast chemistry regimes, we argue that the aforementioned properties
point to the FCHC Lattice Boltzmann as to a potential candidate for the
numerical description of combustion phenomena.
In this paper, we shall focus on two-dimensional reactive flow dynamics
with one passive scalar describing the fuel mixture fraction Z and a second,
active, scalar accounting for temperature dynamics. This picks up the 9s24p
model. It should be noted that within the ”mixed-is-burned” assumption,
the temperature field can be obtained directly from the knowledge of the
mixture fraction Z via the piecewise linear mapping given by eq.(10-11). We
choose however to track the temperature field independently of Z in order
to gain a qualitative feeling for the validity of the adiabatic approximation.
With these preparations, the reactive flow variables are defined as follows:
P P P P
ρ = 24 fi , ρu = 24 fi ci1 , ρv = 24 fi ci2 , ρZ = 24
i=1 fi ci3 , ρT =
P24 i=1 i=1 i=1
i=1 fi ci4
The second scalar T (temperature) is ’activated’ by adding a source term,
say q to all populations projecting upon the fourth dimension. A simple
momentum budget fixes q as a function of the heat released Q̇. This is all
we need to describe fast combustion.

4 Numerical results
As a test-case application, we shall consider a coflow two-dimensional methane-
air laminar flame in a plane domain. The numerical set-up is as follows:
Inlet: Coflowing Fuel (methane)-Oxidizer (air) at T = Tin . The fuel
enters the domain within an inner region centered about the midline and

6
width H/8. The mixture fraction at the inlet is Zin . Outlet: Zero cross-
flow and zero axial derivatives of the streamwise velocity. Upper and lower
boundary are kept at a constant wall temperature Tw . Other simulation
parameters are: inlet mixture fraction profile: Zin = 2.0, Zst G(y), G being a
gaussian profile of width H/8, inlet temperature: Tin = 0.0 (conventionally
800 K), inlet flow: Poiseuille profile, maximum Uin =0.1, aspect ratio L/H =
4.
Our baseline computation is performed on a 256x64 grid and takes about
1µs/step/site on a SUN Sparc1 workstation, ninety percent of which due to
hydrodynamics. Since no turbulence model has been implemented, the flow
in H
viscosity (ν = 0.01) corresponds to a Reynolds number 2U3ν ∼ 420. The
simulation is run over 3-5 axial residency times τR = 1.5L/Uin = 3840.
Chemistry is represented by the following one-step global reaction:
CH4 + 2O2 = 2H2 0 + CO2 (14)
which is supposed to proceed instantly at the stoichiometric surface.
For practical purposes, the reaction is assumed to occurr within a shell
centered about Z = Zst = 0.057 with thickness Zth = 0.2Zst . The heat
source is assumed a (piecewise) constant within this shell.
In Fig.1 we show the equicontours of the mixture fraction at time t=4,000,
corresponding to about one recirculation time. From this picture the typical
shape of a laminar diffusion flame is apparent [11, 12]. Note that the sto-
ichiometric surface is confined to the initial region of the domain (x ≃ 50
lattice units). For a given set of physical-geometrical parameters, this loca-
tion, namely the flame length Lf , depends solely on Zin , i.e the degree of
mixture richness.
In Fig 2 we show the mixture fraction at the midline y = H/2 as a function
of the streamwise coordinate x at two different time snapshots t = 2000, 4000.
The advancement of the mixture front is well visible along with its diffusive
spreading. As expected, after an initial diffusive transient lasting about
τR /Re ≃ 100 time steps, the tip of the front proceeds on a quasi-convective
regime at speed Uf = 2Uin /3. The flame however stabilizes around the
stoichiometric surface long before the tip comes to the outlet of the domain
because Lf is much shorter than the longitudinal span L of the computational
box.
Within the equilibrium chemistry assumption, mid-plane reactions (y =
H/2) should peak around the stoichiometric line x = xst , where Z(xst , y =

7
H/2) = Zst . In our simulation, xst ≃ 52 as highlighted by the horizontal
line in Fig.2. This is the region where a substantial temperature buildup is
expected to occurr. Such a prediction is indeed reflected by Fig 3, where the
median temperature profile is reported as a function of the axial coordinate x.
In line with the behaviour of the mixture fraction, the temperature field shows
a peak centered about the stoichiometric region, with a top temperature of
about 2100 K, slightly below the stoichiometric temperature Tst = 2263K.
For the sake of comparison, the adiabatic profile deduced by the theoretical
”frozen-in” mapping T (Z) computed with Z at t = 4000 is reported right
above (dotted line on top). From this figure, we see that a relatively good
match is obtained in in the post-flame region, while in the pre-flame region
the adiabatic profile provides a gross overestimate of the temperature field.
This is not surprising, since the near-inlet region is characterized by strong
gradients which make the adiabatic assumption rather questionable [13].
Finite-rate chemistry effects are typically highlighted by means of temperature-
mixture fraction scattergrams, such as the one presented in Fig. 4. This
figure, referring to the transversal section x = 40 after 4, 000 time steps,
conveys a clear idea of the cooling effect produced by diffusive processes.
Part of these diffusive processes is likely to result from lattice discreteness
as well, and in particular from the finite-thickness Zth of the reaction region.
These numerical effects can be controlled by a careful tuning of the discretized
version of the Dirac’s delta function δ(Z − Zst ), but no systematic studies
have been conducted in this respect.
In Fig 5 we present the transverse temperature profile at three different
axial locations x = 1 (inlet), x = Lf /2 (mid-flame) and x = 128 (mid-
domain). As expected, a doubly-humped profile is observed within the flame
extension, with the two peaks centered about the stoichiometric surface (com-
pare with Figure 1). Beyond the flame front a typical gaussian profile is re-
covered. The overall structure of the temperature field is finally represented
in Fig 6, where the equicontours of the temperature field are displayed. From
this figure, the transition between double to single humped transverse profiles
is even more apparent.
In conclusion we have shown that minor modifications to the basic FCHC
lattice Boltzmann scheme permit to deal with reactive flows in the limit of
infinitely fast chemistry. Given the fact that chemistry is intrinsically local in
space and Lattice Boltzmann is a proved parallel performer on its own even
for inert flows, excellent performances on parallel machines may be easily

8
anticipated. This said, it is worth emphasizing that this work is still very
preliminary in character and several extensions are required before quan-
titative comparisons with state-of-the art numerical combustion techniques
can be attempted. Among others, important future developments are the
following ones:

1. Two-dimensional laminar to three-dimensional turbulent flames

2. Infinitely fast to finite-rate chemistry

3. Nearly incompressible to compressible flows

The first item does not require any modification to the scheme presented
here, but just higher numerical resolution. To convey a quantitative idea
of how far we can go with direct simulations, let us mention that a purely
hydrodynamic model (18s18p) for isothermal turbulent channel flows is cur-
rently achieving Reynolds numbers of about 3, 000 on the massively parallel
computer Quadrics, with a sustained speed over 10 Gflops [14]. The nu-
merical exploration of this regime should provide valuable informations on
the morphology of turbulent flames. Higher Reynolds numbers will call for
suitable turbulence models. These can easily be incorporated within the
LBE method so long as simple algebraic turbomodels such as Smagorinski,
Baldwin-Lomax and Renormalization Group models are used. All what is
needed is to let the scattering matrix Aij a local function of the stress tensor,
i.e. a pointwise function of the populations [5].
Moving beyond the fast chemistry regime requires the development of
a suitable multispecies transport scheme able to handle several scalars in
three dimensions. Since in the conventional Lattice Boltzmann formalism
each extra-scalar requires six moving populations, a naive extension of the
present scheme does not look very appealing. However it seems possible to
devise modified LBE schemes which are able to track passive scalars using
only population per species [18]. An alternative option is to couple the
present LBE scheme to ’conventional’ grid-based methods.
At variance with the previous two, the third item does involve substantial
changes to the scheme proposed in this work.
The scheme presented in this paper belongs to the simplest class of com-
bustion models; those in which density is not able to keep up dynamically
with significant temperature excursions. So, in principle, the method should

9
only be applied to ’cold’ flames with little heat release in which the main
focus is on flame morphology (i.e. the scenario pertaining to item 1. in the
aforementioned list of developments).
The next class of problems to be adressed is low-Mach reacting flows in
which density is allowed to respond to temperature changes over a signifi-
cant dynamical range of values. The possibility to model this class of flows
by means of appropriate phenomenological extensions of the present LBE
scheme is currently under investigation.
Finally, one should consider the class of fully-compressible, high-Mach
flows in which density is supposed to respond to both temperature and pres-
sure changes. We believe such class is definitely out of reach for the present
LBE model.
In fact, in order to capture full thermo-hydrodynamic behaviour the dis-
crete speeds must necessarily accomodate multienergy levels. A few models
in this class have already been proposed [15, 16], but they still need to be
perfected in many respects. First, they involve many discrete speeds (40
in 3D), let alone the additional degrees of freedom explicitly required to
incorporate passive/active scalar dynamics; second, and more importantly,
it appears like their stability is still matter of concern [17]. Further re-
search work is definitely required to clarify whether discrete speed models
(of reasonable complexity) will ever be able to perform competitively for full
thermo-hydrodynamic combustion applications.
Meanwhile, we believe that the model presented here, although clearly
limited in scope, can provide a useful starting point to explore phenomeno-
logical alternatives to fully-fledged discrete speed models for the numerical
simulation of (some class of) combustion phenomena.

References
[1] Lattice Gas Methods for PDE’s, G. Doolen ed., PhysicaD 47, 1-2, (1991).

[2] R. Benzi, S. Succi, M. Vergassola, Phys. Reports, 222(3), p.147, (1992).

[3] J. Sta. Phys., Special issue on Lattice models, Y. Qian, J. Lebowitz, S.

Orszag ed.s, vol 81,1-2, (1995).

10
 [4] Y.H. Qian, S. Succi, S. Orszag, Ann. Rev. Comp. Phys. III, p. 195,
(1995).

[5] J. Somers, P. Rem, Appl. Sci. Res., 48, 391, (1995).

[6] D. D’Humieres, P. Lallemand, U. Frisch, Europhys. Lett, 2, p.291,

(1986).

[7] S. Chen, Z. Wang, X. Shan, G. Doolen, J. Sta. Phys., 68, 3-4, 379,
(1992).

[8] Y.H. Qian, D. d’Humieres, P. Lallemand, Europhys. Lett., 17, 479,

(1992).

[9] F. Williams, Combustion theory, the Benjamin/Cummings, 2nd edition,

(1985).

[10] N. Peters, Comb. Sci. and Technology, 30, 1, (1983).

[11] A. Linan, Fluid Dyn. App. of Comb. theory, Pitman research notes in
Math. Sci., M. Onofri and A. Tesei ed.s, p.71, (1992).

[12] M. Smooke, in Numerical approaches to combustion modeling, E. Oran,

J. Boris eds, Progress in Astronautics and Aeronautics, vol. 135, Chap.
7, (1991).

[13] F. Williams, in ”The mathematics of combustion’, SIAM Frontiers in

Applied Mathematics, J. Buckmeister ed, 1985

[14] G. Amati, S. Succi, R. Benzi, Fluid Dyn. Res., 1996, in press

[15] F. Alexander, S. Chen, J. Sterling, Phys. Rev. E, 47, 2249, (1993).

[16] Y. Chen, Ph.D Thesis, Dept. of Quantum Engineering and System Sci-
ence, University of Tokyo, (1994).

[17] G. Mc Namara, A. Garcia, B. Alder, in Ref. 3, p.395.

[18] H. Chen, K. Molvig, C. Teixeira, private communication.

11
5 Figure captions
• Fig.1: Isocontours of the mixture fraction field Z(x, y) after 4, 000 time
steps. The stoichiometric line corresponds to Z = 0.057.

• Fig.2: The axial profile of the mixture fraction at the midline y = 128,
after 2, 000 and 4, 000 time steps. The solid line is associated with the
stoichiometric value Zst = 0.057.

• Fig.3: The axial profile of the temperature field at the midline y = 128,
after 2, 000 and 4, 000 time steps. The solid line is the initial profile and
the curve labeled ’Tadiabatic’ represents the adiabatic profile computed
with the mixture field Z at t = 4, 000.

• Fig.4: Temperature versus mixture fraction scatterplot at x = 40 at

time t = 4, 000..

• Fig.5: The transverse temperature profile at three axial locations, x = 1

(inlet), x = 25 (mid-flame) and x = 128 (mid-domain) at t = 4, 000..

• Fig.6: Isocontours of the temperature field after 4, 000 time steps.

12