Nonlinear homogenization: the secant method on a porous plastic solid

Every scheme in MeanFieldHom is linear: it maps phase stiffnesses to an effective stiffness. A nonlinear material can still be treated with those same tools, by replacing each phase with a linear comparison material whose modulus is re-estimated from the strain the phase actually sees. This page walks through one such construction — the modified secant method ([59], [60]) — on the classical test case: a porous solid whose matrix is elastic–perfectly plastic, loaded hydrostatically.

Two ingredients make it work, and both already exist in the package:

1. the composite sphere assemblage (LayeredSphere), which resolves the radial profile of the plastic strain by subdividing the solid into concentric shells, each carrying its own secant modulus;
2. automatic differentiation through a self-consistent scheme, which delivers the second moment of the strain in each shell without ever touching a local field.

The nonlinear closure, and why a derivative gives it

Write the macroscopic strain as ``\boldsymbol{E} = \tfrac{E_v}{3}\boldsymbol{1}

• \boldsymbol{E}d$with$Ev = \mathrm{tr}\,\boldsymbol{E}`` its volumetric

part and $\boldsymbol{E}_d$ its deviator. For an isotropic effective behaviour of moduli $(k^{\hom}, \mu^{\hom})$ the macroscopic elastic energy is

\[W(\boldsymbol{E}) = \tfrac{1}{2}\,\boldsymbol{E}:\mathbb{C}^{\hom}:\boldsymbol{E} = \tfrac{1}{2}\,k^{\hom} E_v^{2} + \mu^{\hom}\,\boldsymbol{E}_d\!:\!\boldsymbol{E}_d .\]

Now let $\mu_i$ be the shear modulus of shell $i$, of volume fraction $f_i$. Differentiating $W$ with respect to $\mu_i$ — at fixed macroscopic strain, so that only the explicit dependence survives — gives the second moment of the deviatoric strain in that shell ([61], [59]):

\[\bigl\langle \boldsymbol{\varepsilon}_d\!:\!\boldsymbol{\varepsilon}_d \bigr\rangle_i \;=\; \frac{1}{f_i}\,\frac{\partial W}{\partial \mu_i} \;=\; \frac{1}{f_i}\left( \tfrac{1}{2}\frac{\partial k^{\hom}}{\partial \mu_i}\,E_v^{2} + \frac{\partial \mu^{\hom}}{\partial \mu_i}\, \boldsymbol{E}_d\!:\!\boldsymbol{E}_d \right).\]

This is the whole point of the method: a quantity that looks local — the quadratic average of the strain inside one shell — is obtained from the derivative of the effective moduli, which ForwardDiff supplies exactly.

Writing $x_i$ for the resulting equivalent deviatoric strain, $x_i = \sqrt{\langle\boldsymbol{\varepsilon}_d\!:\!\boldsymbol{\varepsilon}_d \rangle_i}$, a von Mises solid of elastic shear modulus $\mu_s$ responds with the secant modulus

\[\mu_i = \begin{cases} \mu_s, & x_i \le \varepsilon_0 \quad\text{(still elastic)},\\[2ex] \dfrac{\sqrt{2/3}\;\sigma_0}{2\,x_i}, & x_i > \varepsilon_0 \quad\text{(yielding)} , \end{cases} \qquad \varepsilon_0 = \frac{\sqrt{2/3}\;\sigma_0}{2\mu_s}.\]

The $n$ shells are coupled through $\mathbb{C}^{\hom}$, so this is a fixed point on the vector $(\mu_1,\dots,\mu_n)$ — solved below by direct iteration.

using MeanFieldHom
using TensND
using ForwardDiff
using Printf
using Plots
gr()  # headless backend; GKSwstype is set to "100" before Literate runs

Plots.GRBackend()

§1 The porous material

The solid is elastic–perfectly plastic: bulk modulus $k_s$, shear modulus $\mu_s$, yield stress $\sigma_0$ (a threshold on the deviatoric stress norm $\|\boldsymbol{\sigma}_d\| = 2\mu\|\boldsymbol{\varepsilon}_d\|$). The porosity is $f$. Plasticity only enters through $\sigma_0$; the bulk modulus never changes, since a von Mises solid is plastically incompressible.

