Models API ¶

extract_harmonics(stress, n_points_per_cycle)[source]¶

Extract Fourier harmonics from LAOS stress signal.

For linear model, only I_1 is non-zero.

Parameters:

stress (ndarray) – Stress signal from last cycle
n_points_per_cycle (int) – Points per cycle

Return type:

dict[str, float]

Returns:

Dictionary with I_1, I_3, I_3_I_1 (I_3/I_1 ratio)

Fractional Models¶

Fractional Maxwell Gel¶

rheojax.models.fractional.fractional_maxwell_gel.FractionalMaxwellGel | Handbook: Fractional Maxwell Gel (Fractional) Spring in series with SpringPot for gel-like plateaus.

Parameters¶
Parameter	Units	Bounds	Description
`c_alpha`	Pa·s^α	[0.001, 1e+09]	SpringPot material constant
`alpha`	dimensionless	[0, 1]	Power-law exponent
`eta`	Pa·s	[1e-06, 1e+12]	Dashpot viscosity

class rheojax.models.fractional.fractional_maxwell_gel.FractionalMaxwellGel[source]¶

Bases: BaseModel

Fractional Maxwell Gel model: SpringPot in series with dashpot.

This model describes the rheology of materials transitioning from power-law viscoelastic behavior to terminal flow, such as polymer solutions and gels.

parameters¶: ParameterSet with c_alpha, alpha, eta

Examples

>>> from rheojax.models import FractionalMaxwellGel
>>> from rheojax.core.data import RheoData
>>> import numpy as np
>>>
>>> # Create model with parameters
>>> model = FractionalMaxwellGel()
>>> model.parameters.set_value('c_alpha', 1e5)
>>> model.parameters.set_value('alpha', 0.5)
>>> model.parameters.set_value('eta', 1e3)
>>>
>>> # Predict relaxation modulus
>>> t = np.logspace(-3, 3, 50)
>>> data = RheoData(x=t, y=np.zeros_like(t), domain='time')
>>> data.metadata['test_mode'] = 'relaxation'
>>> G_t = model.predict(data)
>>>
>>> # Predict complex modulus
>>> omega = np.logspace(-2, 2, 50)
>>> data = RheoData(x=omega, y=np.zeros_like(omega), domain='frequency')
>>> G_star = model.predict(data)

__init__()[source]¶: Initialize Fractional Maxwell Gel model.

bayesian_nuts_kwargs()[source]¶

Prefer conservative NUTS settings for the stiff Mittag-Leffler kernel.

Return type:: dict

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [c_alpha, alpha, eta]

Returns:

Model predictions as JAX array

predict_rheodata(rheo_data, test_mode=None)[source]¶

Predict response for RheoData.

Parameters:

rheo_data (RheoData) – Input RheoData with x values
test_mode (str | None) – Test mode (‘relaxation’, ‘creep’, ‘oscillation’) If None, auto-detect from rheo_data

Return type:

Returns:

RheoData with predicted y values

predict(X, test_mode=None, **kwargs)[source]¶

Predict response.

R13-FMG-002: Delegates to BaseModel.predict() for plain-array inputs so that deformation_mode (E*->G* conversion) and test_mode restoration are handled correctly. Only RheoData inputs bypass super() because they return a RheoData wrapper object.

Parameters:

X – RheoData object or array of x-values
test_mode (str | None) – Test mode for prediction (‘relaxation’, ‘creep’, ‘oscillation’) Required when X is a raw array. If None, defaults to ‘relaxation’.
**kwargs – Additional arguments passed to BaseModel.predict()

Returns:

Predicted values (RheoData if input is RheoData, else array)

Fractional Maxwell Liquid¶

rheojax.models.fractional.fractional_maxwell_liquid.FractionalMaxwellLiquid | Handbook: Fractional Maxwell Liquid (Fractional) SpringPot plus dashpot for liquid-like systems with memory.

Parameters¶
Parameter	Units	Bounds	Description
`Gm`	Pa	[0.001, 1e+09]	Maxwell modulus
`alpha`	dimensionless	[0, 1]	Power-law exponent
`tau_alpha`	s^alpha	[1e-06, 1e+06]	Relaxation time

class rheojax.models.fractional.fractional_maxwell_liquid.FractionalMaxwellLiquid[source]¶

Bases: BaseModel

Fractional Maxwell Liquid model: Spring in series with SpringPot.

This model describes materials with elastic response at short times and power-law relaxation at long times, such as polymer melts.

parameters¶: ParameterSet with Gm, alpha, tau_alpha

Examples

>>> from rheojax.models import FractionalMaxwellLiquid
>>> from rheojax.core.data import RheoData
>>> import numpy as np
>>>
>>> # Create model with parameters
>>> model = FractionalMaxwellLiquid()
>>> model.parameters.set_value('Gm', 1e6)
>>> model.parameters.set_value('alpha', 0.7)
>>> model.parameters.set_value('tau_alpha', 1.0)
>>>
>>> # Predict relaxation modulus
>>> t = np.logspace(-3, 3, 50)
>>> data = RheoData(x=t, y=np.zeros_like(t), domain='time')
>>> data.metadata['test_mode'] = 'relaxation'
>>> G_t = model.predict(data)

__init__()[source]¶: Initialize Fractional Maxwell Liquid model.

bayesian_prior_factory(param_name, lower, upper)[source]¶

Provide custom priors that stay near realistic data-informed scales.

bayesian_parameter_bounds(bounds, X, y, test_mode)[source]¶

Tighten tau bounds based on data scale to avoid pathological samples.

Return type:: dict[str, tuple[float | None, float | None]]

bayesian_nuts_kwargs()[source]¶

Prefer conservative NUTS settings for the stiff Mittag-Leffler kernel.

Return type:: dict

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [Gm, alpha, tau_alpha]

Returns:

Model predictions as JAX array

predict_rheodata(rheo_data, test_mode=None)[source]¶

Predict response for RheoData.

Parameters:

rheo_data (RheoData) – Input RheoData with x values
test_mode (str | None) – Test mode (‘relaxation’, ‘creep’, ‘oscillation’) If None, auto-detect from rheo_data

Return type:

Returns:

RheoData with predicted y values

predict(X, test_mode=None, **kwargs)[source]¶

Predict response.

R13-FML-002: Delegates to BaseModel.predict() for plain-array inputs so that deformation_mode (E*->G* conversion) and test_mode restoration are handled correctly. Only RheoData inputs bypass super() because they return a RheoData wrapper object.

Parameters:

X – RheoData object or array of x-values
test_mode (str | None) – Test mode for prediction (‘relaxation’, ‘creep’, ‘oscillation’, ‘flow_curve’). If None, defaults to ‘relaxation’.
**kwargs – Additional arguments passed to BaseModel.predict()

Returns:

Predicted values (RheoData if input is RheoData, else array)

Fractional Maxwell (General)¶

rheojax.models.fractional.fractional_maxwell_model.FractionalMaxwellModel | Handbook: Generalized Fractional Maxwell (Two-Order) Two SpringPots in series for broad relaxation spectra.

Parameters¶
Parameter	Units	Bounds	Description
`c1`	Pa·s^α	[0.001, 1e+09]	Material constant
`alpha`	dimensionless	[0, 1]	First fractional order
`beta`	dimensionless	[0, 1]	Second fractional order
`tau`	s	[1e-06, 1e+06]	Relaxation time

class rheojax.models.fractional.fractional_maxwell_model.FractionalMaxwellModel[source]¶

Bases: BaseModel

Fractional Maxwell Model: Two SpringPots in series with independent orders.

This is the most general fractional Maxwell model, allowing for complex viscoelastic behavior with two independent fractional orders.

parameters¶: ParameterSet with c1, alpha, beta, tau

Examples

>>> from rheojax.models import FractionalMaxwellModel
>>> from rheojax.core.data import RheoData
>>> import numpy as np
>>>
>>> # Create model with parameters
>>> model = FractionalMaxwellModel()
>>> model.parameters.set_value('c1', 1e5)
>>> model.parameters.set_value('alpha', 0.5)
>>> model.parameters.set_value('beta', 0.7)
>>> model.parameters.set_value('tau', 1.0)
>>>
>>> # Predict relaxation modulus
>>> t = np.logspace(-3, 3, 50)
>>> data = RheoData(x=t, y=np.zeros_like(t), domain='time')
>>> data.metadata['test_mode'] = 'relaxation'
>>> G_t = model.predict(data)

__init__()[source]¶: Initialize Fractional Maxwell Model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

CRITICAL: test_mode is now passed as parameter (NOT read from self._test_mode) to ensure correct posteriors in Bayesian inference (v0.4.0 fix).

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [c1, alpha, beta, tau]
test_mode – Explicit test mode for predictions. If None, defaults to ‘relaxation’ for backward compatibility.

Returns:

Model predictions as JAX array

predict_rheodata(rheo_data, test_mode=None)[source]¶

Predict response for RheoData.

Parameters:

rheo_data (RheoData) – Input RheoData with x values
test_mode (str | None) – Test mode (‘relaxation’, ‘creep’, ‘oscillation’) If None, auto-detect from rheo_data

Return type:

Returns:

RheoData with predicted y values

predict(X, test_mode=None, **kwargs)[source]¶

Predict response.

R13-FMM-002: Delegates to BaseModel.predict() for plain-array inputs so that deformation_mode (E*->G* conversion) and test_mode restoration are handled correctly. Only RheoData inputs bypass super() because they return a RheoData wrapper object.

Parameters:

X – RheoData object or array of x-values
test_mode (str | None) – Test mode for prediction (‘relaxation’, ‘creep’, ‘oscillation’) Required when X is a raw array. If None, defaults to ‘relaxation’.
**kwargs – Additional arguments passed to BaseModel.predict()

Returns:

Predicted values (RheoData if input is RheoData, else array)

Fractional Kelvin-Voigt¶

rheojax.models.fractional.fractional_kelvin_voigt.FractionalKelvinVoigt | Handbook: Fractional Kelvin-Voigt (Fractional) Spring plus SpringPot in parallel for solid-like gels.

Parameters¶
Parameter	Units	Bounds	Description
`Ge`	Pa	[0.001, 1e+09]	Equilibrium modulus
`c_alpha`	Pa·s^α	[0.001, 1e+09]	SpringPot constant
`alpha`	dimensionless	[0, 1]	Fractional order

class rheojax.models.fractional.fractional_kelvin_voigt.FractionalKelvinVoigt[source]¶

Bases: BaseModel

Fractional Kelvin-Voigt model: Spring and SpringPot in parallel.

This model describes solid-like materials with power-law creep behavior, typical of filled polymers and soft solids.

parameters¶: ParameterSet with Ge, c_alpha, alpha

Examples

>>> from rheojax.models import FractionalKelvinVoigt
>>> from rheojax.core.data import RheoData
>>> import numpy as np
>>>
>>> # Create model with parameters
>>> model = FractionalKelvinVoigt()
>>> model.parameters.set_value('Ge', 1e6)
>>> model.parameters.set_value('c_alpha', 1e4)
>>> model.parameters.set_value('alpha', 0.5)
>>>
>>> # Predict relaxation modulus
>>> t = np.logspace(-3, 3, 50)
>>> data = RheoData(x=t, y=np.zeros_like(t), domain='time')
>>> data.metadata['test_mode'] = 'relaxation'
>>> G_t = model.predict(data)
>>>
>>> # Predict complex modulus
>>> omega = np.logspace(-2, 2, 50)
>>> data = RheoData(x=omega, y=np.zeros_like(omega), domain='frequency')
>>> G_star = model.predict(data)

__init__()[source]¶: Initialize Fractional Kelvin-Voigt model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [Ge, c_alpha, alpha]

Returns:

Model predictions as JAX array

predict_rheodata(rheo_data, test_mode=None)[source]¶

Predict response for RheoData.

Parameters:

rheo_data (RheoData) – Input RheoData with x values
test_mode (str | None) – Test mode (‘relaxation’, ‘creep’, ‘oscillation’) If None, auto-detect from rheo_data

Return type:

Returns:

RheoData with predicted y values

predict(X, test_mode=None, **kwargs)[source]¶

Predict response.

R13-FKV-002: Delegates to BaseModel.predict() for plain-array inputs so that deformation_mode (E*->G* conversion) and test_mode restoration are handled correctly. Only RheoData inputs bypass super() because they return a RheoData wrapper object.

Parameters:

X – RheoData object or array of x-values
test_mode (str | None) – Test mode for prediction (‘relaxation’, ‘creep’, ‘oscillation’) Required when X is a raw array. If None, defaults to ‘relaxation’.
**kwargs – Additional arguments passed to BaseModel.predict()

Returns:

Predicted values (RheoData if input is RheoData, else array)

Fractional Zener Solid-Liquid¶

rheojax.models.fractional.fractional_zener_sl.FractionalZenerSolidLiquid | Handbook: Fractional Zener Solid-Liquid (Fractional) Equilibrium spring parallel to fractional Maxwell arm.

Parameters¶
Parameter	Units	Bounds	Description
`Ge`	Pa	[0.001, 1e+09]	Equilibrium modulus
`c_alpha`	Pa·s^α	[0.001, 1e+09]	SpringPot constant
`alpha`	dimensionless	[0, 1]	Fractional order
`tau`	s	[1e-06, 1e+06]	Relaxation time

class rheojax.models.fractional.fractional_zener_sl.FractionalZenerSolidLiquid[source]¶

Bases: BaseModel

Fractional Zener Solid-Liquid model.

A fractional viscoelastic model combining equilibrium elasticity with fractional relaxation behavior.

Test Modes¶

Relaxation: Supported
Creep: Supported (via numerical inversion)
Oscillation: Supported
Rotation: Not supported (no steady-state flow)

Examples

>>> import jax.numpy as jnp
>>> from rheojax.models import FractionalZenerSolidLiquid
>>>
>>> # Create model
>>> model = FractionalZenerSolidLiquid()
>>>
>>> # Set parameters
>>> model.set_params(Ge=1000.0, c_alpha=500.0, alpha=0.5, tau=1.0)
>>>
>>> # Predict relaxation modulus
>>> t = jnp.logspace(-2, 2, 50)
>>> G_t = model.predict(t)  # Relaxation mode
>>>
>>> # Predict complex modulus
>>> omega = jnp.logspace(-2, 2, 50)
>>> G_star = model.predict(omega)  # Oscillation mode

__init__()[source]¶: Initialize Fractional Zener Solid-Liquid model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [Ge, c_alpha, alpha, tau]

Returns:

Model predictions as JAX array

Fractional Zener Solid-Solid¶

rheojax.models.fractional.fractional_zener_ss.FractionalZenerSolidSolid | Handbook: Fractional Zener Solid-Solid (Fractional) Two springs plus SpringPot for stiff gels and glasses.

Parameters¶
Parameter	Units	Bounds	Description
`Ge`	Pa	[0.001, 1e+09]	Equilibrium modulus
`Gm`	Pa	[0.001, 1e+09]	Maxwell arm modulus
`alpha`	dimensionless	[0, 1]	Fractional order
`tau_alpha`	s^alpha	[1e-06, 1e+06]	Relaxation time

class rheojax.models.fractional.fractional_zener_ss.FractionalZenerSolidSolid[source]¶

Bases: BaseModel

Fractional Zener Solid-Solid model.

A fractional viscoelastic model with both instantaneous and equilibrium elasticity.

Test Modes¶

Relaxation: Supported
Creep: Supported
Oscillation: Supported
Rotation: Not supported (no steady-state flow)

Examples

>>> import jax.numpy as jnp
>>> from rheojax.models import FractionalZenerSolidSolid
>>>
>>> # Create model
>>> model = FractionalZenerSolidSolid()
>>>
>>> # Set parameters
>>> model.set_params(Ge=1000.0, Gm=500.0, alpha=0.5, tau_alpha=1.0)
>>>
>>> # Predict relaxation modulus
>>> t = jnp.logspace(-2, 2, 50)
>>> G_t = model.predict(t)

__init__()[source]¶: Initialize Fractional Zener Solid-Solid model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [Ge, Gm, alpha, tau_alpha]

Returns:

Model predictions as JAX array

Fractional Zener Liquid-Liquid¶

rheojax.models.fractional.fractional_zener_ll.FractionalZenerLiquidLiquid | Handbook: Fractional Zener Liquid-Liquid (Fractional) Most general fractional Zener with three fractional orders.

Parameters¶
Parameter	Units	Bounds	Description
`c1`	Pa·s^α	[0.001, 1e+09]	First SpringPot constant
`c2`	Pa·s^γ	[0.001, 1e+09]	Second SpringPot constant
`alpha`	dimensionless	[0, 1]	First fractional order
`beta`	dimensionless	[0, 1]	Second fractional order
`gamma`	dimensionless	[0, 1]	Third fractional order
`tau`	s	[1e-06, 1e+06]	Relaxation time

class rheojax.models.fractional.fractional_zener_ll.FractionalZenerLiquidLiquid[source]¶

Bases: BaseModel

Fractional Zener Liquid-Liquid model.

The most general fractional Zener model with three independent fractional orders.

Test Modes¶

Relaxation: Supported (numerical)
Creep: Supported (numerical)
Oscillation: Supported (analytical)
Rotation: Partial support (power-law at high shear rates)

Examples

>>> import jax.numpy as jnp
>>> from rheojax.models import FractionalZenerLiquidLiquid
>>>
>>> # Create model
>>> model = FractionalZenerLiquidLiquid()
>>>
>>> # Set parameters
>>> model.set_params(c1=500.0, c2=100.0, alpha=0.5, beta=0.3, gamma=0.7, tau=1.0)
>>>
>>> # Predict complex modulus
>>> omega = jnp.logspace(-2, 2, 50)
>>> G_star = model.predict(omega)

__init__()[source]¶: Initialize Fractional Zener Liquid-Liquid model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [c1, c2, alpha, beta, gamma, tau]

Returns:

Model predictions as JAX array

Fractional Kelvin-Voigt Zener¶

rheojax.models.fractional.fractional_kv_zener.FractionalKelvinVoigtZener | Handbook: Fractional Kelvin-Voigt-Zener (Fractional) Series spring plus fractional Kelvin-Voigt block.

Parameters¶
Parameter	Units	Bounds	Description
`Ge`	Pa	[0.001, 1e+09]	Series spring modulus
`Gk`	Pa	[0.001, 1e+09]	KV element modulus
`alpha`	dimensionless	[0, 1]	Fractional order
`tau`	s	[1e-06, 1e+06]	Retardation time

class rheojax.models.fractional.fractional_kv_zener.FractionalKelvinVoigtZener[source]¶

Bases: BaseModel

Fractional Kelvin-Voigt Zener model.

A fractional viscoelastic model emphasizing retardation behavior with finite equilibrium compliance.

Test Modes¶

Relaxation: Supported (via inversion)
Creep: Supported (primary mode)
Oscillation: Supported
Rotation: Not supported (no steady-state flow)

Examples

>>> import jax.numpy as jnp
>>> from rheojax.models import FractionalKelvinVoigtZener
>>>
>>> # Create model
>>> model = FractionalKelvinVoigtZener()
>>>
>>> # Set parameters
>>> model.set_params(Ge=1000.0, Gk=500.0, alpha=0.5, tau=1.0)
>>>
>>> # Predict creep compliance
>>> t = jnp.logspace(-2, 2, 50)
>>> J_t = model.predict(t)

__init__()[source]¶: Initialize Fractional Kelvin-Voigt Zener model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [Ge, Gk, alpha, tau]

Returns:

Model predictions as JAX array

Fractional Burgers¶

rheojax.models.fractional.fractional_burgers.FractionalBurgersModel | Handbook: Fractional Burgers Model (Fractional) Combines Maxwell and Kelvin-Voigt elements for asphalt-like flows.

Parameters¶
Parameter	Units	Bounds	Description
`Jg`	1/Pa	[1e-09, 1000]	Glassy compliance
`eta1`	Pa·s	[1e-06, 1e+12]	Viscosity (Maxwell arm)
`Jk`	1/Pa	[1e-09, 1000]	Kelvin compliance
`alpha`	dimensionless	[0, 1]	Fractional order
`tau_k`	s	[1e-06, 1e+06]	Retardation time

class rheojax.models.fractional.fractional_burgers.FractionalBurgersModel[source]¶

Bases: BaseModel

Fractional Burgers model.

A four-parameter fractional viscoelastic model combining instantaneous compliance, viscous flow, and retardation.

Test Modes¶

Relaxation: Supported (via inversion)
Creep: Supported (primary mode)
Oscillation: Supported
Rotation: Partial support (viscous flow at low frequencies)

Examples

>>> import jax.numpy as jnp
>>> from rheojax.models import FractionalBurgersModel
>>>
>>> # Create model
>>> model = FractionalBurgersModel()
>>>
>>> # Set parameters
>>> model.set_params(Jg=1e-6, eta1=1000.0, Jk=5e-6, alpha=0.5, tau_k=1.0)
>>>
>>> # Predict creep compliance
>>> t = jnp.logspace(-2, 2, 50)
>>> J_t = model.predict(t)

__init__()[source]¶: Initialize Fractional Burgers model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [Jg, eta1, Jk, alpha, tau_k]

Returns:

Model predictions as JAX array

Fractional Poynting-Thomson¶

rheojax.models.fractional.fractional_poynting_thomson.FractionalPoyntingThomson | Handbook: Fractional Poynting-Thomson (Fractional) Fractional parallel branch with retardation time.

Parameters¶
Parameter	Units	Bounds	Description
`Ge`	Pa	[0.001, 1e+09]	Instantaneous modulus
`Gk`	Pa	[0.001, 1e+09]	Retarded modulus
`alpha`	dimensionless	[0, 1]	Fractional order
`tau`	s	[1e-06, 1e+06]	Retardation time

class rheojax.models.fractional.fractional_poynting_thomson.FractionalPoyntingThomson[source]¶

Bases: BaseModel

Fractional Poynting-Thomson model.

A fractional viscoelastic model emphasizing instantaneous elastic response with fractional retardation.

Test Modes¶

Relaxation: Supported
Creep: Supported (primary mode)
Oscillation: Supported
Rotation: Not supported (no steady-state flow)

Examples

>>> import jax.numpy as jnp
>>> from rheojax.models import FractionalPoyntingThomson
>>>
>>> # Create model
>>> model = FractionalPoyntingThomson()
>>>
>>> # Set parameters
>>> model.set_params(Ge=1500.0, Gk=500.0, alpha=0.5, tau=1.0)
>>>
>>> # Predict creep compliance
>>> t = jnp.logspace(-2, 2, 50)
>>> J_t = model.predict(t)

__init__()[source]¶: Initialize Fractional Poynting-Thomson model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [Ge, Gk, alpha, tau]

Returns:

Model predictions as JAX array

Fractional Jeffreys¶

rheojax.models.fractional.fractional_jeffreys.FractionalJeffreysModel | Handbook: Fractional Jeffreys Model (Fractional) Two dashpots plus SpringPot for thixotropic-like liquids.

Parameters¶
Parameter	Units	Bounds	Description
`eta1`	Pa·s	[1e-06, 1e+12]	First viscosity
`eta2`	Pa·s	[1e-06, 1e+12]	Second viscosity
`alpha`	dimensionless	[0, 1]	Fractional order
`tau1`	s	[1e-06, 1e+06]	Relaxation time

class rheojax.models.fractional.fractional_jeffreys.FractionalJeffreysModel[source]¶

Bases: BaseModel

Fractional Jeffreys model.

A fractional viscoelastic liquid model combining viscous flow with fractional relaxation behavior.

Test Modes¶

Relaxation: Supported
Creep: Supported
Oscillation: Supported
Rotation: Supported (viscous flow at low frequencies)

Examples

>>> import jax.numpy as jnp
>>> from rheojax.models import FractionalJeffreysModel
>>>
>>> # Create model
>>> model = FractionalJeffreysModel()
>>>
>>> # Set parameters
>>> model.set_params(eta1=1000.0, eta2=500.0, alpha=0.5, tau1=1.0)
>>>
>>> # Predict relaxation modulus
>>> t = jnp.logspace(-2, 2, 50)
>>> G_t = model.predict(t)

__init__()[source]¶: Initialize Fractional Jeffreys model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (time or frequency)
params – Array of parameter values [eta1, eta2, alpha, tau1]

Returns:

Model predictions as JAX array

Non-Newtonian Flow Models¶

Power Law¶

rheojax.models.flow.power_law.PowerLaw | Handbook: Power-Law (Ostwald–de Waele) Two-parameter viscosity curve tau = K*gamma^n.

Parameters¶
Parameter	Units	Bounds	Description
`K`	Pa·sⁿ	[1e-06, 1e+06]	Consistency index
`n`	dimensionless	[0.01, 2]	Flow behavior index

class rheojax.models.flow.power_law.PowerLaw[source]¶

Bases: BaseModel

Power Law model for non-Newtonian flow (ROTATION only).

The Power Law (Ostwald-de Waele) model is a simple empirical relationship that describes shear-thinning (n < 1) and shear-thickening (n > 1) behavior in steady shear flow.

Parameters:

K – Consistency index (Pa·s^n), controls viscosity magnitude
n – Flow behavior index (dimensionless), controls shear-thinning/thickening

Constitutive Equation:: σ(γ̇) = K |γ̇| ^n η(γ̇) = K |γ̇| ^(n-1)
Special Cases:: n = 1: Newtonian fluid with viscosity η = K n < 1: Shear-thinning (viscosity decreases with shear rate) n > 1: Shear-thickening (viscosity increases with shear rate)
Test Mode:: ROTATION (steady shear) only

__init__()[source]¶: Initialize Power Law model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (shear rate γ̇)
params – Array of parameter values [K, n]

Returns:

Model predictions as JAX array (viscosity η)

predict_stress(gamma_dot)[source]¶

Predict shear stress for given shear rates.

Parameters:: gamma_dot (ndarray) – Shear rate data (γ̇)
Return type:: ndarray
Returns:: Predicted shear stress σ(γ̇) = K |γ̇|^n

predict_rheo(rheo_data, test_mode=None, output='viscosity')[source]¶

Predict rheological response for RheoData.

Parameters:

rheo_data (RheoData) – Input rheological data
test_mode (TestModeEnum | None) – Test mode (must be ROTATION)
output (str) – Output type (‘viscosity’ or ‘stress’)

Return type:

Returns:

RheoData with predictions

Raises:

ValueError – If test mode is not ROTATION

__repr__()[source]¶

String representation.

Return type:: str

Carreau¶

rheojax.models.flow.carreau.Carreau | Handbook: Carreau Model Smooth transition from Newtonian plateau to shear thinning.

Parameters¶
Parameter	Units	Bounds	Description
`eta0`	Pa·s	[0.001, 1e+12]	Zero-shear viscosity
`eta_inf`	Pa·s	[1e-06, 1e+06]	Infinite-shear viscosity
`lambda_`	s	[1e-06, 1e+06]	Time constant
`n`	dimensionless	[0.01, 1]	Power-law index

class rheojax.models.flow.carreau.Carreau[source]¶

Bases: BaseModel

Carreau model for non-Newtonian flow (ROTATION only).

The Carreau model describes the smooth transition from a Newtonian plateau at low shear rates (zero-shear viscosity η_0) to power-law shear-thinning behavior at high shear rates, with an optional infinite-shear viscosity plateau.

Parameters:

eta0 – Zero-shear viscosity (Pa·s), Newtonian plateau at low γ̇
eta_inf – Infinite-shear viscosity (Pa·s), Newtonian plateau at high γ̇
lambda – Time constant (s), characteristic relaxation time
n – Power-law index (dimensionless), controls transition steepness

Constitutive Equation:: η(γ̇) = η_∞ + (η_0 - η_∞) [1 + (λ γ̇)²]^((n-1)/2)
Special Cases:: λ → 0: Newtonian fluid with η = η_0 λ → ∞, η_∞ = 0: Power Law behavior n = 1: Newtonian fluid for all shear rates
Test Mode:: ROTATION (steady shear) only

