Some Tips on Numerical Modeling of Airflow and Fires in Road Tunnels

Abstract

The efficiency of tunnels systems is often evaluated using numerical simulations. This concerns both to normal and emergency mode of tunnel systems operation. Therefore the safety level of tunnel users may depend on the quality of numerical models being built. The most often studied areas cover the researches on natural and forced airflows in the normal mode and on fire development and smoke spreading in the emergency mode as well as modeling of fan operation. Thus, many software packages implementing Computational Fluid Dynamics (CFD) are applied here. Despite the available software is recognized as reliable, the problem arises because the built numerical models should be validated at least partially with experimental data. There is a shortage of experimental data from real tunnels due to high costs and many organizational or formal difficulties. Some researchers use data from scaled experiments, but this leads to problems connected with scaling. The paper presents the application of two widely used software packages—Fire Dynamics Simulator (FDS) and ANSYS Fluent to reproduce some scenarios of the operation of a tunnel ventilation system for normal and emergency mode. Most of results were compared with data obtained by own full scale measurements or data available in literature. Some practical issues concerning the application of FDS and ANSYS Fluent were discussed as well.

1. Introduction

Accurate airflow reproducing is the crucial issue while modeling the operation of a tunnel ventilation system at any mode or fire development and smoke spreading in a case of fire [1]. Since commonly used software packages which implement Computational Fluid Dynamics (CFD) are recognized as reliable and soundly validated, the main problem is to built a suitable and faultless numerical model. The two most popular software packages are Fire Dynamics Simulator (FDS) and ANSYS Fluent. The first one is willingly used mainly for simulation of fire development [2,3]. The latter one is a multipurpose package, and can be used for any issue related to fluid and heat flow [4]. There are some design decisions to be made before the start of modeling. All approximations and simplifications have to be justified.

Some of the issues presented here have already been somehow considered, but they may still be a challenge when facing the task of accurate modeling of airflows and fires in tunnels [5,6]. But very often numerical models of tunnels are built without paying sufficient attention to them and eventually it may lead to results which are burden with significant discrepancies. The matter goes worse if there is no possibility to compare numerical results to experimental data and the results become regarded as reliable. The article deals with the following issues:

• modeling of the heat properties of tunnel walls,
• taking into account the drag caused by protruding details of tunnel infrastructure (lamps, traffic signs, cables, fans),
• modeling of jet fans,
• selecting of the turbulence model,
• taking into account the influence of external wind,
• modeling the stack effect,
• monitoring of the solution convergence.

Obviously, the presented above list is not exhaustive. There are some issues not addressed, for instance the influence of vehicles movement on airflow. However, the discussed matters concern the most common problems faced by researchers dealing with flows in road tunnels. Whenever it was reasonable references to experimental data, which validated presented models were provided. The rest of the paper is divided into sections dedicated to successive topics. The last goes the conclusion section, which sums up the most important findings.

2. Modeling of Tunnel Walls

The first step while developing a tunnel model is to build the geometry and specify the properties of the materials the tunnel is built of. This is a very important stage of the work. The tunnel model should accurately reflect the real object. The geometry of the tunnel should be accurately mapped and so the materials it is made of. This will allow for proper modeling of heat flow and reproduction of temperature distributions in case of fire analyzes. The material from which the tunnel vault is made has the ability to accumulate heat, which will have an impact on its distribution in the tunnel space. Only when the considered flow is regarded as isothermal the tunnel walls could be modeled as adiabatic surfaces. Otherwise the heat transfer between air or other gases and tunnel walls must be taken into account. Hence, studies can be found in which the walls of the tunnel are modeled as concrete partitions [7]. Among others, Hua modeled the fire damage in a tunnel structure. He provided a framework to characterize the distributions of fire temperature within a tunnel structure. However, the information on the way of modeling tunnel walls is most often omitted. This can be accepted when the manuscript deals only with the issue of cooled smoke flow in the tunnel, for example due to the evacuation process [8]. However, when modeling the complete picture of smoke flow, it should be remembered that the tunnel wall temperature may affect the results of such analyzes. When there is no information, the tunnel walls were modeled usually as “default”, and in most CFD programs such default setting means the adiabatic boundary conditions. This even applies to manuscripts where the temperature distribution caused by fire was considered [9,10]. Similarly, Kim considered the effect of the tunnel slope on plug-holing phenomenon. He investigated the temperature distribution at different tunnel inclinations, assuming at the same time the adiabatic walls of the tunnel [11].

