The silicates of a CEM I clinker, hydrating end to end
Before this page
Coupling kinetics and equilibrium and Hydration kinetics of a CEM I 52.5 R clinker.
A worked application of Kinetics under partial equilibrium: alite and belite dissolve according to Parrott and Killoh (1984), and the hydrate assemblage that forms is computed by Gibbs minimization rather than imposed from a stoichiometric recipe.
Read the theory page first if you have not — this one assumes the partition
What is different from a stoichiometric model
The usual approach writes reactions like
C₃S + 5.3 H → C₁.₇SH₄ + 1.3 CHand integrates the extent. It is serviceable, but the products and their proportions are decided in advance. Nothing tells you the pore solution pH, and nothing adapts if you change the water/cement ratio, add limestone, or leach the paste.
Here the clinker dissolution rate is prescribed — that is genuinely kinetic — and everything downstream follows from thermodynamics: which hydrates appear, in what amounts, and what the pore solution looks like.
Every block below is executed when this page is built, so the numbers are whatever the code actually produced.
1. Species and system
Take the two silicate clinker phases, the hydrates Cemdata18 (Lothenbach et al., 2019) offers for them, and the aqueous species speciation pulls in from the primaries.
using ChemistryLab, DynamicQuantities, OptimaSolver, OrdinaryDiffEq, Printf, Plots
data = datapath("cemdata18-thermofun.json")
substances = build_species(data)
anhydrous = ["C3S", "C2S"]
hydrates = ["Portlandite", "Jennite"]
sp = speciation(substances, vcat(anhydrous, hydrates); aggregate_state = [AS_AQUEOUS])
cs = ChemicalSystem(sp, CEMDATA_PRIMARIES)
@printf "%d species: %d aqueous, %d crystalline\n" length(cs.species) count(
s -> aggregate_state(s) == AS_AQUEOUS, cs.species
) count(s -> aggregate_state(s) == AS_CRYSTAL, cs.species)17 species: 13 aqueous, 4 crystalline2. The initial state
A paste at w/c = 0.40: one kilogram of clinker, 65 % alite and 11 % belite by mass, and 400 g of water.
state0 = ChemicalState(cs)
set_quantity!(state0, "C3S", 0.65u"kg")
set_quantity!(state0, "C2S", 0.11u"kg")
set_quantity!(state0, "H2O@", 0.40u"kg")
for s in anhydrous
@printf "%-5s %8.4f mol\n" s ustrip(us"mol", moles(state0, s))
endC3S 2.8470 mol
C2S 0.6387 mol3. Dissolution reactions
The clinker does not react into hydrates here — it dissolves into ions, and ChemistryLab balances each reaction from the phase and the primaries. The acid-driven form is what drives the pore solution alkaline.
primaries = ["Ca+2", "SiO2@", "H2O@", "H+"]
reactions = ChemistryLab.Reaction[]
for (phase, pk) in (("C3S", PK84_PARAMS_C3S), ("C2S", PK84_PARAMS_C2S))
rxn = Reaction([cs[phase]], [cs[p] for p in primaries]; symbol = "$phase dissolution")
rxn[:rate] = parrott_killoh_avrami(pk, phase; α_max = 1.0, blaine = 380.0u"m^2/kg")
push!(reactions, rxn)
println(" ", rxn)
end (CaO)₃SiO₂ = 3H₂O@ + 3Ca²⁺ + SiO₂@ + (-6)H⁺
(CaO)₂SiO₂ = 2H₂O@ + 2Ca²⁺ + SiO₂@ + (-4)H⁺The negative H⁺ coefficient is the whole point: dissolving one mole of alite consumes six moles of protons, which is what a pore solution at pH 12.5 and above records.
4. The coupled problem
The reaction list is a positional argument, and equilibrium_solver belongs on the problem. Without it the run is a pure kinetics integration and the aqueous phase never re-speciates; with it, respeciate! solves
The activity model matters. A cement pore solution sits at an ionic strength of 0.1–0.7 mol/kg, where the dilute model is not defensible, so pass HKFActivityModel — to the problem and to the solver, which must agree.
model = HKFActivityModel()
kp = KineticsProblem(
cs, reactions, state0, (0.0, 7 * 86400.0);
activity_model = model,
equilibrium_solver = EquilibriumSolver(cs, model, OptimaOptimizer()),
)
sol = integrate(kp, KineticsSolver(; ode_solver = Rodas5P(), reltol = 1.0e-7, abstol = 1.0e-10))
@printf "%d accepted steps, retcode = %s\n" length(sol.t) sol.retcode150 accepted steps, retcode = SuccessWhat the solvers reported during the integration
5. Reading the result
The composition at a given time is not state_at: that returns the purely kinetic reconstruction speciated_states, which re-solves
times = [1.0, 3.0, 7.0] .* 86400
states = speciated_states(sol, kp; times = times)
@printf "%6s %8s %10s %12s %10s\n" "t [d]" "pH" "Jennite" "Portlandite" "H2O"
for (t, st) in zip(times, states)
@printf "%6.1f %8.2f %10.4f %12.4f %10.3f\n" t / 86400 something(pH(st), NaN) (
ustrip(us"mol", moles(st, "Jennite"))
) ustrip(us"mol", moles(st, "Portlandite")) ustrip(us"mol", moles(st, "H2O@"))
end t [d] pH Jennite Portlandite H2O
1.0 12.56 1.4243 1.7971 17.409
3.0 12.56 2.0191 2.5072 15.451
7.0 12.56 2.3864 2.9254 14.262What the solvers reported during the replay
┌ Warning: equilibrium solve returned `MaxIters`; the composition may not be an equilibrium. Set `ChemistryLab.STRICT_CONVERGENCE[] = true` to raise instead.
└ @ ChemistryLab ~/work/ChemistryLab.jl/ChemistryLab.jl/src/equilibrium/equilibrium_solver.jl:677Check the partition constraint — conservation of matter is one of the three conditions a certificate rests on, and the one with an immediate meaning:
p = sol.prob.p
be = collect(sol(times[end])[1:(p.n_be)])
ne = Float64[ustrip(us"mol", moles(states[end], symbol(cs.species[j]))) for j in kp.idx_equilibrium]
@printf "‖Aₑnₑ − bₑ‖∞ = %.2e mol\n" maximum(abs, p.Ae * ne .- be)‖Aₑnₑ − bₑ‖∞ = 1.78e-15 molDegrees of hydration come from the kinetic amounts against their initial values:
α = degrees_of_hydration(sol, kp; times = times)
for s in anhydrous
@printf "α(%s) = %.3f\n" s α[s][end]
endα(C3S) = 0.750
α(C2S) = 0.393The whole week, the degrees read on the integrated trajectory and the assemblage on its certified replay at twenty instants:
grid = 86400 .* 10 .^ range(-1, log10(7); length = 20)
α_grid = degrees_of_hydration(sol, kp; times = grid)
states_grid = speciated_states(sol, kp; times = grid)
d = grid ./ 86400
p_α = plot(d, [α_grid["C3S"] α_grid["C2S"]]; lw = 2, label = ["alite (C3S)" "belite (C2S)"], xscale = :log10,
ylims = (0, 1), xlabel = "age (days)", ylabel = "degree of hydration", legend = :topleft,
title = "w/c 0.40, Blaine 380 m²/kg, 25 °C", titlefontsize = 10)
p_s = plot(; xscale = :log10, xlabel = "age (days)", ylabel = "mol per kg of clinker", legend = :topleft,
title = "Solids, at equilibrium with the solution", titlefontsize = 10)
for s in ("Jennite", "Portlandite", "C3S", "C2S")
plot!(p_s, d, [ustrip(us"mol", moles(st, s)) for st in states_grid]; lw = 2, marker = :circle, ms = 2, label = s)
end
pHs = [something(pH(st), NaN) for st in states_grid]
@printf("pH from %.4f to %.4f over the twenty instants\n", extrema(pHs)...)
plot(p_α, p_s; layout = (1, 2), size = (980, 380), left_margin = 6Plots.mm, bottom_margin = 7Plots.mm)
The pH does not move: portlandite saturates the solution from the first instant shown, and with no alkali in this paste nothing else sets it.
6. What the thermodynamics decided
Nothing above states that portlandite forms, or in what proportion to the C-S-H. The dissolution reactions produce only Ca²⁺, SiO₂ and H⁺; the assemblage is the minimizer of the Gibbs energy under the element totals the kinetics delivered.
st = states[end]
solids = [(symbol(s), ustrip(us"mol", moles(st, symbol(s)))) for s in cs.species
if aggregate_state(s) == AS_CRYSTAL]
filter!(x -> last(x) > 1.0e-6, solids)
sort!(solids; by = last, rev = true)
for (name, amount) in solids
@printf " %-14s %8.4f mol\n" name amount
end Portlandite 2.9254 mol
Jennite 2.3864 mol
C3S 0.7118 mol
C2S 0.3874 molWhat went in, and what is there after a week, in grams:
g_of(st, s) = ustrip(us"g", moles(st, s) * cs[s][:M])
shown_solids = ["C3S", "C2S", "Portlandite", "Jennite"]
xs = collect(eachindex(shown_solids))
fig = bar(xs .- 0.2, [g_of(state0, s) for s in shown_solids]; bar_width = 0.4, label = "at the mixing", color = :grey60,
xticks = (xs, shown_solids), ylabel = "g", legend = :topright, size = (720, 380),
left_margin = 6Plots.mm, bottom_margin = 5Plots.mm, title = "The solids of the paste, at 0 and 7 days")
bar!(fig, xs .+ 0.2, [g_of(st, s) for s in shown_solids]; bar_width = 0.4, label = "after 7 days", color = :steelblue)
plot!(fig; ylims = (0, 1.1 * ylims(fig)[2]))
The pore solution comes with it:
aq = [(symbol(s), ustrip(us"mol", moles(st, symbol(s)))) for s in cs.species
if aggregate_state(s) == AS_AQUEOUS && symbol(s) != "H2O@"]
sort!(aq; by = last, rev = true)
for (name, amount) in aq[1:min(5, end)]
@printf " %-14s %.3e mol\n" name amount
end OH- 9.407e-03 mol
Ca+2 4.149e-03 mol
CaOH+ 1.110e-03 mol
CaSiO3@ 9.469e-06 mol
HSiO3- 2.970e-07 molCaveats worth carrying
The C-S-H is
Jennite, the Cemdata18 end-member normalized per silicon. A real C-S-H is a solid solution of varying Ca/Si; using the CSHQ solid solution inside a coupled run is not exercised by this package's tests.On a cement equilibrium the interior-point optimizer usually stops at its iteration limit, and its return code is not the thing to read.
integratereports the worst element balance instead, separating the accepted steps — the trajectory — from the other solves, which never enter it. The compositions above do not rest on that solver:speciated_statespasses each instant toDualEquilibriumSolver, which solves the KKT system and certifies the result. The problem is convex, so stationarity of the interior species, the component balance, and undersaturation of every absent phase together prove global optimality.A Parrott–Killoh rate reads only its own degree of reaction, so the trajectory here does not depend on the speciation at all; the speciation is what you read out of it. A rate law reading log-activities feeds the speciation back into the trajectory:
integratethen solves the partition at every evaluation of the right-hand side, by the certified route, and the balance is tight at every step, not just at the instants you ask for.
7. The certificate
Nothing above asks you to take the composition on trust. G is convex — an ideal mixing entropy plus terms linear in the amounts of the pure phases, over a polyhedron — so its minimizer is unique and the KKT conditions are sufficient. optimality_certificate checks them.
sub = ChemistryLab._equilibrium_subsystem(kp.system, kp.idx_equilibrium)
des = DualEquilibriumSolver(sub, model)
snames = symbol.(sub.species)
be = collect(sol(times[end])[1:(sol.prob.p.n_be)])
sub_state = ChemicalState(
sub, [ustrip(us"mol", moles(states[end], s)) for s in snames] .* u"mol";
T = sol.prob.p.T_q[], P = sol.prob.p.P_q[],
)
cert = optimality_certificate(des, sub_state; b = be)
@printf "optimal %s\n" cert.optimal
@printf "stationarity %.2e (RT, over %d interior species)\n" cert.stationarity cert.n_interior
@printf "component balance %.2e mol\n" cert.balance
@printf "worst supersaturation %+.2e (negative: every absent phase undersaturated)\n" cert.worst_supersaturationoptimal true
stationarity 1.89e-16 (RT, over 13 interior species)
component balance 1.78e-15 mol
worst supersaturation -Inf (negative: every absent phase undersaturated)Where to go next
The same route applied to a complete CEM I, aluminate and ferrite included, is A complete CEM I 52.5 N, through its pore solution.