Aerodynamic characterization of isothermal swirling flows in combustors

Abstract

Swirl flame stabilization is widespread among burners’ manufacturers, but the complex flow patterns are not yet fully understood. The interaction of two confined swirling jets leads to the formation of two recirculation zones being the flame located on the shear layer between the both zones. In such conditions, the lean mixtures can be burned producing low emissions. In the present study, flow structure and turbulent mixing of two isothermal coaxial jets are investigated using Large Eddy Simulation (LES). This is a challenging tool to achieve accuracy but it requests demanding spatial resolution and special treatment of results. By contrasting time-averaged radial profiles with experimental data of a classical benchmark, the model is validated. Results show that LES is able to reproduce the basic features of the flow pattern. Besides, the spectra analysis of instantaneous flow fields provides not only the energy decay but also the most energetic flow structures.

Introduction

This work is part of a research project focused on improving the efficiency of burners. Self et al. [1] evidences the finite reserves of coal and the importance of gas fuel recovery. The final goals are to achieve the stable behavior of lean flames and minimize the emissions [2]. First one means the saving of fuel; the second is involved with policy of environment friendly.

Actual trends to achieve lean flame stability and ultralow emissions are to use a swirl burner. The main advantage is the lower head losses than the bluff bodies and cross-flows. This simple burner considers the use of two coaxial nozzles: one central with fuel and another annular with air being fixed-vanes the precursor of a swirling jet. Both nozzles discharge in a chamber with an expansion ratio lower than three to promote the Outer Recirculation Zone (ORZ) with torus shape due to the sudden expansion. Bearing in mind the flow patterns, there are two kinds of burners [3], high swirl injectors that have an Inner Recirculation Zone (IRZ) with bulb shape due to the backflow in the central part of the chamber, whereas the low swirl injectors have an inner divergence zone (IDZ). Shear region is controlled by both, ORZ near the nozzle exit and the IRZ or IDZ. The flame has a fixed position on the shear region where mixing occurs, having with a shape surrounding the IRZ and IDZ. As combustion is strongly affected by interaction with recirculation zones, [4, 5] first step of the research is to model the mixing of the two confined coaxial jets.

There is a wide range of experimental or numerical studies that characterize the flow pattern for different swirl numbers. Hence, Palm et al. [6] review the effect of swirl no. 0, 0.7 and 1.2. Ranga Dinesh et al. [7] provide information for a swirl no. 0.225 with two setups being the swirling jet the central or the annular. Ranga Dinesh et al. [8] also analyze the swirl no. 0.72 1.01 and 1.44 providing not only the velocity fields but the passive scalar one. In addition, his work [9] is devoted to swirl no. 0.3, 0.54 and 0.8. More challenging approach is that of Maciel et al. [10] that tries to study a swirl number of 0.5 in the limit between low and high swirl injectors. Finally, García-Villalba and Frohlich [11] study the coherent structures for non-confined jets with swirl numbers: 0, 0.4, 0.55, 0.7, 1 and 1.2. Valera-Medina et al. [12] analyze the interaction between the IRZ and the Precessing vortex core (PVC).

The benchmark is that of Roback and Johnson [13] whose NASA report provides large experimental data sets of velocity and Reynolds stress on different lines of the test chamber for the isothermal case. Despite the simple geometrical set-up, the flow pattern shows complex aerodynamic behavior. The sudden expansion of the flow when entering the chamber of expansion ratio two produces the ORZ. If swirl number is large enough to let the flow turn back into the center, the vortex break down phenomenon appears to form an IRZ. The region between both recirculation zones with high shear is where mixture occurs. This is a classical test case that let verify numerical models without reactions. Swirl number around one is produced by eight flat blades of 62° of pitch angle and 25 mm chord.

