Previous Article |
Table of Contents
| Next Article 
From The Cover
APPLIED MATHEMATICS
Numerical simulation of a laboratory-scale turbulent V-flame
*Center for Computational Science and Engineering and Environmental Energy Technologies Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720
Communicated by Phillip Colella, University of California, Berkeley, CA, May 19, 2005 (received for review March 15, 2005)
 |
Abstract |
|---|
 Top Abstract Experimental Configuration and...Computational ModelExperimental ComparisonsConclusionReferences |
|---|
adaptive mesh refinement | low-Mach-number flow | turbulent premixed combustion
Numerical simulation offers the potential to augment theory and experiment and overcome the limitations of standard approaches in analyzing laboratory-scale flames. The excessive computational costs of incorporating detailed transport and chemical kinetics have necessitated compromises in the fidelity or scope of simulations for premixed turbulent combustion. Simulation of laboratory-scale systems typically involves models for subgrid-scale turbulent fluctuations. Approaches based on large eddy simulation or Reynolds-averaged Navier-Stokes fall into this class. In addition to the turbulence model, these approaches require a model for the speed of flame propagation in a turbulent field or some other model for turbulence-chemistry interaction.
The goal of the present work is to simulate a laboratory-scale flame, 
10-12 cm in length, without a turbulence model, a burning velocity model, or the introduction of some other type of turbulence closure hypothesis. Standard computational tools for this type of simulation are based on high-order (more than four) explicit integration methods for the compressible Navier-Stokes equations and are typically referred to as direct numerical simulations. The computational requirements of direct numerical simulations have limited most simulations to small-scale 2D models. Recent work by Vervisch et al. (7) presents the simulation of a laboratory-scale "turbulent" premixed V-flame in two dimensions and represents the current state-of-the-art in 2D direct numerical simulations. Tanahashi et al. (8, 9) performed compressible direct numerical simulations of turbulent, premixed hydrogen flames in three dimensions with detailed hydrogen chemistry. Their simulations, performed in a domain <1 cm in each linear dimension, are the first 3D simulations of premixed hydrogen combustion with detailed chemistry.
Low-Mach-number models, which analytically remove acoustic waves on the scale of this problem, provide an alternative formulation for simulation of premixed flames with detailed chemistry. This type of low-Mach-number combustion model was first introduced by Rehm and Baum (10) and rigorously derived from low-Mach-number asymptotic analysis by Majda and Sethian (11). Najm et al. (12-14), and Bell et al. (15) have used low-Mach-number models for simulation of vortex flame interaction with detailed chemistry. Bell et al. (16) used an adaptive low-Mach-number model to simulate a turbulent premixed methane flame in three dimensions with detailed chemistry in an idealized 
(1) cm3 domain.
In this work, we scale up the simulation of Bell et al. (16) to a laboratory-scale turbulent rod-stabilized premixed methane V-flame. This simulation, which models a full laboratory-scale flame by using detailed chemistry and transport, treats a domain more than three orders of magnitude larger in volume than that of any previous efforts and represents a major increment in simulation complexity. We discuss the computational algorithm, techniques used to gather flame statistics, details of the simulation setup, and grid resolution requirements for the adaptive low-Mach-number implementation. We then present detailed comparisons with experimental data. Although there are no explicit subgrid models in the simulation, there are a number of issues that must be addressed in performing this type of computation, such as the fidelity of the chemical kinetics and transport model and characterization of boundary conditions. These issues are discussed in the context of the comparison of experiment and computation.
|
Experimental Configuration and Diagnostics |
|---|
|
TopAbstractExperimental Configuration and...Computational ModelExperimental ComparisonsConclusionReferences |
|---|
= 0.7 exits a 5-cm-diameter circular nozzle with a mean axial velocity of 3 m/s. Turbulence is introduced by a perforated plate mounted 9 cm upstream of the nozzle exit. The integral length scale, 
t, of the turbulence measured by particle image velocimetry (PIV) at the nozzle exit is 3.5 mm. The fluctuation intensity is anisotropic at 7.0% and 5.5% in the axial and radial directions, respectively, relative to the mean axial velocity. The flame is stabilized by a 2-mm-diameter rod spanning the nozzle at its exit. The visible flame extends 15 cm or more downstream from the rod.
|
12 x 12 cm, and the pixel resolution was 0.065 mm per pixel. A cyclone-type particle seeder was used to seed the reactant stream with 0.4-0.6 mm Al2O3 particles (AKP-15, Sumitomo Chemical, Osaka). Data acquisition and analysis were performed on 448 image pairs using software developed by Wernet (17). The two-pass, adaptive cross-correlation interrogation regions of 32 x 32 pixels with final 50% overlapping gave a velocity field spatial resolution of 1 mm. |
Computational Model |
|---|
|
TopAbstractExperimental Configuration and...Computational ModelExperimental ComparisonsConclusionReferences |
|---|
consists of a spatially uniform component, p0, and the perturbation, 
p0 x M2, where M is the local Mach number. The methodology treats the fluid as a mixture of perfect gases and uses a mixture-averaged model for differential species diffusion, ignoring Soret, Dufour, and radiation effects. To second order in the small parameter M, the low-Mach-number equations for flow in an unconfined domain at atmospheric pressure are given by
 | [1] |
 | [2] |
 | [3] |
