Microporomechanics study of anisotropy of ASR under loading

Abstract In this article, we introduce a new micromechanical model for alkali–silica reaction. Our idea was to build a model with the following characteristics. First, the model has to be simple enough to be used to compute damage under loading and chemical attack at the level of each element in a structure code. Second, its parameters must be easy to identify on available alkali–silica reaction lab experiments. We have chosen to model the behavior of concrete containing aggregates such that most of the damage occurs in the cement paste. Using micromechanics and an energy criterion, the model remains analytical except for the minimization of the energy. The parameters were identified on Multon's triaxial experiments and good results were obtained for compressive loadings up to 10 MPa.


Introduction
Concrete is vastly used for construction and the life span of buildings is an important parameter of their rentability. Therefore durability issues such as freeze-thaw, carbonation or endogenous reactions need to be addressed. The alkali-silica reaction is an endogenous reaction discovered in the US in the 40's by Stanton [66]. It can develop several years after building. The reaction is externally visible because of a characteristic surface cracking which tends to align with the main compression directions and also through reaction products which sometimes flow from these cracks. The reaction also induces irreversible expansions which can lead to complications when using the structure. Various kinds of buildings can be affected. They all share some characteristics. First, their aggregates must contain some reactive silica, in the sense that it can dissolve in the concrete interstitial solution. Second, the concrete must contain alkali, usually provided by the cement. Third, there must be a water supply to the structure, since water enhances the transport of chemicals inside the concrete and is absorbed by the alkali-silica gels, leading to their swelling.
How to stop this reaction is an open problem. Today, a lot of effort has been made to make sure that no new building will be affected. For affected structures, mechanical techniques exist to release the stresses due to the reaction or try to stop the expansion mechanically. However, these solutions are expensive. It is therefore necessary to have simulation tools available to estimate the remaining life-span of affected structures and the efficiency of repair ASR modeling can roughly be divided in three branches. First, chemistry and transport models studying at the aggregate or structure scale, the transport of chemicals and water, so as to predict the amounts of gel produced at reactive sites, that is close to or in the aggregates. These models directly simulate the chemical reaction with various simplifying assumptions [59,58]. Second, mechanical models at the structure scale which are usually implemented in more general concrete durability codes accounting for various important phenomena such as creep, shrinkage and macroscopic damage [45,2,63,33,17]. Finally, microscopic scale models which can be analytical or numerical and try, by a fine description of the microstructure, to determine the mechanical consequences of the presence of a swelling gel in the concrete porosity close to the reactive sites [44,16,22,64,49,4]. The aim of these models is at the same time to help studying the affected structures and help understanding the physics of the degradation. It seems to us that analytical models are more likely to help for structure computations since their simplicity makes more obvious the influence of their parameters. Also, we can hope that such a model could serve as a way to compute degradation in each element in structure-size codes. That is why we have focused our research and this article on such models.
Therefore, for our model to be satisfactory, is must be such that its parameters are easily identifiable in the lab on concrete samples, some by direct testing, some by fitting of expansion curves. The model must be simple enough to keep short computation times if we want to couple it with a structure-size code. It must also reproduce expansion of the expansion at intermediate load, and reduction of the total (volumetric) expansion at high loads.
The most advanced effort in understanding the effect of an external load on the mechanical consequences of ASR seems to be the numerical model of Dunant and Scrivener, presented in [22]. In his second paper, the focus is placed on anisotropy [23]. The author gives interesting experimental results of expansion under uniaxial loading. The cracking/damage pattern in the aggregate (where the gel is located, in pockets) and in the cement paste is influenced by the loading. However, the model fails at high loads due to an artificial coalescence of cracks in 2d. This model was recently improved by Giorla et al. [31,30] who introduced creep in the cement paste. The effect of particle shape on the anisotropy of free expansion was successfully reproduced as well as many features of ASR. The effects of creep and damage (using different criteria) are investigated. This model is very advanced but few comparisons with experiments are available at the time.
Finally it seems the most complete experimental work concerning the expansion of concrete samples under loading is Multon's [50,51], since a variety of compressive loads are used (0 MPa, 10 MPa, 20 MPa) for the loading direction, in conjunction with a restriction of the radial expansion by steel rings of two different thicknesses (3 mm rings, 5 mm rings). This set-up induces a triaxial stress state which influences cracking in a more complicated way than the usual uniaxial compression [39,23]. Giorla also developed a new device for ASR under triaxial load which is very promising [31].
Our goal is, as stated in our introduction, to improve the micromechanical modeling of ASR. We will present a model based on micromechanics, poromechanics, and an energy fracture criterion. Therefore, we will start our presentation of the model by explaining the thermodynamic framework we used for our model.

Thermodynamic evolution of a concrete undergoing ASR
In this section we will present the basis of our model. We will start by some geometrical simplifications of the description of the concrete. Our model has been built in the idea of studying concrete made of aggregates which are sometimes classified in the literature as fast reacting aggregates [32,56]. For us the important point is that we focus on aggregates in which the expansive products can easily reach the interface with the cement paste. Therefore the most important cracking phenomenon are considered to be the decohesion between the aggregate and the interface first, and in a second stage the propagation of cracks in the cement paste. Some authors focus on slow reactive aggregates, in which the most important cracking mechanism occurs inside the aggregates [21,22,62,61,31]. This category of aggregates is of great interest, and the methodology developed in this article will hopefully be used in the future for this kind of aggregates.

