Linear viscoelastic models

hermes_rheo.lve_models simulates SAOS, creep and stress relaxation for classical and fractional linear viscoelastic (LVE) models. Each experiment returns a piblin Measurement with the same column names and units as real data, so it works with every hermes-rheo transform.

Each model defines its complex modulus \(G^*(s)\):

• SAOS: \(G^*(i\omega) = G' + iG''\)

• Stress relaxation: \(G(t) = \mathcal{L}^{-1}[G^*(s)/s]\)

• Creep: \(J(t) = \mathcal{L}^{-1}[1/(s\,G^*(s))]\)

Closed forms are used where they exist. Otherwise the Laplace inversion is numerical (fixed Talbot). Models with a parallel dashpot (Kelvin–Voigt, Jeffreys) have a \(\delta(t)\) term in \(G(t)\); only the \(t>0\) part is returned.

Models

Class	Parameters
MaxwellModel	G, tau
MultiModeMaxwellModel	g, tau, G_inf
KelvinVoigtModel	G, eta
JeffreysModel	G, tau, eta_s
SpringPotModel	V, alpha
FractionalMaxwellGelModel	G, V, alpha
FractionalMaxwellLiquidModel	G, eta, beta
FractionalMaxwellModel	G, V, alpha, beta (alpha > beta)
FractionalKelvinVoigtModel	G, V, alpha, beta
FractionalZenerModel	G, V, K, alpha, beta, kappa

Fractional orders lie in the open interval (0, 1). Parameter names follow pyRheo.

Example

import numpy as np
from hermes_rheo.lve_models import FractionalMaxwellGelModel, Protocols

model = FractionalMaxwellGelModel(G=1e4, V=5e3, alpha=0.6)
saos = Protocols.saos(model, np.logspace(-2, 2, 40))
creep = Protocols.creep(model, np.logspace(-2, 2, 40))
relax = Protocols.stress_relaxation(model, np.logspace(-2, 2, 40))

Metadata

conditions["method"] names the protocol and model, and conditions["simulated"] is True. details holds model, parameters, protocol, grid (variable, min, max, points), references and software.

Protocols.experiment_set(model, [...]) builds several experiments in one call.

Validation

The tests check closed forms against the inverted \(G^*\), the identity \(\int_0^t G(t-u)J(u)\,du = t\), the Mittag-Leffler series, and the Maxwell limit of the fractional gel as \(\alpha \to 1\). pyRheo is used only as an optional cross-check; it is GPL-3.0 and is not a dependency.

References

• Bonfanti, A., Kaplan, J. L., Charras, G. and Kabla, A., Fractional viscoelastic models for power-law materials, Soft Matter **16**(26), 6002–6020 (2020). DOI: 10.1039/d0sm00354a.

• Schiessel, H., Metzler, R., Blumen, A., Nonnenmacher, T. F. “Generalized viscoelastic models: their fractional equations with solutions.” J. Phys. A 28, 6567 (1995).

• Garrappa, R. “Numerical evaluation of two and three parameter Mittag-Leffler functions.” SIAM J. Numer. Anal. 53, 1350 (2015).

• Talbot, A. “The accurate numerical inversion of Laplace transforms.” J. Inst. Math. Appl. 23, 97 (1979).

• Song, J., Holten-Andersen, N. and McKinley, G. H., Non-Maxwellian viscoelastic stress relaxations in soft matter, Soft Matter **19**(41), 7885–7906 (2023). DOI: 10.1039/d3sm00736g.

• Miranda-Valdez, I. Y., Niinisto, H., Makinen, T. P., Lejon, E., Koivisto, J. and Alava, M. J., pyRheo: an open-source Python package for complex rheology, Digital Discovery **4**(4), 1075–1082 (2025). DOI: 10.1039/d5dd00021a.

• pyRheo source: https://github.com/3mcloud/pyRheo

API

Simulated linear viscoelastic experiments.

Run SAOS (small-amplitude oscillatory shear), creep and stress relaxation on classical and fractional models. Each experiment returns a piblin Measurement with one composite one-dimensional dataset and metadata describing the model, its parameters and the simulated protocol, so it passes through the same hermes-rheo transforms as measured data.

Models

Classical: MaxwellModel, MultiModeMaxwellModel, KelvinVoigtModel, JeffreysModel.

Fractional: SpringPotModel, FractionalMaxwellGelModel, FractionalMaxwellLiquidModel, FractionalMaxwellModel, FractionalKelvinVoigtModel, FractionalZenerModel.

Examples

>>> import numpy as np
>>> from hermes_rheo.lve_models import MultiModeMaxwellModel, Protocols
>>> model = MultiModeMaxwellModel(g=[1e5, 1e4], tau=[0.1, 1.0])
>>> meas = Protocols.stress_relaxation(model, np.geomspace(1e-3, 1e2, 200))
>>> meas.datasets[0].data_array_names
['time', 'relaxation modulus']
>>> meas.conditions["simulated"], meas.details["protocol"]
(True, 'stress_relaxation')

class hermes_rheo.lve_models.FractionalKelvinVoigtModel(G: float, V: float, alpha: float, beta: float)

Springpot \((G, \beta)\) in parallel with springpot \((V, \alpha)\).

\[G^*(s) = G s^{\beta} + V s^{\alpha}, \qquad G(t) = \frac{G t^{-\beta}}{\Gamma(1-\beta)} + \frac{V t^{-\alpha}}{\Gamma(1-\alpha)}\]

Parameters:

• G (float) – Quasi-modulus in Pa s^beta.

• V (float) – Quasi-modulus in Pa s^alpha.

• alpha (float) – Orders in (0, 1).

• beta (float) – Orders in (0, 1).

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

stress_relaxation(t)

Return the relaxation modulus \(G(t)\) in Pa for \(t > 0\).

Any impulsive term at \(t = 0\) (a parallel dashpot) is excluded.

class hermes_rheo.lve_models.FractionalMaxwellGelModel(G: float, V: float, alpha: float)

Spring \(G\) in series with a springpot \((V, \alpha)\).

\[G^*(s) = \frac{G V s^{\alpha}}{G + V s^{\alpha}}\]

The relaxation modulus is \(G\,E_{\alpha}(-(t/\tau_c)^{\alpha})\) with \(\tau_c = (V/G)^{1/\alpha}\) (Mittag-Leffler function).

Parameters:

• G (float) – Modulus in Pa.

• V (float) – Quasi-modulus in Pa s^alpha.

• alpha (float) – Order, in (0, 1).

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

class hermes_rheo.lve_models.FractionalMaxwellLiquidModel(G: float, eta: float, beta: float)

Springpot \((G, \beta)\) in series with a dashpot \(\eta\).

\[G^*(s) = \frac{G s^{\beta}\,\eta s}{G s^{\beta} + \eta s}\]

Parameters:

• G (float) – Quasi-modulus in Pa s^beta.

• eta (float) – Viscosity in Pa s.

• beta (float) – Order, in (0, 1).

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

class hermes_rheo.lve_models.FractionalMaxwellModel(G: float, V: float, alpha: float, beta: float)

Springpot \((G, \beta)\) in series with springpot \((V, \alpha)\).

\[G^*(s) = \frac{G s^{\beta}\,V s^{\alpha}}{G s^{\beta} + V s^{\alpha}}, \qquad \alpha > \beta\]

Parameters:

• G (float) – Quasi-modulus in Pa s^beta.

• V (float) – Quasi-modulus in Pa s^alpha.

• alpha (float) – Orders in (0, 1) with alpha > beta.

• beta (float) – Orders in (0, 1) with alpha > beta.

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

class hermes_rheo.lve_models.FractionalZenerModel(G: float, V: float, K: float, alpha: float, beta: float, kappa: float)

Springpot \((K, \kappa)\) in parallel with a fractional Maxwell model.

\[G^*(s) = K s^{\kappa} + \frac{G s^{\beta}\,V s^{\alpha}}{G s^{\beta} + V s^{\alpha}}, \qquad \alpha > \beta\]

Parameters:

• G (float) – Quasi-modulus in Pa s^beta.

• V (float) – Quasi-modulus in Pa s^alpha.

• K (float) – Quasi-modulus in Pa s^kappa.

• alpha (float) – Orders in (0, 1) with alpha > beta.

• beta (float) – Orders in (0, 1) with alpha > beta.

• kappa (float) – Orders in (0, 1) with alpha > beta.

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

class hermes_rheo.lve_models.JeffreysModel(G: float, tau: float, eta_s: float)

Jeffreys model: a Maxwell mode in parallel with a dashpot.

With \(\eta_p = G\tau\) and \(\eta = \eta_p + \eta_s\):

\[G^*(s) = \eta_s s + \frac{G\tau s}{1 + \tau s}, \qquad J(t) = \frac{t}{\eta} + \frac{\tau\,\eta_p}{\eta^2}\left(1 - e^{-t/\tau_r}\right), \quad \tau_r = \frac{\eta_s \tau}{\eta}\]

The relaxation modulus is \(G e^{-t/\tau} + \eta_s\,\delta(t)\). For \(t > 0\) stress_relaxation() returns the first term.

Parameters:

• G (float) – Modulus of the Maxwell mode in Pa.

• tau (float) – Relaxation time of the Maxwell mode in s.

• eta_s (float) – Solvent viscosity in Pa s.

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

creep(t)

Return the creep compliance \(J(t)\) in 1/Pa for \(t > 0\).

stress_relaxation(t)

Return the relaxation modulus \(G(t)\) in Pa for \(t > 0\).

Any impulsive term at \(t = 0\) (a parallel dashpot) is excluded.

class hermes_rheo.lve_models.KelvinVoigtModel(G: float, eta: float)

Kelvin-Voigt solid: spring \(G\) and dashpot \(\eta\) in parallel.

\[J(t) = \frac{1}{G}\left(1 - e^{-t/\tau_r}\right), \quad \tau_r = \frac{\eta}{G}, \qquad G^*(s) = G + \eta s\]

The relaxation modulus is \(G + \eta\,\delta(t)\). For \(t > 0\) stress_relaxation() returns the constant \(G\).

Parameters:

• G (float) – Modulus in Pa.

• eta (float) – Viscosity in Pa s.

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

creep(t)

Return the creep compliance \(J(t)\) in 1/Pa for \(t > 0\).

stress_relaxation(t)

Return the relaxation modulus \(G(t)\) in Pa for \(t > 0\).

Any impulsive term at \(t = 0\) (a parallel dashpot) is excluded.

class hermes_rheo.lve_models.LVEModel

Abstract linear viscoelastic model.

Subclasses implement complex_modulus(). They may override relaxation_modulus() and creep_compliance() with exact expressions.

Variables:

• name (str) – Short model name, stored in the measurement metadata.

• parameters (dict) – Model parameters, stored in the measurement metadata.

• references (tuple of str) – Literature for the model, stored in the measurement metadata.

abstractmethod complex_modulus(s: ndarray) → ndarray

Return \(G^*(s)\) in Pa for complex Laplace variable s.

creep(t: ndarray) → ndarray

Return the creep compliance \(J(t)\) in 1/Pa for \(t > 0\).

saos(omega: ndarray) → tuple[ndarray, ndarray]

Return the storage and loss moduli \((G', G'')\) in Pa.

Parameters:

omega (np.ndarray) – Angular frequency in rad/s (positive).

stress_relaxation(t: ndarray) → ndarray

Return the relaxation modulus \(G(t)\) in Pa for \(t > 0\).

Any impulsive term at \(t = 0\) (a parallel dashpot) is excluded.

to_measurement(protocol: str, x: ndarray, *, conditions: dict[str, Any] | None = None, details: dict[str, Any] | None = None) → Measurement

Run a simulated experiment and return a piblin Measurement.

Parameters:

• protocol ({"saos", "creep", "stress_relaxation"}) – Experiment to simulate.

• x (np.ndarray) – Angular frequency in rad/s for "saos", time in s otherwise.

• conditions (dict, optional) – Extra entries merged into Measurement.conditions / details.

• details (dict, optional) – Extra entries merged into Measurement.conditions / details.

Returns:

One composite one-dimensional dataset. conditions hold method and simulated. details hold the model name, parameters, protocol, grid, references and software version.

Return type:

piblin.data.data_collections.measurement.Measurement

class hermes_rheo.lve_models.MaxwellModel(G: float, tau: float)

Maxwell fluid: spring \(G\) and dashpot \(\eta = G\tau\) in series.

\[G(t) = G e^{-t/\tau}, \qquad J(t) = \frac{1}{G}\left(1 + \frac{t}{\tau}\right), \qquad G^*(s) = \frac{G \tau s}{1 + \tau s}\]

Parameters:

• G (float) – Modulus in Pa.

• tau (float) – Relaxation time in s.

Examples

>>> m = MaxwellModel(G=1e5, tau=1.0)
>>> gp, gpp = m.saos(np.array([1.0]))
>>> bool(np.isclose(gp[0], 5e4) and np.isclose(gpp[0], 5e4))
True

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

creep(t)

Return the creep compliance \(J(t)\) in 1/Pa for \(t > 0\).

stress_relaxation(t)

Return the relaxation modulus \(G(t)\) in Pa for \(t > 0\).

Any impulsive term at \(t = 0\) (a parallel dashpot) is excluded.

class hermes_rheo.lve_models.MultiModeMaxwellModel(g: Sequence[float], tau: Sequence[float], G_inf: float = 0.0)

Parallel Maxwell modes with an optional equilibrium modulus.

\[G(t) = G_\infty + \sum_i g_i e^{-t/\tau_i}, \qquad G^*(s) = G_\infty + \sum_i \frac{g_i \tau_i s}{1 + \tau_i s}\]

The creep compliance has no simple closed form and is computed by numerical Laplace inversion (\(t > 0\); at \(t = 0\) it is \(J(0) = 1/G(0)\)).

Parameters:

• g (sequence of float) – Mode strengths in Pa.

• tau (sequence of float) – Relaxation times in s, same length as g.

• G_inf (float, optional) – Equilibrium modulus in Pa. 0 (default) gives a fluid.

Examples

>>> m = MultiModeMaxwellModel(g=[1e5, 1e4], tau=[0.1, 1.0])
>>> bool(np.isclose(m.stress_relaxation(np.array([0.0]))[0], 1.1e5))
True

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

creep(t)

Return the creep compliance \(J(t)\) in 1/Pa for \(t > 0\).

stress_relaxation(t)

Return the relaxation modulus \(G(t)\) in Pa for \(t > 0\).

Any impulsive term at \(t = 0\) (a parallel dashpot) is excluded.

class hermes_rheo.lve_models.Protocols

Simulated experiments on an LVE model. All methods are static.

static creep(model: LVEModel, t: ndarray, **kwargs: Any) → Measurement

Creep compliance. t in s.

static experiment_set(model: LVEModel, protocols: Sequence[dict[str, Any]]) → ExperimentSet | MeasurementSet

Run several protocols on one model.

Parameters:

• model (LVEModel)

• protocols (sequence of dict) – Each dict has "protocol" ("saos", "creep" or "stress_relaxation") and "x" (the grid), plus optional "conditions" and "details".

Returns:

An ExperimentSet when piblin can group the measurements, otherwise a MeasurementSet.

Return type:

ExperimentSet or MeasurementSet

static saos(model: LVEModel, omega: ndarray, **kwargs: Any) → Measurement

Small-amplitude oscillatory shear. omega in rad/s.

static stress_relaxation(model: LVEModel, t: ndarray, **kwargs: Any) → Measurement

Stress relaxation modulus. t in s.

class hermes_rheo.lve_models.SpringPotModel(V: float, alpha: float)

Springpot: power-law element.

\[G(t) = \frac{V t^{-\alpha}}{\Gamma(1-\alpha)}, \qquad J(t) = \frac{t^{\alpha}}{V\,\Gamma(1+\alpha)}, \qquad G^*(s) = V s^{\alpha}\]

Parameters:

• V (float) – Quasi-modulus in Pa s^alpha.

• alpha (float) – Order, in (0, 1).

complex_modulus(s)

Return \(G^*(s)\) in Pa for complex Laplace variable s.

creep(t)

Return the creep compliance \(J(t)\) in 1/Pa for \(t > 0\).

stress_relaxation(t)

Return the relaxation modulus \(G(t)\) in Pa for \(t > 0\).

Any impulsive term at \(t = 0\) (a parallel dashpot) is excluded.

hermes_rheo.lve_models.invert_laplace(F: Callable[[ndarray], ndarray], t: ndarray, n_terms: int = 16) → ndarray

Numerically invert a Laplace transform with the fixed Talbot method.

Parameters:

• F (callable) – Transform \(F(s)\). Must accept and return complex arrays.

• t (np.ndarray) – Times in s. All values must be positive.

• n_terms (int, optional) – Number of contour points. 16 gives about 1e-9 relative accuracy in double precision; larger values lose accuracy to round-off.

Returns:

\(f(t)\), same shape as t.

Return type:

np.ndarray