Discrete Configurational Forces

• Discrete configurational forces are energy-based descriptors defined on discrete lattices, meshes, and grids that quantify energy changes from variations in material configuration.
• They are applied in fields like fracture mechanics, topology optimization, and defect transport to predict crack propagation, adaptive mesh refinement, and material behavior.
• Discrete implementations—from finite-element nodal forces to force ensembles in granular materials—provide precise directional insights and energy thresholds for material evolution.

Searching arXiv for the cited papers to ground the article in current research. Discrete configurational forces are energy-based material-force constructs defined on discrete lattices, pixel grids, finite-element meshes, defect lines, or discretized configuration spaces, and they measure how a system’s total energy changes under changes of material configuration rather than ordinary spatial displacement. In recent work, the concept appears in several technically distinct but related forms: equivalent nodal forces assembled from the Eshelby stress in nonlinear fracture and topology optimization, J-type equivalent domain integrals evaluated from measured elastic fields near blocked slip bands, generalized inner-variation forces and stresses in Kohn–Sham density functional theory, Peierls–Nabarro-modulated driving forces for domain walls in multistable media, Noether-current force densities in mesoscopic Cosserat continua with distributed defects, and a discrete counting of admissible force configurations in jammed granular matter (Santarossa et al., 16 Jul 2025, Stankiewicz et al., 21 Jul 2025, Koko, 25 Mar 2026, Motamarri et al., 2017, Jin et al., 23 Dec 2025, Steinberg, 14 Apr 2026, Sartor et al., 2019).

1. Variational basis and material-force definitions

In configurational mechanics, the central object is an energy–momentum tensor of Eshelby type. For finite-strain hyperelasticity, the mixed-mode fracture study defines

$$\Sigma(F)=\Psi\,I-F^{T}\frac{\partial\Psi}{\partial F},$$

with stored energy density $\Psi(F)$ and first Piola–Kirchhoff stress $P=\partial\Psi/\partial F$. In the topology-optimization formulation, the analogous quantity is written

$$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$

and the local configurational traction is $T_0=\Sigma\cdot N$ on $\partial B_0$. In small-strain blocked-slip analysis, the corresponding material tensor is

$$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$

which yields configurational forces by contour or domain integration (Santarossa et al., 16 Jul 2025, Stankiewicz et al., 21 Jul 2025, Koko, 25 Mar 2026).

The defining principle is that configurational forces act in material space. In fracture mechanics, the energy release rate $G$ is obtained from the flux of $\Sigma$ through a surface around the crack tip, with $G=-\,d\Pi/da$. In blocked-slip problems, Rice’s $\Psi(F)$0-integral appears as the $\Psi(F)$1 component of the configurational force, with the convention $\Psi(F)$2. In Kohn–Sham density functional theory, the same variational idea is implemented through inner variations $\Psi(F)$3, so that the derivative of the constrained ground-state energy is written in terms of a generalized Eshelby tensor, nuclear self terms, and, in pseudopotential calculations, additional non-local contributions (Motamarri et al., 2017).

A recurring distinction in these formulations is between spatial equilibrium and configurational balance. The blocked-slip formulation states the usual equilibrium equations $\Psi(F)$4, while the topology-optimization formulation derives the conventional force balance $\Psi(F)$5. Configurational forces are additional energetic descriptors: they govern where cracks advance, where defects are driven, where a mesh should refine, or which extension directions are energetically favored (Koko, 25 Mar 2026, Stankiewicz et al., 21 Jul 2025).

2. Discrete realizations

Recent work uses several discrete realizations of configurational-force concepts.

