Constitutive laws

Constitutive layer

Material laws owned by the package. Every law carries its coefficients as type parameters, so a result can be differentiated with respect to the parameters themselves and not only with respect to the unknowns — the property that makes inverse calibration and sensitivity analysis possible.

Retention curves

PoroMechanics.AbstractRetention Type

AbstractRetention

Supertype of the water retention curves. A concrete curve implements

• saturation(curve, pc) — liquid saturation at capillary pressure pc

• dsaturation_dpc(curve, pc) — its derivative

source

PoroMechanics.VanGenuchten Type

VanGenuchten(a, n, m)
VanGenuchten(a, m)

Van Genuchten retention curve (van Genuchten, 1980):

a is the air-entry pressure [Pa]; n and m are dimensionless exponents.

The two-argument form applies the classical constraint   . The three-argument form leaves the exponents free — which is not pedantry: published parameter sets often quote a rounded n (1.06383 rather than 1.0638297…), and silently recomputing it from m shifts the curve.

source

PoroMechanics.Gardner Type

Gardner(α)

Exponential retention curve (Gardner, 1958):

α [Pa⁻¹] is the inverse of a characteristic capillary pressure.

The exponential model is the one nonlinearity for which the steady Richards equation integrates in closed form, which is what makes it a verification case rather than a fitting curve. Real soils are usually better described by VanGenuchten.

source

PoroMechanics.ExponentialCutoff Type

ExponentialCutoff(raw, p_c3)

Wrap a retention curve so that below p_c3 it is replaced by an exponential branch joined at p_c3:

For Van Genuchten exponents n > 1, the raw slope vanishes as  . The exponential branch retains a nonzero slope near zero while matching the raw curve's value at the junction. It extends to negative pressures, approaching 1 asymptotically. Clamping it to 1 at zero would introduce a jump in stored water and remove the pressure dependence of storage in nearly saturated cells.

source

PoroMechanics.saturation Function

saturation(curve, pc)

Liquid saturation at capillary pressure pc [Pa]. The saturated branch is law-dependent; ExponentialCutoff approaches 1 asymptotically.

source

PoroMechanics.dsaturation_dpc Function

dsaturation_dpc(curve, pc)

Derivative [Pa⁻¹] of the retention curve.

source

Relative permeability

PoroMechanics.AbstractRelativePermeability Type

AbstractRelativePermeability

Supertype of the relative permeability curves. A concrete curve implements relative_permeability.

source

PoroMechanics.Mualem Type

Mualem(a, m)

Mualem relative permeability built on a Van Genuchten effective saturation (Mualem, 1976), (van Genuchten, 1980):

a and m are the parameters of the curve, which need not equal those of the retention curve.

The derivative   diverges as  . In the dry zone is physically zero, so below S_e = 1e-14 the curve returns exactly zero, which keeps the gradient finite for ForwardDiff.

source

PoroMechanics.GardnerKrl Type

GardnerKrl(α)

Exponential relative permeability (Gardner, 1958):

α [Pa⁻¹] need not equal the α of the matching Gardner retention curve.

Substituting this law into the steady Richards equation turns it into a linear ordinary differential equation for , so the steady profile above a water table has a closed form — see the Gardner infiltration benchmark.

source

PoroMechanics.PowerLawKrl Type

PowerLawKrl(a, n, m)

Relative permeability of the same functional form as a Van Genuchten curve, but with independent exponents:

Used by the non-isothermal drying model, where the two materials are fitted with n = 2 and different m, so the Van Genuchten constraint    does not hold.

source

PoroMechanics.relative_permeability Function

relative_permeability(curve, pc)

Relative permeability to the liquid,  , at capillary pressure pc [Pa]. Returns 1 for pc ≤ 0.

source

PoroMechanics.gas_relative_permeability Function

gas_relative_permeability(sl)

Relative permeability to the gas phase as a function of liquid saturation:

Zero at full saturation. sl is clamped just below 1 so the expression stays differentiable at the endpoint.

source

Tortuosity

PoroMechanics.AbstractTortuosity Type

AbstractTortuosity

Supertype of the tortuosity models. A concrete model implements tortuosity(model, φ, S_l).

source

PoroMechanics.OhJang Type

OhJang(; phi_c, n, ds, tau_agg, sat_exponent = 4.5, cap_fraction = 0.5)

Tortuosity of a cementitious material after Oh & Jang (Oh and Jang, 2004), extended by a saturation factor:

with the paste contribution built from the capillary porosity  :

Field	Meaning
phi_c	percolation threshold [-]
n	Oh-Jang exponent [-]
ds	solid-to-pore diffusivity ratio [-]
tau_agg	aggregate tortuosity factor [-], calibrated on the measured
sat_exponent	exponent of the saturation factor, 4.5 by default
cap_fraction	fraction of the total porosity that is capillary, 0.5 by default

For a saturated medium   and the saturation factor is one, which is why the saturated and unsaturated chloride-ingress models share this single implementation.

source

PoroMechanics.tortuosity Function

tortuosity(model, φ, S_l = 1)

Tortuosity factor at porosity φ and liquid saturation S_l. For a saturated medium the saturation argument may be omitted.

source

Poroelasticity

PoroMechanics.AbstractPoroelastic Type

AbstractPoroelastic

Supertype of the poroelastic materials.

source

PoroMechanics.BiotPoroelastic Type

BiotPoroelastic(; E, nu, k, mu_l, b, N)

Linear isotropic Biot poroelasticity.

Field	Meaning
E, nu	drained Young's modulus [Pa] and Poisson's ratio [-]
k	intrinsic permeability [m²]
mu_l	dynamic viscosity of the fluid [Pa·s]
b	Biot coefficient [-]
N	storage modulus at constant strain [Pa⁻¹], the inverse of the Biot modulus M

Permeability sits in the material rather than beside it because Biot's theory couples Darcy flow to the skeleton: consolidation_coefficient is a property of the pair, not of either half.

source

PoroMechanics.lame Function

lame(m) -> (λ, μ)

Lamé coefficients [Pa], drained.

source

PoroMechanics.shear_modulus Function

shear_modulus(m) -> G [Pa]

source

PoroMechanics.bulk_modulus Function

bulk_modulus(m) -> K = λ + 2μ/3 [Pa]

Drained bulk modulus.

source

PoroMechanics.oedometric_modulus Function

oedometric_modulus(m) -> M_o = λ + 2μ = K + 4μ/3 [Pa]

The modulus that governs uniaxial strain, and the one that appears throughout the consolidation constants.

source

PoroMechanics.biot_modulus Function

biot_modulus(m) -> M = 1/N [Pa]

source

PoroMechanics.compaction_coefficient Function

compaction_coefficient(m) -> c_m [Pa⁻¹]

Geertsma's uniaxial compaction coefficient,

the constant of proportionality in    for problems whose displacement field is irrotational.

source

PoroMechanics.storage_coefficient Function

storage_coefficient(m) -> S [Pa⁻¹]

Storage at constant stress,      .

source

PoroMechanics.consolidation_coefficient Function

consolidation_coefficient(m) -> c [m²/s]

Uniaxial diffusivity  , obtained by eliminating the strain between the constitutive law and the storage equation under uniaxial conditions.

source

PoroMechanics.hydraulic_conductivity Function

hydraulic_conductivity(m) -> κ = k/μ_l [m²/(Pa·s)]

source

PoroMechanics.skempton Function

skempton(m) -> B [-]

Skempton's coefficient     : the fraction of an isotropic stress increment carried by the pore fluid under undrained conditions.

source

PoroMechanics.undrained_poisson Function

undrained_poisson(m) -> ν_u [-]

which tends to 1/2 as the fluid and grains become incompressible ( ).

