Modeling in Continuum Mechanics

**Chapter 1**

**Abstract**

algorithms

**1. Introduction**

approaches [8–13].

**3**

Cohesive Elements or Phase-Field

Fracture: Which Method Is Better

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

titative 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

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–3], the extended finite element method [4, 5], eroded finite elements or eigenfracture strategies [6, 7], and phase-field

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 Γ*C*ð Þ*t* ⊂ ∂Ω in a solid of domain Ω undergoing a deformation *χ*ð Þ *x*, *t* : <sup>Ω</sup> � 0, *<sup>t</sup>*total ½ �! <sup>ℝ</sup><sup>3</sup> in a time *<sup>t</sup>*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

for Dynamic Fracture Analyses?

*Tim Dally, Carola Bilgen, Marek Werner*

*and Kerstin Weinberg*

### **Chapter 1**

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

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

### **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–3], the extended finite element method [4, 5], eroded finite elements or eigenfracture strategies [6, 7], and phase-field approaches [8–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 Γ*C*ð Þ*t* ⊂ ∂Ω in a solid of domain Ω undergoing a deformation *χ*ð Þ *x*, *t* : <sup>Ω</sup> � 0, *<sup>t</sup>*total ½ �! <sup>ℝ</sup><sup>3</sup> in a time *<sup>t</sup>*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 *ε*. Here an additional field *s*ð Þ *x*, *t* : Ω � 0, *t*total ½ �! ℝ characterizes the state of the material and marks the intact or broken state. The set of evolving crack surfaces is replaced by a cracksurface 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,

$$\int\_{\Gamma\_{C}(t)} \mathbf{d}\Gamma \approx \int\_{\Omega} \mathbf{\gamma}(t) \, \mathbf{d}\Omega. \tag{1}$$

distribution in the fast (cracking) specimen; and (iii) we need to quantify the

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.

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

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

We consider a body of domain <sup>Ω</sup> <sup>⊂</sup> <sup>ℝ</sup><sup>3</sup> with external boundary <sup>Γ</sup> � <sup>∂</sup>Ω. The body's displacement field at point *x* and time *t* is denoted by *u x*ð Þ , *t* ; its velocity and acceleration fields are *v* ¼ *u*\_ and *a* ¼ *u*€. A crack splits the body into subbodies Ω<sup>þ</sup> ∪ Ω� ¼ Ω and induces a displacement jump on Γ*<sup>C</sup>* as ½ �� ¼ ½*u u*<sup>þ</sup> � *u*�. The displacements satisfy the Dirichlet boundary conditions *u* ¼ *u* at Γ1. The body is

Linear-elastic material is presumed to follow Hooke's law with elastic strain

where the Lamé material parameters *λ* ¼ *Eν=*ð Þ ð Þ 1 þ *ν* ð Þ 1 � 2*ν* and *μ* ¼ *E=*ð Þ 2 1ð Þ þ *ν* are formulated with Young's modulus *E* and Poisson number *ν*. Equivalently we use the Hookean material tensor with components *Cijkl* ¼ *λδijδkl* þ *μ δikδjl* <sup>þ</sup> *<sup>δ</sup>il <sup>δ</sup>jk* . The strain tensor *<sup>ε</sup>* describes small deformations *<sup>ε</sup>*ð Þ¼ *<sup>u</sup>*

<sup>1</sup>*=*<sup>2</sup> **<sup>∇</sup>***<sup>u</sup>* <sup>þ</sup> **<sup>∇</sup>***u<sup>T</sup>* ; the mechanical stress tensor follows as *<sup>σ</sup>* <sup>¼</sup> *<sup>∂</sup>ε*Ψ*<sup>e</sup>* <sup>¼</sup> *ε*. It holds the

where *ρ* is the mass density and *b* a prescribed body force density. Boundary

<sup>2</sup> <sup>þ</sup> *<sup>μ</sup><sup>ε</sup>* : *<sup>ε</sup>* <sup>¼</sup> <sup>1</sup>

2

div*σ* þ *b* ¼ *ρa* in Ω (3)

*ε* : : *ε* (2)

loaded with traction *t* at boundary Γ2; it holds Γ<sup>1</sup> ∪ Γ<sup>2</sup> ⊂ΓnΓ*C*.

<sup>Ψ</sup>*<sup>e</sup>* <sup>¼</sup> <sup>1</sup>

<sup>2</sup> *<sup>λ</sup>*ð Þ tr *<sup>ε</sup>*

fracture energy and the critical energy release rate of the material.

**2. Governing equations**

**Figure 1.**

**2.1 Elasto-dynamics**

balance of linear momentum,

conditions are prescribed as

**5**

energy density,

By *γ* being not only a function of field *s*ð Þ *x*, *t* but 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 wellknown 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–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

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

#### **Figure 1.**

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

ð

*Modeling and Simulation in Engineering - Selected Problems*

Γ*C*ð Þ*t*

performed with concrete specimen in a Hopkinson-Bar (HB) setup.

dΓ≈ ð Ω

By *γ* being not only a function of field *s*ð Þ *x*, *t* but 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 wellknown 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

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

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

HB spalling experiments are test arrangements to determine the failure strength of brittle materials, see [15–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

*γ*ð Þ*t* dΩ*:* (1)

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 *s*ð Þ *x*, *t* : Ω � 0, *t*total ½ �! ℝ characterizes the state of the material and marks the intact or broken state. The set of evolving crack surfaces is replaced by a cracksurface 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

and their further movement.

boundaries over the body's domain,

quantify the mechanics of fracture.

reliability of the obtained material data.

**4**

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

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 <sup>Ω</sup> <sup>⊂</sup> <sup>ℝ</sup><sup>3</sup> with external boundary <sup>Γ</sup> � <sup>∂</sup>Ω. The body's displacement field at point *x* and time *t* is denoted by *u x*ð Þ , *t* ; its velocity and acceleration fields are *v* ¼ *u*\_ and *a* ¼ *u*€. A crack splits the body into subbodies Ω<sup>þ</sup> ∪ Ω� ¼ Ω and induces a displacement jump on Γ*<sup>C</sup>* as ½ �� ¼ ½*u u*<sup>þ</sup> � *u*�. The displacements satisfy the Dirichlet boundary conditions *u* ¼ *u* at Γ1. The body is loaded with traction *t* at boundary Γ2; it holds Γ<sup>1</sup> ∪ Γ<sup>2</sup> ⊂ΓnΓ*C*.

### **2.1 Elasto-dynamics**

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

$$\Psi^{\mathbf{e}} = \frac{1}{2} \left. \lambda(\mathbf{r} \mathbf{e}) \right|^2 + \mu \mathbf{e} : \mathbf{e} = \frac{1}{2} \left. \mathbf{e} : \mathbb{C} : \mathbf{e} \right. \tag{2}$$

where the Lamé material parameters *λ* ¼ *Eν=*ð Þ ð Þ 1 þ *ν* ð Þ 1 � 2*ν* and *μ* ¼ *E=*ð Þ 2 1ð Þ þ *ν* are formulated with Young's modulus *E* and Poisson number *ν*. Equivalently we use the Hookean material tensor with components *Cijkl* ¼ *λδijδkl* þ *μ δikδjl* <sup>þ</sup> *<sup>δ</sup>il <sup>δ</sup>jk* . The strain tensor *<sup>ε</sup>* describes small deformations *<sup>ε</sup>*ð Þ¼ *<sup>u</sup>* <sup>1</sup>*=*<sup>2</sup> **<sup>∇</sup>***<sup>u</sup>* <sup>þ</sup> **<sup>∇</sup>***u<sup>T</sup>* ; the mechanical stress tensor follows as *<sup>σ</sup>* <sup>¼</sup> *<sup>∂</sup>ε*Ψ*<sup>e</sup>* <sup>¼</sup> *ε*. It holds the balance of linear momentum,

$$
\rho \operatorname{div} \sigma + \overline{\mathbf{b}} = \rho \mathfrak{a} \quad \text{in} \ \Omega \tag{3}
$$

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

$$
\mu = \overline{\mathfrak{u}} \operatorname{ on} \Gamma\_1, \qquad \sigma \mathfrak{u} = \overline{\mathfrak{t}} \operatorname{ on} \Gamma\_2, \qquad \left[ \sigma \mathfrak{u} \right] = \left[ \mathfrak{t} \right] = \mathbf{0} \operatorname{ on} \Gamma\_C, \tag{4}
$$

and with *δ*Π ¼ *δW*int � *δW*ext and Eqs. (2)–(4) we obtain

ð Γ*<sup>C</sup>*

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

denotes the Hilbert space of weak square-integrable functions and its first derivatives. Eq. (8) must hold for all times *t*≥*t*0, which leads to the weak momentum

� � <sup>d</sup>*<sup>t</sup>* <sup>¼</sup> 0 (8)

*δ*½ �� � ½*u t* dΓ �

ð Γ2

*i*¼1

*e*¼1

**Nt** dΓ*:*

ð Γ1

Discretization in time is performed by an implicit Euler method, that is, with time step Δ*t* ¼ *tn*þ<sup>1</sup> � *tn* and, thus, velocity and acceleration are approximated by

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

ð

Ω*e*

For discretization the domain Ω is subdivided into a finite set of nonoverlapping elements. We make use of a conforming ansatz and approximate the displacement

*Niu<sup>i</sup>* <sup>¼</sup> **Nu**^, *<sup>δ</sup>u*<sup>≈</sup> <sup>X</sup>*nk*

**M**€

*<sup>ρ</sup>***NN***<sup>T</sup>* <sup>d</sup>Ω*e*, **<sup>K</sup>** <sup>¼</sup> ∇ ⋃*ne*

**u**^ *n*þ1

**Nb** þ **Nt**dΩ þ

where *Ni* are the piecewise ansatz functions collected in a matrix **N** and the vectors **u**^, **w**^ contain all unknown nodal displacements *u<sup>i</sup>* of the *nk* nodes and its nodal variations *wi*, respectively. Plugging Eq. (10) into Eq. (8) and a straightfor-

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 is

ð Γ2

*δu* � *t* dΓ

*δu* � *t* dΓ ∀*δu*∈ U*:* (9)

*Niw<sup>i</sup>* ¼ **Nw**^ , (10)

**u**^ þ **Ku**^ ¼ **f** (11)

**N∇N***<sup>T</sup>* dΩ*e*,

<sup>≈</sup> **<sup>u</sup>**^*<sup>n</sup>*þ<sup>1</sup> � <sup>2</sup>**u**^*<sup>n</sup>* <sup>þ</sup> **<sup>u</sup>**^*<sup>n</sup>*�<sup>1</sup> � �ð Þ <sup>Δ</sup>*<sup>t</sup>* �<sup>2</sup>

*:* (12)

ð Þ <sup>Ω</sup> <sup>∣</sup>*δ<sup>u</sup>* <sup>¼</sup> 0 on <sup>Γ</sup>1g; *<sup>H</sup>*<sup>1</sup>

ð Þ� *u δu* � *b*dΩ �

for all admissible test functions *<sup>δ</sup>u*<sup>∈</sup> <sup>U</sup> with <sup>U</sup> ¼ f*δu*<sup>∈</sup> *<sup>H</sup>*<sup>1</sup>

balance. For a body with stationary crack it is equivalent to

*δu* � *ρa* þ *δε* : *σ* � *δu* � *b* dΩ ¼

ð*t t*0 ð Ω

*<sup>ρ</sup><sup>v</sup>* � *<sup>δ</sup><sup>v</sup>* <sup>þ</sup> *<sup>δ</sup>*Ψ*<sup>e</sup>*

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

ð Ω

field and its variation with

*<sup>u</sup>*<sup>≈</sup> <sup>X</sup>*nk i*¼1

reformulated in a matrix , and we specify

ð

Ω*e*

ð

Ω*e*∩Γ<sup>1</sup>

<sup>≈</sup> **<sup>u</sup>**^*<sup>n</sup>*þ<sup>1</sup> � **<sup>u</sup>**^*<sup>n</sup>* � �ð Þ <sup>Δ</sup>*<sup>t</sup>* �<sup>1</sup> and €

**M** ¼ ⋃ *ne e*¼1

**f** ¼ ⋃ *ne e*¼1

**3. The cohesive element technique**

\_ **u**^ *n*þ1

**7**

ward calculation gives the global system of equations

with normal vector *n* and 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 Γ*C*ð Þ*t* . 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

$$\mathcal{W}\_{\text{int}}(\mathfrak{u},t) = \int\_{\mathfrak{U}} \Psi^{\epsilon}(\mathfrak{u}) \, \text{d}\mathfrak{Q} + \int\_{\Gamma(t)} \mathcal{G}\_{\mathfrak{C}} \mathrm{d}\Gamma,\tag{5}$$

where G*<sup>c</sup>* is commonly known as Griffith's critical energy release rate (Griffith energy). An optimum of Eq. (5) corresponds to crack growth.

The specific energy G corresponds to the energy *W* dissipated per unit newly created surface area *A*, G ¼ d*W=*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, G ≥ G*c*, which is related to the classical stress intensity factors of linearelastic fracture mechanics given by <sup>G</sup> <sup>¼</sup> *<sup>K</sup>*<sup>2</sup> *<sup>I</sup>* <sup>þ</sup> *<sup>K</sup>*<sup>2</sup> *II* � �*=E*<sup>0</sup> <sup>þ</sup> *<sup>K</sup>*<sup>2</sup> *III=*ð Þ 2*μ* with *E*<sup>0</sup> ¼ *E* regarding the plane stress state and otherwise *<sup>E</sup>*<sup>0</sup> <sup>¼</sup> *<sup>E</sup><sup>=</sup>* <sup>1</sup> � *<sup>ν</sup>*<sup>2</sup> ð Þ. Thus, Griffith's criterion can also be written with the critical stress intensity factor *Kc* as <sup>G</sup>*<sup>c</sup>* <sup>¼</sup> *<sup>K</sup>*<sup>2</sup> *c=E*<sup>0</sup> . The single summands of G may also be understood as mode-dependent Griffith energy densities G*I*, G*II*, and G*III*. 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 J *<sup>c</sup>* ¼ G*c*.

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 *t* and 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 *δ* ¼ ∣*δ*∣,

$$\mathcal{G}\_{\mathfrak{c}} = \int\_{\Gamma} \mathfrak{t} \cdot \mathfrak{d} \quad \mathrm{d}\Gamma = \int\_{0}^{\delta\_{\mathfrak{c}}} \mathfrak{t} \cdot \mathrm{d}\delta. \tag{6}$$

#### **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 ½ � *t*0, *t* defines a functional I½ �¼ *<sup>u</sup>* <sup>Ð</sup>*<sup>t</sup> t*0 Ld*t* with the Lagrangian function as the difference of kinetic and potential energy, Lð Þ¼K *u*, *v* ð Þ� *v* Πð Þ *u* . Stationarity demands the first variation to vanish,

$$
\delta \mathcal{L} = \delta \int\_{t\_0}^t \mathcal{L}(\mathfrak{u}, \mathfrak{v}) \, \mathrm{d}t = \int\_{t\_0}^t \delta \mathcal{K}(\mathfrak{v}) - \delta \Pi(\mathfrak{u}) \, \mathrm{d}t = \mathbf{0} \qquad \forall t \ge t\_0,\tag{7}
$$

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

and with *δ*Π ¼ *δW*int � *δW*ext and Eqs. (2)–(4) we obtain

$$\int\_{t\_0}^t \left( \int\_{\Omega} \rho \mathbf{v} \cdot \delta \mathbf{v} + \delta \Psi^\epsilon(\mathbf{u}) - \delta \mathbf{u} \cdot \overline{\mathbf{b}} \, \mathrm{d}\Omega - \int\_{\Gamma\_0} \delta [\mathbf{u}] \cdot \mathbf{t} \, \mathrm{d}\Gamma - \int\_{\Gamma\_2} \delta \mathbf{u} \cdot \overline{\mathbf{t}} \, \mathrm{d}\Gamma \right) \mathrm{d}t = \mathbf{0} \tag{8}$$

for all admissible test functions *<sup>δ</sup>u*<sup>∈</sup> <sup>U</sup> with <sup>U</sup> ¼ f*δu*<sup>∈</sup> *<sup>H</sup>*<sup>1</sup> ð Þ <sup>Ω</sup> <sup>∣</sup>*δ<sup>u</sup>* <sup>¼</sup> 0 on <sup>Γ</sup>1g; *<sup>H</sup>*<sup>1</sup> denotes the Hilbert space of weak square-integrable functions and its first derivatives. Eq. (8) must hold for all times *t*≥*t*0, which leads to the weak momentum balance. For a body with stationary crack it is equivalent to

$$\int\_{\Omega} \delta \boldsymbol{\mathfrak{u}} \cdot \rho \boldsymbol{\mathfrak{a}} + \delta \boldsymbol{\mathfrak{e}} : \boldsymbol{\sigma} - \delta \boldsymbol{\mathfrak{u}} \cdot \overline{\boldsymbol{\mathfrak{b}}} \, \mathrm{d}\Omega = \int\_{\Gamma\_2} \delta \boldsymbol{\mathfrak{u}} \cdot \overline{\boldsymbol{\mathfrak{t}}} \, \mathrm{d}\Gamma \qquad \forall \delta \boldsymbol{\mathfrak{u}} \in \mathsf{U}. \tag{9}$$

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

$$
\mathfrak{u} \approx \sum\_{i=1}^{n\_k} N\_i \mathfrak{u}\_i = \mathbf{N} \hat{\mathbf{u}}, \qquad \delta \mathfrak{u} \approx \sum\_{i=1}^{n\_k} N\_i \mathfrak{w}\_i = \mathbf{N} \hat{\mathbf{w}}, \tag{10}
$$

where *Ni* are the piecewise ansatz functions collected in a matrix **N** and the vectors **u**^, **w**^ contain all unknown nodal displacements *u<sup>i</sup>* of the *nk* nodes and its nodal variations *wi*, respectively. Plugging Eq. (10) into Eq. (8) and a straightforward calculation gives the global system of equations

$$
\overline{\mathbf{M}}\ddot{\mathbf{u}} + \overline{\mathbf{K}}\hat{\mathbf{u}} = \overline{\mathbf{f}} \tag{11}
$$

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 is reformulated in a matrix , and we specify

$$\begin{split} \overline{\mathbf{M}} &= \bigcup\_{\epsilon=1}^{n\_{\epsilon}} \int\_{\Omega\_{\epsilon}} \rho \mathbf{N} \mathbf{N}^{T} \, \mathrm{d}\Omega\_{\epsilon}, \quad \overline{\mathbf{K}} = \nabla \bigcup\_{\epsilon=1}^{n\_{\epsilon}} \int\_{\Omega\_{\epsilon}} \mathbf{N} \overline{\mathbf{C}} \mathbf{V} \mathbf{N}^{T} \, \mathrm{d}\Omega\_{\epsilon}, \\ \overline{\mathbf{f}} &= \bigcup\_{\epsilon=1}^{n\_{\epsilon}} \int\_{\Omega\_{\epsilon} \cap \Gamma\_{1}} \mathbf{N} \overline{\mathbf{b}} + \mathbf{N} \overline{\mathbf{t}} \, \mathrm{d}\Omega + \int\_{\Gamma\_{1}} \mathbf{N} \overline{\mathbf{t}} \, \mathrm{d}\Gamma. \end{split}$$

Discretization in time is performed by an implicit Euler method, that is, with time step Δ*t* ¼ *tn*þ<sup>1</sup> � *tn* and, thus, velocity and acceleration are approximated by

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

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

*u* ¼ *u*onΓ1, *σn* ¼ *t* onΓ2, ½ �� ¼ ½ ½*σn* ½ �� ¼ *t* **0**onΓ*C*, (4)

with normal vector *n* and traction *t*; ½ �� ½� denotes a jump. Eq. (4)3 corresponds to

Let the evolving internal cracks be represented by a set of boundaries Γ*C*ð Þ*t* . 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

ð Þ *u* dΩ þ

*<sup>I</sup>* <sup>þ</sup> *<sup>K</sup>*<sup>2</sup> *II* � �*=E*<sup>0</sup> <sup>þ</sup> *<sup>K</sup>*<sup>2</sup>

where G*<sup>c</sup>* is commonly known as Griffith's critical energy release rate (Griffith

The specific energy G corresponds to the energy *W* dissipated per unit newly created surface area *A*, G ¼ d*W=*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, G ≥ G*c*, which is related to the classical stress intensity factors of linear-

regarding the plane stress state and otherwise *<sup>E</sup>*<sup>0</sup> <sup>¼</sup> *<sup>E</sup><sup>=</sup>* <sup>1</sup> � *<sup>ν</sup>*<sup>2</sup> ð Þ. Thus, Griffith's crite-

single summands of G may also be understood as mode-dependent Griffith energy densities G*I*, G*II*, and G*III*. 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

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 *t* and 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

rion can also be written with the critical stress intensity factor *Kc* as <sup>G</sup>*<sup>c</sup>* <sup>¼</sup> *<sup>K</sup>*<sup>2</sup>

to the strain energy release rate, that is, it holds the identity J *<sup>c</sup>* ¼ G*c*.

surrounds the cohesive zone gives for an effective crack opening *δ* ¼ ∣*δ*∣,

*t* � *δ* dΓ ¼

The motions of a solid can be characterized by recourse to Hamilton's principle

kinetic and potential energy, Lð Þ¼K *u*, *v* ð Þ� *v* Πð Þ *u* . Stationarity demands the first

of stationary action. The action of a motion within a closed time interval ½ � *t*0, *t*

ð*t t*0 ð*δc* 0

Ld*t* with the Lagrangian function as the difference of

*δ*Kð Þ� *v δ*Πð Þ *u* d*t* ¼ 0 ∀*t* ≥*t*0, (7)

G*<sup>c</sup>* ¼ ð Γ

*t*0

Lð Þ *u*, *v* d*t* ¼

defines a functional I½ �¼ *<sup>u</sup>* <sup>Ð</sup>*<sup>t</sup>*

*δ*I ¼ *δ*

ð*t t*0

variation to vanish,

**6**

**2.3 Weak form of the problem and finite element discretization**

ð

Γð Þ*t*

G*<sup>c</sup>* dΓ, (5)

*III=*ð Þ 2*μ* with *E*<sup>0</sup> ¼ *E*

*t* d*δ:* (6)

*c=E*<sup>0</sup> . The

traction-free crack boundaries. Additional initial conditions may apply.

ð Ω Ψ*e*

*W*intð Þ¼ *u*, *t*

*Modeling and Simulation in Engineering - Selected Problems*

energy). An optimum of Eq. (5) corresponds to crack growth.

