Large Eddy Simulation of Vortex Breakdown/flame Interaction: Related Articles

Large eddy simulation of vortex breakdown/flame interaction
Christophe Duwig and Laszlo Fuchs
Citation: Phys. Fluids 19, 075103 (2007); doi: 10.1063/1.2749812

View online: http://dx.doi.org/10.1063/1.2749812
View Table of Contents: http://pof.aip.org/resource/1/PHFLE6/v19/i7
Published by the American Institute of Physics.
Related Articles
Instability of vortex pair leapfrogging
Phys. Fluids 25, 014107 (2013)
Bistability and hysteresis of maximum-entropy states in decaying two-dimensional turbulence
Phys. Fluids 25, 015113 (2013)
Experimental study on side force alleviation of conical forebody with a fluttering flag
Phys. Fluids 24, 124105 (2012)
Effects of passive control rings positioned in the shear layer and potential core of a turbulent round jet
Phys. Fluids 24, 115103 (2012)
Vortex formation and instability in the left ventricle
Phys. Fluids 24, 091110 (2012)
Additional information on Phys. Fluids

Journal Homepage: http://pof.aip.org/
Journal Information: http://pof.aip.org/about/about_the_journal
Top downloads: http://pof.aip.org/features/most_downloaded
Information for Authors: http://pof.aip.org/authors
Downloaded 14 Mar 2013 to 131.170.6.51. Redistribution subject to AIP license or copyright; see http://pof.aip.org/about/rights_and_permissions
PHYSICS OF FLUIDS 19, 075103 共2007兲
Large eddy simulation of vortex breakdown/flame interaction

Christophe Duwiga兲 and Laszlo Fuchs
Division of Fluid Mechanics, Department of Energy Sciences, Lund University, SE-22100 Lund, Sweden
共Received 24 October 2006; accepted 28 April 2007; published online 16 July 2007兲
The dynamics of a swirl-stabilized premixed flame is studied using large eddy simulation 共LES兲. A
filtered flamelet model is used to account for the subgrid combustion. The model provides a
consistent and robust reaction-diffusion expression for simulating the propagation of turbulent
premixed flames correctly. The numerical results were found to be relatively insensitive to small
changes in the inflow boundary conditions and to the numerical mesh employed. Furthermore, the
results were found to agree well with the available experimental data both for velocity and scalar
fields. In addition, unsteady flame features 关i.e., precessing vortex core 共PVC兲兴 were identified and
compared with experimental data. The agreement between LES results and experimental data, in
terms of flame dynamics, was also good. Increasing swirl did not affect the flame strongly but a
decrease of swirl number was shown to change the flame shape and suppress the PVC. The PVC and
flame dynamics were studied using proper orthogonal decomposition 共POD兲 allowing us to identify
and isolate the PVC from smaller-scale turbulence. The POD results indicate that the PVC
corresponds to a helical wave consisting of two counter-rotating helices. A dynamical reduced
model was also derived do describe the flame response to the PVC. © 2007 American Institute of
Physics. 关DOI: 10.1063/1.2749812兴
INTRODUCTION gases that provide heat to the fresh gases and can thereby
maintain the combustion at rather lean conditions.1,2,7,9 How-
Modern gas turbines 共GT兲 utilize lean premixed burners
ever, it adds an additional complexity since thermal expan-
in order to achieve low levels of NOx emissions.1,2 The
sion affects the breakdown, i.e., it increases the velocity of
flame is kept in the combustor by a relatively strong swirling
the backflow. For example, Selle et al.4 and Schneider et al.7
motion of the 共lean兲 fuel/air mixture. The lean, swirling mix-
ture is highly susceptible to perturbations both in terms of reported that the flame can suppress large-scale coherent
flow and flame unsteadiness as well as flame blowout.2 The structures appearing during vortex breakdown while Syred8
dynamics of the flame is the result of the interactions among reported an increase of the coherent structure frequency with
the fluctuations in fuel-air concentration, heat-release, mix- the fuel/air equivalence ratio. It is clearly a nonlinear two-
ture velocity, and pressure.3–5 The main difficulty lies in the way coupling since the swirling motion enhances the shear-
nonlinear behavior of several of the controlling phenomena layer instability and thereby increases the flame surface
and the interaction among them 共e.g., chemistry and turbu- area.6
lence兲. For example, turbulent mixing instabilities may lead Modeling and understanding the vortex breakdown and
to concentration and heat-release fluctuations. Such heat- its interaction with combustion is then a key issue in flame
release fluctuations form an acoustic source that may en- stabilization. The main difficulty is the unsteady behavior of
hance the fluctuations in equivalence ratio. Several authors5,6 this type of flow:7–10 Large structures, resulting from vortex
showed the importance of large-scale flow structures in the breakdown and the swirling shear layers, directly affect the
flame/flow interaction. This mechanism is of relevance for flame stabilization leading to heat-release fluctuations and
the GT combustors as large scales play an important role in possibly to combustion instabilities. An example of the large-
flame stabilization 共e.g., through vortex breakdown兲. scale unsteady motions arising during vortex breakdown is
Vortex breakdown is commonly utilized to anchor and the precessing vortex core 共PVC兲.8 Experimental7,9–12 and
stabilize flames.1,2,7 This vortex breakdown is the result of numerical4,13 works have reported that vortex breakdown
the swirling motion given to a jet. The critical level of swirl, might result in an off-axis precession of the central recircu-
for attaining the vortex breakdown, depends strongly on the lation zone. The recirculation zone undergoes a periodic ro-
burner geometry and on the radial distribution of the velocity tational and translational motion around the axis and the in-
field.8 At vortex breakdown, a low-pressure region appears
stantaneous flow field is far from being axisymmetric.
around the axis region close to the expansion. If the swirling
However, time averaging the flow also maintains the geo-
motion is strong enough, the longitudinal pressure gradient
metrical axisymmetry. If a flame is stabilized by the recircu-
induces an axial backflow and a stagnation point. In many
lation zone, the vitiated gases will be engulfed into the recir-
combustion systems, the axial backflow consists of hot burnt
culation zone and undergo similar precession.9 Consequently,
a兲 the PVC by itself can act as a source of excitation for the
Author to whom correspondence should be addressed. Present address:
HTAS, Lyngby, Denmark. Telephone: ⫹46 46 222 31 84. Fax: ⫹46 46 222 flame front and thereby become a source of thermoacoustic
47 17. Electronic mail: [email protected] fluctuations.
1070-6631/2007/19共7兲/075103/20/$23.00 19, 075103-1 © 2007 American Institute of Physics
075103-2 C. Duwig and L. Fuchs Phys. Fluids 19, 075103 共2007兲
The issue is therefore to model and understand the inter- In LES, the averaging operation corresponds to a low
action between the flame and the vortex. It is necessary to pass spatial filtering and leads to the unclosed set of PDE
capture in time and analyze the 3D flame stabilization pro- presented above. The main role of the unclosed or subgrid
cess. Such an analysis is carried out in this paper using the scale 共SGS兲 terms is to account for the interaction between
large eddy simulation 共LES兲. First, we present the different the resolved and the unresolved scales. In the momentum
components of the computational model, i.e., LES and the equations, the SGS terms should account for the dissipative
combustion model. Second, a sensitivity analysis of the re- character of turbulence on the small 共unresolved兲 scales as
sults to numerical mesh density and inflow boundary condi- well as for the transfer of energy among the resolved and
tions is presented. The analysis will assess the influence of unresolved scales. Corresponding properties may be attrib-
imperfectly known quantities 共inflow conditions兲 and not re- uted to the SGS in the energy- and the species-transport
solved turbulent structures upon the result. Comparisons equations.
with experimental data are also performed in order to relate A computational grid can support only Fourier compo-
the present results to physics. Finally, we seek to analyze the nents that have longer wavelengths than the grid size. Thus,
large-scale turbulent/flame motions. Since these fluctuations a dependent variable represented on a grid together with a
are smoothed out by time averaging, we focus on different discrete approximation for the derivatives leads to an implicit
statistical quantities obtained through proper orthogonal de- filtering. If no explicit SGS terms are added, the numerical
composition 共POD兲 of the flow and density fields. Finally, scheme should account at least for the small-scale dissipa-
the POD analysis is used to derive a reduced model repre- tion. This is attained for any numerically stable scheme and
senting the large-scale flame motions. is referred as implicit LES 共ILES兲 and also monotonically
integrated LES 共MILES兲 in the literature.14–17 However, nu-
PREMIXED COMBUSTION MODELING FOR LES merical schemes that are too dissipative are highly inappro-
Governing equations and momentum priate for LES since in addition to dissipation on small
equations closure scales, they may also be dissipative on the larger resolved
scales. This effect can be avoided by choosing appropriate
The basic partial differential equations 共PDE兲 describing 共higher-order兲 discretizations; for instance, Grinstein and
the motion of a reacting fluid are the conservation of mo- Fureby14 advocate the use of FTC 共flux-corrected transport兲
mentum, mass, and energy. The system of equations has to schemes while Margolin et al.17 recommend high-order
be closed with an equation of state. Since several species are monotonicity preserving schemes. In the following, we use a
involved, one has to add one transport equation for each of fifth-order weighted essentially nonoscillatory 共WENO兲
the mixture species including respective sink/source terms. scheme18 that ensures numerical stability but does not affect
In addition, a turbulent flow is inherently unsteady and one
significantly the resolved scales 共i.e., low dissipation and dis-
seeks to describe the statistical properties of the dependent
persion characteristics兲.
variables by using a suitable averaging method. If the basic
equations are averaged 共independently of the averaging Modeling of premixed combustion
method兲, the equations can be written 14
Large eddy simulation has also been shown to be very
⳵␳
¯
+ ⵜ · 共¯␳ũ兲 = 0, 共1兲 suitable for handling unsteady coherent structures in
⳵t combustors.4,13,19,20 However, there are several open issues
related to LES and there is only a limited theoretical foun-
⳵␳
¯ ũ dation for modeling the SGS terms related to combustion.
+ ⵜ · 共¯␳ũũ兲 = − ⵜP̄ + ⵜ · 共− ␳uu + ¯␳ũũ + ␮ ⵜ ũ兲, 共2兲
⳵t Incorporating combustion chemistry into LES involves find-
ing a suitable reaction mechanism and solving the filtered
⳵␳
¯ Ỹ i species and energy equations 共3兲 and 共4兲. Appropriate reac-
+ ⵜ · 共¯␳ũỸ i兲 = ⵜ · 共− ␳uY i + ¯␳ũỸ i + ␳Di ⵜ Ỹ i兲 + ␻
¯ i, tion mechanisms, either detailed for simple fuels or reduced
⳵t
for practical fuels, may involve tens of species and hundreds
共3兲 of reaction steps resulting in a large system of stiffly coupled
equations that remain to be solved. Another problem is that
⳵␳
¯ T̃ Eqs. 共3兲 and 共4兲 contain the filtered reaction rates, which are
+ ⵜ · 共¯␳ũT̃兲 = ⵜ · 共− ␳uT + ¯␳ũT̃ + ␳DT ⵜ T̃兲 + ␻
¯ T , 共4兲
⳵t highly nonlinear functions of temperature and species con-
centrations. The closure of Eqs. 共3兲 and 共4兲 is still an open
¯␳ = 冉冊
P
RT
, 共5兲
issue and several alternatives have been proposed. Examples
of such approaches includes 共i兲 implicit LES,21 共ii兲 thickened
flame models 共TFM兲,22 共iii兲 linear eddy models 共LEM兲20
where u is the velocity vector, P is the pressure, T is the with embedded 1D grids, 共iv兲 flamelet models 共FM兲,19,23–27
temperature, ␳ is the density, ␮ is the viscosity, R is the and 共v兲 filtered density function 共FDF兲 models;28 see Ref. 14
specific ideal gas constant, Di is the diffusivity of species i, for further details.
Y i is the mass fraction of species i, and ␻i is the reaction Based on experience and computing efficiency, we con-
term of species i. The bar denotes the averaging operation centrate on the FM class of models and consider a single
and the tilde the density weighted averaging. scalar: the so-called progress-variable 共c兲, which is used as
075103-3 LES of a swirling flame Phys. Fluids 19, 075103 共2007兲
FIG. 1. Production terms ⌸C represented in the density-averaged c space.
the unique coordinate for describing the combustion process. ⳵␳

