Repositorio Institucional de Documentos

2012 24 Juan Antonio Ramírez Vázquez A computational fluid dynamics investigation of turbulent swirling burners Departamento Director/es Instituto Universitario de Investigación Mixto CIRCE Cortés Gracia, Cristóbal Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Juan Antonio Ramírez Vázquez A COMPUTATIONAL FLUID DYNAMICS INVESTIGATION OF TURBULENT SWIRLING BURNERS Director/es Instituto Universitario de Investigación Mixto CIRCE Cortés Gracia, Cristóbal Tesis Doctoral Autor 2012 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Ph.D. Thesis A COMPUTATIONAL FLUID DYNAMICS INVESTIGATION OF TURBULENT SWIRLING BURNERS By Juan Antonio Ramírez Vázquez May 2012 Advisor: Cristóbal Cortés Gracia, Ph.D. Instituto CIRCE Escuela de Ingeniería y Arquitectura Universidad de Zaragoza ii iii A mis padres y mi esposa iv v A computational fluid dynamics investigation of turbulent swirling burners Juan Antonio Ramírez Vázquez Thesis submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy University of Zaragoza, Spain Abstract This thesis presents detailed numerical calculations of the Unsteady, ReynoldsAveraged Navier-Stokes (URANS) equations to simulate isothermal, single-phase flow in the geometries of realistic swirl burners at large Reynolds numbers. Simulations are run with two different turbulence closures, viz., the standard k−and Reynolds stresses (RSM) models. The numerical method is validated concerning convergence, grid density and far-field influence. Results describe a flow that is in any case periodic or pseudo-periodic, and exhibits quite convincing time-dependent features: bubbleand spiral-type vortex breakdowns and vortex core precession. Some simulations are validated by comparison with corresponding experiments. Good agreement with the experiments has been obtained for mean flow, and frequency peaks of the power spectral density of vi pressure fluctuations. In order to asses the reliability of URANS methods within this context, calculated time-averaged flow and coherent structures are documented via 2D graphs, spectral analysis, 3D isosurfaces and advanced, vortex-related visualization methods and 2D snapshot proper orthogonal decomposition (SPOD). Differences arising from the nature of the turbulence model (k−vs. RSM) are very relevant indeed, given the cost factor involved and the apparent verisimilitude of the predicted flow; they are thoroughly analyzed. CONTENTS xiii 5 URANS of turbulent unconfined swirling burner 83 5.1 Experimental Configuration and Computational Setup ....... 83 5.1.1 Equipment description and experimental details ....... 84 5.1.2 Numerical Modeling, boundary conditions and mesh . . . . 87 5.2 Grid Independence Analysis ...................... 91 5.3 Results and Discussion ......................... 95 5.3.1 Mean Flow Field ........................ 95 5.3.2 Fluctuating velocity and spectra ................ 96 5.3.3 Instantaneous flow and coherent structures ......... 99 6 Summary and conclusions 111 6.1 General conclusions .......................... 111 6.2 Effect of inlet and outlet boundary conditions ............ 112 6.3 Comparison of the k−and Reynolds stress turbulence models . . 113 6.4 Coherent structures eduction ..................... 114 6.5 Comparison of numerical simulation with S-PIV ........... 115 6.6 Perspectives for future work ...................... 116 7 Conclusiones 117 7.1 Conclusiones Generales ........................ 118 7.2 Efecto de las condiciones de contorno de entrada y de salida . . . . 119 7.3 Comparación de los modelos k−y esfuerzos de Reynolds . . . . . 120 7.4 Educción de estructuras coherentes .................. 121 7.5 Comparación de la simulación numérica con el S-PIV ........ 122 7.6 Perspectivas para el trabajo futuro .................. 123 xiv CONTENTS A Basic concepts of turbulent swirling flows 125 A.1 Characteristic Turbulent Time and Length Scales .......... 125 A.2 Governing equations .......................... 130 A.3 The k−Model ............................ 132 A.4 Reynolds Stress Model ......................... 134 List of Figures 4.1 Geometry and computational domain. a) General view. b) Detail of primary air inlet. c) Detail of secondary swirler and throat. d) Origin of coordinates, random monitoring points and axial stations. 49 4.2 Details of the computational mesh: a) General view of the complete geometry. b) Details of primary air tangential inlet. c) Details of the expansion throat. Grid refinement: d) coarse, e) medium and f) fine grids. .................................. 53 4.3 Velocity magnitude at point P1, case 2. (a) Fabricated initial flow field. (b) Static initial conditions. ................... 56 4.4 Fourier transform of the modulus of velocity at point P1, case 4 under different numerical time steps. (a) ∆t= 10−3s. (b) ∆t= 10−4s................................ 57 4.5 Time-averaged axial velocity at different axial locations (Fig. 4.1) for varying grid density. ........................ 59 xv xvi LIST OF FIGURES 4.6 Discretization error in axial velocity at two axial positions. ..... 60 4.7 Fourier transform of the modulus of velocity at point P1 for varying grid density. (a) Medium grid, case 4. (b) Fine grid, case 5. . . . . 61 4.8 Time-averaged velocity field in the windbox and secondary swirler, case 2. (a) Streamlines (b) Contours of velocity modulus. (c) Detail of vector plot. Random control points CP1-CP4 at the inlet of the swirler are shown. ........................... 63 4.9 Time series of the modulus of velocity in the control points represented in Figure 4.8, case 4. ................... 65 4.10 (a) Time series of axial velocity in point P1, cases 4 and 7. (b) Power spectral density of the signal. ................. 66 4.11 Time-averaged axial velocity at different axial positions (Figure 4.1), cases 4 and 7. ........................... 67 4.12 Time-averaged axial velocity at different stations in the plane x-z, cases 2 and 4. .............................. 69 4.13 Streamlines of the time-averaged flow in the x-z plane (a) Case 2 (k−model). (b) Case 4 (RSM). .................. 70 4.14 Time series of the modulus of velocity in the control points P1-P5, (a) Case 2 (k−model). (b) Case 4 (RSM). ............ 72 4.15 Power spectral density of the modulus of velocity at point P1. (a) Case 2 (k−model). (b) Case 4 (RSM). .............. 73 4.16 Snapshots of flow streamlines in the y-z plane, case 4: a) 0, b) π/3, c) 2π/3, d) π, e) 4π/3and f) 5π/3................... 75 LIST OF FIGURES xvii 4.17 Instantaneous profiles of normalized turbulence kinetic energy at different axial positions (Figure 4.1), case 4 (RSM). ......... 77 4.18 Instantaneous flow structures, case 4. (a) Isosurfaces of λ2/up2= 0. (b) Isosurface of uz/up= 0.1. (c) Color plot of λ2/upin the xy plane, z/do=−0.57, and contour of uz= 0. (d) Color plot of λ2/up2in the x-z plane and lines of uz= 0...................... 79 4.19 Instantaneous flow structures, case 2. (a) Isosurfaces of λ2/up2= 0. (b) Isosurface of uz/up= 0.1. (c) Color plot of λ2/upin the xy plane, z/do=−0.57, and contour of uz= 0. (d) Color plot of λ2/up2in the x-z plane and lines of uz= 0...................... 82 5.1 Geometry detail: (a) plenum combustor [95], (b) nozzle and (c) levels and monitoring points. ..................... 85 5.2 Computational domain. ........................ 89 5.3 Mesh details: (a)(d) coarse grid, (b)(e) medium grid and (c)(f) fine grid. ................................... 92 5.4 Spectral analysis for pressure monitored at the point Pac for three time steps. ............................... 94 5.5 Mean axial velocity for the three grids. ............... 95 5.6 Axial and Tangential velocity of the URANS models vs. experimental data. ................................. 97 5.7 Time series of the axial velocity and static pressure [Pa]for monitoring points P1−P5and Pac.................. 98 xviii LIST OF FIGURES 5.8 Power spectral density of the static pressure (Pac) and the axial velocity (P1) for RSM and high swirl case. .............. 99 5.9 Axial vorticity snapshots in the r−zplane for RSM and high swirl case for: (a) 0, (b) π/3, (c) 2π/3, (d) π, (e) 4π/3and (f) 5π/3. The maximum value of velocity is 1000 s−1(white) while the minimum value is −1000 s−1(black). ...................... 100 5.10 Snapshots POD modes contribution: randomly and phase averaging. ................................... 104 5.11 POD reconstruction of azimuthal vorticity for high swirl case: (a) experimental [94] and (b) numerical. ................ 106 5.12 Isosurfaces for the high swirl case and RSM. (a)-(c) are numerical and (d) is an experimental reconstruction. .............. 107 A.1 Energy Spectrum of Homogeneous, Isotropic Turbulence [117]. . . 126 A.2 The Normalized Two-Point Velocity Correlation Function ...... 129 List of Tables 4.1 Summary of the computational cases ................. 54 xix Nomenclature Asnapshot data matrix Aicross-sectional area of the primary air inlet Cnonnegative Hermitian matrix Cµconstant Ddiameter of the annular pipe for secondary air Dsymmetric part of the velocity gradient tensor Eexpansion ratio of diameters (=Do/di) Funitary M×Mmatrix Fe kthe kth component of the Fourier transform of length N/2formed from the even components xxi xxii NOMENCLATURE Fo kthe kth component of the Fourier transform of length N/2formed from the odd components Fs factor of safety GCIij fine grid convergence index using the ith and jth grid G†unitary Nt×Ntmatrix, conjugate transpose of G Hncomplex amplitudes Iturbulence intensity Mnumber of realizations of u(x, t) Nnumber of samples Ncnumber of cells Ntnumber of intervals of time Q2Dsecond invariant of ∇u Rtwo-point spatial correlation tensor Re Reynolds Number Swgeometric swirl number Wcomplex number h·i average operator Chapter 1 Scope, aims and outline of this Thesis This thesis describes the study of Unsteady Reynolds Averaged Navier-Stokes (URANS) schemes for the simulation of isothermal flow in swirl burners. The aim of the present thesis is to demonstrate that URANS can simulate complex swirl flows with good quality compared to experimental data. The thesis proposes also that the use of URANS can provide enough information to understand the complex dynamics of the coherent structures in these flows. In the Introduction (Chapter 2) the complex swirl flow and its role in the generation of coherent structures and vortex breakdown is reviewed, with special emphasis in numerical simulations. The structure and dynamics of different types of vortex breakdowns are described as well as their link with coherent structures. Furthermore, advances in the field of computational fluid dynamics research in turbulent swirl flows are discussed in the context of non-reactive flow. Chapter 3gives an overview of numerical methods used for pre-processing and post-processing data such as grid convergence index (GCI), fast Fourier transform 1 2A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners (FFT), proper orthogonal decomposition (POD) and λ2coherent structures visualization technique. Grid convergence index (GCI) assures a grid-independent solution based on a study of the convergence of average magnitudes, in other words, GCI is a method for estimating the standard uncertainty associated with numerical errors. Fast Fourier transform is used as a tool of spectral analysis of the oscillations in the flow. Proper orthogonal decomposition (POD) and λ2coherent structures visualization technique are two advanced and specialized criteria to study and identify flow features such as precessing vortex core (PVC), coherent structures and detached vortices. Results of URANS schemes for the simulation of isothermal single-phase flow in a burner of pulverized solids are described in Chapter 4. The geometry of swirl generation is completely realistic, comprising a tangential inlet for primary air (fuel transport) and movable guide vanes for secondary air. To model turbulence, two different schemes are used: second order closure by a Reynolds Stresses Model (RSM) and the k−model. In both cases, the model is kept intentionally simple, with no modifications over the standard version. Special concern is taken in assuring a grid-independent solution, by studying the convergence of average magnitudes, and also the characteristics of the oscillations numerically reproduced. The upstream placement of fluid inlets, which is relevant both for the real equipment and the economy of the calculation, turns out to have a pronounced effect on the oscillations, and then on the mean flow; the effect is documented at length. In Chapter 5, a complete analysis of URANS schemes applied to simulate an atmospheric low swirl burner under isothermal conditions is described. This Chapter 1. Scope, aims and outline of this Thesis 3 burner was designed and investigated by Legrand et al. [94,95]. Reynolds Stress Model (RSM) and the k−model schemes are also used to model turbulence. Grid-independent solution is assured following the same procedure as in Chapter 4for the confined burner. Flow features are compared with stereo particle image velocimetry (S-PIV) measurements published by Legrand [94]. Dominant structures are also studied in Chapters 4and 5using the advanced post-processing techniques described in Chapter 3, with special attention to the PVC, coherent structures and detached vortices. The relation between the oscillations and the coherent structures is analyzed. Finally, Chapter 6presents a summary and a general discussion on the major results. New contributions and recommendations for continuing research are also included. 4A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners Chapter 2 Introduction 2.1 Motivation Combustion of primary energy sources is the main source to produce electrical and mechanical energy. According to the statistical report of the International Energy Agency (IEA), in 2008, the 81.3% of primary energy was consumed for this purpose [77]. In 2030, consumption will remain high: 68.1%. For this reason, the effort must to be focus to study, develop and/or modify combustion devices to burn more efficiently. Actually, swirl is one of the main technologies used in new burners. Their application allows to stabilize lean flames in gas turbine burners or to provide stable and complete combustion of solids in industrial power boilers [64]. Unfortunately, swirl generates very complex patterns in the flow in form of vortex breakdown and coherent structures. Külsheimer and Büchner [86] suggest that the formation of these vortical structures are the main mechanism which is 5 6A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners responsible of the instabilities in a combustion process. For this reason, it is very important to study the formation of vortex breakdown and coherent structures [6]. Although the application centers on reacting flows, an adequate modeling of isothermal flow is the first step to the adequate modeling of coupled variable density, escalar transport and reaction. We therefore limit the study to isothermal flow conditions. Swirl is applied also in these conditions, for instance in cyclones and swirl tubes. 2.2 Swirl Phenomena, Vortex Breakdown and Coherent Structures Swirl is employed in diverse technical applications: as a means to effect separation in cyclones, to enhance heat transfer in heat exchangers, to stabilize lean flames in gas turbine burners or to provide stable and complete combustion of solids in industrial and power boilers [64,170,154]. Whenever substantial swirl is present, even at moderate Reynolds numbers, the flow usually develops vortex breakdown and coherent structures. According to Benjamin [10] and Squire [147], vortex breakdown can be conceived as a critical phenomenon of swirling flows, much like an hydraulic jump in a channel. Escudier and Sehnder [38] identified three basic types: axisymmetric, spiral, and double helix; other intermediate forms have been observed depending upon the particular combination of Reynolds and swirl numbers [99,130]. At high Reynolds numbers, it has been reported that the core initiates a kink, followed by a spiral [111]. At a value large enough, bubble and Chapter 2. Introduction 7 spiral structures are suppressed and the flow transforms into a nominally axisymmetric cone of swirling turbulent flow. Sarpkaya [131] considers this the fourth fundamental type. Reviews of the topic have been given by Sarpkaya [130], Hall [66], Leibovich [97], Escudier [33] and Lucca-Negro and O’Doherty [101]. For swirling flows at high Reynolds numbers, different patterns of coherent structures (CS) have been documented by advanced experimental methods (LDA and PIV measurement and visualization), both at isothermal and non-isothermal conditions: precessing vortex core (PVC), inner and outer recirculation zones, and inner and outer secondary helical vortices [47,20,156]. Fick et al. [42] observed how the PVC continuously changes its shape and appearance many times within a single cycle, rotating clockwise and twisting anticlockwise against the direction of the rotating fluid. The basic explanation of a backflow due to the dissipation of the main vortex (e.g., Syred [151]) continuous to hold, but it is clear that it fails to produce a steady flow, causing instead the formation of secondary, non-axisymmetric, unsteady structures [116]. Under non-isothermal conditions (e.g., combustion), flow changes mainly through density; the effect is generally stabilizing, but oscillations may persist, and also acoustic coupling can appear [68,164,124,132,65,163,107,152]. For reviews of PVC and its influence in isothermal and combustion systems see e.g. Syred [151] and Huang and Yang [75]. 8A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners 2.3 Use of Computational Fluid Dynamics in Swirl Flows The current status in computational studies of swirling flows can be summarized as follows. Direct numerical simulation (DNS) has been used with success to reproduce and investigate the complex features of isothermal flows at moderateto-low Reynolds numbers. Some examples are the study of the dinamycs of a swirling jet at Re = 500 of Guohuiet et al. [63] , the fundamental studies of Ruith et al. [125] and Gallaire et al. [49] at Re = 200, the analysis of scalar transport of Freitag et al. [44,45], at Re = 5 000, the mixing properties enhancement of coaxial jets of Balarac et al. [8,7], and the simulation of (modeled) two-phase swirling jets of Siamas et al. [140] at Re = 2 000. For practical conditions, turbulence has to be modeled to some degree. In the extreme of maximum simplicity and minimum cost, closure of the Reynolds Averaged Navier-Stokes (RANS) equations has been practiced since more than two decades for many realizations of swirling flows, leading to apparently realistic steady-state solutions that exhibit backflow, see e. g. the axisymmetric calculation by Mondal et al. [109] and Kriaa et al. [85]. However, in modern times, a variety has emerged called Unsteady RANS (URANS) that presents a clear advantage for unstable swirling flows. Operationally, a URANS simulation simply consists in retaining the time derivatives of the RANS equations while relying as usual on a standard, steady-state turbulence model. Then, for some flows and depending also on the specific model, a periodic or pseudo-periodic, converged solution is Chapter 2. Introduction 9 found, even with steady boundary conditions [142]. At the beginning, the effect was thought to be a purely numerical artifact, and the scheme considered not physically sound for turbulent flows. But soon it became clear that a URANS scheme is simply an economical way of simulating flows that develop discrete natural frequencies, for which it is very superior to steady RANS [31,76]. Accordingly, URANS simulations have been attempted for a variety of swirling flow geometries. Guo et al. [61] used the standard k−formulation to model a suddenly expanded jet, and successfully reproduced the PVC and detached vortices. Jakirlic et al. [79] analyzed three versions of the second-moment closure and two eddy-viscosity models when simulating swirling and rotating pipe flows. They observed that the standard k−model invariably results in unrealistic steadystate, solid-body rotation, which was attributed to poor rendering of streamline curvature effects. More refined turbulence modeling was needed, RSM giving the best performance. Cortés and Gil [26] arrived at a similar conclusion after the simulation of gas flow in a cyclone separator. The fact that the k−model can lead to excessive dissipation in coarse grids has been also signaled as the ultimate reason for steady-state URANS solutions that are not physically sound. Advanced RSM has been used by Ali and Georgios [4] in confined geometries and by Jochmann et al. [81] in expanded jets to demonstrate that URANS schemes can realistically predict time-dependent features of turbulent swirling flows. On the other hand, Large Eddy Simulation (LES) is used more and more for advanced calculations, and swirling flows are not a exception. Some noteworthy examples are the studies of large-scale coherent structures and scalar mixing of Garcia-Villalba et al. [51] and Fröhlich et al. [46], the simulation of aerodynamic 10 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners noise by Flemming et al. [43], the effects of the confinement in a combustor by Lin Lin [100], Stone and Menon [149] and Grinstein and Fureby [60], and . For pulsating flows not dominated by wall phenomena, LES is obviously well adapted, although approximate modeling is normally needed for wall regions in closed geometries. In spite of this, a rigorous application of the technique requires much higher spatial and temporal resolution than any model of the effect of turbulence on the mean flow, and this acquires a special relevance when considering URANS schemes. Likewise, LES also requires very long integration times to build an ensemble-averaged solution [5], whereas just a few periods are usually enough in URANS, since the solution tends to be of a deterministic nature. For these reasons, the prediction that RANS-based simulations, albeit less reliable, will remain useful and competitive versus LES [67] has been mostly accomplished, specially for the kind of flows we consider here (see also e. g. Wegner et al. [161]). In addition, the necessity of simulating complex processes where different physics merge and interact also puts a limit on the computational cost of the models that can be used. LES has now successfully surpassed the barrier of single-phase, isothermal flow calculations. Some significant contributions related to swirling flows are the works of Wegner et al. [161] on non-premixed combustion with fast chemistry, Derksen et al. [30] on interacting gas-particle flow in cyclone separators and Duwig and Fuchs [32] on vortex/flame interaction in premixed combustion via flamelet models. However, as it is obvious, many real-world situations are much more complex. For instance consider a swirlstabilized combustor of pulverized solids. A complete simulation would need to model transport of mass, species, momentum and energy, in the gas phase Chapter 2. Introduction 17 spiral breakdown at very low Reynolds numbers (300 ≤Re ≤750); bubble and spiral breakdown are rather unstable. Spiral vortex breakdown occurs at a Reynolds number range of 750 ≤Re ≤2000 and a swirl ratio range of 2.5≤S≤18. They observed that the whole spiral structure rotates about the centerline in the same sense as the outer flow rotation. At Reynolds number range of 2000 ≤Re ≤3200 three types of VB can be formed: closed bubble breakdown, open bubble breakdown and spiral breakdown. The open bubble breakdown is the Type 1 breakdown and the closed bubble breakdown is the Type 0 of those reported by Faler et al. [39]. In the range of 3200 ≤Re ≤3600 and 2.9≤S≤9.5, the vortex breakdown observed is conical breakdown, documented by Sarpkaya [131]. Vortex breakdown was also characterized by Billant et al. [14] using a swirling water jet. They identified four distinct forms of vortex breakdown: the well documented bubble state, a cone configuration in which the vortex takes the form of an open conical sheet, and two associated asymmetric bubble and asymmetric cone states. In the experiments, the confinement effects can be assumed as negligible since the jet exhausts into a large water tank. For the case of asymmetric bubble, they claim that this structure correspond to the spiral mode of breakdown; differs from the bubble by the precession of the stagnation point around the jet axis. The asymmetric cone is a variation on the cone in the same way as the asymmetric bubble on the bubble. Both asymmetric states are observed at large Reynolds numbers. The basic features that have emerged from the experiments are: (1) abrupt and drastic structural changes occur in a vortex breakdown, (2) axial flow in the core decelerates, sometimes resulting in stagnation and reversal flow, (3) 18 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners the flow is unsteady within the breakdown structure and turbulent downstream, (4) axisymmetric breakdown is characterized by slow oscillations, and (5) the breakdown itself is not a result of instability but a sudden and finite transition from one state to the other as suggested in Harvey [69]. Patte-Rouland et al. [113] studied the recirculation zone in an annular jet using POD analysis of PIV measurements. They observed the interaction of the layers: one between the layer which separates the jet flow and the recirculation zone and the other one which separates the jet flow with the surrounded air. The same interaction between shear layers were analyzed in [73] for four modes of flow structures: bubble, dual rings, vortex breakdown, and vortex shedding. They observed an off-axis saddle point which induces large turbulence intensities. The same large turbulence intensities were observed for a non-swirling flow by Tim et al. [15]. Some other studies have been focused in the analysis of the interior of the bubble vortex. Giannadakis et al. [57] identified low azimuthal vorticity values in the upstream region, close to the swirling nozzle and higher ones in the region between the vortex ring core and the bubble’s aft in addition to low turbulent dynamics inside the bubble. Similar results were observed by Ivanic et al. [78], Liang et al. [98] and Vanierschot et al. [158]. These shear layers have been also identified in combustion applications where the reaction zones are formed at the outer or inner shear layers [59,74,126]. The size of the recirculation zone increases when a secondary coaxial stream is added [119,93]. Lee et al. [93] found that the recirculation zone increases up to above 36%compared with the case of no secondary stream, depending on the pressure ratio of the secondary stream. The swirl direction of the secondary Chapter 2. Introduction 19 stream does not significantly change the pressure distributions along the jet axis. The secondary stream of counter-swirl reduces the size of the recirculation region compared with the size recirculation of co-swirl. Another important feature which alters the recirculation zone is the type of injection topologies, i. e., co-axial and radial, leading to different mixing mechanisms and, hence, altering the recirculation zone [112]. Shtork et al. [137,139] confirmed that the time-averaged flow field characteristics indicate the usual features of swirling jet breakdown with central reverse flow, while phase averaged analysis shows an asymmetrical flow pattern with the vortex core center shifted away from the nozzle axis. On the other hand, Alekseenko et al. [1,2] show the evolution of the recirculation zone size. Moreover, an intense generation of turbulence was observed in the initial region of the jet leading by the vortex breakdown. The intensity of this turbulence (the Reynolds stress huviand third-order moments) was about 5 times higher than in the rest of the flow domain. The forcing increases the total turbulence kinetic energy, but it does not affect the mean flow. They observed that large-scale structures rotating in the opposite direction to the mean flow are responsible for the mixing enhancement. In addition, the results of Coghe et al. [25] show the evidence of different recirculation regions which influence the main combustion features. A toroidal central recirculation region influences reactants mixing and flame stabilization [27,24,12,146,13,84]. The corner recirculating zone induces entrainment of a large amount of hot burned gases into the outflowing reactant mixture while the recirculation zone acts as a bluff-body [168,153,160,159,148,156]. Another important coherent structure present in swirl flows is the precessing vor- 20 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners tex core (PVC). The PVC is a three-dimensional time-dependent CS developed in the exhaust nozzle [106]. Froud et al. [47] and Selle et al. [135] claim that the PVC is caused by the displacement of the center of the vortex. Moreover, they also observed a similar displacement from the axis of symmetry of the recirculation zone which rotates through a region of forward flow. This mechanism provides stabilization and increases the mixing processes when combustion is present [165]. Schneider et al. [133] used LDA to investigate fluid dynamical features caused by combustion process in an atmospheric burner. They observed that downstream axial velocity is maintained while tangential momentum is passed over to radial momentum. A precession of the IRZ was observed leading to distinct frequencies in the PSD when it is compared with reacting case. Fick et al. [42] visualized the PVC and RZ in a combustion process. They observed that the PVC continuously changes its shape and appearance many times within a single cycle in addition to it rotates clockwise and is twisted anticlockwise against the direction of the rotating flow. Schildmacher [132] linked the presence of a PVC with thermo-acoustic instabilities and their interaction with the periodic fluctuations of the velocity and pressure. In addition, they found a phase lag between the different signals. 2.4.3 Numerical Studies The main advantage in the use of numerical simulation for the analysis of complex flows is that it allows to understand flow features by means of a detailed Chapter 2. Introduction 21 explanation of the structure and dynamics for isothermal and non-isothermal flows [145,167,56,170,129,108,55,171]. Spall et al. [144] compared the topological structure of four different types of vortex breakdown (weak helical, double helix, spiral and bubble-types) with the experimental analysis made by Faler and Leibovich [39]. They identified velocity fluctuations which are responsible of the exchange of fluid between the inner zones and the free stream. Moreover, these types of vortex breakdown exhibit an axial stagnation point which indicates the origin of the location of VB [11]. This location depends on a number of parameters: the core Reynolds number, the flow divergence, the swirl velocity ratio, and the strength of the vortex. The breakdown location moves upstream as the core Reynolds number increases, or the initial adverse pressure gradient increases. Some studies have shown that the transition between different types of vortex breakdown are caused by absolute instabilities inside of the recirculation zone [70]. Reynolds-averaged Navier-Stokes (RANS) and unsteady RANS play an important role in the approach of the computation of turbulent flows and heat transfer especially in industrial applications. Some studies assess the performance of different turbulence modes for predicting isothermal flow in complex combustors [82]. RANS models are capable of predict mean and turbulent flow quantities reasonably well except near the wall [82]. Standard k−model is the most frequently used model over the past three decades. It has been also used in the prediction of the precessing vortex core. Guo et al. [61,62] claim that the turbulence model performs extremely well for swirl flow and conclude that the intensive mixing after the breakdown may be 22 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners produced primarily by the small-scale turbulence rather than the large-scale flow entrainment. However, the application of the standard k−model in swirling flows results in a solid-body rotation flow under some circumstances [79]. Despite this, the standard k−is capable to predict the central recirculation zone with enough detail in terms of size, location and strength [109]. Moreover, it is possible to use it to model swirl effect in combustors with good accuracy in the prediction of the temperature regions in flames [55]. The main disadvantage of the standard k−model is that the effects of severe streamline bending due to swirl are unconsidered. Launder [88,89] proposed transport equations to solve the Reynolds stress which allows to account the effects of swirl in a more rigorous manner than standard k−model. Some numerical studies have shown capacity to predict in good agreement with measured data axial and tangential velocities, temperature and turbulent correlations and rms (root mean square) of fluctuating velocities. When rms of the fluctuating velocity components are compared with those obtained by the k−model, the results obtained by RSM are closest to the experimental data [96,162,169]. It indicates that RSM may be capable of predicting the correlations and the mean quantities of swirl flows with enough accuracy to model industrial applications [143]. Chapter 3 Pre and post processing methods 3.1 Solution verification in numerical simulations Finite volume discretization is used to obtain a discrete approximation of Unsteady Reynolds Averaged Navier-Stokes equations for swirl flows. Due to this practice, there is a difference between a quantity simulated and the exact solution of governing equations; this is the called numerical error, δnum. Also, the numerical error have an associated standard uncertainty,unum, which corresponds conceptually to an estimate of the standard deviation, σ, of the parent distribution from which δnum is a single realization. These two quantities, δnum and unum, are used to verify the numerical solution. The objective of verification is to establish numerical accuracy, independent of the physical accuracy that is the subject of validation. In other words, the purpose of verification is to detect inaccuracies in numerical solution and provide 23 24 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners an estimate of error. The procedure to obtain an estimate of error is by means of systematic grid refinement. Grid refinement is no longer necessary once the solution satisfies the error level criterion. The most widely used method to obtain an error estimate is classical Richardson Extrapolation (RE) [120,121]. RE is the most popular form of error estimation since the method requires solutions of the same problem on two meshes. The basic idea of Richardson extrapolation is to obtain an approximation of the leading term in the truncation error from suitably weighted solution on two meshes with different cell size. Although grid doubling (or halving) is often used with RE, it is not required [123], and the ratio of grid spacing may be any real number greater than 1.3. Before to obtain estimate error, it must be ensured that iterative convergence is achieved. Otherwise, the incomplete iteration will pollute the uncertainty estimation. A residual drop of three orders of magnitude in properly normalized residuals for each equation solved over the entire computational domain is a commonly used criterion. For time-dependent simulations, iterative convergence at every time step should be checked. 3.1.1 Grid convergence index The grid convergence index (GCI) is used to provide an error band on the grid convergence of the solution. The GCI is based upon a grid refinement error estimator derived from the theory of generalized Richardson Extrapolation [122]. The objective is to provide a measure of uncertainty of the grid convergence. The GCI is a measure of the percentage between the computed value and the value of Chapter 3. Pre and post processing methods 25 the asymptotic numerical value. It indicates how much the solution would change with a further refinement of the grid. A small value of GCI indicates that the computation is within the asymptotic range. Estimation of discretization error is as follows: 1. Mesh or grid size his defined. h="1 Nc Nc X i=1 (∆Vi)#1/3 (3.1) where ∆Viis the volume of the ith cell, and Ncis the total number of cells used for the computations. 2. Three different set of grids are selected and simulations are run to determine the values of the variable φ. It is desirable that the grid refinement factor, r=hcoarse/hfine, be greater than 1.3. 3. For h1< h2< h3and r21 =h2/h1= 2,r32 =h3/h2= 2, the apparent order, p, is calculated using the expression p=1 ln(r21)|ln |32/21|+q(p)|(3.2) q= ln rp 21 −s rp 32 −s(3.3) 26 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners s= 1 ·sign(32/21)(3.4) where 32 =φ3−φ2,21 =φ2−φ1. Negative values of 32/21 <0are an indication of oscillatory convergence. It should be noted that if either 32 =φ3−φ2or 21 =φ2−φ1is "very close" to zero, the above procedure does not work. This might be an indication of oscillatory convergence or, in rare situations, it may indicate that the "exact" solution has been attained. 4. The extrapolated values are calculated from φ21 ext = (rp 21φ1−φ2)/(rp 21 −1) (3.5) 5. The approximate relative error, extrapolated relative error and the fine grid convergence index, along with the apparent order p, are calculated. e21 a= φ1−φ2 φ1(3.6) e21 ext = φ12 ext −φ1 φ12 ext (3.7) GCI21 fine =Fs ·e21 a rp 21 −1(3.8) Chapter 3. Pre and post processing methods 33 tion (flow visualization, conditional methods, variable integration time average, pattern recognition analysis), no a priori is needed for the eduction scheme. CS are defined in an objective and unique manner as the flow realization that possesses the largest projection onto the flow field. Secondly, the POD yields an optimal set of basis functions in the sense that no other decomposition of the same order captures an equivalent amount of kinetic energy. Up to now, POD is only presented as a data analysis method that takes as input an ensemble of data, obtained from physical experiments of from detailed numerical simulations, and extracts basis functions optimal in terms of the representativeness of the data. POD can also be used as an efficient procedure to compute low-dimensional dynamical models of the CS. Due to the optimality of convergence in terms of kinetic energy of the POD functions, only a small number of POD modes are necessary to represent the dynamical evolution of the flow correctly. 3.3.2 POD approximation method Suppose a vector-valued function u(x, t)over some domain of interest Ωs. It can be approximate as a finite sum in the separated-variable form: u(x, t)≃ K X k=1 a(k)(t)φ(k)(x)(3.17) xcan be viewed as a spatial coordinate and tas a temporal coordinate. A classic way to solve this approximation problem is to use for the basis 34 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners functions φk(x), functions given a priori, for example Fourier series, Legendre polynomials of Chebyshev polynomials. An alternative approach could be to determine the functions φk(x)that are naturally intrinsic for the approximation of the function u(x, t). An additional difficulty is that a different sequence of time functions a(k)(t) corresponds to each choice of basic functions φk(x). So, given φk(x), the coefficients a(k)(t)can be determined as follows. Suppose we have chosen orthonormal basis functions, i. e., ZΩs φ(k1)(x)φ(k2)(x)d(x) = δk1k2(3.18) where δk1k2=   0for k16=k2 1for k1=k2 (3.19) is the Kronecker delta symbol, then: a(k)(t) = ZΩs u(x, t)φ(k)(x)dx(3.20) Therefore for orthonormal basis functions, a(k)(t)depends only on φ(k)(x)and not on the other φ. So far selecting the function φ(k)(x), it would be useful to use the orthonormality condition. Now consider experimental or numerical data at Ntdifferent instants of time, Mrealizations of u(x, t)at Mdifferent locations x1,x2,. . .,xM. The Chapter 3. Pre and post processing methods 35 approximation problem of Eq. 3.17 is then equivalent to finding the orthonormal functions φ(k)(x)K k=1 with K≤Ntthat solve: min Nt X i=1 u(x, ti)− K X k=1 u(x, ti), φ(k)(x) 2 2 (3.21) where k·k2define the norm associated with the usual L2inner product (., .). The practical method of solving the minimization problem of Eq. 3.21 is to arrange the data set U={u(x, ti), . . . , u(x, tNt)}in an M×Ntmatrix Acalled the snapshot data matrix A=         u(x1, t1)u(x1, t2)··· u(x1, tNt) u(x2, t1)u(x2, t2)··· u(x2, tNt) . . .. . .. . .. . . u(xM, t1)u(xM, t2)··· u(xM, tNt)         , A ∈RM×Nt(3.22) Each column A:,i ∈RMof the snapshot data matrix represents a single snapshot u(x, ti)of the input ensemble U. It is noted that, if the snapshot data are assumed to be linearly independent, the snapshot matrix has full column rank. 3.3.3 POD applied to turbulent flows Based on previous and basic analysis. Let {u(X),X= (x, tn)∈D=R3×R+} denote the set of snapshots obtained at Ntdifferent time steps tnover a spatial 36 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners domain of interest Ωs. These snapshots could be numerical solutions of velocity fields, vorticity fields, etc. taken at different time steps. The underlaying problem is to extract from this ensemble of random vector fields a coherent structure. Defining a coherent structure as the deterministic function which is the best correlated on average with the realizations u(X). In other words, a function Φthat has the largest mean square projection onto the observations |(u,Φ)|2is looked for. There are two methods to find Φ. One where the average h·i is temporal and is evaluated as an ensemble average, based on the assumptions of stationarity and ergodicity. The variable Xis assimilated to the space variable x= (x, y, z)defined over the domain Ωs. This is the direct method or classical POD. The other method is the so-called snapshot POD method which is the exact symmetry of the classical POD. The average operator h·i is evaluated as a space average over the domain Ωs. The snapshots are taken at different times. The time step is usually constant but this is not necessary. The only requirement is that the snapshots are linearly independent. This method is efficient when the spatial domain is higher than the number of observations. Each method has particular characteristics but it is relatively easy to choose the pertinent method for each practical configuration. For example, on the one hand, data obtained by numerical simulations can be highly resolved in space and time but due to cost considerations only a very short time sample is simulated. Conversely, a good spatial resolution can be obtained by particle image velocimetry, but associated with a poor temporal resolution. On the other hand, experimental approaches such as hot-wire anemometry Chapter 3. Pre and post processing methods 37 or laser Doppler anemometry provide a well-defined time description but with limited spatial resolution. These measurement techniques enabled long time histories and moderate spatial resolution. Data issued form an experimental approach will generally be treated using the classical method and data issued from numerical simulations by the snapshots method. An exception is the case of data sets obtained from particle image velocimetry. 3.3.4 Snapshot POD To derive the discrete eigenvalue problem corresponding to the snapshot POD, it is assumed that Φhas a special form in terms of the original data Φ(x) = Nt X k=1 a(tk)u(x, tk)(3.23) where the coefficients a(tk),k= 1, . . . , Ntare to be determined solving the Fredholm integral eigenvalue problem ZΩs R(x, x0)Φ(x0)dx0=λΦ(x)(3.24) The two-point spatial correlation tensor R(x, x0)is estimated under stationary and ergodicity assumptions as: 38 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners R(x, x0) = 1 TZT u(x, t)⊗u∗(x0, t)dt =1 NT Nt X i=1 u(x, ti)⊗u∗(x0, ti)(3.25) Substituting this expression of Rand the decomposition of Φ(Eq. 3.23) into Equation 3.24 Nt X i=1 "Nt X k=1 1 NtZΩs u(x0, tk)·u∗(x0, ti)dx0a(tk)#×u(x, ti) =λ Nt X k=1 a(tk)u(x, tk) (3.26) and concluding that a sufficient condition for the coefficients a(tk)to be a solution of Equation 3.24 is to verify that Nt X k=1 1 Nt [u(x0, tk)·u∗(x0, ti)] a(tk) = λa(ti), i= 1, . . . , Nt. (3.27) This can be rewritten as the eigenvalue problem CV=λV(3.28) Chapter 3. Pre and post processing methods 39 where C=1 NtZΩs u(x0, tk)·u∗(x0, ti)dx(3.29) and V= [a(t1), a(t2), . . . , a(tNt)]T(3.30) Since Cis a nonnegative Hermitian matrix, it has a complete set of orthogonal eigenvectors V(1) =a(1)(t1), a(1)(t2), . . . , a(1)(tNt)T, V(2) =a(2)(t1), a(1)(t2), . . . , a(2)(tNt)T,..., V(Nt)=a(Nt)(t1), a(Nt)(t2), . . . , a(Nt)(tNt)T (3.31) along with a set of eigenvalues λ(1) ≥λ(2) ≥. . . ≥λ(Nt)≥0. Then, the temporal eigenfunctions Vican be normalized by requiring that 1 Nt (Vn,Vm) = 1 Nt Nt X k=1 a(n)(tk)a(m)∗(tk) =λ(n)δnm (3.32) Then, the POD eigenfunctions Φ(n)(x)are estimated as 40 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners Φ(n)(x) = 1 Ntλ(n) Nt X k=1 a(n)(tk)u(x, tk)(3.33) 3.4 Vortex definition Vortices are a special existence form of fluid motion with origin in the rotation of fluid elements. It can be possible to recognize the existence of vortices first by their intuitive streamline patterns, which are however not Galilean invariant and cannot be used to define a vortex. A natural invariant approach could be based on the vorticity, from which one can extract vorticity lines and vorticity magnitude. Saffman and Baker [128] defined a vortex as a connected fluid region with high concentration of vorticity compared with its surrounding. In other words, a vortex is a vorticity tube surrounded by irrotational flow. But the vortex boundary becomes fuzzy in viscous flow without sharp boundary. There are a some cases where vortices are axisymmetric of which the outer boundary, of radius ro, has the maximum value of the circumferential velocity. However, this criterion cannot be generalized to more complex and nonaxisymmetric vortices. A simple alternative to the vortex definition would be identifying the fluid region with |ω| ≥ |ω0|, where |ω0|is a threshold magnitude. But this criterion is also inadequate because the choice of |ω0|is subjective, and the side boundary of a vorticity tube may significantly differ from an isovorticity surface. A natural basis for developing possible rational criteria is the symmetricantisymmetric decomposition of the velocity gradient tensor, ∇u,∇u=D+ Ω, which suggests that a vortex may be defined as a flow region where the Chapter 3. Pre and post processing methods 41 vorticity (symmetric tensor Omega) prevails over the strain rate (symmetric tensor Omega). This requires the calculation of the invariants of the velocity gradient tensor through its representative matrix, say Aλ; which in cylindrical coordinates reads Aλ=     u,rv,rw,r (u,θ−v)/r (u,θ+u)/r w,θ/r u,zv,zw,z     ,(3.34) where subscript ,ris the partial derivate with respect to radius, ,θis partial derivate with respect to angular coordinate and ,zis partial derivate with respect to axial coordinate. The first criterion along this line was proposed by Weiss for two dimensional incompressible flow (u, v) based on the eigenvalues σof ∇u, of which the characteristic equation is σ2+Q2D= 0 (3.35) where Q2D= u,xv,x u,yv,y=1 2(kΩk2−kDk2) = 1 4ω2−1 2kDk2(3.36) is the second invariant of ∇u(and also the negative of the discriminant ∆2D; the first invariant is tr(∇u)=0). Here, it is considered kSk ≡ tr(S·ST)1/2for any tensor S. When Q2D>0 at a point, the flow is called elliptic and we have 42 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners purely imaginary eigenvalues ±˙ıσ˙ı=±√Q2D; for the case of Q2D<0 the flow is called hyperbolic. Thus, a vortex is defined as a connected fluid region with Q2D=−∆2D=σ2 i>0(3.37) known as the Weiss criterion. In particular, if instead of Cartesian we use cylindrical r-phi coordinates (Eq. 3.34 in 2D), for a 2D vortex, there is Q2D=1 4(1 r ∂ ∂r(rv)2 −r∂ ∂r v r2)=v r ∂v ∂r (3.38) and σ2 i>0precisely defines the vortex as a fluid within r=r0where v= max, in consistency with the common concept of vortex core. Controversy on defining a vortex appears once three-dimensional flow is considered. The characteristic equation for the eigenvalues of ∇uis σ3+Qσ −R= 0 (3.39) where Q≡ −1 2ui,juj,i =1 2(kΩk2−kDk2) =1 21 2ω2−kDk2=σ1σ2+σ1σ3+σ2σ3, (3.40) Chapter 4. URANS of a turbulent confined swirling burner 49 Figure 4.1: Geometry and computational domain. a) General view. b) Detail of primary air inlet. c) Detail of secondary swirler and throat. d) Origin of coordinates, random monitoring points and axial stations. 50 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners 4.1.2 Physical models, boundary conditions and numerical methods The incompressible unsteady Reynolds-averaged Navier-Stokes (URANS) equations are solved adopting two different closures for turbulence, the standard k−model [92] and a Reynolds stress model with linear pressure-strain term [58,48,88]. Wall reflection terms are included in the RSM, in order to consider pressure blocking and redistribution of normal stresses [28]. This was deemed necessary to adequately model wall-dominated regions, such as the spaces between the vanes of the secondary swirler. Standard wall functions [92] are used for near-wall modeling. The first grid point is located at a maximum distance of 40 < y+<60 for all meshes. Inflow conditions are idealized by assuming uniform velocity profiles at the inlet sections, up= 15.11 m/s and us= 1.17 m/s for primary and secondary air, respectively. Turbulence intensity is estimated from fully developed flow correlations as I= 0.16Re−1/8where Re is based on hydraulic diameter. Turbulence kinetic energy is then k= 3/2(uI)2and its dissipation rate for the k−model =C3/4 µk2/3l−1, with Cµ= 0.085. Integral length scales are taken as lp= 0.05 m and ls= 0.1m. For the RSM, inlet Reynolds stresses are determined under the assumption of isotropic turbulence, i.e., u02= 2k/3,hu0v0i= 0. Even for swirling flows, a high value of the Reynolds number normally permits to anticipate a small sensitivity to conditions at the outlet boundary (see e.g., Xia et al. [166]). In order to impose adequate conditions in the present simulations, we performed a brief far-field study. The usual expedient of zero axial velocity Chapter 4. URANS of a turbulent confined swirling burner 51 gradients was imposed on three different geometries, each having a different chamber length of 7do,14doand 21do, which corresponds to the real open end in isothermal air flow conditions (14do), and two imaginary open chambers, one shorter and one longer. Values at the exit plane did change slightly between the first and the second geometry, but they did not change appreciably between the second and the third. Accordingly, we extend the computational domain to a length of 14dodownstream of the throat, and use the usual outflow conditions there. Simulations are performed with the CFD solver FLUENT 6.3.26. We employ second order central differences for convection and diffusion terms and an implicit second-order scheme for the time derivatives. The SIMPLE algorithm is used as the pressure-velocity coupling method, taking care in adopting adequate time steps for the unsteady calculation [9]. 4.1.3 Computational mesh The computational domain is divided in three zones, corresponding to the volumes occupied by primary air, secondary air and combustion chamber. For reasons of convenience, the mesh is unstructured in the tangential inlet of primary air and in the windbox. The remaining zones (annular ducts, secondary air inlet zones and swirler, throat and combustion chamber) are represented by structured meshes. For the study of grid independence, we used three progressively finer meshes. Following the recommendations of Celik et al. [21], a grid refinement factor greater than 1.3was chosen to minimize truncation errors, specifically r= 2. 52 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners Thus, the medium and fine grids have respectively 8and 64 times as many points as the coarse grid. Figures 4.2(d), (e) and (f) illustrate the geometric relationship. The number of computational nodes for the full geometry is, successively, 84,600, 676,800 and 5,414,400. 4.1.4 Computational cases and numerical performance Seven different computational cases were ran, as summarized in Table 4.1. Cases 1,2and 3use the k−model in the coarse, medium and fine grids. They served to the grid independence study, based on recommended techniques from the literature, that were applied to the time-averaged flow (see 4.2). Additionally, variation of time-dependent features between the medium and fine grids were studied using RSM, cases 4and 5respectively. The main investigation then considered the medium grid for the comparison of URANS solutions under different turbulence models, cases 2and 4. Aside from the "complete" geometry shown in Figure 4.2 (a), a simplified geometry that omits the windbox was considered relevant for the study. The reasons and the outcome are explained in Sect. 4. Simulations were repeated accordingly, cases 6and 7. The numerical computation doesn’t converge in neither of these cases to a statistically stationary flow, but to oscillatory fields of velocity and pressure. When using the RSM, a representative time per global iteration is 5.19 s; it does not drop very much if the k−model is used instead: 3.03 s. However, the second order closure typically takes 28 iterations to converge in each time step, whereas the twoequation model only requires 6. In conclusion, times needed are approximately in Chapter 4. URANS of a turbulent confined swirling burner 53 d) e) f) Figure 4.2: Details of the computational mesh: a) General view of the complete geometry. b) Details of primary air tangential inlet. c) Details of the expansion throat. Grid refinement: d) coarse, e) medium and f) fine grids. 54 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners Table 4.1: Summary of the computational cases Run Grid Geometry Turbulence model 1 Coarse Complete k− 2 Medium Complete k− 3 Fine Complete k− 4 Medium Complete RSM 5 Fine Complete RSM 6 Medium Without windbox k− 7 Medium Without windbox RSM a proportion of six-fold. Time values refer to a Beowulf-type cluster using 6CPUs of 2,200 MHz. 4.2 Convergence and grid independence We report in this section the studies undertaken to assure the quality of the numerical predictions. Basically, the computational procedure behaves well and converges within the preset tolerance to a solution that is independent of the initial conditions and the size of spatial and temporal increments. However, our case is very special in these respects, since we have an oscillatory flow. As a consequence, not only instantaneous or averaged values of flow magnitudes must be studied, but also their frequency content, for the range of frequencies that can be considered a genuine prediction of the URANS technique. This can be stated from a slightly different perspective. A numerical simulation with steady-state conditions is able to reproduce flow periodicities only by means of an initial amplification of numerical errors, which triggers the natural instability embedded in the flow Chapter 4. URANS of a turbulent confined swirling burner 55 model, Ruith et al. [125]. Therefore, the converged solution must necessarily satisfy an additional condition, viz., that the oscillations behave independently of arbitrary initial values and grid size. Otherwise, it would be clear that a purely numerical artifact has been obtained, and no claim of representation of the real flow physics could be made. 4.2.1 Convergence of the oscillatory flow Independence with respect to the initial condition is assured by repeating the same computational cases starting from widely different flow patterns. Two extreme possibilities are a stationary medium (nil velocity and pressure everywhere), and an artificial flow field purposely fabricated to shorten the transient period and speed up the convergence to the final solution. The first may represent a first approach to the real transient experienced by the physical system (although this is not of interest in this study). There are many possibilities for the second. For instance, aside from simply guessed fields, one can use solutions from simpler models in the same grid, or interpolated values from the solution in a coarser one. Figure 4.3 shows as an example results from case 2of Table 1. The magnitude represented is the modulus of velocity in the control point P1 of Figure 4.1(d). The fabricated flow field comes from a steady (RANS) solution obtained by using a first order scheme for spatial discretization. The attainment of a steady-periodic regime after an initial transient is clearly observed in both cases. We used in both a variable under-relaxation parameter, starting at a small value and increasing it gradually until the solution began to settle down to a stable oscillation. The 56 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners (a) (b) Figure 4.3: Velocity magnitude at point P1, case 2. (a) Fabricated initial flow field. (b) Static initial conditions. change is apparent at 2.5sin the fabricated flow case, Figure 4.3(a); the final oscillatory flow is attained approximately at this point. Logically, the transient is longer for an initially stationary fluid, Figure 4.3(b). As it is apparent in the figure, and can be demonstrated numerically, the two final solutions share the same average and extreme values, waveforms and frequency. (Phase needs not to be equal, for obvious reasons.) A second consideration is that, since the frequency content predictable by the calculation is unknown beforehand, it is convenient to estimate the effect of the numerical time step ∆ton the oscillations. We repeated selected test cases for three values, ∆t= 10−3,5×104and 10−4s. Figure 4.4 shows the results for case 4, monitoring point P1, for the extreme values of t. Here we represent the power spectral density (PSD) of the velocity modulus vs. the Strouhal number based on Chapter 4. URANS of a turbulent confined swirling burner 57 (a) (b) (c) (d) Figure 4.4: Fourier transform of the modulus of velocity at point P1, case 4 under different numerical time steps. (a) ∆t= 10−3s. (b) ∆t= 10−4s. 58 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners axial momentum and outlet diameter of the throat. The URANS calculation produces distinct low frequency peaks, that can be related to coherent structures formed in the flow, as we will explain later. Tuned with these oscillations, also high frequency peaks are obtained, but with considerably less power (note the logarithmic ordinate), and entering into the inertial subrange of “background” turbulence, where fluctuations are supposedly modeled, so that they shouldn’t convey any fundamental information. The decay exhibited at the right of the graph is just due to a Blackman window function used in the spectral analysis. The spectrum is logically much less noisy the lower the time step, but strength and location of the low frequency peaks only suffer minor variations. A value of ∆t= 10−3s was used accordingly for the rest of the study. 4.2.2 Grid independence To estimate the error of the numerical simulation, we use well-established procedures from the CFD literature, that involve repeating the calculation in three progressively finer grids, cases 1-3 of Table 1. Since the procedures only apply to steady-state situations, we firstly proceed as if they were also valid for our time-averaged flow. After finding out that the outcome is acceptable, we address separately the effect of grid size on the oscillations. Figure 4.5 shows axial velocity profiles at six axial stations computed in the three grids. Differences are indeed small, and very similar flow patterns are predicted. In order to quantify the error, we follow the systematic procedure recommended by Celik et al. [21], or grid convergence index (GCI) method, that Chapter 4. URANS of a turbulent confined swirling burner 65 Figure 4.9: Time series of the modulus of velocity in the control points represented in Figure 4.8, case 4. 66 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners (a) (b) Figure 4.10: (a) Time series of axial velocity in point P1, cases 4 and 7. (b) Power spectral density of the signal. with different extreme values; the oscillation is clearly stronger for the complete geometry. In fact, Figure 4.10(a) clearly suggests different average values in point P1. Figure 4.11 shows the time-averaged radial profile of axial velocity at different axial stations. Differences are indeed noticeable. Magnitude of backflow (positive values) is stronger for case 4, and the inner recirculation zone (IRC) is predicted longer than in case 7. Also the mixing seems stronger in the complete geometry: The positive peaks that signal the presence of primary and secondary fluid injection dissipate upstream in case 4 compared to case 7, and the profile evolves more smoothly. As a conclusion, a different flow is predicted if a simplified representation of Chapter 4. URANS of a turbulent confined swirling burner 67 Figure 4.11: Time-averaged axial velocity at different axial positions (Figure 4.1), cases 4 and 7. the entrance is attempted. It is clear in this case that the simulated oscillation propagates itself upstream to the settling chamber, in a manner that makes the ensemble pulsate differently, with important consequences in the flow patterns. The flow is more oscillatory, which intensifies dissipation and the effect of swirl. Accordingly, cases 2 and 4 are retained for the main study, and 6 and 7 discarded. However, the general implication for the significance of URANS simulations of pulsating flows is not positive. If an adequate prediction of the flow needs such an accurate representation of the entrance, generalization of the models is very difficult. For instance, the small windbox of our pilot combustor is easily amenable to numerical rendition, but this won’t be the case of a large, industrial windbox serving a set of individual burners. And vice versa, predictions of ideal cases or laboratory-scale apparatuses could not be extrapolated easily to full-size 68 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners equipment. 4.4 Results and discussion This section examines the oscillating flow solutions obtained by URANS simulations that use the k−and Reynolds stresses models of turbulence, cases 2 and 4 of Table 1, respectively. When representing time-series and computing Fourier transforms, averages and other statistics, we discard the initial transient period, retaining only the periodic part. We have organized the presentation in two subsections: time-averaged flow and time-series and spectra. 4.4.1 Time-averaged flow Figure 4.12 shows profiles of time-averaged axial velocity. We can observe a long inner recirculation zone (positive velocities) and two characteristic peaks that signal primary and secondary air injections. They last up to z/do=−0.093 and are dissipated by z/do=−0.373. Differences between both predictions are negligible from z/do=−0.747 downstream and slight upstream. However, the latter are significant: the simulation with the RSM predicts higher reverse velocities in the near field of the throat. Streamlines of the time-averaged flow in an axial plane are shown in Figure 4.12, with intervals chosen so as to distinguish vortex and reverse flow zones clearly. Flow inside and close to the throat exhibits what appears to be two toroidal vortexes. The larger is centered approximately at the throat outlet and Chapter 4. URANS of a turbulent confined swirling burner 69 Figure 4.12: Time-averaged axial velocity at different stations in the plane x-z, cases 2 and 4. penetrates into the chamber, conforming what is classically described as an Inner Recirculation Zone (IRC), see e.g. Figure 1.2 of Syred [151]. There is addition a small torus attached just at the lip of the annular duct, obviously a product of flow detachment. Downstream at the sides we observe the secondary corner vortex originated by the jet. Finally, a second recirculation torus develop that completely fills the chamber. As for the differences due to the turbulence model, both flows seem to be similar in their basic features and differ only in the details. With the RSM, larger vortices are predicted in the near flow. Their centers are located at z/do=−0.3, z/do=−2.2for case 4 and at z/do=−0.5,z/do=−1.9for case 2. The upper vortex in this latter case (k−) is of a very reduced size, almost indistinguishable in the plot. In accordance with Figure 4.12(d)-(f), the time-averaged flow apparently 70 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners (a) (b) Figure 4.13: Streamlines of the time-averaged flow in the x-z plane (a) Case 2 (k−model). (b) Case 4 (RSM). converges down-stream to an even closer similitude; for instance, the width of recirculation zone at z/do=−3is roughly the same, about 1.9dofor case 4 and 2dofor case 2. This picture is however misleading, because the values and intervals of streamlines are not the same. Rotation in the plane seems stronger and more extended for the k−prediction (case 2), but it is so because we are employing much shorter intervals of stream function values in order to detect it. Thus, what we see with the black areas of Figure 4.13(a) are only weak vortexes, that consequently encompass large areas. The flow predicted with the RSM (case 4), Chapter 4. URANS of a turbulent confined swirling burner 71 Figure 4.13(b) actually rotates much more vigorously, clearly exhibiting vortexes with concentrated rings and of a richer inner structure. Finally, it should be noted that the prediction of case 2 is almost axisymmetric, which is in agreement with boundary conditions, barring the surely minor detail of the primary swirler. In contrast, case 4 predicts an asymmetric flow, specially at the throat, in a manner that cannot be related easily to geometry. 4.4.2 Time-series and spectra of flow magnitudes Figure 4.14 shows time series of the modulus of velocity in the random control points P1-P5 of Figure 4.1(d), as predicted for cases 2 (k−model) and 4 (RSM). Velocity seems to be oscillatory everywhere, with an amplitude surpassing 50% of the average in points located upstream and close to the wall (P1). The amplitude drops as the flow evolves downstream, points P5, P3 and P2. Also the oscillation is less pronounced as we move farther from the wall, point P2 vs. point P3. Velocity in P4 is only residually oscillatory, especially for case 2. All these details clearly suggest that the peripheral coherent structures of Figure 4.13 move downstream with the flow, with a central IRZ of ascending flow that is more stable. Another noteworthy feature is the fundamental difference of the predictions when the RSM is used (case 4) instead of the k−model (case 2). In both cases, non-sinusoidal, frequency-rich waveforms obtain. However, the signal is purely periodic for case 2, whereas it presents long time variations for case 4. These variations are much slower than the significant frequency content of the apparent wave-form; it seems as if the second-order closure was successful in 72 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners (a) (b) Figure 4.14: Time series of the modulus of velocity in the control points P1-P5, (a) Case 2 (k−model). (b) Case 4 (RSM). modeling a partly stochastic behavior. Obviously this is related to the symmetry of the averaged flow, Figure 4.13: a deterministic prediction will always result in a symmetric flow, whereas a stochastic component seems to realize itself also spatially in non-symmetric flow patterns. The power spectral density of the signal in point P1 is shown in Figure 4.15 for both predictions. It has been calculated by FFT, and a Blackman window function has been applied to attenuate high-frequency components and clarify the logarithmic plot. A clear dominant, low frequency is predicted in both Chapter 4. URANS of a turbulent confined swirling burner 73 (a) (b) Figure 4.15: Power spectral density of the modulus of velocity at point P1. (a) Case 2 (k−model). (b) Case 4 (RSM). calculations, although the values differ notably: St = 0.647 for case 2 and St = 1.0 for case 4, which corresponds to frequencies of 15.11 and 23.49 Hz, respectively. The peak is somewhat more definite for case 4, but the strength is similar. The spectra also contain secondary harmonics whose frequencies are approximately in a relation of 2 with the dominant and between themselves. They are clearly more intense and definite for case 4. The first three have been indicated in the figure. As noted above, this behavior continues with a much diminished strength up to the high-frequency part of the spectrum. The same dominant frequencies were obtained for pressure, velocity components and velocity magnitude in all monitored points. 74 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners 4.4.3 Three-dimensional, time-dependent flow and coherent structures To begin the description of the time-dependent flow, we show in Figure 4.16 six snapshots of flow streamlines in the x-z plane, as predicted in case 4 (RSM), during a complete cycle of oscillation at the dominant frequency (St = 1, after Figure 4.15(b)). We observe how the corner vortex seems static and mostly symmetric, though deformed and shaped by the central jet, that is completely dynamical. Apparently, the central ignitor space acts as the rear of a bluff-body from whose surface an unsteady pattern of alternate vortexes detach. But obviously this is an axisymmetric geometry, so that what we actually see is a spiral vortex (SV) permanently attached to the rim of the primary air outlet, that rotates a whole turn. Thus, a vortex break-down of the spiral type, time-dependent and asymmetric, is predicted for our geometry and Reynolds number. The symmetry (or quasi-) and the throat vortexes seen in Figure 4.13(b) are only the timeaveraged footprint of this coherent structure, with a strong mark where the spiral vortex is attached and a composite down-stream. The geometric complication of the burner itself seems to have a minor role. The secondary fluid stream, of roughly double momentum and half the swirl, doesn’t originate additional CS, but limits itself to be unsteadily throttled and mixed with the primary stream. Figure 4.17 illustrates another way of looking at this question. Although modeled by RSM, the second moments of the turbulence (for instance, its turbulence kinetic energy) should indicate regions of high shear. The figure shows radial profiles of k/u2 pthat correspond to the same instant as Chapter 4. URANS of a turbulent confined swirling burner 81 significant, viz.: geometry of throat and swirler, single air entrance and infinite expansion ratio (unconfined flow). This may perfectly explain the differences in flow structure, notably, the absence of a SV in the near-field. In any case, resemblance is noteworthy. On the other hand, we think that the presence of two different air streams has only a limited effect on flow features, due to the high rate of mixing that attains the flow. Finally, Figure 4.19 repeats the same plots for the prediction of case 2. All the features mentioned above are present, but with significant differences. The prediction with the k−model of turbulence appears obviously less rich, and at once more basic, than its counterpart predicted by using the RSM. However, perhaps important features are lost. In particular, all vortical CS are weaker and dissipate earlier, and this affects specially the two secondary vortexes. Accordingly, the IRZ is only distorted by the main spiral vortex inside the throat; it stabilizes and gets symmetric right after entering the chamber. Thus, the flow there can be considered almost steady, contrary to the more elaborated prediction of case 4. The significance of these differences for a combustion system cannot be assessed at this point and must be investigated, i.e., to what extent the loss of intensity of the main vortex changes the near-flow dynamics and whether the practical absence of downstream secondary vortices and PVC is of importance. 82 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners Figure 4.19: Instantaneous flow structures, case 2. (a) Isosurfaces of λ2/up2= 0. (b) Isosurface of uz/up= 0.1. (c) Color plot of λ2/upin the xy plane, z/do=−0.57, and contour of uz= 0. (d) Color plot of λ2/up2in the x-z plane and lines of uz= 0. Chapter 5 URANS of turbulent unconfined swirling burner 5.1 Experimental Configuration and Computational Setup In this chapter, the capacity of prediction of URANS schemes for the simulation of isothermal flow in an atmospheric low swirl burner is studied. A 50 kW atmospheric low swirl burner designed by Legrand et al. [95] is considered. The main characteristic of this burner is that produces a weak recirculation zone which stabilize the flame achieving ultra-low emissions [22,23]. There are some studies with similar burners which use LES schemes for the flow at isothermal and reactive conditions [72,71,52], all of them show advantages for the identification of coherent structures, but we pretend to analyze the performance of the URANS 83 84 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners schemes under isothermal conditions. It is important to mention that numerical results have been compared with the non-reactive experimental measurements made by Legrand et al. [95]. Their experiments are detailed in the following sections. The burner is completely realistic and some S-PIV measurements have been published [95]. As in Chapter 4, the flow has been simulated with two different, basic, turbulence models, namely second order closure by a Reynolds Stresses Model (RSM) and the k−model. In order to assure good numerical results, special emphasis was put in grid-independent solution and the characteristics of the oscillations varying the time step. Flow features are studied and compared with S-PIV measurements. Snapshot POD is used to extract all the information of numerical simulations and link it to physical measurements. Then, the coherent structures are visualized and described with the help of experimental information. 5.1.1 Equipment description and experimental details In the experiments, velocity and vorticity fields for reactive and non-reactive flows were obtained with S-PIV and then CS were reconstructed via proper orthogonal decomposition (POD). Authors used two CCD cameras and a 532 nm wavelength, 400 mJ Quantel (Twin Brilliant B) double Pulsed Nd:YAG laser for the illumination. The image size was approximately 80 ×80 mm2 (2000 ×2000 pixels) in size. For the isothermal case, they used propylene glycol particles, of which 90% have a diameter less than 2µm with a pick in probability density at 1µm. With this size of particle they captured maximum frequencies of 1.2kHz. Moreover, a Chapter 5. URANS of turbulent unconfined swirling burner 85 21 Figures. Figure 1.Geometry detail: (a) plenum combustor [27], (b) nozzle and (c) levels and monitoring points. (a) (b) (c) Figure 5.1: Geometry detail: (a) plenum combustor [95], (b) nozzle and (c) levels and monitoring points. 86 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners 1/2inch B&K condenser microphone was placed at 5.7Dofrom the axisymmetric axis to acquire pressure signals and via fast Fourier transform calculate acoustic spectral power densities. In this chapter, we attempt to simulate two non-reactive cases based on their experiments: low swirl case (SL= 0.58) and high swirl case (SH= 0.64). Both under the same mass flow rate, Re = 12000 [95]. Reynolds number is based on reference length (Do) and ˙maccounts for the total mass flow rate measured at the burner exit. From continuity, a characteristic velocity is calculated by uo= 4 ˙m/ρoπD2 owith the objective of dimensionless velocity profiles. Some arbitrary control levels and monitoring points were considered to compare numerical simulations with experiments, Fig. 5.1c. Moreover, an additional point, located in the same position as the microphone used in experiments to measure acoustics, is also considered. The control level at z/Do= 0.1is used to analyze the jet flow at the nozzle outlet, where the accelerated flow interacts with surrounding static fluid. The other three control levels (z/Do= 0.5, 1and 2) are located in the zone where coherent structures develop. More monitoring points have been sampled in numerical simulations, but only four are shown. The first two monitoring points are located inside the nozzle due to the fact that in some CFD studies with similar geometries, (inner pipe retracted), it has been demonstrated that CS begin to form inside the nozzle [52]. The main reason to choose P3 and P4 positions is the interest in capturing velocity fluctuations of the main structure and also capturing the inner recirculation zone, respectively. The point used to monitor static pressure Pac„ that is not plotted in the figure has cylindrical coordinates of (0,5.76Do,0). It is important to note that the shape Chapter 5. URANS of turbulent unconfined swirling burner 87 shown in Fig. 5.1b corresponds to a part of the zr-plane where the nozzle and pipe walls are depicted; the inner pipe is retracted a distance Dofrom the nozzle exit. 5.1.2 Numerical Modeling, boundary conditions and mesh Before beginning with the study of the performance of URANS models applied to an unconfined swirl burner, some tests were carried out to make sure that results are independent of the position and type of boundary conditions. The same procedure was followed for a confined swirl burner numerical simulation in Chapter 4. A comparison between two geometries was performed. The first case considering that the computational domain contains the cylindrical plenum and tangential pipes (which feeds the secondary flow), the annular jet and the inner axial pipe (which feeds the primary flow) and the nozzle and surroundings (where both flows merge). The second case considers a reduced geometry where the cylindrical plenum and tangential pipes are not included but approximate boundary conditions used instead. The analysis was as follows. First, numerical simulations for the complete geometry using RSM were executed based on the mass flow rate of low swirl case (SL= 0.58 and Re = 12000). Control levels (z/Do = 0.1,0.5,1and 2) were used to obtain mean velocity values. Also, some monitoring points were located along the annular pipe region, all of them before the nozzle zone. The purpose of these additional monitoring points was to identify velocity fluctuations. It was found that instabilities produced by vortex breakdown propagated always downstream 88 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners and never upstream. In other words, fluctuations inside the annular pipe were inexistent. This result makes it possible to construct a reduced case in which the annular pipe is cut in a region where the flow is fully developed. As inlet boundary condition, a previously simulated velocity profile is imposed in order to ensure that swirl generation is maintained for the reduced geometry that lacks cylindrical plenum and tangential pipes, Fig. 5.2. Values of axial, tangential and angular velocities in control levels (z/Do= 0.1,0.5,1and 2) for this reduced geometry were compared with those of the complete case. It was found that the velocity profiles were similar in all control levels and the frequencies of the fluctuations in the vortex breakdown zone were equal in both cases. It is arguably that the results differ from those in Chapter 4where the conclusion was that the swirler must be simulated. However, this only indicates that local disturbances in confined swirl burners are propagated both upstream and downstream so that it is very important to choose a correct position of inlet boundary conditions, and, if the disturbances propagate until the swirler zone, the swirler needs to be part of the domain. In contrast, for unconfined swirl burners, it seems that local disturbances propagate only in the downstream direction, so that it is possible to make numerical simulations without considering the swirler in the domain. Based on this, the computational domain includes only the annular jet, the inner axial pipe and the nozzle, Fig. 5.2. Inflow boundary conditions are imposed assuming uniform and steady velocity profiles at the inlet sections. For primary air, only the axial velocity component is imposed, while for secondary air, tangential Chapter 5. URANS of turbulent unconfined swirling burner 89 Figure 5.2: Computational domain. 90 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners and axial velocities components are considered. For the high swirl case (SL= 0.64 and Re = 12000), numerical simulations begin with arbitrary values of velocity components. When the axial and tangential velocity profiles are obtained at level z/Do= 0.1, they are compared with experimental measurements at the same level. Then, the values of velocity components are adjusted until swirl number and velocity profiles match the experimental values. The methodology has been also validated by García-Villalba et al. [50,51,52] who demonstrated good accuracy in averaged and statistic results without the costly representation of the plenum. With the problem of inlet boundary condition solved, simulations were performed with the Computational Fluid Dynamics solver FLUENT 6.3. The standard k−model [92] and a Reynolds stress model with linear pressure-strain term [48,58,88] are adopted to solve the incompressible unsteady Reynoldsaveraged Navier-Stokes (URANS) equations. Second-order central differences for convection and diffusion terms and an implicit second-order scheme for the time derivatives are employed. The SIMPLE algorithm is used as the pressure-velocity coupling method. Standard wall functions are used for near-wall modeling [92]. The first grid point is located at a distance of y+= 50 for all meshes. A far-field study is performed using a constant value of 1 bar for pressure as boundary condition at the exit. This value corresponds to the ambient pressure. Three different geometries with different lengths (9Do,18Doand 27Do) are compared following the procedure of Chapter 4. Accordingly, the computational domain is extended to a length of 18Do. With respect to the lateral boundary, a constant value of 1 bar for pressure as boundary condition is used. Also, a Chapter 5. URANS of turbulent unconfined swirling burner 97 (a) axial velocity for low swirl case (b) tangential velocity for low swirl case (c) axial velocity for high swirl case (d) tangential velocity for high swirl case Figure 5.6: Axial and Tangential velocity of the URANS models vs. experimental data. 98 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners Figure 5.7: Time series of the axial velocity and static pressure [Pa]for monitoring points P1−P5and Pac. Chapter 5. URANS of turbulent unconfined swirling burner 99 Figure 5.8: Power spectral density of the static pressure (Pac) and the axial velocity (P1) for RSM and high swirl case. velocity P1, the dominant frequency is 750, but a minor peak at a frequency of 478 Hz can be observed too. In experimental results at the same Reynolds and swirl numbers, the peak is located at a frequency of 500 Hz [94]. 5.3.3 Instantaneous flow and coherent structures Before starting with coherent structures analysis, snapshots of the simulations of axial vorticity are presented in Figure 5.9 for the case of RSM and high swirl number. This figure shows the snapshots during a complete cycle of oscillation based on the stronger frequency peak of P1in Figure 5.8. White color indicates positive values of axial vorticity while dark colors indicate negative values. Thickline indicates zero axial velocity. The purpose of drawing the zero axial velocity 100 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners Figure 5.9: Axial vorticity snapshots in the r−zplane for RSM and high swirl case for: (a) 0, (b) π/3, (c) 2π/3, (d) π, (e) 4π/3and (f) 5π/3. The maximum value of velocity is 1000 s−1(white) while the minimum value is −1000 s−1(black). line is to identify the so-called Inner Recirculation Zone (IRZ) of the swirl flow. This zone is represented in Figure 5.9 by the region delimited by the thick-line in the coordinates −0.5< r/Do<0.5. The shape of the IRZ is distorted and asymmetric. The shape and asymmetry of this coherent structure is related to the vortices observed in Figure 5.9, as explained below. In the right side of the Figure 5.9a there are two white color vortices coupled with three black ones. (Black and white indicate the sense of rotation). Chapter 5. URANS of turbulent unconfined swirling burner 101 Some white vortices clearly touch the zero axial velocity line as well as the complementary three black vortices on the opposite side of the graph. The term Inner Vortex (IV) is applied to those structures which touch the zero axial velocity line, whereas Outer Vortex (OV) is considered to comprise those structures coupled with the IV in the zone of positive axial velocity. It can be noted that the vortices evolve similarly as a 2D von Kármán Vortex Street with the IRZ acting as a bluff body. All the vortices propagate downstream and vanish during the first half of the period, Figure 5.9a-d. During the second half of the period, Figures 5.9d-f and back to 5.9a again, new vortices arise with weak intensity compared with those structures in the first half of the period. However, at this instance, it is impossible to know the type of vortices present in the flow. For that reason, it is necessary an advanced analysis of numerical data. The snapshot-POD method is employed with the complementary method of eduction of vortical coherent structures developed by Jeong and Hussain [80]. Legrand et al. [95] collected enough experimental data, (1,000 S-PIV statistically independent snapshots), and used the proper orthogonal decomposition method as a technique to identify coherent structures. The main advantage of using the POD method is that it consists of a linear procedure which processes input data, based on the Fredholm integral eigenvalue problem, that creates an orthogonal basis of a, in this case, non-linear phenomenon. In other words, POD takes an ensemble of the data and extracts basis functions; these functions are optimal in terms of the representativeness of the data. There are two different POD existent approaches: Classical POD and the snapshot POD. In classical POD, the variable Xis assimilated to the space x=(x,y,z) 102 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners defined over the domain ΩSand, based on the assumptions of stationarity and ergodicity, it is evaluated as an ensemble average or, in other words, as a temporal average. On the other hand, the variable Xin snapshot POD is assimilated to the time tand it is evaluated as a space average over the domain ΩS. It is necessary to evaluate which one of the two POD approaches is more adequate for S-PIV and numerical simulations. S-PIV provides a good spatial resolution, but associated with a poor temporal resolution. Numerical simulations are highly resolved in space and time but only a very short time sample can be simulated. Therefore and based on the previous description, snapshot POD is the chosen method for treating experimental and numerical data. Sirovich [141] introduced snapshot POD method as a way to capture a good picture of the large scale behavior in turbulent flows. The principal advantage of this technique is the fact that the autocovariance matrix can be approximated by a summation of snapshots instead of solving a n×neigenvalue problem which is very time-consuming. In their experiments, Legrand et al. [95] used this method to reconstruct a helicoidal vortex by means of azimuthal vorticity 2Dfields. This structure is, presumably, the responsible of the pressure oscillation mechanisms which leads to the acoustic peak measured by the microphone. But, a deep analysis is necessary to understand the interactions between the coherent structures and measurements. Recently, numerical simulation has emerged as a complementary tool for S-PIV measurements due to its capacity to resolve the system with a considerable quality in space and time. The performance of URANS models is evaluated by using these advanced postprocessing tools. For this, a case has been selected with high swirl, Chapter 5. URANS of turbulent unconfined swirling burner 103 RSM and a time step of δt = 2x10−5. A window of 80 ×80 mm2, positioned in an equivalent place as in the experimental S-PIV measurements, has been used to obtain numerical data. The statistical representation was computed from 1,000 snapshot (N). These 1,000 snapshots were collected with the same sampling rate as in experiment, 0.5seconds. The procedure followed is the following. First, a matrix M=ω0ω1··· ωNis defined, where each column corresponds to a snapshot of the azimuthal vorticity (ω). Then, the autocovariance matrix is calculated, C=MTM, and the eigenvalue problem solved, CAi=λiAi. Once eigenvalues are calculated, they are sorted in descending order, λ0> λ1>··· > λN, as well as the eigenvectors (each column of Ai). Finally, the normalized POD modes were found as: φi=PN n=1 Ai nωn PN n=1 Ai nωn , i = 0,··· , N (5.1) Each POD mode, φi, is the representation of the energy contained in the flow. In this sense, the most energetic realization is the mode 0 where the energy contained is about 72.03%; modes 1 and 2 contain 8.22% and 4.65%, respectively. Only three modes have been analyzed due to the fact that this quantity is enough to produce vortices [95]. It is possible to add more modes to improve smaller details of the vortex field, but some studies have demonstrated that the main structures are unchanged [115,114]. Figure 5.10 the contribution of POD modes to each. The random nature of the sampled numerical data is clearly observed. Because of this, the data need to be shorted. In this context, a phase averaging using POD coefficients can be obtained 104 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners Figure 5.10: Snapshots POD modes contribution: randomly and phase averaging. following the same procedure used for experimental data by et al. [95]. Data of azimuthal vorticity (ωθ) can be decomposed as a sum of the timeindependent mean flow, a quasi-periodic fluctuating component and a random fluctuating component: ωθ=ωθ+fωθ+ωθ0; the term ωθ+fωθis the phase averaged vorticity. Writing it as a Fourier expansion, with the use of only the three first POD modes and retaining the fundamental frequency, hωθi∼ = ωθ+P∞ n=1 (Bncos(ϕ) + Cnsin(ϕ)), Legrand [94] approximates the autocovariance matrix, C, in an alternative way and calculate the eigenvalues as: Chapter 5. URANS of turbulent unconfined swirling burner 105 λ0 i≈1 √N λ1 i≈ ±r2 Ncos 2πi N(5.2) λ2 i≈ ±r2 Nsin 2πi N Using this procedure it can be observed that the eigenvalues depend only on the number of snapshots; also, it can be noted that the eigenvalues 1 and 2 are shifted by a quarter of period. Due to this, some authors have related these two POD modes to the convection of the vortices. Figure 5.10 shows the POD mode contribution that has been alternatively calculated, fit lines, and the POD mode contribution based on the autocovariance matrix C=MTMthat has been alternatively calculated and the fit lines. It can be observed a sinusoidal evolution of the random collected data and their respective strong dispersion when the data have been phase averaged. With this, the vorticity fields can be re-evaluated using any phase angle. Figure 5.11 shows the phase averaged POD reconstruction for (a) experimental and (b) numerical data with a value of phase angle ϕ= 0◦; both graphics share the same vorticity scale. It can be observed that numerical simulations are able to predict the same number of cores. Furthermore, structures captured by SPIV of experimental data are of the same size in comparison with vortices in numerical simulations. Or, in other words, RSM seems to be an acceptable way of predicting the strength of the coherent structures. In addition, it is observed 106 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners (a) (b) Figure 5.11: POD reconstruction of azimuthal vorticity for high swirl case: (a) experimental [94] and (b) numerical. that the inner vortices are surrounded by counter-rotating vortices, and they are clearly divided by shear layers. It can be noted that experimental and numerical POD reconstruction lacks the IRZ observed in Figure 5.9. It indicates that axial vorticity dominates over azimuthal vorticity for the IRZ. Since the axial vorticity is capable to predict both the inner vortices and IRZ, we use it for the analysis of the coherent structures. Since the simulation has been validated with experimental measurements, it is possible to perform a 3D analysis of the coherent structures. Three methods are usually employed to visualize vortices: isosurfaces of axial velocity, isosurfaces of axial vorticity and isosurfaces calculated with the λ2technique (method developed in [80]). Figure 5.12 shows the coherent structures obtained following these Chapter 6. Summary and conclusions 113 The results showed that coherent structures were maintained in the two burners. But, only in the atmospheric swirl burner case, the mean of the axial velocity and the frequency peaks obtained by means of the estimation of power spectral density of the static pressure monitors gave the same value independently of the type of inlet condition. This means that the instabilities produced by vortex breakdown in atmospheric burners propagate downstream and never downstream. In the case of the confined swirl burner, the instabilities propagates in all directions. It is worth mentioning that the position of the outlet boundary condition was also analyzed. Although the influence is less compared with the inlet condition, an optimal distance was found where the effect of the outlet boundary condition was minimal. 6.3 Comparison of the k−and Reynolds stress turbulence models The simulated flow has been studied by several methods of post-processing, including advanced techniques for eduction of vortical coherent structures, also with the objective of comparing the performance of the two turbulence closures. Both give realistic predictions; they describe a complex, and quasi-periodic flow with spiral and helical vortexes that conform a precessing vortex core and a pulsating inner recirculation zone. This behaviour is captured by the monitors. For the confined swirl burner, when these signals are processed by FFT, the frequency 114 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners peak predicted by the k−model differs in a 35% compared with the prediction of the RSM, where the frequency peaks are stronger. Results are quite convincing and compare reasonably with the experimental literature. Particularly realistic is the prediction via RSM, for which the flow lacks exact symmetry and periodicity, and exhibits more stronger and persistent vertical motions. In contrast, albeit with essentially the same features, the k−model leads to a flow of more schematic nature. For the atmospheric swirl burner, the mean velocity showed that both models predict values close to the experimental measurements. For the low swirl case, the axial velocity adjusts better for the experimental data with the RSM compared with the k−model while for the tangential velocity the k−model gives a better approximation. For the high swirl case, the axial and tangential velocities adjust better with the k−model. With respect to the time series, they showed a similar behavior present in the confined burner, namely, the signals are periodic. For this case, the frequency peaks were validated with experimental data where the difference of the numerical simulation and experimental data differ only in 4.4%. 6.4 Coherent structures eduction Actually, there are different techniques to identify the coherent structures present in a turbulent flow. For this thesis, two common techniques (velocity and vorticity isosurfaces) and one advanced visualization technique (λ2) were considered. In all techniques, it was possible identify the recirculation zones, spiral vortex, helical Chapter 6. Summary and conclusions 115 vortex and some spiral and toroidal structures that form around the main flow. However, due to the complexity of coherent structures and the interaction each other, only the λ2technique was able to capture with enough clarity all the structures including those weak structures. Velocity and vorticity isosurfaces are incapable to detect these weak structures. More importantly, the coherent structures predicted by the two turbulence models have great differences when they are analyzed by λ2technique. The coherent structures given by k−model are weak compared with the RSM which indicates that they dissipate quickly. By the other hand, the RSM capture many details of the coherent structures presents in the flow, i. e., the inner recirculation zone (IRZ) which is distorted by the precessing vortex core (PVC). 6.5 Comparison of numerical simulation with S-PIV Exhaustive numerical simulations of the swirling turbulent flow experimentally characterized by S-PIV by Legrand et al. [95] have been presented, under the premise of detecting pulsating phenomena, vortex breakdown and coherent structures, and investigating the usefulness of economic URANS methods in this context. Comparing the numerical results obtained by RSM with the experimental S-PIV measurements, it can be concluded that the turbulence model is able to predict the same number of cores. The numerical flow resembles very much the experimental one, and duplicates its features. Average flow, pulsating frequency and vortex intensities are adequately predicted. The utility of URANS computations for simulating complex 116 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners flows is evident, at least as a complementary tool together with experimental measurement, since they produce all the necessary information to construct 3D coherent structures within an acceptable approximation. 6.6 Perspectives for future work Aside from rigorous experimental validation in the model combustor, future work will be centered on elucidating the essential balance of an URANS simulation: economy of the turbulence model vs. quality of the prediction. To this end, the URANS scheme for the flow will be coupled with a second phenomena, such as turbulent dispersion of particles, density variation through temperature gradients or species mixture, or density variation through simple models of partially premixed gas combustion. From a point of view of post-processing techniques, classical spectral analysis generally uses fast Fourier transform (FFT) which is a common tool in practical applications. However, fast Fourier transform fails in the processing of short data sequences. In other words, fast Fourier transform is not capable to capture very low frequency peaks. But now, various other approaches are available: autoregressive (AR), moving average (MA) and autoregressive moving average (ARMA) methods. These methods, developed for radar, sonar of geophysical applications, etc. are not well known for fluid mechanics and combustion applications. The advantage of the use of modern spectral methods is that very low frequency peaks could be obtained, adding frequency resolution and statistical stability to the current information. Chapter 7 Conclusiones Estudios numéricos de flujo inestable, isotérmico, monofásico, turbulento en un quemador piloto de combustible pulverizado de giro inducido y en un quemador atmosférico de bajo giro inducido han sido realizados usando los modelos k−estándar y el modelo de esfuerzos de Reynolds. Se ha tenido especial cuidado en asegurar precisión numérica, independencia de malla y una adecuada representación del campo lejano. También, la exacta representación de los dispositivos de entrada revela que es especialmente importante para éste tipo de simulaciones. Los cálculos numéricos convergen a una solución oscilatoria para todos los casos analizados, ofreciendo así la posibilidad de un modelado avanzado de flujo pulsante e inestable a un costo numérico relativamente bajo. 117 118 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners 7.1 Conclusiones Generales La principal ventaja de los esquemas URANS es el bajo coste computacional, debido a que estos requieren de menos resolución temporal y espacial que el que necesita LES. Además, se ha demostrado que las URANS necesitan de pocas instantáneas para obtener una aproximación numérica confiable de un flujo pulsante. Uno de los aspectos importante en las simulaciones numéricas es la influencia de la malla sobre la precisión de la solución. El índice de convergencia de malla ha sido usado para estimar la incertidumbre asociada con los errores numéricos. Este análisis de error confirma que la mala ha sido construida adecuadamente. Con la ayuda de técnicas avanzadas de post-proceso, las estructuras coherentes presentes en el flujo han sido identificadas. Estas técnicas permiten observar que las estructuras coherentes rotan alrededor del eje axial. Los tipos básicos de estructuras coherentes que han sido observadas para el caso del quemador confinado de giro inducido son: una zona de recirculación del tipo burbuja y una estructura del tipo espiral que la envuelve. Mientras que para el flujo del quemador atmosférico de giro inducido se ha observado una zona de recirculación del tipo burbuja y una estructura del tipo helicoidal que la envuelve. Las inestabilidades producidas por estas estructuras han sido analizadas mediante la transformada rápida de Fourier de los monitores de presión estática y de velocidad. Chapter 7. Conclusiones 119 7.2 Efecto de las condiciones de contorno de entrada y de salida Fue puesto especial énfasis en la influencia de la posición de las condiciones de contorno en ambos quemadores. Para el análisis, simulaciones en 3D fueron ejecutadas bajo condiciones turbulentas e isotérmicas de un flujo monofásico. Fueron consideradas dos configuraciones para las condiciones de entrada: un dominio con el generador de giro y un dominio sin generador de giro pero con sus correspondientes perfiles de velocidades. Los resultados mostraron que las estructuras coherentes fueron mantenidas en ambos quemadores. Pero, solamente en el caso del quemador atmosférico de giro inducido, el promedio de la velocidad axial y los picos de frecuencia obtenidos mediante la estimación de la densidad de potencia espectral de los monitores de presión estática dieron el mismo resultado independientemente del tipo de condición de entrada. Esto significa que las inestabilidades producidas por el rompimiento de vórtice en los quemadores atmosféricos se propagan aguas abajo y nunca aguas arriba, En el caso del quemador confinado de giro inducido, las inestabilidades se propagan en todas direcciones. Es importante mencionar que la posición de la condición de contorno de salida también fue analizada. Aunque su influencia es menor comparada con la condición de entrada, se encontró una distancia óptima donde el efecto de la condición de salida fuera mínimo. 120 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners 7.3 Comparación de los modelos k−y esfuerzos de Reynolds El flujo simulado ha sido estudiado mediante varios métodos de post-proceso, incluyendo técnicas avanzadas para la educción de estructuras vorticiales coherentes, con el objetivo de comparar el desempeño de los dos modelos de turbulencia. Ambos modelos proporcionan predicciones realistas; describen un flujo complejo, casi-periódico con vórtices espirales y helicoidales que confirman una precesión de núcleo de vórtice y una zona de recirculación pulsante. Este comportamiento es capturado por los puntos de monitorización. Para el quemador confinado de giro inducido, cuando las señales son procesadas por FFT, el pico de frecuencia que predice el modelo k−estándar es distinto en un 35% comparado con la estimación del RSM, donde el pico de frecuencia es, además, más fuerte. Los resultados son muy convincentes y se comparan razonablemente bien con la literatura experimental. Resulta particularmente realista la predicción hecha por RSM, para el cual el flujo carece de simetría y periodicidad exacta, y exhibe movimientos verticales más fuertes y persistentes. En contraste, aunque con esencialmente las mismas características, el modelo k−estándar lleva a un flujo de naturaleza más esquemática. Para el quemador atmosférico de giro inducido, el promedio de la velocidad mostró que ambos modelos predicen valores cercanos a las mediciones experimentales. Para el caso de giro inducido bajo, la velocidad axial ajusta mejor para los datos experimentales con el RSM si se compara con el modelo k−estándar, Chapter 7. Conclusiones 121 mientras que para la velocidad tangencial, el modelo k−estándar proporciona una mejor aproximación. Para el caso de giro inducido alto, las velocidades axial y tangencial ajustan mejor con el modelo k−estándar. En lo que respecta a las series temporales, estos muestran un comportamiento similar al que está presente en el quemador de flujo confinado, es decir, las señales son periódicas. Para este caso, los picos de frecuencia fueron validados con datos experimentales donde la diferencia de la simulación numérica y los datos experimentales difieren solamente en un 4.4%. 7.4 Educción de estructuras coherentes Actualmente, existen diferentes técnicas para identificar las estructuras coherentes presentes en un flujo turbulento. Para esta tesis, dos técnicas comunes (isosuperficies de velocidad y vorticidad) y una técnica de visualización avanzada (λ2) fueron consideradas. Con todas las técnicas, fue posible identificar las zonas de recirculación, los vórtices en espiral, los vórtices helicoidales y algunas estructuras espirales y toroidales que forman alrededor del flujo principal. Sin embargo, debido a la complejidad de las estructuras coherentes y a su interacción entre ellas, solamente la técnica λ2fue capaz de mostrar con suficiente claridad todas las estructuras incluyendo las estructuras débiles. Las isosuperficies de velocidad y vorticidad no son capaces de detectar estas estructuras débiles. Las estructuras coherentes que se predicen mediante los dos modelos de turbulencia tienen grandes diferencias cuando son analizadas por la técnica λ2. Las estructuras coherentes dadas por el modelo k−estándar son débiles en comparación con 122 A Computational Fluid Dynamics Investigation of Turbulent Swirling Burners el modelo RSM. Por otro lado, el modelo RSM captura mejor los detalles de las estructuras coherentes presentes en el flujo, por ejemplo, la zona de recirculación interna la cual es deformada por la precesión de núcleo de vórtice (PVC). 7.5 Comparación de la simulación numérica con el S-PIV Simulaciones numéricas exhaustivas del flujo turbulento caracterizado mediante S-PIV por Legrand et al. [95] han sido presentados, bajo la premisa de la detección del fenómeno pulsante, rompimiento de vórtice y estructuras coherentes, y la investigación de la utilidad de los métodos económicos URANS en este contexto. Comparando los resultados numéricos obtenidos por RSM con las mediciones experimentales S-PIV, puede concluirse que el modelo de turbulencia es capaz de predecir el mismo número de núcleos. El flujo obtenido por medios numéricos se asemeja mucho al experimental, duplicando sus características. Flujo promedio, frecuencia pulsante y la intensidad de vórtices son adecuadamente predichos. La utilidad de los cálculos URANS para simular flujos complejos es evidente al menos en el uso como una herramienta complementaria de las mediciones experimentales, ya que producen toda la información necesaria para construir estructuras coherentes en 3D dentro de una aproximación aceptable.