A novel numerical multi-component model for simulating hydration of cement

Computational Materials Science 78 (2013) 12–21

Contents lists available at SciVerse ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A novel numerical multi-component model for simulating hydration of cement Nghi L.B. Le ⇑, Martijn Stroeven, Lambertus J. Sluys, Piet Stroeven Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, 2628 CN Delft, The Netherlands

a r t i c l e

i n f o

Article history: Received 18 January 2013 Received in revised form 8 May 2013 Accepted 14 May 2013 Available online 10 June 2013 Keywords: Hydration simulation Vector approach Multi-component model Microstructure Pozzolanic blending Penetration rate

a b s t r a c t A numerical multi-component model to simulate the microstructural evolution of plain cement as well as pozzolanic blended cement during hydration is presented. The model is an extension of the so-called ‘‘Integrated Particle Kinetics Model’’ (IPKM) developed by Navi and Pignat [8] to simulate hydration of tricalcium silicate grains. The IPKM is extended to a new model that accounts for different major cement compounds and for a pozzolanic admixture (in case of blended paste). Therefore, the proposed model is referred to as the ‘‘Extended Integrated Particle Kinetics Model’’ (XIPKM). New computational techniques used in the simulation are also presented. A numerical method to estimate the basic penetration rate of the particle reaction front based on experimental data is proposed. The results obtained from four simulated samples of plain cement and pozzolanic blended cement are compared with experimental data from different authors. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction The properties of cement-based materials such as stiffness, strength and transport-based durability evidently depend on its complex heterogeneous microstructure. Understanding the microstructure of cement paste is at the basis of the study on properties of cementitious materials. Experimental methods for studying the spatial microstructure, e.g., by means of serial sectioning and 3D reconstruction are unfortunately laborious, time-consuming and thus expensive. Portland cement (PC) production contributes by about 6% to global emissions of CO2. Partial replacement of Portland cement by pozzolanic mineral admixtures has been proven a possible option for reduction of such emissions. Experiments have demonstrated the profitable effects of blending Portland cement with a pozzolanic admixture that is significantly finer [1]. However, experimental research on the microstructural characteristics and on the properties of pozzolanic blended cement is not simple and even more complicated than in the case of plain cement. It is therefore more economical to conduct research on virtual representations of cement-based materials by means of modern computer facilities. Several computer-based models have been developed in the last few decades for simulating hydration and the evolution of microstructure of cement-based materials. A review of such models has been given in [2], whereby the models have been classified into ⇑ Corresponding author. Tel.: +31 15 278 3332; fax: +31 15 278 5767. E-mail address: [email protected] (N.L.B. Le). 0927-0256/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.05.021

two categories. The first is the so-called ‘‘discretization approach’’ where models are based on a digital-image processing approach. The microstructure is represented by a matrix of ‘‘voxels’’ (volumetric pixels), each of which represents a relevant phase. This approach allows a direct representation of a multiphase material with grains of any shape. CEMHYD3D [3–5] is one of the most popular models using the discretization approach for simulating the microstructure of (pozzolanic blended) cement paste during hydration. It can cope with the main compounds of the cement and the pozzolanic admixture. The second is the so-called ‘‘vector approach’’ (continuum approach), since it does not require a subdivision of the microstructure. Each cement particle is described by its position and a set of radii representing an unhydrated spherical core that is covered by spherical shells of hydration products. The vector approach was first used by Jennings and Johnson [6]. The original version of HYMOSTRUC [7], the ‘‘Integrated Particle Kinetics Model’’ (IPKM) [8–11], SPACE [12] and lic [2] are based on this approach. These models all assume particles to be spherical. In addition, these models were restricted to simulating hydration of single phase cementitious grains. CCBM [13] and an extended version of HYMOSTRUC3D [14] are more recent models using the vector approach whereby the main cement compounds are taken into account. The latter one also considers the pozzolanic reaction of rice husk ash (RHA) in case of a blended binder. The cement particles are modelled similarly in these two models as in the aforementioned ones, but the hydration rate of the unhydrated spherical core is homogenized by averaging the hydration rates of the main compounds.

13

N.L.B. Le et al. / Computational Materials Science 78 (2013) 12–21

Nomenclature b, n B Cs Cv Cw fo i j Ko p ro rin rout VCH Vm,0 Vm,reac Vw,0 Vw

Rosin–Rammler parameters Blaine fineness free surface ratio free volume ratio reduction coefficient due to water decrease volume proportion of the compound particle index time step index basic penetration rate of the hydration front volume proportion of the compound initial radius of the particle inner radius of the spherical sector representing the unhydrated compound outer radius of the particle CH volume at time of hydration total volume of cement (and pozzolan) at initial time total volume of the hydrated cement (and pozzolan) water volume at initial time water volume at time of hydration

A numerical model that is developed in this article is based on the vector approach for simulating hydration and microstructural evolution of PC as well as of pozzolanic blended PC with a multicompound concept. This model is denoted the ‘‘Extended Integrated Particle Kinetics Model’’ (XIPKM) because the chemical reactions, the particle model and the hydration kinetics of the IPKM are extended for taking the major cement compounds and the pozzolan into account. In this study it is assumed that the cement grains consist of four main minerals, i.e., tricalcium silicate (C3S), dicalcium silicate (C2S), tricalcium aluminate (C3A) and tetracalcium aluminoferrite (C4AF), whereas the pozzolanic grains are composed of only silica (SiO2) and an inert phase. The particles expand during hydration by the precipitation of hydration products, thereby gradually invading the pore space. Accurate estimation of particle expansion in a system of complex and implicit interfering particles is important for properly describing the microstructural evolution. A numerical method is proposed for this purpose. Besides the vector-based data of particles during hydration, the hydrating microstructure in this model is also represented and visualized by a 3D digital-image-based ‘‘voxel system’’. This voxel data can be used as the input for further studies on microstructure analysis and for predicting properties of hydrating cementitious materials. In addition, this article also proposes a numerical approach to estimate the basic penetration rate of the particle’s ‘‘reaction front’’ (one of the main parameters of the model) based on experimental data. Finally, the model is validated by comparing numerical results with experimental data obtained by other authors. 2. Formation of the hydration model 2.1. Chemical reactions Hydration of hydraulic materials involves a chemical reaction with water. The volume change of each phase (cement components, pozzolan, free water and hydration products) can be computed via the chemical equilibrium of cement reactions. In the present model, the cement particles are assumed to be composed of four major cement compounds, i.e., C3S, C2S, C3A and C4AF. The chemical equilibrium reaction can be expressed according to Bentz et al. [3–5] in relative volume ratios of water and hydration products to reactants, neglecting the gypsum reaction:

Vprod wo w/c w/b

a b d dtr

c k k

m q n

total volume of hydration product initial water to cement ratio water to cement ratio water to binder ratio degree of hydration parameter to control the hydration stage thickness of the hydration product shell transition thickness of the hydration product shell at a change of mechanism parameter to control the hydration rate index for the hydration products, water and CH as reactant compound index relative volume ratio in chemical reaction specific mass of cement penetration depth of the compound

1:0V C3S þ 1:34V H ! 1:521V CSH þ 0:610V CH

ð1Þ

1:0V C2S þ 1:49V H ! 2:077V CSH þ 0:191V CH

ð2Þ

1:0V C3A þ 1:21V H ! 1:69V CAH

ð3Þ

1:0V C4AF þ 1:41V H ! 1:17V CAH þ 0:259V CH þ 0:545V FH

ð4Þ

with CSH, CH, CAH and FH being calcium silicate hydrate, calcium hydroxide, hydrogarnet and iron hydroxide, respectively. Similarly, Bentz et al. [15] give for the silica (SiO2) in the pozzolanic admixture:

1:0V S þ 1:35V CH þ 1:87V H ! 3:77V CSH :

ð5Þ

It is noted that in these equations the standard cement chemistry abbreviations are used: C = CaO, S = SiO2, A = Al2O3, F = Fe2O3 and H = H2O and the corresponding valences are not shown for simplification. For details, see [3–5,15]. It is well know that in practice the PC clinker is usually ground with a minor amount of gypsum. In this case, C3A reacts with gypsum to form ettringite, crystalizing in the needle shape at the surface of the cement grains. When the gypsum has been consumed, the ettringite starts decomposing and thereupon reacts with additional C3A to form the monosulphate phase [3]. Since the needle-shaped ettringite is formed only during the first hours of hydration and the amount of gypsum is small, the gypsum and the ettringite formation are ignored in this model. 2.2. New particle model In XIPKM, the cement and pozzolan grains are modelled by the so-called ‘‘multi-component and multi-layer spherical particle model’’. Before hydration, each fresh particle is assumed to be a multi-component spherical particle, partaken by the constituent components formed by spherical sectors, as illustrated at the top of Fig. 1. During hydration the unhydrated (multi-component) core of each particle reacts with water – thereby decreasing in size – and the reaction products partly precipitates on the outer surface of the grain itself and in the region that was occupied by the just reacted core. Fig. 1 (bottom) shows the model that describes the state of a hydrating particle. The assumption that the diffusing CSH concentrates and locates at the dissolution source due to low mobility of silicate ionic species is in accordance with Bentz [3]. The core’s components are assumed to maintain their spherical

14

N.L.B. Le et al. / Computational Materials Science 78 (2013) 12–21

C2S

ro

CA CAF

unhydrated state (t = 0)

SiO2

C3S ro CSHin

CH ro

rout ro

hydrated state (t > 0) CAH

rinλ

rinλ

δλ

CSHout

cement λ = C3S, C3S, CA, CAF

δ

rout

λ

FH

pozzolan λ = SiO2, inert

Fig. 1. Particle models of cement, pozzolan and hydration products in the unhydrated state (top) and hydrated state (bottom).

sector shape, however, due to the different hydration rates the radii of the spherical sectors describing the surface of the unhydrated product may differ. The surfaces of the spherical sectors are assumed to be the places where the free water (diffusing through the CSH shell) reacts with the components. These sector surfaces are referred to as the ‘‘reaction fronts’’. The penetration rates of the reaction fronts are assumed controlled by a hydration kinetic model presented in the next section. The shell of the CSH product is composed of two distinct layers, the ‘‘inner CSH’’ (CSHin) that substitutes the volume of the hydrated components and thus grows inwards and the ‘‘outer CSH’’ (CSHout) that precipitates spherically on the outer surface of the particle and thus grows outwards. When the outer CSH layer of the particle interferes with those of neighboring particles, the outer CSH no longer precipitates at the overlap area, but is evenly re-distributed onto the remaining free outer surface of the particle. A model particle is configured by geometrical data (i.e., ro, rkin , rout and coordinates) and fractions of the components, where k denotes the different components (Fig. 1). r inert of the spherical sector characterizing the inert phase in of the pozzolanic particles is assumed unchanged during hydration and thus will be neglected in further considerations. CH, CAH and FH diffuse and either nucleate randomly into crystals in pore space or precipitate on surfaces of already existing crystals [3,4,11]. These species are modelled by simple spherical particles, illustrated at the bottom (three particles on the right hand side) of Fig. 1. The decision to nucleate a CH particle is based on the number of expected grains at time t. In [8] Navi and Pignat propose the following rule:

nCH ðtÞ ¼ nCH; max ð1  eat Þ

ð6Þ

where nCH,max is the maximum number of CH grains at complete hydration and a is a constant. The two constants can be calculated following Jennings and Parrott’s work [16] reporting that there is one CH grain for 12.5 grains of C3S after one day of hydration and five grains of C3S at complete hydration. The size of a CH particle either increases during the hydration process due to the precipitation of the CH produced by the cement reaction, or decreases as a result of the pozzolanic reaction of CH. The number of researches on diffusion and nucleation of CAH and FH is scarce; however Bentz [3,4] shows that the probability of nucleation of diffusing CAH and FH is governed by an equation with a form similar to Eq. (6). Therefore, the nucleation of these species can be described by the following general equation