¯ c̃
Among the different candidate closures,24,26,27 we choose the + ⵜ · 共¯␳ũc̃兲 = ␳0D0ⵜ2c̃ + ␲共c̃,D0兲. 共7兲
⳵t
last one. This choice is supported by Fureby,29 who compares
different premixed combustion models in the LES frame- Note that ␲ is defined to be the RHS of 共7兲 minus the diffu-
work. It was found that all tested closures enable reproduc- sion term. The ␲ term contains information related to the
ing the experimental data well in terms of mean and root- nongradient subgrid transport as well as the creation 共source兲
mean-square 共rms兲 fields, but the closure proposed by Duwig of c. Furthermore, the ␲ term 关defined by Eq. 共7兲兴 is distrib-
and Fuchs27 performed best for prediction of the pointwise uted over one to two filter lengths. Assuming that the source
temperature probability density function 共PDF兲. In addition, of c 共i.e., the reaction layer of the premixed flame兲 is thin
this closure was shown to be suitable for simulating turbulent compared to the grid size 共valid in the flamelet and thin
flames.27,29,30 reaction zone regimes兲, Duwig and Fuchs27 suggested to
model ␲ using a Gaussian distribution in space so that ␲
c-equation closure = ␣ / ⌬ · 共6 / ␲兲0.5 · exp共关−6 · 共x / ⌬兲2兴, where x is the distance to
We define the progress variable as c = 共T − Tu兲 / 共Tb − Tu兲, the reaction layer, ␣ is the integral burning rate, and ⌬ is
where Tu and Tb are the unburned- and burned-gas tempera- the LES filter length. The burning rate is related to the lami-
ture, respectively.14 We have c = 0 in the fresh gas and c = 1 in nar flame speed 共SL兲 and the subgrid flame wrinkling 共⌶兲
the burnt gas. If local variations of the fuel/air equivalence so that ␣ = ␳uSL⌶. Further, a nondimensional parameter
ratio are small 共perfectly premixed combustion兲, the density a = 共␳uSL⌶⌬兲 / ␳0D0 emerges in 共7兲. Finally, a relation x共c兲 is
filtered c equation is derived from 共4兲 and reads needed to close Eq. 共7兲. In the limiting case of a quasiplanar
flame, the analysis of Duwig and Fuchs27 provides an expres-
sion for the term ␲⌬ / ␣ = ⌸共c , a兲 that is stored in lookup
⳵␳
¯ c̃ tables.27 The closed c equation reads
+ ⵜ · 共¯␳ũc̃兲 = ⵜ · 共− ␳uc + ¯␳ũc̃兲 + ⵜ · 共␳Dth ⵜ c兲 + ␻c .
⳵t ⳵␳
¯ c̃ ␳uSL⌶⌬ 2 1
+ ⵜ · 共¯␳ũc̃兲 = ⵜ c̃ + ␳uSL⌶ ⌸c共c̃,a兲. 共8兲
共6兲 ⳵t a ⌬
Figure 1 shows the production term ⌸c mapped in a 共c , a兲

One may observe the unclosed terms on the right-hand side space. For high a values, ⌸c tends toward a parabolic shape
共RHS兲. The first term is the subgrid scale transport term while for a ⬃ 1 its maximum shifts toward c ⬃ 0.8. For
共SGSC兲. The second term is related to the heat diffusivity and a ⬃ 1, the filter effect on the flame front is relatively low so
the last term is the filtered source term ␻¯ c. Note that the that ⌸c approaches the laminar reaction rate.
combination of these three terms results in a reaction- Equation 共8兲 requires an estimation of the subgrid wrin-
transport balance leading to the flame propagation, i.e., the kling factor or a subgrid flame speed SL⌶. These quantities
flame dynamics. Following the analysis of Duwig and are widely used in the framework of combustion
Fuchs,27 the three unclosed terms of Eq. 共6兲 are modeled by modeling22–27 and several expressions or correlations have
a reactive-diffusive balance ensuring correct flame propaga- been proposed for modeling them using the subgrid velocity
tion. The RHS of Eq. 共6兲 is modeled by a diffusion term plus vector fluctuation u⌬. Fureby26 tested the influence of differ-
a production term ␲ as follows: ent subgrid wrinkling correlations and concluded that their
length in the present paper. Downstream of the nozzle, the