source

PoroMechanics.undrained_bulk_modulus Function

undrained_bulk_modulus(m) -> K_u = K + b²M [Pa]

source

The stress–strain interface

The mechanical layer deliberately knows nothing about pore pressure: its signature is the one the Ferrite ecosystem already uses, so a constitutive model written against it travels beyond this package. The poromechanical coupling sits on top.

PoroMechanics.AbstractMaterial Type

AbstractMaterial

Supertype of the stress–strain models. A concrete material implements material_response and initial_state.

source

PoroMechanics.AbstractMaterialState Type

AbstractMaterialState

Supertype of the internal-variable containers. States are immutable: a material returns a new state rather than mutating the one it was given, which keeps the response a pure function of its arguments and therefore differentiable.

source

PoroMechanics.NoState Type

NoState()

The state of a material that has none — linear elasticity, for instance.

source

PoroMechanics.initial_state Function

initial_state(mat) -> state

The internal variables a material starts from. Called once per quadrature point.

source

initial_state(m::LogarithmicElastic, σ0) -> LogarithmicElasticState

Start the material from a known in-situ stress. A pressure-dependent modulus is undefined without one, which is why this material cannot use the zero-argument form.

source

initial_state(m::BBM, σ0, pc_star; suction = 0)

Start from an in-situ stress and a preconsolidation pressure. Both are required: the elastic stiffness is undefined without a stress, and the yield surface without pc_star.

source

PoroMechanics.material_response Function

material_response(mat, ε, state, Δt) -> (σ, ∂σ∂ε, state_new)

Stress, consistent tangent and updated internal variables for the strain increment ending at ε over a step Δt.

• ε is a SymmetricTensor{2,dim}, σ likewise, and ∂σ∂ε a SymmetricTensor{4,dim}.

• state is whatever initial_state returned, and state_new replaces it.

For a rate-independent material Δt is unused but kept in the signature: it is what a viscoplastic or a creeping material needs, and changing the signature later would break every model written against it.

source

material_response(m::BBM, ε, state, Δt) -> (σ, ∂σ∂ε, state_new)

Tension-positive strain in, tension-positive stress out. The suction is read from the state, so a wetting or drying step is applied by handing in a state whose suction has changed — which is exactly how a collapse test is driven.

The tangent returned on a plastic step is the algorithmic one — the exact Jacobian of the discrete return map, from algorithmic_tangent. The continuum tangent is kept as elastoplastic_tangent for comparison; the two differ by .

source

material_response(m::BBM, ε, s, state, Δt) -> (σ, ∂σ∂ε, state_new)

The five-argument form, in which the suction is a loading variable alongside the strain.

That is what it is: the BBM has two controls, and a wetting or drying step is imposed, not undergone. Only the increment has meaning for the elastic response, so the state carries the suction it was computed at and this form supplies the new one. The four-argument form reads the suction back out of the state, which leaves the increment at zero — right for a mechanical step, and wrong for a hydraulic one.

source

material_response(m::DruckerPrager, ε, state, Δt) -> (σ, ∂σ∂ε, state)

The stress, the consistent tangent and the updated state.

The tangent is not written out by hand: it is Tensors.gradient of drucker_prager_return, so it is the derivative of the return actually performed rather than of the one that was meant to be. That is the same rule the finite volume models follow — the algorithmic tangent of a non-associated cone is exactly the kind of expression whose sign errors survive every test that only checks the stress.

One branch is treated differently, on purpose. At the apex the returned stress is a constant, so the true algorithmic tangent is zero and a global stiffness assembled from it is singular. The elastic stiffness stands in there. It costs Newton its quadratic rate on steps that hit the tip; saying so is better than shipping a matrix that cannot be factored.

source

PoroMechanics.stress_controlled_response Function

stress_controlled_response(mat, σ_target, [s,] state, Δt; tol, maxiter) -> (ε, σ, state_new, iters)

