Three-dimensional flame displacement speed and flame front curvature measurements using quad-plane PIV

Three-dimensional flame displacement speed and flame front curvature measurements using quad-plane PIV

Combustion and Flame 160 (2013) 2757–2769 Contents lists available at SciVerse ScienceDirect Combustion and Flame j o u r n a l h o m e p a g e : w ...

4MB Sizes 0 Downloads 8 Views

Combustion and Flame 160 (2013) 2757–2769

Contents lists available at SciVerse ScienceDirect

Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e

Three-dimensional flame displacement speed and flame front curvature measurements using quad-plane PIV Johannes Kerl a, Chris Lawn b, Frank Beyrau a,⇑ a b

Department of Mechanical Engineering, Imperial College London, London SW7 2AZ, United Kingdom School of Engineering and Materials Science, Queen Mary, University of London, London E1 4NS, United Kingdom

a r t i c l e

i n f o

Article history: Received 6 October 2012 Received in revised form 29 June 2013 Accepted 2 July 2013 Available online 26 July 2013 Keywords: Laser diagnostics Particle image velocimetry Lean premixed flames Flame displacement speed Multi-plane Burning velocity

a b s t r a c t In this work we demonstrate the use of a quad-plane particle image velocimetry technique in order to investigate the three-dimensional behavior of several important quantities for combustion research such as the flame displacement speed and the flame front curvature. In results from a premixed methane flame stabilized in a diffuser burner, a comparison of three-dimensional and two-dimensional data is made in order to critically analyze the error of the usually performed planar measurements. It is shown that twodimensional measurements can only give an estimate of the real situation under certain circumstances, such as with mainly spherical structures in the flame and the perfect alignment of both the flame propagation and flow direction to the measurement plane. However, in turbulent flames, this alignment can never be achieved due to fluctuations stemming from turbulence. The application of two crossed planes leads to significant improvements and good agreement with the three-dimensional quantities can be observed, although no perfect match is achieved. Flame displacement speeds ranging from 0.4 sL to 4.5 sL with a mean of 1.1 sL were recorded, but were not correlated with the flame curvature or strain rate. Ó 2013 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Lean premixed turbulent combustion is very attractive for many technical applications because it can greatly reduce the emission of pollutants such as CO, NOx and soot particulates, while allowing high power densities. For the stability of lean premixed flames the interaction of reaction chemistry and turbulent flow field plays an important role, but the understanding of this coupling remains incomplete. Well-designed laser-based techniques applied to model combustors can help to provide insight into these fundamental processes, offering non-intrusive measurements of flow properties such as velocities and strain rates, as well as chemistry-related properties such as the flame speed. However, the applied laser techniques are typically only two-dimensional and hence can capture only an incomplete picture of the real, three dimensional nature of the flame. Recently, direct numerical simulations have shown that such planar measurements provide only limited access to the true three-dimensional properties sought after [1], unless assumptions are made which have only limited applicability in practice [2,3]. In this study, the challenge of measuring three-dimensional properties of lean premixed flames was approached by expanding the existing planar measurement technique of particle image ⇑ Corresponding author. E-mail address: [email protected] (F. Beyrau).

velocimetry (PIV) to four light sheets that cross in a common intersection line. The motion of seeded silicon oil droplets (Dow Corning 200/50cs) was measured to yield the velocity field in each of the investigated planes using a conventional PIV approach. Additionally, the evaporation of the seeded droplets at around 600 K can be exploited in premixed flames within the thin flame regime to identify the location of the flame front [3–8] and thus derive other quantities such as the local flame curvature and shape. When the time separation of two consecutive frames in PIV was chosen to be sufficiently high, the displacement of the flame front between the two exposures could also be determined, allowing simultaneous measurement of the flame displacement speed. The data from each of the planes was then combined to yield several three-dimensional parameters of the flame near the intersection line. In the following we explain how these properties were measured in the past and where the new technique uses established concepts or extends them. To investigate the strain rate in a three-dimensional way, the flow field needs to be determined in three dimensions. This can be achieved by applying multiple light sheet techniques such as dual-plane or crossed-plane PIV arrangements which yield a fully three-dimensional description of the flow [9–13]. In this work, a crossed-plane arrangement was favored over parallel plane approaches due to its several advantages described later on. The flame surface normal can be determined by a three dimensional technique such as described in [6]. For this a crossed-plane

0010-2180/$ - see front matter Ó 2013 The Combustion Institute. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.combustflame.2013.07.002

2758

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

arrangement of light sheets needs to be created in which the determination of the flame surface is possible in each of them. This can be achieved by means of Rayleigh thermometry [14], OH-LIF [2,11] or from PIV raw images, via the step-like change in particle number density [15] or via the disappearance of oil droplets at the flame front [6,8]. As the PIV images are readily available from determining the velocity field, this option was chosen to obtain information on the flame front normal too. But two crossed planes are still not capable of giving an accurate enough description of the flame surface. Therefore, the system was extended by another two crossed light sheets at 45° to the original ones. For these, also PIV was used because of the already available seeding particles, the easy separation of signals by color and the redundancy in flow field information. With this, it is possible to fit a surface through all four light sheets to obtain a reasonably accurate estimate of the flame surface structure and find the important properties of three-dimensional flame curvature and shape factor. The flame displacement speed is defined as the movement of the flame perpendicular to its surface within a certain time relative to the unburnt gas and can be determined from the apparent flame displacement and the local convective flow movement. Practically, it is determined by displacing the flame front in the first set of images by the advection acting on it and then calculating the distance this ‘virtual’ flame front has to travel to arrive at the position in the second set of images taken a certain time interval later. This concept was already applied for a two-dimensional case, e.g. by Long and Hargrave [7]. Although it is an important property for numerical simulations using level-set [16] and flame surface density methodologies [17], experimental data are scarce and only available in two dimensions [2,3,7]. The flame stretch rate can be determined as a combination of strain rate, curvature and flame displacement speed. This quantity is used in numerical simulations involving the flame surface density transport equation [18] or G-equation [19], where the assumption is made that all elements of the flame front are propagating at the laminar burning velocity appropriate to the local flame stretch rate. This premixed flamelet assumption has never been rigorously tested because of the absence of three-dimensional data. Throughout the paper we will compare data from our four plane measurement technique with data that were only evaluated from a single or two perpendicular planes. This is done in order to demonstrate the advantages of applying a three dimensional technique compared to planar measurements but also in order to show cases where a simpler, crossed plane arrangement can lead to sufficiently accurate results. 2. Experimental setup 2.1. Optical diagnostics setup Figure 1 shows a schematic of the optical setup for the crossed plane experiment. The four light sheets (LS) intersect at a common

