Numerical analysis of mechanical damage on concrete under high temperatures

: Concrete is a widespread material all over the world. Due to this material’s heterogeneity and structural complexity, predicting the behavior of concrete structures under extreme environmental conditions is a very challenging task. High temperatures lead to microstructural changes which affect the macrostructural performance. In this context, computational tools that allow the simulation of structures may assist the analysis, by reproducing varied situations of thermal and mechanical loading and boundary conditions. In order to contribute to this scenario, this study proposes a numerical methodology to simulate the thermomechanical behavior of concrete under temperature gradients, through inverse analyses and a user subroutine implemented in Abaqus software. Thermal loading effects were considered as loading data for a damage model. Experimental data available in the literature was adopted for adjustment and validation purposes. The preliminary results presented herein encourage further improvements so as to allow realistic simulations of such an important aspect of concrete’s behavior.


INTRODUCTION
According to ASTM (American Society of Tests and Materials), it is possible to define concrete as a composite material, which comprises a binding agent medium in which are incorporated different aggregates. According to Mehta and Monteiro [1], concrete is one of the most widely adopted construction materials, with world consumption of the order of 33 billion tons per year in 2016. In spite of this fact, the prediction of concrete behavior is rather complex, especially under extreme situations such as large displacements and high temperatures.
Exposure to high temperatures may be due to accidental reasons, such as fire situations, or to ordinary service conditions, which is the case of some components of nuclear power plants, blast furnaces and radioactive waste repositories. In either case, the knowledge of the effects of the temperature elevation on concrete's properties is paramount to an adequate executive or corrective design.
The evaluation of the behavior of concrete when subjected to high temperatures has increasingly aroused the interest of the scientific community. In this context, several constitutive models based on the continuum mechanics have been developed with this objective in mind. Mazars [2] proposed an isotropic damage model based on a single scalar variable in which the concrete has a damaged elastic behavior. In this model, the damage is directly associated with the existence of positive deformations. Thelandersson [3] described the thermomechanical response of concrete considering the thermal strain rate as a function of both rate of temperature change and the current state of stress. Damage is accounted by the change in elastic properties due to a temperature change. Simo and Ju [4] developed a cap plasticity model with an isotropic strain-based damage mechanism and proposed a viscoplastic extension for their model. Cervera et al. [5] implemented a rate-dependent isotropic damage model that incorporates stiffness degradation and stiffness recovery upon load reversals and strain-rate sensitivity in order to make a realistic seismic analysis. Pituba and Lacerda [6], on the other hand, admitted concrete as an initially isotropic medium that starts to present plastic deformations, bimodularity and damage-induced anisotropy.
In the field of experimental studies, Lima et al. [7] analyzed the damage caused by the increase in temperature in a reinforced concrete building. Arioz [8] studied the effects of high temperatures on the physical and mechanical properties of various concrete mixtures, determining weight losses and compressive strength after exposure. Arioz [8] concluded that weight losses and compressive strength are directly linked and that both decrease considerably with increasing temperature. Furthermore, in this study, the damage can be observed macroscopically on the concrete surface. Ehrenbring et al. [9] realized the inspection of a prefabricated hollow core slab of an industrial building, which was exposed to high temperatures due a fire in the underground of the same, estimating the loss of strength of the structural element and attesting the safety of the structure after the accident. Bailey and Toh [10] tested forty-eight horizontally unrestrained two-way spanning reinforced concrete slabs at ambient and elevated temperatures comparing the modes of failure observed in each case in order to provide a wealth of data, which can be used to further develop simple design methods. Souza et al. [11] presented the use of the factorial experimental design methodology with a star configuration for the study of the compressive strength of specimens heated in kilns. Dias et al. [12] studied the concrete subjected to high temperatures and concluded that in this situation the material undergoes significant deterioration as a reduction in the modulus of elasticity and resistance to compression, surface displacement and loss of durability.
As for the work involving numerical analysis, Ribeiro [13] developed a computer system for simulating the behavior of structural elements in a fire situation based on the Finite Element Method and obtained results similar to those observed in experimental tests. Subsequently, Padre et al. [14] used the created program, and developed an algorithm to check the resistance of reinforced concrete sections to oblique composite bending at ambient temperature and in fire situation. Izzuddin and Elghazouli [15] presented two analytical models, detailed and simplified, for the nonlinear analysis of axially restrained lightly reinforced concrete members under ambient and fire conditions with emphasis on the evaluation of steel reinforcement failures. Grassl and Pearce [16] used a meso-scale approach through a damageplasticity model -considering concrete as a three phase material composed of aggregates, matrix and interfacial transition zones -to evaluate the transient thermal creep, concluding that this phenomenon results from the mismatch of thermal expansions of the meso-scale constituents. Ferreira [17] developed a numerical model for the simulation of reinforced concrete columns in a fire situation, also based on the Finite Element Method, in order to predict the structural behavior under high temperatures. Nguyen et al. [18] studied the thermal conductivity change of concrete when exposed to both mechanical and thermal loads through a numerical three-phase plane model using lattice discretization, where damage variable is accounted via crack opening. Rodovalho et al. [19] analyzed the mechanical resistance of concrete blocks subjected to compression and in a fire situation. Assis and Neto [20] evaluated the heat transfer in the cavity of the structural masonry blocks. Barbosa and Haach [21] analyzed the influence of fire in a room of a structural masonry building and checked the large reserve of resistant capacity of this type of construction, in addition to the great ability to redistribute efforts on the walls structural. Machado et al. [22] evaluated the stability of reinforced concrete panels under fire conditions, verifying the influence of the time of exposure to fire on the results obtained.
In turn, this work presents a numerical methodology consisting of simulating the material deterioration through the implementation of Mazars' [2] damage model in a finite element software. Experimental data available in the literature was adopted for data adjustment and numerical validation, and the results encourage further improvements so as to allow the reproduction of varied geometrical, loading and boundary conditions.

