Combustion and Flame 150 (2007) 263–276www.elsevier.com/locate/combustﬂame
Flame acceleration in the early stages of burning in tubes
Vitaly Bychkov
a
,
∗
, V’yacheslav Akkerman
a
,
b
, Gordon Fru
a
,Arkady Petchenko
a
, LarsErik Eriksson
c
a
Institute of Physics, Umeå University, S901 87 Umeå, Sweden
b
Nuclear Safety Institute (IBRAE) of Russian Academy of Sciences, B. Tulskaya 52, 115191 Moscow, Russia
c
Department of Applied Mechanics, Chalmers University of Technology, 412 96 Göteborg, Sweden
Received 20 March 2006; received in revised form 13 December 2006; accepted 18 January 2007
Abstract
Acceleration of premixed laminar ﬂames in the early stages of burning in long tubes is considered. The acceleration mechanism was suggested earlier by Clanet and Searby [Combust. Flame 105 (1996) 225]. Accelerationhappens due to the initial ignition geometry at the tube axis when a ﬂame develops to a ﬁngershaped front, withsurface area growing exponentially in time. Flame surface area grows quite fast but only for a short time. Theanalytical theory of ﬂame acceleration is developed, which determines the growth rate, the total acceleration time,and the maximal increase of the ﬂame surface area. Direct numerical simulations of the process are performed forthe complete set of combustion equations. The simulations results and the theory are in good agreement with theprevious experiments. The numerical simulations also demonstrate ﬂame deceleration, which follows acceleration,and the socalled “tulip ﬂames.”
©
2007 Published by Elsevier Inc. on behalf of The Combustion Institute.
Keywords:
Premixed ﬂames; Flame acceleration; Tulip ﬂames; Direct numerical simulations
1. Introduction
When speaking about ﬂame acceleration, researchers typically mean detonationtodeﬂagration transition (DDT) according to the Shelkin scenario; see[1–6]. However, ﬂame may accelerate not only in thescope of the DDT problem. For example, ﬂame velocity increases due to external turbulence or intrinsic ﬂame instabilities. We do not discuss these effects in the present paper because of limited space;one can see the following reviews and papers on
*
Corresponding author. Fax: +46 90 786 66 73.
Email address:
vitaliy.bychkov@physics.umu.se(V. Bychkov).
turbulent burning [7–14] and on ﬂame instabilities[15–19].Another interesting example of ﬂame accelerationhas been suggested and studied experimentally byClanet and Searby [20]. Let us consider a ﬂame prop
agating in a cylindrical tube of radius
R
with ideallyslip adiabatic walls as shown in Fig. 1. One end of the
tube is closed; the ﬂame is ignited at the symmetryaxis at the closed end. In that case the ﬂame front develops from a hemispherical shape at the beginning toa “ﬁnger”shape; see Fig. 1. In the process of burning,
the volume of the burnt gas increases as(1)
dV dt
=
ΘS
w
U
f
,
where
U
f
is the planar ﬂame velocity,
S
w
is the total surface area of the ﬂame front, and
Θ
=
ρ
f
/ρ
b
00102180/$ – see front matter
©
2007 Published by Elsevier Inc. on behalf of The Combustion Institute.doi:10.1016/j.combustﬂame.2007.01.004
264
V. Bychkov et al. / Combustion and Flame 150 (2007) 263–276
Fig. 1. Geometry of a ﬂame accelerating due to the initialconditions.
is the expansion factor, deﬁned as the density ratioof the fuel mixture and the burnt matter. The majorcontribution to the ﬂame surface area
S
w
in Fig. 1comes from the ﬂame “skirt,”
S
w
≈
2
πRZ
tip
, where
Z
tip
is the coordinate of the ﬂame tip, and the volume of the burnt matter may be roughly evaluated as
V
≈
πR
2
Z
tip
. Then Eq. (1) reduces to(2)
dZ
tip
dt
=
2
ΘU
f
RZ
tip
,
which leads to the acceleration of the ﬂame tip as(3)
Z
tip
∝
exp
(
2
ΘU
f
t/R).
It is interesting to compare this effect to the Shelkinacceleration mechanism. The growth rate 2
ΘU
f
/R
is rather large in comparison with that obtained according to the Shelkin scenario for laminar ﬂames[5,6]. However, the Shelkin mechanism is not limitedin time; it takes place until the detonation is triggered. In contrast, the acceleration mechanism [20]works in a short time interval, and it is unlikely toproduce a DDT under normal conditions. The acceleration starts when the ﬂame evolves from a hemispherical kernel to the ﬁngershaped front of Fig. 1(
t >t
sph
); the acceleration stops when the ﬂame skirttouches the wall (
t < t
wall
). According to the experimental measurements [20], one has
t
sph
≈
0
.
1
R/U
f
,
t
wall
≈
0
.
26
R/U
f
, which leaves only a short time interval
t
wall
−
t
sph
≈
0
.
16
R/U
f
for the acceleration.It is interesting how strongly the ﬂame surface areamay increase during such a short time interval, butRef. [20] did not address this question. The purposeof the present paper is to answer this question as wellas to clarify other aspects of the ﬂame accelerationobtained experimentally by Clanet and Searby. Thecurrent study is also conceptually close to the work of Zeldovich [21].In the present paper we consider ﬂame acceleration in the early stages of burning in tubes with slipwalls, as proposed by Clanet and Searby [20]. We develop the analytical theory of ﬂame acceleration. Wederive the formulas for the acceleration time interval,
t
sph
<t <t
wall
, and for the growth rate. We show thatthe ﬂame surface area increases because of acceleration approximately by a factor of 2
Θ
in comparisonwith the tube cross section. We also perform directnumerical simulations of the ﬂame acceleration. Theanalytical and numerical results agree with the experimental data [20]. The simulations also clarify the
effect of the socalled “tulip” ﬂame.The paper is organized as follows. In Section 2we develop the analytical theory of ﬂame accelerationin the early stages of burning. Details of the directnumerical simulations are presented in Section 3. InSection 4 we compare the theory and the simulationresultstotheexperimentaldata[20].Thepaperiscon
cluded by a brief summary.
2. The analytical theory of ﬂame acceleration
We consider a ﬂame propagating in a cylindricaltube of radius
R
with ideally slip adiabatic walls andwith one end closed. The ﬂame is ignited at the symmetry axis at the closed end. In the analytical theorywe will use the dimensionless coordinates
(η
;
ξ)
=
(r
;
z)/R
, velocities
(w
;
v)
=
(u
r
;
u
z
)/U
f
, and time
τ
=
U
f
t/R
. We employ the standard model of an inﬁnitesimally thin ﬂame front, which propagates normally with the velocity
U
f
(or unity in the dimensionless variables). At the beginning, the ﬂame is ignitedat the tube axis at the closed end
(η
;
ξ)
=
(
0
;
0
)
. Initially, the front is hemispherical, but the ﬂame shapechanges as the ﬂame skirt
η
f
moves along the tubeend wall (
ξ
=
0) from the axis
η
=
0 to the side wall
η
=
1 as shown in Fig. 2. We stress that Fig. 2 and
the calculations below consider only the ﬂow inﬁnitesimally close to the wall, at
ξ
→
0. In that limit theﬂamefronttouching the wallmaybe treated aslocallycylindrical not only in the case of a ﬁnger shape, buteven for the hemispherical front. The ﬂame separatesthe ﬂow into two regions of the fresh fuel mixture andthe burntmatter. The incompressible ﬂowis described
Fig. 2. Flow close to the tube end wall.
V. Bychkov et al. / Combustion and Flame 150 (2007) 263–276
265
by the continuity equation(4)1
η∂∂η(ηw)
+
∂v∂ξ
=
0
.
The boundary condition at the endwall
ξ
=
0 is
v
=
0. We are interested in the ﬂow along the wall,in the limit
ξ
→
0. In the fuel mixture (labeled “1”)the ﬂow is potential; we assume(5)
v
1
=
A
1
ξ,
where the factor
A
1
may depend on time. The assumption is consistent with the “porous piston” effectas the ﬂame approaches the tube wall. In the hemispherical regime of ﬂame propagation it works onlysufﬁciently close to the wall as the leading term in
ξ
for
ξ
→
0. Then the radial velocity in the fuel mixtureis calculated from Eq. (4) as(6)
w
1
=
A
1
2
1
η
−
η
.
In Eq. (6) we have taken into account the boundary condition at the sidewall of the tube
w
=
0 at
η
=
1. The velocity distribution in the burnt matter(labeled “2”) takes the form(7)
v
2
=
A
2
ξ,
(8)
w
2
=−
A
2
2
η.
In general, the ﬂow in the burnt matter of Fig. 1 is rotational because of the curved ﬂame shape. However,close to the wall, the ﬂame front is locally cylindrical and the ﬂow (7) and (8) is potentialy similar to (5)
and (6). In Eq. (8) we have also taken into account the
boundary condition at the tube axis
w
=
0 at
η
=
0.To complete the solution we consider the matchingconditions at the ﬂame front
η
=
η
f
,(9)
dη
f
dτ
−
w
1
=
1
,
(10)
w
1
−
w
2
=
Θ
−
1
,
(11)
v
1
=
v
2
.
The condition (9) speciﬁes the ﬁxed propagation velocity
U
f
of the ﬂame front with respect to the fuelmixture (which is unity in scaled variables). The conditions (10) and (11) describe the jump of the normalvelocity and continuity of the tangential velocity atthe front. We stress that Eq. (11) applies only at theﬂame skirt close to the wall and it follows from theirrotational assumption combined with the cylindricalﬂame shape, i.e., no shear across the ﬂame front. Substituting Eqs. (5)–(8) into Eqs. (9)–(11), we obtain
(12)
A
1
=
A
2
=
2
(Θ
−
1
)η
f
,
and the equation for the ﬂame front(13)
dη
f
dτ
−
(Θ
−
1
)
1
−
η
2f
=
1
.
According to Eq. (13), we can separate two oppositeregimes of ﬂame propagation: when the ﬂame skirtis close to the axis
η
f
≪
1, and when it is close to thewall 1
−
η
f
≪
1.Inthelimitof
η
f
≪
1,theﬂamepropagates with the velocity
dη
f
/dτ
=
Θ
(or
˙
R
f
=
ΘU
f
).The same velocity takes place for an expanding hemispherical ﬂame front. In the other case of 1
−
η
f
≪
1, the ﬂame propagation velocity is
dη
f
/dτ
=
1 (or
˙
R
f
=
U
f
). In that limit a locally cylindrical ﬂame skirtapproaches the wall; the radial velocity of the freshfuel mixture tends to zero, and the ﬂame skirt propagates with the planar ﬂame velocity with respect tothe tube end wall. Integrating Eq. (13) with the initialcondition
η
f
=
0 at
τ
=
0, we ﬁnd(14)
τ
=
12
α
ln
Θ
+
αη
f
Θ
−
αη
f
,
or(15)
η
f
=
Θα
exp
(
2
ατ)
−
1exp
(
2
ατ)
+
1
=
Θα
tanh
(ατ),
where(16)
α
=
Θ(Θ
−
1
).
According to Eq. (15), the ﬂame velocity is equal
to the velocity of a hemispherical front close to theaxis, when 2
ατ
≪
1 and
η
f
=
Θτ
(or
R
f
=
ΘU
f
t
).Respectively, one should expect transition to the“ﬁnger”shape at the characteristic time(17)
τ
sph
≈
1
/
2
α,
when the ﬂame skirt is at
η
f
≈
0
.
46
Θ/α
. Of course,there is no exact mathematical deﬁnition of the transition time. Still, the characteristic time, Eq. (17),comes as a parameter into Eq. (15) and may be compared to the experiments [20]. For the expansion fac
tors
Θ
=
6–8, typical for propane ﬂames, we ﬁnd
α
≈
5
.
5–7
.
5 and 0
.
06
< τ
sph
<
0
.
09. The theoreticalevaluation (17) agrees quite well with the experimental estimates
τ
sph
≈
0
.
1 [20] taking into account arather vague deﬁnition of the value. Transition to the“ﬁnger” shape takes place approximately when theﬂame skirt has moved halfway to the tube side wall. Itis much easier to determine the time instant when theﬂame skirt touches the tube wall. Substituting
η
f
=
1into Eq. (14), we obtain
(18)
τ
wall
=
12
α
ln
Θ
+
αΘ
−
α
.
For the expansion factors
Θ
=
6–8 we ﬁnd 0
.
23
<τ
wall
<
0
.
28, which is in very good agreement withthe experimental evaluation
τ
wall
≈
0
.
26 [20]. Thus,
266
V. Bychkov et al. / Combustion and Flame 150 (2007) 263–276
Fig. 3. The time limits of ﬂame acceleration,
τ
sph
and
τ
wall
,versus the expansion factor
Θ
.
the ﬂame surface area grows almost exponentially intime during the interval
τ
sph
<τ <τ
wall
. The time instant
τ
sph
deﬁned by Eq. (17) and the time
τ
wall
whenthe ﬂame skirt touches the wall, Eq. (18), are shownin Fig. 3 versus the expansion factor
Θ
.We can also ﬁnd evolution of the ﬂame tip. If weconsider the ﬂow along the axis
η
=
0, as shown inFig. 4, then the equation for
ξ
tip
becomes(19)
dξ
tip
dτ
−
v
2
=
Θ.
Equation (19) is the condition of a ﬁxed propagationvelocity of a planar ﬂame front written with respect tothe burnt matter. Similarly to Fig. 2, the ﬂame shape islocally planar close to the axis, at
η
→
0. In that limitthe ﬂow may be described as a potential one, with theaxial velocity component
v
determined by a functionsimilar to Eq. (7). Besides, the solution for
v
alongthe axis has to coincide with Eq. (7) at
η
→
0,
ξ
→
0.Thus weobtain the sameformula for
v
2
along the axisas in Eqs. (7) and (12). We stress that such reasoning
does not hold away from the axis in the burnt gas,where the ﬂow is rotational. Still, in the present calculations we have to know only the gas velocity alongthe axis. As a result, we come to the differential equation for the ﬂame tip,(20)
dξ
tip
dτ
−
2
(Θ
−
1
)η
f
(τ)ξ
tip
=
Θ,
with the initial condition
ξ
tip
(
0
)
=
0 and with the solution
ξ
tip
=
Θ
4
α
exp
(
2
ατ)
−
exp
(
−
2
ατ)
(21)
=
Θ
2
α
sinh
(
2
ατ).
Just after ignition, 2
ατ
≪
1, the ﬂame tip moves inthe same way as the ﬂame skirt,
ξ
tip
=
η
f
=
Θτ
; seeEqs. (15) and (21), which are similar to an expand
ing hemispherical ﬂame front. When the ﬂame skirttouches the wall,
τ
=
τ
wall
, the position of the ﬂametip is(22)
ξ
wall
=
ξ
tip
(τ
wall
)
=
Θ
2
α
sinh
(
2
ατ
wall
)
=
Θ,
Fig. 4. Flow close to the tube axis.
or
Z
wall
=
ΘR
inthedimensionalunits.Bytheendof the acceleration we can neglect the decaying term inEq. (21). Then the ﬂame tip accelerates exponentially,
ξ
tip
∝
exp
(στ)
, with the growth rate(23)
σ
=
2
α
=
2
Θ(Θ
−
1
),
which is a little different from the model estimate 2
Θ
of Eq. (3); see [20]. This difference is small, since
α/Θ
=
0
.
91–0
.
94 for propane burning with
Θ
=
6–8.It is also interesting to evaluate total increase of the ﬂame surface area during the ﬂame acceleration.By the end of the acceleration, the ﬂame shape is almost selfsimilar, with
η
f
≈
1,
ξ
tip
≫
1, and
ξ
tip
∝
exp
(
2
ατ)
. We look for the ﬂame shape
ξ
f
=
f(η,τ)
in the form(24)
f(η,τ)
=
ϕ(η)
exp
(
2
ατ).
The accuracy of such an approximation is
Θ
−
1
≪
1,which is acceptable for the propane ﬂames used in theexperiments [20]. Then the equation of ﬂame evolu
tion may be written with respect to the burnt matter as
∂f ∂τ
+
w
2
∂f ∂η
−
v
2
=
Θ
1
+
∂f ∂η
2
1
/
2
(25)
≈−
Θ∂f ∂η.
Of course, the approximation

