C O M B U S T I O N A N D F L A M E 64:321336 (1986)
321
Effect of Large Scale Structures on Turbulent Flame Propagation A H M E D F. G H O N I E M Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139
The effect of laminar burning velocity, Reynolds number, and turbulent fluctuations on the propagation of a turbulent premixed flame stabilized behind a sudden expansion is investigated numerically using a reaction sheet model and line interface calculations of the front motion. Within this framework, flame propagation as an interface between reactants and products is due to the combined action of the advection of products and burning of reactants. Results indicate that turbulent flames have a narrow range of stability when propagating inside channels At lean conditions, blowoff is caused by the low burning velocity and large turbulent fluctuations, At nearstoichiometric conditions, flashback may occur when the flame blocks the flow and propagates along the wall. Turbulent fluctuations establish an extra mass flux across the flame front that enhances burning at high Reynolds numbers and moderate burning velocities.
I. I N T R O D U C T I O N Most flames of practical interest burn under turbulent conditions to enhance the rate of heat release and to meet the loading requirements in small size combustors [1]. This study focuses on the continuous combustion of premixed gases where the fuel and oxidizer are well homogenized prior to burning, Under these circumstances, the burning process is stabilized by inducing flow separation, which maintains a continuous recirculation of products that provides the energy required to sustain the ignition of reactants. Under turbulent conditions, mixing within the wake vortex may dominate over the diffusivekinetic processes and determine whether the flame will blow off, flash back, or propagate at a steady state [2]. Thus, the anchoring of the flame is a strong function of the local burning rate and the Reynolds number, representing the physical properties of combustion and the dynamics of the flow field, respectively. Flame propagation in flow fields where the scales of turbulent eddies are large with respect to the laminar flame thickness can be modeled by thin reaction sheets which exist on the outer Copyright © 1986 by The Combustion Institute Published by Elsevier Science Publishing Co., Inc. 52 Vanderbdt Avenue, New York, NY 10017
edges of the eddies [3]. Turbulent eddies provide the stirring mechanism between reactants and products, while combustion establishes a density field across the flame structure. The process can thus be modeled in terms of the interaction between three fundamental components: thin laminar flame sheets, turbulent eddies, and volumetric expansion, as well as vorticity adjustment across the flame surfaces. The combined effect of these interactions, if properly simulated, can be used to explain a wide variety of phenomena related to the enhancement of burning due to turbulence, flame stability, the effect of confinement and geometrical obstacles on flame acceleration or extinction, and the interaction between combustion and flow instabilities. In the following, the origin of each component of this model is described and its contribution to turbulent combustion is examined. Laminar flames propagate through the diffusion of hot products carrying energy and radical species ahead of the reaction zone. The slow process of molecular diffusion yields low rates of burning; laminar burning velocities are on the order of 0.100.50 m/s. Turbulence augments burning by (1) distorting the flame surface, thus
00102180/86/$03.50
322 generating larger areas of contact betwee.~ re~. tants and products; and (2) stretching layers o! reactants between products, hence reducing their thickness to molecular scales and improving the mixing process. Consequently, the rate of burning increases sharply, and the apparent flame speeds reach an order of magnitude higher than the laminar values [4]. Since turbulent fields can be decomposed into a spectrum of eddies with a distribution of length scales, the interaction between turbulence and combustion can be analyzed in terms of the relative scale of an eddy and the thickness of a laminar flame [5]. Large scale turbulence enhances the burning process by wrinkling, folding, and stretching the flame surface, hence increasing its area while locally the flame still propagates as a laminar wave. A review of some of the experimental evidence for the existence of wrinkled laminar flames in turbulent fields is given by Chomiak [6]. Small scales of turbulence, of a size smaller than the nominal laminar flame thickness, may affect the laminar flame speed by perturbing the mixing within the flame structure. However, this process has a significant effect only when the Kolmogorov microscale is of the order of magnitude of the laminar flame thickness, a situation that is likely to occur in thick flames at low temperatures. For large scale interactions, the thin wrinkled flame model has been successfully postulated for the interpretation of data on turbulent combustion in premixed gases when the flame thickness is smaller than the characteristic size of turbulent eddies; see, e.g., Ballal and Lefebvre [7]. The formation of large scale eddy structures downstream of points of separation and at areas of high shear, and their persistence by coherent pairing, has been established experimentally in both nonreacting flow (e.g., the extensive review of Cantwell [8]) and reacting flow [9]. Accumulation of vorticity in large structures was observed in the numerical computations of nonreacting mixing layers by Ashurst [10], in confined recirculating flow with combustion by Ghoniem et al. [11], and in nonreacting and reacting jets by Ghoniem et al. [12]. In nonreacting flow, vortex structures appear as a
AHMED F. GHONIEM ult of an intrinsic flow instability, the KelvmHelmholtz instability, and emphasize the deterministic nature of turbulent flow that can be exploited in constructing predictive computational algorithms. In reacting flow of premixed gases, eddies are formed of products of combustion and represent the agent of flame stabilization. The recirculating action of eddies keeps reactants and products in contact long enough for molecular processes to propagate through the mixture. Attempts to study turbulent combustion interactions in terms of thin laminar flame sheets that exist in a field of turbulent eddies and cause volumetric expansion have been plagued by the following complications: (1) the representation of chemical reactions in a fluctuating field; (2) the lack of proper general models to compute turbulent flow; and (3) the numerical instabilities that may arise when sudden expansion across the flame front is present in the model. The crux of the problem of numerical solution of turbulent combustion lies in the presence of the Arrhenius exponential term in the kinetic rate equations that describe the chemical reaction. Attempts to incorporate turbulent fluctuations in the flow equations which model the combustion process using statistical decomposition of the thermodynamic variables result in divergent solutions. Meanwhile, large density and velocity gradients associated with the presence of thin reaction zones cause numerical instability and may require the use of numerical algorithms of high resolution [13]. Moreover, the utilization of fixed Eulerian grids to represent the gradients of flow variables imposes a limitation on the resolution and can only be improved if complicated adaptive grids around zones of sharp gradients are implemented. Numerical methods based on the discretization of functions on fixed grids introduce large values of artificial viscosity that inhibit the growth of natural flow instabilities at high Reynolds numbers and "laminarize" the flow. In this work, some of these difficulties are addressed. A numerical implementation of the thin reaction sheet model is investigated and applied to study flame propagation in the recir
TURBULENT FLAME PROPAGATION culation zone behind a rearwardfacing step. While all the dynamic effects of combustion are removed by assuming that the density ratio across the flame is close to one, the effect of the properties of the velocity field and the normal burning velocity on flame propagation and total burning rate is emphasized. The differential equation that governs the motion of a reaction sheet is formulated in Section II, and the properties of the model are discussed. It is shown that the model is not restricted to simply connected reaction fronts; hence its application can be extended to a regime where multiple reaction sheets coexist. In Section III, the numerical scheme used to compute the solution of the equation is introduced and applied to simple cases in which one mechanism of flame propagation prevails. The accuracy of the calculations is assessed using examples of flame propagation where exact solutions exist. Means of improving the accuracy and optimizing the efficiency of the computations, and the formulation of the boundary conditions associated with flame propagation along walls, are also discussed. Results of flame propagation at three values of Reynolds number in the range of laminar to early turbulent flow are presented in Section IV. In each case, two forms of the velocity field are used to advect the flame, namely, a steady and an unsteady velocity distribution; hence the effect of turbulent fluctuations at high Reynolds numbers appears only at the second case. A range of stable flame propagation, bounded by lean blowoff at high Reynolds numbers and low burning velocities and the possibility of rich flashback at high burning velocities, is established. Numerical results suggest mechanisms for both processes. The change of the overall combustion rate with the normal burning velocity is analyzed in terms of changes in the flame surface geometry and length.
323 Reynolds numberDamkohler number chart. The Damkohler number is defined as the ratio of a characteristic flowlength scale and a laminar flame thickness or, equivalently, the ratio of an average flow velocity and a normal burning velocity. In this work, attention is focused on the reaction sheet regime that arises at conditions of large Damkohler number, and, to a less degree of dependence, at high Reynolds numbers. In this regime, combustion is confined to a thin reaction zone that exists in a turbulence field dominated by large scale eddies. In order to extend the range of application of the model to lower values of the Damkohler number, both single and multiple fronts are considered and the formation of islands of reactants within the products, or vice versa, is allowed. Figure la describes schematically the conceptual model of turbulent combustion in premixed gases based on the conditions described above. Turbulence is represented by a set of eddies with a wide range of length scales that interact continuously according to the laws of conservation of vorticity. Combustion is a slow exothermic process that results in the transformation of reactants into products, accompanied by large rates of heat release and volumetric expansion. The link between the two processes is manifested by the deformation of the surface of the flame by the turbulent field and the emergence of an extra velocity field by the expansion across the flame front. Figure lb demonstrates the wrinkling of a flame surface due to a turbulent eddy that results in the stretch and enlargement of the exposed area and the possibility of formation of islands. Ghoniem, Chorin and Oppenheim [15] postulated the following eikonal equation to describe the propagation of a thin flame in a turbulent field: dz =u+S,n, dt
(1)
II. P H Y S I C A L M O D E L
Different regimes of turbulent combustion in premixed gases have been identified in the studies of Bray [3] and Williams [14], on a
where z = (x, y) is the coordinate for a point on the flame front, u = (u, o) is the flow velocity at z, S,, is the normal burning velocity at constant pressure from the reactants side, n is the normal
324
AHMED F. GHONIEM
TURBULENCE~.
A set of interactingeddies of o wide range of length scales
accelerohon 1 pressure rise vortlclty
Expansion
Advectlon
mlx,ng ,~ curvoture
f
L
~ , ~
stretch
COMBUSTi0 N~ / / ~ /
An exothermlc conversion of reactants into products Slow, selfsustained
(a)
Turbulent Combustion Interochons
R
P (b)
The wrinkling and stretch of flame surface due to the action of turbulent eddies
Fig. 1. Schematic diagram representing the fundamental concepts of thin flame modelling of turbulent combustion m premlxed gases: (a) definitions and modes of interaction between the two basic processes; (b) the wrmkhng of a flame surface caused by a turbulent eddy. to the flame front at z, and t is time. If the flame is moving in a quiescent medium then, u = 0 and the front propagates as a wave described by the line z. In Eq. (1), it is assumed that the flame remains infinitely thin compared with the characteristic dimension of the flow, i.e., that the chemistry is infinitely fast so that the reaction zone is a sharp interface between reactants and products. The combined effect of chemical reaction and molecular diffusion within the flame structure is lumped into the laminar burning velocity S,, which is supplied to the model from previous calculations or from experimental measurements. S, is, in general, a function of the local dynamic and thermodynamic properties of the field. The coupling between the flow field and the dynamics of the burning process is implicit in u and n; the first is employed to move the flame front and the second determines the direction of the laminar flame propagation. The composition of the reacting mixture consists of two components, reactants R and products P. Their interaction is governed by the following irreversible
reaction: R ~P.
(2)
Thus, Su is the transformation rate of reactants into products per unit length of the flame front and, within the framework of this model, the propagation of the flame front is due to relative advection of products with respect to the reactants by the velocity field and the burning of reactants at a rate determined by the laminar burning velocity and the direction normal to the flame front. Therefore, the motion of the flame front depends on the properties of the mixture and the dynamics of the flow field. Furthermore, it is assumed that both components flow as incompressible gases with constant properties and that acoustic interactions are negligible, i.e., the Mach number is low and the possibility of transition to detonation is barred. While the pressure can change temporally as a result of the expansion of gases across the flame front, the flame propagates as a slow deflagration without a noticeable change in pressure.
TURBULENT FLAME PROPAGATION Majda and Sethian [16] constructed a rigorous proof to show that Eq. (1) represents the formal limit of the full set of the reactionconvectiondiffusion equations that describe the conservation of species in a chemically reacting system, provided that the Mach number is low and that the reaction zone thickness is small compared with a characteristic flow length. These limits are consistent with the assumptions used in this section. While the application of this equation is restricted to what is known as the thin wrinkled flame model in turbulent combustion of premixed gases, adapting this model under the proper circumstances helps to remove some of the ambiguity associated with this regime of turbulent combustion and allows the solution of the governing equations using welldeveloped numerical techniques. A similar formulation was suggested by Chorin [17] to study the flamevortex interaction with the help of the line interface motion algorithm which will be described in the next section. Sethian [18] employed the model of Ghoniem et al. [15] to study the effect of confinement on the propagation of a flame in a region of intense vorticity. Hsiao et al. [19] showed that the model is capable of accurately predicting the average velocity profile for the nonreacting and the reacting flow behind a rearwardfacing step. III. T H E C O M P U T A T I O N A L ALGORITHM Equation (1) is integrated numerically in two fractional steps. In the first step, the flame front is transported by the advection of products with respect to reactants using a velocity field that is evaluated from the solution of the fluid flow equations. In the second step, the flame front motion is caused by the burning of reactants and is accompanied by the formation of products. The order of execution of the two steps, advection and burning, should be reversed every time step to avoid a possible bias toward a particular process and improve the accuracy of the fractional splitting scheme. If the equation is linear in terms of the operators that are being sepa
325 rated, this "symmetric splitting scheme" results in an improvement in the accuracy by an order of magnitude and renders the algorithm secondorder in time. In the two steps, a line interface calculation algorithm that was developed by Noh and Woodward [20], and further refined for higher resolution by Chorin [17], is used to construct the flame surface from a description of the concentration field. The latter is subsequently transported by the appropriate velocity field. In the following, the numerical algorithm is summarized.
III.1. Advection of Products The change in the concentration of products due to advection is governed by ac
a
a
  +   (uc)   (vc)=O, at ax +oy
