A Numerical Study on Fire Development in a Confined Space Leading to Backdraft Phenomenon

Abstract

This paper presents the results of numerical experiments on fire development and the backdraft phenomenon. The numerical model of fire development built with the use of Ansys Fluent was then validated based on literature data taken from real fire experiments. Some theoretical foundations of airflow and combustion modelling were added. Some features of the numerical model, which allowed for its high accuracy to be achieved, were widely discussed. Since large buoyancy forces were involved, to reproduce the decrease of the atmospheric pressure with height, a variable static pressure was applied using UDF (user-defined function).The results showed good accordance taking into account both the temperature profiles and the distribution of the airflow velocity. Once the model was validated, the research was extended to examine the backdraft phenomenon. The results revealed characteristic phases of the phenomenon and the occurrence of the gravity current as well, which were reported by empirical experiments.

1. Introduction

Due to a tight shell there are specific conditions for fire development in modern buildings. The way the fire develops in a compartment depends on the oxygen supply. If there is an unlimited oxygen supply (open window or door) and sufficient amount of fuel, the fire develops freely until it reaches the fully developed stage [1]. If the air supply is limited, the fire will become extinct. Products of incomplete combustion will begin to accumulate in the compartment. If there is still no fresh air inflow, the fire will dim, the temperature inside will drop, and the room will be filled with thick smoke [2,3]. Such fire conditions are referred to asunder-ventilated. Under-ventilated fires often occur in buildings where the free air supply is very limited. Fire development in a confined space can be analyzed from different points of view, and many experimental and numerical studies have dealt with this issue [4,5,6]. Some phenomena related to under-ventilated fires are dangerous and pose a great threat to rescue teams. First of all, it concerns the phenomenon of backdraft.

Backdraft can occur when a stream of fresh air is suddenly brought into a compartment where an under-ventilated fire develops. Flammable fire gases with unburned fuel mixed up with the inflowing air start to burn and strengthen the fire [2]. The flames move outside towards the inflowing air. As a result, the flames travel outside through the upper part of the newly formed opening andfresh air inflows through its lower part, forming a so-called gravity current.

Studies on this extremely dangerous backdraft phenomenon have been carried out for a long time. Fleischmann examined backdraft both by model tests and numerical research as well [7,8,9]. He found that the mass fraction of the unburned fuel in a compartment greater than 15% was the necessary condition for backdraft occurrence. Some later research showed that backdraft occurrence was also possible for lower values of the fuel mass fraction. Weng carried out scaled model experiments and confirmed that the mass fraction of unburned fuel was the crucial parameter in determining backdraft occurrence. He determined the threshold value at the level of 9.8% [10]. Next, he examined the possibility of mitigating backdraft effects through the application of water mist. He tried to prove that the water mist weakened the backdraft process because it caused a mass fraction decrease, which is the crucial parameter influencing this phenomenon [11]. Gottuk made similar conclusions about the influence of water injection on weakening the backdraft effects. He claimed that the injection of water spray mainly resulted in a decrease of the fuel mass fraction but not in a temperature decrease. However, he reported that the threshold value of the fuel mass fraction necessary for the occurrence of backdraft should be at least 16% [12].

Very interesting scaled model tests were carried out at the BRE (Building Research Establishment) Center for Fire Safety Engineering. They examined the flashover phenomena as well. Flashover consists in a sudden self-re-ignition of a dimmed under-ventilated fire when fresh air starts to be supplied. The conclusions underlined the great importance of the time of fire development until it becomes under-ventilated and the temperature distribution in the compartment volume. The authors argued that the temperature inside is crucial to the occurrence of backdraft: It is possible only if the temperature inside exceeds 350°C, but if a spark of ignition appears, it can even happen at 300°C. The experiments showed that the occurrence of backdraft is possible even in the early phase of fire development, at a lower temperature than that required for flashover occurrence [13]. Polypropylene as a fuel was used in these experiments instead of the usual methane.

Fuel other than methane was examined by Chen as well. He estimated the possibility of backdraft occurrence for a wood fire [14]. The experiments carried out at Lund University focused on backdraft and provided a lot of information [15] and therefore many scientists refer to the results of these experiments. However, their application to the validation of numerical models of backdraft is limited as the experiments were carried out in a container located outdoors. The authors highlighted the significant impact of the natural wind on backdraft and the gravity current velocity.

As it was already mentioned, air inflow into the compartment is necessary for the occurrence of backdraft. Many researchers argued that just the gravity current itself is of crucial importance here, so it has been examined many times [15,16,17,18].

Since the backdraft phenomenon develops dynamically and, in addition, it is very dangerous, numerical analyses is willingly used to study it. Mowrer carried out some theoretical considerations and built a numerical model of an under-ventilated fire [19]. Horvat numerically analyzed backdraft using Ansys CFX [20]. The DES (detached eddy simulation) turbulence model and the EDM (eddy dissipation model) combustion model were also used. His simulations covered the moment of opening the door to the room in which the fire was developing. In the simulation results, both the backdraft and gravity current were visible. The model showed that the velocity of flame propagation reached 40 m/s. Fire dynamics simulator (FDS) software has also been used for numerical studies on backdraft. Ferraris proposed the application of two combustion regimes: Premixed (fuel and oxidant already mixed and then ignition takes place) and non-premixed (fuel and oxidant enter the combustion area independently) [21,22]. Additionally, the authors distinguished five phases of backdraft:

• A developing fire must reach the necessary condition;
• Ignition and a very short period of spherical propagation;
• Forming of a plain flame front;
• Stretching of the flame front through the opening; and
• Fireball outside the compartment.

CFD methods are commonly applied in fire engineering. Since people’s safety often depends on such obtained results, the reliability of numerical models is of great importance. Therefore, the use experimental data to validate numerical models is desirable. Therefore, the validation of numerical models based on experimental data is necessary as it should confirm the adopted assumptions and simplifications as well as prove that the models describe real physical phenomena [4,23,24].

The objective of the presented work was a numerical experiment on the fire development in a confined space, in which under-ventilated conditions progress rapidly and may result in backdraft. The main contribution of the work is a detailed study on the numerical model features, which allows for reliable modeling of the examined phenomena. The paper includes a description of the principles of fluid flow modeling and briefly discusses combustion modeling, which proves important as some model details refer directly to the physical foundations of fluid flow modeling. It also includes a description of the real fire experiment carried out by Betting et al. [25] and a detailed description of the built-up numerical model, with justification of some of the adopted assumptions. The results, split into two parts, model validation and the backdraft phenomenon analysis, are presented in the third section, followed by the conclusions.

2. Principles of Fluid Flow Modeling

Since the fluid flow and its modeling are very complex, it is impossible to present the problem in brief. Therefore, some of the basic notions and matters mentioned give a better understanding of the case. Nowadays, there are numerous software packages that allow for modeling the different phenomena of fluid flow. Yet, the physical basics of all of them remain the same: The differential equations governing the fluid flow are solved numerically. Ansys Fluent is one of the best recognized software tools, which is widely applied in many fields of technology and science. It solves the flow problem with the finite volumes method [26]. The examined domain is divided into a large number of small volumes called cells. For each cell, the balance equations for mass (Equation (1)), momentum (Equation (2)), and energy (Equation (3)) are solved, which allows the velocity field $\overrightarrow{u}$ [27] to be obtained:

$$\frac{\partial }{\partial t}{{\displaystyle \int }}_V\rho dV+{{\displaystyle \oint }}_S\rho \overrightarrow{u}·d\overrightarrow{S}={{\displaystyle \int }}_V{S}_{m}dV,$$

(1)

$$\frac{\partial }{\partial t}{{\displaystyle \int }}_V\rho \overrightarrow{u}dV+{{\displaystyle \oint }}_S\rho \overrightarrow{u}\overrightarrow{u}·d\overrightarrow{S}={{\displaystyle \int }}_V\left(\rho \overrightarrow{F_m}+\nabla \stackrel{=}{T}\right)dV,$$

(2)

$$\frac{\partial }{\partial t}{{\displaystyle \int }}_V\rho edV={{\displaystyle \int }}_V\left(\rho \overrightarrow{F_m}·\overrightarrow{u}+S_e+\nabla \left(\stackrel{=}{T}·\overrightarrow{u}+\overrightarrow{q_e}\right)\right)dV,$$

(3)

where V denotes the volume and S denotes the surrounding surface. $S_m$, $\overrightarrow{F_m}$, $\stackrel{=}{T}$, and $\overrightarrow{q_e}$ are the mass source efficiency, unit mass force, stress tensor, and heat flux, respectively. Quantity $e=e_i+\frac{1}{2}\overrightarrow{u}\overrightarrow{u}$ expresses the specific energy as a sum of the internal specific energy and kinetic energy.

Each equation includes the change in time of the examined quantity, its convectional and diffusion transport, and its possible generation. If a considered case involves some other phenomena, similar equations must be solved for additional quantities specific for the case. For example, if a combustion process is modeled, the balance equations of the mass fraction of the final and intermediate combustion products are solved.

The next group of equations, which constitute a numerical model, are the equations that describe the boundary and initial conditions. Finally, the last group contains equations describing the material properties and adopted simplifications and approximations. Models of turbulence or combustion belong to this group.

2.1. Turbulence Modeling

Only in rare cases of very low fluid velocity is the flow laminar—it occurs in non-disturbed layers of different velocity. More often, turbulence appears; it is defined as a time and spatial random variability of the quantities describing the flow. The nature of the flow depends on the value of the Reynolds number, which determines the ratio of inertia forces to viscosity forces:

$$Re=\frac{\rho uL}{\mu },$$

(4)

where ρ, μ, and u denote respectively the density, dynamic viscosity, and velocity of the fluid. For Re < 2300, the flow can regarded as laminar; for greater Re values, eddies start to appear; and at Re > 5000, the flow is always fully turbulent. It becomes entirely chaotic and its kinetic energy is transferred into the internal energy (heat) of the fluid in a process called a Kolmogorov energy cascade:

• The biggest eddies contain the most kinetic energy coming directly from the kinetic energy of the flow;
• Energy is gradually transferred towards smaller eddies;
• At the end of the cascade, for the smallest eddies (length of order of 0.0001 m), the viscous processes become the most important factor, and finally, the energy is dissipated;
• The structure of the biggest eddies is highly anisotropic and depends mainly on the geometry of the whole system. Processes at the smallest scale are regarded as isotropic and depend mainly on the fluid properties [27].

The flow of any Newtonian fluid (isotropic and the components of the stress tensor $\stackrel{=}{T}$ are linearly dependent on the strains) at any scale regarding the density variability is governed by the Navier–Stokes equation (µ_b denotes the volume viscosity):

$$\frac{\partial \left(\rho \overrightarrow{u}\right)}{\partial t}+\nabla ·\left(\rho \overrightarrow{u}\overrightarrow{u}\right)+\nabla \left(p-{\mu }_b\nabla ·\overrightarrow{u}\right)-\mu {\nabla }^2\overrightarrow{u}=0.$$

(5)

The unknowns in this equation are the three components of the velocity vector, pressure, and fluid density. As it can be noticed, nowadays, it is practically impossible to directly solve the problem of fluid flow. This is because the size of the cells of a computational domain should be the same scale as the smallest eddies, which leads to extremely large numerical problems. Therefore, some approaches allowing for a decrease of computational complexity have been introduced. The two main ideas are large eddy simulation (LES) and Reynolds-averaged Navier–Stokes (RANS) [28]. The LES approach assumes the direct numerical solution for the biggest eddies and analytical modeling of the smaller ones. It is assumed that the analytically modeled eddies should contain up to 20% of the turbulent energy [29]. The RANS approach origins form the observation that from the engineering point of view, full knowledge on all flow details is unnecessary; only the averaged values depicting the flow are of practical meaning. The main assumption for RANS models is the so-called Reynolds approximation, which states that the instantaneous fluid velocity is a sum of the average velocity and fluctuating velocity:

$$u_i={\overline{U}}_i+u_i^{\prime }i=x,y,z.$$

(6)

By definition, the average of the fluctuating component is equal to zero:

$${\overline{u}}_i^{\prime }=0.$$

(7)

After inserting these expressions into the original equations describing the fluid flow, one can obtain another equations system, which allows the sought averaged velocities to be determined, and a new unknown quantity appears. It is the Reynold’s tensor of viscous stresses. As a symmetric tensor, it contains six variables. Such equations system can be directly solved (Reynold’s stress model, RSM), but six additional equations for the tensor’s components must be added as well as the one for turbulence kinetic energy (k), which leads to high requirements for computational resources. Despite this drawback, the RSM model is willingly used for complex anisotropic flows with violent changes of the fluid velocity. It is the most accurate approach in the RANS family [30].

Further simplification leads to models that are commonly used in science and technology practice. According to the Bussinesq’s hypothesis, the components of the tensor of viscous stresses can be expressed by a new quantity, which is called turbulent viscosity [31]. Generally, in this approach, the energy transport due to turbulence is described in a similar way to molecular transport. Depending on some details, two families of models were developed: k-ε, where turbulent viscosity is determined by the turbulence kinetic energy (k) and the dissipation rate of turbulence energy (ε); and k-ω, where turbulent viscosity is determined by the turbulence kinetic energy and turbulence frequency (ω).

Although k-ε models behave better for bulk fluid, k-ω models better perform for flows near walls. Nowadays, two enhanced models are widely applied: k-ε realizable and k-ω SST, which both give reliable results for most cases in an acceptable time [32]. Unfortunately, their common drawback is a low accuracy for complex anisotropic flows [27]. Some other approaches have been developed, for instance, detached eddy simulation (DES), which combines the LES method for large eddies and RANS method for wall flows.

2.2. Combustion Modeling

The process of combustion is a very complex phenomenon. It involves many chemical reactions among the components of the fuel and the oxidant. It produces various intermediate and final products and releases large amounts of energy as well. Additionally, it is usually accompanied by turbulent flows. This is why some simplifications and pre-assumptions are needed to describe it. There are four combustion models particularly useful for fire simulation, provided by Ansys Fluent software [26]:

• Species transport: A combustion process is modeled as a source, which emits predefined compounds according to the assumed chemical reaction. The mass flow rate and temperature of the emitted species should be adjusted to fit the required fire power. In this case, additional transport equations for the mass fraction of each emitted compound (Y_i)must be solved:

$$\frac{\partial \rho Y_i}{\partial t}+\nabla \left(\rho \overrightarrow{u}{Y}_i\right)=-\nabla \overrightarrow{J_i}+S_i.$$

(8)

The term $\overrightarrow{J_i}$ expresses the mass diffusion stream of the i-th component. Since this model by default ignores the details of combustion chemistry, it is suitable for fires controlled by a fuel supply [23]. Optionally, the feature of taking into account chemical reactions can be switched on. This leads to an eddy dissipation model (EDM), which is quick, but the reactions can be at most two steps.

• Premixed combustion: Fuel and oxidizer are mixed at the molecular level prior to the combustion, which occurs as a flame front propagating into the un-burnt reactants. The mixture is not in the equilibrium state. The main concept is the reaction progress variable c, defined as:

$$c=\frac{\sum Y_i}{\sum Y_{i,ad}},$$

