Numerical evaluation of aggregate size influence on concrete mechanical damage under high temperatures

Abstract The influence of aggregate size on the degradation process of material exposed to high temperatures is not a consensus among the scientific community because changes in the microstructure impact the macrostructural performance. To contribute to this investigation this work presents a thermomechanical model to evaluate aggregate size influence on the concrete mechanical damage under high temperatures. The material is considered as two-phase - aggregate and matrix - and three-phase - in which the interfacial transition zone is added. Concerning geometries, models in 2D and 3D are simulated. A finite element software is used with a weak coupling strategy that reduces the computational cost, and a user subroutine is implemented to define the constitutive model. The results show that the aggregate size influences both the average damage and the damage distribution along the synthetic specimen.


INTRODUCTION
Concrete, a heterogeneous mixture basically composed of cement, coarse aggregate, fine aggregate and water, is amongst the most widely adopted construction materials in the world [1]- [3], with an annual consumption of about 4,7 tons per inhabitant per year [4].Its intrinsic complexity, especially regarding microstructure changes under extreme situations, such as high temperature, justifies the need for investigation by the scientific community.
Exposure of the structure to high temperatures may occur accidentally, in fire episodes, or as ordinary service conditions, as in nuclear powers, blast furnaces and radioactive waste repositories.In both situations, it is very important to have information on the material's response to high temperatures, so as to allow an adequate executive or corrective design.In this context, several computational and experimental studies about the subject have been carried out.
Amongst the available experimental studies, some are interested in post-fire mechanical properties.Arioz [5] studied the effects of high temperatures on the physical and mechanical properties in various concrete mixtures, determining weight loss and compressive strength after exposure, and concluded that both are directly linked and decrease significantly with temperature increase.Morales et al. [6] studied cylindrical mortar specimens submitted to temperature rise, determined their residual strength after cooling and concluded that in all analyzed temperatures there was a strength decrease.Teixeira [7] evaluated the influence of initial strength, temperature level and unheating rate on concrete Young's modulus after a fire and concluded that the lower the initial strength, the greater the decrease in stiffness and that abrupt unheating has a stronger effect on the mechanical residual properties.
Dealing with structures, Bailey and Toh [8] tested forty-eight horizontally unrestrained two-way spanning reinforced concrete slabs at ambient and elevated temperatures, comparing the observed failure modes to provide data to be used in the development of simple design methods.Ehrenbring et al. [9] performed an inspection of a prefabricated hollow core slab of an industrial building, which was exposed to high temperatures due to a fire, estimating the strength loss of structural element and attesting structure safety after the accident.
Regarding computational studies, there are some Ph.D. and M.S thesis talking about it.Ribeiro [10] developed a computer system for simulating structural elements behavior in a fire situation, and had good agreement with experimental tests.Santiago Filho [11] carried out an analysis of the effects of temperature rise in reinforced concrete slabs, and the results were validated through experimental data available in the literature.Ferreira [12] developed a model for simulation of reinforced concrete columns under a fire situation, based on the Finite Element Method, to predict structural behavior under high temperatures.
There are also scientific articles interested in this issue.Using the Finite Element Method (FEM) framework, Grassl and Pearce [13] adopted a mesoscale approach through a damage-plasticity model by considering concrete as a three-phase material composed of aggregate, matrix and interfacial transition zone (ITZ) so as to evaluate transient thermal creep, concluding that this phenomenon results from a mismatch of thermal expansion of mesoscale constituents.Mazzucco et al. [14] evaluated numerically the complex mechanism of polypropylene contribution on concrete behavior under thermal conditions through a coupled hygro-thermalmechanical formulation.Srivastava and Prakash [15] developed a novel coupled framework for analysis of reinforced concrete and steel planar frames subjected to fire with three-way coupling between heat transfer, mechanical deformations and pore pressure build-up and used several numerical examples to demonstrate the accuracy and applicability of the framework.Padre et al. [3] implemented an algorithm to check the resistance of reinforced concrete sections to oblique unsymmetrical bending at ambient temperature and in a fire situation.Magisano et al. [16] proposed an automatic procedure for evaluating the axial force-biaxial bending yield surface of reinforced concrete sections in fire and a strategy to determine the limit fire duration, that is, the time of exposure which leads to structural collapse.Schulthess et al. [17] presented a method named hybrid fire testing and validated it with multiple proof-of-concept tests covering the entire temperature range relevant to structural fire engineering.Assis et al. [18] applied a thermomechanical model to assess concrete mechanical damage under high temperatures using Mazars' [19] theory and experimental data for validation.
Using other approaches, Liao and Huang [20] developed a robust finite element procedure for modeling the localized fracture of reinforced concrete beams at elevated temperatures and validated it against previous fire test results on the concrete beams.Nguyen et al. [21] studied the change in thermal conductivity of concrete when exposed to mechanical and thermal loading through a three-phase plane model using lattice discretization, where the damage variable is accounted for via crack width.They used numerical examples to illustrate and validate the proposal.Dias et al. [22] studied concrete under high temperatures and concluded that, in this situation, the material undergoes significant deterioration with spalling and a decrease in Young's modulus, compressive strength and durability.
There are also in the bibliography studies that deal with the influence of aggregate type and size in concrete under high temperatures behavior.Nince and Figueiredo [23], Kong and Sanjayan [24], Pan et al. [25] and Ali et al. [26] studied the relation between the aggregate size and degradation process of concrete under high temperatures, by observing the spalling of the structure superficial layers, and concluded that spalling increase and aggregate size are inversely proportional.On the other hand, Jansson and Bostrom [27] state that spalling extent is proportional to aggregate size.Souza and Moreno [28] investigated the strength decrease of concrete produced with different aggregates and exposed to 573.15K and 873.15K and found a large decay of compressive strength under the higher temperature value.Fanton [29] reviewed reinforced concrete slab behavior in a fire situation through finite elements software and concluded that concrete produced with limestone is better than one with siliceous aggregate in these situations.
Aiming to contribute to this subject, this work proposes a computational analysis of the aggregate size influence on the mechanical behavior of concrete under high temperatures, applying the Mazars' [19] damage model.Analyses were performed by representing the problem's geometry in 2 and 3 dimensions.For the material modeling, two different approaches were adopted: a two-phase medium, composed of coarse aggregate and mortar, and a three-phase medium, by adding the interfacial transition zone.