(3)
where c is the mass of products per unit volume or, for an incompressible flow, the volume of product per unit volume. Equation (3) is an Eulerian representation of the first term of the righthand side of Eq. (I). The divergence form of the equation is used to satisfy the necessary condition for the construction of a conservative scheme. The solution of Eq. (3) is constructed according to a finite volume scheme [21] by dividing the domain of interest into a number of contiguous cells of dimensions (Ax, Ay), and integrating each term over the boundary of the cell. Let f be the volume fraction of products in a cell, where 0 _< f __< 1 ; the result of the integration is d a t (f+o a x ]Xy)+ { u f AY} +,x'+,/z '/2
+ { v f AX} yjY~+l/aO 1/2 v,
(4)
where the two limits are taken at the cell boundaries, and their difference represents the
326
AHMED F. GHONIEM
net flux of products into the cell (for the definition of symbols, see Fig. 2). Equation (4) expresses the conservation of volumetric concentration of products in terms of the time change inside a cell due to the advection flux across its boundaries. A staggered mesh is utilized to define the value o f f and (u, v) at the center of the cell and across its boundaries, respectively, thus resembling a central differencing approximation. A numerical scheme based on Eq. (4) is conservative; hence it is capable of preserving sharp discontinuities and can be applied across material interfaces without causing excessive artificial diffusion. The algorithm used here differs from standard difference methods by the introduction of a pattern recognition algorithm that can identify the most probable location of products inside each cell where 0 < f < 1 (a general exposition of methods of structural pattern recognition is given in Pavlidis [22]). By comparing the value of f in a cell with its value in the four neighboring cells, the scheme recognizes four possible orientations of the local interface between reactants and products and where the latter is located. These are a vertical interface, two vertical interfaces forming a passage, a corner made of a vertical and a horizontal line, or a horizontal interface. The transpose of all possibilities can be easily accommodated in the algorithm. This modification helps to improve the accuracy of the calculations of fluxes across
the boundary in terms of the most probable shape of the local interface in the cell. Thus Eq. (4) is modified to dft,j dt
1
~ x {~uf} x/_ 1/2___1 {3of} ys+ 1/2 l/z Ay
yj i/2 '
(5)
where 0 _<_ a , 3 < 1, both depending on the expected orientation of the interface. The value of o+ or 3 is 1 for a vertical interface or a passage, and is the ratio of the height of the corner or the horizontal interface to the cell size for the last two cases. The introduction of the pattern recognition algorithm helps keep the numerical diffusion bounded since a sharp transition is identified with any change in f . Moreover, the transition between f = 1 to f = 0 is maintained in the minimum number of cells, thereby keeping the body of products coherent. More details are presented in Chorin [17], Ghoniem et al. [15], and Sethian [18]. Since the method of pattern recognition is most accurate when the motion is in one direction, the method of fractional steps is applied to reduce Eq. (5) to two onedimensional equations:
dft,j dt
 Ax I ~~ I x,_ 1/2'
df~,j_ dt
l {3vf} ~/,/z Ay t/2'
~.
l
•
t l f ' ~ xt+ 1/2
(6)
Where f is the intermediate value o f f and time integration is performed using a first order explicit scheme. The stability condition on the time step is (t,j ÷ I/2)
ft, J I t I,/*£, ...m... kl (t~J)
At < ~u I.t ÷ I/'2,jl
(L,J[]~)
Fig. 2 Staggered grid used m the hne mterface calculations.
min
l
lUlmax ,
x
"
(7)
This condition can be relaxed if one time step is much less than the others; e.g., if IOlmax l ul ..... then Aty >> Atx and calculations may be advanced a number of time steps in the xdirection for each step in the y direction. The accuracy of the time integration of Eq. (5) is improved by reversing the order of application of Eqs. (6) every time step to remove the
TURBULENT FLAME PROPAGATION possible directional bias that can result from a consistent integration in one sequence of fractional steps. This also preserves the symmetry of the results if u = v across the entire field. The pattern recognition algorithm, which is applied to identify the interface, produces exact results if all cells h a v e f = 0 or 1. Thus, keeping the number of cells with fractions small should improve the accuracy of the calculations. This can be achieved if the inequality in Eq. (7) is replaced by an equality, which, in turn, improves the resolution of the computations and produces a sharp interface. This situation is similar to the conditions that arise in the solution of hyperbolic equations where the use of the maximum allowable time step results in an exact solution. 111.2. A d v e c t i o n in a Recirculating F l o w The line interface algorithm was applied to identify the boundary of the recirculation zone established behind a rearwardfacing step in a channel. The expansion ratio across the step is 1 : 2. The velocity distribution over the mesh boundaries, (u, o), was obtained from the calculations of Ghoniem and Sethian [23], in which the random vortex method was applied to compute the vorticity field in the channel over a wide range of Reynolds numbers. The evolution of the flow structure in terms of the formation and interaction between vortex eddies was discussed in detail in the same reference. Four values of the Reynolds number, R = 50, 125, 500, 5000, were used to represent four distinct regimes of recirculating flow, namely, viscous, laminar, transition, and turbulent flow. The Reynolds number is based on the step height S, the approach flow velocity U, and the kinematic viscosity of the reactants ~,. A steady velocity field was evaluated by averaging the stationarystate, timedependent velocity over 500 time steps. These results were then used to construct the streamline plots of Fig. 3. The average velocity field was used, in conjunction with Eq. (6), to advect the field of products and evaluate the steadystate mixture distribution. Initially, products were confined to
327 =
125
R =
500
R
R = 5000
Fig. 3. Steadystate streamlines of the flow behind a sudden expansion in a channel at three values of Reynolds number representing the following three regimes: laminar, transition, and early turbulent flows.
the lower half of the channel, f ( 0 , x, y) = 0, y > 1 . O a n d f ( O , x , y ) = 1 , y < 1.0, w h e r e y = 1.0 corresponds to the step height. Since the boundary of the recirculation zone is the streamline that starts at the step and ends at the reattachment point, and fluid cannot cross a streamline, products should be trapped within the recirculation zone at steady state. Figure 4 shows the boundary of the area containing products at steady state, corresponding to Reynolds numbers of 125,500, and 5000. In all cases, an Eulerian mesh that divides the field into 20 × 100 cells was used in a channel with a length five times it width. Thus, Ax = Ay = 0.1, where dimensions are scaled with respect to the step height. The corresponding Couranttime step At = 0.1, if velocities are scaled with respect to the approach velocity of the flow in the inlet channel. By comparing Figs. 3 and 4, it is clear that the boundary of the
Fig 4. The bounding streamline of the recirculation zone, calculated using the velocity field of Fig. 3 and the hne interface algorithm.
328
AHMED F. GHONIEM
recirculation zone is well reproduced by the advection calculations. The numerical error, measured by X Af, j , where the summation is performed over the whole domain and A signifies the change within one time step, reached 0.010.001 for all calculations, representing the error introduced by the algorithm in material conservation.
111.3. B u r n i n g o f R e a c t a n t s
The displacement of the flame surface due to the burning of reactants is calculated by a direct extension of the previous advection algorithm using the following analogy between flame propagation and wave motion. The motion of products due to a velocity that has the same value on both sides of the interface is equivalent to allowing the flame to propagate in the direction of the velocity. If this process is repeated in all possible directions without removing products from their original location, and the maximum advancement of products in all direction is computed, then the expansion of a flame kernel can be evaluated without any explicit reference to the geometry of its front. This amounts to the application of the Huygens wavepropagation principle to compute the flame motion using the advection of products at the burning velocity without the need to determine the local normal to a surface that has a complex geometry. AlgorithmiCally, Chorin [17] implemented this scheme by using Eqs. (6) with (u, v) calculated from
u(x,+l/9=u(x,l/9=S,, cos 0,
v(yj+l/2)~l)(yj_l/2)=Su
sin O,
(8)
where 0=(kTr)/2,
k=0,
1, 2, .  ' , 7,
and the volume fractions at the end of the time step are calculated according to
f ( t + At) = m a x ( f (t), fo, f l , . . . , f7),
(9)
where fk is the field of fractions that corresponds to 0 = krc/2. In twodimensional calculations, one resorts to a fractional step strategy to extend the onedimensional calculations to flame propagation in two perpendicular directions using Eqs. (6). To eliminate the possible directional bias that results from a consistent application of one order of splitting at all time steps, Sethian [18] suggested that the calculations should be performed for the two possibilities(x  y) and ( y  x) sequencesand the maximum of the two iterations should be considered as the total propagation. This algorithm does not require the evaluation of the normal to the flame front, a process that would normally involve the differentiation of unsmooth data provided by the advection step. Moreover, the difficulties associated with the existence of several intersecting fronts within a small area are eliminated. The process of folding of flame fronts and the concomitant formation of islands of one medium within the other is, hence, permitted. This allows the application of the model to regimes at which the Damkohler number is lower than its value for the regime of a single reaction sheet. The significance of this regime was explained by Klimov [24] as the transition between the reaction sheet mode and the distributed region mode of turbulent combustion. The Courant time step which is necessary for the stability of the computations of burning is given by Eq. (7) with u = v = Su. The laminar burning velocity Su can be described as a variable in the computations. Its value may assume an arbitrary function of space coordinates depending on the local flame curvature and stretch [25], as well as the thermodynamic conditions such as the temperature of the mixture and the quenching at the walls. This variation can be important in determining the stable shape of the flame front. (In this regard, the use of the Huygens principle to evaluate the propagation of the flame by burning can provide stable calculations for unstable phenomena. Zeldovich [26] speculated that the interaction
TURBULENT FLAME PROPAGATION between thermodiffusive instability and Landau instability controls the final shape of propagating flames.)
329 the field. Thus /~ = LfS~,
(1 O)
where Lf is the flame surface area per unit 111.4. E x p a n s i o n o f a Flame Kernel
The burning algorithm was applied to study the propagation of a flame kernel that starts from an ignition point located at the closed end of a channel, as shown in Fig. 5. Conditions across the walls were taken as symmetrical in terms of f . The analytical solution of this problem can be obtained in terms of an expanding circle whose center is fixed at the lower corner of the channel where ignition started, while the radius is expanding as R(t) = Sut. The burning rate/9 is defined as the rate of formation of products in
width. If the flame propagates as a perfect circle inside the channel, then 71"
~ SuEt,
B=
(ll)
Su2t s i n  l ( ~ t
Fig, 5. Flame propagation m a channel due to burning only; ignition starts at t = 0 at the lower lefthand corner. The flame front is defined by the boundary between products (cells with diagonal lines) and reactants (clear zones).
) ,
t>H/S~,
where H is the channel height. Figure 6 depicts the comparison between the exact solution for B and the value calculated from the results of the computations. The latter is obtained from the following expression for the rate of burning: B = (Ax)2 ZAf~ j, At
~i~i~!~~~:,.
t
(12)
where A signifies the change over one computational time step. Computations were performed on the same grid which was described in the advection calculations with Ax  0.1. Different values of At, all fractions of the Courant time step Atc as defined by Eq. (7), were used to check on the convergence of the algorithm. Results in Fig. 6 indicate that the Courant time step yields the most accurate solution. In the first phase of propagation, when the burning rate is increasing linearly with time, the refinement of the time step increases the error in the calculations of the burning rate. In the second phase, after the flame reaches the top wall, the effect of the boundary condition becomes more pronounced, causing deviation from the analytical solution that assumes that the flame front remains a perfect circle. Meanwhile, as the time step becomes smaller, the solution exhibits an oscillatory behavior near the asymptotic value of the burning rate. Two factors contribute to this behavior: (1) as the time step decreases, smaller fractions appear in neighboring cells, thus increasing the uncertainty about the interpretation of the interface geometry; and (2) the
330
AHMED F. GHONIEM
Courant time step provides an exact solution for the eikonal equation that describes the propagation of the flame as a wave front, and smaller fractions of the time step increase the integration error.
this work, we apply a constant S, across the whole field and all boundaries are assumed reflective; i.e., all walls are adiabatic.
111.5. Boundary Conditions
Flame propagation by the combined effect of the advection of products and burning of reactants is applied to study turbulent combustion in premixed gases inside a channel with a sudden expansion. The line interface algorithm described in the previous section is employed to follow the flame motion by both mechanisms in a fractional step scheme. The following features of the algorithm are noted:
Flame propagation along solid walls, where the advection velocity u = O, is strongly affected by the boundary conditions in terms of either the temperature or the heat flux across the surface. Since the advection velocity is very small. laminar burning is the primary mechanism of flame propagation at boundary areas, and heat transfer conditions can either accelerate the flame or cause extinction. Ghoniem et al. [15] showed that for flame propagation at nearly isobaric conditions, T = Cf, where C is a constant that depends on the temperature ratio across the flame and T is a reduced temperature. Therefore, for an adiabatic wall, the proper boundary condition is fw+l = fw~, if w is the cell at the wall, and (w + 1) is an imaginary cell outside the physical domain. On the other hand, if the wall is isothermal, then fw+l should be maintained at the value determined by Tw. In
IV. S O L U T I O N S
1. The scheme provides a uniform resolution over the whole field without the need to add extra points at areas where the flame front geometry becomes exceedingly complicated. This is done without using sorting routines to define neighboring points. 2. While the resolution is limited by the cell size, changes over a distance smaller than one cell dimension are identified by the pattern recognition algorithm. 3. Calculations are confined to cells where an
~o
///
///
///
0.5 ?
~At
!!
18 Ate/8
'
Ate/2 ExocfTM 0'5
I'.0
~ ,,,~
Ate/8 ~tC/8
115
2.0
I
Fig. 6. Burning rate versus time for flame propagation In a quiescent medium The analytical solution is indicated by a thick hne, and the numerical solutions at different values of the time step are shown by thin hnes. At, ,s the Courant time step
TURBULENT FLAME PROPAGATION interface exists; thus although the algorithm employs an Eulerian grid, it has the same efficiency as a Lagrangian calculation that tracks the flame using a set of markers along its front. 4. The algorithm allows for the presence of convoluted fronts that form folds and islands since neighboring points are defined in terms of their coordinates instead of mutual links. Markerincell techniques can only be employed until the front folds in one spot and then they require extra care in defining multivalued interface coordinates. 5. The calculations of burning normal to the front surface do not require identifying the normal to the surface. Computations of the normal to a surface defined by a set of points can pose difficulties associated with the inaccuracy of the numerical differentiation. 6. The calculations are stable provided that the Courant condition is not violated. The splitting scheme allows the use of the onedimensional version of the stability requirement where the time step is chosen as the minimum of the set of maximum values permitted by the different onedimensional steps. The same advantage is gained when computation of several modes of transport is performed simultaneously. In this case, one can use the time step defined by the mechanism that causes the fastest change in the concentration field.
331 of similar order of magnitude which amplifies the exchange across the flame surface. To distinguish between the two effects of the Reynolds number in terms of the geometry of the flow field and the size of turbulent fluctuations on the flame propagation, two cases are computed: (1) steady propagation in which the average velocity distribution is used at all time steps; and (2) unsteady propagation where the timedependent stationary velocity field is used to advect the flame front. Results are presented in terms of the flame front geometry at different values of the laminary burning velocity. Experimentally, changes in the laminar burning velocity Su are achieved by varying the equivalence ratio of the airfuel mixture ~b, where ~b is defined as the actual fuelair ratio divided by the stoichiometric fuelair ratio (the correspondence between Su and 4~ is obtained assuming a propaneair mixture and using an approximate fit for the data of Andrews and Bradley [27]). Several fractions of the Courant time step, defined in Section III, were used until the value of B converged to a constant value. At that limit, the concentration field exhibited a sharp transition between f = 0 to f = 1 within one cell in either direction. The geometry of the flame front was then constructed from the interpretation of the data of the concentration field as seen by the pattern recognition algorithm.
IV.2. Propagation in a Steady Field IV.1. Turbulent Flame Propagation Ghoniem and Sethian [23] show that the Reynolds number affects the flow field in two aspects: (1) the shape of the streamlines, where the size of the recirculation zone grows to a maximum around R = 500, corresponding to transition, and then decreases to a smaller value at turbulent conditions; and (2) the relative magnitude of turbulent fluctuations to the average velocity field. Results indicate that turbulent fluctuations generate a flux of momentum that grows with the Reynolds number. Associated with the momentum flux is a transport of mass
At each Reynolds number, the timedependent velocity field was averaged over 500 time steps and the results were used in the propagation of the flame during the advection step. Figure 7 shows the geometry of the flame front as delineated by the interface between reactants and products for Reynolds numbers 125, 500, and 5000, at values of S, in the range of 0.010.2. These lines are the trace of a ragged interface similar to the boundary between shaded areas and clear areas in Figs. 911. In all cases, the flame reaches an equilibrium position when the downstream advection of products at its front is balanced by the upstream transforma
332
AHMED F. GHONIEM Su~O 2
015
r R=125 Su=O 2
0 I j
001 0 0 2 0 0 5
o 15
F Su=O 2
0 15
V R=5000
sin 6 = S , / l u I,
o I
J R=500
,
OOI 0 0 2
half of the channel, the front is almost a straight line away from the anchoring point. The front angle, 8, is given by the kinematic condition u. n = S , , which yields
005
0 t
001 0 0 2 0 0 5
Ftg. 7. The location of the flame front at steady state for different values of the lammar burning velocity at the Reynolds number correspondmg to three v~scous regimes.
where lu[ is the absolute value of the local flow velocity. 5. The flame length within the computational domain remains almost constant, 1011 step heights, until S,  0.1 when the tip of the front touches the top wall of the channel. Within this range, i.e., S, = 0.010.1, the total burning rate increases linearly with the laminar burning velocity according to Eq. (13), as shown in Fig. 12. At S, > 0.1, the flame length decreases as S, increases so that
LfSu h U = constant, =
tion of reactants by burning in the direction normal to its front. Thus, the criterion for steadystate propagation is
(X~Xf,,)adv + (XAf,,)bur = 0, where adv denotes the advection step and bur indicates the burning step. The following observations can be gathered from these plots: 1. The flame approaches the mixing channel as S, increases. 2. Since the flow is steady, the flame propagates out of the recirculation zone at the smallest values of S, and stabilizes itself at a distance away from its boundary. 3. The distance between the flame front and the bounding streamline of the recirculation zone is proportional to S,. The flame stands at a point defined by the vector sum of the local flow velocity and the burning velocity. Therefore, while at low Su the flame is almost locally tangent to the streamline, at high S, it is approximately normal to the mean flow. 4. The front exhibits slight curvature round the point of separation where the velocity of the flow changes rapidly. At high laminar burning velocities when the flame propagates out of the recirculation zone and into the upper
(13)
(14)
where h is the incoming channel height and U is the approach uniform velocity. Within this range, the burning rate remains constant at its maximum value. Using a steady velocity distribution throughout the range of Reynolds numbers implies that turbulence, in particular the fluctuation associated with the growth of perturbation at high Reynolds number, is not allowed to affect the flame propagation. Therefore, all the previous observations are due to the kinematic effects of the flow field as manifested by the velocity profiles and the recirculation zone. The propagation of the flame front toward the mixing chamber at high values of the laminar burning velocity S,, which correspond to equivalence ratios near the stoichiometric conditions, was observed experimentally by Ganji and Sawyer [28]. The process was initiated by a series of events that led to the flashback of the flame into the incoming channel. Experimental cinematographic records show that as long as the entire flame front is contained within the channel and away from the top wall, flashback, or limited upstream propagation, is not possible. However, as soon as the flame touches the top wall at the downstream end by increasing the equivalence ratio to the appropriate value, a sudden
TURBULENT FLAME PROPAGATION
333
transition to front oscillation followed by an upstream propagation was encountered. Using the numerical results shown in Fig. 7, although lacking both the unsteadiness of the actual flow and the dynamic effects of expansion associated with burning of rich mixtures, one can postulate the following mechanism for the interpretation of the experimental records. For lower values of Su, in the range of 0.010.1, changes in the laminar burning velocity do not affect the front location substantially. For values larger than 0.1, the flame angle at the step increases and its other end touches the upper wall, moving closer to the inlet channel at the slight change in S,,. The distance X between the flame and the step can be obtained from satisfying the continuity across the flame described by Eq. (14), coupled with geometrical consideration at the front, f , as presented by Eq. (13). Assuming uniform conditions across the upper half of the channel, i.e., u ( y ) = U, the following relationship can be obtained: Su
(15)
x  I~/TS_s :
The values obtained from this equation agree with the numerical results of Fig. 7 to the extent that u ( y ) is uniform across the upper half of the channel. The propagation of the flame along the top wall causes heating of reactants which, in turn, increases the value of S.. Equation (15) indicates that the attachment length of the flame at the top wall X decreases by this increase in S.. This process is augmented by the presence of a boundary layer along the top wall, which slows down the flow of reactants. The displacement of the flame along the top wall continues till the flame reaches the incoming chamber. The mechanism suggested here is selfsustained, and is initiated by the flame moving out of the recirculation zone, increasing its angle at the step and reaching the top wall. IV.3. P r o p a g a t i o n in a Stationary Field
To investigate the effect of turbulent fluctuations on flame propagation and burning rates,
the unsteady velocity field is used in the advection of products behind the front. According to Ghoniem and Sethian [23], the crossstream velocity component remains without a noticeable change throughout the range of Reynolds numbers from 125 to 5000, while turbulent kinetic energy in the same direction exhibits a substantial increase starting at R = 500. Of particular interest here is the rise in the turbulent shear stress  u ' o ' , which is accompanied by a material flux normal to the main stream. The latter enhances the mixing between products within the recirculation zone and fresh reactants in the main stream. Samples of the instantaneous streamline plots for three values ofR~, depicting typical vortex structures which form at different regimes, are shown in Fig. 8. Plots of the flame front are shown in Figs. 9, 10, and 11 for R = 125, 500, and 5000, respectively. For each Reynolds number, values of the laminar burning velocity are varied from 0.01 to 0.2 to study the effect of S,, on the flame geometry and the stability of combustion. The plots identify the flame as the interface between clear areas, depicting reactants, and shaded area, representing the products, while the interR
=
125
R
=
500
R
=
5000
Fig. S. Plots o f the instantaneous s t r e a m h n e s o f a stationary flow behind a sudden e x p a n s i o n in a channel at three values o f R e y n o l d s n u m b e r . At R = 125. the flow is l a m i n a r with a r e c i r c u l a t i o n zone f o r m e d o f two eddies. At R = 500, the flow is at translUon with eddy s h e d d i n g from the s e p a r a t i n g s h e a r layer. A turbulent flow is established at R = 5000
334
AHMED F. GHONIEM
Su =
l~:.:
"
O.8I
i~ .....
5u
2::::::

So
=
I].01
Su
=
8.02
$u
=
8.~5 =~mmmmlmtl
~::::
= 8.02
So : 8 . 0 5
~
5u =
8.18
So =
8.18
~':::i~i:~=:~:i::~::?!!)~:~,~:4:~: ..... ~=~:(~ ...... i~. ~::: 5u : B.15
:::~~. =========================================================
Su =
..... Jx,~ 5u = 8 . 8 2

= B.85
:
5u
:!~:: =
8.10
_ ~ ' ~  : : : :
S~ : 0 . 1 5
Su =
[email protected]
l0 F l a m e front c o n t o u r s for R = 500, at dsffercnt v a l u e s of S,.
Fig
"~ ......
~;~
::5~ii::~,: . . . . .
So = B.28
5u = 0.81
.......
:
~..~:,:::~~:~_::~~i ,7 i i ; i i i ~ ,:~;,
0.28
Fig. 9. F l a m e front contours, represented by the interface between reactants and products, for R = 125 at different v a l u e s of the laminar burning velocity S,. Plots of the local interface are obtained from the pattern r e c o g n m o n alg o r i t h m Shaded a r e a s are used to depict products.
5u
e, Su = 8 . 1 5
Fig []. Flame front
contours
for R
=
5000,
at
different
values of S,.
face is constructed using the algorithm interpretation of the concentration field. The following characteristics of turbulent flame propagation in premixed gases are noted from the figures: 1. At low laminar burning velocities and high Reynolds numbers, the combustion process is disrupted at the downstream side of the recirculation zone. Without maintaining a pilot flame at the edge of the step, modeled in the computations by a constant source of products, the flame would disappear. The breakup of the recirculation zone into a set of interacting eddies is due to the instability of vorticity within the recirculation zone at high Reynolds numbers. Moreover, the small values of the laminar burning velocity S, do not allow the rate of formation of products to maintain a "hot" zone near the step. The disruption of burning at the downstream side propagates toward the step and leads to a total extinction of the flame; this is the blowoff that has been experimentally observed in [28] at lean conditions. 2. At intermediate values of S,, the wrinkling of the flame surface by large scale eddies,
TURBULENT FLAME PROPAGATION which are shed from the separating layer at the step, provides the mechanism of enlargement of the combustion zone that leads to a substantial increase in the total burning rate. Moreover, the wrinkling of the flame surface remains of the same order of magnitude after the front moves out of the recirculation zone and into the free stream, indicating that the field produced by vorticity persists for long distances outside the large eddies. 3. The total flame length that affects the rate of burning in Eq. (10) consists of the apparent front between the incoming reactants and products, plus the boundary of the islands that are trapped within the products. That provides two different, related mechanisms for explaining the increase in the burning rate at high Re. On the average, the increase in the burning rate is due to the generation of turbulent mass flux in the crossstream direction. From an unsteady point of view, it is the continuous wrinkling of the flame front due to the passage of the large eddies which causes the front to stretch and allows extra flamesurface area. 4. At large laminar burning speeds S,,, corresponding values of ~ near the stoichiometric conditions, the flame angle at the step increases while its other end reaches the top wall. A slight increase in S, causes the flame to propagate along the wall toward the incoming flow. The burning rate reaches its maximum value defined by Eq. (14). The possibility of flashback, as described in the previous section, is now compounded by the unsteadiness of the flow. (The oscillation of the flame surface may generate pressure waves that propagate into the incoming flow, hence slowing it down and increasing the apparent value of S,,. This, however, is not included in the analysis presented here.) These observations support the postulate that the range of stability of turbulent flames in premixed gases is limited by the fluid mechanical processes associated with recirculating flows. If the disruption of the flame is interpreted as a lean blowoff limit, while the flame reaching the top wall is a nearstoichiometric flashback limit, then from the plots in Fig. l l
335 for R e = 5000 the stable range of operation is 0.01 < S,, < 0.05, which corresponds to an equivalence ratio q~ in the range of 0.4 =< ~ < 0.8 based on the data of Andrews and Bradley [27]. The experimental results by Ganji and Sawyer [28] confirm these conclusions about the narrow range of stability of turbulent flames in premixed gases. IV.4. Total Rate of Burning The total rate of burning, defined by Eq. (12) and nondimensionalized with respect to the flow rate, was computed for all cases considered above and plotted in Fig. 12. Points on the curves were obtained for S, = 0.01, 0.02, 0.05, and 0.1. For the unsteady case, an average value was evaluated after reaching a stationary state. The line representing burning under steady conditions is almost straight, indicating that the flame length stays approximately constant until the front touches the top wall when/Y reaches a maximum equal to ( h U ) . Thereafter, the flame length decreases to a value inversely proportional to S,. As the Reynolds number increases, the flame stretches along the large eddies that form the separating shear layer between the incoming stream and the recirculation zone. Higher burning rates are obtained due to the increase in the flame length, which can also be correlated with the rise in the turbulent statistics. I0

~'O °51
O0 0
I
005
0
Su
Fig. 12. Total burning rate B at different values o f laminar burning velocity S,,, corresponding to varmtlons of the equivalence ratio ~b, lbr steady and timedependent flows at three Reynolds numbers R.
336
AHMED F. GHONIEM
V. C L O S U R E The stability limits for the propagation of a turbulent flame in premixed gases at high Damkohler numbers, in a regime where burning can be modeled by a thin reaction sheet separating reactants and products, are evaluated using line interface calculations of advection and burning. To emphasize the Reynolds number effects while limiting the computational effort, the dynamic processes associated with the density jump across the flame surface are neglected. Thus, the interpretation of the results can only provide qualitative trends and should be supplemented by a comprehensive model of these processes if detailed predictions are attempted. Numerical results show that at high Reynolds numbers and low laminar burning velocity, the combustion process is extinguished. The anchoring mechanism is disrupted by the imbalance between the rate of removal of products from the recirculation zone and the rate of burning. At high burning velocities, the flame propagates out of the recirculation zone and along the top wall of the channel. The heating of the gas near the wall accelerates the flame by increasing the laminar burning velocities of the reactants. At moderate values of the laminar burning velocity, turbulence increases the exchange of mass across the recirculation zone, thus providing higher rates of burning.
5. 6 7. 8 9.
10
11.
12.
13. 14. 15.
16.
17, 18. 19.
20.
This work is supported by the National Science Foundation under Grant CPE8404811 and the Air Force Office of Scientific Research under Grant AFOSR840356. Computations were performed by James Sethian and Omar Knio.
21 22. 23.
REFERENCES 24. 1
2.
3.
4.
Mellor, A M., and Ferguson, C. R., in Turbulent Reacting Flow (P. Libby and F. Wilhams, Eds.), Topics m Applied Physics, 4, SpringerVerlag, 1980, pp. 4564. Cheng, S. I., and Kovitz, A A., 7th Symposium (International) on Combustion, Bunerworths, London, 1959, pp. 681691. Bray, K. N. C., m Turbulent Reacting Flow, (P. Libby and F Wllhams, Eds.), Topics in Applied Phystcs, 4, SprmgerVerlag, 1980, pp. 4564. AbdelGayed, R. G., and Bradley, D., Phil. Trans. R. Soc. Lond. 301.125 (1981).
25.
26. 27. 28.
Damkohler, G., Z. Elekrochemie Angewandte Phys. 46:601626 (1940) (NACA TM 1112, 1947). Chomlak, J., Prog. Energy Combust, Sci. 5:207221 (1979). Ballal, D. R., and Lefebvre, A. H., Proc. R, Soc. Lond. A344:217234 (1975). Cantwell, B. J., Ann. Rev. FluidMech. 13:457515 (1981). Yule, A J,, Chlgier, N. A., and Thompson, D., Symposium on Turbulent Shear Flow, Penn. State Umversity, 7, 1977, p. 41. Ashurst, W T., in Proc. 1st Syrup. on Turbulent Shear Flows (Durst et al.. Eds.), SpringerVerlag, Berlin, 1979, pp. 402413. Ghoniem, A. F., Chorm, A. J., and Oppenheim, A. K., 18th Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, 1981, pp. 13751383. Ghomem, A. F., Chen, D. Y., and Oppenheim, A. K., 22nd Aerospace Sciences Meeting, AIAA Paper No. 840572, Reno, Nevada, January 1984. McDonald, H., Prog. Energy Combust. Sci. 5:97122 (1979). Williams, F, A., AIAA840475, 1984. Ghoniem, A. F., Chorin, A. J., and Oppenheim, A. K., Phd. Trans. Royal Soc. Lond. A304:303325 (1982). Majda, A., and Sethian, J. A., Report PAM197, Center for Pure and Applied Mathematics, Umversity of Cahfornia, Berkeley, 1984. Chorin, A. J,, J. Comput, Phys. 35:111 (1980). Sethian, J. A,, J. Comput. Phys. 54:425456 (1984). Hslao, C. C., Ghomem, A. F., Chorin, A. J., and Oppenheim, A. K., 20th Symposium (International) on Combustion, Ann Arbor, MI, 1984. Noh, W. T., and Woodward, P , 15th International Conference on Numerical Methods in Fluid Dynamics (A, Vooren and P. Zandbergen, Eds.), SprlngerVerlag, Berlin, 1976, pp. 330339. Hlrt, C. W., and Nichols, B. D., J, Comp. Phys. 39:201225, (1981). Pavlidis, T., Structural Pattern Recognition, SpringerVerlag, New York, 1977. Ghoniem, A. F., and Sethian, A. J., 23rd AIAA Aerospace Sciences Meeting, AIAA850146, Reno, Nevada, 1985. Klimov, A. M., Dokl. Akad. Nauk., SSSR 221:56 (1975). Williams, F. A., 1984 Technical Meeting, Eastern Section of the Combustion Institute, Clearwater Beach, FL, 1984. Zeldovlch, Ya. B., Combust. Flame 40:225234 (1981). Andrews, G. E., and Bradley, D., Combust. Flame 19:275288 (1972). Gangi, A. R., and Sawyer, R. E., A I A A Journal 18:817824 (1980).
Recetved 22 March 1985; revtsed 5 December 1985.