Chemical Kinetics

Before this page

The tutorial Chemical Equilibrium, since every step of a kinetic run re-equilibrates the solution, and the theory page Rate laws, and every parameter in them for the rate laws themselves.

This tutorial runs a complete kinetic calculation, the hydration of a Portland clinker, in which the dissolution and precipitation of minerals are governed by rate laws while the aqueous speciation, much faster, is kept at equilibrium, following the methodology of Leal et al. (2017).

Background

Each kinetic mineral advances along its reactions at the rates its laws give, written in Rate laws, and every parameter in them; the aqueous speciation, re-equilibrated as the run advances, provides the activities and the saturation ratios those laws may read. The equations of the coupling are in Kinetics under partial equilibrium.

Two routes, and the one followed here

A kinetic calculation advances in time by one of two routes. The implicit step, kinetic_step, performs one Gibbs minimization per time step with the reaction extents among the unknowns, so that the products are decided by thermodynamics; the integration of an ordinary differential equation, integrate, follows the stoichiometry written in the kinetic reactions and leaves the step length to an adaptive solver. The choice between them, their respective limits and the implicit step in full are the subject of Which route: two ways to advance in time. The run below follows the second route on a Portland clinker, and every call it makes is described in Writing a kinetic model.

Full workflow: OPC clinker hydration with KineticsProblem

Complete chain — ChemicalSystem → ChemicalState → KineticReaction → KineticsProblem → integrate:

using ChemistryLab, OrdinaryDiffEq, DynamicQuantities, OrderedCollections, Printf

# ── 1. ChemicalSystem from CEMDATA18 ────────────────────────────────────────
DATA_FILE = datapath("cemdata18-thermofun.json")
substances = build_species(DATA_FILE)

input_species = split(
    "C3S C2S C3A C4AF " *
    "Portlandite Jennite ettringite monosulphate12 C3AH6 C3FH6 H2O@",
)
species = speciation(substances, input_species; aggregate_state = [AS_AQUEOUS])
cs = ChemicalSystem(species, CEMDATA_PRIMARIES)

# ── 2. Initial state: 1 kg OPC (CEM I 52.5 R), w/c = 0.40 ─────────────────
WC          = 0.40
COMPOSITION = (C3S=0.619, C2S=0.165, C3A=0.080, C4AF=0.087)   # assumed, typical of the class

state0 = ChemicalState(cs)
for (name, frac) in pairs(COMPOSITION)
    set_quantity!(state0, string(name), frac * u"kg")
end
set_quantity!(state0, "H2O@", WC * u"kg")

# ── 3. [ParrottKilloh1984](@citet) rate functions with Powers α_max ───────────────────────
α_max   = powers_alpha_max(WC)
BLAINE  = 380.0u"m^2/kg"
pk_C3S  = parrott_killoh_avrami(PK84_PARAMS_C3S,  "C3S";  α_max, blaine = BLAINE)
pk_C2S  = parrott_killoh_avrami(PK84_PARAMS_C2S,  "C2S";  α_max, blaine = BLAINE)
pk_C3A  = parrott_killoh_avrami(PK84_PARAMS_C3A,  "C3A";  α_max, blaine = BLAINE)
pk_C4AF = parrott_killoh_avrami(PK84_PARAMS_C4AF, "C4AF"; α_max, blaine = BLAINE)

# ── 4. Kinetic reactions (reaction-centric) ─────────────────────────────────
# Reactions follow [LothenbachWinnefeld2006](@citet) — Jennite = Ca₉Si₆O₁₈(OH)₆·8H₂O
# Balanced hydration reactions — ΔᵣH⁰ computed from species ΔₐH⁰.
sp(name) = cs[name]

rxn_C3S = Reaction(
    OrderedDict(sp("C3S") => 1.0, sp("H2O@") => 103/30),
    OrderedDict(sp("Jennite") => 1.0, sp("Portlandite") => 4/3);
    symbol = "C₃S hydration",
)
rxn_C3S[:rate] = pk_C3S

