Chapter 4

NUMERICAL MODELING OF DRYING

SHRINKAGE OF CONCRETE SPECIMENS

In this chapter the modeling of drying shrinkage of concrete specimens is addressed.

First, in section 4.1, the moisture diffusion model used in the simulations is described in
detail for the continuous medium as well as for the cracks, and the coupling strategy is
then discussed in section 4.2. Thereafter, the most relevant modeling results obtained are
described in detail and discussed and possible future research directions are explored.
Preliminary studies of the effect of a single crack on the drying process and the effect of
a single inclusion on the drying-induced microcracking have been performed in order to
provide tools to understand the effects of the main ingredients of the mesostructural
approach on drying shrinkage. A large number of HM simulations at the meso-level
have been calculated for evaluating the ability of the present mesostructural model to
capture the main features of drying shrinkage in concrete specimens. Model parameters
have been adjusted with experimental results (Granger, 1996). Finally, a preliminary
study of the effect of drying under simultaneous loading has been carried out.

4.1. Drying shrinkage: model description

In order to extend the applicability of the mechanical model defined in Chapter 2 for
the case of hygro-mechanical (HM) problems, such as drying shrinkage in concrete, it is
necessary to establish a connection between the purely mechanical analysis and the
diffusion-driven drying process, as defined in section 3.4.3. The origin of these coupling
sources has been discussed in detail in Chapter 3 and may be summarized as follows. On
one hand, the moisture loss due to drying induces a volumetric shrinkage which, if
restrained, generates a self-equilibrated state of stresses within a material cross-section
(see figure 3.11 in Chapter 3). It is reasonable to assume that the induced strain level will
have very little impact on the porosity of the sound material (and hence on its
diffusivity), so that this effect can be neglected. On the other hand, an acceleration of the
drying process is to be expected in the case of micro/macro cracking (see section 3.1.5 in
Chapter 3). This last consideration is usually neglected in most of the models proposed
in the literature, mainly due to the complexity of a correct implementation and the
difficulties encountered for validating a theoretical model, as a result of the scarcity of
experimental data and a lack of a standardized procedure for determining moisture
transfer through a crack (see section 3.1.5). In this section, the moisture diffusion model
through (cracked) concrete adopted throughout this thesis is described in detail.

85
4.1.1. Moisture diffusion through the uncracked porous media
Since the early work of Bazant & Najjar (1972), it is generally accepted that, at least
as a first approximation, moisture movement in concrete basically follows a non-linear
diffusion-type equation which may be advantageously written in terms of the relative
humidity (RH) at the point, denoted as H (varying between 0 and 1). This expression has
been derived in Chapter 3 (eq. 3.37) and is rewritten here for convenience

= −div ( Deff ( H ) grad ( H ) )

dH
(4.1)
dt
where the effective diffusion coefficient Deff [cm2/s] strongly depends on H itself. In this
thesis, following a previous work within the research group (Roncero, 1999), this
dependency has been represented as:
Deff ( H ) = D0 + ( D1 − D0 ) f ( β ,H ) (4.2)

In the previous expression, D0 and D1 are constants determining the value of the
coefficient at zero RH and at fully saturated condition, respectively, and the dependence
on RH is introduced by a hyperbolic function written as
e −β ⋅ H
f ( β ,H ) = , with β = shape factor (4.3)
1 + ( e −β − 1 ) ⋅ H
The dependence of the diffusivity on RH is shown in figure 4.1 for different values of
the shape factor. A value of 3.8 was shown to fit experimental data on mortar specimens
(Roncero, 1999). Boundary conditions for the diffusion equations have been already
discussed in Chapter 3 (eqs. 3.38-3.40). A sensitivity study of the diffusion coefficient
model parameters has been performed elsewhere (Roncero, 1999) and will not be
repeated in this work.

0.09
D1
D H (H)

β =0
=1
=2
=3
= 3.8
D0
0
0 RH 1

Figure 4.1. Relation between the diffusion coefficient and RH for different shape factors
β, as proposed by Roncero (1999).

In this thesis, dealing with a mesostructural representation of concrete materials, the

moisture diffusion has been assumed to take place only through the matrix (which in
turn represents mortar plus smaller aggregates, see Chapter 2), since the discretized
larger aggregates are considered to have a negligible diffusivity with respect to the
matrix one (thus remaining fully saturated throughout the simulation). However,
diffusion through the aggregates could readily be considered if needed, as for instance in
the case of lightweight aggregates, having a porosity as high as that of mortar.

86
4.1.2. Moisture diffusion through the cracks
Cracks may affect diffusivity since they represent potential preferential channels for
moisture migration out of the material. In the case that the mechanical analysis
determines the formation and/or propagation of a crack, the diffusion process will thus
be altered. To analyze this problem, the same FE mesh is used for both the diffusion and
the mechanical calculations, thus simplifying the coupled HM analysis, which is
achieved by the use of a staggered strategy (see next section). This has required the
formulation and implementation of interface elements with double nodes also for the
diffusion analysis (Segura & Carol, 2004; Segura, 2007; Segura & Carol, 2008). Such
elements incorporate longitudinal (KL [cm2/s]) as well as transversal (KT [cm2/s])
diffusivities. The existence of a discontinuity may also represent an obstacle to fluid
flow in the transversal direction, due to the transition from a pore system into an open
channel and back into a pore system, or to the existence of infill material (as in the case
of precipitation of species). This resistance complicates the flow in the transverse
direction and results in a localized potential drop across the discontinuity (Segura, 2007).
In the absence of specific information, the transversal diffusivity is given a very high
value, representing the case of no jump in RH across the crack (see eq. 4.9). For the
longitudinal diffusion two situations are distinguished: before the interface starts
opening, KL takes a zero value (optionally a fixed value Kmin may be chosen, with which
the higher porosity of the aggregate-mortar interface could be simulated), and after the
crack has opened an expression similar to eq. 4.2 is used, in which the diffusivity for
saturated flow K1 [cm2/s] is in this case given by the so-called “cubic law” (see section
3.4.5 in Chapter 3) and K0 [cm2/s] is set as a small fraction of K1 [cm2/s]:
K L ( H ,u ) = K0 ( u ) +  K1 ( u ) − K0 ( u )  f ( β K ,H ) (4.4)

K1 ( u ) = η ⋅ u 2 [cm2/s], and K0 ( u ) = α K1 ( u ) (4.5)

In eqs. 4.4 and 4.5, u [cm] is the crack opening (or crack width), η [1/s] is a parameter
of the model relating the square of the crack width with the diffusivity, to be determined,
α and βK are two model parameters and f is a function given by eq. 4.3 (with β replaced
by βK). The cubic law is recovered when plugging in the previous relation into the
framework described in the following.
Considering a local orthogonal coordinate system (x,y) within a discontinuity, the
mass conservation equation for an incompressible fluid filling the discontinuity may be
written as (Segura, 2007)
dĴ l ∂u
− + q− + q+ = (4.6)
dx ∂t
where Ĵ l is the longitudinal local flow rate and q- and q+ are the leak-off terms,
incoming to the discontinuity from the surrounding porous medium. The total
longitudinal flow along the discontinuity is derived from the first Fick’s law and yields
dH mp
Ĵ 1 = u ⋅ J 1 = −u ⋅ K L (4.7)
dx
in which J l is the local moisture flux, and Hmp represents the RH at the mid-plane of the
fracture (see figure 4.2). Plugging in eq. 4.7 into eq. 4.6 leads to the partial differential
equation governing longitudinal flow along the mid-plane of the crack, yielding

87
 d  dH mp  − + ∂u
 u ⋅ KL +q +q = (4.8)
dx  dx  ∂t
Note that in the last expression, the cubic law is recovered when multiplying KL by u,
yielding η ⋅ u 3 . In the case of the transversal flow, a total potential drop across the
discontinuity is considered as
qT = KT ⋅ Hˆ with Ĥ = H − − H + (4.9)
It should be noted that the hypothesis of an incompressible fluid filling the crack may
not hold for the case of RH, and eq. 4.6 might need to consider an extra term for the
compressibility of the vapor phase. Nevertheless, in the context of this thesis, a similar
formulation as in saturated flow has been used for a first evaluation and comparison with
experimental data, and the refinement of the model with extra terms has been left for
future research.
Equations 4.6 to 4.9 are the ones to be discretized through the FEM. The details of
such discretization may be found elsewhere (Segura & Carol, 2004; Segura, 2007).

Figure 4.2. Longitudinal and transversal flow through a differential joint element (from
Segura & Carol, 2004).

4.1.3. Desorption isotherms model (Norling model, 1994)

An important ingredient for the modeling of the drying process are the so-called
desorption isotherms (see section 3.1.3), relating RH (H) to evaporable water content in
the pores we [g/cm3], this last needed for calculating the overall specimen weight loss but
also for predicting the “shrinkage at a point” (i.e. as a material property). In this study,
the local weight losses are calculated as a post-process, after computing the RH field.
However, in the case that the formulation accounting for a nonlinear moisture capacity
matrix is adopted, the desorption isotherm will have an effect on the RH field (Xi et al.,
1994). Throughout this thesis, the desorption isotherm proposed by Norling (1994) has
been adopted for this purpose, as in the work of Roncero (1999). In that model, it has
been proposed to construct the desorption isotherms as a sum of a gel isotherm and a
capillary one, supported by the fact that the physically bound water (to differentiate it
from the chemically bound water) exists in the gel pores and in the capillary pores
(Norling, 1997). This relation is written in equations 4.10 and 4.11 and is shown in
figure 4.3, for typical values of the model parameters.
we ( H ) / c = x ⋅ (1 − e − z ⋅ H ) + y ⋅ ( e z ⋅ H − 1) (4.10)

88
 0 .15 ⋅ α w / c − 0 .33 ⋅ α
with z = − ( w0 / c ) ⋅ f1 + f 2 , x = −z
and y = 0 (4.11)
1− e ez −1
In the above eqs., c is the cement content [g/cm3], x, y and z are functions of the degree
of hydration (α) and the initial water-to-cement ratio (w0/c), and f1 and f2 represent two
shape factors (Roncero, 1999). Again, a sensitivity study of the different parameters on
the overall weight losses can be found elsewhere (Roncero, 1999).
The derivative of the desorption isotherm is given by equation 4.12 and plotted in
figure 4.3. In case that a formulation based on equation 3.36 is preferred rather than
3.37, the moisture capacity matrix C(H) needs to be computed as the derivative of the
desorption isotherm. In section 4.6, results are presented in which a comparison between
the formulations given by eqs. 3.36 and 3.37 has been performed.
dwe ( H )
= c ⋅ ( x ⋅ z ⋅ e − z⋅H + y ⋅ z ⋅ e z⋅H ) (4.12)
dH

Figure 4.3. Desorption isotherm according to the model proposed by Norling (1994) and
corresponding moisture capacity (derivative of the curve with respect to RH) for the
following model parameters: w0/c = 0.5; α = 0.8; c = 0.473; f1 = 5; f2 = 7.

4.1.4. Volumetric strains due to drying

As stated in Chapter 3 (section 3.4.3), the difficulties encountered for estimating the
shrinkage coefficient, relating volumetric strains to weight losses (or RH variations) at a
local level, have lead most researchers to traditionally assume a constant value of this
coefficient, as shown in equation 4.13.
ε shr = α shr ⋅ ∆ we (a) or ε shr = α shr ⋅ ∆ H (b) (4.13)
Physically, it is clear that variations in moisture content affect mechanical behavior
through the material shrinkage at the point which acts similarly as a non-uniform
thermal contraction, more pronounced in the parts of the specimen exposed to drying.
But reality is hardly ever that simple, and a linear relationship between weight losses and
strains may yield only rough approximations. In order to overcome this problem,
Alvaredo (1995) estimated different shrinkage coefficients at various RH levels by
inverse analysis. Later on, based on these findings, van Zijl (1999) proposed to treat the
shrinkage coefficient as a linear or hyperbolic function of RH. In the present work,
shrinkage at the point has been assumed to be related linearly to the water loss per unit
volume (i.e. constant shrinkage coefficient) at each point in a first stage (López et al.,
2005a,b). Calibration of model parameters with experimental results has put in evidence

89
the advantages of considering a nonlinear relation between strains and weight losses, as
will be shown in this chapter (section 4.6). The dependency adopted in this thesis is
based on the work by van Zijl (1999). The difference is that a relation between strain and
weight losses has been preferred here, instead of the RH variation. The local volumetric
strains can be integrated as follows

∫ ε& dt = ∫ α ( w ) w& dt + cst

s s e e (4.14)

in which the shrinkage coefficient, for the cases of constant value, and linear or
quadratic dependency on the weight loss, is respectively expressed as

 ε s∞

 weenv − we0
 2ε s∞
α s ( we ) =  we (4.15)
(
 ew env 2
) ( )
− we
0 2


 3ε s∞
 env 3 we 2
(
 we ) ( )
− we 0 3

In eq. 4.15, ε s∞ represents the final volumetric shrinkage strain corresponding to the
environmental RH considered, weenv and we0 are the moisture contents corresponding to
the environmental RH and the initial internal RH of the sample, respectively (they can
be computed from the desorption isotherm), and we represents the average moisture
content within the considered time interval. These expressions have been derived by
forcing the same (fixed) final drying-induced shrinkage value for all cases (see figures
4.4 and 4.5). In this way, only the velocity at which drying occurs (or drying rate) is
altered when considering different relations.
Figure 4.4 shows the effect, at a constitutive level, of considering a linear or quadratic
dependence of the shrinkage coefficient on the weight loss, as opposed to the constant
shrinkage coefficient case. It can be seen that considering a nonlinear relationship
accelerates the drying-induced shrinkage strains. Note that the final deformation, which
in this case corresponds to that reached at equilibrium at 50% RH, is the same in the
three cases.
When compared to figure 3.4 in Chapter 3, which shows the same relationship
determined experimentally in small specimens (from Bazant, 1988 and Helmuth & Turk,
1967), it can be seen that considering the linear or quadratic dependence of the shrinkage
coefficient on the weight loss does not mean a better fit of the experimental results. On
the contrary, the calculated trends are the opposite of the experimental ones (the model
predicts larger deformations at the beginning of drying). This difference is not so
surprising and can be partly explained by the skin microcracking effect on shrinkage
strains, which could not be completely hindered in the experiments. This effect is
obviously not considered in the model given by eq. 4.15, since cracking phenomena are
explicitly taken into account via zero-thickness interface elements. It will be seen in
section 4.6 that the lower rate of deformation at the beginning of drying is due to skin
microcracking, independently of the adopted relationship.
Figure 4.5 presents the same effect for the case of a simple FE calculation (mesh size:
12.5x1cm2; upper face exposed to HR=0.5, and sealed conditions for the rest), showing

90
the evolution of shrinkage strains with time for the same three cases proposed in eq.
4.15. It can be seen that the strain rate decreases more slowly in the case of parabolic
shrinkage coefficient.

Figure 4.4. (a) Comparison of different relationships, at a local level, between shrinkage
strains and weight losses: cases of constant, linear and quadratic shrinkage coefficients.
The black vertical line represents the weight loss corresponding to a 50% RH
environment. (b) Shrinkage coefficient vs. water content for the same three cases.

Figure 4.5. Influence of considering different relationships between strains and weight
losses in a FE calculation: RH distribution and deformed mesh (left) for t = 200 days
(scale factor = 200); shrinkage strains (measured at the red points in the deformed mesh)
as a function of time for the cases of constant, linear and parabolic shrinkage
coefficients.

A more detailed inspection of the proposed model for shrinkage strains has been
performed by replacing the weight loss variable in the horizontal axis by the RH in the
previous analysis with the use of the desorption isotherms. In this case there is obviously
an additional nonlinearity introduced by the shape of the isotherms. The results of the
model when drying takes place from a saturated condition to an almost zero value of
RH, have been qualitatively compared with experimental results in thin concrete slices
(Baroghel-Bouny et al., 1999). Experiments were performed in 3mm thick concrete
discs (max. aggregate size of 20mm) with a diameter of 90mm, which are suitable for
comparing with the behavior of the matrix (representing mortar plus smaller aggregates)
in the present mesoscale model. Changes in diameter were recorded at different RH
levels. Although the RH gradients are minimized with the use of very thin slices, it

91
should be noted that even for 3mm thick specimens the skin microcracking is not
completely prevented (Hwang & Young, 1984). Thus, experimental results obtained by
this method do not represent exactly the real unrestrained shrinkage strain.

Figure 4.6 presents the results of this comparison. First of all, it can be seen in figure
4.6a that the linear relation between strains and weight losses for the case of constant
shrinkage coefficient in figure 4.4a is represented by a nonlinear one when plotting
strains against RH, depending on the desorption isotherm model. Moreover, the
comparison with the experimental results in figure 4.6b shows that a constant shrinkage
coefficient may be more suited for high strength concretes, whereas ordinary concrete
behavior is better captured with a quadratic shrinkage coefficient.

Figure 4.6. (a) Relation between longitudinal shrinkage strains (taken as 1/3 of the
volumetric strain) and RH obtained with the model when drying from saturation to
almost zero RH (with the use of desorption isotherms) for the cases of constant, linear
and quadratic shrinkage coefficients (a value of 0.01 [cm3/g] is assumed for the constant
shrinkage coefficient). (b) Experimental results on circular concrete slice (3mm thick),
for ordinary concrete (OC) and high performance concrete (HPC) (data by Baroghel-
Bouny et al., 1999; taken from Gawin et al., 2007).

4.2. Coupling strategy: a staggered approach

Performing a coupled analysis implies a mutual interaction between two or more
processes, so that an interrelation between each of these exists. In the case of the hygro-
mechanical (HM) coupling dealt with in this chapter, the deformation of the porous
medium (cracking phenomena and changes in porosity) may affect the moisture
diffusion mechanisms of the system and, on the other hand, the humidity distribution
given by the moisture diffusion problem affects the volumetric strains distribution in the
mechanical analysis. This problem is said to be coupled, since any of the two processes
may induce an effect on the other.
The resulting system of equations is usually as shown in eq. 4.16 for the case of HM
coupling, in which u and H are, respectively, the nodal displacement vector and the
nodal RH values, and F and Q play the role of external forces and imposed RH flux
vectors, respectively. The sub-matrices K12 and K21 represent the coupling sources of
both processes. This approach is similar to the well known u-p formulation (u =
displacements; p = pressures) in Geomechanics (Zienkiewicz et al., 1977).
 K11 K12   u   F 
K = (4.16)
 21 K 22   H  Q 

92
 There exist two different coupling procedures for solving a coupled nonlinear system
