Rate laws, and every parameter in them
A Gibbs minimization answers what is stable. It says nothing about when, and a cement paste is a material whose whole engineering behavior lives in the when. This page is the catalog of what supplies that: every rate law the package ships, the physics each one encodes, the parameters each one carries, their units, and where the numbers come from.
It is written to be read beside the kinetics tutorial rather than instead of it: the tutorial says how to call these, this page says what they mean and what they are worth.
using ChemistryLab
using DynamicQuantities
using Printf
using Plots
default(framestyle = :box, grid = false)0. What a rate law is, here
Every rate in this package is a KineticFunc with one signature:
r = f(T, P, t, n::StateView, lna::StateView, n_initial::StateView) # mol/sSix arguments, and the choice of six is the point. A rate may depend on the temperature and pressure, on the elapsed time, on the current composition, on the activities of any species in the system, and on the initial composition — which is what lets a law be written in terms of a degree of reaction, α = 1 − n/n₀, without the caller having to track it.
Two consequences worth stating plainly:
A rate law may read any species, not only the ones in its own reaction. A dissolution rate catalyzed by
H+, an SCM whose reaction is gated by the portlandite left, a mechanism inhibited by the sulfate in solution — all of these are written by readingnorlnaat the species concerned. The reaction a rate is attached to fixes what the rate consumes; it does not restrict what the rate may depend on. One consequence has its own section in the tutorial, because it bites: a rate that reads a species the reaction itself consumes to exhaustion — see rate laws that depend on a consumed reactant.The evaluation path is differentiable. No
Float64casts, so aForwardDiff.Dualpropagates throughT, throught, through the composition and through the parameters. That is what makes a rate constant fittable (the calibration page) rather than only adjustable by hand.
There are two families below, and they answer different questions. The first is mechanistic: a rate proportional to how far the solution is from equilibrium with the mineral. The second is empirical: a degree of reaction as a function of time, fitted to calorimetry. Clinker and supplementary materials are described by the second, because nobody can write the first for a multiphase glass.
1. Distance from equilibrium: the Palandri–Kharaka form
For a mineral dissolving into or precipitating from a solution, the rate is driven by the saturation ratio
computed by saturation_ratio from the stoichiometry, the log-activities and the standard Gibbs energies the database supplies. Ω < 1 is undersaturation and the mineral dissolves; Ω > 1 and it precipitates; Ω = 1 is equilibrium and the rate is exactly zero — which is the property that makes this form composable with an equilibrium solver rather than in competition with it.
One mechanism (RateMechanism) contributes
and a mineral's total rate is the sum over its mechanisms — classically an acid one, a neutral one and a base one, which is why the products of activities are called catalysts (RateModelCatalyst): each contributes a_j^{n_j}, with n_j = 0.5 on H+ a typical acid mechanism.
| symbol | meaning | unit |
|---|---|---|
rate constant, usually arrhenius_rate_constant | mol m⁻² s⁻¹ | |
| catalyst activity and its exponent | — | |
| saturation exponent inside the bracket | — | |
| outer exponent | — | |
| saturation ratio | — |
The rate constant itself is the Arrhenius law referred to a reference temperature rather than to an absolute prefactor,
so that k₀ is the measured constant at T_ref and carries its own unit, instead of being an extrapolation to infinite temperature. The shipped parameter sets for this family come from (Palandri and Kharaka, 2004).
2. Clinker hydration: Parrot & Killoh
A clinker phase does not dissolve into a solution at a rate set by its own undersaturation — it is consumed behind a growing layer of hydrate, and the controlling step changes as the layer thickens. The 1984 description, as reported by (Lothenbach et al., 2008) and used by (Lavergne et al., 2018), writes three competing mechanisms in the degree of hydration and lets the slowest one limit:
Each has a physical reading.
the three multipliers being the corrections of section 4.
sets = ["C₃S" => PK84_PARAMS_C3S, "C₂S" => PK84_PARAMS_C2S,
"C₃A" => PK84_PARAMS_C3A, "C₄AF" => PK84_PARAMS_C4AF]
@printf("%-6s %8s %6s %10s %8s %6s %10s\n",
"phase", "k1 [1/d]", "n1", "k2 [1/d]", "k3 [1/d]", "n3", "Ea [kJ/mol]")
for (nm, p) in sets
@printf("%-6s %8.2f %6.2f %10.3f %8.2f %6.1f %10.0f\n",
nm, ustrip(us"1/d", p.k₁), p.n₁, ustrip(us"1/d", p.k₂),
ustrip(us"1/d", p.k₃), p.n₃, ustrip(us"J/mol", p.Ea) / 1000)
endphase k1 [1/d] n1 k2 [1/d] k3 [1/d] n3 Ea [kJ/mol]
C₃S 1.50 0.70 0.050 1.10 3.3 42
C₂S 0.50 1.00 0.006 0.20 5.0 21
C₃A 1.00 0.85 0.040 1.00 3.2 54
C₄AF 0.37 0.70 0.015 0.40 3.7 32Two features of that table are artifacts of the 1984 fit rather than chemistry, and both are useful as checks on an implementation. With
Two Parrot–Killoh variants ship, and only one is attributed
parrot_killoh is a different, smoothed variant that predates this one in the package (the two are set side by side in the tutorial) — min(max(r_{NG}, r_I), r_D) with a damped nucleation term — and its PK_PARAMS_* are not transferable to the canonical law. It is deprecated and the attribution to Parrott & Killoh was withdrawn rather than repaired: its nucleation term carries no Avrami logarithm, its shell coefficient sits in the diffusion expression, and no published set matches its parameters. The primary source is a conference proceedings without a DOI that could not be consulted, so the honest action was to stop claiming it. Measured consequence: with PK_PARAMS_* all four phases land on the diffusion branch almost immediately and a CEM I at w/c 0.40 reaches α ≈ 0.234 at seven days against the ≈ 0.61 the literature reports. Use parrot_killoh_avrami.
3. Supplementary materials: the Waller sigmoid
A slag or a fly ash has no phases, no formula and no single dissolution mechanism, so its reaction is described by the shape it is observed to have — a sigmoid in log time:
The second form is what the package evaluates, because writing the rate as a function of the current degree rather than of the clock is what lets the same temperature, fineness and humidity corrections multiply it as they multiply the clinker rate.
waller_sets = ["fly ash" => WALLER_PARAMS_FLY_ASH, "silica fume" => WALLER_PARAMS_SILICA_FUME,
"slag" => WALLER_PARAMS_SLAG]
α_waller(p, t_days) = 1 / (1 + (ustrip(us"d", p.τ) / t_days)^p.n)
@printf("%-12s %8s %6s %12s %10s %10s\n",
"material", "τ [d]", "n", "Ea [kJ/mol]", "α(28 d)", "α(90 d)")
for (nm, p) in waller_sets
@printf("%-12s %8.0f %6.2f %12.2f %10.3f %10.3f\n",
nm, ustrip(us"d", p.τ), p.n, ustrip(us"J/mol", p.Ea) / 1000,
α_waller(p, 28.0), α_waller(p, 90.0))
endmaterial τ [d] n Ea [kJ/mol] α(28 d) α(90 d)
fly ash 80 0.70 83.14 0.324 0.521
silica fume 80 0.70 83.14 0.324 0.521
slag 100 0.70 83.14 0.291 0.482Pozzolanic and latent-hydraulic reactions are markedly more temperature sensitive than the clinker's: 83.14 kJ/mol against 21–54 kJ/mol above. That is why a blended cement gains so much from a warm cure and loses so much in a cold one, and it is a prediction of the parameter and not an extra rule.
These are not the same numbers as a direct measurement of reacted glass
The α(28 d) column above is ≈ 0.29 for slag and ≈ 0.32 for fly ash. The RILEM TC 238-SCM round robin (Durdziński et al., 2017), which measured the reacted glass directly on its own materials rather than inferring it from heat, reports 38–49 % for two slags and about 20 % for a siliceous fly ash at the same age. The two disagree, and in opposite directions.
Neither is wrong. Waller's parameters are a fit to particular materials, and "a slag" is not a substance — its reactivity depends on its glass content, its basicity and its fineness. The round robin's own conclusion is the one to keep: the precision of any determination of an SCM's degree of reaction is "rather low, at best ± 4-5 %". So a reacted fraction is an input to be stated and swept, never a constant to be trusted to two digits — which is exactly how the blended binder pages treat it.
Silica fume carries the same τ and n as fly ash: its much higher reactivity is represented through the fineness, at an effective Blaine of 2000 m²/kg recommended by (Lavergne et al., 2018). Its BET surface, about 20 000 m²/kg, is a different measurement of a different thing and must not be substituted.
4. The three multiplicative corrections
Each of the empirical laws above is multiplied by three dimensionless factors. They are separated because they answer separate questions, and because each can be switched off by leaving its keyword at nothing.
Temperature, the same Arrhenius factor as section 1:
Fineness (blaine_factor): the rate scales linearly with the specific surface, relative to the fineness the parameters were fitted at — 385 m²/kg for the clinker phases, 400 m²/kg for the Waller sets.
Internal humidity (humidity_factor): hydration stops when the pore water is no longer available, which is taken to happen below 80 % RH.
@printf("β_B at 462 m²/kg (clinker ref) : %.3f\n", blaine_factor(462u"m^2/kg"))
@printf("β_B at 2000 m²/kg (Waller ref) : %.3f\n",
blaine_factor(2000u"m^2/kg"; blaine_ref = 400u"m^2/kg"))
for h in (0.99, 0.90, 0.81, 0.801, 0.80, 0.75)
@printf("β_h(%.3f) = %.4f\n", h, humidity_factor(h))
endβ_B at 462 m²/kg (clinker ref) : 1.200
β_B at 2000 m²/kg (Waller ref) : 5.000
β_h(0.990) = 0.9140
β_h(0.900) = 0.3660
β_h(0.810) = 0.1114
β_h(0.801) = 0.0968
β_h(0.800) = 0.0000
β_h(0.750) = 0.0000That cut is a genuine discontinuity, and the page says so rather than smoothing it: the one-sided limit from above is 0.0953 and the value at 0.80 is exactly 0. It is mild in practice — the rate has already fallen by an order of magnitude from β_h(0.99) = 0.914 — but a solver stepping across it will feel it, which is why PoreHumidity exists to supply h from the current saturation rather than from a schedule.
5. How far the reaction can go: α_max
Every law above is written in powers_alpha_max:
@printf("%6s %12s %12s\n", "w/c", "sealed", "under water")
for wc in (0.25, 0.30, 0.36, 0.40, 0.42, 0.50)
@printf("%6.2f %12.3f %12.3f\n", wc,
powers_alpha_max(wc), powers_alpha_max(wc; curing = :saturated))
end w/c sealed under water
0.25 0.595 0.694
0.30 0.714 0.833
0.36 0.857 1.000
0.40 0.952 1.000
0.42 1.000 1.000
0.50 1.000 1.000Read the last two rows: above 0.42 neither convention does anything, because water has stopped being what limits the reaction. Read the first: at w/c = 0.25 a sealed paste cannot pass 60 % hydration however long it is left, and a cured one cannot pass 69 %. The difference between the columns is the chemical shrinkage, refilled from the bath in one case and not in the other.
6. Which numbers are published and which are fitted
| quantity | where it comes from |
|---|---|
PK84_PARAMS_* | Parrott & Killoh (1984) as reported by (Lothenbach et al., 2008), (Lavergne et al., 2018) |
WALLER_PARAMS_* | Waller (1999), as used by (Lavergne et al., 2018) |
| the dissolution rate constants and their exponents | (Palandri and Kharaka, 2004) |
blaine_ref 385 / 400 m²/kg | the finenesses those fits were made at |
| the 0.80 humidity cut and its exponent | Parrot et al., as used by van Breugel |
powers_alpha_max 0.42 / 0.36 | (Powers, 1948) |
CALIBRATED_THETA | fitted here, on one record — the calibration page |
| an SCM's reacted fraction at a given age | an input, measured or assumed — never a constant of the code |
The last two rows are the ones to keep in view. Everything above them is somebody's published fit to somebody's materials; the calibration page shows what happens when those published parameters meet a calorimetry record they were not fitted to, and how much of the gap a five-parameter fit closes.
See also
The water budget of a hydrating paste — where
comes from, and why it is not thermodynamics Calibrating a hydration model — the fit, the holdout, and the sensitivity
A CEM I from its clinker phases — all of this running, over ninety days