How to add a new solution to equilibrate during a kinetic process

[sama.k]

Hello;
I am running a kinetic model that a substance is hydrated in contact with water for 28 days. After 28 days, I want to add a sulfate solution (known as the sulfate attack) to the resulting system (pore solution combination and precipitated phases), and the simulation is continued for 100 days.
Is my approach correct in the below code?

Any help or advice will be greatly appreciated.
Many thanks

Code: [Select]

TITLE <cement hydration- 28 days curing>
PHASES
Anhydrite
    CaSO4 = Ca+2 + SO4-2
    log_k     -4.36
    delta_h   -1.71 kcal
Brucite
    Mg(OH)2 + 2H+ = 2H2O + Mg+2
    log_k     17.07
    delta_h   -115.66 kcal
C3AH6
    Ca3Al2(OH)12 + 12H+ = 2Al+3 + 3Ca+2 + 12H2O
    log_k     82.22
    delta_h   -595.76 kcal
C3FH6
    Ca3Fe2(OH)12 + 12H+ = 3Ca+2 + 2Fe+3 + 12H2O
    log_k     73.65
    delta_h   -516.96 kcal
Calcite
    CaCO3 = CO3-2 + Ca+2
    log_k     -8.48
    delta_h   -2.297 kcal
Ettringite
    Ca6Al2(SO4)3(OH)12:26H2O + 12H+ = 2Al+3 + 6Ca+2 + 38H2O + 3SO4-2
    log_k     57.73
    delta_h   -389.36 kcal
Gypsum
    CaSO4:2H2O = Ca+2 + 2H2O + SO4-2
    log_k     -4.58
    delta_h   -0.109 kcal
Hemicarboaluminate
    Ca4Al2(CO3)0.5(OH)13:5.5H2O + 13H+ = 2Al+3 + 0.5CO3-2 + 4Ca+2 + 18.5H2O
    log_k     87.88
    delta_h   -604.27 kcal
Hydrotalcite
    Mg4Al2(OH)14:3H2O + 14H+ = 2Al+3 + 17H2O + 4Mg+2
    log_k     75.97
    delta_h   -607.91 kcal
Monocarboaluminate
    Ca4Al2(CO3)(OH)12:5H2O + 12H+ = 2Al+3 + CO3-2 + 4Ca+2 + 17H2O
    log_k     71.54
    delta_h   -533.14 kcal
Monosulfoaluminate
    Ca4Al2(SO4)(OH)12:6H2O + 12H+ = 2Al+3 + 4Ca+2 + 18H2O + SO4-2
    log_k     73.68
    delta_h   -553.08 kcal
Portlandite
    Ca(OH)2 + 2H+ = Ca+2 + 2H2O
    log_k     22.79
    delta_h   -129.66 kcal
lime
    CaO + 2H+ = Ca+2 + H2O
    log_k     32.69
    delta_h   -193.91 kJ
periclase
    MgO + 2H+ = H2O + Mg+2
    log_k     21.5841
    delta_h   -151.23 kJ
EQUILIBRIUM_PHASES 1
    Brucite   0 0
    C2S       0 0.0598
    C3A       0 0.0277
    C3AH6     0 0
    C3FH6     0 0
    C3S       0 0.2912
    C4AF      0 0.0174
    Calcite   0 0.0059
    ECSH1-KSH 0 0
    ECSH1-NaSH 0 0
    ECSH1-SH  0 0
    ECSH1-SrSH 0 0
    ECSH1-TobCa 0 0
    Ettringite 0 0
    Gypsum    0 0.018
    Hemicarboaluminate 0 0
    Hydrotalcite 0 0
    Monocarboaluminate 0 0
    Monosulfoaluminate 0 0
    Portlandite 0 0
SOLUTION 1
    temp      20
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    0.04 # kg
REACTION 1
    Calcite    0.0059
    Gypsum     0.018
    K2O        0.0006
    K2SO4      0.0076
    Na2O       0.0053
    Na2SO4     0.0014
    lime       0.0165
    periclase  0.0158
    0.4672 moles in 1 steps
RATES
    C2S_RATE_HYDRATION