Computational Fluid Dynamics (CFD) models are feasible tools to design flow devices and predict flow features [14]. Traditional turbulence models are RANS that provide reasonable results using home resources. Large Eddy Simulation (LES) is a state of the art approach able to provide better accuracy than traditional RANS models. A paradigm of the computer science and CFD progress is the fact that LES tools are growing very fast. Both industrial and academic communities are devoting resources to this challenging treatment of the turbulence. Different research groups have access to high-performance computing networks. LES solves the unsteady large scales, while the small ones are modeled. Basically, a spatial filter is applied; the new unknown variables are the dissipation of the smallest scales that should be modeled to solve the closure problem. Domingo et al. [15] present LES models with different sub-grid scales closures to predict turbulent combustion of V flames on sudden expansions.

The aim of this paper is to analyze the interaction of confined non-reacting coaxial swirling jets. LES model of Smagorinsky was selected to simulate the mixing process. The flow Reynolds number was low, around 3,000, where it is based on bulk mean velocity and chamber diameter. Then the regime is transitional turbulent. Swirl number on annular jet is around unity. Quality of the mesh determines the accuracy of the model.

Methodology

In this paper, a classical approach for LES is used. The governing Navier–Stokes equations must be filtered to separate the large scales and small scales of turbulence. The large scales are solved by the discretized equations, while the small scales are modeled through the sub-grid scales (SGS) models.

LES equations are obtained by applying a low-pass filter with width ∆ to the Navier–Stokes equations. The resulting momentum equation is represented by

$$\frac{\mathit{\partial }\mathit{\rho }{\overline{u}}_i}{\mathit{\partial }t}+\frac{\mathit{\partial }\mathit{\rho }{\overline{u}}_j{\overline{u}}_i}{\mathit{\partial }{x}_j}=-\frac{\mathit{\partial }\overline{p}}{\mathit{\partial }{x}_i}+\frac{\mathit{\partial }{\overline{\mathit{\sigma }}}_{ij}}{\mathit{\partial }{x}_j}+\frac{\mathit{\partial }{\overline{\mathit{\tau }}}_{ij}^{\text{sgs}}}{\mathit{\partial }{x}_j}$$

(1)

The unknown sub-grid stress tensor is predicted by the Smagorinsky model, proposed in 1963 and summarized by Pope [16] as a function of the strain tensor and a sub-grid stress viscosity

$${\overline{\mathit{\tau }}}_{ij}^{\text{sgs}}=2{\mathit{\mu }}_{\text{sgs}}{\overline{S}}_{ij}$$

(2)

Being

$${\overline{S}}_{ij}=0.5\left(\frac{\mathit{\partial }{\overline{u}}_i}{\mathit{\partial }{x}_j}+\frac{\mathit{\partial }{\overline{u}}_j}{\mathit{\partial }{x}_i}\right)$$

(3)

$${\mathit{\mu }}_{\text{sgs}}=\mathit{\rho }{\left(C_s\mathit{\Delta }\right)}^2\sqrt{2{\overline{S}}_{ij}{\overline{S}}_{ij}}$$

(4)

Bearing in mind that LES performs well far from walls where turbulent is isotropic, a requisite of the Smagorinsky model is the use of van Driest’s wall-damping function to cancel the sub-grid stress viscosity on the surrounding of the wall when both are multiplied.

$$f_{\mathit{\mu }}=1-exp\left(-\frac{y^{+}}{25}\right)$$

(5)

The spatial resolution on the wall should be ∆y ⁺ = 1, ∆x ⁺ = 100 ∆y ⁺, ∆z ⁺ = 30 ∆y ⁺ such as proposed by [17]. Besides, the time-averaged flow field is evaluated after the simulation had achieved full convergence.

The boundary conditions in both inlets and fluid properties are summarized in Table 1. For LES modeling, as important as the velocities are their fluctuations. Intensity of turbulence in the entrance of both nozzles, central and annular, comes from measurements from the test case.

Table 1 Summary of the inlet conditions for the nozzles

The conservative equations are solved using the Total Variation Diminishing (TVD) algorithm. The TVD algorithm is a finite difference scheme that ensures that the conserved quantities remain monotone and positive. It uses a limiter applied on the high-order discretization scheme, so that it prevents the spurious oscillations from appearing. The larger is the limiter, better accuracy and worst convergence will be achieved.

