# Early-Age Cracking Analysis of a HVFA Concrete Structure Based on Thermo-Hygro-Mechanical Modeling Combined with XFEM

Authors: Wei Cui; Yanjun Luo; Shuang Peng
Journal: Advances in Materials Science and Engineering
 (2020)

---

## Abstract

Due to the low hydration rate of high-volume fly ash (HVFA) concrete at early age, the temperature gradient between the concrete core and surfaces could be effectively reduced. However, the low hydration rate results in a lack of hydration degree for early-age HVFA concrete. Thus, during curing, compared to the strength of ordinary Portland cement (OPC), a subsequent lower one of HVFA concrete leads to a more sensitive response to inner stresses induced by thermal and moisture loads. Based on ABAQUS, in this paper, user subroutines in the temperature and moisture fields were developed, with regard to the hydration degree, for simulation of the temperature and moisture influences on concrete. Additionally, the Double Power Law (DPL) model was used to depict early-age deformations of concrete in the mechanical field. Combined with the extended finite element method (XFEM), another subroutine for early-age cracking analysis in the mechanical field was then developed. Together with aforementioned subroutines, a thermo-hygro-mechanical model is derived. For evaluation of early-age cracking initiation and propagation of a pier composed of HVFA concrete, the model was implemented with XFEM. The obtained results show that (1) temperature and restraints are the main causes of cracking; (2) moisture loss affects surface cracks on structures at early age; and (3) although the temperature difference between the core and surfaces is not much obvious due to the reduced rate of heat dissipation from hydration, cracking of early-age HVFA concrete is still likely to happen for its low early strength. Thus, timely curing is critical to prevent early cracking.


---

## Body

## 1. Introduction

Fly ash is an inorganic, noncombustible waste product of coal-burning power plants. The incorporation of high-volume fly ash (HVFA) in concrete not only mitigates the environmental pollution of fly ash but also improves the workability of concrete mixtures by suppressing bleeding and reducing the rate of heat dissipation from hydration [1, 2]. Due to the low heat dissipation rate, HVFA concrete is regarded as a better choice to avoid cracks caused by significant temperature gradients between core and surfaces. However, the corresponding low hydration rate also leads to lower early strength. Thus, cracking behavior for HVFA concrete is worthy of study.The driving forces of crack-causing stresses for early-age concrete casts are combined influences of thermal dilation, shrinkage (including drying shrinkage and autogenous shrinkage), and creep at early age under restrained boundary conditions [3–5]. Up to now, numerical efforts have been taken for evaluation of early-age concrete cracking, ranging from constitutive models to coupled methods. The DPL model is recognized as a powerful tool of high-precision to simulate the creep behavior of early-age concrete [6, 7]. Lee and Kim implemented a numerical analysis performed with the hydration-based microplane model and successfully simulated the typical cracking patterns due to edge restraint in a concrete slab [8]. Yuan and Wan completed a numerical procedure based on a micromechanical model and empirical formulas with consideration of creep, hydration, and moisture transport [3]. Yoneda et al. proposed a local tensile strength model treating micropore water pressure as prestress and applied the model to cracking criterion [9].However, to realistically simulate early-age concrete behavior after casting under complicated environmental conditions and stress development under restrained boundary conditions; none of the abovementioned models adequately substitute for the thermo-hygro-mechanical coupling method [7, 10, 11]. Ulm and Coussy modelled the mechanics and hardening procedure of early-age concrete in terms of coupling with thermal and chemical factors [12, 13]. This was also addressed by Cervera et al. with respect to hydration, aging, damage, and creep in terms of a thermo-chemo-mechanical model [14, 15]. Based on those models, a chemo-plastic model was proposed and implemented onto massive concrete structures by Lackner and Mang [16]. The humidity distribution was taken into account by Gawin et al. with all changes of material properties correlated with hydration degree [17, 18]. And this was further looked into by Di Luzio and Cusatis [19, 20], and by one-way coupling, the solidification-microprestress theory, and the microplane model, early-age concrete behavior was retrieved [21]. In this paper, by adopting the equivalent age of concrete and mechanical properties related, as well as the DPL model for creep, the thermo-hygro-mechanical coupling has been achieved for early-age HVFA concrete behavior.Conventional finite element methods are considered powerful tools for the analysis of smeared cracking. However, due to the complication of mesh alteration when cracking, it is very unlikely to model the real propagation of a crack by finite element methods. This is also attributed to that the singularity of the crack tip could only be estimated. Compared with conventional finite element methods, the extended finite element method (XFEM) is known as an effective crack-simulating tool; it is based on the concept of unity partition to enrich the conventional finite element approximation, to take into consideration the effects of singularity or discontinuity around a crack [22]. The main advantage of XFEM is its capability in modeling discontinuities independently, allowing the model to be discontinuous while cracking [23–26].This paper presents a detailed principle description of the thermo-hygro-mechanical model. By use of ABAQUS user subroutines, numerical analyses with XFEM were performed for verification of the model. Different from evaluation of stress development for prediction of the cracking patterns, the combination of the thermo-hygro-mechanical model and XFEM allows cracks to propagate simultaneously as the stress develops. For numerical insights into HVFA concrete structures, a case study of a freshly cured pier was looked into. The subsequent results show a satisfying agreement with the observed cracks.


## 2. Finite Element Analysis Method

### 2.1. Thermo-Hygro-Mechanical Model

#### 2.1.1. Thermal Submodel

As also employed in [27, 28], the Arrhenius law is found applicable to concrete hydration reaction for generation rate of heat against temperature, which is written as the following equation:(1)kT=A⋅e−EaT/RT,where A is the pre-exponential factor, Ea is the activation energy related to temperature, T is the absolute temperature, and R is the Boltzmann’s constant.It is suggested in [8] that for concrete, Ea = 33500 J/mol for T ≥ 293 K and Ea = 33500 + 1470 × (293-T) J/mol for T < 293 K. According to [28], in this paper, an approximate value for Ea/R of 4000 K is used.By adopting the Arrhenius law, Azenha et al. proposed to determine the maturity degree of concrete by its temperature and curing time [27]. As expressed in equation (2), the maturity degree, characterized as the equivalent age te, is a function of realistic age and temperature in history:(2)te=∫0teEa/R1/T0−1/Tdt,where T0 is the reference temperature (generally T0 = 293 K) and t is the actual curing time. The form of integral is discretized into accumulation of small increments for finite element analysis.Given that the adiabatic temperature rise (θ) is linearly correlated to generated hydration heat (Q), the hydration degree is defined in this paper as follows [29]:(3)αte=QteQmax=θteθmax,where θ(te) is determined empirically by the following equation [4]:(4)θte=Tmax×1−0.57e−0.28te+0.43e−0.004te,where Tmax = 39°C is adopted from Figure 1.Figure 1
Adiabatic temperature rise and autogenous shrinkage strain development with time.According to the work of de Schutter, the relationship between thermal conductivity and hydration degree could be expressed linearly [30]. And equation (5) is adopted for evaluation of the current thermal conductivity:(5)λα=λmax×1.33−0.33α,The specific heat (c) of early-age concrete is proposed by van Breugel as the following equation [31]:(6)c=Wcαccef+Wc1−αcc+Waca+Wwcwρ,where c and W are the specific heat and mass with subscripts c, a, and w representing cement, aggregates, and water. ρ is the density of concrete, and ccef is the assumed specific heat by ccef=0.0084Tc+0.339.By substituting the equivalent age-dependent parameters of equations (4)–(6) into the heat conduction differential equation, the equilibrium of heat within unit volume of concrete is retained as equation (7), also widely employed in [14, 17, 19]:(7)ρcα∂T∂t=∇λα∇T+ρcα∂θ∂t,where c, λ, α, and θ are functions of te as mentioned before, ▽ is the Laplace operator, and t is the curing age of concrete.


#### 2.1.2. Moisture Submodel

