Thermomechanical analysis of reinforced concrete columns exposed to fire

ABSTRACT: A three-dimensional (3D) transient numerical model for thermomechanical analysis developed with Finite Element Method (FEM) using the software ANSYS 19.1 is exposed in this paper. The proposed model aims to predict the structural behavior of reinforced concrete columns in a fire situation since it is known that high temperatures significantly reduce their mechanical resistance. For this, the main factors that govern their structural behavior should be considered. Analyses obtained by the proposed model were validated with results from experimental data, evidencing a good correlation between numerical and experimental fields.


INTRODUCTION
Over the years, mainly due to disasters caused by severe fires in buildings, industrial pavilions, tunnels and other civil engineering works, many research groups have been focusing on evaluating the behavior of structural elements exposed to fire.
Specifically, in Brazil, there are various researches conducted in this area, published through dissertations, theses, and several journal articles. To illustrate it, Silva [1] carried out his thesis evaluating steel structures exposed to fires and continued to contribute with many works on fire safety [2], [3]. Carla Costa, for instance, had important cooperation in designing reinforced concrete elements in fire situations [4], [5].
In the state of Rio Grande do Sul, Lemos [6], through his dissertation, created a computational code able to perform a numerical study of reinforced concrete structures under high temperatures by the finite element method. In addition, Kirchhof [7] developed during her Doctorate at the Structural Models and Testing Laboratory (LEME-UFRGS), a numerical and experimental study analyzing how the water content influence on explosive spalling phenomena in concretes exposed to elevated temperatures. Finally, Bolina [8] presented results of the experimental tests carried out at its Performance/UNISINOS, evaluating the influence of durability requirements on fire safety through prototypes of prefabricated reinforced concrete columns.
Moreover, many studies involving numerical models of structures exposed to fire using the software ANSYS, in the national scenario, was developed by researches of USP -São Carlos, including thermomechanical analyses in steel structures [9] and [10], concrete-filled steel columns [11], composite steel and concrete beams [12], composite concrete and timber beams [13], among others.
It is known that columns of a building structure are the main responsible for its integrity, hence the need to predict their structural behavior in harmful situations through analytical, numerical and experimental methods. Experimental testing is the best way to understand the behavior of these structures. However, laboratory procedures are often limited in quantity due to high costs and time demands. Thus, the development of numerical models, calibrated with experimental data, allows expanding the knowledge with more agility and resource-saving.

Objectives
The current study aims to evaluate, through an advanced calculation model using the finite element method (FEM), reinforced concrete columns behavior in fire situations considering geometric nonlinearity and the elastoplastic regime of the materials. This work was elaborated with the aid of ANSYS 19.1 software [14] where columns models were exposed to a standard temperature-fire curve recommended by ISO 834 [15].
The expectation at the end of the study, by validating the thermomechanical model created with consistent experimental data present in the bibliography, is to conclude that the proposed model can offer good responses and enables several applications and analyses for future researches involving this thematic.

METHODOLOGY -PROPOSED NUMERICAL MODEL
In this work, a three-dimensional numerical model was created using ANSYS 19.1 software [14] to carry out transient thermomechanical analyses, that is, evaluations considering the variability of the parameters in accordance with temperature increase. In order to do this, the coupling of two engineering analyses fields was necessary: thermal and mechanical.
It must be observed that the processing time of all models presented in this work will be indicated in their respective subitems related to results. Besides, all simulations were performed with the aid of an Intel Xeon E3-1225 v5 3.30 GHz Processor, 8.00 GB installed memory, and Intel HD Graphics P530.

Thermal analysis
Firstly, through purely thermal analysis, was achieved the thermal gradient evolution in the structure. Such a study was performed by exposing columns to the standard fire curve ISO 834 [15]. This curve, expressed according to Equation 1, represents a typical building fire whose combustible material considered is cellulosic. Figure 1 shows the evolution of temperature as a function of time in accordance with the curve mentioned above.   [19], such consideration is attributed to the fact that steel has a stable and well-defined microstructure at elevated temperatures. For concrete, in turn, the Annex C of ABNT NBR 15200:2012 [17] points out that specific mass is influenced by water loss; thus, it varies with increasing temperature.