Overview
This work employed a set of experimental data provided by Razafinjato [23] who studied the thermomechanical behavior of concretes of varied compositions.
According to Razafinjato [23], the experimental data were obtained from cylindrical specimens with a proportion of 40% of aggregates. These samples were molded covered with a damp cloth and stored in plastic bags for 90 days. Then the specimens were subjected to heating/cooling cycles, in a furnace, at a rate of 0.50K/min till temperatures of 573.15K, 723.15K, 873.15K and 1023.15K, maintained for two hours. After the heating/cooling cycles, at room temperature (293.15K), the specimens were taken to the press so that, from the uniaxial compression test, the residual concrete Young's modulus could be determined. Figure 1 shows the resulting Young's modulus and temperature relation, for 293.15K, 573.15K, 723.15K and 873.15K -in which it is possible to observe the linear trend of the experimental results adjusted via the least-squares method. For the temperature of 1023.15K it was not possible to perform the uniaxial compression test due to the deterioration of the specimen. These points were used to calculate the objective function of the inverse iterative procedures that calculate mechanical properties and damage parameters, which will be presented in detail in the following sections.

Numerical procedure
To evaluate the damage process in the concrete subjected to high temperatures the first step was to get, through the algorithm developed by Bonifácio [24], the geometric model of the cylindrical specimen of concrete. Then, it was developed a linear elastic analysis employing the Abaqus [25] finite element commercial program, to determine the unknown properties -Young's modulus of aggregate and Poisson's ratio of aggregate and mortar -for initial temperature, without damage, through a numerical adjustment. This process can be seen schematically in Figure 2.
Once evaluated the necessary property, the weakly coupled thermomechanical model was implemented in two steps. Firstly, the thermal model receives the values calculated in Flowchart 1 ( Figure 2) and estimates the specimen temperature field. Then, the temperature field is applied as thermal loading and a damage model is adopted so as to evaluate the material's deterioration. To evaluate the concrete deterioration process, opted to use the Mazars' [2] damage model, which is unavailable in the library of the Abaqus [25]. So, a UMAT subroutine [26] was implemented for the software to perform the analysis considering the desired model. However, Mazars' [2] parameters for the material were unknown and, therefore, it was necessary to use another numerical adjustment.
Finally, from the results of the damage model, it was possible to construct a damage curve as a function of temperature, from which we carried out the analysis of the behavior of the concrete under high thermal gradients.
We emphasize that the strategy used, which is summarized in a weak coupling, has the advantage of low computational cost. This process can be seen schematically in Figure 3. Both in the inverse problems employed to determine the unknown mechanical properties and in the necessary to find the parameters of the Mazars' [2] model, it was used experimental data from a thesis developed at the Cergy-Pontoise University [23], in France, and also functions of the package optimize from the SciPy library [27].

Governing equations
Regarding thermal models, in this study, the transfer of heat by conduction was considered, which is described by Equation 1. From the resolution of a thermal problem with this governing equation, it was possible to find the temperature field in the specimen section considered.
where ρ is the density, c is the specific heat, T is the temperature, t is the time, κ is the thermal conductivity.
In relation to linear elastic model and to damage model Equation 2 was adopted, being that in the second case the Mazars' [2] damage model was considered.
where B are the body forces and σ is the normal stress. About the Mazars' [2] damage model, it modifies the modulus of elasticity E of the intact material according the damage begins to develop, as can be seen in Equation 3, where d E is the Young's modulus of the damaged material and d is the damage variable, which is given by Equation