METHODOLOGY
In this work, a thermomechanical model in finite elements was implemented in Abaqus [30] and applied to assess the influence of aggregate size on concrete behavior when exposed to high temperatures.This model was developed using a weak coupling strategy in which thermal and mechanical analyses are performed separately, for the sake of computational cost, which allowed the use of a mesh with a satisfactory refinement level.
Analysis was performed in three groups of computational samples, concerning the problem and the material representations: two-phase in 2D, three-phase in 2D and two-phase in 3D.In each group, different aggregate grading was adopted so that computational samples were significantly distinct regarding aggregate size.
Synthetic cylindric concrete samples were generated via a Python script in which aggregate particles, with or without an interfacial transition zone, are randomly distributed over a concrete area or volume, concerning a two or a three dimensional representation, respectively, in order to reproduce the grading curve and phase proportion.
Damage in the material was evaluated by means of Mazars' [19] theory, which is not available in the Abaqus' [30] library.Thus, it was necessary to develop a user subroutine in Fortran, to describe the adopted constitutive model.

Governing equations
Concerning thermal analysis, it was considered only heat transfer by conduction.Thus, Equation 1gives the temperature field, provided the initial temperature field and boundary conditions: where  is the density,  is the specific heat,  is the temperature,  is the time,  is the thermal conductivity and  is a source or a sink.
Once the temperature field was obtained, it was applied as a thermal loading so as to solve a mechanical problem, being the relationship between temperature and deformation obtained from Equation 2: where  is the thermal expansion coefficient,  0 is the initial temperature,  is the identity matrix and   is the thermal deformation.
For the mechanical problem resolution, the Cauchy equilibrium equation (Equation 3) was solved with adequate boundary conditions: where  are the body forces and  is the stress tensor.
The linear elastic constitutive model was modified by a damage variable for mortar and interfacial transition zone.Classical Hooke's law, given by Equation 4, was considered for the aggregate: where  is the displacement,  is the strain and  and  are the Lamé constants with: and where  is the Young's modulus and  is the Poisson's ratio.The damage variable  affects directly Young's modulus through the relation given by Equation 7: with   being the damaged Young's modulus.In turn, damage variable  is obtained from linear combination of tension and compression components (  and   ), with weights   and   according to Equation 8: Tension and compression damage variables are both described by Equation 9: where the subscripts  and  refer to tension and compression, respectively,  0 ,  , and  , are the Mazars' [19] model parameters -extracted from uniaxial stress x strain curves via geometrical fitting procedures -and ̃ is the equivalent deformation given by Equation 10: (10) in which  ()+ is the positive principal deformation.

