ChemistryLab.jl

Appearance

The hydrating paste, end to end ​

A worked application of Coupling kinetics and equilibrium: alite and belite dissolve according to (Parrott and Killoh, 1984), and the hydrate assemblage that forms is computed by Gibbs minimization rather than imposed from a stoichiometric recipe.

Read the theory page first if you have not — this one assumes the partition and the map  .

What is different from a stoichiometric model ​

The usual approach writes reactions like

C₃S + 5.3 H → C₁.₇SH₄ + 1.3 CH

and integrates the extent. It is serviceable, but the products and their proportions are decided in advance. Nothing tells you the pore solution pH, and nothing adapts if you change the water/cement ratio, add limestone, or leach the paste.

Here the clinker dissolution rate is prescribed — that is genuinely kinetic — and everything downstream follows from thermodynamics: which hydrates appear, in what amounts, and what the pore solution looks like.

Every block below is executed when this page is built, so the numbers are whatever the code actually produced.

1. Species and system ​

Take the two silicate clinker phases, the hydrates Cemdata18 offers for them, and the aqueous species speciation pulls in from the primaries.

julia
using ChemistryLab, DynamicQuantities, OptimaSolver, OrdinaryDiffEq, Printf

data = datapath("cemdata18-thermofun.json")
substances = build_species(data)

anhydrous = ["C3S", "C2S"]
hydrates = ["Portlandite", "Jennite"]

sp = speciation(substances, vcat(anhydrous, hydrates); aggregate_state = [AS_AQUEOUS])
cs = ChemicalSystem(sp, CEMDATA_PRIMARIES)

@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:  82%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████                             |  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:00
17 species: 13 aqueous, 4 crystalline

2. The initial state ​

A paste at w/c = 0.40: one kilogram of clinker, 65 % alite and 11 % belite by mass, and 400 g of water.

julia
state0 = ChemicalState(cs)
set_quantity!(state0, "C3S", 0.65u"kg")
set_quantity!(state0, "C2S", 0.11u"kg")
set_quantity!(state0, "H2O@", 0.40u"kg")

for s in anhydrous
    @printf "%-5s %8.4f mol\n" s ustrip(us"mol", moles(state0, s))
end
C3S     2.8470 mol
C2S     0.6387 mol

3. Dissolution reactions ​

The clinker does not react into hydrates here — it dissolves into ions, and ChemistryLab balances each reaction from the phase and the primaries. The acid-driven form is what drives the pore solution alkaline.

julia
primaries = ["Ca+2", "SiO2@", "H2O@", "H+"]
reactions = ChemistryLab.Reaction[]
for (phase, pk) in (("C3S", PK84_PARAMS_C3S), ("C2S", PK84_PARAMS_C2S))
    rxn = Reaction([cs[phase]], [cs[p] for p in primaries]; symbol = "$phase dissolution")
    rxn[:rate] = parrot_killoh_avrami(pk, phase; α_max = 1.0, blaine = 380.0u"m^2/kg")
    push!(reactions, rxn)
    println("  ", rxn)
end
  (CaO)₃SiO₂ = 3H₂O@ + 3Ca²⁺ + SiO₂@ + (-6)H⁺
  (CaO)₂SiO₂ = 2H₂O@ + 2Ca²⁺ + SiO₂@ + (-4)H⁺

The negative H⁺ coefficient is the whole point: dissolving one mole of alite consumes six moles of protons, which is what a pore solution at pH 12.5 and above records.

4. The coupled problem ​

The reaction list is a positional argument, and equilibrium_solver belongs on the problem. Without it the run is a pure kinetics integration and the aqueous phase never re-speciates; with it, respeciate! solves once per accepted step.

The activity model matters. A cement pore solution sits at an ionic strength of 0.1–0.7 mol/kg, where the dilute model is not defensible, so pass HKFActivityModel — to the problem and to the solver, which must agree.

julia
model = HKFActivityModel()
kp = KineticsProblem(
    cs, reactions, state0, (0.0, 7 * 86400.0);
    activity_model = model,
    equilibrium_solver = EquilibriumSolver(cs, model, OptimaOptimizer()),
)
sol = integrate(kp, KineticsSolver(; ode_solver = Rodas5P(), reltol = 1.0e-7, abstol = 1.0e-10))
@printf "%d accepted steps, retcode = %s\n" length(sol.t) sol.retcode
┌ 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: 220 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.87e-7 mol ON THE ACCEPTED STEPS, i.e. on the trajectory itself; 9.87e-7 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
149 accepted steps, retcode = Success

5. Reading the result ​

The composition at a given time is not state_at: that returns the purely kinetic reconstruction  , and the redistribution performed by the equilibrium solve cannot be recovered from the stoichiometry. Use speciated_states, which re-solves from the element totals the run carried, walking the instants in order.

julia
times = [1.0, 3.0, 7.0] .* 86400
states = speciated_states(sol, kp; times = times)

