Building materials in Python (`modular`)

The simcoon.modular package is the Python builder for the composable MODUL engine: a constitutive model is assembled from an elasticity definition plus any combination of mechanisms (plasticity, viscoelasticity, damage), and serialized into the self-describing props stream that the C++ engine consumes. The same object drives the in-memory solver, the file solver, or any FEA coupling — see the UMAT catalog for how MODUL fits among the other material names. Units are MPa; Voigt order is [11, 22, 33, 12, 13, 23].

Quick start

from simcoon.modular import (
    ModularMaterial, IsotropicElasticity,
    Plasticity, VonMisesYield, VoceHardening, ArmstrongFrederickHardening,
)
from simcoon import solver

mat = ModularMaterial(
    elasticity=IsotropicElasticity(C1=210000.0, C2=0.3, alpha=1.2e-5,
                                   convention="Enu"),
    mechanisms=[
        Plasticity(
            sigma_Y=300.0,
            yield_criterion=VonMisesYield(),
            isotropic_hardening=VoceHardening(Q=200.0, b=10.0),
            kinematic_hardening=ArmstrongFrederickHardening(C=20000.0, D=100.0),
        ),
    ],
)
print(mat.summary())

step = solver.StepMeca(control=['strain'] + ['stress'] * 5,
                       value=[0.02, 0, 0, 0, 0, 0], ninc=100)
res = solver.solve(step, "MODUL", mat.props, mat.nstatev)

mat.props and mat.nstatev are all any caller needs — they work identically with simcoon.umat() point evaluation, mean-field micromechanics (whose phases are passed in memory through the phases argument of solve()) and FEA couplings such as fedoo.

Elasticity

Exactly one elasticity definition per material. The elastic constants are ordinal slots C1..Cn whose meaning is fixed by the convention argument (enum, int, or string aliases):

Class	Parameters	Conventions (string aliases)
`IsotropicElasticity`	`C1, C2, alpha`	`"Enu"` (C1=E, C2=nu, default), `"nuE"`, `"Kmu"`/`"KG"`, `"muK"`, `"lambdamu"`, `"mulambda"`
`CubicElasticity`	`C1, C2, C3, alpha`	`"EnuG"` (default), `"Cii"` (C11, C12, C44)
`TransverseIsotropicElasticity`	`EL, ET, nuTL, nuTT, GLT, alpha_L, alpha_T, axis`	`"EnuG"` only (axis = isotropy axis, default 3)
`OrthotropicElasticity`	`C1..C9, alpha1..alpha3`	`"EnuG"` (E1,E2,E3,nu12,nu13,nu23,G12,G13,G23; default), `"Cii"`
`NeoHookeanElasticity`, `MooneyRivlinElasticity`, `YeohElasticity`, `IsiharaElasticity`, `GentThomasElasticity`, `SwansonElasticity`	the potential’s parameters (e.g. `C10, C20, C30, kappa` for Yeoh), `volumetric` and `alpha` keyword-only	none: the potentials of the `NEOHC`, `MOORI`, `YEOHH`, `ISHAH`, `GETHH` and `SWANH` UMATs
`HolzapfelElasticity`	`C10, k1, k2, kappa_d, fibres, kappa`, `volumetric` and `alpha` keyword-only; `fibres` is a `simcoon.Rotation`, one entry per fibre family, applied to \(\mathbf{e}_1\)	the `HOLZA` UMAT’s potential (Gasser-Ogden-Holzapfel). The only anisotropic block: see the caveats in Constitutive model (UMAT) catalog before composing it with plasticity or viscoelasticity

alpha are the thermal-expansion coefficients (per direction where applicable).

The hyperelastic blocks are not a constant stiffness: the potential is evaluated at the elastic strain it is handed, bridged by \(\mathbf{b}^{el} = \exp(2\boldsymbol{\varepsilon}^{el})\). Under finite strain that strain is the elastic logarithmic strain, and the mechanisms act additively on it. Like every other law, the block defines its stored energy per reference volume as a function of the elastic strain, \(\psi(\boldsymbol{\varepsilon}^{el}) = W(\mathbf{b}^{el}) + U(J^{el})\), and returns the Kirchhoff stress of \(W\) at \(\mathbf{F} := \mathbf{V}^{el}\). It equals \(\partial \psi / \partial \boldsymbol{\varepsilon}^{el}\) for an isotropic potential, and for the anisotropic HolzapfelElasticity only while it stays coaxial with \(\mathbf{V}^{el}\); the Stress output is \(\boldsymbol{\tau}/J\) with the total \(J\). The volumetric term only sees the elastic volume change \(J^{el} = \exp(\mathrm{tr}\,\boldsymbol{\varepsilon}^{el})\), so a free thermal expansion stays stress-free. Every potential shares its volumetric term \(U(J)\), chosen by the volumetric keyword: "log" (default) for \(U = \kappa (J \ln J - J + 1)\), "quadratic" for \(U = \frac{\kappa}{2} (J - 1)^2\). Both have \(U''(1) = \kappa\), so kappa is the ground-state bulk modulus either way; they differ away from \(J = 1\) (the quadratic one gives a pressure linear in \(J - 1\), the logarithmic one stiffens in compression and blows up as \(J \to 0\)).