swirling jet discharges into a combustion chamber of diam-
eter 0.11 m. The sudden expansion leads to vortex break-
down providing flame stabilization. The frame consists of
three spatial directions with the streamwise direction being
denoted by X. The origin of the X axis has been chosen at the
expansion plane so that X ⬎ 0 is located in the combustion
chamber while X ⬍ 0 lies in the premixing tube 共i.e., X = 0 is
FIG. 2. Schematic of the swirl generator, premixing tube, contraction, sud-
at the expansion兲.
den expansion, and combustion chamber 共Ref. 10兲. All dimensions are in Experimental investigations10–12 were performed using
mm. this setup both for reacting and nonreacting cases. Anacleto
et al.10 studied extensively this swirling flame and reported
velocity as well as temperature measurements. We choose to
effect on the solution is relatively limited. He advocates the concentrate on a particular case exhibiting a strong PVC for
use of a relatively simple power law and uses the fractal which detailed measurements are available. The operating
concept. It reads conditions are as follows:
⌶ = MAX 冋冉冊册
⌫·
兩u⌬兩
SL
D−2
,1 , 共9兲
•
•
Propane/air equivalence ratio: 0.5.
Unburnt gas temperature: Tu = 573 K.
• Bulk velocity in the nozzle: U0 ⬃ 59 m / s.
where D is the fractal dimension and ⌫ is an efficiency func- • Reynolds number: Re⬃ 81000; Mach number: M ⬃ 0.12.
tion presented by Meneveau and Poinsot.31 The fractal di- • Blade angle: 50°.
mension needs to be modeled, and following Fureby,26 we • Swirl number downstream of the swirler 共as defined in
use Ref. 10兲: S = 1.05.
2.05 2.35
D= + . 共10兲 The experimental measurements are summarized in Ref.
1 + 兩u⌬兩/SL 1 + SL/兩u⌬兩
10. They include laser Doppler anemometry 共LDA兲 measure-
Further, several authors19,24–26 use a transport equation to ments and propane mass fraction and temperature measure-
compute the SGS turbulent kinetic energy and the subgrid ments using intrusive techniques.10 The present study fo-
velocity fluctuation. Another alternative is to estimate the cuses on the interaction between the flame and the PVC.
subgrid velocity fluctuation from the resolved velocity field. Therefore, we should present the fuel mass fraction as well
As pointed out by Colin et al.,22 the velocity fluctuation as temperature fields as progress variables that can be com-
should not include thermal expansion effects. They consider pared to the simulations.
the rotational part of the resolved velocity field only and
relate the subgrid fluctuation to the smallest resolved scales.
Consequently, the subgrid velocity vector fluctuation u⌬ is
given by
NUMERICAL METHODS
u⌬ = C2⌬3ⵜ2共ⵜ ⫻ ũ兲, 共11兲
LES code
where C2 is a constant. Colin et al.22 recommend taking
C2 ⬃ 2. The present flamelet formulation has been implemented
Finally, the local filtered density is computed from the in a Cartesian finite-difference LES code solving the low-
temperature using an ideal gas law in the following way: Mach number formulation of the Navier-Stokes
equations27,32 on staggered grids. The spatial derivatives are
␳T ␳ uT u ␳ uT u discretized with a fourth-order centered scheme except for
¯␳ = ⬇ ⬇ , 共12兲
T̃ T̃ Tu + 共Tb − Tu兲c̃ the convective terms that are treated with a fifth-order
WENO scheme.18 The time integration is done by a fully
where the subscripts u and b denote the unburnt and burnt implicit scheme using a second-order upwind finite-
states, respectively. difference approximation for the time derivatives. The com-
putational grid is made of a system of locally refined grids.32
LISBON SWIRLING COMBUSTOR The local refinements are added gradually in regions with
expected 共or computed兲 large gradients.32 The implicit solver
Geometry and available experimental data
uses a multigrid procedure for enhanced efficiency. The 1D
Figure 2 shows the swirl combustor located at the laminar flame speed and burnt/unburnt gas properties are
Lisbon Technical University.10 The swirler’s blade angle can computed using the software Cantera33 together with a pro-
be adjusted to fix the tangential/axial momentum ratio. The pane oxidation scheme.34
swirling annular jet enters a contraction and discharges into a The time step used is 1.5⫻ 10−6 s chosen such that the
premixing pipe of diameter 0.05 m. The pipe ends with a inlet Courant number is ⬃0.1. The mean values have been
converging/diverging nozzle. The minimum diameter within computed over 15000 to 25000 time steps corresponding to
the nozzle is D = 0.04 m and will be used as a characteristic ⬃4 flow-through times the whole computational domain.
TABLE I. Summary of the simulations.
Run Grid Re S Utang / Uax Note
1 Medium 81000 1.05 ⬃0.9 LES ref-Base case

2 Coarse 81000 1.05 ⬃0.9 —
3 Fine 81000 1.05 ⬃0.9 —
4 Coarse 81000 ⬃0.95 ⬃0.8 LES ⫺
5 Coarse 81000 ⬃1.15 ⬃1.0 LES ⫹
6 Medium 81000 0.45 ⬃0.4 Low swirl
7 Medium 81000 1.35 ⬃1.2 High swirl
POD analysis diameter 1.25D and length 4.14D discharging into a cylindri-
cal combustor of diameter 2.75D and length 8.4D. The pre-
Extracting information from the LES data is an essential
mixing pipe has a converging-diverging nozzle 共the diameter
issue for characterizing the turbulent flame and to enable one
at the contraction is 1D兲 before the sudden expansion. The
to understand the mechanisms controlling the specific prob-
flow enters the premixing pipe through an annulus of inner
lem. Traditionally, a turbulent flow is characterized by its
diameter 0.75D and outer diameter 1.25D. A conical bluff
mean and fluctuation 共rms兲 quantities, pointwise spectra, and
body is located inside the annulus. The length of the body
PDFs. An additional tool for data reduction that can be used
is ⬃0.6D. It is worth noting that the combustor is axisym-
to gain improved insight into the problem is proper orthogo-
metric but that the symmetry is not utilized in the LES com-
nal decomposition 共POD兲.35–38 The POD projects the turbu-
putations since we seek to capture spatially 3D turbulent
lent flow field on an orthonormal vector base that maximizes
structures.
the turbulent kinetic energy content for any subset of the
The combustion chamber is considered to operate in a
base. It allows an accurate description of the turbulence data
perfectly premixed mode. There are no available measure-
using only a limited number of modes.35 Given a vector Q
ments of the homogeneity of the fuel/air mixture, but since
containing the field variables and an orthonormal base ⌽, the
the fuel is injected far upstream of the flame, it is reasonable
expansion reads
to approximate the fuel/air mixture by an ideally mixed sys-
N tem. The propane/air equivalence ratio was 0.5. The fresh
Q共x,t兲 ⬇ Q 共x,t兲 = 兺 ai共t兲⌽i共x兲.
N
共13兲 gas temperature is 573 K and the outlet temperature about
i=0 1724 K. The Reynolds number based on the pipe diameter at
the contraction is Re⬃ 81 000. The swirl number down-
Note that the approximation QN of the turbulence data set Q stream of the guiding vanes was varied between S ⬃ 0.45 and
converges to Q when N goes to infinity and that i = 0 corre- 1.35. The global Karlowitz number14 is estimated to be
sponds to the averaged field. The base vectors are computed Ka⬃ 2 – 4 in the flame region. The flame lies in the so-called
so that they satisfy the eigenvalue problem,35,36 thin-reaction regime; see, e.g., Ref. 14. However, Duwig and
Fuchs27,30 showed that the present model can approximate
具Q共x,t兲 · QT共x,t兲典⌽共x兲 = ␭⌽共x兲, 共14兲
successfully a flame in the thin-reaction regime. In addition,
where the superscript T denotes the transposed of the vector the numerical results are compared with the available experi-
and 具.典 is the time-averaging operator. The eigenvalue ␭i mental data to assess their accuracy.
characterizes the turbulent kinetic energy content of the
mode i. For practical reasons, Sirovich’s method of
snapshots36,37 is used to perform the POD. It consists of
transforming the problem into an equivalent eigenvalue
problem of much smaller size. In this paper, we used the 3D
data of the density field plus the three velocity components
and some 400 snapshots for performing the POD. It has been
found that increasing the number of snapshots beyond ⬃300
does not influence the first five POD modes significantly. The
extracted POD modes can be used further to derive dynami-
cal models based on ordinary differential equations 共ODEs兲
enabling capturing the flow/flame dynamics without solving
the full PDE system 共e.g., Ref. 39兲.
Computational geometry
The geometry, corresponding to a portion of the experi- FIG. 3. Systematic error due to time averaging for density variable flows.
mental rig presented in Fig. 2, contains a premixing pipe of Streamwise evolution of the L⬁ norm of SE along radii.
FIG. 4. Grid sensitivity of the normalized mean axial velocity: radial profiles at x / D = 0.25, 0.5, 1.25, and 2.5.
Three computational grids have been used in the present Boundary conditions
work as presented in Table I. For the present flow, we esti-
The boundary conditions are essential for all computa-
mate the Taylor turbulent microscale to be ␭Taylor ⬃ D / 50,
tional results and in particular for LES. Swirling flows in
which gives an idea for a suitable mesh spacing. However,
general are highly sensitive even to small variations in the
this global estimate does not assess on the grid influence
upon the results, which has to be done via a proper sensitiv- temporal or the spatial form of the inlet conditions; see, e.g.,
ity analysis. Three different grids, coarse, medium, and fine, Ref. 27. One may approximate the flow using a fully devel-
are used, with 1.6, 1.9, and 3.4 million computational cells, oped turbulent pipe flow or any other well defined 共mean兲
respectively. The coarse, medium, and fine grids have a mesh velocity profile and adding to that the 共artificially generated兲
size in the flame region of h ⬃ D / 40, h ⬃ D / 46, and instantaneous turbulent fluctuations. In fact, due to the un-
h ⬃ D / 52, respectively. All grids utilize the option of having certainty in experimental conditions, it is not very construc-
local refinements in order to have higher resolution in the tive to strive for the exact reconstruction of the experimental
premixing pipe and in the flame region while using a coarser conditions. Instead, it is more important to assess the sensi-
mesh further downstream. Thus, the mesh spacing close to tivity of the results to the uncertainty in the boundary condi-
the outlet of the domain is four times larger than in the flame tions. In the present study, we start the LES computation at
region. the beginning of the premixing tube. Steady inflow bound-
FIG. 5. Grid sensitivity of the normalized axial velocity fluctuation rms: radial profiles at x / D = 0.25, 0.5, 1.25, and 2.5.
FIG. 6. Grid sensitivity of the Favre-averaged progress variable: radial profiles at x / D = 0.25, 0.5, 1.25, and 2.5.
aries are enforced and the corresponding mean profiles are pressure boundary conditions are needed since we use a stag-
obtained from a separate Reynolds averaged Navier-Stokes gered grid formulation with boundaries located at the cell
共RANS兲 computation computed on an extended domain. The faces.32
use of a steady inflow boundary is supported by the fact that The different computational cases are summarized in
strong shear in the premixing pipe generates turbulence so Table I.
that the flow is strongly turbulent before it reaches the flame.
We assess the influence of the uncertainty related to the mod- RESULTS
eled inflow condition on the results. This is done by chang-
ing the inflow tangential velocity 共Utang兲 by ±10% keeping In the following, quantitative and qualitative results are
the axial velocity Uax 共i.e., the mass flow兲 constant. These presented. The quantitative results presented below are
cases are referred to as LES+ and LES−. normalized using the bulk velocity at the contraction
At the outlet, all variables are assumed to have a zero 共U0 = Ubulk = 59 m / s兲 and the pipe diameter at the contrac-
gradient with a correction 共if needed兲 to ensure global mass tion, i.e., D = 0.04 m. The vortex core has been visualized
conservation. Since the heat losses at the wall are not known, using a criterion based on the second largest eigenvalue of
solid adiabatic walls and nonslip conditions are imposed. the second invariant of the velocity derivative tensor pro-
The impact of the adiabatic walls is expected to be low since posed by Jeong and Hussain40 共the so-called ␭2 technique兲.
we investigate the flame stabilization in the central recircu- In the following figures, the vortex core corresponds to a
lation zone, which is far away from the wall. No explicit region where the eigenvalue is negative.
FIG. 7. Sensitivity of the normalized mean axial velocity to the inlet boundary condition: radial profiles at x / D = 0.225, 0.5, 1.125, and 3.625.
FIG. 8. Sensitivity of the normalized axial velocity fluctuation rms to the inlet boundary condition: radial profiles at x / D = 0.225, 0.5, 1.125, and 3.625.
Before presenting the LES results and comparing them

