Modeling the transport of ions in unsaturated cement-based materials

Available online at www.sciencedirect.com

Computers and Structures 85 (2007) 1740–1756 www.elsevier.com/locate/compstruc

Modeling the transport of ions in unsaturated cement-based materials E. Samson *, J. Marchand CRIB – De´partement de ge´nie civil, Universite´ Laval, Ste-Foy, Que´bec, Canada G1K 7P4 SIMCO Technologies Inc., 1400, boul. du Parc Technologique, Que´bec, Canada G1P 4R7 Received 26 January 2007; accepted 20 April 2007 Available online 18 June 2007

Abstract This paper gives the details of the multiionic transport model STADIUM intended to describe the degradation of cement-based materials exposed to aggressive environments. The main algorithm is based on an operator splitting approach. The first part of the calculations is dedicated to solving the transport equations without considering the chemical reactions. The concentration profile of each ionic species is calculated by taking into account diffusion, electrical coupling between the ions, chemical activity effects and advection caused by a capillary suction flow. In the second part, a chemical module corrects the concentration profiles to enforce the equilibrium between the pore solution and the various solid phases of an hydrated cement paste. The model is compared to experimental results of hydrated cement pastes exposed to pure water and sodium sulfate solutions. Long term simulations were also performed to analyze the behavior of the algorithm.  2007 Elsevier Ltd. All rights reserved. Keywords: Multiionic diffusion; Reactive transport; Nernst–Planck; SNIA; Operator splitting; Concrete; Unsaturated; Sulfate attack

1. Introduction Concrete and other cementitious materials are very reactive media that may be significantly altered from the contact with an aggressive environment. These materials are made of three main phases in equilibrium: aqueous, solid, and gaseous. The aqueous phase, occupying a portion of the porous space, is a highly charged ionic solution containing mainly the following species: OH, Na+, K+, SO2 and Ca2+. The solid phase is a composite mixture 4 of ill-crystallized hydrated calcium silicates (C-S-H) and other crystalline phases like portlandite (Ca(OH)2). The aqueous and solid phases are in equilibrium; when the pore solution is disturbed, an amount of one or more solid phases will be either dissolved or precipitated in order to

* Corresponding author. Address: SIMCO Technologies Inc., 1400, boul. du Parc Technologique, Que´bec, Canada G1P 4R7. Tel.: +1 418 656 1003. E-mail address: [email protected] (E. Samson).

0045-7949/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2007.04.008

reach back the equilibrium state. The solid phase can also incorporate aggregates of different sizes. These aggregates are generally much less reactive than the hydrated cement paste. Consequently, their effect on the equilibrium between the paste and the aqueous solution is neglected. The gaseous phase is a mixture of dry air and water vapor. It contains no ions. Upon drying of the materials, the pressure difference between the gaseous and the aqueous phase, called the capillary pressure, will give rise to capillary suction. The resulting movement of water will result in an advection effect on the ions of the pore solution. The previous description could be applied to most porous materials. What distinguishes cement-based materials is first and foremost the high level of concentration of the ions in the pore solution, which means that the electrical coupling between the ions and chemical activity effects can hardly be neglected. Furthermore, it was already mentioned that the solid phase is very reactive; some phases are very soluble. Consequently, the sustained contact of a

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

concrete structure with water can lead to a fast deterioration of the material, caused by the leaching of the ions in the pore solution. Also, the presence of aluminum in the hydrated cement paste combined with the continuous penetration of external ions like sulfate or chloride will result in the precipitation of aluminum-based solid phases that can be very detrimental to structures during their service-life. This paper presents a multiionic transport model for saturated/unsaturated materials called STADIUM dedicated specifically to cementitious materials. Beside accounting for the classical Fick’s diffusion mechanism, it also considers the electrical coupling between the various ions as well as chemical activity effects. Furthermore, as mentioned earlier, the advection phenomenon due to capillary suction is also considered. Finally, several chemical reactions typical to cement-based materials are taken into account. The model was written under the hypothesis of a constant temperature. The solid phase is assumed non-deformable. External mechanical forces are thus not considered. Hydration of the cement is not considered in the calculations. It implies that the composition of the hydrated cement paste provided to the model as an input parameter will only change through time as a result of chemical reactions caused by the external environment. Finally, the hydrated products forming the solid skeleton are assumed to be uniformly distributed throughout the material. The model is presented for the 1D case. All examples presented in the paper are related to external sulfate attack or calcium and hydroxide leaching. 2. Ionic transport model To model the transport of ions occurring in the liquid (aqueous) phase, the equations were first written at the microscopic scale. They were then integrated over a representative elementary volume (REV, see Fig. 1) using the homogenization (averaging) technique, to yield the equations at the macroscopic scale. Details of this technique can be found in Refs. [1,2]. The macroscopic equation for the transport of ionic species i is based on the extended Nernst–Planck equation with an advection term [3]. Once integrated over the REV, it gives [4]:

Solid

Solid

Vapor Water Solid

Solid

Fig. 1. The representative elementary volume.

1741

oðws csi Þ oðwci Þ þ ot 0 ot 

1

C oB BwDi oci þ w Di zi F ci ow þ wDi ci o ln ci  ci V x C @ |{z} A ox |fflfflffl{zfflfflffl} RTffl{zfflfflfflfflfflfflfflffl ox ox |fflfflfflfflfflfflfflffl ox ffl} ffl} |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflffl diffusion

þ wri ¼ 0

electrical coupling

chemical activity

advection

ð1Þ

where ci is the concentration of the species i in solution (mmol/L), csi is the concentration in solid phase (mol/m3), ws is the volumetric solid content (m3/m3), w is the volumetric water content (m3/m3), Di is the diffusion coefficient (m2/s), zi is the valence number of the species, F is the Faraday constant (96,488.46 C/mol), R is the ideal gas constant (8.3143 J/mol/K), T is the temperature of the material (K), w is the electrical potential (V), ci (–) is the chemical activity coefficient, Vx is the average velocity of the fluid in the pore system under the action of capillary suction (m/s), and ri is a source/sink term accounting for the creation of the ion i in solution as a result of homogeneous chemical reactions (mol/m3/s). The diffusion coefficient Di is expressed as:  7=3  0 w Di ¼ sDi ð2Þ /7=3 where s (–) is the tortuosity of the liquid phase, D0i is the diffusion coefficient in freewater for species i (m2/s) and / is the porosity (m3/m3). The terms sD0i in Eq. (2) consists in the diffusion coefficient for a saturated material. When the material is unsaturated, the reduction in aqueous phase volume decreases the diffusion properties. This is taken into account by the term in parenthesis in Eq. (2). It is based on the power relationship derived by Millington and Quirk [5]. Eq. (1) is different from the classical diffusion model where only the term bearing the ‘‘diffusion’’ label in Eq. (1) is usually considered. This term models the movement of the ions as a result of their thermal agitation and is best known as Fick’s law. But ions in solution bear electrical charges. As ions have different drifting velocities, the fastest ions tend to separate from the slower ones. However, since charges of opposite signs mutually attract each other, the faster ions are slowed down and the slowest ones are accelerated, in order to bring the system near the electroneutral state. This creates an electrical potential w in the material. This term is labeled ‘‘electrical coupling’’ in Eq. (1). Measurements of diffusion potential in cementitious materials were performed by Zhang and Buenfeld [6]. To take this potential into account, Poisson’s equation is added to the model [3]. It is given here after being averaged over the REV [2]: !   N X d dw F sw zi c i ¼ 0 þ w ð3Þ dx dx  i¼1 where  = r0 is the permittivity of the medium, given by the dielectric constant times the permittivity of the vacuum (C2/N/m2) and N is the total number of ionic species. The

1742

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

permittivity of the vacuum is 8.854 · 1012 C2/N/m2. All the calculations are made with the dielectric constant of water, for which 0 = 80 at 20 C. The validity of Eq. (3) is based on the assumption that the electromagnetic signal travels much more rapidly than ions in solutions. This allows considering an equation from electrostatics in a transient problems. The dielectric permittivity of pure water is considered in the calculations. The next term in Eq. (1) that differs from a classical diffusion model is related to chemical activity. To calculate the chemical activity coefficients cis, several models are available. However, well-known models such as Debye– Hu¨ckel or Davies [7] are not adapted to the specific case of cementitious materials, which bear a highly charged pore solution. A modification of Davies’ relationship was found to yield good results [8]: pffiffi Az2i I ð0:2  4:17  105 IÞAz2i I pffiffiffiffiffiffiffiffiffiffi pffiffi þ ln ci ¼  ð4Þ 1 þ ai B I 1000 where I is the ionic strength of the solution (mmol/L): I ¼ 12