Setting	Discrete quantity	Construction
Nonlinear FE fracture	$\Psi(F)$6	$\Psi(F)$7
SIMP topology optimization	$\Psi(F)$8	$\Psi(F)$9
HR-EBSD blocked slip	$P=\partial\Psi/\partial F$0, $P=\partial\Psi/\partial F$1	J-type equivalent domain integral with $P=\partial\Psi/\partial F$2
Axisymmetric domain wall	$P=\partial\Psi/\partial F$3	$P=\partial\Psi/\partial F$4 with periodic modulation
Cosserat line defect	$P=\partial\Psi/\partial F$5	$P=\partial\Psi/\partial F$6
Granular force ensemble	$P=\partial\Psi/\partial F$7	discretized force-space multiplicity

In finite-element fracture, the discrete configurational-force vector at node $P=\partial\Psi/\partial F$8 is assembled elementwise as

$P=\partial\Psi/\partial F$9

The algorithm is explicitly post-processing based: after the nonlinear FE solve, one computes $$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$0, $$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$1, $$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$2, and $$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$3 at Gauss points; integrates element contributions $$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$4; assembles the global nodal vector; identifies nodes on the crack-surface edge plus a narrow “spurious” neighbor band; transforms forces to local cylindrical or Frenet frames; and sums the forces over crack-front segments per unit arc length $$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$5 (Santarossa et al., 16 Jul 2025).

In topology optimization, the discrete force is modified by the design field. Each element carries a pseudo-density $$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$6, from which one defines

$$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$7

with

$$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$8

The nodal configurational force then becomes

$$\Sigma=U_0\,I-F^T P, \qquad U_0=W_0(F;X)+V_0(\varphi;X),$$9

The scalar magnitude $T_0=\Sigma\cdot N$0 is used directly as an adaptive-mesh indicator (Stankiewicz et al., 21 Jul 2025).

In HR-EBSD slip-band analysis, a contour integral is converted to the J-type equivalent domain integral

$T_0=\Sigma\cdot N$1

where $T_0=\Sigma\cdot N$2 is a plateau function equal to $T_0=\Sigma\cdot N$3 inside an inner domain, $T_0=\Sigma\cdot N$4 outside an outer domain, and varying smoothly between them. On the pixel grid, the integral is evaluated by a Riemann sum using HR-EBSD displacement-gradient data, anisotropic elasticity for $T_0=\Sigma\cdot N$5-Ti, and finite differences of $T_0=\Sigma\cdot N$6 (Koko, 25 Mar 2026).

These constructions share a discrete evaluation strategy, but the source of discreteness differs by field. In FE applications, discreteness enters through the mesh and shape functions. In blocked-slip analysis, it enters through the measured pixel grid. In multistable lattices, it enters through the lattice spacing $T_0=\Sigma\cdot N$7 and the associated Peierls–Nabarro modulation. In granular materials, it enters through the subdivision of continuous force space into $T_0=\Sigma\cdot N$8-cells (Jin et al., 23 Dec 2025, Sartor et al., 2019).

3. Crack-front driving forces in soft fracture

For mixed-mode I + III fracture in soft, highly deformable solids, the discrete configurational-force method is implemented as a post-processing algorithm on finite-strain, neo-Hookean FE simulations with no body forces. The material model uses

$T_0=\Sigma\cdot N$9

with $\partial B_0$0 and $\partial B_0$1, together with $\partial B_0$2, $\partial B_0$3, and $\partial B_0$4. Loading is imposed by Dirichlet conditions: first an increasing vertical separation $\partial B_0$5 for Mode I, then an increasing bottom-to-top twist angle $\partial B_0$6 for Mode I + III (Santarossa et al., 16 Jul 2025).

The front is partitioned into planar-segment nodes and facet-edge nodes. For planar segments, the force is expressed in a cylindrical basis $\partial B_0$7; for facets, in the local Frenet triad $\partial B_0$8. The segmentwise measure is

$\partial B_0$9

and its magnitude and orientation are then analyzed (Santarossa et al., 16 Jul 2025).

