# Steps Towards a Real Time Solution of Fire in Tunnels

Bernhard Schrefler, Francesco Pesavento, Francisco Chinesta, Adrien Leygue

► To cite this version:

Bernhard Schrefler, Francesco Pesavento, Francisco Chinesta, Adrien Leygue. Steps Towards a Real Time Solution of Fire in Tunnels. 2nd International Conference on Computational Methods for Thermal Problems, Sep 2011, Dalian, China. hal-01008952

HAL Id: hal-01008952

https://hal.science/hal-01008952

Submitted on 3 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Steps towards a real time solution of fire in tunnels

B.A Schrefler, F. Pesavento

Department of Structural and Transportation Engineering, University of Padua, Italy

F. Chinesta, and A. Leygue

Ecole Centrale, Nantes, France

## ABSTRACT

We show a solution for a realistic tunnel fire simulation using a 2D-3D coupling strategy and a full multiphysics model for concrete (2D sections), known also as three fluids model (water, vapour, dry air). The computing times are too high for a real time solution. We investigate the applicability of the proper generalized decomposition (PGD) to obtain two goals: achieve a full 3D solution for the solid domain and reduce simultaneously the computing time. The first results are encouraging.

## 1. 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 to 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, see Figure 1. 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 one day on a PC.

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 of an extremely fast equation solver which can achieve a speed-up of up to 3600 times [4] and the adoption of the Proper Generalized Decomposition PGD [5,6] for a fast 3D solution of the problem of heated concrete. Steps in this direction as well as the general model will be shown.
## 2. Full 3D solution of heated concrete with PGD

In the case of the complex behaviour of the tunnel vault (transient and non-linear coupled multi-physic models) we would like to avoid the above mentioned 3D to 2D dimensionality reduction. Because of the richness of the thickness description due to many coupled physics with strong and fast evolutions in the thickness direction, a full 3D descriptions may involve millions of degrees of freedom in the solid domain that should be solved many times because of the history dependent chemo-hydro-thermo-mechanical behavior.

Today, the solution of such fully 3D models remains intractable for large size problems despite the impressive progress made in mechanical modeling, numerical analysis, discretization techniques and computer science during the last decade. New numerical techniques are needed for approaching such complex scenarios, able to proceed to the solution of fully 3D multiphysics models in geometrically complex parts. The well established mesh-based discretization techniques fail because of the excessive number of degrees of freedom involved in the fully 3D discretizations where very fine meshes are required in the thickness direction (despite its reduced dimension) and also in the in-plane directions to avoid too distorted meshes. A way to solve this problem is the adoption of the Proper Generalized Decomposition (PGD).

In what follows we illustrate the construction of the Proper Generalized Decomposition of a model defined in tunnel shell $\Xi = \Omega \times I$ with $\mathbf{x} = (x, y) \in \Omega \subset \mathbb{R}^2$ (the tunnel transverse section) and $z \in I = [0, L]$ its axis. We consider a generic transient partial differential equation, linear for the sake of simplicity:

$$ \frac{\partial u}{\partial t} - \nabla \cdot (\mathbf{K} \cdot \nabla u) + S = 0 $$

whose weak form writes:

$$ \int_{\Xi \times I} u^* \cdot \left( \frac{\partial u}{\partial t} - \nabla \cdot (\mathbf{K} \cdot \nabla u) - S \right) d\Xi \cdot dt = 0 $$

with $T$ the time interval in which the model is defined.

The solution $u(x, y, z, t) = u(\mathbf{x}, z, t)$ is searched under the separated form:

$$ u(\mathbf{x}, z, t) \approx \sum_{i=1}^{N} X_i(\mathbf{x}) Z_i(z) \Theta_i(t) $$

In what follows we are illustrating the construction of one such decomposition. For this purpose we assume that at iteration $n < N$ the solution $u^n(\mathbf{x}, z, t)$ is already known:

$$ u^n(\mathbf{x}, z, t) = \sum_{i=1}^{n} X_i(\mathbf{x}) Z_i(z) \Theta_i(t) $$

and that at the present iteration we look for the solution enrichment $R(\mathbf{x}) \cdot S(z) \cdot \theta(t)$:

$$ u^{n+1}(\mathbf{x}, z, t) = u^n(\mathbf{x}, z, t) + R(\mathbf{x}) \cdot S(z) \cdot \theta(t) $$

The test function involved in the weak form is searched under the form:

$$ u^*(\mathbf{x}, z, t) = R^*(\mathbf{x}) \cdot S(z) \cdot \theta(t) + R(\mathbf{x}) \cdot S^*(z) \cdot \theta(t) + R(\mathbf{x}) \cdot S(z) \cdot \theta^*(t) $$
By introducing the trial and test function in the weak form of the problem, we obtain the problem to be solved for computing the functions involved in the enrichment $(R(\mathbf{x}), S(z), \theta(t))$. Now, as the enrichment process is non-linear we propose to search the three functions $(R(\mathbf{x}), S(z), \theta(t))$ by applying an alternating direction fixed point algorithm. Thus, assuming $S(z)$ and $\theta(t)$ known, we compute $R(\mathbf{x})$, then we update $S(z)$ from $\theta(t)$ and the just updated $R(\mathbf{x})$. Finally we update $\theta(t)$ from the just computed $R(\mathbf{x})$ and $S(z)$. The process continues until reaching convergence. The converged solutions allow defining the next term in the finite sums decomposition: $R(\mathbf{x}) \to X_{n+1}(\mathbf{x}) S(z) \to Z_{n+1}(z)$ and $\theta(t) \to \Theta_{n+1}(t)$.

We can notice that this procedure involves a 2D solution for computing functions $R(\mathbf{x})$, and two 1D solutions for computing functions $S(z)$ and $\theta(t)$. Moreover, we can notice that the resulting algorithm is non-incremental because at each iteration we are looking the whole field history.

In the following we show an application of the 2D-3D coupling strategy for a 20MW fire in a tunnel and a PGD solution for the heat transfer problem only.

### 3. 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 Fig.s 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 $\mu=1.8 \times 10^{-5}$ kg/ms, specific heat $c_p=1006$ J/kgK, thermal conductivity $\lambda_f=0.026$, W/mK, density $\rho=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 \times 10^{-18}$ m$^2$, solid density $\rho_s=2564$ kg/m$^3$, solid thermal conductivity $\lambda_s=1.92$ W/mK, solid specific heat $c_{ps}=855.52$ J/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.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Global geometry of the tunnel.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2. Geometry of the tunnel: (A) 3D fluid domain, (B) 2D solid domain (i.e. tunnel cross section)</figcaption></figure>