nk ¼ nk;max ð1  eakt Þ k ¼ CH; CAH; FH:

ð7Þ

It is noted that with higher mobility of calcium and aluminate ions, CH and CAH are distributed somewhat uniformly, whereas with lower mobility of ferrite ions, FH crystallizes near to a dissolution source [3]. CH and CAH that are produced by hydration of a particle are therefore assumed uniformly at random distributed in pore space in the present model, whereas FH is assumed distributed nearby the particle specified by a surrounding sphere with a given radius. 2.3. Hydration kinetic model Three different mechanisms were introduced in the IPKM model by Navi and Pignat [8] for controlling the kinetics of the reaction front of a C3S particle, i.e., ‘‘nucleation and growth’’, ‘‘phase-boundary’’ and ‘‘diffusion’’. The first mechanism is neglected, since it only controls the initial reaction until a relative small degree of hydration of 1–2% [8]. The penetration rate of the reaction front of a C3S particle i for the other two mechanisms can be given by a single formula: i

drin ðtÞ ¼ K o ðdtr =di ðtÞÞb C is ðtÞC w ðtÞ dt

ð8Þ

where riin ðtÞ is the radius of the fresh (unhydrated) core at time t, Ko (lm/h) is the basic rate constant and di is the thickness of the product shell that is equal to rout–rin, where rout is the outer radius of the hydration product. dtr is the transition thickness equal to di at the time that the controlling mechanism changes from the ‘‘phase-boundary’’ stage (b = 0 and dtr/di 6 1) to the ‘‘diffusion’’ stage (b = 1 and dtr/di > 1). C is is a reduction coefficient related to the decrease in free surface (parts in contact with the water or pore) and calculated as the area ratio between the free surface and the total surface (assuming no contact) of the particle. Cw is a global reduction coefficient given by Rohling and cited by van Breugel [7] related to the decrease in water content and calculated by

C w ðtÞ ¼ 1 

v aðtÞ qwo þ aðtÞ

ð9Þ

where v is a constant defined as the volume ratio of the total hydration products (i.e., CSH and CH) and the hydrated amount of C3S, a is the degree of hydration defined as the volume ratio of the hydrated amount and the initial amount of C3S, q (g/cm3) is the specific mass of the fresh material and wo is the initial water-to-binder ratio. The hydration rate in XIPKM is established separately for the reaction front of each mineral component by the following formula

15

N.L.B. Le et al. / Computational Materials Science 78 (2013) 12–21

"

k;i

drin ðtÞ ¼ K ko dt

dktr

#c !b C is ðtÞC w ðtÞ

k;i

d ðtÞ

k ¼ C3 S;C2 S;C3 A;C4 AF;SiO2

ð10Þ

in which the meaning of the symbols is the same as in Eq. (8), except that k is attached to denote the respective components and c is added to have one more parameter to control the hydration rate. dk;i is the thickness of the hydration product layer covering component k that equals riout  rk;i in (Fig. 1). As different components are included, the volume ratio between the total hydration product to the total hydrated amount (v in Eq. (9)) is no longer a constant and thus Eq. (9) for CW is not valid in this case. wo in Eq. (9) equals Vw,0/ (qVm,0) and v and a can be defined as Vprod/Vm,reac and Vm,react/Vm,0, respectively. Vw,0, Vm,0, Vm,reac and Vprod represent the initial water volume, the initial binder volume, the hydrated binder volume, and the total hydration product volume, respectively. This renders possible rewriting Cw in the more general form

C w ðtÞ ¼ 1 

V prod ðtÞ V w;0 þ V m;reac ðtÞ

ð11Þ

When the time of hydration is discretized in a number of time steps and the hydration rate is assumed to be constant during each step, the penetration depth of the reaction front of component k of particle i during time step j can be estimated by

" k Dr k;i in;j ¼ K o

dktr

dk;i j

#c !b C is;j C w;j Dtj k ¼ C3 S;C2 S;C3 A;C4 AF;SiO2 :

ð12Þ

The hydration simulation during a time step is in many vector approach models [2,8,12] an order-based process in which hydration proceeds successively from particle to particle. The hydration rate of a particle is proportional to the remaining water volume and the free surface ratio and, as a result, depends on water consumption and volume expansion by hydration of previous particles. So, the hydration rate of a particle depends on the order of the particles. When the system is out of water for the hydration of a particle, for example, the hydration no longer proceeds to pending particles. To avoid this order-based process in the present study, hence to simulate hydration concurrently, the coefficients in Eq. (12) are estimated at the beginning of step j for all particles in the presented model. The amount of water available for each particle’s hydration is pre-defined as well; the total remaining water is divided into parts, the number of which equals the number of particles and each of which is only available for hydration of a single particle. It is assumed that the proportion between the volumes of these water parts equals that the between water consumption rates of the particles. Hence, the water amount available for hydration of a particle can be estimated by

P

i;k k ¼ C3 S;C2 S;C3 A;C4 AFðcementÞ k K w;j V iw;j ¼ P P i;k V w;j k ¼ SiO2 ðpozzolanÞ K i k w;j

ð13Þ

where Vw,j is the total volume of water available at the beginning of time step j. K i;k w;j represents the water consumption rate of a component that is obtained by multiplying the surface area, the penetration depth and the relative water ratio (in the chemical equations) of the particular component, divided by the time increment. Hence, 0

2 k k;i k i K i;k w;j ¼ 4pðr out;j Þ fo Dr in;j v w =Dtj

k ¼ C3 S;C2 S;C3 A;C4 AF;SiO2

ð14Þ

where f0k is the volume fraction as well as the surface fraction (of the 0 spherical sector) of component k at the initial state, Dri;k in;j is the penetration depth calculated by the right-hand side of Eqs. (12) and v kw is the relative volume ratio of water in the reaction equation

