Skip to content

Nonlinear Materials and Maturity

Design Boundary

AgentFEM separates three things that commercial input decks often present together:

  • material parameters;
  • the local constitutive update;
  • the global finite-element solution procedure.

This matters for path-dependent materials. A correct radial return at one integration point is necessary, but it is not a working elastoplastic FEM solver until quadrature state, consistent tangents, increment control, convergence evidence, and restart are present.

Use constitutive.capabilities() before selecting a nonlinear law.

Compressible Neo-Hookean Hyperelasticity

constitutive.NeoHookeanProperties and model.step(...) is the canonical public entry point. It dispatches a nonlinear_static Study through an inspectable step-provider registry; the older model.hyperelastic_step(...) remains a lower-level compatibility route, not a pattern to duplicate for every material model. New constitutive families register a provider that declares which analysis/material protocol they can lower. The public model language does not acquire one method per material. The formulation uses

F = I + grad(u)
psi = mu/2 (tr(F^T F) - d) - mu ln(J) + lambda/2 ln(J)^2

and solves the first variation of total potential. Ordinary Dirichlet problems use PETSc SNES. Abaqus periodic equations use exact affine reduction and an incremental reduced Newton path. In both cases the tangent is differentiated from the same residual. This follows the same compressible Neo-Hookean energy and automatic differentiation structure as the official DOLFINx hyperelasticity demo.

constitutive.kinematics(u) exposes the standard F, C, J, and Green--Lagrange measures to reusable internals and expert workflows. Ordinary model scripts do not need to import UFL for these quantities. solvers.newton(...) supplies one public nonlinear policy for both SNES and affine-reduction paths; constraint implementation no longer changes the top-level solver language.

neo_hookean(...) in two dimensions is plane strain. neo_hookean_plane_stress(...) is a separate finite-strain membrane material: it solves a positive out-of-plane stretch satisfying P33=0 and condenses the same energy and tangent. The Study assumption must match the material; the provider rejects a mismatch before assembly.

The 3D Abaqus periodic-cell example additionally verifies quadratic tetrahedral geometry, macro-deformation control, equation mismatch, averaged stress, and quadrature-point det(F) bounds. It substitutes Neo-Hookean behavior for the unavailable Abaqus user material and does not claim constitutive equivalence.

Mooney--Rivlin Hyperelasticity

constitutive.mooney_rivlin(...) provides the three-dimensional compressible energy

psi = C10 (I1_bar - 3) + C01 (I2_bar - 3) + K/2 (J - 1)^2

constitutive.mooney_rivlin_plane_stress(...) provides the exact incompressible thin-sheet reduction reported as Eq. (17) by Wang, Fineberg, and Needleman. Both materials are consumed through the same model.step(...) language as Neo-Hookean solids. The finite-element residual, recoverable energy, first-Piola/Cauchy stress, material tangent, Explicit stable-increment estimate, and small-on-large wave analysis all use the declared energy.

The JMPS-inspired weak-interface benchmark accepts this material directly via bulk_material=...; changing the constitutive law does not copy its mesh, cohesive state, preload transfer, time integration, energy ledger, or crack observer. The three-dimensional penalty form and the exact two-dimensional incompressible reduction remain distinct capabilities; neither is presented as a general locking-free three-dimensional Explicit formulation.

J2 Plasticity

constitutive.J2LinearIsotropicHardening implements a small-strain, rate-independent Mises material-point update:

trial elastic predictor
f_trial = q_trial - (sigma_y0 + H p_old)
Delta gamma = f_trial / (3G + H)
radial correction of deviatoric stress

The implementation verifies that the corrected stress lies on the hardened yield surface, includes an exact uniaxial update, and returns the analytical algorithmic consistent tangent. The algorithmic family is the standard closest-point radial return described in the MOOSE radial-return documentation and Abaqus isotropic elastoplasticity theory.

For a three-dimensional nonlinear_static study, model.step(...) now lowers this material to a global DOLFINx path. PE and PEEQ are committed at Basix quadrature points; S and DDSDDE are trial fields updated during Newton. Failed attempts restore displacement and committed material state before automatic cutback. Complete named CellRegion assignments may dispatch different J2 parameter sets without changing the Step API. An otherwise converged attempt is also rejected when its equivalent plastic-strain increment exceeds maximum_inelastic_increment. A portable checkpoint contains displacement, accepted step coordinate, the adaptive next-increment proposal, plastic state, energy/work history, amplitude identity, and schema version.

