HVNM Advanced Theory & Numerical Methods

This page documents the thermodynamic foundations, interphase physics, damage mechanics, and numerical methods underlying the HVNM. For the constitutive equations, see HVNMLocal — Full Model Reference. For protocol derivations, see HVNM Protocol Equations & Derivations. For the unfilled HVM theory, see HVM Advanced Theory & Numerical Methods.

Thermodynamic Framework

Helmholtz Free Energy

The total Helmholtz free energy density extends the HVM (Thermodynamic Framework) with the interphase contribution and interfacial damage:

\[\Psi_{tot} = (1-D)\,\Psi_P(\mathbf{B}_{amp}) + \Psi_E[\boldsymbol{\mu}^E, \boldsymbol{\mu}^E_{nat}] + \Psi_D[\boldsymbol{\mu}^D] + (1-D_{int})\,\Psi_I[\boldsymbol{\mu}^I, \boldsymbol{\mu}^I_{nat}] + p(\det\mathbf{F} - 1)\]

Each subnetwork energy is a Gaussian chain model (neo-Hookean):

Permanent network (hydrodynamically amplified):

\[\Psi_P = \frac{G_P X(\phi)}{2}\left(\text{tr}(\mathbf{B}) - 3\right)\]

where \(X(\phi) = 1 + 2.5\phi + 14.1\phi^2\) is the Guth-Gold strain amplification factor.

Exchangeable (vitrimer) network:

\[\Psi_E = \frac{G_E}{2}\,\text{tr}\!\left(\boldsymbol{\mu}^E - \boldsymbol{\mu}^E_{nat}\right)\]

Dissociative (physical) network:

\[\Psi_D = \frac{G_D}{2}\,\text{tr}(\boldsymbol{\mu}^D - \mathbf{I})\]

Interphase network (NEW):

\[\Psi_I = \frac{G_{I,eff}}{2}\,\text{tr}\!\left(\boldsymbol{\mu}^I - \boldsymbol{\mu}^I_{nat}\right)\]

The new scalar damage variable \(D_{int} \in [0,1]\) couples to \(\Psi_I\), representing interfacial debonding. \(D\) and \(D_{int}\) are independent: permanent network scission and interfacial failure are distinct mechanisms.

Clausius-Duhem Derivation

Applying the Clausius-Duhem procedure to \(\Psi_{tot}\), the total Cauchy stress is:

\[\boldsymbol{\sigma}_{tot} = \underbrace{(1-D)\,\tilde{G}_P\,(\mathbf{B}_{amp} - \mathbf{I})}_{\text{Permanent}} + \underbrace{G_E (\boldsymbol{\mu}^E - \boldsymbol{\mu}^E_{nat})}_{\text{Exchangeable}} + \underbrace{G_D (\boldsymbol{\mu}^D - \mathbf{I})}_{\text{Dissociative}} + \underbrace{(1-D_{int})\,G_{I,eff} (\boldsymbol{\mu}^I - \boldsymbol{\mu}^I_{nat})}_{\text{Interphase}} - p\mathbf{I}\]

where \(\tilde{G}_P = G_P X(\phi)\) is the amplified permanent modulus.

The remaining terms yield four dissipation contributions, each individually non-negative:

Exchangeable network dissipation:

\[\mathcal{D}_{exch} = \frac{G_E}{2} k_{BER}^{mat} \text{tr}\!\left[(\boldsymbol{\mu}^E - \boldsymbol{\mu}^E_{nat})^2 \cdot (\boldsymbol{\mu}^E_{nat})^{-1}\right] \geq 0\]

Dissociative network dissipation:

\[\mathcal{D}_{diss} = \frac{G_D}{2} k_d^D \text{tr}(\boldsymbol{\mu}^D - \mathbf{I})^2 \geq 0\]

Interphase network dissipation (NEW):

\[\mathcal{D}_{int} = \frac{G_{I,eff}}{2} k_{BER}^{int} \text{tr}\!\left[(\boldsymbol{\mu}^I - \boldsymbol{\mu}^I_{nat})^2 \cdot (\boldsymbol{\mu}^I_{nat})^{-1}\right] \geq 0\]

Damage dissipation (dual):

\[\mathcal{D}_{dam} = \Psi_P \dot{D} + \Psi_I \dot{D}_{int} \geq 0\]