(9)

where Y_i,ad denotes the mass fraction of the i-th compound after complete adiabatic combustion. Thus c expresses the progress of the combustion process and it takes a value of 0 for an un-burnt mixture and 1 for a completely burnt mixture. The propagation of the flame front is modeled by solving a transport equation for c:

$$\frac{\partial \rho c}{\partial t}+\nabla ·\left(\rho \overrightarrow{u}c\right)=\nabla \left(\frac{{\mu }_T}{S{c}_T}\nabla c\right)+\rho S_c,$$

(10)

where the terms µ_T, Sc_T, and S_c denote the turbulent viscosity, turbulent Schmidt number, and sources of the reaction progress, respectively. In this approach, the flows are “faster” than chemical reactions [33].

• Non-premixed combustion: The final and intermediate combustion products and released energy depend only on the local composition of the gas mixture and the temperature [34,35]. The mixture reaches equilibrium or is close to it, which can be expressed as “what is mixed is burnt”. This approach assumes that the chemical reactions of combustion are regarded as immediate in comparison with the flows. A simplifying assumption adopted is that a quantity called the mean mixture fraction is kept for all elements. This quantity expresses the mass fraction of the considered elements, which comes from fuel:

$$f_i=\frac{Z_i-Z_{i,ox}}{Z_{i,fuel}-Z_{i,ox}},$$

(11)

where the symbol Z_i denotes the mass fraction of a given element; the indexes describe its origin (fuel or oxidant). If a flow is turbulent, the diffusion coefficients for all compounds are equal and then the values of f_i are the same for each element and they can be replaced by a single value f. Therefore, solving one additional equation describing the transport of the mean mixture fraction $\overline{f}$ is adequate:

$$\frac{\partial (\rho \overline{f)}}{\partial t}+\nabla ·\left(\rho \overrightarrow{u}\overline{f}\right)=\nabla ·\left(\left(\frac{\kappa }{C_p}+\frac{{\mu }_T}{{\sigma }_T}\right)\nabla \overline{f}\right)+S_{\overline{f}},$$

(12)

where the terms 𝜅, C_p, and σ_T denote the laminar thermal conductivity of the mixture, mixture specific heat, and Prandtl number, respectively; and $S_{\overline{f}}$ denotes any sources. Since this approach includes combustion chemistry, it is suitable for fires controlled by ventilation. The possible states of the mixture can be described in terms of the probability density function (PDF), which expresses the fraction of the time that the fluid spends at each state of the species, temperature, and pressure. To speed up the calculation, some preliminary calculations are done and a so-called PDF table is created. It is a loop-up table, which contains the discrete values of PDF and allows for quick access to the chemistry details depending on the actual mixture state [36].

• Partially premixed combustion is a combination of the premixed and non-premixed approach. It can be used to model partially premixed flames ranging from the non-premixed state to the completely premixed mixture. For the first extreme state, the reaction progress variable (c) is equal to 1; for the latter one, the mixture fraction (f) is constant in the domain. In general, the transport equations for the progress variable (Equation (10)) and mixture fraction (Equation (12)) are solved. The average values ($\overline{\Phi }$) of different scalars ($\Phi$) depicting any mixture state are calculated using the probability given by PDF:

  $$\overline{\Phi }={{\displaystyle \int }}_0^1{{\displaystyle \int }}_0^1\Phi \left(f,c\right)p\left(f,c\right)df dc.$$

  (13)

Both non-premixed and partially premixed models can be extended by introducing the idea of a flamelet. This concept postulates that a turbulent flame is an ensemble of thin laminar flames called flamelets. A flamelet has an internal structure that is not significantly altered by turbulence. The basic flamelet model assumes two counterflow jets of a fuel and oxidizer. A new quantity called scalar dissipation is introduced (D is a diffusion coefficient):

$$\chi =2D\left|\nabla f\right|.$$

(14)

The chemistry calculations can be simplified to one dimension along the model axis. They are completely determined by two variables f and χ and can additionally be pre-processed. The flamelet approach allows for the introduction of a detailed combustion chemistry mechanism for all involved reactions.

2.3. Radiation Modeling

An important mode of heat transport is the radiation phenomenon. A full analysis of the radiative heat transport should take into account the surface emission and absorption and the interaction with the fluid, which involves refraction, scattering, emission, and absorption. The fact that radiation of various wavelengths can behave differently should also be considered. The overall equation governing the radiative transport (RTE) for an absorbing, emitting, and scattering medium at position $\overrightarrow{r}$ in the direction $\overrightarrow{s}$ for the radiation of intensity I passing a layer that is ds thick is as follows [37]:

$$\frac{dI\left(\overrightarrow{r},\overrightarrow{s}\right)}{ds}+\left(a+{\sigma }_s\right)I\left(\overrightarrow{r},\overrightarrow{s}\right)=a{n}^2\frac{{\sigma }_B{T}^4}{\mathsf{\pi }}+\frac{{\sigma }_s}{4\mathsf{\pi }}{{\displaystyle \int }}_0^{4\pi }I\left(\overrightarrow{r},{\overrightarrow{s}}^{\prime }\right)\Phi \left(\overrightarrow{s},{\overrightarrow{s}}^{\prime }\right)d{\Omega }^{\prime },$$