where 
is the density, U is the velocity, Ym is the mass fraction of species m, h is the mass-weighted enthalpy of the gas mixture, T is the temperature, and 
is the net destruction rate for species m due to chemical reactions. Also, 
is the thermal conductivity, 
is the stress tensor, cp is the specific heat of the mixture, hm(T) and Dm are the enthalpy and species mixture-averaged diffusion coefficients of species m, respectively, and g
is the gravitational acceleration. These evolution equations are supplemented by an equation of state for a perfect gas mixture,
 | [4] |
where Wm is the molecular weight of species m, and R is the universal gas constant.
Note that Systems 1-3 do not admit acoustic waves because the thermodynamic pressure field is essentially constant. By differentiating the equation of state (Eq. 4) in the frame of the fluid, and by using the conservation equations to replace advective derivatives, we obtain an elliptic constraint on the evolving velocity field
 | [5] |
where W is the mean molecular weight.
The chemical kinetics are modeled by using the DRM-19 subset of the GRI-Mech 1.2 methane mechanism (18). DRM-19 is a detailed mechanism containing 20 chemical species and 84 fundamental reactions. Transport and thermodynamic properties are from the GRI-1.2 databases. Our basic discretization algorithm combines a symmetric operator-split coupling of chemistry and diffusion processes with a density-weighted approximate projection method for incorporating the velocity divergence constraint arising from the low-Mach-number formulation. This basic integration scheme is embedded in a parallel adaptive mesh refinement algorithm. Our approach to adaptive refinement is based on a block-structured hierarchical grid system composed of nested rectangular grid patches (19, 20). The adaptive algorithm is second-order accurate in space and time and discretely conserves species mass and enthalpy. See ref. 21 for details of the low-Mach-number model and its numerical implementation and ref. 16 for previous applications of this methodology to the simulation of premixed turbulent flames.
Nozzle Characterization. The experimental characterization of the turbulent fluctuations consists of measures for the integral scale length and intensity and does not include sufficient information to completely specify the inlet flow. To generate realistic inlet turbulence, we performed two auxiliary nonreacting simulations. The first incorporated the curved boundary of the nozzle walls, the inlet jets, and turbulence generation plate by using an embedded boundary algorithm for 3D time-dependent compressible gas dynamics (see ref. 22). The integration procedure in this code is time-explicit and, in the present case, is severely constrained by Courant-Friedrichs-Levy stability limits due to acoustic waves. Despite this limitation, the solution was evolved to a statistically steady state. Statistics then were collected on the time-dependent fluctuation profiles at the nozzle exit. In the boundary layer region 1.25 mm from the nozzle wall, the rms fluctuation level in the axial velocity increased to 100% over its mean value. The radial fluctuations exhibited a corresponding 80% decrease. Away from the nozzle walls, the turbulent intensities of the resulting flow matched the experimental data taken from the central region of the nozzle exit plane, including the observed anisotropy. Using this procedure to generate the time-dependent boundary data for the reacting flow simulation would require running the two codes simultaneously and managing the coupling of the resulting data sets for the duration of the reacting flow simulation.
Alternatively, we may assume that the fluid in the nozzle is an isothermal homogeneous gas and that the turbulence generation in the central region of the nozzle is not strongly influenced by the side walls of the nozzle. Under these conditions, the flow through the nozzle may then be approximated in a Lagrangian frame by using an temporal evolution of an array of perturbed jets. The diameter and spacing of the jets correspond to those of the turbulence generation plate used in the experiment. This second model was carried out over a triply periodic domain, 5 cm on a side, by using the incompressible flow algorithm discussed in ref. 19. The flowfield was initialized with a hexagonal array of 3.2-mm jets spaced 4.8 mm apart, and the initial velocity magnitudes were chosen so that they average over the domain to the experimentally measured mean flow rate. The simulation was performed on a 2563 grid by using a kinematic viscosity of 1.6 x 10-5 m2/s, corresponding to a methane-air mixture at = 0.7 and T = 300 K, and the flow was evolved for 0.03 s (i.e., the mean residence time of the flow in the nozzle). At t = 0.03 s, the integral length scale, turbulence intensity, and anisotropy of the resulting simulated flow agreed with the limited experimental data measured in the central region of the exit plane. This periodic model for the nozzle flow has a distinct practical advantage over the first approach, namely, that we can cycle through the fluctuation data as often as necessary. To properly reflect the presence of the nozzle wall, the periodic data were shaped using radial fits of the axial and radial fluctuations from the flowfield computed by using the embedded boundary algorithm.
Reacting Flow Simulation Setup. The computational domain for the reacting flow simulation is a cube, 12 cm on a side, with the nozzle exit centered on the lower face. The top and sides of this domain are approximated as outflow boundary conditions that express hydrostatic balance with the vertical gravitational force. The reactant stream inflow profile is formed by superimposing a top-hat velocity profile with thin boundary layers at the nozzle wall on the shaped turbulent fluctuations obtained from the nonreacting computations discussed above. The flame stabilization rod is modeled as a 2-mm-wide no-flow strip at the nozzle exit with smooth increase to mean conditions away from the rod. The mean velocity is scaled to match the experimental mean of 3m/s. We specify an air coflow of 1.5 m/s into the bottom of the domain outside the nozzle to control the early time-transient vortex ring that forms between the ambient gas and the outer layer of the reactant flow downstream of the inflow. This coflow provides a relatively crude approximation to the laboratory environment of the flame. The solution in the reacting flow region is initialized with room-temperature stagnant air throughout the domain and a small hot region just above the rod. As the flow evolves, the heated air ignites a flame near the inlet, and the flame surface propagates downstream. The initial evolution was carried out with a two-level adaptive grid hierarchy, where a factor-of-two refinement from the base grid of 963 dynamically tracks regions of high vorticity and chemical activity (the flame front) until the flame is fully established. Then, an additional factor-of-two refinement was placed at the flame surface and high vorticity regions for an effective fine-grid resolution of 
x = 312.5 µm. The refinement strategy resulted in 12% of the 1,728 cm3 domain being refined to the highest level. The resulting flame was evolved until it reached statistical equilibrium.
It is not feasible, with currently available computer hardware, to fully resolve to high accuracy all of the scales of a laboratory-scale flame of the size considered here. The relevant question is whether the resolution we can obtain is sufficient to capture the interaction between the flame and turbulence as manifested in the basic flame morphology and mean velocity fields. The two significant ingredients needed to meet this requirement are that we accurately represent the burning velocity and accurately represent turbulent inflow conditions. For a premixed methane flame at a fuel equivalence ratio of = 0.7, the thermal thickness is 600 µm. High-resolution 1D simulations using GRI-Mech 3.0 (18), the best available chemical mechanism for methane combustion, predict a laminar burning velocity of 19.00 cm/s. In general, a numerical burning velocity is a function of not only spatial and temporal resolution but also depends on the chemical mechanism and transport properties and details of the numerical algorithm. At a spatial resolution of 312.5 µm and a temporal resolution of 19 µs, our methodology predicts a burning velocity based on DRM-19 of 19.04 cm/s, representing an error of <0.25%. We note that, at this resolution, there are errors of as much as 10% in some of the details of the chemical pathways; nevertheless, the burning velocity is accurately predicted.
The other key scale that must be resolved is that of the inflowing turbulence. Standard estimates show that the Kolmogorov scale of the inflowing turbulent reactants is 220 µm. Although the Kolmogorov scale is often quoted as the smallest scale that needs to be resolved, Moin and Mahesh (23) note that this requirement is probably too stringent and suggest that multiples of 4-8 are adequate to obtain good statistical matches using spectral discretizations. Their work suggests that the x = 312.5 µm needed to accurately predict the burning velocity is sufficient to resolve the inflow turbulence. To further test the resolution, we performed a limited evolution at an additional factor-of-two refinement (x = 156.25 µm). Although we were not able to run at this resolution long enough to compute reasonable statistics, comparison of the instantaneous data with the coarser evolution shows no appreciable difference, suggesting that x = 312.5 µm provided adequate resolution of the flow.
Once the flow reached statistical equilibrium at x = 312.5 µm, we evolved the flame for a total of 43.3 ms, and the analysis was based on simulation data at 0.44-ms intervals. Central processing unit requirements depended on the refinement strategy, but the computation progressed at 38 µs per hour on 256 processors of a parallel computer (IBM SP RS/6000 with 375-MHz processors). The total run, including the refinement study, generated 6 terabytes of data for analysis. An image of the instantaneous flame surface at 312.5-µm resolution appears in Fig. 2.
|
Experimental Comparisons |
|---|
|
TopAbstractExperimental Configuration and...Computational ModelExperimental ComparisonsConclusionReferences |
|---|
Scalar Flame Structure. A number of markers for the position of the flame surface have been suggested in the combustion literature, including specific isopleths of fuel or product mass fractions, temperature, or density. On the scale of the 12-cm viewing window for the experimental data and the 1-mm resolution of the PIV imaging technique, all of these measures are identical. A typical centerline slice of the methane concentration obtained from the simulation is shown in Fig. 3 Left. Experimentally, the instantaneous flame location is determined by using the large differences in Mie scattering intensities from the reactants and products to clearly outline the flame (Fig. 3 Right). Comparing the figures, we see that the wrinkling of the flame in the experiment and the computation is of similar size and structure. The different fine-scale structure on the outer edge of the fuel stream may be related to the difference in the dynamics of fuel versus particles or may be a result of under-resolution in the region of the flow that is not refined to the finest level. To characterize the flame brush, which represents the mean flame location, the position of the flame fronts were obtained from 100 PIV images (see Fig. 3) by an edge-finding algorithm for rendering binarized images. Their average produces a map of the mean reaction progress, c
, where c = 0 in reactants and c = 1 in the products. Contours of c for the simulation data were computed by extracting centerline slices of the fluid density and by using the same imaging processing software used to analyze the experimental data to binarize and average the computational data.
|
are depicted in Fig. 4a. Given the uncertainties in specifying the problem (discussed below), the data show excellent agreement in the mean flame structure on the left and reasonably good agreement on the right. The asymmetries in the computational data suggest that the sample size is somewhat too small to obtain fully converged statistics. To characterize the flame-spreading angle more precisely, we averaged the left and right c curves and fit them to a line. The angle between those fits and the vertical is plotted in Fig. 4b as a function of c. For c = 0.5, the simulation predicts a flame angle of 13.4° compared with 12.2° for the experiment, representing a 10% overestimate of the flame angle. For the central range of c, the simulation and experiment data exhibit similar spreading trends. To provide a more detailed comparison of the flame brush growth, in Fig. 4c, we plot the thickness of the flame brush for the experiment and the simulation as a function of height. The data show that computation underpredicts the growth of the flame brush thickness by 25%.
|
|
W
) and transverse (U) velocity components in a vertical plane centered on the circular nozzle for the simulation and the experiment. Below each of the color plates, horizontal cuts of the averaged data show direct comparisons between the simulated and measured data at height intervals of 20 mm. Profiles are shown only where experimentally meaningful values were available; in the outer shear layer, the PIV diagnostic detects only intermittent seeding due to mixing with the unseeded laboratory air. Data from the simulation was averaged over 43.3 ms. The c = 0.1 and c = 0.9 contours (dashed lines) are superimposed on the 2D images to indicate the mean flame zone location. Away from the rod, the computed and experimental velocity data show good qualitative agreement. The simulation captures with remarkable fidelity the major features of the experimental data, i.e., flame-generated outward deflection in the unburned gases, inward flow convergence, and a centerline flow acceleration in the burned gases. A notable difference in the velocity fields is that the experimental measurements just above the nozzle exhibit a depression in axial velocity at the centerline not apparent in the simulation data, which conversely shows a localized acceleration there. In the experiment, the flow near the rod involves complex vortex shedding, leading to a turbulent wake, which is not represented accurately by the simplistic treatment of the rod boundary condition we used in the computation. As a result, the simulated flame is anchored to the top of the rod, whereas in the experiment the flame is stabilized in the shear layer between the recirculation zone behind the rod and the reactant flow (24). The rod remains at a temperature well below that of the combustion products. We speculate that this effect leads to an enhanced centerline acceleration in the simulated results relative to the experiment where the flame structure is delayed somewhat downstream. Nevertheless, the bulk of flame is entirely outside the wake so that its dynamics are determined predominantly by its interaction with the reactant flow turbulence.
|
Discussion. The discrepancies between the computation and the experiment underscore the formulation issues that arise in performing these types of simulations. One issue concerns the characterization of boundary conditions. The importance of properly characterizing the inflow turbulence is obvious. Also, the flame stabilization region has not been simulated here in detail, which appears to affect the growth of the flame brush close to the rod. Finally, our simple approximation of the complex rod stabilization process has generated an anomalous jet in the center of the simulated flow. Although dramatic improvements are likely achievable by extending the geometrical complexity of simulation algorithm to allow modeling of the stabilization process in detail, an alternative approach would involve modifying the experiment directly in an attempt to simplify the flame stabilization mechanism itself. Initial experiments have been performed using a thin heated wire as a flame stabilizer. The results from this configuration will provide a better basis of comparison with simulations by effectively eliminating the recirculation zone and hence the need to model it. Another more subtle issue arises because the experiment we considered was an open flame in a large laboratory. It is neither possible nor desirable to simulate the entire laboratory, and it is impossible to have anything but a simplistic guess of boundary conditions to represent the laboratory's effect. As with flame stabilization, the best approach to address this problem may instead be to modify the experiment so that the flame is more completely isolated from the laboratory. Thus, one element to improving the comparison between experiment and computation is to design the experiment to minimize uncertainties associated with its characterization.
|
|
Conclusion |
|---|
|
TopAbstractExperimental Configuration and...Computational ModelExperimental ComparisonsConclusionReferences |
|---|
The computation documented here represents an advance in the tools available for studying reacting flow. The ability to perform these types of computations has the potential to have a substantive impact on the study of turbulent combustion. In particular, by designing experiments that are well-characterized and specifically designed as companions to such simulations, we can not only provide a much more comprehensive view of a turbulent flame; we can also establish fundamental linkages between turbulent flame experiments and basic combustion chemistry.
The computations were performed at the National Energy Research Scientific Computing Center. This work was supported by the Office of Science through the Scientific Discovery through Advanced Computing program by the Office of Advanced Scientific Computing Research, Mathematical, Information, and Computational Sciences Division, and through the Office of Basic Energy Sciences, Chemical Sciences Division, under U.S. Department of Energy Contract DE-AC03-76SF00098.
|
Footnotes |
|---|
Present address: Department of Mechanical Engineering, University of Ottawa, 161 Louis Pasteur, Room A215, Ottawa, ON, Canada K1N 6N5. 
To whom correspondence should be addressed at: Center for Computational Science and Engineering, Lawrence Berkeley National Laboratory, Mail Stop 50A-1148, 1 Cyclotron Road, Berkeley, CA 94720-8142. E-mail: jbbell{at}lbl.gov.
|
References |
|---|
|
TopAbstractExperimental Configuration and...Computational ModelExperimental ComparisonsConclusionReferences |
|---|
CiteULike 
Complore 
Connotea 
Del.icio.us 
Digg What's this?
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
T|| = 103 K/mm.
t, is computed by fitting to c