to the experimental data, we consider some of the effects of
averaging. An important but often not recognized shortcom-
SE ⬇ 冓冔 ␳Q
¯␳
−
具␳Q典
具¯␳典
. 共15兲
ing in the use of experimental data for comparison with nu- In the present case, we used the LES data to estimate the SE.
merical results lies in the averaging procedure. Turbulent Figure 3 shows the estimate for three quantities: mean axial
flow variables are meaningful only if statistical quantities velocity, rms of the axial velocity fluctuations 共both being
共usually mean or rms兲 are considered. The set of equations normalized by the bulk velocity U0兲, and the mean progress
共1兲–共5兲 resolves density weighted filtered fields that lead to variable. SE is zero in most of the domain since the flame
Favre-averaged/filtered quantities after time averaging. How- region is relatively confined 共0 ⬍ x / D ⬍ 2兲. SE reaches a
ever, experimental data have been obtained by a different peak at x / D ⬃ 0.5 corresponding to the location of large den-
averaging technique. sity fluctuations. SE might be up to ±0.18 for c and
The uncertainty resulting from the time-averaging proce- ±0.10 U0, ± / −0.04 U0 for the mean axial velocity and its
dure should be estimated if the flows are subject to high- fluctuation, respectively.
density gradients. Denoting the time or the ensemble averag- Another estimate of the uncertainty is carried out by
ing operation by 具.典, the systematic error 共SE兲 on a variable Q computing the mass flux from the LDA data. This procedure
is estimated as is only applicable when density is constant along the radius.
FIG. 9. Comparisons with experimental data 共LDA兲共Ref. 10兲: normalized mean axial velocity radial profiles at x / D = 0.0, 0.225, 1.125, and 3.625.
FIG. 10. Comparisons with experimental data 共LDA兲共Ref. 10兲: normalized mean radial velocity radial profiles at x / D = 0.0, 0.225, 1.125, and 3.625.
In the present case, temperature measurements indicate that size. Figure 5 shows the axial velocity fluctuation rms along
it is restricted to the region x / D ⬎ ⬃ 2. For instance, at the the same lines. Again the discrepancies due to the grids are
location x / D = 3.625, the mass flow computed from the LDA small. The situation is similar also for the more sensitive
data is about 15% lower than the one reported in the original Favre-averaged progress variable along radial lines, as
paper.10 shown in Fig. 6. The difference between the predictions is
marginal, except at the location x / D ⬃ 0.5. At the centerline,
Sensitivity analysis the fine grid prediction departs strongly from the others. It
LES resolves only turbulent scales that are larger than predicts a thinner flame brush so that c reaches ⬃0.8 at the
grid spacing. Hence it is essential to assess the effect of this centerline while it reaches only ⬃0.3 in the other cases.
incomplete resolution not only by using high-order numeri- However, the overall influence of the computational grid was
cal schemes and estimating the resolved scale size, but also found to be low and all three grids predict a very similar flow
by grid-sensitivity study. For the class of problems consid- pattern 共vortex breakdown and central recirculation zone兲
ered in this study, it is difficult to make large grid variations and turbulence levels. Consequently, the coarsest grid can be
and the focus has been on changing the swirling jet resolu- used to conduct sensitivity analysis of different parameters
tion. Figure 4 shows the mean axial velocity radial profiles in since the sensitivity to parameters is not expected to be
the combustor. The differences between the three grids are higher with the finest grid.
not significant and all three computations predict the same The sensitivity to the axial/tangential momentum ratio
central recirculation zone both in terms of magnitude and imposed at the computational inlet has been assessed. Fig-
FIG. 11. Comparisons with experimental data 共LDA兲共Ref. 10兲: normalized axial velocity fluctuation rms radial profiles at x / D = 0.0, 0.225, 1.125, and 3.625.
FIG. 12. Comparisons with experimental data 共LDA兲共Ref. 10兲: normalized radial velocity fluctuation rms radial profiles at x / D = 0.0, 0.225, 1.125, and 3.625.
ures 7 and 8 depict the mean axial velocity and rms of the Comparisons with experimental data: Case S = 1.05
axial velocity fluctuations along radial lines in the combustor
LDA data9 are used for comparison with the present LES
for the different inflow conditions. The effect of variations in
data. Figures 9 and 10 depict mean axial and radial velocity
the inflow tangential velocity 共i.e., swirl兲 is found mainly in
profiles at several axial locations. The overall agreement be-
the rms values and mostly near the inlet. In addition, the
tween the LES results and the LDA data is good and both
sensitivity of the mean progress variable field to the inflow
boundary 共not presented here兲 was not found to be signifi- data sets describe the same flow patterns. The agreement
cant. This behavior is expected since the increase in shear close to the expansion 共x / D ⬍ 0.5兲 is below experimental
increases turbulence production near the inlet, which mani- uncertainty, showing that the LES captures the central recir-
fests in higher rms levels in the swirling jet shear layer. The culation zone as well as the weak outer recirculation accu-
effect of the inlet conditions reduces further downstream as rately. Further downstream 共locations x / D ⬃ 1.125 and
the local generation of turbulence dominates and the contri- 3.625兲, the LES seems to underestimate the back flow and
bution of turbulent transport reduces. Thus, by placing the exhibits somewhat smaller spreading of the swirling jet.
inflow boundary condition at more than four diameters up- However, the uncertainty of the velocity measurements in
stream of the expansion, the uncertainty in the inlet condi- this region 共estimated previously to be ⬃15%兲 suggests that
tions is shown to have a minor effect on the results. the discrepancy is less significant. Figures 11 and 12 show
FIG. 13. Comparisons with experimental data: progress variable radial profiles at x / D = 0.0, 0.25, 0.5, 0.75, 1, and 1.5. The experimental profiles are extracted
from temperature 共circle兲 and from the fuel mass fraction Xf 共triangle兲共Ref. 10兲.
FIG. 15. Visualizations of the vortex core as ␭2 isolevel 共Ref. 40兲共solid
gray surface兲 together with the flame surface 共isolevel c = 0.5/transparent
surface兲.
progress variable fields兲 with the LES results. The overall

