Plasticity and material models

Materials

Barcelona Basic Model

Elastoplasticity for unsaturated soils, in which suction is a second loading variable: drying expands the yield surface, wetting shrinks it, and a soil wetted under constant stress can therefore be pushed to yield and collapse.

PoroMechanics.BBM Type

BBM(; κ, λ0, r, β, M, k_s, nu, e0, p_ref, p_min)

Field	Meaning
κ	elastic swelling index [-]
λ0	virgin compression index at zero suction [-]
r	ratio [-]
β	rate at which the LC curve saturates [Pa⁻¹]
M	slope of the critical state line [-]
k_s	cohesion–suction coefficient [-]
nu	Poisson's ratio [-]
e0	initial void ratio [-]
p_ref	reference pressure of the LC law [Pa]
p_min	floor on the mean compressive stress [Pa]

Defaults are the parameter set of the reference test case.

source

PoroMechanics.BBMState Type

BBMState(σ, ε, pc_star, εv_p, suction)

Internal variables: the stress and strain carried from the previous step (the elastic law is pressure dependent), the preconsolidation pressure at zero suction pc_star — the hardening variable — the accumulated volumetric plastic strain, and the suction the point currently sits at.

source

PoroMechanics.compression_index Function

compression_index(m, s) -> λ(s)

 : the soil compresses less readily as it dries, tending to at high suction.

source

PoroMechanics.preconsolidation Function

preconsolidation(m, s, pc_star) -> p₀(s)

The Loading–Collapse curve, which is the model's central idea:

grows with suction, so drying expands the elastic domain and wetting contracts it. A stress state left untouched while the soil is wetted can therefore end up outside the yield surface — which is collapse.

source

PoroMechanics.yield_function Function

yield_function(m, p, q, s, pc_star) -> f

An ellipse in the plane running from a tensile intercept — the apparent cohesion suction confers — to the preconsolidation pressure . Negative inside, zero on the surface.

source

yield_function(m::DruckerPrager, σ) -> f

   . Negative inside the cone, zero on it.

source

PoroMechanics.mean_pressure Function

Mean compressive net stress    [Pa].

source

PoroMechanics.equivalent_stress Function

Von Mises equivalent stress   [Pa].

source

PoroMechanics.bbm_moduli Function

bbm_moduli(m, p) -> (K, G)

   and    , with p floored at m.p_min because a vanishing modulus would break the solve and a soil in real tension is outside this law.

source

PoroMechanics.suction_stress_increment Function

suction_stress_increment(m, p, s, s_n) -> Δσ_m

The elastic mean-stress change caused by a change of suction at constant strain,

the second half of the BBM's elastic law: volumetric strain responds to suction as well as to stress, ``d\varepsilon_v^e = \frac{\kappa}{1+e_0}\frac{d\bar p}{\bar p}

• \frac{\kappa_s}{1+e_0}\frac{ds}{s + p_{atm}}``.

The atmospheric offset is what keeps the logarithm finite at zero suction, where a saturated soil must still have a defined stiffness. Drying raises s, so at constant net stress the sample shrinks — the reversible half of the shrinkage, the irreversible half being what the LC curve produces.

source

PoroMechanics.log_mean Function

log_mean(a, b)

The logarithmic mean  , with the series    ,  , taking over near   where the quotient is .

It is the exact average of a quantity that varies exponentially between a and b, which is what a pressure obeying    with   does over a step of constant strain rate. That is why it, and not the arithmetic mean, is the right modulus for the deviatoric part of the step.

source

PoroMechanics.step_shear_modulus Function

step_shear_modulus(m, p_n, p_tr, scheme) -> G

The shear modulus to use across a step, from the pressures at its two ends.

Under Val(:exact) it is evaluated at the logarithmic mean pressure, which integrates      exactly when varies exponentially across the step — second order, against the first order of freezing G at the incoming state. Under Val(:explicit) it is frozen at p_n, which is what Bil does.

A constant G_const short-circuits both.

source

PoroMechanics.trial_stress Function

trial_stress(m, ε, s, state, scheme)

The elastic predictor: the stress the material would carry at strain ε and suction s if the step were entirely elastic.

Under Val(:exact) the volumetric part is integrated in closed form. The BBM's elastic law is

which separates and integrates to

with no step-size error at all. Expanding the exponential to first order recovers  , the explicit update, so the two schemes agree in the limit and differ by — which over a path that loads and unloads several times is what accumulates into a per-cent-level error.

Under Val(:explicit) the incremental form is used instead, to reproduce Bil.

source

PoroMechanics.hardening_modulus Function