It is significantly important for modeling of fire development and hot smoke spreading. The commonly accepted scenario assumes that due to the heat transfer from hot gases to tunnel ceiling and sucking cold air the smoke layer is cooled down and therefore it lowers and the initial stratification disappears.

To show the importance of the discussed matter a numerical model was prepared using FDS software. A tunnel of cross-section of 8 m × 6.4 m and 300 m long was modeled. A fire of HRR equal to 4 MW was placed at the middle of the tunnel length. The fire grew linearly during 60 s. Additionally a dynamic pressure of 1Pa was applied to one of the tunnel portals to simulate a weak natural airflow. Figure 1 shows temperatures recorded 0.05 m beneath the ceiling at selected distances downstream the fire.

These outcomes confirmed the expectations: for adiabatic ceiling the temperatures were higher and the difference was more clearly visible with increasing distance from the fire source. In the presented above numerical experiment so called layer zoning devices to measure smoke layer height were also applied. Although the temperature differences are significant, the recorded heights of the smoke layer didn’t differ in such a unequivocal way. These measured differences were of order of fluctuations. It can be explained by the fact, that in an early fire stage for both wall models, despite a clear temperature difference the upper layers were hot enough to keep the similar smoke stratification.

3. Modeling the Additional Drag Caused by Details of Tunnel Infrastructure

Tunnels walls and ceiling are richly equipped with many items (lamps, cables, fans, sensors, road signs), which generate drag against airflow. This can significantly change a perpendicular distribution of air velocity and eventually affect the magnitude of average longitudinal velocity. This drag need to be taken into account when modeling airflow in tunnels. One should have in mind that modeling inaccuracies concern especially regions close to walls and ceiling. Therefore they may impact on the obtained picture of such phenomena as backlayering [12,13] and plug-holing [11], which are essential for assessing the efficiency of ventilation systems.

It is difficult to model all of these details accurately due to very high workload. Fortunately, this additional drag can be modeled by introducing so called equivalent roughness. It is the same way as the natural wall roughness is accounted for [14]. Thanks to the carried out measurements in a 678 m long road tunnel [15] it was possible to depict distributions of air velocity at selected cross-sections of the examined tunnel for different configurations of switched on jet fans. A set of 5 tall and 6 short stand-poles with mounted thermo-anemometers was used for precise airflow velocity measurement. Then using ANSYS-Fluent the flows were reproduced for different heights of wall roughness (h_r). The roughness of the road surface was adopted as 0.002 mm, which is the natural roughness of tarmac. The obtained results were compared with the measurement data. An example for two jet fans operating in the normal mode far from the measurement cross-section is shown in Figure 2.

As a measure of the modeling accuracy the root mean square deviation (RMSD) was adopted. The relation between the value of RMSD and the height of wall roughness (h_r) can be well approximated by a parabola (determination coefficient R² = 0.954), which is shown in Figure 3.

As it can be seen the best fit was obtained for the equivalent wall roughness equal to 0.1 m. Although, the velocity values predicted by Fluent for regions close to walls were still higher than real ones, the overall accuracy was the highest. If the equivalent wall roughness increased the value of RMSD would rapidly increase and the average airflow velocity would be much too low. Since the infrastructure of the examined tunnel can be regarded as a typical one, it is a premise that this value of equivalent wall roughness can be applied commonly when modeling airflows instead of taking into account all equipment details.

4. Modeling of Jet Fans

Jet fans are obviously the most important components of a longitudinal ventilation system. Therefore the accurate airflow modeling strongly depends on the way the fans are introduced into a numerical model. The fundamental formula governing a fan operation describes the relation between outlet velocity (u_out) and unit force (S_a, ratio of the fan thrust to its volume), which corresponds to the pressure jump is given below (ρ denotes fluid density and l denotes the fan length) [16]:

$$S_a=\frac{\rho u_{out}^2}{l}$$

(1)

A very popular approach used in ANSYS Fluent is to set a certain value of the fan outlet velocity (denoted here as IO model) [17,18]. But except the simplicity this fan model has no other advantages. The main drawbacks are unrealistic radial distribution of air velocity and faked composition of the exhaust air (in a case of fire modeling). The summary of jet fan models available in ANSYS Fluent is shown in Figure 4 and briefly described in

Table 1. A brief description of fan models, based on [19].

Model	Description	Input Data
IO	A pair of surfaces of type ‘inlet’ and ‘outlet’	Outlet velocity
A	A region with the volume momentum source	Unit force corresponding to the pressure jump
B	Infinite thin rotor surface	Pressure jump or fan characteristics, optionally axial and radial velocity profiles
C	A cylindrical volume corresponding to a rotor of real dimensions	Pressure jump or fan characteristics

[19].

Models B and C are capable to reproduce accurately the distribution of all components of airflow velocity in a fan jet at almost every distance from an outlet. However, they require more computational resources than other models and they are not very popular. Model B was used to analyze the pressure loss due to traffic jam conditions in a road tunnel [20]. Model A can be regarded as a reasonable tradeoff between accuracy and complexity. This model is encountered in numerical analyzes of the reduced-scale tunnel with longitudinal ventilation system where jet fans that have the inlet/outlet sections inclined at a fixed pitch angle (α = 6°) toward the tunnel floor [21]. Excluding the region just behind the fan outlet the jet is reproduced very accurately. What is important, with a given pressure jump such a fan operates properly for variable air density (for instance due to variable temperature or variable gas composition).

The axial velocity distributions along fan axis for a jet for all discussed fan models are shown in Figure 5. Additionally the experimental data obtained by Giesien team in 2011 [22] are shown. As Giesen reported the diameter of the fan was 0.29 m and its length was 2.6 m. To strictly keep the experiment conditions the following parameters were adopted: outlet velocity of 18 m/s for IO model, unit force of 165.1 N/m³ for A model and pressure jump of 448 Pa for models B and C.

Analyzing the presented data, one can see that all curves but the one of IO model almost converge at a distance shorter than 5 m. Models B and C imitate the flow in the region just behind the fan outlet with a higher accuracy. Outputs of A model are the worst in this region, but its curve very quickly tends to the others. Therefore, having in mind other advantages it is widely applied for modeling of flows in entire tunnels space.

In FDS jet fans are modeled in the way similar to A model: two HVAC components acting as an inlet and an outlet are connected with a duct (Figure 6). However, the fan capacity is determined by a set value of volume flow rate (not thrust), which makes it vulnerable to variable air density [23].

There are manuscripts analyzing the phenomena taking place in tunnels, in which fan modeling is completely omitted. A constant flow at the tunnel inlet is assumed as velocity inlet [24,25]. This distorts the results of the analyzes because the stream generated by the axial fan never covers the entire cross-section of the tunnel. This means that the velocity distribution at a tunnel cross-section is not the same, it is related to the shape of the jet and the roughness of the tunnel walls.

5. Selecting of the Turbulence Model

Fluid flow modeling is a very complex problem mainly due to turbulences. Generally, in CFD approach the entire computational domain is divided into small cells. For each such cell, balance equations should be solved for all quantities being considered, at least for mass, energy and momentum. Since apart from simplest cases the direct numerical solution of these equations, which govern the fluid flow is nowadays still generally impossible due to technical reasons, many different simplifications are introduced to cope with this issue in reasonable time. There are two main approaches used here: Large Eddy Simulation (LES) and Reynolds Averaged Navier-Stokes (RANS).

In the first one only the largest eddies are directly solved and the smaller eddies are modeled analytically [26]. It is that the predominant part of eddies is smaller than cell size, so a low-pass Favre filter is applied to the transport equations for mass, momentum and energy. Hence, just mean values of considered quantities are calculated for each cell [27]. Because of its nature LES method requires transient calculation mode.