Material properties
The thermal and mechanical properties of each concrete constituent phase were obtained from specific codes or literature references.Concerning thermal properties, aggregate and mortar thermal expansion coefficients were considered as temperature linear function, derived from experimental data approximation available in Razafinjato [31].The specific heat of granite and the initial value of specific heat of mortar were obtained from NBR 15220-2 [32].However, the value for granite was considered constant and the value for mortar was adopted as a temperature linear function.Aggregate and mortar thermal conductivity were considered constants and obtained from NBR 15220-2 [32].Aggregate and mortar density were also considered constant but obtained from Razafinjato [31].In the absence of available data for the interfacial transition zone, the same thermal properties used for mortar were considered for this phase.Thermal properties for the three phases are shown in Table 1.Concrete and mortar Young's modulus were known experimentally from Razafinjato [31].Aggregate Young's modulus and aggregate and mortar Poisson's ratio were obtained from an inverse method.For this, a python script solved a mechanical model in which unknown properties were estimated, by trial and error, until the computational sample Young's modulus was adjusted to experimental values.For the interfacial transition zone, it was considered 50% of mortar Young's modulus and the same Poisson's ratio of mortar according to Ramesh et al. [33].Mechanical properties for the three phases are shown in Table 2. From experimental data of concrete Young's modulus for some temperatures [31], the parameters of Mazars' [19] damage model were obtained applying a method described in previous works [18] [34] and then applied to the computational sample to obtain Young's modulus for each evaluated temperature.
In the studied problem the temperature was monotonically increasing, that is, cooling was not considered.Thus, the observed compression level was low and, consequently, the estimated value for   was insignificant when compared to   .So   was neglected and parameters were reduced to  0 ,   and   .Obtained parameters for the three computational sample groups are shown in Table 3.

Geometric modeling
The analysis comprised three different geometric approaches to represent a cylindrical (150mm x 300mm) concrete sample: two-phase in 2D, three-phase in 2D and two-phase in 3D.Aggregate particles were considered as spherical, for simplification purposes.The synthetic samples were generated with the help of an algorithm developed by Bonifacio [35] in Python.The grading curves shown in Figure 1 and an aggregate's proportion of 40% were adopted as input data so as to randomly determine the particles' distribution in the mortar medium.Figure 2 and Figure 3 show, respectively, the 2D and 3D geometries adopted for the samples -in the 2D model, it was considered a quarter of the plane section while the 3D model consisted of one-eighth of the cylinder.Six grading curves, shown in Figure 1, were strategically chosen so as to generate significantly different synthetic concrete samples: curves 1, 2, and 3 were adopted for 2D samples and curves 4, 5 and 6 were considered as input data for 3D samples.The 2D geometries shown in Figure 2 adopted grading curves 1, 2 and 3, while Figure 3 shows 3D synthetic samples resulting from grading curves 4, 5 and 6.
The 2D geometries shown in Figure 2 were applied to simulate concrete as two-phase (mortar and aggregate) and three-phase (mortar, aggregate and interfacial transition zone).The interfacial transition zone (ITZ) was included around every aggregate particle of diameter 2 by means of a concentric circle of diameter 2 + 2 where  is the thickness of the interfacial transition zone and  is the radius.According to Mehta and Monteiro [36], the ITZ thickness varies from 0.01mm to 0.05mm.In the present work, it was adopted  = 0.05mm, for the 2D three-phase model.The 3D counterpart, however, was not performed due to the occurrence of mesh distortions related to the small thickness of the ITZ.
For the finite elements modeling, 2D samples were discretized with DC2D4 (4-node linear heat transfer quadrilateral) and DC2D3 (3-node linear heat transfer triangle) elements for thermal analysis and CPS4 (4-node bilinear plane stress quadrilateral) and CPS3 (3-node linear plane stress triangle) elements were used for mechanical analysis.DC3D4 (4-node linear heat transfer tetrahedron) elements were applied for 3D thermal analysis and C3D4 (4-node linear tetrahedron) elements were adopted for 3D mechanical analysis.