N X

z2i ci

ð5Þ

i¼1

In Eq. (4), A and B are temperature dependent parameters, given by: pffiffiffi 2 2F e 0 A¼ ð6Þ 3=2 8pðRT Þ sffiffiffiffiffiffiffiffi 2F 2 ð7Þ B¼ RT where e0 (1.602 · 1019 C) is the electrical charge of one electron. Finally, the parameter ai (m) in Eq. (4) depends on the ionic species. Its value is 3 · 1010 for OH, 3 · 1010 for Na+, 3.3 · 1010 for K+, 1 · 1010 for 10 SO2 for Cl and 1 · 1013 for Ca2+ [8]. 4 , 2 · 10 Modeling of the movement of water in the pore network under the effect of capillary suction is performed with Richard’s equation [4,9,10]. This equation is developed under the following main assumptions: isothermal conditions, isotropic material, non-deformable solid matrix, negligible gravitational effect, and water movement slow enough to have equilibrium between the liquid and the gaseous phase. The equation is given by:   ow o ow  Dw ¼0 ð8Þ ot ox ox where Dw is the moisture diffusivity coefficient (m2/s). It is the sum of the diffusivity coefficient of vapor and water. An expression is also needed for the average speed of the liquid phase, which appears in Eq. (1). It is given [9] by: V x ¼ DL

ow ox

ð9Þ

where DL is the water diffusivity coefficient (m2/s). Substituting Eq. (9) in (1) gives:

oðws csi Þ oðwci Þ þ ot  ot  o oci Di zi F ow o ln ci ow  wDi ci þ wDi ci þw þ DL ci ox RT ox ox ox ox þ wri ¼ 0

ð10Þ

There are two terms remaining in the model that need to be discussed. The first one on the left-hand side of Eq. (10) involves csi , the concentration of species i in solid phase. This term is used to describe the exchange between the solid and the aqueous phase [1]. The other term is ri, which is a source/sink term in the aqueous to model the creation or removal of ion i. They will be considered in the next section. 3. Modeling the chemical reactions There are several types of chemical reactions (see Ref. [11] for a comprehensive review). The homogeneous reactions are those occurring solely in the aqueous phase, like the formation or the dissociation of acids. The heterogeneous reactions involve more than one phase, like the solid and aqueous ones, for example. They can be divided in two groups: surface and classical reactions. The surface reactions occur at the solid/liquid interface. Adsorption of ions, i.e. the capture of ions by the surface of the solid as a result of electrostatic forces, falls in this category. Finally the dissolution/precipitation reactions are part of the classical category. Dissolution and/or precipitation occur when a solid is no longer in equilibrium with the solution with which it is in contact. If the activity product of the ions in solution involved in the reaction is above the equilibrium constant, some of the ions will precipitate in order to reach equilibrium. In the inverse situation, dissolution will occur. In Eq. (10), the term oðws csi Þ=ot corresponds to heterogeneous chemical reactions whereas the term ri corresponds to homogeneous ones. 3.1. Operator splitting approaches Different techniques have been used to include chemical reactions in ionic transport problems. In early models (see for instance Miller and Benson [12]) the equations that model the chemical reactions were solved simultaneously with the transport relationships. The current trend in reactive transport modeling is to split the transport part of the process from the chemical reactions. This method gained a lot of popularity after the publication of a paper by Yeh and Tripathi in 1989 [13] which showed that important CPU times could be saved by using such an uncoupled approach. Models by Grove and Wood [14], Yeh and Tripathi [15], Walter et al. [16] and Xu et al. [17] are all based on this operator splitting approach. The operator splitting approach is divided in two main classes. First, it is possible to iterate between the transport and chemical reaction set of equations within a time step until convergence is reached. This is called the sequential

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

iterative approach (SIA). Another algorithm consists in solving the transport and chemical reaction part of the model without iterations between these two modules. It is called the sequential non-iterative approach (SNIA). The papers cited previously [14–17] are all based on the SNIA. The operator splitting procedure introduces a numerical error in the resolution of a reactive transport problem. Systematic analyses were performed on simple models to evaluate this error. Herzer and Kinzelbach [18] studied the equation:   ws oc o2 c oc D 2þv ¼0 ð11Þ kd þ 1 ot ox ox / where ws is the solid volumetric content, / is the porosity, c is the concentration, D is the diffusion coefficient, v is the fluid velocity and kd relates c and the concentration cs in solid phase: cs = kdc. The difference between the exact solution and the solution calculated by separating transport and chemistry is given by: err ¼

ðR  1Þ vDt R 2

ð12Þ

where R = 1 + kdws// and Dt is the time step. The error is thus proportional to Dt and tends to zero as the time step is decreased. Valocchi and Malmstead [19] studied the equation: oc o2 c oc  D 2 þ v þ kc ¼ 0 ot ox ox

ð13Þ

where k characterizes a radioactive decay process. The numerical error introduced in the solution due to operator splitting is given by:   1 err ¼ c0 v ð1  ekDt Þ  Dt ekDt ð14Þ k where c0 is related to the flux boundary condition as: (vc  Dc,x)jx=0 = vc0. Again, the error goes to zero as the time step is reduced. Kaluarachchi and Morshed [20] studied the same case as Valocchi and Malmstead but with the following boundary condition:   D oc  þc ¼ ekt ð15Þ vL ox x¼0 where k = 0 and k ! 1 correspond to continuous and pulse boundary conditions, respectively. The conclusions

1743

are the same, the error being reduced as Dt is reduced. Barry et al. [21] studied the more general case: oc os o2 c oc þ D 2þv ¼0 ot ot ox ox

with

os ¼ kc  rs þ g ot

ð16Þ

with both homogeneous and inhomogeneous boundary conditions. The results are summarized in Table 1. The results are in agreement with those presented previously. Except for the unusual homogeneous boundary condition case, the error introduced by the operator splitting procedure is of the order of the time step value. Reducing the time step leads to an increase in calculation time. On the opposite, the SIA does not require time steps as small as in the SNIA, because the iterations between transport and chemistry reduce the splitting error [18]. But iterations within a time step increase CPU time. Walter et al. [16] studied both approaches to determine which is the most efficient. They modeled a test-case that illustrates the controlling geochemical processes occurring in a carbonate aquifer under the impact of acidic mine tailings effluent. Simulations were performed with both SIA and SNIA algorithms. They involved complexation, dissolution/precipitation, redox and gas exchange chemical reactions. The various profiles calculated with both algorithms showed less than a 2% difference in position, calculated with regard to the position of the source of contaminant. However, the CPU load was 3.5 times more important for the SIA model. Following these results, we elected to work with the SNIA algorithm. 3.2. Coupling multiionic transport with chemistry The transport model considered in Refs. [14–17] is a simplified version of Eq. (10). It assumes that the fluid flow is strong enough to neglect diffusion to the profit of dispersion effects (see Ref. [1]). Other phenomena like electrical coupling between the ions and chemical activity (see Eq. (10)) are also neglected, which leaves the linear transport operator: Lðci Þ ¼ D gradðci Þ þ ci v

ð17Þ

where D is the dispersion tensor and v is the fluid velocity, both independent on the ionic species i. This allows to rewrite the transport equations in order to eliminate the terms ri (this topic is discussed in paper [22]). Accordingly the treatment of homogeneous reactions is greatly simpli-

Table 1 Analysis of the accuracy of the operator splitting scheme on various simple models (from Ref. [21]) Case

Model

Accuracy

A B

s_ ¼ g; k ¼ 0; r ¼ 0 s_ ¼ kc; r ¼ 0; g ¼ 0

C

s_ ¼ r½ðR  1Þc  s; g ¼ 0; k = r(R  1) s = (R  1)c, r ! 1 in case C

O(Dt2) O(Dt) O(Dt2) O(Dt)

D

O(Dt) O(Dt2)

Notes

Homogeneous boundary conditions

Homogeneous boundary conditions

1744

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

fied. These simplifications are common to almost all ionic transport models applied to geochemical processes. Such simplifications are not possible for a nonlinear transport model like Eq. (10), which complicates the chemical reaction modeling. However, given the highly reactive nature of the hydrated cement paste, it can be assumed that heterogeneous reactions are more important than homogeneous one in cementitious materials. For example the penetration of sulfate ions may lead to the formation of ettringite and gypsum in concrete in such a quantity that the material sustains physical damage. Portlandite, which occupies more than 5% of the total volume of the material, is highly soluble; its rapid dissolution can have important effect on the overall strength of a concrete structure. Accordingly, the homogeneous reactions are neglected in the present model, which allows to eliminate the term ri from Eq. (10). Only dissolution/precipitation reactions will be considered. Within the framework of the SNIA algorithm, the term involving csi is removed from Eq. (10), which leaves the following equation to solve during the transport step:   oðwci Þ o oci Di zi F ow olnci ow þ D L ci  wDi þ w ci þ wDi ci ¼0 ox ox ot ox RT ox ox