elastic fracture mechanics given by <sup>G</sup> <sup>¼</sup> *<sup>K</sup>*<sup>2</sup>

**2.2 Fracture mechanics**

is composed of

**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* Γ� *<sup>C</sup> and* Γ<sup>þ</sup> *<sup>C</sup> coincide if the element is closed.*

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

### **3.1 Fundamentals**

In cohesive theories, the displacement jump across a cohesive surface Γ*<sup>C</sup>*

$$\delta = [\mathfrak{u}] = \mathfrak{u}^+ - \mathfrak{u}^- \tag{13}$$

^*t*eff <sup>¼</sup> ^*tn* <sup>þ</sup> *<sup>β</sup>*�<sup>2</sup>

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

ized by the conditions *<sup>δ</sup>* <sup>¼</sup> *<sup>δ</sup>*max and \_

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

monly assumed to be linear in *δ* to the origin,

determine via Eq. (6) the specific fracture energy

Ferrante [33].

modeled here.

**Figure 3.**

**9**

^*t*ð Þ¼ *<sup>δ</sup>* ^*t*max *δ*max <sup>j</sup>*tt*<sup>j</sup> <sup>1</sup>*=*<sup>2</sup>

*δ*, if *δ*< *δ*max or \_

*δc*

An appropriate choice of cohesive variable is the maximum attained (effective) crack opening displacement *δ*max. Loading of the cohesive surface is then character-

**Figure 3** shows different common loading envelopes, whereas unloading is com-

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

^*t*ð Þ¼ *<sup>δ</sup> <sup>σ</sup><sup>c</sup>* <sup>1</sup> � *<sup>δ</sup>*

The two parameter cohesive strength *σ<sup>c</sup>* and critical opening displacement *δ<sup>c</sup>*

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

exp 1 � *<sup>δ</sup>*

*δ*0

*δ δ*0

where *δ*<sup>0</sup> is 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 <sup>G</sup>*<sup>c</sup>* <sup>¼</sup> *<sup>σ</sup>cδ*0ð Þ exp 1ð Þ� ð Þ <sup>1</sup> <sup>þ</sup> *<sup>δ</sup>c=δ*<sup>0</sup> exp 1ð Þ � *<sup>δ</sup>c=δ*<sup>0</sup> . Unloading at ^*t*<^*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

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

^*t*ð Þ¼ *<sup>δ</sup> <sup>σ</sup><sup>c</sup>*

*:* (15)

*δ*<0*:* (16)

*δ*≥0; all other states correspond to unloading.

*:* (17)

, (19)

G*<sup>c</sup>* ¼ 1*=*2 *σ<sup>c</sup> δc:* (18)

plays the role of a deformation measure while the tractions *t* furnish the workconjugate 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 Γ*<sup>C</sup>* have to be treated differently. Given a vector field *u* over Γ*C*, its normal and tangential components are *un* ¼ *u* � *n*, and *ut* ¼ ∣*u*Γ∣ ¼ ∣*u* � *unn*∣ ¼ ∣ð Þ *I* � *n* ⊗ *n u*∣, and, the corresponding jump components *δ<sup>n</sup>* ¼ *u*<sup>þ</sup> *<sup>n</sup>* � *u*� *<sup>n</sup>* and *δ<sup>t</sup>* ¼ *u*þ *<sup>t</sup>* � *u*� *<sup>t</sup>* follow from Eq. (13). To further simplify the formulation of mixed-mode cohesive laws, we follow [25] and introduce an effective opening displacement

$$\delta = |\mathfrak{G}| = \left(\beta^2 \delta\_t + \delta\_n\right)^{1/2}.\tag{14}$$

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 *<sup>t</sup>* <sup>¼</sup> *<sup>∂</sup>δ*G, where <sup>G</sup> 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 G has the form G ¼ Gð Þ¼ *δn*, *δ* � *u* Gð Þ *δn*, *δ<sup>t</sup>* and 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 **t** to a single scalar function. Then an effective traction <sup>∣</sup>*t*<sup>∣</sup> <sup>¼</sup> ^*<sup>t</sup>* <sup>¼</sup> ^*t*eff follows,

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

$$\hat{\mathbf{t}}\_{\rm eff} = \left(\hat{\mathbf{t}}\_n + \boldsymbol{\beta}^{-2} |\mathbf{t}\_l|\right)^{1/2}. \tag{15}$$

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 *<sup>δ</sup>* <sup>¼</sup> *<sup>δ</sup>*max and \_ *δ*≥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,

$$
\hat{t}(\delta) = \frac{\hat{t}\_{\text{max}}}{\delta\_{\text{max}}} \delta, \quad \text{if} \ \delta < \delta\_{\text{max}} \text{ or } \dot{\delta} < 0. \tag{16}
$$

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

$$
\hat{t}(\delta) = \sigma\_{\hat{\mathfrak{e}}} \left( \mathbf{1} - \frac{\delta}{\delta\_{\mathfrak{e}}} \right). \tag{17}
$$

The two parameter cohesive strength *σ<sup>c</sup>* and critical opening displacement *δ<sup>c</sup>* determine via Eq. (6) the specific fracture energy

$$
\mathcal{G}\_{\mathfrak{c}} = \mathbf{1}/\mathfrak{2} \; \; \sigma\_{\mathfrak{c}} \; \; \delta\_{\mathfrak{c}}.\tag{18}
$$

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

$$
\hat{t}(\delta) = \sigma\_\varepsilon \frac{\delta}{\delta\_0} \exp\left(\mathbf{1} - \frac{\delta}{\delta\_0}\right),
\tag{19}
$$

where *δ*<sup>0</sup> is 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 <sup>G</sup>*<sup>c</sup>* <sup>¼</sup> *<sup>σ</sup>cδ*0ð Þ exp 1ð Þ� ð Þ <sup>1</sup> <sup>þ</sup> *<sup>δ</sup>c=δ*<sup>0</sup> exp 1ð Þ � *<sup>δ</sup>c=δ*<sup>0</sup> . Unloading at ^*t*<^*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.

#### **Figure 3.**

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

approach combined with an automatic fragmentation and cohesive surface insertion procedure. The method has proven to be reliable and efficient for numerous

*Discretization of a solid with the cohesive element technique (left) and geometry of a 2D and 3D cohesive*

*<sup>C</sup> and* Γ<sup>þ</sup>

In cohesive theories, the displacement jump across a cohesive surface Γ*<sup>C</sup>*

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 Γ*<sup>C</sup>* have to be treated differently. Given a vector field *u* over Γ*C*, its normal and tangential components are *un* ¼ *u* � *n*, and *ut* ¼ ∣*u*Γ∣ ¼ ∣*u* � *unn*∣ ¼

∣ð Þ *I* � *n* ⊗ *n u*∣, and, the corresponding jump components *δ<sup>n</sup>* ¼ *u*<sup>þ</sup>

scalar function. Then an effective traction <sup>∣</sup>*t*<sup>∣</sup> <sup>¼</sup> ^*<sup>t</sup>* <sup>¼</sup> ^*t*eff follows,

plays the role of a deformation measure while the tractions *t* furnish the workconjugate 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

*<sup>t</sup>* follow from Eq. (13). To further simplify the formulation of mixed-mode