of equations. These are the monolithic or fully coupled approach and the staggered
approach (see Segura & Carol, 2008 for a detailed comparison between them in the
framework of a u-p formulation, and introducing zero-thickness interface elements). In
the first case, the HM behavior is described by a group of equations that is solved
simultaneously and that incorporates all the relevant physics of the problem. In the case
of a staggered approach, the solution is obtained, for each time step, from the results of
each problem separately, but in which the information is iteratively transferred between
the two problems (generally, two separate codes), until a certain tolerance is satisfied
(the results obtained from one problem are the input for the other).
It is generally accepted that the monolithic approach is required in cases of strong
coupling such as for instance hydraulic fracture analysis, while staggered approach can
also lead to a solution in cases of relatively weaker coupling [Segura and Carol, 2008],
such as the drying shrinkage or sulfate attack problems studied in this thesis. On the
other hand, the staggered approach also exhibits some advantages over the fully coupled
one. From a practical point of view, it is possible to use existent standalone codes which,
although originally designed for uncoupled analyses, are often powerful tools for
tackling each problem separately (in the present case, the mechanical and the diffusion
analyses). It also permits the separation of the time scale for each problem, making it a
very suitable tool for cases in which the time scale for the diffusion problem is very
different from that required for the mechanical simulation. This is in contrast with the
monolithic approach, in which the solution of the overall problem is governed by the
most critical time scale.
In the proposed model, the hygro-mechanical (HM) coupling is achieved through a
staggered approach, as shown schematically in figure 4.7. One code (DRACFLOW)
performs the nonlinear moisture diffusion analysis, and the results in terms of volumetric
strains at the local level (i.e. at the integration points) serve as input to the second code
(DRAC, Prat et al., 1993), i.e. for the mechanical problem. This latter code obtains,
among other results, an updated displacement field (nodal variables) from which new
crack openings may be derived, which in turn will alter the diffusion analysis. This loop
is successively repeated within each time step until a certain tolerance is satisfied (in
terms of water losses), before passing to the next time interval. Both codes have been
developed within the research group of Mechanics of Materials at UPC (Prof. Ignacio
Carol and Prof. Pere Prat).
As already mentioned, the staggered approach is well-known to sometimes yield
convergence problems in the case of the u-p formulation. However, for the present case
of u-H, this coupling strategy works perfectly and there has been no single case of non-
convergence of the coupling iterations during this work. This could be mainly due to two
reasons. First, the ‘coupling intensity’, given by the magnitude of the off-diagonal terms
in matrix K, has not been significant in the calculations performed so far. Note that these
off-diagonal terms destabilize the system, since an increase of these last causes a
reduction of the determinant of the matrix. Secondly, in the u-p formulation the results
from the hydraulic analysis are expressed in terms of pressures, which enter the
mechanical problem as imposed forces. This could bring convergence problems in the
case of elevated pressure values. On the other hand, the output from the diffusion model
adopted in the present work is expressed in terms of imposed volumetric shrinkage
strains, which is expected to favor the convergence of the system.

93
Figure 4.7. Schematic representation of the coupled hygro-mechanical (HM) system
through a staggered approach.

4.3. Preliminary study of the effect of a single crack on the drying

process
In Chapter 3, a detailed description of the experimental evidence on the effect of
microcracks on the drying process was addressed (see section 3.1.5). It has been shown
that, although it remains an open issue, there is a non-negligible contribution in many
practical and experimental cases. In this work, it has been decided to analyze first the
influence of a single crack in the drying process. More specifically, a sensitivity study of
two key parameters have been performed, namely the crack width and the crack depth.
This section presents the main results of this preliminary study, obtained with the
diffusion model described in the previous section and implemented into the finite
element code DRACFLOW (no mechanical analysis is performed in this stage).
The finite element mesh used in all the calculations is shown in figure 4.8. One single
central line of zero-thickness interface elements has been introduced (then, depending on
the diffusive properties of the individual interface elements along this line, different
crack depths and widths can be readily simulated). This configuration corresponds to a
crack spacing of 3cm, which is a reasonable value to adopt given the crack penetrations
studied in this section (see section 3.1.5.4. in Chapter 3). Boundary conditions consist of
sealing all faces. Additionally, a RH of 50% is imposed at point P in the figure. In this
way, moisture is forced to escape only through this point, either through the continuous
medium or through the active crack. Initial condition is HR = 100% throughout the
sample. The desorption isotherm and diffusivity dependence on RH used in this study
are shown in figure 4.9.

94
Figure 4.8. (a) FE mesh and boundary conditions used throughout the simulations: HR =
0.5 at point P. The middle line (in red) represents the joint elements with diffusive
parameters varying depending on the case studied (in terms of crack width and crack
depth). (b) Configuration for the case of 2cm crack depth, showing two sets of interface
elements.

Figure 4.9. (a) Diffusivity vs. RH and (b) desorption isotherm used in the simulations of
the influence of a single crack on the drying process.

Two series of simulations have been considered: on one hand, the influence of the
crack depth at constant width (50 microns in this study) has been addressed, and on the
other, the influence of the crack width for a fixed crack depth (chosen here as 2cm).
Both series have been compared to a ‘reference’ case, in which there is no crack and
therefore moisture escape occurs only through a point (point P) due to the flux coming
from the surrounding medium (as in figure 4.10a). In all the cases the drying process has
been quantified in time by the total weight loss of the sample and the RH distribution. In
order to consider different crack depths over the same FE mesh, two sets of material
properties for interface elements are defined: one with the desired crack opening and the
second one with zero crack width (therefore nil diffusivity). The different crack widths
automatically fix the diffusivity value of the cracks in each case via the cubic law (eqs.
4.4-4.5).
The model considers an instantaneous moisture loss once it has reached the crack. In
this way, water losses are computed when moisture escape out of the surrounding
continuum elements. This simplification is a good approximation in the cases in which
the crack pattern has a direct connectivity with the drying surface, but special attention
must be paid in the cases where there is no such connectivity. In this last case, moisture
loss should be zero. In any case, and as it will be shown in this chapter, the crack system
in the drying shrinkage simulations at the meso-level is mostly connected. This fact has
been confirmed also with recent experimental observations (Bisschop & van Mier,
2002a).

95
4.3.1. Influence of the crack depth on the drying process
In figure 4.10 the results in terms of RH distribution are presented for the case of
constant crack opening (50 microns). It is clearly observed that the effect of a larger
crack depth consists in increasingly facilitate the moisture migration out of the sample.
Note the semi-circular shape of the RH distribution in figure 4.9a, for the reference case
(zero crack opening), indicating that there is a sink of moisture at point P (shown in
figure 4.8a).
Figure 4.11 presents the results obtained in terms of the total weight losses (we)
evolution in time for the series of constant crack width (50 microns) and variable depth
(from zero, in the reference case, to 5cm, i.e. the total thickness of the sample). Figure
4.11a shows the results in absolute values, while the same results relative to the
reference case (subtracting the weight losses of reference case from each case) are
presented in figure 4.11b. A notorious effect on the rate of drying can be observed,
increasing considerably with crack penetration. It is also observed that, as one would
expect, all cases tend to the same asymptotic value, equal to the moisture total content
between 100 and 50% RH.

Figure 4.10. RH distribution for a drying period of 317.8days for the cases of constant
crack width (50 microns) and different crack depths of: (a) 0.0cm (reference case), (b)
0.25cm, (c) 0.5cm, (d) 1.0cm, (e) 1.5cm, (f) 2.0cm and (g) 5.0cm. Note that in the
limiting case (g) the crack is extended to the total thickness of the sample.

96
Figure 4.11. Evolution of total weight losses (we[g]) for the series of simulations
corresponding to a crack width of 50 microns and different crack depths (from zero, in
the reference case, to 5cm, i.e. the total thickness of the sample): (a) absolute values and
(b) values resulting from the difference between each case and the reference case.

4.3.2. Influence of the crack opening on the drying process

The same conclusions as in the previous subsection can be drawn from figures 4.12
and 4.12, in which the same results are presented for the series of simulations with
constant crack depth (2cm) and variable crack width (from zero, in the reference case, to
1mm). The effect of a higher crack width consists also in increasingly facilitate the
moisture escape. In this case, it can be observed in both figures that for crack widths of
approximately 50 microns or more, the increase in the drying rate due to an increase in
the crack opening is negligible, and that there is a notorious increase of drying rate when
passing from 1µm to 5µm of crack opening.

97
Figure 4.12. RH distribution for a drying period of 317.8days for the cases of constant
crack depth (2cm) and different crack openings of: (a) 0.0µm (reference case), (b) 1µm,
(c) 5µm, (d) 10µm, (e) 50µm and (f) 1,000µm (1mm).

Figure 4.13. Evolution of total weight losses (we[g]) for the series of simulations
corresponding to a crack depth of 2cm and different crack widths (from zero, in the
reference case, to 1mm): (a) absolute values and (b) values resulting from the difference
between each case and the reference case.

4.4. Coupled hygro-mechanical (HM) analysis at the meso-scale

In this section, the first coupled HM calculations at the meso-scale are presented.
More specifically, the simulation of drying shrinkage in concrete specimens is

98
addressed. A staggered approach has been chosen to perform the coupled analysis, as
discussed in section 4.2. Some verification examples of the mechanical behavior and the
moisture distribution are described, and the effect of coupling is evaluated.
In order to evaluate the coupled behavior of concrete specimens when subjected to
drying, a series of analyses has been performed over the same mesostructural FE mesh.
It consists of a 14x14cm2 concrete specimen with an arrangement of 6x6 aggregates
(volume fraction of 28%). As initial condition, RH = 1 throughout the domain is
assumed (full saturation is assumed). At t0= 28 days (not to be confused with the time at
which the mechanical properties are referred to, defined in Chapter 2), boundary
conditions become RH = 0.5 on the left and right edges, and no moisture flow is allowed
through the top and bottom faces. The specimen is assumed to be simply supported, so
that deformation can occur freely. The parameters for the matrix behavior in the
diffusion model have been adopted from the work by Roncero (1999), where they have
been fitted to experimental results in mortar specimens. Reasonable values have been
considered for the mechanical parameters, in order to represent a normal strength
concrete (compressive strength of around 40MPa). Calculations have been repeated for
uncoupled and coupled cases, in order to assess the coupling effect. In the former case,
two independent calculations are performed, the first one being the diffusion analysis
(from which the volumetric strain field entering the mechanical problem is obtained).
Additionally, different constitutive laws have been evaluated in both cases. Elastic and
viscoelastic (with aging) behavior have been compared, as well as the effect of
considering aging phenomena in the zero-thickness interface elements, as described in
Chapter 2.
Figure 4.14 presents the results in terms of RH distribution and total weight losses at
three different drying durations for the following cases: uncoupled (a), coupled with
non-aging joint elements (b), and coupled with interface elements considering the aging
effect (c). In the three cases, linear viscoelasticity with aging has been considered for the
matrix behavior. The aggregates are assumed to behave linear elastically and to remain
fully saturated throughout the simulation (the diffusion coefficient of natural aggregates
is usually negligible as compared to the matrix one). Results show that the effect of
coupling is a higher degree of drying, quantified by the weight losses and the RH
distributions. Accordingly, in the coupled calculations it can be observed that drying is
slightly more pronounced in the case where non-aging interfaces are considered as
compared to the case of aging interface elements. This is due to the fact that more
microcracking is expected to occur in the first case. The RH profiles at different ages for
the uncoupled case, calculated at the middle cross-section of the mesh, can be observed
in figure 4.15 (the discontinuities correspond to the presence of aggregates in the cross-
section).

99
Figure 4.14. RH distributions and corresponding total weight losses for a 14x14cm2
concrete specimen at different drying times (t-t0) of 20 days (left), 200 days (middle) and
10,000 days (right): uncoupled case (a), coupled with joints without aging effect (b), and
coupled with joints accounting for aging effect (c).

Figure 4.15. RH profiles calculated for the uncoupled case at different drying periods
(the discontinuities in each profile correspond to the presence of aggregates in the
middle cross-section analyzed).

The evolution in time of the total weight losses of the above three cases is shown in
figure 4.16. It is confirmed that the effect of coupling yields a higher degree of drying,
although it is not very pronounced in this particular case. In fact, the largest difference
between the coupled and uncoupled cases does not exceed 10% at any drying time.

100
Figure 4.16. Evolution in time of the total weight losses [g/cm] for the uncoupled case
(in red), coupled case with joints without aging (in green) and coupled case with aging
effect (in blue).

The results obtained from the mechanical analysis of the previous cases are presented
in figures 4.17 and 4.18. A sequence of the evolution in time of the energy spent in
fracture processes for the coupled case with aging viscoelasticity and aging interface
elements is depicted in figure 4.17. Initially, as it would be expected, cracks
perpendicular to the two drying surface develop. The red color indicates that the crack is
loaded. As the drying front penetrates into the specimen, the cracks left behind unload
(in blue). Cracks are also observed, at an advanced drying stage, in the interior of the
sample, on the aggregate-matrix interfaces but most importantly at perpendicular
directions to the aggregates. This is due to the restraining effect caused by embedded
(more rigid) particles in a shrinking viscoelastic matrix. Aggregates have a considerably
higher elastic modulus than the matrix one (stiffness of the aggregates is of the order of
3 times that of the matrix). In fact, these results qualitatively agree with recent
experimental observations (Bisschop & van Mier, 2002a,b), as already discussed in
Chapter 3 (see figure 3.15 in section 3.1.5.2.). These findings have been the motivation
of a subsequent deeper study on the effect of aggregates on the drying-induced
microcracking, which is presented in the next section.
In order to evaluate the effect of considering an aging viscoelastic behavior of the
material, another simulation has been performed with a linear elastic constitutive law for
the matrix, with the mechanical properties at an age of 28 days. It should be noticed that
in the linear elastic case, non-aging interface elements have been considered, whereas
the viscoelastic analysis includes joint elements with aging effect. In fact, the elastic
case represents a simplification in comparison with the viscoelastic behavior, since in
the latter the material is capable of internally redistribute stresses and transfer them from
the matrix towards the more rigid inclusions. In figure 4.18 the deformed meshes of both
coupled cases at 2,000 days of drying (scale factor of 150) are compared. It can be
immediately observed that cracking is more pronounced in the linear elastic matrix case,
due to the fact that in this situation there is no aging effect (i.e. strength and elastic
modulus do not increase with time) and stress relaxation is prevented. On the contrary,
these two effects play an important role in the simulation considering a viscoelastic
matrix behavior. The comparison between elastic and viscoelastic behavior is presented
in more detail in section 4.5.2.5.

101
 It should be noted that in the proposed model the aging effect is decoupled from the
moisture diffusion analysis. In fact, it is known that aging occur only as long as moisture
conditions are close to fully saturation (with a threshold of around 80% RH, according to
Bazant & Najjar, 1972). For lower RH, the aging of the material stops. At present, this is
not accounted for in the model and its introduction would certainly represent an
improvement in the model formulation.

Figure 4.17. Dissipated fracture energy in the interface elements for 20, 200, 2,000 and
10,000 days (from left to right), for the coupled case with aging viscoelasticity for the
matrix, and aging interface elements.

Figure 4.18. Deformed meshes at 2,000 days of drying (scale factor of 150): coupled
cases with (a) viscoelastic matrix with aging (and aging interface elements) and (b),
linear elastic matrix with non-aging joint elements. The blue box corresponds to the
dimensions of the undeformed initial mesh. Note the barrel shape of the obtained
deformed meshes.
The previous simulations presented show a barrel type deformed configuration,
meaning that the cross-sections do not remain planar after deformation. A way to
remedy this is to use more slender specimens and measure strains at a distance from the
specimen’s ends, large enough so as to avoid any effect of the edges, which is actually
usual practice in drying shrinkage experiments. From a numerical viewpoint, another
option would be to use the same numerical specimen as before but forcing the end faces
to remain planar throughout the drying process. In this way, the central part of a larger
specimen can be simulated, thus saving computer time. This is shown in figure 4.19, in
which the results of the simulation of the two coupled cases just described are presented
in terms of evolution of longitudinal strains in time. The meshes used in the simulations

102
have dimensions of 15x15cm2 and 15x45cm2 both with 27,6% of aggregate volume
content (rounded shape). Drying is identical in both cases (RH = 0.5 on the lateral faces
and no moisture flow is allowed on top and bottom). For the mechanical analysis, simple
support conditions are considered for the larger specimen (as in the previous examples).
On the other hand, the top and bottom faces of the smaller specimen are forced to remain
planar (materialized on the top face with the addition of a rigid plate and an interface
between them with no resistance to slip, thus only avoiding penetration). In this way, the
central part of a more slender specimen is simulated. It can be seen that, as expected, the
averaged longitudinal strain evolution at the edges is similar in both cases, whereas the
evolution in the central part of the larger specimen is slightly lower (although with the
same asymptotic value). Thus, the validity of the simplified case is confirmed.

Figure 4.19. (a) Evolution of longitudinal strains [mm/m] with time [days] for the
15x15cm2 and 15x45cm2 meshes: average of the strains measured at the left and right
edges, and strain at the center of the 15x45cm2 case. (b) Meshes used in the simulations.

Figure 4.20 depicts a sequence of the drying process for the central part (20cm
length) of the slender specimen (15x45cm2). Moisture distributions together with the
corresponding crack patterns are shown for different drying times, in order to give a
better idea of the coupled behavior. The model predicts that surface microcracks start
unloading before 200 days drying. Also notice that only at an advanced stage of drying
the restraining effect of the aggregates is noticeable. In figure 4.21, the RH profiles are
shown together with the corresponding stress profiles in the same cross section for
different drying periods. The selected cross section corresponds to the upper face (in
contact with the rigid platen) of the 15x15cm2 mesh, since in this particular case the
stress profile can be extracted as the contact normal stresses between the platen and the
specimen (elastic interface elements were inserted in between with very large normal
stiffness and a negligible value for the tangential one). It can be seen in these profiles
that the elevated RH gradient at the beginning of drying (already at 2 days) translates
into important tensile stresses near the drying surface (and compression ones in the
interior). As the drying process evolves with time, the tensile state moves into the
specimen, and the outer layers unload (starting at around 130 days). It is worthy noting
the four interior compression intervals obtained in the simulation at advanced states of
drying. This is due to the presence of aggregates at these same depths and at a short
distance from the cross section (see figure 4.19b), which are known to be under
compression stresses when immersed in a shrinking matrix (see next section).

103
Figure 4.20. Sequence of the drying process for the central part (20cm) of the 15x45cm2
mesh: RH distribution and corresponding crack patterns (represented by the energy spent
in crack formation) for different drying times (16, 200, 500, 3,850, 7,000 and 10,000
days).

104
Figure 4.21. Sequence of the drying process for the case of the 15x45cm2 mesh: RH
profiles and corresponding stress profiles for different drying times (2, 6.5, 21, 130,
2,000, 5,000 and 10,000 days).

105
 Useful information can be obtained by plotting the longitudinal strains against total
