Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture Analyses?

Abstract

Numerical techniques to simulate crack propagation can roughly be divided into sharp and diffuse interface methods. Two prominent approaches to quantitative dynamic fracture analysis are compared here. Specifically, an adaptive cohesive element technique and a phase-field fracture approach are applied to simulate Hopkinson bar experiments on the fracture toughness of high-performance concrete. The experimental results are validated numerically in the sense of an inverse analysis. Both methods allow predictive numerical simulations of crack growth with an a priori unknown path and determine the related material parameter in a quantitative manner. Reliability, precision, and numerical costs differ however.

1. Introduction

One of the main challenges in computational mechanics is the prediction of cracks and fragmentation in dynamic fracture. There are high demands on the modeling side, but mainly the complicated structure and the nonregular behavior of the cracks turn numerical simulations into a difficult task. Every crack in a solid forms a new surface of a priori unknown position, which needs to be identified. Different discretization techniques have been developed to solve such problems, for example the cohesive element technique [1, 2, 3], the extended finite element method [4, 5], eroded finite elements or eigenfracture strategies [6, 7], and phase-field approaches [8, 9, 10, 11, 12, 13].

The numerical techniques to treat the moving boundary problem of crack propagation can roughly be divided into two different strategies: sharp interface and diffuse interface modeling. The sharp interface approach describes a crack as a new boundary ${\mathrm{\Gamma }}_C\left(t\right)\subset \mathrm{\partial \Omega }$ in a solid of domain $\mathrm{\Omega }$ undergoing a deformation $\mathit{\chi }\left(\mathit{x},t\right):\mathrm{\Omega }\times \left[0,t_{\text{total}}\right]\to {\mathrm{ℝ}}^3$ in a time $t_{\text{total}}$ . For a known crack path this is the natural way to capture the mechanics of fracture. In dynamic fracture, however, with a priori unknown ways of crack propagation, kinks, and branching, sophisticated tracking methods need to be employed to localize the boundaries by the position and to enforce (free) boundary conditions along the moving crack. Unfortunately, such surface tracking methods require a high numerical effort and tend to fail for great changes in the topology of the solid, such as the fragmentation into small particles and their further movement.

An alternative way to describe moving boundaries are diffuse interface models where the cracks are smeared over a small but finite length $\epsilon$ . Here an additional field $s\left(\mathit{x},t\right):\mathrm{\Omega }\times \left[0,t_{\text{total}}\right]\to \mathrm{ℝ}$ characterizes the state of the material and marks the intact or broken state. The set of evolving crack surfaces is replaced by a crack-surface density $\gamma$ , which is typically a function of the marker field s and its gradient. This crack-surface density function allows an approximation of the moving crack boundaries over the body’s domain,

$${\int }_{{\mathrm{\Gamma }}_C\left(t\right)}\mathrm{d\Gamma }\approx {\int }_{\mathrm{\Omega }}\gamma \left(t\right)\mathrm{d\Omega }.$$ E1

By $\gamma$ being not only a function of field $s\left(\mathit{x},t\right)$ but also its gradient $\nabla s$ , it regularizes (or diffuses) the local jump of a crack. The net effect of this regularization is to eliminate spurious mesh-dependencies that afflict naive damage schemes. The potential energy of such a regularized cracking body corresponds to the well-known Ambrosio-Tortorelli functional of continuum damage mechanics [14]. Therefore, the diffuse interface approach can be seen as a gradient damage model with the major difference that the order parameter s indicates the material to be either intact ( $s=1$ ) or broken ( $s=0$ ). Intermediate states are not physically meaningful. However, it is still an open question how reliable such a diffuse approach can quantify the mechanics of fracture.

Here we compare a sharp interface method with crack tracking algorithm and a diffuse interface method for its usability in material identification. Background for our comparison are our experimental investigations on the fracture toughness of ultra-high performance concrete (UHPC). Specifically, we use the cohesive element technique and the phase-field fracture approach to simulate spalling experiments performed with concrete specimen in a Hopkinson-Bar (HB) setup.

UHPC is a class of advanced cementitious-based composites whose mechanical strength and durability surpass classical concrete. Typically, UHPC composites are fine grained, almost homogeneous mixtures of small aggregates of cement, a certain amount of silica, other supplements, and a low water content—and so they are more similar to brittle ceramics than to construction concrete. UHPCs are still under development and in order to optimize their composition mechanical tests have to provide material data. Hereby classical experiments determine the concrete’s elasticity as well as its compressive and flexural strength under static loading conditions. For the dynamic properties, however, such as dynamic tensile resistance and fracture energy, it is more complicated to ensure reproducible test conditions. Here numerical simulations in the sense of an inverse analysis are helpful to evaluate the reliability of the obtained material data.

HB spalling experiments are test arrangements to determine the failure strength of brittle materials, see [15, 16, 17, 18, 19]. In these tests the experimental setup of a classical HB is modified in such a way, that the induced pressure impulse is transmitted via an incident bar into the specimen, see Figure 1. Within the specimen a superposition of transmitted and reflected waves determines the stress state. For details of the experimental work we refer to another work [20], here we just use the experimental setup to compare two numerical techniques employed for quantitative analysis. Specifically, for fracture parameter identification we need: (i) numerical methods that are able to find the crack position dependent on the external load and the material parameter of the specimen; (ii) the pressure wave and the stress distribution in the fast (cracking) specimen; and (iii) we need to quantify the fracture energy and the critical energy release rate of the material.

Figure 1.

Illustration of the HB-spallation test setup where the pressure tank accelerates the impactor and the impact-induced wave propagates trough the specimen (left) and UHPC specimen cracked after wave reflection (right).

The remaining paper is organized as follows. In the next section we provide shortly the governing equations of elasto-dynamics and fracture mechanics. Then we introduce the cohesive element technique in Section 3 and the phase-field fracture method in Section 4. Both sections conclude with a short study on the influence of the relevant model parameters. In Section 5 the simulations of the HB spalling experiment are described in detail and a range of values for the fracture parameters is derived. The inverse analysis is presented in Section 6. Here we provide several numerical simulations and evaluate both methods. Such a quantitative comparison is new and has not yet been presented before. In particular, predictive applications of the phase-field approach to fracture are not common by now. A summary of the pros and cons of both methods in Section 7 concludes the paper.

2. Governing equations

We consider a body of domain $\mathrm{\Omega }\subset {\mathrm{ℝ}}^3$ with external boundary $\mathrm{\Gamma }\equiv \mathrm{\partial \Omega }$ . The body’s displacement field at point $\mathit{x}$ and time $t$ is denoted by $\mathit{u}\left(\mathit{x},t\right)$ ; its velocity and acceleration fields are $\mathit{v}=\dot{u}$ and $\mathit{a}=\ddot{\mathit{u}}$ . A crack splits the body into subbodies ${\mathrm{\Omega }}^{+}\cup {\mathrm{\Omega }}^{-}=\mathrm{\Omega }$ and induces a displacement jump on ${\mathrm{\Gamma }}_C$ as $〚\mathit{u}〛={\mathit{u}}^{+}-{\mathit{u}}^{-}$ . The displacements satisfy the Dirichlet boundary conditions $\mathit{u}=\overline{\mathit{u}}$ at ${\mathrm{\Gamma }}_1$ . The body is loaded with traction $\overline{\mathit{t}}$ at boundary ${\mathrm{\Gamma }}_2$ ; it holds ${\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2\subset \mathrm{\Gamma }\backslash {\mathrm{\Gamma }}_C$ .

2.1 Elasto-dynamics

Linear-elastic material is presumed to follow Hooke’s law with elastic strain energy density,

$${\mathrm{\Psi }}^e=\frac{1}{2}\lambda {\left(\mathrm{tr}\mathit{\epsilon }\right)}^2+\mu \mathit{\epsilon }:\mathit{\epsilon }=\frac{1}{2}\mathit{\epsilon }:\mathbb{C}:\mathit{\epsilon }$$ E2

where the Lamé material parameters $\lambda =\mathit{E\nu }/\left(\left(1+\nu \right)\left(1-2\nu \right)\right)$ and $\mu =E/\left(2\left(1+\nu \right)\right)$ are formulated with Young’s modulus $E$ and Poisson number $\nu$ . Equivalently we use the Hookean material tensor $\mathbb{C}$ with components $C_{\mathit{\text{ijkl}}}={\mathit{\lambda \delta }}_{\mathit{ij}}{\delta }_{\mathit{kl}}+\mu \left({\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}+{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}}\right)$ . The strain tensor $\mathit{\epsilon }$ describes small deformations $\mathit{\epsilon }\left(\mathit{u}\right)=1/2\left(\mathbf{\text{∇}}\mathit{u}+\mathbf{\text{∇}}{\mathit{u}}^T\right)$ ; the mechanical stress tensor follows as $\mathit{\sigma }=\partial \mathit{\epsilon }{\mathrm{\Psi }}^e=\mathbb{C}\mathit{\epsilon }$ . It holds the balance of linear momentum,

$$div\mathit{\sigma }+\overline{\mathit{b}}=\rho \mathit{a}\text{in}\mathrm{\Omega }$$ E3

where $\rho$ is the mass density and $\overline{\mathit{b}}$ a prescribed body force density. Boundary conditions are prescribed as

$$\mathit{u}=\overline{\mathit{u}}\mathrm{on}{\mathrm{\Gamma }}_1,\mathit{\sigma n}=\overline{\mathit{t}}\mathrm{on}{\mathrm{\Gamma }}_2,〚\mathit{\sigma n}〛=〚\mathit{t}〛=\mathbf{0}\mathrm{on}{\mathrm{\Gamma }}_C,$$ E4

with normal vector $\mathit{n}$ and traction $\mathit{t}$ ; $〚\cdot 〛$ denotes a jump. Eq. (4) ₃ corresponds to traction-free crack boundaries. Additional initial conditions may apply.

