A comparative numerical study of turbulence models for the simulation of fire incidents: Application in ventilated tunnel fires

Abstract

The objective of this paper is to compare the overall performance of two turbulence models used for the simulation of fire scenarios in ventilated tunnels. Two Reynolds Averaged Navier–Stokes turbulence models were used; the low-Re k–ω SST and the standard k–ε model with wall functions treatment. Comparison was conducted on two different fire scenarios. The varied parameters were the heat release rate and the ventilation rate. Results predicted by the two turbulence models were also compared to the results produced from the commercial package Ansys Fluent. Quite faster simulations were performed using the k–ε turbulence model with wall functions and our findings, as to the basic characteristics of smoke movement, were in good agreement with Ansys Fluent ones.

Keywords:

• tunnel fire
• ventilation
• turbulence models comparison
• numerical study

Public Interest Statement

Much attention has been paid on tunnel fires because of their negative consequences in human lives, repair cost, and cost due to the stop of tunnels operation. The methods used for the simulation of a tunnel fire are experimental, theoretical, and computational. Computational simulations cost less than experimental ones and they are more accurate compared to theoretical approaches.

Modeling of turbulence is a serious aspect of a CFD solver. The objective of this paper is the comparison of the overall performance of two turbulence models in two widely studied fire scenarios in a tunnel. Results predicted are also compared to the results of the general purpose commercial package Ansys Fluent. Both turbulence models produced accurate results and in good agreement to those of Ansys Fluent. Greater temperatures downstream from the heat source were calculated compared to experimental results, showing the importance of radiation and wall heat conduction modeling.

1. Introduction

Much attention has been paid on tunnel fires because of their negative consequences. A case of a fire accident in a tunnel constitutes an extremely dangerous situation for people who are inside it, as hot toxic gases are produced. In addition to this, the repair cost and the cost due to the closing of tunnels operation are often huge.

The methods used for the simulation of a tunnel fire are experimental, theoretical, and computational. Full-scale experimental investigation is often prohibitive due to the necessary high cost, but it provides large amounts of reliable data. Many full-scale experiments for the Memorial tunnel are reported in the Memorial tunnel fire ventilation test program-test report (Massachusetts Highway Department, 1995). Small-scale experiments cost less and require careful choice of the scaling factors between the prototype and small scale (Lee & Ryou, 2005, 2006). Analytical methods and Computational Fluid Dynamics (CFD) simulations are the most affordable and permit the investigation of various alternative cases with proper modifications on the model. CFD is the most accurate between them, but a realistic CFD simulation of a fire scenario requires fluid dynamics, turbulence, radiation, wall conduction, and combustion numerical modeling.

Plenty CFD simulations of tunnel fires could be found in the literature, using commercial, open source and fewer of them research codes. When CFD simulation is applied, the selection of the appropriate turbulence model constitutes a serious aspect. Abanto, Reggio, Barrero, and Petro (2006) used Fluent and a research code to study smoke movement in a case of fire in an underwater tunnel. Turbulence was modeled with the standard k–ε model. Hui, Li, and Lixin (2009) have modeled the longitudinal ventilation of the 4th Beijing subway line using CFX and have compared their results for the critical velocity to previous formulations proposed. They preferred the standard k–ε model rather than LES because of computer load reasons. Lee and Ryou (2006) have studied the aspect ratio effect on smoke movement using Fire Dynamics Simulator (FDS) of the National Institute of Standards and Technology and compared their results to reduced-scale experimental ones. Turbulence was modeled using LES model. Hu, Peng, and Huo (2008) have studied the effect of the place of the fire in the critical velocity using FDS and LES turbulence model. Hu, Huo, Wang, and Yang (2007) have compared FDS with LES turbulence model results to full-scale experimental data with promising conclusions for the validity of FDS solver. Wu and Bakar (2000) applied an experimental and CFD investigation on the critical velocity formulations. For the CFD simulation, Fluent was used with the standard k–ε turbulence model.

