DMT Thixotropic Models¶

Quick Reference¶

Model Summary¶
Model Classes	`DMTLocal`, `DMTNonlocal`
Physics	de Souza Mendes-Thompson thixotropic viscoelasticity
Viscosity Closures	`"exponential"`, `"herschel_bulkley"`
Elasticity	Optional Maxwell backbone (`include_elasticity=True`)
Protocols	FLOW_CURVE, CREEP, RELAXATION, STARTUP, OSCILLATION, LAOS
Key Features	Structure kinetics, stress overshoot, delayed yielding, shear banding

Import:

from rheojax.models import DMTLocal, DMTNonlocal

Basic Usage:

# Exponential closure with Maxwell elasticity
model = DMTLocal(closure="exponential", include_elasticity=True)

# Herschel-Bulkley closure for yield-stress materials
model = DMTLocal(closure="herschel_bulkley", include_elasticity=True)

# Nonlocal variant for shear banding
model = DMTNonlocal(closure="exponential", n_points=51, gap_width=1e-3)

Notation Guide¶

Symbol	Description	Units
\(\lambda\)	Structure parameter (0 = broken, 1 = fully structured)	dimensionless
\(\eta_0\)	Zero-shear viscosity (at \(\lambda = 1\))	Pa·s
\(\eta_\infty\)	Infinite-shear viscosity (at \(\lambda = 0\))	Pa·s
\(\tau_y\)	Yield stress	Pa
\(K\)	Consistency index	Pa·sⁿ
\(n\)	Flow index	dimensionless
\(G\)	Elastic modulus	Pa
\(\theta\)	Relaxation time (\(= \eta/G\))	s
\(t_{eq}\)	Equilibrium (buildup) timescale	s
\(a\)	Breakdown rate coefficient	dimensionless
\(c\)	Breakdown rate exponent	dimensionless

Overview¶

The de Souza Mendes-Thompson (DMT) model [deSouzaMendes2009] [Mendes2012] is a structural-kinetics based thixotropic model that captures time-dependent rheological behavior through a scalar structure parameter \(\lambda \in [0, 1]\).

Key Features¶

Structure-dependent viscosity: Material properties depend on microstructural state tracked by \(\lambda\), with fully structured (\(\lambda = 1\)) giving high viscosity and fully broken (\(\lambda = 0\)) giving low viscosity.
Structure kinetics: The structure evolves through competing buildup (aging at rest) and breakdown (shear-induced destruction) processes.
Multiple viscosity closures: Either smooth exponential dependence or Herschel-Bulkley form with explicit yield stress.
Optional viscoelasticity: Maxwell backbone enables stress overshoot in startup and elastic recoil.
Spatial extension: Nonlocal variant captures shear banding through structure diffusion.

Historical Context¶

The DMT model represents a sophisticated synthesis of several theoretical traditions in thixotropic rheology. Understanding this lineage helps clarify the model’s capabilities and limitations.

Jeffreys/Oldroyd-B Backbone¶

When include_elasticity=True, the DMT model uses a structure-dependent Maxwell element in parallel with a Newtonian element. This yields the classical Jeffreys (three-element) viscoelastic form:

\[\tau_1(\lambda)\dot{\sigma} + \sigma = \eta(\lambda)\left(\dot{\gamma} + \tau_2(\lambda)\ddot{\gamma}\right)\]

where the relaxation time \(\tau_1(\lambda)\) and retardation time \(\tau_2(\lambda)\) depend on the structure parameter. This formulation is numerically equivalent to the Maxwell-in-parallel representation used in RheoJAX:

\[\sigma(t) = \sigma_M(t) + \eta_s(\lambda)\dot{\gamma}(t)\]

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda)} = G(\lambda)\dot{\gamma}\]

This form avoids \(\ddot{\gamma}\) and is ideal for protocol replay and fitting.

Evolution from Simple Thixotropy¶

The DMT model evolved from simpler structural-kinetics models:

Coussot model (early 2000s): Introduced the avalanche effect and viscosity bifurcation for yield-stress fluids. Pure viscosity-structure coupling without elasticity.
Houska model: Added Bingham-like yield stress with separate thixotropic contributions to both yield stress and consistency.
Mujumdar model (2002): Introduced elastic strain as a state variable with a yield function, enabling stress overshoot predictions.
DMT model (2009-2014): Unified these elements into a comprehensive framework with explicit yield stress, optional Maxwell viscoelasticity, and structure-dependent material functions.

Physical Foundations¶

Structure Parameter¶

The structure parameter \(\lambda \in [0, 1]\) represents the degree of microstructural organization:

\(\lambda = 1\): Fully structured (at rest, aged)
\(\lambda = 0\): Fully broken (high shear, rejuvenated)

Physical interpretation includes:

Colloidal networks: Bond connectivity between particles
Polymer solutions: Entanglement density
Emulsions/foams: Droplet/bubble deformation state

Structure Kinetics¶

The structure evolves according to:

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}} - \frac{a \lambda |\dot{\gamma}|^c}{t_{eq}}\]

where:

First term: Buildup (aging) drives \(\lambda \to 1\) at rate \(1/t_{eq}\)
Second term: Breakdown destroys structure at rate proportional to \(|\dot{\gamma}|^c\)

At equilibrium (\(d\lambda/dt = 0\)):

\[\lambda_{eq} = \frac{1}{1 + a|\dot{\gamma}|^c}\]

This gives \(\lambda_{eq} \to 1\) as \(\dot{\gamma} \to 0\) and \(\lambda_{eq} \to 0\) as \(\dot{\gamma} \to \infty\).

Fluidity Interpretation¶

The structure parameter \(\lambda\) has an alternative interpretation in terms of fluidity \(\phi = 1/\eta\), the inverse of viscosity:

High \(\lambda\) (structured) → low \(\phi\) (viscous, solid-like)
Low \(\lambda\) (broken) → high \(\phi\) (fluid, liquid-like)

Fluidity-based kinetics provide an equivalent formulation:

\[\frac{d\phi}{dt} = -\frac{\phi - \phi_{min}}{t_{age}} + c|\dot{\gamma}|^m(\phi_{max} - \phi)\]

This is mathematically equivalent to the \(\lambda\)-form after a monotone transformation. The fluidity approach is particularly natural for:

Spatial extensions (nonlocal fluidity models)
Connection to soft glassy rheology (SGR)
Shear banding analysis

Cooperativity length for nonlocal models:

\[\xi \sim \sqrt{D_\lambda \cdot t_{eq}}\]

where \(D_\lambda\) is the structure diffusion coefficient. This length scale sets the minimum shear band width and prevents ill-posed sharp band interfaces.

Viscosity Closures¶

Exponential Closure¶

A smooth, monotonic relationship between structure and viscosity:

\[\eta(\lambda) = \eta_\infty \left(\frac{\eta_0}{\eta_\infty}\right)^\lambda\]

Properties:

\(\eta(1) = \eta_0\) (zero-shear viscosity)
\(\eta(0) = \eta_\infty\) (infinite-shear viscosity)
No explicit yield stress (power-law-like flow curve)

Herschel-Bulkley Closure¶

Structure-dependent yield stress and consistency:

\[\sigma = \tau_y(\lambda) + K(\lambda) |\dot{\gamma}|^n + \eta_\infty \dot{\gamma}\]

where:

\[\begin{split}\tau_y(\lambda) &= \tau_{y0} \lambda^{m_1} \\ K(\lambda) &= K_0 \lambda^{m_2}\end{split}\]

Properties:

Explicit yield stress \(\tau_y\) controlled by \(\lambda\)
True yield stress behavior (regularized with Papanastasiou)
Structure-dependent flow index contribution

Maxwell Viscoelasticity¶

When include_elasticity=True, a Maxwell element adds elastic response:

\[\frac{d\sigma}{dt} = G(\lambda) \dot{\gamma} - \frac{G(\lambda)}{\eta(\lambda)} \sigma\]

where the elastic modulus depends on structure:

\[G(\lambda) = G_0 \lambda^{m_G}\]

This gives:

Relaxation time: \(\theta(\lambda) = \eta(\lambda) / G(\lambda)\)
Stress overshoot: In startup, stress overshoots before reaching steady state
Stress relaxation: After cessation of flow, stress decays exponentially
SAOS: Storage (\(G'\)) and loss (\(G''\)) moduli from linear response

Steady-State Flow Curve¶

At equilibrium, the structure and stress are uniquely determined by shear rate.

Exponential Closure¶

\[\sigma_{ss}(\dot{\gamma}) = \eta(\lambda_{eq}(\dot{\gamma})) \cdot \dot{\gamma}\]

where \(\eta\) depends on \(\lambda_{eq} = 1/(1 + a|\dot{\gamma}|^c)\).

Herschel-Bulkley Closure¶

\[\sigma_{ss} = \tau_{y0}\lambda_{eq}^{m_1} + K_0\lambda_{eq}^{m_2}|\dot{\gamma}|^n + \eta_\infty\dot{\gamma}\]

Viscosity Bifurcation¶

A defining feature of yield-stress materials is the viscosity bifurcation: a discontinuous transition between solid-like and liquid-like behavior at a critical stress.

Critical Stress¶

For the Herschel-Bulkley closure, the viscosity bifurcation occurs at the yield stress:

\(\sigma < \tau_y(\lambda)\): Effective viscosity diverges (\(\eta \to \infty\))
\(\sigma > \tau_y(\lambda)\): Finite viscosity, flow occurs

At equilibrium structure \(\lambda_{eq}\):

\[\sigma_c = \tau_{y0} \lambda_{eq}^{m_1} = \frac{\tau_{y0}}{(1 + a|\dot{\gamma}|^c)^{m_1}}\]

Avalanche Effect¶

Near the critical stress, thixotropic materials exhibit the avalanche effect [Coussot2002]:

Below yield: \(\sigma_0 < \tau_y\), structure builds up, viscosity increases
Near yield: \(\sigma_0 \approx \tau_y\), metastable state with slow creep
Delayed yielding: Structure slowly breaks down, then catastrophic flow onset
Above yield: Immediate structure breakdown and steady flow

This manifests in creep experiments as an initial slow strain accumulation followed by rapid acceleration when the structure can no longer support the applied stress.

Connection to Herschel-Bulkley¶

In the limit of fast thixotropic kinetics (\(t_{eq} \to 0\)), the DMT model with Herschel-Bulkley closure recovers the classical Herschel-Bulkley constitutive equation:

\[\sigma = \tau_y + K|\dot{\gamma}|^n \quad \text{for } \sigma > \tau_y\]

The key difference is that DMT predicts time-dependent transitions between solid and liquid states, while HB assumes instantaneous equilibrium.

Related Models Comparison¶

The DMT model belongs to a family of structural-kinetics thixotropic models. The table below compares key features:

Model	Yield Stress	Elasticity	Thixotropy	Params	Complexity	Best For
Coussot	Implicit	No	Yes	~4	Low	Avalanche effect, simple thixotropy
Houska	Explicit	No	Yes	~7	Medium	Drilling fluids, concrete
Mujumdar	Explicit	Yes	Yes	~6	Medium	Stress overshoot, transients
Dullaert-Mewis	Implicit	Partial	Yes	~8	Medium	Colloidal suspensions
DMT	Explicit	Optional	Yes	~10	High	Full protocol matching

Coussot Model¶

Simple thixotropy without explicit yield stress:

\[\sigma = \eta(\lambda)\dot{\gamma}, \quad \eta(\lambda) = \eta_0(1 + \lambda^n)\]

\[\frac{d\lambda}{dt} = \frac{1}{\theta_0} - \alpha\lambda|\dot{\gamma}|\]

Strengths: Simple, captures avalanche effect. Limitations: No elasticity, no stress overshoot.

Houska Model¶

Bingham-like with thixotropic contributions:

\[\sigma = (\tau_{y0} + \tau_{yt}\lambda) + (K_0 + K_t\lambda)|\dot{\gamma}|^n\]

Strengths: Explicit yield stress, separate thixotropic contributions. Limitations: No elasticity, many parameters.

Mujumdar Model¶

Elastic strain with yield function:

\[\sigma = G\gamma_e + \eta(\lambda)\dot{\gamma}\]

\[\frac{d\gamma_e}{dt} = \dot{\gamma} - k_{rel}\gamma_e f(\sigma)\]

where \(f(\sigma)\) is a yield function. Strengths: Stress overshoot, elastic recovery. Limitations: Step-function yield behavior.

When to Use Which Model¶

Situation	Recommended Model	Why
Quick thixotropy fitting	Coussot → DMT	Start simple, upgrade if needed
Need yield + thixotropy	Houska or DMT	Explicit yield stress
Stress overshoot important	DMT or Mujumdar	Requires elasticity
Creep with delayed yield	DMT	Best bifurcation physics
LAOS analysis	DMT (with elasticity)	Full nonlinear response
Full protocol matching	DMT	Most comprehensive

Parameters¶

Core Viscosity Parameters¶

Parameter	Description	Units	Default	Bounds
`eta_0`	Zero-shear viscosity	Pa·s	1e5	[1e2, 1e8]
`eta_inf`	Infinite-shear viscosity	Pa·s	0.1	[1e-3, 1e2]

Herschel-Bulkley Parameters (`closure="herschel_bulkley"` only)¶

Parameter	Description	Units	Default	Bounds
`tau_y0`	Fully-structured yield stress	Pa	10.0	[0.1, 1e4]
`K0`	Fully-structured consistency	Pa·sⁿ	5.0	[0.1, 1e3]
`n_flow`	Flow index	—	0.5	[0.1, 1.0]
`m1`	Yield stress exponent	—	1.0	[0.5, 2.0]
`m2`	Consistency exponent	—	1.0	[0.5, 2.0]

Elastic Parameters (`include_elasticity=True` only)¶

Parameter	Description	Units	Default	Bounds
`G0`	Elastic modulus at \(\lambda = 1\)	Pa	100.0	[1e0, 1e6]
`m_G`	Modulus structure exponent	—	1.0	[0.5, 2.0]

Structure Kinetics Parameters¶

Parameter	Description	Units	Default	Bounds
`t_eq`	Equilibrium (buildup) timescale	s	100.0	[0.1, 1e4]
`a`	Breakdown rate coefficient	—	1.0	[1e-3, 1e2]
`c`	Breakdown rate exponent	—	1.0	[0.1, 2.0]

Nonlocal Parameters (`DMTNonlocal` only)¶

Parameter	Description	Units	Default	Bounds
`D_lambda`	Structure diffusion coefficient	m^2/s	1e-9	[1e-12, 1e-6]

Dimensionless Groups¶

Several dimensionless numbers characterize DMT model behavior and help classify rheological regimes.

Deborah Number (De)¶

The ratio of thixotropic timescale to experimental timescale:

\[De = \frac{t_{eq}}{t_{exp}}\]

where \(t_{exp}\) is a characteristic experimental time (e.g., period \(2\pi/\omega\) for oscillation, \(1/\dot{\gamma}\) for steady shear).

\(De \gg 1\): Thixotropy dominates; structure cannot equilibrate during experiment
\(De \ll 1\): Quasi-steady behavior; structure rapidly equilibrates

Weissenberg Number (Wi)¶

Characterizes the structure breakdown rate:

\[Wi = \dot{\gamma} \cdot t_{eq}\]

\(Wi \gg 1\): Strong breakdown, \(\lambda \to 0\)
\(Wi \ll 1\): Structure preserved, \(\lambda \to 1\)

Structure Number (Sn)¶

The breakdown efficiency parameter:

\[Sn = a \cdot |\dot{\gamma}|^c\]

This appears directly in the equilibrium structure:

\[\lambda_{eq} = \frac{1}{1 + Sn}\]

Regime Classification¶

Regime	De	Wi	Behavior
Linear viscoelastic	Low	Low	Standard Maxwell, \(\lambda \approx 1\)
Thixotropic	High	Variable	Time-dependent, hysteresis loops
Shear-thinning	Low	High	Steady power-law, \(\lambda = \lambda_{eq}\)
Glass-like	High	High	Aging + strong breakdown, complex transients

These regimes map to different experimental signatures:

Linear VE: Standard \(G'\), \(G''\) from SAOS
Thixotropic: Hysteresis in flow curves, recovery experiments
Shear-thinning: Rate-dependent steady viscosity
Glass-like: Aging effects, bifurcation, delayed yielding

Protocol-Specific Equations¶

This section provides complete mathematical derivations for each rheological protocol. These equations form the basis for numerical simulations and analytical predictions.

General Governing System¶

For all protocols, the DMT model is governed by a coupled system of differential equations.

State variables:

\(\lambda(t)\): Structure parameter
\(\gamma_e(t)\): Elastic strain (when include_elasticity=True)
\(\sigma(t)\): Stress

Constitutive equation (generalized Newtonian form):

\[\sigma(t) = \eta(\lambda(t), \dot{\gamma}(t)) \cdot \dot{\gamma}(t)\]

Structure kinetics:

\[\frac{d\lambda}{dt} = \underbrace{\frac{1 - \lambda}{t_{eq}}}_{\text{buildup (aging)}} - \underbrace{\frac{a \lambda |\dot{\gamma}|^c}{t_{eq}}}_{\text{breakdown (shear)}}\]

Maxwell backbone (when include_elasticity=True):

\[\sigma(t) = \sigma_M(t) + \eta_s(\lambda) \dot{\gamma}(t)\]

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda)} = G(\lambda) \dot{\gamma}\]

