Goal-oriented adaptive meshing for phase-field fracture in heterogeneous solids

Introduction

Variational fracture and phase-field fracture models have become standard tools for simulating complex crack initiation, propagation, coalescence, and branching without prescribing a crack path. The energetic formulation of brittle fracture introduced by Francfort and Marigo and regularised numerically by Bourdin, Francfort, and Marigo provides the conceptual foundation [1,2,3]. Later thermodynamically consistent finite-element formulations made the method accessible for multiphysics and engineering simulations [4]. Reviews of brittle and cohesive phase-field fracture now cover a broad range of formulations, degradation functions, split energies, and implementation choices [7,8].

The computational cost remains a central obstacle. The phase-field length scale must be resolved by several elements, and heterogeneous solids add sharp stiffness contrasts, interface toughness changes, and competing crack paths. Uniformly resolving the entire domain at process-zone scale is wasteful, but naive refinement near high damage can be too late: by the time damage localises, the preceding stress and energy fields may already have biased crack direction. Adaptive mesh refinement has therefore become an active area in phase-field fracture, with residual, damage-gradient, Nitsche-based, and problem-specific strategies all showing benefits [12,13,14,15,16,17].

Goal-oriented adaptivity is attractive because engineers often care about a quantity of interest rather than a global norm of displacement or damage error. A bridge notch-load capacity, a ceramic joint crack-arrival time, or a composite coupon energy release may be more important than uniform error reduction everywhere. Dual-weighted residual methods provide a principled way to weight local residuals by their influence on a chosen functional [9]. Wick demonstrated goal-functional evaluation for phase-field fracture using partition-of-unity dual-weighted residual mesh adaptivity, showing that adjoint information can guide refinement beyond the immediate crack band [10].

This article develops a practical goal-oriented adaptive meshing workflow for heterogeneous phase-field fracture. The contribution is algorithmic and diagnostic: we combine a dual problem for selected engineering goals, phase-field length-scale safeguards, interface-aware marking, and controlled coarsening. The method is evaluated against uniformly fine reference solutions and against simpler adaptive strategies across heterogeneous benchmarks.

Phase-field fracture model

The baseline model is a quasi-static small-strain phase-field fracture formulation with displacement u and scalar damage d. The stored elastic energy is degraded by a quadratic degradation function with residual stiffness, and fracture energy is regularised by a length scale ell. Unless otherwise stated, we use the AT2 crack-density functional with tension-compression energy split to reduce damage growth under compression, following the unilateral-contact motivation in regularised brittle fracture [5]. The irreversibility constraint is enforced by a history-field active-set procedure and checked against a primal-dual active-set solve on selected benchmarks.

Heterogeneity enters through spatially varying elastic moduli, fracture toughness, and, in interface benchmarks, cohesive-like reduced fracture resistance near material boundaries. This does not make the model interface-width insensitive by itself. Heterogeneous phase-field fracture is sensitive to the relation between length scale, inclusion size, and interface transition width, as shown in recent cohesive and stochastic heterogeneous formulations [18,19]. We therefore report all geometric dimensions in units of ell and avoid interpreting sub-ell path details as physical predictions.

The nonlinear problem is solved by a staggered alternate minimisation scheme with line search. Monolithic solves were tested for two benchmarks but did not change the reference responses enough to justify the additional implementation complexity in this study. Load stepping is adaptive: increments are reduced when active-set changes or energy residuals exceed prescribed tolerances. The displacement field uses quadratic triangles or tetrahedra, and the phase field uses linear elements on the same mesh. This mixed-order choice improves stress resolution without over-resolving damage variables far from cracks.

The crack path is extracted from the d = 0.85 contour for visual comparison and from the ridge of the phase-field gradient for quantitative Hausdorff-distance calculations. The ridge measure is less sensitive to contour threshold when damage saturates behind a crack. Energy dissipation is computed from the regularised fracture term, not from contour length alone.

Goal-oriented estimator

For a chosen quantity of interest J(u,d), the dual-weighted residual estimator approximates the error in J by weighting primal residuals with the solution of an adjoint problem. In a nonlinear phase-field fracture problem, this adjoint is path-dependent because the tangent operator depends on active damage constraints and the current crack state. We linearise about each converged load step and solve the adjoint on an enriched patch space. This follows the optimal-control view of a posteriori finite-element error estimation while keeping the implementation compatible with an existing fracture solver [9].

Four goal functionals are used. J1 is the reaction force at the loading boundary near peak load. J2 is crack-mouth opening displacement at a prescribed force level. J3 is dissipated fracture energy inside a target window around an inclusion or interface. J4 is an interface-arrival functional that integrates damage-gradient magnitude against a smooth indicator around a material boundary. J4 is not a probability by itself; it is a differentiable surrogate for whether the crack reaches a protected region.

