Stoichiometric matrices
Every calculation in this package rests on one object: a matrix that says how each species is built out of a chosen set of components. Equilibrium is the minimization of Gibbs energy subject to A n = b, and A is that matrix. This page shows it — over atoms, over primary species, and in cement notation — and ends with a matrix whose entries are symbols rather than numbers, which is what lets a phase of variable composition be carried through the algebra.
What the matrix is
Take four clinker phases and ask what each is made of. Answer it in atoms and you get the canonical stoichiometric matrix: one row per element, one column per species, the entry being how many of that element the species contains.
using ChemistryLab
using DynamicQuantities
using Symbolics
using LinearAlgebra
clinker = CemSpecies.(split("C3S C2S C3A C4AF"))
csm_atoms = CanonicalStoichMatrix(clinker)
pprint(csm_atoms.A, csm_atoms.primaries, clinker; label = :name)┌───┬─────┬─────┬─────┬──────┐
│ │ C3S │ C2S │ C3A │ C4AF │
├───┼─────┼─────┼─────┼──────┤
│ C │ 3 │ 2 │ 3 │ 4 │
│ S │ 1 │ 1 │ │ │
│ A │ │ │ 1 │ 1 │
│ F │ │ │ │ 1 │
└───┴─────┴─────┴─────┴──────┘Read a column: C3S is three calcium, one silicon, five oxygen — which is
The same phases over oxide components
Atoms are not the only possible components. A cement chemist decomposes the same four phases over oxides, and the matrix changes shape accordingly — four rows instead of five, because there are four oxides and five elements:
oxides = CemSpecies.(split("C S A F"))
cs_ox = ChemicalSystem(clinker, oxides)
pprint(cs_ox.SM.A, oxides, clinker; label = :name)┌───┬─────┬─────┬─────┬──────┐
│ │ C3S │ C2S │ C3A │ C4AF │
├───┼─────┼─────┼─────┼──────┤
│ C │ 3 │ 2 │ 3 │ 4 │
│ S │ 1 │ 1 │ │ │
│ A │ │ │ 1 │ 1 │
│ F │ │ │ │ 1 │
└───┴─────┴─────┴─────┴──────┘That is the same chemistry in a different basis, and the column for C3S now reads what its name says: three lime, one silica. Neither basis is more correct; the choice is made by what the problem is posed in. A datasheet reports oxides, so the oxide basis is the one the Bogue Calculation inverts.
C is lime here, not carbon
The components above are CemSpecies, so C is S is Species the same letters mean carbon and sulfur, and the matrix would be a different one — see the trap.
Over primary species, from a database
For an aqueous system the natural components are not atoms but a set of primary species — the ions a reaction is conventionally written over. A ChemicalSystem computes that matrix at construction:
substances = build_species(datapath("cemdata18-thermofun.json"); verbose = false)
species = speciation(substances, split("C3S Portlandite H2O@");
aggregate_state = [AS_AQUEOUS])
byname = Dict(symbol(s) => s for s in species)
primaries = [byname[s] for s in CEMDATA_PRIMARIES if haskey(byname, s)]
cs = ChemicalSystem(species, primaries)
length(cs.species), length(cs.SM.primaries)(15, 5)The matrix itself, over the first few species so the table stays readable:
pprint(cs.SM.A[:, 1:8], cs.SM.primaries, cs.species[1:8]; label = :symbol)┌───────┬──────┬─────────┬────────┬──────────┬───────┬──────┬────┬───────┐
│ │ H2O@ │ CaSiO3@ │ SiO3-2 │ Si4O10-4 │ SiO2@ │ Ca+2 │ H+ │ CaOH+ │
├───────┼──────┼─────────┼────────┼──────────┼───────┼──────┼────┼───────┤
│ H2O@ │ 1 │ 1 │ 1 │ 2 │ │ │ │ 1 │
│ Ca+2 │ │ 1 │ │ │ │ 1 │ │ 1 │
│ H+ │ │ -2 │ -2 │ -4 │ │ │ 1 │ -1 │
│ SiO2@ │ │ 1 │ 1 │ 4 │ 1 │ │ │ │
│ Zz │ │ │ │ │ │ │ │ │
└───────┴──────┴─────────┴────────┴──────────┴───────┴──────┴────┴───────┘Each column is a formation reaction: the species in the header, written over the primaries. reactions turns those columns back into reactions:
Every independent reaction of this system, reconstructed from the matrix
reactions(cs.SM)10-element Vector{Reaction}:
Si₄O₁₀⁴⁻ = 2SiO₃²⁻ + 2SiO₂@
SiO₃²⁻ + Ca²⁺ = CaSiO₃@
4SiO₃²⁻ + 4H⁺ = 2H₂O@ + Si₄O₁₀⁴⁻
Si₄O₁₀⁴⁻ + 4Ca(OH)⁺ = 2H₂O@ + 4CaSiO₃@
Si₄O₁₀⁴⁻ + 4OH⁻ = 2H₂O@ + 4SiO₃²⁻
4SiO₃²⁻ + 4Ca(HSiO₃)⁺ = 2H₂O@ + 4CaSiO₃@ + Si₄O₁₀⁴⁻
2H₂@ + O₂@ = 2H₂O@
4HSiO₃⁻ = 2H₂O@ + Si₄O₁₀⁴⁻
Si₄O₁₀⁴⁻ + (CaO)₃SiO₂ = 3CaSiO₃@ + 2SiO₃²⁻
Si₄O₁₀⁴⁻ + 2Ca(OH)₂ = 2H₂O@ + 2CaSiO₃@ + 2SiO₃²⁻A matrix with symbols in it
C-S-H has no fixed formula. Its lime-to-silica ratio, its aluminum uptake and its water content all vary with the cement and with age, so a matrix that describes a hydrated paste has to carry those as parameters rather than numbers.
This is the decomposition of (Chen and Brouwers, 2007) and (Chen and Brouwers, 2007), who write a blended-cement paste over its oxides with a parameterized C-S-H. Take
@variables â b̂ ĝ
CSH = CemSpecies(Dict(:C => â, :S => 1, :A => b̂, :H => ĝ))
hydrates = [
CSH,
CemSpecies("M5AH13"), # hydrotalcite
CemSpecies("C6AFS2H8"), # siliceous hydrogarnet
CemSpecies("C6AS̄3H32"), # ettringite
CemSpecies("C2ASH8"), # strätlingite
CemSpecies("C4AH13"), # hydroxy-AFm
]
csm_h = CanonicalStoichMatrix(hydrates)
pprint(csm_h.A, csm_h.primaries, hydrates; label = :name)┌───┬─────────┬────────┬──────────┬──────────┬────────┬────────┐
│ │ CâSAb̂Hĝ │ M5AH13 │ C6AFS2H8 │ C6AS̄3H32 │ C2ASH8 │ C4AH13 │
├───┼─────────┼────────┼──────────┼──────────┼────────┼────────┤
│ C │ â │ │ 6 │ 6 │ 2 │ 4 │
│ M │ │ 5 │ │ │ │ │
│ S │ 1 │ │ 2 │ │ 1 │ │
│ A │ b │ 1 │ 1 │ 1 │ 1 │ 1 │
│ F │ │ │ 1 │ │ │ │
│ S̄ │ │ │ │ 3 │ │ │
│ H │ ĝ │ 13 │ 8 │ 32 │ 8 │ 13 │
└───┴─────────┴────────┴──────────┴──────────┴────────┴────────┘The â, b̂ and ĝ sit in the matrix exactly where numbers sit for the other phases, and every operation below carries them through.
Inverting it
Seven oxides and six hydrates is not square, so one row must go. Dropping H expresses the hydrates over the anhydrous oxides, water following as a consequence rather than as an input — which is how a hydration calculation is posed:
keep = findall(!=(:H), csm_h.primaries)
A_sq = Symbolics.Num.(csm_h.A[keep, :])
anhydrous = CemSpecies.(string.(csm_h.primaries[keep]))
A_inv = simplify.(inv(A_sq))
pprint(A_inv, hydrates, anhydrous; label = :name)┌──────────┬────────────────────────────────┬────────────────────────────────┬───────────────────────────────┬────────────────────────────────┬────────────────────────────────┬────────────────────────
│ │ C │ M │ S │ A │ F │ ⋯
├──────────┼────────────────────────────────┼────────────────────────────────┼───────────────────────────────┼────────────────────────────────┼────────────────────────────────┼────────────────────────
│ CâSAb̂Hĝ │ 15 / (-30 + 20(3 - 3b̂) + 15â) │ 12 / (-30 + 20(3 - 3b̂) + 15â) │ 30 / (-30 + 20(3 - 3b̂) + 15â) │ -60 / (-30 + 20(3 - 3b̂) + 15â) │ -90 / (-30 + 20(3 - 3b̂) + 15â) │ -10 / (-30 + 20(3 - 3 ⋯
│ M5AH13 │ │ 1//5 │ │ │ │ ⋯
│ C6AFS2H8 │ │ │ │ │ 1 │ ⋯
│ C6AS̄3H32 │ │ │ │ │ │ ⋯
│ C2ASH8 │ -15 / (-30 + 20(3 - 3b̂) + 15â) │ -12 / (-30 + 20(3 - 3b̂) + 15â) │ (-4b̂ + â) / (2 - 4b̂ + â) │ 60 / (-30 + 20(3 - 3b̂) + 15â) │ (2 + 8b̂ - 2â) / (2 - 4b̂ + â) │ 10 / (-30 + 20(3 - 3 ⋯
│ C4AH13 │ (1 - b̂) / (2 - 4b̂ + â) │ (2 - â) / (10 - 20b̂ + 5â) │ (2b̂ - â) / (2 - 4b̂ + â) │ (-2 + â) / (2 - 4b̂ + â) │ (-4 + 2b̂ + â) / (2 - 4b̂ + â) │ (-4 + 6b̂ - â) / (6 - ⋯
└──────────┴────────────────────────────────┴────────────────────────────────┴───────────────────────────────┴────────────────────────────────┴────────────────────────────────┴────────────────────────
1 column omittedEach row is one hydrate written over the anhydrous oxides, as a function of the C-S-H composition. Substituting a composition collapses it to numbers — here a C-S-H at C/S = 1.8 taking a little aluminum:
at = Dict(â => 1.8, b̂ => 0.1, ĝ => 4.0)
A_num = Float64.(Symbolics.value.(substitute.(A_inv, (at,))))
pprint(round.(A_num; digits = 4), hydrates, anhydrous; label = :name)┌──────────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────┐
│ │ C │ M │ S │ A │ F │ S̄ │
├──────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┤
│ CâSAb̂Hĝ │ 0.2941 │ 0.2353 │ 0.5882 │ -1.1765 │ -1.7647 │ -0.1961 │
│ M5AH13 │ │ 0.2 │ │ │ │ │
│ C6AFS2H8 │ │ │ │ │ 1.0 │ │
│ C6AS̄3H32 │ │ │ │ │ │ 0.3333 │
│ C2ASH8 │ -0.2941 │ -0.2353 │ 0.4118 │ 1.1765 │ -0.2353 │ 0.1961 │
│ C4AH13 │ 0.2647 │ 0.0118 │ -0.4706 │ -0.0588 │ -0.5882 │ -0.5098 │
└──────────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────┘Why carry the symbols at all
Because the derivative is then free. How much the ettringite content responds to the C-S-H lime ratio is a partial derivative of an entry of this matrix, and with the symbol still in place it is obtained by differentiating rather than by re-solving at two nearby compositions:
∂ = Symbolics.derivative(A_inv[1, 1], â)
simplify(∂)What this matrix does not decide
It is worth being explicit, because the matrix is necessary and nowhere near sufficient. A fixes what is conservable: any composition n with A n = b respects the element balance. It says nothing about which of those compositions is the stable one — that is the Gibbs energy's business, and it needs the thermodynamic data the matrix knows nothing about. A stoichiometric matrix built from formulas alone, as on this page, can balance a reaction and cannot predict whether it proceeds.