Satisfied because \(\Psi_P, \Psi_I \geq 0\) and \(\dot{D}, \dot{D}_{int} \geq 0\) (damage is irreversible or controlled by self-healing).

Nanoparticle Interphase Model

Three-Layer Interphase Structure

Following Papon et al. (2012) and confirmed by NMR and scattering studies (Berriot et al. 2002, Kim et al. 2024), the interphase around each nanoparticle consists of three concentric layers:

Glassy layer (thickness \(\delta_g \sim 1\)–2 nm): Chains strongly adsorbed/bonded to the NP surface with dynamics \(>100\times\) slower than bulk. In the HVNM, this layer is treated as part of the rigid NP inclusion, increasing the effective filler fraction:

\[\phi_{eff} = \phi\left(1 + \frac{\delta_g}{R_{NP}}\right)^3\]

Mobile interphase layer (thickness \(\delta_m \sim 5\)–20 nm): Chains with intermediate dynamics – slower than bulk but not glassy. This is the interphase subnetwork (I) in the HVNM, with volume fraction:

\[\phi_I = \phi_{eff}\left[\left(1 + \frac{\delta_m}{R_{NP} + \delta_g} \right)^3 - 1\right]\]

For small \(\delta_m / R_{NP}\), this reduces to \(\phi_I \approx 3\phi_{eff}\,\delta_m / (R_{NP} + \delta_g)\).

Bulk matrix (remainder): Unperturbed dynamics, described by the P, E, D subnetworks.

Note

The RheoJAX implementation uses the simplified formula \(\phi_I = \phi[(R_{NP} + \delta_m)^3/R_{NP}^3 - 1]\), which absorbs \(\delta_g\) into \(R_{NP}\).

Temperature-Dependent Interphase Thickness

The mobile interphase thickness decreases with temperature following a WLF-type dependence:

\[\delta_m(T) = \delta_m^0 \cdot \left( \frac{T_g^{int} - T_\infty}{T - T_\infty}\right)^{1/\alpha_{int}}\]

where \(T_g^{int}\) is the elevated glass transition in the interphase, \(T_\infty = T_g^{int} - C_2^{int}\) is the Vogel temperature, and \(\alpha_{int} \sim 0.5\)–1 is a scaling exponent.

Physical consequence: For \(T \gg T_g^{int}\), \(\delta_m \to 0\) and the interphase contribution vanishes, recovering the unfilled HVM. Near \(T_g^{int}\), the interphase is thick and NP reinforcement is maximized.

Interphase Modulus

The interphase modulus is expressed as:

\[G_I = \beta_I \cdot G_E\]

where \(\beta_I \in [1.5, 10]\) is the reinforcement ratio, dependent on:

• Surface chemistry: Covalently grafted NPs (\(\beta_I \sim 3\)–10) vs. physically adsorbed (\(\beta_I \sim 1.5\)–3)

• Chain-NP interaction strength: Quantified by the Flory-Huggins interaction parameter \(\chi_{NP}\)

• Temperature: \(\beta_I\) decreases toward 1 as \(T\) increases well above \(T_g^{int}\)

The effective interphase modulus is:

\[G_{I,eff} = \beta_I G_E \phi_I\]

Interphase Percolation

When \(\phi_I\) exceeds the percolation threshold \(\phi_I^{perc} \approx 0.15\)–0.30, a connected interphase network forms spanning the entire sample. This produces dramatic stiffening and slowing of dynamics. The percolation-enhanced modulus is:

\[G_I^{eff}(\phi_I) = G_I \cdot \phi_I \cdot \left[1 + \kappa \left(\frac{\phi_I - \phi_I^{perc}}{\phi_I^{perc}}\right)^+\right]\]

where \(\kappa\) is the percolation enhancement factor and \((\cdot)^+\) denotes the Macaulay bracket. Below percolation (\(\phi_I < \phi_I^{perc}\)), the interphase acts as isolated shells around NPs. Above percolation, the connected network adds geometric stiffening.

Strain Amplification in the Interphase

The interphase chains experience amplified strain because the rigid NP cores do not deform. The effective velocity gradient experienced by the interphase is:

\[\mathbf{L}^I = X_I(\phi)\,\mathbf{L}\]