-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 C2S, j/Mmol
 40 REM   PARM(4) = R,ideal gas constant, L·Pa·K-1·mol-1
 50 REM   PARM(5) = T0,the reference temperature, K
 60 REM   M = the hydration degree
 70 a = M
 80 r1 = 0.5*(1-a)
 90 r2 = (((1-a)^0.67)*0.006)/(1-((1-a)^0.34))
100 r3 = 0.2*((1-a)^5)
110 IF(a > 0.54) THEN b = (((0.54-a)*3.333)+1)^4 ELSE b = 1
120 min = r1
130 IF min > r2 THEN min = r2
140 IF min > r3 THEN min = r3
150 a_rate = b*min*2.0736*(PARM(1)/PARM(2))*EXP((PARM(3)/PARM(4))*((1/PARM(5))-(1/TK)))
160 a_delta = a_rate*TIME/86400
170 SAVE-a_delta
180 PUT(a_delta, 1)
190 PRINT TIME, a_delta, r1, r2, r3, b
-end
    C2S
-start
10 rate = GET(1) * M0
20 SAVE rate
-end
    C3S_RATE_HYDRATION
-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/Mmol
 40 REM   PARM(4) = R,ideal gas constant, L·Pa·K-1·mol-1
 50 REM   PARM(5) = T0,the reference temperature, K
 60 REM   M = the hydration degree
 70 a = M
 80 r1 = 2.14*(1-a)*((-1*LOG(1-a))^0.3)
 90 r2 = (((1-a)^0.67)*0.05)/(1-((1-a)^0.34))
100 r3 = 1.1*((1-a)^3.3)
110 IF(a > 0.72) THEN b = (((0.72-a)*3.333)+1)^4 ELSE b = 1
120 min = r1
130 IF min > r2 THEN min = r2
140 IF min > r3 THEN min = r3
150 a_rate = b*min*2.0736*(PARM(1)/PARM(2))*EXP((PARM(3)/PARM(4))*((1/PARM(5))-(1/TK)))
160 a_delta = a_rate*TIME/86400
170 SAVE-a_delta
180 PUT(a_delta, 2)
-end
    C3S
-start
10 rate = GET(2) * M0
20 SAVE rate
-end
    C3A_RATE_HYDRATION
-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 C3A, j/Mmol
 40 REM   PARM(4) = R,ideal gas constant, L·Pa·K-1·mol-1
 50 REM   PARM(5) = T0,the reference temperature, K
 60 REM   M = the hydration degree
 70 a = M
 80 r1 = 1.18*(1-a)*((-1*LOG(1-a))^0.15)
 90 r2 = (((1-a)^0.67)*0.04)/(1-((1-a)^0.34))
100 r3 = ((1-a)^3.2)
110 IF(a > 0.64) THEN b = (((0.64-a)*3.333)+1)^4 ELSE b = 1
120 min = r1
130 IF min > r2 THEN min = r2
140 IF min > r3 THEN min = r3
150 a_rate = b*min*2.0736*(PARM(1)/PARM(2))*EXP((PARM(3)/PARM(4))*((1/PARM(5))-(1/TK)))
160 a_delta = a_rate*TIME/86400
170 SAVE-a_delta
180 PUT(a_delta, 3)
-end
    C3A
-start
10 rate = GET(3) * M0
20 SAVE rate
-end
    C4AF_RATE_HYDRATION
-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 C4AF, j/Mmol
 40 REM   PARM(4) = R,ideal gas constant, L·Pa·K-1·mol-1
 50 REM   PARM(5) = T0,the reference temperature, K
 60 REM   M = the hydration degree
 70 a = M
 80 r1 = 0.53*(1-a)*((-1*LOG(1-a))^0.3)
 90 r2 = (((1-a)^0.67)*0.015)/(1-((1-a)^0.34))
100 r3 = (0.4*(1-a)^3.7)
110 IF(a > 0.58) THEN b = (((0.58-a)*3.333)+1)^4 ELSE b = 1
120 min = r1
130 IF min > r2 THEN min = r2
140 IF min > r3 THEN min = r3
150 a_rate = b*min*2.0736*(PARM(1)/PARM(2))*EXP((PARM(3)/PARM(4))*((1/PARM(5))-(1/TK)))
160 a_delta = a_rate*TIME/86400
170 SAVE-a_delta
180 PUT(a_delta, 4)
-end
    C4AF
