VLB Protocol Equations & Derivations

This page provides detailed step-by-step derivations of the VLB model predictions for each rheological protocol, including the multi-network variants. All results assume constant dissociation rate \(k_d\).

For the governing equations and notation, see VLB Transient Network Models.

Simple Shear Geometry

In simple shear the velocity field is \(v_x = \dot{\gamma} y\), \(v_y = v_z = 0\). The velocity gradient tensor is:

\[\begin{split}\mathbf{L} = \begin{pmatrix} 0 & \dot{\gamma} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\end{split}\]

The distribution tensor \(\boldsymbol{\mu}\) is symmetric, so we track four independent components: \(\mu_{xx}, \mu_{yy}, \mu_{zz}, \mu_{xy}\).

Component-wise ODE system:

\[\begin{split}\dot{\mu}_{xx} &= k_d(1 - \mu_{xx}) + 2\dot{\gamma}\mu_{xy} \\ \dot{\mu}_{yy} &= k_d(1 - \mu_{yy}) \\ \dot{\mu}_{zz} &= k_d(1 - \mu_{zz}) \\ \dot{\mu}_{xy} &= -k_d \mu_{xy} + \dot{\gamma}\mu_{yy}\end{split}\]

Observations:

• \(\mu_{yy}\) and \(\mu_{zz}\) decouple from the shear components

• Starting from \(\boldsymbol{\mu}(0) = \mathbf{I}\), we have \(\mu_{yy}(t) = \mu_{zz}(t) = 1\) for all \(t\)

• \(\mu_{xy}\) is driven by \(\dot{\gamma}\mu_{yy}\)

• \(\mu_{xx}\) is driven by \(2\dot{\gamma}\mu_{xy}\)

Flow Curve Derivation

Objective: Find steady-state \(\sigma_{12}(\dot{\gamma})\) and \(N_1(\dot{\gamma})\).

Step 1: Set \(\dot{\boldsymbol{\mu}} = 0\):

\[\begin{split}0 &= k_d(1 - \mu_{xx}^{ss}) + 2\dot{\gamma}\mu_{xy}^{ss} \\ 0 &= -k_d \mu_{xy}^{ss} + \dot{\gamma}\end{split}\]

Step 2: Solve for \(\mu_{xy}^{ss}\) from the second equation:

\[\mu_{xy}^{ss} = \frac{\dot{\gamma}}{k_d}\]

Step 3: Substitute into the first equation:

\[\mu_{xx}^{ss} = 1 + \frac{2\dot{\gamma}^2}{k_d^2}\]

Step 4: Since \(\mu_{yy}^{ss} = 1\):

\[\begin{split}\sigma_{12} &= G_0 \mu_{xy}^{ss} = G_0 \frac{\dot{\gamma}}{k_d} = \eta_0 \dot{\gamma} \\ N_1 &= G_0(\mu_{xx}^{ss} - \mu_{yy}^{ss}) = \frac{2G_0\dot{\gamma}^2}{k_d^2}\end{split}\]

Verification conditions:

• \(\sigma_{12}\) is linear in \(\dot{\gamma}\) (Newtonian)

• \(N_1 \propto \dot{\gamma}^2\) (quadratic, typical of Maxwell-type models)

• \(\eta_0 = G_0/k_d\) (Green-Kubo relation)

• \(\Psi_1 = N_1/\dot{\gamma}^2 = 2G_0/k_d^2 = 2\eta_0 t_R\) (first normal stress coefficient)

Startup Shear Derivation

Objective: Find \(\sigma_{12}(t)\) and \(N_1(t)\) for constant \(\dot{\gamma}\) starting from equilibrium.

Initial conditions: \(\mu_{xy}(0) = 0\), \(\mu_{xx}(0) = 1\), \(\mu_{yy}(0) = 1\).

Step 1: Solve the \(\mu_{xy}\) ODE (linear first-order with constant coefficients, \(\mu_{yy} = 1\)):

\[\dot{\mu}_{xy} + k_d \mu_{xy} = \dot{\gamma}\]

Solution:

\[\mu_{xy}(t) = \frac{\dot{\gamma}}{k_d}\left(1 - e^{-k_d t}\right)\]

Hence:

\[\sigma_{12}(t) = \frac{G_0\dot{\gamma}}{k_d}\left(1 - e^{-k_d t}\right)\]

Step 2: Solve the \(\mu_{xx}\) ODE with the known \(\mu_{xy}(t)\):

\[\dot{\mu}_{xx} + k_d \mu_{xx} = k_d + \frac{2\dot{\gamma}^2}{k_d}\left(1 - e^{-k_d t}\right)\]

