Advanced computer model for analysis of reinforced concrete and composite structures at elevated temperatures

: This paper is concerned on the development of a computational model based on finite element procedures for advanced analysis capable of estimating the behavior of reinforced concrete and composite steel-concrete plane structures exposed to fire. The program implemented is called NASEN, the specific thermo-structural module is used to analyze structures under fire conditions. The effects of geometric nonlinearity, material nonlinearity and nonlinear thermal gradients are incorporated into the model, as well as the degradation of material properties with increased temperature. The methodology applied for the solution is based on the unidirectional coupling of the thermal and mechanical solutions. The cross-sections of the structural elements are discretized with two-dimensional meshes for thermal analysis, while the structural analysis uses a one-dimensional beam-column element. Numerical examples are presented to demonstrate the accuracy of the computational code developed in relation to the numerical and experimental solutions found in the literature. In summary, the program adequately predicted the response of the studied structures.


INTRODUCTION
Reinforced concrete and composite steel-concrete structures are widely used in civil construction as load-bearing mechanism in industrial and resistance buildings, bridges and among other engineering applications. In addition to conventional service loads, these structures may be exposed to fire during their design life. Therefore, the verification of these elements in response to fire is essential for designs and analysis of physical behaviors.
In the design context, the thermal and mechanical response and fire resistance of concrete and composite structural elements are usually characterized by the use of simplified calculation methods and tables based on numerical and experimental research [1]. The simplified calculation methods have a limited scope due to their ability to simulate only simple structural components, for example, Eurocode 2 part 1-2 [2] prescribes a methodology to evaluate the fire resistance of structural members using the isotherm method 500ºC.
In general, the geometric effects, support conditions and advanced physical characteristics of the systems are not considered in simplified methods. In this way, the use of advanced calculation methods based on computational tools allows a investigation on complex and highly nonlinear behavior of structures exposed to fire due to large strains and material degradation.
Currently, several engineering software packages are available in academia and industry, such as the generic solutions ANSYS [3] and ABAQUS [4], and specialist packages in the area of structures under fire conditions, such as SAFIR [5] and VULCAN [6]. In addition to finite element (FE) software packages, the development of computational codes is a methodology used in engineering. This line of research requires a knowledge of the numerical solution procedures, physical-mathematical models and characteristics of the studied physical phenomenon, enabling a wide understanding of the physical behavior of the problem [7]. Computational models able to describe the behavior of structures in fire have been developed considering different physical-mathematical approaches, such as the use of the direct stiffness method based on the stability and bowing functions [8], the approach of the embedded and extended finite element model [9], [10], the use of customized beam finite element with skeletal structures [11] and among other methodology and analysis models [12], [13], [14].
In Brazil, there is few scientific researches in the area of developing computer programs for advanced thermomechanical analysis of reinforced concrete or composite steel-concrete structures subjected to fire. In recent decades, research by Caldas [15], Mouço [16], Ribeiro [17] and Maximiano [18] stands out. Recently, Neves [19] presents the computational model NASEN (Numerical Analysis System for Engineering), a program aimed at solving initial engineering problems in the acoustic, thermal, structural and thermo-structural areas. This study presents the specific module for thermo-structural analysis of concrete or composite plane elements under fire condition, NASEN/TSA-FIRE (Thermal-Structural Analysis-Fire). Among the general characteristics of the code, an updated Lagrangian formulation of a beam-column plane element associated with the fiber model is considered, the effects of thermal gradients on the cross-section, degradation of physical and mechanical properties as a function of temperature, large displacements and the actions of the equivalent thermal forces.

BASIC THEORY OF THERMO-STRUCTURAL ANALYSIS
The numerical solution of the thermomechanical problem of structures under high temperatures is based on the knowledge of aspects and parameters related to the modeling of the gases of the environment in fire, the knowledge of the temperature distribution in the cross section of the elements and the behavior of the mechanical response of the structural system.

Heat transfer analysis
The prediction of the thermal field is computed at the level of the cross-section of the structures under fire condition, using a discretization based on plane elements, as shown in Figure 1. The governing equation of the physical phenomenon is described by the transient differential model of conduction of two-dimensional heat [1], [20], as given by Equation 1.
Where θ is the temperature of the structure, t is the time, f is the internal heat generation, ρ is the density of the material, c is the specific heat and = λ D I is the thermal conductivity matrix for an isotropic material. In the area of structures in fire, the physical properties of the materials vary with temperature, making Equation 1 nonlinear. The properties of concrete and steel follow the recommendations prescribed in Eurocode 2 part 1-2 [2] and Eurocode 3 part 1-2 [21]. The effects of convection and radiation caused during exposure to fire are considered by the boundary conditions of the problem, according to Equation 2.
Since g θ is the gas temperature in the proximity of the fire exposed member, c h and r h are called convection and radiation coefficients respectively. Temperature-time curves are used to simplify the evolution of fire behavior [22], for example, the ISO 834 [23] and ASTM E119 [24] curves, as shown in Equation 3 and Where 0 θ is the initial ambient temperature, t and h t are the times of exposure to fire in min and hours respectively. In this way, the numerical solution of Equation 1 is based on Galerkin finite element procedures. The approximation of the temperature field is constructed based on the interpolation functions [7], [25], as shown in Equation 5. . Thus, applying weighted residual sentences, theorems of differential-integral calculus and the approximation given in Equation 5, the resulting system of equations in matrix form is given by Equation 6.
The total capacitance matrix θ K is composed of the contribution of the thermal conductivity matrix and the matrix due to the combined effects of convection-radiation, as given by Equation 7.
The temperature gradient matrix B is written according to the finite element interpolation functions, ∇θ = ∇ =