__init__()[source]¶: Initialize Carreau model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (shear rate γ̇)
params – Array of parameter values [eta0, eta_inf, lambda_, n]

Returns:

Model predictions as JAX array

predict_stress(gamma_dot)[source]¶

Predict shear stress for given shear rates.

Parameters:: gamma_dot (ndarray) – Shear rate data (γ̇)
Return type:: ndarray
Returns:: Predicted shear stress σ(γ̇)

predict_rheo(rheo_data, test_mode=None, output='viscosity')[source]¶

Predict rheological response for RheoData.

Parameters:

rheo_data (RheoData) – Input rheological data
test_mode (TestModeEnum | None) – Test mode (must be ROTATION)
output (str) – Output type (‘viscosity’ or ‘stress’)

Return type:

Returns:

RheoData with predictions

Raises:

ValueError – If test mode is not ROTATION

__repr__()[source]¶

String representation.

Return type:: str

Carreau-Yasuda¶

rheojax.models.flow.carreau_yasuda.CarreauYasuda | Handbook: Carreau–Yasuda Model Adds the Yasuda exponent to tune transition sharpness.

Parameters¶
Parameter	Units	Bounds	Description
`eta0`	Pa·s	[0.001, 1e+12]	Zero-shear viscosity
`eta_inf`	Pa·s	[1e-06, 1e+06]	Infinite-shear viscosity
`lambda_`	s	[1e-06, 1e+06]	Time constant
`n`	dimensionless	[0.01, 1]	Power-law index
`a`	dimensionless	[0.1, 2]	Transition parameter

class rheojax.models.flow.carreau_yasuda.CarreauYasuda[source]¶

Bases: BaseModel

Carreau-Yasuda model for non-Newtonian flow (ROTATION only).

The Carreau-Yasuda model extends the Carreau model with an additional parameter ‘a’ that provides better control over the transition region between Newtonian and power-law behavior. When a=2, it reduces to the standard Carreau model.

Parameters:

eta0 – Zero-shear viscosity (Pa·s), Newtonian plateau at low γ̇
eta_inf – Infinite-shear viscosity (Pa·s), Newtonian plateau at high γ̇
lambda – Time constant (s), characteristic relaxation time
n – Power-law index (dimensionless), controls asymptotic behavior
a – Transition parameter (dimensionless), controls transition width

Constitutive Equation:: η(γ̇) = η_∞ + (η_0 - η_∞) [1 + (λ γ̇)^a]^((n-1)/a)
Special Cases:: a = 2: Reduces to Carreau model λ → 0: Newtonian fluid with η = η_0 n = 1: Newtonian fluid for all shear rates
Test Mode:: ROTATION (steady shear) only

__init__()[source]¶: Initialize Carreau-Yasuda model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (shear rate γ̇)
params – Array of parameter values [eta0, eta_inf, lambda_, n, a]

Returns:

Model predictions as JAX array

predict_stress(gamma_dot)[source]¶

Predict shear stress for given shear rates.

Parameters:: gamma_dot (ndarray) – Shear rate data (γ̇)
Return type:: ndarray
Returns:: Predicted shear stress σ(γ̇)

predict_rheo(rheo_data, test_mode=None, output='viscosity')[source]¶

Predict rheological response for RheoData.

Parameters:

rheo_data (RheoData) – Input rheological data
test_mode (TestModeEnum | None) – Test mode (must be ROTATION)
output (str) – Output type (‘viscosity’ or ‘stress’)

Return type:

Returns:

RheoData with predictions

Raises:

ValueError – If test mode is not ROTATION

__repr__()[source]¶

String representation.

Return type:: str

Cross¶

rheojax.models.flow.cross.Cross | Handbook: Cross Model Alternative viscosity curve with rate exponent m.

Parameters¶
Parameter	Units	Bounds	Description
`eta0`	Pa·s	[0.001, 1e+12]	Zero-shear viscosity
`eta_inf`	Pa·s	[1e-06, 1e+06]	Infinite-shear viscosity
`lambda_`	s	[1e-06, 1e+06]	Time constant
`m`	dimensionless	[0.1, 2]	Rate constant

class rheojax.models.flow.cross.Cross[source]¶

Bases: BaseModel

Cross model for non-Newtonian flow (ROTATION only).

The Cross model describes the transition from a Newtonian plateau at low shear rates to power-law shear-thinning behavior at high shear rates. It uses a different functional form than the Carreau model and is often better suited for polymer solutions.

Parameters:

eta0 – Zero-shear viscosity (Pa·s), Newtonian plateau at low γ̇
eta_inf – Infinite-shear viscosity (Pa·s), Newtonian plateau at high γ̇
lambda – Time constant (s), characteristic relaxation time
m – Rate constant (dimensionless), controls transition steepness

Constitutive Equation:: η(γ̇) = η_∞ + (η_0 - η_∞) / [1 + (λγ̇)^m]
Special Cases:: λ → 0: Newtonian fluid with η = η_0 m → 0: Newtonian fluid for all shear rates λ → ∞: η approaches η_∞
Test Mode:: ROTATION (steady shear) only

__init__()[source]¶: Initialize Cross model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (shear rate γ̇)
params – Array of parameter values [eta0, eta_inf, lambda_, m]

Returns:

Model predictions as JAX array

predict_stress(gamma_dot)[source]¶

Predict shear stress for given shear rates.

Parameters:: gamma_dot (ndarray) – Shear rate data (γ̇)
Return type:: ndarray
Returns:: Predicted shear stress σ(γ̇)

predict_rheo(rheo_data, test_mode=None, output='viscosity')[source]¶

Predict rheological response for RheoData.

Parameters:

rheo_data (RheoData) – Input rheological data
test_mode (TestModeEnum | None) – Test mode (must be ROTATION)
output (str) – Output type (‘viscosity’ or ‘stress’)

Return type:

Returns:

RheoData with predictions

Raises:

ValueError – If test mode is not ROTATION

__repr__()[source]¶

String representation.

Return type:: str

Herschel-Bulkley¶

rheojax.models.flow.herschel_bulkley.HerschelBulkley | Handbook: Herschel-Bulkley Model Yield stress plus power-law viscosity.

Parameters¶
Parameter	Units	Bounds	Description
`sigma_y`	Pa	[0, 1e+06]	Yield stress
`K`	Pa·sⁿ	[1e-06, 1e+06]	Consistency index
`n`	dimensionless	[0.01, 2]	Flow behavior index

class rheojax.models.flow.herschel_bulkley.HerschelBulkley[source]¶

Bases: BaseModel

Herschel-Bulkley model for viscoplastic flow (ROTATION only).

The Herschel-Bulkley model describes materials that require a minimum stress (yield stress σ_y) to flow, and then exhibit power-law behavior. This is widely used for pastes, slurries, and suspensions.

Parameters:

sigma_y – Yield stress (Pa), minimum stress required for flow
K – Consistency index (Pa·s^n), controls viscosity magnitude
n – Flow behavior index (dimensionless), power-law exponent

Constitutive Equation:: σ(γ̇) = σ_y + K |γ̇| ^n for |σ| > σ_y γ̇ = 0 for |σ| ≤ σ_y
Special Cases:: σ_y = 0: Reduces to Power Law model n = 1: Reduces to Bingham model (linear viscoplastic) σ_y = 0, n = 1: Newtonian fluid
Test Mode:: ROTATION (steady shear) only

__init__()[source]¶: Initialize Herschel-Bulkley model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (shear rate γ̇)
params – Array of parameter values [sigma_y, K, n]

Returns:

Model predictions as JAX array

predict_viscosity(gamma_dot)[source]¶

Predict apparent viscosity for given shear rates.

Parameters:: gamma_dot (ndarray) – Shear rate data (γ̇)
Return type:: ndarray
Returns:: Predicted apparent viscosity η_app(γ̇)

predict_rheo(rheo_data, test_mode=None, output='stress')[source]¶

Predict rheological response for RheoData.

Parameters:

rheo_data (RheoData) – Input rheological data
test_mode (TestModeEnum | None) – Test mode (must be ROTATION)
output (str) – Output type (‘stress’ or ‘viscosity’)

Return type:

Returns:

RheoData with predictions

Raises:

ValueError – If test mode is not ROTATION

__repr__()[source]¶

String representation.

Return type:: str

Bingham¶

rheojax.models.flow.bingham.Bingham | Handbook: Bingham Plastic Yield stress and constant plastic viscosity.

Parameters¶
Parameter	Units	Bounds	Description
`sigma_y`	Pa	[0, 1e+06]	Yield stress
`eta_p`	Pa·s	[1e-06, 1e+12]	Plastic viscosity

class rheojax.models.flow.bingham.Bingham[source]¶

Bases: BaseModel

Bingham model for linear viscoplastic flow (ROTATION only).

The Bingham model describes materials that require a minimum stress (yield stress σ_y) to flow, and then exhibit constant (Newtonian) viscosity. This is the simplest viscoplastic model, used for materials like toothpaste, mayonnaise, and drilling mud.

Parameters:

sigma_y – Yield stress (Pa), minimum stress required for flow
eta_p – Plastic viscosity (Pa·s), constant viscosity above yield

Constitutive Equation:: σ(γ̇) = σ_y + η_p |γ̇| for |σ| > σ_y γ̇ = 0 for |σ| ≤ σ_y
Special Cases:: σ_y = 0: Reduces to Newtonian fluid with η = η_p This is a special case of Herschel-Bulkley with n = 1
Test Mode:: ROTATION (steady shear) only

__init__()[source]¶: Initialize Bingham model.

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

Model function for Bayesian inference.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes predictions given input X and a parameter array.

Parameters:

X – Independent variable (shear rate γ̇)
params – Array of parameter values [sigma_y, eta_p]

Returns:

Model predictions as JAX array (shear stress σ)

predict_viscosity(gamma_dot)[source]¶

Predict apparent viscosity for given shear rates.

Parameters:: gamma_dot (ndarray) – Shear rate data (γ̇)
Return type:: ndarray
Returns:: Predicted apparent viscosity η_app(γ̇)

predict_rheo(rheo_data, test_mode=None, output='stress')[source]¶

Predict rheological response for RheoData.

Parameters:

rheo_data (RheoData) – Input rheological data
test_mode (TestModeEnum | None) – Test mode (must be ROTATION)
output (str) – Output type (‘stress’ or ‘viscosity’)

Return type:

Returns:

RheoData with predictions

Raises:

ValueError – If test mode is not ROTATION

__repr__()[source]¶

String representation.

Return type:: str

Soft Glassy Rheology (SGR) Models¶

SGR Conventional¶

rheojax.models.sgr.sgr_conventional.SGRConventional | Handbook: SGR Conventional (Soft Glassy Rheology) — Handbook Statistical mechanics model for soft glassy materials (Sollich 1998).

Parameters¶
Parameter	Units	Bounds	Description
`x`	dimensionless	[0.5, 3.0]	Noise temperature (\(x < 1\): glass, \(1 < x < 2\): SGM, \(x \geq 2\): Newtonian)
`G0`	Pa	[0.001, 1e+09]	Characteristic modulus
`tau0`	s	[1e-09, 1e+06]	Microscopic attempt time

Material Classification by Noise Temperature:

\(x < 1\): Glass (aging, non-ergodic, yield stress)
\(1 < x < 2\): Soft Glassy Material (power-law rheology)
\(x \geq 2\): Newtonian liquid (no memory effects)

class rheojax.models.sgr.sgr_conventional.SGRConventional(dynamic_x=False)[source]¶

Bases: BaseModel

Soft Glassy Rheology (SGR) Conventional Model.

Statistical mechanics model for soft glassy materials (foams, emulsions, pastes, colloidal suspensions) based on Sollich 1998. The model describes rheological behavior through a population of mesoscopic elements trapped in energy wells with exponential distribution rho(E) ~ exp(-E), undergoing thermally activated hopping.

The effective noise temperature x controls material phase: - x < 1: Glass with yield stress (solid-like) - 1 < x < 2: Power-law viscoelastic fluid (G’ ~ G’’ ~ omega^(x-1)) - x >= 2: Newtonian liquid

Parameters:

x – Effective noise temperature (dimensionless), controls phase transition
G0 – Modulus scale (Pa), sets absolute magnitude of elastic response
tau0 – Attempt time (s), characteristic microscopic relaxation timescale

parameters¶: ParameterSet containing x, G0, tau0

Example

>>> from rheojax.models.sgr_conventional import SGRConventional
>>> import numpy as np
>>> model = SGRConventional()
>>> omega = np.logspace(-2, 2, 50)
>>> # Fitting oscillation data
>>> model.fit(omega, G_star, test_mode='oscillation')
>>> G_star_pred = model.predict(omega)
>>> # Bayesian inference
>>> result = model.fit_bayesian(omega, G_star, num_samples=2000)
>>> intervals = model.get_credible_intervals(result.posterior_samples)

Notes

Inherits from BaseModel (includes BayesianMixin for NumPyro NUTS)
Mode-aware Bayesian inference (stores test_mode for correct predictions)
JAX-accelerated kernel functions for GPU computation
Float64 precision critical for numerical stability near x=1

__init__(dynamic_x=False)[source]¶

Initialize SGR Conventional Model.

Creates ParameterSet with: - x (noise temperature): bounds (0.5, 3.0), default 1.5 - G0 (modulus scale): bounds (1e-3, 1e9), default 1e3 - tau0 (attempt time): bounds (1e-9, 1e3), default 1e-3

When dynamic_x=True, additional parameters for x(t) evolution: - x_eq: Equilibrium noise temperature at rest - alpha_aging: Aging rate (drives x toward x_eq) - beta_rejuv: Rejuvenation rate (shear-induced increase in x) - x_ss_A: Steady-state amplitude parameter - x_ss_n: Steady-state power-law exponent

Parameters:: dynamic_x (bool) – If True, x evolves via dx/dt equation. If False, x is constant.

The model is ready for fitting after instantiation.

get_phase_regime()[source]¶

Determine material phase regime from noise temperature x.

The SGR model exhibits three distinct phase regimes based on x: - Glass phase (x < 1): Solid-like with yield stress - Power-law fluid (1 <= x < 2): Viscoelastic with G’ ~ G’’ ~ omega^(x-1) - Newtonian liquid (x >= 2): Viscous with constant viscosity

Returns:: ‘glass’, ‘power-law’, or ‘newtonian’
Return type:: str

Example

>>> model = SGRConventional()
>>> model.parameters.set_value('x', 0.8)
>>> model.get_phase_regime()
'glass'

evolve_x(t, gamma_dot, x_initial)[source]¶

Evolve effective temperature x(t) via ODE integration.

Integrates the evolution equation:: dx/dt = -alpha_aging * (x - x_eq) + beta_rejuv * gamma_dot(t) * (x_ss(t) - x)

Uses JAX ODE integration (jax.experimental.ode.odeint) for stability and compatibility with JAX transformations.

Parameters:

t (ndarray) – Time array (s)
gamma_dot (ndarray) – Shear rate array (1/s), same length as t
x_initial (float) – Initial effective temperature

Returns:

Effective temperature x(t) at each time point

Return type:

Raises:

ValueError – If dynamic_x mode is not enabled or array shapes mismatch

Example

>>> model = SGRConventional(dynamic_x=True)
>>> t = np.linspace(0, 10, 100)
>>> gamma_dot = np.ones_like(t) * 5.0  # Constant shear
>>> x_t = model.evolve_x(t, gamma_dot, x_initial=1.0)

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

Model function for Bayesian inference with NumPyro NUTS.

This method is required by BayesianMixin for NumPyro NUTS sampling. It computes SGR predictions given input X and a parameter array.

Parameters:

X – Independent variable (frequency for oscillation, time for relaxation)
params – Array of parameter values [x, G0, tau0] Length: 3
test_mode – Optional test mode override. If None, uses stored self._test_mode
**kwargs – Protocol-specific arguments (gamma_dot, sigma_applied, etc.)

Returns:

Model predictions as JAX array - For oscillation: Complex modulus [G’, G’’] with shape (M, 2) - For relaxation: Relaxation modulus with shape (M,) - For creep: Creep compliance with shape (M,) - For steady_shear: Viscosity with shape (M,) - For startup: Stress growth coefficient with shape (M,)

Note

Uses stored test_mode from last fit() call for mode-aware inference. This ensures correct kernel functions are used during NUTS sampling.

Example

>>> # During Bayesian inference, NumPyro calls:
>>> predictions = model.model_function(omega, params_sample)
>>> # Where params_sample = [x_sample, G0_sample, tau0_sample]

simulate_laos(gamma_0, omega, n_cycles=2, n_points_per_cycle=256)[source]¶

Simulate LAOS response for given strain amplitude and frequency.

Generates time-domain stress response to sinusoidal strain input:: gamma(t) = gamma_0 * sin(omega * t)

For SGR model, the stress response is computed using the complex modulus in the linear viscoelastic approximation, with nonlinearity arising from strain-dependent softening at large amplitudes.

Parameters:

gamma_0 (float) – Strain amplitude (dimensionless)
omega (float) – Angular frequency (rad/s)
n_cycles (int) – Number of oscillation cycles to simulate
n_points_per_cycle (int) – Number of time points per cycle

Returns:

Strain array gamma(t) stress: Stress array sigma(t)

Return type:

Example

>>> model = SGRConventional()
>>> model.parameters.set_value("x", 1.5)
>>> strain, stress = model.simulate_laos(gamma_0=0.1, omega=1.0)

extract_laos_harmonics(stress, n_points_per_cycle=256)[source]¶

Extract Fourier harmonics from LAOS stress response.

Performs FFT analysis to extract harmonic amplitudes and phases:: sigma(t) = sum_n I_n * sin(n*omega*t + phi_n)

For LAOS, odd harmonics (n = 1, 3, 5, …) dominate due to symmetry.

Parameters:

stress (ndarray) – Stress time series from simulate_laos()
n_points_per_cycle (int) – Points per oscillation cycle

Returns:

I_1, I_3, I_5, …: Harmonic amplitudes
phi_1, phi_3, phi_5, …: Phase angles
I_3_I_1, I_5_I_1, …: Relative intensities

Return type:

Example

>>> strain, stress = model.simulate_laos(gamma_0=0.5, omega=1.0)
>>> harmonics = model.extract_laos_harmonics(stress)
>>> print(f"Third harmonic ratio: {harmonics['I_3_I_1']:.4f}")

compute_chebyshev_coefficients(strain, stress, gamma_0, omega, n_points_per_cycle=256)[source]¶

Compute Chebyshev decomposition of LAOS response.

Decomposes stress into elastic and viscous Chebyshev contributions:

sigma(gamma, gamma_dot) = sum_n e_n * T_n(gamma/gamma_0)

sum_n v_n * T_n(gamma_dot/gamma_dot_0)

where T_n are Chebyshev polynomials of the first kind.

Physical interpretation: - e_n: Elastic (in-phase with strain) Chebyshev coefficients - v_n: Viscous (out-of-phase with strain) Chebyshev coefficients - e_3/e_1 > 0: Strain stiffening - e_3/e_1 < 0: Strain softening - v_3/v_1 > 0: Shear thickening - v_3/v_1 < 0: Shear thinning

Parameters:

strain (ndarray) – Strain array from simulate_laos()
stress (ndarray) – Stress array from simulate_laos()
gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency
n_points_per_cycle (int) – Points per oscillation cycle

Returns:

e_1, e_3, e_5: Elastic Chebyshev coefficients
v_1, v_3, v_5: Viscous Chebyshev coefficients
e_3_e_1, v_3_v_1: Normalized coefficients

Return type:

Example

>>> strain, stress = model.simulate_laos(gamma_0=0.5, omega=1.0)
>>> chebyshev = model.compute_chebyshev_coefficients(
...     strain, stress, gamma_0=0.5, omega=1.0
... )
>>> print(f"Strain stiffening ratio: {chebyshev['e_3_e_1']:.4f}")

get_lissajous_curve(gamma_0, omega, n_points=256, normalized=False)[source]¶

Generate Lissajous curve (stress vs strain) for LAOS.

Parameters:

gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)
n_points (int) – Number of points in curve
normalized (bool) – If True, normalize strain and stress

Returns:

Strain array (one period) stress: Stress array (one period)

Return type:

Example

>>> strain, stress = model.get_lissajous_curve(gamma_0=0.1, omega=1.0)
>>> plt.plot(strain, stress)  # Elastic Lissajous

enable_thixotropy(k_build=0.1, k_break=0.5, n_struct=2.0)[source]¶

Enable thixotropy modeling with structural parameter lambda(t).

Adds thixotropy kinetics parameters to the model. The structural parameter lambda represents the state of internal microstructure: - lambda = 1: Fully built structure - lambda = 0: Fully broken structure

Evolution equation:: d(lambda)/dt = k_build * (1 - lambda) - k_break * gamma_dot * lambda
The effective modulus is coupled to lambda:: G_eff(t) = G0 * lambda(t)^n_struct

Parameters:

k_build (float) – Structure build-up rate (1/s), default 0.1
k_break (float) – Structure breakdown rate (dimensionless), default 0.5
n_struct (float) – Structural coupling exponent, default 2.0

Return type:

Example

>>> model = SGRConventional()
>>> model.enable_thixotropy(k_build=0.1, k_break=0.5, n_struct=2.0)
>>> # Now model can predict stress transients with thixotropy

evolve_lambda(t, gamma_dot, lambda_initial=1.0)[source]¶

Evolve structural parameter lambda(t) for given shear history.

Integrates the thixotropy kinetics equation:: d(lambda)/dt = k_build * (1 - lambda) - k_break * gamma_dot * lambda

Uses JAX lax.scan for vectorized time-stepping (replaces Python for-loop).

Parameters:

t (ndarray) – Time array (s)
gamma_dot (ndarray) – Shear rate array (1/s), same shape as t
lambda_initial (float) – Initial structural parameter [0, 1], default 1.0

Returns:

Structural parameter evolution, same shape as t

Return type:

Raises:

ValueError – If thixotropy not enabled or array shapes mismatch

Example

>>> model = SGRConventional()
>>> model.enable_thixotropy()
>>> t = np.linspace(0, 10, 100)
>>> gamma_dot = np.ones_like(t) * 10.0  # Constant shear
>>> lambda_t = model.evolve_lambda(t, gamma_dot, lambda_initial=1.0)

predict_thixotropic_stress(t, gamma_dot, lambda_t=None, lambda_initial=1.0)[source]¶

Predict stress response with thixotropic modulus.

The effective modulus is coupled to the structural parameter:: G_eff(t) = G0 * lambda(t)^n_struct

Parameters:

t (ndarray) – Time array (s)
gamma_dot (ndarray) – Shear rate array (1/s)
lambda_t (ndarray | None) – Pre-computed lambda trajectory, or None to compute
lambda_initial (float) – Initial lambda if computing [0, 1], default 1.0

Returns:

Stress response (Pa)

Return type:

Example

>>> model = SGRConventional()
>>> model.enable_thixotropy()
>>> t = np.linspace(0, 10, 100)
>>> gamma_dot = np.ones_like(t) * 10.0
>>> sigma = model.predict_thixotropic_stress(t, gamma_dot)

predict_stress_transient(t, gamma_dot, lambda_initial=1.0)[source]¶

Predict stress transient (overshoot/undershoot) for shear step protocol.

For step-up in shear rate: Initially high stress (intact structure) followed by decay as structure breaks down (overshoot).

For step-down in shear rate: Initially low stress (broken structure) followed by increase as structure rebuilds (undershoot).

Parameters:

t (ndarray) – Time array (s)
gamma_dot (ndarray) – Shear rate array (1/s), can include steps
lambda_initial (float) – Initial structural parameter [0, 1]

Returns:

Stress response (Pa) lambda_t: Structural parameter evolution

Return type:

Example

>>> model = SGRConventional()
>>> model.enable_thixotropy()
>>> t = np.linspace(0, 10, 100)
>>> gamma_dot = np.ones_like(t)
>>> gamma_dot[t >= 5] = 10.0  # Step up at t=5
>>> sigma, lambda_t = model.predict_stress_transient(t, gamma_dot)

detect_shear_banding(gamma_dot=None, sigma=None, n_points=100, gamma_dot_range=(0.01, 100.0))[source]¶

Detect shear banding from constitutive curve.

Computes the steady-state flow curve and checks for non-monotonicity (d sigma / d gamma_dot < 0) which indicates shear banding instability.

Parameters:

gamma_dot (ndarray | None) – Shear rate array (1/s). If None, uses gamma_dot_range.
sigma (ndarray | None) – Stress array (Pa). If None, computes from model.
n_points (int) – Number of points if computing flow curve
gamma_dot_range (tuple[float, float]) – Range for computing flow curve if gamma_dot is None

Returns:

True if shear banding detected banding_info: Dict with banding region info, or None

Return type:

tuple[bool, dict | None]

Example

>>> model = SGRConventional()
>>> model.parameters.set_value("x", 0.8)  # Glass regime
>>> is_banding, info = model.detect_shear_banding()

predict_banded_flow(gamma_dot_applied, gamma_dot=None, sigma=None, n_points=100)[source]¶

Predict flow in shear banding regime with lever rule.

When shear banding occurs, the material splits into bands with different local shear rates. This method computes the band fractions and the composite stress.

Parameters:

gamma_dot_applied (float) – Applied average shear rate (1/s)
gamma_dot (ndarray | None) – Shear rate array for flow curve. If None, computed.
sigma (ndarray | None) – Stress array for flow curve. If None, computed.
n_points (int) – Number of points if computing flow curve

Returns:

Dict with band coexistence info, or None

Return type:

dict | None

Example

>>> model = SGRConventional()
>>> model.parameters.set_value("x", 0.8)
>>> coex = model.predict_banded_flow(gamma_dot_applied=1.0)
>>> if coex:
...     print(f"Low band: {coex['fraction_low']:.2%}")
...     print(f"High band: {coex['fraction_high']:.2%}")

SGR Generic (GENERIC Framework)¶

rheojax.models.sgr.sgr_generic.SGRGeneric | Handbook: SGR GENERIC (Thermodynamically Consistent) Thermodynamically consistent SGR using GENERIC framework (Fuereder & Ilg 2013).

Parameters¶
Parameter	Units	Bounds	Description
`x`	dimensionless	[0.5, 3.0]	Noise temperature
`G0`	Pa	[0.001, 1e+09]	Characteristic modulus
`tau0`	s	[1e-09, 1e+06]	Microscopic attempt time

Advantages over Conventional SGR:

Satisfies Onsager reciprocal relations
Enhanced numerical stability near glass transition (x → 1)
Consistent thermodynamic framework

class rheojax.models.sgr.sgr_generic.SGRGeneric(dynamic_x=False)[source]¶

Bases: BaseModel

Soft Glassy Rheology (SGR) GENERIC Thermodynamic Framework Model.

This model implements the GENERIC (General Equation for Non-Equilibrium Reversible-Irreversible Coupling) thermodynamic framework for SGR, ensuring thermodynamic consistency via explicit entropy production tracking.

The GENERIC formulation splits dynamics into: - Reversible (Hamiltonian): dz/dt = L * dF/dz (Poisson bracket L antisymmetric) - Irreversible (dissipative): dz/dt = M * dS/dz (friction M symmetric PSD)

Parameters:

x – Effective noise temperature (dimensionless), controls phase transition
G0 – Modulus scale (Pa), sets absolute magnitude of elastic response
tau0 – Attempt time (s), characteristic microscopic relaxation timescale

State Variables:

z = [sigma, lambda] where: - sigma: Macroscopic stress (Pa) - lambda: Structural parameter [0, 1] representing trap occupation

Thermodynamic Functions:

F(z): Helmholtz free energy = U(z) - T*S(z)
U(z): Internal energy from elastic storage
S(z): Entropy from trap distribution
W: Entropy production rate = (dF/dz)^T M (dF/dz) >= 0

Example

