Cement Hydration Kinetics

[sama.k]

Hello;
I am running a model using the cement hydration RATE and KINETIC blocks, based on the model proposed by Parrot and Killoh. Parrot and Killoh derived a set of empirical equations to describe the hydration rate, R(t,m), of an individual clinker mineral m at time t (m = C3S, C2S, C3A, C4AF):

Nucleation and growth:     R1(t,m)= (K1/N1) * (1-a(t,m)) * (-ln(1-a(t,m)))^(1-N1)
Diffusion:     R2(t,m)= (K2 * (1-a(t,m)) ^ (2/3)) /  (1-a(t,m)) ^ (1/3)
Hydration shell formation:       R3(t,m)= K3 * (1-a(t,m)) ^ N3


The lowest value of R at curing time t is considered as the controlling step for the nucleation and growth rate. Based on the prediction of these equations, the degree of hydration (a(t)) at time t (day) could then be calculated as:


a(t,m)= a(t-1,m) + ∆t * min(R1(t-1,m), R2(t-1,m), R3(t-1,m)) * Beta * Landa * (A/A0) * exp((Ea(m)/R) * ((1/T) - (1/T0)))


where,  Beta= [1 + 3.333 * (H(m) * (w/c) - a(t-1)]^4          if   a(t) > H(m) * (w/c)
            Beta= 1                                                                 if   a(t) <= H(m) * (w/c)


Landa= ((RH-0.55)/0.45)^4


where H(m) is the critical degree of the clinker mineral m, w/c is the water to cement ratio (0.4), a(t-1,m) is the total hydration degree of  clinker mineral m at the previous time step, A is the Blaine surface area of cement (413 m2/kg), A0 is the reference surface area of cement (385 m2/kg), Ea(m) is the apparent activation energy of clinker mineral m (J/mol), T0 is the reference temperature (293.15 K), and RH is the relative humidity (1%). Other parameters such as K1, N1, K2, K3, N3, H(m) and Ea(m) of each cement clinker mineral are determined.
Here is the input data for my model:

RATES
    C2S
-start
 10 REM   PARM(1) = specific surface area of C2S, m^2/kg C2S
 20 REM   PARM(2) = assumed first rate hydration
 30 a = PARM(2)
 40 FOR a = PARM(2) TO TOTAL_TIME=259200  STEP TIME= TIME + 25920
 50 NEXT
 60 r1 = 0.5*(1-a)
 70 r2 = (((1-a)^0.67)*0.006)/(1-((1-a)^0.34))
 80 r3 = 0.2*((1-a)^5)
 90 b = (((0.54-a)*3.333)+ 1)^ 4
100 min = r1
110 IF min > r2 THEN min = r2
120 IF min > r3 THEN min = r3
130 a1 = a+(TIME/3600*24)*min*b*2.23*EXP(0.008-(2.5/TK))
140 a2 = a+(TIME/3600*24)*min*2.23*EXP(0.008-(2.5/TK))
150 IF a1 > 0.54 THEN k = a1 AND a1 = a1 ELSE k = a2 AND a1 = a2
160 rate = k*(1-SR("C2S"))
170 moles = rate*(TIME/3600*24)
180 SAVE moles
-end
    C3S
-start
 10 REM   PARM(1) = specific surface area of C2S, m^2/kg C3S
 20 REM   PARM(2) = assumed first rate hydration
 30 a = PARM(2)
 40 FOR a = PARM(2) TO TOTAL_TIME=259200  STEP TIME= TIME + 25920
 50 NEXT
 60 r1 = 2.1428*(1-a)* ((-1*log(1-a))^0.3)
 70 r2 = (((1-a)^0.6666)*0.05)/(1-((1-a)^0.3333))
 80 r3 = 1.1*((1-a)^3.3)
 90 b = (((0.72-a)*3.333)+1)^4
100 min = r1
110 IF min > r2 THEN min = r2
120 IF min > r3 THEN min = r3
130 a1 = a+(TIME/3600*24)*min*b*2.23*EXP(0.017-(4.9996/TK))
140 a2 = a+(TIME/3600*24)*min*2.23*EXP(0.017-(4.9996/TK))
150 IF a1 > 0.72 THEN k = a1 AND a1 = a1 ELSE k = a2 AND a1 = a2
160 rate = k*(1-SR("C3S"))
170 moles = rate*(TIME/3600*24)
180 SAVE moles
-end
    C3A
-start
 10 REM   PARM(1) = specific surface area of C3A, m^2/kg C3S
 20 REM   PARM(2) = assumed first rate hydration
 30 a = PARM(2)
 40 FOR a = PARM(2) TO TOTAL_TIME=259200  STEP TIME= TIME + 25920
 50 NEXT
 60 r1 = 1.1764*(1-a)* ((-1*log(1-a))^0.15)
 70 r2 = (((1-a)^0.6666)*0.04)/(1-((1-a)^0.3333))
 80 r3 = ((1-a)^3.3)
 90 b = (((0.64-a)*3.333)+1)^4
100 min = r1
110 IF min > r2 THEN min = r2
120 IF min > r3 THEN min = r3
130 a1 = a+(TIME/3600*24)*min*b*2.23*EXP(0.0221-(6.4994/TK))
140 a2 = a+(TIME/3600*24)*min*2.23*EXP(0.0221-(6.4994/TK))
150 IF a1 > 0.64 THEN k = a1 AND a1 = a1 ELSE k = a2 AND a1 = a2
160 rate = k*(1-SR("C3A"))
170 moles = rate*(TIME/3600*24)
180 SAVE moles
-end
    C4AF
-start
 10 REM   PARM(1) = specific surface area of C4AF, m^2/kg C4AF
 20 REM   PARM(2) = assumed first rate hydration
 30 a = PARM(2)
 40 FOR a = PARM(2) TO TOTAL_TIME=259200  STEP TIME= TIME + 25920
 50 NEXT
 60 r1 = 0.5285 *(1-a)*((-1*LOG(1-a))^0.3)
 70 r2 = (((1-a)^0.6666)*0.015)/(1-((1-a)^0.3333))
 80 r3 = (0.4*(1-a)^3.7)
 90 b = (((0.58-a)*3.333)+1)^4
100 min = r1
110 IF min > r2 THEN min = r2
120 IF min > r3 THEN min = r3
130 a1 = a+(TIME/3600*24)*min*b*2.23*EXP(0.0139-(4.0997/TK))
140 a2 = a+(TIME/3600*24)*min*2.23*EXP(0.0139-(4.0997/TK))
150 IF a1 > 0.58 THEN k = a1 AND a1 = a1 ELSE k = a2 AND a1 = a2
160 rate = k*(1-SR("C4AF"))
170 moles = rate*(TIME/3600*24)
180 SAVE moles
-end
KINETICS 1
C2S
    -formula  (CaO)2SiO2  0.0598
    -m        0.0598
    -m0       0.0598
    -parms    385 0.001
    -tol      0.001
C3A
    -formula  (CaO)3Al2O3  0.0277
    -m        0.0277
    -m0       0.0277
    -parms    385 0.001
    -tol      0.0001
C3S
    -formula  (CaO)3SiO2  0.2912
    -m        0.2912
    -m0       0.2912
    -parms    385 0.001
    -tol      0.0001
C4AF
    -formula  (CaO)4(Al2O3)(Fe2O3)  0.0174
    -m        0.0174
    -m0       0.0174
    -parms    385 0.001
    -tol      0.0001
-steps       259200 in 10 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 15000


Can you please tell me what's wrong with my script for formula of cement hydration rate expression.
Any help or advice will be greatly appreciated.
Many thanks

[dlparkhurst]

You can use PRINT statements to help debug your RATES definitions. For example,

190 PRINT moles, TIME, rate, SR("C2S")

I am surprised this is not an error, but it definitely is not doing what you think.

40 FOR a = PARM(2) TO TOTAL_TIME=259200  STEP TIME= TIME + 25920

TOTAL_TIME and TIME are functions, and cannot be set as variables, so, I think what this statement does is consider (TOTAL_TIME=259200) as a logical value (is TOTAL_TIME equal to 259200, yes 1, no 0), same with (TIME=TIME + 25920), so, I think this statement is probably equivalent to

40 for a PARM(2) TO 0  STEP 1 + 25920

But, to get back to the point, you should not be trying to adjust TIME or TOTAL_TIME, or looping through them. The rate will be integrated over the time step(s) defined in KINETICS, and TIME gives the length of the current subinterval being integrated and TOTAL_TIME is the time at the beginning of that subinterval. I don't understand the point of your loop in the first place (probably should be 50 next a, but you aren't doing anything within the loop anyway.)