(15)

where the terms a, σ_s, n, and σ_B denote the absorption coefficient, scattering coefficient, refractive index, and Stefan–Boltzmann constant, respectively. The phase function Φ describes the scattering of radiation incoming at direction $\overrightarrow{s}$ into direction ${\overrightarrow{s}}^{\prime }$; the scattered radiation is integrated over the full solid angle Ω’. This function depends on the fluid nature and its properties.

Depending on the additional assumptions and simplifications, there are four main radiation models available in Ansys Fluent [38]. The main criterion when selecting the proper radiation model is the optical density of the fluid, which is a product of the absorption coefficient and typical beam length (a∙L).

• P1—for optically dense fluids; takes into account the absorption by the fluid (a∙L> 1). Any wall is assumed to be a diffuse gray surface.
• S2S—describing only the radiative energy exchange between solid surfaces, and omitting the absorption (a∙L= 0). At the initial stage of the simulation, the so-called view factors are calculated. These variables keep the data on the surfaces’ arrangement and are then used to speed up the calculations.
• DTRM—radiation emitted by a surface element in a certain range of solid angles is approximated by a single ray. Then, the history of each ray that passes through the fluid volume is traced (any value of a∙L). The ray paths are calculated and stored prior to the fluid flow calculations. An important drawback is that it cannot be used in the parallel calculation mode.
• DO—instead of single-ray tracing, the full equation governing the radiation is solved for each finite volume cell. Since the model takes into account all processes, it is the most accurate but also the most computationally complex (any value of a∙L).

3. Description of the Modeled Real Experiment

The experiment was carried out by Betting et al. [25]. It was a series of test fires in a confined space with a limited fresh air supply. The test chamber was composed of a standard maritime container. The walls and the ceiling were insulated by a thick layer of mineral wool. The floor was made of concrete. A grid of 36 propane burners was mounted at one side of the container and a roof opening was placed on the opposite side. The burner’s grid formed a square with an edge of 1 m, mounted at a height of 0.4 m above the floor. The air velocity distribution was monitored using particle image velocimetry (PIV) [39]. For this purpose, alumina particles with a diameter of 50 μm were injected through an air inlet positioned below the burner. The additional amount of air was as low as 10% of the stoichiometric amount needed for combustion and did not disturb the overall air movement in the chamber. A rack of 10 thermocouples type-K was mounted at the centerline of the opening, 0.1 m below it. The thermocouples’ spacing was 0.1 m. A simplified sketch of the experimental setup is shown in Figure 1.

As the roof opening was the only connection with the environment, it was the only possible source of fresh air (besides the small amount of air for the alumina particles’ injection). Therefore, depending on the adopted fire power, it allowed for testing of both well-ventilated and severely under-ventilated conditions.

4. Numerical Model

The numerical model was prepared using Ansys Fluent software and it accurately reproduced the experimental set-up. The size of the container, locations of the equipment, and the structure and properties of the walls, ceiling, and floor were exactly mapped. The roof opening was modeled as a ‘pressure outlet’. Heat losses through the container walls, roof, and floor were modeled using the Ansys Fluent feature ‘shell conduction’. It allowed for analytical modeling of the heat flow in bulkheads without an extension of the computational domain.

A non-premixed combustion model was selected for the first phase of the simulation (fire development) as the stream of propane was introduced into a volume filled by air. As the fire developed, oxygen burnt out, which might lead to severe under-ventilated conditions for high-power fires. Since the backdraft phenomenon is a case of deflagration—a quick combustion—the equilibrium assumptions for non-premixed combustion might not be met. This was the reason why partially premixed combustion was checked for this phase as well.

The burner consisted of the cuboid space of real dimensions to which a ‘source term’ property was assigned. It allowed the generation of a suitable amount of propane in this space to achieve the desired heat release rate (HRR). The generated amount of propane brought the relevant upward momentum as well, which imitated the operation of a real burner.

Since our intension was to examine the phenomenon of backdraft, our numerical research went further than the real experiment. An additional door was modeled in the gable wall opposite the burner and a room adjacent to the door was added. The wall of the added room, opposite the door, was modeled as the ‘pressure outlet’. Initially, the door was closed. Then, after some time, when an amount of unburned propane and intermediate flammable combustion products had accumulated in the container, the door was opened. It was implemented by changing the ‘wall’ boundary condition into ‘interior’ for the door surface. A sketch of the numerical model is shown in Figure 2.

The modeled phenomenon involves buoyancy forces; therefore, the value of static pressure at both boundary conditions of the ‘pressure outlet’ type was controlled by the user-defined function (UDF) to reproduce the dependence of the atmospheric pressure on the height:

$$p\left(z\right)=p_0 – z\rho g.$$

(16)

A UDF is a plug-in, which can be used for the free modeling of many features of a numerical model built using Ansys Fluent. For instance, UDFs can alter the standard boundary conditions or material properties. The code of the used UDF is presented below (Box 1). It is rather simple and contains a loop, which iterates over a list of all faces for a given boundary and calculates the atmospheric pressure for each face according to Equation(14).