Discretization of the structure into finite elements
In order to simulate the thermal response, reinforced concrete columns are discretized using three different types of elements available in ANSYS 19.1 internal library [14]. To represent solid materials, the elements SOLID70 and LINK33 were used. To enable the application of fire boundary conditions, SURF152 was used overlaid onto a heated surface of 3D thermal solid elements.
SOLID70 is an eight nodes element used for three-dimensional applications and LINK33 is a bar element with only two nodes, both having the temperature as a single degree of freedom, at each node. It allows the element capability to simulate thermal conductivity, being the solid element applicable to concrete simulation and the bar element applicable to steel. SURF152, in its turn, is a thermal surface effect element that is defined by four to nine nodes with the possibility of using an extra node. The function of the extra node, located away from the base element, is to be a logical point for the fire curve application. The objective of surface elements is to simulate heat transfers by convection and radiation. For that, material properties such as emissivity and convection heat transfer coefficient should be assigned.
Thus, with the mentioned discretization, it is possible to encompass the three main heat transfer mechanisms. Fire exposed surfaces contribute to the temperature rise by convection and radiation through the SURF152 element, while internal temperature conduction capability is attributed by SOLID70 and LINK33 elements.
In addition, depending on the structural arrangement of each model, if possible, the principle of symmetry is used, reducing computational costs and analysis time. For better understanding, Figure 2 presents a generic model created as described above.
The thermal analysis results consist of determining temperatures at each node of the mesh. It is also noteworthy that, according to the analyses realized, the temperature remains uniform longitudinally while in the structure cross-section, the temperature variation was observed. It is important to emphasize that materials properties are assigned according to the average temperature of elements, calculated from the temperature of each element node. This approach to average the nodal values is used in both thermal and thermomechanical analysis defined in the next subitem.

Thermomechanical analysis
The thermomechanical simulation is performed using coupled fields analysis offered by ANSYS 19.1 [14]. Such work consists of using thermal outputs (nodal temperatures) as inputs in the structural analysis through the command LREAD. The nodal temperatures are applied as body forces on the nodes, as showed in the flowchart presented in Figure 3. It is important to note that, for the right operation in the interest of field coupling, the mesh of the second analysis (in this case, thermomechanical) must be identical to the mesh defined in the first process (thermal). This methodology ensures the correct imposition of input data on each node of the model.
To create the thermostructural model is required that at first, a simpler analysis must be performed at room temperature. The objective is to quickly certify that the criteria established in the model design are working well. After that, it is possible to proceed safely for structural analysis in a fire situation, where the material mechanical properties, as well as the discretization of the model, must be appropriately assigned.

Thermomechanical properties
The temperature-dependent properties required are specific thermal strain, compressive strength, tensile strength, modulus of elasticity, and the stress-strain diagram of the materials. The specific thermal strains of steel and concrete are determined by Annex E of ABNT NBR 14323: 2013 [16] and Annex C of ABNT NBR 15200: 2012 [17] respectively.
For steel, the variation of the characteristic yield strength and the modulus of elasticity with the increase in temperature are obtained by decreasing the values corresponding to ambient temperature by reduction factors, , and , respectively. The reduction coefficients, distinguished for steel CA-50 and CA-60, can be found in ABNT NBR 15200: 2012 [17].
For concrete, in its turns, the mean tensile strength ( ) and the modulus of elasticity at room temperature ( ) are calculated from the mean compressive strength obtained in the laboratory ( ) according to Equations 2 and 3 suggested by the Fib Model Code [20]. The variation of the properties with the increase in temperature, as in the case of the steel properties, are obtained by decreasing the value at room temperature through the reduction factors , , , and , present in tabular form in ABNT NBR 15200:2012 [17]. According to Eurocode [18], the values of , must be obtained from , , as indicated in Equation 4. It is noteworthy that the Brazilian standards do not inform how to proceed with the , calculation, which justifies the adoption of the recommendation made by European code. Moreover, differently from national technical standardization, Eurocode [18] suggests different values for concretes consisting predominantly of siliceous aggregates from those composed of calcareous aggregates. Where: f ∆ is equal to 8 MPa; E α is a coefficient, being 1.0 for silica concrete and 0.9 calcareous. Finally, stress-strain relations to be used for both materials were defined. For the steel was attributed a perfect elastoplastic model with bilinear constitutive law, respecting Hooke's law to the point of deformation corresponding to the yield strength, followed by a straight line to the point of the ultimate tensile strength ( Figure 4). To simulate the behavior of concrete, several studies ( [21], [22], [23], among others) used the diagram proposed by Eurocode [18] (Equation 5). However, Lemos [6], opts for the model indicated by fib [20]. Moreover, as pointed out by Araújo [24], the relation of fib [20] represents in a better way the behavior of reinforced concrete. Therefore, in the current study, it was decided to analyze the structures at both room temperature and fire situation through the two mentioned diagrams. It was noteworthy that the curve suggested by fib [20] was idealized for room temperature and it was necessary to adjust the parameters as a function of temperature. The original formulation can be found in the cited reference, and the adapted formulation can be observed in Equation 6 below.   , , , Where: The schematic representation of the stress-strain relation proposed by fib [20] can be visualized in a general way in Figure 5. The strain values used for steel and concrete diagrams calculations were extracted from ABNT NBR 15200:2012 [17] and Eurocode [18].