line near the middle axis of the diffuser burner with an angle of 45° between each of them. To create the light sheets, three of the four lines of a Thales multi-channel Nd:YAG laser cluster are utilized in double-pulse operation with a time separation of 150 ls. Light sheet 1 is created from the frequency tripled emission of one line of the multi-Nd:YAG system at a wavelength of 355 nm. The respective camera is equipped with a 105 mm f/2.8 UV camera lens to capture the signal and a Schott UG11 filter to block interfering signals from the other three light sheets. Cameras 2 and 3 are equipped with Nikon 105 mm f/2.0 and Sigma 105 mm f/2.8 camera lenses, respectively, and band-pass filters for 532 nm to collect the scattered light from the corresponding light sheets. The polarizations of the two 532 nm laser beams are set perpendicular to each other and the signals from the two sheets are separated through different polarization filters. Finally, the last light sheet has an emission wavelength of 563 nm, produced by an external Raman laser [20] pumped with one of the frequency-doubled lines and then separated from the pump beam by a dichroic mirror. The signal is captured by camera 4, which is equipped with a band-pass filter for 560 nm. All cameras are equipped with Scheimfplug mounts to facilitate focusing the cameras when looking from a slight angle onto the light sheets. This is not depicted in Fig. 1 for clarity reasons. For the detection of the Mie scattering from light sheets 1 and 3, Lavision Imager ProX 4 M (4 megapixel resolution, 24.7 mm  24.7 mm field of view) cameras have been used. Light sheets 2 and 4 have been captured using Lavision Imager Intense cameras (1 megapixel resolution, 24.7 mm  20.2 mm field of view). To identify the position of intersection in each of the camera images, the four laser light sheets were aligned to a suspended thin filament with an attached weight. The measurements were taken with the aligned light sheets and after that, images of the filament were taken with each of the cameras. These were treated with the same dewarping algorithms as the measurement data of the flames to guarantee the precise determination of the position of intersection. 2.2. Flame geometry The investigated diffuser burner has been previously characterized using planar techniques [21,22]. It is made of quartz glass and has a diffuser half-angle of 12.6° from an inlet diameter of 58 mm. It was chosen for this study because of the perpendicularity of the flame surface to the mean flow direction which can be achieved if a high momentum annular wall jet is used to avoid separation from the diffuser wall. Aligning the four light sheets perpendicular to the flame surface yields highly accurate flame displacement measurements. Furthermore this arrangement allows the long interframing time required, while not causing too many out-of-plane losses of seeding particles. Because of the decrease of the mean flow velocity in the quartz glass diffuser, the flame stabilized at a specific height due to the instantaneous equilibrium of turbulent flame speed and mean flow velocity. For turbulence generation, a perforated plate with a blockage ratio of 37% was installed 25 mm upstream of the burner nozzle. 2.3. Image processing

Fig. 1. Schematic of the optical setup (Cam: Camera, LS: light sheet, LSFO: light sheet forming optics, BPx: band pass filter for wavelength x, PF: polarization filter, UG11: Schott UG11 filter, k/2: half wave plate). Scheimpflug mounts were used for all cameras but are not depicted here.

Processing of the raw data is a crucial step in the analysis of the results. The raw images of the four cameras were dewarped into their respective imaging planes by using a dot-target. After this step, particle image velocimetry software was used to calculate the velocity fields. Then, the dewarped raw images were filtered and utilized to extract the flame front location from each image. The information from each of the light sheets was then combined

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

to yield quantities such as the three-dimensional flame surface, the three-dimensional curvature, and the flame displacement speed. 2.3.1. Particle image velocimetry The dewarped particle images were fed into the commercial software DaVis from Lavision to calculate the velocity fields. Here, a multi-stage cross-correlation algorithm (see e.g., [23]) was used in order to gain high resolution vector fields. For cameras 1 and 3, interrogation window sizes of 128  128 pixels down to 32  32 pixels were used. The respective interrogation window sizes for cameras 2 and 4 were 64  64 to 16  16. All velocity calculations were done with a 50% overlap of the interrogation windows. This yielded very good results despite the high interframing times necessary for the determination of the flame displacement. It should be noted here that with an increase in interframing time, the flame movement can be determined more accurately due to a larger pixel displacement. But if the interframing time is chosen too long, the vector yield and the accuracy of vector determination is decreased due to a loss of seeding particles via out-of-plane motion. Hence an optimum needs to be determined for the given experimental conditions which was found here to be 150 ls. For our experiment, the seeding particles move by a maximum of 7.5 + 22.5 lm into the radial direction, which guarantees that around 70% of the particles stay in a 100 lm light sheet, resulting in a high vector yield of more than 97%. The in-plane displacement of droplets near the flame front was around 4.4 pixels for cameras 1 and 3 and 2.7 pixels for cameras 2 and 4, respectively. The resulting vector spacing is 192 lm for cameras 1 and 3 and 157 lm for cameras 2 and 4. A total of 126  126 vectors for cameras 1 and 3 and 159  130 vectors for cameras 2 and 4 have been calculated in this way. To remove erroneous vectors, groups of five or less vectors were removed. This filter removes patches of vectors in the burnt region of the image that are falsely found by the algorithm due to camera noise. Additionally, a filter removed all vectors which deviated by more than twice the rootmean-squared-value of the neighboring vectors. Finally, a 3  3 linear smoothing filter was applied. To evaluate the data at specific positions, linear interpolation of the vector information was used. The three-dimensional velocity components at the intersection line were determined by a linear combination of the information from the four light sheets. The axial component was calculated by a mean of all four light sheet’s axial component of velocity. For the radial components, values from light sheets 1 and 3 were used as reference, because they coincide with the x- and y-axes of the selected coordinate system. The radial components from light sheets 2 and 4 were then projected onto them via a coordinate transformation and the mean of the respective radial components was calculated. 2.3.2. Image filtering The possibility of extracting the flame front position from raw particle images has been used in several studies. This can be done using the step in particle number density when non-combustible ceramic particles are used, see e.g., Pfadler et al. [15,24] or Steinberg et al. [25,26], or from the disappearing of oil droplets, which evaporate close to the flame front [27]. Hereby the accuracy with which the flame front can be located depends more on the particle seeding density than on the image resolution. In the present work, a combination of filtering algorithms and contour extraction is used to determine the flame front position and characteristics with sub-pixel accuracy; the procedure is described in the following. A local standard deviation filter is used to remove any flame luminescence from the images. This is especially important for the second frame which has an inherently longer exposure for interline transfer CCD cameras. In a second processing step, peak

