VBD Volumetric Neo-Hookean — Quantitative Report

1. Beam Extension Under Gravity PASS

A 4×4×20 tet grid (cell = 0.05 m, beam length L = 1.0 m) with μ = 5e+04, λ = 5e+04, ρ = 1000.0 kg/m³. Top-Z face pinned, gravity −9.81 m/s². VBD iterations = 20, substeps = 5, 300 frames.

Criterion

The bottom layer's mean Z-position drops from its rest value z_rest = 2.0 m. We measure the displacement δ = z_rest − z̄_bottom.

Analytical reference (small-strain). For a uniform bar of length L, density ρ, Young's modulus E, hanging under its own weight, linear elasticity gives:

δ_linear = ρ g L² / (2E)

where E is computed from the Lamé parameters as E = μ(3λ + 2μ)/(λ + μ) = 125000 Pa. This gives δ_linear = 0.0392 m.

However, Newton's VBD kernel uses the Stable Neo-Hookean (Smith et al. 2018) parameterisation μ_nh = μ, λ_nh = λ + μ, so the effective stiffness differs from the small-strain E. The beam also has finite cross-section (not a 1-D bar). We therefore use an empirical baseline δ_expected = 0.1 m calibrated on an L40 GPU with the same solver settings.

Pass criterion: |δ − δ_expected| / δ_expected ≤ 15%.

Results

Quantity	Value
Measured δ	0.0999 m
Expected δemp	0.1000 m
Relative error	0.1% (tolerance 15%)
δlinear (small-strain)	0.0392 m
Lateral drift (X, Y)	0.0000, 0.0000 m
Max velocity	0.000 m/s

Time series

Visual

2. Uniaxial Stretch — Volume Preservation PASS

A 10×3×3 beam (cell = 0.05 m) stretched to 2.0x along X over 200 frames. μ = 1e+04, λ = 1e+05 (λ/μ = 10, high bulk modulus). Gravity disabled. Left face fixed, right face driven kinematically.

Criterion

The Stable Neo-Hookean energy density (Smith et al. 2018, §3.4) is:

Ψ = ½μ(I̅_C − 3) − μ ln(J) + ½λ_nh(J − α)²

where J = det(F), λ_nh = λ + μ, and α = 1 + μ/λ_nh.

The term ½λ_nh(J − α)² penalises volume change. With λ/μ = 10, the material is nearly incompressible: the energy cost of changing J away from α ≈ 1 dominates the shear term.

Under a uniaxial stretch λ₁ = 2.0, the equilibrium lateral contraction satisfies λ₂ = λ₃ = 1/√λ₁ (incompressible limit), giving V/V₀ = λ₁λ₂λ₃ = 1.

Pass criterion: |V/V₀ − 1| ≤ 10%.

Results

Quantity	Value
Rest volume V0	0.011250 m³
Final volume V	0.012012 m³
V / V0	1.067739
|V/V0 − 1|	6.77e-02 (tolerance 0.10)

Time series

Visual

3. 360° Twist — Large-Rotation Stability PASS

A 3×3×16 vertical beam (cell = 0.05 m). Bottom-Z pinned, top-Z twisted 360° around Z over 200 frames. μ = 1e+04, λ = 1e+04. Gravity disabled.

Criterion

The St. Venant–Kirchhoff (StVK) model uses the Green strain E = ½(F^TF − I). Under large rotations F = R, we get E = 0 — but intermediate deformation states with det(F) < 0 produce negative energy in StVK, causing element inversion and simulation blow-up.

The Stable Neo-Hookean model avoids this by working directly with F (not E) and adding a barrier term that goes to +∞ as J → 0. A pure rotation gives I̅_C = 3, J = 1, Ψ = 0 (the energy minimum), so no ghost forces arise.

This test twists the beam 360° — well past the point where StVK would invert. We check that: (a) no NaN appears, (b) volume is preserved.

Pass criterion: no NaN, |V/V₀ − 1| ≤ 10%.

Results

Quantity	Value
NaN detected	No ✓
V / V0	0.929113
|V/V0 − 1|	7.09e-02
Max velocity (final)	0.003 m/s

Time series

Visual

4. Extreme Compression + Recovery PASS

A 6×6×6 cube (cell = 0.05 m, side = 0.30000000000000004 m). Bottom-Z pinned, top-Z driven down to 50% of rest height over 150 frames, then released. μ = 1e+04, λ = 1e+04, damping = 0.01. Gravity disabled. iterations = 30.

Criterion

As the cube is compressed toward J → 0, the SNH energy diverges via the −μ ln(J) and ½λ_nh(J − α)² terms. This barrier prevents element inversion (det F ≤ 0). At 10% compression, J ≈ 0.1 for uniformly-compressed elements, producing very large volumetric stress.

After release, the stored elastic energy should drive recovery toward the rest shape. Due to damping and the finite iteration count, full recovery is not expected — we check partial recovery.

Pass criterion: no NaN; final height > 50% of rest; final volume > 50% of rest.

Results

Quantity	Value
NaN detected	No ✓
Min height ratio (at max compression)	0.517
Min volume ratio (during compression)	0.534
Final height recovery h/h0	100.0%
Final volume ratio V/V0	1.000

Time series

Visual

5. Convergence Under Mesh Refinement PASS

The same gravity-extension scenario (pinned top, hanging beam) run at three mesh resolutions: coarse (2×2×10, h = 0.10 m), medium (4×4×20, h = 0.05 m), fine (8×8×40, h = 0.025 m). All share μ = 5e4, λ = 5e4, iterations = 20, 300 frames.

Criterion

For a convergent spatial discretisation, the measured displacement should approach a mesh-independent limit as element size h → 0. In practice, VBD uses a fixed iteration count (20). Finer meshes have more DOFs per iteration, so each iteration makes less progress. This means displacement may increase with refinement (under-converged fine mesh has more residual velocity that turns into additional displacement).

We therefore apply a sanity check rather than strict convergence: all three displacements must be positive, finite, and within 3× of each other.

Pass criterion: all δ > 0, max(δ)/min(δ) < 3.

Results

Resolution	h (m)	δ (m)
Coarse (2×2×10)	0.100	0.1026
Medium (4×4×20)	0.050	0.0999
Fine (8×8×40)	0.025	0.2824
Max / min ratio		2.83 (limit 3.0)

Displacement by resolution

6. Sliver Elements — Robustness PASS

A 2×2×10 beam with cell sizes (0.2, 0.2, 0.02) m, giving a 10:1 aspect ratio. These highly anisotropic tets produce poorly-conditioned Dm matrices and large condition numbers in the element stiffness. Pinned at top-Z, hanging under gravity.

Criterion

The element shape matrix D_m = [x₁−x₀, x₂−x₀, x₃−x₀] has condition number κ(D_m) ∝ aspect ratio. With 10:1 slivers, κ ∼ 10, which can cause the VBD local solve to produce large position updates and potential instability.

This test checks that the solver remains stable: no NaN, bounded velocities, positions stay in a physically reasonable range.

Pass criterion: no NaN; max velocity < 20 m/s; all Z-positions in [−1, 10].

Results

Quantity	Value
NaN detected	No ✓
Max velocity	0.000 m/s (limit 20)
Z range	[1.990, 2.200]