The resulting interpretation is explicitly predictive. The scalar magnitude $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$0 is proportional to the local energy-release rate $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$1, and crack advance is expected when $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$2 exceeds a material-specific threshold $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$3, equivalently $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$4. The orientation of $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$5 indicates the preferred growth direction. On tilted facets, the dominant component lies along the binormal $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$6, showing that cracks advance while maintaining their tilt angle $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$7. A nonzero component along the facet normal $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$8 indicates a slight tendency for facet rotation, and under mixed-mode I + III the angle $$P_{kj}=W\,\delta_{kj}-\sigma_{ij}\,u_{i,k}, \qquad W=\tfrac12\,\sigma_{mn}\,\varepsilon_{mn},$$9 between $G$0 and the facet-plane normal can change sign for small $G$1, revealing shear-induced reorientation away from the original tilt (Santarossa et al., 16 Jul 2025).

The same discrete force field clarifies elastic interactions along complex crack fronts. On planar segments, $G$2 often exceeds the values on facets, especially under mixed-mode loading, which predicts the later growth of “type B” bridging cracks. Varying the facet spacing $G$3 shows that close spacing amplifies $G$4 on planar regions through constructive elastic interaction, whereas very small $G$5 leads to shielding. In facet-coalescence studies, mesh relaxation by DynaMesh-R is used to join facets smoothly, after which configurational-force redistribution shows a drop near the new junction and amplification farther away, producing an asymmetric driving-force distribution that controls non-uniform growth and the final echelon morphology (Santarossa et al., 16 Jul 2025).

4. Interface motion, pinning, and blocked-slip energetics

In elastically coupled multistable metamaterials, an axisymmetric domain wall of normalized radius $G$6 is described by a reduced-order model in which the total potential energy is the sum of plate bending and foundation contributions. In normalized form,

$G$7

$G$8

and

$G$9

Here $\Sigma$0 is the asymmetry parameter and $\Sigma$1 is the discreteness parameter, with $\Sigma$2 and $\Sigma$3 favoring expansion (Jin et al., 23 Dec 2025).

The configurational driving force is

$\Sigma$4

In the continuum limit $\Sigma$5, the force follows from differentiating $\Sigma$6. For finite discreteness, the model superposes a periodic Peierls–Nabarro correction of period $\Sigma$7 and amplitude $\Sigma$8,

$\Sigma$9

which yields a modulated configurational force with local minima whenever $G=-\,d\Pi/da$0 and $G=-\,d\Pi/da$1. These minima create pinning wells. The paper identifies three regimes—expansion, shrinking, and metastable pinning—and states that local minima exist only for $G=-\,d\Pi/da$2, within a finite band $G=-\,d\Pi/da$3 (Jin et al., 23 Dec 2025).

The continuum nucleation radius $G=-\,d\Pi/da$4 is determined by the global maximum of $G=-\,d\Pi/da$5. If the initial wall radius satisfies $G=-\,d\Pi/da$6, then $G=-\,d\Pi/da$7 and the wall shrinks; if $G=-\,d\Pi/da$8, then $G=-\,d\Pi/da$9 and the wall expands. The same axisymmetric reduced-order model extends to non-axisymmetric polygonal walls. For convex $\Psi(F)$00 polygons, a regular hexagon of side length $\Psi(F)$01 is metastable whenever $\Psi(F)$02 lies in the discrete stability set $\Psi(F)$03, and an irregular convex polygon is fully pinned if every edge length $\Psi(F)$04 satisfies $\Psi(F)$05. For concave $\Psi(F)$06 corners, pinning requires that the complementary regular hexagons be stable under the complementary asymmetry $\Psi(F)$07 (Jin et al., 23 Dec 2025).

A different interface-localization problem appears in blocked slip bands at grain boundaries in $\Psi(F)$08-Ti. There, the configurational force vector

$\Psi(F)$09