>>> from rheojax.models.sgr_generic import SGRGeneric
>>> import numpy as np
>>> model = SGRGeneric()
>>> omega = np.logspace(-2, 2, 50)
>>> model._test_mode = 'oscillation'
>>> G_star = model.predict(omega)
>>> # Check thermodynamic consistency
>>> state = np.array([100.0, 0.5])
>>> W = model.compute_entropy_production(state)
>>> assert W >= 0, "Second law violated!"

Notes

Inherits from BaseModel (includes BayesianMixin for NumPyro NUTS)
Predictions match SGRConventional in linear viscoelastic regime
GENERIC structure guarantees thermodynamic consistency
Reference: Fuereder & Ilg 2013 PRE 88, 042134

__init__(dynamic_x=False)[source]¶

Initialize SGR GENERIC Model.

Creates ParameterSet with: - x (noise temperature): bounds (0.5, 3.0), default 1.5 - G0 (modulus scale): bounds (1e-3, 1e9), default 1e3 - tau0 (attempt time): bounds (1e-9, 1e3), default 1e-3

Parameters:: dynamic_x (bool) – If True, enable dynamic noise temperature evolution with 3D state [sigma, lambda, x]. Default False for backward compatibility.

free_energy(state)[source]¶

Compute Helmholtz free energy F(z) = U(z) - T*S(z).

The free energy functional for SGR combines: - Elastic energy storage from stressed trap elements - Entropic contribution from trap occupation distribution

Parameters:: state (ndarray) – State vector [sigma, lambda] where sigma is stress (Pa) and lambda is structural parameter [0, 1]
Return type:: float
Returns:: Free energy F (J/m^3 or Pa, depending on normalization)

Notes

F = U - T*S where T is the noise temperature x (in units of trap depth)

internal_energy(state)[source]¶

Compute internal energy U(z) from elastic storage.

The internal energy represents energy stored in elastically deformed trap elements. For SGR with stress sigma and structural parameter lambda:

U = (1/2) * (sigma^2 / (G0 * lambda^n))

where the effective modulus depends on structure.

Parameters:: state (ndarray) – State vector [sigma, lambda]
Return type:: float
Returns:: Internal energy U (J/m^3)

entropy(state)[source]¶

Compute entropy S(z) from trap occupation distribution.

The entropy represents the configurational entropy of trap occupation. For the structural parameter lambda in [0, 1], we use a mixing entropy form:

S = -k * [lambda * ln(lambda) + (1-lambda) * ln(1-lambda)]

This captures the entropy associated with the distribution of elements between trapped (structured) and free (unstructured) states.

Parameters:: state (ndarray) – State vector [sigma, lambda]
Return type:: float
Returns:: Entropy S (dimensionless, normalized by kB)

poisson_bracket(state)[source]¶

Compute Poisson bracket operator L(z).

The Poisson bracket generates reversible (Hamiltonian) dynamics. It must be antisymmetric: L = -L^T.

For SGR, the Poisson bracket couples stress sigma to strain rate:: L = [[0, L_12], [-L_12, 0]]

where L_12 encodes the stress-strain rate coupling from the constitutive relation.

Parameters:: state (ndarray) – State vector [sigma, lambda]
Return type:: ndarray
Returns:: 2x2 antisymmetric Poisson bracket matrix L

Notes

L is state-dependent in general
Antisymmetry ensures energy conservation: dE/dt = 0 for reversible part

friction_matrix(state)[source]¶

Compute friction matrix M(z).

The friction matrix generates irreversible (dissipative) dynamics. It must be symmetric and positive semi-definite: M = M^T, M >= 0.

For SGR, the friction matrix encodes: - Viscous dissipation (stress relaxation) - Structural evolution (trap hopping)

Parameters:: state (ndarray) – State vector [sigma, lambda]
Return type:: ndarray
Returns:: 2x2 symmetric positive semi-definite friction matrix M

Notes

M is state-dependent
Positive semi-definiteness ensures entropy production W >= 0
The noise temperature x appears in M controlling dissipation rate

reversible_dynamics(state)[source]¶

Compute reversible (Hamiltonian) part of dynamics.

dz/dt|_rev = L(z) * dF/dz

This represents the energy-conserving part of the dynamics, encoding the reversible coupling between variables.

Parameters:: state (ndarray) – State vector [sigma, lambda]
Return type:: ndarray
Returns:: Time derivative dz/dt from reversible dynamics

irreversible_dynamics(state)[source]¶

Compute irreversible (dissipative) part of dynamics.

dz/dt|_irrev = M(z) * dS/dz

where dS/dz = (1/T) * dF/dz for systems at effective temperature T.

This represents the entropy-producing part of the dynamics, encoding irreversible relaxation processes.

Parameters:: state (ndarray) – State vector [sigma, lambda]
Return type:: ndarray
Returns:: Time derivative dz/dt from irreversible dynamics

full_dynamics(state)[source]¶

Compute full GENERIC dynamics.

dz/dt = L(z) * dF/dz + M(z) * dS/dz

The total dynamics combines reversible (Hamiltonian) and irreversible (dissipative) contributions.

Parameters:: state (ndarray) – State vector [sigma, lambda]
Return type:: ndarray
Returns:: Total time derivative dz/dt

entropy_production_rate(state)[source]¶

Compute entropy production rate dS/dt.

This is equivalent to compute_entropy_production() but expressed in terms of entropy rather than free energy.

Parameters:: state (ndarray) – State vector [sigma, lambda]
Return type:: float
Returns:: Entropy production rate dS/dt >= 0

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

Model function for Bayesian inference with NumPyro NUTS.

Required by BayesianMixin for NumPyro NUTS sampling.

Parameters:

X – Independent variable (frequency or time)
params – Array of parameter values [x, G0, tau0]
test_mode – Optional test mode override
**kwargs – Protocol-specific arguments (gamma_dot, sigma_applied, etc.)

Returns:

Model predictions as JAX array

get_phase_regime()[source]¶

Determine material phase regime from noise temperature x.

Returns:: ‘glass’, ‘power-law’, or ‘newtonian’
Return type:: str

enable_thixotropy(k_build=0.1, k_break=0.5, n_struct=2.0)[source]¶

Enable thixotropy modeling with structural parameter lambda(t).

Evolution equation:: d(lambda)/dt = k_build * (1 - lambda) - k_break * gamma_dot * lambda
The effective modulus is coupled to lambda:: G_eff(t) = G0 * lambda(t)^n_struct

Parameters:

k_build (float) – Structure build-up rate (1/s), default 0.1
k_break (float) – Structure breakdown rate (dimensionless), default 0.5
n_struct (float) – Structural coupling exponent, default 2.0

Return type:

Example

>>> model = SGRGeneric()
>>> model.enable_thixotropy(k_build=0.1, k_break=0.5, n_struct=2.0)
>>> # Now model can predict stress transients with thixotropy

evolve_lambda(t, gamma_dot, lambda_initial=1.0)[source]¶

Evolve structural parameter lambda(t) for given shear history.

Integrates the thixotropy kinetics equation:: d(lambda)/dt = k_build * (1 - lambda) - k_break * gamma_dot * lambda

Uses JAX lax.scan for vectorized time-stepping (replaces Python for-loop).

Parameters:

t (ndarray) – Time array (s)
gamma_dot (ndarray) – Shear rate array (1/s), same shape as t
lambda_initial (float) – Initial structural parameter [0, 1], default 1.0

Returns:

Structural parameter evolution, same shape as t

Return type:

Raises:

ValueError – If thixotropy not enabled or array shapes mismatch

Example

>>> model = SGRGeneric()
>>> model.enable_thixotropy()
>>> t = np.linspace(0, 10, 100)
>>> gamma_dot = np.ones_like(t) * 10.0  # Constant shear
>>> lambda_t = model.evolve_lambda(t, gamma_dot, lambda_initial=1.0)

predict_thixotropic_stress(t, gamma_dot, lambda_t=None, lambda_initial=1.0)[source]¶

Predict stress response with thixotropic modulus.

The effective modulus is coupled to the structural parameter:: G_eff(t) = G0 * lambda(t)^n_struct

Parameters:

t (ndarray) – Time array (s)
gamma_dot (ndarray) – Shear rate array (1/s)
lambda_t (ndarray | None) – Pre-computed lambda trajectory, or None to compute
lambda_initial (float) – Initial lambda if computing [0, 1], default 1.0

Returns:

Stress response (Pa)

Return type:

Example

>>> model = SGRGeneric()
>>> model.enable_thixotropy()
>>> t = np.linspace(0, 10, 100)
>>> gamma_dot = np.ones_like(t) * 10.0
>>> sigma = model.predict_thixotropic_stress(t, gamma_dot)

predict_stress_transient(t, gamma_dot, lambda_initial=1.0)[source]¶

Predict stress transient (overshoot/undershoot) for shear step protocol.

For step-up in shear rate: Initially high stress (intact structure) followed by decay as structure breaks down (overshoot).

For step-down in shear rate: Initially low stress (broken structure) followed by increase as structure rebuilds (undershoot).

Parameters:

t (ndarray) – Time array (s)
gamma_dot (ndarray) – Shear rate array (1/s), can include steps
lambda_initial (float) – Initial structural parameter [0, 1]

Returns:

Stress response (Pa) lambda_t: Structural parameter evolution

Return type:

Example

>>> model = SGRGeneric()
>>> model.enable_thixotropy()
>>> t = np.linspace(0, 10, 100)
>>> gamma_dot = np.ones_like(t)
>>> gamma_dot[t >= 5] = 10.0  # Step up at t=5
>>> sigma, lambda_t = model.predict_stress_transient(t, gamma_dot)

detect_shear_banding(gamma_dot=None, sigma=None, n_points=100, gamma_dot_range=(0.01, 100.0))[source]¶

Detect shear banding from constitutive curve.

Computes the steady-state flow curve and checks for non-monotonicity (d sigma / d gamma_dot < 0) which indicates shear banding instability.

Parameters:

gamma_dot (ndarray | None) – Shear rate array (1/s). If None, uses gamma_dot_range.
sigma (ndarray | None) – Stress array (Pa). If None, computes from model.
n_points (int) – Number of points if computing flow curve
gamma_dot_range (tuple[float, float]) – Range for computing flow curve if gamma_dot is None

Returns:

True if shear banding detected banding_info: Dict with banding region info, or None

Return type:

tuple[bool, dict | None]

Example

>>> model = SGRGeneric()
>>> model.parameters.set_value("x", 0.8)  # Glass regime
>>> is_banding, info = model.detect_shear_banding()

predict_banded_flow(gamma_dot_applied, gamma_dot=None, sigma=None, n_points=100)[source]¶

Predict flow in shear banding regime with lever rule.

When shear banding occurs, the material splits into bands with different local shear rates. This method computes the band fractions and the composite stress.

Parameters:

gamma_dot_applied (float) – Applied average shear rate (1/s)
gamma_dot (ndarray | None) – Shear rate array for flow curve. If None, computed.
sigma (ndarray | None) – Stress array for flow curve. If None, computed.
n_points (int) – Number of points if computing flow curve

Returns:

Dict with band coexistence info, or None

Return type:

dict | None

Example

>>> model = SGRGeneric()
>>> model.parameters.set_value("x", 0.8)
>>> coex = model.predict_banded_flow(gamma_dot_applied=1.0)
>>> if coex:
...     print(f"Low band: {coex['fraction_low']:.2%}")
...     print(f"High band: {coex['fraction_high']:.2%}")

simulate_laos(gamma_0, omega, n_cycles=2, n_points_per_cycle=256)[source]¶

Simulate LAOS response for given strain amplitude and frequency.

Generates time-domain stress response to sinusoidal strain input:: gamma(t) = gamma_0 * sin(omega * t)

For SGR model, the stress response is computed using the complex modulus in the linear viscoelastic approximation, with nonlinearity arising from strain-dependent softening at large amplitudes.

Parameters:

gamma_0 (float) – Strain amplitude (dimensionless)
omega (float) – Angular frequency (rad/s)
n_cycles (int) – Number of oscillation cycles to simulate
n_points_per_cycle (int) – Number of time points per cycle

Returns:

Strain array gamma(t) stress: Stress array sigma(t)

Return type:

Example

>>> model = SGRGeneric()
>>> model.parameters.set_value("x", 1.5)
>>> strain, stress = model.simulate_laos(gamma_0=0.1, omega=1.0)

extract_laos_harmonics(stress, n_points_per_cycle=256)[source]¶

Extract Fourier harmonics from LAOS stress response.

Performs FFT analysis to extract harmonic amplitudes and phases:: sigma(t) = sum_n I_n * sin(n*omega*t + phi_n)

For LAOS, odd harmonics (n = 1, 3, 5, …) dominate due to symmetry.

Parameters:

stress (ndarray) – Stress time series from simulate_laos()
n_points_per_cycle (int) – Points per oscillation cycle

Returns:

I_1, I_3, I_5, …: Harmonic amplitudes
phi_1, phi_3, phi_5, …: Phase angles
I_3_I_1, I_5_I_1, …: Relative intensities

Return type:

Example

>>> strain, stress = model.simulate_laos(gamma_0=0.5, omega=1.0)
>>> harmonics = model.extract_laos_harmonics(stress)
>>> print(f"Third harmonic ratio: {harmonics['I_3_I_1']:.4f}")

compute_chebyshev_coefficients(strain, stress, gamma_0, omega, n_points_per_cycle=256)[source]¶

Compute Chebyshev decomposition of LAOS response.

Decomposes stress into elastic and viscous Chebyshev contributions:

sigma(gamma, gamma_dot) = sum_n e_n * T_n(gamma/gamma_0)

sum_n v_n * T_n(gamma_dot/gamma_dot_0)

where T_n are Chebyshev polynomials of the first kind.

Parameters:

strain (ndarray) – Strain array from simulate_laos()
stress (ndarray) – Stress array from simulate_laos()
gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency
n_points_per_cycle (int) – Points per oscillation cycle

Returns:

e_1, e_3, e_5: Elastic Chebyshev coefficients
v_1, v_3, v_5: Viscous Chebyshev coefficients
e_3_e_1, v_3_v_1: Normalized coefficients

Return type:

Example

>>> strain, stress = model.simulate_laos(gamma_0=0.5, omega=1.0)
>>> chebyshev = model.compute_chebyshev_coefficients(
...     strain, stress, gamma_0=0.5, omega=1.0
... )
>>> print(f"Strain stiffening ratio: {chebyshev['e_3_e_1']:.4f}")

get_lissajous_curve(gamma_0, omega, n_points=256, normalized=False)[source]¶

Generate Lissajous curve (stress vs strain) for LAOS.

Parameters:

gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)
n_points (int) – Number of points in curve
normalized (bool) – If True, normalize strain and stress

Returns:

Strain array (one period) stress: Stress array (one period)

Return type:

Example

>>> strain, stress = model.get_lissajous_curve(gamma_0=0.1, omega=1.0)
>>> plt.plot(strain, stress)  # Elastic Lissajous

evolve_x(t, gamma_dot, x0=None)[source]¶

Evolve noise temperature x(t) for aging/rejuvenation dynamics.

Evolution equation:

dx/dt = alpha_aging * (x_eq - x) + beta_rejuv * abs(gamma_dot) * (x_ss - x)

where x_ss = x_eq + x_ss_A * abs(gamma_dot)^x_ss_n is the steady-state value under shear.

Parameters:

t (ndarray) – Time array (s)
gamma_dot (ndarray) – Shear rate array (1/s), same shape as t
x0 (float | None) – Initial noise temperature. If None, uses current x value.

Returns:

Noise temperature evolution, same shape as t

Return type:

Example

>>> model = SGRGeneric(dynamic_x=True)
>>> t = np.linspace(0, 100, 1000)
>>> gamma_dot = np.where(t < 50, 10.0, 0.0)  # Shear then rest
>>> x_t = model.evolve_x(t, gamma_dot, x0=1.0)

free_energy_gradient(state)[source]¶

Compute gradient dF/dz of free energy.

The gradient components are: - dF/d(sigma): Conjugate to stress (strain-like) - dF/d(lambda): Conjugate to structure (chemical potential-like) - dF/d(x) = -S: Conjugate to temperature (dynamic x mode only)

Parameters:: state (ndarray) – State vector [sigma, lambda] or [sigma, lambda, x]
Return type:: ndarray
Returns:: Gradient [dF/d(sigma), dF/d(lambda)] or [dF/d(sigma), dF/d(lambda), dF/d(x)]

compute_entropy_production(state)[source]¶

Compute entropy production rate W at given state.

The entropy production is:: W = (dF/dz)^T * M(z) * (dF/dz) >= 0

This must be non-negative (second law of thermodynamics).

Parameters:: state (ndarray) – State vector [sigma, lambda] or [sigma, lambda, x]
Return type:: float
Returns:: Entropy production rate W (must be >= 0)
Raises:: Warning if W < 0 due to numerical errors –

verify_thermodynamic_consistency(state, tol=1e-10)[source]¶

Verify all GENERIC thermodynamic consistency conditions.

Checks: 1. Poisson bracket antisymmetry: L = -L^T 2. Friction matrix symmetry: M = M^T 3. Friction matrix positive semi-definiteness: eigenvalues >= 0 4. Entropy production non-negativity: W >= 0

Parameters:

state (ndarray) – State vector [sigma, lambda] or [sigma, lambda, x]
tol (float) – Numerical tolerance for consistency checks

Return type:

Returns:

Dictionary with consistency check results

STZ Models¶

STZ Conventional¶

rheojax.models.stz.conventional.STZConventional | Handbook: Shear Transformation Zone (STZ) Shear Transformation Zone model for amorphous plasticity (Langer 2008).

Parameters¶
Parameter	Units	Bounds	Description
`G0`	Pa	[1e6, 1e12]	Elastic shear modulus
`sigma_y`	Pa	[1e3, 1e9]	Yield stress scale
`chi_inf`		[0.01, 0.5]	Steady-state effective temperature
`tau0`	s	[1e-14, 1e-9]	Molecular attempt time
`epsilon0`		[0.01, 1.0]	Strain increment per flip
`c0`		[0.1, 100]	Specific heat
`ez`		[0.1, 5.0]	STZ formation energy

class rheojax.models.stz.conventional.STZConventional(variant='standard')[source]¶

Bases: STZBase

Conventional Shear Transformation Zone (STZ) Model.

Implements STZ plasticity with Langer (2008) formulation. Supports Minimal, Standard, and Full complexity variants.

Protocols: - Steady-State Flow: Algebraic solution for flow curve - Transient: Diffrax ODE integration for creep/relaxation/startup - SAOS/LAOS: Diffrax ODE integration + FFT for harmonic analysis

__init__(variant='standard')[source]¶

Initialize STZ Conventional Model.

Parameters:: variant (Literal['minimal', 'standard', 'full']) – Model variant (‘minimal’, ‘standard’, ‘full’)

simulate_laos(gamma_0, omega, n_cycles=2, n_points_per_cycle=256)[source]¶

Simulate LAOS response.

Parameters:

gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)
n_cycles (int) – Number of oscillation cycles
n_points_per_cycle (int) – Points per cycle

Return type:

Returns:

(strain, stress) arrays

extract_harmonics(stress, n_points_per_cycle=256)[source]¶

Extract Fourier harmonics from LAOS stress response.

Parameters:

stress (ndarray) – Stress array from simulate_laos
n_points_per_cycle (int) – Points per cycle

Return type:

Returns:

Dictionary with I_1, I_3, I_5 amplitudes and ratios

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

NumPyro/BayesianMixin model function.

Routes to appropriate prediction based on test_mode.

Elasto-Plastic Models (EPM)¶

Lattice EPM¶

rheojax.models.epm.lattice.LatticeEPM | Handbook: Elasto-Plastic Models (EPM) Mesoscopic lattice model for amorphous solids with spatial heterogeneity and avalanches.

Parameters¶
Parameter	Units	Bounds	Description
`mu`	Pa	[0.1, 100.0]	Shear modulus
`tau_pl`	s	[0.01, 100.0]	Plastic relaxation timescale
`sigma_c_mean`	Pa	[0.1, 10.0]	Mean yield threshold
`sigma_c_std`	Pa	[0.0, 1.0]	Disorder strength (std dev of thresholds)
`smoothing_width`	Pa	[0.01, 1.0]	Width for smooth yielding approx (inference)

class rheojax.models.epm.lattice.LatticeEPM(L=64, dt=0.01, mu=1.0, tau_pl=1.0, sigma_c_mean=1.0, sigma_c_std=0.1, n_fluid=1.0, n_bayesian_steps=200, fluidity_form='overstress')[source]¶

Bases: EPMBase

2D Lattice Elasto-Plastic Model (EPM).

A mesoscopic model for amorphous solids (glasses, gels) that explicitly resolves spatial heterogeneity, plastic avalanches, and stress redistribution.

Physics:

Lattice of elastoplastic blocks.
Elastic loading (affine).
Local yielding when stress > threshold.
Long-range stress redistribution via quadrupolar Eshelby propagator.
Structural renewal (disorder).

Parameters:

mu (float) – Shear modulus. Default 1.0.
tau_pl (float) – Plastic relaxation timescale. Default 1.0.
sigma_c_mean (float) – Mean yield threshold. Default 1.0.
sigma_c_std (float) – Disorder strength (std dev of thresholds). Default 0.1.
smoothing_width (float) – Width for smooth yielding approx (inference only). Default 0.1.

Configuration:: L (int): Lattice size (LxL). Default 64. dt (float): Time step. Default 0.01.

__init__(L=64, dt=0.01, mu=1.0, tau_pl=1.0, sigma_c_mean=1.0, sigma_c_std=0.1, n_fluid=1.0, n_bayesian_steps=200, fluidity_form='overstress')[source]¶

Initialize the Lattice EPM.

See EPMBase.__init__ for the description of fluidity_form.

Tensorial EPM¶

rheojax.models.epm.tensor.TensorialEPM | Handbook: Tensorial Elasto-Plastic Model (EPM) Full stress tensor implementation with normal stress difference predictions.

Parameters¶
Parameter	Units	Bounds	Description
`mu`	Pa	[0.1, 100.0]	Shear modulus
`nu`	dimensionless	[0.3, 0.5]	Poisson’s ratio (plane strain)
`tau_pl_shear`	s	[0.01, 100.0]	Plastic relaxation time for shear
`tau_pl_normal`	s	[0.01, 100.0]	Plastic relaxation time for normal stresses
`sigma_c_mean`	Pa	[0.1, 10.0]	Mean yield threshold
`sigma_c_std`	Pa	[0.0, 1.0]	Disorder strength (std dev)
`w_N1`	dimensionless	[0.1, 10.0]	Weight for N₁ in combined fitting
`hill_H`	dimensionless	[0.1, 5.0]	Hill anisotropy parameter H
`hill_N`	dimensionless	[0.1, 5.0]	Hill anisotropy parameter N

Configuration: L=64 (lattice size), dt=0.01 (timestep), yield_criterion="von_mises" or "hill"

Key Features: - Tracks full stress tensor [σ_xx, σ_yy, σ_xy] - Predicts normal stress differences N₁, N₂ - Von Mises (isotropic) or Hill (anisotropic) yield criteria - Flexible fitting: shear-only or combined [σ_xy, N₁]

class rheojax.models.epm.tensor.TensorialEPM(L=64, dt=0.01, mu=1.0, nu=0.48, tau_pl=1.0, tau_pl_shear=1.0, tau_pl_normal=1.0, sigma_c_mean=1.0, sigma_c_std=0.1, n_fluid=1.0, yield_criterion='von_mises', n_bayesian_steps=200, fluidity_form='overstress')[source]¶

Bases: EPMBase

3-Component Tensorial Lattice EPM.

A mesoscopic model for amorphous solids that explicitly tracks the full stress tensor, enabling predictions of normal stress differences and anisotropic flow.

Physics:

Lattice of elastoplastic blocks with tensorial stress state.
Elastic loading (affine) for all stress components.
Von Mises or Hill yield criterion for anisotropic materials.
Component-wise plastic flow rule (Prandtl-Reuss).
Tensorial Eshelby propagator for stress redistribution.

Parameters:

mu (float) – Shear modulus. Default 1.0.
nu (float) – Poisson’s ratio for plane strain. Default 0.48 (avoid 0.5 singularity).
tau_pl (float) – Base plastic relaxation timescale (legacy). Default 1.0.
tau_pl_shear (float) – Plastic relaxation time for shear. Default 1.0.
tau_pl_normal (float) – Plastic relaxation time for normal stresses. Default 1.0.
sigma_c_mean (float) – Mean yield threshold. Default 1.0.
sigma_c_std (float) – Disorder strength (std dev of thresholds). Default 0.1.
smoothing_width (float) – Width for smooth yielding approx (inference only). Default 0.1.
w_N1 (float) – Weight for N₁ in combined fitting loss. Default 1.0.
hill_H (float) – Hill anisotropy parameter H. Default 0.5.
hill_N (float) – Hill anisotropy parameter N. Default 1.5.

Configuration:: L (int): Lattice size (LxL). Default 64. dt (float): Time step. Default 0.01. yield_criterion (str): “von_mises” or “hill”. Default “von_mises”.

__init__(L=64, dt=0.01, mu=1.0, nu=0.48, tau_pl=1.0, tau_pl_shear=1.0, tau_pl_normal=1.0, sigma_c_mean=1.0, sigma_c_std=0.1, n_fluid=1.0, yield_criterion='von_mises', n_bayesian_steps=200, fluidity_form='overstress')[source]¶

Initialize the Tensorial EPM.

Parameters:

L (int) – Lattice size (LxL grid).
dt (float) – Time step for integration.
mu (float) – Shear modulus.
nu (float) – Poisson’s ratio for plane strain constraint.
tau_pl (float) – Base plastic relaxation time (for compatibility).
tau_pl_shear (float) – Plastic relaxation time for shear components.
tau_pl_normal (float) – Plastic relaxation time for normal stress components.
sigma_c_mean (float) – Mean yield threshold.
sigma_c_std (float) – Standard deviation of yield thresholds (disorder).
yield_criterion (str) – Yield criterion name (“von_mises” or “hill”).
n_bayesian_steps (int) – Number of time steps for Bayesian inference. Default 200.

get_shear_stress(stress)[source]¶

Extract shear stress component from stress tensor.

Parameters:: stress (Array) – Stress tensor of shape (3, L, L) or (…, 3, L, L).
Return type:: Array
Returns:: Shear stress σ_xy, shape (L, L) or (…, L, L).

get_normal_stress_differences(stress, nu=None)[source]¶

Compute normal stress differences from stress tensor.

For plane strain: σ_zz = ν(σ_xx + σ_yy)

Normal stress differences: - N₁ = σ_xx - σ_yy - N₂ = σ_yy - σ_zz

Parameters:

stress (Array) – Stress tensor of shape (3, L, L).
nu (float | None) – Poisson’s ratio. If None, uses parameter value.

Return type:

tuple[Array, Array]

Returns:

Tuple (N₁, N₂), each of shape (L, L).

predict_normal_stresses(data, **kwargs)[source]¶

Convenience method to predict normal stress differences.

Runs the simulation and extracts N₁ and N₂ spatial averages over time.

Parameters:

data (RheoData) – RheoData with protocol specification.
**kwargs – Additional arguments passed to predict().

Return type:

tuple[Array, Array]

Returns:

Tuple (N₁_array, N₂_array) with time-averaged values.

Raises:

NotImplementedError – Not yet implemented (future feature).

SPP LAOS Models¶

SPP Yield Stress¶

rheojax.models.spp.spp_yield_stress.SPPYieldStress | Handbook: SPP Yield Stress Model Yield stress model for Sequence of Physical Processes (SPP) LAOS analysis.

Parameters¶
Parameter	Units	Bounds	Description
`sigma_y`	Pa	[0, 1e+06]	Static yield stress
`G_cage`	Pa	[0.001, 1e+09]	Cage modulus
`n`	dimensionless	[0.01, 2]	Power-law exponent

Note: Use with rheojax.transforms.spp_decomposer.SPPDecomposer for complete LAOS analysis. See SPP Analysis API for the full SPP API reference.

class rheojax.models.spp.spp_yield_stress.SPPYieldStress[source]¶

