TNT Tanaka-Edwards (Basic Transient Network) — Handbook

Quick Reference

Use when:

  • Simple associating polymers with reversible physical crosslinks

  • Telechelic networks (e.g., PEG-PEO telechelics)

  • Dilute wormlike micelles (initial approximation)

  • Initial model for any transient network system

  • Maxwell-like behavior expected (single relaxation time)

Parameters:

3 parameters: \(G\) (network modulus, Pa), \(\tau_b\) (bond lifetime, s), \(\eta_s\) (solvent viscosity, Pa·s)

Key equation:
\[\frac{d\mathbf{S}}{dt} = \boldsymbol{\kappa} \cdot \mathbf{S} + \mathbf{S} \cdot \boldsymbol{\kappa}^T - \frac{\mathbf{S} - \mathbf{I}}{\tau_b}\]
Test modes:

All 6 protocols (FLOW_CURVE, OSCILLATION, RELAXATION, STARTUP, CREEP, LAOS)

Material examples:

  • Hydrophobically modified ethoxylated urethane (HEUR) solutions

  • PVA-borax gels

  • Telechelic PEG-PEO associative polymers

  • Reversible gelatin networks

  • Casein micelle dispersions (dilute)

Key characteristics:

  • Newtonian steady-state flow (no shear thinning without variants)

  • Quadratic normal stress difference: \(N_1 \propto \dot{\gamma}^2\)

  • Single exponential stress relaxation

  • Stress overshoot in startup flow

  • Maxwell-type SAOS response

Notation Guide

Symbol

Units

Description

\(\mathbf{S}\)

dimensionless

Conformation tensor (3×3 symmetric, normalized second moment of end-to-end vector)

\(G\)

Pa

Network elastic modulus (related to chain density via \(G \approx n_{chains} k_B T\))

\(\tau_b\)

s

Bond lifetime — mean time between bond breaking and reformation events

\(\eta_s\)

Pa·s

Solvent viscosity (non-network contribution to total viscosity)

\(\eta_0\)

Pa·s

Zero-shear viscosity: \(\eta_0 = G \tau_b + \eta_s\)

\(\boldsymbol{\kappa}\)

1/s

Velocity gradient tensor: \(\kappa_{ij} = \partial v_i / \partial x_j\)

\(\mathbf{D}\)

1/s

Rate of deformation tensor: \(\mathbf{D} = (\boldsymbol{\kappa} + \boldsymbol{\kappa}^T)/2\)

\(\dot{\gamma}\)

1/s

Shear rate (scalar, for simple shear flow)

\(Wi\)

dimensionless

Weissenberg number: \(Wi = \tau_b \dot{\gamma}\) (ratio of relaxation to flow timescales)

\(\mathbf{I}\)

dimensionless

Identity tensor (equilibrium conformation)

\(\mathbf{R}\)

m

End-to-end vector of network chain between crosslinks

\(R_0\)

m

Equilibrium root-mean-square end-to-end distance

\(N_1\)

Pa

First normal stress difference: \(N_1 = \sigma_{xx} - \sigma_{yy}\)

\(N_2\)

Pa

Second normal stress difference: \(N_2 = \sigma_{yy} - \sigma_{zz}\)

Overview

The Tanaka-Edwards model is the foundational transient network theory for viscoelastic materials where polymer chains are connected by reversible physical crosslinks. Originally proposed by Green and Tobolsky (1946) in their seminal work on elastomers with interchanging bonds, and later formalized by Tanaka and Edwards (1992) using a conformation tensor approach, this model describes materials where bonds break and reform with a constant rate \(1/\tau_b\), independent of chain stretch or applied force.

At equilibrium, bonds break and reform at equal rates, maintaining a constant network structure. Under deformation, chains between crosslinks stretch and orient, generating stress proportional to the deviation of the conformation tensor \(\mathbf{S}\) from its equilibrium value (the identity tensor \(\mathbf{I}\)). When a bond breaks, the chain relaxes instantly to its equilibrium conformation; when a new bond forms, the chain joins the network with the current equilibrium state. This kinetic balance gives rise to a single relaxation mode with Maxwell-like linear viscoelastic behavior.