weight losses. In this way, results from the diffusion and mechanical analyses are
contrasted at each time step, which gives us a better insight to the coupled problem.
These results are presented in figure 4.22 for the 15x15cm2 mesh (forcing the end faces
to remain planar), being the longitudinal strain the average of the measures at the left
and right edges. Some degree of nonlinearity can be observed, which at a first glance
could be unexpected due to the initially adopted linear relationship between volumetric
strain and weight loss at the local level (i.e. at the constitutive level, for a constant
shrinkage coefficient). But in fact, these results qualitatively agree with similar curves
that were experimentally determined by Granger (1996) and analyzed in Granger et al.
(1997a,b). A more detailed analysis of the results shows that three well-defined phases
may be differentiated. In the first phase, shrinkage strains evolve much slower than total
weight losses, which in our calculations correspond to the drying-induced microcracking
effect at the specimen skin. These microcracks relax the compressive strains, since a
crack opening displacement results in a tensile deformation, contrary to the measured
shrinkage longitudinal strains. The second phase of the curve corresponds to a fairly
linear relationship between strains and weight losses, which could be interpreted as the
adopted (linear) intrinsic relationship of the material. Finally, in the last portion of the
curve, a second kink is present that correspond to an advanced drying stage in which
relatively large weight losses produce little deformations (the specimen is reaching
equilibrium with the environment).

Figure 4.22. Longitudinal strain [mm/m] vs. total weight loss (in %) for the 15x15cm2
mesh, in which the top and bottom faces remain planar throughout the simulation.

Finally, another interesting way of analyzing the results in terms of RH distribution

consists of plotting the evolution in time of drying for fixed points in the FE mesh. This
allows observing the RH gradient with respect to time (or drying rate). This is presented
in figure 4.23, for fixed points in the mesh at different depths (distance to the drying
surface) and considering the effect of cracks on the drying process. These results
correspond to a coupled drying shrinkage simulation with aging of a 15x30cm2 mesh
(obtained at the middle cross-section) with 22% aggregate content (6x12 arrangement).
It may be observed that the RH gradient decreases with the increase of the distance to
the drying surface. Moreover, for two points situated at the same depth (in the figure
0.3cm and 1.3cm), drying is more intense in the case in which the point lies on the crack
wall.

106
Figure 4.23. Evolution of internal RH in time of exposure to drying, corresponding to
mesh points at different depths (distance to the drying surface): influence of cracks on
the drying rate for depths of 0.3cm and 0.7cm.

4.5. Effect of the aggregates on the drying-induced microcracking

Homogeneous materials, as hardened cement paste, subjected to drying shrinkage,
typically show cracks perpendicular to the drying surface, with limited penetration, as
shown in the previous section and discussed in Chapter 3. When adding inclusions in a
shrinking matrix, as in the case of concrete, shrinkage strains are restrained by their
presence, thus inducing a stress state in the system like the one presented in figure 4.24.
The stress level is very much dependent on the stiffness ratio between aggregates and
matrix, being higher for large values of this ratio. If the induced stresses exceed the
tensile strength of the matrix, cracks are generated in radial as well as tangential
direction to the aggregates (Goltermann, 1995; Hearn, 1999; Hsu, 1963; Koufopoulos &
Theocaris, 1969). These cracks may be important from the durability point of view,
since they can significantly influence the mechanical and transport properties of the
material (Yurtdas et al., 2005). Recent experimental results (Bisschop & van Mier,
2002b, Bisschop, 2002) have shown the importance of the effect of inclusions (stiffer
than the matrix), regarding both the quantity and their size, in the drying-induced crack
patterns in cementitious systems with glass spheres. To this end, they prepared a
cementitious composite with mono-sized glass spheres as inclusions (and a HCP matrix),
as well as concrete specimens with irregular aggregates of controlled-size, from which
they extracted similar conclusions (Bisschop, 2002).
The effect of the inclusions on the drying-induced microcracking has been
numerically studied with the present model. First, some simulations with only one
central inclusion have been performed in order to address the effect of the size of the
inclusion. Thereafter, an extended series of calculations has been carried out to study the
influence of volume fraction and size of aggregates in concrete specimens. In this
section, the most relevant results are presented, together with a discussion of the main
observations.

107
Figure 4.24. Stress distribution in the matrix and the inclusion/s obtained by photo-
elasticity, for the case of (a) a single inclusion, and (b) a squared array of inclusions
(from Koufopoulos & Theocaris, 1969).

4.5.1. Microcracking around one single inclusion

Prior to the study of the effect of polygonal aggregates (generated through Voronoï-
Delaunay theory) on the drying-induced microcracking, which is presented in the next
section, a series of simulations have been carried out in order to evaluate the effect of a
single circular inclusion on the matrix crack pattern. To this end, the program DRAC has
been used for the mechanical analysis, in which a controlled shrinkage profile has been
assumed (local volumetric strains are artificially generated so that homogeneous
contraction of the matrix, as well as a constant gradient of shrinkage strains through the
specimen thickness, are obtained). No diffusion analysis is performed in these cases.
Unstructured meshes with circular inclusion/s have been generated with the program
PARSIFAL (Particle Simulation for Analysis), developed by Saouma and coworkers
from the University of Colorado at Boulder (Puatatsananon et al., 2008). Zero-thickness
interface elements have been added a posteriori in all matrix-matrix and inclusion-
matrix contacts (inclusions are not allowed to crack). Figure 4.25 shows the meshes used
in the single-inclusion simulation series and their dimensions. The same series has been
repeated with the same material parameters but in which the size of the meshes is
reduced by one order of magnitude, in order to analyze 2mm, 4mm and 6mm in diameter
inclusions.

Figure 4.25. Meshes used in the simulations with a single circular inclusion with a
diameter of 2, 4 and 6cm.

108
 The simulations entail inducing a controlled homogeneous volumetric strain
throughout the matrix so that the resulting crack pattern of the composite and its
dependence on inclusion size can be evaluated. Mechanical properties are similar to
those used for the simulations of the preceding section, although a linear elastic matrix
with non-aging interface elements have been adopted in order to highlight the obtained
crack patterns (Ematrix = 25,000MPa; Einclusion = 72,000MPa; i.e. a ratio of almost 3).
Figure 4.26 shows the principal stresses distribution in the matrix and the inclusion for
the case of a 4cm inclusion. It can be observed that the aggregate is always subjected to
compression, whereas the matrix is subjected to circumferential tensile stresses and a
compression radial state, which agrees with theoretical considerations (Goltermann,
1994).

Figure 4.26. Vectorial distribution of principal stresses for the case of an inclusion with a
diameter of 4cm: (a) σI, (b) σII. σI varies between -10 and 15MPa and σII between -7,7
and 0MPa). The circumferential tensile state in the matrix is readily observed (in red).
The case corresponds to a homogeneous shrinkage strain of the matrix (εεvol = 0,75mm/m
at all points in the matrix).

In figures 4.27 and 4.28, the results are presented for the two series of simulated
meshes (i.e. 2, 4 and 6mm in figure 4.27, and 20, 40 and 60mm in figure 4.28) in terms
of the energy dissipated in fracture processes, over the deformed configurations. In all
the cases, cracks with perpendicular directions, as well as tangential, to the aggregates
are observed. It may be noted that, for the same level of contraction of the matrix,
microcracking is more pronounced as the inclusion size increases, for the whole range of
sizes, i.e. for diameters varying between 2 and 60mm (note the different scales used for
the dissipated energy, with respect to GFI, between the two series). Maximum crack
openings in the first series do not reach 0.5µm, whereas in the second series they vary
between 1 and 10µm. These results agree well with experimental observations and
analytical results found in the literature (Bisschop & van Mier, 2002b; Goltermann,
1995). Some of the simulations have been repeated with a constant gradient of shrinkage
strains through the specimen thickness (minimum at the top face and maximum at the
bottom one), which is a closer situation to reality, obtaining similar microcrack patterns.

109
Figure 4.27. Energy dissipated in fracture processes (Wcr maximum of 0.024.GFI) for the
cases of 2 (a), 4 (b) and 6mm (c) inclusions. Results over deformed configuration (scale
factor of 150).

Figure 4.28. Energy dissipated in fracture processes (Wcr maximum of 0.534.GFI) for the
cases of 20 (a), 40 (b) and 60mm (c) inclusions. Results over deformed configuration
(scale factor of 150).

4.5.2. Microcracking in concrete specimens with multiple inclusions

Recent experimental observations highlight the influence of the size and quantity of
the aggregates in the microcrack-pattern of concrete subjected to drying (Bisschop &
van Mier, 2002b). When non-shrinking inclusions are inserted in a shrinking matrix,
microcracks may appear in radial as well as tangential directions to the aggregates,
depending on the size and number of inclusions, and the behavior of the matrix.
In order to evaluate the ability of the HM mesostructural model to capture the effect
of the aggregate distribution, number and size on the drying-induced microcracking, a
significant number of mesostructural meshes of 120x40mm2 has been generated and
numerically tested. Polygonal aggregates inscribed in circumferences (of constant
diameter in each case, see figure 4.29) have been adopted, in the spirit of recent
experimental tests performed by van Mier and coworkers (Bisschop & van Mier, 2002b),
in which different sizes (for each case) of mono-sized glass spheres have been used as
aggregates (see figure 3.15 in Chapter 3 for an example), with three different volume
fractions. First, aggregates of 2mm, 4mm and 6mm in diameter have been used, as in the
experiments (Bisschop, 2002), with three different volume fractions of 20%, 30% and
40% (different volume fractions to the experimental ones, which were of 10, 21 and
35%, were used because a mesh with 10% volume fraction would not be realistic with
the present mesh generation procedure), yielding a total of nine cases. Additionally, in
some of these cases (namely the series of 4 and 6mm with 20, 30 and 40% volume
fraction) four different meshes have been generated and simulated in order to study the
effect of random generation of geometries.

110
 In a second step, the same series of simulations have been repeated for larger
specimens (300x100mm2) with polygonal aggregates of 10, 15 and 20mm and the same
volume fractions (20, 30 and 40%). In fact, the same meshes as the 120x40mm2 series
have been used by adopting a scaling factor of 2.5.
The geometry and boundary conditions are shown in figure 4.29 for the diffusion and
mechanical analyses. Samples are assumed to be initially fully saturated (RH = 1). At
the age of 28 days a RH of 0.3 (as in the experiments) is imposed on the top face of the
specimen, and no moisture flux is allowed on the remaining ones. For the mechanical
analysis a simply supported beam case is considered. Aging viscoelasticity is assumed
for the matrix, with aging interface elements. Additionally, some simulations have been
performed assuming linear elastic matrix behavior with non-aging interface elements, in
order to assess the effect of creep on drying-induced microcracking. Similarly to the
experimental campaign (Bisschop & van Mier, 2002b), only the central part of the
specimen is used for quantification of the crack patterns (i.e. an area of 40x40mm2), so
that the boundary effects are eliminated. In order to be able to compare simulations
made with different matrix volume fractions, and thus different initial contents of
evaporable water, the degree of drying is defined as the weight loss at time t with respect
to the initial (before drying starts) water content (Bisschop & van Mier, 2002b).

Figure 4.29. Geometry and boundary conditions adopted for the simulations (mechanical
and diffusion analysis), for the case with 30% volume fraction (4mm aggregates).

With the purpose of quantifying the crack patterns and determine the preferential
orientation of the microcracks, rosettes (polar diagrams) have been constructed as a post-
process for all the simulated cases, in the spirit of (Sicard et al., 1992; Bisschop & van
Mier, 2002b), as shown in figures 4.30 and 4.31. Polar diagrams are constructed by
adding up the length of all (opened) crack-segments of the matrix with the same
orientation (the angular space has been divided in six groups, every π/6 radians). The
aggregate-matrix cracked interfaces have not been considered in the quantification, to
make results more comparable to the experimental ones, in which these were also
neglected. In the cases of 4mm and 6mm aggregates (for all the volume fractions), the
rosettes have been computed as the average of four meshes (i.e. 4 different simulations)
with random aggregate distribution, but keeping at a constant value both the size and the
volume fraction.
Polar diagrams help us determining the mechanisms acting on drying shrinkage
microcracking (Bisschop, 2002). As stated in previous paragraphs, a homogeneous
material will mostly produce cracks perpendicular to the drying surface when subjected
to a low RH environment, thus resulting in rosettes with elongated shapes. On the other
hand, when cracking due to the effect of aggregate restraining is predominant, cracks
will radiate from the inclusions, thus producing (at least in theory) cracks along all

111
possible directions, with approximately the same probability. In this case, it is expected
that the rosettes have a rounded shape, due to this random orientation of cracks. Finally,
an intermediate case, in which both mechanisms are present (and which is generally the
case), will show rosettes with a shape varying between these two limiting cases
(Bisschop & van Mier, 2002b). Figure 4.30 shows the results obtained by Bisschop and
van Mier on the experimental drying shrinkage tests. The influence of the aggregate
volume fraction at constant aggregate size and of the aggregate size at constant volume
fraction is evaluated for cementitious composites with glass spheres. These two cases
clearly show that an increase in any of these two variables yield an increase of the
degree of microcracking and of the aggregate restraining effect (although this was not
the case in all of the experimental series, see Bisschop, 2002).

Figure 4.30. Experimental results of drying shrinkage microcracking in cementitious

composites with glass spheres: crack patterns of the central part of the samples
(superposed maps of 4 cross-sections) and corresponding rosettes. (a) Effect of the
aggregate size (0.5, 1, 2 & 4mm) at constant volume fraction (35%), and (b) effect of the
aggregate volume fraction (0, 10, 21 & 35%) at constant size of 4mm (from Bisschop &
van Mier, 2002b).

Note that the construction of a polar diagram of the entire system of interface
elements in a given mesh, gives us information on the quality of the representation.
Accordingly, an appropriate mesh should have a fairly equal distribution of interface
elements (regarded as potential cracks) in each direction. On the other hand, a more
elongated rosette (in any direction) means that the mesh has a preferential cracking
direction, thus loosing objectivity of the results. An example of the initial distribution of
interface elements (excluding the aggregate-matrix interfaces) for the case of 4mm
aggregates with 20%(+/- 1%) volume fraction is given in figure 4.31, together with the

112
effect of randomness in the geometry generation of 4 different meshes of the same case.
In this case, the axis represents the percentage of all interfaces in a specific direction.
Since the angular space has been divided in six groups, a perfect distribution would have
around 16% of the interface elements in each direction. The average rosette has also
been calculated and is shown in figure 4.31e. It may be seen that the initial distribution is
acceptable in all cases, and that considering the average of four different cases yields a
much more regular distribution.

Figure 4.31. Initial rosettes or polar diagrams for four different meshes with 4mm
aggregates and 20% volume, and average rosette of the same four cases, showing a fairly
rounded shape in all cases (values in %).

4.5.2.1.Effect of the degree of drying

Figure 4.32 presents some typical results of the evolution of drying-induced
microcracking, for the cases of 6mm and 2mm aggregates with 40% volume fraction, in
terms of the energy spent in fracture processes vs. percentage of weight loss (10 and
30% of the initial water content). It can be seen in both cases that at the beginning of
drying microcracks are mainly produced perpendicular to the drying surface, thus
presenting an elongated shape of the polar diagram. On the other hand, for an advance
stage of drying (30% in this example), the crack penetration front is more pronounced,
which is characterized by a more ‘rounded’ polar diagram, with microcracks radiating
from the aggregates. Most of the simulations performed in this work show the same
trend.

113
Figure 4.32. Effect of the degree of drying on the drying-induced microcracking, for the
cases of 6mm (top) and 2mm (bottom) aggregates with 40% volume fraction: crack
patterns for 10% and 30% drying (of the central 40x40mm2 part) and corresponding
rosettes (10% in white and 30% in brown). The average rosette of 4 random geometries
is shown for the case of 6mm aggregates.

4.5.2.2.Effect of aggregate volume fraction

Figure 4.33 presents the results regarding the microcrack maps of the drying
shrinkage simulations for 3 different volume fractions (20, 30 and 40%) in the case of
4mm aggregates (Idiart et al., 2007). The upper row shows the microcrack patterns,
represented by the fracture energy, obtained at 30% drying for the 3 cases. No clear
trend of the variation of the microcracking depth of penetration with increasing volume
fraction is observed, which is in agreement with experimental observations (figure
4.30b). In the lower row, the average polar diagrams corresponding to the same drying
situation for each case are presented (average of 4 different calculations for each volume
fraction). It can be clearly seen that an increase in volume fraction results in much less
elongated polar diagrams (note also the increase of the aspect ratio, calculated as the
relation between horizontal and vertical axis of the polar diagram). This means that the
aggregate restraining effect (causing radial cracks to the inclusions in any direction)
becomes more important when a larger number of inclusions are present. For the lowest
volume fraction, an elongated diagram is found, suggesting that aggregate restraining is
not very important in this case: microcracking occurs mainly perpendicular to the drying
surface, due to the fact that the inner layers are drying much slower than the outer ones.
These results are in agreement with the experimental findings by Bisschop & van Mier
(2002b).
The results obtained for the 2mm aggregates are presented in figure 4.34, for the three
volume fractions (20, 30 & 40%) and 3 different degrees of drying (10, 20 & 30%). As
in the case of 4mm aggregates, crack patterns and polar diagrams are included, although
in this case only one simulation for each volume fraction has been performed (due to the
computational cost). Overall, the same trends as in the previous example are observed.

114
Figure 4.33. Effect of aggregate volume fraction (20, 30 & 40%) in drying shrinkage
microcracking at 30% drying, for the case of 4mm inclusions: (upper row) microcrack-
maps of mid-specimens (40x40mm2) and (lower row) average rosettes and aspect ratios.

Figure 4.34. Effect of aggregate volume fraction (20, 30 and 40%) in drying shrinkage
microcracking at different degrees of drying (10, 20 and 30%), for the series of 2mm
inclusions: microcrack-maps of mid-specimens (40x40mm2) and corresponding rosettes
(rightmost). Note the small difference between rosettes for 10 and 20% drying.

The aspect ratios at 30% drying are 0.16, 0.205 and 0.333 for 20, 30 and 40% volume
fraction, respectively, showing again that the influence of volume fraction on drying-

115
induced microcracking is an increase of the aggregate restraining effect.
4.5.2.3.Effect of aggregate size
The second studied variable, apart from the volume fraction, is the effect of the size
of the aggregates at constant volume fraction. Results by Bisschop and van Mier are
shown in figure 4.30a. It can be seen that slightly more rounded rosettes are found with
increasing aggregate size, but also considerably higher crack lengths are observed for the
cases with large aggregates as compared to the small ones. The results obtained with the
present model from the study of the aggregate size influence are shown in figure 4.35, in
terms of the polar diagrams and corresponding aspect ratios. At a first glance, the
comparison does not look very encouraging (although considerations explained later
change totally this first impression). Although the aspect ratio and the micro-cracks
depth shows the correct trend (an increase in the aspect ratio and the microcrack depth
with aggregate size indicates a more important aggregate restraining effect), the rosettes
expressed in cm (figure 4.35a) show decreasing values for increasing size, contrary to
what has been experimentally observed (figure 4.30a).