150 IF a1 > 0.54 THEN k = a1 AND a1 = a1 ELSE k = a2 AND a1 = a2

Again, (a1 AND a1 = a1) will be interpreted as a logical value, same with a2 AND a1 = a2. I think what you mean is

150 IF a1 > 0.54 THEN k = a1 ELSE k = a2
155 IF a1 > 0.54 THEN a1 = a1 ELSE  a1 = a2

This statement is evaluated as

170 moles = rate*(TIME/3600)*24

That is a start. I think you need a lot of PRINT statements to check your logic.

[sama.k]

Thanks for your help. I revised the model based on your suggestions and comments and also, I used PRINT statements to help debug my RATES definitions.
First, my model takes a long time to execute and even an extended period simulation (such as 28 days or 2419200 seconds) caused the program to hang. I want to know if this is a laptop problem or my model has a problem.
Second, I checked calculations of moles and rates by PRINT statements and I found that moles and rates of cement clinker mineral were calculated correctly before 86400 seconds (1 day) but moles and rates of cement clinker mineral weren't true after 86400 seconds. I'm so confused about what is the problem.

 
For example for C2S:   mole                time                 rate             
                                    1.6379              86400             27.389               it is wrong 
                                 8.0369e-5           258.57           0.001344             it is true
                                 1.9477                  61.193            0.0003257          it is true