where \(X_I = X(\phi_{eff})\) uses the effective NP volume fraction (including glassy layer). This amplification appears in the affine terms of the interphase evolution equation:

\[\dot{\boldsymbol{\mu}}^I = X_I\,(\mathbf{L}\boldsymbol{\mu}^I + \boldsymbol{\mu}^I\mathbf{L}^T) + k_{BER}^{int}\,(\boldsymbol{\mu}^I_{nat} - \boldsymbol{\mu}^I)\]

Physical consequence: Strain amplification means interphase chains reach large stretch (and damage threshold) earlier than bulk matrix chains at the same macroscopic strain. This explains the early onset of nonlinearity in filled systems.

Non-Affine Interphase Dynamics

An alternative to strain amplification introduces a monomer-particle friction parameter \(\xi_{NP}\) that produces non-affine, partially suppressed deformation:

\[\dot{\boldsymbol{\mu}}^I = (1 - \xi_{NP})\, (\mathbf{L}\boldsymbol{\mu}^I + \boldsymbol{\mu}^I\mathbf{L}^T) + k_{BER}^{int}\,(\boldsymbol{\mu}^I_{nat} - \boldsymbol{\mu}^I)\]

where \(\xi_{NP} \in [0, 1]\):

• \(\xi_{NP} = 0\): Full affine deformation (all strain transmitted to interphase)

• \(\xi_{NP} = 1\): Completely pinned layer (no deformation)

The two approaches (strain amplification and friction) are complementary: \(X_I > 1\) models the geometric concentration of strain; \(\xi_{NP}\) models the dynamical suppression. They can be combined as \((1 - \xi_{NP}) X_I\), but this introduces a degeneracy. The recommended default is to use \(X_I\) alone.

Enhanced Damage Mechanics

Two Damage Variables

The HVNM carries two independent scalar damage variables:

Permanent network damage \(D\) (chain scission, identical to HVM):

\[\dot{D} = \Gamma_0\,\bigl(\lambda_{eff}^{perm} - \lambda_{crit}\bigr)^+ \cdot (1 - D)\]

Interfacial damage \(D_{int}\) (debonding/desorption/interfacial bond rupture):

\[\dot{D}_{int} = \Gamma_0^{int}\, \bigl(\lambda_{chain}^{int} - \lambda_{crit}^{int}\bigr)^+ \cdot (1 - D_{int}) - h_{int}(T)\,(D_{int})^{n_h}\]

where the interfacial chain stretch is:

\[\lambda_{chain}^{int} = \sqrt{\frac{\text{tr}(\boldsymbol{\mu}^I)}{3}}\]

Key differences from permanent damage:

• Lower critical stretch \(\lambda_{crit}^{int} < \lambda_{crit}\): Confined interphase chains have less extensibility, and the interface concentrates stress due to modulus mismatch.

• Self-healing term \(-h_{int}(T)\,(D_{int})^{n_h}\): Above \(T_v^{int}\), interfacial BER can reform broken bonds, reducing \(D_{int}\) over time. The healing rate follows TST:

  \[h_{int}(T) = h_0 \exp\!\left(-\frac{E_a^{heal}}{k_B T}\right)\]

  The exponent \(n_h \in [0.5, 1]\) controls healing kinetics shape (\(n_h = 1\) for first-order, \(n_h = 0.5\) for diffusion-limited).

This makes interfacial damage reversible above \(T_v^{int}\) but irreversible below it – the hallmark of vitrimer nanocomposites.

Weissenberg Number Fracture Criterion

The competition between loading rate and BER rate determines whether the material flows or fractures, captured by dual Weissenberg numbers:

\[\text{Wi}^{mat} = \frac{\dot{\epsilon}}{k_{BER}^{mat}}, \qquad \text{Wi}^{int} = \frac{\dot{\epsilon}}{k_{BER}^{int}}\]

Three fracture regimes:

Regime	Condition	Behavior
I: Flow	Wimat \(\ll 1\), Wiint \(\ll 1\)	All networks relax faster than load accumulates. Viscous flow, no damage.
II: Partial	Wimat \(\ll 1\), Wiint \(\gtrsim 1\)	Matrix flows but interphase cannot relax. \(D_{int}\) grows. Ductile.
III: Brittle	Wimat \(\gg 1\), Wiint \(\gg 1\)	Neither network relaxes. Both \(D\) and \(D_{int}\) grow. Thermoset-like.

