PaperlayerView exact revision on arXiv

Source-linked AI summary

A mesh-free multiresolution deep energy method with phase-field modeling of brittle fracture

Han Zhang, Mehrisadat Makki Alamdari, Babak Shahbodagh, Mohammad Vahab, Cosmin Anitescu, Timon Rabczuk, Elena Atroshchenko

arXiv:2608.24126v1cs.LGmath.NA

TL;DR

Brittle-fracture phase-field simulations must resolve a localization band whose unknown path challenges mesh-based discretization. This paper uses a mesh-free neural representation with multiresolution C1 encoding and per-iteration resampling, matching finite-element responses and correctly classifying 90% of seeded cracks.

  • Problem

    Phase-field fracture simulations require resolving a narrow, unknown-path localization band, creating demanding discretization and integration requirements.

  • Method

    A single neural network minimizes incremental energy using multiresolution C1 B-spline features and integration points redrawn each iteration.

  • Results

    90% of seeded cracks were correctly classified as active or dormant, while single-edge-notched peak loads agreed with finite-element references within about 1%.

  • Takeaways & Limitations

    The encoding–resampling pairing supports crack propagation and lets second- and fourth-order fracture models share one discretization.

  • Takeaways & Limitations

    Boundary-condition corners can trigger artificial damage, so some runs require a grip-smoothing device whose post-failure side effects are documented.

Abstract

from arXiv · show

Phase-field modeling of brittle fracture removes the need to track cracks explicitly by recasting their evolution as the minimization of an energy functional. In return it requires a discretization dense enough to resolve a localization band whose width is set by a regularization length and whose path is not known in advance. We propose a mesh-free discretization in which a single neural network represents the displacement and phase fields and is trained by minimizing the incremental energy directly. The coordinates enter the network through a multiresolution feature encoding built from $C^1$ quadratic B-spline grids, so the finest scale the representation can express is set by choice rather than reached through slow training, and the energy is estimated by stratified Monte Carlo integration on points redrawn at every optimizer iteration. This pairing proves critical, since the crack fails to advance both when the integration points are held fixed and when the encoding is too coarse to represent the band, while each ingredient tolerates a wide range of settings once the other is in place. Because the representation is globally $C^1$, the second- and the fourth-order fracture energy densities run on the identical discretization. Across six problems, from single-edge-notched tension and shear to a thick-walled ring on a single spline patch, the computed load-displacement curves follow staggered finite element references at matched regularization length, with peak loads within about 1% on the single-edge-notched tests and within 8% where the crack pattern changes topology. On a public benchmark dataset of random multi-crack configurations the method classifies the active or dormant state of 90% of the seeded cracks in twenty zero-shot runs, where the deep Ritz baseline of the dataset authors fails.

1 Introduction

The work addresses phase-field fracture’s need to resolve an unknown crack band without explicit crack tracking by combining multiresolution B-spline neural representations with resampled energy integration. Its discretization makes the band representable from the first iteration and supports both second- and fourth-order fracture energies on the same globally C1 representation.

  • Motivation: Phase-field regularization removes explicit crack bookkeeping by recasting nucleation, propagation, kinking, branching, and merging as outcomes of energy minimization.The approach replaces the unknown crack set with a regularized phase field while retaining the energetic formulation.
  • Motivation: The discretization must resolve a crack band whose width is proportional to l, although the crack path is unknown in advance.This creates a quantitative resolution challenge for methods that would otherwise need global or adaptive refinement.
  • Motivation: A fixed collocation set caused phase-field cracks to advance faster than finite element references, with the error depending on activation function.Piecewise linear activation reproduced the reference, whereas smooth activations failed in the reported single-edge-notched tension and shear tests.
  • Contributions: One network uses multiresolution C1 quadratic B-spline grids whose finest level sets the smallest representable scale, making crack-band resolution a design choice rather than a training outcome.The method retains mesh-free automatic differentiation while encoding coordinates with structured multiresolution features.
  • Contributions: The globally C1 field representation lets second- and fourth-order fracture densities use the same discretization, with the higher-order form requiring one added differentiation pass.The fourth-order density therefore avoids the interface conditions otherwise introduced by the formulation.
  • Contributions: Redrawing integration points at every optimizer iteration and choosing sufficiently fine encoding levels are both required, because removing either device stops crack advance.The resampling prevents a fine representation from deteriorating between fixed integration points, while the encoding supplies the needed expressiveness.