2.2 Fracture mechanics

Let the evolving internal cracks be represented by a set of boundaries ${\mathrm{\Gamma }}_C\left(t\right)$ . According to the linear-elastic fracture theory of Griffith and Irwin [21, 22], a material fails upon attainment of a critical surface-energy density. The crack growth corresponds to the creation of new surfaces and hence the internal work of the body is composed of

$$W_{\mathrm{int}}\left(\mathit{u},t\right)={\int }_{\mathrm{\Omega }}{\mathrm{\Psi }}^e\left(\mathit{u}\right)\mathrm{d\Omega }+{\int }_{\mathrm{\Gamma }\left(t\right)}{\mathcal{G}}_c\mathrm{d\Gamma },$$ E5

where ${\mathcal{G}}_c$ is commonly known as Griffith’s critical energy release rate (Griffith energy). An optimum of Eq. (5) corresponds to crack growth.

The specific energy $\mathcal{G}$ corresponds to the energy W dissipated per unit newly created surface area $A$ , $\mathcal{G}=\mathrm{d}W/\mathrm{d}A$ . The Griffith fracture criterion states that a crack will grow when the available energy release rate is greater than or equal to a critical value, $\mathcal{G}\ge {\mathcal{G}}_c$ , which is related to the classical stress intensity factors of linear-elastic fracture mechanics given by $\mathcal{G}=\left(K_I^2+K_{\mathit{II}}^2\right)/E^{\prime }+K_{\mathit{III}}^2/\left(2\mu \right)$ with $E^{\prime }=E$ regarding the plane stress state and otherwise $E^{\prime }=E/\left(1-{\nu }^2\right)$ . Thus, Griffith’s criterion can also be written with the critical stress intensity factor $K_c$ as ${\mathcal{G}}_c=K_c^2/E^{\prime }$ . The single summands of $\mathcal{G}$ may also be understood as mode-dependent Griffith energy densities ${\mathcal{G}}_I$ , ${\mathcal{G}}_{\mathit{II}}$ , and ${\mathcal{G}}_{\mathit{III}}$ . The indices denote the fracture modes of crack opening, in-plane and out-of-plane sliding. For a brittle elastic material the $\mathcal{J}$ -integral is equal to the strain energy release rate, that is, it holds the identity ${\mathcal{J}}_c={\mathcal{G}}_c$ .

Another fracture criterion is the crack tip opening displacement with critical value ${\delta }_c$ . The underlying theory of Dugdale [23] and Barenblatt [24] explains crack growth as loss of cohesion in a cohesive zone. If the decohesion is modeled in a nonlinear way, that is, for a given relation between the vector of cohesive traction $\mathit{t}$ and crack opening $\mathit{\delta }$ , a standard application of the $\mathcal{J}$ -integral will establish a link between the critical energy release rate and the critical crack opening displacement. Starting with the relation $\mathit{t}\left(\mathit{\delta }\right)$ and choosing a contour $\mathrm{\Gamma }$ for the evaluation of the $\mathcal{J}$ -integral that surrounds the cohesive zone gives for an effective crack opening $\delta =\mid \mathit{\delta }\mid$ ,

$${\mathcal{G}}_c={\int }_{\mathrm{\Gamma }}\mathit{t}\cdot \mathit{\delta }\mathrm{d\Gamma }={\int }_0^{{\delta }_c}\mathit{t}\mathrm{d}\delta .$$ E6

2.3 Weak form of the problem and finite element discretization

The motions of a solid can be characterized by recourse to Hamilton’s principle of stationary action. The action of a motion within a closed time interval $\left[t_0,t\right]$ defines a functional $\mathcal{I}\left[\mathit{u}\right]={\int }_{t_0}^t\mathcal{L}\mathrm{d}t$ with the Lagrangian function as the difference of kinetic and potential energy, $\mathcal{L}\left(\mathit{u},\mathit{v}\right)=\mathcal{K}\left(\mathit{v}\right)-\mathrm{\Pi }\left(\mathit{u}\right)$ . Stationarity demands the first variation to vanish,

$$\delta \mathcal{I}=\delta {\int }_{t_0}^t\mathcal{L}\left(\mathit{u},\mathit{v}\right)\mathrm{d}t={\int }_{t_0}^t\delta \mathcal{K}\left(\mathit{v}\right)-\delta \mathrm{\Pi }\left(\mathit{u}\right)\mathrm{d}t=0\forall t\ge t_0,$$ E7

and with $\delta \mathrm{\Pi }=\delta W_{\mathrm{int}}-\delta W_{\mathrm{ext}}$ and Eqs. (2)–(4) we obtain

$${\int }_{t_0}^t\left({\int }_{\mathrm{\Omega }}\rho \mathit{v}\cdot \delta \mathit{v}+\delta {\mathrm{\Psi }}^e\left(\mathit{u}\right)-\delta \mathit{u}\cdot \overline{\mathit{b}}\mathrm{d\Omega }-{\int }_{{\mathrm{\Gamma }}_C}\delta 〚\mathit{u}〛\cdot \mathit{t}\mathrm{d\Gamma }-{\int }_{{\mathrm{\Gamma }}_2}\delta \mathit{u}\cdot \overline{\mathit{t}}\mathrm{d\Gamma }\right)\mathrm{d}t=0$$ E8

for all admissible test functions $\delta \mathit{u}\in \mathfrak{U}$ with $\mathfrak{U}=\{\delta \mathit{u}\in H^1\left(\mathrm{\Omega }\right)\mid \delta \mathit{u}=0$ on ${\mathrm{\Gamma }}_1\}$ ; $H^1$ denotes the Hilbert space of weak square-integrable functions and its first derivatives. Eq. (8) must hold for all times $t\ge t_0$ , which leads to the weak momentum balance. For a body with stationary crack it is equivalent to

$${\int }_{\mathrm{\Omega }}\delta \mathit{u}\cdot \rho \mathit{a}+\delta \mathit{\epsilon }:\mathit{\sigma }-\delta \mathit{u}\cdot \overline{\mathit{b}}\mathrm{d\Omega }={\int }_{{\mathrm{\Gamma }}_2}\delta \mathit{u}\cdot \overline{\mathit{t}}\mathrm{d\Gamma }\forall \delta \mathit{u}\in \mathfrak{U}.$$ E9

For discretization the domain $\mathrm{\Omega }$ is subdivided into a finite set of nonoverlapping elements. We make use of a conforming ansatz and approximate the displacement field and its variation with

$$\mathit{u}\approx \sum _{i=1}^{n_k}{N}_i{\mathit{u}}_i=\mathbf{N}\widehat{\mathbf{u}},\delta \mathit{u}\approx \sum _{i=1}^{n_k}{N}_i{\mathit{w}}_i=\mathbf{N}\widehat{\mathbf{w}},$$ E10

where $N_i$ are the piecewise ansatz functions collected in a matrix $\mathbf{N}$ and the vectors $\widehat{\mathbf{u}}$ , $\widehat{\mathbf{w}}$ contain all unknown nodal displacements ${\mathit{u}}_i$ of the $n_k$ nodes and its nodal variations ${\mathit{w}}_i$ , respectively. Plugging Eq. (10) into Eq. (8) and a straightforward calculation gives the global system of equations

$$\overline{\mathbf{M}}\ddot{\widehat{\mathbf{u}}}+\overline{\mathbf{K}}\widehat{\mathbf{u}}=\overline{\mathbf{f}}$$ E11

where $\overline{\mathbf{M}}$ and $\overline{\mathbf{K}}$ denote the mass and stiffness matrix of the finite element discretization, and $\overline{\mathbf{f}}$ is the vector of external forces. The Hookean tensor $\mathbb{C}$ is reformulated in a matrix $\overline{\mathbb{C}}$ , and we specify

$$\begin{array}{l}\overline{\mathbf{M}}=\bigcup _{e=1}^{n_e}\underset{{\mathrm{\Omega }}_e}{\int }\rho {\mathbf{NN}}^T\mathrm{d}{\mathrm{\Omega }}_e,\overline{\mathbf{K}}=\nabla \bigcup _{e=1}^{n_e}\underset{{\mathrm{\Omega }}_e}{\int }\mathit{\text{∇}}\mathbf{N}\overline{\mathbb{C}}{\mathbf{\text{∇N}}}^T\mathrm{d}{\mathrm{\Omega }}_e,\\ \overline{\mathbf{f}}=\bigcup _{e=1}^{n_e}\underset{{\mathrm{\Omega }}_e\cap {\mathrm{\Gamma }}_1}{\int }\mathbf{N}\overline{\mathbf{b}}+\mathbf{N}\overline{\mathbf{t}}\mathrm{d\Omega }+{\int }_{{\mathrm{\Gamma }}_1}\mathbf{N}\overline{\mathbf{t}}\mathrm{d\Gamma }.\end{array}$$

Discretization in time is performed by an implicit Euler method, that is, with time step $\mathrm{\Delta }t=t_{n+1}-t_n$ and, thus, velocity and acceleration are approximated by

$${\dot{\widehat{\mathbf{u}}}}^{n+1}\approx \left({\widehat{\mathbf{u}}}^{n+1}-{\widehat{\mathbf{u}}}^n\right){\left(\mathrm{\Delta }t\right)}^{-1}\text{and}{\ddot{\widehat{\mathbf{u}}}}^{n+1}\approx \left({\widehat{\mathbf{u}}}^{n+1}-2{\widehat{\mathbf{u}}}^n+{\widehat{\mathbf{u}}}^{n-1}\right){\left(\mathrm{\Delta }t\right)}^{-2}.$$ E12

3. The cohesive element technique

