Extended LES-PaSR model for simulation of turbulent combustion

fingersfieldΜηχανική

22 Φεβ 2014 (πριν από 3 χρόνια και 8 μήνες)

162 εμφανίσεις

4
TH

EUROPEAN CONFERENCE
FOR AEROSPACE SCIENC
ES (EUCASS)

Copyright


2011 by First Author

and
Second Author
. Published by the EUCASS association with permission.

E
xtended LES
-
PaSR model

for

simulation of turbulent
combustion


V. Sabelnikov


& C. Fureby

,


Fundamental and Applied Energetics Department, ONERA,

Chemin de la Hunière, 91761 Palaiseau Cedex, France


Defense Security Systems Technology, The Swedish Defens
e Research Agency


FOI,

SE 147 25. Tumba Stockholm, Sw
e
den



Abstract

In this work we propose a novel model for Large Eddy Simulations (LES) of pr
e
mixed, non
-
premixed and partially premixed turbulent flames. The develo
p
ment is motivated by the need for
m
ore accurate and versatile LES combustion models that can be used to simulate a range of
engineering applications such as airbreathing engines, IC
-
engines in order to assist in developing
improved combustion devices. The model is based on the f
i
nite chemis
try approach in which species
transport equations are solved with chemistry represented by a (typically reduced) reaction
mechanism. Provided that the filtered rea
c
tion rates are modeled sufficiently accurately this approach
has the adva
n
tage of being able

to represent also different types of flames under different conditions.
The mo
d
eling of the filtered reaction rate provides the challenge: As most of the chemical acti
v
ity,
and thus also most of the exothermicity occurs on the subgrid scales, this model n
eeds to be based on
the properties of fine
-
scale turbulence and mixing and A
r
rhenius chemistry. The model developed
here makes use of the similarities with the mathemat
i
cal treatment of multi
-
phase flows together with
the knowledge of fine
-
scale turbulence

and chemistry o
b
tained by Direct Numerical Simulation
(DNS) and experiments. In the model deve
l
oped, equations are proposed for the fine
-
structure
composition and vo
l
ume fraction that are solved together with the LES equations for the resolved
scales. If
su
b
grid convection can be neglected the proposed model simplifies to the well
-
known
Partially Stirred R
e
actor (PaSR) model. To validate the proposed LES model comparisons are made
with DNS data of a planar turbulent flame in homogenous isotropic turb
u
lence
.


1. Introduction

Predictive modeling of turbulent combustion is becoming increasingly important for the d
e
vel
opment of air
-
, land
-

and marine gas turbines, airbreathing engines, IC
-
engines and for studies of e.g. ammunition and waste destruction as
wel
l as fire spread.
Reynolds Averaged Navier
-
Stokes (RANS) based combustion models, [1], have been successful
in predicting some gross features of com
bus
tion, such as combustor exit temperatures, whereas they are unable to
pr
e
dict tran
sient phenomena such

as combustion insta
bilities, cycle
-
to
-
cycle vari
a
tions, self
-
ignition and emission
formation. Large Eddy Simu
lation (LES), [2
-
3], based combustion models is a promising alter
na
tive, [4
-
6], that can
provide accurate data for improvement of combustion d
evices. The idea of LES is to compute the large (energetic)
scales of the flow, affected by the boundary co
n
di
tions, whilst mod
eling the small (less energetic) scales of the flow.
For non
-
reactive flows this translat
es into providing closure models for
the subgrid transport terms, [2, 4], whilst for
re
ac
t
ing flows also the filtered mixing and (highly non
-
linear) reaction rate terms requires proper modeling, [4
-
6], as
the reaction zone is not often well r
e
solved on the grid.


To overcome these difficult
ies several methods have been developed, such as (i)
Flam
e
let Progress
Variable (LES
-
FPV) models
, including G
-
equation, [7], and flame surface de
n
sity, [8
-
9], models in which the
reactions are assumed to take place in thin layers, separating the reactant m
ixture from the product mixture, wrinkled
by the turbulence; (ii)
Probability Density Function

mo
d
els

(LES
-
PDF), [10], employing presumed shape or
transported probability density functions to model the filtered reaction rates; (iii)
Finite Rate Chemistry m
odels

(FRC
-
LES) in
cluding the Thic
k
ened Flame Model (LES
-
TFM), [11], Partially Stirred Reactor (LES
-
PaSR) model,
[12], as well as the Eddy Dissipation Concept (LES
-
EDC) model, [13], using different mathematical and
SESSION NUMBER & NAM
E


2

phenomenolog
i
cal models for the filtered

reaction rates; (iv)
Conditional Moment Closure

(LES
-
CMC) mo
d
els, [14],
in which the species equations are conditionally averaged on a few variables on which the reaction rates are critically
known to depend on, and (v)
Linear Eddy Models

(LES
-
LEM), [15],

using a grid
-
within
-
the
-
grid approach to solve
1D species equations with full resol
u
tion. All these models have features that limit their usefulness in different ways,
and thus it is d
e
sirable to develop more vers
a
tile LES combustion models.


This paper p
roposes a novel LES finite rate chemistry model for turb
u
lent combustion. The idea of this
model is based on the ideas of Chomiak, [16
-
17], which in turn are based on the non
-
reactive experimental data

of
Batchelor & Townsend, [18], indicating that
at high

Re
y
nolds (Re) numbers turbulent fine structure is not uniformly
distributed but concentrated into small is
o
lated regions, whose entire volume is a small fraction of the total vo
l
ume.
This has more recently been verified by Direct Numerical Simulations (DN
S) in non
-
reactive, [19], as well as
rea
c
tive, [20], flows. Later, Magnussen, [13, 21], starting from the Chomiak model, developed the well
-
known EDC
model.
The core of the model proposed here is the use of mathematical treatment of multi
-
phase flows for t
he
modeling of the high
-
intensity fine
-
scale structures embedded in an environ
-
ment of low
-
intensity turbulence, with
slower mixing. The proposed model framework incorp
o
rates the PaSR model, [12], as limiting case when subgrid
convection is n
e
glected, and
thus the model will be referred to as the Expended PaSR (EPaSR) model.
Here, we
develop the EPaSR model and test it against in
-
house DNS data of a premixed planar flame in homogeneous
is
o
tropic turbulence
.

To compare the proposed model with conventional LE
S combustion model compar
i
sons are
also shown with a few ordinary LES
-
FPV and LES
-
FRC models.
The LES
-
EPaSR (and PaSR) model reproduces the
experimental and DNS data well and typically shows improved agre
e
ment with the reference data than the
conventional
LES
-
FPV and FRC mo
d
els.


1.1
A conceptual representation of turbulent combustion


To develop successful physically
-
based subgrid Turbulence Chemistry Interaction (TCI) mo
d
els for LES
-
FRC we
have to take into account that the spatial and temporal distributi
ons of velocity gradients have a major impact on the
spatial and temporal structure of the chemical reactions r
e
gions in turbulent flows. As was observed first
(in non
-
reactive case)

by
Batchelor & Tow
n
send, [18
], velo
c
ity gradients are highly intermittent
, which means that the fine
structures, possessing most of the small
-
scale turbulence, are distributed non
-
uniformly in space, with the
intermittency increasing with increasing Re
y
nolds (Re) number. This picture of small
-
scale turbulence is quite
different