Nθ Bθ
. The thermal capacity matrix θ C is given by Equation 8, considering the interpolation functions and the product between density and specific heat.
The vector of thermal actions θ F is given by Equation 9.
The finite difference technique is applied to approximate the first order temporal differential operator [19]. The effective algebraic system computed at n 1 + is given by Equation 10.
Where ∆t is the time step and β is the time integration parameter equal to 2/3, characterizing an unconditionally stable scheme.

Mechanical analysis
The structural nonlinear analysis is constructed based on the plane beam-column element with six degrees of freedom (6DF), two translations and one rotation at each element node, as shown in Figure 1. Thus, the differential deformation relationships-displacement based on the Green-Lagrange tensor, accounting for axial and shear deformations, as shown below.
The first terms represent the linear parts and the other terms correspond to the non-linear parts, according to Equation 11 and Equation 12. The Euler-Bernoulli hypothesis is also adopted, which establishes that the plane cross-sections remain normal to the axis of element in flexion. In this way, the displacements x u and y u at a generic point can be associated with the displacements u and v of the beam, as given by Equation 13.
Thus, following the principle of virtual displacements and performing some mathematical operations, the integral sentence resulting from the nonlinear structural problem, disregarding the high order terms [26] is given below by Equation 14.
Based on finite element procedures, the field of axial and vertical displacement is approached with the aid of interpolation functions, as shown in Equation 15 and Equation 16 respectively.
The terms 1 N and 3 N represent the vectors of the interpolation functions, linear for axial displacement and cubic for transversal displacement respectively. The approximations presented in Equation 15 and Equation 16 are applied in the integral formulation of the problem, as shown in Equation 14. As the virtual displacements are arbitrary at the natural boundary, consequently, the resulting algebraic system is given by Equation 17.
Where ∆d and ∆f are the displacement and force incremental vectors, respectively. The matrix of elastic stiffness e K is obtained by the numerical evaluation of the first terms of Equation 14, as shown below in Equation 18.
The nonlinear term is computed in the geometric matrix g K , given by Equation 19, where a simplified theory is adopted that disregards the combined effects of axial and flexural behavior [27].
Having defined the elementary elastic and geometric matrices for each element of the structural system, the global stiffness matrix is constructed by the composition of the local matrices of each element. In addition, the equivalent axial (EA) and flexural (EI) stiffnesses are calculated in relation to the fibers of the cross-section, as shown in the following.
In the computational model developed, the reduction factor of the elasticity modulus of the concrete as a function of temperature, based on the experimental results obtained by Hertz [28] at the Technical University of Denmark, is approximately equal to the square of the reduction factor of the compressive strength of the concrete. This simplification is a particular characteristic of the zone method described in Eurocode 2 part 1-2 [2].
The equivalent stiffnesses are computed based on the values of the elasticity modulus of the materials. In a fire condition, there are numerous models available in the literature to describe the degradation of the modulus of elasticity of steel [21], [29]- [31] and concrete [2], [32]- [36] as a function of temperature rise, as can be seen in Figure 2. The effects of temperature in the cross-section of the structure are computed in the beam-column model based on the vector thermal fixed-end forces th f , consisting of the contributions of the effects of thermal expansion and curvature due to the gradients temperature [16]. The components of the vector are shown in the following.
Where th ε is coefficient of thermal expansion of the materials due to the increase in temperature in the structural element, following the behavior described in European codes [2], [21]. In addition, the vector thermal fixed-end forces is included in the global force vector of the structural system.

