Ageing creep of solidifying cementitious materials

The ageing-creep model of [56]: one phase solidifies progressively — as C-S-H does during hydration — so $\mathbb R^{\rm hom}(t, t')$ depends on the observation time $t$ and the loading time $t'$ independently. Laplace–Carson no longer applies; the homogenization runs directly in the time domain, through homogenize_alv.

The composite has three phase types:

Both viscoelastic phases obey a Maxwell relaxation law with separate bulk and shear characteristic times.

using MeanFieldHom
using TensND
using LinearAlgebra
using Plots
gr()  # headless backend; GKSwstype is set to "100" in make.jl

# Matrix
const E0, ν0, f0 = 1.0, 0.2, 0.6
const k0, μ0 = E0 / (3(1 - 2ν0)), E0 / (2(1 + ν0))
const η0, γ0 = 0.2, 0.133          # bulk / shear relaxation times
# Solidifying phase
const E1, ν1, finf = 5.0, 0.3, 0.3
const k1, μ1 = E1 / (3(1 - 2ν1)), E1 / (2(1 + ν1))
const η1, γ1 = 1.0, 1.67
# Pore (elastic, near-zero)
const Ep, νp = 1.0e-8, 0.2
const kp, μp = Ep / (3(1 - 2νp)), Ep / (2(1 + νp))
const fp = 1 - f0 - finf
const C_p = TensISO{3}(3kp, 2μp)

make_R0() = maxwell_iso(k0, μ0, η0, γ0)
make_R1() = maxwell_iso(k1, μ1, η1, γ1)

Solidification kinetics

The solidified fraction grows as $f(t) = f_\infty\, t^\alpha/(1+t^\alpha)$; the setting time of the layer carrying midpoint fraction $f_k = (k+\tfrac12) f_\infty/N$ is $t_k = (f_k/(f_\infty-f_k))^{1/\alpha}$.

function setting_times(N, α)
    F = [(i + 0.5) * finf / N for i in 0:(N - 1)]
    return [(f / (finf - f))^(1 / α) for f in F]
end

Per-layer relaxation law: history-dependent vs frozen

A newly formed layer is deposited stress-free and creeps only from its setting time on. History-dependent (fixed = false): layer $i$ responds as a solid only if it had set at the loading time $t'$. Frozen (fixed = true): the decision is made once, at the start of the observation window $t_0$ — a cheaper but physically approximate model.

function inclusion_law(t_set, t0; fixed)
    if fixed
        t0 ≥ t_set && return make_R1()
        return ViscoLaw((t, tp) -> (t < tp ? zero(C_p) : C_p), :relaxation)
    else
        R1 = make_R1()
        return ViscoLaw(
            function (t, tp)
                t < tp && return zero(C_p)
                tp ≥ t_set ? R1.eval_fun(t, tp) : C_p
            end, :relaxation)
    end
end

Two equivalent RVE topologies

[56]'s key contribution is that the $N$ solidifying shells and the pore can be packed into a single composite sphere instead of $N+1$ separate inclusions — reducing $N+1$ Eshelby problems to one. MeanFieldHom supports both: :whole_pores ($N$ separate spherical inclusions) and :layers (one LayeredSphere whose per-layer moduli are ageing relaxation laws).

function build_rve_whole_pores(N, α, t0; fixed)
    rve = RVE(:M)
    add_matrix!(rve, Ellipsoid(1.0), Dict(:C => make_R0()))
    add_phase!(rve, :PORE, Ellipsoid(1.0), Dict(:C => heaviside_law(C_p)); fraction = fp)
    t_sets = setting_times(N, α)
    for i in 1:N
        add_phase!(rve, Symbol(:INC_, i), Ellipsoid(1.0),
            Dict(:C => inclusion_law(t_sets[i], t0; fixed = fixed)); fraction = finf / N)
    end
    return rve
end

function build_rve_layers(N, α, t0; fixed)
    rve = RVE(:M)
    add_matrix!(rve, Ellipsoid(1.0), Dict(:C => make_R0()))
    t_sets = setting_times(N, α)
    f_layers = vcat([fp], fill(finf / N, N))        # pore innermost, shells outward
    cumulative = cumsum(f_layers)
    radii = ntuple(k -> cumulative[k]^(1 / 3), N + 1)
    moduli = ntuple(N + 1) do k
        k == 1 ? heaviside_law(C_p) : inclusion_law(t_sets[N - k + 2], t0; fixed = fixed)
    end
    sphere = LayeredSphere(radii, moduli)
    add_phase!(rve, :INCLUSION, sphere, Dict(:C => heaviside_law(C_p)); fraction = fp + finf)
    return rve
end

build_rve(N, α, t0, model; fixed) =
    model === :layers ? build_rve_layers(N, α, t0; fixed = fixed) :
    build_rve_whole_pores(N, α, t0; fixed = fixed)

Time-domain homogenization and effective creep

homogenize_alv returns the $6n\times6n$ block relaxation matrix over the time grid; its Volterra inverse (volterra_inverse) is the creep-compliance matrix, from which the uniaxial creep $E_0 J^E_{\rm eff}(t,t_0)$ follows.

function uniaxial_creep(R)
    J = volterra_inverse(R; block_size = 6)
    n = size(J, 1) ÷ 6
    return [sum(J[6(i - 1) + 1, 6(j - 1) + 1] for j in 1:n) for i in 1:n]
end

function creep_curve(N, α, t0, T, model; fixed)
    R = homogenize_alv(build_rve(N, α, t0, model; fixed = fixed), MoriTanaka(), :C; times = T)
    return uniaxial_creep(R)
end

Results

Following [56], the effective creep is computed for five loading ages $t_0$ (history-dependent, solid +; frozen, dashed) with both RVE topologies, side by side as in the Echoes book. N = 100 layers are used.

The elastic reference $1/E^{\rm hom}(t)$ (black dotted) is the instantaneous compliance, obtained from the same topology by evaluating the ageing homogenization over a single time (the relaxation at zero elapsed time is the glassy/elastic modulus). Because it uses the same morphology, every creep curve starts exactly on this reference at $t = t_0$.

const N, α_solid, t_max = 100, 4.0, 10 / 3
loading_ages = (1 / 3, 2 / 3, 4 / 3, 2.0, 8 / 3)
cmap = palette(:viridis, length(loading_ages))

# Instantaneous elastic compliance of a given topology at time t: one-time ALV.
elastic_ref(t, model) =
    uniaxial_creep(homogenize_alv(build_rve(N, α_solid, t, model; fixed = false),
        MoriTanaka(), :C; times = [t]))[1]

function creep_panel(model, title)
    p = plot(; xlabel = "t", ylabel = "E₀ · J^E_eff(t, t₀)", legend = :topleft,
        framestyle = :box, xlims = (0, t_max), ylims = (0, 20), title = title)
    for (k, t0) in enumerate(loading_ages)
        T = collect(range(t0, t_max; length = 31))
        Jh = creep_curve(N, α_solid, t0, T, model; fixed = false)
        Jf = creep_curve(N, α_solid, t0, T, model; fixed = true)
        plot!(p, T, E0 .* Jh; lw = 2, color = cmap[k], marker = :+, ms = 2,
            label = "history t₀=$(round(t0, digits = 2))")
        plot!(p, T, E0 .* Jf; lw = 1.5, color = cmap[k], ls = :dash, label = "")
    end
    T_ref = vcat([1e-3], filter(≤(t_max), setting_times(N, α_solid)), [t_max])
    plot!(p, T_ref, [elastic_ref(t, model) for t in T_ref]; lw = 2, color = :black,
        ls = :dot, label = "1/E^hom(t)")
    return p
end

plot(creep_panel(:layers, "model = :layers"),
    creep_panel(:whole_pores, "model = :whole_pores");
    layout = (2, 1), size = (760, 820))

Three observations, all reproducing [56]:

1. Ageing — early loading ages ($t_0$ small) give much larger creep because many layers have not yet solidified; the compliance decreases toward the elastic limit as $t_0$ grows.
2. History vs frozen — the frozen approach overestimates creep at early ages (it ignores solidification before $t_0$) and converges with the history-dependent result at late ages. Each curve starts on the dotted elastic reference of its own panel.
3. Morphology — the two panels do not coincide: :whole_pores is systematically more compliant than :layers.

A note on the two topologies

The :layers composite sphere and the :whole_pores collection of $N+1$ separate inclusions are different morphologies — the first packs the pore and the solidifying shells concentrically (as hydrates deposit around a pore), the second scatters them independently in the matrix.