rxn_C2S = Reaction(
    OrderedDict(sp("C2S") => 1.0, sp("H2O@") => 73/30),
    OrderedDict(sp("Jennite") => 1.0, sp("Portlandite") => 1/3);
    symbol = "C₂S hydration",
)
rxn_C2S[:rate] = pk_C2S

rxn_C3A = Reaction(
    OrderedDict(sp("C3A") => 1.0, sp("H2O@") => 6.0),
    OrderedDict(sp("C3AH6") => 1.0);
    symbol = "C₃A hydration",
)
rxn_C3A[:rate] = pk_C3A

rxn_C4AF = Reaction(
    OrderedDict(sp("C4AF") => 1.0, sp("Portlandite") => 2.0, sp("H2O@") => 10.0),
    OrderedDict(sp("C3AH6") => 1.0, sp("C3FH6") => 1.0);
    symbol = "C₄AF hydration",
)
rxn_C4AF[:rate] = pk_C4AF

# ── 5. Problem + semi-adiabatic calorimeter ─────────────────────────────────
# `Cp` is the vessel alone: the heat capacity of the paste, Σᵢ nᵢ Cp°ᵢ(T), is
# summed over its composition at every step and must not be counted again.
cal = SemiAdiabaticCalorimeter(;
    Cp        = 900.0u"J/K",
    T_env     = 293.15u"K",
    heat_loss = ΔT -> 0.30 * ΔT + 0.003 * ΔT^2,
    T0        = 293.15u"K",
)

kp = KineticsProblem(
    cs, [rxn_C3S, rxn_C2S, rxn_C3A, rxn_C4AF], state0, (0.0, 7.0 * 86400.0);
    calorimeter = cal,
    equilibrium_solver = nothing,
)

# ── 6. Integrate and post-process ───────────────────────────────────────────
ks  = KineticsSolver(; ode_solver = Rodas5P(), reltol = 1e-6, abstol = 1e-9)
sol = integrate(kp, ks)

_, T_vec = temperature_profile(sol, cal)
_, Q_vec = cumulative_heat(sol, cal)

α_C3S = degrees_of_hydration(sol, kp)["C3S"]

@printf("ΔT_max = %.2f °C   Q_7d = %.1f kJ/kg   α(C3S) = %.4f\n",
        maximum(T_vec .- 273.15) - 20.0, Q_vec[end] / 1000, α_C3S[end])

Choosing equilibrium_solver

Setting equilibrium_solver = nothing skips the Gibbs minimization, which is appropriate for the (Parrott and Killoh, 1984) model: its rate closure ignores the lna argument, so re-speciation cannot change α(t) or the calorimetry. It does change what the products are — with nothing, the hydrate assemblage is whatever the hand-written reaction stoichiometry says, not what thermodynamics gives. For transition_state models, whose rate depends on the saturation ratio Ω = IAP/K, an EquilibriumSolver is required for the rates themselves to mean anything.

The solver may be passed on the KineticsProblem or on the KineticsSolver; both work. If it is set on both, the one on the problem is used and the conflict is reported.

The equilibrium–kinetics coupling

With an equilibrium_solver on the problem, the species are split into a kinetic partition, the minerals carrying a rate law, and an equilibrium partition, everything else, which is re-equilibrated as the run advances. The ODE state is (bₑ, nₖ): the element amounts of the equilibrium partition and the moles of the kinetic minerals. Why it carries element amounts, the equations it advances and when the partition may be frozen within a step are in Kinetics under partial equilibrium; with the rate laws of this page, which read only the kinetic amounts, it is solved once per accepted step. The partition of a problem can be read off it:

[symbol(kp.system.species[i]) for i in kp.idx_kinetic]       # the kinetic minerals
[symbol(kp.system.species[i]) for i in kp.idx_equilibrium]   # everything re-equilibrated

Each kinetic reaction is oriented so that its controlling mineral carries ν = −1: a positive rate is a dissolution.

A failed re-speciation is reported, not hidden

