Hydro-Thermo-Mechanical Analysis of an Existing Gravity Dam Undergoing Alkali–Silica Reaction

Abstract

The alkali–silica reaction is a chemical phenomenon that, by inducing expansion and the formation of cracks in concrete, can have a severe impact on the safety and functioning of existing concrete dams. Starting from a phenomenological two-phase isotropic damage model describing the degradation of concrete, the effects of alkali-silica reaction in an existing concrete gravity dam are evaluated and compared with real monitoring data. Considering the real temperature and humidity variations, the influence of both temperature and humidity are considered through two uncoupled diffusion analyses: a heat diffusion analysis and a moisture diffusion analysis. The numerical analyses performed with the two-phase damage model allow for prediction of the structural behaviour, both in terms of reaction extent and increase of crest displacements. The crest displacements are compared with the real monitoring data, where reasonably good agreement is obtained.

1. Introduction

The alkali–silica reaction (ASR) is a chemical phenomenon that may occur over time in concrete between the alkaline cement paste and the reactive silica present in some kinds of aggregates. This reaction causes expansion of the altered aggregate by the formation of a viscous gel, which increases in volume when absorbing water. This gel exerts an expansive pressure inside the siliceous aggregate and/or at the interface between aggregates and the cement paste, causes spalling and loss of stiffness of the concrete and, finally, leads to the formation of cracks [1]. While ASR can be avoided in new concrete structures by the selection of non-reactive aggregates and/or by the chemical composition of cement (which determines the alkali ion concentration in the pore solution [2]) and additions [3,4,5], ASR can lead to serious cracking in concrete structures built several decades ago in very wet environments (e.g., dams), resulting in critical structural problems. It takes many years for the reaction to develop and the symptoms to become visible, in terms of an increase in displacements and of a pattern of diffuse localized micro- and macro-cracks [6]. Several examples of important concrete structures, such as dams or shielding in nuclear plants, subject to ASR have been reported in the literature (see, e.g., [7,8,9]).

The expansion of concrete as a result of the reaction between cement and aggregates was first studied in [1]. Since then, much research effort has been devoted to describing the reactions occurring between aggregates and cement paste, in order to experimentally characterize the consequent expansion (as in [10]) and to develop material models to describe the mechanical effects at different scales (see, e.g., [11,12]) but, very often, there is a lack of validation of the proposed models at a structural level, as has been pointed out in [13]. Significant examples of analyses conducted on dams affected by ASR can be found in [14,15,16,17].

The purpose of this paper is to further validate the bi-phase damage model proposed in [18] by simulating the experimental campaign reported in [19] and to present a case study concerning the analysis, using that model, of a real concrete dam affected by ASR. As the kinetics of ASR strongly depend on the evolution of the temperature and humidity, as described in [20,21], we use real monitoring data of temperature variations and air humidity to compute the temperature and degree of saturation fields through preliminary heat and moisture diffusion analyses. Then, following a weakly coupled strategy, we perform an elasto-damage analysis, which allows for computation of the displacement evolution and the damage patterns.

2. Bi-Phase Damage Model

In this section, a brief summary of the chemo-damage model is given (for details, see [18]). At the meso-scale, concrete affected by ASR is described as a two-phase heterogeneous material constituted by a solid skeleton s and a wet gel g produced by ASR. The total volume V is, hence, divided into $V_s$, the solid volume, and $V_g$, the volume occupied by the dry gel $dg$, liquid water w, and the vapor and air phase $va$ (see Figure 1).

Assuming small strains and quasi-static conditions, the compatibility and equilibrium equations for the bi-phase solid are

$$\mathit{\epsilon }=\frac{1}{2}(grad\mathbf{u}+gra{d}^T\mathbf{u}),$$

(1)

$$div\mathit{\sigma }+\rho \mathbf{b}=\mathbf{0},$$

(2)

where $\mathit{\epsilon }$ is the small strain tensor of the skeleton; $\mathbf{u}$ is the skeleton displacement; $\mathit{\sigma }$ is the Cauchy stress; $\rho ={\rho }_s(1-\varphi )+{\rho }_g\varphi$ is the mass density (${\rho }_s$ and ${\rho }_g$ being the intrinsic solid and fluid mass densities, respectively); and $\rho \mathbf{b}$ is the body force of the combination of solid and fluid. In this model, the variable which accounts for the gel formation and expansion (considered simultaneously) is the volumetric fraction of the wet gel,