Geometric model and mesh
So as to simulate the laboratory experiments developed by Razafinjato [23], it was necessary to reproduce the specimen geometry, by taking into account the adopted aggregate gradation and proportion. For reasons of symmetry, it was decided to consider 1/4 of the longitudinal section of a cylindrical specimen with 150mm in diameter and 300mm in height, composed of two phases: mortar and aggregate.
The center coordinate and the aggregates radius were generated using an algorithm developed by Bonifácio [24], in Python language, able to simulate the distribution of aggregate in a specimen, having as input data its dimensions, aggregates proportion and grading.
The model was discretized in 12175 linear elements, after convergence test. For the elastic model (Flowchart 1 - Figure 2) and in damage model (Flowchart 2 - Figure 3), elements of CPS4 type (4-node bilinear plane stress quadrilateral) and CPS3 type (3-node linear plane stress triangle) were employed. For the thermal model (Flowchart 2 - Figure 3), elements of the DC2D4 type (4-node linear heat transfer quadrilateral) and DC2D3 type (3-node linear heat transfer triangle) were adopted [25]. The adopted mesh is presented in Figure 4. . Mesh adopted to discretize the model.

Linear elastic model mesh
From mechanical properties, by Razafinjato [23], we had the Young's modulus of mortar and concrete. However, the Young's modulus of aggregate as well as the Poisson's ratio of the aggregate and the mortar were unknown and, to determine them, an inverse procedure was used.
The numerical procedure consisted of simulating a 1/4 of the longitudinal section of a cylindrical specimen with 150mm in diameter and 300mm in height, composed by two phases, mortar and aggregate, under increasing axial loading, through a displacement applied at the top of the sample as represented in Figure 5. Flowchart 1, in Figure 2 represents the applied inverse methodology, based on a trial and error procedure implemented in Python language, applying the differential_evolution function, from the optimize package of the SciPy library [27], assuming as objective function the error between the experimental Young's modulus and numerical. The latter was considered as the ratio between the vertical reactions in response to the applied displacement. Table 1 shows the obtained values.

Thermal model
A transient thermal analysis was performed via Abaqus [25], which adopts the Finite Element Method so as to evaluate the temperature field resulting from exposing a concrete specimen to thermal loading. The loading, initial and boundary conditions are indicated in Figure 6. It is emphasized that the initial and boundary conditions were chosen so to represent the experimental process developed by Razafinjato [23].  Figure 7 shows the data provided by Razafinjato [23] for the thermal expansion coefficient of aggregate and mortar. The dashed lines are the simplified curves adopted herein. The curve is limited to 803.15K, since higher temperatures lead to strong non-linearities. Due to the lack of experimental data regarding the specific heat of the mortar, it was considered as an extrapolation of the concrete curve [23], as shown in Figure 8, and therefore it was assumed that the variation of this property in the mortar is proportional to the variation in the concrete. As a simplification, the specific heat of the granite was considered constant. The initial value for the specific heat of the mortar, the constant value for the specific heat of the granite and the thermal conductivity for both phases were provided by the Brazilian Standard NBR 15220-2: 2005 [28]. The density was adopted according to Razafinjato [23]. Table 2 presents a summary of the assumed properties of the materials.

Damage model
Concrete deterioration was then quantified through the Mazars' [2] damage model, which is not available in the Abaqus [25] libraries. For this reason, the model was incorporated to the analysis via a user UMAT subroutine [26], schematically described in the Flowchart 3 shown in Figure 9. We emphasize that, in this work, it was considered that only the mortar and, consequently, the concrete, suffer damage. For the aggregate, the linear elastic behavior was adopted, without damage. Furthermore, in the proposed model, temperature variation is of the increasing monotonic type, that is, no cooling is considered. Thus, the level of compression identified in the concrete was low and, consequently, the estimated c α values were lower than the t α values. Thus, in Equation 5 it was adopted c α = 0, which reduced the parameters of interest to t A , t B and 0 d e , in Equation 4.
These parameters were obtained from a second inverse problem using the function f_min of the package optimize from the SciPy library [27] in Python language, having as objective function the error between experimental and numerical mortar Young's modulus. This process can be seen in Flowchart 2 ( Figure 3) and the obtained Mazars' [2] parameters may be seen in Table 3.

