Cohesive Phase-Field Model Advances

• Cohesive phase-field models are variational formulations that regularize cracks via a scalar damage field to recover finite strength and traction–separation behavior.
• They employ tailored degradation functions and eigenstrain-based techniques to capture elastic, pre-fractured, and fully fractured states under static and dynamic loading.
• Recent developments include inverse mapping of cohesive laws, anisotropic and multi-cohesive frameworks, and advanced finite-element implementations for multiphase materials.

Searching arXiv for recent and foundational papers on cohesive phase-field fracture and related models. arXivSearch: query="cohesive phase-field fracture", max_results=10 Cohesive phase-field model denotes a class of variational fracture formulations in which a crack is regularized by a scalar damage or phase-field variable, while the limiting fracture energy depends on the displacement jump through a cohesive law rather than through a purely Griffith-type surface term. In these models, intact material is typically represented by $v\approx 1$ or $d\approx 0$, fully broken material by $v\approx 0$ or $d\approx 1$, and the fracture process zone is replaced by a diffused band of width controlled by an internal length. The central objective is to recover finite strength, gradual energy release, and traction–separation behavior within a phase-field or gradient-damage framework, while retaining a variational structure amenable to $\Gamma$-convergence, equilibrium analysis, and finite-element implementation (Conti et al., 2014, Alessi et al., 16 Jul 2025).

1. Variational architecture

A cohesive phase-field model is usually formulated in terms of a displacement field $u$ and a scalar phase-field $v$ or $d$. In one broad family of models, the regularized energy has the form

$$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$