Microscopic description of a concrete undergoing ASR
In this section we introduce the geometrical description of our attacked concrete at a given time t. We will call Ω the actual volume of concrete that we study.
3.1.1. Description of the concrete under attack before cracking Our concrete, at the scale at which we consider it, is composed of two main components: grains and cement paste matrix, whose properties vary due to the attack. Let us detail their properties, beginning with the grains. They are an ensemble of aggregates and sand particles (which size are extremely variable, from microns to a few centimeters). Let us call their number N g in the volume Ω (g stands for grain). All grains are considered spherical. The grain i ∈ N g is of radius R i , and its current attack degree is called α i (t). The attack degree represents the proportion of the radius of the grain which has undergone attack. Therefore the shell between radii (1 − α i (t))R i and R i will be called the attacked zone, while the sphere of radius (1 − α i (t))R i will be called the sound zone, at time t (Fig. 1). The sound zone has the mechanical properties of  a sound aggregate. It is therefore taken as a linear elastic isotropic material, which stiffness tensor is called C a . The attacked zone has been partly dissolved by the attack of ions coming from the cement paste. We do not describe this attack but only its mechanical consequences. We assume that it transforms the sound aggregate into a porous material of porosity ρ i , whose poroelastic properties are a tensor of elasticity C i p , a Biot coefficient b i p and a Biot compliance M i p , which are defined according to the habits in poromechanics as presented by Coussy [18], except for the Biot compliance, which is taken as the inverse of the usual Biot modulus. We get the following constitutive law for the attacked zone of grain i: Where σ is the stress tensor,φ i is the deformed (due to strains ε and pressure p) porosity of the attacked zone. We assume that the attack keeps the isotropy of the aggregate, therefore the Biot coefficient b i p is a scalar. In our model the porosity ρ i can be different in different aggregates, but is constant in time. The attack progresses by increasing of the size of the attacked zone (described by α i (t)) only. The poroelastic coefficients of this zone are obtained using the Mori-Tanaka estimate [5], assuming the porosity is composed of spherical cavities. This simple homogenization step is not described here.
The grain is surrounded by an interfacial transition zone (ITZ ). We think it is important to take it into account for two reasons. First, attempts to predict the stiffness of concretes have shown that if micromechanics models are used considering concrete as a two-phase (aggregates and cement paste) material, the mechanical properties were overestimated. Nielsen and Monteiro [54] argued that it is better to consider three phases (aggregates, ITZ, and cement paste) to predict the mechanical properties. It confirmed mechanically observations that were made long before by Farran for example [25] about the modification of the packing of cement grains close to the aggregates leading to different composition and porosities in this zone. There is discussions about the dependence of the properties of this zone on the total aggregate volume fraction [35] based on generalizations of the self-consistent scheme of Christensen and Lo [60,15], or on the aggregate size [36], and whether it is possible to model it as a uniform zone or not [36,10,53]. The second reason why we need to take the ITZ into account is related to ASR. The ITZ plays the role of a reservoir for gel escaping from the aggregate, limiting the pressure increase. Not willing to put too much detail in the ITZ which is not the focus of our study, we choose to model it as a homogeneous porous medium of constant thickness l c and porosity ρ itz (both independent on the grain considered). The properties are written in the same manner as for the attacked zone: a tensor of elasticity C t , a Biot coefficient b t and a Biot compliance M t in the following constitutive law: Whereφ itz is the deformed (due to strains ε and pressure p) porosity of the ITZ. The poroelastic coefficients relative to the ITZ are also determined using the Mori-Tanaka estimate. Now that we have described in detail the geometry of our concrete during attack, we will proceed with the description of cracking.

Description of cracking
The concrete described in the previous paragraph is going to change morphology due to the high pressures which will develop in the pores of the attacked zone and ITZ. Let us describe the changing microstructure by a number of damage parameters for each grain i. First, a decohesion parameter d i (t) which can take values of 0 and 1 only and is an increasing function of time. At the beginning of the attack and as long as the aggregate remains bounded to the cement paste, d i (t) = 0. When the interface at radius R i , that is between the attacked zone and ITZ breaks, d i (t) = 1. Second, crack size parameters we note collectively x i (t) which describe the size of the cracks developing in the cement paste. The choice of crack shape and exact meaning of these parameters will be explained in § 4.2, since some approximations will be related to specificities of the micromechanical description. The crack sizes are also increasing functions of time.
The geometrical properties of our sample has been described and we have defined the damage parameters accounting for the evolution of the microstructure. Let us sum things up defining the problem we want to solve: (P): Find the evolution of the damage parameters d i (t) and x i (t) due to the ASR, represented by the set of attack degrees α i (t).
To simplify the resolution of this problem, let us introduce aggregate families.

Definition of the aggregate families
An aggregate family is a group of aggregates which are assumed to have identical evolutions. Usually it will correspond to a given aggregate size, but a given aggregate size can also be divided in families. Let us call their number N f . We assume that the grains of family i are of size R i . The volume fraction of grains of site i is called f i , and the number of grains per unit volume in family i writes Aggregates which are in the same family will follow the same evolution in terms of attack and cracking. The vicinity of each aggregate contain various kinds of pores. First, pores in the aggregate which appear due to the chemical attack. As explained in § 3.1.1, they are located in the attacked zone between radius 1 − α i R i and radius R i , and represent a volume fraction ρ i of the region where they are located. That means that in the family i, there is a total porous volume in each aggregate equal to ρ i 4π The overall constitutive law has been written in Eq. 1. Second, pores in the ITZ which, due to its small thickness l c occupy a volume ρ itz 4πl c (R i ) 2 for each aggregate. This volume is available to the gel, which we think has threshold behavior and because of that, cannot penetrate farther in smaller pores of the cement paste [14]. The constitutive law of this zone was given in Eq. 2. Third, if decohesion has occurred, the gel is also able to occupy the space between the aggregate and the ITZ in each site, and if cracking has occurred, it can go in the cracks. These zones have no undeformed volume, but gel under pressure can invade them.
In our model, we assume that all pores belonging to the same family are at the same pressure that we call p i , which means that we consider that the alkali-silica reaction develops slowly enough so that the gel can flow to every part of the porosity without dissipating large amounts of energy.
To solve problem (P) in various cases (for example for different macroscopic loading conditions on the sample), let us introduce a virtual problem.