The strain that produces a prescribed stress, by Newton on the material response.

Laboratory paths are prescribed in stress, not in strain: an isotropic compression test holds a cell pressure and measures the volume change. A model whose interface is   cannot follow one directly, so this inverts it, using the same consistent tangent the global solve uses — which is the cheapest possible check that the tangent is right, since a wrong one shows up immediately as a failure to converge.

The optional s is the second loading variable of a model that has one, the suction of the Barcelona Basic Model in particular.

tol is relative to the target stress, so it means the same thing at 1 kPa and at 1 MPa.

source

PoroMechanics.LinearElastic Type

LinearElastic(λ, μ)
LinearElastic(; E, nu)

Isotropic linear elasticity, ``\sigma = \lambda,\mathrm{tr}(\varepsilon)\mathbf

• 2\mu,\varepsilon``.

Type-parameterized like every other law here, so a response can be differentiated with respect to λ or μ and not only with respect to the strain.

source

PoroMechanics.elastic_stiffness Function

elastic_stiffness(mat, ::Val{dim}) -> SymmetricTensor{4,dim}

The isotropic stiffness ``C_{ijkl} = \lambda\delta_{ij}\delta_

• \mu(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})``.

source

PoroMechanics.skeleton Function

skeleton(m::BiotPoroelastic) -> LinearElastic

The drained skeleton of a Biot material, as a stress–strain model in its own right. This is the hinge between the two layers: everything below it is ordinary solid mechanics, and everything above it is poromechanics.

source

PoroMechanics.total_stress Function

total_stress(b, σ_eff, p) -> σ

Total Cauchy stress of a poroelastic medium,     , from the effective (skeleton) stress and the pore pressure. b is the Biot coefficient.

Sign convention: tension positive, so a positive pore pressure relieves the skeleton.

source

PoroMechanics.poro_response Function

poro_response(m::BiotPoroelastic, ε, p, state, Δt) -> (σ, ∂σ∂ε, ∂σ∂p, state_new)

Total stress of a Biot medium and its tangents, obtained by asking the skeleton for its mechanical response and adding the pressure coupling.

∂σ∂p = -b I is constant for linear Biot, but it is returned rather than assumed so that a model with a stress-dependent Biot coefficient can be dropped in without changing the callers.

source

poro_response(m::BiotPlastic, ε, p, state, Δt) -> (σ, ∂σ∂ε, ∂σ∂p, state)

Total stress of the medium and its tangents, from the skeleton's own response.

Identical in shape to the BiotPoroelastic method — ask the skeleton, then add the pressure coupling — which is the point: the poromechanical layer does not care whether the skeleton yielded, and a plastic skeleton needs no new coupling code.

source

Pressure-dependent elasticity

PoroMechanics.LogarithmicElastic Type

LogarithmicElastic(; κ, nu, e0, p_min)

Elasticity with a bulk modulus proportional to the mean effective compressive stress,

the elastic law of critical-state soil mechanics, and the elastic part of the Barcelona Basic Model. κ is the swelling index, e0 the initial void ratio.

Tension positive throughout, so   in compression. Under tension the modulus would go to zero or negative, which is unphysical and would break the solve, so p_min floors it: a soil that is genuinely in tension is outside this law's range and the floor says so rather than producing a plausible number.

Unlike LinearElastic, this material has a state and its tangent changes with it — which is what makes it the honest rehearsal for a return mapping.

source

PoroMechanics.LogarithmicElasticState Type

LogarithmicElasticState(σ, ε)

Internal variables of LogarithmicElastic: the stress and strain the material carries from the previous step. A hypoelastic law needs them because its stiffness depends on where it currently sits, so the response is an increment rather than a function of the total strain.

source

PoroMechanics.tangent_moduli Function

tangent_moduli(m, σ) -> (K, G)

Bulk and shear moduli at the current stress.

source

PoroMechanics.mean_compressive_stress Function