Discretization of the structure into finite elements
To proceed with the modeling, elements compatible with those defined in item 2.1 were adopted. In this case, SOLID65 and LINK180 elements were used and there was no need to maintain the surface element for mechanical boundary conditions. The element SOLID65 used to represent concrete as well as the solid element previously described, was a three-dimensional element with eight nodes, but with three degrees of freedom per node. The LINK180 element was a two-node bar element with a degree of freedom per node capable of simulating the compressive and tensile stresses on reinforcing steel bars.

Concrete failure modes
The material model predicts concrete failure considering cracking and crushing effects [14]. These phenomena, being crushing the most important for a central compression column, imply a significant strength loss of the structural member, leading it to the failure. In the numerical model developed with ANSYS 19.1, through the SOLID65 element, according to ANSYS [14], the failure envelope proposed by Willam and Warnke [25] could be defined. The envelope aims to consider the crushing and cracking effects of concrete, performing a fundamental role in defining its behavior. The failure criterion due to a multiaxial stress state could be expressed according to Equation 7.
Where: S is a continuous failure surface.
In order to define the failure surface, five parameters are required: ultimate uniaxial compressive strength, ultimate uniaxial tensile strength ( ), ultimate biaxial compressive strength ( ) and ultimate compressive strength for a state of biaxial and uniaxial compression superimposed on hydrostatic stress state ( 1 and 2 ). The relative magnitudes of the principal stresses in the octahedral plane are described using the angle of similarity ( ) and functions 1 and 2 , as shown in Figure 6. The three-dimensional failure surface resulting from the main stress state is presented in Figure 6 and the detailed mathematical formulation can be consulted on [25]. According to ANSYS [14], it is possible to inform only two inputs (f t and f t ) for the failure envelope calculation if sufficient experimental data is not available. The other parameters are determined according to software default values: .725 c f . In the case of an element cracking, a plane of weakness was introduced in a direction normal of the crack face that meets the criterion. In order to incorporate the element stiffness reduction due to cracking, additional parameters are required. These parameters consist of shear transfer coefficients for an open ( t β ) and a closed crack ( c β ). The values range from 0 to 1, with 0 representing a smooth crack (complete loss of shear transfer) and 1 representing a rough crack (no loss of shear transfer). For the present work, as in Kumar and Kodur [23], the values of 0.53 and 0.98 for open and closed cracks were adopted.
The ANSYS documentation [14] informs the difficulty of solution convergence in cases where concrete cracking and crushing are simultaneously considered. Thus, it is possible to adopt only one failure mode, assigning the value of -1 for which it is intended to be disregarded. Based on this statement and on the experience acquired from the analyses performed during the development of this work, it was decided to consider only the concrete cracking through the Willam and Warnke model [25].
Crushing, in turn, is related to concrete strength imposed on the material stress-strain diagram. Therefore, concrete crushing was considered when the element was not able to resist the stress state originated from mechanical and thermal loading anymore.

VALIDATION OF NUMERICAL MODEL AT ROOM TEMPERATURE
In order to create a consistent numerical model, room temperature analyses were performed for reference before proceeding with coupled fields. To this end, the study elaborated by Ramos and Giongo [26] was chosen to evaluate the criteria adopted in the model developed in the current study.

Brief description of the experimental program
Ramos and Giongo [26] analyzed 16 normal strength columns considering variation in the following parameters: geometric dimensions of the structural element, longitudinal reinforcement ratio, transverse reinforcement spacing and configuration. Tests were performed on a servo-controlled hydraulic machine which allows the application of load under displacements control. The monitoring of the reinforcement deformations was conducted through electric strain gauges fixed on the steel bars and on the column faces. The deflectometers that were positioned on the column faces measured displacements derived by the compressive force. Further details on the procedures adopted in this experimental program can be obtained by consulting the above reference.