ð18Þ

This equation still has to be coupled to Eqs. (3) and (8). The transport part of the process is thus closed. The unknowns are w, w, and N times ci, i.e. a concentration for each ionic species taken into account. There is accordingly N equations (18), one for each ionic species, Eq. (3) to solve the electrical potential and Eq. (8) for the water content. The other parameters are either physical constants, like T, F, R, zi, e0, , or material parameters that have to be determined experimentally: Di, s, Dw, DL. At the end of this calculation step, ionic profiles are known for each ionic species. After the transport step, the concentrations may not be in equilibrium with the solids. They will be corrected in the chemical reaction step, by adjusting the solid phases properly. This is achieved by enforcing the algebraic relationship describing the equilibrium state of each solid phase present at a given location. This relationship is written in a general form as [17]:

Km ¼

N Y

cmi mi cmj mi

with m ¼ 1; . . . ; M

ð19Þ

i¼1

where M is the number of solid phases, N is the number of ions, Km is the equilibrium constant (or solubility constant) of the solid m, ci is the concentration of the ionic species i, ci is its chemical activity coefficient, and mmi is the stoichiometric coefficient of the ith ionic species in the mth mineral. The right-hand side of Eq. (19) is called the ion activity product. To illustrate the chemical modeling, we consider a problem where portlandite (CH), ettringite (AFt) and gypsum (G) are involved (see Table 2). The chemical equilibrium equation for each solid is: K CH ¼ cCa c2OH ½Ca½OH K AFt ¼

2

ð20Þ

6 4 3 2 c6Ca c4OH c3SO4 c2AlðOHÞ4 ½Ca ½OH ½SO4  ½AlðOHÞ4 

K G ¼ cCa cSO4 ½Ca½SO4 

ð21Þ ð22Þ

where the square brackets represent the concentrations. If the concentrations in solution do not respect the equilibrium equations (20)–(22), solids will either be dissolved or precipitated until these algebraic relationships are satisfied. Let us assume that the concentrations c0i (e.g. [Ca0], [OH0], 0 ½SO04 , ½AlðOHÞ4 ) are the output of the previous transport calculation step and do not respect equilibrium. The equilibrium concentrations can be expressed as: ci ¼ c0i þ

M X

ð23Þ

mmi X m

m¼1

where the Xms represent the amount of a given solid that has to dissolve in order to reach the equilibrium state. This is similar to the approach proposed in Ref. [23]. Following this, the system of Eqs. (20)–(22) is rewritten as: K CH ¼ cCa c2OH ½Ca0 þ X CH þ 6X AFt þ X G ½OH0 þ 2X CH þ 4X AFt 2

ð24Þ K AFt ¼

c6Ca c4OH c3SO4 c2AlðOHÞ4 ½Ca0

þ X CH þ 6X AFt þ X G  4

0

6 3

 ½OH þ X CH þ 4X AFt  ½SO04 þ X G þ 3X AFt  0

2

ð25Þ

 ½AlðOHÞ4 þ 2X AFt  0

K G ¼ cCa cSO4 ½Ca þ X CH þ 6X AFt þ

X G ½SO04

þ X G þ 3X AFt 

ð26Þ

Table 2 Solid phases in hydrated cement Solid phase

Chemical description

Equilibrium relationship

log(Ksp)

Portlandite C-S-H Ettringite Monosulfates Gypsum Mirabilite

Ca(OH)2 1.65CaO Æ SiO2 Æ (2.45)H2Oa 3CaO Æ Al2O3 Æ 3CaSO4 Æ 32H2O 3CaO Æ Al2O3 Æ CaSO4 Æ 12H2O CaSO4 Æ 2H2O Na2SO4 Æ 10H2O

Ksp = {Ca}{OH}2 Ksp = {Ca}{OH}2 b Ksp = {Ca}6{OH}4{SO4}3{Al(OH)4}2 Ksp = {Ca}4{OH}4{SO4}{Al(OH)4}2 Ksp = {Ca}{SO4} Ksp = {Na}2{SO4}

5.2 6.2 44.0 29.1 4.6 1.2

{  } indicates chemical activity. a A C/S of 1.65 is assumed for the C-S-H. b The C-S-H decalcification is modeled as portlandite dissolution with a lower Ksp.

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

This nonlinear system of equations is solved for the Xms. Following the convention adopted previously, a positive Xm means dissolution and a negative one indicates the precipitation of solid m. At the end of the chemical equilibrium calculation, each solid phase S is adjusted according to the following relationship: S tm ¼ S mt1  wX m Cm =q with m ¼ 1; . . . ; M

ð27Þ

where Sm is the amount of a given solid phase (g/kg of material), t indicates the time step, Cm is the molar mass of the solid m (g/mol) and q is the density of the saturated, wet material (kg/m3). The value of the density q is about 2400 kg/m3 for concrete and 1700 kg/m3 for neat cement paste. The ionic concentrations are also adjusted following the chemical equilibrium calculation, according to Eq. (23). The solids considered in this paper are listed in Table 2, along with their equilibrium constant. They are portlandite, C-S-H, ettringite, monosulfates, gypsum and mirabilite. The C-S-H is treated in a particular way. This solid phase exhibits an incongruent behavior [24], which means that the dissolution rate of the ions composing this mineral does not follow the stoichiometric coefficient. As a consequence, the composition of the C-S-H varies as the dissolution progresses. According to Ref. [24], the calcium to silicate ratio of the C-S-H (C/S) starts at about 1.65 in the sound hydrated cement paste and decreases as calcium is lost. Berner modeled the incongruent behavior of the CS-H by assuming that it is composed of two different solids, Ca(OH)2 and CaH2SiO4, each having its own equilibrium constant that depends on the C/S ratio. The part associated to portlandite dissolves first. The remaining CaH2SiO4 dissolves very slowly [24]. Studies have shown that C-S-H exposed to pure water for a long time keeps a C/S ratio between 0.8 and one in the decalcified zone [25]. Considering this, we model the dissolution of C-S-H by considering only the Ca(OH)2 part. Berner’s approach is further simplified by assuming that this fraction of the C-S-H dissolves with a constant equilibrium value. Its value is an average of the solubility constant of Ca(OH)2 at C/S = 1.65 and C/S = 1.0. It is given in Table 2. The chemical reactions, beside capturing or releasing ions as solid phases are precipitated or dissolved, will have an effect on the transport properties by affecting the porosity of the material. If for instance gypsum is formed at some location, the local porosity will decrease, thus reducing the area across which ions are able to diffuse. This will reflect on their diffusion coefficient. Models for ionic transport in groundwater often use a modified version of the Kozeny–Carman relationship [26] to take this effect into account:  H D ð/Þ ¼

/ /0

3 

1  /0 1/

2 ð28Þ

where /0 is the initial porosity of the cementitious material (m3/m3). For the present model, a relationship was devised

1745

based on migration and porosity test results obtained on a wide range of material [27]. The relationship is given as: H D ð/Þ ¼

e4:3/=V p e4:3/0 =V p

ð29Þ

where Vp is the paste volume (m3/m3) of the material. Upon chemical modification to the material, the porosity of the paste is calculated as: / ¼ /0 þ

M X ðV init m  V mÞ

ð30Þ

m¼1

where Vm is the volume of a given solid phase, per unit volume of material, and M is the total number of solid phases. The function HD multiplies the diffusion coefficient of each ionic species. It is also assumed that water diffusivity is affected similarly by chemical reactions. 3.3. Validity of the local equilibrium assumption The use of equilibrium relationships like Eq. (19) to model chemical reactions is valid under the local equilibrium assumption (LEA). This hypothesis implies that chemical reaction time scales are shorter than the time scales of the transport processes. On the contrary, if transport processes are fast compared to chemical reactions, the equilibrium relationships are no longer valid and must be replaced by kinetic expressions. The validity of the LEA is checked by performing an analysis based on the dimensionless Damko¨hler number, given by [28]: tD DaD ¼ ð31Þ tr where tD and tr are the characteristic times for diffusion and chemical reactions respectively. The parameter tD is calculated as: tD ¼

l2D Deff

ð32Þ