for $(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$, where $d\approx 0$0 governs stiffness degradation, $d\approx 0$1 encodes the cohesive scaling, and the terms $d\approx 0$2 and $d\approx 0$3 regularize damage localization (Alessi et al., 16 Jul 2025). A more specialized Ambrosio–Tortorelli-type version used for one-dimensional cohesive fracture is

$d\approx 0$4

with $d\approx 0$5 and $d\approx 0$6 a.e. (Bonacini et al., 2023).

The degradation mechanism is the constitutive core of the model. In variational cohesive formulations, $d\approx 0$7 is chosen so that damaged regions carry reduced stiffness but still induce a nontrivial limiting traction–separation response. In the foundational phase-field approximation of cohesive fracture, the elastic coefficient is computed from the damage variable $d\approx 0$8 through a function $d\approx 0$9 of the form $v\approx 0$0, with $v\approx 0$1 diverging near the undamaged state (Conti et al., 2014). This scaling differs qualitatively from brittle AT1/AT2 approximations because it is designed to generate an opening-dependent surface energy in the sharp limit.

Several later formulations preserve the same variational logic but modify the constitutive ingredients. Cohesive gradient-damage models for dynamics use a crack density functional linear in the damage field together with non-polynomial degradation functions, in order to obtain a linear elastic regime prior to damage onset and controlled strain-softening thereafter (Geelen et al., 2018). In eigenstrain-based cohesive phase-field models, cohesive behavior is not introduced by modifying the elastic degradation $v\approx 0$2; instead, fracture eigenstrains $v\approx 0$3 represent crack opening and sliding, and a strength potential $v\approx 0$4 governs crack nucleation, while the phase-field still regularizes propagation (Hageman, 23 Mar 2026). A further variational branch introduces an eigenstrain-like internal variable $v\approx 0$5 directly into the free energy and constructs a convex strength surface $v\approx 0$6 through its support function $v\approx 0$7, with degradation $v\approx 0$8 acting on the strength domain rather than on a conventional stiffness split (Vicentini et al., 13 Jun 2025).

These formulations are therefore best regarded as a family rather than a single equation. What unifies them is the replacement of a sharp discontinuity by a regularized field and the emergence, either in the $v\approx 0$9-limit or in the local constitutive response, of a cohesive law that depends on crack opening.

2. Cohesive law and sharp-crack limit

The defining property of a cohesive phase-field model is the emergence of a limiting free-discontinuity functional with jump energy $d\approx 1$0. In the foundational scalar model, the $d\approx 1$1-limit of the phase-field energies is

$d\approx 1$2

with $d\approx 1$3 determined by a one-dimensional optimal profile problem in coupled displacement–damage variables $d\approx 1$4 (Conti et al., 2014). In the generalized one-dimensional framework, the limit is

$d\approx 1$5

for $d\approx 1$6, with $d\approx 1$7 in the limit (Alessi et al., 16 Jul 2025).

The surface density $d\approx 1$8 is defined by a cell problem. In the generalized framework,

$d\approx 1$9

where $\Gamma$0 is a microscopic displacement profile with total jump $\Gamma$1, and $\Gamma$2 is a rescaled damage profile satisfying $\Gamma$3 (Alessi et al., 16 Jul 2025). In the earlier Ambrosio–Tortorelli-type construction, the analogous optimal profile problem yields a cohesive surface energy that is linear in the opening $\Gamma$4 at small values of $\Gamma$5 and has a finite limit as $\Gamma$6 (Conti et al., 2014).

This sharp-limit structure distinguishes cohesive from Griffith-type phase-field fracture. In brittle scaling, the same general family may instead converge to

$\Gamma$7

which is Griffith-type because the surface term is constant per jump and independent of opening (Alessi et al., 16 Jul 2025). By contrast, in cohesive scaling the jump cost is opening-dependent, typically nondecreasing, subadditive, and saturating.

In the one-dimensional cohesive approximation analyzed in detail in (Bonacini et al., 2023), the limiting sharp cohesive functional is

$\Gamma$8

with

$\Gamma$9

and $u$0, $u$1, $u$2 nondecreasing, Lipschitz with Lipschitz constant $u$3, and $u$4 as $u$5. In that notation, the threshold opening for complete decohesion is

$u$6

with $u$7 if $u$8 (Bonacini et al., 2023).

A common simplification is that phase-field fracture is intrinsically Griffith-like. The variational results above show that this is false in general: by choosing the degradation and damage terms appropriately, one obtains a cohesive free-discontinuity limit with a prescribed opening-dependent surface density.

3. Equilibria, critical points, and irreversible evolution

The sharp cohesive functional admits nontrivial critical configurations. In the one-dimensional setting with fixed elongation $u$9, critical points of the relaxed cohesive energy are completely classified by a constant global stress $v$0. They consist of elastic states, pre-fractured cohesive states, and fractured states (Bonacini et al., 2023). Elastic states have

$v$1

Pre-fractured cohesive states with $v$2 are piecewise affine with constant bulk strain $v$3 and finitely many jumps of equal opening $v$4 such that $v$5. Fractured states with $v$6 are piecewise constant with fully opened cracks satisfying $v$7 and $v$8 (Bonacini et al., 2023).

The phase-field approximations do not converge to arbitrary critical points of the limit. For critical points $v$9 of the phase-field energy with uniformly bounded energy and Dirichlet data $d$0, $d$1, $d$2, the asymptotic limit is selected by the minimum $d$3. Three regimes occur: $d$4 gives an elastic limit; $d$5 gives a single cohesive jump at $d$6; $d$7 gives a fully fractured limit with $d$8 (Bonacini et al., 2023). In particular, the regularized one-dimensional theory selects $d$9 critical points with at most one crack at the midpoint, excluding overstressed elastic states, diffuse singular strain, and multiple cracks.

Irreversible cohesive evolution requires additional internal variables. In the one-dimensional quasi-static theory of cohesive fracture derived from static phase-field models, the fracture state is described by a crack set $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$0 and a maximal opening history $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$1, and the energy uses a two-variable density $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$2, where $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$3 is the current opening and $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$4 is the maximal previous opening (Bonacini et al., 2020). This permits different loading and unloading responses while remaining within a variational framework. The quasi-static evolution is characterized by global stability and exact energy balance, with $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$5 increasing and $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$6 nondecreasing.

A related compactness issue arises in one-dimensional phase-field damage models with cohesive interface. There, the physical maximal opening is replaced by a fictitious variable $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$7 in order to prove existence of energetic evolutions, and a subsequent argument based on temporal regularity recovers the equivalence between the fictitious and real variables under general loading–unloading regimes (Bonetti et al., 2020). This shows that cohesive phase-field models often require an explicit treatment of history variables if unloading hysteresis is part of the constitutive target.

4. Reconstruction of prescribed cohesive laws and tunable strength surfaces

A major development in recent work is the inversion of the phase-field-to-cohesive mapping. Instead of merely deriving $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$8 from a chosen phase-field model, one may prescribe a target cohesive law $$\mathcal{F}_\varepsilon(u,v,A):= \int_A \left( \varphi( f_\varepsilon^2(v))\, |u'|^2 + \frac{(1-v)}{4} + \varepsilon\, |v'|^2 \right)\,dx,$$9 and reconstruct phase-field ingredients that reproduce it (Alessi et al., 16 Jul 2025). In the second part of the 2025 three-part framework, the one-dimensional cell problem is reduced to a scalar minimization over the minimum phase-field value $(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$0, encoded in a function $(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$1, together with an auxiliary function $(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$2 that governs whether the minimizing profile is absolutely continuous or contains a jump. When $(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$3 is strictly decreasing, the resulting cohesive law is $(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$4, concave, and satisfies

$(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$5

which yields an Abel integral equation linking the target cohesive law to the damage potential and degradation functions (Alessi et al., 16 Jul 2025).

This reconstruction framework yields multiple phase-field models with the same macroscopic $(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$6 but different localized phase-field profiles. Explicit constructions are worked out for Dugdale, linear softening, bilinear softening, hyperbolic softening, quadratic hyperbolic softening, exponential softening, and logarithmic softening laws (Alessi et al., 16 Jul 2025). A plausible implication is that the macroscopic cohesive response and the microscopic diffused-crack structure are not uniquely coupled: the same traction–separation law may correspond to distinct regularized localization patterns.

A second line of development makes the strength surface itself an explicit constitutive object. In the variational model with reversible eigenstrain $(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$7 and scalar damage $(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$8, one chooses an initial elastic domain

$(u,v)\in H^1(\Omega,\mathbb{R}\times[0,1])$9

and defines its support function

$d\approx 0$00

The degraded eigenstrain potential is

$d\approx 0$01

This construction accommodates an arbitrary convex strength surface, independent of the regularization length scale, and naturally enforces a sharp non-interpenetration condition (Vicentini et al., 13 Jun 2025). In the AT2 case, the corresponding one-dimensional traction–separation law is

$d\approx 0$02

so the cohesive law is Barenblatt-type and independent of the phase-field regularization length (Vicentini et al., 13 Jun 2025).

The eigenstrain-based return-mapping formulation introduces tunable strength surfaces directly at the constitutive level. Two explicit criteria are considered: a non-smooth criterion with independent tensile and shear strengths,

$d\approx 0$03

and a smooth Drucker–Prager-like criterion with pressure-dependent strengthening under compression (Hageman, 23 Mar 2026). In that framework, crack nucleation is controlled by the strength potential $d\approx 0$04 and material strengths such as $d\approx 0$05 and $d\approx 0$06, whereas crack propagation is controlled by $d\approx 0$07. The model therefore decouples strength from toughness rather than identifying both with a single regularization length.

5. Computational formulations and implementation strategies

The computational literature on cohesive phase-field models is heterogeneous because the constitutive target varies. The main variants can be organized as follows.

Model class	Distinctive ingredient	Representative paper
Variational cohesive $d\approx 0$08-limit	$d\approx 0$09 from a cell problem	(Alessi et al., 16 Jul 2025)
Critical-point approximation	ODE reduction and midpoint crack selection in 1D	(Bonacini et al., 2023)
Eigenstrain-based cohesive PF	Local return mapping with fracture eigenstrains	(Hageman, 23 Mar 2026)
Interface cohesive PF	Energy-consistent $d\approx 0$10 across diffuse interfaces	(Chen et al., 2024)
Multi-cohesive anisotropic PF	Directional degradation with multiple cohesive lengths	(Lacave et al., 13 Oct 2025)

From a numerical viewpoint, one attractive feature of variational cohesive phase-field energies is that they are smooth functionals on Hilbert spaces, and common algorithms such as alternate minimization typically compute stationary points rather than global minimizers (Bonacini et al., 2023). This is one reason why the convergence of critical points, not only minimizers, has become important in the mathematical analysis.

The eigenstrain-based framework reformulates cohesive evolution as a local constitutive problem, analogous to plasticity, because the fracture eigenstrains require no spatial gradients (Hageman, 23 Mar 2026). At each integration point, one updates the stress and eigenstrain multipliers by return mapping, while the global problem still contains only displacement and phase-field degrees of freedom. Consistent tangent operators are derived for both the non-smooth and the Drucker–Prager-like criteria, and the model can therefore be inserted into standard finite-element solvers. The paper reports implementation in FEniCSx with Numba and provides open-source code (Hageman, 23 Mar 2026).

In heterogeneous solids with interfaces, cohesive phase-field fracture can be compared directly with intrinsic cohesive-zone elements and hybrid models. A comparative study of fracture in multiphase materials found that the cohesive phase-field method is in agreement with a hybrid model combining CPFM and CZM when the interface zone thickness is not excessively small, which suggests that CPFM may serve as a unified fracture approach for multiphase materials provided the interface zone thickness is comparable to that of the other phases (Koopas et al., 2023). The same study reports that CPFM is length-scale insensitive with respect to $d\approx 0$11, whereas standard brittle phase-field required finer meshes and exhibited less satisfactory behavior on the benchmark considered (Koopas et al., 2023).

Length-scale-insensitive interface formulations make the same point from another direction. In a diffuse-interface cohesive phase-field model, the local fracture energy is interpolated as

$d\approx 0$12

and $d\approx 0$13 is chosen so that the integrated fracture energy across the diffuse interface equals the sharp cohesive interface energy $d\approx 0$14 (Chen et al., 2024). This energy-consistency condition is what makes the model insensitive to both the regularized interface thickness $d\approx 0$15 and the regularized fracture surface thickness $d\approx 0$16. The same paper extends the formulation to thermo-chemo-mechanical fracture in lithium-ion battery materials (Chen et al., 2024).

6. Anisotropy, dynamics, scope, and unresolved issues

Cohesive phase-field modeling has expanded beyond isotropic scalar settings. One route introduces anisotropy through the fracture energy and strength as arbitrary functions of crack direction, while preserving the length-scale-insensitive constitutive structure of isotropic cohesive phase-field damage models (Rezaei et al., 2021). In this formulation, the crack direction is extracted from $d\approx 0$17, the toughness becomes $d\approx 0$18, and the damage-function parameter $d\approx 0$19 is set by

$d\approx 0$20

so that direction-dependent strength is introduced in addition to direction-dependent fracture energy (Rezaei et al., 2021). The reported advantage is that one can increase the mesh size and reduce computational time without severe change in predicted crack path and load–displacement curves. The same work also states that these models still lack to capture mode-dependent fracture properties, identifying a limitation rather than a closed theory (Rezaei et al., 2021).

A different anisotropic approach is the multi-cohesive model for orthotropic materials. There a single damage variable $d\approx 0$21 is combined with a tensorial degradation of the orthotropic stiffness,

$d\approx 0$22

where the scalar functions

$d\approx 0$23

introduce distinct cohesive lengths along principal material directions (Lacave et al., 13 Oct 2025). This yields direct formulas for directional strengths such as

$d\approx 0$24

thereby decoupling nucleation strength from propagation properties with only one scalar damage field (Lacave et al., 13 Oct 2025).

Dynamic cohesive phase-field fracture is another major branch. A dynamic extension of the cohesive gradient-damage formulation uses a regularized fracture energy linear in the damage field together with non-polynomial degradation functions, and derives the governing equations from macro- and microforce balance theories while accounting for irreversible microstructural changes (Geelen et al., 2018). Numerical examples in that work include crack branching, Kalthoff–Winkler experiments, and fragmentation, with an augmented Lagrangian staggered scheme used to enforce irreversibility (Geelen et al., 2018). The eigenstrain-based 2026 model also treats dynamic loading and reports that complex phenomena such as crack branching under dynamic loading are naturally captured (Hageman, 23 Mar 2026).

Despite these advances, several restrictions remain explicit in the current literature. The convergence of critical points in (Bonacini et al., 2023) is one-dimensional, and extension to higher dimensions is stated to be challenging. The anisotropic cohesive formulation of (Rezaei et al., 2021) argues that mode-dependent fracture properties are still not captured satisfactorily. The eigenstrain-based model of (Hageman, 23 Mar 2026) is formulated in small strain and linear elasticity, although extensions to finite strain, anisotropy, and multiphysics are identified as natural directions. The variational model with flexibly tunable strength surface (Vicentini et al., 13 Jun 2025) requires a balance between the regularization length and the cohesive length: the ratio must be sufficiently small to ensure strain hardening, but also large enough to destabilize the homogeneous damaged state if crack nucleation is desired.

These limitations clarify the present state of the subject. Cohesive phase-field modeling is no longer confined to ad hoc softening modifications of brittle phase-field fracture. It now encompasses rigorously derived $d\approx 0$25-limits, critical-point convergence, inverse construction of prescribed cohesive laws, explicit strength surfaces, interface-specific regularizations, anisotropic multi-cohesive constitutive structures, and dynamic return-mapping formulations. At the same time, the choice between these variants remains constitutively meaningful: it determines whether the model is intended primarily as a regularization of a cohesive free-discontinuity functional, as a constitutive surrogate for a traction–separation law, or as a computational framework for crack nucleation, propagation, and branching under complex loading.