Numerical model aspects
For validation, columns models P1-10-120 and P1-12,5-100 extracted from Ramos and Giongo [26] were chosen, which will be renamed in this article as Column 1 and Column 2 respectively. The information required to perform the mechanical analysis is shown in Table 1. With the purpose of simulating the same pillar conditions found in the experimental test, the structure was modeled respecting actual dimensions. However, since the columns present double symmetry, only a quarter of their volume was modeled according to the generic model shown in Figure 7. To ensure an appropriate simulation, nodes present in the symmetrical plane had their movement restricted in the perpendicular direction.
Also, rigid body motions were restricted at nodes belonging to the structure base. For nodes belonging to the top, only translation in the load application direction is allowed. The distributed load was applied to the nodes of the upper structure section. Furthermore, extremely rigid plates were placed at the top and bottom of the model allowing uniform loading on the nodes and avoiding localized plastic deformations due to different displacement between them.
The mesh used for analysis was done manually, dividing the cross-sectional region of the covering and concrete core in a manner such that the element longest side does not exceed 1 cm. These details can be seen in Figure 7. About the solution options adopted, it is important to point out that loading was applied incrementally, dividing the total load into sub-steps. Each load sub-step starts with the last complete sub-step value and proceeds linearly until the end. This way of load application is called [14] as ramped load and must be enabled with KBC command, 0. Additionally, geometric nonlinearity was considered through NLGEOM, ON command. The convergence criteria adopted for force and displacement were 5% of the vector norm. The solution was given through the full Newton-Raphson iterative method.

Validation results
The first analysis performed with Column 1 was the behavior of different stress-strain diagrams. Figure 8 shows a good difference between the curves proposed by European Committee For Standardization [18] and Federation Internationale du Béton [20], and the latter showed a behavior closer to that found in the experimental field at room temperature. Figure 9 indicates an even greater agreement on the reinforced concrete column behavior when subjected to axial compression. As for the rupture force, Column 1 through numerical simulation presented an ultimate force of 1049 kN, representing an error of only 2.19% when compared to experimental test. The numerical model for Column 2, in its turn, presented an ultimate force equivalent to 1376.25 kN, which represents an error of 6.52%. In the face of it, the strategy adopted and described in item 2 can be regarded as satisfactory, obtaining consistent results for performed analyses and allowing to advance to the thermomechanical model. It should be noted that the processing time for Column 1 was 6 min 10 sec, while Column 2 required 6 min 43 sec.

VALIDATION OF THE THERMOMECHANICAL NUMERICAL MODEL
The numerical model was compared with results obtained by experimental tests conducted by Wu and Lie [27] and by analytical results presented in a later study published by Zhu and Lie [28] using the same columns. Thus, it is verified, with both experimental and analytical results, the proposed model validity and effectiveness in predicting the reinforced concrete columns thermomechanical behavior.

Brief description of the experimental program
Wu and Lie [27] conducted an experimental program including seven fire-exposed reinforced concrete columns, the aim of this study was to provide practical and real data to validate calculation methods that were being developed at the time. The columns were built with siliceous and calcareous aggregates, and the parameter studied were structural elements geometric dimensions and load types (centered and eccentric).
The experiments were carried out by exposing the columns to fire, fixed at the top and bottom, in a specially built furnace to evaluate preloaded columns. For centered loading, the prescribed load is applied by a hydraulic jack located under the column. After stabilization of deformation caused by loading, the structural element is exposed to fire in a controlled manner following the standard curve recommended by ISO 834 [15] which remains unchanged until the present days.
Concrete and steel temperatures were obtained using Type K Thermocouples (Cromel/Alumel) located at points of interest previously defined. The column axial displacements were measured by transducers that monitored the hydraulic jack motion during the test with an accuracy of 0.002 mm. More information about the procedures and instruments adopted in the experimental program can be obtained at Wu and Lie [27].

Brief description of the analytical method
The analytical method defined by Zhu and Lie [28] consists of different steps, involving, besides temperature calculation, calculation of mean deformations, and resistance capacity presented by columns in a fire situation.
First, temperatures were defined by subdividing the cross-sectional structure domain using the Finite Difference Method. The number of elements was set to allow that they can be arranged in a triangular mesh within the section boundaries. The methods used to derive the heat transfer equation and determine temperatures in reinforced concrete columns have been published in previous work [28].
In a second moment, the authors present the methodology used to determine strains, stresses, and strength of columns. Where the previously used triangular mesh was replaced by a square one. To obtain the column strength, the authors used a method based on a load-deflection analysis, idealizing columns, which are fixed at the ends during the tests, as pin-ended with a reduced length. From the curvature, the column strength is determined for each temperature. Strains were given for concrete and steel as the sum of mechanical strain, thermal expansion properties and elements curvature. With the values of strains, through material stress-strain relation, it was finally possible to determine the stresses acting on the structure. Further information about methods and parameters adopted for analytical calculation could be found in Zhu and Lie [28].