const KS = 2.0e3      # bulk modulus of the solid phase
const MUS = 1.0e3     # shear modulus of the solid phase
const SIG0 = 1.0      # yield stress
const FPORE = 0.1     # porosity

0.1

The RVE is a composite sphere: a void core of volume fraction $f$, surrounded by $n$ shells of equal volume fraction $f_i = (1-f)/n$. Radii follow from cumulative volume fractions, $r_i = (\text{cumulated } f)^{1/3}$, the outer radius being 1.

function porous_sphere(μs, f)
    n = length(μs)
    T = eltype(μs)
    tiny = T(1.0e-8) * KS
    # Radii: the core carries f, then each shell adds (1-f)/n.
    radii = ntuple(i -> i == 1 ? f^(1 / 3) : (f + (1 - f) * (i - 1) / n)^(1 / 3), n + 1)
    # Every layer must share one element type, or ForwardDiff cannot promote
    # the tuple when only the shells carry Duals.
    moduli = (
        TensISO{3}(3tiny, 2tiny),
        ntuple(i -> TensISO{3}(3 * T(KS), 2 * μs[i]), n)...,
    )
    return LayeredSphere(radii, moduli)
end

porous_sphere (generic function with 1 method)

The composite sphere fills the whole RVE (volume fraction 1) and is homogenized self-consistently: the reference medium is the effective medium itself, which is what makes the assemblage a model of the porous material rather than of dilute voids in a matrix.

The property declared for the phase is irrelevant here — for a heterogeneous inclusion MeanFieldHom reads the moduli from the inclusion, not from the phase dictionary — but one must be supplied, so the solid stiffness is used.

function kmu_hom(μs, f)
    C_ref = TensISO{3}(3 * eltype(μs)(KS), 2 * μs[1])
    rve = RVE(:M)
    add_matrix!(rve, Ellipsoid(1.0, 1.0, 1.0), Dict(:C => C_ref))
    add_phase!(rve, :I, porous_sphere(μs, f), Dict(:C => C_ref); fraction = 1.0)
    return collect(k_mu(homogenize(rve, SelfConsistent(), :C)))
end

kmu_hom (generic function with 1 method)

§2 The secant closure by automatic differentiation

One ForwardDiff.jacobian call gives the two rows $\partial k^{\hom}/\partial\mu_i$ and $\partial\mu^{\hom}/\partial\mu_i$ at once — the whole nonlinear closure, in one line.

function secant_update(μs, f, Ev, Ed2)
    n = length(μs)
    J = ForwardDiff.jacobian(m -> kmu_hom(m, f), μs)   # row 1: k, row 2: μ
    fi = (1 - f) / n
    ε₀ = sqrt(2/3) * SIG0 / (2 * MUS)
    return map(1:n) do i
        x = sqrt(max(0.0, (0.5 * J[1, i] * Ev^2 + J[2, i] * Ed2) / fi))
        x ≤ ε₀ ? MUS : sqrt(2/3) * SIG0 / (2x)
    end
end

secant_update (generic function with 1 method)

Direct iteration converges in a few tens of steps; converged is returned so that a silently unconverged point can never be plotted as if it were a result.

function solve_secant(μ0, f, Ev, Ed2; tol = 1.0e-9, maxiter = 200)
    μ = copy(μ0)
    for _ in 1:maxiter
        μnew = secant_update(μ, f, Ev, Ed2)
        Δ = maximum(abs.(μnew .- μ) ./ μ)
        μ = μnew
        Δ < tol && return μ, true
    end
    return μ, false
end

solve_secant (generic function with 1 method)

§3 The stress–strain response

Loading is purely hydrostatic: $\boldsymbol{E} = \tfrac{E_v}{3}\boldsymbol{1}$, so $\boldsymbol{E}_d = 0$ and the response is read on the mean stress $\Sigma_m = k^{\hom} E_v$. Note that plasticity develops all the same: the void makes the local strain deviatoric even under a purely volumetric macroscopic strain, and it is that local deviator which yields.

Each strain step starts from the previous converged $\mu$ (continuation), which keeps the iteration count low along the whole curve.

