VLB Advanced Theory & Numerical Methods¶

This page documents the theoretical foundations for VLB extensions and the numerical methods used in their implementation. Four extensions are implemented in RheoJAX; one (Langevin chains) remains as a theory reference for future development.

For usage, see VLBVariant — Bell, FENE-P & Temperature (Bell, FENE-P, Temperature) and VLBNonlocal — Shear Banding PDE (spatial PDE with banding).

Implementation Status

Extension	Status	Classes
Bell breakage	Implemented	`VLBVariant`, `VLBNonlocal`
FENE-P stress	Implemented	`VLBVariant`, `VLBNonlocal`
Temperature dependence	Implemented	`VLBVariant`
Nonlocal PDE	Implemented	`VLBNonlocal`
Langevin chains	Theory reference	(not yet implemented)

Chain Distribution Decomposition¶

The chain distribution function \(\varphi(\mathbf{r},t)\) admits a factorization into density and shape:

\[\varphi(\mathbf{r},t) = c(t) \cdot P(\mathbf{r},t)\]

where \(c(t) = \int \varphi \, d^3r\) is the total chain density (number of elastically active chains per unit volume) and \(P(\mathbf{r},t)\) is the normalized probability density of the end-to-end vector.

Angle-bracket operator: For any quantity \(f(\mathbf{r})\), the ensemble average over the distribution is:

\[\langle f \rangle = \frac{1}{c} \int f(\mathbf{r}) \, \varphi(\mathbf{r},t) \, d^3r = \int f(\mathbf{r}) \, P(\mathbf{r},t) \, d^3r\]

The distribution tensor is therefore:

\[\boldsymbol{\mu} = \frac{3}{\langle r_0^2 \rangle} \langle \mathbf{r} \otimes \mathbf{r} \rangle\]

Master evolution equation: The Smoluchowski equation for \(\varphi\) under affine convection with source/sink kinetics:

\[\frac{\partial \varphi}{\partial t} + \nabla_r \cdot (\mathbf{L} \cdot \mathbf{r} \, \varphi) = k_a \varphi_0 - k_d \varphi\]

Taking the second moment of this equation (multiplying by \(\mathbf{r} \otimes \mathbf{r}\) and integrating) directly yields the distribution tensor evolution equation used in the VLB model.

At equilibrium \(k_a = k_d\) (detailed balance), and the equilibrium distribution \(\varphi_0\) is isotropic Gaussian with \(\boldsymbol{\mu}_{eq} = \mathbf{I}\).

Numerical Methods¶

Semi-Implicit Exponential Integrator¶

For the ODE-based protocols (LAOS, startup/relaxation with nonlinear \(k_d\)), RheoJAX uses the diffrax.Tsit5 explicit Runge-Kutta solver with adaptive stepping. For finite-element or custom integration contexts, a semi-implicit exponential integrator is recommended:

Step 1 — Explicit convective update:

\[\boldsymbol{\mu}^* = \boldsymbol{\mu}^n + \Delta t \bigl(\mathbf{D} \cdot \boldsymbol{\mu}^n + \boldsymbol{\mu}^n \cdot \mathbf{D}\bigr)\]

Step 2 — Exponential kinetic decay:

\[\boldsymbol{\mu}^{n+1} = \mathbf{I} + (\boldsymbol{\mu}^* - \mathbf{I}) e^{-k_d \Delta t}\]

This scheme is exact for the kinetic part (exponential decay toward \(\mathbf{I}\)) and explicit for the convective part, providing unconditional stability for the relaxation term regardless of \(k_d \Delta t\).

Symmetry Enforcement¶

After each time step, enforce the symmetry of the distribution tensor:

\[\boldsymbol{\mu} \leftarrow \frac{1}{2}(\boldsymbol{\mu} + \boldsymbol{\mu}^T)\]

This is necessary because floating-point arithmetic can introduce asymmetric errors, particularly in the off-diagonal components. In JAX:

mu = 0.5 * (mu + mu.T)

Time Step Selection¶

For stability and accuracy, the time step should satisfy:

\[\Delta t < \min\!\left(\frac{1}{|\mathbf{D}|}, \quad \frac{\pi}{10 \omega}\right)\]

