1 Introduction↩︎

Many real-world structures undergo fracture under dynamic loading, either deliberately (e.g.mining, fracking) or accidentally (e.g.crashes). For fracture under high strain rates, complex mechanisms emerge, including limiting crack-tip speeds [1], [2], crack-tip instabilities [3]–[6], and (micro-)branching phenomena [7]–[9]. The phase-field approach naturally lends itself to capturing intricate crack patterns without explicit branching or propagation criteria, tracking algorithms, or re-meshing. Originally proposed in the quasi-static regime as a regularization of Griffith’s crack propagation criterion [10] and later shown to amount to a gradient damage model [11], the phase-field approach to brittle fracture [12] formulates crack evolution as the minimization of a total energy functional [13]. Extensions to dynamics [14]–[16] retain the same regularization and supplement it with inertia effects; their governing equations can be derived from an extension of Hamilton’s principle to systems with dissipation [17]–[19].

Dynamic phase-field fracture models have reproduced at least qualitatively a wide range of experimentally observed phenomena, including the dependence of the energy release rate on crack-tip velocity [20], velocity toughening [21], a limiting mode-I crack speed of approximately 60% of the Rayleigh wave speed [22], inter- and transonic propagation [23], oscillatory crack-tip instabilities [24], micro- and macro-branching [25], and fragmentation [26], [27]. Substantial work has addressed mathematical aspects such as existence of solutions and convergence [28], stability [19], as well as numerical aspects [29]–[34].

However, several works have reported anomalous behavior in the interaction of elastic waves with regularized cracks. Steinke et al. [35] observe high-frequency oscillations following the interaction of an elastic wave with a pre-existing phase-field crack in a quasi-1D setting. They interpret these oscillations as a numerical artifact which can be alleviated through dissipative time integration and a crack region spanning at least one fully damaged element. Schlüter et al. [30] report similar oscillations in a 1D bar and note that, while some numerical solution approaches yield ‘no or minor oscillations’ [30], with all tested approaches the shape of the elastic wave is altered by the interaction with the pre-existing phase-field crack. More recently, Durussel et al. [36] interpret the phenomenon as a ‘trapping of elastic waves’ within the damaged zone, leading to spurious effects and a widening of the phase-field profile; they argue that these effects vanish as the phase-field regularization length goes to zero and classify them as numerical artifacts.

In classical models, the coupling between the displacement and the phase-field (or damage) variable is realized by multiplying the stiffness of the material by a degradation function depending on the phase-field variable. In the following, we refer to these models as phase-field models of brittle fracture with stiffness degradation. In dynamics, other models have been proposed which degrade both the stiffness and the mass density [24]. Their motivation is that the mismatch between degraded stiffness and undegraded density spuriously reduces the local wave speed within the damaged zone. Subsequent works in the large-strain setting [4], [37], [38] adopt this variant and report on its ability to capture higher crack-tip velocities and the experimentally observed oscillatory crack-tip instability. Tian et al. [38] further argue that, without density degradation, energy accumulates at the crack tip which may yield unrealistic crack-tip splitting. In this paper, we refer to these models as phase-field models of brittle fracture with stiffness+density degradation.

The phase-field models of cohesive fracture recently proposed in [39], [40] preserve the elastic stiffness and degrade only the material strength; hence, we denote them as phase-field models of cohesive fracture with strength degradation. These models have been proposed for quasi-static fracture and have not yet been extended to the dynamic regime. In this paper, we perform this extension and show that these models, besides their advantages in the quasi-static framework, are also promising alternatives for dynamics. In summary, this work makes the following contributions:

• For the classical model with stiffness degradation, we show that the high-frequency oscillations and the widening of the phase-field profile reported in the literature are not purely numerical artifacts, but consequences of the phase-field regularization in the dynamic setting. By a scattering analysis and computations with the transfer matrix method, we show that the wave-crack interaction is governed by the ratio of the regularization length to the wavelength. By finite element simulations, we explain the widening of the phase-field profile. We demonstrate that degrading the mass density in addition to the stiffness does not resolve the above issues.

• We extend the recently proposed phase-field model of cohesive fracture [39], [40] to the dynamic setting. We study under which conditions this new formulation solves the above issues and recovers the desired elastodynamic response. Under these conditions, we further derive analytically a dynamic cohesive zone law and identify the ratio of the wave speed to the Irwin cohesive length as parameter governing the temporal response.

The remainder of this work is organized as follows. Section 2 analyzes the phase-field models of dynamic brittle fracture in both variants with stiffness and stiffness+density degradation. Section 3 extends the cohesive phase-field formulation with strength degradation in [40] to the dynamic setting. In both sections, after presenting the governing equations in the multidimensional setting, we analyze the 1D problem of the interaction between an elastic wave and a phase-field crack both analytically and numerically. Section 4 compares the numerical results from all studied formulations on a 2D benchmark of dynamic fracture. Section 5 concludes with a summary and outlook.

2 Phase-field modeling of dynamic brittle fracture with stiffness and stiffness+density degradation↩︎

We start by analyzing phase-field models of dynamic brittle fracture based on stiffness degradation [18], [19], [41]. Then, we shift our attention to the model proposed by [24], [38] with stiffness and density degradation.

2.1 Brittle model with stiffness degradation↩︎

In this section, we first recapitulate the main aspects of phase-field models of dynamic brittle fracture with stiffness degradation. Subsequently, we analyze the interaction of elastic waves with a pre-existing phase-field crack in a 1D bar.

2.1.1 Energetic quantities↩︎

We consider a homogeneous, isotropic, linear elastic body occupying a portion of a \(d\)-dimensional space \(\Omega \subset \mathbb{R}^d\) with boundary \(\partial\Omega\). The state at a point \(\boldsymbol{x} \in \Omega\) and time \(t \in [0,T]\) is described by the displacement field \(\boldsymbol{u}(\boldsymbol{x}, t) : \Omega \times [0,T] \rightarrow \mathbb{R}^d\) and the phase (or damage) field \(\alpha(\boldsymbol{x}, t) : \Omega \times [0,T] \rightarrow [0,1]\), where \(\alpha=0\) denotes the pristine material and \(\alpha=1\) the fully damaged material. The infinitesimal strain tensor is defined as \(\boldsymbol{\varepsilon}(\boldsymbol{x}, t) = \nabla_{\text{sym}} \boldsymbol{u}(\boldsymbol{x}, t)\), with \(\nabla_{\text{sym}}\) as the symmetric spatial gradient operator. We denote with \(\boldsymbol{\bar u}(\boldsymbol{x},t)\) the prescribed displacement on the Dirichlet part of the boundary \(\partial\Omega_D\) and with \(\boldsymbol{f}(\boldsymbol{x},t)\) the prescribed traction on the Neumann part of the boundary \(\partial\Omega_N\). Unless stated otherwise, we assume the domain to be initially undamaged and at rest, namely \(\alpha(\boldsymbol{x}, 0)=0\) and \(\dot{\boldsymbol{u}}(\boldsymbol{x}, 0)=\boldsymbol{0}\) for all \(\boldsymbol{x} \in \Omega\), where \(\dot{(\bullet)}\) stands for the time derivative of \((\bullet)\).

The total potential energy of the body reads \[\label{eq:P} \mathcal{P} (\boldsymbol{u}, \alpha) = \mathcal{P}^\text{e} (\boldsymbol{u},\alpha) + \mathcal{P}^\text{f} (\alpha) - \mathcal{W} (\boldsymbol{u}) = \int_{\Omega} \psi (\boldsymbol{\varepsilon}(\boldsymbol{u}),\alpha) \,\mathrm{d}\boldsymbol{x} + \int_{\Omega} \frac{G_{\text{c}}}{c_w} \left( \frac{w(\alpha)}{\ell} + \ell |\nabla \alpha|^2 \right) \mathrm{d}\boldsymbol{x} - \int_{\partial \Omega_{\text{N}}} \boldsymbol{f} \cdot \boldsymbol{u} \,\mathrm{d}\boldsymbol{x}\tag{1}\] comprising the elastic energy \(\mathcal{P}^\text{e}\), the fracture energy \(\mathcal{P}^\text{f}\), and the work of the external forces \(\mathcal{W}\).

Following [42]–[46], the elastic strain energy density \(\psi\) is expressed as \[\label{eq:psi95brittle} \psi(\boldsymbol{\varepsilon}, \alpha) = g(\alpha) \psi_{\text{D}} (\boldsymbol{\varepsilon}) + \psi_{\text{R}} (\boldsymbol{\varepsilon}) \quad\text{,}\tag{2}\] where \(\psi_{\text{D}}\) is the degraded contribution driving damage evolution and \(\psi_{\text{R}}\) is the residual part. In this work, we employ the volumetric-deviatoric split [42], \[\label{eq:vol95dev95decomposition} \psi_{\text{D}} (\boldsymbol{\varepsilon}) = \frac{\kappa_0}{2} \langle \text{tr} (\boldsymbol{\varepsilon}) \rangle_+^2 + \mu_0 \left|\boldsymbol{\varepsilon}_{\text{dev}}\right|^2 \quad\text{,}\quad \psi_{\text{R}} (\boldsymbol{\varepsilon}) = \frac{\kappa_0}{2} \langle \text{tr} (\boldsymbol{\varepsilon}) \rangle_-^2 \quad\text{, with}\quad \psi_{\text{D}}(\boldsymbol{\varepsilon}) + \psi_{\text{R}}(\boldsymbol{\varepsilon}) = \psi_0 (\boldsymbol{\varepsilon})\quad\text{,}\quad\tag{3}\] where \(\kappa_0\) and \(\mu_0\) are the undamaged bulk and shear moduli, while \(\psi_0 (\boldsymbol{\varepsilon}) = \frac{\kappa_0}{2} \text{tr}(\boldsymbol{\varepsilon})^2 + \mu_0 \left|\boldsymbol{\varepsilon}_{\text{dev}}\right|^2\) is the undamaged elastic strain energy density. Also, \(\text{tr} (\boldsymbol{\varepsilon})\) is the trace of \(\boldsymbol{\varepsilon}\) and \(\boldsymbol{\varepsilon}_{\text{dev}}=\boldsymbol{\varepsilon}-\frac{1}{d}\text{tr} (\boldsymbol{\varepsilon})\boldsymbol{I}\) is its deviatoric part with \(\boldsymbol{I}\) as the second-order identity tensor, while we use \(\langle(\bullet)\rangle_\pm = \tfrac{1}{2} ((\bullet) \pm |(\bullet)|)\). The degradation function \(g(\alpha) = (1-\alpha)^2 + g_0\) couples the phase field and the displacement field and contains the constant \(0 < g_0 \ll 1\) (\(g_0=o(\ell)\)) leading to a small residual stiffness in fully damaged conditions. In the quasi-static case, the residual stiffness is used to ensure strict convexity of the total energy functional in \(\boldsymbol{u}\). In the dynamic case, we retain it for the reasons clarified in Section 2.1.5.

The second contribution in 1 is the Ambrosio-Tortorelli (AT) regularization of the surface energy [47], with fracture toughness \(G_{\text{c}}\), regularization length \(\ell\), and normalization constant \(c_w = 4 \int_0^1 \sqrt{w(\beta)} \mathrm{d}\beta\) depending on the choice of \(w(\alpha)\). The two most common choices are the AT1 model with \(w(\alpha) = \alpha\) (\(c_w = 8/3\)) and the AT2 model with \(w(\alpha) = \alpha^2\) (\(c_w=2\)). In this paper, for the brittle case we adopt the AT1 model, since AT2 leads to a vanishing elastic domain. The implications of adopting the AT2 model in the dynamic setting are briefly discussed in Section 2.3.

Finally, the kinetic energy reads \[\label{eq:K} \mathcal{K} (\dot{\boldsymbol{u}}) = \int_\Omega \frac{1}{2} \rho_0 \left|\dot{\boldsymbol{u}}\right|^2 \mathrm{d}\boldsymbol{x}\tag{4}\] with the mass density \(\rho_0\).

2.1.2 Governing equations↩︎

Following [18], [19], [41], the governing equations are obtained imposing the principle of stationary action, the energy balance, and the irreversibility condition. The latter reads \[\label{eq:irr} \dot{\alpha} \geq 0 \qquad\text{.}\tag{5}\] Along with the initial condition \(\alpha(\boldsymbol{x}, 0) \geq 0\), 5 ensures the non-negativity of the damage parameter \(\alpha \geq 0\). The space-time action functional between two time instants \(t_2>t_1\) is defined as \[\mathcal{A}(\boldsymbol{u}, \alpha) = \int_{t_1}^{t_2} \mathcal{L} (\boldsymbol{u}, \dot{\boldsymbol{u}}, \alpha) \,\mathrm{d}t \qquad\text{with}\qquad \mathcal{L} (\boldsymbol{u}, \boldsymbol{v}, \alpha) = \mathcal{P} (\boldsymbol{u}, \alpha) - \mathcal{K} (\boldsymbol{v})\,,\] where \(\mathcal{L}\) is the Lagrangian of the system. The principle of stationary action requires the first variation of \(\mathcal{A}\) to be non-negative for all admissible variations of the state variables over any interval \([t_1, t_2]\) [17], [48], namely \[\label{eq:PSA} \mathcal{A}^\prime (\boldsymbol{z}) (\hat{\boldsymbol{z}} - \boldsymbol{z}) \geq 0\,, \qquad \forall \hat{\boldsymbol{z}} \in \mathcal{Z}\,,\tag{6}\] where \(\boldsymbol{z} = (\boldsymbol{u}, \alpha)\) is the state vector, \(\mathcal{A}^\prime (\boldsymbol{z}) (\hat{\boldsymbol{z}} - \boldsymbol{z})\) denotes the Gâteaux derivative of \(\mathcal{A}\) in the direction \(\hat{\boldsymbol{z}} - \boldsymbol{z}\), and \(\mathcal{Z}\) is a sufficiently regular functional space incorporating the Dirichlet boundary conditions. Finally, the energy balance requires conservation of the total energy, \[\label{eq:EB} \frac{\mathrm{d}}{\mathrm{d}t} \mathcal{E} (\boldsymbol{u}, \dot{\boldsymbol{u}}, \alpha) = - \int_{\partial \Omega_{\text{N}}}\!\! \dot{\boldsymbol{f}} \cdot \boldsymbol{u} \,\mathrm{d}\boldsymbol{x} + \int_{\partial \Omega_{\text{D}}}\!\! \boldsymbol{\sigma} \boldsymbol{n} \cdot \dot{\bar{\boldsymbol{u}}} \,\mathrm{d}\boldsymbol{x} \qquad\text{with}\qquad \mathcal{E} (\boldsymbol{u}, \dot{\boldsymbol{u}}, \alpha) = \mathcal{P} (\boldsymbol{u},\alpha) + \mathcal{K} (\dot{\boldsymbol{u}}) \qquad\text{.}\tag{7}\]

Defining the stress \[\label{eq:sigma95brittle} \boldsymbol{\sigma} (\boldsymbol{\varepsilon}, \alpha) = \frac{\partial \psi (\boldsymbol{\varepsilon}, \alpha)}{\partial \boldsymbol{\varepsilon}} = g(\alpha) \left[ \kappa_0 \langle \text{tr} (\boldsymbol{\varepsilon}) \rangle_+ \boldsymbol{I} + 2 \mu_0 \boldsymbol{\varepsilon}_{\text{dev}} \right] + \kappa_0 \langle \text{tr} (\boldsymbol{\varepsilon}) \rangle_- \boldsymbol{I}\tag{8}\] and the damage energy release rate \[\label{eq:standard95err} Y (\boldsymbol{\varepsilon}, \alpha) = - \frac{\partial \psi (\boldsymbol{\varepsilon}, \alpha)}{\partial \alpha} = - g^\prime(\alpha) \psi_{\text{D}} (\boldsymbol{\varepsilon})\qquad\text{,}\tag{9}\] 5 , 6 and 7 yield the dynamic equilibrium equations \[\tag{10} \begin{equation}\tag{11} - \nabla \cdot \boldsymbol{\sigma} + \rho_0 \ddot{\boldsymbol{u}} = \boldsymbol{0}\,, \qquad \forall (\boldsymbol{x},t) \in \Omega \times [0,T]\,, \end{equation} \begin{equation}\tag{12} \boldsymbol{\sigma} \boldsymbol{n} = \boldsymbol{f}\,, \qquad \forall (\boldsymbol{x},t) \in \partial \Omega_{\text{N}} \times [0,T]\,,\\ \end{equation}\] which differ from those of the quasi-static case only by the presence of the inertia term \(\rho_0 \ddot{\boldsymbol{u}}\), and the Karush-Kuhn-Tucker (KKT) conditions for damage evolution \[\tag{13} \begin{equation} - Y (\boldsymbol{z}) + \frac{G_{\text{c}}}{c_w} \left( \frac{w^\prime(\alpha)}{\ell} - 2 \ell \Delta \alpha \right)\geq 0\,, \dot{\alpha} \geq 0\,, \left[ - Y (\boldsymbol{z}) + \frac{G_{\text{c}}}{c_w} \left( \frac{w^\prime(\alpha)}{\ell} - 2 \ell \Delta \alpha \right) \right] \dot{\alpha} = 0 \,\,\forall (\boldsymbol{x},t) \in \Omega \times [0,T] \tag{14} \end{equation} \begin{equation} \nabla \alpha \cdot \boldsymbol{n} \geq 0 \qquad \dot{\alpha} \geq 0 \qquad (\nabla \alpha \cdot \boldsymbol{n}) \dot{\alpha} = 0 \qquad \forall (\boldsymbol{x},t) \in \partial \Omega \times [0,T] \qquad\text{,} \tag{15} \end{equation}\] which coincide with those of the quasi-static case.