The nucleation and the propagation of cracks are efficiently modeled through the cohesive zone model where fracture is assumed to happen along an extended crack tip triggered by tractions on the crack flanks, [23, 24]. A particularly appealing aspect of the cohesive zone model is that it fits naturally in the framework of finite element analysis and leads directly to the cohesive element technique introduced by Needleman, Ortiz, and co-workers [25, 26, 27]. The main idea of this approach is to add cohesive interfaces between the continuum elements that are able to model crack growth, see Figure 2. We employ this classical cohesive element approach combined with an automatic fragmentation and cohesive surface insertion procedure. The method has proven to be reliable and efficient for numerous applications, see among others [2, 28, 29, 30].

Figure 2.

Discretization of a solid with the cohesive element technique (left) and geometry of a 2D and 3D cohesive element, respectively (right). The surfaces ${\mathrm{\Gamma }}_C^{-}$ and ${\mathrm{\Gamma }}_C^{+}$ coincide if the element is closed.

3.1 Fundamentals

In cohesive theories, the displacement jump across a cohesive surface ${\mathrm{\Gamma }}_C$

$$\mathit{\delta }=〚\mathit{u}〛={\mathit{u}}^{+}-{\mathit{u}}^{-}$$ E13

plays the role of a deformation measure while the tractions $\mathit{t}$ furnish the work-conjugate stress measure. They follow a traction separation law—the cohesive law—which models locally the loss of material resistance during cracking. If the cohesive element has attained a critical opening displacement ${\delta }_c$ , no tractions can be transfered and the adjacent continuum elements are de-facto disconnected.

Typically, isotropic and anisotropic materials behave differently in crack opening (mode-I separation) and sliding (mode-II and mode-III separation) and, therefore, normal and tangential components of the displacement jump across the surface ${\mathrm{\Gamma }}_C$ have to be treated differently. Given a vector field $\mathit{u}$ over ${\mathrm{\Gamma }}_C$ , its normal and tangential components are $u_n=\mathit{u}\cdot \mathit{n},$ and $u_t=\mid {\mathit{u}}_{\mathrm{\Gamma }}\mid =\mid \mathit{u}-u_n\mathit{n}\mid =\mid \left(\mathit{I}-\mathit{n}\otimes \mathit{n}\right)\mathit{u}\mid$ , and, the corresponding jump components ${\delta }_n=u_n^{+}-u_n^{-}$ and ${\delta }_t=u_t^{+}-u_t^{-}$ follow from Eq. (13). To further simplify the formulation of mixed-mode cohesive laws, we follow [25] and introduce an effective opening displacement

$$\delta =\mid \mathit{\delta }\mid ={\left({\beta }^2{\delta }_t+{\delta }_n\right)}^{1/2}.$$ E14

Here the parameter $\beta$ assigns different weights to the sliding and normal opening components. This allows us to formulate the traction as a function of the effective opening displacement $\delta$ only.

3.2 Cohesive laws

A cohesive law defines the relation between crack opening displacements d and tractions on the crack flanks $\mathit{t}$ . Generally, a cohesive law can be stated in the general form $\mathit{t}={\partial }_{\mathit{\delta }}\mathcal{G}$ , where $\mathcal{G}$ is the specific fracture energy describing the dissipation in the cohesive zone. It is subject to the restrictions imposed by material frame indifference, material symmetry, and the isotropy of sliding. The most general dependence of $\mathcal{G}$ has the form $\mathcal{G}=\mathcal{G}\left({\delta }_n,\mathit{\delta }\cdot \mathit{u}\right)=\mathcal{G}\left({\delta }_n,{\delta }_t\right)$ and with Eq. (14) follows $\mathcal{G}\left(\delta \right)$ . Thus, a key benefit of the potential structure of the cohesive law is that it reduces the identification of the cohesive law from the three components of $\mathbf{t}$ to a single scalar function. Then an effective traction $\mid \mathit{t}\mid =\widehat{t}={\widehat{t}}_{\mathrm{eff}}$ follows,

$${\widehat{t}}_{\mathrm{eff}}={\left({\widehat{t}}_n+{\beta }^{-2}|{\mathit{t}}_t|\right)}^{1/2}.$$ E15

An appropriate choice of cohesive variable is the maximum attained (effective) crack opening displacement ${\delta }_{\max }$ . Loading of the cohesive surface is then characterized by the conditions $\delta ={\delta }_{\max }$ and $\dot{\delta }\ge 0$ ; all other states correspond to unloading. Figure 3 shows different common loading envelopes, whereas unloading is commonly assumed to be linear in $\delta$ to the origin,

Figure 3.

Typical cohesive laws: (left) linear cohesive envelope (blue) of Eq. (17); concave bilinear envelope (red) with ${\delta }_1={\delta }_c/4,{\sigma }_1={\sigma }_c/4$ ; and (right) convex trilinear envelope (blue), modification with smooth transitions (red), both with ${\delta }_1={\delta }_c/4,{\delta }_2={\delta }_c/2$ , and exponential law (green) of Eq. (19) with ${\delta }_0=3{\delta }_c/20$ .

$$\widehat{t}\left(\delta \right)=\frac{{\widehat{t}}_{\max }}{{\delta }_{\max }}\delta ,\text{if}\delta <{\delta }_{\max }\text{or}\dot{\delta }<0.$$ E16

The simplest cohesive law for brittle materials has a linear loading envelope

$$\widehat{t}\left(\delta \right)={\sigma }_c\left(1-\frac{\delta }{{\delta }_c}\right).$$ E17

The two parameter cohesive strength ${\sigma }_c$ and critical opening displacement ${\delta }_c$ determine via Eq. (6) the specific fracture energy

$${\mathcal{G}}_c=1/2{\sigma }_c{\delta }_c.$$ E18

There are several modifications of Eq. (17), for example bilinear laws for concrete [31] or convex cohesive laws for ductile materials [32]. A cohesive law that can be adapted to brittle and ductile behavior is the universal binding law of Smith and Ferrante [33].

$$\widehat{t}\left(\delta \right)={\sigma }_c\frac{\delta }{{\delta }_0}exp\left(1-\frac{\delta }{{\delta }_0}\right),$$ E19

where ${\delta }_0$ is the position of the maximal traction. Note that here the cohesive stress does not vanish at the critical separation ${\delta }_c$ . The corresponding fracture energy is ${\mathcal{G}}_c={\sigma }_c{\delta }_0\left(exp\left(1\right)-\left(1+{\delta }_c/{\delta }_0\right)exp\left(1-{\delta }_c/{\delta }_0\right)\right)$ . Unloading at $\widehat{t}<{\widehat{t}}_{\max }$ follows a linear relation of the form Eq. (16). Upon closure, the cohesive surfaces are subject to unilateral contact constraints, including friction. Therefore, suitable contact conditions have to be applied under compression. Friction in the cohesive zone is regarded as an independent phenomena, which is not modeled here.

3.3 Finite element implementation

The class of cohesive elements considered here consists of two surface elements, which coincide in the reference configuration of the solid. Each surface element has $n_{\mathrm{coh}}$ nodes; the total number of nodes of the cohesive element is $2n_{\mathrm{coh}}$ . The particular triangulation depicted in Figure 2 is compatible with edge or line elements of four nodes.

Basis of the finite element implementation is Hamilton’s principle given in Eq. (7). Inserting the balance of linear momentum and the static boundary conditions gives the deformation power with variation

$$\delta W_{\mathrm{int}}\left(\mathit{u}\right)={\int }_{\mathrm{\Omega }}\delta \mathit{\epsilon }:\mathbb{C}:\mathit{\epsilon }\mathrm{d\Omega }+\delta W_{\mathrm{int}}^{\mathrm{coh}}\left(\mathit{u}\right)\text{with}\delta W_{\mathrm{int}}^{\mathrm{coh}}\left(\mathit{u}\right)≔{\int }_{{\mathrm{\Gamma }}_{\mathrm{coh}}}\delta 〚\mathit{u}〛\cdot \mathit{t}\mathrm{d\Gamma }$$ E20

where ${\mathrm{\Gamma }}_{\mathrm{coh}}$ are the total cohesive surfaces and $\delta 〚\mathit{u}〛$ is the variation of the separation vector given in Eq. (13). The first term of Eq. (20) ₁ as well as the remaining energy contributions correspond to standard finite element forms and will not be repeated here. The variation of the cohesive energy leads with ansatz Eq. (10) in Eq. (20) ₂ for one cohesive element to

$$\delta W_{int,e}^{\mathrm{coh}}\left(\mathit{u}\right)={\int }_{{\mathrm{\Gamma }}_{\mathrm{coh},e}}\delta {\left(\mathbf{N}{〚\widehat{\mathbf{u}}〛}_e\right)}^T\mathbf{t}\mathrm{d\Gamma }=\delta {{〚\widehat{\mathbf{u}}〛}_e}^T{\int }_{{\mathrm{\Gamma }}_{\mathrm{coh},e}}{\mathbf{N}}^T\mathbf{t}\mathrm{d\Gamma }.$$ E21

The kinetic energy does not have any support in the cohesive element and only the external virtual work has to be determined. For one cohesive element it is

$$\delta W_{ext,e}^{\mathrm{coh}}\left(\mathit{u}\right)=\delta W_{int,e}^{\mathrm{coh}}=\delta {{〚\widehat{\mathbf{u}}〛}_e}^T{\mathbf{f}}_{\mathrm{coh},e}^{\pm }\text{with}{\mathbf{f}}_{\mathrm{coh},e}^{\pm }=\mp {\int }_{{\mathrm{\Gamma }}_{\mathrm{coh},e}}{\mathbf{N}}^T\mathbf{t}\mathrm{d\Gamma }.$$ E22

