Dynamical behaviour of a premixed turbulent open V-flame

The level-set approach of Osher & Sethian to tracking interfaces is successfully adapted to the simulation of a premixed turbulent open V-flame including the effects of exothermicity and baroclinicity. In accord with experimental observations this algorithm, along with a flame anchoring scheme, predicts flame cusping for a case in which a strong vortex pair interacts with the flame front. The computed velocity and scalar statistics obtained for the turbulent V-flame compare reasonably well with experimental results by Cheng & Shepherd, and demonstrate the importance of flame-generated vorticity in the determination of flame dynamics and product velocity characteristics.


Introduction
The propagation of a flame and its interaction with a surrounding turbulent flow field in a premixed medium is one of the fundamental and challenging problems of combustion phenomena. The coupling of the chemical reaction with the hydrodynamic flow field makes the problem difficult to analyse, although many models are now available and moderately successful calculations (Ghoniem, Chorin & Oppenheim 1982;Sethian 1984;Ashurst 1987) have been conducted. Some models have relied on statistical methods using a probability density function (e.g. Pope 1976). This p.d.f. is specified empirically or may be derived from an evolution equation which requires certain modelling assumptions. Alternative descriptions have been based on the laminar flamelet concept (Peters 1986), wherein the flow field consists of a collection of flame elements embedded in the turbulent stream. The structure of these flamelets (Clavin 1985) is analysed separately from the flow-field description so that the complicated chemistry problems are decoupled from the simulation of the turbulent flow field. In many practical situations, the flame thickness is smaller that the smallest turbulent length scale and flamelet descriptions are relevant.
Here, we simplify one step further beyond the flamelet model and take the zero limit for the flame thickness and solve only for the geometry of the flame front. This is commonly referred to as the reaction sheet limit. The flame front acts as an infinitesimally thin boundary which separates regions of constant density and propagates into the fresh mixture at a prescribed flame speed. The speed is assumed to be a function of the local curvature of the flame front (though not of local strainsee later discussion), which avoids hydrodynamic and diffusive-thermal instability (Darrieus 1938;Landau 1944). A Lewis number greater than unity is implied in the analysis. Fluid elements at the flame front undergo an increase in volume as they burn, creating a jump in the normal component of velocity and a concomitant decrease in density, assuming constancy of pressure across the flame front. This density drop combined with a non-uniform velocity field along the flame creates vorticity and this new vorticity field contributes to the advection of the flame.
In the reaction sheet model, which is discussed in detail in the recent book by Liiian &Williams (1993), the flame is treated as a discontinuity without internal structure and consequently there is no intrinsic length scale, namely the laminar flame thickness. In this infinite Damkohler number limit chemistry factors, such as the rate equation, transport effects and energy balance are perforce excluded from the model. The density ratio across the flame and the normal laminar flame speed are chosen a priori as input parameters based on experimental data appropriate to the situation to be simulated. Quenching is likewise excluded from the model, since it involves the balance between reaction rate and strain rate, as represented by the Karlovitz number which is given by the ratio of the laminar flame thickness to the normal laminar flame speed multiplied by the flame stretch. The reaction sheet model is therefore the zero Karlovitz number limit. Quenching is customarily expected for Karlovitz number exceeding unity, the Klimov-Williams criterion, although recent work by Poinsot, Vegnante & Candel (1991) indicates that local quenching of premixed adiabatic flames is an unlikely occurrence and that this criterion may underestimate the Karlovitz number quenching boundary in combustion diagrams by as much as an order of magnitude. Our model simulation may therefore have applicability to a wide range of practical premixed combustion situations, such as are enumerated by Liiiin & Williams. Pindera & Talbot (1986) proposed a computational model for premixed turbulent combustion and concluded that production of vorticity at the flame front forms an integral part of the overall velocity field. Their model was based on solving for the turbulent flow field using the discrete vortex method developed by Chorin (1973). The flame front was represented by a set of marker particles. A regridding procedure was employed when marker particles come together in regions where the curvature of the propagating front is large (Hyman 1984). This regridding procedure, however, resembles diffusion and dominates the actual effects of curvature. In view of the fact that propagating flames can develop cusps, we adopt a level-set algorithm, introduced by Osher & Sethian (1988), for following fronts propagating with curvature-dependent speed (Sethian 1984). The algorithm handles cusping of the flame front naturally and provides an accurate calculation of the curvature which is needed for the evaluation of the flame speed and the estimation of flame-generated vorticity. The level-set field equation which governs the flame motion is itself not new, and has been employed by others (Kerstein, Ashurst & Williams 1988) and is often referred to in the literature as the 'G equation'.
Many experimental studies have been performed to examine the interactions between fluid turbulence and the combustion process. The open V-flame stabilized by a rod is one of the simple configurations for which various velocity and scalar statistics have been reported (Cheng 1984;Cheng & Shepherd 1986;Shepherd, Cheng & Goix 1990). Their results showed that increases in the unconditioned velocity fluctuations and the sharp increase in the Reynolds stress in the flame zone are caused by intermittency contributions associated with the difference in the mean velocities in the reactant and the product regions. They also obtained conditioned velocities and Reynolds stresses in premixed V-flames using a conditional sampling technique to investigate the true nature of flame-generated turbulence and these conditioned statistics have been used for comparison with numerical results by vortex dynamics models (Pindera & Talbot 1986) and the BML model (Bray & Moss 1977;Bray, Libby & Moss 1984).
The primary objective of the present study is to adapt the Osher & Sethian level-set algorithm to the case of a rod-stabilized premixed V-shaped open flame. To adapt the level-set algorithm to this problem, several extensions to the basic numerical schemes are required. They include : ( a ) the effect of exothermic volume generation on the flow field, due to the temperature rise across the flame front in the constant-pressure flow; (b) the effect on the flow field of baroclinic vorticity generation, described by the Hayes (1959) formulation of the vorticity jump condition, thus extending the earlier work of Pindera & Talbot on this effect; (c) the effect of the free-stream turbulence on the flame dynamics, which is modelled by application of Chorin's random vortex technique.
The numerical studies to be described allow the separate evaluation of the importance of each of these effects, and then their combined effect on the overall flame dynamics. The case chosen to be studied in most detail corresponds to the experimental conditions of the investigations by Cheng (1984) and Cheng & Shepherd (1986), and comparisons between their experimental results and the predictions of the present numerical model are presented. A number of computed flow-field properties and turbulence statistical quantities are compared with experimental results. In addition, a separate numerical study is presented describing the interaction of two strong discrete vortices with the flame front and which specifically exhibits the capability of the algorithm to predict the flame cusping phenomena which have been observed experimentally.
In $2 we describe the theoretical aspects of our model. Detailed numerical implementation is described in $3. We apply the above ideas to a rod-stabilized open V-shaped premixed flame in $4. These numerical predictions are compared with experimental results of Cheng & Shepherd in $5. We summarize the results in $6.

Physical model
In this Section we present our model for flame propagation and its accompanying flow field. The main assumptions of the model are: (a) the flow is two-dimensional and inviscid; (b) the Mach number is vanishingly small, hence the flow is dynamically incompressible ; (c) the flame acts as an interface which separates two regions of different but constant density and propagates into the unburnt gas at a prescribed flame speed S, which is dependent on the local curvature.

2.1, General background
A common approach to the problem of the interaction between the flow field and the flame consists of separating the fluid dynamic treatment from that of the flame zone by treating the flame front as a surface of discontinuity separating two fluids of different densities. Darrieus (1938) and Landau (1944) found independently that perturbations along a density interface that propagates at a constant speed grow rapidly without bounds and that such flames are absolutely unstable to any perturbation. This is the so-called hydrodynamic instability. Heat conduction improves flame front stability by changing the local value of the burning speed. The concave/convex part of the perturbed flame heats up the reactants more/less than an unperturbed planar flame, thus increasing/decreasing the burning speed. In this work, we model a flame under the stabilizing effects of heat conduction associated with curvature along the flame. Markstein (1964) showed that if the burning speed varies with the curvature of the flame front such that concave/convex flames propagate faster/slower than planar flames, a stabilizing mechanism is realized for shortwavelength perturbations, and which is verified experimentally. According to this formulation, the laminar burning speed for weak curvature is given by where Si is the laminar burning speed of a planar flame, often denoted in the literature as S,, K is the curvature of the flame front, R, = 1 / K the radius of curvature of the front taken positive in the direction of the products, and L is the Markstein length scale. An expression for L was deduced by Garcia-Ybarra, Nicoli & Clavin (1984) using highactivation-energy asymptotics. Two observations are in order with respect to the propagation relation (2.1), the first concerning the value of the Markstein length L and the second regarding the curvature factor. The analysis of Garcia-Ybarra et al. yielded the result that L, divided by the laminar flame thickness, ranges from 2 (rich) to 6 (lean) for ethylene flames. Since in the reaction sheet formulation the laminar flame thickness is undefined, a characteristic length has to be chosen apriori to represent the physical flame configuration simulated, and the only such length that can be unambiguously defined in the computational model is the integral scale of the reactant turbulence which is computed a posteriori. Hence an initial 'guess' has to be made concerning the relationship between the numerical grid size and the physical scale of the flame being simulated, which has then to be later validated. In the present computations, it was assumed that the grid size (= 0.02 in non-dimensional units) corresponded to a physical turbulence integral scale of 1 mm, and that the thickness of the laminar flame being simulated was about 0.5 mm. With the Markstein ratio taken to be the intermediate value of 4, these values together yield the non-dimensional numerical value L = 0.04, which in fact is the same value employed by Ashurst. The 'guess' of the grid size equivalence to 1 mm turned out from the results of the computations to be quite good, and iteration was not deemed necessary.
In the propagation equation only the curvature has been included, when in fact the stretch K (see (2.6~)) which includes strain effects as well should be employed. If K is replaced by K in (2.1), the computation at each time step becomes an implicit one because the volumetric source term depends on the flame speed S,, which has to be evaluated recursively until convergence is achieved. Ashurst used the full stretch term in his computations, and describes a procedure for performing the iteration. However, for the weak turbulence of the present computations, u'/Sz N O(l), where u' denotes the turbulence r.m.s. velocity, the contribution of strain to the overall stretch was at maximum roughly only 10% and was neglected. This finding is consistent with the DNS results of TrouvC & Poinsot (1993) who found that for u'/Si = 10 the two terms were of the same order of magnitude. Since the strain term scales with u'/l, where 1 is a characteristic length of the free-stream turbulence, it is expected that a ten-fold reduction in u ' / S i would be accompanied by a similar reduction in the instantaneous values of the strain at the flame interface. Ashurst likewise reports that the strain term made only a minor contribution to the velocity statistics.
Also, in the case of weak combined curvature and strain stretch, the flame speed is linearly dependent on the stretch, which is the case considered here. However, for moderate or strong stretch ( K of order unity or greater) strain and curvature produce qualitatively different effects. It was confirmed from the results that the computations conformed to the weak stretch regime.

Flow-field description
As was done by Pindera & Talbot (1986), we decompose the velocity field Uinto three components : with the individual components satisfying the following conditions : u= u,y+ u,+ up, is the velocity, Us the velocity field due to volume expansion across the flame front, U, the rotational velocity field due to vorticity, U, the potential velocity field, m the volume source strength per unit length of the flame, x f the coordinate describing the position vector of the flame front, w ( x ) the vorticity field, $ the velocity potential of the incident flow, and 8 ( . ) the two-dimensional Dirac delta function.
We derive an expression for the effect of volume expansion on the velocity field, Us, using the conservation of mass. Let S , and S, be the relative normal fluid velocities with respect to the flame on the unburnt and burnt sides, respectively, and let p, and pb be the fluid densities of the unburnt and burnt mixtures. Mass continuity yields As reactants are converted into products, volume m is generated in the expansion process which per unit length of flame is given by The vorticity transport equation is The last term on the right of (2.5) is the baroclinic torque term and is a source of vorticity through the interaction of gradients of density and pressure. Pressure gradients tangential to the flame cause different accelerations in the light and heavy gases and hence vorticity will be produced at the flame by the mean density gradient across the flame and the pressure gradient tangential to the flame. It should be noted that although the baroclinic torque term involves the pressure gradient tangential to the flame, pressure gradient effects are of second order as regards the flow field due to volumetric sources, and are neglected in the solution of the Poisson equation governing the velocity potential associated with the volumetric sources and the irrotational flow field. Hayes (1959) derived a vorticity production equation in the inviscid limit that involves quantities which are continuous across the flame or which have known jump conditions, and Lagrangian time derivatives of quantities in the tangential direction of where here Us is the flow velocity at the flame in the tangential direction, V s its gradient along the flame, and V, the absolute normal flame speed; d/dt denotes the time derivative taken at a point which always lies on the front and which moves in a direction normal to the discontinuity as it moves. Figure 1 shows a sketch of the normal and tangential directions denoted by the subscripts n and s. Defining the stretch K as where A is an elemental flame front area, we have The unsteady term, dUS/dt, was not included in the analysis of Pindera & Talbot.

Equation for flame propagation
It is convenient to formulate the flame propagation problem in terms of an equation for a scalar field Y ( x , y , t ) which describes the flame location. Following Osher & Sethian (1988), we define a continuous initial distance function Y(x, y , t ) such that Y > 0 in the unburnt region, Y < 0 in the burnt region and the zero-level surface Y = 0 represents a flame front configuration as is depicted in figure 1.
To find an equation for the evolution of Y which corresponds to the propagating flame, consider the motion of some level set Y ( x , t ) = C. Let x(t) be the trajectory of a particle located on this level set, so Y ( x ( t ) , t ) = c.
(2.8) The particle speed ax/at in the direction -n normal to the flame is given by the flame speed S,(K). Thus, with the convention adopted here that n points in the direction of the burnt gas (2.11) 2.4. Flame propagation with advection A flame propagates itself with its own burning speed S,(K) and is also advected by the accompanying flow field. The flame front affects the flow field by volume generation and vorticity production and this flow field influences the flame location by advection. Hence there is a mutual interaction between flame and flow field. In this case the field equation becomes q+ss,(K)lvYl+ U -V Y = 0. (2.12) The second term represents propagation and the third term denotes advection of the scalar field Y. The velocity U is the convection velocity of the unburned fluid just upstream of the flame sheet.

Numerical methods
In this Section, we present the numerical details for an approximation to the equation of motion. To determine the flame motion, a numerical scheme which incorporates the effects of flame propagation, exothermicity and baroclinicity in the level-set algorithm is described. The algorithm we employ uses the method of fractional steps to decompose the motion into propagation and advection. First we propagate the flame as a result of burning. Then we locate the volume sources due to exothermicity and the flame-induced vortices due to baroclinicity to obtain the velocity field for the advection part. A numerical scheme for the flameholder is developed to stabilize the flame against blowout.

Flame propagation with advection
We begin by briefly reviewing the level-set approach to tracking propagating interfaces. Complete details may be found in Osher  and Sethian (1989).
As shown above, the equation of motion for the propagating interface without advection is given by where K is the curvature of the front. We refer to this equation as a 'Hamilton-Jacobi' level-set formulation. Strictly speaking, it is only a Hamilton-Jacobi equation in the case when S , is constant, but the flavour of Hamilton-Jacobi equations is present.
This yields an equation of motion for a higher-dimensional function Y for which a particular level set always corresponds to the motion of the original front. Another way to say this is that we have transformed the Lagrangian equation which would have resulted from a parameterization of the moving interface (the approach taken by Pindera & Talbot and by Ashurst in their V-flame simulations) into an Eulerian equation on a fixed grid of one higher dimension. Thus, if the interface is an (n-1)dimensional hypersurface, we have traded it in for an n-dimensional problem.
Fortunately, the advantages of this exchange far outweigh the additional computational energy required by the extra dimension. To begin, we observe that the function Y(x,t) always remains a function, even if the level surface Y = 0 corresponding to the front changes topology, breaks, or merges. In such cases, parameterizations of the front often break down. As an example, consider two circles in R2 whose distance between centres is initially greater than the sum of their radii expanding outwards with unit normal velocity. The initial function Y is doublehumped. As Y evolves under the Hamilton-Jacobi equation of motion, the topology of the front Y = 0 changes. When the two circles expand, they meet and merge into a single closed curve with two corners. This is reflected in the change of topology of the level set Y = 0.
Thus, the level-set approach avoids the complex bookkeeping that plagues discrete parameterization techniques when the interface changes topology. Another advantage is that the technique is applicable in any number of space dimensions; calculations of interfaces propagating in three space dimensions are discussed in detail in Osher & Sethian (1988) and Chopp & Sethian (1994).
Finally, a crucial advantage of this approach is that, because we have posed an Eulerian problem for the motion of the propagating interfaces, fixed-grid finite differences may be used to approximate the equations of motion. While care must be taken to choose differencing schemes that satisfy an entropy condition for propagating fronts, the most basic versions of the schemes presented in Osher & Sethian are extremely straightforward and simple to program.
It is tempting to use a central difference approximation to the gradient in the Yequation and thus produce the obvious explicit scheme (central difference in space, forward difference in time) for the update Yn+l. Unfortunately, such an approximation is unworkable, for reasons which we now explain. For details, see Sethian (1985), Osher & Sethian (1988) and Sethian (1989).
Consider the simple case of a front propagating with speed function S,  Sethian (1985), followed by a proof in Osher c > 0, the right-hand side diffuses sharp gradients and forces the Y to stay smooth for all time. Conversely, for Y = 0, corners develop, and a singularity develops in the curvature. This situation is analogous to solutions of hyperbolic conservation laws, in which the absence of viscosity on the right-hand side allows the development of shock discontinuities in the propagating solution. Indeed, an entropy condition is required to force the correct solution for propagating interfaces which is equivalent to the one required for hyperbolic conservation laws. A full description of this entropy condition and the parallel between propagating interfaces and hyperbolic conservation laws is given in Sethian (1989).
Thus, an accurate numerical approximation to the equation for a propagating interface must pick out the correct entropy-satisfying solution and avoid excessive smearing at sharp discontinuities. This leads quite naturally to the use of schemes borrowed from the numerical solutions of hyperbolic conservation laws, where stable, consistent, entropy-satisfying schemes have a rich history.
A complete explanation of the use of shock schemes for approximating the level-set equation may be found in Osher & Sethian. Briefly, consider a one-dimensional version of the level-set equation for a front advancing with unit speed, and define H(YJ = (Y;)lI2. Then a forward time-discrete version of the equation may be written as (3.2) y n + 1 = Let g be an appropriate numerical flux function approximating H . Then we may directly approximate the spatial term and write where DZ (D;) is the forwards (backwards) difference operator. In multiple space dimensions and the special case where H(u) = u2, a particularly straightforward numerical flux function was given in Osher & Sethian, namely y n + l = !Pn+Atg(D;q,Diq), (3.4) This conservative monotone scheme is an upwind method, in that it differences in the direction of propagating characteristics. Equation (3.4) completely specifies the numerical approximation to (3.1). Details may be found in Osher & Sethian, and Sethian (1 989).
Addition of the curvature correction term to the flame speed S,(K) = Sz + S: (K) requires the separation of the two effects. We approximate the advection part, St with upwind differences as in the above, and the parabolic curvature term S~( K ) using central differences. Note that the curvature term K may be easily computed in this levelset setting by observing that the curvature is the divergence of the unit normal, and hence The above derivatives are easily approximated using central differences.
3.2. Velocity field due to exothermicity At each time step, !Pis updated to new values on an Eulerian grid by propagation and advection. A new flame position is obtained by passing field data Ti to a contouring routine and finding a set of flame segments corresponding to the level set Y = 0. The contouring routine works as follows. We locate the cell with four grid values of different signs, then we bisect the cell by the line separating the burnt from the unburnt region. The position and length of the flame segment within the cell determines the location and the strength of the volumetric source within the cell, m, which is given by where A1 is the elemental flame length. This volume source strength becomes the source term in the Poisson equation for the velocity potential computed on an Eulerian grid: where Q, is the velocity potential due to volume sources and f is the source strength. In order to use the Poisson solver, the volume source of strength m located in the middle of the flame segment within a cell must be distributed to the four nearest grid points with values Lj using an area weighting scheme.
The computational domain is a rectangle with no-flow boundary conditions across y = 0,1, and inflow-outflow conditions at x = 0,2. The flame propagates in the negative x-direction and volume is generated inside the burnt region owing to the density drop across the flame front. The boundary conditions for the Poisson equation should conform to the Neumann compatibility condition which requires that the total volume generated inside the computational domain matches the net outflow across the boundaries of the domain.
Depending on the boundary conditions prescribed, this expansion velocity field affects either both the incoming and outgoing flows or one of them. Since we want the incoming flow sufficiently far from the flame front to remain intact, the boundary condition is set such that the additional volume generated inside the domain goes out in the positive x-direction only.
The Poisson equation is solved on a discrete grid using a fast Poisson solver (see Swarztrauber 1974). The corresponding vector velocity field due to volumetric sources, Us, is computed by central differences on @ l ' j . The irrotational velocity field on the unburnt side of the flame is used for the vorticity generation in the rotational part of the velocity field.
The flame acts as a source of vorticity via the baroclinic torque term in the vorticity equation. The vorticity jump [w] across the flame front is given by (2.6).
The normal and tangential directions on the flame front are as indicated in figure 1 and the corresponding unit normal and tangential vectors are expressed in terms of

Velocity field due to baroclinicity y by
At the four grid points of the cell containing a flame segment we can compute the normal and tangential velocity components of the composite velocity field U (see (2.2)) : The normal and tangential velocity components at the midpoint of the flame segment within the cell are obtained by bilinear interpolation using the velocities at the four surrounding grid points. The absolute flame speed, V, , , is given by Once all the values along the flame front are calculated, the space derivative terms in (2.6) are approximated in the usual fashion. Special care is required for the time derivative term dUJdt, defined by The values of Y and U at time t -At are used to compute Us at the midpoint of a flame segment at time t-At. At time t, we obtain Us and V, at the midpoint of a flame segment using values of Y and U at time t. The first term on the right-hand side of (3.8) is computed using the values of Us at time steps t and t -At. To compute the second term on the right-hand side, we trace back a distance of V,At normal to the flame segment into the unburnt region to get Us and calculate the spatial derivative.
Flame-induced vortices are injected on the burnt side at each time step such that These vortices will occupy an area given by where AI is the length of the flame segment in a given cell, and the corresponding circulation can be expressed as 4 = w b S, A1 At, where S, is the relative velocity of the unburnt gas with respect to the flame front. This vorticity field located behind the flame on the burnt side becomes the source term in the Poisson equation for the vorticity stream function y?: V".. a? = -qj, (3.9) where the wij at the cell corners on an Eulerian grid are evaluated by area weighting the circulation contained within the cell. We use the vortex-in-cell method (Leonard 1980) to decrease the computation time. As in the case of the velocity potential due to volume sources, the Neumann compatibility condition should be satisfied such that the negative of the total circulation inside the domain is equal to the contour integral of the normal derivative of the vorticity stream function; given the impermeability condition along the boundaries in the y-direction of the computational domain, this normal component of the vorticity-induced velocity along those boundaries is cancelled by solving the Laplace equation V 2 0 i j = 0 with the boundary condition that the normal derivative of O along the impermeable boundaries is equal but opposite in sign to the normal velocity component of the vorticity-induced velocity. Also, the Neumann compatibility condition is satisfied by specifying the boundary conditions along the other boundaries such that the resulting contour integral of the normal derivative of the velocity potential O along the boundary should be zero, since the source term is zero in this case.
The addition of the two velocity components, one due to the vorticity stream function and the other due to the potential solver to accommodate the impermeability condition, yields the rotational velocity field due to the vorticity field. This rotational velocity combined with the volume-source-induced irrotational velocity and the uniform incident velocity specify the total velocity which moves the field data by advection.
We have the equation for flame propagation with advection as is given in (2.12),
A first-order upwind scheme is used for the advection part, U -V Y U l D ? qj+ UtTDT qj+ F+D? !J$+ 5-D: qj.

3.4.
Flame anchoring algorithm Power-production devices often require that a flame be retained in the system. Gas velocities in combustors usually exceed the laminar flame speed and hence the flames must be stabilized against blowout, a condition at which flames are transported through the exit of the burner and thus combustion ceases. In order for a point on a flame front to be stationary, the local condition for stabilization of the flame is that the normal velocity of the unburnt gas and the normal flame speed must be equal. This can be enforced by the presence of a flameholder or retention point (or line) or a small retention region.
Various objects may serve as flameholders: a continuous electrical discharge, a heated metal wire, a blunt body (so that hot combustion products lie in the recirculating region behind it), an auxiliary flame, or other means for continuous ignition of the fresh mixture. The ignition is transmitted to neighbouring portions of gas at the velocity component of the flow tangent to the flame.
To incorporate the above idea of flame stabilization into the level-set scheme, the flame is attached to the flameholder by laying down an initial ignition field Th on an Eulerian grid and letting lu,, act as a source of an ignition impulse. As the flame propagates and is advected by the accompanying flow field it is continually reignited at the flameholder by superposing 5, onto the existing Y. Further details can be found in Rhee (1992). We emphasize that our interest in the present work is in the dynamic far-field behaviour of open V-shaped flames, and that no attempt is made to model the detailed flow-field structure of any physical flameholder which would be needed to produce such a flame. The flame-retention algorithm described is used to ensure that our modelled flames, regardless of their initial assumed shapes, always remain within the computational domain and begin at the retention point (or, more accurately, at the retention cell).
3.5. Initial and boundary conditions The initial conditions for the level-set approach are straightforward. Given the initial position of the front, the function P i s initialized by using the signed distance function, that is where d is the minimum distance from a point x to the initial front, with the positive (negative) distance if the point is in the unburnt (burnt) region. In the case of a flameholder, we initialize a small set of cells around the flameholder; that is, we assume that fluid inside that small set of cells is burnt, and the remaining fluid is unburnt. The boundary conditions depend, of course, on the region under study. In the calculations that follow, we consider a flameholder placed between two walls. On the solid walls, the boundary conditions for Yare mirror conditions; this matches the fluid requirement of a Neumann boundary no flow condition. At the inlet, we simply choose a fixed boundary condition; the periodic reinitialization of the level-set function Y ensures that these upstream boundary conditions for the level-set function do not affect the downstream calculation. The downstream boundary conditions for Y are outflow conditions, and the upward differencing precludes needing boundary conditions along the boundary.
For the fluid mechanics, the boundary conditions are more subtle. As discussed earlier, a potential solver is used to construct the potential solution. On the sidewalls, the impermeability condition is satisfied. At the inlet, a uniform incoming velocity is prescribed. Downstream at the outlet, the outflow boundary condition is a cause for some concern. In Sethian & Ghoniem (1988) an extensive numerical study of downstream boundary conditions for vortex methods was performed, where flow over a backwards facing step was considered. The results of these calculations indicate that a uniform downstream boundary condition (Neumann condition) more than ten step heights downstream did not seriously affect the formation of vortices and vortex shedding near the step. Motivated by this, we choose a uniform downstream condition, although as will be seen, some computations were made with a ' top-hat ' exit velocity profile.
The choice of value for this uniform exit velocity is difficult. The presence of sidewalls accelerates the unburnt flow between the burnt flame portion and the walls themselves. Thus, by enforcing a uniform exit velocity, conservation of mass is slightly violated. We believe that this does not significantly affect the flame profile upstream. We note that for the physical unbounded flame which we model, there is no conservation-of-mass requirement.

Some Jinal numerical comments
There are several ways to improve the above algorithm. To begin, the use of a vortexin-cell method to compute the velocity field from the vortex-vortex interaction means that the vorticity distributions on a scale smaller than the grid size are necessarily lost; this adds numerical smoothing to the flow field which lowers the effective Reynolds number of the computed flow. The representation in the next Section of upstream turbulence by means of a distribution of vortex elements of constant strength is a very simple way to mimic the effects of turbulence; a more accurate way to capture the real statistical flow variations in a turbulent field might come from adjusting both the strengths and the core sizes of the upstream vortex elements. The representation of the flame volume generation on the corners of the grid cell smears somewhat the width of the reaction zone, although as the grid size goes to zero the discretization becomes more accurate. Finally, a more accurate calculation could be performed using more time steps and a larger computational domain by taking advantage of fast N-body solvers, adaptive meshes to track the flame zone more finely, and higher-order approximates in the Hamilton-Jacobi level-set solvers. Nonetheless, the above algorithm is a good first approximation to the essential physics.

Results
The numerical method for flame propagation and flow field description developed above is applied to the dynamic behaviour of a V-shaped open premixed flame in an incident turbulent flow field. Although we are simulating an open flame, numerical work requires a finite computational domain. Hence, the solution domain is truncated to a rectangle of axial (x) length equal to 2 and transverse ( y ) width equal to 1. Incoming fresh mixture enters the computational domain at x = 0 with velocity equal to 1. The transverse width which is set to 1 corresponds to the region from y = 0 to y = 1 in the computational domain. All velocities and lengths are scaled with incoming free-stream velocity U , and transverse width respectively; the combination of the two gives a time scale. The flameholder is located at (x, y ) = (0.5,0.5). Flame speed S z / U , and density ratio p u / p b are set to 0.08 and 6.0 respectively, unless otherwise specified. The Markstein length scale L in the flame-speed equation is set to 0.04. The grid size is equal to 0.02 and the time step is 0.004 which is based on the Courant condition; 500 time steps were used in computing statistical averages.
Both the kinematic and dynamic responses of the flame are investigated. The kinematic part is simply the case without exothermicity and baroclinicity. The dynamic part is composed of exothermicity and baroclinicity. Also the turbulent flow field of the incoming flow is simulated by injecting vortices at the domain entrance x = 0 whose transverse locations y between 0 and 1 are given in a random fashion. The development of flame cusps in response to a somewhat strong vortex pair is depicted. Various statistical data associated with the turbulence modelling are presented.

Flame response without exothermicity and baroclinicity
In this case the V-flame simply adjusts kinematically to the incoming flow. The sine of the angle between the flame and the centreline is equal to the laminar flame speed divided by the free-stream velocity. Since in the absence of exothermicity the incoming

Flame response with exothermicity and without baroclinicity
As the fresh mixture undergoes combustion, volume is generated owing to the temperature increase at constant pressure. There are several possible outflow boundary conditions to accommodate this volume production. The Neumann compatibility condition for the Poisson solver requires that the volume generated within the solution domain matches the integral of the net outflow across the boundaries. Since the two lateral boundaries y = 0 and y = 1 act as impermeable walls, this generated volume is prescribed along the exit boundary x = 2. If the flame does not encompass the whole width of the solution domain, the outflow boundary condition at x = 2 can be given in various forms as long as the Neumann compatibility condition is met. Figure 3(a) shows the flame response for the uniform outflow boundary condition at x = 2. Compared to the case for the same value of S l / U , but without exothermicity, the flame with exothermicity at equilibrium (i.e. after sufficiently many time steps) reaches out to the unburnt side due to the change of the velocity field in the unburnt side, because the reactants are deflected from the centreline toward the sidewalls and accelerated in the axial direction. The flame near the flame stabilizer reaches out to the unburnt side and the flame near the exit is mostly convected downstream by the dominant axial flow and hence is positioned almost in line with the flow velocity. It is perhaps worth noting here the danger inherent in attempting to determine experimentally the laminar flame speed S: from the observed flame angle, as has often been done. Without exothermicity, the equilibrium angle relationship O,, = sin-' (Sl/ U,) holds, but with exothermicity (the physical case) flow divergence and acceleration of the reactants invalidates this relationship, even for an unbounded Vflame (see below). For further discussion on this point, see Cheng, Shepherd & Talbot (1988). Figure 3(b) shows the results for the 'top-hat' boundary condition in which the excess volume manifests itself only on the burnt portion at x = 2. This boundary condition involves an abrupt jump of velocity across the flame location at x = 2. Since the Poisson equation is elliptic, the change in boundary condition is expected to influence the whole solution domain. But the flame location and velocity field up to x = 1.5 for this case is almost identical to case (a) for uniform exit boundary condition. From that location to the exit, the axial velocity on the unburnt side in case (b) decelerates and hence the flow converges towards the centreline to adjust to the 'tophat' boundary condition, causing the contraction of the burnt region of the flame at the exit. This unphysical flame response near the exit is strictly due to the boundary condition. Since our main interest lies away from the exit, the uniform boundary condition will be henceforth employed.
An increased lateral extent of the computational domain would reduce the acceleration of the reactants exterior to the flame, and would also allow increased excursions of the instantaneous flame interface due to turbulence interactions modelled by random discrete vortices, as described later. It would also make more viable a nonuniform exit boundary condition, which is more physically realistic than the uniformvelocity exit condition employed. However, the acceleration of the reactants exterior to the flame interface is not unphysical. It is observed in experiments, though not nearly as pronounced as that obtained from the computations. At the axial station x = 1.0 where we will make our comparisons between computed and experimental results, the experiments show at most about a 20 YO increase in reactant axial velocity, whereas the computations yield about a 60 % increase. This comparison evidently provides a clue to the lateral extent of the computational domain required for more faithful modelling, but owing to computational expense, such an inquiry was not pursued.
Although we are simulating an open flame at constant pressure, small pressure gradients in fact are present owing to exothermic flow deflection and acceleration, which we neglect except in the calculation of the flame-induced vorticity, since they are of second order as regards the flow field as a whole. In an open flame at essentially constant pressure, the far-field boundary condition due to volume expansion is determined at each time step by (3.6) summed over the instantaneous total length of the flame. The boundary conditions along the truncated domain of the numerical simulation are such that the volume generated due to exothermicity is prescribed at the exit (x = 2) and the impermeability condition is imposed at the lateral boundaries ( y = 0) and ( y = l), similar to the flow in a channel. In this confined computational domain, overall mass conservation would require that the mass flow at the inlet (x = 0) be equal to the mass flow at the exit (x = 2). However, for the cases considered the computed values of the outflow turned out to be greater by about 15 % than the inflow. Since we are modelling the far-field boundary condition of an open flame for which no mass conservation requirements pertain with boundary conditions for a finite domain, the check on the mass conservation though of interest computationally is in fact physically irrelevant.

Response to two strong vortices
Before addressing the more general question of how our model predits turbulent Vflame dynamics, including the effects of baroclinicity, it is of interest to analyse one of the constituents of the model, namely the interaction of a single line vortex pair with a flame interface. This problem has been studied experimentally by Namer et al. (1984) and Hertzberg, Namazian & Talbot (1984) in the context of a Karman vortex street interacting with a V-flame, and by Roberts & Driscoll (1991) who investigated the response of a flame sheet to a laminar vortex ring. The latter work was subsequently followed by a parametric numerical study of the interaction process (Wu & Driscoll 1992) using the SLIC algorithm, although the effects of curvature-dependent flame speed, volume generation and baroclinicity were not incorporated. Here we do not present such a parametric study, but rather just a single example which demonstrates how the level-set algorithm naturally reproduces the flame cusping phenomenon which is well-established experimentally.
The test case we have analysed is the release adjacent to the flameholder at x = 0.5 and t = 0 of two line vortices of non-dimensional circulation r = f0.2, the positive vortex being on the right and the negative one on the left. The value Irl = 0.2, with our chosen vortex core radius of 0.02, corresponds to a maximum azimuthal velocity-flame speed ratio ( Uoma,J/SE of 20, which places it according to the parametric study of Wu & Driscoll in the 'severely wrinkled' flame regime. Figure 4(a-d) shows the sequence of events as the vortex pair passes through the flame. In (a) the vortex pair is seen to have moved outward from its original release points ( y = 0.44,0.56 at x = 0.5) as it has been advected by the divergent reactant flow field, but at the same time there is a necking-in of the flame interface due to the strong circulation of the vortices. As the vortices pass through the flame in (b), (c) and (d), cusps develop, which later smooth out as the vortices traverse the product region behind the flame. The computations include volume production at the flame interface, which is reflected in the fact that the interface is curved, but do not include baroclinic vorticity generation. Other calculations have been carried out which incorporate baroclinicity, but the results are not qualitatively different from those presented, since the baroclinic vorticity is considerably weaker than the vorticity associated with the interacting vortex pair. While this single computation does not provide the extent of information given by the investigation of Wu & Driscoll, it does demonstrate the capability of the level-set algorithm to deal effectively with vortex-flame interaction, including exothermicity and curvature-dependent flame speed, which is an essential ingredient in the modelling of the dynamical behaviour of turbulent flames.

The complete model including exothermicity, baroclinicity
and free-stream turbulence The modelling of the effect of exothermicity on the dynamics of the flame interface has been described in $4.2. We now add to our model the additional effects on the flame dynamics of turbulence in the oncoming reactant flow and baroclinic vorticity generation at the flame interface. Additionally, computations were carried out in which free-stream turbulence and exothermicity were incorporated, but baroclinicity omitted. In the interest of conciseness these results, although referred to, will not be given here. Complete details can be found in Rhee (1992).
To simulate computationally the effect of free-stream turbulence, vortices of circulation r = k 0.01 are injected at x = 0, this value of r having been found to yield a free-stream turbulence level of about 8 %. The transverse ( y ) locations of the injected vortices are distributed between y = 0 and 1 using a random number generator. The same number of positive and negative vortices are injected, so that the total free-stream circulation remains zero at all times. The time step interval between successive injections of the vortices is chosen such that the mean transverse and axial distances between adjacent vortices are almost equal, as required for the simulation of isotropic turbulence.
The injection of vortex elements to mimic free-stream turbulence is delicate. In reality, the problem is three-dimensional, in which vortex stretching would contribute to turbulence decay and cross-correlation mixing effects (see Sethian & Ghoniem). Thus, our assumption of two-dimensional flow is a significant restriction in our analysis of the various competing factors of exothermicity, baroclinicity, and freestream turbulence in the flame motion and wrinkling. Some features of turbulence decay could have been incorporated in the model by adding a random walk component to the vortex displacement at each time step, but this was not done. In future calculations, we hope to construct a three-dimensional and improved version of the algorithm described here.
To model the flame-generated vorticity, the formulation by Hayes (1959) which is composed of steady and unsteady terms is used to give the strengths of the vortices. Since the unsteady term, which was omitted in the simulation by Pindera & Talbot (1986), is calculated to be of the same order of magnitude as the steady terms, both are included in the computation. The baroclinic vortices are injected on the burnt side of the flame at each time step. They are advected by the local flow field and modify the flame motion and the velocity fields of both reactants and products.
In making comparisons between computed and experimental results, it is desirable to discriminate between the passage of unburnt and burnt fluid at a point by means of conditional sampling. To this end, we define an instantaneous progress variable c(x, t ) such that c is zero in the reactants and unity in the products. While in the limit of a vanishingly thin flame interface, c can have only 0 , l values, its meadvalue (c(x)) can range between zero and unity, depending on the fraction of time the observation point is occupied by either reactants or products. This mean value, also termed the intermittency factor, which is the computations is obtained by sampling the time record of a variable and evaluating N , / N where N , is the number of samples associated with product fluid and N the total number of samples, forms the basis for conditional sampling. For example, the unconditional Eulerian mean axial velocity at a point is given by u = (1 -(c(x))) u,+ (C(X)> where U, and U p are reactant and product axial velocities respectively. Plots of (c(x)), which do not reveal any exceptional behaviour, are contained in Rhee (1992). Figure 5 shows the superposition of instantaneous flame shapes at successive time intervals. The brush-like structure of the flame region is evident. The individual flames were found to be generally smooth, with cusping rarely observed. Flame excursions were somewhat greater than what was found in computations in which baroclinic vorticity production was omitted, but the flame shapes were similar. The mean flame included angle for 0.5 d x d 1.0 is about 25". The velocity field and vortex distribution of one instantaneous flame at t = 1.4 are shown in figure 6. One can observe from the velocity field, figure 6(a), that there is an axial acceleration of the flow in both the reactants and the products and that in the products this acceleration is most pronounced at the centreline y = 0.5. We will return to this point later. Figure 6(b) shows the distribution of both incident and flamegenerated vortices within the flow field for this instantaneous flame. Vortices of positive and negative circulations are represented by filled and open circles, respectively. About 8500 flame-generated vortices are present at this time step (and a comparable number in the reactant flow) so for plotting purposes they are combined on an Eulerian grid. Of most interest is the predominance of negative vorticity in the products behind the right-hand (lower) flame sheet, and of positive vorticity behind the left-hand (upper) sheet, as is required from the sign of the baroclinic term (l/p2) V p x VP. The majority of the vortices in the product region are baroclinic, as computations show that when baroclinicity is omitted the vorticity concentration in the products associated with passage of free-stream vortices through the flame is reduced owing to the dilatation of the flow. Note that in the coordinate system employed, which adopts the same sign convention used by Ashurst, the right-hand sheet of the flame occupies the region 0 < y d 0.5, and the left-hand sheet the region 0.5 < y < 1.0.
Positive values of V are from right to left. This sign convention is opposite to that employed by Cheng (1984) in presenting his experimental results. figure 9. The acceleration of the flow in the products near the centreline is evident. This acceleration is a baroclinic effect, since it was not observed in the absence of flamegenerated vortices. The conditioned-average transverse velocity at x = 1 .O is shown in figure 10. As seen, the mean V velocity is away from the flame in the reactants, and towards the centreline in the products, owing to the reversal of flow direction across the flame brush.
Values of unconditioned and conditioned kinematic Reynolds stress -UV are plotted in figures l l ( a ) and l l ( b ) (u and u are the fluctuating components of U and V , respectively). While the distributions are quite noisy, one feature of interest is that the reversal spikes at the flame locations y = 0.3,0.7 present in the unconditioned data are absent in the conditioned data. This again is an intermittency effect. Notice also that the Reynolds stress within the products is clearly positive on the left-hand side and negative on the right-hand side. This is a manifestation of the predominance of positive and negative flame-generated vorticity in the products on the left-and right-hand sides respectively, since the u-fluctuations have no preferred sign. A comment is in order regarding the noise level in the statistics of the fluctuating velocity components. Since we were mainly interested in the baroclinic effect, we used a large number of vortices in the flow-field computations, as was noted earlier. However, because of limitations in computational resources, we employed only 500 time steps in computing our statistical averages. By way of comparison Ashurst employed only 400 vortices in his simulation of the turbulent flow field, but averaged over 6000 time steps, and consequently his averages are far smoother. The experimental data of Cheng (1984), with which comparisons will be made, are likewise less noisy than our computed results, because each data point represents the average of 4096 validated velocity data pairs. Clearly, we could have reduced the noise level in our statistics by increasing significantly the number of time steps used in the computations, but this would have become prohibitively expensive, since even with the limited number used a typical run required about two hours of Cray I1 CPU time. Evidently, a major objective of future studies will be to optimize the apportionment of computational resources between vortex population size and number of time steps, as well as the required size of the computational domain which was discussed earlier with respect to the outflow boundary conditions. The probability distribution of values of flame curvature is shown in figure 12. The

Discussion and comparison with experiment
In a typical experimental premixed turbulent V-flame configuration (Cheng 1984;Cheng & Shepherd 1984;Cheng et al. 1988) which is suitable for comparison with this numerical simulation, an inner core of fuel/air mixture 5.0 cm in diameter is surrounded by an outer co-flow of air to reduce shear layer effects in the reactants. The turbulent flame is stabilized by a 1.0 mm diameter rod placed at the exit of the flow nozzle. The incident free-stream turbulence is generated by a square mesh grid or a perforated plate. The free-stream velocity is about 5 m s-l, the turbulence level is between 5 and 8 % of the free stream and the laminar flame speed is 40 to 50 cm s-'. Flame statistics such as turbulence intensity, Reynolds stress, length scales, and flame brush thickness are reported, and provide the means for direct comparison with the present numerical simulation.
The values of St, pJpb and the strength, number and frequency of injected vortices for the numerical simulation of the incident-flow turbulence level were specifically chosen with the intention to compare with the above experimental results. Since we are simulating an open flame and do not model the detailed features of the flameholder region, there is no intrinsic length scale present. One way of establishing a scaling relationship is to relate the computed integral length scale to the physical integral length scale. The reactant integral length scale was calculated to be 1.2 times the grid size. This is the transverse length scale I, . The measured integral scale was obtained to be about 3 mm (Cheng et al., 1988). This is, however, the streamwise length scale, l,, determined from a single-point time series using the Taylor hypothesis. Since the freestream turbulence in the numerical simulation was very nearly isotropic, it is assumed that the condition of 1, = 21, applies, and this allows us to relate the numerical scale of the computations to the physical scale of the experiments. This means that the grid size of 0.02 in the computations corresponds to about 1.2 mm in physical scale. All numerical results can then be converted to the physical scale using this relation for direct comparison with experimental results. Most of statistical data in the product region of this work were obtained at x = 1 .O which corresponds to 30 mm downstream of the flameholder in the physical scale and these will be compared with the experimental results reported at the location 50 mm downstream of the flame-holder in the work of Cheng (1984).
We first consider the computed flame response to two strong vortices, in comparison with the experimental observations of Hertzberg et al. (1984) on the interaction of a Karman vortex street with a V-flame. Flame cusping was observed in the experiment, though not as strong as was obtained in the numerical simulation. The flame configuration produced in the cusping process clearly depends on the ratio of the maximum vortex-induced velocity to the free-stream velocity as the vortex interacts with the flame. For a Karman vortex street at rod-generator Re = 240 of the experimental conditions, this ratio is in the vicinity of 0.3 (Blake 1986), whereas this value was about 1.6 in the numerical simulation. This rough estimate suggests that the cusping observed in the experiment might be expected to be less severe than that found in the numerical simulation, although it is clear that they share the same general structure. For the weak incident turbulence level chosen for the complete simulation, essentially no flame cusping was observed.
The mean velocities in the flame region of the complete simulation show that the products move faster than the reactants and that flow direction is towards the centreline for products and is away from the centreline towards the sidewalls for reactants, which is also observed experimentally. There is a significant increase in the centreline axial velocity owing to flame vorticity as compared with the case of exothermicity only, a characteristic of the flow field clearly evident in Cheng's experimental results, which are plotted in figure 9 in terms of the product velocity distribution relative to the adjacent reactant velocity at the same station. This effect was also found by Pindera & Talbot and is a striking feature of the effect of baroclinicity. The conditioned values of U,,,, V,,, and the Reynolds stress in the burnt region with baroclinicity included show an overall increase, compared with exothermicity only. These conditioned data indicate that the flame-generated vorticity plays an important part in determining the turbulent flow field. The flame-generated vortices are responsible for the increase in the conditioned turbulent intensities in the products compared to the case without flame-generated vortices and are therefore considered to be the true source for the so-called flame-generated turbulence since the conditioned results do not include the effects of intermittent flame motions.
This conclusion is supported by Ashurst's (1987) results, in which no increase of the r.m.s. fluctuations in the products was found. Ashurst notes that the inclusion of baroclinic vorticity production would probably have produced better agreement with experiment. In the experimental results of Cheng (1984), the conditioned values of U,,,, V,,, within the flame brush are about two times higher than those in the reactant region. But these values within the flame brush decay to the level in the reactant in the post flame region. Since the turbulence cascade and viscous dissipation are not included in the numerical model, this decay of turbulence level was not observed in the numerical simulation. But we can at least estimate how much the turbulence level decays in the product region. According to the experimental results by Batchelor & Townsend (1948) on the decay of isotropic turbulence in the initial period, about a 50 % decay in the turbulence level from x = 1 .O to 2.0 (which is equivalent to 60 mm in distance) is estimated. This suggests that after the turbulence level within the flame brush region is increased by flame-generated vortices, viscous dissipation takes over within the product region to reduce the turbulence level significantly. The kinematic Reynolds stress within the flame brush in the experiment had a maximum absolute value of about 0.08 (m s-')'. The absolute magnitude of the Reynolds stress in the product region is calculated to be roughly about 0.004 which is equal to 0.1 (m s-l)' in physical dimension, in reasonable agreement with experiment.
Two experimental values of flame brush thickness reported by Cheng are plotted in figure 13, scaled to non-dimensional units. It is seen that agreement with the computational results is good. In the work of Pindera & Talbot (1986), the flame brush thickness including baroclinicity was found to be slightly smaller than that without baroclinicity, whereas the opposite is true according to our simulation. While the reason for this is unclear, it may be associated with the neglect by Pindera & Talbot of the unsteady term in the vorticity jump relation, the relatively sparse number of flame vortices (about 400) employed by them as compared with the present computations which involved as many as 8500, and possibly most importantly, the superiority of the level-set algorithm over the marker particle method employed by Pindera & Talbot to track the flame motion. A qualitative comparison between a computed flame interface and one observed experimentally is displayed in DOE Lawrence Berkeley Laboratory Report LBL 33445, UC 350, July 1993. The comparison reveals that the experimentally observed flame excursions are considerably greater than those computed, undoubtedly in part because of the constraint of the finite lateral extent of the computational domain, and the rather coarse grid used in the computations. Both of these limitations could be relaxed in future work, at the expense of increased computational cost. However, despite these limitations, it is of interest to note that the computed mean flame included angle of 25" compares quite well with Cheng's experimental value of 28". In contrast, Ashurst obtained a flame angle of 65", and Pindera & Talbot a value of 45" in their unbounded flow calculations. While the reason for the discrepancies is unclear, one conjecture is that it may be related to the fact that the latter two works employed two-sided volume sources, while the present computation used an unidirectional formulation.
The distance between flame crossings at maximum probability is found to be about 0.2 which in physical dimension means that the instantaneous flames meet with the average location of the flame at 12 mm intervals at maximum probability. This translates into a mean flame crossing frequency of about 500Hz whereas the experimental value was found to be about 600Hz (Cheng et al. 1988), in good agreement. This agreement between computed and measured flame crossing frequencies further implies that the flame wrinkle scales produced in our simulation are consistent with those observed experimentally.
The probability distribution of curvature along the flame front compares well with experimental results (Shepherd & Ashurst 1992) although the experiments were done on a stagnation flame. Positive and negative curvatures are about equally present with slight negative skewness for a weakly turbulent flame.

Conclusions
The level-set algorithm, with volume generation and baroclinic vorticity included, provides a simulation of turbulent premixed flame dynamics in reasonably good agreement with experimental results for input parameters chosen to agree with experimental conditions. No adjustable parameters other than the Markstein length scale are involved. Some of the important features predicted for the case investigated are: (a) included flame angle, (6) flame brush thickness, (c) velocity r.m.s. and Reynolds stress levels, ( d ) product axial velocity acceleration, (e) flame crossing frequency, ( f ) flame curvature distribution.
Flame-generated vorticity is found to have an important effect in the prediction of velocity r.m.s. values, Reynolds stresses and particularly the axial flow acceleration in the product region behind the flame which is observed experimentally. The unsteady and steady terms in the vorticity jump relationship contribute about equally to the total flame-generated vorticity.
The level-set algorithm is shown to predict flame cusping for a situation in which a strong vortex interacts with the flame front, in accord with experimental observations of the passage of a Karman vortex street through a flame. In addition, for the relatively low turbulence level chosen for the full numerical simulations, cusping was rarely found, and the statistics of flame curvature exhibit essentially a symmetrical distribution about a zero mean, similar to what has been observed experimentally for a stagnation-point flame. Volume generation at the flame front is shown to play a decisive role in determining the flame angle, and it is shown clearly that the experimental evaluation of the 'turbulent burning velocity' through observation of the flame angle leads to an erroneous result, if the outward deflection and acceleration of the incoming flow are not taken into account, as has also been established experimentally by Cheng & Shepherd (1986). The authors would like to thank Drs Robert Cheng and Ian Shepherd for their valuable discussions and suggestions.