A complete CEM I 52.5 N, through its pore solution
Before this page
The silicates of a CEM I clinker, hydrating end to end runs the two silicate clinker phases this way: prescribe the dissolution, let the thermodynamics decide the hydrates. This page does the same for a complete CEM I — alite, belite, aluminate, ferrite, gypsum and limestone filler — and then reads its calorimetry off the result.
The aluminates are what make the difference. A stoichiometric model has to state, in advance and by hand, that the aluminate goes to ettringite while sulfate lasts, then to monocarboaluminate if carbonate is available, then to monosulphate, then to hydrogarnet. Here none of that is written anywhere. The clinker only dissolves, into Ca²⁺, SiO₂, AlO₂⁻, FeO₂⁻, SO₄²⁻, CO₃²⁻ and H⁺, and a Gibbs minimization at every accepted step decides what is stable. §4 runs the same paste twice, with and without limestone, and gets two different aluminate histories out of one model.
The model lives in scripts/ionic_hydration.jl, which this page includes, so the page and the script cannot drift apart.
Every composition below is proved optimal
speciated_states passes each instant to DualEquilibriumSolver, which solves the KKT system and returns a certificate. The Gibbs problem is convex — an ideal mixing entropy plus terms linear in the amounts of the pure phases, over a polyhedron — so stationarity of the interior species, the component balance, and undersaturation of every absent phase together prove global optimality. All forty instants below are certified, with element balances between 1e-15 and 1e-11 mol.
The interior-point solve alone would not support that claim: on this package's own calcite reference it returns pH 6.96 where the certified solve gives 9.90, and on these problems its return code does not distinguish the two cases.
1. The system and the dissolution reactions
using ChemistryLab, DynamicQuantities, OptimaSolver, OrdinaryDiffEq
using OrderedCollections, Printf, Plots
gr()
include(joinpath(pkgdir(ChemistryLab), "scripts", "ionic_hydration.jl"))
cs = build_ionic_system(:opc)
@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)67 species: 54 aqueous, 13 crystallineThe dissolution reactions are balanced by ChemistryLab from the phase and the primaries, and come out in the expected acid-driven form:
for r in ionic_reactions(cs; wb = 0.5, blaine = 380.0u"m^2/kg")[1:3]
println(" ", r.reaction)
end (CaO)₃SiO₂ = 3H₂O@ + 3Ca²⁺ + SiO₂@ + (-6)H⁺
(CaO)₂SiO₂ = 2H₂O@ + 2Ca²⁺ + SiO₂@ + (-4)H⁺
(CaO)₃Al₂O₃ = 2H₂O@ + 3Ca²⁺ + 2AlO₂⁻ + (-4)H⁺Note the negative H⁺ coefficients: dissolving alite consumes six protons per mole, which is what drives the pore solution to pH 12.5 and above.
2. Running the coupling
The formulation is the CEM I 52.5 N of Lavergne et al. (2018), Table 9: Bogue composition C₃S 65 / C₂S 11 / C₃A 11 / C₄AF 8, gypsum 4.6 %, calcite 3.5 %, Blaine 380 m²/kg, w/b = 0.50, one kilogram of binder so every extensive result is per kilogram.
# Table 9 of Lavergne et al. (2018), read by scripts/ionic_hydration.jl from
# data/literature/Lavergne2018.json
IONIC_CEMENT.clinker # (C3S = 0.65, C2S = 0.11, C3A = 0.11, C4AF = 0.08)
TEND = 28 * 86400.0
run_cal = run_ionic_hydration(;
wb = 0.5, IONIC_CEMENT.clinker, IONIC_CEMENT.gypsum, IONIC_CEMENT.filler, tend = TEND,
)One Gibbs minimization per accepted step, so this is the expensive part — 216 accepted steps, some four minutes. The activity model is HKFActivityModel on both halves of the coupling: a cement pore solution sits at I ≈ 0.1–0.7 mol/kg, where a dilute model is not defensible.
This run is performed by the build, and it is the expensive part
The call above is shown and not executed here; it is made by scripts/precomputed.jl, which the next block includes. That script memoizes per process, so this page's phase history and its calorimetry come from one integration rather than two — Documenter runs the whole site in a single process, which is usually a hazard and here is the thing that makes this affordable.
Nothing is approximated: the trajectory is integrated, then replayed and certified instant by instant, by this build. It used to be read from a stored file — a coupled equilibrium cost 583 ms and the site called for thousands of them — and two solver fixes later it costs 17 ms, which makes computing it affordable again. That is worth the build time: a stored result is a claim about code that may since have changed, and keeping the two in step needed a guard, a procedure, and a list of traps.
include(joinpath(pkgdir(ChemistryLab), "scripts", "precomputed.jl"))
# The two pastes this page compares are two independent 28-day integrations,
# each followed by a certified replay at every reported instant. Nothing is
# shared between them, so they are computed together rather than one after the
# other; on a single-threaded session this is exactly the same work in the same
# order. Every `read_precomputed` below is then a cache hit.
warm_precomputed(["ionic_opc_phases", "ionic_nolimestone_phases"])
phases_c = read_precomputed("ionic_opc_phases")
heat_c = read_precomputed("ionic_opc_heat")
for line in phases_c.provenance
println(" ", line)
end CEM I with 3.5 % limestone filler -- volume fractions of the phase families, pH, porosity
computed by this documentation build -- no stored result is read
ChemistryLab version: 0.42.0, OptimaSolver version: 0.10.0
commit: 73c59df23
40 log-spaced instants over 28.00 days
w/b 0.50, Blaine default, 20 C, 1 kg of binder
clinker (mass fractions): (C3S = 0.65, C2S = 0.11, C3A = 0.11, C4AF = 0.08)
certified: 40 of 40 replayed instants proved optimal against the KKT conditionsWhat the solvers reported while computing these trajectories
┌ Warning: 35 equilibrium solve(s) stopped short of the optimizer's tolerance and were used anyway. Judge them on the element balance, not on the retcode: worst |Aₑn − bₑ|∞ over the run was 1.62e-11 mol ON THE ACCEPTED STEPS, i.e. on the trajectory itself; 1.62e-11 mol counting also the Jacobian probes and rejected steps, which never enter the solution. Read the first figure: 1e-10 mol is machine precision whatever the system, while 1e-2 mol against a 0.3 mol sulfate budget is not. How much it matters depends on the RATE LAWS: `bₑ` is integrated from the rates alone, so a law that reads only its own degree of reaction (Parrott-Killoh, Waller) gives a trajectory independent of the speciation, and this figure then bears on the reported composition only — recover that with `speciated_states`, which certifies each instant against the KKT conditions. A law reading log-activities (a saturation ratio) does feed the speciation back into the trajectory, and there this figure is a direct measure of the error. Do NOT simply loosen the optimizer tolerance — on the calcite reference case that degrades the speciation from 4 % to 250 % against Reaktoro.
└ @ KineticsOrdinaryDiffEqExt ~/work/ChemistryLab.jl/ChemistryLab.jl/ext/KineticsOrdinaryDiffEqExt.jl:216The columns come back as the page's own variables, so everything downstream is the same calculation it always was:
families = phase_families(phases_c)
TIMES = phases_c.columns["time_s"]
times = TIMES
pore_pH = phases_c.columns["pore_pH"]
poro = [
(
liquid = phases_c.columns["poro_liquid"][i],
void = phases_c.columns["poro_void"][i],
total = phases_c.columns["poro_total"][i],
) for i in eachindex(TIMES)
]
fracs = [Dict(f => phases_c.columns[f][i] for f in families) for i in eachindex(TIMES)]3. The pore solution
Certifying forty instants is the dominant cost of the whole page, and the calorimetry of §5 needs the very same compositions, so the replay is done once here and handed on.
p_pH = plot(
times ./ 86400, pore_pH; xscale = :log10, lw = 2, legend = false,
xlabel = "time [days]", ylabel = "pore-solution pH",
title = "Pore solution — a computed quantity",
ylims = (11.5, 13.5), size = (760, 380),
)
The solution holds at pH 12.7 throughout, which is where a Portland cement pore solution belongs. Nothing in the input fixes it.
The sequence of instants is not optional
Each speciation warm-starts from the previous one, and the guess is first clipped to the element budget and projected back into {Aₑn = bₑ, n ≥ 0}. Both matter: the warm start is the equilibrium of the previous bₑ, so once the sulfate is spent it demands more of it than now exists, and the interior-point solve would start outside its own feasible set. Solved instead from a cold guess, the map returns no hydrates at all and a pore solution at pH 6, while the run itself computed 2.2 mol of C-S-H.
speciated_states also walks the chain up to the first instant requested, through a few earlier times whose compositions are discarded. That instant has no predecessor of its own, and without the run-up its interior-point answer held 56 interior species where the answer has 25 — every candidate hydrate present, four of them at 1e-5 to 1e-6 mol — which no certifying Newton recovers from, since it inherits the start.
3bis. The paste, from cast to 28 days
The figure a cement chemist reads first: what the volume of the specimen is made of, at every instant, stacked from the anhydrous grains at the bottom to the pore water and the empty porosity at the top. Each band is a volume fraction of the fresh specimen, so the total stays at one and the Le Chatelier contraction appears as the void band opening rather than as a shrinking total.
Nothing in this figure was imposed. The clinker was told only to dissolve; that the C-S-H and the portlandite grow together, that the ettringite forms early and then gives way, that the empty porosity opens after set — all of it is what the Gibbs minimization chose at each of the eighty instants.
cols = family_colors()
shown = [f for f in families if maximum(phases_c.columns[f]) > 1.0e-4]
p_stack = plot(;
xscale = :log10, xlabel = "time [days]",
ylabel = "volume fraction of the fresh paste",
title = "CEM I, w/b = 0.50 — the paste from cast to 28 days",
legend = :outerright, size = (900, 440), ylims = (0, 1),
left_margin = 8Plots.mm, bottom_margin = 8Plots.mm,
)
areaplot!(
p_stack, times ./ 86400,
hcat((phases_c.columns[f] for f in shown)...);
label = permutedims(shown), color = permutedims([cols[f] for f in shown]),
fillalpha = 0.85, linewidth = 0,
)
The same information at the two ends, which is what a mix design is actually compared on — what went in, and what it became:
first_i, last_i = 1, length(times)
bars = [f for f in families if max(phases_c.columns[f][first_i],
phases_c.columns[f][last_i]) > 1.0e-3]
# Two series side by side, with `Plots` alone: the x positions are shifted by
# half a bar width rather than reaching for a grouped-bar package.
x = collect(eachindex(bars))
w = 0.38
p_bar = bar(
x .- w / 2, [phases_c.columns[f][first_i] for f in bars];
bar_width = w, label = "at 1.2 h", color = :grey70,
xticks = (x, bars), xrotation = 30,
ylabel = "volume fraction of the fresh paste",
title = "What went in, and what it became",
size = (860, 400), legend = :topright,
left_margin = 8Plots.mm, bottom_margin = 12Plots.mm,
)
bar!(
p_bar, x .+ w / 2, [phases_c.columns[f][last_i] for f in bars];
bar_width = w, label = "at 28 days", color = :steelblue,
)
@printf("%-12s %10s %10s\n", "family", "1.2 h", "28 d")
for f in bars
@printf("%-12s %10.4f %10.4f\n", f,
phases_c.columns[f][first_i], phases_c.columns[f][last_i])
end
@printf("%-12s %10.4f %10.4f\n", "TOTAL",
sum(phases_c.columns[f][first_i] for f in families),
sum(phases_c.columns[f][last_i] for f in families))family 1.2 h 28 d
anhydrous 0.3387 0.0676
gypsum 0.0160 0.0000
calcite 0.0104 0.0000
C-S-H 0.0000 0.3516
CH 0.0000 0.1288
AFt 0.0098 0.0778
AFm 0.0000 0.1089
FH3 0.0000 0.0092
water 0.6129 0.1867
void 0.0122 0.0693
TOTAL 1.0000 1.0000The total holding at one to four decimals is not a normalization: the volume fractions are computed independently against the fresh specimen, so their sum closing is a check that nothing was created or lost.
4. The aluminate sequence, and how to deplete the ettringite
The same paste, run a second time with the limestone removed and nothing else changed.
phases_n = read_precomputed("ionic_nolimestone_phases")
heat_n = read_precomputed("ionic_nolimestone_heat")
families_n = phase_families(phases_n)
fracs_nl = [
Dict(f => phases_n.columns[f][i] for f in families_n) for i in eachindex(TIMES)
]
td = times ./ 86400
p_seq = plot(;
xscale = :log10, xlabel = "time [days]", ylabel = "volume fraction",
title = "Aluminates: a result, not an input", legend = :topleft, size = (760, 420),
)
plot!(p_seq, td, [get(f, "AFt", 0.0) for f in fracs]; lw = 2, color = 1, label = "AFt — 3.5 % calcite")
plot!(p_seq, td, [get(f, "AFm", 0.0) for f in fracs]; lw = 2, color = 2, label = "AFm — 3.5 % calcite")
plot!(p_seq, td, [get(f, "AFt", 0.0) for f in fracs_nl]; lw = 2, ls = :dash, color = 1, label = "AFt — no limestone")
plot!(p_seq, td, [get(f, "AFm", 0.0) for f in fracs_nl]; lw = 2, ls = :dash, color = 2, label = "AFm — no limestone")
for (lbl, f) in ("with 3.5 % calcite" => fracs, "no limestone" => fracs_nl)
aft = [get(fi, "AFt", 0.0) for fi in f]
j = argmax(aft)
@printf "%-20s AFt peak %.4f at %5.2f d AFt(28 d) %.4f AFm(28 d) %.4f\n" lbl aft[j] (
times[j] / 86400
) aft[end] get(f[end], "AFm", 0.0)
endwith 3.5 % calcite AFt peak 0.0778 at 28.00 d AFt(28 d) 0.0778 AFm(28 d) 0.1089
no limestone AFt peak 0.0717 at 0.35 d AFt(28 d) 0.0000 AFm(28 d) 0.1028Two histories out of one model:
With limestone the ettringite forms and survives to 28 days. The carbonate reacts with the aluminate to form monocarboaluminate, so the sulfate is never called upon to feed an AFm and the AFt is stabilized. This is the well-known limestone effect.
Without limestone the ettringite peaks within hours and is then depleted, reaching zero by 28 days. Once the gypsum is exhausted the AFt is the only sulfate reservoir left, and the remaining aluminate converts it to monosulphate.
Nothing in the code distinguishes the two cases. The free energy does.
5. Calorimetry
What is computed, and why it is not the heat of the reactions
Under partial equilibrium the kinetic reactions only dissolve the clinker into ions and the hydrates are precipitated by the minimization, so the heat is the fall of the enthalpy of the whole composition, Eqs. (17)–(21) of (Lavergne et al., 2018), not the heat of the kinetic reactions (What the package counts as heat, The calorimeters). heat_release reads it on the certified speciations of §3, which also makes the curves monotone:
t_cal = heat_c.columns["time_s"]
# The stored curves are already per gram of binder; the page below works in
# joules per kilogram, so they are scaled back to it rather than the reverse.
BINDER_G = 1000.0 # the runs simulate 1 kg of binder
Q_c = heat_c.columns["Q_J_per_g"] .* BINDER_G
qd_c = heat_c.columns["heat_flow_W_per_g"] .* BINDER_G
Q_n = heat_n.columns["Q_J_per_g"] .* BINDER_G
qd_n = heat_n.columns["heat_flow_W_per_g"] .* BINDER_G
@printf "monotone: with limestone %s, without %s\n" all(diff(Q_c) .>= -1.0e-9) all(
diff(Q_n) .>= -1.0e-9
)monotone: with limestone true, without trueIsothermal calorimetry at 20 °C
p_Q = plot(;
xscale = :log10, xlabel = "time [days]", ylabel = "Q [J / g of binder]",
title = "Heat released, isothermal at 20 °C", legend = :topleft, size = (760, 420),
)
plot!(p_Q, t_cal ./ 86400, Q_c ./ BINDER_G; lw = 2, color = 1, label = "with 3.5 % calcite")
plot!(p_Q, t_cal ./ 86400, Q_n ./ BINDER_G; lw = 2, color = 2, ls = :dash, label = "no limestone")
p_q = plot(;
xscale = :log10, xlabel = "time [days]", ylabel = "q̇ [mW / g of binder]",
title = "Heat rate", legend = :topright, size = (760, 420),
)
plot!(p_q, t_cal ./ 86400, qd_c ./ BINDER_G .* 1000; lw = 2, color = 1, label = "with 3.5 % calcite")
plot!(p_q, t_cal ./ 86400, qd_n ./ BINDER_G .* 1000; lw = 2, color = 2, ls = :dash, label = "no limestone")
Both axes above are logarithmic in time, which is the scale a calorimetrist reads: it opens out the induction period and the peak, which occupy the first decade and would be a single vertical rise otherwise. It also flatters the late curve, where nothing much happens and a decade of time is a centimeter of paper.
The linear axis answers the other question — how much of the total is already released at a date on the calendar, and how flat the curve has become. Same arrays, no recomputation:
p_Qlin = plot(;
xlabel = "time [days]", ylabel = "Q [J / g of binder]",
title = "Heat released, linear time", legend = :bottomright, size = (760, 420),
)
plot!(p_Qlin, t_cal ./ 86400, Q_c ./ BINDER_G; lw = 2, color = 1, label = "with 3.5 % calcite")
plot!(p_Qlin, t_cal ./ 86400, Q_n ./ BINDER_G; lw = 2, color = 2, ls = :dash, label = "no limestone")
vline!(p_Qlin, [1, 7, 28]; ls = :dot, color = :gray, label = "1, 7, 28 days")
And the pairing matters here, because the two curves separate late, not early. From the figures printed above:
| 1 day | 7 days | 28 days | peak | |
|---|---|---|---|---|
| with 3.5 % calcite | 179.1 | 354.0 | 428.3 | 3.89 mW/g at 7.15 h |
| no limestone | 177.6 | 342.7 | 413.1 | 3.93 mW/g at 7.15 h |
1.5 J/g apart at one day — eight tenths of a percent, invisible — and 15.2 J/g apart at twenty-eight, which is 3.7 % and plainly visible. The peaks sit at the same time, 7.15 h, and at nearly the same height.
So the limestone does not shift when the heat comes out; it adds to how much, and it does so gradually, over the weeks the logarithmic axis compresses into its last centimeter. A reader shown only the log plot would conclude that the calcite does nothing, because everything it does happens where that axis has no room left.
for (lbl, Q, qd) in (("with 3.5 % calcite", Q_c, qd_c), ("no limestone", Q_n, qd_n))
j = argmax(qd)
@printf "%-20s Q: %5.1f (1 d) %5.1f (7 d) %5.1f (28 d) J/g peak %.2f mW/g at %.2f h\n" lbl (
Q[argmin(abs.(t_cal .- 86400))] / BINDER_G
) (Q[argmin(abs.(t_cal .- 7 * 86400))] / BINDER_G) (Q[end] / BINDER_G) (
qd[j] / BINDER_G * 1000
) (t_cal[j] / 3600)
endwith 3.5 % calcite Q: 178.8 (1 d) 354.1 (7 d) 428.5 (28 d) J/g peak 4.06 mW/g at 7.15 h
no limestone Q: 177.6 (1 d) 342.7 (7 d) 413.1 (28 d) J/g peak 4.23 mW/g at 7.15 hThe 28-day figures, about 420 J/g with limestone against 405 J/g without, are the ordinary range for a CEM I. The limestone raises the heat slightly rather than diluting it, because the carbonate is not inert here: it converts the aluminate to monocarboaluminate and stabilizes the ettringite (§4), and both reactions are exothermic. Substituting more limestone would eventually reverse the sign of that effect, which is the trade the LC³ literature is about.
The semi-adiabatic cell
A semi-adiabatic cell is not read off a trajectory like the curves above: the heat raises the temperature, the temperature raises the rates, and the cell belongs in the kinetic problem. A semi-adiabatic calorimeter, inside the kinetics runs this model in the NF EN 196-9 cell of (Lavergne et al., 2018) and compares the temperature with theirs.
6. Porosity, and what it is referred to
The porosity of a setting binder is not V_liquid / V_total: the denominator shrinks with the reactions, while a sealed specimen keeps the volume it was cast with, and the empty porosity left by the Le Chatelier contraction is not a species at all. ionic_phase_history returns the two-argument porosity, referred to the fresh paste and counting the chemical shrinkage as void.
@printf "%6s %8s %8s %8s\n" "t [d]" "φ" "S" "void"
for t in (0.25, 1.0, 3.0, 7.0, 28.0)
i = argmin(abs.(td .- t))
@printf "%6.2f %8.4f %8.4f %8.4f\n" td[i] poro[i].total (
poro[i].liquid / poro[i].total
) get(fracs[i], "void", 0.0)
end t [d] φ S void
0.25 0.5971 0.9483 0.0309
0.93 0.5135 0.9082 0.0471
2.89 0.4468 0.8715 0.0574
6.50 0.4092 0.8464 0.0629
28.00 0.3630 0.8091 0.0693A sealed paste desaturates as it hydrates though no water ever leaves it: the saturation falls from 0.94 to 0.81 while the chemical shrinkage grows to about 7 % of the fresh volume.
Where the mechanics goes
The volume fractions above are exactly what a micromechanical estimate of the elastic modulus needs. That extension — the same chemistry, feeding a four-scale self-consistent/Mori–Tanaka scheme, and the setting threshold that comes with it — is the chapter Hydration through the pore solution of MeanFieldHomogenization.jl, which duplicates the model of this page and adds the homogenization.
Where to go next
The kinetic parameters used here are published values; fitting them to a measured calorimetry curve, and what such a fit can and cannot identify, is Calibrating the hydration kinetics of a CEM I on measured calorimetry.