Figure 4.35. Effect of aggregate size for a constant volume fraction of 20% (2, 4 and
6mm aggregates), at 30% degree of drying: polar diagrams and aspect ratios for each
case. Values are represented in cm (a) and also as a percentage (%) of the total interface
elements in each direction (b). The cases of 4mm and 6mm correspond to the average of
4 different geometries. (c) Microcrack maps for each case.

116
 After a first study of the results, this lack of agreement with experimental
observations could be attributed to a number of factors. First of all, the fact of
performing 2D calculations is a simplification by which the in-plane cracks due to out-
of-plane aggregates in real 3D cases are not considered. Also, as discussed in Chapter 2
(section 2.2.4), the adoption of a correct 2D geometry from a 3D composite is not an
easy task, and in this work it has been decided to simulate mono-sized aggregates in 2D
sections, whereas in the real specimens, generally none of the cross-sections fulfill this
hypothesis. Finally, the structured discretization of the mesh generation procedure forces
larger total interface element length with decreasing aggregate size (at constant volume
fraction). This means that meshes with smaller aggregates have more potential crack
paths available. If the same results are normalized with the total interface element length
in each direction, the rosettes presented in figure 4.35b are obtained. Note that in this
case the correct trend is observed, showing increasing rosettes, both in magnitude and
aspect ratio, for increasing aggregate size.
But a better interpretation of the numerical results, may be obtained by filtering (i.e.
eliminating from the plotted results) microcracks with crack openings (or, equivalently,
spent fracture energies) that are below a certain threshold. In fact, in Bisschop (2002) it
is underlined that microcracks with openings much smaller than 1 micron are not always
captured by the experimental technique (although the exact threshold is not clear). Thus,
it seems reasonable to consider a threshold of e.g. 0.1 microns (corresponding to 0.3% of
GFI with the present parameters). Results obtained by applying such a filter as a post
process to the previous cases are shown in figure 4.36, in terms of the polar diagrams
and crack patterns. It can be seen that the trend is now totally changed even in terms of
absolute crack lengths, and that larger rosettes are obtained by increasing the aggregate
size, which is in agreement with the experimental observations.

Figure 4.36. Effect of aggregate size (2, 4 and 6mm aggregates) for 20 % volume
fraction, at 30% degree of drying: results obtained by filtering the crack patterns with a
lower bound for the fracture energy of 0.3% of GFI. (a) Polar diagrams and aspect ratios
for each case (values represented in cm). (b) Microcrack maps for each case.

117
 Note also that the difference of the crack patterns is more pronounced, showing
deeper crack fronts with increasing aggregate size. It is concluded that, although the
polar diagrams seem to be very sensitive to the filtering of cracks with very small
widths, the model is indeed capable of capturing the aggregate size effect on drying-
induced microcracking.
4.5.2.4. Randomness effect
As already mentioned, the series of 4 and 6mm aggregates, with 20, 30 and 40%
volume fraction, have been used to study the randomness effect on the microcrack
patterns by generating four different geometries with the same parameters. In figure
4.37, the results of the 4mm aggregates with 40% volume fraction are presented for the
four simulations in terms of polar diagrams (left, including the average rosette) and
fracture energy of the mid-part of the specimen. Good agreement between the results for
the different geometries is found: crack patterns show similar crack penetration fronts
and the polar diagrams present a reasonable variation. Only one of the cases presents an
aspect ratio with a 40% deviation from the average value of 0.449. Similar results have
been obtained for the rest of the studied cases (i.e. 4mm aggregates with 20 and 30%
volume fraction, and 6mm with 20 and 30% volume fraction), with the exception of the
case with 6mm aggregates and 40% volume fraction, for which more important
variations for each case have been observed (one of the meshes has shown 74% variation
from the average behavior).

Figure 4.37. Randomness effect of four generated geometries, for the case of 4mm
aggregate with 40% volume fraction, on drying shrinkage microcracking at 30% degree
of drying: polar diagrams for each case and average rosette (left), and microcrack-maps
of mid-specimens (40x40mm2) for the same cases (right).

4.5.2.5. Effect of creep

An interesting factor is the behavior of the matrix and its influence in microcracking.
The effect of creep on the drying-induced microcracking has been addressed by
considering two cases with aging viscoelasticity (and aging interface elements) and

118
elasticity (with non-aging interfaces) for the matrix, but otherwise the same material
parameters (López et al., 2007), similar to the cases in section 4.4.1. The case with 4mm
aggregates with 30% volume fraction has been selected. As previously stated, the effect
of considering linear elasticity for the matrix is mainly to prevent any stress
redistribution. Moreover, since there is no aging effect, microcracking is expected to be
slightly higher at advanced drying states. Results of the two simulations are presented in
figure 4.38 for two drying states (10 and 20%). It can clearly be observed that the effect
of creep consists of considerably decreasing the penetration depth of the cracking front,
and thus the microcracking due to the aggregate restraining effect. This is also supported
by the obtained rosettes, being more elongated for the case of aging viscoelasticity, as
well as presenting a much more reduced total crack length.

Figure 4.38. Effect of creep on the drying-induced microcracking for the case of 4mm
aggregates with 30% volume fraction for 10 and 20% drying: crack pattern (quantified
by fracture energy) of the central part (40x40mm2) and corresponding rosettes for aging
viscoelastic matrix with aging interface elements (top), and elastic matrix with non-
aging interfaces (bottom).

4.5.2.6. Effect of aggregate shape

Finally, a number of simulations considering circular aggregates have been
performed, in order to assess the effect of the shape of aggregates. One of the meshes
generated for this purpose is shown in figure 4.39a, together with the imposed boundary
conditions (which were similar to the previous cases). The software PARSIFAL
(Puatatsananon et al., 2008) has been used to generate the meshes, with the addition a
posteriori of the zero-thickness interface elements. Mono-size circular aggregates have
been generated with a diameter of 6mm (aggregate volume fraction equal to 27% of the
total volume). Aging interface elements are inserted along all the contacts between
matrix continuum elements and around all the aggregates, with different mechanical
properties for the latter (no interfaces within the aggregates), in order to allow the most
relevant failure mechanisms to develop. Note that the matrix polar diagram in figure
4.39b has a very reasonable round shape, indicating an appropriate initial distribution of
the total matrix interface elements. Linear visco-elasticity with aging is adopted for the

119
matrix while aggregates are considered to remain linear elastic. Material parameters
were kept the same as in the previous examples.
The main advantage in this case is that the adoption of circular inclusions makes the
simulations more comparable to the experiments (Bisschop & van Mier, 2002b).
Nonetheless, it is not the objective of this work to repeat the entire series of simulations,
but to show the potential of the present coupled model with different mesostructural
representations. In figure 4.40, results of a wider study are presented, in terms of the
fracture energy spent at 30% drying (with respect to the initial water content). Typical
microcracks perpendicular to the drying surface can be observed (maximum crack width
of around 18 microns), but also radial and circumferential microcracks around the
aggregates appear as a result of drying (in this case, only microcracks with a crack width
larger than 0.5 microns have been included, in order to obtain more clear crack patterns).
In figure 4.41, the fracture energy spent at 5% (a) and 30% drying (b) is presented for
the central third of the specimen, and the corresponding rosettes are shown (5% drying
in white; 30% drying in grey; aggregate-matrix interface cracks have not been included).
At 5% drying, the rosette has an elongated shape in the vertical direction, meaning that
at the beginning of the drying process microcracks occur perpendicular to the drying
surface. At 30% drying, aggregate restraining is the predominant effect causing
microcracking in all directions, since cracks develop radial to the inclusions. Thus, the
polar diagram obtained shows the same trends as with the polygonal shape aggregates.

Figure 4.39. (a) Mesh (120x40mm2) and boundary conditions used in the simulation with
6mm aggregates, and (b) initial rosette for the matrix joint elements distribution.

Figure 4.40. Energy spent in fracture processes (red = loaded joints; blue = unloaded
joints), at a drying state corresponding to 30% of the total initial water content.

120
Figure 4.41. Energy spent in fracture processes (red = loaded; blue = unloaded joints), at
a drying state corresponding to 5% (a) and 30% (b), and corresponding rosettes (c).

4.6. Simulation of the experiments by Granger (1996)

4.6.1. Description of the tests
This section presents the main results of the analysis of the drying shrinkage
experiments by Granger (1996). In that work, the drying of cylindrical specimens (16cm
in diameter and 100cm height) of different concrete types was studied. All specimens
were subjected to a 50% (+/- 5%) RH environment, allowing a radial drying (top and
bottom faces were sealed). In this thesis, a normal strength concrete has been chosen for
the simulations (the so-called Penly concrete in the series by Granger). Details of the test
and material properties are summarized in table 4.1.

Penly concrete - concrete mixture details & summary of experimental setup

Specimens size φ 16cm, 100cm high Density (kg/m3) 2270
Aggregate Silico-limestone E 28d (GPa) 36,2
3
G: 12.5/25 (kg/m ) 682 fc 28d (MPa) 34,3
g: 5/12.5 (kg/m3) 330 ft 28d (MPa) 3,4
3
s: 0/5 (kg/m ) 702 Relative humidity 50% +/- 5%
3
Filler (kg/m ) 50 Temperature 20º +/- 1ºC
3
Cement (kg/m ) 350 Drying Radial
3
Water (kg/m ) 202 Measurement basis 50cm central
3
Plasticiser (kg/m ) 1.15 Curing period (days) 28
Table 4.1. Composition and mechanical properties at 28 days of the normal strength
concrete (Penly concrete), and summary of test conditions (Granger, 1996).

These experimental results are of particular interest in the sense that this study
includes the relationship between longitudinal strains and total weight loss, which is a
key feature in order to construct a coupled hygro-mechanical model. Most of the
experimental campaigns found in the literature focus their attention on either the strains
or the weight losses. This relationship allows us to link the mechanical behavior (via the
longitudinal strains) with the drying process (via weight losses), and is thus a measure of
the coupled behavior. In this sense, availability of RH profiles would have been optimal
in order to obtain a more robust fit, although this has been unfortunately not the case. In
spite of that, the experiments by Granger (1996) remain one of the most complete
experimental campaigns concerning drying shrinkage in concrete, which has been the

121
reason to adopt his experiments in this work, even if not exempt from specific
difficulties.
Perhaps the main of such difficulties is the cylindrical shape of the specimens, which
is a disadvantage for the present 2D mesostructural model, since an axisymmetric
analysis is not suitable for this kind of representations. Indeed, an axisymmetric analysis
in this case would imply two main inconsistencies. First, the symmetry axis should not
cut aggregates in order to avoid the generation of irregular shape of aggregates. But even
more important, in the axisymmetric representation, the aggregates would represent in
fact bodies of revolution, i.e. rigid ‘donuts’ within a more flexible matrix. This may
totally alter the analysis, yielding undesired consequences in the results. Nonetheless, the
lack of similar experimental results obtained on prismatic specimens has forced the
selection of these specimens. In order to eliminate to some extent the above mentioned
undesired features, it has been decided to perform a ‘semi-axisymmetric’ analysis for the
moisture diffusion simulation, coupled with a 2D plane stress analysis for the
mechanical problem, over the same FE mesh (which is a necessary condition in the
present model). In the case of moisture diffusion, the axis of revolution is positioned at
the center of the mesh (i.e. at D/2, with D being the diameter of the cylinder), and the
spin is limited to a π value (instead of 2π). In this way, an undesired behavior due to
‘donut shape’ of the aggregates in the mechanical simulations is avoided and the
cylindrical shape for the moisture diffusion problem, which is an essential characteristic,
is taken into account.
The FE mesh used in the simulations has the dimensions 16x50cm2. The height of the
mesh has been chosen to be equal to the strain measurement basis used in the
experiments (Granger et al., 1997b). The aggregate volume fraction is 26% and the
maximum and minimum sizes are 25.9mm and 10mm, respectively (the shape is
polygonal with aggregates inscribed in circumferences of variable diameter). In order to
simulate the central part of a larger specimen (the real specimens are 16x100cm2),
similar boundary conditions as in figure 4.19 (for the 15x15cm2 mesh) have been
adopted. Material parameters finally adopted for the best fit are summarized in table 4.2.

Material parameters adopted

Diffusion analysis (matrix) Mechanical analysis (continuum)
Initial humidity 100% E matrix aging Maxwell chain
2 -5
D0 (cm /day) 5x10 E aggr (MPa) 70000
2
D1 (cm /day) 2x10-1 ν matrix ; ν aggr 0.2 ; 0.2
β (−) 3.0 Mechanical analysis (interface elements)*
3
C (g. cem/cm mat) 0.473 χ (MPa) 2.0 ; 4.0
α hydration 0.90 c (MPa) 7.0 ; 14.0
w 0 /c 0.50 tan φ 0.6 ; 0.6
f1 ; f2 5.0 ; 8.0 tan φ residual 0.2 ; 0.2
α shr 0.01 GFI 0.03 ; 0.06
Diffusion analysis (interface elements) GFIIa 10x GFI
η (1/day) 100.0x10
6
σdil (MPa) 40
K 0 /K 1 0.01 p χ , p c , p GF 0.4, 0.5, 0.8
β k (−) 0.0 K χ , K c , K GF 1.0, 1.0, 1.0

Table 4.2. Adopted parameters for the diffusion and mechanical models that yield the
best fit to experimental results.

122
4.6.2. Simulation results
First, the results in terms of weight losses as a function of time are compared to the
experimental values in figure 4.42, for both coupled and uncoupled analyses, with the
same material parameters (from table 6.2). It can be observed that the effect of coupling
is small and may be neglected in this case, and that the numerical results agree well with
the experimental ones. However, when the same weight losses are contrasted with the
longitudinal strains (measured in the simulation as the average of left and right strains,
with the same experimental base of 50cm), as shown in figure 4.43, an important gap
exists between numerical (coupled and uncoupled cases) and experimental values.
Moreover, the effect of coupling yields smaller longitudinal strains, as compared to the
uncoupled case, for the same state of drying, although again the effect is small. This
result is to be expected, since the coupling effect manifests itself in a higher degree of
microcracking, thus decreasing longitudinal shrinkage strains.
In order to explain the above-mentioned discrepancies, the possible effect of
considering a formulation with the calculation of the moisture capacity matrix based on
the derivative of the desorption isotherm (see section 4.1.3.) has also been studied. The
same uncoupled simulation has been repeated with the introduction of this extra source
of nonlinearity, as shown in figure 4.44. The diffusion coefficient is in this case
considered as 0.25 times the values in table 4.2 (D0 and D1), in order to introduce a sort
of mean value of the moisture capacity between 0.5 and 1.0 RH (figure 4.44a). The
comparison in terms of longitudinal strains and weight losses is shown in figure 4.44b,
in which it can be observed that the difference is negligible, provided that the effective
diffusion coefficient is corrected. This result was not to be expected, since a highly
nonlinear desorption isotherm has been used with a moisture capacity varying between
1.1 and 0.10 [g/cm3], i.e. one order of magnitude.

Figure 4.42. Comparison of numerical (coupled and uncoupled analyses) and

experimental results in terms of weight loss (in %) evolution in time (days), and mesh
used throughout the simulations (the red line represents the spin axis).

123
Figure 4.43. Comparison of numerical and experimental results in terms of weight losses
(%) and longitudinal strains [mm/m] for coupled and uncoupled simulations with
constant shrinkage coefficient (of 0.01[cm3/g]), and RH distribution at 1,000 days.

Although the skin microcracking effect is well captured (large weight losses induce
small longitudinal strains at the beginning of drying), two important deficiencies arise
from the simulation: the numerical curve seems to be shifted in the weight loss axis with
respect to the experimental one, and the second curved part of the experimental curve, in
which strains grow more slowly than weight losses, is not well captured by the model.
The reasons why the numerical fit is not entirely satisfactory may be twofold.
On one side, the elasto-plastic nature of the model used for the zero-thickness
interface elements implies that the crack closure effect is not taken into account.
Accordingly, as with any elasto-plastic model, the deformation (in this case relative
displacements) is irreversible on unloading (elastic unloading is with the initial stiffness,
which is assigned a very high value in the present model, see Chapter 2). This fact may
be of importance, since the skin microcracks eventually unload and are even subjected to
compression stresses when the drying front advances towards the interior of the sample,
as shown in section 4.4. In this particular case, most of the skin microcracks unload at an
age of 7 days of drying, corresponding to a weight loss of 0.8%. If crack closure were
not prevented, the longitudinal strains would thus increase (crack closure may be
regarded as extra shrinkage strains), and the numerical curve would be shifted in the
weight loss axis (to the left). It is to be expected that, if the first non-linearity in the
strain vs. weight loss curve is due to the effect of microcracks opening, then the
consideration of the crack closure effect should be non-negligible in the global
relationship. Other authors have adopted continuum-based isotropic damage models in
order to account for cracking in drying shrinkage studies, thus allowing the crack closure
(Benboudjema, 2002; Benboudjema et al., 2005c; Bazant & Xi, 1994).
The second factor having an influence on the results is obviously the consideration of
a constant shrinkage coefficient (i.e. a linear relation between local shrinkage strains and
local weight losses, see section 4.1.4.). This could be the cause of the lack of a second
curved part in figure 4.43. Other authors have already proposed to adopt shrinkage
coefficients depending on the moisture content (van Zijl, 1999; Alvaredo, 1995).

124
 1.2

We/We(H=1)
(a)

dWe/dH
1 sorption isotherm
derivative
0.8

0.6

0.4
0.25
0.2

0
(b) 0 0.2 0.4 0.6 0.8
H 1
0.5
ε (mm/m)
0.4
C_const UNCOUPLED SIM.

0.3 C_var UNCOUPLED SIM.

0.2

0.1

w (%)
0
0 0.5 1 1.5 2 2.5 3

Figure 4.44. Influence of the moisture capacity matrix in the drying shrinkage
simulations. (a) Desorption isotherm used in the analysis and its derivative (0.25 is the
value adopted in the constant moisture capacity case). (b) Comparison of simulations
performed by assuming a constant and a RH-dependent moisture capacity matrix, in
terms of weight loss and longitudinal strains [mm/m] for an uncoupled analysis with
constant shrinkage coefficient (of 0.01[cm3/g]).

Because the introduction of the above mentioned effects would possibly result in a
better approximation to the experimental results, it was decided to repeat the simulations
with the consideration of a linear and a parabolic dependence of the shrinkage
coefficient on the RH (as shown in section 4.1.4.). Otherwise, all other material
parameters have been kept the same. For the shrinkage coefficient, a value has been
adopted such that the final local shrinkage is the one corresponding to a constant
shrinkage coefficient of 9x10-3[cm3/g]. The case of constant shrinkage coefficient has
been repeated with the new value. The Dirichlet boundary conditions have been replaced
by convective type boundary conditions with a film coefficient of 5[mm/day] (see
section 3.4.2.). The only reason for this change is that, in order to integrate incrementally
the nonlinear relationship between strains and weight losses in a more accurate way, a
smoother variation of the RH near the drying surface, at the beginning of the drying
process, is preferred. The results in terms of longitudinal strains vs. total weight losses
are presented in figure 4.45 for the uncoupled cases. A considerable improvement of the
global behavior can be observed, especially with respect to the second kink in the curve,
which is now satisfactorily captured. From the obtained final fit one may still suspect
that the fact of not considering the crack closure effect may also be of importance, and
that the introduction of this effect could perhaps improve the approximation, as
discussed above. Further investigation of this effect, however, would require a new