The current boundary is explicit: 3D small strain, complete nonoverlapping material regions, natural or strong-displacement loading, structurally benchmarked MPI global equilibrium, and full-Step restart across MPI rank counts. Committed quadrature state can additionally be saved collectively and restored across MPI rank counts using physical-cell, quadrature-rule, mesh, material, and state-schema identity. A named tabular amplitude may load, unload, and reverse while the internal step coordinate remains monotone. Strong prescribed-displacement paths record generalized reaction, external work, internal energy, and balance histories. The result retains S/PE/PEEQ/MISES at constitutive integration points and adds separately named *_CELL weighted DG0 recovery fields. The recovered fields preserve element and material boundaries and are never labeled as raw integration-point values or smoothed nodal contours. Plane stress, kinematic hardening, external distributed structural benchmarks, finite-strain plasticity, and a general UMAT path remain future work.

Power-Law Creep

constitutive.PowerLawCreep provides a normalized Mises time-hardening law, exact constant-stress integration, a relaxation solution, and associative tensor increments. The constant-stress and relaxation formulas are checked against the equations used in the Abaqus creep-integration verification.

constitutive.integrate_stress_history(...) is the next reusable layer: it integrates piecewise-constant scalar or tensor stress intervals with the exact time-hardening increment and returns a named CreepHistory. It is useful for material tests and prescribed stress paths, but is explicitly not a global FE creep solver.

constitutive.ArrheniusPowerLawCreep adds a normalized temperature factor whose coefficient is calibrated at a declared reference temperature. This is the appropriate local basis for high-temperature component workflows, but it does not yet promote creep to a global coupled solver.

The global 3D route reuses the J2 quadrature transaction, implements a backward-Euler local update with analytical consistent tangent, and passes relaxation, forced-cutback, restart, and weighted field-recovery contracts. It also reproduces the official Abaqus creep_usr_creep.inp held-stress case: the published time-hardening constants, 20,000 psi stress, and 100,000 s duration give the closed-form equivalent creep strain 0.1. The versioned 100-increment contract declares its expected backward-Euler error instead of freezing a solver-specific number as physical truth.

The next milestone is nonuniform multi-element and high-temperature component evidence, followed by accepted temperature-field coupling and portable quadrature state. A new creep-law name alone is not that milestone.

Creep Damage, Sinh Flow, and Modified Theta

KachanovRabotnovCreep couples effective-stress creep and scalar damage:

epsilon_dot = A (q / sigma_ref)^n / (1 - omega)^n
omega_dot   = B (q / sigma_ref)^m / (1 - omega)^phi

For a piecewise-constant stress interval, AgentFEM integrates both equations analytically. One interval and any subdivision therefore recover the same material-point state up to floating-point tolerance. SinhCreep provides the separate hyperbolic-sine Mises rate family used when a power law is too rigid over a wide stress range.

ModifiedThetaProjection represents

epsilon = epsilon_0 + A1 (1 - exp(-alpha t))
                    + B1 (exp(alpha t) - 1)

and supplies deterministic nonnegative fitting, strain/rate projection, and a time-to-strain criterion without adding SciPy to the core. It is classified as a curve/life assessment, not a global FE stress update.

The release demo examples/creep_hot_wall_assessment.py connects an implicit heat-transfer solve, sequential thermoelastic stress, the governing sampled equivalent stress, K-R damage history, and modified-theta projection through one SimulationResult. Its material constants are explicitly illustrative. It demonstrates the power-component workflow and data contract while the global quadrature creep step remains a visible next gate.

Stress-Life Fatigue

The fatigue module is a result postprocessor:

  • Basquin and tabulated log-log S-N curves;
  • turning-point extraction and rainflow cycle counting;
  • optional linear Goodman mean-stress correction;
  • Palmgren-Miner cumulative damage.

NASA fatigue guidance describes rainflow counting as the bridge from a stress response history to cumulative S-N damage; see the NASA spectral fatigue report. AgentFEM currently accepts a scalar/equivalent stress history. Multiaxial critical-plane fatigue is a separate future capability.

constitutive.assess_history(...) returns counted cycles, Miner damage, and the life in repeated copies of that history. assess_result_history(...) accepts a named SimulationResult history and preserves its source in the assessment. This connects analysis results to fatigue without making fatigue a solver step or hiding which scalar history was used.

Verification Inventory

benchmarks.list_benchmarks() returns the current test-linked obligations. Examples teach use; benchmarks carry a criterion and automated evidence.