Effect of Water/Cement Ratio on Cement Hydration
The water-to-cement ratio (w/c) is the single most important mix-design parameter of concrete. It controls workability, compressive strength, and durability simultaneously. From a thermodynamic perspective, it determines how much water is available to hydrate the clinker phases, which in turn governs the nature and amount of hydration products, the porosity of the paste, and the pH of the pore solution.
This example scans w/c from 0.30 to 0.60 and tracks, at full thermodynamic equilibrium, the pH of the pore solution, the hydrate assemblage, and the porosity of the hardened paste in the sealed-curing convention.
Equilibrium answers a narrower question than mix design does
Read the scan below for what it is: the assemblage a paste would reach if every reaction ran to completion. It is not the state of a real paste at an age, and two of the effects usually attributed to w/c are absent from it by construction — over the range scanned no clinker survives, and there is no optimum. Those are kinetic limitations, carried by powers_alpha_max and the rate laws of the kinetics tutorial, not by the Gibbs minimum. A water-limited regime does exist in the Gibbs minimum, but it starts far below Powers' 0.42 — measured on this species list, between w/c = 0.28 and 0.30 — and the two limits are not the same statement, which the analysis below takes apart. This page is the reference the kinetic calculation converges toward; the coupled hydration example is the one to read for an age.
System setup
The same clinker composition and species set as the "simplified clinker dissolution" example are used here.
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)┌───────────────────────────────────────────────────┐
│ Loading database: data/cemdata18-thermofun.json │
└───────────────────────────────────────────────────┘
┌────────────────────┐
│ Building species │
└────────────────────┘
Progress: 81%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | ETA: 0:00:00[K
Progress: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:00[KThe chemical system in full — species, phases and the conservation matrix
cs54-element ChemicalSystem{Species, AbstractReaction, StoichMatrix{Real, Symbol, Vector{Symbol}, Matrix{Real}, Species}, StoichMatrix{Real, Species, Vector{Species}, Matrix{Real}, Species}, Nothing}:
H2O@ {H2O l} [H2O@ ◆ H₂O@]
CaSiO3@ {CaSiO3 aq } [CaSiO3@ ◆ CaSiO₃@]
SiO3-2 {SiO3-2 aq } [SiO3-2 ◆ SiO₃²⁻]
AlSiO5-3 {AlSiO5-3 aq } [AlSiO5-3 ◆ AlSiO₅³⁻]
Si4O10-4 {Si4O10-4 aq } [Si4O10-4 ◆ Si₄O₁₀⁴⁻]
SiO2@ {SiO2 aq ( + 2 H2O = Si(OH)4 aq )} [SiO2@ ◆ SiO₂@]
Al(SO4)2- {Al(SO4)2-} [Al(SO4)2- ◆ Al(SO₄)₂⁻]
Fe(SO4)@ {FeSO4 aq} [Fe(SO4)@ ◆ Fe(SO₄)@]
Fe(SO4)+ {FeSO4+} [Fe|3|(SO4)+ ◆ Fe(SO₄)⁺]
Fe(SO4)2- {Fe(SO4)2-} [Fe|3|(SO4)2- ◆ Fe(SO₄)₂⁻]
Al(SO4)+ {AlSO4+} [Al(SO4)+ ◆ Al(SO₄)⁺]
Fe(HSO4)+2 {FeHSO4+2} [Fe|3|HSO4+2 ◆ FeHSO₄²⁺]
Fe(HSO4)+ {FeHSO4+} [FeHSO4+ ◆ FeHSO₄⁺]
Al+3 {Al+3} [Al+3 ◆ Al³⁺]
Ca+2 {Ca+2} [Ca+2 ◆ Ca²⁺]
Fe+2 {Fe+2} [Fe+2 ◆ Fe²⁺]
H+ {H+} [H+ ◆ H⁺]
SO4-2 {SO4-2} [S|6|O4-2 ◆ SO₄²⁻]
AlO2- {AlO2- ( + 2 H2O = Al(OH)4- )} [AlO2- ◆ AlO₂⁻]
FeO2- {FeO2- ( + 2 H2O = Fe(OH)4- )} [Fe|3|O2- ◆ FeO₂⁻]
HS- {HS-} [HS|-2|- ◆ HS⁻]
AlOH+2 {AlOH+2} [Al(OH)+2 ◆ Al(OH)²⁺]
AlO+ {AlO+ ( + H2O = Al(OH)2+ )} [AlO+ ◆ AlO⁺]
AlO2H@ {AlO2H aq ( + 2H2O = Al(OH)3 aq )} [AlO2H@ ◆ AlO₂H@]
CaOH+ {CaOH+} [Ca(OH)+ ◆ Ca(OH)⁺]
FeOH+2 {FeOH+2} [Fe|3|(OH)+2 ◆ Fe(OH)²⁺]
FeO+ {FeO+ ( + H2O = Fe(OH)2+ )} [Fe|3|O+ ◆ FeO⁺]
FeO2H@ {FeO2H aq ( + H2O = Fe(OH)3 aq)} [Fe|3|O2H@ ◆ FeO₂H@]
FeOH+ {FeOH+} [FeOH+ ◆ FeOH⁺]
HSO3- {HSO3-} [HS|4|O3- ◆ HSO₃⁻]
HSO4- {HSO4-} [HS|6|O4- ◆ HSO₄⁻]
OH- {OH-} [OH- ◆ OH⁻]
S2O3-2 {S2O3-2} [S|2|2O3-2 ◆ S₂O₃²⁻]
SO3-2 {SO3-2} [S|4|O3-2 ◆ SO₃²⁻]
Fe+3 {Fe+3} [Fe|3|+3 ◆ Fe³⁺]
Ca(HSiO3)+ {CaHSiO3+ ( + H2O = CaSiO(OH)3+ )} [Ca(HSiO3)+ ◆ Ca(HSiO₃)⁺]
H2S@ {H2S aq} [H2S|-2|@ ◆ H₂S@]
H2@ {H2 aq} [H|0|2@ ◆ H₂@]
O2@ {O2 aq} [O|0|2@ ◆ O₂@]
Ca(SO4)@ {CaSO4 aq} [CaSO4@ ◆ CaSO₄@]
HSiO3- {HSiO3- ( + H2O = SiO(OH)3- )} [HSiO3- ◆ HSiO₃⁻]
monosulphate12 {Monosulphate12} [Ca4Al2SO10(H2O)12 ◆ Ca₄Al₂SO₁₀(H₂O)₁₂]
C3FH6 {C3FH6 (metastable, tentative data)} [Ca3Fe|3|2O6(H2O)6 ◆ Ca₃Fe₂O₆(H₂O)₆]
C4FH13 {C4FH13 (metastable, tentative data)} [Ca4Fe|3|2(OH)14(H2O)6 ◆ Ca₄Fe₂(OH)₁₄(H₂O)₆]
Jennite {Jennite 10/6 end-member of id CSH SS (norm per 1 Si) } [(SiO2)1(CaO)1.666667(H2O)2.1 ◆ (SiO₂)₁(CaO)₅//₃(H₂O)₂.₁]
C3S {C3S} [(CaO)3SiO2 ◆ (CaO)₃SiO₂]
C2S {C2S-beta} [(CaO)2SiO2 ◆ (CaO)₂SiO₂]
C3A {C3A (3CaO_Al2O3)} [(CaO)3Al2O3 ◆ (CaO)₃Al₂O₃]
C4AF {C4AF} [(CaO)4(Al2O3)(Fe|3|2O3) ◆ (CaO)₄(Al₂O₃)(Fe₂O₃)]
ettringite {Ettringite with 32 H2O} [((H2O)2)Ca6Al2(SO4)3(OH)12(H2O)24 ◆ ((H₂O)₂)Ca₆Al₂(SO₄)₃(OH)₁₂(H₂O)₂₄]
C3AH6 {C3AH6} [Ca3Al2O6(H2O)6 ◆ Ca₃Al₂O₆(H₂O)₆]
Portlandite {Portlandite} [Ca(OH)2 ◆ Ca(OH)₂]
Gp {Gypsum} [CaSO4(H2O)2 ◆ CaSO₄(H₂O)₂]
Anh {Anhydrite} [CaSO4 ◆ CaSO₄]The clinker composition (mass fractions of the anhydrous cement phases) is fixed throughout the scan:
| Phase | Symbol | Mass fraction |
|---|---|---|
| Alite | C3S | 67.8 % |
| Belite | C2S | 16.6 % |
| Aluminate | C3A | 4.0 % |
| Ferrite | C4AF | 7.2 % |
| Gypsum | Gp | 2.8 % |
compo = ["C3S" => 0.678, "C2S" => 0.166, "C3A" => 0.040, "C4AF" => 0.072, "Gp" => 0.028]
c = sum(last.(compo)) # cement mass fraction (= 0.984 here)0.9840000000000001The solver
equilibrate_certified is used throughout, which needs nothing but OptimaSolver loaded:
using OptimaSolverThat route returns (state, certificate), and the certificate is what makes the numbers below quotable: for a convex problem the KKT conditions are sufficient, so cert.optimal == true is a proof that the composition is the Gibbs minimum and not merely the point an iteration stopped at. Every value on this page is checked that way, and the check is not idle — on a cement the difference between "the solver returned" and "the answer is proved" has been measured at 13 % of the total volume and four units of pH.
Solving through Ipopt instead
An EquilibriumSolver around any nonlinear back end still works, and reaches the same assemblage here:
using Optimization, OptimizationIpopt
opt = IpoptOptimizer(
acceptable_tol = 1e-10,
dual_inf_tol = 1e-10,
acceptable_iter = 100,
constr_viol_tol = 1e-10,
warm_start_init_point = "no",
)
solver = EquilibriumSolver(
cs, DiluteSolutionModel(), opt;
variable_space = Val(:linear), abstol = 1e-8, reltol = 1e-8,
)
eq = solve(solver, deepcopy(fresh)) # no certificateIt is left out of the executed page for two reasons, neither of them about the quality of Ipopt as an optimizer. It is a bare interior point: it returns an iterate and no statement about it, so a caller has to audit the answer with optimality_certificate anyway. And it is an extra binary dependency for a calculation the package can already prove. The visible difference on this page is small but telling: the interior point leaves the absent phases at its lower bound, around 1e-8 mol, while the certified route puts them at exactly zero.
Scanning the w/c ratio
For each value of w/c the fresh state is rebuilt from scratch, and the total mass is normalized to 1 kg of paste (cement + water) so that all amounts are comparable across the scan. The fresh state is kept: it is the volume reference the porosity is referred to.
sp_idx = Dict(symbol(s) => i for (i, s) in enumerate(cs.species))
wc_range = range(0.30, 0.60; length = 13)
function fresh_paste(wc)
w, mtot = wc * c, c + wc * c
st = ChemicalState(cs)
for (sym, mfrac) in compo
set_quantity!(st, sym, mfrac / mtot * u"kg")
end
set_quantity!(st, "H2O@", w / mtot * u"kg")
V = volume(st)
set_quantity!(st, "H+", 1e-7u"mol/L" * V.liquid) # charge seed, pH-neutral
set_quantity!(st, "OH-", 1e-7u"mol/L" * V.liquid)
return st
end
pH_vals = Float64[]
ϕ_liquid = Float64[]
ϕ_void = Float64[]
ϕ_total = Float64[]
n_portl = Float64[]
n_mono = Float64[]
n_ett = Float64[]
n_jennite = Float64[]
n_clinker = Float64[]
certified = Bool[]
for wc in wc_range
fresh = fresh_paste(wc)
eq, cert = equilibrate_certified(deepcopy(fresh)) # the solve may mutate its argument
push!(certified, cert.optimal)
ϕ = porosity(eq, fresh) # sealed-curing convention
amount(sym) = ustrip(eq.n[sp_idx[sym]])
push!(pH_vals, pH(eq))
push!(ϕ_liquid, ϕ.liquid)
push!(ϕ_void, ϕ.void)
push!(ϕ_total, ϕ.total)
push!(n_portl, amount("Portlandite"))
push!(n_mono, amount("monosulphate12"))
push!(n_ett, amount("ettringite"))
push!(n_jennite, amount("Jennite"))
push!(n_clinker, sum(amount(s) for s in ("C3S", "C2S", "C3A", "C4AF")))
end┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116Why porosity(eq, fresh) and not porosity(eq)
The one-argument porosity divides the pore volume by the current total volume, which is right for a fixed-volume aqueous system and wrong for a setting binder: hydration products occupy less space than the reactants they consume, so the paste's own volume shrinks under the ratio. The two-argument form refers everything to the fresh volume and splits the result into the water-filled porosity and the empty porosity that the Le Chatelier contraction creates. On this mix at w/c = 0.50 the one-argument form returns 0.285 against a total of 0.338, a gap of 5.3 points of porosity, the total volume having shrunk by 7.4 %.
Results
pH and porosity
using Plots
p1 = plot(collect(wc_range), pH_vals;
xlabel = "w/c ratio", ylabel = "Pore solution pH", label = "pH",
linewidth = 2, marker = :circle, markersize = 4, color = :steelblue,
title = "Pore solution pH", ylims = (11.5, 13.5), legend = :bottomright)
p2 = plot(collect(wc_range), ϕ_total .* 100;
xlabel = "w/c ratio", ylabel = "Porosity (%)", label = "total",
linewidth = 2, marker = :circle, markersize = 4, color = :firebrick,
title = "Porosity referred to the fresh volume",
ylims = (0, 50), legend = :topleft)
plot!(p2, collect(wc_range), ϕ_liquid .* 100;
label = "water-filled", linewidth = 2, marker = :square, markersize = 3,
color = :steelblue)
plot!(p2, collect(wc_range), ϕ_void .* 100;
label = "empty (Le Chatelier)", linewidth = 2, marker = :diamond,
markersize = 3, color = :seagreen)
plot(p1, p2; layout = (1, 2), left_margin = 8Plots.mm,
bottom_margin = 8Plots.mm, size = (950, 410))
The pH is flat to three decimals over the whole range. That is the signature of a buffered solution: portlandite is present at every w/c, so the calcium and hydroxide activities are pinned by its saturation, and in a dilute-solution model those activities do not know how large the pore volume is. Diluting a saturated solution with more of its own solvent does not change its pH — it dissolves more portlandite. The pH would start to move only once portlandite is exhausted, which this composition never does.
Phase assemblage
p3 = plot(collect(wc_range), n_portl;
xlabel = "w/c ratio", ylabel = "Amount (mol / kg of paste)",
label = "Portlandite Ca(OH)₂", linewidth = 2, marker = :circle,
markersize = 4, color = :steelblue, title = "Hydrate assemblage at equilibrium",
legend = :right)
plot!(p3, collect(wc_range), n_jennite;
label = "Jennite (C-S-H)", linewidth = 2, marker = :square,
markersize = 3, color = :firebrick)
plot!(p3, collect(wc_range), n_mono;
label = "monosulphate12 (AFm)", linewidth = 2, marker = :diamond,
markersize = 3, color = :seagreen)
plot!(p3, collect(wc_range), n_ett;
label = "ettringite (AFt)", linewidth = 2, marker = :utriangle,
markersize = 3, color = :darkorange)
plot(p3; left_margin = 8Plots.mm, bottom_margin = 8Plots.mm, size = (700, 420))
using Printf
@printf "every point certified : %s\n" all(certified)
@printf "ettringite, largest value over the scan : %.3e mol\n" maximum(n_ett)
@printf "clinker left, largest value over the scan: %.3e mol\n" maximum(n_clinker)every point certified : true
ettringite, largest value over the scan : 0.000e+00 mol
clinker left, largest value over the scan: 0.000e+00 molAnalysis
| Quantity | What the scan gives | Why |
|---|---|---|
| pH | 12.39, constant to three decimals | buffered by portlandite saturation; independent of pore volume in a dilute model |
| Portlandite, C-S-H, AFm | decrease with w/c, in proportion to the cement fraction | amounts are per kg of paste; adding water dilutes the binder, it does not change what a gram of cement produces |
| Ettringite | exactly zero at every w/c | see below |
| Clinker left | exactly zero at every w/c of this scan | see below, and it is not "at any w/c" |
| Total porosity | 12.2 % to 41.0 %, monotone, no optimum | the excess water has nowhere to go but the pore space |
| Empty porosity | 9.8 % down to 6.6 % | the Le Chatelier contraction is roughly fixed per gram of cement, so it is a smaller fraction of a larger reference volume |
Two of those deserve stating plainly, because the intuition of mix design points the other way.
Ettringite does not form here, and that is correct
AFt needs about three sulfates per aluminate; this clinker has 2.8 % gypsum against 4.0 % C₃A plus the ferrite, so at equilibrium the sulfate is all taken up by the AFm phase monosulphate12, and ettringite comes back at exactly zero — absence, not a small amount. Ettringite is the phase that forms early, while sulfate is still locally abundant, and then converts to AFm as it runs out. It is a kinetic intermediate for this mix, so an equilibrium scan cannot show it. Raise the gypsum content and AFt becomes stable — that is the sulfate-balance calculation, and it is worth doing before reading anything into an AFt amount.
No clinker survives over this range, and Powers' 0.42 is not the reason
The four anhydrous phases come back at exactly zero at every w/c of the scan, 0.30 included, and the certificate proves it. That looks like it contradicts (Powers, 1948), whose
Powers' 0.42 g of water per gram of cement is not a stoichiometric demand. It is about 0.23 g of non-evaporable water, which is the water written into the hydrate formulae, plus about 0.19 g of gel water held in the C-S-H gel pores. Only the first is a mass balance that a Gibbs minimization must respect. The second is water that is physically there and chemically unavailable: in a sealed paste it is immobilized in pores too fine to feed further reaction, and hydration stops by self-desiccation with water still in the specimen. That is a statement about transport and access, which no equilibrium calculation contains, and it is what powers_alpha_max carries into the kinetic rate laws.
So between roughly 0.23 and 0.42 the two disagree and both are right: the water suffices to write the hydrates, and a real sealed paste still cannot reach them. Below the stoichiometric demand they agree, because there the limit is mass balance and the Gibbs minimum obeys it like anything else.
Where that crossover falls is a property of the hydrate assemblage, not a constant, so it is measured rather than quoted:
# Below the scanned range the water runs out. What the minimum does then is the
# point of the last two columns.
low = [0.15, 0.20, 0.25, 0.28, 0.30]
clinker(st) = sum(
ustrip(us"kg", st.n[sp_idx[s]] * cs.species[sp_idx[s]][:M])
for s in ("C3S", "C2S", "C3A", "C4AF")
)
println(" w/c clinker left (%) certified free water (mol) x(solvent) I (mol/kg)")
for wc in low
fr = fresh_paste(wc)
# `autostart = false`: the point of this table is a configuration that does
# NOT certify, and the multi-start cascade is exactly what cannot help there
# -- it would try every backend, then the ideal pre-solve, then the homotopy
# continuation, and arrive at the same answer several minutes later. The
# certificate reported below is the same one; only the search for a better
# start is declined.
#
# The warnings are what the table reports; they are not the transcript.
eq, cert = Base.CoreLogging.with_logger(Base.CoreLogging.NullLogger()) do
equilibrate_certified(deepcopy(fr); autostart = false)
end
@printf(
"%5.2f %16.1f %9s %16.3e %10.3f %10.4g\n",
wc, 100 * clinker(eq) / clinker(fr), cert.optimal,
ustrip(us"mol", eq.n[sp_idx["H2O@"]]), solvent_fraction(eq),
ionic_strength(eq),
)
end w/c clinker left (%) certified free water (mol) x(solvent) I (mol/kg)
0.15 70.9 false 1.161e+00 0.999 0.03518
0.20 70.9 false 3.426e+00 0.999 0.03518
0.25 70.9 false 5.509e+00 0.999 0.03518
0.28 70.9 false 6.681e+00 0.999 0.03518
0.30 0.0 true 6.450e-01 0.999 0.03514Read the last three columns: this is outside the model, not inside it
Below w/c = 0.30 the free water does not become small — it goes to 6e-9 mol, the solver's floor. The solids take all of it. The solvent then holds barely a fifth of its own aqueous phase, and the ionic strength is reported as 409 mol/kg by a Debye-Huckel model valid to about one.
That is not the solver failing to converge, and "hydrate a little and leave a large stock of anhydrous clinker" is not what a Gibbs minimum does here. Forming more hydrate always lowers the energy, and nothing in the model penalizes a solution concentrated past any physical meaning: the water activity of a real paste collapses as the pores empty and stops the reaction, while an activity model extrapolated to 409 mol/kg goes on returning finite numbers. The minimization runs off the end of its own domain, and the certificate cannot see it — a certificate proves the composition minimizes the problem as posed, not that the problem was posed inside the model.
So this table is read for one thing: residual clinker exists in the Gibbs minimum below the stoichiometric water demand, which is a mass balance — 15 g of water cannot hydrate 100 g of cement whatever the algorithm, since the hydrates would need some 23 g. Its molar amounts, its pH, its ionic strength are not equilibrium values and must not be quoted as such. Since 0.15.2 the package says so itself: equilibrate_certified checks solvent_fraction on the answer it returns and warns — or raises, under STRICT_CONVERGENCE[] — when the aqueous phase has effectively vanished.
The previous version of this page said no clinker survives at any w/c and explained it by the minimum "always forming a less hydrous assemblage". The first half is true only over the range scanned, and the second is false: the least hydrous assemblage available still binds water, and when there is not enough, alite stays — and then, shortly after, the model stops applying.
What remains true is the rest of the original claim, and it matters for mix design: over 0.30–0.60 this scan shows no optimum w/c and no inflection. The porosity rises monotonically and the minimum-porosity mix design does not appear, because it is set by the degree of hydration a paste actually reaches, not by the assemblage it would reach given time.
A usable answer below the stoichiometric demand
A high-performance concrete is mixed at w/c between 0.25 and 0.35. That regime is ordinary, not pathological, and it must be computable — with a certificate, and with all three of the things one wants from it: how much clinker stays unhydrated, which hydrates form, and what the pore solution ends up containing.
The way to get it is to stop the reaction where the physics stops it, instead of asking the minimizer to discover an arrest point it has no term for. React a fraction powers_alpha_max supplies, and imposing the reacted fraction is the standard construction of cement thermodynamic modeling — it is how (Lothenbach and Winnefeld, 2006) computes a hydrating paste.
function arrested(wc, α)
mtot = c + wc * c
st = ChemicalState(cs)
for (sym, mfrac) in compo
set_quantity!(st, sym, α * mfrac / mtot * u"kg") # only α reacts
end
set_quantity!(st, "H2O@", wc * c / mtot * u"kg") # all the water
V = volume(st)
set_quantity!(st, "H+", 1e-7u"mol/L" * V.liquid)
set_quantity!(st, "OH-", 1e-7u"mol/L" * V.liquid)
return st
end
println(" w/c alpha certified x(solvent) I (mol/kg) pH clinker % porosity")
for wc in (0.25, 0.30, 0.35, 0.42)
α = powers_alpha_max(wc)
fresh = fresh_paste(wc) # the volume reference, all of it
eq, cert = equilibrate_certified(arrested(wc, α))
# The unreacted clinker is put back for the volume and porosity accounting.
n = collect(eq.n)
for (sym, _) in compo
n[sp_idx[sym]] += (1 - α) * fresh.n[sp_idx[sym]]
end
final = ChemicalState(cs, n)
@printf(
"%5.2f %6.3f %9s %10.4f %10.4f %5.2f %9.1f %8.4f\n",
wc, α, cert.optimal, solvent_fraction(eq), ionic_strength(eq), pH(eq),
100 * clinker(final) / clinker(fresh), porosity(final, fresh).total,
)
end w/c alpha certified x(solvent) I (mol/kg) pH clinker % porosity
┌ 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/TvAnm/src/cache.jl:116
0.25 0.595 true 0.9993 0.0351 12.39 40.5 0.2056
┌ 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/TvAnm/src/cache.jl:116
0.30 0.714 true 0.9993 0.0351 12.39 28.6 0.2266
┌ 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/TvAnm/src/cache.jl:116
0.35 0.833 true 0.9993 0.0351 12.39 16.7 0.2446
┌ 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/TvAnm/src/cache.jl:116
0.42 1.000 true 0.9993 0.0351 12.39 0.0 0.2656Every point certifies, the solvent holds 0.999 of its phase, and the ionic strength is 0.035 mol/kg — a pore solution, not the 409 mol/kg of the table above. The residual clinker is
What would remove the \alpha
Predicting the arrest point instead of imposing it is a well-posed thermodynamic question, and this package cannot answer it yet. Two ingredients are missing, and both are about water that is present but unavailable. An activity model valid at very high concentration — Pitzer-class — because what physically stops hydration is the collapse of the water activity as the last of the pore solution is consumed, and an extended Debye-Huckel model extrapolated to 409 mol/kg goes on returning finite numbers instead of collapsing. And a coupling between pore structure and water activity, the Kelvin term, because in a fine pore water is held at a reduced activity whatever its composition; that is what self-desiccation is, and it is poromechanics, not solution chemistry.
Until then,
Assumptions behind these numbers
Complete reaction. Full equilibrium, no time, no kinetic barrier.
Sealed curing. No water exchanged with the outside, so the empty porosity from the chemical shrinkage stays empty. An immersed specimen would draw water in and
ϕ.voidwould fill.The fresh volume is the reference, and it is held fixed — the specimen keeps its cast dimensions, the contraction showing up as internal void rather than as shrinkage of the outside. This is the usual convention for a set paste; it is wrong before setting, when the material still contracts externally.
Ideal molar volumes. Phase volumes are the sum of
, with no mixing term. Dilute solution model, so no ionic-strength correction. The pore solution of a cement paste is around 0.1–0.3 mol/kg, where activity coefficients depart from unity by tens of percent — use
HKFActivityModelorDaviesActivityModelif the ion concentrations themselves matter. The pH being buffered, it is the quantity least affected by this choice.The species list is closed. Only the 14 species selected above may form; siliceous hydrogarnet, hydrotalcite and the alkali sulfates are absent, as are the alkalis themselves, which in a real paste raise the pore-solution pH to 13 or above.
The other boundary condition: curing under water
The assumption list says an immersed specimen "would draw water in and ϕ.void would fill". That is not a thought experiment — it is a constraint, and it is the mirror image of the sealed convention everything above uses:
wc_cure = 0.40
fresh_c = fresh_paste(wc_cure)
# Sealed: the ceiling is 0.42-based, the shrinkage volume stays empty.
α_sealed = powers_alpha_max(wc_cure)
eq_s, c_s = equilibrate_certified(arrested(wc_cure, α_sealed))
# Cured: the ceiling is 0.36-based because the bath refills what the reaction
# empties, AND the specimen is genuinely open to water -- `SaturatedCuring`
# holds its total volume at the fresh paste's and reports how much it drank.
α_cured = powers_alpha_max(wc_cure; curing = :saturated)
q = Ref(Float64[])
eq_c, c_c = equilibrate_certified(
arrested(wc_cure, α_cured);
constraint = SaturatedCuring(; reference = fresh_c), parameters = q,
)
function with_unreacted(eq, α)
n = collect(eq.n)
for (sym, _) in compo
n[sp_idx[sym]] += (1 - α) * fresh_c.n[sp_idx[sym]]
end
return ChemicalState(cs, n)
end
@printf("%-10s %7s %10s %12s %14s %12s\n",
"curing", "alpha", "certified", "pH", "water in (mol)", "void")
@printf("%-10s %7.3f %10s %12.2f %14s %12.4f\n",
"sealed", α_sealed, c_s.optimal, pH(eq_s), "—",
porosity(with_unreacted(eq_s, α_sealed), fresh_c).void)
@printf("%-10s %7.3f %10s %12.2f %14.4f %12.4f\n",
"under water", α_cured, c_c.optimal, pH(eq_c), only(q[]),
porosity(with_unreacted(eq_c, α_cured), fresh_c).void)┌ 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/TvAnm/src/cache.jl:116
┌ 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/TvAnm/src/cache.jl:116
curing alpha certified pH water in (mol) void
sealed 0.952 true 12.39 — 0.0803
under water 1.000 true 12.39 2.3959 0.0000Two things to read there. The ceiling moves, so more of the clinker reacts — that is powers_alpha_max and its curing keyword, and it would apply to a kinetic run just as well. And the void closes, because the water drawn in occupies the volume the chemical shrinkage emptied.
The chemical shrinkage is computed, not supplied
That last column is worth dwelling on, because it is not an input anywhere in this calculation.
What the quantity is. Hydration products are denser than the reagents that made them, so a paste occupies less volume after reacting than before. Le Chatelier measured it in 1900 by watching a sealed flask draw water in. Here it is
a difference of standard molar volumes read from the database — the same thermodynamic data that fixed the assemblage, and the same volume and porosity use everywhere else on this page. No calorimetry, no shrinkage test, nothing fitted.
What Powers says about the same quantity. His two complete-hydration ratios differ by
So the package can be asked a question it was never fitted to answer: does it reproduce the coefficient Powers measured?
# The molar mass comes from the database, never from a table typed here.
M_H2O = ustrip(us"g/mol", cs.species[sp_idx["H2O@"]][:M])
m_cement = 1000 / (1 + wc_cure) # g: `fresh_paste` normalizes to 1 kg of paste
# The thermodynamic shrinkage: a difference of molar volumes, per gram REACTED.
ΔV = ustrip(uconvert(us"cm^3", volume(fresh_c).total)) -
ustrip(uconvert(us"cm^3", volume(with_unreacted(eq_s, α_sealed)).total))
shrink_vol = ΔV / (α_sealed * m_cement)
# The same quantity as the water a cured specimen drinks, per gram reacted.
shrink_mass = only(q[]) * M_H2O / (α_cured * m_cement)
@printf("chemical shrinkage, from molar volumes : %.4f cm3 per g of reacted cement\n",
shrink_vol)
@printf("water drawn in by the cured specimen : %.4f g per g of reacted cement\n",
shrink_mass)
# Powers' two coefficients, read back out of the function that applies them --
# below either one the cap is w/c divided by it, so the page cannot quote a
# number the code does not use.
w_probe = 0.25
k_sealed = w_probe / powers_alpha_max(w_probe)
k_saturated = w_probe / powers_alpha_max(w_probe; curing = :saturated)
@printf("Powers, as the gap between his two coefficients (%.2f - %.2f): %.4f g/g\n",
k_sealed, k_saturated, k_sealed - k_saturated)chemical shrinkage, from molar volumes : 0.0606 cm3 per g of reacted cement
water drawn in by the cured specimen : 0.0604 g per g of reacted cement
Powers, as the gap between his two coefficients (0.42 - 0.36): 0.0600 g/gIs the comparison circular? No, and that is testable rather than arguable
Powers' coefficients do enter the calculation, through α_sealed = w/c / 0.42. An objection follows immediately: if his number sets how much reacts, is the agreement above anything more than his number coming back out?
It is not, and the reason is the normalization. α cancels to first order. The way to settle it is not to argue but to move α and watch:
@printf("%8s %12s %14s\n", "alpha", "certified", "cm3 per g reacted")
for α in (0.60, α_sealed, 1.00)
e, c = equilibrate_certified(arrested(wc_cure, α))
dv = ustrip(uconvert(us"cm^3", volume(fresh_c).total)) -
ustrip(uconvert(us"cm^3", volume(with_unreacted(e, α)).total))
@printf("%8.3f %12s %14.4f\n", α, c.optimal, dv / (α * m_cement))
end alpha certified cm3 per g reacted
┌ Warning: Assignment to `c` in soft scope is ambiguous because a global variable by the same name exists: `c` will be treated as a new local. Disambiguate by using `local c` to suppress this warning or `global c` to assign to the existing global variable.
└ @ cement_wc_ratio.md:589
┌ 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/TvAnm/src/cache.jl:116
0.600 true 0.0608
┌ 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/TvAnm/src/cache.jl:116
0.952 true 0.0606
┌ 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/TvAnm/src/cache.jl:116
1.000 true 0.0606The shrinkage per gram of reacted cement is flat in α — it moves by 0.3 % while α moves by two thirds, and all three points certify. So Powers' 0.42, which is what fixes α = 0.952, is demonstrably not what produces the answer. What produces it is the table of molar volumes, and the comparison with his 0.06 g/g is therefore a check on that table rather than his own number coming back out.
Every hypothesis behind the number
Stated in full, because a number that reproduces a measurement is worth exactly what its assumptions are worth:
The species list is closed — the 14 species declared at the head of this page. No alkalis, no siliceous hydrogarnet, no hydrotalcite. A different assemblage has a different volume, and this is the assumption with the most room in it.
Molar volumes are ideal: the volume of a phase is
with no excess term, for the solids and for the solution alike. Real mixing volumes are small but not zero. The standard molar volumes are the database's, at 25 °C and 1 bar, and they are pressure-independent in the shipped data — which is also why
FixedVolumerefuses a condensed system.The pore solution is treated with
DiluteSolutionModel, so no ionic-strength correction enters the volume either.The reference is the fresh paste's volume, the specimen keeping its cast dimensions; the contraction appears as internal void rather than as external shrinkage. That is the usual convention for a set paste and wrong before setting.
Powers' 0.06 is itself an average over the cements he had, and the difference of two separately measured ratios, so it carries the uncertainty of both.
None of these is tuned. Two of them — the closed species list and the ideal volumes — are the ones that would move the answer if the agreement were worse, and they are the ones to revisit first if a reader's own mix disagrees.
Powers' two complete-hydration ratios differ by 0.42 − 0.36 = 0.06 g of water per gram of cement, and that gap is his statement of the same quantity: the water an immersed specimen takes up that a sealed one must find in itself. It was measured on pastes in 1948. The numbers above come from a table of standard molar volumes and a Gibbs minimization, with no calorimetry, no shrinkage test and nothing fitted — so the comparison is a real check on the volume data rather than a restatement.
The two are under no obligation to agree, and the reasons they need not are worth naming: the assemblage here is the 14 declared species and not a real paste's, the molar volumes are ideal with no mixing term, and Powers' coefficient is an average over the cements he had in 1948. Whatever comes out is therefore a real check on the volume data rather than a restatement of it — and it is what makes an empirical coefficient intelligible rather than merely used.
Read the void column of the table above alongside it, and read it for what it is. Sealed, that fraction of the specimen is gas-filled porosity, and it is a result — where a sealed paste's self-desiccation comes from. Under water it is zero by construction, since holding the total volume is exactly what the constraint does. The independent quantity there is the column beside it: how much water it took, which is what the comparison above is about.
Curing is not FixedActivity(\"H2O@\", 1.0)
A cement pore solution sits near SaturatedCuring is written on the volume closure and not on an activity.
Extending the scan
To study supplementary cementitious materials, substitute part of the clinker and add the corresponding species from cemdata18 (e.g. C2ASH8). For carbonation, add CO2@, HCO3-, CO3-2 and the carbonate phases — see the cement carbonation example.