Physics models and solvers

Models

Physics models that the package ships, as opposed to those a script defines for itself. What earns a model its place here is that someone else would write it again: the equation is fixed and only the data changes, so shipping it turns the next script into a case rather than a re-implementation.

That is not the same as being complicated. Fick diffusion is a storage term and a flux, and Writing a model measures that writing such a model through the package costs three lines more than writing it directly against VoronoiFVM — the abstraction pays for itself in what it shares, not in what it saves per model. A model whose equation is particular to one study is still better written where it is used.

Every model carries its boundary data in a dirichlet field rather than in a method, because an imposed value is what distinguishes one case from another and not one model from another. A value may be a number or a function of time.

PoroMechanics.dirichlet_value Function

dirichlet_value(value, t) -> Real

Resolve one entry of a model's dirichlet data at time t. A number is imposed as it stands; anything else is called, so a boundary that ramps or cycles is written (region, t -> p_top * min(1, t / t_c)) and needs no method of its own.

The distinction is deliberate — a number cannot be called and a schedule cannot be added, so dispatch separates them with no run-time test and no allocation. The value is returned as given, which keeps a ForwardDiff.Dual boundary value dual and lets a result be differentiated with respect to what is imposed on the boundary.

source

PoroMechanics.apply_dirichlet! Function

apply_dirichlet!(f, u, bnode, dirichlet; species = 1)

Impose every entry of a model's dirichlet data at boundary node bnode, each value resolved at the current time by PoroMechanics.dirichlet_value.

A tuple is consumed head-and-tail rather than in a loop, which is not a stylistic choice: ((1, 0.0), (2, t -> ramp(t))) is a heterogeneous tuple, a plain for over it has no single element type, and the boxing that follows shows up as allocations inside the assembly loop — the one place in a finite volume code where they are paid per facet and per Newton iteration. Recursing on Base.tail is unrolled at compile time, so each entry is specialized on its own type. Anything that is not a tuple falls back to the loop.

source

Fickian diffusion

PoroMechanics.FickModel Type

FickModel(; phi, D, dirichlet)

Diffusion of one solute in a saturated porous medium. The unknown is the concentration [mol/m³].

field	meaning
phi	porosity [-] — a scalar, or one value per cell region
D	effective diffusion coefficient [m²/s] — a scalar, or one value per cell region
dirichlet	imposed concentrations, as ((region, value), …)

Boundaries not named are sealed: zero flux is VoronoiFVM's default and the usual meaning of an unlisted boundary in a transport problem.

dirichlet carries the boundary data rather than a method, because that is what distinguishes one case from another and not one model from another — see RichardsModel for the same argument at greater length. A value may be a number or a function of time, resolved by PoroMechanics.dirichlet_value.

A layered barrier is likewise a case, not a second model: pass a collection for phi and D, indexed by the cell region.

FickModel(; phi = [0.30, 0.12], D = [1.0e-10, 2.0e-12], dirichlet = ((1, 1.0),))

The coefficients are type parameters, so a profile can be differentiated with respect to D and not only with respect to the unknowns — which is what an inverse identification of a diffusion coefficient from a measured profile needs.

source

PoroMechanics.diffusivity Function

diffusivity(m::FickModel, region) -> D

Effective diffusion coefficient of cell region region, in m²/s.

source

Darcy flow

PoroMechanics.DarcyModel Type

DarcyModel(; k_int, mu_l, storativity, dirichlet)

Linear single-phase flow in a saturated porous medium. The unknown is the pore pressure [Pa].

field	meaning
k_int	intrinsic permeability [m²] — a scalar, or one value per cell region
mu_l	dynamic viscosity [Pa·s]
storativity	storage coefficient [Pa⁻¹] — a scalar, or one value per cell region
dirichlet	imposed pressures, as ((region, value), …)

Boundaries not named are impermeable. A value in dirichlet may be a number, or a function of time — which is how a pressure is ramped rather than stepped:

DarcyModel(; dirichlet = ((1, 0.0), (2, t -> p_top * min(1, t / t_c))))