2 Phase-field modeling of brittle fracture

The paper formulates brittle fracture as phase-field minimization of a regularized energy, with the length l controlling diffuse-crack width. It compares second- and fourth-order fracture densities, tensile-driven damage, and practical numerical treatments for higher regularity and corner singularities.

  • Regularized energy: The regularized fracture energy recovers Gc per unit crack length and Γ-converges to the Griffith functional as l approaches zero.At finite l, the model is interpreted as gradient damage with an internal length, and comparisons use the same l.
  • Second- and fourth-order models: The fourth-order density produces a C1, faster-decaying, more compact crack band than the second-order density at equal l, while both dissipate Gc per unit crack length.The fourth-order model also provides smoother nearby stress fields and more accurate surface-energy approximation, but requires H2 regularity.
  • Boundary conditions: Energy minimization avoids explicitly imposing the higher-order natural boundary conditions required by residual-based formulations for the fourth-order phase-field equation.The second-order condition is ∇ϕ · n = 0, whereas the fourth-order case also involves the normal derivative of ∆ϕ.
  • Damage driving force: The default hybrid formulation combines isotropic mechanical degradation with phase-field evolution driven only by spectral tensile energy, preventing damage in compressed material.Its single staggered energy formulation advances displacement and phase-field updates together, with gradients matching an explicitly staggered implementation to machine precision.
  • Corner treatment: Corner stress singularities can create unbounded fracture driving forces when l is resolved near changes from fixed to loaded boundaries.A smooth grip bump suppresses this numerical damage without changing the incremental problem structure, though its post-failure side effects are documented separately.

3 Mesh-free deep energy method

The mesh-free deep energy method represents displacement and phase fields with one neural network and minimizes incremental energy using Monte Carlo integration. C1 multiresolution features and fresh stratified sampling address phase-field localization and quadrature-induced failure modes.

  • Field representation: One shared neural network supplies raw displacement and phase outputs, which closed-form constructions convert into fields satisfying essential boundary conditions and carrying the pre-existing crack exactly.The incremental problems are solved by first-order optimization until the estimated energy plateaus.
  • Multiresolution encoding: C1 quadratic spline interpolation keeps strains and phase gradients continuous while making second derivatives piecewise constant, allowing second- and fourth-order fracture terms on one discretization.C0 linear interpolation would make the fourth-order term non-integrable at cell boundaries.
  • Multiresolution encoding: h = 0.26l sets the finest feature spacing, giving about four features across one profile decay length and roughly fifteen across its 4l half-width.The finest level is prescribed as a fraction of the regularization length so the localized damage band remains representable.
  • Monte Carlo integration: nr = 1 redraws stratified integration points every optimizer iteration within each load step, preserving the incremental minimization structure while preventing fixed-point zero-energy modes.The density is frozen within each load step, and stratification controls the variance introduced by fresh sampling.
  • Failure-mode controls: Neither resolution capping nor resampling suffices alone: the cap cannot prevent fixed-quadrature artifacts, while resampling cannot remove derivative oscillations from overly fine representations.The method therefore pairs a function-space resolution cap with iteration-wise point redrawing to address the two characteristic failure modes.

4 Numerical examples · 4.1 Single-edge-notched tension and shear