Finally, PISO (Pressure Implicit with Splitting of Operators) was the algorithm for pressure–velocity coupling since it was the option that provided faster convergence, besides, it does not require under relaxation. Basically, the procedure consists on solving the momentum equations to get provisional velocity fields. Next, the pressure correction equation is solved to get corrected pressure field. Following the velocity correction equation is solved and updated pressure and velocity fields are calculated.

Results and discussion

Parallel processing using the library of subroutines OpenFoam v2.11 was performed using the geometric model and boundary conditions from experimental measurements performed by Roback and Johnson [13], see a review in Table 1. The numerical error of convective transport is controlled using TVD with a limiter that corresponds with good accuracy.

One important concern was the optimum method for spatial decomposition. The criterion was minimizing the number of interface cells since time of information transfer is important.

LES is expected to describe mixing more accurately than traditional RANS approaches. However, the use of a structured mesh with a suitable spatial resolution is a prerequisite of LES [18, 19]. Uniform hexahedral mesh was built and the criterion was enough spatial and time resolution to predict correctly the energy spectrum. The LES models run on parallelization using 128 or 256 nodes for 4 or 2 weeks.

As for postprocess tasks, time-averaged and instantaneous flow fields have been analyzed. Only after full convergence has been achieved, the time average can start. Figure 1 shows the stream lines projected on a plane inside the test chamber over the contours of axial velocity. It is hard to predict flow patterns from instantaneous fields, Fig. 1a, since backflow with negative averaged velocities are often positive in instantaneous fields. Besides, streamlines can be far from being comprehensible. In Fig. 1b, the IRZ has two counter-rotating vortices and is limited by two stagnation points on the axis of the chamber. Its longitudinal length (0.158 m) is determined by the distance between these stagnation points, whereas its transversal length (0.103 m) is the sum of the vortices’ diameters. The centers of the vortices are located at 0.057 m from the chamber’s entrance.

Fig. 1

Velocity field (m/s) and streamlines. Axial thick markers are diameters of the test chamber. a Instantaneous flow. b Averaged flow

A detailed analysis of the influence of the mesh is compulsory. In LES models, the spatial resolution plays the role of a filter of scales smaller than the cell size. Figure 2 shows the comparison with experimental data for two different meshes of 2.5 and 10 million grid points. Details of spatial resolution and y ⁺ values are given in Table 2. Mesh was built using SnappyHexMesh and temporal resolution is ∆t = 10⁻⁵ s

Fig. 2

Radial distribution of averaged axial velocity (left) non-dimensional axial position = 0.04 (right) non-dimensional axial position = 0.42. Axial and radial locations are made non-dimensional with the chamber diameter. Experimental results from Roback and Johnson [13]

Table 2 Comparison of two mesh models

Bearing in mind that both models have the same boundary conditions, it is possible to conclude from Fig. 2a that a fine mesh manages to extract more energy from the main flow generating more turbulence on the macroscopic scales. In Fig. 2b, we can conclude that the mesh of 10 million cells captures better the main characteristics of the experimental profile than the mesh of 2.5 million cells.

The criteria of vortices associated with local low-pressure or high-vorticity make it difficult to visualize 3-D turbulent fields. Alternative criterion to visualize the vortex kernel may be lambda 2 (λ ₂) lower than 0 as shown in Fig. 3. Lambda 2 is defined as the intermediate eigenvalue of the tensor S ² + Ω² (calculated from the strain and rotation part of the gradient velocity).

Fig. 3

Isovolumes of Lambda 2 inside the test chamber representing the core of vortex

The energy spectrum E(κ) is the energy contained in eddies of size l and wavenumber κ, defined as κ = 2π/l. The determination of the energy spectrum requires simultaneous measurements of velocity fluctuations at different points. It is common to assume Taylor’s hypothesis of frozen turbulence [16], which let to measure one velocity component at one point over a certain period of time and then convert the time signal to a spatial signal using x = Ut with U being the time-averaged velocity. It is only valid for u′/U ≪ 1, which is not always the case.

