A dissolution model of alite coupling surface topography and ions transport under different hydrodynamics conditions at microscale

Portland cement is the most produced material in the world. The hydration process of cement consists of a group of complex chemical reactions. In order to investigate the mechanism of cement hydration, it is vital to study the hydration of each phase separately. An integrated model is proposed in this paper to simulate the dissolution of alite under different hydrodynamic conditions at microscale, coupling Kinetic Monte Carlo model (KMC), Lattice Boltzmann method (LBM) and diffusion boundary layer (DBL). The dissolution of alite is initialised with KMC. Two Multiple-relaxation-time (MRT) LB models are used to simulate the fluid flow and transport of ions, respectively. For solid-liquid interface, DBL is adapted to calculate the concentration gradient and dissolution flux. The model is validated with experiment from literature. The simulation results show good agreements with the results published in the literature.


Introduction
Portland cement is the most produced material in the world.The hydration process of cement consists of a group of complex chemical reactions.In order to investigate the mechanism of cement hydration, it is vital to study the hydration of each phase separately.
Portland cement contains four main components: alite, belite, calcium aluminate and calcium aluminoferrite.Alite (tricalcium silicate) is the major and characteristic mineral phase.The thermodynamics and kinetics of alite hydration still attracts many attentions [9][10][11][12][13][14], after it has being studied for more than 100 years.Some models has been proposed in order to predict the dissolution of alite.μic [7] was used by Kumar et al. [15] to study the dissolution and precipitation mechanisms of alite.μic takes the vector approach, similar to [1,3].μic uses arbitrarily chosen hydration kinetics to calculate the hydration rate and converts the amount of hydration product into volumetric term [7], thermodynamics and ions transport are also not taken into accounts.
Nicoleau and Bertolim [16] proposed an analytical model to simulate the dissolution of alite.This model described the etch pit formation on the surface of alite, and calculated the corresponding dissolution rate.However, diffusion and chemistry was not considered in this model.
In previous study [17], Martin et al. proposed a Kinetic Monte Carlo model (KMC) at nanoscale, which could simulate the dissolution of crystal not only in dilute solutions but also at close to saturated conditions.It essentially simulated all dissolution mechanisms at nano-microscales.An integrated model is proposed in this paper to simulate the dissolution of alite under different hydrodynamic conditions at microscale, coupling KMC [17], Lattice Boltzmann method (LBM) [18] and diffusion boundary layer (DBL) [19].
The dissolution of alite is initialised with KMC.Two Multiplerelaxation-time (MRT) LB models are used to simulate the fluid flow and transport of ions, respectively.For solid-liquid interface, diffusion boundary layer is adapted to calculate the concentration gradient and dissolution flux.The model is validated with experiment from literature [11].

Methodology
The dissolution of alite is initialised by KMC.Each voxel on the surface of alite is simulated individually.The ion concentration of Ca 2+ at the surface of each voxel is converted to Gibbs free energy ΔG with the function of ions concentration [20].KMC takes ΔG as input, and the dissolved volume is given as output.The dissolution rate at nanoscale is calculated with time interval, and it is upscaled to microscale directly.The dissolved volume of each voxel on the surface of alite is calculated with upscaled dissolution rate, followed by the calculation of dissolved thicknesses Δz.The overall dissolution rate of alite is calculated with the mean Δz [11].