With regard to concrete structures that are exposed to environmental air at early age, the moisture content in concrete decreases due to the diffusion of moisture. According to Fick’s second law, the rate of moisture transfer per unit area in a certain direction is proportional to the gradient of the moisture concentration in that direction [32]. In analogy to equation (7), by taking advantage of the mass conservation differential equation, the moisture diffusion is controlled by equation (8), also addressed for in [17, 19]:(8)∂h∂t=∂∂xDxh∂h∂x+∂∂yDyh∂h∂y+∂∂zDzh∂h∂z+∂QT∂t+∂Qh∂t+∂QC∂t,where t is the time, h is the relative humidity in concrete, and Dx(h), Dy(h), and Dz(h) are moisture diffusion coefficients in the direction of the x, y, and z-axes, respectively. The subscripts T, h, and C denote diffusion flow induced by temperature gradient, hydration reaction, and cement carbonation. With the consideration that early-age concrete is often composed of high moisture degree in a relatively stable environmental at a constant temperature. The effects of temperature gradient and hydration reaction on moisture distribution are overlooked in this paper. Also given that significant cement carbonization occurs after quite a long time, the last term is omitted. Thus, above all, it is assumed in this paper that the moisture distribution is primarily described by the moisture diffusion coefficients and that the moisture diffusion coefficients are uniformly identical to D(h).By adopting the S-shaped form by Bažant and Najjar and in tandem with the principle of equivalent agete, the moisture diffusion coefficient D(h) is defined by the following equation [33]:(9)Dh=D1⋅α+1−α1+1−h/1−hcn⋅fte,where α=D0/D1, D1 is the maximum of D(h) for h = 100%, D0 is the minimum of D(h) for h = 0%, hc is the pore relative humidity, denoting the humidity drops half-way between D1 and D0, and f(te) is the function of equivalent age. In this paper, α = 0.05, hc = 0.8, and n = 15 are employed. And D1 is estimated from equation (7):(10)D1=D1.0fck/fck0,where D1.0 = 3.6 × 10−6 m2/h, fck0 = 10.0 MPa, and the characteristic compressive strength fck is estimated by mean compressive strength fcm through fck = fcm−8.0 MPa. In analogy to the definition of te, f(te) is characterized as follows:(11)fte=expEadR1T0−1T,where Ead = 35 kJ/mol is the activation energy for moisture diffusion.For boundary conditions to be covered below, the modified Menzel’s equation gives the surface moisture diffusion coefficient (Dsurf) which is as follows [3]:(12)Dsurf=A0.253+0.06Vah−he,where A is an empirical coefficient, determined by the w/c ratios, Va is the average wind speed, and h and he are the relative humidity of the surface and environment, respectively.


#### 2.1.3. Initial and Boundary Conditions

For temperature analysis, the initial condition isT(x, y, z, 0) = Tcast, in which Tcast is the casting temperature. As illustrated in [34, 35], the boundary conditions are chosen as follows:(13)−λ∂T∂n=βT−Te,where n is the normal direction of the surface, T and Te are the absolute temperature of the surface and environment, and β is the heat emission coefficient.For humidity analysis,h(x, y, z, 0) = 100% is chosen as the initial condition. For surfaces confronted with air, the boundary condition is shown as follows:(14)−D∂h∂n=Dsurfh−he,where D and Dsurf are coped by equations (9) and (12).


#### 2.1.4. Mechanical Submodel

The development of the short-term mechanical properties of early-age concrete, such as tensile strength, compressive strength, and Young’s modulus, is strongly related to the degree of hydration, which depends on both curing age and temperature [36]. The mechanical analyses are conducted with respect to the maturity concept, by adopting the inputs from the thermo-hygro submodels, i.e., the equivalent age te and the potential thermal and shrinkage strains.The total strain of the concrete structure at timeτ is expressed as(15)ετ=εeτ+εcτ+εTτ+εsτ,where εeτ is the linear elastic strain caused by stress, εcτ is the creep strain depending on stress and concrete age at loading, εTτ is the deformation due to temperature fluctuation, and εsτ is the deformation due to shrinkage including drying shrinkage and autogenous shrinkage.The basic creep of concrete is accounted for by the DPL, which has reasonably good performance for either early-age or long-term spans [6]. Note that the DPL formulation is expanded into the Taylor series for discretized analysis. Detailed derivations are found in [37]:(16)εcτ=∫0tJt,τdστ,where(17)Jt,τ=1Eτ+qEττ−dt−τp,where t is the curing age, τ is the concrete age at loading, E−1τ is the modulus of elasticity at loading age, t−τ is the continuous time under load, and p, d, and q are material parameters calibrated by experimental creep tests. According to [37], q = 0.6, d = 0.2, and p = 0.2.The increments of thermal and shrinkage strains are evaluated by the following equations:(18)ΔεTτ=αTΔT,(19)Δεsτ=ksΔh.As addressed for in [38], ks=−ε0sEt0/Ete. ε0s is the final shrinkage strain, and E(t0) and E(te) are the elastic modulus at the initial time and the equivalent age te.To account for the evolution of concrete elastic modulus, tensile strength, and compressive strength during early cement hydration, the following expressions are adopted [6]:(20)Ete=E28⋅es1−672/te−t0nE,ftte=ft28⋅es1−672/te−t0nt,fcte=fc28⋅es1−672/te−t0,where E28, ft28, and fc28 are the elastic modulus, tensile strength, and compressive strength at the age of 28 days, respectively, and 27.45 GPa, 1.78 MPa, and 17.83 MPa are used for E28, ft28, and fc28 (obtained from Figure 2). t0 = 10 hrs is the time when the strength and stiffness are defined to be zero; s, nt, and nE are parameters that account for the cement type used and are fitted by experimental data. 0.211, 0. 59, and 0.37 are selected here, respectively.Figure 2
Development of concrete strength with time.


### 2.2. Extended Finite Element Method

The formulation of the standard displacement function is based on the finite element method and is enriched by the Heaviside function and the crack tip enrichment functions [23]:(21)ux=∑i∈NNixui+∑j∈NdisNjHxai+∑k∈NasyNk∑a=14Φaxbka,where j is the set of nodes whose shape function support is cut by a crack, i is the set of all nodes of the domain, Ni is the shape function associated with the node i, K is the set of nodes whose shape function support contains the crack front, and ui, ai, and bka are the displacements associated with the corresponding nodes.The Heaviside jump function,Hx, is a discontinuous function across the crack surface and is constant on each side of the crack [26]:(22)Hx+1above the crack−1beneath the crack.The crack tip enrichment functions,Φax, are obtained from the asymptotic displacement fields [22, 23]:(23)Φar,θ=rsinθ2,rcosθ2,rsinθsinθ2,rsinθcosθ2.The cohesive property of the crack tip is considered in this paper, and that means the stresses at the crack tip are not singular. [39, 40] carried out an asymptotic analysis for very large structure of the mechanical fields in a cohesive zone. It is better to consider the following functions [41]:(24)Φar,θ=rsinθ2orr3/2sinθ2orr2sinθ2.Constitutive formulation about the cohesive crack model can be seen in [24, 41, 42] while formulation information of XFEM can be found in [22, 23, 43].Gauss integral is employed for the numerical integration. For simplification, the stresses of Gauss integral points are used to control the crack propagation. Elements fail when maximum principal stresses of the Gauss integral points exceed the tensile strength.


### 2.3. Submodel Coupling