Across six examples, the method is quantitatively assessed against finite element and other reference solutions, with single-edge-notched tension and shear providing detailed benchmarks. In these tests, both model orders reproduce crack paths and key load–displacement features, while the fourth-order bands are more compact.

  • 4 Numerical examples: Six examples assess quantitative agreement, crack branching, three-crack coalescence, elastic response against a closed-form solution, nucleation, and random multi-crack configurations.The single-edge-notched tests provide quantitative comparisons for both second- and fourth-order models.
  • 4 Numerical examples: The finite element references use a 512 × 512 bilinear-element mesh, element size l/5, and 789 507 displacement-and-phase-field degrees of freedom.The references were generated with an in-house staggered solver distributed with the paper’s code.
  • 4.1 Single-edge-notched tension and shear: Single-edge-notched tension and shear use a pre-existing mid-height crack, displacement-driven loading, and identical discretizations for second- and fourth-order models.The models differ only in the density ψc,n entering Eq. (1), and each test compares both solutions with the finite element reference.
  • 4.1.1 Tension: Under tension, the crack crosses the ligament horizontally, and the second-order peak of 97.6 N reproduces the reference peak of 98.7 N within 1.1%.The corresponding peak displacements are δ = 4.90 × 10^-4 mm and 5.09 × 10^-4 mm, respectively.
  • 4.1.2 Shear: Under shear, the crack curves toward the bottom edge, and the second-order peak of 66.42 N agrees with the reference value of 66.35 N within 0.1%.The computed response also reproduces the reference’s softening dip and second rise, while retaining a single band without a spurious secondary branch.
  • 4.1.2 Shear: In shear, computed intermediate and final fields agree with finite element references, the full curved crack path is reproduced, and the fourth-order band is more compact.The two specimen parts separate across the fully developed band in agreement with the reference fields.
  • 4.1.2 Shear: Four random-seed shear runs complete the loading program and retain the response features, with a peak load of 66.3±0.3 N.They produce the same single curved crack band while varying initialization and sampling sequence.

4.2 Crack branching · 4.3 Coalescence of en-echelon cracks · 4.4 Plate with a circular hole

The topology tests show that the method reproduces branching and coalescence patterns without explicit crack tracking, while the masked circular-hole example nucleates cracks at the analytically predicted stress-concentration sites. Across these cases, computed responses and crack paths agree with corresponding finite element or closed-form references where available.

  • 4.2 Crack branching: The undecomposed-energy branching test produces a symmetric Y-shaped crack through incremental-minimization bifurcation rather than dynamic crack branching.The comparison is therefore made with a finite element solution using the same functional, not with experiment.
  • 4.2 Crack branching: 59.3 N is the branching peak load, about 8% below the finite element reference computed with the same undecomposed energy.The branches curve smoothly apart, arrest near the right edge, and remain symmetric within the band width.
  • 4.3 Coalescence of en-echelon cracks: The three 45° en-echelon cracks link through curving inner tips and form a single staircase crack under vertical tension.The multi-segment seeded profile is formed by the minimum distance to the three crack segments, with no other method ingredient changed.
  • 4.3 Coalescence of en-echelon cracks: 195.4 N at 6.4 × 10^-4 mm is the coalescence peak, 5.8% above the finite element reference of 184.8 N at 6.1 × 10^-4 mm.Both solution orders link the three cracks along the same staircase pattern.
  • 4.4 Plate with a circular hole: The circular hole is represented by an indicator mask in the estimator, without boundary fitting, for a unit-square plate with a 0.1 mm central hole.The specimen is pulled horizontally in uniaxial tension.
  • 4.4 Plate with a circular hole: The elastic-stage hoop stress matches the Kirsch solution, including a stress concentration factor of three at θ = ±90°.The comparison verifies the combined domain mask, sampling, quadrature, and automatically differentiated stresses before fracture is released.
  • 4.4 Plate with a circular hole: Damage nucleates at θ = ±90° and grows into two symmetric cracks toward the top and bottom edges.The nucleation site follows the elastic maximum-hoop-stress concentration, and the load–displacement curve shows a single sharp failure event.

4.5 Thick-walled ring · 4.6 Random multi-crack configurations