RESULTS AND DISCUSSIONS
A thermal model was developed on the Abaqus [25] to evaluate the thermomechanical behavior of concrete through the Mazars' [2] damage model, obtaining the temperature field active in the entire section of the specimen. Figure 10 shows the uniform distribution of the initial temperature T 0 = 293.15K and the temperature range obtained for T = 573.15K, T = 723.15K and T = 803.15K. The mesh used in the problem was shown at the initial temperature so that it was possible to identify the problem geometry and hidden in the other temperature levels for better visualization of the results. The maximum temperature obtained in the section is found at the top and on the right side, which are the faces in contact with the thermal flow. Consequently, in the center of the specimen, the temperature was the minimum of the section. A damage model was implemented on the Abaqus [25] using the temperature field as imposed loading. As already explained, at this moment, a UMAT subroutine [26] was used for the model Mazars' [2] to be considered and to find its parameters an inverse problem was solved, from which the values of interest have been determined. Figures 11 and 12 show, respectively, the damage and Young's modulus map for temperatures T 0 = 293.15K, T = 573.15K, T = 723.15K and T = 803.15K. The mesh used in the problem was shown at the initial temperature so that it was possible to identify the problem geometry and hidden in the other temperature levels for better visualization of the results. It is possible to identify that the greatest damages, and consequently, the smallest Young's modulus, are present in the areas of mortar that interconnect the aggregates.   Figure 13 shows the mortar and concrete damage evolution. It is possible to observe that the beginning of the damage in the section occurs around 350K, reaching expressive values from 380K. In addition, it is observed that the final damage of the mortar is 0.84. As for concrete, there is a slightly lower value of 0.74, justified by the fact that the aggregate has null damage and, thus, contributes to the decrease in the homogenized Young's modulus being smaller. Therefore, at the final temperature of 803.15K Young's modulus of the aggregate, mortar, and concrete are equal to 39486MPa, 4920MPa, and 9000MPa respectively. It is important to notice that since the aggregate does not present damage, its Young's modulus keeps constant, as shown in Table 1. In turn, the Young's modulus of mortar is estimated by considering the weighted average of this property in each element of the mesh belonging to this phase. Finally, the Young's modulus of concrete is obtained through the simulation of the uniaxial compression test assuming elastic behavior, as indicated in Figure 5, and by considering the previously estimated values for the mortar and the aggregate. Figure 14 shows the evolution of Young's modulus obtained experimentally [23] and numerically -normalized from 0 to 1 -and the concrete damage variation with temperature. It is observed that the damage is inversely proportional to the Young's modulus. It is also verified that Young's modulus is reduced by half at 650K, the temperature under which the damage has a value of 0.5. In addition, the relative error between the numerical and experimental values of concrete Young's modulus under 293.15K, 573.15K, 723.15K and 803.15 K was approximately 0.03%. Then, considering the Mazars parameters and the mechanical properties obtained previously, the damage model was applied to other specimens in order to analyze the influence of aggregate grading on the damage evolution. In this sense, specimens were created with 40% of aggregates of fixed dimensions, of 4mm, 8mm, 16mm and 32mm to assess which diameter would result in the greatest damage. As a result, the graph in Figure 15 was obtained, showing that the larger the diameter of the aggregate, the greater the damage, although the difference is not so significant. Specimens were also created with different proportions of aggregate, 40%, 30%, 20% and 10%, keeping the aggregate grading curve constant. As a result, the graph in Figure 16 was obtained, showing that the lower the percentage of aggregate, slower is the initial increase of the damage. Regarding the damage to higher temperatures, there are three simultaneous phenomena. The first concerns the damage to the mortar, which is directly proportional to the amount of aggregate. The second is about reducing the initial Young's modulus of concrete by reducing the relative volume of aggregate that has a Young's modulus greater than that of the mortar. The third, on the other hand, is related to the increase in the damaged portion of the model, by reducing the percentage of aggregates that represents the nondamaged portion of the model. These three phenomena, with opposite effects, mean that there is no well-defined pattern for the entire temperature domain, although at the final temperature of 803.15K the highest percentage of aggregates results in slightly higher damage.

CONCLUSIONS
The present work achieved the proposed objectives: the Mazars' [2] damage model was implemented in the Abaqus [25] software through the UMAT user subroutine, making it possible to observe how the increase in temperature influences the concrete damage and, consequently, in the material Young's modulus, obtaining results consistent with the theory.
Also, the proposed computational model was able to represent the results obtained experimentally by Razafinjato [23] with an error of the order of 0.03%, which is relatively small. Besides, the data available were referents to only four temperatures, while by the curve obtained numerically it was possible to evaluate Young's modulus at any point, between the maximum (803.15K) and minimum (293.15K) temperatures considered.
It is noteworthy the performance of the proposed model, which, with little experimental data, is able to contribute to the elucidation of the mechanism of damage generated in concrete subjected to high temperatures.
In this study, the Mazars' [2] model was implemented for a bidimensional domain. Thus, the implementation of a subroutine in which the model is considered in its tridimensional geometry, would be an improvement on the work done so far.