Open access peer-reviewed chapter

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

By Tim Dally, Carola Bilgen, Marek Werner and Kerstin Weinberg

Submitted: February 13th 2020Reviewed: March 18th 2020Published: May 6th 2020

DOI: 10.5772/intechopen.92180

Downloaded: 130


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.


  • 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 ΓCt∂Ωin a solid of domain Ωundergoing a deformation χxt:Ω×0ttotal3in a time ttotal. 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 ε. Here an additional field sxt:Ω×0ttotalcharacterizes 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 γ, 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,


By γbeing not only a function of field sxtbut also its gradient 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 Ω3with external boundary Γ∂Ω. The body’s displacement field at point xand time tis denoted by uxt; its velocity and acceleration fields are v=u̇and a=u¨. A crack splits the body into subbodies Ω+Ω=Ωand induces a displacement jump on ΓCas u=u+u. The displacements satisfy the Dirichlet boundary conditions u=u¯at Γ1. The body is loaded with traction t¯at boundary Γ2; it holds Γ1Γ2Γ\ΓC.

2.1 Elasto-dynamics

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


where the Lamé material parameters λ=/1+ν12νand μ=E/21+νare formulated with Young’s modulus Eand Poisson number ν. Equivalently we use the Hookean material tensor Cwith components Cijkl=λδijδkl+μδikδjl+δilδjk. The strain tensor εdescribes small deformations εu=1/2u+uT; the mechanical stress tensor follows as σ=εΨe=Cε. It holds the balance of linear momentum,


where ρis the mass density and b¯a prescribed body force density. Boundary conditions are prescribed as


with normal vector nand traction t; denotes a jump. Eq. (4) 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 ΓCt. 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


where Gcis commonly known as Griffith’s critical energy release rate (Griffith energy). An optimum of Eq. (5) corresponds to crack growth.

The specific energy Gcorresponds to the energy W dissipated per unit newly created surface area A, G=dW/dA. 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, GGc, which is related to the classical stress intensity factors of linear-elastic fracture mechanics given by G=KI2+KII2/E+KIII2/2μwith E=Eregarding the plane stress state and otherwise E=E/1ν2. Thus, Griffith’s criterion can also be written with the critical stress intensity factor Kcas Gc=Kc2/E. The single summands of Gmay also be understood as mode-dependent Griffith energy densities GI, GII, and GIII. The indices denote the fracture modes of crack opening, in-plane and out-of-plane sliding. For a brittle elastic material the J-integral is equal to the strain energy release rate, that is, it holds the identity Jc=Gc.

Another fracture criterion is the crack tip opening displacement with critical value δ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 tand crack opening δ, a standard application of the J-integral will establish a link between the critical energy release rate and the critical crack opening displacement. Starting with the relation tδand choosing a contour Γfor the evaluation of the J-integral that surrounds the cohesive zone gives for an effective crack opening δ=δ,


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 t0tdefines a functional Iu=t0tLdtwith the Lagrangian function as the difference of kinetic and potential energy, Luv=KvΠu. Stationarity demands the first variation to vanish,


and with δΠ=δWintδWextand Eqs. (2)(4) we obtain


for all admissible test functions δuUwith U={δuH1Ωδu=0on Γ1}; H1denotes the Hilbert space of weak square-integrable functions and its first derivatives. Eq. (8) must hold for all times tt0, which leads to the weak momentum balance. For a body with stationary crack it is equivalent to


For discretization the domain Ω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


where Niare the piecewise ansatz functions collected in a matrix Nand the vectors û, ŵcontain all unknown nodal displacements uiof the nknodes and its nodal variations wi, respectively. Plugging Eq. (10) into Eq. (8) and a straightforward calculation gives the global system of equations


where M¯and K¯denote the mass and stiffness matrix of the finite element discretization, and f¯is the vector of external forces. The Hookean tensor Cis reformulated in a matrix C¯, and we specify


Discretization in time is performed by an implicit Euler method, that is, with time step Δt=tn+1tnand, thus, velocity and acceleration are approximated by


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 Γ C − and Γ C + coincide if the element is closed.