-start
10 rate = GET(4) * M0
20 SAVE rate
-end
INCREMENTAL_REACTIONS True
KINETICS 1
C2S
    -formula  (CaO)2SiO2  1
    -m        0.0598
    -m0       0.0598
    -tol      1e-008
C3S
    -formula  (CaO)3SiO2  1
    -m        0.2912
    -m0       0.2912
    -tol      1e-008
C2S_RATE_HYDRATION
    -formula  H  0
    -m        0.001
    -m0       0.001
    -parms    413 385 20.785 8.3144 293.15
    -tol      1e-008
C3S_RATE_HYDRATION
    -formula  H  0
    -m        0.001
    -m0       0.001
    -parms    413 385 41.57 8.3144 293.15
    -tol      1e-008
C3A
    -formula  (CaO)3Al2O3  1
    -m        0.0277
    -m0       0.0277
    -tol      1e-008
C3A_RATE_HYDRATION
    -formula  H  0
    -m        0.001
    -m0       0.001
    -parms    413 385 54.04 8.3144 293.15
    -tol      1e-008
C4AF
    -formula  (CaO)4(Al2O3)(Fe2O3)  1
    -m        0.0174
    -m0       0.0174
    -tol      1e-008
C4AF_RATE_HYDRATION
    -formula  H  0
    -m        0.001
    -m0       0.001
    -parms    413 385 34.087 8.3144 293.15
    -tol      1e-008
-steps       2419200 in 28 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 5000
SAVE solution 2
END
TITLE sulfate attack
      sulfate solution= Na2SO4
      "solution" concentration= 44 g/L
      TIME= 100 DAYS
SOLUTION 3
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Na        14.54 g/kgw
    S(6)      30.36 g/kgw
    -water    0.98 # kg
RATES
    Portlandite
-start
 10 REM M = current number of moles of portlandite
 20 REM M0 = number of moles of portlandite initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_port = SI("portlandite")
 70 IF(M <= 0 AND si_port < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_port))
120 moles = rate*TIME
130 SAVE moles
-end
    Monosulfoaluminate
-start
 10 REM M = current number of moles of Monosulfoaluminate
 20 REM M0 = number of moles of Monosulfoaluminate initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_ms = SI("Monosulfoaluminate")
 70 IF(M <= 0 AND si_ms < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_ms))
120 moles = rate*TIME
130 SAVE moles
-end
    Ettringite
-start
 10 REM M = current number of moles of Ettringite
 20 REM M0 = number of moles of Ettringite initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_Et = SI("Ettringite")
 70 IF(M <= 0 AND si_Et < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_Et))
120 moles = rate*TIME
130 SAVE moles
-end
    Hydrotalcite
-start
 10 REM M = current number of moles of Hydrotalcite
 20 REM M0 = number of moles of Hydrotalcite initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_Ht = SI("Hydrotalcite")
 70 IF(M <= 0 AND si_Ht < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_Ht))
120 moles = rate*TIME
130 SAVE moles
-end
    Gypsum
-start
 10 REM M = current number of moles of Gypsum
 20 REM M0 = number of moles of Gypsum initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_Gyp = SI("Gypsum")
 70 IF(M <= 0 AND si_Gyp < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_Gyp))
120 moles = rate*TIME
130 SAVE moles
-end
KINETICS 2
Portlandite
    -formula  Ca(OH)2  -1
    -m        1
    -m0       1
    -parms    16.5 0.67 2.24e-008
    -tol      1e-008
Gypsum
    -formula  CaSO4:2H2O  -1
    -m        1
    -m0       1
    -parms    9.8 0.67 0.0016
    -tol      1e-008
Hydrotalcite
    -formula  Mg4Al2(OH)14:3H2O  -1
    -m        1
    -m0       1
    -parms    9.8 0.67 1e-009
    -tol      1e-008
