Breakage Mechanics

• Breakage mechanics is defined by coupling state variables like energy, stress, and mass density to predict rupture initiation, propagation, and completion across scales.
• It integrates energetic, stochastic, and threshold formulations with computational and constitutive methods to accurately model crack paths, fragment distributions, and failure kinetics.
• The framework spans applications from nanoscale thermal ruptures to macroscopic granular and turbulent systems, providing actionable insights into diverse breakage phenomena.

Breakage mechanics concerns the initiation, propagation, completion, and statistical consequences of rupture events in systems ranging from molecular bonds and nanowires to granular particles, colloidal networks, droplets, and particulate populations. In the literature summarized here, breakage is not restricted to externally driven fracture: it also includes thermally activated rupture in quasi-one-dimensional systems, collision-induced fragmentation in population balances, hydrodynamically activated failure in turbulence and porous media, and constitutive descriptions of distributed grain crushing in geomaterials. The unifying theme is the coupling of a breakage criterion to a state description—free energy, stress, mass density, porosity, daughter distributions, or internal geometry—so that rupture becomes a calculable part of mechanics rather than a purely phenomenological endpoint (Nisoli et al., 2013, Volokh, 2016, Cil et al., 2019, Giri et al., 2020).

1. Conceptual scope and state variables

A central feature of breakage mechanics is that the “object” that breaks depends on the model class. In quasi-one-dimensional thermal rupture, a site variable $\sigma_n=+1$ denotes solid and $\sigma_n=-1$ gas, with breakages located at domain-walls; in particle fragmentation, the relevant field is a sub-particle stress or energy landscape; in brittle-continuum failure, fracture is interpreted as a local material sink through the referential mass density $\rho(x,t)$; in granular constitutive theory, breakage is encoded by a scalar breakage index $B\in[0,1]$; and in population-balance models, the state variable is the particle-size density $f(x,t)$ or $u(t,x)$ together with collision and daughter kernels (Nisoli et al., 2013, Jiang et al., 2020, Volokh, 2016, Cil et al., 2019, Jaiswal et al., 2024).

These formulations differ in representation, but they share two structural requirements. First, a local or mesoscopic condition must determine when rupture is activated: examples include $y(x)\ge \theta$ for energy-landscape damage, $\sigma(t)\ge \sigma_{cr}$ for hydrodynamic activation, $\Delta x_j>\Delta_c$ for spring failure, or $\sigma_t=T_0$ and $\sigma_n=-1$0 for pore-scale bond failure (Jiang et al., 2020, Marchioli et al., 2015, Baek, 2018, Hashemi et al., 2023). Second, breakage must feed back into geometry, stress redistribution, fragment size, or macroscopic response. That feedback is explicit in models that predict fragment meshes, daughter distributions, porosity evolution, or cluster statistics (Jiang et al., 2018, Bariwal et al., 2022, Tripathi et al., 2023).

A plausible implication is that breakage mechanics is best understood not as a single constitutive law but as a family of tightly coupled state-evolution problems. The choice of state variable determines which aspects of rupture become tractable: domain lengths in nanowires, crack paths in individual grains, grading evolution in sands, or scaling profiles in collision-induced fragmentation.

2. Energetic, stochastic, and threshold formulations

Energetic descriptions are prominent across the literature. In quasi-one-dimensional thermal breakage, the rupture cost is encoded by $\sigma_n=-1$1, the chemical-potential bias by $\sigma_n=-1$2, and the free-energy density by $\sigma_n=-1$3, where $\sigma_n=-1$4 is the larger effective fugacity. This framework yields two distinct regimes. In the “ferromagnetic” limit $\sigma_n=-1$5, the crossover width is $\sigma_n=-1$6 near $\sigma_n=-1$7; in the “antiferromagnetic” limit $\sigma_n=-1$8, a plateau appears around $\sigma_n=-1$9 for $\rho(x,t)$0 between $\rho(x,t)$1, corresponding to alternating single-site solid-gas clusters. The same model predicts narrow peaks in fluctuations and specific heat, and at $\rho(x,t)$2 gives $\rho(x,t)$3 for $\rho(x,t)$4, linking rupture probability directly to nanocluster spacing (Nisoli et al., 2013).

