NUMERICAL PREDICTION MODEL FOR TEMPERATURE DISTRIBUTIONS IN CONCRETE AT EARLY AGES

A finite element-finite difference numerical model is developed for predicting non-uniform temperature development in hydrating concrete with respect to t ime and space. The results obtained from this model can be used by structural and construction engineers to predict critical thermal stresses induced due to differential temperatures between the core and the surface of the concrete at early ages and between t he zero-stress temperatures and the minimum equilibrat ing mbient temperatures that the concrete experien ces during its service life. The prediction of zero-str e s temperatures also enables to quantify the exten t of builtin curl developed in concrete structures. The finit e element is used to space discretization while the finite difference is used to obtain transient solutions of the model. The numerical formulations are then programmed in Matlab. The numerical results were co mpared with experimental results found in literatur e and demonstrated very good agreement.


INTRODUCTION
Non-uniform temperature caused by the hydration of cement, coupled with low conductivity of concrete can induce non-linear thermal gradients in the interior and exterior surfaces of concrete. This behavior is followed by non-uniform dimensional changes that will result in internal stresses (Branco et al., 1992). When such internal stresses exceed the concrete strength, which is relatively low at early ages, cracking may occur. The resulting cracks strongly affect the long term concrete performance as concrete will have to survive the early-age problems in order to guarantee a good long-term performance.
Cracking of concrete may also occur due to longterm critical thermal stresses, which are induced by the difference between the zero-stress temperature and the minimum concrete temperature that the concrete will be exposed to during its service life. Hence, proper determination of zero-stress temperatures also helps for a better design of reinforcement associated with thermal stresses and prediction of the extent of built-in curl.
The effects of non-uniform and zero-stress temperatures are highly pronounced on massive concrete structures such as dams, foundation slabs, concrete pavements and bridge decks. In these structures, temperature developments around the core follow nearly adiabatic scenarioswhile around the surfaces are highly influenced by the environmental conditions. NCHRP (2003) considers only monthly ambient temperature and cement content to compute the average zero stress temperature in the concrete cross-section. It also uses a constant with default value of-10°F for the permanent curl/warp effective temperature difference, which is the zero-stress temperature gradient. Lui and Fwa (2003) presented non-linear temperature development across concrete cross section is best described by a quadratic equation with unknown coefficients that can be derived from measured temperature data at top, middle and bottom of the concrete. Schindler et al. (2004) extensively researched on concrete temperature predictions and presented that the maximum concrete temperatures and the zero-stress temperatures appear to be similarly impacted by a change in the input variables affecting temperature development in hydrating concrete. Fig. 1. Factors affecting temperature prediction models along with the long-term and short-term critical temperature difference (Modified from Schindler et al., 2004) They also showed that the zero-stress temperatures are about 94 and 92% of the corresponding maximum temperatures under normal and hot weather concreting, respectively.

AJEAS
In this study, a finite element-finite difference temperature development prediction model is developed that accounts for the cement composition, the supplementary cementitious materials, the amount of cement, the fineness of cement, the water-cement ratio, the concrete placing temperature, the presence of mineral and chemical admixtures, the ambient air temperature, the wind speed, solar radiation, solar irradiation, concrete geometry, concrete size, types, lengths and area of the underlying and overlying materials. The finite element-finite difference numerical model is used to determine the temperature distributions in space and time. The zero-stress temperatures are determined from the corresponding maximum temperatures developed during the early ages. The numerical formulations are programmed using Matlab. The program results were then compared with experimental results published by Ballim (2004). Figure 1 summarizes the parameters involving in the temperature prediction models along with the critical long-term and short-term temperature differences.