The latter approach assumes that from engineering point of view just averaged values are important. Thus the instantaneous fluid velocity (u) is regarded as a sum of the average velocity ($\overline{U}$) and fluctuating velocity (u′):

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

(2)

By definition the average of fluctuating component is equal to zero ${\overline{u}}_i^{\prime }=0$. By substituting this to the basic Navier-Stokes equations, a new quantity called Reynold’s tensor of viscous stresses is achieved. The values of the tensor’s components can be solved directly, what leads to the Reynolds Stress Model (RSM). The RSM model is accurate even for complex and anisotropic flows with strong gradients but is computationally demanding. Further simplifications adopt Bussinesq’s hypothesis, which assumes that turbulent stresses are transported in a way similar to viscous stresses. Hence the components of the Reynold’s tensor are expressed by a new single variable—turbulent viscosity. Depending on some details there are two models developed: k-ε and k-ω. Term k corresponds to the kinetic turbulence energy, ε is dissipation rate of turbulence energy and ω is specific dissipation rate of turbulence energy (ω = ε/k). Both of them are nowadays somewhat like industry standards, for instance Chen et al. and Zhao et al. used k-ε model [14,28], meanwhile Król and Król used k-ω model [29,30]. However using them one must be aware of their drawbacks and limitations. Very simplifying the issues, the k-ε model better describes the bulk flows, meanwhile the k-ω model performs better for flows near walls.

FDS uses LES turbulence model by default, thus a list of researchers, who use LES is long [31,32,33]. It is dedicated for low-Mach buoyancy driven flows, which are specific for the fire phenomenon. Meanwhile ANSYS Fluent supports a majority of known models.

The performances of different turbulence models were examined again using experimental data obtained by Giesien team [22]. The fan model of type A was applied and different turbulence models were checked for free and wall jets [5,19]. The parameters of the A fan model were adopted as earlier in Section 4.

Since LES can work only in the transient mode the data shown for it came from averaging of 10 series recorded successively by 2 s. As a quantity being a comparison criterion the distribution of the axial velocity for free and a wall jets was selected. It was considered over a distance of 5 to 15 m behind the fan outlet (Figure 7).

The exact comparison of LES with other turbulence models is difficult due to different nature of them. But it is clearly visible that for free jets all models of the RANS family gave almost the same predictions, which were in accordance with the measurement. Meanwhile LES predicted much higher velocity, especially close to the fan outlet. At larger distances LES outputs fit with the measurements in a higher degree.

The situation was different for wall jets: the accordance between the RSM prediction and the measurement was almost perfect, other RANS models overestimated the air velocity and the LES predictions oscillated between both curves.

The observed divergences decreased as the distance from fan outlet increased. Thus from practical point of view, when fully developed flows in long tunnels are considered the differences are meaningless and this is why all mentioned turbulence models can be and actually are widely applied. But if phenomena close to a fan outlet are examined the selection of a turbulence model should be preceded by a deep analysis.

6. Impact of the External Wind

When considering the phenomena occurring in a road tunnel equipped with a ventilation system, it should be taken into account that flows in the tunnel under real conditions are affected by the external wind. Obviously, the wind has a particular importance when natural ventilation is designed in a tunnel. Here, the wind is the main driving force affecting the airflow next to the stack effect, which will be discussed later and the impact of the road traffic. Meanwhile, analyzes often ignore the effect of wind on the flows inside the tunnel. Simply, a specific flow velocity in the tunnel is assumed, without specifying whether this velocity results only from the operation of fans or also from the influence of the external wind [9]. Sometimes researchers consciously assume that the pressure inlet at the tunnel portal is 0 Pa. This means that the influence of external wind is not taken into account in the analyzes and this is clearly indicated [20].

The most accurate way to take into account the impact of the external wind is to accurately model the whole vicinity of tunnel portals. The added volumes should be of size up to several times greater than respectively tunnel portal width and height. Such approach leads in an obvious way to an abnormal extension of the computational domain. Figure 8 presents the numerical model of the influence of the external wind, which contains the outer volume neighboring to the windward portal.

This can be avoided by an application of a certain value of dynamic pressure (Δp) directly to the tunnel portals according to well known Bernoulli’s formula (u_wind denotes the external wind velocity):