The local error indicator combines three terms: a dual-weighted displacement residual, a dual-weighted phase-field residual, and a length-scale safeguard that enforces h <= ell/3 wherever d exceeds 0.15 or the history field has a strong gradient. The safeguard prevents the adjoint from under-refining the future process zone when the current goal functional is insensitive to a region that will soon localise. Interface-aware marking adds refinement along material boundaries where toughness or modulus jumps coincide with high adjoint weight.

Coarsening is allowed only when d > 0.98, the crack ridge has moved more than 4ell away, and the local contribution to all active goal indicators is below 1% of the marked-element median. This conservative rule prevents the solver from erasing residual compliance behind a crack too early. The mesh is regenerated by local bisection in two dimensions and newest-vertex subdivision in three dimensions, with interpolation of u, d, and history variables followed by a small energy-projection correction.

Adaptive algorithm

Each load step follows a solve-estimate-mark-refine loop. First, the nonlinear fracture problem is solved to the current load. Second, goal functionals and primal residuals are evaluated. Third, the adjoint problems for active goals are solved on the enriched patch space. Fourth, elements are marked by a combined DWR and length-scale criterion using a bulk parameter theta = 0.55. Fifth, the mesh is refined and, if allowed, coarsened. The load step is then repeated on the new mesh if the estimated goal error exceeds the tolerance; otherwise the algorithm advances.

Multiple goals are combined by normalising each indicator by the current magnitude or an engineering scale for that goal and then taking a weighted maximum. The default weights are equal, but the interface-arrival functional can be prioritised in protected-inclusion problems. This is important for heterogeneous solids: a mesh optimised for peak force may under-resolve a low-energy secondary crack that controls whether damage reaches an interface.

The reference comparison uses three other mesh strategies. Uniform refinement resolves the entire domain with h approx ell/4. Damage-gradient refinement marks elements where d, grad d, or the history field exceed thresholds. Residual-only refinement uses primal residuals without adjoint weights. All methods use the same nonlinear tolerances and load-step controller. This isolates mesh strategy from solver settings, a frequent source of misleading comparisons in adaptive fracture simulations [16].

The method is implemented in a finite-element code with PETSc linear solvers and MUMPS direct factorisation for two-dimensional benchmarks. The three-dimensional aggregate coupon uses algebraic multigrid preconditioning for displacement blocks and direct local solves for phase-field patches. Wall-clock timings are reported for a 32-core AMD EPYC node. Because implementation details matter, we report both degrees of freedom and wavefront-adjusted linear-solver time rather than element counts alone.

Benchmarks

The first benchmark is a notched three-point-bend beam containing 140 circular stiff inclusions with stiffness ratio 5:1 and fracture toughness ratio 1.8:1 relative to the matrix. Inclusion diameters range from 3ell to 12ell. The target functional is peak reaction force and dissipated energy in a 20ell window around the expected crack-inclusion interaction zone. The benchmark tests whether the mesh refines ahead of the crack near inclusions that may deflect the path.

The second benchmark is a fibre-reinforced plate with a central notch under tension. Fibres are represented as elongated high-stiffness, high-toughness regions embedded in a matrix. Two crack paths have similar energy: one cuts through the matrix between fibres, and one follows a weak interfacial band. The goal functional is crack arrival at a protected fibre cluster. This case is designed to expose the limitation of local refinement strategies when small stress differences decide path selection.

The third benchmark is a layered ceramic joint under mixed-mode loading with alternating layers of modulus and toughness. The model includes reduced interface fracture resistance over a transition width of 1.5ell. The goal is the mode-mix-sensitive energy dissipated in the joint layer and crack-mouth opening displacement. This benchmark tests interface-aware marking and the relation between phase-field length scale and material-layer thickness.

The fourth benchmark is a three-dimensional aggregate coupon with 312 randomly packed ellipsoidal inclusions. The mesh begins with 1.2 million degrees of freedom and reaches a maximum of 5.8 million in the adaptive run. A uniformly fine h = ell/4 mesh would exceed 22 million degrees of freedom and was used only on a cropped reference window. The goal is dissipated energy in a central fracture process volume and the average crack-plane normal.

Accuracy and cost

Across the first three two-dimensional benchmarks, the goal-oriented adaptive method predicted peak load within 1.9% of the uniformly fine reference solution, compared with 4.8% for damage-gradient refinement and 3.6% for residual-only refinement at similar element budgets. Dissipated fracture energy in target windows was within 2.7% for the goal-oriented method, 6.4% for damage-gradient refinement, and 5.1% for residual-only refinement. Crack-path Hausdorff distance was below 1.6ell for all goal-oriented runs after path bifurcation was excluded from the distance calculation.