2759

maxima are equalized. This is accomplished by an analysis of the image’s histogram. The maximum intensity threshold is set to the intensity at which the number of pixels showing this value goes down to only a hundredth of the maximum number of pixels showing the same value. To illustrate this, a schematic of a histogram is shown in Fig. 2. Although only around 2% of the image pixels exceed this threshold, their removal still leads to a significant homogenization of the high intensity part of the image. The next step consists of a nonlinear sliding average filtering of the image to smooth the flame contour, avoiding the detection of artificial high frequency wrinkles produced by the particle pattern. To prepare the image for further processing, a background subtraction and normalization is applied, using histogram analysis again. A non-linear diffusion filter algorithm published by Malm et al. [28], using a diffusivity function originally proposed by Weickert [29], was then applied to achieve the goal of enhancing the flame front gradients and smoothing other gradients, e.g. from particles. This leads to a step function at the flame front position. After this, another non-linear sliding average filter flattens any remaining peaks near the flame front that are disregarded by the previous algorithm. After this, a Gaussian convolution of the image was used to extract the gradient map of the image. This was fed into an algorithm proposed by Steger [30] for the sub-pixel accurate detection of curvilinear structures. This algorithm has been used in combustion diagnostics by Sweeney and Hochgreb [31] to improve the results from a Canny algorithm [32] but without utilizing the advantage of sub-pixel accurate line detection. With this algorithm, the results are not discrete anymore, in contrast to the Canny algorithm, and therefore do not artificially produce the high curvature values that usually need to be smoothed by applying a polynomial fit function (see [33]). The result is parameterized contour information for each of the two frames of the four camera images. Figure 3 shows an example of a raw image with superimposed extracted contour information. In the image one can see that the contrast between burnt and unburnt part is very low, but still the flame contour is extracted reliably (only about 1% of the images have been discarded due to flame contours which were unrecognizable by the software). This image was chosen intentionally to show the worst contrast that occurred in the raw images of all cameras. 2.3.3. Three-dimensional evaluations The information from each of the light sheets is then combined in the following way to yield flame curvature, flame displacement speed, strain and stretch rate in three dimensions.

Fig. 2. Schematic of an image histogram illustrating the procedure to find the threshold intensity for the homogenization of the high intensity part of the image.

2760

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

sh ¼

Fig. 3. Example dewarped raw image with extracted contour superimposed. Also shown are excerpts of the intensity information from a horizontal and vertical line depicted in red. A 24.7 mm by 20.2 mm area of the second frame of light sheet 4 is shown. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The three-dimensional mean curvature jm is defined as the mean of curvatures j1 and j2 of the principal directions of the given surface, see the following equation:

jm ¼

ðj1 þ j2 Þ 2

ð1Þ

j1 and j2 are determined by calculating their reciprocals, the principal radii R1 and R2, by solving the following equation[34]: 4

ðrt  s2 ÞR2 þ h½2pqs  ð1  p2 Þt  ð1  q2 ÞrR þ h ¼ 0

ð2Þ

With



@z ; @x



@z ; @y



@2z ; @x2



@2z ; @x@y



@2z @y2

ð3Þ

where z denotes the downstream position of the flame and x and y are the horizontal axes of planes 1 and 3, which are perpendicular to each other and the intersection line. Additionally,

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h ¼ 1 þ p2 þ q2

ð4Þ

In the formulation of the variable s in Eq. (3), the requirement of four light sheets at 45° to each other (see Fig. 1) is evidenced by the combined derivative in x and y which could not be calculated from two-plane data alone. The resulting equation for the mean curvature is given as

jm ¼

ðj1 þ j2 Þ rð1 þ q2 Þ  2pqs þ tð1 þ p2 Þ ¼ 3 2 2ð1 þ p2 þ q2 Þ2

ð5Þ

The data from the four measurement planes can be combined to yield the discretized function z(x, y), which can be used in Eqs. (3) and (4) to calculate the derivatives. It is obvious that the flame surface must not be parallel to one or more of the measurement planes in order to avoid ambiguous descriptions of z(x, y). The first and second order derivatives in the above equations were calculated via the central differences methods, with a spacing of 0.5 mm between the data points. To evaluate the data at these specific positions, linear interpolation was used. To analyze the flame shape, a shape factor can be defined [35,36]

j1 ; with jj1 j < jj2 j j2

ð6Þ

This value gives information on the three-dimensional structure of the surface, being greater than 0 for elliptic, ‘‘sphere-like’’ and smaller than 0 for hyperbolic, ‘‘saddle-like’’ structures. For a value of 0, the flame is cylindrically shaped (parabolic). To illustrate these different flame shapes, schematics are given in Fig. 4. This is a crucial distinction as it can give a hint as to whether a twodimensional projection can be viable for the flame in question. For mainly hyperbolic structures, it can be doubted that a projection will correctly capture the real three-dimensional structure of the flame surface. Knowing the flow field, the three-dimensional curvatures and the respective normal vector to the flame surface at the intersection, the next possible step is to determine the three-dimensional flame displacement speed. This is done by translating the intersection point of the flame front in the first frame by the velocity measured 0.4 mm upstream in the unburnt, as it is not possible to evaluate the velocity exactly at the intersection point due to a lack of seeding droplets in the burnt region. A side effect of this choice of location is that the flame displacement speed does not need to be corrected for reduced gas densities, which would be the case if the velocities were taken at a location closer to the flame front or even in the burnt region [2,10]. Additionally, thermophoretic effects negatively influence the accuracy of PIV near the flame front [37]. Balusamy et al. [38] investigate the behavior of the fluid near the flame front further. The suggested methods should be taken into account if velocity information closer to the flame is to be gained. Following the flame surface normal from this advected virtual point, one can calculate the distance to the flame surface of the second frame by applying a cubic interpolation spline to it and solving the set of linear equations leading to the distance. The procedure is explained for a two-dimensional measurement in [7]. The accuracy of determining the flame front position is assumed to be 0.25 pixels due to the sub-pixel nature of the information gained. It is furthermore assumed to be the biggest factor of uncertainty. Thus, the error in determining the flame displacement speed is estimated to be around 10%. The flame stretch rate is defined as the temporal rate of change of flame area per unit area [39]