from the Kolmog
orov 1941
turbulences theory, [22
-

23
]. In this K41 turbulence theory, the fine
-
structures,
and associated quantities such as dissipation and acceleration, are distri
b
uted quasi
-
uniformly in space and in time.
The sp
a
tial and temporal stat
istics of fine
-
structures in the K41 theory is determined solely by the Kolmogorov
length and time scales, i.e.
4
/
1
3
)
/
(



K


and
2
/
1
)
/
(




K
, respectively, [22
-
23
], in which


is the mean
dissipation. The experimental study of Kuo &
Corrsin, [2
4
], suggests that fine
-
structure regions consist of
topological stru
c
tures such as vortex sheets, ri
b
bons and tubes folded in regions of the fl
ow. More recently, DNS,
e.g. [19
-
20
], provided support to this view, revealing how high intensity vort
ices merge in co
m
plex shaped filaments
that are embedded into sheets or arcs of low(er) intensity vo
r
ticity.


The intermittency is also apparent by long
-
range spatio
-
temporal correlations (co
m
parable with integral
length and time scales,
I


and
I

, respectively) of quantities such as dissipation and acceler
a
tion. Kolmogorov,
[25], and Oboukhov, [26
], proposed 1962 a modified theory of small
-
scale turbulence, where intermi
t
tency is taken
into account, following the w
ell
-
known r
e
mar
k of Landau, [27
]. The most exciting result of the K62 theory is that the
spatial and temporal st
a
tistics of the fine
-
stru
c
tures is determined not solely by
K


and
K

, as was supposed in the
K41 theo
ry, but also by
I


and
I

. This interrelation of the Kolmogov and larger scales (up to integral scales) can
be illustrated by considering a tubular region in a turbulent flow with length
I


and

diameter
K


(<<
I

). This
tubular region is thus subjected to the stretch produced by the velocity grad
i
ent of larger vortices,
I
v

/

, where
v


represents the fluctuating velo
city. Since the velocity is correlated to a length of order
I

, it may be assumed that
the characteristic time scale of this stretch is of order
v
I

/

, being applied over a length of order
I

. Thi
s situation
is similar to Burgers vo
r
tic
es, and time of life of this tube will be of the order of integral time scale
v
I

/

. A more
d
e
tailed discussion of inte
r
mittency in high Re turbulent flows may

be found in Monin & Yaglom, [28], Frisc
h, [29
],
Ts
i
nober, [3
0], and other review papers, [31
-
32
].


Recent DNS of combustion (a
t moderate Re numbers), e.g. [20
], show a similar organiz
a
tion of the flow
with the fine
-
structure vortices at the flame being essentially parallel to the flame whereas
those behind the flame are
predominantly perpendicular to the flame. Regions of high exothermicity and volumetric expa
n
sion are found in
tubular stru
c
tures distributed among the fine
-
structure vortices whereas regions of low exothermicity may be
distribute
d more randoml
y. [16
-
17
]. Thus we can conclude that the geometry of heat release zones (more generally,
First Author, Second Author
.
INSTRUCTIONS FOR THE

PREPARATION OF PAPER
S



3

zones, where the chemistry reactions take place) and their spatial and temporal distributions are highly correlated
with fine turbulent structures. From

physical viewpoint this conclusion is not surpri
s
ing as these small
-
scale
structures contain the molecular fluxes on which the mixing and chem
i
cal reactions, are so heavily depending on.
The consequences of the last conclusion are very i
m
portant for the f
ormulation of TCI models. If the spectrum of
life
-
times of the fine structures can vary from the Kolmogorov scale to the integral scale the micromixing times that
inevitably must e
n
ter any subgrid TCI model have also to respect this wide range of scales. T
he proposed model that
will be described next is formulated upon the cartoon sketched above.


2
.
LES combustion modeling using a multi
-
phase analogy

The reactive flow equations are the bal
ance equations of mass, momentum and en
ergy describing convection,

diffusion and chemical reactions, [4
-
6, 10].
In LES, [2], all variables are deco
m
pos
ed into resolved and unresolved
(subgrid) components by a spatial filter such that
f
f
f




~
, where


/
~
f
f


is the Favre filtered variable.
Filtering the r
e
active flow equations yields,


















































).
(
~
)
(
~
)
(
~
)
~
~
(
)
~
(
),
(
)
~
~
(
)
~
(
,
)
(
)
~
~
(
)
~
(
,
0
)
~
(
)
(
,
1
















i
f
i
N
i
h
t
s
s
t
t
i
i
i
i
i
t
t
h
w
p
p
h
h
p
w
Y
Y


b
h
v
v
S
v
B
S
v
v
v
b
j
v
v

(1)


Here,


is the density,
v

the velocity, p the pressure,
S

the viscous stress tensor,



T
T
i
p
i
i
s
dT
C
Y
h
0
,

the sensible
ent
halpy, T the temperature,
h

the heat flux ve
c
tor, Y
i

the sp
e
cies mass
-
fraction,
i
w


the species re
action rate,
j
i

th
e
species mass
-
flux,

i
f
h
,

the species form
a
tion enthalpies and


the radiative heat lo
ss which will be neglected.
The
subgrid flow phys
ics is concealed in the su
b
grid stress tensor
)
~
~
(
~
v
v
v
v
B






and flux vectors
)
~
~
(
~
i
i
i
Y
Y
v
v
b




and
)
~
~
(
~
h
h
h
v
v
b



, which r
e
sults from filtering the convective terms.
Following [6] we
post
u
late that the gas mi
x
ture b
e
haves as a linear viscous fluid with Fourier heat conduction and Fickian species
diff
u
sion so that
i
i
i
Y
D
~


j
,
T
R
p
~


,
D
D
S
~
2


,
T
~



h
, respectively. Here,
i
i
Sc
D
/



is the
species diffusivities, R the composition dependent gas constant,


the viscosity,
I
v
v
v
D
)
~
(
)
~
~
(
~
3
1
2
1







T
D

the deviatoric part of the rate of strain tensor,
Pr
/




the thermal diffusivity, in which Sc
i

and Pr are the
Schmidt and Prandtl numbers, respe
c
tively. Finally, the filter
ed reactions rates are de
fined as
j
ij
M
j
i
i
w
P
M
w


1



,
in which M
i

is the molar mass of specie i, P
ij

the stoichiometric matrix and
ij
j
ij
N
i
P
i
N
i
j
P
j
Y
k
T
w
1
1









the
reaction rate of the j
th

reaction step, with

k
j

b
e
ing the Arrhenius rate.

2.1

Subgrid flow model
ing

To close the LES equations (1) we need to provide models for
B
,
b
i
,
b
E

and
i
w

, and considering first the subgrid
stress tensor and flux ve
c
tors, these are not unique to reactive flows and closure models can be acquired from the
ple
thora of su
b
grid models for non
-
reactive flows, [2
-
3].
Here,
we use the Mixed Model (MM), [33
], with
D
k
D
v
v
v
v
B
~
2
)
~
~
~
~
~
~
(
~







,
i
Sc
i
i
i
Y
Y
Y
t
k
~
)
~
~
~
~
~
~
(
~






v
v
b

and
s
s
s
h
h
h
h
t
k
~
)
~
~
~
~
~
~
(
Pr
~






v
v
b
, in which
the subgrid viscosity,
2
/
1
k
c
k
k



, is provided by an eq
uation for the

sub
grid kinetic energy, k, [33
], and where
Sc
t

and Pr
t

are the turbulent Schmidt and Prandtl numbers. The model coefficients are calculated u
n
der the
assum
p
tion of an infinite inertial sub
-
range, [2]. To reduce the computational cost we use

wall
-
modeled

LES, in
which a model is used to handle

the near
-
wall flow physics, [34
]. This model is based on r
e
pla
c
ing the viscosity
SESSION NUMBER & NAM
E


4

(

+

k
), and thermal and species diffusivities (

/Pr+

k
/Pr
T

and

/Sc
i
+

k
/Sc
T
, respe
c
tively) in the first grid point at the

wall, with va
l
ues co
n
sistent with the log
a
rithmic law of the wall.

2
.2 A multi
-
phase framework for subgrid combustion modeling

Considering next the filtered reaction rates
i
w

, incorporating the turbulence chemistry intera
c
tions, being

notoriously difficult to model because of their non
-
linear nature, [1]. Using the ph
e
nomenological models, [16
-
17,
35
], based on experim
ental data, and DNS results, [19
-
20
], a ca
r
toon of turbulent mixing and combustion may be
drawn. In this, turbulent r
e
a
cting flows may be viewed as a muddle of vortex structures of different topological
character, sheets, ri
b
bons and tubes, in which the tubes and ribbons carry most of the high
-
intensity vorticity and
dissipation. This implies that the fine
-
structure region
s, denoted by (
*
), are embedded in a su
r
rounding fluid, here
denoted by (
0
), will be responsible for most of the molecular mixing, chem
i
cal reactions (if the temperature is
sufficiently high) and heat release. This cartoon is sim
i
lar to that of the Eddy Di
ssipation Concept (EDC), [
13, 21,
36
], but here a different strategy, based on an analogy with mode
l
ing multiphase flows, will be used to model the
filtered reaction rates
i
w

.


For the derivation of the proposed model we let
]
,
[
s
i
h
Y

Ψ

be the composition space that satisfies the
equation
i
i
i
i
t













k
v
)
(
)
(
, in which
i
i
Y


,
i
i
j
k


and
i
i
w





for 1<i<N, where N
denotes the number of species in the species equation (1
3
)
whereas


N

1

h
s
,
h
k


1
N
, and
)
(
~
~
,
1
1



i
f
i
N
i
t
N
h
w
p
p













v
v
S

are the terms associated with the e
n
ergy equation (1
4
). As for the
theor
y of multi
-
phase flows, e.g. [37
], a phase indicator function, d
e
noted

I

, wi
th

=1 in the fine structures (*) and

=0 in the surroundings (
0
) is next i
n
troduced to di
f
ferentiate betwe
en these regions so that the local balance
equations for the comp
o
sition space becomes,




,
)
(
)
(
)
(
,
,
,
,
,

















i
i
i
i
i
t
I
I
I
I










k
v



(2)


in which

,
i


denote the exchange terms at the immaterial surface separating the fine structures and

=0 in
the
surroundings. By summing over


in (2), provided that the exchange terms sati
s
fies
0
)
(
,
1
0






i
, the
comp
o
sition space equations are recovered.

Although achievable, velo
c
ity differences between the fine
-
structures
and su
r
roundings are not included. In the frame
work of LES the species mass fraction equations are filtered over
space to obtain equ
a
tions for the large energetic scales of motion. If
we take
the filtering to correspond to a box
-
filter, cover
ing the cell volume, ∆V, as often done in LES, [2], the vo
l
ume fraction of phase


is,




,
1
with
/
)
,
(
1
0
1



















I
V
V
dV
t
I
V
V
x



(3)


so that for products



,
i
I

and





I
i
,

it can easi
ly be shown that







I
I
i
i
/
)
(
,
,


an
d






I
i
,















)
~
(
)
(
)
(
,
,
i
i

. Filtering the equations (2) then yields,




,
)
(
)
)
(
(
)
~
)
~
(
)
(
(
)
)
~
(
)
(
(
,
,
,
,
,





























i
i
i
i
i
i
t











b
k
v
(4)


in which
v
v
b
~
)
~
(
)
(
,
,















i
i
i
I



denotes the subgrid transport term, which will typ
i
cally b
e
modeled by a conventional subgrid viscosity model. Summing over the filtered equ
a
tions (4) results in that the
filtered exchange terms satisfy
0
)
(
,
1
0






i

so that the LES composition space equations (1
2
) and (1
4
) are
recovered. Returning next

to the old notions, with (*) d
e
noting the fine structures and (
0
) the surroundings, such that

*
i


1
1
,
)
~
(
i


and
2
2
,
0
)
~
(
i
i



, and



0
0
*
*
~
i
i
i






0
*
*
)
1
(
*
i
i






, we then obtain the set of
equat
ions,

First Author, Second Author
.
INSTRUCTIONS FOR THE

PREPARATION OF PAPER
S



5



























.
))
(
(
)
~
(
)
(
,
))
(
(
)
~
(
)
(
0
0
0
0
0
0
0
0
0
0
*
*
*
*
*
*
*
*
*
*
i
i
i
i
i
i
t
i
i
i
i
i
i
t






















b
k
v
b
k
v


(5)


The above equations can be solved together for


i
*

and


i
0

after which the LES quantity
i

~

can be constructed
from
0
*
*
)
1
(
~
*
i
i
i









once the fine structu
re volume fraction,
*

, is known.
A more convenient and
straightforward approach is however to solve the balance equ
a
tions for the fine structure fractions,
*
i

, (5
1
) together
with the LES balance equations for
i

~
, or more pr
e
cise
ly the LES species equ
a
tions (1
2
) and energy equation (1
4
).


Noticing that
0
*
*
*
)
1
(
i
i
i











, and more specifically that
0
*
*
*
)
1
(
i
i
i
w
w
w








, and taking
a
dvantage of the theoretical, [16
-
17], experimental, [35
], an
d

computational, [20
], o
b
servations that (fast exothermal)
reactions mainly occur in the fine structures, (1
2
) becomes,




,
)
(
)
~
~
(
)
~
(
*
*
i
i
i
i
i
t
w
Y
Y













b
j
v





(6)


emphasizing the importance of the volume fraction of the fine structures. Equation (6) further implie
s that the
reaction rate c
an be expressed as
*
0
0
*
*
*
)
1
(
)
(
)
(
i
i
i
i
i
w
w
w
d
w
w
















Ψ
Ψ
Ψ
, in which
)
(
Ψ


is
hence a probability de
n
sity function of the form
)
(
)
1
(
)
(
)
(
0
*
*
*
Ψ
Ψ
Ψ
Ψ
Ψ










. Moreover, by
summing over the N species equations in (5) and taking specific
ally into a
c
count that
1
)
(
)
(
0
1
1
*






i
N
i
i
N
i



and
0
)
(
)
(
0
1
*
1






i
N
i
i
N
i




, it immediately fo
l
lows that,




),
(
)
~
(
)
(

and

)
(
)
~
(
)
(
0
1
0
0
0
0
*
1
*
*
*
*
i
N
i
t
i
N
i
t














v
v











(7)


the sum of which must satisfy the continuity equation (1
1
). This constraint subsequently results in that
)
(
)
(
0
*
i
i
i
i
m








, where we denoted by
m


the exchange rate of mass between the fine structures and the
surroundings. This exchange rate of mass between the fine structures and su
r
roundings is a key quantity for accurate
predictions of species m
ixing and reactions.

2.3

Subgrid combustion modeling

To provide a closure of the exchange rate of mass,
m

, through the immaterial surface separating the fine structures
and surroundings, we suppose that if there is a dynamic equilibr
ium state b
e
tween the fine structures and the
su
r
roundings, which is characterized by the value
*
eq


(which is to be discussed later), then mass exchange rate
between the fine structures and the su
r
roundings is equal to zero, i.e.
0
)
(
*

eq
m


. For this equilibrium state, the
surface separating the surroun
d
ings and fine structures can be considered material. Developing
m


in the vicinity of
*
eq

, whilst retaining only the first
-
order li
n
e
ar term, assumed proportional to
)
(
*
*
eq



, results in that,




,
/
)
(
*
*
*




eq
m











(8)


where
*


is the fine structure residence time, and where the term
*
/



follows from dim
e
n
sional consid
erations.
The hypothesis underlying (8) refers to high Re
-
number flows, in which combustion takes place in
concentrated
regions, whose entire volume is a small fraction of the total vo
l
ume. The application of this approach for lower Re
-
numbers, i.e. in the

flamelet r
e
gime, is not straightforward.
We conclude from (8), that if
0

m

,
*
*
eq



, then
the exchange rate of mass is directed from the surroundings to the fine structures. Otherwise, if
0

m

,
*
*
eq



,
then the exchange rate of mass is directed from the fine structures to the surroun
d
ings. We also note, that this mass
SESSION NUMBER & NAM
E


6

exchange rate is due to convection through the interface separating fine structures and surroundings. Using the
closure

model (8),
*


can be o
b
tained from (7
1
),




.
/
)
(
)
~
(
)
(
*
*
*
*
*
*
*









eq
t






v





(9)


Next we consider the necessary
modelling

of the exchange terms
*
i


and
0
i


in equation (
5). As follows from the
core
physical co
n
siderations,
*
i


and
0
i


contain two kinds of terms. The first type of term, here denoted by
*
i


and
0
i

, is due to the exchange rate of mass between the fine struc
tures and the surroundings (through the surface
separating the fine structures and surroun
d
ings) as discussed above. If the exchange rate of mass is absent (i.e. in the
dynamic equilibrium state, when the surface separating the fine structures and surround
ings can be considered as
m
a
terial one),
,
0

m


as there is no transport through the interface between the fine structures and the surroundings,
these two terms both turn to zero. The second type of term, here d
e
noted by
*
i


and
0
i

, is due to molecular
diffusion through the interface between the fine structures and surroundings. Indeed, even if the exchange rate of
mass is absent (i.e. in dynamic equ
i
librium), i.e.
,
0

m


there is ex
change through the inte
r
face due to molecular
diffusion. If
0
*
i
i



, these terms however turn to zero, since in this case there is no molecular diffusion. Based on
the above physical considerations we pr
o
pose to model these exchange terms a
s,




,
,
)
(
)
(
,
/
)
(
,
/
)
(
,
,
*
0
2
1
0
2
1
*
0
*
*
0
0
*
*
0
0
0
*
*
*
*
*
*
*
*
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
m
m
m
m















































(10)


thereby satisfying
0
)
(
,
1
0






i

(in fact, separately for both two types of terms,
0
)
(
,
1
0






i

and
0
)
(
,
1
0






i
, as has to be). We note, that in correspondence with the underlying physics
, if mass exchange is
directed from the surroun
d
ings to the fine structures,
0

m

, one obtain from (10) that
0
*
i
i
m





and
0
*
0
i
i
i
m








. Similarly, if the mass exchange is directed from the fine structures to the

su
r
roundings, i.e.
for
0

m

, one has that
*
*
i
i
m





and
*
*
0
i
i
i
m








.


By assuming that the convective and diffusive terms in (5
1
) can be neglected we have that
0
*
*
*



i
i



. Summing in this relation o
ver i

[1,…, N], and taking into account that
0
)
(
*


i
i


, as follows
from the law of mass action, it follows that
0
)
(
*



i
i
, which, by the definition of
m


as
)
(
*
i
i
m




, results
in that
0

m

. By using the linear model introduced in (8) this results in that
*
*
eq



. The fact that
0

m


when
convection and diffusion can be n
e
glected at the subgrid level in (5
1
) further i
m
plies that
0
*


i
, which in turn
implies that
*
*
/
)
(
0
*
*
*





i
i
i
i






. From the relation
0
*
*
*



i
i




it then follows that
*
*
/
)
(
0
*





i
i
i



, which may be further rearranged, u
s
ing the afor
e
mentioned definition
0
*
*
)
1
(
~
*
i
i
i








, to the more appropria
te form,




,
)
1
/(
/
)
~
(
*
*
*
*
i
i
i
w
















(11)


which can be recognized as the standard LES
-
PaSR model, [12]. Note also that the EDC model, [13, 21], bears some
resemblance of
the algebraic model (11). Equation (11) can be a
r
rived at and understood also fr
om simple mass
balance considerations as done in the conve
n
tional PaSR and EDC mo
d
els but is here (i) put on a more solid
mathematical foundation and (ii) extended to non
-
equilibrium situations prevailing in some combustion applic
a
tions.


The framework dev
eloped so far can be used to provide a more rigo
r
ous foundation for the LES
-
PaSR
co
m
bustion model, [12
], and to propose an extension of this model, referred to as the Extended PaSR (EPaSR) LES
First Author, Second Author
.
INSTRUCTIONS FOR THE

PREPARATION OF PAPER
S



7

m
odel. In the LES
-
PaSR model, [12
], the continuity, m
o
mentum an
d energy equations, (1
1
), (1
3
) and (1
4
),
respectively, are solved togeth
er with the species equ
a
ti
-
on (6), in which the filtered reaction rates are modeled by
)
(
*
*
*
*
Ψ
i
i
i
w
w
w







, in which
*
i


r
e
sults from solving (11) in which
*


is provided by
*
eq

, yet to be
specified. The PaSR
-
LES model has been proven
reliable and accurate, [12, 38
], but a better model is needed for
more complex flows. In the proposed LES
-
EPaSR model the con
ti
nuity, mom
entum and energy equ
a
tions, (1
1
), (1
3
)
and (1
4
), respectively, are solved simultaneously with the species equ
a
tion (6) in which
)
(
*
*
Ψ
i
i
w
w




, where
*
i


results from solving the transport equations (5
1
) with source terms
*
i
M
, provided by (10
1
) and


*

provided by
solving (9). To avoid calculating the fine
-
structure de
n
sity,
*

, equations (5
1
) and (9) will here be simplified by
replacing
*


by


so that in the LES
-
EPaSR model the fine structure composition and enthalpy (or temper
a
ture) is
given by,



























,
/
)
(
)
~
(
)
(
,
)
1
/(
/
)
~
(
))
(
(
)
~
(
)
(
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*


























eq
t
i
i
i
i
i
i
i
i
t
m
v
b
k
v



(12)


for 1≤i≤N+1, in which
0
*
*
)
1
(
~
*
i
i
i









has been used to simplify the rhs of (12
1
) and in wh
ich
*
/
)
(
*
*




eq
m




.
To close (12) the value of
*
eq


is required which will be discussed next.

2.4

Modeling the subgrid time scale and the equilibrium reacting volume fraction

To finally close the EPaSR and PaSR models, the subgr
id time,
*

, and reacting volume fra
c
tion,
*

, needs to be
provided. In what follows we propose sub
-
models for both these quantities based on the physical description
provided in Section 2, and the references cited
.


To estimate the subgrid time,
*

, we first observe that the fine scale structures in a volume,
*
V

, are
generally anisotropic (sheets, ribbons and tubes) and influenced by the subgrid velocity stretch,


/
v
, the
characteristic time of which is
v


/
. The smallest quasi
-
equilibrium scale of the fine structures,
D

, is controlled
by molecular viscosity and the characteristic time of the stretch


/

v
, such that
2
/
1
)
/
(
v
D





. Its similarity
with the Taylor length scale,
2
/
1
)
/
(
v
I
T





, defin
e
d from the molecular viscosity and the integral time scale
v
I

/

, su
g
gests that
D


is the dissip
a
tiv
e length scale of the smallest r
e
solved scales. We further assume that area
-
volume ratio,
*
*
/
V
S


, is solely determined by
D


so that
D
V
S

/
1
/
*
*



. For small deviations from the
equilibrium state, a simple rela
xation equation can be formulated for the vari
a
tion of
*
V


so that,




),
(
)
/
(
/
)
(
*
*
*
*
*
eq
eq
K
V
V
V
S
v
dt
V
d












(13)


in which we take into account that the change of the volume
*
V


is due to the propag
a
tion of the surface
*
S


with
relative velocity on the order of the Kolmogorov velocity,
K
v
, [27
]. From (13) it also becomes evident that the
subgrid time scale may be expressed as,





.
/
)
)
/
(
(
1
*
*
*
K
D
eq
K
v
V
S
v













(14)


Note also that (13) a
nd (14) degenerates into (9) when divided by ∆V. It should however also be noted that if it
assumed that a
rea
-
volume density,
*
*
/
V
S


, is instead determined by the Ko
l
mo
gorov scale,
K

, we obtain
K



*
, which contradicts the physical description of the fine stru
c
tur
es di
s
cussed above. Using the definition of
the Kolmogorov length
-

and time scales,
4
/
1
3
)
/
(



K


and
K
K
K
v
/



, re
spectively, this finally r
e
sults in
SESSION NUMBER & NAM
E


8

that





K
*

in which
v




/


is the time of the subgrid velocity stretch. This represents the ge
o
metrical
mean of the (short) Kolmogorov time and the (long) time scale associated with the subgrid v
e
locity stretch.


To model the equilibrium r
eacting fine
-
structure fraction,
*
eq

, we use dimensional con
side
-
rations and
analysis of two limiting cases (fast and slow chemistry in comparison with mi
x
ing). More precisely, following from
the dimensional considera
tions we write,




).
/
,
/
,
/
(
*
K
K
m
c
eq
f















(15)


Since we consider the equilibrium state (
0

m

), the parameters that characterize un
steady and convective ef
fects
are not present in (15). Let us for simplicity denote the charac
teristic time and len
gth scales ratios by
*
/


c
x

,
K
y


/



and
K
z

/


, respectively, and by neglecting intermittency we have that
3
/
2
z
y

. Thus we
consider further that
*
eq


depends solely of x

and y,




).
,
(
*
y
x
f
eq











(16)


There are two limiting cases for th
e parameter y, namely y>>1 and y≈1. The first case is indeed LES since from
above it follows
that
K



, whereas the second case is indeed DNS since all scales are r
e
solved, i.e.
K





and
K



, and th
us
1
)
,
(
*


y
x
f
eq


irrespectively of x. For the first case we consider two limiting
situations: (i) x<<1 (or
*



c
), which corresponds to fast chemistry, in which combustion is limited by mixing,
and (ii) x>>1 (or
*



c
), which corresponds to slow chemistry, in which combustion is limited by the chemical
reaction rate. In the first situation (x<<1 or
*



c
), when combustion is limited by mixing, the pro
d
uct
*
*
i
eq




has to be of th
e order
*
/
1

, and if we take into account that
c
i


/
1
*


,




.
1


1

,
)
,
(
*




y
and
x
x
y
x
f
eq






(17)



In the second situation (x>>1 or
*



c
), when combustion is limited by the reaction rate, the prod
uct
*
*
i
eq




has to be of the order
c

/
1
.
Thus, if we take into account that
c
i


/
1
*


,




.
1

,
1

,
1
)
,
(
*




y
x
y
x
f
eq







(18)


In order to provide a reliable model for
*
eq

, respecting the underlying
cases and situations di
s
cussed abov
e we
propose to model
*
eq


in the following manner,


),
/
(
)]
/(
|
))
/
(
1
(
)
(
)]
1
/(
[
))
(
1
(
)
,
(
*
*
K
c
c
K
eq
F
F
y
F
x
x
y
F
y
x
f




















(19)


where
)
(
y
F

is a monotonous function satisfying
0
)
1
(

F

and
1
)
(


F
.
In this paper solely the case

y

1
,
(
K




) is considered resulting in that (19) reduces to,




.
)
/(
)
1
/(
)
,
(
*
*









c
c
eq
x
x
y
x
f





(20)


This model for the reacting volume fraction is identical to the phenomenologically derive
d mod
el successfully used
in previous LES, in which the chemical time scale is estimated as


u
u
c
s
/



2
/
u
s


and the subgrid time scale
First Author, Second Author
.
INSTRUCTIONS FOR THE

PREPARATION OF PAPER
S



9

is estimated as
K





*
. Introducing
v




/


and
2
/
1
)
/
(




K
, in which



/
)
(
3
v

,
(20)
becomes
)
/(
)
(
4
/
5
4
/
3
4
/
3
2
4
/
5
4
/
3
*
v
s
v
u
eq








.

2.5

Relations with other turbulence chemistry interaction models

The proposed LES
-
EPaSR model is next placed in perspective by comparison with contemporary LES combustion
models such as e.g. the LES
-
TFM, EDC, PaSR and FPV models.
The co
m
parison is performed by first assessing the
theoretical formulations of these models and secondly by comparing predictions of these models with reference data,
as will be discussed in Sectio
ns 7 and 8. In addition to the continuity, momentum and energy equations the LES
-
TFM, EDC, PaSR and EPaSR models rely on species equations of the form
i
i
k
i
i
i
t
w
Y
D
D
Y
Y










)
~
)
((
)
~
~
(
)
~
(
v



, whereas the LES
-
FPV family of models makes use of an
equation for a progress vari
able, c, (and a mixture fraction equation if appropriate) such that
c
k
t
w
c
D
c
c









)
~
(
)
~
~
(
)
~
(
v



, in which D
k

is the subgrid diffusivity, and
i
w


and
c
w


the reaction rate
terms. Table 1 details D
k

and reaction rate terms fo
r the models tested here. A
l
though many LES
-
FPV models are
available, we have here used the model of Boger
et al
., [39], as a representative candidate. In these models F and E
represent the flame thickening and eff
i
ciency factors, respectively,
*
PaSR

,
*
EPaSR


and
*
EDC


are the reacting
volume fractions for the PaSR, EPaSR and EDC models, respectively, whereas
*
k
Y

are
*
T

are the fine
-
structure
species mass fractions and t
emperatures. For further details we refer to the references and additional re
f
erences
therein. Regarding the expressions for the reacting vo
l
ume fractions both PaSR and EPaSR use the same model,


PaSR
*


eq
*


c
/(

c


*
)
, as introduced ea
r
lier, which may be r
eform
u
lated as


PaSR
*


eq
*

(

3
/
4

v
5
/
4
)
/(
s
u
2

3
/
4


3
/
4

v
5
/
4
)
. For the EDC model
3
*
*
)
/
(
v
v
EDC



, which may, according to [6],
be reformulated as
4
/
3
*
)
/
(
02
.
1
v
EDC





. In the LES
-
TFM model an expre
s
sion similar to the reacting volume
fraction, E/F, occurs, a
l
though havin
g a slightly different i
n
terpretation, which may, according to the approximation
in [11] be reform
u
lated as
F
F
E
/
/



, in which
)
/
,
/
min(
1
u
u
s
v









is the resolved flame
wrinkling and


the efficiency fun
c
tion which may be approximated as
3
/
2
3
.
0
)
/
)(
)
/
/(
2
.
1
exp(
75
.
0
u
u
s
v






, [11]. To compare the models in Table 1 we apply an analogy with
the laminar flame speed a
p
proximation
2
/

w
D
s
L



such that
2
/
)
(

w
D
D
s
k
T



, which typically
re
sult in that
L
T
s
s


, where


is the flame wrinkling factor as included in Table 1 for the different models. This
comparison can only be approximate due to the necessary assum
p
tions, and a more co
m
prehensive comparison is
thus provided in Se
ctions 7 and 8 where predictions using the models listed in Table 1 are compared for (i) a planar
flame propagating in h
o
mogeneous isotropic tu
r
bulence for which in
-
house DNS data is available for detailed
compar
i
son and (ii) for an ax
i
symmetric dump combu
stor for which high
-
quality experimental data is available for
compar
i
son.


Table 1.
LES combustion model comparison


Model


Diffusivity

Reaction rate models

Flame wrinkling

TFM, [11]

i
D
FE
)
1
(


)
~
,
~
,
(
T
Y
w
k
i
F
E



)
/
,
/
min(
1
u
u
s
v






P
aSR, [12]

T
Sc
/


)
,
,
~
,
~
,
(
*
*
*
*
T
Y
T
Y
w
k
k
i
PaSR




)
/
1
(
*
D
D
k
PaSR



EPaSR

T
Sc
/


)
,
,
~
,
~
,
(
*
*
*
*
T
Y
T
Y
w
k
k
i
eq




)
/
1
(
*
D
D
k
eq



EDC, [6
, 13
]

T
Sc
/


)
,
,
~
,
~
,
(
*
*
*
*
T
Y
T
Y
w
k
k
i
EDC




)
/
1
(
*
D
D
k
EDC



FPV, [39
]

T
Sc
/






/
)
1
(
4
c
c
s
s
u
u
u
u



|)
|
/(
)
1
(
4
|
|
/
c
c
c
c











SESSION NUMBER & NAM
E


10

3.

Global hydrocarbon combustion chemistry


Accurate predictions of temperatures and species concentrations require the use of carefully s
e
lected reduced
schemes, and here reduced three
-
step

mechanisms for C
n
H
m
-
air combustion are used.
The mechanis
ms are detailed in
Table 2 for CH
4
, for which the Arrhenius coefficients have been o
p
timiz
ed to match the laminar flame speeds and
tem
pera
tures of the Gri
-
Mech 3.0, in lean pr
e
mixed flames, [40].

The third reaction step, modeling NO formation, is
the sum of two rates, one for thermal NO (only dependent on temper
a
ture) and one for prompt NO (depend
ing on the
hydrocarbon fuel and H
2
O). Most NO
X

reacts to NO
2

after combustion, however, close to the
flame, NO is dominant,
and therefore, it is suffi
cient to limit the mechanism to NO.
Assoc
i
ated with the reaction mechanism is the treatment
of molecular diffusion, that is here simpl
i
fied and based on species Schmidt numbers so that
i
i
i
Y
Sc
~
)
/
(



j
,
where the values of the individ
ual species Schmidt numbers
i
Sc

are

:
CH
4



0.68

Sc
, C
3
H
8



0.72, O
2



0.76,
CO


0.76, CO
2



0.98, H
2
O


0.60, N
2



0.75
, respectively
. For a critical evaluation and discussion of this

approach
to model the molecular diffusion we refer to Giacomazzi
et al
, [41].


Table 2. Rate parameters for the reduced three
-
step CH
4
-
air mech
a
nisms


No

Reaction

A [m,kg,mol,s]

n
CH4

n
O2

n
CO

T
a
[K]

1

CH
4
+1.5O
2


C伫2H
2
O

7.20∙10
11

0.9

1.1


17387

2

CO+0.5

O
2
CO
2

4.00∙10
8

0

0.5

1.0

6047

3

O
2
+N
2

2乏

1.50·10
16

-
0.3



38440




4.

Numerical method


For this study we use the C++ library OpenFoam, [42], as the computational platf
H
orm, which has previously been
used fo
r applications of varying complexity.
The code employs an unstructured collocate
d Finite Volume (FV)
method, [43
], in which the dis
cre
ti
zation is based on Gauss the
o
rem together with a semi
-
implicit time
-
inte
gration
scheme. Given the vector of unknowns
,

u

T
i
E
Y
]
~

,
~

,
~

,
[




v
, the semi
-
discretized equ
a
tions are,




),
,
(
)]
,
(
)
(
)
(
[
)
(
1
u
u
u
u
F
u
F
u
F
u
P
B
f
D
f
C
f
f
V
P
t
s
P









(21)

where


f
T
i
C
f
d
Y
E
A
v
v
v
v
v
u
F



]
~
,
~
,
~
,
~
[
)
(







,

f
T
i
D
f
d
p
E
p
A
j
h
v
S
v
v
S
I
u
F







]
,
~
~
~
~
,
,
0
[
)
(



,
,
0
[
)
,
(

u
u
F
B
f

f
T
i
d
A
b
b
B

]
,
,



and
T
i
P
w
s
]
,
~
,
,
0
[
)
,
(






0
u
u


ar
e the convective, diffusive, su
b
grid fluxes
and source terms, respe
c
tively. For the convective fluxes,
)
(
u
F
C
f
, a monotonicity preserving reco
n
struction of the
form
)]
(
)
(
))[
(
1
(
)
(
)
(
,
,
,
u
F
u
F
u
u
F
u
F
L
C
f
H
C
f
H
C
f
C
f






is used, [43
], in which
)
(
,
u
F
H
C
f

den
otes a 2
nd

order linear reconstruction,
)
(
,
u
F
L
C
f

a 1
st

order upwind b
i
ased reconstruction and
)
(
u




a non
-
linear flux
limiter. The non
-
linear flux limiter is used to switch between the two underlying reconstruction algorith
ms, and here
the
MC limiter, [44
], is used for the momentum, e
n
ergy and species equations whereas only the higher order scheme
is used for the continuity equ
a
tion. To mini
mize the non
-
orthogonal
ity errors in the viscous and subgrid fluxes,
)
(
u
F
D
f

and
)
,
(
u
u
F
B
f
, respectively, these are split into or
thogonal

and non
-
or
tho
gonal parts, [43
]. Central
difference approximation and gradient face interpolation are used for the orthogonal and non
-
orthogonal parts
respectively. A 2
n
d

order semi
-
implicit Crank
-
Nichol
son time
-
integration scheme is used together with a PISO
-
like
algorithm, in the spirit of Rhie & Chow, for the cell
-
centered data storage stru
c
ture, [45
]. Stability is imposed using
compact sten
cil
s and by en
forc
ing c
onservation of k
i
netic en
er
gy. The equa
tions are solved sequentially, with
iteration over the non
-
linear source terms to o
b
tain rapid conve
r
gence, with a fixed CFL number of about 0.4.


First Author, Second Author
.
INSTRUCTIONS FOR THE

PREPARATION OF PAPER
S



11

5
.
Results from LES of a planar flame in homogeneous isotropic turb
ulence


The initial computational tests of the LES
-
EPaSR model are carried out for a CH
4
-
air turbulent premixed flame
propagating in homogeneous isotropic turbulence. The results from this simul
a
tion are gauged against corresponding
simulations using LES
-
P
aSR, TFM, EDC and LES
-
FPV models. To provide an independent reference against which
all LES predictions can be compared a DNS was undertaken. The computational domain chosen is a box with a side
length of 18 mm, and is initialized with an isotropic turbule
nt velocity onto which a planar flame, separating
rea
c
tants and products, is superimposed. Periodic conditions are applied to the domain boundaries orthogonal to the
flame and w
ave transmissive conditions, [46
], are used at the inlet and ou
t
let, with an in
let velocity of
)
,
,
(
)
,
,
(
t
z
y
s
t
z
y
x
u
in
w
e
v



, in which
e
x

is the unit vector in the x d
i
rection and
w


a

synthetic turbulence
field, [47
]. To avoid the flame propagating towards the i
n
let and outlet boundaries the entire computational domai
n is
allowed to move with the mean flame such that the flame is always well within the computational box. The target
rms
-
velocity fluctu
a
tions, laminar flame speed, flame thickness, integral, Tay
lor and Ko
l
mo
gorov scales are
v

=3.12
m/s, s
u
=0.38 m/s,

u
=0.04 mm,

I
=1.81 mm,

T
=0.46 mm and

K
=0.029 mm, re
spectively, such that v´/s
u
≈8.2,

I
/

u
≈45, placing the flame in the lower part of the distributed r
e
action region.
A
l
ternatively, can the case be
characterized by Re
I
≈374, Re
T
≈44, Ka≈3.4 and Da=5.5 as
summariz
ed in Table 3
.
Here,

u

is based on the ratio
between D and s
u
, usually pla
c
ing the estimate of

u

on the thinner side. For the DNS a grid of 640
3

cells is
e
m
ployed, corr
e
sponding to a grid re
so
lution of 0.028 mm, and for the LES a grid of 64
3

ce
lls is used, corr
e
sponding
to a grid resol
u
tion of 0.281 mm. This implies that each LES cell is resolved by 10
3

cells in the DNS. In this DNS,
K


is resolved by one cell across whereas

u

is resolved by two cells across. This is not id
eal, but is reasonable
considering that no detailed rea
c
tion. All LES and DNS computations were run for three eddy turnover times,
)
/
(
3
v
T


≈0.43
ms, as implied by Han & Huh, [48
], to be sufficient for the statistics of the turbulent consumpti
on
velocity to have converged.



Table 3
. Characteristic properties for Planar flame propagating in h
omogeneous isotropic turbulence


Case

v


[m⽳]

s
u

[m/s]


u

[m]


I

[m]


T

[m]


K

[m]

v


u


I
/

u


T





1

3.12

0.38

0.04

0.00181

0.00046

0.00029

8.2



㌷3

3.4

5.5



䙩Fur攠 1 show r敳e汴l from 䑎匬 䑎D fi汴敲敤 on瑯 瑨e LE匠 gr楤, 慮d LES
-
EP慓a mod敬e pr敤楣瑩ins,
r敳e散瑩t敬y, of 瑨攠f污m攠prop慧慴楮g from 瑨攠ho琠produ捴cmix瑵r攠(r楧h琩 瑯 瑨攠捯汤 r敡捴慮琠m楸瑵r攠(汥l琩 慦瑥t
two 敤dy 瑵rn
-
ov敲 瑩m敳
楮 t敲ms of vor瑩捩ty 捯r敳 (b汵攩, 瑥
m
p敲慴ar攠捯n瑯urs 慮d iso
-
surf慣es of 瑨e CH
4

mass
fraction (gray) and heat release (ye
l
low). The DNS results in figure 1a show that the vortices in the cold reactant
mixture distribute almost ra
n
domly in space, where
as when passing through the flame these vortices d
e
cay in
strength due to the increase of viscosity and volumetric expansion due to the higher temperature of the pro
d
uct
mixture. The strongest vortices survive passing through the flame, and exit into the h
ot pro
d
uct mixture with
preferred directions either perpendicular or parallel to the flame. The vortices tra
v
ersing the flame transport cold
reactants into the hot product mixture where pockets of cold rea
c
tants are formed, [49
]. In addition, these penetra
ting
vortices also wri
n
kle the flame, creating a multiply
-
wrinkled flame. At the cusps of the flame it is observed that the
heat release is locally decreased due pr
i
marily to large strain rate tangential to the flame. Regions of high heat release
(>80% of
the laminar heat release) are found more frequently between penetrating vortices pr
o
truding into the hot
reactant mixture. The maximum heat release rate of this turbulent pr
e
mixed flame is about 120% of the laminar heat
release. Regions of low heat release

(<20% of the lam
i
nar heat release) occur in between regions of high he
at release,
as also found by [20
]. The heat release topology observed and described here supports the aforementioned
description of fine structures and surroundings underlying the class

of models put forward in this study. The filtered
DNS data is in this investigation obtained by filtering the DNS data using a conventional box fi
l
ter having the same
width as the LES grid. The filtered DNS results essentially reveal all fe
a
tures observed

in the DNS, but in a coarse
-
grained manner, with the smallest structures and strongest events b
e
ing eliminated or smoothened out. The LES
models tested (LES
-
EPaSR, PaSR, EDC, TFM and FPV) manage to emulate this b
e
havior to a variable degree of
success. Bo
th the LES
-
EPaSR (figure 1c) and PaSR models are quite successful in qualitatively ca
p
turing this,
whereas the LES
-
TFM and FPV models are not, as will be elaborated on next.


SESSION NUMBER & NAM
E


12

(a)
(b)
(c)


Figure 1.

Perspective view from a planar flame propagating in homo
geneous isotropic turbulence of vorticity (blue),
reaction pr
o
gress variable (gray), heat release (yellow) and temperature contours ran
-
ging from cold (blue) to hot
(red) from (a) DNS, (b) filtered DNS and (c) LES
-
EPaSR.



Next, quantitative comparisons be
tween the different LES model predictions with the DNS and filtered
DNS data are presented. To achieve this, isomorphisms from R
3

R to R
2

are defined so that
)
~
(
~
)
,
(
~
c
f
t
f
i

x
, in
which
c
~

is the progress variable
)
/(
)
~
(
~
u
fu
b
fu
u
fu
fu
Y
Y
Y
Y
c



, with
b
fu
Y

and
u
fu
Y

being the burnt and unburnt
fuel mass fractions, r
e
spectively and

˜
f

an LES field. In fi
g
ures 2a to 2c isomorphism maps (i.e. scatter plots) of
temperature,
T
~
, CO mass fraction,
CO
Y
~
, and vorticity magnitude,
|
~
|
ω
, respectively, are shown. Figure 2a indicate
that the LES
-
EPaSR, PaSR, EDC and TFM models all predict
T
~

in good accordance with
the DNS and filtered
DNS data, whereas the LES
-
FPV model overpredicts
T
~

in the hot products (for
6
.
0
~

c
). Figure 2b shows a
larger spreading of the
CO
Y
~

data as compared with the
T
~

data as a consequence of
CO
Y
~

being an intermed
i
ate
specie, whereas the almost linear dependence of
T
~

on
c
~

is a consequence of the overall thermodynamic effects. No
direct comparison for
CO
Y
~

can be made for the LES
-
FPV model since this model is formulated explicitly in terms of
c
~
, [39
]. The LES
-
EPaSR and PaSR models shows best agreement with the DNS and filtered DNS data, whereas
both the LES
-
TFM and EDC m
odels predict slightly to high values of
CO
Y
~
. Figure 2c reveals that the DNS data
show higher values of
|
~
|
ω

on the reactant side (
c
~
<0.5) than any other data sets as a cons
e
quence of all vortex
str
uctures b
e
ing included in the DNS data. On the product side (
c
~
>0.5) this difference is less evident due to the
aforementioned decay of vorticity across the flame. The LES
-
EPaSR, PaSR, EDC and TFM models all show good
agreement with the

filtered DNS data whereas the LES
-
FPV model shows lower values of
|
~
|
ω

across the entire
range of
c
~
. So far we conclude that the LES
-
EPaSR, PaSR models show good agreement with the DNS data,
closely followed by the
LES
-
EDC and TFM models, whereas the LES
-
FPV model show less satisfactory agreement.



(a)
(b)
(c)


Figure 2.

Isomorphism maps (s
catter
-
plots) of (a) temperature, (b) CO mass fraction and (c) vorti
c
ity magnitude.
Legend: (

) DNS, (

) fi
l
tered DNS, (

) LES
-
PaSR, (

) LES
-
EDC, (

) LES
-
TFM, (

) LES
-
EPaSR and (

) LES
-
FPV.


First Author, Second Author
.
INSTRUCTIONS FOR THE

PREPARATION OF PAPER
S



13


Next we attempt to compare some key quantities of the different LES combustion models with their
definitions as obtained from the DNS data. The quantities chosen are the flame surface density
,

, the filtered heat
release,
Q

, and the reacting volume fraction,
*

. Based on the DNS data the flame surface density is evaluated as
)
(
|
|
*
c
c
c





, with


being the Dirac fun
c
tion,
the overbar denoting box
-
filtering from the 640
3

DNS grid
to the 64
3

LES grid, and
8
.
0
*

c

to match the peak reaction rate of CH
4
-
air combustion. Similarly, the filtered heat
r
e
lease is computed as
)
)
,
,
(
(
,
1


i
f
k
j
ij
i
N
i
h
Y
T
w
P
M
Q





, whereas the reactin
g volume fraction is computed as
3
*
/
0



Q
Q
V



, in which
0
Q
Q
V




is the volume in which
Q


exceeds the threshold value
0
Q

, selected to be 5% of
the peak laminar heat release. The flame surface

density is used in this comparison is obtained from its modeling,




/
)
1
(
4
c
c

, in the LES
-
FPV model, and from the estimate
)
~
(
|
~
|
*
c
c
c






in the remaining LES
models. Similarly, the heat
-
release is obtained from



h
s
Q
u
u



i
n the LE
S
-
FPV model, where ∆h is the heat of
combustion, and from
)
)
,
,
(
(
,
*
*
*
1



i
f
k
j
ij
i
N
i
h
Y
T
w
P
M
Q






in the LES
-
EPaSR, PaSR and EDC models, and
from
)
)
~
,
~
,
(
(
,
1


i
f
k
j
F
E
ij
i
N
i
h
Y
T
w
P
M
Q






in the LES
-
TFM model. The reaction volume fraction is pr
o
vided by
directly by the LES
-
EPaSR, PaSR and ED
C models, and as estimated as E/F for the LES
-
TFM model. Figure 3a
reveals that all LES predi
c
tions of


result in
parabolic
-
like


di
s
tributions, with


from the LES
-
FPV model being
sy
m
metric around
5
.
0
~

c

(by co
n
stru
ction), whereas the


distributions from LES
-
EPaSR, EDC, PaSR and TFM
are all slightly skewed t
o
wards higher values of

˜
c
. The filtered DNS data also appears to be slightly skewed towards
higher values of

˜
c
, as a
re
recent experimental data, [50
-
52
], thus suggesting that the LES
-
EPaSR, EDC, PaSR and
TFM models provide a more accurate description of the flame kinema
t
ics than does the LES
-
FPV model. Figure 3b
shows that all LES predictions, except that of the LES
-
FPV mo
del
result in
Q


distributions skewed towards
8
.
0
*

c

as expected from CH
4
-
air combustion, whereas the results from the LES
-
FPV model show a symmetric
Q


distr
i
bution centered around
5
.
0
~

c

(by construction). The filtered DNS data also show a
Q


distribution that
is skewed t
o
wards
8
.
0
*

c
, having similar fluctuation levels as the
LES
-
EPaSR, PaSR and EDC mo
d
els. The heat
release of the LES
-
TFM model

is here somewhat underpredicted, with its peak values only reaching 80% of the peak
laminar reaction rate, whereas the peak reaction rate of the LES
-
EPaSR, PaSR and EDC models are closer to 120%
as indicated by the filtered DNS data. In fi
g
ure 3c we compa
re the different estimates of
*

. The filtered DNS and
the
LES
-
EPaSR, PaSR and EDC models all result in a rather narrow band of data around
8
.
0
*


, whereas the
LES
-
TFM model result in an almost constant value of
4
.
0
*


. The LES
-
EPaSR model pr
o
vides best agre
e
ment
with the filtered DNS data, closely followed by the PaSR model, whereas the EDC model results in somewhat higher
values of
*

. Consequently, these compar
i
sons suggest that the
LES
-
EPaSR and PaSR models show the best
agreement with the filtered DNS data.



(a)
(b)
(c)


Figure 3.

Isomorphism maps (s
catter
-
plots) of (a) flame surface density, (b) heat release and (c) reac
t
ing fra
c
tion.
Legend: (

) DNS, (

) filtered DNS, (

) LES
-
P
aSR, (

) LES
-
EDC, (

) LES
-
TFM, (

) LES
-
EPaSR and (

) LES
-
FPV.



SESSION NUMBER & NAM
E


14

5
.
Summary and concluding remarks


This paper proposes a novel LES finite rate chemistry combustion model for turbulent flows. We start developing
this model from a simple cartoon of turbulent

flow (subsequently also involving combustion) that has evolved ever
since Kolmogorov proposed his 1941 th
e
ory, [22]. Batchelor & Townsend, [18
], made a strong case for non
-
uniform
spa
tial distrib
u
tion of the turbulent fine
-
structures based on experimenta
l evidence in high Re number flows and thus
also of molecular mixing and diss
i
pation. Kuo & Corrsin, [24
], showed that such fine
-
structure regions consist of
vor
tex sheets, rib
bons and tubes folded in certain r
egions of the flow. Chomiak, [16
-
17
], r
e
fine
d these theories and
extended them into the field of combustion, proposing that
at high Re nu
m
bers turbulent fine structure is not
uniformly distributed but concentrated into small isolated r
e
gions, whose entire volume is a small fraction of the
whole volu
me.
DNS, e.g. [19
], have pro
vided su
p
port to this view, revealing that high in
tensity vortices amal
ga
mate
in fila
ments em
bedded into sheets of low intensity vor
ticity. For combustion, a similar orga
nization of the flow has
been o
b
served from DNS, e
.g. [20
] and the results in this paper, supporting t
he earlier ideas of Chomiak, [16
-
17
].
Such DNS suggests that the fine
-
structure vortices at the flame are essentially para
l
lel to the flame whereas those
behind the flame are mostly perpendicular to the f
lame. Re
gions of heat release and volumetric expansion are found
to be regularly distributed among the fine
-
structure vor
tices in a background of low heat release and volumetric
expansion. Ma
g
nussen, [13
], used a similar cartoon of turbulent combustion t
o develop the Eddy Dissipation
Concept (EDC)


the first he
t
erogeneous multi
-
scale turbulent combustion model ever used.


The basic idea of the proposed model is to use the math
e
matical treatment of multi
-
phase flows for the
description of fine
-
scale struc
tures dissolved in the background turbulence, chara
c
terized by low intensity mixing.
The model results in closing the filtered reaction rates (typically described by Arrhenius rate expressions),
)
,
,
(
k
i
Y
T
w


, using the same rate expressions but
evaluated at the fine
-
structure conditions and taking the fine
structure volume fraction into a
c
count so that
)
,
,
(
)
,
,
(
*
*
*
k
i
k
i
Y
T
w
Y
T
w






. In the proposed LES
-
EPaSR
model, the fine stru
c
ture conditions
}
,
{
*
*

T
Y
i

are computed by sol
v
ing a set of t
ransport equations, thus being more
general than the LES
-
PaSR model in which the fine structure conditions were obtained from l
o
cal algebraic
equations. For both the LES
-
PaSR and EPaRR models, the reacting fine structure volume fraction,
*

, and the
subgrid time scale,
*

, needs to be provided by submodels. Based on the multi
-
phase analogy, a separate transport
equation can be form
u
lated for the reacting vo
l
ume fraction such that
*
*
*
*
*
/
)
(
)
~
(
)
(









eq
t






v
, in which
*
eq


and
*

, needs to be modeled. A model for the
subgrid time scale is here also proposed, based essentially on geome
t
rical considerations, resulting in that





K
*
, in which
K


is

the Kolmogorov time scale and



is the time scale of the resolved shear.
This
particular time scale was recently found by Y
e
ung
et al
, [53
], to be the time scale that best characterized the
Lagrangian statistics of dissip
a
tion, and t
herefore appears as a reasonable choice for representing the subgrid
processes of inte
r
est to mixing and combustion at high Re
-
numbers.
A model for the equili
b
rium reacting volume
fra
c
tion,
*
eq

, is proposed based on dimensional argument
s and analysis of two cases (fast and slow chemistry in
co
m
parison with mi
x
ing) that takes the form
)
/(
*
*






c
c
eq
, in which

c

is a chemical time scale. When
subgrid convection is negligible the transport equation for


*

d
e
ge
nerates into the analyt
i
cal expression
)
/(
*
*






c
c

that is e
m
ployed in the LES
-
PaSR model, [12].


The proposed model, hereafter referred to as the LES
-
EPaSR model, was then successfully evaluated
against other well
-
known LES combustion models
such as the LES
-
PaSR, EDC, TFM and FPV models for a planar
flame propagating in h
omogeneous isotropic turbulence.

For the planar flame in homogeneous is
o
tropic turbulence
DNS was also carried out to provide a detailed reference against which all LES models

could be bechmarked.
Comparisons were carried out both against the raw DNS data and against the filtered DNS data. Best agreement was
obtained by the proposed LES
-
EPaSR model closely followed by the PaSR model.


Acknowledgement
s


Financial support for th
is work was provided by the Office National d'Etudes et de Recherches Aero
spatiales (VS)
and by the Swedish Armed Forces (CF). Helpful discussion with A. Mura and his co
m
ments are gratefully
acknowledged.


First Author, Second Author
.
INSTRUCTIONS FOR THE

PREPARATION OF PAPER
S



15

References


[1]
T. Poinsot and D. Veynante, Theo
retical and Numerical Combustion, R. T. Edwards, ed., Phil
a
delphia, 2001.

[2]
P. Sagaut, Large Eddy Simulation for Incom
pressible Flows, Springer Verlag, Heide
l
berg, 2001.

[3]
F.F. Grinstein, L. Margolin and B. Rider (Eds.), Implicit Large Eddy Simulati
on: Compu
t
ing Turbu
lent Fluid
Dynamics, Cambridge Un
i
versity Press, New York, 2007.

[4]
J. Janicka and A. Sadiki, Large Eddy Simulation of Turbulent Combustion Systems, Proc.
Comb. Inst. 30 (2005)
pp. 537
-
565.

[5]
H. Pitsch, Large Eddy Simulation of Tur
bulent Combustion, Annu.
Rev Fluid Mech., 38 (2006) pp. 453
-
483.

[6]
C. Fureby, LES Modeling of Combustion for Propulsion Applications, Phil.
Trans. R. Soc. A, 367 (2008) pp.
2957
-
2969.

[7]
P. Wang and X.S. Bai, Large Eddy Simulation of Turbulent Premixe
d Flames using Level
-
set G
-
equ
a
tion, Proc.
Comb. Inst. 30 (2005) pp. 583
-
591.

[8]
E.R. Hawkes and R.S. Cant, A Flame Surface Density Approach to Large Eddy Simul
a
tion of Pr
e
mixed
Turbulent Combustion, Proc.
Comb.. Inst. 28 (2000) pp. 51
-
59.

[9]
R. Knikk
er and D. Veynante, Experimental Study of the Filtered Progress Variable A
p
proach for LES of
Premixed Combustion, in Advances in LES of Complex Flows, (Eds.) R. Fri
e
drich and W. Rodi, Kluwer, Dordrecht,
2000 pp. 353
-
366.

[10]
P Givi, Filtered Density Func
tion for Subgrid Scale Modeling of Turbulent Combu
s
tion, AIAA.J., 44 (2006)
pp. 16
-
23.

[11]
O. Colin, F. Ducros, D. Veynante and T. Poinsot, A Thickened Flame Model for Large Eddy Simul
a
tions of
Turbulent Premixed Combustion, Phys.
Fluids 12 (2000) pp.
1843
-
1863.

[12]
M. Berglund, E. Fedina, J. Tegnér, C. Fureby and V. Sabelnikov, Finite Rate Chemistry LES of Self I
g
nition in
a Supersonic Combustion Ramjet, AIAA.J. 48 (2010) pp. 540
-
550.

[13]
B.F. Magnussen, On the Structure of Turbulence and a General
ised Eddy Dissipation Co
n
cept for Chemical
Reactions in Turbulent Flow, 19
th
AIAA Sc. meeting, St. Louis, 1981, USA

[14]
R.W. Bilger, Conditional Moment Closure for Turbulent Reacting Flow, Phys.
Fluids A 5 (1993), pp. 436
-
444.

[15]
V. Sankaran and S.

Menon, Subgrid Combustion Modeling of 3D Premixed Flames in the Thin
-
Reac
tion
-
Zone
Regime, Proc.
Comb Inst. 30 (2005), p 575
-
583.

[16]

J. Chomiak, A Possible Propagation Mechanism of Turbulent Flames at High Re
y
nolds Numbers, Comb.
Flame, 15 (1970) pp.
319
-
321.

[17]
J. Chomiak, Basic Considerations in the Turbulent Flame Propagation in Premixed Gases, Progress Energy
Comb. Sci., 5 (1979) pp. 207
-
221.

[18]
G.K. Batchelor and A.A. Townsend, The Nature of Turbulent Motion at Large Wave
-
numbers, Proc. Roy.

Soc.
London A, 199 (1949) pp. 238
-
255.

[19]
P.R. Woodward, D.H. Porter, I. Sytine, S.E. Anderson, A.A. Mirin, B.C. Curtis, R.H. Cohen, W.P. Dannevik,
A.M. Dimits, D.E. Eliason, K.
-
H. Winkler and S.W. Hodson, Very High Resolution Simulations of Compressib
le
Turbulent Flows, Computational Fluid Dynamics, Pr
o
ceedings of the 4th UNAM Supercomputing Conference
Mexico City, June 2000, R
a
mos E., Cisneros G., Fernandez
-
Flores A., & Santillan
-
Gonzalez A., eds., World
Scie
n
tific, p. 3..

[20]
M. Tanahashi, M. Fujim
ura and T. Miyauchi, Coherent Fine Scale Eddies in Turb
u
lent Premixed Flames, Proc.
Comb.. Inst. 28 (2000), pp. 579
-
587.

[21]
I.S. Ertesvåg and B.F. Magnussen, The Eddy Dissipation Turbulence Energy Ca
s
cade Model, Comb. Sci.
Tech.
159 (2000), pp. 213
-
236.

[22]
A.N. Kolmogorov, Local Structure of Turbulence in an Incompressible Fluid for Very Large Re
y
nolds
Numbers, Comptes rendus (Doklady) de l’Academie des Sciences de l’U.R.S.S., 31 (1941) pp. 301
-
305. Reprinted in
Friedlander S.K. & Topper L., eds., 196
1, “Turbulence: Classic Papers on Statistical Theory”, NewYork:
Inte
r
science Publishers

[23]
H. Tennekes and J.L. Lumley, A First Course of Turbulence, MIT Press, Cambridge, MA, USA, 1972.

[24]
Y.S. Kuo and S. Corrsin, Experiments on Internal Intermitten
cy and Fine Structures Di
s
tribution Func
tions in
Fully Turbulent Fluid, J. Fluid Mech., 50 (1971) pp. 285
-
319.

[25]

A.N. Kolmogorov, A Refinement of Previous Hypotheses Concerning the Local Structure of Turb
u
lence in a
Viscous Incompressible Fluid at Hig
h Reynolds Number, J. Fluid Mech., 13 (1962), pp. 82
-
85.

[26]
A.M. Obukhov, Some Specific Features of Atmospheric Turbulence, J. Fluid Mech. 13 (1962) pp. 77
-
81.

[27]
L.D. Landau and E.M. Lifshitz, Fluid Mechanics, (1959) Footnote on p. 126, Perg
a
mon, Lo
ndon

[28]

A.S. Monin and A.M Yaglom, Statistical Fluid Mechanics: Mechanics of Tu
r
bulence, 1971, Vol. 2, MIT Press,
Ca
m
bridge, MA, USA.

[29]
U. Frisch. Turbulence: The Legacy of A.N. Kolmogorov, Cambridge and NewYork: Ca
m
bridge Unive
r
sity
Press, 1995.

SESSION NUMBER & NAM
E


16

[3
0]
A. Tsinober, An Informal Conceptual Introduction to Turbulence, Springer, 2009.

[31]
K.R. Sreenivasan and R.A. Antonia, The Phenomenology of Small
-
scale Turb
u
lence, Ann.
Rev. Fluid Mech.,
29 (1997) pp. 435
-
472.

[32]
K.R. Sreenivasan and C. Meneveau,

The Fractal Facets of Turbulence, J. Fluid Mech. 173 (1986) pp. 356
-
386.

[33]
R. Bensow and C. Fureby, On the Justification and Extension of Mixed Models in LES, J. Turb. 8(54) (2007)
pp 1
-
12.

[34]
C. Fureby, On LES and DES of Wall Bounded Flows, Ercoft
ac Bulletin No 72, 2007, Marsh I
s
sue.

[35]
F. Tichy, Laser
-
Based Measurements in Non
-
premixed Jet Flames, Ph.D. Th
e
sis, Dept. of Appl.
Mech.,
Thermodynamics and Fluid Dynamics, Norwegian University of Sc
i
ence and Technology, Trondheim, No
r
way
(1997).

[36]

B.F. Magnussen, The Eddy Dissipation Concept, ECCOMAS Thematic Confe
r
ence on Comput
a
tional
Combustion, Lisbon, Portugal, 21
-
24 June (2005).

[37]
D.A. Drew, Mathematical Modeling of Two
-
Phase Flow, Ann.
Rev. Fluid Mech., 15 (1983), pp. 261
-
291.

[38]
E.
Fedina and C. Fureby, A Comparative Study of Flamelet and Finite Rate Chemistry LES for an
Ax
i
symmetric Dump Combustor, Submitted to J. Turb.
(2010).

[39]
M. Boger, D. Veynante, H. Boughanem and A. Trouvé, Direct Numerical Sim
u
lation Analysis of Flame
Sur
face Density Concept for Large Eddy Simulation of Turbulent Pr
e
mixed Combustion, Proc.
Comb. Inst. 27
(1998), pp. 917
-
925.

[40]
G.P. Smith, D.M. Golden, M. Frenklach, N.W. Moriarty, B. Eiteneer, M. Goldenberg, C.T. Bowman, R.K.
Ha
n
son, S. Song W.C. Gardin
er, V.V. Lissianski and Z. Qin, http://www.me.ber
-
keley.edu/gri_mech/

[41]
E. Giacomazzi, V. Battaglia and C. Bruno C, The Coupling of Turbulence and Chemi
s
try in a Premixed Bluff
-
body Flame as Studied by LES, Comb.
Flame, 138 (2004), pp 320
-
335.

[42]
H.
G. Weller, G. Tabor, H. Jasak and C. Fureby, A Tensorial Approach to CFD using O
b
ject Oriented
Techniques, Comp. in Physics, 12 (1997), pp. 620
-
631.

[43]
D. Drikakis, C. Fureby, F.F. Grinstein & M. Liefendahl, ILES with Limiting Alg
o
rithms, In Implicit La
rge
Eddy Simulation: Computing Turbulent Fluid Dynamics, Eds.
Gri
n
stein F.F., Margolin L. & Rider B., Cambridge
University Press, 2007, pp 94
-
129.

[44]
B. van Leer, Towards the Ultimate Conservative Difference Scheme III.
U
p
stream
-
Center
-
ed Finite
-
Differe
nce
Schemes for Ideal Compressible Flow, J. Comp. Phys., 23 (1977), pp 263
-
275.

[45]

C.M Rhie and W.L. Chow, Numerical Study of the Turbulent Flow Past an Airfoil with Trailing Edge
Separation, AIAA J., 21 (1983), pp. 1525
-
1532.

[46]
T.J. Poinsot and S.K
. Lele, Boundary Conditions for Direct Simulation of Co
m
pressible Viscous Reac
t
ing
Flows, J. Comp.
Phys., 101 (1992), p 104
-
129.

[47]
N. Jarrin, S. Benhamadouche, D. Laurence and R. Prosser, A Synthetic
-
Eddy
-
Method for Generating Inflow
Conditions for Lar
ge Eddy Simulations, Int. J. Heat and Fluid Flow, 27 (2006), pp. 585
-
593.

[48]
I. Han and K.Y. Huh, Roles of Displacement Speed on Evolution of Flame Surface De
n
sity for Different
Turbulent Intensities and Lewis Numbers in Turbulent Premixed Combu
s
tion, C
omb.
Flame, 152 (2008), pp 194.

[49]
J.H. Chen, T. Eckekki and W. Kollmann, The Mechanism of Two
-
Dimensional Pocket Formation in Lean
Premixed Methane
-
Air Flames with Implications to Turbulent Combustion, Comb.
Flame, 116 (1998), pp. 15
-
48.


[50]
J. Hul
t, S. Gashi, N. Chakraborty, M. Klein, K.W. Jenkins, S. Cant and C.F. Kaminski, Measurement of Flame
Surface Density for Turbulent Premixed Flames using PLIF and DNS, Proc. Inst. Comb, 30 (2007), pp. 1319
-
1326.

[51]
Z.S. Li, B. Li, Z.W. Sun, X.
-
S. Bai and

M. Aldén, Turbulence and Combustion Intera
c
tion: High Resolution
Local Flame Front Structure Visualization using Simultaneous Single
-
Shot PLIF Imaging of CH, OH, and CH
2
O in a
Piloted Premixed Jet Flame, Comb.
Flame, 157 (2010), pp. 1087
-
1096.

[52]
D. Ve
ynante, G. Lodato, P. Domingo, L Vervisch and E. Hawkes, Estimation of Three
-
Dimensional Flame
Surface Densities from Planar Images in Turbulent Premixed Co
m
bustion, Exp. Fluids, 49 (2010), pp. 267
-
278.

[53
]
P.K. Yeung, S.B. Pope and B.L. Sawford, Reynold
s Number Dependence of Lagra
n
gian Statistics in Large
Numerical Simulations of Isotropic Turbulence, J. Turb., 7 (2006), pp. 1
-
12.