Large eddy simulations of flame extinction with infinitely fast chemistry

Large eddy simulations of flame extinction with N2 as extinguish agent are performed focusing on combustion and radiation modelling with infinitely fast chemistry. The use of a EDC combustion model with dynamically determined coefficients, an enthalpy-based flame extinction model based on a locally variable critical flame temperature and the use of the WSGGM for radiation are employed in order to predict flame extinction in a turbulent CH4 line burner. The numerical predictions of mean temperatures, combustion efficiencies and radiative fractions with different grid sizes are compared to the experiments by White et al. (2015). Overall, the results from the numerical simulations agree well qualitatively and, to some extent, quantitatively with the experimental data when small grid sizes are employed. More specifically, the maximum values and the profile widths of the mean temperatures at two axial locations examined are reasonably well predicted. The decrease in the combustion efficiencies as the extinction limit is approached is reproduced by the numerical simulations. The decreasing trend in the radiative fractions as the oxidizer stream is diluted with N2 is also captured by the simulations using a WSGGM model for radiation with a dynamically determined beam length which is calculated based on the local heat release rate. Nevertheless, the resulting radiative fractions are still over-predicted as the extinction limit is approached. Limitations of some of the typically used approaches regarding radiation modelling in flame extinction scenarios are outlined. Additionally, possible deficiencies on the use of a fixed critical flame temperature and/or re-ignition temperature in the flame extinction/re-ignition modelling are discussed and possible extensions of the modelling approaches are presented.


INTRODUCTION
Flame extinction is a complex physical process which has received a lot of attention from the scientific community in the past and is inherently important when it comes to fire safety engineering since it deals with the protection of material property (i.e., industrial and/or residential) but also the life safety of the public. The use of Computational Fluid Dynamics (CFD) in this area of research is crucial and its advantages are multifold. The use of CFD can be used towards evaluating the performance of current fire detection systems, the design and testing of new ones but also the virtual simulation of a wide range of hypothetical fire scenarios that in reality would either be too expensive to conduct or would not be permitted for fire safety reasons.
Extinction can be broadly described as the reduction of the chemical reaction rate below a critical value for which self-sustained combustion cannot occur. In general, three different mechanisms of extinction can be defined: aerodynamic quenching (i.e., reduction of the flame residence times through flow perturbations), thermal quenching (i.e., radiative heat losses and/or water droplet evaporation) and quenching by dilution (dilution of the fuel and/or the oxidizer stream). Given the type of scenarios usually involved in fire applications, thermal quenching and quenching due to dilution are expected to be more the predominant mechanisms for flame extinction as opposed to aerodynamic quenching. CFD models should, in practice, be applicable and produce reliable and accurate predictions for a wide range of conditions/scenarios if they are to be used for practical applications of fire safety engineering. Currently, fire-related CFD codes mainly rely on the concept of a critical flame temperature or critical Damköhler number while other models have been recently presented as well such as, e.g., the reactive volume fraction model.
The work presented here aims to continue the authors' previous work [1] towards predictive fire modelling. The main objective of the paper is the application of advanced modelling approaches related to turbulence, combustion and radiation modelling along with the use of an enthalpy-based extinction model using a variable critical flame temperature and their evaluation when it comes to predicting flame extinction with infinitely fast chemistry. For this purpose, Large Eddy Simulations (LES) of a turbulent CH 4 line burner with N 2 as extinguish agent are presented and the results are compared to the experiments by White et al. [2]. The chosen experiments are very relevant for the fire research community since they are one of the target test cases of the newly established Measurement and Computation of Fire Phenomena (MaCFP) workshop (https://github.com/MaCFP), focusing on advancing predicting fire modelling, making them suitable for evaluation/validation purposes of CFD codes.