The most dangerous factor for human lives in tunnel fires is the excess of the concentration limits of the combustion products and not the extreme temperatures. Therefore, the appropriate ventilation system is required for the smoke control. One of the most important aims of a numerical model is the accurate prediction of the back-layering length and critical velocity produced by the current ventilation system. The distance of the smoke front from the heat source is the back-layering length. The critical velocity refers to longitudinal ventilation systems and is the lowest ventilation velocity that could prevent smoke back-layering. Longitudinal ventilation systems drive smoke to tunnel exit, ensuring the safe escape of passengers through the tunnel entrance and/or emergency exits. The ventilation system should prevent back-layering, but high ventilation velocities feed the fire with more oxygen, augmenting the heat release rate (Chow, 1998) and increase the resistance to the passengers, reducing the escaping rate (Hui et al., 2009). For the calculation of the critical velocity, mostly, Froude number-based analytical formulae have been used. Despite their simplicity, they do not account for some specific characteristics of each tunnel such as the existence of lateral evacuation hallways (Banjac & Nikolic, 2008), or the conditions of each fire scenario such as the place where accident happened or the obstructions that may exist. It was found that obstructions affect significantly the critical velocity (Kang, 2006; Oka & Atkinson, 1995). In case of not taking into account all these facts, huge investment costs for the ventilation system may be produced, or inadequate safety measures may be adopted. This fact makes CFD simulation the ideal method for the design of ventilation systems.

Our purpose is the comparison of the overall performance of the low-Re k–ω SST and the standard k–ε with wall functions treatment turbulence models. The developed numerical solver is based on the incompressible Navier–Stokes and energy equations. Buoyancy force, due to temperature differences, was approximated by means of Boussinesq approximation. For both turbulence models, the necessary buoyancy source terms are included. Our results are also compared to those produced by the commercial package Ansys Fluent, using the standard k–ε turbulence model with wall functions under the same simplifications.

2. Mathematical model and numerical procedure

2.1. Test case and governing equations

The test case presented below is a widely studied one. Apte, Green, and Kent (1991) have carried out the experimental investigation. Fletcher, Kent, Apte, and Green (1994) have also presented experimental results and a numerical investigation using a steady-state approach and k–ε model with wall functions. Gao, Liu, Chow, and Fong (2004) have presented a numerical investigation using an unsteady approach and an LES turbulence model. Miloua, Azzi, and Wang (2011) also numerically studied this test case using FDS for the comparison of combustion models and wall boundary conditions, using LES turbulence model.

The tunnel geometry is described in Figure 1. A pool fire exists at 59.5 m from the entrance of the tunnel. The pool fire was assumed to be a cubic volumetric heat source, with the heat release rate being constant and having its maximum value from the beginning till the end of the simulation. The flow field in all cases tested was regarded as incompressible, because Mach number remained below 0.3. For the buoyancy forces, as a result of the temperature differences, Boussinesq approximation was adopted. Moreover, the viscous dissipation term in the energy equation was neglected, because the thermal energy due to viscous shear at low velocities is small (Gao et al., 2004). For the closure of the Reynolds Averaged Navier–Stokes (RANS) mean flow equations, two alternative turbulence models were used and compared; the low-Re k–ω SST and the less computationally demanding standard k–ε with wall functions treatment.

Figure 1. Sketch of the tunnel and heat source.

The general form of the conservation equations is:(1) $$\frac{\mathit{\partial }\mathit{\varphi }}{\mathit{\partial }t}+\mathrm{\nabla }·\left(\overrightarrow{u}\mathit{\varphi }\right)=\mathrm{\nabla }·\left({\mathrm{\Gamma }}_{\mathit{\varphi }}\mathrm{\nabla }\mathit{\varphi }\right)+S_{\mathit{\varphi }}$$(1)

The conserved quantity “ϕ”, the diffusion coefficient “Γ_ϕ” and the source term “S_ϕ” for each one of the conservation equations are given in Table 1 in tensorial notations.

Table 1. Coefficients of Equation 1 and constants values