The model predicts:

  • Linear viscoelasticity: Single-mode Maxwell response in SAOS with \(G' \sim \omega^2\) and \(G'' \sim \omega\) at low frequencies

  • Steady-state flow: Newtonian behavior (constant viscosity \(\eta_0\))

  • Normal stresses: Quadratic in shear rate (\(N_1 = 2G(\tau_b \dot{\gamma})^2\))

  • Transient response: Stress overshoot in startup flow, single exponential relaxation

  • LAOS: Nonlinear harmonics from conformation tensor evolution

Despite its simplicity (only 3 parameters), the Tanaka-Edwards model provides the molecular foundation for understanding more complex transient network materials and serves as a starting point for extensions incorporating force-dependent breakage, finite extensibility, and multiple relaxation modes.

Historical Context

The development of transient network theory spans several decades:

  1. Green & Tobolsky (1946): First transient network model for rubber-like materials with reversible crosslinks. They introduced the concept of a network with bonds that can break and reform, leading to stress relaxation without permanent deformation.

  2. Yamamoto (1956): Provided a more rigorous kinetic theory based on statistical mechanics, deriving the evolution equations for network structure under deformation.

  3. Lodge (1956): Developed the network model framework emphasizing the role of entanglements and temporary junctions in polymer melts and concentrated solutions.

  4. Tanaka & Edwards (1992): Formulated the conformation tensor approach used in modern TNT models, providing a unified framework for transient networks with various kinetic mechanisms. Their work established the connection between molecular-scale bond dynamics and macroscopic rheological response.

  5. Modern developments: The TNT framework has been extended to incorporate force-dependent breakage (Bell model), finite chain extensibility (FENE-P), multiple species, and Rouse dynamics for chain relaxation.

The Tanaka-Edwards model is mathematically equivalent to the upper-convected Maxwell (UCM) model, but derives from molecular network kinetics rather than continuum mechanical arguments. This molecular foundation provides physical interpretation of the parameters and enables systematic extensions to more complex network behaviors.

Physical Foundations

Mechanical Analogue

The Tanaka-Edwards model can be represented as a simple mechanical analogue:

┌─────────┐
│ Spring  │──── Network contribution (modulus G)
│  (G)    │
└────┬────┘
     │
┌────┴────┐
│ Dashpot │──── Network relaxation (viscosity η_p = G·τ_b)
│ (η_p)   │
└─────────┘
     ║
     ║  (parallel)
     ║
┌─────────┐
│ Dashpot │──── Solvent contribution (viscosity η_s)
│ (η_s)   │
└─────────┘

This is equivalent to the Jeffreys model or Oldroyd-B model (when solvent viscosity is included). The spring represents the elastic response of stretched network chains. The series dashpot represents stress relaxation due to bond breaking (at rate \(1/\tau_b\)). The parallel dashpot represents the contribution of the solvent and any free (non-network) chains.

Physical interpretation:

  • Spring (G): When chains are stretched, they store elastic energy. The modulus \(G \approx n_{chains} k_B T\) where \(n_{chains}\) is the number density of elastically active chains.

  • Series dashpot (η_p = G·τ_b): Bonds break randomly at rate \(1/\tau_b\), releasing stored elastic energy and allowing chains to relax. The longer the bond lifetime, the higher the viscosity contribution from the network.

  • Parallel dashpot (η_s): Solvent and free chains provide immediate viscous resistance to flow, independent of network dynamics.

Network Kinetics

The transient network is characterized by the following kinetic processes:

Bond breakage and reformation:

  • Breakage rate: \(1/\tau_b\) (constant, independent of chain conformation)

  • Reformation rate: \(1/\tau_b\) (maintains equilibrium at zero stress)

  • Active chain fraction: \(\phi = 1\) (all chains are connected, constant)

Kinetic balance at equilibrium:

At equilibrium (zero applied stress), bonds break and reform at equal rates:

\[\text{breakage rate} = \text{reformation rate} = \frac{1}{\tau_b}\]

This maintains a constant network structure with isotropic conformation \(\mathbf{S} = \mathbf{I}\).

Under deformation:

When the material is deformed, the conformation tensor \(\mathbf{S}\) deviates from \(\mathbf{I}\). Chains become stretched and oriented:

  1. Bond breaks: The chain instantly relaxes to equilibrium conformation \(\mathbf{S} = \mathbf{I}\) and contributes \(-(\mathbf{S} - \mathbf{I})/\tau_b\) to the rate of change of the average conformation.

  2. Bond forms: The chain joins the network with the current equilibrium conformation \(\mathbf{I}\) and contributes \(+(\mathbf{S} - \mathbf{I})/\tau_b\) to restore equilibrium.

  3. Deformation: The velocity gradient \(\boldsymbol{\kappa}\) continuously stretches and rotates the conformation tensor, contributing \(\boldsymbol{\kappa} \cdot \mathbf{S} + \mathbf{S} \cdot \boldsymbol{\kappa}^T\).

The competition between deformation (which stretches chains) and bond breaking (which relaxes chains) determines the steady-state network structure and stress.

Constant breakage assumption:

The key simplification of the Tanaka-Edwards model is that the breakage rate \(1/\tau_b\) is independent of chain stretch. This is appropriate when:

  • Bonds are weak enough that thermal fluctuations dominate over mechanical stress

  • The activation energy for bond breaking is independent of chain conformation

  • The applied stress is much smaller than the bond dissociation energy

For materials where stress significantly affects bond lifetime (e.g., strong shear thinning), force-dependent breakage models (Bell, catch-slip) are required.

Conformation Tensor

The conformation tensor \(\mathbf{S}\) is the normalized second moment of the end-to-end vector \(\mathbf{R}\) of network chains:

\[\mathbf{S} = \frac{\langle \mathbf{R} \otimes \mathbf{R} \rangle}{R_0^2}\]

where \(R_0 = \sqrt{\langle R^2 \rangle_0}\) is the equilibrium root-mean-square end-to-end distance (related to the number of Kuhn segments per chain).

Physical meaning:

  • \(\mathbf{S} = \mathbf{I}\): Chains are unstretched and isotropically oriented (equilibrium state)

  • \(S_{ii} > 1\): Chains are stretched along direction \(i\)

  • \(S_{ii} < 1\): Chains are compressed along direction \(i\)

  • \(S_{ij} \neq 0\) (off-diagonal): Chains are oriented at an angle to the principal axes

In simple shear flow:

For shear in the \(xy\)-plane with velocity \(v_x = \dot{\gamma} y\):

  • \(S_{xx}\): Chain extension in flow direction (increases with \(Wi = \tau_b \dot{\gamma}\))

  • \(S_{yy}\): Chain extension in gradient direction (decreases slightly)

  • \(S_{zz}\): Chain extension in vorticity direction (unchanged for 2D flow, \(S_{zz} = 1\))

  • \(S_{xy}\): Chain orientation (proportional to shear stress)

Stress-conformation relation:

The stress tensor is proportional to the deviation of \(\mathbf{S}\) from equilibrium:

\[\boldsymbol{\sigma}_{network} = G (\mathbf{S} - \mathbf{I})\]

This is the Giesekus-type constitutive relation, derived from entropic elasticity of network chains. The total stress includes the solvent contribution:

\[\boldsymbol{\sigma}_{total} = G (\mathbf{S} - \mathbf{I}) + 2 \eta_s \mathbf{D}\]

where \(\mathbf{D} = (\boldsymbol{\kappa} + \boldsymbol{\kappa}^T)/2\) is the rate of deformation tensor.

Governing Equations

Constitutive Equation

The evolution of the conformation tensor is governed by the upper-convected derivative:

\[\frac{d\mathbf{S}}{dt} = \boldsymbol{\kappa} \cdot \mathbf{S} + \mathbf{S} \cdot \boldsymbol{\kappa}^T - \frac{\mathbf{S} - \mathbf{I}}{\tau_b}\]

where:

  • \(\frac{d\mathbf{S}}{dt}\) is the material derivative (rate of change following the flow)

  • \(\boldsymbol{\kappa} \cdot \mathbf{S} + \mathbf{S} \cdot \boldsymbol{\kappa}^T\) is the upper-convected term (describes affine deformation and rotation of the conformation tensor with the flow)

  • \(-(\mathbf{S} - \mathbf{I})/\tau_b\) is the relaxation term (bond breaking returns chains to equilibrium at rate \(1/\tau_b\))

Affine deformation assumption:

The upper-convected derivative assumes that network chains deform affinely with the bulk flow, i.e., the end-to-end vector \(\mathbf{R}\) transforms as a material line element. This is valid when:

  • Chains are much smaller than any length scale of flow variation

  • Chains are well-connected to the surrounding network

  • Chain relaxation is much slower than local flow rearrangements

2D simple shear flow:

For simple shear with \(\boldsymbol{\kappa} = \dot{\gamma} \mathbf{e}_x \otimes \mathbf{e}_y\), the conformation tensor evolution reduces to:

\[ \begin{align}\begin{aligned}\frac{dS_{xx}}{dt} &= 2 \dot{\gamma} S_{xy} - \frac{S_{xx} - 1}{\tau_b}\\\frac{dS_{yy}}{dt} &= -\frac{S_{yy} - 1}{\tau_b}\\\frac{dS_{zz}}{dt} &= -\frac{S_{zz} - 1}{\tau_b}\\\frac{dS_{xy}}{dt} &= \dot{\gamma} S_{yy} - \frac{S_{xy}}{\tau_b}\end{aligned}\end{align} \]

Total stress:

The total shear stress and normal stress differences are:

\[ \begin{align}\begin{aligned}\sigma_{xy} &= G S_{xy} + \eta_s \dot{\gamma}\\N_1 &= \sigma_{xx} - \sigma_{yy} = G (S_{xx} - S_{yy})\\N_2 &= \sigma_{yy} - \sigma_{zz} = G (S_{yy} - S_{zz})\end{aligned}\end{align} \]

For the Tanaka-Edwards model with constant breakage, \(N_2 = 0\) in simple shear (since \(S_{yy} = S_{zz}\)).

Steady-State Solutions

At steady state, \(d\mathbf{S}/dt = \mathbf{0}\). Solving the ODEs for simple shear:

Conformation tensor components:

\[ \begin{align}\begin{aligned}S_{xy} &= \tau_b \dot{\gamma}\\S_{yy} &= 1\\S_{zz} &= 1\\S_{xx} &= 1 + 2 (\tau_b \dot{\gamma})^2\end{aligned}\end{align} \]

Shear stress (flow curve):

\[\sigma_{xy} = G \tau_b \dot{\gamma} + \eta_s \dot{\gamma} = \eta_0 \dot{\gamma}\]

where \(\eta_0 = G \tau_b + \eta_s\) is the zero-shear viscosity. The flow curve is Newtonian (constant viscosity).

First normal stress difference:

\[N_1 = G (S_{xx} - S_{yy}) = 2 G (\tau_b \dot{\gamma})^2\]

This is quadratic in shear rate, consistent with Maxwell-type viscoelasticity. The normal stress coefficient \(\Psi_1 = N_1 / \dot{\gamma}^2 = 2 G \tau_b^2\) is constant.

Second normal stress difference:

\[N_2 = G (S_{yy} - S_{zz}) = 0\]

The Tanaka-Edwards model predicts zero second normal stress difference in simple shear.

Physical Insight: \(N_2 = 0\)

The prediction \(N_2 = 0\) in simple shear is a direct consequence of the upper-convected derivative (affine kinematics). The gradient direction (\(y\)) and vorticity direction (\(z\)) are equivalent because the affine deformation acts only through \(\boldsymbol{\kappa}\), which has no \(yy\) or \(zz\) components in simple shear. Therefore \(S_{yy}(t) = S_{zz}(t)\) for all \(t\), giving \(N_2 = G(S_{yy} - S_{zz}) = 0\).

Experimental observation of \(N_2 \neq 0\) is a clear indicator that non-affine deformation is present. See TNT Non-Affine (Gordon-Schowalter) — Handbook for the Gordon-Schowalter extension that breaks this symmetry.

Stress Relaxation Non-Exponentiality

For the base Tanaka-Edwards model, stress relaxation is a perfect single exponential \(G(t) = G e^{-t/\tau_b}\). However, even within the TNT framework, non-exponential relaxation arises naturally from several extensions:

  • Bell breakage: Initially stretched chains have shorter \(\tau_b^{\text{eff}}\), producing fast initial decay that slows as chains relax toward equilibrium. This mimics a broad spectrum without requiring multiple modes.

  • Multi-mode (StickyRouse, MultiSpecies): Explicit multi-exponential decay from discrete \(\{G_k, \tau_k\}\) spectrum.

  • Cates model: \(G(t) \sim \exp(-\sqrt{2t/\tau_b})\) in the fast-breaking regime — a stretched exponential from scission/recombination dynamics.

Weissenberg number dependence:

The dimensionless Weissenberg number \(Wi = \tau_b \dot{\gamma}\) controls the relative magnitude of elastic and viscous stresses:

  • \(Wi \ll 1\): Viscous-dominated, \(\sigma \approx \eta_0 \dot{\gamma}\)

  • \(Wi \sim 1\): Comparable elastic and viscous contributions

  • \(Wi \gg 1\): Elastic-dominated, \(N_1 \gg \sigma_{xy}\)

SAOS Response

For small-amplitude oscillatory shear (SAOS), the conformation tensor is linearized around equilibrium \(\mathbf{S} = \mathbf{I}\):

\[\mathbf{S} = \mathbf{I} + \delta \mathbf{S} e^{i \omega t}\]

Substituting into the evolution equation and solving for the complex modulus:

\[G^*(\omega) = G' + i G''\]

where the storage and loss moduli are:

\[ \begin{align}\begin{aligned}G'(\omega) &= G \frac{(\omega \tau_b)^2}{1 + (\omega \tau_b)^2}\\G''(\omega) &= G \frac{\omega \tau_b}{1 + (\omega \tau_b)^2} + \eta_s \omega\end{aligned}\end{align} \]

Low-frequency behavior (\(\omega \tau_b \ll 1\)):

\[ \begin{align}\begin{aligned}G'(\omega) &\approx G (\omega \tau_b)^2 \quad \text{(quadratic)}\\G''(\omega) &\approx G \omega \tau_b + \eta_s \omega = \eta_0 \omega \quad \text{(linear)}\end{aligned}\end{align} \]

This is characteristic of Maxwell-type viscoelasticity with a single relaxation mode.

High-frequency behavior (\(\omega \tau_b \gg 1\)):

\[ \begin{align}\begin{aligned}G'(\omega) &\approx G \quad \text{(plateau)}\\G''(\omega) &\approx \frac{G}{\omega \tau_b} + \eta_s \omega \quad \text{(decreasing + solvent)}\end{aligned}\end{align} \]

The storage modulus approaches the network modulus \(G\), representing the elastic response of the network before bonds have time to break.

Crossover frequency:

The loss tangent \(\tan \delta = G'' / G'\) equals 1 at the crossover frequency:

\[\omega_c = \frac{1}{\tau_b}\]

This provides a direct experimental measure of the bond lifetime.

Complex viscosity:

\[\eta^*(\omega) = \frac{G^*(\omega)}{i \omega} = \eta' - i \eta''\]

where:

\[ \begin{align}\begin{aligned}\eta'(\omega) &= \frac{G'(\omega)}{\omega} = G \frac{\omega \tau_b}{1 + (\omega \tau_b)^2} + \eta_s\\\eta''(\omega) &= \frac{G''(\omega)}{\omega} = G \frac{(\omega \tau_b)^2}{\omega [1 + (\omega \tau_b)^2]}\end{aligned}\end{align} \]

At low frequencies, \(\eta'(\omega \to 0) = \eta_0\).

Relaxation

For stress relaxation after a step shear strain \(\gamma_0\) applied at \(t = 0\):

Initial condition:

\[\mathbf{S}(t=0) = \mathbf{I} + \gamma_0 (\mathbf{e}_x \otimes \mathbf{e}_y + \mathbf{e}_y \otimes \mathbf{e}_x)\]

so \(S_{xy}(0) = \gamma_0\) and all other components are at equilibrium.

Evolution:

For \(t > 0\), with \(\boldsymbol{\kappa} = \mathbf{0}\) (no further deformation):

\[\frac{dS_{xy}}{dt} = -\frac{S_{xy}}{\tau_b}\]

Solution:

\[S_{xy}(t) = \gamma_0 e^{-t/\tau_b}\]

Stress relaxation:

\[\sigma(t) = G S_{xy}(t) = G \gamma_0 e^{-t/\tau_b}\]

This is a single exponential decay with relaxation time \(\tau_b\). The relaxation modulus is:

\[G(t) = \frac{\sigma(t)}{\gamma_0} = G e^{-t/\tau_b}\]

Characteristic times:

  • Time to relax to \(1/e\) of initial stress: \(t = \tau_b\)

  • Time to relax to 5% of initial stress: \(t \approx 3 \tau_b\)

  • Time to relax to 1% of initial stress: \(t \approx 4.6 \tau_b\)

Startup Flow

For startup of steady shear at constant \(\dot{\gamma}\) starting from equilibrium at \(t = 0\):

Initial condition:

\[\mathbf{S}(t=0) = \mathbf{I}\]

Evolution:

The conformation tensor evolves according to the ODEs derived above. For shear stress, the analytical solution is:

\[\sigma(t) = G \tau_b \dot{\gamma} [1 - e^{-t/\tau_b}] + \eta_s \dot{\gamma}\]

This can be rewritten as:

\[\sigma(t) = \eta_0 \dot{\gamma} [1 - \frac{G \tau_b}{\eta_0} e^{-t/\tau_b}]\]

Limiting behavior:

  • Initial slope (\(t \to 0\)): \(\frac{d\sigma}{dt}\bigg|_{t=0} = G \dot{\gamma}\) (elastic response)

  • Steady state (\(t \to \infty\)): \(\sigma_{ss} = \eta_0 \dot{\gamma}\) (Newtonian flow)

  • Time to steady state: \(t \approx 5 \tau_b\) (>99% of steady state)

Stress overshoot:

For the Tanaka-Edwards model with constant breakage, there is no stress overshoot in startup shear. The stress increases monotonically from zero to the steady-state value. Stress overshoot requires either:

  • Force-dependent breakage (shear thinning → overshoot)

  • Finite extensibility (strain hardening → overshoot)

  • Nonlinear damping (Giesekus-type nonlinearity)

Startup Phase Decomposition

The startup transient can be understood in three phases:

  1. Phase 1 — Elastic loading (\(t \ll \tau_b\)): Chains stretch affinely with the flow. \(S_{xy} \approx \dot{\gamma} t\) grows linearly. Stress rises as \(\sigma \approx G \dot{\gamma} t\) (purely elastic). No bonds have yet broken.

  2. Phase 2 — Bond turnover onset (\(t \sim \tau_b\)): Bonds begin breaking at rate \(1/\tau_b\). For force-dependent variants (Bell), \(\text{tr}(\mathbf{S})\) rises and \(\tau_b^{\text{eff}}\) decreases, accelerating destruction beyond creation. This produces a stress overshoot in Bell/FENE variants (but not in the base Tanaka-Edwards model with constant breakage).

  3. Phase 3 — Steady state (\(t \gg \tau_b\)): Formation-destruction balance is reached. The conformation tensor settles to its steady-state value. For the base model, this is \(S_{xy} = \tau_b \dot{\gamma}\).

Note

The base Tanaka-Edwards model reaches steady state monotonically (no overshoot). A stress overshoot requires force-dependent breakage (TNT Bell (Force-Dependent Breakage) — Handbook), finite extensibility (TNT FENE-P (Finite Extensibility) — Handbook), or non-affine slip (TNT Non-Affine (Gordon-Schowalter) — Handbook).

First normal stress difference in startup:

The analytical solution for \(N_1(t)\) involves both \(S_{xx}\) and \(S_{yy}\):

\[N_1(t) = 2 G (\tau_b \dot{\gamma})^2 [1 - e^{-t/\tau_b} - \frac{t}{\tau_b} e^{-t/\tau_b}]\]

The normal stress exhibits a transient overshoot before approaching the steady-state value \(N_1^{ss} = 2 G (\tau_b \dot{\gamma})^2\). The overshoot occurs at \(t_{max} \approx 2 \tau_b\).

Creep

For creep under constant applied stress \(\sigma_0\), the shear rate \(\dot{\gamma}(t)\) and strain \(\gamma(t)\) must be solved from the coupled ODEs:

Governing equations:

\[ \begin{align}\begin{aligned}\frac{dS_{xy}}{dt} &= \dot{\gamma} S_{yy} - \frac{S_{xy}}{\tau_b}\\\frac{dS_{yy}}{dt} &= -\frac{S_{yy} - 1}{\tau_b}\\\sigma_0 &= G S_{xy} + \eta_s \dot{\gamma}\end{aligned}\end{align} \]

Shear rate from stress balance:

\[\dot{\gamma}(t) = \frac{\sigma_0 - G S_{xy}(t)}{\eta_s}\]

Strain:

\[\gamma(t) = \int_0^t \dot{\gamma}(t') dt'\]

Analytical solution:

The creep compliance \(J(t) = \gamma(t) / \sigma_0\) is:

\[J(t) = \frac{1}{G} [1 - e^{-t/\tau_b}] + \frac{t}{\eta_0}\]

This consists of:

  1. Elastic contribution: \(J_e = 1/G\) (instantaneous elastic deformation, approached at \(t \approx 5\tau_b\))

  2. Viscous contribution: \(J_v(t) = t/\eta_0\) (linear in time, steady-state flow)

Limiting behavior:

  • Short time (\(t \ll \tau_b\)): \(J(t) \approx t/G \tau_b = t/\eta_p\) (network relaxation)

  • Long time (\(t \gg \tau_b\)): \(J(t) \approx 1/G + t/\eta_0\) (steady flow with elastic offset)

Creep recovery:

If stress is removed at \(t = t_1\), the recoverable strain is:

\[\gamma_{rec} = \frac{\sigma_0}{G} [1 - e^{-t_1/\tau_b}]\]

The non-recoverable (viscous) strain is:

\[\gamma_{visc} = \frac{\sigma_0 t_1}{\eta_0}\]

Creep Rupture (Force-Dependent Extension)

While the base Tanaka-Edwards model with constant breakage always reaches a finite steady-state creep rate \(\dot{\gamma}_{ss} = \sigma_0/\eta_0\), the introduction of force-dependent breakage (Bell variant) creates the possibility of creep rupture:

  1. Applied stress \(\sigma_0\) stretches chains → \(\text{tr}(\mathbf{S})\) increases

  2. Stretch increases \(k_{\text{off}}\) (Bell) → effective \(\tau_b\) decreases

  3. Reduced \(\tau_b\) means lower effective viscosity → \(\dot{\gamma}\) increases

  4. Higher \(\dot{\gamma}\) stretches chains further → positive feedback loop

  5. Above a critical stress \(\sigma_c\), the system undergoes delayed yielding (creep rate accelerates without bound)

This viscosity bifurcation is analogous to the yielding transition in thixotropic fluids. See TNT Bell (Force-Dependent Breakage) — Handbook for the full treatment.

LAOS

For large-amplitude oscillatory shear (LAOS) with strain \(\gamma(t) = \gamma_0 \sin(\omega t)\):

Strain rate:

\[\dot{\gamma}(t) = \gamma_0 \omega \cos(\omega t)\]

Governing ODEs:

The conformation tensor components evolve according to:

\[ \begin{align}\begin{aligned}\frac{dS_{xx}}{dt} &= 2 \gamma_0 \omega \cos(\omega t) S_{xy} - \frac{S_{xx} - 1}{\tau_b}\\\frac{dS_{yy}}{dt} &= -\frac{S_{yy} - 1}{\tau_b}\\\frac{dS_{xy}}{dt} &= \gamma_0 \omega \cos(\omega t) S_{yy} - \frac{S_{xy}}{\tau_b}\end{aligned}\end{align} \]

Stress response:

\[\sigma(t) = G S_{xy}(t) + \eta_s \gamma_0 \omega \cos(\omega t)\]

Fourier decomposition:

The stress can be decomposed into odd harmonics:

\[\sigma(t) = \sum_{n=1,3,5,\ldots} [G'_n(\omega, \gamma_0) \sin(n \omega t) + G''_n(\omega, \gamma_0) \cos(n \omega t)]\]

where \(G'_1\) and \(G''_1\) are the linear viscoelastic moduli (at small \(\gamma_0\)), and higher harmonics \(G'_3, G''_3, G'_5, \ldots\) characterize nonlinearity.

Lissajous curves:

Plots of \(\sigma\) vs. \(\gamma\) (elastic Lissajous) and \(\sigma\) vs. \(\dot{\gamma}\) (viscous Lissajous) reveal:

  • Linear regime (\(\gamma_0 \ll 1\)): Elliptical curves

  • Nonlinear regime (\(\gamma_0 \sim 1\)): Distorted ellipses with higher harmonics

Numerical solution:

LAOS requires numerical integration of the ODEs over multiple cycles until a periodic steady state is reached. The RheoJAX implementation uses Diffrax with adaptive time stepping.

LAOS: Linear Stress Response at All Amplitudes

A remarkable property of the base Tanaka-Edwards model (constant breakage, Hookean chains) is that the stress response in LAOS is perfectly sinusoidal regardless of strain amplitude \(\gamma_0\). This means:

  • \(I_3/I_1 = 0\) (no third harmonic) at all \(\gamma_0\)

  • The elastic Lissajous curve is always a perfect ellipse

  • The Fourier spectrum contains only the fundamental frequency

This occurs because the constitutive equation is linear in \(\mathbf{S}\) (the upper-convected Maxwell model). The time-varying shear rate \(\dot{\gamma}(t) = \gamma_0 \omega \cos(\omega t)\) enters linearly through the velocity gradient \(\boldsymbol{\kappa}\), and the relaxation term is also linear.

Consequence for model identification: If LAOS data show any higher harmonics (\(I_3/I_1 > 0\)), the base Tanaka-Edwards model is immediately ruled out. One must invoke nonlinear extensions:

  • Bell breakage → odd harmonics from \(\exp[\nu(\text{tr}\mathbf{S} - 3)]\)

  • FENE stress → harmonics from \(L^2/(L^2 - \text{tr}\mathbf{S})\) nonlinearity

  • Non-affine slip → harmonics from Gordon-Schowalter coupling

Parameters

Tanaka-Edwards Model Parameters

Parameter

Symbol

Default

Bounds

Units

Description

G

\(G\)

1000

(1, 108)

Pa

Network elastic modulus (related to chain density via \(G \approx n_{chains} k_B T\))

tau_b

\(\tau_b\)

1.0

(10-6, 104)

s

Mean bond lifetime (characteristic relaxation time of the network)

\(\eta_s\)

\(\eta_s\)

0.0

(0, 104)

Pa·s

Solvent viscosity (non-network contribution, can be zero for entangled systems)

Parameter Interpretation

Network modulus (G):

The network modulus is related to the number density of elastically active chains:

\[G \approx n_{chains} k_B T\]

where:

  • \(n_{chains}\) is the number of chains per unit volume (m-3)

  • \(k_B = 1.38 \times 10^{-23}\) J/K is Boltzmann’s constant

  • \(T\) is absolute temperature (K)

For a polymer solution with concentration \(c\) (g/mL) and molecular weight \(M_w\) (g/mol):

\[n_{chains} \approx \frac{c N_A}{M_w}\]

where \(N_A = 6.02 \times 10^{23}\) mol-1 is Avogadro’s number.

Typical values:

  • Dilute polymer solutions: \(G \sim 1-100\) Pa

  • Semi-dilute solutions: \(G \sim 100-10^4\) Pa

  • Entangled melts and gels: \(G \sim 10^4-10^6\) Pa

Bond lifetime (τ_b):

The bond lifetime \(\tau_b\) is the mean time between bond breaking and reformation events. It is related to the activation energy \(E_a\) for bond dissociation via Arrhenius kinetics:

\[\tau_b = \tau_0 \exp\left(\frac{E_a}{k_B T}\right)\]

where \(\tau_0\) is a pre-exponential attempt time (typically 10-9 to 10-12 s).

Experimental determination:

  • SAOS: Crossover frequency \(\omega_c = 1/\tau_b\)

  • Stress relaxation: Time to decay to \(1/e\) of initial stress

  • DLS: Characteristic time for structural relaxation

  • Fluorescence recovery: Bond exchange timescale

Typical values:

  • Weak physical gels: \(\tau_b \sim 10^{-3}-1\) s

  • Associating polymers: \(\tau_b \sim 0.1-100\) s

  • Strong supramolecular networks: \(\tau_b \sim 10-10^4\) s

Solvent viscosity (η_s):

The solvent viscosity \(\eta_s\) accounts for viscous dissipation not associated with network relaxation. It includes:

  • Viscosity of the solvent (e.g., water ~ 0.001 Pa·s at 20°C)

  • Contribution from free (non-network) chains

  • Hydrodynamic interactions

For entangled polymer melts without added solvent, \(\eta_s \approx 0\) is a reasonable approximation.

Zero-shear viscosity:

The zero-shear viscosity is the sum of network and solvent contributions:

\[\eta_0 = G \tau_b + \eta_s\]

This can be measured directly from steady-state flow curves (Newtonian plateau at low \(\dot{\gamma}\)).

Parameter correlations:

  • \(G\) and \(\tau_b\) are strongly correlated when fitting SAOS or relaxation data (both affect the magnitude and timescale of the response)

  • \(\eta_s\) is most easily determined from high-frequency (or high shear rate) behavior where network relaxation is fast

Bayesian priors:

For Bayesian inference, recommended priors are:

\[ \begin{align}\begin{aligned}G &\sim \text{LogNormal}(\log(1000), 2) \quad \text{(broad, uninformative)}\\\tau_b &\sim \text{LogNormal}(\log(1), 2) \quad \text{(centered at 1 s)}\\\eta_s &\sim \text{HalfNormal}(10) \quad \text{(weakly informative, allows zero)}\end{aligned}\end{align} \]

These can be refined based on prior knowledge of the material system.

Validity and Assumptions

The Tanaka-Edwards model is valid under the following assumptions:

1. Constant breakage rate:

  • Bond lifetime \(\tau_b\) is independent of chain stretch or applied force

  • Appropriate when thermal fluctuations dominate over mechanical stress

  • Fails for materials with strong force-dependent bond kinetics (use Bell or catch-slip variants instead)

2. Affine deformation:

  • Network chains deform affinely with the bulk flow (no slip)

  • Valid when chains are well-connected and much smaller than flow length scales

  • Fails for poorly connected networks or spatially heterogeneous deformation

3. Instant reformation at equilibrium:

  • When a bond breaks, the chain instantly relaxes to \(\mathbf{S} = \mathbf{I}\)

  • When a bond forms, the chain joins with equilibrium conformation

  • Neglects finite Rouse relaxation time (use TNTStickyRouse for chain dynamics)

4. Linear springs (Hookean elasticity):

  • No finite extensibility limit for network chains

  • Stress is proportional to \(\mathbf{S} - \mathbf{I}\) at all strain levels

  • Fails at high strains where chains approach full extension (use FENE-P variant)

5. Monodisperse network:

  • Single bond lifetime \(\tau_b\) (single relaxation mode)

  • All chains have the same \(G\) and \(R_0\)

  • For polydisperse systems with multiple relaxation times, use TNTMultiSpecies or sum of Maxwell modes

6. No entanglements beyond crosslinks:

  • Chain relaxation occurs only via bond breaking, not reptation or constraint release

  • Appropriate for dilute to semi-dilute solutions

  • For entangled melts, additional relaxation mechanisms may be important

7. Homogeneous deformation:

  • Assumes no shear banding or spatial gradients in structure

  • Fails when flow induces heterogeneous microstructure (use 1D nonlocal model if needed)

8. Incompressibility:

  • Volume is conserved (\(\text{tr}(\boldsymbol{\kappa}) = 0\))

  • Appropriate for polymer solutions and melts

  • May fail for compressible systems (e.g., foams)

Regimes and Behavior

The Tanaka-Edwards model exhibits distinct rheological regimes depending on the dimensionless Weissenberg number \(Wi = \tau_b \dot{\gamma}\) (for flow) or dimensionless frequency \(\omega \tau_b\) (for oscillatory shear).

Linear viscoelastic regime (\(Wi \ll 1\) or \(\omega \tau_b \ll 1\)):

  • SAOS: \(G' \sim \omega^2\), \(G'' \sim \omega\) (terminal regime)

  • Flow: Newtonian behavior, \(\sigma = \eta_0 \dot{\gamma}\)

  • Normal stress: Negligible, \(N_1 \ll \sigma\)

  • Physical picture: Bonds break and reform faster than chains can be significantly stretched

Crossover regime (\(Wi \sim 1\) or \(\omega \tau_b \sim 1\)):

  • SAOS: \(G' \approx G'' \approx G/2\), \(\tan \delta = 1\)

  • Flow: Comparable elastic and viscous stresses

  • Normal stress: \(N_1 \sim \sigma^2 / G\)

  • Physical picture: Deformation and relaxation timescales are matched

Elastic-dominated regime (\(Wi \gg 1\) or \(\omega \tau_b \gg 1\)):

  • SAOS: \(G' \to G\) (rubber-like plateau), \(G'' \ll G'\)

  • Flow: Normal stress dominates, \(N_1 \gg \sigma\)

  • Physical picture: Chains are highly stretched before bonds break, elastic energy storage dominates

Key behavioral features:

  1. Newtonian flow curve: The Tanaka-Edwards model always predicts constant viscosity \(\eta_0\) at steady state. This is a key limitation for materials that exhibit shear thinning. For shear-thinning behavior, use force-dependent breakage variants (Bell, catch-slip).

  2. Quadratic normal stress: \(N_1 \propto \dot{\gamma}^2\) is a universal prediction of single-mode Maxwell models. Experimental data showing different scaling (e.g., \(N_1 \propto \dot{\gamma}^{1.5}\) for some wormlike micelles) indicates more complex network dynamics.

  3. Single relaxation time: All transient responses (relaxation, startup, creep) are characterized by the single timescale \(\tau_b\). Materials with multiple relaxation times require multi-mode extensions.

  4. No stress overshoot in startup shear: Unlike shear-thinning models, the Tanaka-Edwards model predicts monotonic stress increase. Experimental overshoot indicates force-dependent kinetics or finite extensibility effects.

  5. Zero second normal stress difference: \(N_2 = 0\) in simple shear. Experimental observation of \(N_2 < 0\) indicates non-affine deformation or additional microstructural effects.

What You Can Learn

Fitting the Tanaka-Edwards model to experimental data provides the following physical insights:

From SAOS (:math:`G’`, :math:`G’’` vs. ω):

  1. Network modulus (G): High-frequency plateau of \(G'\) gives \(G\) directly. This is related to chain density via \(G \approx n_{chains} k_B T\), allowing estimation of the number of elastically active chains per unit volume.

  2. Bond lifetime (τ_b): Crossover frequency \(\omega_c\) where \(G' = G''\) gives \(\tau_b = 1/\omega_c\). This is the characteristic timescale for bond breaking and network relaxation.

  3. Zero-shear viscosity (η_0): Low-frequency limit \(\eta_0 = G'' / \omega\) as \(\omega \to 0\) gives \(\eta_0 = G \tau_b + \eta_s\).

  4. Solvent viscosity (η_s): High-frequency viscous contribution \(\eta_s \omega\) determines the non-network viscosity.

From stress relaxation (σ vs. t):

  1. Bond lifetime (τ_b): Time to decay to \(1/e\) of initial stress is \(\tau_b\). This is the most direct measure of bond dynamics.

  2. Network modulus (G): Initial stress \(\sigma(0) = G \gamma_0\) after step strain \(\gamma_0\).

From steady-state flow (σ vs. γ̇):

  1. Zero-shear viscosity (η_0): Slope of the Newtonian flow curve.

  2. Model validity check: If the flow curve shows shear thinning (decreasing viscosity with \(\dot{\gamma}\)), the Tanaka-Edwards model is not appropriate. Switch to force-dependent breakage models (Bell, catch-slip).

From normal stress (N₁ vs. γ̇):

  1. Bond lifetime (τ_b): From \(N_1 = 2 G (\tau_b \dot{\gamma})^2\), the ratio \(N_1 / \sigma^2 = 2 \tau_b / \eta_0\) gives a second estimate of \(\tau_b\) that is independent of SAOS.

  2. Consistency check: If \(N_1 / \dot{\gamma}^2\) is not constant, the model assumptions are violated (force-dependent kinetics or finite extensibility).

From creep (γ vs. t):

  1. Compliance (J): Long-time slope \(1/\eta_0\) and elastic offset \(1/G\).

  2. Recoverable strain: Fraction of strain that recovers after stress removal is \(\gamma_{rec} / \gamma_{total} = (1 - e^{-t_1/\tau_b}) / [1 + G t_1 / (\eta_0 \tau_b)]\).

From LAOS (harmonics vs. γ_0):

  1. Onset of nonlinearity: Strain amplitude \(\gamma_0\) at which higher harmonics become significant. For the Tanaka-Edwards model, nonlinearity appears at \(\gamma_0 \sim 1\) (order-unity strain).

  2. Nonlinear parameter space: Odd harmonics \(G'_3, G''_3, \ldots\) characterize the shape of the stress-strain response.

Temperature dependence:

If data are available at multiple temperatures, Arrhenius analysis of \(\tau_b(T)\) gives the activation energy \(E_a\) for bond dissociation:

\[\ln(\tau_b) = \ln(\tau_0) + \frac{E_a}{k_B T}\]

This provides insight into the molecular mechanism of bond breaking (hydrogen bonding, hydrophobic interactions, electrostatic interactions, etc.).

Experimental Design

Recommended protocols for parameter extraction:

  1. SAOS (essential):

    • Frequency range: At least 2 decades around \(\omega_c = 1/\tau_b\)

    • Recommended: \(10^{-2}\) to \(10^2\) rad/s for typical \(\tau_b \sim 1\) s

    • Strain amplitude: \(\gamma_0 < 0.1\) (linear regime)

    • Output: \(G\), \(\tau_b\), \(\eta_s\)

  2. Stress relaxation (highly informative):

    • Step strain: \(\gamma_0 \sim 0.1-0.5\) (linear regime)

    • Time range: At least \(5 \tau_b\) to reach baseline

    • Output: \(\tau_b\) (most direct), \(G\)

  3. Steady-state flow (model validation):

    • Shear rate range: \(10^{-3}\) to \(10^2\) s-1 (to cover \(Wi < 1\) and \(Wi > 1\))

    • Check: Flow curve should be Newtonian (constant \(\eta_0\))

    • If shear thinning is observed: model is inappropriate, use Bell variant

  4. Normal stress (optional, for consistency check):

    • Same shear rate range as flow curve

    • Check: \(N_1 / \dot{\gamma}^2 = 2 G \tau_b^2\) should be constant

    • Provides independent estimate of \(\tau_b\)

  5. Startup flow (optional):

    • Useful for validating transient behavior

    • Apply constant \(\dot{\gamma}\) from equilibrium

    • For Tanaka-Edwards: no stress overshoot expected (monotonic rise)

  6. Creep (alternative to relaxation):

    • Apply constant \(\sigma_0 < G\) (linear regime)

    • Time range: At least \(5 \tau_b\) to observe steady flow

    • Output: \(G\), \(\eta_0\), \(\tau_b\)

  7. LAOS (for advanced characterization):

    • Strain amplitude sweep: \(\gamma_0 = 0.01\) to \(10\)

    • Fixed frequency: \(\omega \sim 1/\tau_b\) (crossover)

    • Output: Nonlinear parameters, Lissajous curves

Experimental considerations:

  • Sample equilibration: Allow sufficient rest time (\(\sim 10 \tau_b\)) between measurements to restore equilibrium structure

  • Temperature control: Maintain constant temperature (±0.1°C) as \(\tau_b\) is strongly temperature-dependent

  • Wall slip: Check for slip (especially at low \(\dot{\gamma}\)) using parallel plate gap variation or roughened surfaces

  • Edge fracture: At high \(Wi\), normal stresses can cause edge instabilities in parallel plate geometry

  • Instrument compliance: Subtract instrument inertia and compliance (especially important for low-modulus materials \(G < 100\) Pa)

Minimal experimental set for fitting:

For quick parameter estimation with limited instrument time:

  1. SAOS frequency sweep (\(\gamma_0 = 0.05\), \(10^{-2}\) to \(10^2\) rad/s)

  2. Stress relaxation (step \(\gamma_0 = 0.2\), measure for \(5\tau_b\))

These two experiments are sufficient to extract all three parameters and validate the model assumptions.

Computational Implementation

RheoJAX TNTSingleMode implementation:

The Tanaka-Edwards model is implemented as TNTSingleMode with breakage="constant":

from rheojax.models import TNTSingleMode

# Create model with constant breakage (Tanaka-Edwards)
model = TNTSingleMode(breakage="constant")

# Default parameters: G=1000 Pa, tau_b=1.0 s, eta_s=0.0 Pa·s
# Bounds: G in [1, 1e8], tau_b in [1e-6, 1e4], eta_s in [0, 1e4]

Test modes and predictions:

  1. OSCILLATION: Analytical solution for \(G'(\omega)\) and \(G''(\omega)\)

  2. RELAXATION: Analytical solution \(\sigma(t) = G \gamma_0 e^{-t/\tau_b}\)

  3. STARTUP: Analytical solution \(\sigma(t) = \eta_0 \dot{\gamma} [1 - (G \tau_b / \eta_0) e^{-t/\tau_b}]\)

  4. FLOW_CURVE: Analytical solution \(\sigma = \eta_0 \dot{\gamma}\)

  5. CREEP: ODE integration for \(\gamma(t)\)

  6. LAOS: ODE integration for \(\mathbf{S}(t)\) over multiple cycles

ODE solver details:

For protocols requiring numerical integration (CREEP, LAOS):

  • Solver: Diffrax Tsit5 (adaptive 5th-order Runge-Kutta)

  • Tolerances: rtol=1e-6, atol=1e-8 (configurable)

  • Time stepping: Adaptive with safety factor 0.9

  • JIT compilation: All functions are JIT-compiled for GPU acceleration

Conformation tensor representation:

The 3×3 symmetric conformation tensor \(\mathbf{S}\) is stored as a 6-component vector: [S_xx, S_yy, S_zz, S_xy, S_xz, S_yz]. For 2D simple shear, only the first 4 components are active.

Numerical stability:

  • At very high \(Wi\) (\(> 10^3\)), the conformation tensor can become stiff. The adaptive solver automatically reduces time step size to maintain accuracy.

  • For LAOS at high \(\gamma_0\) (\(> 10\)), the model may predict unphysically large chain extension (since there is no finite extensibility limit). Use the FENE-P variant if strain stiffening is observed experimentally.

Performance:

  • SAOS, RELAXATION, STARTUP, FLOW_CURVE: Analytical solutions, no ODE integration required. Prediction time: ~0.1-1 ms for 100 points (GPU).

  • CREEP, LAOS: ODE integration required. Prediction time: ~10-100 ms for 1000 time points (GPU, depends on stiffness and number of cycles).

Fitting workflow:

The recommended workflow uses NLSQ for fast point estimation followed by Bayesian inference for uncertainty quantification:

# Step 1: NLSQ fit (fast, deterministic)
model.fit(omega, G_star, test_mode='oscillation')
print(f"NLSQ: G={model.G}, tau_b={model.tau_b}, eta_s={model.eta_s}")

# Step 2: Bayesian fit with NLSQ warm-start
result = model.fit_bayesian(omega, G_star, test_mode='oscillation',
                             num_warmup=1000, num_samples=2000, num_chains=4)

# Step 3: Extract credible intervals
intervals = model.get_credible_intervals(result.posterior_samples, credibility=0.95)
print(f"95% CI for G: {intervals['G']}")

Custom parameter bounds:

For materials with known parameter ranges, custom bounds improve fitting:

from rheojax.core import Parameter

model = TNTSingleMode(breakage="constant")
model.parameters['G'].bounds = (100, 10000)  # Pa, for dilute solutions
model.parameters['tau_b'].bounds = (0.01, 100)  # s, for fast networks
model.fit(omega, G_star, test_mode='oscillation')

Fitting Guidance

SAOS fitting:

For small-amplitude oscillatory shear data:

import numpy as np
from rheojax.models import TNTSingleMode

# Prepare data: frequency, complex modulus
omega = np.logspace(-2, 2, 50)  # 50 points, 10^-2 to 10^2 rad/s
G_star = G_prime + 1j * G_double_prime  # Complex modulus

# Option 1: Fit to complex modulus (recommended)
model = TNTSingleMode(breakage="constant")
model.fit(omega, G_star, test_mode='oscillation')

# Option 2: Fit to G' and G'' separately (if complex data not available)
# Stack G' and G'' as [G'_1, ..., G'_N, G''_1, ..., G''_N]
y_data = np.concatenate([G_prime, G_double_prime])
model.fit(omega, y_data, test_mode='oscillation')

Tips for SAOS fitting:

  • Ensure frequency range covers at least 1 decade below and above \(\omega_c = 1/\tau_b\)

  • If \(\tau_b\) is unknown, use \(\omega = 10^{-3}\) to \(10^3\) rad/s

  • Weight low-frequency points more heavily if zero-shear viscosity is important

  • Check for instrument compliance (spurious upturn in \(G''\) at high \(\omega\))

Stress relaxation fitting:

For step strain relaxation:

# Prepare data: time, shear stress
t = np.linspace(0, 10, 200)  # 10 s, 200 points
sigma = ...  # Measured stress after step strain gamma_0

# Fit (note: gamma_0 must be provided as a keyword argument)
model = TNTSingleMode(breakage="constant")
model.fit(t, sigma, test_mode='relaxation', gamma_0=0.2)

# Extract parameters
print(f"tau_b = {model.tau_b:.3f} s")
print(f"G = {model.G:.0f} Pa")

Tips for relaxation fitting:

  • Measure for at least \(5\tau_b\) to reach baseline (stress decays to <1% of initial)

  • If \(\tau_b\) is unknown, measure for at least 100 s

  • Subtract residual stress (instrument drift, sample thixotropy) if non-zero at long times

  • The fit is relatively insensitive to \(\eta_s\) (relaxation is dominated by network)

Flow curve fitting:

For steady-state shear stress vs. shear rate:

# Prepare data: shear rate, shear stress
gamma_dot = np.logspace(-3, 2, 50)  # 10^-3 to 10^2 s^-1
sigma = ...  # Measured steady-state stress

# Fit
model = TNTSingleMode(breakage="constant")
model.fit(gamma_dot, sigma, test_mode='flow_curve')

# Expected: constant viscosity (Newtonian)
eta_0 = model.G * model.tau_b + model.eta_s
print(f"Zero-shear viscosity: {eta_0:.2f} Pa·s")

Important: If the flow curve shows shear thinning (decreasing \(\sigma / \dot{\gamma}\) with increasing \(\dot{\gamma}\)), the Tanaka-Edwards model is not appropriate. Switch to force-dependent breakage:

# For shear-thinning materials
model = TNTSingleMode(breakage="bell")  # Bell force-dependent breakage
model.fit(gamma_dot, sigma, test_mode='flow_curve')

Startup flow fitting:

For transient shear stress during startup:

# Prepare data: time, shear stress during startup at constant gamma_dot
t = np.linspace(0, 10, 200)
sigma = ...  # Measured stress during startup

# Fit (gamma_dot must be provided)
model = TNTSingleMode(breakage="constant")
model.fit(t, sigma, test_mode='startup', gamma_dot=1.0)

Creep fitting:

For strain vs. time under constant stress:

# Prepare data: time, shear strain
t = np.linspace(0, 100, 500)
gamma = ...  # Measured strain under constant sigma_0

# Fit (sigma_applied must be provided)
model = TNTSingleMode(breakage="constant")
model.fit(t, gamma, test_mode='creep', sigma_applied=100.0)

LAOS fitting:

For large-amplitude oscillatory shear:

# Prepare data: time, shear stress (multiple cycles)
t = np.linspace(0, 10*T, 1000)  # 10 cycles, T = 2*pi/omega
sigma = ...  # Measured stress

# Fit (gamma_0 and omega must be provided)
model = TNTSingleMode(breakage="constant")
model.fit(t, sigma, test_mode='laos', gamma_0=0.5, omega=1.0)

Bayesian inference with NLSQ warm-start:

For parameter uncertainty quantification:

# Step 1: NLSQ for initial guess
model = TNTSingleMode(breakage="constant")
model.fit(omega, G_star, test_mode='oscillation')

# Step 2: Bayesian with warm-start (4 chains for robust diagnostics)
result = model.fit_bayesian(omega, G_star, test_mode='oscillation',
                             num_warmup=1000,    # Burn-in
                             num_samples=2000,   # Samples per chain
                             num_chains=4,       # Parallel chains (default)
                             seed=42)            # Reproducibility

# Step 3: Check diagnostics
print(f"R-hat: {result.diagnostics['r_hat']}")  # Should be < 1.01
print(f"ESS: {result.diagnostics['ess']}")      # Should be > 400

# Step 4: Extract credible intervals
intervals = model.get_credible_intervals(result.posterior_samples, credibility=0.95)
print(f"G: {intervals['G'][0]:.0f} - {intervals['G'][1]:.0f} Pa (95% CI)")
print(f"tau_b: {intervals['tau_b'][0]:.3f} - {intervals['tau_b'][1]:.3f} s (95% CI)")

Parameter initialization:

For models that are difficult to fit (e.g., multiple relaxation modes, high noise), smart initialization improves convergence:

from rheojax.utils.initialization import initialize_from_saos

# Extract initial guess from SAOS data
initial_params = initialize_from_saos(omega, G_star, model_type='maxwell')

# Set initial values
model = TNTSingleMode(breakage="constant")
model.G = initial_params['G']
model.tau_b = initial_params['tau_b']

# Fit with improved initial guess
model.fit(omega, G_star, test_mode='oscillation')

Common fitting issues:

  1. Poor fit quality (low R²):

    • Check data quality (noise, instrument artifacts)

    • Verify model assumptions (Newtonian flow curve, single relaxation time)

    • Try multi-mode model if data show multiple relaxation times

  2. Parameter bounds hit:

    • If \(G\) hits upper bound (10⁸ Pa): check for instrument compliance or dimensionless data

    • If \(\tau_b\) hits lower bound (10⁻⁶ s): relaxation is too fast for instrument or model is inappropriate

    • If \(\eta_s\) hits upper bound: network relaxation may be too slow, or solvent dominates

  3. Bayesian convergence issues (R-hat > 1.1, low ESS):

    • Increase num_warmup (try 2000-5000 for difficult posteriors)

    • Check for multimodal posteriors (plot pair plots)

    • Tighten parameter bounds if physically justified

    • Ensure NLSQ warm-start is used (critical for convergence)

  4. Numerical instability in ODE integration (CREEP, LAOS):

    • Reduce rtol and atol (try 1e-8 and 1e-10)

    • Check for unphysical parameter values (negative \(\eta_s\), too large \(Wi\))

    • Use implicit solver for stiff systems (not yet implemented in RheoJAX)

Usage

Basic SAOS prediction and fitting:

from rheojax.models import TNTSingleMode
import numpy as np
import matplotlib.pyplot as plt

# Create model with default parameters
model = TNTSingleMode(breakage="constant")
print(f"Default: G={model.G} Pa, tau_b={model.tau_b} s, eta_s={model.eta_s} Pa·s")

# Generate SAOS prediction
omega = np.logspace(-2, 2, 100)  # 0.01 to 100 rad/s
G_star = model.predict(omega, test_mode='oscillation')

# Extract G' and G''
G_prime = G_star.real
G_double_prime = G_star.imag

# Plot
plt.figure(figsize=(8, 6))
plt.loglog(omega, G_prime, 'o-', label="G' (storage)")
plt.loglog(omega, G_double_prime, 's-', label="G'' (loss)")
plt.xlabel('Angular frequency (rad/s)')
plt.ylabel('Modulus (Pa)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.title('Tanaka-Edwards SAOS Response')
plt.show()

# Fit to synthetic data (with noise)
np.random.seed(42)
noise = 1 + 0.05 * np.random.randn(len(omega))  # 5% noise
G_star_data = G_star * noise

model_fit = TNTSingleMode(breakage="constant")
model_fit.fit(omega, G_star_data, test_mode='oscillation')
print(f"Fitted: G={model_fit.G:.0f} Pa, tau_b={model_fit.tau_b:.3f} s, eta_s={model_fit.eta_s:.3f} Pa·s")

Stress relaxation:

# Predict stress relaxation after step strain
t = np.linspace(0, 10, 200)  # 0 to 10 s
gamma_0 = 0.2  # Step strain (20%)

sigma = model.predict(t, test_mode='relaxation', gamma_0=gamma_0)

# Plot
plt.figure(figsize=(8, 6))
plt.semilogy(t, sigma, 'o-')
plt.axhline(sigma[0] / np.e, color='red', linestyle='--',
            label=f'τ_b = {model.tau_b:.2f} s')
plt.xlabel('Time (s)')
plt.ylabel('Stress (Pa)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.title('Stress Relaxation (Tanaka-Edwards)')
plt.show()

Startup flow:

# Predict stress during startup at constant shear rate
t = np.linspace(0, 5, 200)  # 0 to 5 s
gamma_dot = 1.0  # 1 s^-1

sigma = model.predict(t, test_mode='startup', gamma_dot=gamma_dot)

# Analytical steady-state
eta_0 = model.G * model.tau_b + model.eta_s
sigma_ss = eta_0 * gamma_dot

# Plot
plt.figure(figsize=(8, 6))
plt.plot(t, sigma, 'o-', label='Transient')
plt.axhline(sigma_ss, color='red', linestyle='--', label=f'Steady state ({sigma_ss:.1f} Pa)')
plt.xlabel('Time (s)')
plt.ylabel('Stress (Pa)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.title(f'Startup Flow (γ̇ = {gamma_dot} s⁻¹)')
plt.show()

Flow curve:

# Predict steady-state flow curve
gamma_dot = np.logspace(-3, 2, 50)  # 0.001 to 100 s^-1
sigma = model.predict(gamma_dot, test_mode='flow_curve')

# Calculate viscosity
eta = sigma / gamma_dot

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

ax1.loglog(gamma_dot, sigma, 'o-')
ax1.set_xlabel('Shear rate (s⁻¹)')
ax1.set_ylabel('Shear stress (Pa)')
ax1.grid(True, alpha=0.3)
ax1.set_title('Flow Curve (Tanaka-Edwards)')

ax2.semilogx(gamma_dot, eta, 'o-')
ax2.axhline(eta_0, color='red', linestyle='--', label=f'η₀ = {eta_0:.1f} Pa·s')
ax2.set_xlabel('Shear rate (s⁻¹)')
ax2.set_ylabel('Viscosity (Pa·s)')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_title('Viscosity (Newtonian)')

plt.tight_layout()
plt.show()

Normal stress difference:

# Predict first normal stress difference (requires custom calculation)
# For Tanaka-Edwards: N_1 = 2 * G * (tau_b * gamma_dot)^2

gamma_dot = np.logspace(-2, 2, 50)
N_1 = 2 * model.G * (model.tau_b * gamma_dot)**2

# Normal stress coefficient
Psi_1 = N_1 / gamma_dot**2

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

ax1.loglog(gamma_dot, N_1, 'o-')
ax1.set_xlabel('Shear rate (s⁻¹)')
ax1.set_ylabel('N₁ (Pa)')
ax1.grid(True, alpha=0.3)
ax1.set_title('First Normal Stress Difference')

ax2.semilogx(gamma_dot, Psi_1, 'o-')
ax2.set_xlabel('Shear rate (s⁻¹)')
ax2.set_ylabel('Ψ₁ (Pa·s²)')
ax2.grid(True, alpha=0.3)
ax2.set_title('Normal Stress Coefficient (Constant)')

plt.tight_layout()
plt.show()

Creep simulation:

# Predict creep under constant stress
t = np.linspace(0, 50, 500)  # 0 to 50 s
sigma_0 = 100.0  # Applied stress (Pa)

gamma = model.predict(t, test_mode='creep', sigma_applied=sigma_0)

# Analytical creep compliance
J_t = (1/model.G) * (1 - np.exp(-t/model.tau_b)) + t / eta_0
gamma_analytical = sigma_0 * J_t

# Plot
plt.figure(figsize=(8, 6))
plt.plot(t, gamma, 'o-', label='Numerical', markersize=3)
plt.plot(t, gamma_analytical, 'r--', label='Analytical')
plt.xlabel('Time (s)')
plt.ylabel('Strain')
plt.legend()
plt.grid(True, alpha=0.3)
plt.title(f'Creep (σ₀ = {sigma_0} Pa)')
plt.show()

LAOS simulation:

# Predict LAOS response
omega = 1.0  # rad/s
gamma_0 = 0.5  # Strain amplitude
n_cycles = 10
T = 2 * np.pi / omega
t = np.linspace(0, n_cycles * T, 1000)

# Note: For LAOS, predict() returns stress vs. time
# Strain must be computed separately
gamma_t = gamma_0 * np.sin(omega * t)
sigma = model.predict(t, test_mode='laos', gamma_0=gamma_0, omega=omega)

# Plot last 2 cycles (steady periodic state)
t_plot = t[t > (n_cycles - 2) * T]
gamma_plot = gamma_t[t > (n_cycles - 2) * T]
sigma_plot = sigma[t > (n_cycles - 2) * T]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Time series
ax1.plot(t_plot, sigma_plot, 'o-', markersize=3, label='Stress')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Stress (Pa)')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_title('LAOS Time Series')

# Lissajous (elastic)
ax2.plot(gamma_plot, sigma_plot, 'o-', markersize=3)
ax2.set_xlabel('Strain')
ax2.set_ylabel('Stress (Pa)')
ax2.grid(True, alpha=0.3)
ax2.set_title('Elastic Lissajous Curve')

plt.tight_layout()
plt.show()

Bayesian inference workflow:

from rheojax.pipeline.bayesian import BayesianPipeline

# Create synthetic SAOS data
omega_true = np.logspace(-2, 2, 30)
G_true, tau_b_true, eta_s_true = 1000, 1.0, 0.1
model_true = TNTSingleMode(breakage="constant")
model_true.G = G_true
model_true.tau_b = tau_b_true
model_true.eta_s = eta_s_true
G_star_true = model_true.predict(omega_true, test_mode='oscillation')

# Add noise
np.random.seed(42)
noise_factor = 1 + 0.1 * np.random.randn(len(omega_true))
G_star_data = G_star_true * noise_factor

# Step 1: NLSQ fit for initial guess
model_nlsq = TNTSingleMode(breakage="constant")
model_nlsq.fit(omega_true, G_star_data, test_mode='oscillation')
print(f"NLSQ: G={model_nlsq.G:.0f}, tau_b={model_nlsq.tau_b:.3f}, eta_s={model_nlsq.eta_s:.3f}")

# Step 2: Bayesian inference with NLSQ warm-start
result = model_nlsq.fit_bayesian(omega_true, G_star_data, test_mode='oscillation',
                                  num_warmup=1000, num_samples=2000, num_chains=4, seed=42)

# Step 3: Check diagnostics
print(f"R-hat (G): {result.diagnostics['r_hat']['G']:.4f}")  # Should be < 1.01
print(f"ESS (G): {result.diagnostics['ess']['G']:.0f}")      # Should be > 400

# Step 4: Credible intervals
intervals = model_nlsq.get_credible_intervals(result.posterior_samples, credibility=0.95)
print(f"G: [{intervals['G'][0]:.0f}, {intervals['G'][1]:.0f}] Pa (true: {G_true})")
print(f"tau_b: [{intervals['tau_b'][0]:.3f}, {intervals['tau_b'][1]:.3f}] s (true: {tau_b_true})")
print(f"eta_s: [{intervals['eta_s'][0]:.3f}, {intervals['eta_s'][1]:.3f}] Pa·s (true: {eta_s_true})")

# Step 5: Posterior predictive check
# (Sample from posterior and compare predictions to data)
posterior_samples = result.posterior_samples
n_samples = min(100, len(posterior_samples['G']))

plt.figure(figsize=(10, 6))
for i in range(n_samples):
    model_sample = TNTSingleMode(breakage="constant")
    model_sample.G = posterior_samples['G'][i]
    model_sample.tau_b = posterior_samples['tau_b'][i]
    model_sample.eta_s = posterior_samples['eta_s'][i]
    G_star_sample = model_sample.predict(omega_true, test_mode='oscillation')
    plt.loglog(omega_true, np.abs(G_star_sample), 'gray', alpha=0.1)

plt.loglog(omega_true, np.abs(G_star_data), 'ro', label='Data')
plt.loglog(omega_true, np.abs(G_star_true), 'b-', linewidth=2, label='True')
plt.xlabel('Angular frequency (rad/s)')
plt.ylabel('|G*| (Pa)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.title('Posterior Predictive Check (100 samples)')
plt.show()

Using BayesianPipeline for streamlined workflow:

from rheojax.pipeline.bayesian import BayesianPipeline

# Streamlined workflow with diagnostics and plotting
pipeline = BayesianPipeline()
(pipeline.load_data(omega_true, G_star_data, test_mode='oscillation')
         .fit_nlsq('tnt_tanaka_edwards')  # Automatic model selection
         .fit_bayesian(num_warmup=1000, num_samples=2000, num_chains=4)
         .plot_trace()      # MCMC trace plots
         .plot_pair()       # Parameter correlations
         .plot_forest()     # Credible intervals
         .save('tanaka_edwards_results.hdf5'))

See Also

Related TNT model variants:

Alternative constitutive models:

  • Giesekus Nonlinear Viscoelastic Models — Giesekus model family (nonlinear damping)

  • /models/ptt/index — Phan-Thien-Tanner model (shear thinning via trace of stress)

  • /models/rolie_poly/index — Rolie-Poly model (entangled polymer melts)

Theoretical background:

API Reference

class rheojax.models.TNTSingleMode(breakage='constant', stress_type='linear', xi=0.0)[source]

Single-mode Transient Network Theory model.

Implements the Green-Tobolsky / Tanaka-Edwards transient network model with composable physics variants. The conformation tensor S tracks the average chain configuration between reversible crosslinks.

The constitutive equation is:

dS/dt = L·S + S·L^T + g₀·I - β(S)·S

Stress is computed from S via σ = G·f(S)·(S - I) + η_s·γ̇.

Parameters:

  • breakage (Literal['constant', 'bell', 'power_law', 'stretch_creation']) – Breakage rate function: - “constant”: β = 1/τ_b (Tanaka-Edwards, UCM-like) - “bell”: β = (1/τ_b)·exp(ν·(stretch-1)) (force-dependent) - “power_law”: β = (1/τ_b)·stretch^m - “stretch_creation”: β = (1/τ_b), g₀ = (1+κ·stretch)/τ_b

  • stress_type (Literal['linear', 'fene']) – Stress formula: - “linear”: σ = G·(S - I) (Gaussian chains) - “fene”: σ = G·f(tr(S))·(S - I) (finitely extensible)

  • xi (float) – Gordon-Schowalter slip parameter (0=upper-convected, 1=corotational)

parameters

Model parameters

Type:

ParameterSet

fitted_

Whether the model has been fitted

Type:

bool

Examples

Basic Tanaka-Edwards model:

>>> model = TNTSingleMode()
>>> gamma_dot = np.logspace(-2, 2, 50)
>>> sigma = model.predict(gamma_dot, test_mode='flow_curve')

Bell force-dependent breakage:

>>> model = TNTSingleMode(breakage="bell")
>>> # Now has additional parameter 'nu' (force sensitivity)

See also

TNTLoopBridge

Two-species loop-bridge kinetics

TNTCates

Living polymer (wormlike micelle) model

__init__(breakage='constant', stress_type='linear', xi=0.0)[source]

Initialize single-mode TNT model.

Parameters:

  • breakage (Literal['constant', 'bell', 'power_law', 'stretch_creation']) – Breakage rate function type

  • stress_type (Literal['linear', 'fene']) – Stress formula type

  • xi (float) – Slip parameter for Gordon-Schowalter derivative

property G: float

Get network modulus G (Pa).

property tau_b: float

Get bond lifetime τ_b (s).

property eta_s: float

Get solvent viscosity η_s (Pa·s).

property eta_0: float

Get zero-shear viscosity η₀ = G·τ_b + η_s (Pa·s).

property breakage: str

Get breakage type.

property stress_type: str

Get stress type.

property xi: float

Get slip parameter ξ.

model_function(X, params, test_mode=None, **kwargs)[source]

NumPyro/BayesianMixin model function.

Routes to appropriate prediction based on test_mode. This is the stateless function used for both NLSQ optimization and Bayesian inference.

Parameters:

  • X (array-like) – Independent variable

  • params (array-like) – Parameter values in ParameterSet order: [G, tau_b, eta_s, …]

  • test_mode (str, optional) – Override stored test mode

  • **kwargs – Protocol-specific parameters: gamma_dot, sigma_applied, gamma_0, omega

Returns:

Predicted response

Return type:

jnp.ndarray

predict_flow_curve(gamma_dot, return_components=False)[source]

Predict steady shear stress and viscosity.

For constant breakage: analytical (UCM-like, no shear thinning). For non-constant breakage: ODE-to-steady-state (shear thinning).

Parameters:

  • gamma_dot (ndarray) – Shear rate array (1/s)

  • return_components (bool) – If True, return (sigma, eta, N1)

Returns:

Shear stress σ (Pa), or (σ, η, N₁) if return_components=True

Return type:

ndarray | tuple[ndarray, ndarray, ndarray]

predict_saos(omega, return_components=True)[source]

Predict SAOS storage and loss moduli.

In the linear regime, TNT reduces to single-mode Maxwell: G’(ω) = G·(ωτ_b)²/(1+(ωτ_b)²) G’’(ω) = G·(ωτ_b)/(1+(ωτ_b)²) + η_s·ω

Parameters:

  • omega (ndarray) – Angular frequency array (rad/s)

  • return_components (bool) – If True, return (G’, G’’)

Returns:

(G’, G’’) if return_components=True, else |G*|

Return type:

tuple[ndarray, ndarray] | ndarray

predict_normal_stresses(gamma_dot)[source]

Predict first and second normal stress differences.

Basic TNT (constant/linear/UC): N₁ = 2G·(τ_b·γ̇)², N₂ = 0. Gordon-Schowalter (ξ > 0): N₂ ≠ 0. FENE-P: N₁ enhanced by Peterlin factor f(trS).

Parameters:

gamma_dot (ndarray) – Shear rate array (1/s)

Returns:

(N₁, N₂) in Pa

Return type:

tuple[ndarray, ndarray]

simulate_startup(t, gamma_dot, return_full=False)[source]

Simulate startup flow at constant shear rate.

Parameters:

  • t (ndarray) – Time array (s)

  • gamma_dot (float) – Applied shear rate (1/s)

  • return_full (bool) – If True, return full conformation tensor components

Returns:

Shear stress σ(t), or (S_xx, S_yy, S_xy, S_zz) if return_full

Return type:

ndarray | tuple[ndarray, ndarray, ndarray, ndarray]

simulate_relaxation(t, gamma_dot_preshear, return_full=False)[source]

Simulate stress relaxation after cessation of steady shear.

For constant breakage + linear stress, relaxation is analytical: σ(t) = G·S_xy(0)·exp(-t/τ_b). For non-constant breakage or FENE stress, ODE integration is used.

Parameters:

  • t (ndarray) – Time array (s), starting from t=0 (cessation)

  • gamma_dot_preshear (float) – Shear rate before cessation (1/s)

  • return_full (bool) – If True, return full conformation tensor components

Returns:

Relaxing stress σ(t), or (S_xx, S_yy, S_xy, S_zz) if return_full

Return type:

ndarray | tuple[ndarray, ndarray, ndarray, ndarray]

simulate_creep(t, sigma_applied, return_rate=False)[source]

Simulate creep deformation under constant stress.

Parameters:

  • t (ndarray) – Time array (s)

  • sigma_applied (float) – Applied constant stress (Pa)

  • return_rate (bool) – If True, also return shear rate γ̇(t)

Returns:

Strain γ(t), or (γ, γ̇) if return_rate=True

Return type:

ndarray | tuple[ndarray, ndarray]

simulate_laos(t, gamma_0, omega, n_cycles=None)[source]

Simulate Large-Amplitude Oscillatory Shear (LAOS).

Parameters:

  • t (ndarray) – Time array (s), or None to auto-generate

  • gamma_0 (float) – Strain amplitude (dimensionless)

  • omega (float) – Angular frequency (rad/s)

  • n_cycles (int | None) – Number of oscillation cycles (overrides t)

Returns:

Dictionary with keys: ‘t’, ‘strain’, ‘stress’, ‘strain_rate’

Return type:

dict[str, ndarray]

get_overshoot_ratio(gamma_dot, t_max=None)[source]

Compute stress overshoot ratio in startup flow.

For constant breakage, there is no overshoot (UCM behavior). Overshoot requires Bell or stretch-dependent breakage.

Parameters:

  • gamma_dot (float) – Shear rate (1/s)

  • t_max (float | None) – Maximum simulation time (default: 20·τ_b)

Returns:

(overshoot_ratio, strain_at_overshoot)

Return type:

tuple[float, float]

get_relaxation_spectrum(t=None, n_points=100)[source]

Get relaxation modulus G(t).

For single-mode TNT: G(t) = G·exp(-t/τ_b)

Parameters:

  • t (ndarray | None) – Time array (default: logspace from 0.01·τ_b to 100·τ_b)

  • n_points (int) – Number of points if t not provided

Returns:

(t, G(t))

Return type:

tuple[ndarray, ndarray]

extract_laos_harmonics(laos_result, n_harmonics=5)[source]

Extract Fourier harmonics from LAOS stress response.

Parameters:

  • laos_result (dict[str, ndarray]) – Result from simulate_laos()

  • n_harmonics (int) – Number of harmonics to extract

Returns:

Dictionary with ‘n’, ‘sigma_prime’, ‘sigma_double_prime’, ‘intensity’, ‘I3_I1’

Return type:

dict[str, ndarray]

__repr__()

Return string representation.

Return type:

str

static compute_stretch(S_xx, S_yy, S_zz)

Compute stretch ratio from conformation tensor.

stretch = sqrt(tr(S)/3)

At equilibrium (S=I): stretch = 1. stretch > 1 indicates chains are extended beyond equilibrium.

Parameters:

  • S_xx (float) – Diagonal conformation tensor components

  • S_yy (float) – Diagonal conformation tensor components

  • S_zz (float) – Diagonal conformation tensor components

Returns:

Stretch ratio (dimensionless, ≥ 0)

Return type:

float

deborah_number(omega)

Compute Deborah number De = τ_b·ω.

Parameters:

omega (float) – Angular frequency (rad/s)

Returns:

Deborah number (dimensionless)

Return type:

float

fit(X, y, method='nlsq', check_compatibility=False, use_log_residuals=None, use_multi_start=None, n_starts=5, perturb_factor=0.3, deformation_mode=None, poisson_ratio=0.5, auto_init=False, return_result=False, check_physics=False, uncertainty=None, **kwargs)

Fit the model to data using NLSQ optimization.

This method uses NLSQ (GPU-accelerated nonlinear least squares) by default for fast point estimation. The optimization result is stored for potential warm-starting of Bayesian inference.

For very wide frequency ranges (>10 decades), multi-start optimization is automatically enabled to escape local minima.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – Target values

  • method (str) – Optimization method (‘nlsq’ by default for compatibility)

  • check_compatibility (bool) – Whether to check model-data compatibility before fitting. If True, warns when model may not be appropriate for data. Default is False for backward compatibility.

  • use_log_residuals (bool | None) – Whether to use log-space residuals for fitting. Recommended for wide frequency ranges (>8 decades) to prevent optimizer bias. If None (default), automatically detected based on data range. Explicit True/False overrides auto-detection.

  • use_multi_start (bool | None) – Whether to use multi-start optimization to escape local minima. Recommended for very wide ranges (>10 decades). If None (default), automatically enabled for >10 decades.

  • n_starts (int) – Number of random starts for multi-start optimization (default: 5)

  • perturb_factor (float) – Perturbation magnitude for multi-start random starts (default: 0.3). Parameters are perturbed by ± perturb_factor * (value or range). Larger values (0.7-0.9) explore wider parameter space.

  • auto_init (bool) – If True, calls auto_p0() to estimate initial parameters from data before running the optimizer (default: False).

  • return_result (bool) – If True, returns a FitResult instead of self. This intentionally breaks method chaining for workflows that need structured result objects (default: False).

  • check_physics (bool) – If True, runs post-fit physics validation and emits RheoJaxPhysicsWarning for any violations (default: False).

  • uncertainty (str | None) – Post-fit uncertainty method. "hessian" for fast Cramér-Rao bounds, "bootstrap" for residual bootstrap CIs, or None to skip (default: None).

  • **kwargs – Additional fitting options passed to _fit()

Return type:

BaseModel | Any

Returns:

self for method chaining (default), or FitResult if return_result=True.

Example

>>> model = Maxwell()
>>> model.fit(t, G_data)  # Uses NLSQ by default
>>> model.fit(t, G_data, method='nlsq', max_iter=1000)
>>> model.fit(t, G_data, check_compatibility=True)  # Check compatibility
>>> model.fit(omega, G_star, use_log_residuals=True)  # Force log-residuals
>>> model.fit(mastercurve, None, use_multi_start=True, n_starts=10)  # Multi-start
>>> result = model.fit(t, G_data, return_result=True)  # Structured result
>>> result = model.fit(t, G_data, auto_init=True, check_physics=True,
...                    return_result=True)  # Full pipeline
fit_bayesian(X, y=None, num_warmup=1000, num_samples=2000, num_chains=4, initial_values=None, test_mode=None, seed=None, deformation_mode=None, poisson_ratio=0.5, **nuts_kwargs)

Perform Bayesian inference using NumPyro NUTS sampler.

This method delegates to BayesianMixin.fit_bayesian() to run NUTS sampling for Bayesian parameter estimation. If initial_values is not provided and the model has been previously fitted with fit(), the NLSQ point estimates are automatically used for warm-starting.

Multi-chain sampling is enabled by default (num_chains=4) to provide reliable convergence diagnostics (R-hat, ESS) and parallel execution on multi-GPU systems.

Parameters:

  • X (numpy.typing.ArrayLike) – Independent variable data (input features) or RheoData object

  • y (Optional[numpy.typing.ArrayLike]) – Dependent variable data (observations to fit). If X is RheoData, y is ignored and extracted from X.

  • num_warmup (int) – Number of warmup/burn-in iterations (default: 1000)

  • num_samples (int) – Number of posterior samples per chain (default: 2000)

  • num_chains (int) – Number of MCMC chains (default: 4). Multiple chains enable proper R-hat computation and parallel execution. Chain method is auto-selected: ‘parallel’ on multi-GPU, ‘vectorized’ on single GPU/CPU.

  • initial_values (dict[str, float] | None) – Optional dict of initial parameter values for warm-start. If None and model is fitted, uses NLSQ estimates.

  • test_mode (str | None) – Explicit test mode (e.g., ‘relaxation’, ‘creep’, ‘oscillation’). If None, inferred from RheoData.metadata[‘test_mode’] or defaults to ‘relaxation’. Overrides RheoData metadata if provided.

  • seed (int | None) – Random seed for reproducibility. If None, uses seed=0 for deterministic results. Set to different values for independent runs.

  • **nuts_kwargs – Additional arguments passed to NUTS sampler (e.g., target_accept_prob, chain_method)

Return type:

BayesianResult

Returns:

BayesianResult containing posterior samples, summary statistics, and convergence diagnostics (R-hat, ESS, divergences)

Example

>>> model = Maxwell()
>>> # Warm-start from NLSQ with explicit mode
>>> model.fit(t, G_data, test_mode='relaxation')  # NLSQ optimization
>>> result = model.fit_bayesian(t, G_data, test_mode='relaxation')
>>>
>>> # RheoData with embedded mode (recommended)
>>> rheo_data = RheoData(x=omega, y=G_star, metadata={'test_mode': 'oscillation'})
>>> result = model.fit_bayesian(rheo_data)
>>>
>>> # Or provide explicit initial values
>>> result = model.fit_bayesian(
...     t, G_data,
...     initial_values={'G0': 1e5, 'eta': 1e3},
...     test_mode='creep'
... )
fit_predict(X, y, **kwargs)

Fit model and return predictions.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – Target values

  • **kwargs – Additional fitting options

Return type:

numpy.typing.ArrayLike

Returns:

Model predictions on training data

classmethod from_dict(data)

Create model from dictionary.

Parameters:

data (dict[str, Any]) – Dictionary representation

Return type:

BaseModel

Returns:

Model instance

get_bayesian_result()

Get stored Bayesian inference result.

Return type:

BayesianResult | None

Returns:

BayesianResult from fit_bayesian(), or None if not run

Example

>>> model.fit_bayesian(t, G_data)
>>> result = model.get_bayesian_result()
>>> print(result.diagnostics['r_hat'])
get_credible_intervals(posterior_samples, credibility=0.95)

Compute highest density intervals (HDI) for posterior samples.

Computes the highest posterior density interval for each parameter, which is the shortest interval containing the specified probability mass. This is preferred over equal-tailed intervals for skewed distributions.

Parameters:

  • posterior_samples (dict[str, ndarray]) – Dictionary mapping parameter names to posterior sample arrays (from BayesianResult.posterior_samples)

  • credibility (float) – Probability mass to include in interval (default: 0.95) Common values: 0.68 (1 sigma), 0.95 (2 sigma), 0.997 (3 sigma)

Return type:

dict[str, tuple[float, float]]

Returns:

Dictionary mapping parameter names to (lower, upper) credible interval tuples. All values are float64.

Example

>>> result = model.fit_bayesian(X, y)
>>> intervals_95 = model.get_credible_intervals(
...     result.posterior_samples, credibility=0.95
... )
>>> print(f"95% CI for a: {intervals_95['a']}")
static get_equilibrium_conformation()

Return equilibrium conformation tensor S_eq = I.

In the absence of flow, chains adopt their equilibrium configuration: S = I (isotropic, unstretched).

Returns:

Equilibrium state [S_xx, S_yy, S_zz, S_xy] = [1, 1, 1, 0]

Return type:

Array

get_nlsq_result()

Get stored NLSQ optimization result.

Returns:

OptimizationResult from NLSQ fit, or None if not fitted

Example

>>> model.fit(t, G_data)
>>> result = model.get_nlsq_result()
>>> if result:
...     print(f"Converged: {result.success}")
get_parameter_dict()

Get all parameters as a dictionary.

Returns:

Dictionary of parameter name -> value

Return type:

dict[str, float]

get_parameter_uncertainties()

Get standard errors for fitted parameters from NLSQ covariance.

Returns:

std_error}, or None if covariance unavailable

Return type:

dict of {param_name

get_params(deep=True)

Get model parameters.

Parameters:

deep (bool) – If True, return parameters of sub-objects

Return type:

dict[str, Any]

Returns:

Dictionary of parameter names and values

initialize_from_flow_curve(gamma_dot, sigma)

Initialize parameters from flow curve data.

Parameters:

  • gamma_dot (ndarray) – Shear rate array (1/s)

  • sigma (ndarray) – Stress array (Pa)

Return type:

None

initialize_from_saos(omega, G_prime, G_double_prime)

Initialize parameters from SAOS data.

Uses the crossover frequency and modulus to estimate τ_b and G. In the linear regime, TNT reduces to Maxwell: crossover at ω_c·τ_b = 1.

Parameters:

  • omega (ndarray) – Angular frequency array (rad/s)

  • G_prime (ndarray) – Storage modulus G’ (Pa)

  • G_double_prime (ndarray) – Loss modulus G’’ (Pa)

Return type:

None

property pcov_

Parameter covariance matrix from NLSQ fit.

Returns:

ndarray of shape (n_params, n_params), or None if not fitted

property popt_

Optimal parameter values from NLSQ fit.

Returns:

ndarray of shape (n_params,), or None if not fitted

precompile(test_mode='relaxation', X=None, y=None)

Precompile NLSQ residual functions to eliminate JIT cold-start.

Triggers JIT compilation by running a minimal fit (max_iter=1) with dummy data. The model parameters are reset afterwards so the model is left in its original state.

This is useful for interactive sessions or benchmarks where the ~870ms first-fit JIT overhead should be excluded.

Parameters:

  • test_mode (str) – Test mode to precompile for (default: ‘relaxation’).

  • X (Optional[numpy.typing.ArrayLike]) – Optional input data for shape inference. If None, uses a 10-point logspace array.

  • y (Optional[numpy.typing.ArrayLike]) – Optional output data. If None, generates ones matching X.

Return type:

float

Returns:

Compilation time in seconds.

Example

>>> model = Maxwell()
>>> t = model.precompile(test_mode='relaxation')
>>> print(f"Compiled in {t:.2f}s")
>>> model.fit(X, y)  # No JIT overhead
precompile_bayesian(X=None, y=None, test_mode=None, num_chains=4)

Precompile NUTS kernel to eliminate JIT overhead in subsequent calls.

Triggers JIT compilation of the NumPyro model by running a minimal sampling (1 warmup, 1 sample). This caches the compiled kernel so that subsequent fit_bayesian() calls are 2-5x faster.

Parameters:

  • X (ndarray | RheoData | None) – Sample input data for determining array shapes. If None, uses a default 10-point linspace [0.01, 100].

  • y (ndarray | None) – Sample output data. If None, generates dummy data.

  • test_mode (str | TestModeEnum | None) – Test mode to precompile for. If None, defaults to ‘relaxation’.

Returns:

Compilation time in seconds.

Return type:

float

Example

>>> model = Maxwell()
>>> compile_time = model.precompile_bayesian(test_mode='relaxation')
>>> print(f"Compiled in {compile_time:.1f}s")
>>> # Now fit_bayesian() will be faster
>>> result = model.fit_bayesian(X, y)  # No compilation overhead
predict(X, test_mode=None, deformation_mode=None, poisson_ratio=None, **kwargs)

Make predictions.

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • test_mode (str | None) – Optional test mode (‘oscillation’, ‘relaxation’, ‘creep’, ‘flow’). If provided, sets model’s test_mode before prediction. Useful for data generation without fitting.

  • deformation_mode (str | DeformationMode | None) – Optional deformation mode for output conversion. If None, uses the mode stored from fit(). If tensile, converts G* predictions to E* space.

  • poisson_ratio (float | None) – Poisson’s ratio for conversion. If None, uses value stored from fit() (default 0.5).

  • **kwargs – Additional arguments passed to the internal _predict method.

Return type:

numpy.typing.ArrayLike

Returns:

Model predictions (in E* space if deformation_mode is tensile)

sample_prior(num_samples=1000, seed=None)

Sample from prior distributions over parameter bounds.

Samples from uniform prior distributions defined by parameter bounds. This is useful for prior predictive checks and understanding the prior’s influence on the posterior.

Parameters:

  • num_samples (int) – Number of samples to draw from prior (default: 1000)

  • seed (int | None) – Random seed for reproducibility (default: None)

Return type:

dict[str, ndarray]

Returns:

Dictionary mapping parameter names to arrays of prior samples. Each array has shape (num_samples,) and dtype float64.

Raises:

Example

>>> model = MyModel()
>>> prior_samples = model.sample_prior(num_samples=500, seed=42)
>>> print(prior_samples["a"].shape)  # (500,)
score(X, y)

Compute model score (R² by default).

Parameters:

  • X (numpy.typing.ArrayLike) – Input features

  • y (numpy.typing.ArrayLike) – True target values

Return type:

float

Returns:

Model score (R² coefficient)

set_parameter_dict(params)

Set parameters from a dictionary.

Parameters:

params (dict[str, float]) – Dictionary of parameter name -> value

Return type:

None

set_params(**params)

Set model parameters.

Parameters:

**params – Parameter names and values

Return type:

BaseModel

Returns:

self for method chaining

to_dict()

Serialize model to dictionary.

Return type:

dict[str, Any]

Returns:

Dictionary representation of model

weissenberg_number(gamma_dot)

Compute Weissenberg number Wi = τ_b·γ̇.

Parameters:

gamma_dot (float) – Shear rate (1/s)

Returns:

Weissenberg number (dimensionless)

Return type:

float

Key methods:

  • fit(x, y, test_mode, **kwargs): NLSQ fit to data

  • predict(x, test_mode, **kwargs): Model prediction

  • fit_bayesian(x, y, test_mode, **kwargs): Bayesian inference with NUTS

  • get_credible_intervals(samples, credibility): Extract credible intervals from posterior

Parameters (attributes):

  • G: Network elastic modulus (Pa), default 1000, bounds (1, 1e8)

  • tau_b: Bond lifetime (s), default 1.0, bounds (1e-6, 1e4)

  • eta_s: Solvent viscosity (Pa·s), default 0.0, bounds (0, 1e4)

Test modes:

  • 'oscillation': Small-amplitude oscillatory shear (SAOS)

  • 'relaxation': Stress relaxation after step strain

  • 'startup': Startup of steady shear

  • 'flow_curve': Steady-state flow curve

  • 'creep': Creep under constant stress

  • 'laos': Large-amplitude oscillatory shear

References