s_ ¼

1 dA A dt

ð7Þ

It can also be related to the strain, local flame displacement speed and curvature [39,40]

s_ ¼ at þ 2sD jm

ð8Þ

jm being the mean of the principal curvatures of the flame and sD the local flame displacement speed. at is described as the surface divergence due to the velocity gradients tangent to the flame front. It is defined by [40,41] as at ¼ ðdij  ni nj Þ

@ui @xj

ð9Þ

Here ui and xi are the velocity and flame normal vector component, respectively. dij denotes the delta function, and ni and nj are the unit normal vector components of the flame surface. This information was extracted from the vector fields of all four planes. The flame stretch rate is therefore a combination of flow field and reaction information which can only be obtained by simultaneously measuring flame structure and flow field in a three-dimensional manner. In the graphs shown in Section 3, these quantities are sometimes normalized. This is done for the curvature by multiplication with the laminar flame thickness dL calculated as the ratio of the thermal diffusivity a evaluated at the mean flame temperature

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

2761

Fig. 4. Schematics of different possible flame shapes. From left to right: hyperbolic (sh < 0), parabolic (sh = 0), elliptic (sh > 0).

[42] for the given stoichiometry and the laminar burning velocity sL [43]. For the operating condition here, these values were sL = 0.35 m/s, a = 1.7  104 m2/s, and thus dL = 0.4 mm. Stretch rates are normalized by multiplying by a timescale, dL/sL = 1.1 ms, to form one possible definition of the Karlovitz number:

Ka  s_

dL a  s_ 2 sL sL

ð10Þ

The flame displacement speeds are normalized by the laminar burning velocity. They may be analyzed in terms of the components due to the reaction, due to normal diffusion, and due to tangential diffusion [1]:

sD  sr þ sn þ st

ð11Þ

where st  2Djm , with D being the mass diffusivity. Thus, the term sr + sn can also be analyzed by a subtraction of st from sD. The normalized expression for Eq. (11) under the assumption of unity Lewis number reads

sD sr þ sn ¼  2dL jm sL sL

Fig. 5. Axial velocity (s), fluctuation () and turbulence intensity (h) profiles along the intersection line of the four light sheets (for clarity reasons, only every 10th data point is represented by a marker).

ð12Þ

The value for sD is measured directly in the way described above (and, in more detail, in Ref. [7]) and jm is determined via Eq. (5) and the measured values for z(x, y). With this information, st and the term sr + sn can be determined through Eq. (12). 3. Discussion of results 3.1. Isothermal measurements The isothermal measurements were done with an outlet bulk velocity of 4 m/s and wall jet velocity of 6 m/s, resulting in a momentum ratio of 0.167. These flow conditions match the ones used in the later reactive measurements that are discussed in Section 3.2. In Fig. 5 the axial velocity along the centerline for all four light sheets is presented, along with the root-mean-square fluctuations and the resulting turbulence intensity. The graph shows a very good agreement of the axial velocity component in all four light sheets, indicating that the light sheet calibration procedure is working properly. It can be seen that the mean axial flow rate decreases downstream as expected due to the area expansion of the glass diffuser [44], and some boundary layer growth. The area increase of 23% between 53 mm and 73 mm accounts almost exactly for the measured axial velocity decrease. At the same time, the fluctuations in the flow are slightly increasing. Together this leads to a rise in turbulence intensity within the measured region from 0.14 to 0.24, with a value of 0.18 at the position that turns out to be the mean flame height. This behavior of turbulence intensities is typical for a conical diffuser [44].

A comparison of the axial velocity component in the radial direction at the lower end of the measured region, 53 mm downstream of the burner exit, can be seen in Fig. 6. Here, the difficulty in homogenizing the flow can be seen. Although the resulting profiles are not perfectly flat, the deviations are only around 10%. Given the fact that data are only gained close to or directly at the intersection position, this deviation can be neglected. The root-mean-square fluctuations of the axial velocity component are similar in all four light sheets, indicating a rather homogeneous turbulence distribution. It can be seen from Fig. 7 that the fluctuations of the radial velocity component, which is depicted together with the radial mean velocity, are only half as great as the fluctuations of the axial component. This indicates an anisotropic turbulence behavior of the flow, as can be expected in such a diffuser configuration [44], which is similar to the turbulence profile in a pipe flow. As the radial components are only projections of the real flow velocity into the respective investigated plane, it is not expected to see them precisely zero on the axis. The values for the radial velocity are only 3% of the axial component and calculation based on continuity shows that the gradient on the axis is such that only a 2 mm asymmetry is needed to account for this. The root-meansquare fluctuations show good agreement for all four investigated light sheets. 3.2. Reactive measurements A planar premixed methane air flame was stabilized within the diffuser burner. A depiction of it is given in Fig. 8. The operating conditions were chosen to be / = 0.95 with a bulk velocity of

2762

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

Fig. 6. Axial velocity (s) and fluctuation () profiles for all four light sheets at a height of 53 mm above the burner exit (for clarity reasons, only every 10th data point is represented by a marker in the right figure).

Fig. 7. Radial velocity (s) and fluctuation () profiles along the centerline in all four light sheets (for clarity reasons, only every 10th data point is represented by a marker).