In the quasi-static case, the 1D counterpart of 14 written for the local damage model (i.e., for \(\alpha^\prime \equiv 0\)) delivers the 1D elastic domain [11] which, for the AT1 model, reads \[\label{eq:elastic95domain95brittle951D} \mathcal{S} (\alpha) = (-\infty, \hat{\sigma}_{\text{c}}(\alpha)] \qquad\text{with}\qquad \hat{\sigma}_{\text{c}}(\alpha) = \left( (1-\alpha)^2 +g_0 \right)\sqrt{\frac{3 G_{\text{c}} E_0}{8 \ell (1-\alpha)}} \qquad\text{.}\tag{16}\] The critical stress at the initial undamaged state \[\label{eq:tens95str} \hat{\sigma}_{\text{c}}(0) = \left(1+g_0 \right)\sqrt{\frac{3 G_{\text{c}} E_0}{8 \ell}}\tag{17}\] coincides with the tensile strength of the material.

2.1.3 Interaction of a tensile wave with a fixed phase-field crack in a 1D bar↩︎

We now consider a 1D bar which contains a phase-field crack at its center \(x=0\), represented by the optimal AT1 phase-field profile [49] \[\label{eq:profile95AT1} \alpha (x) = \begin{cases} \left( 1 - \frac{|x|}{2\ell} \right)^2 &x \in (-2\ell, 2\ell)\\ 0 &\text{else}\\ \end{cases} \qquad\text{.}\tag{18}\] In this section, this phase-field profile is held constant in time and denoted as a fixed phase-field crack (whereas in later sections we will let it evolve and denote it as an evolving phase-field crack). Hence, the damage field is known and only the elastodynamic problem 10 with the appropriate boundary and initial conditions needs to be solved. For simplicity, when studying the interaction of a tensile wave with a fixed phase-field crack, we ignore the volumetric-deviatoric energy decomposition and apply the degradation function to the whole elastic energy density, i.e. we adopt \(\psi_{\text{D}} (\varepsilon) = \psi(\varepsilon)\) and \(\psi_{\text{R}} (\varepsilon) = 0\). Combining the kinematics \(\varepsilon = u^\prime\) and the constitutive law \(\sigma = g(\alpha) E_0 \varepsilon\) with the balance of linear momentum 10 , we obtain the governing equation \[\label{eq:wave95equation95standard951D95nosplit} - \frac{\partial}{\partial x} \left[ g(\alpha) E_0 u^\prime \right] + \rho_0 \ddot{u} = 0\tag{19}\] which is the wave equation in a 1D medium with smoothly varying material properties.