3.1 Fundamentals

In cohesive theories, the displacement jump across a cohesive surface ΓC


plays the role of a deformation measure while the tractions tfurnish 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 δ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 ΓChave to be treated differently. Given a vector field uover ΓC, its normal and tangential components are un=un,and ut=uΓ=uunn=Innu, and, the corresponding jump components δn=un+unand δt=ut+utfollow from Eq. (13). To further simplify the formulation of mixed-mode cohesive laws, we follow [25] and introduce an effective opening displacement


Here the parameter β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 δonly.

3.2 Cohesive laws

A cohesive law defines the relation between crack opening displacements d and tractions on the crack flanks t. Generally, a cohesive law can be stated in the general form t=δG, where Gis 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 Ghas the form G=Gδnδu=Gδnδtand with Eq. (14) follows Gδ. 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 tto a single scalar function. Then an effective traction t=t̂=t̂efffollows,


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

Figure 3.

Typical cohesive laws: (left) linear cohesive envelope (blue) of Eq. (17); concave bilinear envelope (red) with δ 1 = δ c / 4 , σ 1 = σ c / 4 ; and (right) convex trilinear envelope (blue), modification with smooth transitions (red), both with δ 1 = δ c / 4 , δ 2 = δ c / 2 , and exponential law (green) of Eq. (19) with δ 0 = 3 δ c / 20 .


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


The two parameter cohesive strength σcand critical opening displacement δcdetermine via Eq. (6) the specific fracture energy


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 δ0is the position of the maximal traction. Note that here the cohesive stress does not vanish at the critical separation δc. The corresponding fracture energy is Gc=σcδ0exp11+δc/δ0exp1δc/δ0. Unloading at t̂<t̂maxfollows 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 ncohnodes; the total number of nodes of the cohesive element is 2ncoh. 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


where Γcohare the total cohesive surfaces and δuis the variation of the separation vector given in Eq. (13). The first term of Eq. (20) 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 Kcoh,e=Γcoh,eNTt/uN. To avoid singularities at δ0in 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 KcohΔû=Δfcohwith Kcoh=e=1ncohKcoh,e,fcoh=e=1ncohfcoh,e, where Δ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 1 and k 2 is inside of the body, the edge with all nodes will be duplicated: k1k2¯∂ΩΓt=; case 2: one of the marked nodes is on the boundary, one node will be duplicated: k1∂ΩΓtk2∂ΩΓt; and case 3: both nodes are on the boundary, two nodes (all nodes) will be duplicated: k1∂ΩΓtk2∂ΩΓt, 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 σcand critical opening displacement δcon 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=60GPa,ν=0.2has a unit domain 01×01; 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=0and y=1, a traction boundary condition is prescribed, t¯y=15MPa. A crack is predefined in the area x0,0.2, 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 σc=20MPaand the critical crack opening displacement δc=104. Following Eq. (18) this corresponds to a critical energy release rate of Gc=1N/mm.

During the simulation a cohesive element will be added if the effective traction Eq. (15), here t̂=t̂effβ, exceeds a critical value, t̂σi. We identify the insertion criterium σiwith the cohesive strength σ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 σyin 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 x y ∈ 0 1 × 0 1 m2 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 σ y [MPa] is plotted.

At next we apply the traction t¯y=15MPa as a linearly increasing load within 200 load steps. A cohesive element will be added if the effective tensile stress exceeds σi=σc. The crack length is determined via amaxx01:δnx>δc0.2. Figure 6 shows the crack growth normalized to the crack growth at σc=5MPa. 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 σishould correspond to the cohesive strength σcto 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 sxt01, xΩ,t. It has a value of s=1in the intact material and indicates for s=0the 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 ΓCtis replaced by a regularizing crack-surface density function γtand 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


The parameter ε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.


which has been used, for example, in [38, 39]. However, in a finite element discretization ansatz given in Eq. (24) requires C1-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