$${\zeta }_g=\frac{V_g}{V}.$$

(3)

To take into account the influence of humidity and temperature, the degree of saturation $S_w=V_w/(V-V_s-V_{dg})$ and the temperature T are included as additional state variables and proper conservation and transport laws are necessary. In particular, by combining the mass conservation law for liquid water and the transport law for moisture in its liquid form, as proposed in [22], and the relation between the water pressure and the degree of saturation given by the experimental capillarity curve (see, e.g., [23]), one obtains the non-linear transport law for moisture (in its liquid form) as

$$div\left(D_{w}grad{S}_w\right)+{\varphi }_w\frac{\partial S_w}{\partial t}=0,$$

(4)

where $D_w$ is the permeability of concrete, dependent on the degree of saturation $S_w$

$$D_w\left(S_w\right)=m_1{m}_2\left(1-\frac{1}{m_2}\right)k\frac{{[1-{(1-S_w^{1/m_3})}^{m_3}]}^2}{\left[{\eta }_w{S}_w^{(1/2+m_2)}{(S_w^{-m_2}-1)}^{1/m_2}\right]},$$

(5)

where $m_1$, $m_2$, and $m_3$ are material parameters; k is the intrinsic permeability; and ${\eta }_w$ is the dynamic viscosity.

Regarding the temperature, the heat conservation law, combined with linear and isotropic Fourier conduction laws, gives the following heat transport law

$$C\frac{\partial T}{\partial t}=div\left(D_{T}gradT\right),$$

(6)

where C is the volumetric heat capacity of concrete and $D_T$ is the isotropic heat conductivity coefficient.

In a bi-phase model, the stress $\mathit{\sigma }$ is the sum of the effective stress ${\mathit{\sigma }}^{\mathbf{\prime }}$ (acting on the skeleton) and of the stress -$bp\mathbf{1}$ (acting on the fluid phase), where b is Biot’s coefficient and p is the fluid pressure (see [24]). In the bi-phase model proposed in [18,25], the fluid phase is constituted of the wet gel and the material behaviour is assumed to be elastic with isotropic damage. The state equations, hence, read

$$\mathit{\sigma }=(1-D)[2G\mathbf{e}+K(tr\mathit{\epsilon }-\alpha \theta \left)\mathbf{1}\right]-bp\mathbf{1},$$

(7)

$$p=-(1-D)M(btr\mathit{\epsilon }-{\zeta }_g-{\alpha }_g\theta ),$$

(8)

where G and K are the shear and bulk moduli of the homogenized concrete skeleton, respectively; $\mathbf{e}$ is the deviatoric strain tensor; M is the Biot modulus; $\alpha$ and ${\alpha }_g$ are the volumetric coefficients of thermal expansion of the skeleton and of the gel, respectively; $\theta =T-T_0$ is the temperature variation, with respect to a reference temperature $T_0$; and D is the isotropic damage variable.

The constitutive model is completed by the definition of the evolution law for the gel volume content, assumed to be proportional to the reaction extent $\xi$, a phenomenological internal variable varying from 0 to 1,

$${\zeta }_g=\frac{K+M{b}^2}{Mb}{\mathit{\epsilon }}_{ASR}^{\infty }\xi ,$$

(9)

where ${ϵ}_{ASR}^{\infty }$ is the free asymptotic volumetric expansion due to ASR in the isothermal fully saturated case ($S_w$ = 1).

The following form for the reaction rate is considered:

$$\dot{\xi }=\frac{\langle f_{Sw}-\xi \rangle }{\tilde{t}},$$

(10)

with

$$f_{Sw}=\frac{1+b_{1}exp(-b_2)}{1+b_{1}exp(-b_2{S}_w)},$$

(11)

where $b_1$ and $b_2$ are material parameters, which are calibrated with experimental data.

In the present work, the following law for the intrinsic time $\tilde{t}$ is assumed, as in [18]

