The material is explicit and typed. There is no integer constitutive selector and no shared bag of parameters: every law names its own arguments, and every argument carries its units.
An InstantaneousMaterial can contain an elastic law, a viscous law, or both:
from pyimr import CarreauYasuda, Gent, InstantaneousMaterial
material = InstantaneousMaterial(
elastic=Gent(shear_modulus_pa=2500.0, extensibility=250.0),
viscous=CarreauYasuda(
zero_shear_viscosity_pa_s=0.5,
infinite_shear_viscosity_pa_s=0.02,
time_constant_s=20e-6,
transition_exponent=2.0,
power_index=0.45,
),
)Elastic: NeoHookean, MooneyRivlin, Yeoh, Fung, Gent, ArrudaBoyce,
Ogden.
Generalized-Newtonian: Newtonian, PowerLaw, CarreauYasuda, Cross,
PowellEyring, ModifiedPowellEyring, HerschelBulkley, Bingham.
Carreau is CarreauYasuda(transition_exponent=2); simplified Cross is
Cross(transition_exponent=1).
Ogden takes matched tuples and is the only elastic law here that depends on
the principal stretches rather than on I1 alone, so exponents may be negative
or fractional:
material = Ogden(shear_moduli_pa=(1800.0, 600.0, -300.0), exponents=(1.3, 4.0, -2.0))A single term with exponents=(2.0,) is neo-Hookean, and reduces to it
exactly rather than asymptotically. The small-strain shear modulus is
sum(shear_moduli_pa * exponents) / 2, which must be positive; individual
moduli may be negative, as above.
The standard Blatz-Ko strain energy is distinguished by its dependence on I3,
and this solver assumes an incompressible spherical deformation with stretches
(l^-2, l, l), so I3 = 1 identically. In that limit Blatz-Ko is
MooneyRivlin(c10=0.0, c01=mu/2) -- verified equal to machine precision -- so a
separate class would be an alias, not a new capability. A genuinely
compressible Blatz-Ko needs the incompressibility assumption relaxed throughout
the radial dynamics.
The pair uses the standard laws, eta_inf + (eta_0 - eta_inf) * asinh(x)/x and
its log1p(x)/x variant with x = lambda*|gdot|. Both reduce exactly to
Newtonian(eta_0) as lambda -> 0. IMRv2's f_viscosity.m instead uses
sinh(x)/x^nc, which is shear-thickening and diverges exponentially, and a
log(1+x)/x^nc variant with no finite zero-shear limit unless nc == 1;
neither was copied. See upstream.md.
PowerLaw, HerschelBulkley and Bingham require a positive
regularization_rate_per_s. The latter two use a smooth yield-stress
regularization so the implicit radial equations retain a finite tangent at zero
strain rate.
Gent lock-up is a material-domain error: if a trajectory reaches
SimulationError rather than continuing with
nonphysical stress.
Prepared Gauss-Legendre rules evaluate the finite-interval stress integrals. The
solver evaluates the stress-rate terms and acceleration coefficient
analytically, including the viscosity tangent -- no finite-difference derivative
is used inside the radial dynamics. The specialized NeoHookeanKelvinVoigt path
and the equivalent composable material agree to solver tolerance.
The finite-dimensional hot paths are Zener, QuadraticZener, CubicZener,
QuadraticKelvinVoigt, OldroydB and LinearMaxwell:
material = Zener(
shear_modulus_pa=2500.0,
viscosity_pa_s=0.1,
relaxation_time_s=40e-6,
retardation_time_s=8e-6,
)QuadraticKelvinVoigt is NeoHookeanKelvinVoigt with a stiffening term
quadratic in I1 - 3, which is Yeoh truncated after c2; QuadraticZener
adds the Maxwell arm to it, and Zener is QuadraticZener at zero stiffening.
CubicZener carries the next term of the same expansion, with cubic = 0
reducing to QuadraticZener exactly. On measured collapses I1 - 3 reaches
24 to 119, where the cubic term is not a small correction. LinearMaxwell is
Zener without the parallel spring: no modulus, no retardation, an elastic
target of zero.
Each Zener above carries its own hand-derived stress integral, so relaxation
was available to the neo-Hookean family and to nothing else. Gent, Fung,
Arruda-Boyce, Mooney-Rivlin, Yeoh and Ogden existed only as
InstantaneousMaterial, elastic with no memory, which is the wrong comparison
to make against a relaxing model. RelaxingMaterial is the Zener construction
freed from one potential: any ElasticModel as the equilibrium target, with a
Maxwell arm on top.
from pyimr import RelaxingMaterial, Yeoh
material = RelaxingMaterial(
elastic=Yeoh(c1_pa=1250.0, c2_pa=100.0, c3_pa=10.0),
viscosity_pa_s=0.1,
relaxation_time_s=40e-6,
retardation_time_s=8e-6,
)The equilibrium target is taken by the same quadrature InstantaneousMaterial
uses rather than by a closed form derived per law, so one class covers every
elastic model in the package and any added later. That costs the quadrature at
each step (quadrature_points, default 32) and buys a controlled comparison.
With elastic=NeoHookean(...) it agrees with Zener; the closed form is the
one to use where it applies.
A single relaxation time is a strong assumption for a crosslinked biopolymer, and the
one-mode residual on the gelatin records is correlated at lag one rather than white. Two
laws enrich it in different directions, and both reduce to QuadraticZener exactly rather
than approximately -- which is what their tests assert, since a reduction that merely held
closely would leave room for an error of its own size.
TwoModeQuadraticZener adds a second Maxwell arm: two timescales, with second_share
splitting both the elastic target and the viscous forcing between them. At
second_share=0 the second memory is driven by nothing, decays from zero, and the law is
QuadraticZener.
material = TwoModeQuadraticZener(
shear_modulus_pa=204.3,
viscosity_pa_s=0.047,
relaxation_time_s=2.0e-7,
retardation_time_s=0.0,
stiffening=5.3,
second_relaxation_time_s=2.0e-6,
second_share=0.3,
)CarreauZener keeps one arm and lets its timescale move. A Maxwell arm is a spring and a
dashpot in series, so lambda = eta/G; if the dashpot obeys Carreau rather than Newton,
the relaxation time falls as the medium is sheared harder, and a collapse spans decades of
shear rate. Only the dashpot thins -- the elastic target and the solvent term are the
one-mode law's. At power_index=1 the thinning factor is exactly one at every shear rate,
so the reduction is bit-identical.
material = CarreauZener(
shear_modulus_pa=204.3,
viscosity_pa_s=0.047,
relaxation_time_s=2.0e-7,
retardation_time_s=0.0,
stiffening=5.3,
thinning_time_s=2.0e-6,
power_index=0.6,
)Both are comparison candidates, in EXTENDED_MODELS rather than STANDARD_MODELS: at six
free parameters the grid quadrature in pyimr.selection costs count**6, so they are
scored by candidate_log_evidence instead. See model selection.
Giesekus and LinearPTT evolve radial and hoop stress on a prepared,
wall-clustered Lagrangian grid:
material = Giesekus(
viscosity_pa_s=0.1,
relaxation_time_s=40e-6,
retardation_time_s=8e-6,
mobility=0.2,
)result.stress_state contains radial stress followed by hoop stress on
result.stress_reference_radius_ratio.
The constitutive equations at each material point are ordinary differential
equations -- there are no spatial derivatives -- so the only spatial
approximation is the quadrature for the stress integral. quadrature="gauss"
(the default, 240 points) places the material points at Gauss-Legendre nodes and
converges spectrally: about five orders of magnitude more accurate than the
trapezoid grid it replaced, and cheaper. The table is in
discretization.md.
At zero mobility or zero extensibility these converge to the analytic Oldroyd-B
solution. Use OldroydB directly when that closure applies -- it is
substantially smaller and faster.
Memory materials can initialize from a resolved equilibrium-to-maximum-radius precursor rather than an assumed unstressed state:
from pyimr import CollapseInitialization
config = SimulationConfig(
R0=225e-6, Req=37.5e-6,
material=Zener(2500.0, 0.1, 40e-6, 8e-6),
collapse=CollapseInitialization(),
)
problem = prepare(config)
problem.collapse_statsPreparation brackets the precursor velocity, shoots to R/R0 == 1, and retains
the complete memory state at the maximum. CollapseStats records the root,
achieved maximum, integration work, and immutable stress state. Sensitivities
differentiate the event time and shooting root implicitly. Oldroyd-B and
distributed Giesekus/PTT states use the same mechanism; the Zener precursor
retains the upstream IMRv2 formulation.