Fundamentally, concrete structures are under complex states of stress after being placed at construction sites due to the stress-inducing mechanisms of temperature effects, moisture effects, mechanical effects, and their combinations. For simplification, several assumptions shall be assumed between submodels. Further information regarding coupling methods is addressed for in this section.As shown in Figure3, the temperature, moisture, and mechanical fields are connected in a certain order via one-way coupling. Temperature analysis is conducted first, in which the hydration procedure for HVFA concrete is simulated according to equation (4); both element temperature and the corresponding hydration degree are obtained.Figure 3
Flowchart of the modeling process.Even though moisture diffusion affects temperature analysis to some extent, it is small enough to be omitted. Moisture diffusion properties are, however, strongly dependent on temperature. Accordingly, temperature analysis is initially conducted before moisture computation. Then, the thermal results are used as inputs for calculation of the moisture distribution under temperature effects.With respect to mechanical fields, temperature evolution will not only lead to volumetric changes but also influence the development of the concrete properties based on the equivalent age. Moreover, drying shrinkage induced by obvious surface moisture gradients via moisture transport to the environment will lead to shrinkage stress under restraint boundary conditions, thereby significantly affecting the stress responses in the mechanical field.The mechanical field also influences the temperature and the moisture fields upon two reasons: (1) the inner transportation of temperature and moisture are subsequently affected by the existence and propagation of cracks in the concrete structure and (2) the variations of concrete volume will slightly influence the boundary conditions of the temperature and moisture fields. Though it may lead to inaccuracy, it is assumed to be omitted the effects of cracks on temperature and moisture fields, which is characterized in Figure3 as a one-way coupling method. For more realistic results, a two-way coupling method could be expected to be implemented afterwards.Consequently, the mechanical field is the last step in the one-way coupling method, where the results of the previous two fields (temperature and moisture) are used as inputs to analyze the stress development and corresponding crack propagation. In this manner, relatively accurate results of early-age concrete behavior are obtained. By use of ABAQUS user subroutines, the specific implementation is as follows:Temperature Analysis. User subroutine UMATHT is firstly called by use of adiabatic temperature rise equation and other temperature parameters, such as thermal conductivity and specific heat. At the meantime, the corresponding initial temperature is imported together with reasonable boundary conditions (ambient temperature and heat transfer coefficient) in terms of user subroutine FILM. The obtained results are saved and named in text documents in the order of the corresponding element numbers at the end of the computation.Moisture Analysis. Similar user subroutines are used herein for the humidity distribution. UMATHT is called to input controlling equations to simulate variations of the moisture diffusion coefficient, and FILM is used to define the surface moisture diffusion coefficient and ambient relative humidity. Since moisture diffusion coefficients are dependent on temperature effects, a self-defined procedure is adopted to read the temperature data, and then the texts are saved in the form of vectors for moisture distribution calculation. In the end, the results are collected in the same manner as the temperature analysis.Mechanical-Cracking Analysis. A user subroutine UMAT was coded for mechanical analysis. As expressed by equation (15), the changing distribution of temperature and humidity is translated as additional strains, in tandem with the creep strain, into the stress field after the former coupling analysis. Thermal expansion coefficient is adopted for coupling of the temperature-mechanical field, while taking into consideration the development of mechanical properties, such as the modulus of elasticity, compressive strength, and tensile strength, over the equivalent age. As indicated by equation (19), drying shrinkage strains are gained from moisture distribution, which is superposed onto the stress field together with the thermal strains. Hence, the thermal, moisture, and creep effects are considered altogether, and the corresponding stress is computed.The XFEM method employed in this paper has been built in ABAQUS with a user-defined cohesive law. Prior to cracking, the initial response is controlled by the preceding UMAT. The damage criterion and evolution are, respectively, determined by the current tensile strength as equations (20) and (25) and fracture energy controlled traction-separation response as equation (26).(25)f=σmaxftte,where  is the Macaulay bracket and σmax represents the maximum principal stress. Damage is assumed to initiate when the value of the principal stress ratio, i.e., f, reaches a value of one.(26)tn=1−DTn,Tn≥0,Tn,otherwise,tJ=1−DTJ,J=s,t,where D represents the averaged overall damage at the intersection between the crack surfaces and the edges of cracked elements, evolving initially from 0 to 1. D is calculated by ABAQUS in the form of an exponential softening procedure with the Power Law.In practice, the UEL subroutine is adopted to combine XFEM with the DPL model to simulate crack initiation, evolution, and the corresponding stress development. The fracture energy was empirically estimated at 40 N/m for the normal mode and 10 N/m for the shear mode, in line with the values suggested by the ABAQUS documentation. And the contact properties between crack surfaces are specified as Penalty for tangential behavior and “Hard” Contact for normal behavior.


## 2.1. Thermo-Hygro-Mechanical Model

### 2.1.1. Thermal Submodel

As also employed in [27, 28], the Arrhenius law is found applicable to concrete hydration reaction for generation rate of heat against temperature, which is written as the following equation:(1)kT=A⋅e−EaT/RT,where A is the pre-exponential factor, Ea is the activation energy related to temperature, T is the absolute temperature, and R is the Boltzmann’s constant.It is suggested in [8] that for concrete, Ea = 33500 J/mol for T ≥ 293 K and Ea = 33500 + 1470 × (293-T) J/mol for T < 293 K. According to [28], in this paper, an approximate value for Ea/R of 4000 K is used.By adopting the Arrhenius law, Azenha et al. proposed to determine the maturity degree of concrete by its temperature and curing time [27]. As expressed in equation (2), the maturity degree, characterized as the equivalent age te, is a function of realistic age and temperature in history:(2)te=∫0teEa/R1/T0−1/Tdt,where T0 is the reference temperature (generally T0 = 293 K) and t is the actual curing time. The form of integral is discretized into accumulation of small increments for finite element analysis.Given that the adiabatic temperature rise (θ) is linearly correlated to generated hydration heat (Q), the hydration degree is defined in this paper as follows [29]:(3)αte=QteQmax=θteθmax,where θ(te) is determined empirically by the following equation [4]:(4)θte=Tmax×1−0.57e−0.28te+0.43e−0.004te,where Tmax = 39°C is adopted from Figure 1.Figure 1
Adiabatic temperature rise and autogenous shrinkage strain development with time.According to the work of de Schutter, the relationship between thermal conductivity and hydration degree could be expressed linearly [30]. And equation (5) is adopted for evaluation of the current thermal conductivity:(5)λα=λmax×1.33−0.33α,The specific heat (c) of early-age concrete is proposed by van Breugel as the following equation [31]:(6)c=Wcαccef+Wc1−αcc+Waca+Wwcwρ,where c and W are the specific heat and mass with subscripts c, a, and w representing cement, aggregates, and water. ρ is the density of concrete, and ccef is the assumed specific heat by ccef=0.0084Tc+0.339.By substituting the equivalent age-dependent parameters of equations (4)–(6) into the heat conduction differential equation, the equilibrium of heat within unit volume of concrete is retained as equation (7), also widely employed in [14, 17, 19]:(7)ρcα∂T∂t=∇λα∇T+ρcα∂θ∂t,where c, λ, α, and θ are functions of te as mentioned before, ▽ is the Laplace operator, and t is the curing age of concrete.


### 2.1.2. Moisture Submodel

With regard to concrete structures that are exposed to environmental air at early age, the moisture content in concrete decreases due to the diffusion of moisture. According to Fick’s second law, the rate of moisture transfer per unit area in a certain direction is proportional to the gradient of the moisture concentration in that direction [32]. In analogy to equation (7), by taking advantage of the mass conservation differential equation, the moisture diffusion is controlled by equation (8), also addressed for in [17, 19]:(8)∂h∂t=∂∂xDxh∂h∂x+∂∂yDyh∂h∂y+∂∂zDzh∂h∂z+∂QT∂t+∂Qh∂t+∂QC∂t,where t is the time, h is the relative humidity in concrete, and Dx(h), Dy(h), and Dz(h) are moisture diffusion coefficients in the direction of the x, y, and z-axes, respectively. The subscripts T, h, and C denote diffusion flow induced by temperature gradient, hydration reaction, and cement carbonation. With the consideration that early-age concrete is often composed of high moisture degree in a relatively stable environmental at a constant temperature. The effects of temperature gradient and hydration reaction on moisture distribution are overlooked in this paper. Also given that significant cement carbonization occurs after quite a long time, the last term is omitted. Thus, above all, it is assumed in this paper that the moisture distribution is primarily described by the moisture diffusion coefficients and that the moisture diffusion coefficients are uniformly identical to D(h).By adopting the S-shaped form by Bažant and Najjar and in tandem with the principle of equivalent agete, the moisture diffusion coefficient D(h) is defined by the following equation [33]:(9)Dh=D1⋅α+1−α1+1−h/1−hcn⋅fte,where α=D0/D1, D1 is the maximum of D(h) for h = 100%, D0 is the minimum of D(h) for h = 0%, hc is the pore relative humidity, denoting the humidity drops half-way between D1 and D0, and f(te) is the function of equivalent age. In this paper, α = 0.05, hc = 0.8, and n = 15 are employed. And D1 is estimated from equation (7):(10)D1=D1.0fck/fck0,where D1.0 = 3.6 × 10−6 m2/h, fck0 = 10.0 MPa, and the characteristic compressive strength fck is estimated by mean compressive strength fcm through fck = fcm−8.0 MPa. In analogy to the definition of te, f(te) is characterized as follows:(11)fte=expEadR1T0−1T,where Ead = 35 kJ/mol is the activation energy for moisture diffusion.For boundary conditions to be covered below, the modified Menzel’s equation gives the surface moisture diffusion coefficient (Dsurf) which is as follows [3]:(12)Dsurf=A0.253+0.06Vah−he,where A is an empirical coefficient, determined by the w/c ratios, Va is the average wind speed, and h and he are the relative humidity of the surface and environment, respectively.