The cost reduction was substantial. In the notched beam, the adaptive mesh used 71% fewer peak degrees of freedom than uniform refinement and reduced wall-clock time by 64%. In the fibre plate, degrees of freedom fell by 58% and wall-clock time by 49%. In the layered joint, degrees of freedom fell by 81% and wall-clock time by 76%. The three-dimensional coupon achieved a 74% reduction in degrees of freedom relative to the estimated uniformly fine model for the same target region, although solver time reduction was smaller because adaptive meshes had less favourable sparsity patterns.

Goal-oriented marking refined regions that damage-gradient marking missed. In the notched beam, the adjoint field highlighted inclusions ahead of the notch before damage localised there. In the layered joint, it refined a stiffness-contrast corner that controlled mode mix but had low damage at early load. In the fibre plate, the interface-arrival goal forced refinement along a weak fibre-matrix band that residual-only marking reached too late. These cases explain why goal error dropped faster than global energy-norm error.

The effectivity index, defined as estimated goal error divided by observed reference error, ranged from 0.72 to 1.38 for smooth response segments. Near snapback and crack branching, it widened to 0.41-2.7. This is not surprising. The adjoint linearisation assumes a locally smooth map from solution to goal, while crack branching can make the map nonsmooth. We therefore use the estimator to guide refinement, not as a certified bound on all fracture events.

Comparison with existing adaptive strategies

The results are consistent with recent adaptive phase-field fracture studies showing that local refinement can reduce cost without losing crack geometry [13,14,15,16,17]. The difference is that the present method asks why refinement is needed. Damage-gradient refinement is reliable once a crack exists, but it can under-resolve future bifurcation zones. Residual indicators are broader but do not distinguish residuals that matter for the chosen engineering output from residuals that do not. Goal-oriented indicators add this discrimination.

The method also inherits known complications. A posteriori estimators for phase-field fracture with irreversibility constraints must account for active sets and nonsmooth evolution [11]. We handle this by solving adjoints on the converged active set and by enforcing length-scale safeguards independent of adjoint weight. This pragmatic choice improves robustness but means the estimator is not a rigorous bound for all constrained nonlinear steps.

For heterogeneous materials, interface treatment is as important as crack-band treatment. Cohesive and interface-width-insensitive formulations show that length-scale choices can bias crack-interface competition [18]. Stochastic heterogeneous phase-field studies further show that small material variations can change crack path and variance [19]. Our adaptive method cannot remove those modelling sensitivities. It can, however, make them visible by refining the locations where the goal functional is sensitive to interface choices.

Limitations

The method is intended for quasi-static fracture in solids where a phase-field regularisation is already acceptable. Dynamic branching, fatigue crack growth, ductile plasticity, and hydraulic fracture introduce additional time scales and state variables. Adaptive algorithms exist for dynamic and fatigue phase-field settings, but goal-oriented adjoints must then account for time history and possibly discontinuous events [6,13].

The benchmarks use synthetic heterogeneous microstructures with known material maps. Real materials may have uncertain inclusions, imperfect interface toughness, residual stress, and evolving process-zone mechanisms. A mesh-adaptive numerical method cannot compensate for wrong material parameters. In practice, adaptive fracture simulations should be paired with calibration, imaging uncertainty, and mesh-convergence studies on the quantities of interest.

The adjoint problems add computational overhead. In the two-dimensional cases, estimator and adjoint solves added 11-18% to total runtime before mesh savings. In the three-dimensional coupon, they added 24%. The overhead is justified when the goal-oriented mesh substantially reduces primal solve cost or when the selected goal would otherwise be inaccurate. For exploratory simulations where only qualitative crack pattern is needed, simpler damage-gradient refinement may be sufficient.

Finally, crack-path ambiguity remains. When two paths differ by less than approximately 1% in total energy, small changes in mesh, load step, or material sampling can switch the selected path. The adaptive algorithm can reduce discretisation bias, but it cannot make a physically ambiguous path unique. We recommend reporting alternative-path indicators and path-sensitivity tests whenever adaptive phase-field fracture is used for design decisions.

Conclusion

We developed a goal-oriented adaptive finite-element strategy for phase-field fracture in heterogeneous solids. The method combines dual-weighted residual indicators, phase-field length-scale safeguards, interface-aware marking, and conservative coarsening behind saturated cracks. Across two-dimensional and three-dimensional benchmarks, it preserved peak load, target fracture energy, and crack-path metrics while reducing degrees of freedom by 58-81% relative to uniformly fine meshes.

