A damage law written from scratch

The whole page as one script

Every block of this page, in order, is scripts/material_laws/02_damage.jl. From a clone of the repository, julia --project=docs -e 'using Pkg; Pkg.instantiate()' once, then julia scripts/material_laws/02_damage.jl, runs it end to end.

This example writes a law that is not shipped with the package, from its equations to a cyclic test, in about fifteen lines. It is a scalar isotropic damage law of the simplest kind, chosen for being easy to check by hand: the point is the procedure, which is the same for any law with internal variables.

The model

The material is linear elastic, of stiffness , degraded by a scalar damage  :

The damage is driven by the elastic energy of the undamaged material,     . Below a threshold nothing happens; beyond it the damage takes the value   , for which the stored energy    stays equal to . Damage never heals: when the strain decreases, keeps its largest value so far.

Over a step this reads: if   the step is elastic,  ; otherwise solves   .

The law in code

using MeanFieldHomogenization, StaticArrays, Plots

k, μ = 20e3, 12e3                                     # bulk and shear moduli, MPa
J = identity_mandel() * identity_mandel()' / 3        # the projector on volume changes
ℂ = 3k * J + 2μ * (one(SMatrix{6, 6, Float64}) - J)   # isotropic stiffness, Mandel form

Y(ε, p) = transpose(ε) * (p.ℂ * ε) / 2
d_star(ε, p) = 1 - sqrt(p.Y₀ / Y(ε, p))

damage = ResidualLaw(;
    name = :ScalarDamage,
    stress = (ε, y, p) -> (1 - y.d) * (p.ℂ * ε),
    residual = (y, g, prev, Δt, p) -> (d = y.d - d_star(g.ε, p),),
    elastic = (g, prev, Δt, p) -> d_star(g.ε, p) <= prev.y.d,
    layout = (d = 1,),
    p = (ℂ = ℂ, Y₀ = 1e-3),                            # Y₀ in MPa (energy per volume)
)

ScalarDamage (residual route; ε → σ; internal variables d)

Each line is one piece of the model:

code	model
stress
residual	on a damaging step
elastic	: below the threshold, or unloading
layout	one internal variable, the scalar d
p	the stiffness and the threshold

The derivative of the stress with respect to the strain, the consistent tangent of a damaging step, is not written: it comes from the implicit-function theorem applied to the residual. The contract checker confirms it against a finite difference:

check_material_interface(damage; verbose = false)

true

A test with unloading

Uniaxial strain (  ), loaded, unloaded to zero, then loaded further:

ε_axial(t) = t < 0.4 ? 5e-4 * t / 0.4 : t < 0.6 ? 5e-4 * (0.6 - t) / 0.2 : 1e-3 * (t - 0.6) / 0.4
path = LoadingPath(range(0, 1; length = 101); ε11 = ε_axial, ε22 = 0, ε33 = 0)
r = point_response(damage, path)
plot(r[:ε11] * 1e3, r[:σ11]; lw = 2, label = "ScalarDamage", legend = :topleft,
    xlabel = "axial strain ε₁₁ (‰)", ylabel = "axial stress σ₁₁ (MPa)")

The unloading goes back to the origin along a straight line of reduced slope   , and the reloading follows it up to the previous maximum before damage resumes.

Checking it by hand

In uniaxial strain    . The damage is the running maximum of , and the stress follows:

C11 = ℂ[1, 1]
d_ref = accumulate(max, [e == 0 ? 0.0 : max(0.0, 1 - sqrt(1e-3 / (C11 * e^2 / 2))) for e in r[:ε11]])
maximum(abs, r[:d] - d_ref), maximum(abs, r[:σ11] - (1 .- d_ref) .* C11 .* r[:ε11])

(1.1102230246251565e-16, 1.7763568394002505e-15)

Both differences are at the rounding of the numbers.

---

This page was generated using Literate.jl.