Localizing Gradient Damage Method (LGDM)

• The paper introduces LGDM as a gradient-based damage formulation that mitigates mesh dependence by confining damage to a finite-width band.
• It couples local damage evolution with nonlocal variables such as micromorphic fields or equivalent strains to represent diffuse network scission and fracture localization.
• LGDM provides a thermodynamically consistent framework that predicts fracture toughness emergently, making it effective in both elastomer and metallic fracture models.

Searching arXiv for papers on the Localizing Gradient Damage Method and closely related gradient-damage formulations. Localizing Gradient Damage Method (LGDM) denotes, in the literature cited here, a class of gradient-regularized damage formulations designed to suppress pathological mesh dependence while permitting damage to localize into a finite-width band. Across these formulations, the central idea is the introduction of a nonlocal field—directly the damage variable, or an auxiliary micromorphic quantity such as nonlocal chain stretch, nonlocal equivalent strain, or regularized damage invariants—together with an internal length scale that controls the width of the damage process zone. A distinctive strand of recent work, particularly for elastomers, frames LGDM as a thermodynamically consistent alternative to phase-field fracture in which fracture localization and diffuse network scission are both represented and fracture toughness may emerge as an output rather than an input (Mousavi et al., 5 Feb 2025). Other works use the same label for phase-field-like gradient damage, modified non-local damage laws with fixed-width bands, mixed finite-element regularizations, and micromorphic tensorial damage extensions (Saji et al., 30 Jun 2025, Sarkar, 2023, Velden et al., 2023).

1. Terminology and conceptual scope

The nomenclature surrounding LGDM is not uniform. In the elastomer formulation of Mousavi et al., the method is a stretch-based gradient-enhanced damage model with a local damage variable $d(X)\in[0,1]$, a perturbed-Lagrangian pressure $p(X)$, and a nonlocal chain-stretch measure $\bar\lambda(X)$. Damage is not assigned its own evolution PDE; instead, it is obtained from a constitutive damage function driven by the historic nonlocal stretch, while $\bar\lambda$ satisfies a micromorphic balance (Mousavi et al., 5 Feb 2025). In the modified non-local damage model of Saji et al., LGDM is realized through a nonlocal equivalent strain $\bar\varepsilon$ governed by a Helmholtz-type equation, while the local damage $d$ is updated through a chosen damage law $d=d(\bar\varepsilon)$ (Saji et al., 30 Jun 2025).

In other papers summarized under the same label, LGDM is closer to gradient-damage or phase-field practice. Alkhoury et al. formulate a regularized total energy in terms of displacement $\mathbf u$ and scalar damage $d$, with degradation $g(d)=(1-d)^2$, fracture toughness $p(X)$0, and internal length $p(X)$1, leading to a canonical Helmholtz-type damage PDE (Alkhoury et al., 3 Jul 2025). Riesselmann and Balzani employ a mixed finite-element formulation for finite strains with displacement $p(X)$2, damage $p(X)$3, and a Lagrange multiplier $p(X)$4 enforcing irreversibility exactly, while Holthusen et al. regularize anisotropic damage by introducing nonlocal invariants $p(X)$5 associated with a second-order damage tensor $p(X)$6 (Riesselmann et al., 2022, Velden et al., 2023).

This variation in usage suggests that LGDM is best understood not as a single constitutive model, but as a methodological family centered on localization-capable gradient regularization.

2. Variational structure and governing fields

The formulations collected under LGDM are variationally structured, but they differ in which field is regularized. In the elastomer model of Mousavi et al., the reference configuration is $p(X)$7, with deformation map $p(X)$8, deformation gradient $p(X)$9, Jacobian $\bar\lambda(X)$0, and right Cauchy–Green tensor $\bar\lambda(X)$1. The local elastic energy is compressible neo-Hookean,

$\bar\lambda(X)$2

and near-incompressibility is enforced by the perturbed-Lagrangian term $\bar\lambda(X)$3 with $\bar\lambda(X)$4. The nonlocal contribution introduces a stiffness $\bar\lambda(X)$5, an internal length $\bar\lambda(X)$6, and a damage-weighted relaxation function $\bar\lambda(X)$7 (Mousavi et al., 5 Feb 2025).

In that framework, damage is introduced through a statistical-mechanics-inspired conjugate force $\bar\lambda(X)$8. A critical chain stretch $\bar\lambda(X)$9, steepness $\bar\lambda$0, maximum damage parameter $\bar\lambda$1, and a history operator $\bar\lambda$2 define a constitutive law such that no damage occurs for $\bar\lambda$3, whereas for $\bar\lambda$4, $\bar\lambda$5 increases with the maximum historic nonlocal stretch. Solving $\bar\lambda$6 yields the damage function $\bar\lambda$7, so the only genuine additional PDE is the one for $\bar\lambda$8 (Mousavi et al., 5 Feb 2025).