Unidirectional coupling of thermal and mechanical processes
The two-way unidirectional (or sequential) coupling is based on thermal and mechanical solvers in a similar strategy used in the research by Barros et al. [37], Caldas et al. [38] and Silva and Landesmann [39]. Figure 3 illustrates the general process of the solution steps adopted in the advanced computational model developed. The program starts by performing a structural nonlinear analysis at ambient temperature, considering external loads. In the occurrence of fire, the gas temperature is determined by following the standard temperature-time curves or experimental data from test furnaces. Then, the temperature distribution in the structural cross-section is computed and the material properties and equivalent thermal forces are calculated due to fire exposure for each element. In the last solution stage, structural analysis is performed considering the effects of the thermal gradients of the cross-section, obtaining the displacement field. The solution of each system of equations is based on direct methods and iterative processes are implemented due to the nonlinear behavior of the analyzed problems. At the end of the computer simulation an investigation of the results obtained is carried out.

NUMERICAL APPLICATION TESTS
The performance evaluation of the developed program is based on the analysis of reinforced concrete and composite steel-concrete structures. The thermophysical properties of the materials and the parameters used in the simulations follow the recommendations of European codes [2], [21], [22]. For concrete, the variation of density with temperature is defined by Equation 22 and the value of 2300 kg/m 3 at 20ºC.
The specific heat of dry concrete as a function of temperature is determined by Equation 23. In addition, the specific heat of the concrete is influenced by the moisture content of the concrete. According to the recommendations of Eurocode 2 part 1-2 [2], this behavior is modeled by a constant peak value between 100ºC and 115ºC, and followed by a linear decrease between 115ºC and 200ºC. The peak value, for 3% moisture content by weight, is equal to 1470 J/kg K.
The upper limit of the thermal conductivity of normal weight concrete varying with increasing temperature is described by a quadratic expression, as given by Equation 24.
For steel elements, the specific mass is constant and equal to 7850 kg/m 3 . The specific heat of steel as a function of temperature is determined by Equation 25.
The variation of the thermal conductivity of the steel with the temperature is modeled by a simple mathematical expression, as given by Equation 26.
In addition, the parameters related to the convection and radiation heat transfer mechanisms, for the face exposed to fire, a convection coefficient equal to 25 W/m 2 ºC is used, while for the face unexposed to fire the value of the coefficient is taken as 9 W/m 2 ºC. The emissivity factor, related to the effect of heat transfer by radiation, is equal to 0,7.

Reinforced concrete beam subjected to fire
A numerical analysis of a simply supported reinforced concrete beam with an overhanging subjected to fire is performed in this section. The structure was subjected to a fire modeled with the ASTM E119 curve [24]. This example is part of the experimental tests carried out at the Construction Technology Laboratories of the Portland Cement Association by Ellingwood and Lin [40]. The beam has an overhang of 1,83 m and length between supports of 6,1 m. The fire acts only on the central span, while the overhang remains at ambient temperature. Figure 4 shows the geometric configuration of the problem. MPa, respectively. This example has been numerically analyzed by several authors [10], [15], [18], [41].
The first step of the analysis seeks to validate the temperature increase in the reinforcement steel bars of the section. Results are shown in Figure 5. It is noted that the temperature change in the top and bottom bars of the reinforced concrete section presents acceptable results in relation to reference solutions. Top bars present slightly lower temperature values than bottom ones, due to the upper face of the concrete section not being directly exposed to fire. Although temperature-time curves are commonly used to evaluate the performance of a numerical model, the thermal response can also be visualized by a two-dimensional temperature field in the cross section of the structure, allowing identification of critical regions and a qualitative estimate of temperature change within the domain. As such, Figure 4 shows the temperature distribution in the reinforced concrete cross-section at 30, 90 and 180 min of fire exposure.
Finally, the evolution of the vertical displacement at mid-span between supports is measured as a function of time of heat exposure. Figure 6 shows a result comparison between NASEN and reference solutions from the literature. The results obtained with NASEN present reasonably good values, and the best agreement with the experimental tests occur between 45 and 200 min. Outside of this time range, the values diverge slightly. Furthermore, the solution extracted from NASEN is similar to the results obtained by Caldas [15].

