Abstract
Thermal cracking of massive concrete structures occurs as a result of stresses caused by hydration in real environment conditions. The extended finite element method that combines thermal fields and creep is used in this study to analyze the thermal cracking of massive concrete structures. The temperature field is accurately simulated through an equivalent equation of heat conduction that considers the effect of a cooling pipe system. The time-dependent creep behavior of massive concrete is determined by the viscoelastic constitutive model with Prony series. Based on the degree of hydration, we consider the main properties related to cracking evolving with time. Numerical simulations of a real massive concrete structure are conducted. Results show that the developed method is efficient for numerical calculations of thermal cracks on massive concrete. Further analyses indicate that a cooling system and appropriate heat preservation measures can efficiently prevent the occurrence of thermal cracks.
1. Introduction
In massive hardening concrete, the thermal gradients under the heat of hydration during hardening can cause thermal cracking. The thermal cracking of concrete structures is a serious concern in construction. Some cracks with wide openings can particularly result in durability problems.
To address this concern, various numerical models have been proposed to investigate the nature of thermal cracks for massive concrete by means of finite element calculations. de Borst and van den Boogaard [1] studied deformation and cracking in early-age concrete modeled by finite element method (FEM). Xiang et al. [2] conducted thermoelastic analysis by examining time-dependent evolutions of heat release caused by hydration and material properties. Mazars and Bournazel [3], Waller et al. [4], and Lackner and Mang [5] performed thermal-elastic (or plastic) analyses of massive concrete without considering creep strain. de Schutter [6] simulated crack initiation and propagation by employing a simple stress-based cracking criterion and a smeared crack model based on the heat of hydration in hardening massive concrete elements. Benboudjema and Torrenti [7] developed a numerical model to predict early-age cracking of massive concrete through Kelvin-Voigt elements and an elastic damage model that considered autogenous and thermal shrinkage. Briffaut et al. [8, 9] examined the early-age cracking of massive concrete structures caused by thermal restrained shrinkage by using the thermal active restrained shrinkage ring test and FE methods. Buffo-Lacarrière et al. [10] proposed FE modeling to predict thermal cracking through the evolution of the damage variable law, which considered creep and damage behavior based on nonlinear viscoelasticity.
Massive concrete has been widely used in civil engineering, such as in the construction of dams, buildings, and bridges. Traditional finite element methods have been widely applied in the thermal analysis of massive concrete. However, conventional FEM is formulated with continuous media, so additional remeshing is necessary to accurately predict irregular crack propagation [11]. The extended finite element method (XFEM) [12–15] enhances conventional FEM capabilities by excluding the mesh requirement to conform to discontinuities. XFEM has been widely used in numerous fields with discontinuous problems, particularly in fracture mechanics, because XFEM is an excellent method of addressing discrete crack propagation in various types of materials. XFEM was applied to thermal problems in [16] and to shear band problems with thermal effects in [17]. It was applied by Duflot for the first time in steady-state thermoelastic fracture mechanics [18]. Zamani and Eslami [19, 20] investigated thermoelastic fractures by implementing XFEM for dynamic fractures and by improving accuracy through high-order crack tip enrichments. In these papers, the authors employed XFEM to characterize both the displacement discontinuity and the temperature discontinuity. In addition to XFEM, the embedded finite element method (EFEM) [21–23], meshfree methods (MMs) [24–27], and boundary elements (BEMs) [28–31] are also efficient computational methods based on the continuum-mechanics for discrete cracks.
An alternative pathway to simulate crack initiation and propagation offers methods based on the discrete-mechanics such as the particle simulation method [32–37] and the particle element method [38, 39]. The methods have been recently developed to discretize the continuum into particles so that crack generation at a contact between two particles is a natural process. Zhao et al. [40] simulated magma intrusion problems using solid particles for the surrounding rock and fluid particles for the intruded magma, respectively. The particle simulation method can efficiently simulate hydraulic fracturing because the fluid particles are allowed to move through the cracks between the solid particles. In addition, the particle simulation method has also been used to solve a broad range of scientific and engineering problems [41–43].
In accordance with the literature review, this paper aims to develop a numerical method based on XFEM to analyze thermal cracks in massive concrete structures. A practical thermal model that considers an artificial water-cooling system is adopted to simulate the temperature field. Time-dependent viscoelastic behavior is implemented for creep. As input parameters for XFEM that combines the temperature field and the viscoelastic constitutive model, the main properties develop over time when thermal cracking is calculated. The reliability of the developed method is verified by a case study on a lift of an arch dam, which underwent a cold wave during its construction. Subsequently, further analysis on thermal cracking is conducted.
2. XFEM for Thermal Cracking of Massive Concrete
2.1. Briefing on XFEM
Compared with the classical FEM, XFEM provides significant benefits in modeling crack propagation. The crack geometry in XFEM does not need to be aligned with the element edges. XFEM is based on the partition of unity concept and introduces additional degrees of freedom which are tied to the nodes of the elements intersected by the crack [12, 13].
XFEM displacement approximation around the crack is enriched with the step function to model the crack surface and with asymptotic displacement fields to model crack tips. The displacement approximation takes the following form [13]: where the first term is the classical FE approximation, is the FE shape function, and is the vector of the regular degrees of nodal freedom. The other terms are the enriched terms, where and stand for the set of nodes that have crack faces and crack tips in their support domain, respectively. and denote the vectors of additional degrees of nodal freedom for modeling crack faces and crack tips, respectively.
Step function simulates the displacement at the crack location with the help of the signed distance function to the crack, which is also called the normal level set [13, 15]:
The four following branch functions, which define the discontinuity near the crack tip, enrich the nodes with a distance to the front inferior to a prescribed enrichment radius: where and are polar coordinates in the local crack-tip coordinate system.
Duflot [18] applied the same procedure to temperature enrichment. The step function can account for the temperature jump. The leading term of the asymptotic field for the temperature of an adiabatic crack can be written as where coefficient provides the strength of flux singularity at the crack tip. The temperature field is discretized similar to the displacement field but with only the second branch function (4), which is the only discontinuous branch function. The approximation of the temperature field is then written as
As in standard FEM, performing numerical integration over the element domain is necessary to compute the element stiffness matrix. The subdomain method [44], in which the crack is one of the subdomain boundaries to carry out the numerical integration, is adopted for the elements that are cut by the crack.
2.2. Equivalent Equation of Heat Conduction
An isotropic homogenous domain bounded by boundary is considered. The thermal properties of this domain are independent of temperature, so the equation of thermal conduction is where is the heat flux, is the thermal conductivity, is the temperature, is the density, is the specific heat, and is the heat source.
For massive concrete structures, the hydration heat of concrete, which varies with age, is the main heat source. Temperature control is one of the critical problems in the construction of massive concrete, and embedding water-cooling pipes in massive concrete is one of the key methods to address this problem. Zhu [45] suggested that heat removal through a cooling system can be regarded as a negative heat source. Therefore, the heat source consists of the following two parts in consideration of the cooling pipe system: In the equation, the first part is the heat source caused by hydration, and the second part is the “cool source” caused by the water cooling system. is the adiabatic temperature rise of concrete at time , is the increment step, is the temperature of the cooling water inlet, and is the cooling function illustrated by Yang et al. [46]. Yang et al. [46] emphasized that although the temperature gradient near the cooling pipes cannot be obtained, uniformly dispersing heat throughout the solution domain is feasible.
For the adiabatic temperature rise of concrete, Zhu [47] proposed a compound exponential formula that is convenient for mathematical operation and is in accordance with experimental results. The hydration model can be written as follows: where and are the parameters that can be obtained through the optimization method or the trial method whose details can be found in the literature [47].
2.3. Thermal Boundary Conditions
When the concrete is in contact with air, we assume that the heat flux on the concrete surface is proportional to the difference between concrete surface temperature and air temperature . Therefore, the thermal boundary conditions at the concrete surface can be written as where is the surface heat convection coefficient and represents the normal direction of outer surfaces.
2.4. Viscoelastic Creep Model
Early-age concrete exhibits high creep, which reduces the stresses significantly more than any other influencing parameter, as well as assisting in mitigating thermal cracking by reducing thermal stress. Bažant [48] briefly summarized the method of calculating concrete creep.
The mechanism of concrete creep is not entirely clear, but the viscoelastic constitutive model composed of spring and dashpots can describe its performance. In this study, creep is included in time-dependent viscoelastic behavior. The constitutive relationship of viscoelastic materials can be represented in an integral form as follows [49]: where and are the shear and bulk moduli, respectively, which are functions of reduced time . Superscript dots denote differentiation with respect to time , and and are the mechanical deviatory and volumetric strains, respectively. The bulk () and shear () moduli can be defined individually through Prony series representation, which is expressed as where , , , and are material constants. Variable numbers of Prony series parameters end in and . and are the instantaneous relaxation moduli obtained by the instantaneous Young’s moduli as follows: where is the Poisson’s ratio.
To obtain the relaxation modulus, the normalized shear and bulk moduli can be defined as functions of the creep coefficient, which is expressed as follows [50]: where is the creep coefficient, which can be calculated as follows [51]: with where , , , , , and represent a series of correction coefficients that were described in detail in the literature [51]. Through nonlinear least-squares fit, the Prony series parameters in (11) are calculated.
The developed XFEM method that combines the equivalent equation of heat conduction and the viscoelastic creep model, as previously described, can be used in the numerical calculation of the thermal cracking of massive concrete. Therefore, considering hydration and creep, we can simulate a few discrete thermal cracks on a massive concrete structure.
3. Fracture Criteria and Evolution of Material Properties
The temperature gradient in massive concrete leads to tensile stress, which in turn causes crack initiation and evolution. In modeling thermal cracking evolution, an appropriate law on crack initiation and damage evolution should be defined. The mechanical properties of concrete are essential to predict the development of thermal stress in massive concrete elements.
3.1. Criterion for Initiation
As a possible first approximation, the maximum principal stress criterion can be applied to simulate the onset of thermal cracks. This criterion can be defined as follows: where is the maximum principle stress ratio and and denote the maximum principle stress and maximum allowable principle stress (tensile strength), respectively. The symbol represents the so-called Macaulay brackets, which enable to have an alternative value of 0 (when ) and (when ). The crack initiates when .
3.2. Criterion for Evolution
In this study, the analysis of crack growth is based on the approach of linear elastic fracture mechanics and is conducted with the framework of the extended-FEM method. The power law criterion proposed by Wu and Reuter [52] can be extended to 3D problems as follows [53]: where , , and are parameters with values of 1 to 2, respectively. , , and denote the strain energy release rate parameters of mode , mode , and mode components, respectively; , , and are the critical values of the strain energy release rate (the fracture toughness) for the three fracture modes. The law forms an energy envelope surface, and if the energy release rate exceeds this surface, the crack propagates.
3.3. Evolution of Material Properties
During the construction period in which concrete changes from an almost liquid state to a solid state, most of the mechanical properties rapidly vary in relation to the age of concrete. Among these properties, modulus of elasticity, tensile strength, and fracture toughness are the key parameters used in analyzing thermal cracking.
Degree of hydration is a fundamental parameter that describes concrete behavior during hardening. de Schutter and Taerwe [54, 55] proposed a degree-of-hydration-based description of these material properties as follows: where , , and are the Young’s modulus, tensile strength, and critical values of strain energy release rate at degree of reaction , respectively; and , , and correspond to these values when . , , , and are parameters that depend on the concrete composition.
A practical estimation of the degree of hydration at each moment can be expressed as the ratio between the accumulated heat of hydration at time and that at the end of reaction, which can be calculated through the adiabatic temperature rise: where is the degree of hydration at time and is the time at end of hydration reaction.
The coefficient of thermal expansion of concrete is another key parameter for thermal cracking analysis. Its measured values tend to be stable within two days [56], so a constant for it is adopted in this paper.
4. Numerical Examples
4.1. Introduction
A 286 m high parabolic double-curvature arch dam under construction is located in Southwest China at an altitude of 610 m above sea level (level of crest). The dam consists of 31 monoliths, each containing a series of lifts. The 10th lift of the 13th monolith was poured on September 5, 2009 and underwent a cold wave before the next lift was poured. This lift was selected for the thermal cracking analysis. The 9th lift was also used in our calculation as the lower boundary of the 10th lift.
Figure 1 illustrates the finite element mesh of the two lifts with cast corners (they were used to facilitate the construction and increase the rigidity of the bottom of the dam) at downstream; each lift (1.5 m high) contains four-layer elements. Each element in Figure 1 has an edge length of about 0.984 m and a thickness of 0.375 m. The numerical model includes a total of 15488 elements, in which 5544 elements (the 10th lift) are used in the XFEM mesh.
4.2. Boundary Conditions
During construction, the 10th lift of the monolith had been standing for 86 days before the next lift was poured. The increasing upstream and downstream surface of the poured concrete slabs and the top surface of the latest slab were assumed to be directly exposed to the environment, where heat transfer by convention occurred. Because the other lifts under the 9th lift and in the other monoliths were not considered, the boundary condition of the lower surface of the 9th lift was assumed to be adiabatic with fixed vertical displacement.
The calculations correspond to a period of 128 days (the duration between starting the 9th lift and starting the 11th lift), and the 10th lift started to be poured 42 days after the 9th lift was poured. The daily mean air temperature curve at the dam site, as shown in Figure 2, was used in this analysis.
4.3. Properties of Materials
According to [47], the hydration model in the engineering case analyzed in this paper may be written as follows: The equation is in accordance with the experimental data. Each concrete slab contains one-layer staggered heterogeneous iron cooling pipes in the middle for artificial cooling; related property values are illustrated in [46, 57].
Equations (18) to (19) show that the Young’s modulus, the maximum allowable principle stress, and the critical values of strain energy release rate can be written as time-dependent functions (, , and ) when (20) is used. In the simplest case, the fracture process can be isotropic, and parameters and . The main properties of concrete are presented in Table 1. Figure 3 shows the inclusion of the creep coefficient in viscoelastic material behavior during the cracking analysis for the 10th lift.
4.4. Results and Discussions
Accurate temperature fields are required to calculate thermal cracking. Figure 4 compares the temperature histories of the concrete in the calculated results and the actual measured temperature from thermometers buried at middle height of the lift. Overall, the calculated results obtained from the methods used in this study agree with the measured data. A sudden temperature drop occurred in the position of the buried thermometers five days after the cold wave, and an inflection point was observed in the computed curve two days in advance. These conditions can be related to the positions of the thermometers and the thermal conductivity.
Both elastic and viscoelastic models were applied in the stress analysis. Four elements of different heights at the center of the slab (see Figure 1) were selected, with the time-dependent maximum principal stress at the centroid of these elements shown in Figure 5. Overall creep can reduce the stress levels by approximately 20%. At the beginning, the bottom element bears low tensile stress levels because of the effects from the lower lift. The three other elements bear compressive stresses because of the restraint of thermal dilatation and the rise to their peaks seven days after the lift casting. After the temperature decreases in the core as a result of artificial water cooling, compressive stresses decrease gradually, and then, tensile stress occurs. Affected by the ambient temperature, the stresses fluctuate over time and increase suddenly as a result of the cold wave. The stress of the surface element that contacts with air has the largest varied amplitude and exceeds the tensile strength eventually.
The cold wave results in a thermal gradient (between the core and the skin of the lift), which causes tensile stress that is particularly serious on the surface. The cracks occur at discrete locations where the maximum principle stress exceeds the allowable values and propagates based on the criterion. The cracks simulated by XFEM and observed on the slab are presented in Figure 6. This figure illustrates that the numerical results are generally consistent with real crack occurrence. Most of the main cracks are initiated at the edges of the lift where high-temperature gradients occur, and this condition favors crack initiation. The evolution paths of the simulated cracks are almost close to the lines because of the isotropic assumption for the material and the crack evolution law. The discrepancy between the simulated and real crack directions in the middle of the lift surface is a result of excluding the other lifts (except the 9th and 10th) and the foundation in the calculation. In addition, crack crossing cannot be simulated because of the limitations of XFEM. The numerical simulation indicates that ambient temperature is among the main reasons for the thermal cracking of massive concrete structures.
(a)
(b)
The casting temperature obviously affects the cracking risk if no cooling system is used [58]. The concrete pouring temperature in the preceding analysis is 6.7°C, and then, we gradually increase the casting temperature to the ambient one (30.3°C). The maximum principle stresses with different pouring temperatures are presented in Figure 7. The result indicates that an increase in casting temperature does not exert significant effects on the stresses during the early stage because the embedded water-cooling system decrease the concrete temperature effectively. A negligible discrepancy is observed between these stresses at a later stage.
Heat preservation measures can be adopted during the intermission between the pouring of two lifts to reduce the thermal gradient of concrete. In this study, we used different heat convection coefficients to simulate different levels of thermal insulation effects. The stress evolution with a series of heat convection coefficients is illustrated in Figure 8. The stress levels decrease significantly with a decrease in the coefficients, particularly from 400 kJ/m·day °C to 200 kJ/m·day °C. In contrast to the case during the early stage, a small convection coefficient assumes a critical function when the cold wave arrives. The results demonstrate that appropriate heat insulation measures can effectively reduce the risk of thermal cracking.
5. Conclusions
This study presents an effective XFEM method to simulate the thermal cracking of massive concrete. The temperature field in the slab can be obtained by the equivalent equation of heat conduction, which includes the hydration heat of concrete and the effect of cooling pipe systems. The time-dependent viscoelastic behavior presented by the Prony series for concrete is included in the prediction of stresses to describe the stress relaxation that results from creep. The Young’s modulus, tensile strength, and critical values of strain energy release rate which evolve over time are the key parameters that affect the likelihood of thermal cracking in massive concrete. Thermal cracking of massive concrete can be described through a combination of XFEM with the temperature field and the viscoelastic creep model.
Compared with existing numerical methods to simulate thermal cracks on massive concrete, the developed XFEM method in this study has advantages in efficiently simulating multiple discrete thermal crack propagations without presupposing crack path and remeshing. Therefore, this method can be conveniently applied to analyze thermal cracking for massive concrete structures during the design and construction phases.
The precision of the developed method was verified by an engineering application. The simulated and actual results were in agreement, as shown in two comparisons: the comparison between the calculated time-dependent temperature curves and the measured ones and the comparison between the numerical simulated cracks and the actual observed crack occurrence. If the structure is embedded with a cooling system and the ambient temperature is appropriate, the casting temperature for concrete exerts a limited effect on thermal cracking risk that is different for the case without a cooling system. However, heat preservation measures can exert a significant effect on thermal crack control for massive concrete, particularly when a cold wave occurs. These measures deserve attention during the construction of massive concrete structures.
Acknowledgments
This study was financially supported by National Basic Research Program of China (2013CB035902), Research Project of State Key Laboratory of Hydroscience and Engineering of Tsinghua University (2011-KY-4), and the National Natural Science Foundation of China (51279087).