Saji et al. formulate the free-energy density as

$\bar\lambda$9

with

$\bar\varepsilon$0

Here $\bar\varepsilon$1 carries the regularization, the stress is degraded by $\bar\varepsilon$2, and the nonlocal field is driven by $\bar\varepsilon$3 rather than by $\bar\varepsilon$4 alone (Saji et al., 30 Jun 2025).

By contrast, the phase-field-like version summarized by Alkhoury et al. uses the regularized total energy

$\bar\varepsilon$5

which yields mechanical equilibrium together with the damage balance

$\bar\varepsilon$6

Irreversibility is imposed by replacing $\bar\varepsilon$7 with the history field $\bar\varepsilon$8 (Alkhoury et al., 3 Jul 2025).

Micromorphic anisotropic extensions retain the same general logic. Holthusen et al. introduce a symmetric damage tensor $\bar\varepsilon$9, a tuple of local invariants $d$0, and their nonlocal counterparts $d$1, with micromorphic energy

$d$2

The regularized variables may be the full six independent tensor components, three trace invariants, or a reduced volumetric-deviatoric pair (Velden et al., 2023).

3. Localization mechanism and suppression of damage-band broadening

The principal theoretical motivation for LGDM is the avoidance of spurious localization pathologies. Classical local damage models are mesh dependent after loss of ellipticity. Conventional gradient-damage formulations remove that mesh dependence by introducing an internal length, but several papers emphasize that they may still exhibit unrealistic widening of the damage band as damage progresses (Saji et al., 30 Jun 2025). Mousavi et al. make the same criticism of classical gradient-enhanced damage models for elastomers, noting that they often fail to capture the cascade from damage to fracture and instead produce damage-zone broadening (Mousavi et al., 5 Feb 2025).

In the stretch-based elastomer formulation, two specific modeling features are introduced to overcome this deficiency. First, the local chain-stretch driver $d$3 is bounded by replacing it with

$d$4

so that once $d$5, the local driver saturates. Second, nonlocal coupling is relaxed inside highly damaged regions through

$d$6

According to the authors, these two features permit a physically meaningful diffuse damage field together with simultaneous fracture localization (Mousavi et al., 5 Feb 2025).

The modified non-local damage model achieves a related objective through two different thermodynamic conditions. The auxiliary degradation term

$d$7

is added so that $d$8 as $d$9, and hence the thermodynamic damage driving force $d=d(\bar\varepsilon)$0 vanishes at full damage. In parallel, the forcing term in the nonlocal equation is attenuated with

$d=d(\bar\varepsilon)$1

so that $d=d(\bar\varepsilon)$2 as $d=d(\bar\varepsilon)$3. The stated consequence is a non-growing constant-width damage band (Saji et al., 30 Jun 2025).

Earlier analytical localization studies provide a useful backdrop. In a one-dimensional rate-independent regularized damage model, the internal length $d=d(\bar\varepsilon)$4 controls the width of the localized zone, and linearization around the onset of damage gives a critical wavelength

$d=d(\bar\varepsilon)$5

Modified regularization with a variable length $d=d(\bar\varepsilon)$6 yields a constant-zone-width response for $d=d(\bar\varepsilon)$7, whereas $d=d(\bar\varepsilon)$8 produces expansion and $d=d(\bar\varepsilon)$9 contraction (Jirásek et al., 2014). This suggests that a central issue in LGDM is not only regularization per se, but the precise asymptotic behavior of the nonlocal coupling as the damage variable approaches its terminal state.

4. Irreversibility, fracture measures, and the distinction from phase-field fracture

Irreversibility is enforced in several ways. In the stretch-based elastomer model, the constitutive structure is governed by $\mathbf u$0 together with the Clausius–Planck condition $\mathbf u$1; because $\mathbf u$2 appears only algebraically, irreversibility is handled pointwise through the history-dependent constitutive relation $\mathbf u$3 (Mousavi et al., 5 Feb 2025). In phase-field-like variants, a history field $\mathbf u$4 stores the maximum previous energy density and enters the damage PDE directly (Alkhoury et al., 3 Jul 2025). In the mixed finite-element formulation of Riesselmann and Balzani, irreversibility is written as a KKT system,

$\mathbf u$5

with $\mathbf u$6 and exact enforcement by a Lagrange multiplier (Riesselmann et al., 2022).