Conservation equation	ϕ	Γϕ	Sϕ
Continuity	1	0	0
x-momentum	u	${\mathit{\nu }}_{\text{eff}}$	$-\frac{1}{\mathit{\rho }}\frac{\mathit{\partial }p}{\mathit{\partial }x}+\frac{\mathit{\partial }}{\mathit{\partial }{x}_j}\left({\mathit{\nu }}_{\text{eff}}\frac{\mathit{\partial }{u}_j}{\mathit{\partial }x}\right)$
y-momentum	v	${\mathit{\nu }}_{\text{eff}}$	$-\frac{1}{\mathit{\rho }}\frac{\mathit{\partial }p}{\mathit{\partial }y}+\frac{\mathit{\partial }}{\mathit{\partial }{x}_j}\left({\mathit{\nu }}_{\text{eff}}\frac{\mathit{\partial }{u}_j}{\mathit{\partial }y}\right)$
z-momentum	w	${\mathit{\nu }}_{\text{eff}}$	$-\frac{1}{\mathit{\rho }}\frac{\mathit{\partial }p}{\mathit{\partial }z}+\frac{\mathit{\partial }}{\mathit{\partial }{x}_j}\left({\mathit{\nu }}_{\text{eff}}\frac{\mathit{\partial }{u}_j}{\mathit{\partial }z}\right)+g{\mathit{\beta }}_{\mathrm{T}}\mathrm{\Delta }T$
Energy	T	$\frac{\mathit{\nu }}{Pr}+\frac{{\mathit{\nu }}_t}{P{r}_t}$	$\frac{{\dot{q}}_C}{\mathit{\rho }{C}_p}$
k-equation^a	k	$\mathit{\nu }+\frac{{\mathit{\nu }}_t}{{\mathit{\sigma }}_k}$	$p_k+G_k-\mathit{\epsilon }$
ε-equation	ε	$\mathit{\nu }+\frac{{\mathit{\nu }}_t}{{\mathit{\sigma }}_{\mathit{\epsilon }}}$	$[C_{\mathit{\epsilon }1}(P_k+C_{\mathit{\epsilon }3}{G}_k)-C_{\mathit{\epsilon }2}\mathit{\epsilon }]\frac{\mathit{\epsilon }}{k}$
k-equation^b	k	$\mathit{\nu }+\frac{{\mathit{\nu }}_t}{{\mathit{\sigma }}_k}$	$P_k+G_k-{\mathit{\beta }}^{*}\mathit{\omega }k$
ω-equation	ω	$\mathit{\nu }+\frac{{\mathit{\nu }}_t}{{\mathit{\sigma }}_{\mathit{\omega }}}$	$\mathit{\gamma }\frac{\mathit{\omega }}{k}{p}_k+C{D}_{\mathit{\omega }}+\mathit{\gamma }{C}_{\mathit{\epsilon }3}\frac{\mathit{\omega }}{k}{G}_k-\mathit{\beta }{\mathit{\omega }}^2$
${\mathit{\nu }}_{\text{eff}}=\mathit{\nu }+{\mathit{\nu }}_t,P_k={\mathit{\nu }}_t\left(\frac{\mathit{\partial }{u}_i}{\mathit{\partial }{x}_j}+\frac{\mathit{\partial }{u}_j}{\mathit{\partial }{x}_i}\right)\frac{\mathit{\partial }{u}_i}{\mathit{\partial }{x}_j},G_k=-\frac{{\mathit{\nu }}_t}{P{r}_t}g{\mathit{\beta }}_{\mathrm{T}}\frac{\mathit{\partial }\mathrm{T}}{\mathit{\partial }z},C{D}_{\mathit{\omega }}=2(1-f_1){\mathit{\sigma }}_{\mathit{\omega }2}\frac{1}{\mathit{\omega }}\frac{\mathit{\partial }k}{\mathit{\partial }{x}_j}\frac{\mathit{\partial }\mathit{\omega }}{\mathit{\partial }{x}_j}$
$\mathit{\rho }=1.18\frac{\text{kg}}{{\text{m}}^3},g=9.81\frac{\text{m}}{{\text{s}}^2},{\mathit{\beta }}_{\text{t}}=3.34·{10}^{-3}\frac{1}{\text{k}},cp=1,006\frac{\text{J}}{\text{kg}\text{K}},Pr=0.71,{Pr}_t=0.9,C_{\mathit{\epsilon }1}=1.44,C_{\mathit{\epsilon }2}=1.92,C_{\mathit{\epsilon }3}=1.0$

^a Standard k–ε turbulence model is used.