hardening_modulus(m) -> (1+e₀)/(λ(0) − κ)

The coefficient of the hardening law    .

source

PoroMechanics.dyield_dp Function

   .

source

PoroMechanics.dyield_dq Function

 .

source

PoroMechanics.deviator Function

Deviatoric stress tensor.

source

PoroMechanics.deviatoric_tolerance Function

deviatoric_tolerance(m, p_tr) -> q_tol

Below what deviatoric stress a state counts as hydrostatic.

This is not defensive padding: at   the invariant   is the apex of a cone and is not differentiable, so   is a unit tensor pointing in the direction of e — whatever that direction happens to be. Under isotropic loading e is pure round-off, so the expression yields a full-magnitude tensor with a random orientation, and it poisons the tangent rather than vanishing from it.

The tolerance is relative to the stress the material actually carries, because an absolute one would be meaningless across the range of pressures a soil sees.

source

PoroMechanics.return_residual Function

return_residual(m, x, p_tr, q_tr, s, pc_star_n, K, G) -> (r1, r2)

Residual of the return map, in the two unknowns x = (p, Δγ).

Given and , the deviatoric return and the hardening follow in closed form, so the system reduces to two equations: the volumetric return, and consistency.

The hardening is integrated exponentially,   , which keeps it positive whatever the step size — a plain forward increment can drive it negative and destroy the yield surface.

source

PoroMechanics.solve_return_map Function

solve_return_map(m, p_tr, q_tr, s, pc_star_n, p_n, scheme; tol, maxiter)

Newton on the two-equation residual, with the Jacobian from ForwardDiff. Returns (p, Δγ, q, pc_star, converged).

source

PoroMechanics.elastoplastic_tangent Function

elastoplastic_tangent(m, C_e, p, q, s, pc_star, Δγ, K, G, dev_tr, q_tr, ::Val{dim})

The continuum elastoplastic tangent, from the consistency condition differentiated at the converged point:

with  . It is obtained from  , not by differentiating through the Newton loop of the return map — which is the right instinct, and also the only practical one, since nested automatic differentiation through that loop returns NaN.

What this is not. The algorithmic (consistent) tangent, which is the exact Jacobian of the discrete return map, differs from this one by terms of order . Measured on a plastic isotropic step, the gap against a finite-difference Jacobian falls from 65 % at    to 15 % at   — correct in the limit, and enough to drive a global Newton, but at a linear rather than quadratic rate. algorithmic_tangent is what material_response actually returns; this one is kept because the comparison between the two is the measurement that justifies it.

source

PoroMechanics.algorithmic_tangent Function

algorithmic_tangent(m, C_e, p, Δγ, p_tr, q_tr, s, pc_star_n, K, G, dev_tr, ::Val{dim})

The algorithmic tangent: the exact Jacobian of the discrete return map, as opposed to the continuum tangent of elastoplastic_tangent.

The returned stress is a function of the trial state through the converged unknowns,

and depend on implicitly, through  . The implicit function theorem gives those sensitivities from the same 2×2 Jacobian Newton already assembled,

so the tangent costs two extra back-substitutions, not a second solve. Chaining with   closes it.

This is what makes a global Newton converge quadratically. It is also where Julia earns its place on this model: ForwardDiff supplies and exactly, where Bil differentiates the same residual by hand.

source

PoroMechanics.ContinuumTangent Type

ContinuumTangent(model)

A material that answers exactly like model but returns the continuum tangent instead of the algorithmic one.

It exists to be measured against, not to be used: it is how the claim that the algorithmic tangent is worth its derivation gets tested rather than asserted. Wrapping rather than adding a flag to BBM keeps the choice out of the model, where it does not belong — a constitutive law should not carry a switch that changes only how fast the solver reaches the answer it would have given anyway.

source

PoroMechanics.ExplicitPredictor Type

ExplicitPredictor(model)

A material that answers like model but integrates the elastic law incrementally, freezing the moduli at the incoming state, instead of in closed form.

That is the scheme Bil uses, and it is first-order accurate in the step: over a path that loads and unloads several times the error accumulates to the per-cent level. This wrapper exists so that the comparison against Bil can be made on Bil's own terms, and so that the improvement from integrating exactly is measured rather than asserted. It is not the recommended way to run the model.

source

Drucker-Prager

Perfect plasticity with a non-associated flow rule: friction sets the yield cone, dilatancy sets the plastic flow direction, and keeping the two apart is what lets the model shear without inventing volume.

PoroMechanics.DruckerPrager Type

DruckerPrager(; E, nu, cohesion, friction, dilatancy)

Perfect Drucker-Prager plasticity on isotropic linear elasticity.