*δ<sup>t</sup>* þ *δ<sup>n</sup>* <sup>1</sup>*=*<sup>2</sup>

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

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 *<sup>t</sup>* <sup>¼</sup> *<sup>∂</sup>δ*G, where <sup>G</sup> 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 G has the form G ¼ Gð Þ¼ *δn*, *δ* � *u* Gð Þ *δn*, *δ<sup>t</sup>* and 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 **t** to a single

cohesive laws, we follow [25] and introduce an effective opening displacement

*<sup>δ</sup>* <sup>¼</sup> <sup>∣</sup>*δ*<sup>∣</sup> <sup>¼</sup> *<sup>β</sup>*<sup>2</sup>

*δ* ¼ ½½ �� ¼ *u u*<sup>þ</sup> � *u*� (13)

*<sup>C</sup> coincide if the element is closed.*

*<sup>n</sup>* � *u*�

*:* (14)

*<sup>n</sup>* and *δ<sup>t</sup>* ¼

applications, see among others [2, 28–30].

*Modeling and Simulation in Engineering - Selected Problems*

*element, respectively (right). The surfaces* Γ�

**3.1 Fundamentals**

**Figure 2.**

*u*þ *<sup>t</sup>* � *u*�

**8**

opening displacement *δ* only.

**3.2 Cohesive laws**

### **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*coh nodes; the total number of nodes of the cohesive element is 2*n*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 \mathcal{S} W\_{\text{int}}(\mathfrak{u}) = \int\_{\mathfrak{U}} \delta \mathfrak{e} : \mathbb{C} : \mathfrak{e} \,\mathrm{d}\Omega + \delta W\_{\text{int}}^{\text{coh}}(\mathfrak{u}) \text{ with } \delta W\_{\text{int}}^{\text{coh}}(\mathfrak{u}) := \int\_{\Gamma\_{\text{coh}}} \delta[\mathfrak{u}] \cdot \mathfrak{t} \,\mathrm{d}\Gamma \tag{20}$$

where Γcoh are the total cohesive surfaces and *δ*½ �� ½*u* is 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

$$\delta \mathcal{W}\_{\text{int},\varepsilon}^{\text{coh}}(\mathfrak{u}) = \int\_{\Gamma\_{\text{coh},\varepsilon}} \delta \Big(\mathbf{N} \, \|\hat{\mathbf{u}}\|\_{\varepsilon}\Big)^{T} \mathbf{t} \, \mathrm{d}\Gamma = \delta \|\hat{\mathbf{u}}\|\_{\varepsilon} \, ^{T} \Big[\, \underset{\Gamma\_{\text{coh},\varepsilon}}{\mathbf{N}^{T}\mathbf{t}} \, \mathrm{d}\Gamma. \tag{21}$$

of the body, the edge with all nodes will be duplicated: *k*1*k*2∩ ∂ð Þ¼ Ω ∪ Γð Þ*t* ∅; case 2:

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

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

Here we illustrate the influence of the cohesive strength *σ<sup>c</sup>* and critical opening displacement *δ<sup>c</sup>* on the cohesive element simulations. Exemplarily we investigate a plane mode-I tension test under quasistatic conditions; however, the observations

The test specimen with material data *E* ¼ 60 GPa, *ν* ¼ 0*:*2 has a unit domain ½ �� 0, 1 ½ � 0, 1 ; 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, *ty* ¼ 15 MPa. A crack is predefined in the area *x*∈½ � 0, 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 *σ<sup>c</sup>* ¼ 20 MPa and the critical crack opening displacement *<sup>δ</sup><sup>c</sup>* <sup>¼</sup> <sup>10</sup>�4. Following Eq. (18) this corresponds to a critical energy release rate of

During the simulation a cohesive element will be added if the effective traction

*Crack propagation of the mode-I tension test with x*ð Þ , *<sup>y</sup>* <sup>∈</sup> ½ �� 0, 1 ½ � 0, 1 *<sup>m</sup><sup>2</sup> 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 σ<sup>y</sup> [MPa] is plotted.*

Eq. (15), here ^*<sup>t</sup>* <sup>¼</sup> ^*t*effð Þ *<sup>β</sup>* ! <sup>∞</sup> , exceeds a critical value, ^*<sup>t</sup>* <sup>≥</sup>*σi*. We identify the insertion criterium *σ<sup>i</sup>* with 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 *σ<sup>y</sup>* in **Figure 5** correspond very well to the analytical solution.

one of the marked nodes is on the boundary, one node will be duplicated: *k*<sup>1</sup> ∈ ∂Ω ∪ Γð Þ*t* ⊻ *k*<sup>2</sup> ∈ ∂Ω ∪ Γð Þ*t* ; and case 3: both nodes are on the boundary, two nodes (all nodes) will be duplicated: *k*<sup>1</sup> ∈ ∂Ω ∪ Γð Þ*t* ∧ *k*<sup>2</sup> ∈ ∂Ω ∪ Γð Þ*t* , for the twodimensional 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

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

have been confirmed also in other models.

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

algorithms, cf. [35].

**Figure 4.**

*on the boundary.*

G*<sup>c</sup>* ¼ 1 N*=*mm.

**Figure 5.**

**11**

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 \mathcal{W}\_{\text{ext},\varepsilon}^{\text{coh}}(\mathfrak{u}) = \delta \mathcal{W}\_{\text{int},\varepsilon}^{\text{coh}} = \delta \left[ \hat{\mathbf{u}} \right]\_{\epsilon}^{T} \mathbf{f}\_{\text{coh},\varepsilon}^{\pm} \quad \text{with} \quad \mathbf{f}\_{\text{coh},\varepsilon}^{\pm} = \mp \int\_{\Gamma\_{\text{coh},\varepsilon}} \mathbf{N}^{T} \mathbf{t} \, d\Gamma. \tag{22}$$

The tangent stiffness matrix follows by its consistent linearization, with the result **<sup>K</sup>**coh,*<sup>e</sup>* <sup>¼</sup> <sup>Ð</sup> Γcoh,*<sup>e</sup>* **<sup>N</sup>***<sup>T</sup>∂***t***=∂*½ �� <sup>½</sup>*<sup>u</sup>* **<sup>N</sup>** <sup>d</sup>Γ. To avoid singularities at *<sup>δ</sup>* ! 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 **<sup>K</sup>**coh � <sup>Δ</sup>½ �� ¼ <sup>½</sup>**u**^ <sup>Δ</sup>**<sup>f</sup>** coh with **<sup>K</sup>**coh <sup>¼</sup> <sup>⋃</sup>*<sup>n</sup>*coh *<sup>e</sup>*¼<sup>1</sup>**K**coh,*<sup>e</sup>*, **<sup>f</sup>** coh <sup>¼</sup> <sup>⋃</sup>*<sup>n</sup>*coh *<sup>e</sup>*¼<sup>1</sup> **<sup>f</sup>** coh,*<sup>e</sup>*, where <sup>Δ</sup> symbolizes the incremental solution procedure, which is required for the nonlinear crackopening 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*<sup>1</sup> and *k*<sup>2</sup> is inside *Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

**Figure 4.**

**3.3 Finite element implementation**

conditions gives the deformation power with variation

*Modeling and Simulation in Engineering - Selected Problems*

*<sup>δ</sup><sup>ε</sup>* : : *<sup>ε</sup>* <sup>d</sup><sup>Ω</sup> <sup>þ</sup> *<sup>δ</sup>W*coh

Eq. (10) in Eq. (20)2 for one cohesive element to

ð

int,*<sup>e</sup>* ¼ *δ*½ �� ½**u**^ *<sup>e</sup>*

Γcoh,*<sup>e</sup>*

*δ* **N**½ �� ½**u**^ *<sup>e</sup>* � �*<sup>T</sup>*

*<sup>T</sup>***f**�

the external virtual work has to be determined. For one cohesive element it is

int,*<sup>e</sup>*ð Þ¼ *u*

*δW*coh

ext,*<sup>e</sup>*ð Þ¼ *<sup>u</sup> <sup>δ</sup>W*coh

Γcoh,*<sup>e</sup>*

as **<sup>K</sup>**coh � <sup>Δ</sup>½ �� ¼ <sup>½</sup>**u**^ <sup>Δ</sup>**<sup>f</sup>** coh with **<sup>K</sup>**coh <sup>¼</sup> <sup>⋃</sup>*<sup>n</sup>*coh

four nodes.

*δW*intð Þ¼ *u*

*δW*coh

result **<sup>K</sup>**coh,*<sup>e</sup>* <sup>¼</sup> <sup>Ð</sup>

opening problem.

**10**

**3.4 Adaptive meshing**

ð Ω

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*coh nodes; the total number of nodes of the cohesive element is 2*n*coh. The particular triangulation depicted in **Figure 2** is compatible with edge or line elements of

int ð Þ *<sup>u</sup>* with *<sup>δ</sup>W*coh

**t** dΓ ¼ *δ*½ �� ½**u**^ *<sup>e</sup>*

� coh,*<sup>e</sup>* ¼ ∓

**<sup>N</sup>***<sup>T</sup>∂***t***=∂*½ �� <sup>½</sup>*<sup>u</sup>* **<sup>N</sup>** <sup>d</sup>Γ. To avoid singularities at *<sup>δ</sup>* ! 0 in the deriva-

*<sup>e</sup>*¼<sup>1</sup>**K**coh,*<sup>e</sup>*, **<sup>f</sup>** coh <sup>¼</sup> <sup>⋃</sup>*<sup>n</sup>*coh

The kinetic energy does not have any support in the cohesive element and only

coh,*<sup>e</sup>* with **f**

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

tives 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

izes the incremental solution procedure, which is required for the nonlinear crack-

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*<sup>1</sup> and *k*<sup>2</sup> is inside

int ð Þ *u* ≔

*T* ð

Γcoh,*<sup>e</sup>*

ð

Γcoh,*<sup>e</sup>*

ð

Γcoh

*δ*½ �� � ½*u t*dΓ (20)

**N***<sup>T</sup>***t** dΓ*:* (21)

**N***<sup>T</sup>***t** dΓ*:* (22)

*<sup>e</sup>*¼<sup>1</sup> **<sup>f</sup>** coh,*<sup>e</sup>*, where <sup>Δ</sup> symbol-

Basis of the finite element implementation is Hamilton's principle given in Eq. (7). Inserting the balance of linear momentum and the static boundary

where Γcoh are the total cohesive surfaces and *δ*½ �� ½*u* is 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

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

of the body, the edge with all nodes will be duplicated: *k*1*k*2∩ ∂ð Þ¼ Ω ∪ Γð Þ*t* ∅; case 2: one of the marked nodes is on the boundary, one node will be duplicated: *k*<sup>1</sup> ∈ ∂Ω ∪ Γð Þ*t* ⊻ *k*<sup>2</sup> ∈ ∂Ω ∪ Γð Þ*t* ; and case 3: both nodes are on the boundary, two nodes (all nodes) will be duplicated: *k*<sup>1</sup> ∈ ∂Ω ∪ Γð Þ*t* ∧ *k*<sup>2</sup> ∈ ∂Ω ∪ Γð Þ*t* , for the twodimensional 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].

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

Here we illustrate the influence of the cohesive strength *σ<sup>c</sup>* and critical opening displacement *δ<sup>c</sup>* 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 GPa, *ν* ¼ 0*:*2 has a unit domain ½ �� 0, 1 ½ � 0, 1 ; 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, *ty* ¼ 15 MPa. A crack is predefined in the area *x*∈½ � 0, 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 *σ<sup>c</sup>* ¼ 20 MPa and the critical crack opening displacement *<sup>δ</sup><sup>c</sup>* <sup>¼</sup> <sup>10</sup>�4. Following Eq. (18) this corresponds to a critical energy release rate of G*<sup>c</sup>* ¼ 1 N*=*mm.

During the simulation a cohesive element will be added if the effective traction Eq. (15), here ^*<sup>t</sup>* <sup>¼</sup> ^*t*effð Þ *<sup>β</sup>* ! <sup>∞</sup> , exceeds a critical value, ^*<sup>t</sup>* <sup>≥</sup>*σi*. We identify the insertion criterium *σ<sup>i</sup>* with 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 *σ<sup>y</sup>* in **Figure 5** correspond very well to the analytical solution.

#### **Figure 5.**

*Crack propagation of the mode-I tension test with x*ð Þ , *<sup>y</sup>* <sup>∈</sup> ½ �� 0, 1 ½ � 0, 1 *<sup>m</sup><sup>2</sup> 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 σ<sup>y</sup> [MPa] is plotted.*

**Figure 6.** *Normalized crack growth versus the critical cohesive stress.*

They disappear at the crack flanks and are maximal at the crack tip with the typical butterfly shape.

*<sup>γ</sup>*ð Þ¼ *<sup>s</sup>*, **<sup>∇</sup>***s*,**<sup>∇</sup>** � **<sup>∇</sup>***<sup>s</sup>* ð Þ <sup>1</sup> � *<sup>s</sup>*

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

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

discretization ansatz given in Eq. (24) requires <sup>C</sup><sup>1</sup>

ð Þþ *u* G*cγ*ð Þ*s* dΩ ¼

function, which is here of the form *g s*ðÞ¼ *s*

be achieved.

*W*ð Þ¼ *u*, *s*

ð Ω Ψ*e*

general Allen-Cahn form,

works [42, 43].

**13**

**4.1 Elastic strain energy split**

reads

**Figure 7.**

crack.

2 <sup>4</sup>*<sup>ε</sup>* <sup>þ</sup> *<sup>ε</sup>*j j **<sup>∇</sup>***<sup>s</sup>*

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

which has been used, for example, in [38, 39]. However, in a finite element

[40] we investigated this formulation in more detail, applied quadratic NURBSansatz functions, and found that a smoother crack and a better convergence rate can

> ð Ω 1

The tensor <sup>∗</sup> in the elastic strain energy density is derived from Eq. (2) by means of a substitute-material approach, <sup>∗</sup> <sup>¼</sup> *g s*ð Þ, where *g s*ð Þ is a degradation

eter introduced for numerical reasons only. With *<sup>s</sup>* <sup>¼</sup> 0 tensor <sup>∗</sup> models the *empty*

where the nonnegative function *M* is the mobility with unit [mm<sup>2</sup>*=*Ns]; *δ<sup>s</sup>* denotes 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

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

*<sup>s</sup>*\_ ¼ �*Mδs*<sup>Ψ</sup> ¼ �*<sup>M</sup> <sup>ε</sup>* : : *<sup>ε</sup><sup>s</sup>* � <sup>G</sup>*<sup>c</sup>*

At this point the evolution equation for the phase-field parameter *s* is stated in a

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

<sup>2</sup> *<sup>ε</sup>* : <sup>∗</sup> : *<sup>ε</sup>* <sup>þ</sup> <sup>G</sup>*<sup>c</sup>*

2

<sup>2</sup> <sup>þ</sup> *<sup>ε</sup>*3ð Þ **<sup>∇</sup>** � **<sup>∇</sup>***<sup>s</sup>*

2


ð Þ 1 � *s* 2 <sup>4</sup>*<sup>ε</sup>* <sup>þ</sup> *<sup>ε</sup>*j j **<sup>∇</sup>***<sup>s</sup>*