125
constitutive interface model with secant unloading, and it is therefore out of the scope of
the present thesis.

0.6
ε (mm/m)
0.5
PENLY CONCRETE
SIMULATION (ALPHA PARABOLIC)
0.4 SIMULATION (ALPHA LINEAR)
SIMULATION (ALPHA CONSTANT)
0.3

0.2

0.1

0
0 0.5 1 1.5 2 2.5 w(%) 3

Figure 4.45. Comparison of numerical and experimental results in terms of weight losses
(%) and longitudinal strains [mm/m]: influence of the shrinkage coefficient (constant,
linear and parabolic functions) on the uncoupled analysis of drying shrinkage.

4.7. Drying shrinkage under an external compression load

Even though a detailed systematic study of drying creep and the Pickett effect is not
part of the present work, some preliminary calculations have been made with the present
model to evaluate the behavior under simultaneous drying and mechanical loading. A
compression load has been applied to a numerical specimen in a direction perpendicular
to the drying direction. In this section, results obtained from simulations made on
15x15cm2 (27.6% aggregate volume fraction) concrete specimens subjected to drying
and different compression levels are presented. Initial and boundary conditions are
identical to the case in figure 4.19. Mechanical material properties are Ematrix(28days):
18,220MPa (Maxwell chain model), Eaggr: 68,000MPa (yielding an elastic modulus at 28
days of 25,000MPa for the overall concrete specimen), and rest of parameters equal to
table 4.2. For the diffusion analysis the parameters shown in table 4.3 have been
adopted.

Material parameters adopted

Diffusion analysis (matrix) Diffusion analysis (interface elements)
Initial humidity 100% η (1/day) 100.0x10
6

-5
D0 (cm2/day) 5x10 α hydration 0.70 K 0 /K 1 0.01
-2
D1 (cm2/day) 8x10 w 0 /c 0.40 β k (−) 0.0
β (−) 3.8 f1 ; f2 8.0 ; 6.0
3
C (g. cem/cm mat) 0.715 α shr 0.01

Table 4.3. Parameters adopted for the diffusion model in the simulations of simultaneous
drying and loading under compression.

126
 If during the drying process a certain level of compression load is applied
perpendicular to the direction of drying, it is to be expected that microcracking
perpendicular to the drying surface will be reduced, as shown schematically in figure
4.46. In fact, this reduction has been observed experimentally (see Sicard et al., 1992,
and section 3.1.5.1.) and by numerical simulations considering linear elastic behavior
(Sadouki & Wittmann, 2001). This reduction of cracking induces an increase of the
longitudinal effective shrinkage strain. The effective strain is defined here by subtracting
the time-dependent strain measured on the same test, but in which drying is prevented
(basic creep strain), from the total longitudinal strain. In the limiting case, when the
compression level is high enough so that microcracking is totally prevented, the
effective shrinkage strain would be maximized. This case can be simulated with the
present model by considering an elastic matrix (or viscoelastic) and without the addition
of interface elements (cracks are not allowed to form).

Figure 4.46. Schematic effect of simultaneous drying and mechanical compression load
on the microcracks perpendicular to the drying surface: (a) effect of drying only, (b)
combined effect at a low level of compression, and (c) at a higher compression level,
showing the reduction of microcracking.

As a first verification of the effect of simultaneous drying and loading, two identical
coupled simulations of free drying shrinkage have been performed, for which the only
difference is the constitutive law of the interface elements. In one case, a linear elastic
constitutive law for the interfaces is adopted (with KN and KT regarded as penalty
coefficients, with very high values, see Chapter 2) so as to simulate the case of no-
cracking, while the second simulation considers an aging elasto-plastic law (with
mechanical parameters adopted from table 4.2), for which cracking is not prevented at
all. Additionally, two levels of compression of 2MPa and 5MPa (peak load of around
30MPa) have been tested for both sealed (as in a standard basic creep test) and non-
sealed (i.e. drying) conditions. Such low levels of pre-compression have been chosen in
order to remain within the elastic regime of the material. In this way the subtraction of
basic creep strains from the total strain in the case of simultaneous drying and loading
would yield the deformations due to drying only.
Figure 4.47 presents the results obtained in terms of the evolution of longitudinal
strains in time for the different cases. It can be seen that the lower-bound limit is given
by the case of free drying shrinkage with microcrack formation. The upper-bound
represents the case of free-load drying shrinkage in which cracking is prevented (thus,
shrinkage strains are maximal). The intermediate curves, for loads of 2MPa and 5MPa,
are the ‘effective’ shrinkage strains defined above. It may be observed that the effective
shrinkage strain approximates the maximal shrinkage strain (for the case of no crack
formation) with the increase of the compression level from 0 to 5MPa. These results
agree with the assumptions made above and with observations made in the literature

127
(Sicard et al., 1992) and confirm previous numerical simulations made with linear
elasticity (Sadouki & Wittmann, 2001) also for the case of aging viscoelasticity.

Figure 4.47. Evolution of longitudinal strains [mm/m] with time (days) for the following
cases: free drying shrinkage with microcracking, drying shrinkage under a 2MPa
compression load, idem with 5MPa, and free drying shrinkage with prevented cracking.

As a second step in this study, a number of simulations were also performed with the
aim of reproducing the experimental series needed to study the Pickett effect (see
Chapter 3, section 3.2.2.), at three different compression load levels (5MPa, 10MPa and
20MPa). The goal is to assess the model response under simultaneous drying and
loading and to determine whether the Pickett effect can be well captured by the present
model or not. Using the same FE mesh as in the previous example (15x15cm2), drying
shrinkage, basic creep, and drying with simultaneous compression load simulations have
been performed. The results are presented in figure 4.48 for the three load levels. The
right trend is obtained in all cases, i.e. creep strain under drying conditions is larger than
the sum of basic creep and drying shrinkage strains.
However, the magnitude of the difference varies between 5 and 10%, which remain
still far from the magnitude of the observed Pickett effect. This conclusion agrees with
previous studies (e.g. Bazant & Xi, 1994) and supports the idea that the largest part of
the drying creep strains is due to an additional intrinsic physical mechanism and not to
the simple superposition of drying shrinkage and crack closure due to superimposed
loading, as was proposed in the past (see Wittmann & Roelfstra, 1980). Nevertheless, it
may be noted that the consideration in the model of the coupling between aging and
moisture diffusion (briefly discussed in Chapter 2), which has been left out of the
present work, could have an influence in the skin microcracking effect, yielding more
severe cracking due to a less hydrated material in the outer layers near the exposed
surfaces, as compared to the results obtained in this thesis. It can be concluded that the
present model may be certainly used to study drying shrinkage in concrete, but in order
to analyze drying creep tests an improvement of the formulation would need to be
introduced, as briefly discussed in Chapter 3 (section 3.5.2.). This improvement,
however, has been left out of the present work.

128
Figure 4.48. Longitudinal strains [mm/m] vs. time [days] curves for the following cases:
drying shrinkage, basic creep, drying plus basic creep, and drying under simultaneous
load, for a uniaxial compression load of 5MPa (a), 10MPa (b) and 20MPa (c).

4.8. Partial conclusions on HM modeling of drying shrinkage

A hygro-mechanical (HM) coupled model for the analysis of drying shrinkage of
concrete at the meso-level has been presented that is able to represent well the essential
features of this complex phenomenon. Non-uniform moisture distribution due to the

129
presence of aggregates and cracks, strain vs. weight loss relationship, stress profiles and
crack patterns are well captured by the proposed model. Cracking is introduced via zero-
thickness interface elements equipped with an elasto-plastic constitutive law based on
NLFM, and the additional moisture diffusion through open cracks is explicitly accounted
for.
The simulations presented in this chapter have shown that the effect of coupling is a
slight increase of the drying state but is small in most practical cases. In fact, the
difference (typically less than 15%) could as well be considered to fall in the range of
scatter of results. This is mainly due to the fact that the maximum crack openings found
in the simulations are in the order of 20µm, with typical values of 5µm (varying between
2 and 8µm) near the drying surface and 1µm for the inner microcracks. Thus, it is to be
expected that the effect of coupling is noticeable only at the beginning of drying, when
the drying through surface microcracks (of maximum crack openings) is most important.
This has been confirmed in the simulations, in which the iterations needed by the
staggered approach to converge were of the order of three times larger at the beginning
of drying than at more advanced stages. This is an important conclusion saying that
uncoupled analyses can be performed without major loss of consistency of the obtained
results and with significant savings in computational cost. Moreover, it should be noted
that the roughness (at the micro level) of the crack faces has not been considered in the
simulations (traditionally done by introducing a ‘hydraulic crack aperture’, with a
smaller value than the real crack width, see e.g. Segura, 2007). Thus, the results obtained
in this thesis should be considered as an upper bound of the effect of cracks on drying
(provided that the crack pattern is accurate enough), supporting the previous conclusion.
The effects of the aggregates on the drying-induced microcracking have been studied
in some detail. The performance of the model in this regard has proven to be satisfactory
in most cases, with the only observation that the effect of the aggregate size for a
constant volume fraction required filtering cracks with very small crack widths from the
data post-process. This could be due to the simplified mesostructure adopted in this
work: using mono-sized inclusions in a real 3D specimen does not mean that a 2D
section would contain mono-sized sections of these aggregates. In fact not even the
volume fraction is an exact value, since stereology only tell us that the average area
fraction of a (large) number of cross-sections will be equal to the real 3D volume
fraction, but not each single cross-section area fraction. Finally, the simulated 2D cross-
sections can not represent the microcracking due to out-of-plane aggregates, and thus a
lower degree of microcracking is expected. All these deficiencies can only be overcome
with a full 3D diffusion/mechanical analysis, as included in the final chapter on
perspectives for future research.
A series of simulations has also been performed in order to adjust model parameters
to experimental results on concrete specimens obtained by Granger (1996). Numerical
results agree well with experimental measurements. In addition, it has been shown that
the consideration of a nonlinear local relationship between shrinkage strains and weight
losses can be more accurate for simulating drying shrinkage experiments, which is an
important observation. Results have also underlined the possible need of introducing a
constitutive law for zero-thickness interface elements that accounts for the crack closure
effect. In this regard, a new damage-plastic model for interface elements would certainly
be an improvement and future (and on-going) work is aimed at this direction. On the
other hand, it has been found that considering the moisture capacity matrix as the
derivative of the desorption isotherm, although more consistent from a theoretical point
of view, has a negligible effect on the overall response.

130
 The analysis of drying under simultaneous loading has confirmed that, although the
correct trend is obtained, the well-known Pickett effect, or drying creep, cannot be
quantitatively represented by the skin microcracking effect. It is concluded that the
model presented in this thesis, as it is, cannot be used to study drying creep experiments,
and that a correct representation of this phenomenon would require the introduction of a
new intrinsic mechanism, as discussed in Chapter 3.

131
132
Chapter 5

EXTERNAL SULFATE ATTACK ON

SATURATED CONCRETE

The durability assessment and evaluation of the long-term performance of underground

nuclear waste containments, among other types of structures, are of vital importance due to
the huge consequences that their failure would have. Different forms of chemical attack,
sometimes with various types of ions present simultaneously, may be at the origin of the
possible degradation of concrete. In this chapter we focus on one of these degradation
processes, namely the external sulfate attack on concrete. Although it is a classical type of
deterioration widely described in concrete textbooks (see for instance Mehta & Monteiro,
1994, Neville, 2002, or Soroka, 1993), important issues regarding the degradation
mechanisms remain still not well understood. There is nowadays a renewed interest in
rationally describing the mechanisms behind expansive processes leading to cracking and
spalling of concrete structures exposed to sulfate solutions. Recent advances in the
experimental and also in the modeling fields have shown encouraging results towards this
goal. Thus, in this chapter a review of the most relevant progress in the sulfate attack field
will be addressed. In the first section of this chapter a summary of the main experimental
findings existing in the literature will be attempted. In the second part some proposed
empirical relations will be addressed, and the most relevant models developed to evaluate the
degradation due to external sulfate attack will be critically described and compared, trying to
highlight the main achievements and limitations. Finally, the need for a correct representation
of the cracking phenomena will be underlined.

5.1. Some experimental evidence of sulfate attack

Several forms of sulfate attack on concrete have widely been recognized (Hime & Matter,
1999; Skalny et al., 2002; Collepardi, 2003; Neville, 2004). There is a large amount of
experimental evidence on the different kinds of chemical or physical sulfate attack. It has
been argued that sulfate attack is not a decisive issue in concrete durability, due to the scarcity
of real cases where severe damage of concrete structures in the field have occurred (Neville,
2004; Skalny et al., 2002). In fact, it is very difficult to find in the literature real structures in
field conditions attacked by sulfate solutions (see Skalny et al., 2002 for some case studies).
However, there is a renewed interest in understanding chemical and physical origins of sulfate
attack and the main reason is the need to evaluate the long-term performance of underground
nuclear waste containments (Shah & Hookham, 1998; Gallé et al., 2006), where water
tightness is of major importance (Planel et al., 2006; Lothenbach & Wieland, 2006; Bary,
2008). Early attempts to model this long-term behavior were based on empirical relations of
limited applicability (Atkinson & Hearne, 1984), so that reliable predictions could not be
attained.

133
5.1.1. Fundamentals of external sulfate attack
In this thesis we restrict our attention to the study of external sulfate attack (ESA), which
may be defined as the mechanical deterioration of cementitious materials due to the ingress of
sulfate ions present in the surrounding sulfate-rich environment. By mechanical deterioration
it is meant that concrete may suffer cracking, scaling (erosion-type appearance of the concrete
surface), spalling, loss of overall strength, and even complete disintegration under severe
attack conditions (see for instance the theses by Schmidt, 2007 and Akpinar, 2006). Typical
expansion curves and spalling and degradation phenomena in mortar and concrete specimens
are shown in figures 5.1 to 5.3, under sodium sulfate as well as magnesium sulfate solutions.
These figures are illustrative of the extent of the degradation process in sulfate solutions with
high concentrations (although these conditions are rarely found in field cases). The main
sources of sulfates are the sodium, magnesium, calcium and potassium sulfates that may be
present in soils, ground water, solid or liquid industrial wastes, fertilizers or atmospheric SO3
(Skalny et al., 2002; Escadeillas & Hornain, 2008). This last source, often coming from acid
rain, is often studied separately as acid attack of concrete. It has been experimentally
determined that sulfates enter the specimen through a diffusion-driven process (Planel, 2002;
Le Bescop & Solet, 2006), as in the leaching process of cement pastes. The main conditions
that have to be fulfilled in order for ESA to occur are a sulfate-rich environment, a high
permeability (or diffusivity) of concrete and the presence of a moist environment, favoring
diffusion of sulfates (Collepardi, 2003). In this way, ESA may be described according to three
processes (Planel, 2002): 1) transport of the sulfate ions through the porous network, mainly
controlled by the permeability and diffusivity of concrete (the w/c ratio is the key parameter
in this regard), as well as through cracks present in the material; 2) chemical reactions
between the hardened cement paste (hcp) components and the sulfate ions (once these ions
have entered the material, the type of cement and the aluminates content will mainly
determine the degree of reactions occurring); 3) expansion phenomena as a consequence of
the formation of new crystalline phases.

Figure 5.1. (a) Crack pattern of concrete specimens attacked by sulfate solutions (sections
exposed on four edges), showing spalling and cracking phenomena (from Al-Amoudi, 2002).
(b) Expansions of mortar specimens with different w/c ratios as a function of the time of
exposure, showing the advantage of reducing the mixing water (from Naik et al., 2006).

134
Figure 5.2. 5mm cube mortar specimens (s/c = 2.0, w/c = 0.45, Type I cement, different
amounts of limestone filler replacement: (a): 0%, (b): 10%, (c): 20%, (d): 30%) exposed to
sodium sulfate NaSO4 solution (33,800ppm SO42-) during 1 year (from Lee et al., 2008).

Figure 5.3. 5mm cube mortar specimens (s/c = 2.0, w/c = 0.45, Type I cement, different
amounts of limestone filler replacement: (a): 0%, (b): 10%, (c): 20%, (d): 30%) exposed to
magnesium sulfate MgSO4 solution (33,800ppm SO42-) during 1 year (from Lee et al., 2008).

Mechanisms behind expansions

In general, the presence of sulfates from external sources results in the formation of new
phases like ettringite (3CaO·Al2O3·3CaSO4·32H2O, or C6A S 3H32 in cement chemistry
notation) and gypsum (C S H2). Most of the experimental evidence has shown that secondary
ettringite formation is the main factor behind expansions1 (Odler & Colán-Subauste, 1999;
Santhanam et al., 2001; Brown & Hooton, 2002; Al-Amoudi, 2002; Irassar et al., 2003). On
the other hand, gypsum’s effect on overall expansion is still an open question (Neville, 2004).
Some authors suggest that gypsum formation may lead to expansion (Wang, 1994; Tian &
Cohen, 2000a and references therein), based on specimens where ettringite formation was
eliminated (Tian & Cohen, 2000b; Santhanam et al., 2003a). In fact, it seems that gypsum is
the primary reaction product of sulfate attack at high sulfate-ion concentrations (> 8,000ppm
SO42-) (Santhanam et al., 2001). The formation of gypsum is related to the formation of other
products in sulfate attack, since it may be combined with other products to form ettringite.
The mechanism by which ettringite causes macroscopic or microscopic expansion is still
not well documented and there is no general consensus. Several proposals can be found in the
literature, of which none is able to directly correlate ettringite concentration to observed
expansions (Brown & Taylor, 1998; Skalny et al., 2002; Planel, 2002). In fact, this is
confirmed by several experimental observations, in which this correlation is usually also
missing. The mere exercise of reviewing the different nature of proposed mechanisms is a
good indicator of the level of controversy that exists around ettringite-related expansion. Such
proposals include (non extensively) that a) expansions are the result of an increase of the solid

1
Primary ettringite forms in fresh concrete with no harmful consequences and then disappears during hydration
to form sulfate and monosulfate hydrates, see Brown & Taylor, 1998.