### 2.1.3. Initial and Boundary Conditions

For temperature analysis, the initial condition isT(x, y, z, 0) = Tcast, in which Tcast is the casting temperature. As illustrated in [34, 35], the boundary conditions are chosen as follows:(13)−λ∂T∂n=βT−Te,where n is the normal direction of the surface, T and Te are the absolute temperature of the surface and environment, and β is the heat emission coefficient.For humidity analysis,h(x, y, z, 0) = 100% is chosen as the initial condition. For surfaces confronted with air, the boundary condition is shown as follows:(14)−D∂h∂n=Dsurfh−he,where D and Dsurf are coped by equations (9) and (12).


### 2.1.4. Mechanical Submodel

The development of the short-term mechanical properties of early-age concrete, such as tensile strength, compressive strength, and Young’s modulus, is strongly related to the degree of hydration, which depends on both curing age and temperature [36]. The mechanical analyses are conducted with respect to the maturity concept, by adopting the inputs from the thermo-hygro submodels, i.e., the equivalent age te and the potential thermal and shrinkage strains.The total strain of the concrete structure at timeτ is expressed as(15)ετ=εeτ+εcτ+εTτ+εsτ,where εeτ is the linear elastic strain caused by stress, εcτ is the creep strain depending on stress and concrete age at loading, εTτ is the deformation due to temperature fluctuation, and εsτ is the deformation due to shrinkage including drying shrinkage and autogenous shrinkage.The basic creep of concrete is accounted for by the DPL, which has reasonably good performance for either early-age or long-term spans [6]. Note that the DPL formulation is expanded into the Taylor series for discretized analysis. Detailed derivations are found in [37]:(16)εcτ=∫0tJt,τdστ,where(17)Jt,τ=1Eτ+qEττ−dt−τp,where t is the curing age, τ is the concrete age at loading, E−1τ is the modulus of elasticity at loading age, t−τ is the continuous time under load, and p, d, and q are material parameters calibrated by experimental creep tests. According to [37], q = 0.6, d = 0.2, and p = 0.2.The increments of thermal and shrinkage strains are evaluated by the following equations:(18)ΔεTτ=αTΔT,(19)Δεsτ=ksΔh.As addressed for in [38], ks=−ε0sEt0/Ete. ε0s is the final shrinkage strain, and E(t0) and E(te) are the elastic modulus at the initial time and the equivalent age te.To account for the evolution of concrete elastic modulus, tensile strength, and compressive strength during early cement hydration, the following expressions are adopted [6]:(20)Ete=E28⋅es1−672/te−t0nE,ftte=ft28⋅es1−672/te−t0nt,fcte=fc28⋅es1−672/te−t0,where E28, ft28, and fc28 are the elastic modulus, tensile strength, and compressive strength at the age of 28 days, respectively, and 27.45 GPa, 1.78 MPa, and 17.83 MPa are used for E28, ft28, and fc28 (obtained from Figure 2). t0 = 10 hrs is the time when the strength and stiffness are defined to be zero; s, nt, and nE are parameters that account for the cement type used and are fitted by experimental data. 0.211, 0. 59, and 0.37 are selected here, respectively.Figure 2
Development of concrete strength with time.


## 2.1.1. Thermal Submodel

As also employed in [27, 28], the Arrhenius law is found applicable to concrete hydration reaction for generation rate of heat against temperature, which is written as the following equation:(1)kT=A⋅e−EaT/RT,where A is the pre-exponential factor, Ea is the activation energy related to temperature, T is the absolute temperature, and R is the Boltzmann’s constant.It is suggested in [8] that for concrete, Ea = 33500 J/mol for T ≥ 293 K and Ea = 33500 + 1470 × (293-T) J/mol for T < 293 K. According to [28], in this paper, an approximate value for Ea/R of 4000 K is used.By adopting the Arrhenius law, Azenha et al. proposed to determine the maturity degree of concrete by its temperature and curing time [27]. As expressed in equation (2), the maturity degree, characterized as the equivalent age te, is a function of realistic age and temperature in history:(2)te=∫0teEa/R1/T0−1/Tdt,where T0 is the reference temperature (generally T0 = 293 K) and t is the actual curing time. The form of integral is discretized into accumulation of small increments for finite element analysis.Given that the adiabatic temperature rise (θ) is linearly correlated to generated hydration heat (Q), the hydration degree is defined in this paper as follows [29]:(3)αte=QteQmax=θteθmax,where θ(te) is determined empirically by the following equation [4]:(4)θte=Tmax×1−0.57e−0.28te+0.43e−0.004te,where Tmax = 39°C is adopted from Figure 1.Figure 1
Adiabatic temperature rise and autogenous shrinkage strain development with time.According to the work of de Schutter, the relationship between thermal conductivity and hydration degree could be expressed linearly [30]. And equation (5) is adopted for evaluation of the current thermal conductivity:(5)λα=λmax×1.33−0.33α,The specific heat (c) of early-age concrete is proposed by van Breugel as the following equation [31]:(6)c=Wcαccef+Wc1−αcc+Waca+Wwcwρ,where c and W are the specific heat and mass with subscripts c, a, and w representing cement, aggregates, and water. ρ is the density of concrete, and ccef is the assumed specific heat by ccef=0.0084Tc+0.339.By substituting the equivalent age-dependent parameters of equations (4)–(6) into the heat conduction differential equation, the equilibrium of heat within unit volume of concrete is retained as equation (7), also widely employed in [14, 17, 19]:(7)ρcα∂T∂t=∇λα∇T+ρcα∂θ∂t,where c, λ, α, and θ are functions of te as mentioned before, ▽ is the Laplace operator, and t is the curing age of concrete.


## 2.1.2. Moisture Submodel

With regard to concrete structures that are exposed to environmental air at early age, the moisture content in concrete decreases due to the diffusion of moisture. According to Fick’s second law, the rate of moisture transfer per unit area in a certain direction is proportional to the gradient of the moisture concentration in that direction [32]. In analogy to equation (7), by taking advantage of the mass conservation differential equation, the moisture diffusion is controlled by equation (8), also addressed for in [17, 19]:(8)∂h∂t=∂∂xDxh∂h∂x+∂∂yDyh∂h∂y+∂∂zDzh∂h∂z+∂QT∂t+∂Qh∂t+∂QC∂t,where t is the time, h is the relative humidity in concrete, and Dx(h), Dy(h), and Dz(h) are moisture diffusion coefficients in the direction of the x, y, and z-axes, respectively. The subscripts T, h, and C denote diffusion flow induced by temperature gradient, hydration reaction, and cement carbonation. With the consideration that early-age concrete is often composed of high moisture degree in a relatively stable environmental at a constant temperature. The effects of temperature gradient and hydration reaction on moisture distribution are overlooked in this paper. Also given that significant cement carbonization occurs after quite a long time, the last term is omitted. Thus, above all, it is assumed in this paper that the moisture distribution is primarily described by the moisture diffusion coefficients and that the moisture diffusion coefficients are uniformly identical to D(h).By adopting the S-shaped form by Bažant and Najjar and in tandem with the principle of equivalent agete, the moisture diffusion coefficient D(h) is defined by the following equation [33]:(9)Dh=D1⋅α+1−α1+1−h/1−hcn⋅fte,where α=D0/D1, D1 is the maximum of D(h) for h = 100%, D0 is the minimum of D(h) for h = 0%, hc is the pore relative humidity, denoting the humidity drops half-way between D1 and D0, and f(te) is the function of equivalent age. In this paper, α = 0.05, hc = 0.8, and n = 15 are employed. And D1 is estimated from equation (7):(10)D1=D1.0fck/fck0,where D1.0 = 3.6 × 10−6 m2/h, fck0 = 10.0 MPa, and the characteristic compressive strength fck is estimated by mean compressive strength fcm through fck = fcm−8.0 MPa. In analogy to the definition of te, f(te) is characterized as follows:(11)fte=expEadR1T0−1T,where Ead = 35 kJ/mol is the activation energy for moisture diffusion.For boundary conditions to be covered below, the modified Menzel’s equation gives the surface moisture diffusion coefficient (Dsurf) which is as follows [3]:(12)Dsurf=A0.253+0.06Vah−he,where A is an empirical coefficient, determined by the w/c ratios, Va is the average wind speed, and h and he are the relative humidity of the surface and environment, respectively.