Monosulfoaluminate
    -formula  Ca4Al2(SO4)(OH)12:6H2O  -1
    -m        1
    -m0       1
    -parms    5.7 0.67 6.76e-012
    -tol      1e-008
Ettringite
    -formula  Ca6Al2(SO4)3(OH)12:26H2O  -1
    -m        1
    -m0       1
    -parms    9.8 0.67 7.08e-013
    -tol      1e-008
C2S
    -formula  (CaO)2SiO2  1
    -m        0.023991
    -m0       0.023991
    -tol      1e-008
C2S_RATE_HYDRATION
    -formula  H  0
    -m        0.60056
    -m0       0.60056
    -parms    413 385 20.785 8.3144 293.15
    -tol      1e-008
C3A
    -formula  (CaO)3Al2O3  1
    -m        0.0059445
    -m0       0.0059445
    -tol      1e-008
C3A_RATE_HYDRATION
    -formula  H  0
    -m        0.78734
    -m0       0.78734
    -parms    413 385 54.04 8.3144 293.15
    -tol      1e-008
C3S
    -formula  (CaO)3SiO2  1
    -m        0.049792
    -m0       0.049792
    -tol      1e-008
C3S_RATE_HYDRATION
    -formula  H  0
    -m        0.83101
    -m0       0.83101
    -parms    413 385 41.57 8.3144 293.15
    -tol      1e-008
C4AF
    -formula  (CaO)4(Al2O3)(Fe2O3)  1
    -m        0.0052423
    -m0       0.0052423
    -tol      1e-008
C4AF_RATE_HYDRATION
    -formula  H  0
    -m        0.70057
    -m0       0.70057
    -parms    413 385 34.078 8.3144 293.15
    -tol      1e-008
-steps       8640000 in 100 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
USER_GRAPH 1
    -headings               Time Portlandite Monosulfoaluminate Ettringite Hydrotalcite Gypsum
    -axis_titles            "Time (days)" "phase weight (percent)" ""
    -chart_title            "phase weight (percent)"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X total_time/86400
20 GRAPH_Y ((equi("Portlandite") * 74) / ((equi("Portlandite") * 74)+( TOT("water")*1000))) * 100, ((equi("Monosulfoaluminate") * 74) / ((equi("Monosulfoaluminate") * 74)+( TOT("water")*1000))) * 100, ((equi("Ettringite ") * 74) / ((equi("Ettringite") * 74)+( TOT("water")*1000))) * 100, ((equi("Hydrotalcite") * 74) / ((equi("Hydrotalcite") * 74)+( TOT("water")*1000))) * 100, ((equi("Gypsum") * 74) / ((equi("Gypsum") * 74)+( TOT("water")*1000))) * 100
  -end
    -active                 true
EQUILIBRIUM_PHASES 2
    Brucite   0 0
    C2S       0 0
    C3A       0 0
    C3AH6     0 0
    C3FH6     0 0
    C3S       0 0
    C4AF      0 0
    Calcite   0 0
    ECSH1-KSH 0 0
    ECSH1-NaSH 0 0
    ECSH1-SH  0 0
    ECSH1-SrSH 0 0
    ECSH1-TobCa 0 0
    Ettringite 0 0
    Gypsum    0 0
    Hemicarboaluminate 0 0
    Hydrotalcite 0 0
    Monocarboaluminate 0 0
    Monosulfoaluminate 0 0
    Portlandite 0 0
END




[dlparkhurst]

I don't think you should define minerals in both KINETICS and EQUILIBRIUM_PHASES. By using EQUILIBRIUM_PHASES, a mineral will instantly dissolve or precipitate to equilibrium (or dissolve entirely); in this case, what is the point of kinetics?

[sama.k]

I modified the code according to your suggestion (Phases in the both KINETICS and EQUILIBRIUM_PHASES or REACTION and EQUILIBRIUM_PHASES were modified).
I used RATES and KINETICS keyword for defining hydration and precipitation of phases dependence to time, and the KINETICS keyword was applied for dissolution of phase dissolved completely as the cement comes into contact with water according to the assumptions. Finally, I need the weight percentage of the phases, pH, and the concentration of the elements relative to the time as output.