4 m/s. For these conditions, the combustion of the seeded silicone oil droplets contributes only 0.5% to the total heat release of the flame. As has been seen in Fig. 5 turbulence intensities of around 20% are generated in the investigated flame stabilization region at a height of around 63 mm above the burner exit. This places the flame in the ‘corrugated’ regime of the Borghi diagram [22]. The Karlovitz and Damkoehler numbers are 0.3 and 8.1, respectively, calculated from the definitions given in [42]. 3.2.1. Three-dimensional curvatures The probability density function (pdf) of the three-dimensional curvature is compared with the two-dimensional curvature measured in one plane in Fig. 9. The two-dimensional curvatures were extracted from light sheet 3 (see Fig. 1). Curvatures statistics from the other light sheets are virtually identical. A set of 1104 recordings was used for the compilation of the probability density functions. Both pdfs show a peak at slightly positive values, but the twodimensional curvatures form a broader probability density function. Other experimental investigations [5,6,45] have shown such a behavior for two-dimensional curvatures of flat flames, too, where the probability density function is centered around zero with a skew to the positive. To illustrate the relationship of twoand three-dimensional curvature values, Fig. 10 shows a comparison in the form of a scatter plot with one dot per measurement.

Furthermore, a distinction according to the shape of the threedimensional flame surface is shown in differently colored points, where blue1 depicts an elliptic (sphere-like) and red a hyperbolic (saddle-like) shape of the flame front. A correlation between two- and three-dimensional curvatures is evident but the scatter is significant, indicating the inability of the planar measurements to capture the 3D curvature. It is also notable that for elliptic shapes the correlation seems better, as they occur mainly near the intersection line, whereas hyperbolic shapes mainly occur off the bisecting line. A linear regression of the conditioned values shows the same slope, but the root-mean square-error of the regression is with a value of 0.096 for hyperbolic shapes nearly double the value for elliptic shapes, which is 0.052. This behavior is quite obvious, as for elliptic shapes the two principal curvatures have the same sign. As the two principal curvatures are by definition the highest and lowest curvatures of the surface, every twodimensional curvature chosen from a random plane through the same point needs to have the sign of the principal curvatures. This leads to a rather good positive correlation. On the other hand, for hyperbolic shapes, the principal curvatures have different sign and randomly chosen curvatures through the same point can have either positive or negative sign, which can lead to a negative correlation or no correlation at all. In order to improve curvature measurements, a crossed plane setup involving two planes has been be used in the past. This approach yields ‘‘quasi-three-dimensional’’ data in the form of the mean curvature between the curvatures calculated in the two respective planes. This is shown in Fig. 11, where a comparison with the true three-dimensional curvature from four planes is made. The information for the two perpendicular planes was taken from light sheets 1 and 3 (see Fig. 1). A very good correlation can be observed between the two approaches and there is only limited scatter which indicates that a less complex set-up of only two crossed light sheets might already yield a very good estimate of the three-dimensional curvature. The distinction between elliptic and hyperbolic shapes in the plot only signifies that hyperbolic shapes mainly occur for curvatures around zero, whereas elliptic shapes are seen more often for higher values of curvature. The flame surface can be further investigated by considering the shape factor of curvatures. For this, the pdf of said shape factor is plotted in Fig. 12. Positive values correspond to spherically shaped elliptic structures, whereas negative values correspond to hyper-

1 For interpretation of color in Fig. 10, the reader is referred to the web version of this article.

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

2763

Fig. 8. Photograph (left) and short exposure chemiluminescence image (right) of the flame at operating conditions (/ = 0.95, ubulk = 4 m/s).

Fig. 9. Probability density function of three-dimensional curvature calculated from the four light sheets compared to two-dimensional curvature evaluated from light sheet 3.

Fig. 10. Normalized three-dimensional curvature calculated from four light sheets in comparison to the normalized two-dimensional curvature in light sheet 3.

bolic saddle structures. Values of zero are calculated for cylindrical structures. It can be seen that there is a greater probability of finding spherically-shaped flame fronts than hyperbolically-shaped ones

Fig. 11. Comparison of normalized three-dimensional curvatures calculated from four planes (from Eq. (5)) versus the normalized mean value of the curvatures from two perpendicular planes (light sheets 1 and 3).

in this burner. A scatter plot of the shape factor, comparing it to the mean curvature, is shown in Fig. 13. The trend to positive, spherically-shaped curvatures can also be observed here. Elliptic shapes can be observed throughout the spectrum of mean curvatures, whereas hyperbolically-shaped flame fronts are scarce for high absolute mean curvatures and mainly occur from 0.25 mm1 to 0.25 mm1. This was already observed in previous plots and can be explained by the definition of both shape factor and mean curvature. For a hyperbolic structure, the shape factor is smaller than zero and the principal curvatures therefore have different signs, resulting in low mean curvatures. To achieve high mean curvatures on the other hand, the principal curvatures need to be high and of the same order of magnitude, which leads to a positive shape factor representing elliptic shapes. It could be that the nature of the burner system leads to such a distribution, favoring elliptic shapes, because of the tendency to mean curvature of the flame. Comparison with the results of DNS calculations [1], where the peak of the pdf was at zero indicating cylindrical structures, reinforces this conclusion. We can also compare the shape factor derived from only two crossed planes with the ‘‘true’’ shape factor from the four plane

2764

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

Fig. 12. Probability density function of the shape factor evaluated from four planes.

Fig. 13. Scatter plot of the shape factor versus the three-dimensional curvature evaluated from four planes.

measurements. For the two-plane case, the shape factor is not calculated by the principal curvatures, but by a set of curvatures measured in perpendicular, but arbitrary planes. Here, light sheets 1 and 3 are chosen again. A comparison to the real shape factor is given in Fig. 14. Although a clear correlation of the two values can be seen, it is also evident that the scatter increases for high absolute values of shape factor. Analyzing the data conditioned by the mean cur-

Fig. 14. Comparison of the shape factor calculated from four planes and the shape factor determined from only two perpendicular planes (light sheets 1 and 3).