Crack tip analysis: Near a propagating crack tip, the local strain rate diverges as \(\dot{\epsilon} \sim v_{crack}/r\). The fracture process zone size is:

\[r_{fz}^{mat} = \frac{v_{crack}}{k_{BER}^{mat}}, \qquad r_{fz}^{int} = \frac{v_{crack}}{k_{BER}^{int}}\]

Since \(k_{BER}^{int} < k_{BER}^{mat}\), the interphase fracture zone is larger (\(r_{fz}^{int} > r_{fz}^{mat}\)), meaning interfacial damage extends further ahead of the crack tip. This dissipates energy and is the micromechanical origin of NP toughening in vitrimer nanocomposites.

Cooperative Shielding (4-Network)

The HVM’s cooperative shielding concept is extended to include the interphase. The effective permanent chain stretch is:

\[\lambda_{eff}^{perm} = \lambda_{chain}^{perm} \cdot \frac{(1-D)\,\tilde{G}_P} {(1-D)\,\tilde{G}_P + G_E(t) + G_D(t) + (1-D_{int})\,G_{I,eff}(t)}\]

where \(G_\alpha(t)\) represents each network’s current load share.

NP toughening mechanism: The interphase carries load that would otherwise stress the permanent network, delaying \(D > 0\). When \(D_{int}\) grows (interface fails), shielding decreases and load transfers to the permanent network, causing delayed permanent damage – a cascading failure that produces gradual, ductile failure.

Payne Effect

The Payne effect (strain-amplitude-dependent modulus drop in filled rubbers) emerges naturally from the HVNM:

1. Small \(\gamma_0\): Interphase intact (\(D_{int} \approx 0\)), full modulus \(G' = G_P X + G_E + G_D + G_{I,eff} X_I\)

2. Increasing \(\gamma_0\): Interphase natural state begins tracking deformation via BER, reducing \(\sigma_I\)

3. Large \(\gamma_0\): \(\sigma_I \to 0\) at steady state, \(G'\) drops to unfilled level

The critical strain for Payne onset is:

\[\gamma_c \approx \frac{\lambda_{crit}^{int} - 1}{X_I(\phi)}\]

Higher \(\phi\) lowers the onset strain through the \(X_I\) amplification factor. See LAOS Derivation for detailed LAOS analysis.

Mullins Effect

Under cyclic loading, the HVNM predicts stress softening through three mechanisms:

• Irreversible (from \(D\)): Permanent chain scission softens the elastic response permanently.

• Partially reversible (from \(D_{int}\)): Interfacial damage reduces interphase stress, but self-healing above \(T_v^{int}\) partially restores it between cycles.

• Fully reversible (from BER): Bond exchange in transient networks relaxes stress during each cycle.

The Mullins effect in vitrimers is temperature-dependent: above \(T_v^{int}\), softening is partially recovered between cycles (due to self-healing); below \(T_v^{int}\), it is permanent. See Cyclic Loading & Mullins Effect for the full cyclic loading analysis.

Diffusion-Limited Slow Mode

Karim, Vernerey & Sain (Macromolecules, 2025) identified that the long-term mechanical response of vitrimers requires a constant-rate kinetic term representing diffusion-driven chain dynamics slower than the BER timescale. In a nanocomposite, this is physically motivated by:

• Reptation-like chain diffusion through the network (\(\tau_{diff} \gg \tau_{BER}\))

• Slow chain extraction from the interphase as chains desorb from NP surfaces

• Long-range topological reorganization requiring multiple sequential BER events

The HVNM adds a constant background rate \(k_{diff}\) to both kinetics:

\[k_{eff}^{mat} = k_{BER}^{mat}(T,\boldsymbol{\sigma}^E) + k_{diff}^{mat}(T)\]

\[k_{eff}^{int} = k_{BER}^{int}(T,\boldsymbol{\sigma}^I) + k_{diff}^{int}(T)\]

The diffusion rates follow simple Arrhenius:

\[k_{diff}^{mat} = k_{diff,0}^{mat}\exp\!\left( -\frac{E_a^{diff}}{k_B T}\right), \qquad k_{diff}^{int} = k_{diff,0}^{int}\exp\!\left( -\frac{E_a^{diff,int}}{k_B T}\right)\]

Key properties:

• The diffusion activation energy \(E_a^{diff}\) is typically 1.5–3\(\times\) larger than \(E_a^{mat}\), reflecting the higher barrier for coordinated multi-bond rearrangement.

• The diffusion rate is stress-independent (not TST-activated) because it represents thermally driven Brownian motion.

• \(k_{diff}^{int} \ll k_{diff}^{mat}\) because chain extraction from NP surfaces is entropically penalized.

When to enable: The include_diffusion=True flag should be used when stress relaxation data shows a long-time tail beyond the BER relaxation, or when creep data shows continuing slow deformation at \(t \gg 1/k_{BER}\).

Numerical Implementation

ODE state vector: 17 components in simple shear (18 with interfacial damage):

Index	Component	Description
0–2	\(\mu^E_{xx}, \mu^E_{xy}, \mu^E_{yy}\)	Exchangeable distribution tensor
3–5	\(\mu^{E,nat}_{xx}, \mu^{E,nat}_{xy}, \mu^{E,nat}_{yy}\)	Exchangeable natural state
6–8	\(\mu^D_{xx}, \mu^D_{xy}, \mu^D_{yy}\)	Dissociative distribution tensor
9–11	\(\mu^I_{xx}, \mu^I_{xy}, \mu^I_{yy}\)	Interphase distribution tensor
12–14	\(\mu^{I,nat}_{xx}, \mu^{I,nat}_{xy}, \mu^{I,nat}_{yy}\)	Interphase natural state
15	\(\gamma\)	Accumulated strain
16	\(D\)	Permanent network damage
17	\(D_{int}\)	Interfacial damage (if include_interfacial_damage=True)

ODE solver: diffrax Tsit5 (explicit 5th-order Runge-Kutta) with PIDController adaptive stepping (rtol=1e-8, atol=1e-10).

Note

Implicit solvers (e.g., Kvaerno5) were tested but produce TracerBoolConversionError due to lineax LU transpose checks during JAX tracing. Tsit5 is the recommended solver.

Stiffness at high \(\phi\): High \(\phi\) amplifies the affine deformation (\(X_I \dot{\gamma}\)), creating stiff ODEs. If the solver diverges, increase max_steps or reduce the shear rate.

Square-root guard: The BER rate computation uses:

safe_stretch = jnp.sqrt(jnp.maximum(stretch_invariant, 0.0) + 1e-30)

This prevents infinite gradients at \(\sigma_{VM} = 0\).

Initial conditions: All tensors at identity (\(\mu_{xx} = \mu_{yy} = 1\), \(\mu_{xy} = 0\)), \(\gamma = 0\), \(D = D_{int} = 0\).

\(\phi = 0\) Recovery Verification

When \(\phi = 0\), the HVNM must recover the HVM exactly. Mathematically:

1. \(X(\phi=0) = 1\) and \(X_I = X(\phi_{eff}=0) = 1\): No strain amplification.

2. \(\phi_I = 0 \cdot [(\ldots)^3 - 1] = 0\): Zero interphase volume fraction.

3. \(G_{I,eff} = \beta_I G_E \cdot 0 = 0\): No interphase modulus contribution.

4. The \(D_{int}\) equation decouples (no interphase stress to drive damage).

5. The \(\boldsymbol{\mu}^I\) equations decouple (zero prefactor in stress).

The remaining equations are identical to HVM with \(D = 0\), \(X = 1\). This is verified numerically to machine precision (relative error \(< 10^{-14}\)) across all six protocols in the test suite.

References

1. Vernerey, F.J., Long, R. & Brighenti, R. (2017). “A statistically-based continuum theory for polymers with transient networks.” J. Mech. Phys. Solids, 107, 1–20. https://doi.org/10.1016/j.jmps.2017.05.016

2. 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

3. Vernerey, F.J., Brighenti, R., Long, R. & Shen, T. (2018). “Statistical Damage Mechanics of Polymer Networks.” Macromolecules, 51(17), 6609–6622. https://doi.org/10.1021/acs.macromol.8b01052

4. Meng, F., Saed, M.O. & Terentjev, E.M. (2019). “Elasticity and Relaxation in Full and Partial Vitrimer Networks.” Macromolecules, 52(19), 7423–7429. https://doi.org/10.1021/acs.macromol.9b01123

5. Shen, T., Song, Z., Cai, S. & Vernerey, F.J. (2021). “Nonsteady fracture of transient networks: The case of vitrimer.” PNAS, 118(29), e2105974118. https://doi.org/10.1073/pnas.2105974118

6. Song, Z., Wang, Z. & Cai, S. (2021). “Mechanics of vitrimer with hybrid networks.” Mech. Mater., 153, 103687. https://doi.org/10.1016/j.mechmat.2020.103687

7. Papon, A., Montes, H., Lequeux, F. et al. (2012). “Glass-transition temperature gradient in nanocomposites: Evidence from nuclear magnetic resonance and differential scanning calorimetry.” Phys. Rev. Lett., 108, 065702. https://doi.org/10.1103/PhysRevLett.108.065702

8. Berriot, J., Montes, H., Lequeux, F. et al. (2002). “Filler-elastomer interaction in model filled rubbers, a ¹H NMR study.” J. Non-Crystalline Solids, 307–310, 719–724. https://doi.org/10.1016/S0022-3093(02)01552-1

9. Duan, P., Zhao, H., Chen, Q. et al. (2023). “Insights into Uniaxial Tension and Relaxation of Nanorod-Filled Polymer Vitrimer Nanocomposites: A Molecular Dynamics Simulation.” Macromolecules, 56(12), 4468–4481. https://doi.org/10.1021/acs.macromol.3c00296

10. Li, Z., Zhao, H., Duan, P., Zhang, L. & Liu, J. (2024). “Manipulating the Properties of Polymer Vitrimer Nanocomposites by Designing Dual Dynamic Covalent Bonds.” Langmuir, 40(14), 7769–7780. https://doi.org/10.1021/acs.langmuir.4c00699

11. Duan, P., Zhao, H., Liu, M. et al. (2024). “Molecular Insights into the Topological Transition, Fracture, and Self-Healing Behavior of Vitrimer Composites with Exchangeable Interfaces.” Macromolecules, 57, 9725–9736. https://doi.org/10.1021/acs.macromol.4c01541

12. Kim, J., Thompson, B.R., Tominaga, T. et al. (2024). “Suppression of Segmental Chain Dynamics on a Particle’s Surface in Well-Dispersed Polymer Nanocomposites.” ACS Macro Lett., 13(6), 720–725. https://doi.org/10.1021/acsmacrolett.4c00168

13. Karim, M.R., Vernerey, F. & Sain, T. (2025). “Constitutive Modeling of Vitrimers and Their Nanocomposites Based on Transient Network Theory.” Macromolecules, 58(10), 4899–4912. DOI: 10.1021/acs.macromol.4c02872 PDF

14. Li, Z., Zhao, H., Zhang, L., Liu, J. et al. (2025). “Designing All-Vitrimer Nanocomposites to Combine Low Energy Consumption, Mechanical Robust and Recyclability.” Nano Energy, 142, 111215.

15. Alkhoury, K. & Chester, S.A. (2025). “A chemo-thermo-mechanically coupled theory of photo-reacting polymers: Application to modeling photo-degradation with irradiation-driven heat transfer.” J. Mech. Phys. Solids, 197, 106050. https://doi.org/10.1016/j.jmps.2025.106050

16. Wagner, R.J. & Silberstein, M.N. (2025). “A foundational framework for the mesoscale modeling of dynamic elastomers and gels.” J. Mech. Phys. Solids, 194, 105914. https://doi.org/10.1016/j.jmps.2024.105914

17. Hayashi, M. & Ricarte, R.G. (2025). “Towards the next development of vitrimers: Recent key topics for the practical application and understanding of the fundamental physics.” Prog. Polym. Sci., 170, 102026. DOI: 10.1016/j.progpolymsci.2025.102026 PDF

18. Zhao, H., Wei, X., Fang, Y. et al. (2022). “Molecular dynamics simulation of the structural, mechanical, and reprocessing properties of vitrimers based on a dynamic covalent polymer network.” Macromolecules, 55, 1091–1103. https://doi.org/10.1021/acs.macromol.1c02034