Thermomechanical model
This work adopts a weak coupling strategy to perform a thermomechanical analysis which consists of two models: a transient thermal one, in which the temperature field is obtained for a thermal boundary condition, followed by a mechanical one, with suitable boundary conditions, where the mechanical damage generated by thermal loading is evaluated.Figure 4 shows the boundary conditions in each face for the 2D and 3D cases.In the first model,  = 573.15K,723.15K, 803.15K was applied on exterior faces and  = 0 on internal faces, with the initial condition  0 = 293.15K.In the second model, a displacement restriction was applied on internal faces.Material damage was incorporated via the Mazars' [19] constitutive model, implemented in a Fortran subroutine named UMAT.In this case, software accomplishes the pre-processing and post-processing normally while the algorithm developed externally for the user is used for the processing.This subroutine is schematically shown in Figure 5.

RESULTS
Figure 6 shows the temperature field results referring to the maximal imposed temperature (803.15K) for the twophase in 2D model, which were quite similar to those related to the three-phase in 2D model.The results for the twophase in 3D geometry are shown in Figure 7, for a vertical section of the computational sample.It is noted in all cases that the highest temperatures are located on the faces in contact with the heat flow, as expected.Furthermore, Figure 7 indicates that aggregate diameters do not seem to have significantly affected temperature distribution.However, there is a higher average temperature in 3D geometries (approximately 788K), which is evidenced by the significantly higher minimum temperature compared to those verified for 2D geometries.Subsequently, this thermal loading was applied to the mechanical model with damage from which it was possible to obtain the damage map for the 2D cases, as shown in Figures 8 and 9.For the two-phase in 3D case, the results for a vertical section are shown in Figure 10.In all geometries, greatest damage is located close to the aggregates.In addition, for all cases, there is a high degree of damage in practically the entire section.It is also observed that for larger aggregates there are more regions with high damage values (red color) between the particles, although least damaged regions (yellow color) are observed in the computational sample as a whole.Graphs in Figure 11 present the average damage values obtained for each computational sample considered.Through them it is possible to quantitatively assess that the average damage is inversely proportional to the diameter of the aggregate.By considering that spalling is an important damage process, the results are in agreement with studies developed by Nince and Figueiredo [23], Kong and Sanjayan [24], Pan et al. [25] and Ali et al. [26].Still in Figure 11, more significant differences may be identified in the 2D cases and it is possible to notice that the temperature at which the damage starts is not highly influenced by the particle's size.
According to Figure 11, for the two-phase in 2D geometries, the highest average damage was 0.729 for curve 1 and the lowest was 0.711 for curve 3.For the three-phase in 2D geometries the highest average damage was 0.728 for curve 1 and the lowest was 0.711 for curve 3.For the two-phase in 3D geometries, the highest average damage was 0.740 for curve 4 and the lowest was 0.734 for curve 6.
Figures 12 and 13 show the Young's modulus map obtained from the mechanical model with damage for 2D computational samples.Figure 14 shows these results for the two-phase in 3D geometries for a vertical section.Damage is quantified by means of  decay.There is also a more uniform distribution of Young's modulus values in the sections with smaller aggregates, while in the sections with larger aggregates regions with more discrepant values are found.Finally, graphs in Figure 15 present the average Young's modulus values obtained for each computational sample considered.It is possible to quantitatively assess that this property is higher for larger aggregates, while it becomes relatively lower for smaller aggregates, which is explained by the inverse relationship between this property and damage.The graphs also show the experimental reference [31] applied for calibration purposes.It is possible to notice that results related to 2D computational samples (curve 1) are in better agreement with the experimental counterpart.Results obtained for the 3D analyses were practically coincident.Considering the maximal imposed temperature (T=803.15K),results denote that for the two-phase in 2D computational samples the highest average Young's modulus was 9868MPa for curve 3 and the lowest was 9286MPa for curve 1.For the three-phase in 2D computational samples, the highest average Young's modulus was 9906MPa for curve 3 and the lowest was 9350MPa for curve 1.Finally, for the two-phase in 3D computational samples, the highest average Young's modulus was 9073MPa for curve 6 and the lowest was 8867MPa for curve 4. It should be noted that the experimental reference value is 9000MPa.
As observed computationally, the thermal model was not influenced by the aggregate grading.The mechanical model with damage denotes that smaller aggregates lead to a higher average damage.Nevertheless, more elements with higher damage are observed in samples with large aggregates.In contrast, for small aggregates the damage is better distributed throughout the section.In relation to Young's modulus, due to its inverse relationship with damage, it is noticed that higher mean values are observed for large aggregates.
Regarding the sensitivity of the computational model, the 2D cases are apparently more sensitive to aggregate grading, while 3D geometry results did not indicate a significant influence of this parameter.
Concerning the 2D representation of concrete composition, results for the two-phase and three-phase geometries were quite similar.It is worth mentioning that the only aspect considered distinct between the samples of the same group was the size of the aggregates.Thus, the observed sensitivity was considered satisfactory and future adjustments in the thermal and mechanical properties, which undergo changes according to the particle size, will contribute to more distinct results.
The influence of aggregate grading on concrete damage is a controversial issue in the technical literature.A number of experimental researches lead to divergent results, while there is a lack of numerical studies on this topic.In such a context, the present work aims to contribute for a better knowledge on the subject, by means of a numerical tool capable of relating aggregate diameters to damage evolution in a concrete medium.