$$\frac{1}{\tilde{t}}=\frac{\frac{\xi }{f_{Sw}}+exp(-\frac{{\tau }_{lat}}{{\tau }_{ch}})}{{\tau }_{ch}(1+exp(-\frac{{\tau }_{lat}}{{\tau }_{ch}}))},$$

(12)

where ${\tau }_{lat}$ and ${\tau }_{ch}$ are defined as (with i = $lat$, $ch$)

$${\tau }_i(T,S_w)=\left[{\tau }_i(\overline{T},1)+\frac{{\tau }_i(\overline{T},0)-{\tau }_i(\overline{T},1)}{1+c_{1i}exp\left[-\frac{c_{2i}(1-S_w)}{S_w(1-S_w)}\right]}\right]exp\left[U_i\left(\frac{1}{T}-\frac{1}{\overline{T}}\right)\right],$$

(13)

where $U_{ch}$ and $U_{lat}$ [K] denote the Arrhenius activation energies, representing the energy that must be reached for a chemical reaction to occur.

For simplicity, and in view of the specific application of this paper, we only consider the activation of damage for prevaling tensile stress state, instead of the more general bi-dissipative damage model [26]. The evolution of the tensile damage variable D is governed by the loading–unloading conditions, defined in terms of $\mathit{\sigma }$ and p through the inelastic effective stress ${\mathit{\sigma }}^{\mathbf{\prime }\mathbf{\prime }}$ (with $0<\beta \le b$):

$${\mathit{\sigma }}^{\mathbf{\prime }\mathbf{\prime }}=\mathit{\sigma }+\beta p\mathbf{1},$$

(14)

$$f=\frac{1}{2}\mathbf{s}:\mathbf{s}-a_t{\left(tr{\mathit{\sigma }}^{\mathbf{\prime }\mathbf{\prime }}\right)}^2+b_t\left(tr{\mathit{\sigma }}^{\mathbf{\prime }\mathbf{\prime }}\right)h-k_t{h}^2,$$

(15)

$$f\le 0\dot{D}\ge 0f\dot{D}=0,$$

(16)

where f is the damage activation function; $\mathbf{s}$ is the deviatoric stress; $a_t$, $b_t$, and $k_t$ are material parameters governing the shape and dimensions of the elastic domain (see [26] for details); and h is the softening function, defined as

$$h\left(D\right)={[1-{\left(D\right)}^c]}^{0.75}.$$

(17)

The parameter c in Equation (17) tunes the specific fracture energy of the model.

3. Model Validation

The model used in this work, already employed in [18,25] to simulate the experimental campaigns in [21,27], is here further validated by simulating the recent experimental campaign of [19,28]. In these works, the expansion of cubic specimens (of 254 mm side) of concrete affected by ASR and subject to multi-axial stresses was measured. The post-tensioning method of pre-stressing was chosen as the technique to apply and sustain the desired compressive stresses in the concrete specimens. The stress application in one direction required four high-strength bolts to be inserted into four rigid sheaths embedded in the concrete during casting. Figure 2a shows the experimental configuration for a three-axial stress-state application. Each bolt was stretched to a desired stress level by producing a reaction against a pair of bearing plates placed symmetrically on the opposite surfaces of a concrete specimen. Different dimensions of bolts allow for obtaining the average nominal compression stress in the concrete, in a range of 0–9.6 MPa.

Four (three reactive and one non-reactive to serve as the control specimen) cubic specimens were cast for each of the different stress states. Expansion for reactive specimens due to ASR can be decoupled from the creep and shrinkage strains by taking the algebraic difference between the expansion of a reactive specimen and the corresponding expansion of the control specimen (non-reactive concrete) at an identical stress state. We will refer to these data as ASR-only results. The stress was applied for 52–56 days after casting at room temperature. After 180 days, the specimens were conditioned inside an acceleration chamber maintained at 50 ${}^{\circ }$C and close to 100% relative humidity to promote and accelerate ASR. During these tests, the strains were monitored in the three directions as the average of the four measurements.