Box 1. The code of the applied UDF.

#include "udf.h"

#define GACC 9.81

 

DEFINE_PROFILE(phstat, thread, i)

{ face_t f;

 real z, vec[ND_ND];

 

begin_f_loop(f, thread)

 { F_CENTROID(vec, f ,thread);

 z = vec[2];

 F_PROFILE (f, thread, i) = -C_R(f, thread) * GACC* z;

 }end_f_loop(f, thread)

}

The numerical mesh was un-structural—it contained mainly tetrahedron elements. The basic size of the mesh cells was 0.2 m, but the mesh was significantly condensed in areas, where complex flows were expected. The volume of the burner was divided into rectangular elements of a size of 0.05 m, and the neighborhood of the alumina particles’ inlet was divided in the same manner. The size of cells in the entire container volume, which contained the burner, was 0.07 m. Additionally, the space close to the roof opening and the door was divided into cells of 0.1 m in size. The total number of cells was about 418,000. Such a meshing scheme was used on the basis of some preliminary calculations, which were carried out for different meshes.

The preliminary calculations showed a highly instable flow through the roof opening. At the beginning of the simulation, the air and the combustion products flowed out through the entire cross-section of the opening; later, both the outflow and inflow of fresh air were observed at the same time. It appeared that both the k-ε and k-ω models led to the results, which were a bit inconsistent with the data reported by Betting [24]. The application of the RSM model allowed sufficient accuracy to be obtained. Furthermore, it was especially important in the modeling of the backdraft stage, when complex anisotropic flows occurred. The preliminary calculations also proved that the generated soot amount was negligible, thus the fluid in the container was transparent at the beginning of the simulation. Next, when three-atomic molecules were generated in the combustion process (mainly carbon dioxide and water vapor), absorption and scattering became important. Since the calculations were in the parallel mode, the discrete ordinate (DO) radiation model was selected, despite its complexity. The absorption coefficient was calculated automatically according to the weighted-sum-of-gray-gases model (WSGGM).

The last task was to adjust the time step, which had to meet the Courant–Friedrichs–Lewy criterion, which shows the relation among the characteristic velocity of the phenomena, time step, and typical element size (C is Courant number, which must not be significantly larger than one):

$$\frac{u\Delta t}{\Delta x}\le C.$$

(17)

Since the course of phenomena for the examined case varied significantly during the simulation’s progress, it was difficult to adopt a single value for the time step duration. Therefore, the feature of an adaptive time step was enabled, which allowed Ansys Fluent to select the most suitable time step for the current state of the process. A summary of the numerical model details is presented in

Table 1. Numerical model summary.

Feature		Value
Turbulence model		RSM
Fluid material		Mixture of air and combustion products—ideal gas
Combustion model	burning phase	Non-premixed
	backdraft phase	Non-premixed/partially premixed
Fuel		Propane (C3H8)
Radiation model		DO
Operating pressure		1025 hPa
Gravitational acceleration		9.81 m/s2
Solver		Pressure based
Pressure spatial discretization		Body force weighted
Calculation mode		Transient
Time step	burning phase	Adaptive 0.001 s—0.2 s
	backdraft phase	Adaptive 0.00001s—0.02 s
Pressure/velocity coupling		SIMPLE
Under-relaxation factors	Pressure	0.5
	Momentum	0.5
	Energy	0.95
	Species	0.8
	Others	Default
Heat flow in bulkheads		Shell conduction

.

5. Results

The first step was the validation of the built-up numerical model. Therefore, the fire power was set to 0.2 MW, which corresponds to a propane mass flow rate equal to 0.0043 kg/s. The fire development can be divided into two phases. During the first phase, the air inside the container was warmed up and underwent thermal expansion. The buoyancy forces began to force air movement. The roof opening acted only as an outlet. In the second phase, when the air movement in the container was fully developed, a certain amount of fresh air was sucked in. In this way, a part of the roof opening still acted as an outlet, but the rest of it played the role of an inlet. Such a velocity pattern remained almost stable to the end of the simulation. Figure 3 shows the velocity vectors at the roof opening in both phases.

Figure 4 shows the temperature profiles for the experimental and numerical results. The temperature was recorded by the set of thermocouples mounted just beneath the roof opening (see Figure 1). The experimental values were averaged in time, and the numerical ones were the instantaneous values recorded every one second.

As can be observed, the accordance between both curves is acceptable despite some fluctuations. The discrepancy is especially visible in the transient region in which the air started to be sucked into the container (about the 30th second). Significant flow direction and velocity fluctuations occurred in this region, which influenced the recorded temperature. The fresh air supply mitigated the inside temperature increase, thus a flattening of the temperature profile can be observed. The temperature distribution in Figure 5 is in the form of a temperature field at the vertical symmetry plane caught at the 300th second of the fire development. The influence of the inflowing cold air is clearly visible in the region below the roof opening: The temperature is much lower there.

It is also possible to compare the velocity fields from the PIV experiment and numerical simulations. Figure 6 presents a square region of the vertical symmetry plane going in through the middle of the roof opening. The square is positioned just beneath the ceiling (marked in Figure 7) and corresponds to the region where the PIV measurements were carried out [25]. The presented snapshot was also taken at the 300th second of the fire development(exactly as in the real experiment).