where lD is a characteristic diffusion length and Deff is the effective diffusion coefficient. The characteristic time for the chemical reactions is calculated according to: tr ¼

1 k eff

ð33Þ

where keff is the effective rate constant. If Da  1, the time scales of the reactions are much shorter than the time scales associated with ionic diffusion. In that case, the local equilibrium assumption is valid [28]. Assuming a value of 5 · 1011 m2/s for Deff (see Ref. [27]) and 2 mm for lD (which corresponds to a typical element size), one gets tD = 8 · 104 s. The calculation of tr is based on data found in Ref. [29]. In this paper, the highly reactive nature of hydrated cement paste is taken into account through reaction rate constants k that are approximately six orders of magnitude higher than for rock minerals. In the case of ettringite and port-

1746

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

landite for instance, k = 1 · 108 mol/m2/s. Values of surface area for hydrated cement pastes are given in Ref. [30]. A typical value of 100 m2/gdry paste is used for the present calculations. The amount of dry paste in the material is calculated using the following mixture proportions of a 0.45 w/c concrete: 380 kg/m3 of cement (C), 170 kg/m3 of water (W), 2400 kg/m3 density, and 0.12 m3/m3 porosity. The amount of dry paste Mp in the concrete is given by: M p ¼ C þ W  1000/ ¼ 380 þ 170  1000  0:12 ¼ 430 kgdry paste Using the methodology presented in Section 5, the amount of ettringite in the paste is estimated at 20 g/kgconcrete. The reactive area Ar is thus given by: Ar ¼ 100

m2 gdry paste



430 gdry paste 1000 gconcrete 1255 gettringite   2400 gconcrete 20 gettringite mol

 1  106 m2 =mol

4. Numerical model

Using this result, the effective rate constant is given by keff = kAr = 1 · 108 · 1 · 106 = 0.01 s1. According to Eq. (33), the characteristic time for the chemical reactions is tr = 100 s. With these values for tD and tr, the Damko¨hler number is 800. It thus meets the validity criteria for the LEA. Similar calculations have been performed with lower surface area values and higher diffusion coefficients. In all cases, the criteria Da  1 was respected. The low diffusion coefficients, high reactive rate of the paste and high surface area typical of cementitious materials thus contribute to validate the use of the chemical equilibrium algebraic expressions to model reactions coupled with ionic transport. The same analysis can be performed based on a characteristic advection time ta instead of tD. In that case: ta ¼

la v

ð34Þ

where v (m3/m2/s) is the moisture flow in the material and la is a characteristic advection length. The corresponding Damko¨hler number is given by: Daa ¼

ta tr

2 mm for la, the characteristic advection time ta is 4.3 · 106 s for the 0.65 mixture. With the value of tr = 100 s calculated in the previous paragraphs, the advection Damko¨hler number Daa is 4 · 104. This again confirms the validity of the local equilibrium assumption for ionic transport in cementitious materials. The previous results are in agreement with analyses published previously. Barbarulo et al. [31] analyzed the validity of the LEA based on a dimensional analysis. They concluded that the assumption is valid in most practical cases of ionic diffusion in cementitious materials. Knapp [32] published a study where the limit of validity of both kinetic and equilibrium approaches is evaluated. The analysis is based on two dimensionless parameters: the Peclet number Pe = vl/D and the Damko¨hler number Da. Using the parameters of the previous paragraphs with Knapp methodology also confirms the validity of the local equilibrium assumption.

ð35Þ

Typical values of v were obtained by performing ISO 12572:2001(E) vapor transmission tests on ordinary concrete mixtures prepared with an ASTM Type V cement. After 28 days of curing in a fog room, samples were exposed to a 100–50% relative humidity gradient until steady-state is reached. The moisture flow for a 0.45 water to cement ratio concrete was 0.023 kg/m2/d (2.7 · 1010 m3/ m2/s). A flow of 0.040 kg/m2/d (4.6 · 1010 m3/m2/s) was measured on a 0.65 concrete. Using again a value of

This section is divided in two parts. First, the numerical algorithm for the transport model, i.e. Eqs. (18), (3) and (8), is presented. Then, the details on the resolution of the equilibrium problem in the chemical part of the model are given.

4.1. Finite element discretization of the transport equations The numerical algorithms for the Nernst–Planck/Poisson model, i.e. the diffusion of ions taking into account the electrical coupling but not the chemical activity nor the capillary suction, has already been studied [34]. The authors tested two algorithms using the finite element method. The first one, with the equations uncoupled and solved by successive iterations, was found to converge only for very small concentrations. This algorithm is not suitable for cement-based materials since their porous network is filled with an highly charged aqueous solution. In the other algorithm, equations were solved simultaneously with the Newton–Raphson method. The results in this case were not limited to a given concentration level. The same algorithm was later used to solve the extended Nernst– Planck model [35], which takes into account chemical activity gradients. Following these results, the coupled algorithm will be used with the current model, which adds a convection term to the extended Nernst–Planck equation. Numerical tests will be performed to assess the validity of this approach. To ease the writing, the numerical model will be shown for a two ion case in 1D. But the model can easily be extended to take into account more species. The numerical examples that will be presented later involve six ions. The weighted residual form of the set of Eqs. (18), (3) and (8), upon integration by parts, is given as:

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

38 9 w 0 c1 0 > > > c_1 > > > > 7> 6 Z = < c_2 > 6 0 w c2 0 7 > 7 6 dx W ¼ hdWi6 7 6 0 0 1 0 7> L w_ > > > > 5> 4 > > > ; : _ > 0 0 0 0 w 9 38 2 1 z1 F c1;x > wD1 0 DL c1 DRT wc1 > > > > > > 7> > > 6 Z  D2 z2 F = < 7 6 c 2;x 0 wD D c wc odW 6 2 L 2 27 RT dx þ 7 6 > 7> ox 6 0 w L > > ;x 0 D 0 w > > 5> 4 > > > ; : w;x 0 0 0 ws 38 9 2 c1 > wD1 oðlnoxc1 Þ 0 0 0 > > > > > > > 7 6 > >   Z oðln c2 Þ < 7 6 0 wD2 ox 0 0 7 c2 = odW 6 þ dx 7 6 >w> ox 6 L > 0 0 0 07 > > 5> 4 > > > ; : > w 0 0 0 0 38 9 2 0 0 0 0 > > > c1 > > > > 7> 6 > > Z = < 7 6 0 0 0 0 c 2 7 6 þ hdWi6 dx þ W ext ¼ 0 7 > 7> 6 0 0 0 0 L w > > > 5> 4 > > > ; : > F F   z1 w   z2 w 0 0 w 2

ð36Þ ext

where W includes the boundary integrals associated to the flux boundary conditions and hdWi is the vector of the weighting function, defined as: hdWi  hdc1 dc2 dw dwi

ð37Þ

The weak form is discretized using the Galerkin method with a standard linear two-node element (see Ref. [36] for a complete text on the method). The approximation of the solution on each element is written as: 8 9 c1 > > > > > > > = < c2 > ¼ ½N fU n g ð38Þ > > w > > > > > > : ; w 2 3 0 0 N2 0 0 0 N1 0 6 7 0 0 N2 0 0 7 6 0 N1 0 7 ð39Þ ½N  ¼ 6 6 0 0 0 N2 0 7 0 N1 0 4 5 0

0

0

N1

0

0

hU n i ¼ hc11 c21 w1 w1 c12 c22 w2 w2 i

0

N2 ð40Þ

where N1 and N2 are the shape functions. The subscripts i and j for the concentrations in the vector hUni (Eq. (40)) designate the species i at the node j of one element. The elementary matrices are thus expressed as: Z ½K e  ¼ ð½BT ½D1 ½B þ ½BT ½D2 ½N  þ ½N T ½D3 ½N Þdx ð41Þ e Zl ð42Þ ½M e  ¼ ½N T ½C½N dx le

1747

with 2

N 1;x 6 0 6 ½B ¼ 6 4 0 0 2

0 N 1;x 0 0 0

wD1

0 0 N 1;x 0

0 0 0 N 1;x

D L c1

N 2;x 0 0 0

D 1 z1 F RT D 2 z2 F RT

6 wD2 DL c2 6 0 ½D1  ¼ 6 4 0 0 Dw 0 0 0 2 oðln c1 Þ wD1 ox 0 6 6 0 wD2 oðlnoxc2 Þ ½D2  ¼ 6 6 0 0 4 0 0 2 0 0 0 6 0 0 0 6 ½D3  ¼ 6 4 0 0 0 F F   z1 w   z2 w 0 2 3 w 0 c1 0 6 0 w c 07 2 6 7 ½C ¼ 6 7 4 0 0 1 05 0 0 0 0

