VLB Model Equation Verification

Date: 2026-03-30 Sources: Vernerey, Long & Brighenti (2017) JMPS 107:1-20; Vernerey (2018) JMPS 115:217-247 (PMC6824477); Vernerey et al. (2018) Polymers 10(8):848 (PMC6403683); Bell (1978) Science 200:618-627.

---

1. Distribution Tensor Evolution Equation

Paper (Eq. 3b in Polymers 2018; Eq. 17 in JMPS 2018)

The full evolution equation from Vernerey (2018 JMPS, Eq. 17) is:

D(c·μ)/Dt = k_a·(C - c)·I - k_d·(c·μ) + c·(L·μ) + c·(L·μ)^T

where c = attached chain concentration, C = total concentration, k_a = association rate, k_d = dissociation rate, L = velocity gradient tensor.

For incompressible material with constant concentration (c = const, C = const), and using the steady-state relation c = C·k_a/(k_a + k_d), this simplifies. In the simplest case where k_a = k_d (equal rates) and c = C/2, or more generally absorbing concentration into the effective rates:

dμ/dt = k_a·((C-c)/c)·I - k_d·μ + L·μ + μ·L^T

At dynamic equilibrium of concentrations, k_a·(C - c)/c = k_d, giving:

dμ/dt = k_d·(I - μ) + L·μ + μ·L^T

RheoJAX Implementation

dmu/dt = k_d*(I - mu) + L·mu + mu·L^T

VERDICT: CORRECT. The implementation uses the full velocity gradient L (not just symmetric D), matching the paper. The docstring at the top of _kernels.py incorrectly states D·mu + mu·D but the actual code uses L·mu + mu·L^T which is correct.

Note: The docstring header says D·mu + mu·D (symmetric part only) but the derivation in the comments and the actual code at lines 127-131 correctly use the full velocity gradient. The docstring should be corrected to say L·mu + mu·L^T.

---

2. Cauchy Stress

Paper (Eq. 28 in JMPS 2018, Gaussian limit N→∞)

For Gaussian chains (N → ∞), the stress from Eq. 28 reduces to:

σ = c·k_B·T·(μ - I) + p·I

where G₀ = c·k_B·T is the network modulus.

RheoJAX Implementation

sigma_xy = G0 * mu_xy           # (since I_xy = 0)
N1 = G0 * (mu_xx - mu_yy)       # = sigma_xx - sigma_yy

VERDICT: CORRECT. The shear stress σ_xy = G₀·μ_xy and N1 = G₀·(μ_xx - μ_yy) are correct for the Gaussian limit.

---

3. Free Energy

Paper (Eq. 5 in Polymers 2018; Eq. 26/27 in JMPS 2018)

For Gaussian chains, the stored elastic energy (from the leading term of Eq. 5 in Polymers 2018):

ΔΨ_e = (1/2)·c·k_B·T·[tr(μ) - 3 - ln det(μ)]
     = (1/2)·G₀·[tr(μ) - 3 - ln det(μ)]

This is the Neo-Hookean form. The quantity [tr(μ) - 3 - ln det(μ)] ≥ 0 with equality only at μ = I.

RheoJAX Implementation

NOT IMPLEMENTED. No vlb_free_energy function exists in _kernels.py.

VERDICT: MISSING (not a bug — the free energy is not needed for stress/strain predictions, only for thermodynamic consistency checks).

---

4. Dissipation

Paper (from Eqs. 19-21 in JMPS 2018)

For Gaussian chains, the dissipation is:

D = G₀·k_d·[tr(μ) - 3 - ln det(μ)] ≥ 0

This is 2·k_d·ΔΨ_e, confirming non-negative dissipation.

RheoJAX Implementation

NOT IMPLEMENTED. No vlb_dissipation function exists.

VERDICT: MISSING (same rationale as free energy — not needed for predictions).

---

5. Simple Shear Flow Components

Paper derivation

With L = [[0, γ̇, 0], [0, 0, 0], [0, 0, 0]] and L^T = [[0, 0, 0], [γ̇, 0, 0], [0, 0, 0]]:

(L·μ)_xx = γ̇·μ_xy,     (μ·L^T)_xx = μ_xy·γ̇     → 2γ̇·μ_xy
(L·μ)_yy = 0,           (μ·L^T)_yy = 0           → 0
(L·μ)_zz = 0,           (μ·L^T)_zz = 0           → 0
(L·μ)_xy = γ̇·μ_yy,     (μ·L^T)_xy = 0           → γ̇·μ_yy

Therefore:

dμ_xx/dt = k_d·(1 - μ_xx) + 2·γ̇·μ_xy
dμ_yy/dt = k_d·(1 - μ_yy)
dμ_zz/dt = k_d·(1 - μ_zz)
dμ_xy/dt = -k_d·μ_xy + γ̇·μ_yy

RheoJAX Implementation (lines 127-131)

dmu_xx = k_d * (1.0 - mu_xx) + 2.0 * gamma_dot * mu_xy
dmu_yy = k_d * (1.0 - mu_yy)
dmu_zz = k_d * (1.0 - mu_zz)
dmu_xy = -k_d * mu_xy + gamma_dot * mu_yy

VERDICT: CORRECT. Matches exactly. This is equivalent to the upper-convected Maxwell model.

---

6. Shear Stress and N1

Paper (Gaussian limit)

σ_xy = G₀·(μ_xy - I_xy) = G₀·μ_xy      (since I_xy = 0)
N1 = G₀·(μ_xx - μ_yy)

RheoJAX Implementation

def vlb_shear_stress(mu_xy, G0):
    return G0 * mu_xy

def vlb_normal_stress_1(mu_xx, mu_yy, G0):
    return G0 * (mu_xx - mu_yy)

VERDICT: CORRECT.

---

7. Steady-State Shear

Paper derivation

At steady state (dμ/dt = 0):

• μ_yy = 1 (from k_d·(1-μ_yy) = 0)

• μ_xy = γ̇·μ_yy/k_d = γ̇/k_d

• σ = G₀·γ̇/k_d = η₀·γ̇ where η₀ = G₀/k_d

This is Newtonian with zero-shear viscosity η₀ = G₀/k_d.

RheoJAX Implementation

def vlb_steady_shear(gamma_dot, G0, k_d):
    return G0 * gamma_dot / k_d

VERDICT: CORRECT. Newtonian response for constant k_d.

---

8. SAOS (Small Amplitude Oscillatory Shear)

Paper derivation

The VLB model in the Gaussian limit reduces to the Maxwell model. For a Maxwell element with modulus G₀ and relaxation time τ = 1/k_d:

G'(ω) = G₀·ω²τ²/(1 + ω²τ²) = G₀·ω²/(ω² + k_d²)
G''(ω) = G₀·ωτ/(1 + ω²τ²) = G₀·ω·k_d/(ω² + k_d²)

RheoJAX Implementation

t_R = 1.0 / k_d
wt = omega * t_R
wt2 = wt * wt
denom = 1.0 + wt2
G_prime = G0 * wt2 / denom        # = G0·ω²/(ω² + k_d²)
G_double_prime = G0 * wt / denom   # = G0·ω·k_d/(ω² + k_d²)  [since wt/denom = ωτ/(1+ω²τ²)]

Wait — let me verify: wt/denom = (ω/k_d)/(1 + ω²/k_d²) = ω·k_d/(k_d² + ω²). So G'' = G₀·ω·k_d/(ω² + k_d²). But we can also write this as G₀·ω/(k_d + ω²/k_d).

The formula G''(ω) = G₀·ω·k_d/(ω² + k_d²) matches the standard Maxwell form.

VERDICT: CORRECT.

---

9. Stress Relaxation

Paper

Single exponential decay (Maxwell element):

G(t) = G₀·exp(-k_d·t) = G₀·exp(-t/τ)

where τ = 1/k_d.

RheoJAX Implementation

def vlb_relaxation_modulus(t, G0, k_d):
    return G0 * jnp.exp(-k_d * t)

VERDICT: CORRECT.

---

10. Creep Compliance

Analysis

For a Maxwell element (spring G₀ in series with dashpot η = G₀/k_d), the creep compliance under constant stress σ₀ is:

J(t) = 1/G₀ + t/η = 1/G₀ + k_d·t/G₀ = (1 + k_d·t)/G₀

This is not J(t) = (1/G₀)(1 - exp(-k_d·t)) + t·k_d/G₀ — that would be the Kelvin-Voigt or SLS form.

The Maxwell creep: instantaneous elastic strain 1/G₀ followed by linear viscous flow t/η.

RheoJAX Implementation

def vlb_creep_compliance_single(t, G0, k_d):
    return (1.0 + k_d * t) / G0

VERDICT: CORRECT. This is the standard Maxwell creep compliance. The question mark in the user’s query is resolved: J(t) = (1 + k_d·t)/G₀ is correct for Maxwell. The exponential form would apply to a Kelvin-Voigt element or Standard Linear Solid, which is implemented separately as vlb_creep_compliance_dual for the case with a permanent elastic network (G_e > 0).

---

11. Bell Model Force-Dependent Dissociation