The tangent stiffness matrix follows by its consistent linearization, with the result ${\mathbf{K}}_{\mathrm{coh},e}={\int }_{{\mathrm{\Gamma }}_{\mathrm{coh},e}}{\mathbf{N}}^T\partial \mathbf{t}/\partial 〚\mathit{u}〛\mathbf{N}\mathrm{d\Gamma }$ . To avoid singularities at $\delta \to 0$ in the derivatives of Eq. (17) some numerical modifications are required, cf. [34]. Finally, the equations added by the cohesive elements to the system Eq. (11) can be formulated as ${\mathbf{K}}_{\mathrm{coh}}\cdot \mathrm{\Delta }〚\widehat{\mathbf{u}}〛=\mathrm{\Delta }{\mathbf{f}}_{\mathrm{coh}}$ with ${\mathbf{K}}_{\mathrm{coh}}={\bigcup }_{e=1}^{n_{\mathrm{coh}}}{\mathbf{K}}_{\mathrm{coh},e},{\mathbf{f}}_{\mathrm{coh}}={\bigcup }_{e=1}^{n_{\mathrm{coh}}}{\mathbf{f}}_{\mathrm{coh},e}$ , where $\mathrm{\Delta }$ symbolizes the incremental solution procedure, which is required for the nonlinear crack-opening problem.

3.4 Adaptive meshing

Since in most problems the expected crack path is not known the decision where the cohesive elements should be inserted has to be made during the simulation. The analysis proceeds incrementally in time. Our decision criterion is based on the effective tensile stress given in Eq. (15), which has to exceed a threshold. This means, in every time step of the calculations, this condition is checked for each internal face. The faces that met the criterion are flagged for subsequent processing. A cohesive element will be inserted at the flagged face and in this manner, the shape and location of a successive crack front is itself an outcome of the calculations.

Within the finite element mesh the insertion of cohesive elements requires topological changes. The local sequential numbering of the corner-nodes defines the orientation; the mid-side node is subsequently duplicated. Owing to the variable environment of the edges in the triangulation, the data structure has to be adapted as illustrated in Figure 4 with case 1: the marked edge with nodes k ₁ and k ₂ is inside of the body, the edge with all nodes will be duplicated: $\overline{k_1{k}_2}\cap \left(\mathrm{\partial \Omega }\cup \mathrm{\Gamma }\left(t\right)\right)=\varnothing$ ; case 2: one of the marked nodes is on the boundary, one node will be duplicated: $k_1\in \mathrm{\partial \Omega }\cup \mathrm{\Gamma }\left(t\right)⊻k_2\in \mathrm{\partial \Omega }\cup \mathrm{\Gamma }\left(t\right)$ ; and case 3: both nodes are on the boundary, two nodes (all nodes) will be duplicated: $k_1\in \mathrm{\partial \Omega }\cup \mathrm{\Gamma }\left(t\right)\wedge k_2\in \mathrm{\partial \Omega }\cup \mathrm{\Gamma }\left(t\right)$ , for the two-dimensional situation. In three dimensions the situation consists of four cases but the main idea is the same. However, loops over all facets, edges, and nodes of the mesh are computationally expensive and require very efficient search and insertion algorithms, cf. [35].

Figure 4.

Classification according to whether the fractured segment has zero (case 1), one (case 2), or two (case 3) nodes on the boundary.

3.5 Effect of the crack initiation criteria in the cohesive model

Here we illustrate the influence of the cohesive strength ${\sigma }_c$ and critical opening displacement ${\delta }_c$ on the cohesive element simulations. Exemplarily we investigate a plane mode-I tension test under quasistatic conditions; however, the observations have been confirmed also in other models.

The test specimen with material data $E=60\mathrm{GPa},\nu =0.2$ has a unit domain $\left[0,1\right]\times \left[0,1\right]$ ; all lengths are given in meter. On the vertical boundaries the displacements in $x$ -direction are constrained and on the lower and upper side $y=0$ and $y=1$ , a traction boundary condition is prescribed, ${\overline{t}}_y=15\mathrm{MPa}$ . A crack is predefined in the area $x\in \left[\mathrm{0,0.2}\right]$ , $y=0.5$ , by including cohesive elements that are completely open. We employ a cohesive law with linear envelope given in Eq. (17), set the cohesive strength ${\sigma }_c=20\mathrm{MPa}$ and the critical crack opening displacement ${\delta }_c={10}^{-4}$ . Following Eq. (18) this corresponds to a critical energy release rate of ${\mathcal{G}}_c=1\mathrm{N}/\mathrm{mm}$ .

During the simulation a cohesive element will be added if the effective traction Eq. (15), here $\widehat{t}={\widehat{t}}_{\mathrm{eff}}\left(\beta \to \infty \right)$ , exceeds a critical value, $\widehat{t}\ge {\sigma }_i$ . We identify the insertion criterium ${\sigma }_i$ with the cohesive strength ${\sigma }_c$ . Because the specimen is pulled with constant load, the computation is stopped at a crack length of 0.7 m.

Figure 5 demonstrates the crack evolution in a mesh of 25 × 25 squares, each divided into eight triangular finite elements with linear shape functions. The computed stresses ${\sigma }_y$ in Figure 5 correspond very well to the analytical solution. They disappear at the crack flanks and are maximal at the crack tip with the typical butterfly shape.

Figure 5.

Crack propagation of the mode-I tension test with $\left(x,y\right)\in \left[0,1\right]\times \left[0,1\right]$ m² and an initial crack of length 0.2 m. the meshes show the initial configuration and the crack in the final configuration magnified by a factor of 25. On the right the corresponding stress component ${\sigma }_y$ [MPa] is plotted.

At next we apply the traction ${\overline{t}}_y=15$ MPa as a linearly increasing load within 200 load steps. A cohesive element will be added if the effective tensile stress exceeds ${\sigma }_i={\sigma }_c$ . The crack length is determined via $a≔\left(max\left\{x\in \left[0,1\right]:{\delta }_n\left(x\right)>{\delta }_c\right\}-0.2\right)$ . Figure 6 shows the crack growth normalized to the crack growth at ${\sigma }_c=5$ MPa. Obviously, the insertion criteria has an influence on the appearance of cohesive elements and so—indirectly—on the crack propagation. In consequence, the insertion stress ${\sigma }_i$ should correspond to the cohesive strength ${\sigma }_c$ to avoid artificial numerical stiffness.

Figure 6.

Normalized crack growth versus the critical cohesive stress.

4. The phase-field fracture approach

The evolving crack in a solid with potential energy of Eq. (5) is represented in the phase-field fracture approach by an additional continuous field $s\left(\mathit{x},t\right)\in \left[0,1\right]$ , $\forall \mathit{x}\in \mathrm{\Omega },t\in \mathrm{ℝ}$ . It has a value of $s=1$ in the intact material and indicates for $s=0$ the cracked zones. Continuity requires a transition zone between both phases. Such a transition zone cannot reflect the sharp boundary of a crack but models a diffuse interface instead. Thus, in a phase-field approach the crack set ${\mathrm{\Gamma }}_C\left(t\right)$ is replaced by a regularizing crack-surface density function $\gamma \left(t\right)$ and the corresponding boundary integral in Eq. (5) is approximated by Eq. (1).

The crack-surface density function can be chosen in different ways. If a second-order phase-field approach is considered (like originally proposed in [36, 37]) it has the form

$$\gamma \left(s,\mathbf{\text{∇}}s\right)=\frac{{\left(1-s\right)}^2}{4\epsilon }+\epsilon {\left|\mathbf{\text{∇}}s\right|}^2.$$ E23

The parameter $\epsilon$ has the unit of a length and is a measure for the width of the diffuse interface zone, see Figure 7. Moreover, this parameter weights the influence of the linear and the gradient term whereby the gradient enforces a regularization of the sharp interface. Another ansatz for the crack-surface density function is the fourth-order form

Figure 7.

Analytic solution of a crack in a second- and a fourth-order phase-field.

$$\gamma \left(s,\mathbf{\text{∇}}s,\mathbf{\text{∇}}\cdot \mathbf{\text{∇}}s\right)=\frac{{\left(1-s\right)}^2}{4\epsilon }+\frac{\epsilon {\left|\mathbf{\text{∇}}s\right|}^2}{2}+\frac{{\epsilon }^3{\left(\mathbf{\text{∇}}\cdot \mathbf{\text{∇}}s\right)}^2}{4}$$ E24

which has been used, for example, in [38, 39]. However, in a finite element discretization ansatz given in Eq. (24) requires ${\mathcal{C}}^1$ -continuous basis functions. In [40] we investigated this formulation in more detail, applied quadratic NURBS-ansatz functions, and found that a smoother crack and a better convergence rate can be achieved.

Here we use the ansatz of Eq. (23) and the corresponding total potential energy reads

$$W\left(\mathit{u},s\right)={\int }_{\mathrm{\Omega }}{\mathrm{\Psi }}^e\left(\mathit{u}\right)+{\mathcal{G}}_c\gamma \left(s\right)\mathrm{d\Omega }={\int }_{\mathrm{\Omega }}\frac{1}{2}\mathit{\epsilon }:{\mathbb{C}}^{\ast }:\mathit{\epsilon }+{\mathcal{G}}_c\left(\frac{{\left(1-s\right)}^2}{4\epsilon }+\epsilon {\left|\mathbf{\text{∇}}s\right|}^2\right)\mathrm{d\Omega }$$ E25

The tensor ${\mathbb{C}}^{\ast }$ in the elastic strain energy density is derived from Eq. (2) by means of a substitute-material approach, ${\mathbb{C}}^{\ast }=g\left(s\right)\mathbb{C}$ , where $g\left(s\right)$ is a degradation function, which is here of the form $g\left(s\right)=s^2+s_0$ , with $s_0\ll 1$ being a small parameter introduced for numerical reasons only. With $s=0$ tensor ${\mathbb{C}}^{\ast }$ models the empty crack.

At this point the evolution equation for the phase-field parameter s is stated in a general Allen-Cahn form,