wc1

0 N 2;x 0 0 3

7 wc2 7 7 5 0 ws 3 0 0 7 0 07 7 7 0 05 0 0 3 0 07 7 7 05 0

0 0 N 2;x 0

3 0 0 7 7 7 0 5 N 2;x ð43Þ

ð44Þ

ð45Þ

ð46Þ

ð47Þ

The assembly of elementary matrices [Ke] and [Me] leads to the following system of equations: ½MfU_ g þ ½KfU g ¼ fF ext g

ð48Þ

The time discretization is performed using the implicit Euler scheme:  U t  U tDt ð49Þ ½M t  þ ½K t fU t g ¼ fF ext t g Dt where Dt is the time step. The subscript t stands for the actual time step and t  Dt the previous one. Defining the matrices ½K ¼ ½M t  þ Dt½K t 

ð50Þ

fF ext t g

ð51Þ

fF g ¼

þ ½M t fU tDt g

the system of equations can be written as: ½KfU t g ¼ fF g

ð52Þ

This nonlinear system of equations is solved at each time step with the modified Newton–Raphson method. Numerical simulations have shown that the convergence rate with a tangent matrix calculated without the nonlinear terms arising from the coupling between the various variables in the model is almost the same as the one with a complete tangent matrix. However, the calculation time is reduced since fewer terms need to be calculated. Furthermore, its wider radius of convergence makes this algorithm even more interesting. The elementary tangent matrix is thus given as:

1748

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

½K eT  ¼ ½M e  þ Dt½K e 

ð53Þ

The elementary matrices in Eqs. (41) and (42) are evaluated through Gauss numerical integration method. Accordingly, the variables appearing in the matrices [D1], [D2], [D3], and [C] are calculated at the integration points. To calculate the terms o(ln ci)/ox in (45), the ionic strength (Eq. (5)) is calculated at each nodes. Then the (ln ci)s are calculated for each species and at each nodes with Eq. (4). The value of o(ln ci)/ox can then be evaluated at every integration points. 4.2. Solving the chemical equilibrium equations After solving the transport equation, the concentrations are used as input in the chemical equilibrium module. The numerical method is illustrated using the nonlinear system of Eqs. (24)–(26). Newton’s algorithm is used to solve the system of equations: k k1 ½K k1 Þg T fDX g ¼ fRðX k

X ¼X

k1

þ DX

ð54Þ

k

ð55Þ

All terms in the Jacobian matrix are calculated with values found at the previous iteration level. The system of equation is solved by iterating over Eqs. (54) and (55). This calculation procedure is repeated two times. The first time, the chemical equilibrium calculation is made only for solids that are present at a given node. During the second step, the calculation is made for solid present at the node and also for the solids for which the ion activity product is greater than the equilibrium constant. This procedure allows avoiding the impossible situation where dissolution of a non-existing solid occurs at a given location. 5. Experimental validation Cement pastes were prepared at a water/cement ratio of 0.6 using an ordinary Portland cement (Canadian CSA Type 10). The chemical composition of the cement and the mixture proportions are listed in Table 3. The mixture was prepared using deionized water and without any chemical admixture. It was batched under vacuum (at 10 mbar)

where [KT] is the tangent (or Jacobian) matrix, DX is a variation calculated at each iteration, {R} is the residual vector, and k indicates the iteration level. The residuals are given by:

Table 3 Parameters for the numerical simulations on the cement paste submitted to the degradation tests Properties

Values

RCH ¼ cCa c2OH ½Ca0 þ X CH þ 6X AFt þ X G ½OH0 þ 2X CH þ 4X AFt 2  K CH

Cement type w/c

CSA, T10 0.6

Mixture proportions (kg/m3) Cement Water Coarse aggregates Fine aggregates Density

1088.8 653.3 – – 1742.1

Cement composition (% mass) CaO SiO2 Al2O3 SO3

62.04 19.78 4.39 3.2

Initial solid phases (g/kg) Portlandite C-S-H Ettringite Monosulfate

198.8 375.2 34.1 104.8

Porosity (%)

52.20

Diffusion coefficients (1011 m2/s) OH Na+ K+ SO2 4 Ca2+ AlðOHÞ 4

16.1 4.1 6.0 3.3 2.4 1.7

Initial pore solution (mmol/L) OH Na+ K+ SO2 4 Ca2+ AlðOHÞ 4

429.3 111.1 327.0 5.6 1.3 0.2

ð56Þ RAFt ¼ c6Ca c4OH c3SO4 c2AlðOHÞ4 ½Ca0 þ X CH þ 6X AFt þ X G 6 0

 ½OH þ 2X CH þ 4X AFt  

½AlðOHÞ04

4

½SO04

þ X G þ 3X AFt 

3

2

ð57Þ

þ 2X AFt   K AFt

0

RG ¼ cCa cSO4 ½Ca þ X CH þ 6X AFt þ

X G ½SO04

þ X G þ 3X AFt   K G

ð58Þ

The Jacobian matrix is given by: 2 oR 3 oR oR CH

CH

CH

6 CH 6 AFt K T ¼ 6 oR 4 oX CH

oX AFt

oX G

oRAFt oX AFt

oRAFt oX G

oRG oX AFt

oRG oX G

oX

oRG oX CH

7 7 7 5

ð59Þ

The calculation of each term in this matrix is made by considering the activity coefficient constant. The latter are updated after each iteration from the updated concentrations oRCH (see Eq. (23)). For example the term oX in the Jacobian CH matrix is given by:

 oRCH oRCH o½Ca oRCH o½OH 2 ¼ cCa cOH þ ð60Þ oX CH o½Ca oX CH o½OH oX CH where the concentrations [Ca] and [OH] are calculated according to Eq. (23). The calculation of Eq. (60) thus gives: oRCH 2 ¼ cCa c2OH ½½OH þ 4½Ca½OH oX CH

ð61Þ

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

tiplied by the appropriate factor, i.e. multiply each quantity by 0.9 for a 90% hydration degree. Given the high water to cement ratio of the paste, complete hydration is assumed. The remaining disks were submitted to a degradation test. It consists in exposing the pastes to different controlled environments. The 20 mm disks were first sealed on the side and on one face with silicon. They were then immersed in 30 L of aggressive solution for a minimum of 3 months. Two different solutions were used: deionized water and a 50 mmol/L Na2SO4 solution. In the case of the Na2SO4 solution, disks were also exposed during 6 months and 1 year. During the exposure period, the solution was renewed on a weekly basis in order to maintain uniform conditions. After the test, small prisms (20 mm · 3 mm · 12 mm) were cut from the disk and immersed in isopropylic alcohol for 2 weeks. They were then dried under vacuum for a week. The samples were finally impregnated with an epoxy resin, polished and coated with carbon. After the conditioning was completed, the samples were analyzed with a microprobe (Cameca SX-100 operating at 15 kV and 20 nA), which allows to determine their degradation state. The analysis was performed over 1000 points along the length of the sample. Profiles of calcium and sulfur were measured and compared to the numerical model. The microprobe profiles are expressed in count/s. The simulations were performed with the parameters listed in Table 3. The boundary conditions are of the Dirichlet type. For the deionized water case, all concentrations are set to zero at x = 0. The water content is set to the porosity value of 0.522 at x = 0. Finally, the potential w is also set to zero at x = 0. No boundary conditions are applied at x = L. The spatial discretization was performed with a uniform 100 element mesh over a 20 mm domain. Time steps of 1800 s were used. The output of the model consists in the profiles of the variables, i.e. the N concentrations along with the water content and the electrical potential, as well as the profiles of the M solid phases. The total calcium content is

CCa CCa CCa CCa CCa S CH þ 1:65 S csh þ 6 S AFt þ 4 S AFm ¼ 10CaO CCH Ccsh CAFt CAFm CCaO CSi CCa S csh ¼ 10SiO2 ð63Þ Ccsh CSiO2 CAl CAl CAl 2CAl 2 S AFt þ 2 S AFm þ 2nAl2 O3 S csh ¼ 10Al2 O3 CAFt CAFm CAl2 O3 CAl2 O3

ð64Þ 3

CS CS CS S AFt þ S AFm ¼ 10SO3 CAFt CAFm CSO3

ð65Þ

Microprobe (count/sec)

ð62Þ

1600

400

1400

350

1200

300

1000

250

800

200

600

150

400

100 Microprobe measurement

200