The ramp is a property of the case, not of Darcy's law, so it belongs in the data and not in a method. Its practical use is to remove the discontinuity between a zero initial condition and an imposed pressure, which otherwise forces the step-size controller down to Δt_min at startup: a Dirichlet condition is applied unconditionally, so the jump it creates does not shrink when Δt does.

A layered column is a case as well: pass a collection for k_int or storativity, indexed by the cell region.

source

PoroMechanics.storativity Function

storativity(m::DarcyModel, region) -> S

Storage coefficient of cell region region, in Pa⁻¹.

source

PoroMechanics.mobility Function

mobility(m::DarcyModel, region) -> k/μ

The coefficient in front of in Darcy's law, in m²/(Pa·s).

source

mobility(m::PoroplastModel, p) -> k_h

, the coefficient in front of in the mass flux.

source

Richards' equation

PoroMechanics.RichardsModel Type

RichardsModel(; phi, rho_l, k_int, mu_l, p_g, gravity, gravity_axis, retention, rel_perm, dirichlet)

Unsaturated single-phase flow. The unknown is the liquid pressure [Pa].

field	meaning
phi	porosity [-]
rho_l	liquid density [kg/m³]
k_int	intrinsic permeability [m²] — a scalar, or one value per cell region
mu_l	dynamic viscosity [Pa·s]
p_g	gas pressure, held constant [Pa]
gravity	signed component of gravity along gravity_axis [m/s²]
gravity_axis	coordinate axis gravity acts along: 1 for , 2 for , 3 for
retention
rel_perm
dirichlet	imposed pressures, as ((region, value), …)

gravity_axis defaults to 1, which is the only orientation a 1D column along can have. It has to be said out loud in two or more dimensions: a vertical column meshed in the plane is driven by gravity along , and projecting it on silently removes the drainage rather than failing.

dirichlet carries the boundary data because that is what distinguishes one case from another, not one model from another: a Richards model with the pressure imposed on the right is the same physics as one with it imposed on the left. Hard-wiring a region number into the model would make the model unusable for the next problem, which is exactly what happens when a physics model is written inside a script.

Boundaries not named are no-flow, which is VoronoiFVM's default and the usual meaning of an unlisted boundary in a flow problem. A value may be a number or a function of time — PoroMechanics.dirichlet_value resolves it — so a pressure that ramps or cycles is data like any other.

A layered or composite medium is likewise a case, not a second model: pass a collection for k_int, indexed by the cell region, and the same struct covers it.

RichardsModel(; k_int = [8.9e-12, 8.9e-13], gravity = -9.81, gravity_axis = 2, …)

source

PoroMechanics.liquid_conductivity Function

liquid_conductivity(m::RichardsModel, pc, region = 1) -> ρ_l k_int k_rl(p_c) / μ_l

The coefficient in front of in the mass flux, in kg/(m·s·Pa). region selects the permeability when the model carries one per cell region.

source

Saturation from a solution

For Richards flow, saturation follows from the liquid pressure. For drying, it also depends on the air pressure, temperature and material region.

PoroMechanics.liquid_saturation Function

liquid_saturation(m::RichardsModel, pc)

from the model's retention curve.

source

liquid_saturation(model::DryingModel, u, region)
liquid_saturation(model::DryingModel, profiles, node, region)

Liquid saturation from a local (p_l, p_a, T) state, or column node of a profile matrix. Applies the Kelvin relation and thermal shift before evaluating the region's retention curve. Pressures can be continuous across a material interface while saturation differs between its two sides.

source

Non-isothermal drying

Configure the material laws, shared fluid properties and boundary values, then pass the model and your grid to fvm_system. Each cell region indexes materials. The three rows contain liquid pressure, dry-air pressure and temperature; the balances store water mass, dry-air mass and entropy.

PoroMechanics.DryingMaterial Type

DryingMaterial(; phi, k_int, lam_s, C_s, retention, rel_perm)

Properties of one rigid porous material: porosity [-], intrinsic permeability [m²], solid thermal conductivity [W/(m·K)], volumetric solid heat capacity [J/(m³·K)], retention and liquid relative-permeability laws. Numerical coefficients are promoted so that material parameters can carry ForwardDiff.Dual values.