<sup>2</sup> <sup>þ</sup> *<sup>s</sup>*0, with *<sup>s</sup>*<sup>0</sup> <sup>≪</sup> 1 being a small param-

*s*\_ ¼ �*Mδs*Ψ*:* (26)

<sup>2</sup>*<sup>ε</sup>* <sup>1</sup> � *<sup>s</sup>* <sup>þ</sup> <sup>4</sup>*ε*<sup>2</sup> **<sup>∇</sup>** � **<sup>∇</sup>***<sup>s</sup>* � � � �*:* (27)

!

<sup>4</sup> (24)

2

dΩ

(25)

At next we apply the traction *ty* ¼ 15MPa as a linearly increasing load within 200 load steps. A cohesive element will be added if the effective tensile stress exceeds *σ<sup>i</sup>* ¼ *σc*. The crack length is determined via *a* ≔ ð Þ max f g *x*∈½ � 0, 1 : *δn*ð Þ *x* >*δ<sup>c</sup>* � 0*:*2 . **Figure 6** shows the crack growth normalized to the crack growth at *σ<sup>c</sup>* ¼ 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 *σ<sup>i</sup>* should correspond to the cohesive strength *σ<sup>c</sup>* to avoid artificial numerical stiffness.

### **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*ð Þ *x*, *t* ∈½ � 0, 1 , ∀*x*∈ Ω, *t*∈ ℝ. 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 Γ*C*ð Þ*t* is replaced by a regularizing crack-surface density function *γ*ð Þ*t* 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 secondorder phase-field approach is considered (like originally proposed in [36, 37]) it has the form

$$\chi(\mathfrak{s}, \nabla \mathfrak{s}) = \frac{\left(\mathfrak{1} - \mathfrak{s}\right)^2}{4\varepsilon} + \mathfrak{e} \left|\nabla \mathfrak{s}\right|^2. \tag{23}$$

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

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

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

$$\chi(\mathbf{s}, \nabla \mathbf{s}, \nabla \cdot \nabla \mathbf{s}) = \frac{\left(\mathbf{1} - \mathbf{s}\right)^2}{4\varepsilon} + \frac{\varepsilon \left|\nabla \mathbf{s}\right|^2}{2} + \frac{\varepsilon^3 \left(\nabla \cdot \nabla \mathbf{s}\right)^2}{4} \tag{24}$$

which has been used, for example, in [38, 39]. However, in a finite element discretization ansatz given in Eq. (24) requires <sup>C</sup><sup>1</sup> -continuous basis functions. In [40] we investigated this formulation in more detail, applied quadratic NURBSansatz 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

$$\mathcal{W}(\mathfrak{u}, \mathfrak{s}) = \int\_{\mathfrak{U}} \Psi^{\varepsilon}(\mathfrak{u}) + \mathcal{G}\_{\varepsilon} \boldsymbol{\chi}(\mathfrak{s}) \, \mathrm{d}\mathfrak{L} = \int\_{\mathfrak{U}} \frac{1}{2} \left. \mathfrak{e} : \mathbb{C}^{\*} : \mathfrak{e} + \mathcal{G}\_{\varepsilon} \left( \frac{\left(\mathfrak{1} - \mathfrak{s}\right)^{2}}{4\varepsilon} + \varepsilon |\nabla \mathfrak{s}|^{2} \right) \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm{d}\mathfrak{L} \, \mathrm$$

The tensor <sup>∗</sup> in the elastic strain energy density is derived from Eq. (2) by means of a substitute-material approach, <sup>∗</sup> <sup>¼</sup> *g s*ð Þ, where *g s*ð Þ is a degradation function, which is here of the form *g s*ðÞ¼ *s* <sup>2</sup> <sup>þ</sup> *<sup>s</sup>*0, with *<sup>s</sup>*<sup>0</sup> <sup>≪</sup> 1 being a small parameter introduced for numerical reasons only. With *<sup>s</sup>* <sup>¼</sup> 0 tensor <sup>∗</sup> 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\_\circ \Psi. \tag{26}
$$

where the nonnegative function *M* is the mobility with unit [mm<sup>2</sup>*=*Ns]; *δ<sup>s</sup>* denotes the variational or total derivative of Ψ wrt. *s*. With Eq. (25) it follows

$$\dot{\mathcal{S}} = -\overline{\mathcal{M}} \delta\_{\mathfrak{s}} \Psi = -\overline{\mathcal{M}} \left[ \mathbf{e} : \mathbb{C} : \mathfrak{e} \, \mathfrak{s} - \frac{\mathcal{G}\_{\mathfrak{c}}}{2\varepsilon} \left( \mathbb{1} - \mathfrak{s} + 4\mathfrak{e}^{2} \, \nabla \cdot \nabla \mathfrak{s} \right) \right]. \tag{27}$$

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

They disappear at the crack flanks and are maximal at the crack tip with the typical

At next we apply the traction *ty* ¼ 15MPa as a linearly increasing load within 200 load steps. A cohesive element will be added if the effective tensile stress exceeds *σ<sup>i</sup>* ¼ *σc*. The crack length is determined via *a* ≔ ð Þ max f g *x*∈½ � 0, 1 : *δn*ð Þ *x* >*δ<sup>c</sup>* � 0*:*2 . **Figure 6** shows the crack growth normalized to the crack growth at *σ<sup>c</sup>* ¼ 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 *σ<sup>i</sup>* should correspond to the cohesive strength *σ<sup>c</sup>* to avoid artificial

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*ð Þ *x*, *t* ∈½ � 0, 1 , ∀*x*∈ Ω, *t*∈ ℝ. 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 Γ*C*ð Þ*t* is replaced by a regularizing crack-surface density function *γ*ð Þ*t* and the corresponding boundary

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

> 2 <sup>4</sup>*<sup>ε</sup>* <sup>þ</sup> *<sup>ε</sup>*j j **<sup>∇</sup>***<sup>s</sup>*

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

2

*:* (23)

*<sup>γ</sup>*ð Þ¼ *<sup>s</sup>*,**∇***<sup>s</sup>* ð Þ <sup>1</sup> � *<sup>s</sup>*

butterfly shape.

**Figure 6.**

numerical stiffness.

the form

fourth-order form

**12**

**4. The phase-field fracture approach**

*Normalized crack growth versus the critical cohesive stress.*

*Modeling and Simulation in Engineering - Selected Problems*

integral in Eq. (5) is approximated by Eq. (1).

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 s*ð Þ and state <sup>Ψ</sup>*<sup>e</sup>* ð Þ¼ *<sup>ε</sup>*, *<sup>s</sup> g s*ð ÞΨ*e*<sup>þ</sup> <sup>þ</sup> <sup>Ψ</sup>*e*�. This split is based on the spectral decomposition of the strain tensor *<sup>ε</sup>* <sup>¼</sup> <sup>P</sup><sup>3</sup> *<sup>a</sup>*¼1*εan<sup>a</sup>* <sup>⊗</sup> *<sup>n</sup>a*, where *ε<sup>a</sup>* denote the principal strains and *n<sup>a</sup>* the corresponding principal directions, *a* ¼ 1, 2, 3. Based on this representations and using the Macaulay brackets

**A** ¼ ⋃ *ne e*¼1 4*ε*<sup>2</sup> ð Ω

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

equations is obtained:

**<sup>∇</sup>N**ð Þ **<sup>∇</sup><sup>N</sup>** *<sup>T</sup>* <sup>d</sup>Ω, **<sup>D</sup>** <sup>¼</sup> <sup>⋃</sup>

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

**<sup>M</sup>**ð Þ <sup>Δ</sup>*<sup>t</sup>* �<sup>1</sup> <sup>þ</sup> **<sup>A</sup>** <sup>þ</sup> **<sup>C</sup>** <sup>þ</sup> **<sup>D</sup>**

explicit time discretization, which simplifies the solution.

linearly increasing displacement up to *<sup>u</sup>* ¼ �1*:*<sup>2</sup> � <sup>10</sup>�3.

**4.3 Effect of the coefficients** *M* **and** *ε*

only limited by the numerics.

**Figure 8.**

**15**

Here we assumed the same ansatz functions for *u* and *s* for the ease of writing. After the discretization in time by using Eq. (12) the following system of global

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

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 0, 1 ½ �� ½ � 0, 1 , with a predefined crack at *x*∈ ½ � 0, 0*:*5 , *y* ¼ 0*:*5, and boundary conditions as depicted in **Figure 8**. The finite element mesh of 100 � 100 � 4 triangular elements has a uniform size of *h* ¼ 0*:*005. The simulation is quasistatic with a prescribed,

The kinematic mobility parameter *M* ¼ *M*G*c=*ð Þ 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 *<sup>M</sup>* <sup>¼</sup> 1 is demonstrated. We see a steady state for *<sup>M</sup>* <sup>≥</sup><sup>50</sup> mm2sN. In this case the crack length is defined as *a* ≔ max f g *x*∈½ � 0, 1 : *s x*ð Þ , 0*:*5 ≤0*:*05 � 0*:*5*:* The mobility *M* [s] can be seen as the inverse of a kinematic viscosity *η* [s�<sup>1</sup>

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

*Mode-I tension-test with x*ð Þ , *y* ∈ ½ �� 0, 1 ½ � 0, 1 *, an initial crack of length 0.5, and a prescribed total displacement of <sup>u</sup>* ¼ �1*:*<sup>2</sup> � <sup>10</sup>�<sup>3</sup>*. Plotted is the normalized crack length* <sup>a</sup> *versus the kinematic mobility* <sup>M</sup>*.*

h i**s***<sup>n</sup>*þ<sup>1</sup> <sup>¼</sup> **<sup>f</sup>** <sup>þ</sup> ð Þ <sup>Δ</sup>*<sup>t</sup>* �<sup>1</sup>

*ne e*¼1 4*ε* G*c* ð Ω

<sup>Ψ</sup>*<sup>e</sup>*þð Þ *<sup>u</sup>*, *<sup>s</sup>* **NN***<sup>T</sup>* <sup>d</sup>Ω*:*

**Ms***<sup>n</sup>* (30)

], *M* ¼

$$\left\langle \mathbf{x} \right\rangle^{+} = \begin{cases} \mathbf{0}, & \mathbf{x} < \mathbf{0} \\ \mathbf{x}, & \mathbf{x} \ge \mathbf{0}, \end{cases} \qquad \left\langle \mathbf{x} \right\rangle^{-} = \begin{cases} \mathbf{x}, & \mathbf{x} < \mathbf{0} \\ \mathbf{0}, & \mathbf{x} \ge \mathbf{0}, \end{cases}$$

we define the positive and the negative parts of the strain tensor as *<sup>ε</sup>*� <sup>¼</sup> <sup>P</sup><sup>3</sup> *<sup>a</sup>*¼1h i *<sup>ε</sup><sup>a</sup>* � *n<sup>a</sup>* ⊗ *na*. 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, <sup>Ψ</sup>*<sup>e</sup>*<sup>þ</sup> <sup>¼</sup> <sup>1</sup>*=*2h i *<sup>σ</sup>* <sup>þ</sup> : h i*<sup>ε</sup>* <sup>þ</sup> <sup>¼</sup> <sup>1</sup>*=*2h i*<sup>ε</sup>* <sup>þ</sup> : <sup>∗</sup> : h i*<sup>ε</sup>* <sup>þ</sup>.

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 *s*\_ ¼ 0 if *s* ¼ 0. The same effect has a product of *s* with the mobility parameter *M* ¼ *M*G*c=*2*ε*, where we additionally formulate the evolution Eq. (26) in a dimensionless form, *<sup>s</sup>*\_ ¼ �*Mδs*<sup>Ψ</sup> ¼ �*<sup>M</sup>* <sup>Ψ</sup>*<sup>e</sup>*<sup>þ</sup> � ð Þ� <sup>1</sup> � *<sup>s</sup>* <sup>4</sup>*ε*<sup>2</sup> ½ � **<sup>∇</sup>** � **<sup>∇</sup>***<sup>s</sup>* .

### **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\_{\mathfrak{U}} \left( \frac{\dot{\mathfrak{s}}}{\overline{M}} \delta \mathfrak{s} + \partial\_{\mathfrak{s}} \Psi^{e+} \delta \mathfrak{s} + 4 \mathcal{G}\_{\mathfrak{c}} \varepsilon \nabla \mathfrak{s} \nabla \delta \mathfrak{s} - \frac{\mathcal{G}\_{\mathfrak{c}}}{2 \mathfrak{c}} (\mathfrak{1} - \mathfrak{s}) \delta \mathfrak{s} \right) d \mathfrak{Q} = \mathbf{0} \quad \forall \delta \mathfrak{s} \in \mathfrak{V} \tag{28}$$

with test functions *<sup>δ</sup><sup>s</sup>* � *<sup>v</sup>* and <sup>V</sup> <sup>¼</sup> *<sup>v</sup>*<sup>∈</sup> *<sup>H</sup>*<sup>1</sup> ð Þj Ω **∇***v* � *n* ¼ 0 on ∂Ω, *v* ¼ 0 on Γ<sup>1</sup> � �.

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

$$\boldsymbol{s} \approx \sum\_{i=1}^{n\_k} N\_i \boldsymbol{s}\_i = \mathbf{N} \hat{\mathbf{s}}, \qquad \boldsymbol{\delta} \approx \sum\_{i=1}^{n\_k} N\_i \boldsymbol{v}\_i = \mathbf{N} \hat{\mathbf{v}} \tag{29}$$

with ansatz functions *N*1, … , *Nnk* for *nk* 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 *<sup>∂</sup>s*Ψ*<sup>e</sup>* <sup>¼</sup> <sup>2</sup>*s*Ψ*<sup>e</sup>*<sup>þ</sup> <sup>¼</sup> *<sup>s</sup>*h i *<sup>ε</sup>*ð Þ *<sup>u</sup>* <sup>þ</sup> : : h i *ε*ð Þ *u* <sup>þ</sup>.

In short hand notation we get for the mechanical problem the matrix Eq. (11), whereby now the Hookean matrix <sup>∗</sup> ð Þ *<sup>u</sup>*, *<sup>s</sup>* depends on the phase-field. For Eq. (28) the approximation leads to **<sup>M</sup>**s\_ <sup>þ</sup> **<sup>A</sup>** <sup>þ</sup> **<sup>C</sup>** <sup>þ</sup> **<sup>D</sup>** � �**<sup>s</sup>** <sup>¼</sup> **<sup>f</sup>** with **<sup>s</sup>** <sup>¼</sup> *si* ½ � and

$$\overline{\mathbf{M}} = \bigcup\_{\epsilon=1}^{n\_{\epsilon}} \frac{1}{M} \int\_{\Omega} \mathbf{N} \mathbf{N}^{T} \, \mathrm{d}\Omega,\qquad \overline{\mathbf{C}} = \bigcup\_{\epsilon=1}^{n\_{\epsilon}} \int\_{\Omega} \mathbf{N} \mathbf{N}^{T} \, \mathrm{d}\Omega,\qquad \overline{\mathbf{f}} = \bigcup\_{\epsilon=1}^{n\_{\epsilon}} \int\_{\Omega} \mathbf{N} \mathbf{N}^{T} \mathbf{1} \, \mathrm{d}\Omega.$$

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