^b k–ω SST turbulence model is used.

The flow variables are the pressure “p,” the Cartesian velocity components “u,” “v,” and “w” in the “x,” “y,” and “z” directions, respectively, and the temperature “T.” Turbulence modeling inserts two additional dependent variables. The kinetic energy of turbulence “k” and turbulence dissipation “ε” when standard k–ε model is used and the turbulence kinetic energy “k” and specific dissipation rate “ω” when k–ω SST is used. The term “gβ_ΤΔΤ” in the “z” direction momentum equation is the buoyancy force term, where “g” is the gravitational acceleration, “β_Τ” is the thermal expansion coefficient, and “ΔΤ” is the subtraction of the ambient temperature Τ_ο (which is equal to 26°C) from the local temperature. “ρ” is the fluid density at the ambient temperature and “C_p” is the specific heat at constant pressure. “${\dot{q}}_c$” is the heat release rate per unit volume. Effective kinematic viscosity “ν_eff” is the sum of the molecular kinematic viscosity “ν” and the turbulence kinematic viscosity “ν_t”. Turbulence kinematic viscosity is calculated by ν_t = c_μk²/ɛ when standard k–ε is used and ${\mathit{\nu }}_t=\frac{{\mathit{\alpha }}_{1}k}{max({\mathit{\alpha }}_1\mathit{\omega },\mathrm{\Omega }{F}_2)}$ when k–ω SST is used. For more information on the k–ω SST turbulence model readers are referred to (Menter, 1993). “Pr” is the Prandtl number and “Pr_t” is the turbulence Prandtl number for temperature. “σ_k”, “σ_ε”, and “σ_ω” are turbulence model constants (Barakos, Mitsoulis, & Assimacopoulos, 1994; Menter, 1993). Terms “P_k” and “G_k” correspond to shear and buoyancy production rates of the turbulence kinetic energy respectively. The inclusion of the “G_k” term in the turbulence conservative equations is of vital importance, since without this extremely false back-layering lengths may be calculated (Fletcher et al., 1994). “CD_ω” is the cross-diffusion term (Menter, 1993).

2.2. Boundary and initial conditions

Two different test cases were simulated, varying in the total heat release rate of the heat source and the ventilation rate. The total heat release rate and the ventilation velocity for these two cases are given in Table 2.

Table 2. Cases data and calculated results

Case	Ventilation velocity (m/s)	Total heat release rate (MW)	Numerical approach	Flame angle (°)	Back-layering length
1	U1 = 0.85	${\dot{Q}}_1$ = 2.57	Present solver	8	Till entrance
			Ansys Fluent	9	Till entrance
2	U2 = 2	${\dot{Q}}_2$= 2.29	Present solver	58	3.2 m
			Ansys Fluent	58	3.7 m

At the inlet of the tunnel uniform velocity profile was prescribed, with the “u” velocity being equal to the ventilation velocity while “v” and “w” velocities were equal to zero (v = w = 0). At the outlet of the tunnel, pressure was prescribed and set equal to the ambient pressure. At the tunnel walls no-slip boundary conditions were applied and they were assumed to be adiabatic, which corresponds to the worst situation. When low-Re k–ω SST turbulence model was applied, turbulence variables at the walls were(2) $$k=0,\mathit{\omega }=10\frac{6\mathit{\nu }}{{\mathit{\beta }}_1{(\mathrm{\Delta }y)}^2}$$(2)