## 2.1.3. Initial and Boundary Conditions

For temperature analysis, the initial condition isT(x, y, z, 0) = Tcast, in which Tcast is the casting temperature. As illustrated in [34, 35], the boundary conditions are chosen as follows:(13)−λ∂T∂n=βT−Te,where n is the normal direction of the surface, T and Te are the absolute temperature of the surface and environment, and β is the heat emission coefficient.For humidity analysis,h(x, y, z, 0) = 100% is chosen as the initial condition. For surfaces confronted with air, the boundary condition is shown as follows:(14)−D∂h∂n=Dsurfh−he,where D and Dsurf are coped by equations (9) and (12).


## 2.1.4. Mechanical Submodel

The development of the short-term mechanical properties of early-age concrete, such as tensile strength, compressive strength, and Young’s modulus, is strongly related to the degree of hydration, which depends on both curing age and temperature [36]. The mechanical analyses are conducted with respect to the maturity concept, by adopting the inputs from the thermo-hygro submodels, i.e., the equivalent age te and the potential thermal and shrinkage strains.The total strain of the concrete structure at timeτ is expressed as(15)ετ=εeτ+εcτ+εTτ+εsτ,where εeτ is the linear elastic strain caused by stress, εcτ is the creep strain depending on stress and concrete age at loading, εTτ is the deformation due to temperature fluctuation, and εsτ is the deformation due to shrinkage including drying shrinkage and autogenous shrinkage.The basic creep of concrete is accounted for by the DPL, which has reasonably good performance for either early-age or long-term spans [6]. Note that the DPL formulation is expanded into the Taylor series for discretized analysis. Detailed derivations are found in [37]:(16)εcτ=∫0tJt,τdστ,where(17)Jt,τ=1Eτ+qEττ−dt−τp,where t is the curing age, τ is the concrete age at loading, E−1τ is the modulus of elasticity at loading age, t−τ is the continuous time under load, and p, d, and q are material parameters calibrated by experimental creep tests. According to [37], q = 0.6, d = 0.2, and p = 0.2.The increments of thermal and shrinkage strains are evaluated by the following equations:(18)ΔεTτ=αTΔT,(19)Δεsτ=ksΔh.As addressed for in [38], ks=−ε0sEt0/Ete. ε0s is the final shrinkage strain, and E(t0) and E(te) are the elastic modulus at the initial time and the equivalent age te.To account for the evolution of concrete elastic modulus, tensile strength, and compressive strength during early cement hydration, the following expressions are adopted [6]:(20)Ete=E28⋅es1−672/te−t0nE,ftte=ft28⋅es1−672/te−t0nt,fcte=fc28⋅es1−672/te−t0,where E28, ft28, and fc28 are the elastic modulus, tensile strength, and compressive strength at the age of 28 days, respectively, and 27.45 GPa, 1.78 MPa, and 17.83 MPa are used for E28, ft28, and fc28 (obtained from Figure 2). t0 = 10 hrs is the time when the strength and stiffness are defined to be zero; s, nt, and nE are parameters that account for the cement type used and are fitted by experimental data. 0.211, 0. 59, and 0.37 are selected here, respectively.Figure 2
Development of concrete strength with time.


## 2.2. Extended Finite Element Method

The formulation of the standard displacement function is based on the finite element method and is enriched by the Heaviside function and the crack tip enrichment functions [23]:(21)ux=∑i∈NNixui+∑j∈NdisNjHxai+∑k∈NasyNk∑a=14Φaxbka,where j is the set of nodes whose shape function support is cut by a crack, i is the set of all nodes of the domain, Ni is the shape function associated with the node i, K is the set of nodes whose shape function support contains the crack front, and ui, ai, and bka are the displacements associated with the corresponding nodes.The Heaviside jump function,Hx, is a discontinuous function across the crack surface and is constant on each side of the crack [26]:(22)Hx+1above the crack−1beneath the crack.The crack tip enrichment functions,Φax, are obtained from the asymptotic displacement fields [22, 23]:(23)Φar,θ=rsinθ2,rcosθ2,rsinθsinθ2,rsinθcosθ2.The cohesive property of the crack tip is considered in this paper, and that means the stresses at the crack tip are not singular. [39, 40] carried out an asymptotic analysis for very large structure of the mechanical fields in a cohesive zone. It is better to consider the following functions [41]:(24)Φar,θ=rsinθ2orr3/2sinθ2orr2sinθ2.Constitutive formulation about the cohesive crack model can be seen in [24, 41, 42] while formulation information of XFEM can be found in [22, 23, 43].Gauss integral is employed for the numerical integration. For simplification, the stresses of Gauss integral points are used to control the crack propagation. Elements fail when maximum principal stresses of the Gauss integral points exceed the tensile strength.


## 2.3. Submodel Coupling

Fundamentally, concrete structures are under complex states of stress after being placed at construction sites due to the stress-inducing mechanisms of temperature effects, moisture effects, mechanical effects, and their combinations. For simplification, several assumptions shall be assumed between submodels. Further information regarding coupling methods is addressed for in this section.As shown in Figure3, the temperature, moisture, and mechanical fields are connected in a certain order via one-way coupling. Temperature analysis is conducted first, in which the hydration procedure for HVFA concrete is simulated according to equation (4); both element temperature and the corresponding hydration degree are obtained.Figure 3
Flowchart of the modeling process.Even though moisture diffusion affects temperature analysis to some extent, it is small enough to be omitted. Moisture diffusion properties are, however, strongly dependent on temperature. Accordingly, temperature analysis is initially conducted before moisture computation. Then, the thermal results are used as inputs for calculation of the moisture distribution under temperature effects.With respect to mechanical fields, temperature evolution will not only lead to volumetric changes but also influence the development of the concrete properties based on the equivalent age. Moreover, drying shrinkage induced by obvious surface moisture gradients via moisture transport to the environment will lead to shrinkage stress under restraint boundary conditions, thereby significantly affecting the stress responses in the mechanical field.The mechanical field also influences the temperature and the moisture fields upon two reasons: (1) the inner transportation of temperature and moisture are subsequently affected by the existence and propagation of cracks in the concrete structure and (2) the variations of concrete volume will slightly influence the boundary conditions of the temperature and moisture fields. Though it may lead to inaccuracy, it is assumed to be omitted the effects of cracks on temperature and moisture fields, which is characterized in Figure3 as a one-way coupling method. For more realistic results, a two-way coupling method could be expected to be implemented afterwards.Consequently, the mechanical field is the last step in the one-way coupling method, where the results of the previous two fields (temperature and moisture) are used as inputs to analyze the stress development and corresponding crack propagation. In this manner, relatively accurate results of early-age concrete behavior are obtained. By use of ABAQUS user subroutines, the specific implementation is as follows:Temperature Analysis. User subroutine UMATHT is firstly called by use of adiabatic temperature rise equation and other temperature parameters, such as thermal conductivity and specific heat. At the meantime, the corresponding initial temperature is imported together with reasonable boundary conditions (ambient temperature and heat transfer coefficient) in terms of user subroutine FILM. The obtained results are saved and named in text documents in the order of the corresponding element numbers at the end of the computation.Moisture Analysis. Similar user subroutines are used herein for the humidity distribution. UMATHT is called to input controlling equations to simulate variations of the moisture diffusion coefficient, and FILM is used to define the surface moisture diffusion coefficient and ambient relative humidity. Since moisture diffusion coefficients are dependent on temperature effects, a self-defined procedure is adopted to read the temperature data, and then the texts are saved in the form of vectors for moisture distribution calculation. In the end, the results are collected in the same manner as the temperature analysis.Mechanical-Cracking Analysis. A user subroutine UMAT was coded for mechanical analysis. As expressed by equation (15), the changing distribution of temperature and humidity is translated as additional strains, in tandem with the creep strain, into the stress field after the former coupling analysis. Thermal expansion coefficient is adopted for coupling of the temperature-mechanical field, while taking into consideration the development of mechanical properties, such as the modulus of elasticity, compressive strength, and tensile strength, over the equivalent age. As indicated by equation (19), drying shrinkage strains are gained from moisture distribution, which is superposed onto the stress field together with the thermal strains. Hence, the thermal, moisture, and creep effects are considered altogether, and the corresponding stress is computed.The XFEM method employed in this paper has been built in ABAQUS with a user-defined cohesive law. Prior to cracking, the initial response is controlled by the preceding UMAT. The damage criterion and evolution are, respectively, determined by the current tensile strength as equations (20) and (25) and fracture energy controlled traction-separation response as equation (26).(25)f=σmaxftte,where  is the Macaulay bracket and σmax represents the maximum principal stress. Damage is assumed to initiate when the value of the principal stress ratio, i.e., f, reaches a value of one.(26)tn=1−DTn,Tn≥0,Tn,otherwise,tJ=1−DTJ,J=s,t,where D represents the averaged overall damage at the intersection between the crack surfaces and the edges of cracked elements, evolving initially from 0 to 1. D is calculated by ABAQUS in the form of an exponential softening procedure with the Power Law.In practice, the UEL subroutine is adopted to combine XFEM with the DPL model to simulate crack initiation, evolution, and the corresponding stress development. The fracture energy was empirically estimated at 40 N/m for the normal mode and 10 N/m for the shear mode, in line with the values suggested by the ABAQUS documentation. And the contact properties between crack surfaces are specified as Penalty for tangential behavior and “Hard” Contact for normal behavior.


