- Email: [email protected]
Simulation of a turbulent premixed open V-shaped flame using contour advection with surgery
Applied Mathematical Modelling 30 (2006) 1103–1117 www.elsevier.com/locate/apm
Simulation of a turbulent premixed open V-shaped flame using contour advection with surgery B.H.Y. Tang
a,*
, C.K. Chan a, J.S.L. Lam
b
a
b
Department of Applied Mathematics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong Environmental Protection Department, The Government of Hong Kong Administrative Region, 24/F Southern Centre, 130 Hennessy Road, Wanchai, Hong Kong Received 17 January 2005; received in revised form 24 June 2005; accepted 24 August 2005 Available online 25 October 2005
Abstract Despite its capability of high spatial resolution, simulation of turbulent flows with traditional Lagrangian (front tracking) scheme is often discouraged by numerical instability caused by clustering of marker nodes and topological changes of fronts. Contour advection surgery (CAS), being a robust front tracking scheme, can limit the growth of front complexity during simulation without jeopardizing accuracy or efficiency. This allows it to open up an advantage over traditional front-tracking schemes. It has already been demonstrated that CAS, with incorporation of the reaction sheet model, can accurately simulate the propagation and advection of a turbulent premixed V-shaped flame. In this study, it is further tested with 10 values of vortex circulation. A range of upstream turbulence levels of 1.8–19.8% was obtained. Results indicate that upstream turbulence increase the average flame length, flame zone area and the overall burning rate. Flame surface density R was also estimated. Maximum values of R obtained lie in the range 0.1–1.4 mm1. Skewness towards the burnt region was observed in all profiles of R. Similar to results from laboratory experiments, it was found that R values decreases with upstream turbulence. From this study, the ability of CAS to cope with intense turbulence is demonstrated and a better quantitative understanding on the scheme has also been acquired. 2005 Elsevier Inc. All rights reserved. Keywords: Contour advection with surgery; V-shaped flame; Flame surface density; Burning rate
1. Introduction Turbulent premixed combustion is widely used in a range of engineering devices and has attracted many research interests over the years. In the field of turbulent premixed combustion, of primary research interest is the accurate simulation of a premixed flame propagating in an ambient turbulent flow. Due to the non-linear coupling of mechanical turbulence and combustion process, devising an accurate numerical method which is capable of capturing details of the propagation without masking the nature of physics, still remains a challenge today. *
Corresponding author. E-mail address: [email protected] (B.H.Y. Tang).
0307-904X/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2005.08.016
1104
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
In this paper, a two-dimensional rod-anchored turbulent V-shaped premixed flame is considered. This is one of the most common configurations employed for studying turbulent flame in laboratory experiments. Amongst the many theoretical models developed for simulating such flame numerically, the reaction sheet model [1] appears to be the most frequently invoked. This model allows the internal structure of the flame to be neglected by using the fact that, when the reaction rate is high, the flame front can be approximated as an infinitesimally thin boundary separating burnt and unburnt regions. With the exclusion of the flame frontÕs internal structure, the geometry of the flame front becomes the sole element that governs the evolution of the flow field. In turbulent premixed flames, two common occurrences namely, the topological changes (merger and breaking of the flame front) and the development of cusps, are the major influences on the geometry of the flame front. There are in general two approaches for capturing the flame front movement accurately. The first one is called the front-capturing approach which was proved to be very successful in numerous works [2–4]. However, being Eulerian in nature, the spatial resolution is somehow limited by the grid size of the computational domain. Thus, any sub-grid features would not be resolved by this grid-based approach. The second one, termed the front-tracking approach, being Lagrangian in nature, does not have such spatial resolution limitation. This approach has received some attention [5,6]. However, in general applications, it remains to be less popular than the Eulerian approach. This is because numerical instability is often caused by the two common occurrences in turbulent premixed flames mentioned earlier. Clearly, in order to exploit its advantage of high spatial resolution and thus capturing the sub-grid features, this numerical instability problem must first be solved. In view of this, Lam et al. [7] employed a front-tracking scheme known as contour advection with surgery to treat the topological changes of flame fronts and development of cusps. This front-tracking method originates from a numerical technique called contour dynamics (hereinafter CD) meaning the dynamics of boundaries of vorticity discontinuity (i.e. contours). Contour dynamics is commonly used in simulating inviscid fluid motions, particularly in vortex motion. It is based on the ‘‘water-bag’’ model for solving the Vlasov equation that appears in plasma physics [8] and was first introduced by Zabusky [9] in his simulation of two-dimensional fluid motions governed by the barotropic vorticity equation in a domain of infinite extent. One major drawback of CD that prevents its further application is the inability to deal with the rapid increase in the number of nodes used to represent the contours. To maintain a suitable spatial resolution, it is often necessary to continually insert extra nodes at regions of fast-growing curvature. This inevitably slows down the computation or in the worse case scenario, brings the calculation to a virtual standstill. Another drawback is the unmanageable growth of contour perimeters. Due to continual evolution of the contours, the contour perimeters will grow in time and very often resulting in a flow too complex for CD to resolve. In view of this, Dritschel [10,11] extended CD to ‘‘Contour Surgery’’ (hereinafter CS). It was introduced to especially treat the development of fine-scale features and thus limiting the growth in the number of nodes in a consistent manner. This is accomplished by truncating features that are formed under a predefined scale, d. It also topologically reconnects contours/divide a contour if the contours/parts of a contour become closer than d. The front-tracking scheme CAS was developed [12] after the success of CS and it was mainly used in geophysical research both as a numerical tool to advect fronts and a diagnostic tool to investigate fine-scale features [13–15]. Recently Lam et al. [7] have adopted the concept of CAS and applied it to a two-dimensional turbulent premixed open V-flame. Their computed velocity statistics and estimated flame surface density are in good agreement with corresponding laboratory measurements [16]. The present work, as well as being part of a series of extensive tests on CAS, is a continuation of the work done by Lam et al. [7] with the primary objective of acquiring a better quantitative understanding of this method. By varying the circulation of the vortices injected at the computational domain entrance, a wide range of upstream turbulence levels was obtained. Analysis was made on the resulting average flame length, flame area, flame surface density R and burning rate. 2. Basic equations The main assumptions of the model are: (i) the reaction is high so that the flame front is infinitesimally thin,
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1105
(ii) burnt and unburnt regions have distinct uniform densities, (iii) Mach number is very small so that the flow can be regarded as incompressible on either side of the flame, (iv) the mechanism of vorticity production is inviscid. At low Mach number, the velocity field, u for the combustion process can be decomposed into three components namely us the solenoidal component due to volume expansion across the flame front, ur the rotational component due to vorticity distribution x(x) and up the potential velocity field u ¼ us þ u r þ u p .
ð1Þ
The following three equations describe the conditions which have to be satisfied by individual components: r us ¼ mdðx xflm Þ;
r us ¼ 0;
ð2Þ
where m is the volume source strength along the flame front, x is the position vector, with subscript flm denoting the flame position and d( ) is the two-dimensional Dirac delta function r ur ¼ 0; up ¼ r/;
r ur ¼ xðxÞ;
ð3Þ ð4Þ
where / is the velocity potential of the incident flow. An expression for m is now derived. Let Su and Sb be the relative normal fluid velocities with respect to the flame front on the unburnt and burnt sides, respectively. The volume source strength is the net volume flux generated in the burnt products and consumed in the reactants, then m ¼ Sb Su.
ð5Þ
Denoting qu and qb as the fluid densities of the unburnt and burnt regions, respectively, the mass conservation law yields qu S u ¼ qb S b .
ð6Þ
From Eqs. (5) and (6), m can be expressed as m ¼ ðqu =qb 1ÞS u .
ð7Þ
According to Markstein [17], Su varies with curvature in such a way that a perturbed flame is stabilized as long as the curvature along the flame remains weak. Hence for weak curvature, Su can be approximated as S u ¼ S L ð1 þ ejÞ;
ð8Þ
where e is the Markstein length scale and j is the local curvature which is taken as positive/negative when the centre of the circular arc lies left/right of the flame front. Since Eq. (8) does not involve the use of any correction term, the effect of strain on the flame propagation velocity is not modelled. According to Ashurst [18], such exclusion of strain term is possible. He modelled a turbulent V-flame under conditions similar to some of that employed in the present study with the full stretch term and found that the velocity statistics did not show any dependence on the strain term. In the present simulation, the vorticity x(x) has two sources, namely the prescribed upstream turbulence and the flame front itself. The vorticity generated from the latter source is also known as flame-induced vorticity. Pressure gradients tangential to the flame result in different accelerations in the unburnt and burnt sides with distinct densities. Vorticity is thus generated at the flame by the mean density gradient across the flame and the pressure gradient tangential to the flame. The upstream source is simulated by injecting small uniform circular vortices at the domain entrance. To account for the flame-induced vorticity, the flame front is divided into small segments of equal length Ds with one small uniform circular vortex being introduced immediatelypbehind the mid-point of each segment in each time step Dt. The radius of this flame-induced vortex is given by (SuDsDt/p) so as to conserve the volume flux across the flame front in each time step. The determination of the flame-induced vorticity involves the use of an expression proposed by Hayes [19] for the vorticity jump across the flame front 1 1 q qu ½x ¼ fDut þ ut rt ut ut un j un rt un g; ð9Þ rt ðqu S u Þ þ b qb qu qu S u
1106
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
where ut and un are respectively the relative tangential and absolute normal components with respect to the flame front, $t is the gradient along the flame front and D is the material time derivative taken at a point which always lies on the flame front and moves in a direction normal to the discontinuity as it moves. Both the turbulence vortices and the flame-induced vortices are treated by the Random Vortex Method [20]. It can simulate the diffusive effect of viscosity on vorticity distribution and prevent the clumping of vortices about the flame front. 3. Numerical algorithm—CAS Since CAS was not originally devised for turbulent reacting flow, it is necessary to make modifications from the original. Following Lam et al. [7], the essential details of steps involved in CAS are presented in the following subsections. 3.1. Geometric configuration of front As with typical Lagrangian schemes, the front is discretized into an orderly set of connected marker nodes {xi} as illustrated in Fig. 1. Arrows on the front signify its orientation. In the present study, the flame front is oriented such that the unburnt region (reactants) always lies on the left while the burnt region (fully burnt products) always lies on the right. The displacement between two consecutive marker nodes xi and xi+1 is denoted by ei, i.e. ei ðai ; bi Þ ¼ xiþ1 xi .
ð10Þ
For ei, an outward normal vector fi is defined as fi ¼ ðbi ; ai Þ;
ð11Þ
which points towards the reactant as illustrated in Fig. 1. At any marker node xi, the continuity of tangent slope, ti is ensured by using the approximation ei1 ei ti ¼ þ 2. 2 jei1 j jei j Similar to ei, a corresponding outward normal ni can be defined for ti ^ ti ; ni ¼ k jti j
Flame front
right
left UNBURNT ti
xi+1 ei
fi
BURNT xi ni ei-1 xi-1
Fig. 1. Geometric configuration of a discretized front {xi}.
ð12Þ
ð13Þ
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1107
ti ni ti xi
ni
xi (a) κ i > 0
(b) κ i < 0
Fig. 2. Sign convention of curvature ji.
^ is the upward unit vector normal to the plane of the flow domain. In general, the two orthogonal pairs where k {ei, fi} and {ti, ni} are not identical. At any marker node xi, the curvature is determined by fitting an arc through its preceding and following nodes i.e. nodes xi1, xi and xi+1. It is then given as ji ¼
2ðai1 bi bi1 ai Þ 2
jti jjei j jei1 j
.
ð14Þ
Evidently, the curvature defined in Eq. (14) is signed. As illustrated in Fig. 2, ji takes a positive/negative value if the centre of the arc lies to the left/right of the front. 3.2. Contour surgery In every time step, each of the marker node is advected using the fourth-order Runge–Kutta scheme. Once all the marker nodes are advected and propagated to a new position, contour surgery is performed if necessary. At the end of each time step, the distance between any two non-adjacent marker nodes is computed. If this distance is less than the predefined threshold value, d (note that d no longer represents the delta Dirac function as in the preceding chapter), a so-called surgery will be carried out. There are three kinds of CS, namely fission, fusion and shorten surgery. The decision of which type of surgery to perform is based on the positions of the two non-adjacent marker nodes. If they are along the same front, a fission surgery is performed and the front is separated into two disconnected fronts as illustrated in Fig. 3(a). If the two non-adjacent marker nodes are on two disconnected front, a fusion surgery is performed as depicted in Fig. 3(b). In the simulation of flame front propagation, overshooting can cause an undesirable effect of unburnt reactants appearing as burnt product. Both types of aforementioned CS are effective in preventing overshooting of fronts and the associated numerical errors. In addition, even when overshooting does not occur, it is possible for fronts/parts of a front to coincide if the distance between them is small enough. This results in a redundant and dynamically ineffective front, as illustrated in Fig. 4. With surgery, this is avoided because any two fronts/ parts of a front which have a distance less than d between them is treated immediately. There exists a special type of fission surgery namely, shorten surgery. When the fission surgery is invoked, two disconnected fronts are form. If one of these disconnected fronts has a spatial extent that is less than the threshold scale d, it is eradicated as depicted in Fig. 3(c) and this is the shorten surgery. With this surgery, the problem of numerical instability triggered by the clustering of marker nodes at regions of high curvature can be resolved. Both fission and shorten surgeries are in fact the so-called Ôde-loopingÕ techniques [21] which remove certain marker nodes on the front so as to prevent the formation any needless loops (filamentary cusps). 3.3. Redistribution of marker nodes To maintain a suitable resolution and prevent clustering of marker nodes at a region of high curvature, it is essential to redistribute the marker nodes after surgery. Before the redistribution, each front is first split into
1108
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
Fig. 3. Three types of contour surgery (a) fission surgery, (b) fusion surgery and (c) shorten surgery.
Fig. 4. Schematic evolution of fronts following Fig. 3(a) and (b) without (left) and with (right) contour surgery.
curved segments delimited by extremities and/or corners. A marker node xi is labelled as a corner if the displacement vector ei and ei1 form an acute angle. For each curved segment, an appropriate average node density, ki, between each pair of adjacent marker nodes is calculated as follows: 1 1 ~ 1 ~ ~i þ ; j ki ¼ min ; ð15Þ d 2l L ~ ~ ~iþ1 ~i þ j j ~ ~ ~i ¼ ; ð16Þ j 2 ð~ ji1 =jei1 jÞ þ ð~ ji =jei jÞ ~ ~i ¼ j ; ð17Þ ð1=jei1 jÞ þ ð1=jei jÞ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ji þ jiþ1 2 ~i ¼ ; ð18Þ j þ 2 L2 where L is the characteristic length scale which depends on the characteristic radius of curvature of the initial flame front and l 6 1 is a positive non-dimensional input parameter for overall density of marker nodes. The physical significance of these two parameters will be discussed later in this sub-section. Examining Eq. (15) more closely, one can deduce that ki is an increasing function of ji and must be less than or equal to 1/d. The former condition ensures adequate resolution by placing more marker nodes at region of high curvature while the latter condition prevents clustering of marker nodes by limiting the minimum distance
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1109
between any two adjacent marker nodes to nearly d. Eqs. (16)–(18) give immediate values of averaged curvatures (weighted or not) and are all positive (P1/L). As commented by Dritschel [10,11], to determine the optimal placement of nodes, one requires more than the local instantaneous curvature of the contour. This is because any approaching small-scale features will not be detected with the use of a local method of node adjustment. Hence Eqs. (16)–(18) are formulated such that they depend on not only the relative positions of three consecutive marker nodes but also on the geometry of a broader neighbourhood. To explain the physical significance of L and d as mentioned previously, one may begin with the minimum density of marker nodes. Clearly, minimum density of marker nodes is obtained with a straight front. From Eqs. (14) and (16)–(18), the local and non-local curvatures are respectively zero and 1/L. Note that the non-local curvature is non-zero even for a straight front. In this way, marker nodes are ÔreservedÕ and further addition is possible as the straight front evolves. Rewritting Eq. (15) for a straight front, the minimum density of marker nodes is 1 1 ki ¼ min ; . ð19Þ d lL Since 1/d is the upper limit for ki, one must select L and d such that d < lL.
ð20Þ
In fact, the relation described by Eq. (20) holds not only for a straight but also for any curved front. This is because for a curved front, ki is bounded as shown below 1 1 6 ki 6 . lL d
ð21Þ
Returning to the straight front case, setting the non-dimensional parameter l as unity i.e. its maximum value, would give the minimum density of marker nodes as 1 ki ¼ . ð22Þ L This means that at least one marker node must be inserted per arc-length of L. Decreasing the value of l simply means that the number of marker nodes per arc-length of L multiplies. It follows that l is in fact an inverse magnification factor that reserves extra marker nodes for the subsequent evolution of the turbulent flame front. After redistributing marker nodes with the results from Eqs. (15)–(18) for each curved segments, interpolation of marker nodes is carried out. Following Lam et al. [7], consecutive marker nodes xi and xi+1 are interpolated using piecewise cubic spline as follows: xðsÞ ¼ xi þ sei þ gðsÞf i ;
ð23Þ
where s is the fractional distance between marker nodes xi and xi+1 and gðsÞ ¼ ai s þ bi s2 þ ci s3
ð24Þ
whose coefficients are given by 1 ai ¼ ð2ji þ jiþ1 Þjei j; 6 1 bi ¼ ji jei j; 2 1 ci ¼ ðjiþ1 ji Þjei j. 6
ð25Þ ð26Þ ð27Þ
4. Applying CAS to a turbulent premixed V-shaped flame 4.1. Experimental set-up Similar to Chan et al. [2,3], the numerical conditions in the present study are based on a series of V-flame laboratory experiments performed by Cheng [22–25]. Fig. 5 depicts the experimental set-up employed in
1110
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117 Turbulent Flame 150 mm
Su
U
Flameholder
100
V 50
Turbulence
x
Generator
y 0
Air
C 2H 4
Air
+Air
Fig. 5. Experimental set-up for investigating a rod-stabilized V-shaped flame.
ChengÕs investigation. It consists of a coaxial nozzle providing a circular jet with an inner core of ethene and air. This premixed fuel is shielded from the ambient environment with an outer layer of air. A turbulence grid is placed 50 mm upstream of the nozzle exit. When the premixed fuel passes the grid, turbulence is generated within the flow. A circular rod with a diameter of 1 mm is placed at the nozzle exit. When the premixed fuel is ignited, a turbulent V-shaped flame anchoring at the circular rod would be observed. This turbulent V-shaped flame travels with a normal velocity Su with respect to the reactants. While the V-flame problem is threedimensional in nature, the third velocity component in the z-direction is ignored in ChengÕs experiments. Thus it is reasonable to assume the flow is essentially two-dimensional. 4.2. Computational domain With CAS, the computational domain is unbounded. However, for a proper treatment of the flame front, one must define a spatial extent because an infinite flame cannot be achieved physically or computationally. Based on the size of the V-flame obtained in ChengÕs experiment, the solution domain is truncated to a rectangle of axial length (x-direction) equal to 150 mm and transverse width (y-direction) equal to 120 mm. For convenience, the present problem is transformed into a non-dimensional one. The non-dimensionalising length scale is selected to match the diameter of the inner coaxial cylinder of the nozzle. While the non-dimensionalising time scale is chosen such that the non-dimensional inflow velocity, U0 can be scaled to unity. Thus the following length and time scales are chosen for scaling all parameters: L ¼ 50 mm; 1 s. T ¼ 110
ð28Þ ð29Þ
Consequently, the computational domain now becomes non-dimensional with axial length, Lx equal to 3 and transverse width, Ly equal to 1.2. This is illustrated in Fig. 6. There is one boundary condition imposed at the lower boundary of the computational domain or domain entrance at Lx = 0. The value of up is varied in such a way that the axial component of total flow at the domain entrance is maintained identical to the prescribed axial component of incoming reactants U0.
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1111
Lx
u v xanchor
(flame anchor) y Ly
0
-L y
premixed reactants Fig. 6. Computational domain for turbulent premixed open V-flame.
4.3. Numerical experiment Following ChengÕs experimental set-up, the reactants under consideration in this study is an ethene (C2H4)/ air mixture with an equivalence ratio of 0.7. The density ratio of reactants to products is s = 6.7. The velocity of the flow entering the domain U0 is set as 5.5 ms1. Using U0 and the diameter of the inner coaxial cylinder used in ChengÕs experiment as the characteristic length (which is 50 mm), a Reynolds number, Re of 28,000 is yielded. The laminar flame velocity, SL is given as 0.44 ms1 and e is chosen as 1 mm in accordance with the range of values computed by Garcia-Ybarra et al. [26] for thermodiffusively stable hydrocarbon flames. The values of three important quantities relating to CAS must also be defined. They are the prescribed threshold scale, d, the characteristic length, L and the inverse magnification factor, l. The first quantity is required for deciding if contour surgery should be performed while the last two are primarily for computing the desired density of marker node distribution. Since d is closely linked to the connectivity of a flame front, it is therefore natural to select the value of laminar flame front thickness dl = 0.2 mm as d. It is clear from the description of L in the preceding section that it represents the typical distance between consecutive marker nodes which are able to characterize the flame when it is at its initial state. Similar to Lam et al., the nondimensional values of L is taken as 0.25 in the present study. The optimal value of l has to be determined by trial and error. This can be accomplished by first running the simulation with l set as say unity. This is then followed by repeating the run but with a smaller l value. The shape of the flame fronts resulting from the two simulations are then compared. In particular, the average distance between the two fronts are compared at various downstream locations. If sufficient convergence is acquired, then the lower l value is taken as optimal. In the case when convergence is not shown, the whole procedure is repeated but with progressively smaller values of l. In the present study, the optimal value of l used is 0.1. Upstream turbulence is incorporated by injecting 24 vortices at random along the domain entrance of y = ±60 mm. Care is taken to ensure that the same number of positive and negative vortices are introduced at each injection. The vortex radius is fixed at 1 mm throughout this study. Ten simulations are performed with the non-dimensional vortex circulation starting at 0.002 and ascending in step of 0.002 to a vortex circulation of 0.02. In all simulations, the non-dimensional time step dt is 0.005 based on the CFL condition. Vortices are introduced for every 10 time steps. A total of 5000 time steps are performed. All statistic results are obtained by averaging the instantaneous values after the 2000th time step.
1112
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
5. Results and discussion 5.1. Turbulence intensity Resulted upstream mean turbulence intensities are 1.8%, 3.5%, 4.4%, 7.0%, 8.3%, 11.0%, 12.5%, 15.2%, 16.9% and 19.8%. Fig. 7 shows profiles of rms axial velocity fluctuation u 0 upstream (x = 25 mm) and downstream (x = 100 mm) of the flame holder which is located at x = 50 mm. Profiles for three selected cases with turbulence levels of 7.0%, 12.5% and 19.8% (with vortex circulation of 0.008, 0.014 and 0.02, respectively) are presented. On the left hand side, profiles of upstream turbulence are shown. Clearly, upstream turbulence level increases with the circulation of injected vortices. Although the circulation was raised linearly, the increase of upstream turbulence does not follow in the same fashion. On the right hand side, profiles of downstream turbulence are presented. All three profiles begin to display an obvious increase in turbulence intensity at y = 10 mm. The profiles peak at a position just after y = 20 mm. When the flow passes through the flame, density decreases rapidly while temperature increases across the flame resulting in local amplifications of velocity fluctuations. Hence the peak on each profile represents the position of the flame front. From the figure, it can be deduced that the effect of increasing vortex circulation (hence turbulence intensity) on flame front position is negligible. It is also true that the peak (flame front position) is not as sharp at higher levels of turbulence intensities. This is because the extent of flame wrinkling is greater. Consequently, the flame has a thicker appearance and the corresponding peak becomes wider and less sharp. Profiles of rms transverse velocity fluctuation v 0 appear similar to that of u 0 and are shown in Fig. 8. 5.2. Flame surface density The flame surface density can be defined as the ratio of mean surface area to unit volume. For two-dimensional cases, such as the V-flame considered in the present study, R can be redefined as the ratio of the average flame length DL to the flame zone area DA RðhciÞ ¼
DLðhciÞ . DAðhciÞ
ð30Þ
All three terms in Eq. (30) are expressed as a function of the mean reaction progress variable hci. The mean reaction progress variable can be regarded as a parameter that indicates the extent of burning. It takes a value of zero in the unburnt region (reactants) and a value of unity in the fully burnt region (products). DL(hci) and DA(hci) correspond to the average flame length and the flame zone area for regions having mean reaction progress variable lying between hci ± Dhci. In this study, Dhci is chosen to be 1/800. 0.25
0.2
19.8% at 25 mm 12.5% at 25 mm 7.0% at 25 mm 19.8% at 100 mm 12.5% at 100 mm 7.0% at 100 mm
u' / U
0.15
0.1
0.05
0 -40
-30
-20
-10
0
10
20
30
40
y , mm Fig. 7. Profiles of rms axial velocity fluctuation u 0 at x = 25 mm (filled symbols) and at x = 100 mm (empty symbols).
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1113
0.25
0.2
19.8% at 25 mm 12.5% at 25 mm 7.0% at 25 mm 19.8% at 100 mm 12.5% at 100 mm 7.0% at 100 mm
v' / U
0.15
0.1
0.05
-40
-30
-20
-10
0
0
10
20
30
40
y, mm Fig. 8. Profiles of rms transverse velocity fluctuation v 0 at x = 25 mm (filled symbols) and at x = 100 mm (empty symbols).
Before proceeding, it is appropriate to elaborate further the concept of hci. To determine the reaction progress variable hci, the computational domain is first divided into a mesh with a suitable grid size. In the present study, a grid size of 0.1 mm · 0.1 mm is used. Under the flamelet regime, the reaction is assumed to take place at the flame front which is infinitely thin. Taking the shape of the flame front at a particular time step and superimposing it onto the mesh results in the mesh being segregated into burnt and unburnt regions. For the burnt region (i.e. area enclosed by the flame), each of the gird is assigned a value of unity. The remaining grids, which fall into the unburnt region, a value of zero is assigned. In situation when the grid actually contains a flame, it is treated as unburnt and thus a value of zero is allocated. This is repeated for all time steps. The assigned value on each of the grid is accumulated for all time steps and the time average gives the value of hci of that particular grid. A map of hci is thus obtained. Fig. 9 displays half the hci map for the case with simulated upstream turbulence level of 7%. Evenly spaced hci contours can be seen from the figure. Since hci equals unity in fully burnt products and zero in unburnt reactants, any value of hci between zero and unity
Fig. 9. The hci map obtained for the case with simulated upstream turbulence level of 7.0% (only the left side is shown).
1114
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
indicates the position of the mean flame. The streamwise increase in the thickness of the mean flame can be realized by inspecting the mean reaction progress variable hci. Using results from the same three simulations as before, Fig. 10 is constructed to display the variations of DL and DA with hci. Fig. 10 confirms that both DL and DA increase with turbulence. As mentioned in the preceding sub-section, the extent of flame wrinkling is greater at higher level of turbulence. More wrinkling means longer flame length, thus explaining the trend observed in Fig. 10. From individual subplots, the increase of DL with hci becomes more apparent as turbulence intensity gets higher. Since there is greater cusping of the flame front on the burnt side, DL is longer at higher hci. Also, increased flame wrinkling causes the flame to thicken, which is displayed in Fig. 10. All three DA profiles appear roughly symmetrical with a min1.4
30
1.2
25
Length Area
1 2
10
0.4
5
0.2 0
0 0.2
(a)
0.4 0.6 Reaction progress variable,
0.8
1.0
1.4
30
1.2
25
Length, mm
1
15 0.6 10
0.4
5
0.2
0 0.2
(b)
0.4
0.6
0.8
1.0
Reaction progress variable,
1.4
30
1.2
25
1 Length, mm
Area
20
0.8
0 0.0
Length
Area, mm2
0.0
15 0.6
Area
10 5
0.2 0 0.0
Length
20
0.8
0.4
(c)
Area, mm
15
0.6
Area, mm2
Length, mm
20
0.8
0 0.2
0.4
0.6
0.8
1.0
Reaction progress variable,
Fig. 10. Variations of mean flame length DL and mean flame area DA with reaction progress variation hci for three different levels of upstream turbulence (a) 7.0%, (b) 12.5% and (c) 19.8%.
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1115
imum at the centre of the flame zone hci = 0.5. As anticipated, DA increases sharply as hci tends to 0 (unburnt) and 1 (fully burnt). Using values of DL and DA, the flame surface density R is estimated and results are illustrated in Fig. 11. Asymmetry is observed for all profiles and are comparable in shape to those obtained by other researchers [16,27] for a two-dimensional V-shaped flame, under different turbulence conditions. Skewing towards the burnt side is also observed. Results from some laboratory experiments [28] indicate that the peak of the profile moves further away from hci = 0.5 as turbulence increases. However, this is not apparent in the present simulations. 5.3. Overall burning rate Following Shepherd [16], the overall burning rate is estimated as Z gb hW i ¼ R dg;
ð31Þ
gu
where g is an integration path which is normal to the hci contours and takes a value of zero at hci = 0.5. Fig. 11 illustrates that R is significantly less at higher turbulence intensity. One might expect that lower R constitute lower W. However, experimental data of Shepherd [16] illustrate that, W is in fact greater at higher turbulence intensity. Hence a higher R does not necessary imply a higher W. As shown in Fig. 10, the flame area increases with flame length and since W is estimated by integrating R through the flame zone, it is possible for cases having low R but a higher W. The overall burning rates for all 10 simulations are displayed in Fig. 12. 7.00%
0.25
Flame surface density , Σ
12.50% 0.2
19.80%
0.15
0.1
0.05
0 0
0.2
0.4
0.6
0.8
1
Reaction progress variable Fig. 11. Variation of flame surface density R with reaction progress variable hci.
Overall burning rate, W
4 3.5 3 2.5 2 1.5 1 0.5 0
0
5
10
15
Upstream turbulence % Fig. 12. Variation of overall burning rate, W with upstream turbulence.
20
1116
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
The minimum and maximum burning rates obtained are 1.26 and 3.36, respectively. It can be observed that the more intense the vortex circulation (and hence a higher upstream turbulence intensity), the greater the overall burning rate. A second-order polynomial can be fitted through the data points indicating that the rate of increase of W gets progressively higher with upstream turbulence. Using this figure to estimate W for the Vshaped flame studied by Shepherd [16] gives a value of 1.47, which is very close to the values shown in Table 2 of his paper. 6. Conclusions Simulations were carried out to test a numerical algorithm termed contour advection with surgery (CAS) when applied to turbulent premixed combustion. In order to acquire a better quantitative understanding of the numerical scheme, it was tested by varying the vortex circulation and this resulted in a wide range of upstream turbulence (1.8–19.8%). It was demonstrated that even at intense upstream turbulence (<15%), CAS can capture the flame front evolution with ease. Results confirm that both average flame length and flame zone area increase with turbulent intensity. Flame length is also longer in regions having higher values of mean progress variable hci. The flame surface density, R displayed a similar trend as those found in similar laboratory experiments: asymmetrical and peaked towards the burnt region. Acknowledgements This work was partially supported by a studentship of The Hong Kong Polytechnic University and a grant from the Research Committed of The Hong Kong Polytechnic University (Grant No. A-PD99). References [1] A. Linan, F.A. Williams, Fundamental Aspects of Combustion, Oxford University Press, Oxford, 1993. [2] C.K. Chan, K.S. Lau, B.L. Zhang, Simulation of a premixed turbulent flame with the discrete vortex method, Int. J. Numer. Methods Eng. 48 (2000) 613–627. [3] C.K. Chan, H.Y. Wang, H.Y. Tang, Effect of intense turbulence on turbulent premixed V-flame, Int. J. Eng. Sci. 41 (2003) 903–916. [4] C.W. Rhee, L. Talbot, J.A. Sethian, Dynamical behaviour of a premixed turbulent open V-flame, J. Fluid Mech. 300 (1995) 87–115. [5] M.Z. Pindera, L. Talbot, Some fluid dynamic considerations in the modeling of flames, Combust. Flame. 73 (1988) 111–125. [6] J. Qian, G. Tryggvason, C.K. Law, A front tracking method for the motion of premixed flames, J. Comput. Phys. 144 (1998) 52–69. [7] J.S.L. Lam, C.K. Chan, L. Talbot, I.G. Shepherd, On the high-resolution modelling of a turbulent premixed open V-flame, Combust. Theory Modell. 7 (2003) 1–28. [8] H.L. Berk, K.V. Roberts, The water-bag model, Meth. Comput. Phys. 9 (1967) 87–134. [9] N.J. Zabusky, M.H. Hughes, K.V. Roberts, Contour dynamics for the Euler equations in two dimensions, J. Comput. Phys. 30 (1979) 96–106. [10] D.G. Dritschel, CS-a topological reconnection scheme for extended integrations using contour dynamics, J. Comput. Phys. 77 (1988) 240–266. [11] D.G. Dritschel, Contour dynamics and CS: numerical algorithm for extended, high-resolution modelling of vortex dynamics in twodimensional, inviscid, incompressible flows, Comput. Phys. Rep. 10 (1989) 78–146. [12] D.W. Waugh, R.A. Plumb, Contour advection with surgery—a technique for investigating finescale structure in tracer transport, J. Atmos. Sci. 51 (1994) 530–540. [13] W.A. Norton, Breaking Rossby waves in a model stratosphere diagnosed by a vortex-following coordinate system and a technique for advecting material contours, J. Atmos. Sci. 51 (1994) 654–673. [14] A. Mariotti, M. Moustaoui, B. Legras, H. Teitelbaum, Comparison between vertical ozone soundings and reconstructed potential vorticity maps by contour advection with surgery, J. Geophys. Res. 102 (1997) 6131–6142. [15] A. Mariotti, C.R. Mechoso, B. Legras, V. Daniel, The evolution of the ozone ‘‘collar’’ in the Antarctic lower stratosphere during early August 1994, J. Atmos. Sci. 57 (2000) 402–414. [16] I.G. Shepherd, Flame surface density and burning rate in premixed turbulent flames, Proc. Combust. Inst. 26 (1996) 373–379. [17] G.H. Markstein, Experimental and theoretical studies of flame front stability, J. Aeronaut. Sci. 18 (1951) 199–209. [18] W.T. Ashurst, Vortex Simulation of Unsteady Wrinkled Laminar Flames, Combust. Sci. Technol. 52 (1987) 325–351. [19] H.D. Hayes, The vorticity jump across a gas dynamic discontinuity, J. Fluid Mech. 2 (1959) 595–600. [20] A.J. Chorin, Numerical study of slightly viscous flow, J. Fluid Mech. 57 (1973) 786–796. [21] J.A. Sethain, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Geometry, Fluid Mechanics, Computer Vision and Materials Science, Cambridge University Press, Cambridge, UK, 1999. [22] R.K. Cheng, Conditioned sampling of turbulence intensities and Reynolds stress in premixed turbulent flame, Combust. Sci. Technol. 41 (1984) 109–142.
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1117
[23] R.K. Cheng, T.T. Ng, Velocity statistics in premixed turbulent flames, Combust. Flame 52 (1983) 185–202. [24] R.K. Cheng, T.T. Ng, On defining the turbulent burning velocity in premixed V-shaped turbulent flames, Combust. Flame 57 (1984) 155–167. [25] R.K. Cheng, L. Talbot, F. Robben, Conditional velocity statistics in mixed CH4–AIR and C2H4–AIR turbulent flames, Proc. Combust. Inst. 20 (1984) 453–461. [26] P. Garcia-Ybarra, C. Nicoli, P. Clavin, Soret and dilution effects on premixed flames, Combust. Sci. Technol. 42 (1984) 87–109. [27] D. Veynante, J. Piana, J.M. Duclos, C. Martel, Experimental analysis of flame surface density model for premixed turbulent combustion, Proc. Combust. Inst. 26 (1996) 413–420. [28] G.G. Lee, K.Y. Huh, H. Kobayashi, Measurement and analysis of flame surface density for turbulent premixed combustion on a nozzle-type burner, Combust. Flame 122 (2000) 43–57.
Simulation of a turbulent premixed open V-shaped flame using contour advection with surgery B.H.Y. Tang
a,*
, C.K. Chan a, J.S.L. Lam
b
a
b
Department of Applied Mathematics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong Environmental Protection Department, The Government of Hong Kong Administrative Region, 24/F Southern Centre, 130 Hennessy Road, Wanchai, Hong Kong Received 17 January 2005; received in revised form 24 June 2005; accepted 24 August 2005 Available online 25 October 2005
Abstract Despite its capability of high spatial resolution, simulation of turbulent flows with traditional Lagrangian (front tracking) scheme is often discouraged by numerical instability caused by clustering of marker nodes and topological changes of fronts. Contour advection surgery (CAS), being a robust front tracking scheme, can limit the growth of front complexity during simulation without jeopardizing accuracy or efficiency. This allows it to open up an advantage over traditional front-tracking schemes. It has already been demonstrated that CAS, with incorporation of the reaction sheet model, can accurately simulate the propagation and advection of a turbulent premixed V-shaped flame. In this study, it is further tested with 10 values of vortex circulation. A range of upstream turbulence levels of 1.8–19.8% was obtained. Results indicate that upstream turbulence increase the average flame length, flame zone area and the overall burning rate. Flame surface density R was also estimated. Maximum values of R obtained lie in the range 0.1–1.4 mm1. Skewness towards the burnt region was observed in all profiles of R. Similar to results from laboratory experiments, it was found that R values decreases with upstream turbulence. From this study, the ability of CAS to cope with intense turbulence is demonstrated and a better quantitative understanding on the scheme has also been acquired. 2005 Elsevier Inc. All rights reserved. Keywords: Contour advection with surgery; V-shaped flame; Flame surface density; Burning rate
1. Introduction Turbulent premixed combustion is widely used in a range of engineering devices and has attracted many research interests over the years. In the field of turbulent premixed combustion, of primary research interest is the accurate simulation of a premixed flame propagating in an ambient turbulent flow. Due to the non-linear coupling of mechanical turbulence and combustion process, devising an accurate numerical method which is capable of capturing details of the propagation without masking the nature of physics, still remains a challenge today. *
Corresponding author. E-mail address: [email protected] (B.H.Y. Tang).
0307-904X/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2005.08.016
1104
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
In this paper, a two-dimensional rod-anchored turbulent V-shaped premixed flame is considered. This is one of the most common configurations employed for studying turbulent flame in laboratory experiments. Amongst the many theoretical models developed for simulating such flame numerically, the reaction sheet model [1] appears to be the most frequently invoked. This model allows the internal structure of the flame to be neglected by using the fact that, when the reaction rate is high, the flame front can be approximated as an infinitesimally thin boundary separating burnt and unburnt regions. With the exclusion of the flame frontÕs internal structure, the geometry of the flame front becomes the sole element that governs the evolution of the flow field. In turbulent premixed flames, two common occurrences namely, the topological changes (merger and breaking of the flame front) and the development of cusps, are the major influences on the geometry of the flame front. There are in general two approaches for capturing the flame front movement accurately. The first one is called the front-capturing approach which was proved to be very successful in numerous works [2–4]. However, being Eulerian in nature, the spatial resolution is somehow limited by the grid size of the computational domain. Thus, any sub-grid features would not be resolved by this grid-based approach. The second one, termed the front-tracking approach, being Lagrangian in nature, does not have such spatial resolution limitation. This approach has received some attention [5,6]. However, in general applications, it remains to be less popular than the Eulerian approach. This is because numerical instability is often caused by the two common occurrences in turbulent premixed flames mentioned earlier. Clearly, in order to exploit its advantage of high spatial resolution and thus capturing the sub-grid features, this numerical instability problem must first be solved. In view of this, Lam et al. [7] employed a front-tracking scheme known as contour advection with surgery to treat the topological changes of flame fronts and development of cusps. This front-tracking method originates from a numerical technique called contour dynamics (hereinafter CD) meaning the dynamics of boundaries of vorticity discontinuity (i.e. contours). Contour dynamics is commonly used in simulating inviscid fluid motions, particularly in vortex motion. It is based on the ‘‘water-bag’’ model for solving the Vlasov equation that appears in plasma physics [8] and was first introduced by Zabusky [9] in his simulation of two-dimensional fluid motions governed by the barotropic vorticity equation in a domain of infinite extent. One major drawback of CD that prevents its further application is the inability to deal with the rapid increase in the number of nodes used to represent the contours. To maintain a suitable spatial resolution, it is often necessary to continually insert extra nodes at regions of fast-growing curvature. This inevitably slows down the computation or in the worse case scenario, brings the calculation to a virtual standstill. Another drawback is the unmanageable growth of contour perimeters. Due to continual evolution of the contours, the contour perimeters will grow in time and very often resulting in a flow too complex for CD to resolve. In view of this, Dritschel [10,11] extended CD to ‘‘Contour Surgery’’ (hereinafter CS). It was introduced to especially treat the development of fine-scale features and thus limiting the growth in the number of nodes in a consistent manner. This is accomplished by truncating features that are formed under a predefined scale, d. It also topologically reconnects contours/divide a contour if the contours/parts of a contour become closer than d. The front-tracking scheme CAS was developed [12] after the success of CS and it was mainly used in geophysical research both as a numerical tool to advect fronts and a diagnostic tool to investigate fine-scale features [13–15]. Recently Lam et al. [7] have adopted the concept of CAS and applied it to a two-dimensional turbulent premixed open V-flame. Their computed velocity statistics and estimated flame surface density are in good agreement with corresponding laboratory measurements [16]. The present work, as well as being part of a series of extensive tests on CAS, is a continuation of the work done by Lam et al. [7] with the primary objective of acquiring a better quantitative understanding of this method. By varying the circulation of the vortices injected at the computational domain entrance, a wide range of upstream turbulence levels was obtained. Analysis was made on the resulting average flame length, flame area, flame surface density R and burning rate. 2. Basic equations The main assumptions of the model are: (i) the reaction is high so that the flame front is infinitesimally thin,
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1105
(ii) burnt and unburnt regions have distinct uniform densities, (iii) Mach number is very small so that the flow can be regarded as incompressible on either side of the flame, (iv) the mechanism of vorticity production is inviscid. At low Mach number, the velocity field, u for the combustion process can be decomposed into three components namely us the solenoidal component due to volume expansion across the flame front, ur the rotational component due to vorticity distribution x(x) and up the potential velocity field u ¼ us þ u r þ u p .
ð1Þ
The following three equations describe the conditions which have to be satisfied by individual components: r us ¼ mdðx xflm Þ;
r us ¼ 0;
ð2Þ
where m is the volume source strength along the flame front, x is the position vector, with subscript flm denoting the flame position and d( ) is the two-dimensional Dirac delta function r ur ¼ 0; up ¼ r/;
r ur ¼ xðxÞ;
ð3Þ ð4Þ
where / is the velocity potential of the incident flow. An expression for m is now derived. Let Su and Sb be the relative normal fluid velocities with respect to the flame front on the unburnt and burnt sides, respectively. The volume source strength is the net volume flux generated in the burnt products and consumed in the reactants, then m ¼ Sb Su.
ð5Þ
Denoting qu and qb as the fluid densities of the unburnt and burnt regions, respectively, the mass conservation law yields qu S u ¼ qb S b .
ð6Þ
From Eqs. (5) and (6), m can be expressed as m ¼ ðqu =qb 1ÞS u .
ð7Þ
According to Markstein [17], Su varies with curvature in such a way that a perturbed flame is stabilized as long as the curvature along the flame remains weak. Hence for weak curvature, Su can be approximated as S u ¼ S L ð1 þ ejÞ;
ð8Þ
where e is the Markstein length scale and j is the local curvature which is taken as positive/negative when the centre of the circular arc lies left/right of the flame front. Since Eq. (8) does not involve the use of any correction term, the effect of strain on the flame propagation velocity is not modelled. According to Ashurst [18], such exclusion of strain term is possible. He modelled a turbulent V-flame under conditions similar to some of that employed in the present study with the full stretch term and found that the velocity statistics did not show any dependence on the strain term. In the present simulation, the vorticity x(x) has two sources, namely the prescribed upstream turbulence and the flame front itself. The vorticity generated from the latter source is also known as flame-induced vorticity. Pressure gradients tangential to the flame result in different accelerations in the unburnt and burnt sides with distinct densities. Vorticity is thus generated at the flame by the mean density gradient across the flame and the pressure gradient tangential to the flame. The upstream source is simulated by injecting small uniform circular vortices at the domain entrance. To account for the flame-induced vorticity, the flame front is divided into small segments of equal length Ds with one small uniform circular vortex being introduced immediatelypbehind the mid-point of each segment in each time step Dt. The radius of this flame-induced vortex is given by (SuDsDt/p) so as to conserve the volume flux across the flame front in each time step. The determination of the flame-induced vorticity involves the use of an expression proposed by Hayes [19] for the vorticity jump across the flame front 1 1 q qu ½x ¼ fDut þ ut rt ut ut un j un rt un g; ð9Þ rt ðqu S u Þ þ b qb qu qu S u
1106
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
where ut and un are respectively the relative tangential and absolute normal components with respect to the flame front, $t is the gradient along the flame front and D is the material time derivative taken at a point which always lies on the flame front and moves in a direction normal to the discontinuity as it moves. Both the turbulence vortices and the flame-induced vortices are treated by the Random Vortex Method [20]. It can simulate the diffusive effect of viscosity on vorticity distribution and prevent the clumping of vortices about the flame front. 3. Numerical algorithm—CAS Since CAS was not originally devised for turbulent reacting flow, it is necessary to make modifications from the original. Following Lam et al. [7], the essential details of steps involved in CAS are presented in the following subsections. 3.1. Geometric configuration of front As with typical Lagrangian schemes, the front is discretized into an orderly set of connected marker nodes {xi} as illustrated in Fig. 1. Arrows on the front signify its orientation. In the present study, the flame front is oriented such that the unburnt region (reactants) always lies on the left while the burnt region (fully burnt products) always lies on the right. The displacement between two consecutive marker nodes xi and xi+1 is denoted by ei, i.e. ei ðai ; bi Þ ¼ xiþ1 xi .
ð10Þ
For ei, an outward normal vector fi is defined as fi ¼ ðbi ; ai Þ;
ð11Þ
which points towards the reactant as illustrated in Fig. 1. At any marker node xi, the continuity of tangent slope, ti is ensured by using the approximation ei1 ei ti ¼ þ 2. 2 jei1 j jei j Similar to ei, a corresponding outward normal ni can be defined for ti ^ ti ; ni ¼ k jti j
Flame front
right
left UNBURNT ti
xi+1 ei
fi
BURNT xi ni ei-1 xi-1
Fig. 1. Geometric configuration of a discretized front {xi}.
ð12Þ
ð13Þ
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1107
ti ni ti xi
ni
xi (a) κ i > 0
(b) κ i < 0
Fig. 2. Sign convention of curvature ji.
^ is the upward unit vector normal to the plane of the flow domain. In general, the two orthogonal pairs where k {ei, fi} and {ti, ni} are not identical. At any marker node xi, the curvature is determined by fitting an arc through its preceding and following nodes i.e. nodes xi1, xi and xi+1. It is then given as ji ¼
2ðai1 bi bi1 ai Þ 2
jti jjei j jei1 j
.
ð14Þ
Evidently, the curvature defined in Eq. (14) is signed. As illustrated in Fig. 2, ji takes a positive/negative value if the centre of the arc lies to the left/right of the front. 3.2. Contour surgery In every time step, each of the marker node is advected using the fourth-order Runge–Kutta scheme. Once all the marker nodes are advected and propagated to a new position, contour surgery is performed if necessary. At the end of each time step, the distance between any two non-adjacent marker nodes is computed. If this distance is less than the predefined threshold value, d (note that d no longer represents the delta Dirac function as in the preceding chapter), a so-called surgery will be carried out. There are three kinds of CS, namely fission, fusion and shorten surgery. The decision of which type of surgery to perform is based on the positions of the two non-adjacent marker nodes. If they are along the same front, a fission surgery is performed and the front is separated into two disconnected fronts as illustrated in Fig. 3(a). If the two non-adjacent marker nodes are on two disconnected front, a fusion surgery is performed as depicted in Fig. 3(b). In the simulation of flame front propagation, overshooting can cause an undesirable effect of unburnt reactants appearing as burnt product. Both types of aforementioned CS are effective in preventing overshooting of fronts and the associated numerical errors. In addition, even when overshooting does not occur, it is possible for fronts/parts of a front to coincide if the distance between them is small enough. This results in a redundant and dynamically ineffective front, as illustrated in Fig. 4. With surgery, this is avoided because any two fronts/ parts of a front which have a distance less than d between them is treated immediately. There exists a special type of fission surgery namely, shorten surgery. When the fission surgery is invoked, two disconnected fronts are form. If one of these disconnected fronts has a spatial extent that is less than the threshold scale d, it is eradicated as depicted in Fig. 3(c) and this is the shorten surgery. With this surgery, the problem of numerical instability triggered by the clustering of marker nodes at regions of high curvature can be resolved. Both fission and shorten surgeries are in fact the so-called Ôde-loopingÕ techniques [21] which remove certain marker nodes on the front so as to prevent the formation any needless loops (filamentary cusps). 3.3. Redistribution of marker nodes To maintain a suitable resolution and prevent clustering of marker nodes at a region of high curvature, it is essential to redistribute the marker nodes after surgery. Before the redistribution, each front is first split into
1108
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
Fig. 3. Three types of contour surgery (a) fission surgery, (b) fusion surgery and (c) shorten surgery.
Fig. 4. Schematic evolution of fronts following Fig. 3(a) and (b) without (left) and with (right) contour surgery.
curved segments delimited by extremities and/or corners. A marker node xi is labelled as a corner if the displacement vector ei and ei1 form an acute angle. For each curved segment, an appropriate average node density, ki, between each pair of adjacent marker nodes is calculated as follows: 1 1 ~ 1 ~ ~i þ ; j ki ¼ min ; ð15Þ d 2l L ~ ~ ~iþ1 ~i þ j j ~ ~ ~i ¼ ; ð16Þ j 2 ð~ ji1 =jei1 jÞ þ ð~ ji =jei jÞ ~ ~i ¼ j ; ð17Þ ð1=jei1 jÞ þ ð1=jei jÞ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ji þ jiþ1 2 ~i ¼ ; ð18Þ j þ 2 L2 where L is the characteristic length scale which depends on the characteristic radius of curvature of the initial flame front and l 6 1 is a positive non-dimensional input parameter for overall density of marker nodes. The physical significance of these two parameters will be discussed later in this sub-section. Examining Eq. (15) more closely, one can deduce that ki is an increasing function of ji and must be less than or equal to 1/d. The former condition ensures adequate resolution by placing more marker nodes at region of high curvature while the latter condition prevents clustering of marker nodes by limiting the minimum distance
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1109
between any two adjacent marker nodes to nearly d. Eqs. (16)–(18) give immediate values of averaged curvatures (weighted or not) and are all positive (P1/L). As commented by Dritschel [10,11], to determine the optimal placement of nodes, one requires more than the local instantaneous curvature of the contour. This is because any approaching small-scale features will not be detected with the use of a local method of node adjustment. Hence Eqs. (16)–(18) are formulated such that they depend on not only the relative positions of three consecutive marker nodes but also on the geometry of a broader neighbourhood. To explain the physical significance of L and d as mentioned previously, one may begin with the minimum density of marker nodes. Clearly, minimum density of marker nodes is obtained with a straight front. From Eqs. (14) and (16)–(18), the local and non-local curvatures are respectively zero and 1/L. Note that the non-local curvature is non-zero even for a straight front. In this way, marker nodes are ÔreservedÕ and further addition is possible as the straight front evolves. Rewritting Eq. (15) for a straight front, the minimum density of marker nodes is 1 1 ki ¼ min ; . ð19Þ d lL Since 1/d is the upper limit for ki, one must select L and d such that d < lL.
ð20Þ
In fact, the relation described by Eq. (20) holds not only for a straight but also for any curved front. This is because for a curved front, ki is bounded as shown below 1 1 6 ki 6 . lL d
ð21Þ
Returning to the straight front case, setting the non-dimensional parameter l as unity i.e. its maximum value, would give the minimum density of marker nodes as 1 ki ¼ . ð22Þ L This means that at least one marker node must be inserted per arc-length of L. Decreasing the value of l simply means that the number of marker nodes per arc-length of L multiplies. It follows that l is in fact an inverse magnification factor that reserves extra marker nodes for the subsequent evolution of the turbulent flame front. After redistributing marker nodes with the results from Eqs. (15)–(18) for each curved segments, interpolation of marker nodes is carried out. Following Lam et al. [7], consecutive marker nodes xi and xi+1 are interpolated using piecewise cubic spline as follows: xðsÞ ¼ xi þ sei þ gðsÞf i ;
ð23Þ
where s is the fractional distance between marker nodes xi and xi+1 and gðsÞ ¼ ai s þ bi s2 þ ci s3
ð24Þ
whose coefficients are given by 1 ai ¼ ð2ji þ jiþ1 Þjei j; 6 1 bi ¼ ji jei j; 2 1 ci ¼ ðjiþ1 ji Þjei j. 6
ð25Þ ð26Þ ð27Þ
4. Applying CAS to a turbulent premixed V-shaped flame 4.1. Experimental set-up Similar to Chan et al. [2,3], the numerical conditions in the present study are based on a series of V-flame laboratory experiments performed by Cheng [22–25]. Fig. 5 depicts the experimental set-up employed in
1110
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117 Turbulent Flame 150 mm
Su
U
Flameholder
100
V 50
Turbulence
x
Generator
y 0
Air
C 2H 4
Air
+Air
Fig. 5. Experimental set-up for investigating a rod-stabilized V-shaped flame.
ChengÕs investigation. It consists of a coaxial nozzle providing a circular jet with an inner core of ethene and air. This premixed fuel is shielded from the ambient environment with an outer layer of air. A turbulence grid is placed 50 mm upstream of the nozzle exit. When the premixed fuel passes the grid, turbulence is generated within the flow. A circular rod with a diameter of 1 mm is placed at the nozzle exit. When the premixed fuel is ignited, a turbulent V-shaped flame anchoring at the circular rod would be observed. This turbulent V-shaped flame travels with a normal velocity Su with respect to the reactants. While the V-flame problem is threedimensional in nature, the third velocity component in the z-direction is ignored in ChengÕs experiments. Thus it is reasonable to assume the flow is essentially two-dimensional. 4.2. Computational domain With CAS, the computational domain is unbounded. However, for a proper treatment of the flame front, one must define a spatial extent because an infinite flame cannot be achieved physically or computationally. Based on the size of the V-flame obtained in ChengÕs experiment, the solution domain is truncated to a rectangle of axial length (x-direction) equal to 150 mm and transverse width (y-direction) equal to 120 mm. For convenience, the present problem is transformed into a non-dimensional one. The non-dimensionalising length scale is selected to match the diameter of the inner coaxial cylinder of the nozzle. While the non-dimensionalising time scale is chosen such that the non-dimensional inflow velocity, U0 can be scaled to unity. Thus the following length and time scales are chosen for scaling all parameters: L ¼ 50 mm; 1 s. T ¼ 110
ð28Þ ð29Þ
Consequently, the computational domain now becomes non-dimensional with axial length, Lx equal to 3 and transverse width, Ly equal to 1.2. This is illustrated in Fig. 6. There is one boundary condition imposed at the lower boundary of the computational domain or domain entrance at Lx = 0. The value of up is varied in such a way that the axial component of total flow at the domain entrance is maintained identical to the prescribed axial component of incoming reactants U0.
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1111
Lx
u v xanchor
(flame anchor) y Ly
0
-L y
premixed reactants Fig. 6. Computational domain for turbulent premixed open V-flame.
4.3. Numerical experiment Following ChengÕs experimental set-up, the reactants under consideration in this study is an ethene (C2H4)/ air mixture with an equivalence ratio of 0.7. The density ratio of reactants to products is s = 6.7. The velocity of the flow entering the domain U0 is set as 5.5 ms1. Using U0 and the diameter of the inner coaxial cylinder used in ChengÕs experiment as the characteristic length (which is 50 mm), a Reynolds number, Re of 28,000 is yielded. The laminar flame velocity, SL is given as 0.44 ms1 and e is chosen as 1 mm in accordance with the range of values computed by Garcia-Ybarra et al. [26] for thermodiffusively stable hydrocarbon flames. The values of three important quantities relating to CAS must also be defined. They are the prescribed threshold scale, d, the characteristic length, L and the inverse magnification factor, l. The first quantity is required for deciding if contour surgery should be performed while the last two are primarily for computing the desired density of marker node distribution. Since d is closely linked to the connectivity of a flame front, it is therefore natural to select the value of laminar flame front thickness dl = 0.2 mm as d. It is clear from the description of L in the preceding section that it represents the typical distance between consecutive marker nodes which are able to characterize the flame when it is at its initial state. Similar to Lam et al., the nondimensional values of L is taken as 0.25 in the present study. The optimal value of l has to be determined by trial and error. This can be accomplished by first running the simulation with l set as say unity. This is then followed by repeating the run but with a smaller l value. The shape of the flame fronts resulting from the two simulations are then compared. In particular, the average distance between the two fronts are compared at various downstream locations. If sufficient convergence is acquired, then the lower l value is taken as optimal. In the case when convergence is not shown, the whole procedure is repeated but with progressively smaller values of l. In the present study, the optimal value of l used is 0.1. Upstream turbulence is incorporated by injecting 24 vortices at random along the domain entrance of y = ±60 mm. Care is taken to ensure that the same number of positive and negative vortices are introduced at each injection. The vortex radius is fixed at 1 mm throughout this study. Ten simulations are performed with the non-dimensional vortex circulation starting at 0.002 and ascending in step of 0.002 to a vortex circulation of 0.02. In all simulations, the non-dimensional time step dt is 0.005 based on the CFL condition. Vortices are introduced for every 10 time steps. A total of 5000 time steps are performed. All statistic results are obtained by averaging the instantaneous values after the 2000th time step.
1112
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
5. Results and discussion 5.1. Turbulence intensity Resulted upstream mean turbulence intensities are 1.8%, 3.5%, 4.4%, 7.0%, 8.3%, 11.0%, 12.5%, 15.2%, 16.9% and 19.8%. Fig. 7 shows profiles of rms axial velocity fluctuation u 0 upstream (x = 25 mm) and downstream (x = 100 mm) of the flame holder which is located at x = 50 mm. Profiles for three selected cases with turbulence levels of 7.0%, 12.5% and 19.8% (with vortex circulation of 0.008, 0.014 and 0.02, respectively) are presented. On the left hand side, profiles of upstream turbulence are shown. Clearly, upstream turbulence level increases with the circulation of injected vortices. Although the circulation was raised linearly, the increase of upstream turbulence does not follow in the same fashion. On the right hand side, profiles of downstream turbulence are presented. All three profiles begin to display an obvious increase in turbulence intensity at y = 10 mm. The profiles peak at a position just after y = 20 mm. When the flow passes through the flame, density decreases rapidly while temperature increases across the flame resulting in local amplifications of velocity fluctuations. Hence the peak on each profile represents the position of the flame front. From the figure, it can be deduced that the effect of increasing vortex circulation (hence turbulence intensity) on flame front position is negligible. It is also true that the peak (flame front position) is not as sharp at higher levels of turbulence intensities. This is because the extent of flame wrinkling is greater. Consequently, the flame has a thicker appearance and the corresponding peak becomes wider and less sharp. Profiles of rms transverse velocity fluctuation v 0 appear similar to that of u 0 and are shown in Fig. 8. 5.2. Flame surface density The flame surface density can be defined as the ratio of mean surface area to unit volume. For two-dimensional cases, such as the V-flame considered in the present study, R can be redefined as the ratio of the average flame length DL to the flame zone area DA RðhciÞ ¼
DLðhciÞ . DAðhciÞ
ð30Þ
All three terms in Eq. (30) are expressed as a function of the mean reaction progress variable hci. The mean reaction progress variable can be regarded as a parameter that indicates the extent of burning. It takes a value of zero in the unburnt region (reactants) and a value of unity in the fully burnt region (products). DL(hci) and DA(hci) correspond to the average flame length and the flame zone area for regions having mean reaction progress variable lying between hci ± Dhci. In this study, Dhci is chosen to be 1/800. 0.25
0.2
19.8% at 25 mm 12.5% at 25 mm 7.0% at 25 mm 19.8% at 100 mm 12.5% at 100 mm 7.0% at 100 mm
u' / U
0.15
0.1
0.05
0 -40
-30
-20
-10
0
10
20
30
40
y , mm Fig. 7. Profiles of rms axial velocity fluctuation u 0 at x = 25 mm (filled symbols) and at x = 100 mm (empty symbols).
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1113
0.25
0.2
19.8% at 25 mm 12.5% at 25 mm 7.0% at 25 mm 19.8% at 100 mm 12.5% at 100 mm 7.0% at 100 mm
v' / U
0.15
0.1
0.05
-40
-30
-20
-10
0
0
10
20
30
40
y, mm Fig. 8. Profiles of rms transverse velocity fluctuation v 0 at x = 25 mm (filled symbols) and at x = 100 mm (empty symbols).
Before proceeding, it is appropriate to elaborate further the concept of hci. To determine the reaction progress variable hci, the computational domain is first divided into a mesh with a suitable grid size. In the present study, a grid size of 0.1 mm · 0.1 mm is used. Under the flamelet regime, the reaction is assumed to take place at the flame front which is infinitely thin. Taking the shape of the flame front at a particular time step and superimposing it onto the mesh results in the mesh being segregated into burnt and unburnt regions. For the burnt region (i.e. area enclosed by the flame), each of the gird is assigned a value of unity. The remaining grids, which fall into the unburnt region, a value of zero is assigned. In situation when the grid actually contains a flame, it is treated as unburnt and thus a value of zero is allocated. This is repeated for all time steps. The assigned value on each of the grid is accumulated for all time steps and the time average gives the value of hci of that particular grid. A map of hci is thus obtained. Fig. 9 displays half the hci map for the case with simulated upstream turbulence level of 7%. Evenly spaced hci contours can be seen from the figure. Since hci equals unity in fully burnt products and zero in unburnt reactants, any value of hci between zero and unity
Fig. 9. The hci map obtained for the case with simulated upstream turbulence level of 7.0% (only the left side is shown).
1114
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
indicates the position of the mean flame. The streamwise increase in the thickness of the mean flame can be realized by inspecting the mean reaction progress variable hci. Using results from the same three simulations as before, Fig. 10 is constructed to display the variations of DL and DA with hci. Fig. 10 confirms that both DL and DA increase with turbulence. As mentioned in the preceding sub-section, the extent of flame wrinkling is greater at higher level of turbulence. More wrinkling means longer flame length, thus explaining the trend observed in Fig. 10. From individual subplots, the increase of DL with hci becomes more apparent as turbulence intensity gets higher. Since there is greater cusping of the flame front on the burnt side, DL is longer at higher hci. Also, increased flame wrinkling causes the flame to thicken, which is displayed in Fig. 10. All three DA profiles appear roughly symmetrical with a min1.4
30
1.2
25
Length Area
1 2
10
0.4
5
0.2 0
0 0.2
(a)
0.4 0.6 Reaction progress variable,
0.8
1.0
1.4
30
1.2
25
Length, mm
1
15 0.6 10
0.4
5
0.2
0 0.2
(b)
0.4
0.6
0.8
1.0
Reaction progress variable,
1.4
30
1.2
25
1 Length, mm
Area
20
0.8
0 0.0
Length
Area, mm2
0.0
15 0.6
Area
10 5
0.2 0 0.0
Length
20
0.8
0.4
(c)
Area, mm
15
0.6
Area, mm2
Length, mm
20
0.8
0 0.2
0.4
0.6
0.8
1.0
Reaction progress variable,
Fig. 10. Variations of mean flame length DL and mean flame area DA with reaction progress variation hci for three different levels of upstream turbulence (a) 7.0%, (b) 12.5% and (c) 19.8%.
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1115
imum at the centre of the flame zone hci = 0.5. As anticipated, DA increases sharply as hci tends to 0 (unburnt) and 1 (fully burnt). Using values of DL and DA, the flame surface density R is estimated and results are illustrated in Fig. 11. Asymmetry is observed for all profiles and are comparable in shape to those obtained by other researchers [16,27] for a two-dimensional V-shaped flame, under different turbulence conditions. Skewing towards the burnt side is also observed. Results from some laboratory experiments [28] indicate that the peak of the profile moves further away from hci = 0.5 as turbulence increases. However, this is not apparent in the present simulations. 5.3. Overall burning rate Following Shepherd [16], the overall burning rate is estimated as Z gb hW i ¼ R dg;
ð31Þ
gu
where g is an integration path which is normal to the hci contours and takes a value of zero at hci = 0.5. Fig. 11 illustrates that R is significantly less at higher turbulence intensity. One might expect that lower R constitute lower W. However, experimental data of Shepherd [16] illustrate that, W is in fact greater at higher turbulence intensity. Hence a higher R does not necessary imply a higher W. As shown in Fig. 10, the flame area increases with flame length and since W is estimated by integrating R through the flame zone, it is possible for cases having low R but a higher W. The overall burning rates for all 10 simulations are displayed in Fig. 12. 7.00%
0.25
Flame surface density , Σ
12.50% 0.2
19.80%
0.15
0.1
0.05
0 0
0.2
0.4
0.6
0.8
1
Reaction progress variable
Overall burning rate, W
4 3.5 3 2.5 2 1.5 1 0.5 0
0
5
10
15
Upstream turbulence % Fig. 12. Variation of overall burning rate, W with upstream turbulence.
20
1116
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
The minimum and maximum burning rates obtained are 1.26 and 3.36, respectively. It can be observed that the more intense the vortex circulation (and hence a higher upstream turbulence intensity), the greater the overall burning rate. A second-order polynomial can be fitted through the data points indicating that the rate of increase of W gets progressively higher with upstream turbulence. Using this figure to estimate W for the Vshaped flame studied by Shepherd [16] gives a value of 1.47, which is very close to the values shown in Table 2 of his paper. 6. Conclusions Simulations were carried out to test a numerical algorithm termed contour advection with surgery (CAS) when applied to turbulent premixed combustion. In order to acquire a better quantitative understanding of the numerical scheme, it was tested by varying the vortex circulation and this resulted in a wide range of upstream turbulence (1.8–19.8%). It was demonstrated that even at intense upstream turbulence (<15%), CAS can capture the flame front evolution with ease. Results confirm that both average flame length and flame zone area increase with turbulent intensity. Flame length is also longer in regions having higher values of mean progress variable hci. The flame surface density, R displayed a similar trend as those found in similar laboratory experiments: asymmetrical and peaked towards the burnt region. Acknowledgements This work was partially supported by a studentship of The Hong Kong Polytechnic University and a grant from the Research Committed of The Hong Kong Polytechnic University (Grant No. A-PD99). References [1] A. Linan, F.A. Williams, Fundamental Aspects of Combustion, Oxford University Press, Oxford, 1993. [2] C.K. Chan, K.S. Lau, B.L. Zhang, Simulation of a premixed turbulent flame with the discrete vortex method, Int. J. Numer. Methods Eng. 48 (2000) 613–627. [3] C.K. Chan, H.Y. Wang, H.Y. Tang, Effect of intense turbulence on turbulent premixed V-flame, Int. J. Eng. Sci. 41 (2003) 903–916. [4] C.W. Rhee, L. Talbot, J.A. Sethian, Dynamical behaviour of a premixed turbulent open V-flame, J. Fluid Mech. 300 (1995) 87–115. [5] M.Z. Pindera, L. Talbot, Some fluid dynamic considerations in the modeling of flames, Combust. Flame. 73 (1988) 111–125. [6] J. Qian, G. Tryggvason, C.K. Law, A front tracking method for the motion of premixed flames, J. Comput. Phys. 144 (1998) 52–69. [7] J.S.L. Lam, C.K. Chan, L. Talbot, I.G. Shepherd, On the high-resolution modelling of a turbulent premixed open V-flame, Combust. Theory Modell. 7 (2003) 1–28. [8] H.L. Berk, K.V. Roberts, The water-bag model, Meth. Comput. Phys. 9 (1967) 87–134. [9] N.J. Zabusky, M.H. Hughes, K.V. Roberts, Contour dynamics for the Euler equations in two dimensions, J. Comput. Phys. 30 (1979) 96–106. [10] D.G. Dritschel, CS-a topological reconnection scheme for extended integrations using contour dynamics, J. Comput. Phys. 77 (1988) 240–266. [11] D.G. Dritschel, Contour dynamics and CS: numerical algorithm for extended, high-resolution modelling of vortex dynamics in twodimensional, inviscid, incompressible flows, Comput. Phys. Rep. 10 (1989) 78–146. [12] D.W. Waugh, R.A. Plumb, Contour advection with surgery—a technique for investigating finescale structure in tracer transport, J. Atmos. Sci. 51 (1994) 530–540. [13] W.A. Norton, Breaking Rossby waves in a model stratosphere diagnosed by a vortex-following coordinate system and a technique for advecting material contours, J. Atmos. Sci. 51 (1994) 654–673. [14] A. Mariotti, M. Moustaoui, B. Legras, H. Teitelbaum, Comparison between vertical ozone soundings and reconstructed potential vorticity maps by contour advection with surgery, J. Geophys. Res. 102 (1997) 6131–6142. [15] A. Mariotti, C.R. Mechoso, B. Legras, V. Daniel, The evolution of the ozone ‘‘collar’’ in the Antarctic lower stratosphere during early August 1994, J. Atmos. Sci. 57 (2000) 402–414. [16] I.G. Shepherd, Flame surface density and burning rate in premixed turbulent flames, Proc. Combust. Inst. 26 (1996) 373–379. [17] G.H. Markstein, Experimental and theoretical studies of flame front stability, J. Aeronaut. Sci. 18 (1951) 199–209. [18] W.T. Ashurst, Vortex Simulation of Unsteady Wrinkled Laminar Flames, Combust. Sci. Technol. 52 (1987) 325–351. [19] H.D. Hayes, The vorticity jump across a gas dynamic discontinuity, J. Fluid Mech. 2 (1959) 595–600. [20] A.J. Chorin, Numerical study of slightly viscous flow, J. Fluid Mech. 57 (1973) 786–796. [21] J.A. Sethain, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Geometry, Fluid Mechanics, Computer Vision and Materials Science, Cambridge University Press, Cambridge, UK, 1999. [22] R.K. Cheng, Conditioned sampling of turbulence intensities and Reynolds stress in premixed turbulent flame, Combust. Sci. Technol. 41 (1984) 109–142.
B.H.Y. Tang et al. / Applied Mathematical Modelling 30 (2006) 1103–1117
1117
[23] R.K. Cheng, T.T. Ng, Velocity statistics in premixed turbulent flames, Combust. Flame 52 (1983) 185–202. [24] R.K. Cheng, T.T. Ng, On defining the turbulent burning velocity in premixed V-shaped turbulent flames, Combust. Flame 57 (1984) 155–167. [25] R.K. Cheng, L. Talbot, F. Robben, Conditional velocity statistics in mixed CH4–AIR and C2H4–AIR turbulent flames, Proc. Combust. Inst. 20 (1984) 453–461. [26] P. Garcia-Ybarra, C. Nicoli, P. Clavin, Soret and dilution effects on premixed flames, Combust. Sci. Technol. 42 (1984) 87–109. [27] D. Veynante, J. Piana, J.M. Duclos, C. Martel, Experimental analysis of flame surface density model for premixed turbulent combustion, Proc. Combust. Inst. 26 (1996) 413–420. [28] G.G. Lee, K.Y. Huh, H. Kobayashi, Measurement and analysis of flame surface density for turbulent premixed combustion on a nozzle-type burner, Combust. Flame 122 (2000) 43–57.