A major conceptual distinction concerns fracture toughness. In standard phase-field formulations, $\mathbf u$7 is a prescribed material input appearing explicitly in the crack-density part of the free energy (Alkhoury et al., 3 Jul 2025). Mousavi et al. emphasize that their LGDM does not prescribe an input $\mathbf u$8. Instead, one loads a cracked or notched specimen in displacement control, computes the domain $\mathbf u$9-integral, and identifies the plateau in

$d$0

as the predicted fracture toughness. In their $d$1 plane-strain panel with a pre-damage zone, the plateau gives $d$2 (Mousavi et al., 5 Feb 2025).

That predicted value was then used as the input toughness in a phase-field benchmark with the same geometry, $d$3, $d$4, and $d$5. The reported force-displacement curves, $d$6-integral versus displacement curves, and damage-field versus phase-field contours all matched to within a few percent, while LGDM retained two stated differences: $d$7 had a micromechanical interpretation as a chain-statistical nonlocal length, and the damage field $d$8 followed chain-stretch physics rather than serving as a crack interpolant (Mousavi et al., 5 Feb 2025).

The comparison is not uniformly favorable to all gradient-enhanced damage variants. In an earlier elastomer study, the measured energy release rate during crack propagation in a phase-field model with artificial viscosity did not comply with the imposed critical energy release rate and showed non-monotonic behavior. The stretch-based GED alternative again produced $d$9 as an output, but was reported to suffer from numerical issues to which strain-based GED formulations are susceptible, including damage-zone broadening and high sensitivity of $g(d)=(1-d)^2$0 to artificial viscosity $g(d)=(1-d)^2$1 (Mousavi et al., 2024). This is one of the clearest controversies in the area: some damage-based GED formulations were said to be identical to phase-field methods and to inherit their limitations, whereas other LGDM variants were explicitly designed to avoid that equivalence (Mousavi et al., 5 Feb 2025).

5. Discretization strategies and computational implementations

LGDM has been implemented in several computational settings, and the implementation strategy depends strongly on which field is regularized.

Implementation line	Core discretization	Notable claim
Stretch-based elastomer LGDM	Taylor–Hood type $g(d)=(1-d)^2$2 for $g(d)=(1-d)^2$3	Fracture toughness extracted from $g(d)=(1-d)^2$4-integral plateau
Abaqus UMAT + UMATHT	Damage PDE mapped to heat equation	Applicable in large deformation
Mixed FE with Lagrange multiplier	$g(d)=(1-d)^2$5-$g(d)=(1-d)^2$6, bubble-enriched $g(d)=(1-d)^2$7-$g(d)=(1-d)^2$8, $g(d)=(1-d)^2$9-$p(X)$00	No penalty or stabilization
FFT-based nonlocal ductile failure	FFT-Galerkin mechanics + Krylov Helmholtz solve	Millions of DOFs in 3D RVEs
JAX-FEM differentiable GED	HEX8 with automatic differentiation	Mesh-independent, data-driven extension
Vectorized MATLAB LGDM	Sparse vectorized assembly	Significant speedup on consumer hardware

In the elastomer formulation of Mousavi et al., the finite-element spaces are $p(X)$01, $p(X)$02, and $p(X)$03, discretized with vector-valued quadratics for $p(X)$04, linears for $p(X)$05, and linears for $p(X)$06. A staggered quasi-static algorithm is used: first solve the mechanics block for $p(X)$07 via Newton–SNES, then solve the nonlocal PDE for $p(X)$08 with bound constraints using vinewtonssls, then update $p(X)$09 algebraically until $p(X)$10. Near incompressibility is enforced through $p(X)$11 (Mousavi et al., 5 Feb 2025).

Alkhoury et al. translate the large-deformation damage PDE into the spatial frame to exploit Abaqus’s heat-equation solver. With the identification $p(X)$12, the transient heat equation is matched term by term via

$p(X)$13

UMATHT prescribes $p(X)$14 and $p(X)$15, while UMAT uses the RPL mechanism with

$p(X)$16

The implementation was validated against an independent UEL code, with coincidence to within numerical precision in the reported decoupled tests (Alkhoury et al., 3 Jul 2025).

Riesselmann and Balzani propose a mixed finite element for gradient damage in finite strains based on a Lagrange multiplier and a bubble-enriched damage space. The damage field uses $p(X)$17-Lagrange interpolation on vertices plus one quartic bubble at the element center, and the multiplier is element-wise constant. Because the multiplier and bubble degree of freedom are local to each element, both are eliminated by static condensation, leaving a global system with only the $p(X)$18-displacement and vertex-based $p(X)$19-damage unknowns. The formulation is reported to require no penalty or stabilization (Riesselmann et al., 2022).