[dlparkhurst]

You will have to decide if your input and output match your conceptual model.

Output can be achieved with SELECTED_OUTPUT/USER_PUNCH and USER_GRAPH. Look at the manual for these data blocks, and the section on the Basic Interpreter for Basic functions that provide the information that you want; for instance TOTAL_TIME is the Basic function for the current time in the simulation.

[sama.k]

Thanks for your guidance
My question is about the concept of my model, whether the sulfate solution exactly reacts with the cementitious system produced by 28 days of cement hydration?
When running the bellow codes, no results or errors are shown. the running takes days with no results. So, my model definitely has a problem. Do you have any idea that which part of the model could have a problem and I should focus on it more?

Code: [Select]

TITLE <cement hydration- 28 days curing>
PHASES
Anhydrite
    CaSO4 = Ca+2 + SO4-2
    log_k     -4.36
    delta_h   -1.71 kcal
Brucite
    Mg(OH)2 + 2H+ = 2H2O + Mg+2
    log_k     17.07
    delta_h   -115.66 kcal
C3AH6
    Ca3Al2(OH)12 + 12H+ = 2Al+3 + 3Ca+2 + 12H2O
    log_k     82.22
    delta_h   -595.76 kcal
C3FH6
    Ca3Fe2(OH)12 + 12H+ = 3Ca+2 + 2Fe+3 + 12H2O
    log_k     73.65
    delta_h   -516.96 kcal
Calcite
    CaCO3 = CO3-2 + Ca+2
    log_k     -8.48
    delta_h   -2.297 kcal
Ettringite
    Ca6Al2(SO4)3(OH)12:26H2O + 12H+ = 2Al+3 + 6Ca+2 + 38H2O + 3SO4-2
    log_k     57.73
    delta_h   -389.36 kcal
Gypsum
    CaSO4:2H2O = Ca+2 + 2H2O + SO4-2
    log_k     -4.58
    delta_h   -0.109 kcal
Hemicarboaluminate
    Ca4Al2(CO3)0.5(OH)13:5.5H2O + 13H+ = 2Al+3 + 0.5CO3-2 + 4Ca+2 + 18.5H2O
    log_k     87.88
    delta_h   -604.27 kcal
Hydrotalcite
    Mg4Al2(OH)14:3H2O + 14H+ = 2Al+3 + 17H2O + 4Mg+2
    log_k     75.97
    delta_h   -607.91 kcal
Monocarboaluminate
    Ca4Al2(CO3)(OH)12:5H2O + 12H+ = 2Al+3 + CO3-2 + 4Ca+2 + 17H2O
    log_k     71.54
    delta_h   -533.14 kcal
Monosulfoaluminate
    Ca4Al2(SO4)(OH)12:6H2O + 12H+ = 2Al+3 + 4Ca+2 + 18H2O + SO4-2
    log_k     73.68
    delta_h   -553.08 kcal
Portlandite
    Ca(OH)2 + 2H+ = Ca+2 + 2H2O
    log_k     22.79
    delta_h   -129.66 kcal
lime
    CaO + 2H+ = Ca+2 + H2O
    log_k     32.69
    delta_h   -193.91 kJ
periclase
    MgO + 2H+ = H2O + Mg+2
    log_k     21.5841
    delta_h   -151.23 kJ
EQUILIBRIUM_PHASES 1
    Ettringite 0 0
    Hemicarboaluminate 0 0
    Hydrotalcite 0 0
    Monocarboaluminate 0 0
    Monosulfoaluminate 0 0
    Portlandite 0 0
    CSHQ-JenD 0 0
    CSHQ-JenH 0 0
    CSHQ-TobD 0 0
    CSHQ-TobH 0 0
SOLUTION 1
    temp      20
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    0.04 # kg
REACTION 1
    Calcite    0.0059
    Gypsum     0.018
    K2O        0.0006
    K2SO4      0.0076
    Na2O       0.0053
    Na2SO4     0.0014
    lime       0.0165
    periclase  0.0158
    0.4672 moles in 1 steps
RATES
    C2S_RATE_HYDRATION