The tensor Cin the elastic strain energy density is derived from Eq. (2) by means of a substitute-material approach, C=gsC, where gsis a degradation function, which is here of the form gs=s2+s0, with s01being a small parameter introduced for numerical reasons only. With s=0tensor Cmodels the 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 M¯is the mobility with unit [mm2/Ns]; δsdenotes the variational or total derivative of Ψwrt. 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 gsand state Ψeεs=gsΨe++Ψe. This split is based on the spectral decomposition of the strain tensor ε=a=13εanana, where εadenote the principal strains and nathe corresponding principal directions, a=1,2,3. Based on this representations and using the Macaulay brackets


we define the positive and the negative parts of the strain tensor as ε±=a=13εa±nana. 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 σ. This finally leads to an elastic energy density function, which only accounts for tension, Ψe+=1/2σ+:ε+=1/2ε+:C:ε+.

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 ṡ=0if s=0. The same effect has a product of swith the mobility parameter M=M¯Gc/2ε, where we additionally formulate the evolution Eq. (26) in a dimensionless form, ṡ=MδsΨ=MΨe+1s4ε2s.

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 δsvand V=vH1Ωvn=0on∂Ωv=0onΓ1.

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 N1,,Nnkfor nknodes. 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 sΨe=2sΨe+=sεu+:C:εu+.

In short hand notation we get for the mechanical problem the matrix Eq. (11), whereby now the Hookean matrix Cusdepends on the phase-field. For Eq. (28) the approximation leads to M¯ṡ+A¯+C¯+D¯s=f¯with s=siand


Here we assumed the same ansatz functions for uand 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 εin the phase-field fracture model. To this end we again investigate a plane mode-I tension test of unit length 01×01, with a predefined crack at x0,0.5, y=0.5, and boundary conditions as depicted in Figure 8. The finite element mesh of 100×100×4triangular elements has a uniform size of h=0.005. The simulation is quasistatic with a prescribed, linearly increasing displacement up to u¯=±1.2103.

Figure 8.

Mode-I tension-test with x y ∈ 0 1 × 0 1 , an initial crack of length 0.5, and a prescribed total displacement of u ¯ = ± 1.2 ⋅ 10 − 3 . Plotted is the normalized crack length a versus the kinematic mobility M.

The kinematic mobility parameter M=M¯Gc/2ε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=1is demonstrated. We see a steady state for M50 mm2sN. In this case the crack length is defined as amaxx01:sx, mobility M [s] can be seen as the inverse of a kinematic viscosity η[s−1], M=1/η, 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 ε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


whereby the Lebesgue-measure of a set is denoted with eand the diffuse interface zone is included in the relevant set.

Please note, that the effect of the length-scale parameter ε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 Gc/εcompetes with the elastic energy density Ψe+determined by E,νor combinations thereof. Here εhas the effect of a material parameter. On the other hand, εis determined by the mesh size hbecause it has to be large enough to enable the approximation of a diffuse interface. This clearly requires an adaption on the mesh size, ε>h.

In Figure 9 the crack width b is plotted over different values of ε. The computation starts with ε=0.01, which is the smallest possible parameter for the chosen mesh size, that is εh. A clear proportionality between the length-scale parameter εand the width of the crack can be seen, and we set b/2=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 εis smaller than the diffuse crack width, see also Figure 7. Therefore, to approximate a sharp crack the length-scale parameter εshould be as small as possible. Based on some parameter studies we propose to choose ε=13hfor linear shape functions, which is similar to other author’s suggestion of lc=2ε>2h, [12].

Figure 9.

Influence of the length-scale parameter ε 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 ε , 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=59GPa; 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 σTis 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   μs .

The aim of spalling experiments is to determine a material’s resistance to fracture, specifically its fracture energy Wcor specific energy Gcand/or its tensile resistance Rmt. 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 q¯=σmaxftwith f0tt1=t/t1, ft1<tt2=1, and, ft2<tt3=t3t/t3t2, where t1=12μsand t2=47μsand t3=62μsand σmax=17MPa. Within the 80μsof 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 vspec=6.8m/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 Rmtis 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 Rmtwith the cohesive strength σcand 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 δc=10μm for σc=15MPa, Gc=75N/m, β=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 σc. In tension, the elements can subsequently open, in compression contact conditions apply. Once the opening has exceeded δc, a crack has formed. The simulation stops after a spallation plane has built or after 120 μs.

Figure 11 shows one symmetry half of the specimen at the end of the simulation for different values of σc. At first, with σc=20MPa, the cohesive stress is obviously too high, no elements are inserted and no crack can grow. For σc=17MPaa small, localized crack zone develops. Lowering σc=148MPa gives a wider zone and, moreover, the cracked zone moves toward the free end. This follows from the fact that the superposed pulse σR+σIreaches the value of σcearlier. However, this position does not correspond to the measured crack position. In all our experiments the crack appeared at zc=120130mm, with coordinate zstarting 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 σcplays 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 δc. Here we observe for higher values of δca longer time of crack opening whereas δchas 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 σcand δ0, cf. Figure 12.

Figure 12.

Influence of the critical crack opening displacement on the time of crack initiation t 1 and of total crack opening t 2 for σ c = 15 MPa with (left) δ c in the linear cohesive law Eq. (17) and (right) δ 0 in the exponential form with δ c = 30 μ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 εis responsible for the width of the diffuse interface. Because of the fact, that ε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/6mm. We set ε=3hand M˜=200,000s1and define regions with phase-field parameter s<0.2as cracked.

For a large critical energy release rate, Gc>150N/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 Gc=90N/ma sharp crack zone appears at the expected position. Smaller values, for example Gc=30N/min 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 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 G c , approximated with R-coefficient R 2 ≈ 0.9894 .


and the cracked zone is given by s0z0.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 Gc. Similar is the situation for an decrease of parameter ε. A value of ε<6hincreases the time of crack formation at constant Gc.

Figure 14.

Influence of the specific fracture energy G c is shown in left plot and the length parameter ε on the time of crack initiation t 1 and of total crack opening t 2 is shown in the right plot, computed with M = 0.2 μs − 1 and ε = 3 h (left), G c = 90 N / 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 Rmt, the critical opening displacement δc, and the specific fracture energy Gcare in the ranges 12<RmtMPa<18,5<δcμm17.5,and 40GcN/m105.

6. Inverse analysis: determination of the Gc

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 Wccorresponds to the amount of work necessary to form such a new surface. In order to determine Wcwe balance the energy before and after crack initiation, that is, at time 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 mfra,1, mfra,2, and mspec=mfra,1+mfra,2it follows for the difference of kinetic energy:

Figure 15.

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


where vfra,1,vfra,2refer 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 ΔK=Wcto be the fracture energy of the specimen.

In order to deduce Gcthe fracture energy Eq. (33) needs to be referred to the fractured surface Acof the specimen. Ideally, the crack is smooth and perpendicular to the specimen’s axis. Then the specific fracture energy simply follows as Gc=Wc/Acwith Ac=2πrspec2. In practice, there may be spalled and rough crack surfaces and so we replaced Acby 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 Gcinp, simulate the spallation experiment, and from the computed velocities and masses we deduce via Eq. (33) the “measured” value Gcsim. To this end we determine the velocities of the fragments by vspec=2eKe/mspec1/2, vfra,1=2efra1Ke/mfra,11/2, and, vfra,2=2efra2Ke/mfra,21/2, with the kinetic energy of the finite elements, Ke=12meve2.

6.2 Cohesive element technique

We mesh the specimen uniformly with 2560 elements and employ the linear cohesive law from Eq. (17) with parameter β=1. The cohesive stress is σc=15MPa to ensure crack formation at the expected position and we vary the critical opening displacement δcto obtain values of Gcbetween 20 N/m and 150 N/m. All computations result in similar fragments with masses mfra,1895,905gand mfra,2572,582g. Also the fractured surfaces show little variations, Ac9.5,11.5cm2.

In Figure 16 the computed specific fracture energy Gcsimis shown. In the optimal case is Gcinp=Gcsim, this is marked with a dotted line. The computed values show a clear proportionality between Gcinpand its simulated counterpart Gcsim. This generally validates the chosen method of energy balance. However, the value of Gcdeduced 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 Gcinp<Gcsimwith 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 G c sim derived from the velocity and mass data obtained in the simulation versus the input parameter G c inp of the experiment. The dotted line marks the identity G c sim = G c inp ; the dashed line is approximated with R 2 ≈ 0.9928 and corresponds roughly to G c sim = 1.5 G c inp .

6.3 Phase-field fracture

For phase-field simulations we use a finer mesh of 5760 elements, h=5/6 mmand set ε=2.5 mmand M=0.6 μs1. At first we accelerate the specimen to vspec=6.8m/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 Gcinp, 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 st1and st2. Here we consider st1and set for symmetry st2=1st1. A lower “crack initiation” point st1clearly needs less time and energy to open the crack. On the one hand, if crack initiation equals crack opening, st1=st2=0.5, that is the crack appears immediately after initiation and the dissipated energy Gcsimwill be close to zero. The input value Gcinpwould 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 st1>0.95, st2>0.05the crack may not fully open, which gives inaccurate fragment velocities. A stable state we found for st10.75,0.9and so we set st1=0.8and st2=0.2.

In Figure 17 the computed specific fracture energy Gcsimover the input value Gcinpis shown. Again, the optimal case Gcinp=Gcsimis marked with a dotted line and the computed values show a clear proportionality between Gcinpand Gcsim. The value of Gcsimis now lower than the input value Gcinpwith 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 Gcinp<20N/m the specimen cracks basically without dissipation whereas too high values of Gcinpprevent any crack. In the range of interest, Gcinp=30100N/m both values Gcsimand Gcinpcorrelate very well with an R-coefficient of 0.9969. Specifically we have Gcsim0.9Gcinp. 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 G c sim derived from the velocity and mass data obtained simulations with ε = 5 / 2   mm , s t 1 = 0.8 , s t 2 = 0.2 versus the input parameter G c inp of the experiment. The dotted line marks the identity G c sim = G c inp ; the dashed line is approximated with R 2 ≈ 0.9969 and corresponds well to G c sim = 0.9 G c 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, tδ, 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 σc, the critical separation δc, and the weight β, which relates shear and tension. These parameters depend on the specific material and can be determined experimentally whereby δcis implicitly given via the specific crack energy Gc. Further specifications of the cohesive law may require additional material parameters, for example, in the classical exponential Rose-Ferrante law an additional parameter δ0needs 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 Ψeusof the material. The remaining relevant parameters are the mobility M, the specific crack energy Gc, and the length-scale parameter ε. The mobility has only numerical character and controls the rate of phase-field decrease during crack formation. The critical length-scale parameter ε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 Gc/εin the crack resistance and therefore εneeds to be calibrated carefully. The only material parameter in the phase-field model with a clear physical background is the Griffith energy Gc.

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 ε>h.

For numerical computation of the coupled fields u,sthere exist two different types of solution. On the one hand it is possible to determine the displacements uand 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<ε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 δ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 ε  δmodels the crack. The values 0<s<1can not be interpreted physically so that threshold values have to be chosen. Crack closure needs to be prevented by explicit local constraints, for example ṡ<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.

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Tim Dally, Carola Bilgen, Marek Werner and Kerstin Weinberg (May 6th 2020). Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture Analyses?, Modeling and Simulation in Engineering - Selected Problems, Jan Valdman and Leszek Marcinkowski, IntechOpen, DOI: 10.5772/intechopen.92180. Available from:

chapter statistics

130total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Methods of Nonequilibrium Statistical Mechanics in Models for Mixing Bulk Components

By Anna Kapranova, Daria Bahaeva, Dmitry Stenko, Ivan Verloka, Anton Lebedev and Mikhail Tarshis

Related Book

First chapter

Digital Image Processing with MATLAB

By Mahmut Sinecen

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us