where \(\tau_1(\lambda) = \eta_M(\lambda) / G(\lambda)\) is the structure-dependent relaxation time.

Rotation / Steady-State Flow Curve¶

Protocol: Constant shear rate \(\dot{\gamma} = \text{const}\), wait for equilibrium.

Equilibrium Condition¶

At steady state, \(d\lambda/dt = 0\):

\[0 = \frac{1 - \lambda_{eq}}{t_{eq}} - \frac{a \lambda_{eq} |\dot{\gamma}|^c}{t_{eq}}\]

Solving for equilibrium structure:

\[\boxed{\lambda_{eq}(\dot{\gamma}) = \frac{1}{1 + a|\dot{\gamma}|^c}}\]

Limiting behaviors:

\(\dot{\gamma} \to 0\): \(\lambda_{eq} \to 1\) (fully structured)
\(\dot{\gamma} \to \infty\): \(\lambda_{eq} \to 0\) (fully broken)

Steady-State Stress¶

Exponential closure:

\[\sigma_{ss}(\dot{\gamma}) = \eta(\lambda_{eq}) \cdot \dot{\gamma} = \eta_\infty \left(\frac{\eta_0}{\eta_\infty}\right)^{\lambda_{eq}} \dot{\gamma}\]

Herschel-Bulkley closure:

\[\sigma_{ss}(\dot{\gamma}) = \tau_{y0} \lambda_{eq}^{m_1} + K_0 \lambda_{eq}^{m_2} |\dot{\gamma}|^n + \eta_\infty \dot{\gamma}\]

Maxwell backbone (steady state, \(\dot{\sigma}_M = 0\)):

\[\sigma_M^{ss} = G(\lambda_{eq}) \tau_1(\lambda_{eq}) \dot{\gamma}\]

\[\sigma_{ss} = \sigma_M^{ss} + \eta_s(\lambda_{eq}) \dot{\gamma} = \left[ G(\lambda_{eq}) \tau_1(\lambda_{eq}) + \eta_s(\lambda_{eq}) \right] \dot{\gamma}\]

Controlled-Stress Mode¶

When the rheometer controls stress \(\sigma\) rather than rate \(\dot{\gamma}\), solve the implicit equation:

\[\sigma = \eta(\lambda_{ss}, \dot{\gamma}) \cdot \dot{\gamma}\]

together with:

\[\lambda_{ss} = \frac{1}{1 + a|\dot{\gamma}|^c}\]

This can yield multiple solutions (S-shaped flow curve), representing a route to shear banding in local models when the middle branch is unstable.

Start-up of Steady Shear¶

Protocol: \(\dot{\gamma}(t) = \dot{\gamma}_0 H(t)\) from rest (\(\lambda_0 = 1\)).

Governing ODEs¶

For \(t > 0\):

Structure evolution:

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}} - \frac{a \lambda |\dot{\gamma}_0|^c}{t_{eq}}\]

Closed-form solution: Let \(r = \frac{1}{t_{eq}} + \frac{a|\dot{\gamma}_0|^c}{t_{eq}}\). Then:

\[\lambda(t) = \lambda_{ss} + (\lambda_0 - \lambda_{ss}) e^{-rt}\]

where \(\lambda_{ss} = \frac{1}{1 + a|\dot{\gamma}_0|^c}\) is the target equilibrium.

Maxwell stress evolution (with elasticity):

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda(t))} = G(\lambda(t)) \dot{\gamma}_0\]

This ODE has time-varying coefficients through \(\lambda(t)\), requiring numerical integration.

Stress Overshoot Mechanism¶

The stress overshoot characteristic of thixotropic materials arises from the interplay of elastic loading and structure breakdown:

Initial elastic rise (\(t \ll t_{eq}\)): \(\sigma \approx G_0 \cdot \gamma = G_0 \dot{\gamma}_0 t\) (linear elastic loading while \(\lambda \approx 1\))
Structure breakdown begins: \(\lambda\) starts decreasing toward \(\lambda_{ss}\)
Viscosity drop: As \(\lambda\) decreases, \(\eta(\lambda)\) drops exponentially (exponential closure) or polynomially (HB closure)
Peak stress: Maximum occurs when the rate of viscosity decrease exceeds the rate of strain increase
Approach to steady state: \(\sigma \to \sigma_{ss}\), \(\lambda \to \lambda_{ss}\)

Peak strain estimate (order of magnitude):

\[\gamma_{peak} \sim \frac{\tau_y}{G_0} \quad \text{(yield strain)}\]

Overshoot ratio: \(\sigma_{peak}/\sigma_{ss}\) increases with shear rate and initial structure level \(\lambda_0\).

Purely Generalized Newtonian (No Elasticity)¶

Without elasticity (include_elasticity=False), the stress follows viscosity directly:

\[\sigma(t) = \eta(\lambda(t), \dot{\gamma}_0) \cdot \dot{\gamma}_0\]

In this case, there is no elastic stress overshoot. However, transient “viscosity overshoot/undershoot” can occur depending on how \(\eta(\lambda, \dot{\gamma})\) depends on both arguments. For dense glasses with true overshoots, the Maxwell backbone is required.

Stress Relaxation (Cessation of Flow)¶

Protocol: Shear at \(\dot{\gamma}_0\) until \(t = 0\), then \(\dot{\gamma}(t > 0) = 0\).

Structure Recovery¶

For \(t > 0\), with \(\dot{\gamma} = 0\):

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}}\]

(pure buildup, no breakdown). Solution:

\[\lambda(t) = 1 - (1 - \lambda_0) e^{-t/t_{eq}}\]

where \(\lambda_0 = \lambda_{ss}(\dot{\gamma}_0)\) is the structure at cessation.

As \(t \to \infty\), \(\lambda \to 1\) (full recovery).

Stress Relaxation (Generalized Newtonian)¶

Without elasticity:

\[\sigma(t > 0) = \eta(\lambda(t), 0) \cdot 0 = 0\]

There is no stress relaxation in the strict sense because there is no stored elastic stress. What is predicted is structural recovery (thixotropic rebuild).

Stress Relaxation (Maxwell Backbone)¶

With elasticity, the Maxwell stress relaxes:

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda(t))} = 0\]

Solution:

\[\sigma_M(t) = \sigma_M(0^+) \exp\left( -\int_0^t \frac{ds}{\tau_1(\lambda(s))} \right)\]

Key coupling effect: As structure rebuilds (\(\lambda \to 1\)), the viscosity \(\eta(\lambda)\) increases exponentially. This causes the relaxation time \(\tau_1 = \eta/G\) to become very large, effectively arresting the relaxation.

This is the signature of yielding behavior: a partially relaxed stress becomes “frozen in” as the material re-solidifies.

Two Timescales¶

The relaxation exhibits two timescales:

Fast elastic relaxation: \(\tau_1(\lambda_0) = \eta(\lambda_0)/G(\lambda_0)\) (initial, while structure is still broken)
Slow structural recovery: \(t_{eq}\) (controls rebuilding)

The observed relaxation is an interplay of these, with early fast decay transitioning to arrested relaxation as the material re-gels.

Creep (Step Stress)¶

Protocol: At \(t = 0\), apply constant stress \(\sigma(t) = \sigma_0 H(t)\).