This is a linear ODE whose solution gives:

\[N_1(t) = G_0(\mu_{xx}(t) - 1) = \frac{2G_0\dot{\gamma}^2}{k_d^2} \left(1 - e^{-k_d t}\right) - \frac{2G_0\dot{\gamma}^2}{k_d} t \, e^{-k_d t}\]

Verification conditions:

• \(\sigma_{12}(0) = 0\) (starts from rest)

• \(\sigma_{12}(\infty) = G_0\dot{\gamma}/k_d = \eta_0\dot{\gamma}\) (steady state)

• \(\sigma_{12}(t)\) is monotonically increasing (no overshoot for constant \(k_d\))

• Time to reach 63% of steady state: \(t_{63\%} = 1/k_d = t_R\)

• \(N_1(t)\) can be non-monotonic for \(\text{Wi} = \dot{\gamma}/k_d > 1\) (the \(t \, e^{-k_d t}\) term creates a transient undershoot)

Stress Relaxation Derivation

Objective: Find \(G(t)\) after step strain \(\gamma_0\).

Approach: At \(t = 0^+\), the material has been instantaneously strained. For \(t > 0\), \(\dot{\gamma} = 0\).

Step 1: With \(\dot{\gamma} = 0\), the ODE simplifies:

\[\dot{\mu}_{xy} = -k_d \mu_{xy}\]

Solution from initial condition \(\mu_{xy}(0^+) = \gamma_0\) (affine deformation in the linear regime):

\[\mu_{xy}(t) = \gamma_0 \, e^{-k_d t}\]

Step 2: The relaxation modulus:

\[G(t) = \frac{\sigma_{12}(t)}{\gamma_0} = \frac{G_0 \mu_{xy}(t)}{\gamma_0} = G_0 \, e^{-k_d t}\]

Verification conditions:

• \(G(0) = G_0\) (instantaneous modulus)

• \(G(\infty) = 0\) (liquid, complete relaxation)

• \(\ln G(t)\) is linear with slope \(-k_d\)

• Area under \(G(t)\): \(\int_0^\infty G(t)\,dt = G_0/k_d = \eta_0\)

Multi-network:

\[G(t) = G_e + \sum_{I=0}^{M-1} G_I \, e^{-k_{d,I} t}\]

Creep Derivation

Objective: Find \(J(t) = \gamma(t)/\sigma_0\) under constant stress \(\sigma_0\).

Single network (Maxwell creep):

Since the VLB single network with constant \(k_d\) is exactly Maxwell:

\[J(t) = \frac{1}{G_0} + \frac{t}{\eta_0} = \frac{1 + k_d t}{G_0}\]

This can also be derived from the stress constraint \(\sigma_0 = G_0 \mu_{xy}\), giving \(\mu_{xy} = \sigma_0/G_0\), and integrating \(\dot{\gamma} = k_d \mu_{xy}\) to obtain \(\gamma(t) = \sigma_0/G_0 + (\sigma_0/\eta_0)t\).

Verification conditions:

• \(J(0) = 1/G_0\) (instantaneous elastic compliance)

• \(dJ/dt = 1/\eta_0 = k_d/G_0\) (constant viscous flow rate)

• \(J(t)\) is linear — hallmark of Maxwell creep

Dual-network (1 transient + permanent) — Standard Linear Solid creep:

\[J(t) = \frac{1}{G_0 + G_e} + \frac{G_0}{G_e(G_0 + G_e)} \left(1 - e^{-t/\tau_{ret}}\right)\]

where the retardation time is:

\[\tau_{ret} = \frac{G_0 + G_e}{G_e \cdot k_d}\]

Verification conditions:

• \(J(0) = 1/(G_0 + G_e)\) (unrelaxed compliance)

• \(J(\infty) = 1/G_e\) (relaxed compliance)

• Bounded creep (solid-like behavior due to \(G_e > 0\))

General multi-network creep requires ODE integration because the stress constraint \(\sigma_0 = \sum G_I \mu_{xy,I} + \eta_s \dot{\gamma}\) creates implicit coupling between modes. RheoJAX solves this via diffrax with fixed-point stress balance at each time step.

SAOS Derivation