Bases: BaseModel

SPP-based yield stress model for LAOS analysis.

This model extracts yield stress parameters from amplitude-sweep LAOS data using the SPP framework. It parameterizes the nonlinear response in terms of physically meaningful quantities: cage modulus, static/dynamic yield stresses, and post-yield flow parameters.

The model supports both OSCILLATION mode (LAOS amplitude sweeps) and ROTATION mode (steady shear flow curves) for comprehensive yield stress characterization.

Parameters:

G_cage (float) – Apparent cage modulus (Pa), elastic stiffness of the cage structure
sigma_sy_scale (float) – Static yield stress scale factor (Pa)
sigma_sy_exp (float) – Static yield stress exponent for amplitude scaling
sigma_dy_scale (float) – Dynamic yield stress scale factor (Pa)
sigma_dy_exp (float) – Dynamic yield stress exponent for amplitude scaling
eta_inf (float) – Infinite-shear viscosity (Pa·s)
n_power_law (float) – Power-law flow index (dimensionless)
noise (float) – Observation noise scale (Pa), for Bayesian inference
Equations (Constitutive)
----------------------
γ_0) (For OSCILLATION (LAOS at amplitude) – σ_sy(γ_0) = sigma_sy_scale * γ_0^sigma_sy_exp σ_dy(γ_0) = sigma_dy_scale * γ_0^sigma_dy_exp σ_max(γ_0) = G_cage * γ_0 + η_inf * ω * γ_0
γ̇) (For ROTATION (steady shear at rate) – σ(γ̇) = σ_dy + η_inf * γ̇^n_power_law (Herschel-Bulkley-like)
Modes (Test)
----------
OSCILLATION (-)
ROTATION (-)

Examples

Fit to amplitude sweep data:

>>> from rheojax.models import SPPYieldStress
>>> from rheojax.transforms import SPPDecomposer
>>>
>>> # Amplitude sweep data (multiple γ_0 values)
>>> gamma_0_values = jnp.array([0.01, 0.1, 0.5, 1.0, 2.0, 5.0])
>>> sigma_sy_data = jnp.array([1.0, 10.0, 40.0, 80.0, 140.0, 280.0])
>>>
>>> model = SPPYieldStress()
>>> model.fit(gamma_0_values, sigma_sy_data, test_mode='oscillation')
>>> print(model)

Bayesian inference:

>>> result = model.fit_bayesian(gamma_0_values, sigma_sy_data,
...                             test_mode='oscillation')
>>> print(result.summary)

__init__()[source]¶: Initialize SPP yield stress model.

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

Model function for predictions and Bayesian inference.

This is the NumPyro-facing model used by BayesianMixin. It maps inputs and parameter vectors to stress predictions.

Parameters:

X (Array) – Independent variable - OSCILLATION: strain amplitudes (gamma_0) - ROTATION: shear rates (gamma_dot)
params (Array) –

[G_cage, sigma_sy_scale, sigma_sy_exp, sigma_dy_scale,
sigma_dy_exp, eta_inf, n_power_law, noise]
test_mode (TestModeEnum | None) – Mode selector; defaults to the model’s stored test_mode.

Returns:

Predicted stress values.

Return type:

bayesian_prior_factory(name, lower, upper)[source]¶

Create NumPyro-friendly priors for model parameters.

Uses physically-motivated prior distributions: - LogNormal for scale parameters (positive, often log-distributed) - Beta for bounded exponents - HalfCauchy for noise scale

Parameters:

name (str) – Parameter name
lower (float | None) – Lower bound
upper (float | None) – Upper bound

Returns:

NumPyro distribution, or None to use default uniform

Return type:

Distribution | None

predict_amplitude_sweep(gamma_0_array, omega=1.0, yield_type='static')[source]¶

Predict full amplitude sweep response.

Computes both static and dynamic yield stresses across amplitudes.

Parameters:

gamma_0_array (ndarray) – Strain amplitudes to evaluate
omega (float) – Angular frequency (rad/s)
yield_type (str) – ‘static’, ‘dynamic’, or ‘both’

Returns:

Dictionary with: - gamma_0: strain amplitudes - sigma_sy: static yield stresses (if requested) - sigma_dy: dynamic yield stresses (if requested)

Return type:

predict_flow_curve(gamma_dot_array)[source]¶

Predict steady shear flow curve.

Parameters:: gamma_dot_array (ndarray) – Shear rates to evaluate
Returns:: Dictionary with: - gamma_dot: shear rates - sigma: shear stress - eta_app: apparent viscosity
Return type:: dict

predict_rheo(rheo_data, test_mode=None)[source]¶

Predict rheological response for RheoData.

Parameters:

rheo_data (RheoData) – Input rheological data
test_mode (TestModeEnum | None) – Test mode override

Returns:

Predicted response

Return type:

__repr__()[source]¶

String representation.

Return type:: str

Fluidity Models¶

FluidityLocal¶

rheojax.models.fluidity.local.FluidityLocal | Handbook: Fluidity Local (Homogeneous Fluidity Model) — Handbook Local (0D) fluidity model with aging/rejuvenation kinetics for thixotropic yield-stress fluids.

Parameters¶
Parameter	Units	Bounds	Description
`G`	Pa	[1e+03, 1e+09]	Elastic modulus
`tau_y`	Pa	[10, 1e+06]	Yield stress
`K`	Pa·sⁿ	[1, 1e+06]	Flow consistency (HB K parameter)
`n_flow`	–	[0.1, 2]	Flow exponent (HB n parameter)
`f_eq`	1/(Pa·s)	[1e-12, 0.001]	Equilibrium fluidity (aging limit)
`f_inf`	1/(Pa·s)	[1e-06, 1]	High-shear fluidity (rejuvenation limit)
`theta`	s	[0.1, 1e+04]	Structural relaxation time (aging timescale)
`a`	–	[0, 100]	Rejuvenation amplitude
`n_rejuv`	–	[0, 2]	Rejuvenation exponent

class rheojax.models.fluidity.local.FluidityLocal[source]¶

Bases: FluidityBase

Local (0D) Fluidity Model for yield-stress fluids.

Implements a homogeneous fluidity model where the material state is characterized by a single fluidity value f(t) that evolves via:

df/dt = (f_eq - f)/θ + a|γ̇|^n(f_inf - f)

This captures: - Aging: structural build-up at rest, f → f_eq - Rejuvenation: flow-induced breakdown, f → f_inf

The stress evolves as a viscoelastic solid with plastic flow: dσ/dt = G(γ̇ - σf)

Protocols: - Flow Curve: Algebraic steady-state solution - Startup: ODE integration with constant γ̇ - Relaxation: ODE integration with γ̇=0, stress decays - Creep: ODE integration with constant σ - SAOS/LAOS: ODE integration + FFT for moduli

__init__()[source]¶: Initialize Local Fluidity Model.

simulate_laos(gamma_0, omega, n_cycles=2, n_points_per_cycle=256)[source]¶

Simulate LAOS response.

Parameters:

gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)
n_cycles (int) – Number of oscillation cycles
n_points_per_cycle (int) – Points per cycle

Return type:

Returns:

(strain, stress) arrays

extract_harmonics(stress, n_points_per_cycle=256)[source]¶

Extract Fourier harmonics from LAOS stress response.

Parameters:

stress (ndarray) – Stress array from simulate_laos
n_points_per_cycle (int) – Points per cycle

Return type:

Returns:

Dictionary with I_1, I_3, I_5 amplitudes and ratios

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

NumPyro/BayesianMixin model function.

Routes to appropriate prediction based on test_mode. Accepts protocol-specific kwargs (gamma_dot, sigma_applied, etc.).

FluidityNonlocal¶

rheojax.models.fluidity.nonlocal_model.FluidityNonlocal | Handbook: Fluidity Nonlocal (Coussot-Ovarlez Cooperative Model) — Handbook Nonlocal (1D) fluidity model with cooperativity length for shear banding prediction.

Parameters¶
Parameter	Units	Bounds	Description
`G`	Pa	[1e+03, 1e+09]	Elastic modulus
`tau_y`	Pa	[10, 1e+06]	Yield stress
`K`	Pa·sⁿ	[1, 1e+06]	Flow consistency (HB K parameter)
`n_flow`	–	[0.1, 2]	Flow exponent (HB n parameter)
`f_eq`	1/(Pa·s)	[1e-12, 0.001]	Equilibrium fluidity (aging limit)
`f_inf`	1/(Pa·s)	[1e-06, 1]	High-shear fluidity (rejuvenation limit)
`theta`	s	[0.1, 1e+04]	Structural relaxation time
`a`	–	[0, 100]	Rejuvenation amplitude
`n_rejuv`	–	[0, 2]	Rejuvenation exponent
`xi`	m	[1e-09, 0.001]	Cooperativity length (non-local diffusion scale)

class rheojax.models.fluidity.nonlocal_model.FluidityNonlocal(N_y=64, gap_width=0.001)[source]¶

Bases: FluidityBase

Non-Local (1D PDE) Fluidity Model for yield-stress fluids.

Implements the Coussot-Ovarlez non-local fluidity model where the fluidity field f(y,t) evolves across the gap (y-direction) via:

∂f/∂t = (f_loc(σ) - f)/θ + ξ²∂²f/∂y²

where: - f_loc(σ) is the local equilibrium fluidity from HB flow curve - θ is the relaxation time - ξ is the cooperativity length (non-local diffusion)

This captures shear banding: localized flow in yield-stress fluids where the cooperativity length ξ determines band width.

Key features: - 1D Couette gap discretization (N_y points) - Neumann (zero-flux) boundary conditions at walls - Diffrax Dopri5 solver (explicit, robust) for PDE - Shear banding metrics: CV and max/min ratio

N_y¶: Number of grid points across gap

gap_width¶: Physical gap width (m)

__init__(N_y=64, gap_width=0.001)[source]¶

Initialize Non-Local Fluidity Model.

Parameters:

N_y (int) – Number of spatial grid points (default 64)
gap_width (float) – Physical gap width in meters (default 1 mm)

get_fluidity_profile(time_idx=-1)[source]¶

Get fluidity profile at specified time index.

Parameters:: time_idx (int) – Time index (default -1 for final time)
Return type:: ndarray
Returns:: Fluidity field across gap, shape (N_y,)

get_shear_banding_metric(f_field=None)[source]¶

Compute coefficient of variation as shear banding metric.

CV = std(f) / mean(f) CV > 0.3 typically indicates significant shear banding.

Parameters:: f_field (ndarray | None) – Fluidity field. If None, uses final simulation state.
Return type:: float
Returns:: Coefficient of variation (dimensionless)

get_banding_ratio(f_field=None)[source]¶

Compute max/min fluidity ratio as banding metric.

ratio > 10 indicates strong localization.

Parameters:: f_field (ndarray | None) – Fluidity field. If None, uses final simulation state.
Return type:: float
Returns:: Banding ratio (dimensionless)

is_banding(f_field=None, cv_threshold=0.3)[source]¶

Check if shear banding is occurring.

Parameters:

f_field (ndarray | None) – Fluidity field. If None, uses final simulation state.
cv_threshold (float) – CV threshold for banding (default 0.3)

Return type:

bool

Returns:

True if CV > threshold

simulate_laos(gamma_0, omega, n_cycles=2, n_points_per_cycle=256, f_init=None)[source]¶

Simulate LAOS response.

Parameters:

gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)
n_cycles (int) – Number of oscillation cycles
n_points_per_cycle (int) – Points per cycle

Return type:

Returns:

(strain, stress) arrays

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

NumPyro/BayesianMixin model function.

Accepts protocol-specific kwargs (gamma_dot, sigma_applied, etc.).

Fluidity-Saramito EVP Models¶

FluiditySaramitoLocal¶

rheojax.models.fluidity.saramito.local.FluiditySaramitoLocal | Handbook: Fluidity-Saramito EVP Model Local (0D) elastoviscoplastic model combining Saramito tensorial viscoelasticity with thixotropic fluidity evolution and Von Mises yielding.

Parameters¶
Parameter	Units	Bounds	Description
`G`	Pa	[10, 1e+08]	Elastic modulus
`eta_s`	Pa·s	[0, 1e+03]	Solvent viscosity
`tau_y0`	Pa	[0.1, 1e+05]	Base yield stress (Von Mises threshold)
`K_HB`	Pa·sⁿ	[0.01, 1e+05]	Herschel-Bulkley consistency index
`n_HB`	–	[0.1, 1.5]	Herschel-Bulkley flow exponent
`f_age`	1/(Pa·s)	[1e-12, 0.01]	Aging fluidity limit
`f_flow`	1/(Pa·s)	[1e-06, 1]	Flow fluidity limit (rejuvenation)
`t_a`	s	[0.01, 1e+05]	Aging timescale (structural build-up)
`b`	–	[0, 1e+03]	Rejuvenation amplitude
`n_rej`	–	[0.1, 3]	Rejuvenation rate exponent

class rheojax.models.fluidity.saramito.local.FluiditySaramitoLocal(coupling='minimal')[source]¶

Bases: FluiditySaramitoBase

Local (0D) Fluidity-Saramito Model for elastoviscoplastic fluids.

Implements a homogeneous Saramito EVP model where the material state is characterized by tensorial stress [τ_xx, τ_yy, τ_xy] and scalar fluidity f(t) that evolve via coupled ODEs.

The constitutive equation (Saramito 2009): λ(∇̂τ) + α(τ)τ = 2η_p D

where:

λ = 1/f: Fluidity-dependent relaxation time
∇̂τ: Upper-convected derivative
α = max(0, 1 - τ_y/|τ|): Von Mises plasticity
η_p = G/f: Polymeric viscosity

Fluidity evolves via: df/dt = (f_age - f)/t_a + b|γ̇|^n(f_flow - f)

Parameters:: coupling (Literal['minimal', 'full']) – Coupling mode for yield stress evolution: - “minimal”: τ_y = τ_y0 (constant) - “full”: τ_y(f) = τ_y0 + a_y/f^m (aging yield stress)

Examples

Basic flow curve fitting:

>>> model = FluiditySaramitoLocal(coupling="minimal")
>>> model.fit(gamma_dot, sigma, test_mode="flow_curve")
>>> sigma_pred = model.predict(gamma_dot)

Startup with stress overshoot:

>>> model = FluiditySaramitoLocal()
>>> model.fit(t, sigma, test_mode="startup", gamma_dot=1.0)
>>> strain, stress, f = model.simulate_startup(t, gamma_dot=1.0)

Bayesian inference:

>>> result = model.fit_bayesian(gamma_dot, sigma, test_mode="flow_curve",
...                             num_warmup=1000, num_samples=2000)
>>> intervals = model.get_credible_intervals(result.posterior_samples)

__init__(coupling='minimal')[source]¶

Initialize Local Fluidity-Saramito Model.

Parameters:: coupling (Literal['minimal', 'full']) – Coupling mode for yield stress evolution

simulate_startup(t, gamma_dot, t_wait=0.0)[source]¶

Simulate startup response with full trajectory.

Parameters:

t (ndarray) – Time array (s)
gamma_dot (float) – Applied shear rate (1/s)
t_wait (float) – Waiting time before startup (s)

Return type:

Returns:

strain (np.ndarray) – Accumulated strain γ(t)
stress (np.ndarray) – Shear stress τ_xy(t) (Pa)
fluidity (np.ndarray) – Fluidity f(t) (1/(Pa·s))

simulate_creep(t, sigma_applied, t_wait=0.0)[source]¶

Simulate creep response.

Parameters:

t (ndarray) – Time array (s)
sigma_applied (float) – Applied stress (Pa)
t_wait (float) – Waiting time before creep (s)

Return type:

Returns:

strain (np.ndarray) – Accumulated strain γ(t)
fluidity (np.ndarray) – Fluidity f(t) (1/(Pa·s))

simulate_relaxation(t, gamma_0=None, sigma_0=None, t_wait=0.0)[source]¶

Simulate stress relaxation after a step strain.

After an instantaneous step strain γ_0, the stress jumps elastically to σ_0 = G·γ_0 and then decays under zero applied rate. This is the analogue of simulate_creep() / simulate_startup() for the relaxation protocol — without it the only public path is predict(), which silently falls back to sigma_init = tau_y0 (the von Mises plasticity parameter is then exactly zero and stress never decays).

Parameters:

t (ndarray) – Time array (s).
gamma_0 (float | None) – Step strain amplitude. The initial stress is taken as G * gamma_0. Mutually exclusive with sigma_0.
sigma_0 (float | None) – Initial shear stress (Pa) immediately after the step. Takes precedence over gamma_0 if both are given.
t_wait (float) – Pre-relaxation aging time (s). Sets the initial fluidity via the same exponential schedule used by the other simulate_* helpers.

Return type:

Returns:

stress (np.ndarray) – Shear stress τ_xy(t) (Pa).
fluidity (np.ndarray) – Fluidity f(t) (1/(Pa·s)).

simulate_laos(gamma_0, omega, n_cycles=3, n_points_per_cycle=256)[source]¶

Simulate LAOS response.

Parameters:

gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)
n_cycles (int) – Number of oscillation cycles
n_points_per_cycle (int) – Points per cycle

Return type:

Returns:

t (np.ndarray) – Time array (s)
strain (np.ndarray) – Strain array
stress (np.ndarray) – Stress array (Pa)

extract_harmonics(stress, n_points_per_cycle=256)[source]¶

Extract Fourier harmonics from LAOS stress response.

Parameters:

stress (ndarray) – Stress array from simulate_laos
n_points_per_cycle (int) – Points per cycle

Returns:

Dictionary with I_1, I_3, I_5 amplitudes and ratios

Return type:

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

NumPyro/BayesianMixin model function.

Routes to appropriate prediction based on test_mode. Accepts protocol-specific kwargs (gamma_dot, sigma_applied, etc.).

Parameters:

X (array-like) – Independent variable
params (array-like) – Parameter values in order
test_mode (str, optional) – Override stored test mode
**kwargs – Protocol-specific arguments (gamma_dot, sigma_applied, gamma_0, omega)

Returns:

Predicted response

Return type:

jnp.ndarray

get_overshoot_ratio(gamma_dot, t_max=100.0, n_points=1000)[source]¶

Compute stress overshoot ratio σ_max / σ_ss.

Parameters:

gamma_dot (float) – Shear rate (1/s)
t_max (float) – Maximum simulation time (s)
n_points (int) – Number of time points

Returns:

Overshoot ratio (1.0 = no overshoot)

Return type:

float

get_critical_stress()[source]¶

Get critical stress for creep bifurcation.

Returns:: Critical stress σ_c ≈ τ_y (Pa)
Return type:: float

FluiditySaramitoNonlocal¶

rheojax.models.fluidity.saramito.nonlocal_model.FluiditySaramitoNonlocal | Handbook: Fluidity-Saramito EVP Model Nonlocal (1D) Saramito EVP model with spatial cooperativity for shear banding.

Parameters¶
Parameter	Units	Bounds	Description
`G`	Pa	[10, 1e+08]	Elastic modulus
`eta_s`	Pa·s	[0, 1e+03]	Solvent viscosity
`tau_y0`	Pa	[0.1, 1e+05]	Base yield stress (Von Mises threshold)
`K_HB`	Pa·sⁿ	[0.01, 1e+05]	Herschel-Bulkley consistency index
`n_HB`	–	[0.1, 1.5]	Herschel-Bulkley flow exponent
`f_age`	1/(Pa·s)	[1e-12, 0.01]	Aging fluidity limit
`f_flow`	1/(Pa·s)	[1e-06, 1]	Flow fluidity limit
`t_a`	s	[0.01, 1e+05]	Aging timescale
`b`	–	[0, 1e+03]	Rejuvenation amplitude
`n_rej`	–	[0.1, 3]	Rejuvenation rate exponent
`xi`	m	[1e-07, 0.01]	Cooperativity length (interface width)

class rheojax.models.fluidity.saramito.nonlocal_model.FluiditySaramitoNonlocal(coupling='minimal', N_y=51, H=0.001, xi=1e-05)[source]¶

Bases: FluiditySaramitoBase

Nonlocal (1D) Fluidity-Saramito Model with spatial diffusion.

Implements a spatially-resolved Saramito EVP model where fluidity varies across a Couette gap and can form shear bands.

The fluidity evolution includes a diffusion term: ∂f/∂t = (f_loc - f)/t_a + D_f * ∇²f

where D_f = ξ²/t_a is the fluidity diffusivity and ξ is the cooperativity length (interface width parameter).

Parameters:

coupling (Literal['minimal', 'full']) – Coupling mode for yield stress evolution
N_y (int) – Number of spatial grid points
H (float) – Gap width (m)
xi (float) – Cooperativity length (m)

Notes

The model solves a coupled PDE system for [Σ, f(y)] where Σ is the bulk (gap-averaged) stress. In Couette geometry, stress is approximately uniform across the gap, enabling this simplification.

Shear bands appear when the fluidity profile develops large gradients, with high-fluidity (flowing) bands coexisting with low-fluidity (jammed) regions.

Examples

Basic flow curve with banding check:

>>> model = FluiditySaramitoNonlocal(N_y=51, H=1e-3, xi=1e-5)
>>> model.fit(gamma_dot, sigma, test_mode="flow_curve")
>>> is_banded, cv, ratio = model.detect_shear_bands()
>>> print(f"Shear banding detected: {is_banded}")

Startup transient with spatial profile:

>>> model = FluiditySaramitoNonlocal()
>>> t, sigma, f_field = model.simulate_startup(t, gamma_dot=1.0)
>>> model.plot_fluidity_profile()  # Shows spatial variation

__init__(coupling='minimal', N_y=51, H=0.001, xi=1e-05)[source]¶

Initialize Nonlocal Fluidity-Saramito Model.

Parameters:

coupling (Literal['minimal', 'full']) – Coupling mode for yield stress evolution
N_y (int) – Number of spatial grid points
H (float) – Gap width (m)
xi (float) – Cooperativity length (m)

simulate_startup(t, gamma_dot)[source]¶

Simulate startup with spatial resolution.

Parameters:

t (ndarray) – Time array (s)
gamma_dot (float) – Applied shear rate (1/s)

Return type:

Returns:

t (np.ndarray) – Time array
sigma (np.ndarray) – Bulk stress (Pa)
f_field (np.ndarray) – Final fluidity profile, shape (N_y,)

simulate_creep(t, sigma_applied)[source]¶

Simulate creep with spatial resolution.

Parameters:

t (ndarray) – Time array (s)
sigma_applied (float) – Applied stress (Pa)

Return type:

Returns:

gamma (np.ndarray) – Bulk strain
f_field (np.ndarray) – Final fluidity profile

detect_shear_bands(f_profile=None, threshold=0.3)[source]¶

Detect shear banding from fluidity profile.

Parameters:

f_profile (ndarray | None) – Fluidity field. If None, uses stored profile.
threshold (float) – CV threshold for banding detection

Return type:

tuple[bool, float, float]

Returns:

is_banded (bool) – True if shear banding detected
cv (float) – Coefficient of variation
ratio (float) – Max/min fluidity ratio

get_banding_metrics(f_profile=None)[source]¶

Get detailed shear banding metrics.

Parameters:: f_profile (ndarray | None) – Fluidity field. If None, uses stored profile.
Returns:: Metrics including cv, ratio, band_width, etc.
Return type:: dict[str, float]

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

NumPyro/BayesianMixin model function.

Accepts protocol-specific kwargs (gamma_dot, sigma_applied, etc.).

Parameters:

X (array-like) – Independent variable
params (array-like) – Parameter values
test_mode (str, optional) – Override stored test mode
**kwargs – Protocol-specific arguments (gamma_dot, sigma_applied)

Returns:

Predicted response

Return type:

jnp.ndarray

property y_grid: ndarray¶

Get spatial grid across gap.

Returns:: Position array (m)
Return type:: np.ndarray

__repr__()[source]¶

Return string representation.

Return type:: str

DMT Thixotropic Models¶

DMTLocal¶

rheojax.models.dmt.local.DMTLocal | Handbook: DMT Thixotropic Models Local (0D) de Souza Mendes-Thompson thixotropic model with scalar structure parameter and exponential or Herschel-Bulkley viscosity closure.

Parameters¶
Parameter	Units	Bounds	Description
`eta_0`	Pa·s	[100, 1e+08]	Zero-shear viscosity (fully structured, \(\lambda=1\))
`eta_inf`	Pa·s	[0.001, 100]	Infinite-shear viscosity (fully broken, \(\lambda=0\))
`t_eq`	s	[0.1, 1e+04]	Structural equilibrium (buildup) timescale
`a`	–	[0.001, 100]	Breakdown rate coefficient
`c`	–	[0.1, 2]	Breakdown rate exponent (shear rate sensitivity)

class rheojax.models.dmt.local.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)

See also

DMTNonlocal: Nonlocal (1D) variant with shear banding
FluidityLocal: Simpler fluidity-based thixotropic model

References

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.” Rheol. Acta 52, 673-694.

__init__(closure='exponential', include_elasticity=True)[source]¶

Initialize DMTLocal model.

simulate_startup(gamma_dot, t_end, dt=0.01, lam_init=1.0)[source]¶

Simulate startup of steady shear from rest.

Parameters:

gamma_dot (float) – Applied constant shear rate [1/s]
t_end (float) – Simulation end time [s]
dt (float) – Time step [s]
lam_init (float) – Initial structure parameter (default: 1.0, fully structured)

Return type:

Returns:

t (array) – Time array [s]
stress (array) – Stress response [Pa]
lam (array) – Structure parameter evolution

simulate_relaxation(t_end, dt=0.01, sigma_init=100.0, lam_init=0.5)[source]¶

Simulate stress relaxation after cessation of shear.

Requires include_elasticity=True.

Parameters:

t_end (float) – Simulation end time [s]
dt (float) – Time step [s]
sigma_init (float) – Initial stress at cessation [Pa]
lam_init (float) – Initial structure at cessation

Return type:

Returns:

t (array) – Time array [s]
stress (array) – Relaxing stress [Pa]
lam (array) – Recovering structure

simulate_creep(sigma_0, t_end, dt=0.01, lam_init=1.0)[source]¶

Simulate creep under constant applied stress.

For the Maxwell variant (include_elasticity=True), the total strain includes both elastic and viscous contributions:

γ(t) = γ_e(t) + γ_v(t)

where: - γ_e(t) = σ₀/G(λ(t)) is the elastic strain (changes with structure) - γ_v(t) = ∫₀ᵗ σ₀/η(λ(s)) ds is the viscous strain

This correctly captures: - Initial elastic jump: γ(0⁺) = σ₀/G(λ_init) - Elastic strain recovery/growth as structure evolves - Viscous flow accumulation

Parameters:

sigma_0 (float) – Applied constant stress [Pa]
t_end (float) – Simulation end time [s]
dt (float) – Time step [s]
lam_init (float) – Initial structure parameter

Return type:

tuple[ndarray, ndarray, ndarray, ndarray]

Returns:

t (array) – Time array [s]
gamma (array) – Total accumulated strain (elastic + viscous for Maxwell variant)
gamma_dot (array) – Total shear rate evolution [1/s]
lam (array) – Structure parameter evolution

predict_saos(omega, lam_0=1.0)[source]¶

Predict SAOS moduli G’(ω) and G’’(ω).

Requires include_elasticity=True. Assumes small amplitude so structure remains at λ₀.

Parameters:

omega (ndarray) – Angular frequency [rad/s]
lam_0 (float) – Reference structure level (default: 1.0, fully structured)

Return type:

Returns:

G_prime (array) – Storage modulus G’(ω) [Pa]
G_double_prime (array) – Loss modulus G’’(ω) [Pa]

simulate_laos(gamma_0, omega, n_cycles=10, points_per_cycle=128, lam_init=1.0)[source]¶

Simulate LAOS (Large Amplitude Oscillatory Shear).

Parameters:

gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency [rad/s]
n_cycles (int) – Number of cycles to simulate
points_per_cycle (int) – Points per cycle for discretization
lam_init (float) – Initial structure parameter