vatures gives two main results. One of them can be seen already from Fig. 13, i.e. that for high mean curvatures mainly elliptic shapes are present, and for negative shape factors, one finds mainly small mean curvatures around zero. But for these small mean curvatures the scatter in the shape factor plot is the highest. This can be explained by the vicinity of the two-dimensional curvatures to zero. A slight change in orientation of the coordinate system can lead to a change in sign for one of the curvatures, which results in an anti-correlating behavior of the shape factor. Bearing this in mind, it can also be explained why there are no strongly anti-correlating points for high positive values of shape factor. As these values are mainly corresponding to high absolute mean curvatures and both principal curvatures are of comparable magnitude, a change in orientation does not result in a change in sign so that the shape factor is more or less correlated. This suggests that use of a set-up with only two crossed planes will lead to unreliable information on the shape of the flame, especially for low absolute mean curvatures. Another interesting structural feature is the orientation of the flame surface. As the flame is propagating perpendicular to its surface, the orientation of the said surface determines the main propagation direction. To investigate this, the unit flame normal vector is split up into its components N1, N2 and N3, N1 being the component parallel to the burner axis and N2 the component aligned with the x-axis, which is found in light sheet 1 (c.f. Fig. 1). N3 corresponds to the component aligned with the y-axis and therefore is within light sheet 3. Figure 15 shows the probability density functions of ðN21 þ N22 Þ and ðN21 þ N23 Þ. As both distributions involving N1 show a high probability around 1, a predominant propagation of the flame front in the mean flow direction is evident. This also means that it is well aligned with the intersection line of the four light sheets leading to the maximum possible sensitivity for flame displacement measurements, which are evaluated in Section 3.2.2. However, the two probability density functions do not perfectly match, which indicates a slight misalignment of the flame normal to the centerline of the burner. This is also found and discussed in Fig. 16. The orientation and propagation direction of the flame can be further investigated by creating scatter plots for the relations between N1, N2 and N3. These can be found in Fig. 16. As observed before, the predominant propagation direction for the flame front is along the N1 direction opposite to the flow, with arbitrary scatter in N2 and N3. The few positive sample points for N1 indicate a propagation of the flame in the mean flow direction. These are directly correlated to the cases of negative flame displacement speeds which are discussed later on in the context of Figs. 18 and 21. Given that the scatter cloud for N2 and N3 is not completely homogeneous, it can be deduced that the flame surface orientation was not perfectly parallel to the intersection line. This behavior can also be seen in DNS simulations of statistically planar flames, e.g., in [1]. But although the main propagation direction of the flame is aligned with the axis of the burner, for the comparison of twodimensional and three-dimensional flame displacement speeds, one other parameter needs to be found to describe whether the flame is mainly propagating within the investigated light sheet or whether the out-of-plane component is significant. Therefore, a probability density function of a- defined as the angle between the two-dimensional plane investigated, which is light sheet 3 in this case, and the flame normal vector – is shown in Fig. 17. One can see that, in good agreement with the data found in [1] for statistically planar flames, the propagation of the flame does align well with the detection plane. A plot of the angle for the orthogonal plane shows a similar trend but is not depicted here. Although the main propagation direction is aligned with the axis of the burner, the movement of the flame is significant and

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

2765

Fig. 15. Probability density function of ðN21 þ N22 Þ and ðN21 þ N23 Þ showing an orientation of the flame surface normal along the centerline.

Fig. 16. Scatter plots of the relations between N1, N2, and N3 showing a predominant propagation of the flame against the flow.

therefore a projection of three-dimensional displacement speed into one plane is expected to give higher values of displacement speed. This will be further analyzed in the next section. 3.2.2. Determining the flame displacement speed As for the curvatures in Fig. 10, a comparison of the threedimensional flame displacement speed with the two-dimensional one calculated from light sheet 3 can be found in Fig. 18. A positive correlation can be observed but nonetheless there is a large scatter, which indicates the unreliability of the reconstruction of threedimensional flame displacement speeds from two-dimensional data.

This point can be further corroborated by the comparison of probability density functions from three- and two-dimensional flame displacement speeds, shown in Fig. 19. The probability of finding high values is larger in the two-dimensional case and the distribution for the three-dimensional data is narrower. The mean value of the normalized three-dimensional flame displacement speed is 1.1 and therefore the mean displacement speed is indeed approximately equal to the laminar flame speed, but there is a very broad distribution with a standard deviation of 0.71 in sD/sL. This deviation is very similar to those in the simulations of [1], and this figure constitutes a major result of this investigation, which will be discussed further in Section 3.2.3.

2766

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

Fig. 17. Probability density function of the angle a between the two-dimensional plane of light sheet 3 and the flame surface normal.

Fig. 19. Probability density function of the normalized three-dimensional flame displacement speed calculated from four planes compared to the normalized twodimensional flame displacement speed evaluated from light sheet 3.

Fig. 18. Normalized three-dimensional flame displacement speed calculated from four planes in comparison to the normalized two-dimensional flame displacement speed evaluated from light sheet 3.

Fig. 20. Probability density function of the normalized three-dimensional flame displacement speed calculated from four planes compared to the angle corrected normalized two-dimensional flame displacement speed evaluated with information from light sheets 1 and 3 (2D corrected).

The two-dimensional mean value, however, shows a shift to higher values of flame displacement speed. Bearing in mind the findings from Fig. 17 it can be explained by the misalignment of the flame propagation direction with the investigated light sheet. This can be demonstrated by correcting the two dimensional values with the angle a and the out-of-plane component u3 of the velocity by using the following equation[2]:

  s2Dcorr ¼ s2D D D  u3 tan a  cos a

ð13Þ

The resulting probability density function is given in Fig. 20. With this correction applied, the statistics of two-dimensional and three-dimensional flame displacement speeds do overlap almost perfectly. This shows that even with two orthogonal planes, the statistics of three-dimensional flame displacement speed can be derived satisfactorily by the two-dimensional flame speed determined from one plane and the angle a and the out-of-plane component of the velocity determined in the orthogonal plane. To investigate the point-by-point correlation of the threedimensional flame displacement speed and the corrected twodimensional one, a scatter plot of the two values, similar to Fig. 18, is shown in Fig. 21. It can be clearly seen that even if the statistics of three-dimensional flame displacement speed can be deduced from the twodimensional data in two crossed planes, the individual values still