## 3. Numerical Simulation

### 3.1. Material

This section introduces the results of a contrast test to show the differences between HVFA concrete and ordinary Portland cement (OPC) in heat dissipation, autogenous shrinkage, and strength properties. Some input data for numerical simulation are obtained from this experiment. The replacement proportion level of fly ash is 60% by weight for the paste at the water-to-binder (w/b) ratio of 0.41 in the present paper. Details of the mix proportions of the two pastes are shown in Table1, where ① represents HVFA concrete, and ② represents OPC.Table 1
Concrete-mixed proportion.
w/b
Sand/aggregates
Weight per unit volume (kg/m3)
Water
Cement
Fly ash
Sand
Gravel
0.41
0.41
① 164
160
240
736
1060
② 164
400
0
736
1060As shown in Figure2, the thermal and mechanical properties were gained by experiments for data fitting and verification. The replacement of cement with fly ash significantly reduced the early-age compressive strength, tensile strength, and elastic modulus. It could be deduced that the relatively lower early strength increases the likelihood of early-age cracking.Adiabatic temperature rise and autogenous shrinkage history for the proposed OPC and HVFA concrete are shown in Figure1. Clearly, the influence of hydration heat and autogenous shrinkage are efficiently reduced by fly ash. Autogenous shrinkage is added into the total strain according to equation (15).


### 3.2. Outline of the Analysis

The YD sluice is situated in the northeast region of Tianjin located on the unconsolidated quaternary deposits of the alluvial plain. The lock floor is a monolithic raft-plain structure, measuring 35 meters width, with twenty lock chambers in total, two of which are connected together as one independent part. 60% HVFA concrete was used for pier construction, and the mixture proportion is listed in Table1. According to the design, piers were cast five days after the floor had been placed at the construction site, and the casting for each pier lasted for 2 days during November 2008. The formwork wall for casting consists of 5 cm thick wooden plates which were removed after 28 days. After removing the formwork wall, cracks were found in several piers as illustrated in Figure 4.Figure 4
Shapes and positions of the objective cracks.For simulation of the crack propagation, a pier that underwent the lowest temperature in history and held the most representative cracks was evaluated herein. As shown in Figure4, a typical pattern of early-age cracking was found after the formwork wall had been removed. Cracks embedded at the pier appeared at the distance of 3.2 m, 6 m, and 9 m from its round cut point and stopped their propagation at the length of 1 m, 0.9 m, and 1 m; the width is around 0.2 mm.For the numerical analysis presented here, a massive concrete pier was modelled, with dimensions of 12.4 m in height along thez-axis, 28 m length along the y-axis, and 2 m width along the x-axis. Additionally, the pier was placed on a concrete foundation base with dimensions of 35 m length along the x-axis, 36 m in width along the y-axis, and 2–3.5 m varying thickness along the z-axis for simulation of heat transfer and constraint of the concrete structures. The shape and configurations of the structure are shown in Figure 5. For simplification, the formwork wall was not modelled in this numerical simulation, but its influence was considered in the form of temperature boundary conditions. The casting temperature for the pier was 15°C at the beginning of the calculation, and the foundation was naturally cooled to zero. Concrete structures were meshed with eight-node isoperimetric solid elements. Input parameters used in the analysis are shown in Table 2. The effects of night cooling and solar radiation were properly neglected in this study.Figure 5
Structural configuration and data acquisition points.Table 2
Input data in numerical analysis.
Parameters
Pier
Foundation
Adiabatic temperature rise
Maximum temperature rise (°C)
39
—
Thermal properties
Thermal conductivity [kJ/(m·d·°C)]
10
10
Specific heat [kJ/(kg·°C)]
0.87
0.87
Unit density (kg/m3)
2360
2360
Thermal expansion coefficient (/°C)
0.9 × 10−6
1 × 10−5
Convective heat transfer coefficient [W/(m2·°C)]
45
45
Moisture properties
Surface moisture diffusion coefficient (m2/h)
2.1 × 10−4
2.1 × 10−4
Moisture diffusion coefficienta (m2/h)
3.6 × 10−6
3.6 × 10−6
Mechanical properties
Final setting time (h)
10
10
Poisson’s ratio
0.167
0.167
Tensile strengthb (MPa)
1.78
1.78
Modulus of elasticityb (GPa)
20
20
Comprehensive strengthb (MPa)
24.5
24.5
Fracture energy (normal mode) (N/m)
40
−
Fracture energy (shear mode) (N/m)
10
−
Environmental conditions
Temperaturec (°C)
−1
−1
Humidity (%)
60
60
aMaximum value of moisture diffusion coefficient according to equation (10); bstrength of 28 days; cinfluenced by the existence of the formwork wall.The presence of autogenous shrinkage leads to extracompression force imposed on the reinforcement; therefore, reducing its effect on the control of early-age cracking [8]. Moreover, the reinforcement had little impact on stress distribution in the computation. Considering the similarity between the present analysis and the one conducted in [8], it is reasonable to neglect the reinforcement influence in this paper.To compare the variation of temperature and moisture in different regions, three representative output points were separately selected, as shown in Figure5: the pier core (point A), the core of the bottom (point B), and the lateral surface (point C).


### 3.3. Temperature Analysis

Figure6 reveals the temperature distribution 100 hours and 672 hours after casting; the corresponding typical temperature history curves around the three reference points are depicted in Figure 7. As a result of cement hydration, a significant temperature rise was observed during the early stage. However, because of constant heat dissipation into the environment, the concrete underwent nonuniform cooling, which occurred faster on the surface and slower in the core. Thus, a significant temperature gradient between the core and the surface was captured (Figure 6(a)). As the analysis proceeded, the temperatures in both the core and the surface continuously decreased. At the final step, the computed temperature reached a relatively uniform state and was close to the surrounding temperature (−1°C).Figure 6
Temperature distributions at 100 hours (a) and 672 hours (b) after casting.
(a)
(b)Figure 7
Temperature development in OPC and HVFA concrete.Figure7 shows the numerical temperature history curves of three marked points in HVFA concrete and a core point in OPC (Table 1). Obvious temperature differences were observed between the OPC point and point A. This phenomenon led to a lower temperature gradient between the core and the surface in HVFA concrete compared to OPC if consistent environmental temperature conditions were considered. In other words, the stress effects of the HVFA concrete caused by thermal contraction in conjunction with structural restraint are weakened, and thus, the cracking potential of the concrete is reduced compared to OPC.For 60% HVFA concrete, as observed in Figure7, the temperature development trends of the inner points (A, B) were almost consistent, i.e., a rapid temperature rise during the initial period followed by a slow decrease. For the surface point (C), a steady decrease followed by a relatively rapid one was captured. The inner portion reached its maximum temperature of 35.03°C after 14 hours. Within the framework of the temperature submodel, the degree of hydration at point A was the largest among the three points. The temperature at point A underwent the highest temperature rise history, whereas the temperature at point C underwent the lowest. The figure clearly shows that the surface region of the pier suffered more temperature losses and greater temperature gradient emerged between the core and the surface.


### 3.4. Moisture Analysis