$$\overline{\mathbf{A}} = \bigcup\_{\epsilon=1}^{n\_{\epsilon}} 4\epsilon^{2} \int\_{\Omega} \mathbf{V} \mathbf{N} (\mathbf{V} \mathbf{N})^{T} \, \mathrm{d}\Omega,\qquad \overline{\mathbf{D}} = \bigcup\_{\epsilon=1}^{n\_{\epsilon}} \frac{4\epsilon}{\mathcal{G}\_{\epsilon}} \int\_{\Omega} \Psi^{\epsilon+}(\mathbf{u}, \boldsymbol{\epsilon}) \mathbf{N} \mathbf{N}^{T} \, \mathrm{d}\Omega.$$

Here we assumed the same ansatz functions for *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}}(\Delta t)^{-1} + \overline{\mathbf{A}} + \overline{\mathbf{C}} + \overline{\mathbf{D}}\right] \mathbf{s}^{n+1} = \overline{\mathbf{f}} + (\Delta t)^{-1} \overline{\mathbf{M}} \mathbf{s}^{n} \tag{30}$$

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

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.

*ε<sup>a</sup>* denote the principal strains and *n<sup>a</sup>* the corresponding principal directions, *a* ¼

h i *<sup>x</sup>* � <sup>¼</sup> *<sup>x</sup>*, *<sup>x</sup>*<<sup>0</sup>

�

0, *x*≥0,

ð Þ¼ *<sup>ε</sup>*, *<sup>s</sup> g s*ð ÞΨ*e*<sup>þ</sup> <sup>þ</sup> <sup>Ψ</sup>*e*�. This split is

dΩ ¼ 0 ∀*δs*∈ V (28)

*Nivi* ¼ **Nv**^ (29)

ð Þj Ω **∇***v* � *n* ¼ 0 on ∂Ω, *v* ¼ 0 on Γ<sup>1</sup> � �.

> *ne e*¼1

ð Ω

**NN***<sup>T</sup>***1** dΩ

*<sup>a</sup>*¼1*εan<sup>a</sup>* <sup>⊗</sup> *<sup>n</sup>a*, where

We use the degradation function *g s*ð Þ and state <sup>Ψ</sup>*<sup>e</sup>*

*Modeling and Simulation in Engineering - Selected Problems*

based on the spectral decomposition of the strain tensor *<sup>ε</sup>* <sup>¼</sup> <sup>P</sup><sup>3</sup>

h i *<sup>x</sup>* <sup>þ</sup> <sup>¼</sup> 0, *<sup>x</sup>*<sup>&</sup>lt; <sup>0</sup>

�

*<sup>s</sup>*\_ ¼ �*Mδs*<sup>Ψ</sup> ¼ �*<sup>M</sup>* <sup>Ψ</sup>*<sup>e</sup>*<sup>þ</sup> � ð Þ� <sup>1</sup> � *<sup>s</sup>* <sup>4</sup>*ε*<sup>2</sup> ½ � **<sup>∇</sup>** � **<sup>∇</sup>***<sup>s</sup>* .

*<sup>a</sup>*¼1h i *<sup>ε</sup><sup>a</sup>* �

**4.2 Discretization**

ð Ω

h i *ε*ð Þ *u* <sup>þ</sup>.

**14**

**M** ¼ ⋃ *ne e*¼1

1 *M* ð Ω

*s*\_ *M*

1, 2, 3. Based on this representations and using the Macaulay brackets

*x*, *x*≥0,

we define the positive and the negative parts of the strain tensor as *<sup>ε</sup>*� <sup>¼</sup> <sup>P</sup><sup>3</sup>

only accounts for tension, <sup>Ψ</sup>*<sup>e</sup>*<sup>þ</sup> <sup>¼</sup> <sup>1</sup>*=*2h i *<sup>σ</sup>* <sup>þ</sup> : h i*<sup>ε</sup>* <sup>þ</sup> <sup>¼</sup> <sup>1</sup>*=*2h i*<sup>ε</sup>* <sup>þ</sup> : <sup>∗</sup> : h i*<sup>ε</sup>* <sup>þ</sup>.

weak formulation of the phase-field equation is set up analogously

*<sup>δ</sup><sup>s</sup>* <sup>þ</sup> *<sup>∂</sup>s*Ψ*<sup>e</sup>*þ*δ<sup>s</sup>* <sup>þ</sup> <sup>4</sup>G*cε***∇***s***∇***δ<sup>s</sup>* � <sup>G</sup>*<sup>c</sup>*

with test functions *<sup>δ</sup><sup>s</sup>* � *<sup>v</sup>* and <sup>V</sup> <sup>¼</sup> *<sup>v</sup>*<sup>∈</sup> *<sup>H</sup>*<sup>1</sup>

*<sup>s</sup>* <sup>≈</sup> <sup>X</sup>*nk i*¼1

� �

tion 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

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 *s*\_ ¼ 0 if *s* ¼ 0. The same effect has a product of *s* with the mobility parameter *M* ¼ *M*G*c=*2*ε*, where we additionally formulate the evolution Eq. (26) in a dimensionless form,

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

<sup>2</sup>*<sup>ε</sup>* ð Þ <sup>1</sup> � *<sup>s</sup> <sup>δ</sup><sup>s</sup>*

*i*¼1

**NN***<sup>T</sup>* <sup>d</sup>Ω, **<sup>f</sup>** <sup>¼</sup> <sup>⋃</sup>

The shape and test functions given in Eq. (10) are inserted in Eq. (9) to obtain the

with ansatz functions *N*1, … , *Nnk* for *nk* 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 *<sup>∂</sup>s*Ψ*<sup>e</sup>* <sup>¼</sup> <sup>2</sup>*s*Ψ*<sup>e</sup>*<sup>þ</sup> <sup>¼</sup> *<sup>s</sup>*h i *<sup>ε</sup>*ð Þ *<sup>u</sup>* <sup>þ</sup> : :

In short hand notation we get for the mechanical problem the matrix Eq. (11), whereby now the Hookean matrix <sup>∗</sup> ð Þ *<sup>u</sup>*, *<sup>s</sup>* depends on the phase-field. For Eq. (28)

> *ne e*¼1

ð Ω

the approximation leads to **<sup>M</sup>**s\_ <sup>þ</sup> **<sup>A</sup>** <sup>þ</sup> **<sup>C</sup>** <sup>þ</sup> **<sup>D</sup>** � �**<sup>s</sup>** <sup>¼</sup> **<sup>f</sup>** with **<sup>s</sup>** <sup>¼</sup> *si* ½ � and

**NN***<sup>T</sup>* <sup>d</sup>Ω, **<sup>C</sup>** <sup>¼</sup> <sup>⋃</sup>

discrete mechanical problem of Eq. (11). For the phase-field we approximate

*Nisi* <sup>¼</sup> **<sup>N</sup>**^**s**, *<sup>δ</sup>s*<sup>≈</sup> <sup>X</sup>*nk*

*n<sup>a</sup>* ⊗ *na*. The positive parts contain contributions due to positive dilata-

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 0, 1 ½ �� ½ � 0, 1 , with a predefined crack at *x*∈ ½ � 0, 0*:*5 , *y* ¼ 0*:*5, and boundary conditions as depicted in **Figure 8**. The finite element mesh of 100 � 100 � 4 triangular elements has a uniform size of *h* ¼ 0*:*005. The simulation is quasistatic with a prescribed, linearly increasing displacement up to *<sup>u</sup>* ¼ �1*:*<sup>2</sup> � <sup>10</sup>�3.

The kinematic mobility parameter *M* ¼ *M*G*c=*ð Þ 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 *<sup>M</sup>* <sup>¼</sup> 1 is demonstrated. We see a steady state for *<sup>M</sup>* <sup>≥</sup><sup>50</sup> mm2sN. In this case the crack length is defined as *a* ≔ max f g *x*∈½ � 0, 1 : *s x*ð Þ , 0*:*5 ≤0*:*05 � 0*:*5*:* The mobility *M* [s] can be seen as the inverse of a kinematic viscosity *η* [s�<sup>1</sup> ], *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.

#### **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 <sup>u</sup>* ¼ �1*:*<sup>2</sup> � <sup>10</sup>�<sup>3</sup>*. Plotted is the normalized crack length* <sup>a</sup> *versus the kinematic mobility* <sup>M</sup>*.*

**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 <sup>R</sup>*<sup>2</sup> <sup>¼</sup> <sup>0</sup>*:*9986*, and the right plot demonstrates the phase-field* <sup>s</sup> *in the diffuse interface zone to the crack in corresponding colors.*

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

$$b \coloneqq \left| \{ y \in [0, 1] : s(0.9, y) \le 0.95 \} \right|\_{\epsilon}, \tag{31}$$

transmitted into the specimen. This stress pulse *σ<sup>T</sup>* 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

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

The aim of spalling experiments is to determine a material's resistance to fracture, specifically its fracture energy *Wc* or specific energy G*<sup>c</sup>* and/or its tensile

pression but fails under tension. Plastic deformations do not matter, neither do

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

A first challenge was the correct reproduction of the incident and reflected stress

*<sup>m</sup>* is usually defined as the

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* ¼ *σ*max*f t*ð Þ with *f*ð Þ¼ 0≤ *t*≤*t*<sup>1</sup> *t=t*1, *f t*ð Þ¼ <sup>1</sup> <*t*≤ *t*<sup>2</sup> 1, and, *f t*ð Þ¼ <sup>2</sup> <*t*≤*t*<sup>3</sup> *t*<sup>3</sup> � *t=t*<sup>3</sup> � *t*2, where *t*<sup>1</sup> ¼ 12μs and *t*<sup>2</sup> ¼ 47 μs and *t*<sup>3</sup> ¼ 62μs and *σ*max ¼ �17MPa. Within the 80μ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*spec ¼ 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

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

*Schematics of the HB-spalling experiment and typical incoming stress wave in the middle of the specimen with*

temperature effects; all tests are conducted under ambient conditions.

*m*. The spalling test presumes brittle materials which can sustain com-

specimen of **Figure 1**.

**5.1 Finite element discretization**

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

numerical details we refer to [34, 45].

The dynamic tensile resistance of a brittle material *Rt*

**5.2 Fracture parameter**

**Figure 10.**

**17**

*t*<sup>1</sup> ¼ 0 *and t*<sup>2</sup> ¼ *t*<sup>3</sup> ¼ 30 μs*.*

resistance *Rt*

studies.

whereby the Lebesgue-measure of a set is denoted with j j� *<sup>e</sup>* and 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 G*c=ε* competes with the elastic energy density Ψ*<sup>e</sup>*<sup>þ</sup> determined by *E*, *ν* or combinations thereof. Here *ε* has the effect of a material parameter. On the other hand, *ε* 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, *ε*>*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 ¼ *cε* 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 *ε* ¼ 1 … 3*h* for linear shape functions, which is similar to other author's suggestion of *lc* ¼ 2*ε* >2*h*, [12].

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

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

transmitted into the specimen. This stress pulse *σ<sup>T</sup>* 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**.

The aim of spalling experiments is to determine a material's resistance to fracture, specifically its fracture energy *Wc* or specific energy G*<sup>c</sup>* and/or its tensile resistance *Rt m*. 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 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

(31)

*b* ≔ f*y*∈ ½ � 0, 1 : *s*ð0*:*9, *y*Þ≤0*:*95gj *<sup>e</sup>*,

*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 <sup>R</sup>*<sup>2</sup> <sup>¼</sup> <sup>0</sup>*:*9986*, and the right plot demonstrates the phase-field* <sup>s</sup> *in the diffuse interface zone to the crack in*

whereby the Lebesgue-measure of a set is denoted with j j� *<sup>e</sup>* and the diffuse

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 G*c=ε* competes with the elastic energy density Ψ*<sup>e</sup>*<sup>þ</sup> determined by *E*, *ν* or combinations thereof. Here *ε* has the effect of a material parameter. On the other hand, *ε* is determined by the mesh size *h* because it has to be large enough to enable the approximation of a diffuse interface. This clearly

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 ¼ *cε* 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 *ε* ¼ 1 … 3*h* for linear shape functions, which is similar to other author's

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

called crack-free. Hence the width of the crack is defined as

*Modeling and Simulation in Engineering - Selected Problems*

interface zone is included in the relevant set.

**Figure 9.**

*corresponding colors.*

requires an adaption on the mesh size, *ε*>*h*.

**5. Simulation of the HB-spalling experiment**

suggestion of *lc* ¼ 2*ε* >2*h*, [12].

**16**

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* ¼ *σ*max*f t*ð Þ with *f*ð Þ¼ 0≤ *t*≤*t*<sup>1</sup> *t=t*1, *f t*ð Þ¼ <sup>1</sup> <*t*≤ *t*<sup>2</sup> 1, and, *f t*ð Þ¼ <sup>2</sup> <*t*≤*t*<sup>3</sup> *t*<sup>3</sup> � *t=t*<sup>3</sup> � *t*2, where *t*<sup>1</sup> ¼ 12μs and *t*<sup>2</sup> ¼ 47 μs and *t*<sup>3</sup> ¼ 62μs and *σ*max ¼ �17MPa. Within the 80μ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*spec ¼ 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 *Rt <sup>m</sup>* 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

#### **Figure 10.**

*Schematics of the HB-spalling experiment and typical incoming stress wave in the middle of the specimen with t*<sup>1</sup> ¼ 0 *and t*<sup>2</sup> ¼ *t*<sup>3</sup> ¼ 30 μs*.*

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<sup>t</sup> <sup>m</sup>* with the cohesive strength *σ<sup>c</sup>* 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 *δ<sup>c</sup>* ¼ 10μm for *σ<sup>c</sup>* ¼ 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 *σ<sup>c</sup>* ¼ 20MPa, the cohesive stress is obviously too high, no elements are inserted and no crack can grow. For *σ<sup>c</sup>* ¼ 17MPa a small, localized crack zone develops. Lowering *σ<sup>c</sup>* ¼ 14 … 8MPa gives a wider zone and, moreover, the cracked zone moves toward the free end. This follows from the fact that the superposed pulse *σ<sup>R</sup>* þ *σ<sup>I</sup>* reaches the value of *σ<sup>c</sup>* earlier. However, this position does not correspond to the measured crack position. In all our experiments the crack appeared at *zc* ¼ 120 … 130mm, with coordinate *z* starting at the free end.

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 *<sup>ε</sup>* <sup>¼</sup> <sup>3</sup>*<sup>h</sup>* and *<sup>M</sup>*<sup>~</sup> <sup>¼</sup> 200, 000 s�<sup>1</sup> and define regions with phase-field parameter

*Influence of the critical crack opening displacement on the time of crack initiation* t*<sup>1</sup> and of total crack opening* t*<sup>2</sup> for σ<sup>c</sup>* ¼ 15*MPa with (left) δ<sup>c</sup> in the linear cohesive law Eq. (17) and (right) δ*<sup>0</sup> *in the exponential form*

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