Returns:

‘t’: time array ‘strain’: strain γ(t) ‘strain_rate’: strain rate γ̇(t) ‘stress’: stress σ(t) ‘lam’: structure λ(t)

Return type:

extract_harmonics(laos_result, n_harmonics=5)[source]¶

Extract Fourier harmonics from LAOS stress response.

Parameters:

laos_result (dict) – Result from simulate_laos()
n_harmonics (int) – Number of harmonics to extract

Returns:

‘sigma_prime’: in-phase coefficients (odd harmonics) ‘sigma_double_prime’: out-of-phase coefficients ‘I_n_1’: normalized harmonic intensities

Return type:

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

NumPyro/BayesianMixin model function for DMT.

Routes to appropriate JAX-traceable prediction based on test_mode. Required by BayesianMixin for NumPyro NUTS sampling.

Parameters:

X (array-like) – Independent variable (shear rate, time, or frequency)
params (array-like) – Parameter values in ParameterSet order
test_mode (str, optional) – Override stored test mode

Returns:

Predicted response (stress, strain, or complex modulus)

Return type:

jnp.ndarray

DMTNonlocal¶

rheojax.models.dmt.nonlocal_model.DMTNonlocal | Handbook: DMT Thixotropic Models Nonlocal (1D) DMT model with structure diffusion for shear banding prediction.

Parameters¶
Parameter	Units	Bounds	Description
`eta_0`	Pa·s	[100, 1e+08]	Zero-shear viscosity (fully structured)
`eta_inf`	Pa·s	[0.001, 100]	Infinite-shear viscosity (fully broken)
`t_eq`	s	[0.1, 1e+04]	Structural equilibrium timescale
`a`	–	[0.001, 100]	Breakdown rate coefficient
`c`	–	[0.1, 2]	Breakdown rate exponent
`D_lambda`	m²/s	[1e-10, 0.01]	Structure diffusion coefficient (cooperativity)

class rheojax.models.dmt.nonlocal_model.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)

See also

DMTLocal: Local (0D) variant for homogeneous flow
FluidityNonlocal: Simpler nonlocal fluidity model

References

Coussot, P. et al. (2002). “Viscosity bifurcation in thixotropic,: yielding fluids.” J. Rheol. 46, 573-589.

__init__(closure='exponential', include_elasticity=True, n_points=51, gap_width=0.001)[source]¶

Initialize DMTNonlocal model.

get_cooperativity_length()[source]¶

Compute cooperativity length scale.

ξ = √(D_λ · t_eq)

Returns:: Cooperativity length [m]
Return type:: float

simulate_steady_shear(gamma_dot_avg, t_end, dt=0.1, lam_init=1.0)[source]¶

Simulate approach to steady state under imposed average shear rate.

For planar Couette: v(0) = 0, v(H) = V_wall = γ̇_avg · H

The stress is uniform (σ(y) = Σ) at low Reynolds number. The local shear rate γ̇(y) varies to satisfy the local constitutive relation with the uniform stress.

Parameters:

gamma_dot_avg (float) – Average (imposed) shear rate [1/s]
t_end (float) – Simulation end time [s]
dt (float) – Time step [s]
lam_init (float | ndarray) – Initial structure (scalar for uniform, array for profile)

Returns:

‘t’: time array ‘lam’: structure profiles λ(y, t_i) ‘gamma_dot’: shear rate profiles γ̇(y, t_i) ‘velocity’: velocity profiles v(y, t_i) ‘stress’: wall stress Σ(t)

Return type:

detect_banding(result, threshold=0.1)[source]¶

Detect shear banding from steady-state profiles.

A shear band is detected when the shear rate profile shows significant spatial variation (standard deviation / mean > threshold).

Parameters:

result (dict) – Result from simulate_steady_shear()
threshold (float) – Relative variation threshold for banding detection

Returns:

‘is_banding’: bool ‘band_contrast’: max/min shear rate ratio ‘band_width’: approximate band width [m] ‘band_location’: center of high-shear band [m] ‘gamma_dot_profile’: final shear rate profile ‘lam_profile’: final structure profile

Return type:

plot_profiles(result, time_indices=None, figsize=(12, 4))[source]¶

Plot structure and shear rate profiles.

Parameters:

result (dict) – Result from simulate_steady_shear()
time_indices (list[int] | None) – Indices of time points to plot (default: [0, -1])
figsize (tuple) – Figure size

Returns:

Matplotlib figure and axes

Return type:

fig, axes

IKH Models (Isotropic-Kinematic Hardening)¶

MIKH¶

rheojax.models.ikh.mikh.MIKH | Handbook: Maxwell-Isotropic-Kinematic Hardening (MIKH) Single-mode thixotropic elasto-viscoplastic model with Maxwell viscoelasticity, Armstrong-Frederick kinematic hardening, and thixotropic structure evolution.

Parameters¶
Parameter	Units	Bounds	Description
`G`	Pa	[0.1, 1e+09]	Shear modulus
`eta`	Pa·s	[0.001, 1e+12]	Maxwell viscosity (relaxation time \(= \eta/G\))
`C`	Pa	[0, 1e+09]	Kinematic hardening modulus
`gamma_dyn`	–	[0, 1e+04]	Dynamic recovery parameter
`m`	–	[0.5, 3]	AF recovery exponent
`sigma_y0`	Pa	[0, 1e+09]	Minimal yield stress (destructured)
`delta_sigma_y`	Pa	[0, 1e+09]	Structural yield stress contribution
`tau_thix`	s	[1e-06, 1e+12]	Thixotropic rebuilding timescale
`Gamma`	–	[0, 1e+04]	Structural breakdown coefficient
`eta_inf`	Pa·s	[0, 1e+09]	High-shear viscosity
`mu_p`	Pa·s	[0.001, 1e+12]	Plastic viscosity (Perzyna regularization)

class rheojax.models.ikh.mikh.MIKH[source]¶

Bases: IKHBase

Maxwell-Isotropic-Kinematic Hardening (MIKH) Model.

A thixotropic elasto-viscoplastic model combining: 1. Armstrong-Frederick kinematic hardening (backstress evolution). 2. Isotropic hardening/softening via structural parameter lambda (thixotropy). 3. Maxwell viscoelastic element for proper relaxation behavior. 4. Viscous background solvent.

Two Formulations:

Maxwell ODE (via Diffrax): For creep/relaxation protocols
Return Mapping: For startup/LAOS protocols (incremental)

Governing Equations:

σ_total = σ + η_inf * γ̇

Stress Evolution (ODE formulation): dσ/dt = G(γ̇ - γ̇ᵖ) - (G/η)σ

Yield Surface: |σ - α| ≤ σ_y(λ) σ_y(λ) = σ_y0 + Δσ_y * λ

Structure Evolution: dλ/dt = (1-λ)/τ_thix - Γ*λ*|γ̇ᵖ|

Backstress Evolution (Armstrong-Frederick): dα = C*dγ_p - γ_dyn*|α|^(m-1)*α*|dγ_p|

Parameters:

G – Shear modulus [Pa]
eta – Maxwell viscosity [Pa·s] (controls relaxation time τ = η/G)
C – Kinematic hardening modulus [Pa]
gamma_dyn – Dynamic recovery parameter for backstress [-]
m – AF recovery exponent [-] (typically 1.0)
sigma_y0 – Minimal (destructured) yield stress [Pa]
delta_sigma_y – Yield stress increment (structured - destructured) [Pa]
tau_thix – Thixotropic rebuilding time scale [s]
Gamma – Structural breakdown coefficient [-]
eta_inf – High-shear viscosity [Pa·s]
mu_p – Plastic viscosity for Perzyna regularization [Pa·s]

__init__()[source]¶: Initialize base model.

predict_flow_curve(gamma_dot)[source]¶

Predict steady-state flow curve.

Return type:: ndarray

predict_startup(t, gamma_dot=1.0)[source]¶

Predict startup shear response.

Parameters:

t (ndarray) – Time array
gamma_dot (float) – Constant shear rate

Return type:

Returns:

Stress vs time

predict_relaxation(t, sigma_0=100.0)[source]¶

Predict stress relaxation.

Parameters:

t (ndarray) – Time array
sigma_0 (float) – Initial stress

Return type:

Returns:

Stress vs time

predict_creep(t, sigma_applied=50.0)[source]¶

Predict creep response.

Parameters:

t (ndarray) – Time array
sigma_applied (float) – Applied constant stress

Return type:

Returns:

Strain vs time

predict_laos(t, gamma_0=1.0, omega=1.0)[source]¶

Predict LAOS response.

Parameters:

t (ndarray) – Time array
gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency

Return type:

Returns:

Stress vs time

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

NumPyro model function for Bayesian inference.

Accepts protocol-specific kwargs (gamma_dot, sigma_applied, sigma_0). Falls back to values cached during _fit() if not provided.

MLIKH¶

rheojax.models.ikh.ml_ikh.MLIKH | Handbook: Multi-Lambda Isotropic-Kinematic Hardening (ML-IKH) Multi-layer IKH model with multiple viscoelastic-thixotropic modes for broad relaxation spectra in complex yield-stress fluids.

Parameters (per mode `i`)¶
Parameter	Units	Bounds	Description
`G_i`	Pa	[0, 1e+09]	Mode i shear modulus
`eta_i`	Pa·s	[0.001, 1e+12]	Mode i Maxwell viscosity
`C_i`	Pa	[0, 1e+09]	Mode i kinematic hardening modulus
`gamma_dyn_i`	–	[0, 1e+04]	Mode i dynamic recovery parameter
`sigma_y0_i`	Pa	[0, 1e+09]	Mode i minimal yield stress
`delta_sigma_y_i`	Pa	[0, 1e+09]	Mode i structural yield stress
`tau_thix_i`	s	[1e-06, 1e+12]	Mode i thixotropic timescale
`Gamma_i`	–	[0, 1e+04]	Mode i breakdown coefficient

Global Parameters¶
Parameter	Units	Bounds	Description
`eta_inf`	Pa·s	[0, 1e+09]	High-shear viscosity (shared)

class rheojax.models.ikh.ml_ikh.MLIKH(n_modes=2, yield_mode='per_mode')[source]¶

Bases: IKHBase

Multi-Lambda Isotropic-Kinematic Hardening (ML-IKH) Model.

Extends MIKH to N modes connected in parallel. Each mode evolves its own internal variables (stress, backstress, structural lambda) with distinct timescales (tau_thix_i) and properties.

Two Yield Mode Options:

per_mode (default): Each mode has independent yield surface. Total stress is sum of mode stresses.
weighted_sum: Single global yield surface with structure contribution from all modes: σ_y = σ_y0 + k3·Σ(wᵢ·λᵢ)

Per-Mode Parameters (for each mode i=1..N):

G_i: Shear modulus C_i: Backstress modulus gamma_dyn_i: Dynamic recovery sigma_y0_i: Minimal yield stress delta_sigma_y_i: Structural yield stress tau_thix_i: Rebuilding timescale Gamma_i: Breakdown coefficient

Weighted-Sum Parameters:

G: Global shear modulus C: Global hardening modulus gamma_dyn: Global dynamic recovery sigma_y0: Base yield stress k3: Structure-yield coupling tau_thix_i: Per-mode rebuilding timescales Gamma_i: Per-mode breakdown coefficients w_i: Per-mode structure weights

Global Parameters (both yield modes):

eta_inf: High-shear viscosity mu_p: Plastic viscosity (Perzyna regularization) — controls creep/flow rate

Supported Protocols:

FLOW_CURVE: Steady-state stress vs shear rate (analytical solution)
STARTUP: Transient stress growth at constant shear rate (return mapping)
RELAXATION: Stress decay at constant strain (ODE formulation via Diffrax)
CREEP: Strain evolution at constant stress (ODE formulation via Diffrax)
OSCILLATION: Small amplitude oscillatory shear response
LAOS: Large amplitude oscillatory shear (return mapping with sinusoidal strain)

Note

Both yield modes (per_mode, weighted_sum) support all protocols. ODE protocols (creep, relaxation) use Diffrax for numerical integration. Return mapping protocols (startup, LAOS) use JAX scan for time stepping.

Parameters:

n_modes (int) – Number of structural modes (default: 2)
yield_mode (Literal['per_mode', 'weighted_sum']) – Yield formulation (‘per_mode’ or ‘weighted_sum’)

__init__(n_modes=2, yield_mode='per_mode')[source]¶

Initialize base model.

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

NumPyro model function with protocol-aware dispatch.

Accepts protocol-specific kwargs (gamma_dot, sigma_applied, sigma_0). Falls back to values cached during _fit() if not provided.

Parameters:

X – Input data
params – Parameter array or dict from NumPyro
test_mode – Optional protocol override
**kwargs – Protocol-specific arguments

Returns:

Predicted response

property n_modes: int¶: Number of structural modes.

property yield_mode: str¶: Yield formulation mode (‘per_mode’ or ‘weighted_sum’).

predict_flow_curve(gamma_dot)[source]¶

Predict steady-state flow curve.

Parameters:: gamma_dot (ndarray) – Array of shear rates
Return type:: ndarray
Returns:: Array of steady-state stresses

predict_startup(t, gamma_dot=1.0, strain=None)[source]¶

Predict startup shear response.

Parameters:

t (ndarray) – Time array
gamma_dot (float) – Applied shear rate (default: 1.0)
strain (ndarray | None) – Optional strain array (if None, uses gamma_dot * t)

Return type:

Returns:

Array of stresses

predict_relaxation(t, sigma_0=100.0)[source]¶

Predict stress relaxation after step strain.

Parameters:

t (ndarray) – Time array
sigma_0 (float) – Initial stress (default: 100.0)

Return type:

Returns:

Array of decaying stresses

predict_creep(t, sigma_applied=50.0)[source]¶

Predict creep response under constant stress.

Parameters:

t (ndarray) – Time array
sigma_applied (float) – Applied constant stress (default: 50.0)

Return type:

Returns:

Array of strains

predict_laos(t, gamma_0=1.0, omega=1.0)[source]¶

Predict large amplitude oscillatory shear response.

Parameters:

t (ndarray) – Time array
gamma_0 (float) – Strain amplitude (default: 1.0)
omega (float) – Angular frequency in rad/s (default: 1.0)

Return type:

Returns:

Array of stresses

FIKH Models (Fractional IKH)¶

FIKH¶

rheojax.models.fikh.fikh.FIKH | Handbook: Fractional Isotropic-Kinematic Hardening (FIKH) Fractional IKH model with Caputo fractional derivative for power-law memory in structure evolution (single-mode).

Parameters¶
Parameter	Units	Bounds	Description
`G`	Pa	[0.1, 1e+09]	Shear modulus
`eta`	Pa·s	[0.001, 1e+12]	Maxwell viscosity
`C`	Pa	[0, 1e+09]	Kinematic hardening modulus
`gamma_dyn`	–	[0, 1e+04]	Dynamic recovery parameter
`m`	–	[0.5, 3]	AF recovery exponent
`sigma_y0`	Pa	[0, 1e+09]	Minimal yield stress (destructured)
`delta_sigma_y`	Pa	[0, 1e+09]	Structural yield stress contribution
`tau_thix`	s	[1e-06, 1e+12]	Thixotropic rebuilding timescale
`Gamma`	–	[0, 1e+04]	Structural breakdown coefficient
`alpha_structure`	–	[0.01, 0.99]	Fractional order for structure evolution
`eta_inf`	Pa·s	[0, 1e+09]	High-shear viscosity
`mu_p`	Pa·s	[0.001, 1e+12]	Plastic viscosity (Perzyna regularization)

class rheojax.models.fikh.fikh.FIKH(include_thermal=True, include_isotropic_hardening=False, alpha_structure=0.5, n_history=100, stable_dt=0.01)[source]¶

Bases: FIKHBase

Fractional Isotropic-Kinematic Hardening (FIKH) Model.

A thixotropic elasto-viscoplastic model extending MIKH with: 1. Caputo fractional derivative for structure evolution (power-law memory). 2. Full thermokinematic coupling (Arrhenius + viscous heating).

The fractional derivative captures memory effects in thixotropic recovery, where the structure remembers its history via a power-law kernel rather than simple exponential decay.

Governing Equations:

σ_total = σ + η_inf·γ̇

Stress Evolution (ODE):: dσ/dt = G(γ̇ - γ̇ᵖ) - (G/η)σ
Yield Surface:: |σ - α| ≤ σ_y(λ, T) σ_y = σ_y0 + Δσ_y·λ^m_y · exp(E_y/R·(1/T - 1/T_ref))
Fractional Structure Evolution (Caputo):: D^α_C λ = (1-λ)/τ_thix - Γ·λ·|γ̇ᵖ|
Backstress (Armstrong-Frederick):: dα = C·dγᵖ - γ_dyn·|α|^(m-1)·α·|dγᵖ|
Temperature:: ρc_p·dT/dt = χ·σ·γ̇ᵖ - h·(T - T_env)

Parameters (22 with thermal):

G: Shear modulus [Pa] eta: Maxwell viscosity [Pa·s] C: Kinematic hardening modulus [Pa] gamma_dyn: Dynamic recovery parameter [-] m: AF recovery exponent [-] sigma_y0: Minimal yield stress [Pa] delta_sigma_y: Structural yield contribution [Pa] tau_thix: Thixotropic time scale [s] Gamma: Breakdown coefficient [-] alpha_structure: Fractional order (0 < α < 1) [-] eta_inf: High-shear viscosity [Pa·s] mu_p: Plastic viscosity [Pa·s] T_ref: Reference temperature [K] E_a: Viscosity activation energy [J/mol] E_y: Yield stress activation energy [J/mol] m_y: Structure exponent for yield [-] rho_cp: Volumetric heat capacity [J/(m³·K)] chi: Taylor-Quinney coefficient [-] h: Heat transfer coefficient [W/(m²·K)] T_env: Environmental temperature [K]

Limiting Behavior:

α → 1: Recovers classical IKH/MIKH (exponential structure relaxation) E_a = E_y = 0: Isothermal behavior (temperature-independent)

Example

>>> # Isothermal FIKH
>>> model = FIKH(include_thermal=False, alpha_structure=0.7)
>>> model.fit(omega, G_star, test_mode='oscillation')

>>> # Thermal FIKH with startup
>>> model = FIKH(include_thermal=True)
>>> result = model.fit(t, stress, test_mode='startup', strain=strain)
>>> sigma_pred = model.predict_startup(t_new, gamma_dot=1.0)

__init__(include_thermal=True, include_isotropic_hardening=False, alpha_structure=0.5, n_history=100, stable_dt=0.01)[source]¶

Initialize FIKH model.

Parameters:

include_thermal (bool) – Enable thermokinematic coupling (Arrhenius + heating).
include_isotropic_hardening (bool) – Enable isotropic hardening R.
alpha_structure (float) – Fractional order for structure (0 < α < 1). - α → 0: Strong memory (slow recovery) - α → 1: Weak memory (fast, exponential recovery)
n_history (int) – History buffer size for Caputo derivative.
stable_dt (float) – Internal integration substep (seconds) for startup / LAOS. See FIKHBase for the full explanation. Coarse user grids are densified to this step before the explicit return-mapping kernel runs. Set to 0 to disable. Default 0.02 s.

predict_flow_curve(gamma_dot, T=None)[source]¶

Predict steady-state flow curve.

Parameters:

gamma_dot (numpy.typing.ArrayLike) – Shear rate array.
T (float | None) – Temperature (if thermal enabled).

Return type:

numpy.typing.ArrayLike

Returns:

Stress array.

predict_startup(t, gamma_dot=1.0, T_init=None)[source]¶

Predict startup shear response.

Parameters:

t (numpy.typing.ArrayLike) – Time array.
gamma_dot (float) – Constant shear rate.
T_init (float | None) – Initial temperature.

Return type:

numpy.typing.ArrayLike

Returns:

Stress vs time.

predict_relaxation(t, sigma_0=100.0, T_init=None)[source]¶

Predict stress relaxation.

Parameters:

t (numpy.typing.ArrayLike) – Time array.
sigma_0 (float) – Initial stress.
T_init (float | None) – Initial temperature.

Return type:

numpy.typing.ArrayLike

Returns:

Stress vs time.

predict_creep(t, sigma_applied=50.0, T_init=None)[source]¶

Predict creep response.

Parameters:

t (numpy.typing.ArrayLike) – Time array.
sigma_applied (float) – Applied constant stress.
T_init (float | None) – Initial temperature.

Return type:

numpy.typing.ArrayLike

Returns:

Strain vs time.

predict_oscillation(omega, gamma_0=0.01, n_cycles=5)[source]¶

Predict oscillatory response (SAOS G*, G’, G’’).

For small amplitudes, this uses the linearized response. For accurate nonlinear response, use predict_laos().

Parameters:

omega (numpy.typing.ArrayLike) – Angular frequency array.
gamma_0 (float) – Strain amplitude (should be small).
n_cycles (int) – Number of cycles to simulate.

Return type:

Returns:

Complex modulus G* = G’ + i·G’’ for each frequency.

predict_laos(t, gamma_0=1.0, omega=1.0, T_init=None, strain=None)[source]¶

Predict LAOS (Large Amplitude Oscillatory Shear) response.

Integrates on the densified stable_dt grid (same as the fit path) so fit and predict stay in the linearization regime of the explicit return mapping. If strain is omitted, a clean sinusoid gamma_0·sin(omega·t) is synthesized; pass the measured strain array for fit/predict consistency on experimental data.

Parameters:

t (numpy.typing.ArrayLike) – Time array at which to report the response.
gamma_0 (float) – Strain amplitude used if strain is not given.
omega (float) – Angular frequency used if strain is not given.
T_init (float | None) – Initial temperature (thermal models only).
strain (Optional[numpy.typing.ArrayLike]) – Optional measured strain array aligned with t. When supplied, gamma_0 / omega are ignored for the simulation and the measured trace drives the return mapping.

Return type:

dict[str, Array]

Returns:

Dictionary with ‘time’, ‘strain’, ‘stress’, and optionally ‘temperature’.

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

Model function for NumPyro Bayesian inference.

This method provides a pure function interface for Bayesian sampling, capturing the test_mode via closure for correct mode-aware inference.

Parameters:

X (numpy.typing.ArrayLike) – Input data.
params (Union[numpy.typing.ArrayLike, dict[str, Any]]) – Parameter array or dictionary.
test_mode (str | None) – Protocol (uses stored _test_mode if None).
**kwargs – Protocol-specific arguments (e.g., strain, sigma_0).

Return type:

Returns:

Predicted values.

get_limiting_behavior()[source]¶

Get limiting behavior diagnostics.

Return type:: dict[str, Any]
Returns:: Dictionary with limiting cases and expected behavior.

precompile(test_mode='relaxation', X=None, y=None, *, n_points=100, verbose=True)[source]¶

Precompile JIT kernels for faster subsequent predictions.

Triggers JAX JIT compilation of the core FIKH kernels by running a small dummy prediction. This is useful when you want to avoid the compilation overhead on first real prediction.

Parameters:

test_mode (str) – Accepted for parent compatibility (unused).
X – Accepted for parent compatibility (unused).
y – Accepted for parent compatibility (unused).
n_points (int) – Number of time points for dummy data.
verbose (bool) – Print compilation time if True.

Return type:

float

Returns:

Compilation time in seconds.

Example

>>> model = FIKH(include_thermal=True)
>>> compile_time = model.precompile()  # Triggers JIT
>>> # Now predictions will be fast
>>> sigma = model.predict_startup(t_real, gamma_dot=1.0)

__repr__()[source]¶

String representation.

Return type:: str

FMLIKH¶

rheojax.models.fikh.fmlikh.FMLIKH | Handbook: Fractional Multi-Lambda IKH (FMLIKH) Fractional multi-layer IKH model combining multi-mode viscoelasticity with Caputo fractional structure kinetics.

Parameters (per mode `i`)¶
Parameter	Units	Bounds	Description
`G_i`	Pa	[0.1, 1e+09]	Mode i shear modulus
`eta_i`	Pa·s	[0.001, 1e+12]	Mode i Maxwell viscosity
`C_i`	Pa	[0, 1e+09]	Mode i kinematic hardening modulus
`gamma_dyn_i`	–	[0, 1e+04]	Mode i dynamic recovery parameter

Shared Parameters¶
Parameter	Units	Bounds	Description
`sigma_y0`	Pa	[0, 1e+09]	Minimal yield stress
`delta_sigma_y`	Pa	[0, 1e+09]	Structural yield stress contribution
`tau_thix`	s	[1e-06, 1e+12]	Thixotropic rebuilding timescale
`Gamma`	–	[0, 1e+04]	Structural breakdown coefficient
`alpha_structure`	–	[0.01, 0.99]	Fractional order for structure evolution
`eta_inf`	Pa·s	[0, 1e+09]	High-shear viscosity
`mu_p`	Pa·s	[0.001, 1e+12]	Plastic viscosity

class rheojax.models.fikh.fmlikh.FMLIKH(n_modes=3, include_thermal=True, include_isotropic_hardening=False, shared_alpha=True, alpha_structure=0.5, n_history=100, stable_dt=0.01)[source]¶

Bases: FIKHBase

Fractional Multi-Layer Isotropic-Kinematic Hardening (FMLIKH) Model.

A multi-mode extension of FIKH supporting multiple viscoelastic relaxation processes with shared yield and thixotropy behavior.

The total stress is the sum of contributions from N modes:

σ_total = Σᵢ σᵢ + η_inf·γ̇

Each mode i has its own:

Gᵢ: Shear modulus
ηᵢ: Viscosity (defines τᵢ = ηᵢ/Gᵢ)
Cᵢ: Kinematic hardening modulus
γ_dyn,ᵢ: Dynamic recovery parameter

Shared across all modes:

σ_y0, δσ_y: Yield stress parameters
τ_thix, Γ: Thixotropy parameters
α: Fractional order (or per-mode if shared_alpha=False)
Thermal parameters

This structure allows capturing materials with:

Multiple relaxation time scales
Complex SAOS signatures (wide frequency range)
Non-trivial startup overshoot dynamics

Parameters:

n_modes (int) – Number of viscoelastic modes.
shared_alpha (bool) – If True, use single α for all modes. If False, α_i per mode.
FIKHBase. (Other parameters inherited from)

Example

>>> model = FMLIKH(n_modes=3, include_thermal=True, shared_alpha=True)
>>> model.fit(t, stress, test_mode='startup', strain=strain)

>>> # Access per-mode parameters
>>> G_values = [model.parameters.get_value(f'G_{i}') for i in range(3)]

__init__(n_modes=3, include_thermal=True, include_isotropic_hardening=False, shared_alpha=True, alpha_structure=0.5, n_history=100, stable_dt=0.01)[source]¶

Initialize FMLIKH model.

Parameters:

n_modes (int) – Number of viscoelastic modes (≥1).
include_thermal (bool) – Enable thermokinematic coupling.
include_isotropic_hardening (bool) – Enable isotropic hardening R.
shared_alpha (bool) – Use single fractional order (True) or per-mode (False).
alpha_structure (float) – Fractional order (used if shared_alpha=True).
n_history (int) – History buffer size.
stable_dt (float) – Internal integration substep for startup / LAOS return mapping. See FIKHBase for details.

property n_modes: int¶: Number of viscoelastic modes.

predict_oscillation(omega, gamma_0=0.01, n_cycles=5)[source]¶

Predict oscillatory response (SAOS G*) for multi-mode model.

Parameters:

omega (numpy.typing.ArrayLike) – Angular frequency array.
gamma_0 (float) – Strain amplitude.
n_cycles (int) – Number of cycles to simulate.

Return type:

Returns:

Complex modulus G* = G’ + i·G’’ for each frequency.

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

Model function for Bayesian inference.

Return type:: Array

get_mode_info()[source]¶

Get information about each mode.

Return type:: dict[str, Any]
Returns:: Dictionary with per-mode parameters and derived quantities.

__repr__()[source]¶

String representation.

Return type:: str

Hébraud-Lequeux Model¶

HebraudLequeux¶