Despite visible differences, which are a result of the different imaging techniques, both velocity fields seem to be very similar, which proves the simulation accuracy. The analyzed velocity distributions correspond to the fully developed air movement in the container. The stream of fresh air coming through the opening is easily visible. Figure 7 shows the velocity vectors recorded at the whole symmetry plane. The data come from the simulation only, as the PIV measurements do not cover the entire area.

As can be seen, a swirl that covered all the container volume was formed. It forced the fire plume to wave around the vertical axis. This effect was much weaker for greater fire powers, which resulted from greater buoyancy forces in comparison with the inertia forces of the moving fluid.

The next step of the numerical research involved study on the backdraft phenomena. For this purpose the simulation was performed in the same configuration as previously until the mass fraction of flammable gases in the container reached the level required for backdraft occurrence. Since a single simulation of the burning phase can last several dozens of hours, the high-power fire was adopted. The assumed HRR was equal to 800 kW. Although the violent course of the combustion phenomenon forced a very short time step, it allowed the under-ventilation conditions to be obtained in a relatively short time (about 57 s of the real process). Additionally, to speed up the calculation, the Ansys Fluent feature ‘symmetry plane’ was enabled, which allowed for the computational domain to be diminished twice. The gases mixture contained large amounts of carbon monoxide and unburned propane. At the chosen time, the door was opened. Figure 8 showsthe calculated temperature distribution in the container at the moment before the door was opened. The fact that the temperature of the gases in the container mostly did not exceed the value of about 1700 K is significant; locally, slightly higher temperatures were observed.

Figure 9 shows the calculated distributions of the mass fraction of carbon dioxide and unburned propane at the same moment. As can be observed, a value of 15% for CO is above the threshold and the value of 20% for C₃H₈ is locally high as well. However, the share of unburned propane in the mixture of flammable gases is rather low.

To monitor the backdraft development, three sets of thermocouples were located at the symmetry plane: Two inside the container and one in the adjacent room. This is shown in Figure 10.

Figure 11, Figure 12 and Figure 13 present the time profiles of the recorded temperature after the door was opened. The profiles are depicted according to Figure 10, and the time of the door being opened (about 57.2 s of simulation) is marked by a dashed line.

The profiles inside the container at the distance of 1m from the door (Figure 12) show significant instabilities caused by the rapid inflow of cold air through the roof opening. This was the result of a sudden pressure drop after opening the door. All profiles are evidence of the deflagration that occurred when the gases mixture with a significant amount of flammable carbon monoxide met the fresh air rich in oxygen. The temperature rose up to 2300 K. This concerns both the inside of the container and the adjacent room. Figure 11 shows a rapid temperature jump in the adjacent room, where the temperature under the ceiling remained high, but after about 1 s, the temperature at heights under 1 m dropped to more or less the initial value. The inside profiles in Figure 12 and Figure 13 are generally shifted in time—the difference is about 0.3 s. This gives an image of the flame front moving towards the inside of the container. The above observations can be confirmed after analyzing the snap shots of the temperature distribution and mass fraction of carbon dioxide in successive time moments. This is shown in Figure 14.

Since according to Figure 9 propane was mainly accumulated just above the burner and the overall amount of unburned propane is significantly smaller than the amount of carbon monoxide, the analyses were focused on the latter one. As can be seen, the regions of the highest temperatures overlap with the regions where hot gases containing large amounts of carbon dioxide meet fresh air. This triggers the ignition of carbon dioxide (and other flammable components, which are not shown), which is the essence of the backdraft. As a result, carbon dioxide burns out at this interface (while the region of high propane abundance is located far from this interface). When the door is opened, the process starts in its whole area, which results in a sudden temperature jump even at low heights. Then, the cold fresh air begins to inflow into the container through the lower part of the door and the temperature drops at low heights. This is in accordance with Figure 11. A stream of fresh air incoming through the roof opening is visible as well. It strengthens the backdraft inside the container. It takes only 0.3 s to spread the blast over the entire volume of the adjacent room. The air movement at 1 s after opening the door is shown in Figure 15.

The velocity distribution coincides with previous observations: Hot burning gases flow out through the upper part of the door and the fresh air inflows through its lower part, forming the gravity current. In this way, the backdraft occurs in both parts—at the interface of both masses of gases. Additionally, the outflow through the roof opening is visible again. This is because the pressure is almost equalized in both rooms and a blast due to backdraft occurs.

The observed variability of pressure also met the expectations. Figure 16 presents the distributions of the total pressure at the symmetry plane for selected moments. The total pressure is the sum of the static and dynamic pressure expressed in reference to the operating pressure. As can be seen, the pressure changes were rather small and ranged from −25 to 15 Pa. The hydrostatic stratification of pressure is clearly visible for all distributions. The first pressure distribution (30th second of the simulation) corresponds to the phase of oxygen abundance. The whole amount of propane was burned, the temperature was high, and the pressure inside was significantly higher than outside. The second one was taken just before the door was opened. The pressure drop can be observed here; it was caused by the fire dimming due to a lack of oxygen. As a result, the temperature inside and the pressure decreased. The third pressure distribution corresponds to the moment 0.3 s after the door was opened. The pressure inside dropped significantly due to the large opening. The last pressure distribution expresses the process of the pressure equalization in both parts.

