dirt_bond

Bonded Particle Model (BPM) forces for DIRT — elastic bonds, breakage, and plasticity.

What it does

dirt_bond adds the DemBondPlugin, which treats each bond between two particles as a solid cylindrical beam resisting four independent deformation channels: axial stretch/compression, transverse shear, twist (torsion about the bond axis), and bending (rotation perpendicular to the axis). Stiffness can be derived from material properties (Young's / shear modulus + bond radius) or supplied directly; damping is set as a critical-damping ratio; bonds can break under a configurable failure criterion and yield under configurable plasticity models.

Bonds are created at setup either by auto-bonding initially touching pairs or by loading a LAMMPS Bonds section. Periodic minimum-image distances are used so closed-loop / wrapping bonds get the correct rest length. Across MPI ranks, ghost_cutoff is automatically extended to keep bonded partners (and their 1-3 shared neighbours) visible as ghosts. Granular contact between bonded pairs is suppressed via soil_core's BondStore, so the bond is the sole interaction until it breaks.

Bond geometry

Per bond the cross-section is a solid cylinder of radius r_b = bond_radius_ratio · min(R_i, R_j), giving area A = π r_b², polar second moment J = ½ π r_b⁴ (torsion), bending second moment I = ¼ π r_b⁴ = ½ J, and length L = r₀.

Four-channel force model

The shear force acts at the bond mid-point, producing a lever-arm torque τ = (L/2) n̂ × F_t on both particles. Per-channel damping uses a critical-damping ratio β ∈ [0, 1] (γ = 2 β √(m*·K) for forces, 2 β √(I*·K) for moments) with reduced mass / reduced MOI of the pair; each channel accepts an optional raw-γ override.

Breakage

A bond breaks when the active criterion (breakage::BreakageCriterion) trips. The menu (breakage::BreakageConfig) covers three families × force/stress/strain measures plus an Unbreakable default:

• axial_force / axial_stress / axial_strain — independent tensile + shear branches
• combined_stress / combined_strain — extreme-fibre beam criterion (Potyondy–Cundall / Guo), tensile and shear branches
• interaction_linear_force / _stress / _strain — single linear damage envelope Σ |X_i| / X_{i,c} ≥ 1 over the four channels (LAMMPS bpm/rotational)

Each threshold is a ThresholdDistribution (kind = …), one of two:

• constant — the same value for every bond.
• weibull — length-scaled 2-parameter Weibull (mean, m, l_calib, optional l_min) giving the weakest-link size effect.

Per-bond thresholds are drawn once at creation from the bond's tag pair plus seed, so the breakage pattern is MPI-decomposition independent and bit-reproducible.

Plasticity

The bending and axial channels are independently optional (plasticity::PlasticityConfig); omit a channel to keep it purely elastic.

• bending: guo_bending (elastic–perfectly-plastic cap at M^p = (4/3) σ_0 r_b³), guo_trilinear (Guo 2018 Eq. 32: elastic → K_ep = K_e/2 elasto-plastic → perfectly-plastic trilinear envelope; expands at setup to a two-breakpoint piecewise and needs [bonds].youngs_modulus), or piecewise
• axial: piecewise (breakpoints in axial strain, slope multipliers)

Plastic anchors and max-strain history are carried per bond in BondHistoryStore and survive MPI communication and atom reordering.

Key types

Configuration

[bonds]
auto_bond = true
bond_tolerance = 1.001
bond_radius_ratio = 1.0
ghost_cutoff_multiplier = 2.5   # MPI: ghost skin to cover bond + 1-3 reach

# Material mode (paper-standard beam theory):
youngs_modulus = 1.0e9     # E (Pa) → K_n, K_bend
shear_modulus  = 4.0e8     # G (Pa) → K_t, K_tor
# Direct overrides used when the matching modulus is absent:
# normal_stiffness / shear_stiffness / twist_stiffness / bending_stiffness

# Critical-damping ratios (raw *_damping overrides also available):
beta_normal = 0.05
beta_shear  = 0.05
beta_twist  = 0.05
beta_bending = 0.05

seed = 0   # per-bond threshold RNG seed

[bonds.breakage]
kind    = "combined_stress"
tensile = { kind = "weibull", mean = 5.0e7, m = 8.0, l_calib = 0.020 }
shear   = { kind = "constant", value = 3.0e7 }
# Threshold kinds: "constant" or "weibull".

[bonds.plasticity.bending]
kind         = "guo_bending"
yield_stress = 1.23e8
# Or the Guo 2018 trilinear envelope (needs [bonds].youngs_modulus):
# kind         = "guo_trilinear"
# yield_stress = 1.23e8

[bonds.plasticity.axial]
kind               = "piecewise"
breakpoint_strains = [0.01, 0.02, 0.03]
slope_multipliers  = [0.5, 0.1, 0.0]

# Load explicit bonds instead of auto-bonding:
# file = "bonds.lammps"
# format = "lammps_data"

Usage

This crate exports DemBondPlugin (and BondConfig, the breakage/plasticity menus, the history stores). It depends only on dirt_atom, soil_core, and the GRASS app/scheduler — not on dirt_core. A full DEM application is assembled through the dirt_core umbrella crate, which re-exports DemBondPlugin from its prelude alongside the core/granular plugin groups:

use dirt_core::prelude::*;          // re-exports DemBondPlugin, CorePlugins, GranularDefaultPlugins

let mut app = App::new();
app.add_plugins(CorePlugins)        // atom data, neighbor list, dump/thermo I/O
    .add_plugins(GranularDefaultPlugins) // Hertz-Mindlin contact + Verlet integration
    .add_plugins(DemBondPlugin);    // bond forces (suppresses contact on bonded pairs)
app.start();

If you depend on dirt_bond directly (without the umbrella), import the plugin from this crate instead — use dirt_bond::DemBondPlugin; — and supply your own core/granular plugins. DemBondPlugin needs the granular contact plugin present so non-bonded pairs still interact, and contact on bonded pairs is suppressed via soil_core's BondStore.

Thermo output

BondMetrics publishes three keys each thermo step:

Concepts

Bond vs. contact. A bonded pair is registered in soil_core's BondStore, and the granular contact plugin skips any pair found there — so a bond is the sole interaction between its two particles until it breaks. On breakage the pair is removed from BondStore and normal Hertz/Hooke contact resumes automatically (relevant for overlapping bonded spheres, which would otherwise double-count force).

Sign convention. Each bond is owned by the lower-tag atom. Forces are applied as +F on atom i (lower tag) and −F on atom j (higher tag) — Newton's third law by construction. Moments follow the same ownership.

Per-bond history lifecycle. A bond is created at setup (auto-bond of touching pairs, or a loaded Bonds section). At creation its breakage thresholds are drawn once from (tag_i, tag_j, seed), so the failure pattern is decomposition-independent and bit-reproducible. Each step the shear Δs and twist/bend Δθ accumulators (and any plastic anchors / max-strain history) are carried in BondHistoryStore, which survives MPI communication and atom reordering. When the active criterion trips, the bond is removed.

bond_type is parsed and stored on each bond entry but is not currently consumed by the force model — it is reserved for future per-type bond parameter sets.

Examples

In the repo's examples/: bond_fiber_tensile (recovers the input Young's modulus from a pulled fiber), bond_fiber_tensile_overlap (contact-suppression check on overlapping bonded spheres), bond_cantilever, and bond_mpi_drift.

cargo run --release --example bond_fiber_tensile --no-default-features -- \
    examples/bond_fiber_tensile/config.toml

License

MIT OR Apache-2.0