To simulate these experimental results, we assumed the mechanical parameters of concrete given in [19] and calibrated the material parameters governing the ASR evolution with the difference between the free expansion curves of reactive and non-reactive specimens (ASR-only data). Other material parameters were defined according to reference values, such as the porosity $\varphi =0.16$ (from which it is possible to evaluate Biot’s coefficient $b=0.41$) and Biot’s modulus M= 5000 MPa. A comparison between experiments and the model result is shown in Figure 3, in terms of mean strain evolution. It should be noted that, in these tests, the temperature conditions changed while the degree of saturation was constant and equal to 1; therefore, according to Equation (13), the latency and characteristic times changed. The following values were identified at 23 ${}^{\circ }$C: ${\tau }_{lat}$ = 478 days and ${\tau }_{ch}$ = 274 days while, at 50 ${}^{\circ }$C, we have ${\tau }_{lat}$ = 34 days and ${\tau }_{ch}$ = 60 days. The identified asymptotic volumetric strain is ${\mathit{\epsilon }}_{ASR}^{\infty }$ = $4.9\times {10}^{-3}$.

In order to simulate the loaded tests in a realistic way, a finite element (FE) mechanical analysis was carried out. In fact, the stress state within the specimens induced by the pre-stressed bolts was non-homogeneous. In Figure 2b, the finite element model used for the simulation of the experimental tests is shown (as the problem is symmetric, only one eighth of the cubic specimen was modelled), in which the plates were modelled and the load was applied in the circular regions representing the bolts. With the above parameters calibrated with the free exansion test, we simulated the tests with uniaxial compression and biaxial compression. In Figure 4, comparisons between the experimental data and numerical results are represented in terms of strain evolutions in the three directions. Even though the material model was isotropic, the load-induced anisotropy of the expansion curves was well predicted.

The difference between the expansion in the X and Y directions in the bi-axial loaded case (Figure 4b) was due to the different locations of the loading plates on the two faces of the cube, which induced a stress state far from being equi-biaxial at each point. The numerical simulation was able to predict this complex stress state and the consequent non-uniform damage pattern (see Figure 5).

4. Results and Discussion. Structural Analysis of a Concrete Dam Affected by Alkali–Silica Reaction

4.1. Numerical Approach

The model presented above allows us to perform numerical simulations, using the finite element method, of the mechanical consequences of ASR on real concrete structures.

The approach we followed is weakly coupled. First, the histories of degree of saturation and temperature are computed by solving the moisture diffusion problem (Equation (4)) and the heat diffusion problem (Equation (6)) with corresponding boundary conditions, possibly varying in time. Then, the reaction extent evolution can be evaluated, through Equation (10). As this evolution depends on temperature and humidity, the development of the reaction is non-uniform inside the structure; in particular, the reaction is faster in warmer parts of the dam and expansion is higher in parts where the degree of saturation is higher.

The reaction extent, together with the external loads and temperature variations, is the input of the elasto-damage analysis, which allows us to compute the evolution of displacements and the damage patterns leading to the formation of cracks.

4.2. Case Study: Gravity Dam in Canada

Reference is made to the left wing of the gravity dam of the Beauharnois power plant (Québec, Canada), shown in Figure 6, courtesy of Hydro-Québec. The construction of the whole plant begun in 1932 and the gravity dam, considered here, was built between 1956 and 1960. It is 21 m high (h) and it is characterized by a base width of b = 15 m and a crest dimension of a = 3 m (Figure 6b). The mechanical properties of the concrete and foundation provided by Hydro-Québec are collected in

Table 1. Mechanical and thermal parameters of concrete and rock for the Beauharnois dam.

Parameter	Value	Unit
$E_{concrete}$	26,000	MPa
${\nu }_{concrete}$	0.21	-
${\rho }_{concrete}$	2365	kg/m${}^3$
$f_{ct}$	3	MPa
${\alpha }_{concrete}$	1 $\times {10}^{-5}$	1/${}^{\circ }$C
$E_{rock}$	30,000	MPa
${\nu }_{rock}$	0.2	-
${\alpha }_{rock}$	1 $\times {10}^{-5}$	1/${}^{\circ }$C

. The fracture energy in tension typical of concrete characterized by a large aggregate size was also given: G${}_{ft}$ = 350 N/m.

