Skip to content

API — MaterialLaws ​

MeanFieldHomogenization.MaterialLaws Module
julia
MeanFieldHomogenization.MaterialLaws

Constitutive laws at a material point: writing a law, linear or not, in a few lines, integrating it over a step with its consistent tangent, driving it along a loading path, and plugging it into any caller of the Gauss-point contract, a finite-element code included.

Four independent pieces

piecewhat it sayshow
the lawwhat the material isResidualLaw (local equations) or PotentialLaw (free energy and dissipation), parameters in p
the local algorithmhow a step is integratedLocalNewton, with the tangent by the implicit-function theorem
the loadingalong which pathLoadingPath, each component driven in strain or in stress
the couplingwhere the law is evaluatedpoint_response at one point, material_response in an element loop

Any law runs with any algorithm, along any path, in any caller. A law is a AbstractMFHMaterial like a homogenized material, so the coupling code of the package serves both.

julia
using MeanFieldHomogenization

law = J2Plasticity(; E = 200e3, ν = 0.3, σ₀ = 250.0, H = 2e3)
path = LoadingPath(range(0, 1; length = 101); ε11 = t -> 5e-3 * t)   # uniaxial tension
r = point_response(law, path)
r[:σ11], r[:p]                        # stress and cumulated plastic strain

Numbers

Strains and stresses are Mandel vectors (SVector{6}), tangents SMatrix{6, 6}: a step allocates nothing. Every number may be a ForwardDiff.Dual, nested or not, in the strain, the state or the parameters: the response then carries its derivatives, exact (no finite differences). The small-strain kinematics (SmallStrain) is the only one so far.

Writing a law ​

MeanFieldHomogenization.MaterialLaws.AbstractPointLaw Type
julia
AbstractPointLaw <: AbstractMFHMaterial

Supertype of the constitutive laws of MaterialLaws: a law written with ResidualLaw or PotentialLaw, integrated over a step by point_update.

Being an AbstractMFHMaterial, a law answers material_response and plugs into any caller of that contract, a finite-element element loop included, exactly as a homogenized material does.

MeanFieldHomogenization.MaterialLaws.ResidualLaw Type
julia
ResidualLaw(; fluxes, residual, layout, gradients = (ε = 6,), flux_names = (σ = 6,),
              p = (;), elastic = nothing, y₀ = zeros, name = :ResidualLaw)
ResidualLaw(; stress, residual, layout, …)

A law given by the equations of its implicit integration over a step.

  • fluxes(g, y, p) returns the fluxes, a NamedTuple by the flux names: (σ = …,) for a mechanical law, (σ = …, φ = …) for a poroelastic one. A mechanical law may give stress(ε, y, p) instead.

  • residual(y, g, prev, Δt, p) returns the residual of the local equations, a NamedTuple by the names of the internal variables (or a flat SVector), whose root is the state at the end of the step, in any units.

  • layout names the internal variables and their sizes, e.g. (εᵖ = 6, p = 1), or (d = 1,) for a scalar damage.

  • gradients, flux_names: the primal fields and their conjugate fluxes, with their sizes, in the same order. The default is the strain and the stress.

  • p holds every parameter the functions read.

  • elastic(g, prev, Δt, p), optional, returns true when the step leaves the internal variables unchanged: the elastic predictor of a rate-independent law, or the threshold of a damage law. Without it every step is solved.

  • y₀, optional, the initial internal variables (zero by default).

The functions receive g, y and the previous state prev (fields prev.y, prev.g) as NamedTuples: g.ε, g.p, y.εᵖ, prev.y.d. A size-6 entry is a symmetric tensor in Mandel components, a size-1 entry a scalar.

point_update solves the residual by Newton's method with a Jacobian by automatic differentiation, and returns the consistent tangent of the fluxes with respect to the gradients by the implicit-function theorem,

every coupling block included, so that no derivative is written by hand.

julia
# A viscous shear flow in series with a spring: σ = 2μ (ε − εᵛ), ε̇ᵛ = dev σ / (2η)
law = ResidualLaw(;
    stress = (ε, y, p) -> 2p.μ * (ε - y.εᵛ),
    residual = (y, g, prev, Δt, p) -> (εᵛ = y.εᵛ - prev.y.εᵛ - Δt * deviator(2p.μ * (g.ε - y.εᵛ)) / (2p.η),),
    layout = (εᵛ = 6,), p = (μ = 1.0, η = 10.0),
)
MeanFieldHomogenization.MaterialLaws.PotentialLaw Type
julia
PotentialLaw(; energy, dissipation = nothing, layout = (;), gradients = (ε = 6,),
               flux_names = (σ = 6,), signs = nothing, p = (;), y₀ = zeros,
               name = :PotentialLaw)

A law given by two potentials, in the framework of generalized standard materials: the free energy of the gradients and the internal variables, and the dissipation potential .

  • energy(g, α, p) returns ;

  • dissipation(α̇, α, p) returns , α being the internal variables at the beginning of the step; it is required as soon as the law has internal variables;

  • each flux is the derivative of with respect to its gradient, times the sign signs gives it (+1 by default): a poroelastic potential written in strain and pore pressure gives   and   , hence signs = (σ = 1, φ = -1);

  • layout, gradients, flux_names, p, y₀ and name as for ResidualLaw.

