Computing flare dynamics using large eddy simulations

Computing the dynamics of flares is motivated by the increased need for efficient and safe flaring of unwanted gases during hydrocarbon and petrochemical processing. To understand the unsteady flame shape dynamics, Large Eddy Simulations (LES) are used to study natural gas flares in the presence of crosswind. Different flare regimes are created by the competing forces of the jet inertia, flare buoyancy, and crosswind inertia. These regimes are studied by performing LES simulations which include different flare jet velocities, crosswind velocities and flare tip diameters. As the ratio between the flare jet to crosswind momentum decreases for constant flare tip diameters, the computed flame shape is changed from a buoyancy dominated regime to a fully wake stabilized regime. These simulation results are in agreement with the flame shapes observed in wind tunnel experiments reported in the open literature. With this confidence in the simulation tool, a parametric study is carried out to understand the effect of flare tip diameters and velocities typical of large-scale industrial flares. In addition to the jet to crosswind momentum ratio, the individual velocity magnitudes also play a significant role in changing the flame regimes. This study shows the capability of the LES method to predict dynamics of flare flame shape.

C o m p u tin g F la r e D y n a m ic s U sin g L a r g e E d d y S im u la tio n s Padmabhushana R. Desam1, Philip J. Smith, Stainslav G. Borodai and Seshadri Kumar Combustion and Reaction Simulations (CRSIM) Research Group University of Utah Salt Lake City, UT-84112 ABSTRACT Computing the dynamics of flares is motivated by the increased need for efficient and safe flaring of unwanted gases during hydrocarbon and petrochemical processing. To understand the unsteady flame shape dynamics, Large Eddy Simulations (LES) are used to study natural gas flares in the presence of crosswind. Different flare regimes are created by the competing forces of the jet inertia, flare buoyancy, and crosswind inertia. These regimes are studied by performing LES simulations which include different flare jet velocities, crosswind velocities and flare tip diameters. As the ratio between the flare jet to crosswind momentum decreases for constant flare tip diameters, the computed flame shape is changed from a buoyancy dominated regime to a fully wake stabilized regime. These simulation results are in agreement with the flame shapes observed in wind tunnel experiments reported in the open literature. With this confidence in the simulation tool, a parametric study is carried out to understand the effect o f flare tip diameters and velocities typical of large-scale industrial flares. In addition to the jet to crosswind momentum ratio, the individual velocity magnitudes also play a significant role in changing the flame regimes. This study shows the capability of the LES method to predict dynamics of flare flame shape. f Corresponding Author: [e-mail] IN T R O D U C T IO N The primary objective of flare systems is to burn unwanted combustible gases released during hydrocarbon and petrochemical processing operations in an effective and safe manner in the open atmosphere [1]. The need for effective disposal of gases stems to meet the environmental regulations since methane is 20 times more effective in trapping heat in the atmosphere, thus results in potential global warming, compared to the combustion products such as CO2 [2]. In operating the flares, radiation threat to the surroundings is a major safety concern. The amount of radiation to the surroundings depends on the location of the flame shape, size and amount of soot. Experimental studies using photographic techniques showed that the flame shape and size depend on factors like ambient wind conditions, flare velocities and fuel composition [3]. A priori knowledge of the flame shape and size for different conditions facilitate the design of flare stack heights and address several safety related issues. A common measure to describe the effectiveness of combustion in flares is combustion efficiency, which is defined as the ratio between the carbon present in unburned flare gases to the carbon in products in the form of CO2 [3]. Previous studies on measuring the combustion efficiency of laboratory scale flares in closed-loop wind tunnels showed that the wind conditions play a significant role [3-5]. The reason for this is the dynamics of turbulent combustion resulted by the local mixing between fuel and oxidizer. Understanding the unsteady mixing phenomena of flares is inevitable to design the flare equipment for achieving higher carbon efficiency. Experimental studies on measuring carbon efficiency, flame shape and size are only limited to laboratory scale flares [3-6]. However, the wide variety of fuel composition, fuel velocities, ambient wind conditions and the size of the flare equipment make experimental studies not feasible for designing industrial flare equipment. Recent advances in computational combustion aid to tackle this problem cost-effectively using numerical simulations. In order to use the simulations with confidence in design and operation of flares, the simulation results must be compared with the experimental results to determine whether the physics of the flares are accurately represented. 2 Industrial flares operate under turbulent flow conditions, which involve wide ranges of length and time scales. The biggest of these scales can be on the order of diameter of the flare tip and the smallest is determined by the viscosity. In general, the fuel from the flare tip and ambient air as oxidizer come together in separate streams before mixing takes place. The chemical reactions take place after the fuel and oxidizer is mixed at molecular level. This kind of configuration is referred as nonpremixed combustion where the reaction rates are limited by mixing rates. Detailed chemical mechanisms of combustion reactions involve several thousand elementary reactions steps among hundreds of species with a wide range of time scales. These chemical reactions are exothermic, resulting in both convective and radiative heat transfer. Radiation is the dominant mode of heat transfer at high temperatures, whose transfer depends on the absorptive, emissive, and scattering properties of the gas mixture. The presence of soot in the flame enhances the rate of radiative heat transfer. All of these processes are highly coupled. For example, turbulence enhances mixing and thus chemical reactions. Chemical reaction changes the temperature through the amount of heat generated, and this changes the density and thus the intensity of mixing via turbulence. Resolving all the length and time scales in practical turbulent combustion applications is not possible even on supercomputers for the reasons mentioned in the previous section. In general, time-averaged or spatially-filtered governing equations are solved with subgrid models to represent the unresolved scales and interactions among them. By time-averaging the equations (Reynolds Averaged Navier-Stokes-RANS), unsteady information such as instantaneous mixing and flame shape can not be captured. Instead, important features of the flame can be captured by resolving large length and time scales responsible for controlling the dynamics, with models for more homogenous smaller scales. This method is called Large Eddy Simulations (LES) and a more detailed description of this method follows this section. However, resolving even the large length and time scales requires large amounts of computational resources. So, this task has to be done in a massively parallel environment with the aid of computational tools that are scalable. Over the past 5 years, the Combustion and Reaction Simulations (CRSIM) research group at the University of Utah has been developing a massively-parallel LES code (ARCHES), which is linearly scalable up to 1000 processors, in a component-based 3 framework called Uintah Computational Framework (UCF). The UCF provides the framework for large-scale parallelization of different applications [7]. In this study, LES of natural gas flares are performed to understand the dynamics of flame shape and size. Different flare jet velocities, crosswind velocities and flare tip diameters are considered to simulate different regimes of laboratory scale flares as described in the reference [3]. The computations are compared with the corresponding experimental data. To understand the dynamics of large industrial scale flares, computations are performed with large flare tip diameters and velocities representative of industrial-scale flares. The following section describes the computational procedure used in the simulations. C O M P U T A T IO N A L P R O C E D U R E LES is employed to model the fluid dynamics and the convection-diffusion scalar transport in the flares. This method successfully captures the transient nature of the coherent vortical structures responsible for the dynamic features of mixing and flame shape. The LES approach is based on solving the set of density-weighted, filtered, timedependent, coupled conservation equations for mass, momentum, mixture fraction, and enthalpy in a Cartesian coordinate system as shown below [8]. I p + V ^ (P J) = 0 (1) dpV + V • (p J J ) = - V p + V • a - V • (p(VV - J J )) dt (2) dt P ^ dt dt + V • (PJ Z ) = V • (pz~z VZ) - V • (P(uZ - DZ)) (3) + V • (PJA ) = V • (P rV J ) - V • (P(Vh - » ) ) - V • (q „d ) (4) Here, p is the density, t is the time, v is the velocity vector, p is the pressure, Z is the mixture fraction, DZ is the molecular diffusion coefficient of the mixture fraction, a is the viscous stress, h is the enthalpy, r is the thermal diffusion coefficient, and qrad is the radiative flux vector. 4 These filtered equations are discretized on a 3-dimentional, structured, Cartesian staggered grid using a second-order differencing scheme. An explicit second order strong stability preserving Runge-Kutta time integration scheme [9] is used for advancing the variables in time. It is reasonable to assume that acoustic waves do not play a significant role in the dynamics of the flares. Therefore, the low-Mach number variable density formulation analogous to one by Najm et al. [10] is used for the present calculations. In this approach, the intermediate velocity is computed from the pressure-free momentum equation. Poisson equation for pressure is then obtained from mass conservation considerations. After pressure is solved for, the intermediate velocity is corrected to obtain velocity field satisfying mass conservation. One of the important characteristics of the flares addressed here is that it occurs in an open domain. Because the open domain must be modeled with a finite numerical domain, pressure-based open domain boundary conditions are employed on all the boundaries except for inlet and crosswind boundary. Boundary conditions for the open domains available in literature [11, 12] proved to be inadequate, so new formulation for the open domain boundary conditions have been developed. At the inlet boundary, the boundary condition is a constant specified mass flux with initially constant velocity profile. At the crosswind boundary, a constant crosswind velocity across the boundary is specified. The variables of interest in radiative heat transfer analysis are the distributions of radiative heat flux vectors ( q(r) ) and the radiative source terms (V • q rad). The radiative source term describes the conservation of radiative energy within a control volume and goes into the total energy equation (Eq.4), thereby coupling radiation with the other physical processes that occur in a multi-physics application. In order to determine the distributions of these quantities, a method is first required for the solution of the radiative transfer equation (RTE). The solution to the RTE gives distributions of spectral intensities. The heat fluxes and the source terms are readily computed by performing integrations over the spectral intensities. Discrete Ordinates (DO) method is employed to solve the RTE. More details on the equations solved and solution procedure can be found 5 in Ref [13]. Soot is included in calculating the radiative properties along with the gases participate in radiative heat transfer such as CO2 and H2O. The unresolved subgrid scale Reynolds stress (momentum equation), the scalar mass flux (mixture fraction equation), and the turbulent energy flux (enthalpy equation) are modeled using the Smagorinsky approach [14]. In this formulation, the subgrid scale fluxes in the momentum and mixture fraction transport equations are given by p (u u - u u ) = - 2 p v tS (5) ~ 1 ~ ~ T S = _((V ~) + (V ~ )T (6) p(uZ - ~Z) = -p D t VZ (7) with and where vt is the subgrid kinematic eddy viscosity and Dt is the subgrid eddy diffusivity. The eddy viscosity vt is given by the Smagorinsky model as vt = CA2|S (8) where C is the Smagorinsky constant and A is the filter width. The subgrid diffusivity Dtt is determined from Dt = vt / Sct (9) where Sct is the turbulent Schmidt number and assumed to be a constant in the current calculation. Similarly, turbulent energy diffusivity is computed assuming a constant turbulent Prandtl number. The Smagorinsky constant is set to 0.17 for all of the calculations. A value of 0.4 is used for both the turbulent Prandtl and Schmidt number constants. In nonpremixed combustion, the chemical reactions associated with combustion are fast compared to the subgrid scale mixing of fuel and oxidizer, thus limits the rate of reactions. Under this assumption and the assumption of equal diffusivity of all species, 6 state space variables such as temperature, density and species concentrations can be expressed as functions of mixture fraction (Z), a conserved scalar [15]. In order to include the effect of enthalpy loss by radiation on the state space variables, an additional parameter, heat loss is also included. Heat loss is a normalized variable, describes the amount of sensible enthalpy lost/gained by heat transfer. A non-adiabatic equilibrium model is applied to relate the instantaneous state space variables with the mixture fraction and heat loss [16]. Subgrid scale mixing is modeled by assuming a P-PDF shape in the mixture fraction. Filtered values of the state space variables are computed by convoluting the variables over this prescribed PDF shape. Computing the shape of the P-PDF requires the mean and the variance of the mixture fraction values. Since no transport equation for the subgrid scale variance is solved, it is modeled using the assumption that the subgrid scale variance production and dissipation rates are in equilibrium [17]. ARCHES is developed to address gaseous combustion problems. The presence of the flare pipe in the domain is accounted for by the Material Point Method (MPM). [18] MPM uses material points (i.e., point masses) to describe the state (i.e., velocities, temperatures, stresses) of the solid. In this paper, the state of the solid is kept fixed, and so the only effect of MPM, as far as the solid is concerned, is in describing its geometry. The momentum and heat transfer exchange due to the fire-structure interaction are accounted for in the fluid phase, but not in the solid phase. In other words, the temperature rise in the solid pipe is neglected. This is acceptable since the pipe has a high heat capacity, and therefore it would take a long time for the temperature of the pipe to rise significantly. The pipe does, however, radiate heat to the fire, through its specification as a boundary at 298 K. The fire radiates energy to the pipe, and this heat loss is calculated, along with the convective heat transfer implicit in the presence of the wall cell. R E S U L T S & D IS C U S S IO N The computational grid in most cases consists of 2,625,000 cells, with uniform spacing in each direction. Computations are performed with 48 processors on a parallel 7 Linux cluster consists of Pentium IV 2.4Ghz processors. Typically, the computational time varied between a minimum of 72 hours to a maximum of 148 hours. Natural gas with the composition given in Ref [3] is used as fuel. As a part of the flame shape validation, three cases are considered that match the experimental conditions as described in Ref [3]. These three cases demonstrated distinctively different flame shape and size due to the jet and crosswind velocities. For comparison purposes, a non-dimensional parameter to describe the relative magnitudes of density-weighed momentum flux ratio as defined in Ref [3] is also adopted in this study, which is Pj V R =j r paU a (1 0 ) In the above equation, pj is the density of fuel, Vj is the velocity of flare jet, p„ is the crosswind density and U„ is the crosswind velocity. Figures 1 to 3 show the flame shape from computations for different values of ‘R ', with a constant flare tip diameter. The flame shapes shown in these figure are extracted after the flame is fully developed and reached a quasi-steady state. For a value of R=5.78 (Fig. 1), the effect of crosswind is only limited to the near tip region of the flame. Beyond that, buoyancy effects dominate because of the low density properties of the combustible gases. Even, the jet momentum is insufficient to overcome the buoyancy forces and results in a lazy buoyant flame. As the value of ‘R ' is decreased to 0.15, as shown in Fig. 2, the crosswind affects the flame in different ways. A wake is formed behind the flare pipe which intensifies the mixing of fuel and oxidizer due to the low pressure region. This results in part of the flame attaching to the flare pipe. Beyond this region, it can also be observed that the crosswind momentum washes out the flame because of its high momentum. The effect of buoyancy results in the flame being lifted from its horizontal position. A further decrease in ‘R ' value to 0.021 (Fig. 3) results in the flame being fully attached to the wake region and is referred to as being wake stabilized. In this case, the buoyancy has negligible effect on the flame shape and the flame is aligned horizontally. As ‘R ' decreases in these cases, the flame length is initially increased and then decreased for a very low value of ‘R '. All of the above observations are in agreement with the experimental results, also shown in the figures. From this study, it is concluded that simulation tool successfully 8 Figure 1. Flame shape comparison between simulations and experiments for R=5.79; (Left) Volume rendered temperature field from the LES simulation; (Middle) Planar image of the temperature field from the LES simulation through the center of the flame; (Right) Photographic image from experiments [3]. Figure 2. Flame shape comparison between simulations and experiments for R=0.15. (Top) Volume rendered temperature field from the LES simulation; (Middle) Planar image of the temperature field from the LES simulation through the center of the flame; (Bottom) Photographic image from experiments [3]. 9 (c) (d) Figure 3. Flame shape comparison between simulations and experiments for R=0.021; (a) Volume rendered temperature field from the LES simulation; (b) Planar image of the temperature field from the LES simulation through the center of the flame; (c) Photographic image from experiments (short exposure); (d) Photographic image from experiments (long exposure) [3]. captures the dynamics of the flame shape and size. Apparently, the flame shape is dominated by large scale structures in the flow. LES is able to capture this dynamic phenomena. The advantage of the simulation tools is its ability to simulate different conditions by changing the boundary conditions and the computational domain. Thus, further studies were carried out to study the effect of large flare tip diameters and velocities on flame shape and size. To achieve the same task with the experiments would be a difficult task since an increase in flare tip diameter needs a large wind tunnel and a powerful fan to create the wind conditions. The remainder of this section discusses the results obtained for large flare tip diameters. 10 To study the effect of scaling flare equipment on flame shape, the diameter of flare tip is increased 10 times. A crosswind velocity of 7.39 m/sec and a jet velocity of 3.78 m/sec are considered, which corresponds to an ‘R ' value of 0.15. These parameters represent typical conditions for industrial scale flaring equipment [5]. Figure 4 shows the flame shape for this configuration. The crosswind momentum dominates the near tip region of the flame, without creating any flame in the wake region of the flare pipe. On moving away from the flare tip, buoyancy forces play a more significant role and dominate the tail region of the flame. In contrast to the partly wake attached flame shown in the Fig. 2 for the same ‘R ' value, this flare fully resides above the flare tip. This result indicates that the flame shape does not scale with ‘R ' alone but shows a diameter dependency. The effect of increasing both the crosswind and jet velocities without changing the value of ‘R ' is also studied. In this case, a crosswind velocity of 12.0 m/sec Figure 4. Flame shape for larger diameter flare tip with a jet velocity of 3.78 m/sec and crosswind velocity of 7.39 m/sec (R=0.15). (Top) Volume rendered temperature field from the LES simulation; (Bottom) Planar image of the temperature field from the LES simulation through the center of the flame; 11 Figure 5. Flame shape for larger diameter flare tip with a jet velocity of 6.13 m/sec and crosswind velocity of 12.0 m/sec (R=0.15); (Top) Volume rendered temperature field from the LES simulation; (Bottom) Planar image of the temperature field from the LES simulation through the center of the flame; and a jet velocity of 6.13 m/sec are used. The flame shape for these conditions is shown in Figure 5. In spite of the constant ‘R ' value, the flame is distinctively different from the previous case shown in Fig 4. These conditions create a wake behind the flare tip with part of the flame attached to the pipe. The effect of buoyancy forces is limited to the tail region of the flare. Also, the length of the flame is increased compared to the previous case. This flame shape is qualitatively similar to the one shown in Fig.2 for R=0.15 with less shedding frequency. C O N C L U S IO N S Natural gas flares with different conditions are studied using LES to understand the different flame shape regimes such as buoyancy dominated, wake attached and wake stabilized. A set of calculations is performed to validate the computed flame shape by matching experimental conditions. There is a good agreement between the simulations 12 and experiments, which gives confidence for subsequent parametric studies. The effect of large flare tip diameters and different magnitudes of crosswind and jet velocities are studied using the simulation tool. An increase in flare tip diameter changes the flame from partly wake attached to fully reside over the tip. Increasing crosswind and jet velocities resulted in having part of the flame develop in the wake region. These parametric studies show evidence that the flame shape also depends on the diameter of the flare tip and individual velocity magnitudes of flare jet and crosswind. ACKNOW LEDGEM EN T Our sincere thanks to Dr. Larry Kostuik of University of Alberta for sharing the photographic images of flares in Ref [3]. REFEREN CES 1. Charles, E. Baukal Jr., "The John Zink Combustion Hand Book", CRC Press, p. 589­ 634, 2001. 2. U.S. Environmental Protection Agency, "Greenhouse Gases and Global Warming Potential Values", http://www.epa.gov. 3. Johnson, M. R., and Kostiuk, L. W., "Efficiencies of Low-Momentum Jet Diffusion Flames in Crosswinds", Combustion and Flame, Vol. 123, pp.189-200, 2000. 4. Bourguignon, E., Johnson, M. R., and Kostiuk, L. W., "The Use of Closed-Loop Wind Tunnel for Measuring the Combustion Efficiency of Flames in a Cross Flow", Combustion and Flame, Vol. 119, pp.319-334, 1999. 5. Johnson, M. R., and Kostuik, L. W., "A Parametric Model for the Efficiency of a Flare in Crosswind", Proceedings of the Combustion Institute, Vol. 29, pp. 1943­ 1950, 2002. 6. Majeski, A. J., Wilson, D. J., and Kostuik, L. W., "Size and Trajectory of a Flare in a Cross Flow", Combustion Canada 1999, Calgary Alberta, May 1999. 7. Davison de St. Germain, J., McCorquodale, J., Parker, S. G., Johnson, C. R., "Uintah: A Massively Parallel Problem Solving Environment," Ninth IEEE International Symposium on High Performance and Distributed Computing (HPDC) , August 2000. 13 8. Rawat, R., Spinti, J. P., Yee, W., and Smith, P. J., "Simulation of Large Hydrocarbon Pool Fires on Massively Parallel Computers", 7th International Symposium on Fire Safety Science, Massachusetts, June, 2002. 9. Gottlieb, S., Shu, C. W., and Tadmor, E., "Strong Stability-Preserving High-Order Time Discretization Methods", SIAM Review, Vol. 43, No.1, pp.89-112, 2001. 10. Najm, H. N., Wyckoff, P. S., Knio, O. M., "A semi-implicit numerical scheme for reacting flow," Journal of Computational Physics, Vol. 143, No. 2, pp. 381, 1998. 11. Boersma, B. J., Brethouwer, G., Nieuwstadt, F. T. M., "A numerical investigation on the effect of the inflow conditions on the self-similar region of a round jet," Physics of Fluids, Vol. 10, No. 4, pp. 899-909, 1998. 12. Akselvoll, K. and Moin, P., "Large-Eddy Simulation of Turbulent Confined Coannular Jets," Journal of Fluid Mechanics, Vol. 315, pp. 387, 1996. 13. Gautham Krishnamoorthy, Rajesh Rawat and Philip J. Smith, Parallel Computations of Radiative Heat Transfer Using the Discrete Ordinates Method, Numerical Heat Transfer, Part B: Fundamentals (accepted for publication). 14. Mason, P. J., "Large-eddy simulation: A critical review of the technique," Q. J. R. Meteorol. Soc, Vol. 120, pp. 1-26, 1994. 15. Bilger, R. W., "Turbulent Flows with Nonpremixed Reactants," Topics in Applied Physics: Turbulent Reacting Flows, Libby, P. A. and Williams, F. A., Eds., Springer Verlag: New York, 1980. 16. Sivathanu, Y. R., and Faeth, G. M., 1990, "Generalized State Relationships for Scalar Properties in Nonpremixed Hydrocarbon/Air Flames," Combustion and Flame, Vol. 82, pp.211-230. 17. Pierce, H. and Moin, P., "A dynamic model for subgrid-scale variance and dissipation rate of a conserved scalar," Physics of Fluids, Vol. 10, pp. 3041, 1998. 18. Guilkey, J.E., and Weiss, J.A., "Implicit Time Integration for the Material Point Method: Quantitative and Algorithmic Comparisons with the Finite Element Method," Int.J.Num.Meth.Engg.,Vol 1., p.1 (2002). 14