Virtual problem: loading with imposed strains and attack degree
As explained in the previous section, the chemical attack is described by a number of attacked degrees. Let us now take external loadings into account, since one of the goals of our model is to be able to estimate expansions of an attacked concrete under loading.

Total energy function
Let us now describe problem (P * ) which will give an energetic framework for solving problem (P). Let us first be clearer about the loading parameters of our problem. Our problem is driven by a set of time-dependent parameters α i (t) which represent the attack degree for each site present in the sample. The second loading parameter from now on will be the macroscopic strain E. That means the displacements E.x are prescribed on the external boundary of Ω. This is a particular external loading which is convenient for the presentation, but our model is suited for any macroscopic loading as will be shown in the examples section at the end of this article (section 5). Now that our loading parameters are defined, let us choose a given loading state E, α i and write formally the total energy of our sample for a set of virtual damage parameters d i * , x i * . This total energy is the sum of the elastic energy stored in all components of the system, plus the surface (or dissipated) energy due to crack creation: Where the loading parameters appear as parameters for the energy functions, while the damage parameters appear as arguments. The values of these energy functions will be specified later. Let us now explain the energy criterion used to compute the crack state of the sample for a given loading.

Postulate: law of evolution of the damage (Francfort-Marigo criterion)
The damage law we choose is as follows: the real damage parameters are found by minimizing the total energy over all possible virtual damage parameters, as proposed by Francfort and Marigo, Fedelich and Ehrlacher, or Mielke [28,26,47]: Where A is the set of admissible damage states. Without entering too much details about this, it requires that damage is increasing, and that cracks surfaces are not interpenetrating.
These definitions of the energy and evolution law form problem (P * ). We now need to specify the shape of the total energy, which means mainly of the elastic energy. However it is a difficult problem since the attack degrees influences the energy in a complex way. Let us introduce a simpler problem where we assume pressure is controlled.

Auxiliary problem: loading with imposed strain and pressures
We do so by means of an auxiliary problem (P aux ) to the virtual problem (P * ). In this problem, the loading variables are changed. In problem (P * ), the loading parameters are E, α i , corresponding to the physical parameters imposed in our representation of alkali-silica reaction. While we apply a strain and the chemical attack progresses, all regions filled with gel (pores in the attacked zone, pores in the ITZ, cracked aggregates interfaces, cracks in the cement paste) in family i are under pressure p i .

New loading parameters and resulting strain
For the purpose of writing explicitly the elastic energy, we decide to use E, p i as loading parameters for the mechanical auxiliary problem. Following the description we have proposed of the morphology of the attacked concrete, there are two type of sites we need to distinguish, which we will call Type I and Type II, as shown on Fig. 2 of Type I are those where decohesion has not occurred yet. Therefore the porosity is only that of the attacked aggregate and of the ITZ. The sites of Type II have undergone decohesion. Therefore, from the point of view of the solid skeleton, the aggregate is not visible anymore. Only the pressure on the cavity and the crack lips remain. The N f families are therefore divided into N I families of Type I and N II families of Type II.

. The sites
In this new mechanical problem the displacement and hence the strain is, everywhere in the structure, proportional to the loading parameters because all materials are linear elastic. Let us define two sorts of localization tensors. First, a strain localization tensor A(x) such that if all pressures p i are equal to zero, the strain in the structure writes: It is a fourth order tensor, with the minor symmetry but not necessarily major symmetry. Second, As a result our strain writes, for any set of loading parameters in the auxiliary problem: 3.3.2. Elastic energy of the solid skeleton and constitutive law We are now able to write, at least formally, the expression of the elastic energy of the skeleton by unit volume of the porous medium in the auxiliary problem for a virtual damage state, using a poromechanical description at the scale of the REV. We make the strong assumption that the matrix of Biot compliances is diagonal, that is when pressure is applied in a given porous zone, the average strain on this zone is much larger than anywhere else in the porous medium. It makes the expressions simpler and avoids keeping terms that cannot be approximated easily by micromechanics tools. Let us illustrate this assumption by showing a numerical evaluation of the value of the crossed Biot compliances in a very simple case : the case where we have 2d plane strain elastic medium with circular pores, gathered in two families, each with volume fraction f . These results are unpublished but were established by the authors while comparing finite element results with micromechanical estimated. This work is already published in [13], where all details concerning the simulations can be found. The materials parameters chosen were E = 1 for Young's modulus and ν = 0.25 for Poisson's ratio. As one can see on Figure 3, the proportion between diagonal Biot compliances M 11 = M 22 and the extra-diagonal ones M 12 = M 21 greatly varies with volume fraction. When reaching f = 0.2 for both families, which means a total porosity of 0.4 and is close to what can be reached with the random sequential addition algorithm we used for these simulations, crossed and diagonal Biot compliances are close, with opposite signs. This means that at zero macroscopic strain, the pressure in one family induces an increase of its volume of just a little larger amount than the decrease of volume of the other family. This fact is surely to be considered to precisely study the coupled evolution of damage around different reactive sites. We underline that the use of the strain localization tensor to determine the Biot compliance which is made in this paper induces the replacement of M ii by j M ij and M ij by zero ( § 4.1). However we stick to this assumption which is exact if the various families undergo exactly the same evolution, even if obviously we loose a coupling phenomenon by doing this. We keep in mind that this could be an idea for an improvement of the model, which will be however hard to implement due to the lack of simple micromechanical formula.
While the constitutive equations of the overall porous medium write: Where B i is the overall Biot coefficient of family i, describing how a pressure in a pore creates a macroscopic stress Σ in a macroscopically restrained sample, or equivalently how a macroscopic strain induces a pore volume change in a pore free of pressure, M ii is the overall Biot compliance of family i, linking the pressure in a pore family to its volume change when the macroscopic strain is zero, and C hom is the macroscopic stiffness tensor. The computation of these poroelastic properties will be detailed in section 4. The dual quantities of the imposed deformation E and the imposed pressures p i are the average stress on the porous medium Σ and φ i −f i which is the difference between the volume fraction of pores in site i relatively to the volume of porous medium in deformed configuration φ i and the volume fraction of pores in site i relatively to the volume of porous medium in undeformed configuratioñ f i , that means the volume variation of the porous space of site i relatively to the volume of porous medium.
We have completed the presentation of the auxiliary problem, let us explain how it is related to the virtual problem.