function response(n, f, Evs)
    μ = fill(MUS, n)
    Σm = Float64[]
    allok = true
    for Ev in Evs
        μ, ok = solve_secant(μ, f, Ev, 0.0)
        allok &= ok
        push!(Σm, kmu_hom(μ, f)[1] * Ev)
    end
    return Σm, allok
end

Evs = range(0.0, 4.0e-3; length = 51)[2:end]
shells = (1, 2, 5, 20)

curves = Dict{Int, Vector{Float64}}()
for n in shells
    Σm, ok = response(n, FPORE, Evs)
    ok || @warn "n = $n: fixed point did not converge at every strain step"
    curves[n] = Σm
    @printf("n = %2d shells : Σm plateau = %.4f σ₀\n", n, Σm[end] / SIG0)
end

n =  1 shells : Σm plateau = 1.8974 σ₀
n =  2 shells : Σm plateau = 1.6837 σ₀
n =  5 shells : Σm plateau = 1.5708 σ₀
n = 20 shells : Σm plateau = 1.5378 σ₀

The curves show the two regimes: a common elastic branch of slope $k^{\hom}$, then a plateau once the whole solid has yielded. Refining the radial discretization lowers the plateau and converges.

plt = plot(
    xlabel = "volumetric strain Eᵥ", ylabel = "mean stress Σₘ / σ₀",
    legend = :bottomright, title = "Porous plastic solid, hydrostatic loading (f = $FPORE)",
)
for n in shells
    plot!(plt, Evs, curves[n] ./ SIG0; marker = :circle, markersize = 2,
        label = "n = $n shell" * (n > 1 ? "s" : ""))
end
hline!(plt, [(2 / 3) * log(1 / FPORE)]; linestyle = :dash, color = :black,
    label = "hollow sphere, (2/3)ln(1/f)")
plt

§4 What the plateau converges to

For a rigid–perfectly plastic hollow sphere under hydrostatic loading, limit analysis gives the exact collapse stress $\Sigma_m = \tfrac{2}{3} \sigma_0\ln(1/f)$ — the hydrostatic point of the Gurson criterion ([62]). Refining the radial discretization drives the secant estimate onto it:

@printf("\nexact rigid-plastic hollow sphere : %.4f σ₀\n", (2 / 3) * log(1 / FPORE))
for n in shells
    @printf(
        "  n = %2d : %.4f σ₀   (%+.1f %%)\n", n, curves[n][end] / SIG0,
        100 * (curves[n][end] / SIG0 / ((2 / 3) * log(1 / FPORE)) - 1)
    )
end

exact rigid-plastic hollow sphere : 1.5351 σ₀
  n =  1 : 1.8974 σ₀   (+23.6 %)
  n =  2 : 1.6837 σ₀   (+9.7 %)
  n =  5 : 1.5708 σ₀   (+2.3 %)
  n = 20 : 1.5378 σ₀   (+0.2 %)

So the number of shells is not a numerical detail — it is the model, and the two ends of the sweep are both meaningful.

A single shell reproduces the classical variational estimate. With $n = 1$ the second moment is taken over the whole solid at once, which is exactly the linear-comparison construction of the variational / modified secant method ([60], [59]): one uniform secant modulus for the entire matrix. The resulting plateau matches $\tfrac{2}{3}\sigma_0(1-f)/\sqrt{f}$ to five digits — verified here at $f = 0.05$, $0.1$ and $0.3$:

@printf("\n  n = 1 plateau                    : %.5f σ₀\n", curves[1][end] / SIG0)
@printf("  (2/3)·(1-f)/√f                   : %.5f σ₀\n", (2 / 3) * (1 - FPORE) / sqrt(FPORE))

  n = 1 plateau                    : 1.89737 σ₀
  (2/3)·(1-f)/√f                   : 1.89737 σ₀

That estimate sits about 24 % above the exact collapse load at this porosity: a uniform secant modulus cannot represent a solid that yields near the void long before it yields far from it, and the quadratic average of $\boldsymbol{\varepsilon}_d$ over the whole shell overstates the dissipation.

Resolving the radial profile removes that gap. Each shell then carries its own secant modulus, the plastic front is free to progress outwards from the void, and the estimate converges to the exact limit — within 0.2 % at $n = 20$. Over the range computed here the approach is monotone and always from above.

This page was generated using Literate.jl.