Abstract
A computational model using Large Eddy Simulation (LES) for turbulence modelling was implemented, by means of the Eddy Dissipation Concept (EDC) combustion model using the fireFoam solver. A small methanol pool fire experiment was simulated in order to validate and compare the numerical results, hence trying to validate the effectiveness of the solver. A detailed convergence analysis is performed showing that a mesh of approximately two million elements is sufficient to achieve satisfactory numerical results (including chemical kinetics). A good agreement was achieved with some of the experimental and previous computational results, especially in the prediction of the flame height and the average temperature contours.
1. Introduction
Correct prediction and description of a fire have become one of the main concerns in safety engineering and risk analysis. Studies on fire safety have been developed mainly with emphasis on fire detection, heating of structures, and smoke-filling rates [1]. Usually this phenomenon has been analyzed through different experimental techniques. However, due to the destructive nature of fires, experimentation can be highly expensive and due to its randomness, it can be nearly impossible to replicate. Taking this into account, in recent years there has been an increasing interest in computational and numerical modelling of fires.
Nevertheless, there are several difficulties that arise when trying to fully develop models for fire dynamics simulations. First and foremost, the fact that fires consider several different physical phenomena that take place simultaneously such as turbulent flow, turbulent mixture processes, thermodynamics, heat transfer (especially through radiation, which in turn allows pyrolysis), and chemical kinetics [2]. Another difficulty is the coupling of the several analyses considered since these include a long range of length and time scales (for example the turbulent and chemical time scales). This is where different assumptions are applied on the combustion, chemical kinetics, and fluid dynamics processes [3].
These simplifications lead to the development of different models. When it comes to conflagrations (defined as an uncontrolled fire spread [2]), there are two main types of models: zone and field models. The first one basically divides the space in which the fire is taking place into two major zones: one in which predominantly remain the products of the combustion process and another one where the reacting air (oxygen) remains, until it is consumed by the reaction. On the other hand, field models are developed using Computational Fluid Dynamics (CFD) for reacting flows, in order to solve the Navier-Stokes equations coupled with the chemical kinetics solution (i.e., mixture fraction solution) [4].
While zone models are fairly straightforward, they have a major limitation in that these are only able to consider fires in an enclosure, therefore restricting the geometries and cases which can be explained through them. In contrast, field models have a higher mathematical complexity, while being able to adjust to almost any geometric domain and constraints. Therefore, there has been an increasing interest in this field to develop reliable field models which are accurate and exact enough to predict fires, complementing the experimental analysis. Some of the most notable field models are the Flame Surface Density model developed by Trouvé and Poinsot [5], the Partially Stirred Reactor model of Chen [6], or the Laminar Flamelet Model initially proposed by Peters [7].
One of the main combustion models used within field models is the Eddy Dissipation Concept (EDC), developed by Magnussen. In it, the author suggests a way to “relate the rate of combustion to the rate of dissipation of eddies” [8] assuming that the rate of reaction is a function of the mean concentration of a reacting specie, turbulent kinetic energy and its dissipation rate [9]. The EDC started as a model capable of considering both turbulent and momentum mixing, while considering the chemical kinetics solution and particularly soot formation. Most recently the interest in this model has shifted to develop an EDC that can account for Large Eddy Simulation (LES) turbulence models instead of the Reynolds Averaged Navier-Stokes (RANS) models traditionally used [9–12].
The biggest concern with the use of LES models is the fact that these filter the turbulent properties according to the different length scales for the eddies created (Integral, Taylor, Kolmogorov scale). Therefore, it becomes complicated to discuss a certain mesh convergence in computational models which use LES, since any change in the size of the cells of the mesh would also lead to a change in the size of the filter for the turbulent model, rendering different results which not necessarily tend to convergence [10].
Regarding the simulation of pool fires, there have been several studies in recent years. One of the most studied cases is the experiment proposed by Tieszen et al. in 2002 [13]. This experiment has been simulated by different authors, using mainly the Fire Dynamics Simulator (FDS). Such is the case of the simulations by Xin et al. in 2008 [14] and Cheung and Yeoh in 2009 [15]. Their results mainly show an interest in the flow’s velocity with great details of the fluctuating velocity field [14, 15]. Another example is the study of Novozhilov and Koseki [16] who use the FIRE software to simulate the results obtained in different experimental analysis including the experiment by Weckman and Strong [17]. Thus, they report results for the evolution of temperature with time, as well as normalized flame heights and burning rates for pool fires of different fuels and sizes. Finally, in 2014 Chen et al. [18] also recreate computationally the results obtained by Weckman and Strong [17], in an attempt to test the capabilities of the LES models developed in the solver fireFoam. Chen developed a new model for radiative emission, obtaining different results for the flame heights, heat generation, temperature, and emissive contours. However, some of these results do not agree completely with the experimental ones, particularly flame heights and the temperature root-mean-square contours.
Hence, the main objective of this work is to predict the dynamics of a small-scale pool fire and study the influence of the computational domain’s discretization especially in the vertical (flame propagation) direction, when using LES turbulence modelling with the EDC as the combustion model. As test case, the experiment proposed by Weckman and Strong in 1996 [17] was selected and simulated. The idea is to validate the numerical results obtained, by comparing them to the experimental results as well as other computational results found in literature. Numerical results show the influence of the mesh resolution on some variable contours and the flame height.
2. Methods
2.1. Experimental Setup
The base case used in the present computational study is the experiment developed by Weckman and Strong which consisted of a 30.5 cm methanol pool fire [17]. Weckman and Strong used different experimental techniques (for example Doppler anemometry) in order to measure a broad range of data: velocities, centerline and contour temperatures, mean and root-mean squares (rms) of fluctuations in velocities and temperatures, and flame heights [17]. This particular case has been well documented in literature, as different authors have developed computational analysis to validate numerical results (FDS [19], FIRE [16], SOFIE [20], and fireFoam [18, 21]).
The experimental facility consists of a pan burner of 30.5 cm in diameter where the methanol was injected at a rate of 1.35 cm3/s to assure a heat release rate of 24.6 kW. The burner is designed to minimize the obstruction of ambient air flow into the fire, such that it is guaranteed that the pan is at least one pool diameter above the floor. It is worth mentioning that a natural draft hood was used to exhaust the combustion products.
2.2. Computational Setup
In order to carry out the simulations, the open source code OpenFOAM–2.4.0 (OF) was used. This OF version was compiled with the fireFoam solver developed between CFD Direct and FMGlobal. This way the Navier-Stokes equations for compressible flows can be solved using the finite volume methods available as part of the OF library for turbulent flows. FireFoam uses a PIMPLE algorithm to iterate the calculation of the reacting flow properties. This iteration allows combining both the SIMPLE and PISO algorithm, such that it can use up to three correction equations for the pressure and velocity before solving the other transport equations [9].
2.3. Computational Domain, Mesh, and Boundary Conditions
The computational domain is a cylinder of 1.8 m in height and 1.8 m in diameter (grossly 6 pool diameters) to make sure that the boundary conditions are far enough from the pool base. This way the flame can develop freely, such that the pulsating phenomena are due only to buoyancy and not an effect of interacting boundary conditions. To generate the mesh, OF’s utility blockMesh was used. Hexahedral elements within the domain were used and a total of 90 cells across the diameter of the pool fire were used to run the simulations [18]. The mesh generation was restrained by the vertices of the geometry, which in turn forces a total of 270 cells across the domain’s diameter. However, there is no major study on the influence of the number of divisions in the vertical direction (flame propagation) in the numerical results. Four different meshes were generated in order to pursue this study and Table 1 shows the most important details of these meshes. All the meshes were refined towards the bottom face, so that there are more elements near the burner’s exit.
The minimum size of the cube root volume (CRV) for the coarse mesh is nearly 7 mm, which is smaller than the experimental value of the Taylor length scale [17]. As shown in Figure 1(c), in the region close to the burner, the cell volume is small enough to capture this length scale. Finally, Figure 2 shows a visualization and details of the coarse mesh (1.62 million cells).
(a)
(b)
(c)
(a)
(b)
Four different boundaries were defined: the base, sides, outlet, and an inlet patch on the “base” face. To make sure that the phenomenon could be modelled correctly, the boundary conditions were set according to Table 2.
The zeroGradient condition determines that the variable does not have gradients in the direction normal to that boundary. The inletOutlet condition acts differently depending on the flow’s direction. When the flow is going out of the domain it acts as a zeroGradient condition; however when the flow is entering the domain (backflow) it acts as a constant value given by the user. The inlet values used for the inletOutlet boundary condition were = 0 m/s, = 300 K, m2/s2, CH3OH mass fraction = 0, and O2 mass fraction = 0. The pressureInletOutletVelocity performs as an inletOutlet condition, but for backflow, instead of having a constant value, this is calculated from the pressure obtained from the internal cell value [22]. The kinetic turbulent energy was calculated according to the relationship: where is the turbulence intensity () and assuming isotropic turbulence. Typical values of turbulence intensity for small-to-medium scale pool fires can be between 5 and 15%. At the “Outlet” boundary, also has a value of m2/s2 when performing as an inlet condition.
2.4. Computational Models
The selected turbulence LES model for the present work was the one proposed by Menon et al. as the Spectral Eddy Viscosity model, which “assumes that the energy and energy transfer are in phase with each other in wavenumber space” [23]. The great advantage of this model is entirely computational; since it consists of only one equation then less computational cost is required. Even though this assumption is big enough to significantly deviate the results from the exact equations [23], previous computational studies [18] have proven that the computational results obtained using this model show a good agreement with experimental data.
In order to model the combustion, the Eddy Dissipation Concept proposed by Magnussen was chosen. This model basically expresses the turbulent mix as a function of the turbulent kinetic energy and its dissipation, the latter which is equal to the reaction rate. This leads to an expression of the mass fraction in terms of the two turbulent properties [24]. This way, the model only needs to solve one additional transport equation for the fuel’s mass fraction [9]. The EDC model is implemented in the fireFoam version that is to be compiled with OF. However, since it is based on Favre-averaged reaction rates [9], it has been mostly used with turbulent RANS models, therefore drawing a recent interest in coupling it to LES models.
Regarding the chemical kinetics solution, the irreversibleinfiniteReaction model was used which is available in the OF library. This model takes into account only one reaction (see (1)) to solve the chemistry of the case at hand, basically a one equation equilibrium model. By using a one reaction model, the chemical kinetics solution is greatly simplified since it limits the number of species (five in this case) thus reducing the amount of iterations and reaction rate equations the solver has to solve, therefore reducing the global computational cost of the simulation.
The finite volume Discrete Ordinates Model (fvDOM) was also included in the simulation in order to model the radiation heat transfer. In the present work, the model was implemented considering the gaseous phase as a gray gas participating media, while the Radiative Transfer Equation (RTE) is solved through a “loose coupling” method. This means that the RTE is solved using the solution of the previous time step (lagged in time). Then, the divergence of the radiative heat field is used as a source term for the energy equation which is solved in the PIMPLE (PISO) algorithm [25].
The simulations were set to have an adjustable time step, such that the Courant-Friedrich-Lewy number (CFL) would fluctuate around 0.5. In other words, on average the simulation would have a CFL of 0.5. This way the entire model depends on the specified CFL rather than on a fixed time step, which will allow solving the problem depending on the velocity of the fluid flow. This is relevant considering that the velocities in reactive turbulent flows are not constant; therefore the CFL would change in both, space and time, as the flow evolves. Finally, the PISO algorithm will be solved using a Gauss Linear scheme, with an implicit Euler time discretization.
3. Results
Figures 3 and 4 show instantaneous temperature contours (in the midplane and 3D isosurfaces) for the coarse (C) mesh. Particularly, these images show the buoyant phenomenon as the principal transport mechanism for the reactive flow, therefore, suggesting the correct implementation of a pool fire model.
(a)
(b)
(c)
(d)
(a)
(b)
Since the instantaneous results vary from one simulation to another, it was necessary to look at the time averaged results for the different variables of interest, that is, temperature, velocity and mass fractions as shown in Figures 5 and 6 for meshes C and SF, respectively. The results shown are based on averages taken in a time frame window of 20 s.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
By comparing Figures 5 and 6, it is clear that the finer mesh is able to obtain more detailed results as expected. For example, by comparing images (c) and (d) from both figures, the SF mesh is able to detail the chemical kinetics solution for the higher part of the domain. This also occurs with the velocity profiles (image (b) in both figures). In this case, the coarse mesh only obtains velocities up to 2.87 m/s while the SF mesh goes up to 3.18 m/s. Also, it is worth noting that the top velocities cover a higher range in the SF mesh than the C mesh (i.e., the red contour is larger in Figure 6 than in Figure 5). On the other hand, the mass fraction contours show in general the same values, with the exception that these also have a better resolution in the higher part of the domain. It is noteworthy that the profiles for oxygen (reactant) are a bit smaller in comparison to the carbon dioxide (product). This can give insight to the outermost points at which the reaction is taking place (those points where both species show a contour), and at which location begins to govern the hot gases from the combustion process. Looking at the temperature contours, the higher temperatures also cover a larger spatial region in the SF mesh than the C mesh; however the latter shows a value for the maximum temperature 13% higher than the finer mesh. Finally, it is interesting to notice that the maximum values for each contour (the lowest values for the oxygen contours, since it is a reactant) along the axis of the domain have similar location.
The flame height for the simulations was determined using four different approximations: image analysis, mixture fraction calculation, fireFoam’s flame height tool, and Hekestad’s empirical method [1]. The first one was performed by visually locating the highest point at which the products of combustion had the highest mass fraction, based on the average mass fraction contours (images (d) in Figures 5 and 6). This would in turn show the point at which the reaction was still taking place, therefore representing the visible flame.
The second method estimates the average mixture fraction () given by (2) [3]:Figure 7 shows the time averaged values for the mixture fraction along the axial position of the flame. From this result, the flame height can be defined as the position at which , which corresponds to the stoichiometric value.
The third method used to estimate the flame height is a tool included in the fireFoam solver. In this case the flame height is calculated using (3) [18]:where is the vertical coordinate (flame propagation direction), while where the mass fractions are density-weighted [18].
The last method used to estimate the flame height is entirely empirical (see (4)) in which the flame height is related to the total heat generated by the fire and the pool fire’s diameter [1]Finally, Weckman and Strong determine the flame height experimentally as the highest point in which the flame is still visible, using instantaneous photography. The experimentally measured flame height is not higher than 0.3 m [17]. However, it is worth mentioning that a methanol flame burns blue (no soot formation); thus it can be difficult to accurately determine the height visually [18].
3.1. Mesh Convergence
Flame height is then the criterion used for mesh convergence, as it is shown in Figure 8 and Table 3 where the flame heights for the different methods are collected.
It is clear that the flame height has an asymptotic behavior showing a convergence. For the visual contour method, the differences between the experimental and numerical results are only 8.6% and 5.3% for the coarse and superfine meshes, respectively, with a minimum of 3% for the medium mesh. When looking at the results obtained by the mixture fraction method, a difference of 3% with respect to the reported experimental data is observed for the (SF) mesh.
The fireFoam’s tool for flame height calculation (mass fraction method) highly overpredicts the experimental result of 0.3 m. This results is expected since the definition of includes all the region in the computational domain in which mass fractions of fuel and oxidant exist but it does not consider the mass fractions of the products. The boundary of this region is larger than the region in which the stoichiometric value of the mixture fraction (second method) is observed.
Figure 8 suggests that there is a convergence to the model proposed. Based on these results, the medium (M) mesh was selected as the one that better compromises numerical precision with computational cost. This is also shown when looking at the normalized results, which according to McCaffrey characterize the behavior of a pool fire. Normalized results as shown in Figures 9 and 10 present a better agreement with the experimental results reported, even better than other computational results reported in literature [18].
It is evident that there is actually a behavior that tends to convergence when the computational mesh is refined especially in the region far from the burner exit. Figure 9 also shows the expected behavior as suggested by Heskestad [1], when reaching the region in which the slope of the curve far from the burner exit approaches −5/3. Figures 11–14 show the average axial velocity as function of the radial position at different vertical position (6, 10, 18, and 30 cm). Better agreements between the computational and experimental results are observed far from the burner exit (Figure 14)
The convergent behavior of the numerical results can also be appreciated when looking at the average temperature contours. To do this, 12 contours were selected ranging between 300 K and 1092 K (the highest average temperature obtained for the finer mesh); therefore, from the inside out, the contours increase every 66 K. In Figure 15, the colored contours show the superfine mesh (2.84 M cells) while the black contours show the coarser meshes.
(a)
(b)
(c)
However, when comparing the results obtained in the simulations, with the ones obtained experimentally by Weckman and Strong [17], it is clear that the random behavior of the fire phenomena cannot be replicated to its full extent. The contours are different not only geometrically, but also numerically. The experimental results show values up to 1300 K, while the only simulation to arrive to such a value is the coarsest mesh. However, there are some observations that can be found from this comparison. Although the geometry of the contours is not similar, the general behavior is. For example, the contours near the burner’s exit tend to flatten as they approach the burner exit. Also, consider the fact that the high temperature areas are similar; that is, the 1000 K contour in both the experimental and numerical results achieves a height of about 25 cm and radial position of 6 cm. Therefore, although the numerical results are somewhat different from the experimental ones, the general behavior of the pool fire seems to be correct. Particularly, since the high temperature (above 1000 K) area is contained in the same regions (below 25 cm high and 6 cm wide). It is evident that the medium (M) mesh behaves similar to the finer meshes. While the coarse (C) mesh does not show any bending of the temperature contours close to the burner exit, the medium (M) mesh captures part of this behavior.
One of the main advantages of the fireFoam solver is that it is able to solve the chemical kinetics of this particular combustion process. These results are shown in Figure 16.
The highest consumption of methanol (the fuel) happens at the same time as the biggest increase in water vapor and carbon dioxide (dotted line on Figure 16). Also, consider the fact that the reaction effectively reaches an equilibrium, due to the one reaction system used as part of the Irreversible Infinite Reaction model for the chemical kinetics solution of this case.
Finally, Table 4 shows that the simulation on the coarse (C) mesh took approximately 180 hours of CPU time running 12 parallel processes in an Intel Xeon processor E5-2695 v2, while the finer (SF) mesh took 267 hours of CPU with the same number of parallel processes. This again proves the necessity of the different models used, considering the complexity of the combustion process. Especially the turbulent and chemical kinetics models used are necessary to significantly reduce the computational cost.
4. Conclusions
In the present work, simulations were carried out in order to replicate the experimental results obtained by Weckman and Strong [17] with a small-scale methanol pool fire, through the use of the open source code fireFoam. The numerical results show a good agreement with results such as flame heights and normalized results. Nevertheless, when looking at time averaged results and particularly temperature contours the results show some discrepancy, once again proving the random nature of fires and turbulent reactive flows. Again, this is proved by comparing the numerical results of the contours of average temperature in the midplane with the experimental ones; this shows not only numerical but also geometrical discrepancies. However, the general behavior of these results shows that the models that were implemented can simulate effectively a pool fire. The regions with high temperatures are similar in both cases (experimental and numerical), or the fact that the chemical reaction takes place in similar regions (flame height results). Also, mesh convergence was analyzed by running simulations with four different mesh sizes. This proved that effectively the combustion and turbulent models have an asymptotic behavior as a function of the number of elements that compromise the computational domain. Based on the convergence analysis, it is appropriate to affirm that a mesh with 90 divisions across the pool fire’s diameter and 50 divisions in the axial direction (medium size mesh) is fine enough to obtain reliable results; an error of 3% is obtained with respect to flame height while there is a negligible difference when looking at Hekestad’s normalized results (Figure 9) and the temperature contours (Figure 15). Taking all of this into account, the fireFoam tool has proven to be an effective solver for the numerical simulations of pool fires. This can in turn lead to suggesting the use of the computational tool for large scale diffusive turbulent flames. For example, in the prediction, prevention, or a posteriori analyses of large scale fires, it is capable of reliably calculating different variables interest such as temperatures, velocities, and species mass fractions.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors would like to thank the engineering schools of both Universidad de los Andes and Universidad del Valle (Grant CI2868) for providing the resources and support for the development of the present research and, also, the HPC support team at the Universidad de los Andes, for providing the cluster computing service in which the simulations took place.