- Email: [email protected]
Large eddy simulation of a premixed flame with approximate deconvolution modeling
Proceedings of the Combustion Institute, Volume 29, 2002/pp. 1995–2000
LARGE EDDY SIMULATION OF A PREMIXED FLAME WITH APPROXIMATE DECONVOLUTION MODELING JOSEPH MATHEW Department of Aerospace Engineering Indian Institute of Science Bangalore 560 012, India
Approximate deconvolution modeling is a very recent approach to large eddy simulation of turbulent flows. It has been applied to compressible flows with success. Here, a premixed flame which forms in the wake of a flameholder has been selected to examine the subgrid-scale modeling of reaction rate by this new method because a previous plane two-dimensional simulation of this wake flame, using a wrinkling function and artificial flame thickening, had revealed discrepancies when compared with experiment. The present simulation is of the temporal evolution of a round wakelike flow at two Reynolds numbers, Re ⳱ 2000 and 10,000, based on wake defect velocity and wake diameter. A Fourier-spectral code has been used. The reaction is single-step and irreversible, and the rate follows an Arrhenius law. The reference simulation at the lower Reynolds number is fully resolved. At Re ⳱ 10,000, subgrid-scale contributions are significant. It was found that subgrid-scale modeling in the present simulation agrees more closely with unresolved subgrid-scale effects observed in experiment. Specifically, the highest contributions appeared in thin folded regions created by vortex convection. The wrinkling function approach had not selected subgrid-scale effects in these regions.
Introduction Large eddy simulation (LES) is well-established for incompressible flows. The dynamic subgrid-scale model [1], which continuously estimates a model parameter during the course of the simulation, has proved to be the adequate development to take LES to analyses of incompressible, engineering flows [2]. This success has encouraged searches for methods to make LES more widely applicable. Turbulent reacting flows present both the need and the opportunity for applications of LES. Such flows exhibit large-scale organization which determines mixing and, in turn, the progress of reactions. The associated large-scale unsteadiness places severe demands on modeling for direct computations of the mean field (Reynolds-averaged Navier-Stokes simulations). LES is a natural choice for such flows because there is no modeling of these large-scale structures. But, LES for reacting flow is still developing because, in addition to special difficulties, the heuristic modeling allowed by incompressible flow cannot be extended to general flows. LES requires models for the effects of unresolved small scales on the large scales arising from every nonlinear term. These are often termed subgridscale (SGS) models denoting the modeling of scales which are smaller than grid spacing. For incompressible flows, the nonlinear terms are the convection terms, and SGS effects appear as a tensor, like the stress tensor. Heuristic SGS modeling, extended
from physical models, takes this tensor to depend on the computed strain rate similar to the models for Reynolds stress. In both cases, the basic physical model is the stress-strain rate relationship for Newtonian fluids. For general flows, like compressible or reacting flows, there are other nonlinear terms. Then, modeling must be done in some other way because physical models to guide heuristic modeling may not be available for every term. Often, several terms are either ignored or assumed to be modeled together. A recent development in LES which circumvents this difficulty is the approximate deconvolution model (ADM) proposed by Adams and Leonard [3]. Several papers reporting excellent results from a priori and a posteriori tests of LES with ADM for non-reacting, compressible flows have appeared [4–6]. The present study was undertaken to examine the effectiveness of ADM for reacting flow. Other simulations of compressible flows have already shown that ADM is able to provide good modeling for convection. So, the specific nonlinearity considered here is the Arrhenius reaction rate term. To add focus, a premixed flame which forms in the wake of a flameholder has been considered because an LES of this configuration, studied by Nottin et al. [7] using other techniques, was reportedly found wanting: there were qualitative differences in the SGS contributions between laboratory measurements and LES computations. They had used a Smagorinsky-type model for the flow and a wrinkling function for the reaction rate term. They concluded
1995
1996
TURBULENT COMBUSTION—Large Eddy and Direct Numerical Simulations
that the sensor for the wrinkling function was inadequate and also that some kind of dynamic procedure was needed. The ADM is potentially better because it combines beneficial effects of the Smagorinsky and scale similarity models and is also a dynamic procedure. LES of realistic premixed flames faces an essential difficulty: The flame thickness is much smaller than the smallest turbulence eddy scale we wish to resolve in an LES. One approach, in different manifestations, has been to track the flame front (G equation) or find the flame surface density [8,9] or the flame wrinkling factor [10] and thence obtain the low-pass filtered reaction rate. Another approach has been to obtain a partial resolution of flame structure on the LES grid by thickening the flame artificially and also use a wrinkling function to compensate for the modified flame-turbulence interactions of the thickened flame. Nottin et al. [7] used the thickened flame method to simulate an experiment in which a propane/air mixture burns in the wake of a flame holder with upstream acoustic excitation. A two-dimensional simulation was considered adequate for the flow near the flame holder. They found global statistical quantities from simulation and experiment to agree very well, but local quantities showed discrepancies. Specifically, the SGS contributions to reaction rate in the simulation were largest in stretched regions around shear layer vortices, whereas the maximum of the unresolved part in experiments occurred in the folds of concentration between a vortex and braid. Although ADM may offer benefits when incorporated into other approaches to reacting flow LES, the present study has the limited aim of assessing its potential for the thickened flame method. A flow which is qualitatively similar to that considered by Nottin et al. [7] has been simulated to determine whether the discrepancy which they observed can be overcome. ADM is described in the following section. Characteristics of the underlying direct numerical simulation (DNS) code have been discussed previously. Accordingly, only the specifics for premixed reacting flow are presented. Next, overall characteristics of the simulations at two Reynolds numbers (fully resolved at 2000 and requiring ADM at 10,000) are presented. A detailed examination then shows that ADM provides SGS contributions to the reaction rate which are closer to experiments. Approximate Deconvolution Modeling The essential aspects of this method can be understood by applying it to the generic, one-dimensional transport equation for u(x,t), u f(u) Ⳮ ⳱0 (1) t x where f(u) is a nonlinear function. Computing a
large wavenumber solution, as in an LES, implies a filtering of equation 1: u¯ f(u) Ⳮ G* ⳱0 t x
(2)
Here, G is a low-pass filter and filtered quantities are obtained by convolution: u¯ ⳱ G * u ⳱ G (x ⳮ x⬘) u(x⬘) dx⬘. The nonlinear term can also be written in terms of the filtered variable. Then equation 2 reads u¯ f(u¯) Ⳮ ⳱R t x
(3)
with a remainder R⳱
f(u¯) f(u) ⳮG* x x
For an LES of incompressible flow, the differential equation 1 is replaced by the Navier-Stokes equations, and the remainder is the SGS stress term, ⳮsij/xj, where sij ⳱ uiuj ⳮ u¯iu¯j and ui are velocity components. Subgrid-scale modeling finds models Rm(u¯) for the remainder as a function of the computed solution u¯. While eddy viscosity models replace this term with a strain rate, the scale similarity model of Bardina [11] uses the computed field to estimate this term. The approximate deconvolution model of Adams and Leonard [3] sets Rm ⳱
f(u¯) f(u*) ⳮG* x x
(4)
where u*(x,t) ⳱ Q * u¯ is an approximation to u(x,t) obtained by approximate deconvolution. On substituting equation 4 into equation 3, and assuming, consistently with equation 4, that G * u* G * u, the LES problem reads G*
u*
冤 t
Ⳮ
f(u*) ⳱0 x
冥
(5)
For a numerical solution, equation 5 can be split into an integration of the original differential equation 1 without any explicit SGS term, followed by a filtering. For example, with superscript n indicating timestepping, we can find u*(n) ⳱ Q * u¯(n) by approximate deconvolution of the LES field u¯(n), advance through an integration step u*(nⳭ1) ⳱ u*(n) Ⳮ Dt
u* Ⳮ O(Dt)2 t
(6)
and filter to obtain the LES field u¯(nⳭ1) ⳱ G * u*(nⳭ1). This three-stage procedure ensures that the SGS modeling prescribed by the ADM (equation 4) has been incorporated. When executed as a sequence, the two filtering steps between each integration step (equation 6) are equivalent to the single
LES WITH APPROXIMATE DECONVOLUTION
filtering Q * G * u*(nⳭ1). If Q were the exact inverse of G, the resultant filtering Q * G * u*(nⳭ1) would not alter the field. But, as an approximate inverse, QG removes content at the high wave number end. This equivalence between filtering and the approximate deconvolution model implies that implementation in DNS codes involves just an explicit filtering step. In addition to the ADM procedure discussed above, Adams and Leonard [3] added an ad hoc regularization term v(G2 ⳮ I) * u¯ to their filtered equations. They took G2 ⳱ QG. It is straightforward to show that the effect of adding this term is equivalent to filtering u¯ with G2 every 1/(vDt) time steps within the temporal truncation error of the integration method. Such an interpretation allows the use of existing DNS codes easily: integration procedures are not modified. Instead, there is an additional, or secondary, filtering of the variables between time-steps. It has also been shown by Mathew [12] that the form of the ad hoc term is close to that which can be obtained by a formal procedure leading to such terms for every nonlinear term. But, using the ad hoc term is simpler and seems to have been adequate.
Reacting Flow Modeling The flow is governed by the Navier-Stokes equations. For premixed flames at unit Lewis number, a progress variable c(x, t) determines species concentrations (mass fractions Yi; i ⳱ F, O, P for fuel, oxidizer, and product) and temperature T from the relations: c ⳱ YP/YP,b ⳱ (T ⳮ Tu)/(Tb ⳮ Tu), YF ⳱ (1 ⳮ YP)/(1 Ⳮ c), and YO ⳱ cYF. Subscripts b and u denote burned and unburned states, respectively. The reaction is single-step and irreversible. For simplicity, the unburned mixture was taken to be at the stoichiometric ratio; c ⳱ 1 and YP,b ⳱ 1. The transport equation for the progress variable is Dc ⳱
1 2c Ⳮ w ˙P ReSc
(7)
where D /t Ⳮ u • and u is the velocity. Two values of the Reynolds number Re were used, but the Schmidt number Sc ⳱ 1 in both simulations. The Arrhenius type reaction rate model considered here is w˙P ⳱ AYFYO exp(ⳮTa/T), and w˙P ⳱ w˙P(c) only. In the simulations, A ⳱ 105, Ta/Tb ⳱ 8, and Tb/Tu ⳱ 2. On filtering equation 7 to perform an LES, remainder terms of the form introduced in equation 3 appear. For the convection and reaction rate terms, they are Rt ⳱ u¯ • c¯ ⳮ u • c and Rw ⳱ w ˙ P(c) ⳮ w˙P(c¯), respectively. These are the subgrid contributions which were obtained from ADM.
1997
Simulations and Results The experiment considered by Nottin et al. [7] is of a premixed flame which is held in the wake of a triangular prism flameholder in a rectangular duct. Their simulation was a plane two-dimensional approximation of the experiment. The flow was acoustically forced using loudspeakers placed upstream. The present simulations are of the temporal evolution of unforced, round, wake-like flows. Flow evolution at successive times in the simulation can be compared to flow development at successive downstream stations of laboratory flows. A notable approximation is that these are incompressible flow simulations, but it should be acceptable because the study is of the modeling of the Arrhenius reaction rate term. Also, compressibility effects were not important for Nottin et al. [7] either in their low Mach number flow. However, significantly, the present three-dimensional simulations preserve essential turbulence dynamics. Simulations were performed using a code which had been developed for DNS. Characteristics of the algorithm and DNS solutions have been reported before [13,14]. It is a Fourier-spectral algorithm in Cartesian coordinates with physical space evaluation of nonlinear terms. The initial wake-like velocity profile was constructed from a thin cylindrical vortex sheet. The streamwise component rises from zero in the cylindrical core to a level across the vortex sheet (tanh velocity profile). The velocity scale was U and the mean diameter D of the initial vortex sheet was the length scale. The Reynolds number Re is based on these two scales. The initial flame lies in the annular shear layer enclosing product and surrounded by premixed reactant. The product distribution (YP ⳱ 1 ⳮ u) mimics the region immediately downstream of flame holders where burned gases accumulate. Two simulations are presented below. The initial fields were the same. Both solutions had periods 3D ⳯ 3D ⳯ 3D on a grid of 64 ⳯ 64 ⳯ 64 points. The first simulation is at Re ⳱ 2000. At this relatively low Reynolds number, the flow is essentially laminar and there is no loss of resolution on this grid. Subgrid-scale modeling is not required. The second is at Re ⳱ 10,000. Now, the initial evolution through vortex ring dynamics is faster and is followed by a breakdown to turbulent flow. ADM provides SGS modeling. The filter G (equation 2) was a Gaussianlike, Pade-type filter specified by a parameter ⳮ0.5 ␣ 0.5. When ␣ is increased, filter width increases and the cutoff of the effective filter QG moves closer to the high wave number end (see Adams and Leonard [3] for details). Here, a high value ␣ ⳱ 0.4 could be used because the underlying scheme is spectral. The initial vortex sheet suffers a Kelvin-Helmholtz-like instability. Vorticity grows and rolls up into a sequence of vortex rings. Since the most amplified
1998
TURBULENT COMBUSTION—Large Eddy and Direct Numerical Simulations
Fig. 2. Distributions of reaction rate w˙P(c) at successive times (t ⳱ 2,3,4,5,7, and 9 from left to right). One-half of streamwise period is shown in each frame. Gray level varies from light to dark as reaction rate increases. (a) Re ⳱ 2000; (b) Re ⳱ 10,000. Fig. 1. Distributions on sections from simulation at Re ⳱ 2000 (t ⳱ 5). Gray level varies from light to dark as the value of the quantity increases. (a) Vorticity magnitude on longitudinal section | ⳯ u|; (b) product mass fraction YP ⳱ c.
wavelength scales with sheet thickness, small perturbations were added to select this rollup and a subsequent pairing [13]. As an example of the solution, Fig. 1a shows vorticity distribution on a longitudinal section at t ⳱ 5 (Re ⳱ 2000). A sequence of vortex rings has formed. Product distribution is shown in Fig. 1b. The action of the vortex rings bringing in surrounding fluid can be seen as tongues of lower concentration along braids between neighboring rings. Reaction rates (part of Fig. 2a, discussed below) are highest around ring vortices and in the braids where unburned fluid has reached and temperatures are high. In temporal simulations, flow structures which form at successive times correspond to those which develop at successive downstream stations in laboratory wakes. For a comprehensive picture, reaction rate fields at successive times (t ⳱ 2,3,4,5,7,9) have
been arranged in a sequence in Fig. 2. One-half of the streamwise extent at each time has been shown. Fig. 2a is for Re ⳱ 2000 and Fig. 2b is for Re ⳱ 10,000. Although the initial evolution looks similar, soon differences appear. At the lower Reynolds number the reaction zone eventually becomes a nearly uniform annular sheet, insensitive to the vortices which dissipate. At the higher Reynolds number, the reaction zone fills the wake due to strong (turbulent) mixing. The total mass of product scaled with the initial mass is plotted as a function of time in Fig. 3. The greater mixing in a turbulent flow can be seen in the significantly higher product formation rate following breakdown of the flow beyond t 7. The discrepancies in SGS contributions observed by of Nottin et al. [7] were close to neighboring vortices. Fig. 2 shows that the data at t ⳱ 4 (Re ⳱ 10,000), when a sequence of large vortices with welldefined braids exist, can be used to examine this issue in the present simulation. Fig. 4 shows distri˙ P(c*), butions of product, filtered reaction rate w ˙ P(c*) ⳮ w˙P(c¯). vorticity, and the SGS reaction rate w Four parts of a vortex have been marked in Fig. 4d. We may term A the outer rear, B the outer head, C the inner fold, and D the braid. In experiments, Nottin et al. [7] found unresolved SGS values to be
LES WITH APPROXIMATE DECONVOLUTION
Fig. 3. Growth of mass of product scaled with initial mass. (—) Re ⳱ 2000; (—) Re ⳱ 10,000.
1999
in their experiments and here. But, in the experiments there was little contribution in the braids (D) unlike here. In their simulations [7], subgrid values were highest in region B, unlike both their experiment and the present simulation. High levels were seen in regions A and D but not in C. So we find that ADM is closer to experiment with regard to the folded regions, but differs in the contributions at the braids. Also, unlike the wrinkling function which determined from the vorticity field, here the SGS modeling in ADM depends directly on the term it originates from: The vorticity has peaks within the ring vortex cores in Fig. 4c but the SGS effect on reaction rate vanishes in the cores. It must be noted that there were other differences between simulation and experiment: In experiments, reaction zones had become thicker and peak levels had diminished with downstream distance. In their simulations, reaction rate intensity persisted downstream, and the flame remained thinner. The absence of SGS effects in the braids in experiment could be due to this thickening. Then, the reaction rate can be resolved and SGS effects disappear. This argument is further supported by the presence of SGS effects in the thin, first braid (next to flame holder; see Fig. 4d in Ref. [7]). So SGS effects seen in the braids in the present ADM simulation and theirs need not be a discrepancy in modeling; rather, the cause is a difference in underlying flow structure. Taken together, there is strong support that ADM selects the right SGS contributions to reaction rate. Conclusions
Fig. 4. Distributions on sections from simulation at Re ⳱ 10,000 (t ⳱ 4). (a) Product mass fraction YP ⳱ c; (b) filtered reaction rate w ˙ P(c*). (c) vorticity magnitude on longitudinal section | ⳯ u|; (d) subgrid-scale contribution ˙ P(c*) ⳮ w˙P(c¯). Quantities in (a)–(c) are to reaction rate w positive, while SGS term in (d) takes on positive and negative values. Gray level is darker for larger magnitudes.
An LES of a premixed flame forming in the shear layer of a round wake has been performed using the approximate deconvolution model. Other simulations of compressible flows have already shown that ADM is able to provide good modeling for convection. A close examination of the SGS effects provided by ADM for an Arrhenius reaction rate term shows that it is close to that observed as unresolved small-scale effects in experiments of a similar flow. This is in contrast to using the wrinkling function with a thickened flame. It remains to be determined whether artificial thickening is the best-suited for LES of reacting flow but present results suggest that ADM may be the better complement to the thickened flame method. Acknowledgments
strongest in region C where the flame is folded close together and has high curvature. Subgrid contributions from the present ADM simulations are also highest in region C. The reaction rate is also highest here. Fig. 4b shows reaction rate reach peak levels where fresh reactants have been carried into the thin folds between vortex and braid (see also Fig. 4a). Subgrid values are also significant in region A, both
It is a pleasure to thank Prof. R. Friedrich of TU-Munich who introduced me to the ADM for LES. I also thank Dr. Amit J. Basu who developed the underlying DNS code. REFERENCES 1. Germano, M., Piomelli, U., Moin, P., and Cabot, W. H., Phys. Fluids A 3(7):1760–1771 (1991).
2000
TURBULENT COMBUSTION—Large Eddy and Direct Numerical Simulations
2. Moin, P., in Turbulence and Shear Flow Phenomena-2 (E. Lindberg et al., ed.), Royal Institute of Technology, Stockholm, 2001, pp. 1–10. 3. Adams, N. A., and Leonard, A., in Direct and LargeEddy Simulation III (P. Voke et al., ed.), Kluwer Academic Publishers, Dordrecht, 1999, pp. 201–212. 4. Stolz, S., and Adams, N. A., Phys. Fluids 11:1699–1701 (1999). 5. Stolz, S., Adams, N. A., and Kleiser, L., Phys. Fluids 13:997–1015 (2001). 6. Stolz, S., Adams, N. A., and Kleiser, L., Phys. Fluids 13:2985 (2001). 7. Nottin, C., Knikker, R., Boger, M., and Veynante, D., Proc. Combust. Inst. 28:67–73 (2000). 8. Boger, M., Veynante, D., Boughanem, H., and Trouve´, A., Proc. Combust. Inst. 27:917–925 (1998).
9. Hawkes, E. R., and Cant, R. S., Proc. Combust. Inst. 28:51–58 (2000). 10. Weller, H. G., Tabor, G., Gosman, A. D., and Fureby, C., Proc. Combust. Inst. 27:899–907 (1998). 11. Bardina, J., ‘‘Improved Turbulence Models Based on Large Eddy Simulations of Homogeneous, Incompressible, Turbulent Flows,’’ Ph.D. thesis, Department of Mechanical Engineering, Stanford University, Stanford, CA, 1983. 12. Mathew, J., The Deconvolution Method for Large Eddy Simulation, report 2000-FM-7, Department of Aerospace Engineering, Indian Institute of Science, Bangalore, 2000. 13. Mathew, J., and Basu, A. J., Flow Turbulence Combust. 64:71–93 (2000). 14. Mathew, J., Proc. Combust. Inst. 27:1207–1212 (1998).
COMMENTS Stephen B. Pope, Cornell University, USA. There is a substantial difference between the subgrid modeling of the stress and the reaction rate. For the stress, the dominant contribution arises from motions that are just smaller than the smallest resolved scales, whereas for the reaction rate, the contributions may arise from much smaller scales. Can the ADM procedure accurately represent the subgrid reaction rate in such circumstances? Author’s Reply. ADM is a mathematically consistent treatment of small scales close to the smallest represented scales in a calculation. For ADM to work satisfactorily, unrepresentative small scales must have a small effect on the solution. Smaller scales must not be an essential part of the problem. These are general expectations in LES. Therefore, ADM can be useful for the reaction rate for thick, or thickened, flames. When flame structure is thin, other suitable modeling will be required for the reaction rate term even when ADM is adequate for other nonlinearities. This is analogous to wall layer treatments in LES that either resolve wall layer structure with a fine grid or use an additional model. ● Sebastie´n Candel, Ecole Centrale Paris, France. You have nicely illustrated the ADM method. However, the calculations are rather far from conditions prevailing in the experiment. In the case of a ducted flame the fresh gases are initially moving faster than the burned reactants but further downstream the situation is reversed because the
hot products accelerate. This strongly influences the shear layer development. Author’s Reply. A model flow whose flow structure resembles the near wake flow in the experiment was chosen. The temporal, constant density simulation was not expected to closely model a ducted flame. It was pointed out that in the present ADM simulation, SGS contributions to the reaction rate appeared in the folded regions as well. A wrinkling function may still be needed as an adjunct to ADM for quantitatively close results, and such a quantitative assessment must use a close model of the flame. ● Anaias Tomboulides, ETH Zurich, Switzerland. In the ADM procedure, there is a relaxation term accounting for energy transfer from large to small scales. It would be interesting to learn how you dealt with that term, particularly in the species concentration equation that includes the reaction rate. Author’s Reply. The ad hoc term which is used with ADM and its interpretation as a secondary filtering has been discussed in the paper. This is higher-order modeling, and its effectiveness has not been established yet. Primary and secondary filtering result in a single equivalent filter, and the solution depends on this equivalent filter only and not separately on the two constituents. The present simulations were conducted with a single equivalent filter, with a very high cutoff wave number, for all variables including concentration.
LARGE EDDY SIMULATION OF A PREMIXED FLAME WITH APPROXIMATE DECONVOLUTION MODELING JOSEPH MATHEW Department of Aerospace Engineering Indian Institute of Science Bangalore 560 012, India
Approximate deconvolution modeling is a very recent approach to large eddy simulation of turbulent flows. It has been applied to compressible flows with success. Here, a premixed flame which forms in the wake of a flameholder has been selected to examine the subgrid-scale modeling of reaction rate by this new method because a previous plane two-dimensional simulation of this wake flame, using a wrinkling function and artificial flame thickening, had revealed discrepancies when compared with experiment. The present simulation is of the temporal evolution of a round wakelike flow at two Reynolds numbers, Re ⳱ 2000 and 10,000, based on wake defect velocity and wake diameter. A Fourier-spectral code has been used. The reaction is single-step and irreversible, and the rate follows an Arrhenius law. The reference simulation at the lower Reynolds number is fully resolved. At Re ⳱ 10,000, subgrid-scale contributions are significant. It was found that subgrid-scale modeling in the present simulation agrees more closely with unresolved subgrid-scale effects observed in experiment. Specifically, the highest contributions appeared in thin folded regions created by vortex convection. The wrinkling function approach had not selected subgrid-scale effects in these regions.
Introduction Large eddy simulation (LES) is well-established for incompressible flows. The dynamic subgrid-scale model [1], which continuously estimates a model parameter during the course of the simulation, has proved to be the adequate development to take LES to analyses of incompressible, engineering flows [2]. This success has encouraged searches for methods to make LES more widely applicable. Turbulent reacting flows present both the need and the opportunity for applications of LES. Such flows exhibit large-scale organization which determines mixing and, in turn, the progress of reactions. The associated large-scale unsteadiness places severe demands on modeling for direct computations of the mean field (Reynolds-averaged Navier-Stokes simulations). LES is a natural choice for such flows because there is no modeling of these large-scale structures. But, LES for reacting flow is still developing because, in addition to special difficulties, the heuristic modeling allowed by incompressible flow cannot be extended to general flows. LES requires models for the effects of unresolved small scales on the large scales arising from every nonlinear term. These are often termed subgridscale (SGS) models denoting the modeling of scales which are smaller than grid spacing. For incompressible flows, the nonlinear terms are the convection terms, and SGS effects appear as a tensor, like the stress tensor. Heuristic SGS modeling, extended
from physical models, takes this tensor to depend on the computed strain rate similar to the models for Reynolds stress. In both cases, the basic physical model is the stress-strain rate relationship for Newtonian fluids. For general flows, like compressible or reacting flows, there are other nonlinear terms. Then, modeling must be done in some other way because physical models to guide heuristic modeling may not be available for every term. Often, several terms are either ignored or assumed to be modeled together. A recent development in LES which circumvents this difficulty is the approximate deconvolution model (ADM) proposed by Adams and Leonard [3]. Several papers reporting excellent results from a priori and a posteriori tests of LES with ADM for non-reacting, compressible flows have appeared [4–6]. The present study was undertaken to examine the effectiveness of ADM for reacting flow. Other simulations of compressible flows have already shown that ADM is able to provide good modeling for convection. So, the specific nonlinearity considered here is the Arrhenius reaction rate term. To add focus, a premixed flame which forms in the wake of a flameholder has been considered because an LES of this configuration, studied by Nottin et al. [7] using other techniques, was reportedly found wanting: there were qualitative differences in the SGS contributions between laboratory measurements and LES computations. They had used a Smagorinsky-type model for the flow and a wrinkling function for the reaction rate term. They concluded
1995
1996
TURBULENT COMBUSTION—Large Eddy and Direct Numerical Simulations
that the sensor for the wrinkling function was inadequate and also that some kind of dynamic procedure was needed. The ADM is potentially better because it combines beneficial effects of the Smagorinsky and scale similarity models and is also a dynamic procedure. LES of realistic premixed flames faces an essential difficulty: The flame thickness is much smaller than the smallest turbulence eddy scale we wish to resolve in an LES. One approach, in different manifestations, has been to track the flame front (G equation) or find the flame surface density [8,9] or the flame wrinkling factor [10] and thence obtain the low-pass filtered reaction rate. Another approach has been to obtain a partial resolution of flame structure on the LES grid by thickening the flame artificially and also use a wrinkling function to compensate for the modified flame-turbulence interactions of the thickened flame. Nottin et al. [7] used the thickened flame method to simulate an experiment in which a propane/air mixture burns in the wake of a flame holder with upstream acoustic excitation. A two-dimensional simulation was considered adequate for the flow near the flame holder. They found global statistical quantities from simulation and experiment to agree very well, but local quantities showed discrepancies. Specifically, the SGS contributions to reaction rate in the simulation were largest in stretched regions around shear layer vortices, whereas the maximum of the unresolved part in experiments occurred in the folds of concentration between a vortex and braid. Although ADM may offer benefits when incorporated into other approaches to reacting flow LES, the present study has the limited aim of assessing its potential for the thickened flame method. A flow which is qualitatively similar to that considered by Nottin et al. [7] has been simulated to determine whether the discrepancy which they observed can be overcome. ADM is described in the following section. Characteristics of the underlying direct numerical simulation (DNS) code have been discussed previously. Accordingly, only the specifics for premixed reacting flow are presented. Next, overall characteristics of the simulations at two Reynolds numbers (fully resolved at 2000 and requiring ADM at 10,000) are presented. A detailed examination then shows that ADM provides SGS contributions to the reaction rate which are closer to experiments. Approximate Deconvolution Modeling The essential aspects of this method can be understood by applying it to the generic, one-dimensional transport equation for u(x,t), u f(u) Ⳮ ⳱0 (1) t x where f(u) is a nonlinear function. Computing a
large wavenumber solution, as in an LES, implies a filtering of equation 1: u¯ f(u) Ⳮ G* ⳱0 t x
(2)
Here, G is a low-pass filter and filtered quantities are obtained by convolution: u¯ ⳱ G * u ⳱ G (x ⳮ x⬘) u(x⬘) dx⬘. The nonlinear term can also be written in terms of the filtered variable. Then equation 2 reads u¯ f(u¯) Ⳮ ⳱R t x
(3)
with a remainder R⳱
f(u¯) f(u) ⳮG* x x
For an LES of incompressible flow, the differential equation 1 is replaced by the Navier-Stokes equations, and the remainder is the SGS stress term, ⳮsij/xj, where sij ⳱ uiuj ⳮ u¯iu¯j and ui are velocity components. Subgrid-scale modeling finds models Rm(u¯) for the remainder as a function of the computed solution u¯. While eddy viscosity models replace this term with a strain rate, the scale similarity model of Bardina [11] uses the computed field to estimate this term. The approximate deconvolution model of Adams and Leonard [3] sets Rm ⳱
f(u¯) f(u*) ⳮG* x x
(4)
where u*(x,t) ⳱ Q * u¯ is an approximation to u(x,t) obtained by approximate deconvolution. On substituting equation 4 into equation 3, and assuming, consistently with equation 4, that G * u* G * u, the LES problem reads G*
u*
冤 t
Ⳮ
f(u*) ⳱0 x
冥
(5)
For a numerical solution, equation 5 can be split into an integration of the original differential equation 1 without any explicit SGS term, followed by a filtering. For example, with superscript n indicating timestepping, we can find u*(n) ⳱ Q * u¯(n) by approximate deconvolution of the LES field u¯(n), advance through an integration step u*(nⳭ1) ⳱ u*(n) Ⳮ Dt
u* Ⳮ O(Dt)2 t
(6)
and filter to obtain the LES field u¯(nⳭ1) ⳱ G * u*(nⳭ1). This three-stage procedure ensures that the SGS modeling prescribed by the ADM (equation 4) has been incorporated. When executed as a sequence, the two filtering steps between each integration step (equation 6) are equivalent to the single
LES WITH APPROXIMATE DECONVOLUTION
filtering Q * G * u*(nⳭ1). If Q were the exact inverse of G, the resultant filtering Q * G * u*(nⳭ1) would not alter the field. But, as an approximate inverse, QG removes content at the high wave number end. This equivalence between filtering and the approximate deconvolution model implies that implementation in DNS codes involves just an explicit filtering step. In addition to the ADM procedure discussed above, Adams and Leonard [3] added an ad hoc regularization term v(G2 ⳮ I) * u¯ to their filtered equations. They took G2 ⳱ QG. It is straightforward to show that the effect of adding this term is equivalent to filtering u¯ with G2 every 1/(vDt) time steps within the temporal truncation error of the integration method. Such an interpretation allows the use of existing DNS codes easily: integration procedures are not modified. Instead, there is an additional, or secondary, filtering of the variables between time-steps. It has also been shown by Mathew [12] that the form of the ad hoc term is close to that which can be obtained by a formal procedure leading to such terms for every nonlinear term. But, using the ad hoc term is simpler and seems to have been adequate.
Reacting Flow Modeling The flow is governed by the Navier-Stokes equations. For premixed flames at unit Lewis number, a progress variable c(x, t) determines species concentrations (mass fractions Yi; i ⳱ F, O, P for fuel, oxidizer, and product) and temperature T from the relations: c ⳱ YP/YP,b ⳱ (T ⳮ Tu)/(Tb ⳮ Tu), YF ⳱ (1 ⳮ YP)/(1 Ⳮ c), and YO ⳱ cYF. Subscripts b and u denote burned and unburned states, respectively. The reaction is single-step and irreversible. For simplicity, the unburned mixture was taken to be at the stoichiometric ratio; c ⳱ 1 and YP,b ⳱ 1. The transport equation for the progress variable is Dc ⳱
1 2c Ⳮ w ˙P ReSc
(7)
where D /t Ⳮ u • and u is the velocity. Two values of the Reynolds number Re were used, but the Schmidt number Sc ⳱ 1 in both simulations. The Arrhenius type reaction rate model considered here is w˙P ⳱ AYFYO exp(ⳮTa/T), and w˙P ⳱ w˙P(c) only. In the simulations, A ⳱ 105, Ta/Tb ⳱ 8, and Tb/Tu ⳱ 2. On filtering equation 7 to perform an LES, remainder terms of the form introduced in equation 3 appear. For the convection and reaction rate terms, they are Rt ⳱ u¯ • c¯ ⳮ u • c and Rw ⳱ w ˙ P(c) ⳮ w˙P(c¯), respectively. These are the subgrid contributions which were obtained from ADM.
1997
Simulations and Results The experiment considered by Nottin et al. [7] is of a premixed flame which is held in the wake of a triangular prism flameholder in a rectangular duct. Their simulation was a plane two-dimensional approximation of the experiment. The flow was acoustically forced using loudspeakers placed upstream. The present simulations are of the temporal evolution of unforced, round, wake-like flows. Flow evolution at successive times in the simulation can be compared to flow development at successive downstream stations of laboratory flows. A notable approximation is that these are incompressible flow simulations, but it should be acceptable because the study is of the modeling of the Arrhenius reaction rate term. Also, compressibility effects were not important for Nottin et al. [7] either in their low Mach number flow. However, significantly, the present three-dimensional simulations preserve essential turbulence dynamics. Simulations were performed using a code which had been developed for DNS. Characteristics of the algorithm and DNS solutions have been reported before [13,14]. It is a Fourier-spectral algorithm in Cartesian coordinates with physical space evaluation of nonlinear terms. The initial wake-like velocity profile was constructed from a thin cylindrical vortex sheet. The streamwise component rises from zero in the cylindrical core to a level across the vortex sheet (tanh velocity profile). The velocity scale was U and the mean diameter D of the initial vortex sheet was the length scale. The Reynolds number Re is based on these two scales. The initial flame lies in the annular shear layer enclosing product and surrounded by premixed reactant. The product distribution (YP ⳱ 1 ⳮ u) mimics the region immediately downstream of flame holders where burned gases accumulate. Two simulations are presented below. The initial fields were the same. Both solutions had periods 3D ⳯ 3D ⳯ 3D on a grid of 64 ⳯ 64 ⳯ 64 points. The first simulation is at Re ⳱ 2000. At this relatively low Reynolds number, the flow is essentially laminar and there is no loss of resolution on this grid. Subgrid-scale modeling is not required. The second is at Re ⳱ 10,000. Now, the initial evolution through vortex ring dynamics is faster and is followed by a breakdown to turbulent flow. ADM provides SGS modeling. The filter G (equation 2) was a Gaussianlike, Pade-type filter specified by a parameter ⳮ0.5 ␣ 0.5. When ␣ is increased, filter width increases and the cutoff of the effective filter QG moves closer to the high wave number end (see Adams and Leonard [3] for details). Here, a high value ␣ ⳱ 0.4 could be used because the underlying scheme is spectral. The initial vortex sheet suffers a Kelvin-Helmholtz-like instability. Vorticity grows and rolls up into a sequence of vortex rings. Since the most amplified
1998
TURBULENT COMBUSTION—Large Eddy and Direct Numerical Simulations
Fig. 2. Distributions of reaction rate w˙P(c) at successive times (t ⳱ 2,3,4,5,7, and 9 from left to right). One-half of streamwise period is shown in each frame. Gray level varies from light to dark as reaction rate increases. (a) Re ⳱ 2000; (b) Re ⳱ 10,000. Fig. 1. Distributions on sections from simulation at Re ⳱ 2000 (t ⳱ 5). Gray level varies from light to dark as the value of the quantity increases. (a) Vorticity magnitude on longitudinal section | ⳯ u|; (b) product mass fraction YP ⳱ c.
wavelength scales with sheet thickness, small perturbations were added to select this rollup and a subsequent pairing [13]. As an example of the solution, Fig. 1a shows vorticity distribution on a longitudinal section at t ⳱ 5 (Re ⳱ 2000). A sequence of vortex rings has formed. Product distribution is shown in Fig. 1b. The action of the vortex rings bringing in surrounding fluid can be seen as tongues of lower concentration along braids between neighboring rings. Reaction rates (part of Fig. 2a, discussed below) are highest around ring vortices and in the braids where unburned fluid has reached and temperatures are high. In temporal simulations, flow structures which form at successive times correspond to those which develop at successive downstream stations in laboratory wakes. For a comprehensive picture, reaction rate fields at successive times (t ⳱ 2,3,4,5,7,9) have
been arranged in a sequence in Fig. 2. One-half of the streamwise extent at each time has been shown. Fig. 2a is for Re ⳱ 2000 and Fig. 2b is for Re ⳱ 10,000. Although the initial evolution looks similar, soon differences appear. At the lower Reynolds number the reaction zone eventually becomes a nearly uniform annular sheet, insensitive to the vortices which dissipate. At the higher Reynolds number, the reaction zone fills the wake due to strong (turbulent) mixing. The total mass of product scaled with the initial mass is plotted as a function of time in Fig. 3. The greater mixing in a turbulent flow can be seen in the significantly higher product formation rate following breakdown of the flow beyond t 7. The discrepancies in SGS contributions observed by of Nottin et al. [7] were close to neighboring vortices. Fig. 2 shows that the data at t ⳱ 4 (Re ⳱ 10,000), when a sequence of large vortices with welldefined braids exist, can be used to examine this issue in the present simulation. Fig. 4 shows distri˙ P(c*), butions of product, filtered reaction rate w ˙ P(c*) ⳮ w˙P(c¯). vorticity, and the SGS reaction rate w Four parts of a vortex have been marked in Fig. 4d. We may term A the outer rear, B the outer head, C the inner fold, and D the braid. In experiments, Nottin et al. [7] found unresolved SGS values to be
LES WITH APPROXIMATE DECONVOLUTION
Fig. 3. Growth of mass of product scaled with initial mass. (—) Re ⳱ 2000; (—) Re ⳱ 10,000.
1999
in their experiments and here. But, in the experiments there was little contribution in the braids (D) unlike here. In their simulations [7], subgrid values were highest in region B, unlike both their experiment and the present simulation. High levels were seen in regions A and D but not in C. So we find that ADM is closer to experiment with regard to the folded regions, but differs in the contributions at the braids. Also, unlike the wrinkling function which determined from the vorticity field, here the SGS modeling in ADM depends directly on the term it originates from: The vorticity has peaks within the ring vortex cores in Fig. 4c but the SGS effect on reaction rate vanishes in the cores. It must be noted that there were other differences between simulation and experiment: In experiments, reaction zones had become thicker and peak levels had diminished with downstream distance. In their simulations, reaction rate intensity persisted downstream, and the flame remained thinner. The absence of SGS effects in the braids in experiment could be due to this thickening. Then, the reaction rate can be resolved and SGS effects disappear. This argument is further supported by the presence of SGS effects in the thin, first braid (next to flame holder; see Fig. 4d in Ref. [7]). So SGS effects seen in the braids in the present ADM simulation and theirs need not be a discrepancy in modeling; rather, the cause is a difference in underlying flow structure. Taken together, there is strong support that ADM selects the right SGS contributions to reaction rate. Conclusions
Fig. 4. Distributions on sections from simulation at Re ⳱ 10,000 (t ⳱ 4). (a) Product mass fraction YP ⳱ c; (b) filtered reaction rate w ˙ P(c*). (c) vorticity magnitude on longitudinal section | ⳯ u|; (d) subgrid-scale contribution ˙ P(c*) ⳮ w˙P(c¯). Quantities in (a)–(c) are to reaction rate w positive, while SGS term in (d) takes on positive and negative values. Gray level is darker for larger magnitudes.
An LES of a premixed flame forming in the shear layer of a round wake has been performed using the approximate deconvolution model. Other simulations of compressible flows have already shown that ADM is able to provide good modeling for convection. A close examination of the SGS effects provided by ADM for an Arrhenius reaction rate term shows that it is close to that observed as unresolved small-scale effects in experiments of a similar flow. This is in contrast to using the wrinkling function with a thickened flame. It remains to be determined whether artificial thickening is the best-suited for LES of reacting flow but present results suggest that ADM may be the better complement to the thickened flame method. Acknowledgments
strongest in region C where the flame is folded close together and has high curvature. Subgrid contributions from the present ADM simulations are also highest in region C. The reaction rate is also highest here. Fig. 4b shows reaction rate reach peak levels where fresh reactants have been carried into the thin folds between vortex and braid (see also Fig. 4a). Subgrid values are also significant in region A, both
It is a pleasure to thank Prof. R. Friedrich of TU-Munich who introduced me to the ADM for LES. I also thank Dr. Amit J. Basu who developed the underlying DNS code. REFERENCES 1. Germano, M., Piomelli, U., Moin, P., and Cabot, W. H., Phys. Fluids A 3(7):1760–1771 (1991).
2000
TURBULENT COMBUSTION—Large Eddy and Direct Numerical Simulations
2. Moin, P., in Turbulence and Shear Flow Phenomena-2 (E. Lindberg et al., ed.), Royal Institute of Technology, Stockholm, 2001, pp. 1–10. 3. Adams, N. A., and Leonard, A., in Direct and LargeEddy Simulation III (P. Voke et al., ed.), Kluwer Academic Publishers, Dordrecht, 1999, pp. 201–212. 4. Stolz, S., and Adams, N. A., Phys. Fluids 11:1699–1701 (1999). 5. Stolz, S., Adams, N. A., and Kleiser, L., Phys. Fluids 13:997–1015 (2001). 6. Stolz, S., Adams, N. A., and Kleiser, L., Phys. Fluids 13:2985 (2001). 7. Nottin, C., Knikker, R., Boger, M., and Veynante, D., Proc. Combust. Inst. 28:67–73 (2000). 8. Boger, M., Veynante, D., Boughanem, H., and Trouve´, A., Proc. Combust. Inst. 27:917–925 (1998).
9. Hawkes, E. R., and Cant, R. S., Proc. Combust. Inst. 28:51–58 (2000). 10. Weller, H. G., Tabor, G., Gosman, A. D., and Fureby, C., Proc. Combust. Inst. 27:899–907 (1998). 11. Bardina, J., ‘‘Improved Turbulence Models Based on Large Eddy Simulations of Homogeneous, Incompressible, Turbulent Flows,’’ Ph.D. thesis, Department of Mechanical Engineering, Stanford University, Stanford, CA, 1983. 12. Mathew, J., The Deconvolution Method for Large Eddy Simulation, report 2000-FM-7, Department of Aerospace Engineering, Indian Institute of Science, Bangalore, 2000. 13. Mathew, J., and Basu, A. J., Flow Turbulence Combust. 64:71–93 (2000). 14. Mathew, J., Proc. Combust. Inst. 27:1207–1212 (1998).
COMMENTS Stephen B. Pope, Cornell University, USA. There is a substantial difference between the subgrid modeling of the stress and the reaction rate. For the stress, the dominant contribution arises from motions that are just smaller than the smallest resolved scales, whereas for the reaction rate, the contributions may arise from much smaller scales. Can the ADM procedure accurately represent the subgrid reaction rate in such circumstances? Author’s Reply. ADM is a mathematically consistent treatment of small scales close to the smallest represented scales in a calculation. For ADM to work satisfactorily, unrepresentative small scales must have a small effect on the solution. Smaller scales must not be an essential part of the problem. These are general expectations in LES. Therefore, ADM can be useful for the reaction rate for thick, or thickened, flames. When flame structure is thin, other suitable modeling will be required for the reaction rate term even when ADM is adequate for other nonlinearities. This is analogous to wall layer treatments in LES that either resolve wall layer structure with a fine grid or use an additional model. ● Sebastie´n Candel, Ecole Centrale Paris, France. You have nicely illustrated the ADM method. However, the calculations are rather far from conditions prevailing in the experiment. In the case of a ducted flame the fresh gases are initially moving faster than the burned reactants but further downstream the situation is reversed because the
hot products accelerate. This strongly influences the shear layer development. Author’s Reply. A model flow whose flow structure resembles the near wake flow in the experiment was chosen. The temporal, constant density simulation was not expected to closely model a ducted flame. It was pointed out that in the present ADM simulation, SGS contributions to the reaction rate appeared in the folded regions as well. A wrinkling function may still be needed as an adjunct to ADM for quantitatively close results, and such a quantitative assessment must use a close model of the flame. ● Anaias Tomboulides, ETH Zurich, Switzerland. In the ADM procedure, there is a relaxation term accounting for energy transfer from large to small scales. It would be interesting to learn how you dealt with that term, particularly in the species concentration equation that includes the reaction rate. Author’s Reply. The ad hoc term which is used with ADM and its interpretation as a secondary filtering has been discussed in the paper. This is higher-order modeling, and its effectiveness has not been established yet. Primary and secondary filtering result in a single equivalent filter, and the solution depends on this equivalent filter only and not separately on the two constituents. The present simulations were conducted with a single equivalent filter, with a very high cutoff wave number, for all variables including concentration.