source

PoroMechanics.DryingParameters Type

DryingParameters(; kwargs...)

Shared water, vapor, dry-air and thermal-retention coefficients for DryingModel. Defaults describe the existing non-isothermal drying model. All pressures are in Pa, temperature in K, densities in kg/m³, viscosities in Pa·s, conductivities in W/(m·K), heat capacities in J/(kg·K), and latent heat in J/kg. M_vsR and M_asR are molar mass divided by the gas constant [kg·K/J]. D_av0 is the reference air-vapor diffusivity [m²/s], H_a the Henry constant [Pa], and alpha_T the thermal-retention coefficient [1/K].

Use positive dry-air pressure and temperature. Thermal retention requires 1 - alpha_T * (T - T_0) > 0. These are constitutive-domain restrictions, not a guarantee that Newton trial states will remain admissible.

source

PoroMechanics.DryingModel Type

DryingModel(; materials, parameters=DryingParameters(),
            dirichlet=((), (), ()), heat_flux=())

Non-isothermal water and air transport in a rigid porous medium. The unknowns are liquid pressure p_l [Pa], dry-air pressure p_a [Pa], and temperature T [K]. Storage contains total water mass, total dry-air mass (including dissolved air), and volumetric entropy. Fluxes combine Darcy flow, vapor diffusion, dissolved-air advection, and heat conduction. Gravity and mechanical deformation are omitted. The entropy balance retains the reference model's approximation without a separate volumetric entropy-production term.

materials[region] supplies a DryingMaterial for each cell region; a tuple, vector or indexable mapping can be used. Geometry, initial conditions and time controls are supplied separately to the grid and solver.

dirichlet contains three boundary tuples, one per unknown, with entries (boundary_region, value); a value is a number or t -> value. heat_flux contains (boundary_region, Q) entries with incoming heat flux Q in W/m², also constant or time-dependent. Heat is imposed as entropy flux Q/T. The denominator is floored at 200 K during trial evaluations, as in the original model; this is a numerical guard, not a physical extension to low temperatures. Unspecified boundaries have zero flux. Use disjoint temperature-Dirichlet and heat-flux boundaries. The time integrator must stop at discontinuities of Q.

source

PoroMechanics.drying_material Function

drying_material(model, region)

Material on the specified cell region. At an interface, query each side separately: VoronoiFVM assembles its storage contributions with their respective region labels.

source

PoroMechanics.vapor_pressure Function

vapor_pressure(parameters::DryingParameters, p_l, T)

Vapor pressure [Pa] from the modified Kelvin relation, including the temperature variation of latent heat. At (p_l0, T_0) it returns p_v0.

source

vapor_pressure(model::DryingModel, p_l, T)

Evaluate the Kelvin relation using the model's shared DryingParameters.

source

The drying tutorial supplies a clay/rock case. Its drying_case factory keeps geometry, initial conditions and heating history explicit, and run_drying(case = custom_case) runs a modified experiment without redefining any balance callback.

Reactive transport

NernstPlanck transports ion concentrations and an electric potential. EquilibratedTransport transports conserved component totals and obtains local aqueous concentrations from ChemistryLab. These are distinct formulations: the equilibrium-coupled model currently omits electromigration. See Reactive transport for units, boundary data and optional chemistry setup.

PoroMechanics.NernstPlanck Type

NernstPlanck(; phi, D, z, tortuosity=nothing, q_background=0, dirichlet=())

Multi-ionic transport with a zero-current closure.

field	meaning
phi	porosity, scalar
D	one diffusivity per ion; effective in bulk unless multiplied by tortuosity
z	charge numbers, one per ion, in a tuple
tortuosity	an AbstractTortuosity; τ here is D_eff/D⁰, not a geometric factor
dirichlet	one boundary tuple per unknown, ions first and the potential last
q_background	prescribed Σ zᵢ cᵢ [mol/m³ of solution], opposite to the untransported charge

