Viscoelasticity

Overview

Viscoelastic models describe time-dependent stress-strain behavior that is partly elastic (recoverable) and partly viscous (rate-dependent). A viscoelastic material carries memory of past loading through internal strain variables, but the constitutive response stays linear — there is no yield surface and no consistency parameter.

NEML2 builds viscoelastic models from two primitive rheological elements — a linear spring with modulus \(E\) and a Newtonian dashpot with viscosity \(\eta\) — combined in series and parallel. Standard textbook configurations are shipped as pre-assembled element models for convenience (KelvinVoigtElement, ZenerElement, WiechertElement, BurgersElement), while LinearDashpot is exposed as a leaf for assembling arbitrary networks from primitives. Each model treats the spring and dashpot relations component-wise on \(\boldsymbol{\sigma}\) and \(\boldsymbol{\varepsilon}\), so the same parameters describe both the deviatoric and volumetric response — a common simplification when only one effective modulus is available from experiment.

Governing equations

The two primitives are

\[ \boldsymbol{\sigma} = E \boldsymbol{\varepsilon} \quad\text{(spring)}, \qquad \boldsymbol{\sigma} = \eta \dot{\boldsymbol{\varepsilon}} \quad\text{(dashpot)}. \]

A series connection adds strains under a shared stress; a parallel connection adds stresses under a shared strain. The four named topologies in the catalog follow from these two rules. Taking \(\boldsymbol{\varepsilon}_v\) as the dashpot (viscous) strain in each branch:

For KelvinVoigtElement — spring and dashpot in parallel — the total stress is the sum of the two branches and there is no separate viscous-strain unknown:

\[ \boldsymbol{\sigma} = E \boldsymbol{\varepsilon} + \eta \dot{\boldsymbol{\varepsilon}}. \]

For ZenerElement — equilibrium spring in parallel with a Maxwell branch (Standard Linear Solid):

(13)\[\begin{align} \boldsymbol{\sigma} &= (E_\infty + E_M) \boldsymbol{\varepsilon} - E_M \boldsymbol{\varepsilon}_v, \\ \dot{\boldsymbol{\varepsilon}}_v &= \frac{E_M}{\eta_M} (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_v). \end{align}\]

For WiechertElement — an equilibrium spring in parallel with two Maxwell branches (generalized Maxwell), one viscous strain \(\boldsymbol{\varepsilon}_{v,i}\) per branch:

(14)\[\begin{align} \boldsymbol{\sigma} &= E_\infty \boldsymbol{\varepsilon} + \sum_{i=1}^{2} E_i \left( \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_{v,i} \right), \\ \dot{\boldsymbol{\varepsilon}}_{v,i} &= \frac{E_i}{\eta_i} \left( \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_{v,i} \right), \quad i = 1, 2. \end{align}\]

For BurgersElement — a Maxwell branch in series with a Kelvin-Voigt block, with \(\boldsymbol{\varepsilon}_v^M\) the Maxwell dashpot strain and \(\boldsymbol{\varepsilon}^K\) the Kelvin-Voigt block strain:

(15)\[\begin{align} \boldsymbol{\sigma} &= E_M \left( \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_v^M - \boldsymbol{\varepsilon}^K \right), \\ \dot{\boldsymbol{\varepsilon}}_v^M &= \boldsymbol{\sigma} / \eta_M, \\ \dot{\boldsymbol{\varepsilon}}^K &= \left( \boldsymbol{\sigma} - E_K \boldsymbol{\varepsilon}^K \right) / \eta_K. \end{align}\]

In every case the constitutive model exposes the rate \(\dot{\boldsymbol{\varepsilon}}_v\) of each internal strain variable as a derived output. SR2BackwardEulerTimeIntegration turns each rate into an implicit residual, and an ImplicitUpdate wrapping a Newton solver advances the unknowns one time step at a time.

Example: Zener element (Standard Linear Solid)

The simplest non-trivial case is the Zener element. The single viscous strain \(\boldsymbol{\varepsilon}_v\) is the only Newton unknown, and the canned ZenerElement closure provides both the stress and the viscous-strain rate.

# neml2
# Stress relaxation of a Zener (Standard Linear Solid) viscoelastic element under a ramped strain
# loading. With a single Maxwell branch in parallel with an equilibrium spring, the stress
# overshoots during loading and relaxes to the equilibrium value E_inf * strain when the strain is
# held fixed.
[Tensors]
  [end_time]
    type = Python
    expr = 'Scalar(torch.logspace(-1, 3, 10, dtype=torch.float64))'
  []
  [times]
    type = Python
    expr = 'Scalar(torch.stack([torch.linspace(0, t.item(), 100, dtype=torch.float64) for t in end_time.data], dim=-1))'
  []
  [strains]
    # Mandel-fill (10, 6) max_strain = (exx, eyy, ezz, 0, 0, 0) = (0.01, -0.005, -0.005, 0, 0, 0),
    # then linspace ramp from 0 to max across 100 timesteps -> shape (100, 10, 6).
    type = Python
    expr = 'SR2(torch.tensor([0.01, -0.005, -0.005, 0.0, 0.0, 0.0], dtype=torch.float64).reshape(1, 1, 6).expand(100, 10, 6) * torch.linspace(0.0, 1.0, 100, dtype=torch.float64).reshape(100, 1, 1))'
  []