where \(|\mathbf{D}|\) is the largest eigenvalue of the deformation rate tensor. The second condition ensures at least 20 time points per oscillation cycle in LAOS/SAOS simulations.

For multi-network models with widely separated \(k_d\) values, the stiffest mode (largest \(k_d\)) sets the stability limit. The semi-implicit scheme avoids this restriction for the kinetic term.

Steady-State Detection¶

For flow curve and other steady-state protocols, convergence is declared when:

\[\frac{|\dot{\boldsymbol{\mu}}|}{|\boldsymbol{\mu}|} < \varepsilon\]

where \(\varepsilon = 10^{-8}\) is the default tolerance and \(|\cdot|\) is the Frobenius norm. For multi-network models, all modes must satisfy this criterion.

Multiple Time-Scale Handling¶

Multi-network VLB models with \(k_d\) values spanning several decades create stiff ODE systems. RheoJAX handles this via:

Adaptive stepping (diffrax.PIDController): automatic step size reduction near fast transients
Mode-ordered initialization: log-spaced \(k_d\) values ensure well-separated time scales
Separate analytical/ODE split: fast modes that have reached steady state can be treated analytically while slow modes continue integrating

Force-Dependent Dissociation Rate¶

Physical Motivation¶

In real networks, the bond breaking rate depends on the force applied to the chain. Under tension, bonds break faster (slip-bond behavior); under compression, they may break slower (catch-bond behavior). This force-dependence introduces shear thinning, stress overshoot, and nonlinear LAOS response.

Bell Model¶

The simplest force-dependent model (Bell, 1978):

\[k_d(\boldsymbol{\mu}) = k_d^0 \exp\!\left(\nu \cdot (\lambda_c - 1)\right)\]

where:

\(k_d^0\) is the unstressed dissociation rate
\(\nu\) is the force sensitivity parameter (dimensionless)
\(\lambda_c = \sqrt{\text{tr}(\boldsymbol{\mu})/3}\) is the normalized average chain stretch

Effect on predictions:

Protocol	Constant \(k_d\)	Bell \(k_d\)
Flow curve	Newtonian	Shear thinning (\(\eta \propto \dot{\gamma}^{n-1}\))
Startup	Monotonic	Stress overshoot
LAOS \(\sigma_{12}\)	Linear (no harmonics)	Nonlinear (\(I_3 > 0\))
Extension	Singularity at \(\dot{\varepsilon} = k_d/2\)	Softening before singularity

New parameter: \(\nu\) (force sensitivity).

\(\nu = 0\): recovers constant \(k_d\)
\(\nu \sim 1\) - 5: moderate force sensitivity
\(\nu > 10\): strong force weakening (catch-bond to slip-bond transition)

Evolution equation becomes nonlinear:

\[\dot{\boldsymbol{\mu}} = k_d(\boldsymbol{\mu})(\mathbf{I} - \boldsymbol{\mu}) + \mathbf{L} \cdot \boldsymbol{\mu} + \boldsymbol{\mu} \cdot \mathbf{L}^T\]

This requires ODE integration for all protocols (no more analytical solutions).

Comparison with TNT Bell:

The VLB Bell model is mathematically equivalent to the TNTSingleMode with breakage="bell" at the Gaussian chain level. The VLB formulation is preferred when extending to Langevin chains (see below).

Slip-Bond and Catch-Bond¶

Slip-bond (\(k_d\) increases with force):

\[k_d(F) = k_d^0 \exp\!\left(\frac{F \cdot \delta}{k_B T}\right)\]

where \(\delta\) is the distance to the transition state.

Catch-bond (\(k_d\) first decreases then increases with force):

\[k_d(F) = k_{slip} \exp\!\left(\frac{F \delta_s}{k_B T}\right) + k_{catch} \exp\!\left(-\frac{F \delta_c}{k_B T}\right)\]

Catch-bond behavior has been observed in biological systems (selectin-ligand, fibrin) and in some synthetic networks with dual-lock mechanisms.

Langevin Chain Finite Extensibility¶

Physical Motivation¶