The unknowns are the n ion concentrations followed by the potential Ψ, so nspecies == length(D) + 1.

Coefficients retain their numeric types, so a ForwardDiff.Dual can enter a parameter. Use fvm_system(model, grid; reaction = true) to assemble electroneutrality. Initial and boundary concentrations must satisfy that constraint. Set a Dirichlet reference for the potential to remove its constant nullspace.

The tortuosity convention

OhJang returns D_eff/D⁰, porosity included. The flux carries no extra φ, while the storage does. Mixing that up with FickModel, whose D is a pore diffusivity, costs a factor 1/φ.

source

PoroMechanics.nions Function

Number of transported ions, excluding the electric potential.

source

PoroMechanics.ipot Function

Index of the dimensionless electric potential in the unknown vector.

source

PoroMechanics.effective_diffusivity Function

effective_diffusivity(m, i)

Dᵢ τ(φ), the coefficient the flux actually carries.

source

PoroMechanics.net_charge Function

net_charge(m, u) -> Vector

Σᵢ zᵢ cᵢ at every node, to be compared with q_background. The algebraic constraint preserves this value when initial and boundary data are consistent.

source

PoroMechanics.edge_current Function

edge_current(m, u, dx) -> (current, scale)

Σᵢ zᵢ Jᵢ on each edge of a uniformly spaced one-dimensional grid, and the sum of |zᵢ Jᵢ| to scale it by. dx is the node spacing. At zero-current conditions the first vector should vanish to solver accuracy.

source

PoroMechanics.ComponentSet Type

ComponentSet(; names, z, D, fixed)

Transported primaries, charges and diffusivities. fixed is a named tuple of all remaining primary totals (mol/m³ of medium). The charge row Zz defaults to zero.

source

PoroMechanics.ncomp Function

Number of transported components, excluding fixed totals.

source

PoroMechanics.EquilibratedTransport Type

EquilibratedTransport

Transport of conserved component inventories [mol/m³ of medium]. Construct it with equilibrated_transport, which checks the complete chemical basis and seed. Storage is the inventory itself. Fluxes use aqueous component concentrations from certified local equilibrium, with one Fickian diffusivity per component.

This reduced model omits electromigration; components.z records charge numbers but does not close a current balance. Use NernstPlanck for explicit ionic migration. No mineral or surface term is added outside the conserved inventories.

source

PoroMechanics.equilibrated_transport Function

equilibrated_transport(components, system; initial_state, phi=0.121,
                       tortuosity=nothing, dirichlet=())

Build a model with a complete basis and a nonnegative starting composition for one m³ of medium. initial_state supplies T, P and a starting guess, never the imposed inventories. Specify water and aluminum through fixed when held constant.

Requires using ChemistryLab, DynamicQuantities to activate the optional extension. Also load OptimaSolver, which provides ChemistryLab's certified solver. Use the patched ChemistryLab environment prepared by scripts/prepare_chemistrylab.jl. ChemistryLab owns the equilibrium solve, its certificate and implicit sensitivities. D denotes a bulk effective diffusivity, or a free-water diffusivity multiplied by tortuosity = D_eff/D⁰; no extra porosity factor is applied to the flux.

source

PoroMechanics.component_totals Function

Complete signed totals, including fixed components and the charge row.

source

PoroMechanics.equilibrium_state Function

equilibrium_state(m, totals) -> (state, certificate)

Solve with explicit signed totals. Dual-valued seed amounts select ChemistryLab's implicit differentiation route; their values remain a physical starting composition. ChemistryLab evaluates the certificate on the primal solution.

source

PoroMechanics.speciate Function

speciate(m, totals) -> (aqueous_component_concentrations, ok)

Concentrations are in mol/m³ of solution, summed over all aqueous species. Uncertified answers return ok=false; flux! rejects them. Solver exceptions propagate to VoronoiFVM, whose transient controller can retry a smaller time step.

source

PoroMechanics.LocalEquilibriumError Type

Raised when a transport flux would use an uncertified local equilibrium.

source

Linear Biot assembly and time integration

