LARGE EDDY SIMULATION OF BLUFF-BODY STABILIZED
FLAME BY USING FLAME SURFACE DENSITY
APPROACH
Rajani Kumar Akula*, Amsini Sadiki† and Johannes Janicka
*Darmstadt University of Technology, Faculty EKT, Petersonstr. 30,
64287 Darmstadt, Germany e-mail: akula@ekt.tu-darmstadt.de
†
Corresponding author: e-mail: sadiki@ekt.tu-darmstadt.de
Key words: Large eddy simulation (LES), premixed combustion, flame surface density
Abstract. Despite significant advances in understanding and modelling of combustion-LES, LES of premixed combustion remains a challenging issue. This problem is more challenging when an incompressible formulation is used where the density coupling is weak, and an improper density coupling can lead to unrealistic flame speeds and thickness. For this purpose the flame surface density approach is used along with an appropriate numerical technique. By incorporating the proper flame wrinkling model, the LES of Bluff-body stabilized premixed flame are carried out and numerical results are compared with experimental data. A good agreement testifies the accuracy of the numerical method.
1 INTRODUCTION
flame wrinkling models available in the literature. Colin’s9 model is a simple one but requires input of turbulent Reynolds number. Charlette et al.10, 11 proposed two versions of power-law flame wrinkling model. In a first version they took power-law exponent as a constant value10 and in the second one they calculated it by using dynamic formulation11. Nevertheless one can use simple turbulent flame speed models like the flame speed model proposed by Pocheau 12 to account for the wrinkling effects. All these models require accurate description of subgrid scale velocity fluctuation u∆′ , which in turn demands proper subgrid scale modelling. In this work, the model by Boger et al. is used in the frame of the flame surface density13.
For the simulations we follow a dynamic SGS model formulation. Germano et al.14 developed a dynamic Smagorinsky SGS model (DSM) which overcomes some of the drawbacks of the Smagorinsky model. This model has many desirable features in that it requires only one input parameter, i.e., ratio of test filter width to grid filter width. DSM has been successfully applied to LES of transitional and turbulent channel flows14 and both incompressible and compressible isotropic turbulence15. However, the dynamically computed model coefficient was averaged either in the global volume of the domain or in a homogeneous plane. To overcome this problem Meneveau et al.16 proposed a Lagrangian dynamic model that we used in this present work. It is worth mentioning that for better predictions the use of one equation dynamic models based on subgrid scale kinetic energy17 is encouraged.
With regard to numerical implementation, an important problem is how to include the variable density in an incompressible flow solver where the weak density coupling exists. The procedure for this coupling may differ from solver to solver based on the type of numerical procedure used. Improper inclusion of the variable density can lead to false results. Different procedures are proposed for the explicit procedure based numerical solvers18, 19 and semi-implicit ones20. In the present work we used an implicit solver, and therefore developed and applied an adequate variable density approach. Some details of the implementation of this technique are presented in the following section. Handling of outlet convective boundary condition is also presented along with flame-wall interaction strategies.
To measure the influence of this methodology to capture the behaviour of the premixed flame the so-called VOLVO test case21 referred as Validation Rig 1 (VR 1) is investigated, as some experimental results are available. This case was already investigated by using large eddy simulation in the context of the G-equation approach by Ping22. We provide and discuss some comparisons between experimental and simulation results.
2 NUMERICAL PROCEDURE AND MODELLING
the Favre filtered continuity, Navier-Stokes and the progress variable equations. The two first equations are i 0 i i u t x ρ ρ ∂ ∂ + = ∂ ∂ (1) i i i
(
)
i i 2 j 3 j i i k i j ij ij j i j i j k u u p u u u u T t x x x x x x ρ ρ ⎡ρυ⎛∂ δ ⎞ ρ ⎤ ∂ + ∂ = −∂ + ∂ ⎢ ⎜ +∂ − ∂ ⎟ + ⎥ ⎜ ⎟ ∂ ∂ ∂ ∂ ⎢⎣ ⎝ ∂ ∂ ∂ ⎠ ⎥⎦ (2) where tildas ( )⋅ refer to Favre-filtered variables, whereas overbars ( )⋅ refer to simple spatial filtering (u=ρ ρu/ ). The effect of the unresolved subgrid scales is represented by the SGS stressi j ij u ui j u u
τ = − (3)
In the Smagorinsky model, the anisotropic part of the SGS-turbulent stress, a ij
τ , is related to the resolved strain-rate tensor Sij by
(
)
1/ 2 2 2( ) 2 a ij Cs S Smn mn Sij τ = − ∆ (4)where Cs represents the Smagorinsky coefficient. To compute this coefficient for complex geometries we introduce a test scale filter represented by a tilde in accordance with the Germano procedure. The purpose of doing this is to utilize the information between the grid- and test-scale filters to determine the characteristics of the SGS motion. The Smagorinsky coefficient can be then calculated dynamically using the following expression:
2 ij ij s ij ij L M C M M = (5) where 2 2 2 2 ij ij ij M = − ∆ S S + ∆ S S (6) k i i i i j ij i j ij ij L =u u −u u or =T −τ (7) k i i ij i j i j T =u u −u u (8)
∆ is a test filter width and S =
(
)
1/ 2
2S Smn mn
Different dynamic techniques can be applied to compute Eq. (5). Among these the Lagrangian dynamic model from Meneveau et al.16, is geometry independent and hence will be used for the complex geometries, considered in this work.
i i
( )
. .i 4 6 i(1 i) 16 6 t L u u L t S C C C uC C S t Sc µ ρ ρ ρ ρ π π ⎡⎛ ∆ ⎞ ⎤ ∂ − + ∇ ⋅ = ∇ ⎢⎜⎜ + ⎟∇ ⎥+ Ξ ∂ ⎢⎣⎝ ⎠ ⎦ ∆ (9) i i 1 u u or C C ρ ρ ρ ρ ρ τ τ − = = + (10)In equation (9) the last term represents the reaction term which is based on algebraic expression of flame surface density. The extra diffusion term (first term) in the r.h.s is added to preserve the correct flame propagation speed when turbulence is small. Ξ denotes flame wrinkling factor. In the present work we used the Charlette formulation10.
1 min , l l u S β δ ∆ ⎡ ⎛∆ ′ ⎞⎤ Ξ = +⎢ ⎜ Γ ⎟⎥ ⎝ ⎠ ⎣ ⎦ (11)
where Γ is an efficiency function, the laminar flame thickness, S is the laminar flame speed l and β is the power-law exponent. More details about this wrinkling model can be found in10.
By using equations (1), (9) and (10) one can find that the following relation for the divergence of velocity field
i i(1 i) . . 4 6 16 6 i L t u u L i u t u S C C C S x Sc µ τ ρ ρ π ρ π ⎡ ⎡⎛ ⎞ ⎤ ⎤ ∂ = ⎢∇ ∆ + ∇ + Ξ − ⎥ ⎢⎜⎜ ⎟ ⎥ ∂ ⎢⎣ ⎢⎣⎝ ⎠ ⎦ ∆ ⎥⎦ (12)
holds. This relation has been derived by replacing the iρCfrom the equation (9) by using the relation from equation (10) and using continuity equation (1). In the numerical procedure, the density was updated by using the relation equation (12) and the continuity equation (1). More details of this variable density formulation are going to be presented elsewhere23. In addition Equation (12) is also used for the calculation of convective outlet velocity. In general convective outlet boundary condition is defined by using
i i i . . 16 6 *( _ ) (1 ) 4 6 t L u t C S inlet v u u L S C Sc U outflow surface u ds dv C C S µ ρ τ π ρ ρ π ⎡ ⎡⎛ ∆ ⎞ ⎤⎤ ⎢∇ ⎢⎜⎜ + ⎟∇ ⎥⎥ ⎢ ⎢⎣⎝ ⎠ ⎦⎥ − = ⎢ ⎥ − ⎢+ Ξ ⎥ ⎢ ⎥ ⎣ ∆ ⎦
∫
∫
(14)The above relation allows then to calculate the convective outlet velocity U . In order to C fulfil the mass conservation through out the system, equation (12) will be further used. For this, two types of methods are used. One is additive correction and the other is a multiplicative one. For open boundaries additive correction ingeneral provides good results where as for the wall bounded flows multiplicative correction is known to give better results23.
3 RESULTS AND DISCUSSION
For the present investigation, the VOLVO configuration (VR1)21 has been considered as sketched in Fig.1.
Fig. 1 Bluff body test configuration
It is a bluff body stabilized flame configuration, and consists of a rectilinear channel with rectangular cross section, divided into an inlet section and a combustion section equipped with a two-dimensional triangular shape flame holder. The channel is 0.24m width and 0.12m high. The blockage ratio of the obstacle was 1/3. The primary intention of this rig is to investigate phenomena occurring in afterburners for jet engines.
In figure 2 mean axial velocity and stream lines contour have been showed, while the flame front can be seen in Figure 3. It has been observed that flame will be detached from the bluff-body if any wrinkling model is not used. Furthermore it has been observed that the wrinkling model plays a main role near the bluff-body region. It has been observed that without flame wrinkling model flame will detach from the flame holder. Figure 4 shows the axial velocity prediction along the centreline in the stream-wise direction. From this figure it is clear that velocity is over-predicted and the recirculation length is under-predicted. Ping22 also reported similar results.
Figure 5, 6 and 7 show the mean temperature along three sections behind the flame holder, (z=0.15m, z=0.35m and z=0.55m). The agreement between experiments and the numerical simulations is satisfactory. The discrepancies near the wall are mainly due to the numerical quenching of the flame near wall. However a proper quenching near the wall is possible if the wall temperature is provided. Unfortunately for the present case wall temperatures are not available. Normally quenching near the wall can be done by two ways. One is side wall quenching in which flame travels along the wall and the other one is head on quenching in which flame travels towards wall. In the present case, due to the unavailability of the wall temperature we have performed the side wall type quenching.
4 CONCLUSION & OUTLOOK
In the present paper we have performed the large eddy simulations of bluff-body stabilized premixed flame by using the flame surface density model with the help of an implicit numerical solver. For this purpose a simple technique for the density coupling has been used for the simulations along with a proper handling of the outlet convective boundary conditions. Comparisons between the experimental and numerical results have showed that the numerical method implemented is able to provide satisfactory agreement.
As future work species concentration such as CO2 and O2 are going to be predicted by including proper chemistry for the propane. Furthermore an extension of present work to partially premixed flames by including mixture fraction equation is undergoing.
ACKNOWLEDGEMENTS
For financial support we gratefully acknowledge the DFG (German Research Foundation). REFERENCES
1. Janicka J., and Sadiki. A., 2004, Proc. Combust. Institute, 30: 537-547.
2. Kerstein, A. R., Ashurst, W., and Williams, F. A., 1988, Phys. Rev. A, 37(7):2728– 2731.
3. Menon, S., and Jou, W. H., 1991, Combust. Sci. Technolo., 75:53–72.
5. Boger, M., Veynante, D., Boughanem, H., and Trouve´, 1998, A., Proc. Combust. Institute, 27:917–926.
6. Butler, T., and O’Rourke, P., 1977, Proc. Combust. Institute, 16:1503–1516. 7. O’Rourke, P., and Bracco, F. V., 1979, J. Comp. Phys. 33:185–203.
8. Angelberger, C., Veynante, D., Egolfopoulos, F., and Poinsot, T., 1998, Annual Research Briefs 1998, Center for Turbulent Research, NASA Ames/Stanford University, Stanford, CA.
9. Colin, O., Ducros, F., Veynante, D., and Poinsot, T., 2000, Phys. Fluids, 12(7):1843– 1863.
10. Charlette F., Meneveau C. and Veynante D., 2002, Combustion and Flame, 131:159– 180.
11. Charlette F., Meneveau C. and Veynante D., 2002, Combustion and Flame, 131:181– 197.
12. Pocheau A., 1994, Physical Review E, 49:1109-1122.
13. Boger M and Veynante D., 2000, Advances in Turbulence VIII, 449-451.
14. Germano M., Piomelli U., Moin P., and Cabot W. H., 1991, Phys. Fluids A, 3: 1760-1765.
15. Moin P., Squires K., Cabot W., and Lee S., 1991, Phys. Fluids A, 3:2746-2752.
16. Meneveau C., Lund T. S., and Cabot W. H., 1996, Journal of Fluid Mechanics, 319: 353-385.
17. Akula R. A., Sadiki A., And Janicka J., 2005, Turbulent Shear Flow Phenomena IV, 2: 595-600.
18. McMurtry PA, Jou WH, Riley JJ, and Metcalfe RW, 1986, AIAA Journal, 24(6): 962-970.
19. Boersma BJ, 1998, Annula Research Briefs 1998, Centre of turbulence research, Stanford University, 47-56.
20. Najm HN, Wyckoff PS, and Knio OM, 1998, J. Computational Physics, 143: 381-402. 21. Sjunnesson A., Olovsson S., and Sjoblom B., 1991, Volvo Flygmotor internal report,
VFA9370-308.
22. Ping W, 2005, Large eddy simulation of turbulent swirling flows and turbulent premixed combustion, Lund University.
Fig. 2 Mean axial velocity and stream lines contours
Fig. 3 Flame front
Fig. 4 Axial velocity prediction along the
centreline in the stream-wise direction. Fig. 5 temperature behind the bluff-body at z=150mm
Fig. 6 temperature behind the bluff-body at z=350mm
