The AFm of a CEM I 52.5 N and its miscibility gap: three ways to declare a binary that can unmix
Before this page
The AFm and AFt phases of a Portland cement are binaries — sulfate against hydroxide, sulfate against carbonate — and the Redlich-Kister parameters (Redlich and Kister, 1948) published for them in Cemdata18 (Lothenbach et al., 2019) are strong enough that the mixing energy is concave over an interval. Where an energy is concave the Gibbs minimum is not one composition but two: the phase unmixes, and the equilibrium is a miscibility gap.
This page runs the same cement three ways on the AFm sulfate/hydroxide binary, because the three are genuinely different claims:
ideal mixing — an answer, certified, to a question that was changed;
the published parameters with one composition — the question kept; the answer is a minimum only if the composition falls outside the gap, and the certificate tests it;
the published parameters with two compositions — room for the pair, which only a composition inside the gap needs.
The theory is in the solid-solution chapter; this page is the measurement.
using ChemistryLab
using DynamicQuantities
using Logging
using OptimaSolver
using OrderedCollections
using Printf
using Plots
default(framestyle = :box, grid = false)
substances = build_species(datapath("cemdata18-thermofun.json"); verbose = false)
byname = Dict(symbol(s) => s for s in substances)
molar_mass_g(n) = ustrip(us"g/mol", byname[n][:M])1. Where the gap is, before any cement is involved
spinodal_interval scans the second derivative of the molar mixing energy and returns the interval over which it is negative, or nothing when the model is convex. It needs no solve and no system:
RT = R_GAS * 298.15
# Cemdata18's dimensionless Guggenheim parameters, from
# data/literature/Lothenbach2019.json: a₀ = α₀RT and a₁ = α₁RT.
function published_rk(pair)
p = literature_row("Lothenbach2019", "guggenheim_parameters", pair)
return RedlichKisterModel(a0 = p.alpha0 * RT, a1 = p.alpha1 * RT)
end
models = OrderedDict(
"ideal" => IdealSolidSolutionModel(),
"AFm SO4/OH (published)" => published_rk("AFm SO4/OH"),
"AFt SO4/CO3 (published)" => published_rk("AFt SO4/CO3"),
"regular, W = 1.9 RT" => RegularSolutionModel([0.0 1.9RT; 1.9RT 0.0]),
"regular, W = 2.1 RT" => RegularSolutionModel([0.0 2.1RT; 2.1RT 0.0]),
)
for (label, m) in models
gap = spinodal_interval(m, 2)
if gap === nothing
@printf(" %-26s convex everywhere\n", label)
else
@printf(" %-26s concave for x in [%.3f, %.3f]\n", label, gap[1], gap[2])
end
end ideal convex everywhere
AFm SO4/OH (published) concave for x in [0.631, 0.914]
AFt SO4/CO3 (published) concave for x in [0.564, 0.830]
regular, W = 1.9 RT convex everywhere
regular, W = 2.1 RT concave for x in [0.391, 0.609]The last two are the classical check: a symmetric regular solution unmixes above
The curve itself makes the shape of the problem visible. A convex energy has one minimum; a concave stretch means the straight line between two compositions lies below the curve, and the system takes the line — which is to say, it takes the two compositions and not any single one between them.
xs = range(0.001, 0.999; length = 400)
g(m, x) = begin
# molar Gibbs energy of mixing, in units of RT
ideal = x * log(x) + (1 - x) * log(1 - x)
excess = if m isa RedlichKisterModel
x * (1 - x) * (m.a0 + m.a1 * (2x - 1) + m.a2 * (2x - 1)^2) / RT
else
0.0
end
ideal + excess
end
afm = models["AFm SO4/OH (published)"]
gap = spinodal_interval(afm, 2)
fig = plot(xs, [g(afm, x) for x in xs]; label = "published Redlich-Kister",
color = :firebrick, linewidth = 2,
xlabel = "x (C4AH13, the first end-member)", ylabel = "g / RT",
title = "The mixing energy of the AFm sulfate/hydroxide binary",
size = (760, 420), bottom_margin = 8Plots.mm, left_margin = 8Plots.mm)
plot!(fig, xs, [g(IdealSolidSolutionModel(), x) for x in xs];
label = "ideal mixing", color = :steelblue, linewidth = 2, linestyle = :dash)
vspan!(fig, [gap[1], gap[2]]; color = :firebrick, alpha = 0.12, label = "spinodal")
2. The cement, declared three ways
A CEM I paste, with the AFm binary as the only phase whose model changes between the three runs:
# The Bogue composition of the CEM I 52.5 N of [Lavergne2018](@citet), Table 9,
# and its gypsum.
bogue_table = literature_table("Lavergne2018", "cement_bogue")
CLINKER = OrderedDict(zip(bogue_table.phase, bogue_table.percent ./ 100))
GYPSUM = literature_value("Lavergne2018", "gypsum_percent") / 100
WB = 0.50
BINDER_G = 100.0
pure = split(
"C3S C2S C3A C4AF Gp Anh Cal Portlandite ettringite monocarbonate " *
"hemicarbonate C3AH6 C3FH6 straetlingite FeOOHmic AlOHmic Amor-Sl Brc"
)
CSHQ = ["CSHQ-JenD", "CSHQ-JenH", "CSHQ-TobD", "CSHQ-TobH", "KSiOH", "NaSiOH"]
AFM = ["C4AH13", "monosulphate12"]
extra_aqueous = ["SO4-2", "CO2@"]
# 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). This paste carries no alkali, so neither set describes
# its calcium hydroxide and sulfate solution more closely; the KOH set is used,
# as on the other cement pages.
model = cemdata18_activity_model(:KOH)
"""
One system, differing only in how the AFm binary is declared.
`autostart` is a keyword because two of the three cases below are *meant* not to
certify. For those the multi-start cascade -- every backend, then the ideal
pre-solve, then the homotopy continuation -- cannot help: there is no admissible
single-composition minimum for it to find, so it spends minutes arriving at the
answer the first route already gave. The certificate reported is the same one
either way; declining the cascade declines only the search for a better start.
"""
function run_case(afm_phase; autostart = true)
sp = speciation(substances, vcat(pure, CSHQ, AFM, extra_aqueous);
aggregate_state = [AS_AQUEOUS])
ss = [SolidSolutionPhase("CSHQ", [byname[m] for m in CSHQ]), afm_phase]
cs = ChemicalSystem(sp, CEMDATA_PRIMARIES; solid_solutions = ss)
st = ChemicalState(cs)
for (phase, frac) in CLINKER
set_quantity!(st, phase,
BINDER_G * (1 - GYPSUM) * frac / molar_mass_g(phase) * u"mol")
end
set_quantity!(st, "Gp", BINDER_G * GYPSUM / molar_mass_g("Gp") * u"mol")
set_quantity!(st, "H2O@", BINDER_G * WB / molar_mass_g("H2O@") * u"mol")
b = budget(st)
# Every case prints its certificate; the warning of a refusal would only repeat it.
eq, cert = with_logger(NullLogger()) do
equilibrate_certified(st; model = model, b = b, autostart = autostart)
end
return cs, eq, cert
end
afm_em = [byname[m] for m in AFM]
published = published_rk("AFm SO4/OH")Case 1 — ideal mixing
cs1, eq1, c1 = run_case(SolidSolutionPhase("AFm_SO4_OH", afm_em))
@printf("optimal=%-5s worst SI=%+.2e balance=%.1e pH=%.4f\n",
c1.optimal, c1.worst_supersaturation, c1.balance, pH(eq1, model))optimal=true worst SI=-8.60e-14 balance=2.3e-12 pH=12.4768Case 2 — the published parameters, one composition
SolidSolutionPhase refuses this declaration outright, and names the interval:
err = try
SolidSolutionPhase("AFm_SO4_OH", afm_em; model = published)
nothing
catch e
e
end
println(err.msg)SolidSolutionPhase "AFm_SO4_OH": the mixing energy of this model is CONCAVE for x in [0.631, 0.914] at T = 298.15 K, so the phase unmixes there: the Gibbs minimum inside that interval is two coexisting compositions, not one, and this formulation has a single amount per species to describe it with. Declare it with `instances = 2` to give it two, which is what a miscibility gap needs; or use ideal mixing, or parameters that keep the energy convex; or pass `check_convexity = false` to proceed with one composition anyway — in which case the optimality certificate no longer proves anything, its sufficiency resting on convexity.Waiving the refusal is what a caller does who believes the answer stays outside the gap. The certificate decides it, testing the present phase against splitting as well as for stationarity:
cs2, eq2, c2 = run_case(SolidSolutionPhase("AFm_SO4_OH", afm_em;
model = published, check_convexity = false);
autostart = false)
@printf("optimal=%-5s worst SI=%+.2e balance=%.1e pH=%.4f\n",
c2.optimal, c2.worst_supersaturation, c2.balance, pH(eq2, model))
if hasproperty(c2, :worst_violation_split)
@printf("phases reported as wanting to split: %d worst split measure %+.2e\n",
length(c2.split_phases), c2.worst_violation_split)
endoptimal=true worst SI=-1.40e-13 balance=1.6e-14 pH=12.4768
phases reported as wanting to split: 0 worst split measure -1.40e-13Here the belief is right: the AFm of this paste has a C4AH13 fraction of 0.28, outside the pair section 3 computes, and the certificate finds no phase that wants to split. The answer is one composition, and a minimum. Until a phase with an excess term was inverted by Newton's method, this solve did not converge, and the tangent-plane test, applied to the point it stopped at, reported the phase as wanting to split: a failure of the search read as one of the chemistry.
Case 3 — the published parameters, two compositions
cs3, eq3, c3 = run_case(SolidSolutionPhase("AFm_SO4_OH", afm_em;
model = published, instances = 2);
autostart = false)
@printf("optimal=%-5s worst SI=%+.2e balance=%.1e pH=%.4f\n",
c3.optimal, c3.worst_supersaturation, c3.balance, pH(eq3, model))
n3 = ustrip.(us"mol", eq3.n)
for (grp, ph) in zip(cs3.ss_groups, cs3.solid_solutions)
startswith(name(ph), "AFm") || continue
tot = sum(n3[i] for i in grp)
x = tot > 1.0e-12 ? n3[grp[2]] / tot : NaN
@printf(" %-14s total %8.5f mol x(sulfate) = %s\n",
name(ph), tot, isnan(x) ? "—" : @sprintf("%.4f", x))
endoptimal=true worst SI=+3.36e-14 balance=4.2e-11 pH=12.4768
AFm_SO4_OH total 0.03695 mol x(sulfate) = 0.7231
AFm_SO4_OH#2 total 0.00000 mol x(sulfate) = —With a second instance offered, the paste splits nothing: both instances hold the composition of case 2, at the same pH, and the answer is certified. How the AFm is shared between them is not a result. Two instances of a phase outside its gap are degenerate: any division of one composition between them has the same Gibbs energy, so only their sum is determined, and the division printed is wherever the search stopped. That is why instances = :auto exists: it gives a phase its second instance only when the certificate of the one-composition answer asks for it.
3. The pair itself is computable, from the model alone
The two compositions are not unknown. They follow from the mixing model alone, by the construction Glynn and Reardon (1990) set out and that PHREEQC (Parkhurst and Appelo, 2013) uses for a binary solid solution: the pair at which a single straight line is tangent to
ct = common_tangent(published, 2)
sp = spinodal_interval(published, 2)
@printf("common tangent : x = %.4f and %.4f\n", ct[1], ct[2])
@printf("spinodal : x = %.4f to %.4f (contained in it, as it must be)\n",
sp[1], sp[2])common tangent : x = 0.4999 and 0.9700
spinodal : x = 0.6312 to 0.9142 (contained in it, as it must be)common_tangent solves two equations in two unknowns by Newton and costs microseconds — it involves no chemical system at all. Its answer is checkable without trusting the code: for a symmetric model
A = 3.0 # W/RT, well past the threshold of 2
sym = RegularSolutionModel([0.0 A*RT; A*RT 0.0])
a, b = common_tangent(sym, 2)
@printf("symmetric model: x = %.6f and %.6f, and b = 1 - a to %.1e\n",
a, b, abs(b - (1 - a)))
@printf("residual of ln(x/(1-x)) + A(1-2x) = 0 : %+.2e and %+.2e\n",
log(a / (1 - a)) + A * (1 - 2a), log(b / (1 - b)) + A * (1 - 2b))symmetric model: x = 0.070720 and 0.929280, and b = 1 - a to 0.0e+00
residual of ln(x/(1-x)) + A(1-2x) = 0 : -2.26e-13 and +2.27e-13So the endpoints of the gap are available to a user who needs them. What does not work is handing them to the minimization as a starting point and expecting it to keep the two lobes apart:
And the split follows, by mass balance
Inside a gap only the proportions move: the two compositions are the same whatever the overall composition is. So once the pair is known, how a given overall composition miscibility_split returns, together with the Gibbs energy the separation releases:
# The composition of the AFm of this paste, as the single-composition answer of
# case 2 returns it: x is the fraction of C4AH13, the first end-member.
n2 = ustrip.(us"mol", eq2.n)
grp2 = only(g for (g, ph) in zip(cs2.ss_groups, cs2.solid_solutions) if startswith(name(ph), "AFm"))
x̄_paste = n2[grp2[1]] / sum(n2[i] for i in grp2)
for x̄ in (0.20, 0.40, x̄_paste, 0.80, 0.96)
r = miscibility_split(published, x̄, 2)
if r.f_beta == 0
@printf("x̄ = %.4f : homogeneous (outside the pair)\n", x̄)
else
# Two calls rather than one: `@printf` takes a LITERAL format string,
# and `"a" * "b"` is an expression -- `ArgumentError: First argument to
# @printf after io must be a format string`, raised while the block is
# lowered, so it kills the build rather than the line.
@printf("x̄ = %.4f : %.1f %% at x=%.4f and %.1f %% at x=%.4f",
x̄, 100r.f_alpha, r.x_alpha, 100r.f_beta, r.x_beta)
@printf(" — releases %6.1f J/mol (mass balance %+.0e)\n",
r.Δg, r.f_alpha * r.x_alpha + r.f_beta * r.x_beta - x̄)
end
endx̄ = 0.2000 : homogeneous (outside the pair)
x̄ = 0.4000 : homogeneous (outside the pair)
x̄ = 0.2769 : homogeneous (outside the pair)
x̄ = 0.8000 : 36.2 % at x=0.4999 and 63.8 % at x=0.9700 — releases 102.8 J/mol (mass balance +0e+00)
x̄ = 0.9600 : 2.1 % at x=0.4999 and 97.9 % at x=0.9700 — releases 2.1 J/mol (mass balance +0e+00)Inside the gap
Everything above costs microseconds and needs no solver. What it does not give is instances = :auto gives the phase its second instance and the split passes look for the pair, with the aqueous solution iterated along; the theory chapter runs that on a binary its element budget holds in the gap (Solid solutions, section 6).
Where this leaves the three questions
| detect a gap | yes — the certificate refuses a single composition inside one, and names the phase |
| locate it | yes — spinodal_interval and common_tangent, from the model alone |
| represent it | yes — instances = 2, or instances = :auto, which adds the second when the certificate asks for it; the species exist, the groups stay disjoint, conservation is untouched |
| split a given overall composition | yes — miscibility_split, exact, with the energy it releases |
| iterate that back through the aqueous equilibrium | yes, with instances = :auto where the budget holds the composition inside the gap (Solid solutions, section 6); this paste does not need it |
The first four are what a user needs to answer "is this phase homogeneous, and if not, into what?". The last is the coupling with the rest of the paste.
Stated plainly because a reader deciding whether to trust a number here deserves to know which of these five they are relying on.
Why the second instance is refused for a convex model
Two instances of a convex phase are degenerate: every way of splitting the amount between them has the same energy, so the minimum is a flat manifold and the optimizer is asked to choose a point on it for no reason. Inside a spinodal the common-tangent pair is unique, and the degeneracy does not arise. SolidSolutionPhase therefore admits instances > 1 exactly where spinodal_interval reports a gap.
4. The three answers side by side
labels = ["ideal\nmixing", "published,\none composition", "published,\ntwo compositions"]
certs = [c1, c2, c3]
eqs = [eq1, eq2, eq3]
vols = [ustrip(uconvert(us"cm^3", volume(e).total)) for e in eqs]
phs = [pH(e, model) for e in eqs]
ok = [c.optimal for c in certs]
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", title = "Pore solution pH",
color = [o ? :seagreen : :firebrick for o in ok],
ylims = (12.0, 13.5))
fig3 = plot(p1, p2; layout = (1, 2), size = (900, 420),
bottom_margin = 16Plots.mm, left_margin = 8Plots.mm)
Green is a proof and red is its absence — not a claim that a red bar is far from the truth, which nothing here establishes. The three give the same pH to four figures. The middle one is proved with the published parameters kept, because this paste's AFm lies outside the gap; the third, which offers a second instance that nothing needs, is the one the search cannot close.
What this does and does not settle
optimal = true in case 2 is a proof of a KKT point at which no present phase wants to split and no absent phase is supersaturated. For a convex problem those conditions are sufficient for a global minimum, and that is why the certificate means what it means on every other page. With a concave mixing model the problem is no longer convex, so global optimality is not implied by stationarity alone — what the certificate adds here is the tangent-plane test, which is exactly the condition that separates a one-composition minimum from a spurious stationary point; its scope does not claim a global minimum.
Where to go next
A cement declared with every solid solution of its database, and the phase list treated as a modeling decision, is A CEM I at equilibrium, with every solid solution declared. The aqueous counterpart of this chapter continues with The Pitzer model, and what it can be used with.