Keep the Ferrite mesh, spaces and constraints explicit, select a material for each cell, then assemble once and integrate at the requested times. A step callback can record fields or diagnostics without implementing the time loop.

PoroMechanics.assemble_biot_matrices Function

assemble_biot_matrices(dh, cv_u, cv_p, material; constraints=nothing, valuetype=...)
    -> (K1, K2)

Assemble the stationary and storage matrices of linear Cartesian Biot poroelasticity. material is a BiotPoroelastic, or a function cell -> material for a heterogeneous medium. The closed Ferrite DofHandler must contain exactly :u followed by :p, matching the supplied displacement and pressure CellValues. Use one cell type and interpolation pair throughout the mesh.

Pass the closed ConstraintHandler as constraints when using affine constraints: its master/slave couplings must be present in the sparsity pattern before assembly. The returned matrices are unconstrained and share that pattern.

valuetype defaults to the material's numeric type for a homogeneous medium, and Float64 for a material selector. For parameter differentiation through a selector, pass the promoted coefficient type explicitly. The CellValues are reinitialized in place; the material data and constraints are not changed.

source

PoroMechanics.assemble_biot_load Function

assemble_biot_load(dh, facets, fv_u, model; valuetype=Float64) -> Vector

Integrate mechanical surface loads on facets using facet_load!. The callback receives a zeroed local displacement vector and reinitialized Ferrite FacetValues; it adds the case's traction contributions. They are scattered into field :u of the global load vector, leaving the pressure entries zero.

For several loaded boundaries, sum the returned vectors. A time-dependent load can be rebuilt inside the load(t) callback of solve_biot. Specify valuetype when the loading carries a numeric type other than Float64.

source

PoroMechanics.solve_biot Function

solve_biot(K1, K2, constraints; inival, times, load=nothing, on_step=nothing, linsolve=\)
    -> Vector

Solve K2 * dx/dt + K1 * x = load(t) by backward Euler at each supplied time. times includes the initial time and must be finite and strictly increasing. The matrices are constant, unconstrained sparse matrices with the same sparsity pattern, including any affine-constraint couplings. assemble_biot_matrices builds this pattern when given constraints.

load is a constant vector, a function t -> vector evaluated at each new time, or nothing for zero loading. Allocate the matrices and initial values in a type that can hold all coefficients and loads. The default linear solver is \; linsolve(A, rhs) can supply a different solver (for example, a dense solve for small systems carrying ForwardDiff.Dual values).

Initial Dirichlet and affine constraints are applied to a copy of inival. At every step the constraint values are updated, the matrix is rebuilt as K1 + K2 / dt, and the constrained system is solved. Affine slave values are reconstructed after solving. Factorizations are not cached.

on_step(x, t, step) runs after each solve, with step starting at 1. It can record diagnostics or save copy(x) for a time history; treat x as read-only, since its buffer is reused. Only the final vector is returned. The input matrices, initial vector, and load vectors are preserved; the constraint handler's time and values are updated in place.

source

Poroplasticity

One-dimensional axisymmetric poroplasticity: a Richards-like liquid balance coupled to a skeleton that may yield. The state carried between steps is the material's own, so any AbstractMaterial can be the skeleton.

PoroMechanics.PoroplastModel Type

PoroplastModel(; material, phi0, rho_l0, k_l, p_l0)

field	meaning
material	a BiotPlastic: the skeleton, b, beta, N, k, mu_l
phi0	porosity at the reference state [-]
rho_l0	fluid density at p_l0 [kg/m³]
k_l	fluid bulk modulus [Pa]
p_l0	reference pore pressure [Pa]

The fluid is compressible,     , which is what makes the storage term more than the porosity change.

source

PoroMechanics.PoroplastState Type

PoroplastState(solid, p)

What one quadrature point carries between steps: the skeleton's own state — stress, strain, plastic strain, prestress — and the pore pressure it was last converged at.

The pressure is here rather than derived because two terms need the previous one: the storage term needs the mass it implies, and the mobility is evaluated there rather than at the current iterate. That second choice is Bil's, not a simplification — k_h = val_n.Permeability_liquid in Poroplast.cpp, exactly as in its Richards.

