ODE-Based Constitutive Models

Overview

This tutorial covers constitutive models that solve systems of ordinary differential equations (ODEs) to predict transient rheological behavior. Unlike analytical models (Maxwell, fractional), these models use numerical ODE integration via diffrax to compute time-dependent responses in startup shear, creep, stress relaxation, and large-amplitude oscillatory shear (LAOS) protocols.

ODE-based models capture rich physics including:

  • Stress overshoot in startup shear (thixotropic/structural evolution)

  • Normal stress differences (\(N_1\), \(N_2\)) in shear flow

  • Yield stress behavior with Von Mises criterion

  • Thixotropic loops (buildup/breakdown cycles)

  • Nonlinear viscoelasticity beyond linear response theory

RheoJAX implements three major ODE-based model families:

  1. Giesekus (nonlinear polymer viscoelasticity)

  2. IKH/FIKH (isotropic-kinematic hardening with optional fractional memory)

  3. Fluidity-Saramito EVP (elastoviscoplastic fluids with thixotropic structure evolution)

Key Insight

ODE-based models excel at transient protocols where analytical solutions don’t exist. The trade-off is computational cost: each prediction requires numerical integration, making fits 10-100× slower than analytical models. Use NLSQ warm-start before Bayesian inference.

When to Use ODE-Based Models

Choose ODE-based models when:

  • Transient phenomena dominate: Stress overshoot, creep acceleration, thixotropic loops

  • Normal stress predictions needed: Giesekus predicts \(N_1/N_2\) from shear flow

  • Yield stress with elasticity: Saramito EVP combines Bingham yield with viscoelastic backbone

  • Fractional memory in thixotropy: FIKH/FMLIKH for long-time structural memory

  • Shear banding (nonlocal variants): Spatial heterogeneity in flow

Don’t use ODE models when:

  • Linear viscoelasticity suffices (use Maxwell, Fractional Maxwell instead)

  • Only steady-state flow curves needed (use Carreau-Yasuda, power-law)

  • Fast fitting required (>1000 datasets in batch)

Warning

ODE models compile slowly on first call (~30-90s for diffrax JIT). Pre-compile with model.precompile() if running repeated predictions. Each NUTS leapfrog step requires forward+backward ODE solves, limiting Bayesian sample counts.

Theoretical Foundations

Constitutive Equations as ODEs

Classical analytical models have closed-form stress-strain relations:

\[\sigma(t) = \int_0^t G(t-t') \dot{\gamma}(t') dt' \quad \text{(Maxwell)}\]

ODE models instead evolve internal state variables \(\mathbf{s}(t)\) with differential equations:

\[\frac{d\mathbf{s}}{dt} = \mathcal{F}(\mathbf{s}, \dot{\gamma}, T, \ldots)\]

where \(\mathcal{F}\) is a nonlinear operator. The stress \(\sigma(t)\) is computed from the current state \(\sigma = \mathcal{G}(\mathbf{s})\).

Three Key State Variables

  1. Stress/conformation tensor \(\boldsymbol{\tau}\) or \(\boldsymbol{\mu}\) (Giesekus, Saramito)

    • Tracks anisotropic polymer/network deformation

    • Upper-convected derivative: \(\nabla \boldsymbol{\tau} = \partial_t \boldsymbol{\tau} + \mathbf{v} \cdot \nabla \boldsymbol{\tau} - (\nabla \mathbf{v})^T \cdot \boldsymbol{\tau} - \boldsymbol{\tau} \cdot \nabla \mathbf{v}\)

  2. Structure parameter \(\lambda \in [0,1]\) (IKH, Saramito, DMT)

    • Dimensionless measure of microstructural state (0 = broken, 1 = intact)

    • Evolves via buildup/breakdown kinetics

  3. Backstress tensor \(\boldsymbol{\alpha}\) (kinematic hardening, IKH only)

    • Internal stress arising from structural rearrangements

    • Causes thixotropic hysteresis in loading/unloading

Key Insight

The dimensionality of the ODE system determines computational cost. Giesekus (5 components: 3 tensor + 2 symmetry) is faster than Saramito (4 components + fluidity). Multi-mode models multiply the ODE size by the number of modes.

Model Family 1: Giesekus

Physical Basis

The Giesekus model describes nonlinear viscoelastic polymers (solutions, melts) with:

  • Anisotropic drag coefficient (mobility parameter \(\alpha\))

  • Shear thinning via stress-dependent relaxation

  • Normal stress differences (\(N_1 > 0\), \(N_2 < 0\))

Constitutive equation:

\[\boldsymbol{\tau} + \lambda \left( \nabla \boldsymbol{\tau} + \frac{\alpha}{G} \boldsymbol{\tau} \cdot \boldsymbol{\tau} \right) = 2\eta_0 \mathbf{D}\]

where \(\mathbf{D}\) is the rate-of-deformation tensor, \(\lambda\) is relaxation time, \(G = \eta_0/\lambda\) is modulus, and \(\alpha \in [0, 0.5]\) controls nonlinearity.

Special case \(\alpha = 0\): upper-convected Maxwell (UCM) model.

Parameters and Protocols

Single-mode parameters (3):

  • G: Elastic modulus (Pa)

  • lambda_1: Relaxation time (s)

  • alpha: Mobility parameter (0-0.5, dimensionless)

Multi-mode parameters (per mode: 3):

  • G_i, lambda_i, alpha_i for mode \(i=1,\ldots,N\)

Supported protocols (6):

  • FLOW_CURVE: Steady-state shear stress vs. shear rate

  • OSCILLATION: SAOS (\(G'\), \(G''\) vs. frequency)

  • STARTUP: Transient stress growth at constant shear rate

  • RELAXATION: Stress decay after step strain

  • CREEP: Strain growth under constant stress

  • LAOS: Large-amplitude oscillatory shear (nonlinear harmonics)

Practical Implementation

Basic SAOS Fitting

from rheojax.models import GiesekusSingleMode
import numpy as np

# Load experimental data
omega = np.logspace(-2, 2, 50)  # rad/s
G_prime = ...  # Pa
G_double_prime = ...  # Pa
G_star = G_prime + 1j * G_double_prime

# Fit Giesekus model
model = GiesekusSingleMode()
model.fit(omega, G_star, test_mode='oscillation')

# Extract parameters
G = model.parameters.get_value('G')
lambda_1 = model.parameters.get_value('lambda_1')
alpha = model.parameters.get_value('alpha')

print(f"Modulus: {G:.1f} Pa")
print(f"Relaxation time: {lambda_1:.3f} s")
print(f"Mobility: {alpha:.3f}")

Startup Shear with Overshoot

# Predict stress growth at constant shear rate
t = np.linspace(0.01, 50, 300)  # s
gamma_dot = 1.0  # s⁻¹

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

# Find overshoot
idx_max = np.argmax(stress)
t_max = t[idx_max]
sigma_max = stress[idx_max]
sigma_ss = stress[-1]

print(f"Overshoot at t={t_max:.2f}s: {sigma_max:.1f} Pa")
print(f"Steady-state: {sigma_ss:.1f} Pa")
print(f"Overshoot ratio: {sigma_max/sigma_ss:.2f}")

Normal Stress Differences

# Predict N₁ and N₂ in steady shear
gamma_dot_range = np.logspace(-2, 2, 50)
N1, N2 = model.predict_normal_stresses(gamma_dot_range)

# First normal stress coefficient
Psi_1 = N1 / gamma_dot_range**2

# Check N₁/|N₂| ratio (theory: 1/(2α) - 1)
ratio = -N1[25] / N2[25]  # Midpoint
alpha_effective = 1 / (2 * (ratio + 1))
print(f"N₁/|N₂| ratio: {ratio:.2f}")
print(f"Implied alpha: {alpha_effective:.3f}")

Multi-Mode Fitting

from rheojax.models import GiesekusMultiMode

# Fit with 3 relaxation modes
model = GiesekusMultiMode(n_modes=3)
model.fit(omega, G_star, test_mode='oscillation')

# Extract mode parameters
for i in range(3):
    G_i = model.parameters.get_value(f'G_{i+1}')
    lambda_i = model.parameters.get_value(f'lambda_{i+1}')
    alpha_i = model.parameters.get_value(f'alpha_{i+1}')
    print(f"Mode {i+1}: G={G_i:.1f} Pa, λ={lambda_i:.3f} s, α={alpha_i:.3f}")

Note

Giesekus startup predictions require ODE integration of the conformation tensor \(\boldsymbol{\mu}\). First call compiles (~30s). Subsequent calls are fast (0.1-1s per curve).

Model Family 2: IKH/FIKH (Isotropic-Kinematic Hardening)

Physical Basis

IKH models describe thixotropic fluids with evolving microstructure:

  • Isotropic hardening: Structure parameter \(\lambda(t)\) modulates viscosity/modulus

  • Kinematic hardening: Backstress tensor \(\boldsymbol{\alpha}(t)\) tracks internal stress

  • Fractional variants (FIKH/FMLIKH): Power-law memory via Caputo derivative

Four model variants:

  1. MIKH (Maxwell Isotropic-Kinematic Hardening): Liquid-like, single relaxation time

  2. MLIKH (Maxwell Liquid Isotropic-Kinematic Hardening): Extended liquid formulation

  3. FIKH (Fractional Isotropic-Kinematic Hardening): MIKH + fractional derivative

  4. FMLIKH (Fractional Maxwell Liquid Isotropic-Kinematic Hardening): MLIKH + fractional derivative

Structure Evolution

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

  • \(\lambda = 1\): Fully structured (rest state)

  • \(\lambda = 0\): Fully broken (high shear)

  • \(t_{\text{eq}}\): Equilibrium buildup time

  • \(a, c\): Breakdown rate and exponent

Parameters and Protocols

MIKH/MLIKH parameters (5-6):

  • G_0: Modulus at full structure (Pa)

  • eta_inf: Viscosity at infinite shear (Pa·s)

  • lambda_0: Initial structure parameter (0-1)

  • a: Breakdown rate coefficient (s^c)

  • c: Breakdown exponent (0.5-2.0)

  • t_eq: Buildup time (s, MLIKH only)

FIKH/FMLIKH parameters (6-7):

  • All MIKH/MLIKH parameters plus:

  • alpha: Fractional order (0-1, dimensionless)

Supported protocols (6 each):

  • FLOW_CURVE, OSCILLATION, STARTUP, RELAXATION, CREEP, LAOS

Practical Implementation

Thixotropic Loop in Flow Curves

from rheojax.models import MIKHLocal

# Forward sweep: low → high shear rate
gamma_dot_forward = np.logspace(-2, 2, 30)

model = MIKHLocal()
model.parameters.set_value('lambda_0', 1.0)  # Start fully structured
stress_forward = model.predict(gamma_dot_forward, test_mode='flow_curve')

# Reverse sweep: high → low shear rate (requires re-fitting)
gamma_dot_reverse = gamma_dot_forward[::-1]
model.parameters.set_value('lambda_0', 0.5)  # Start partially broken
stress_reverse = model.predict(gamma_dot_reverse, test_mode='flow_curve')

# Hysteresis area
import matplotlib.pyplot as plt
plt.loglog(gamma_dot_forward, stress_forward, 'b-', label='Up sweep')
plt.loglog(gamma_dot_reverse, stress_reverse, 'r--', label='Down sweep')
plt.xlabel('Shear rate (s⁻¹)')
plt.ylabel('Stress (Pa)')
plt.legend()

Fractional Memory Effects

from rheojax.models import FMLIKHLocal

# Fit FMLIKH with fractional order α
model = FMLIKHLocal()
model.fit(t, G_t, test_mode='relaxation')

alpha = model.parameters.get_value('alpha')
print(f"Fractional order: {alpha:.3f}")

# Compare with integer-order MLIKH
from rheojax.models import MLIKHLocal
model_integer = MLIKHLocal()
model_integer.fit(t, G_t, test_mode='relaxation')

# Fractional model captures power-law relaxation better
# α ≈ 0.3-0.7 for complex fluids

Startup Overshoot and Structure Decay

# Simulate startup at constant shear rate
t = np.linspace(0.01, 100, 400)
gamma_dot = 1.0

# Get full ODE state (stress + structure)
result = model.simulate_startup(t, gamma_dot, return_full=True)
stress = result['stress']
lambda_t = result['lambda']

# Plot stress and structure evolution
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(t, stress)
ax1.set_ylabel('Stress (Pa)')
ax1.axhline(stress[-1], ls='--', color='gray', label='Steady state')

ax2.plot(t, lambda_t)
ax2.set_ylabel('Structure λ')
ax2.set_xlabel('Time (s)')
ax2.set_ylim([0, 1])

Key Insight

Fractional IKH models (FIKH/FMLIKH) add power-law memory to thixotropic evolution. The fractional order \(\alpha\) quantifies non-exponential structural relaxation, critical for soft matter with broad timescale distributions (clays, emulsions, blood).

Model Family 3: Fluidity-Saramito EVP

Physical Basis

Fluidity-Saramito models combine:

  • Saramito tensorial viscoelasticity: Upper-convected Maxwell (UCM) stress tensor

  • Von Mises yield criterion: \(\alpha_y = \max(0, 1 - \tau_y/|\boldsymbol{\tau}|)\)

  • Thixotropic fluidity evolution: \(f(t)\) controls both yield stress and relaxation time

Two coupling modes:

  1. Minimal coupling: Only relaxation time \(\lambda = 1/f\)

  2. Full coupling: Relaxation time + aging yield stress \(\tau_y(f) = \tau_{y,0} + k_{\text{age}} (1 - f)\)

Structure Evolution (Fluidity)

\[\frac{df}{dt} = -\underbrace{\frac{f - f_{\text{eq}}}{t_{\text{th}}}}_{\text{aging}} + \underbrace{b |\dot{\gamma}|^n}_{\text{rejuvenation}}\]

  • \(f = 1\): High fluidity (flowing)

  • \(f = 0\): Low fluidity (jammed/aging)

  • \(t_{\text{th}}\): Thixotropic timescale

  • \(b, n\): Rejuvenation rate and exponent

Parameters and Protocols

Local model parameters (7-8):

  • G: Elastic modulus (Pa)

  • k_d_0: Baseline relaxation rate (s⁻¹)

  • tau_y: Baseline yield stress (Pa)

  • f_eq: Equilibrium fluidity (0-1)

  • t_th: Thixotropic time (s)

  • b: Rejuvenation rate (s^(1-n))

  • n: Rejuvenation exponent (0.5-2.0)

  • k_age: Aging yield increment (Pa, full coupling only)

Nonlocal model adds:

  • n_points: Spatial discretization (51-201)

  • gap_width: Gap size (m)

  • D_f: Fluidity diffusion coefficient (m²/s)

Supported protocols (6):

  • FLOW_CURVE, STARTUP, CREEP, RELAXATION, OSCILLATION, LAOS

Practical Implementation

Yield Stress with Elasticity

from rheojax.models.fluidity.saramito import FluiditySaramitoLocal

# Minimal coupling (λ = 1/f only)
model = FluiditySaramitoLocal(coupling="minimal")

# Fit to flow curve
gamma_dot = np.logspace(-3, 2, 50)
sigma = ...  # Experimental stress (Pa)

model.fit(gamma_dot, sigma, test_mode='flow_curve')

# Extract yield stress
tau_y = model.parameters.get_value('tau_y')
print(f"Yield stress: {tau_y:.1f} Pa")

# Verify Herschel-Bulkley behavior at low shear rates
# σ = τ_y + K·γ̇^n

Creep with Viscosity Bifurcation

# Apply constant stress above/below yield
t = np.linspace(0, 500, 1000)

# Below yield: creep compliance saturates
sigma_low = 0.8 * tau_y
gamma_low, f_low = model.simulate_creep(t, sigma_low)

# Above yield: unbounded creep (viscous flow)
sigma_high = 1.2 * tau_y
gamma_high, f_high = model.simulate_creep(t, sigma_high)

# Plot bifurcation
plt.figure()
plt.plot(t, gamma_low, label=f'σ = {sigma_low:.1f} Pa (< τ_y)')
plt.plot(t, gamma_high, label=f'σ = {sigma_high:.1f} Pa (> τ_y)')
plt.xlabel('Time (s)')
plt.ylabel('Strain')
plt.legend()

Normal Stress from UCM Tensor

# Saramito uses tensorial stress → N₁ prediction
gamma_dot_range = np.logspace(-2, 2, 50)
N1, N2 = model.predict_normal_stresses(gamma_dot_range)

# UCM limit: N₁ = 2λG·γ̇², N₂ = 0
# (Saramito modifies λ → 1/f(γ̇))

print(f"First normal stress coefficient: {N1[25]/gamma_dot_range[25]**2:.1f} Pa·s²")

Shear Banding with Nonlocal Model

from rheojax.models.fluidity.saramito import FluiditySaramitoNonlocal

# Nonlocal model with structure diffusion
model_nl = FluiditySaramitoNonlocal(coupling="full", n_points=51)
model_nl.parameters.set_value('D_f', 1e-9)  # m²/s
model_nl.parameters.set_value('gap_width', 1e-3)  # 1 mm gap

# Predict steady-state profile
gamma_dot_avg = 10.0
result = model_nl.simulate_steady_shear(gamma_dot_avg, t_end=500)

# Extract spatial profiles
y_positions = result['y']
gamma_dot_profile = result['gamma_dot']
fluidity_profile = result['fluidity']

# Detect banding (high/low shear regions)
banding = model_nl.detect_banding(result, threshold=0.1)
if banding['has_banding']:
    print(f"Band interface at y = {banding['interface_position']:.4f} mm")
    print(f"Shear rate ratio: {banding['shear_ratio']:.1f}")

Warning

Nonlocal models solve PDEs on spatial grids (51-201 points). Memory scales as O(n_points × n_timesteps). Use n_points=51 for prototyping, increase to 101-201 only for publication-quality profiles.

Bayesian Inference for ODE Models

Workflow: NLSQ Warm-Start → NUTS

ODE models are computationally expensive for MCMC. Always use NLSQ to initialize:

from rheojax.models import GiesekusSingleMode

# Step 1: NLSQ point estimate (fast)
model = GiesekusSingleMode()
model.fit(omega, G_star, test_mode='oscillation', max_iter=1000)

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

# Step 3: Check convergence
import arviz as az
rhat = az.rhat(result.posterior_samples)
ess = az.ess(result.posterior_samples)

print(f"R-hat (should be <1.01): {rhat.max():.3f}")
print(f"ESS (should be >400): {ess.min():.0f}")

Memory Considerations

# For ODE models, reduce sample counts to avoid OOM

# Standard protocols (SAOS, flow curves): 1000 warmup + 2000 samples OK
result = model.fit_bayesian(omega, G_star, num_warmup=1000, num_samples=2000)

# Transient protocols (startup, LAOS): reduce to 25-50 samples in FAST_MODE
import os
if os.environ.get("FAST_MODE", "0") == "1":
    num_warmup, num_samples = 25, 50
else:
    num_warmup, num_samples = 500, 1000

result = model.fit_bayesian(t, stress, test_mode='startup',
                            num_warmup=num_warmup, num_samples=num_samples)

Prior Selection

# Tighten priors for ODE models (improves convergence)
from rheojax.core.parameter import Prior

model = MIKHLocal()

# Structure parameter: informative beta prior
model.parameters['lambda_0'].prior = Prior('beta', alpha=2.0, beta=2.0)

# Breakdown exponent: narrow around physical values
model.parameters['c'].prior = Prior('truncated_normal', loc=1.0, scale=0.2)

# Log-normal for positive-only parameters (G, eta_inf)
model.parameters['G_0'].prior = Prior('lognormal', loc=np.log(1000), scale=1.0)

Note

ODE models with >8 parameters may require tighter priors or reparameterization. Consider fixing less-sensitive parameters (e.g., c, n) based on literature values.

Visualization and Diagnostics

Stress Overshoot Visualization

from rheojax.visualization import plot_startup_overshoot

# Automatic detection of peak stress
fig, axes = plot_startup_overshoot(
    model,
    gamma_dot_values=[0.1, 1.0, 10.0],
    t_max=50,
    mark_overshoot=True
)
plt.savefig('giesekus_startup.png', dpi=300)

Normal Stress Comparison

from rheojax.visualization import plot_normal_stresses

fig, ax = plot_normal_stresses(
    model,
    gamma_dot_range=np.logspace(-2, 2, 50),
    plot_ratio=True  # Add N₁/|N₂| ratio
)

Thixotropic Loop

# Custom visualization for hysteresis
fig, ax = plt.subplots()

# Forward sweep
ax.loglog(gamma_dot_forward, stress_forward, 'b-o',
         label='Increasing $\dot{\gamma}$', markersize=4)

# Reverse sweep
ax.loglog(gamma_dot_reverse, stress_reverse, 'r--s',
         label='Decreasing $\dot{\gamma}$', markersize=4)

# Shade hysteresis region
ax.fill_between(gamma_dot_forward, stress_forward, stress_reverse,
                alpha=0.2, color='gray')

ax.set_xlabel('Shear rate $\dot{\gamma}$ (s$^{-1}$)')
ax.set_ylabel('Stress $\sigma$ (Pa)')
ax.legend()
ax.grid(alpha=0.3)

Bayesian Posterior Plots

# Pair plot for parameter correlations
from rheojax.pipeline.bayesian import BayesianPipeline

pipeline = BayesianPipeline()
(pipeline.load_data(omega, G_star, test_mode='oscillation')
         .fit_nlsq('giesekus_single')
         .fit_bayesian(num_samples=2000)
         .plot_pair(divergences=True)  # Mark divergent samples
         .plot_forest(hdi_prob=0.95)   # Credible intervals
         .plot_trace()                  # MCMC chains
         .save_results('giesekus_bayes.hdf5'))

Model Selection and Comparison

Comparison Table

ODE-Based Model Selection Guide

Model

Best For

Key Feature

Parameters

ODE Size

Fit Time

Giesekus (single)

Polymer solutions

Shear thinning + \(N_1\)

3

5

~30s (startup)

Giesekus (multi)

Polydisperse polymers

Multi-mode relaxation

3N

5N

~60s (startup)

MIKH

Thixotropic fluids

Structure buildup/breakdown

5

6

~40s (startup)

MLIKH

Complex thixotropy

Extended liquid formulation

6

6

~40s (startup)

FIKH

Fractional thixotropy

Power-law memory

6

6

~50s (startup)

FMLIKH

Advanced thixotropy

Fractional + liquid

7

6

~50s (startup)

Saramito (local)

EVP fluids

Yield + thixotropy

7-8

4

~35s (startup)

Saramito (nonlocal)

Banding EVP

Spatial heterogeneity

8-9

4×N_pts

~300s (startup)

Performance Benchmarks

Typical fit times on CPU (Apple M1, 8 cores):

# NLSQ optimization (100 data points)
# Giesekus SAOS: 2-5s
# Giesekus startup: 20-30s
# IKH SAOS: 3-6s
# IKH startup: 30-40s
# Saramito SAOS: 2-4s
# Saramito startup: 25-35s

# Bayesian NUTS (1000 warmup + 2000 samples, 4 chains)
# Giesekus SAOS: 3-5 min
# Giesekus startup: 15-25 min
# IKH SAOS: 4-6 min
# IKH startup: 20-30 min
# Saramito SAOS: 3-5 min
# Saramito startup: 18-28 min

Key Insight

Startup/LAOS protocols are 5-10× slower than SAOS due to fine temporal resolution required for stress overshoot capture. Use adaptive timestepping (diffrax default) and compile models before batch processing.

Limitations and Pitfalls

Numerical Stability

  • Stiff ODEs at high shear rates: Use implicit solvers (Kvaerno3/5) or reduce rtol/atol

  • Tensor symmetry breaking: Giesekus/Saramito enforce \(\tau_{xy} = \tau_{yx}\) but errors accumulate

  • Negative moduli: Structure parameter \(\lambda \to 0\) can cause \(G(\lambda) < 0\) → add lower bound

Computational Cost

  • First call compilation: 30-90s for diffrax JIT → use model.precompile()

  • Bayesian OOM: LAOS with 1000 cycles × 4 chains × 2000 samples exceeds 16GB RAM → reduce num_samples

  • Batch processing: ODE models don’t vectorize well → use joblib parallelism over datasets

Physical Validity

  • Mobility bounds: Giesekus \(\alpha > 0.5\) violates thermodynamics → enforce alpha.bounds = (0, 0.5)

  • Fractional order: FIKH/FMLIKH with \(\alpha > 0.9\) behaves like integer model → check if fractional is needed

  • Yield stress artifacts: Von Mises criterion at \(|\tau| \approx \tau_y\) causes stress discontinuities

Warning

ODE models can fit noise if over-parameterized. Always compare to simpler analytical models (Maxwell, Carreau-Yasuda) and verify that overshoot/thixotropy are genuine physical effects, not measurement artifacts.

Tutorial Notebooks

Each model family has 6-12 tutorial notebooks in examples/:

Giesekus Notebooks

  • examples/giesekus/01_giesekus_flow_curve.ipynb: Steady shear viscosity

  • examples/giesekus/02_giesekus_saos.ipynb: Linear viscoelasticity (\(G'\), \(G''\))

  • examples/giesekus/03_giesekus_startup.ipynb: Stress overshoot prediction

  • examples/giesekus/04_giesekus_normal_stresses.ipynb: \(N_1\), \(N_2\) in shear flow

  • examples/giesekus/05_giesekus_creep.ipynb: Creep compliance

  • examples/giesekus/06_giesekus_relaxation.ipynb: Stress relaxation

  • examples/giesekus/07_giesekus_laos.ipynb: Nonlinear oscillatory response

IKH/FIKH Notebooks

  • examples/ikh/01-06_mikh_*.ipynb: MIKH model (6 protocols)

  • examples/ikh/07-12_mlikh_*.ipynb: MLIKH model (6 protocols)

  • examples/fikh/01-06_fikh_*.ipynb: FIKH fractional model (6 protocols)

  • examples/fikh/07-12_fmlikh_*.ipynb: FMLIKH fractional model (6 protocols)

Fluidity-Saramito Notebooks

  • examples/fluidity/13-18_saramito_local_*.ipynb: Local EVP (6 protocols)

  • examples/fluidity/19-24_saramito_nonlocal_*.ipynb: Nonlocal with banding (6 protocols)

Note

All notebooks include FAST_MODE flag for quick testing (reduced samples/points). Set FAST_MODE=0 for publication-quality fits with full Bayesian diagnostics.

References

Giesekus Model

  • Giesekus, H. (1982). “A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility.” Journal of Non-Newtonian Fluid Mechanics, 11(1-2), 69-109. https://doi.org/10.1016/0377-0257(82)85016-7

  • Bird, R. B., Armstrong, R. C., & Hassager, O. (1987). Dynamics of Polymeric Liquids, Vol. 1: Fluid Mechanics (2nd ed.). Wiley. ISBN: 978-0-471-80245-7

IKH Models

  • 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

  • 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

Fractional Extensions

  • Bagley, R. L., & Torvik, P. J. (1983). “A theoretical basis for the application of fractional calculus to viscoelasticity.” Journal of Rheology, 27(3), 201-210. https://doi.org/10.1122/1.549724

  • Palade, L. I., Verney, V., & Attané, P. (1996). “A modified fractional model to describe the entire viscoelastic behavior of polybutadienes from flow to glassy regime.” Rheologica Acta, 35(3), 265-273. https://doi.org/10.1007/BF00366913

Saramito EVP Model

Fluidity Approach

  • Derec, C., Ducouret, G., Ajdari, A., & Lequeux, F. (2003). “Aging and nonlinear rheology in suspensions of polyethylene oxide–protected silica particles.” Physical Review E, 67(6), 061403. https://doi.org/10.1103/PhysRevE.67.061403

  • Moorcroft, R. L., & Fielding, S. M. (2013). “Criteria for shear banding in time-dependent flows of complex fluids.” Physical Review Letters, 110(8), 086001. https://doi.org/10.1103/PhysRevLett.110.086001

See Also

Related tutorials in this documentation:

Model handbooks:

Performance optimization:

  • /user_guide/03_advanced_topics/performance_optimization — JAX compilation, batching

  • /user_guide/04_reference/api_core — Core API reference