where the Cis are the different molar masses. The term nAl2 O3 is an amount of Al2O3 that substitutes in the C-SH. According to Ref. [38], that substitution rate varies between 2% and 5%. The value of 2% is used for the calculations. To account for the hydration degree of the material, the quantity of CaO, SiO2, Al2O3 and SO3 should be mul-

Total calcium content (g/kg)

to prevent the formation of air void during mixing. It was then cast in plastic cylinders (diameter = 7 cm; height = 20 cm). Molds were sealed and rotated for the first 24 h to prevent bleeding. Cylinders were then demolded and sealed with an adhesive aluminum foil for at least 12 months at room temperature. After the hydration period, the cylinders were sawn in disks for testing, in order to determine the different parameters needed for the ionic transport simulations. Prior to testing, the disks were vacuum saturated in a 300 mmol/L NaOH solution during 24 h. The porosity of the paste was determined according to the ASTM C642 standard procedure. A value of 52.2% was measured. Since the ionic transport during the validation tests occurs in saturated condition, the water content is assumed constant in the material and has the same value as the porosity. The pore solution was extracted following the experimental procedure described in Ref. [37]. It allows determining the concentration of these ions: OH, Na+, K+, SO2 4 , and Ca2+. Because of experimental error, the solution did not respect electroneutrality. The concentration in OH was adjusted to make sure that the positive and negative charges are equal. The concentration of AlðOHÞ 4 is evaluated with the chemical equilibrium code described in the previous section in order to have a solution in equilibrium with portlandite, ettringite and monosulfates. The ionic diffusion coefficients were determined from migration tests. A complete description of the testing method is given in Ref. [33]. It consists in accelerating the transport of ions in a saturated concrete sample by applying an external electrical potential. The procedure is similar to the ASTM C1202 test method. The test was performed on two 20 mm disks. The measured diffusion coefficients parameters are listed in Table 3. The table also includes the initial amount of each solid phase in the material: portlandite (CH), C-S-H, ettringite (AFt) and monosulfates (AFm). It is calculated from the cement composition (CaO, SiO2, Al2O3 and SO3) by solving the following set of equations:

1749

50

Simulation

0

0 0

2

4 6 Position (mm)

8

10

Fig. 2. Validation test – comparison of the model with the calcium profile for the paste exposed to water for 3 months. The boundary x = 0 corresponds to the exposed surface.

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

obtained from the addition of the calcium in each solid phase. It is then converted in grams of calcium per kg of dry material. The same procedure is followed for sulfur. The comparison between the numerical simulations and the microprobe measurements is shown in Figs. 2–4. Fig. 2 shows the comparison of the model with the experiment for the case of the cement paste exposed 3 months to water. The experimental profile show high peaks of calcium starting at about 2 mm and going deeper toward x = L. They correspond to portlandite. At around 2 mm, these peaks disappear, which indicates the dissolution of portlandite. A plateau is observed between 0.8 and 2 mm, followed by a drop of calcium near the surface of the material. This drop in calcium content is associated with the decalcification of C-S-H. The plateau indicates that the decalcification of C-S-H does not start immediately after the dissolution of portlandite. The model was able to reproduce with accuracy the main features of the microprobe profile. As shown in Fig. 2, the predicted dissolution of portlandite begins at around 2 mm. The simulations show the presence of a plateau and matches very well with the experiment for the calcium drop near the surface. The model is then compared to the sulfur profiles measured on the paste exposed to sodium sulfate at three differ-

3500

ent time intervals as shown in Fig. 3. Following the x-axis from 10 mm to the surface, the measurements show a low level of sulfate, corresponding to the level in the sound paste. At some point, the level of sulfate increases. This increase is located roughly at 4.5 mm after 3 months, 5.5 mm after 6 months and 8 mm after a year. It corresponds to the penetration of an ettringite front in the material, as a result of the exposure to sodium sulfate. Finally, an important peak of sulfur is observed near the surface, corresponding to gypsum formation. Fig. 3 shows that the model reproduces very well the experimental profiles. The penetration of the ettringite front is accurately followed, as well as the position of the gypsum peak. Furthermore, the model also predicts the widening of the gypsum peak through time. The same analysis is done for the calcium profiles at 3, 6 and 12 months for the paste exposed to sodium sulfate. The results are shown in Fig. 4. In all cases, the model predicts accurately the dissolution of portlandite. In Fig. 4b and c, it was observed that the plateau of calcium between the dissolution of portlandite and the decalcification of C-S-H is higher than in Fig. 2. This can be attributed to the gypsum peak, which forms right next to the portlandite dissolution front. The model predicts this higher calcium plateau. The

120

3500

120

Microprobe measurement

Microprobe measurement

2000 60 1500 40 1000 20

500 0

0 2

4 6 Position (mm)

8

Microprobe (count/sec)

80

Totalsulfur content (g/kg)

2500

100

Simulation

2500 80 2000 60 1500 40 1000 20

500 0

10

0 0

2

(a) Three month exposure

4 6 Position (mm)

8

(b) Six month exposure

3500

120 Microprobe measurement

3000

100 Simulation

2500 80 2000 60 1500 40 1000

Total sulfur content (g/kg)

0

3000

100

Simulation

Microprobe (count/sec)

Microprobe (count/sec)

3000

20

500 0

0 0

2

4 6 Position (mm)

8

10

(c) Twelve month exposure

Fig. 3. Validation test – comparison of the model with the sulfur profiles for the paste exposed to sodium sulfate.

10

Total sulfur content (g/kg)

1750

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756 1800

1800

450

450

1600

1600

400

1000

250

800

200

600

150

400

100

350 1200

300

1000

250

800

200

600

150

400 200

50 0 0

2

4 6 Position (mm)

8

50

Simulation

Simulation 0

100

Microprobe measurement

Microprobe measurement 200

0

10

Total calcium content (g/kg)

300

1400 Microprobe (count/sec)

350 1200

400

Total calcium content (g/kg)

Microprobe (count/sec)

1400

1751

0 0

2

(a) Three month exposure

4 6 Position (mm)

8

10

(b) Six month exposure

1800 450 1600

Microprobe (count/sec)

350 1200

300

1000

250

800

200

600

150

400

100

Total calcium content (g/kg)

400 1400

Microprobe measurement 200

50

Simulation 0

0 0

2

4 6 Position (mm)

8

10

(c) Twelve month exposure

Fig. 4. Validation test – comparison of the model with the calcium profiles for the paste exposed to sodium sulfate.

only discrepancy between the model and the measurements is observed in Fig. 4a. The model predicts the presence of a calcium plateau occurring after the dissolution of portlandite. But it cannot be observed on the experimental measurement. So far, no explanation could be found for this. Another series of simulations was performed to evaluate the impact of the electrical coupling effect on the numerical results. The case of the cement paste exposed to sodium sulfate was considered. The first simulation was performed

with the materials properties listed in Table 3, but with the electrical coupling terms neglected. This gives concentration profiles that do not respect the electroneutrality requirement. Another simulation was performed using the average of the diffusion coefficients in Table 3 for each ionic species. This allows maintaining a neutral solution throughout the calculations since the initial concentrations and boundary conditions are neutral. Fig. 5 shows that neglecting the electrical coupling leads to results that do

160

350 Simulation with electrical coupling

300

Simulation without coupling

Total sulfure content (g/kg)

Total sulfure content (g/kg)

140

Simulation without coupling - Average D

120 100 80 60 40

250 200 150 100

20

50

0

0

Simulation with electrical coupling Simulation without coupling Simulation without coupling - Average D

0

2

6 4 Position (mm)

(a) Sulfur profiles

8

10

0

2

4 6 Position (mm)

8

10

(b) Calcium profiles

Fig. 5. Simulations performed without electrical coupling and chemical activity terms for the paste exposed to sodium sulfate for 1 year.

1752

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

not match with the experimental measurements. This is especially true in Fig. 5a, where the coupling strongly affects the position of the gypsum peak. Simulations made while neglecting chemical activity showed that this term has very little influence on the solution. 6. Long-term simulation Other simulations are performed to investigate the behavior of the SNIA algorithm used in the model, as well as the effect of a gradient in water content on the ionic transport. The case considered consists in a 15 cm concrete slab lying on a soil bearing a high concentration in sodium sulfate. The concrete has a water/cement ratio of 0.65 and is made with a ASTM type V cement. It is assumed that the soil is saturated (100% RH). The top of the slab is exposed to air at a 65% relative humidity. This gradient in relative humidity will create a flow of water directed upward that will contribute, along with diffusion, to move the SO2 4 ions into the material. The simulations were all performed at a 25 C temperature. The parameters needed to perform simulations were evaluated from laboratory concrete samples. The ordinary concrete was prepared without any admixtures. Table 4 shows the mixture proportions. After the mixing sequence, the material was cast in cylindrical molds (10 cm diameter, 20 cm length). The cylinders were demolded 24 h after casting and stored in a fog room (100% RH) for 28 days before testing. Table 4 also shows slump test measurements, air content and compressive strength measurements performed at 7, 28 and 91 days of hydration. After the curing period, cylinders were taken from the fog room and cut to perform the transport properties tests. The properties are listed in Table 5. The procedures performed to measure the porosity, initial pore solution composition, initial solid phase content and diffusion coefficients have already been described in the previous secTable 4 Properties of the 0.65 w/c concrete mixture Properties