Theoretical energy spectrum as a function of the wavenumber covering the production, inertial and dissipation sub-ranges can be represented by E(κ) = f _L f _η 1.5 ε ^2/3 κ ^−5/3

where the inertial range depends on the wavenumber and the turbulent energy dissipation from dimensional analysis. The production range is governed by f _L , which goes to unity for large eddy scales; whereas the dissipation range is governed by f _η which goes to unity for small scales.

For a turbulence Reynolds number, Re _L  = k ² ε ⁻¹ ν ⁻¹ ≅ 2,400 equivalent to a Taylor-scale Reynolds number of Re _λ  = (20 Re _L /3)^1/2 ≅ 125, the theoretical energy spectrum is represented by the red line in Fig. 4. One probe located on the free turbulence region, (0.12, 0, 0.16) times the test chamber diameter, provide information of the energy decay, line blue. The FFT (Fast Fourier Transformation) analysis of instantaneous register of probes located on the free turbulence region let to confirm that the grid resolution is fine enough to capture the slope −5/3 of the inertial range but the turbulence production associated to the energy containing in the macro-scales is not clear.

Fig. 4

Energy (m²/s²) decay versus wave number (m⁻¹) of a probe at (0.12, 0, 0.16) D, 20,000 samples Δt = 50 μs using Taylor’s hypothesis of frozen turbulence

Proper Orthogonal Decomposition (POD) is a powerful tool for analyzing data since it let identify patterns in data, and express the data in such a way as to highlight their similarities and differences [20]. The POD makes possible to reduce the complete high-dimensional turbulent field to a low-dimensional subspace where different modes are identified and let reconstruct in a simple way a high percentage of the variance of the field. Rotating modes are instabilities which are commonly observed in swirling flows. The cold flow rotating mode is essentially hydrodynamic and corresponds to the well-known PVC (precessing vortex core) observed in many swirled confined flows.

The procedure to decompose the velocity fields into modes is:

• 1. Get as many samples as possible of an instantaneous variable in a cross-section.

• 2. Subtract the mean field from each sample to get the fluctuation of the variable.

• 3. Calculate the covariance matrix.

• 4. Calculate the eigenvectors and eigenvalues of the covariance matrix.

• 5. The eigenvector with the highest eigenvalue is the principle mode that represents the highest energy of the flow.

• 6. Projecting snapshots on the POD modes to form new data set for the principal modes.

Having captured 400 snapshot samples with Δt = 10⁻⁴ s, the eigenvalue spectrum is in Fig. 5 (top-left) where it is possible to establish the number of modes with significant energy contain. Mode 0, Fig. 5 (top-right), shows the streamlines for averaged flow, hence a necessary condition of the number of samples is its symmetrical behavior. Because of the swirling nature of the flow, the modes are coupled by pairs of nearly equal energy, hence mode 1 and 2 are similar, as well as 3 and 4, Fig. 5 (down).

Fig. 5

POD at x = 50 mm (top, left) Energy distribution for each mode (top, right) Mode 0 corresponded to the averaged field (down, left) Modes 1–2 (down, right) Modes 3–4

Conclusions

LES is a very challenging approach to model flows. The accuracy of well-performed models is worthy. The models provide huge amount of information that should be treated on the frequency domain to achieve further details that provide important insight to transport phenomena. LES has been applied to the well-known test case of Roback and Johnson.

Comparison with experimental data evidences that the finest mesh generates more energy decay than the poor resolution mesh. The FFT analysis of probes located on the free turbulence region manages to capture the slope −5/3 of the inertial range. The Proper Orthogonal Decomposition (POD) has been performed on the provisional results. A sensibility analysis to the number of samples and period of sampling is necessary to trust that it is large enough to perform a realistic identification of the most energetic flow structures.

Future works involve other LES models such as implicit LES that could produce better results for transitional flows.