At the molecular scale, breakage kinetics can be decomposed into free-energy and friction contributions. For water-bond breakage, the one-dimensional Fokker–Planck description uses the separation coordinate $\rho(x,t)$5, free energy $\rho(x,t)$6, and a radial diffusivity $\rho(x,t)$7 extracted from mean first-passage times. The reported minimum $\rho(x,t)$8 in the first hydration shell represents a nearly six-fold reduction relative to $\rho(x,t)$9, while the free-energy barrier is only of order $B\in[0,1]$0 near $B\in[0,1]$1 nm. Ignoring the spatial variation of $B\in[0,1]$2 would underestimate the mean first-passage times by roughly a factor of two. The interpretation given is that the true reaction path has substantial excursions in orthogonal coordinates, so breakage kinetics cannot be reduced to a scalar barrier alone (Hansen et al., 2011).

Threshold models also distinguish between activation and completion. In turbulent breakage of ductile aggregates, activation begins when $B\in[0,1]$3, equivalently $B\in[0,1]$4, but rupture is completed only when the accumulated energy

$B\in[0,1]$5

reaches $B\in[0,1]$6. The corresponding rate $B\in[0,1]$7 is systematically smaller than the brittle rate $B\in[0,1]$8, and for weak aggregates the brittle law exhibits $B\in[0,1]$9 with $f(x,t)$0, whereas the ductile rate flattens and can plateau as $f(x,t)$1 (Marchioli et al., 2015).

A related probabilistic picture appears in colloidal gels. There, strand breaking is detected through jumps $f(x,t)$2 in chemical distance, and for $f(x,t)$3 the metrics $f(x,t)$4, $f(x,t)$5, and $f(x,t)$6 decrease, $f(x,t)$7 and $f(x,t)$8 increase, and $f(x,t)$9 decreases before rupture. The data are summarized by an Arrhenius-like law $u(t,x)$0, with $u(t,x)$1, indicating exponential sensitivity of breakage to local stress or stored strain energy (Bhaumik et al., 2024).

3. Continuum and constitutive breakage mechanics

A distinctive continuum formulation is Volokh’s coupled mass-momenta balance, which interprets fracture as a local material sink. In referential form, the exact local mass balance is

$u(t,x)$2

with mass flux $u(t,x)$3 and sink

$u(t,x)$4

The resulting two-field equation contains the characteristic length $u(t,x)$5, and the momentum balance retains the standard finite-strain form $u(t,x)$6, with $u(t,x)$7. In this theory the role often played by a damage or phase field is assumed by the mass-density field $u(t,x)$8, so regularization enters through mass diffusion rather than an ad hoc internal variable (Volokh, 2016).

For brittle granular materials, breakage mechanics has been reformulated into a constitutive model that couples breakage and dilation. The internal variable $u(t,x)$9 interpolates between the initial and ultimate particle-size distributions through

$y(x)\ge \theta$0

while porosity $y(x)\ge \theta$1 evolves alongside breakage. The model introduces breakage-dependent porosity bounds $y(x)\ge \theta$2 and $y(x)\ge \theta$3, a breakage-energy conjugate $y(x)\ge \theta$4, and a mixed stress-energy yield function

$y(x)\ge \theta$5

According to the calibration and comparisons reported for Kurnell and Cambria sands, the framework captures strain softening in dense specimens, the reduction of dilatancy and peak strength with increasing confinement, the transition from dilatant to compactive response, and the evolution of particle grading due to distributed breakage events (Cil et al., 2019).

Tripathi et al. extend the constitutive viewpoint by integrating abrasion and breakage in a “particle geometry space” with axes $y(x)\ge \theta$6 and $y(x)\ge \theta$7. In this space, the breakage line, sphere line, and average-shape-conserving line partition the plane into five zones: Zone A, Zone $y(x)\ge \theta$8, Zone $y(x)\ge \theta$9, and the impossible zones $\sigma(t)\ge \sigma_{cr}$0 and $\sigma(t)\ge \sigma_{cr}$1. Pure breakage follows the inherited regression $\sigma(t)\ge \sigma_{cr}$2, pure abrasion tends toward the sphere line $\sigma(t)\ge \sigma_{cr}$3, and equal occurrence follows the slope $\sigma(t)\ge \sigma_{cr}$4 line through $\sigma(t)\ge \sigma_{cr}$5. This framework makes shape, size, area, and volume simultaneous state variables rather than auxiliary descriptors (Tripathi et al., 2023).

At the pore scale in porous media, authigenic-fine detachment by breakage is modeled by combining Timoshenko beam theory, creeping-flow drag, and strength criteria for the particle-rock bond. The explicit breakage velocity

$\sigma(t)\ge \sigma_{cr}$6