agreement between the LES and the measurements is good
considering the systematic error induced by the averaging
procedure. In addition, the discrepancy between the two sets
of experimental data provides an estimate for the experimen-
tal uncertainty. However, the two experimental progress vari-
able fields agree reasonably well and show similar tends.
These trends are well captured by the LES tool. In particular,
comparing Figs. 9 and 13, it is observed that the central
recirculation zone extends upstream of the expansion plane
共x / D ⬍ 0兲 whereas the flame brush starts at x / D ⬃ 0.2. This
FIG. 14. Snapshots of the flame surface captured as isosurfaces of the is surprising since the negative velocity at the centerline is
progress variable c. Top, isolevel c = 0.6; bottom, isolevel c = 0.2. expected to drag the flame upstream of the expansion plane.
Actually, one should stress that averaged fields are only sta-
tistical representations of the flow. For instance, the averaged
the rms of the axial and radial velocity fluctuations along data are axisymmetric while the real flame is definitely not
several lines in the combustion chamber. The LES predic- symmetric, as shown in Fig. 14. In particular, the visualiza-
tions appear to overestimate, by a factor 2, the axial velocity tions of the flame surface exhibit strong large-scale wrinkling
fluctuation levels close to the expansion but capture the ra- with a helical tail of the flame rotating around the axis.
dial velocity fluctuations well. Since the mean and rms fields Figure 15 shows the flame surface as well as the vortex
are different statistical views of the same turbulent flow
fields, it is expected that an overestimation of the turbulent
levels would affect all components. As pointed out previ-
ously, the systematic error due to the averaging is large in the
flame region where both intermittency and density gradients
are important. It explains why the discrepancies between the
LDA and the LES results are observed for one component
only. Thus, further downstream, the LES predicts the axial
velocity fluctuation levels correctly. Considering the radial
velocity fluctuations, the LES results perform well capturing
both the amplitude and the spatial distribution of the fluctua-
tions. In particular, close to the expansion 共x / D ⬍ 0.5兲, both
the LES results and the LDA measurements capture a strong
peak on the axis. Fernandes et al.12 stressed that this particu-
lar peak corresponds to an off-axis motion, i.e., to the PVC.
The agreement between the LES and the LDA below the
experimental uncertainty suggests that the PVC is captured
accurately.
Next, the flame brush predicted by LES is compared to
the measurements. Figure 13 shows the mean progress vari-
able fields obtained from normalizing the temperature and FIG. 16. Power density spectrum of the turbulent kinetic energy vs Strouhal
fuel mass fraction measurements 共referred to as experimental number at the location 共x / D = 0.25, r / D = 0.5兲. The line shows a −5 / 3 slope.
FIG. 17. Effect of swirl on the normalized mean axial velocity: radial profiles at x / D = −3.6, −2, 0, and 1.
core. The vortex core is aligned with the centerline in the upstream. Since the flame does not penetrate the premixing
premixing tube but undergoes a strong off-axis rotation close tube 共it starts at x / D ⬃ 0.25兲, the initiation of the PVC is not
to the expansion where vortex breakdown occurs. This mo- affected by thermal expansion. Further characterization of
tion is referred to as PVC and includes a strong precession of the PVC is performed below using POD. Figure 16 also in-
the central recirculation zone. The recirculation zone con- dicates the presence of a smaller peak at St⬃ 0.4.
taining burnt gases is not located instantaneously on the axis
and therefore the flame is not allowed to propagate upstream. Swirl number modulation
It is worth noting that the averaged field offers a poor de-
scription of the PVC and therefore RANS-based techniques The PVC plays an essential role in terms of flame sta-
fail to predict this type of feature. Figure 16 shows a power bility. Hence it is important to study the sensitivity of the
density spectrum sampled in the PVC region. The spectrum PVC and flame stabilization to the value of the swirl number.
exhibits a typical −5 / 3 slope and a well defined peak at the Two additional flow cases with different values of swirl num-
Strouhal number St= f · D / U0 ⬃ 0.6. The corresponding fre- bers have been computed. These cases are referred to as the
quency f ⬃ 880 Hz is in good agreement with the LDA low and the high swirl cases, corresponding to swirl numbers
measurements10 f ⬃ 800 Hz and well in line with PVC fre- of ⬃0.45 and ⬃1.35 共compared to S = 1.05 for the base case兲,
quencies that are reported in the literature, e.g., Ref. 4. respectively. Figures 17 and 18 show the averaged radial
Note that unlike previous studies,7,9,13 we 共as well as and tangential velocity profiles both in the premixing tube
Ref. 10兲 do not witness a strong effect of combustion on the 共x / D ⬍ 0兲 and in the combustor 共x / D ⬎ 0兲. Close to the com-
PVC in this particular configuration. The reason is that the putational inlet 共x / D ⬃ −3.6兲, the differences are small for
PVC initiates upstream from the expansion as pointed out by the three cases. The tangential velocity scales with the swirl
Syred.9 As a consequence, the rms of the velocity fluctua- number. At x / D ⬃ 2, the high swirl and base cases are very
tions 共Figs. 11 and 12兲 shows large values at x / D ⬃ 0 and similar and the profiles match. The axial flow is distributed
FIG. 18. Effect of swirl on the normalized mean tangential velocity: radial profiles at x / D = −3.6, −2, 0, and 1.
FIG. 19. Effect of swirl on the Favre-averaged progress variable: radial profiles at x / D = 0.25, 0.5, 1, and 2.
toward the large radius with a low-velocity region at the axis. at x / D ⬃ 0.5 where the progress variable at the centerline
The low swirl case is still characterized by lower tangential increases with swirl. The low swirl case resembles a short jet
velocity and the axial mass flow is distributed relatively uni- flame due to the swirling motion and the low-velocity region
formly along the radius. Further downstream, vortex break- along the axis. Visualization of the instantaneous flame front
down occurs in the high swirl and base cases as described in and vortex core is presented in Fig. 20 共the low swirl case兲.
the preceding section. The low swirl case exhibits instead a
straight-jet-like behavior with a small low- 共but positive兲 ve-
locity region around the axis. Figure 19 shows the progress
variable profiles along lines in the combustor. The high swirl
and base cases are again similar with a noticeable difference
FIG. 20. Visualizations of the vortex core and flame surface in the low swirl
case: Top: flame surface 共isolevel c = 0.6兲. Bottom: ␭2 isolevel 共Ref. 40兲 FIG. 21. Power density spectra of the turbulent kinetic energy vs Strouhal
共solid gray surface兲 together with the flame surface 共isolevel number at the location 共x / D = 0.25, r / D = 0.5兲. Bottom: high-swirl case. Top:
c = 0.5/transparent surface兲. low-swirl case.
TABLE II. Characteristics of the POD modes.
Normalized Fraction of the turbulent kinetic

Mode number eigenvalue of the energy contained in the N first
N mode N modes
0 0.80048 0.00000
1 0.01470 0.07369
2 0.01454 0.14655
3 0.00576 0.17541
4 0.00529 0.20190
5 0.00404 0.22216
6 0.00377 0.24104
7 0.00367 0.25945
8 0.00346 0.27678
9 0.00340 0.29383
10 0.00319 0.30980
11 0.00317 0.32568
12 0.00305 0.34098
13 0.00294 0.35571
14 0.00286 0.37005
15 0.00265 0.38333
The instantaneous flame front has a plug-type shape; the