The method reproduces the prescribed crack path and force response in a thick-walled ring on one NURBS patch, then evaluates zero-shot behavior on random multi-crack benchmarks under tension and shear.

  • 4.5 Thick-walled ring: The ring uses a symmetric half-domain with two outer notches, vertical loading on the upper semicircle, and a fixed lower semicircle.The ring has inner radius 5 mm, outer radius 20 mm, and 3 mm notches; cracks grow horizontally toward the center.
  • 4.5 Thick-walled ring: A single quadratic NURBS patch represents the curved specimen, while fields and encoding live on the parametric square and sampling uses physical measure.The same patch supplies the geometry, control net, basis functions, pullback strains, and Jacobian-weighted integration.
  • 4.5 Thick-walled ring: 3.91 kN is the computed reaction-force maximum at δ = 2.67 × 10−2 mm before softening as the crack crosses the wall.The crack follows the reported horizontal path and severs the full ligament down to the inner bore.
  • 4.5 Thick-walled ring: The displacement fields satisfy the prescribed upper and lower boundary motions and jump across the severed ligament.The pulled upper semicircle moves vertically, the lower semicircle remains at rest, and both displacement components reflect the fracture separation.
  • 4.6 Random multi-crack configurations: The benchmark contains 1,000 random configurations with 10–20 interior cracks, evaluated under biaxial tension and shear using staggered finite elements on an 800 × 800 mesh.The method tests ten published configurations under both loadings, with identical solver settings, no training data, and no per-sample adjustment.
  • 4.6 Random multi-crack configurations: Per-crack active/dormant classification is the discriminating benchmark measure because reference solutions leave several seeded cracks dormant while neighboring cracks grow, link, and percolate.Dice similarity is also reported, but classification directly tests whether the seeded crack states match the reference.
  • 4.6 Random multi-crack configurations: 0.5% is the maximum force disagreement at the reference peak displacement under both loadings, although computed failure occurs later than the reference.For configuration 106244, the tension peak is 4418 N versus a reference peak of 4223 N, with the discrepancy attributed to different load increments and metastable-branch tracking.
  • 4.6 Random multi-crack configurations: 87 of 88 per-crack decisions remain identical across eight initialization runs, while peak loads vary by at most 2.0% in tension and 1.1% in shear.The deep Ritz baseline produced no repeated crack pattern across its initializations, contrasting with the stability attributed to resampled integration points.

5 Conclusions

The proposed mesh-free phase-field framework combines a single network, multiresolution C1 encoding, and per-iteration resampling to minimize incremental fracture energy without fixed meshes. It matches finite-element benchmarks and identifies most seeded cracks, while its conclusions remain limited to fixed regularization lengths and a non-variational hybrid formulation.

  • Contributions: The framework represents displacement and phase fields with one network, sets the finest resolvable scale through multiresolution C1 encoding, and redraws integration points every iteration.It treats second- and fourth-order fracture energy densities on the same discretization and imposes essential boundary conditions exactly through lifts.
  • Benchmark results: Within about 1% at the peaks, single-edge-notched load–displacement curves agree with staggered finite element references at the same regularization length.The curves also reproduce the softening valley and second local maximum of the shear response.
  • Benchmark results: 90% of seeded cracks are classified as active or dormant correctly in twenty zero-shot runs, while the dataset authors’ deep Ritz baseline fails on these configurations.The reported mean Dice scores are 0.739 in tension and 0.824 in shear, against 0.680 and 0.733 for the authors’ surrogates.
  • Ablation findings: Removing either per-iteration resampling or sufficiently fine encoding stops crack advancement, whereas the method tolerates a wide range of settings when the other ingredient is present.Coarsening the encoding beyond the band width or holding integration points fixed prevents the crack from advancing.
  • Limitations: The study reports results only at a fixed regularization length per example, makes no sharp-crack-limit claim, and uses a hybrid formulation that is not variationally consistent.The authors accept this inconsistency because spectral decomposition stalls the mode II shear test.
  • Extensions: The smooth representation and unfixed quadrature extend to coupled multi-field problems, allowing staggered subproblems to integrate on one freshly drawn point set without interpolation between discretizations.The estimator is indifferent to the number of fields and whether they use one network or several.

Appendix A Convention map and surface-energy check for the fourth-order model