field	meaning
elastic	the elastic law, a LinearElastic
cohesion	[Pa]
friction	[rad]
dilatancy	[rad]

Angles are stored in radians. Deck files usually quote degrees — convert at the call site rather than inside the model, so the units of a stored parameter are never in doubt.

Every coefficient is type-parameterized, so a ForwardDiff.Dual can enter the cohesion or the friction angle and a result can be differentiated with respect to them.

source

PoroMechanics.DruckerPragerState Type

DruckerPragerState(σ, ε, εp, γp)

Stress, total strain, plastic strain and the cumulative plastic shear strain  .

γp is carried even though perfect plasticity does not use it: it is the measure a softening or hardening cohesion is a function of, and a state that cannot report how much plastic straining it has been through is useless for anything but the perfect case.

Each tensor carries its own type parameter. That is not pedantry: differentiating with respect to a material parameter — the cohesion, say — makes the stress and the plastic strain Dual while the imposed strain stays Float64, and a state that demanded one shared type would reject exactly the case this package exists to support.

source

PoroMechanics.drucker_prager_return Function

drucker_prager_return(m, ε, εp_n, σ0 = zero(ε)) -> (σ, εp, Δγ, at_apex)

The elastic predictor and its plastic corrector, as a pure function of the strain.

σ0 is a prestress: the stress the material already carries at zero strain. Geotechnical cases are posed this way — base/Poroplast starts at −11.5 MPa isotropic with the mesh undeformed — and folding it into the predictor is the only place it belongs, because the yield check has to see the total stress and not the increment.

Perfect plasticity with a linear elastic predictor makes the smooth return closed-form: consistency on after the return gives

so the returned state follows in one step. No inner Newton is needed, which is the whole difference in cost between this and BBM.

Two returns exist and both are needed. The smooth one slides the trial stress back onto the side of the cone; it is valid only while the corrected stays positive. When the trial state is so tensile that   would go negative, the correct projection is onto the apex, where the cone has no unique normal and the deviator vanishes entirely. Omitting the apex return is the classic way to get a Drucker-Prager implementation that passes every test until a tensile corner appears, and then returns a negative .

source

PoroMechanics.friction_coefficient Function

friction_coefficient(m::DruckerPrager) -> f_f

 , the slope of the cone in the plane.

source

PoroMechanics.cohesion_intercept Function

cohesion_intercept(m::DruckerPrager) -> c_c

  , the value of at which the cone crosses  .

source

PoroMechanics.dilatancy_coefficient Function

dilatancy_coefficient(m::DruckerPrager) -> d_d

 , the volumetric part of the flow direction. Equal to friction_coefficient exactly when the rule is associated.

source

PoroMechanics.apex_pressure Function

apex_pressure(m::DruckerPrager) -> p_apex

The mean stress at which the cone closes, . A stress path that would leave the cone through its tip returns here instead of to its side; see material_response.

Infinite for a frictionless material, which is the von Mises limit and has no apex.

source

A Biot medium with an arbitrary skeleton

PoroMechanics.BiotPlastic Type

BiotPlastic(; skeleton, b, beta, N, k, mu_l)

A Biot medium whose skeleton is an arbitrary elastoplastic material.

BiotPoroelastic hard-wires a LinearElastic skeleton, which is what makes poro_response short and what stops it short of base/Poroplast: that case is a Drucker-Prager skeleton under a Biot coupling, and no shipped date of it is elastic. This struct is the same hinge with the skeleton left open.

field	meaning
skeleton	the drained stress–strain law — any AbstractMaterial
b	Biot coefficient of the stress coupling and of the elastic porosity change [-]
beta	coefficient with which the plastic volume change enters the porosity [-]
N	storage modulus at constant strain [Pa⁻¹]
k	intrinsic permeability [m²]
mu_l	dynamic viscosity [Pa·s]

b and beta are separate because Bil's Poroplast keeps them separate: its porosity is

so a plastic volume change opens pore space at a different rate than an elastic one. They coincide in base/Poroplast, where both are 0.8, and a comparison that assumed they always do would pass there and be wrong everywhere else.

source

PoroMechanics.porosity Function

porosity(m::BiotPlastic, phi0, ε, εp, p, p0) -> φ

Lagrangian porosity, following Bil's Poroplast:

Splitting the volumetric strain into its elastic and plastic parts, each with its own coefficient, is what b and beta are for. With beta == b this collapses to the usual      .

source

porosity(m::FickModel, region) -> φ

Porosity of cell region region. Dispatch rather than a run-time branch, so a single-region model pays nothing for the generality.

source