Figure 1
Factors affecting temperature prediction models along with the long-term and short-term critical temperature difference (Modified from Schindler et al., 2004).
The Fourier equation governing the heat generation and temperature distribution for two dimensional elements with isotropic solid environment is expressed by (Aurich et al., 2009;Wang and Tian, 2005) Equation 1: where, Q • is the rate of the internal heat generated by cement hydration per unit volume, Tis the concrete temperature, k x and k y are the concrete thermal conductivity in the x-and y-directions, ρ is the density of the concrete and c p is the specific heat of the concrete. Thermal conductivity, specific heat and density of concrete have relatively constant values during and after the cement hydration but they are highly dependent on the type of aggregate used (Faria et al., 2006). Initial values, boundary conditions and rate of heat generation are to be determined before solving the governing equation numerically.
Concrete placement temperatures are considered as initial temperatures and given by Equation 2: where, T 0 (x, y) is the initial or placement temperature in the concrete region.
Mathematically, the concrete boundaries should satisfy either Drichlet or Cauchy boundary conditions. Drichlet, prescribed surface temperature, boundary conditions are given by Equation 3: Cauchy, heat flux, boundary conditionsare given by Equation 4: where, T P is the known value of the nodal points of the temperatures of the boundaries; q is flowing heat from surface or to the surface and its value is zero at insulated boundaries; h is the equivalent over all heat loss or gain coefficient including conduction, convection, radiation and irradiation; T S is the unknown temperatures at the boundary nodal points; T A is the ambient temperature; and n is the normal vector to the boundary.
Heat transfer between concrete and its surroundings takes place by thermal conduction, convection, radiation and irradiation. Thermal conduction is defined as heat transport in a material by transfer of heat between portions of the material that are in direct contact with each other. Conduction involves the exchange of heat between the concrete and underlying or overlying materials. The overall equivalent two dimensional conduction heat transfer coefficient for multiple layers of underlying and overlying materials is given by Equation 5: where, A is the unit area through which heat transfer is occurring and for each i th layer, L i is the thickness and k i is the conduction coefficient. The convective heat transfer can be expressed in terms of the Newton's cooling law Equation 6: where, q is the convective heat flux per unit area, h C is the convection coefficient, T S and T A are the surface and the air temperatures, respectively. The rate of heat flow from a horizontal surface is controlled by the magnitude of the temperature difference, the speed of the air flow and also the surface texture of the member. Based on the wind speed w (m/s), the convection coefficient, h C (kJ/m 2 /h/°C), can be estimated using the following Equation 7 The amount of solar radiation that reaches to the concrete surface can be estimated from the structure's geographical location, orientation, altitude, atmospheric conditions, time of the day and day of the year. Total radiation absorbed by the concrete surface can be determined from (Isgor and Razaqpur, 2004) where, α is the absorptivity of concrete, which varies between 0.5 and 1.0 and depends on the color and the texture of the surface and In is the direct solar radiation intensity which can be approximated using Equation 9: where, C N is the clearness number, which depends on the geographic location of the structure and can be obtained from meteorological map. The quantity β is the solar zenith angle, defined as the angle between the solar radiation and the normal to the horizontal plane at the location of the structure and is a function of the time of the day. G Z and τ are parameters quantifying the effect of the Julian day on the amount of radiation. Irradiation, also known as re-radiation, from the concrete surface can be obtained by using Stefan-Boltzman law Equation 10: where, ε is the emissivity of concrete, which is equal to the ratio of emission from a gray concrete surface to that from perfect radiator at the same temperature, (usually within the range of 0.85-0.95 and σ is the Stefan-Boltzmann constant (~5.67*10-8 W m -2 K -4 ). The non-linearity problem of the thermal irradiation heat transfer through the boundary can be bypassed by converting the non-linear irradiation heat transfer into the quasi-linear "irradiant convection" through a single convection irradiation coefficient given by (Faria et al., 2006) Because the heat transfer at the boundaries of a concrete can occur due to all convection, conduction and irradiation simultaneously, the overall heat transfer coefficient that includes all these mechanisms can be computed as Equation 13: Modeling of the ambient temperature that captures its diurnal variations is one of the challenges in developing the mathematical model for prediction of temperature development in concrete. At design stage of a construction project, it is unlikely that reliable and continuous ambient temperature values are available for the area in which the construction is to take place. However, daily ambient maximum and minimum temperatures are usually available from the local meteorological office. Hence, an approximate model that varies sinusoidally about the mean temperature was assumed to predict the ambient temperature as shown below (Ballim, 2004) where, t d is the clock time of day at which the prediction is being made (0 to 24 h); t m is the time at which the minimum overnight temperature occurs; T max and T min are the maximum and minimum temperatures for the day under considerations, respectively. In this study, adiabatic hydration model is used to simulate the heat of hydration generated during hardening of fresh concrete. When concrete is placed in adiabatic condition, the heat of hydration is completely converted into temperature and the adiabatic temperature rise of concrete is given by (Jafaar et al., 2007) where, K t maximum temperature rise of the concrete under an adiabatic condition is, α is a parameter accounting for the rate of heat generation and is the time (h). The cumulative heat generated per volume due to hydration up to time t is given by (Aurich et al., 2009): The expression for the rate of heat of hydration can be obtained by deriving Equation 16 with respect to time as Equation 17: Fairbairn et al. (2003) extensively presented adiabatic test results corresponding to different concrete mixes. Table  1 and 2 show the composition of the concrete mixes and their thermo-chemo-mechanical parameters, respectively. Figure 2 demonstrates the adiabatic temperature rise for the different concrete mixes which are used in this study as reference to take appropriate adiabatic parameters required for the numerical analysis. Each mix has unique combination of the adiabatic parameters depending on the water cement ratio, cement content, cement fineness, cement compositions, type and amount of supplementary cementitious material and type and amount of admixtures used.

MATERIALS AND METHODS
The numerical solution scheme used in this study is based on the Taylor-Galerkin finite element method in combination with finite difference method. A linear element is used for the one dimensional elements at the boundaries of the concrete and a bilinear element is used for the two dimensional elements inside the boundaries. Upon applying the numerical methods at the governing heat transfer equation, the following system of equation is derived: where, [K] is heat stiffness (conduction, convection, irradiation) matrix, [C] is the capacitance matrix, {f} is the total heat load vector due to heat generations and heat fluxes actions, {T} is the nodal temperature vector and is the rate of change of nodal temperature vector.
For time dependence problems, where temporal discretization is needed, it is effective to employ numerical solutions in the time domain to get the solution of the differential equation (Wang and Tian, 2005). From the mean value theorem for the differentiation Equation 19: where, ∆t is the time step and {T} t+∆t and {T} t are the vectors containing nodal temperature values at time t+∆t and t, respectively.
The vector {T} at time t = t * within the time step ∆t is given as: where, t t t * − θ = ∆ , the vector {f} is equal at t and t+∆t.

Re-writing Equation 18 in terms of {T} t and {T} t+∆t results in:
Equation 21 gives nodal values {T} t+∆t in terms of a set of known values, {T} t and the ratio,θ. There are four different types of finite difference schemes depending on the value of θ. In this study, the backward finite difference method, for which θ equals to 1, is used because a Fourier stability analysis showed that the numerical solutions of Equation 21 are unconditionally Science Publications AJEAS stable for θ ≥1/2 (Lewis et al., 1996). Hence, the Equation 21 reduces to: A Matlabprogram is developed to obtain numerical results of Equation 22.

RESULTS AND DISCUSSION
The main output of the temperature prediction program is the temperature of hardening concrete at each node of the finite element mesh at each time step. Figure  3 shows part of the mesh and the selected 21 nodes along the depth of the concrete section to evaluate temperature development with depth versus time. Table 3 shows the general material and geometrical input parameters used for evaluating the numerical model results based on a theoretical scenario. The scenario assumes that the concrete slab is placed under an ambient temperature varying sinusoidally between maximum and minimum daily temperatures of 35 and 15°C, respectively. The placement temperature is 25°C. The Cauchy's boundary coefficient for the top surface is taken as 5.15 whereas at the bottom of the concrete is about 2.15. The total heat loss/gain Cauchy's boundary coefficients are obtained by dividing the total heat loss/gain coefficients, h top and h bottom , by the conductivity of the concrete. This scenario emulates concrete pavements or slabs on grade. Figure 4 shows the plots of temperature development at thirteenodd number nodes along the depth of the concrete. Figure 4 confirms the expectation that the temperature developments at nodes around the core appeared to be higher but less affected by the ambient temperature variation than those nodes close to the top and bottom surfaces during the first couple of daysafter casting. Due to a higher heat transfer coefficient at the top surface than at the bottom surface, the effect of the ambient temperature appeared to be higher at the top nodes than at the bottom nodes. Referring to Fig. 4, the maximum temperature difference critical for early age cracking was found to be about 6°C, which was between nodes 9 (close to mid depth) and 21(at the top surface) and occurred around 13 h after casting time. It can also be seen that the maximum temperatures at the nodes occurred at different times in such a way that the lowest of the maximum temperatures occur first and the highest occur in the last. After about 72 h from the time the concrete was placed, the effect of heat of hydration became almost zero and the variations with depth are more or less followed the ambient temperature, which was assumed as varying sinusoidally in a 24 h period.   Figure 5 shows the zero stress temperature distributions along the depth of the concrete based on the findings that they are about 94% of the corresponding maximum temperatures shown in Fig. 4. It is worth noting that the zero stress temperature distribution is non-uniform and non-linear.
Finally, temperature results measured in the laboratory by Ballim (2004) were compared with the corresponding results from the numerical model developed in this study. Ballim (2004) measured the temperatures at different location after casting a 0.5 m 3 concrete block with a nominal compressive strength of 25 MPa. Figure 6 Shows the concrete block used by Ballim which was insulated on the two opposite 700 by1000 mm faces.

AJEAS
With the exception of the adiabatic parameters, the rest of the environmental and material parameters for executing the present model are taken from the literature including the thermal conductivity of the concrete (3.5 W/mK), heat capacity (1228 J/kgK), the heat transfer coefficients (30 W/m 2 K for exposed concrete surface at top and 5 W/m 2 K for the 700×700 mm end surfaces). The ambient temperature varied between 21-16°C during the five days when the test was conducted. Having related the properties of the concrete used for the test block and the adiabatric curves documented for different concrete mixes by Fairbairn et al. (2003) as shown in Fig. 2, it was reasonably determined that the adiabatic maximum temperature parameters is about 45°C and the adiabatic rate parameters is about 0.0375. The following three plots show the comparisons of the temperatures measured by Ballim (2004) and modeled in this study at three locations on the xy plane at the center of the block along the depth, which are at the bottom, at the mid depth and at 200 mm from the top of the block. As seen from Fig. 7, the overall results from the numerical model were in very good agreement to the experimental results with a maximum discrepancy of about 2°C at the mid depth at around 40 h after the concrete block was cast. In general, the maximum temperatures from the model at all the three locations are almost equal to the measured values. However, with regard to when these maximum temperatures occurred, the comparison plots show that results from the numerical model lag their corresponding from the experimental results. The maximum lag was about 6 h which was observed at the mid depth of the block.

CONCLUSION
This study presents a numerical model predicting early-age temperature development of concrete which incorporates the effects of cement type, cement content, water cement ratio, types and amount of admixtures (in terms of adiabatic parameters), environmental conditions, concrete thickness and curing conditions. This provides structural and construction engineers a better tool and understanding to estimate the potential of thermal cracking at early-ages during hydration as well as in the long-term when the concrete is exposed to cold temperatures. As expected, the numerical model shows that temperature developments of the nodes around the core are much less influenced to environmental variations than those nodes close to the surface. It was also observed that the zero-stress temperature distribution along the depth of the concrete is far from being linear.