$$\Delta p=c\frac{\rho u_{wind}^2}{2}$$

(3)

The value of c coefficient should be adjusted to obtain the same velocity of natural airflow as for real external wind. It depends on the geometry of the vicinity of tunnel entrance and may require the method of trial and error.

This issue was examined by comparing two numerical models: the first was the full one (as in Figure 8), the second contained just the tunnel interior, but the relevant dynamic pressure was applied at the entrance portal.

The value of the dynamic pressure was calculated using the former numerical model at a vertical plane just behind the tunnel entrance. Figure 9 presents the relation between the additional pressure caused by the wind and the external wind velocity. A perfect accordance with Formula (3) can be observed.

Then the calculated values of dynamic pressure were applied to the second model. Additionally, the airflow velocity inside the tunnel was calculated using well known semi-empirical Darcy-Weisbach formula with wall roughness adopted as above (D_h denotes a tunnel hydraulic diameter, L denotes its length, λ is Darcy drag coefficient depending on Reynolds number and wall roughness):

$$\Delta p=\lambda \frac{\rho u_{air}^{2}L}{2D_h}$$

(4)

Using Moody’s diagram and fitting the relative roughness height h_r/D_h the value of Darcy drag coefficient can be determined as λ = 0.042. Figure 10 shows the comparison of the obtained results. As it can be seen almost exact accordance was obtained. However, one must have in mind that the relevant value of dynamic pressure representing the wind impact must be determined in advance in some way.

7. The Stack Effect Modeling

A stack effect is the second cause of the natural airflow occurrence. This appears when colder air enters the inclined tunnel whose walls are warmer. Such temperature difference could be caused by the circadian rhythm or the occurrence of geothermal sources. The heat is transferred from walls to air and warms it, hence air density decreases and due to buoyancy forces it starts to move upwards. In such a way an inclined tunnel acts as a heat engine. Its cycle runs between the temperature of the tunnel wall (heat source) and the ambient temperature (cooler).

There are two ways to reproduce this phenomenon. The simplest one is to introduce additional dynamic pressure at the lower tunnel portal in a similar way as for external wind [34,35]. Assuming air is an ideal gas, this thermally induced pressure can be expressed by the following formula [36,37]:

$$\Delta p={\rho }_0{T}_{amb}\left(\frac{1}{T_{amb}}-\frac{1}{T_i}\right)g\Delta h$$

(5)

In this formula T_i and T_amb denote the average air temperature inside the tunnel and the ambient temperature respectively, Δh denotes the difference of portal heights.

The more complex approach requires to model tunnel walls as a boundary condition with given temperature or even to model the whole tunnel shell of concrete with given temperature deep in the bulk. To show this some simulations were carried out for two tunnels inclined by 3% and 6%. The ambient temperature was 7 °C, a the airflow was examined for selected values of tunnel walls temperature (the phenomena inside walls bulk were not taken into account). The numerical model contained also the volumes in the vicinity of both tunnel portals. Its structure and the applied boundary conditions are shown in Figure 11. The external wind velocity (at ‘velocity inlet’ boundary condition) was set to 0 for this case.

The average airflow longitudinal velocity at the cross-section in the middle of tunnel length was regarded as the result. It is presented in Figure 12. Additionally values used in Equation (5) are shown. A crucial issue was to specify the temperature of air inside the tunnel, because it is not equal to the tunnel walls temperature. Eventually the average value coming from the mentioned full simulation was taken. Then, as previously the flow velocity was determined using Darcy-Weisbach formula for calculated pressure value.

The first observation is quite obvious: the simulation results proved that the naturally induced airflow could be significant enough to impact on the operation of a ventilation system. In turn this suggests that the contribution of the stack effect must be taken into account when modeling airflows in an inclined tunnel. The approach which uses Equation (5) gives the accurate predictions on airflow velocity but the real air temperature inside the tunnel must be known.