rheojax.models.hl.hebraud_lequeux.HebraudLequeux | Handbook: Hébraud–Lequeux (HL) Model — Handbook Mean-field stochastic model for soft glassy materials with stress-induced yielding and mechanical noise coupling.

Parameters¶
Parameter	Units	Bounds	Description
`alpha`	–	[0, 1]	Coupling parameter (\(\alpha < 0.5\) implies yield stress)
`tau`	s	[1e-06, 1e+04]	Microscopic yield timescale
`sigma_c`	Pa	[0.001, 1e+06]	Critical yield stress threshold

class rheojax.models.hl.hebraud_lequeux.HebraudLequeux[source]¶

Bases: BaseModel

Hébraud–Lequeux (HL) Model for Soft Glassy Materials.

The HL model (1998) is a mean-field description of yield-stress fluids where mesoscopic blocks of stress evolve via elastic loading, plastic yielding, and stress diffusion (mechanical noise) generated by yielding events elsewhere.

It predicts: - Finite yield stress for coupling parameter alpha < 0.5 - Herschel-Bulkley flow curves (stress ~ gdot^0.5) near yield - Creep and delayed yielding - Stress overshoots in startup flow - Non-linear LAOS response

Parameters:

alpha – Coupling parameter (dimensionless). Controls phase state. alpha < 0.5: Glassy (yield stress) alpha >= 0.5: Fluid (no yield stress)
tau – Microscopic yield timescale (s).
sigma_c – Critical yield stress threshold (Pa).

parameters¶: ParameterSet containing alpha, tau, sigma_c.

__init__()[source]¶: Initialize Hébraud–Lequeux Model.

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

Model function for Bayesian inference (NumPyro NUTS).

Parameters:

X (ndarray) – Input array
params (ndarray) – Parameter values [alpha, tau, sigma_c]
test_mode (str | None) – Override test mode
**kwargs – Protocol kwargs forwarded by BayesianMixin. Falls back to _last_fit_kwargs for values not provided.

get_phase_state()[source]¶

Return the phase state based on alpha.

Return type:: str

Giesekus Models¶

GiesekusSingleMode¶

rheojax.models.giesekus.single_mode.GiesekusSingleMode | Handbook: Giesekus Model — Handbook Nonlinear differential constitutive model with anisotropic mobility for shear-thinning polymer solutions and melts.

Parameters¶
Parameter	Units	Bounds	Description
`eta_p`	Pa·s	[0.001, 1e+06]	Polymer viscosity (zero-shear polymer contribution)
`lambda_1`	s	[1e-06, 1e+04]	Relaxation time
`alpha`	–	[0, 0.5]	Mobility factor (0 = UCM, 0.5 = max anisotropy)
`eta_s`	Pa·s	[0, 1e+04]	Solvent viscosity (Newtonian contribution)

class rheojax.models.giesekus.single_mode.GiesekusSingleMode[source]¶

Bases: GiesekusBase

Single-mode Giesekus nonlinear viscoelastic model.

The Giesekus model extends the Upper-Convected Maxwell framework with a quadratic stress term representing anisotropic molecular mobility:

τ + λ∇̂τ + (αλ/η_p)τ·τ = 2η_p D

This captures:

Shear-thinning: Viscosity decreases with increasing shear rate
Normal stresses: Both N₁ > 0 and N₂ < 0 (with N₂/N₁ = -α/2)
Stress overshoot: Peak stress in startup flow
Nonlinear LAOS: Higher harmonics in large-amplitude oscillation

Parameters:

eta_p (float) – Polymer viscosity η_p (Pa·s). Default: 100.0
lambda_1 (float) – Relaxation time λ (s). Default: 1.0
alpha (float) – Mobility factor α (dimensionless, 0 ≤ α ≤ 0.5). Default: 0.3
eta_s (float) – Solvent viscosity η_s (Pa·s). Default: 0.0

parameters¶

Model parameters for fitting

Type:: ParameterSet

fitted_¶

Whether the model has been fitted

Type:: bool

Examples

Basic fitting and prediction:

>>> model = GiesekusSingleMode()
>>> # Generate synthetic data
>>> gamma_dot = np.logspace(-2, 2, 50)
>>> sigma_data = model.predict(gamma_dot, test_mode='flow_curve')
>>> # Fit to data
>>> model.fit(gamma_dot, sigma_data, test_mode='flow_curve')

Bayesian inference:

>>> result = model.fit_bayesian(gamma_dot, sigma_data, test_mode='flow_curve')
>>> print(result.diagnostics)

See also

GiesekusMultiMode: Multi-mode extension with N relaxation times

__init__()[source]¶: Initialize single-mode Giesekus model.

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

NumPyro/BayesianMixin model function.

Routes to appropriate prediction based on test_mode.

Parameters:

X (array-like) – Independent variable
params (array-like) – Parameter values [eta_p, lambda_1, alpha, eta_s, beta_cc]
test_mode (str, optional) – Override stored test mode
**kwargs (dict) – Protocol-specific arguments (gamma_dot, sigma_applied, gamma_0, etc.)

Returns:

Predicted response

Return type:

jnp.ndarray

predict_flow_curve(gamma_dot, return_components=False)[source]¶

Predict steady shear stress and viscosity.

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 viscoelastic regime, the Giesekus model reduces to Maxwell behavior (α-independent).

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:

predict_normal_stresses(gamma_dot)[source]¶

Predict first and second normal stress differences.

The Giesekus model predicts:

N₁ = τ_xx - τ_yy > 0  (first normal stress difference)
N₂ = τ_yy - τ_zz < 0  (second normal stress difference)

with the diagnostic ratio N₂/N₁ = -α/2.

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.

The Giesekus model predicts stress overshoot in startup flow, where the stress first exceeds then relaxes to its steady-state value.

Parameters:

t (ndarray) – Time array (s)
gamma_dot (float) – Applied shear rate (1/s)
return_full (bool) – If True, return full stress tensor components

Returns:

Shear stress τ_xy(t), or (τ_xx, τ_yy, τ_xy, τ_zz) if return_full=True

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.

The Giesekus model predicts faster-than-exponential relaxation due to the quadratic τ·τ term.

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 stress tensor components

Returns:

Relaxing stress τ_xy(t), or (τ_xx, τ_yy, τ_xy, τ_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:

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

Simulate Large-Amplitude Oscillatory Shear (LAOS).

The Giesekus model produces nonlinear stress response in LAOS, characterized by higher harmonics (I₃, I₅, …).

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’: Time array - ‘strain’: γ(t) = γ₀·sin(ωt) - ‘stress’: τ_xy(t) - ‘strain_rate’: γ̇(t) = γ₀·ω·cos(ωt)

Return type:

extract_laos_harmonics(laos_result, n_harmonics=5)[source]¶

Extract Fourier harmonics from LAOS stress response.

The nonlinear stress response can be decomposed as:

σ(t) = Σ [σ'_n·sin(nωt) + σ''_n·cos(nωt)]

Only odd harmonics (n=1,3,5,…) are expected for symmetric response.

Parameters:

laos_result (dict[str, ndarray]) – Result from simulate_laos()
n_harmonics (int) – Number of harmonics to extract

Returns:

Dictionary with: - ‘n’: Harmonic indices [1, 3, 5, …] - ‘sigma_prime’: In-phase (elastic) components - ‘sigma_double_prime’: Out-of-phase (viscous) components - ‘intensity’: |σ_n| = sqrt(σ’_n² + σ’’_n²)

Return type:

get_overshoot_ratio(gamma_dot, t_max=None)[source]¶

Compute stress overshoot ratio in startup flow.

The overshoot ratio is defined as:

overshoot_ratio = σ_max / σ_ss

where σ_max is the peak stress and σ_ss is the steady-state stress.

Parameters:

gamma_dot (float) – Shear rate (1/s)
t_max (float | None) – Maximum simulation time (default: 20λ)

Returns:

(overshoot_ratio, strain_at_overshoot)

Return type:

tuple[float, float]

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

Get effective relaxation spectrum G(t).

Note: For single-mode Giesekus, this is single-exponential in the linear regime but deviates due to nonlinearity.

Parameters:

t (ndarray | None) – Time array (default: logspace from 0.01λ to 100λ)
n_points (int) – Number of points if t not provided

Returns:

(t, G(t))

Return type:

GiesekusMultiMode¶

rheojax.models.giesekus.multi_mode.GiesekusMultiMode | Handbook: Giesekus Model — Handbook Multi-mode Giesekus model for broad relaxation spectra in polymer systems.

Parameters (per mode `i`)¶
Parameter	Units	Bounds	Description
`eta_p_i`	Pa·s	[0.001, 1e+06]	Mode i polymer viscosity
`lambda_1_i`	s	[1e-06, 1e+04]	Mode i relaxation time
`alpha_i`	–	[0, 0.5]	Mode i mobility factor

Global Parameters¶
Parameter	Units	Bounds	Description
`eta_s`	Pa·s	[0, 1e+04]	Solvent viscosity (shared)

class rheojax.models.giesekus.multi_mode.GiesekusMultiMode(n_modes=3)[source]¶

Bases: BaseModel

Multi-mode Giesekus nonlinear viscoelastic model.

This model extends the single-mode Giesekus with N parallel Maxwell modes, each with its own relaxation time, viscosity, and mobility factor.

The constitutive equation for each mode is:

τᵢ + λᵢ∇̂τᵢ + (αᵢλᵢ/η_p,i)τᵢ·τᵢ = 2η_p,i D

Total stress: σ = η_s·γ̇ + Σᵢ τᵢ

Parameters:: n_modes (int) – Number of relaxation modes (N ≥ 1). Default: 3

parameters¶

Model parameters including per-mode values

Type:: ParameterSet

fitted_¶

Whether the model has been fitted

Type:: bool

Notes

The multi-mode model is particularly useful for:

Fitting broad SAOS spectra that single-mode cannot capture
Representing polydisperse polymer systems
Capturing multiple relaxation processes

Each mode can have different α_i values, allowing different molecular weight fractions to exhibit different anisotropy.

See also

GiesekusSingleMode: Single relaxation time variant
GeneralizedMaxwell: Linear multi-mode Maxwell model

__init__(n_modes=3)[source]¶

Initialize multi-mode Giesekus model.

Parameters:: n_modes (int) – Number of relaxation modes (must be ≥ 1)
Raises:: ValueError – If n_modes < 1

property n_modes: int¶: Get number of modes.

property eta_s: float¶: Get solvent viscosity η_s (Pa·s).

property eta_0: float¶: Get zero-shear viscosity η₀ = η_s + Σ η_p,i (Pa·s).

get_mode_params(mode_idx)[source]¶

Get parameters for a specific mode.

Parameters:: mode_idx (int) – Mode index (0 to n_modes-1)
Returns:: Dictionary with keys ‘eta_p’, ‘lambda_1’, ‘alpha’
Return type:: dict[str, float]

set_mode_params(mode_idx, eta_p=None, lambda_1=None, alpha=None)[source]¶

Set parameters for a specific mode.

Parameters:

mode_idx (int) – Mode index (0 to n_modes-1)
eta_p (float | None) – Polymer viscosity (Pa·s)
lambda_1 (float | None) – Relaxation time (s)
alpha (float | None) – Mobility factor (0 ≤ α ≤ 0.5)

Return type:

get_mode_arrays()[source]¶

Get all mode parameters as JAX arrays.

Uses vectorized extraction via get_values() + slicing for ~3x speedup over N individual get_value() calls.

Returns:: (eta_p_modes, lambda_modes, alpha_modes), each shape (n_modes,)
Return type:: tuple[Array, Array, Array]

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

NumPyro/BayesianMixin model function.

Parameters:

X (array-like) – Independent variable
params (array-like) – All parameter values in order
test_mode (str, optional) – Override stored test mode
**kwargs (dict) – Protocol-specific arguments (gamma_dot, sigma_applied, etc.)

Returns:

Predicted response

Return type:

jnp.ndarray

predict_saos(omega, return_components=True)[source]¶

Predict SAOS storage and loss moduli.

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:

predict_flow_curve(gamma_dot, return_components=False)[source]¶

Predict steady shear stress.

Parameters:

gamma_dot (ndarray) – Shear rate array (1/s)
return_components (bool) – If True, return (sigma, eta)

Returns:

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

Return type:

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 per-mode stresses

Returns:

Total shear stress, or dict with per-mode stresses

Return type:

get_relaxation_spectrum()[source]¶

Get discrete relaxation spectrum.

Returns:: (lambda_i, G_i) where G_i = η_p,i / λ_i
Return type:: tuple[ndarray, ndarray]

get_continuous_spectrum(t=None, n_points=200)[source]¶

Get continuous relaxation modulus G(t).

Parameters:

t (ndarray | None) – Time array
n_points (int) – Number of points if t not provided

Returns:

(t, G(t))

Return type:

__repr__()[source]¶

Return string representation.

Return type:: str

ITT-MCT Models (Integration Through Transients MCT)¶

ITTMCTSchematic¶

rheojax.models.itt_mct.schematic.ITTMCTSchematic | Handbook: ITT-MCT Schematic (F_1_2) F₁₂ schematic MCT model for dense colloidal suspensions and glassy materials with memory kernel \(m(\Phi) = v_1\Phi + v_2\Phi^2\).

Parameters¶
Parameter	Units	Bounds	Description
`v1`	–	[0, 5]	Linear vertex coefficient (typically 0 for F₁₂)
`v2`	–	[0.5, 10]	Quadratic vertex coefficient (glass at \(v_2 > 4\))
`Gamma`	1/s	[1e-06, 1e+06]	Bare relaxation rate
`gamma_c`	–	[0.01, 0.5]	Critical strain for cage breaking
`G_inf`	Pa	[1, 1e+12]	High-frequency elastic modulus

class rheojax.models.itt_mct.schematic.ITTMCTSchematic(epsilon=None, v2=None, integration_method='volterra', n_prony_modes=10, decorrelation_form='gaussian', memory_form='simplified', stress_form='schematic', phi_volume=None, k_BT=1.0)[source]¶

Bases: ITTMCTBase

ITT-MCT Schematic F₁₂ Model.

The F₁₂ model uses a quadratic memory kernel m(Φ) = v₁Φ + v₂Φ² to describe the cage effect in dense colloidal suspensions.

The glass transition occurs when the non-ergodicity parameter f (long-time limit of Φ) becomes non-zero, which happens at v₂ = v₂_c = 4 for v₁ = 0.

Parameters:

epsilon (float | None) – Separation parameter. If provided, v₂ is set to achieve this ε. ε < 0: fluid state ε = 0: critical point ε > 0: glass state
v2 (float | None) – Quadratic vertex coefficient. Alternative to epsilon.
integration_method (Literal['volterra', 'history']) – Integration method for memory kernel
n_prony_modes (int) – Number of Prony modes for Volterra integration
decorrelation_form (Literal['gaussian', 'lorentzian']) – Strain decorrelation function form: - “gaussian”: h(γ) = exp(-(γ/γ_c)²) - faster exponential decay - “lorentzian”: h(γ) = 1/(1+(γ/γ_c)²) - slower algebraic decay
memory_form (Literal['simplified', 'full']) – Memory kernel form: - “simplified”: single decorrelation m(Φ) = h[γ_acc] × (v₁Φ + v₂Φ²) - “full”: two-time decorrelation m(t,s,t₀) = h[γ(t,t₀)] × h[γ(t,s)] × (v₁Φ + v₂Φ²)
stress_form (Literal['schematic', 'microscopic']) – Stress computation form: - “schematic”: σ = G_∞ × γ̇ × ∫ Φ² × h(γ) dt (standard schematic) - “microscopic”: σ = (k_BT/60π²) × ∫dk k⁴ [S’/S²]² Φ² (structure factor weighted)
phi_volume (float | None) – Volume fraction for Percus-Yevick S(k). Required if stress_form=”microscopic”.
k_BT (float) – Thermal energy k_B × T in Joules. Default 1.0 gives dimensionless stress.

parameters¶

Model parameters with the following: - v1: Linear vertex (default 0) - v2: Quadratic vertex (default 2.0, fluid state) - Gamma: Bare relaxation rate (default 1.0 s⁻¹) - gamma_c: Critical strain (default 0.1) - G_inf: High-frequency modulus (default 1e6 Pa)

Type:: ParameterSet

Examples

>>> model = ITTMCTSchematic(epsilon=-0.1)  # Fluid state
>>> model.get_glass_transition_info()
{'is_glass': False, 'epsilon': -0.1, ...}

>>> model = ITTMCTSchematic(epsilon=0.05)  # Glass state
>>> sigma = model.predict(np.logspace(-3, 2, 50), test_mode='flow_curve')
>>> # Shows yield stress at low rates

>>> # Use Lorentzian decorrelation for materials with extended yielding
>>> model = ITTMCTSchematic(epsilon=0.05, decorrelation_form="lorentzian")

>>> # Use full two-time memory kernel (Fuchs & Cates 2002)
>>> model = ITTMCTSchematic(epsilon=0.05, memory_form="full")

>>> # Use microscopic stress with Percus-Yevick S(k)
>>> model = ITTMCTSchematic(
...     epsilon=0.05,
...     stress_form="microscopic",
...     phi_volume=0.5,
...     k_BT=4.11e-21,  # Room temperature
... )

__init__(epsilon=None, v2=None, integration_method='volterra', n_prony_modes=10, decorrelation_form='gaussian', memory_form='simplified', stress_form='schematic', phi_volume=None, k_BT=1.0)[source]¶

Initialize F₁₂ Schematic Model.

Parameters:

epsilon (float | None) – Separation parameter ε = (v₂ - v₂_c)/v₂_c. Mutually exclusive with v2.
v2 (float | None) – Direct vertex coefficient. Mutually exclusive with epsilon.
integration_method (Literal['volterra', 'history']) – Integration method for memory kernel
n_prony_modes (int) – Number of Prony modes
decorrelation_form (Literal['gaussian', 'lorentzian']) – Form of the strain decorrelation function h(γ): - “gaussian”: h(γ) = exp(-(γ/γ_c)²) - faster decay (default, Fuchs & Cates 2002) - “lorentzian”: h(γ) = 1/(1+(γ/γ_c)²) - slower algebraic decay (Brader et al. 2008)
memory_form (Literal['simplified', 'full']) – Memory kernel form: - “simplified”: single decorrelation m(Φ) = h[γ_acc] × (v₁Φ + v₂Φ²) - “full”: two-time decorrelation m(t,s,t₀) = h[γ(t,t₀)] × h[γ(t,s)] × (v₁Φ + v₂Φ²)
stress_form (Literal['schematic', 'microscopic']) – Stress computation form: - “schematic”: σ = G_∞ × γ̇ × ∫ Φ² × h(γ) dt (standard schematic) - “microscopic”: σ = (k_BT/60π²) × ∫dk k⁴ [S’/S²]² Φ² (structure factor weighted)
phi_volume (float | None) – Volume fraction for Percus-Yevick S(k). Required if stress_form=”microscopic”.
k_BT (float) – Thermal energy k_B × T in Joules. Default 1.0 gives dimensionless stress. Use 4.11e-21 J for T=298K with real units.

get_glass_transition_info()[source]¶

Get information about the glass transition state.

Returns:: Glass transition properties: - is_glass: bool - epsilon: separation parameter - v2_critical: critical v₂ value - f_neq: non-ergodicity parameter - lambda_exponent: MCT exponent parameter
Return type:: dict[str, Any]

property epsilon: float¶: Get separation parameter ε = (v₂ - v₂_c)/v₂_c.

get_laos_harmonics(t, gamma_0=0.1, omega=1.0, n_harmonics=5)[source]¶

Extract Fourier harmonics from LAOS response.

Parameters:

t (ndarray) – Time array covering at least one full period
gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency
n_harmonics (int) – Number of odd harmonics to extract

Return type:

Returns:

sigma_prime_n (np.ndarray) – In-phase coefficients [σ’₁, σ’₃, σ’₅, …]
sigma_double_prime_n (np.ndarray) – Out-of-phase coefficients [σ’’₁, σ’’₃, σ’’₅, …]

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

Static model function for Bayesian inference.

NOTE: Bayesian inference is not yet supported for ITT-MCT models. The Prony decomposition depends on parameters (v1, v2) and would need to be recomputed for each MCMC sample, which is computationally prohibitive.

Use NLSQ fitting with bootstrap resampling for uncertainty quantification.

Parameters:

X (ndarray) – Independent variable
params (ndarray) – Parameter array [v1, v2, Gamma, gamma_c, G_inf] in parameter order
test_mode (str) – Protocol type
**kwargs – Additional protocol-specific parameters

Raises:

NotImplementedError – Always raised - Bayesian inference not supported for ITT-MCT models

Return type:

precompile(test_mode='relaxation', X=None, y=None, **kwargs)[source]¶

Pre-compile the diffrax ODE solver for fast subsequent calls.

Triggers JIT compilation with dummy data so the first real prediction doesn’t incur the compilation cost. Useful when predictable timing is important (e.g., in interactive applications or benchmarks).

Returns:: Compilation time in seconds (0.0 if diffrax not available)
Return type:: float

Examples

>>> model = ITTMCTSchematic(epsilon=0.05)
>>> compile_time = model.precompile()
>>> print(f"Compilation took {compile_time:.1f}s")
>>> # Now flow curve predictions will be fast
>>> sigma = model.predict(gamma_dot, test_mode='flow_curve')

Notes

First call to flow curve prediction triggers JIT compilation which can take 30-90 seconds. This method triggers that compilation upfront.

Only affects diffrax-based flow curve solver. Other protocols (oscillation, startup, etc.) use scipy and don’t need precompilation.

property decorrelation_form: str¶: Get the strain decorrelation function form.

property memory_form: str¶

Get the memory kernel form.

Returns:: “simplified” for single decorrelation m(Φ) = h[γ_acc] × (v₁Φ + v₂Φ²) “full” for two-time decorrelation m(t,s,t₀) = h[γ(t,t₀)] × h[γ(t,s)] × (v₁Φ + v₂Φ²)
Return type:: str

property stress_form: str¶

Get the stress computation form.

Returns:: “schematic” for σ = G_∞ × γ̇ × ∫ Φ² × h(γ) dt “microscopic” for σ = (k_BT/60π²) × ∫dk k⁴ [S’/S²]² Φ²
Return type:: str

__repr__()[source]¶

Return string representation.

Return type:: str

ITTMCTIsotropic¶

rheojax.models.itt_mct.isotropic.ITTMCTIsotropic | Handbook: ITT-MCT Isotropic (ISM) Isotropic Schematic Model (ISM) with Percus-Yevick structure factor \(S(k)\) for quantitative predictions from microscopic parameters.

Parameters¶
Parameter	Units	Bounds	Description
`phi`	–	[0.1, 0.64]	Volume fraction (glass at \(\phi \approx 0.516\))
`sigma_d`	m	[1e-09, 0.001]	Particle diameter
`D0`	m²/s	[1e-18, 1e-06]	Bare short-time diffusion coefficient
`kBT`	J	[1e-24, 1e-18]	Thermal energy \(k_B T\)
`gamma_c`	–	[0.01, 0.5]	Critical strain for cage breaking

class rheojax.models.itt_mct.isotropic.ITTMCTIsotropic(phi=None, sk_source='percus_yevick', k_data=None, sk_data=None, n_k=100, integration_method='volterra', n_prony_modes=10)[source]¶

Bases: ITTMCTBase

ITT-MCT Isotropically Sheared Model with k-resolved correlators.

The ISM computes density correlators Φ(k,t) for an array of wave vectors, using the static structure factor S(k) to compute the MCT memory kernel.

The model can use: - Built-in Percus-Yevick S(k) for hard spheres (default) - User-provided S(k) data

Parameters:

phi (float | None) – Volume fraction. If provided with Percus-Yevick, determines S(k).
sk_source (Literal['percus_yevick', 'user_provided']) – Source of structure factor data
k_data (ndarray | None) – Wave vectors for user-provided S(k)
sk_data (ndarray | None) – Structure factor values for user-provided S(k)
n_k (int) – Number of k-grid points
integration_method (Literal['volterra', 'history']) – Integration method for memory kernel

k_grid¶

Wave vector array (1/m or dimensionless)

Type:: np.ndarray

S_k¶

Structure factor at k_grid points

Type:: np.ndarray

vertex¶

MCT vertex matrix V(k,q)

Type:: np.ndarray

Examples

>>> # Using Percus-Yevick for hard spheres
>>> model = ITTMCTIsotropic(phi=0.55)
>>> model.get_glass_transition_info()
{'is_glass': True, 'phi': 0.55, 'phi_mct': 0.516, ...}

>>> # Using user-provided S(k)
>>> model = ITTMCTIsotropic(
...     sk_source="user_provided",
...     k_data=k_experimental,
...     sk_data=sk_experimental
... )

__init__(phi=None, sk_source='percus_yevick', k_data=None, sk_data=None, n_k=100, integration_method='volterra', n_prony_modes=10)[source]¶

Initialize ISM model.

Parameters:

phi (float | None) – Volume fraction for Percus-Yevick S(k)
sk_source (Literal['percus_yevick', 'user_provided']) – Source of structure factor
k_data (ndarray | None) – User-provided structure factor data
sk_data (ndarray | None) – User-provided structure factor data
n_k (int) – Number of k-grid points
integration_method (Literal['volterra', 'history']) – Integration method
n_prony_modes (int) – Number of Prony modes

update_structure_factor(phi=None, k_data=None, sk_data=None)[source]¶

Update structure factor (e.g., after parameter change).

Parameters:

phi (float | None) – New volume fraction (for Percus-Yevick)
k_data (ndarray | None) – New user-provided S(k) data
sk_data (ndarray | None) – New user-provided S(k) data

Return type:

get_glass_transition_info()[source]¶

Get information about the glass transition state.

Returns:: Glass transition properties including: - is_glass: bool - phi: current volume fraction - phi_mct: MCT glass transition (≈0.516) - S_max: peak of S(k)
Return type:: dict[str, Any]

get_sk_info()[source]¶

Get information about current S(k).

Returns:: S(k) properties
Return type:: dict[str, Any]

get_laos_harmonics(t, gamma_0=0.1, omega=1.0, n_harmonics=5)[source]¶

Extract Fourier harmonics from LAOS response.

Parameters:

t (ndarray) – Time array covering at least one full period
gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency
n_harmonics (int) – Number of odd harmonics to extract

Return type:

Returns:

sigma_prime_n (np.ndarray) – In-phase coefficients [sigma’_1, sigma’_3, sigma’_5, …]
sigma_double_prime_n (np.ndarray) – Out-of-phase coefficients [sigma’’_1, sigma’’_3, sigma’’_5, …]

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

Static model function for Bayesian inference.

NOTE: Bayesian inference is not yet supported for ITT-MCT Isotropic models. The ISM kernel requires full k-grid integration and Prony decomposition for each MCMC sample, making NUTS sampling computationally prohibitive.

Use NLSQ fitting with bootstrap resampling for uncertainty quantification.

Raises:: NotImplementedError – Always raised - Bayesian inference not supported for ITT-MCT Isotropic models
Return type:: ndarray

__repr__()[source]¶

Return string representation.

Return type:: str

TNT Models (Transient Network Theory)¶

TNTSingleMode¶

rheojax.models.tnt.single_mode.TNTSingleMode | Handbook: TNT Bell (Force-Dependent Breakage) — Handbook Single-mode transient network model with configurable breakage kinetics (Bell, power-law, stretch-creation) and optional FENE-P springs.

Parameters¶
Parameter	Units	Bounds	Description
`G`	Pa	[1, 1e+08]	Network modulus (active chain density)
`tau_b`	s	[1e-06, 1e+04]	Bond lifetime (mean detachment time)
`eta_s`	Pa·s	[0, 1e+04]	Solvent viscosity (Newtonian background)
`nu`	–	[0.01, 20]	Bell force sensitivity (higher = more shear-thinning)
`L_max`	–	[2, 100]	Maximum chain extensibility (FENE-P, if enabled)

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

Bases: TNTBase

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:

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:

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:

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:

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:

TNTLoopBridge¶

rheojax.models.tnt.loop_bridge.TNTLoopBridge | Handbook: TNT Loop-Bridge (Two-Species Kinetics) — Handbook Loop-bridge transient network model with stress-dependent loop-to-bridge conversion.

Parameters¶
Parameter	Units	Bounds	Description
`G`	Pa	[1, 1e+08]	Network modulus (fully bridged state)
`tau_b`	s	[1e-06, 1e+04]	Bridge detachment time
`tau_a`	s	[1e-06, 1e+04]	Loop attachment time
`nu`	–	[0.01, 20]	Bell force sensitivity
`f_B_eq`	–	[0.01, 0.99]	Equilibrium bridge fraction at rest
`eta_s`	Pa·s	[0, 1e+04]	Solvent viscosity

class rheojax.models.tnt.loop_bridge.TNTLoopBridge[source]¶

Bases: TNTBase

Loop-bridge kinetics model for telechelic polymer networks.

Implements reversible loop-bridge kinetics for telechelic polymers where chains can exist as load-bearing bridges (both ends on different junctions) or non-load-bearing loops (both ends on same junction).

The bridge fraction f_B evolves dynamically via attachment (loops → bridges) and force-activated dissociation (bridges → loops via Bell kinetics).

State Variables¶

f_B: Bridge fraction (0 to 1)
S: Conformation tensor of bridges [S_xx, S_yy, S_zz, S_xy]

Key Equations¶

Bridge fraction kinetics:

df_B/dt = (1 - f_B)/τ_a - f_B·β(stretch)
β(stretch) = (1/τ_b)·exp(ν·(stretch - 1))

Conformation tensor (bridges only):

dS/dt = L·S + S·L^T + g_0·I - β(stretch)·S

Stress (only bridges carry load):

σ = f_B·G·S_xy + η_s·γ̇

param G:: Network modulus (fully bridged, Pa), bounds (1e0, 1e8)
type G:: float, default 1e3
param tau_b:: Bridge detachment time (s), bounds (1e-6, 1e4)
type tau_b:: float, default 1.0
param tau_a:: Loop attachment time (s), bounds (1e-6, 1e4)
type tau_a:: float, default 5.0
param nu:: Bell force sensitivity (dimensionless), bounds (0.01, 20)
type nu:: float, default 1.0
param f_B_eq:: Equilibrium bridge fraction, bounds (0.01, 0.99)
type f_B_eq:: float, default 0.5
param eta_s:: Solvent viscosity (Pa·s), bounds (0.0, 1e4)
type eta_s:: float, default 0.0

parameters¶

Model parameters

Type:: ParameterSet

fitted_¶

Whether the model has been fitted

Type:: bool

Examples

Basic usage:

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

Startup with bridge fraction tracking:

>>> t = np.linspace(0, 100, 500)
>>> stress, f_B = model.simulate_startup(
...     t, gamma_dot=10.0, return_bridge_fraction=True
... )

See also

TNTSingleMode: Basic single-mode TNT (constant breakage)
TNTCates: Living polymer (wormlike micelle) model

__init__()[source]¶: Initialize TNT loop-bridge model.

initialize_from_creep(t, strain, sigma_applied=1.0)[source]¶

Initialize parameters from creep data for numerical stability.

Key insight: For creep simulation, eta_s must be non-zero to prevent infinite initial shear rate when S_xy starts at 0.

Parameters:

t (ndarray) – Time array (s)
strain (ndarray) – Strain array
sigma_applied (float) – Applied constant stress (Pa)

Return type:

initialize_from_relaxation(t, modulus)[source]¶

Initialize parameters from stress relaxation data.

Uses conservative tau_b estimate to ensure numerical stability with typical pre-shear rates (Wi = gamma_dot * tau_b should be < ~100).

Parameters:

t (ndarray) – Time array (s)
modulus (ndarray) – Relaxation modulus G(t) (Pa)

Return type:

property G: float¶: Get network modulus G (Pa).

property tau_b: float¶: Get bridge detachment time τ_b (s).

property tau_a: float¶: Get loop attachment time τ_a (s).

property nu: float¶: Get Bell force sensitivity ν (dimensionless).

property f_B_eq: float¶: Get equilibrium bridge fraction f_B_eq (dimensionless).

property eta_s: float¶: Get solvent viscosity η_s (Pa·s).

property G_eff: float¶

Get effective modulus G_eff = f_B_eq·G (Pa).

This is the linearized modulus at equilibrium.

property eta_0: float¶: Get zero-shear viscosity η₀ = f_B_eq·G·τ_b + η_s (Pa·s).

get_equilibrium_state()[source]¶

Return equilibrium state [f_B_eq, 1, 1, 1, 0].

At rest: f_B = f_B_eq, S = I (unstretched, isotropic)

Returns:: Equilibrium state [f_B, S_xx, S_yy, S_zz, S_xy]
Return type:: Array

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: [G, tau_b, tau_a, nu, f_B_eq, eta_s]
test_mode (str, optional) – Override stored test mode

Returns:

Predicted response

Return type:

jnp.ndarray

predict_flow_curve(gamma_dot, return_components=False)[source]¶

Predict steady shear stress and viscosity.

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, loop-bridge model reduces to effective Maxwell: G_eff = f_B_eq·G, τ_eff = τ_b

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:

predict_normal_stresses(gamma_dot)[source]¶

Predict first and second normal stress differences.

N₁ = f_B·G·(S_xx - S_yy) N₂ = 0 (upper-convected derivative)

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_bridge_fraction=False)[source]¶

Simulate startup flow at constant shear rate.

Parameters:

t (ndarray) – Time array (s)
gamma_dot (float) – Applied shear rate (1/s)
return_bridge_fraction (bool) – If True, also return f_B(t)

Returns:

Shear stress σ(t), or (σ, f_B) if return_bridge_fraction=True

Return type:

simulate_relaxation(t, gamma_dot_preshear, return_bridge_fraction=False)[source]¶

Simulate stress relaxation after cessation of steady shear.

Parameters:

t (ndarray) – Time array (s), starting from t=0 (cessation)
gamma_dot_preshear (float) – Shear rate before cessation (1/s)
return_bridge_fraction (bool) – If True, also return f_B(t)

Returns:

Relaxing stress σ(t), or (σ, f_B) if return_bridge_fraction=True

Return type:

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:

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’, ‘f_B’

Return type:

get_bridge_fraction_vs_rate(gamma_dot)[source]¶

Compute steady-state bridge fraction vs shear rate.

Parameters:: gamma_dot (ndarray) – Shear rate array (1/s)
Returns:: (gamma_dot, f_B_steady)
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:

TNTStickyRouse¶

rheojax.models.tnt.sticky_rouse.TNTStickyRouse | Handbook: TNT Sticky Rouse (Multi-Mode Sticker Dynamics) — Handbook Sticky Rouse model combining Rouse chain dynamics with reversible sticker groups.

Parameters (per mode `k`)¶
Parameter	Units	Bounds	Description
`G_k`	Pa	[1, 1e+08]	Mode k modulus
`tau_R_k`	s	[1e-06, 1e+04]	Mode k Rouse relaxation time

Global Parameters¶
Parameter	Units	Bounds	Description
`tau_s`	s	[1e-06, 1e+04]	Sticker lifetime
`eta_s`	Pa·s	[0, 1e+04]	Solvent viscosity

class rheojax.models.tnt.sticky_rouse.TNTStickyRouse(n_modes=3)[source]¶

Bases: TNTBase

Sticky Rouse model for associative polymers.

Multi-mode Maxwell model where sticker dynamics impose a relaxation time floor: τ_eff_k = max(τ_R_k, τ_s).

Creates a plateau in G(t) at intermediate times (sticker-dominated regime) before terminal relaxation (slowest Rouse mode).

Parameters:: n_modes (int) – Number of Rouse modes

parameters¶

Model parameters: - G_0, G_1, …, G_{N-1}: Mode moduli (Pa) - tau_R_0, tau_R_1, …, tau_R_{N-1}: Rouse relaxation times (s) - tau_s: Sticker lifetime (s) - eta_s: Solvent viscosity (Pa·s)

Type:: ParameterSet

Notes

The model reduces to standard multi-mode Maxwell when tau_s → 0. For tau_s → ∞, all modes relax at tau_s (single network behavior).

Examples

>>> # 3-mode sticky Rouse
>>> model = TNTStickyRouse(n_modes=3)
>>> model.fit(omega, G_star, test_mode='oscillation')
>>>
>>> # Predict plateau modulus
>>> G_plateau = model.predict_plateau_modulus()
>>>
>>> # Predict startup with stress overshoot
>>> t = np.linspace(0, 10, 200)
>>> sigma = model.predict(t, test_mode='startup', gamma_dot=1.0)
>>>
>>> # Extract effective relaxation times
>>> tau_eff = model.get_effective_times()

__init__(n_modes=3)[source]¶

Initialize Sticky Rouse model.

Parameters:: n_modes (int) – Number of Rouse modes (must be >= 1)

property n_modes: int¶: Number of Rouse modes.

property tau_s: float¶: Sticker lifetime (s).

property eta_s: float¶: Solvent viscosity (Pa·s).

get_effective_times()[source]¶

Get effective relaxation times for all modes.

Returns:: Effective times τ_eff_k = τ_R_k + τ_s, shape (N,)
Return type:: ndarray

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

Compute model prediction for given parameters.

Parameters:

X (Array) – Independent variable (time, frequency, or shear rate)
params (Array) – Parameter array [G_0, tau_R_0, G_1, tau_R_1, …, tau_s, eta_s] Length: 2*N + 2
test_mode (str | None) – Protocol: ‘oscillation’, ‘relaxation’, ‘flow_curve’, ‘startup’, ‘creep’, or ‘laos’

Returns:

Predicted response (protocol-dependent)

Return type:

predict_plateau_modulus()[source]¶

Compute plateau modulus G_N = Σ G_k for modes with τ_R_k < τ_s.

The plateau modulus is the sum of mode moduli for modes dominated by sticker lifetime (fast Rouse modes).

Returns:: Plateau modulus G_N (Pa)
Return type:: float

predict_zero_shear_viscosity()[source]¶

Compute zero-shear viscosity η₀ = Σ G_k·τ_eff_k + η_s.

Returns:: Zero-shear viscosity η₀ (Pa·s)
Return type:: float

predict_saos(omega, return_components=True)[source]¶

Predict SAOS storage and loss moduli.

Analytical superposition for multi-mode Maxwell:: G’(ω) = Σ G_k·(ωτ_eff_k)² / (1 + (ωτ_eff_k)²) G’’(ω) = Σ G_k·(ωτ_eff_k) / (1 + (ωτ_eff_k)²) + η_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:

predict_terminal_time()[source]¶

Return longest effective relaxation time (terminal mode).

Returns:: Terminal time τ_terminal = max(τ_eff_k) (s)
Return type:: float

predict_normal_stress_difference(gamma_dot)[source]¶

Predict first normal stress difference N₁(γ̇).

N₁ = Σ 2·G_k·τ_eff_k²·γ̇² / (1 + (τ_eff_k·γ̇)²)

Parameters:: gamma_dot (float | ndarray) – Shear rate (1/s)
Returns:: N₁ (Pa)
Return type:: ndarray

get_trajectory()[source]¶

Get stored ODE trajectory from last prediction.

Returns:: Dictionary with keys like ‘time’, ‘stress’, ‘strain’, ‘S_xy’
Return type:: dict[str, ndarray] | None

initialize_from_saos(omega, G_prime, G_double_prime)[source]¶

Initialize parameters from SAOS data.

Uses crossover frequency to estimate sticker lifetime and plateau modulus to distribute mode strengths.

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:

__repr__()[source]¶

Return string representation.

Return type:: str

TNTCates¶

rheojax.models.tnt.cates.TNTCates | Handbook: TNT Cates (Living Polymers / Wormlike Micelles) — Handbook Cates living polymer model with coupled reptation and reversible scission for worm-like micelle solutions.

Parameters¶
Parameter	Units	Bounds	Description
`G_0`	Pa	[1, 1e+08]	Plateau modulus
`tau_rep`	s	[1e-04, 1e+06]	Reptation time (curvilinear diffusion)
`tau_break`	s	[1e-06, 1e+04]	Average breaking time (scission events)
`eta_s`	Pa·s	[0, 1e+04]	Solvent viscosity

class rheojax.models.tnt.cates.TNTCates[source]¶

Bases: TNTBase

Cates living polymer (wormlike micelle) model.

Implements the Cates theory for living polymers with reversible scission. In the fast-breaking limit, the system behaves as a single Maxwell mode with effective relaxation time τ_d = √(τ_rep · τ_break).

The constitutive equation is identical to the basic TNT model (constant breakage, linear stress, upper-convected derivative), but with τ_d replacing the single bond lifetime τ_b:

dS/dt = L·S + S·L^T + (1/τ_d)·I - (1/τ_d)·S σ = G₀·S_xy + η_s·γ̇

Parameters:

G_0 (float, default 1e3) – Plateau modulus (Pa). Network elastic modulus.
tau_rep (float, default 10.0) – Reptation time (s). Curvilinear diffusion time along the tube.
tau_break (float, default 0.1) – Average breaking time (s). Mean time between scission events.
eta_s (float, default 0.0) – Solvent viscosity (Pa·s). Newtonian background contribution.
Derived
-------
tau_d (float) – Effective relaxation time τ_d = √(τ_rep · τ_break)
eta_0 (float) – Zero-shear viscosity η₀ = G₀·τ_d + η_s

parameters¶

Model parameters for fitting

Type:: ParameterSet

fitted_¶

Whether the model has been fitted

Type:: bool

Examples

Basic usage with default parameters:

>>> model = TNTCates()
>>> print(model.tau_d)  # Effective time
1.0

Fit to SAOS data:

>>> omega = np.logspace(-2, 2, 50)
>>> G_star = generate_synthetic_data(omega)
>>> model.fit(omega, G_star, test_mode='oscillation')

Predict flow curve:

>>> gamma_dot = np.logspace(-2, 2, 50)
>>> sigma = model.predict_flow_curve(gamma_dot)

See also

TNTSingleMode: Single-mode TNT with variants
TNTLoopBridge: Two-species loop-bridge kinetics

__init__()[source]¶: Initialize Cates living polymer model.

property G_0: float¶: Get plateau modulus G₀ (Pa).

property tau_rep: float¶: Get reptation time τ_rep (s).

property tau_break: float¶: Get breaking time τ_break (s).

property eta_s: float¶: Get solvent viscosity η_s (Pa·s).

property tau_d: float¶: Get effective relaxation time τ_d = √(τ_rep · τ_break) (s).

property eta_0: float¶: Get zero-shear viscosity η₀ = G₀·τ_d + η_s (Pa·s).

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: [G_0, tau_rep, tau_break, eta_s]
test_mode (str, optional) – Override stored test mode

Returns:

Predicted response

Return type:

jnp.ndarray

predict_flow_curve(gamma_dot, return_components=False)[source]¶

Predict steady shear stress and viscosity.

For Cates model with constant breakage: σ = G₀·τ_d·γ̇ + η_s·γ̇ (UCM-like, no 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.

Cates model reduces to single-mode Maxwell with τ_d: G’(ω) = G₀·(ωτ_d)²/(1+(ωτ_d)²) G’’(ω) = G₀·(ωτ_d)/(1+(ωτ_d)²) + η_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:

predict_normal_stresses(gamma_dot)[source]¶

Predict first and second normal stress differences.

Cates model (UCM-like): N₁ = 2G₀·(τ_d·γ̇)² N₂ = 0

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)[source]¶

Simulate stress relaxation after cessation of steady shear.

Analytical single-exponential decay: σ(t) = G₀·τ_d·γ̇·exp(-t/τ_d)

Parameters:

t (ndarray) – Time array (s), starting from t=0 (cessation)
gamma_dot_preshear (float) – Shear rate before cessation (1/s)

Returns:

Relaxing stress σ(t)

Return type:

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:

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:

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

Get relaxation modulus G(t).

For Cates model: G(t) = G₀·exp(-t/τ_d)

Parameters:

t (ndarray | None) – Time array (default: logspace from 0.01·τ_d to 100·τ_d)
n_points (int) – Number of points if t not provided

Returns:

(t, G(t))

Return type:

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:

__repr__()[source]¶

Return string representation.

Return type:: str

TNTMultiSpecies¶

rheojax.models.tnt.multi_species.TNTMultiSpecies | Handbook: TNT Multi-Species (Multiple Bond Types) — Handbook Multi-species transient network model with multiple independent bond types of distinct lifetimes and moduli.

Parameters (per species `i`)¶
Parameter	Units	Bounds	Description
`G_i`	Pa	[1, 1e+08]	Network modulus for bond species i
`tau_b_i`	s	[1e-06, 1e+04]	Bond lifetime for species i

Global Parameters¶
Parameter	Units	Bounds	Description
`eta_s`	Pa·s	[0, 1e+04]	Solvent viscosity

class rheojax.models.tnt.multi_species.TNTMultiSpecies(n_species=2)[source]¶

Bases: TNTBase

Multi-species Transient Network Theory model.

Implements a network with N independent bond populations, each with its own modulus G_i and lifetime τ_b_i. The total stress is a superposition of N Maxwell modes.

This is equivalent to a generalized Maxwell model where each mode represents a distinct physical crosslink species rather than a mathematical decomposition.

Parameters:: n_species (int) – Number of distinct bond species (N ≥ 1)

parameters¶

Model parameters: [G_0, tau_b_0, G_1, tau_b_1, …, G_{N-1}, tau_b_{N-1}, eta_s]

Type:: ParameterSet

fitted_¶

Whether the model has been fitted

Type:: bool

_n_species¶

Number of species

Type:: int

Notes

The state vector has 4*N components:: [S_xx_0, S_yy_0, S_zz_0, S_xy_0, …, S_xy_{N-1}]

Each species evolves independently via the upper-convected derivative with constant breakage/creation rates.

See also

TNTSingleMode: Single-mode TNT with variant breakage rates
GeneralizedMaxwell: Mathematical multi-mode decomposition

__init__(n_species=2)[source]¶

Initialize multi-species TNT model.

Parameters:: n_species (int) – Number of distinct bond species (must be ≥ 1)

initialize_from_creep(t, strain, sigma_applied=1.0)[source]¶

Initialize parameters from creep data for numerical stability.

For multi-species, sets eta_s to prevent infinite initial shear rate.

Parameters:

t (ndarray) – Time array (s)
strain (ndarray) – Strain array
sigma_applied (float) – Applied constant stress (Pa)

Return type:

initialize_from_relaxation(t, modulus)[source]¶

Initialize parameters from stress relaxation data.

Parameters:

t (ndarray) – Time array (s)
modulus (ndarray) – Relaxation modulus G(t) (Pa)

Return type:

property n_species: int¶: Get number of bond species N.

property eta_s: float¶: Get solvent viscosity η_s (Pa·s).

property G_total: float¶: Get total modulus G_total = Σ G_i (Pa).

property eta_0: float¶: Get zero-shear viscosity η₀ = Σ G_i·τ_b_i + η_s (Pa·s).

get_equilibrium_conformation_multimode()[source]¶

Return equilibrium conformation for all N modes.

Returns:: Equilibrium state [1, 1, 1, 0, …, 1, 1, 1, 0] with shape (4N,)
Return type:: Array

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: [G_0, tau_b_0, G_1, tau_b_1, …, G_{N-1}, tau_b_{N-1}, eta_s] Total length: 2*N + 1
test_mode (str, optional) – Override stored test mode

Returns:

Predicted response

Return type:

jnp.ndarray

predict_flow_curve(gamma_dot, return_components=False)[source]¶

Predict steady shear stress and viscosity.

Analytical superposition: σ = Σ G_i·τ_b_i·γ̇ / (1 + (τ_b_i·γ̇)²) + η_s·γ̇

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.

Analytical superposition:: G’(ω) = Σ G_i·(ωτ_i)² / (1 + (ωτ_i)²) G’’(ω) = Σ G_i·(ωτ_i) / (1 + (ωτ_i)²) + η_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:

predict_normal_stresses(gamma_dot)[source]¶

Predict first and second normal stress differences.

For multi-mode UCM (conformation tensor):: N₁ = Σ 2·G_i·(τ_b_i·γ̇)² N₂ = 0 (upper-convected derivative)

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 for all modes

Returns:

Shear stress σ(t), or dict with ‘S_xx’, ‘S_yy’, ‘S_xy’, ‘S_zz’ (each shape (T, N)) if return_full=True

Return type:

simulate_relaxation(t, gamma_dot_preshear)[source]¶

Simulate stress relaxation after cessation of steady shear.

Analytical multi-mode relaxation:: σ(t) = Σ σ₀_i·exp(-t/τ_b_i)

Parameters:

t (ndarray) – Time array (s), starting from t=0 (cessation)
gamma_dot_preshear (float) – Shear rate before cessation (1/s)

Returns:

Relaxing stress σ(t) (Pa)

Return type:

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:

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:

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

Get relaxation modulus G(t).

For multi-species TNT: G(t) = Σ G_i·exp(-t/τ_b_i)

Parameters:

t (ndarray | None) – Time array (default: logspace from 0.01·min(τ) to 100·max(τ))
n_points (int) – Number of points if t not provided

Returns:

(t, G(t))

Return type:

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:

__repr__()[source]¶

Return string representation.

Return type:: str

VLB Models (Vernerey-Long-Brighenti)¶

VLBLocal¶

rheojax.models.vlb.local.VLBLocal | Handbook: VLB Transient Network Models Local VLB transient network model with chain-length-resolved bond kinetics and non-affine deformation.

Parameters¶
Parameter	Units	Bounds	Description
`G0`	Pa	[1, 1e+08]	Network modulus (\(n k_B T\))
`k_d`	1/s	[1e-06, 1e+06]	Dissociation rate (inverse relaxation time)

class rheojax.models.vlb.local.VLBLocal[source]¶

Bases: VLBBase

Single transient network VLB model (2 params: G0, k_d).

Implements the VLB framework for a single population of dynamic crosslinks with constant dissociation rate. Analytically equivalent to the Maxwell model but with molecular-statistical foundations.

The distribution tensor mu has equilibrium mu = I and evolves via:

dmu/dt = k_d*(I - mu) + L·mu + mu·L^T

Cauchy stress: sigma = G0*(mu - I)

Parameters:

G0 (float) – Network modulus (Pa), physically G0 = c*k_B*T where c is chain density
k_d (float) – Dissociation rate (1/s), inverse of relaxation time t_R = 1/k_d

parameters¶

Model parameters (G0, k_d)

Type:: ParameterSet

fitted_¶

Whether the model has been fitted

Type:: bool

See also

VLBMultiNetwork: Multi-network VLB with N transient + permanent + solvent

__init__()[source]¶: Initialize single-network VLB model.

property G0: float¶: Get network modulus G0 (Pa).

property k_d: float¶: Get dissociation rate k_d (1/s).

property relaxation_time: float¶: Get relaxation time t_R = 1/k_d (s).

property viscosity: float¶: Get zero-shear viscosity eta_0 = G0/k_d (Pa*s).

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: [G0, k_d]
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 optionally viscosity, N1).

Newtonian: sigma = G0*gamma_dot/k_d.

Parameters:

gamma_dot (ndarray) – Shear rate array (1/s)
return_components (bool) – If True, return (sigma, eta, N1)

Returns:

Stress, or (stress, viscosity, N1) if return_components=True

Return type:

ndarray | tuple[ndarray, ndarray, ndarray]

predict_saos(omega, return_components=True)[source]¶

Predict SAOS storage and loss moduli.

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:

predict_normal_stresses(gamma_dot)[source]¶

Predict steady-state first and second normal stress differences.

N1 = 2*G0*(gamma_dot/k_d)^2, N2 = 0 (upper-convected Maxwell).

Parameters:: gamma_dot (ndarray) – Shear rate array (1/s)
Returns:: (N1, N2) in Pa. N2 is always zero for upper-convected VLB.
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 (sigma, N1, strain)

Returns:

Stress, or (stress, N1, strain) if return_full=True

Return type:

ndarray | tuple[ndarray, ndarray, ndarray]

simulate_relaxation(t)[source]¶

Simulate stress relaxation G(t) = G0*exp(-k_d*t).

Parameters:: t (ndarray) – Time after cessation of flow (s)
Returns:: Relaxation modulus G(t) (Pa)
Return type:: ndarray

simulate_creep(t, sigma_0, return_full=False)[source]¶

Simulate creep under constant applied stress.

Parameters:

t (ndarray) – Time array (s)
sigma_0 (float) – Applied stress (Pa)
return_full (bool) – If True, return (strain, compliance)

Returns:

Strain gamma(t), or (gamma, J) if return_full=True

Return type:

simulate_laos(t, gamma_0, omega)[source]¶

Simulate large-amplitude oscillatory shear (LAOS).

Parameters:

t (ndarray) – Time array (s), should span at least 3-5 full cycles
gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)

Returns:

Keys: ‘time’, ‘strain’, ‘stress’, ‘N1’, ‘gamma_dot’

Return type:

predict_uniaxial_extension(epsilon_dot, return_trouton=False)[source]¶

Predict steady-state uniaxial extensional stress.

Parameters:

epsilon_dot (ndarray) – Extensional strain rate (1/s)
return_trouton (bool) – If True, also return Trouton ratio

Returns:

Extensional stress, or (stress, Trouton_ratio)

Return type:

simulate_uniaxial_extension(t, epsilon_dot)[source]¶

Simulate transient uniaxial extension.

Parameters:

t (ndarray) – Time array (s)
epsilon_dot (float) – Extensional strain rate (1/s)

Returns:

(extensional_stress, extensional_viscosity) as functions of time

Return type:

get_relaxation_spectrum()[source]¶

Get relaxation spectrum as list of (G, tau) pairs.

Returns:: [(G0, 1/k_d)] — single mode
Return type:: list[tuple[float, float]]

extract_laos_harmonics(laos_result, n_harmonics=5)[source]¶

Extract Fourier harmonics from LAOS results.

Parameters:

laos_result (dict) – Output from simulate_laos()
n_harmonics (int) – Number of harmonics to extract

Returns:

Keys: ‘harmonic_index’, ‘sigma_harmonics’, ‘N1_harmonics’

Return type:

VLBMultiNetwork¶

rheojax.models.vlb.multi_network.VLBMultiNetwork | Handbook: VLB Transient Network Models Multi-network VLB model with multiple independent transient sub-networks of distinct moduli and kinetics.

Parameters (per network `i`)¶
Parameter	Units	Bounds	Description
`G0_i`	Pa	[1, 1e+08]	Network i modulus
`k_d_i`	1/s	[1e-06, 1e+06]	Network i dissociation rate

class rheojax.models.vlb.multi_network.VLBMultiNetwork(n_modes=2, include_permanent=False)[source]¶

Bases: VLBBase

Multi-network VLB model: M transient + optional permanent + solvent.

Implements a network with N independent transient crosslink populations, each with modulus G_i and dissociation rate k_d_i. The total stress is a superposition of N Maxwell modes.

Parameters:

n_modes (int) – Number of distinct transient network populations (N >= 1)
include_permanent (bool) – Whether to include a permanent (elastic) network (G_e)

parameters¶

Model parameters: [G_0, k_d_0, G_1, k_d_1, …, eta_s, (G_e)]

Type:: ParameterSet

fitted_¶

Whether the model has been fitted

Type:: bool

_n_modes¶

Number of transient modes

Type:: int

Notes

Parameter ordering: [G_0, k_d_0, G_1, k_d_1, …, G_{N-1}, k_d_{N-1}, eta_s, (G_e)] Total parameter count: 2N + 1 (without permanent) or 2N + 2 (with permanent)

See also

VLBLocal: Single transient network (2 parameters)

__init__(n_modes=2, include_permanent=False)[source]¶

Initialize multi-network VLB model.

Parameters:

n_modes (int) – Number of transient network populations (must be >= 1)
include_permanent (bool) – Include permanent elastic network