Gaussian chain statistics assume unlimited extensibility. Real polymer chains have a finite contour length \(L_{max} = N_K b\) (product of Kuhn segments and Kuhn length). As chains approach full extension, the restoring force diverges.

Langevin Force Law¶

For a freely jointed chain with \(N_K\) segments, the force-extension relation is given by the inverse Langevin function:

\[\mathbf{F} = \frac{k_B T}{b} \mathcal{L}^{-1}\!\left(\frac{r}{N_K b}\right) \hat{\mathbf{r}}\]

where \(\mathcal{L}(x) = \coth(x) - 1/x\) is the Langevin function.

Stress modification:

\[\boldsymbol{\sigma} = G_0 \, f\!\left(\text{tr}(\boldsymbol{\mu})\right) (\boldsymbol{\mu} - \mathbf{I})\]

where \(f\) is a strain-amplification (FENE-like) function:

\[f(\text{tr}(\boldsymbol{\mu})) = \frac{L_{max}^2}{L_{max}^2 - \text{tr}(\boldsymbol{\mu}) + 3}\]

Effects:

Bounds the extensional stress (regularizes the \(\dot{\varepsilon} = k_d/2\) singularity)
Strain hardening at large deformations
Affects \(N_1\) predictions and LAOS harmonics

New parameter: \(L_{max}\) (maximum chain extensibility, typically 3-100).

Comparison with TNT FENE-P:

The Langevin/FENE-P stress modification is identical to the FENE variant of TNTSingleMode. The VLB formulation with the distribution tensor provides a more natural route to combining FENE with force-dependent \(k_d\).

Padé Approximation¶

The inverse Langevin function is typically approximated by a Padé form for computational efficiency:

\[\mathcal{L}^{-1}(x) \approx x \frac{3 - x^2}{1 - x^2}\]

This is exact at \(x = 0\) and \(x \to 1\), and accurate to within 4% for all \(x \in [0, 1)\).

Nonlocal Formulation¶

Physical Motivation¶

In spatially heterogeneous flows (e.g., Couette, pipe, channel), the distribution tensor can vary across the gap. Shear banding — the coexistence of high-shear and low-shear bands — is a common phenomenon in transient networks with force-dependent \(k_d\).

Nonlocal Evolution Equation¶

The distribution tensor acquires a spatial diffusion term:

\[\frac{\partial \boldsymbol{\mu}}{\partial t} = k_d(\boldsymbol{\mu})(\mathbf{I} - \boldsymbol{\mu}) + \mathbf{L} \cdot \boldsymbol{\mu} + \boldsymbol{\mu} \cdot \mathbf{L}^T + D_\mu \nabla^2 \boldsymbol{\mu}\]

where \(D_\mu\) is a diffusivity with dimensions m²/s.

Cooperativity length: \(\xi = \sqrt{D_\mu / k_d}\) sets the characteristic width of shear band interfaces.

Coupled 1D system (gap direction \(y\) ):

\[\begin{split}\frac{\partial \mu_{xy}}{\partial t} &= -k_d \mu_{xy} + \dot{\gamma}(y)\mu_{yy} + D_\mu \frac{\partial^2 \mu_{xy}}{\partial y^2} \\ \frac{\partial \dot{\gamma}}{\partial y} &= \text{(momentum balance)}\end{split}\]

This is a PDE system similar to the existing FluidityNonlocal and DMTNonlocal.

New parameters: \(D_\mu\) (diffusivity), \(n_{points}\) (spatial grid), \(L_{gap}\) (gap width).

Temperature Dependence¶

Arrhenius Bond Kinetics¶

The dissociation rate follows Arrhenius kinetics:

\[k_d(T) = k_d^0 \exp\!\left(-\frac{E_a}{k_B T}\right)\]

where \(E_a\) is the activation energy for bond dissociation.

Modulus temperature dependence:

\[G_0(T) = c(T) k_B T\]

where \(c(T)\) is the chain density (may depend on temperature if bonds have different association constants at different temperatures).

Application: Rheological master curves. By fitting \(k_d(T)\) at multiple temperatures, one obtains \(E_a\), which can be compared with calorimetric data (DSC) or molecular simulations.