Appendix A maps the phase-field and length conventions to Borden et al. and verifies that both second- and fourth-order fracture densities recover the prescribed surface energy. The implementation check returns 0.042467 per unit crack length for both models, with a −0.007% error.

  • Convention map: The convention map identifies Borden et al.’s phase field c and length ℓ0, with the present decay-length parameter related by l = 2ℓ0.The mapping follows by equating the decay lengths of the second-order optimal profiles.
  • Surface-energy check: Each optimal fracture profile integrates to exactly Gc per unit crack length; second-order bulk and gradient terms contribute equally, while fourth-order terms split 5 : 2 : 1.For the fourth-order model, the bulk, gradient, and Laplacian contributions sum to Gc.
  • Implementation check: 0.042467 per unit crack length is returned by both models for l = 0.01 and Gc = 0.04247, corresponding to an error of −0.007%.The check evaluates closed-form profiles through the solver’s energy routines rather than solving a boundary-value problem.
  • Implementation check: 0.026544, 0.010617, and 0.005306 are the implemented fourth-order term integrals, compared with 0.026544, 0.010618, and 0.005309 analytically.The numerical check uses automatic-differentiation derivatives and trapezoidal integration on a half-width 30l window with 4×10^5 points.

Appendix B Hyperparameters

Production runs used shared network and optimization settings, with resolution controlled mainly by multiresolution feature-grid parameters. Runs used a single V100 GPU, while fourth-order discretizations were unchanged from second-order counterparts and cost 1.4–1.6 times more.

  • Computational cost: 1.4 to 1.6 times as much compute was required by fourth-order runs as by second-order counterparts on identical discretizations.Every production run used a single NVIDIA V100 GPU with 12 CPU cores.
  • Parameterization: 380 803, 1 309 059 and 2 287 491 total parameters were used for unit-square specimens, multi-crack configurations and the ring, respectively.The perceptron contributed 51 075 parameters, while the feature grids determined most of the total count.
  • Fracture and scaling parameters: γir = 103, κ = 10−6 and τ = 0 were shared settings, with τ = 0.02 used only for multi-crack shear runs.Displacement scales varied by specimen, including Uref = 10−3 mm for unit-square specimens and 3.3 × 10−2 mm for the ring.
  • Reference-solver comparison: 500 to 800 load steps were used by finite element references versus 80 to 126 deep energy steps, making the finite element solver roughly an order of magnitude cheaper per increment.The references used four CPU cores, compared with the production runs’ single V100 GPU and 12 CPU cores.
  • Shared training settings: The production network used four hidden layers of width 128, learning rates 5 × 10−4 for the network and 2 × 10−3 for feature grids, and two feature channels per level.Fourth-order runs used the unchanged discretization of their second-order counterparts.

Appendix C Failure modes of alternative configurations

Appendix C evaluates plausible alternatives to the production configuration on the shear test, changing one ingredient at a time while holding the others fixed. The failures reproduce documented phase-field or deep-Ritz failure modes, while successful variations measure each ingredient’s tolerance.

  • Experimental setup: Alternative configurations were examined on the shear test with one ingredient changed at a time and all others held fixed.The comparison targets the proposed method’s production configuration from Section 3.
  • Failure modes: The alternatives that fail reproduce failure modes documented in the phase-field or deep-Ritz literature.These failures provide the appendix’s diagnosis of why plausible configurations are unsuitable.
  • Configuration tolerance: Variations around the production setting measure the range tolerated by each ingredient.The alternatives were used not only to identify failures but also to assess robustness around the selected configuration.

C.1 Reuse of the integration points

Reusing one integration-point set throughout a load step decouples the estimated-energy descent from the true energy, preventing fracture despite continued optimization. Per-iteration resampling restores the 68.0 N peak and complete failure, while reuse for 10 iterations remains unchanged because the representation cannot learn the point set.

  • C.1 Reuse of the integration points: 161.5 N at δ = 1.2 × 10−3 mm is reached without failure when integration points remain fixed through a load step.The optimizer still lowers the estimated energy, but the reaction force climbs monotonically to more than twice the production peak.
  • C.1 Reuse of the integration points: Per-iteration resampling restores the response to a 68.0 N peak and complete failure by preventing fine encoding levels from lowering energy between points.With fixed points, the estimator cannot see fields deteriorating between points, so estimated-energy descent decouples from true-energy descent.
  • C.1 Reuse of the integration points: 74 is the median within-step reduction factor for the frozen energy estimate, while redraws undo it by a median factor of 81 and as much as 218.The apparent descent reflects worsening discrepancy between the fields and one particular point set, not a corresponding reduction in energy.
  • C.1 Reuse of the integration points: nr = 10 iterations of point reuse leaves the response unchanged, showing that the estimator needs a point set no representation can learn rather than simply more points.The seeded band becomes visibly speckled as reuse increases, indicating a continuous progression from the frozen-point failure.