If the equilibrium solve fails, that step keeps its frozen composition, the first failure is reported with its exception, and integrate states how many steps were affected. A run in which re-speciation never succeeded must not look like a healthy one.

Calorimetry and ΔᵣH⁰

Without an equilibrium partition the heat is that of the kinetic reactions, from their reaction enthalpies, which are built from the ΔₐH⁰ of the species — so reactions must be mass-balanced (The calorimeters). For custom species that lack a ΔₐH⁰ entry (GGBS, MK, …), set :ΔᵣH⁰ directly on the reaction: rxn[:ΔᵣH⁰] = NumericFunc((T,) -> -36_100.0, (:T,), u"J/mol").

Post-processing a kinetics run

The ODE state vector carries only the kinetic species. Every other amount is held in a single buffer that the integrator mutates in place, so after a run it reflects the last accepted step and nothing else. Recovering the composition at an arbitrary time therefore means replaying the stoichiometry, which is what reaction_extents and state_at do.

sol = integrate(kp, ks)

α  = degrees_of_hydration(sol, kp)          # Dict(symbol => α(t))
ᾱ  = mean_degree_of_hydration(sol, kp)      # mass-weighted binder average

ξ  = reaction_extents(sol, kp)              # extent of each reaction, mol
st = state_at(sol, kp, 28 * 86400.0)        # full ChemicalState at 28 days

state_at rebuilds every species from n = n₀ + νᵀξ, so the result satisfies A(n - n₀) = (Aνᵀ)ξ to machine precision — whatever the quadrature error. How well the elements themselves balance is a property of the reactions you wrote, not of the reconstruction; check it with maximum(abs, A * ν').

Re-speciation is not replayed

When the run used an equilibrium solver, the aqueous partition was re-speciated at every accepted step, and that redistribution cannot be recovered from the stoichiometry. state_at then returns the purely kinetic reconstruction; call equilibrate on it to speciate.

From moles to volume fractions

volume_fractions turns a state into the input a mean-field homogenization scheme consumes, using the standard molar volumes V⁰ that CEMDATA18 (Lothenbach et al., 2019) supplies for every cement phase.

groups = [
    "anhydrous" => ["C3S", "C2S", "C3A", "C4AF"],
    "C-S-H"     => "Jennite",
    "CH"        => "Portlandite",
    "AFt"       => "ettringite",
    "water"     => "H2O@",
]
f = volume_fractions(st, groups; reference = state0)

Passing reference selects the sealed-curing convention of (Lavergne et al., 2018): fractions are referred to the initial volume, held fixed, and the deficit left by the reactions appears as a "void" phase. That void is the chemical shrinkage — hydration products occupy less space than the reactants they consume — and it is a genuine phase of the microstructure, not a rounding error. Without reference, fractions are referred to the current volume and the void does not exist.

sum(values(f))                              # 1.0, void included
f["void"]                                   # Le Chatelier contraction
chemical_shrinkage(st, state0)              # the same thing, as a volume

Species without a molar volume are invisible

A species carrying no V⁰ contributes nothing to volume, porosity or volume_fractions — a custom species added by hand, for instance. Call missing_molar_volumes before trusting a volume balance:

isempty(missing_molar_volumes(st)) || @warn "incomplete volume balance"

Individual fractions may be slightly negative for aqueous solutes, whose standard partial molar volumes are negative (electrostriction contracts the solvent around an ion). They are kept so that volume_fractions and volume always agree; grouping folds them back into the liquid.

Where to go next

The syntax used above — rate functions, rate constants, kinetic reactions, calorimeters — is described call by call in Writing a kinetic model. The coupling between kinetics and equilibrium is derived in Coupling kinetics and equilibrium, and the application pages carry it further: Hydration kinetics of a CEM I 52.5 R clinker on the same clinker, The silicates of a CEM I clinker, hydrating end to end with a hydrate assemblage computed rather than imposed, and Calibrating the hydration kinetics of a CEM I on measured calorimetry when the parameters themselves are the question.