Reinforced concrete column with eccentric load
At the Technical University of Braunschweig, Hass [42] carried out a series of full-scale experimental tests on reinforced concrete columns exposed to fire. Among the tests performed, three of these cases are studied numerically in Bamonte and Lo Monte [43]. For analysis, the case known as Hass 16 is considered, characterized by a square column of reinforced concrete heated according to the ISO 834 curve [23], length equal to 4,76 m, axial load of 460 kN, eccentricity equal to 9 cm and the dimensions and details of the cross section are shown in Figure 7. Regarding the physical characteristics of the problem, the modulus of elasticity and the yield strength of steel equal to 210 GPa and 462 MPa respectively, while the characteristic value of compressive strength of concrete is equal to 30,7 MPa. The concrete column is discretized into 7 one-dimensional beam-column elements and 618 triangular elements are used to mesh the cross section. The results are based on the analysis of the evolution of the axial displacement measured in the middle of the concrete column, as shown in Figure 8. The numerical solutions obtained in Maximiano [18] and Bamonte and Lo Monte [43] are used to calibrate the performance of the program NASEN. The results are shown to be in good agreement with the previsions in the literature. In addition, usually in structural design, simplified calculation methods are applied, for example, the 500ºC isotherm method for concrete structures under fire condition [2]. Using the characteristics of the studied physical problem, the NASEN program can be applied to compute the effective section representative of a concrete region of the inner crosssection the isotherm of 500ºC for 30, 60, 90 and 120 min of exposure to fire, assuming that the concrete with temperature higher than 500ºC is completely neglected, as shown in Figure 9. In a simplified way, the concrete contained in the interior region is not affected by fire, considering the properties of materials at ambient temperature, however, exceptional coefficients are adopted in the elaboration of the structural design [44]. Furthermore, another important analysis in the context of fire design of structures is the investigation of the influence of the concrete column cover. Table 1 shows the temperature values measured in the upper steel reinforcing located in the region that is most heated by fire and unfavorable for safety. In order to make comparisons and discussions of the results obtained with the computer program, the simplified value of 550ºC is adopted as the reference temperature of the steel reinforcement of the concrete cross section. This value is not prescribed in standards and is used only as a reference value for the present study, since for temperatures above this value, steel shows significant reductions in stiffness and strength. As can be seen in Table 1, the steel reinforcements of the concrete column are exposed to high temperatures when the cover values are lower, which can provide severe risks to the safety of the structure. In contrast, higher values of the concrete column cover, the temperature levels decrease. However, it is important to check the accuracy of the details of the steel reinforcement and the requirements in the preparation of the structural design, in order to guarantee the requirements of resistance and fire safety.

Steel-concrete composite beam
The thermomechanical behavior of simply supported steel-concrete beams subjected to concentrated loads in a fire situation is studied using an advanced computational model. The system is exposed to the standard temperature-time curve ISO 834 [23]. Two configurations are analyzed, the first case is shown in Figure 10a, where the steel profile is partially encased by the concrete and subject to a symmetrical load equal to 36 kN.  Dotreppe et al. [45] presented numerical results and experimental tests for first case, carried out at the Technical University of Braunschweig. Results from both approaches are adopted as references to evaluate the numerical results obtained with NASEN. A mesh with 6 one-dimensional elements is adopted to discretize the structural model. The validation of results for temperature evolution at the cross section and the deflection of the beam as a function of the time of exposure to fire is shown in Figure 11. It is observed that the results obtained with NASEN once again present satisfactory agreement with the experimental data. Figure 12 shows the two-dimensional mesh with 477 nodes and 879 three-node linear triangular elements used for discretization of the cross-section, as well as the temperature distribution at 105 min of fire exposure. The second case analyzed is shown in Figure 10b, where the steel profile is only in contact with concrete slab. The beam is subjected to four concentrated loads equal to P with an intensity of 62,36 kN. The support conditions adopted in the structural model allow rotation on the left end and horizontal and rotational displacement are free on the right support. The proposed example is reported, initially, to the tests performed by Wainman and Kirby [46], later, studied numerically by Huang et al. [47]. The numerical meshes adopted for the discretization of the cross section and the structural model are, respectively, equal to 413 linear triangular elements and 6 one-dimensional finite elements. Figure 13 shows the temperature profile along the center line of symmetry of the cross section. It is observed that the temperature levels in the steel profile remain practically the same, while in the concrete the thermal gradients of temperature are noticed.  Figure 14 shows the maximum vertical deflection in the beam in relation to the temperature in the bottom flange of the steel profile. The results obtained with the NASEN program indicate a good agreement between the experimental tests, as well as in relation to the numerical model in the literature. It should be noted that the numerical solution with the developed program shows a more conservative behavior over the temperature increase in the structure compared to the experimental results.

CONCLUSIONS
This work presents the nonlinear numerical procedures of finite elements used in the development of the NASEN computational tool for thermomechanical analysis of reinforced concrete and composite steel-concrete plane structures at high temperatures. Different loading configurations, cross-section characteristics and boundary conditions are studied.
The results obtained with the developed program were calibrated with data from computer simulations and experimental tests available in the literature. The predictions of thermal behavior were well evaluated by the developed program, since the solution of scalar problems are already well known in the literature. The presence of non-linearity of the material requires some additional numerical procedures, which does not present difficulties in the implementation of the code. The mechanical responses present results consistent with the predictions in the literature. However, these structural models are more sensitive and complex to simulate in relation to the thermal model. It can be observed in each test performed, the NASEN program and the models in the literature present different behaviors, due to the physical hypotheses adopted in the development of each model.
In general, the results show satisfactory behavior in all cases studied, indicating the good performance of the NASEN program. Hence, this advanced numerical model allows to investigate and simulate adequately the behavior of structures in fire situations. In future projects and research, updates are made to the models and implementation of new computational modules in the program.