MeanFieldHomogenization.jl

Appearance

Ageing creep of solidifying cementitious materials ​

The ageing-creep model of [60]: one phase solidifies progressively — as C-S-H does during hydration — so depends on the observation time and the loading time independently. Laplace–Carson no longer applies; the homogenization runs directly in the time domain, through homogenize_alv.

The composite has three phase types:

PhaseFractionStiffnessRheology
Matrix   Maxwell
Solidifying inclusions   Maxwell (per-layer setting time)
Pore   elastic

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

julia
using MeanFieldHomogenization
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    ; the setting time of the layer carrying midpoint fraction    is   .

julia
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 responds as a solid only if it had set at the loading time . Frozen (fixed = true): the decision is made once, at the start of the observation window — a cheaper but physically approximate model.

julia
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 ​

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

julia
function build_rve_whole_pores(N, α, t0; fixed)
    rve = RVE()
    add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => make_R0()); fraction = :rest)
    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()
    add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => make_R0()); fraction = :rest)
    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   block relaxation matrix over the time grid; its Volterra inverse (volterra_inverse) is the creep-compliance matrix, from which the uniaxial creep follows.

julia
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 [60], the effective creep is computed for five loading ages (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 (black dotted) is the instantaneous (glassy) compliance of the microstructure frozen at time : every layer past its setting time carries its elastic stiffness , the others are still pores. It is built from the same topology but computed by the purely elastic pipeline — homogenize, no Volterra algebra anywhere — so its agreement with the start of each creep curve is an independent cross-check of homogenize_alv, not a restatement of it.

julia
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))
const C_0_el, C_1_el = TensISO{3}(3k0, 2μ0), TensISO{3}(3k1, 2μ1)

# Elastic RVE frozen at time t — same two topologies, elastic moduli only.
layer_stiffness(t, t_set) = t ≥ t_set ? C_1_el : C_p

function build_elastic_rve(N, α, t, model)
    rve = RVE()
    add_phase!(rve, :M, Ellipsoid(1.0), Dict(:C => C_0_el); fraction = :rest)
    t_sets = setting_times(N, α)
    if model === :layers
        cumulative = cumsum(vcat([fp], fill(finf / N, N)))
        radii = ntuple(k -> cumulative[k]^(1 / 3), N + 1)
        moduli = ntuple(k -> k == 1 ? C_p : layer_stiffness(t, t_sets[N - k + 2]), N + 1)
        add_phase!(rve, :INCLUSION, LayeredSphere(radii, moduli), Dict(:C => C_p);
            fraction = fp + finf)
    else
        add_phase!(rve, :PORE, Ellipsoid(1.0), Dict(:C => C_p); fraction = fp)
        for i in 1:N
            add_phase!(rve, Symbol(:INC_, i), Ellipsoid(1.0),
                Dict(:C => layer_stiffness(t, t_sets[i])); fraction = finf / N)
        end
    end
    return rve
end

function elastic_ref(t, model)
    C_hom = homogenize(build_elastic_rve(N, α_solid, t, model), MoriTanaka(), :C)
    K_hom, μ_hom = TensND.get_data(C_hom)[1] / 3, TensND.get_data(C_hom)[2] / 2
    return E0 / max(9K_hom * μ_hom / (3K_hom + μ_hom), 1.0e-12)
end

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
    # Sample the setting times (where 1/E^hom jumps) *and* the loading ages,
    # so the start of each curve is read off the reference without
    # interpolation error.
    T_ref = sort(vcat([1e-3], filter(≤(t_max), setting_times(N, α_solid)),
        collect(loading_ages), [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), left_margin = 8Plots.mm, bottom_margin = 8Plots.mm, size = (860, 880))

Three observations, all reproducing [60]:

  1. Ageing — early loading ages ( small) give much larger creep because many layers have not yet solidified; the compliance decreases toward the elastic limit as grows.

  2. History vs frozen — the frozen approach overestimates creep at early ages (it ignores solidification before ) 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.

Elastic cross-check ​

  must hold exactly: the trapezoidal block of every phase kernel is , its glassy modulus, and the Volterra products, inverses and layered recurrences all preserve that block. The two sides below come from disjoint code paths — the time-domain ALV pipeline and the elastic Mori–Tanaka estimate — so the agreement validates one against the other.

julia
using Printf
for model in (:layers, :whole_pores), t0 in (2 / 3, 2.0)
    alv = E0 * creep_curve(N, α_solid, t0, [t0, t_max], model; fixed = false)[1]
    ela = elastic_ref(t0, model)
    @printf("%-13s t₀ = %.3f   ALV %9.6f   elastic %9.6f   rel. err. %.1e\n",
        model, t0, alv, ela, abs(alv - ela) / ela)
end
layers        t₀ = 0.667   ALV  1.565784   elastic  1.565784   rel. err. 7.1e-16
layers        t₀ = 2.000   ALV  0.708016   elastic  0.708016   rel. err. 6.3e-16
whole_pores   t₀ = 0.667   ALV  1.940932   elastic  1.940932   rel. err. 4.6e-16
whole_pores   t₀ = 2.000   ALV  0.868795   elastic  0.868795   rel. err. 0.0e+00

A note on the two topologies ​

The :layers composite sphere and the :whole_pores collection of   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.

Reproducing Echoes: the two topologies genuinely differ

MeanFieldHomogenization reproduces Echoes for both topologies to better than 1 % (:layers ≈ 1.60 → 11.06; :whole_pores ≈ 1.96 → 17.54 at  ), and the two differ by in .

That gap is a modeling result, not a discrepancy: the composite-sphere packing of [60] is an efficient model — one Eshelby problem instead of   — not an exact reformulation of the separate-inclusion RVE. Choosing between them is choosing a morphology.