Cement Hydration Kinetics and Rate

[MaryamZamzami]

Dear PHREEQC users!

I'm trying to run the kinetics for cement hydration based on equations by Parrot and Killoh (1984). I already found the same question on the PhreeqcUsers website but it does not give me the correct graph that I want. So I tried to change the code based on a paper by Holmes et. al. (2022). Right now my problem is that I don't think "TIME" is exactly doing what I want. I need to define "a_t" as "a+R_hyd * delta_t" where delta_t is the difference between the current and previous time step. Could anyone please help me what should I do?
Thank you for your time.

Code: [Select]

DATABASE D:\Software\PHREEQC\GroutDegradation\concrete_3T_V07_02.dat
TITLE Cement Hydration # CEM I 42.5 N - SR 3 MH/LA with a w/c ratio of 0.47

SOLUTION 1
    temp      25
    pH        7    charge
    -water    0.047 # kg
   
EQUILIBRIUM_PHASES 1
    Portlandite 0 0
    Brucite 0 0
    HydrotaliceOH 0 0
    HydrotalciteC 0 0
    Syngenite 0 0
    CAH10 0 0
    Al(OH)3am    0 0
    Fe(OH)3mic 0 0

SOLID_SOLUTIONS 1
CSH-I
-comp SiO2am 0
-comp CSHtob1 0
CSH-II
-comp CSHjen 0
-comp CSHtob2 0
AFt
-comp Ettringite 0
-comp Tricarboaluminate 0
-comp Fe-ettringite 0
AFm-OH-SO4
-comp C4AH13 0
-comp C4FH13 0
-comp C2AH8 0
-comp C2FH8 0
-comp stratlingite 0
-comp Fe-stratlingite 0
-comp Monosulfoaluminate 0
-comp Fe-monosulfate 0
AFm-monocarbo
-comp Monocarboaluminate 0
-comp Fe-monocarbonate 0
AFm-hemicarbo
-comp Hemicarboaluminate 0
-comp Fe-hemicarbonate 0
Hydrogarnet
-comp HydrogarnetOH 0
-comp HydrogarnetFe 0
   
REACTION 1
CaO 0.011412268 
SiO2 0.00349534
Al2O3 0.000343272
Fe2O3 0.000288058
MgO 0.000173697
K2O 6.58174E-05
Na2O 1.1294E-05 
SO3 0.000274794
0.47
RATES
    C3S
-start
 10 REM   PARM(1) = A,Blaine surface area of cement, m^2/kg
 20 REM   PARM(2) = A0,reference surface area of cement, m^2/kg
 30 REM   PARM(3) = E,the apparent activation energy of C3S, J/mol
 40 REM   PARM(4) = R,ideal gas constant
 50 REM   PARM(5) = T0,the reference temperature, K
 60 a = M
 70 r1 = (1.5/0.7)*(1-a)*((-LOG(1-a))^0.3)
 80 r2 = (((1-a)^(2/3))*0.05)/(1-((1-a)^(1/3)))
 90 r3 = (1.1*(1-a)^3.3)
110 IF(a > (1.8*0.5)) THEN f = ((((1.8*0.5)-a)*3.333)+1)^4 ELSE f = 1
130 min = r1
140 IF min > r2 THEN min = r2
150 IF min > r3 THEN min = r3
160 R_hyd = f*min*(PARM(1)/PARM(2))*EXP((PARM(3)/PARM(4))*((1/PARM(5))-(1/TK)))
170 a_t = a + R_hyd * TIME / 86400
180 PUT(a_t, 1)
190 unhydrated_C3S = Get(1) * M0
200 SAVE unhydrated_C3S
300 print r1, r2, r3, f, a, a_t
-end

KINETICS 1
C3S
    -formula  (CaO)3(SiO2)
    -m0       0.2621
    -m        0.2621
    -parms    386 385 41570 8.3144 293.15
    -tol      1e-008
   
-steps      864000 in 100 steps

-step_divide  1
-runge_kutta  3
-bad_step_max 5000
# -cvode true
INCREMENTAL_REACTIONS true

USER_GRAPH 1
    -headings               Time C3S
    -axis_titles            "Time (days)" "C3S (%)"
    -chart_title            "C3S"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X total_time/86400
20 GRAPH_Y KIN("C3S") * 228.32
  -end
    -active                 true

END

[dlparkhurst]

TIME is an internally generated length of the current sub-time step used for RATES definitions.

TOTAL_TIME is the current total elapsed time in seconds.

KIN_TIME is the length of the current time step for KINETICS.

[MaryamZamzami]

Thank you for your answer. I fixed it using IPHREEQC and Python language.
 
I am just curious to know if writing a loop in PHREEQC is possible. So, what I want is to calculate the unhydrated amount of each phase (C3S, C2S, C3A, and C4AF) in the KINETICS block, which the current code can do, and then put the mole of hydrated phases inside the EQUILIBRIUM_PHASES.
 
Generally, is it possible (in PHREEQC) to have a variable as the initial mole in EQUILIBRIUM_PHASES that can be updated at each time-step?

[dlparkhurst]

I think the short answer is no.

Perhaps you could treat the mineral as a KINETIC reactant instead of an EQUILIBRIUM_PHASE? You can keep track of initial moles within the RATES definition through the integrated moles, which updates M, or with GET and PUT statements. It will be slower to have more KINETIC reactants, and there is more chance of problems with more kinetic reactants.

Another option is simply to add the moles to the solution through the unhydrated KINETICS reaction. The equilibrium phases will precipitate the appropriate phases. The initial moles won't change, but the current moles will. I think this would be the best way to ensure mole balance.

It is possible to change the initial moles of an equilibrium phase with EQUILIBRIUM_PHASES_MODIFY. The changes defined by this keyword are done before any other calculations, such as initial solution, reaction, or TRANSPORT calculations. However, it is done once per simulation (END to END statement), not at every time step of a TRANSPORT calculation. Note that the amount of an equilibrium phase does not affect the calculations except when the phase dissolves completely.