is then upscaled into a Maximum Retention Function $\sigma(t)\ge \sigma_{cr}$7, decomposed into DLVO and breakage contributions. Coreflood validation on Castlegate sandstone yielded $\sigma(t)\ge \sigma_{cr}$8 for breakthrough-concentration curves and $\sigma(t)\ge \sigma_{cr}$9 for pressure drop, with bimodal effluent size distributions supporting the two-population interpretation (Hashemi et al., 2023).

4. Particle-scale stress analysis and computational breakage criteria

Particle-resolved computation has become a major route for turning breakage mechanics into a predictive methodology. One approach constructs a scalar energy-release density

$\Delta x_j>\Delta_c$0

derived from the sub-particle stress field $\Delta x_j>\Delta_c$1 obtained by boundary-element analysis. A point is deemed damaged when $\Delta x_j>\Delta_c$2, and ridge points satisfy the Lindeberg conditions $\Delta x_j>\Delta_c$3 and $\Delta x_j>\Delta_c$4. Fragmentation is declared when a connected damaged path separates the particle into subdomains, equivalently when there exists a ridge point $\Delta x_j>\Delta_c$5 with $\Delta x_j>\Delta_c$6. In a 2D Brazilian-disc study with coordination numbers $\Delta x_j>\Delta_c$7, the ridge-based method reproduced the correct number of fragments and approximate crack orientations within $\Delta x_j>\Delta_c$8 in most cases, and predicted critical force within $\Delta x_j>\Delta_c$9 of experiments (Jiang et al., 2020).

A second route is the Boundary-Spheropolygon Element Method, a one-way coupling of a spheropolygon discrete element method with a boundary element method. In BSEM, contact forces obtained from SDEM are mapped as Neumann boundary data for BEM, which then returns a continuous sub-particle stress field $\sigma_t=T_0$0. Tensile fracture initiation is predicted where the maximum principal stress satisfies $\sigma_t=T_0$1, and propagation follows contours of highest $\sigma_t=T_0$2 gradients. Validation against the Brazilian-disc benchmark with 223 linear boundary elements and 200 interior sampling points showed agreement within $\sigma_t=T_0$3 over $\sigma_t=T_0$4, while comparison with ABAQUS indicated that for errors $\sigma_t=T_0$5 BEM is $\sigma_t=T_0$6–$\sigma_t=T_0$7 faster than FEM and requires $\sigma_t=T_0$8–$\sigma_t=T_0$9 fewer degrees of freedom (Jiang et al., 2018).

These methods share a move away from average-stress substitutes toward resolved internal fields. In the BSEM study, the claimed advantages include avoiding agglomerate clusters, volume meshes, and remeshing after fracture; in the energy-landscape method, fragment geometry follows maximal-$\sigma_n=-1$00 ridges rather than an imposed crack template (Jiang et al., 2018, Jiang et al., 2020). This suggests that contemporary breakage mechanics increasingly treats crack initiation, crack path, and fragment shape as outcomes of the same internal field computation.

Dynamic loading can also determine where breakage occurs, even in minimal models. In a homogeneous harmonic chain pulled by a force $\sigma_n=-1$01, the spring extension becomes piecewise linear in the continuum limit, and the first breaking spring is found by minimizing the threshold-crossing time over all springs. For $\sigma_n=-1$02, the chain breaks at the pulled end $\sigma_n=-1$03; for sufficiently small $\sigma_n=-1$04, the earliest failure can occur near the wall $\sigma_n=-1$05. As $\sigma_n=-1$06 varies, the breaking index moves back and forth, the “anomalous” breaking phenomenon explained by competition between intrinsic wave-transit times and the extrinsic loading time $\sigma_n=-1$07 (Baek, 2018).

5. Population-balance models, analysis, and numerical approximation

At the ensemble level, breakage mechanics is formulated through nonlinear integro-differential equations for the particle-size density. In one version, the collision-induced breakage equation combines gain and loss terms associated with binary collisions, daughter production, and secondary destruction, with mass conservation in each breakage event enforced by

$\sigma_n=-1$08

Another equivalent weak formulation uses

$\sigma_n=-1$09

together with the symmetric power-law collision kernel $\sigma_n=-1$10 (Bariwal et al., 2024, Giri et al., 2020).

The rigorous theory shows a sharp dependence on kernel homogeneity. When $\sigma_n=-1$11, there exists at least one global-in-time weak mass-conserving solution for initial data in $\sigma_n=-1$12, and higher moments remain bounded if initially finite. When $\sigma_n=-1$13 with $\sigma_n=-1$14, only local-in-time mass-conserving solutions are obtained, up to

$\sigma_n=-1$15