mean_compressive_stress(m, σ) -> p̄ [Pa]

  , floored at m.p_min.

source

Unsaturated effective stress

When two fluid phases share the pore space there is no single pore pressure, and the question of what replaces it is a modeling choice rather than a formula — so the Bishop coefficient is a model of its own.

PoroMechanics.AbstractBishop Type

AbstractBishop

Supertype of the Bishop coefficient models. A concrete model implements bishop_coefficient(model, pc).

source

PoroMechanics.SaturationBishop Type

SaturationBishop(retention)

 , so the equivalent pore pressure is the saturation-weighted average of the two phase pressures. Coussy's form with the interfacial energy dropped, and the usual default.

source

PoroMechanics.PowerBishop Type

PowerBishop(retention, n)

 , the fitted generalisation used when   overestimates the contribution of the liquid — as it does in fine-grained soils, where part of the water is held in menisci that transmit little stress. n = 1 recovers SaturationBishop.

source

PoroMechanics.bishop_coefficient Function

bishop_coefficient(model, pc) -> χ ∈ [0, 1]

Weight given to the liquid pressure in the equivalent pore pressure at capillary pressure pc [Pa].

source

PoroMechanics.equivalent_pore_pressure Function

equivalent_pore_pressure(model, p_l, p_g = zero(p_l)) -> π

The single pressure that plays the role of in the saturated theory,

At full saturation ( , so  ) this returns exactly, which is the check that the unsaturated theory contains the saturated one.

source

PoroMechanics.unsaturated_total_stress Function

unsaturated_total_stress(b, χmodel, σ_eff, p_l, p_g = 0) -> σ

Total stress of an unsaturated poroelastic medium,     , with the equivalent pore pressure.

Tension positive, as everywhere else here.

source

PoroMechanics.suction_stress Function

suction_stress(χmodel, p_l, p_g = 0) -> χ p_c

The isotropic tension that suction exerts on the skeleton,  : what pulls the grains together and gives an unsaturated soil its apparent cohesion. It is the part of the effective stress that vanishes on wetting, and therefore the mechanism behind collapse on saturation.

source

Tabulated material curves

A retention or relative-permeability curve given as data rather than as a formula — measured points, or a curve digitized from a reference. Interpolation is linear in the tabulated variable and the coefficients are parameterized by their own type, so a table can carry ForwardDiff.Dual values like any closed-form law.

PoroMechanics.Tabulated Type

Tabulated(pc, sl)

A retention curve given as a table rather than a formula: sl[k] is the saturation at capillary pressure pc[k], with linear interpolation between the points and clamping to the end values outside the range.

Measured curves arrive this way, and so do the curves another code has already discretized — Bil writes the table it interpolated next to its deck, and reading that table back is what lets a comparison separate a disagreement about the physics from a disagreement about how a curve was sampled.

pc must be sorted and strictly increasing; the constructor checks it, because a table that is merely nearly sorted produces a plausible curve with a wrong branch in it.

Both vectors are type parameters, so a table of Duals differentiates like any other coefficient — the values of a measured curve are material parameters too.

pc, sl, krl = eachcol(readdlm("billes"))
model = RichardsModel(; retention = Tabulated(pc, sl), rel_perm = TabulatedKrl(pc, krl))

source

PoroMechanics.TabulatedKrl Type

TabulatedKrl(pc, krl)

Relative permeability given as a table, krl[k] at capillary pressure pc[k], linearly interpolated and clamped outside the range. The counterpart of Tabulated for retention; see there for why a table is worth having as a first-class law.

source

PoroMechanics.interpolate_table Function

interpolate_table(x, y, q)

Linear interpolation of the table (x, y) at q, clamped to the end values outside the range.

The + zero(q) on the clamped branches is not decoration: without it a Dual argument falls back to a Float64 return and the derivative is silently lost — the same trap as returning the literal 1.0 from a saturated branch.

source