is evaluated from HR-EBSD measurements. Directionality is resolved by defining a Virtual Extension Direction and rotating the displacement-gradient field into the frame of each candidate slip-system trace in the neighboring grain. The polar function $\Psi(F)$10 exhibits peaks at the most energetically favorable extension directions. In the reported example, the maximum $\Psi(F)$11 occurs on a pyramidal $\Psi(F)$12 trace that is geometrically admissible but not the one with the highest Schmid factor; the highest-Schmid prismatic variant has only moderate $\Psi(F)$13, and a basal variant with $\Psi(F)$14 still produces non-zero $\Psi(F)$15. The paper therefore states that Schmid factor, $\Psi(F)$16, and residual Burgers vector do not uniquely predict the local energetic driving force $\Psi(F)$17 (Koko, 25 Mar 2026).

5. Adaptive computation and electronic-structure forces

In SIMP-based topology optimization, discrete configurational forces are used as a mesh-adaptivity criterion because, with the relaxed Eshelby-stress interpolation, they localize both in grey transition regions and in highly stressed regions. The localization mechanism is explicit in the paper: fully void elements have $\Psi(F)$18 and therefore $\Psi(F)$19; fully solid elements have $\Psi(F)$20, but $\Psi(F)$21 is often only moderate if energy density and stress vary smoothly; by contrast, grey transition elements generate large values because the density gradients align with $\Psi(F)$22, while high-stress solid regions also generate large values through the local Eshelby stress (Stankiewicz et al., 21 Jul 2025).

The marking rule uses the normalized maximum

$\Psi(F)$23

and thresholds

$\Psi(F)$24

with the common choice

$\Psi(F)$25

Refinement is aggressive, whereas coarsening is conservative: an element is coarsened only if all its vertices satisfy the coarsening threshold. The algorithm stores the refinement tree, allowing multilevel coarsening up to a user-defined maximum refinement level $\Psi(F)$26 (Stankiewicz et al., 21 Jul 2025).

The numerical results reported in the paper are problem-specific. For a cantilever compliance-minimization example with a $\Psi(F)$27 base mesh and max-level $\Psi(F)$28, configurational-force-based adaptivity achieves approximately $\Psi(F)$29 of the total CPU time of density-based adaptivity with identical compliance. For a U-beam example with compliance, volume, and $\Psi(F)$30-norm stress constraint, configurational-force-based adaptivity requires approximately $\Psi(F)$31 of the CPU time of density-based adaptivity for $\Psi(F)$32, while yielding equal or better enforcement of stress constraints. The paper contrasts this with density-based refinement, which refines everywhere grey, and von-Mises-based refinement, which refines stress-critical corners but can miss geometric boundaries (Stankiewicz et al., 21 Jul 2025).

In Kohn–Sham density functional theory, configurational forces are derived as generalized variational forces obtained from inner variations of the Kohn–Sham energy functional with respect to the position of a material point $\Psi(F)$33. The derivative of the ground-state energy has the form

$\Psi(F)$34

where $\Psi(F)$35 is the generalized Eshelby tensor, $\Psi(F)$36 the nuclear-self contribution, and $\Psi(F)$37 the non-local pseudopotential term. By choosing $\Psi(F)$38 as an atom-centered compact generator $\Psi(F)$39, one obtains atomic forces; by choosing $\Psi(F)$40, one obtains the stress tensor in a periodic cell (Motamarri et al., 2017).

The method is variationally consistent and the paper states that Pulay corrections are inherently included, so no separate Pulay-correction step is needed. It also treats pseudopotential and all-electron calculations in a single framework. In the reported higher-order FE benchmarks, forces and stresses show $\Psi(F)$41 convergence for CO, CH$\Psi(F)$42, N$\Psi(F)$43, SiF$\Psi(F)$44, Al fcc, and Li bcc, and the paper reports errors typically below $\Psi(F)$45 in forces and below $\Psi(F)$46 in stresses when compared with reference calculations (Motamarri et al., 2017).

6. Defects, Noether currents, and discrete line-force limits

A geometrically richer generalization appears in the mesoscopic Cosserat theory with distributed defects. The basic fields are the coframe $\Psi(F)$47 and an independent $\Psi(F)$48-connection $\Psi(F)$49, with torsion and curvature

$\Psi(F)$50

