Finite Element Modeling of Concrete Behavior at Early Age Modelagem por Elementos Finitos do Comportamento do Concreto nas Primeiras Idades

© 2009 IBRACON a Pontifícia Universidade Católica do Rio Grande do Sul, maurenaurich@pucrs.br, Avenida Ipiranga, 6681, Prédio 30, Bloco B, Sala 101, 90619-900, Porto Alegre, RS, Brasil; b Universidade Federal do Rio Grande do Sul, americo@ufrgs.br, Av. Osvaldo Aranha, 99, 3o. Andar, 90035-150, Porto Alegre, RS, Brasil; c Escola Politécnica da Universidade de São Paulo, tulio.bittencourt@poli.usp.br, Av. Prof. Almeida Prado, Trav. 02, número 271, Butantã, 05508-900, São Paulo, SP, Brasil; d Center for Advanced Cement-Based Materials, Northwestern University, s-shah@ northwestern.edu, 2145 Sheridan Road, Evanston, Illinois


Introduction
Over the past several years, many numerical tools have been developed to evaluate the measures adopted to reduce the cracking risk at early age (Lura, 2000;Morabito, 1998;Shah, 1994). These tools are based on mathematical models capable to describe the coupling of thermal, moisture diffusion, chemical and mechanical analysis. In this work, according to Aurich (2008), a methodology for evaluate the potential cracking risk in concrete structures at early age is described. It will be shown numerical models to describe the main aspects of the chemical, thermal, moisture diffusion and mechanical phenomena that occur during the first days after the concrete cast and a computational implementation of these models through the finite element method. In the chemical analysis, the heat released due to the exothermic reactions of the cement hydration is determined by a concrete adiabatic temperature rise curve. This heat released due to the temperature rise inside the body is the input for the next step. In the thermal analysis, besides the heat released due to the exothermic reactions of the cement hydration, the heat flux on account of the difference of temperature between the body and the environment is considered. The nodal temperatures are evaluated from the body thermal and geometrical properties. These temperatures will be input data for the last step, the mechanical analysis. The follow analysis is the moisture diffusion. Taking into account the similarity between the differential equation that governs the heat transfer analysis and the moisture diffusion, the same procedures used to evaluate the nodal temperatures in the previous step are used to determine the nodal values of relative pore humidity. These values will also be used as input data in the mechanical analysis. In the mechanical analysis, the stress states due to the temperature and humidity changes evaluated in the previous steps and due to concrete shrinkage and creep are computed. When the stress state in any sample point reaches the failure surface, this point is considered cracked. The software considers concrete cracking through a smeared model with fixed crack. The developed software results were compared with experimental values found in the literature, demonstrating an excellent approach for all the implemented analysis.

Chemical Analisys
In this analysis, the heat released due to the exothermic reactions of the cement hydration is determined by a concrete adiabatic temperature rise curve. In this study, the curve proposed by the Japan Society of Civil Engineers (JSCE, 1999), as shown in figure 1, or a curve determined from experimental data could be used. According to JSCE (1999), the adiabatic temperature value, during the hydration reactions development, can be estimated according to equation 1.
where: • ∆T max : maximum temperature achieved by cement in the calorimetric test; • t: time (between 0 and 70 days). This analysis generates a volumetric load, in a perfect analogy with the structure self-weight, but instead of a force, it generates heat. This heat released due to the temperature rise inside the body is the input for the next step, the thermal analysis.

Thermal Analysis
According with Bathe (1996), for a heat transfer analysis in a twodimensional body, assuming that the material obeys the Fourier's Law of heat conduction, it can be written that: where: • q x , q y : heat flows conducted per unit area; • θ: body temperature; • k x , k y : thermal conductivities corresponding to the principal axes x,y. Considering the heat flow equilibrium inside the body, it is obtained that: where q B is the rate of heat generated per unit volume. On the surface of the body, the following conditions must be satisfied: a term that takes account of the rate at which heat is stored within the material. The rate of heat absorption may be considered as: where c is the material specific heat capacity. The variable q c can be understood as a part of the heat generated, which must be subtracted from the otherwise generated heat q B in the expression 7. This effect leads to a transient response solution.

Moisture Diffusion Analysis
In the diffusion most general case, according to CEB-FIP Model Code 1990 (1993), Kwon (2008)  where: • θ = θ(x, y, t): temperature as a function of location variables x and y and the time variable t; • k: thermal conductivity; • c: specific heat; • ρ: density. Similarly, rewriting the differential equation that governs the moisture diffusion phenomenon, it is transformed to: where: • H = H (x, y, t): pore relative humidity as a function location variables x and y and the time variable t; where: • θ S : known temperature at the surface S θ ; • k n : thermal conductivity of the body; • n: unit normal vector to the surface of the body; • q S : prescribed heat flux input on the body surface S q . Three boundary conditions may be considered in the heat transfer analysis: n Temperature conditions: the temperature may be prescribed or fixed on certain points or surfaces of the body, denoted by S θ in equation 4; n Heat flow conditions: the heat flow input may be prescribed in certain points or surfaces of the body. These heat flow boundary conditions S q are specified in equation 5; n Convection boundary conditions: the convection may be included as boundary condition, adding the following term in equation 5: where h is the convection coefficient. Here at environmental temperature θ e is known, but the surface temperature θ S is unknown. For the finite element solution of the heat transfer problem, using the principle of virtual temperatures, given as: where: and Q i is the concentrated heat flow input. Each Q i is equivalent to a surface heat flow, over a very small area. The bar over the temperature θ indicates that a virtual temperature distribution is being considered.
In the heat transfer problem considered above only it is assumed steady state conditions. However, when significant heat flow input changes are specified over a short time period (due to a change of any boundary conditions or the heat generation in the body), it is important to include where: • n ∂ θ ∂ : thermal gradient at the surface; • h: convective coefficient; • θ e : temperature at the surface of the solid; • θ S : temperature of the ambient space. Also, considering the boundary conditions in the moisture diffusion analysis, it can be written that: where: • n H ∂ ∂ : moisture gradient at the surface; • α m : mass transfer factor; • H e : relative pore humidity at the surface of the porous medium; • H S : relative humidity of the ambient space.
Comparing the equations that describe the moisture diffusion and heat transfer analysis, it is clear that exists an analogy between both problems. The correspondence between the parameters in the analyses can be seen in table 1.