The total thermal power involved by the fire is increasing in 10 minutes up to 20MW 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 $\alpha_c=500 \text{ W/m}^2\text{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.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 $\alpha_c$. As far as the mass exchange between the surface of concrete and the surrounding environment is concerned, a water vapor 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 \text{ Pa}$, $p^c=7 \times 10^{-7} \text{ Pa}$, $T=298.15 \text{ K}$. 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.

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 \text{ 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.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Velocity field (m/s) at $t= 600$ s (C), (A) cross-section and (B) longitudinal section.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Temperature distribution (K) at $t = 3600$ s in the fluid domain and in sections S2, S3 (fire), S4.</figcaption></figure>

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.
<figure><img src="image_6.png" /></figure>

A

<figure><img src="image_7.png" /></figure>

B

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Temperature and rel. humidity distribution (A) and gas pressure and total damage distribution (B) in the top of the central section S3 (fire) at $t = 3600$ s.</figcaption></figure>

A

<figure><img src="image_9.png" /></figure>

B

<caption>**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).</caption>

**4. Heat transfer simulation in the tunnel vault with PGD**

In this section we present some preliminary results of the simulation of the tunnel fire described above obtained by applying the Proper Generalized Decomposition techniques illustrated in section 2. The main aim of this calculation is to demonstrate the validity of the approach considering, as initial stage of the on-going research, a simplified model limited to heat transfer through conduction.

The mesh used in the same as in the previous example in the transversal section, but using nine-noded elements., while in the longitudinal direction 8000 elements have been used. The duration is 1 hour with a time step of 1 sec.
As already pointed out the problem is solved assuming some simplifications as far as both the physical problem and the boundary conditions are concerned. The concrete vault is treated as a solid domain for which only the pure conductive thermal problem is considered, neglecting convective transfer of energy, phase changes and the mass transport inside the pores of the material. The boundary conditions consist in the definition of a stationary normalized thermal flux on the inner surface of the tunnel vault according to the following formula:

$$norm. heat flux = F_x(x, y) \times F_z(z)$$

in which $F_x(x, y) = y$ and the distribution function along the z-coordinate is shown in Figure 7B:

<figure><img src="image_10.png" /><figcaption>Figure 7. Sketch of the mesh (A) and of the bc.s (B) used in the case of tunnel fire solved by PGD.</figcaption></figure>

The whole 20 MW power is applied directly to the heated surface without taking in to account the role played by air, i.e. by using Neumann condition instead of Robin condition defined trough a heat exchange coefficient $\alpha$ between the surface of concrete structure and the surrounding fluids. Because of this we expect a large difference in terms of temperature field between the solution obtained by means of PGD and the original one obtained by using the full model. The results of PGD analysis are depicted in Figure 8.

<figure><img src="image_11.png" /><figcaption>Figure 8. Temperature distribution along z-axis (A) and in the thickness of the central cross section of the tunnel (B) obtained from PGD analysis.</figcaption></figure>

The results shown in Figure 8 have been obtained by using 13 modes and took 2 min and 30 sec of computational time on a standard INTEL cpu based computer using MATLAB. A fast check has shown that the computed stored energy in the concrete vault after one hour equals the heat input. The time span is too short for heat flux to reach the boundary to the surrounding rocks, see the temperature gradient in Figure 8B.
The surface temperatures reached are comparable to the ones obtained with the full model, [1], in the case of pure CFD simulations for which the solid domain is neglected together with the related heat exchange fluxes. In that case the temperature values were much higher than those obtained with the model of section 3 and approximately equal to 4000 K [1], that is the same order of magnitude of thermal field attained by using PGD in this work.

**5. Conclusions**

We have shown a 3D-2D coupled solution for a tunnel fire simulation, taking the three fluids concrete model into account. For a fire of a duration of one hour the calculation time on a PC is well over one day. For the thermal problem in the tunnel vault a PGD approach on a full 3D solid domain model has been shown. The computing time was reduces to two and half minutes.

Considering the short computational time and the fact that it is possible to obtain physically reasonable results just by introducing in PGD analysis more sophisticated physical models and physically correct bc.s (i.e. taking into account the role of the air), this set of preliminary computations are extremely encouraging.

**References**

1. B.A. Schrefler, R. Codina, F. Pesavento, J. Principe, (2010), Thermal coupling of fluid flow and structural response of a tunnel induced by fire, *IJNME*, DOI 10.1002/nme.3077.

2. D. Gawin, F. Pesavento, B.A. Schrefler, (2011) What physical phenomena can be neglected when modelling concrete at high temperature? A comparative study: Part 1: physical phenomena and mathematical model. Int. J. Solids Struct., doi:10.1016/j.ijsolstr.2011.03.004.

3. D. Gawin, F. Pesavento, B.A. Schrefler, (2011) What physical phenomena can be neglected when modelling concrete at high temperature? A comparative study: Part 2: comparison between models. Int. J. Solids Struct., doi:10.1016/j.ijsolstr.2011.03.003.

4. A. Bertoldo, M. Bianco, and G. Pucci, (2004) A fast multifrontal solver for non-linear multi-physics problems, *Computational Science, ICCS 2004*, PT 2, Book Series: *Lecture Notes in Computer Science* Vol. 3037, 614-617.

5. D. Dureisseix, P. Ladeveze, B.A. Schrefler (2003), A computational strategy for multiphysics problems – application to poroelasticity, *Int. J. Numer.Methods Engng.*, 56, 1489-1510.

6. F. Chinesta, A. Ammar, E. Cueto (2010), Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models, *Archives of Computational Methods in Engineering*, 17, 327-350.

7. Brite EURAM Research Project “HITECO - High Temperature Concrete” (1999), final report.