Here is the input data for my model:
RATES
    C2S
-start
 10 REM   PARM(1) = specific surface area of C2S, m^2/kg C2S
 20 REM   PARM(2) = assumed first rate hydration
 30 a = PARM(2)
 40 r1 = 0.5*(1-a)
 50 r2 = (((1-a)^0.67)*0.006)/(1-((1-a)^0.34))
 60 r3 = 0.2*((1-a)^5)
 70 b = (((0.54-a)*3.333)+1)^4
 80 min = r1
 90 IF min > r2 THEN min = r2
100 IF min > r3 THEN min = r3
110 a1 = a+(TIME/(3600*24))*min*b*2.23*EXP(0.008-(2.5/TK))
120 a2 = a+(TIME/(3600*24))*min*2.23*EXP(0.008-(2.5/TK))
130 IF a1 > 0.54 THEN k = a1 ELSE k = a2
140 IF a1 > 0.54 THEN a = a1 ELSE a = a2
150 rate = k
160 moles = (rate*M0)
170 SAVE moles
180 PRINT moles, TIME, rate, SR("C2S"), k
-end
    C3S
-start
 10 REM   PARM(1) = specific surface area of C2S, m^2/kg C3S
 20 REM   PARM(2) = assumed first rate hydration
 30 a = PARM(2)
 40 r1 = 2.1428*(1-a)*((-1*LOG(1-a))^0.3)
 50 r2 = (((1-a)^0.6666)*0.05)/(1-((1-a)^0.3333))
 60 r3 = 1.1*((1-a)^3.3)
 70 b = (((0.72-a)*3.333)+1)^4
 80 min = r1
 90 IF min > r2 THEN min = r2
100 IF min > r3 THEN min = r3
110 a1 = a+(TIME/(3600*24))*min*b*2.23*EXP(0.017-(4.9996/TK))
120 a2 = a+(TIME/(3600*24))*min*2.23*EXP(0.017-(4.9996/TK))
130 IF a1 > 0.54 THEN k = a1 ELSE k = a2
140 IF a1 > 0.54 THEN a = a1 ELSE a = a2
150 rate = k
160 moles = (rate*M0)
170 SAVE moles
180 PRINT moles, TIME, rate, SR("C3S"), k
-end
    C3A