The described above approaches were also validated with experimental data: using data on wind velocity and pressure at the tunnel portals, temperature at the portals and inside the tunnel both additional pressure components were calculated according to Equations (3) and (5), then the flow velocity was determined using Equation (4) and compared with the measured natural airflow velocity. The accuracy about 13% was achieved (1.11 m/s measured and 1.28 m/s calculated) [35].

An additional issue must be considered when modeling in ANSYS Fluent flows which are driven by the buoyancy in inclined tunnels. The gravity must be switched on (while, it is FDS default). It is recommended to set operating density to 0 and to use ‘Body Force Weighted’ method for pressure spatial discretization. If any pressure boundary conditions are applied, the pressure value must be set accounting for the pressure drop along with the height (p₀ is the reference pressure):

$$p=p_0-\rho gh$$

(6)

This will generate the natural atmospheric pressure profile (again it is FDS default) It can be done with the use of a User Defined Function (UDF), which controls the given boundary conditions. For the discussed simplified models the contributions of wind influence and stack effect can be taken into account here together as additional pressure components as well.

8. Monitoring of the Solution Convergence

The issue considered in this section concerns mainly steady mode numerical modeling in ANSYS Fluent(or other similar codes). It may seem to be of less importance than others, but neglecting it might lead to significantly inaccurate results, which is noticed by some researchers [38,39]. Despite its apparent simplicity (just a tube) tunnels pose specific challenges in modeling. It is a result of the aspect ratio of the computational domain. Tunnels are extremely long in comparison to their width or height. One must be aware that the common way of monitoring the convergence of solution could be inadequate for some calculations, especially in steady mode. This can be illustrated in a simple numerical model of jet fan operation, which was built up in ANSYS Fluent. In the presented example a steady flow in a long tube (length 300 m, diameter 6 m) for a given fan capacity (pressure jump 450 Pa) was sought using k-ω turbulence model. The fan of A type was placed at the centerline of the tube.

The default convergence criteria (10⁻⁶ for energy and 10⁻³ for others residuals) were switched off. Figure 13 presents the solution progress for different initial values of air velocity. The points, for which the solution would converge according to the default convergence criteria settings are marked.

As it is easily visible there is a difference of 0.35 m/s between both solutions. Therefore it is advisable to switch off the default convergence criteria and introduce additional monitored values and use them as a convergence index.

9. Conclusions

The article deals with some details concerning modeling of airflow and fire development in road tunnels. The considered issues might seem trivial for those who are experienced in this area, but the literature review proved that some authors seem to be quite unaware of the importance of the matter discussed here. Meanwhile, similar papers on good practices have appeared for other domains, in which CFD methods are also involved [38,39]. The problem arises because the calculated flows almost never fit accurately the measured ones, and this is an obvious fact to be reconciled with. But being aware of the pitfalls one can prepare his numerical model as precisely as possible. On the other hand some necessary tradeoffs, which may result in inaccuracy would be adopted with the full awareness.

The presented work is a result of plenty of time spent when trying to fit numerical models of tunnel airflows and fires to experimental data, which came from own or literature measurements. The carried out numerical experiments showed that:

• in cases involving tunnel fires walls must not be adiabatic,
• additional flow drag caused by tunnel infrastructure can be replaced by the equivalent wall roughness equal to 0.1 m,
• the fan model of A type (with given unit volume force) is a reasonable solution in most cases modeled in ANSYS Fluent,
• all examined turbulence models give almost similar and reliable results for fully developed flows, however the RSM model is the most accurate and the LES model will not provide an average steady picture of a flow,
• in some situations, when phenomena close to a fan outlet or complex flows are considered a deep preliminary analysis is recommended to avoid inappropriate turbulence model selection,
• the influence of the external wind and the stack effect on flows in tunnels is significant and can be modeled by an additional pressure at the windward portal, however determination of the exact value of this pressure is not an obvious issue,
• it is recommended to monitor the convergence process to avoid its premature termination (this concerns mainly steady calculations in ANSYS Fluent).

Basic level issues such as mesh density determination or mesh sensitivity analysis are omitted in the paper. One should have in mind that from the first moment when building a numerical model every assumption may have an impact on modeling accuracy. Hence, the authors hope the presented information will appear helpful for those, who deal with modeling of flows and fires in tunnels.