front is wrinkled but no recirculation is visible. The vortex
core is aligned with the centerline and enters 2D into the
combustion chamber. The vortex core undergoes a low am-
plitude 共compared to the PVC兲 spiraling motion along the
axis, indicating the presence of coherent motions. Figure 21
shows the power density spectra for the low and high swirl FIG. 22. 共Color online兲 Visualizations of the vortex core as ␭2 共Ref. 40兲
cases. The high swirl case is characterized by a strong well isolevel: Modes 1 共top兲 and 2 共bottom兲.
defined peak at St⬃ 0.55. This case is very similar to the
base case in terms of combustor aerodynamics. The low swirl
number case, on the other hand, does not exhibit any strong
off-axis motion and its power density spectrum is relatively
number itself. Our low swirl case does not show a well de-
flat. Since the sampling location is not on the axis, it might
fined PVC but is relatively close to the critical swirl level as
not capture the small spiraling motion seen in Fig. 20.
the axial velocity at the centerline drops significantly at the
Our results suggest that a well defined off-axis motion
expansion 共see Fig. 17兲.
共PVC兲 appears at some critical swirl number when the cen-
tral recirculation zone is well developed and penetrates into
POD analysis of the PVC
the premixing pipe. Anacleto et al.10 pointed out that the
vortex breakdown occurred at S ⬃ 0.5. They also noted that As pointed out above, averaged quantities are inadequate
the Strouhal number decreases with increasing swirl number to characterize the turbulent flow and the flame under con-
up to S ⬃ 0.9 and decreases when increasing S further. It sideration. In order to understand the PVC, one seeks to iso-
follows from our analysis that the range 0.5⬍ S ⬍ 0.9 corre- late the strong coherent motion from the turbulent fluctua-
sponds to a recirculation zone starting downstream of the tions. One way is to use proper orthogonal decomposition
expansion plane. As S increases, the recirculation moves up- 共POD兲 in order to sort out the different structures. Table II
stream and enters the premixing pipe for S ⬎ 0.9. The PVC shows the eigenvalues of the first 16 modes. One may note
occurs in the expansion region of the premixing pipe and the that mode 0 共averaged field兲 contains about 80% of the total
recirculation is prevented from moving upstream by the con- kinetic energy. Modes 1 and 2 each contain ⬃7% of the
traction. The expansion region of the premixing pipe con- turbulent kinetic energy. The higher 共i.e., nonspiraling兲
strains the annular jet opening angle explaining why an in- modes contain less turbulent kinetic energy. However, it is
crease in swirl number has a small effect on the flame observed that the 15 first modes only contain 38% of the
stabilization. total turbulent kinetic energy. The slow convergence of the
It has to be added that the critical swirl number is not POD series limits the possibilities of efficient reconstructions
universal in any way. The onset of PVC depends on the of the turbulent field from the POD modes. Here we inves-
geometrical details and the radial distribution of the azi- tigate the PVC and therefore we are interested in a relatively
muthal and axial velocity distributions. Thus, the “critical” small fraction of the turbulent kinetic energy only. Actually,
value also depends on the definition used to define the swirl the POD analysis shows that only modes 1 and 2 are relevant
FIG. 23. Vector plots of the modes 1 共top兲 and 2 共bottom兲: longitudinal cut
at a symmetry plane.
to study the PVC. Higher modes have characteristic frequen-

cies 共not presented here兲 that do not match St⬃ 0.6 or any of
its harmonics.
Figure 22 shows a visualization of the vortex core for
mode 1 and 2. Both modes consist of two spiraling vortices.
Modes 1 and 2 are orthogonal by construction such that they
are almost identical but shifted by a rotation around the axis
of ␲ / 2. Figure 23 shows a longitudinal cut representing the
vector plots of modes 1 and 2, where the vortices are local-
ized in the shear layer around the central recirculation zone.
The combination of modes 1 and 2 results in the propagation FIG. 24. Vector plots of the modes 1 共top兲 and 2 共bottom兲: cross section at
of the vortices toward the outlet of the combustor. Figure 24 x / D = 0.35.
shows vector plots of modes 1 and 2 in a cross-sectional
plane, showing that both modes consist of two counter-
rotating vortices. Figure 22 shows a 3D picture of the corre-
sponding vortices but does not give any information about
where x, ␪, and t represent the axial coordinate, the polar
the rotation sign. These findings are well in line with Wang
angle, and the time, respectively. k, m are the axial and an-
and Yang38 POD analysis of an isothermal PVC. Figure 25
gular dimensionless wave numbers while ␻ is the dimension-
shows that time coefficients are presented in Fourier space.
less frequency. From Figs. 23–25, one obtains that for m =
The two modes exhibit a strong and clear peak at the PVC
+ 1, k = 2␲ · 0.96 and ␻ = 2␲ · 0.6. Furthermore, the dimension-
frequency 共St⬃ 0.6兲. The phase angle between the two
less axial phase velocity, cx = ␻ / k ⬃ 0.6, is computed corre-
modes is ␲ / 2, indicating that the modes and time coeffi-
sponding to the speed at which the vortices move toward the
cients are orthogonal both in space and time. In addition, the
outlet in an axial plane 共e.g., Fig. 23兲. The angular phase
combination of modes 1 and 2 characterizes the PVC as a
velocity corresponding to the rotation of the double helix is
double helix 共with counter-rotating branches兲 rotating around
found to be c␪ = ␻ / m ⬃ 1.2␲.
the axis.
Figure 26 shows the central recirculation zone dynamics
Fernandes et al.12 described the PVC as a helical wave,
during one PVC cycle obtained by reconstruction using the
i.e., a harmonic oscillation. It can be characterized by a phase
POD modes 0, 1, and 2. The dynamics consists of a large
function of the form
off-axis rotation so that isolevels characterizing the recircu-
冉冊
lation bubble take a spiral form. Note that the recirculation
x t · U0 zone stabilizes the flame; a similar dynamics for the flame
␸共x, ␪,t兲 = ␸1 k · +m·␪−␻· , 共16兲
D D front is expected.
FIG. 25. Fourier analysis of the POD time coefficients ai共t兲. Top: power
density spectra 共PDS兲 vs St for mode 1. Center: power density spectra vs St
for mode 2. Bottom: temporal phase angle ␾ between modes 1 and 2 vs St.
Study and derivation of a reduced model representing

the flame dynamics
Consider the flow variables vector Q described as FIG. 26. Reconstruction of the PVC using modes 0, 1, and 2 at different
冋册
phase angle during one cycle. Visualization of the central recirculation zone
¯␳ 共isolevel U / U0 = −0.1兲.
Q共x,t兲 = , 共17兲
␶ · ⵜ␳u
N
where the parameter ␶ 共with units of seconds兲 is introduced ␶ · ⵜ␳u共x,t兲 ⬇ 共␶ · ⵜ␳u兲共x,t兲 = 兺 ai共t兲⌽i␶·ⵜ␳u共x兲.
N
共19兲
such that Q will be expressed in units of kg/ m3. Practically, i=1
we choose ␶ so that the terms of the Q vector will be of the
same order of magnitude: here ␶ = D / U0. Figure 27 shows the contribution of each mode to the
We apply the POD on Q in order to get the POD vectors overall density fluctuation. Modes 1 and 2 clearly dominate
for the variables in 共17兲, over the others and have almost the same contributions
⬃8%. Mode 3 contributes some ⬃3% while modes 4 and 5
N have the same contribution of about 2%. The five first modes
¯␳共x,t兲 ⬇ ¯␳N共x,t兲 = 兺 ai共t兲⌽i␳共x兲, 共18兲 together represent only ⬃23% of the total density fluctua-
i=1 tion. Consequently, the higher modes have a global contribu-
FIG. 27. Contribution to the density fluctuation for each of the POD modes.
tion of ⬃77% but this contribution can be attributed to a

very large range of modes and frequencies. As for the POD
analysis of the velocity field, we focus on the modes that are
related to the PVC. Figure 28 shows the first five modes
represented as density isosurfaces. Modes 1 and 2 are rather
similar. Both consist of a pair of central pockets surrounded
by helixes. It should be noted that modes 1 and 2 are or-
thogonal and shifted by a rotation around the axis by ␲ / 2.
Modes 1 and 2 are each antisymmetric as one can see by
rotation of an angle of ␲ around the axis. Consequently,
the corresponding angular dimensionless wave number is
m = + 1 for both modes 1 and 2. One should stress that these
two modes are very similar to the modes extracted from the
POD on the velocity field.
Unlike the previous modes, mode 3 is seemingly axi-
symmetric, it does not belong to a pair and corresponds to a
longitudinal motion of the flame. Modes 4 and 5 are again a
pair of modes and are similar to each other both in term of
contribution and shape. Again it is noted that these modes are
shifted by a rotation around the axis of ␲ / 2. These modes
are invariant if rotated by an angle of ␲ around the axis and
the corresponding angular dimensionless wave number is
m = + 2. Figure 29 shows the time coefficients represented in
Fourier space. Modes 1 and 2 are exhibiting a strong and
clear peak at St⬃ 0.6 which corresponds the PVC frequency.
Mode 3 shows two peaks around St⬃ 0.1 and St⬃ 0.2, indi-
cating a slow longitudinal motion of the flame which is in-
dependent from the PVC. As depicted in Fig. 21, both the
low and medium swirl cases exhibit peaks at relatively low
frequency, suggesting the slow longitudinal motion of the
flame originates in the premixing pipe. In this particular case,
the slow longitudinal motion has low amplitude. It does not
enter in resonance with the modes of the cavity since the first FIG. 28. Isodensity level visualizations of the five first POD modes: light
longitudinal acoustic mode of the cavity is at St⬃ 0.45. surface ⬃ + 0.1共␳u − ␳b兲, dark surface ⬃−0.1共␳u − ␳b兲. Modes 1 to 5 from top
Modes 4 and 5 exhibits the same well defined peaks at to bottom.
St⬃ 1.2 indicating that they correspond to a harmonic of the
mode pair 1 and 2. The angular phase velocity corresponding phase velocity. These modes represent ⬃20% of the total
to the rotation of the modes 1, 2, 4, and 5 is c␪ ⬃ 1.2␲. The density fluctuation. Actually, 80% of these fluctuations are
effect of the PVC or helical wave on the flame is therefore concentrated in modes 1 and 2. It should also be noted that
represented as the sum of four modes with the same angular the first tangential resonant mode of the combustor is at St
TABLE III. Coefficients in Eq. 共20兲.
bi,j / c j i=1 i=2
j=1 −0.11 −3.71