Over a step of length the internal variables minimize the incremental energy

the stationarity condition being solved as the residual of a ResidualLaw: the gradient and the Hessian come from automatic differentiation, so a law is written as its two potentials and nothing else. Newton's method finds a stationary point; with convex potentials it is the minimum. A step of zero length leaves the internal variables unchanged.

julia
# Nonlinear elasticity, no internal variable
law = PotentialLaw(; energy = (g, α, p) -> p.k / 2 * trace_mandel(g.ε)^2 + p.μ * sum(abs2, deviator(g.ε))^2,
                   p = (k = 10.0, μ = 3.0))
MeanFieldHomogenization.MaterialLaws.PointState Type
julia
PointState(y, g) <: AbstractMaterialState

The state a material point carries from one step to the next: its internal variables y and its gradients g (the strain, and any other primal field of the law) at the end of the step, both flat SVectors laid out as the law's state_layout and gradient_layout say.

The functions of a law read the state of the previous step by name, as prev.y.εᵖ or prev.g.ε. A state is never mutated: point_update returns a new one, so that a step a caller rejects is retried from the old one.

MeanFieldHomogenization.MaterialLaws.kinematics Function
julia
kinematics(law) -> SmallStrain()

The kinematics a law is written in. SmallStrain is the only one so far.

MeanFieldHomogenization.MaterialLaws.SmallStrain Type
julia
SmallStrain()

The kinematics of a law written in infinitesimal strain: the law receives the symmetric strain tensor and returns the Cauchy stress . It is the only kinematics so far; the trait kinematics exists so that a finite-strain one can be added without changing the laws written for this one.

MeanFieldHomogenization.MaterialLaws.state_layout Function
julia
state_layout(law) -> NamedTuple

The names of a law's internal variables and the size of each, e.g. (εᵖ = 6, p = 1) for von Mises plasticity: what a post-processing reads, and what PointResponse indexes by name.

MeanFieldHomogenization.MaterialLaws.gradient_layout Function
julia
gradient_layout(law) -> NamedTuple

The primal fields a law consumes and the fluxes it returns, with their sizes: (ε = 6,) and (σ = 6,) for a mechanical law, (ε = 6, p = 1) and (σ = 6, φ = 1) for a poroelastic one. A size-6 entry is a symmetric tensor in Mandel components, a size-1 entry a scalar.

MeanFieldHomogenization.MaterialLaws.flux_layout Function
julia
flux_layout(law) -> NamedTuple

The fluxes a law returns, with their sizes, in the order of the gradients they are conjugate to (see gradient_layout).

MeanFieldHomogenization.MaterialLaws.law_parameters Function
julia
law_parameters(law) -> NamedTuple

The parameters p the law's functions receive. Building the law again with ForwardDiff.Dual values in their place differentiates its response with respect to them.

Kelvin–Mandel algebra ​

MeanFieldHomogenization.MaterialLaws.Mandel Type
julia
Mandel{T}

A symmetric second-order tensor of dimension 3 in Kelvin–Mandel components, an SVector{6, T} ordered

The basis is orthonormal:     (dot(σ, ε)), the norm of the tensor is norm(v), and a fourth-order tensor with the minor symmetries is an SMatrix{6, 6} acting by the matrix product. The laws of MaterialLaws are written on these vectors.

MeanFieldHomogenization.MaterialLaws.identity_mandel Function
julia
identity_mandel(T = Float64) -> Mandel{T}

The second-order identity , in Kelvin–Mandel components.

MeanFieldHomogenization.MaterialLaws.trace_mandel Function
julia
trace_mandel(v::Mandel)

The trace   .

MeanFieldHomogenization.MaterialLaws.deviator Function
julia
deviator(v::Mandel) -> Mandel

The deviatoric part   .

MeanFieldHomogenization.MaterialLaws.von_mises Function
julia
von_mises(σ::Mandel)

The von Mises equivalent stress  , being the deviator of : equal to in uniaxial tension.

Not differentiable at zero

The square root has no derivative at a vanishing deviator. A law whose response must be differentiated there writes its flow rule with instead, as NortonViscoplasticity does.

MeanFieldHomogenization.MaterialLaws.equivalent_strain Function
julia
equivalent_strain(ε::Mandel)

The equivalent strain  , being the deviator of : work-conjugate to von_mises, and equal to for an incompressible uniaxial strain.

Integrating a step ​

MeanFieldHomogenization.MaterialLaws.AbstractLocalAlgorithm Type
julia
AbstractLocalAlgorithm

Supertype of the algorithms that integrate a law over a step. LocalNewton is the default; a new algorithm subtypes this and adds a method of MaterialLaws._solve_local.

MeanFieldHomogenization.MaterialLaws.LocalNewton Type
julia
LocalNewton(; abstol = 1e-12, reltol = 1e-12, maxiters = 50)

Newton's method on the local equations of a law, with the Jacobian by automatic differentiation.

The step is converged when the Newton correction of the internal variables satisfies   abstol reltol , the residual being no larger than at the start of the step or the correction small at the scale of the starting state: an iterate that has run away has a small correction against its own size only, and is not accepted. The correction has the units of the internal variables whatever the units of the residual, so the residual of a law may be written in any units: the equations of a ResidualLaw in strain units and the stationarity condition of a PotentialLaw, in units of stress, converge alike. After maxiters iterations without convergence, or at a singular Jacobian, a LocalIntegrationError is thrown, which a driver catches to cut the step.

A correction that leads the internal variables where the residual is not finite (a law returns NaN outside its domain, a negative crack density for instance) is halved until the residual is finite again, at most 30 times.

MeanFieldHomogenization.MaterialLaws.LocalIntegrationError Type
julia
LocalIntegrationError(law, iterations, correction, reason)

Thrown by point_update when the local equations of a law could not be solved: the law's name, the iterations done, the last Newton correction (largest entry) and the reason. A driver that catches it can retry the step in smaller pieces, as point_response does; any other error is a bug and is not caught.

MeanFieldHomogenization.MaterialLaws.point_update Function
julia
point_update(law, g, prev, Δt; algorithm = LocalNewton(), external = (;)) -> (F, 𝕂, state)

Integrate law over a step of length Δt that brings its gradients to g, from the PointState prev. Returns the fluxes F, the consistent tangent 𝕂 = dF/dg and the new state.

g and F are flat SVectors laid out as gradient_layout and flux_layout say: for a mechanical law, the strain and the stress as Mandel vectors and 𝕂 the 6×6 tangent stiffness; for a coupled law, the fields one after the other and 𝕂 with every coupling block.

This is the function every caller ends up in: material_response wraps it for the TensND contract, point_response drives it along a loading path, an element loop calls it at each quadrature point. Its inputs and outputs are static arrays and it allocates nothing.

external gives, by name, the variables the law reads without returning their flux at the end of the step: a temperature (T = 350.0,), say. The law reads them in p, beside its parameters.

The numbers may be ForwardDiff.Dual, in g, in prev, in Δt, in the law's parameters or in external, nested or not: the fluxes, the tangent and the state then carry the derivatives, exactly.

Laws shipped ​

MeanFieldHomogenization.MaterialLaws.LinearElasticity Function
julia
LinearElasticity(ℂ)
LinearElasticity(; E, ν)

Linear elasticity,     , from a stiffness tensor (TensND, any symmetry, any basis) or, isotropic, from Young's modulus and Poisson's ratio. The control case of the potential route: its response is   with the tangent , exactly.

MeanFieldHomogenization.MaterialLaws.PowerLawElasticity Function
julia
PowerLawElasticity(; k, σ₀, ε₀, m)

Nonlinear elasticity of power-law type (deformation theory of plasticity):

linear in the volume change (bulk modulus k) and of exponent m in the deviatoric strain: the equivalent stress is   ,   being linear elasticity of shear modulus and   rigid–perfect plasticity.

Exponents below one

For   the shear stiffness is infinite at zero deviatoric strain, so the tangent there is not finite; any nonzero strain is fine.

MeanFieldHomogenization.MaterialLaws.MaxwellViscoelasticity Function
julia
MaxwellViscoelasticity(; E, ν, τ)

A Maxwell viscoelastic solid written as two potentials: a spring of isotropic stiffness (E, ν) in series with a dashpot of viscosity  , so that every mode relaxes with the same time τ:

hence    : held at a fixed strain, the stress relaxes as . The example of the potential route with an internal variable.

MeanFieldHomogenization.MaterialLaws.BiotPoroelasticity Function
julia
BiotPoroelasticity(; E, ν, b, M)

Saturated linear poroelasticity, a law of two coupled primal fields: the strain and the pore pressure , conjugate to the stress and to the change of fluid content (volume of fluid per unit volume). With the drained isotropic stiffness (E, ν), the Biot coefficient b and the Biot modulus M (see Coussy (2004)),

derived from the potential        ,   and   . Gradients (ε = 6, p = 1), fluxes (σ = 6, φ = 1); the tangent blocks are those of the Gauss-point contract, σε , σp   , φε  , φp . The example of a coupled law.

MeanFieldHomogenization.MaterialLaws.NortonViscoplasticity Function
julia
NortonViscoplasticity(; E, ν, K, n)

Isotropic elasticity (E, ν) with Norton's power-law creep of the plastic strain :

integrated by the backward Euler scheme. Under a constant uniaxial stress the creep rate is . Requires  .

MeanFieldHomogenization.MaterialLaws.J2Plasticity Function
julia
J2Plasticity(; E, ν, σ₀, H = 0, Hkin = 0, hardening = nothing)

Rate-independent von Mises plasticity with isotropic elasticity (E, ν), isotropic and linear kinematic hardening:

being the cumulated plastic strain and     the yield stress, or hardening(p) when a function is given (any law of isotropic hardening). The state is (εᵖ = 6, p = 1). An elastic predictor decides whether the step is plastic; a plastic step is solved by Newton on the backward Euler equations, which for linear hardening reproduce the radial return.

MeanFieldHomogenization.MaterialLaws.CrystalViscoplasticity Function
julia
CrystalViscoplasticity(; C₁₁, C₁₂, C₄₄, lattice = :fcc, frame = CanonicalBasis, K, n, τ₀,
                         Q = 0, b = 0, latent = 1)

Viscoplasticity of a single crystal by slip. The elasticity is cubic (C₁₁, C₁₂, C₄₄), the slip systems those of lattice (see slip_systems), and the crystal is oriented by frame, the TensND basis whose vectors are its cube axes. On each system of normal and direction , with the Schmid tensor   ,

 . A system slips only beyond its critical resolved shear stress , at a rate of power in the excess; is the slip it has accumulated, which hardens every system through the interaction matrix (latent = 1: all systems harden alike). The state is (γ = 12, v = 12); the steps are integrated by the backward Euler scheme, a step on which no system exceeds its critical stress being elastic. The flow and hardening rules have the form of the single-crystal model of Méric et al. (1991), without its kinematic hardening.

MeanFieldHomogenization.MaterialLaws.slip_systems Function
julia
slip_systems(lattice::Symbol) -> Vector{Tuple{SVector{3}, SVector{3}}}

The slip systems of a cubic lattice, as pairs (n, m) of the unit normal to the slip plane and the unit slip direction, in the frame of the cube:

  • :fcc, face-centered cubic: the 12 systems ;

  • :bcc, body-centered cubic: the 12 systems .

They are generated from their definition, every direction of the second family that lies in a plane of the first.

Nonlinear homogenization ​

MeanFieldHomogenization.MaterialLaws.NonlinearHomogenized Function
julia
NonlinearHomogenized(rve, laws; scheme = MoriTanaka(), method = ModifiedSecant())

The law of a composite whose phases follow nonlinear laws, homogenized by a secant method. rve is a cell of the package whose phases carry an isotropic stiffness :C (any value for a phase given a law: it is replaced), laws maps phase names to laws of MaterialLaws without internal variable (nonlinear elastic, e.g. PowerLawElasticity); a phase without a law keeps its linear stiffness. scheme is any linear homogenization scheme of the package, method a SecantMethod.

At the macroscopic strain , every phase with a law takes its secant moduli at its effective strain , the linear composite they form is homogenized into , and

a fixed point in the effective strains. The result is a ResidualLaw whose internal variables are those effective strains, one per phase with a law, named after the phases: the consistent tangent of the composite follows from the implicit-function theorem, and the law runs in point_response and in a finite-element computation like any other.

MeanFieldHomogenization.MaterialLaws.SecantMethod Type
julia
SecantMethod

How the effective strain of a phase is measured in a secant method: ClassicalSecant by its average, ModifiedSecant by its second moment.

MeanFieldHomogenization.MaterialLaws.ClassicalSecant Type
julia
ClassicalSecant()

The effective strain of a phase is the equivalent strain of its average strain,   , being the strain concentration tensor of the phase in the linear comparison composite.

MeanFieldHomogenization.MaterialLaws.ModifiedSecant Type
julia
ModifiedSecant()

The effective strain of a phase is the square root of the second moment of its equivalent strain,

the derivative of the effective stiffness of the linear comparison composite with respect to the shear modulus of the phase giving the energy the phase stores (Kreher, 1990). It is the method whose estimates coincide with the variational ones of Ponte Castañeda (1991) for power-law phases (see the theory page).

MeanFieldHomogenization.MaterialLaws.secant_moduli Function
julia
secant_moduli(law, εeq) -> (k, μ)

The isotropic secant moduli of a phase law without internal variable at the equivalent strain εeq, read on its own response: the shear modulus from a deviatoric strain    , of equivalent strain , and the bulk modulus from a volume change of the same size.

Damage by the growth of microcracks ​

MeanFieldHomogenization.MaterialLaws.CrackGrowthDamage Function
julia
CrackGrowthDamage(rve; scheme = MoriTanaka(), resistance, τ = nothing, families)

The law of a solid weakened by families of microcracks that grow. rve is a cell of the package whose crack phases (declared with a CrackDensity) are the families; their densities in rve are the initial ones. scheme is any linear scheme that handles cracks. resistance gives the resistance of a family as a function of its density: a CrackGrowthResistance, any function d -> R, or a NamedTuple of them keyed by family. families defaults to every crack phase of rve.

At the strain and the densities ,

This is the rate-independent law, τ = nothing. It suits a resistance that increases with the density. Under a CrackGrowthResistance, which decreases, the growth is unstable: the densities jump at the threshold, and Newton's method cannot follow the jump from the densities of the last step. A time τ > 0 then regularizes the growth by a viscosity,

 , which tends to the rate-independent law as becomes small against the duration of the loading.

The result is a ResidualLaw whose internal variables are the densities, named after the families: its consistent tangent includes the growth of the cracks, and it runs in point_response and in a finite-element computation like any other law. The cracks stay open (no unilateral contact). The approach is that of Dormieux and Kondo (2016); the derivation is in the theory page.

MeanFieldHomogenization.MaterialLaws.CrackGrowthResistance Type
julia
CrackGrowthResistance(; Gc, a₀, d₀)

The resistance to the growth of a family of penny-shaped cracks per unit volume whose radius grows at fixed number, the crack density being  . Griffith's criterion (Griffith, 1921) on each crack front, of toughness (an energy per unit area), reads in terms of the density

and being the radius and the density in a reference state (the initial one, say), which fix  . The resistance decreases as the cracks grow: the growth is unstable. Units: is a stress, e.g. Gc in N/mm and a₀ in mm for stresses in MPa.

Uniform fields ​

MeanFieldHomogenization.MaterialLaws.AbstractPhaseComposite Type
julia
AbstractPhaseComposite

Supertype of the laws of MaterialLaws built from the laws of several phases: UniformStrain and UniformStress.

MeanFieldHomogenization.MaterialLaws.UniformStrain Type
julia
UniformStrain(laws, fractions)

The composite of the phase laws in which every phase undergoes the macroscopic gradient (the strain): the Voigt bound for linear elastic phases, the Taylor model for the polycrystals. laws and fractions are NamedTuples by phase name, e.g. UniformStrain((M = matrix, I = inclusion), (M = 0.7, I = 0.3)). The phases may carry internal variables and several fields; they share their fields.

The internal variables are those of the phases, renamed after them (M_εᵖ, M_p, …). See the theory page.

MeanFieldHomogenization.MaterialLaws.UniformStress Type
julia
UniformStress(laws, fractions)

The composite of the phase laws in which every phase carries the macroscopic flux (the stress): the Reuss bound for linear elastic phases, the Sachs model for the polycrystals. Same arguments as UniformStrain.

The gradient of each phase and the macroscopic flux are unknowns of the step: the internal variables are the macroscopic flux (Σ_σ, …), the gradients of the phases (M_ε, …) and their internal variables (M_εᵖ, …). See the theory page.

The incremental tangent approach ​

MeanFieldHomogenization.MaterialLaws.IncrementalTangent Type
julia
IncrementalTangent(rve, laws; scheme = MoriTanaka())

The composite of the phase laws (a NamedTuple by phase name, one law per phase of rve) by the incremental tangent approach of Hill (1965): over a step, every phase is replaced by a linear material of its consistent tangent at the end of the step, projected on the isotropic tensors as Doghri and Ouaar (2003) do, and the linear scheme concentrates the increment of the macroscopic strain:

The geometry and the fractions come from rve; the stiffnesses it declares are replaced. The internal variables are the strains of the phases (M_ε, …) and their own internal variables (M_εᵖ, …). See the theory page.

The affine formulation ​

MeanFieldHomogenization.MaterialLaws.AffineComposite Type
julia
AffineComposite(cell, laws; scheme = MoriTanaka(), options...)

The composite of two phases of nonlinear elastic (or viscous) laws (a NamedTuple by phase name, laws of the strain alone such as PowerLawElasticity) by the affine formulation of Masson et al. (2000): each phase is linearized about its average strain,    , with the tangent and the prestress    , and the average strains are those of the linear thermoelastic comparison composite they form, which Levin's theorem gives for two phases:

cell is an RVE or a Laminate of the two phases, scheme a linear scheme for it (Laminated() for a laminate); the tangents, anisotropic in general, replace the stiffnesses the cell declares; options go to the computation of the concentration tensors (reltol of the Hill tensors of an anisotropic medium, for instance). The internal variable is the average strain of the first phase (M_ε, after its name). See the theory page.

Polycrystals ​

MeanFieldHomogenization.MaterialLaws.PolycrystalModel Type
julia
PolycrystalModel

How the grains of a Polycrystal share the macroscopic fields: TaylorModel, KronerModel, BerveillerZaouiModel.

MeanFieldHomogenization.MaterialLaws.TaylorModel Type
julia
TaylorModel()

Every grain undergoes the macroscopic strain (Bishop and Hill, 1951): the UniformStrain composite of the grains.

MeanFieldHomogenization.MaterialLaws.KronerModel Type
julia
KronerModel(; E, ν)

Every grain is a spherical inclusion, of plastic strain , in a matrix of plastic strain , the average of the grains', grains and matrix sharing the isotropic elasticity (E, ν) (Kröner, 1961). Eshelby's solution gives the interaction law      for deviatoric plastic strains,    . See the theory page.

MeanFieldHomogenization.MaterialLaws.BerveillerZaouiModel Type
julia
BerveillerZaouiModel(; E, ν)

KronerModel with a matrix that flows plastically: its stiffness in the interaction is the secant one of an isotropic elastoplastic medium,    with    ( the deviator of ) (Berveiller and Zaoui, 1978), which tempers the interaction by the plastic accommodation factor  :     . Out of proportional loading, is read as the ratio of the norms of and . See the theory page.

MeanFieldHomogenization.MaterialLaws.TangentModel Type
julia
TangentModel()

The tangent self-consistent model of Molinari et al. (1987), with the isotropic matrix of their computations: each spherical grain interacts with a matrix of the tangent viscosity of the aggregate,  , its secant viscosity. See the theory page.

MeanFieldHomogenization.MaterialLaws.Polycrystal Type
julia
Polycrystal(grains, weights = equal; model = TaylorModel())

An aggregate of grains, a vector of laws of MaterialLaws of one type (the same law with different orientations, typically a CrystalViscoplasticity per grain), with the volume weights (equal by default), sharing the macroscopic fields as model says. It runs in point_response and in the finite-element coupling like a law; its state holds the state of every grain. See the theory page.

MeanFieldHomogenization.MaterialLaws.PolycrystalState Type
julia
PolycrystalState

The state of a Polycrystal: the state of every grain (grains), and the macroscopic strain and stress it was reached at.

MeanFieldHomogenization.MaterialLaws.grain_states Function
julia
grain_states(st::PolycrystalState) -> Vector

The states of the grains, in the order of the grains of the polycrystal.

MeanFieldHomogenization.MaterialLaws.random_orientations Function
julia
random_orientations(n; rng = Random.default_rng()) -> Vector

n orientations drawn uniformly over the rotations (a random texture), as TensND bases: Euler angles in the convention of TensND, with uniform on   and , uniform on , which is the uniform measure on the rotations.

MeanFieldHomogenization.MaterialLaws.PowerLawCrystal Type
julia
PowerLawCrystal(; γ̇₀, τ₀, n, lattice = :fcc, frame = CanonicalBasis)

A rigid–viscoplastic single crystal: on each slip system of lattice (see slip_systems), oriented by frame, the slip rate is a power of the resolved shear stress,

the Schmid tensors,   . The strain rate is deviatoric and homogeneous of degree in the stress. The grain of a ViscoplasticPolycrystal.

MeanFieldHomogenization.MaterialLaws.crystal_strain_rate Function
julia
crystal_strain_rate(c::PowerLawCrystal, σ) -> SVector{6}

The strain rate of the crystal under the stress σ (Kelvin–Mandel components).

MeanFieldHomogenization.MaterialLaws.ViscoplasticPolycrystal Type
julia
ViscoplasticPolycrystal(grains, weights = equal)

An aggregate of PowerLawCrystal grains of one exponent, with the volume weights, in steady viscoplastic flow: flow_stress gives its stress under a strain rate.

MeanFieldHomogenization.MaterialLaws.flow_stress Function
julia
flow_stress(p::ViscoplasticPolycrystal, D; model = TangentModel(), algorithm = LocalNewton())
    -> (σ = Σ, grains = σ_g)

The stress deviator of the aggregate under the deviatoric strain rate D (Kelvin–Mandel components), and the stress of every grain, by the TaylorModel (every grain at the strain rate D) or the TangentModel. Iterative: numbers that can be compared only.

Strength of a frictional matrix ​

MeanFieldHomogenization.MaterialLaws.DruckerPrager Type
julia
DruckerPrager(; h, friction, dilatancy = friction)

The strength of a Drucker–Prager matrix,     , with the mean stress,   the norm of the stress deviator, h the isotropic tensile strength and friction the friction coefficient . Its plastic flow follows the potential   , dilatancy being , with   :   is the associated (normality) rule, a smaller a less dilatant flow (Maghous et al., 2009).

MeanFieldHomogenization.MaterialLaws.Heterogeneities Type
julia
Heterogeneities

What the phases other than the matrix are in a FrictionalComposite: Pores or RigidInclusions.

MeanFieldHomogenization.MaterialLaws.Pores Type

Empty pores: the phases other than the matrix carry no stress.

MeanFieldHomogenization.MaterialLaws.RigidInclusions Type

Rigid inclusions perfectly bonded to the matrix: no strain rate inside them.

MeanFieldHomogenization.MaterialLaws.FrictionalComposite Type
julia
FrictionalComposite(rve, matrix::DruckerPrager; matrix_phase = :M, heterogeneities,
                    scheme = MoriTanaka(), regularization = 1e-3)

A composite of a Drucker–Prager matrix (phase matrix_phase of rve) and heterogeneities, Pores or RigidInclusions, which every other phase of rve is. The geometry and the fractions come from rve, the linear estimates from scheme; the stiffnesses declared in rve are replaced. regularization is the parameter of the viscous regularization, relative to  , the dilatancy times the size of the macroscopic strain rate: the limit states are the limits  , which the computed ones approach at first order in .

limit_stress gives the macroscopic limit stress state reached under a macroscopic strain rate. The method and its limits are in the theory page.

MeanFieldHomogenization.MaterialLaws.limit_stress Function
julia
limit_stress(fc::FrictionalComposite, D; algorithm = LocalNewton()) -> (; σ, dv, dd)

The macroscopic stress σ (Kelvin–Mandel components) of the limit state that the composite reaches under the macroscopic strain rate D, with the effective volume and deviatoric strain rates of the matrix, dv and dd. It depends on the direction of D only. For Pores every direction gives a point of the boundary of the strength domain. For RigidInclusions the limit states form a cone whose points share one flow direction: D must satisfy the macroscopic flow rule to within the regularization, its distance to it setting the point on the cone (see the example); a strain rate that admits no limit state is refused (an ArgumentError, or a LocalIntegrationError when the iteration finds no fixed point at all). The iteration compares values, so that symbolic numbers are refused.

Nonlinear poroelasticity of cracked rocks ​

MeanFieldHomogenization.MaterialLaws.ClosingCrackPoroelasticity Function
julia
ClosingCrackPoroelasticity(rve; scheme = MoriTanaka(), closure, porosity = 0, families)

The law of a saturated rock whose microcracks close progressively under the effective stress. rve is a cell of the package whose crack phases, each with a given orientation, are the families; their densities in rve are those of the cracks all open. closure gives the open fraction of a family as a function of the normal effective traction on its plane: an ExponentialClosure, any function σn -> f, or a NamedTuple of them keyed by family. porosity is the porosity occupied by the fluid (the volume of the cracks, small for thin cracks), which sets the Biot modulus; it is held at its given value.

The microstructure sees the Terzaghi effective stress    : family has the open density    . The drained stiffness , the Biot tensor and the Biot modulus of the solid with those open densities give the rates

integrated by a backward Euler step. Gradients (ε = 6, p = 1), fluxes (σ = 6, φ = 1), internal variables the effective stress Σ′ and the fluid content φ. The derivation is in the theory page.

MeanFieldHomogenization.MaterialLaws.ExponentialClosure Type
julia
ExponentialClosure(σc)

The open fraction of a family of cracks whose closure stresses follow an exponential distribution of mean σc: under the normal effective traction on the plane of the family (negative in compression),

every crack is open under tension, and a fraction is still open under a compression equal to . The closure stress of a crack being proportional to its aspect ratio, this is an exponential distribution of the aspect ratios. Any function σn -> f, between 0 and 1, can be given instead.

Driving a point ​

MeanFieldHomogenization.MaterialLaws.LoadingPath Type
julia
LoadingPath(times; ε11 = …, σ22 = …, p = …, …)

A loading path for a material point: the times at which the response is computed, and, for each component of each field of the law, whether it is driven by its gradient or by its flux. A value is a number (held constant) or a function of time.

  • A tensor field is driven component by component: ε11, …, ε12 for the strain, σ11, …, σ12 for the stress (tensor components: ε12 is , not the engineering shear ).

  • A scalar field by its name: p for the pore pressure, φ for the change of fluid content of a BiotPoroelasticity law.

  • A component named in neither form has a zero flux: a free face.

  • external gives the variables the law reads without returning their flux, by name, as numbers or functions of time: external = (T = t -> 293 + 0.5t,) for a temperature (see point_update).

julia
tension   = LoadingPath(range(0, 1; length = 101); ε11 = t -> 2e-3 * t)   # uniaxial stress
creep     = LoadingPath(range(0, 100; length = 201); σ11 = 150.0)          # constant stress
drained   = LoadingPath(0:0.1:1; ε11 = t -> 1e-3 * t, p = 0)              # pore pressure held at zero
undrained = LoadingPath(0:0.1:1; ε11 = t -> 1e-3 * t)                     # φ = 0: no fluid exchanged

The first time is the initial instant: the response there is computed from the initial state with a step of zero length. The names are checked against the law's fields when the path is run.

MeanFieldHomogenization.MaterialLaws.point_response Function
julia
point_response(material, path::LoadingPath; algorithm = LocalNewton(), state = initial_state(material),
               rtol = 1e-10, maxiters = 25, max_cuts = 8, on_step = nothing) -> PointResponse

Drive a material point along path and return its response at every time of the path.

material is a law of MaterialLaws or any material of the Gauss-point contract (AbstractMFHMaterial). At each time the components the path drives by their flux are matched to a relative residual rtol by Newton's method on the free gradients, with the consistent tangent; when the local integration or this iteration fails, the step is retried in two halves, down to max_cuts halvings.

on_step(i, t, g, F, state), if given, is called after each time of the path, for instance to write the response to a file as it is computed, so that an interrupted run keeps what it had.

julia
law = J2Plasticity(; E = 200e3, ν = 0.3, σ₀ = 250.0, H = 1e3)
r = point_response(law, LoadingPath(range(0, 1; length = 51); ε11 = t -> 5e-3 * t))
r[:σ11], r[:εᵖ11]
MeanFieldHomogenization.MaterialLaws.PointResponse Type
julia
PointResponse

What point_response returns: the times t, and at each of them the gradients g, the fluxes F (flat vectors, laid out as the law's fields) and the state. Indexed by name:

  • r[:t], the times;

  • a component of a field, by the names LoadingPath uses: r[:ε11], r[:σ12] (tensor components, the √2 of the shears removed), r[:p], r[:φ];

  • an internal variable, by the name the law's state_layout gives it: r[:p] would be shadowed by a field of that name, so a plastic law calls its cumulated strain p only when no field does; a six-entry one also by component, r[:εᵖ11].

MeanFieldHomogenization.MaterialLaws.PointDriverError Type
julia
PointDriverError(time, iterations, residual)

Thrown by point_response when the components driven by their flux could not be matched at time, even after the step was cut as far as allowed.

Coupling to a finite-element code ​

MeanFieldHomogenization.MaterialLaws.QuadratureStates Type
julia
QuadratureStates(material, ncells, nqp)

The states of the quadrature points of a mesh, nqp per cell, kept twice:

  • states.old[cell][qp], the state at the last converged step, which the material reads;

  • states.new[cell][qp], the state the current iteration writes.

commit! makes the new states the old ones once a step has converged; rollback! discards them when it is rejected. A material never mutates a state, so nothing else is needed to retry a step.

MeanFieldHomogenization.MaterialLaws.commit! Function
julia
commit!(states::QuadratureStates)

The step has converged: the states it wrote become those the next step starts from.

MeanFieldHomogenization.MaterialLaws.rollback! Function
julia
rollback!(states::QuadratureStates)

The step is rejected: the states it wrote are discarded.

MeanFieldHomogenization.MaterialLaws.ModelingHypothesis Type
julia
ModelingHypothesis

How the strain of a finite element is a strain of the material:

hypothesiselementstrain of the material
ThreeD()3-D, coordinates the element's
PlaneStrain()2-D,    
Axisymmetric()2-D, the hoop strain   as third component

A hypothesis converts the strain of an element, and the virtual strain of each shape function, into the Mandel strain of the material, so that the element residual is   and its Jacobian    whatever the hypothesis.

MeanFieldHomogenization.MaterialLaws.ThreeD Type
julia
ThreeD()

Three-dimensional elements: the strain of the element is that of the material.

MeanFieldHomogenization.MaterialLaws.PlaneStrain Type
julia
PlaneStrain()

Two-dimensional elements in plane strain: the out-of-plane strains vanish, and the out-of-plane stress is whatever the material returns. Exact, unlike plane stress, which would have to condense the law at every point.

MeanFieldHomogenization.MaterialLaws.Axisymmetric Type
julia
Axisymmetric()

Two-dimensional elements of an axisymmetric body, in the coordinates : the axes of the material are , the hoop strain being the third component. Volumes are measured per radian, with the weight .

MeanFieldHomogenization.MaterialLaws.mandel_strain Function
julia
mandel_strain(hypothesis, ε, u = nothing, x = nothing) -> Mandel

The strain of the material for the strain ε of an element (a Tensors.SymmetricTensor{2, dim}), in Mandel components. In Axisymmetric the displacement u and the position x at the point give the hoop strain . The same function gives the virtual strain of a shape function, from its symmetric gradient and its value.

MeanFieldHomogenization.MaterialLaws.volume_weight Function
julia
volume_weight(hypothesis, x)

The factor of the volume element at the position x: in Axisymmetric (per radian), 1 otherwise (true, the unit of any number type).

MeanFieldHomogenization.MaterialLaws.quasistatic! Function
julia
quasistatic!(system!, prescribe!, u, states, times; rtol = 1e-8, maxiters = 25,
             max_cuts = 8, on_step = nothing) -> history

Drive a structure through the load steps times, by Newton's method on its equilibrium with the consistent tangent of its materials, independently of the finite-element library.

  • prescribe!(u, t) sets the prescribed degrees of freedom of u at time t;

  • system!(u, t, Δt) assembles, at the displacement u, the Jacobian K and the residual r (internal minus external forces) of the step ending at t, the rows of the prescribed degrees of freedom set to those of the identity and to zero; it reads the states from states.old and writes them into states.new, and returns (K, r);

  • u is updated in place, states a QuadratureStates.

A step is converged when   rtol times the norm of its first residual. It is then committed; otherwise, or when a material fails to integrate (LocalIntegrationError), the step is rolled back and retried in two halves, down to max_cuts halvings. on_step(i, t, u, states) is called after each time of times. Returns, for each step taken (cut steps included), the time and the norms of its successive residuals, from which the quadratic convergence of a consistent tangent can be read.

MeanFieldHomogenization.MaterialLaws.GlobalNewtonError Type
julia
GlobalNewtonError(time, iterations, residual)

Thrown by quasistatic! when the global Newton iteration did not converge at time, even after the step was cut as far as allowed.

MeanFieldHomogenization.MaterialLaws.material_element! Function
julia
material_element!(Ke, re, cv, material, hypothesis, ue, states, cell, Δt;
                  algorithm = LocalNewton(), external = (;), cache = nothing)

Assemble the element stiffness Ke and the internal-force vector re of one cell of a small-strain mechanical problem, whose material is any material of the Gauss-point contract (a law of MaterialLaws or a homogenized material), under the ModelingHypothesis hypothesis; the states of the cell's quadrature points are read from states.old[cell] and written into states.new[cell].

Methods are added by the adapter of a finite-element library; with Ferrite, using Ferrite provides one, cv being its CellValues reinitialized on the cell and ue the element's displacements.

MeanFieldHomogenization.MaterialLaws.material_system Function
julia
material_system(dh, ch, cv, material, hypothesis, states; external_forces! = nothing, …)
    -> (system!, prescribe!)

The two functions quasistatic! calls, for a finite-element library: the assembly of the Jacobian and the residual over the whole mesh with material_element!, and the application of the prescribed displacements. Methods are added by the adapter of the library; with Ferrite, dh is the DofHandler, ch the ConstraintHandler and cv the CellValues.

MeanFieldHomogenization.MaterialLaws.write_material_vtk Function
julia
write_material_vtk(path, dh, u, law, states)

Write to a VTK file the displacement u and, cell by cell, the average over the quadrature points of the committed states of a law of MaterialLaws: its fluxes (the stress, by component σ11, …), its gradients (ε11, …) and its internal variables by name (p, εᵖ11, …). Methods are added by the adapter of a finite-element library; with Ferrite, dh is the DofHandler.