$$\dot{s}=-\overline{M}{\delta }_s\mathrm{\Psi }.$$ E26

where the nonnegative function $\overline{M}$ is the mobility with unit [ ${\mathrm{mm}}^2/\mathrm{Ns}$ ]; ${\delta }_s$ denotes the variational or total derivative of $\mathrm{\Psi }$ wrt. s. With Eq. (25) it follows

$$\dot{s}=-\overline{M}{\delta }_s\mathrm{\Psi }=-\overline{M}\left[\mathit{\epsilon }:\mathbb{C}:\mathit{\epsilon }s-\frac{{\mathcal{G}}_c}{2\epsilon }\left(1-s+4{\epsilon }^2\mathbf{\text{∇}}\cdot \mathbf{\text{∇}}s\right)\right].$$ E27

This evolution equation can also be deduced in a fully variational manner from energy dissipation, cf. [41]. For more theoretical details we also refer to our recent works [42, 43].

4.1 Elastic strain energy split

The quadratic form of the elastic strain energy density does not distinguish between tensile and pressure states in the material. A direct use of the formulations Eq. (25) or Eq. (27) would allow a crack to grow also in a compressive regime, which clearly contradicts the physics of the underlying problem. For that reason a split of the elastic strain energy into a tensile and a compressive part is necessary. We use the degradation function $g\left(s\right)$ and state ${\mathrm{\Psi }}^e\left(\mathit{\epsilon },s\right)=g\left(s\right){\mathrm{\Psi }}^{e+}+{\mathrm{\Psi }}^{e-}$ . This split is based on the spectral decomposition of the strain tensor $\mathit{\epsilon }={\sum }_{a=1}^3{\epsilon }_a{\mathit{n}}_a\otimes {\mathit{n}}_a$ , where ${\epsilon }_a$ denote the principal strains and ${\mathit{n}}_a$ the corresponding principal directions, $a=\mathrm{1,2,3}$ . Based on this representations and using the Macaulay brackets