j=2 2.9 −0.27
dynamics. Using the POD modes presented in the previous

paragraphs, the N time coefficients ai are computed enabling
to solve a system of N ODE 关Eq. 共20兲兴. In the present case,
we consider only the two dominant modes, i.e., mode 1 and
2. Equation 共20兲 gives
da1 da2
c1 · ␶ = − a1b1,1 − a2b2,1, c2 · ␶
dt dt
= − a1b1,2 − a2b2,2 . 共21兲
The discriminant of the equation system 共21兲 is

D2 = 共b1,1 / c1 − b2,2 / c2兲2 + 4b1,2 / c2 · b2,1 / c1.. In the case
D2 ⬍ 0, it gives a general solution,
a1共t兲 = A exp − ␣ 冉冊冋冉 t
␶
t
· sin ␻ + B
␶
冊册 , 共22兲
where A and B are constants. The damping coefficient is

␣ = −1 / 2共b1,1 / c1 + b2,2 / c2兲 and the frequency is given by
␻ = 冑−共b1,2 / c2 · b2,1 / c1兲 − 1 / 4共b1,1 / c1 − b2,2 / c2兲2. An arbitrary
time t = 0 is chosen so that a1共t = 0兲 = 0, i.e., B = 0. The ampli-
tude A is obtained by projecting the POD vector ⌽␳1 on the
data set.
The system is further simplified by making use of the
symmetries. For instance, we notice that modes 1 and 2 are
identical but shifted in space by an angle of ␲ / 2. Thus, one
obtains that c1 = c2 and b1,1 = b2,2. In addition, the modes are
found to be antisymmetric by rotation of an angle ␲ leading
FIG. 29. Fourier analysis of the POD time coefficients ai共t兲, PDS vs St.
Modes 1 to 5 from top to bottom. to b1,2 = −b2,1.
From Eq. 共22兲, it is shown that having self-sustained
oscillations 共such as the PVC兲 requires D2 ⬍ 0 and
b1,1 = b2,2 = 0. In addition, the nondimensional oscillation
⬃ 2.9 and therefore it does not enter in resonance with the
PVC. The strongest density fluctuation can be attributed to frequency St is given by St= 冑−共b1,2 / 2␲ · c2 · b2,1 / 2␲ · c1兲
combustor aerodynamics. = 兩b2,1兩 / 2␲ · c1.
We seek to describe the effect of the PVC on the flame We might assess the numerical uncertainty linked to the
as a dynamic system. We insert Eqs. 共18兲 and 共19兲 into the prediction of the ci and bi,j coefficients that are presented in
continuity equation 共1兲, use the orthogonality of the base Table III. First, the coefficients do not exactly satisfy the
vectors, and find the following ODE after projection on the symmetry conditions. The uncertainty is approximately 25%.
POD vector ⌽␳ j: From Table III, the damping coefficient is estimated to be
N
␣ ⬃ 0.19 and the characteristic frequency St⬃ 0.53. The un-
da j ␳ ␳ da j 1 certainty of the characteristic frequency is approximately
关⌽ j ,⌽ j 兴 = c j = − 兺 ai关⌽i␶·ⵜ␳u,⌽␳j 兴
dt dt ␶ i=1 12%. The damping coefficient significantly affects the simu-
lation after about 2 – 3 cycles. The present estimation of the
N
1 uncertainty is relatively high and can be explained by the
= − 兺 aibi,j , 共20兲 limited number of snapshots used for performing the POD.
␶ i=1
However, increasing significantly the number of snapshots
where 关,兴 denotes the scalar product between two vectors. 共ideally one would need statistically independent snapshots兲
Equation 共20兲 describes the dynamical system, i.e., the PVC is not feasible for practical reasons.
1
SUMMARY A. H. Lefevbre, “The role of fuel preparation in low-emission combus-
tion,” J. Eng. Gas Turbines Power 117, 617 共1995兲.
2
S. M. Correa, “Power generation and aeropropulsion gas turbines: From
The present paper presents a numerical study of a swirl- combustion science to combustion technology,” Sym. 共Int.兲 Combust.,
ing flame. The results show that the numerical predictions 关Proc.兴 27, 1793 共1998兲.
only depend weakly on the 共uncertainty in the兲 inflow bound-
3
C. Stone and S. Menon, “Open-loop control of combustion instabilities in
a model gas turbine combustor,” J. Turbomach. 4, 1 共2003兲.
ary conditions and on grid refinement. Comparisons of the 4
L. Selle, G. Lartigue, T. Poinsot, R. Koch, K.-U. Schildmacher, W. Krebs,
present numerical results with experimental data was shown B. Prade, P. Kaufmann, and D. Veynante, “Compressible large eddy simu-
to be limited by uncertainty regarding the averaging proce- lation of turbulent combustion in complex geometry on unstructured
dure as well as by uncertainties related to the experimental meshes,” Combust. Flame 137, 489 共2004兲.
5
C. O. Paschereit, E. Gutmark, and W. Weisenstein, “Control of thermoa-
techniques. Given these limitations, our results indicate coustic instabilities and emissions in an industrial-type gas turbine com-
clearly that large eddy simulation supplemented by a filtered bustor,” Sym. 共Int.兲 Combust., 关Proc.兴 27, 1817 共1998兲.
6
flamelet model captures the flow and the flame well. In par- C. M. Coats, “Coherent structures in combustion,” Prog. Energy Combust.
Sci. 22, 427 共1996兲.
ticular, the behavior of the PVC is captured accurately both 7
C. Schneider, A. Dreizler, and J. Janicka, “Fluid dynamical analysis of
in terms of precessing frequency and radial velocity fluctua- atmospheric reacting and isothermal swirling flows,” Flow, Turbul.
tion amplitude along the centerline. Combust. 74, 103 共2005兲.
8
The numerical results indicate that the PVC is sensitive O. Lucca-Negro and T. O’Doherty, “Vortex breakdown: A review,” Prog.
Energy Combust. Sci. 27, 431 共2001兲.
to the swirl and disappears if the swirl number is reduced 9
N. Syred, “A review of oscillation mechanisms and the role of the pre-
below S ⬃ 0.45. We also found that the complex swirling cessing vortex core 共PVC兲 in swirl combustion systems,” Prog. Energy
flames dynamics are inadequately represented by averaged Combust. Sci. 32, 93 共2006兲.
10
P. M. Anacleto, E. C. Fernandes, M. V. Heitor, and S. I. Shtork, “Swirl
fields that limit the understanding of the flow physics. In flow structure and flame characteristics in a model lean premixed combus-
particular, the PVC induces a nontrivial flame dynamics. As tor,” Combust. Sci. Technol. 175, 1369 共2003兲.
11
a consequence, the averaged flame front does not stabilize C. E. Cala, E. C. Fernandes, M. V. Heitor, and S. I. Shtork, “Coherent
structures in unsteady swirling jet flow,” Exp. Fluids 40, 267 共2006兲.
close to the averaged stagnation point. Instead, the flame 12
E. C. Fernandes, M. V. Heitor, and S. I. Shtork, “An analysis of unsteady
stabilizes further downstream within the averaged central re- highly turbulent swirling flow in a model combustor,” Exp. Fluids 40, 177
circulation zone. It follows that the PVC is initiated upstream 共2006兲.
13
of the flame and does not affect the precessing frequency S. Roux, G. Lartigue, T. Poinsot, U. Meier, and C. Bérat, “Studies of mean
and unsteady flow in a swirled combustor using experiments, acoustic
significantly compared to the nonreacting case. analysis, and large eddy simulations,” Combust. Flame 141, 40 共2005兲.
Performing a POD analysis on the LES data allows iso- 14
T. Poinsot and D. Veynante, Theoretical and Numerical Combustion 共R. T.
lating and identifying large-scale coherent structures. The Edwards, Philadelphia, 2001兲.
15
C. Fureby and F. F. Grinstein, “Large eddy simulation of high-Reynolds-
POD analysis complements averaged data fields with statis-
number free and wall-bounded flows,” J. Comput. Phys. 181, 68 共2002兲.
tical quantities 共i.e., POD modes兲 leading to an accurate rep- 16
F. F. Grinstein and C. Fureby, “Recent progress on MILES for high Rey-
resentation of the flow dynamics. In the present study, a he- nolds number flows,” J. Fluids Eng. 124, 848 共2002兲.
17
lical wave corresponding to the PVC is identified and L. G. Margolin, P. K. Smolarkiewicz, and A. A. Wyszgrodzki, “Implicit
turbulence modeling for high Reynolds number flows,” J. Fluids Eng.
represented by a pair of POD modes or counter-rotating he- 124, 862 共2002兲.
lical vortices. It has also been demonstrated that POD can be 18
G. S. Jiang and C. W. Shu, “Efficient implementation of weighted ENO
applied to the flame front 共density field兲 in order to study the schemes,” J. Comput. Phys. 126, 202 共1996兲.
19
W.-W. Kim, S. Menon, and H. C. Mongia, ”Large-eddy simulation of a
flame dynamics. The effect of the PVC on the flame can be gas turbine combustor flow,” Combust. Sci. Technol. 143, 25 共1999兲.
described using two pairs of modes and an axial mode. The 20
V. Chakravarthy and S. Menon, “Large-eddy simulation of turbulent pre-
second pair is a harmonic of the first pair 共double frequency兲 mixed flames in the flamelet regime,” Combust. Sci. Technol. 162, 1
with the same angular phase velocity corresponding to the 共2000兲.
21
F. F. Grinstein and K. K. Kailasanath, “Three dimensional numerical simu-
helical wave. The POD has also been used to represent the lations of unsteady reactive square jets,” Combust. Flame 100, 2 共1994兲.
22
flame motion and heat release fluctuation using reduced dy- O. Colin, F. Ducros, D. Veynante, and T. Poinsot, “A thickened flame
namical models. Such reduced models can be used advanta- model for large eddy simulation of turbulent premixed combustion,” Phys.
Fluids 12, 1843 共2000兲.
geously to derive control strategies for combustion devices. 23
H. Pitsch and L. Duchamp de Lageneste, “Large-eddy simulation of pre-
mixed turbulent combustion using a level-set approach,” Proc. Combust.
Inst. 29, 2001 共2003兲.
24
ACKNOWLEDGMENTS H. G. Weller, G. Tabor, A. D. Gosman, and C. Fureby, “Application of a
flame wrinkling LES combustion model to a turbulent mixing layer,” Sym.
共Int.兲 Combust., 关Proc.兴 27, 899 共1998兲.
The authors thank Dr. Sergei Shtork 共Lisbon Tech. Uni- 25
C. Fureby, “A computational study of combustion instabilities due to vor-
versity兲 for kindly providing the experimental data and an- tex shedding,” Proc. Combust. Inst. 28, 783 共2000兲.
26
swering many questions about the Lisbon burner. The au- C. Fureby, “A fractal flame-wrinkling large eddy simulation model for
premixed turbulent combustion,” Proc. Combust. Inst. 30, 593 共2005兲.
thors also want to thank Dr. Mirko Salewski and Lisa Prahl 27
C. Duwig and L. Fuchs, “Study of flame stabilization in a swirling com-
for contributing to the paper. The computations were run on bustor using a new flamelet formulation,” Combust. Sci. Technol. 177,
LUNAC and HPC2N computing facilities within the SNIC 1485 共2005兲.
28
P. Domingo, L. Vervisch, S. Payet, and R. Hauguel, “DNS of a premixed
program. This work was supported by the Swedish Energy
turbulent V flame and LES of a ducted flame using a FSD-PDF subgrid
Agency 共STEM兲 and by the Swedish Research Council scale closure with FPI-tabulated chemistry,” Combust. Flame 143, 566
共VR兲. 共2005兲.
29
C. Fureby, “A comparison of flamelet LES models for premixed turbulent Peters and B. Rogg, Lecture Notes in Physics Vol. m-15 共Springer, Berlin,
combustion,” 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, 1993兲, p. 3.
35
NV, AIAA Pap. 2006-155, 2006, pp. 9–12. G. Berkooz, P. Holmes, and J. L. Lumley, “The proper orthogonal decom-
30
C. Duwig, L. Fuchs, P. Griebel, P. Siewert, and E. Boschek, ”Study of a position in the analysis of turbulent flows,” Annu. Rev. Fluid Mech. 25,
confined turbulent jet: Influence of combustion and pressure on the flow,” 539 共1993兲.
AIAA J. 45, 624 共2007兲. 36
T. Smith, J. Moehlis, and P. Holmes, “Low-dimensional modeling of tur-
31
C. Meneveau and T. Poinsot, “Stretching and quenching of flamelets in bulence using the proper orthogonal decomposition: A tutorial,” Nonlinear
premixed turbulent combustion,” Combust. Flame 86, 311 共1991兲. Dyn. 41, 275 共2005兲.
32 37
J. Gullbrand, X. S. Bai, and L. Fuchs, “High-order Cartesian grid method L. Sirovich, “Turbulence and the dynamics of coherent structures. Part 1:
for calculation of incompressible turbulent flows,” Int. J. Numer. Methods Coherent structures,” Q. Appl. Math. 45, 561 共1987兲.
Fluids 36, 687 共2001兲. 38
S. Wang and V. Yang, “Unsteady flow evolution in swirl injectors with
33
D. G. Goodwin, “An open-source, extensible software suite for CVD pro- radial entry. II. External excitations,” Phys. Fluids 17, 045107 共2005兲.
39
cess simulation,” in Proceedings of CVD XVI and EuroCVD Fourteen, M. Couplet, C. Basdevant, and P. Sagaut, “Calibrated reduced-order POD-
edited by M. Allendorf, F. Maury, and F. Teyssandier 共Electrochemical Galerkin system for fluid flow modeling,” J. Comput. Phys. 207, 192
Society, Pennington, NJ, 2003兲, p. 155. 共2005兲.
34 40
N. Peters, “Flame calculation with reduced mechanisms—An outline,” in J. Jeong and F. Hussain, “On the identification of a vortex,” J. Fluid Mech.
Reduced Mechanism for Applications in Combustion Systems, edited by N. 285, 69 共1995兲.