Figure 6a shows the locations of the various instruments for monitoring. The monitoring system for the dam-foundation is represented by topographic measurements along dam crest through the reference targets shown by black dots in Figure 6a. Displacements in the three directions (X, Y, and Z) are measured at all reference targets. Due to ASR, an increase in the mean annual crest displacements has been observed since 1990.

Furthermore, a set of thermometers determined the temperature of the concrete in contact with air. These values, recorded over six decades, have been interpolated by a periodic function with a period of one year. Figure 7a shows the annual profile thus obtained: $T\left(t\right)=\overline{T}+\Delta Tg\left(t\right)$, where $\Delta T=T_{max}-T_{min}$ and $g\left(t\right)$ is the time function with period of one year.

The reservoir-level oscillations were very limited (with a maximum oscillation of 0.2 m, from a maximum value of 46.1 m to a minimum of 45.9 m). For this reason, a mean value of 46 m was considered in all analyses.

The degree of saturation inside the dam was, however, non-constant due to the variation of relative humidity of the air. The annual variation of relative humidity of the air is given in Figure 7b.

In view of the geometry, plane strain conditions were assumed and a 2D mesh with triangular (3 nodes) elements was employed, as depicted in Figure 8 (for the dam foundation system). Firstly, the temperature and humidity effects were evaluated using the diffusion analyses described previously. Then, the mechanical analysis of the studied dam was carried out, in order to compare the structural response obtained by the model with the real monitoring data.

To determine the history of the temperature field within the dam, a heat diffusion analysis (by solving Equation (6) with proper initial conditions and varying boundary conditions) should be performed. The boundary conditions (BC) were set as follows (see Figure 9). The temperature of the dam surface in contact with air (the downstream face and the upstream face above water) was assumed equal to the air temperature, whose annual variation is known (as previously discussed) and shown in Figure 7a. The upstream surface below the reservoir level, assumed constant and equal to 46 m, was subject to a temperature linearly varying with the depth (as proposed in [29]), with an annual oscillation similar to that of the air, $T\left(t\right)=\overline{T}+\Delta Tg\left(t\right)$, but with different maximum and minimum values (from 1 ${}^{\circ }$C to 15 ${}^{\circ }$C at the water surface and from 3 ${}^{\circ }$C to 7 ${}^{\circ }$C at the dam foot, as depicted in Figure 9). The rock temperature at the external boundary was assumed constant and equal to 4 ${}^{\circ }$C. The specific heat capacity and thermal conductivity of concrete, given by Hydro-Québeq, were C= 917 J/kg${}^{\circ }$C and $D_T$ = $2.9$ W/m${}^{\circ }$C, respectively. These values were all within the quite large range of variability of concrete thermal properties [30,31].

The procedure followed to determine the temperature at each point of the dam, at each time of the year, and then periodically repeating, $T(\mathbf{x},t)$, is schematically shown in Figure 10.

In the first step, fixing the temperature at different portions of the boundary of the dam as equal to the mean value over one year, $\overline{T}$, in air, water, and rock, a steady state analysis was used to evaluate the stabilized temperature in the internal nodes of the mesh $T(\mathbf{x},0)$, which was used as initial conditions for the following transient analysis.

In Figure 11a, the contour plot of the stabilized initial temperature $T(\mathbf{x},0)$ is shown.

The second step consisted of a thermal transient analysis, in which the real periodic annual variation of air and water temperature was applied at the buondary $T=\overline{T}+\Delta Tg\left(t\right)$. The analysis was performed until a stabilized periodic pattern was obtained throughout the whole body of the dam. This stabilized periodic temperature field was used as input for the following chemo-mechanical analysis. In Figure 12a, the yearly histories of temperature at different points of the dam (A = in contact with air, B = close to the reservoir level, and C and D = internal nodes; see Figure 9) obtained through this second transient analysis are depicted.

The evaluation of the yearly evolution of the degree of saturation field followed the same procedure. As in the case of thermal analyses, a preliminary moisture diffusion analysis was performed before the transient moisture diffusion analysis.

