Thermal analysis of two-dimensional structures in fire

Resumo The structural materials, as reinforced concrete, steel, wood and aluminum, when heated have their mechanical proprieties degraded. In fire, the structures are subject to elevated temperatures and consequently the load capacity of the structural elements is reduced. The Brazilian and European standards show the minimal dimensions for the structural elements had an adequate bearing capacity in fire. However, several structural checks are not contemplated in methods provided by the standards. In these situations, the knowledge of the temperature distributions inside of structural elements as function of time of exposition is required. The aim of this paper is present software developed by the authors called ATERM. The software performs the thermal transient analysis of two-dimensional structures. The structure may be formed of any material and heating is provided by means of a curve of temperature versus time. The data input and the visualization of the results is performed thought the GiD software. Several examples are compared with software Super TempCalc and ANSYS. Some conclusions and recommendations about the thermal analysis are presented.

Thermal analysis of two-dimensional structures in fire

Introduction
In several engineering areas, it is very important to know the temperature distribution within a solid body.In civil engineering, especially in the structures fire design and in the design of roll compacted concrete dams, the interest in the study of temperature field evolution within structural elements have significantly grown in recent years.Concrete structure fire design is ruled by ABNT NBR 15200:2012 [1].This standard presents simplified methods and allows using advanced methods for checking structural elements in fire, besides the use of tests.The simplified methods are applicable in several situations.Nonetheless, there are design situations that cannot be represented by them, asking for a more refined analysis.The use of advanced methods for verifying the fire safety of structures requires models for determining the development and distribution of temperatures in structural parts.The thermal analysis may be based on the heat transfer principles and hypotheses.The model adopted may take into consideration the relevant thermal actions and the variation in thermal properties.When relevant, the effects of non-uniform thermal exposition and heat transfer through components of adjacent constructions should be included.This paper presents a formulation for the thermal analysis of twodimensional structures using finite elements.This formulation was implemented in a piece of software named ATERM [2], developed by the authors using FORTRAN 90 language.For data input and visualization of the field of temperatures, GiD [3] software is used.For the validation of results, comparisons were made using Super Tempcalc (STC) [4], developed in Lund, Sweden.STC began being developed in 1985 for the thermal analysis of two-dimensional structures.The software was checked against several experimental tests since its first version and was used for elaborating Eurocode 2 part 1.2 [5].The results obtained by means of the ATERM [2] software will also be compared to that obtained by ANSYS [6].SOLID70 elements will be used in the models for the thermal analysis.Convection and radiation effects will be simulated by means of the SURF152 element using one extra node.Models constructed in ANSYS software are tridimensional, but they can be compared to the two-dimensional models obtained by ATERM because there is no temperature variation along the longitudinal axis of the structural element.ATERM was also used by SILVA [7] in the study of reinforced concrete columns in fire.In that study, a collection of isotherms for several cross sections was created.

Finite elements method
In its beginning, the Finite Elements Method was applied to mechanics problems in structural engineering.The thermal analysis was the first non-structural area to use it for modelling engineering problems.The determination of the field of temperatures within the structural element is the first step for structure fire analysis.Temperature distribution affects the stress distribution within the structural element.Heat propagates within concrete by conduction.This phenomenon is governed by Poisson's equation, which, in its two-dimensional form, is given by equation (1). (1)