135
volume, b) expansions are due to topochemical reactions (the reaction product is formed in
the space originally occupied by one of the reactants) or c) through-solution reactions (the
reaction product precipitates randomly from the liquid phase), d) expansions are due to crystal
growth pressure, e) or that expansions are caused by swelling phenomena. A complete review
on several mechanisms can be found elsewhere (Brown & Taylor, 1998; Planel, 2002; Skalny
et al., 2002). Moreover, the possibility of simultaneous action of more than one of these
mechanisms may be the cause of the overall expansions (Planel, 2002). In fact, it has been
detected that the presence of ettringite is not necessarily followed by expansions. The reason
is that the presence of pores of different sizes and cracks present in the cement paste and the
paste-aggregates interfaces provide space for ettringite precipitation. It has been suggested, on
the basis of numerical simulations, that the pore size distribution plays a vital role in the
resulting overall expansions (Schmidt-Döhl & Rostásy, 1999). Brown & Taylor (1998) stated
that any theory of ettringite expansion must account for the facts that ettringite forms by a
through-solution mechanism, exhibits a true solubility and must occupy more space than that
which was initially available to it.
In the present study, following the work by Krajcinovic et al. (1992), Clifton &
Pommersheim (1994) and Tixier & Mobasher (2003a,b), the idea that the solid reaction
product volume is larger than the solid reactants volume has been adopted to explain and
quantify the local expansion. Surprisingly, the review by Skalny and coworkers (2002) points
out that most of the reactions in which sulfate ions are involved are associated with a
chemical shrinkage (i.e. the volume of reaction product is smaller than the sum of the volume
of the reactants). Figure 5.4 shows the relation between the expansions observed in different
cement paste specimens under external sulfate attack (with water curing and humid air curing)
and the measured amount of formed ettringite (Odler & Colán-Subauste, 1999). It may me
clearly observed that there exist a threshold amount of ettringite that has to be formed before
significant expansions are measured, and that high levels of expansion were accompanied
with high ettringite concentrations. The comparison of the mass change (sulfate-induced mass
increase) and the overall expansion of mortar samples (Schmidt, 2007), as presented in figure
5.5(a), shows a moderate expansion at the beginning of attack, while mass change rapidly
increases. These results also support the idea that during the first period the sulfates uptake
mainly results in the filling of the pores and air voids initially present. Once this space is filled
a further increase in the mass change is followed by substantial global expansions. The pore
filling process with reaction products is well documented from microstructural experimental
observations.

Figure 5.4. Expansions of different cement paste specimens as a function of the amount of
formed ettringite for water curing (a), and humid air curing (b) conditions (from Odler &
Colán-Subauste, 1999).

136
 Figures 5.5(b) (Schmidt, 2007) and 5.6 show the precipitation of ettringite in air voids and
pores (Brown & Hooton, 2002). It may be observed that the crystals do not necessarily fill
them entirely (Brown & Hooton, 2002) but there exists a ‘macro porosity’ of the reaction
product. These observations are in favor of the idea that expansion starts well before the pore
space is filled by reaction products, which is usually used as a modeling hypothesis, as it will
be shown in the next section.

Figure 5.5. (a) Mass change versus expansion of different mortar specimens immersed in
sodium sulfate solution at 8ºC and 20ºC. In these cases, ettringite as well as thaumasite were
formed. (b) Ettringite formation in pores, voids and cracks of the same mortar after 180 days
of exposure at a depth of 0.5mm (from Schmidt, 2007).

Figure 5.6. (a) Ettringite-filled air void and typical tiger stripe morphology of ettringite.
Cracks radiating from the ettringite propagate through the paste, denoting expansion and
distress induced by reaction product formation. In these cases, ettringite as well as thaumasite
and gypsum were formed. (b) Ettringite partially filling air voids in the paste at 20mm below
the exposed surface, and reduced resulting crack density (from Brown & Hooton, 2002).

Main chemical reactions involved

The most relevant chemical reactions that potentially take place during the exposure to a
sodium sulfate solution are well-known and are summarized below (Cohen, 1983; Neville,
2002; Planel, 2002; Tixier, 2000; Skalny et al., 2002; Schmidt, 2007; Escadeillas & Hornain,
2008).

137
- Formation of gypsum from portlandite (CH2, eq. 5.1a) and CSH (eq. 5.1b):
CH + Na 2SO4 + H 2 O → CSH 2 + 2NaOH
(5.1)
CSH + SO24 − + H 2 O → CSH 2

- Formation of ettringite from monosulfoaluminate (eq. 5.2a), unreacted tricalcium aluminate

grains (C3A) (eq. 5.2b), tetracalcium aluminate hydrate (C4AH13) (eq. 5.2c), and alumino-
ferrite phase (C4AF, eq. 5.2d), although this last reaction seems to be of little importance due
to its low reaction kinetics (see Schmidt, 2007):
C4 ASH12 + 2CSH 2 + 16H → C6 AS3 H32
C3 A + 3CSH 2 + 26H → C6 AS3H 32
(5.2)
C4 AH13 + 2CH + 3S + 17H → C6 AS3H32
3C4 AF + 12CSH 2 + xH → 4 ( C6 AS3H 32 ) + 2 ( A, F ) H 3 

The formation of gypsum is mostly important for high sulfate concentrations, as in the case
of laboratory samples (Santhanam, 2001; Santhanam et al., 2001; Monteiro & Kurtis, 2003),
but its presence in field conditions, under low sulfate concentrations, is rather questionable
(Marchand et al., 2002). Precipitation of gypsum in the cracks has been detected by many
authors as gypsum veins (see Glasser et al., 2008; Planel, 2002 and references therein),
sometimes filling them completely, although its presence is thought not to have caused these
cracks. Gypsum formation is usually considered to be the cause of loss of cohesion or strength
of the cement paste as a result of decalcification of the CSH phase (Planel, 2002; Shazali et
al., 2006). Its participation in the observed overall expansion is, as stated above, still not clear
(Tian & Cohen, 2000b) and will not be considered in the present study.
5.1.2. Factors affecting sulfate attack
The intensity of the attack (or degradation extent) is very difficult to predict in real
structures and depends on many factors, although they may be divided in two groups:
- concrete quality: within this category, we may cite the type of cement, the w/c ratio,
placing and curing conditions, concrete deterioration before the attack, mineral admixtures
types and content, etc.;
- environmental conditions: concentration, distribution and type (sodium, magnesium, etc.)
of sulfates, humidity and temperature, groundwater flow (defining the pH of the solution),
permeability of the soil in contact with the structure, combined effect of different degradation
processes, etc.
a) Effect of the type of sulfate
It is now well-known that the type of sulfate (sulfate anion plus a cation) has a great
influence on the degradation process, being the most aggressive the magnesium sulfate
(MgSO4), followed by sodium (Na2SO4) or potassium (K2SO4) sulfates (alkali sulfates) and
calcium sulfate (CaSO4) (Neville, 2004; Escadeillas & Hornain, 2008). The last two types of
sulfates (i.e. alkali and calcium sulfates) act in a similar way (Schmidt, 2007). Briefly, alkali
sulfates react with calcium ions provided by the calcium hydroxide (CH or portlandite) or
eventually from decomposition of calcium silicate hydrates (CSH) to form gypsum
(CaSO4·2H2O, or C S H2 in cement chemistry notation). This process involves calcium

2
From now on the standard cement chemistry notation will be used: CaO2 = C, Al2O3 = A,
SO3 = S , SiO2 = S, CO2 = C , Na2O = N and H2O = H

138
leaching of the cement paste to a degree that depends on the concentration of sulfate solution,
among other factors. Consecutively, the formed gypsum reacts with the aluminate phase,
which is mostly present as monosulfoaluminates (C4A S H12, often referred to as an AFm
phase), but also as unreacted tricalcium aluminate (C3A) or tetracalcium aluminate (C4AH13),
to form ettringite (C6A S 3H32, often referred to as an AFt phase) (Brown & Taylor, 1998;
Skalny et al., 2002). In the case of high concentration of sulfates in the solution, when the
aluminates are depleted, gypsum starts to be formed rather than ettringite (Skalny et al., 2002;
Planel, 2002), although the formation of gypsum and its effects are a matter of controversy
(Neville, 2004). The difference between an alkali sulfate and a calcium sulfate attack is that in
this last case a calcium source is not needed, so that portlandite depletion is reduced, leading
to a lower amount of leached calcium from the cement paste (Neville, 2004; Schmidt, 2007).
In contrast, under magnesium sulfate attack, depletion of portlandite is accelerated, forming
magnesium hydroxide (brucite) and gypsum. Gradual decomposition of the CSH phase takes
place (significantly faster than with other sulfates), converting it into an amorphous hydrous
silica (Brown & Taylor, 1998) without mechanical properties. The amount of ettringite
formed is in this case rather low because the hcp tends to disintegrate, due to degradation of
CSH phase, before significant amounts of ettringite may be formed (Skalny et al., 2002).
Thus, magnesium sulfate attack is characterized by loss of strength and disintegration of
concrete, rather than scaling, spalling and expansion, which are the characteristic symptoms
in the case of sodium sulfate attack (Skalny et al., 2002). In the present study attention will be
exclusively focused on sodium sulfate attack (Na2SO4), mainly because it is the best
documented type of attack (and perhaps better understood), at least experimentally, and
because it has been suggested that experiments with magnesium sulfate solutions are not valid
to characterize field conditions (Neville, 2004). In the field, however, the general scenario is
given by a mixture of ions, and several types of sulfates may be present simultaneously (as
well as other type of ions, like chlorides), which is out of the scope of this thesis.
Recently, an interesting mechanism for sodium sulfate attack of mortar specimens has been
proposed, based on a vast experimental campaign, that considers chemical (microstructural)
as well as mechanical (physical) aspects of the deterioration process (Santhanam et al., 2002,
2003b) and may be summarized as follows:
1) Ettringite and gypsum formation near the exposed surface of the mortar result in an
attempt of expansion of the specimen’s skin, generating compressive stresses in this
layer and tensile stresses (self-equilibrated) in the inner unaltered layers;
2) These tensile stresses result in cracks appearing in the interior of the mortar;
3) The surface zone continues to deteriorate and the solution reaches the cracked interior
zone, reacting with hydration products and leading to their deposition within the cracks
and cement paste;
4) This new region tries to expand, causing cracking in the interior zone, so that three
distinct zones can be identified: disintegrated surface, reaction product’s deposition
zone, and chemically unaltered interior cracked zone;
5) The attack progresses at a steady rate until complete disintegration of the specimen.
This mechanism is based on their experimental results (Santhanam, 2001) and is in
agreement with other experiments found in the literature (Planel, 2002) as well as with what is
expected from purely mechanical considerations. Contrary to the case of drying shrinkage,
where tensile stresses in the outer layers and compressive stresses in the inner ones act in a
self-equilibrated manner at the beginning of drying, under an external sulfate attack scenario

139
these stresses change their sign, being compressive in the outer layers and tensile in the inner
ones.
b) Effect of the pH of the solution
Another important factor influencing the extent of external sulfate attack is the pH value of
the attacking solution and whether it is controlled (i.e. at a fixed value) or not. In fact, the
ASTM standard test method (ASTM C1012, 1995) does not consider a controlled pH in the
solution, which will vary in time due to the relative low volume ratio between attacking
solution and specimen. Instead, a controlled pH is rather what is found in field conditions due
to underground water flow and a large volume of the solution. The ettringite is stable at a pH
of 11.5 and if it changes gypsum rather than ettringite is formed, so that the influence in the
degradation process is notorious. This was first documented in the early 80’s (Brown, 1981)
and is now well recognized that accelerated attack conditions can be attained under controlled
pH (the pH value does not seem to affect expansions as much as the fact of not controlling it),
as shown in figure 5.7(a).
c) Effect of the initial C3A content
It is well-known that the initial C3A content of cement plays a very important role in
expansions due to sulfate attack, especially in the case of exposure to a sodium sulfate
solution (Ouyang et al., 1988; Akpinar, 2006). In this last case, sodium sulfates react with
portlandite to form gypsum and the latter reacts mainly with monosulfoaluminates to form
ettringite. A low C3A content minimizes the monosulfoaluminate content and thus the
formation of ettringite, leading to much more reduced expansions (Ouyang et al., 1988;
Skalny et al., 2002), as shown in figure 5.7(b). This is why Type V cements, with a very low
C3A content, were developed. It should be noted that expansions have also been observed in
pastes where the initial C3A content was almost eliminated (González & Irassar, 1997). In
these experiments the C3S content of the unhydrated cement (before mixing) was the key
factor determining expansions.

Figure 5.7. Evolution of overall expansions of mortar specimens. (a) effect of the pH of the
0.35M sodium sulfate solution, showing higher expansions for controlled pH solutions (from
Brown, 1981); (b) effect of reducing the C3A content on expansions under external sulfate
attack (from Ouyang et al., 1988).

d) Effect of combined leaching and sulfate attack

Recently, some authors have studied the effect of combined leaching and sulfate attack of
cement pastes or mortars (Planel et al., 2006; Marchand et al., 1998; Skalny et al., 2002;
Samson & Marchand, 2007; Bary, 2008). In fact, leaching and sulfate attack cannot be

140
generally separated (Schmidt, 2007; Planel, 2002). Calcium leaching, or simply leaching, is
the process by which cement paste progressively decalcifies (dissolution of calcium
hydroxide, or portlandite in cement chemistry, which could also result in decalcification of
calcium silicate hydrates) when subjected to pure weakly deionized water. It is also a
diffusion-driven phenomenon between the highly charged interstitial water in the pore
solution of the specimen and the aggressive solution (weakly charged), producing the
transport of ions by diffusion and causing dissolution of some of the cement paste constituents
(see e.g. Planel, 2002; Gérard et al., 2002; Le Bescop & Solet, 2006). An interesting state-of-
the-art of the leaching process and the main factors affecting calcium leaching may be found
in (Nguyen, 2005). It has been shown that this process negatively influences the mechanical
(local decrease in strength) and diffusive properties (increase in diffusivity) of cement-based
materials (Marchand et al., 1998; Haga et al., 2005). Recent studies have shown that under
combined leaching and external sulfate attack, the leaching process is not greatly affected and
follows approximately the same kinetics than that in deionized water (Planel et al., 2006),
thus depleting portlandite and increasing the diffusivity near the surface. However, it has been
suggested that the quantity of calcium ions released under combined attack is significantly
reduced as compared to exposure to pure water (Le Bescop & Solet, 2006). It is emphasized
that in this work a combined analysis of leaching and external sulfate attack has not been
performed. Instead, following the work by Mobasher and coworkers (Tixier & Mobasher,
2003a,b; Mobasher & Ferraris, 2004), it is considered that the penetrating sulfate ions
instantaneously react with the available portlandite (CH) and calcium silicate hydrates (CSH)
to form gypsum (CaSO4·2H2O). In this sense, it could be argued that the depletion of
portlandite is considered. However, the implicit assumption is made that calcium supply
remains constant. Finally, neither the decrease of mechanical properties (which is expected
not to be very important) nor the increase in porosity due to leaching (competing with the
decrease due to pore filling effect) have been considered in the present work.
5.1.3. Other kinds of sulfate attack
Other kinds of sulfate attack, which have been left out of the present study, are
summarized below. It is important to differentiate between these kinds of attack in order to
completely delimit the problem studied. The model which will be presented in the next
chapters is, at present, able to simulate the external sulfate attack as defined above and for
saturated and isothermal conditions, although its adaption to some other kinds of sulfate
attack presented herein would be straightforward provided the chemistry is well understood,
especially in the case of internal sulfate attack.
a) The internal sulfate attack (ISA), often named delayed ettringite formation (DEF), has
been defined as “the formation of ettringite in a cementitious material by a process that begins
after hardening is substantially complete and in which none of the sulfate comes from outside
the cement paste” (Taylor et al., 2001). It is also referred to as heat-induced internal sulfate
attack, because of its dependence on high curing temperatures (typically in steam-curing
conditions, but also in mass concrete). The necessary (but not sufficient) conditions for
expansion to occur from DEF are that the internal temperature must be above 70°C
(approximately) for a sufficient period of time, and that after temperature returns to normal
values the moisture content must be sufficiently high (Glasser, 1996; Taylor et al., 2001;
Collepardi, 2003). In contrast with external sulfate attack, expansions over time usually
present an S-shape type, so that expansions are bounded. ISA has been first detected around
twenty years ago in concrete railway ties, so that literature is rather limited and DEF origins
are not yet well understood. Indeed, some authors argue that ISA may occur even at normal
temperatures and propose to intimately relate it to pre-existing microcracks in concrete
(Collepardi, 2003; Diamond, 1996). The main hypotheses for explaining the mechanism of

141
expansion are the crystal growth at aggregate interfaces (large crystals) and crystal growth of
much smaller size formed within the cement paste, between CSH layers (Taylor et al., 2001).
Finally, it should be noted that although in recent literature ISA and DEF have become
synonyms, ISA also includes the case of attack by excess (with respect to the clinker
aluminate phase) of cement sulfate (Skalny et al., 2002).
b) The thaumasite form of sulfate attack (TSA) may occur in principle either in an ESA or
ISA scenario (although most of the reported cases involve ESA). TSA deterioration may be
much more severe than conventional sulfate attack as it is associated with depletion of the
CSH, which represent the main binding phase of the hcp (Macphee & Diamond, 2003).
Thaumasite (CaSiO3·CaCO3·CaSO4·15H2O, or C3S C S H15 in cement chemistry notation) is
formed by a general reaction involving calcium ions, silicate, sulfate, carbonate and sufficient
water to permit both transport of the potentially reactive species and to form thaumasite where
these species can interact under the low temperature conditions prevailing (Bensted, 1999),
being carbonate the key ingredient in this case. It is also a temperature dependent reaction,
being its formation favored by low temperatures (< 10°C), although it has been found to be
stable also at ambient temperature (Brown & Hooton, 2002; Schmidt, 2007). The carbonate
ions may be supplied by carbonate rocks constituting the aggregates, limestone addition,
atmospheric CO2, or ground water with high concentrations of dissolved CO2 (Skalny et al.,
2002; Schmidt, 2007). One important difference with conventional ettringite form of sulfate is
that expansions in TSA may be unbounded, since thaumasite can form as long as external
sulfates and carbonates are available. Indeed, the needed silica is supplied by the CSH, so that
all this phase could in theory be destroyed in an indefinitely sustained attack (Macphee &
Diamond, 2003). Contrary to this, formation of ettringite in conventional sulfate attack is
associated and limited with the alumina content of the cement. Once the alumina phase is
depleted, expansion stops or, at least, it is strongly reduced. Finally, it should be emphasized
that TSA may lead to significant spalling, loss of strength, and eventually converts concrete
into a structureless mass (Schmidt, 2007).
c) It has been proposed in the literature to consider the physical sulfate attack (often called
sulfate salt crystallization) as a special form of attack, in the sense that it should not be treated
as a chemical sulfate attack. Accordingly, a difference should be made between physical and
chemical attack, the former involving crystallization of salts (of which one example is sulfate)
exerting pressure over the pore volume (thus causing expansion), and the latter necessarily
involving the sulfate ion (Neville, 2004). Briefly, physical sulfate attack occurs when
crystallization takes place in the pores and the pressure built-up provokes damage of the
material (Al-Amoudi, 2002). Other authors claim that this division into chemical and physical
sulfate attack is not so clear, as reality shows “chemical processes with physical
consequences” (Skalny et al., 2002). In any case, this discussion is out of the scope of this
study, where focus will be exclusively made on external sulfate attack (involving the sulfate
ions).
5.1.4. Final remarks on experimental evidence of sulfate attack
In the preceding paragraphs a summary of the most relevant features of sulfate attack in
cementitious materials has been presented, as well as a brief discussion of the main factors
affecting the degradation of concrete in a sulfated environment. It is clear from this review
that a complete treatment of the problem should involve chemical as well as mechanical
aspects of sulfate ingress and their consequences in the overall behavior, in order to reliably
predict the durability of concrete structures under sulfate attack. There is a high level of
complexity due to the large number of the intervening factors and the extent to which they
affect the macroscopic response. The most important factors are on one side the sulfate