Paper (Bell 1978, applied in Vernerey framework)

Bell (1978) proposed for bond dissociation under force:

k_off(f) = k_off_0 · exp(f·γ / k_B·T)

In the Vernerey VLB context, this is adapted using chain stretch rather than force directly. The mean chain stretch is λ̄ = sqrt(tr(μ)/3), and the Bell-type relation becomes:

k_d(μ) = k_d_0 · exp(ν · (λ̄ - 1))

where ν absorbs the force sensitivity. At equilibrium (tr(μ) = 3, λ̄ = 1), k_d = k_d_0.

RheoJAX Implementation

def vlb_breakage_bell(mu_xx, mu_yy, mu_zz, k_d_0, nu):
    tr_mu = mu_xx + mu_yy + mu_zz
    stretch = jnp.sqrt(jnp.maximum(tr_mu / 3.0, 1e-30))
    return k_d_0 * jnp.exp(nu * (stretch - 1.0))

VERDICT: CORRECT. The stretch definition √(tr(μ)/3) and the Bell-type exponential are consistent with the VLB framework. The 1e-30 floor prevents sqrt(0) gradient issues.

---

12. FENE-P Finite Extensibility

Paper context

The FENE-P (Peterlin) closure introduces a nonlinear spring factor that diverges as chain extension approaches its limit. For a distribution tensor with equilibrium tr(μ) = 3:

f(μ) = L²/(L² - (tr(μ) - 3))  = L²/(L² - tr(μ) + 3)

The stress becomes: σ = G₀·f(μ)·(μ - I) + p·I

At equilibrium: f(tr=3) = L²/L² = 1 (recovers linear stress). As tr(μ) → L² + 3: f → ∞ (chains reach maximum extensibility).

RheoJAX Implementation

def vlb_fene_factor(mu_xx, mu_yy, mu_zz, L_max):
    tr_mu = mu_xx + mu_yy + mu_zz
    L2 = L_max * L_max
    return L2 / jnp.maximum(L2 - tr_mu + 3.0, 1e-10)

def vlb_stress_fene_xy(mu_xx, mu_yy, mu_zz, mu_xy, G0, L_max):
    f = vlb_fene_factor(mu_xx, mu_yy, mu_zz, L_max)
    return G0 * f * mu_xy

VERDICT: CORRECT. The offset of +3 in the denominator correctly accounts for tr(I) = 3.

---

Summary Table

#	Equation	Paper Reference	RheoJAX	Status
1	Distribution tensor evolution	JMPS 2018 Eq.17, Polymers 2018 Eq.3b	_kernels.py:127-131	CORRECT
2	Cauchy stress	JMPS 2018 Eq.28 (N→∞)	_kernels.py:176-216	CORRECT
3	Free energy	Polymers 2018 Eq.5	—	NOT IMPLEMENTED
4	Dissipation	JMPS 2018 Eqs.19-21	—	NOT IMPLEMENTED
5	Simple shear components	Derived from Eq.17	_kernels.py:127-131	CORRECT
6	σ_xy and N1	Gaussian limit of Eq.28	_kernels.py:176-216	CORRECT
7	Steady-state shear	Polymers 2018 §3.2	_kernels.py:225-250	CORRECT
8	SAOS G’(ω), G’’(ω)	Maxwell limit	_kernels.py:387-415	CORRECT
9	Stress relaxation G(t)	Maxwell limit	_kernels.py:357-378	CORRECT
10	Creep compliance J(t)	Maxwell element	_kernels.py:462-483	CORRECT
11	Bell breakage	Bell (1978) + VLB	_kernels.py:1055-1081	CORRECT
12	FENE-P factor	FENE-P closure	_kernels.py:1090-1115	CORRECT

Issues Found

1. Docstring inconsistency (minor): Line 17 of _kernels.py states dmu/dt = k_d*(I - mu) + D·mu + mu·D using the symmetric part D, but the actual model uses the full velocity gradient L·mu + mu·L^T. The code is correct; only the docstring needs updating.

2. Free energy and dissipation not implemented: These are thermodynamic quantities useful for consistency checks but not required for mechanical predictions. Consider adding for completeness.

---

Key Clarification: L vs D in the Evolution Equation

The VLB paper uses the full velocity gradient L = ∇v, not just its symmetric part D = (L + L^T)/2. This is critical because:

• With L·μ + μ·L^T: the model is equivalent to the upper-convected Maxwell model

• With D·μ + μ·D: the model would be the corotational form

The implementation correctly uses L·μ + μ·L^T (verified by the simple shear components where dμ_yy/dt has NO deformation contribution, which is only true with the full L).