The main benefit is not simply fewer elements. Goal-oriented adaptivity refines regions that influence the engineering output before visible damage arrives there, which is especially important near inclusions and weak interfaces. The approach is most reliable for smooth crack-growth segments and less reliable near branching or snapback, where the adjoint linearisation loses predictive power. Future work will extend the estimator to time-dependent fatigue, probabilistic material maps, and coupled thermo-mechanical fracture with certified bounds on selected goal functionals.

Data and code availability

All benchmark meshes, material maps, solver input files, adaptive refinement logs, reference-solution summaries, crack-path extraction scripts, and post-processing notebooks are included in the supplementary archive. The finite-element code was compiled with GCC 12.3, PETSc 3.20, MUMPS 5.6, and METIS 5.1. Analysis scripts use Python 3.11, NumPy 1.26, SciPy 1.11, meshio 5.3, and PyVista 0.43.

References

1. Francfort, G. A. & Marigo, J.-J. Revisiting brittle fracture as an energy minimization problem. J. Mech. Phys. Solids 46, 1319-1342 (1998).
2. Bourdin, B., Francfort, G. A. & Marigo, J.-J. Numerical experiments in revisited brittle fracture. J. Mech. Phys. Solids 48, 797-826 (2000).
3. Bourdin, B., Francfort, G. A. & Marigo, J.-J. The variational approach to fracture. J. Elast. 91, 5-148 (2008).
4. Miehe, C., Welschinger, F. & Hofacker, M. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations. Int. J. Numer. Methods Eng. 83, 1273-1311 (2010).
5. Amor, H., Marigo, J.-J. & Maurini, C. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J. Mech. Phys. Solids 57, 1209-1229 (2009).
6. Borden, M. J., Verhoosel, C. V., Scott, M. A., Hughes, T. J. R. & Landis, C. M. A phase-field description of dynamic brittle fracture. Comput. Methods Appl. Mech. Eng. 217-220, 77-95 (2012).
7. Ambati, M., Gerasimov, T. & De Lorenzis, L. A review on phase-field models of brittle fracture and a new fast hybrid formulation. Comput. Mech. 55, 383-405 (2015).
8. Vignollet, J., May, S., de Borst, R. & Verhoosel, C. V. Phase-field models for brittle and cohesive fracture. Meccanica 49, 2587-2601 (2014).
9. Becker, R. & Rannacher, R. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numer. 10, 1-102 (2001).
10. Wick, T. Goal functional evaluations for phase-field fracture using PU-based DWR mesh adaptivity. Comput. Mech. 57, 1017-1035 (2016).
11. Walloth, M. & Wollner, W. A posteriori estimator for the adaptive solution of a quasi-static fracture phase-field model with irreversibility constraints. SIAM J. Sci. Comput. 44, B479-B505 (2022).
12. Gupta, A., Krishnan, U. M., Mandal, T. K., Chowdhury, R. & Nguyen, V. P. An adaptive mesh refinement algorithm for phase-field fracture models: application to brittle, cohesive, and dynamic fracture. Comput. Methods Appl. Mech. Eng. 399, 115347 (2022).
13. Freddi, F. & Mingazzi, L. Mesh refinement procedures for the phase field approach to brittle fracture. Comput. Methods Appl. Mech. Eng. 388, 114214 (2022).
14. Muixi, A., Fernandez-Mendez, S. & Rodriguez-Ferran, A. Adaptive refinement for phase-field models of brittle fracture based on Nitsche's method. Comput. Mech. 66, 69-85 (2020).
15. Rohracker, M., Kumar, P. & Mergheim, J. A comparative assessment of different adaptive spatial refinement strategies in phase-field fracture models for brittle fracture. Forces Mech. 10, 100157 (2023).
16. Xu, W. et al. An adaptive mesh refinement strategy for 3D phase modeling of brittle fracture. Eng. Fract. Mech. 284, 109241 (2023).
17. Badnava, H., Msekh, M. A., Etemadi, E. & Rabczuk, T. An h-adaptive thermo-mechanical phase field model for fracture. Finite Elem. Anal. Des. 138, 31-47 (2018).
18. Zhou, Q. Q., Wei, Y. G., Zhou, Y. C. & Yang, L. An interface-width-insensitive cohesive phase-field model for fracture evolution in heterogeneous materials. Int. J. Solids Struct. 256, 111980 (2022).
19. Wu, J.-Y., Yao, J.-R. & Le, J.-L. Phase-field modeling of stochastic fracture in heterogeneous quasi-brittle solids. Comput. Methods Appl. Mech. Eng. 416, 116332 (2023).