3830
J. Phys. Chem. C 2010, 114, 3830–3836
Master Equation Approach to Gypsum Needle Crystallization G. Dumazer,† A. Smith,‡ and A. Lemarchand*,† CNRS, UniVersite´ Pierre et Marie Curie - Paris 6, Laboratoire de Physique The´orique de la Matie`re Condense´e, UMR 7600, 4 Place Jussieu, Case Courrier 121, 75252 Paris Cedex 05, France, Ecole Nationale Supe´rieure de Ce´ramique Industrielle, Groupe d’Etude des Mate´riaux He´te´rogènes, E. A. 3178, 47-73, AVenue Albert Thomas, 87065 Limoges, France ReceiVed: NoVember 12, 2009; ReVised Manuscript ReceiVed: January 22, 2010
We write and solve a master equation which governs the stochastic growth of interacting gypsum needles. Needle blocking by steric hindrance is reproduced in an effective way. The comparison of the results with submicrometric simulations enables us to choose a small number of fitting parameters to optimize the master equation approach. The stochastic description qualitatively predicts the main features of the dynamics, including the transition toward incomplete dissolution of plaster for increasing values of the number of nuclei per plaster grain. The master equation quantitatively reproduces some morphological properties of the final material, such as mean value and standard deviation of the probability distribution of needle length. 1. Introduction Gypsum is widely used as a building material for its insulating and fireproofing properties, as molds in the industry of ceramics for its porosity, or as high strength material for dental care. The variety of applications gave rise to many experiments which describe the effect of a given additive on the resulting material and its crystallographic structure.1–13 The main existing kinetics models rely on differential equations for macroscopic variables such as the concentrations of some chemical species.14–18 The results of these models are fitted to specific experimental data, so that the models are often unable to predict the effect of a new additive. Only few studies were devoted to the design of a theoretical, predictive tool, able to make the link between the dynamics of gypsum crystallization and the properties of the final material. We recently developed a two-dimensional submicrometric-scale simulation model of gypsum crystallization in order to capture the main steps of the mechanism and to extract a minimal number of control parameters.19 We chose stochastic rules for the growth without overlap of gysum needles, whereas the dissolution of plaster grains was supposed to follow adiabatically the crystallization process. The simulation model has been validated by the comparison of the results with experiments, according to which an appropriate treatment of the plaster grains, aimed at increasing the number of nucleation sites, could lead to smaller gypsum needles with a faster kinetics. We believe that such a description of gypsum crystallization on an intermediate scale, which omits the details of the atomistic scale but includes the fluctuations, is well suited to the optimization of the material properties. Nevertheless, obtaining simulation results is time-consuming and our goal in this paper is to build a more tractable analytical model on the same submicrometric scale. As far as the simulation results compared well with experiments,19 they are regarded here as reference properties which we would like to find thanks to a faster procedure and for the same given parameter values. If successful, we will possess an efficient analytical tool, allowing us to predict * Corresponding author. Tel: +33 (0)1 44 27 72 90. Fax: +33 (0)1 44 27 72 87. E-mail:
[email protected]. † CNRS, Universite´ Pierre et Marie Curie. ‡ Ecole Nationale Supe´rieure de Ce´ramique Industrielle.
material properties in the entire range of parameter values. We propose a stochastic description by a master equation which governs the evolution of the state of interacting needles of gypsum, each needle being described by its length and its blocked or free nature. To build a minimal model and obtain a master equation which accounts for the simulation results as simply as possible, we start from an oversimplified, deterministic model and complexify it progressively. In addition, this procedure has the advantage to offer a better interpretation of the simulation results. The paper is organized as follows. We first recall the simulation model viewed as a benchmark and build step by step an analytical, stochastic description of gypsum crystallization aiming at reproducing the main simulation results. We begin with a simple deterministic analysis based on a differential equation for the length of free gypsum needles that grow without interactions. Then, the model is improved by the description of needle length fluctuations, and we write a master equation which governs the evolution of the probability of the length of a single needle. This simple stochastic approach gives already access to analytical expressions for the mean value and the variance of the reaction time. The comparison with the simulation results allows us to make precise the conditions for which the interaction between needles can be neglected. Finally, we show that the main simulation results about dynamics and morphology can be captured by a multivariate master equation for interacting needles. We prove that it is sufficient to account for steric hindrance in an effective way to recover the statistical results on needle length, without needing the time-consuming generation of spatial configurations of plaster grains and gypsum needles. 2. Submicrometric Simulation Model In the simulation,19 we consider N grains of plaster or calcium sulfate hemihydrate, which are modeled by disks of radius R, that are randomly spread in a square box of side L. We do not describe the formation of the nuclei and G nuclei are randomly placed on the perimeter of each grain at the beginnig of the simulation. They are supposed to already have the shape of embryonic needles of gypsum or calcium sulfate dihydrate: The
10.1021/jp910775w 2010 American Chemical Society Published on Web 02/17/2010
Master Equation for Gypsum Needle Crystallization nuclei are rectangles of length l0 and width l0/R, where R, defined as the ratio of the length and the width, is the aspect ratio of the needles. At any time, the solution is supposed to be supersaturated with respect to dihydrate.15 During the time step ∆t, the growth of each needle is considered and dissolution of hemihydrate is supposed to follow adiabatically. We assume that, during a time step, each needle reaches a state of local equilibrium, so that the growth follows thermodynamic rules and obeys Wulff theorem.20,21 Hence, the aspect ratio R of the needles remains constant. In the initial autocatalytic growth regime, a possible growth is accepted with perimeter-dependent probability pl/lc, where lc is a critical length, from which growth is controlled by the diffusion of the reactants in the liquid phase and occurs with constant probability p. Nevertheless, before accepting the actual growth of a needle, it is necessary to check if its sides are not blocked by steric hindrance, i.e., by the contact with other solids, such as hemihydrate grains and other needles. A needle with four blocked sides does not grow. The growth of a needle with two opposite blocked sides or three blocked sides, which would modify the aspect ratio, is not allowed either. Provided that its growth has been accepted with probability pl/ lc (if l < lc) or p (if l g lc), a needle actually increases its length by ∆l ) 2∆l0 if all four sides are free and by only ∆l ) ∆l0 if one side or two perpendicular sides are blocked (∆l0 is a constant). Simultaneously, the width is increased by ∆l/R, so that the aspect ratio of the needle is preserved. On the simulation time scale, the growth process is irreversible: The needles never shrink. The total volume which precipitated during ∆t determines the width ∆R(t) of the external crown which is lost by each hemihydrate grain. After a succession of growth and dissolution steps, the simulation ends either when the total amount of hemihydrate has dissolved or when all the needles are blocked by steric hindrance. The simulations are performed at a fixed, low density which corresponds to a volume ratio of hemihydrate to water of H/W ) 0.15. The parameters are set to the following values: diffusion-controlled growth probability p ) 0.8, needle aspect ratio R ) 20, length scale set to ∆x ) 1 and time step set to ∆t ) 1 for the sake of simplicity, hemihydrate disk radius R ) 100, initial length l0 ) 3, critical length lc ) 76.2, and ∆l0 ) 2. Note that we derived the order of magnitude of the space and time steps thanks to the comparison of the simulation results with experiments in the same dilute conditions.19 By fitting the time evolution of hemihydrate percentage during a simulation to experimental conductivity curves during plaster setting, we found ∆t = 2 s. The computation of the mean needle length in final configurations obtained by simulation and the analysis of optical microscopy images of gypsum led to ∆x = 50 nm, which justifies the use of submicrometric for the scale of the description. In the simulations, the speed at which a needle increases obeys 0 e ∆l/∆t e 4. The values of ∆t and ∆x define the temporal and spatial resolutions of the description under which the model is not supposed to be valid. A series of simulation has been performed for different values of the number G of nuclei per grain. For each set of parameter values, we have determined averages over 10 realizations of the initial configuration which counts each 600 nuclei. In the following, the statistics is given for 6000 needles whatever the value of G. Such averaged quantities q are denoted qj. 3. Deterministic Approach We first oversimplify the simulation model in order to isolate the characteristics of the dynamics when both the length fluctuations and the interactions between needles are ignored.
J. Phys. Chem. C, Vol. 114, No. 9, 2010 3831 In this section, we derive a differential equation for the length of a single gypsum needle. We consider as initial condition a single hemihydrate disk of radius R with G “nuclei” of length l0. Then, we admit that G needles grow without interaction. If we neglect the initial surface l02/R of the nuclei and the difference of density between hemihydrate and dihydrate, conservation of matter imposes
G
l1002 ) πR2 R
(1)
where l100 is the final needle length that is reached when all of the available plaster has been consumed. We consider the mean growth dynamics of a single needle or, equivalently, several needles without interaction. In the frame of such a simple, deterministic approach, the evolution is governed by a differential equation for the needle length l
dl p∆l ) l, dt lc∆t dl p∆l ) , dt ∆t
l e lc
(2)
lc < t < l100
(3)
where ∆t is the time step and ∆l the length increase. The quantity ∆l/∆t appears as the needle growth speed. For l ) l0 at initial time t ) 0, the solution of eqs 2 and 3 is
( )
l(t) ) l0 exp
l(t) ) lc +
p∆l t, lc∆t
p∆l(t - tc) , ∆t
(4)
t e tc
tc < t < t100
(5)
where tc ) (lc∆t/p∆l) ln(lc/l0) is the critical time at which the length reaches the critical value lc, and t100 is the final time at which the length reaches l100. At time t, the percentage of remaining hemihydrate is given by
%HH ) 100 - z,
with z ) 100
Gl(t)2 RπR2
(6)
where the needle length l(t) obeys eqs 4 and 5. The time tz, at which z% of hemihydrate has been dissolved is then
tz ) tc +
∆t p∆l
(
)
zRπ R - lc 100G
(for z such that tz > tc)
(7)
In particular, the expresion of the final time
t100 )
( ()
)
lc ∆t l ln + l100 - lc p∆l c l0
(8)
is recovered by choosing z ) 100 in eq 7. For the parameter values chosen, tc is supposed to be smaller than t100. In the following, for an easier comparison with the simulation results, we focus on time t10, at which 10% of hemihydrate has been
3832
J. Phys. Chem. C, Vol. 114, No. 9, 2010
dissolved, to characterize the induction period and on time t95, at which 95% of hemihydrate has been dissolved, to evaluate the duration of the reaction. The low density at which the simulations have been performed and the choice of a single nucleation site G ) 1 are favorable conditions for the comparison with the deterministic approach, which ignores needle interactions and steric hindrance. In the deterministic description, the parameters p, R, R, lc, and ∆x are set to the same values as in the simulations. The initial length l0 and the needle growth speed ∆l/∆t are used as fitting parameters, so that the deterministic predictions for the times t10 and t95 that are deduced from eq 7 with G ) 1, coincide with the corresponding simulation results, jt10/∆t ) 200 and jt95/ ∆t ) 537. We find l0 ) 12.849 in place of l0 ) 3 in the simulation and ∆l/∆t ) 1.936. In the simulations, the intantaneous needle growth speed varies in the range 0 e ∆l/∆t e 4, depending on needle environment. The growth of one of the tips of the needle is slowed down by the presence of the grain. The two opposite sides of the needle are rarely allowed to grow simultaneously and the growth speed rarely reaches 4. For a small number of nuclei per grain, here G ) 1, and a low density, the growth speed is larger at the beginning of the simulation, since the needles are far the one from the others, the typical distance between two needles being of the order of the distance between two plaster grains. Then, the needle growth speed decreases due to the strengthening of steric hindrance between different needles. The deterministic approach which ignores these geometrical constraints is only able to produce a meanfield picture through the values of the two parameters l0 and ∆l/∆t. The fitted value of initial length l0 is found larger that in the simulation, in order to account for the larger growth speed at the beginning of the simulation. The fitted growth speed ∆l/ ∆t is found close to the middle of the allowed interval [0, 4]: The mean-field approach gives the picture of a needle with one nearly systematically blocked side and allowed growth on the opposite side. We use eq 6, where l is given by eqs 4 and 5 with l0 ) 12.849 and ∆l/∆t ) 1.936, to draw Figure 1a. The averaged evolution of hemihydrate percentage versus time deduced from the simulations is given in Figure 1c. The simulation model reproduces the typical S-shape of explosive or autocatalytic phenomena with three steps: an induction period with a slow dynamics, a fast regime, and a relaxation toward the final state at which reaction ends. As expected, the deterministic approach fails in mimicing the final slowing down of the reaction due to steric hindrance and the evolution abruptly stops in Figure 1a, simply because the available hemihydrate has been consumed. The evolution of times t10 and t95 with the number G of nucleation sites is deduced from eq 7. According to Figures 2 and 3, the deterministic approach predicts that times t10 and t95 monotonically decrease as the number G of nucleation sites increases. By construction, the results of the deterministic approach are identical to the simulation results for G ) 1. As G increases, the discrepancies between the two approaches increase, revealing the effect of steric hindrance in the simulations. The deterministic predictions are always smaller that the simulation results: t10 e jt10, t95 e jt95. As G increases in the simulations, needle originating from nuclei on the same grains may overlap even at short times: steric hindrance appears from the very beginning of the reaction and increases the duration of the induction period jt10 with respect to the mean-field description. The largest disagreement between the two approaches is obtained for t95, for which the simulations reveal a nonmonotonic behavior. At large G (G g 10), for which needle
Dumazer et al.
Figure 1. Time evolution of hemihydrate percentage %HH for different numbers G of nuclei per grain, G ) 1 (solid line), 10 (dotted line), 30 (a,b), or 25 (c) (dashed line). (a) Predictions of the deterministic approach, (b) Results of master equation integration by Gillespie method. (c) Simulation results. In cases (b) and (c), the bar over %HH refers to averages over 10 realizations of the initial condition and dissolution is not complete for the dashed line. See text for the parameter values.
Figure 2. Averaged time jt10 at which 10% of hemihydrate has been dissolved versus number G of nuclei per grain. The solid line gives the prediction of the deterministic approach, squares correspond to master equation results and triangles to simulation results. The solid symbols are associated with uncomplete dissolution for at least one of the 10 realizations.
entanglement becomes limiting, jt95 increases as G increases. Clearly, the mean-field approach is unable to predict such a refined effect.
Master Equation for Gypsum Needle Crystallization
J. Phys. Chem. C, Vol. 114, No. 9, 2010 3833 evolution of the conditional probability P(l′, t′|l, t) that the needle length is l′ at time t′ knowing that it was equal to l at time t
∂tP(l', t'|l, t) ) γ(l)[P(l', t'|l, t) - P(l', t'|l + ∆l, t)]
(9) with the following transition rate:
γ(l) )
pl p H(l - l) + H(l - lc) lc∆t c ∆t
(10)
where H is the Heaviside function. The probability that the length reaches the boundary l100 at a time larger than t′, knowing l100 P(l′, t′|l, t). that she starts from l at time t, is G(l, t′ - t) ) Σl′)l 0 Consequently, the mean time 〈T(l f l100)〉 necessary for the needle to reach the final length l100 knowing that she starts from length l is given by22,23 Figure 3. Same caption as Figure 2 for the averaged time jt95 at which 95% of hemihydrate has been dissolved.
〈T(l f l100)〉 )
∫0∞ t'∂t'G(l, t') dt') - ∫0∞ G(l, t') dt' (11)
Similarly, the second moment of time T obeys
〈T2(l f l100)〉 ) -
∫0∞ 2t'G(l, t') dt'
(12)
Summing eq 9 over l′ and integrating with respect to t′, we find
〈T(l f l100)〉 - 〈T(l + ∆l f l100)〉 )
1 γ(l)
(13)
A straighforward recurrence leads to 〈T(l0 f l100)〉 ) l100-∆l [1/γ(l)]. Using the expression of γ(l) and the first-order ∑l)l 0 approximation of the Euler-Maclaurin formula24 between sums and integrals, we find Figure 4. Same caption as Figure 2 for the averaged mean needle length 〈lj〉. The solid line gives the prediction of the deterministic approach for the final length l100.
According to eq 1, conservation of matter leads to the expression of the final needle length, l100. The variation of l100 versus the number of nuclei per grain G is given in Figure 4. Contrary to the times t10 and t95, no fit is possible for the final length, since its value is entirely imposed by the initial amount of the available hemihydrate. The value of l100 is independent of the dynamics. As seen in Figure 4, the deterministic prediction of the final needle length is always larger than the mean length deduced from the simulations: l100 g 〈lj〉. The simulations reveal another consequence of steric hindrance: in average, the simulated needles are shorter than in the mean-field description. 4. Stochastic Approach for a Single Needle We now go beyond the deterministic approach and include the description of the length fluctuations, but still ignore the interactions between needles. The length of a needle is seen as a discrete random variable associated with jumps ∆l. We consider the backward master equation,22,23 which governs the
〈T(l0 f l100)〉 =
( ()
lc l100 - lc lc ∆t lc + + ln p ∆l l0 2l0 ∆l
)
(14)
where we assume that lc . l0. The result of the stochastic approach for the mean reaction time 〈T(l0 f l100)〉 introduces the correction (∆t/p)(lc/2l0) with respect to the result t100 of the deterministic approach which is given in eq 8. We follow an analogous procedure to evaluate the second moment of time T. Multiplying eq 9 by (t′ - t), summing over l′ and integrating with respect to t′, we obtain the following recurrence relation: 〈T2(l f l100)〉 - 〈T2(l + ∆l f l100)〉 ) [2〈T(lfl100)〉]/γ(l) or equivalently 〈T2(l0 f l100)〉 ) l100-∆l Σl)l [2〈T(l f l100)〉/γ(l)]. To the leading order, the variance 0 2 σ (T) reads
σ2(T) ) 〈T2(l0 f l100)〉 - 〈T(l0 f l100)〉2
( ) (
2 ∆t 2 lc 1 1 = + p l0 4l0 ∆l
)
(15)
3834
J. Phys. Chem. C, Vol. 114, No. 9, 2010
Dumazer et al. 5. Stochastic Approach for Interacting Needles We now propose a stochastic description for many interacting needles and write a master equation for the probability of observing n interacting needles in a given state at time t. The state of needle i, with 1 e i e n, is defined by its length li and its free or blocked nature bi, where bi ) 0 if needle growth is allowed and bi ) 1 if the needle is blocked. For N hemihydrate disks of radius R, the maximum length that a needle can reached is ltot2 ) RNπR2. The initial state corresponds to n free embryonic needles, such that li ) l0, bi ) 0 for 1 e i e n. The master equation reads n
∂tP(l1, b1, ..., li, bi, ..., ln, bn) )
∑ {γ(li - ∆l)P(li -
i)1 n
∆l, bi)δ(bi)H(ltot2 Figure 5. Simulation results for the time evolution of hemihydrate percentage %HH for a single nucleation site per grain G ) 1. The two dashed lines give two individual realizations for different initial configurations of hemihydrate disks and nucleation site positions. The solid line gives the averaged result %HH over 10 realizations of the initial condition.
Note that the variance is independent of l100: The fluctuations mainly develop during the autocatalytic regime at the beginning of the growth. The scaled standard deviation is given by
σ(T)/〈T(l0 f l100)〉 =
(
()
ln
lc l0
)
∆l ∆l +1 l0 4l0 l100 - lc ∆l + + 2l0 lc
(16)
Using eq 1 for a single nucleation site per hemihydrate grain (G ) 1), we obtain l100 ) 792.7. For the values l0 ) 3, ∆l ) 2, and lc ) 76.2, we find that the scaled standard deviation reaches 6.8%. This prediction agrees well with the low-density simulation results obtained for G ) 1. Figure 5 shows two representative realizations of hemihydrate percentage evolution for two different initial configurations of plaster grains and nucleation sites: large fluctuations of the reaction time are observed. Note the abrupt end of the reaction for the individual realizations in dashed lines: the reaction suddenly stops because a solid phase (hemihydrate) disappears and not because chemical equilibrium has been reached. Whatever the size of the simulated system, we have checked that nonvanishing slopes are observed as the percentage of hemihydrate %HH tends to zero. The smoothing of the leading edge of the solid line curve %HH(t) in Figure 5 is a consequence of the averaging procedure over different realizations. At the end of the reaction, the effect of the fluctuations is nontrivial and the shape of the averaged curve %HH(t) differs from the shape of the individual realizations. The simulation results lead to a scaled standard deviation of time t95 of the order of 7%. This good agreement between the stochastic approach for a single needle and the simulations enables us to conclude that, for a low density (small volume ratio of hemihydrate to water H/W ) 0.15) and for a single nucleation site per grain (G ) 1), the interaction between needles can be neglected.
∑
lj2 - (li - ∆l)2) +
j)1,j*i
[-γ(li)Pδ(bi) + k(li)P(li, bi - 1)δ(bi - 1) n
k(li)Pδ(bi)]H(ltot2 -
∑ lj2)}
(17)
j)1
with γ(li) ) (pli/lc∆t)H(lc - li) + (p/∆t)H(li - lc), k(li) ) kG + kdH(li - d), in which H denotes the Heaviside step function. kG and kd are effective rate constants that accounts for needle blocking due to steric hindrance. In the right-hand side, the needle lengths and natures, at which the probability P has to be evaluated, are omitted, when they are identical to those of the left-hand side. The rate constant kG aims at reproducing that a needle can be blocked from the beginning, in particular, due to the presence of other needles on the same disk. The value of kG is supposed to increase as the number G of nucleation sites increases. The rate constant kd characterizes the blocking that arises after the needle has reached a cutoff length d. The cutoff length varies with the mean distance between hemihydrate grains and increases with density. The rate constant kd reproduces the interaction of needles that originate from different discs. The master equation governs the evolution of the probability of observing a given state l1, b1, ..., li, bi, ..., ln, bn for the n needles. The state evolves due to jump processes. To write the master equation, we have identified all the processes that start from the state l1, b1,.. ., li, bi,.. ., ln, bn and those that lead to it.22,23 Such processes can be seen as chemical reactions and the probability that a process occurs is controlled by a chemical rate constant. Let us consider as an example a process in which only the variable bi changes, i.e., a process which modifies the free needle i into a blocked needle, all the other variables remaining unchanged kG
l1, b1...li, bi ) 0, ..., ln, bn f l1, b1, ..., li, bi ) 1, ..., ln, bn (18) where the needles are in the state l1, b1, ..., li, bi ) 0, ..., ln, bn before the process occurs and in the state l1, b1, ..., li, bi ) 1, ..., ln, bn after. This process leads to the term -kGP(l1, b1, ..., li, n bi, ..., ln, bn)δ(bi)H(ltot2 - Σj)1 lj2) in the master equation. This expression is contained in the term -k(li)Pδ(bi)H(ltot2 - Σnj)1lj2) of eq 17 where k(li) ) kG + kdH(li - d). The contribution is negative since the process destroys the state l1, b1, ..., li, bi, ..., ln, bn provided that bi ) 0, which justifies the Dirac distribution δ(bi): The contribution vanishes if the initial state is l1, b1, ...,
Master Equation for Gypsum Needle Crystallization
Figure 6. Effective blocking rate constant kG versus number G of nuclei per grain. Values used in the master equation approach. The solid symbol is associated with uncomplete dissolution for at least one of the 10 realizations.
Figure 7. Same caption as Figure 4 for the averaged standard deviation σ j (lj) of the distribution function of needle length.
li, bi ) 1, ..., ln, bn, which agrees with the hypothesis of irreversible blocking. The process given in eq 18 modifies the probability that the needles are in the state l1, b1, ..., li, bi, ..., ln, bn proportionally to the probability P(l1, b1, ..., li, bi, ...., ln, bn) to be in this state before the process occurs and proportionnally to the rate constant kG. The Heaviside24 step function H(ltot2 n lj2) ensures that the process does not occur when the total Σj)1 mass of available gypsum (proportional to ltot2 in a 2-dimensional model) has already been precipitated. We now explain how we choose the parameter values and begin with the cutoff length d introduced in the expression of the rate constant kd to mimic the interaction of needles which grow on different hemihydrate disks. A rough estimate of the distance between the centers of two discs, for N discs placed in a square of side L, is given by L/N. For the density chosen in the simulations, we have L/N ) 460. We choose a value d of the cutoff length close to L/N - 2R and impose d ) 300. The parameters p, R, lc, R, ∆x and ∆t are fixed at the same values as in the simulations: p ) 0.8, R ) 20, lc ) 76.2, R )
J. Phys. Chem. C, Vol. 114, No. 9, 2010 3835 100, ∆x ) 1, ∆t ) 1. The values of the rate constants kG)1 and kd, the initial needle length l0 and the length increase ∆l are chosen such that 〈lj〉, σ j (l), jt10, and jt95 are the same for the master equation and simulation results for G ) 1. For the rate constants, we find kG)1 ) 6.8 × 10-4 and kd ) 10-4. The nonvanishing value of kG)1 reproduces the fact that even a small needle without neighboring needles can already be blocked because of the hemihydrate grains. For the lengths, we find l0 ) 8 and ∆l ) 2.255. Note that the values of 〈lj〉 and σ j (lj) are mainly determined by the choice of kG)1 and kd, whereas the values of jt10 and jt95 essentially depend on l0 and ∆l. For the sake of simplicity, the different parameters do not vary with G except kG. The values of kG for the considered values of G > 1 are found by matching the results of the master equation for the averaged mean needle length 〈lj〉 with the corresponding simulation results. To solve the master equation, we use the Gillespie algorithm, which generates the stochastic processes according to the correct statistics.25 The time t is discretized and the time step ∆t is sampled from an exponential distribution. The growth or blocking of a randomly chosen needle is accepted according to the instantaneous value of its growing or blocking rate. If the growth process is accepted, the length of the needle is increased by ∆l and the total amount of available hemihydrate is decreased by the part that just precipitated. If the blocking process is accepted, the length of the needle is irreversibly fixed at the value reached at time t. After a process is accepted, the corresponding rate is updated and time is increased by ∆t. The procedure is repeated until hemihydrate has been completely consumed or all the needles are blocked. In order to compare with the simulation results we generate initial conditions with the same number (600) of embryonic needles and perform averages over 10 realizations. As shown in Figures 1b, 1c, the results of the master equation about kinetics qualitatively agree with the corresponding simulation results. In particular, the transition toward the incomplete dissolution of hemihydrate, due to the blocking of all the needles, is recovered in the stochastic approach. As shown in Figures 2 and 3, the agreement between the results of the simulations and the master equation is excellent for the induction period jt10 but is only qualitative for the reaction time jt95. Let us recall that only the values of jt10 and jt95 deduced from the master equation for G ) 1 are fixed to the corresponding simulation results. Only the rate constant kG is adjusted as G varies, in order to ensure the equality of the mean needle lengths deduced from the two descriptions, but not to optimize dynamics. Thus, the prediction of the master equation for the variation of jt10 and jt95 with the number of nuclei G can be considered as a test for the stochastic approach. The effective way of taking steric hindrance in the case of the master equation is able to predict a minimum for the variation of jt95 with G but this minimum is reached for a larger value of G, about 20, whereas the simulations lead to a minimum for G ) 10. The results of the master equation concerning the morphology of the material are good. However, the excellent agreement between the simulations and the master equation for the mean needle length 〈lj〉 observed in Figure 4 is simply inherent to the fitting procedure. As shown in Figure 6, this requirement leads to a linear variation of kG with G in the explored range. According to Figure 7, the good agreement between the standard deviations σ j (l) of the length deduced from the simulations and the master equation constitutes a satisfying test of the master equation and in particular of its manner to account for the interactions between needles.
3836
J. Phys. Chem. C, Vol. 114, No. 9, 2010
Dumazer et al.
6. Conclusion
References and Notes
The multivariate master equation that we propose in this work offers a stochastic description of gypsum needle growth in which needle blocking is reproduced in an effective way. The comparison with submicrometric simulations, that give access to the spatial organization of the needles, enables us to choose a small number of fitting parameters and to optimize the results of the master equation. The master equation approach accounts for the main qualitative properties of dynamics and morphology that were exibited by the simulations: The stochastic dynamics reproduces (i) the autocatalytic behavior with an induction period and a sudden acceleration, (ii) the nonmonotonic variation of the reaction time as the number G of nuclei per plaster grain varies, (iii) the incomplete dissolution of plaster for large values of G. The two last results prove that the linear dependence of the blocking rate constant kG with G, deduced from the fitting procedure, correctly mimics spatial effects such as steric hindrance between needles. Some important morphological properties of gypsum are captured by the stochastic description: The mean value and the standard deviation of the needle length are in good quantitative agreement with the simulations. In particular, the mean needle length is smaller than the rough prediction of the deterministic approach which ignores the interactions between needles. We conclude that the master equation offers a satisfying alternative to more costly simulations. However, the submicrometric simulations remain a powerful tool to investigate more deeply the spatial properties of the resulting material. The analysis of the final needle configurations should enable us to define morphological criteria that could be related to the mechanical properties of the material. Some relevant quantities, that require the resort to simulations, are the mean number of contact per needle to characterize global resistance, the needle radial distribution function around the initial plaster grains to check the homogeneity of the structure, or the pore distribution function. The comparison with experiments should help in choosing appropriate criteria.
(1) Rinaudo, C.; Robert, M. C.; Lefaucheux, F. J. Cryst. Growth 1985, 71, 803–806. (2) Amathieu, L.; Boistelle, R. J. Cryst. Growth 1986, 79, 169–177. (3) Rinaudo, C.; Franchini-Angela, M.; Boistelle, R. J. Cryst. Growth 1988, 89, 257–266. (4) De Vreugd, C. H.; Witkamp, G. J.; van Rosmalen, G. M. J. Cryst. Growth 1994, 144, 70–78. (5) Rinaudo, C.; Lanfranco, A. M.; Boistelle, R. J. Cryst. Growth 1996, 158, 316–321. (6) Bosbach, D.; Junta-Rosso, J. L.; Becker, U.; Hochella, M. F., Jr. Geochim. Cosmochim. Acta 1996, 60, 3295–3304. (7) Badens, E.; Veesler, S.; Boistelle, R. J. Cryst. Growth 1999, 198/ 199, 704–709. (8) Boisvert, J.-P.; Domenech, M.; Foissy, A.; Persello, J.; Mutin, J.C. J. Cryst. Growth 2000, 220, 579–591. (9) Brandt, F.; Bosbach, D. J. Cryst. Growth 2001, 233, 837–845. (10) Hill, J.-R.; Plank, J. J. Comput. Chem. 2004, 25, 1438–1448. (11) Rashad, M. M.; Mahmoud, M. H. H.; Ibrahim, I. A.; Abdel-Aal, E. A. J. Cryst. Growth 2004, 267, 372–379. (12) Reynaud, P.; Saadaoui, M.; Meille, S.; Fantozzi, G. Mater. Sci. Eng., A 2006, 442, 500–503. (13) Ersen, A.; Smith, A.; Chotard, T. J. Mater. Sci. 2006, 41, 7210– 7217. (14) Rosmalen, G. M.; Daudey, P. J.; Marche´e, W. G. J. J. Cryst. Growth 1981, 52, 801–811. (15) Amathieu, L.; Boistelle, R. J. Cryst. Growth 1988, 88, 183–192. (16) Beretka, J.; van der Touw, J. W. J. Chem. Tech. Biotechnol. 1989, 44, 19–30. (17) Hand, R. J. Cem. Concr. Res. 1994, 24, 885–895. (18) Hernandez, A.; La Rocca, A.; Power, H.; Graupner, U.; Ziegenbalg, G. J. Cryst. Growth 2006, 295, 217–230. (19) Dumazer, G.; Narayan, V.; Smith, A.; Lemarchand, A. J. Phys. Chem. C 2009, 113, 1189–1195. (20) Venables, J. A. Introduction to Surface and Thin Film Processes; Cambridge University Press: Cambridge, U.K., 2000. (21) Herring, C. Phys. ReV. 1951, 82, 87–93. (22) Van Kampen, N. G. Stochastic processes in physics and chemistry; North Holland: Amsterdam, 1983. (23) Gardiner, C. W. Handbook of stochastic methods; Springer: Berlin, 1985. (24) Abramowitz, M.; Stegun, I. A. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables; Dover: New York, 1972. (25) Gillespie, D. T. J. Phys. Chem. 1977, 81, 2340–2361.
JP910775W
