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) |
|---|---|---|
|
|
|
|
|
|
|
|
|
|
|
|
|
the potential’s parameters (e.g. |
none: the potentials of the |
|
|
the |
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 |
|---|---|---|
|
— |
\(\sqrt{\tfrac{3}{2}\,\mathbf{s}:\mathbf{s}}\) |
|
— |
Maximum shear stress |
|
|
Drucker \(J_2\)–\(J_3\) criterion |
|
|
Hill 1948 quadratic anisotropy |
|
|
Deshpande–Fleck–Ashby (pressure-sensitive) |
|
|
Full quadratic form; P admissibility requirements: see the |
Isotropic hardening
The isotropic hardening law \(R(p)\) of the accumulated plastic strain
\(p\). \(R\) is the hardening stress, the growth of the yield stress
(\(\sigma_Y + R(p)\), \(R(0) = 0\)); in the notation of the theory manual it is
\(F_p = \partial F/\partial p\), the derivative of the hardening function of the free
energy. Its work \(\int R\,\mathrm{d}p\) is the irrecoverable work Wm_ir:
Class |
Parameters |
Law |
|---|---|---|
|
— |
\(R = 0\) (perfect plasticity, default) |
|
|
\(R = H\,p\) |
|
|
\(R = k\,p^m\) — for \(m < 1\) the singular onset slope is C1-regularized below \(p = 10^{-6}\) (exact above) |
|
|
\(\sigma_Y + R = (\sigma_Y + B p^n)(1 + C\ln(\dot p/\dot p_0))(1 - \theta^{*m})\),
the |
|
|
\(R = Q\,(1 - e^{-b\,p})\) |
|
|
\(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 \(\mathbf{a}\) (thermodynamic variable); \(\mathbf{X} = \tfrac{2}{3} C \mathbf{a}\):
Class |
Parameters |
Law |
|---|---|---|
|
— |
\(\mathbf{X} = 0\) (default) |
|
|
Linear: \(\dot{\mathbf{a}} = \dot{p}\,\mathbf{n}\) |
|
|
\(\dot{\mathbf{X}} = \tfrac{2}{3}C\,\dot{\boldsymbol{\varepsilon}}^p - D\,\mathbf{X}\dot{p}\) |
|
|
\(\mathbf{X} = \sum_i \mathbf{X}_i\), each branch Armstrong–Frederick |
Mechanisms
Class |
Parameters |
Physics |
|---|---|---|
|
|
Rate-independent plasticity (Fischer–Burmeister return mapping) |
|
|
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 |
|
|
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). |
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.
Thermomechanical face
The same composition serves a thermomechanical block (thermomechanical=True in
simcoon.solver.solve(), sim.umat_T, the Abaqus umat_singleT) with the name
MODUL and the props [rho, c_p] + mat.props: the density and the heat capacity in
front of the stream, as the dedicated kernels. The heat of the increment and its
derivatives are written once on the mechanisms (ModularUMAT::run_T): the entropy
\(\eta = c_0 \ln(\theta/\theta_0) + \boldsymbol{\alpha}:\boldsymbol{\sigma}\), the dissipation
\(\Delta\gamma_{loc}\) of the increment, the thermodynamic split of the plastic work (the
isotropic hardening work \(\tfrac12 (R_n + R_{n+1})\,\Delta p\) stored in Wm_ir; the
back-stress energy \(\sum_i \tfrac12 (\mathbf{X}_i^n + \mathbf{X}_i^{n+1}) : \Delta\mathbf{a}_i\)
in Wm_r, since a reverse loading recovers it; the rest dissipated, Wm_d) and
\(r = \Delta\gamma_{loc} - \theta_{n+1}\,\Delta\eta\) (see Use the solver). A
composition equivalent to EPICP or EPKCP is that kernel, stress, temperature and
works to round-off, under every tangent mode and integrator and under both finite-strain
forms. The derivatives of the heat are exact for isotropic hardening (the turning of the
flow direction is in them) and hold the direction of the back stresses fixed over the
increment (Prager as EPKCP, 2e-6 on a non-proportional probe, tested at 1e-5;
Armstrong-Frederick and Chaboche 1e-4 under the cutting plane, whose explicit iterates leave the
secant \(\Delta\mathbf{a}/\Delta p\) as the derivative of the back-strain, 3e-7 under the
closest point, where the backward-Euler update has a consistent one). A viscoelastic mechanism dissipates
\(\sum_i \tfrac12 (\mathbf{A}_i^n + \mathbf{A}_i^{n+1}) : \Delta\mathbf{E}^v_i\) on its branch
forces \(\mathbf{A}_i = \mathbf{L}_i (\boldsymbol{\varepsilon} - \boldsymbol{\alpha}(\theta - \theta_0) -
\mathbf{E}^v_i)\), its heat derivatives exact (linear branches, nothing held fixed): a Maxwell
composition is PRONK to round-off, and ZENER (a Kelvin branch in series with the spring)
is the one-branch Maxwell composition with, per isotropic mode \(X \in \{K, G\}\),
\(X_1^M = X_0^2/(X_0 + X_1)\) and \(\eta^M = \eta^K X_0^2/(X_0 + X_1)^2\) (the map
its adapter applies; simcoon.modular.zener_as_maxwell() gives the branch).
A damage mechanism dissipates the trapezoid \(\tfrac12 (Y_n + Y_{n+1})\,\Delta D\) on its
driving force (exact for the linear law), its heat derivatives written on the effective stress.
A hyperelastic block carries its temperature dependence through the omega keyword of the
block, \(\Psi = f(\theta)\,\bar{W} + U(J)\) with
\(f(\theta) = 1 - \omega + \omega\,\theta/\theta_0\): omega=0 (default) is energetic
thermoelasticity, the limit of the linear blocks (the entropy
\(c_0 \ln(\theta/\theta_0) + \boldsymbol{\alpha}:\boldsymbol{\sigma}\), cooling under adiabatic
stretch); omega=1 the entropic elasticity of a rubber (\(\mu \propto \theta\), the
entropy gains \(-\bar{W}/\theta_0\), Gough-Joule heating beyond the thermoelastic inversion at
\(\lambda - 1 \simeq \alpha\theta_0\)). The volumetric term stays energetic, the heat capacity
at fixed deformation stays \(c_0\), and the mechanical block scales its isochoric stress by
\(f(\theta)\) too when the temperature moves (T_init must then be an absolute
temperature). “Entropic” names where the elastic energy is stored, not a dissipation: the block is
hyperelastic, \(\gamma_{loc} = 0\), its heat \(-\theta\,\Delta\eta\) is reversible
(Wt_r) and Wm_ir = Wm_d = 0. A load-unload cycle therefore closes without hysteresis
under isothermal and under adiabatic conditions (the heat expelled on loading is taken back on
unloading); a loop only opens under a finite heat exchange with the surroundings (thermoelastic
damping), where its area is the net heat given to the surroundings and is booked in Wm_r,
the closure \(W_m - W_m^{ir} - W_m^d\), which is not a state function when the temperature
moves. ncp="min" reaches the face through simcoon.solver.solve() and the
ncp keyword of sim.umat_T. The face serves every block and every mechanism,
omega != 0 composed with mechanisms included: the entropy then carries
\(-\bar{W}(\mathbf{H}_e)/\theta_0\) on the elastic strain
\(\mathbf{H}_e = \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{in} - \boldsymbol{\alpha}\,\Delta\theta\),
its derivatives through the chain
\(d\mathbf{H}_e = \mathbf{L}^{-1}\,(d\boldsymbol{\sigma}_{eff} - (\omega/\theta_0)\,\mathbf{T}_{iso}\,d\theta)\)
on the tangents of the effective stress, \(\mathbf{T}_{iso}\) the isochoric stress of the block
at \(f = 1\) (its stress is \(f(\theta)\,\mathbf{T}_{iso} + p_{vol}\mathbf{I}\)) (exact against finite differences with plasticity,
viscoelasticity and damage, to the accuracy of those tangents).
See also
Constitutive model (UMAT) catalog — where
MODULand the adapter-served legacy names meet, props streams and statev layoutsHands-on: writing a UMAT with simcoon — writing a dedicated UMAT by hand instead
examples/mechanical/MODUL.py— runnable gallery example (elasto-plasticity)examples/mechanical/MODUL_finite.py— the same composition under NLGEOMexamples/mechanical/MODUL_hyper_visco.py— Yeoh block + Prony branches, finite-strain viscoelasticity (rate sweep, relaxation, energy split)