Values 3

Mixture proportions (kg/m ) Cement Water Sand Coarse aggregates Density

280 182 850 1055 2367

Slump tests (mm) 10 min 30 min

45 40

Compressive strength (MPa) 7 days 28 days 91 days

26.0 33.3 42.5

Air content (%) 10 min 30 min

1.6 1.7

Table 5 Transport properties of the 0.65 w/c concrete Properties

Values

Cement type w/c

ASTM Type V 0.65

Mixture proportions (kg/m3) Cement Water Coarse aggregates Fine aggregates Density

280 182 1055 850 2367

Cement composition (% mass) CaO SiO2 Al2O3 SO3

65.3 21.6 3.6 2.1

Initial solid phases (g/kg) Portlandite C-S-H Monosulfate Ettringite

41.9 77.5 15.1 2.8

Porosity (%)

12.7 11

Diffusion coefficients (10 OH Na+ K+ SO2 4 Ca2+ AlðOHÞ 4

2

m /s) 40.6 10.3 15.1 8.2 6.1 4.2

Water diffusivity A (1011 m2/s) B (–)

1.02 34.8

Initial pore solution (mmol/L) OH Na+ K+ SO2 4 Ca2+ AlðOHÞ 4

171.1 83.4 82.7 1.1 3.7 0.2

tion. In the case of the migration tests, they were performed on 50 mm concrete samples. The transport parameters Dw (Eq. (8)) and DL (Eq. (18)) are needed to model the transport of water in the concrete slab. It was shown in Ref. [39] that for high values of relative humidity, Dw  DL. This parameter is nonlinear and is generally expressed as an exponential function for construction materials [40]: Dw = AeBw. The parameter was evaluated from a drying test, which consists in exposing disks of initially saturated material to two different relative humidities. In this case, the relative humidities were controlled by saturated NaCl (75% RH) and MgNO3 (52% RH) salt solutions in tightly closed plastic boxes. The samples are sealed on one end and on the side to enforce 1D water transport. The mass of the samples is measured over approximately 1 month. Richards’ Eq. (8) is used to reproduce the two mass loss curves by adjusting A and B. The analysis gave A = 1.02 · 1011 m2/s and B = 34.8. The

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

Before showing any simulation results, a few words should be written on the ability of the algorithm to handle the water flowing through the material as a result of capillary suction. The relative importance of the water flow compared to diffusion can be evaluated with the Peclet number (Pe). It is given by V L/D, where V is the velocity of the flow, L is the length of the domain and D is the diffusion coefficient. In the present case, V can be calculated with Eq. (9), with the gradient of water content evaluated from the analytical solution of Eq. (8) in steady-state. It gives 4.6 · 1010 m3/m2/s. D is taken as the OH diffusion coefficient, 40.6 · 1011 m2/s. This results in a low Peclet number of 0.06. The numerical problems commonly associated with convection-dominated problems, like strong oscillations, are not likely to occur at such low values of Pe. All the simulations were performed over a 20 year period. The first simulation was performed with a constant time steps of day over a uniform 75 element mesh. The solid phase profiles are shown in Fig. 6. It shows that the presence of sulfate in the groundwater leads to the penetra-

120 Portlandite

Solid phase content (g/kg)

values are reported in Table 5. Parallel to the drying tests, ISO12572 vapor transmission tests [41] were performed on two 5-cm samples. This test consists in exposing the two faces of a disk to different humidities until steady-state is reached. The moisture flux across the samples is calculated from the mass variations. For this paper, a reversed setup was used to better control the boundary conditions: the top surface of the sample was maintained in contact with water while the other face was maintained in a 50% RH chamber. The test was performed on two samples. The average steady flux for the two samples is 0.040 kg/m2/d. Using the values of A and B evaluated from the drying test to reproduce the ISO experiment gives 0.0395 kg/m2/d, which confirms the reliability of the parameters. The boundary conditions at x = 0 corresponds to the ionic concentration in the groundwater. For this particular simulation, the groundwater bears a 52 mmol/L (5000 ppm) Na2SO4 concentration. Accordingly, concentrations in Na+ of 104 mmol/L and of 52 mmol/L for SO2 are imposed through Dirichlet-type conditions. All 4 other concentrations are set to zero. At x = 0.15, the material is in contact with air. A null flux is imposed for all the species at this location. As mentioned previously, the bottom part of the slab is in contact with a saturated soil. The water content at that location is set (Dirichlet boundary condition) to the porosity value. The upper surface is in contact with a constant 65% relative humidity environment. The relative humidity is translated into water content at x = L using the desorption isotherm method proposed by Xi et al. [42]. This corresponding water content w1 is then imposed as a convection flux boundary condition: qn = hw(w  w1), with hw = 3 · 108 m/s. Values for hw are proposed in Ref. [43]. For the electrical potential w, a value of zero is set at x = 0.

1753

C-S-H

100

Ettringite Monosulfate

80

Gypsum

60 40 20 0 0

2

4

6

8

10

12

Position (cm)

Fig. 6. Solid phase distribution after 20 years.

40

Time step: 1/4 day Time step: 1/2 day

35

Time step: 1 day

Ettringite content (g/kg)

30

Time step: 2 days Time step: 5 days

25

20

15

10

5

0 0

2

4

6

8 Position (cm)

10

12

14

Fig. 7. Effect of time steps on the position of the ettringite front after 20 years. The simulations were made with 75 elements.

14

1754

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756 40

30 elements 35

75 elements 150 elements

Ettringite content (g/kg)

30

300 elements

25

20

15

10

5

0 0

2

4

6

8 Position (cm)

10

12

14

Fig. 8. Effect of the mesh density on the position of the ettringite front after 20 years. The simulations were made with time steps of 43,200 s (0.5 day).

tion of an ettringite front in the concrete slab as well as the propagation of a gypsum peak toward the top of the structure. These features were also observed experimentally on the microprobe analyses of the previous section. The position of the ettringite front matches with the dissolution front of monosulfates. The model also predicts the dissolution of portlandite and the decalcification of C-S-H, as a result of OH and Ca2+ leaching in the groundwater. This phenomenon was also observed on the hydrated cement paste submitted to the degradation tests. Finally, it is interesting to observe the presence of ettringite formation on the top of the slab (x = 15 cm). This is caused by the drying process, which leads to an increase of all concentrations as the volume of liquid water is reduced. The sensitivity of the model to time discretization was tested by doing the same simulation with different time steps: 1/4, 1/2, 1, 2, and 5 days. The results are shown in

Fig. 7. All simulations were made with 75 elements. It clearly shows the influence of reducing the time steps on the solution. The ettringite front is spread over a wide area for large time steps. As the latter decreases, the front gets sharper. The figure also shows the convergence of the front as the time steps gets smaller. Using a 1/2-day time step offers a good balance between precision and calculation time. A similar sensitivity analysis was performed for the spatial discretization. Fig. 8 shows the ettringite front location calculated with various mesh densities. A constant time step of 0.5 day was used. It clearly shows that this parameter has not a large influence on the numerical solution, contrary to the time step. The simulation results are similar with 75, 150 and 300 elements. However, for the 30-element case, some oscillations are observed. The last series of simulations was made to show the influence of the water transport on the chemical alteration

50 45

Solid phase content (g/kg)

40

Portlandite - unsaturated

35

Portlandite - saturated

30

Ettringite - unsaturated 25

Ettringite - saturated

20 15 10 5 0 0

2

4

6

8

10

12

14

Position (cm)

Fig. 9. Effect of the water flow on the ettringite and portlandite fronts after 20 years.

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