The so-called scattering analysis, i.e. the analysis of the frequency-dependent reflection and transmission of an elastic wave in a heterogeneous elastic medium, is a classical subject [50], [51] and its main concepts and results are summarized in Appendix 5.1. As noted in [24], [36], [38], the spatially variable stiffness leads to a wave speed \(c(\alpha)= \sqrt{g(\alpha)E_0/\rho_0}=\sqrt{g(\alpha)}c_0\) which decreases from \(c_0\) in the intact material to a very small value (depending on \(g_0\)) for \(\alpha=1\) (Fig. 1). As a consequence, in absence of a residual stiffness no mechanical signal reaching the crack can propagate away from it. For this reason, we argue that \(g_0 > 0\) is needed not only in the quasi-static but also in the dynamic brittle case. However, in contrast to the claims in literature [24], [36], [38], the reduction of the local wave speed alone does not characterize the wave-crack interaction. As shown in Appendix 5.1, the full characterization of the scattering behavior is given by two profiles [51], [52]: the local wave speed \(c(\alpha)\), and the acoustic impedance \[\label{eq:impedance} Z(\alpha) = \sqrt{g(\alpha)E_0\,\rho_0} = \sqrt{g(\alpha)}\,Z_0 \quad\text{with}\quad Z_0 = \sqrt{E_0\rho_0} \quad\text{,}\tag{20}\] both reported in Fig. 1. In the model with stiffness degradation, both profiles degrade simultaneously, so that the regularized crack acts on the wave not only as a region of reduced speed, but also as a smooth impedance well. According to the Wentzel–Kramers–Brillouin (WKB) criterion for smoothly graded media [51], [52], a harmonic wave of angular frequency \(\omega\) traverses a heterogeneity without appreciable reflection wherever the relative impedance variation per local wavelength is small, i.e., wherever \[\label{eq:WKB} \delta(x) = \frac{|Z'(x)|}{Z(x)} \frac{c(x)}{\omega} \ll 1 \quad\text{,}\tag{21}\] while values of \(\delta\) of order one mark the onset of reflection. Note that, with a little abuse of notation, we are denoting \(c(x)=c(\alpha(x))\) as the composition of \(c(\alpha)\) with the phase-field profile 18 , and similarly for \(Z(x)\). For the optimal AT1 profile 18 we have \(1-\alpha(x) \simeq |x|/\ell\), therefore \(\sqrt{g(\alpha(x))} \simeq |x|/\ell\) to leading order near the center of the crack. The relative impedance gradient then diverges toward the center of the crack, \(|Z'|/Z \simeq 2/|x|\), but the local wavelength shrinks at exactly the compensating rate, \(c(x)/\omega \simeq (\lambda/2\pi)\,|x|/\ell\), so that \[\label{eq:WKB-AT1} \delta \simeq \frac{\lambda}{\pi\ell}\tag{22}\] to leading order near the crack center. Hence, the wave-crack interaction is controlled by the single global parameter \(\ell/\lambda\), to which \(\delta\) is inversely proportional: a small \(\ell/\lambda\) corresponds to the non-adiabatic regime \(\delta \gg 1\), in which the wave is reflected, while a large \(\ell/\lambda\) corresponds to the adiabatic regime \(\delta \ll 1\), in which the wave is transmitted as through a homogeneous medium; a transition region separates the two. Note that this analysis identifies the governing parameter and the two limit regimes, but not the quantitative location of the transition.

A quantitative analysis can be performed using the transfer matrix method (TMM), which approximates a heterogeneous medium with a stack of homogeneous layers with constant material properties (Fig. 2). Its main features are recalled in Appendix 5.2. The TMM makes the separation of roles between wave speed and impedance explicit. The global transfer matrix 65 is assembled from two types of factors with disjoint dependencies: the discontinuity matrices 63 depend only on the impedance ratios of adjacent layers, while the propagation matrices 64 depend only on the local wavenumber \(k(\alpha) = \omega/c(\alpha)\), i.e., on the local wave speed, and accumulate phase without generating reflection. Consequently, the reflected power is produced solely by the impedance profile, whereas the wave speed enters the scattering problem only indirectly, by dictating the local wavelength against which the impedance variation is measured. The wave speed reduction within the phase-field support is thus not per se the origin of the anomalous reflection and transmission; it is the accompanying impedance variation that governs them. Fig. 2a shows the obtained reflection and transmission coefficients in dependence of \(\ell/\lambda\) for the optimal AT1 profile with \(g_0 = 10^{-6}\). Three distinct regimes can be distinguished. If the wavelength is sufficiently larger than the phase-field support, in our case for \(\ell/\lambda \lessapprox 0.08\) (at this ratio, the reflection power coefficient first drops below \(99\)%), the quantitative behavior resembles that of a sharp crack, leading to a full reflection. The opposite happens when the wavelength is smaller than (or comparable to) the phase-field support, in our setup for \(\ell/\lambda \gtrapprox 0.22\)¹ (from then onwards, the reflection coefficient is below \(1\)%). In this case the behavior is similar to that of a homogeneous material and the wave is fully transmitted. In between these two extremes there is a regime where the wave is partially transmitted and partially reflected. The TMM results thus suggest that the brittle model with stiffness degradation yields the intended wave-crack interaction only for \(\lambda \gtrapprox \ell/0.08\). This condition is difficult to enforce in realistic scenarios where crack propagation can induce high-frequency waves. Moreover, as shown by 17 , in the present model the regularization length plays the role of a material parameter if the tensile strength of the material has to be quantitatively reproduced.

2.1.4 Interaction of a compressive wave with a fixed phase-field crack in a 1D bar↩︎

For a compressive wave, the expected full transmission (with both crack faces in contact) is recovered in the 1D case through the strain energy decomposition. With \(\psi_{\text{D}} (\varepsilon) = \tfrac{1}{2} E_0 \langle \varepsilon \rangle_+^2\) and \(\psi_{\text{R}} (\varepsilon) = \tfrac{1}{2} E_0 \langle \varepsilon \rangle_-^2\), the balance of linear momentum becomes \[\label{eq:wave95equation95standard951D95split} - g^\prime(\alpha) \alpha^\prime E_0 \langle u^\prime \rangle_+ - g(\alpha) E_0 H (u^\prime) u^{\prime\prime} - E_0 H (- u^\prime) u^{\prime\prime} + \rho_0 \ddot{u} = 0\tag{23}\] which collapses to the bulk wave equation for \(u^\prime < 0\). In 23 , \(H(\bullet)\) is the Heaviside function that takes the value \(1\) if \((\bullet) \geq 0\) and \(0\) if \((\bullet) < 0\).

2.1.5 Numerical examples in 1D↩︎

As follows, we report the results of a numerical investigation in which a sample pulse interacts with a phase-field crack.

2.1.5.1 Problem setup

To study the interaction between phase-field regularization and elastodynamics, we consider a 1D bar of length \(L\) (Fig. 3a) subjected to an imposed displacement \[\label{eq:1D95pulse95loading} \bar{u} (t) = \begin{cases} \frac{\tilde{\sigma}}{\rho_0 c_0} \frac{\tilde{T}}{2\pi} \left[ \cos \left( 2 \pi \frac{t}{\tilde{T}} \right) - 1 \right] &t \in \left[0, \tfrac{\tilde{T}}{2} \right]\\ - \frac{\tilde{\sigma}}{\rho_0 c_0} \frac{\tilde{T}}{\pi} &\text{else}\\ \end{cases} \qquad\text{with}\qquad \tilde{T} = \frac{\tilde{\lambda}}{c_0}\tag{24}\] at \(x=-L/2\) (Fig. 3b). This generates a sinusoidal tensile stress half-wave with wavelength \(\tilde{\lambda}\), full period \(\tilde{T}\) (or angular frequency \(\tilde{\omega}=2\pi/\tilde{T}\)) and amplitude \(\tilde{\sigma}\) traveling at a wave speed \(c_0 = \sqrt{E_0 / \rho_0}\) in the undamaged material, with \(E_0\) as the undamaged Young’s modulus (Fig. 3c). Note that, as a truncated sinusoid, the pulse is not monochromatic; we therefore distinguish its nominal wavelength \(\tilde{\lambda}\) from the wavelength \(\lambda\) of a generic harmonic component. The results expected from the interaction of this stress pulse with a sharp crack are shown in Fig. 4 by means of space-time diagrams (4a,c) and snapshots of the stress field at various time instants (4b,d). The first column shows the result for a tensile wave, where the crack acts as a free boundary, reflecting the wave as undistorted compressive wave without energy transmission to the right part of the bar. The second column shows the results for a compressive wave, which is transmitted without any distortion to the other side of the crack, as the crack lips are in contact.

In the numerical simulations, the problem is solved by the finite element method (FEM) with linear elements and a mesh size \(\Delta x \approx \ell/5\) (unless specified otherwise), known to be sufficiently fine in both quasi-static [53], [54] and dynamic [22], [27] cases. For time integration, we use the implicit Newmark-\(\beta\) scheme [55] with \(c_0\Delta t/\Delta x = 0.1\), \(\gamma = 1/2\), and \(\beta = 1/4\), providing unconditional stability, second-order accuracy, and no numerical dissipation in the undamaged elastodynamic case. We stop the computations as soon as the wave has reached the free end of the bar \(x= L/2\), hence there are no boundary reflections. Details of our numerical implementation are reported in Appendix 5.3.

In the following analyses, we adopt for the bar \(L=1000\) mm, for the material \(E_0 = 30\,000\) MPa, \(G_{\text{c}} = 0.1\) N/mm, and \(\rho_0 = 2\,400\) kg/m\(^3\), and for the imposed stress wave \(\tilde{\sigma} = 4.5\) MPa and \(\tilde{\lambda}=600\) mm. We test two values of the regularization length, namely \(\ell = 15\) mm and \(\ell = 45\) mm; both ensure a small phase-field support compared to domain size (\(\ell/L=0.015\) and \(0.045\)) and wavelength of the stress wave (\(\ell/\tilde{\lambda} = 0.025\) and \(0.075\)). With the above values, it is \(\hat{\sigma}_{\text{c}}(0)=8.7\) MPa and \(5\) MPa, hence \(\tilde{\sigma} / \hat{\sigma}_{\text{c}}(0) = 0.52\) and \(0.9\). Finally, we set \(g_0 = 10^{-6}\).

2.1.5.2 Numerical results for a tensile wave and a fixed phase-field crack

Fig. 5 reports a space-time diagram of the stress (Fig. 5a-c), the fixed phase-field profile (Fig. 5d-f), and snapshots of the stress field at selected time steps (Fig. 5g-i). The left column (Fig. 5b,d,g) shows the results for \(\ell/\tilde{\lambda} = 0.075\). Once the wave enters the region of the phase-field crack, its speed is reduced and its shape is distorted. Upon arriving at the center of the crack at \(x=0\), part of the wave is reflected with an altered shape and high-frequency oscillations, while the remaining portion is transmitted to the right part of the bar. In [30], [35] this behavior is interpreted as a numerical artifact due to a coarse spatial and temporal discretization. To show that this is not the only reason, we illustrate in the middle column (Fig. 5b,e,h) the results obtained with a significantly finer discretization of \(\Delta x \approx \ell/50\) (keeping \(c_0\Delta t/\Delta x = 0.1\)). There, the high-frequency distortions of the wave are even more pronounced, and part of the wave energy is still transmitted through the crack. In the last row, we overlay the results with further ones obtained using the generalized-\(\alpha\) time integration scheme with parameters \(\alpha_m = 0.2\), \(\alpha_f = 0.4\), \(\gamma = \tfrac{1}{2} - \alpha_m + \alpha_f\) and \(\beta= \tfrac{1}{4} (1 - \alpha_m + \alpha_f)^2\). These parameters correspond to a high-frequency spectral radius \(\rho_\infty = 2/3\), so as to introduce controlled numerical dissipation of high-frequency modes while retaining second-order accuracy and unconditional stability. While the numerical damping slightly lowers the amplitudes of the oscillations, the same behavior obtained with Newmark-\(\beta\) persists. From now on, we exclusively employ the generalized-\(\alpha\) scheme with the mentioned parameters. Finally, the last column of Fig. 5 (c,f,i) shows the results for a smaller ratio \(\ell/\tilde{\lambda} = 0.025\), but keeping \(\Delta x \approx \ell/50\). Confirming the theoretical results, the interaction of the elastic wave and the phase-field crack resembles more closely that of the sharp crack (Fig. 4), while the shape of the reflected wave is still significantly altered.

2.1.5.3 Numerical results for a compressive wave and a fixed phase-field crack

In the FEM simulations a compressive wave is transmitted with no distortion, as shown in Fig. 6, matching the expected behavior (Fig. 4b,d).

2.1.5.4 Numerical results for a tensile wave and an evolving phase-field crack

So far we held the phase field fixed (\(\dot{\alpha} = 0\)). Next, we repeat the computation while letting both displacement and phase field evolve. We retain the setup of Fig. 3 with \(\ell/\tilde{\lambda} = 0.075\), \(\ell/\Delta x \approx 50\), and \(c_0 \Delta t / \Delta x = 0.1\), and adopt the decomposition of the strain energy. The results, illustrated in Fig. 7, show a widening of the phase-field profile as soon as the wave interacts with the regularized crack, indicating the fulfillment of the damage evolution criterion 13 and a spurious increase of the dissipated energy that no longer matches the fracture toughness \(G_{\text{c}}\).

This phenomenon can be explained considering the 1D elastic domain 16 , from which the profile of the material strength along the bar is reported in Fig. 7.

As the wave enters the phase-field support, damage evolution is triggered since the stress exceeds the local \(\hat{\sigma}_{\text{c}}(\alpha)\), which becomes extremely small (depending on \(g_0\)) at the center of the phase-field crack (Fig. 7c,d). As the damage reduces \(\hat{\sigma}_{\text{c}}(\alpha)\), further points are exposed to damage evolution (Fig. 7e-h), which widens the damaged band. In addition, the local reduction of the wave speed forces the incoming stresses to pile up as the wave stalls, further facilitating damage. The phase-field evolution along with the wave–stiffness interaction discussed earlier involving partial wave reflection ultimately leads to high-frequency oscillations as illustrated in Fig. 7i.

Note that, for a fixed phase field, the reflection and transmission behavior quantified with the TMM is independent of the stress wave amplitude. Instead, the widening behavior (as well as the accompanying reflection and transmission) does depend on the amplitude of the stress, as this amplitude directly influences the attainment of the local material strength.

2.2 Brittle model with stiffness+density degradation↩︎

To avoid a spatially varying wave speed, Chen et al. [24] and Tian et al. [38] propose to modify the kinetic energy by degrading also the mass density. In this section, we first briefly recapitulate the main equations of this model and then we again analyze the interaction of elastic waves with a pre-existing phase-field crack in a 1D bar as predicted by this model.

2.2.1 Energetic quantities and governing equations↩︎

With the brittle model with stiffness+density degradation, the kinetic energy reads \[\label{eq:kinetic95energy95degraded} \mathcal{K} (\dot{\boldsymbol{u}}, \alpha) = \int_\Omega \frac{1}{2} h(\alpha) \rho_0 \left|\dot{\boldsymbol{u}}\right|^2 \mathrm{d}\boldsymbol{x} \quad\text{,}\tag{25}\] where \(h(\alpha)\) is a monotonically decreasing density degradation function, which, along with 5 , leads to the violation of the mass conservation condition, namely \[\dot{m} = \frac{\mathrm{d}}{\mathrm{d} t} \int_{\Omega} h(\alpha) \rho_0 \mathrm{d}\boldsymbol{x} = \int_{\Omega} h^\prime(\alpha) \dot{\alpha} \rho_0 \mathrm{d}\boldsymbol{x} \leq 0 \qquad\text{for}\qquad \dot{\alpha} \geq 0 \qquad\text{.}\]

By reapplying the principles in Section 2.1.2 with 25 in place of 4 , the previous governing equation 10 becomes \[\label{eq:BF95strong95form95u95degraded} - \nabla \cdot \boldsymbol{\sigma} + \underbrace{h^\prime(\alpha) \dot{\alpha} \rho_0 \dot{\boldsymbol{u}}}_{(\dagger)} + h(\alpha) \rho_0 \ddot{\boldsymbol{u}} = \boldsymbol{0} \qquad \forall (\boldsymbol{x},t) \in \Omega \times [0,T]\quad,\tag{26}\] still accompanied by boundary conditions 12 , where \((\dagger)\) is an additional term not present in the sharp crack case nor in the brittle model with stiffness degradation. The KKT conditions 13 for the phase field remain unchanged, but with a modified damage energy release rate \[Y(\boldsymbol{\varepsilon}, \dot{\boldsymbol{u}}, \alpha) = - g^\prime(\alpha) \psi_{\text{D}} (\boldsymbol{\varepsilon}) + \underbrace{\frac{1}{2} h^\prime(\alpha) \rho_0 \left|\dot{\boldsymbol{u}}\right|^2}_{(\diamond)} \qquad\text{,}\] which now includes the contribution \((\diamond)\) depending on the velocity \(\dot{\boldsymbol{u}}\). Since \(h^\prime(\alpha) \leq 0\) and \(\rho_0 \left|\dot{\boldsymbol{u}}\right|^2 \geq 0\), this term leads to a lower energy release rate compared to 9 .

In the 1D case, the critical stress reads \[\label{eq:elastic95domain95brittle951D95densdegrad} \hat{\sigma}_{\text{c}}(\alpha,\dot{u}) = \left( (1-\alpha)^2 +g_0 \right) \sqrt{\frac{3 G_{\text{c}}E_0}{8 \ell (1-\alpha)} + E_0 \rho_0 \dot{u}^2 } \qquad\text{,}\tag{27}\] i.e. it now depends also on \(\dot{u}\) and grows for increasing velocities.

2.2.2 Interaction of a tensile wave with a fixed phase-field crack in a 1D bar↩︎

We now consider the same problem of Section 2.1.3, i.e. we want to study in 1D the interaction of a tensile wave with a phase-field crack represented by the optimal profile 18 which is held constant in time. As in Section 2.1.3, we ignore the volumetric-deviatoric decomposition. The equation of motion becomes \[\label{eq:wave95equation95density} - \frac{\partial}{\partial x} \left[ g(\alpha) E_0 u^\prime \right] + \frac{\partial}{\partial t}\left[h(\alpha)\rho_0 \dot{u}\right] = 0\qquad .\tag{28}\] For a fixed phase field (\(\dot{\alpha} = 0\)), the elastodynamic problem reduces again to a heterogeneous wave equation, now with independently degraded stiffness \(E(\alpha) = g(\alpha)E_0\) and density \(\rho(\alpha) = h(\alpha)\rho_0\). Thus, the scattering analysis in Appendix 5.1 continues to apply, now with \[\label{eq:wavespeed95density95deg} c (\alpha) = \sqrt{\frac{g(\alpha)}{h(\alpha)}}c_0\qquad Z(\alpha) = \sqrt{g(\alpha)\,h(\alpha)}Z_0\tag{29}\] The choice \(h(\alpha) = g(\alpha)\) leads to a constant wave speed \(c(\alpha)=c_0\), whereas the acoustic impedance becomes \[Z(\alpha) \;=\; g(\alpha)\,Z_0 \quad\text{.} \label{eq:impedance-dd}\tag{30}\] Thus, the density degradation does not remove the impedance well accompanying the regularized crack, but rather makes it stronger. Let us now analyze the adiabatic criterion 21 . With the wave speed equal to \(c_0\), the local wavenumber is \(k(\alpha) = \omega/c_0 = 2\pi/\lambda\). The local wavelength no longer contracts toward the crack center, and nothing compensates the diverging relative impedance gradient, hence \[\label{eq:WKB-dd} \delta(x) = \frac{|Z'(x)|}{Z(x)}\,\frac{c_0}{\omega} \simeq\; \frac{\lambda}{2\pi\,|x|}\, \xrightarrow[\;x\to 0\;]{}\; \infty \quad\text{.}\tag{31}\] Thus, a non-adiabatic region of size \(|x| \lesssim \lambda/(2\pi)\) exists at every frequency and its extent is set by the wavelength, not by \(\ell\). Across this region the impedance drops by several orders of magnitude, down to its residual value at the crack center. Total reflection therefore occurs independently of the \(\ell/\lambda\) ratio. In the TMM, the propagation matrices 64 accumulate a spatially uniform phase and the entire scattering response is produced by the discontinuity matrices 63 through the impedance ratios \(Z^{l+1}/Z^{l}\) (with \(l\) denoting the layer index), whose profile is reported in Fig. 8b. This behavior is shown in Fig. 8a, for which we set \(h(\alpha) = (1-\alpha)^2 + h_0\) with a residual density \(h_0\) and investigate two values of the residual density \(h_0\), one equal to \(g_0=10^{-6}\) and a much larger one \(h_0=10^{-2}\). While TMM results are unaffected by this choice, the reason for using a larger value will become apparent later. Unlike in the stiffness-degraded model, no adiabatic (transmissive) regime can be reached by tuning \(\ell/\lambda\).

2.2.3 Interaction of a compressive wave with a fixed phase-field crack in a 1D bar↩︎

To study the interaction of a compressive wave, we use the volumetric-deviatoric decomposition as in Section 2.1.4, which yields the governing equation \[\label{eq:wave95equation95stiffdens} \begin{align} &- g^\prime(\alpha) \alpha^\prime E_0 \langle u^\prime \rangle_+ - g(\alpha) E_0 H (u^\prime) u^{\prime\prime} - E_0 H (- u^\prime) u^{\prime\prime} + h(\alpha) \rho_0 \ddot{u} = 0 \quad\text{.} \end{align}\tag{32}\] In compression, the stiffness is not degraded but the density is, leading to \[\label{eq:wavespeed95density95deg95compr} c (\alpha) = \frac{c_0}{\sqrt{h(\alpha)}}\qquad Z(\alpha) = \sqrt{h(\alpha)}Z_0 \quad\text{,}\tag{33}\] thus the impedance still decreases within the phase-field support while the local wavelength now grows toward the crack center, giving \(\delta(x) \simeq \lambda\ell/(2\pi x^2)\) and hence an even more strongly non-adiabatic region.

2.2.4 Numerical results↩︎

As follows we report numerical results obtained with the same problem setup of Section 2.1.5, this time using the brittle model with stiffness and density degradation.

2.2.4.1 Numerical results for a tensile wave and a fixed phase-field crack

We repeat the numerical experiment of Fig. 5 for the model with stiffness+density degradation and report the results in Fig. 9. We set \(h_0 = 10^{-2}\), a value significantly larger than \(g_0 = 10^{-6}\). This choice stems from the fact that the degraded mass induces high frequency oscillations of the stress and numerical instabilities once the wave interacts with the phase-field profile. To limit oscillations, alternative strategies are available; however, these modify the material behavior by introducing viscosity in the constitutive law [24], [38] or give up variational consistency by neglecting the concave term in \(Y\) as in [36]. Apart from the numerical difficulties, another issue with this model is that the velocity-dependent term \(\tfrac{1}{2} h^\prime(\alpha) \rho_0 \left|\dot{\boldsymbol{u}}\right|^2\) in \(Y\) potentially leads to a negative energy release rate in case of high velocities. A negative energy release rate can lead to the loss of convexity of the damage subproblem as explained in [56].

Figs. 9a,d,g involve a tensile stress wave and use \(\ell/\Delta x \approx 5\), while Figs. 9b,e,h report the results with \(\ell/\Delta x \approx 50\). The wave is fully reflected (apart from minor numerical artifacts for \(\Delta x \approx \ell/5\) which are visible in Fig. 9g). Unlike with the model with stiffness degradation, here the results are independent of the discretization.

2.2.4.2 Numerical results for a compressive wave and a fixed phase-field crack

Let us now study the interaction of a compressive wave (\(|\tilde{\sigma}| = 4.5\) MPa and remaining parameters same as before) with a phase-field crack, see Figs. 9c,f,i for \(\ell/\Delta x \approx 50\). For a compressive wave, the combination of undegraded stiffness and degraded density results in an increase of the wave speed approaching \(\alpha=1\) (diverging to \(+\infty\) if \(h_0\) is neglected), see 33 . In the results of Fig. 9c,i, the transmitted wave is reduced in amplitude and altered in shape, while a significant reflection occurs at the regularized crack, highlighting that the strain energy decomposition does not yield the expected behavior. Moreover, in Fig. 9c a faster motion of the wave is visible around \(t/T \approx 0.6\) at the regularized crack, which reflects the increase of wave speed indicated by 33 in the compressive case. For the compressive wave, we obtained an even more strongly non-adiabatic region, consistent with the significant reflection and waveform alteration observed in Fig. 9c,i.

2.2.4.3 Numerical results for a tensile wave and an evolving phase-field crack

We now repeat the experiment of Fig. 7, where the phase field is allowed to evolve, for the model with stiffness+density degradation. The results in terms of \(x\)-\(t\) diagram are illustrated in Fig. 10a, while the phase-field profile and the stress field along with the material strength are reported in Fig. 10b and Fig. 10c-i, respectively. Unlike in Fig. 7, now the damage does not evolve, hence the phase-field profile does not widen and the tensile stress wave keeps being completely reflected. In Fig. 10c-i we plot the material strength given by 27 for \(\alpha=0\) along with its counterpart without the velocity-dependent contribution, 17 , to show that, for the case analyzed here, the dependence of the strength on the velocity is negligible.

2.3 Summary of the results with brittle models↩︎

From the results reported in Sections 2.1 and 2.2, we conclude that phase-field formulations of dynamic brittle fracture cannot consistently represent the interaction of stress waves with a pre-existing crack, regardless of whether they involve degrading the stiffness or both the stiffness and the mass density. Phase-field models with degradation functions different from the ones chosen here should display similar issues if the same degradation mechanisms and coupling between damage and mechanical fields are retained. This is the case for the models recovering cohesive-like behavior proposed in [57]–[59] for the quasi-static case and extended to dynamics in [27], [60], [61]. We therefore do not consider them separately here. Concerning the AT2 model, its vanishing elastic domain implies that any non-zero stress wave triggers local damage evolution and hence a local change of material properties, as well as a mix of reflections and transmissions, ultimately leading to a distorted wave regardless of the wave amplitude.

Hence, a different modeling approach is required. A potential candidate is the cohesive model proposed for the quasi-static case in [39], [40] where the material stiffness is not degraded. In the following section, we extend this model to dynamics and analyze its predictions.

3 Phase-field modeling of dynamic cohesive fracture with strength degradation↩︎

In this section, we extend the recently proposed phase-field regularization of cohesive fracture [39], [40] to dynamics. Following the structure of the previous section on brittle models, we first derive the governing equations and subsequently analyze the model behavior in 1D, focusing on the interaction of a stress wave with a pre-existing crack.

3.1 Energetic quantities and governing equations↩︎

The main modification introduced by the cohesive model in [39], [40] lies in the elastic strain energy density and the regularity of the fields. In particular, the cohesive formulation introduces an additional reversible kinematic field, the eigenstrain \(\boldsymbol{\eta} : \Omega \times [0,T] \rightarrow \mathbb{M}^{d}_{\text{sym}}\), where \(\mathbb{M}^{d}_{\text{sym}}\) is the set of symmetric second-order tensors. Specifically, the elastic strain energy density is written as² [40] \[\label{eq:psi95cohesive} \psi (\boldsymbol{\varepsilon}, \boldsymbol{\eta}, \alpha) = \frac{\kappa_0}{2} \left( \text{tr}(\boldsymbol{\varepsilon}) - \text{tr}(\boldsymbol{\eta}) \right)^2+ \mu_0 \left( \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| - \left| \boldsymbol{\eta}_{\text{dev}} \right| \right)^2+ \pi (\boldsymbol{\eta}, \alpha) \quad\text{,}\tag{34}\] where the eigenstrain potential \(\pi(\boldsymbol{\eta},\alpha)\) is defined as \[\label{eq:eigenstrain95potential} \pi(\boldsymbol{\eta},\alpha) = \begin{cases} a(\alpha)\pi_0(\text{tr}(\boldsymbol{\eta}), \left|\boldsymbol{\eta}_{\text{dev}}\right|) \qquad&\text{if}\;\text{tr}(\boldsymbol{\eta}) \geq 0\\ +\infty &\text{otherwise} \end{cases}\quad\text{.}\tag{35}\] Here, \(\pi_0\) is the support function of the initial elastic domain of the material \(\mathcal{S}_0\), which is considered a material property. The eigenstrain potential 35 is the support function of the damaged elastic domain \(\mathcal{S} (\alpha)\), which is assumed to homothetically shrink with increasing \(\alpha\) from \(\mathcal{S}_0=\mathcal{S} (0)\) through the degradation function \(a(\alpha)\). This means that in this model \(\alpha\) leads to the degradation of the local strength rather than of the stiffness. Further assumptions involve \(\mathcal{S} (\alpha)\) being convex and bounded in any direction apart from the purely hydrostatic compressive stresses, i.e. for \(\boldsymbol{\sigma}_{\text{dev}}=\boldsymbol{\sigma}-\tfrac{1}{d}\,\text{tr}(\boldsymbol{\sigma})\boldsymbol{I}=\boldsymbol{0}\) and \(\text{tr}(\boldsymbol{\sigma})<0\) [40]. Apart from these properties, the selection of \(\mathcal{S} (\alpha)\) is arbitrary and, for the case at hand, is specified in Section 3.2 along with the choice of \(a(\alpha)\).

Assuming sufficient temporal regularity and adopting the spatial regularity framework of the quasi-static setting [40], the state vector \(\boldsymbol{z} = (\boldsymbol{u}, \boldsymbol{\eta}, \alpha)\) at each time instant is piecewise smooth with a singular part localized on a jump set \(J(\boldsymbol{z}) \subset \Omega\) of co-dimension \(1\). The displacement \(\boldsymbol{u}\) is continuously differentiable on \(\Omega \setminus J(\boldsymbol{z})\) and admits jumps only across \(J(\boldsymbol{z})\), so the strain decomposes into a regular and a singular part as \[\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}_{\text{R}} + \boldsymbol{\varepsilon}_{\text{S}} \quad \text{, with} \quad \boldsymbol{\varepsilon}_{\text{R}} = \nabla_{\text{sym}} \boldsymbol{u} \quad \text{, and} \quad \boldsymbol{\varepsilon}_{\text{S}} = (\llbracket \boldsymbol{u} \rrbracket \otimes_{\text{sym}} \boldsymbol{m}) \, \delta_{J(\boldsymbol{z})} \qquad\text{,} \label{eq:CF95strain95jump}\tag{36}\] where \(\delta_{J(\boldsymbol{z})}\) is the Dirac surface measure concentrated on \(J(\boldsymbol{z})\), \(\otimes_{\text{sym}}\) denotes the symmetrized outer product and \(\llbracket \boldsymbol{u} \rrbracket(\boldsymbol{x},t) = \boldsymbol{u}^+(\boldsymbol{x},t) - \boldsymbol{u}^-(\boldsymbol{x},t)\) is the displacement jump across \(J(\boldsymbol{z})\), with \(\boldsymbol{u}^+(\boldsymbol{x},t)\) the limit of \(\boldsymbol{u}(\boldsymbol{x},t)\) approaching \(\boldsymbol{x}\) from the direction of the unit normal \(\boldsymbol{m}\) to the jump set at any given \(t\). Analogously to 36 , we have \(\boldsymbol{\eta} = \boldsymbol{\eta}_{\text{R}} + \boldsymbol{\eta}_{\text{S}}\) and, enforcing \(\boldsymbol{\varepsilon} - \boldsymbol{\eta} \in L^2(\Omega)\) to guarantee a finite energy, we obtain \[\boldsymbol{\eta}_{\text{S}} = \boldsymbol{\varepsilon}_{\text{S}}= \left( \llbracket \boldsymbol{u} \rrbracket \otimes_{\text{sym}} \boldsymbol{m} \right) \, \delta_{J(\boldsymbol{z})} \qquad\text{,} \label{eq:CF95eta95jump}\tag{37}\] which gives the link between the singular part of the eigenstrain and the displacement jump. The condition \(\text{tr}(\boldsymbol{\eta}) \geq 0\) in 35 naturally enforces non-interpenetration at the crack lips, \(\llbracket \boldsymbol{u} \rrbracket \cdot \boldsymbol{m} \geq 0\). For further details, we refer to [39], [40], [62].

In contrast to the brittle model with stiffness and density degradation of Section 2.2, the cohesive model does not degrade the mass, so the bulk density remains \(\rho_0\); the only question is whether the crack itself carries mass. Writing the total mass as \(m = \int_{\Omega \setminus J(\boldsymbol{z})} \rho_0 \, \mathrm{d}x + m_J\), where \(m_J\) denotes any mass concentrated on the jump set, mass conservation between an undamaged configuration (\(J(\boldsymbol{z}) = \emptyset\)) at time \(t_A\) and a cracked configuration at \(t_B\) requires \(m_J = 0\): the bulk integral is unchanged since the density is undegraded, so the added surface set \(J(\boldsymbol{z})\), of co-dimension 1, must carry no mass. Hence the density admits no singular part and the kinetic energy is \[\label{eq:kin95cohesive} \mathcal{K} (\dot{\boldsymbol{u}}) = \int_\Omega \frac{1}{2} \left|\dot{\boldsymbol{u}}\right|^2 \mathrm{d}m = \int_{\Omega\setminus J(\boldsymbol{z})} \frac{1}{2} \rho_0 \left|\dot{\boldsymbol{u}}\right|^2 \mathrm{d}\boldsymbol{x} \quad\text{.}\tag{38}\] The fracture energy remains identical to the brittle case of 1 .

The governing equations are obtained as usual from the principle of stationary action with irreversibility and energy balance. The momentum balance is obtained as \[\label{eq:u95strongform95cohesive} - \nabla \cdot \boldsymbol{\sigma} + \rho_0 \ddot{\boldsymbol{u}} = \boldsymbol{0} \qquad \forall (\boldsymbol{x},t) \in \Omega \times [0,T] \qquad\text{,}\quad \boldsymbol{\sigma} \cdot \boldsymbol{n} = \boldsymbol{f} \qquad \forall (\boldsymbol{x},t) \in \partial \Omega_{\text{N}} \times [0,T]\tag{39}\] with the cohesive stress tensor given by \[\label{eq:stress95cohesive95pressure95shear} \boldsymbol{\sigma} (\boldsymbol{\varepsilon}, \boldsymbol{\eta}, \alpha) =\frac{\partial \psi(\boldsymbol{\varepsilon}, \boldsymbol{\eta}, \alpha)}{\partial\boldsymbol{\varepsilon}}= p \boldsymbol{I} + \tau \frac{\boldsymbol{\varepsilon}_{\text{dev}}}{\left|\boldsymbol{\varepsilon}_{\text{dev}}\right|}\tag{40}\] where \[p = \frac{\partial \psi(\boldsymbol{\varepsilon}, \boldsymbol{\eta}, \alpha)}{\partial \text{tr}(\boldsymbol{\varepsilon})} = \kappa_0 \left( \text{tr}(\boldsymbol{\varepsilon}) - \text{tr}(\boldsymbol{\eta}) \right) \quad\text{,}\quad \tau = \frac{\partial \psi(\boldsymbol{\varepsilon}, \boldsymbol{\eta}, \alpha)}{\partial \left| \boldsymbol{\varepsilon}_{\text{dev}}\right|} = 2 \mu_0 \left( \left|\boldsymbol{\varepsilon}_{\text{dev}} \right| - \left| \boldsymbol{\eta}_{\text{dev}} \right| \right)\] are the hydrostatic and the shear stress, respectively.

The eigenstrain evolution criterion [40] reads \[\label{eq:eigenstrain95evolution95criterion} p \zeta + \tau \xi \leq a(\alpha) \pi_0^\prime (\text{tr}(\boldsymbol{\eta}), \left|\boldsymbol{\eta}_{\text{dev}}\right|) (\zeta, \xi)\tag{41}\] for any admissible variations \(\zeta\) and \(\xi\). Since \(\nabla\alpha\) may be discontinuous across \(J(\boldsymbol{z})\), for the phase field we obtain two sets of KKT conditions \[\begin{equation} -Y (\boldsymbol{\eta}, \alpha) + \frac{G_{\text{c}}}{c_w}\!\left(\!\frac{w^\prime(\alpha)}{\ell}\!-\!2 \ell \Delta \alpha\!\right) \geq 0 \text{,}\, \dot{\alpha} \geq 0 \text{,}\, \left[-Y (\boldsymbol{\eta}, \alpha) + \frac{G_{\text{c}}}{c_w}\!\left(\!\frac{w^\prime(\alpha)}{\ell}\!-\!2 \ell \Delta \alpha\!\right) \right] \dot{\alpha} = 0 \quad \forall (\boldsymbol{x},t) \in \left( \Omega \setminus J(\boldsymbol{z}) \right) \times [0,T] \,\text{,} \end{equation} \begin{equation} -Y (\boldsymbol{\eta}, \alpha) - 2 \ell \frac{G_{\text{c}}}{c_w} \llbracket \nabla \alpha \rrbracket \cdot \boldsymbol{n} \geq 0 \quad\text{,}\quad \dot{\alpha} \geq 0 \quad\text{,}\quad \left[-Y (\boldsymbol{\eta}, \alpha) - 2 \ell \frac{G_{\text{c}}}{c_w} \llbracket \nabla \alpha \rrbracket \cdot \boldsymbol{n} \right] \dot{\alpha} = 0 \quad \forall (\boldsymbol{x},t) \in J(\boldsymbol{z})\times [0,T] \,\text{,} \end{equation}\] in addition to 15 at the boundary. The bulk conditions mirror those of the brittle model 13 with an energy release rate \[Y (\boldsymbol{\eta}, \alpha) = -\frac{\partial \psi(\boldsymbol{\varepsilon}, \boldsymbol{\eta}, \alpha)}{\partial \alpha}= -\frac{\partial \pi(\boldsymbol{\eta}, \alpha)}{\partial \alpha}=-a^\prime(\alpha) \pi_0 (\text{tr}(\boldsymbol{\eta}), \left|\boldsymbol{\eta}_{\text{dev}}\right|) \qquad\text{.}\]

3.2 Adopted model and compact reformulation↩︎

We adopt \(w(\alpha)=\alpha^2\) and \(c_w=2\) as in [40], and the linear degradation function \(a(\alpha) = 1-\alpha\) from [39]. Note that no residual value is necessary for this model, i.e. we can have \(a(\alpha=1) = 0\). For \(\pi_0\) we adopt the \(2\)-norm model in [40] which reads \[\label{eq:strength95pot} \pi_0(\text{tr}(\boldsymbol{\eta}), \left|\boldsymbol{\eta}_{\text{dev}}\right|) = \sqrt{p_{\text{c}}^2 \text{tr}(\boldsymbol{\eta})^2 + \tau_{\text{c}}^2 \left| \boldsymbol{\eta}_{\text{dev}} \right|^2}\,,\tag{42}\] where \(p_{\text{c}}\) and \(\tau_{\text{c}}\) are the critical pressure and shear strength, respectively. As clarified later in Section 4, this leads to the same initial elastic limit as for the brittle case with the volumetric-deviatoric split 3 .

Before analyzing the model behavior, we introduce a more compact reformulation for the special case \(p_{\text{c}}^2 / \tau_{\text{c}}^2 = \kappa_0 / (2 \mu_0)\). For the considered special case, the strength potential reads \[\label{eq:cohesive95strength95potential} \pi_0 (\text{tr}(\boldsymbol{\eta}), \left|\boldsymbol{\eta}_{\text{dev}}\right|) = \sqrt{w_{\text{c}}} \sqrt{\kappa_0 \text{tr}(\boldsymbol{\eta})^2 + 2 \mu_0 \left| \boldsymbol{\eta}_{\text{dev}} \right|^2} \quad \text{with} \quad w_{\text{c}} = \frac{\tau_{\text{c}}^2}{2\mu_0} = \frac{p_{\text{c}}^2}{\kappa_0} \quad\text{,}\tag{43}\] hence it is \(\pi(\boldsymbol{\eta}, \alpha) = a(\alpha) \sqrt{w_{\text{c}}} \sqrt{ \mathbb{C}_0 \boldsymbol{\eta} \cdot \boldsymbol{\eta}}\), with \(\mathbb{C}_0\) as the undamaged stiffness tensor and the parameter \(w_{\text{c}}\) controlling the size of the initial elastic domain \(\mathcal{S}_0\). Substituting 42 in 34 , it is possible to optimize the action functional with respect to the eigenstrain, yielding a strain energy density depending only on the displacement and the phase field. This introduces a piece-wise defined elastic strain energy density but eliminates the need to solve for \(\boldsymbol{\eta}\) as done in [40]. Instead, the eigenstrain is evaluated at the Gauss points during assembly, analogously as with a return mapping algorithm in plasticity. Although additional nonlinear iterations are introduced whenever the eigenstrain is non-zero, the resulting problem for the displacement is unconstrained, thus improving the efficiency of the numerical computations. This approach shares some similarities with the one in [63]; however, the present case only requires the evaluation of an algebraic expression rather than the solution of a system of equations at each Gauss point.

Stationarity with respect to \(\boldsymbol{\eta}\) under the constraint \(\text{tr}(\boldsymbol{\eta})\geq 0\) yields \[\label{eq:psi95cohesive95condensed} \hat{\psi} (\boldsymbol{\varepsilon}, \alpha) = \begin{cases} a(\alpha) \sqrt{w_{\text{c}}} \sqrt{ 2\psi_0 (\boldsymbol{\varepsilon}) } - \frac{1}{2} a(\alpha)^2 w_{\text{c}} + \epsilon \psi_0 (\boldsymbol{\varepsilon}) &\text{if}\,\,\, \text{tr}(\boldsymbol{\varepsilon}) \geq 0, \psi_0 (\boldsymbol{\varepsilon}) \geq \frac{1}{2} a(\alpha)^2 w_{\text{c}}\\ \frac{\kappa_0}{2} \left( \text{tr}(\boldsymbol{\varepsilon}) \right)^2 - \frac{1}{2} a(\alpha)^2 w_{\text{c}} + a(\alpha) \sqrt{2 \mu_0 w_{\text{c}}} \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| + \epsilon \psi_0 (\boldsymbol{\varepsilon}) &\text{if}\,\,\, \text{tr}(\boldsymbol{\varepsilon}) < 0, \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| \geq a(\alpha) \sqrt{\frac{w_{\text{c}}}{2 \mu_0}}\\ (1+\epsilon) \psi_0(\boldsymbol{\varepsilon}) &\text{else}\\ \end{cases} \quad \text{.}\tag{44}\] The term \(\epsilon \psi_0\) with \(\epsilon = 10^{-7}\) is introduced, analogously to the brittle case, as a residual energy density to avoid numerical issues which would result from the linearity of \(\hat{\psi}\) in \(\boldsymbol{\varepsilon}\) in the first two branches (in the second one linearity emerges only for purely deviatoric deformations). We add this term in all branches for numerical stability. Having \(a(\alpha=1)>0\) as in the brittle models would not fix the rank-deficiency of the stiffness matrix with this compact reformulation. The detailed derivation as well as the optima for the eigenstrain and the condensed form of the stress are given in Appendix 5.4.

As in the quasi-static case, we require that the strain-hardening condition be fulfilled, which yields [40] \[\label{eq:strain95hardening} \ell \leq \ell_{\text{ch}} \qquad\text{with}\qquad \ell_{\text{ch}} = \min_{\boldsymbol{\sigma}_{\text{c}} \in \partial \mathcal{S}_0} \frac{G_{\text{c}}}{\mathbb{S}_0 \boldsymbol{\sigma}_{\text{c}} \cdot \boldsymbol{\sigma}_{\text{c}}}\tag{45}\] where \(\ell_{\text{ch}}\) is the cohesive (or Irwin) length of the material and \(\mathbb{S}_0 = \mathbb{C}_0^{-1}\) is the undamaged compliance tensor. In the present case, this condition simplifies to \(\ell \leq \ell_{\text{ch}} = G_{\text{c}} / w_{\text{c}}\).

3.3 1D cohesive model↩︎

Following the same procedure as for the brittle models, we now analyze the cohesive model in the 1D setting. The 1D elastic domain is \(\mathcal{S}_0 = \{ \sigma \in \mathbb{R} : \sigma \leq \sigma_{\text{c}} \}\), where \(\sigma_{\text{c}}>0\) denotes the tensile strength, so that the eigenstrain potential reads \[\pi (\eta, \alpha) = \begin{cases} a(\alpha) \sigma_{\text{c}} \eta &\text{if } \, \eta \geq 0\\ +\infty &\text{else} \end{cases} \quad\text{.}\] In this case, 36 and 37 simplify to \[\varepsilon = u' + \llbracket u \rrbracket \, \delta_{J(\boldsymbol{z})},\eta = \eta_R + \llbracket u \rrbracket \, \delta_{J(\boldsymbol{z})}, \label{eq:1Djump}\tag{46}\] where we see that at the location of a cohesive crack the eigenstrain coincides with the displacement jump (which we also refer to as crack opening in the following). The eigenstrain evolution criterion 41 takes the explicit form \[\label{eq:strength95criterion951D} \sigma (\varepsilon, \eta) \leq a(\alpha) \sigma_{\text{c}} \qquad \eta\geq0 \qquad \left(\sigma (\varepsilon, \eta) - a(\alpha) \sigma_{\text{c}}\right) \eta = 0 \qquad \forall (x,t) \in \Omega \times [0, T]\tag{47}\] with \(\sigma = E_0\,(\varepsilon-\eta)\). Stationarity of \(\psi\) with respect to \(\eta \geq 0\) yields the optimal eigenstrain \[\eta^\star (\varepsilon, \alpha) = \arg\,\text{stat}_{\eta \geq 0} \psi (\varepsilon, \eta, \alpha) = \begin{cases} 0 &\varepsilon < (1-\alpha) \frac{\sigma_{\text{c}}}{E_0}\\ \varepsilon - (1-\alpha) \frac{\sigma_{\text{c}}}{E_0} &\text{else}\\ \end{cases}\,.\] The resulting reduced form of the elastic strain energy density is represented in Fig. 11a and reads \[\hat{\psi} (\varepsilon, \alpha) = \psi (\varepsilon, \eta^\star(\varepsilon, \alpha), \alpha) = \begin{cases} \psi_0 (\varepsilon) + \epsilon \psi_0 (\varepsilon) &\varepsilon < (1-\alpha) \frac{\sigma_{\text{c}}}{E_0}\\ (1-\alpha) \sigma_{\text{c}} \varepsilon - \frac{1}{2} (1-\alpha)^2 \frac{\sigma_{\text{c}}^2}{E_0} + \epsilon \psi_0 (\varepsilon) &\text{else}\\ \end{cases}\qquad\text{.}\]

The KKT conditions for the phase field on the jump set read \[\label{eq:phase95field95KKT951D} a^\prime(\alpha) \sigma_{\text{c}} \llbracket u \rrbracket - \frac{2 G_{\text{c}} \ell}{c_w} \llbracket \alpha^\prime \rrbracket \geq 0 \qquad \dot{\alpha} \geq 0 \qquad \left[ a^\prime(\alpha) \sigma_{\text{c}} \llbracket u \rrbracket - \frac{2 G_{\text{c}} \ell}{c_w} \llbracket \alpha^\prime \rrbracket \right] \dot{\alpha} = 0 \qquad \forall (x,t) \in J(\boldsymbol{z})\times[0,T] \qquad\text{.}\tag{48}\] The complementarity condition in case of evolving damage (\(\dot{\alpha} > 0\)) and using the optimal AT2 profile, for which \(\llbracket \alpha^\prime \rrbracket = -2\alpha/\ell\), yields \[\label{eq:lindeg95alpha95jump} \llbracket u \rrbracket = \frac{2G_{\text{c}}}{\sigma_{\text{c}}} \alpha\quad\text{.}\tag{49}\] Combined with the complementarity condition of 47 for \(\eta > 0\), this gives the cohesive law \[\label{eq:cohesive95stress95jump} \sigma = \left(1 - \frac{\sigma_{\text{c}}}{2G_{\text{c}}} \llbracket u \rrbracket \right) \sigma_{\text{c}} \qquad\text{,}\tag{50}\] hence the model exhibits linear softening (Fig. 11b) [39]. For a displacement jump value \(\llbracket u \rrbracket_{\text{ult}}={2G_{\text{c}}}/{\sigma_{\text{c}}}\), complete failure is reached. In Fig. 11c we observe that, upon unloading from the softening branch, the displacement jump gradually closes at constant stress \(\sigma = (1-\alpha) \sigma_{\text{c}}\). Once the jump is fully closed, the response follows the initial linear elastic branch until complete unloading. The initial phase of reloading follows again the initial elastic branch up to \(\sigma = (1-\alpha) \sigma_{\text{c}}\), after which the jump reopens at constant stress until the maximum value reached previously. Then, the phase field continues to evolve with the response governed by the softening branch.

3.4 Interaction of a stress wave with a pre-existing fully developed cohesive crack in a 1D bar↩︎

We now investigate the interaction of an elastic wave with a pre-existing crack. Note that a cohesive crack can transmit cohesive forces if it is not fully developed, i.e. if the displacement jump across its faces has never reached \(\llbracket u \rrbracket_{\text{ult}}\) (and the maximum value of the phase field has never reached \(1\)). In the following, we begin with the simpler case of a fully developed crack and will treat the case of a partially developed crack in the next subsection.

Thus, at the center of our bar we place a crack that has been previously opened beyond \(\llbracket u \rrbracket_{\text{ult}}\), and whose faces have been then brought back in contact, i.e. \(\alpha=1\) at \(x=0\) and \(\eta = 0\) everywhere. From now on, we denote the value of the phase-field at the central crack as \(\breve{\alpha}(t)\) (here, \(\breve{\alpha}(0)=1\)). The initial crack is introduced by imposing the optimal AT2 profile \(\alpha(x,0) = \exp(-|x|/\ell)\); here, however, we directly allow for the phase field to evolve.

Apart from these conditions, the setup is identical to the brittle case (Section 2.1.5, Fig. 3) with a bar of length \(L=1000\) mm and material parameters \(E_0 = 30\,000\) MPa, \(G_{\text{c}} = 0.1\) N/mm, and \(\rho_0 = 2\,400\) kg/m\(^3\). As the optimal profile of the AT2 model is wider, we now set \(\ell = 20\) mm. The strength is set to \(\sigma_{\text{c}} = 5\) MPa, matching one of the values of \(\hat{\sigma}_{\text{c}}(0)\) for the brittle phase-field model 16 . We again impose a half-sine pulse at the left boundary via 24 , either compressive with \(\ell/\tilde{\lambda} = 0.05\), or tensile with \(\ell/\tilde{\lambda} = 0.05\) or \(0.15\). The pulse amplitude is again chosen such that \(\tilde{\sigma} / \sigma_{\text{c}} = 0.9\). The FEM discretization uses \(\ell / \Delta x \approx 20\) and \(c_0\Delta t / \Delta x = 0.1\).

The results are summarized in Fig. 12, where the four rows show the stress field in \(x\)-\(t\) space, the evolution of the displacement jump, the phase-field profile, and the stress along with the local strength at selected time instants. Since the stiffness is not degraded, both the wave speed and the acoustic impedance within the support of the phase field remain constant and equal to that of the pristine material, while no energy reflection takes place. The original wave equation is fully recovered, hence, the compressive wave (Fig. 12a,d,g,j) is completely transmitted across the crack, with the condition \(\eta \geq 0\) preventing interpenetration. Since the stress is always negative, the eigenstrain evolution criterion 47 remains inactive and the phase field does not evolve. This holds independently of the wave shape and frequency, as long as the wave is fully compressive. The negligible evolution of the eigenstrain up to values still below \(10^{-3}\) visible in Fig. 12d is due to numerical artifacts.

For the tensile case with \(\ell/\tilde{\lambda} = 0.05\) (Fig. 12b,e,h,k), the wave behaves as expected in the sharp-crack case. The results in Fig. 12e can be explained using 47 . In this case the stress remains always lower than the local strength of the material, i.e. \(\sigma(x,t)=E_0 u^\prime (x,t) < \left(1-\alpha(x,t)\right) \sigma_{\text{c}}\) in all points and at all times, therefore neither the eigenstrain nor the damage evolve and the behavior remains linear elastic. Also, since the fully developed cohesive crack can no longer transmit tensile cohesive stresses, it acts as a free end, hence the incoming wave is completely reflected without distortion but with reversed sign (i.e., as a compressive wave). The jump opens up exclusively at the two Gauss points of the element with \(\alpha = 1\), while neighboring points with \(\alpha < 1\) remain unaffected.

For the tensile wave with the shorter wavelength \(\ell/\tilde{\lambda} = 0.15\) (Fig. 12c,f,i,l) the behavior is different. Here, the stress gradient of the incoming wave is steeper than the spatial variation of the material strength \(a(\alpha(x,t))\sigma_{\text{c}}\) (reported as a dashed line in Fig. 12l), leading to the evolution of the eigenstrain and, hence, of the damage at several locations in the neighborhood of the crack center. Although much less severely than in Fig. 7, this triggers the widening of the damage band and an increase of the dissipated energy. Each of the jumps reflects and transmits a portion of the wave, and the subsequent closure of the jumps introduces stress discontinuities in time or shock waves (better illustrated in Section 3.5), causing the observed oscillations.

The above results suggest that the interaction of stress waves with a fully developed cohesive crack depends on the \(\ell/\tilde{\lambda}\) ratio, with smaller ratios leading to the expected sharp-crack response and larger ratios yielding diffuse jumps accompanied by shock waves and high-frequency oscillations. To achieve a better quantitative understanding of this phenomenon, in the following we formulate a general condition for the emergence of diffuse jumps.

Preventing the occurrence of diffuse jumps requires that the local strength of the material is reached only at the center of the crack, i.e. where \(\alpha=1\) and the strength vanishes. At any other position, the stress must be strictly below the strength profile at all times, namely \[\label{eq:stress95strengthprofile} \sigma (x, t) < a(\alpha (x,t)) \sigma_{\text{c}} \qquad \forall x \in \Omega \setminus \{0\}\quad\text{,}\quad\forall t \in [0,T]\text{.}\tag{51}\] To check whether this condition is fulfilled, we trace the maximum stress reached during the reflection process at each material point and compare it to the strength profile. For an incoming wave 24 , the temporal maximum at each spatial point reads \[\label{eq:max95stress} \max_{t} \sigma (x, t) = \begin{cases} \tilde{\sigma} \sin (4 \pi \tfrac{|x|}{\tilde{\lambda}}) &\text{for} \, |x| \leq \frac{\tilde{\lambda}}{8}\\ \tilde{\sigma} &\text{else}\\ \end{cases}\quad\text{with} \quad \tilde{\sigma}\le\sigma_c\quad\text{.}\tag{52}\] The strength profile is given by the optimal AT2 profile \[\label{eq:opt95profile95AT2} a(\alpha(x,t))\sigma_{\text{c}} = \left(1 - \alpha(x,t)\right) \sigma_{\text{c}} \qquad\text{with}\qquad \alpha (x,t) = \breve{\alpha}(t) \exp \left( - \frac{|x|}{\ell} \right) \qquad\text{.}\tag{53}\] Both 52 and 53 are illustrated in Fig. 13a for the case studied here, i.e. with \(\breve{\alpha}(0)=1\) and \(\ell/\tilde{\lambda} = 0.05\) or \(0.15\). Using 52 and 53 in 51 gives the condition \[\label{eq:nodiffusecrack95cond} \frac{\tilde{\sigma}}{\sigma_{\text{c}}} \sin \left(4 \pi \frac{\ell}{\tilde{\lambda}} \frac{|x|}{\ell}\right) < \left(1 - \breve{\alpha}(0) \exp \left(-\frac{|x|}{\ell}\right)\right) \quad \tilde{\sigma}\le\sigma_c\qquad \forall x \in [-L/2, L/2] \setminus \{0\}\quad\text{,}\tag{54}\] which depends on the \(\tilde{\sigma}/\sigma_{\text{c}}\) and \({\ell}/\tilde{\lambda}\) ratios.

Evaluating 54 with an equal sign gives the limit which separates the cases where no diffuse jumps are expected from those where the boundary of the elastic domain is reached, triggering the evolution of eigenstrain and damage. The resulting domains for the case at hand are illustrated in Fig. 13b. The red line represents the limit case; for \(\tilde{\sigma}/\sigma_{\text{c}}\) and \({\ell}/\tilde{\lambda}\) ratios falling within the lower white area, we expect the correct wave-crack interaction, i.e. a sharp-crack behavior, whereas \(\tilde{\sigma}/\sigma_{\text{c}}\) and \({\ell}/\tilde{\lambda}\) ratios within the upper gray-shaded area are expected to lead to diffuse jumps accompanied by shock waves and high-frequency oscillations. The tensile cases illustrated in Fig. 12 are marked and, as expected, the waves with \(\ell/\tilde{\lambda} = 0.05\) and \(0.15\) lie respectively below and above the red limit line, further confirming the obtained results. Although the criterion 54 is derived for harmonic waves, we can draw some general observations. For sufficiently small stress wave amplitudes compared to the initial material strength (\(\tilde{\sigma}/\sigma_{\text{c}} \rightarrow 0\)), the whole range of \(\ell/\tilde{\lambda}\) ratios recovers the expected sharp-crack behavior. The same occurs for the whole range of amplitudes \(\tilde{\sigma}/\sigma_{\text{c}}\) if the regularization length is sufficiently small compared to the wavelength (\(\ell/\tilde{\lambda} \rightarrow 0\)). In this context, an advantage of the cohesive model lies in the interpretation of the regularization length \(\ell\) as a pure numerical parameter, which can be adapted to the expected wavelength. We remark that for the special case of the applied wave, these results do not only hold for the half-wave of width \(\tilde{\lambda}\), but also for the full wave of wavelength \(\lambda\). Note that additionally the strain-hardening condition 45 must be satisfied, independently of \(\ell/\lambda\). However, in case of high-frequency waves (\(\lambda\to0\)) or shock waves, avoiding diffuse jumps calls for \(\ell \to 0\).

Note that, although the dependence on \(\ell/\lambda\) resembles that of the brittle model with stiffness degradation of Section 2.1.5 (both recovering the sharp crack for small ratios), the cohesive model admits two additional degrees of freedom to reach it: the sharp-crack behavior is restored for any \(\ell/\lambda\) at small enough \(\tilde{\sigma}/\sigma_{\text{c}}\), and \(\ell\) – being decoupled from the strength unlike in the brittle case – can be freely reduced to accommodate high-frequency waves.

3.5 Interaction of a stress wave with a pre-existing partially developed cohesive crack in a 1D bar↩︎

We now place at the center of our bar a cohesive crack that has experienced a maximum damage \(\breve{\alpha} = 0.5\), hence a maximum displacement jump \(2 G_{\text{c}} \breve{\alpha} / \sigma_{\text{c}}\) (see 49 ). This crack can still transmit a cohesive stress \(a(\breve{\alpha})\sigma_{\text{c}}\) across its faces, and we further assume that its faces have been then brought back in contact, i.e. \(\eta = 0\) everywhere. We adopt the same setup as before, but we set \(\ell = 10\) mm. Moreover, we choose \(\ell/\Delta x = 20\), \(c_0 \Delta t / \Delta x = 0.1\), \(\ell / \tilde{\lambda} = 0.015\) and \(\tilde{\sigma} / \sigma_{\text{c}} = 0.9\). To test different cohesive responses, we vary the fracture toughness as \(G_{\text{c}} \in \{0.01, 0.025, 0.1\}\) N/mm, yielding different Irwin lengths \(\ell_{\text{ch}} = E_0 G_{\text{c}} / \sigma_{\text{c}}^2\) following 45 .

The numerical results are shown in Fig. 14, while the time evolution of the displacement jump \(\llbracket u \rrbracket(t)\), the phase-field value at the crack \(\breve{\alpha}(t)\), and the cohesive stress \(\sigma(0,t)\) at the center of the crack are illustrated in Fig. 15. Once the wave reaches the crack, the response shows four different phases:

1. As long as the strength criterion at the crack is not met, i.e. \(\sigma(0,t) < a(\breve{\alpha}(0))\sigma_{\text{c}}\), the wave is fully transmitted and \(\eta = 0\), \(\dot{\eta} = 0\), \(\dot{\alpha} = 0\).

2. When at time \(t_{I}\) the cohesive strength is reached, i.e. \(\sigma(0,t_{I}) = a(\breve{\alpha}(0))\sigma_{\text{c}}\), the crack starts reopening, \(\dot{\eta} > 0\), with a transmitted stress equal to \(a(\breve{\alpha}(0))\sigma_{\text{c}}\) while the damage remains constant and equal to \(\breve{\alpha}(0)\). This phase continues until the displacement jump reaches the maximum value attained in the past, namely \(2 G_{\text{c}} \breve{\alpha}(0) / \sigma_{\text{c}}\), at time \(t_{II}\). This corresponds to the horizontal branch of the reloading path in Fig. 11c.

3. Once this value of the opening is reached, the damage evolution criterion is fulfilled and the phase field starts to evolve along with the eigenstrain, i.e. \(\dot{\eta} > 0\) and \(\dot{\alpha} > 0\) (softening branch in Fig. 11c). Provided that the incoming wave supplies sufficient energy, this phase continues until at the crack center \(\alpha = 1\) and accordingly \(\llbracket u \rrbracket_{\text{ult}}=2 G_{\text{c}} / \sigma_{\text{c}}\) is reached at time \(t_{III}\).

4. From this point on, no cohesive stresses can be transmitted. Therefore, the crack behaves as a free end and any further portion of wave is reflected while the displacement jump keeps increasing.

As Fig. 14i shows, for the case with \(G_{\text{c}} = 0.1\) N/mm the energy carried by the wave is not enough to reach \(\llbracket u \rrbracket=2 G_{\text{c}} \breve{\alpha}(0) / \sigma_{\text{c}}\) needed to activate the damage evolution criterion. Therefore, after the initial linear elastic reloading phase, the transmitted stress remains constant and the phase field does not evolve. The displacement jump first increases and, when the energy carried by the incoming wave starts decreasing, it begins decreasing again following the horizontal unloading branch in Fig. 11c. Once the stress of the incoming wave drops below the cohesive strength, the jump closes completely, triggering a shock wave propagating through the domain (Fig. 14l). This behavior can be explained considering that the constant stress during the jump closure forces also the velocities of the crack faces to be constant due to the balance of linear momentum. When the jump finally closes, the faces velocities drop abruptly from a constant non-zero value to zero, producing a discontinuity in velocity, hence, a shock wave. As discussed in Section 3.4, this creates high-frequency oscillations as visible in Fig. 12l and also, albeit to a more limited extent, in Fig. 14l. This explains the scatter visible in Fig. 12c: since the strength criterion is reached at multiple points, multiple cracks open and close again subsequently, which triggers multiple smaller shock waves. With the emergence of shock waves, the assumption of ‘sufficiently regular’ temporal derivatives made in the model derivation (Section 3.1) is no longer satisfied. The behavior, directly stemming from the reversibility of the eigenstrain, might be reduced or avoided by modifying the reversibility condition so as to involve a different unloading/reloading behavior as also noted in [40] and [64]. This is, however, outside of the scope of the present work.

From Fig. 15 we note that a smaller value of \(G_{\text{c}}\) leads to a faster rate of change of the displacement jump. To make this dependence quantitative, we analytically derive the evolution in time of the cohesive crack opening from the kinematic conditions, stress continuity, and the complementarity conditions for the eigenstrain and the phase field. The full derivation is given in Appendix 5.5, while here we report only the essential results. Under the condition that the strength criterion is met at a single point only (i.e., for sufficiently small \(\ell/\lambda\)), the opening evolution is governed by the ordinary differential equation (ODE) \[\label{eq:opening95evolution95ODE95mainbody} \dot{v}_L^{\text{ref}}(0^-,t) - \frac{c_0}{\ell_{\text{ch}}} v_L^{\text{ref}}(0^-,t) = \dot{v}_L^{\text{inc}}(0^-,t) \qquad \forall t_{\text{II}} \leq t \leq t_{\text{III}} \qquad\text{,}\tag{55}\] where \(v_L^{\text{ref}}\), \(v_L^{\text{inc}}\) are respectively the reflected and the incoming velocities of the wave traveling toward the left at the crack position. The ODE admits the closed-form solution \[v_L^{\text{ref}}(0^-,t) = \exp \left( \tfrac{c_0}{\ell_{\text{ch}}} (t-t_{\text{II}}) \right) v_L^{\text{ref}}(0^-,t_{\text{II}}) + \int_{t_{\text{II}}}^{t} \exp \left( \tfrac{c_0}{\ell_{\text{ch}}} (t-\tau) \right) \dot{v}_L^{\text{inc}}(0^-,\tau) \mathrm{d}\tau \qquad \forall t_{\text{II}} \leq t \leq t_{\text{III}} \quad\text{,}\] from which the crack opening follows as \[\label{eq:dynamic95opening95evolution} \llbracket u \rrbracket (t) = - 2 \tfrac{\ell_{\text{ch}}}{c_0} \Bigl[ v_L^{\text{ref}}(0^-,t) - v_L^{\text{ref}}(0^-,t_{\text{II}}) - \left(v_L^{\text{inc}}(0^-,t) - v_L^{\text{inc}}(0^-,t_{\text{II}})\right) \Bigr] - 2 u_L^{\text{ref}}(0^-,t_{\text{II}}) \quad\text{.}\tag{56}\] From 56 we deduce that the evolution is independent of \(\ell\) and governed only by the \(c_0/\ell_{\text{ch}}\) ratio and by the incoming wave, which determines the initial condition at \(t_{\text{II}}\). An interactive visualization of the theoretical solution is available at https://github.com/jonas-heinzmann/phase_field_dynamics in the folder /math/crackopening.py. In Appendix 5.5.2, we also evaluate these expressions for the specific case of the half-sine pulse used in the numerical computations. There, we also show that the theoretical results and the FEM results from Fig. 15a coincide.

Analyzing 55 the following two limit cases can be obtained:

• For \(c_0/\ell_{\text{ch}} \rightarrow 0\) (small wave speed or very large Irwin length), the second term on the left-hand side vanishes and \(\dot{v}_L^{\text{ref}}(0^-,t) = \dot{v}_L^{\text{inc}}(0^-,t)\), the reflected and incoming wave velocities can differ only by a constant offset, while \(\breve{\alpha}(t)\) does not evolve and a constant cohesive stress \(a(\breve{\alpha}(t))\sigma_{\text{c}}\) is transmitted.

• For \(c_0/\ell_{\text{ch}} \rightarrow \infty\), the ODE becomes singularly perturbed and the phase field jumps instantaneously to \(\alpha = 1\) since the time range \(t_{III}-t_{II}\to0\). This corresponds to the brittle limit \(G_{\text{c}}/\sigma_{\text{c}} \rightarrow 0\) (or \(\ell_{\text{ch}} \rightarrow 0\)), for which no cohesive softening is present.

Fig. 16, which illustrates the time evolution of \(\llbracket u \rrbracket(t)\) and \(\sigma/\sigma_c\) for different \(c_0/\ell_{\text{ch}}\) ratios for the harmonic wave used here, confirms these observations. The derivation of the curves in Fig. 16 is reported in Appendix 5.5.2.

3.6 Summary of the 1D analyses: Stiffness and stiffness+density degradation vs. strength degradation↩︎

We close the 1D analyses by comparing in Fig. 17 how the phase-field regularization affects the local material properties in the three formulations and briefly summarizing the results. The figure uses AT1 for the brittle models and AT2 for the cohesive model, consistently with the dissipation functions adopted in Sections 2 and 3; the qualitative comparison is independent of this choice. In the brittle model with stiffness degradation (Fig. 17a), the wave speed and the acoustic impedance are reduced within the phase-field support, and elastic waves are partially reflected and partially transmitted at the regularized crack. Combined with the locally reduced critical stress \(\hat{\sigma}_{\text{c}}(\alpha)\), this drives the widening of the phase-field profile (Fig. 7). Compressive waves are transmitted correctly. Degrading the density alongside the stiffness restores a constant wave speed and leads to an even more strongly reduced acoustic impedance for tensile loading, which yields full reflection (Fig. 17b). However, mass balance is not fulfilled; moreover, in compression the stiffness remains undegraded while the density is degraded, leading to partial reflection and transmission. Both brittle variants recover the sharp-crack response – including undistorted waveforms and correct reflection or transmission – only in the limit \(\ell/\lambda \rightarrow 0\) (Section 2.1.5).

In the cohesive model (Fig. 17c), neither the stiffness nor the density is degraded, so both the bulk wave speed and the acoustic impedance are preserved while fulfilling mass balance. The phase-field regularization reduces the local strength in the vicinity of the crack, while the sharp crack itself is represented by the eigenstrain at the jump set. No strain energy decomposition is needed since tension-compression asymmetry follows from the eigenstrain constraint \(\text{tr}(\boldsymbol{\eta}) \geq 0\). The only restriction is that elastic waves must reach the strength criterion exclusively on the jump set and not in its vicinity, leading to limits in terms of \(\ell/\lambda\) and \(\tilde{\sigma}/\sigma_c\) ratios for which the sharp-crack behavior is fully recovered (Section 3.4).

4 Comparison of the models in the multi-dimensional setting↩︎

We now compare the three formulations on a benchmark problem in 2D plane-strain conditions (\(d=3\)). Preliminarily, we compare the formulations of the elastic domain for the three models and introduce a closeness measure to the strength surface.

4.1 Elastic domains and a closeness measure to the strength surface for the brittle and cohesive models↩︎

The equations of motion in 3D for the brittle models with stiffness and stiffness+density degradation and for the cohesive model are compared in Table 1, while the expressions for the stress and energy release rate are presented in Table 2. We remark that, in contrast to pure linear elastodynamics, the dilational and deviatoric waves are no longer decoupled as soon as the phase field is non-zero for the brittle models in the multidimensional setting. This comes from the degradation of the elastic properties, and is expected for heterogeneous materials [65].

We analyze now the elastic domains obtained for the brittle and cohesive local models (i.e., for \(\nabla\alpha \equiv \boldsymbol{0}\)). For the brittle models the elastic domain is a consequence of the adopted strain energy decomposition [40] and of the selected fracture toughness, regularization length and elastic parameters [46], [53]. It is obtained from the KKT conditions 13 ; for the AT1 dissipation function, quadratic degradation function, and volumetric-deviatoric split 3 it reads [46] \[\label{eq:elastic95domain95brittle} \frac{\langle \text{tr}(\boldsymbol{\sigma}) \rangle_+^2}{d^2 \kappa_0} + \frac{\left|\boldsymbol{\sigma}_{\text{dev}}\right|^2}{2 \mu_0} \leq \begin{cases} \frac{3 (1-\alpha)^3 G_{\text{c}}}{8\ell} &\text{stiffness degradation}\\ \frac{3 (1-\alpha)^3 G_{\text{c}}}{8\ell} - h^\prime (\alpha) (1-\alpha)^3 \frac{\rho_0}{2}\left|\dot{\boldsymbol{u}}\right|^2 &\text{stiffness+density degradation}\\ \end{cases}\tag{57}\] where the residual stiffness is neglected and \(d\) is the spatial dimension. Although the elastic domain is affected by the damage gradient, looking at the elastic domains for the local models is useful to compare their behavior. Moreover, the initial elastic domain obtained for zero damage, i.e. for an intact domain, is exact and its boundary provides the elastic limit for the pristine material. Evaluating 57 along the pure hydrostatic pressure (\(\tau=0\)) and pure shear (\(p=0\)) directions yields the damage-dependent critical pressure \(\hat{p}_{\text{c}}(\alpha)\) and shear \(\hat{\tau}_{\text{c}}(\alpha)\) given in Table 3. As observed in Section 2.2 for the 1D case, also in 3D the model with stiffness+density degradation leads to a strength surface that depends on the local velocity \(\dot{\boldsymbol{u}}\). For both models the initial elastic domain is an ellipse for positive pressure and an unbounded cylinder in the negative pressure direction and it shrinks homothetically with damage, so the elastic domain \(\mathcal{S}(\alpha)\) can be written as \[\label{eq:elastic95domain} \begin{cases} \tau \leq \hat{\tau}_{\text{c}} (\alpha) &p < 0\\ \left( \frac{p}{\hat{p}_{\text{c}} (\alpha)} \right)^2 + \left( \frac{\tau}{\hat{\tau}_{\text{c}} (\alpha)} \right)^2 \leq 1 &p\geq0\\ \end{cases} \quad\text{.}\tag{58}\]

For the cohesive model, the elastic domain is governed by the strength potential \(\pi_0(\boldsymbol{\eta})\) and by the degradation function \(a(\alpha)\), which homothetically scales \(\pi_0(\boldsymbol{\eta})\). In particular, the strength potential is independent of \(\ell\) and \(\nabla\alpha\). Also, it is parametrized through \(p_{\text{c}}\) and \(\tau_{\text{c}}\), which are independent input parameters. Although its choice is flexible [40], with the adopted strength potential 43 the elastic domain can still be written as 58 , with the damaged critical pressure and shear given in Table 3 for \(a(\alpha)=1-\alpha\).

From Fig. 18, where the strength surfaces of the brittle model with stiffness degradation and of the cohesive model are compared, we can observe that the initial elastic domain is the same for both models. However, with our choice of degradation functions, the brittle elastic domain shrinks at a faster rate compared to the cohesive one. The brittle model with stiffness+density degradation is not reported since it leads to an elastic domain that depends on the velocity, whereas it gives the same elastic domain of the brittle model with stiffness degradation in the quasi-static case. Exploiting 58 and Table 3, we finally define a closeness measure \(s \in [0,1]\) of a stress state to the damaged strength surface as follows: \[s(p, \tau, \alpha) = \begin{cases} \frac{\tau}{\hat{\tau}_{\text{c}} (\alpha)} &p<0\\ \sqrt{\left( \frac{p}{\hat{p}_{\text{c}} (\alpha)} \right)^2 + \left( \frac{\tau}{\hat{\tau}_{\text{c}} (\alpha)} \right)^2} &p\geq0\\ \end{cases}\] illustrated with green arrows in Fig. 18. The measure is devised so that when \(s=1\) a stress state lies at the boundary of the elastic domain. We recall that, for the cohesive model, \(s\) is exact by virtue of the explicit strength criterion, while for the brittle models it is exact only under homogeneous damage conditions.

4.2 Comparison in the brittle limit↩︎

To compare the various models, let us consider the benchmark problem of a brittle notched plate under tension, adopting the setup from [22] and [66]. The setup is presented in Fig. 19a and involves a rectangular 2D domain of length \(L=100\) mm and height \(H=40\) mm with an initial crack of length \(a_0 = 50\) mm starting from the left edge and placed at \(y=H/2\). We describe this initial crack through an initial phase field, enabling us to evaluate the interaction of the waves with the initial phase-field crack. Following [22], we consider plane-strain conditions, Young’s modulus \(E_0 = 32\,000\) MPa, Poisson’s ratio \(\nu_0 = 0.2\), mass density \(\rho_0 = 2\,450\) kg/m\(^3\), fracture toughness \(G_{\text{c}} = 0.003\) N/mm, and regularization length \(\ell = 0.25\) mm. According to Table 3, the resulting initial critical pressure and shear for the brittle model with stiffness degradation are \(p_{\text{c}} = 8.94\) MPa and \(\tau_{\text{c}} = 10.95\) MPa. The same values are also adopted for the cohesive model to allow for a fair comparison. The maximum of the ultimate displacement jumps in case of pure hydrostatic or shear loading, i.e. \(2G_{\text{c}} / \min(p_{\text{c}}, \tau_{\text{c}}) = 0.6\;\mu\text{m}\), is five orders of magnitude smaller than the smallest structural dimension, which indicates that we are close to the brittle limit. Following Section 3.5, we expect the ratios of the wave speeds to the cohesive length to govern the dynamic crack opening evolution also in the multi-dimensional case. Given the value of \(3.5 \cdot 10^6~\text{s}^{-1}\) obtained here, we expect that the period during which the crack transmits cohesive forces (phase III in Section 3.5) is close to zero.

Unlike in [22] where the tensile loading is applied abruptly, we impose a traction on the top and bottom edges via a smooth ramp with \(T_{\text{pulse}}=1.6\) \(\mu\)s (Fig. 19b) to avoid introducing a shock wave. We consider a traction per unit thickness of final magnitude \(\hat{f}_y \in \{1, 2\}\) N/mm to promote more or less dominant crack branching. The time domain is \(t \in [0, T]\) with \(T=80\) \(\mu\)s and \(\Delta t = 0.01\) \(\mu\)s. The spatial discretization employs a structured mesh of linear quadrilateral elements with \(\ell / \Delta x = \ell / \Delta y \approx 5\), yielding \(1\,604\,802\) nodes and \(1\,602\,000\) elements. This gives \(c_{\text{S}_0} \Delta t / \Delta x = 0.47\), with \(c_{\text{S}_0}\) denoting the shear wave speed (lower than the compressive wave speed).

4.2.1 Single branching case↩︎

We start with \(\hat{f}_y = 1\) N/mm, for which we observe a single crack branching event. Fig. 20 shows the phase field before first branching at \(t=30\;\mu\text{s}\), shortly after the branching at \(t=40\;\mu\text{s}\), with fully developed branches at \(t=60\;\mu\text{s}\), and at the final time \(t=80\;\mu\text{s}\) for the three models.

While the brittle model with stiffness degradation and the cohesive model both branch, the brittle model with stiffness and density degradation produces straight crack propagation up to complete separation of the domain in two parts. In [38] the crack tip splitting for the brittle case is attributed to the wave speed reduction within the damaged regions, which causes an accumulation of energy at the crack tip. With this interpretation, the different behavior of the brittle model with stiffness+density degradation is likely due to the constant wave speed under tensile loading which avoids any energy accumulation at the crack tip. Consistently with the closeness to the brittle limit and with the large wave speeds to cohesive length ratios of the adopted cohesive model, in the cohesive results the phase field directly evolves to \(\alpha=1\) during propagation. The plots of \(s(\boldsymbol{x})\) further reveal that for the brittle model with stiffness degradation and for the cohesive model the crack tip emits ripples into the domain alongside the elastic energy release. For the brittle model with stiffness+density degradation no ripples are observed, likely due to the additional velocity-dependent term which decreases the value of the energy release rate. The obtained results are further confirmed in Fig. 21, where vertical slices of \(s\) and \(\alpha\) are shown. High-frequency oscillations of \(s\) related to the ripples are visible for the brittle model with stiffness degradation and for the cohesive model but not for the brittle model with stiffness+density degradation. For the brittle model with stiffness degradation, the oscillations stem from the interaction of the elastic waves with the regularized crack leading to partial reflection and transmission and from the widening of the phase-field support discussed in Section 2.1.5. Instead, for the cohesive model they correspond to the occurrence of diffuse jumps due to the strength criterion being met in the vicinity of the crack (Section 3.4).

The brittle model with stiffness degradation and the cohesive model produce very similar crack patterns, which we attribute to their identical initial elastic domain. Also, the different rates of homothetic shrinkage with damage evolution (Fig. 18) do not play a big role since \(\alpha\) rapidly reaches \(1\). The cohesive model has a slightly higher tendency to branch, especially as the cracks approach the boundaries. We ascribe this to the brittle model dissipating additional energy by spuriously widening the phase-field profile – as a result, as discussed earlier, the energy per unit crack extension increases over \(G_{\text{c}}\) – whereas in the cohesive model this excess energy is available for branching.

Further information comes from the elastic and fracture energy contributions (Fig. 22a,b). The brittle model with stiffness degradation and the cohesive model show very similar fracture energy evolutions, mirroring their similar crack patterns. Their elastic energies also agree up to roughly \(t=40\) \(\mu\)s, after which both brittle models show a higher elastic energy storage (Fig. 22b). Such difference can be traced back to the distinct representations of a fully developed crack. In the brittle models, the crack opening is smeared as an elastic strain over the localization band; since the degradation is complete only on the centerline of the phase-field profile, the elastic energy density \(g(\alpha) \psi_{\text{D}}(\boldsymbol{\varepsilon})\) does not vanish in the band, where small values of \(g(\alpha)\) meet very large strains. As the crack faces continue to separate under the sustained loading, the strains in the band — and with them this stored energy — keep growing, which explains the increasing surplus of elastic energy after \(t\approx40\) \(\mu\)s. In the cohesive model, by contrast, the opening is accommodated by the eigenstrain \(\boldsymbol{\eta}\), and beyond the ultimate displacement jump the stored energy is insensitive to further opening; the elastic energy therefore reflects only the bulk response and saturates accordingly. The system mass \(m=\int_\Omega h(\alpha) \rho_0\, \mathrm{d}\boldsymbol{x}\), normalized by the undegraded mass \(m_0 = \int_\Omega \rho_0 \mathrm{d}\boldsymbol{x}\), is shown in Fig. 22c. The brittle model with stiffness+density degradation exhibits a mass loss of roughly \(1.2\) % during crack propagation.

The crack tip speeds \(\dot{a}\) (Fig. 22d) are computed assuming the presence of, at most, two crack tips as the moving averages of the temporal derivative of the crack tip positions. Namely, we track the \(x\)-coordinates corresponding to minimum and maximum \(y\)-coordinate where \(\alpha \geq 0.95\). When no branching occurs, both of them coincide, while in the event of branching or widening, their vertical distance increases. In particular, we identify the start of widening as the instant at which upper and lower tip coordinates are separated by more than \(\Delta y\) for the first time. The instant at which the interpolation of the phase-field variable on the line connecting the two tips shows values below \(0.95\) for the first time marks the end of the widening, and hence the occurrence of branching.

As expected for stable dynamic crack propagation [2], [67], the brittle model with stiffness degradation and the cohesive model both remain below \(60\) % of the Rayleigh wave speed \(c_{\text{R}_0}\) before branching, after which the individual branches slow down and then re-accelerate. The cohesive model attains slightly slower speeds and yields an additional branching event around \(t=65\) \(\mu\)s accompanied by spikes in the measured crack tip speeds. In the final part of the test, the crack tip velocity decreases for all models, due to the reflection of the tensile wave at the boundary that introduces compressive stresses. For \(t \lesssim 25\) \(\mu\)s, the brittle model with stiffness+density degradation shows a similar behavior to the other two; at later times, the crack tip speed reaches values above \(60\) % of \(c_{\text{R}_0}\) without branching. It then decreases around \(t \approx 35\) \(\mu\)s due to the boundary wave reflections.

In Fig. 22d we have also marked with dashed lines the time instants at which widening of the phase-field profile starts and ends. Accordingly, the second dashed line of a respective color marks the branching event. Compared to the cohesive model, the brittle model with stiffness degradation exhibits a longer region of widening prior to branching, as well as a more extensive widening resulting in a thicker damage band (Fig. 23). This behavior is consistent with the observations made in Section 3.4 for the 1D model. The separated branches themselves do not show a significantly widened phase-field profile in either model.

4.2.2 Multiple branching case↩︎

Fig. 24 shows the phase field and \(s(\boldsymbol{x})\) for \(\hat{f}_y = 2\) N/mm. The brittle model with stiffness+density degradation again does not branch, producing a single straight crack up to the right edge. Conversely, the other two models branch significantly more compared to the case in Section 4.2.1, with the cohesive model showing a higher number of tip splitting events. This is consistent with the observation in Fig. 23. The brittle model with stiffness degradation triggers a longer widening period of the phase-field profile, which spuriously dissipates more energy; this energy is no longer available for further tip splittings. In addition, in this case the crack branches generated after a crack tip splitting event display a wider phase-field profile (Fig. 24). The cohesive model both branches more and yields a less symmetric pattern, features we attribute to its higher sensitivity to small perturbations.

To examine the branching mechanism of the cohesive model in more detail, Fig. 25 shows a close-up of \(\alpha(\boldsymbol{x})\), \(s(\boldsymbol{x})\), and the volumetric and deviatoric parts of the eigenstrain, \(\text{tr}(\boldsymbol{\eta})\) and \(\left|\boldsymbol{\eta}_{\text{dev}}\right|\) at different instants. At a branching event, the eigenstrain is not localized in a single element but is instead diffuse over a region close to the crack: the smooth transition from a single crack to two branches necessarily widens the damaged band, and with it the eigenstrain. A diffuse jump band persists, which we attribute to the interactions between the ripples emitted by the crack tip, as well as the interaction with the boundaries. The resulting opening and closing of these diffuse jumps then generate shock waves and with them high-frequency oscillations, which in turn foster the further widening of the crack band.

5 Conclusions↩︎

In this work, we analyzed the influence of phase-field regularization on the dynamic extensions of both brittle and cohesive phase-field fracture formulations. The following main insights were obtained.

For the brittle models, the interaction of a phase-field crack with an incoming elastic wave is dictated by the profiles of the wave speed and of the acoustic impedance, which depend on the choice of the degradation functions for the stiffness and for the mass density. For the model with stiffness degradation, we demonstrated the occurrence of spurious partial reflection and transmission of tensile elastic waves. The overall wave-crack interaction is controlled by the ratio \(\ell/\lambda\) of the phase-field regularization length to the elastic wavelength, and the sharp-crack response is recovered fully only in the asymptotic limit where \(\ell/\lambda \rightarrow 0\). Incidentally, this limit is not compatible with the interpretation of \(\ell\) as a material parameter, which is needed with this model to calibrate the desired nucleation behavior under quasi-static conditions. Furthermore, wave-crack interactions directly trigger a widening of the phase-field profile, causing the energy dissipated per unit crack extension to spuriously increase over the material’s critical energy release rate \(G_{\text{c}}\). However, with this model a phase-field crack correctly transmits a compressive wave if combined with a standard energy decomposition. The brittle model with stiffness and density degradation, while leading to a constant wave speed and to full reflection of a tensile wave if both degradation functions coincide, inherently violates mass conservation and cannot correctly transmit a compressive wave with a standard strain energy decomposition.

To attempt to resolve the intrinsic issues of brittle regularization in the dynamic setting, we extended our recently proposed cohesive phase-field fracture model with strength degradation to elastodynamics. Since it does not modify the stiffness nor the density as a result of damage, this model leaves the bulk wave equation unchanged with respect to the one of the pure linearly elastic problem. Moreover, with this model a phase-field crack correctly transmits a compressive wave with no need for an energy decomposition; it also correctly (i.e. completely) reflects a tensile wave under the condition that this wave does not exceed the locally degraded strength adjacent to the crack, which is controlled by the combination of the \(\ell/\lambda\) ratio and of the ratio \(\tilde{\sigma}/\sigma_c\) between the amplitude of the stress wave and the tensile strength of the material. Incidentally, with this model the value of \(\ell\) is completely independent from the material strength and can be chosen based on purely numerical considerations. For the one-dimensional case, we further derived an analytical crack-opening evolution law and demonstrated that, for a combination of \(\ell/\lambda\) and \(\tilde{\sigma}/\sigma_c\) leading to the full reflection of a tensile wave, the dynamic response in time is fundamentally governed by the ratio \(c_0/\ell_{\text{ch}}\) of the bulk wave speed to the Irwin cohesive characteristic length.

Finally, a comparative study on a 2D pre-notched plate under tensile loading highlighted critical differences in the predictions of the three models. While the brittle model with stiffness degradation and the cohesive model yielded remarkably similar crack branching patterns, the brittle model with stiffness and density degradation failed to capture crack-tip branching instabilities, resulting in purely straight crack propagation.

Taken together, these results identify the cohesive formulation as a promising starting point for further work on dynamic phase-field fracture.

Acknowledgements↩︎

We gratefully acknowledge funding from the Swiss National Science Foundation through Grant No. 200021-219407 ‘Phase-field modeling of fracture and fatigue: from rigorous theory to fast predictive simulations’.

Data and code availability↩︎

Data will be made available upon request. The implementations used in this work are publicly available at https://github.com/jonas-heinzmann/phase_field_dynamics. Animations of all numerical results can be found at https://github.com/jonas-heinzmann/phase_field_dynamics in the folder animations.

Declaration of generative AI in the writing process↩︎

During the preparation of this work, the authors used Anthropic Claude in order to improve the readability of the manuscript. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

Appendices↩︎

5.1 Scattering analysis of the 1D heterogeneous wave equation↩︎

In this section, we analyze the frequency-dependent reflection and transmission of an elastic wave in a 1D medium with smoothly varying Young’s modulus \(E(x) = g(x)E_0\) and density \(\rho(x) = h(x)\rho_0\), as obtained from both brittle models for a fixed phase field (whereby the model with stiffness degradation takes \(h(x)\equiv 1\)). We first derive an exact reformulation of the wave equation that separates propagation from reflection [50], from which we can obtain the adiabatic criterion 21 ; the transfer matrix method (TMM) then arises as its exact solution for layerwise-constant properties.

Non-dimensionalizing the elastic wave equation by \(\check{x}=x/\ell\) and \(\check{t}=c_0 t/\ell\), as well as \(\check{u} = u / u_0\) with some reference displacement magnitude \(u_0\), yields \[\label{wave} \frac{\partial}{\partial \check{x}}\!\left(g(\check{x})\,\frac{\partial \check{u}}{\partial \check{x}}\right) = h(\check{x})\,\frac{\partial^{2}\check{u}}{\partial \check{t}^{\,2}} \quad\text{,}\tag{59}\] where, with a little abuse of notation, we are denoting \(g(\check{x})=g(\alpha(\check{x}))\) and \(h(\check{x})=h(\alpha(\check{x}))\) as the compositions of the degradation functions for stiffness and density, respectively, with the phase-field profile 18 . Since, for a fixed phase field, the properties of the medium are time-independent, we work in the frequency domain: for a time-harmonic field \(\check u (\check x, \check t)=\text{Re}[U(\check x)e^{-i\check \omega \check t}]\) at dimensionless angular frequency \(\check\omega = \omega\ell/c_0\), and denoting with uppercase letters the complex amplitudes of the corresponding lower-case fields, we introduce the velocity and stress amplitudes \(V = -i\check\omega\,U\) and \(\Sigma = g\,U'\), by which 59 can be reformulated as the first-order system \[V' = -\frac{i\check k}{\check Z}\,\Sigma , \qquad \Sigma' = -i\check k\,\check Z\,V , \qquad\text{with}\qquad \check k(\check x) = \check\omega\sqrt{\frac{h(\check x)}{g(\check x)}} , \quad \check Z(\check x) = \sqrt{g(\check x)\,h(\check x)}. \label{eq:B-VS}\tag{60}\] Note that this system is fully characterized by two acoustic profiles: the local wavenumber \(\check k\) (or equivalently the local wave speed \(\check c(\check x) = \sqrt{g(\check x)/h(\check x)}\)) and the local impedance \(\check Z\). In a homogeneous medium, right- and left-going waves satisfy \(\Sigma = \mp\check Z V\); accordingly, we define the local wave amplitudes \[a = \tfrac12\Bigl( \sqrt{\check Z}\, V - \Sigma/\sqrt{\check Z} \Bigr) , \qquad b = \tfrac12\Bigl( \sqrt{\check Z}\, V + \Sigma/\sqrt{\check Z} \Bigr) , \label{eq:B-ab}\tag{61}\] normalized such that \(|a|^2 - |b|^2\) equals (twice) the time-averaged rightward power flux. Differentiating 61 and using 60 yields the first-order system \[a' - i\check k\,a = \frac{\check Z'}{2\check Z}\,b , \qquad b' + i\check k\,b = \frac{\check Z'}{2\check Z}\,a . \label{eq:B-coupled}\tag{62}\] which is still exact and fully equivalent to 59 . In each equation of 62 , the term proportional to \(i\check k\) rotates the phase of the corresponding amplitude at the local rate \(\check k(\check x)\) while conserving its modulus, and describes transmission; the right-hand sides, with coefficient \(\check Z'/2\check Z\), couple the two amplitudes and are the only mechanism converting right- into left-going waves, i.e., generating reflection. Moreover, \(|a|^2 - |b|^2\) is exactly conserved. Since the coupling vanishes for \(\check Z'=0\), a medium with constant impedance is reflectionless regardless of its wave-speed profile.

The WKB (adiabatic) approximation consists in neglecting the coupling wherever it is slow compared to the phase rotation: factoring out the accumulated phase, the reflected contributions generated over one local wavelength then cancel by destructive interference [51]. The validity ratio of the two rates, \(\bigl(|\check Z'(\check x)|/2\check Z(\check x)\bigr)/\check k(\check x)\), is (up to the factor of two) the quantity \(\delta\) introduced in 21 .

5.2 Transfer matrix method↩︎

For quantitative computations we solve 61 exactly by the TMM, which approximates the heterogeneous medium with a stack of \(n\) homogeneous layers with constant material properties \(g_l\), \(h_l\); for more details on the TMM applied to elastic waves see e.g. [68]. Within each layer \(l\) the coupling vanishes and the amplitudes only accumulate phase, so that the solution reads \[U_l(\check{x})=A_l\,e^{i\check{k}_l\check{x}}+B_l\,e^{-i\check{k}_l\check{x}}, \qquad \check{k}_l=\check{\omega}\sqrt{h_l/g_l} \quad\text{,}\] where \(A_l\) and \(B_l\) are the amplitudes of the rightward and leftward waves, respectively. By means of continuity of displacement and stress across the interface between layers \(l\) and \(l{+}1\), the discontinuity matrix \[\label{discont95matrix} \mathbf{D}(\check{Z}_l,\check{Z}_{l+1})=\frac{1}{2} \begin{bmatrix} 1+\check{Z}_{l+1}/\check{Z}_l & 1-\check{Z}_{l+1}/\check{Z}_l\\[2pt] 1-\check{Z}_{l+1}/\check{Z}_l & 1+\check{Z}_{l+1}/\check{Z}_l \end{bmatrix}\tag{63}\] relates the amplitudes of the two layers, with the reflection and transmission at the interface being governed by the acoustic impedances \(\check{Z}_l=\sqrt{g_l h_l}\); for a small impedance contrast, its off-diagonal entries linearize to \(\Delta\check Z/2\check Z\), recovering the coupling terms of 62 . Propagation across a layer of thickness \(\Delta\check{x}_l\) is described by the propagation matrix \[\label{propag95matrix} \mathbf{P}(\check{k}_l,\Delta\check{x}_l)= \begin{bmatrix} e^{i\check{k}_l\Delta\check{x}_l} & 0\\ 0 & e^{-i\check{k}_l\Delta\check{x}_l} \end{bmatrix} \quad\text{,}\tag{64}\] which integrates the propagation terms of 62 within a layer. The global transfer matrix is assembled as \[\label{global95transfer95matrix} \mathbf{T}(\check{\omega})=\mathbf{D}(\check{Z}_0,\check{Z}_1)\, \prod_{l=1}^{n-1}\mathbf{P}(\check{k}_l,\Delta\check{x}_l)\, \mathbf{D}(\check{Z}_l,\check{Z}_{l+1}) \quad\text{,}\tag{65}\] and links the wave amplitudes in the two semi-infinite homogeneous media bounding the layer stack – the undamaged material on either side of the regularized crack – via \[\begin{bmatrix}A_\text{inc}\\ A_\text{ref}\end{bmatrix} =\mathbf{T}(\check{\omega})\begin{bmatrix}A_\text{tra}\\ 0\end{bmatrix} \quad\text{,}\] with \(A_\text{inc}\), \(A_\text{ref}\) and \(A_\text{tra}\) being the amplitudes of the incoming, reflected and transmitted waves, respectively. Accounting for the impedance mismatch between the incident and transmitted media, the frequency-dependent power coefficients follow as \[t(\check{\omega})=\frac{\check{Z}_n}{\check{Z}_0}\left|\frac{1}{T_{11}(\check{\omega})}\right|^{2} \quad\text{, and}\quad r(\check{\omega})=\left|\frac{T_{21}(\check{\omega})}{T_{11}(\check{\omega})}\right|^{2} \quad\text{.}\] The results in Figs. 2 and 8 were obtained with \(n=10^{5}\) layers and \(10^{3}\) frequencies, expressed through the dimensionless ratio \(\ell/\lambda\).

5.3 Finite element implementation↩︎

We use FEniCSx v0.9.0 [69]–[71] and PETSc release 3.24.5 [72]–[74] for the implementation of the models in this work. The starting point is the respective action functional, with the temporal approximations according to either Newmark-\(\beta\) or generalized-\(\alpha\) inserted already. We use UFL [75] to automatically derive the residuals and their linearization, which are passed to PETSc through the petsc4py interface [76].

While some previous contributions use a monolithic solution scheme to solve the time-discrete equations of dynamic phase-field fracture in a fully coupled way [15], [22], [38], we adopt a staggered solution scheme as also done in [35], [41], [77]. As described in [56], this means that we solve alternately the \(\boldsymbol{u}\)-problem while keeping \(\alpha\) fixed, and the \(\alpha\)-problem while keeping the last converged \(\boldsymbol{u}\) fixed. We solve both subproblems with Newton’s method until convergence, while we employ the reduced-space active set strategy based on Newton’s method [78] for the damage problem to enforce the constraints. For the linear solve in each Newton iteration, we use the MUMPS package [79], [80] with the Cholesky factorization. For the model with density degradation, we use the \(LU\) factorization instead. We determine convergence of both subproblems as \(\|\mathbf{R}_{\mathbf{u}} (\mathbf{u}, \boldsymbol{\alpha}^{i-1})\|_2 \leq 10^{-7}\) and \(\|\mathbf{R}_{\boldsymbol{\alpha}} (\mathbf{u}^i, \boldsymbol{\alpha})\|_2 \leq 10^{-7}\), while we check the convergence of the staggered scheme as \(\|\mathbf{R}_{\mathbf{u}} (\mathbf{u}^i, \boldsymbol{\alpha}^{i})\|_2 \leq 10^{-5}\), where \(i\) is the staggered iteration. All our implementations are available at https://github.com/jonas-heinzmann/phase_field_dynamics.

5.4 Compact reformulation of the cohesive strain energy density↩︎

To derive a compact reformulation for the cohesive model, we seek a stationary point of the action functional with respect to \(\boldsymbol{\eta}\), respecting the constraint \(\text{tr}(\boldsymbol{\eta}) \geq 0\). Since \(\boldsymbol{\eta}\) carries no inertia, we can optimize the strain energy functional directly. By the normality rule [40], the sign of \(\text{tr}(\boldsymbol{\eta})\) coincides with the sign of the pressure \(p = \tfrac{1}{d} \text{tr}(\boldsymbol{\sigma})\), allowing us to distinguish the following two cases.

5.4.0.1 Non-negative pressure

For \(p \geq 0\), stationarity of \(\psi^{\text{c}}\) with respect to \(\boldsymbol{\eta}\) yields \[\label{eq:CF95psi95voldev95postr} \frac{\partial}{\partial \boldsymbol{\eta}} \left( \frac{1}{2} \mathbb{C}_0 (\boldsymbol{\varepsilon} - \boldsymbol{\eta}^\star) \cdot (\boldsymbol{\varepsilon} - \boldsymbol{\eta}^\star) + a(\alpha) \sqrt{w_{\text{c}}} \sqrt{ \mathbb{C}_0 \boldsymbol{\eta}^\star \cdot \boldsymbol{\eta}^\star} \right) = \boldsymbol{0} \quad\Rightarrow\quad \boldsymbol{\varepsilon} - \boldsymbol{\eta}^\star = a(\alpha) \sqrt{w_{\text{c}}} \frac{\boldsymbol{\eta}^\star}{\sqrt{ \mathbb{C}_0 \boldsymbol{\eta}^\star \cdot \boldsymbol{\eta}^\star}} \quad\text{,}\tag{66}\] where \(\boldsymbol{\eta}^\star\) is the optimum. Noting the coaxiality between \(\boldsymbol{\eta}^\star\) and \(\boldsymbol{\varepsilon}\), we can write \(\boldsymbol{\eta}^\star = q(\alpha) \boldsymbol{\varepsilon}\) with \[\boldsymbol{\varepsilon} \left( 1 - q(\alpha) - \frac{a(\alpha) \sqrt{w_{\text{c}}}}{\sqrt{ \mathbb{C}_0 \boldsymbol{\varepsilon} \cdot \boldsymbol{\varepsilon}}} \right) = \boldsymbol{0} \qquad\Rightarrow\qquad q(\alpha) = 1 - \frac{a(\alpha) \sqrt{w_{\text{c}}}}{\sqrt{ \mathbb{C}_0 \boldsymbol{\varepsilon} \cdot \boldsymbol{\varepsilon}}} \qquad\text{.}\] Requiring \(q(\alpha) \geq 0\) yields \(\psi_0 (\boldsymbol{\varepsilon}) \geq \frac{1}{2} a(\alpha)^2 w_{\text{c}}\). Substituting in the cohesive strain energy density gives \[\hat{\psi}^{\text{c}} (\boldsymbol{\varepsilon},\alpha) = a(\alpha) \sqrt{w_{\text{c}}} \sqrt{2\psi_0(\boldsymbol{\varepsilon})} - \frac{1}{2} a(\alpha)^2 w_{\text{c}} \qquad \text{if}\,\,\, \text{tr}(\boldsymbol{\varepsilon}) \geq 0, \psi_0 (\boldsymbol{\varepsilon}) \geq \frac{1}{2} a(\alpha)^2 w_{\text{c}}.\]

5.4.0.2 Negative pressure

For \(p < 0\), the eigenstrain potential 35 enforces \(\text{tr}(\boldsymbol{\eta}) = 0\), reducing 34 to \[\label{CF95psi95voldev95treta950} \psi^{\text{c}} (\boldsymbol{\varepsilon}, \boldsymbol{\eta}, \alpha) |_{\text{tr}(\boldsymbol{\eta}) = 0} = \frac{\kappa_0}{2} \left( \text{tr}(\boldsymbol{\varepsilon}) \right)^2 + 2 \mu_0 \left( \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| - \left| \boldsymbol{\eta}_{\text{dev}} \right| \right)^2 + a(\alpha) \sqrt{w_{\text{c}}} \sqrt{2 \mu_0} \left| \boldsymbol{\eta}_{\text{dev}} \right| \quad\text{.}\tag{67}\] Stationarity with respect to \(\left| \boldsymbol{\eta}_{\text{dev}} \right|\) then gives the optimum \[\label{eq:CF95compact95trneg95opt} \left| \boldsymbol{\eta}_{\text{dev}} \right|^\star = \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| - a(\alpha) \sqrt{\frac{w_{\text{c}}}{2\mu_0}} \quad\text{.}\tag{68}\] Non-negativity of the norm requires \(\left| \boldsymbol{\varepsilon}_{\text{dev}} \right| \geq a(\alpha) \sqrt{w_{\text{c}}/(2\mu)}\) for this branch. Combining the two cases and adding the residual energy density \(\epsilon \psi_0\) to avert numerical issues yields the condensed form of the strain energy density reported in the main body, 44 .

From it, the stress 40 can now be evaluated directly from \(\hat{\psi}^{\text{c}}\): \[\begin{equation} p (\boldsymbol{\varepsilon}, \alpha) = \frac{\partial \hat{\psi}^{\text{c}} (\boldsymbol{\varepsilon}, \alpha)}{\partial \text{tr}(\boldsymbol{\varepsilon})} = \begin{cases} \kappa_0 \left( \frac{a(\alpha) \sqrt{w_{\text{c}}}}{\sqrt{2\psi_0 (\boldsymbol{\varepsilon})}} + \epsilon \right) \text{tr}(\boldsymbol{\varepsilon}) &\text{if}\,\,\, \text{tr}(\boldsymbol{\varepsilon}) \geq 0, \psi_0 (\boldsymbol{\varepsilon}) \geq \frac{1}{2} a(\alpha)^2 w_{\text{c}}\\ \kappa_0 (1 + \epsilon) \text{tr}(\boldsymbol{\varepsilon}) &\text{else} \end{cases} \qquad\text{and} \end{equation} \begin{equation} \tau (\boldsymbol{\varepsilon}, \alpha) = \frac{\partial \hat{\psi}^{\text{c}} (\boldsymbol{\varepsilon}, \alpha)}{\partial \left| \boldsymbol{\varepsilon}_{\text{dev}} \right|} = \begin{cases} 2\mu_0 \left( \frac{ a(\alpha) \sqrt{w_{\text{c}}}}{\sqrt{2\psi_0 (\boldsymbol{\varepsilon})}} + \epsilon \right) \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| &\text{if}\,\,\, \text{tr}(\boldsymbol{\varepsilon}) \geq 0, \psi_0 (\boldsymbol{\varepsilon}) \geq \frac{1}{2} a(\alpha)^2 w_{\text{c}}\\ a(\alpha) \sqrt{2 \mu_0 w_{\text{c}}} + 2 \epsilon \mu_0 \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| &\text{if}\,\,\, \text{tr}(\boldsymbol{\varepsilon}) < 0, \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| \geq a(\alpha) \sqrt{\frac{w_{\text{c}}}{2\mu_0}}\\ 2\mu_0 (1 + \epsilon) \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| &\text{else} \end{cases} \qquad\text{.} \end{equation}\] The eigenstrain can be post-processed from \(\boldsymbol{\varepsilon}\) as \[\begin{equation} \text{tr}(\boldsymbol{\eta}) = \begin{cases} q(\alpha) \text{tr}(\boldsymbol{\varepsilon}) &\text{if}\,\,\, \text{tr}(\boldsymbol{\varepsilon}) \geq 0, \psi_0 (\boldsymbol{\varepsilon}) \geq \frac{1}{2} a(\alpha)^2 w_{\text{c}}\\ 0 &\text{else} \end{cases} \qquad\text{and}\qquad \end{equation} \begin{equation} \left| \boldsymbol{\eta}_{\text{dev}} \right| = \begin{cases} q(\alpha) \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| &\text{if}\,\,\, \text{tr}(\boldsymbol{\varepsilon}) \geq 0, \psi_0 (\boldsymbol{\varepsilon}) \geq \frac{1}{2} a(\alpha)^2 w_{\text{c}}\\ \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| - a(\alpha) \sqrt{\frac{w_{\text{c}}}{2\mu_0}} &\text{if}\,\,\, \text{tr}(\boldsymbol{\varepsilon}) < 0, \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| \geq a(\alpha) \sqrt{\frac{w_{\text{c}}}{2\mu_0}}\\ 0 &\text{else} \end{cases} \end{equation}\] where the residual energy density does not affect the definitions since the stationarity conditions with respect to \(\boldsymbol{\eta}\), 66 and 68 , are not affected by \(\epsilon \psi_0(\boldsymbol{\varepsilon})\).

5.5 The dynamic cohesive law↩︎

In the following, we derive the crack opening evolution for the dynamic cohesive phase-field model.

5.5.1 General solution↩︎

We consider an isolated material point at \(x=0\) at which a crack will open, adjacent to two semi-infinite half-lines representing the bulk on the left and right sides of the crack. We assume that the bulk is linearly elastic (\(\eta \equiv 0\) for all \(x \neq 0\)), i.e., the wave interacts only with a single crack, see Section 3.4 and Fig. 26a. We allow for a pre-existing phase-field crack at \(x=0\), with corresponding maximum phase-field value \(\breve{\alpha}(t)\), while the jump is initially closed, \(\llbracket u \rrbracket(t=0)=0\). The tensile loading is assumed large enough that the stress at \(x=0\) reaches the critical stress \(\sigma_{\text{c}}\) at some time, and that it provides enough energy for the crack to fully open. On the left side of the crack (\(x_L = \{x : x < 0\}\)) we consider an incoming, right-traveling wave \(u_L^{\text{inc}}\) and a potentially reflected, left-traveling wave \(u_L^{\text{ref}}\); on the right side of the crack (\(x_R = \{x : x > 0\}\)) we consider only a transmitted, right-traveling wave \(u_R^{\text{tra}}\).

In both half-spaces, the standard d’Alembert solution holds along the characteristics \(t \pm x/c_0 = \text{const}\): \[\label{eq:u95DAlambert} u_L(x_L,t) = \underbrace{f_L(t - \tfrac{x_L}{c_0})}_{u_L^{\text{inc}}(x_L,t)} + \underbrace{g_L(t + \tfrac{x_L}{c_0})}_{u_L^{\text{ref}}(x_L,t)} \qquad\text{and}\qquad u_R(x_R,t) = \underbrace{f_R(t - \tfrac{x_R}{c_0})}_{u_R^{\text{tra}}(x_R,t)} \qquad\text{,}\tag{69}\] with the arbitrary perturbations \(f_L,\,g_L,\,f_R\) determined by the initial conditions. The stress follows directly from the velocity, \[\label{eq:DAlambert95stess95velocity} \sigma_L(x_L,t) = Z_0 \left( - v_L^{\text{inc}}(x_L,t) + v_L^{\text{ref}}(x_L,t) \right) \qquad\text{and}\qquad \sigma_R(x_R,t) = - Z_0\, v_R^{\text{tra}}(x_R,t)\quad\text{,}\tag{70}\] with the (constant) bulk acoustic impedance \(Z_0 = E_0/c_0 = \rho_0 c_0 = \sqrt{E_0 \rho_0}\).

At the crack, the kinematic compatibility on the jump and the stress continuity (with the cohesive law 50 ) read \[\begin{align} \llbracket u \rrbracket(t) &= u_R^{\text{tra}}(0^+,t) - u_L^{\text{inc}}(0^-,t) - u_L^{\text{ref}}(0^-,t) \quad\text{,}\tag{71}\\ \sigma_L(0^-,t) &= \sigma_R(0^+,t) = \left(1 - \breve{\alpha}(t)\right)\sigma_{\text{c}} = \left(1 - \frac{\sigma_{\text{c}}}{2G_{\text{c}}} \llbracket u \rrbracket(t)\right)\sigma_{\text{c}} \quad\text{.}\tag{72} \end{align}\] Using 70 in the stress continuity gives \(v_R^{\text{tra}}(0^+,t) = v_L^{\text{inc}}(0^-,t) - v_L^{\text{ref}}(0^-,t)\). Inserting this into the time-integrated form of 71 eliminates the incident and transmitted contributions and yields \[\label{eq:cohesive95response95raw95jump} \llbracket u \rrbracket(t) = -2 \int_{0}^{t} v_L^{\text{ref}}(0^-,\tau)\,\mathrm{d}\tau \quad\text{.}\tag{73}\] Owing to the reversibility of the eigenstrain, we have to distinguish the four phases of Section 3.5, illustrated in Fig. 26b-i.

5.5.1.1 Phase I: \(\dot{\eta} = 0\), \(\dot{\alpha} = 0\)

For \(0 \leq t \leq t_{\text{I}}\), the stress at \(x=0\) remains below \(\sigma_{\text{c}}\), so the incoming wave is fully transmitted with no reflection: \[u_L^{\text{ref}}(x_L,t) = 0 \quad\text{,}\quad \llbracket u \rrbracket(t) = 0 \quad\text{,}\quad u_R^{\text{tra}}(0^+,t) = u_L^{\text{inc}}(0^-,t)\quad\text{,}\] with the maximum phase field still at its initial value, \(\breve{\alpha}(t_{\text{I}}) = \breve{\alpha}(0)\).

5.5.1.2 Phase II: \(\dot{\eta} > 0\), \(\dot{\alpha} = 0\)

For \(t_{\text{I}} \leq t \leq t_{\text{II}}\), the jump reopens to the value associated with the (irreversible) phase-field value \(\breve{\alpha}(t_{\text{I}})\). The strength criterion thus holds with \(\sigma_L(0^-,t) = (1 - \breve{\alpha}(t_{\text{I}}))\sigma_{\text{c}}\), and 70 together with the initial condition \(u_L^{\text{ref}}(0^-,t_{\text{I}}) = 0\) yields \[\begin{align} v_L^{\text{ref}}(0^-,t) &= v_L^{\text{inc}}(0^-,t) + \frac{1 - \breve{\alpha}(t_{\text{I}})}{Z_0}\,\sigma_{\text{c}}\quad\text{,} \\ u_L^{\text{ref}}(0^-,t) &= u_L^{\text{inc}}(0^-,t) - u_L^{\text{inc}}(0^-,t_{\text{I}}) + \frac{1 - \breve{\alpha}(t_{\text{I}})}{Z_0}\,\sigma_{\text{c}}\,(t - t_{\text{I}})\quad\text{.} \end{align}\] The jump and the transmitted wave then follow from 73 and 71 as \[\begin{align} \llbracket u \rrbracket(t) &= -2\,\frac{1 - \breve{\alpha}(t_{\text{I}})}{Z_0}\,\sigma_{\text{c}}\,(t - t_{\text{I}}) - 2 u_L^{\text{inc}}(0^-,t) + 2 u_L^{\text{inc}}(0^-,t_{\text{I}})\quad\text{,} \\ u_R^{\text{tra}}(0^+,t) &= -\frac{1 - \breve{\alpha}(t_{\text{I}})}{Z_0}\,\sigma_{\text{c}}\,(t - t_{\text{I}}) + u_L^{\text{inc}}(0^-,t_{\text{I}})\quad\text{.} \end{align}\] Phase II ends at \(t_{\text{II}}\), defined by the jump reaching the value at which the phase field is about to evolve, \(\llbracket u \rrbracket(t_{\text{II}}) = \tfrac{2G_{\text{c}}}{\sigma_{\text{c}}}\breve{\alpha}(t_{\text{I}})\).

5.5.1.3 Phase III: \(\dot{\eta} > 0\), \(\dot{\alpha} > 0\)

For \(t_{\text{II}} \leq t \leq t_{\text{III}}\), the cohesive stress decreases as the jump grows. Substituting 73 into 72 and differentiating in time yields, after introducing the Irwin length \(\ell_{\text{ch}} = E_0 G_{\text{c}}/\sigma_{\text{c}}^2\) via \(\sigma_{\text{c}}^2/(Z_0 G_{\text{c}}) = c_0/\ell_{\text{ch}}\), the linear ODE \[\label{eq:opening95evolution95ODE} \dot{v}_L^{\text{ref}}(0^-,t) - \frac{c_0}{\ell_{\text{ch}}}\,v_L^{\text{ref}}(0^-,t) = \dot{v}_L^{\text{inc}}(0^-,t) \quad\text{,}\tag{74}\] which reveals a dependence of the solution on the ratio of wave speed to Irwin length and on the incoming wave. With the initial condition at \(t_{\text{II}}\), its solution is \[\label{eq:cohesive95response95ODE95solution} v_L^{\text{ref}}(0^-,t) = \exp\!\left(\tfrac{c_0}{\ell_{\text{ch}}}(t-t_{\text{II}})\right) v_L^{\text{ref}}(0^-,t_{\text{II}}) + \int_{t_{\text{II}}}^{t} \exp\!\left(\tfrac{c_0}{\ell_{\text{ch}}}(t-\tau)\right) \dot{v}_L^{\text{inc}}(0^-,\tau)\,\mathrm{d}\tau\quad\text{.}\tag{75}\] Phase III ends at \(t_{\text{III}}\) when \(\breve{\alpha}\) reaches one, marking full opening. All other quantities follow from 75 . To avoid the explicit evaluation of the integral, we directly integrate the ODE 74 between \(t_{\text{II}}\) and \(t\), which together with 73 yields \[\begin{align} u_L^{\text{ref}}(0^-,t) &= u_L^{\text{ref}}(0^-,t_{\text{II}}) + \tfrac{\ell_{\text{ch}}}{c_0}\Bigl[v_L^{\text{ref}}(0^-,t) - v_L^{\text{ref}}(0^-,t_{\text{II}}) - \bigl(v_L^{\text{inc}}(0^-,t) - v_L^{\text{inc}}(0^-,t_{\text{II}})\bigr)\Bigr]\quad\text{,} \\ \llbracket u \rrbracket(t) &= -2 u_L^{\text{ref}}(0^-,t_{\text{II}}) - 2\tfrac{\ell_{\text{ch}}}{c_0}\Bigl[v_L^{\text{ref}}(0^-,t) - v_L^{\text{ref}}(0^-,t_{\text{II}}) - \bigl(v_L^{\text{inc}}(0^-,t) - v_L^{\text{inc}}(0^-,t_{\text{II}})\bigr)\Bigr]\quad\text{,} \end{align}\] the latter being the expression 56 reported in the main body.

5.5.1.4 Phase IV: \(\dot{\eta} \gtreqless 0\), \(\dot{\alpha} = 0\)

For \(t \geq t_{\text{III}}\), the cohesive stress is zero and the crack behaves as a free surface, \[v_L^{\text{ref}}(x_L,t) = v_L^{\text{inc}}(x_L,t)\quad\text{,}\qquad u_L^{\text{ref}}(x_L,t) = u_L^{\text{ref}}(x_L,t_{\text{III}}) + u_L^{\text{inc}}(x_L,t) - u_L^{\text{inc}}(x_L,t_{\text{III}})\quad\text{.}\]

In all phases, the reflected and transmitted wave solutions in the bulk are obtained from these crack-located expressions by propagation along the characteristics.

5.5.2 Solution for a harmonic half-sine wave↩︎

We finally evaluate the above for the sinusoidal half-pulse 24 used in our numerical examples. With \(\tilde{T} = \lambda/c_0\), the incoming wave and its velocity are \[u_L^{\text{inc}}(x_L,t) = \begin{cases} 0 & t - \tfrac{x_L}{c_0} < 0\\ \frac{\tilde{\sigma}}{Z_0}\frac{\tilde{T}}{2\pi}\!\left(\cos\!\left(2\pi \frac{t - x_L/c_0}{\tilde{T}}\right) - 1\right) & t - \tfrac{x_L}{c_0} \in [0,\tfrac{\tilde{T}}{2}]\\ -\frac{\tilde{\sigma}}{Z_0}\frac{\tilde{T}}{\pi} & t - \tfrac{x_L}{c_0} > \tfrac{\tilde{T}}{2} \end{cases} \text{,}\quad v_L^{\text{inc}}(x_L,t) = \begin{cases} -\frac{\tilde{\sigma}}{Z_0}\sin\!\left(2\pi \frac{t - x_L/c_0}{\tilde{T}}\right) & t - \tfrac{x_L}{c_0} \in [0,\tfrac{\tilde{T}}{2}]\\ 0 & \text{else} \end{cases}\quad\text{.}\] The reflected wave velocity and displacement at any \(x_L\) then follow from the general derivation as \[v_L^{\text{ref}}(x_L,t) = \begin{cases} 0 & t < t_{\text{I}}\\ v_L^{\text{inc}}\!\left(0,t+\tfrac{x_L}{c_0}\right) + \bigl(1-\breve{\alpha}(t_{\text{I}})\bigr)\frac{\sigma_{\text{c}}}{Z_0} & t_{\text{I}} \leq t < t_{\text{II}}\\ \exp\!\left(\tfrac{c_0}{\ell_{\text{ch}}}\bigl(t+\tfrac{x_L}{c_0}-t_{\text{II}}\bigr)\right)\left[v_L^{\text{ref}}(0,t_{\text{II}}) - \tilde{v}^{\text{ref}}(t_{\text{II}})\right] + \tilde{v}^{\text{ref}}\!\left(t+\tfrac{x_L}{c_0}\right) & t_{\text{II}} \leq t < t_{\text{III}}\\ v_L^{\text{inc}}\!\left(0,t+\tfrac{x_L}{c_0}\right) & t_{\text{III}} \leq t \end{cases}\] with \[\tilde{v}^{\text{ref}}(\tau) = \frac{\dfrac{2\pi \tilde{\sigma}}{\tilde{T} Z_0}}{\left(\dfrac{c_0}{\ell_{\text{ch}}}\right)^{2} + \left(\dfrac{2\pi}{\tilde{T}}\right)^{2}}\left[\frac{c_0}{\ell_{\text{ch}}}\cos\!\left(2\pi\tfrac{\tau}{\tilde{T}}\right) - \frac{2\pi}{\tilde{T}}\sin\!\left(2\pi\tfrac{\tau}{\tilde{T}}\right)\right]\quad\text{,}\] and \[u_L^{\text{ref}}(x_L,t) = \begin{cases} 0 & t < t_{\text{I}}\\ u_L^{\text{inc}}(0,t+\tfrac{x_L}{c_0}) - u_L^{\text{inc}}(0,t_{\text{I}}) + (1-\breve{\alpha}(t_{\text{I}}))\frac{\sigma_{\text{c}}}{Z_0}\bigl(t+\tfrac{x_L}{c_0}-t_{\text{I}}\bigr) & t_{\text{I}} \leq t < t_{\text{II}}\\ u_L^{\text{ref}}(0,t_{\text{II}}) + \frac{\ell_{\text{ch}}}{c_0}\Bigl[v_L^{\text{ref}}(0,t+\tfrac{x_L}{c_0}) - v_L^{\text{ref}}(0,t_{\text{II}}) - \bigl(v_L^{\text{inc}}(0,t+\tfrac{x_L}{c_0}) - v_L^{\text{inc}}(0,t_{\text{II}})\bigr)\Bigr] & t_{\text{II}} \leq t < t_{\text{III}}\\ u_L^{\text{ref}}(0,t_{\text{III}}) + u_L^{\text{inc}}(0,t+\tfrac{x_L}{c_0}) - u_L^{\text{inc}}(0,t_{\text{III}}) & t_{\text{III}} \leq t \end{cases}\] from which all remaining quantities of interest can be obtained. The end of Phase I is given in closed form by \(t_{\text{I}} = \tfrac{\tilde{T}}{2\pi}\sin^{-1}\!\left((1-\breve{\alpha}(0))\,\sigma_{\text{c}}/\tilde{\sigma}\right)\), while \(t_{\text{II}}\) and \(t_{\text{III}}\) must be obtained numerically.

In Fig. 27, we show a comparison between these theoretical results and the FEM results from Fig. 15a: the chosen fracture toughnesses yield \(c_0/\ell_{\text{ch}} = 2.95 \cdot 10^5\) s\(^{-1}\) for \(G_{\text{c}} = 0.01\) N/mm and \(c_0/\ell_{\text{ch}} = 1.18 \cdot 10^5\) s\(^{-1}\) for \(G_{\text{c}} = 0.025\) N/mm. Clearly, analytical and FEM results coincide; the case producing a shock wave is excluded from this comparison.

References↩︎

---

1. Recall from 18 that the length of the phase-field support is \(4\ell\).↩︎

2. Originally, the deviatoric contribution is written as \(\mu_0 \left| \boldsymbol{\varepsilon}_{\text{dev}} - \boldsymbol{\eta}_{\text{dev}} \right|^2\). However, since the elastic domain depends only on the spherical and deviatoric components of \(\boldsymbol{\eta}\) and is axisymmetric about the purely volumetric direction, it can be replaced by its lower bound \(\mu_0 \left( \left| \boldsymbol{\varepsilon}_{\text{dev}} \right| - \left| \boldsymbol{\eta}_{\text{dev}} \right| \right)^2\) which is attained when \(\boldsymbol{\varepsilon}_{\text{dev}}\) and \(\boldsymbol{\eta}_{\text{dev}}\) are co-aligned (see [40]). This substitution, originally derived from energy minimization in the quasi-static setting, remains valid in dynamics: the original \(\psi\) (with the deviatoric term \(\mu_0 \left| \boldsymbol{\varepsilon}_{\text{dev}} - \boldsymbol{\eta}_{\text{dev}} \right|^2\)) is convex in \(\boldsymbol{\eta}\) and the eigenstrain carries no inertia, so the point-wise stationarity condition of the action functional with respect to \(\boldsymbol{\eta}\) coincides with the first-order minimality condition of the quasi-static problem – and, by convexity, this stationary point is in fact a minimum.↩︎