-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 C2S, j/Mmol
 40 REM   PARM(4) = R,ideal gas constant, L·Pa·K-1·mol-1
 50 REM   PARM(5) = T0,the reference temperature, K
 60 REM   M = the hydration degree
 70 a = M
 80 r1 = 0.5*(1-a)
 90 r2 = (((1-a)^0.67)*0.006)/(1-((1-a)^0.34))
100 r3 = 0.2*((1-a)^5)
110 IF(a > 0.54) THEN b = (((0.54-a)*3.333)+1)^4 ELSE b = 1
120 min = r1
130 IF min > r2 THEN min = r2
140 IF min > r3 THEN min = r3
150 a_rate = b*min*2.0736*(PARM(1)/PARM(2))*EXP((PARM(3)/PARM(4))*((1/PARM(5))-(1/TK)))
160 a_delta = a_rate*TIME/86400
170 SAVE-a_delta
180 PUT(a_delta, 1)
190 PRINT TIME, a_delta, r1, r2, r3, b
-end
    C2S
-start
10 rate = GET(1) * M0
20 SAVE rate
-end
    C3S_RATE_HYDRATION
-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/Mmol
 40 REM   PARM(4) = R,ideal gas constant, L·Pa·K-1·mol-1
 50 REM   PARM(5) = T0,the reference temperature, K
 60 REM   M = the hydration degree
 70 a = M
 80 r1 = 2.14*(1-a)*((-1*LOG(1-a))^0.3)
 90 r2 = (((1-a)^0.67)*0.05)/(1-((1-a)^0.34))
100 r3 = 1.1*((1-a)^3.3)
110 IF(a > 0.72) THEN b = (((0.72-a)*3.333)+1)^4 ELSE b = 1
120 min = r1
130 IF min > r2 THEN min = r2
140 IF min > r3 THEN min = r3
150 a_rate = b*min*2.0736*(PARM(1)/PARM(2))*EXP((PARM(3)/PARM(4))*((1/PARM(5))-(1/TK)))
160 a_delta = a_rate*TIME/86400
170 SAVE-a_delta
180 PUT(a_delta, 2)
-end
    C3S
-start
10 rate = GET(2) * M0
20 SAVE rate
-end
    C3A_RATE_HYDRATION
-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 C3A, j/Mmol
 40 REM   PARM(4) = R,ideal gas constant, L·Pa·K-1·mol-1
 50 REM   PARM(5) = T0,the reference temperature, K
 60 REM   M = the hydration degree
 70 a = M
 80 r1 = 1.18*(1-a)*((-1*LOG(1-a))^0.15)
 90 r2 = (((1-a)^0.67)*0.04)/(1-((1-a)^0.34))
100 r3 = ((1-a)^3.2)
110 IF(a > 0.64) THEN b = (((0.64-a)*3.333)+1)^4 ELSE b = 1
120 min = r1
130 IF min > r2 THEN min = r2
140 IF min > r3 THEN min = r3
150 a_rate = b*min*2.0736*(PARM(1)/PARM(2))*EXP((PARM(3)/PARM(4))*((1/PARM(5))-(1/TK)))
160 a_delta = a_rate*TIME/86400
170 SAVE-a_delta
180 PUT(a_delta, 3)
-end
    C3A
-start
10 rate = GET(3) * M0
20 SAVE rate
-end
    C4AF_RATE_HYDRATION
-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 C4AF, j/Mmol
 40 REM   PARM(4) = R,ideal gas constant, L·Pa·K-1·mol-1
 50 REM   PARM(5) = T0,the reference temperature, K
 60 REM   M = the hydration degree
 70 a = M
 80 r1 = 0.53*(1-a)*((-1*LOG(1-a))^0.3)
 90 r2 = (((1-a)^0.67)*0.015)/(1-((1-a)^0.34))
100 r3 = (0.4*(1-a)^3.7)
110 IF(a > 0.58) THEN b = (((0.58-a)*3.333)+1)^4 ELSE b = 1
120 min = r1
130 IF min > r2 THEN min = r2
140 IF min > r3 THEN min = r3
150 a_rate = b*min*2.0736*(PARM(1)/PARM(2))*EXP((PARM(3)/PARM(4))*((1/PARM(5))-(1/TK)))
160 a_delta = a_rate*TIME/86400
170 SAVE-a_delta
180 PUT(a_delta, 4)
-end
    C4AF