142
concentration, pH and the type of cation in the solution (sodium, magnesium, etc.), and on the
other hand the w/c ratio of the material as well as its C3A content, etc. The precise way in
which chemical reactions produce changes in the mechanical properties of the specimens
remains an open question. Ettringite is believed to be the key factor behind expansions,
although a widely accepted mechanism to explain expansion and degradation is still missing.
Instead, a number of proposed mechanisms coexist. The role of gypsum formation is also
controversial. Whether some authors believe that gypsum is the responsible for loss of
strength and cohesion, others affirm that its presence or formation does not cause significant
changes in mechanical properties. It can be concluded that there is still room for new
experimental campaigns aiming at precisely explaining, among other issues, the exact origin
of expansions, and that advanced models based on a complete treatment of the chemo-
transport-mechanical problem should provide new insight in order to design effective
experiments.

5.2. Modeling of external sulfate attack

In the previous section it has been shown that external sulfate attack mechanisms are not
yet clearly understood. The mechanism by which the presence (or formation) of ettringite
causes overall expansion, and the participation of gypsum in the degradation process are
perhaps the most important issues in sulfate attack, and yet the most controversial ones. As a
consequence, discrepancies and lack of consensus on the origin of resulting expansions and
cracking phenomena have obviously been translated to the computational modeling field. In
fact, as will be shown in the following paragraphs, only a few models dealing with external
sulfate attack can be found in the literature, some of them very recent, and even fewer able to
perform a chemo-transport-mechanical analysis, either under coupled or uncoupled
conditions. Instead, given the high level of complexity of the problem, unilateral efforts have
been made in different fields, such as: 1) performing chemical equilibrium analyses of initial
compounds and reaction products, basically dealing with chemical speciation and
composition of binder systems at equilibrium (without any contribution from mechanics or
transport processes); 2) performing single or multi-ionic transport analyses of chemical
species, where different chemical species and compounds profiles, such as sulfate or ettringite
concentrations, are sought (also disregarding mechanical consequences and their interaction,
although in this case some of the models proposed consider the reactive transport problem,
thus performing chemical equilibrium calculations); 3) elaborating empirical and
phenomenological models to estimate e.g. the spalling depth or the evolution of expansions of
specimens or structures subjected to sulfated environments, lacking a sound description of the
chemical processes or rigorous mechanical foundations; 4) and finally some advanced
mechanical models for cementitious materials have been proposed recently that incorporate a
simplified chemical approach to sulfate attack and the transport of ions considered in the
chemical analysis. It may be noted that the model proposed in this thesis lies in this last
category. However, the model described in this work is the only one capable of performing
simulations at the meso-scale, i.e. explicitly considering the main heterogeneities of the
material and the effect of cracking in the degradation process.
5.2.1. Chemo-transport models for sulfate ions
It should be emphasized the distinction made in the literature between purely chemical
analysis and transport processes. Traditionally, these fields have been considered as separated
scientific disciplines, with little contact between them. However, recent advances, mostly in
the field of geochemistry in hydrogeological systems, have closed this gap and broaden the
applicability of existent geochemical models with the development of so-called reactive

143
transport models (for a comprehensive review on this topic, see van der Lee & De Windt,
2001, van der Lee et al., 2003; Samson et al., 2000 or Samson & Marchand, 2007). The
resolution of this type of analysis is often called a mixed problem, since an algebraic system
of equations must be solved for the chemical analysis simultaneously with a partial
differential system for the transport problem (Samson et al., 2000; Samson & Marchand,
2007). It is only in the late 90’s that reactive transport models have started to be adapted to
study cementitious materials like concrete, with increasing success (Damidot & Le Bescop,
2008).
5.2.1.1. Modeling of the composition at chemical equilibrium
Within the field of chemical analysis, there are mainly two approaches that have been
mostly followed in the literature to numerically calculate the composition of a mixture at
equilibrium: by satisfying the equations of mass balance through the well known law of mass
action (or LMA) approach (Planel, 2002; van der Lee et al., 2003; Samson & Marchand,
2007), or by minimizing the Gibbs free energy of the chemical system (GEM approach),
techniques which are often referred to as thermodynamic modeling (Lothenbach &
Winnefeld, 2006; Schmidt et al., 2008). It should be noted that these two approaches are
thermodynamically equivalent (Guimarães et al., 2007) and yield almost the same results.
Nowadays, with the advent of readily accessible computational modeling, thermodynamic
modeling is becoming a widely used approach to calculate a composition at chemical
equilibrium (Gordon & McBride, 1994; Guimarães et al., 2007), providing equilibrium exists,
and has emerged as a powerful tool to analyze cementitious materials (Ayora et al., 1998;
Sugiyama & Fujita, 2006; Lothenbach & Winnefeld, 2006; Matschei, 2007). Recently, it has
been proposed to analyze sulfate attack on cement paste through the use of thermodynamic
modeling (Schmidt, 2007; Schmidt et al., 2008). This technique is based on the minimization
of the Gibbs free energy of a system to find its composition at chemical equilibrium supported
by consistent thermodynamic databases, when available, and assuming thermodynamic
equilibrium between different phases at a given time (see e.g. Lothenbach & Winnefeld, 2006
and references therein). The method has become a powerful and robust tool to analyze
chemical speciation due in part to the improvement and updating of the thermodynamic
databases, which are a key ingredient in this type of calculations. Several computer codes, of
which some are freely available, have been developed over the last years (see van der Lee &
De Windt, 2001 for a list of available programs). However, it is important to note that this
type of simulations do not provide any information on the kinetics of reaction, which is an
essential feature to study the evolution of a system (Lothenbach & Winnefeld, 2006; Schmidt,
2007), nor on the transport processes of the ions involved. Indeed, the minimization is
performed at fixed time and space, so that if the evolution of the system is to be determined,
additional input data concerning kinetics of reactions have to be provided, and equilibrated
compositions have to be calculated at each time increment. This drawback has been often
remedied by the addition of empirical relations giving information on kinetics (Lothenbach &
Winnefeld, 2006), with the consequent loss of consistency, since kinetic databases are much
smaller compared to thermodynamic ones and are more dependent on chemical conditions
(van der Lee & De Windt, 2001), thus adopting a more empirical character. It should be
emphasized that until a general consensus is reached on what is happening at the chemical
level, it will be very difficult to obtain fully reliable predictive models of the concrete
mechanical response under sulfate attack.
5.2.1.2. Modeling of the transport processes including chemical reactions
The transport processes modeling field is much more familiar to civil engineers. A
diffusion-driven phenomenon such as drying of concrete is a good example and has been

144
already presented in detail in the preceding chapters. The first proposals of ions migration
through cementitious materials dealt with the transport of a single ion. The selection of the ion
considered obviously depends on the problem studied. In this way, for the case of sulfate
attack in concrete, the sulfate ions have been mostly used (Tumidajski et al., 1995;
Gospodinov et al., 1996; Gospodinov et al., 1999; Gospodinov, 2005; Tixier, 2000). Other
examples are the calcium ions in the case of leaching (Gérard et al., 2002), chlorides in the
classical case of corrosion (Andrade, 1993; Oh & Jang, 2007; Glasser et al., 2008), or carbon
dioxide in the case of carbonation (Saetta & Vitaliani, 2005; Thiery, 2005; Ferreti & Bazant,
2006). In this approach, chemical reactions are included as a sink term (Gospodinov et al.,
1996) in the mass balance equation. This term explicitly contains the kinetics of reaction,
which typically takes the form of a zero, first or second order reaction. The chemical reaction
itself is implicitly taken into account by directly defining the amount of the reaction product
that is formed at the expense of the diffusing ion depletion. More sophisticated models
consider the transport of several ions simultaneously, as in the multi-ionic diffusion analyses
(Samson et al., 1999; Truc et al., 2000; Revil & Jougnot, 2008; Krabbenhoft & Krabbenhoft,
2008), restricting themselves to the consideration of the most important ions and their
electrical interaction, but disregarding chemical reactions. These models consider one mass
balance equation for each ion plus an additional equation that couples all of the previous ones,
regarding the electrical field between ions (either by the electroneutrality condition or the
Poisson’s equation). The complexity of these models is increased (and so is uncertainty) when
unsaturated conditions are introduced (Revil & Jougnot, 2008). In this case, some additional
features have to be taken into account, such as the transport of ions by convection and the
transport properties in gaseous and liquid phases. Finally, the last generation of transport
models in cementitious materials includes a sound description of the chemical reactions,
through the so called reactive transport approach (Marchand et al., 2002; Samson et al.,
2000; Samson & Marchand, 2007; Glasser et al., 2008; Guillon, 2004; Galíndez et al., 2006;
Damidot & Le Bescop, 2008). In these models, chemical reactions are studied simultaneously
(either in a fully coupled way or in a staggered way) with the transport equations, by
considering either the LMA (Marchand et al., 2002; Planel, 2002) or the GEM (Guimarães et
al., 2007) approach. Recent efforts have been made to simulate the chemical reactions and
transport of ions of specimens subjected to external sulfate attack, with encouraging results
(Marchand et al., 2002; Planel, 2002; Glasser et al., 2008). This type of simulations should
definitely serve to improve our understanding of the underlying phenomena behind
macroscopic degradation due to sulfate ions ingress. However, in order to achieve this, the
mechanical properties of the individual components (e.g. CSH, CH or ettringite) should be
known and more information on the internal pore structure of the hcp would be needed. At
present, this data is still rather scarce, although recent developments have permitted to obtain,
e.g., ettringite, CH and CSH elastic properties (see e.g. Speziale et al., 2008 or Bary, 2008
and references therein). In this thesis, focus is made on the mechanical behavior of concrete at
the meso-level, trying to gradually incorporate the fundamentals of chemical degradation,
reason by which these analyses at the micro-scale are out of the scope of this study.
5.2.2. Models for the degradation of cementitious materials under sulfate attack
5.2.2.1. Empirical and phenomenological models
From a purely civil engineering point of view, efforts have been focused on developing
empirical relations that may serve to quantify the degradation of concrete structures when
exposed to sulfate attack, based on experience. Their main drawback is that neither of the
proposed models is capable of accounting for more than two or three parameters (typically the
C3A content of the cement, the sulfate concentration of the solution and/or the w/c ratio),
whereas for the particular case of sulfate attack it is well known that there is a large number of

145
factors affecting the overall response of concrete exposed to sulfate solutions. Thus, their use
for extrapolating to cases that have not been taken into account in the formulation is limited.
The most complete empirical model, based on an extended experimental campaign for over
forty years on more than 100 concrete specimens, is probably that one proposed by Monteiro
and coworkers (Kurtis et al., 2000). In this case, it is the expansion evolution of concrete
specimens that is calculated in terms of the C3A content and the w/c ratio of concrete, instead
of the degradation depth or spalling depth of previous models (see for instance Atkinson and
Hearne, 1984, later improved in Atkinson & Hearne, 1990). The expressions proposed for two
different ranges of C3A content are as follows,
EXP(%) = α1 + α 2 ⋅ t(years) ⋅ w / c + α 3 ⋅ t(years) ⋅ C3 A(%), C3 A < 8%
(5.3)
ln ( EXP(%) ) = α1 + α 2 ⋅ t(years) + α 3 ⋅ ln ( t(years) ⋅ C3 A(%) ) , C3 A > 10%

where α i are model parameters, w/c is the water cement ratio, and C3 A(%) is the
concentration of calcium aluminates. The difficulty of translating expansions of concrete
specimens to a prediction of the degradation process of concrete structures is an important
disadvantage of this empirical model, although it does permit an easy evaluation of the
influence of two key parameters as are the w/c ratio and the C3A content.
Another empirical relation, proposed more than a decade earlier than that of Kurtis, was
developed by Atkinson & Hearne (1990). This model has some remarkable features, as it is
based on a more mechanistic point of view of the sulfate attack problem and on the fact that it
can be treated as a diffusion-driven phenomenon. They considered that ettringite is the only
cause of expansion (assuming that only 5% of the volume of ettringite formed is translated
into actual expansions), and that due to these imposed volumetric strains, the solid is
subjected to a stress field able to produce a crack when the elastic energy reaches a
predetermined crack surface energy. With these hypotheses they established the rate of
degradation (R) as (from which the degradation or spalling depth X can be extracted)
X E ⋅ β 2 ⋅ c 0 ⋅ CE ⋅ Di
R= = (5.4)
t α ⋅ γ ⋅ (1 − υ )
where t is the time of exposure α ⋅ γ accounts for the fracture energy of concrete, E and ν are
respectively the Young modulus and Poisson’s ratio, c0 is the external concentration of sulfate
solution (mol/m3), CE is the formed ettringite (mol/m3), β the coefficient of proportionality
between ettringite formed and expansions (β = 1.8x10-6) and Di is the diffusion coefficient of
concrete. The most important parameters in this model are the diffusion coefficient and the
fracture energy. The predictions made for two types of concrete were rather crude, with a big
gap between model and experimental results, but some of the ideas introduced in this model
have been of great value for subsequent proposals.
Additionally, some phenomenological models have been proposed in the literature to
predict the expansion, based on a somewhat more profound understanding of the sulfate
attack mechanism. The model by Clifton & Pommersheim (1994) is within the most cited
work within this category, and it will be briefly presented in this section because it is the base
of some of the more recent models in the literature, including the one presented in this thesis.
It relies on the assumption that expansion occurs when the reaction products occupy a volume
larger than the capillary porosity, which they estimate through a Powers’ type model.
Accordingly, expansion may be calculated as

146
 X = h(X p − Φ c ), Xp > Φc
(5.5)
X = 0, Xp ≤ Φc

w / c − 0.39α
with Φ c = f c ⋅ , w / c > 0.39α (and 0 otherwise) (5.6)
w / c + 0.32
where X is the fractional expansion, Φ c is the capillary porosity of the concrete (fc=volume
fraction of cement; w/c=water/cement ratio; α=degree of hydration), X p is the net fractional
increase of solid volume on reaction, and h is a constant introduced to account for the degree
to which the potential expansive volume, measured by (X p − Φ c ) , is translated into actual
expansion (which they estimated to be around 0.05, from observations by Atkinson & Hearne,
1990). The model implicitly assumes that reaction of sulfates to form ettringite, rather than
diffusion of ions, is the predominant process (i.e., if diffusion is much faster than reactions,
then the formation of products do not depend on the availability of the reactants but only on
the rate of reaction). This is in contrast to the model proposed by Atkinson & Hearne (1990).
Note that expansion occurs in this model only when the reaction products (they considered
formation of ettringite as well as gypsum to be the reaction products, although they concluded
that gypsum does not contribute much to expansions) fill the entire volume of capillary
porosity, in contrast with some other proposals (Tixier & Mobasher, 2003a). To compute the
net fractional increase of solid volume X p they proposed the following relation

X p = b ⋅ r ⋅ ys ⋅ v ⋅ α r /(1− ∈) (5.7)
In this expression, ys is the volume of concrete that can participate in expansion (typically the
volume of C3A), v is the expansion factor (expressing how much larger is the solid product
volume than that of the solid reactants), α r is the degree of reaction (from 0 to 1 for complete
reaction), r is the ratio between volume of capillary pores and total pore space (including gel
pores), ∈ is the porosity of the product phase (e.g. space between ettringite crystals) and b is
the reactant ratio, defined as the total volume of solid reactant to the volume of the rate
limiting solid reactant. Note that already in this early model, the ‘porosity’ of the reaction
product was considered to be an important factor. If figures 5.5 and 5.6 are observed, it can be
seen that the reaction products filling the pores, with tiger stripe morphology, have a large
portion of free space, even in the cases where air voids have been totally filled with ettringite.
This model was later included into a simplified computer program (CONCLIFE) from NIST
(freely available), to predict the deterioration by sulfate attack and freeze-thaw of bridge
decks and pavements, assuming that ions enter only by sorption and not by diffusion (see e.g.
Bentz et at., 2001).
A complete and critical review of empirical and phenomenological models dealing with
sulfate attack may be found elsewhere (Skalny et al., 2002, Chapter 7).
5.2.2.2. Advanced chemo-transport-mechanical models
In recent years, some advanced chemo-transport-mechanical models have been proposed to
study the problem of external sulfate attack in cementitious materials (Schmidt-Döhl &
Rostásy, 1999; Planel, 2002; Tixier & Mobasher, 2003a,b; Shazali, 2006; Bary, 2008). The
common feature in all of these models is that they attempt to quantify the mechanical
consequences of the degradation processes on the basis of a chemo-transfer calculation of the
main species involved and the relevant (i.e. expansive) reaction products formed (i.e.
ettringite, gypsum, or a combination of the two). However, only some of them perform a full

147
tenso-deformational mechanical analysis (Planel, 2002; Tixier & Mobasher, 2003a; Bary,
2008).
a) Model by Krajcinovic et al. (1992)
The theoretical model proposed in the early 90’s by Krajcinovic and coworkers
(Krajcinovic et al., 1992) based on micromechanics of brittle solids is worth to be mentioned.
Some of the hypotheses introduced in this early model have been adopted in more recent
proposals (e.g. Tixier & Mobasher, 2003a). The chemo-transport analysis of external
magnesium sulfate attack in 1D (implemented with the FEM) is performed with the use of a
single-ion diffusion-reaction equation in terms of the ingressing sulfate ions (denoted as
c[mol/m3]) that react with the unhydrated C3A phase to form ettringite (denoted as
ce[mol/m3]), as follows
∂c
= ∇ ⋅ (D eff ∇c)- kc(c0a − ce ) (5.8)
∂t
where Deff [m2/s] is the effective diffusion coefficient (accounting for damage and
microcracking via micromechanical considerations), k is the (phenomenological) rate constant
of the reaction (taken as 1.35x10-5 m3/mol/s for the experimental fit), and c0a is the initial
concentration of C3A in the hardened cement paste. The assumption is made that ettringite
causes volumetric expansions ( ε t ), through a simple model accounting for ettringite growth
at the surface of spherical unhydrated C3A inclusions, which is also derived from
micromechanical considerations, yielding
ε t = f t β (ω )ε r (5.9)
in which f.t is the inclusion (ettringite) volume density, β (ω ) is a factor accounting for the
elastic properties of the inclusion (ettringite) and the microcrack density ( ω ), and ε r is the
volumetric expansion of the sphere (given by an expression obtained with the simple model
mentioned above). For the mechanical analysis, they adopted elastic-perfectly-brittle behavior
in tension for simplicity and introduced some concepts of LEFM and damage mechanics.
Satisfactory agreement with experimental expansions on mortar specimens with different C3A
contents (by Ouyang et al., 1988) were obtained with this model, although the extraction of
axial strains from a 1D analysis across the thickness is rather questionable. Its extension to
sodium sulfate attack and to 2D or 3D analysis would be of great value. Results obtained in
terms of axial expansions and sulfate concentration profiles are shown in figure 5.8. Note the
plateau of the sulfate profile at 240 days of exposure near the surface, which reflects the effect
of the cracking front (spalling) on the effective diffusion coefficient (included in the model).