The Palatini-type action leads to Euler–Lagrange equations

$\Psi(F)$51

which combine the standard force and couple balances with defect-excitation terms $\Psi(F)$52 and $\Psi(F)$53 (Steinberg, 14 Apr 2026).

Material invariance then yields configurational currents as Noether quantities. Under infinitesimal material translations, the theory gives

$\Psi(F)$54

where $\Psi(F)$55 is the configurational stress current, $\Psi(F)$56 the configurational momentum, and $\Psi(F)$57 the configurational force density. Under material rotations, an analogous identity yields configurational moment density. The same framework also provides dynamic Bianchi transport laws,

$\Psi(F)$58

which connect defect transport directly to configurational forces and moments (Steinberg, 14 Apr 2026).

The discrete line-defect limit is especially explicit. If a defect line $\Psi(F)$59 carries

$\Psi(F)$60

then the net configurational force per unit length along $\Psi(F)$61 is

$\Psi(F)$62

The first term reproduces the classical Peach–Koehler force, and the second is a couple-stress correction. The worked example in the paper shows a time-dependent micro-rotation field generating torsion and curvature, which then produce defect excitations, force stress, and finally a localized configurational force that oscillates in space and decays as $\Psi(F)$63 (Steinberg, 14 Apr 2026).

This formulation explicitly separates defect densities from classical compatibility. Torsion and curvature are treated as independent primary fields rather than constrained to vanish. The theory therefore places discrete configurational forces within a broader geometric setting in which they are not only post-processed outputs but also Noether currents linked to defect transport (Steinberg, 14 Apr 2026).

7. Statistical-mechanical interpretation and common distinctions

In jammed granular materials, “discrete configurational forces” take a different form. The Force Network Ensemble fixes particle positions and treats the contact forces $\Psi(F)$64 as the degrees of freedom, subject to force balance on every grain and fixed global stress trace. The allowed set of force configurations is a convex polytope for frictionless packings, and the phase-space volume is

$\Psi(F)$65

The configurational entropy is then

$\Psi(F)$66

where $\Psi(F)$67 is the elementary force-cell volume, directly analogous to Planck’s constant in ordinary phase space (Sartor et al., 2019).

The direct-measurement protocol proceeds by generating jammed packings, constructing the rigidity matrix, computing the left nullspace to obtain the $\Psi(F)$68 states of self stress, intersecting the positivity cone with the unit sphere to obtain a convex spherical polytope of force solutions, and measuring its $\Psi(F)$69-dimensional volume $\Psi(F)$70. The paper reports

$\Psi(F)$71

and matches the resulting microscopic entropy with a macroscopic entropy obtained from angoricity via overlapping histograms (Sartor et al., 2019).

Taken together, these results suggest that discrete configurational forces are not a single formula but a family of energy-based descriptors tied to changes in material configuration. One common misconception is to identify them with ordinary force balance; the cited works instead distinguish spatial equilibrium from material driving forces. Another misconception is that geometric admissibility alone determines evolution: the blocked-slip study shows a marked decoupling between geometric metrics and the configurational-force response, and the domain-wall study shows that lattice discreteness can create metastable pinning even when the smooth continuum force would predict monotonic expansion or shrinking (Koko, 25 Mar 2026, Jin et al., 23 Dec 2025).

A further common feature is that thresholds and localization bands are problem dependent. In fracture, $\Psi(F)$72 must exceed a material-specific threshold $\Psi(F)$73. In multistable metamaterials, pinning requires $\Psi(F)$74 and occurs only over a finite interval of wall sizes. In blocked slip, the paper proposes a future critical threshold $\Psi(F)$75 for slip transfer or crack nucleation. This suggests a general pattern: discrete configurational forces provide the energetic ranking and directional information, while the onset of actual evolution requires constitutive or material-specific criteria supplied by the particular application (Santarossa et al., 16 Jul 2025, Jin et al., 23 Dec 2025, Koko, 25 Mar 2026).