- Email: [email protected]
Global stability analysis of gasless flames propagating in a cylindrical sample of energetic material: Influence of radiative heat-losses
Combustion and Flame 162 (2015) 1996–2005
Contents lists available at ScienceDirect
Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e
Global stability analysis of gasless flames propagating in a cylindrical sample of energetic material: Influence of radiative heat-losses Vadim N. Kurdyumov a,⇑, Carmen Jiménez a, Vladimir V. Gubernov b,c, Andrei V. Kolobov b a
Department of Energy, CIEMAT, Avda. Complutense 22, 28040 Madrid, Spain I.E. Tamm Theory Department, P.N. Lebedev Physical Institute of RAS, Leninskii prosp. 53, Moscow 119991, Russian Federation c ICE Lab., Far Eastern Federal University, Sukhanova str. 8, Vladivostok 690950, Russian Federation b
a r t i c l e
i n f o
Article history: Received 9 December 2014 Received in revised form 23 December 2014 Accepted 24 December 2014 Available online 24 January 2015 Keywords: Combustion waves Combustion instabilities Flame oscillations Energetic materials Gasless flame
a b s t r a c t Axisymmetric and non-axisymmetric gasless flames propagating in a narrow cylindrical sample of an energetic material are considered. The work is focused on the impact of radiative heat-losses from the outer surface of the slab on the flame structure and stability. A typical C-shaped response curve is found with two different flame velocities in the case of axisymmetric flames where the slowest velocity represents an unstable solution. Global linear stability analysis of such flames is performed. It is demonstrated that an increase in heatlosses promotes the flame instabilities for values of the Zel’dovich number lower than in the cases with the heat-losses neglected. The results of the stability analysis are successfully compared with the results of direct three-dimensional calculations showing oscillatory and spinning propagation of the gasless combustion waves. Ó 2015 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1. Introduction In 1973 in experiments with self-propagating high temperature synthesis [1] it was shown that in solid gasless combustion of cylindrical slabs the flame front can propagate in auto-oscillatory and spinning wave regimes. Flame front oscillations in such configurations had been observed as early as 1950 [2], whereas more complex dynamical regimes such as spinning waves [1] and other periodic or aperiodic combustion waves in cylindrical samples were found much later [3,4]. The nature of combustion wave pulsations was explained in [5] based on the numerical analysis of a one-dimensional model and was attributed to the diffusive-thermal instability of flame fronts. The dynamics of one-dimensional gasless flames was again investigated in [6] analytically by using large activation energy asymptotics. Cylindrical geometry was employed in [7], where the adiabatic front propagation in an infinite cylindrical slab was also studied in terms of the reaction sheet model. It was demonstrated that the basic solution, a planar flame front, may lose stability due to Hopf-type bifurcation. The stability analysis was carried out and the dispersion relation was obtained. This allowed to qualitatively explain that as the diameter of the sample was increased either ⇑ Corresponding author. Fax: +34 91 346 6269. E-mail address: [email protected] (V.N. Kurdyumov).
traveling or one-dimensional pulsating instabilities occurred, which could result in the emergence of spinning or pulsating regimes of combustion wave propagation. The inclusion of heat losses to the surroundings through the surface of the cylindrical sample makes the planar one dimensional traveling wave solution nonexistent in a cylindrical slab. As shown in both qualitative analytical and numerical investigations in [8], the steady flame front becomes curved in this case. As the heat loss is increased, the angle between the vector normal to the flame surface at the point of front attachment to the outer boundary of the slab and the cylinder axis grows, until for a certain heat transfer rate the flame front lifts off from the cylinder walls and local quenching occurs near the surface. The stability analysis for the case of conductive heat losses was carried out in [9] using an infinitely thin reaction zone model and assuming a small Biot number, which served as the asymptotically small expansion parameter. The neutral stability boundaries were found in the plane of the dimensionless activation energy versus radius for both adiabatic and nonadiabatic cases for different angular and radial wave numbers. It was demonstrated that depending on the diameter of the sample, instabilities with different cylindrical symmetries occur as the activation energy is increased. Different nonstationary and/or multidimensional regimes of flame propagation were the subject of extensive numerical analysis [10–14]. Classification of the solid-flame propagation modes was suggested in [15].
http://dx.doi.org/10.1016/j.combustflame.2014.12.018 0010-2180/Ó 2015 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1997
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
Recent advances in the design and fabrication of micro- and nanostructured composite energetic materials of shell-core type [16–22] require the development of further qualitative and quantitative understanding of the stability of flame fronts propagation in such micro-sized samples of cylindrical geometry. Due to the large surface-to-volume ratio of these micro samples, the interfacial heat transport effects become increasingly important. In view of this, in our current work we focus on the numerical analysis of the flame front propagation in cylindrical samples with radiative heat losses to the surroundings. In general terms, the previous investigations of gasless flames in cylindrical samples can be classified in two groups. The first group comprises the investigations carried out by means of the direct numerical simulations of the time-dependent equations with Arrhenius type kinetics, such as [10–14]. In this kind of studies, the accurate quantitative determination of the stability boundaries lacks details and precision. The second group of studies is based on the reaction sheet model where the reaction zone is approximated by an infinitely thin layer. When this model is used the problem becomes amenable to an analytical analysis, e.g., [6,7,9]. The main aim of the current study is to provide the global stability analysis in order to test the stability properties of the axisymmetric steady-state solutions. We consider an Arrhenius reaction model assuming that the thickness of the reaction zone is comparable to the diameter of the mixture slab, as is the case in micro-sized systems. To the best of our knowledge, such analysis has not been carried out previously. The article is arranged as follows: in Section 2, we give the general formulation; the numerical details are described in Section 3; the results describing the effect of heat-losses on the axisymmetric steady-state solutions are outlined in Section 4; the global linear stability problem is formulated in Section 5; the stability results are presented in Section 6, where the adiabatic and non-adiabatic cases are considered separately. The results of time-dependent three dimensional calculations are compared with the stability analysis in Section 7. Finally, conclusions are drawn in the last section.
where r is the Stefan–Boltzmann constant, is the emittance and k is the thermal conductivity. We assume in Eq. (1) that the environment temperature is equal to the initial temperature of the sample T0. In order to describe the flame propagation a reference frame attached to the flame is used. Following the temperature distribution along the axis of the sample, starting from the unburned side, we choose the first point z0 ¼ z0 where the temperature is equal to some value T ¼ T (the reference temperature below), see [23–25]. In the following, the reference frame is attached to this point. It is common to characterize the flame properties by the adiabatic planar burning velocity, SL , and the adiabatic flame temperature, T a ¼ T 0 þ QY 0 =c, where Y 0 is the initial fuel mass fraction and c is the heat capacity. The value of SL provides the corresponding thermal flame thickness defined as dT ¼ DT =SL , where DT ¼ k=qc is the thermal diffusivity. In the following, the combustible mass fraction is normalized with respect to Y 0 , and a non-dimensional temperature h ¼ ðT T 0 Þ=ðT a T 0 Þ is introduced. We shall use dT and R to measure the longitudinal and transverse coordinates, z ¼ z0 =dT and r ¼ r 0 =R; dT =SL and SL are used as scales for the time and velocity. Accordingly, we have 0
2. General formulation
b ¼ EðT a T 0 Þ=RT 2a , and the heat release parameter, c ¼ ðT a T 0 Þ=T a . The factor up ¼ SL =U L appears in Eq. (4) if we take the planar deflagration velocity, SL , as a scale for the velocity, where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U L ¼ DT b1 B expðE=2RT a Þ is the asymptotic value of the flame
The deflagration process in a cylindrical sample of energetic material at initial temperature T 0 is considered. In the following, z0 is the longitudinal coordinate along the sample axis, r 0 and u are the radial and azimuthal coordinates, respectively, and R denotes the radius of the sample. Primes here and thereafter denote dimensional quantities if the same notation is used for dimensional and non-dimensional variables. In all cases considered below the fresh cold side is situated at the left of the flame (at z0 ! 1), see also Fig. 3. The combustion process is modeled by an irreversible reaction of the form F ! P þ Q , where F denotes the combustible substance, P the product, and Q the heat released per unit mass of fuel. The combustion rate, X, defined as the mass of fuel consumed per unit volume and unit time, is assumed to follow an Arrhenius law of the form
X ¼ qBY expðE=RTÞ; where q; B; Y; E; R and T represent, respectively, the density, the pre-exponential factor, the mass fraction of combustible substance, the activation energy, the universal gas constant and the temperature. The heat losses considered in the present study correspond to radiation from the sample surface. They are modeled by equating the flux of energy radiating from and the conductive heat flux to the sample surface:
k
@T ¼ rðT 4 T 40 Þ 0 ; 0 @r r0 ¼R r ¼R
ð1Þ
@h @h @ 2 h 1 þ uf ¼ þ Dru h þ x; @t @z @z2 a2 @Y @Y þ uf ¼ x; @t @z
ð2Þ ð3Þ
where Dru ¼ @ 2 =@r 2 þ r 1 @=@r þ r 2 @ 2 =@ u2 and
x¼
b bðh 1Þ : Y exp 2 up 1 þ cðh 1Þ
ð4Þ
The following non-dimensional parameters appear in Eqs. (2)–(4): the non-dimensional radius of the sample measured in terms of the thermal flame thickness, a ¼ R=dT , the Zel’dovich number,
speed calculated in the limit b ! 1. The factor up ensures that for a given b and c the non-dimensional speed of a planar adiabatic flame equals one. Adequate calculation of up requires the solution of the following eigenvalue problem 2
dh=dn ¼ d h=dn2 þ x; h ¼ Y 1 ¼ 0;
dY=dn ¼ x;
n ! þ1 :
n ! 1 :
h 1 ¼ Y ¼ 0;
where x is given by Eq. (4). The numerical values of up are plotted in Fig. 1 as a function of the Zeldovich number b for various values of the heat release parameter c. We emphasize that up naturally appears in Eq. (4) when the planar burning velocity SL is used in the scale definitions. The instantaneous values of uf ðtÞ ¼ U f ðt0 Þ=SL are determined by the additional condition
hðz ; r ¼ 0; tÞ ¼ h ;
ð5Þ
where h ¼ ðT T 0 Þ=ðT a T 0 Þ is the non-dimensional reference temperature and z ¼ z0 =dT . The numerical values of the reference temperature used in the calculations were h ¼ 0:5 0:6. It should be noted that the numerical results are independent on the choice of h due to translation invariance of the problem along the direction of motion, z ! z þ const.
1998
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
3. Numerical treatment
1.08
up
γ=0.5
1.06
γ=0.6
1.04
γ=0.7
γ=0.85
1.02
1
γ=0.9
0.98
0.96
γ=0.98 0
5
10
15
20
25
30
β
35
Fig. 1. Numerical values of the factor up ¼ SL =U L appearing in Eq. (4) plotted as a function of the Zel’dovich number, b ¼ EðT a T 0 Þ=RT 2a , for various values of the heat release parameter c ¼ ðT a T 0 Þ=T a .
Eqs. (2) and (3) are to be solved subject to the following boundary conditions. The functions h and Y are 2p-periodic functions of u. The boundary condition (1) takes the following dimensionless form:
r¼1:
h i @h ¼ a b ð1 þ ch=ð1 cÞÞ4 1 ; @r
ð6Þ
where the non-dimensional heat-loss coefficient is b ¼ rT 40 dT =k ðT a T 0 Þ ¼ rT 40 =qSL QY 0 . For axisymmetric calculations, imposing @=@ u ¼ 0, the standard symmetry condition is required
r¼0:
@h=@r ¼ 0;
ð7Þ
while for three dimension calculations the lack of singularity is demanded for the temperature field at the axis, see the Appendix A. The mass conservation Eq. (3) does not contain derivatives in the radial direction r and no boundary condition is required for the fuel mass fraction at r ¼ 0 or r ¼ 1. The temperature and the fuel mass fraction take their prescribed upstream values
z ! 1 :
h ¼ 0;
Y ¼ 1:
ð8Þ
Steady as well as time-dependent computations were carried out in a finite domain, zmin < z < zmax . The typical values were zmin ¼ 20 and zmax ¼ 20 for two-dimensional calculations and zmin ¼ 10 and zmax ¼ 10 for three dimensional calculations. The spatial derivatives were discretized on a uniform grid using second order, three-point central differences for the temperature Eq. (2) and three-point upwind differences for the convection term of the mass fraction Eq. (3). The typical number of grid points was 251 51 for two dimensional calculations and 161 51 51 for three dimensional ones. The number of grid points was doubled in some cases without any significant differences in the results. For unsteady calculations an explicit marching procedure was used for the temperature equation and an implicit finite differences scheme for the mass fraction equations, both with first order discretization in time. The typical time step was s ¼ 105 . The absence of a diffusion term allowed to organize calculations of Eq. (3) also in explicit form. No significant differences were found in the results when s was halved. In order to determine steady (but not necessarily stable) solutions, the steady counterpart (@=@t ¼ 0) of Eqs. (2) and (3) were solved using a Gauss–Seidel method with over-relaxation.
4. Steady-state axisymmetric solutions The influence of volumetric heat losses on a freely propagating premixed flame has been known from the work by Spalding [26]. He found a typical C-shaped response curve with two different burning velocities where the slowest one represents an unstable solution which is difficult to observe in experiments. The same trend was found for flames with radiative heat-losses stabilized near a porous-plug burner [27]. Consider the influence of the radiative heat-losses on the steady axisymmetric gasless flames in cylindrical samples. The dimensionless velocity uf ¼ U f =SL is plotted in Fig. 2 as a function of the non-dimensional radius a ¼ R=dT calculated for various values
b=0
uf 1
b=0.0001
0.8
For the temperature outlet boundary condition we require
z!1:
2
2
@ h=@z ¼ 0:
b=0.0003 ð9Þ
The fuel mass fraction Eq. (3) contains only the first spatial derivative in z direction and no outlet boundary condition is required. This mild outlet boundary condition for the temperature field replaces a more severe zero-temperature condition which should be imposed far downstream when the heat losses are not negligible. The numerical simulations reported below showed that the influence of the downstream boundary condition becomes negligible, as it should be, if the size of the computational domain is reasonably large downstream the flame. For the non-dimensional description given above, it is convenient to separate geometric parameters of a sample, e.g. the radius in the case of cylindrical samples, which can be easily modified in experiments, from parameters dependent on the intrinsic characteristics of the energetic material. Parameters b, c and b represent the intrinsic characteristics of the energetic material and the radiating property of the external sample surface.
0.6
b=0.0005 0.4
0.2
0
0
2
4
6
8
a
10
Fig. 2. Computed steady velocity uf ¼ U f =SL as a function of the sample radius a ¼ R=dT calculated for c ¼ 0:7; b ¼ 7 (solid lines), b ¼ 6:8 (dashed lines), b ¼ 6:5 (dash-dotted lines) and several values of the heat-loss parameter b. The critical values of a ¼ amin below which propagation is not longer possible are shown with the symbol . The steady-state solutions corresponding to the points marked with the D-symbols are shown in Fig. 3.
1999
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
of the heat-loss parameter b and the Zel’dovich number b. The typical temperature and reaction rate distributions corresponding to the upper and lower branches of Fig. 2 are illustrated in Fig. 3 calculated with a ¼ 3, b ¼ 6:8; c ¼ 0:7 and b ¼ 0:0003. In both solutions the temperature decreases near the wall due to heat-losses. Anticipating the results of the stability analysis, the lower-branch solutions are unconditionally unstable and have the reaction rate significantly thicker than the upper-branch one. Notice that all solutions calculated with b ¼ 0 (no heat losses) constitute planar flame fronts. They are shown in Fig. 2 by the unique horizontal line uf ¼ 1 by virtue of up incorporated in Eq. (4). One can see that for non-zero heat-loss intensity, b > 0, the response curves take a typical C-shaped form with a critical radius, a ¼ ac , below which the propagation is not longer possible. These points are marked in Fig. 2 by open circles. It is interesting to remark that the curves calculated with different Zel’dovich numbers (solid, dashed and dash-dotted lines) are close together. The weak dependence of uf on b emerges also by virtue of the factor up included in Eq. (4).
( ) @F @ 2 F 1 @ 2 F 1 @F n2 kF ¼ uf þ 2 þ 2 þ 2 F þ AF þ BG; @z @z a @r 2 r @r r kG ¼ uf
@G AF BG; @z
ð12Þ
where
b2 Y 0
A¼
exp
u2p ½1 þ cðh0 1Þ2 b bðh0 1Þ ¼ 2 exp up 1 þ cðh0 1Þ
bðh0 1Þ ; 1 þ cðh0 1Þ
The axisymmetric steady-state solutions reported in the previous section are used in the following to examine the stability of the steady state. Two-dimensional axisymmetric distributions of the steady-state temperature and mass fraction, all now denoted by the subindex ‘‘0’’, are perturbed as usual with small perturbations
h ¼ h0 ðz; rÞ þ F ðz; rÞ expðkt þ i nuÞ;
ð10Þ
Y ¼ Y 0 ðz; rÞ þ Gðz; rÞ expðkt þ i nuÞ;
where k is a complex number the real part of which gives the growth rate, n ¼ 0; 1; 2; . . . is the azimuthal wave number and 1 is the perturbation amplitude. The mode with n ¼ 0 represents axisymmetric perturbations. As shown in [24], the corresponding infinitesimal perturbation of the flame propagation velocity can be excluded from the analysis without loss of generality. The linearized eigenvalue problem obtained when substituting Eq. (10) into Eqs. (2) and (3) reduces to finding the non-trivial solutions of the following system
ð13Þ
are both functions of x and r. The linearized temperature boundary condition at the surface becomes
r¼1:
@F ¼ CF ; @r
ð14Þ
where
is a function of z determined by the steady-state temperature h0 at the surface. Inspecting the asymptotic behavior of the temperature perturbation near the axis shows that F rn at r ! 0. Then, the following conditions should be imposed at the axis
r¼0:
@F =@r ¼ 0;
for n ¼ 0;
F ¼ 0;
for n > 0:
z ! 1 :
F ¼ G ¼ 0;
z ! þ1 :
ð16Þ
0.95 0.8 0.65 0.5 0.35 0.2 0.05
burned
unburned 0
0
@ 2 F =@z2 ¼ 0;
where again no boundary condition for the perturbed mass fraction is required downstream. The eigenvalue problem given by Eqs. (11)–(16) is not amenable to analytical analysis and numerical calculations are required. The numerical method described in [24] was applied to calculate the
a=3, β=6.8, γ=0.7, b=0.0003, uf ≈0.862 (the upper branch)
-5
ð15Þ
Taking into account that only the first z-derivative appears in Eq. (3), no boundary conditions are required for the perturbed mass fraction G neither at the axis nor at the sample surface. Far upstream and downstream we require
3
-10
B
3 4abc c C¼ 1þ h0 ð1 cÞ 1 c r¼1
5. Global linear stability analysis formulation
-3 -15
ð11Þ
5
10
15
1.3 1.1 0.9 0.7 0.5 0.3 0.1
a=3, β=6.8, γ=0.7, b=0.0003, uf ≈0.155 (the lower branch) 0.55 0.45 0.35 0.25 0.15 0.05
3
0
-3
0.05 0.04 0.03 0.02 0.01 0
-15
-10
-5
0
5
10
15
Fig. 3. Structure of axisymmetric flames with radiative heat-losses for the upper (upper plot) and lower (lower plot) branches of the C-shaped curve, calculated for b ¼ 6:8, c ¼ 0:7; b ¼ 0:0003 and a ¼ 3; upper half: temperature distribution, lower half: reaction rate distribution; the steady-states correspond to points marked with D-symbols in Fig. 2.
2000
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
eigenvalue with the greatest real part, or the main eigenvalue. This eigenvalue determines completely the stability properties of a given steady-state solution. If the real part of this eigenvalue is positive kR > 0, then the steady-state is unstable, and, conversely, if the real part is not positive, kR 6 0, then the steady-state is linearly stable.
λR β=7.25 β=7
0.1
6.1. Flame stability analysis in adiabatic cylindrical samples
F ðz; rÞ ¼ f ðzÞJ n ðl rÞ;
Gðz; rÞ ¼ gðzÞJ n ðl rÞ
dJn ðlÞ=dl ¼ 0:
ð18Þ ðnÞ j ,
In the following, we denote the roots of Eq. (18) by l where n indicates the order of the Bessel function and j is the root number. For the sake of completeness, the first values of l are given in Table 1. Notice that only for n ¼ 0 the root l ¼ 0 represent a nontrivial solution. After substituting Eq. (17) into the appropriate equations, the functions f ðzÞ and gðzÞ are determined by a non-trivial solution of the eigenvalue problem 2
kf ¼ df =dz þ d f =dz þ Af þ Bg k f ;
ð19Þ
kg ¼ dg=dz Af Bg; f ¼ g ¼ 0;
z ! þ1 :
2
β=6.954
-0.1
-0.2
β=6.5 -0.3
0
0.2
0.4
0.6
0.8
1
k
1.2
Fig. 4. Planar flame stability results: the dependence of kR on the wavenumber k calculated for c ¼ 0:7 and various values of b; for dashed segments kI ¼ 0 and for solid segments kI – 0. The critical point is shown with an open circle.
ð17Þ
where Jn ð Þ is the Bessel function of the first kind. The adiabatic boundary condition at the sample surface requires for l to satisfy the equation
2
β=6.8
0
The stability analysis of the gasless flames propagation in cylindrical thermally insulated (adiabatic) samples of a solid fuel was carried out in [7] using the reaction-sheet model. In the frame of that model the spatially distributed reaction rate given by Eq. (4) is replaced by an infinitely narrow reaction sheet modeled by a Dirac d-function. For the sake of completeness we present here the stability analysis obtained using the Arrhenius kinetics model. When b ¼ 0 is imposed in Eq. (6), the steady-state temperature and mass fraction profiles become independent on the radial and angular coordinates. Consequently, the functions A and B appearing in Eqs. (11) and (12) are also functions of z only. Following [7], eigenfunctions F and G appearing in Eqs. (11) and (12) are sought in the form
z ! 1 :
β=7.5
0.2
6. Stability results
2
γ=0.7
0.3
2
d f =dz ¼ 0;
ð20Þ
where k ¼ l=a. Notice that uf ¼ 1 is used in Eq. (19) by virtue of up introduced in Eq. (4). The eigenvalue problem given by Eqs. (19) and (20) is completely identical to the problem of determining the stability of a planar flame front in an unbounded environment where k is the transverse wavenumber. The heat release parameter c remains the only parameter in the problem. Here it was solved using the method described in [24]. The computed growth rate kR is exemplified in Fig. 4 as a function of k for different b and for c ¼ 0:7. The
solid segments represent the eigenvalues with non-zero imaginary part while for the dashed-line segment kI ¼ 0. The dashed-line segments match almost exactly for different b and the transition from a complex-conjugate eigenvalue to a purely real eigenvalue is indicated with dark circles. One can see that for values of the Zel’dovich number below a critical value, b < bc (e.g. bc 6:954 was calculated for c ¼ 0:7), the planar flame front is stable to any perturbations. With increasing values of b above bc an interval of the wavenumber appears, kmin < k < kmax , inside which the growth rate kR is positive. For ever increasing values of the Zel’dovich numbers the minimum terminal point of the unstable wavenumber interval, kmin , becomes zero at b ¼ bc (e.g. bc 7:25 was calculated for c ¼ 0:7) and for b > bc the growth rate kR is positive for 0 6 k < kmax . The critical values of the Zel’dovich number bc above which the planar flame becomes unstable are shown in Fig. 5 with a solid line as a function of c. We plot in this figure with a dashed line the imaginary part kIc in the critical point. The corresponding critical wavenumber kc indicated in Fig. 4 with an open circle is shown in the inset of the figure. Using the stability curves calculated for a freely propagating planar flame front, the stability properties of flames propagating in a thermally insulated sample of a given radius can be obtained in the following way. Consider, for example, the curve corresponding to b ¼ 7 in Fig. 4 where the growth rate kR becomes positive for ðnÞ 0:333 K k K 0:596. Using the relation a ¼ lj =k for different n and j along this curve, we plot in Fig. 6 the dependence of the growth rate kR as a function of the non-dimensional radius a for different
Table 1 The roots of Eq. (18). n
0 1 2 3 4
j 1
2
3
4
5
0 1.841183781 3.054236928 4.201188941 5.317553126
3.831705970 5.331442774 6.706133194 8.015236598 9.282396285
7.015586670 8.536316366 9.969467823 11.34592431 12.68190844
10.17346814 11.70600490 13.17037086 14.58584829 15.96410704
13.32369194 14.86358863 16.34752232 17.78874787 19.19602880
2001
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
β=6.5, γγ=0.7, b=0.0003, n=0
0.48
kc
7.5
1.68
0.47
βc
λIc
0.46
0.6
λR
0.45
0.45 0.5
0.6
0.7
0.8
0.9
7
γ
the lower branch of Fig. 2
λ
1.66 1
0.4
0.4
1.64
the upper branch of Fig. 2
a**
1.62 6.5
0.35 1.28
1.29
1.3
1.31
0.2
a 1.32
1.6
ac 6 0.5
0.6
0.7
0.8
0.9
γ
1.58 1
0
Fig. 5. Critical values of the Zel’dovich number bc (solid line), the imaginary part of the main eigenvalue kIc (dashed line) and the wave number kc (the inset of the figure) in the critical point (exemplified with an open circle for b ¼ 6:954 in Fig. 4) as a function of the heat release parameter c for a gasless planar flame front.
β=7, γ=0.7 0.02
n=2
n=1
1
15
n=1,4
2
25
a
3
Fig. 7. The growth rate kR calculated for the n ¼ 0 mode as a function of a near the extinction radius; the inset illustrates in detail the dependence of kR along the upper and lower branches of the C-shaped curve shown in Fig. 2.
β=6.5, γ=0.7, b=0.0003
0.8
n=0 n=3
a*
λI
λR
0.6
0.01
0.4
0
0.2
-0.01
ac 0
-0.02
0
5
10
a
15
Fig. 6. The growth rate kR for b ¼ 7; c ¼ 0:7 for various modes calculated for a thermally insulated slab, b ¼ 0.
n. All these curves have the same maximum value equal to that of the curve shown in Fig. 4. It is obvious that for b < bc the flame propagation in the thermally insulated slab is stable at any a and the flame shape is planar. On the other hand the flame behavior is unstable in any, even very narrow, adiabatic slabs when the Zel’dovich number exceeds some value (namely bc 7:25 for c ¼ 0:7) for which kmin is equal to zero.
6.2. Stability of flames with heat-losses It should be noted that the eigenvalue k ¼ 0 calculated for the axisymmetric mode n ¼ 0 always exists for a traveling wave solution. These perturbations correspond to translation invariance along the
1
a** 1.5
2
2.5
a
3
Fig. 8. The frequency of oscillations kI as a function of a calculated for the axisymmetric mode n ¼ 0 near the extinction radius.
direction of motion, z ! z þ const. The numerical method used in the current study, see [24] for details, yields the eigenvalue with greatest real part. Consequently, only segments with non-negative main eigenvalue can be obtained for the n ¼ 0 mode in this study. In order to gain insight into the heat-losses influence on stability, let us consider the case with the Zel’dovich number b smaller than bc shown in Fig. 5. Figure 7 shows the growth rate kR plotted versus the non-dimensional radius a for the axisymmetric mode n ¼ 0 calculated for b ¼ 6:5, c ¼ 0:7 and b ¼ 0:0003. The vertical dashed line indicates the critical radius of the slab, ac , below which the flame propagation is not longer possible. The solid segments represent the real part of k when the imaginary part is nonzero, kI – 0, and the dashed segments display the purely real eigenvalues. The inset of the figure shows the eigenvalue behavior close to the extinction value ac for two branches of the C-shaped response curve of Fig. 2. One can see that k remains real and
2002
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
β=6.5, γ=0.7, b=0.0003
0.15
λR
(a)
n=0
0.8
n=1
0
0.6
n=2
ac
n=3
-0.05
(b)
-0.1
0
200
1.2
400
600
t
800
t
800
a=7, β=6.8, γ=0.7, b=0.0003
uf 1
-0.15
-0.2
a=5, β=6.8, γ=0.7, b=0.0003
uf 1
0.1
0.05
1.2
0.8
0
5
a
10
15
Fig. 9. The growth rate kR as a function of the dimensionless radius a ¼ R=dR calculated for various modes and b ¼ 6:5; c ¼ 0:7 and b ¼ 0:0003. The modes with n ¼ 0 and n ¼ 1 are unstable.
0.6
(c)
0
200
1.2
400
600
a=9.5, β=6.8, γ=0.7, b=0.0003
uf 1
β=6.8, γ=0.7,b=0.0003
0.2
λR
n=0
0.15
0.8
n=1
0.6 0
n=2 (a) n=3
0.1
n=4
(b)
0.05
(c)
n=5
-0.05
n=0
n=7 n=6
-0.15 -0.2
0
5
10
15
400
600
800
t
1000
Fig. 11. Three time histories of uf calculated as a velocity of a point along the slab axis with a fixed temperature, see Eq. (5). The cases correspond to the points (a), (b) and (c) indicated with open circles in Fig. 10.
0
-0.1
200
a
20
Fig. 10. The growth rate kR as a function of the dimensionless radius a ¼ R=dR calculated for various modes and b ¼ 6:8; c ¼ 0:7 and b ¼ 0:0003. The modes with n P 7 are stable.
positive along the lower branch of the dual solution, which is unconditionally unstable. The curve for kR corresponding to the upper branch indicates that the axisymmetric solution losses stability to the n ¼ 0 perturbation mode for values of the radius a below some critical value a marked in Fig. 7 with an open circle. Nevertheless, the main eigenvalue has a nonzero imaginary part inside the interval a < a < a , suggesting that axisymmetric flame oscillations occur. On further decreasing a the main eigenvalue k becomes real and positive inside the interval ac < a < a indicating that the upper branch solution becomes also unconditionally unstable. Figure 8 shows the dependence of the imaginary part of the main eigenvalue for the n ¼ 0 mode near the extinction radius ac . The growth rate kR calculated along the upper branch for b ¼ 6:5; c ¼ 0:7 and various wavenumbers n (including the n ¼ 0
mode) is shown in Fig. 9. From this figure we notice that a small interval appears for the n ¼ 1 mode inside which the growth rate becomes slightly positive indicating that unstable flames exist for this value of b for slabs with radius within this interval. All modes with n > 1 remain stable for this value of the Zel’dovich number. For a somewhat higher value of the Zel’dovich number, b ¼ 6:8, still below the critical Zel’dovich number equal to bc ¼ 6:954 (calculated for c ¼ 0:7), the influence of the heat-losses on flame destabilization is still more visible. Figure 10 shows the growth rate kR plotted versus a for b ¼ 6:8. One can see that the modes with 0 6 n 6 6 become unstable for this set of parameters within certains intervals of a. It is notable that a maximum radius of the slab appears, amax 17:2, above which the gasless flame turns to be stable. This could be related to the fact that for large radii the heat-losses effect is confined to a relatively narrow region near the surface while the core of the slab, where the flame is nearly planar, remains unaffected. It should be noted that the linear stability results are reported here for the first time for both adiabatic and nonadiabatic models with distributed Arrhenius reaction rates. These results agree qualitatively with the predictions of the reaction sheet model [7,9]. 7. Unsteady flame dynamics Evidently, the results of the linear stability analysis can not predict with certainty the structure of an unstable flame. Nevertheless one can expect that if for a given set of parameters the mode with a certain n is the most unstable, namely, the growth rate of this mode exceeds the growth rates of the rest of modes, an unstable flame behavior corresponding to this symmetry mode should develop.
2003
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
a=5, β=6.8, γ=0.7, b=0.0003
a=5, β=6.8, γ=0.7, b=0.0003
θ*=0.94
z*=1.35 5
5
0
0
-10 -5
-5
z 0
5 10
0
-5
-10 -5
-5
5
z 0
5 10
5
0
-5
Fig. 12. Three-dimensional instantaneous iso-surface h ¼ 0:94 (left plot) and the temperature contour plot in the z ¼ 1:35 slice (right plot) illustrating a one-headed spinning regime. The case corresponds to a point indicated with the open circle (a) in Fig. 10 where the mode with n ¼ 1 is the most unstable.
a=7, β=6.8, γ=0.7,b=0.0003
a=7, β=6.8, γ=0.7,b=0.0003
θ*=0.98
z*=1
5
5
0
0
-5
-10 -5
z
-5
-10 -5
0 5 10
-5
0
z
5
0
5 10
5
0
-5
Fig. 13. Three-dimensional instantaneous iso-surface h ¼ 0:98 (left plot) and the temperature contour plot in the z ¼ 1 slice (right plot) illustrating a two-headed spinning regime. The case corresponds to a point indicated with the open circle (b) in Fig. 10 where the mode with n ¼ 2 is the most unstable.
a=9.5, β=6.8, γ=0.7, b=0.0003
a=9.5, β=6.8, γ=0.7, b=0.0003
θ*=0.98
z*=1.4 10
10
5
5
0
-15
0
-5 -10
-5
z
0
5
10
15
-10
-5
0
5
-10 10
-5
-15
-10
-5
z
0
5
10
15 -10
-5
0
5
-10 10
Fig. 14. Three-dimensional instantaneous iso-surface h ¼ 0:98 (left plot) and the temperature contour plot in the z ¼ 1:4 slice (right plot) illustrating a three-headed spinning regime. The case corresponds to a point indicated with the open circle (c) in Fig. 10 where the mode with n ¼ 3 is the most unstable.
In order to check this conjecture three-dimensional unsteady calculations were carried out for b ¼ 6:8; c ¼ 0:7, b ¼ 0:0003 and the slab radii marked with open circles (a), (b) and (c) in Fig. 10. The stability curves given in Fig. 10 indicate that for these operating points the modes with n ¼ 1; 2 and 3, respectively, are the most unstable.
The flame velocity uf is plotted in Fig. 11 as a function of time for these three cases. One can see that after a relatively long transition time the function uf , the velocity of the central part of the flame front, approaches a constant value in all cases. It should be noted that uf defined by Eq. (5) is the velocity of a point along the axis where the temperature takes a given value. Figures 12–14 present
2004
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
the flame shapes for the three cases marked (a), (b) and (c) in Fig. 10. The figures demonstrate that, despite the constancy of uf , the flames adopt finally multi-headed structures with a spinning behavior in the clockwise or anti clockwise direction depending on initial conditions. The spinning flame dynamics corresponding to these cases are illustrated in the Supplementary material. 8. Conclusion In this work a problem of propagation of solid gassless combustion waves in a cylindrical sample of infinite length is numerically investigated, based on a model with a single step Arrhenius reaction, for both thermally insulated and radiative heat losses boundary conditions on the cylinder wall. In the adiabatic limit, the steady-state combustion front is a planar traveling wave which exists for all parameter values. It is demonstrated that in the nonadiabatic case the flame speed is a C-shaped function of the cylinder radius i.e. there exists a minimal quenching radius of the cylindrical sample below which traveling combustion wave solutions do not exist. It is remarkable that in the case of cylindrical sample the flame speed at the quenching conditions is smaller than the critical flame speed for the case of planar combustion waves propagating in infinite media, which is usually considered to be e1=2 times slower than the adiabatic flame speed [26]. This is due to the flame front deformation mechanism which is qualitatively described in [8]. The linear stability problem for the traveling combustion waves in the infinite solid fuel cylinder under both adiabatic and nonadiabatic boundary conditions is solved numerically in three dimensional geometry, for the first time without any simplifying assumptions (like asymptotically narrow reaction zone or asymptotically small Biot number conditions). This is achieved by the use of a method proposed in [24], which allows to directly determine the mode of the linear stability problem with the largest rate of exponential growth rate as well as the corresponding eigenvalue. In the adiabatic case we have found the neutral stability boundary as a curve in the Zeldovich number vs thermal expansion parameter plane, which qualitatively agrees with the results of the one-dimensional analysis in [5], where this curve is phenomenologically fitted by a line. In the case of cylindrical geometry the traveling instability, however, occurs earlier in the parameter space than the one-dimensional pulsating instability studied numerically in [5]. In terms of the cylinder radius the following conclusions are drawn: if the sample radius is smaller or comparable to the flame thickness, only planar one-dimensional instabilities can emerge, as was the case for example in the micro and nanosized cylindrical samples studied in [22]. As the radius is increased unstable modes with angular number n ¼ 1 (single-head spinning mode), n ¼ 2 (two-head spinning mode), n ¼ 0 (drum mode), n ¼ 3 (three-head spinning mode) and so on emerge sequentially. Depending on the sample radius unstable modes with any arbitrary angular number can occur. The situation is different for the nonadiabatic case. For the range of radii smaller or comparable to the flame thickness the combustion wave does not exist due to the fold bifurcation responsible for flame quenching. The lower solution branch always exhibits the uniform instability. As the radius is increased above the critical value for the upper solution branch, unstable modes can emerge with cylindrical angular numbers in the natural order: n ¼ 0; n ¼ 1; n ¼ 2 and so on. This is different from the results described in [9], where a reaction sheet model is used and asymptotically small heat losses are assumed. The order of occurrence of unstable modes with various angular numbers in [9] follows the sequence of appearance for the adiabatic model. Thus with the increase of heat losses the instability characteristics for the traveling flame front are modified. In contrast to the adiabatic limit,
when the combustion wave stability is lost, not all cylindrical modes become unstable, but only modes with n smaller than some critical angular number and only for a certain range of radius of the cylinder. Once the traveling wave becomes unstable, nonsteady spinning modes of flame propagation with different cylindrical symmetries emerge. The detailed map of instability emergence obtained here allows us to compare the type of dominating instability and the symmetry of the spinning regime which appears as the traveling combustion wave loses stability. We have demonstrated that in the cases when there is a dominant unstable mode of the linear stability problem with angular number n, that is a mode with the largest growth rate, an exactly n-headed spinning regime emerges once the traveling wave solution becomes unstable. It is interesting to compare briefly the spinning flames observed in cases of gasless combustion of cylindrical slabs and the gaseous spinning premixed flames in narrow circular channels [28–31], see also the review [32]. It should be underlined that, in a certain sense, the dynamics of gaseous flames in narrow channels is a different problem, strongly affected by gas velocity characteristics. Moreover, in most studies, the gaseous flames are investigated either when the fresh gaseous mixture is injected from a porous head-wall [28] or the flame is located in a sudden expansion tube [29] and in a divergent channel [30] (but not in a freely propagating configuration as considered in [31]). In all cases reported in the literature, to the best of our knowledge, only one-headed spinning gaseous flames were reported which corresponds to excitation of the n ¼ 1 mode. Apparently, the difference between gasless and gaseous cases can be attributed to the infinite Lewis number of the gasless flames. Acknowledgments VNK and CJ acknowledge the support of Spanish MEC under Projects #ENE2011-27686-C02-01 and #CSD2011-00011. VVG would like to acknowledge the financial support from the Ministry of Education and Science of Russian Federation (Project 14.Y26.31.0003).
Appendix A. Numerical treatment of the axial singularity in the cylindrical coordinates The numerical treatment of the axial singularity appearing in Eq. (2), namely in the term Dru h, is based on the fact that there is not any physical singularity at the axis and the temperature field is an analytic function of coordinates everywhere. This procedure was used in cases of the three dimensional calculations. Using the conventional cylindrical coordinates ðx; yÞ ¼ ðr cos u; r sin uÞ, the lack of singularity at the axis allows the Taylor series expansion of the temperature field, hðx; y; zÞ, as follows
hðr cos u; r sin u;zÞ ¼ h0 ðzÞ þ fhx cos u þ hy sin ugr 1 1 2 hxx cos2 u þ hxy sin u cos u þ hyy sin u r 2 þ 2 2 1 1 3 2 þ hxxx cos u þ hxxy cos u sin u 6 2 1 1 2 3 þ hxyy cos u sin u þ hyyy sin u r3 2 6 þ Cðz; uÞr4 þ Oðr 6 Þ;
ðA1Þ
valid for r 1, where h0 ðzÞ represents the temperature profile along the axis. All x- and y-derivatives appearing in the righthand side of Eq. (A1) are evaluated at the axis. The factor Cðz; uÞ
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
containing the fourth derivatives is not written out in full for the sake of brevity. Averaging of Eq. (A1) with respect to the angle u gives 2 6 4 hðr; zÞ ¼ h0 ðzÞ þ 1 ðDxy hj ð0;0;zÞ Þ r þ Cr þ Oðr Þ; 4
ðA2Þ
R 2p where f ¼ ð2pÞ1 0 f du. It is easily seen that terms containing odd powers of r become equal to zero in Eq. (A2). Evidently, the value of the Laplace operator remains independent on the type of coordinates, namely Dru h ¼ Dxy h should be everywhere including the axis. By applying Eq. (A2) to the rings of grid points situated at r ¼ hr and we have 2hr and eliminating C
Dru hjr¼0 ¼
15 h0 ðzÞ þ 16 hðhr ; zÞ hð2hr ; zÞ 2
3hr
4
þ Oðhr Þ;
ðA3Þ
where hr is the radial grid step. The approximation given by Eq. (A3) 4 is accurate to within Oðhr Þ. The averaging by u procedure was carried out using the Simpson’s rule for numerical integrations. Appendix B. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.combustflame. 2014.12.018. References [1] A.G. Merzhanov, A.K. Filonenko, I.P. Borovinskaya, Sov. Phys. Dokl. 208 (4) (1973) 892–894. [2] A.F. Belyaev, L.D. Komkova, Zhur. Fiz. Khim. Akademia Nauk SSSR 24 (11) (1950) 1302–1311. [3] Yu.M. Maksimov, A.G. Merzhanov, A.T. Pak, M.N. Kuchkin, Combust. Explos. Shock Waves 17 (1981) 393–400.
2005
[4] A.V. Dvoryankin, A.G. Strunina, Combust. Explos. Shock Waves 27 (1991) 168– 172. [5] G. Shkadinskii, I. Khaikin, A.G. Merzhanov, Combust. Explos. Shock Waves 7 (1971) 15–22. [6] B.J. Matkovsky, G. Sivashinsky, SIAM J. Appl. Math. 35 (1978) 465–478. [7] G.I. Sivashinsky, SIAM J. Appl. Math. 40 (1981) 432–438. [8] V.V. Aleksandrov, A.A. Davydenko, Yu.A. Kovalenko, N.P. Poddubnyi, Combust. Explos. Shock Waves 23 (1987) 182–191. [9] J.J. Thiart, H.J. Viljoen, N.F.J. Van Rensburg, J.E. Gatica, V. Hlavacek, Combust. Sci. Technol. 82 (1992) 185–204. [10] A. Bayliss, B.J. Matkowsky, Physica D 128 (1999) 18–40. [11] A. Bayliss, B.J. Matkowsky, A.P. Aldushin, Physica D 166 (2002) 104–130. [12] T.P. Ivleva, A.G. Merzhanov, Phys. Rev. E 64 (2001) 036218. [13] T.P. Ivleva, A.G. Merzhanov, Chaos 13 (2003) 80–86. [14] T.P. Ivleva, A.G. Merzhanov, Combust. Explos. Shock Waves 39 (2003) 300– 3008. [15] T.P. Ivleva, A.G. Merzhanov, Dokl. Phys. Chem. 391 (1-3) (2003) 171–173. [16] J.T. Abrahamson, N. Nair, M.S. Strano, Nanotechnology 19 (2008) 195701. [17] S. Shimizu, W. Choi, J.T. Abrahamson, M.S. Strano, Phys. Status Solidi B 248 (2011) 2445–2448. [18] J.T. Abrahamson, C. Song, J.H. Hu, J.M. Forman, S.G. Mahajan, N. Nair, W. Choi, E.-J. Lee, M.S. Strano, Chem. Mater. 23 (2011) 4557–4562. [19] J.T. Abrahamson, M.S. Strano, J. Phys. Chem. Lett. 1 (2010) 3514–3519. [20] W. Choi, J.T. Abrahamson, J.M. Strano, M.S. Strano, Mater. Today 13 (2010) 22– 33. [21] W. Choi, S. Hong, J.T. Abrahamson, J.-H. Han, C. Song, N. Nair, S. Baik, M.S. Strano, Nat. Mater. 9 (2010) 423–429. [22] V.V. Gubernov, V.N. Kurdyumov, A.V. Kolobov, Combust. Flame 161 (2014) 2209–2214. [23] V.N. Kurdyumov, E. Fernández-Tarrazo, Combust. Flame 128 (2002) 382–394. [24] V.N. Kurdyumov, Combust. Flame 158 (2011) 1307–1317. [25] V.N. Kurdyumov, C. Jiménez, Combust. Flame 161 (2014) 927–936. [26] D.B. Spalding, Proc. R. Soc. London Ser. A 240 (1957) 83–100. [27] V.N. Kurdyumov, M. Sánchez-Sans, Proc. Combust. Inst. 34 (2013) 989–996. [28] G. Pizza, C.E. Frouzakis, J. Mantzaras, A.G. Tomboulides, K. Boulouchos, J. Fluid Mech. 658 (2010) 463–491. [29] M.J. Kwon, B.J. Lee, S.H. Chung, Combust. Flame 105 (1996) 180–188. [30] B. Xu, Y. Lu, Proc. Combust. Int. 31 (2007) 3285–3292. [31] V.N. Kurdyumov, in: Abstracts of Ginzburg Conference on Physics, Moscow, 2012, pp. 186–187. [32] Y. Lu, K. Maruta, Prog. Energy Combust. Sci. 37 (2011) 669–715.
Contents lists available at ScienceDirect
Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e
Global stability analysis of gasless flames propagating in a cylindrical sample of energetic material: Influence of radiative heat-losses Vadim N. Kurdyumov a,⇑, Carmen Jiménez a, Vladimir V. Gubernov b,c, Andrei V. Kolobov b a
Department of Energy, CIEMAT, Avda. Complutense 22, 28040 Madrid, Spain I.E. Tamm Theory Department, P.N. Lebedev Physical Institute of RAS, Leninskii prosp. 53, Moscow 119991, Russian Federation c ICE Lab., Far Eastern Federal University, Sukhanova str. 8, Vladivostok 690950, Russian Federation b
a r t i c l e
i n f o
Article history: Received 9 December 2014 Received in revised form 23 December 2014 Accepted 24 December 2014 Available online 24 January 2015 Keywords: Combustion waves Combustion instabilities Flame oscillations Energetic materials Gasless flame
a b s t r a c t Axisymmetric and non-axisymmetric gasless flames propagating in a narrow cylindrical sample of an energetic material are considered. The work is focused on the impact of radiative heat-losses from the outer surface of the slab on the flame structure and stability. A typical C-shaped response curve is found with two different flame velocities in the case of axisymmetric flames where the slowest velocity represents an unstable solution. Global linear stability analysis of such flames is performed. It is demonstrated that an increase in heatlosses promotes the flame instabilities for values of the Zel’dovich number lower than in the cases with the heat-losses neglected. The results of the stability analysis are successfully compared with the results of direct three-dimensional calculations showing oscillatory and spinning propagation of the gasless combustion waves. Ó 2015 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1. Introduction In 1973 in experiments with self-propagating high temperature synthesis [1] it was shown that in solid gasless combustion of cylindrical slabs the flame front can propagate in auto-oscillatory and spinning wave regimes. Flame front oscillations in such configurations had been observed as early as 1950 [2], whereas more complex dynamical regimes such as spinning waves [1] and other periodic or aperiodic combustion waves in cylindrical samples were found much later [3,4]. The nature of combustion wave pulsations was explained in [5] based on the numerical analysis of a one-dimensional model and was attributed to the diffusive-thermal instability of flame fronts. The dynamics of one-dimensional gasless flames was again investigated in [6] analytically by using large activation energy asymptotics. Cylindrical geometry was employed in [7], where the adiabatic front propagation in an infinite cylindrical slab was also studied in terms of the reaction sheet model. It was demonstrated that the basic solution, a planar flame front, may lose stability due to Hopf-type bifurcation. The stability analysis was carried out and the dispersion relation was obtained. This allowed to qualitatively explain that as the diameter of the sample was increased either ⇑ Corresponding author. Fax: +34 91 346 6269. E-mail address: [email protected] (V.N. Kurdyumov).
traveling or one-dimensional pulsating instabilities occurred, which could result in the emergence of spinning or pulsating regimes of combustion wave propagation. The inclusion of heat losses to the surroundings through the surface of the cylindrical sample makes the planar one dimensional traveling wave solution nonexistent in a cylindrical slab. As shown in both qualitative analytical and numerical investigations in [8], the steady flame front becomes curved in this case. As the heat loss is increased, the angle between the vector normal to the flame surface at the point of front attachment to the outer boundary of the slab and the cylinder axis grows, until for a certain heat transfer rate the flame front lifts off from the cylinder walls and local quenching occurs near the surface. The stability analysis for the case of conductive heat losses was carried out in [9] using an infinitely thin reaction zone model and assuming a small Biot number, which served as the asymptotically small expansion parameter. The neutral stability boundaries were found in the plane of the dimensionless activation energy versus radius for both adiabatic and nonadiabatic cases for different angular and radial wave numbers. It was demonstrated that depending on the diameter of the sample, instabilities with different cylindrical symmetries occur as the activation energy is increased. Different nonstationary and/or multidimensional regimes of flame propagation were the subject of extensive numerical analysis [10–14]. Classification of the solid-flame propagation modes was suggested in [15].
http://dx.doi.org/10.1016/j.combustflame.2014.12.018 0010-2180/Ó 2015 The Combustion Institute. Published by Elsevier Inc. All rights reserved.
1997
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
Recent advances in the design and fabrication of micro- and nanostructured composite energetic materials of shell-core type [16–22] require the development of further qualitative and quantitative understanding of the stability of flame fronts propagation in such micro-sized samples of cylindrical geometry. Due to the large surface-to-volume ratio of these micro samples, the interfacial heat transport effects become increasingly important. In view of this, in our current work we focus on the numerical analysis of the flame front propagation in cylindrical samples with radiative heat losses to the surroundings. In general terms, the previous investigations of gasless flames in cylindrical samples can be classified in two groups. The first group comprises the investigations carried out by means of the direct numerical simulations of the time-dependent equations with Arrhenius type kinetics, such as [10–14]. In this kind of studies, the accurate quantitative determination of the stability boundaries lacks details and precision. The second group of studies is based on the reaction sheet model where the reaction zone is approximated by an infinitely thin layer. When this model is used the problem becomes amenable to an analytical analysis, e.g., [6,7,9]. The main aim of the current study is to provide the global stability analysis in order to test the stability properties of the axisymmetric steady-state solutions. We consider an Arrhenius reaction model assuming that the thickness of the reaction zone is comparable to the diameter of the mixture slab, as is the case in micro-sized systems. To the best of our knowledge, such analysis has not been carried out previously. The article is arranged as follows: in Section 2, we give the general formulation; the numerical details are described in Section 3; the results describing the effect of heat-losses on the axisymmetric steady-state solutions are outlined in Section 4; the global linear stability problem is formulated in Section 5; the stability results are presented in Section 6, where the adiabatic and non-adiabatic cases are considered separately. The results of time-dependent three dimensional calculations are compared with the stability analysis in Section 7. Finally, conclusions are drawn in the last section.
where r is the Stefan–Boltzmann constant, is the emittance and k is the thermal conductivity. We assume in Eq. (1) that the environment temperature is equal to the initial temperature of the sample T0. In order to describe the flame propagation a reference frame attached to the flame is used. Following the temperature distribution along the axis of the sample, starting from the unburned side, we choose the first point z0 ¼ z0 where the temperature is equal to some value T ¼ T (the reference temperature below), see [23–25]. In the following, the reference frame is attached to this point. It is common to characterize the flame properties by the adiabatic planar burning velocity, SL , and the adiabatic flame temperature, T a ¼ T 0 þ QY 0 =c, where Y 0 is the initial fuel mass fraction and c is the heat capacity. The value of SL provides the corresponding thermal flame thickness defined as dT ¼ DT =SL , where DT ¼ k=qc is the thermal diffusivity. In the following, the combustible mass fraction is normalized with respect to Y 0 , and a non-dimensional temperature h ¼ ðT T 0 Þ=ðT a T 0 Þ is introduced. We shall use dT and R to measure the longitudinal and transverse coordinates, z ¼ z0 =dT and r ¼ r 0 =R; dT =SL and SL are used as scales for the time and velocity. Accordingly, we have 0
2. General formulation
b ¼ EðT a T 0 Þ=RT 2a , and the heat release parameter, c ¼ ðT a T 0 Þ=T a . The factor up ¼ SL =U L appears in Eq. (4) if we take the planar deflagration velocity, SL , as a scale for the velocity, where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U L ¼ DT b1 B expðE=2RT a Þ is the asymptotic value of the flame
The deflagration process in a cylindrical sample of energetic material at initial temperature T 0 is considered. In the following, z0 is the longitudinal coordinate along the sample axis, r 0 and u are the radial and azimuthal coordinates, respectively, and R denotes the radius of the sample. Primes here and thereafter denote dimensional quantities if the same notation is used for dimensional and non-dimensional variables. In all cases considered below the fresh cold side is situated at the left of the flame (at z0 ! 1), see also Fig. 3. The combustion process is modeled by an irreversible reaction of the form F ! P þ Q , where F denotes the combustible substance, P the product, and Q the heat released per unit mass of fuel. The combustion rate, X, defined as the mass of fuel consumed per unit volume and unit time, is assumed to follow an Arrhenius law of the form
X ¼ qBY expðE=RTÞ; where q; B; Y; E; R and T represent, respectively, the density, the pre-exponential factor, the mass fraction of combustible substance, the activation energy, the universal gas constant and the temperature. The heat losses considered in the present study correspond to radiation from the sample surface. They are modeled by equating the flux of energy radiating from and the conductive heat flux to the sample surface:
k
@T ¼ rðT 4 T 40 Þ 0 ; 0 @r r0 ¼R r ¼R
ð1Þ
@h @h @ 2 h 1 þ uf ¼ þ Dru h þ x; @t @z @z2 a2 @Y @Y þ uf ¼ x; @t @z
ð2Þ ð3Þ
where Dru ¼ @ 2 =@r 2 þ r 1 @=@r þ r 2 @ 2 =@ u2 and
x¼
b bðh 1Þ : Y exp 2 up 1 þ cðh 1Þ
ð4Þ
The following non-dimensional parameters appear in Eqs. (2)–(4): the non-dimensional radius of the sample measured in terms of the thermal flame thickness, a ¼ R=dT , the Zel’dovich number,
speed calculated in the limit b ! 1. The factor up ensures that for a given b and c the non-dimensional speed of a planar adiabatic flame equals one. Adequate calculation of up requires the solution of the following eigenvalue problem 2
dh=dn ¼ d h=dn2 þ x; h ¼ Y 1 ¼ 0;
dY=dn ¼ x;
n ! þ1 :
n ! 1 :
h 1 ¼ Y ¼ 0;
where x is given by Eq. (4). The numerical values of up are plotted in Fig. 1 as a function of the Zeldovich number b for various values of the heat release parameter c. We emphasize that up naturally appears in Eq. (4) when the planar burning velocity SL is used in the scale definitions. The instantaneous values of uf ðtÞ ¼ U f ðt0 Þ=SL are determined by the additional condition
hðz ; r ¼ 0; tÞ ¼ h ;
ð5Þ
where h ¼ ðT T 0 Þ=ðT a T 0 Þ is the non-dimensional reference temperature and z ¼ z0 =dT . The numerical values of the reference temperature used in the calculations were h ¼ 0:5 0:6. It should be noted that the numerical results are independent on the choice of h due to translation invariance of the problem along the direction of motion, z ! z þ const.
1998
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
3. Numerical treatment
1.08
up
γ=0.5
1.06
γ=0.6
1.04
γ=0.7
γ=0.85
1.02
1
γ=0.9
0.98
0.96
γ=0.98 0
5
10
15
20
25
30
β
35
Fig. 1. Numerical values of the factor up ¼ SL =U L appearing in Eq. (4) plotted as a function of the Zel’dovich number, b ¼ EðT a T 0 Þ=RT 2a , for various values of the heat release parameter c ¼ ðT a T 0 Þ=T a .
Eqs. (2) and (3) are to be solved subject to the following boundary conditions. The functions h and Y are 2p-periodic functions of u. The boundary condition (1) takes the following dimensionless form:
r¼1:
h i @h ¼ a b ð1 þ ch=ð1 cÞÞ4 1 ; @r
ð6Þ
where the non-dimensional heat-loss coefficient is b ¼ rT 40 dT =k ðT a T 0 Þ ¼ rT 40 =qSL QY 0 . For axisymmetric calculations, imposing @=@ u ¼ 0, the standard symmetry condition is required
r¼0:
@h=@r ¼ 0;
ð7Þ
while for three dimension calculations the lack of singularity is demanded for the temperature field at the axis, see the Appendix A. The mass conservation Eq. (3) does not contain derivatives in the radial direction r and no boundary condition is required for the fuel mass fraction at r ¼ 0 or r ¼ 1. The temperature and the fuel mass fraction take their prescribed upstream values
z ! 1 :
h ¼ 0;
Y ¼ 1:
ð8Þ
Steady as well as time-dependent computations were carried out in a finite domain, zmin < z < zmax . The typical values were zmin ¼ 20 and zmax ¼ 20 for two-dimensional calculations and zmin ¼ 10 and zmax ¼ 10 for three dimensional calculations. The spatial derivatives were discretized on a uniform grid using second order, three-point central differences for the temperature Eq. (2) and three-point upwind differences for the convection term of the mass fraction Eq. (3). The typical number of grid points was 251 51 for two dimensional calculations and 161 51 51 for three dimensional ones. The number of grid points was doubled in some cases without any significant differences in the results. For unsteady calculations an explicit marching procedure was used for the temperature equation and an implicit finite differences scheme for the mass fraction equations, both with first order discretization in time. The typical time step was s ¼ 105 . The absence of a diffusion term allowed to organize calculations of Eq. (3) also in explicit form. No significant differences were found in the results when s was halved. In order to determine steady (but not necessarily stable) solutions, the steady counterpart (@=@t ¼ 0) of Eqs. (2) and (3) were solved using a Gauss–Seidel method with over-relaxation.
4. Steady-state axisymmetric solutions The influence of volumetric heat losses on a freely propagating premixed flame has been known from the work by Spalding [26]. He found a typical C-shaped response curve with two different burning velocities where the slowest one represents an unstable solution which is difficult to observe in experiments. The same trend was found for flames with radiative heat-losses stabilized near a porous-plug burner [27]. Consider the influence of the radiative heat-losses on the steady axisymmetric gasless flames in cylindrical samples. The dimensionless velocity uf ¼ U f =SL is plotted in Fig. 2 as a function of the non-dimensional radius a ¼ R=dT calculated for various values
b=0
uf 1
b=0.0001
0.8
For the temperature outlet boundary condition we require
z!1:
2
2
@ h=@z ¼ 0:
b=0.0003 ð9Þ
The fuel mass fraction Eq. (3) contains only the first spatial derivative in z direction and no outlet boundary condition is required. This mild outlet boundary condition for the temperature field replaces a more severe zero-temperature condition which should be imposed far downstream when the heat losses are not negligible. The numerical simulations reported below showed that the influence of the downstream boundary condition becomes negligible, as it should be, if the size of the computational domain is reasonably large downstream the flame. For the non-dimensional description given above, it is convenient to separate geometric parameters of a sample, e.g. the radius in the case of cylindrical samples, which can be easily modified in experiments, from parameters dependent on the intrinsic characteristics of the energetic material. Parameters b, c and b represent the intrinsic characteristics of the energetic material and the radiating property of the external sample surface.
0.6
b=0.0005 0.4
0.2
0
0
2
4
6
8
a
10
Fig. 2. Computed steady velocity uf ¼ U f =SL as a function of the sample radius a ¼ R=dT calculated for c ¼ 0:7; b ¼ 7 (solid lines), b ¼ 6:8 (dashed lines), b ¼ 6:5 (dash-dotted lines) and several values of the heat-loss parameter b. The critical values of a ¼ amin below which propagation is not longer possible are shown with the symbol . The steady-state solutions corresponding to the points marked with the D-symbols are shown in Fig. 3.
1999
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
of the heat-loss parameter b and the Zel’dovich number b. The typical temperature and reaction rate distributions corresponding to the upper and lower branches of Fig. 2 are illustrated in Fig. 3 calculated with a ¼ 3, b ¼ 6:8; c ¼ 0:7 and b ¼ 0:0003. In both solutions the temperature decreases near the wall due to heat-losses. Anticipating the results of the stability analysis, the lower-branch solutions are unconditionally unstable and have the reaction rate significantly thicker than the upper-branch one. Notice that all solutions calculated with b ¼ 0 (no heat losses) constitute planar flame fronts. They are shown in Fig. 2 by the unique horizontal line uf ¼ 1 by virtue of up incorporated in Eq. (4). One can see that for non-zero heat-loss intensity, b > 0, the response curves take a typical C-shaped form with a critical radius, a ¼ ac , below which the propagation is not longer possible. These points are marked in Fig. 2 by open circles. It is interesting to remark that the curves calculated with different Zel’dovich numbers (solid, dashed and dash-dotted lines) are close together. The weak dependence of uf on b emerges also by virtue of the factor up included in Eq. (4).
( ) @F @ 2 F 1 @ 2 F 1 @F n2 kF ¼ uf þ 2 þ 2 þ 2 F þ AF þ BG; @z @z a @r 2 r @r r kG ¼ uf
@G AF BG; @z
ð12Þ
where
b2 Y 0
A¼
exp
u2p ½1 þ cðh0 1Þ2 b bðh0 1Þ ¼ 2 exp up 1 þ cðh0 1Þ
bðh0 1Þ ; 1 þ cðh0 1Þ
The axisymmetric steady-state solutions reported in the previous section are used in the following to examine the stability of the steady state. Two-dimensional axisymmetric distributions of the steady-state temperature and mass fraction, all now denoted by the subindex ‘‘0’’, are perturbed as usual with small perturbations
h ¼ h0 ðz; rÞ þ F ðz; rÞ expðkt þ i nuÞ;
ð10Þ
Y ¼ Y 0 ðz; rÞ þ Gðz; rÞ expðkt þ i nuÞ;
where k is a complex number the real part of which gives the growth rate, n ¼ 0; 1; 2; . . . is the azimuthal wave number and 1 is the perturbation amplitude. The mode with n ¼ 0 represents axisymmetric perturbations. As shown in [24], the corresponding infinitesimal perturbation of the flame propagation velocity can be excluded from the analysis without loss of generality. The linearized eigenvalue problem obtained when substituting Eq. (10) into Eqs. (2) and (3) reduces to finding the non-trivial solutions of the following system
ð13Þ
are both functions of x and r. The linearized temperature boundary condition at the surface becomes
r¼1:
@F ¼ CF ; @r
ð14Þ
where
is a function of z determined by the steady-state temperature h0 at the surface. Inspecting the asymptotic behavior of the temperature perturbation near the axis shows that F rn at r ! 0. Then, the following conditions should be imposed at the axis
r¼0:
@F =@r ¼ 0;
for n ¼ 0;
F ¼ 0;
for n > 0:
z ! 1 :
F ¼ G ¼ 0;
z ! þ1 :
ð16Þ
0.95 0.8 0.65 0.5 0.35 0.2 0.05
burned
unburned 0
0
@ 2 F =@z2 ¼ 0;
where again no boundary condition for the perturbed mass fraction is required downstream. The eigenvalue problem given by Eqs. (11)–(16) is not amenable to analytical analysis and numerical calculations are required. The numerical method described in [24] was applied to calculate the
a=3, β=6.8, γ=0.7, b=0.0003, uf ≈0.862 (the upper branch)
-5
ð15Þ
Taking into account that only the first z-derivative appears in Eq. (3), no boundary conditions are required for the perturbed mass fraction G neither at the axis nor at the sample surface. Far upstream and downstream we require
3
-10
B
3 4abc c C¼ 1þ h0 ð1 cÞ 1 c r¼1
5. Global linear stability analysis formulation
-3 -15
ð11Þ
5
10
15
1.3 1.1 0.9 0.7 0.5 0.3 0.1
a=3, β=6.8, γ=0.7, b=0.0003, uf ≈0.155 (the lower branch) 0.55 0.45 0.35 0.25 0.15 0.05
3
0
-3
0.05 0.04 0.03 0.02 0.01 0
-15
-10
-5
0
5
10
15
Fig. 3. Structure of axisymmetric flames with radiative heat-losses for the upper (upper plot) and lower (lower plot) branches of the C-shaped curve, calculated for b ¼ 6:8, c ¼ 0:7; b ¼ 0:0003 and a ¼ 3; upper half: temperature distribution, lower half: reaction rate distribution; the steady-states correspond to points marked with D-symbols in Fig. 2.
2000
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
eigenvalue with the greatest real part, or the main eigenvalue. This eigenvalue determines completely the stability properties of a given steady-state solution. If the real part of this eigenvalue is positive kR > 0, then the steady-state is unstable, and, conversely, if the real part is not positive, kR 6 0, then the steady-state is linearly stable.
λR β=7.25 β=7
0.1
6.1. Flame stability analysis in adiabatic cylindrical samples
F ðz; rÞ ¼ f ðzÞJ n ðl rÞ;
Gðz; rÞ ¼ gðzÞJ n ðl rÞ
dJn ðlÞ=dl ¼ 0:
ð18Þ ðnÞ j ,
In the following, we denote the roots of Eq. (18) by l where n indicates the order of the Bessel function and j is the root number. For the sake of completeness, the first values of l are given in Table 1. Notice that only for n ¼ 0 the root l ¼ 0 represent a nontrivial solution. After substituting Eq. (17) into the appropriate equations, the functions f ðzÞ and gðzÞ are determined by a non-trivial solution of the eigenvalue problem 2
kf ¼ df =dz þ d f =dz þ Af þ Bg k f ;
ð19Þ
kg ¼ dg=dz Af Bg; f ¼ g ¼ 0;
z ! þ1 :
2
β=6.954
-0.1
-0.2
β=6.5 -0.3
0
0.2
0.4
0.6
0.8
1
k
1.2
Fig. 4. Planar flame stability results: the dependence of kR on the wavenumber k calculated for c ¼ 0:7 and various values of b; for dashed segments kI ¼ 0 and for solid segments kI – 0. The critical point is shown with an open circle.
ð17Þ
where Jn ð Þ is the Bessel function of the first kind. The adiabatic boundary condition at the sample surface requires for l to satisfy the equation
2
β=6.8
0
The stability analysis of the gasless flames propagation in cylindrical thermally insulated (adiabatic) samples of a solid fuel was carried out in [7] using the reaction-sheet model. In the frame of that model the spatially distributed reaction rate given by Eq. (4) is replaced by an infinitely narrow reaction sheet modeled by a Dirac d-function. For the sake of completeness we present here the stability analysis obtained using the Arrhenius kinetics model. When b ¼ 0 is imposed in Eq. (6), the steady-state temperature and mass fraction profiles become independent on the radial and angular coordinates. Consequently, the functions A and B appearing in Eqs. (11) and (12) are also functions of z only. Following [7], eigenfunctions F and G appearing in Eqs. (11) and (12) are sought in the form
z ! 1 :
β=7.5
0.2
6. Stability results
2
γ=0.7
0.3
2
d f =dz ¼ 0;
ð20Þ
where k ¼ l=a. Notice that uf ¼ 1 is used in Eq. (19) by virtue of up introduced in Eq. (4). The eigenvalue problem given by Eqs. (19) and (20) is completely identical to the problem of determining the stability of a planar flame front in an unbounded environment where k is the transverse wavenumber. The heat release parameter c remains the only parameter in the problem. Here it was solved using the method described in [24]. The computed growth rate kR is exemplified in Fig. 4 as a function of k for different b and for c ¼ 0:7. The
solid segments represent the eigenvalues with non-zero imaginary part while for the dashed-line segment kI ¼ 0. The dashed-line segments match almost exactly for different b and the transition from a complex-conjugate eigenvalue to a purely real eigenvalue is indicated with dark circles. One can see that for values of the Zel’dovich number below a critical value, b < bc (e.g. bc 6:954 was calculated for c ¼ 0:7), the planar flame front is stable to any perturbations. With increasing values of b above bc an interval of the wavenumber appears, kmin < k < kmax , inside which the growth rate kR is positive. For ever increasing values of the Zel’dovich numbers the minimum terminal point of the unstable wavenumber interval, kmin , becomes zero at b ¼ bc (e.g. bc 7:25 was calculated for c ¼ 0:7) and for b > bc the growth rate kR is positive for 0 6 k < kmax . The critical values of the Zel’dovich number bc above which the planar flame becomes unstable are shown in Fig. 5 with a solid line as a function of c. We plot in this figure with a dashed line the imaginary part kIc in the critical point. The corresponding critical wavenumber kc indicated in Fig. 4 with an open circle is shown in the inset of the figure. Using the stability curves calculated for a freely propagating planar flame front, the stability properties of flames propagating in a thermally insulated sample of a given radius can be obtained in the following way. Consider, for example, the curve corresponding to b ¼ 7 in Fig. 4 where the growth rate kR becomes positive for ðnÞ 0:333 K k K 0:596. Using the relation a ¼ lj =k for different n and j along this curve, we plot in Fig. 6 the dependence of the growth rate kR as a function of the non-dimensional radius a for different
Table 1 The roots of Eq. (18). n
0 1 2 3 4
j 1
2
3
4
5
0 1.841183781 3.054236928 4.201188941 5.317553126
3.831705970 5.331442774 6.706133194 8.015236598 9.282396285
7.015586670 8.536316366 9.969467823 11.34592431 12.68190844
10.17346814 11.70600490 13.17037086 14.58584829 15.96410704
13.32369194 14.86358863 16.34752232 17.78874787 19.19602880
2001
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
β=6.5, γγ=0.7, b=0.0003, n=0
0.48
kc
7.5
1.68
0.47
βc
λIc
0.46
0.6
λR
0.45
0.45 0.5
0.6
0.7
0.8
0.9
7
γ
the lower branch of Fig. 2
λ
1.66 1
0.4
0.4
1.64
the upper branch of Fig. 2
a**
1.62 6.5
0.35 1.28
1.29
1.3
1.31
0.2
a 1.32
1.6
ac 6 0.5
0.6
0.7
0.8
0.9
γ
1.58 1
0
Fig. 5. Critical values of the Zel’dovich number bc (solid line), the imaginary part of the main eigenvalue kIc (dashed line) and the wave number kc (the inset of the figure) in the critical point (exemplified with an open circle for b ¼ 6:954 in Fig. 4) as a function of the heat release parameter c for a gasless planar flame front.
β=7, γ=0.7 0.02
n=2
n=1
1
15
n=1,4
2
25
a
3
Fig. 7. The growth rate kR calculated for the n ¼ 0 mode as a function of a near the extinction radius; the inset illustrates in detail the dependence of kR along the upper and lower branches of the C-shaped curve shown in Fig. 2.
β=6.5, γ=0.7, b=0.0003
0.8
n=0 n=3
a*
λI
λR
0.6
0.01
0.4
0
0.2
-0.01
ac 0
-0.02
0
5
10
a
15
Fig. 6. The growth rate kR for b ¼ 7; c ¼ 0:7 for various modes calculated for a thermally insulated slab, b ¼ 0.
n. All these curves have the same maximum value equal to that of the curve shown in Fig. 4. It is obvious that for b < bc the flame propagation in the thermally insulated slab is stable at any a and the flame shape is planar. On the other hand the flame behavior is unstable in any, even very narrow, adiabatic slabs when the Zel’dovich number exceeds some value (namely bc 7:25 for c ¼ 0:7) for which kmin is equal to zero.
6.2. Stability of flames with heat-losses It should be noted that the eigenvalue k ¼ 0 calculated for the axisymmetric mode n ¼ 0 always exists for a traveling wave solution. These perturbations correspond to translation invariance along the
1
a** 1.5
2
2.5
a
3
Fig. 8. The frequency of oscillations kI as a function of a calculated for the axisymmetric mode n ¼ 0 near the extinction radius.
direction of motion, z ! z þ const. The numerical method used in the current study, see [24] for details, yields the eigenvalue with greatest real part. Consequently, only segments with non-negative main eigenvalue can be obtained for the n ¼ 0 mode in this study. In order to gain insight into the heat-losses influence on stability, let us consider the case with the Zel’dovich number b smaller than bc shown in Fig. 5. Figure 7 shows the growth rate kR plotted versus the non-dimensional radius a for the axisymmetric mode n ¼ 0 calculated for b ¼ 6:5, c ¼ 0:7 and b ¼ 0:0003. The vertical dashed line indicates the critical radius of the slab, ac , below which the flame propagation is not longer possible. The solid segments represent the real part of k when the imaginary part is nonzero, kI – 0, and the dashed segments display the purely real eigenvalues. The inset of the figure shows the eigenvalue behavior close to the extinction value ac for two branches of the C-shaped response curve of Fig. 2. One can see that k remains real and
2002
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
β=6.5, γ=0.7, b=0.0003
0.15
λR
(a)
n=0
0.8
n=1
0
0.6
n=2
ac
n=3
-0.05
(b)
-0.1
0
200
1.2
400
600
t
800
t
800
a=7, β=6.8, γ=0.7, b=0.0003
uf 1
-0.15
-0.2
a=5, β=6.8, γ=0.7, b=0.0003
uf 1
0.1
0.05
1.2
0.8
0
5
a
10
15
Fig. 9. The growth rate kR as a function of the dimensionless radius a ¼ R=dR calculated for various modes and b ¼ 6:5; c ¼ 0:7 and b ¼ 0:0003. The modes with n ¼ 0 and n ¼ 1 are unstable.
0.6
(c)
0
200
1.2
400
600
a=9.5, β=6.8, γ=0.7, b=0.0003
uf 1
β=6.8, γ=0.7,b=0.0003
0.2
λR
n=0
0.15
0.8
n=1
0.6 0
n=2 (a) n=3
0.1
n=4
(b)
0.05
(c)
n=5
-0.05
n=0
n=7 n=6
-0.15 -0.2
0
5
10
15
400
600
800
t
1000
Fig. 11. Three time histories of uf calculated as a velocity of a point along the slab axis with a fixed temperature, see Eq. (5). The cases correspond to the points (a), (b) and (c) indicated with open circles in Fig. 10.
0
-0.1
200
a
20
Fig. 10. The growth rate kR as a function of the dimensionless radius a ¼ R=dR calculated for various modes and b ¼ 6:8; c ¼ 0:7 and b ¼ 0:0003. The modes with n P 7 are stable.
positive along the lower branch of the dual solution, which is unconditionally unstable. The curve for kR corresponding to the upper branch indicates that the axisymmetric solution losses stability to the n ¼ 0 perturbation mode for values of the radius a below some critical value a marked in Fig. 7 with an open circle. Nevertheless, the main eigenvalue has a nonzero imaginary part inside the interval a < a < a , suggesting that axisymmetric flame oscillations occur. On further decreasing a the main eigenvalue k becomes real and positive inside the interval ac < a < a indicating that the upper branch solution becomes also unconditionally unstable. Figure 8 shows the dependence of the imaginary part of the main eigenvalue for the n ¼ 0 mode near the extinction radius ac . The growth rate kR calculated along the upper branch for b ¼ 6:5; c ¼ 0:7 and various wavenumbers n (including the n ¼ 0
mode) is shown in Fig. 9. From this figure we notice that a small interval appears for the n ¼ 1 mode inside which the growth rate becomes slightly positive indicating that unstable flames exist for this value of b for slabs with radius within this interval. All modes with n > 1 remain stable for this value of the Zel’dovich number. For a somewhat higher value of the Zel’dovich number, b ¼ 6:8, still below the critical Zel’dovich number equal to bc ¼ 6:954 (calculated for c ¼ 0:7), the influence of the heat-losses on flame destabilization is still more visible. Figure 10 shows the growth rate kR plotted versus a for b ¼ 6:8. One can see that the modes with 0 6 n 6 6 become unstable for this set of parameters within certains intervals of a. It is notable that a maximum radius of the slab appears, amax 17:2, above which the gasless flame turns to be stable. This could be related to the fact that for large radii the heat-losses effect is confined to a relatively narrow region near the surface while the core of the slab, where the flame is nearly planar, remains unaffected. It should be noted that the linear stability results are reported here for the first time for both adiabatic and nonadiabatic models with distributed Arrhenius reaction rates. These results agree qualitatively with the predictions of the reaction sheet model [7,9]. 7. Unsteady flame dynamics Evidently, the results of the linear stability analysis can not predict with certainty the structure of an unstable flame. Nevertheless one can expect that if for a given set of parameters the mode with a certain n is the most unstable, namely, the growth rate of this mode exceeds the growth rates of the rest of modes, an unstable flame behavior corresponding to this symmetry mode should develop.
2003
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
a=5, β=6.8, γ=0.7, b=0.0003
a=5, β=6.8, γ=0.7, b=0.0003
θ*=0.94
z*=1.35 5
5
0
0
-10 -5
-5
z 0
5 10
0
-5
-10 -5
-5
5
z 0
5 10
5
0
-5
Fig. 12. Three-dimensional instantaneous iso-surface h ¼ 0:94 (left plot) and the temperature contour plot in the z ¼ 1:35 slice (right plot) illustrating a one-headed spinning regime. The case corresponds to a point indicated with the open circle (a) in Fig. 10 where the mode with n ¼ 1 is the most unstable.
a=7, β=6.8, γ=0.7,b=0.0003
a=7, β=6.8, γ=0.7,b=0.0003
θ*=0.98
z*=1
5
5
0
0
-5
-10 -5
z
-5
-10 -5
0 5 10
-5
0
z
5
0
5 10
5
0
-5
Fig. 13. Three-dimensional instantaneous iso-surface h ¼ 0:98 (left plot) and the temperature contour plot in the z ¼ 1 slice (right plot) illustrating a two-headed spinning regime. The case corresponds to a point indicated with the open circle (b) in Fig. 10 where the mode with n ¼ 2 is the most unstable.
a=9.5, β=6.8, γ=0.7, b=0.0003
a=9.5, β=6.8, γ=0.7, b=0.0003
θ*=0.98
z*=1.4 10
10
5
5
0
-15
0
-5 -10
-5
z
0
5
10
15
-10
-5
0
5
-10 10
-5
-15
-10
-5
z
0
5
10
15 -10
-5
0
5
-10 10
Fig. 14. Three-dimensional instantaneous iso-surface h ¼ 0:98 (left plot) and the temperature contour plot in the z ¼ 1:4 slice (right plot) illustrating a three-headed spinning regime. The case corresponds to a point indicated with the open circle (c) in Fig. 10 where the mode with n ¼ 3 is the most unstable.
In order to check this conjecture three-dimensional unsteady calculations were carried out for b ¼ 6:8; c ¼ 0:7, b ¼ 0:0003 and the slab radii marked with open circles (a), (b) and (c) in Fig. 10. The stability curves given in Fig. 10 indicate that for these operating points the modes with n ¼ 1; 2 and 3, respectively, are the most unstable.
The flame velocity uf is plotted in Fig. 11 as a function of time for these three cases. One can see that after a relatively long transition time the function uf , the velocity of the central part of the flame front, approaches a constant value in all cases. It should be noted that uf defined by Eq. (5) is the velocity of a point along the axis where the temperature takes a given value. Figures 12–14 present
2004
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
the flame shapes for the three cases marked (a), (b) and (c) in Fig. 10. The figures demonstrate that, despite the constancy of uf , the flames adopt finally multi-headed structures with a spinning behavior in the clockwise or anti clockwise direction depending on initial conditions. The spinning flame dynamics corresponding to these cases are illustrated in the Supplementary material. 8. Conclusion In this work a problem of propagation of solid gassless combustion waves in a cylindrical sample of infinite length is numerically investigated, based on a model with a single step Arrhenius reaction, for both thermally insulated and radiative heat losses boundary conditions on the cylinder wall. In the adiabatic limit, the steady-state combustion front is a planar traveling wave which exists for all parameter values. It is demonstrated that in the nonadiabatic case the flame speed is a C-shaped function of the cylinder radius i.e. there exists a minimal quenching radius of the cylindrical sample below which traveling combustion wave solutions do not exist. It is remarkable that in the case of cylindrical sample the flame speed at the quenching conditions is smaller than the critical flame speed for the case of planar combustion waves propagating in infinite media, which is usually considered to be e1=2 times slower than the adiabatic flame speed [26]. This is due to the flame front deformation mechanism which is qualitatively described in [8]. The linear stability problem for the traveling combustion waves in the infinite solid fuel cylinder under both adiabatic and nonadiabatic boundary conditions is solved numerically in three dimensional geometry, for the first time without any simplifying assumptions (like asymptotically narrow reaction zone or asymptotically small Biot number conditions). This is achieved by the use of a method proposed in [24], which allows to directly determine the mode of the linear stability problem with the largest rate of exponential growth rate as well as the corresponding eigenvalue. In the adiabatic case we have found the neutral stability boundary as a curve in the Zeldovich number vs thermal expansion parameter plane, which qualitatively agrees with the results of the one-dimensional analysis in [5], where this curve is phenomenologically fitted by a line. In the case of cylindrical geometry the traveling instability, however, occurs earlier in the parameter space than the one-dimensional pulsating instability studied numerically in [5]. In terms of the cylinder radius the following conclusions are drawn: if the sample radius is smaller or comparable to the flame thickness, only planar one-dimensional instabilities can emerge, as was the case for example in the micro and nanosized cylindrical samples studied in [22]. As the radius is increased unstable modes with angular number n ¼ 1 (single-head spinning mode), n ¼ 2 (two-head spinning mode), n ¼ 0 (drum mode), n ¼ 3 (three-head spinning mode) and so on emerge sequentially. Depending on the sample radius unstable modes with any arbitrary angular number can occur. The situation is different for the nonadiabatic case. For the range of radii smaller or comparable to the flame thickness the combustion wave does not exist due to the fold bifurcation responsible for flame quenching. The lower solution branch always exhibits the uniform instability. As the radius is increased above the critical value for the upper solution branch, unstable modes can emerge with cylindrical angular numbers in the natural order: n ¼ 0; n ¼ 1; n ¼ 2 and so on. This is different from the results described in [9], where a reaction sheet model is used and asymptotically small heat losses are assumed. The order of occurrence of unstable modes with various angular numbers in [9] follows the sequence of appearance for the adiabatic model. Thus with the increase of heat losses the instability characteristics for the traveling flame front are modified. In contrast to the adiabatic limit,
when the combustion wave stability is lost, not all cylindrical modes become unstable, but only modes with n smaller than some critical angular number and only for a certain range of radius of the cylinder. Once the traveling wave becomes unstable, nonsteady spinning modes of flame propagation with different cylindrical symmetries emerge. The detailed map of instability emergence obtained here allows us to compare the type of dominating instability and the symmetry of the spinning regime which appears as the traveling combustion wave loses stability. We have demonstrated that in the cases when there is a dominant unstable mode of the linear stability problem with angular number n, that is a mode with the largest growth rate, an exactly n-headed spinning regime emerges once the traveling wave solution becomes unstable. It is interesting to compare briefly the spinning flames observed in cases of gasless combustion of cylindrical slabs and the gaseous spinning premixed flames in narrow circular channels [28–31], see also the review [32]. It should be underlined that, in a certain sense, the dynamics of gaseous flames in narrow channels is a different problem, strongly affected by gas velocity characteristics. Moreover, in most studies, the gaseous flames are investigated either when the fresh gaseous mixture is injected from a porous head-wall [28] or the flame is located in a sudden expansion tube [29] and in a divergent channel [30] (but not in a freely propagating configuration as considered in [31]). In all cases reported in the literature, to the best of our knowledge, only one-headed spinning gaseous flames were reported which corresponds to excitation of the n ¼ 1 mode. Apparently, the difference between gasless and gaseous cases can be attributed to the infinite Lewis number of the gasless flames. Acknowledgments VNK and CJ acknowledge the support of Spanish MEC under Projects #ENE2011-27686-C02-01 and #CSD2011-00011. VVG would like to acknowledge the financial support from the Ministry of Education and Science of Russian Federation (Project 14.Y26.31.0003).
Appendix A. Numerical treatment of the axial singularity in the cylindrical coordinates The numerical treatment of the axial singularity appearing in Eq. (2), namely in the term Dru h, is based on the fact that there is not any physical singularity at the axis and the temperature field is an analytic function of coordinates everywhere. This procedure was used in cases of the three dimensional calculations. Using the conventional cylindrical coordinates ðx; yÞ ¼ ðr cos u; r sin uÞ, the lack of singularity at the axis allows the Taylor series expansion of the temperature field, hðx; y; zÞ, as follows
hðr cos u; r sin u;zÞ ¼ h0 ðzÞ þ fhx cos u þ hy sin ugr 1 1 2 hxx cos2 u þ hxy sin u cos u þ hyy sin u r 2 þ 2 2 1 1 3 2 þ hxxx cos u þ hxxy cos u sin u 6 2 1 1 2 3 þ hxyy cos u sin u þ hyyy sin u r3 2 6 þ Cðz; uÞr4 þ Oðr 6 Þ;
ðA1Þ
valid for r 1, where h0 ðzÞ represents the temperature profile along the axis. All x- and y-derivatives appearing in the righthand side of Eq. (A1) are evaluated at the axis. The factor Cðz; uÞ
V.N. Kurdyumov et al. / Combustion and Flame 162 (2015) 1996–2005
containing the fourth derivatives is not written out in full for the sake of brevity. Averaging of Eq. (A1) with respect to the angle u gives 2 6 4 hðr; zÞ ¼ h0 ðzÞ þ 1 ðDxy hj ð0;0;zÞ Þ r þ Cr þ Oðr Þ; 4
ðA2Þ
R 2p where f ¼ ð2pÞ1 0 f du. It is easily seen that terms containing odd powers of r become equal to zero in Eq. (A2). Evidently, the value of the Laplace operator remains independent on the type of coordinates, namely Dru h ¼ Dxy h should be everywhere including the axis. By applying Eq. (A2) to the rings of grid points situated at r ¼ hr and we have 2hr and eliminating C
Dru hjr¼0 ¼
15 h0 ðzÞ þ 16 hðhr ; zÞ hð2hr ; zÞ 2
3hr
4
þ Oðhr Þ;
ðA3Þ
where hr is the radial grid step. The approximation given by Eq. (A3) 4 is accurate to within Oðhr Þ. The averaging by u procedure was carried out using the Simpson’s rule for numerical integrations. Appendix B. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.combustflame. 2014.12.018. References [1] A.G. Merzhanov, A.K. Filonenko, I.P. Borovinskaya, Sov. Phys. Dokl. 208 (4) (1973) 892–894. [2] A.F. Belyaev, L.D. Komkova, Zhur. Fiz. Khim. Akademia Nauk SSSR 24 (11) (1950) 1302–1311. [3] Yu.M. Maksimov, A.G. Merzhanov, A.T. Pak, M.N. Kuchkin, Combust. Explos. Shock Waves 17 (1981) 393–400.
2005
[4] A.V. Dvoryankin, A.G. Strunina, Combust. Explos. Shock Waves 27 (1991) 168– 172. [5] G. Shkadinskii, I. Khaikin, A.G. Merzhanov, Combust. Explos. Shock Waves 7 (1971) 15–22. [6] B.J. Matkovsky, G. Sivashinsky, SIAM J. Appl. Math. 35 (1978) 465–478. [7] G.I. Sivashinsky, SIAM J. Appl. Math. 40 (1981) 432–438. [8] V.V. Aleksandrov, A.A. Davydenko, Yu.A. Kovalenko, N.P. Poddubnyi, Combust. Explos. Shock Waves 23 (1987) 182–191. [9] J.J. Thiart, H.J. Viljoen, N.F.J. Van Rensburg, J.E. Gatica, V. Hlavacek, Combust. Sci. Technol. 82 (1992) 185–204. [10] A. Bayliss, B.J. Matkowsky, Physica D 128 (1999) 18–40. [11] A. Bayliss, B.J. Matkowsky, A.P. Aldushin, Physica D 166 (2002) 104–130. [12] T.P. Ivleva, A.G. Merzhanov, Phys. Rev. E 64 (2001) 036218. [13] T.P. Ivleva, A.G. Merzhanov, Chaos 13 (2003) 80–86. [14] T.P. Ivleva, A.G. Merzhanov, Combust. Explos. Shock Waves 39 (2003) 300– 3008. [15] T.P. Ivleva, A.G. Merzhanov, Dokl. Phys. Chem. 391 (1-3) (2003) 171–173. [16] J.T. Abrahamson, N. Nair, M.S. Strano, Nanotechnology 19 (2008) 195701. [17] S. Shimizu, W. Choi, J.T. Abrahamson, M.S. Strano, Phys. Status Solidi B 248 (2011) 2445–2448. [18] J.T. Abrahamson, C. Song, J.H. Hu, J.M. Forman, S.G. Mahajan, N. Nair, W. Choi, E.-J. Lee, M.S. Strano, Chem. Mater. 23 (2011) 4557–4562. [19] J.T. Abrahamson, M.S. Strano, J. Phys. Chem. Lett. 1 (2010) 3514–3519. [20] W. Choi, J.T. Abrahamson, J.M. Strano, M.S. Strano, Mater. Today 13 (2010) 22– 33. [21] W. Choi, S. Hong, J.T. Abrahamson, J.-H. Han, C. Song, N. Nair, S. Baik, M.S. Strano, Nat. Mater. 9 (2010) 423–429. [22] V.V. Gubernov, V.N. Kurdyumov, A.V. Kolobov, Combust. Flame 161 (2014) 2209–2214. [23] V.N. Kurdyumov, E. Fernández-Tarrazo, Combust. Flame 128 (2002) 382–394. [24] V.N. Kurdyumov, Combust. Flame 158 (2011) 1307–1317. [25] V.N. Kurdyumov, C. Jiménez, Combust. Flame 161 (2014) 927–936. [26] D.B. Spalding, Proc. R. Soc. London Ser. A 240 (1957) 83–100. [27] V.N. Kurdyumov, M. Sánchez-Sans, Proc. Combust. Inst. 34 (2013) 989–996. [28] G. Pizza, C.E. Frouzakis, J. Mantzaras, A.G. Tomboulides, K. Boulouchos, J. Fluid Mech. 658 (2010) 463–491. [29] M.J. Kwon, B.J. Lee, S.H. Chung, Combust. Flame 105 (1996) 180–188. [30] B. Xu, Y. Lu, Proc. Combust. Int. 31 (2007) 3285–3292. [31] V.N. Kurdyumov, in: Abstracts of Ginzburg Conference on Physics, Moscow, 2012, pp. 186–187. [32] Y. Lu, K. Maruta, Prog. Energy Combust. Sci. 37 (2011) 669–715.