Simplified Clinker Dissolution
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, total = 1)
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)┌ Warning: Verbosity toggle: missing_second_order_ad
│ The selected optimization algorithm requires second order derivatives, but `SecondOrder` ADtype was not provided. So a `SecondOrder` with AutoForwardDiff() for both inner and outer will be created, this can be suboptimal and not work in some cases so an explicit `SecondOrder` ADtype is recommended.
└ @ OptimizationBase ~/.julia/packages/OptimizationBase/BEgWO/src/cache.jl:116The solved state in full — every species, with its amount
state_eqChemicalState{Species, AbstractReaction, DynamicQuantities.Quantity{Float64, DynamicQuantities.SymbolicDimensions{DynamicQuantities.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.4392│ │
├┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┤
│ H2O@│ 4.5639│ 82.2186│ 82.4623│ 55.3608│
│ OH-│ 0.00203441│ 0.0345992│ -0.00957768│ 0.0246777│
│ Ca+2│ 0.000846972│ 0.033945│ -0.0156171│ 0.0102739│
│ CaOH+│ 0.000347807│ 0.0198546│ 0.00200423│ 0.00421896│
│ AlO2-│ 7.2335e-6│ 0.000426628│ 6.8482e-5│ 8.77434e-5│
│ CaSiO3@│ 3.0664e-6│ 0.000356192│ 0.0│ 3.71958e-5│
│ HSiO3-│ 6.4551e-8│ 4.97624e-6│ 2.92149e-7│ 7.83014e-7│
│ Ca(SO4)@│ 3.03929e-8│ 4.13751e-6│ 1.42866e-7│ 3.68671e-7│
│ SO4-2│ 1.47856e-8│ 1.42025e-6│ 1.90996e-7│ 1.79352e-7│
│ Ca(HSiO3)+│ 1.05399e-8│ 1.23493e-6│ -7.10022e-8│ 1.2785e-7│
│ FeO2-│ 9.10608e-9│ 7.99905e-7│ 4.11187e-9│ 1.10458e-7│
│ SiO3-2│ 7.47692e-9│ 5.68859e-7│ 0.0│ 9.06962e-8│
│ SiO2@│ 1.68406e-10│ 1.01184e-8│ 2.70516e-9│ 2.0428e-9│
│ AlO2H@│ 8.16864e-12│ 4.90017e-10│ 1.06269e-10│ 9.90869e-11│
│ FeO2H@│ 4.03028e-12│ 3.58095e-10│ 2.90556e-11│ 4.88879e-11│
│ AlSiO5-3│ 2.22885e-12│ 3.01033e-10│ -7.60746e-11│ 2.70363e-11│
│ H+│ 3.31947e-14│ 3.34602e-14│ 0.0│ 4.02656e-13│
│ FeO+│ 2.16722e-292│ 1.55702e-290│ -9.10683e-291│2.6288800000000002e-291│
│ Si4O10-4│9.859679999999999e-305│ 2.68509e-302│ 0.0│ 1.19599e-303│
│ Al(SO4)2-│9.859679999999999e-305│ 2.16019e-302│ 3.06753e-303│ 1.19599e-303│
│ Fe(SO4)@│9.859679999999999e-305│ 1.49769e-302│ 1.64841e-304│ 1.19599e-303│
│ Fe(SO4)+│9.859679999999999e-305│ 1.49769e-302│-2.5984500000000003e-304│ 1.19599e-303│
│ Fe(SO4)2-│9.859679999999999e-305│ 2.44478e-302│3.0057499999999996e-303│ 1.19599e-303│
│ Al(SO4)+│9.859679999999999e-305│ 1.21311e-302│ -5.93457e-304│ 1.19599e-303│
│ Fe(HSO4)+2│9.859679999999999e-305│ 1.50763e-302│2.2873100000000004e-304│ 1.19599e-303│
│ Fe(HSO4)+│9.859679999999999e-305│ 1.50763e-302│ 1.85425e-303│ 1.19599e-303│
│ Al+3│9.859679999999999e-305│ 2.66029e-303│ -4.4608e-303│ 1.19599e-303│
│ Fe+2│9.859679999999999e-305│ 5.50614e-303│ -2.23227e-303│ 1.19599e-303│
│ HS-│9.859679999999999e-305│ 3.2604e-303│1.9925899999999998e-303│ 1.19599e-303│
│ AlOH+2│9.859679999999999e-305│ 4.33713e-303│ -2.68924e-304│ 1.19599e-303│
│ AlO+│9.859679999999999e-305│ 4.23774e-303│3.0302300000000003e-305│ 1.19599e-303│
│ FeOH+2│9.859679999999999e-305│ 7.18297e-303│ -2.49849e-303│ 1.19599e-303│
│ FeOH+│9.859679999999999e-305│ 7.18297e-303│ -1.64784e-303│ 1.19599e-303│
│ HSO3-│9.859679999999999e-305│ 7.99275e-303│ 3.24943e-303│ 1.19599e-303│
│ HSO4-│9.859679999999999e-305│ 9.5702e-303│ 3.43522e-303│ 1.19599e-303│
│ S2O3-2│9.859679999999999e-305│ 1.10544e-302│ 2.72051e-303│ 1.19599e-303│
│ SO3-2│9.859679999999999e-305│ 7.89336e-303│-4.058230000000001e-304│ 1.19599e-303│
│ Fe+3│9.859679999999999e-305│ 5.50614e-303│ -3.72599e-303│ 1.19599e-303│
│ H2S@│9.859679999999999e-305│ 3.35978e-303│ 3.44608e-303│ 1.19599e-303│
│ H2@│9.859679999999999e-305│1.9877100000000003e-304│ 2.49098e-303│ 1.19599e-303│
│ O2@│9.859679999999999e-305│ 3.1549e-303│ 3.00729e-303│ 1.19599e-303│
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│ # 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.695│ │
╞══════════════════════════════════════════════════════════════════════════════════════════════════╡
│ pH : 12.3924 │
│ pOH : 1.6077 │
│ porosity : 0.175517 │
│ saturation : 1.0 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘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
equilibrate accepts a variable_space keyword that selects the optimization variable space:
variable_space | Variables | Recommended when |
|---|---|---|
Val(:linear) | mole amounts nᵢ ≥ 0 | most systems, default |
Val(:log) | log nᵢ | systems spanning many orders of magnitude |
state_eq_log = equilibrate(state; variable_space=Val(: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
Tighter tolerances are passed directly as keyword arguments and forwarded to the underlying Ipopt solver:
state_eq_tight = equilibrate(state; abstol=1e-12, reltol=1e-12)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",
)
solver = EquilibriumSolver(
cs,
DiluteSolutionModel(),
opt;
variable_space = Val(:linear),
abstol = 1e-10,
reltol = 1e-10,
)Once built, solver is called with any compatible ChemicalState:
using Optimization, OptimizationIpopt #hide
opt = IpoptOptimizer( #hide
acceptable_tol = 1e-12, #hide
dual_inf_tol = 1e-12, #hide
acceptable_iter = 1000, #hide
constr_viol_tol = 1e-12, #hide
warm_start_init_point = "no", #hide
) #hide
solver = EquilibriumSolver( #hide
cs, #hide
DiluteSolutionModel(), #hide
opt; #hide
variable_space = Val(:linear), #hide
abstol = 1e-10, #hide
reltol = 1e-10, #hide
) #hide
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ᵢ |
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 Water/Cement Ratio on Cement Hydration. A clinker described by its four main phases and its sulfate carrier, rather than simplified, is A CEM I from its clinker phases.