&
In equation (1), λ is the material thermal conductivity, Q  is the internally generated heat per volume unit per time, ρ is the den- sity, c is the specific heat, θ is the temperature and t is the time.
To solve equation (1), it is necessary to impose the boundary and initial conditions of the mathematical model.(2) One adiabatic, or thermally insulated, surface can be simulated by imposing zero heat flow ( The convection and radiation phenomena are included into the numerical model by means of Neumann's boundary condition, which can be written by equation (3).
(3) ( ) ( ) α is the convection coefficient, θ and θ ∞ are, respectively, the temperatures in the structure and out of it, ε is the emissivity, σ is the Stefan-Boltzmann constant, T and T ∞ are, respectively, the absolute temperatures in the structure and out of it.By making the part due to radiation linear, equation (3) can be rewritten by means of equation ( 4), (4) ( ) ( ) ( ) c r q a q q a q q a q q In equation (4), r α is the coefficient for heat transfer by radiation provided by equation (5).
(5) As the coefficient for heat transfer by radiation depends on the structure temperature, the process becomes iterative.

( )( )
The solution for the differential equation (1) can be numerically obtained by the finite elements method (FEM).For the use of FEM, it is necessary to write the weak problem formulation, which is obtained by equation ( 6). (6) In equation ( 6), w is an arbitrary function named weight function and Ω is the problem dominium.
Integrating the two first terms of equation ( 6) by parts and applying Gauss' Theorem, we obtain equation ( 7). (7) In equation ( 7), S represents the problem boundary region and l and m are the director cosines.Substituting equation ( 7) in (6), we have equation ( 8). ( Discretizing the problem dominium in a finite number of elements and using Galerkin Method, the set of interpolation functions N, as weight functions, the temperature in any point of the finite element can be approximated by equation ( 9).
(9) ( ) , , In equation (9), n is the number of nodes of the element, ( ) are the interpolation functions and i θ are the nodal temperatures of the element.Substituting equation ( 9) in ( 8), we get to the equations system (10), which represents the thermal balance in each finite element. (10) In equations (10), A and S represent the area and the surface of the element dominium.
The derivative of the temperature field in relation to the normal to the boundary surface, n, can be written in function of the director cosines, l and m, by means of equation ( 11). ( Substituting equations (2), ( 4) and (11) in equation ( 10), we obtain expression (12), (12) & Defining the interpolation function vector N and the nodal temperatures vector θ by means of expressions (13). (13 Where n is the number of nodes of the element and N i are the interpolation functions of the finite element.Thus, we can rewrite equation (9) in the matrix form by means of equation ( 14). (

Nq
We also define the conductivity matrix of the material, for the case of isotropic material, by means of equation ( 15) and the gradient vector by equation ( 16).
Thermal analysis of two-dimensional structures in fire (15) 13) and ( 16) in equation ( 12) and rearranging the terms, we obtain the following set of equations ( 17).
The set of equations ( 17) can be rewritten in the matrix form by means of equation ( 18). (18) In equation ( 18),  θ is the first derivative of the temperatures field in relation to time.
The total capacitance matrix ( K ) of the element is given by equation (19). ( The element conductivity matrix is given by equation ( 20).We define the thermal capacity of the material as the product density ( ρ ) x specific heat ( c ).Thus, the thermal capacity matrix ( C ) of the element is given by equation ( 23). ( The vector of the consistent thermal actions ( F ) is given by expression (24). ( In equation ( 24), θ ∞ is the temperature outside the structure.
Matrixes N and B depend on the type of element used in struc- ture discretizing.
In ATERM, we employ rectangular and triangular planes of four and three nodes, respectively, and two nodes special bar elements for considering convection and radiation heat transfer effects, which can be superimposed in any face of the plane elements.These elements are briefly described in the following items.Further details on the used finite elements and the computational implementation can be found in Pierin [2].

Four nodes rectangular element
The four-node plane rectangular finite element with sides 2a by 2b used herein is represented in Figure 1.The interpolation functions i N are defined by means of equation ( 25). (25) By substituting equations (15), ( 21) and (25) in equation (20), we obtain equation (26), which represents the matrix of conductivity of the four-node plane rectangular element. (26)

Triangular element
The triangular finite element used in ATERM software has three nodes identified by 1, 2 and 3 and is represented in Figure 2. (28) ( )

N a b x c y i A
In equation (28), A represents the area of the finite element and the coefficients i a , i b and i c are provided by expressions (29). (29) In expressions (29), i X and i Y represent the coordinates of node i.

Two nodes bar special element
The Neumann's boundary conditions can be imposed to the finite elements model by means of a two-node linear element with (32)

N
Note that the convection linear element is fully compatible with the previously formulated four-node rectangular element because it uses linear functions for temperatures interpolation along the element side.Substituting equation (32) in expression ( 22), the convection matrix and the vector of consistent thermal actions due to the convective heat flow are provided by equations ( 33) and (34), respectively,

Cavities
In some structural elements, such as alveolar, ceramic blocks and ribbed slabs filled with EPS blocks, there is the presence of enclosed air inside cavities.The existence of cavities within the structure generates heat transfer by convection and radiation, due to the heating of the enclosed air.
Discretizing the analyzed structure boundary in plane finite elements and the cavity contour in special two-node finite elements, as per Figure 4, we can obtain the temperature of the air inside the cavity from the nodal temperatures by means of equation ( 35). ( In equation ( 35

Pre-processing
The data input to ATERM software is accomplished by text files that contain information on the two-dimensional structure to be thermally analyzed.This information comprehends the nodes coordinates, the elements connectivity and thermal properties of the materials.
In order to ease the creation of the input data file, some geometry models, such as rectangular and T sections, were created based on GiD software [3], as shown in Figure 5.In these models, it is only necessary to inform the dimensions of the sections, of the finite elements and the type of material.
The exposition to fire can be configured directly in the model, by the indication of the surfaces exposed to standard fire and the surfaces exposed to a constant temperature, as shown in Figure 6.In addition, emissivity and convection coefficients for the consideration of heat transfer by convection and radiation should be informed.
The thermal properties of the materials, such as conductivity, density and coefficient of strength reduction, are informed by means of tables in function of the temperature as shown in Figure 7.We can also define materials the thermal properties of which are independent of temperature.
Other data also necessary to the thermal analysis, such as total time of fire exposition, time increment, initial temperature of the structure, α coefficient for the choice of the time integration method, tolerance for temperatures convergence, time increase for printing results, should be informed as shown in Figure 8. Coefficient α varies from 0 (Euler Implicit Method) to 1 (Euler Explicit Method) and defines the numerical stability for the integration along time.STC software uses α=2/3, known as Galerkin Method.

Numerical simulation
In order to validate the software developed in this work, three simulations were performed, the results of which were compared to those obtained by means of STC [4] and ANSYS [6] software.
The thermal properties adopted for concrete were based on ABNT NBR 15200.2012[1].In the examples herein, we used siliceous concrete humidity content equal to 1.5% in weight.In all of the examples, the thermal action is determined according to ISO 834 [8] standard fire curve.

Concrete beam
As example # 1, the heating of a 20 x 50 cm rectangular  were adopted equal to 0.7 and 25 W/m².ºC,respectively.For the initial temperature of the structure, the value 20ºC is adopted.
The structure was discretized in square elements.For ATERM software, 3 meshes were used: in Mesh 1, we used elements with 5 cm sides, in Mesh 2, with 2.5 cm sides and, in Mesh 3, with 1 cm sides.For modeling in ANSYS and STC, we only used Mesh 3.
Figure 10 shows the influence of the mesh refinement on the evolution of temperature in function of time in the node in the middle of the smaller heated side.With mesh refinement, the temperatures is observed to rise a little up to instant t = 28 minutes and after this instant, the temperatures obtained for the different meshes are virtually equal.The solution is verified to tend to converge with mesh refinement and that, for Mesh 3, the results from ATERM, ANSYS and STC are virtually equal.

Ribbed slab
Figure 11 shows the geometry of the thermally analyzed ribbed slab (dimensions in cm).The slab underside is exposed to fire.The adopted values for emissivity factor and convection coefficient were 0.7 and 25 W/m².ºC,respectively, as per Eurocode 2 part 1.2 [5] recommendations.In the upside of the slab, not exposed to fire, the combined phenomena of convection and ra- Thermal analysis of two-dimensional structures in fire diation were simulated by α c = 9 W/m² °C [10].The adopted initial structure temperature was 20 ºC.Due to symmetry, only half of the slab was modeled and the symmetry axis was considered to be adiabatic.
Figure 13 shows the variation in the nodal temperature in function of time obtained by ATERM and STC for three distinct points located in the slab symmetry axis (as shown in Figure 12): point Ain the base of the rib.Point B -at the level of the junction between the rib and the flange.Point C -in the upside of the slab.The slab was discretized by triangular elements with 0.5 cm sides, 1653 elements being generated in STC and 1809 elements in ATERM.The difference between the numbers of elements generated by the two pieces of software is due to the distinct mesh generators used by them.The results obtained from the two pieces of software are virtually equal.It is also noted that there is good correlation between the results from ATERM and ANSYS.
Figure 14 shows the field of temperatures after 60 minutes of fire, obtained by ATERM and STC.

Design of padded ribbed slabs
In Brazil, ribbed slabs filled with EPS on cement plates, as shown in Figure 15, are also used.When heated, EPS decomposes quickly, resulting in a large space filled with trapped air.This trapped air contributes to the heat transfer between the cement plate, which is at higher temperatures, and the flange.This heat transfer can be taken into account in the numerical model by means of the procedure indicated in item 3 herein.
In the face directly heated by the fire, we applied the ISO 834 [7] standard curve and, in the face not directly exposed to the heat, a combination of convection and radiation, simulated by a c = 9 W/ m² °C, was used.
The thermal properties of the cement plates used in the analyses were: thermal conductivity equal to 2.22 W/mºC, specific heat equal to 840 J/kg K and density equal to 1200 kg/m³ [11].The variation in temperature was analyzed in two points, as shown in Figure 16 The consideration of the trapped air in the model is observed to cause heating in the upside of the slab in the region between ribs.
The fields of temperature at 60 minutes of standard fire on the downside of the slab, with and without the trapped air consideration, are represented in Figures 18 and 19, respectively.Disregarding the air trapped inside the cavity does not generate the heating of the upside of the slab in the region between ribs, i.e., after 60 minutes of fire, the upside of the slab is at room temperature.There is only a slight rise of the upside temperature in the region of the ribs.
In reality, the trapped air generates heat transfer by convection and radiation between the heated surface (the cement plate) and the face not exposed to fire (flange).This example shows that this heat transfer can be modeled by means of cavities.

Conclusion
Within this work, we presented the formulation of the Finite Element Method applied to non-linear thermal analysis of two-dimensional structures, employed in the elaboration of the ATERM software.This kind of analysis is critical to the study of structures in fire.
The results were validated with STC [4] and by means of ANSYS [6] software.It was observed that the results obtained by both types of software present an excellent correlation in all the examples performed.However, depending on the mesh used, the results were observed not to correspond to the physical behavior.Thus, a study of geometrical and time meshes is recommended before undertaking a thermal analysis.
In some situations, such as in ribbed slabs and ribbed slabs filled with inert material, cavities do exist within the structural element.Such cavities are filled with trapped air that, when heated, transfer heat by convection and radiation to the neighboring regions.
In the finite element analysis, the heat transfer by means of the trapped air can be considered, both by ATERM and STC softwares.
The performance of ATERM showed that the software is very effective in relation to STC, having a much shorter processing time.By means of the GiD software, several sections were created for data input, such as rectangular and T sections, which facilitate the creation of the geometry and the mesh of finite elements.
And the matrix due to the combined effects of convection and radiation of the element is given by equation (22).

FunctionsiN
used for interpolation of the nodal temperatures within the triangular finite element are defined by equation (28).
Substituting equation (28) in equation (23), we obtain the matrix of thermal capacity of this element, expressed by equation (31).

FigureL
Figure 1 -Four nodes rectangular element c and ar ρ are the specific heat and density of air, respectively.i L is the length of the special finite element.Equation (35) can be rewritten by means of equation (36).
Figure 4 -Mesh of finite elements

Figure 5 -
Figure 5 -Models for geometry creation

Figure 8 -Figure 10 -Figure 11 -
Figure 8 -Other data Figure 12 -Analyzed model Figure 13 -Temperature evolution in function of time

Figure 15 -
Figure 15 -Model with EPS/cement plate (dimensions in mm) Figure 16 -Control points for the analysis

Figure 17 -
Figure 17 -Evolution of temperature in function of time for points A and B ), i α is the convection coefficient of element i,