For the first analysis, we assumed an initial uniform field for the degree of saturation, $S_w$=$0.85$, and boundary conditions corresponding to a water level of 46 m, such that the submerged nodes were characterized by a complete saturation, while, for the nodes in contact with air, the mean of the sinusoidal annual variation (equal to 0.66 and corresponding to RH = 84%) was applied. The parameters $m_1$, $m_2$, and $m_3$ of the permeability were assumed, as in [22]; namely $m_1$ = $18.62$ MPa, $m_2$ = $2.27$, and $m_3$ = $0.44$. The obtained pattern of saturation degree is shown in Figure 11b. Then, the transient moisture diffusion analysis was performed, considering the annual oscillation of Figure 7b. Figure 12b reports the obtained stabilized yearly variations of the saturation degree at different nodes. It should be noted that the moisture diffusive process affects the external boundary of the structure only to a depth of about 20 cm (i.e., the depth of one element) and the degree of saturation is constant elsewhere, equal to the imposed initial value of 0.85. Therefore, the effect of humidity on the ASR evolution is expected to have a marginal role in the structural response.

The reaction extent evolution was then computed by integrating Equation (10). No specific tests of ASR expansion on concrete specimens of this dam were available. Therefore, we conducted several analyses by changing the latency, the characteristic times, and the asymptotic ASR expansion within the ranges previously established (e.g., in [18]) and, finally, we arrived at the values given in

Table 2. Parameters for the alkali–silica reaction (ASR) model used for the simulation of Beauharnois dam.

Parameter	Value	Unit
${\tau }_{lat}(\overline{T},1)$	500	days
${\tau }_{lat}(\overline{T},0)$	1100	days
${\tau }_{ch}(\overline{T},1)$	350	days
${\tau }_{ch}(\overline{T},0)$	700	days
${ϵ}_{ASR}^{\infty }$	0.01	-

.

Figure 13 shows the evolution of the reaction extent for the points A, B, C, and D of Figure 9. The reaction extent developed more quickly at nodes A and B, located in the upper part of the dam, where the temperature was higher. The influence of the different degrees of saturation between points A and B was very limited. The temperature and the reaction extent fields were then the input into the chemo-thermo-mechanical analysis of the dam, aimed at evaluating the response due to ASR.

The development of ASR produced an overall expansion of the dam, with a permanent variation of upstream–downstream displacements, which was superposed to the normal seasonal oscillation. Furthermore, ASR induced damage in the material. Figure 14 shows the damage at two different time steps. At the beginning, the damage appeared only on the external skin of the structure; the only part affected by the external humidity conditions. This was the first manifestation of ASR development, due to the different expansions at points with different degrees of saturation. The first macroscopic crack in the body of the dam and visible on the surface of the structure in contact with the rock developed 30 years after construction. One should observe that, as expected in a real massive structure at moderate environmental temperatures, the time scale involved in the damage process due to ASR was completely different from that of one typical in a standard test on mortar specimens accelerated by high temperature (see, e.g., [32]). Figure 15 shows a comparison between the monitoring data (mean values at points $P_1$ and $P_2$ of Figure 6a) and the numerical results, in terms of horizontal crest displacement, positive in the downstream–upstream direction. A good agreement between them was obtained.

5. Conclusions

In this work, a bi-phase damage model for concrete subject to the alkali-silica reaction present in the literature was used to simulate a recent experimental campaign, as well as to perform a mechanical analysis of an existing concrete gravity dam. In the adopted model, based on Biot’s theory of multi-phase materials, the progressive expansion of the alkali–silica gel is described by a phenomenological internal variable, depending both on temperature and on degree of saturation.

The structural analysis of an existing concrete gravity dam was aimed at assessing the extent of structural damage caused by the reaction. The material parameters used in the model to define the mechanical behaviour and the reaction development were calibrated and fixed using the monitoring information available, provided by Hydro-Quebec.

Firstly, the effects of temperature and of humidity on the hydraulic structure were evaluated through a heat diffusion analysis and a moisture diffusion analysis, respectively. The moisture gradient was present only near the dowstream face of the dam and, so, its effect on the ASR evolution was limited. Then, the solutions of these two analyses were used as input for a consequent chemo-thermo-mechanical analysis, used to define the response due to ASR. The results, in terms of horizontal crest displacement, were in good agreement with the available real structural data.