-start
10 rate = GET(4) * M0
20 SAVE rate
-end
INCREMENTAL_REACTIONS True
KINETICS 1
C2S
    -formula  (CaO)2SiO2  1
    -m        0.0598
    -m0       0.0598
    -tol      1e-008
C3S
    -formula  (CaO)3SiO2  1
    -m        0.2912
    -m0       0.2912
    -tol      1e-008
C2S_RATE_HYDRATION
    -formula  H  0
    -m        0.001
    -m0       0.001
    -parms    413 385 20.785 8.3144 293.15
    -tol      1e-008
C3S_RATE_HYDRATION
    -formula  H  0
    -m        0.001
    -m0       0.001
    -parms    413 385 41.57 8.3144 293.15
    -tol      1e-008
C3A
    -formula  (CaO)3Al2O3  1
    -m        0.0277
    -m0       0.0277
    -tol      1e-008
C3A_RATE_HYDRATION
    -formula  H  0
    -m        0.001
    -m0       0.001
    -parms    413 385 54.04 8.3144 293.15
    -tol      1e-008
C4AF
    -formula  (CaO)4(Al2O3)(Fe2O3)  1
    -m        0.0174
    -m0       0.0174
    -tol      1e-008
C4AF_RATE_HYDRATION
    -formula  H  0
    -m        0.001
    -m0       0.001
    -parms    413 385 34.087 8.3144 293.15
    -tol      1e-008
-steps       2419200 in 28 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 5000
SAVE solution 2
END
TITLE sulfate attack
      sulfate solution= Na2SO4
      "solution" concentration= 44 g/L
      TIME= 100 DAYS
SOLUTION 3
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Na        14.54 g/kgw
    S(6)      30.36 g/kgw
    -water    0.98 # kg
RATES
    Portlandite
-start
 10 REM M = current number of moles of portlandite
 20 REM M0 = number of moles of portlandite initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_port = SI("portlandite")
 70 IF(M <= 0 AND si_port < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_port))
120 moles = rate*TIME
130 SAVE moles
-end
    Monosulfoaluminate
-start
 10 REM M = current number of moles of Monosulfoaluminate
 20 REM M0 = number of moles of Monosulfoaluminate initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_ms = SI("Monosulfoaluminate")
 70 IF(M <= 0 AND si_ms < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_ms))
120 moles = rate*TIME
130 SAVE moles
-end
    Ettringite
-start
 10 REM M = current number of moles of Ettringite
 20 REM M0 = number of moles of Ettringite initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_Et = SI("Ettringite")
 70 IF(M <= 0 AND si_Et < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_Et))
120 moles = rate*TIME
130 SAVE moles
-end
    Hydrotalcite
-start
 10 REM M = current number of moles of Hydrotalcite
 20 REM M0 = number of moles of Hydrotalcite initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_Ht = SI("Hydrotalcite")
 70 IF(M <= 0 AND si_Ht < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_Ht))
120 moles = rate*TIME
130 SAVE moles
-end
    Gypsum
-start
 10 REM M = current number of moles of Gypsum
 20 REM M0 = number of moles of Gypsum initially present
 30 REM parm(1) = A/V, m^2/L
 40 REM parm (2) = exponent for M/M0
 50 REM parm (3) = rate constant (rf), (mol.m-2.s-1)
 60 si_Gyp = SI("Gypsum")
 70 IF(M <= 0 AND si_Gyp < 0) THEN GOTO 130
 80 IF M0 > 0 THEN t = M/M0
 90 IF t = 0 THEN t = 1
100 area = PARM(1)*(t)^PARM(2)
110 rate = area*1e-3*PARM(3)*(1-10^(si_Gyp))
120 moles = rate*TIME
130 SAVE moles
-end
KINETICS 2
Portlandite
    -formula  Ca(OH)2  -1
    -m        1
    -m0       1
    -parms    16.5 0.67 2.24e-008
    -tol      1e-008