where “Δy” is the distance to the next point away from the wall. When high-Re k–ε turbulence model was applied, appropriate wall functions were used leading to coarser numerical meshes and less time demanding calculations. The values for “u”, “v”, “w”, “T”, “k”, and “ε” were calculated using the following wall functions.(3) $$u^{+}=\left\{\begin{array}{cc}\frac{1}{\mathit{\kappa }}ln{y}^{+}+5.5,\hfill & y^{+}\ge 11.6\\ y^{+},\hfill & y^{+}\le 11.6\end{array}\right.$$(3) (4a) $$T^{+}=Pr{y}^{+}{e}^{-\mathrm{\Gamma }}+\left[2.12ln\left(1+y^{+}\right)+\mathit{\beta }(Pr)\right]e^{-1/\mathrm{\Gamma }}$$(4a) (4b) $$\mathit{\beta }(Pr)={\left(3.85{Pr}^{1/3}-1.3\right)}^2+2.12ln(Pr)\mathrm{\Gamma }=\frac{0.01{\left(Pr{y}^{+}\right)}^4}{1+5{Pr}^3{y}^{+}}$$(4b) (5) $$y^{+}=\frac{y{u}_{\mathit{\tau }}}{\mathit{\nu }}$$(5) (6) $$u^{+}=\frac{u}{u_{\mathit{\tau }}}$$(6) (7) $$T^{+}=\frac{T_w-T}{T_{\mathit{\tau }}}$$(7) (8) $$u_{\mathit{\tau }}=\sqrt{\mathit{\rho }{\mathit{\tau }}_w}$$(8) (9) $$T_{\mathit{\tau }}=\frac{q_w}{\mathit{\rho }{C}_p{u}_{\mathit{\tau }}}$$(9) (10) $$k=\frac{u_{\mathit{\tau }}^2}{\sqrt{C_{\mathit{\mu }}}}$$(10) (11) $$\mathit{\epsilon }=\frac{u_{\mathit{\tau }}^3}{{\mathit{\kappa }}_y}$$(11)

For y⁺ ≤ 11.6 the first grid point is in the viscous sub-layer and for y⁺ > 11.6 in the logarithmic region. It is desirable that the first grid point belongs to the logarithmic region of the boundary layer, where the viscous effects are weaker compared to the turbulent ones (Launder & Spalding, 1974). For the temperature near the wall, the formulation proposed by Kader was used (Kader, 1981), which is valid in all zones of the boundary layer and for a large range of Prandtl numbers.

As initial condition, the converged to steady-state flow field of the isothermal case (without the heat source) was used. The temperature was equal to the ambient temperature.

2.3. Numerical methodology

The solution procedure of the mean flow equations has been analyzed in details by Stokos, Vrahliotis, Pappou, and Tsangaris (2012, 2013) and Vrahliotis, Pappou, and Tsangaris (2012). A node centered finite volume discretization technique has been developed for hybrid numerical meshes. The solution of the incompressible flow equations was based on the artificial compressibility method. A Roe’s approximate Riemann solver was used for the evaluation of the convective fluxes. For the discretization of the viscous fluxes a central scheme was used. For the evaluation of the convective fluxes of the turbulence equations a first-order upwind scheme was applied, as proposed by Menter (1993). For the calculation of turbulence equations viscous fluxes, the same method used for the mean flow equations was applied. For the temporal discretization, a dual-time stepping approach was adopted. Specifically, an implicit second-order backward difference scheme was implemented for the physical-time marching and an implicit first-order scheme for the pseudo-time marching. To pass from one physical time step to another convergence in pseudo-time needed to be achieved. The mean flow and turbulence equations were loosely coupled. The mean flow equations were first solved with the turbulence kinematic viscosity fixed. Then, the turbulence equations were solved with the velocity and temperature fields fixed. The resulting non-linear algebraic system of equations was linearized using Newton’s method. The arising linear system was solved using the Jacobi iterative method with suitable under-relaxation. Parallel processing was used for the acceleration of the computations.

2.4. Numerical mesh

The computing domain was composed of 481,950 nodes and 449,000 hexahedrons when standard k–ε turbulence model with wall functions was used (Figure 2). The first layer thickness was approximately equal to 0.01. When low-Re k–ω SST was used the computing mesh was composed of 1,094,252 nodes and 1,075,400 hexahedrons (Figure 3). The first layer thickness was of the order of 10⁻⁵ so that the y⁺ < 2 condition is fulfilled and the growing factor was equal to 1.2. Both numerical meshes were denser near the heat source. The physical time step was set equal to 0.01 s.

Figure 2. Numerical mesh used for the standard k–ε with wall functions simulations.

Notes: Section y = 2.7 m near the vicinity of the heat source (top figure). Cross-section x = 60 m (bottom figure).

Figure 3. Numerical mesh used for the low-Re k–ω SST simulations.

Notes: Section y = 2.7 m near the vicinity of the heat source (top figure). Cross-section x = 60 m (bottom figure).

3. Results and discussion

3.1. Transient results

When no concentrations of smoke are calculated and radiative heat transfer is not considered, it could be assumed that smoke movement is analogous to the temperature field. In Figures 4 and 5, the evolution of the flame and temperature field for cases 1 and 2, respectively, for the first 10 s are presented. It is obvious that the predicted smoke movement through time by the two turbulence models is similar.

Figure 4. Temperature isolines through time for case 1.

Notes: (a) Standard k–ε at 0.5 s; (b) Low-Re k–ω SST at 0.5 s; (c) Standard k–ε at 2 s; (d) Low-Re k–ω SST at 2 s; (e) Standard k–ε at 5 s; (f) Low-Re k–ω SST at 5 s; (g) Standard k–ε at 10 s; and (h) Low-Re k–ω SST at 10 s.

Figure 5. Temperature isolines through time for case 2.

Notes: (a) Standard k–ε at 0.5 s; (b) Low-Re k–ω SST at 0.5 s; (c) Standard k–ε at 2 s; (d) Low-Re k–ω SST at 2 s; (e) Standard k–ε at 5 s; (f) Low-Re k–ω SST at 5 s; (g) Standard k–ε at 10 s; and (h) Low-Re k–ω SST at 10 s.

In a tunnel fire case surrounding medium (air and smoke) in the vicinity of the heat source is heated and raises up till the ceiling of the tunnel. Then having reached the ceiling it moves to the side walls and along the ceiling to the tunnel exit and the tunnel inlet forming the back-layering length. Reaching the side walls, smoke moves downward to the ground. In Figures 6 and 7 velocity vectors predicted by both turbulence models are given at characteristic sections and moments. Velocity vectors reveal smoke movement. It is obvious that smoke requires less than 1 s to reach the ceiling and less than 5 s to reach the side walls for both cases. Similar flow patterns were predicted by the turbulence models.

Figure 6. Velocity vectors at characteristic sections and moments for case 1 predicted by the standard k–ε model (left column) and the low-Re k–ω SST model (right column).

Notes: Section y = 2.7 m near the heat source at 0.5 s (top row), cross-section x = 60 m at 1 s (middle row), and cross-section x = 60 m at 5 s (bottom row).

Figure 7. Velocity vectors at characteristic sections and moments for case 2 predicted by the standard k–ε model (left column) and the low-Re k–ω SST model (right column).

Notes: Section y = 2.7 m near the heat source at 0.5 s (top row), cross-section x = 62 m at 1 s (middle row) and cross-section x = 62 m at 5 s (bottom row).

In Figure 8, we present the predicted using both turbulence models temperature vertical profiles at 18 and 30 m downstream from the heat source for case 2, 20 s after fire breaking. Small differences (less than 6%) are observed for the vertical temperature profiles.

Figure 8. Case 2, 20 s after fire breaking.

Notes: Vertical temperature profiles at a distance of 18 m (top figure) and 30 m (bottom figure) downstream from the heat source.

In Figure 9, we present for case 2 the convergence history for the temperature corrections in pseudo-time for three physical time steps, using the standard k–ε turbulence model. The starting point is t = 2 s. It is obvious that less than 250 pseudo-time steps are required for each physical time step. However, it should be noted that the first 2 s are the most numerically crucial and more pseudo-time steps are required for pseudo-time convergence compared to the rest of the simulation. Same smooth convergence occurred and for the other dependent variables.

Figure 9. Convergence curve.

Notes: Temperature corrections in pseudo-time for three physical time steps with starting point t = 2 s.

3.2. Steady state results

Gao et al. (2004) claim that flame shape is defined by the maximum temperature gradients. In the literature, many definitions for the flame angle have been found (Anderson et al., 2006). In Figures 10 and 11, flame shapes and the definition used for the calculation of the flame angle are given. Flame angle β_t is defined by the vertical line passing through the core of the heat source and the line connecting the core of the heat source with the upper point of the flame. In the same figures, temperature fields in the vicinity of the heat source are compared to those predicted by Ansys Fluent. The greatest temperature values, flame shapes, flame angles and back-layering lengths seem to agree satisfactorily. Flame angles and back-layering lengths predicted by our solver and Ansys Fluent are given in Table 2. It is evident that the backflow is less and flame tilt greater for higher ventilation velocities. The time needed for the smoke front to reach the steady-state back-layering length was about 50 and 15 s for case 1 and case 2, respectively.

Figure 10. Flame shape comparison for case 1.

Notes: Present solver with the k–ω model (top figure), Ansys Fluent with the k–ε model (bottom figure).

Figure 11. Flame shape comparison for case 2.

Notes: Present solver with the k–ε model (top figure), Ansys Fluent with k–ε model (bottom figure). Definition of flame angle.

Radiation plays a significant role in a case of fire in a tunnel. The fraction of the heat release from a heat source in a tunnel fire in form of radiation is in the order of 20–50% (Bettelini, 2001; Hostikka, 2008; Grant, Jagger, & Lea, 1998). However, the solution of the radiative heat transfer equation is time consuming because a great number of radiation intensities have to be calculated at each computational node. More than 100 solid angles are utilized for each node at a finite volume-based method. Consequently, heat transfer due to radiation is often not taken into account, or the assumed amount of radiation loss is subtracted from the heat release rate (Kang, 2006). In Figure 12, temperature profiles are given 18 and 40 m downstream of the heat source. Greater values for the temperature are computed compared to the experimental values of Fletcher et al. (1994). This discrepancy is attributed to the omission of radiation modeling, according to the aforementioned role of radiation, and the omission of heat conduction inside tunnel wall. Similar discrepancies were calculated by Gao et al. (2004) who applied the same simplifications. However, the curves calculated are of the same “s” form. Temperature increases with increasing height and reaches its maximum values inside the smoke plume. Both turbulence models calculated approximate temperature profiles 18 and 40 m downstream from the heat source. Experimental investigation also showed approximate temperature profiles 18 and 40 m from the heat source.

Figure 12. Comparison of temperature profiles along height when steady state was reached.

Notes: Case 1 (top figure), and case 2 (bottom figure).

Velocity vectors for case 2 with the standard k–ε model when steady state was achieved are presented in Figure 13. The stagnation point which defines the back-layering length is marked with a red circle.

Figure 13. Velocity vectors at section y = 2.7 m and stagnation point (red circle).

Finally, as far as computation time is concerned, a speedup approximately equal to 2.5, at each pseudo-time step, was gained when k–ε turbulence model with wall functions was used compared to the low-Re k–ω SST turbulence model. This fact is mostly a result of the less dense grid required from the k–ε turbulence model due to the use of wall functions.

4. Conclusions

The performance of the standard k–ε with wall functions and the low-Re k–ω SST turbulence models was compared in terms of accuracy, stability, and calculation time. Flow pattern and temperature fields were in good agreement. Our results were also compared to the results produced by the commercial general purpose code Ansys Fluent. Back-layering lengths and flame angles for both cases agreed satisfactorily with Ansys Fluent ones. However, the omission of radiation and wall conduction modeling led to higher temperatures compared to experimental results. Simulations including both of these models will be presented in a future article with improved calculated temperature fields.

Acknowledgments

Simulations were performed on the high-performance platform “VELOS” of the Laboratory of Thermal Turbomachines, Parallel CFD and Optimization Unit, School of Mechanical Engineering of the National Technical University of Athens.

Additional information

Funding

Funding. This work was funded under the EUROSTARS Project (E!5292) “Structural and Aerodynamic Design of TUNnels under Fire Emergency Conditions” from: the Greek General Secretariat for Research and Technology, Bundesministerium fur Bildung and Forschung (BMBF), FiDES DV-Partner GmbH, Germany and ELXIS Engineering Consultants SA. It was also funded from the Greek State Scholarships Foundation.

Notes on contributors

Konstantinos G. Stokos

Our group’s research and working areas relate to the development and application of software for computational mechanics. Specifically, the development and application for heat transfer problems (conduction, convection, and radiation), structural problems, aerodynamics, fluid–structure interaction, biofluid mechanics, numerical mesh generation, and optimization constitutes our field. Our team is involved in several transnational research projects in collaboration with organizations and companies. Many papers of our work have been published in scientific journals and international conferences.