WLF/VFT Kinetics¶

For glass-forming systems, the Arrhenius model breaks down and one uses the Williams-Landel-Ferry (WLF) or Vogel-Fulcher-Tammann (VFT) form:

\[\log \frac{k_d(T)}{k_d(T_r)} = \frac{-C_1(T - T_r)}{C_2 + T - T_r}\]

Connection to Finite Elements¶

The VLB evolution equation is a first-order tensorial ODE at each material point. In a finite element (FE) context:

Gauss-point state: Each integration point stores \(\boldsymbol{\mu}\) as internal variables
Stress computation: Given \(\boldsymbol{\mu}\) and \(\mathbf{L}\), compute \(\boldsymbol{\sigma} = G_0(\boldsymbol{\mu} - \mathbf{I})\)
State update: Integrate \(\dot{\boldsymbol{\mu}}\) over the time step using implicit or semi-implicit schemes
Algorithmic tangent: The material tangent \(\partial \boldsymbol{\sigma}/\partial \boldsymbol{\varepsilon}\) is needed for Newton-Raphson iteration in the FE solver

The VLB formulation is particularly well-suited for FE implementation because \(\boldsymbol{\mu}\) has a clear physical interpretation and the evolution equation is relatively simple.

Integration schemes:

Exponential map (exact for constant \(\mathbf{L}\)): \(\boldsymbol{\mu}^{n+1} = \mathbf{I} + (\boldsymbol{\mu}^n - \mathbf{I}) e^{-k_d \Delta t} + \ldots\)
Semi-implicit (stable for large \(\Delta t\)): kinetic term implicit, convective term explicit
Fully implicit (most stable, requires tangent linearization)

The JAX automatic differentiation framework in RheoJAX provides the algorithmic tangent “for free” via jax.jacfwd or jax.jacrev.

Summary of Extensions¶

Extension	New Params	New Physics	Status
Bell \(k_d\)	\(\nu\)	Shear thinning, overshoot, LAOS harmonics	Implemented (`VLBVariant`)
Langevin/FENE	\(L_{max}\)	Extensional hardening, bounded stress	Implemented (`VLBVariant`)
Nonlocal	\(D_\mu\), \(n_{pts}\), \(L_{gap}\)	Shear banding	Implemented (`VLBNonlocal`)
Temperature	\(E_a\) or \(C_1, C_2, T_r\)	TTS, Arrhenius/WLF	Implemented (`VLBVariant`)
Catch-bond	\(k_s, k_c, \delta_s, \delta_c\)	Non-monotonic \(k_d(F)\)	Future
FE coupling	(mesh)	Spatial mechanics	External (via JAX tangent)

References¶

Vernerey, F.J. (2018). “Transient response of nonlinear polymer networks: A kinetic theory.” J. Mech. Phys. Solids, 115, 230-247. https://doi.org/10.1016/j.jmps.2018.02.018 PDF
Bell, G.I. (1978). “Models for the specific adhesion of cells to cells.” Science, 200(4342), 618-627. https://doi.org/10.1126/science.347575
Rubinstein, M. & Semenov, A.N. (2001). “Dynamics of entangled solutions of associating polymers.” Macromolecules, 34(4), 1058-1068. https://doi.org/10.1021/ma0013049
Meng, F. & Pritchard, R.H. & Terentjev, E.M. (2016). “Stress relaxation, dynamics, and plasticity of transient polymer networks.” Macromolecules, 49(7), 2843-2852. https://doi.org/10.1021/acs.macromol.5b02667
Long, R., Qi, H.J. & Dunn, M.L. (2013). “Modeling the mechanics of covalently adaptable polymer networks with temperature-dependent bond exchange reactions.” Soft Matter, 9, 4083-4096. https://doi.org/10.1039/C3SM27945F
Hui, C.-Y., Cui, F., Zehnder, A. & Vernerey, F.J. (2021). “Physically motivated models of polymer networks with dynamic cross-links: comparative study and future outlook.” Proc. R. Soc. A, 477, 20210608. DOI: 10.1098/rspa.2021.0608 PDF