and, under an additional positivity assumption on fragmentation, no global mass-conserving solution exists. For $\sigma_n=-1$16 and Beta-type daughter distributions $\sigma_n=-1$17, the non-existence result is even stronger: no nontrivial mass-conserving weak solution exists, not even locally in time (Giri et al., 2020).

Self-similar scaling complements this existence theory. Under homogeneous kernels and breakage functions, a mass-conserving self-similar solution takes the form

$\sigma_n=-1$18

so the typical particle size decays like $\sigma_n=-1$19. The scaling profile $\sigma_n=-1$20 is positive, has an upper power-law bound near zero, and a super-exponential lower bound at large size. These results identify a dynamical-scaling regime in which the rescaled size distribution becomes stationary while mass is conserved (Jaiswal et al., 2024).

Finite-volume discretization makes these equations computationally usable. For the pure collisional breakage equation, a nonuniform-mesh finite-volume scheme with explicit Euler time stepping is shown to converge weakly under a CFL-type stability condition $\sigma_n=-1$21, with piecewise-constant reconstructions weakly compact in $\sigma_n=-1$22. On a uniform mesh, if $\sigma_n=-1$23, then

$\sigma_n=-1$24

yielding first-order convergence when $\sigma_n=-1$25; the reported numerical tests gave $\sigma_n=-1$26 (Bariwal et al., 2022). For the non-conservative approximation of the collisional-induced breakage equation, the finite-volume method likewise produces weak convergence under a feasible stability condition and an $\sigma_n=-1$27 error estimate

$\sigma_n=-1$28

again with first-order convergence on uniform meshes and numerical confirmation by $\sigma_n=-1$29 (Bariwal et al., 2024).

6. Regimes, applications, and recurring misconceptions

Several application domains show that breakage does not always manifest as a localized catastrophic crack. In depletion colloidal gels under constant shear stress, strand failures occur at a constant rate per particle per unit macroscopic strain, $\sigma_n=-1$30, independent of system size; the event-event radial distribution is flat, $\sigma_n=-1$31, and temporal correlations are negligible, so breakages are Poisson-like in accrued strain. At the macroscopic level, yielding is identified by a critical strain $\sigma_n=-1$32, the failure time obeys $\sigma_n=-1$33 with $\sigma_n=-1$34, and aging increases resistance through $\sigma_n=-1$35. The corresponding narrative is homogeneous, ductile failure rather than fracture localization or shear banding (Bhaumik et al., 2024).

Droplet impact provides a different inversion of the usual breakage logic: an interfacial instability can suppress fragmentation. For particle-coated droplets on superhydrophobic surfaces, fingering begins as early as $\sigma_n=-1$36, compared with $\sigma_n=-1$37 for clean droplets, and the relevant rim-Bond threshold is $\sigma_n=-1$38–$\sigma_n=-1$39. The instability increases viscous losses through finger formation, so the rebound stretching ratio can remain below the Rayleigh–Plateau threshold $\sigma_n=-1$40; moreover, finger pinch-off does not occur because $\sigma_n=-1$41 stays below the threshold $\sigma_n=-1$42. The result is an “island of stability” in the $\sigma_n=-1$43–ML map where droplet fragmentation is fully suppressed despite high impact energy (Lathia et al., 2021).

Turbulent agglomerates show yet another regime structure. In homogeneous isotropic turbulence laden with adhesive particles, the adhesion parameter

$\sigma_n=-1$44

controls sticking, rebound, and breakage. For doublets, the normalized breakage rate obeys

$\sigma_n=-1$45

while for agglomerates with primary-particle count $\sigma_n=-1$46, the normalized rate scales linearly as $\sigma_n=-1$47, with $\sigma_n=-1$48. Violent collisions and breakages are associated with particles ejected from strong vortices into straining sheets, emphasizing that flow topology can shape breakage statistics even when the criterion is formulated at collision (Chen et al., 2020).

Three recurring misconceptions are corrected by this body of work. First, breakage is not exclusively mechanical: quasi-one-dimensional nanostructures can rupture via thermal activation (Nisoli et al., 2013). Second, threshold crossing does not always imply immediate failure: ductile rupture requires subsequent energy accumulation, and local precursors can persist for many characteristic times before actual rupture (Marchioli et al., 2015, Bhaumik et al., 2024). Third, the first break is not determined solely by the largest local load: in the pull-or-jerk chain, changing the loading rate moves the breaking point from the pulled end toward the wall and back again because intrinsic and extrinsic timescales compete (Baek, 2018). Taken together, these results define breakage mechanics as a multiscale discipline in which rupture location, completion time, fragment geometry, and macroscopic response emerge from the interaction between constitutive thresholds, transport processes, and internal structural evolution.