Simplified dissolution of a CEM I clinker
Before this page
The tutorial Chemical Equilibrium and Cement chemist notation.
This example illustrates ChemistryLab's equilibrium workflow applied to cement chemistry: computing the hydration state of a simplified Portland clinker at a given water-to-cement ratio. ChemistryLab computes thermodynamic equilibrium by minimizing the Gibbs free energy of the system subject to element-conservation constraints. The workflow always follows the same four steps:
Build a
ChemicalSystem(species + stoichiometric matrix).Create an initial
ChemicalState(temperature, pressure, initial amounts).Call
equilibrate(or useEquilibriumSolverexplicitly).Inspect the resulting
ChemicalState.
Minimal workflow
The convenience function equilibrate handles everything with sensible defaults. The example below reproduces a simplified clinker dissolution calculation.
using ChemistryLab
using DynamicQuantities
substances = build_species(datapath("cemdata18-thermofun.json"))
input_species = split("C3S C2S C3A C4AF Gp Anh Portlandite Jennite H2O@ ettringite monosulphate12 C3AH6 C3FH6 C4FH13")
species = speciation(substances, input_species; aggregate_state=[AS_AQUEOUS])
cs = ChemicalSystem(species, CEMDATA_PRIMARIES)state = ChemicalState(cs)
# Clinker + gypsum composition, mass fractions summing to 0.984: an assumed
# CEM I composition, typical of the class; its source is not recorded.
compo = ["C3S" => 0.678, "C2S" => 0.166, "C3A" => 0.04, "C4AF" => 0.072, "Gp" => 0.028]
c = sum(last.(compo))
wc = 0.4 # water-to-cement ratio
w = wc * c
mtot = c + w
for x in compo
set_quantity!(state, x.first, x.second / mtot * u"kg")
end
set_quantity!(state, "H2O@", w / mtot * u"kg")
# pH-neutral seed for H⁺ and OH⁻
V = volume(state)
set_quantity!(state, "H+", 1e-7u"mol/L" * V.liquid)
set_quantity!(state, "OH-", 1e-7u"mol/L" * V.liquid)
state_eq = equilibrate(state)The solved state in full — every species, with its amount
state_eqChemicalState{Species, AbstractReaction, Quantity{Float64, SymbolicDimensions{FRInt32}}, Float64}
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ T : 298.15 K │
│ P : 1.0 bar │
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│ # liquid #│ n [mol]│ m [g]│ V [cm³]│ c [mol/L]│
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ tot. liquid│ 4.56714│ 82.3078│ 82.4378│ │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ H2O@│ 4.5639│ 82.2186│ 82.4609│ 55.3617│
│ OH-│ 0.00203441│ 0.0345992│ -0.00957768│ 0.0246781│
│ Ca+2│ 0.000846972│ 0.033945│ -0.0156171│ 0.0102741│
│ CaOH+│ 0.000347807│ 0.0198546│ 0.00200423│ 0.00421903│
│ AlO2-│ 7.2335e-6│ 0.000426629│ 6.8482e-5│ 8.77449e-5│
│ CaSiO3@│ 3.0664e-6│ 0.000356193│ 0.0│ 3.71965e-5│
│ HSiO3-│ 6.4551e-8│ 4.97624e-6│ 2.92149e-7│ 7.83027e-7│
│ Ca(SO4)@│ 3.0393e-8│ 4.13752e-6│ 1.42866e-7│ 3.68677e-7│
│ SO4-2│ 1.47857e-8│ 1.42025e-6│ 1.90996e-7│ 1.79355e-7│
│ Ca(HSiO3)+│ 1.05399e-8│ 1.23493e-6│ -7.10022e-8│ 1.27852e-7│
│ FeO2-│ 9.10608e-9│ 7.99906e-7│ 4.11187e-9│ 1.1046e-7│
│ SiO3-2│ 7.47692e-9│ 5.68859e-7│ 0.0│ 9.06977e-8│
│ SiO2@│ 1.68406e-10│ 1.01184e-8│ 2.70516e-9│ 2.04283e-9│
│ H2@│ 1.26873e-11│ 2.55775e-11│ 3.20535e-10│ 1.53901e-10│
│ FeOH+│ 8.69277e-12│ 6.33286e-10│ -1.45281e-10│ 1.05446e-10│
│ AlO2H@│ 8.16864e-12│ 4.90017e-10│ 1.06269e-10│ 9.90885e-11│
│ FeO2H@│ 4.03028e-12│ 3.58095e-10│ 2.90556e-11│ 4.88888e-11│
│ AlSiO5-3│ 2.22886e-12│ 3.01033e-10│ -7.60747e-11│ 2.70368e-11│
│ HS-│ 4.41466e-14│ 1.45984e-12│ 8.9218e-13│ 5.35513e-13│
│ H+│ 3.31947e-14│ 3.34602e-14│ 0.0│ 4.02663e-13│
│ Fe+2│ 1.11061e-14│ 6.20218e-13│ -2.51446e-13│ 1.3472e-13│
│ FeO+│ 1.26443e-17│ 9.0842e-16│ -5.31324e-16│ 1.5338e-16│
│ AlO+│2.2752800000000003e-18│ 9.77926e-17│ 6.99274e-19│ 2.75999e-17│
│ SO3-2│1.1683100000000001e-18│ 9.35313e-17│ -4.80874e-18│ 1.4172e-17│
│ HSO4-│ 5.80285e-19│ 5.63248e-17│ 2.02178e-17│7.039060000000001e-18│
│ Fe(SO4)@│ 3.55104e-19│ 5.39406e-17│ 5.93689e-19│ 4.30753e-18│
│ H2S@│1.7417000000000001e-19│5.935020000000001e-18│6.0874600000000004e-18│ 2.11274e-18│
│ S-2│ 1.09345e-20│ 3.5056e-19│1.2816899999999999e-25│ 1.32639e-19│
│ Si4O10-4│ 2.72466e-23│ 7.42008e-21│ 0.0│ 3.30511e-22│
│ HSO3-│7.828020000000001e-24│ 6.34578e-22│ 2.57986e-22│9.495660000000001e-23│
│ AlOH+2│ 3.98071e-25│ 1.75105e-23│ -1.08574e-24│ 4.82874e-24│
│ FeOH+2│ 1.5418e-26│1.1232300000000001e-24│ -3.90699e-25│ 1.87026e-25│
│ S2O3-2│6.764689999999999e-29│7.584369999999999e-27│1.8665299999999998e-27│ 8.20581e-28│
│ Fe(HSO4)+│ 9.42834e-31│1.4416799999999998e-28│ 1.77314e-29│ 1.14369e-29│
│ AlHSiO3+2│ 2.87448e-31│ 2.99152e-29│ 9.9485e-36│ 3.48685e-30│
│ Al+3│ 1.45756e-32│ 3.93273e-31│-6.594419999999999e-31│1.7680699999999998e-31│
│ FeHSiO3+2│ 3.79568e-33│ 5.04579e-31│1.2190299999999999e-37│4.6042899999999996e-32│
│ Al(SO4)+│ 2.08237e-35│ 2.5621e-33│ -1.25339e-34│2.5259899999999996e-34│
│ Fe+3│ 9.64621e-37│ 5.38693e-35│ -3.64533e-35│ 1.17012e-35│
│ Fe(SO4)+│ 1.90225e-39│ 2.88953e-37│ -5.01325e-39│ 2.30749e-38│
│ Al(SO4)2-│3.744479999999999e-40│ 8.20391e-38│ 1.16498e-38│4.542179999999999e-39│
│ Fe(SO4)2-│7.484279999999999e-45│ 1.85578e-42│ 2.2816e-43│ 9.07869e-44│
│ Fe2(OH)2+4│ 7.77913e-50│1.1334500000000001e-47│ -4.45016e-53│ 9.43636e-49│
│ Fe(HSO4)+2│ 2.05766e-51│3.1463499999999998e-49│4.7734999999999996e-51│ 2.49602e-50│
│ Fe3(OH)4+5│ 2.49749e-63│ 5.88317e-61│ -2.77411e-66│ 3.02955e-62│
│ O2@│ 2.18289e-74│ 6.9848e-73│ 6.65799e-73│ 2.64792e-73│
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│ # solid #│ n [mol]│ m [g]│ V [cm³]│ │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ tot. solid│ 6.0689│ 917.692│ 387.255│ │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ Portlandite│ 2.89109│ 214.207│ 95.5796│ │
│ Jennite│ 2.85524│ 546.425│ 223.851│ │
│monosulphate12│ 0.118057│ 73.4912│ 36.6095│ │
│ C3FH6│ 0.107551│ 46.893│ 16.7012│ │
│ C3AH6│ 0.0969546│ 36.6761│ 14.5143│ │
│ C4FH13│ 0.0│ 0.0│ 0.0│ │
│ C3S│ 0.0│ 0.0│ 0.0│ │
│ C2S│ 0.0│ 0.0│ 0.0│ │
│ C3A│ 0.0│ 0.0│ 0.0│ │
│ C4AF│ 0.0│ 0.0│ 0.0│ │
│ ettringite│ 0.0│ 0.0│ 0.0│ │
│ Gp│ 0.0│ 0.0│ 0.0│ │
│ Anh│ 0.0│ 0.0│ 0.0│ │
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│ # TOTAL #│ n [mol]│ m [g]│ V [cm³]│ │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ │ 10.636│ 1000.0│ 469.693│ │
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│ pH : 12.3924 │
│ pOH : 1.6077 │
│ porosity : 0.175514 │
│ saturation : 1.0 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘What went in and what came out, the solids and the free water, in grams per kilogram of paste:
using Plots
grams(st, i) = ustrip(us"g", st.n[i] * st.system.species[i][:M])
function inventory(st)
out = Dict(symbol(st.system.species[i]) => grams(st, i) for i in st.system.idx_crystal if grams(st, i) > 0.1)
out["free water"] = grams(st, findfirst(sp -> symbol(sp) == "H2O@", st.system.species))
return out
end
before, after = inventory(state), inventory(state_eq)
names_ba = sort(collect(union(keys(before), keys(after))); by = k -> -max(get(before, k, 0.0), get(after, k, 0.0)))
xs = collect(eachindex(names_ba))
fig = bar(xs .- 0.2, [get(before, k, 0.0) for k in names_ba]; bar_width = 0.4, label = "before",
color = :grey60, xticks = (xs, names_ba), xrotation = 30, ylabel = "g per kg of paste",
ylims = (0, 1.1 * maximum(max(get(before, k, 0.0), get(after, k, 0.0)) for k in names_ba)),
size = (820, 400), left_margin = 6Plots.mm, bottom_margin = 10Plots.mm, titlefontsize = 11,
title = "Clinker and gypsum at w/c 0.4, then at equilibrium, 25 °C")
bar!(fig, xs .+ 0.2, [get(after, k, 0.0) for k in names_ba]; bar_width = 0.4, label = "at equilibrium",
color = :steelblue)
fig
Quick shortcut
If your ChemicalState was constructed from a ChemicalSystem built with CEMDATA_PRIMARIES, calling equilibrate(state) with no other arguments is usually sufficient for cement-chemistry problems.
Inspecting the equilibrium state
The returned ChemicalState carries all derived thermodynamic quantities:
println("pH = ", pH(state_eq))
println("pOH = ", pOH(state_eq))
println("porosity = ", porosity(state_eq))
println("saturation = ", saturation(state_eq))Phase volumes and mole amounts are accessible via named tuples:
v = volume(state_eq)
println("V liquid = ", v.liquid)
println("V solid = ", v.solid)
println("V total = ", v.total)
m = moles(state_eq)
println("n liquid = ", m.liquid)
println("n solid = ", m.solid)Individual species amounts (in mol):
cs_eq = state_eq.system
for (i, sp) in enumerate(cs_eq.species)
n_i = state_eq.n[i]
println(rpad(symbol(sp), 20), ustrip(n_i), " mol")
endControlling the solver
Variable space: :linear vs :log
The choice between the two variable spaces, and why the logarithmic one is a refinement rather than a start, is explained in Variable space: :linear vs :log.
Convergence
Solving a system of equations in chemistry can be a difficult undertaking. The orders of magnitude can vary greatly, and convergence is not guaranteed.
Tolerances
Solver options are passed as keyword arguments and forwarded to the optimizer; which ones each back end reads is in the manual.
Using EquilibriumSolver explicitly
For batch calculations where many different initial states share the same system and activity model, construct an EquilibriumSolver once and reuse it:
using Optimization, OptimizationIpopt
opt = IpoptOptimizer(
acceptable_tol = 1e-12,
dual_inf_tol = 1e-12,
acceptable_iter = 1000,
constr_viol_tol = 1e-12,
warm_start_init_point = "no",
additional_options = Dict{String, Any}("print_level" => 0, "sb" => "yes"),
)
solver = EquilibriumSolver(
cs,
DiluteSolutionModel(),
opt;
variable_space = Val(:linear),
reltol = 1e-10,
)Ipopt prints an iteration table for every solve unless told otherwise; print_level = 0 and sb = "yes" (no banner) keep it silent, and the default solver of equilibrate is built the same way. reltol becomes Ipopt's convergence tolerance tol.
Once built, solver is called with any compatible ChemicalState:
state_eq2 = solve(solver, state)Performance
The potential function μ(n, p) is compiled once during EquilibriumSolver construction. Repeated calls to solve(solver, ...) with different states reuse it, avoiding redundant compilation overhead.
Activity models
All activity models inherit from AbstractActivityModel. There are three built-in ones — DiluteSolutionModel, HKFActivityModel and DaviesActivityModel; see Activity models for the other two and for Reading the aqueous properties back. This page uses the first, which implements:
| Phase | Law | Expression |
|---|---|---|
| Solvent (H₂O) | Raoult | ln a = ln xₛ |
| Aqueous solutes | Henry | ln a = ln(cᵢ / c°), c° = 1 mol/L |
| Crystals | Pure solid | ln a = 0 |
| Gas | Ideal mixture | ln a = ln xᵢ + ln(P/P°) |
To implement a custom activity model, define a new subtype and extend activity_model:
struct MyModel <: AbstractActivityModel
# model parameters
end
function ChemistryLab.activity_model(::ChemicalSystem, ::MyModel)
# return a closure lna(n, p) -> Vector{Float64}
return (_, _) -> begin
# compute and return log-activities
end
endPass your model to equilibrate or EquilibriumSolver:
state_eq = equilibrate(state; model=MyModel(...))Where to go next
The same system is swept over the water-to-cement ratio in Effect of the water/cement ratio on the hydration of a CEM I. A clinker described by its four main phases and its sulfate carrier, rather than simplified, is A CEM I from its clinker phases.