Mechanical Analysis
Finite elements for axis-symmetrical, plane stress and strain analy-ses have been implemented. According to figure 2, the concrete is represented through a viscoelastic model, corresponding to a chain composed by five Maxwell´s elements in parallel. Bazant and Wu (1974) presented an algorithm to automatically determinate the parameters E µ (t) and η µ (t) for each age, t, of concrete, from experimental results or standard concrete values. According to Aurich (2008)  where α F is a coefficient that depends on the maximum aggregate size, also given in table 2.

Numerical Applications
Numerical applications of the described methodologies throughout Finite Element Modeling of Concrete Behavior at Early Age this work are presented in this item. The experimental results of concrete rings tests will be compared with the developed software results for the time at which the cracks occur and the increase of the crack opening. Grzybowski and Shah (1990) developed a test using a ring specimen to study the concrete cracking with restrained shrinkage. The test consists in a concrete ring internally limited by a steel ring, in which are placed strain-gages for strain measuring, and an external microscope to measure the crack width.
The tests proposed by the authors were conducted at ACBM laboratory (Center for Advanced Cement-Based Materials at Northwestern University, Evanston, Illinois) with the purpose to compare different types of shrinkage-reducing admixtures performances at early age. The dimensions of the specimens are given in figure 5. The specimens were cured for 4 hours at 20 o C and 100% RH, and then after demoulding exposed to drying in the humidity room at the ambient temperature of 20 o C and 40% RH. Figure 6 illustrates the test, the position of the strain-gages in the steel rings and concrete ring prepared for testing.
To analysis the algorithm results, this paper studies two rings tested by Shah, Karaguler and Sarigaphuti (1993) and Shah et al (1994). The results for temperature rising, moisture diffusion and stress development are presented. The comparison between the results for the time at which the cracks occur and the increase of the crack width is also showed. The mesh used for the finite element analysis has 168 elements and can be seen in figure 7. The boundary conditions adopted for the axis-symmetrical analysis (thermal, diffusion and mechanical) are shown in figures 8, 9 and 10. The parameters considered in the thermal, moisture diffusion and mechanical analysis for both studied rings are shown in tables 3, 4 e 5. Figure 11 shows the temperature rising for two nodes, P1, situated in the central zone of the concrete ring, and P2, situated in the external surface of the concrete ring. The temperature development maps are shown in figure 12. It can be seen that the maximum temperatures, due to the heat generated by the cement hydration reactions, occur in the central zone of the concrete ring. This situation was already verified in figure  11, where the node P1, situated in the central zone of the concrete ring, reached a higher temperature than the node P2, situated in the external surface of the concrete ring. In the fifth day, the concrete ring temperature had already equalized with the environment temperature. Figure 13 shows the moisture diffusion development for two nodes, P1, situated in the central zone of the concrete ring, and P2, situated in the external surface of the concrete ring. The moisture dif-Finite Element Modeling of Concrete Behavior at Early Age fusion development maps are shown in figure 14. It can be seen that the drying process, due to the environmental relative humidity is lesser than the initial relative humidity of the concrete ring, occur form the external surface to the central zone of the concrete ring. This situation was already verified in figure 13. The mechanical analysis results are shown for both rings tested, which properties were observed in table 5. The ring 2 (Shah, Karaguler and Sarigaphuti, 1993) were made with a concrete which compressive strength and Young modulus higher than ring 1 (Shah et al, 1994). Figure 15 shows, for ring 1 (Shah et al, 1994), the comparison between experimental and numerical results for the time at which the cracks occur and the increase of the crack width. Figure 16 shows the same comparisons for ring 2 (Shah, Karaguler and Sarigaphuti, 1993). The stresses development maps for ring 1 (Shah et al, 1994) are shown in figure 17. In these maps it could be observed that the crack opening occurs from the external surface to the interior of the concrete ring, to support what was observed in the test, according to figure 18. The agreement between the numerical and experimental results, for the time at which the crack occurs and the increase of the crack width, is excellent. The numerical analysis shows why ring 2, with higher strength, cracks earlier than ring 1. The concrete shrinkage prescribed deformations generate stresses proportional to Young modulus, and the tensile strength development at early age is not enough to avoid the cracking formation in the stiffener ring. It can be seen that, in the stresses development maps showed in figure 17, in the tenth day the red color appears in the external surface of the ring. This color corresponds to the stress level equivalent to the concrete tensile strength. So, the crack opens on the external surface, reducing the stress value. At the fourteenth day, the red stress level is reached again in a further internal position (the stress on external surface is close to zero). In sequence, with the crack propagation, the stress level falls down. At the nineteenth day, the stress peak is reached in a more internal position, being followed by the stress reducing. This circumstance successively repeats each time closer to the steel ring surface.

Conclusions
The results achieved, even considering the concrete behavior inherent variability, demonstrated that the adopted models could simulate the actual concrete behavior at early age. The implemented finite element model proved to be an accurate and appropriate tool to predict the concrete cracking in different concrete structures.