do not correlate perfectly. This makes three-dimensional measurements essential in order to yield data of significance. A decomposition of the flame displacement speed into its tangential-diffusion (curvature) component and the combined normal-diffusion and reaction component (derived by subtracting the curvature component from the total in Eq. (12)) has been shown to be important in the context of level-set modeling [17,46]. It yields the pdfs in Figs. 22 and 23, respectively. In the plot of the tangential-diffusion component, predominantly negative values can be found. This can be explained by the definition of the component and the fact that the mean curvature of the investigated flame is positive. The distribution is therefore more or less a mirrored graph of the mean curvature histogram in Fig. 9. Thus, the two-dimensional distribution is here also broader than the three-dimensional one, but the means of both are similar, with st/sL = 0.1. The pdf of the combined normal-diffusion and reaction components shows a similar trend to the original flame displacement speed pdf, but there is a smaller chance of finding values in the negative part of the histogram. This leads to the conclusion that the tangential-diffusion component is mainly responsible for displacement speeds in the same direction as the mean flow, whereas the combined normal-diffusion and reaction components tend to

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

Fig. 21. Normalized three-dimensional flame displacement speed calculated from four planes in comparison to the angle corrected normalized two-dimensional flame displacement speed evaluated with information only from light sheets 1 and 3.

Fig. 22. Probability density function of the tangential diffusion component of the normalized three-dimensional flame displacement speed calculated from four planes compared to the tangential diffusion component of the normalized twodimensional flame displacement speed evaluated from light sheet 3.

move the flame against the flow direction. Similarly to Fig. 19, the two-dimensional distribution is overestimating the combined normal-diffusion and reaction components of the flame displacement speed.

3.2.3. Flame stretch relations It is now of interest to see whether these deviations in the flame displacement speed from the unstretched laminar burning velocity can be related to the flame stretch. To analyze this, first the relations of the quantities contributing to the flame stretch have been plotted against the flame stretch. A scatter plot of the normalized flame stretch, the Karlovitz number, versus the normalized mean curvature of the flame can be seen as Fig. 24. A positive correlation can clearly be seen for both hyperbolic and elliptic flame structures, but it is in combination with the flame displacement speed that the contribution to flame stretch is made (Eq. (8)). The contribution of the tangential component of strain rate can be seen in Fig. 25. A positive correlation can be seen here too. Furthermore, negative values for tangential strain rate are relatively scarce. The strain rate was calculated 0.4 mm upstream of the flame position. The distance is about the same as the laminar flame

2767

Fig. 23. Probability density function of combined normal-diffusion and reaction components of the normalized three-dimensional flame displacement speed calculated from four planes in comparison to the combined normal-diffusion and reaction components of the normalized two-dimensional flame displacement speed evaluated from light sheet 3.

Fig. 24. Normalized flame stretch versus normalized mean curvature, both calculated from four planes.

thickness, and about a fifth of the smallest radius of flame curvature. This was chosen to be close enough to the flame to be of relevance but far enough to have sufficient unburnt velocity data in the vicinity. In summary, Figs. 24 and 25 show that the flame stretch rate is almost always positive for the measured condition which indicates that the flame is constantly stretched, not only by the convective contribution of the flow field but also by the chemical reaction. From the slope of the positive correlation with strain rate, it can be concluded that about 80% of the flame stretch is due to motion of the fluid. The flame displacement speed is displayed as a function of the flame stretch in Fig. 26, but no substantial correlation can be seen. Even more surprisingly in view of the DNS results in Ref. [1], there was no correlation between the flame displacement speed and the three dimensional curvature or the strain rate when these contributions to flame stretch were isolated and examined independently. To determine the reason for the variability in the displacement speed may require data for the small-scale turbulence characteristics in the preheat zone of the flame. Such turbulence, while not perturbing the shape of the flame front, might have a substantial influence on the diffusion of species and heat

2768

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769

However, a simplification of the measurement technique using only two crossed planes was explored and found to deliver sufficiently accurate data for several quantities in comparison to the fully three-dimensional ones. With the instantaneous local Karlovitz numbers seldom greater than unity, the investigated flame was in the corrugated regime. Speeds ranging from 0.4 sL to 4.0 sL with a mean of 1.1 sL were recorded in a burner designed to generate a mean turbulent flame front that was planar. These major deviations from the unstretched laminar burning velocity bore no relationship to the flame stretch, nor to the contributions to stretch due to strain or curvature. These are believed to be the first simultaneous measurements of the three-dimensional flame shape, curvature and displacement speed, and refinement of the technique to yield velocity data closer to the flame front should be possible. Acknowledgments Fig. 25. Normalized flame stretch versus normalized strain rate tangential to the flame surface, both calculated from four planes.

The authors gratefully acknowledge the invaluable help from Stephen Ramsey with the manufacturing of the glass cone. References

Fig. 26. Normalized flame stretch versus normalized flame displacement speed, both calculated from four planes.

in its vicinity, and hence on its rate of propagation. Gulder and Smallwood [47] discuss the possibility that transport of heat and species by turbulence with scales smaller than the laminar flame front thickness may have a fundamental effect on the local propagation rate of the flame. However, such measurements clearly present a considerable challenge. 4. Conclusions A quad plane particle image velocimetry experiment for the three dimensional investigation of turbulent premixed flame structures was presented here. The described technique not only gives access to the flame displacement speed, but also to flow characteristics such as the strain rate. A combination of these can yield the flame stretch rate or a decomposition of the flame displacement speed into its tangential-diffusion and combined normal-diffusion and reaction components. These very relevant quantities can all be gathered three dimensionally, and not merely in a two-dimensional light sheet, as in many earlier studies. Thus, a comparison of two- and three-dimensional data was possible, showing that for certain flame geometries, two-dimensional data may be sufficient to represent the processes occurring. But in most real flames, a two-dimensional measurement cannot reproduce the three-dimensional information necessary for model validation or a deeper understanding of the underlying processes.