-start
 10 REM   PARM(1) = specific surface area of C3A, m^2/kg C3S
 20 REM   PARM(2) = assumed first rate hydration
 30 a = PARM(2)
 40 r1 = 1.1764*(1-a)*((-1*LOG(1-a))^0.15)
 50 r2 = (((1-a)^0.6666)*0.04)/(1-((1-a)^0.3333))
 60 r3 = ((1-a)^3.3)
 70 b = (((0.64-a)*3.333)+1)^4
 80 min = r1
 90 IF min > r2 THEN min = r2
100 IF min > r3 THEN min = r3
110 a1 = a+(TIME/(3600*24))*min*b*2.23*EXP(0.0221-(6.4994/TK))
120 a2 = a+(TIME/(3600*24))*min*2.23*EXP(0.0221-(6.4994/TK))
130 IF a1 > 0.54 THEN k = a1 ELSE k = a2
140 IF a1 > 0.54 THEN a = a1 ELSE a = a2
150 rate = k
160 moles = (rate*M0)
170 SAVE moles
180 PRINT moles, TIME, rate, SR("C3A"), k
-end
    C4AF
-start
 10 REM   PARM(1) = specific surface area of C4AF, m^2/kg C4AF
 20 REM   PARM(2) = assumed first rate hydration
 30 a = PARM(2)
 40 r1 = 0.5285*(1-a)*((-1*LOG(1-a))^0.3)
 50 r2 = (((1-a)^0.6666)*0.015)/(1-((1-a)^0.3333))
 60 r3 = (0.4*(1-a)^3.7)
 70 b = (((0.58-a)*3.333)+1)^4
 80 min = r1
 90 IF min > r2 THEN min = r2
100 IF min > r3 THEN min = r3
110 a1 = a+(TIME/(3600*24))*min*b*2.23*EXP(0.0139-(4.0997/TK))
120 a2 = a+(TIME/(3600*24))*min*2.23*EXP(0.0139-(4.0997/TK))
130 IF a1 > 0.54 THEN k = a1 ELSE k = a2
140 IF a1 > 0.54 THEN a = a1 ELSE a = a2
150 rate = k
160 moles = (rate*M0)
170 SAVE moles
180 PRINT moles, TIME, rate, SR("C4AF"), k
-end
KINETICS 1
C2S
    -formula  (CaO)2SiO2  0.0598
    -m        0.0598
    -m0       0.0598
    -parms    385 1e-005
    -tol      0.001
C3A
    -formula  (CaO)3Al2O3  0.0277
    -m        0.0277
    -m0       0.0277
    -parms    385 0.001
    -tol      0.0001
C3S
    -formula  (CaO)3SiO2  0.2912
    -m        0.2912
    -m0       0.2912
    -parms    385 0.001
    -tol      0.0001
C4AF
    -formula  (CaO)4(Al2O3)(Fe2O3)  0.0174
    -m        0.0174
    -m0       0.0174
    -parms    385 0.001
    -tol      0.0001
-steps       259200 in 3 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 5000

Can you please tell me what's wrong with my script for the formula of cement hydration rate expression.
Any help or advice will be greatly appreciated.
Many thanks

[dlparkhurst]

If you know the rate is wrong, then you should be able to debug it with more PRINT statements. I doubt it is a problem with PHREEQC, the program is simply doing the calculation that you specify.