C.2 Resolution of the encoding

The encoding must resolve the fracture band for crack propagation, but performance remains stable across a broad range of finer and coarser resolutions. A level four times finer than the regularization length preserves the response, whereas an encoding unable to localize the band prevents crack advance.

  • C.2 Resolution of the encoding: 71.0 N peak load, identical crack path, and identical failure event persist when the finest encoding level increases from 384 to 1536.The 1536 cap is four times finer than the regularization length requires.
  • C.2 Resolution of the encoding: 65.9 N still reproduces the response with levels 32/128 and spacing h = 0.79l.This coarser encoding remains effective with resampling active and the point budget unchanged.
  • C.2 Resolution of the encoding: 82.0 N at the final increment occurs without failure when a single resolution level of 4 sets h at 50l and cannot localize the fracture band.The seeded crack remains stationary, with only faint diffuse damage beyond its tip.
  • C.2 Resolution of the encoding: Neither encoding resolution nor point resampling suffices alone: removing either ingredient stops crack propagation.The tested threshold is broad, including a cap four times finer, a cap three times coarser, and a redraw interval ten times longer than production.

C.3 Spectral decomposition under shear

Under spectral decomposition, the mode II shear test stalls because only tensile principal strains drive fracture, leaving the crack at its seeded tip and the reaction increasing without failure.

  • C.3 Spectral decomposition under shear: The crack never advances beyond the seeded tip, while the reaction rises monotonically to 104.8 N at δ = 1.2 × 10−3 mm without failure.This behavior occurs when spectral decomposition replaces the hybrid driving energy.
  • C.3 Spectral decomposition under shear: In pure shear, opposite principal-strain pairs mean only the tensile component drives fracture, giving the band about half the available driving force.The resulting deficit prevents crack advance at the loads reached in the test; the formulation, rather than the solver, causes the stall.

C.4 Corner damage under the volumetric–deviatoric decomposition

Without corner toughening, the volumetric–deviatoric decomposition produces spurious damage at the fixed-edge corner despite otherwise correct crack growth. The hybrid formulation avoids this corner damage on the identical loading program while producing a clean single fracture band.

  • Corner damage: Without corner toughening, the volumetric–deviatoric decomposition drives the fixed-edge corner damage near one, creating a spurious zone absent from the reference solution.The program used 36 increments of 5×10−5 mm.
  • Crack response: 64.4 N peak load and full ligament severance by δ = 1.8×10−3 mm show that the crack itself advances correctly despite the corner artifact.The corner damage is therefore a localized decomposition-related failure rather than a failure of crack propagation in the test.
  • Corner damage: All deviatoric deformation contributes to damage, so the shear-rich stress concentration at the fixed-edge corner supplies damage that a discretization resolving l everywhere cannot ignore.This mechanism explains the spurious corner zone under the volumetric–deviatoric decomposition.
  • Hybrid control: 68.5 N peak load and a clean single band with no corner damage are obtained by the hybrid formulation on the identical program without toughening.The hybrid control also records a reaction of 47.4 N at the reported displacement.

Appendix D Seed ensembles of the random multi-crack configurations

Appendix D examines four random-seed runs for configuration 106244 under biaxial tension and shear, comparing their responses with finite element references. The ensemble shows limited peak variation despite the baseline producing different crack patterns across initializations.

  • Seed ensembles: 2.0% and 1.1% are the peak spreads across four random seeds in tension and shear, respectively.The curves are overlaid on the dataset’s reference solutions.
  • Baseline comparison: The deep Ritz baseline produced a different crack pattern from every network initialization.This contrast motivates examining the proposed method’s seed ensembles.
  • Load–displacement response: Four random seeds of configuration 106244 are compared under biaxial tension and shear against the finite element reference.Figure 34 overlays the eight load–displacement curves with the reference solutions.
  • End-state fields: The end-state fields report displacement components u and v and the phase field for all four tension seeds on shared color scales.The corresponding shear ensemble is arranged analogously.
Loading 2608.24126v1…