Yield criteria

Used inside Plasticity; the criterion defines the equivalent stress \(\sigma_{eq}\) of the yield function \(f = \sigma_{eq}(\boldsymbol{\sigma} - \mathbf{X}) - (\sigma_Y + R(p))\):

Class	Parameters	Criterion
`VonMisesYield`	—	\(\sqrt{\tfrac{3}{2}\,\mathbf{s}:\mathbf{s}}\)
`TrescaYield`	—	Maximum shear stress
`DruckerYield`	`b, n`	Drucker \(J_2\)–\(J_3\) criterion
`HillYield`	`F, G, H, L, M, N`	Hill 1948 quadratic anisotropy
`DFAYield`	`F, G, H, L, M, N, K`	Deshpande–Fleck–Ashby (pressure-sensitive)
`AnisotropicYield`	`P11, P22, P33, P12, P13, P23, P44, P55, P66`	Full quadratic form; P admissibility requirements: see the `EPANI` row of the catalog.

Isotropic hardening

The isotropic hardening law \(R(p)\) of the accumulated plastic strain \(p\):

Class	Parameters	Law
`NoIsotropicHardening`	—	\(R = 0\) (perfect plasticity, default)
`LinearIsotropicHardening`	`H`	\(R = H\,p\)
`PowerLawHardening`	`k, m`	\(R = k\,p^m\) — for \(m < 1\) the singular onset slope is C¹-regularized below \(p = 10^{-6}\) (exact above)
`VoceHardening`	`Q, b`	\(R = Q\,(1 - e^{-b\,p})\)
`CombinedVoceHardening`	`terms = ((Q_1, b_1), ...)`	\(R = \sum_i Q_i\,(1 - e^{-b_i p})\) (standard independent sum — note this differs from the removed legacy EPCHG coupling, see the catalog)

Kinematic hardening

The back stress \(\mathbf{X}\) shifting the yield surface. The stored internal variable is the back-strain \(\boldsymbol{\alpha}\) (thermodynamic variable); \(\mathbf{X} = \tfrac{2}{3} C \boldsymbol{\alpha}\):

Class	Parameters	Law
`NoKinematicHardening`	—	\(\mathbf{X} = 0\) (default)
`PragerHardening`	`C`	Linear: \(\dot{\boldsymbol{\alpha}} = \dot{p}\,\mathbf{n}\)
`ArmstrongFrederickHardening`	`C, D`	\(\dot{\mathbf{X}} = \tfrac{2}{3}C\,\dot{\boldsymbol{\varepsilon}}^p - D\,\mathbf{X}\dot{p}\)
`ChabocheHardening`	`terms = ((C_1, D_1), ...)`	\(\mathbf{X} = \sum_i \mathbf{X}_i\), each branch Armstrong–Frederick

Mechanisms

Class	Parameters	Physics
`Plasticity`	`sigma_Y, yield_criterion, isotropic_hardening, kinematic_hardening`	Rate-independent plasticity (Fischer–Burmeister return mapping)
`Viscoelasticity`	`terms = ((E_i, nu_i, etaB_i, etaS_i), ...)`	Generalized Maxwell (Prony) branches: per branch a spring (\(E_i, \nu_i\)) in series with bulk/shear dashpots (\(\eta_B, \eta_S\)) — same rheology, layout and closed-form backward-Euler step as the kept `PRONK` kernel (identical results); the branches follow the total strain and enter the tangent exactly
`Damage`	`Y_0, Y_c, damage_type, A, n`	Scalar stiffness-degradation damage, \(\boldsymbol{\sigma} = (1-D)\,\boldsymbol{\sigma}_{eff}\), driven by the undamaged energy \(Y = \psi_0\); the other mechanisms (plastic yield included) act on the effective stress (strain equivalence). `damage_type` selects the evolution law and its extra parameters: `LINEAR` (none), `EXPONENTIAL` (`A`), `POWER_LAW` (`n`) or `WEIBULL` (`A, n`)

Multiple mechanisms compose additively on the inelastic strain; the registration order defines the statev layout (see the catalog statev section).

Tangent operator and finite strain

MODUL honors the solver’s tangent_mode (continuum or algorithmic, algorithmic being the 2.0 default — Use the solver). The algorithmic tangent of a composition is the exact derivative of its discrete update whenever the plasticity it contains is exact on its own (von Mises with isotropic hardening): the viscoelastic branches take a closed-form step and enter by the chain rule, and damage scales the effective response (plasticity yields on the effective stress \(\boldsymbol{\sigma}/(1-D)\)). With kinematic hardening or another yield criterion the plastic cutting-plane update makes it approximate (see the note on the scope of the algorithmic tangent in Use the solver). Under the finite-strain control types the composition acts as a Hencky hyperelastic law on the logarithmic strain and requires corate_type = 3 (log_R) — which is the simcoon.solver.solve() default, so nothing needs to be passed. Any other corate (e.g. corate="logarithmic", the XBM rate, code 2) is rejected under NLGEOM.

Building materials in Python (modular)