This is the most diagnostic protocol for thixotropic yield-stress materials, revealing viscosity bifurcation and delayed yielding.

Inversion for Shear Rate¶

The constitutive equation must be inverted to find \(\dot{\gamma}(t)\):

\[\sigma_0 = \eta(\lambda(t), \dot{\gamma}(t)) \cdot \dot{\gamma}(t)\]

If \(\eta\) does not depend on \(\dot{\gamma}\) (pure structure-dependent viscosity):

\[\dot{\gamma}(t) = \frac{\sigma_0}{\eta(\lambda(t))}\]

If \(\eta\) depends on \(\dot{\gamma}\) (HB or other shear-thinning): solve the algebraic equation for \(\dot{\gamma}(t)\) at each time given \(\lambda(t)\).

Structure Evolution Under Creep¶

The kinetics become:

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}} - \frac{a \lambda |\dot{\gamma}(t)|^c}{t_{eq}}\]

with \(\dot{\gamma}(t)\) determined from the stress constraint.

Creep compliance:

\[J(t) = \frac{\gamma(t)}{\sigma_0} = \frac{1}{\sigma_0} \int_0^t \dot{\gamma}(s) \, ds\]

Viscosity Bifurcation¶

The creep response exhibits viscosity bifurcation [Coussot2002]:

Stress Regime	Behavior
\(\sigma_0 < \tau_y(\lambda = 1)\)	Arrested creep: Structure builds faster than it breaks. \(\lambda \to 1\), \(\eta \to \eta_0\), flow stops. \(J(t) \to \text{const}\) (solid-like).
\(\sigma_0 \approx \tau_y\)	Delayed yielding: Metastable state with slow creep. Eventually, catastrophic flow onset.
\(\sigma_0 > \tau_y\)	Immediate flow: Breakdown dominates. \(\lambda\) drops, \(\eta\) drops, \(\dot{\gamma}\) accelerates. \(J(t) \sim t\) at long times (liquid-like).

This bifurcation is the hallmark of yield-stress fluids: below a critical stress, \(\eta \to \infty\) (solid); above it, \(\eta \to\) finite (liquid).

Maxwell Variant Creep¶

With the Maxwell backbone, solve the coupled system:

Stress constraint:

\[\sigma_0 = \sigma_M(t) + \eta_s(\lambda(t)) \dot{\gamma}(t)\]

Rearranging:

\[\dot{\gamma}(t) = \frac{\sigma_0 - \sigma_M(t)}{\eta_s(\lambda(t))}\]

Maxwell stress evolution:

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda)} = G(\lambda) \dot{\gamma}(t)\]

Structure kinetics:

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}} - \frac{a \lambda |\dot{\gamma}(t)|^c}{t_{eq}}\]

The Maxwell variant adds an initial elastic jump:

\[\gamma(0^+) = \gamma_e(0) = \frac{\sigma_0}{G(\lambda_0)}\]

followed by combined elastic and viscous creep.

SAOS (Small Amplitude Oscillatory Shear)¶

Protocol: \(\gamma(t) = \gamma_0 \sin(\omega t)\) with \(\gamma_0 \ll 1\).

Limitations of Pure DMT¶

The pure DMT model (generalized Newtonian + structure kinetics) is not a linear viscoelastic model. It cannot produce genuine \(G'(\omega)\) and \(G''(\omega)\) in the standard sense because there is no elastic energy storage.

What pure DMT can predict:

Time-dependent apparent viscosity under oscillatory shear
Phase shifts only through structural kinetics (not through elastic storage)

For proper SAOS response, the Maxwell backbone (include_elasticity=True) is required.

Governing Equations (Full Nonlinear)¶

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

\[\sigma(t) = \eta(\lambda(t), \dot{\gamma}(t)) \cdot \dot{\gamma}(t)\]

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}} - \frac{a \lambda |\dot{\gamma}(t)|^c}{t_{eq}}\]

Linear Regime (Small Amplitude)¶

For small \(\gamma_0\) where breakdown is weak, linearize around \(\lambda \approx \lambda_0\):

\(|\dot{\gamma}(t)|^c \sim (\gamma_0 \omega)^c |\cos(\omega t)|^c\)

This is not purely sinusoidal, so even “small amplitude” produces harmonics in \(\lambda\) unless further approximations are made.

Maxwell Variant SAOS¶

With the Maxwell backbone at a fixed structure level \(\lambda_0\):

Complex modulus:

\[G^*(\omega) = \frac{i \omega G(\lambda_0) \tau_1(\lambda_0)}{1 + i \omega \tau_1(\lambda_0)}\]

Storage modulus:

\[G'(\omega) = G(\lambda_0) \frac{(\omega \tau_1)^2}{1 + (\omega \tau_1)^2}\]

Loss modulus:

\[G''(\omega) = G(\lambda_0) \frac{\omega \tau_1}{1 + (\omega \tau_1)^2} + \eta_s \omega\]

where \(\tau_1 = \tau_1(\lambda_0)\).

Limiting behaviors:

Low frequency: \(G' \sim \omega^2\), \(G'' \sim \omega\) (Maxwell liquid)
High frequency: \(G' \to G(\lambda_0)\), \(G'' \to \eta_s \omega\) (elastic solid + viscous)
Crossover: \(\omega_c = 1/\tau_1(\lambda_0)\)

For fully structured material (\(\lambda_0 = 1\)):

\(G'(\omega \to 0) \to G_0\) (solid-like plateau)
\(G''\) shows a minimum (liquid-like at very low \(\omega\), solid-like at intermediate \(\omega\))

LAOS (Large Amplitude Oscillatory Shear)¶

Protocol: \(\gamma(t) = \gamma_0 \sin(\omega t)\) with \(\gamma_0\) finite (nonlinear).

LAOS is the natural regime for DMT because large amplitude drives substantial breakdown/rebuild within cycles, making the structure kinetics dominant.

Full Governing System¶

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

Maxwell stress evolution:

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda(t))} = G(\lambda(t)) \dot{\gamma}(t)\]

Total stress:

\[\sigma(t) = \sigma_M(t) + \eta_s(\lambda(t)) \dot{\gamma}(t)\]

Structure kinetics:

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}} - \frac{a \lambda |\gamma_0 \omega \cos(\omega t)|^c}{t_{eq}}\]

Fourier Decomposition¶

The stress response is periodic but non-sinusoidal:

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

Only odd harmonics appear due to symmetry under \(\gamma \to -\gamma\).

First harmonic moduli:

\[G'_1 = \frac{\sigma'_1}{\gamma_0}, \quad G''_1 = \frac{\sigma''_1}{\gamma_0}\]

Third harmonic ratio (nonlinearity measure):

\[I_{3/1} = \frac{\sqrt{(\sigma'_3)^2 + (\sigma''_3)^2}}{\sqrt{(\sigma'_1)^2 + (\sigma''_1)^2}}\]

Chebyshev Decomposition¶

An alternative representation using Chebyshev polynomials of the first kind:

\[\sigma'(\gamma, \dot{\gamma}) = \gamma_0 \sum_n \left[ e_n T_n(x) + v_n T_n(y) \right]\]

where \(x = \gamma/\gamma_0\) and \(y = \dot{\gamma}/\dot{\gamma}_{max}\).

The coefficients \(e_n\) (elastic Chebyshev) and \(v_n\) (viscous Chebyshev) provide physical interpretation of the nonlinear response.

Intra-Cycle Structure Evolution¶

In LAOS, the structure parameter \(\lambda\) oscillates within each cycle:

Near \(|\dot{\gamma}| = \gamma_0 \omega\) (maximum shear rate): Strong breakdown, \(\lambda\) decreases
Near \(\dot{\gamma} = 0\) (strain extrema): Structure recovery, \(\lambda\) increases (if \(t_{eq}\) is comparable to period)

The amplitude of \(\lambda\) oscillation depends on the ratio \(2\pi / (\omega \cdot t_{eq})\):

\(\omega t_{eq} \gg 1\): Structure cannot follow, averages to intermediate value
\(\omega t_{eq} \ll 1\): Structure tracks instantaneous rate closely

LAOS Signatures from Thixotropy¶

DMT-type models predict characteristic LAOS features:

Feature	Physical Origin
Strain softening (\(G'_1\) decreases with \(\gamma_0\))	Structure breakdown at large \(\gamma_0\)
Higher harmonics (\(I_{3/1} > 0\))	Nonlinear structure kinetics
Asymmetric Lissajous curves	Different buildup/breakdown rates (\(t_{eq}\) vs \(t_{br}\))
Intra-cycle yielding	Stress peak before strain peak
Secondary loops in Lissajous	Elastic recoil combined with thixotropy

Pipkin Diagram¶

The Pipkin diagram maps behavior in (Wi, De) space:

Wi = \(\gamma_0 \omega \cdot t_{eq}\) (Weissenberg number, strain amplitude effect)
De = \(\omega \cdot t_{eq}\) (Deborah number, frequency effect)

Region	\(\gamma_0\)	\(\omega\)	Behavior
Linear viscoelastic	Low	Any	\(\lambda \approx 1\), standard \(G'\), \(G''\)
Quasi-steady thixotropic	High	Low	Structure equilibrates within cycle
Nonlinear viscoelastic	High	High	Viscoelastic nonlinearity, limited thixotropy
Thixotropic LAOS	Intermediate	Intermediate	Complex coupling of all effects

Fluidity-Maxwell Extension¶

This section describes the theoretical foundation for combining DMT structural kinetics with Maxwell-type viscoelastic stress dynamics. This extension is essential for capturing true stress relaxation, SAOS moduli, and elastic overshoots.

Jeffreys Form (Original Formulation)¶

The original DMT-style formulation uses a structure-dependent Jeffreys (three-element) viscoelastic backbone:

\[\tau_1(\lambda) \dot{\sigma} + \sigma = \eta(\lambda) \left( \dot{\gamma} + \tau_2(\lambda) \ddot{\gamma} \right)\]

where:

\(\tau_1(\lambda)\): Structure-dependent relaxation time
\(\tau_2(\lambda)\): Structure-dependent retardation time
\(\eta(\lambda)\): Structure-dependent viscosity scale

This arises from “a structure-dependent Maxwell element in parallel with a Newtonian element,” yielding the Jeffreys/Oldroyd-B form with times depending on structure.

Maxwell-in-Parallel Representation¶

For numerical implementation, the Maxwell-in-parallel form is preferred:

Decomposition:

\[\sigma(t) = \sigma_M(t) + \eta_s(\lambda) \dot{\gamma}(t)\]

where \(\sigma_M\) is the Maxwell branch stress and \(\eta_s\) is the parallel Newtonian viscosity.

Maxwell branch ODE:

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda)} = G(\lambda) \dot{\gamma}\]

Advantages:

Avoids computing \(\ddot{\gamma}\) (difficult in rate-controlled experiments)
Natural for time-stepping numerical schemes
Clear separation of elastic and viscous contributions
Easy initialization: \(\sigma_M(0) = 0\) for stress-free initial state

Material Functions¶

The structure-dependent material functions are:

Relaxation time:

\[\tau_1(\lambda) = \frac{\eta_M(\lambda)}{G(\lambda)}\]

Common choices:

Linear: \(\tau_1(\lambda) = \tau_0 + \tau_{slope} \cdot \lambda\)
Reciprocal: \(\tau_1(\lambda) = \tau_0 / \lambda\) (diverges as \(\lambda \to 0\))
Power-law: \(\tau_1(\lambda) = \tau_0 \lambda^{m_\tau}\)

Elastic modulus:

\[G(\lambda) = G_0 \lambda^{m_G}\]

Parallel viscosity (often constant):

\[\eta_s(\lambda) = \eta_{s0}\]

Protocol Equations with Maxwell Backbone¶

Flow curve (steady state, \(\dot{\sigma}_M = 0\)):

\[\sigma_{ss} = G(\lambda_{eq}) \tau_1(\lambda_{eq}) \dot{\gamma} + \eta_s(\lambda_{eq}) \dot{\gamma}\]

where \(\lambda_{eq} = 1/(1 + a|\dot{\gamma}|^c)\).

Start-up (constant \(\dot{\gamma}_0\)):

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda(t))} = G(\lambda(t)) \dot{\gamma}_0\]

\[\lambda(t) = \lambda_{ss} + (\lambda_0 - \lambda_{ss}) e^{-rt}\]

This can produce overshoot-like behavior depending on how fast \(\tau_1\) collapses under shear (rejuvenation).

Cessation (stop shear at \(t = 0\)):

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda(t))} = 0, \quad t > 0\]

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}}\]

As \(\lambda\) rebuilds, \(\tau_1\) increases, causing stress relaxation to slow down and eventually arrest.

Creep (constant \(\sigma_0\)):

Solve for \(\dot{\gamma}\) from the Maxwell relation:

\[\dot{\gamma}(t) = \frac{1}{G(\lambda)} \left[ \dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda)} \right]\]

With \(\dot{\sigma} = 0\) in creep:

\[\dot{\gamma}(t) = \frac{\sigma_0}{G(\lambda(t)) \tau_1(\lambda(t))}\]

Then:

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}} - a \lambda \left| \frac{\sigma_0}{G(\lambda) \tau_1(\lambda)} \right|^c \frac{1}{t_{eq}}\]

SAOS (small amplitude, fixed \(\lambda \approx \lambda_*\)):

\[G^*(\omega) = \frac{i \omega G(\lambda_*) \tau_1(\lambda_*)}{1 + i \omega \tau_1(\lambda_*)} + i \omega \eta_s\]

Giving:

\[G'(\omega) = G(\lambda_*) \frac{(\omega \tau_1)^2}{1 + (\omega \tau_1)^2}\]

\[G''(\omega) = G(\lambda_*) \frac{\omega \tau_1}{1 + (\omega \tau_1)^2} + \eta_s \omega\]

Aging enters through slow time-dependence of \(\lambda_*\).

LAOS (full nonlinear):

\[\dot{\sigma}_M + \frac{\sigma_M}{\tau_1(\lambda(t))} = G(\lambda(t)) \gamma_0 \omega \cos(\omega t)\]

\[\frac{d\lambda}{dt} = \frac{1 - \lambda}{t_{eq}} - a \lambda |\gamma_0 \omega \cos(\omega t)|^c \frac{1}{t_{eq}}\]

Nonlocal Extension for Shear Banding¶

This section describes the spatial extension of the DMT model for capturing shear banding through fluidity/structure diffusion.

Physical Motivation¶

Shear banding occurs when the material separates into coexisting regions of different shear rates under uniform stress. Local DMT models can predict S-shaped flow curves (multiple solutions), but the sharp band interfaces are ill-posed without spatial regularization.

The nonlocal fluidity model introduces a diffusive coupling that:

Sets a minimum shear band width (cooperativity length)
Regularizes the mathematical problem
Captures the spatial structure evolution

Governing Equations¶

Spatial extension: Promote \(\lambda\) to a field \(\lambda(y, t)\) where \(y\) is the gap position:

\[\frac{\partial \lambda}{\partial t} = \underbrace{\frac{1 - \lambda}{t_{eq}} - \frac{a \lambda |\dot{\gamma}|^c}{t_{eq}}}_{\text{local kinetics}} + \underbrace{D_\lambda \frac{\partial^2 \lambda}{\partial y^2}}_{\text{structure diffusion}}\]

Stress balance (planar Couette, low Reynolds number):

\[\sigma(y, t) = \Sigma(t) \quad \text{(uniform across gap)}\]

The stress is uniform, but \(\dot{\gamma}(y, t)\) varies spatially according to the local constitutive relation:

\[\Sigma(t) = \eta(\lambda(y, t), \dot{\gamma}(y, t)) \cdot \dot{\gamma}(y, t)\]

Constraint: The spatially-averaged shear rate equals the imposed value:

\[\frac{1}{H} \int_0^H \dot{\gamma}(y, t) \, dy = \dot{\gamma}_{avg}\]

where \(H\) is the gap width.

Cooperativity Length¶

The structure diffusion introduces a characteristic length scale:

\[\xi = \sqrt{D_\lambda \cdot t_{eq}}\]

This cooperativity length sets the minimum shear band width. Typical values:

\(D_\lambda \sim 10^{-9}\) to \(10^{-12}\) m^2/s
\(t_{eq} \sim 1\) to \(1000\) s
\(\xi \sim 1~\mu\text{m}\) to \(1\) mm

Physical interpretation: Structure rearrangements are not purely local but involve cooperative motion of neighboring material elements over distance \(\xi\).

Boundary Conditions¶

No-flux (common for rheometer walls):

\[\frac{\partial \lambda}{\partial y}\bigg|_{y=0} = \frac{\partial \lambda}{\partial y}\bigg|_{y=H} = 0\]

Fixed structure (idealized smooth/rough walls):

\[\lambda(0, t) = \lambda_{wall}, \quad \lambda(H, t) = \lambda_{wall}\]

The boundary condition choice affects band nucleation location.

Shear Banding Criteria¶

Shear banding occurs when:

The steady-state flow curve is non-monotonic (S-shaped)
The applied stress lies in the unstable region (negative slope)
The diffusion length \(\xi\) is much smaller than the gap \(H\)

Band contrast: The ratio of shear rates in high-shear vs low-shear bands:

\[C = \frac{\dot{\gamma}_{high}}{\dot{\gamma}_{low}}\]

For strong banding, \(C \gg 1\).

Lever rule: At steady state, the volume fractions of bands satisfy:

\[f_{high} \dot{\gamma}_{high} + (1 - f_{high}) \dot{\gamma}_{low} = \dot{\gamma}_{avg}\]

Transient Banding¶

Shear bands can be:

Transient: Appear during startup but disappear at steady state
Steady-state: Persist indefinitely at fixed conditions
Oscillatory: Bands move or oscillate under steady forcing

The DMT nonlocal model can capture all these behaviors depending on parameters.

Connection to Fluidity Models¶

The structure parameter \(\lambda\) and fluidity \(\phi\) are related by:

\[\phi = \frac{1}{\eta(\lambda)}\]

For the exponential closure:

\[\phi = \frac{1}{\eta_\infty} \left( \frac{\eta_\infty}{\eta_0} \right)^\lambda\]

The fluidity-based kinetics:

\[\frac{\partial \phi}{\partial t} = -\frac{\phi - \phi_{min}}{t_{age}} + c |\dot{\gamma}|^m (\phi_{max} - \phi) + D_\phi \nabla^2 \phi\]

is mathematically equivalent to the \(\lambda\)-form after appropriate transformation.

Industrial Applications¶

The DMT model family is particularly well-suited for industrially relevant materials that exhibit combined thixotropy, yield stress, and viscoelasticity.

Waxy Crude Oils¶

Waxy crude oils form gel structures at temperatures below the wax appearance temperature. DMT models capture:

Gel strength (yield stress) from wax crystal network
Thixotropic breakdown during pipeline startup
Structure recovery during shut-in periods
Restart pressure prediction for pipeline design

Key parameters: High \(\tau_{y0}\), large \(t_{eq}\) (slow recovery), temperature-dependent.

Drilling Fluids and Muds¶

Water-based and oil-based drilling muds exhibit complex thixotropic behavior:

Gel strength for cuttings suspension during circulation stops
Low viscosity under high shear for efficient pumping
Progressive gel development over time

The HB closure is preferred for explicit yield stress modeling.

Concrete and Cement Slurries¶

Fresh concrete and cement pastes show pronounced thixotropy:

Formwork pressure prediction requires thixotropic modeling
Pumpability assessment
Self-compacting concrete (SCC) flow design

Cementitious materials often show strong structure recovery (\(t_{eq} \sim\) minutes).

Food Products¶

Many food materials are thixotropic gels:

Mayonnaise: Thixotropic emulsion, important for dispensing
Yogurt: Structure breakdown affects mouthfeel
Ketchup: Classic yield-stress thixotropic material

Cosmetics and Personal Care¶

Lotions and creams: Must flow during application, then stop
Toothpaste: Yield stress for tube stability, thixotropy for spreading
Hair gel: Structure recovery for styling hold

Connection to Other Model Families¶

Soft Glassy Rheology (SGR)¶

The SGR model [Sollich1997] and DMT model address similar physics from different perspectives:

Aspect	DMT	SGR
State variable	Structure \(\lambda\) (scalar)	Trap depth distribution \(P(E)\)
Kinetics	Phenomenological buildup/breakdown	Activated hopping with effective temperature
Yield stress	Explicit via closure	Emergent from trap distribution
Aging	\(\lambda \to 1\) as \(t \to \infty\)	Mean trap depth increases
SAOS	Requires Maxwell backbone	Natural from trap dynamics

Both models predict similar phenomena (aging, rejuvenation, yield stress) through different mechanisms.

Kinetic Elasto-Viscoplastic (KEVP) Models¶

The DMT model is a member of the broader KEVP class:

Saramito model: Combines Oldroyd-B viscoelasticity with Herschel-Bulkley plasticity
Bautista-Manero-Puig (BMP): Structure-dependent Maxwell model
Mujumdar model: Elastic strain with yield function

The DMT model is distinguished by its explicit structure-dependent yield stress and comprehensive treatment of structure kinetics.

API Reference¶

DMTLocal¶

class rheojax.models.dmt.DMTLocal(closure='exponential', include_elasticity=True)[source]

Bases: DMTBase

Local (0D) DMT model for homogeneous thixotropic flow.

This model assumes spatially homogeneous flow (no shear banding). For shear banding analysis, use DMTNonlocal.

The model captures: - Yielding: Stress plateau at low shear rates (HB closure) - Thixotropy: Time-dependent viscosity via structure kinetics - Viscoelasticity: Optional Maxwell backbone for overshoot/SAOS

Two viscosity closures: - “exponential”: η(λ) = η_∞·(η_0/η_∞)^λ (smooth, original DMT) - “herschel_bulkley”: Explicit yield stress τ_y(λ) + K(λ)|γ̇|^n

Parameters:

closure (Literal['exponential', 'herschel_bulkley']) – Viscosity closure type.
include_elasticity (bool) – Include Maxwell viscoelastic backbone for stress overshoot and SAOS.

Examples

>>> from rheojax.models.dmt import DMTLocal
>>>
>>> # Create model with Herschel-Bulkley closure
>>> model = DMTLocal(closure="herschel_bulkley", include_elasticity=True)
>>>
>>> # Fit to flow curve data
>>> model.fit(gamma_dot, stress, test_mode="flow_curve")
>>>
>>> # Simulate startup shear
>>> t, stress, lam = model.simulate_startup(gamma_dot=10.0, t_end=100.0)

DMTNonlocal¶

class rheojax.models.dmt.DMTNonlocal(closure='exponential', include_elasticity=True, n_points=51, gap_width=0.001)[source]

Bases: DMTBase

Nonlocal (1D) DMT model for shear banding analysis.

Extends the local DMT model with spatial structure diffusion:

∂λ/∂t = (1-λ)/t_eq - a·λ·|γ̇|^c/t_eq + D_λ·∂²λ/∂y²

The diffusion term introduces a cooperativity length scale: ξ ~ √(D_λ · t_eq)

which regularizes the problem and sets the minimum width of shear bands.

This model solves for: - λ(y, t): Structure field across the gap - v(y, t): Velocity profile (from momentum balance) - γ̇(y, t): Local shear rate

Parameters:

closure (Literal['exponential', 'herschel_bulkley']) – Viscosity closure type.
include_elasticity (bool) – Include Maxwell viscoelastic backbone.
n_points (int) – Number of spatial grid points across the gap.
gap_width (float) – Gap width H [m] (e.g., for Couette cell).

n_points

Spatial grid resolution

Type:: int

gap_width

Physical gap width [m]

Type:: float

y

Spatial coordinate array [m]

Type:: array

Examples

>>> from rheojax.models.dmt import DMTNonlocal
>>>
>>> # Create nonlocal model for banding analysis
>>> model = DMTNonlocal(
...     closure="herschel_bulkley",
...     n_points=101,
...     gap_width=1e-3
... )
>>>
>>> # Simulate steady shear with banding
>>> result = model.simulate_steady_shear(
...     gamma_dot_avg=10.0, t_end=1000.0
... )
>>>
>>> # Check for banding
>>> banding_info = model.detect_banding(result)

Parameter Estimation Guidance¶

Systematic parameter estimation improves model fits and convergence. Below are recommended strategies for extracting DMT parameters from different experiments.

From Steady-State Flow Curve¶

The flow curve \(\sigma(\dot{\gamma})\) provides several key parameters:

Parameter	Estimation Method
\(\tau_{y0}\)	Extrapolate low-shear plateau; fit Herschel-Bulkley to low \(\dot{\gamma}\) data
\(K_0, n\)	Power-law fit to intermediate shear rate region above yield
\(\eta_\infty\)	High-shear slope: \(\sigma/\dot{\gamma} \to \eta_\infty\) as \(\dot{\gamma} \to \infty\)
\(a, c\)	Fit equilibrium structure: \(\lambda_{eq} = 1/(1 + a\|\dot{\gamma}\|^c)\)

The shape of the shear-thinning transition (log-log curvature) constrains \(a\) and \(c\).

From Transient Tests¶

Startup and cessation experiments provide kinetic and elastic parameters:

Parameter	Estimation Method
\(G_0\)	Initial slope of startup stress: \(\sigma(t \ll t_{eq}) \approx G_0 \cdot \gamma\)
\(t_{eq}\)	Recovery time constant from structure rebuild after cessation
\(a, c\)	Overshoot decay rate; rate-dependence of peak strain

The stress overshoot ratio \(\sigma_{peak}/\sigma_{ss}\) increases with shear rate and initial structure level.

From SAOS¶

For the Maxwell variant, the linear viscoelastic moduli provide:

Parameter	Estimation Method
\(G_0\)	High-frequency plateau: \(G'(\omega \to \infty) \to G_0\)
\(\theta_0\)	Crossover frequency: \(\omega_c\) where \(G' = G''\); \(\theta_0 = 1/\omega_c\)
\(\eta_0\)	Low-frequency loss: \(G''(\omega)/\omega \to \eta_0\) as \(\omega \to 0\)

Typical Parameter Ranges¶

Parameter	Typical Range	Unit	Notes
\(\tau_{y0}\)	1 - 1000	Pa	Higher for strongly structured gels
\(K_0\)	0.1 - 100	Pa·sⁿ	Consistency at full structure
\(n\)	0.2 - 0.8	—	Shear-thinning index
\(\eta_\infty\)	0.001 - 1	Pa·s	Often close to solvent viscosity
\(G_0\)	10 - 10000	Pa	Plateau modulus
\(m_1, m_2, m_G\)	0.5 - 2.0	—	Structure exponents (often ~1)
\(t_{eq}\)	1 - 10000	s	Longer for strong gels
\(a\)	0.01 - 100	—	Breakdown intensity
\(c\)	0.5 - 2.0	—	Often ~1 for linear kinetics

Fitting Strategy¶

Recommended order for parameter estimation:

Flow curve first: Fix \(\tau_{y0}\), \(K_0\), \(n\), \(\eta_\infty\) from steady-state data
Add transients: Fit \(t_{eq}\), \(a\), \(c\) from startup/cessation
Add elasticity: Fit \(G_0\), \(m_G\) from SAOS or overshoot magnitude
Refine jointly: Use Bayesian inference with all protocols simultaneously

Usage Examples¶

Flow Curve Fitting¶

import numpy as np
from rheojax.models import DMTLocal

# Experimental flow curve data
gamma_dot = np.logspace(-2, 2, 20)  # 1/s
stress = np.array([...])  # Pa

# Fit with exponential closure
model = DMTLocal(closure="exponential", include_elasticity=True)
model.fit(gamma_dot, stress, test_mode='flow_curve')

# Predict flow curve
gamma_dot_pred = np.logspace(-3, 3, 100)
stress_pred = model.predict(gamma_dot_pred, test_mode='flow_curve')

Startup Shear with Stress Overshoot¶

from rheojax.models import DMTLocal

model = DMTLocal(closure="exponential", include_elasticity=True)

# Startup at γ̇ = 10 s⁻¹ from fully-structured state
t, stress, lam = model.simulate_startup(
    gamma_dot=10.0,
    t_end=100.0,
    dt=0.01,
    lam_init=1.0  # Aged state
)

# Find stress overshoot
peak_idx = np.argmax(stress)
overshoot_ratio = stress[peak_idx] / stress[-1]

Creep with Delayed Yielding¶

The creep response differs significantly between viscous and Maxwell variants.

Viscous Variant (include_elasticity=False):

Pure viscous flow: \(\gamma(t) = \int_0^t \sigma_0 / \eta(\lambda(s)) \, ds\)

from rheojax.models import DMTLocal

model = DMTLocal(closure="herschel_bulkley", include_elasticity=False)

# Apply constant stress
t, gamma, gamma_dot, lam = model.simulate_creep(
    sigma_0=50.0,  # Applied stress (Pa)
    t_end=1000.0,
    dt=0.1,
    lam_init=1.0  # Start from aged state
)

# Observe delayed yielding: initial slow creep, then acceleration
# as structure breaks down

Maxwell Variant (include_elasticity=True):

Total strain includes both elastic and viscous contributions:

\[\gamma(t) = \underbrace{\frac{\sigma_0}{G(\lambda(t))}}_{\gamma_e(t)} + \underbrace{\int_0^t \frac{\sigma_0}{\eta(\lambda(s))} \, ds}_{\gamma_v(t)}\]

Key features:

Initial elastic jump: \(\gamma(0^+) = \sigma_0 / G(\lambda_0)\) — instantaneous response
Elastic strain evolution: As structure breaks down (\(\lambda \downarrow\)), \(G \downarrow\), so \(\gamma_e \uparrow\)
Viscous flow: Accumulates continuously via \(\dot{\gamma}_v = \sigma_0 / \eta(\lambda)\)

from rheojax.models import DMTLocal
import numpy as np

model = DMTLocal(closure="exponential", include_elasticity=True)

# Parameters for initial elastic strain estimate
G0 = model.parameters.get_value("G0")
sigma_0 = 100.0  # Pa

# Expected initial elastic strain: γ_e(0) = σ₀/G₀
gamma_e_expected = sigma_0 / G0
print(f"Expected initial elastic strain: {gamma_e_expected:.4f}")

# Simulate creep
t, gamma, gamma_dot, lam = model.simulate_creep(
    sigma_0=sigma_0,
    t_end=500.0,
    dt=0.1,
    lam_init=1.0
)

# Verify initial strain includes elastic contribution
print(f"Actual initial strain: {gamma[0]:.4f}")

# As structure breaks (λ decreases), elastic strain increases
# because G(λ) = G₀·λ^m_G decreases
print(f"Structure: {lam[0]:.3f} → {lam[-1]:.3f}")
print(f"Final strain: {gamma[-1]:.4f}")

SAOS Predictions (Maxwell Variant)¶

import numpy as np
from rheojax.models import DMTLocal

model = DMTLocal(closure="exponential", include_elasticity=True)

omega = np.logspace(-3, 3, 50)  # rad/s
G_prime, G_double_prime = model.predict_saos(omega, lam_0=1.0)

# Crossover frequency ω_c where G' = G''
# Related to relaxation time θ = η₀/G₀

LAOS Analysis¶

Pipkin Diagram (Wi-De space):

The nonlinear oscillatory response can be classified using the Pipkin diagram, which maps behavior in the (strain amplitude, frequency) space:

Region	Conditions	Behavior
Linear viscoelastic	Low \(\gamma_0\), any \(\omega\)	\(\lambda \approx 1\), standard \(G'\), \(G''\)
Quasi-steady thixotropic	High \(\gamma_0\), low \(\omega\)	Structure equilibrates within cycle
Nonlinear viscoelastic	High \(\gamma_0\), high \(\omega\)	Viscoelastic nonlinearity, limited thixotropy
Thixotropic LAOS	Intermediate	Complex coupling of all effects

Intra-cycle structure evolution:

In LAOS, the structure parameter \(\lambda\) oscillates within each cycle:

Near \(|\dot{\gamma}| = \gamma_0\omega\) (maximum shear rate): Strong breakdown, \(\lambda\) decreases
Near \(\dot{\gamma} = 0\) (strain extrema): Structure recovery, \(\lambda\) increases

This intra-cycle variation produces:

Strain softening: \(G'_1\) decreases with \(\gamma_0\)
Higher harmonics: \(I_{3/1} > 0\) from nonlinear structure kinetics
Asymmetric Lissajous curves: Different buildup/breakdown rates
Secondary loops: Elastic recoil combined with thixotropy

from rheojax.models import DMTLocal

model = DMTLocal(closure="exponential", include_elasticity=True)

# Simulate LAOS
result = model.simulate_laos(
    gamma_0=0.5,  # Strain amplitude
    omega=1.0,    # Angular frequency (rad/s)
    n_cycles=10,
    points_per_cycle=128
)

# Extract harmonics
harmonics = model.extract_harmonics(result, n_harmonics=5)
# harmonics["G1_prime"], harmonics["G3_prime"], etc.

Shear Banding with Nonlocal Model¶

from rheojax.models import DMTNonlocal

model = DMTNonlocal(
    closure="exponential",
    include_elasticity=True,
    n_points=101,
    gap_width=1e-3  # 1 mm gap
)

# Simulate steady shear
result = model.simulate_steady_shear(
    gamma_dot_avg=10.0,  # Average shear rate
    t_end=500.0,
    dt=1.0
)

# Detect banding
banding_info = model.detect_banding(result, threshold=0.1)
print(f"Shear banding: {banding_info['is_banding']}")
print(f"Band contrast: {banding_info['band_contrast']:.2f}")

Numerical Implementation¶

ODE Integration¶

Time-stepping simulations use jax.lax.scan for efficient compilation:

def step(state, _):
    lam, sigma = state
    # Update structure
    dlam = structure_evolution(lam, gamma_dot, t_eq, a, c)
    lam_new = clip(lam + dt * dlam, 0, 1)
    # Update stress (Maxwell)
    dsigma = G(lam) * gamma_dot - sigma / theta(lam)
    sigma_new = sigma + dt * dsigma
    return (lam_new, sigma_new), (sigma_new, lam_new)

_, (stress, lam) = jax.lax.scan(step, init_state, None, length=n_steps)

JIT Compilation¶

All core kernels are JIT-compiled for performance:

@jax.jit
def equilibrium_structure(gamma_dot, a, c):
    return 1.0 / (1.0 + a * jnp.abs(gamma_dot) ** c)

@jax.jit
def viscosity_exponential(lam, eta_0, eta_inf):
    return eta_inf * jnp.power(eta_0 / eta_inf, lam)

First compilation may take 1-2 seconds; subsequent calls are fast.

Papanastasiou Regularization¶

For numerical stability, the Herschel-Bulkley yield stress is regularized:

\[\sigma = \tau_y \left(1 - e^{-m|\dot{\gamma}|}\right) + K|\dot{\gamma}|^n + \eta_\infty\dot{\gamma}\]

where \(m\) is a large regularization parameter (default: 1000).

Literature References¶

[deSouzaMendes2009]

de Souza Mendes, P. R. (2009). “Modeling the thixotropic behavior of structured fluids.” Journal of Non-Newtonian Fluid Mechanics, 164(1-3), 66-75. https://doi.org/10.1016/j.jnnfm.2009.08.005

[Mendes2012]

de Souza Mendes, P. R., & Thompson, R. L. (2012). “A critical overview of elasto-viscoplastic thixotropic modeling.” Journal of Non-Newtonian Fluid Mechanics, 187-188, 8-15. https://doi.org/10.1016/j.jnnfm.2012.08.006

[Thompson2014]

de Souza Mendes, P. R., & Thompson, R. L. (2013). “A unified approach to model elasto-viscoplastic thixotropic yield-stress materials and apparent yield-stress fluids.” Rheologica Acta, 52(7), 673-694. https://doi.org/10.1007/s00397-013-0699-1

[Coussot2002] (1,2)

Coussot, P., Nguyen, Q. D., Huynh, H. T., & Bonn, D. (2002). “Avalanche behavior in yield stress fluids.” Physical Review Letters, 88(17), 175501. https://doi.org/10.1103/PhysRevLett.88.175501

[Mujumdar2002]

Mujumdar, A., Beris, A. N., & Metzner, A. B. (2002). “Transient phenomena in thixotropic systems.” Journal of Non-Newtonian Fluid Mechanics, 102(2), 157-178. https://doi.org/10.1016/S0377-0257(01)00176-8

[Dullaert2006]

Dullaert, K., & Mewis, J. (2006). “A structural kinetics model for thixotropy.” Journal of Non-Newtonian Fluid Mechanics, 139(1-2), 21-30. https://doi.org/10.1016/j.jnnfm.2006.06.002

[Larson2019]

Larson, R. G., & Wei, Y. (2019). “A review of thixotropy and its rheological modeling.” Journal of Rheology, 63(3), 477-501. https://doi.org/10.1122/1.5055031

[MewisWagner2009]

Mewis, J., & Wagner, N. J. (2009). “Thixotropy.” Advances in Colloid and Interface Science, 147-148, 214-227. https://doi.org/10.1016/j.cis.2008.09.005

[Sollich1997]

Sollich, P., Lequeux, F., Hébraud, P., & Cates, M. E. (1997). “Rheology of soft glassy materials.” Physical Review Letters, 78(10), 2020-2023. https://doi.org/10.1103/PhysRevLett.78.2020

[SoftMatter2011]

de Souza Mendes, P. R. (2011). “Thixotropic elasto-viscoplastic model for structured fluids.” Soft Matter, 7(6), 2471-2483. https://doi.org/10.1039/C0SM01021A

DMT Thixotropic Models¶

Quick Reference¶

Notation Guide¶

Overview¶

Key Features¶

Historical Context¶

Jeffreys/Oldroyd-B Backbone¶

Evolution from Simple Thixotropy¶

Physical Foundations¶

Structure Parameter¶

Structure Kinetics¶

Fluidity Interpretation¶

Viscosity Closures¶

Exponential Closure¶

Herschel-Bulkley Closure¶

Maxwell Viscoelasticity¶

Steady-State Flow Curve¶

Exponential Closure¶

Herschel-Bulkley Closure¶

Viscosity Bifurcation¶

Critical Stress¶

Avalanche Effect¶

Connection to Herschel-Bulkley¶

Related Models Comparison¶

Coussot Model¶

Houska Model¶

Mujumdar Model¶

When to Use Which Model¶

Parameters¶

Core Viscosity Parameters¶

Herschel-Bulkley Parameters (closure="herschel_bulkley" only)¶

Elastic Parameters (include_elasticity=True only)¶

Structure Kinetics Parameters¶

Nonlocal Parameters (DMTNonlocal only)¶

Dimensionless Groups¶

Deborah Number (De)¶

Weissenberg Number (Wi)¶

Structure Number (Sn)¶

Regime Classification¶

Protocol-Specific Equations¶

General Governing System¶

Rotation / Steady-State Flow Curve¶

Equilibrium Condition¶

Steady-State Stress¶

Controlled-Stress Mode¶

Start-up of Steady Shear¶

Governing ODEs¶

Stress Overshoot Mechanism¶

Purely Generalized Newtonian (No Elasticity)¶

Stress Relaxation (Cessation of Flow)¶

Structure Recovery¶

Stress Relaxation (Generalized Newtonian)¶

Stress Relaxation (Maxwell Backbone)¶

Two Timescales¶

Creep (Step Stress)¶

Inversion for Shear Rate¶

Structure Evolution Under Creep¶

Viscosity Bifurcation¶

Maxwell Variant Creep¶

SAOS (Small Amplitude Oscillatory Shear)¶

Limitations of Pure DMT¶

Governing Equations (Full Nonlinear)¶

Linear Regime (Small Amplitude)¶

Maxwell Variant SAOS¶

LAOS (Large Amplitude Oscillatory Shear)¶

Full Governing System¶

Fourier Decomposition¶

Chebyshev Decomposition¶

Intra-Cycle Structure Evolution¶

LAOS Signatures from Thixotropy¶

Pipkin Diagram¶

Fluidity-Maxwell Extension¶

Jeffreys Form (Original Formulation)¶

Maxwell-in-Parallel Representation¶

Material Functions¶

Protocol Equations with Maxwell Backbone¶

Nonlocal Extension for Shear Banding¶

Physical Motivation¶

Governing Equations¶

Cooperativity Length¶

Herschel-Bulkley Parameters (`closure="herschel_bulkley"` only)¶

Elastic Parameters (`include_elasticity=True` only)¶

Nonlocal Parameters (`DMTNonlocal` only)¶