property n_modes: int¶: Number of transient network modes.

property include_permanent: bool¶: Whether a permanent network is included.

property G_e: float¶: Permanent network modulus (Pa). 0 if not included.

property eta_s: float¶: Solvent viscosity (Pa*s).

property G_total: float¶

sum G_i + G_e.

Type:: Total modulus

property eta_0: float¶

sum G_i/k_d_i + eta_s.

Type:: Zero-shear viscosity

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

NumPyro/BayesianMixin model function.

Routes to appropriate prediction based on test_mode.

Parameters:

X (array-like) – Independent variable
params (array-like) – Parameter values: [G_0, k_d_0, …, G_{N-1}, k_d_{N-1}, eta_s, (G_e)]
test_mode (str, optional) – Override stored test mode
**kwargs – Protocol-specific parameters

Returns:

Predicted response

Return type:

jnp.ndarray

predict_saos(omega, return_components=True)[source]¶

Predict multi-mode SAOS moduli.

Parameters:

omega (ndarray) – Angular frequency (rad/s)
return_components (bool) – If True, return (G’, G’’)

Returns:

(G’, G’’) or |G*|

Return type:

get_relaxation_spectrum()[source]¶

Get relaxation spectrum as list of (G, tau) pairs.

Returns:: [(G_i, 1/k_d_i)] for each transient mode
Return type:: list[tuple[float, float]]

__repr__()[source]¶

Return string representation.

Return type:: str

VLBVariant¶

rheojax.models.vlb.variant.VLBVariant | Handbook: VLB Advanced Theory & Numerical Methods VLB model variant with modified kinetics (Bell, FENE-P, stretch-creation) for enhanced nonlinear rheology.

class rheojax.models.vlb.variant.VLBVariant(breakage='constant', stress_type='linear', temperature=False)[source]¶

Bases: VLBBase

VLB with Bell breakage, FENE-P stress, and/or temperature dependence.

This is the composable variant class for VLB models. It supports all 6 protocols via ODE integration (required when k_d depends on state).

When breakage=”constant” and stress_type=”linear”, the model matches VLBLocal exactly (regression verified).

Parameters:

breakage (Literal['constant', 'bell']) – Breakage rate function: “constant” or “bell”
stress_type (Literal['linear', 'fene']) – Stress formula: “linear” or “fene”
temperature (bool) – If True, enable Arrhenius temperature dependence

__init__(breakage='constant', stress_type='linear', temperature=False)[source]¶

Initialize VLBVariant model.

property G0: float¶: Network modulus (Pa).

property k_d_0: float¶: Unstressed dissociation rate (1/s).

property nu: float | None¶: Force sensitivity parameter (Bell only).

property L_max: float | None¶: Maximum extensibility (FENE only).

property relaxation_time: float¶: Equilibrium relaxation time t_R = 1/k_d_0 (s).

property viscosity: float¶: Zero-shear viscosity eta_0 = G0/k_d_0 (Pa·s).

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

NumPyro/BayesianMixin model function.

Routes to appropriate prediction based on test_mode.

Parameters:

X (array-like) – Independent variable
params (array-like) – Parameter values in ParameterSet order: [G0, k_d_0, eta_s, (nu), (L_max), (E_a, T_ref)]
test_mode (str, optional) – Override stored test mode
**kwargs – Protocol-specific: gamma_dot, sigma_applied, gamma_0, omega, T

Returns:

Predicted response

Return type:

jnp.ndarray

predict_flow_curve(gamma_dot)[source]¶

Predict steady-state flow curve.

Parameters:

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

Return type:

Returns:

sigma (np.ndarray) – Steady-state shear stress (Pa)
eta (np.ndarray) – Apparent viscosity (Pa·s)

predict_saos(omega)[source]¶

Predict SAOS moduli (Maxwell, analytical).

In the linear regime, Bell reduces to constant k_d = k_d_0.

Parameters:

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

Return type:

Returns:

G_prime (np.ndarray) – Storage modulus G’ (Pa)
G_double_prime (np.ndarray) – Loss modulus G’’ (Pa)

predict_normal_stresses(gamma_dot)[source]¶

Predict steady-state first normal stress difference N1.

For Bell breakage, this requires ODE integration.

Parameters:: gamma_dot (ndarray) – Shear rate array (1/s)
Returns:: N1 values (Pa)
Return type:: ndarray

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

Simulate startup shear.

Parameters:

t (ndarray) – Time array (s)
gamma_dot (float) – Applied shear rate (1/s)
return_full (bool) – If True, return dict with stress, N1, strain

Returns:

Shear stress sigma(t), or dict with full trajectory

Return type:

ndarray | dict

simulate_relaxation(t, gamma_dot_preshear=10.0)[source]¶

Simulate stress relaxation after cessation of flow.

Parameters:

t (ndarray) – Time array (s)
gamma_dot_preshear (float) – Pre-shear rate (1/s)

Returns:

Relaxing stress sigma(t)

Return type:

simulate_creep(t, sigma_applied)[source]¶

Simulate creep (stress-controlled).

Parameters:

t (ndarray) – Time array (s)
sigma_applied (float) – Applied stress (Pa)

Returns:

Strain gamma(t)

Return type:

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

Simulate Large Amplitude Oscillatory Shear (LAOS).

Parameters:

t (ndarray | None) – Time array (if None, auto-generated from n_cycles)
gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)
n_cycles (int) – Number of cycles (if t is None)

Returns:

‘t’, ‘strain’, ‘stress’, ‘gamma_dot’

Return type:

predict_uniaxial_extension(eps_dot)[source]¶

Predict steady-state extensional stress.

For FENE-P, extensional stress is bounded (no singularity).

Parameters:: eps_dot (ndarray) – Extension rate array (1/s)
Returns:: Extensional stress sigma_E (Pa)
Return type:: ndarray

extract_laos_harmonics(result=None, n_harmonics=5)[source]¶

Extract Fourier harmonics from LAOS data.

Parameters:

result (dict | None) – LAOS result dict (uses stored _trajectory if None)
n_harmonics (int) – Number of harmonics to extract

Returns:

‘harmonics’: array of harmonic intensities ‘I3_I1’: third-to-first harmonic ratio

Return type:

weissenberg_number(gamma_dot)[source]¶

Compute Weissenberg number Wi = t_R * gamma_dot.

Return type:: float

deborah_number(omega)[source]¶

Compute Deborah number De = t_R * omega.

Return type:: float

VLBNonlocal¶

rheojax.models.vlb.nonlocal_model.VLBNonlocal | Handbook: VLB Advanced Theory & Numerical Methods Nonlocal VLB model with spatial diffusion for shear banding and inhomogeneous deformation prediction.

class rheojax.models.vlb.nonlocal_model.VLBNonlocal(breakage='constant', stress_type='linear', n_points=51, gap_width=0.001)[source]¶

Bases: VLBBase

Nonlocal VLB with tensor diffusion for shear banding.

Solves a 1D PDE across the gap of a Couette cell. The state at each spatial point is (mu_xx, mu_yy, mu_zz, mu_xy), plus a single wall stress Sigma (spatially uniform at low Reynolds number).

Shear banding occurs when the Bell breakage rate creates a non-monotonic flow curve. The diffusion term D_mu * nabla^2(mu) regularizes the interface with width ~ xi = sqrt(D_mu / k_d_0).

Parameters:

breakage (Literal['constant', 'bell']) – “constant” or “bell”
stress_type (Literal['linear', 'fene']) – “linear” or “fene”
n_points (int) – Spatial grid resolution
gap_width (float) – Gap width (m)

__init__(breakage='constant', stress_type='linear', n_points=51, gap_width=0.001)[source]¶

Initialize VLBNonlocal model.

property G0: float¶

property k_d_0: float¶

get_cooperativity_length()[source]¶

Cooperativity length xi = sqrt(D_mu / k_d_0).

This sets the shear band interface width.

Returns:: Cooperativity length (m)
Return type:: float

simulate_steady_shear(gamma_dot_avg, t_end=100.0, dt=0.1, perturbation=0.01)[source]¶

Simulate approach to steady state under imposed average shear rate.

Parameters:

gamma_dot_avg (float) – Imposed average shear rate (1/s)
t_end (float) – Simulation end time (s)
dt (float) – Output time step (s)
perturbation (float) – Initial spatial noise amplitude

Returns:

‘t’: time array ‘y’: spatial grid ‘mu_xy’: mu_xy profiles (N_t, N_y) ‘gamma_dot’: shear rate profiles (N_t, N_y) ‘stress’: wall stress Sigma(t)

Return type:

simulate_startup(gamma_dot_avg, t_end=50.0, dt=0.05)[source]¶

Simulate startup from rest with banding evolution.

Parameters:

gamma_dot_avg (float) – Imposed average shear rate (1/s)
t_end (float) – End time (s)
dt (float) – Output time step (s)

Returns:

Same format as simulate_steady_shear

Return type:

simulate_creep(sigma_0, t_end=100.0, dt=0.1)[source]¶

Simulate stress-controlled creep with spatial resolution.

In creep, the stress Sigma is held fixed (no feedback).

Parameters:

sigma_0 (float) – Applied stress (Pa)
t_end (float) – End time (s)
dt (float) – Output time step (s)

Returns:

‘t’, ‘y’, ‘gamma_dot’, ‘mu_xy’, ‘velocity’

Return type:

detect_banding(result, threshold=0.1)[source]¶

Detect shear banding from steady-state profiles.

Parameters:

result (dict) – Result from simulate_steady_shear()
threshold (float) – Relative variation threshold for banding detection

Returns:

‘is_banding’: bool ‘band_contrast’: max/min shear rate ratio ‘band_width’: approximate band width (m) ‘band_location’: center of high-shear band (m)

Return type:

get_velocity_profile(result)[source]¶

Compute velocity profile from final shear rate profile.

v(y) = integral_0^y gamma_dot(y’) dy’

Parameters:: result (dict) – Result from simulate_steady_shear()
Returns:: Velocity profile v(y)
Return type:: ndarray

plot_profiles(result, ax=None)[source]¶

Plot spatial profiles (shear rate and mu_xy).

Parameters:

result (dict) – Result from simulate_steady_shear()
ax (matplotlib axes, optional) – If None, creates new figure

Return type:

matplotlib figure

HVM Model (Hybrid Vitrimer)¶

HVMLocal¶

rheojax.models.hvm.local.HVMLocal | Handbook: HVM Advanced Theory & Numerical Methods Hybrid vitrimer constitutive model with permanent covalent crosslinks, associative exchangeable bonds (BER/TST kinetics), and optional dissociative physical bonds.

Parameters¶
Parameter	Units	Bounds	Description
`G_P`	Pa	[1, 1e+08]	Permanent network modulus (covalent crosslinks)
`G_E`	Pa	[1, 1e+08]	Exchangeable network modulus (vitrimer bonds)
`G_D`	Pa	[0, 1e+08]	Dissociative network modulus (physical bonds, optional)
`nu_0`	1/s	[1e+06, 1e+14]	BER attempt frequency (Eyring pre-factor)
`E_a`	J/mol	[4e+04, 2.5e+05]	BER activation energy
`V_act`	m³	[1e-08, 0.01]	BER activation volume
`k_d_D`	1/s	[1e-06, 1e+06]	Dissociation rate for D-network (optional)
`T`	K	[250, 350]	Temperature

class rheojax.models.hvm.local.HVMLocal(kinetics='stress', include_damage=False, include_dissociative=True)[source]¶

Bases: HVMBase

Local (0D) Hybrid Vitrimer Model.

A constitutive model for vitrimers combining: - Permanent network (P): covalent crosslinks, elastic - Exchangeable network (E): vitrimer bonds with TST kinetics - Dissociative network (D): physical bonds, Maxwell relaxation

Parameters:

kinetics (Literal['stress', 'stretch']) – TST coupling mechanism for bond exchange rate
include_damage (bool) – Enable cooperative damage shielding
include_dissociative (bool) – Include dissociative (D) network

Examples

>>> from rheojax.models import HVMLocal
>>> model = HVMLocal()
>>> omega = np.logspace(-2, 2, 50)
>>> G_prime, G_double_prime = model.predict_saos(omega)

>>> # Partial vitrimer (Meng 2019)
>>> model = HVMLocal(include_dissociative=False)

>>> # With TST stress feedback
>>> model = HVMLocal(kinetics='stress')
>>> t = np.linspace(0, 10, 200)
>>> result = model.simulate_startup(t, gamma_dot=1.0, return_full=True)

__init__(kinetics='stress', include_damage=False, include_dissociative=True)[source]¶

Initialize VLB base model.

predict_flow_curve(gamma_dot, return_components=False)[source]¶

Predict steady-state flow curve.

At steady state, mu^E → mu^E_nat, so sigma_E → 0. Only the D-network contributes viscous stress.

Parameters:

gamma_dot (ndarray) – Shear rate array (1/s)
return_components (bool) – If True, return dict with subnetwork contributions

Returns:

Steady-state stress (Pa) or component dict

Return type:

predict_saos(omega, return_components=True)[source]¶

Predict SAOS storage and loss moduli.

Two Maxwell modes (E, D) plus permanent plateau (P).

Parameters:

omega (ndarray) – Angular frequency array (rad/s)
return_components (bool) – If True, return (G’, G’’); if False, return |G*|

Returns:

(G’, G’’) or |G*|

Return type:

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

Simulate startup shear flow.

Uses analytical solution for constant-rate case, or ODE with TST feedback for nonlinear regime.

Parameters:

t (ndarray) – Time array (s)
gamma_dot (float) – Applied shear rate (1/s)
return_full (bool) – If True, return dict with all trajectories

Returns:

Stress array or full trajectory dict

Return type:

simulate_relaxation(t, gamma_step=1.0, return_full=False)[source]¶

Simulate stress relaxation after step strain.

Parameters:

t (ndarray) – Time array after step (s)
gamma_step (float) – Applied step strain
return_full (bool) – If True, return full trajectory dict

Returns:

G(t) relaxation modulus or trajectory dict

Return type:

simulate_creep(t, sigma_0, return_full=False)[source]¶

Simulate creep under constant stress.

Parameters:

t (ndarray) – Time array (s)
sigma_0 (float) – Applied constant stress (Pa)
return_full (bool) – If True, return full trajectory dict

Returns:

Strain gamma(t) or trajectory dict

Return type:

simulate_laos(t, gamma_0, omega)[source]¶

Simulate LAOS (Large Amplitude Oscillatory Shear).

Parameters:

t (ndarray) – Time array (s)
gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)

Returns:

Keys: time, strain, stress, gamma_dot, N1, mu_E_xy, mu_E_nat_xy, mu_D_xy, damage

Return type:

predict_normal_stresses(gamma_dot)[source]¶

Predict steady-state normal stress differences.

At steady state, E-network contributes zero normal stress (mu^E → mu^E_nat). Only D-network contributes N1.

N1 = 2 * G_D * (gamma_dot / k_d_D)^2 N2 = 0 (upper-convected Maxwell)

Parameters:: gamma_dot (ndarray) – Shear rate array (1/s)
Returns:: (N1, N2) arrays (Pa)
Return type:: tuple[ndarray, ndarray]

extract_laos_harmonics(laos_result, n_harmonics=5)[source]¶

Extract Fourier harmonics from LAOS simulation.

Parameters:

laos_result (dict[str, ndarray]) – Output from simulate_laos()
n_harmonics (int) – Number of harmonics to extract

Returns:

Keys: harmonic_index (1, 3, 5, …), sigma_harmonics, N1_harmonics

Return type:

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

NumPyro/BayesianMixin model function for HVM.

Routes to appropriate JAX-traceable prediction based on test_mode. Required by BayesianMixin for NumPyro NUTS sampling.

Parameters:

X (array-like) – Independent variable
params (array-like) – Parameter values in ParameterSet order
test_mode (str, optional) – Override stored test mode
**kwargs – Protocol-specific: gamma_dot, sigma_applied, gamma_0, omega

Returns:

Predicted response

Return type:

jnp.ndarray

classmethod neo_hookean(G_P=10000.0)[source]¶

Create neo-Hookean elastic solid (G_E=0, G_D=0).

Parameters:: G_P (float) – Permanent network modulus (Pa)
Returns:: Model with only P-network active
Return type:: HVMLocal

classmethod maxwell(G_D=10000.0, k_d_D=1.0)[source]¶

Create single Maxwell fluid (G_P=0, G_E=0).

Parameters:

G_D (float) – Network modulus (Pa)
k_d_D (float) – Dissociation rate (1/s)

Returns:

Model with only D-network active

Return type:

classmethod zener(G_P=10000.0, G_D=10000.0, k_d_D=1.0)[source]¶

Create Zener/SLS model (G_E=0).

Parameters:

G_P (float) – Permanent network modulus (Pa)
G_D (float) – Dissociative network modulus (Pa)
k_d_D (float) – Dissociation rate (1/s)

Returns:

Model with P + D networks (no vitrimer exchange)

Return type:

classmethod pure_vitrimer(G_E=10000.0, nu_0=10000000000.0, E_a=80000.0, V_act=1e-05, T=300.0)[source]¶

Create pure vitrimer (G_P=0, G_D=0).

Parameters:

G_E (float) – Exchangeable network modulus (Pa)
nu_0 (float) – TST parameters
E_a (float) – TST parameters
V_act (float) – TST parameters
T (float) – TST parameters

Returns:

Model with only E-network active

Return type:

classmethod partial_vitrimer(G_P=10000.0, G_E=10000.0, nu_0=10000000000.0, E_a=80000.0, V_act=1e-05, T=300.0)[source]¶

Create partial vitrimer (G_D=0, Meng 2019).

Parameters:

G_P (float) – Permanent network modulus (Pa)
G_E (float) – Exchangeable network modulus (Pa)
nu_0 (float) – TST parameters
E_a (float) – TST parameters
V_act (float) – TST parameters
T (float) – TST parameters

Returns:

Model with P + E networks (no dissociative bonds)

Return type:

HVNM Model (Hybrid Vitrimer Nanocomposite)¶

HVNMLocal¶

rheojax.models.hvnm.local.HVNMLocal | Handbook: HVNM Advanced Theory & Numerical Methods Hybrid vitrimer nanocomposite model extending HVM with a 4th interphase sub-network around nanoparticles. Guth-Gold strain amplification \(X(\phi)\).

Parameters¶
Parameter	Units	Bounds	Description
`G_P`	Pa	[1, 1e+08]	Permanent network modulus (amplified by \(X(\phi)\))
`G_E`	Pa	[1, 1e+08]	Exchangeable network modulus (matrix TST kinetics)
`G_D`	Pa	[0, 1e+08]	Dissociative network modulus (optional)
`G_I`	Pa	[1, 1e+08]	Interphase network modulus (optional)
`nu_0`	1/s	[1e+06, 1e+14]	Matrix BER attempt frequency
`E_a`	J/mol	[4e+04, 2.5e+05]	Matrix BER activation energy
`nu_0_int`	1/s	[1e+06, 1e+14]	Interfacial BER attempt frequency (optional)
`E_a_int`	J/mol	[4e+04, 2.5e+05]	Interfacial BER activation energy (optional)
`phi`	–	[0, 0.4]	NP volume fraction
`R_NP`	m	[1e-09, 1e-06]	NP radius
`beta_I`	–	[1, 10]	Interphase reinforcement factor
`T`	K	[250, 350]	Temperature

class rheojax.models.hvnm.local.HVNMLocal(kinetics='stress', include_damage=False, include_dissociative=True, include_interfacial_damage=False, include_diffusion=False)[source]¶

Bases: HVNMBase

Local (0D) Hybrid Vitrimer Nanocomposite Model.

A constitutive model for NP-filled vitrimers combining: - Permanent network (P): covalent crosslinks, amplified by X(phi) - Exchangeable network (E): vitrimer bonds with matrix TST kinetics - Dissociative network (D): physical bonds, Maxwell relaxation - Interphase network (I): NP-bound chains with interfacial TST kinetics

Parameters:

kinetics (Literal['stress', 'stretch']) – TST coupling mechanism for bond exchange rates
include_damage (bool) – Enable matrix cooperative damage shielding
include_dissociative (bool) – Include dissociative (D) network
include_interfacial_damage (bool) – Enable interfacial damage with self-healing
include_diffusion (bool) – Enable diffusion-limited relaxation tails

Examples

>>> from rheojax.models import HVNMLocal
>>> model = HVNMLocal()
>>> model.parameters.set_value("phi", 0.1)
>>> omega = np.logspace(-2, 2, 50)
>>> G_prime, G_double_prime = model.predict_saos(omega)

>>> # Unfilled limit (recovers HVM)
>>> model = HVNMLocal()
>>> model.parameters.set_value("phi", 0.0)

__init__(kinetics='stress', include_damage=False, include_dissociative=True, include_interfacial_damage=False, include_diffusion=False)[source]¶

Initialize VLB base model.

predict_flow_curve(gamma_dot, return_components=False)[source]¶

Predict steady-state flow curve.

At steady state, mu^E -> mu^E_nat and mu^I -> mu^I_nat, so sigma_E -> 0 and sigma_I -> 0. Only the D-network contributes viscous stress: sigma_D = eta_D * gamma_dot.

Parameters:

gamma_dot (ndarray) – Shear rate array (1/s)
return_components (bool) – If True, return dict with subnetwork contributions

Returns:

Steady-state stress (Pa) or component dict

Return type:

predict_saos(omega, return_components=True)[source]¶

Predict SAOS storage and loss moduli.

Three Maxwell modes (E, D, I) plus amplified permanent plateau (P).

Parameters:

omega (ndarray) – Angular frequency array (rad/s)
return_components (bool) – If True, return (G’, G’’); if False, return |G*|

Returns:

(G’, G’’) or |G*|

Return type:

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

Simulate startup shear flow with dual TST feedback.

Parameters:

t (ndarray) – Time array (s)
gamma_dot (float) – Applied shear rate (1/s)
return_full (bool) – If True, return dict with all trajectories

Returns:

Stress array or full trajectory dict

Return type:

simulate_relaxation(t, gamma_step=1.0, return_full=False)[source]¶

Simulate stress relaxation after step strain.

Parameters:

t (ndarray) – Time array after step (s)
gamma_step (float) – Applied step strain
return_full (bool) – If True, return full trajectory dict

Returns:

G(t) relaxation modulus or trajectory dict

Return type:

simulate_creep(t, sigma_0, return_full=False)[source]¶

Simulate creep under constant stress.

Parameters:

t (ndarray) – Time array (s)
sigma_0 (float) – Applied constant stress (Pa)
return_full (bool) – If True, return full trajectory dict

Returns:

Strain gamma(t) or trajectory dict

Return type:

simulate_laos(t, gamma_0, omega)[source]¶

Simulate LAOS (Large Amplitude Oscillatory Shear).

Parameters:

t (ndarray) – Time array (s)
gamma_0 (float) – Strain amplitude
omega (float) – Angular frequency (rad/s)

Returns:

Keys: time, strain, stress, gamma_dot, N1, mu_E_xy, mu_E_nat_xy, mu_D_xy, mu_I_xy, mu_I_nat_xy, damage, damage_int

Return type:

predict_normal_stresses(gamma_dot)[source]¶

Predict steady-state normal stress differences.

At steady state, E and I networks contribute zero normal stress. Only D-network contributes N1.

N1 = 2 * G_D * (gamma_dot / k_d_D)^2 N2 = 0

Parameters:: gamma_dot (ndarray) – Shear rate array (1/s)
Returns:: (N1, N2) arrays (Pa)
Return type:: tuple[ndarray, ndarray]

extract_laos_harmonics(laos_result, n_harmonics=5)[source]¶

Extract Fourier harmonics from LAOS simulation.

Parameters:

laos_result (dict[str, ndarray]) – Output from simulate_laos()
n_harmonics (int) – Number of harmonics to extract

Returns:

Keys: harmonic_index (1, 3, 5, …), sigma_harmonics

Return type:

get_payne_parameters()[source]¶

Extract Payne effect parameters from model.

The Payne effect manifests as modulus drop with increasing strain amplitude, driven by interphase softening.

Returns:: G_0: zero-strain modulus G_inf: high-strain modulus (X*G_P only) gamma_c: approximate critical strain (1/X_I)
Return type:: dict[str, float]

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

NumPyro/BayesianMixin model function for HVNM.

Routes to appropriate JAX-traceable prediction based on test_mode. Required by BayesianMixin for NumPyro NUTS sampling.

Parameters:

X (array-like) – Independent variable
params (array-like) – Parameter values in ParameterSet order
test_mode (str, optional) – Override stored test mode
**kwargs – Protocol-specific: gamma_dot, sigma_applied, gamma_0, omega

Returns:

Predicted response

Return type:

jnp.ndarray

classmethod unfilled_vitrimer(G_P=10000.0, G_E=10000.0, G_D=1000.0, nu_0=10000000000.0, E_a=80000.0, V_act=1e-05, T=300.0, k_d_D=1.0)[source]¶

Create unfilled vitrimer (phi=0, recovers HVM exactly).

Parameters:

G_P (float) – Subnetwork moduli (Pa)
G_E (float) – Subnetwork moduli (Pa)
G_D (float) – Subnetwork moduli (Pa)
nu_0 (float) – TST parameters
E_a (float) – TST parameters
V_act (float) – TST parameters
T (float) – TST parameters
k_d_D (float) – Dissociative rate (1/s)

Returns:

Model with phi=0 (no interphase contribution)

Return type:

classmethod filled_elastomer(G_P=10000.0, phi=0.1, R_NP=2e-08, delta_m=1e-08)[source]¶

Create filled elastomer (no exchange networks).

Parameters:

G_P (float) – Permanent network modulus (Pa)
phi (float) – NP volume fraction
R_NP (float) – NP radius (m)
delta_m (float) – Mobile interphase thickness (m)

Returns:

Model with only amplified P-network (no E, D, or active I)

Return type:

classmethod partial_vitrimer_nc(G_P=10000.0, G_E=10000.0, phi=0.1, nu_0=10000000000.0, E_a=80000.0, V_act=1e-05, T=300.0, **nc_kwargs)[source]¶

Create partial vitrimer nanocomposite (G_D=0).

Parameters:

G_P (float) – Network moduli (Pa)
G_E (float) – Network moduli (Pa)
phi (float) – NP volume fraction
nu_0 (float) – TST parameters
E_a (float) – TST parameters
V_act (float) – TST parameters
T (float) – TST parameters
**nc_kwargs – NP geometry: R_NP, delta_m, beta_I, nu_0_int, E_a_int, V_act_int

Returns:

Model with P + E + I networks (no D)

Return type:

classmethod conventional_filled_rubber(G_P=10000.0, phi=0.1, R_NP=2e-08, delta_m=1e-08, G_D=1000.0, k_d_D=1.0)[source]¶

Create conventional filled rubber (no E-network, frozen interphase).

Parameters:

G_P (float) – Permanent network modulus (Pa)
phi (float) – NP volume fraction
R_NP (float) – NP geometry (m)
delta_m (float) – NP geometry (m)
G_D (float) – Dissociative modulus (Pa)
k_d_D (float) – Dissociative rate (1/s)

Returns:

Model with P + D + frozen I (no exchange)

Return type:

classmethod matrix_only_exchange(G_P=10000.0, G_E=10000.0, phi=0.1, nu_0=10000000000.0, E_a=80000.0, V_act=1e-05, T=300.0, **nc_kwargs)[source]¶

Create model with frozen interphase (k_BER^int=0).

The interphase acts as a dead (non-exchanging) reinforcement layer.

Parameters:

G_P (float) – Network moduli (Pa)
G_E (float) – Network moduli (Pa)
phi (float) – NP volume fraction
nu_0 (float) – Matrix TST parameters
E_a (float) – Matrix TST parameters
V_act (float) – Matrix TST parameters
T (float) – Matrix TST parameters
**nc_kwargs – NP geometry: R_NP, delta_m, beta_I

Returns:

Model with active matrix exchange, frozen interphase

Return type: