Weller B Xi Combustion Model PDF

The Development of a New Flame Area
Combustion Model
Using Conditional Averaging
H.G. Weller
Thermo-Fluids Section Report TF/9307
Department of Mechanical Engineering
Imperial College of Science Technology and Medicine
23th March 1993
Contents
1 Introduction 6
1.0.1 Analysis of EBU based Combustion Models : : : : : : : : 7
2 Development of Flame Area Models 11
2.1 Conditional Averaging : : : : : : : : : : : : : : : : : : : : : : : : 11
2.1.1 Formulae for Conditional Averaging : : : : : : : : : : : : : 12
2.1.1.1 Conditional Averaging Dierential Operations : : 12
2.1.1.2 Decomposition Operations : : : : : : : : : : : : : 13
2.1.2 Conditional Averaging Transport Equations : : : : : : : : 13
2.1.3 Interface Transport Equation : : : : : : : : : : : : : : : : 14
2.1.4 A General Denition for Wrinkle Factor : : : : : : : : : 15
2.1.5 Transport Equation for b in terms of : : : : : : : : : : : 16
2.1.6 Transport Equation : : : : : : : : : : : : : : : : : : : : 17
2.1.6.1 jrbj Transport Equation : : : : : : : : : : : : : : 17
2.1.6.2 Transport Equation : : : : : : : : : : : : : : : 19
2.1.6.3 Transport Equation Assembly : : : : : : : : : 19
2.1.7 Modelling the b Transport Equation : : : : : : : : : : : : : 22
2.1.8 Modelling the Transport Equation : : : : : : : : : : : : 25
2.2 Conditionally Averaged Velocities : : : : : : : : : : : : : : : : : : 28
2.2.1 Conditionally Averaged Continuity Equations : : : : : : : 29
2.2.2 Conditionally Averaged Momentum equations : : : : : : : 29
2.2.3 Solution of the Conditionally Averaged Equations : : : : : 31
2.3 Ensemble Averaged Equations : : : : : : : : : : : : : : : : : : : : 32
2.3.1 Transport Equation for eb : : : : : : : : : : : : : : : : : : : 32
2.3.2 Transport Equation for : : : : : : : : : : : : : : : : : : 35
2.3.3 A One Equation Model : : : : : : : : : : : : : : : : : : : : 35
3 Assessment of the eb { Models 37
3.1 Constant Density : : : : : : : : : : : : : : : : : : : : : : : : : : : 39
3.2 Variable Density : : : : : : : : : : : : : : : : : : : : : : : : : : : 42
3.3 Boundary Layer Flame Propagation : : : : : : : : : : : : : : : : : 47
4 Conclusion 52
Contents 2
5 Acknowledgement 54
List of Figures
2.1 Flame Interface : : : : : : : : : : : : : : : : : : : : : : : : : : : : 12
2.2 Ensemble of Flame Interfaces : : : : : : : : : : : : : : : : : : : : 16
2.3 Construction of Term III for a Gaussian Flame : : : : : : : : : : 21
2.4 \Cusp" Formation for a Laminar Flame in a Taylor Vortex : : : : 22
2.5 Modelling Surface Averages : : : : : : : : : : : : : : : : : : : : : 25
3.1 One Dimensional Open Duct Test Case : : : : : : : : : : : : : : : 37
3.2 One Dimensional Semi-Closed Duct Test Case : : : : : : : : : : : 38
3.3 eb { Model b Prole : : : : : : : : : : : : : : : : : : : : : : : : : 39
3.4 eb { Model Prole : : : : : : : : : : : : : : : : : : : : : : : : : 40
3.5 eb { Model Prole : : : : : : : : : : : : : : : : : : : : : : : : : 40
3.6 \Cusp" Formation and Dissipation : : : : : : : : : : : : : : : : : 41
3.7 Eect of Doubling \Cusp" Dissipation : : : : : : : : : : : : : : : 42
3.8 \eddy break-up" Model b Prole : : : : : : : : : : : : : : : : : : : 43
3.9 eb { Model b Prole : : : : : : : : : : : : : : : : : : : : : : : : : 43
3.10 \eddy break-up" Model Flame Propagation : : : : : : : : : : : : : 44
3.11 eb { Model Flame Propagation : : : : : : : : : : : : : : : : : : : 45
3.12 eb { Model Development : : : : : : : : : : : : : : : : : : : : : 45
3.13 1 Equation Model Flame Propagation : : : : : : : : : : : : : : : : 46
3.14 eb { and \eddy break-up" Model Exit Velocity : : : : : : : : : : 47
3.15 Boundary Layer Test Case : : : : : : : : : : : : : : : : : : : : : 48
3.16 eb { Model Flame Propagation : : : : : : : : : : : : : : : : : : : 49
3.17 eb { 1-Equation Model Flame Propagation : : : : : : : : : : : : 50
3.18 \Eddy break-up" Model Flame Propagation : : : : : : : : : : : : 51
Nomenclature
Roman Symbols
A Area
A Combustion model coecient
b Regress variable (normalised fuel concentration) (b = 1 ? c)
C Constant
c Progress variable (normalised product concentration) (c = 1 ? b)
D Diusive transport term
D Diusion coecient
f fuel mass fraction
G Generation coecient
G Generation rate
h Specic enthalph
I Indicator function
I Identity matrix
k Turbulence kinetic energy
l Turbulence integral length scale
l Spectral length scale
n Unit normal
p Pressure
Q Tensor property of any rank
R Gas constant
R Source term
S Surface
S Laminar ame speed
T Temperature
t Time
U Velocity
V Volume
x Spatial position
Nomenclature 5
Greek Symbols
Phase volume fraction
Total ame \wrinkling"
(l) Flame \wrinkle" density function
Length scale
Flame area per unit volume
Stress tensor
" Turbulence dissipation rate
Laminar kinematic viscosity
Equivalence ratio
Density
Time scale
Subscripts
Qb Burnt gas
Qeq Equilibrium
QI Interface
Qs Surface
Qt Turbulent
Qu Unburnt gas
Q Zone in which ! 1
Q Kolmogorov
Superscripts
QT Transpose
Q000 Fluctuation with respect to ensemble average
Q Fluctuation with respect to density weighted ensemble average
Q] Fluctuation with respect to surface average
Overlines
Q Ensemble average
Q Conditional ensemble average
Q
f Density weighted ensemble average
Q
f Conditional density weighted ensemble average
Q Surface average
z}|{
Q_ Rate of change
Chapter 1
Introduction
The need for turbulent combustion modelling became apparent through the study
of the computational cost of handling ame kinetics reactions schemes directly.
These problems are two-fold, rstly there is the the complexity of the reaction
schemes, although this problem may be overcome by the use of simplied schemes
[1], and secondly there is the problem of averaging, the average reaction rates
cannot simply be expressed in terms of the average temperature [2]. With the
advent of \transported joint PDF" techniques for handling the coupling between
the kinetics and the ow [3], the need for turbulent combustion models appears
to have lessened. It must, however, be remembered that \transported PDF"
methods for turbulent combustion closure are computationally expensive, the
extent of which is only now becoming apparent. Also there are numerical and
modelling diculties, particularly with respect to the turbulence and mixing,
with these schemes which require further research. Meanwhile there is still much
that can be gained from the more traditional approaches to turbulent combustion
modelling, and under some condition, for example in internal combustion engines,
these methods may prove superior.
There are several traditional approaches to combustion modelling, although
there are marked similarities between many of them, even between some derived
using fundamentally dierent methods. Under some combustion conditions the
chemical kinetics may be assumed innitely fast in which case the combustion
process is controlled by mixing. Using this idea, Spalding [4] proposed the \eddy
break-up" (EBU) model which has spawned many derivatives [5, 6, 7, 8]. One
form of these were developed in an attempt to extend the range of applicability of
the basic model by including chemical kinetics eects, often by the introduction
of the laminar ame speed to represent the combustion process occurring between
the unburnt and burnt gases as they mix [9, 10, 11, 12, 13].
Alternative methods for combining both turbulent mixing and chemical ki-
netics eects in turbulent combustion models include the \presumed PDF" [14,
15, 16, 17] and the \laminar- amelet" [18, 19, 20, 21, 22, 23, 24, 25, 26, 20]
approaches. In the former, the reaction rate probability density function is no
Introduction 7
longer assumed to be a double function, implicit in the EBU models, but some
other simple alternative is chosen, and kinetics eects are then included directly
using some simplied kinetics scheme. The \laminar amelet" model is funda-
mentally similar, in that a PDF is assumed, usually a double function with or
without embellishments, the dierence being in the way the kinetics information
is introduced. It is supposed that the turbulent ame may be represented as a
distribution of reaction zones, \ amelets", each of which is considered to be a
laminar ame evaluated at the local conditions. In simple forms of this modelling
the laminar amelets are assumed to be undisturbed by the turbulent ow eld
but in more recent work the local stretch and curvature eects are incorporated.
Many of these models bare a marked similarity to the extended EBU models in
the form of the reaction expression used [23, 24], or, in the limit of fast kinetics
reduce to the EBU expression [14, 15, 16, 17], and are expected to exhibit the
same limitation as the original models under such conditions. The behaviour of
the EBU type models will now be analysed in order to clarify its limitations and
to give a reference point for comparison of the present work.
1.0.1 Analysis of EBU based Combustion Models

It is now clear from numerous tests that combustion models based on the \eddy
break-up" expression fail to reproduce the qualitative behaviour of ames near
walls, in complex ow elds and in large geometries. This is apart from the
inability of such models to reproduce the quantitative behaviour of ames without
\tuning" numerous coecients.
The deciencies of such models near walls are often attributed to the lack of
chemical kinetics input. It is thought that such a dependency would compensate
the ame acceleration due to the turbulence eect by the reduction in the chemical
reaction rate in the thermal boundary layer near the wall. It has been shown [13]
that it is possible to force the kinetics dependency to fulll the promise of such
compensation, but only if extreme modelling is applied. It is not clear from such
eort that the lack of kinetics input is the major contributing factor to the poor
model performance near walls, although it cannot be neglected.
In the EBU models, the burning zone consists of alternating fragments of
unburnt and fully burnt gas and the rate of fuel consumption is proportional to
the turbulent mixing of the unburnt and burnt gases. Thus,
we_ / cg
00 2
(1.1)
t
where cg 00 2
is the density weighted ensemble average progress variable uctuation
and t is a characteristic turbulent time scale. A consequence of the rst assump-
tion is that the probability density function of the progress variable is a double
Introduction 8
function, in which case cg00 2

= ce(1 ? ce) hence
we_ / c(1? c)
e e
(1.2)
t
In most simple combustion models, t is considered to be the \eddy turn-over"
time, and evaluated from the turbulence kinetic energy ke and dissipation rate "e
as
t = k"e
e
(1.3)
Consider the solution of the simple, modelled, density weighted ensemble aver-
aged convective-diusive transport equation for a reaction progress variable ce,
incorporating the simple \eddy break-up" combustion model :
@ce + r (U fce) ? r (D e c rce) = "e ce(1 ? ce) (1.4)
@t Ae
t
k
The form of the reaction rate expression in the above equation may be rein-
terpreted through the use of the asymptotic analysis proposed by Kolmogorov,
Petrovsky and Piscounov (KPP) [27], and applied to simple combustion models
by Hakberg [28], which suggests that the turbulent ame speed is proportional
to u0. Thus the reaction rate expression may be rewritten by decomposing the
time scale thus
"e = u0 ; (1.5)
ke l
where the turbulence integral length scale l is dened as
0
l u"k : (1.6)
Thus the ame length scale is eectively being imposed rather than it being
obtained through the solution of the transport equation. This lies at the heart
of the instability problem associated with such models, as will be demonstrated
in the following simple one dimensional analysis.
Consider the solution of the above equation for a one dimensional steady state
ame propagation through static homogeneous turbulence. Assuming the ame
reaches a stable wave-like propagation, the turbulent ame speed for this problem
can be obtained by integrating the volumetric reaction rate across the ame from
the unburnt gas to the fully burnt gas [8]:
Z1
St = R_ dl (1.7)
?1 ec
Introduction 9
where
0
R_ ec = A ul ce(1 ? ce) (1.8)
For such a simple case the above integral may be decomposed into a characteristic
length scale for the ame f and the integral of R_ ec in ce space
Z ec=1
St = f R_ ec dc (1.9)
ec=0
Clearly if the length scale of the ame increases the speed will increase pro-
portionally. So, if the initial hypothesis of the ame scale being proportional to
the turbulent length scale is correct, then the ame speed will be stable. How-
ever, if there is an instability in the equation set such that the ame thickness
increases, then the speed will also increase.
To test this point numerical solutions of Eqn. (1.4) have been obtained for
the aforementioned 1D circumstances, using a ne mesh and small time step to
eliminate numerical errors. Initially, the heat release is set to zero so that no
combustion-induced ow is generated. The result, is, as expected, a stable wave-
like ame which propagates at a constant speed proportional to the prescribed
u0 and has a width of approximately one integral length scale. This result is
consistent with KPP analysis, and supports the premise that the ame thickness
is known and may be incorporated directly into the model. However, with heat
release and the associated ow changes due to gas expansion, the behaviour is
quite dierent. Now, the ame has a tendency to expand due to the divergent
velocity eld, which increases the ame speed. The increase in ame thickness
increases the integral heat release, which in turn maintains the velocity diver-
gence. Hence the coupling between the equation for ce and the density weighted
ensemble average momentum equations is unstable and unphysical. This result
is inconsistent with KPP analysis due to it's reliance on the assumption that the
ame propagation is governed by the behaviour of the ame tip rather than the
integral heat release. An alternative method for analysing the ame propaga-
tion, which make no assumption about ame behaviour is eigenvalue analysis [8]
which, for this case, does not produce a single ame speed eigenvalue indicating
that the model is unstable.
It is now time to review the fundamental basis of the simple models by com-
parison with recent multi-equation models derived from a more physical under-
standing of ame behaviour. Through such analysis it may be possible to assess
the deciencies of the simple models and propose more qualitatively reasonable
ones by careful reduction of the more complex schemes.
An extension of the \laminar amelet" approach proposed by Marble and
Broadwell [29], and further developed by many others [30, 31, 32], considers the
ame as a sheet propagating locally as a laminar ame and increasing in area
due to the formation of distortions by its interaction with the turbulent ow
Introduction 10
[33]. These models comprise transport equations for the progress variable and
the ame area per unit volume , and do not resemble the EBU models in any
way. The physical approach used to derive such models makes them appealing
and enables development in many directions.
The present work is an alternative approach to the traditional based models
in that the ame distribution is represented by the ame \wrinkle" density ,
which is the ame area per unit area resolved in the mean direction of propa-
gation. This enables more physically justiable modelling of the velocity eld
generated by the combustion, rather than the traditional density weighted en-
semble averaging and gradient transport modelling. Also, may be decomposed
as a spectral density function in length scale phase space, enabling detailed rep-
resentation of the turbulence- ame interactions processes [33].
This paper contains a detailed derivation of the \wrinkle" density trans-
port equation, and the corresponding equation for the regress variable (which
may be considered a normalised fuel concentration), from its origins as a ame
area model to a closed pair of coupled transport equations, and a simple one
equation model derived there-from. This derivation relies heavily on the use of
conditional-averaging techniques developed by Dopazo [34]. As a result of this
analysis the resulting transport equations involve conditional average properties,
such as the unburnt gas conditional velocity, which must be obtained for their so-
lution. The conditional averaged momentum equations are proposed as a method
of obtaining the required information, but their solution is considered unneces-
sarily complicated for the present work, and simple models for all the conditional
properties required are developed. The resulting models are analysed for stability
and compared with the simple EBU models. Finally, the models are solved in
simple one-dimensional and boundary layer test cases, the results again compared
with the EBU model. It is found that both the one and two-equation models are
inherently stable to both width and speed perturbations, and produce plausible
solutions of ame behaviour for all cases tested, unlike the EBU model which
failed in all cases except one dimensional propagation without heat release.
Chapter 2
Development of Flame Area
Models
The transport equation for ce is usually derived by ensemble or time averaging
the transport equation for a reacting species mass fraction [21, 2]. The Reynolds
ux term in the resulting equation is then modelled by applying the Boussi-
nesq approximation (\gradient transport") and reaction rate modelled either by
heuristic argument or on the basis of a presumed PDF. In the case of homoge-
neous turbulent combustion of a highly reactive species (high Damk}ohler number)
at moderate Reynolds number, the PDF may be assumed to be a double delta
function, i.e. the ame reaction zone is considered innitesimally thin, which
allows for various mathematical and modelling approaches to be applied to the
ame. For example the ame may be considered a sheet, the dynamics of which
may be represented and modelled through a transport equation for the ame area
per unit volume [29, 31, 32, 30]. Alternatively the PDF may be substituted into
the Reynolds averaged equations to represent combustion eects, together with
a representation of the burning rate [22, 23]. Both of these approaches may be
considered traditional in that they are well represented in the literature and are
becoming widely accepted as a basis for combustion modelling. However, if the
ame may be considered innitesimally thin, then a dierent approach to the
derivation of the averaged transport equations can be taken { that of conditional
averaging [35], as will now be described.
2.1 Conditional Averaging

Conditional averaging is applied to combusting systems by considering separately
the unburnt and burnt gases and the interface, i.e. reaction zone, see Fig. (2.1)
Equations are conditionally averaged by multiplying by an indicator function
I , which is 1 in one zone and 0 in the other, and then applying conventional
averaging techniques, density weighted or otherwise. This approach is a simple
2.1 Conditional Averaging 12
Unburnt
Burnt Su I=1
I=0
Figure 2.1: Flame Interface
extension of that applied to intermittent turbulent ows by Dopazo [34]. With

the use of the indicator function rather than the intermittency function, and by
allowing the phase interface to propagate, the analytical methods developed by
Dopazo can be applied to the present problem [34].
The terms resulting from conditional averaging are then manipulated by ap-
plying formulae presented in [34], and outlined below, replacing terms involving
the indicator function with equivalent terms involving the phase volume fraction
. A subscript will be used to denote the conditional average of a property in
the zone for which I = 1, i.e. for which is the volume fraction.
2.1.1 Formulae for Conditional Averaging

We dene the indicator function I (x; t) as:
(
if point (x; t) is in zone 1
I (x; t) = 01 if (2.1)
point (x; t) is in zone 0
The zone volume fraction is calculated as the probability of point (x; t) being
in zone 1 :
= I (x; t) (2.2)
where the overbar denotes the ensemble average.
2.1.1.1 Conditional Averaging Dierential Operations
Let Q(x; t) be any uid-mechanical property, scalar or tensor of any rank, then
rQI = r(QI ) = r( Q ) (2.3)
where Q denotes the conditional average Q in zone 1, denoted by subscript .
Also
@ QI = @ (QI ) = @ ( Q ) (2.4)
@t @t @t
2.1.1.2 Decomposition Operations

Consider the mean of I rQ = r(QI ) ? QrI over the control volume V . Note
rI is non-zero only at the interface, where it has the absolute value of the Dirac
delta function and the direction of the unit normal n to the interface, pointing
into zone 1. Then
Z
I rQ = r(QI ) ? Vlim 1 Qn (x; t) dS (2.5)
!0 V S (x; t)
where S (x; t) = 0 is the equation for the interface. We dene the surface average
Q of a property Q as the surface integral per unit volume divided by the surface
z}|{
area per unit volume , thus

R
limV !0 V1 S (x; t) Q(x; t) dS
Q (2.6)
z}|{

where
1 Z
= Vlim
!0 V S(x; t)
dS: (2.7)
Equation Eqn. (2.5) may then be re-written as
I rQ = r( Q ) ? Qn : (2.8)
z }| {
Also
I r Q = r ( Q ) ? Q:n (2.9)
z }| {
and
@ ( Q )
I @@tQ =
@t + QUs:n (2.10)
z }| {
where Us is the velocity of the interface.
2.1.2 Conditional Averaging Transport Equations

In what follows, the conditional averaging techniques outlined above are ap-
plied to the basic transport equations representing turbulent combustion at high
Reynolds and Damk}ohler numbers, by choosing zone 1 as the unburnt gas or
burnt gas and zone 0 as the other. In this way two transport equation are pro-
duced, for each one averaged, one for each zone. In addition to the transport
equations for mass, momentum, energy and turbulence properties, an equation
representing the average position of the interface is required. This is dealt with
rst.
2.1.3 Interface Transport Equation

Substitution of Q 1 into Eqn. (2.10) yields
@ ( ) = ?U :n z }| {
(2.11)
s
@t
which represents the rate of change of the unburnt gas volume fraction or the
burnt gas volume fraction depending on the choice of I and hence of . Clearly the
two equations are not independent due to the fact that the two volume fractions
sum to 1; it is therefore only necessary to solve one. Although the choice may be
made arbitrarily, the unburnt gas is chosen as the reference zone, this choice being
more convenient as will be shown later. Having made this choice the unburnt
zone volume fraction will be denoted as the regress variable b, following Pope [31],
and the direction of propagation nu will be denoted as n. Thus the transport
equation for b may be written
@b = ?U :n: z }| {
(2.12)
s
@t
The surface velocity Us may be decomposed into Uu the velocity in the unburnt
gas at the interface and Sunu, the velocity of advancement of an element of the
interface, relative to the unburnt gas, with nu pointing into the latter. Thus
Us = Uu + Sunu (2.13)
or equivalently in terms of the burnt gas properties
Us = Ub + Sbnb (2.14)
where
b Sb = ?uSu: (2.15)
and
nb = ?nu : (2.16)
At this stage the two decompositions are equivalent, but one or other, or a com-
bination, may prove advantageous at the modelling stage. For generality, the
choice will not be made now, and the surface velocity will be decomposed into
an interface uid velocity and propagation velocity in the direction of the normal
pointing towards the unburnt gas,
Us = UI + SI n (2.17)
Decomposing Eqn. (2.12) using Eqn. (2.17) gives
@b + U :n = ?S
z }| { z }| {
(2.18)
I I
@t
Decomposing UI :n into surface averages and surface uctuation correlations :

z }| {
UI :n = UI : n + UI ] :n]
z }| {
(2.19)
z }| { z }| {
z}|{
and substituting into Eqn. (2.18) results in

@b + U : n + U ] :n] = ?S :
z }| {
z}|{
z }| {
(2.20)
z }| {
@t I I I
Further simplication may be made to Eqn. (2.20) by substituting Q 1 into
Eqn. (2.8) resulting in
n = rb z}|{
(2.21)
which may be substituted into Eqn. (2.20) yielding
@b + U :rb + U ] :n] = ?S
z }| { z }| {
(2.22)
z }| {
I I I
@t
In the expectation that the surface uctuation correlation term will be modelled
as some kind of diusion term, it will be denoted Db, in which case the transport
equation for b may be written
@b + U :rb + D = ?S
z }| {
(2.23)
z }| {
I b I
@t
where
Db = UI ]:n]
z }| {
(2.24)
2.1.4 A General Denition for Wrinkle Factor

The ame area per unit volume may be decomposed into the ame wrinkle
factor , the ame surface area per unit projected area, and a simple function
of b, using conditional averaging techniques. The utility of this approach will
become clear as the analysis proceeds.
Consider an ensemble of turbulent ames created in an ensemble of experi-
ments all conducted under the same conditions, Fig. (2.2). Consider a volume
element V containing facets of the ensemble of ames. Then the ame wrinkle
factor is dened as
Af
z }| {
(2.25)
As
z }| {
where Af is the average ame area per unit volume in the control volume V , i.e.
z }| {
1 Z
Af = Vlim dS
z }| {
!0 V S(x; t)
= (2.26)
δV
Contours of b
Figure 2.2: Ensemble of Flame Interfaces
and As is the average ame area projected onto the mean propagation direction
z }| {
n^ per unit volume in the control volume, i.e.

1 Z
As = Vlim n^ :n dS
z }| {
!0 V S(x; t)
= j n j: z}|{
(2.27)
From Eqs. (2.26, 2.27 and 2.21) may be expressed as
= 1 (2.28)
jnj
z}|{
= (2.29)
jrbj
and
rb
n^ n^ = jr
z}|{
(2.30)
bj
2.1.5 Transport Equation for b in terms of
Decomposing the source term in Eqn. (2.23) using Eqn. (2.29) gives
@b + U :rb + D = ?S jrbj:
z }| { z }| {
(2.31)
I b I
@t
The signicance of this equation does not become apparent until jrbj is expressed
as n^ :rb, in which case
@b + U :rb + D = ?S ^n:rb;
z }| {
(2.32)
z }| {
@t I b I
the right hand side of which is clearly a propagation term, with SI the turbulent
z }| {
ame propagation speed and n^ the direction of propagation. The meaning of the
ame wrinkle factor is now apparent: it is the turbulent/laminar ame speed
ratio. For convenience a total transport velocity is dened as
Ut = UI + SI ^n (2.33)
z }| { z }| { z }| {
in which case Eqn. (2.31) may be rewritten as

@b + U :rb + D = 0
z }| {
(2.34)
@t t b
2.1.6 Transport Equation

A transport equation for may be derived by taking the time derivative of Eqn.
(2.29) at position x :
! ! !
@ = 1 @ ? 2 @ jrbj
@t x jrbj @t !x jrbj @t ! x
= 1 @@t ? @ jr bj (2.35)
jrbj x jrbj @t x
Further development of this equation requires expressions for the time derivatives
of jrbj and .
2.1.6.1 jrbj Transport Equation
jrbj may be decomposed and dierentiated with respect to time thus
q
jrbj = rb:rb (2.36)
! !
@ jrbj = p1 @ rb:rb
@t x 2 rb:rb ! @t x
= rb :r @b
jrbj !@t x
= n^ :r @b
@t x (2.37)
Taking the gradient of Eqn. (2.34) gives

!
@b
r @t = ?r(Ut :rb) ? rDb
z }| {
x
= ?Ut :rrb ? (rUt):rb ? rDb (2.38)
z }| { z }| {
Combining Eqn. (2.37) with Eqn. (2.38) and substituting Eqn. (2.30) produces
!
@ jrbj = ?U :rrb:n^ ?n^ :rU :rb ?n^ :rD :
z }| { z }| {
t t b (2.39)
@t x
I II III IV
Noting that
r(rb:rb) = 2(rrb):rb (2.40)
term II in Eqn. (2.39) may be rewritten
II = ? 21 Ut : r(rb:rb) : (2.41)
z }| {
jrbj
From Eqn. (2.29)
2
rb:rb = (2.42)
in which case
2
r = 2 r
= 2 jrbj r ? jrbjr :

(2.43)
Combining Eqn. (2.39), Eqn. (2.41) and Eqn. (2.43) gives
!
@ jrbj = ? Ut : r ? jrbjr ? n^ :rUt :rb ? n^ :rDb
z }| {
z }| {
@t x (2.44)
which concludes the derivation of the transport equation for jrbj. Substitution
of Eqn. (2.44) into Eqn. (2.35) yields
!
@ + U :r = ^n:rU :n^
z }| { z }| {
@t x t t
!
1 @ + Ut :r
z }| {
+
jrbj @t x jrbj
+ ^n :rDb : (2.45)
jrbj
2.1.6.2 Transport Equation

The derivation of the transport equation for has been published in several
papers [31, 32, 30] and need not be reproduced here. The result, here stated in
non-conservative form is
@ + U :r = ?n:rU :n:
z }| {
(2.46)
z }| {
@t s s
2.1.6.3 Transport Equation Assembly

Substitution Eqn. (2.46) into Eqn. (2.45) produces
@ + U :r = ^n:rU :n^ ? n:rU :n
z }| { z }| { z }| {
t t s
@t
+ Ut ? Us : r
z }| { z }| {
jrbj
+ ^n :rDb : (2.47)
jrbj
After substituting for from Eqn. (2.29) and rearranging, there results
@ + U :r = ? n:rU :n
z }| {
I
z }| {
s s
@t
+ ^n:rUt :n^ II
z }| {
rjrbj
+ Ut ? Us :
jrbj III
z }| { z }| {
+ ^n :rDb IV:

jrbj
Term I on the r.h.s. may be decomposed into a contribution involving the surface
averages of Us and n and residual terms, representing all the eects of turbulence
and the distribution of interface direction on . Performing this decomposition,
noting from Eqn. (2.28) that
n = j n jn^ = n^
z}|{ z}|{
(2.48)
yields
n:rUs :n = 1 n^ :rUs:n^ ? G (2.49)
z }| { z }| {
where
z
G = ? n] :rUs:n] + n] :rUs] : n + n :rUs ]:n]
}| { z }| { z }| {
z}|{ z}|{
(2.50)
Substitution of Eqn. (2.49) into Eqn. (2.48), and renaming term IV as D , in

the expectation that it will form a diusion like term for , yields
!
@ + U :r + D =
z }| {
G I
@t x s
+ ^n:rUt :n^ ? 1 n^ :rUs:n^ II

z }| { z }| {
+ Ut ? Us : rjrbj III :
z }| { z }| {
jrbj
The rate of change and transport terms on the l.h.s. of the above are self ex-
planatory, but note the velocity transporting is the surface-averaged interface
value. The analysis of D is deferred to the modelling stage of the equation
development.
Term I on the r.h.s. represents turbulence and interface uctuation eects,
resulting in an increase or decrease of the interface area. This term may be
identied with the modelling of ame interface behaviour presented in an earlier
paper by the author [33].
Term II on the r.h.s. of Eqn. (2.51) represents the mean stretch and prop-
agation eects on , the rst part of which is the generation due to meam ow
stretch and the interaction of propagation with mean curvature and the second
part correcting for these eects on the projected area. If this correction were not
present the ame \wrinkling" could be reduced even if the surface were smooth.
The signicance of term III is not immediately obvious, however a clue to the
behaviour is given by the velocity part, which represents the dierence between
the overall propagation velocity and the average interface velocity, which increases
with interface distortion. The term in b may be best interpreted by constructing
the distribution through a hypothetical Gaussian ame Fig. (2.3). The choice of
a Gaussian is considered reasonable on the basis of experimental results for rod
stabilised \V-shaped" ames [36, 37]. From Fig. (2.3) it can be seen that the
function of b in term III tends to +1 at the back of the ame and ?1 at the
front. Thus at the front of the ame term III has the eect of smoothing out
the ame by forcing ! 1. Note cannot become less than 1, by virtue of
the fact that the dierence between the two propagation velocities is zero for a
smooth ame. This point becomes clear if the velocity dierence is evaluated in
terms of using Eqs. (2.33 and 2.17). Surface averaging of Eqn. (2.17) creates
a correlation between the laminar ame speed and the local interface normal.
Neglecting this correlation and applying Eqn. (2.48) produces
Us UI + SI n^ (2.51)
z }| { z }| { z }| {
in which case the velocity dierence becomes

z }| { z }| { 1
Ut ? Us = ? SI n^
z }| {
(2.52)
0
n.x
b
(Gaussian)
b
b
Figure 2.3: Construction of Term III for a Gaussian Flame

which ! 0 as ! 1. At the back of the ame there is no similar eect to

counter-balance the rise in term III , in which case ! 1 at the back of
the ame. From this description it is clear that the term being described has
the properties of \cusp" formation, i.e. smoothing the front of the ame and
forming points of innite curvature at the back. The cusp formation process is
illustrated in Fig. (2.4), which is a simulation of an initially plane ame interface
propagating into a single vortex (which has a Taylor type velocity distribution).
The smoothing of the front of distortion and the formation of the \cusp" can be
clearly seen late into the interaction.
Sur face Distr ibution
2.0
1.6
time
1.2
0.8
y
0.4
0.0
- 0.4
- 6.0 - 4.0 - 2.0 0.0 2.0 4.0 6.0
x
Figure 2.4:
\Cusp" Formation for a Laminar Flame in a Taylor
Vortex
2.1.7 Modelling the b Transport Equation

Although the b and transport equations developed in the previous section are
largely closed, there remain a few terms which require modelling, these being
terms involving surface average velocities and correlations between the interface
normal and the local velocity eld or the propagation speed. In the present anal-
ysis it is additionally assumed that conditionally averaged transport equations
are solved in order to obtain the properties of the unburnt and burnt gases. Possi-
ble modelling of these properties in terms of density weighted ensemble averaged

values is deferred to a later section.
The problem which remains is the relationship between the surface average
interface uid velocity and the conditional averages. Such a relationship must
involve the interface distribution, as this must be correlated with the velocity eld
in that the distortions are caused by the turbulence content of the surrounding
gas. As already mentioned, there is experimental evidence that the distribution
of jrbj within the ame is close to Gaussian, in which case it is expected that
\gradient transport" approximations may be used to represent the generation of
the ame distribution. Thus, extending the work of Cant & Pope [32] on constant-
density ames, the surface-average interface uid velocity is decomposed into the
interface conditional average, UI , and a surface correlation term
UI UI ? DI r : (2.53)
z }| {
where the diusion coecient DI will be modelled in terms of the conditionally

averaged turbulence properties.
The interface conditionally averaged velocity | interface normal correlation
term in Db, if treated similarly, results in
UI ] :n] ?DI r n
z }| {
z}|{
(2.54)
which may be decomposed using Eqs. (2.28 and 2.29) to give
DI r n = DI r n^
z}|{

r
= ?DI :rb + DI (r n^ )^n:rb: (2.55)

The last term of Eqn. (2.55) represents the correlation between the velocity uc-
tuations and the ame mean curvature. Such an eect is likely to be negligible in
all cases except perhaps early kernel growth. For clarity this term is neglected in
the present work. Combining Eqs. (2.53 and 2.55) to form the complete modelled
form of Db and substituting into Eqn. (2.34) provides the transport equation for
b in closed form
@b + U :rb = ?S jrbj: z }| {
(2.56)
I I
@t
or equivalently, dening a suitable total transport velocity UT ,
@b + U :rb = 0 (2.57)
T
@t
where
UT UI + SI ^n (2.58)
z }| {
It is clear from Eqn. (2.57) that the use of \gradient transport" modelling
for the surface uctuation correlations has provided an equation which does not
contain explicitly a dispersion process for the interface: however such information
is contained in the conditional averaged velocity. Also Eqn. (2.57) has the form
of a wave equation and hence is bounded and non-conservative, which is to be
expected since the transported property b need not be conserved but must be
bounded between 0 and 1. Another important feature of Eqn. (2.57) is that the
ame burning velocity eigenvalue is simply SI , i.e. the turbulent burning veloc-
z }| {
ity resulting from the modelling of . Unfortunately the ame width eigenvalue
cannot be evaluated simply from the form of the equation as it results from the
interaction of the combustion with the unburnt gas.
There remains the problem of choosing the zone or zone averaging used in
specifying UI and SI in terms of the known conditional average velocities. There
z }| {
is evidence from direct simulation studies that the creation of ame distortion is
dominated by the turbulence in the unburnt gas, due to the density change across
the interface \ejecting" the vortices from the burnt side preventing interaction
with the ame [38]. On this basis one might expect that DI Du. However, it
seem unlikely that the interface behaviour can be represented purely in terms of
the unburnt gas properties: at least the interface uid velocity must be depen-
dent on both the unburnt and burnt gases. An alternative approach proposed by
Dancy [39] and developed by Chen et al. [35] is to postulate that uctuations of
the interface properties, about the interfacial point average, are uncorrelated and
the interface average can be approximated by the conventional average. Although
this may be a reasonable approximation for the case of no density change across
the interface, as in [35], the application to the general combustion case is unrea-
sonable. Consider a \simple" ame interface containing cusps, Fig. (2.5). At the
back of the ame the unburnt gas is partially trapped by the \cusps" and one
would expect that the local interface properties will be highly correlated with the
unburnt gas conditional properties. Conversely, the burnt gas is predominantly
in the \free" gas surrounding the ame and unlikely to be strongly correlated
with the interface properties. A similar condition exists at the front of the ame,
except that the interface has a lower curvature than at the back and hence is less
contained.
These arguments suggest that a form of inverse averaging (in which the value
of the interface average, as the unburnt gas fraction approaches 1, is dominated
by the burnt gas conditional properties and vice versa) should be applied to the
conditional averages in order to construct the interface transport conditional av-
erage. It is considered appropriate to average the momentum (density weighted
velocity) at the interface rather than the velocity, as the former varies continu-
ously across the interface and the later discontinuously. A second advantage in
burnt
gas
Ub
Ub
partially Uu
trapped
unburnt gas
unburnt
Su gas
Ub Ub
Figure 2.5: Modelling Surface Averages
the use of density weighting is that the unburnt gas properties will dominate the
interface behaviour as suggested by direct simulation [38]. The density weighted
interface uid velocity is constructed thus
I UI (1 ? b)uU b Ub
gu + b g (2.59)
where
I = (1 ? b)u + bb: (2.60)
Similarly
I SI = uSu
z }| { z }| {
(1 ? b)uSu ? bbSb (2.61)

z }| { z}|{
and
I DI (1 ? b)uD
gu + b D
f
b b (2.62)
2.1.8 Modelling the Transport Equation

The modelling methods developed in the previous section may equally well be
applied to the transport equation for , Eqn. (2.51). The transport velocities
Ut and Us are modelled as in the previous section, leaving the diusion term,
z }| { z }| {
D, and the turbulence interaction term G . The former comprises a combination
of the diusion term of the b equation and the uctuation of the interface uid
velocity part of the convective transport term, the models for which may be
combined to create the diusion term for the transport equation

D = Us ? Us :r ? ^n :rDb

z }| {
jrbj !
r
= ?DI :r + n
^ n
^
jrbj :r jrbjDI r : (2.63)
where
Us = UI + SI n^ (2.64)
z }| {
If, for the purposes of the modelling of diusion in the transport equation,
mean surface curvature eects are neglected, the second term on the r.h.s. of
Eqn. (2.66) becomes
!
n
^
r jrbjDI r = jrbj r jrbjDI r 1 (2.65)
Combining Eqn. (2.63) and Eqn. (2.65) forms
D = ?r DI r ? DI r: rjrbj (2.66)

jrbj
the rst term of which is a \gradient transport" term in standard conservative
form and the second term has a form similar to the \cusp formation" term de-
scribed earlier. The latter is of particular interest in that at the back of the
ame the term in b is large and positive and the gradient of is expected to be
negative throughout the ame, thus the contribution of this term is such as to
counter-balance the eect of the \cusp" formation, i.e. \cusp" dissipation due
to turbulence interaction. It is important to note that the ame width is ex-
pected to be dictated by the balance between the eects of \cusp" formation and
dissipation.
As mentioned earlier, the turbulence interaction and ame propagation eects
on the ame \wrinkling", collected into G , are the subject of a previous study,
[33], in which the surface is represented as a spectrum of distortions, and all the
processes aecting the spectrum are modelled individually and combined into a
set of integro-dierential equations in length scale space. The solution of this
equation set has been accomplished for pseudo-homogeneous ames, producing
turbulent ame speeds and spectral distributions of many ame statistics. The
length scale integral of the net \wrinkling" generation rate may be identied with
the term G in the transport equation for i.e.
Z 1 @ (l) !
G = 0 @t dl: (2.67)
x
It has not proved possible to integrate analytically the spectral equations to

create a simple model for G , however, it is possible to perform the opposite, that
is, dierentiate the transport equation for with respect to length scale. The
resulting transport equation for the ame area density function (l) would be
extremely costly to solve and this development is left to future work. For present
purposes a simple model is sought which has the correct qualitative behaviour
and is numerically stable. From the spectral work [33] and direct simulation
[40, 41, 42] it is found that the rate of generation of area per unit area of a
passive surface in homogeneous turbulence is approximately 0:28= . Also, it is
found that the rate of removal of the area per unit area due to propagation is
proportional to the turbulent ame speed. Combining these eects the following
model is proposed
G = G ? R2 (2.68)
where
G = 0:28 (2.69)

and the removal rate coecient R is as yet unknown. Consider a case where the
turbulent ame speed is known, and the ame speed is measured as the speed of
the point in the ame where rjrbj = 0. At this point the dominant terms in the
equation, \cusp" formation and dissipation, are zero; also the other terms are
small due to the low gradient of . Under these conditions the generation and
removal rate must balance at equilibrium, the resulting equilibrium , eq , being
equal to the turbulent to laminar ame speed ratio. Thus,
R = G=eq (2.70)
where
eq = St : (2.71)
Su
z }| {
eq may be specied in a wide variety of ways, i.e. from the solution of the spectral
equation, from correlations of experimental data, or from simple models of ame
propagation e.g. fractal models. For the purposes of the present investigation, a
simple conventional linear relation with the turbulence intensity u0 will be used:
0
eq = 1 + u (2.72)
Su
z }| {
2.2 Conditionally Averaged Velocities 28
Collecting the terms, the nal modelled equation is

!
@ + U :r ? r D r =
@t x s I
+ G ? R2
+ ^n:rUt :n^ ? 1 n^ :rUs:n^
z }| { z }| {

+ SI ? + DI r : rjrbj :
1 (2.73)
z }| {
jrbj
It is interesting to consider now, the dierences between the development
of the transport equation and the modelling of the transport equation by
the various developers [29, 43, 32]. The present analysis does not add any new
information to the original form of the transport equation [31], however, the
manipulation to convert into has separated out a term which represents ame
annihilation by \cusps", a term which does not require modelling. Whereas, in the
modelling of the equation, all the work to date has involved the introduction
of a modelled term representing this important feature of ame behaviour. It
is concluded that while it is necessary to introduce a term to represent ame
area removal by propagation (R2 in Eqn. (2.73)), it is unnecessary to give this
term spatial dependency (as in the work of Marble and Broadwell [29] and all
subsequent developments), such that a stable ame width may be obtained. This
function is fullled by the dierence in propagation behaviour of , and hence of
, manifesting itself as a \cusp" term.
A comparison may also be draw between the transport equation for b, Eqn.
(2.23), and that for G in the work of Peters [44]. Clearly, if G is considered an
indicator function, rather than a continuous scalar eld, the G and b are identical.
However, in the work of Peters, G is not given such a concrete denition, and the
interpretation of G as an indicator function is open to question. The use of condi-
tional averaging in the present work introduces the concept of surface averaged as
well as conditionally averaged properties, which requires a form of closure quite
unlike the more traditional moment closure, as used in the G equation. Finally,
the way in which surface eects are treated in the two approaches are quite dier-
ent. In the G equation, surface stretch and curvature eects are treated directly
through the mean and variance of G, whereas in the present work these eects are
represented through the ame \wrinkling" and the physical processes involved in
it's development. Thus, although the starting point of the equation derivation
may be considered similar, the modelling and closure are quite dierent.
2.2 Conditionally Averaged Velocities

Although a ame area based combustion model comprising of Eqs. (2.15 and
2.73) has now been formulated, its utility still relies on the availability of the
conditional average velocities. Although these may be obtained by solving the

conditional average momentum and continuity equations, together with the zone
equations of state, this approach is complicated and expensive; and a model
based on the density-weighted ensemble-averaged properties may be more appro-
priate, if it proves to be suciently accurate. However, before any such modelling
can be undertaken, it is necessary to develop and close the complete conditional
transport equation system, comprising of continuity, momentum, energy, turbu-
lence and state equations, using the conditional averaging techniques developed
in Section 2.1 from which simple models may be derived.
2.2.1 Conditionally Averaged Continuity Equations

The continuity equation
@ + r U = 0 (2.74)
@t
is conditionally averaged by multiplying by the indicator function I , averaging,
using a conventional method, and decomposing using the expressions presented in
Section 2.1. The conventional technique chosen here is that of density-weighted
averaging, as this avoids the creation of all the density uctuation correlations.
Density weighting in conditional averaging seems more reasonable than when
used to average across the ame structure, as the density variations within each
zone are an order of magnitude less than that across the ame interface. Also
the use of density weighting is common in the development of turbulence models
for variable density ows [45] on which it is intended to base the present work.
Thus applying density-weighted conditional averaging to Eqn. (2.74) and making
use of Eqs. (2.3, 2.4, 2.8, 2.9 and 2.10) there results
@
+ r ( U ) = ? S (2.75)
f
z }| {
@t
where the subscript represents conditional average for the zone and the tilde
represents the use of density weighting.
2.2.2 Conditionally Averaged Momentum equations

The momentum equation
@U + r UU + r = ?rp (2.76)
@t
where

= ? rU + rU T (2.77)
z
@t
e
@ U
}| {
z }| {
^
is conditionally averaged in the same manner as the continuity equation, yielding
+r ( U U ) + r = ?r p
?US +n : ? r (Un + n U) + pn
z }| {
z }| {
(2.78)
where

= ? r U + r U T : (2.79)
The rst line of Eqn. (2.78) resembles the density-weighted ensemble averaged
momentum equation except for the presence of the zone volume fraction in
all the terms. The second line however does not have any counterpart in the
^
conventional averaged form and represents the eect of the interface on the zone
momentum. Apart from the various terms in the laminar viscosity there remains
the problem of decomposing and modelling the velocity correlation U U and
the surface average terms involving velocity and pressure. The rst of these may
be decomposed into terms involving the average and a uctuation correlation in
^
the same manner as ensemble averaging
U U = U
f U ^
f + U00 U00
(2.80)
the second term of which requires either modelling or the solution of a suitable
transport equation. It is considered that, in these preliminary studies of the use
of conditional averaging for combusting systems, the solution of Reynolds stress
equations is unnecessary and the use of \gradient transport" models for the zone
property correlations is adequate. This judgement is based on the earlier obser-
vation that the density variations within the zones are an order of magnitude
less than those across the interface, in which case the eect of \counter-gradient
transport" observed in ames might be adequately represented through condi-
tional averaging without recourse to second order closure. However, it is noted
that for complex ow problems involving curvature and adverse pressure gra-
dients \gradient transport" approximations have proved inadequate, for which
cases, with combustion, it is also expected that this approximation will prove
inadequate. However, such analysis is best deferred until the principles presented
here are proven. The application of \gradient transport" to Eqn. (2.80) results
in
^
U00 U00 = ?D rUf + rUfT ? 2 (r U
3

f )I + 2 kI
3 (2.81)
where I is the identity matrix.
The surface-average term in Eqn. (2.78) involving velocity may be decom-
posed, noting that S is constant across the interface (see Eqn. (2.15)), in which
case uctuations of this property may be considered small compared with uctu-
ations in the interface velocity which have the order of magnitude of the density
change across the interface, thus
I UI SI I SI UI : (2.82)
z }| { z }| {z }| {
However, it is dicult to dene unambiguously the uid velocity at the ame in-

terface, due to it's change across the interface; and it seems reasonable to suppose
such a denition should be linked to the denition of the interface pressure, which
drives the change in momentum on either side. Using the approach developed in
Section 2.1.7, UI is decomposed into a conditionally averaged part and a term
z }| {
representing surface uctuations, in which case the consistent choice of average

interface pressure is
pI = (1 ? b)pu + bpb: (2.83)
Decomposing the interface pressure | interface normal correlation
pn = p n^ + p] n ]
z }| {
z }| { z}|{
(2.84)
and neglecting all pressure uctuation correlations yields
pn pI n^ :
z }| {
(2.85)
Combining the above models and neglecting terms in the laminar viscosity,
the assumption being that the Reynolds number is high enough for turbulent
transport to dominate, results in the following closed equation for the conditional
average momentum:
e
@ U h i
@t + r ( U
f U
f ) ? r D rU

f + rUfT ? 2 (r U
3
f )I + 2 kI =
3
? rp ? (p ? pI )r ? I SI UI (2.86)
z }| {z }| {
The pressure terms have been decomposed and combined to form the pressure
driving force within the zone and that due to the pressure dierence between
the interface and the zone. The last term on the right hand side represents
the transfer of momentum between the zones due to the gas transfer by ame
propagation.
2.2.3 Solution of the Conditionally Averaged Equations

In the previous two sections the conditionally averaged continuity and momentum
equations were derived and modelled in order to obtain a closed system. These
together with the zone equations of state
p = R T (2.87)
2.3 Ensemble Averaged Equations 32
and the transport equations for the zone energies, from which the zone temper-
atures are obtained, may be solved to obtain the conditional average properties.
These are then used to construct the interface average properties which in turn
are used to solve for the interface structure, b and . This equation set should be
capable of describing much of behaviour of premixed turbulent ames, but is far
more complex than the ensemble-average equations set presently being used to
describe combustion in engineering calculations [46]. It may be that such com-
plexity is unnecessary for most engineering purposes, and that one need only solve
for the ensemble average properties in order to adequately represent turbulent
ame behaviour. This possibility will now be explored.
2.3 Ensemble Averaged Equations

2.3.1 Transport Equation for be
The proposed simplied system of equations should rstly preserve fuel mass
conservation, as this is the property governing the heat release from the ame,
secondly be bounded, and lastly have the same transport form as the density
weighted ensemble averaged momentum and enthalpy equations. This last re-
quirement results from the study of coupling between the fuel, enthalpy and
momentum equations: any dierences in the transport form of which creates
instabilities in the ame structure and can lead to unphysical propagation.
One approach would be to solve a conservation equation for the density-
weighted regress variable eb which is directly related to the fuel mass fraction
through
fe ? fb
eb = (2.88)
fu ? f b
and may be expressed in terms of the unweighted value through [23]

eb = ub
1 ? eb = b 1 ? b

= ub + b 1 ? b (2.89)
Eqn. (2.75) may then be written
@eb + r U eb = ? S jrbj
z }| {
(2.90)
u u u
@t
which is in conservative form with respect to the unburnt gas but not with respect
to the ensemble average properties. Also, due to the presence of the unburnt gas
density and velocity in this equation, it cannot be expressed in bounded form by
combination with either the ensemble averaged or conditional averaged continu-

ity equations, a technique commonly used for the ensemble averaged transport
equations.
Another approach would be to solve the equation for the unburnt gas volume
fraction b, Eqn. (2.56)
@b + U :rb = ?S z }| {
(2.91)
I I
@t
which is bounded but requires the average interface uid velocity UI and is not
in conservative form.
Considering the importance of conservation and the uncertainty in modelling
the average interface uid velocity, it is considered preferable to begin modelling
from the conditionally averaged unburnt gas continuity equation, Eqn. (2.90),
and construct a simple model for the conditional average velocity Uu.
The ensemble average velocity may be represented in terms of the conditional
averages through [23]
U
f = bUu + (1 ? b)Ub
u b (2.92)
which may be rearranged in terms of the slip velocity Uub
Uu = U
f + (1 ? eb)Uub (2.93)
where
Uub = Uu ? Ub: (2.94)
By analogy with the properties of laminar ames, the velocity dierence across a
turbulent ame is a function of the turbulent ame speed and the density ratio,
in which case performing a simple one dimensional analysis of the turbulent ame
brush there results
!

Uub ? 1 S u^n:
u z}|{
(2.95)
b
However, this does not correctly account for the distribution of , the assump-
tion being that is constant through the ame, nor curvature eects or turbulent
uctuations in the phase velocities. Noting that required transport form of the
resulting equation must be the same as that of the rest of the ensemble averaged
transport equations, a model for the turbulence eect is postulated, for consis-
tency, as a simple gradient diusion form, producing the following result for the
slip velocity:
!
Uub u ? 1 S u^n ? De e rb e :
e
(2.96)
z}|{
b b(1 ? b)
where the turbulent diusion coecient De is the same as that used in the enthalpy
transport equation. Substitution of Eqs. (2.96 and 2.93) into Eqn. (2.90) yields
" ! #
@eb + r (U e b reb) = ? Su jrbj ? r (1 ? eb) u ? 1 S u ^
feb) ? r (D
z }| {
n
z}|{
@t t u b (2.97)
which may be simplied by rearranging the second term on the r.h.s using Eqn.
(2.89), in the form
!

(1 ? eb) u ? 1 = u ? (2.98)
b
resulting in
@eb + r (U e b reb) = ? Su jrbj
feb) ? r (D
z }| {
@t t u

? r (u ? ) eb S u^n : (2.99)
z}|{
Decomposing the last term of Eqn. (2.99) using Eqn. (2.89) produces

r (u ? ) eb S u^n = S u^n:ru eb ? b + u eb ? b r S u^n
z}|{ z}|{ z}|{
(2.100)
which, when combined with Eqn. (2.99) nally results in
@eb + r (U e b reb) = ? Su jrebj
feb) ? r (D
z }| {
@t t u

? u eb ? b r S u^n
z}|{

? eb ? b S u^n:ru (2.101)
z}|{
The last term on the r.h.s of the above will be small in comparison with the
rst due to the gradients of the unburnt gas properties being small compared to
the ensemble-averaged properties within the ame. However the second term on
the r.h.s will be signicant and dominated by the gradient of at the back of
the ame produced by \cusp" formation. This term creates a problem with the
equation due to the requirements for conservation and boundedness. If removed
the equation will be bounded but non-conservative, and if retained the equation
will be conservative and unbounded. The boundedness requirement is crucial
due to the need for a smooth distribution of b in order to evaluate the second
derivatives in the equation, apart from the realisability requirement. If conser-
vation is violated in this equation but corrected for by constructing a consistent
fuel consumption rate for the energy equation, the ame will not propagate at
the correct speed. In the expectation that such an error in speed will be small
or correctable, the bounded approach is chosen. However, the removal of the

gradient of term will destabilise the ame as b ! 0, by allowing the back of
the ame to accelerate as increases due to \cusp" formation. It is therefore
desirable to compensate for this eect by adjusting at the back of the ame.
This adjustment is only required with respect to the density change: for a ame
with no heat release no adjustment would be required. It is thus proposed to set
an upper limit of at the equilibrium value, see Eqn. (2.71), which is simple,
stable and allows the ame to propagate at the correct speed. A problem with
this method occurs in complex ow systems where there may be area generation
due to mean ow eects, in which case eq as specied would not be correct. For
the present analysis this limitation will not be addressed. The nal, modelled,
transport equation for eb and corresponding heat release term for the enthalpy
equation are then:
@eb + r (U e b reb) = ?( + ( ? ) min(; eq ))Su jrebj
feb) ? r (D
z }| {
@t t u
(2.102)
R_ h = ?c H (fu ? fb ) ( + (u ? ) min(; eq ))Sujrebj

z }| {
(2.103)
2.3.2 Transport Equation for
The transport equation requires little alteration, but the interface uid velocity
UI must be expressed in terms of ensemble averaged properties, consistent with
the modelling of the eb equation. The interface uid velocity may be expressed
in terms of the density weighted ensemble average and the \slip" velocities by
manipulating Eqs. (2.59, 2.60 and 2.89), resulting in
!
f + (1 ? eb) ? b eb Uub :
UI = U (2.104)
I u
where
!
I = u + eb b ? 1 (2.105)
u
and the \slip" velocity is modelled as Eqn. (2.96).
2.3.3 A One Equation Model

For some engineering purposes it may not be necessary to reproduce all the
details of the ame structure or the transient response of the turbulent ame
speed. Under these conditions the transport equation need not be solved at
all; instead a simple analytical expression for can be used in the eb equation.
For example substituting = eq (assuming local equilibrium) into Eqn. (2.102)
results in the following simple transport equation for eb:
@eb + r (U
feb) ? r (De b reb) = ? eq Su jrebj
z }| {
(2.106)
@t t u
Such a model will not represent \cusp" formation and dissipation processes which
are responsible for maintaining a stable ame width: thus the ame may increase
in width indenitely, and hence the model may only be applicable to small ge-
ometries, for example internal combustion engines. The lack of transient response
could be remedied by obtaining eq from analytical expressions for the turbulent
burning speed, for example those developed to represent ame kernel growth in
internal combustion engines [47, 48]. However, in the present work, such com-
plexity is not warranted and the following simple linear expression for eq will be
used in both the one and two equation models:
0
eq = 1 + A u : (2.107)
Su
z }| {
Chapter 3
Assessment of the b { Models f
In order to study the behaviour of the one and two equation models described,
Eqs. (2.102, 2.73 and 2.106), simple one-dimensional test cases are used, from
which information about the ame structure and stability is obtained. The results
are compared and contrasted with those of the simple \eddy break-up" model
Eqn. (1.4) analysed in the introduction. Initial tests are performed in an open
duct, see Fig. (3.1), with the ame stabilised in the centre by adjusting the inlet
velocity to balance the propagation, which allows for good spatial resolution
and cheap computation. The conditions of the case are give in Table 3.1. For
conditions under which a stable ame could not be obtained in this geometry, a
semi-closed duct is used, see Fig. (3.2) the ame being initiated at the left-hand
end, which is closed, and then propagates towards the open end. This allows the
form of the instability to be analysed in order that the range of applicability of
the model may be ascertained. The conditions for the semi-closed duct case are
also as in Table 3.1 except that the unburnt gas velocity need not be specied.
Flame
Ub Position Uu
Constant Specified
Pressure
Figure 3.1: One Dimensional Open Duct Test Case
Assessment of the eb { Models 38
Property Symbol Value Dimension

Initial pressure p 1.0 bar
Initial temperature T 300.0 K
Fuel Methane
Equivalance Ratio 0.8
Laminar Flame Speed Su 0.3 m=s
Density Ratio u=b 6.72
Turbulence Intensity (xed) u0 1.2 m=s
Integral Length Scale (xed) l = u0k=" 0.01 m
Laminar Viscosity 15:6 10?6 m2 =s
Derived Properties
Turbulent Kinetic Energy k 2.16 m2 =s2
Turbulent Dissipation Rate " 259.2 m2 =s3
u0=Su 4
St=Su = 1 + u0=Su 5
Inlet velocity Uu = ?St -1.5 m=s
Equilibrium Wrinkling eq = St=Su 5
Table 3.1: One Dimensional Test Conditions.
ρu
S
ρb u
Constant
Pressure
Figure 3.2: One Dimensional Semi-Closed Duct Test Case
3.1 Constant Density 39
3.1 Constant Density

With no heat release, the eb { combustion model, Eqs. (2.102 and 2.73), produces
a ame which stabilises with a width of approximately 0.5l, dened as the distance
between the peaks in the Laplacian of b, and propagates close to the expected
turbulent ame speed Table 3.1, as specied through the equilibrium \wrinkling"
relationship, Eqn. (2.107), see Fig. (3.3). The reason why the speed is not exactly
1.0
0.8
0.6
b
0.4
0.2
0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
x/l
Figure 3.3: eb { Model b Prole
as specied is that the point in the ame where the \cusp" terms are zero does
not correspond to the point where the diusion term in the \slip" velocity is zero,
see Eqs. (2.73 and 2.96). This dierence causes \wrinkling" from the back of the
ame to diuse into the ame increasing the speed of propagation. This problem
could be remedied by reformulating the diusion term as in the work of Cant
and Pope [32], which in this case, where the surface correlation term is already
accounted for, has the form
?De rjr
jrbj
bj : (3.1)
However, it is dicult to extend this form of diusion modelling to account for the
density eects at the ame interface, accompanying combustion with heat release.
Alternative diusion representations are possible, which may be formulated to
allow for the correct propagation speed but their development and testing is left
to future work. Suce it to say that the error in the propagation speed is of the
order of 1% and hence is not considered critical.
The prole, Fig. (3.4), and ame width are very similar to those produced
by Cant and Pope [32], which exhibits a smooth Gaussian like shape at the front
of the ame and a much steeper prole at the back. This is indicative of the
9.0
8.0
7.0
6.0
2
5.0
Σ /10
4.0
3.0
2.0
1.0
0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
x/l
Figure 3.4: eb { Model Prole
presence of \cusps" at the back of the ame which can be seen more clearly
in the distribution, plotted against b in Fig. (3.5), which is monotonic and
increasing rapidly at the back tending to innity. The long leading edge in the
20.0
16.0
12.0
Ξ
8.0
4.0
0.0
0.0 0.2 0.4 0.6 0.8 1.0
b
Figure 3.5: eb { Model Prole
prole can be attributed to the gradual decrease in . This quantity does tend
towards 1 at the front of the ame but only at the very tip, behaviour which can
be attributed to the smoothing of the \wrinkles" due to propagation. The balance
between the \cusp" formation and dissipation terms appearing in Eqn. (2.73) is
presented in Fig. (3.6) If the diusion coecient in the \cusp" dissipation term
is doubled the ame becomes correspondingly wider see Fig. (3.7), which lends
weight to the physical interpretation of the term as opposing \cusp" formation
2.0
1.0 Generation
Dissipation
3
Cusp Terms /10
0.0
-1.0
-2.0
-3.0
0.0 0.2 0.4 0.6 0.8 1.0
b
Figure 3.6: \Cusp" Formation and Dissipation
and hence an important factor in determining the ame width. If the dissipation
term is removed the ame becomes unstable with the back propagating into ame
which becomes innitesimally thin.
It is dicult to make judgements about the correctness of the ame width due
to the dissparity in published results. As already mentioned the present results
are in good agreement with those of Cant and Pope [32] but are considerably
dierent to those obtained by Pope using a transported joint PDF method [3], in
which a ame, under the present conditions, would have a width of approximately
2l i.e. four time the value obtained here. Experimental data is not available for
ames under these conditions and so the issue of the ame width cannot be
immediately resolved. However, data is available for various types of burner
and stabilised ame, which suggest the ame width is of the order of the integral
length scale, so it is not expected that the present model is excessively inaccurate.
It may be possible to resolve the issue by comparison with direct simulation
results, work which is left for the future.
The one-equation model, Eqn. (2.102), does not produce a completely stable
solution in the open duct case due to the need for a prole to counter-balance
the diusion eect. The width increases continuously at a decreasing rate, al-
though the ame speed remains constant at the correct value corresponding to
the equilibrium wrinkling eq .
The \eddy break-up" model also produces a stable ame for the constant
3.2 Variable Density 42
1.0
0.8 Cusp Dissipation x2

Standard Model
0.6
b
0.4
0.2
0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
x/l
Figure 3.7: Eect of Doubling \Cusp" Dissipation
density case, with a near-symmetric prole for b (see Fig. (3.8)). However, in order
to obtain this stable condition the coecient A in Eqn. (1.4) required adjustment
until the turbulent ame speed reached the expected value and counter-balanced
the inlet unburnt gas velocity, this coecient value being 3.1. The ame is wider
than that for the eb { model, at approximately 1:2l, and the shape is quite
dierent, not exhibiting the strong asymmetry attributed to \cusps".
3.2 Variable Density

With heat release the eb { combustion model still produces a stable ame,
Fig. (3.3), and again the ame propagates very close to the expected speed with
approximately 4% error, which may be attributed to removal of the gradient
term in the eb equation (Eqn. (2.101) ! Eqn. (2.102)) as well as the diusion
error described in the last section.
Contrary to expectation, the ame is thinner with heat release than without,
0:38l as opposed to 0:5l in the previous case, the expected result being a linear
increase in ame width with density ratio [3]. This behaviour is attributed to
taking the \slip" velocity to be proportional to (Eqn. (2.96)) which causes
unphysical results at the back of the ame. Both experiment [49] and second order
closure models [26, 3] show that the \slip" velocity reaches a maximum towards
1.0
0.8
0.6
b
0.4
0.2
0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
x/l
Figure 3.8: \eddy break-up" Model b Prole
1.0
0.8
0.6
b
0.4
0.2
0.0
0.0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2
x/l
Figure 3.9: eb { Model b Prole

the back of the ame and then decreases. It is expected that this behaviour will
be corrected if the full conditionally averaged system is solved.
The \eddy break-up" fails to produce a result in this case because of the
instability induced by the coupling between the reaction rate and the velocity
divergence via the density, resulting in unphysical ame spreading, as described
in Section 1.0.1. A transient solution can however be obtained for the semi-closed
duct by tracking the ame propagation from the ignition plane to the exit, the
results of which are compared with the equivalent obtained from the eb { model
and the one-equation model. eb proles at time intervals of 0.01s are plotted in
Figs. (3.10, 3.11 and 3.13) for each of the three models. In addition, the prole
at time intervals of 0.01s is plotted in Fig. (3.12) for the eb { model.
1.0
0.8
0.6
b
0.4
0.2
0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
x/l
Figure 3.10: \eddy break-up" Model Flame Propagation
The \eddy break-up" results show the characteristic ame spreading, and it
is apparent that the the ame propagation is so slow at the back, that within the
time of calculation it does not release from the wall, and the resulting increase
in ame width produces a corresponding increase in the integral heat release and
propagation speed (see Section 1.0.1). Clearly this behaviour is unrealistic and
unacceptable. Eigenvalue analysis supports this nding in that for such condi-
tions a steady solution cannot be obtained [8]. In contrast, the eb { model
produces a smooth development of ame width (Fig. (3.11)), \wrinkling" (Fig.
(3.12)), and speed. The one-equation model produces qualitatively similar re-
sults, Fig. (3.13), the ame width developing slowly and appears to approach a
steady level, although it is not expected that a stable value would be obtained for
an arbitrarily large duct. Clearly the ame assumes a physically realistic prole,
1.0
0.8
0.6
b
0.4
0.2
0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
x/l
Figure 3.11: eb { Model Flame Propagation
8.0
7.0
6.0
5.0
Ξ
4.0
3.0
2.0
1.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
x/l
Figure 3.12: eb { Model Development

1.0
0.8
0.6
b
0.4
0.2
0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
x/l
Figure 3.13: 1 Equation Model Flame Propagation
Fig. (3.13), compared to the \eddy break-up" model, gebuPropb, and within
the space of 10l does not become unreasonably wide. It would thus seem that
the one-equation model may be applied with condence to small and moderately
sized combustion systems such as internal combustion engines, furnaces, turbines
etc., but may not be reasonable for unconned vapour cloud explosions.
In order to study the approach to equilibrium of the eb { and one-equation
models, and the speed instability of the \eddy break-up" model in more detail,
the normalised duct exit velocity vs normalised time is plotted in Fig. (3.14)
for a 50l long duct. The continuous acceleration of the \eddy break-up" model
is clear, whereas both the eb { and one-equation models approach a state of
steady propagation. The eb { model stabilises with a speed speed 30% lower
than for the open duct case, due to the dierence in ame structure, caused by
the compression of the back of the ame by the expanding burnt gas driving the
ame forward. This compression eect decreases the ame wrinkling through the
mean stretch eect term in the equation, Eqn. (2.73).
The ame speed for the one-equation model develops faster than for the two-
equation eb { model, Fig. (3.14), due to the use of the equilibrium rather
than solving for it's time development. This has an unfortunate side eect of
causing large oscillations in the velocity eld about the correct value which damp
out with time. This is a feature of the ignition process being a plane ame,
the sudden introduction of which causes pressure and hence velocity waves. Such
behaviour is not obtained for a spherically propagating ame in which the ignition
kernel is small. The ame does, however, attain the expected speed, Sueq , the
3.3 Boundary Layer Flame Propagation 47
2.0
1.5
t
Exit U/(R-1)S
1.0
0.5 Two-Equation Model

One-Equation Model
EBU Model
0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
t/τ
t
Figure 3.14: eb { and \eddy break-up" Model Exit Velocity
compression eects present in the equation not being represented through the
equilibrium assumption.
3.3 Boundary Layer Flame Propagation

The behaviour of ames in boundary layers is of great importance in many com-
bustion devices, particularly in internal combustion engines where much of the
combustion occurs near-walls. It is, however, found that many simple combustion
models are unable to reproduces, even qualitatively, near-wall ame behaviour.
It is the purpose of this section to construct a simple boundary layer test case
in which the eb { , one-equation models and \eddy break-up" models may be
studied with moderate computational eort.
Clearly, \wall functions" are inappropriate for such a study due to the many
assumptions used, and their inability to represent the important features of the
turbulence structure in the near-wall region. Alternative approaches are one-
equation models, using simple length-scale prescriptions in the near-wall region,
or low-Reynolds number models, in which damping functions are applied to the
"-equation to represent length scale changes in the near wall region. Because
neither of these schemes have been constructed with variable density systems in
mind, the use of density weighted averaging will be assumed to allow for such
eects. On this basis the low-Reynolds approach is considered more appropriate.
This choice being made, there is left the decision as to which low-Reynolds number
model is the most appropriate. The various models have been tailored to specic
ranges of ow conditions, the important features to be represented here are the
position and magnitude of the peaks in the k and " distributions near the wall.
On the basis of the review by Patel et.al. [50] the Lam-Bremhorst model [51] is
considered to represent best the important ow features of the present study.
There are various ways in which a turbulent boundary layer, through which a
ame propagates, may be generated. For example, it may be self induced as in a
ame tube experiment, or, it may be a driven ow, in which case the boundary
layer exists before the ame ignition. The problem with the former case is that a
transition to turbulence must be modelled, and that the computational domain
must be large enough for both this development period and and a period of
ame propagation through the resulting turbulent boundary layer. The latter
case results in the propagation of two opposed ames, one of which propagates
towards the gas inlet. To avoid both these sets of problems a rather articial
case is constructed, which is more computationally convenient. It consists of
two parallel plates, innitely large in one direction, Fig. (3.15). A premixed
Initial flame "Fully-developed"
Profile Symmetry Plane Profile
Ub U u
Fully developed Adiabatic Wall Constant

inlet flow Pressure
Figure 3.15: Boundary Layer Test Case
fuel-air mixture is injected at one end with a velocity and turbulence prole
corresponding to the fully-developed state. The ow is maintained until a steady
state is reached. The ame is initiated by imposing a simple ramp prole in eb
near the inlet, and setting the inlet conditions to the equivalent fully-developed
ow for the burnt gas state. The ame then propagates along the duct, the
computation ceasing when the ame reaches the outlet. The dimensions of the
duct and the ow conditions are given in Table 3.2.
The eb { model produces a ame that appears to reach a steady state of
propagation, although the duct is not long enough to be certain that such a state
is attained. The velocity vectors, eb and contours are shown in Fig. (3.16).
Property Symbol Value Dimension
Distance Between Plates h 20.0 mm
Initial mean pressure p 1.0 bar
Initial temperature T 300.0 K
Initial Velocity U Fully Developed m=s
Inlet Mean Velocity U 30.0 m=s
Initial Flame Prole
ramp Length 8.0 mm
Fuel Methane
Equivalance Ratio 0.8
Laminar Flame Speed Su 0.3 m=s
Density Ratio u=b 6.72
Laminar Viscosity 15:6 10?6 m2 =s
Table 3.2: Boundary Layer Test Conditions.
The inlet center line velocity being approximately 30m=s, the black contours of
ranging from 0:05 ! 0:95 and the coloured contours of ranging from 1 ! 6.
eb
The ame shape is clearly plausible and similar to schlieren photographs of ames
Figure 3.16: eb { Model Flame Propagation
in ducts [52]. More than this cannot be said until more detailed experiments
of ames near-walls are made. There are, however, many interesting features
in these results which require closer examination, and may shed light on the
behaviour of ames in boundary layers. It is interesting to see that the ame
\wrinkling" at the ame tip is low, indicating that the gas expansion in the
combustion process accelerates and stretches the ame tip, forcing it to become
smooth. The ame \wrinkling", with it's slower propagation rate (see Eqs. (2.64
and 2.58)), becomes dominant towards the back of the ame, and peaks near the
peak in the near-wall turbulence kinetic energy distribution. Another important
feature of the ame is that the ame width does not appear to scale with the
local turbulence integral length scale, a consequence of the propagation of ame

from the bulk towards the wall. Finally, it is interesting to note that within the
ame the contours of eb and are close to being normal to each other, although
it is not obvious why this should be the case.
The 1-equation model produces a similar result to the two-equation eb {
model in that a near fully-developed ame prole is obtained, although the eb
prole is shorter and less curved, Fig. (3.17). Again the inlet centre-line velocity
is approximately 30m=s, the black contours of eb range from 0:05 ! 0:95,[M
but the coloured contours of range from 6 ! 14. Clearly, the behaviour is
Figure 3.17: eb { 1-Equation Model Flame Propagation
profoundly dierent to the two-equation eb { model, the equilibrium assumption

creating a strong peak in the distribution corresponding to the point where k
peaks, and no ame smoothing eect is apparent at the front of the ame. In
this case it might be expected that the ame thickness would scale with the
local integral length scale, but this is not obtained because the local equilibrium
assumption is applied to the ame speed, not the thickness; which is represented
through the eb distribution. Thus ame scale information is still propagated from
the bulk towards the wall.
Considering the dierences in the distribution between the two models it
may be concluded that the local equilibrium assumption is inappropriate for
this boundary layer case, in which the transport of ame information is very
important.
The \eddy break-up" model behaves in a ridiculous manner in the region
where " reaches a maximum. Due to the coupling between the diusion term,
the combustion source term and the turbulence generation, the ame tip speed
in the near-wall region exceeds 108m=s, Fig. (3.18). In Fig. (3.18) the coloured
contours are those of eb, the range being 0:05 ! 0:95. The rapid propagation of
the ame in the near-wall region causes a shock wave to propagate away from
the wall. The ame tip speed is not limited by the speed of sound because
the transport term driving the ame tip is the diusion term, rather that the
convective transport, and thus is not limited by any physical speed constraints.
The speed is not limited either by the speed of light, due to the lack of relativistic
Figure 3.18: \Eddy break-up" Model Flame Propagation
constraints on the equation system!. Clearly, this behaviour is unphysical and

renders the \eddy break-up" model useless for combustion simulations in which
the ame approaches or propagates along a wall.
Chapter 4
Conclusion
Two new premixed combustion models have been developed by the application
of conditional averaging techniques to the turbulent ame interface. These \lam-
inar amelet" based models require conditional average properties of both the
unburnt and burnt gases. To this end, conditional averaged momentum and con-
tinuity equations have been developed with the view to solve the whole coupled
equation set. For engineering purposes it is considered unnecessary to represent
the ow eld in such detail and models are presented, relating the conditional
averaged properties to density weighted ensemble averaged values, such that the
combustion model may be incorporated into existing turbulent combustion codes.
There results two models, a two equation model involving the regress variable and
the ame \wrinkling" factor, and a one equation derivative in which the \wrinkle
factor" is represented through an analytical expression. Both models are simple
in form and do not pose any great diculty in their numerical representation.
Testing has been performed for a set of simple one-dimensional duct and
two-dimensional boundary layer cases, in which the expected ame behaviour
is known. In order to demonstrate the advantages of these models, comparison
runs for each of the test cases has been performed with a simple \eddy break-up"
model, commonly used in engineering applications.
In a open-ended one-dimensional duct it is found that two-equation model
produces a ame of constant width and speed, the width scaling with the integral
length scale of the turbulence, and the speed resulting directly from the turbulent
burning speed expression used in the model. The one-equation model is unable
to produce a ame of stable width, although the ame spreading is very slow.
The \eddy break-up" model creates a ame that spreads so rapidly that this
geometry cannot be used to assess the models capabilities.
In order to compare the eb { and one-equation models with the \eddy break-
up" model in a case for which all give a solution, test calculations have been
performed in a semi-closed duct with ignition at the closed end. In this case the
two-equation eb ? model produces a ame which approaches a steady width
and speed, the speed again corresponding closely to the turbulent burning speed
Conclusion 53
relation used in the model. As in the previous case, the one-equation model
produces a ame which is qualitatively similar to the two-equation model, except
that it's width increases continuously, but slowly. This model does, however,
reproduce the correct speed behaviour and thus is considered applicable to ame
propagation in moderate sized geometries such as internal combustion engines.
The \eddy break-up" model produces a ame which increases in speed and width
continuously and rapidly, and is therefore considered unphysical.
The dierence between the models is even more extreme in the boundary layer
calculation. Both the eb ? and one-equation models produce plausible ame
shapes and speeds, although there is a marked dierence in the ame structure
due to the local equilibrium assumption used in the one-equation model. The
\eddy break-up" model produced a ame in which a needle like protrusion of the
front propagates between the laminar sub-layer and the buer layer at phenome-
nal speeds (> 108m=s) and creates a shock wave parallel to the wall. Clearly, this
result is unphysical and demonstrates that such simple models cannot be used in
near-wall combustion.
It is thus concluded that the two models developed in this report are simple,
stable and produce encouraging results for a range of simple test cases. It remains
to apply these models to more complex combusting systems such as spark ignition
internal combustion engines and explosions in complex geometries, for which the
models are considered particularly suited.
Chapter 5
Acknowledgement
I acknowledge with gratitude, guidance and tuition freely given by Dr. C.J.
Marooney, without whose help this work would not have been possible. I would
also like to thank Prof. A.D. Gosman and Dr. R.P. Lindstedt for all the help and
encouragement they have given me throughout my work. This work has been
partly funded through a research grant given by Renault.
References
[1] Jones, W.P. and Linstedt, R.P.: Combustion and Flame, 73, 1988.
[2] Libby, P.A. and Williams, F.A.: \Fundamental Aspects", In Libby, P.A.
and Williams, F.A., editors, Turbulent Reacting Flows, Topics in Applied
Physics, volume 44 of Lecture Notes in Physics. Springer-Verlag, 1980.
[3] Anand, M.S. and Pope, S.B.: \Calculation of Premixed Turbulent Flames
by PDF Methods", Combustion and Flame, 67:127{142, 1987.
[4] Spalding, D.B.: \Mixing and Chemical Reation in Steady Conned Tur-
bulent Flames", In Thirteenth Symposium (International) on Combustion,
pages 649{657. The Combustion Institute, 1970.
[5] Spalding, D.B.: \Development of the Eddy-Break-Up Model of Turbulent
Combustion", In Sixteenth Symposium (International) on Combustion, pages
1657{1663. The Combustion Institute, 1976.
[6] Mason, H.B. and Spalding, D.B.: \Prediction of Reaction-Rates in Turbu-
lent Premixed Boundary-Layer Flows", In Combustion Institute European
Symposium, pages 601{606. Academic Press, 1973.
[7] Magnussen, B.F. and Hjertager, B.H.: \On Mathematical Modelling of Tur-
bulent Combustion with Special Emphasis on Soot Formation and Combus-
tion", In Sixteenth Symposium (International) on Combustion, pages 719{
729. The Combustion Institute, 1976.
[8] Catlin, C.A. and Lindstedt, R.P.: \Premixed Turbulent Burning Velocities
Derived from Mixing Controlled Reaction Models With Cold Front Quench-
ing", Combustion and Flame, 85:427{439, 1991.
[9] Delamare, L., Borghi, R., and Gonzalez, M.: \The Modelling and Cal-
culation of a Turbulent Premixed Flame Propagation in a Closed Vessel:
Comparison of Three Models with Experiments", Technical Report 910265,
SAE, 1991.
References 56
[10] Borghi, R., Argueyrolles, B., Gaue, S., and Souhaite, P.: \Application
of a \Presumed p.d.f." Model of Turbulent Combustion to Reciprocating
Engines", In Twenty-First Symposium (International) on Combustion, pages
1591{1599. The Combustion Institute, 1986.
[11] Said, R. and Borghi, R.: \A Simulation with a \Cellular Automaton" for
Turbulent Combustion Modelling", In Twenty-Second Symposium (Interna-
tional) on Combustion, pages 569{577. The Combustion Institute, 1988.
[12] Abraham, J., Bracco, F.V., and Reitz, R.D.: \Comparison of Computed and
Measured Premixed Charge Engine Combustion", Combustion and Flame,
60:309{322, 1985.
[13] Jennings, M.J.: \Multi-Dimensional Modelling of Turbulent Premixed
Charge Combustion", Technical Report 920589, SAE, 1992.
[14] Borghi, R. and Dutoya, D.: \On the Scales of the Fluctuations in Turbulent
Combustion", In Seventeenth Symposium (International) on Combustion,
pages 235{244. The Combustion Institute, 1978.
[15] Delamare, L.: Modelisation Numerique de la Propagation d'une Flamme
Turbulente en Milieu Conne, PhD thesis, Universite de Rouen, 1992.
[16] Mantel, T. and Borghi, R.: \Modelling
of a Premixed Turbulent Flame
in a Closed Vessel Using a \ Y0 ? "fY " Model", In Proceedings of the 3rd
g02
International Conference Innovation and Reliability in Automotive Design
and Testing, pages 701{710. ATA, 1992.
[17] Mantel, T.: \Three Dimensional Study of Flame Kernel Formation Around
a Spark Plug": SAE paper 920587.
[18] Bray, K.N.C.: \Turbulent Flows with Premixed Reactants", In P.A., Libby
and F.A., Williams, editors, Turbulent Reacting Flows, Topics in Applied
Physics, volume 44. Springer-Verlag, 1980.
[19] Peters, N.: \Laminar Flamelet Concepts in Turbulent Combustion", In
Twenty-First Symposium (International) on Combustion, pages 1231{1250.
The Combustion Institute, 1986.
[20] Bray, K.N.C.: \Methods of Including Realistic Chemical Reaction Mecha-
nisms in Turbulent Combustion Models", In J., Warnatz, editor, Complex
Chemical Reaction Systems. Springer-Verlag, 1987.
[21] Bray, K.N.C. and Moss, J.B.: \A Unied Statistical Model of the Premixed
Turbulent Flame", Acta Astronautica, 4:291{319, 1977.
References 57
[22] Libby, P.A. and Bray, K.N.C.: \Implications of the Laminar Flamelet
Model in Premixed Turbulent Combustion", Combustion and Flame, 39:33{
41, 1980.
[23] Bray, K.N.C., Libby, P.A., and Moss, J.B.: \Unied Modelling Approach
for Premixed Turbulent Combustion | Part I: General Formulation", Com-
bustion and Flame, 61:87{102, 1985.
[24] Cant, R.S. and Bray, K.N.C.: \Strained Laminar Flamelet Calculations
of Premixed Turbulent Combustion in a Closed Vessel", In Twenty-Second
Symposium (International) on Combustion, pages 791{799. The Combustion
Institute, 1988.
[25] El Tahry, S.H.: \A Turbulent Combustion Model for Premixed Charge
Engines": GM Research Report FM-133, USA 1988.
[26] Masuya, G.: \In uence of Laminar Flame Speed on Turbulent Premixed
Combustion", Combustion and Flame, 64:353{367, 1986.
[27] Kolmogorov, A., Petrovsky, I., and Piskunov, N.: section A 1(6), MG Mosaic
State University, USSR, 1937.
[28] Hakberg, B. and Gosman, A.D.: \Analytical Determination of Turbulent
Flame Speed from Combustion Models", In Twentieth Symposium (Interna-
tional) on Combustion, pages 225{232. The Combustion Institute, 1984.
[29] Marble, F.E. and Broadwell, J.E.: \The Coherent Flame Model of Chemical
Reactions", Technical Report TRW-9-PU, Project Squib Rep., 1977.
[30] Candel, S.M. and Poinsot, T.J.: \Flame Stretch and the Balance Equation
for Flame Area", Combust. Sci. and Tech, 70:1{15, 1990.
[31] Pope, S.B.: \The Evolution of Surfaces in Turbulence", Int. J. Engng Sci.,
26(5):445{469, 1988.
[32] Cant, R.S., Pope, S.B., and Bray, K.N.C.: \Modelling of Flamelet Surface-
to-Volume Ratio in Turbulent Premixed Combustion", In Twenty-Third
Symposium (International) on Combustion, pages 809{815. The Combus-
tion Institute, 1990.
[33] Weller, H.G., Marooney, C.J., and Gosman, A.D.: \A New Spectral Method
for Calculation of the Time-Varying Area of a Laminar Flame in Homoge-
neous Turbulence", In Twenty-Third Symposium (International) on Com-
bustion, pages 629{636. The Combustion Institute, 1990.
[34] Dopazo, C.: \On conditional averages for intermittent turbulent ows", J.
Fluid Mech., 81:433{438, 1977.
References 58
[35] Chen, J.Y., Lumley, J.L., and Gouldin, F.C.: \Modelling of Wrinkled Lam-
inar Flames with Intermitency and Conditional Statistics", In Twenty-First
Symposium (International) on Combustion, pages 1483{1491. The Combus-
tion Institute, 1986.
[36] Namazian, M., Talbot, L., and Robben, F.: \Density Fluctuations in Pre-
mixed Turbulent Flames", In Twentieth Symposium (International) on Com-
bustion, pages 411{419. The Combustion Institute, 1984.
[37] Namazian, M., Shephard, I.G., and Talbot, L.: \Charaterization of the
Density Fluctuations in Turbulent V-Shaped Premixed Flames", Combus-
tion and Flame, 64:299{308, 1986.
[38] Poinsot, T.J.: \Private communication", 1992.
[39] Dancy, C.L.: A Study of Second-Order Modelling of Scalar Mixing with
Scalar Intemittency, PhD thesis, Cornell University, 1983.
[40] Yeung, P.K. and Pope, S.B.: \Lagrangian Statistics from Direct Numerical
Simulation of Isotropic Turbulence", J. Fluid Mech., 207:531{586, 1989.
[41] Yeung, P.K., Girimaji, S.S., and Pope, S.B.: \Straining and Scalar Dis-
sipation on Material Surfaces in Turbulence: Implications for Flamelets",
Combustion and Flame, 79:340{365, 1990.
[42] Cant, R.S., Rutland, C.J., and Trouve, A.: \Statistics for Laminar Flamelet
Modelling", Proc. 1990 Summer Program, Ctr. Turb. Res., Stanford U.,
1990.
[43] Candel, S.M., Veynate, D., Lucas, F., Maistret, E., Darabiha, N., and
Poinsot, T.J.: \in Recent Advances in Combustion Modelling", In Lar-
routurou, B.E., editor, Series on Advances in Mathematics for Applied Sci-
ences. World Scientic, Singapore, 1990.
[44] Peters, N.: \A Spectral Closure for Premixed Turbulent Combustion in the
Flamelet Regime", J. Fluid Mech., 242:611{629, 1992.
[45] Jones, W.P.: \Models for Turbulent Flows with Variable Density and Com-
bustion", In Kollman, W., editor, Prediction Methods for Turbulent Flows,
page 379. Hemisphere Pub. Co., 1980.
[46] Issa, R.I., Ahmadi-Befrui, B., Beshay, K.R., and Gosman, A.D.: \Solution of
the Implicitly Discretised Reacting Flow Equations By Operator-Splitting",
J. Comp. Phys., 93:388{410, 1991.
References 59
[47] Herweg, R: Die Ent ammung brennbarer, turbulenter Gemishe durch elek-
trische Zundanlagen - Bildung von Flammenkernen, PhD thesis, Von der
Fakultat Elektrotechnik der Universitat Stuttgart, 1992.
[48] Santavicca, D.A., Liou, D., and North, G.L.: \A Fractal Model of Turbulent
Flame Kernel Growth", Technical Report 900024, SAE, 1990.
[49] Moss, J.B.: \Simultaneous Measurements of Concentration and Velocity in
an Open Premixed Turbulent Flame", Combustion Science and Technology,
22:119{129, 1980.
[50] Patel, V.C., Rodi, W., and Scheuerer: \Turbulence Models for Near-Wall
and Low Reynolds Number Flows: A Review", AIAA Journal, 23(9):1308{
1319, September 1985.
[51] Lam, C.K.G. and Bremhorst, K.: \A Modied Form of the k ? " Model
for Predicting Wall Turbulence", Transactions of the ASME, 103:456{460,
September 1981.
[52] Schmidt, E., Steinicke, H., and Neubert, U.: \Aufnahmen der Verbren-
nung von Gasgemischen in Rohren mit dem Eigenlicht der Flamme und bei
Schlierenbeleuchtung", VDI-Forchungsheft, 413, 1951.

Weller B Xi Combustion Model PDF

Uploaded by

Copyright:

Available Formats

Weller B Xi Combustion Model PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Weller B Xi Combustion Model PDF

Uploaded by

Copyright:

Available Formats

The Development of a New Flame Area

1.0.1 Analysis of EBU based Combustion Models

 function, in which case cg00 2

2.1 Conditional Averaging

Figure 2.1: Flame Interface

extension of that applied to intermittent turbulent ows by Dopazo [34]. With

2.1.1 Formulae for Conditional Averaging

2.1.1.2 Decomposition Operations

area per unit volume , thus

where Us is the velocity of the interface.

2.1.2 Conditional Averaging Transport Equations

2.1.3 Interface Transport Equation

Decomposing UI :n into surface averages and surface uctuation correlations :

and substituting into Eqn. (2.18) results in

2.1.4 A General De nition for Wrinkle Factor

Figure 2.2: Ensemble of Flame Interfaces

n^ per unit volume in the control volume, i.e.

in which case Eqn. (2.31) may be rewritten as

2.1.6 Transport Equation

Taking the gradient of Eqn. (2.34) gives

2.1.6.2  Transport Equation

2.1.6.3 Transport Equation Assembly

+ ^n :rDb IV:

Substitution of Eqn. (2.49) into Eqn. (2.48), and renaming term IV as D , in

+ ^n:rUt :n^ ? 1 n^ :rUs:n^ II

in which case the velocity di erence becomes

Figure 2.3: Construction of Term III for a Gaussian Flame

which ! 0 as ! 1. At the back of the ame there is no similar e ect to

2.1.7 Modelling the b Transport Equation

ble modelling of these properties in terms of density weighted ensemble averaged

where the di usion coecient DI will be modelled in terms of the conditionally

Figure 2.5: Modelling Surface Averages

 (1 ? b)uSu ? bbSb (2.61)

2.1.8 Modelling the Transport Equation

combined to create the di usion term for the transport equation

Combining Eqn. (2.63) and Eqn. (2.65) forms

D = ?r DI r ? DI r: rjrbj (2.66)

It has not proved possible to integrate analytically the spectral equations to

Collecting the terms, the nal modelled equation is

2.2 Conditionally Averaged Velocities

conditional average velocities. Although these may be obtained by solving the

2.2.1 Conditionally Averaged Continuity Equations

2.2.2 Conditionally Averaged Momentum equations

However, it is dicult to de ne unambiguously the uid velocity at the ame in-

representing surface uctuations, in which case the consistent choice of average

2.2.3 Solution of the Conditionally Averaged Equations

2.3 Ensemble Averaged Equations

combination with either the ensemble averaged or conditional averaged continu-

or correctable, the bounded approach is chosen. However, the removal of the

R_ h = ?c H (fu ? fb ) ( + (u ? ) min(; eq ))Sujrebj

2.3.3 A One Equation Model

Property Symbol Value Dimension

Table 3.1: One Dimensional Test Conditions.

3.1 Constant Density

Figure 3.3: eb { Model b Pro le

Figure 3.4: eb { Model  Pro le

Figure 3.5: eb { Model Pro le

Figure 3.6: \Cusp" Formation and Dissipation

0.8 Cusp Dissipation x2

function, in which case cg00 2

area per unit volume , thus

2.1.4 A General Denition for Wrinkle Factor

2.1.6.2 Transport Equation

in which case the velocity dierence becomes

which ! 0 as ! 1. At the back of the ame there is no similar eect to

where the diusion coecient DI will be modelled in terms of the conditionally

(1 ? b)uSu ? bbSb (2.61)

combined to create the diusion term for the transport equation

D = ?r DI r ? DI r: rjrbj (2.66)

However, it is dicult to dene unambiguously the uid velocity at the ame in-

R_ h = ?c H (fu ? fb ) ( + (u ? ) min(; eq ))Sujrebj

Figure 3.3: eb { Model b Prole

Figure 3.4: eb { Model Prole

Figure 3.5: eb { Model Prole

Figure 3.7: Eect of Doubling \Cusp" Dissipation

Figure 3.8: \eddy break-up" Model b Prole

Figure 3.9: eb { Model b Prole

profoundly dierent to the two-equation eb { model, the equilibrium assumption