of component k (Eqs. (1)–(5)). It is noted that the hydration of the components of a particle also need to be simulated simultaneously. Therefore, the amount of water available for a particle is subdivided into sub-parts, the number of which equals the number of components and each of which is only available for hydration of a single component. Hence, the water required for hydration of a component of a particle during a single time step cannot exceed the sub-part available for this component. This is expressed by i K i;k k ¼ C3 S;C2 S;C3 A;C4 AF;ðcementÞ 4p h k;i 3 w;j k;i 3 k k ðr in;j Þ  ðr k;i V iw;j in;j  Dr in;j Þ f0 v w 6 P i;k 3 k ¼ SiO2 ðpozzolanÞ K k w;j ð15Þ

which yields a constraint for the penetration depth: k;i k;i 3 Dr k;i in;j 6 r in;j  ððr in;j Þ 

3V iw;j K i;k k ¼ C3 S;C2 S;C3 A;C4 AFðcementÞ w;j Þ1=3 : P k ¼ SiO2 ðpozzolanÞ 4pf0k v kw k K i;k w;j ð16Þ

Since CH is also a reactant in the pozzolanic reaction, the available CH volume also constrains the penetration depth. Hence, 3 k;i Drk;i ðr k;i in;j 6 r in;j  in;j Þ 

3V iCH;j

!1=3 k ¼ SiO2

4pf0k v kCH

ð17Þ

2 where v SiO CH is the relative volume ratio of CH in the reaction equation of SiO2 (Eq. (5)) and V iCH;j is the CH volume available for the pozzolanic reaction of particle i that can be estimated by

K i;k CH;j V iCH;j ¼ P i;k V CH;j i K CH;j

k ¼ SiO2

ð18Þ

where VCH,j is the total CH volume at the beginning of time step j and K i;k CH;j is the CH consumption rate of a component of the particle that can be estimated by 0

2 k k;i k i K i;k CH;j ¼ 4pðr out;j Þ fo Dr in;j v CH =Dtj

k ¼ SiO2 :

ð19Þ Drk;i in;j

Besides the above constraints, the penetration depth cannot exceed the current radius rk;i in;j . In summary, the penetration depth of the reaction front of component k of particle i during time step j is given by 0

00

000

k;i k;i k;i k;i Dr k;i in;j ¼ minfDr in;j ; Dr in;j ; Dr in;j ; r in;j g

k ¼ C3 S;C2 S;C3 A;C4 AFðcementÞ k ¼ SiO2 ðpozzolanÞ ð20Þ

0 Drk;i in;j ,

00 Dr k;i in;j

000 Dr k;i in;j

where and are estimated by the right-hand sides of Eqs. (12), (16), and (17), respectively. 2.4. Calculation of volume change The volumes of the consumed water, of CSH, CAH, FH, of the produced CH (CHprod) (cement hydration) and of the reacted CH (CHreact) (pozzolanic hydration) caused by the hydration of particle i during time step j can be calculated by DV ij;j ¼

i k ¼ C3 S;C2 S;C3 A;C4 AFðcementÞ 4p Xh k;i 3 k;i 3 k k k ¼ SiO2 ðpozzolanÞ ððr in;j Þ  ðr k;i in;j  Dr in;j Þ Þf0 v j 3 k j ¼ water;CSH;CAH;FH;CHprod ;CHreact ð21Þ

where v kj is the relative volume ratio corresponding to phase j by hydration of componentk, derived from the chemical equations (Eqs. (1)–(5)). Since part of the CSH volume (CSHin) substitutes the total volume of the components that has reacted in this step, the CSH volume (CSHout) precipitating on the surface of the particle is given by

16

N.L.B. Le et al. / Computational Materials Science 78 (2013) 12–21

DV iCSH;out;j ¼

i k ¼ C S;C S;C A;C AFðcementÞ 4p Xh k;i 3 3 2 3 4 k;i k;i 3 k ððr in;j Þ  ðrin;j  Dr in;j Þ Þf0 ðv kCSH  1Þ : 3 k k ¼ SiO2 ðpozzolanÞ ð22Þ

3. Algorithmic aspects 3.1. Estimation of free surface In the present model, the effect of the free surface on the hydration rates of the reaction fronts is explicitly taken into account by the area ratio of the free surface and the total surface, C is , determined for individual particles. To estimate the free surface ratio C is , the ‘‘sampling point method’’ [12] is used. A number of sampling points are ‘‘evenly’’ dispersed on the total surface of a particle. Then, the free surface ratio is numerically estimated as follows

C is ¼ 1  n=ntot

ð23Þ

where ntot is the total number of the sampling points and n is the number of the points inside other particles. The ‘‘general spiral points’’ method, which has been demonstrated to be one of the most effective and simplest ways for dispersing points evenly on a spherical surface [17,18], is used in this model. The coordinates of point k of N spiral points on a unit sphere’s surface is given by

qffiffiffiffiffiffiffiffiffiffiffiffiffi 1  z2k cosð/k Þ qffiffiffiffiffiffiffiffiffiffiffiffiffi yk ¼ 1  z2k sinð/k Þ

xk ¼

k1 zk ¼ 1 þ 2 ;1 6 k 6 N N 1 0 1 3:6 1 C B /k ¼ @/k1 þ pffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiAðmod2pÞ; 2 6 k 6 N  1; /1 ¼ /N ¼ 0: N 1  z2 k ð24Þ 3.2. Estimation of particle expansion The increment of the outer radius Driout of the particle i resulting from a new product volume that has precipitated on the particle’s surface governs the expansion of the particle. The relationship between the precipitating volume DV iprec and the corresponding radius increment Driout is given by

DV iprec ¼

4p i ððr out þ Dr iout Þ3  ðr iout Þ3 ÞC iv ðDr iout Þ 3

ð25Þ

where C iv is the ratio between the free volume (the part that has no contact with other particles) and the total volume of a spherical shell with the thickness equal to Driout (Fig. 2). Eq. (25) is an implicit

equation and cannot be solved directly for the unknown Dr iout . If Cv can be estimated from a value of Driout , nonetheless, the solution of Eq. (25) can be obtained by a standard numerical method for solving a nonlinear equation. Still, estimation of C iv ðDriout Þ could be impossible by analytical means due to complex interferences of the spherical shell with the surrounding particles. Therefore, the ‘‘sampling point method’’ is used again:

C iv ðDr iout Þ ¼ 1  n=ntot

where ntot is the total number of sampling points uniformly at random dispersed in the spherical shell with thickness Driout , and n is the number of the points situated inside another particle as illustrated in Fig. 2. Herein, the coordinates of a sampling point k are given by

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  ð2u  1Þ2 cosð2pv Þððriout Þ3 þ wððr iout þ Driout Þ3  ðr iout Þ3 ÞÞ1=3 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi yk ¼ 1  ð2u  1Þ2 sinð2pv Þððr iout Þ3 þ wððr iout þ Dr iout Þ3  ðriout Þ3 ÞÞ1=3 xk ¼

zk ¼ ð2u  1Þððr iout Þ3 þ wððr iout þ Driout Þ3  ðr iout Þ3 ÞÞ1=3 ð27Þ where u, v and w e (0, 1) are the random variables. 3.3. 3D microstructure by the voxel system The microstructure simulated by the hydration model is represented and visualized by the ‘‘digital-image-based’’ technique. The simulated space is subdivided into small cubes, called ‘‘voxels’’ (3D equivalent of image pixels), each of which is used to represent a certain phase. Herein, the microstructure of the hydrating paste is composed of the unhydrated cement, the unhydrated pozzolan, CSHin, CSHout, CH, CAH, FH and the pore space. To simplify the visualization, the constituent phases of the unhydrated core of each particle are merged into a single ‘‘average’’ phase, represented by a spherical core of which the volume equals the total volume of the individual phases. Fig. 3 is an illustration of a microstructure visualized by the voxel system, whereby the red and yellow represent the unhydrated cement and unhydrated pozzolan, respectively. Some voxels can change phases during the cement hydration process; unhydrated cement or unhydrated pozzolan voxels can change into CSHin ones, pore voxels can change into CSHout, CH, CAH or FH ones, and CH voxels can change into pore ones (pozzolanic hydration). It is noted that this voxel system is only used for representing the microstructure of the hydrating paste and thus does not affect the hydration simulation process or calculation speed. An advantage of this voxel system is its ability to precisely identify the phases that are located at the overlap of the growing particles (of different phases) during the hydration evolution, which is

uniform at random sampling points surrounding particles

spherical shell with i thickness Δrout

ð26Þ

free points points in other particles

i rout

Fig. 2. Sampling point method for estimation of the free volume ratio of a spherical shell.

17

N.L.B. Le et al. / Computational Materials Science 78 (2013) 12–21 n X wi

unhyrated cement unhydrated pozzolan

i¼1

ðr io Þ3

ðnk Þ3  3

n X wi i¼1

ðr io Þ2

! ðnk Þ2 þ 3

n X wi i¼1

r io

! nk  ak ¼ 0: ð29Þ

When the DOH of a particular compound, e.g. aC3S ref , at time tref is C3S known, taking the integral of drin ðtÞ in Eq. (10) from t = 0 to t = tref yields the following equation

CSH in CSH out

nC3S ref ¼

CH CAH FH

Fig. 3. Visualization of the hydrating paste with the phase voxel system. Each unhydrated multi-phase core is represented by a sphere of which the volume equals the total volume of the phases.

very difficult in the vector approaches. Indeed, with a vector approach, only the boundaries between phases are accurately represented, not the phases itself. Phase representation at the layer-bylayer overlap among neighboring expanding particles cannot be uniquely represented, as illustrated in Fig. 4a. On the contrary, with the recording-phase voxel system, the phases at the complex overlap zone can be simply characterized as in Fig. 4b. Additionally, the voxel data are also useful for the analysis of the microstructural characteristics as well as for the transformation to a finite element method mesh for evaluation of mechanical properties. 3.4. Numerical estimation of the basic penetration rate of hydration front Model coefficient K ko reflects the penetration rate of the hydration front of component k without considering any of the reduction coefficients. A numerical calculation is proposed in this study to obtain the basic rate constants from the hydration data of physical experiments. For that purpose a real cement paste is modelled according to the multi-component particle model presented in Section 2.2. It should be noted that the particle size range and the particle size distribution (PSD) of the model cement do not have to resemble those of the real cement. Under the assumption that, at early hydration, the penetration rate of each component is the same for all particles, the relationship between the degree of hydration (DOH) ak and the penetration depth nk of component k is given by the following equation [19]

ak ¼ 1 

!

 3 n X nk wi 1  i ro i¼1

ð28Þ

where n is the number of particles, r io is the initial radius of particle i and wi is the weight fraction derived from the PSD function corresponding to size di ðdi ¼ 2r io Þ. nk can be obtained by solving the following equation that is obtained after rearranging Eq. (28)

Z

tref

0

K C3S o C s ðtÞC w ðtÞdt

k C3S C3S where nC3S ref is calculated by solving Eq. (29) with a ¼ aref . aref and tref can be taken from the experimental hydration curve of C3S. It is noted that this calculation is taken at early hydration, so the mechanism change is neglected (b in Eq. (10) equals 0). The penetration rate of each component is assumed to be the same for all particles, so a unique global coefficient Cs, accounting for the free surface reduction should be included instead of individual values C is . This overall coefficient can be obtained from the area ratio of the total free surfaces and the total spherical surface of all particles:

Pn Cs ¼

2 i i i¼1 ðr out Þ C s ðtÞ : Pn 2 i i¼1 ðr out Þ

ð31Þ

When tref is discretized into a number of even steps, Eq. (30) transforms into

nC3S ref ¼

X C3S K o C w;j C s;j Dt

ð32Þ

j

where j denotes the step index and Cw,j can be calculated by Eq. (11) P PP with Vm,react equal to k akj fok V c;0 and Vprod equal to k j akj fok V c;0 v jk , k where fo , Vc,0 and v jk are the initial volume proportion of component k, the initial volume of cement and the chemical ratio of product j (see Eq. (21)), respectively. akj denotes the corresponding DOH of the components at time tj. akj can be taken from experimental hydration curves of the considered cement. To obtain Cs,j, a cubic volume of the model paste is taken, whereupon hydration simulation is implemented with the same penetration rate for all components. Then, the values of Cs (by Eq. (31)) and the overall DOH of the model paste at all hydration steps are recorded as a data set for interpolation. Cs,j is then obtained by interpolation at the overall DOH P aglobal that equals k akj fok . With transposition of some terms, Eq. (32) j becomes

K C3S ¼ o

NnC3S P ref : tref j C w;j C s;j

ð33Þ

where N is the number of time step (tref = NDt). The basic penetration of the other components (i.e., C2S, C3A and C4AF) is estimated in a similar way. 4. Examples For the hydration simulation of a cementitious paste, an initial packing is required of (blended) cement particles representing the fresh state of the paste. An advanced dynamic DEM system (HADES) [20–22] is utilized for this purpose. HADES allows for si-

OR

(a) vector approach

ð30Þ

(b) voxel system

Fig. 4. Phase representation at layer-by-layer overlaps of two neighboring particles during the expansion due to hydration of the particles.

18

N.L.B. Le et al. / Computational Materials Science 78 (2013) 12–21

Table 1 Input data for estimating the basic rate Ko. Refs.

El-Hemaly Maruyama et al. [13] et al. [23]

Maruyama et al. [13]

Escalante-García and Sharp [24]

Odler and Odler and Odler and Schüppstuhl [25] Schüppstuhl [25] Schüppstuhl [25]

w/c Blaine fineness (cm2/g) fo-C3S:C2S:C3A:C4AF Component to estimate aref-C3S;C2S;C3A;C4AF tref (h)

0.25 3500 0.94:0:0:0 C3S 0.36;–;–;– 24

0.5 3430 0.298:0.594:0.023:0.085 C3S;C2S 0.62;0.117;0.135;0.056 24

0.5 3765 0.605:0.218:0.084:0.093 C3S;C3A;C4AF 0.541;0.375;0.43;0.11 24

0.7 3000 1:0:0:0 C3S 0.133;–;–;– 12

0.5 3300 0.696:0.111:0.106:0.087 C3S;C2S 0.530;0.12;0.236;0 24

0.7 4700 1:0:0:0 C3S 0.194;–;–;– 12

0.7 7000 1: 0:0: 0 C3S 0.223;–;–;– 12

Size range (lm): 3–30; Rosin–Rammler coefficients: n = 1.107 and b = 0.023.

mulating particle packing under the influence of external forces (drag force or gravity). Mechanical interaction in HADES is based on a contact mechanism algorithm. The hydration model is validated by comparing the hydration curves of plain PC and of rice husk ash (RHA) blended PC specimens of simulations with those of experiments. Before proceeding to the hydration simulations, however, the basic penetration rates of C3S, C2S, C3A and C4AF components have to be defined. 4.1. Estimation of the basic penetration rates In this part, the basic penetration rates (Ko) of mineral components are obtained from a series of experimental hydration data. The input data for estimating Ko of the components is given in Table 1, whereby fo is the initial volume proportion of the components. aref and tref denote the reference values of the DOH and the corresponding time (Section 3.4). All samples are represented by model pastes with the same particle size distribution (PSD) function (Rosin–Rammler) and the same particle size range (PSR) as in the last row of Table 1. The values of Ko of the components are calculated by means of the procedure presented in Section 3.4. Ko-values for C3S are plotted as functions of component proportions and Blaine fineness in Fig. 5. The left plot shows Ko-values as a function of the C3S proportion by volume (0.298, 0.605, 0.696 and 0.94) taken from the samples of which the fineness is close to 3500 cm2/g. The right plot presents Ko-values as a function of (Blaine) fineness (3000, 4700 and 7000 cm2/g) taken from the samples having the same C3S proportion equal to 1. It is noted that the samples with different w/c ratios are modelled by different model pastes to account for the difference in Cw and Cs at the same DOH. Therefore, these Ko-value functions encompass the influence of the w/c ratio. When the experimental samples are modelled with PSR as well as PSD functions similar to those of the experimental ones – so, the fineness of the experimental and model pastes are the same – the Ko-value function is expected to be independent of the fineness of the real paste. In the simulation, nevertheless, the fineness (number of particles) of the model paste is usually lower than that of real paste to reduce computational effort. It can be seen in the samples with a fineness of 3000, 4700 and 7000 cm2/g that a higher fineness results in a higher hydration rate (Table 1). Note that these samples are modelled by the same model paste, so the Ko-value increases with higher fineness to adjust the hydration rate of the model pastes to those of the experiments. The functions for Ko of C3S can be estimated through linear regression analysis by means of

K C3S o ðpÞ ¼ 0:06485 þ 0:05168  ð1  pÞ K C3S o ðBÞ ¼ 0:05194 þ 0:01249  ðB=1000  3:5Þ where B is the Blaine fineness (cm2/g) and p is the proportion of C3S. Note that the function of K C3S o ðpÞ is determined from the samples of which the fineness is close to 3500 cm2/g. A question put forward here is how K C3S o ðpÞ is estimated for a sample of which B is different from 3500 cm2/g. By normalizing the function K C3S o ðBÞ at B equals to

3500 cm2/g and then coupling the normalized function with K C3S o ðpÞ, a unique function K C3S o ðp; BÞ can be obtained by

K oC3S ðp; BÞ ¼ ½0:06485 þ 0:05168  ð1  pÞ  ½1 þ 0:01249=0:05194  ðB=1000  3:5Þ The Ko-values for C2S, C3A and C4AF are computed as 0.0127, 0.0639 and 0.0130, respectively. Because of a lack of data these Ko-values cannot be determined as a function of p. Note that this value is computed from the sample having B nearly equal to 3500 cm2/g. In the similar way to take the B into account as in the function for Ko of C3S, the functions for Ko of C2S, C3A and C4AF are given by, respectively

K oC2S ðBÞ ¼ 0:0127  ½1 þ 0:01249=0:05194  ðB=1000  3:5Þ K oC3A ðBÞ ¼ 0:0639  ½1 þ 0:01249=0:05194  ðB=1000  3:5Þ K oC4AF ðBÞ ¼ 0:0130  ½1 þ 0:01249=0:05194  ðB=1000  3:5Þ: Hereafter, these functions can be applied to calculate basic penetration rates for the model cement paste (with size range and PSD as in Table 1) that is used to simulate a real cement paste with specific values of phase proportions and fineness. 4.2. Comparison between model and experiment results To validate the model, the experimental data of four plain PC and RHA-blended PC samples are used as input data for the hydration simulation. The input data of the samples are detailed in Table 2 with the values Ko calculated by the procedure presented in Section 4.1. Fig. 6 depicts 2D sections of the hydrated samples revealing their simulated microstructure. The sections at the top demonstrate that at the same time of hydration (28 days), the DOH of the sample with w/c = 0.4 exceeds that of the sample with w/ c = 0.2, so that the unhydrated cement fraction (color1 coded in red) is more significantly reduced. Further, porosity of the sample (i.e., the pore fraction that is color coded in black) with w/c = 0.25 is much lower than that of the sample with w/c = 0.4. The sections at the bottom of Fig. 6 shows that the CH fraction of the sample (color coded in green) is reduced at higher blending dosage as a result of the higher amount of available SiO2. Comparisons of the degree of hydration of plain PC samples shown in Fig. 7 indicate good agreement between the experimental results and the simulations. Note that trial values are used for the model parameters, i.e., dktr and c (Eq. (12)), until the simulation curves fit well to experimental data as in Fig. 7. Although a calibration process is still required, the calibration for Ko-values as the most significant model parameter can be eliminated as a result of the proposed estimation of Ko-values. The total porosities of the RHA-blended PC samples obtained by simulation and from experimental results of Nguyen [14] are com1 For interpretation of color in Fig. 6, the reader is referred to the web version of this article.

19

N.L.B. Le et al. / Computational Materials Science 78 (2013) 12–21

Fig. 5. Calculated values of penetration rate Ko of C3S particles.

Table 2 Input data of samples used for hydration simulation. Sample Refs.

W40Pc Frigione and Marra [26]

W25Pc Danielson [27]

W40Po10 Nguyen [14]

W40Po20 Nguyen [14]

w/b Cement Blaine (cm2/g) % blended RHA fo-C3S:C2S:C3A:C4AF fo-SiO2 Model cement PSD & size range

0.40 3200 0 0.542:0.278:0.064:0.117 – n = 1.107 b = 0.023 3–30 lm – 2027 0 0.0707;0.0109;0.0544;0.0111 –

0.25 3120 0 0.649:0.193:0.081:0.078 – n = 1.107 b = 0.023 3–30 lm – 2562 0 0.0821;0.0118;0.0593;0.0120 –

0.40 4500 10 0.714:0.101:0.097:0.088 0.875 * 3–30 lm

0.40 4500 20 0.714:0.101:0.097:0.088 0.875 * 3–30 lm

*

*

Model RHA PSD Number of cement particles Number of RHA particles Ko-C3S;C2S;C3A;C4AF (lm/h) Ko-RHA (lm/h)

3–30 lm 2050 883 0.0925;0.0148;0.0742;0.0151 0.0060

3–30 lm 1784 1726 0.0690;0.0084;0.0596;0.0169 0.0060

n and b are Rosin–Rammler’s parameters. PSD functions of cement and RHA are experimentally measured and plotted in [14].

*

unhydrated cement unhydrated pozzolan

CSH in W40Pc at 28 days

W25Pc at 28 days

CSH out CH CAH FH pore

W40Po10 at 90days

W40Po20 at 90 days Fig. 6. 2D sections of the samples.

20

N.L.B. Le et al. / Computational Materials Science 78 (2013) 12–21

Fig. 7. Hydration curves of plain cement samples. Fig. 10. Error in numerically obtained porosity as compared to an exact calculation.

to exact values for porosity as shown in Fig. 10. The numerical porosity is derived from the voxel system and equals the fraction of the number of pore voxels on the total number of voxels. The exact porosity is calculated by (Vtotal  Vsolid(t))/Vtotal, where Vtotal is the total volume of the sample and Vsolid(t) is the total volume of the solid phases. Vsolid(t) is the summation of all volumes of the solid phases determined by Eq. (21). It can be seen that the error increases at higher degree of hydration. This is caused by the increasing complexity in the interaction between particles during hydration. However, the error is still small at a high degree of hydration (smaller than 5  105), which demonstrates the efficiency of the proposed procedure for calculating expansion.

Fig. 8. Porosity curves of RHA-blended PC samples for prolonged hydration.

Fig. 9. CH fraction curves of the RHA-blended PC samples for prolonged hydration.

pared in Fig. 8, while a comparison in the CH content is shown in Fig. 9. There are several reasons why the simulated values are different from the experimental data. The number of samples (Section 4.1) is too small to obtain reliable estimates for the functions of Ko-values. Moreover, the experimental data are resulting from an indirect method, which leads to experimental biases. The proposed numerical estimation of outward particle expansion is validated indirectly on the basis of the error in numerical data on porosity as compared

5. Discussion and conclusions In this article, the IPKM for simulating the hydration of C3S has been extended to a multi-component integrated particle kinetic model (XIPKM) for simulating the hydration of cement and pozzolanic admixtures. The unhydrated core of a model particle consists of four main components, i.e., C3S, C2S, C3A and C4AF for cement or two components for the pozzolanic admixture, i.e., SiO2 and inert material. The microstructure a model paste is visualized by a digital-based ‘‘voxel system’’. By forming the constituent components of the unhydrated part of a particle by different spherical sectors, the hydration rate of a particle is controlled by the different penetration rates of the phases. As a result, the presented model is able to account for the changes in the constituent proportions of a particle during hydration. Contrary, in many multi-component hydration models, the constituent proportions are assumed unchanged during hydration, while the unhydrated core is described by a homogenized single phase sphere [13,14]. This assumption also means that the hydration rate is the same for all components, which is unrealistic since in practice different components hydrate at different rates. Although based on the multi-component concept involving the four main cement compounds, CCBM [13] assumes the hydration products represented by a homogenized hydration product that adheres to the cement particles. Further, the extended version of HYMOSTRUC3D [14] assumes the hydration product to include only CSH and CH. In the proposed model, four distinct hydration products of the four major cement compounds are involved, i.e., CSH, CH, CAH and FH. CSH precipitates on the cement particles, CH and CAH diffuses and crystallizes randomly in the pore space and FH diffuses and crystallizes near the dissolution source. This

N.L.B. Le et al. / Computational Materials Science 78 (2013) 12–21

is in agreement with the assumptions by Bentz [4] based on the mobility of diffusing ions. By subdividing the remaining contents of water and CH (for pozzolanic reaction) according to the number of particles and then according to the number of components, the simulated hydration process can proceed concurrently, thereby avoiding the orderbased procedure that is popular in vector-based approaches. However, only the hydration reactions proceed concurrently; the deposition of the hydration products and the resulting expansion evolution of the articles are still governed by the order-based procedure. This concurrent process makes the hydration simulation of the individual particles as well as the formation of the hydrate structure more realistic, of course. In the hydration simulation systems based on the vector approach, the estimation of outward growth of the hydrating particles plays a crucial role in the simulation of microstructural evolution. Indeed, since the particles expand by the precipitation of the hydration products and thus gradually occupy part of the pore space, the precision in estimating the particle expansion obviously affects the reliability of measures for geometry and topology (pore continuity!) of pore space. Estimating the size increment of a spherical particle caused by a new volume of precipitated hydration products is a complicated task due to the growing complexity of interferences among particles. In the presented model, the particle expansion is estimated by a numerical procedure that accounts for such implicit interferences. In many vector approach models, the basic rate constant of a phase is commonly determined through a calibration process [7,13,14] in which many trial values of the basic rate constants are attempted until the hydration curves (DOH vs. time curves) fit well with experimental results. Instead of a laborious calibration process, herein, a numerical procedure is proposed to obtain the basic penetration rates of different minerals from the hydration data of the experiments.

21

References [1] D.D. Bui, Rice husk ash as a mineral admixture for high performance concrete, PhD Thesis, Delft University of Technology, DUP Science, Delft, 2001. [2] S. Bishnoi, K.L. Scrivener, Cement Concrete Res. 39 (2009) 266–274. [3] D.P. Bentz, A Three-Dimensional Cement Hydration and Microstructure Program. I. Hydration Rate, Heat of Hydration, and Chemical Shrinkage, NISTIR 5756 U.S. Department of Commerce, 1995. [4] D.P. Bentz, J. Am. Ceram. Soc. 80 (1997) 3–21. [5] D.P. Bentz, P.V. Coveney, E.J. Garboczi, M.F. Kleyn, P.E. Stutzman, Model. Simul. Mater. Sci. Eng. 2 (1994) 738–808. [6] H.M. Jennings, S.K. Johnson, J. Am. Ceram. Soc. 69 (1986) 790–795. [7] K. van Breugel, Cement Concrete Res. 25 (1995) 319–331. [8] P. Navi, C. Pignat, Adv. Cement Based Mater. 4 (1996) 58–67. [9] P. Navi, C. Pignat, Comput. Mater. Sci. 16 (1999) 285–293. [10] P. Navi, C. Pignat, Cement Concrete Res. 29 (1999) 507–514. [11] C. Pignat, P. Navi, K. Scrivener, Mater. Struct. 38 (2005) 459–466. [12] P. Stroeven, M. Stroeven, Cement Concrete Compos. 23 (2001) 189–200. [13] I. Maruyama, T. Matsushita, T. Noguchi, Numerical modeling of Portland cement hydration based on particle kinetic model and multi-component concept, in: Proceedings of the 12th International Congress of the Chemistry of Cement (ICCC), Montréal, Canada, 2007. [14] V.T. Nguyen, Rice Husk Ash as a Mineral Admixture for Ultra High Performance Concrete, PhD Thesis, Delft University of Technology, Ipskamp Drukkers, Delft, 2011. [15] D.P. Bentz, O.M. Jensen, A.M. Coats, F.P. Glasser, Cement Concrete Res. 30 (2000) 953–962. [16] H.M. Jennings, L.J. Parrott, J. Mater. Sci. 21 (1986) 4053–4059. [17] E.A. Rakhmanov, E.B. Saff, Y.M. Zhou, Math. Res. Lett. 1 (1994) 647–662. [18] E.B. Saff, A.B.J. Kuijlaars, Math. Intell. 19 (1997) 5–11. [19] M.D. Cohen, R.D. Cohen, J. Mater. Sci. 23 (1988) 3816–3820. [20] P. Stroeven, H. He, M. Stroeven, J. Zhejiang Univ. Sci. A 12 (2011) 335–344. [21] P. Stroeven, J. Hu, M. Stroeven, Comput. Concr. 6 (2009) 133–153. [22] P. Stroeven, L.B.N. Le, L.J. Sluys, H. He, Image Anal. Stereol. 31 (2012) 55–63. [23] S.A.S. El-Hemaly, R. El-Sheikh, F.H. Mosalamy, H. El-Didamony, Thermochim. Acta 78 (1984) 219–225. [24] J.I. Escalante-García, J.H. Sharp, Cement Concrete Res. 28 (1998) 1245–1257. [25] I. Odler, J. Schüppstuhl, Cement Concrete Res. 11 (1981) 765–774. [26] G. Frigione, S. Marra, Cement Concrete Res. 6 (1976) 113–127. [27] U. Danielson, Heat of hydration of cement as affected by water-cement ratio, in: The 4th International Symposium on the Chemistry of Cement, Washington, 1962, pp. 519–526.