Back to P * , introduction of the gel
In this section we wish to compute the pressures in problem P aux , and introduce them in the energies given in the previous paragraph to yield energies of problem P * .

Undeformed and deformed gel volume
First let us write the undeformed gel volume, from the attack degree α i , of site i. We define the parameter δ that we call the expansion factor of the gel. It is the volume of gel created by unit volume of dissolved aggregate. The volume fraction of new porosity due to the attack at site i per unit volume of concrete writes: While the undeformed volume fraction of the ITZ of family i writes: Since we assume that the gel has a volume δ times that of the dissolved aggregate, the total undeformed volume of gel of site i writes: Then assuming our gel behaves elastically and is characterized by a bulk modulus K g , the deformed volume of gel in sites i per unit volume of porous medium V i writes: This gel volume has to be compatible with the volume available.

Available volume
This volume of gel has to be compared to the available volume in deformed configuration for site i by unit volume of porous medium. For families of Type I, it simply is equal to the current volume fraction of the pores φ i , which according to the constitutive law at the macroscopic scale writes: Where the undeformed volumef i = V i 0 + V i itz is the sum of the newly created porosity and the porosity of the ITZ and the deformed porosity is occupied by gel if there is enough gel φ i = V i . Therefore for families of Type I, the equation we obtain is: Which can be solved for the pressure as: In families of Type II, the undeformed porosity in the constitutive equation But the deformed porosity is occupied both by gel and a deformed a , and the pressure writes: Where K a is the bulk modulus of the sound aggregate. In Eqs. 16 and 17, F i , B i and M ii depend on the attack degree α i , and B i and M ii depend on the damage state. If Eqs. 16 and 17 predict negative pressures, they are taken as zero instead.
3.5. Total energy in problem P * The combination of energies of the problem P aux which are expressed in terms of E, (p i ) i=1:N f for a damage state d i * , x i * , combined with the expression of the pressure under the loading E, (α i ) i=1:N f at the same damage state, will now yield the energy in terms of imposed deformation and attack degrees, which corresponds to the situation of problem P * .

Contribution of the solid skeleton
Now that we know how to compute the elastic energy of the skeleton at given macroscopic loading and pressure (Eq. 8) and how to compute the pressure for a given macroscopic loading and given attack degrees (Eqs. 16 and 17), we can write the elastic energy of the skeleton as a function of our real loading parameters, for a given virtual damage state d i * , x i * :

Contribution of the gel
The elastic energy of the gel is written as follows for site i per unit volume of porous material :

Contribution of the aggregates
While for Type II families, in which the aggregates are no more accounted for as a part of the skeleton, the energy stored in the aggregates of site i writes:

Dissipated energy
It is the sum of crack surfaces weighted by surface fracture energies, which we take equal to G dec for decohesion and G f iss for cracks in the cement paste. Therefore it writes: We have shown the way we want to solve this evolution problem. One information is missing: how are the homogenized properties of the concrete related to the material properties of the various regions for a given damage state. We address this question in the next section.

Resolution technique: microporomechanics
We have all the tools to write the elastic energy of our attacked concrete for any virtual damage state d i * , x i * and any imposed strain and attack degrees E, (α i ) i=1:N f . We need to compute the poromechanical properties of the concrete appearing in the overall constitutive equations (Eq. 9). We will not detail this part but give final expressions. Details about this derivation can be found in [19,11,55].

Poroelastic coefficients
The expressions write: Introducing the volume fractions of each zone, α referring to Type I families, β to Type II families : Fig. 2). The homogenized stiffness tensor writes: Where, as shown on Fig. 2, C c is the elasticity tensor of the cement paste, C a is the elasticity tensor of the sound aggregate on domain Ω α a , (C p , b p , M p ) and (C t , b t , M t ) are the elasticity tensor, Biot coefficient and Biot compliance of the attacked aggregate and ITZ, respectively on domains Ω α p and Ω α t . < A > D denotes the volume average of the field A on domain D. The Biot coefficients of the various zones containing fluid write: Let us now use micromechanics to give approximations of the averages of the various localization tensors.