Objective: Find \(G'(\omega)\) and \(G''(\omega)\).

Step 1: Assume \(\gamma(t) = \gamma_0 e^{i\omega t}\), so \(\dot{\gamma}(t) = i\omega \gamma_0 e^{i\omega t}\).

Step 2: Seek solution \(\mu_{xy}(t) = \hat{\mu}_{xy} e^{i\omega t}\):

\[i\omega \hat{\mu}_{xy} = -k_d \hat{\mu}_{xy} + i\omega \gamma_0\]

\[\hat{\mu}_{xy} = \frac{i\omega \gamma_0}{k_d + i\omega}\]

Step 3: Complex modulus:

\[G^*(\omega) = \frac{G_0 \hat{\mu}_{xy}}{\gamma_0} = \frac{G_0 i\omega}{k_d + i\omega} = G_0 \frac{i\omega t_R}{1 + i\omega t_R}\]

Step 4: Separate real and imaginary parts:

\[\begin{split}G'(\omega) &= G_0 \frac{\omega^2 t_R^2}{1 + \omega^2 t_R^2} \\ G''(\omega) &= G_0 \frac{\omega t_R}{1 + \omega^2 t_R^2}\end{split}\]

Limiting behaviors:

Regime	\(G'\)	\(G''\)
\(\omega \ll k_d\)	\(G_0 \omega^2/k_d^2\) (slope 2 on log-log)	\(G_0 \omega/k_d\) (slope 1 on log-log)
\(\omega = k_d\)	\(G_0/2\) (crossover)	\(G_0/2\) (crossover)
\(\omega \gg k_d\)	\(G_0\) (plateau)	\(G_0 k_d/\omega\) (slope -1)

Cole-Cole plot:

In the \((G', G'')\) plane, the Maxwell model traces a semicircle:

\[\left(G' - \frac{G_0}{2}\right)^2 + G''^2 = \left(\frac{G_0}{2}\right)^2\]

This provides a visual diagnostic: deviation from a semicircle indicates non-Maxwell behavior.

LAOS Implementation

Under \(\gamma(t) = \gamma_0 \sin(\omega t)\), the full 4-component ODE is integrated numerically via diffrax.Tsit5 with adaptive step size:

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

ODE system (state = [μ_xx, μ_yy, μ_zz, μ_xy]):

\[\begin{split}\dot{\mu}_{xx} &= k_d(1 - \mu_{xx}) + 2\gamma_0\omega\cos(\omega t)\mu_{xy} \\ \dot{\mu}_{yy} &= k_d(1 - \mu_{yy}) \\ \dot{\mu}_{zz} &= k_d(1 - \mu_{zz}) \\ \dot{\mu}_{xy} &= -k_d\mu_{xy} + \gamma_0\omega\cos(\omega t)\mu_{yy}\end{split}\]

Since \(\mu_{yy} = 1\) and \(\mu_{xy}\) satisfies a linear ODE, the shear stress \(\sigma_{12} = G_0\mu_{xy}\) is purely linear (fundamental frequency only, no higher harmonics).

Normal stress \(N_1 = G_0(\mu_{xx} - 1)\) is driven by \(2\dot{\gamma}\mu_{xy}\), which is a product of two oscillatory quantities and therefore contains \(2\omega\) harmonics:

\[\mu_{xx}(t) - 1 \approx A_0 + A_2 \cos(2\omega t) + B_2 \sin(2\omega t) + \ldots\]

RheoJAX implementation:

• Simulates \(n\) complete cycles (default 10)

• Discards first 5 cycles (transient decay)

• Returns \(\sigma_{12}(t)\), \(N_1(t)\), and FFT harmonics

• Verifies linearity via \(I_3/I_1\) ratio for \(\sigma_{12}\)

Uniaxial Extension

Velocity gradient:

\[\begin{split}\mathbf{L} = \begin{pmatrix} \dot{\varepsilon} & 0 & 0 \\ 0 & -\dot{\varepsilon}/2 & 0 \\ 0 & 0 & -\dot{\varepsilon}/2 \end{pmatrix}\end{split}\]

Component equations (axial + transverse only, by symmetry):

\[\begin{split}\dot{\mu}_{11} &= k_d(1 - \mu_{11}) + 2\dot{\varepsilon}\mu_{11} \\ \dot{\mu}_{22} &= k_d(1 - \mu_{22}) - \dot{\varepsilon}\mu_{22}\end{split}\]

Steady state:

\[\begin{split}\mu_{11}^{ss} &= \frac{k_d}{k_d - 2\dot{\varepsilon}} \quad (\dot{\varepsilon} < k_d/2) \\ \mu_{22}^{ss} &= \frac{k_d}{k_d + \dot{\varepsilon}}\end{split}\]

Extensional stress:

\[\sigma_E = G_0(\mu_{11}^{ss} - \mu_{22}^{ss}) = G_0\dot{\varepsilon}\left(\frac{1}{k_d - 2\dot{\varepsilon}} + \frac{1}{k_d + \dot{\varepsilon}}\right)\]

Transient response from equilibrium:

\[\begin{split}\mu_{11}(t) &= \frac{k_d}{k_d - 2\dot{\varepsilon}} \left(1 - e^{-(k_d - 2\dot{\varepsilon})t}\right) + e^{-(k_d - 2\dot{\varepsilon})t} \\ \mu_{22}(t) &= \frac{k_d}{k_d + \dot{\varepsilon}} \left(1 - e^{-(k_d + \dot{\varepsilon})t}\right) + e^{-(k_d + \dot{\varepsilon})t}\end{split}\]

Singularity at \(\dot{\varepsilon} = k_d/2\):

The axial component \(\mu_{11} \to \infty\) — chains cannot relax fast enough to compensate the extensional stretching. This divergence is:

• Regularized by finite extensibility (Langevin chains)

• Analogous to the coil-stretch transition in dilute solutions

• A useful diagnostic: if data shows extensional hardening, VLB with constant \(k_d\) is insufficient

General Analytical Solution

The full tensorial evolution equation (1) (see VLB Transient Network Models) admits a compact integral-form solution through the substitution:

\[\mathbf{A} = \mathbf{F}^{-1} \cdot \boldsymbol{\mu} \cdot \mathbf{F}^{-T}\]

where \(\mathbf{F}(t)\) is the deformation gradient. Under this change of variables, the convective terms cancel and the evolution reduces to:

\[\dot{\mathbf{A}} = k_d\!\left(\mathbf{F}^{-1}\mathbf{F}^{-T} - \mathbf{A}\right)\]

This is a first-order linear ODE in \(\mathbf{A}\) with solution:

\[\mathbf{A}(t) = \mathbf{A}(0) e^{-k_d t} + k_d \int_0^t e^{-k_d(t-s)} \mathbf{F}^{-1}(s) \mathbf{F}^{-T}(s) \, ds\]

Transforming back:

\[\boldsymbol{\mu}(t) = \mathbf{F}(t) \cdot \mathbf{A}(t) \cdot \mathbf{F}^T(t)\]

This integral form is the basis for all closed-form solutions in the VLB model. For constant velocity gradient \(\mathbf{L}\), the deformation gradient is \(\mathbf{F}(t) = e^{\mathbf{L} t}\) and the integral can be evaluated analytically for each protocol.

LAOS: Exact Sinusoidal Proof

Theorem: For constant \(k_d\), the VLB shear stress \(\sigma_{12}(t)\) under LAOS is purely sinusoidal at the fundamental frequency \(\omega\), with zero higher harmonics.

Proof: The \(\mu_{xy}\) component satisfies (with \(\mu_{yy} = 1\)):

\[\dot{\mu}_{xy} + k_d \mu_{xy} = \dot{\gamma}(t) = \gamma_0 \omega \cos(\omega t)\]

This is a linear first-order ODE with constant coefficients. The steady-state (post-transient) particular solution is:

\[\mu_{xy}^{ss}(t) = \frac{\gamma_0 \omega}{k_d^2 + \omega^2} \bigl[k_d \cos(\omega t) + \omega \sin(\omega t)\bigr]\]

which is a pure sinusoid at frequency \(\omega\). Since \(\sigma_{12} = G_0 \mu_{xy}\), the shear stress contains only the fundamental frequency. No harmonic generation occurs because the ODE is linear — there is no mechanism to couple \(\omega\) to \(n\omega\).

Corollary: The intrinsic nonlinearity ratio \(I_3/I_1 = 0\) exactly for the VLB model with constant \(k_d\). Any measured \(I_3/I_1 > 0\) in experimental data indicates that \(k_d\) depends on deformation.

Note

This linearity is a property of the constant-\(k_d\) model only. The VLBVariant with Bell breakage (\(k_d = k_d(\boldsymbol{\mu})\)) produces nonlinear coupling and \(I_3/I_1 > 0\).

LAOS: \(N_1\) Harmonic Structure

While \(\sigma_{12}\) is purely sinusoidal, the first normal stress difference \(N_1 = G_0(\mu_{xx} - \mu_{yy}) = G_0(\mu_{xx} - 1)\) has a richer harmonic content.

The \(\mu_{xx}\) equation is driven by \(2\dot{\gamma}\mu_{xy}\), which is a product of the driving frequency and the fundamental response:

\[2\dot{\gamma}(t) \cdot \mu_{xy}^{ss}(t) \propto \cos(\omega t) \cdot [\cos(\omega t) + \ldots]\]

This product generates a DC offset (constant term) and a \(2\omega\) component via the trigonometric identity \(\cos^2(\omega t) = \frac{1}{2}(1 + \cos(2\omega t))\).

Explicit solution (steady state):

\[N_1(t) = \overline{N}_1 + N_1^{(2c)} \cos(2\omega t) + N_1^{(2s)} \sin(2\omega t)\]

where the nonzero mean \(\overline{N}_1 > 0\) reflects the time-averaged normal stress, and the \(2\omega\) oscillation captures the stretch-compression asymmetry within each cycle.

Key features:

• \(N_1\) oscillates at \(2\omega\) (not \(\omega\))

• Nonzero mean \(\overline{N}_1 = 2G_0 \gamma_0^2 \omega^2 / (k_d^2 + \omega^2)^2 \cdot k_d^2\)

• No higher harmonics (\(4\omega, 6\omega, \ldots\)) for constant \(k_d\)

• Amplitude grows as \(\gamma_0^2\) — quadratic in strain amplitude

LAOS Linearity Validation

In practice, verifying \(I_3/I_1 \approx 0\) for the VLB model serves as both a correctness check on the numerical implementation and a diagnostic for the constant-\(k_d\) assumption.

Validation methodology:

1. Numerical integration: Simulate LAOS via the ODE system for 10+ cycles using diffrax.Tsit5 with tight tolerances (rtol=1e-8, atol=1e-10)

2. Transient removal: Discard the first 5 cycles to eliminate initial transients

3. FFT analysis: Compute the Fourier spectrum of \(\sigma_{12}(t)\) over the last 5 cycles

4. Harmonic ratio: \(I_3/I_1 = |a_3|/|a_1|\) where \(a_n\) are Fourier coefficients

Expected results:

• \(I_3/I_1 < 10^{-10}\) (machine precision) for constant \(k_d\)

• \(I_3/I_1 \sim 10^{-3}\) to \(10^{-1}\) for Bell \(k_d\) with moderate \(\nu\)

Diagnostic use: If \(I_3/I_1 > 10^{-6}\) is measured experimentally, the constant-\(k_d\) VLB is insufficient and force-dependent breakage (VLBVariant with breakage="bell") should be considered.

Extensional Singularity and Safety

The steady-state axial component \(\mu_{11}^{ss} = k_d/(k_d - 2\dot{\varepsilon})\) diverges at the critical Weissenberg number:

\[\text{Wi}_{ext} = \frac{\dot{\varepsilon}}{k_d} = \frac{1}{2}\]

This is the extensional catastrophe: chains are stretched faster than they can relax, and the Gaussian chain assumption breaks down.

Physical interpretation:

• For \(\text{Wi}_{ext} < 0.5\): steady state exists, Trouton ratio is finite

• At \(\text{Wi}_{ext} = 0.5\): transition to unbounded stretch

• For \(\text{Wi}_{ext} > 0.5\): no steady state; transient stress grows without bound

Safety clamping in RheoJAX:

The implementation clamps \(\mu_{11}\) to prevent numerical overflow:

# In _kernels.py: safety clamp for extensional flows
mu_11 = jnp.minimum(mu_11, L_max_sq)  # FENE-P bound
# Or for Gaussian chains:
mu_11 = jnp.minimum(mu_11, 1e6)  # numerical safety

Regularization options:

1. FENE-P (VLBVariant(stress_type="fene")): Bounds stress via \(f = L_{\max}^2/(L_{\max}^2 - \text{tr}(\boldsymbol{\mu}) + 3)\), which diverges at finite \(\text{tr}(\boldsymbol{\mu}) = L_{\max}^2 + 3\) rather than at \(\text{Wi}_{ext} = 0.5\)

2. Bell breakage (VLBVariant(breakage="bell")): Enhanced breakage at high stretch softens the singularity but does not fully remove it

For extensional rheology applications, the FENE-P variant is strongly recommended.

Multi-Network Protocol Summary

Protocol	Method	Formula
Flow curve	Analytical	\(\sigma = (\sum G_I/k_{d,I} + \eta_s)\dot{\gamma}\)
Startup	Analytical	\(\sigma(t) = \sum \frac{G_I\dot{\gamma}}{k_{d,I}}(1-e^{-k_{d,I}t}) + \eta_s\dot{\gamma}\)
Relaxation	Analytical	\(G(t) = G_e + \sum G_I e^{-k_{d,I}t}\)
Creep	ODE / analytical (SLS)	General: ODE with stress balance; 1-mode + \(G_e\): SLS formula
SAOS	Analytical	\(G^* = i\omega\eta_s + G_e + \sum G_I\frac{i\omega/k_{d,I}}{1+i\omega/k_{d,I}}\)
LAOS	ODE (diffrax)	Multi-mode ODE with oscillatory driving