## 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.

### Keywords

- Split-Hopkinson bar experiment
- UHPC
- cohesive elements
- phase-field fracture
- inverse analysis
- dynamic fracture
- crack propagation
- crack tracking algorithms

## 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

An alternative way to describe moving boundaries are diffuse interface models where the cracks are smeared over a small but finite length *s* and its gradient. This crack-surface density function allows an approximation of the moving crack boundaries over the body’s domain,

By *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 (

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.

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

### 2.1 Elasto-dynamics

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

where the Lamé material parameters

where

with normal vector _{3} 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

where

The specific energy *W* dissipated per unit newly created surface area

Another fracture criterion is the crack tip opening displacement with critical value

### 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

and with

for all admissible test functions

For discretization the domain

where

where

Discretization in time is performed by an implicit Euler method, that is, with time step

## 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].

### 3.1 Fundamentals

In cohesive theories, the displacement jump across a cohesive surface

plays the role of a deformation measure while the tractions

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

Here the parameter

### 3.2 Cohesive laws

A cohesive law defines the relation between crack opening displacements d and tractions on the crack flanks

An appropriate choice of cohesive variable is the maximum attained (effective) crack opening displacement

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

The two parameter cohesive strength

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].

where

### 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

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

where _{1} 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) _{2} for one cohesive element to

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

The tangent stiffness matrix follows by its consistent linearization, with the result

### 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* _{1} and *k* _{2} is inside of the body, the edge with all nodes will be duplicated:

### 3.5 Effect of the crack initiation criteria in the cohesive model

Here we illustrate the influence of the cohesive strength

The test specimen with material data

During the simulation a cohesive element will be added if the effective traction Eq. (15), here

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

At next we apply the traction

## 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 *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

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

The parameter

which has been used, for example, in [38, 39]. However, in a finite element discretization ansatz given in Eq. (24) requires

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

The tensor *empty* crack.

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

where the nonnegative function *s*. With Eq. (25) it follows

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

we define the positive and the negative parts of the strain tensor as

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

### 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

with test functions

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

with ansatz functions

In short hand notation we get for the mechanical problem the matrix Eq. (11), whereby now the Hookean matrix

Here we assumed the same ansatz functions for *s* for the ease of writing.

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

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 ε

With the problem formulation at hand we now illustrate the influence of the parameter mobility *M* and regularization length

The kinematic mobility parameter *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* [s] can be seen as the inverse of a kinematic viscosity ^{−1}], *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 *s* is larger than 0.95 is called crack-free. Hence the width of the crack is defined as

whereby the Lebesgue-measure of a set is denoted with

Please note, that the effect of the length-scale parameter

In Figure 9 the crack width *b* is plotted over different values of _{,} which is the smallest possible parameter for the chosen mesh size, that is *c* depends on the model via definition of Eq. (31). At any case it is greater than one and so it is obvious that

## 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 *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

The aim of spalling experiments is to determine a material’s resistance to fracture, specifically its fracture energy

### 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

### 5.2 Fracture parameter

The dynamic tensile resistance of a brittle material *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

#### 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 _{.}

Figure 11 shows one symmetry half of the specimen at the end of the simulation for different values of

Obviously, the cohesive stress

We further studied the influence of the critical crack opening displacement

#### 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

For a large critical energy release rate,

and the cracked zone is given by

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

Further studies have been performed to simulate the experiments, cf. [34], and with the collected knowledge we conclude, that for UHPC the tensile strength _{,} and the specific fracture energy

## 6. Inverse analysis: determination of the 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 *t* _{1} right before cracking and at time *t* _{2} immediately after the crack opened. In the simple case of two fragments, cf. Figure 15, with masses _{,} and

where

In order to deduce

For the inverse analysis we proceed as follows: we use the data obtained in our previous simulations to define a range of input data

### 6.2 Cohesive element technique

We mesh the specimen uniformly with 2560 elements and employ the linear cohesive law from Eq. (17) with parameter ^{2}.

In Figure 16 the computed specific fracture energy

### 6.3 Phase-field fracture

For phase-field simulations we use a finer mesh of 5760 elements,

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

In Figure 17 the computed specific fracture energy

## 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, _{,} which relates shear and tension. These parameters depend on the specific material and can be determined experimentally whereby

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 _{,} and the length-scale parameter

### 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

For numerical computation of the coupled fields *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

### 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

In phase-field fracture by definition a continuous function *s* and a diffuse interface with width

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.