Nonlinear solid mechanics

Overview

The solid simulations so far assumed a linear spring: stress proportional to strain, displacement proportional to load. Two physical effects break that assumption, and each is important enough to need its own model. This chapter covers both.

Hyperelasticity is about large deformation. Rubber and soft tissue stretch so far that the small-strain approximation fails, yet they spring back to shape. They stay elastic — no permanent deformation — but the stress–strain law is nonlinear.

Plasticity is about permanent deformation. A metal is elastic up to a stress threshold, then yields: beyond it the material flows and does not return. The threshold is the yield stress, and the permanent motion is plastic strain.

Finite-strain hyperelasticity

Linear elasticity measures strain by the symmetric gradient of $u$ and assumes it is small. That is fine for steel, but wrong for rubber stretched to twice its length. The finite-strain description replaces the small-strain measure with the deformation gradient $F = I + \nabla u$, which captures large rotations and stretches exactly.

A hyperelastic material then goes one step further than a stress law: it stores a scalar strain-energy density $W(F)$ and derives the stress from it. The energy is a function of the invariants of $F$ — quantities that do not change when the block is rigidly rotated, so the material behaves the same regardless of orientation. The default model is the compressible Neo-Hookean solid, whose energy is

\[W = \frac{\mu}{2}(I_1 - 3) - \mu\ln J + \frac{\lambda}{2}(\ln J)^2 ,\]

where $\mu$ and $\lambda$ are the Lamé parameters, $J = \det F$ is the volume ratio ($J = 1$ means no volume change), and $I_1 = \mathrm{tr}(F^{T}F)$. The two terms give the material its two behaviours: the $I_1$ term resists shearing and stretching, and the $\ln J$ terms resist volume change.

HyperElasticity

W = \frac{\mu}{2}(I_1 - 3) - \mu\ln J + \frac{\lambda}{2}(\ln J)^2,

julia> HyperElasticity() isa AbstractPhysics
true

A worked example: stretching a rubber block

Clamp a soft block at $x = 0$ and stretch the $x = 1$ face by $2\%$. The prescribed displacement is given through the dof indices fixed and the values g (a stretch is a Dirichlet condition, not a force):

m = tetrahedralize_box(1.0, 1.0, 1.0, 6, 6, 6)
n = size(m.nodes, 2)

rubber = Material(youngs_modulus=1e6, poisson_ratio=0.49)

fixed = Int[]; g = Float64[]
for i in 1:n
    x = m.nodes[1, i]
    if x < 1e-9
        append!(fixed, 3i-2, 3i-1, 3i); append!(g, 0.0, 0.0, 0.0)
    elseif x > 1 - 1e-9
        append!(fixed, 3i-2, 3i-1, 3i); append!(g, 0.02, 0.0, 0.0)
    end
end

u = solve(HyperElasticity(), m, rubber; fixed=fixed, g=g)
maximum(u[1:3:end])

0.02 m

The $x$ displacement ramps linearly from $0$ at the clamp to the prescribed $0.02$ at the free end — a uniform stretch, as a block pulled from one face should have:

vtk(m; field=u[1:3:end], fieldname="u_x", edges=false)

J2 plasticity

A metal is not springlike forever. Pull it hard enough and it yields: at the yield stress $\sigma_y$ it stops storing energy elastically and starts to flow, accumulating permanent (plastic) strain that remains when the load is removed.

J2 plasticity (von Mises) decides when yielding happens by a single scalar — the von Mises (or deviatoric) stress $q = \sqrt{3 J_2}$, where $J_2$ is the second invariant of the deviatoric stress (the part of the stress that shears, with the hydrostatic pressure subtracted). The rule is: elastic while $q < \sigma_y$, yielding when $q = \sigma_y$. With linear hardening the yield stress grows with the accumulated plastic strain, ``\sigmay \leftarrow \sigmay

• H\,\varepsilon_p``, so the material stiffens as it hardens.

The numerical core is the radial return map. Each Newton step first computes an elastic predictor stress; if it lies outside the yield surface, the stress is projected straight back onto the surface — the shortest (radial) return — which is exactly the J2 flow rule. The consistent tangent from that projection then drives Newton's method.

Plasticity

julia> Plasticity() isa AbstractPhysics
true

A worked example: a bar below yield

A load below the yield stress produces a purely elastic response, identical to the linear elasticity solve. Here a steel bar is pulled with $10^6$ N over its $1\ \text{m}^2$ end face — $1$ MPa, far below the $250$ MPa yield point:

steel = Material(youngs_modulus=200e9, poisson_ratio=0.3,
                 yield_strength=250e6, hardening_modulus=2e9)

right = [i for i in 1:n if m.nodes[1, i] > 1 - 1e-9]
left  = Int[]
for i in 1:n
    m.nodes[1, i] < 1e-9 && append!(left, 3i-2, 3i-1, 3i)
end

loads = zeros(3n)
for i in right
    loads[3i-2] = 1e6 / length(right)
end

up = solve(Plasticity(), m, steel; fixed=left, loads=loads)
maximum(abs.(up))

7.529832156726544e-6 m

The displacement is tiny — $\approx 7.5 \times 10^{-6}$ m — and exactly the linear-elastic prediction, because no element reached yield. The return map only changes the answer once $q$ reaches $\sigma_y$; below that, J2 plasticity is linear elasticity.

vtk(m; field=up[1:3:end], fieldname="u_x", edges=false)

Next steps

Nonlinear solids share their solver with the rest of the nonlinear physics. The Newton iteration and its variants live in Nonlinear solvers, and the flow counterpart of this chapter — the nonlinear Navier–Stokes equations — is in the fluid-flow chapter.