Fully coupled numerical simulation of fire in tunnels : From fire scenario to structural response

In this paper we present an efficient tool for simulation of a fire scenario in a tunnel. The strategy adopted is based on a 3D-2D coupling technique between the fluid domain and the solid one. So, the thermally driven CFD part is solved in a three dimensional cavity i.e. the tunnel, and the concrete part is solved on 2D sections normal to the tunnel axis, at appropriate intervals. The heat flux and temperature values, which serve as coupling terms between the fluid and the structural problem, are interpolated between the sections. Between the solid and the fluid domain an interface layer is created for the calculation of the heat flux exchange based on a “wall law”. In the analysis of the concrete structures, concrete is treated as a multiphase porous material. Some examples of application of this fully coupled tool will be shown.


INTRODUCTION
The availability of an efficient tool for simulation of a fire scenario in a tunnel is of paramount importance for fire safety management in emergency situations, for training of fire brigades prior to emergency cases in order to be able take the right decisions when needed and to evaluate measures geared to increase the resistance of existing tunnel vaults against spalling.We have developed such a tool, which takes the thermal fluid-structural coupling in a tunnel fire fully into account [1].It appears as one of the largest coupled problems actually solved in the community of computational interaction problems.
The simulation of a realistic fire scenario is still a time consuming task and the tool is not yet completely ready for the first of the above mentioned three goals.One of the bottlenecks is the heavy computational burden linked with the three fluids model for concrete.
It is not possible to disregard the enormous heat sink the tunnel vault represents with the phase changes and chemical reactions going on in heated concrete.Such an omission can yield temperature fields also some 1000 • C above measured ones in an experiment.On the other hand simplifications of these phenomena are not possible as highlighted in two recent companion papers [2,3].In the existing model [1] we have chosen a 3D-2D coupling strategy where the thermally driven CFD part is solved in a three dimensional cavity i.e. the tunnel, and the concrete part is solved on 2D sections normal to the tunnel axis, at appropriate intervals.The heat flux and temperature values, which serve as coupling terms between the fluid and the structural problem, are interpolated between the sections.With this approximation the heat transfer in the tunnel vault in the direction of the tunnel axis is disregarded.As an example, with such an approach the fully coupled simulation on a realistic tunnel for a fire of 20 MW of the duration of one hour lasts more than three days on a high-end workstation.
The aim of our current research effort is twofold: realize a true 3D-3D coupling on one hand and reduce drastically the computing time on the other hand.The way for achieving this is through adoption MATEC Web of Conferences of an extremely fast equation solver which can achieve a speed-up of up to 3600 times [4].At the same time the on-going research aims also at the assessment of various methods to limit the consequences of a fire in the tunnel.In particular in the framework of HYDROFIRESHOCK project the research team is evaluating the effectiveness of an innovative fire suppression system.