Gypsum
    -formula  CaSO4:2H2O  -1
    -m        1
    -m0       1
    -parms    9.8 0.67 0.0016
    -tol      1e-008
Hydrotalcite
    -formula  Mg4Al2(OH)14:3H2O  -1
    -m        1
    -m0       1
    -parms    9.8 0.67 1e-009
    -tol      1e-008
Monosulfoaluminate
    -formula  Ca4Al2(SO4)(OH)12:6H2O  -1
    -m        1
    -m0       1
    -parms    5.7 0.67 6.76e-012
    -tol      1e-008
Ettringite
    -formula  Ca6Al2(SO4)3(OH)12:26H2O  -1
    -m        1
    -m0       1
    -parms    9.8 0.67 7.08e-013
    -tol      1e-008
C2S
    -formula  (CaO)2SiO2  1
    -m        0.023991
    -m0       0.023991
    -tol      1e-008
C2S_RATE_HYDRATION
    -formula  H  0
    -m        0.60056
    -m0       0.60056
    -parms    413 385 20.785 8.3144 293.15
    -tol      1e-008
C3A
    -formula  (CaO)3Al2O3  1
    -m        0.0059445
    -m0       0.0059445
    -tol      1e-008
C3A_RATE_HYDRATION
    -formula  H  0
    -m        0.78734
    -m0       0.78734
    -parms    413 385 54.04 8.3144 293.15
    -tol      1e-008
C3S
    -formula  (CaO)3SiO2  1
    -m        0.049792
    -m0       0.049792
    -tol      1e-008
C3S_RATE_HYDRATION
    -formula  H  0
    -m        0.83101
    -m0       0.83101
    -parms    413 385 41.57 8.3144 293.15
    -tol      1e-008
C4AF
    -formula  (CaO)4(Al2O3)(Fe2O3)  1
    -m        0.0052423
    -m0       0.0052423
    -tol      1e-008
C4AF_RATE_HYDRATION
    -formula  H  0
    -m        0.70057
    -m0       0.70057
    -parms    413 385 34.078 8.3144 293.15
    -tol      1e-008
-steps       8640000 in 100 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
USER_GRAPH 1
    -headings               Time Portlandite Monosulfoaluminate Ettringite Hydrotalcite Gypsum
    -axis_titles            "Time (days)" "phase weight (percent)" ""
    -chart_title            "phase weight (percent)"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X total_time/86400
20 GRAPH_Y ((equi("Portlandite") * 74) / ((equi("Portlandite") * 74)+( TOT("water")*1000))) * 100, ((equi("Monosulfoaluminate") * 74) / ((equi("Monosulfoaluminate") * 74)+( TOT("water")*1000))) * 100, ((equi("Ettringite ") * 74) / ((equi("Ettringite") * 74)+( TOT("water")*1000))) * 100, ((equi("Hydrotalcite") * 74) / ((equi("Hydrotalcite") * 74)+( TOT("water")*1000))) * 100, ((equi("Gypsum") * 74) / ((equi("Gypsum") * 74)+( TOT("water")*1000))) * 100
  -end
    -active                 true
EQUILIBRIUM_PHASES 2
    Ettringite 0 0
    Gypsum    0 0
    Hemicarboaluminate 0 0
    Hydrotalcite 0 0
    Monocarboaluminate 0 0
    Monosulfoaluminate 0 0
    Portlandite 0 0
    CSHQ-JenD 0 0
    CSHQ-JenH 0 0
    CSHQ-TobD 0 0
    CSHQ-TobH 0 0
END

[dlparkhurst]

Try using -cvode in the KINETICS definition.

However, let's start with your first simulation. I don't see any sink for Si, it just continuously increases to over 15 molal during the first KINETIC reaction. Shouldn't some Ca-Si (or even something like chalcedony) start to form? You can add the following to start looking at the evolution of concentrations.

Code: [Select]

USER_GRAPH 1
    -headings               Days Si
    -axis_titles            "Days" "Moles per kilogram water" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X total_time/86400
20 GRAPH_Y TOT("Si")
  -end