The full Portland cement, through its pore solution
The hydrating paste, 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 against a certified 9.90, and it rarely reports convergence at all, so its return code cannot tell the two cases apart.
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)┌───────────────────────────────────────────────────┐
│ Loading database: data/cemdata18-thermofun.json │
└───────────────────────────────────────────────────┘
┌────────────────────┐
│ Building species │
└────────────────────┘
Progress: 86%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | ETA: 0:00:00[K
Progress: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:00[K
62 species: 49 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.
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, clinker = CLINKER, gypsum = 0.046, filler = 0.035, 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┌ Warning: Replacing docs for `Main.IONIC_DEFAULT_SYSTEM :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.ionic_water_group :: Tuple{Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.GEL_WATER_PER_CSH :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.build_ionic_system :: Union{Tuple{}, Tuple{Symbol}}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_CALIBRATION :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_INDUCTION_TAU :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_INDUCTION_M :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_INDUCTION_PHASES :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.ionic_induction :: Union{Tuple{}, Tuple{Any}, Tuple{Any, Any}}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_PK84 :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.ionic_reactions :: Tuple{Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.run_ionic_hydration :: Tuple{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.CALORIMETRY_MIX_C100 :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.CALORIMETRY_LOSS_A :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.CALORIMETRY_VESSEL_CP :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.sand_heat_capacity :: Tuple{Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.semiadiabatic_cell :: Tuple{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.langavant_temperature :: Tuple{Any, Any, Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.ionic_phase_history :: Tuple{Any, Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_DEFAULT_SYSTEM :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.ionic_water_group :: Tuple{Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.GEL_WATER_PER_CSH :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.build_ionic_system :: Union{Tuple{}, Tuple{Symbol}}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_CALIBRATION :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_INDUCTION_TAU :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_INDUCTION_M :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_INDUCTION_PHASES :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.ionic_induction :: Union{Tuple{}, Tuple{Any}, Tuple{Any, Any}}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.IONIC_PK84 :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.ionic_reactions :: Tuple{Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.run_ionic_hydration :: Tuple{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.CALORIMETRY_MIX_C100 :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.CALORIMETRY_LOSS_A :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.CALORIMETRY_VESSEL_CP :: Union{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.sand_heat_capacity :: Tuple{Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.semiadiabatic_cell :: Tuple{}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.langavant_temperature :: Tuple{Any, Any, Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌ Warning: Replacing docs for `Main.ionic_phase_history :: Tuple{Any, Any}` in module `Main`
└ @ Base.Docs docs/Docs.jl:253
┌───────────────────────────────────────────────────┐
│ Loading database: data/cemdata18-thermofun.json │
┌───────────────────────────────────────────────────┐
│ Loading database: data/cemdata18-thermofun.json │
└───────────────────────────────────────────────────┘
└───────────────────────────────────────────────────┘
┌────────────────────┐
│ Building species │
└────────────────────┘
┌────────────────────┐
│ Building species │
└────────────────────┘
Progress: 67%|███████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | ETA: 0:00:00[K
Progress: 71%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | ETA: 0:00:00[K
Progress: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:00[K
Progress: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:00[K
┌ Warning: Using arrays or dicts to store parameters of different types can hurt performance.
│ Consider using tuples instead.
└ @ SciMLBase ~/.julia/packages/SciMLBase/1DcFZ/src/performance_warnings.jl:33
┌ 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:398
┌ Warning: the dual equilibrium solve did not certify optimality; audit it with `optimality_certificate`.
└ @ ChemistryLab ~/work/ChemistryLab.jl/ChemistryLab.jl/src/equilibrium/dual_solver.jl:283
┌ Warning: 256 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 9.42e-6 mol ON THE ACCEPTED STEPS, i.e. on the trajectory itself; 9.42e-6 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 (Parrot-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:168
┌ Warning: 586 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 8.85e-6 mol ON THE ACCEPTED STEPS, i.e. on the trajectory itself; 8.85e-6 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 (Parrot-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:168
┌───────────────────────────────────────────────────┐
│ Loading database: data/cemdata18-thermofun.json │
└───────────────────────────────────────────────────┘
┌────────────────────┐
│ Building species │
└────────────────────┘
Progress: 81%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | ETA: 0:00:00[K
Progress: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:00[K
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.18.0, OptimaSolver version: 0.6.0
commit: 1fcbb4a5
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 conditionsThe 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.3385 0.0676
gypsum 0.0160 0.0000
calcite 0.0104 0.0000
C-S-H 0.0000 0.3515
CH 0.0000 0.1287
AFt 0.0098 0.0778
AFm 0.0000 0.1088
FH3 0.0000 0.0092
water 0.6129 0.1867
void 0.0125 0.0697
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.1088
no limestone AFt peak 0.0716 at 0.35 d AFt(28 d) 0.0000 AFm(28 d) 0.1027Two 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
The heat comes from the enthalpy of the whole system, Eqs. (17)–(21) of (Lavergne et al., 2018): for a quasi-static isobaric process the released heat balances the change of enthalpy, and the enthalpy is a sum over species of molar enthalpies of formation,
At fixed temperature this collapses to Q(t) = H(t₀) − H(t). Enthalpy is a state function, so reactants, ions and hydrates are each counted once and no reaction stoichiometry has to be written down — which matters here, because the hydrates are not produced by any reaction the model declares.
That is the whole difficulty of doing calorimetry on this model, and it is worth stating plainly. heat_rate sums rᵢ(−Δ_r H⁰ᵢ) over the kinetic reactions, which is right for a stoichiometric model whose reactions produce the hydrates directly. Here the kinetic reactions only dissolve the clinker into ions; the hydrates are precipitated by the Gibbs minimization, whose heat that sum cannot see. Driving a semi-adiabatic cell from it put the temperature rise at 207 K.
Nor can the enthalpy be read from the composition the integrator carries: under partial equilibrium that composition comes from an in-run, warm-started minimization which is not certified, and a single hydrate is worth hundreds of kilojoules. Read that way the curve came out at 12.7, 145, 1174, 936 and 631 J/g at 1 h, 6 h, 12 h, 1 d and 2 d — heat that rises and then falls, which no calorimeter has ever measured. heat_release therefore reads the certified speciations of §3.
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: 179.1 (1 d) 354.0 (7 d) 428.3 (28 d) J/g peak 3.89 mW/g at 7.15 h
no limestone Q: 177.6 (1 d) 342.7 (7 d) 413.1 (28 d) J/g peak 3.93 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 Langavant test (NF EN 196-9) lets the heat raise the temperature of the sample against the losses of the vessel. (Lavergne et al., 2018) write the loss as their Eq. (23),
and the numbers used below are theirs, for the plain-cement mix C100 of their Table 11 at w/b = 0.5:
| quantity | value | source |
|---|---|---|
| binder / dry sand / water | 371 g / 1113 g / 196 g | Table 11, C100 |
calorimeter vessel C_vessel | 380 J/K | §4.1 — see the note below |
| sand heat capacity | 812 J/K | Qtz of CEMDATA18, 0.73 J/(g·K) |
loss coefficient a | 75 J/(h·K) = 0.0208 W/K | Eq. (23), NF EN 196-9 calibration |
loss coefficient b | 0.260 J/(h·K²) = 7.22e-5 W/K² | Eq. (23) |
The sand takes no part in the chemistry; it is there, as the paper says, "to avoid large temperatures", and enters only through its heat capacity. The paste's own Σᵢ nᵢ C°_{p,i}(T) — about 900 J/K at 28 days — comes from the database at each instant, so it is not counted twice.
The vessel heat capacity is read as 380 J/K, not 380 kJ/K
The paper prints "about 380 kJ/K", and that cannot be the figure its own results correspond to. Its Table 11 mix holds 371 g of binder releasing some 420 J/g, i.e. about 156 kJ; against 380 kJ/K the temperature would rise by 0.4 K, where the test reports tens of kelvin. The rest of the setup is consistent with joules — sand and water alone contribute roughly 1.6 kJ/K — so 380 J/K puts the total near 2.1 kJ/K and the adiabatic rise near 75 K, which is the order the measurements show. It is read as 380 J/K here, and this note is deliberate: the alternative is to change a published number in silence.
# Computed alongside the run, in `precomputed.jl`: the cell temperature
# needs the heat capacity of the paste at each instant, so it needs the states
# themselves rather than the heat curve alone.
T_c = heat_c.columns["T_semiadiabatic_K"]
T_n = heat_n.columns["T_semiadiabatic_K"]
p_T = plot(;
xscale = :log10, xlabel = "time [days]", ylabel = "T − T_env [K]",
title = "Semi-adiabatic cell (NF EN 196-9)", legend = :topleft, size = (760, 420),
)
plot!(p_T, t_cal ./ 86400, T_c .- 293.15; lw = 2, color = 1, label = "with 3.5 % calcite")
plot!(p_T, t_cal ./ 86400, T_n .- 293.15; lw = 2, color = 2, ls = :dash, label = "no limestone")
# The states themselves, not just the columns: a heat capacity is a property of
# the whole paste, so it needs the speciation and not a plotted curve. They come
# from the same memoized run as the tables above.
states_c = coupled_states("ionic_opc")
states_n = coupled_states("ionic_nolimestone")
m_binder_g = ustrip(us"kg", CALORIMETRY_MIX_C100.binder) * 1000
C_fixed = CALORIMETRY_VESSEL_CP + sand_heat_capacity(CALORIMETRY_MIX_C100.sand)
for (lbl, T, Q, st) in (("with 3.5 % calcite", T_c, Q_c, states_c),
("no limestone", T_n, Q_n, states_n))
j = argmax(T)
C_tot = C_fixed + ustrip(us"J/K", heat_capacity(st[end])) * m_binder_g / 1000
@printf "%-20s ΔT max %.1f K at %.1f h adiabatic ΔT(28 d) %.1f K\n" lbl (
T[j] - 293.15
) (t_cal[j] / 3600) (Q[end] / BINDER_G * m_binder_g / C_tot)
endwith 3.5 % calcite ΔT max 19.8 K at 22.3 h adiabatic ΔT(28 d) 75.9 K
no limestone ΔT max 19.3 K at 22.3 h adiabatic ΔT(28 d) 72.9 KA rise of about 19 K at roughly one day, against an adiabatic 75 K: the sand and the losses absorb three quarters of the heat, which is what the test is designed to do.
One approximation, and it is in the direction you would expect
The heat rate above was computed at 20 °C. The temperature reached in the cell accelerates the reactions — Parrot–Killoh carries activation energies of 42, 21, 54 and 32 kJ/mol for C₃S, C₂S, C₃A and C₄AF — and that feedback is not included, so the true peak comes earlier and higher. Closing the loop needs the heat source inside the ODE, which under partial equilibrium requires differentiating the equilibrium map; KineticsProblem refuses that combination with a warning rather than returning a number it cannot support. For a stoichiometric model, where the reactions do produce the hydrates, the fully coupled version is scripts/opc_semiadiabatic_calorimetry.jl.
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.5975 0.9470 0.0317
0.93 0.5139 0.9071 0.0478
2.89 0.4473 0.8705 0.0579
6.50 0.4097 0.8456 0.0633
28.00 0.3635 0.8083 0.0697A 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.