source

PoroMechanics.poroplast_initial_states Function

poroplast_initial_states(m, nodes, σ0_total) -> Vector{Vector{PoroplastState}}

Two quadrature states per element, from the total initial stress a deck quotes.

The conversion lives here because it is the step that is easy to get wrong and impossible to notice: the skeleton's prestress is the effective stress

and base/Poroplast quotes sigma_0 = -11.5 MPa with p_0 = 4.7 MPa, so the skeleton starts at -7.74 MPa, not -11.5. Handing the total stress straight to the skeleton leaves the initial state out of equilibrium by beta*p_0 — 3.76 MPa here — and the first step then produces a perfectly smooth, entirely spurious wave of displacement and pressure.

source

PoroMechanics.poroplast_element_residual Function

poroplast_element_residual(m, dofs, dofs_n, r1, r2, states_n, Δt) -> (residual, states)

Residual of one two-node element, and the quadrature states it leaves behind.

dofs is [u₁, p₁, u₂, p₂]. Two-point Gauss quadrature, with the axisymmetric weight r folded into each contribution.

The returned residual is the internal one; boundary tractions are added by the caller, which is where they belong — a traction is a property of the problem, not of an element.

source

PoroMechanics.poroplast_step! Function

poroplast_step!(m, nodes, states, dofs, Δt; σ_inner, σ_outer, p_inner, p_outer) -> dofs

One implicit step of the coupled problem, by Newton.

The element Jacobian is ForwardDiff.jacobian of poroplast_element_residual — four degrees of freedom per element, so differentiating it costs almost nothing and it cannot disagree with the residual it linearizes. The quadrature states are frozen at their converged values during the differentiation, which is what makes the tangent consistent rather than continuum.

The global matrix is dense. With two unknowns on a hundred elements that is a 202×202 factorization, far cheaper than the sparse bookkeeping it would replace, and this model is one-dimensional by construction.

σ_inner and σ_outer are the radial stresses applied at the two ends, compression negative. They enter as boundary work -σ r δu at the inner face, whose outward normal points inward, and +σ r δu at the outer one.

Convergence is measured on the Newton increment, scaled per field, which is Bil's own criterion — its Objective Variations block names exactly these two scales, u_1 = 1e-3 and p_l = 1e5 for this deck.

Neither obvious alternative works, and both were tried. A single norm of the residual tests whichever block carries the larger units: the mechanical one is a force per radian of order σ r ≈ 10⁸ while the hydraulic one is a mass rate of order 10⁻⁷, so ‖R‖ < 10⁻⁹ is not strict, it is a mechanical test with the hydraulic block along for the ride. Scaling that norm by the applied traction inverts the failure rather than removing it: after the cavity pressure is released the hydraulic residual is fifteen orders of magnitude below the mechanical scale, every step converges at the first iteration, and the pressure field simply stops evolving — a smooth, plausible, entirely frozen answer.

source

PoroMechanics.fluid_density Function

fluid_density(m::PoroplastModel, p) -> ρ_l

   .

source

PoroMechanics.liquid_mass Function

liquid_mass(m::PoroplastModel, ε, εp, p) -> m_l

 , with the porosity from porosity — so an elastic and a plastic volume change enter it with their own coefficients.

source

PoroMechanics.intrinsic_permeability Function

intrinsic_permeability(m::RichardsModel, region) -> k

Permeability of cell region region. Dispatch rather than a runtime branch, so a single-region model pays nothing for the generality and a Vector{<:Dual} still differentiates.

source

PoroMechanics.axisymmetric_strain_1d Function

axisymmetric_strain_1d(du_dr, u, r) -> SymmetricTensor{2,3}

   — the strain of a radial displacement field under plane-strain conditions along the axis.

source

Computational homogenization

A periodic cell solved under an imposed macroscopic strain, or under an imposed macroscopic stress by an outer Newton loop on the strain that produces it. The cells may be plastic, in which case the tangent is the algorithmic one and the state is carried between macroscopic steps.

PoroMechanics.PeriodicCell Type