Note that you do not use parm(1), nor do you use SR in the rate calculation. You probably need a factor like (1 - SR(mineral) to allow the minerals to approach equilibrium.

[sama.k]

The total hydration degree of clinker mineral at each time step (a(t,m)) is calculated based on the total hydration degree of clinker mineral at the previous time step (a(t-1,m)) . could you please assist me how I can update hydration degree parameter at each time step by RATE keyword. The way I've tried to solve it is by using the STEP_NO statement (as below) but it doesn't work. do you have any other suggestions?
And one more question, does the numbering order or numbering steps (i.e. 10, 20, 30, 40 or 10, 20, 100, 110) have any effect on the calculations?

Hydration degree formula: a(t,m)= a(t-1,m) + ∆t * min(R1(t-1,m), R2(t-1,m), R3(t-1,m)) * Beta * Landa * (A/A0) * exp((Ea(m)/R) * ((1/T) - (1/T0)))

Here is the input data for my model but it doesn't work correctly:
10   REM   PARM(1) = specific surface area of C2S, m^2/kg C2S
20   REM   PARM(2) = the first-rate hydration degree 
30   a = PARM(2)
40   n = STEP_NO 
50   counter = 0
60   r1 = 0.5*(1-a)
70   r2 = (((1-a)^0.67)*0.006)/(1-((1-a)^0.34))
80   r3 = 0.2*((1-a)^5)
90   b = (((0.54-a)*3.333)+1)^4
100   min = r1
110   IF min > r2 THEN min = r2
120   IF min > r3 THEN min = r3
130   a1 = a+((TIME/(3600*24))*min*b*0.00538*PARM(1)*EXP(0.008-(2.5/TK))) 
140   a2 = a+(TIME/(3600*24))*min*0.00538*PARM(1)*EXP(0.008-(2.5/TK)) 
150   IF a1 > 0.54 THEN k = a1 ELSE k = a2
160   IF a1 > 0.54 THEN a = a1 ELSE a = a2
170   counter = counter+1
180   IF counter <= n THEN GOTO 60
190   rate = k*(1 – SR("C2S"))
200   moles = rate
210   SAVE moles
220   PRINT moles, TIME, rate, SR("C2S")

Thanks in advance for any help

[dlparkhurst]

You certainly have a way of taking a tangled mess and making it exponentially worse.

My thought is to define two rates, one for the evolution of hydration and one for mineral reaction. It has the advantage that you can work first on the hydration equations before you worry about mineral reaction.

The approach integrates the hydration from a starting value of 0.001. The moles of reactant (M or KIN("C2S_hydration")) represent the hydration degree. In order for it to increase, the SAVEd value must be negative. The C2S_hydration reaction has no effect on solution because the stoichiometric coefficient of H is zero.

I have done my best to interpret your rates, but you should definitely work carefully to get the equations translated to Basic correctly. There seemed to be a problem with rate 2, and rates 1 and 3 are orders of magnitude faster than rate 2. Check also the logic for calculating b.

Once you have the hydration working reasonably, you can add C2S to KINETICS and use the value from C2S_hydration in the C2S RATES definition. For C2S, positive SAVEd values represent dissolution and negative SAVEd values represent precipitation; think of it relative to the solution--positive increases concentrations, negative decreases concentrations.

The line numbers determine the order of the Basic calculations, the statements are executed from lowest to highest, regardless of the order of definition (10 will be executed before 20 even if 20 is defined first). Note that if two statements have the same line number, then only the last definition is used.

Code: [Select]

RATES
C2S_Hydration
-start
10   REM   PARM(1) = specific surface area of C2S, m^2/kg C2S
20   REM   M = the hydration degree
30   a = M
60   r1 = 0.5*(1-a)
70   r2 = (((1-a)^0.67)*0.006)/((1-a)^0.34)
80   r3 = 0.2*((1-a)^5)
90   if (a > 0.54) then b = (((0.54-a)*3.333)+1)^4 else b = 1
95 print r1, r2, r3, b
100   min = r1
110   IF min > r2 THEN min = r2
120   IF min > r3 THEN min = r3
140   a_rate = b*min*0.00538*PARM(1)*EXP(0.008-(2.5/TK))
150   a_delta = a_rate*TIME/86400
210   SAVE -a_delta
220 print a_delta
-end
C2
-start
10 rate = KIN("C2S_Hydration)*(1 - SR("C2S")
20 SAVE rate * time / 86400
-end
SOLUTION
KINETICS
C2S_Hydration
-formula H 0
-M 0.001
-step 8640000 in 10
-parm 385
USER_GRAPH 1
    -axis_titles            "Days" "Hydration degree" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x total_time/86400
20 graph_y KIN("C2S_hydration")
  -end
    -active                 true
END