For large heterogeneous RVEs, an FFT framework solves the fully coupled problem in staggered form: a conventional mechanical problem by FFT-Galerkin under mixed macroscopic loading control, and a generalized Helmholtz equation for each nonlocal variable by preconditioned conjugate gradient. In three-dimensional particle-reinforced composites, $p(X)$20 runs with roughly $p(X)$21 voxels completed in a few hours on a moderate multicore workstation, and localized-band thickness scaled linearly with the regularization length $p(X)$22 (Magri et al., 2021).

Two implementation studies highlight opposite ends of the software spectrum. A vectorized MATLAB code reported runtime reductions from $p(X)$23 s to $p(X)$24 s for a one-dimensional bar with $p(X)$25 elements, from $p(X)$26 s to $p(X)$27 s for a $p(X)$28 two-dimensional SEN problem, and from $p(X)$29 s to $p(X)$30 s for a $p(X)$31 three-dimensional SEN benchmark, relative to non-vectorized baselines on a Ryzen 5 5600H laptop (Sarkar, 2023). At the other end, a differentiable JAX-FEM framework combines gradient-enhanced damage with physics-augmented neural networks; the nonlocal field $p(X)$32 obeys a PDE of the form

$p(X)$33

and convergence studies showed virtually identical damage evolution across three meshes in uniform tension and notch-tip damage $p(X)$34 across meshes of $p(X)$35, $p(X)$36, and $p(X)$37 elements in a notched-plate simulation (Wilkinson et al., 3 Apr 2026).

6. Applications, limitations, and evolving research directions

Applications of LGDM span quasi-brittle elastomers, metallic fracture, heterogeneous ductile failure, and anisotropic damage at finite strains. In elastomers, the method is used to represent large deformation, near incompressibility, diffuse chain scission, and crack propagation. Mousavi et al. explicitly compare LGDM with phase-field fracture and argue that both can deliver nearly identical macroscopic toughness and load response, while LGDM additionally carries a built-in model of diffuse network scission (Mousavi et al., 5 Feb 2025). A later elastomer study extends this perspective by relating the gradient-enhanced damage zone to a fractocohesive length

$p(X)$38

defined as the width of the region where irreversible bond scission and energy dissipation occur; the reported implication is that the full damage-zone width is typically of order $p(X)$39 (Mousavi et al., 30 Aug 2025).

For metallic fracture, the gradient-damage framework is coupled to mechanism-based strain-gradient plasticity, introducing separate length scales for fracture regularization $p(X)$40 and plastic hardening $p(X)$41. The reported effect of $p(X)$42 is twofold: it elevates crack-tip stresses and facilitates brittle-style failure for sharp defects, but delays localization and increases peak load for blunt defects (Shishvan et al., 2021). In multiphase microstructures, phase-specific characteristic lengths $p(X)$43 allow nonlocal variables to diffuse only within their host phase, producing mesh-objective stress-strain responses where local Gurson or Lemaitre damage models diverge strongly with grid size (Magri et al., 2021). For anisotropic damage, a reduced two-degree-of-freedom volumetric-deviatoric regularization reproduces the full six-degree-of-freedom response in post-peak softening and damage pattern, whereas a three-invariant alternative yields slightly wider damage bands and higher dissipation (Velden et al., 2023).

Several limitations recur. One is the longstanding tendency of classical gradient-damage models toward damage-zone broadening, which motivates the localizing modifications of both Mousavi et al. and Saji et al. (Mousavi et al., 5 Feb 2025, Saji et al., 30 Jun 2025). Another is numerical sensitivity under unstable crack growth, especially when artificial viscosity is introduced; in the elastomer comparison of phase-field and GED, the reported effect on $p(X)$44 was particularly pronounced for GED (Mousavi et al., 2024). A further source of confusion is conceptual rather than numerical: some “damage-based GED” formulations are described as identical to phase-field fracture and therefore as inheriting phase-field limitations, while other authors reserve LGDM for formulations that regularize auxiliary micromorphic drivers and derive damage algebraically (Mousavi et al., 5 Feb 2025).

Taken together, these works position LGDM as a broad research program rather than a single closed formalism. Its unifying objective is stable, mesh-independent softening with a finite process-zone width; its differentiating questions concern which internal variable should be regularized, how localization should be arrested at full damage, whether $p(X)$45 is an input or an emergent output, and how closely the regularized damage field should be tied to micromechanical physics.