MODELLING
The standard version of FireFOAM 2.2.x (https://github.com/fireFoam-dev), developed by FM Global, is employed here. The code uses a fully compressible flow formulation and solves the Navier-Stokes equations, along with transport equations for species mass fractions and sensible enthalpy, using Favre-filtered quantities, assuming a unity Lewis number. The dynamic Smagorinsky model is used to model turbulence with a variable Prandtl number formulation [1]. Combustion is considered to be infinitely fast and it is based on the Eddy Dissipation Concept (EDC) [1] calculating the fuel reaction rate as: where ߛ = ‫ܥ‬ ఊ ൫ߥߝ ௦௦ /݇ ௦௦ ଶ ൯ ଵ/ସ is the size of the fine structures, ߯ is the reacting part of the fine structures and ߬ = min (߬ , ߬ ௧௨ ) is a time scale which considers mixing under both laminar (߬ = ߂ ଶ ‫ܥ(/‬ ߥ)) and turbulent conditions (߬ ௧௨ = ‫ܥ‬ ఛ ൫ߥ/ߝ ௦௦ ൯ ଵ/ଶ ). The constants ‫ܥ‬ ఊ and ‫ܥ‬ ఛ are determined dynamically based on the local Reynolds and Damköhler number similar to the work by Parente et al. [3]. The Damköhler number is defined as ‫ܽܦ‬ = ߬/߬ with the chemical time scale given by an Arrhenius equation and the constants for methane combustion taken from [4]. For radiation modelling, the radiative transfer equation is solved using the finite volume discrete ordinates model with the absorption/emission term modelled through the weightedsum-of-gray-gases model (WSGGM). The path length, L, is calculated as ‫ܮ‬ = ‫ܣ/ܸ6.3‬ with the volume, V, calculated by summing all the cell volumes where reaction takes place. Assuming a rectangular shape (due to the geometry of the line burner), then the corresponding surface area, A, can be calculated. This approach bypasses any modelling uncertainties related to prescribing a fixed path length (which will exist as the flame extinction limit is approached) and allows it to vary dynamically during the numerical simulations. The modelling approach employed to account for local extinction considers an enthalpy balance based on the concept of a critical flame temperature. The criterion, examines whether a reactant mixture within a 3rd European Symposium on Fire Safety Science IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1107 (2018) 062003 doi :10.1088/1742-6596/1107/6/062003 computational cell will have enough energy to raise its temperature above a critical flame temperature to sustain combustion and is expressed as [5]: where ܻ is the mass fraction, ℎ is the sensible enthalpy, the sub-scripts F, O and D denote the fuel, oxidizer (i.e., air) and diluents (i.e., inert gases or products of combustion) present in the reacting mixture, ܶ ෨ is the (filtered) initial temperature cell temperature, ‫ܪ∆‬ ,ி is the heat of combustion of the fuel and ܶ ி் is a critical flame temperature. The composition of the reactant mixture within a computational cell that can react is calculated based on the composition of the species present in that cell as where ܻ ෨ ி , ܻ ෨ ை , ܻ ෨ are the (filtered) mass fractions of fuel, oxidizer and diluents in the computational cell and s is the oxidizer to fuel mass ratio. The enthalpy balance is based on the initial composition of the reactant mixture before any combustion is computed (i.e., after scalar transport but before any reaction). If the inequality is false then reaction within a computational cell is inhibited for the current time step. Based on the formulation of the model, any excess fuel within a cell acts as a diluent while any excess air does not. In general, the model is applicable to scenarios involving thermal quenching and quenching by dilution, since both mechanisms can be accounted for in the enthalpy balance.
For the evaluation of the flame extinction criterion, a critical flame temperature, ܶ ி் , is needed. A fixed value of 1700 K [6] is suggested by the CFD code Fire Dynamics Simulator (FDS). A sensitivity study with values between ܶ ி் = 1500 − 1700 K in [7] revealed some sensitivity in the numerical predictions when simulating the same test case. A different approach is employed here which allows the critical flame temperature, ܶ ி் , to vary depending on the local composition and residence time in every computational cell.
The extinction temperatures as a function of the residence times, W, are taken the from Perfectly Stirred Reactor simulations for different diluent agents (N 2 , CO 2 , H 2 O) with varying concentrations reported in [8].
The local critical flame temperature on every computational cell is calculated as a mass-weighted average of the critical flame temperature based on the local inert (N 2 , CO 2 ), given by a power law in the form of ܶ ி், (ଶ,େଶ) = 1407.4߬ ି.ଷ , and the critical flame temperature based on the non-inert diluents (H 2 O) given by a power law in the form of ܶ ி், (ୌଶ) = 1517.8߬ ି.ଷସ . The resulting critical flame temperatures are bounded between 1480 K and 1780 K, common minimum and maximum limits for most hydrocarbon fuels. In general, different approaches can be employed for determining ܶ ி் (e.g., from PSR simulations, opposed flow calculations with detailed chemistry, experiments) for different fuels when data are not available. Alternatively to what is employed in this study, ܶ ி் values as a function of strain rate can also be considered bringing this way some sub-grid scale effects in the flame extinction modelling. Nevertheless, these two approaches are not completely different from one another as small residence times in PSR calculations will result in ܶ ி் values similar to what you would get from e.g. opposed flow calculations with large strain rates and the other way around. Flame re-ignition is also an important aspect when it comes modelling flame extinction with infinitely fast chemistry. At any point in space and time, any unburned fuel in the computational domain can mix with the oxidizer and react regardless of the temperature of the mixture. This behavior is not realistic and can cause erroneous phenomena of re-ignition further downstream from the fire source where they should not be present. For this reason, an ignition temperature, ܶ , is specified in the numerical simulations which controls when any re-ignition phenomena can occur. This criterion essentially only allows any unburned fuel to react with the oxidizer if the temperature in the computational cell is above the specified ܶ value (and the enthalpy balance in Eq. (2) is satisfied). However, this approach cannot be globally applied within all the computational domain without affecting the primary ignition and combustion process. Following the work performed by White et al. [6], a small pilot flame area just above the fuel inlet, 25 mm in height, is specified in which combustion is allowed to take place without allowing for any flame extinction or re-ignition to occur (i.e., both phenomena essentially cannot occur there). Inside this region it is ܶ = 0 K. For the region outside the pilot flame area, ܶ = 1200 K is used, taken from extrapolated values of ignition temperatures of methane combustion at low strain rates, typically found in fires, based on [9]. This value is also close to the range of re-ignition temperatures suggested in [10] for the same test case. It is worth noting that this value is above the auto-ignition temperature of methane (i.e., 910 K). The use of this value is motivated by the fact that any re-ignition phenomena will occur for a lean methane mixture, diluted by hot air and combustion products, which will raise the ignition temperature of the mixture above the auto ignition temperature of methane. Nevertheless, a more advanced re-ignition modelling approach, that is not based on a user input parameter, is worth considering in future numerical simulations. A possible approach could be the use of a re-ignition temperature as a function of strain rate for different methane dilution ratios based on [9].

NUMERICAL SETUP
A computational domain of 1.6 m x 2 m x 1 m is used. The base mesh consists of 25 mm cells stretched towards the top boundary with an expansion ratio of 1.5. A grid refinement strategy is then employed up to the experimentally reported flame height. A box having dimension of 0.8 m x 0.8 m x 0.6 m refines the mesh to 12.5 mm while a second box of dimensions 0.6 m x 0.6 m x 0.4 m refines the mesh to 6.25 mm. The total number of cells is then approximately 886000. The methane flow rate is set to 1 g/s while the co-flowing oxidizer is provided at 85 g/s, according to the experiments. For the cases involving flame extinction, the oxygen mass fraction in the oxidizer port is initially maintained at ܻ ைଶ = 0.233 for 4 s and is subsequently decreased by 0.005/s (followed by equivalent increase in ܻ ேଶ ). For the solution of the RTE, 48 solid angles are used for angular discretization. All simulations were set to run for 40 sec with a varying time step, limited by a maximum Courant number of 0.9. Time averaging over the last 35 s was used to produce time-averaged quantities. A second order backward scheme is used for time discretization. A second order filtered linear scheme is used for the convective terms which introduces small amounts of upwind if the solution is unbounded. For scalar transport, a second order TVD scheme using a Sweby limiter is applied.

RESULTS
The centerline mean and rms temperatures for different mesh sizes are first presented in Fig. 1(a) for a test case with ܺ ைଶ = 0.18 in the oxidizer stream. The extinction model was not employed as the experimentally measured combustion efficiency was about 99%. Maximum flame temperatures of approximately 1240±40 K are predicted with all three mesh sizes, showing reasonable grid convergence, with an increase in the rms temperatures as the grid size decreases. The slightly lower maximum flame temperatures with the finer grid are attributed to a combined effect of turbulence, combustion and radiation. Even though the integrated reaction rates for all cases were the same (i.e., 1 g/s), the reaction rates right at the fuel source decrease and shift slightly downstream with decreasing grid size. Additionally, the predicted radiative fractions for the 25, 12.5 and 6.25 mm cases were 0.172, 0.196 and 0.212, respectively, while the experimental was about 0.195. The lack of influence of grid resolution on the mean and rms temperature profiles is attributed to the use of dynamic models that are able to self adjust their coefficients depending the local properties of the flow. Nevertheless, further investigation is needed in order to make more concrete statements regarding this aspect. The radial profiles of the mean temperatures at two different heights for the ܺ ைଶ = 0.18 case are presented in Fig. 2. Overall, the numerical predictions with the finest grid size employed (i.e. 6.25 mm) are able to reproduce quite reasonably both the flame width and the maximum flame temperatures at these two heights.
The predicted combustion efficiencies from the numerical simulations are presented in Fig. 3(a). It is interesting to note that the numerical predictions with the coarsest grid (i.e., 25 mm) are strongly affected by the extinction/re-ignition model and cannot accurately simulate the scenario. As the grid size is refined this influence diminishes with minimal influence on the 12.5 mm grid size (i.e., combustion efficiency slightly lower than 1 in the initial stage when ܺ ைଶ = 0.21 in the oxidizer stream). Nevertheless, the predictions with both the 12.5 mm and 6.25 mm grid sizes are able to predict the oxygen concentration where extinction takes place (i.e., approximately ܺ ைଶ = 0.122 based on the experimental data). The predictions with the 6.25 mm grid are in better agreement with the experiments, when compared to the 12.5 cm grid case, and are able to reproduce the experimental trends towards flame extinction even though the decay phase somewhat starts faster. Based on the results shown above, it is clear that with the current modelling approach a relatively fine grid is needed in the numerical simulations in order to accurately predict flame extinction. This can perhaps pose some limitations in cases where large-scale fire scenarios need to be simulated. It is also worth noting that the re-ignition temperature, ܶ , can influence the predictions and the resulting combustion efficiency (not shown here). Effectively an increase in the ܶ values results in flame extinction earlier than expected (i.e., for concentrations ܺ ைଶ > 0.12). The start of the decay phase of the combustion efficiency is also partially affected by the choice of this value. Equally important can be the choice for the critical flame temperature, ܶ ி் , particularly in the case where a constant value is chosen and does not vary locally in every computational cell. The combination of different constant ܶ ி் and ܶ values can lead to equally accurate predictions of flame extinctions for perhaps the wrong reasons (i.e., compensating effects). It is obvious that accurate predictions of flame extinction with infinitely fast chemistry are directly linked to the combined modelling of combustion and radiation since it is the resulting (filtered) temperatures that will determine the degree of flame extinction (through comparison with ܶ ி் in the enthalpy balance and with ܶ in order to determine re-ignition). The influence of radiation/turbulence interaction should also be explored in the future.
Similar observations can be made by examining the predicted radiative fractions, ߯ , shown in Fig 3(b). Reasonably good predictions are obtained for the 12.5 mm and 6.25 mm cases with the simulations being able to reproduce the decrease in the heat released due to radiation as the flame extinction limit is approached. However, some over-predictions on the resulting radiative fractions on the finest grid case are still evident. It is worth noting that capturing this decrease in ߯ can be very important when trying to accurately predict flame extinction in any random fire scenario. For example, the use of a constant radiative approach cannot be directly employed, since it will lead to flame extinction faster than expected, nor can an a-priori tabulation of the radiative fractions with decreasing ܺ ைଶ values be expected for any given scenario. Additionally, the use of the simpler gray-gas model is known to over-predict the resulting radiative fractions [11] unless a calibration constant is introduced which also puts some limitations in the application of the model to random scenarios. The accuracy of the WSGGM with a dynamically determined beam length in this scenario was satisfactory in terms of predicting the resulting radiative fractions. However, application of the model in more scenarios involving different fuels is needed before making concrete statements regarding its accuracy. Equally important, the use of a dynamically determined beam length will be more challenging when it comes to scenarios involving flame spread where a distinct flame shape cannot be so easily determined.   [2]. Use of advanced modelling approaches related to turbulence, combustion and radiation was made along with the application of an enthalpy-based extinction model with a locally varying critical flame temperature. The predictive capabilities of the models was analyzed by comparing the numerical predictions for different grid sizes to the experimental thermocouple temperatures, combustion efficiencies and radiative fractions.
The numerical predictions were, overall, in good agreement with the experimental data when fine grid sizes were used (i.e., in the order of 6.25 mm). More specifically, the resulting maximum flame temperatures and flame widths matched reasonably well the experimental profiles at two different heights examined. The numerical simulations were also able to reproduce the decreasing trend of the combustion efficiencies and radiative fractions with some discrepancies when the extinction limit was reached.
From the results of the current work, it is suggested that caution must be taken when modelling flame extinction with infinitely fast chemistry using an enthalpy-based extinction model based on the concept of a critical flame temperature. Accurate modelling of combustion (i.e., in terms of resulting mean and rms flame temperatures) and radiation (i.e., in terms of decreasing radiative fractions), along with relatively small grid sizes, is a pre-requisite for predicting flame extinction. Re-ignition modelling can also be important. The choice of a constant re-ignition temperature approach might be worth revisiting in the future by employing a more advanced re-ignition model based e.g. the local composition and strain.