of the material. The simulation made with 75 elements and a time step of 0.5 day (Fig. 6) is compared to the same simulation without any gradient of water content. In that case, the material remains saturated. The comparison for the ettringite and portlandite fronts is showed in Fig. 9. The ettringite front penetrates faster in the material under the influence of a water content gradient, even if the diffusion coefficients are lowered according to Eq. (2). The water flow thus has a stronger influence than diffusion on the penetration of the front. On the contrary, the dissolution of portlandite remains virtually unchanged. In that case, the water flow is opposed to the diffusion flux of OH and Ca2+ ions which tend to leach out of the material. The results indicate that the flow compensates the reduction in diffusion of the fast OH ion. 7. Conclusion The numerical model presented is based on the SNIA algorithm, in which the transport and chemical parts of the ionic diffusion process are solved separately without iterations between the two. The transport equations take into account the electrical coupling between the ions and the chemical activity effects. These effects are considered because cement-based materials bear a highly charged pore solution. The model was compared to experimental results, showing a good match with the measurements. The model was able to reproduce simultaneous chemical effect like the precipitation of ettringite and gypsum, as well as the dissolution of portlandite and decalcification of C-S-H. Comparison with the experimental results emphasized the importance of electrical coupling. Simulations made without considering this term in the model showed large discrepancies with the measurements. Another series of simulations were used to test the sensitivity of the algorithm to various numerical parameters, namely mesh density and time steps. While the mesh density clearly showed no major influence on the numerical solution, the situation is different for the time step. The simulations showed that a small time step is needed in order to have a reliable result. Large time steps tend to produce wide fronts. As the time step gets smaller, the solution converges to a unique sharp front. Acknowledgement The authors would like to acknowledge the work of Dr. Yannick Maltais, who provided all the experimental results presented in this paper. References [1] Bear J, Bachmat Y. Introduction to modeling of transport phenomena in porous media. Netherlands: Kluwer Academic Publishers; 1991.

1755

[2] Samson E, Marchand J, Beaudoin JJ. Describing ion diffusion mechanisms in cement-based materials using the homogenization technique. Cem Concr Res 1999;29(8):1341–5. [3] Helfferich F. Ion exchange. USA: McGraw-Hill; 1961. [4] Samson E, Marchand J, Snyder KA, Beaudoin JJ. Modeling ion and fluid transport in unsaturated cement systems in isothermal conditions. Cem Concr Res 2005;35:141–53. [5] Millington RJ, Quirk JP. Permeability of porous solids. Trans Faraday Soc 1961;57:1200–7. [6] Zhang JZ, Buenfeld NR. Presence and possible implications of a membrane potential in concrete exposed to chloride solution. Cem Concr Res 1997;27(6):853–9. [7] Pankow JF. Aquatic chemistry concepts. USA: Lewis Publishers; 1994. [8] Samson E, Lemaire G, Marchand J, Beaudoin JJ. Modeling chemical activity effects in strong ionic solutions. Comput Mater Sci 1999;15(3):285–94. [9] Daian J-F. Condensation and isothermal water transfer in cement mortar: 1. Pore size distribution, equilibrium water condensation and imbibition. Transp Porous Media 1988;3:563–89. [10] Whitaker S. Simultaneous heat, mass, and momentum transfer in porous media: a theory of drying. Adv Heat Transfer 1977;13:119–203. [11] Rubin J. Transport of reacting solutes in porous media: relation between mathematical nature of problem formulation and chemical nature of reactions. Water Resour Res 1983;19(5):1231–52. [12] Miller CW, Benson LV. Simulation of solute transport in a chemically reactive heterogeneous system: model development and application. Water Resour Res 1983;19(2):381–91. [13] Yeh GT, Tripathi VS. A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components. Water Resour Res 1989;25(1):93–108. [14] Grove DB, Wood WW. Prediction and field verification of subsurface-water quality during artificial recharge, Lubbock, Texas. Groundwater 1979;17(3):250–7. [15] Yeh GT, Tripathi VS. A model for simulating transport of reactive multispecies components: model development and demonstration. Water Resour Res 1991;27(12):3075–94. [16] Walter AL, Frind EO, Blowes DW, Ptacek CJ, Molson JW. Modeling of multicomponent reactive transport in groundwater. 1. Model development and evaluation. Water Resour Res 1994;30(11):3137–48. [17] Xu T, Samper J, Ayora C, Manzano M, Custodio E. Modeling of non-isothermal multi-component reactive transport in field scale porous media flow systems. J Hydrol 1999;214:144–64. [18] Herzer J, Kinzelbach W. Coupling of transport and chemical processes in numerical transport models. Geoderma 1989;44:115–27. [19] Valocchi AJ, Malmstead M. Accuracy of operator splitting for advection–dispersion–reaction problems. Water Resour Res 1992;28(5):1471–6. [20] Kaluarachchi JJ, Morshed J. Critical assessment of the operatorsplitting technique in solving the advection–dispersion–reaction equation: 1. First-order reaction. Adv Water Resour 1995;18(2):89–100. [21] Barry DA, Miller CT, Culligan-Hensley PJ. Temporal discretisation errors in non-iterative split-operator approaches to solving chemical reaction/groundwater transport models. J Contam Hydrol 1996;22:1–17. [22] Samson E, Marchand J, Beaudoin JJ. Modeling the influence of chemical reactions on the mechanisms of ionic transport in porous materials: an overview. Cem Concr Res 2000;30:1895–902. [23] Liu CW, Narasimhan TN. Redox-controlled multiple-species reactive chemical transport – 1. Model development. Water Resour Res 1989;25(5):869–82. [24] Berner UR. Modelling the incongruent dissolution of hydrated cement minerals. Radiochim Acta 1988;44/45:387–93. [25] Henocq P. Modeling ionic interaction on the calcium silicate hydrate surface. Phd thesis, Laval University, Canada; 2005 [in French].

1756

E. Samson, J. Marchand / Computers and Structures 85 (2007) 1740–1756

[26] Xu T, Sonnenthal E, Spycher N, Pruess K. TOUGHREACT – A simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media: applications to geothermal injectivity and CO2 geological sequestration. Comput Geosci 2006;32:145–65. [27] Zhang T, Samson E, Marchand J. Effect of temperature on ionic transport properties of concrete. In: Proceedings of the ConMAT conference, Vancouver, Canada, August 2004. [28] MacQuarrie KTB, Mayer KU. Reactive transport modeling in fractured rock: a state-of-the-science review. Earth-Sci Rev 2005;72:189–227. [29] Lichtner PC, Pabalan RT, Steefel CI. Model calculations of porosity reduction resulting from cement-tuff diffusive interaction. Mater Res Soc Sympos 1998;56:709–18. [30] Thomas JJ, Jennings HM, Allen AJ. The surface area of hardened cement paste as measured by various techniques. Concr Sci Eng 1999;1:45–64. [31] Barbarulo R, Marchand J, Snyder KA, Prene´ S. Dimensional analysis of ionic transport problems in hydrated cement systems Part 1. Theoretical considerations. Cem Concr Res 2000;30:1955–60. [32] Knapp RB. Spatial and temporal scales for local equilibrium in dynamic fluid-rock systems. Geochim Cosmochim Acta 1989;53:1955–64. [33] Samson E, Marchand J, Snyder KA. Calculation of ionic diffusion coefficients on the basis of migration test results. Mater Struct 2003;36:156–65.

[34] Samson E, Marchand J, Robert J-L, Bournazel J-P. Modelling ion diffusion mechanisms in porous media. Int J Num Methods Eng 1999;46:2043–60. [35] Samson E, Marchand J. Numerical solution of the extended Nernst– Planck model. J Colloid Interface Sci 1999;215:1–8. [36] Reddy JN, Gartling DK. The finite element method in heat transfer and fluid dynamics. USA: CRC Press; 1994. [37] Barneyback RS, Diamond S. Expression and analysis of pore fluid from hardened cement pastes and mortars. Cem Concr Res 1981;11:279–85. [38] Planel D. Les effets couple´s de la pre´cipitation d’espe`ces secondaires sur le comportement me´canique et la de´gradation chimique des be´tons. Phd thesis, Universite´ de Marne-la-Valle´e, France; 2002 [in French]. [39] Crausse P, Bacon G, Bories S. Etude fondamentale des transferts couple´s chaleur-masse en milieu poreux. Int J Heat Mass Transfer 1981;24(6):991–1004. [40] Hall C. Barrier performance on concrete: a review of fluid transport theory. Mater Struct 1994;27:291–306. [41] ISO 12572 – Hygrothermal performance of building materials and products – Determination of water vapour transmission properties. Ref. number: ISO 12572:2001(E); 2001. [42] Xi Y, Bazˇant ZP, Jenning HM. Moisture diffusion in cementitious materials – adsorption isotherms. Adv Cem Based Mater 1994;1:248–57. [43] Sakata K. A study on moisture diffusion in drying and drying shrinkage of concrete. Cem Concr Res 1983;13:216–24.