INTERFACE CONDITIONS AND COUPLING STRATEGY
To solve the thermo-mechanical problem on the domain we consider a geometric domain decomposition of the problem by means of a non-overlapping subdomain approach.We split the domain into the fluid and structural subdomains as = S ∪ F .The conditions to be satisfied at the interface are the continuity of the temperatures and velocities as well as the normal components of heat fluxes and tractions.If we denote by SF (t) the surface of the tunnel walls, these conditions read where it is explicitly noted that SF (t) is actually a moving surface and its position depends on the solution of the structural problem.The approximations proposed in [1,5] are based on two assumptions: 1.The velocities of the solid medium are small (compared to the dimensions of the tunnel).
2. The mechanical traction produced by the fluid on the solid is small.
These assumptions result in uS ≈ 0 and n • F | SF ≈ 0 and therefore the boundary SF (t) remains fixed and the interface conditions become: ( Thus the coupling between the solid and the fluid is due to the thermal problem only.A further approximation was developed in [5] to consider the problem of the strong boundary layers present in a turbulent flow.This second approximation is based on a non-overlapping domain decomposition of the problem in three subdomains, one in the solid region and two in the fluid region.One of the fluid subdomains will be a thin region of thickness near the solid surface and the other will be the rest of the fluid domain.Using the well-known wall function approach [6], the problem in boundary layer is approximately solved and an iteration strategy between the remaining subdomains is proposed.Therefore this second approach also involves two subdomains but the interface conditions now read: and where t is now the tangential stress acting on the fluid and and c are parameters that depend on the wall function coefficients.The first one is computed from the universal velocity profile assumed in the wall function approach, as usual [6].The second one could also be computed in a similar fashion but it is estimated from experimental values in this work.It is also important to remark that we do not consider any mass exchange between the solid and the fluid.Considering a geometric domain decomposition of the problem by means of a non-overlapping subdomain approach, we expect to construct the solution problem from the solution of local problems for the fluid and the structure using the interface conditions above described.This is carried out by iteratively solving local problems on each domain until convergence on the interface conditions is satisfied, that is to say we use an iteration-by-subdomain strategy.The description of the local problems for the concrete structure (i.e. the solid domain) and the fluid domain is beyond the scope of this paper.In the following we just describe briefly the coupling strategy between the local problems.
If L denotes the differential operator on the domain and B the differential operator on the boundary the situation is as follows.On the fluid region we solve at the iteration i of each time step: where boundary conditions on SF depend on the solution of the structural problem at the previous iterate.On the structural domain we solve at the i-th iteration where we can take k = i or k = i−1.In the first case the solution of this problems is sequential whereas in the second one it can be parallel.The choice of the boundary conditions of the local problems should be such that presented interface conditions are satisfied when convergence is achieved.It is well known from the theory of domain decomposition methods that in the case of non-overlapping subdomains we can choose Dirichlet-Neumann(Robin); Neumann(Robin)-Dirichlet or Robin-Robin.In this work we apply the interface conditions already described.One important point of this strategy is that we already have programs that solve the fluid dynamics problem and the structural problem.Then a master/slave algorithm was implemented by developing a third code (the master code) in order to control the iterative process.The MPICH2 library, an implementation of the MPI-2 standard, provides functions for process communications that are used to interchange the data needed to apply boundary conditions on each dedicated (slave) code.Some minor modifications on these codes are needed in order to exchange data with the master.In order to perform a calculation, input data for each subproblem needs to be generated and the master code starts the calculation by starting the slave process (this is only possible under MPI-2 standard).During the calculation the master code needs to define the boundary conditions to be applied on each subproblem.
A particular aspect of this implementation is the need of coupling a two dimensional code with a three dimensional one.On the one hand, due to the high number of state variables of the structural model, only two dimensional calculations are performed using the structural code.On the other hand three dimensional calculations can be performed using the fluid dynamics code.Therefore, the three dimensional concrete tunnel vault is approximated by several two dimensional sections and variables are linearly interpolated to generate boundary conditions on the fluid.This task is also performed by the master code.The values of the temperature or heat flux to be applied as boundary condition on an interface node need to be calculated from the temperature on the other subdomain.This is done by finding the host element and interpolating.The element search strategy used is based on the quad-tree algorithm.

Tunnel fire: Example of the 3D-2D coupling strategy
The structure under consideration is the tunnel of Virgolo close to Bolzano (Italy) that has been also used for an experimental test in the framework of UPTUN project [7].We have considered the central part of the tunnel, 80 m long.Its geometry is decomposed in the fluid and the solid domains, see Figs. 1,2.The solid domain consists in the cross section of the tunnel vault.In the simulations five cross sections are considered at 0, 30, 40, 50, 80 meters along the longitudinal axis z.The location of fire is the section at 40 m.The fluid is considered as an ideal gas and has the following properties: dynamic viscosity = 1.8 × 10 −5 kg/ms, specific heat c p = 1006 J/kgK, thermal conductivity f = 0.026, W/mK, density = 1.225 kg/m 3 .Concrete used for the solid domain (i.e. the sections of the concrete tunnel) is C60 concrete (with a final compressive strength equal to 60 MPa) and has the following main properties at ambient temperature (20 • C): elastic modulus E = 40000 MPa, porosity n = 0.082, intrinsic permeability k = 2 × 10 −18 m 2 , solid density s = 2564 kg/m 3 , solid thermal conductivity s = 1.92 W/mK, solid specific heat c ps = 855.52J/kgK.The volumetric heat source corresponding to the fire is located at the coordinates (x,y,z) = (1.0,0.5, 40.0) and has a volume equal to 8 m 3 .This means that the fire is located in the central section of the tunnel at 0.5 m from the longitudinal axis and at a height of 1 m from the road pavement.
The total thermal power involved by the fire is increasing in 10 minutes up to 20 MW following a linear law and then kept constant.For this analysis 15300 hexaedral elements are used in the fluid domain, while each cross section is discretized with 640 quadrilateral elements with eight nodes.The initial and the boundary conditions selected for this case are: • For the fluid domain: the atmospheric pressure is imposed at the ends of the tunnel, the initial fluid velocity is equal to zero, close to the tunnel vault the fluid can exchange heat with the concrete structure surface according to the universal profiles ("wall law") described in [1], with a heat exchange coefficient equal to c = 500 W/m 2 K.The initial temperature is set to 298.15 K for the whole fluid domain.• For the solid domain: on the inner side of the cross section, i.The total time of simulation is 1 hour.The case under consideration corresponds to a real fire case in terms of the total thermal power involved, the duration of the fire and the value of the heat exchange coefficient selected.Figure 3 shows that the velocity of the ascensional flux close to the fire source is higher than 9 m/s, while the horizontal fluxes flowing toward the ends of the tunnel have a velocity equal to 3 m/s (t = 600 s).The temperature distribution in the fluid domain and in the top part of the sections S2, S3 (the section of the fire) and S4 are shown in Figure 4.The central section, that is the most stressed one, is the most exposed to spalling risk.Indeed, the peak of gas pressure (2.4 MPa) in the external layer of the concrete vault (10 cm thick), the formation of the "moisture-clog" (the relative humidity reaches in this zone a value higher than 90%) and, finally, the total damage ( ∼ = 55%) can lead to a progressive spalling starting from that layer, see Fig. 5. Finally, Figure 6 shows the distribution of maximum temperature in the cross sections along the longitudinal axis and a comparison between the heat source temperature evolution in time and the temperature profiles most used in literature.The case considered corresponds to a fire with a temperature profile between an ISO-Fire and a Hydro-carbon Fire.

FIRE MITIGATION STRATEGIES
In this section we present an application of the FIT tool described in section 2 within the HYDROFIRESHOCK Project.The project aims at the development, testing and deployment of an innovative automatic system for the suppression of fire in tunnels.The system consists of (see Fig. 7): • intervention and control stations at regular intervals (typically every 42 m); • each control station is equipped with two cameras in the visible range and two in the I.R. (Infra Red) range; • the control stations are interconnected via a pipe under pressure ( 10 bar) for supply of water or anti-fire foam; • mobile structures (i.e.overhead trolley located every 800 m) connected to the control/docking stations (42 m); • each mobile structures has two extinguishing devices (remotely controlled) with a capacity of 1000 lt/min.
The work is in progress and the numerical tool will be employed for the analysis of the effectiveness of the extinguishing system.In particular, various configurations of the system will be studied, with
e. the vault surface in contact with the fluid, two convective (i.e.Robin) boundary conditions are imposed.The convective heat exchange is governed by the same universal profiles described for the fluid domain with the same exchange coefficient c .As far as the mass exchange between the surface of concrete and the surrounding environment is concerned, a water vapour pressure equal to 1300 Pa and an exchange coefficient of 0.02 m/s are set.The initial condition for the concrete structure are p g = 101325 Pa, p c = 7 × 10 −7 Pa, T =298.15K.This set of values corresponds to an initial relative humidity equal to 58%.On the outer side of the cross sections the values of gas pressure, capillary pressure and temperature are fixed (i.e.Dirichlet bc.s) to the initial ones.

Figure 6 .BFigure 7 .
Figure 6.Distribution of the maximum temperature of the cross-sections along the z-axis (A) and comparison of the source temperature evolution with the main heating profiles available in the literature (B).