In the final stage of the work, the partially premixed combustion model was used for the backdraft phase. This was in accordance with the work of Ferraris, who also examined both combustion regimes [21,22].Two sub-models were examined: Chemical equilibrium (the default one) and steady diffusion flamelet. For the latter one, the detailed chemistry was imported from a CHEMKIN file, which was generated at the GRI-Mech 3.0 website [40]. Snapshots of the temperature distribution at successive time moments are shown in Figure 17.

As can be seen, the course of the phenomenon is similar to the non-premixed model. The expansion of hot burning gases is a bit slower, but at 0.3 s, the blast reaches the leftmost edge of the room as well. However, the highest temperature (about 2300 K) is observed just at the interface between the hot gases and the fresh air, which is rich in oxygen. In general, the steady diffusion flamelet sub-model predicts slightly higher temperatures. The distribution of areas with the highest temperature is more similar to the one for the non-premixed combustion model. The differences result from different assumptions of the chemical equilibrium for these combustion models. Since the differences are so subtle, it is difficult to settle which model is the most suitable. This requires accurate experimental data.

6. Conclusions

This work presented a numerical simulation of fire development in a confined space and the backdraft phenomenon. Detailed studies on numerical model features were carried out to achieve a high level of accuracy. Since the work was split in two stages, both of them are summarized separately.

6.1. Validation of the Adopted Assumptions

The validation of the built-up numerical model was possible thanks to the experimental data published by Betting et al. [25]. Figure 4 and Figure 6 prove the high reliability of the numerical model—the results are in good accordance with the reported data. Obtaining such convincing results was possible thanks to a deep study on the numerical modeling of fires and then adopting the relevant assumptions.

The non-premixed combustion model was applied for the fire development stage of the simulation and it was a natural choice for a case in which a stream of fuel entered into a volume filled with an oxidizer.

One of the important observations is that for CFD modeling of the fire development in confined spaces, the RSM turbulence model is the most suitable from the whole RANS family. This is due its capability to model complex anisotropic flows with violent changes of the fluid velocity.

Since the mixture of gases contained optical dense gases like carbon dioxide and water vapor, the DO radiation model was used, but no significant differences were observed compared to the S2S model. This was probably due to the very strong convection, whose contribution to the heat exchange was prevailing.

6.2. Modeling the Backdraft Phenomena

Once the model was validated, the research was extended to the backdraft phenomena. The advantage of the presented approach is the full simulation of fire development—from the ignition moment to the fully developed backdraft. It is different from some previous studies, where the preliminary conditions for backdraft were prepared in advance [20,21]. It is difficult to indicate the phase of spherical propagation probably because it lasts at most a few milliseconds. However, as can be seen in Figure 14, the obtained results showed all other phases:

• The developing fire reached the deeply under-ventilated state;
• Then, after the door was opened, the plane flame front was formed at the interface between hot flammable gases and incoming fresh air (0.004 s);
• The plane flame front stretched through the opening (0.014–0.030 s); and
• Finally, a blast outside could be observed (0.300 s and later).

The intensive gravity current was also reproduced (Figure 15). The pressure variability during both stages of the process was as expected (Figure 16).

In the backdraft stage, both non-premixed and partially premixed models were examined. This approach is in accordance with the work of Ferraris, who examined both combustion models as well [21,22]. Although the course of the simulated phenomena was similar for both models, some differences were observed. The default partially premixed model resulted in a generally lower temperature of the hot gases, which form the blast. In turn, the temperature distribution predicted by the steady diffusion flamelet sub-model was similar to that obtained by the non-premixed model. Hence, the mentioned similarity may suggest that although the backdraft is a quick combustion, the burning conditions were close to the equilibrium state. Despite the observed differences, the overall image of the backdraft predicted by all applied models was coherent.

Ferraris used the LES turbulence model, whereas we used RSM. Although RSM has higher demands on computational resources, it is more accurate for complex flows, especially flows near walls, anisotropic flows, and flows with great gradients. Such kinds of flows can be expected in the case of backdraft, which is a very violent phenomenon. In general, Ansys Fluent allows for more precise modeling of real systems than FDS, which uses LES as the default.

6.3. Further Work

The presented work is just the first stage of the authors’ research. Despite the fact that only one case was considered, the results are promising. Due to the very high fire power adopted, the temperature inside the modeled container was very high (about 1700 K), which made the occurrence of backdraft certain. Once the numerical model is built, it allows for deeper analyses of this dangerous phenomenon. Further research will cover cases of a lower initial temperature and different compounds of the gas mixture. As mentioned above, the use of Ansys Fluent gives hope to examine the course of backdraft in systems with a complex layout as well.

Since the authors are convinced that numerical models should be validated whenever it is possible, many efforts will be dedicated to the acquisition of experimental data by their own experiments or by literature reviews. Therefore, the next stages of our work will focus on examining the backdraft phenomenon under different conditions, while taking into account available experimental data. Accurate data on the temperature values in the blast region would be particularly valuable.