@printf "%6s %8s %10s %12s %10s\n" "t [d]" "pH" "Jennite" "Portlandite" "H2O"
for (t, st) in zip(times, states)
    @printf "%6.1f %8.2f %10.4f %12.4f %10.3f\n" t / 86400 something(pH(st), NaN) (
        ustrip(us"mol", moles(st, "Jennite"))
    ) ustrip(us"mol", moles(st, "Portlandite")) ustrip(us"mol", moles(st, "H2O@"))
end
┌ 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
 t [d]       pH    Jennite  Portlandite        H2O
   1.0    12.59     1.4243       1.7968     17.409
   3.0    12.59     2.0191       2.5069     15.451
   7.0    12.59     2.3864       2.9252     14.262

Check the partition constraint — conservation of matter is one of the three conditions a certificate rests on, and the one with an immediate meaning:

julia
p = sol.prob.p
be = collect(sol(times[end])[1:(p.n_be)])
ne = Float64[ustrip(us"mol", moles(states[end], symbol(cs.species[j]))) for j in kp.idx_equilibrium]
@printf "‖Aₑnₑ − bₑ‖∞ = %.2e mol\n" maximum(abs, p.Ae * ne .- be)
‖Aₑnₑ − bₑ‖∞ = 1.03e-13 mol

Degrees of hydration come from the kinetic amounts against their initial values:

julia
α = degrees_of_hydration(sol, kp; times = times)
for s in anhydrous
    @printf "α(%s) = %.3f\n" s α[s][end]
end
α(C3S) = 0.750
α(C2S) = 0.393

6. What the thermodynamics decided ​

Nothing above states that portlandite forms, or in what proportion to the C-S-H. The dissolution reactions produce only Ca²⁺, SiO₂ and H⁺; the assemblage is the minimizer of the Gibbs energy under the element totals the kinetics delivered.

julia
st = states[end]
solids = [(symbol(s), ustrip(us"mol", moles(st, symbol(s)))) for s in cs.species
    if aggregate_state(s) == AS_CRYSTAL]
filter!(x -> last(x) > 1.0e-6, solids)
sort!(solids; by = last, rev = true)
for (name, amount) in solids
    @printf "  %-14s %8.4f mol\n" name amount
end
  Portlandite      2.9252 mol
  Jennite          2.3864 mol
  C3S              0.7118 mol
  C2S              0.3874 mol

The pore solution comes with it:

julia
aq = [(symbol(s), ustrip(us"mol", moles(st, symbol(s)))) for s in cs.species
    if aggregate_state(s) == AS_AQUEOUS && symbol(s) != "H2O@"]
sort!(aq; by = last, rev = true)
for (name, amount) in aq[1:min(5, end)]
    @printf "  %-14s %.3e mol\n" name amount
end
  OH-            9.909e-03 mol
  Ca+2           4.393e-03 mol
  CaOH+          1.122e-03 mol
  CaSiO3@        9.463e-06 mol
  HSiO3-         3.122e-07 mol

Caveats worth carrying ​

  • The C-S-H is Jennite, the Cemdata18 end-member normalized per silicon. A real C-S-H is a solid solution of varying Ca/Si; using the CSHQ solid solution inside a coupled run is not exercised by this package's tests.

  • The interior-point optimizer rarely reports convergence on a cement equilibrium, and its return code is not the thing to read. integrate reports the worst element balance instead, separating the accepted steps — the trajectory — from the Jacobian probes that never enter it. The compositions above do not rest on that solver: speciated_states passes each instant to DualEquilibriumSolver, which solves the KKT system and certifies the result. The problem is convex, so stationarity of the interior species, the component balance, and undersaturation of every absent phase together prove global optimality.

  • A Parrot–Killoh rate reads only its own degree of reaction, so the trajectory here does not depend on the speciation at all; the speciation is what you read out of it. A rate law reading log-activities would feed the speciation back into the trajectory, and would need the balance above to be tight at every step, not just at the instants you ask for.

7. The certificate ​

Nothing above asks you to take the composition on trust. G is convex — an ideal mixing entropy plus terms linear in the amounts of the pure phases, over a polyhedron — so its minimizer is unique and the KKT conditions are sufficient. optimality_certificate checks them.

julia
sub = ChemistryLab._equilibrium_subsystem(kp.system, kp.idx_equilibrium)
des = DualEquilibriumSolver(sub, model)
snames = symbol.(sub.species)
be = collect(sol(times[end])[1:(sol.prob.p.n_be)])
sub_state = ChemicalState(
    sub, [ustrip(us"mol", moles(states[end], s)) for s in snames] .* u"mol";
    T = sol.prob.p.T_q[], P = sol.prob.p.P_q[],
)
cert = optimality_certificate(des, sub_state; b = be)
@printf "optimal          %s\n" cert.optimal
@printf "stationarity     %.2e  (RT, over %d interior species)\n" cert.stationarity cert.n_interior
@printf "component balance %.2e mol\n" cert.balance
@printf "worst supersaturation %+.2e  (negative: every absent phase undersaturated)\n" cert.worst_supersaturation
optimal          true
stationarity     2.16e-15  (RT, over 12 interior species)
component balance 1.03e-13 mol
worst supersaturation -Inf  (negative: every absent phase undersaturated)