Figure8 presents the analysis results for the relative humidity distribution within 28 days. The moisture variation was going at a much slower rate than the temperature was, with only 0.043 descent for humidity observed on the surface after the concrete had been cast for 100 hours (Figure 8(a)). After 672 hours, an obvious moisture gradient appeared, and the relative humidity in the core remained unchanged. In the bottom of the pier, the section that is in contact with the foundation (where humidity is transferred more directly than other places) exhibited the earliest apparent humidity gradient. At the end of the analysis (Figure 8(b)), 40 cm beneath the concrete surface showed the greatest moisture gradient. The lowest relative moisture (which was close to atmospheric humidity) was present at the bottom of the pier and was approximately 0.62. The shrinkage stress caused by the loss of moisture strongly affects the surface. Thus, for young concrete, it is essential to attach equal importance to both moisture and temperature curing to prevent large moisture gradients.Figure 8
Relative humidity distributions at 100 hours (a) and 672 hours (b) after casting.
(a)
(b)Figure9 shows the relative humidity trend of the corresponding reference points. A point in the center maintained a constant humidity, whereas the other points suffered monotonic moisture losses. The curve of point B decreased quickly first and then turned slowly when they were close to the ambient humidity. By comparing the computed relative humidity among all three points, the bottom of the pier clearly suffered the most severe humidity losses and therefore was most sensitive to shrinkage stress and subsequent failures.Figure 9
Relative humidity development with time.


### 3.5. Evolution of Cracking

Taking into consideration temperature strain, shrinkage strain (including drying shrinkage and autogenous shrinkage), and creep, the maximum principal strain evolution along with the corresponding crack propagation are illustrated in Figure10.Figure 10
Cracking extension and the corresponding mechanical strain development. (a) 300 hours after casting. (b) 672 hours after casting.
(a)
(b)The figure shows that strain near the three cracks is much larger than elements at other locations. As Figure10(a) shows, 300 hours after casting, the maximum principal strain appeared on the bottom surface around the tips of cracks ② and ③, whereas a rather uniform strain distribution was observed around the tip of crack ①. These differences in strain distribution implied that cracks ② and ③ still had the potential to propagate. At the end of the calculation (672 hours), as shown in Figure 10(b), cracks ② and ③ continued to propagate, while the maximum strain appeared at the initial position of crack ①, and that indicates crack ① remained unchanged afterwards.Figure11 reveals the maximum principal stress evolution of the six reference points marked in Figure 10(b) along with the tensile strength of 60% HVFA concrete used in this paper. Due to the cohesive zone existing around the crack tip, stresses were no longer singular. Cracks propagated when the maximum principle stresses of Gauss integral points exceeded the corresponding tensile strength. Tensile stress produced at the concrete surface during the early-age stage was caused by the rapid temperature variation inside the concrete structure and the corresponding restraint effect. The figure shows that the pier experienced a rapid early stress development, and 100 hours after casting, points 1, 2, and 3 attained the tensile stress of 0.80 MPa. At this point, the tensile strength of the pier was still at a rather low level; however, the achieved principal stress of 0.80 MPa had been large enough to form a crack. The same phenomenon also occurred at points 4, 5, and 6. Despite the relatively slower stress development and the simultaneously developing concrete strength, the stress development of the concrete still exceeded its tensile strength and resulted in new crack occurrence.Figure 11
Tensile strength evolution against time.During the cracking procedure, crack ① in Figure10 led throughout the experiment and stopped propagation first, whereas cracks ② and ③ achieved their final form slightly later. This particular cracking patterns occurred because the effect of shrinkage in the downstream pier round cut was larger than in the middle part and formed a gradual decreasing trend of tensile stress from the downstream part to the middle of the pier.Though the addition of fly ash reduced the temperature differences between the core and the surface to some degree and slowed down the consequent stress development, the relatively lower early stress was still large enough to cause cracking because of the low early strength of HVFA concrete. Additionally, the great humidity gradient that emerged from the rapid moisture losses at the surface (or the bottom) after exposure also contributed to the higher surface (or the bottom) tensile stress development and the consequent surface cracking.According to the numerical analysis, typical patterns of early-age cracking for a HVFA concrete pier caused by thermal, hydric, and mechanical conditions are shown in Figure10(b). Three computed cracks appeared at distances of 2.9 m, 5.7 m, and 8.5 m from the round cut. Though, due to other external factors, such as the chosen mesh element and inadequate environmental and construction conditions, as well as simplified governing equation of moisture, the numerical results did not exactly match the objective results, and the distribution trends of temperature, humidity, and cracks were well captured.


## 3.1. Material

This section introduces the results of a contrast test to show the differences between HVFA concrete and ordinary Portland cement (OPC) in heat dissipation, autogenous shrinkage, and strength properties. Some input data for numerical simulation are obtained from this experiment. The replacement proportion level of fly ash is 60% by weight for the paste at the water-to-binder (w/b) ratio of 0.41 in the present paper. Details of the mix proportions of the two pastes are shown in Table1, where ① represents HVFA concrete, and ② represents OPC.Table 1
Concrete-mixed proportion.
w/b
Sand/aggregates
Weight per unit volume (kg/m3)
Water
Cement
Fly ash
Sand
Gravel
0.41
0.41
① 164
160
240
736
1060
② 164
400
0
736
1060As shown in Figure2, the thermal and mechanical properties were gained by experiments for data fitting and verification. The replacement of cement with fly ash significantly reduced the early-age compressive strength, tensile strength, and elastic modulus. It could be deduced that the relatively lower early strength increases the likelihood of early-age cracking.Adiabatic temperature rise and autogenous shrinkage history for the proposed OPC and HVFA concrete are shown in Figure1. Clearly, the influence of hydration heat and autogenous shrinkage are efficiently reduced by fly ash. Autogenous shrinkage is added into the total strain according to equation (15).


## 3.2. Outline of the Analysis

The YD sluice is situated in the northeast region of Tianjin located on the unconsolidated quaternary deposits of the alluvial plain. The lock floor is a monolithic raft-plain structure, measuring 35 meters width, with twenty lock chambers in total, two of which are connected together as one independent part. 60% HVFA concrete was used for pier construction, and the mixture proportion is listed in Table1. According to the design, piers were cast five days after the floor had been placed at the construction site, and the casting for each pier lasted for 2 days during November 2008. The formwork wall for casting consists of 5 cm thick wooden plates which were removed after 28 days. After removing the formwork wall, cracks were found in several piers as illustrated in Figure 4.Figure 4
Shapes and positions of the objective cracks.For simulation of the crack propagation, a pier that underwent the lowest temperature in history and held the most representative cracks was evaluated herein. As shown in Figure4, a typical pattern of early-age cracking was found after the formwork wall had been removed. Cracks embedded at the pier appeared at the distance of 3.2 m, 6 m, and 9 m from its round cut point and stopped their propagation at the length of 1 m, 0.9 m, and 1 m; the width is around 0.2 mm.For the numerical analysis presented here, a massive concrete pier was modelled, with dimensions of 12.4 m in height along thez-axis, 28 m length along the y-axis, and 2 m width along the x-axis. Additionally, the pier was placed on a concrete foundation base with dimensions of 35 m length along the x-axis, 36 m in width along the y-axis, and 2–3.5 m varying thickness along the z-axis for simulation of heat transfer and constraint of the concrete structures. The shape and configurations of the structure are shown in Figure 5. For simplification, the formwork wall was not modelled in this numerical simulation, but its influence was considered in the form of temperature boundary conditions. The casting temperature for the pier was 15°C at the beginning of the calculation, and the foundation was naturally cooled to zero. Concrete structures were meshed with eight-node isoperimetric solid elements. Input parameters used in the analysis are shown in Table 2. The effects of night cooling and solar radiation were properly neglected in this study.Figure 5
Structural configuration and data acquisition points.Table 2
Input data in numerical analysis.
Parameters
Pier
Foundation
Adiabatic temperature rise
Maximum temperature rise (°C)
39
—
Thermal properties
Thermal conductivity [kJ/(m·d·°C)]
10
10
Specific heat [kJ/(kg·°C)]
0.87
0.87
Unit density (kg/m3)
2360
2360
Thermal expansion coefficient (/°C)
0.9 × 10−6
1 × 10−5
Convective heat transfer coefficient [W/(m2·°C)]
45
45
Moisture properties
Surface moisture diffusion coefficient (m2/h)
2.1 × 10−4
2.1 × 10−4
Moisture diffusion coefficienta (m2/h)
3.6 × 10−6
3.6 × 10−6
Mechanical properties
Final setting time (h)
10
10
Poisson’s ratio
0.167
0.167
Tensile strengthb (MPa)
1.78
1.78
Modulus of elasticityb (GPa)
20
20
Comprehensive strengthb (MPa)
24.5
24.5
Fracture energy (normal mode) (N/m)
40
−
Fracture energy (shear mode) (N/m)
10
−
Environmental conditions
Temperaturec (°C)
−1
−1
Humidity (%)
60
60
aMaximum value of moisture diffusion coefficient according to equation (10); bstrength of 28 days; cinfluenced by the existence of the formwork wall.The presence of autogenous shrinkage leads to extracompression force imposed on the reinforcement; therefore, reducing its effect on the control of early-age cracking [8]. Moreover, the reinforcement had little impact on stress distribution in the computation. Considering the similarity between the present analysis and the one conducted in [8], it is reasonable to neglect the reinforcement influence in this paper.To compare the variation of temperature and moisture in different regions, three representative output points were separately selected, as shown in Figure5: the pier core (point A), the core of the bottom (point B), and the lateral surface (point C).


