A CEM I at equilibrium, with every solid solution declared
Before this page
A thermodynamic model of a cement is only as good as the list of phases it is allowed to form. Leave one out and the calculation cannot report it; put one in that does not belong and it can take mass that belongs elsewhere. The list is therefore a modeling decision, and the honest way to make it is to declare everything the database defines and let the minimization decide.
This page does exactly that on a CEM I at w/c = 0.5, from the oxide analysis on the datasheet through the Bogue conversion to the anhydrous charge, and then to the hydrated assemblage. It also measures what declaring only the two or three solid solutions one expects to find changes: on this cement, the same answer, with a certificate that covers fewer phases.
using ChemistryLab
using DynamicQuantities
using OptimaSolver
using OrderedCollections
using Printf
using Plots
# Frame style only — never the font. Documenter runs all `@example` blocks in
# ONE process, so `default(fontfamily = ...)` here would apply to every page
# built after this one. The font is set once, for the whole site, in
# `docs/make.jl`, and the reason is recorded there.
default(framestyle = :box, grid = false)1. The datasheet, and the Bogue conversion
The input is what a cement works publishes: an oxide analysis, in grams per hundred grams. Nothing else about the cement is assumed. The analysis below is a representative CEM I; where it was first taken from is not recorded, so it stands for a class of cement rather than for a product.
# A representative CEM I analysis, g per 100 g; its source is not recorded.
oxides = OrderedDict(
"CaO" => 65.03, "SiO2" => 21.4, "Al2O3" => 3.84, "Fe2O3" => 4.49,
"MgO" => 1.0, "K2O" => 0.46, "Na2O" => 0.13, "SO3" => 2.3, "CO2" => 0.0,
)
water_g = 50.0 # w/c = 0.5
# At w/c = 0.5 the whole clinker can react, and that is computed rather than
# assumed: `powers_alpha_max(0.5)` is 1, because complete hydration needs about
# 0.42 g of water per gram of cement and there is more than that here. Below
# 0.42 it would not be 1, the paste would arrest with water still in it, and
# reacting the whole clinker would be the wrong question -- which is the subject
# of [the w/c example](@ref sec-wc-ratio).
@assert powers_alpha_max(water_g / 100) == 1.0
sum(values(oxides)) # the analysis does not close at 100 g98.64999999999999It sums to 98.65 g, and the missing 1.35 g is loss on ignition and minor oxides the sheet does not report. Scaling the analysis to 100 g is a choice, and it moves every element by 1.37 %; it is made here so that "100 g of cement" means 100 g of the oxides that are modeled.
The database is loaded now rather than later, because the molar masses used below are read from it. None is written down here: an oxide formula goes through Species, which computes the mass from the formula and the element data, while a cement phase is a database species read by name — Species("C3S") would read that name as the formula C₃S, three carbons and a sulfur.
substances = build_species(datapath("cemdata18-thermofun.json"); verbose = false)
byname = Dict(symbol(s) => s for s in substances)
oxide_mass(ox) = ustrip(us"g/mol", Species(ox)[:M])
phase_mass(name) = ustrip(us"g/mol", byname[name][:M])The Bogue conversion turns those four main oxides into the four clinker phases. It is not reimplemented here: Bogue Calculation obtains it by inverting the mass stoichiometric matrix of the clinker phases over the oxide components, which is built from the species data, so no molar mass and no tabulated coefficient appears anywhere.
clinker = CemSpecies.(split("C3S C2S C3A C4AF"))
oxide4 = CemSpecies.(split("C S A F"))
bogue_inverse = inv(mass_matrix(CanonicalStoichMatrix(clinker)).A)
for (i, sp) in enumerate(clinker)
@printf(" %-5s = %+8.4f CaO %+8.4f SiO2 %+8.4f Al2O3 %+8.4f Fe2O3\n",
symbol(sp), bogue_inverse[i, 1], bogue_inverse[i, 2], bogue_inverse[i, 3], bogue_inverse[i, 4])
end C3S = +4.0714 CaO -7.6000 SiO2 -6.7177 Al2O3 -1.4298 Fe2O3
C2S = -3.0714 CaO +8.6000 SiO2 +5.0678 Al2O3 +1.0786 Fe2O3
C3A = +0.0000 CaO +0.0000 SiO2 +2.6500 Al2O3 -1.6920 Fe2O3
C4AF = +0.0000 CaO +0.0000 SiO2 +0.0000 Al2O3 +3.0432 Fe2O3Those are the classical Bogue coefficients, recovered rather than quoted. The lime combined with the sulfate is not available to the silicates, so it comes off before the inversion:
f = 100 / sum(values(oxides))
CaO_free = oxides["CaO"] -
oxides["SO3"] * oxide_mass("CaO") / oxide_mass("SO3") -
oxides["CO2"] * oxide_mass("CaO") / oxide_mass("CO2")
grams = bogue_inverse * [CaO_free, oxides["SiO2"], oxides["Al2O3"], oxides["Fe2O3"]]
batch = OrderedDict{String, Float64}(
symbol(sp) => grams[i] * f / phase_mass(symbol(sp))
for (i, sp) in enumerate(clinker)
)
# CEMDATA18 has no periclase and no anhydrous SO3, so the magnesia enters as
# brucite and the sulfate as gypsum. The water they bring is part of the cement
# and comes off the mixing water.
n_ox = Dict(k => v * f / oxide_mass(k) for (k, v) in oxides)
batch["Gp"] = n_ox["SO3"]
batch["Brc"] = n_ox["MgO"]
batch["K2O"] = n_ox["K2O"]
batch["Na2O"] = n_ox["Na2O"]
free_water = water_g / oxide_mass("H2O") - 2n_ox["SO3"] - n_ox["MgO"]
for (k, v) in batch
@printf(" %-5s %10.6f mol (%6.2f g)\n", k, v, v * phase_mass(k))
end
@printf(" %-5s %10.6f mol\n", "H2O@", free_water) C3S 0.281272 mol ( 64.22 g)
C2S 0.079776 mol ( 13.74 g)
C3A 0.009675 mol ( 2.61 g)
C4AF 0.028502 mol ( 13.85 g)
Gp 0.029123 mol ( 5.01 g)
Brc 0.025151 mol ( 1.47 g)
K2O 0.004950 mol ( 0.47 g)
Na2O 0.002126 mol ( 0.13 g)
H2O@ 2.692069 molnames_c = collect(keys(batch))
mass_c = [batch[k] * phase_mass(k) for k in names_c]
bar(names_c, mass_c;
legend = false, ylabel = "g per 100 g of oxides",
title = "The anhydrous charge, from the oxide analysis",
color = :steelblue, size = (720, 380), bottom_margin = 6Plots.mm,
left_margin = 6Plots.mm)
2. The phase list, and how many solid solutions there really are
The pure phases are the CEMDATA18 list a Portland cement is given. The solid solutions are every multi-end-member phase the database defines — not a selection.
pure = split(
"AlOHmic Kln Gr C12A7 C2S C3A C3S C4AF CA CA2 C2AH7.5 C3AH6 CAH10 " *
"monosulphate10.5 monosulphate12 monosulphate14 monosulphate16 " *
"monosulphate9 chabazite zeoliteP_Ca straetlingite5.5 monocarbonate9 " *
"hemicarbonat10.5 hemicarbonate hemicarbonate9 monocarbonate " *
"ettringite13 ettringite9 Arg Cal C3FH6 C4FH13 C3FS1.34H3.32 " *
"Fe-hemicarbonate Femonocarbonate Dis-Dol Ord-Dol Lim Portlandite " *
"Anh Gp hemihydrate Fe Sd Mag FeOOHmic Py Tro Melanterite K2SO4 " *
"syngenite K2O hydrotalcite Mgs Brc Na2SO4 natrolite zeoliteX " *
"zeoliteY Na2O Sulfur Amor-Sl"
)
solutions = [
"CSHQ" => ["CSHQ-JenD", "CSHQ-JenH", "CSHQ-TobD", "CSHQ-TobH",
"KSiOH", "NaSiOH"],
"C3(AF)S0.84H" => ["C3AFS0.84H4.32", "C3FS0.84H4.32"],
"AFt_SO4" => ["ettringite", "ettringite30"],
"AFt_SO4_CO3" => ["tricarboalu03", "ettringite03_ss"],
"AFm_SO4_OH" => ["C4AH13", "monosulphate12"],
"straetlingite" => ["straetlingite", "straetlingite7"],
"hydrotalc-pyro" => ["Mg3AlC0.5OH", "Mg3FeC0.5OH"],
"MSH" => ["M075SH", "M15SH"],
]
length(solutions), length(reduce(vcat, last.(solutions)))(8, 20)Eight phases over twenty end-members. A reader coming from GEM-Selektor (Kulik et al., 2013) will count ten phases in its CEMDATA18 setup, and the difference is worth stating because it is a difference of formulation and not of chemistry.
Why eight here and ten there
Two of the ten repeat their end-members with the roles exchanged: SO4_OH_AFm and OH_SO4_AFm are both C4AH13 + monosulphate12, and SO4_CO3_AFt and CO3_SO4_AFt are both tricarboalu03 + ettringite03_ss. Declaring the same binary twice is how that formulation opens a miscibility gap: two immiscible compositions of one binary can then coexist as two phases. It requires a species to belong to two phases at once.
Here the composition vector has one entry per species, so a species belongs to exactly one phase, and a gap is expressed by declaring the phase with instances = 2, which gives it a second copy of each end member (see The AFm of a CEM I 52.5 N and its miscibility gap). The mixing models of this page are ideal and open no gap, so the eight distinct phases are the whole of the chemistry.
Note also that AFt_SO4_CO3 describes the same aluminate sulfate as AFt_SO4 in a different normalization: ettringite03_ss is ettringite divided by three. Both are declared, because deciding in advance which normalization the answer will use is exactly the kind of prior knowledge this page is trying to do without.
members = reduce(vcat, last.(solutions))
species = speciation(
substances, vcat(pure, members, ["CO2@"]); aggregate_state = [AS_AQUEOUS]
)
cs = ChemicalSystem(
species, CEMDATA_PRIMARIES;
solid_solutions = [
SolidSolutionPhase(nm, [byname[m] for m in mem]) for (nm, mem) in solutions
],
)
length(cs.species), length(cs.solid_solutions)(150, 8)3. The solve
The composition is the anhydrous charge and the mixing water; the element vector is fixed once from it and passed to the solve, so whatever starting point is used can only decide whether the optimum is reached, never which optimum it is.
state = ChemicalState(cs)
for (sym, n) in batch
n > 0 && set_quantity!(state, sym, n * u"mol")
end
set_quantity!(state, "H2O@", free_water * u"mol")
set_quantity!(state, "CO2@", 1.0e-9u"mol")
b = budget(state)
# The activity model Cemdata18 prescribes (its Eq. C.1): extended Debye-Hückel,
# with the common ion size and B-dot the paper gives for KOH solutions (it also
# gives them for NaOH). The alkalis of this cement are mostly potassium.
model = cemdata18_activity_model(:KOH)
# The molar K/Na ratio of its alkalis, which is why the KOH set applies.
Mox(ox) = ustrip(us"g/mol", Species(ox)[:M])
println("molar K/Na of the alkalis: ", round((2 * oxides["K2O"] / Mox("K2O")) / (2 * oxides["Na2O"] / Mox("Na2O")); digits = 1))
eq, cert = equilibrate_certified(state; model = model, b = b)
@printf("certificate: optimal=%s worst supersaturation=%.3e balance=%.1e\n",
cert.optimal, cert.worst_supersaturation, cert.balance)
@printf("pH = %.4f total volume = %.3f cm3\n",
pH(eq, model), ustrip(uconvert(us"cm^3", volume(eq).total)))molar K/Na of the alkalis: 2.3
certificate: optimal=true worst supersaturation=-1.712e-13 balance=1.2e-14
pH = 13.1156 total volume = 74.195 cm3optimal = true is a proof, not a report of successful iteration: the Gibbs minimization is convex when the mixing terms are, so the KKT conditions are sufficient and optimality_certificate decides. A negative worst supersaturation says every absent phase is undersaturated at the answer — which is the statement that the phase list was searched, not merely used.
n = ustrip.(us"mol", eq.n)
present = [
(symbol(cs.species[i]), n[i]) for i in cs.idx_crystal if n[i] > 1.0e-8
]
sort!(present; by = last, rev = true)
for (s, x) in present
@printf(" %-18s %10.6f mol\n", s, x)
end Portlandite 0.453629 mol
CSHQ-JenD 0.170159 mol
CSHQ-TobD 0.130473 mol
CSHQ-JenH 0.107997 mol
C3AFS0.84H4.32 0.044371 mol
KSiOH 0.015343 mol
NaSiOH 0.007822 mol
ettringite 0.007743 mol
C3FS0.84H4.32 0.006317 mol
hydrotalcite 0.006288 mol
CSHQ-TobH 0.005419 mol
ettringite30 0.001961 molbar([s for (s, _) in present], [x for (_, x) in present];
legend = false, ylabel = "mol per 100 g of oxides", yscale = :log10,
title = "The hydrated assemblage", color = :seagreen,
xrotation = 45, size = (760, 420),
bottom_margin = 14Plots.mm, left_margin = 8Plots.mm)
4. Which solid solutions the calculation keeps, and why
Three of the eight carry mass. The other five are absent — and the point of this page is that the calculation says so, rather than being told.
The criterion is exact for an ideal solid solution and worth deriving, because it is what makes "declare everything" a sound procedure rather than a hopeful one. Let the phase have end-members
where
supersaturated above and undersaturated below. Two things follow. The first is that
saturation_indices reports the index of an end-member in its phase, which is
si = saturation_indices(eq, model)
lna = log_activities(eq, model)
inv_ln10 = inv(log(10))
Ω = OrderedDict(
name(ss) => sum(
10^(si[symbol(cs.species[i])] + lna[symbol(cs.species[i])] * inv_ln10)
for i in grp
)
for (grp, ss) in zip(cs.ss_groups, cs.solid_solutions)
)
for (nm, w) in Ω
@printf(" %-16s log10 Omega = %+8.3f %s\n", nm, log10(w),
abs(log10(w)) < 1.0e-6 ? "saturated, present" :
log10(w) > 0 ? "supersaturated" : "undersaturated")
end CSHQ log10 Omega = +0.000 saturated, present
C3(AF)S0.84H log10 Omega = +0.000 saturated, present
AFt_SO4 log10 Omega = -0.000 saturated, present
AFt_SO4_CO3 log10 Omega = -0.032 undersaturated
AFm_SO4_OH log10 Omega = -1.236 undersaturated
straetlingite log10 Omega = -2.883 undersaturated
hydrotalc-pyro log10 Omega = -9.843 undersaturated
MSH log10 Omega = -3.510 undersaturatedThe three present phases come out at exactly zero, which is the check on the criterion itself: a phase that the minimization put in must be saturated, and this formula, derived independently of the solve, says it is.
nms = collect(keys(Ω))
vals = [log10(w) for w in values(Ω)]
bar(nms, vals;
legend = false, ylabel = "log₁₀ Ω",
title = "How far each solid solution is from forming",
color = [v > -1.0e-6 ? :seagreen : :steelblue for v in vals],
xrotation = 30, size = (760, 420),
bottom_margin = 12Plots.mm, left_margin = 8Plots.mm)
hline!([0.0]; color = :black, linewidth = 1.5, label = "")
One number on that chart deserves attention. AFt_SO4_CO3 sits at −0.03: thirty thousandths of a log unit from forming. A phase list chosen by expectation would have dropped it as "obviously absent", and on a cement only slightly richer in carbonate it would not be. That is the argument for declaring everything, in one number.
5. What declaring only some of them costs
Declaring only the solid solutions one expects to find, and leaving the rest out, can be done in two ways. On this cement both reproduce the answer, and neither proves what the complete declaration proves.
Leaving the other end-members in the system as pure phases. Their saturation is then still tested, each end-member as a pure solid:
three = solutions[1:3]
cs_pure = ChemicalSystem(
species, CEMDATA_PRIMARIES; # the same species list
solid_solutions = [
SolidSolutionPhase(nm, [byname[m] for m in mem]) for (nm, mem) in three
],
)
st = ChemicalState(cs_pure)
for (sym, x) in batch
x > 0 && set_quantity!(st, sym, x * u"mol")
end
set_quantity!(st, "H2O@", free_water * u"mol")
set_quantity!(st, "CO2@", 1.0e-9u"mol")
b2 = budget(st)
eq2, cert2 = equilibrate_certified(st; model = model, b = b2)
@printf("optimal=%s worst supersaturation=%+.3f pH=%.4f V=%.3f cm3\n",
cert2.optimal, cert2.worst_supersaturation, pH(eq2, model),
ustrip(uconvert(us"cm^3", volume(eq2).total)))optimal=true worst supersaturation=-0.000 pH=13.1156 V=74.195 cm3The answer is that of the eight-phase run, to the printed digit, and it is certified. The assemblage shows why nothing else could have happened:
n2 = ustrip.(us"mol", eq2.n)
bad = [(symbol(cs_pure.species[i]), n2[i])
for i in cs_pure.idx_crystal if n2[i] > 1.0e-8]
sort!(bad; by = last, rev = true)
for (s, x) in bad
@printf(" %-18s %10.6f mol\n", s, x)
end Portlandite 0.453629 mol
CSHQ-JenD 0.170159 mol
CSHQ-TobD 0.130473 mol
CSHQ-JenH 0.107997 mol
C3AFS0.84H4.32 0.044371 mol
KSiOH 0.015343 mol
NaSiOH 0.007822 mol
ettringite 0.007743 mol
C3FS0.84H4.32 0.006317 mol
hydrotalcite 0.006288 mol
CSHQ-TobH 0.005419 mol
ettringite30 0.001961 molNo end-member offered as a pure solid has taken any mass, and none could have. As a pure phase an end-member has unit activity, whereas inside its solution its activity is its mole fraction, below one; its chemical potential as a pure solid is therefore never lower than inside the solution, and the pure phase can form only where the solution itself would consist of that end-member alone. Declaring end-members as pure phases is thus harmless, but it is not a test of the solution they belong to. Each pure end-member may be undersaturated while their mixture is not, since the entropy of mixing lowers the Gibbs energy of the solution below the composition-weighted mean of its end-members, and this certificate tests the five omitted solid solutions only through their end-members.
Keeping the other end-members out of the system entirely. This is what a hand-written species list usually does, and it reproduces the answer as well:
species3 = speciation(
substances, vcat(pure, reduce(vcat, last.(three)), ["CO2@"]);
aggregate_state = [AS_AQUEOUS],
)
cs3 = ChemicalSystem(
species3, CEMDATA_PRIMARIES;
solid_solutions = [
SolidSolutionPhase(nm, [byname[m] for m in mem]) for (nm, mem) in three
],
)
st3 = ChemicalState(cs3)
for (sym, x) in batch
x > 0 && set_quantity!(st3, sym, x * u"mol")
end
set_quantity!(st3, "H2O@", free_water * u"mol")
set_quantity!(st3, "CO2@", 1.0e-9u"mol")
b3 = budget(st3)
eq3, cert3 = equilibrate_certified(st3; model = model, b = b3)
@printf("optimal=%s worst supersaturation=%+.3f pH=%.4f V=%.3f cm3\n",
cert3.optimal, cert3.worst_supersaturation, pH(eq3, model),
ustrip(uconvert(us"cm^3", volume(eq3).total)))optimal=true worst supersaturation=-0.000 pH=13.1156 V=74.195 cm3Same pH, same volume, same assemblage as the eight-phase run, and a certificate again — but a weaker statement, since the five solid solutions it never declared are not among the phases it tested, not even through their end-members. The eight-phase certificate, whose margin on AFt_SO4_CO3 is only 0.03, is the one that actually proves the five are out. A certificate over a small phase list says nothing about the phases left out of it, however comfortable its margin.
labels = ["all eight\ndeclared", "three declared,\nothers pure",
"three declared,\nothers absent"]
vols = [ustrip(uconvert(us"cm^3", volume(e).total)) for e in (eq, eq2, eq3)]
phs = [pH(e, model) for e in (eq, eq2, eq3)]
ok = [cert.optimal, cert2.optimal, cert3.optimal]
p1 = bar(labels, vols; legend = false, ylabel = "total volume (cm³)",
color = [o ? :seagreen : :firebrick for o in ok], title = "Volume")
p2 = bar(labels, phs; legend = false, ylabel = "pH",
color = [o ? :seagreen : :firebrick for o in ok], title = "pH",
ylims = (12.0, 13.4))
plot(p1, p2; layout = (1, 2), size = (860, 400),
bottom_margin = 12Plots.mm, left_margin = 8Plots.mm)
Green is a certified answer, red is not.
6. A second opinion
Reaktoro solves the same problem from the same literature (Leal et al., 2017) and reads the same ThermoFun file this package ships, which makes it a useful independent check — see Validation against Reaktoro for the small case where every knob is matched and the agreement is quantified species by species.
On this cement the cross-check was run with scripts/crosscheck/cem1_solid_solutions_reaktoro.py, over the same database, the same species, the same ideal mixing and the same activity model (the extended Debye-Hückel law of Cemdata18, with its KOH parameters), and — the point of the exercise — the same element vector, transferred as moles rather than as grams of oxide so that the two codes do not each convert the analysis with their own atomic masses. The element vectors agree to eight decimal places, so what follows is a comparison of chemistry.
| ChemistryLab | Reaktoro | |
|---|---|---|
| solid solutions declared | 8 | 3 |
| pH | 13.1156 | 13.1426 |
| total volume | 74.195 cm³ | 74.239 cm³ |
The two agree to 0.027 units of pH and 0.06 % of volume, on the same twelve solid phases, with the majors within a few percent.
Two observations about the run itself:
The iteration count. A full cement assemblage takes several hundred iterations, more than the default limit of the Reaktoro run; with
EquilibriumOptions.optima.maxitersraised, it converges in 306 iterations.The eight-phase configuration needs a regularization. A solid solution whose end-members are all at zero has no mole fractions, so its ideal-mixing term is undefined there and its gradient depends on the direction of approach, a cone singularity at the origin of that phase's subspace. Here it is handled by the
ϵfloor of the activity closure together with the multi-start ofequilibrate_certified. The Reaktoro run uses the three-phase configuration documented for this system, hence the 8 and 3 of the table; on this cement the two reach the same twelve solid phases.
7. What this does not settle
The mixing is ideal, and that is an assumption. All eight phases are declared with ideal mixing here. CEMDATA18 publishes non-ideal parameters for the AFm and AFt binaries (Table 1 of Lothenbach et al. (2019), in data/literature/Lothenbach2019.json), and they are not used on this page. A non-ideal phase changes
The five absent phases are absent on this cement. AFt_SO4_CO3 at −0.03 is the warning: a different sulfate or carbonate content moves it across. The procedure transfers; the result does not.
No kinetics. This is the assemblage the cement would reach given unlimited time and complete reaction. What a paste actually reaches, and why it stops short, is The silicates of a CEM I clinker, hydrating end to end and Self-desiccation.
Where to go next
The equilibrium view continues with the water content in Effect of the water/cement ratio on the hydration of a CEM I and with carbonation in Carbonation of a CEM I paste. The binders that replace part of the clinker begin with CEM II/A-LL and CEM II/B-S.