Large Eddy Simulation of Vortex Breakdown/flame Interaction: Related Articles

Uploaded by

Copyright:

Available Formats

Large Eddy Simulation of Vortex Breakdown/flame Interaction: Related Articles

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Large Eddy Simulation of Vortex Breakdown/flame Interaction: Related Articles

Uploaded by

Copyright:

Available Formats

Large eddy simulation of vortex breakdown/flame interaction

Christophe Duwig and Laszlo Fuchs

Citation: Phys. Fluids 19, 075103 (2007); doi: 10.1063/1.2749812

Additional information on Phys. Fluids

Large eddy simulation of vortex breakdown/flame interaction

1070-6631/2007/19共7兲/075103/20/$23.00 19, 075103-1 © 2007 American Institute of Physics

FIG. 1. Production terms ⌸C represented in the density-averaged c space.

the unique coordinate for describing the combustion process. ⳵␳

Figure 1 shows the production term ⌸c mapped in a 共c , a兲

length in the present paper. Downstream of the nozzle, the

TABLE I. Summary of the simulations.

Run Grid Re S Utang / Uax Note

1 Medium 81000 1.05 ⬃0.9 LES ref-Base case

Before presenting the LES results and comparing them

progress variable fields兲 with the LES results. The overall

TABLE II. Characteristics of the POD modes.

Normalized Fraction of the turbulent kinetic

The instantaneous flame front has a plug-type shape; the

to study the PVC. Higher modes have characteristic frequen-

Study and derivation of a reduced model representing

tion of ⬃77% but this contribution can be attributed to a

TABLE III. Coefficients in Eq. 共20兲.

bi,j / c j i=1 i=2

j=1 −0.11 −3.71

dynamics. Using the POD modes presented in the previous

The discriminant of the equation system 共21兲 is

where A and B are constants. The damping coefficient is

You might also like