Numerical model aspects
In order to proceed with validation, Column 3 presented in Wu and Lie [27] and Zhu and Lie [28] was selected. The necessary information related to the column characteristics and properties can be seen in Table 1. The motivation in choosing such a column is given by analytical results demonstrated greater agreement with experimental ones.
For thermal analysis, the structure faces that will be exposed to fire and the time-temperature curve adopted are the only boundary conditions required. As discussed in section 3.2, Column 3 of [27] also allows the adoption of double symmetry. Therefore, only faces that do not represent symmetry planes are exposed to fire (see example in Figure 2). The material emissivity assigned was 0.7, and the convective heat transfer coefficient adopted was 25 W/(m 2 .k).
The standard fire indicated by ISO 834 [15] was applied incrementally, i.e., every time step of fire, the software calculates the nodal temperature evolution. The thermal analysis uses an explicit default method to solve the problem differential equations (TRNOPT, FULL). Besides, it is noteworthy that the convergence criterion adopted for temperature was 1% of the vector norm.
The load applied to the column top section before starting the temperature rise was 1180 kN, and the other thermomechanical model aspects are identical to those presented in item 3.2. The only difference is the addition of LREAD command after mechanical loading, thus obtaining the influence of temperature rise on the preloaded column (see item 2.2).

Validation results
The present subitem exposes, results obtained by the proposed thermomechanical numerical model. These results include temperature rise in concrete structure according to the depth measured from concrete face at cross-section midline -12 mm, 38 mm, 63 mm and 152 mm - (Figure 10), longitudinal bars temperature ( Figure 11) and, finally, thermomechanical analysis represented by the column axial displacement during the time of fire exposure (Figure 12). Note that, for subtitling purposes, the names EXP were used for experimental results (item 4.1), CALC for those calculated analytically (item 4.2), and ANSYS for those obtained with the software aid (item 4.3). Figures 10 and 11 show the convergence in temperature values during the standard fire. It is noteworthy that the numerical model provides, for steel, values corresponding almost exactly with those found experimentally, evidencing the satisfactory representativity obtained by computational simulation. Figure 12 presents a new comparative analysis between numerical models using stress-strain diagrams suggested by Eurocode [18] and Fib Model Code [20]. In contrast to item 3.3, when dealing with high temperatures, the relation proposed by European code presents a behavior closer to that found in experimental data.
The thermomechanical validation shown in Figure 12, although numerical simulation presents a higher axial displacement, shows especially an important concordance on general column behavior when exposed to fire. Furthermore, it indicates 103 min fire resistance, differing approximately 5.8% of the value found in the laboratory (109 min), showing that the numerical model is in favor of safety. It is important to note that the processing time for Column 3 thermal analysis was 35 min 46 sec, while for thermomechanical analysis, it was necessary 2h 47 min 20 sec.
It should be emphasized that, according to the experimental program description exposed in item 4.1, the columns are preloaded and subsequently exposed to fire. Thus, before observing a structure expansion due to the increase in temperature, there is a negative axial displacement (compression). However, this initial compressive displacement was neglected because the experimental test, as well as numerical model objective, is to analyze the structural behavior during the fire.  Additionally, Figure 12 also reveals that there is a divergence of approximately twice the column displacement when comparing the final values obtained between Eurocode relation [18] and the experimental data. This difference can be attributed to the high concrete theoretical capacity of deformation at high temperatures specified by current standards, according to the thermomechanical properties already mentioned in item 2.2.
To conclude thermomechanical result analysis, the structure specific strain (axial) is used as a parameter. This, in turn, represents the relation between the variation of axial displacement and column length. For Column 3, which length is 3.5 m, a divergence in the axial strain of only 0.0005 m/m or 0.5 mm/m is observed. This difference can be considered not representative when considering all thermomechanical variables involved and the complexity of determining reinforced concrete structures behavior when exposed to high temperatures.

CONCLUSIONS
This paper shows that the concrete constitutive law proposed for fib [20] better represented the behavior of the material at room temperature. However, it was found that increasing temperatures, simulations using the stress-strain relation specified by Eurocode [18] achieved the best results when compared to experimental data.
Moreover, it is concluded at the end of research that the proposed numerical model based on the finite element method, according to validations and comparisons made with experimental results, can predict the reinforced concrete columns thermomechanical behavior when exposed to fire.
It is also noted that the methodology presented here enables a wide range of studies in the building fire area, allowing adaptation of the created model to study other structures such as beams, slabs, and frames. The current study also allows us to expand the knowledge by parametric analyses of temperature influence on structural members with different dimensions, concrete cover, steel bar diameters, aggregate types, moisture content, among others.