$${⟨x⟩}^{+}=\left\{\begin{array}{cc}0,& x<0\\ x,& x\ge 0,\end{array}\right.{⟨x⟩}^{-}=\left\{\begin{array}{cc}x,& x<0\\ 0,& x\ge 0,\end{array}\right.$$

we define the positive and the negative parts of the strain tensor as ${\mathit{\epsilon }}^{\pm }={\sum }_{a=1}^3{⟨{\epsilon }_a⟩}^{\pm }{\mathit{n}}_a\otimes {\mathit{n}}_a$ . The positive parts contain contributions due to positive dilatation and contributions due to positive principal strains. Only this part of the strain energy is responsible for crack growth. A similar decomposition can be deduced for the stress tensor $\mathit{\sigma }$ . This finally leads to an elastic energy density function, which only accounts for tension, ${\mathrm{\Psi }}^{e+}=1/2{⟨\mathit{\sigma }⟩}^{+}:{⟨\mathit{\epsilon }⟩}^{+}=1/2{⟨\mathit{\epsilon }⟩}^{+}:{\mathbb{C}}^{\ast }:{⟨\mathit{\epsilon }⟩}^{+}$ .

Furthermore, the irreversibility of the crack growth has to be considered. This can be done by Dirichlet constraints on the phase-field parameter, that is $\dot{s}=0$ if $s=0$ . The same effect has a product of $s$ with the mobility parameter $M=\overline{M}{\mathcal{G}}_c/2\epsilon$ , where we additionally formulate the evolution Eq. (26) in a dimensionless form, $\dot{s}=-M{\delta }_s\mathrm{\Psi }=-M\left[{\mathrm{\Psi }}^{e+}-\left(1-s\right)-4{\epsilon }^2\mathbf{\text{∇}}\cdot \mathbf{\text{∇}}s\right]$ .

4.2 Discretization

The numerical solution of the phase-field fracture model within the finite element framework leads to a coupled-field problem. The weak formulation of the mechanical field is derived in the usual way with the result given in Eq. (9). The weak formulation of the phase-field equation is set up analogously

$${\int }_{\mathrm{\Omega }}\left(\frac{\dot{s}}{\overline{M}}\mathit{\delta s}+{\partial }_s{\mathrm{\Psi }}^{e+}\mathit{\delta s}+4{\mathcal{G}}_c\epsilon \mathbf{\text{∇}}s\mathbf{\text{∇}}\mathit{\delta s}-\frac{{\mathcal{G}}_c}{2\epsilon }\left(1-s\right)\mathit{\delta s}\right)\mathrm{d\Omega }=0\forall \mathit{\delta s}\in \mathfrak{V}$$ E28

with test functions $\mathit{\delta s}\equiv v$ and $\mathfrak{V}=\left\{v\in H^1\left(\mathrm{\Omega }\right)|\mathbf{\text{∇}}v\cdot \mathit{n}=0\mathrm{on}\mathrm{\partial \Omega },v=0\mathrm{on}{\mathrm{\Gamma }}_1\right\}$ .

The shape and test functions given in Eq. (10) are inserted in Eq. (9) to obtain the discrete mechanical problem of Eq. (11). For the phase-field we approximate

$$s\approx \sum _{i=1}^{n_k}{N}_i{s}_i=\mathbf{N}\widehat{\mathbf{s}},\mathit{\delta s}\approx \sum _{i=1}^{n_k}{N}_i{v}_i=\mathbf{N}\widehat{\mathbf{v}}$$ E29

with ansatz functions $N_1,\dots ,N_{n_k}$ for $n_k$ nodes. Plugging this ansatz into Eq. (28) leads after a straightforward calculation to a finite element system of equations. Additionally we specify the phase-field driving force ${\partial }_s{\mathrm{\Psi }}^e=2s{\mathrm{\Psi }}^{e+}=s{⟨\mathit{\epsilon }\left(\mathit{u}\right)⟩}^{+}:\mathbb{C}:{⟨\mathit{\epsilon }\left(\mathit{u}\right)⟩}^{+}$ .

In short hand notation we get for the mechanical problem the matrix Eq. (11), whereby now the Hookean matrix ${\mathbb{C}}^{\ast }\left(\mathit{u},s\right)$ depends on the phase-field. For Eq. (28) the approximation leads to $\overline{\mathbf{M}}\dot{\mathrm{s}}+\left(\overline{\mathbf{A}}+\overline{\mathbf{C}}+\overline{\mathbf{D}}\right)\mathbf{s}=\overline{\mathbf{f}}$ with $\mathbf{s}=\left[s_i\right]$ and

$$\overline{\mathbf{M}}=\bigcup _{e=1}^{n_e}\frac{1}{M}{\int }_{\mathrm{\Omega }}{\mathbf{NN}}^T\mathrm{d\Omega },\overline{\mathbf{C}}=\bigcup _{e=1}^{n_e}{\int }_{\mathrm{\Omega }}{\mathbf{NN}}^T\mathrm{d\Omega },\overline{\mathbf{f}}=\bigcup _{e=1}^{n_e}{\int }_{\mathrm{\Omega }}{\mathbf{NN}}^T\mathbf{1}\mathrm{d\Omega }$$

$$\overline{\mathbf{A}}=\bigcup _{e=1}^{n_e}4{\epsilon }^2{\int }_{\mathrm{\Omega }}\mathbf{\text{∇N}}{\left(\mathbf{\text{∇N}}\right)}^T\mathrm{d\Omega },\overline{\mathbf{D}}=\bigcup _{e=1}^{n_e}\frac{4\epsilon }{{\mathcal{G}}_c}{\int }_{\mathrm{\Omega }}{\mathrm{\Psi }}^{e+}\left(\mathit{u},s\right){\mathbf{NN}}^T\mathrm{d\Omega }.$$

Here we assumed the same ansatz functions for $\mathit{u}$ and s for the ease of writing.

After the discretization in time by using Eq. (12) the following system of global equations is obtained:

$$\left[\overline{\mathbf{M}}{\left(\mathrm{\Delta }t\right)}^{-1}+\overline{\mathbf{A}}+\overline{\mathbf{C}}+\overline{\mathbf{D}}\right]{\mathbf{s}}^{n+1}=\overline{\mathbf{f}}+{\left(\mathrm{\Delta }t\right)}^{-1}\overline{\mathbf{M}}{\mathbf{s}}^n$$ E30

Note that the coupled problem of Eq. (11) and Eq. (30) is nonlinear. The solution of the implicit problem is obtained with recourse to a Newton-Raphson method. The necessary linearization (tangent stiffness matrix) can be calculated monolithically or by recourse to a staggered scheme. For the HB experiment we employ an explicit time discretization, which simplifies the solution.

4.3 Effect of the coefficients M and $\epsilon$

With the problem formulation at hand we now illustrate the influence of the parameter mobility M and regularization length $\epsilon$ in the phase-field fracture model. To this end we again investigate a plane mode-I tension test of unit length $\left[0,1\right]\times \left[0,1\right]$ , with a predefined crack at $x\in \left[\mathrm{0,0.5}\right]$ , $y=0.5$ , and boundary conditions as depicted in Figure 8. The finite element mesh of $100\times 100\times 4$ triangular elements has a uniform size of $h=0.005$ . The simulation is quasistatic with a prescribed, linearly increasing displacement up to $\overline{u}=\pm 1.2\cdot {10}^{-3}$ .

Figure 8.

Mode-I tension-test with $\left(x,y\right)\in \left[0,1\right]\times \left[0,1\right]$ , an initial crack of length 0.5, and a prescribed total displacement of $\overline{u}=\pm 1.2\cdot {10}^{-3}$ . Plotted is the normalized crack length a versus the kinematic mobility M.

The kinematic mobility parameter $M=\overline{M}{\mathcal{G}}_c/\left(2\epsilon \right)$ is responsible for the rate of phase-field evolution. For higher values of M the crack will show (and propagate) faster. In Figure 8 the crack growth for different values of M normalized with the crack length for $M=1$ is demonstrated. We see a steady state for $M\ge 50 {\mathrm{mm}}^2\mathrm{sN}$ . In this case the crack length is defined as $a≔max\left\{x\in \left[0,1\right]:s\left(x\mathrm{,0.5}\right)\le 0.05\right\}-0.5.$ The mobility M [s] can be seen as the inverse of a kinematic viscosity $\eta$ [s⁻¹], $M=1/\eta$ , which illustrates its effect on the crack evolution within a phase-field fracture computation. Large values of M correspond to a small viscosity, that is, the material will crack fast, whereas a highly viscous material retards crack growth. In this sense, low values of M correspond a high dissipation. However, in the case of quasi-static brittle fracture there exist no local material dissipation because energy is only stored by elastic deformations and in the surface energy of the crack. Consequently, the mobility parameter M has to be sufficiently large for meaningful computations, it is only limited by the numerics.

The influence of the length-scale parameter $\epsilon$ on the crack is now demonstrated. For definiteness, the domain where the phase-field parameter s is larger than 0.95 is called crack-free. Hence the width of the crack is defined as

$$b≔{\left|\{y\in \left[0,1\right]:s(0.9,y)\le 0.95\}\right|}_e,$$ E31

whereby the Lebesgue-measure of a set is denoted with ${\left|\cdot \right|}_e$ and the diffuse interface zone is included in the relevant set.

Please note, that the effect of the length-scale parameter $\epsilon$ is two-fold. On one hand, it enters the material because, in the sense of Griffith’s criteria for crack growth, the fracture energy density ${\mathcal{G}}_c/\epsilon$ competes with the elastic energy density ${\mathrm{\Psi }}^{e+}$ determined by $E,\nu$ or combinations thereof. Here $\epsilon$ has the effect of a material parameter. On the other hand, $\epsilon$ is determined by the mesh size $h$ because it has to be large enough to enable the approximation of a diffuse interface. This clearly requires an adaption on the mesh size, $\epsilon >h$ .

In Figure 9 the crack width b is plotted over different values of $\epsilon$ . The computation starts with $\epsilon =0.01$ _, which is the smallest possible parameter for the chosen mesh size, that is $\epsilon \approx h$ . A clear proportionality between the length-scale parameter $\epsilon$ and the width of the crack can be seen, and we set $b/2=\mathit{c\epsilon }$ with constant $c>1$ . The exact value of c depends on the model via definition of Eq. (31). At any case it is greater than one and so it is obvious that $\epsilon$ is smaller than the diffuse crack width, see also Figure 7. Therefore, to approximate a sharp crack the length-scale parameter $\epsilon$ should be as small as possible. Based on some parameter studies we propose to choose $\epsilon =1\dots 3h$ for linear shape functions, which is similar to other author’s suggestion of $l_c=2\epsilon >2h$ , [12].

Figure 9.

Influence of the length-scale parameter $\epsilon$ on the crack width in the converged state of the mode-I-tension-test of Figure 8 . In the left panel the half crack width $b/2$ is plotted versus $\epsilon$ , the linear regression coefficient is $R^2=0.9986$ , and the right plot demonstrates the phase-field s in the diffuse interface zone to the crack in corresponding colors.

5. Simulation of the HB-spalling experiment

A classical Split-Hopkinson-Pressure Bar consists of a steel projectile (striker), an incident bar and a transmission bar. The specimen is placed between the bars and an analysis of the propagating waves allows to deduce its Young’s modulus. For our UHPC mixture the result is $E=59$ GPa; details of the experiments can be found in [20, 44]. For the spalling experiment the HB setup is modified; we have a striker, an incident bar and a cylindrical specimen but no transmission bar. The impact of the striker generates a compressive stress pulse traveling through the incident bar. When the pulse reaches the end of the incident bar, a part of the stress is transmitted into the specimen. This stress pulse ${\sigma }_T$ is of interest for us. At the free end of the specimen it is reflected as a tensile pulse and spalling occurs if the superposition of both pulses is beyond the cohesive strength, cf. Figure 10. Please note that with the formation of a crack boundary the wave propagation and reflection continue, which may result in several further cracks; this can be seen in the specimen of Figure 1.

Figure 10.

Schematics of the HB-spalling experiment and typical incoming stress wave in the middle of the specimen with $t_1=0$ and $t_2=t_3=30 \mathrm{\mu s}$ .

The aim of spalling experiments is to determine a material’s resistance to fracture, specifically its fracture energy $W_c$ or specific energy ${\mathcal{G}}_c$ and/or its tensile resistance $R_m^t$ . The spalling test presumes brittle materials which can sustain compression but fails under tension. Plastic deformations do not matter, neither do temperature effects; all tests are conducted under ambient conditions.

5.1 Finite element discretization

The UHPC specimen has a length of 200 mm and a diameter of 20 mm. Because of the cylindrical symmetry of the problem we can use an axialsymmetric finite element model. This model maps a fully three-dimensional material behavior with the reduced effort of a plane mesh, which allows us to do extensive parametric studies.

A first challenge was the correct reproduction of the incident and reflected stress pulses in the specimen. From the strain pulse measured in the incident bar, the difference in impedance, and a low-amplitude pulse measured in the specimen we conclude on the shape of the transmitted wave. It is applied on the (left) boundary as a pressure impulse of trapezoidal form $\overline{q}={\sigma }_{\max }f\left(t\right)$ with $f\left(0\le t\le t_1\right)=\mathit{t}/{\mathit{t}}_1$ , $f\left(t_1<t\le t_2\right)=1$ , and, $f\left(t_2<t\le t_3\right)=t_3-t/t_3-t_2$ , where $t_1=12\mathrm{\mu s}$ and $t_2=47\mathrm{\mu s}$ and $t_3=62\mathrm{\mu s}$ and ${\sigma }_{\max }=-17\mathrm{MPa}$ . Within the $80\mathrm{\mu s}$ of simulation time the stress wave travels through the specimen, is reflected and reaches the left end again. Further propagation (and further cracks) are not considered. The average velocity of the specimen before spallation is $v_{\text{spec}}=6.8$ m/s. For time discretization we use a special central difference scheme with a weighted displacement field, which results in stress pulses that largely correspond to the measured data, see Figure 10; for numerical details we refer to [34, 45].

5.2 Fracture parameter

The dynamic tensile resistance of a brittle material $R_m^t$ is usually defined as the maximum tension a material can sustain. A higher stress results in fracture—in our experiments in spallation. There are several attempts but no established way to deduce the tensile resistance directly from the data measured in the experiments, see [46] for a discussion. One way is to measure the incident and reflected waves, “shift” them to the position of fracture z _c, and determine the superposed elastic stress state there. The dynamic tensile strength is then defined as the level of the tensile stress reached at the location of fracture. In this sense we identify $R_m^t$ with the cohesive strength ${\sigma }_c$ and try to capture it by adaptive insertion of cohesive elements.

5.2.1 Cohesive element technique

Our specimen is meshed uniformly with triangular finite elements (2560 elements) and a mixed-mode cohesive law is employed. We start with the linear envelope given in Eq. (17) and an effective opening displacement ${\delta }_c=10\mathrm{\mu }$ m for ${\sigma }_c=15\mathrm{MPa}$ , $G_c=75$ N/m, $\beta =1$ . As outlined in Section 3.5 a cohesive element will be added if the effective tensile stress given in Eq. (15) exceeds the value of ${\sigma }_c$ . In tension, the elements can subsequently open, in compression contact conditions apply. Once the opening has exceeded ${\delta }_c$ , a crack has formed. The simulation stops after a spallation plane has built or after $120 \mathrm{\mu s}$ _.

Figure 11 shows one symmetry half of the specimen at the end of the simulation for different values of ${\sigma }_c$ . At first, with ${\sigma }_c=20\mathrm{MPa}$ , the cohesive stress is obviously too high, no elements are inserted and no crack can grow. For ${\sigma }_c=17\mathrm{MPa}$ a small, localized crack zone develops. Lowering ${\sigma }_c=14\dots 8$ MPa gives a wider zone and, moreover, the cracked zone moves toward the free end. This follows from the fact that the superposed pulse ${\sigma }_R+{\sigma }_I$ reaches the value of ${\sigma }_c$ earlier. However, this position does not correspond to the measured crack position. In all our experiments the crack appeared at $z_c=120\dots 130$ mm, with coordinate $z$ starting at the free end.

Figure 11.

Finite element mesh of the specimen with adaptively inserted cohesive elements. A variation of the critical cohesive stress leads to different positions of the crack. Here the inserted but not yet fully open elements are purple, the cracks (open elements) are black.

Obviously, the cohesive stress ${\sigma }_c$ plays a significant role in the simulation of crack growth with cohesive elements. If it is too low, an unrealistic scattered crack zone will be computed, which does not give fully opened cracks (assuming a fixed fracture energy). However, some energy is dissipated here and the crack plane cannot appear at the right position anymore.

We further studied the influence of the critical crack opening displacement ${\delta }_c$ . Here we observe for higher values of ${\delta }_c$ a longer time of crack opening whereas ${\delta }_c$ has no influence on the position of the crack. Studies with the exponential cohesive law given in Eq. (19) show a very similar behavior in variations of ${\sigma }_c$ and ${\delta }_0$ , cf. Figure 12.

Figure 12.

Influence of the critical crack opening displacement on the time of crack initiation t ₁ and of total crack opening t ₂ for ${\sigma }_c=15$ MPa with (left) ${\delta }_c$ in the linear cohesive law Eq. (17) and (right) ${\delta }_0$ in the exponential form with ${\delta }_c=30\mathrm{\mu m}$ Eq. (19).

5.2.2 Phase-field fracture

In phase-field simulations the fracture energy is the essential parameter for crack growth. Other parameter, like the mobility, are of numerical nature and can be calibrated. Further, the length-scale parameter $\epsilon$ is responsible for the width of the diffuse interface. Because of the fact, that $\epsilon$ depends on the mesh size, the mesh has to be very fine in the domain of a potential crack. To avoid numerical artifacts we chose here a uniform mesh of 5760 triangular elements; the mesh size is $h=5/6$ mm. We set $\epsilon =3h$ and $\tilde{M}=200,000{\mathrm{s}}^{-1}$ and define regions with phase-field parameter $s<0.2$ as cracked.

For a large critical energy release rate, ${\mathcal{G}}_c>150\mathrm{N}/\mathrm{m}$ , no crack is computed. Obviously, here such values are too high. Reducing the value leads to the appearance of a zone with $s<1$ ; however, it depends on the defining threshold if this is already considered to be a crack. For ${\mathcal{G}}_c=90\mathrm{N}/\mathrm{m}$ a sharp crack zone appears at the expected position. Smaller values, for example ${\mathcal{G}}_c=30\mathrm{N}/\mathrm{m}$ in Figure 13(left), show a wider cracked zone but the position of the crack is nearly the same, see Figure 13(right). Specifically we define the diffuse zone as

Figure 13.

Influence of the specific fracture energy ${\mathcal{G}}_c$ on the width given in Eq. (32) of the crack; left plot shows the phase-field parameter at crack position in longitudinal direction and right plot shows width b as function of ${\mathcal{G}}_c$ , approximated with R-coefficient $R^2\approx 0.9894$ .

$$b_{\text{trans}}≔{\left|\{z\in \left[\mathrm{0,200}\right]:0.2\le s(0,z)\le 0.8\}\right|}_e$$ E32

and the cracked zone is given by $s\left(0,z\right)\le 0.8$ .

Figure 14 shows the effect of the specific fracture energy on the time of crack evolution. In opposite to the cohesive model the time difference between crack formation and final state grows linearly with ${\mathcal{G}}_c$ . Similar is the situation for an decrease of parameter $\epsilon$ . A value of $\epsilon <6h$ increases the time of crack formation at constant ${\mathcal{G}}_c$ .

Figure 14.

Influence of the specific fracture energy ${\mathcal{G}}_c$ is shown in left plot and the length parameter $\epsilon$ on the time of crack initiation t ₁ and of total crack opening t ₂ is shown in the right plot, computed with $M=0.2{\mathrm{\mu s}}^{-1}$ and $\epsilon =3h$ (left), $G_c=90\mathrm{N}/\mathrm{m}$ (right).

Further studies have been performed to simulate the experiments, cf. [34], and with the collected knowledge we conclude, that for UHPC the tensile strength $R_m^t$ , the critical opening displacement ${\delta }_c$ _, and the specific fracture energy ${\mathcal{G}}_c$ are in the ranges $12<R_m^t\left[\mathrm{MPa}\right]<18,5<{\delta }_c\left[\mu \mathrm{m}\right]\le 17.5,$ and $40\le {\mathcal{G}}_c\left[\mathrm{N}/\mathrm{m}\right]\le 105$ .

6. Inverse analysis: determination of the ${\mathcal{G}}_c$

The aim of our study was to evaluate the two fracture simulation methods for the use in an inverse analysis where the measured data of the experiment are used to deduce fracture parameters and the obtained results are used to simulate the experiment. The deduced parameters are considered “correct” when the difference between experiment and simulation is small.

6.1 Derivation of the specific fracture energy from spallation

After spallation two fragments result with the crack located at the position where the stress exceeds the tensile resistance first. Depending on the energy of the incoming wave the same process may continue in both fragments with the results of additional cracks. The total fracture energy $W_c$ corresponds to the amount of work necessary to form such a new surface. In order to determine $W_c$ we balance the energy before and after crack initiation, that is, at time t ₁ right before cracking and at time t ₂ immediately after the crack opened. In the simple case of two fragments, cf. Figure 15, with masses $m_{\mathrm{fra},1}$ , $m_{\mathrm{fra},2}$ _, and $m_{\text{spec}}=m_{\mathrm{fra},1}+m_{\mathrm{fra},2}$ it follows for the difference of kinetic energy:

Figure 15.

Cracked specimen in a phase-field simulation and fragments after spallation.

$$\mathrm{\Delta }\mathcal{K}=\frac{1}{2}{m}_{\text{spec}}{\left(v_{\text{spec}}\left(t_{\underset{¯}{1}}\right)\right)}^2-\frac{1}{2}{m}_{\mathrm{fra},1}{\left(v_{\mathrm{fra},1}\left(t_{\underset{¯}{2}}\right)\right)}^2-\frac{1}{2}{m}_{\mathrm{fra},2}{\left(v_{\mathrm{fra},2}\left(t_{\underset{¯}{2}}\right)\right)}^2$$ E33

where $v_{\mathrm{fra},1},v_{\mathrm{fra},2}$ refer to the velocities of the fragments 1 and 2. Because we presume the UHPC to behave linear elastically, the loss of kinetic energy between initial and spalled state will completely be related to the fracturing process. Consequently, we state $\mathrm{\Delta }\mathcal{K}=W_c$ to be the fracture energy of the specimen.

In order to deduce ${\mathcal{G}}_c$ the fracture energy Eq. (33) needs to be referred to the fractured surface $A_c$ of the specimen. Ideally, the crack is smooth and perpendicular to the specimen’s axis. Then the specific fracture energy simply follows as ${\mathcal{G}}_c=W_c/A_c$ with $A_c=2\pi r_{\text{spec}}^2$ . In practice, there may be spalled and rough crack surfaces and so we replaced $A_c$ by the measured surface size.

For the inverse analysis we proceed as follows: we use the data obtained in our previous simulations to define a range of input data ${\mathcal{G}}_c^{\mathrm{inp}}$ , simulate the spallation experiment, and from the computed velocities and masses we deduce via Eq. (33) the “measured” value ${\mathcal{G}}_c^{\mathrm{sim}}$ . To this end we determine the velocities of the fragments by $v_{\text{spec}}={\left(2{\sum }_e{\mathcal{K}}_e/m_{\text{spec}}\right)}^{1/2}$ , $v_{\mathrm{fra},1}={\left(2{\sum }_{e\left(\mathrm{fra},1\right)}\mathcal{K}{e}_{}/m_{\mathrm{fra},1}\right)}^{1/2}$ , and, $v_{\mathrm{fra},2}={\left(2{\sum }_{e\left(\mathrm{fra},2\right)}{\mathcal{K}}_e/m_{\mathrm{fra},2}\right)}^{1/2}$ , with the kinetic energy of the finite elements, ${\mathcal{K}}_e=\frac{1}{2}{m}_e{\left|{\mathit{v}}_e\right|}^2$ .

6.2 Cohesive element technique

We mesh the specimen uniformly with 2560 elements and employ the linear cohesive law from Eq. (17) with parameter $\beta =1$ . The cohesive stress is ${\sigma }_c=15$ MPa to ensure crack formation at the expected position and we vary the critical opening displacement ${\delta }_c$ to obtain values of ${\mathcal{G}}_c$ between 20 N/m and 150 N/m. All computations result in similar fragments with masses $m_{\mathrm{fra},1}\in \left[\mathrm{895,905}\right]\mathrm{g}$ and $m_{\mathrm{fra},2}\in \left[\mathrm{572,582}\right]\mathrm{g}$ . Also the fractured surfaces show little variations, $A_c\in \left[\mathrm{9.5,11.5}\right]$ cm².

In Figure 16 the computed specific fracture energy ${\mathcal{G}}_c^{\mathrm{sim}}$ is shown. In the optimal case is ${\mathcal{G}}_c^{\mathrm{inp}}={\mathcal{G}}_c^{\mathrm{sim}}$ , this is marked with a dotted line. The computed values show a clear proportionality between ${\mathcal{G}}_c^{\mathrm{inp}}$ and its simulated counterpart ${\mathcal{G}}_c^{\mathrm{sim}}$ . This generally validates the chosen method of energy balance. However, the value of ${\mathcal{G}}_c$ deduced from the simulation is higher than the real one. This corresponds to a higher difference in the kinetic energy, that is, for the given initial velocity the velocity of the fragments is underestimated, that is there is too much energy dissipated. This likely results from the fact that some of the adaptively inserted cohesive elements are “useless,” they are inserted in a wrong place and do not fully open. This partial opening costs energy, which is dissipated but does not contribute to spallation. Therefore we observe ${\mathcal{G}}_c^{\mathrm{inp}}<{\mathcal{G}}_c^{\mathrm{sim}}$ with discrepancies between 30 and 60%, see Figure 16.

Figure 16.

Determination of the specific fracture energy with the cohesive element technique: displayed are the values ${\mathcal{G}}_c^{\mathrm{sim}}$ derived from the velocity and mass data obtained in the simulation versus the input parameter ${\mathcal{G}}_c^{\mathrm{inp}}$ of the experiment. The dotted line marks the identity ${\mathcal{G}}_c^{\mathrm{sim}}={\mathcal{G}}_c^{\mathrm{inp}}$ ; the dashed line is approximated with $R^2\approx 0.9928$ and corresponds roughly to ${\mathcal{G}}_c^{\mathrm{sim}}=1.5{\mathcal{G}}_c^{\mathrm{inp}}$ .

6.3 Phase-field fracture

For phase-field simulations we use a finer mesh of 5760 elements, $h=5/6 \mathrm{mm}$ and set $\epsilon =2.5 \mathrm{mm}$ and $M=0.6 {\mathrm{\mu s}}^{-1}$ . At first we accelerate the specimen to $v_{\text{spec}}=6.8\mathrm{m}/\mathrm{s}$ , then the actual simulation starts. Again, during spallation the velocity of the fragments develops differently, the first fragment moves faster then the second one. The time of fracture depends on the input value ${\mathcal{G}}_c^{\mathrm{inp}}$ , the lower it is the easier both fragments split, cf. Figure 14.

Since the phase-field model is a diffuse interface approach, which does give a discrete distinction between the both states, it is necessary to specify tolerances for the states $s\left(t_1\right)$ and $s\left(t_2\right)$ . Here we consider $s\left(t_1\right)$ and set for symmetry $s\left(t_2\right)=1-s\left(t_1\right)$ . A lower “crack initiation” point $s\left(t_1\right)$ clearly needs less time and energy to open the crack. On the one hand, if crack initiation equals crack opening, $s\left(t_1\right)=s\left(t_2\right)=0.5$ , that is the crack appears immediately after initiation and the dissipated energy ${\mathcal{G}}_c^{\mathrm{sim}}$ will be close to zero. The input value ${\mathcal{G}}_c^{\mathrm{inp}}$ would work only as a threshold for the fragments to loose contact and (almost) no dissipation will take place. On the other hand, if we raise $s\left(t_1\right)>0.95$ , $s\left(t_2\right)>0.05$ the crack may not fully open, which gives inaccurate fragment velocities. A stable state we found for $s\left(t_1\right)\in \left[\mathrm{0.75,0.9}\right]$ and so we set $s\left(t_1\right)=0.8$ and $s\left(t_2\right)=0.2$ .

In Figure 17 the computed specific fracture energy ${\mathcal{G}}_c^{\mathrm{sim}}$ over the input value ${\mathcal{G}}_c^{\mathrm{inp}}$ is shown. Again, the optimal case ${\mathcal{G}}_c^{\mathrm{inp}}={\mathcal{G}}_c^{\mathrm{sim}}$ is marked with a dotted line and the computed values show a clear proportionality between ${\mathcal{G}}_c^{\mathrm{inp}}$ and ${\mathcal{G}}_c^{\mathrm{sim}}$ . The value of ${\mathcal{G}}_c^{\mathrm{sim}}$ is now lower than the input value ${\mathcal{G}}_c^{\mathrm{inp}}$ with relative difference less than 10%. Reason for this is basically the definition of crack initiation with $s<1$ , which attributes the beginning of the crack still to the full specimen. Please note that for very small values of ${\mathcal{G}}_c^{\mathrm{inp}}<20$ N/m the specimen cracks basically without dissipation whereas too high values of ${\mathcal{G}}_c^{\mathrm{inp}}$ prevent any crack. In the range of interest, ${\mathcal{G}}_c^{\mathrm{inp}}=30\dots 100$ N/m both values ${\mathcal{G}}_c^{\mathrm{sim}}$ and ${\mathcal{G}}_c^{\mathrm{inp}}$ correlate very well with an R-coefficient of 0.9969. Specifically we have ${\mathcal{G}}_c^{\mathrm{sim}}\approx 0.9{\mathcal{G}}_c^{\mathrm{inp}}$ . Therefore, it is possible and useful to employ a phase-field fracture simulation for a quantitative analysis of the experiments.

Figure 17.

Determination of the specific fracture energy with the phase-field fracture approach: displayed are the values ${\mathcal{G}}_c^{\mathrm{sim}}$ derived from the velocity and mass data obtained simulations with $\epsilon =5/2 \mathrm{mm}$ , $s\left(t_1\right)=0.8$ , $s\left(t_2\right)=0.2$ versus the input parameter ${\mathcal{G}}_c^{\mathrm{inp}}$ of the experiment. The dotted line marks the identity ${\mathcal{G}}_c^{\mathrm{sim}}={\mathcal{G}}_c^{\mathrm{inp}}$ ; the dashed line is approximated with $R^2\approx 0.9969$ and corresponds well to ${\mathcal{G}}_c^{\mathrm{sim}}=0.9{\mathcal{G}}_c^{\mathrm{inp}}$ .

7. Conclusions

In the previous we compared the possibilities of a sharp interface method and a diffuse interface method for crack nucleation and quantitative dynamic fracture analysis. Exemplarily, we validated investigations on the fracture toughness of high-performance concrete in a Hopkinson bar spallation experiment whereby, in particular, the fracture energy values have been determined. Both methods, the cohesive element technique and the phase-field fracture approach, allow numerical simulations of crack growth with an a priori unknown path, and both methods allow to determine the related material parameter in a quantitative manner. Reliability, precision, and numerical costs differ however. Pros and cons of both methods are summarized in the following.

7.1 Model parameters

The core of the cohesive zone model is a cohesive law, $\mathit{t}\left(\mathit{\delta }\right)$ , which describes the forces between the crack flanks as a function of separation. Such cohesive laws allow for pure mode-I cracks in the sense of Griffith as well as for mixed-mode cracks, for example by using the effective traction and separation. Essential cohesive parameters are the critical cohesive stress ${\sigma }_c$ , the critical separation ${\delta }_c$ , and the weight $\beta$ _, which relates shear and tension. These parameters depend on the specific material and can be determined experimentally whereby ${\delta }_c$ is implicitly given via the specific crack energy ${\mathcal{G}}_c$ . Further specifications of the cohesive law may require additional material parameters, for example, in the classical exponential Rose-Ferrante law an additional parameter ${\delta }_0$ needs to be set. All these parameters have a clear physical meaning.

Sensitive for the cohesive element technique is the critical traction for adaptive insertion of the cohesive elements, which has no direct physical background but strongly influences energy dissipation and numerical efficiency. Wrongly inserted elements may dissipate energy but do not contribute to fracture and skew the simulation results.

The phase-field approach to fracture is based on an evolution equation that essentially refers to the elastic strain energy density ${\mathrm{\Psi }}^e\left(\mathit{u},s\right)$ of the material. The remaining relevant parameters are the mobility $M$ , the specific crack energy ${\mathcal{G}}_c$ _, and the length-scale parameter $\epsilon$ . The mobility has only numerical character and controls the rate of phase-field decrease during crack formation. The critical length-scale parameter $\epsilon$ is a measure for the width of the diffuse interface, which relates to both, the finite element mesh size and the material properties. The latter enters the crack growth criteria via the term ${\mathcal{G}}_c/\epsilon$ in the crack resistance and therefore $\epsilon$ needs to be calibrated carefully. The only material parameter in the phase-field model with a clear physical background is the Griffith energy ${\mathcal{G}}_c$ .

7.2 Numerical implementation and computation

In the cohesive zone model cohesive elements are adaptively inserted between the continuum elements to describe the crack opening. The continuum elements themselves are not directly affected and the crack can only propagate along the element boundaries, which results in a certain mesh dependence. The adaptive insertion of cohesive elements require a continuous update of the data structure, which leads to a significant programming effort and also increases the costs of computation. Additionally, the cohesive zone has to be equipped with contact constraints in order to prevent penetration in case of unloading. In total, the numerical implementation of an adaptive cohesive zone model becomes very complex.

In contrast, the structure of the finite element mesh in the phase-field approach remains constant during the simulation. The phase-field parameter can decrease to zero at each node and the crack is able to propagate theoretically everywhere in the whole domain. Essential requirement is a very fine mesh with $\epsilon >h$ .

For numerical computation of the coupled fields $\mathit{u},s$ there exist two different types of solution. On the one hand it is possible to determine the displacements $\mathit{u}$ and the phase-field parameter s successive such that two system of equations have to be built up—this method is called staggered scheme; on the other hand a fully coupled system of equations can be constructed. In general the implementation effort is much larger for the coupled approach and the staggered scheme converges (when implemented correctly) to the same solution so that this method is usually preferred. Much more relevant than the rise in the degrees of freedom by the additional field is, however, the required fine mesh size. The resolution of the mesh needs to be $h<\epsilon$ everywhere where the crack will be able to form. This can be solved by adaptive mesh refinement—which results in a similar effort like the insertion of cohesive elements—or by an a priori fine mesh in the relevant zones, which, in turn, raises the numerical costs tremendously.

7.3 Constraints and driving forces

The major advantage of the cohesive zone model is that the crack properties can be mapped exactly. The local opening is known, the crack width is the separation $\delta$ of the crack flanks, and the corresponding normal and tangential forces follow from the cohesive law. Since the loading and unloading processes are distinguished, compressive forces do not contribute to crack growth. The irreversibility of crack propagation is guaranteed and although the separation can decrease during the simulation the crack will not “heal.” Regarding mixed-mode problems it is positive that the normal and tangential traction components are weighted by a parameter so that different crack opening modes can be realized for each specific material under consideration.

In phase-field fracture by definition a continuous function s and a diffuse interface with width $\epsilon \propto \delta$ models the crack. The values $0<s<1$ can not be interpreted physically so that threshold values have to be chosen. Crack closure needs to be prevented by explicit local constraints, for example $\dot{s}<0$ , or by a modified mobility. In order to ensure that only tensile stresses contribute to crack growth the strain energy function needs to be decomposed in tensile and compressive components, which is by no means unique, cf. [47], and also influences the resulting crack.

Summarizing we state that both methods are mechanically consistent and have a clear variational structure. The cohesive element technique is difficult to implement but provides a strong physical background. For static computations with expected way of crack propagation it is definitely preferred because it allows cohesive laws, which may consider anisotropy, friction, and other material specific properties. In general dynamic applications of unknown crack path, however, its numerical drawbacks, together with the fact that a suboptimal insertion may lead to wrong predictions, dominate. Here the phase-field approach to fracture is clearly the better choice. Unknown crack paths can simply be followed—as long as the mesh resolution is fine enough. The major drawback of phase-field fracture is its parameter sensitivity. Also, extensions to more complex fracture models, which account, for example, for sliding, anisotropy, and interlocking, contradict the original variational derivation and are still an open problem.