Kinetic Monte Carlo model
Kinetic Monte Carlo (KMC) methods use statistics to describe the time evolution of a system as a sequence of possible events with a certain rate to happen [21].Their main advantage is that allow to simulate time scales comparable to experiments with atomic scale resolution of the process [22][23][24][25][26].In order to carry out KMC study the dissolution of any mineral, we have developed a code called KIMERA [27], which has been used in this work to investigate the specific case of alite.
The dissolution process can be described as sequence of events taking place with different rates.Far from equilibrium conditions the rates follow a trend in agreement with the Arrhenius equation [28].
where f f is the fundamental reaction frequency k b T ℏ , equal to 4 ⋅ 10 13 s − 1 (at room T) [28][29][30], E D is the dissolution reaction energy barrier (kJ⋅mol − 1 ), k B is the Boltzmann constant (1.38 ⋅ 10 − 23 J⋅K − 1 ) and T the temperature (K).Close to equilibrium the saturation of the solution in contact with the dissolving solid starts to play an important role.In order to get the dissolution rate at close to equilibrium conditions it is necessary to couple the precipitation term [17].
where ΔG* is the local Gibbs free energy at the water mineral interface, and E P is the precipitation energy barrier.E D and E P are parameters which depends specifically on each mineral, solvent, pH, structure and topography and are possible to be calculated by ab initio calculations [23,28,[31][32][33].In this model, presented in detail in [17], Eq. (1) describes the dissolution reaction rate, and Eq. ( 2) describes the precipitation reaction rate.Note that the precipitation reaction does not involve C-S-H nucleation and growth, but a reformation reaction bonds in the alite surface.Using both equations together, we take into account the microscopic reversibility of the chemical reactions.As physical interpretation, it can be understood as if the dissolution process is intrinsic of the material, and therefore independent of the pore solution chemistry, yet the precipitation does depend on the ionic concentration of the pore solution: the higher the concentration of ions in solution, the higher the precipitation rate.Therefore, the Gibbs free energy difference is included only in the precipitation reaction, reducing the exponential term in the rate and consequently increasing the rate.Far from equilibrium the ΔG* is small and it does not affect the rate.As a model system to reproduce Alite, we used a Kossel crystal, or Terrace Ledge Kink system (TLK).The TLK is a simple mineral model consisting of a cubic structure with six first neighbors [30] and different topological features on the surface (adatoms, terraces, vacancies, dislocations) (shown in Fig. 1).Despite the simplicity of this system, it ensures enough topographical and structural details so as to reproduce the mechanisms attributed to the dissolution process [17].In this scheme, both E D and E P depend on the number of neighbors assuming a lineal contribution [22].
where E da and E pa is the contribution for an atom with only one bond, i.e. an adatom (shown in Fig. 1), for dissolution and precipitation respectively.
It remains necessary to establish a relationship between the local ΔG* and the macroscopic ΔG.Experimentally, ΔG can be calculated by the ions concentration of the mineral within the water [28,34]: where Π is the ion activity product of the dissolved material and K s its thermodynamic solubility product.The ion activity product Π of alite is calculated by: The macroscopic ΔG is usually different form local ΔG* used in Eq. ( 2), since close to the surface the ionic concentration differs from the bulk one due to the local charge gradient in the electrical double layer [35].In equilibrium conditions, when ΔG=0 and the dissolution and precipitation rates described in Eqs.(1) and 2 should be equal, ΔG and ΔG* can be related by equation [17,27]: where n k is the coordination of a kink site.A kink site, also referred to as the "half-crystal position" is a surface site which has half the number of neighboring atoms as an atom in the crystal bulk, regardless of the type of crystal lattice.Kink sites are located in terraces, with three neighbors in a Kossel crystal.It is a special site because its cohesive energy is also half of that in the bulk.In precipitation and dissolution, the energy of attachment and detachment are therefore the same, and the equilibrium condition is that in which kink sites do not dissolve and do not precipitate.
The dissolution of alite is a very complex process, with multiple coupled reactions in a disordered interface [36], and the reaction energy barriers have not been calculated by atomistic simulations.Therefore, we have determined the values of E da and E pa in our model by fitting the experimental dissolution rate from reference [11].In that work, Juilland and Galluci used a flow-through-cell shown in Fig. 3, controlling the calcium hydroxide concentrations and flow rates of solution, to investigate the alite dissolution rates under different hydrodynamics conditions.First, we calculate the E da fitting the rate at far from equilibrium conditions (ΔG ≥ 25 kcal⋅mol − 1 in Fig. 2) when precipitation does not contribute to the overall dissolution rate and only modeldis is relevant.Once E da is known, we use the ΔG crit value to obtain E pa .The ΔG crit is the free energy value when the dissolution rate suddenly increases, as can be seen in Fig. 2. That increase in ΔG crit is associated to the undersaturation concentration at which dislocations starts to open.For a Kossel crystal, the change in the regime implies that sites at dislocations, with coordination 4, cannot dissolve above ΔG crit , and start to dissolve below the ΔG crit .In a previous work [17] we showed that it was possible to obtain an approximately the ΔG crit value with the numerical relationship: where α = − 2.3 ± 0.9 k B T units [17].From the fitting, we obtained E da =7.47 and E pa =4.95 k B T units.
Once the model parameters for alite were determined, they are used to investigate the dissolution rate of a kossel crystal representing the crystal with the different topographical features (terraces, adatoms, vacancies and dislocations).A value of 10 14 m − 2 dislocation density was considered, which lies in the typical values for most of the minerals (10 10 -10 14 m − 2 ) [37].A random number with random of terraces, adatoms and vacancies are constructed initially to represent the variability of the grain surfaces.The limit cases, thought very unlikely, are on the one hand a perfect surface, and on the other hand a perfect terrace.The system size changes with the voxel size later considered.It varies between 25 × 25 and 150 × 150 Å.The simulated time lies in the order of seconds.The results of the dissolution rate with ΔG, showed in Fig. 2 feed the Lattice Boltzmann simulations, as explained below.

Lattice Boltzmann method
The dissolution of alite under hydrodynamics condition is a classic advection-diffusion problem.The advection is the transport of heat, mass, and momentum by bulk motion of fluid particles.The advection-diffusion process is a process where both advection and diffusion take place simultaneously [18].Lattice Boltzmann Method (LBM) is widely used to solve advection-diffusion equation [38][39][40][41][42].
The discrete velocity space in LBM is classified as DnQm scheme, where Dn and Qm stand for n dimensions and m velocities, respectively.For example, D3Q15 is a 3-dimensional Lattice Boltzmann model with 15 discrete velocities on a cubic grid, shown in Fig. 3 [43].Each line with an arrow represents a velocity in both magnitude and direction.
To overcome the numerical instability of single-relaxation-time model (also called the Bhatnagar-Gross-Krook (BGK) model) [44][45][46], Multiple-relaxation-time (MRT) [47] LB models are used to simulate the fluid flow and transport of ions in this paper.

Fluid flow
A D3Q15-MRT model is used to simulate the fluid flow, given in Eqs. 9 f i is the density distribution of fluid.δ t is the time increment.Ω f is collision term in MRT. e i are the discrete velocities.f i eq is the corresponding equilibrium distribution function.The equilibrium distribution for fluid flow f i eq is given in Eqs.10: where c i is unit vector along the streaming direction, c s is the speed of sound and w i is the associated weight coefficients.For D3Q15, and w i is given by: , ω 1− 6 = 1 9 , ω 7− 14 = The macroscopic fluid density ρ and velocity u are defined by To collision term in MRT for fluid flow is: where N is the transformation matrix, given in Eqs. 15.
and the corresponding diagonal relaxation matrix S f is denoted by where 0<λ<2.

Ions transport
A D3Q7-MRT is used to simulate the transport of ions, instead of D3Q15-MRT, to reduce the burden of computation [18].The density distribution of ions g i is given by Eq. ( 17) The equilibrium distribution for ions transport g i eq is given in: For D3Q7, c s = Δx/(2Δt) and w i is given by The collision matrix in MRT for ions transport is where transform matrix P is defined as: and the corresponding diagonal relaxation matrix S g is: where S 0 could be any value such as 0; S 1 = S D =1/τ relates to diffusivity; S 2 =2 − S D , corresponding to second-order momentum.

Diffusion boundary layer
Diffusion boundary layer (DBL) is a thin liquid layer between bulk solution and the surface of solid [19].When flowing fluid is presented, the transport of ions in the bulk solution occurs through convection, while the transport in the diffusion layers occurs only through diffusion.In this layer, a transverse concentration gradient of ions is formed, leading to ions transport.The exchange of ions between bulk solution and surface of solid is significantly limited by this phenomenon [19,48,49].In this paper, DBL thickness δ DBL (unit: mm) is calculated with Eqs.23 where U is mean velocity of fluid, D is diffusion coefficient of ions in the liquid, parameter a = 1686.1 and b = 0.1 [50].
Based on Fick's first law of diffusion, the dissolution flux J from solid surface to bulk solution is calculated with Eq. ( 24), where C l is the mean concentration of ions in the fluid flow.The dissolution flux J equals to the net dissolution rate in KMC, so concentration of ions on the solid surface C s can be deduced and used in next iteration by KMC.

Simulation procedure
The dissolution experiment performed by Juilland and Gallucci [11] is mimicked to validate the model under same hydrodynamic conditions.The experiment was carried out in the flow-through-cell shown in Fig. 4. The cross-section of the channel was 14 mm 2 , the wetted surface of the alite was 15 mm 2 , and the calcium hydroxide concentrations and flow rates of solution flow were controlled to investigate the dissolution rates under different hydrodynamics conditions.
The concentration of calcium hydroxide inlet solution ranged from 0 (deionised water) to 20 mM (saturated), while the flow rate ranged from 9 ml.min − 1 .mm − 2 to 2600 ml.min − 1 .mm − 2 .The absolute height differences Δz were measured with vertical scanning interferometry, to observe the dissolution front.Consequently, the overall rate of alite dissolution r was calculated with the mean Δz over the whole dissolved area, with Eqs. ( 25) and ( 26) [11].
The cross-section of simulation domain is given in Fig. 5.The solution inlet is injected from left boundary, and flows out from right boundary.The other boundaries of the domain is bounce-back wall.Similar with the experiment, the specimen of alite is located at the bottom, close to the outlet.The dimension of domain is 6 × 1.6 × 1.6mm, of which the resolution is 20μm/voxel.
The flow chart of simulation is given in Fig. 6.The simulation is started by the initialisation of liquid with flow and ions concentration.ΔG is calculated based on the ion concentration of Ca 2+ at the surface of alite and ΔG as a function of ion concentration in [20].In every simulation step, the mass changed due to dissolution is calculated by KMC, and the dissolution rate is calculated with Eqs. ( 25) and (26).When the dissolution rate does not change for 1000 cycles, steady state is assumed and the simulation is ended.If steady state is not reached, LBM is used to update the flow and ions concentration in the fluid while DBL is used to calculated the ions concentration on the surface of alite.The simulation cycle continues until steady state is reached and the dissolution rate is compared with the experiment data in the literature [11].

Results and discussion
The dissolution of alite is simulated under different hydrodynamic conditions and different calcium hydroxide concentrations.In Fig. 7, the simulation in saturated portlandite with flux rate 100 ml.min − 1 .mm − 2 is Fig. 4. Experimental setup of dissolution process under hydrodynamic conditions [11].
given as an example.The accumulated dissolved depth of each voxel on the surface of alite is given in Fig. 7a, and instant dissolution rate is calculated by KMC, given in Fig. 7b.The overall dissolution rate is calculated with the averaged Δz over the whole dissolved area, based on Eqs. ( 25) and ( 26) [11].The average dissolution profile of the total alite surface is given in Fig. 7c.
Notably, the accumulated dissolved depth and dissolution rate are not corresponded because KMC calculates the dissolution rate based on the estimation of the probability for atom to truly leave the mineral taking into account both dissolving and precipitation [17], which is an instant value at any moment.However, the dissolution rate fluctuates around the mean value when ΔG is stable, so the dissolved depth converges as the increase of time.
The dissolutions under different hydrodynamic conditions and different calcium hydroxide concentrations are modelled,and the dissolution rates are compared with the experimental data [11].

Influence of the hydrodynamic conditions on the rates of dissolution
The ions dissolution rate from alite surface influenced by the flow rate of fluid is studied from low to high hydrodynamic conditions.The inputs of the simulation are configured as the same as the experiments in reference [11], i.e., the flow rate ranges between 9 ml.min − 1 .mm − 2 and 1300 ml.min − 1 .mm − 2 , and the concentration of calcium hydroxide ranges between 0 (deionised water) and 20 mM (saturated).The simulated dissolution rate of alite as a function of flow rate with varying initial calcium hydroxide concentration is given in Fig. 8 and compared with experiments.
The simulation result shows good agreements with the experiments in general.The dissolution rate increases rapidly at the beginning as the boost of flow rate, until it reaches to a plateau value.The increase of flow rate shows minor impact on the dissolution rate after 200 ml.min − 1 .mm − 2 .This feature is clearly captured by the model.
Specifically, the simulation matches with the experimental data under low initial calcium hydroxide concentration, e.g., deionised water or 5 mM.For higher initial calcium hydroxide concentration, the simulation results have greater dissolution rates, compared with the experiments.The simulated dissolution rate in saturated portlandite condition has the same higher tendency, compared to experiments as shown in Fig. 9.
In Fig. 9, both simulation and experiments show that the dissolution rate rises with the increase of the flow rate under saturated condition, but a plateau is not observed in the experiments.It was explained that this was due to large error bars at high flow rates [11].The simulation shows higher dissolution rate, compared with the experimental data, but the plateau value is reached within the scope of error bar.It is indicated that the simulation model can be used to predict the dissolution rate of alite under extreme hydrodynamic conditions, i.e., high flow rates.

Influence of calcium hydroxide concentration in the solution
The influence of initial concentration in the solution on the dissolution rate is studied by comparing the simulation and experiment results at low flow rates.The rate of dissolution decreases with the increase of the initial concentration in the solution at any flow rate.
The simulated dissolution rates with varying flow rates as a function of the natural logarithm of the ionic product ln(Π) are given in Fig. 10 and compared with experiments.Both simulation and experiment plots have the same S-shape [20], which is typical for pure dissolution of most minerals.
The simulated dissolution rate rises with the increase of flow rates, but the influence of flow rate decreases dramatically when ln(Π) approaches to equilibrium.This feature is also clearly captured by the simulation.
Some results with different initial concentrations are selected to make the matricial evolution, shown in Fig. 11.The slope derived from each set of point in Fig. 11 is a good indicator of the dissolution rate dependency on the flow rate.
In general, the simulation shows higher dissolution rates under same fluid flow rates compared with the experiments.One of the reasons may be attributed to the calculation of diffusion boundary layer.In this paper a simple linear model Eqs.24 is used, while more complex numerical models [51][52][53] may contribute to more accurate calculation of thickness.
In this paper, the dissolution process of alite is simulated under artificial conditions (same as the experiment in [11]), which are quite different than the hydration of cement paste in practical application.Under this artificial condition, hydration products are prevented theoretically, so the dissolution mechanisms and kinetics of alite can be studied.However, dissolution of alite is influenced in twofold in practical application: first, dissolved ions will accumulate in fluid due to fluid with zero flow rate, thus the dissolution rate is decreased; second, precipitations of hydration product consumes the ions, leading to increase of dissolution rate, but they also hinder the transport of ions from the surface of alite to the liquid.The model presented in this paper can simulate the dissolution and transport of ions in practical application, but the precipitation of hydration product still requires further implementation.

Conclusion
A integrated dissolution model of alite is proposed to deal with different hydrodynamics conditions, coupling surface topography and ions transport.In this model, the dissolution of alite is simulated with a Kinetic Monte Carlo model [17], which not only deals with dilute solutions but also close to saturated conditions.Two Multiple-relaxation-   The influence of the hydrodynamic conditions on the rates of dissolution is studied from low to high hydrodynamic conditions.The simulation result shows good agreements with the reference in general.At higher initial concentration, the simulation shows greater dissolution rate, compared with the experimental data.The plateau value is reached at saturated condition, which indicates that the simulation model can be used to predict the dissolution rate of alite under extreme hydrodynamic conditions, i.e., high flow rates.
The influence of initial concentration in the solution on the dissolution of rate is studied by comparing the simulation and experiment results at low flow rates.The simulation plots have the same S-shape with the experimental data, even though they show higher dependency on the flow rate.The cause may be attributed to the choice of model to calculate the thickness of diffusion boundary layer.
In this study, the dissolution of alite is simulated explicitly.This provides a good start point for simulation of other cement components, i. e., belite, calcium aluminate and calcium aluminoferrite.For a complete study on the dissolution of cement particle, it is necessary to involve all of them individually and simultaneously.Additionally, further studies on the diffusion boundary layer are suggested in order to improve the accuracy of the dissolution model proposed in this study.

Fig. 1 .
Fig. 1.Kossel crystal system.Its characteristic topological features on the surface are labeled.

Fig. 2 .Fig. 3 .
Fig. 2. Alite dissolution rate.In brown, experimental dissolution rates versus ΔG in pure water at 300 K from [11].In black, simulations results using the model parameters E da =7.47 and E pa =4.95 k B T units.The dashed line represents the ΔG crit value used for the fitting, 10.7 kcal⋅mol − 1 .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .Fig. 6 .
Fig. 5. Cross-section of simulation domain: mimicking the experiment setup in [11] with solution inlet, solution outlet and specimen (in the red box).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .
Fig. 8. Influences of the flow rate on the dissolution rate of alite for solutions having different initial calcium hydroxide concentrations.

Fig. 9 .
Fig. 9. Influences of the flow rate on the dissolution rate of alite in saturated portlandite conditions.

Fig. 10 .
Fig. 10.Rate of dissolution of alite as function of ln(Π) for three different flow rates.

Fig. 11 .
Fig. 11.Matricial evolution of the alite dissolution rate with the flow rate and the calcium hydroxide concentrations of the wetting solution.