CONCLUSION
This work presents a thermomechanical model applied to the simulation of concrete behavior under imposed temperature raise.The material was represented as an heterogeneous medium, composed of aggregate-mortar and or aggregate-mortar-ITZ, mostly aiming to verify the analysis sensitivity to aggregate grading concerning damage evolution, which was evaluated according to Mazar's Model.
Results denote that aggregates' particle size influences the damage distribution in the medium.Synthetic concrete samples with smaller inclusions show a more homogeneous damage distribution in spite of a higher average damage than that verified in computational samples with aggregates of larger dimensions.
An important aspect of the proposed weak coupling strategy is the fact that it demands low computational cost and a small set of experimental data so as to provide information on the material's degradation under temperature exposition.
In this context, the average processing time of thermal models was 7.0 minutes for the two-phase in 2D, 26.8 minutes for the three-phase in 2D and 68.3 minutes for the two-phase in 3D.For the mechanical model, the average processing time was 0.5 minutes for the two-phase in 2D, 2.5 minutes for the three-phase in 2D and 10.3 minutes for the two-phase in 3D.In both cases a machine with a Intel(R) Core(TM) i7-7500U CPU @ 2.70GHz 2.90 GHz processor was used.
The sensitivity presented by the computational model encourages further improvements so as to contribute to the understanding of concrete's behavior under high temperatures.

Figure 1 .
Figure 1.Grading curves adopted for the 2D and 3D synthetic concrete samples.

Figure 2 .
Figure 2. Two-phase in 2D computational samples for curves 1 to 3, from left to right.

Figure 3 .
Figure 3. Two-phase in 3D computational samples for curves 4 to 6, from left to right.

Figure 4 .
Figure 4. Faces of the 2D and 3D models.

Figure 6 .
Figure 6.Temperature field in Kelvin for curves 1, 2 and 3 (two-phase in 2D), from left to right, at 803.15K.

Figure 7 .
Figure 7. Temperature field in Kelvin for curves 4, 5 and 6 (two-phase in 3D), from left to right, at 803.15K.

Figure 15 .
Figure 15.Average Young's modulus obtained in each computational sample evaluated for the two-phase in 2D case, three-phase in 2D case and two-phase in 3D case.

Table 1 .
Thermal properties, T being the temperature considered.