[]

[Drivers]
  [driver]
    type = TransientDriver
    model = 'model'
    prescribed_time = 'times'
    force_SR2_names = 'strain'
    force_SR2_values = 'strains'
    save_as = 'result.pt'
  []
  [regression]
    type = TransientRegression
    driver = 'driver'
    reference = 'gold/result.pt'
  []
[]

[Models]
  [zener]
    type = ZenerElement
    equilibrium_modulus = 1000
    maxwell_modulus = 5000
    maxwell_viscosity = 100
  []
  [integrate_Ev]
    type = SR2BackwardEulerTimeIntegration
    variable = 'viscous_strain'
  []
  [implicit_rate]
    type = ComposedModel
    models = 'zener integrate_Ev'
  []
[]

[EquationSystems]
  [eq_sys]
    type = NonlinearSystem
    model = 'implicit_rate'
    unknowns = 'viscous_strain'
    residuals = 'viscous_strain_residual'
  []
[]

[Solvers]
  [newton]
    type = Newton
    linear_solver = 'lu'
  []
  [lu]
    type = DenseLU
  []
[]

[Models]
  [predictor]
    type = ConstantExtrapolationPredictor
    unknowns_SR2 = 'viscous_strain'
  []
  [update]
    type = ImplicitUpdate
    equation_system = 'eq_sys'
    solver = 'newton'
    predictor = 'predictor'
  []
  [model]
    type = ComposedModel
    models = 'update zener'
    additional_outputs = 'viscous_strain'
  []
[]

Walking through the wiring:

  1. [zener] is the constitutive closure. It consumes strain (the forcing) and viscous_strain (the unknown) and produces stress and viscous_strain_rate. The three parameters — equilibrium_modulus (\(E_\infty\)), maxwell_modulus (\(E_M\)), and maxwell_viscosity (\(\eta_M\)) — are the only material constants in the model.

  2. [integrate_Ev] is a SR2BackwardEulerTimeIntegration instance. It consumes viscous_strain (current iterate), viscous_strain_rate (from zener), and the old-state pair, and emits a residual viscous_strain_residual that is zero when the backward-Euler identity \(\boldsymbol{\varepsilon}_v - \boldsymbol{\varepsilon}_v^{n-1} = \Delta t \dot{\boldsymbol{\varepsilon}}_v\) holds.

  3. [implicit_rate] is a ComposedModel that glues zener and integrate_Ev together by variable name — there is no explicit edge list, just shared input/output names.

  4. [eq_sys] declares the NonlinearSystem: the unknown is viscous_strain, the residual is viscous_strain_residual, and the model that computes the residual is implicit_rate.

  5. [update] wraps the system in an ImplicitUpdate. On each forward pass, the ConstantExtrapolationPredictor seeds viscous_strain from the previous step, the Newton solver drives the residual to zero, and the converged viscous_strain is fed forward.

  6. [model] composes the implicit update with a fresh evaluation of zener so that the post-solve stress is computed from the converged internal state. The additional_outputs = 'viscous_strain' line keeps the internal variable visible to the driver for time integration of the next step.

Swapping in any other pre-assembled element is a one-line change to the [zener] block — the rest of the wiring (time integration, equation system, implicit update) is identical for any single-unknown viscoelastic model. The WiechertElement and BurgersElement classes extend this pattern to multi-branch generalized-Maxwell and Burgers networks without requiring you to assemble springs and dashpots yourself.

Composing custom networks

When a topology doesn’t match a named element — for example a 5-element generalized Maxwell or a Burgers-plus-extra-dashpot — the network can be assembled from primitives directly. The recipe is purely topological:

  • Use LinearDashpot for each dashpot leaf. It maps \(\boldsymbol{\sigma} \mapsto \dot{\boldsymbol{\varepsilon}}_v = \boldsymbol{\sigma} / \eta\).

  • Use SR2LinearCombination (or a full LinearIsotropicElasticity) as the constitutive law for each spring, and again for the strain decompositions and stress sums that encode series/parallel topology.

  • Wire it all up with ComposedModel. The graph is by variable name: series-chain elements share a stress variable, parallel branches share a strain variable, and each dashpot’s strain is one unknown in the NonlinearSystem.

Worked examples for every named topology, plus a 5-element generalized Maxwell, live under tests/regression/solid_mechanics/viscoelasticity/composed_*/. The Burgers composition is the most instructive — two viscous strains share the chain stress, so it exercises both series and parallel assembly in one input file:

# neml2
# Burgers response built from primitives — Maxwell branch in series with a Kelvin-Voigt block,
# both internal-state dashpots solved together. Same load history and parameters as
# `burgers/model.i`. The trick for any series chain: a single shared stress flows through every
# element, so the Maxwell spring's `stress` output is consumed both by the Maxwell dashpot AND
# (after subtracting the Kelvin spring's stress) by the Kelvin dashpot. Two viscous strains, two
# residuals, one Newton solve.
#
# Math, cross-checked against the closed-form BurgersElement:
#   stress = E_M * (strain - mvs - kvs)
#   maxwell_viscous_strain_rate = stress / eta_M
#   kelvin_dashpot_stress = stress - E_K * kvs
#   kelvin_voigt_strain_rate = kelvin_dashpot_stress / eta_K
[Tensors]
  [end_time]
    type = Python
    expr = 'Scalar(torch.logspace(-1, 3, 10, dtype=torch.float64))'
  []
  [times]
    type = Python
    expr = 'Scalar(torch.stack([torch.linspace(0, t.item(), 100, dtype=torch.float64) for t in end_time.data], dim=-1))'
  []
  [strains]
    # Mandel-fill (10, 6) max_strain = (exx, eyy, ezz, 0, 0, 0) = (0.01, -0.005, -0.005, 0, 0, 0),
    # then linspace ramp from 0 to max across 100 timesteps -> shape (100, 10, 6).
    type = Python
    expr = 'SR2(torch.tensor([0.01, -0.005, -0.005, 0.0, 0.0, 0.0], dtype=torch.float64).reshape(1, 1, 6).expand(100, 10, 6) * torch.linspace(0.0, 1.0, 100, dtype=torch.float64).reshape(100, 1, 1))'
  []
[]

[Drivers]
  [driver]
    type = TransientDriver
    model = 'model'
    prescribed_time = 'times'
    force_SR2_names = 'strain'
    force_SR2_values = 'strains'
    save_as = 'result.pt'
  []
  [regression]
    type = TransientRegression
    driver = 'driver'
    reference = 'gold/result.pt'
  []
[]

[Models]
  # Maxwell-branch elastic strain: strain absorbed by the series chain that isn't already in the
  # Maxwell dashpot or in the Kelvin-Voigt block.
  [elastic_strain]
    type = SR2LinearCombination
    from = 'strain maxwell_viscous_strain kelvin_voigt_strain'
    weights = '1 -1 -1'
    to = 'elastic_strain'
  []
  [maxwell_spring]
    type = SR2LinearCombination
    from = 'elastic_strain'
    weights = '5000'
    to = 'stress'
  []
  [maxwell_dashpot]
    type = LinearDashpot
    stress = 'stress'
    viscous_strain_rate = 'maxwell_viscous_strain_rate'
    viscosity = 200
  []

  # Kelvin-Voigt block (parallel of spring and dashpot, both under stress = chain stress)
  [kelvin_spring_stress]
    type = SR2LinearCombination
    from = 'kelvin_voigt_strain'
    weights = '2000'
    to = 'kelvin_spring_stress'
  []
  [kelvin_dashpot_stress]
    type = SR2LinearCombination
    from = 'stress kelvin_spring_stress'
    weights = '1 -1'
    to = 'kelvin_dashpot_stress'
  []
  [kelvin_dashpot]
    type = LinearDashpot
    stress = 'kelvin_dashpot_stress'
    viscous_strain_rate = 'kelvin_voigt_strain_rate'
    viscosity = 100
  []

  [integrate_EvM]
    type = SR2BackwardEulerTimeIntegration
    variable = 'maxwell_viscous_strain'
  []
  [integrate_EK]
    type = SR2BackwardEulerTimeIntegration
    variable = 'kelvin_voigt_strain'
  []
  [implicit_rate]
    type = ComposedModel
    models = 'elastic_strain maxwell_spring maxwell_dashpot
              kelvin_spring_stress kelvin_dashpot_stress kelvin_dashpot
              integrate_EvM integrate_EK'
  []
[]

[EquationSystems]
  [eq_sys]
    type = NonlinearSystem
    model = 'implicit_rate'
    unknowns = 'maxwell_viscous_strain kelvin_voigt_strain'
    residuals = 'maxwell_viscous_strain_residual kelvin_voigt_strain_residual'
  []
[]

[Solvers]
  [newton]
    type = Newton
    linear_solver = 'lu'
  []
  [lu]
    type = DenseLU
  []
[]

[Models]
  [predictor]
    type = ConstantExtrapolationPredictor
    unknowns_SR2 = 'maxwell_viscous_strain kelvin_voigt_strain'
  []
  [update]
    type = ImplicitUpdate
    equation_system = 'eq_sys'
    solver = 'newton'
    predictor = 'predictor'
  []
  [model]
    type = ComposedModel
    models = 'update elastic_strain maxwell_spring maxwell_dashpot
              kelvin_spring_stress kelvin_dashpot_stress kelvin_dashpot'
    additional_outputs = 'maxwell_viscous_strain kelvin_voigt_strain stress'
  []
[]

The unknowns are now maxwell_viscous_strain and kelvin_voigt_strain — two SR2BackwardEulerTimeIntegration instances, two residuals, one Newton solve. The hand-composed input is mathematically identical to tests/regression/.../burgers/model.i (which uses BurgersElement directly) and produces the same gold result — the difference is that the composed form lets you change the topology by editing the input file, without authoring a new Model class.

See also