Approximation of the localization tensors
In our approach, one of the main objectives is to keep the model very simple. Hence, we have decided to approximate the localization tensors thanks to classical micromechanics schemes. We have decided to use a scheme called the interaction direct derivative (IDD ) scheme, which was introduced by Zheng and Du [68,20]. Let us first present a dilute estimate for our problem, since it is the easiest way to write and also serves as a basis for the IDD scheme.

Dilute approach
In this approach, we assume that the heterogeneities only occupy a very small volume fraction of concrete. Hence, to estimate the localization tensors averages on the various zones, we can use simpler mechanical problems, assuming that each grain lies in an infinite cement paste. Let us first show our approach for Type I families.
Type I sites. In these families which have not undergone decohesion, we have a number of spherical shells. At the center, the sound part of the aggregate in a sphere of radius (1 − α i )R i and modulus C a . Then, the attacked zones between radii 1 − α i R i and R i , and poromechanical properties (C p , b p , M p ). Next, the ITZ between R i and R i + l c with proporties (C t , b t , M t ). Finally, an infinite cement paste of elasticity tensor C c .
As explained in § 4.1, we need to compute the localization tensors under two types of loading: strain imposed at infinity on the one hand and pressure in the porous zones on the other hand. This is conveniently achieved using Love's solution [46] which is also used by Christensen and Lo in [15]. This solution yields the displacement field under an imposed shear strain. For the pressure and the spherical part of the imposed strain, it is much easier due to the spherical symmetry of the solution. We won't give the details of this solution which can be found in the author's PhD thesis and are classical [11].
The dilute estimate of the strain localization tensor is an isotropic tensor which can hence be written for the three regions as: Where J and K are the classical basis tensors for isotropic fourth order tensors and the index z takes values a, p and t. Second, the space averages of the pressure localization tensors can be determined simply using a spherical symmetry solution where the porous zones are under pressure.
Type II sites. For these families we haven't yet specified the shape of the cracked cavities. We need a very simple geometrical description so we can write the solution inspiring ourselves of Eshelby's solution. These sites are made of a spherical cavity on which the pressure of the gel surrounding the aggregate which has undergone decohesion is applied, and cracks in the cement paste which appear due to the high pressures remaining even after decohesion.
First, the cracks are assumed to be penny-shaped (PS) and centered on the center of the cavity. In fact, they have an annular shape. Each crack is only described by its orientation and external diameter (x). With this description it is possible to have N cracks of sizes x i , i ∈ 1 : N around each cavity, and orientations θ i , φ i . We want to reduce the number of cracks. Thanks to finite element calculations in 2d, we have shown that the overall properties of a medium containing randomly oriented cracks are very close to those of a medium containing cracks oriented along 2 orthogonal directions only [11,13]. Hence, we assume that anisotropic crack distributions can be represented by orthogonal cracks of adequate sizes. Therefore, we assume that in 3d, only 3 orthogonal cracks can grow around each cavity (see Fig. 4) This geometrical simplification leads us to try to compute the desired averages of strain localization tensors using a superposition of various Eshelby solutions. Our complex cracked cavity is hence considered as a superposition of a spherical cavity of Eshelby tensor S sph c and volume fraction f i , and 3 PS cracks identical to the cracks around the cavity, of aspect ration X f and sizes x 1 , x 2 , x 3 in directions e 1 , e 2 , e 3 . The associated Eshelby tensors only depend on the orientation, the aspect ratio and the outside medium (here the cement paste). We call these Eshelby tensors S f1 c , S f2 c , S f3 c . This decomposition was studied in detail in 2d in [11].
An important point in the choice of the volume fraction affected to the cracks (here we use volume fraction instead of the more common crack density parameter introduced by Budiansky and O'Connell [8]). If there is one crack in each direction τ for each family i of size x i τ and aspect ratio X f , the volume fraction for each direction and family writes: However, we have noticed that better results are obtained when correcting this volume fraction to account for the fact that the cracks intersect the sphere. Thanks to 2d computations detailed in [11], we came up with a correction volume fraction equal to: The final volume fraction we use for the crack in direction τ of family i is then: All the Eshelby tensors can be found in Mura's book [52]. Finally the dilute localization tensors corresponding to these families can be written: One could think that since the volume fractions of the cracks are much smaller than those of the cavity, the corresponding terms are negligible, but in fact some terms of the Eshelby tensors diverge when the aspect ratio goes to zero, so that their product by the volume fraction converges to a constant value [19].
We now have a dilute approximation of the averages of the localization tensors of Eqs. 22, 23, and 24. Let us explain how we obtain a better estimate with the IDD scheme.

Use of the IDD estimate
The IDD estimate was proposed by Zheng and Du [68,20]. From our point of view, it gives good results with cracks and is easy to use, so it makes it the best candidate to improve our estimate of the localization tensors [11,13]. With this estimate, it is easy to account for the diversity of inclusion shapes and space distribution. In many cases it is identical to classical estimates: the two-phase estimate [65] when we only have one type of ellipsoidal inclusion, Mori-Tanaka [5] when the inclusions are spherical. We have shown its effeciency thanks to 2d finite element computations on randomly generated microstructures, in cases where the interpenetration of cracks was forbidden [13]. Therefore, we don't expect our model to remain valid at high cracking states where percolation occurs. Also, this explains why we have considered the IDD estimate to perform much better than self-consistent estimates. Let us now build the IDD estimate from the dilute estimate.
Strain localization tensor. The IDD localization tensors are obtained from the dilute ones by multiplication by a correction term accounting for the interaction between the heterogeneities. It makes use of a Hill tensor describing the shape of the ellipsoidal inclusions P k c and another one describing the ellipsoidal spatial distribution of the inclusions P Dk c . Let us recall that the Hill tensor of a given inclusion k in medium c is related to the Eshelby tensor of the same inclusion in the same medium by the relation The correction term is the same for all inclusions. We call it T. More precisely, the following sum must be computed: However, this description only works for ellipsoidal inclusions, which is not the case for our Type I families. These cannot be fully described by Hill tensors, that is why we had to use the Love solution. If we go back to the idea behind this correction term, we can write it as a function of the dilute localization tensor: It writes: Then, if the inclusions are not ellipsoids but we know their dilute localization tensor (we do thanks to Love solution) and a Hill tensor approaching their space distribution, we know all the terms needed. Concerning the spatial distribution of the aggregates and the cracks, we keep things simple assuming that it is spherical (which might be wrong if the aggregates are flat, but then a lot of things would have to be adapted). We get: Therefore, we know how to write the IDD estimate for strain localization tensors. Concerning pressure localization tensors, we have not been able to make the same kind of modification and have decided to keep the same expression as for the dilute estimates.

Simple examples
We have explained in the previous sections the theoretical basis of our model. Let us now demonstrate how it works by very simple examples of concrete under attack. When using our model, the most time-consuming operation is the minimization of the total energy, which requires many evaluations of the energies for different crack configurations. Hence, to make things simpler, we need to decrease the number of crack sizes over which the energy is minimized. There are two simple ways to do that. First, reduce the number of aggregate families. Second, assume the symmetry of the crack state not to have 3 independent crack sizes for each aggregate family, but just one or two. In this section and next section we present our results, which were created using a second implementation of our program (compared to [11]), in Python, in the Materials Ageing Platform developed at EDF R&D [40].

Imposed strain
Our theoretical derivation has been performed in the case of an imposed macroscopic strain. In this example we assume that the strain is imposed to the value: Which corresponds to a moderate traction along axis e 1 . The result of this traction is that when the pressure develops around the aggregates, once decohesion has occurred due to high pressures, an anisotropic cracking takes place due to the combination of the pressure and the macroscopic loading, which corresponds precisely to the situation on real structures which we are modeling. Let us detail on this example how the degradation of the attacked concrete develops. We will choose simple val-ues for the parameters of the model: And, in this example, there will only be one aggregate family of size R 1 = 1 mm and volume fraction f 1 = 0.1.
Let us start with the observation of the evolution of the pressure. One can see on time during which no pressure increase is seen. This period corresponds to the filling of the ITZ by the gel. Then the pressure increases very fast to high values (60 MPa) in this case, but here the parameters have not been identified properly, so this value is not to be trusted. Following this increase, there are two pressure drops. These correspond to the decohesion of the aggregates and the beginning of cracking, with a slight crack size jump in this case. This interpretation of the evolution of pressure is confirmed by the observation of the crack size (Fig. 6). We see the crack jump at the same time as pressure drop. The crack size plateau which occurs around attacks depths of 4.10 −4 m is purely due to our model: using the Full Range IDD scheme, at some point the spatial distribution of the cracks is assumed to flatten. We haven't discussed this issue in detail here, but it is explained extensively in [57,3,13,11]. This evolution of the shape of the crack distribution induces this temporary crack size stabilization. It could be avoided by smoothing the relationship between the distribution and the crack volume fraction. Due to this increase in crack size in direction e 1 , there is a drop in mechanical properties in this direction, as can bee seen on Fig. 7. The first property loss is due to decohesion, which is why it affects all three directions identically. The second drop is only in direction 1 and is due to cracking. The Young's modulus of the skeleton diminishes to a very low value, while the undrained modulus, which takes the value of the gel and aggregate in the cavities into account diminishes much less (Fig. 8). We also display the Biot coefficients ( Fig. 9) and Biot compliance (Fig. 10). Their evolution are similar, but the Biot coefficient becomes anisotropic as the stiffness tensor, which is due to the crack orientation. The values of the anisotropic Biot coefficient might seem large (when it reaches 1, it means a stress inside the pore network is fully transmitted at the macroscopic scale, if strain is prevented), but are consistent with the Young's modulus values. When the Biot coefficient reaches 0.7 in direction 1, the drained Young's modulus is divided by almost a factor of 4. This large increase of the Biot coefficient ( Fig. 9) and decrease of the Young's modulus (Fig. 7) are due for the isotropic part to decohesion (spherical cavities of volume fraction 0.1 here yield a Biot coefficient of 0.2 with the IDD estimate), and for the anisotropic part, to a high concentration of cracks normal to direction 1, not in terms of volume fractions which are always negligible for cracks, but in terms of crack density parameter which is not computed here (see Fig. 6), [8]. Biot coefficient Figure 9: Evolution of Biot coefficients under uniaxial traction strain.
Finally, during this pressure build-up and crack creation under fixed displacement, stresses develop (Fig. 11). While at the beginning of attack the largest stress is in the direction of the imposed traction, cracking creates a high compressive stress which is most important in the cracking direction. It is also interesting to take a look at the evolution of energies during the evolution of the concrete (Fig. 12). The total energy is continuous, even at moments where the crack state is discontinuous. The dissipated energy jumps, and the elastic energy drops exactly of the same amount at decohesion and at the crack jump. The potential and elastic energies are identical be-  cause there are no imposed forces in this example. We have then shown that our model is able to determine the evolution of cracking depending on an external loading of the imposed strain type. Let us now give an example with imposed stresses.

Imposed stress
We choose to impose an unidirectional compression stress.
All parameters remain the same as given in Eq. 36, but the aggregate volume fraction is higher in this example (R 1 = 1 mm, f 1 = 0.3). In this example, a compression is imposed to the sample along the first direction. It is therefore more likely to see cracks appear along axes 2 and 3. As expected, this is what happens, as can be seen on Crack size  (Fig. 14). We have hence shown that our model is also able to deal with imposed stresses. The equations relative to this case have not been detailed because they are very close to those of the imposed strain case. The only difference is that when computing the total energy, the potential energy includes a term taking account of the work of the imposed stresses. Finally let us present a more complicated test which was used by Multon in his thesis [50].

Multon's test
In this test, cylindrical concrete samples are submitted to various compressive loads along their longitudinal axis, while their radial expansion is restrained by steel rings [50]. This set-up is interesting because of the complex three-dimensional stress state which is induced. The idea is that because of the longitudinal compressive load, cracking starts in the radial directions. This induces tension in the rings and eventually high stresses. Therefore, at some point radial cracking can be stopped and cracking can appear in the longitudinal direction where compression is still applied. This phenomenon of cracking reorientation is very interesting and a good way to test the qualitative behavior of our model. This behavior is shown here by using the same set of parameters as previously (Eq. 36), again with only one aggregate size (R 1 = 1 mm and f 1 = 0.3). As in Multon's test, the sample's radius is 6.5 cm and the steel ring's Young's modulus is 193 GPa. The rings cover all the height of the sample but are not connected. Hence they only impact the radial stress in the sample. In the present case the steel rings are of radius 3 mm and the compressive load is 10 MPa. Let us take a look at the development of cracking in the sample. It starts with cracks in directions 2 and 3, which are the directions impacted by the rings (Fig. 15). At some point, compressive stresses become much more important in the radial direction that in the longitudinal direction (Fig. 16). Cracking changes orientation and stresses in the longitudinal and radial directions progressively equilibrate.

Conclusion on the qualitative behavior of the model
In this section we have shown some examples of the ability of the model to predict cracking of an attacked concrete under various loadings. However, no comparison with experiments has been shown yet. This will be the focus of next section.

Comparison with experiments
In order to test our model, we need to compare its results to experimental data. We have chosen to use Mul-  ton's data [50,51]. Until more data is available, particularly on triaxial tests in which the radial and axial loads are independently imposed (such as in the testing devise developed by Giorla [31]), it seems to us it is the most complete data set for ASR under loading. Let us give some details about the test.

Multon's test
Multon performed 9 different loading cases on reactive samples, and 3 loading on non reactive sample, in order to correct the strains of the reactive sample from shrinkage and creep. Even if this correction is not perfect, we follow Multon's approach and use the corrected expansion curves. A better solution would be to simulate shrinkage and creep, which is done for example in Grimal's model [33,34] in his phenomenological model and which was developed by Giorla in a microscale model [30,31]. Since creep of cementitous materials is faster at small scales, we think that creep plays a role very locally on crack propagation. This role is neglected here, and could be investigated in a further version of the model where the cement paste would be considered as viscoelastic, as is done by Giorla through finite element simulations.
The samples are vertical cylinders of 13 cm diameter and 24 cm length. They are sealed with aluminum adhesive foils. A creep frame is used to submit them to 10 MPa and 20 MPa loads. The steel rings are either 3 mm or 5 mm thick and 1 cm high, covering all the lateral surface of the samples but not touching each other. The temperature is kept at 38˚C. The 9 tests combine 3 radial restraints: rings of 0 mm, 3 mm, and 5 mm, and three axial loadings: 0 MPa, 10 MPa, and 20 MPa. In our model we assume that the volume fraction of reactive silica in the reactive aggregates is ρ = 0.1 [49].
The particle size distribution is shown on Fig. 17. There is 23 % of non reactive aggregates (represented with a radius of 1 mm, but this radius does not matter since they don't undergo degradation). The rest of the aggregates are reactive. They are limestone crushed aggregates with silica veins, which can react. We think the created gels can easily migrate to the ITZ because of the veins. Therefore, our model seems well suited for this aggregate since most fracture phenomena will take place in the cement paste. All reactive aggregates are rather large, starting from 2.5 mm to 2.8 cm. We here underline that our model is able to deal with such variety of particle sizes, very large computations would be required to deal with a representative elementary volume of such a material.

Comments on the test results
The strain and stress computed from the measurements are given on Figs. 18 and 19. Looking at the strains (Fig. 18), we see that the free expansion is quite anisotropic, which was also mentioned by Comi and Grimal who worked on these tests results as well [34,17]. The axial strain (E multon zz ) decreases with increasing axial load, except when looking at the couples (3 mm, 10 MPa; 3 mm, 20 MPa) and (5 mm, 10 MPa; 5 mm, 20 MPa) and tends to increase with increasing rings thickness, but it is not very pronounced when looking at couples (3 mm, 10 MPa; 5 mm, 10 MPa) and (3 mm, 20 MPa; 5 mm, 20 MPa). The radial strain slightly increases with the axial load and decreases with the rings thickness. It is difficult to interpret these strains, possibly because of measurement errors, as mentioned by Multon concerning the radial strain, or of the correction for creep and shrinkage from experiments with non reactive aggregates and therefore, different stress states. The computation of the stress is more strainghforward (Fig. 19) and it seems easier to interpret: the radial stress increases with increasing axial load and rings thickness.
Let us use these experiments to identify the parameters of our model.

Parameter identification
We have performed the identification of some of the parameters of our model on Multon's data. First, we have considered that some parameters could be fixed because their values can be measured: We have fixed the aspect ratio of the cracks since the model is insensitive to this parameter as long as it is very small. We then assume that the fracture energy of the cement paste is twice that of the interface G f iss = 2G dec , which is consistent with values given in [1]. We also introduce a time scale parameter η, which links the attack depth to the physical time: We assume the attack depth is the same for all aggregates and scales with the square root of time, which was also our choice in [12], following [48,29]. The reason for that it that the progress of the attack is linked to diffusion processes. Therefore, the parameters of our model which need to be identified are K g , l c , δ, η, G dec . The optimization of the parameters is done only on the expansion curves of two experiments: (0 mm, 0 MPa) and (3 mm, 10MPa). Only the strain curves were used to identify the parameters. Then all experiments were simulated, and stress were also computed for comparison with the experimental results. The obtained set of parameters is the following: Where the time scale parameters η is written in meters per square root of days. Let us first comment on this value of η. With this value, 400 days of experiment correspond to an attack depth of 0.53 mm. Compared to the size of the aggregates (2.5 mm to 2.8 cm, see Fig. 17), we see that the smallest aggregates are significantly attacked, while for the larger ones, the attack is rather superficial. The fit of the strains on experiments (0 mm, 0 MPa) and (3 mm, 0 MPa) give satisfactory results (Fig. 18). Our model cannot reproduce the anisotropy of the free expansion. One way to do it would be to use non spherical aggregates, as shown by Giorla [31]. This improvement could be introduced in a more advanced version of the model but was not attempted yet. The other strain curves are also well reproduced for loads of 0 and 10 MPa, but the model obviously fails at 20 MPa. The decohesion does not occur at such loads, because our model currently does not allow for partial decohesion, and also forbids interpenetration of cracked surfaces. In the case of high loads (20 MPa here), the pressure is not sufficient to prevent interpenetration of the aggregate and its surrounding cement paste after decohesion. Therefore decohesion is prevented. This could be improved by modeling partial decohesion of the aggregate but would require important improvements of its micromechanical description. This discrepancy of the model is also seen on the 10 MPa tests where the brutality of decohesion induces important strain jumps. This is also due to the low number of aggregate size and the total absence of variability of the material parameters. One could point out that the expansion seems to slow down at the end of the test period (450 days), which is less visible in the predictions than in the tests results. Indeed, the depletion of reactants such as alkali is not modeled here, while it might have occurred in the tests. Some models at the same scale consider this point [48].
The results are also satisfactory in terms of stress for most cases (Fig. 19), even if stresses are usually overestimated. This can be linked to the fact that creep is not taken into account in our model, but only by correcting the expansion curves using experiments on non-reactive samples.
Finally, let us comment on the identified parameters (Eq. 40). The gel bulk modulus is roughly a tenth of that of water, and two orders of magnitude below that of C-S-H. It is also much lower that the modulus chosen by Giorla in his simulations [31] and those measured by Leeman and Lura on samples from a 45 years old dam [41]. We think that the kind of gel that forms outside the aggregates is very different from the swelling pockets found in more slowreactive concrete as modeled by Giorla. Considering the gel as a porous solid, this bulk modulus corresponds to its drained stiffness, since the experiment is long compared to the time of water transport in gels. Therefore we think the bulk modulus is in the right order of magnitude. We think that the expansion factor δ and the ITZ thickness are also realistic.
In contrast to that, the fracture energies G dec and G f iss are an order of magnitude above the values reported in the literature, for example by Wittmann or Alexander [1,67].
We think that the reason for that is that our model creates fracture patterns that are more oriented than in real concretes, where crack directions are highly influenced by the relative position of aggregates. To compensate for this intrinsic lack, the identification lead to increased values of the fracture energies, allowing cracking to occur at higher pressures and distributing cracks more isotropically.
We think that the model is well suited for loads below 10 MPa. The drawbacks of our model, such as the brutality of decohesion and the high fracture energies, do not invalidate it, but give hints to improve the model. Another drawback which was not underlined is the overestimation of the decrease of the Young's moduli with our model. We haven't shown these results, but these moduli are roughly divided by four, which is too much compared to experimental data. We think the main reason is that in our model, the decohesion removes completely the shear stiffness of the aggregates, while in real concrete this stiffness would contribute to the overall stiffness in compression. Another interesting output of the model would be to compute the resistance of the concrete at different attack stages, as done by Esposito and Hendriks [24]. This would require that the failure mechanism during a resistance test be close to that modeled here for alkali-silica reaction, which needs to be discussed.

Conclusion
In this article our main goal was to show that it is possible to build simple models for alkali-silica reaction based on micromechanics. To do so, we have chosen to represent the behavior of a concrete containing aggregates such that most of the damage occurs at the interface between the aggregates and the cement paste and in the cement paste. We have presented the micromechanical framework as well as the energy criterion that allow computing cracking around the aggregates under macroscopic loading. We have been able to identify the parameters of our model on triaxial experiments by Multon. The model reproduces the interesting phenomenon of cracking reorientation under loading and behaves well up to compressive loads of 10 MPa. The model was implemented in the Materials Ageing Platform at EDF R&D.
Some improvements are needed for the model. For example, the description of decohesion is not refined enough, inducing brutal evolutions of strain in some cases, and forbidding swelling in any direction in cases where the concrete is highly compressed.
To make the model more practical, now that we are confident that that kind of model can be built and that parameters can be identified on experiments, we plan to build a version that would describe degradation for slowreactive aggregates. That would require computing damage inside the aggregates rather than outside. On the one hand, from a micromechanical point of view, such a model seems easier to build since interaction between aggregates will be diminished if damaged is located inside the aggregates. On the other hand, on larger time-scales, the neglect of creep would become a bigger approximation. This path surely needs to be explored.