∂f/∂η
≫
1 holds onlynear the wall, certainly not near the axis. Still, theregion near the wall contributes most to the ﬂame surface area. Taking into account the exponential regimeof ﬂame acceleration (24) and the velocity distribution (7) and (8) we reduce Eq. (25) to
(26)2
(α
−
Θ
+
1
)ϕ
=
(Θ
−
1
)η
−
Θ
dϕdη.
Integrating Eq. (26), we ﬁnd(27)
ϕ
=
C
Θ
−
(Θ
−
1
)η
χ
,
where
C
is the constant of integration and(28)
χ
=
2
αΘ
−
1
−
1
=
2
ΘΘ
−
1
−
1
.
In the solution (27) and (28) we have to keep onlythe terms of the principal order in
Θ
−
1
≪
1, and takeinto account that
f(
0
)
=
ξ
wall
=
Θ
, Eq. (22). Then
V. Bychkov et al. / Combustion and Flame 150 (2007) 263–276
267
Eqs. (24), (27), and (28) reduce to(29)
f
=
Θ(
1
−
η)
1
/Θ
.
We stress that the solution (29) does not satisfy theboundary condition at the axis
∂f/∂η
=
0 at
η
=
0.This is, of course, a consequence of the approximation(30)
1
+
∂f ∂η
2
1
/
2
≈−
∂f ∂η
used in Eq. (25) in the limit of
Θ
−
1
≪
1. Still, theprincipal contribution to the ﬂame surface area comesfrom the ﬂame skirt. Another shortcoming of the evaluation (29) is that we adopted the potential ﬂow (7)
and (8) everywhere in the burnt matter, while in reality it holds only close to the tube end wall and tothe tube axis. Outside of these regions, the ﬂow maybe rotational because of the curved ﬂame shape. Using the evaluation (29), the scaled surface area of the
curved ﬂame front is calculated as
S
w
/πR
2
=
2
1
0
η
1
+
∂f ∂η
2
dη
(31)
≈−
2
1
0
η∂f ∂ηdη
=
2
Θ
2
Θ
+
1
.
Thus, for stoichiometric propane ﬂames with
Θ
=
8,we should expect a maximal ﬂame surface areaas large as
S
w
/πR
2
≈
14
.
2. Of course, orderofmagnitude evaluations may be obtained within aneven simpler and cruder approximation of a cylindrical shape of radius
R
and height
ΘR
with
S
w
/πR
2
≈
2
Θ
. Relative difference between this approximationand Eq. (31) is about 1
/Θ
.
3. Details of the numerical simulations
To assess accuracy of the analytical theory developed in the previous section, we performed directnumerical simulations of the hydrodynamic and combustion equations including chemical kinetics andtransport processes. In this section we are using dimensional variables. We investigated the case of anaxisymmetric ﬂow described by the equations(32)
∂ρ∂t
+
1
r∂∂r(rρu
r
)
+
∂∂z(ρu
z
)
=
0
,∂∂t (ρu
z
)
+
∂∂z
ρu
2
z
−
γ
zz
(33)
+
1
r∂∂r
r(ρu
z
u
r
−
γ
zr
)
=−
∂P ∂z,∂∂t (ρu
r
)
+
∂∂z(ρu
z
u
r
−
γ
zr
)
(34)
+
1
r∂∂r
r
ρu
2
r
−
γ
rr
=−
∂P ∂r
−
1
rγ
θθ
,∂ε∂t
+
∂∂z
(ε
+
P)u
z
−
γ
zz
u
z
−
γ
zr
u
r
+
q
z
+
1
r∂∂r
r
(ε
+
P)u
r
−
γ
rr
u
r
−
γ
zr
u
z
+
q
r
(35)
=
0
,∂∂t (ρY)
+
1
r∂∂r
ρru
r
Y
−
rµ
Sc
∂Y ∂r
+
∂∂z
ρu
z
Y
−
µ
Sc
∂Y ∂z
(36)
=−
ρY τ
R
exp
(
−
E
a
/R
p
T),
where(37)
ε
=
ρ(QY
+
C
V
T)
+
ρ
2
u
2
z
+
u
2
r
is the total energy per unit volume,
Y
is the mass fraction of the fuel mixture,
Q
is the energy release in thereaction, and
C
V
is the heat capacity per unit massat constant volume. We consider a single irreversiblereaction of the ﬁrst order. Temperature dependenceof the reaction rate obeys the Arrhenius law with theactivation energy
E
a
and the constant of time dimension
τ
R
. The stress tensor
γ
ij
is(38)
γ
zz
=
µ
43
∂u
z
∂z
−
23
∂u
r
∂r
−
23
u
r
r
,
(39)
γ
rr
=
µ
43
∂u
r
∂r
−
23
∂u
z
∂z
−
23
u
r
r
,
(40)
γ
θθ
=
µ
43
u
r
r
−
23
∂u
z
∂z
−
23
∂u
r
∂r
,
(41)
γ
zr
=
µ
∂u
z
∂r
+
∂u
r
∂z
,
and the heat diffusion vector
q
i
is given by(42)
q
z
=−
µ
C
P
Pr
∂T ∂z
+
Q
Sc
∂Y ∂z
,
(43)
q
r
=−
µ
C
P
Pr
∂T ∂r
+
Q
Sc
∂Y ∂r
.
Here
µ
=
ρν
is the dynamic viscosity,
C
P
is the heatcapacityperunitmassatconstantpressure,andPrandSc are the Prandtl and Schmidt numbers, respectively.The fuel mixture is assumed to be a perfect gas of constant molecular weight
m
=
2
.
9
×
10
−
2
kg
/
mol,with
C
P
=
7
R
p
/
2
m
and
C
V
=
5
R
p
/
2
m
, where
R
p
=
8
.
31 J
/(
molK
)
is the perfect gas constant. The equation of state is(44)
P
=
ρmR
p
T.