148
Figure 5.8. Experimental (by Ouyang et al., 1988) vs. numerical results obtained with the
model by (Krajcinovic et al., 1992): (a) expansions vs. time for three mortars with different
C3A contents (5.3, 8.8 and 12%); (b) normalized sulfate profiles for three different times of
exposure (note the plateau at 240 days near the surface, reflecting the effect of cracking or
spalling on the effective diffusivity); and simplified model for ettringite growth around C3A
spherical particles.

The model has been recently updated (Basista & Weglewski, 2009), with the addition of
new features from other models, such as a relation between the diffusion coefficient and the
porosity (from Garboczi & Bentz, 1992, see also Chapter 6, section 6.1.2.) and a shift factor
for calculating the volumetric strains considering a fraction of the capillary porosity as a
threshold value (from Tixier & Mobasher, 2003, described at the end of this section). In
another contribution of the same authors (Basista & Weglewski, 2008), two forms of
chemical reaction of the ettringite formation, i.e. through-solution (in the pore solution) and
topochemical (solid-solid type), were considered and compared, leading to the conclusion that
the hypothesis of a through-solution reaction can not explain the overall expansions observed,
although the way in which crystallization pressures are applied to the macroscopic continuous
medium is questionable (see e.g. Flatt & Scherer, 2008).
b) Model by Shazali et al. (2006)
In the case of Shazali and coworkers (Shazali et al., 2006), they consider an empirical
relation between the degree of reaction of sulfate attack and the loss of strength, through a
relative strength loss function. In this way, degradation of the specimen is defined locally in
the finite element mesh, although a mechanical analysis is not performed. Coupling between
the chemo-transport model and mechanical degradation is introduced via a damage parameter,
also related to the degree of reaction through an empirical relation. In their work, the chosen
state variable for sulfate attack is gypsum, focusing on the modeling of the impact of gypsum
as a single form of sulfate attack, assuming that this product is the only responsible for
strength loss in concrete. Additionally, non-saturated conditions are taken into account
through a convection term in gypsum mass balance equation and the introduction of the
volumetric pore water content effect. The result is a diffusion-advection-reaction equation in
terms of sulfates concentration and water content in order to locally quantify the gypsum
concentration. Obviously, the consideration of non-saturated conditions adds a second mass
conservation equation for water, coupled through a staggered strategy to the previous
equation.
c) Model by Schmidt-Döhl & Rostásy (1999)
Schmidt-Döhl & Rostásy (1999) proposed a complete model based on an explicit finite
difference method for the transport processes, accounting for non-saturated conditions and

149
allowing diffusion of various solved species simultaneously. Additionally, two separate
modules account for the chemical equilibrium calculation, based on the minimization of the
Gibbs free energy of the system, and the reaction kinetics. The effect of chemical reactions on
the transport properties are not taken into consideration. The effect of cracks on transport
processes is considered in the case of the diffusion of solved species, although its accuracy in
the context of a finite difference scheme is rather questionable. The initial phase assemblage
considered includes CSH, monosulfoaluminate (AFm), ettringite (Aft), portlandite (CH) and
water. Multi-ionic diffusion is considered in the calculations (although neither the number of
ions nor their types have been specified in their publication). Results obtained in terms of
ettringite profiles of mortar cylindrical specimens (1D transport from an end face) are in
relative good agreement with their own experimental campaign. Expansion (measured as
changes in diameter of the cylinder at different depths) obtained by assuming ettringite as its
main cause agree well with experimental results. The change in material strength is
considered by assigning theoretical strengths to each of the solid components of the cement
paste (e.g. CSH phase is assigned a value of 365MPa, which seems somewhat arbitrary). The
expansions are computed when ettringite volume exceeds the existing pore volume, which is
in agreement with most of the proposals in the literature. Unfortunately, the mechanical
module of the model is not described in any of the available authors’ publications (nor are
other important aspects regarding transport processes). Figure 5.9 shows some results
obtained with the model and compared to experimental data of external sodium sulfate attack
(44 g/L Na2SO4) at fixed times of exposure, in terms of ettringite concentrations (at 154 days)
and expansions (at 154 and 302 days) at different depths (profiles). There is a satisfactory
agreement between experimental and numerical results in this case, although the authors
recognize that the model fails to represent correctly the degraded zone.

Figure 5.9. Experimental vs. numerical results obtained by Schmidt-Döhl & Rostásy (1999):
ettringite profile at 154 days of exposure to sodium sulfate solution (left), and expansions
measured as changes in diameter of the cylindrical mortar specimen at different depths (right).

d) Model by Planel (2002)

The model developed by Planel (2002) in the course of his doctoral thesis (see also Planel
et al., 2001) represents a progress in terms of mechanical modeling. A damage model that
accounts also for plastic strains, previously developed by Sellier and coworkers (Sellier et al.,
2001) has been used in this study. It is an orthotropic damage model with a rotating smeared
crack criterion and accounts for additional damage produced by precipitation of secondary
phases in a special way (Planel, 2002), although some strong assumptions are introduced
regarding the origin of tensile as well as compressive damage variables (tensors in this case).
The chemo-transport analysis has been performed with the commercial code HYTEC (van der
Lee et al., 2003), in which a chemical equilibrium with species and surface code (CHESS),
based on thermodynamic modeling (LMA approach), is coupled with a 1D transport finite
difference model (RT1D). The analysis considers isothermal and saturated conditions. The

150
influence of precipitation of reaction products on the transport properties may be accounted
for in the model. At each time step, the transport analysis is first performed and the results
enter the chemical equilibrium code which alters the transport process, and this same process
is repeated until a certain tolerance is satisfied in terms of concentration of products
(staggered approach). The main constituents of the cement paste considered in their
simulations are portlandite (CH), monosulfoaluminate (AFm), ettringite (Aft) and CSH (with
a ratio Ca/Si=1.65, although this compound is only used to determine the capillary porosity),
which may be expressed as a function of the concentration of calcium, silicon, aluminum and
sulfur, and the single species considered in the calculations are Ca2+, Na+, SO42-, Al+3, H+,
OH- and H2O. Typical results obtained with this model, together with the comparison with
their own experimentally-determined species profiles (gypsum, ettringite and portlandite)
from external sodium sulfate attack, are shown in figure 5.10(a). The same type of
distribution for the main solid phases has been found in another series of experiments for low
concentration sodium sulfate solution (Le Bescop & Solet, 2006). It can be observed that
although the chemo-transport model is more sophisticated than a single-ion transport model, it
fails to capture the distributions of ettringite and gypsum. Moreover, the mechanical analysis
is decoupled from the chemo-transport model. Thus, the effect of cracks in the transport
properties can not be taken into account, limiting its applicability to the onset of cracking
phenomena (Planel, 2002). In this case, ettringite as well as gypsum are assumed to cause
expansion (and damage) once a fixed fraction of the pore space (which in turn evolves) is
filled. Results of the mechanical analysis in terms of damage distribution are included in
figure 5.11(b), together with the experimental crack pattern observed in figure 5.11(a).

Figure 5.10. (a) Experimental (by XRD technique) vs. numerical profiles of gypsum (gypse in
French), ettringite (Aft) and portlandite (CH) obtained by Planel (2002). The material is a
cement paste with w/c ratio of 0.4 and type I cement and the time of exposure is 84 days (12
weeks) and the specimen is cylindrical (φ of 70mm, height of 8mm). It may be seen that the
approximation is rather crude (from Planel, 2002). (b) Results of the simulation of the same
experiment at the same time of exposure (Planel, 2002) with the model by Bary (2008). In this
case, however, the profiles have been numerically fitted, as opposed to the previous model.

e) Model by Bary (2008)

More recently, within the French Atomic Energy Commission (CEA), a different chemo-
transport-mechanical model also for isothermal saturated conditions has been developed to
analyze sulfate attack with combined leaching of cement paste (Bary, 2008). The chemo-
transport model takes into account two ionic diffusion processes (i.e. two mass balance
equations): calcium ions (which are assumed to completely define the degradation state of
cement pastes due to leaching), and sulfate ions (assumed to define the sulfate migration as

151
well as gypsum and ettringite formation in an intelligent way). As input data for the chemical
model, the initial concentrations of the solid phases are required. Similar to the previous
model (Planel, 2002), the phases considered are CSH, portlandite, AFm, AFt and gypsum
(also initial porosity has to be provided). In this way, the model is able to calculate the
distribution of portlandite, gypsum or ettringite, as in the model by Planel (2002), although in
a very different way (Bary, 2008), without performing chemical equilibrium calculations.
Indeed, since the reaction rate coefficients considered are of phenomenological nature, they
have to be determined by back analysis in order to fit gypsum, Aft and CH profiles at a fixed
time to then check the model’s response at any arbitrary time. The results obtained with this
model, in terms of distribution of phases for the fitted time (which should be compared with
results in figure 5.10(a) by Planel), are shown in figure 5.10(b). The effect of damage on the
diffusion process is not taken into account in this model, although reduction in porosity due to
secondary phases’ precipitation (and increasing porosity due to dissolution of initial solid
phases) is accounted for. Thus, as in Planel (2002), simulations are valid only until crack
initiation, when the effect of damage on diffusivity becomes non-negligible. The mechanical
model is based on a poroelasticity analogy, where the traditional pore pressure takes the form
of crystallization pressure exerted by ettringite (resulting from the interaction between
growing AFm crystals from supersaturated pore solution and the surrounding CSH matrix)
and the Biot coefficient translates into an interaction coefficient αAFm (Bary, 2007). The
damage model of Mazars is adopted for crack representation. The input parameters are obtain
from a homogenization technique (with the use of the Mori-Tanaka scheme), considering
CSH, portlandite, AFm, ettringite, gypsum and non-hydrated compounds as main solid
phases. The chemo-transport-mechanical problem is coupled through a staggered strategy in
Cast3m FE code. This model seems promising and the assumptions made therein allow
simplifying the chemo-transport model to two relatively simple mass conservation equations
for calcium and sulfates. However, its application has been tested only for the case of thin
cement paste samples, which is still rather limited. In addition, the simplicity of the adopted
damage model, with its limitations, is in contradiction with the sophisticated framework of
poroelasticity and homogenization technique included in the model. The not very realistic
damage patterns obtained may be related to this fact, as shown in figure 5.11(c).

152
Figure 5.11. Comparison of experimental crack patterns with numerical results obtained with
the models by Planel (2002) and Bary (2008). (a) Crack patterns at 30 (left) and 40 (right)
days for a cement paste specimen of 1.5x4x16cm3 (w/c=0.6) immersed in 15mmol/L sodium
sulfate solution (from Planel, 2002 and Bary, 2008). (b) Distribution of mechanical damage
(left) and chemical damage (right) at 95 days over the deformed mesh for the model by Planel
(from Planel, 2002); (c) damage distribution in the deformed mesh at 36 (left) and 50 (right)
days for the model by Bary, for the simulation of the same experiment (from Bary, 2008).

f) Model by Tixier & Mobasher (2003)

Finally, the single-ion diffusion-reaction model proposed by Tixier and Mobasher
(2003a,b) considers a simplified view of the problem, treating sulfate attack in a more
phenomenological way than the previous two models, as in the case of Krajcinovic et al.,
1992 and Shazali et al., 2006, although with very different underlying mechanisms. Given the
complexity of the problem, they assumed that external sulfate attack can be analyzed with
only one diffusing ion type, represented by the ingressing sulfates, which react with the non-
diffusing alumina phases of the hardened cement paste, present in the form of
monosulfoaluminates, tetracalcium-alumina hydrate (C4AH13) and unreacted C3A, to form
ettringite. This last compound is assumed to be the only reaction product governing
expansions. In this model, the chemo-transport analysis is performed simultaneously via a
diffusion-reaction equation, assuming a constant rate of reaction coefficient. Additionally,
they consider a second equation that completes the chemo-transport model that accounts for
the alumina-phase depletion, which depends also on the rate of reaction and a stoichiometric
coefficient. The model assumes that gypsum is formed first following
Ca(OH) 2 +Na 2SO 4 ⋅10H 2O → CaSO 4 ⋅ 2H 2 O (5.10)

153
 In a second step, the model assumes that the three compounds mentioned above may react
with ingressing sulfates, represented in the form of gypsum, according to the following
reactions
C4 AH13 + 3CSH 2 + 14H → C6 AS3H 32 + CH
C4 ASH12 + 2 CSH 2 + 16H → C6 AS3H 32 (5.11)
C3 A + 3CSH 2 + 26H → C6 AS3H 32
For any of the individual reactions described above, the volumetric change due to the
difference in specific gravity can be calculated using stoichiometric constants. The above
equations on the formation of ettringite may be lumped into one single equation as
CA + qS → C6 AS3H 32 (5.12)

where CA=γ 1C4 AH13 + γ 2 C4 ASH12 +γ 3C3 A , being γi the proportion of each phase, and
q=3γ 1 + 2γ 2 +3γ 3 represents the weighted stoichiometric coefficients of the sulfate phase.
The preceding chemical reactions take place according to the sulfates and calcium
aluminates availability, which is determined in time and space through a diffusion-reaction
equation for the sulfate concentration (denoted as U[mol/m3 of material]) and an additional
equation for the depletion of calcium aluminates (denoted as C[mol/m3 of material]):
∂U ∂2U
= D 2 - kUC (5.13)
∂t ∂x
∂C UC
= -k (5.14)
∂t q
in which D[m2/s] is the diffusion coefficient, assumed to be dependent (increase) on cracking
through a damage variable, k[m3/mol/s] is the rate of take-up of sulfates constant, and t[s] and
x[m] are the time and space coordinate, respectively. The effect of pore filling on the
diffusivity is not taken into account in this model. The volumetric strain εV(t) is obtained from
the amount of reacted calcium aluminate and the volume change associated with it. An
averaging scheme is again used for the three phases, yielding
ε V (t) = α s ⋅ CA reacted − f ⋅ Φ (5.15)
In the previous equation, CAreacted represents the reacted calcium aluminates (difference
between initial and non-reacted calcium aluminates), αs is calculated as a sum of the
volumetric changes of each reaction considered, Φ is the capillary porosity (estimated through
a Powers’-type model) and f is the fraction of capillary porosity that has to be filled before
any expansion occurs.
For the mechanical analysis they adopted a simple damage model previously proposed by
Karihaloo (1995) for the description of concrete under uniaxial tension. The damage variable
is arbitrarily defined as a function of strains. More details on the mechanical model may be
found elsewhere (Tixier, 2000). The resolution of the model’s equations is done through a
finite difference scheme in 1D (with Matlab) and extended to 2D with the use of the
superposition method (yielding only approximate solutions). Moreover, results in terms of
expansions are obtained by extrapolating the strain field in a 2D section to the out of plane
direction (Tixier, 2000). Chemo-mechanical coupling is accomplished with the use of the
moving boundary problem approach (Tixier & Mobasher, 2003a). Unfortunately, the

154
assumption made in order to reduce the system of eqs. 5.13 and 5.14 to one differential
equation is not correct. They proposed to define a new variable Z = U – qC, and since no
diffusion of calcium aluminates is possible they ended up with the simple expression
∂Z ∂2Z
=D 2 (5.16)
∂t ∂x
& = - kUC (dot being derivative with respect to time) into
However, a simple substitution of qC
equation 5.13 yields

& = D ∂ U - qC
2
U & (5.17)
∂x 2
Rearranging terms one finally obtains
∂Z ∂2U ∂ 2 (Z + qC)
=D 2 =D (5.18)
∂t ∂x ∂x 2
Clearly, since the diffusion coefficient of sulfate ions (D) is the one previously defined
(and is certainly not zero), a gradient of calcium aluminates (C) will have an effect on the
variation of Z. In other words, a diffusion of the alumina-phase has been erroneously
introduced.
Typical results in terms of expansions obtained with this model are shown in figure 5.12.
Relatively good agreement is found, although the underlying assumptions introduced to
calculate these expansions from 1D calculations yield only approximate solutions.

Figure 5.12. Comparison of experimental and numerical results obtained with the model by
Tixier and Mobasher (2003a). Expansions evolution for a mortar sample tested by Ferraris et
al. (left) and for two cement pastes with different C3A contents tested by Ouyang et al.
(right). Note the stair-type curves due to the discretization adopted in their simulations (from
Tixier & Mobasher, 2003b).

5.2.2. Some final comments on the modeling of external sulfate attack

A complete review of, at least to the author’s knowledge, all of the relevant numerical
models dealing with sulfate attack has shown the different approaches that have been
followed in the literature in order to assess the degradation of cementitious materials when
exposed to a sulfate-rich environment. All of these models treat the material as a continuous
and homogeneous medium, and some of them propose a more detailed representation of the
chemical reactions occurring. The common feature in all of these models is the
phenomenological treatment of expansions as a result of the formation of reaction products

155
(ettringite and/or gypsum). This is not surprising since it is the most controversial issue also
from an experimental point of view. Only the model by Bary (2008) is capable of performing
fully two-dimensional coupled analysis in a FE environment. But even this model lacks an
explicit introduction of the dependency of the transport processes on the damage or cracking
level (in fact, the analyses made by the author were uncoupled). The rest of the models
consider different techniques to extract expansions or degradation profiles from 1D
calculations. Although the use of chemical equilibrium models in order to obtain the phase
assemblage at equilibrium seems a very promising tool in fully coupled models, results obtain
so far show that the resulting mechanical response does not capture the correct behavior in
terms of spalling and crack patterns, and expansion evolution has not been reported at all in
these cases. Finally, it should be mentioned that none of the models proposed so far are able
to capture the main crack patterns correctly (this is logical since the mechanical analyses in
these proposals are based on simplified damage models), nor are they able to account for the
effect that crack and spalling may have on the transport processes (instead, only some of these
models propose an arbitrary increase in the diffusivity as a function of damage). This is an
important feature, since in the cases where spalling occurs, a drastic change in the boundary
conditions of the chemo-transport analysis should be expected, thus accelerating the
degradation process. The model proposed in this thesis will address these issues within the
framework of a meso-mechanical approach with discrete cracking.

156