## 3.3. Temperature Analysis

Figure6 reveals the temperature distribution 100 hours and 672 hours after casting; the corresponding typical temperature history curves around the three reference points are depicted in Figure 7. As a result of cement hydration, a significant temperature rise was observed during the early stage. However, because of constant heat dissipation into the environment, the concrete underwent nonuniform cooling, which occurred faster on the surface and slower in the core. Thus, a significant temperature gradient between the core and the surface was captured (Figure 6(a)). As the analysis proceeded, the temperatures in both the core and the surface continuously decreased. At the final step, the computed temperature reached a relatively uniform state and was close to the surrounding temperature (−1°C).Figure 6
Temperature distributions at 100 hours (a) and 672 hours (b) after casting.
(a)
(b)Figure 7
Temperature development in OPC and HVFA concrete.Figure7 shows the numerical temperature history curves of three marked points in HVFA concrete and a core point in OPC (Table 1). Obvious temperature differences were observed between the OPC point and point A. This phenomenon led to a lower temperature gradient between the core and the surface in HVFA concrete compared to OPC if consistent environmental temperature conditions were considered. In other words, the stress effects of the HVFA concrete caused by thermal contraction in conjunction with structural restraint are weakened, and thus, the cracking potential of the concrete is reduced compared to OPC.For 60% HVFA concrete, as observed in Figure7, the temperature development trends of the inner points (A, B) were almost consistent, i.e., a rapid temperature rise during the initial period followed by a slow decrease. For the surface point (C), a steady decrease followed by a relatively rapid one was captured. The inner portion reached its maximum temperature of 35.03°C after 14 hours. Within the framework of the temperature submodel, the degree of hydration at point A was the largest among the three points. The temperature at point A underwent the highest temperature rise history, whereas the temperature at point C underwent the lowest. The figure clearly shows that the surface region of the pier suffered more temperature losses and greater temperature gradient emerged between the core and the surface.


## 3.4. Moisture Analysis

Figure8 presents the analysis results for the relative humidity distribution within 28 days. The moisture variation was going at a much slower rate than the temperature was, with only 0.043 descent for humidity observed on the surface after the concrete had been cast for 100 hours (Figure 8(a)). After 672 hours, an obvious moisture gradient appeared, and the relative humidity in the core remained unchanged. In the bottom of the pier, the section that is in contact with the foundation (where humidity is transferred more directly than other places) exhibited the earliest apparent humidity gradient. At the end of the analysis (Figure 8(b)), 40 cm beneath the concrete surface showed the greatest moisture gradient. The lowest relative moisture (which was close to atmospheric humidity) was present at the bottom of the pier and was approximately 0.62. The shrinkage stress caused by the loss of moisture strongly affects the surface. Thus, for young concrete, it is essential to attach equal importance to both moisture and temperature curing to prevent large moisture gradients.Figure 8
Relative humidity distributions at 100 hours (a) and 672 hours (b) after casting.
(a)
(b)Figure9 shows the relative humidity trend of the corresponding reference points. A point in the center maintained a constant humidity, whereas the other points suffered monotonic moisture losses. The curve of point B decreased quickly first and then turned slowly when they were close to the ambient humidity. By comparing the computed relative humidity among all three points, the bottom of the pier clearly suffered the most severe humidity losses and therefore was most sensitive to shrinkage stress and subsequent failures.Figure 9
Relative humidity development with time.


## 3.5. Evolution of Cracking

Taking into consideration temperature strain, shrinkage strain (including drying shrinkage and autogenous shrinkage), and creep, the maximum principal strain evolution along with the corresponding crack propagation are illustrated in Figure10.Figure 10
Cracking extension and the corresponding mechanical strain development. (a) 300 hours after casting. (b) 672 hours after casting.
(a)
(b)The figure shows that strain near the three cracks is much larger than elements at other locations. As Figure10(a) shows, 300 hours after casting, the maximum principal strain appeared on the bottom surface around the tips of cracks ② and ③, whereas a rather uniform strain distribution was observed around the tip of crack ①. These differences in strain distribution implied that cracks ② and ③ still had the potential to propagate. At the end of the calculation (672 hours), as shown in Figure 10(b), cracks ② and ③ continued to propagate, while the maximum strain appeared at the initial position of crack ①, and that indicates crack ① remained unchanged afterwards.Figure11 reveals the maximum principal stress evolution of the six reference points marked in Figure 10(b) along with the tensile strength of 60% HVFA concrete used in this paper. Due to the cohesive zone existing around the crack tip, stresses were no longer singular. Cracks propagated when the maximum principle stresses of Gauss integral points exceeded the corresponding tensile strength. Tensile stress produced at the concrete surface during the early-age stage was caused by the rapid temperature variation inside the concrete structure and the corresponding restraint effect. The figure shows that the pier experienced a rapid early stress development, and 100 hours after casting, points 1, 2, and 3 attained the tensile stress of 0.80 MPa. At this point, the tensile strength of the pier was still at a rather low level; however, the achieved principal stress of 0.80 MPa had been large enough to form a crack. The same phenomenon also occurred at points 4, 5, and 6. Despite the relatively slower stress development and the simultaneously developing concrete strength, the stress development of the concrete still exceeded its tensile strength and resulted in new crack occurrence.Figure 11
Tensile strength evolution against time.During the cracking procedure, crack ① in Figure10 led throughout the experiment and stopped propagation first, whereas cracks ② and ③ achieved their final form slightly later. This particular cracking patterns occurred because the effect of shrinkage in the downstream pier round cut was larger than in the middle part and formed a gradual decreasing trend of tensile stress from the downstream part to the middle of the pier.Though the addition of fly ash reduced the temperature differences between the core and the surface to some degree and slowed down the consequent stress development, the relatively lower early stress was still large enough to cause cracking because of the low early strength of HVFA concrete. Additionally, the great humidity gradient that emerged from the rapid moisture losses at the surface (or the bottom) after exposure also contributed to the higher surface (or the bottom) tensile stress development and the consequent surface cracking.According to the numerical analysis, typical patterns of early-age cracking for a HVFA concrete pier caused by thermal, hydric, and mechanical conditions are shown in Figure10(b). Three computed cracks appeared at distances of 2.9 m, 5.7 m, and 8.5 m from the round cut. Though, due to other external factors, such as the chosen mesh element and inadequate environmental and construction conditions, as well as simplified governing equation of moisture, the numerical results did not exactly match the objective results, and the distribution trends of temperature, humidity, and cracks were well captured.


## 4. Conclusions

From the numerical analysis for HVFA concrete structures, the following conclusions are obtained with focus on its early-age behavior:(1)
The proposed XFEM-combined thermal-hygro-mechanical model has the ability to predict the crack path during the early stage of HVFA concrete structures.(2)
From the results of the numerical analysis for piers containing HVFA concrete, thermal gradient, surface moisture losses, and restraint are the main causes of early-age cracking initiation and propagation.(3)
The computed cracks are roughly consistent with the observed ones in structures, presenting inclined propagation from the middle part of the foundation to the downstream pier round cut.(4)
Although the addition of fly ash could reasonably reduce the rate of heat dissipation from hydration and temperature gradient between the core and the surface of the concrete structure to some extent, significance should still be attached to crack problems due to the low early strength of fly ash. Therefore, timely curing is very important for young HVFA concrete.Though crack tip enrichment functions are considered in the formulation, it will penetrate the entire element once an element cracks. This problem will be further studied in the following research.



---

*Source: 5793902-2020-01-08.xml*