For a large critical energy release rate, G*<sup>c</sup>* >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 G*<sup>c</sup>* ¼ 90N*=*m a sharp crack zone appears at the expected position. Smaller values, for example G*<sup>c</sup>* ¼ 30N*=*m in **Figure 13(left)**, show a wider cracked zone but the position of the crack is nearly the same, see

*b*trans ≔ f*z*∈½ � 0, 200 : 0*:*2 ≤*s*ð0, *z*Þ≤0*:*8gj

**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 G*c*. Similar is the situation for an decrease of parameter *ε*. A value of *ε*< 6*h* increases the time of crack formation at

Further studies have been performed to simulate the experiments, cf. [34], and with the collected knowledge we conclude, that for UHPC the tensile strength *Rt*

*Influence of the specific fracture energy* G*<sup>c</sup> on the width given in Eq. (32) of the crack; left plot shows the phasefield parameter at crack position in longitudinal direction and right plot shows width b as function of* G*c,*

*e*

*m*,

(32)

**Figure 13(right)**. Specifically we define the diffuse zone as

and the cracked zone is given by *s*ð Þ 0, *z* ≤0*:*8.

*s*<0*:*2 as cracked.

*with δ<sup>c</sup>* ¼ 30μm *Eq. (19).*

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

**Figure 12.**

constant G*c*.

**Figure 13.**

**19**

*approximated with* R*-coefficient R*<sup>2</sup> ≈0*:*9894*.*

Obviously, the cohesive stress *σ<sup>c</sup>* 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 *δc*. Here we observe for higher values of *δ<sup>c</sup>* a longer time of crack opening whereas *δ<sup>c</sup>* 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 *σ<sup>c</sup>* and *δ*0, cf. **Figure 12**.

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


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

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

**Figure 12.**

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<sup>t</sup>*

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 *δ<sup>c</sup>* ¼ 10μm for *σ<sup>c</sup>* ¼ 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

**Figure 11** shows one symmetry half of the specimen at the end of the simulation for different values of *σc*. At first, with *σ<sup>c</sup>* ¼ 20MPa, the cohesive stress is obviously too high, no elements are inserted and no crack can grow. For *σ<sup>c</sup>* ¼ 17MPa a small, localized crack zone develops. Lowering *σ<sup>c</sup>* ¼ 14 … 8MPa gives a wider zone and, moreover, the cracked zone moves toward the free end. This follows from the fact that the superposed pulse *σ<sup>R</sup>* þ *σ<sup>I</sup>* reaches the value of *σ<sup>c</sup>* earlier. However, this position does not correspond to the measured crack position. In all our experiments the crack appeared at *zc* ¼ 120 … 130mm, with coordinate *z* starting at the free end. Obviously, the cohesive stress *σ<sup>c</sup>* 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

We further studied the influence of the critical crack opening displacement *δc*. Here we observe for higher values of *δ<sup>c</sup>* a longer time of crack opening whereas *δ<sup>c</sup>* 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 *σ<sup>c</sup>* and *δ*0, cf.

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

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

the cohesive strength *σ<sup>c</sup>* and try to capture it by adaptive insertion of cohesive

elements.

**Figure 12**.

**Figure 11.**

**18**

*purple, the cracks (open elements) are black.*

*5.2.2 Phase-field fracture*

*5.2.1 Cohesive element technique*

after a spallation plane has built or after 120 μs.

*Modeling and Simulation in Engineering - Selected Problems*

cannot appear at the right position anymore.

*<sup>m</sup>* with

*Influence of the critical crack opening displacement on the time of crack initiation* t*<sup>1</sup> and of total crack opening* t*<sup>2</sup> for σ<sup>c</sup>* ¼ 15*MPa with (left) δ<sup>c</sup> in the linear cohesive law Eq. (17) and (right) δ*<sup>0</sup> *in the exponential form with δ<sup>c</sup>* ¼ 30μm *Eq. (19).*

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 *<sup>ε</sup>* <sup>¼</sup> <sup>3</sup>*<sup>h</sup>* and *<sup>M</sup>*<sup>~</sup> <sup>¼</sup> 200, 000 s�<sup>1</sup> and define regions with phase-field parameter *s*<0*:*2 as cracked.

For a large critical energy release rate, G*<sup>c</sup>* >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 G*<sup>c</sup>* ¼ 90N*=*m a sharp crack zone appears at the expected position. Smaller values, for example G*<sup>c</sup>* ¼ 30N*=*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

$$b\_{\text{trans}} \coloneqq \left| \{ z \in [0, 200] : 0.2 \le s(0, z) \le 0.8 \} \right|\_{\epsilon} \tag{32}$$

and the cracked zone is given by *s*ð Þ 0, *z* ≤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 G*c*. Similar is the situation for an decrease of parameter *ε*. A value of *ε*< 6*h* increases the time of crack formation at constant G*c*.

Further studies have been performed to simulate the experiments, cf. [34], and with the collected knowledge we conclude, that for UHPC the tensile strength *Rt m*,

#### **Figure 13.**

*Influence of the specific fracture energy* G*<sup>c</sup> on the width given in Eq. (32) of the crack; left plot shows the phasefield parameter at crack position in longitudinal direction and right plot shows width b as function of* G*c, approximated with* R*-coefficient R*<sup>2</sup> ≈0*:*9894*.*

**Figure 14.**

*Influence of the specific fracture energy* G*<sup>c</sup> is shown in left plot and the length parameter ε on the time of crack initiation* <sup>t</sup>*<sup>1</sup> and of total crack opening* <sup>t</sup>*<sup>2</sup> is shown in the right plot, computed with M* <sup>¼</sup> <sup>0</sup>*:*2μs�<sup>1</sup> *and <sup>ε</sup>* <sup>¼</sup> <sup>3</sup>*<sup>h</sup> (left), Gc* ¼ 90N*=*m *(right).*

For the inverse analysis we proceed as follows: we use the data obtained in our

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

experiment, and from the computed velocities and masses we deduce via Eq. (33) the

We mesh the specimen uniformly with 2560 elements and employ the linear cohesive law from Eq. (17) with parameter *β* ¼ 1. The cohesive stress is *σ<sup>c</sup>* ¼ 15MPa to ensure crack formation at the expected position and we vary the critical opening displacement *δ<sup>c</sup>* to obtain values of G*<sup>c</sup>* between 20 N/m and 150 N/m. All computations result in similar fragments with masses *m*fra,1 ∈ ½ � 895, 905 g and *m*fra,2 ∈½ � 572, 582 g.

validates the chosen method of energy balance. However, the value of G*<sup>c</sup>* 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

*<sup>c</sup>* with discrepancies between 30 and 60%, see **Figure 16**.

and set *<sup>ε</sup>* <sup>¼</sup> <sup>2</sup>*:*<sup>5</sup> mm and *<sup>M</sup>* <sup>¼</sup> <sup>0</sup>*:*<sup>6</sup> <sup>μ</sup>s�1. At first we accelerate the specimen to *v*spec ¼ 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 Ginp

the easier both fragments split, cf. **Figure 14**.

For phase-field simulations we use a finer mesh of 5760 elements, *h* ¼ 5*=*6 mm

*<sup>e</sup>*ð Þ fra,1 <sup>K</sup>*e=m*fra,1 � �<sup>1</sup>*=*<sup>2</sup>

P

*<sup>c</sup>* . To this end we determine the velocities of the fragments by

, with the kinetic energy of the finite elements, <sup>K</sup>*<sup>e</sup>* <sup>¼</sup> <sup>1</sup>

*<sup>c</sup>* , this is marked with a dotted line. The computed values show a

*<sup>c</sup>* and its simulated counterpart <sup>G</sup>sim

*<sup>c</sup>* , simulate the spallation

.

*<sup>c</sup>* is shown. In the optimal

*<sup>c</sup>* . This generally

*<sup>c</sup>* , the lower it is

<sup>2</sup> *me <sup>v</sup><sup>e</sup>* j j<sup>2</sup> .

, and, *v*fra,2 ¼

previous simulations to define a range of input data Ginp

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

, *v*fra,1 ¼ 2

Also the fractured surfaces show little variations, *Ac* <sup>∈</sup> ½ � <sup>9</sup>*:*5, 11*:*<sup>5</sup> cm2

In **Figure 16** the computed specific fracture energy Gsim

"measured" value Gsim

P

*<sup>e</sup>*ð Þ fra,2 <sup>K</sup>*e=m*fra,2 � �<sup>1</sup>*=*<sup>2</sup>

*<sup>e</sup>*K*e=m*spec � �<sup>1</sup>*=*<sup>2</sup>

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

**6.2 Cohesive element technique**

*v*spec ¼ 2

**Figure 15.**

case is Ginp

Ginp *<sup>c</sup>* <sup>&</sup>lt; <sup>G</sup>sim

**21**

*<sup>c</sup>* <sup>¼</sup> <sup>G</sup>sim

**6.3 Phase-field fracture**

clear proportionality between Ginp

2 P

the critical opening displacement *δc*, and the specific fracture energy G*<sup>c</sup>* are in the ranges 12<*R<sup>t</sup> <sup>m</sup>* ½ � MPa <18, 5< *δ<sup>c</sup>* ½ � *μ*m ≤17*:*5, and 40≤ G*<sup>c</sup>* ½ � N*=*m ≤105.

### **6. Inverse analysis: determination of the** G*<sup>c</sup>*

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 *Wc* corresponds to the amount of work necessary to form such a new surface. In order to determine *Wc* we balance the energy before and after crack initiation, that is, at time *t*<sup>1</sup> right before cracking and at time *t*<sup>2</sup> immediately after the crack opened. In the simple case of two fragments, cf. **Figure 15**, with masses *m*fra,1, *m*fra,2, and *m*spec ¼ *m*fra,1 þ *m*fra,2 it follows for the difference of kinetic energy:

$$
\Delta\mathcal{K} = \frac{1}{2} m\_{\rm spec} \left( v\_{\rm spec}(t\_{\underline{2}}) \right)^2 - \frac{1}{2} m\_{\rm fra,1} \left( v\_{\rm fra,1}(t\_{\underline{2}}) \right)^2 - \frac{1}{2} m\_{\rm fra,2} \left( v\_{\rm fra,2}(t\_{\underline{2}}) \right)^2 \tag{33}
$$

where *v*fra,1, *v*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 ΔK ¼ *Wc* to be the fracture energy of the specimen.

In order to deduce G*<sup>c</sup>* the fracture energy Eq. (33) needs to be referred to the fractured surface *Ac* of the specimen. Ideally, the crack is smooth and perpendicular to the specimen's axis. Then the specific fracture energy simply follows as G*<sup>c</sup>* ¼ *Wc=Ac* with *Ac* <sup>¼</sup> <sup>2</sup>*πr*<sup>2</sup> spec. In practice, there may be spalled and rough crack surfaces and so we replaced *Ac* by the measured surface size.

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

**Figure 15.**

the critical opening displacement *δc*, and the specific fracture energy G*<sup>c</sup>* are in the

*Influence of the specific fracture energy* G*<sup>c</sup> is shown in left plot and the length parameter ε on the time of crack initiation* <sup>t</sup>*<sup>1</sup> and of total crack opening* <sup>t</sup>*<sup>2</sup> is shown in the right plot, computed with M* <sup>¼</sup> <sup>0</sup>*:*2μs�<sup>1</sup> *and <sup>ε</sup>* <sup>¼</sup> <sup>3</sup>*<sup>h</sup>*

*<sup>m</sup>* ½ � MPa <18, 5< *δ<sup>c</sup>* ½ � *μ*m ≤17*:*5, and 40≤ G*<sup>c</sup>* ½ � N*=*m ≤105.

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

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 *Wc* corresponds to the amount of work necessary to form such a new surface. In order to determine *Wc* we balance the energy before and after crack initiation, that is, at time *t*<sup>1</sup> right before cracking and at time *t*<sup>2</sup> immediately after the crack opened. In the simple case of two fragments, cf. **Figure 15**, with masses *m*fra,1, *m*fra,2, and *m*spec ¼ *m*fra,1 þ *m*fra,2 it follows for the

*m*fra,1 *v*fra,1 *t*<sup>2</sup>

where *v*fra,1, *v*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. Conse-

In order to deduce G*<sup>c</sup>* the fracture energy Eq. (33) needs to be referred to the fractured surface *Ac* of the specimen. Ideally, the crack is smooth and perpendicular to the specimen's axis. Then the specific fracture energy simply follows as G*<sup>c</sup>* ¼

<sup>2</sup> � <sup>1</sup>

2

spec. In practice, there may be spalled and rough crack surfaces

*m*fra,2 *v*fra,2 *t*<sup>2</sup>

<sup>2</sup> (33)

**6. Inverse analysis: determination of the** G*<sup>c</sup>*

*Modeling and Simulation in Engineering - Selected Problems*

between experiment and simulation is small.

difference of kinetic energy:

*m*spec *v*spec *t*<sup>1</sup>

<sup>2</sup> � <sup>1</sup>

and so we replaced *Ac* by the measured surface size.

2

quently, we state ΔK ¼ *Wc* to be the fracture energy of the specimen.

<sup>Δ</sup>K ¼ <sup>1</sup> 2

*Wc=Ac* with *Ac* <sup>¼</sup> <sup>2</sup>*πr*<sup>2</sup>

**20**

**6.1 Derivation of the specific fracture energy from spallation**

ranges 12<*R<sup>t</sup>*

*(left), Gc* ¼ 90N*=*m *(right).*

**Figure 14.**

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

For the inverse analysis we proceed as follows: we use the data obtained in our previous simulations to define a range of input data Ginp *<sup>c</sup>* , simulate the spallation experiment, and from the computed velocities and masses we deduce via Eq. (33) the "measured" value Gsim *<sup>c</sup>* . To this end we determine the velocities of the fragments by