PeriodicCell(grid, materials; ip_order = 1, qr_order = 2)

A unit cell with one material per cell region, to be homogenized under periodic boundary conditions.

materials maps a cell region index to an AbstractMaterial. Regions are the grid's own cellsets-equivalent — with a mesh read from Gmsh, the compacted elementary tags.

Plane strain in two dimensions: the out-of-plane strain is zero, which is what a cell extruded along a third axis does and what the codes this is checked against assume.

source

PoroMechanics.periodic_cell Function

periodic_cell(nodes, cells, regions, materials; tol = 1e-8) -> PeriodicCell

Build a cell from raw connectivity: nodes is dim × n, cells is nodes_per_cell × m, regions one index per cell.

The periodic facet pairs are collected geometrically, from the coordinates, rather than from named boundary sets. A mesh written by another code carries its own naming convention — Bil's composite0.msh numbers its edges 13, 14, 104, 105, … with no left/right meaning — and matching on position is the one rule that survives crossing that boundary.

source

PoroMechanics.cell_states Function

cell_states(cell) -> Matrix

Fresh material states, one per quadrature point of every element: states[q, e].

A path-dependent cell needs them; an elastic one is handed NoState() and never looks.

source

PoroMechanics.homogenize_stress Function

homogenize_stress(cell, ε_macro, states_n = nothing; maxiter = 30, rtol = 1e-10)
    -> (σ_macro, u, states)

Solve the cell under the imposed mean strain and return the volume-averaged stress, the fluctuation field, and the quadrature states it leaves behind.

Newton on the fluctuation, so a path-dependent phase is handled by the same routine as a linear one — an elastic cell simply converges on the first correction. states_n are the states the cell converged at previously; nothing starts from fresh ones, which is only right for a first step or a rate-independent elastic cell.

Convergence is measured on the residual relative to the force the macroscopic strain alone would generate, not absolutely. The residual of a cell is an integrated stress and its magnitude follows the loading; an absolute tolerance would be a statement about the units of the moduli.

source

PoroMechanics.homogenize_to_stress Function

homogenize_to_stress(cell, σ_target, states_n; ε_guess, maxiter = 25, rtol = 1e-8)
    -> (; ε, σ, states, iterations)

Find the macroscopic strain that puts the cell under a prescribed in-plane stress, in plane strain — the control a single-element two-scale problem reduces to.

The three in-plane components are stress-controlled and  , so comes out as a result. Newton on the three unknowns, with the homogenized tangent obtained by perturbing the cell — three extra cell solves per iteration. That is not a shortcut: once a phase yields, the tangent of the cell is not the volume average of anything, and Bil homogenizes its own by finite differences for the same reason.

states_n are the states at the last converged macroscopic step; every trial re-solves the cell from those, never from the previous trial, so the path stays single-valued.

source

PoroMechanics.homogenized_stiffness Function

homogenized_stiffness(cell::PeriodicCell) -> Matrix

The effective in-plane stiffness of a linear cell: homogenized_tangent evaluated at zero strain, in the (11, 12, 22) convention.

Valid only while every phase is linear. A cell with a plastic phase has no such matrix — its tangent depends on where it is on its path, and it has to be asked for there.

source

PoroMechanics.homogenized_tangent Function

homogenized_tangent(cell, ε_macro, states = nothing; perturbation = 1e-8) -> Matrix

The 3×3 in-plane macroscopic tangent at a given macroscopic strain, by forward differences on the cell — one column per component of (ε₁₁, ε₁₂, ε₂₂), and the same convention on the rows, so no engineering factor of two enters anywhere.

Finite differences and not an assembled quantity, deliberately. While the cell is elastic the tangent could be condensed out of the microscopic stiffness, but the moment a phase yields there is nothing to condense: the macroscopic tangent depends on which quadrature points are on their yield surface, and only a perturbation sees that. Bil homogenizes its own the same way.

source

PoroMechanics.plane_strain Function

plane_strain(ε2::SymmetricTensor{2,2}) -> SymmetricTensor{2,3}

Embed a plane strain tensor in three dimensions, with  .

source