[1] N. Chakraborty, G. Hartung, M. Katragadda, C.F. Kaminski, Combust. Flame 158 (7) (2011) 1372–1390. [2] G. Hartung, J. Hult, R. Balachandran, M.R. Mackley, C.F. Kaminski, Appl. Phys. B: Lasers Opt. 96 (4) (2009) 843–862. [3] B. Renou, A. Boukhalfa, D. Puechberty, M. Trinite, in: Twenty-Seventh Symposium (International) on Combustion, vols. 1 and 2, 1998, pp. 841–847. [4] K.H.H. Goh, P. Geipel, F. Hampp, R.P. Lindstedt, Proc. Combust. Inst. 34 (2013) 3311–3318. [5] D.A. Knaus, F.C. Gouldin, Proc. Combust. Inst. 28 (2000) 367–373. [6] D.A. Knaus, F.C. Gouldin, D.C. Bingham, Combust. Sci. Technol. 174 (1) (2002) 101. [7] E.J. Long, G.K. Hargrave, Flow Turbul. Combust. 86 (3–4) (2011) 455–476. [8] T.D. Upton, D.D. Verhoeven, D.E. Hudgins, Exp. Fluids 50 (1) (2011) 125–134. [9] K. Yamamoto, S. Isii, M. Ohnishi, Proc. Combust. Inst. 33 (2011) 1285–1292. [10] P.J. Trunk, I. Boxx, C. Heeger, W. Meier, B. Böhm, A. Dreizler, Proc. Combust. Inst. 34 (2) (2013) 3565–3572. [11] M. Shimura, T. Ueda, G.-M. Choi, M. Tanahashi, T. Miyauchi, Proc. Combust. Inst. 33 (1) (2011) 775–782. [12] S. Pfadler, J. Kerl, F. Beyrau, A. Leipertz, A. Sadiki, J. Scheuerlein, F. Dinkelacker, Proc. Combust. Inst. 32 (2009) 1723–1730. [13] T. Landenfeld, A. Kremer, E.P. Hassel, J. Janicka, T. Schafer, J. Kazenwadel, C. Schulz, J. Wolfrum, in: Twenty-Seventh Symposium (International) on Combustion, vols. 1 and 2, 1998, pp. 1023–1029. [14] D.A. Knaus, S.S. Sattler, F.C. Gouldin, Combust. Flame 141 (3) (2005) 253–270. [15] S. Pfadler, F. Beyrau, A. Leipertz, Opt. Express 15 (2007) 15444–15456. [16] N. Peters, Turbulent Combustion, Cambridge University Press, Cambridge, 2000. [17] N. Chakraborty, R.S. Cant, Proc. Combust. Inst. 32 (2009) 1445–1453. [18] D. Veynante, L. Vervisch, Prog. Energy Combust. Sci. 28 (3) (2002) 193–266. [19] H. Pitsch, Annu. Rev. Fluid Mech. 38 (2006) 453–482. [20] J. Kerl, T. Sponfeldner, F. Beyrau, Combust. Flame 158 (2011) 1905–1907. [21] V.R. Savarianandam, C.J. Lawn, Combust. Flame 146 (1–2) (2006) 1–18. [22] C.J. Lawn, R.W. Schefer, Combust. Flame 146 (1–2) (2006) 180–199. [23] J. Westerweel, D. Dabiri, M. Gharib, Exp. Fluids 23 (1) (1997) 20–28. [24] S. Pfadler, F. Dinkelacker, F. Beyrau, A. Leipertz, Combust. Flame 156 (8) (2009) 1552–1564. [25] A.M. Steinberg, J.F. Driscoll, S.L. Ceccio, Exp. Fluids 47 (3) (2009) 527–547. [26] A.M. Steinberg, J.F. Driscoll, Combust. Flame 157 (7) (2010) 1422–1435. [27] I.G. Shepherd, E. Bourguignon, Y. Michou, I. Gokalp, in: Twenty-Seventh Symposium (International) on Combustion, vols. 1 and 2, 1998, pp. 909–916. [28] H. Malm, G. Sparr, J. Hult, C.F. Kaminski, J. Opt. Soc. Am. A – Opt. Image Sci. Vision 17 (12) (2000) 2148–2156. [29] J. Weickert, in: Proceedings of the Seventh European Conference on Mathematics in Industry, vol. 9, 1994, pp. 355–362. [30] C. Steger, IEEE Trans. Pattern Anal. Mach. Intell. 20 (2) (1998) 113–125. [31] M. Sweeney, S. Hochgreb, Appl. Opt. 48 (19) (2009) 3866–3877. [32] J. Canny, IEEE Trans. Pattern Anal. Mach. Intell. 8 (6) (1986) 679–698. [33] R.S.M. Chrystie, I.S. Burns, J. Hult, C.F. Kaminski, Meas. Sci. Technol. 19 (12) (2008). [34] I.N. Bronstein, K.A. Semendjajew, G. Musiol, H. Mühlig, Taschenbuch der Mathematik, Verlag Harri Deutsch, Frankfurt am Main, 2001. [35] C.J. Rutland, A. Trouve, Combust. Flame 94 (1–2) (1993) 41–57. [36] N. Chakraborty, R.S. Cant, Numer. Heat Transfer Part A – Appl. 50 (7) (2006) 623–643. [37] C.J. Sung, C.K. Law, R.L. Axelbaum, Combust. Sci. Technol. 99 (1–3) (1994) 119– 132.

J. Kerl et al. / Combustion and Flame 160 (2013) 2757–2769 [38] S. Balusamy, A. Cessou, B. Lecordier, Exp. Fluids 50 (4) (2011) 1109–1121. [39] L.W. Kostiuk, K.N.C. Bray, Combust. Sci. Technol. 95 (1–6) (1994) 193–212. [40] D. Bradley, A.K.C. Lau, M. Lawes, Philos. Trans. Roy. Soc. Lond. Ser. A – Math. Phys. Eng. Sci. 338 (1650) (1992) 359–387. [41] K.W. Jenkins, M. Klein, N. Chakraborty, R.S. Cant, Combust. Flame 145 (1–2) (2006) 415–434. [42] J.F. Driscoll, Prog. Energy Combust. Sci. 34 (1) (2008) 91–134.

2769

[43] T. Tahtouh, F. Halter, C. Mounaim-Rousselle, Combust. Flame 156 (9) (2009) 1735–1743. [44] P.A.C. Okwuobi, R.S. Azad, J. Fluid Mech. 57 (20) (1973) 603–622. [45] I.G. Shepherd, R.K. Cheng, T. Plessing, C. Kortschik, N. Peters, Proc. Combust. Inst. 29 (2002) 1833–1840. [46] E.R. Hawkes, J.H. Chen, Proc. Combust. Inst. 30 (2005) 647–655. [47] O.L. Gulder, G.J. Smallwood, Combust. Sci. Technol. 179 (1–2) (2007) 191–206.