$$v\_{\rm spec} = \left(2\sum\_{\epsilon} \mathcal{K}\_{\epsilon} / m\_{\rm spec}\right)^{1/2}, v\_{\rm frz,1} = \left(2\sum\_{\epsilon(\rm frz,1)} \mathcal{K}\epsilon / m\_{\rm frz,1}\right)^{1/2}, \text{and}, v\_{\rm frz,2} = \left(2\sum\_{\epsilon(\rm frz,2)} \mathcal{K}\_{\epsilon} / m\_{\rm frz,2}\right)^{1/2}, \text{with the kinetic energy of the finite elements, } \mathcal{K}\_{\epsilon} = \frac{1}{2}m\_{\epsilon}|\mathbf{v}\_{\epsilon}|^{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 *β* ¼ 1. The cohesive stress is *σ<sup>c</sup>* ¼ 15MPa to ensure crack formation at the expected position and we vary the critical opening displacement *δ<sup>c</sup>* to obtain values of G*<sup>c</sup>* between 20 N/m and 150 N/m. All computations result in similar fragments with masses *m*fra,1 ∈ ½ � 895, 905 g and *m*fra,2 ∈½ � 572, 582 g. Also the fractured surfaces show little variations, *Ac* <sup>∈</sup> ½ � <sup>9</sup>*:*5, 11*:*<sup>5</sup> cm2 .

In **Figure 16** the computed specific fracture energy Gsim *<sup>c</sup>* is shown. In the optimal case is Ginp *<sup>c</sup>* <sup>¼</sup> <sup>G</sup>sim *<sup>c</sup>* , this is marked with a dotted line. The computed values show a clear proportionality between Ginp *<sup>c</sup>* and its simulated counterpart <sup>G</sup>sim *<sup>c</sup>* . This generally validates the chosen method of energy balance. However, the value of G*<sup>c</sup>* 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 Ginp *<sup>c</sup>* <sup>&</sup>lt; <sup>G</sup>sim *<sup>c</sup>* with discrepancies between 30 and 60%, see **Figure 16**.

#### **6.3 Phase-field fracture**

For phase-field simulations we use a finer mesh of 5760 elements, *h* ¼ 5*=*6 mm and set *<sup>ε</sup>* <sup>¼</sup> <sup>2</sup>*:*<sup>5</sup> mm and *<sup>M</sup>* <sup>¼</sup> <sup>0</sup>*:*<sup>6</sup> <sup>μ</sup>s�1. At first we accelerate the specimen to *v*spec ¼ 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 Ginp *<sup>c</sup>* , the lower it is the easier both fragments split, cf. **Figure 14**.

inaccurate fragment velocities. A stable state we found for *s t*ð Þ<sup>1</sup> ∈½ � 0*:*75, 0*:*9 and so we

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

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

*<sup>c</sup>* <sup>≈</sup>0*:*9Ginp

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

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 *δ<sup>c</sup>* is implicitly given via the specific crack energy 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 *δ*<sup>0</sup> needs to be set. All these parameters have a

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

The phase-field approach to fracture is based on an evolution equation that

remaining relevant parameters are the mobility *M*, the specific crack energy G*c*, and the length-scale parameter *ε*. The mobility has only numerical character and

essentially refers to the elastic strain energy density Ψ*<sup>e</sup>*

to employ a phase-field fracture simulation for a quantitative analysis of the

*<sup>c</sup>* <sup>¼</sup> <sup>G</sup>sim

*<sup>c</sup>* <20N/m the specimen cracks basically without dissipation

*<sup>c</sup>* prevent any crack. In the range of interest, <sup>G</sup>inp

*<sup>c</sup>* over the input value

*<sup>c</sup>* . The value

*c* ¼

*<sup>c</sup>* is marked with a dotted line and

*<sup>c</sup>* with relative difference less than

*<sup>c</sup>* correlate very well with an R-coefficient

*<sup>c</sup>* . Therefore, it is possible and useful

ð Þ *u*, *s* of the material. The

*<sup>c</sup>* and <sup>G</sup>sim

In **Figure 17** the computed specific fracture energy Gsim

the computed values show a clear proportionality between Ginp

*<sup>c</sup>* and <sup>G</sup>inp

set *s t*ð Þ¼ <sup>1</sup> 0*:*8 and *s t*ð Þ¼ <sup>2</sup> 0*:*2.

very small values of Ginp

whereas too high values of Ginp

30 … 100N/m both values Gsim

of 0.9969. Specifically we have Gsim

methods are summarized in the following.

*<sup>c</sup>* is shown. Again, the optimal case <sup>G</sup>inp

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

*<sup>c</sup>* is now lower than the input value <sup>G</sup>inp

Ginp

of Gsim

experiments.

**7. Conclusions**

**7.1 Model parameters**

clear physical meaning.

simulation results.

**23**

#### **Figure 16.**

*Determination of the specific fracture energy with the cohesive element technique: displayed are the values* Gsim *c derived from the velocity and mass data obtained in the simulation versus the input parameter* Ginp *<sup>c</sup> of the experiment. The dotted line marks the identity* Gsim *<sup>c</sup>* <sup>¼</sup> <sup>G</sup>inp *<sup>c</sup> ; the dashed line is approximated with R*<sup>2</sup> ≈0*:*9928 *and corresponds roughly to* Gsim *<sup>c</sup>* <sup>¼</sup> <sup>1</sup>*:*5Ginp *<sup>c</sup> .*

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 t*ð Þ<sup>1</sup> and *s t*ð Þ<sup>2</sup> . Here we consider *s t*ð Þ<sup>1</sup> and set for symmetry *s t*ð Þ¼ <sup>2</sup> 1 � *s t*ð Þ<sup>1</sup> . A lower "crack initiation" point *s t*ð Þ<sup>1</sup> clearly needs less time and energy to open the crack. On the one hand, if crack initiation equals crack opening, *s t*ð Þ¼ <sup>1</sup> *s t*ð Þ¼ <sup>2</sup> 0*:*5, that is the crack appears immediately after initiation and the dissipated energy Gsim *c* will be close to zero. The input value Ginp *<sup>c</sup>* 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 t*ð Þ<sup>1</sup> > 0*:*95, *s t*ð Þ<sup>2</sup> > 0*:*05 the crack may not fully open, which gives

#### **Figure 17.**

*Determination of the specific fracture energy with the phase-field fracture approach: displayed are the values* Gsim *<sup>c</sup> derived from the velocity and mass data obtained simulations with ε* ¼ 5*=*2 mm*,s t*ð Þ¼ <sup>1</sup> 0*:*8*,s t*ð Þ¼ <sup>2</sup> 0*:*2 *versus the input parameter* Ginp *<sup>c</sup> of the experiment. The dotted line marks the identity* <sup>G</sup>sim *<sup>c</sup>* <sup>¼</sup> <sup>G</sup>inp *<sup>c</sup> ; the dashed line is approximated with R*<sup>2</sup> ≈0*:*9969 *and corresponds well to* Gsim *<sup>c</sup>* <sup>¼</sup> <sup>0</sup>*:*9Ginp *<sup>c</sup> .*

inaccurate fragment velocities. A stable state we found for *s t*ð Þ<sup>1</sup> ∈½ � 0*:*75, 0*:*9 and so we set *s t*ð Þ¼ <sup>1</sup> 0*:*8 and *s t*ð Þ¼ <sup>2</sup> 0*:*2.

In **Figure 17** the computed specific fracture energy Gsim *<sup>c</sup>* over the input value Ginp *<sup>c</sup>* is shown. Again, the optimal case <sup>G</sup>inp *<sup>c</sup>* <sup>¼</sup> <sup>G</sup>sim *<sup>c</sup>* is marked with a dotted line and the computed values show a clear proportionality between Ginp *<sup>c</sup>* and <sup>G</sup>sim *<sup>c</sup>* . The value of Gsim *<sup>c</sup>* is now lower than the input value <sup>G</sup>inp *<sup>c</sup>* 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 Ginp *<sup>c</sup>* <20N/m the specimen cracks basically without dissipation whereas too high values of Ginp *<sup>c</sup>* prevent any crack. In the range of interest, <sup>G</sup>inp *c* ¼ 30 … 100N/m both values Gsim *<sup>c</sup>* and <sup>G</sup>inp *<sup>c</sup>* correlate very well with an R-coefficient of 0.9969. Specifically we have Gsim *<sup>c</sup>* <sup>≈</sup>0*:*9Ginp *<sup>c</sup>* . Therefore, it is possible and useful to employ a phase-field fracture simulation for a quantitative analysis of the experiments.

### **7. Conclusions**

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 t*ð Þ<sup>1</sup> and *s t*ð Þ<sup>2</sup> . Here we consider *s t*ð Þ<sup>1</sup> and set for symmetry *s t*ð Þ¼ <sup>2</sup> 1 � *s t*ð Þ<sup>1</sup> . A lower "crack initiation" point *s t*ð Þ<sup>1</sup> clearly needs less time and energy to open the crack. On the one hand, if crack initiation equals crack opening, *s t*ð Þ¼ <sup>1</sup> *s t*ð Þ¼ <sup>2</sup> 0*:*5, that is the crack appears immediately after initiation and the dissipated energy Gsim

*Determination of the specific fracture energy with the cohesive element technique: displayed are the values* Gsim

*<sup>c</sup>* <sup>¼</sup> <sup>G</sup>inp

*derived from the velocity and mass data obtained in the simulation versus the input parameter* Ginp

fragments to loose contact and (almost) no dissipation will take place. On the other hand, if we raise *s t*ð Þ<sup>1</sup> > 0*:*95, *s t*ð Þ<sup>2</sup> > 0*:*05 the crack may not fully open, which gives

*Determination of the specific fracture energy with the phase-field fracture approach: displayed are the values*

*<sup>c</sup> derived from the velocity and mass data obtained simulations with ε* ¼ 5*=*2 mm*,s t*ð Þ¼ <sup>1</sup> 0*:*8*,s t*ð Þ¼ <sup>2</sup> 0*:*2

*<sup>c</sup> of the experiment. The dotted line marks the identity* <sup>G</sup>sim

*<sup>c</sup>* <sup>¼</sup> <sup>0</sup>*:*9Ginp *<sup>c</sup> .*

*<sup>c</sup>* would work only as a threshold for the

*<sup>c</sup> ; the dashed line is approximated with R*<sup>2</sup> ≈0*:*9928

will be close to zero. The input value Ginp

*experiment. The dotted line marks the identity* Gsim

*<sup>c</sup>* <sup>¼</sup> <sup>1</sup>*:*5Ginp *<sup>c</sup> .*

*Modeling and Simulation in Engineering - Selected Problems*

*and corresponds roughly to* Gsim

**Figure 16.**

**Figure 17.**

*versus the input parameter* Ginp

*line is approximated with R*<sup>2</sup> ≈0*:*9969 *and corresponds well to* Gsim

Gsim

**22**

*c*

*<sup>c</sup>* <sup>¼</sup> <sup>G</sup>inp

*<sup>c</sup> ; the dashed*

*c*

*<sup>c</sup> of the*

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 *δ<sup>c</sup>* is implicitly given via the specific crack energy 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 *δ*<sup>0</sup> 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 Ψ*<sup>e</sup>* ð Þ *u*, *s* of the material. The remaining relevant parameters are the mobility *M*, the specific crack energy G*c*, and the length-scale parameter *ε*. The mobility has only numerical character and

controls the rate of phase-field decrease during crack formation. The critical lengthscale 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 G*c=ε* 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 G*c*.

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.

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

variational derivation and are still an open problem.

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

Tim Dally†, Carola Bilgen†, Marek Werner\*† and Kerstin Weinberg†

\*Address all correspondence to: marek.werner@uni-siegen.de

† These authors contributed equally.

provided the original work is properly cited.

Lehrstuhl für Festkörpermechanik, Fakultät IV, Department Maschinenbau, Institut für Mechanik und Regelungstechnik - Mechatronik, Siegen, Germany

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

**Author details**

**25**

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

### **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*, *s* there exist two different types of solution. On the one hand it is possible to determine the displacements *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*<*ε* 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*<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 *s*\_ < 0, or by a modified mobility.

### *Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

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.

### **Author details**

controls the rate of phase-field decrease during crack formation. The critical lengthscale 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 G*c=ε* in the crack resistance and therefore *ε* needs to be calibrated carefully. The only material parameter in the phase-field model

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

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

For numerical computation of the coupled fields *u*, *s* there exist two different types of solution. On the one hand it is possible to determine the displacements *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*<*ε* 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,

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

In phase-field fracture by definition a continuous function *s* and a diffuse interface with width *ε* ∝ *δ* 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 *s*\_ < 0, or by a modified mobility.

whole domain. Essential requirement is a very fine mesh with *ε*>*h*.

which, in turn, raises the numerical costs tremendously.

**7.3 Constraints and driving forces**

consideration.

**24**

with a clear physical background is the Griffith energy G*c*.

**7.2 Numerical implementation and computation**

*Modeling and Simulation in Engineering - Selected Problems*

complex.

Tim Dally†, Carola Bilgen†, Marek Werner\*† and Kerstin Weinberg† Lehrstuhl für Festkörpermechanik, Fakultät IV, Department Maschinenbau, Institut für Mechanik und Regelungstechnik - Mechatronik, Siegen, Germany

\*Address all correspondence to: marek.werner@uni-siegen.de

† These authors contributed equally.

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

## **References**

[1] Pandolfi A, Ortiz M. Solid modeling aspects of three-dimensional fragmentation. Engineering with Computers. 1998;**14**:287-308

[2] Scheider I, Brocks W. Cohesive elements for thin-walled structures. Computational Materials Science. 2006; **37**(1–2):101-109

[3] Xu XP, Needleman A. Numerical simulations of dynamic interfacial crack growth allowing for crack growth away from the bond line. International Journal of Fracture. 1995;**74**:253-275

[4] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics. 2002;**69**(7):813-833

[5] Stazi FL, Budyn E, Chessa J, Belytschko T. An extended finite element method with higher-order elements for curved cracks. Computational Mechanics. 2003;**31**: 38-48

[6] Pandolfi A, Ortiz M. An eigenerosion approach to brittle fracture. International Journal for Numerical Methods in Engineering. 2012;**92**: 694-714

[7] Schmidt B, Fraternali F, Ortiz M. Eigenfracture: An eigen deformation approach to variational fracture. SIAM Journal on Multiscale Modeling and Simulation. 2009;**7**:1237-1266

[8] Borden MJ, Verhoosel CV, Scott MA, Hughes TJR, Landis CM. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering. 2012;**217– 220**:77-95

[9] Bourdin B. The variational formulation of brittle fracture: Numerical implementation and extensions. In: Volume 5 of IUTAM Symposium on Discretization Methods for Evolving Discontinuities, IUTAM Bookseries, Chapter 22; Springer, Netherlands. 2007. pp. 381-393

of the tensile strength and fracture energy of concrete at high strain rates. International Journal of Impact Engineering. 2006;**32**(10):1635-1650

*DOI: http://dx.doi.org/10.5772/intechopen.92180*

International Journal for Numerical Methods in Engineering. 1999;**44**:

[27] Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. Journal of the Mechanics

[28] Ferrara A, Pandolfi A. Numerical modelling of fracture of human arteries. Computer Methods in Biomechanics and Biomedical Engineering. 2008;**11**:

and Physics of Solids. 1994;**42**:

[29] Pandolfi A, Weinberg K. A numerical approach to the analysis of failure modes in anisotropic plates. Engineering Fracture Mechanics. 2011;

Rosakis AJ. Three-dimensional modeling of intersonic shear-crack growth in asymmetrically loaded unidirectional composite plates. International Journal of Solids and Structures. 2002;**39**:6135-6157

[31] Bazant ZP. Concrete fracture models: Testing and practice.

Engineering Fracture Mechanics. 2002;

[32] Tvergaard V, Hutchinson JW. The

relation between crack growth resistance and fracture process parameters in elastic-plastic solids. Journal of the Mechanics and Physics of

[33] Rose JH, Ferrante J, Smith JR. Universal binding energy curves for metals and bimetallic interfaces. Physical Review Letters. 1981;**47**:

Bruchmechanik. Technical Report at

Solids. 1992;**40**:1377-1397

[34] Dally T. Vergleich von Kohäsivelement-Methode und Phasenfeld-Methode anhand ausgewählter Probleme der

Universität Siegen; 2017

[30] Yu C, Pandolfi A, Ortiz M, Coker D,

1267-1282

*Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture…*

1397-1434

553-567

**78**:2052-2069

**69**:165-205

675-678

[18] Weerheijma J, Van Doormaal JCAM. Tensile failure of concrete at high loading rates: New test data on strength and fracture energy from instrumented spalling tests. International Journal of Impact Engineering. 2007;**34**:609

[19] Zhang L, Hu S-S, Chen D-X, Yu Z-Q, Liu F. An experimental technique for spalling of concrete. Experimental Mechanics. 2009;**49**(4):

[20] Khosravani MR, Wagner P, Fröhlich D, Weinberg K. Dynamic fracture investigations of ultra-high performance concrete by spalling tests. Engineering Structures. 2019;**201**:

[21] Griffith AA. The phenomena of rupture and flow in solids. Philosophical Transactions of the Royal Society of

[22] Irwin GR. Elasticity and plasticity:

[23] Dugdale DS. Yielding of steel sheets

Mechanics and Physics of Solids. 1960;**8**:

[24] Barenblatt GI. The mathematical theory of equilibrium of cracks in brittle

Computational modelling of impact

[26] Ortiz M, Pandolfi A. A class of cohesive elements for the simulation of three-dimensional crack propagation.

London. 1921;**221**:163-198

Fracture. In: Függe S, editor. Encyclopedia of Physics. 1958

containing clits. Journal of the

fracture. Advances in Applied Mechanics. 1962;**7**:55-129

[25] Camacho GT, Ortiz M.

damage in brittle materials. International Journal of Solids and Structures. 1996;**33**:2899-2938

523-532

109844

100-104

**27**

[10] Henry H, Levine H. Dynamic instabilities of fracture under biaxial strain using a phase field model. Physics Review Letters. 2004;**93**:105505

[11] Karma A, Kessler DA, Levine H. Phase-field model of mode III dynamic fracture. Physical Review Letters. 2001; **81**:045501

[12] Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering. 2010;**199**:2765-2778

[13] Verhoosel CV, de Borst R. A phasefield model for cohesive fracture. International Journal for Numerical Methods in Engineering. 2013. In press

[14] Ambrosio L, Tortorelli VM. Approximation of functional depending on jumps by elliptic functional via γconvergence. Communications on Pure and Applied Mathematics. 1990;**348– 349**:13-16

[15] Diaz-Rubio FG, Perez JR, Galvez VS. The spalling of long bars as a reliable method of measuring the dynamic tensile strength of ceramics. International Journal of Impact Engineering. 2002;**27**:161-177

[16] Klepaczko JR, Brara A. An experimental method for dynamic tensile testing of concrete by spalling. International Journal of Impact Engineering. 2001;**25**(4):387-409

[17] Schuler H, Mayrhofer C, Thoma K. Spall experiments for the measurement *Cohesive Elements or Phase-Field Fracture: Which Method Is Better for Dynamic Fracture… DOI: http://dx.doi.org/10.5772/intechopen.92180*

of the tensile strength and fracture energy of concrete at high strain rates. International Journal of Impact Engineering. 2006;**32**(10):1635-1650

**References**

**37**(1–2):101-109

38-48

694-714

**220**:77-95

**26**

[1] Pandolfi A, Ortiz M. Solid modeling

*Modeling and Simulation in Engineering - Selected Problems*

Symposium on Discretization Methods for Evolving Discontinuities, IUTAM Bookseries, Chapter 22; Springer, Netherlands. 2007. pp. 381-393

[10] Henry H, Levine H. Dynamic instabilities of fracture under biaxial strain using a phase field model. Physics

Review Letters. 2004;**93**:105505

[12] Miehe C, Hofacker M,

**81**:045501

**349**:13-16

[11] Karma A, Kessler DA, Levine H. Phase-field model of mode III dynamic fracture. Physical Review Letters. 2001;

Welschinger F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering. 2010;**199**:2765-2778

[13] Verhoosel CV, de Borst R. A phasefield model for cohesive fracture. International Journal for Numerical Methods in Engineering. 2013. In press

Approximation of functional depending on jumps by elliptic functional via γconvergence. Communications on Pure and Applied Mathematics. 1990;**348–**

[15] Diaz-Rubio FG, Perez JR, Galvez VS. The spalling of long bars as a reliable method of measuring the dynamic tensile strength of ceramics. International Journal of Impact Engineering. 2002;**27**:161-177

[14] Ambrosio L, Tortorelli VM.

[16] Klepaczko JR, Brara A. An experimental method for dynamic tensile testing of concrete by spalling. International Journal of Impact Engineering. 2001;**25**(4):387-409

[17] Schuler H, Mayrhofer C, Thoma K. Spall experiments for the measurement

[2] Scheider I, Brocks W. Cohesive elements for thin-walled structures. Computational Materials Science. 2006;

[3] Xu XP, Needleman A. Numerical simulations of dynamic interfacial crack growth allowing for crack growth away from the bond line. International Journal of Fracture. 1995;**74**:253-275

[4] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics. 2002;**69**(7):813-833

[5] Stazi FL, Budyn E, Chessa J, Belytschko T. An extended finite element method with higher-order elements for curved cracks.

approach to brittle fracture.

Computational Mechanics. 2003;**31**:

International Journal for Numerical Methods in Engineering. 2012;**92**:

[7] Schmidt B, Fraternali F, Ortiz M. Eigenfracture: An eigen deformation approach to variational fracture. SIAM Journal on Multiscale Modeling and Simulation. 2009;**7**:1237-1266

[8] Borden MJ, Verhoosel CV, Scott MA, Hughes TJR, Landis CM. A phase-field description of dynamic brittle fracture.

Mechanics and Engineering. 2012;**217–**

Computer Methods in Applied

[9] Bourdin B. The variational formulation of brittle fracture: Numerical implementation and extensions. In: Volume 5 of IUTAM

[6] Pandolfi A, Ortiz M. An eigenerosion

aspects of three-dimensional fragmentation. Engineering with Computers. 1998;**14**:287-308

[18] Weerheijma J, Van Doormaal JCAM. Tensile failure of concrete at high loading rates: New test data on strength and fracture energy from instrumented spalling tests. International Journal of Impact Engineering. 2007;**34**:609

[19] Zhang L, Hu S-S, Chen D-X, Yu Z-Q, Liu F. An experimental technique for spalling of concrete. Experimental Mechanics. 2009;**49**(4): 523-532

[20] Khosravani MR, Wagner P, Fröhlich D, Weinberg K. Dynamic fracture investigations of ultra-high performance concrete by spalling tests. Engineering Structures. 2019;**201**: 109844

[21] Griffith AA. The phenomena of rupture and flow in solids. Philosophical Transactions of the Royal Society of London. 1921;**221**:163-198

[22] Irwin GR. Elasticity and plasticity: Fracture. In: Függe S, editor. Encyclopedia of Physics. 1958

[23] Dugdale DS. Yielding of steel sheets containing clits. Journal of the Mechanics and Physics of Solids. 1960;**8**: 100-104

[24] Barenblatt GI. The mathematical theory of equilibrium of cracks in brittle fracture. Advances in Applied Mechanics. 1962;**7**:55-129

[25] Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. International Journal of Solids and Structures. 1996;**33**:2899-2938

[26] Ortiz M, Pandolfi A. A class of cohesive elements for the simulation of three-dimensional crack propagation.

International Journal for Numerical Methods in Engineering. 1999;**44**: 1267-1282

[27] Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. Journal of the Mechanics and Physics of Solids. 1994;**42**: 1397-1434

[28] Ferrara A, Pandolfi A. Numerical modelling of fracture of human arteries. Computer Methods in Biomechanics and Biomedical Engineering. 2008;**11**: 553-567

[29] Pandolfi A, Weinberg K. A numerical approach to the analysis of failure modes in anisotropic plates. Engineering Fracture Mechanics. 2011; **78**:2052-2069

[30] Yu C, Pandolfi A, Ortiz M, Coker D, Rosakis AJ. Three-dimensional modeling of intersonic shear-crack growth in asymmetrically loaded unidirectional composite plates. International Journal of Solids and Structures. 2002;**39**:6135-6157

[31] Bazant ZP. Concrete fracture models: Testing and practice. Engineering Fracture Mechanics. 2002; **69**:165-205

[32] Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. Journal of the Mechanics and Physics of Solids. 1992;**40**:1377-1397

[33] Rose JH, Ferrante J, Smith JR. Universal binding energy curves for metals and bimetallic interfaces. Physical Review Letters. 1981;**47**: 675-678

[34] Dally T. Vergleich von Kohäsivelement-Methode und Phasenfeld-Methode anhand ausgewählter Probleme der Bruchmechanik. Technical Report at Universität Siegen; 2017

[35] Pandolfi A, Ortiz M. An efficient adaptive procedure for threedimensional fragmentation simulations. Engineering with Computers. 2002;**18**: 48-159

[36] Bourdin B, Francfort GA, Marigo JJ. The variational approach to fracture. Journal of Elasticity. 2008;**9**:5-148

[37] Francfort GA, Marigo J-J. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids. 1998; **46**:1319-1342

[38] Borden MJ, Hughes TJR, Landis CM, Verhoosel CV. A higher-order phasefield model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Computer Methods in Applied Mechanics and Engineering. 2014; **273**(0):100-118

[39] Hesch C, Stefan S, Dittmann M, Franke M, Weinberg K. Isogeometric analysis and hierarchical refinement for higher-order phase-field models. Computer Methods in Applied Mechanics and Engineering. 2016;**303**: 185-207

[40] Weinberg K, Hesch C. A high-order finite-deformation phase-field approach to fracture. Continuum Mechanics and Thermodynamics. 2017;**29**:935-945

[41] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering. 2010;**83**: 1273-1311

[42] Bilgen C, Kopanicakova A, Krause R, Weinberg K. A phase-field approach to conchoidal fracture. Meccanica. 2018;**53**:1203-1219

[43] Weinberg K, Dally T, Schuß S, Werner M, Bilgen C. Modeling and numerical simulation of crack growth and damage with a phase field approach. GAMM-Mitteilungen. 2016;**39**(1):55-77

**Chapter 2**

**Abstract**

stochastic model

**1. Introduction**

**29**

ally divided into several types:

(segregation, classification).

Methods of Nonequilibrium

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

Mixing Bulk Components

Statistical Mechanics in Models for

When describing the mechanics of the behavior of bulk materials during their mixing, a theoretical basis for the design of the specified equipment is formed. In recent years, the most well-known methods of modeling this process include the stochastic approach, in the framework of which models of the following types are actively developing: cell, managerial, with time series, energy, etc. Moreover, as a rule, predicting the quality of the finished mixture according to the selected

criterion is achieved by using numerical calculation methods based on the generated cyber system. Of particular interest is the use of the energy method from the statistical mechanics of nonequilibrium processes due to the possibility of obtaining analytical simulation results. The paper describes the motion models of bulk components in rarefied flows, which are built on the basis of the energy method and

Mixing of bulk materials belongs to the category of mechanical processes in the field of chemical technology in a variety of manufacturing enterprises, including paint and varnish, glass, construction, food, pharmaceutical, chemical, etc. The result of this technological operation is a loose mixture with a regulated degree of homogeneity, the quality of which is influenced by many factors. Moreover, factors that have a decisive influence on the quality of the final product can be condition-

• Properties of components—physical, mechanical (particle-size distribution, tendency to adhesion, agglomeration, degree of moisture, etc.), and chemical.

working chamber of the apparatus under the established modes of its operation

• Features of the behavior of different-sized particles when moving in the

take into account the main characteristics of the studied mixing process.

**Keywords:** apparatus, process, mixing, bulk material, rarefied flow,

[44] Weinberg K, Khosravani MR. On the tensile resistance of UHPC at impact. The European Physical Journal Special Topics. 2018;**227**:167-177

[45] Park KC, Lim SJ, Huh H. A method for computation of discontinuous wave propagation in heterogeneous solids: Basic algorithm description and application to one-dimensional problems. International Journal for Numerical Methods in Engineering. 2012;**91**:622-643

[46] Erzar B, Forquin P. An experimental method to determine the tensile strength of concrete at high rates of strain. Experimental Mechanics. 2010;**50**(7):941-955

[47] Bilgen C, Weinberg K. On the crack-driving force of phase-field models in linearized and finite elasticity. Computer Methods in Applied Mechanics and Engineering. 2019; **353**(15):348-372

### **Chapter 2**

[35] Pandolfi A, Ortiz M. An efficient adaptive procedure for three-

48-159

**46**:1319-1342

**273**(0):100-118

185-207

1273-1311

**28**

dimensional fragmentation simulations. Engineering with Computers. 2002;**18**:

*Modeling and Simulation in Engineering - Selected Problems*

numerical simulation of crack growth and damage with a phase field approach. GAMM-Mitteilungen. 2016;**39**(1):55-77

[44] Weinberg K, Khosravani MR. On the tensile resistance of UHPC at impact. The European Physical Journal Special Topics. 2018;**227**:167-177

[45] Park KC, Lim SJ, Huh H. A method for computation of discontinuous wave propagation in heterogeneous solids: Basic algorithm description and application to one-dimensional problems. International Journal for Numerical Methods in Engineering.

experimental method to determine the tensile strength of concrete at high rates of strain. Experimental Mechanics.

[47] Bilgen C, Weinberg K. On the crack-driving force of phase-field models in linearized and finite elasticity.

Computer Methods in Applied Mechanics and Engineering. 2019;

2012;**91**:622-643

2010;**50**(7):941-955

**353**(15):348-372

[46] Erzar B, Forquin P. An

[36] Bourdin B, Francfort GA, Marigo JJ. The variational approach to fracture. Journal of Elasticity. 2008;**9**:5-148

[37] Francfort GA, Marigo J-J. Revisiting

[38] Borden MJ, Hughes TJR, Landis CM, Verhoosel CV. A higher-order phasefield model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Computer Methods in Applied Mechanics and Engineering. 2014;

[39] Hesch C, Stefan S, Dittmann M, Franke M, Weinberg K. Isogeometric analysis and hierarchical refinement for higher-order phase-field models. Computer Methods in Applied

Mechanics and Engineering. 2016;**303**:

[40] Weinberg K, Hesch C. A high-order finite-deformation phase-field approach to fracture. Continuum Mechanics and Thermodynamics. 2017;**29**:935-945

[41] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering. 2010;**83**:

[42] Bilgen C, Kopanicakova A, Krause R, Weinberg K. A phase-field approach to conchoidal fracture. Meccanica. 2018;**53**:1203-1219

[43] Weinberg K, Dally T, Schuß S, Werner M, Bilgen C. Modeling and

minimization problem. Journal of the Mechanics and Physics of Solids. 1998;

brittle fracture as an energy
