Open access peer-reviewed chapter

# Mode-I and Mode-II Crack Tip Fields in Implicit Gradient Elasticity Based on Laplacians of Stress and Strain. Part III: Numerical Simulations

Written By

Carsten Broese, Jan Frischmann and Charalampos Tsakmakis

Submitted: 10 July 2020 Reviewed: 17 August 2020 Published: 21 September 2020

DOI: 10.5772/intechopen.93619

From the Edited Volume

## Nanomechanics - Theory and Application

Edited by Alexander V. Vakhrushev

Chapter metrics overview

View Full Metrics

## Abstract

A two-dimensional formulation of the 3-PG Model of implicit gradient elasticity has been developed in Part I. The predicted near-tip fields for Mode-I and Mode-II crack problems have been derived in Part II. It has been found that both the classical Cauchy stress and the nonclassical double stress are singular with the order of singularity r−12. In the present chapter, the two-dimensional model formulation is implemented in a finite element code. For verification of the resulting finite element model, a square section with a circular hole subjected to displacement-controlled tension loading is considered and discussed. The main concerns of the chapter are, on the one hand, to validate the analytical solutions of Part II. On the other hand, the chapter aims to investigate the effect of nonclassical material parameters on the stress intensity factors.

### Keywords

• finite elements
• square section with a hole
• mode-I and mode-II crack problems
• stress intensity factors
• angular functions

## 1. Introduction

Mode-I and Mode-II crack problems have been discussed analytically in Part II for the 3-PG Model of implicit gradient elasticity. Solutions for the near-tip fields have been obtained by employing the method of asymptotic expansion of Williams’ type (see Williams [1]). It has been proved that both the classical Cauchy stress and the nonclassical double stress are singular with the order of singularity r12. Even more, the first two terms in the asymptotic expansion of the Cauchy stress are identical to those in the context of classical elasticity. The leading terms of the asymptotic expansions of the classical and the nonclassical stresses are represented by the so-called stress intensity factors.

The present chapter deals with numerical simulations that employ the 3-PG Model to solve crack problems of Mode-I and Mode-II types. A finite element model for plane strain is developed in the framework of a weak formulation based on the principle of virtual work. To verify the finite element formulation and implementation, a representative example is considered: a square section with a circular hole subjected to tension loading. The predicted stress concentration factors are compared with the corresponding stress concentration factors predicted by classical elasticity. The main objectives are to confirm the assumptions and the analytical results of Part II as well as to assess the effect of the nonclassical material parameters of the 3-PG Model on the stress intensity factors for Mode-I and Mode-II crack problems. It is perhaps of interest to remark that similar investigations for the case of micropolar elastic continua are provided in Diegele et al. [2].

The scope of the chapter is organized as follows: Section 2 gives some details about the implementation of the 3-PG Model in a finite element code. A square section with a circular hole subjected to tension loading is discussed in Section 3. The numerical simulations verify the finite element model and its ability to predict length scale effects. Further, they provide a first comparison to classical elasticity by calculating the corresponding stress intensity factors. Section 4 is devoted to an analysis of edge-cracked specimens. The analysis comprises, among others, the effect of material parameters on the stress intensity factors. Moreover, it indicates a very good agreement between the numerical and the analytical solutions of the angular functions. This confirms, a posteriori, the assumed symmetry conditions of the micro-deformation. The chapter closes with some concluding remarks in Section 5.

Throughout the chapter, the same notation as introduced in Part I applies.

## 2. Finite element formulation

The following formulations refer to three dimensions and apply especially to two-dimensional cases when conditions for plane strain are imposed. Let us consider once more the equilibrium equations (see Section 3.1 “The 3-PG Model as particular case of micro-strain elasticity” in Part I) to be solved,

iΣij=0,E1
iμijk+σjk=0,E2

and specify the corresponding boundary conditions (see Section 4.7 “Boundary conditions” of Part I) as follows:

ui=ui0onVui,njΣji=Pi0onVPi,E3
Ψij=Ψij0onVΨij,nkμkij=Tij0onVTij,E4

with

VuiVPi=V,VuiVPi=Ø,E5
VΨijVTij=V,VΨijVTij=Ø.E6

These equations reflect Dirichlet boundary conditions for the macro-displacement and the micro-deformation and Neumann boundary conditions for the Cauchy stress and the couple stress.

The first step toward a finite element formulation is to elaborate a weak form of the above boundary value problem.

### 2.1 Weak form of the boundary value problem

Let the fields uix and Ψijx as well as the so-called virtual fields δujx and δΨjkx belong to the following function spaces:

uiSuiuiH1V¯ui=ui0onVui,E7
ΨijTΨijΨijH1V¯Ψij=Ψij0onVΨij,E8
δujVδujδujH1V¯δuj=0onVui,E9
δΨjkWδΨjkδΨjkH1V¯δΨjk=0onVΨij.E10

Now, multiply Eq. (1) by δuj and Eq. (2) by δΨjk and take the integrals over V,

VδujiΣijdv=0,E11
VδΨjkiμijk+σjkdv=0.E12

Next, add up these two equations, use partial integration, the divergence theorem, and the boundary conditions (3) and (4) to receive

ViδujΣijdvVδΨjkσjkdv+ViδΨjkμijkdvVPδujPj0daVTδΨjkTjk0da=0.E13

As usually, in favor of a short notation, we use the integral over VP to indicate the summation of single integrals over VPi, which generally do not coincide. The meaning of the integration over VT is analogous. Eq. (13) is the weak form of the boundary value problem and is the starting point of the finite element formulation.

### 2.2 Discretization

According to the finite element method (see, e.g., Hughes [3]), the domain V is approximated by nel finite elements Ve, so that

VVhe=1nelVe.E14

The elements are connected to each other at selected points called nodal points, or simply nodes, and the following notation holds.

Ku: Set of global node numbers with macro-displacement degree of freedom.

KΨ: Set of global node numbers with micro-deformation degree of freedom.

The exact solutions ui and Ψij are approximated by

uixuihxAKuNAuxuiA,E15
ΨijxΨijhxAKΨNAYxΨijA,E16

where uiA and ΨijA are the unknown values of ui and Ψij at node A. The so-called shape functions NAu and NAΨ belong to finite-dimensional function spaces Sh and Th, which approximate the function spaces S and T, respectively. Similarly, δuj and δΨjk are approximated by

δujxδujhxBKuNBuxδujB,E17
δΨjkxδΨjkhxBKΨNBuxδΨjkB,E18

with δujB and δΨjkB being constants. The functions δujh and δΨjkh are elements of finite-dimensional function spaces Vh and Wh, which approximate the function spaces V and W, respectively.

We may use the above approximations (17) and (18) to rewrite Eq. (13) in the form

e=1nelVeiNBuδujBΣijdvVeNBΨδΨjkBσjkdv+VeiNBΨδΨjkBμijkdvVePNBuδujBPj0daVeTNBΨδΨjkBTjk0da}=0,E19

with the meaning of VeP and VeT being obvious.

We next employ the elasticity laws (see Section 3.1 “The 3-PG Model as particular case of micro-strain elasticity” of Part I)

Σij=c2c1Cijmnεmnc2c1c1CijmnΨmn,E20
σjk=c2c1c1CjkmnεmnCjkmnΨmn,E21
μijk=c2c1iΨmnCmnjk,E22

in order to replace the stresses Σij and σjk as well as the double stresses μijk in Eq. (19):

e=1nelVeiNBuδujBc2c1Cijmnmunc2c1c1CijmnΨmndvc2c1c1VeNBΨδΨjkBCjkmnmunCjkmnΨmndv+c2c1VeiNBΨδΨjkBiΨmnCmnjkdvVePNBuδujBPj0daVeTNBΨδΨjkBTjk0da}=0.E23

Finally, we use the approximations (15) and (16), to find

or equivalently,

e=1nelδujBc2c1VeiNBuCijmnmNAudvunAδujBc2c1c1VeiNBuCijmnNAΨdvΨmnAδΨjkBc2c1c1VeNBΨCjkmnmNAudvunA+δΨjkBc2c1c1VeNBΨCjkmnNAΨdvΨmnA+δΨjkBc2c1VeiNBΨiNAΨCjkmndvΨmnAδujBVePNBuPj0daδΨjkBVeTNBΨTjk0da=0.E25

The form of this equation suggests to define the element stiffness matrices:

KjBnAuu,ec2c1VeiNBuCijmnmNAudv,E26
KjBmnAuΨ,ec2c1c1VeiNBuCijmnNAΨdv,E27
KjkBnAΨu,ec2c1c1VeNBΨCjkmnmNAudv,E28
KjkBmnAΨΨ,ec2c1c1VeNBΨCjkmnNAΨ+c1iNBΨiNAΨCjkmndv,E29

and the element force vectors:

FjBu,eVePNBuPj0da,E30
FjkBΨ,eVeTNBΨTjk0da,E31

and to recast Eq. (25) as

e=1nelδujBKjBnAuu,eunA+KjBmnAuΨ,eΨmnAFjBu,e+δΨjkBKjkBnAΨu,eunA+KjkBmnAΨΨ,eΨmnAFjkBΨ,e=0.E32

As δujB and δΨjkB may be chosen arbitrary, we obtain the system of equations:

e=1nelKjBnAuu,eunA+KjBmnAuΨ,eΨmnAFjbu,e=0,E33
e=1nelKjkBnAΨu,eunA+KjkBmnAΨΨ,eΨmnAFjkbΨ,e=0.E34

This means, that the solution of the boundary value problems (1)(6) is reduced to the solution of the systems (33) and (34).

The proceeding approximation of the weak form, constrained to plane strain states, has been implemented in the finite element code FEAP. Isoparametric elements are employed, that is, the space coordinates are represented by using the shape functions, NAu,

xixihAKuNAuxiA,E35

where xiA are constants.

## 3. Square section with a circular hole

The main objective of this section is to validate the implemented finite element code. To this end, we consider the plane strain problem shown in Figure 1a, where the square section with a circular hole, located at the center of the section, is stretched in the y-direction. The length of the section is b=5mm, while the radius of the hole is r=0.25mm. With respect to the Cartesian coordinate system xy, the boundaries x=±b2 are assumed to be free of classical and nonclassical tractions. At the boundary y=b2, the macro-displacement component uy, the component Px of the classical traction, and the nonclassical traction components Tij are assumed to vanish. At the boundary y=b2, the macro-displacement in the y-direction is given by uy=0.1mm, while Px=0 and Tij=0 are imposed. The whole circular hole is assumed to be free of classical and nonclassical tractions. For a small circular hole, a nearly uniform stress component

Σ0Σyyy=b2E36

will be required to realize the given boundary conditions.

The most simple case in classical elasticity, analogous to the boundary value problem above, is to consider the square section in the context of a plane stress problem subjected to the traction boundary condition Σyyy=b2=Σ0. Attention is focused on the distribution of Σyy along the section A-A (see Figure 1a), as a function of the local coordinate

axr,E37

with xr. A so-called stress concentration factor k is defined by

kΣyyΣ0,ΣyyΣyya=0y=0,E38

and turns out to be k=3 (see, e.g., Gould [4], p. 124).

Now consider the square section in Figure 1a in the context of a plane strain problem with the boundary conditions stated in the first paragraph within classical elasticity. The stress distribution along the section A-A has been determined by employing the standard elastic element of FEAP using the classical material parameters

E=100GPa,ν=0,3.E39

As depicted schematically in Figure 1b, a radial mesh of 16 stripes of quadratic elements with reduced integration, that is, elements with 8 nodes, is used. Every stripe consists of 180 elements and the elements are chosen to decrease in size the closer to the circular hole the elements are placed. In the context of classical elasticity, a mesh consisting of 720 such stripes is used, whereas the mesh for the simulations of the 3-PG Model only needs 72 such stripes of elements.

Figure 2 illustrates the distribution of the dimensionless stress ΣyyΣ0 along the section A-A. The value of ΣyyΣ0 at a=0 represents the stress concentration factor k and turns out to be k=3,14 for classical elasticity.

We consider next corresponding distributions predicted by the 3-PG Model. To this end, we employ the finite element model developed in the last section. The classical material parameters are given by Eq. (39). For the purposes of the present chapter, we find it convenient to use the nonclassical parameters c1 and

c3c2c1c1.E40

Note that since c2>c1>0 (see Section 3.1 “The 3-PG Model as particular case of micro-strain elasticity” of Part I), the constrain c3>0 applies. It becomes apparent, from the elasticity laws (20)(22), that for c30, the 3-PG Model will approach to classical elasticity. Hence, we expect the distributions of the dimensionless stress ΣyyΣ0 along the section A-A to be close to the classical one whenever c3 is sufficiently small. Indeed, Figure 2 confirms this expectation: It can be seen, that for sufficiently large values of a, the two graphs almost coincide. Another way to illustrate this issue is to consider the effect of the nonclassical material parameter c3 on the values of the stress concentration factor k. We expect that for c1=const. and c30, the values of k will approach to the classical value k=3,14. This is exactly what we can observe in Figure 3. Alternatively, we can consider the effect of the material parameter c1 on the values of the stress concentration factor k. For the assumed geometry of the specimen, the imposed boundary conditions and for a fixed value c3=1, the effect of c1 on k is illustrated in Figure 4 by the graph referred to as specimen 1. A convenient way to illustrate the effect of c1 is to consider further boundary value problems. Thus, we consider in addition three further problems, denoted as specimen 4, 20, and 200, which arise by multiplying the given geometry of the specimen and the imposed boundary conditions with the factors n=4,20 and 200, respectively. The corresponding graphs of k as a function of c1 are displayed in Figure 4 and are referred to as specimens 4, 20, and 200, respectively. Keeping in mind, that c1 is an internal material length, we rescale the abscissa c1 by considering the graphs of k as a function of c1n2. We expect, that all distributions should coincide, and in fact this is shown in Figure 5.

Altogether, the calculated responses reflect the expected results, which in turn provides a validation of the developed and implemented finite element approximation of the 3-PG Model.

## 4. Finite element analysis of crack problems

### 4.1 Edge-cracked specimen

The remaining analysis is referred to the edge-cracked specimen shown in Figure 6a. The assumed width and length of the specimen are, respectively, b=11mm and 2h=20mm, while the crack length is chosen to be a=1mm. The origin of the Cartesian and the cylindrical coordinate systems, to which we refer, is located at the crack tip.

The specimen is discretized by two meshes (see Figure 6b), a rectangular mesh for the main part of the specimen, consisting of 320 quadratic eight-node elements (with reduced integration) and a radial mesh around the crack tip. The radial mesh, schematically shown in Figure 6c, consists of 16 stripes of quadratic elements with reduced integration, that is, elements with eight nodes, are used. Every strip consists of 45 elements, which decrease in size the closer to the crack tip they are placed. The elements, which contain the crack tip itself, are singular and are the so-called quarter point elements. That means, that one edge of each quadratic element degenerates into a point and two intermediate nodes of two adjacent edges are repositioned. To be more specific, they are moved to a position only a quarter of the edge length away from the crack tip, as indicated in Figure 6d. Due to this repositioning, the corresponding shape functions differ from those of regular 8–node elements. Specifically, the shape functions of quarter point elements consist of linear shape functions extended by a term proportional to the square root of the corresponding component of the position vector (see Henshell and Shaw [5] or Barsoum [6]).

For all calculations in the remainder of the chapter, the values of the classical material parameters E and ν are given by Eq. (39) and the following Dirichlet boundary conditions are imposed:

urr=0=uφr=0=0,uφr=10mm,φ=0=0,E41

which excludes the possibility of rigid body motions.

The crack faces are subjected to the following boundary conditions for the classical stress: Mode-I near-tip fields are produced by imposing an internal pressure p=100MPa to act, while Mode-II near-tip fields are enforced by subjecting the crack faces to a shear stress loading of 100MPa. All nonclassical traction components are supposed to vanish for both Mode-I and Mode-II crack types.

### 4.2 Stress intensity factors

Finite element simulations of Mode-I and Mode-II crack problems allow to verify, numerically, the order of singularity and the angular functions as well as to determinate the stress intensity factors. First, we verify the order of singularity and show how to determine stress intensity factors numerically. Motivated by the analytical solutions in Section 5 “Discussion of the asymptotic solutions” of Part II, assume that the stress Σ and the double stress μ in the vicinity of the crack tip are given by

Σαβ=2πrp1K˜IfαβIφ+K˜IIfαβIIφ,E42
μαβγ=2πrp1L˜IgαβγIφ+L˜II,1gαβγII,1φ+L˜II,2gαβγII,2φ.E43

Then,

Σφφφ=0=2πrp1K˜I,E44
Σφ=0=2πrp1K˜II,E45
μφrφφ=0=2πrp1L˜I,E46
12μφrr+μφφφφ=0=2πrp1L˜II,1,E47
12μφrrμφφφφ=0=2πrp1L˜II,2,E48

and hence

logΣφφφ=0=p1log2πr+logK˜I,E49
logΣφ=0=p1log2πr+logK˜II,E50
logμφrφφ=0=p1log2πr+logL˜I,E51
log12μφrr+μφφφφ=0=p1log2πr+logL˜II,1,E52
log12μφrrμφφφφ=0=p1log2πr+logL˜II,2.E53

These equations indicate, respectively, a linear respone with slope p1. This means that, as in the case of classical elastic fracture mechanics, one can fit the exponent p on the basis of values of stress components calculated by the finite element method at the nodes ahead of the crack tip (see, e.g., Figure 7). Computed responses with the finite element model have confirmed the value p=12 with great accuracy. Therefore, the value p=12 will be fixed, to avoid errors owing to inaccurate numerical determination of p when discussing predicted responses. All stress intensity factors are determined from Eqs. (49)(53) for the fixed value p=12 by applying the least square method. It can be recognized from Figure 7, that the linear response of the stress applies for a radius r107101. The numerical determination of the angular functions below is referred to a fixed radius r=5106mm. The verification of the angular functions and the effect of material parameters on the stress intensity factor is discussed seperately for Mode-I and Mode-II crack problems.

### 4.3 Results for mode-I

Figures 811 display plots of the angular functions fαβI, hαβI and gαβγI for both the analytical and the finite element solutions. The plots of the latter are constructed by dividing the values of Σαβ, ΨαβΨ¯αβ, and μαβγ at a radius r=5106mm by K˜I2πr, r2πL˜Ic2c12μr2πL˜Ic1c32μ and L˜I2πr, respectively. The stress intensity factors K˜I and L˜I have been determined by least square fitting as described in the last section. In the finite element computations, the values c1=104mm2 and c3=1 for the nonclassical material parameters are chosen. The values of Ψ¯ij are determined to be Ψ¯11=0,001374, Ψ¯22=0,00359, and Ψ¯12=0. The general observation is that there is good agreement between the analytical and the finite element predictions. The nonvanishing values of Ψ¯ijΨ¯ijI verify the existence of this constant terms. The fact that the analytical and the numerical results of the angular functions hαβI agree very well verifies the assumed symmetry conditions for Ψαβ (see Section 4.8 “Symmetry conditions” in Part I).

It is worth remarking that even though the asymptotic solutions of Σαβ have the same form as the ones of classical elasticity, the corresponding values the stress intensity factors K˜I and K˜II of the 3-PG Model will be, in general, different from the stress intensity factors KI and KII of classical elasticity. This difference results from the fact that the 3-PG Model and the classical elasticity imply, in general, different distributions of the components Σαβ for identical classical boundary conditions.

Recall also, that the 3-PG Model includes two material parameters more than the classical elasticity, namely, c1 and c3. The effect of these nonclassical material parameters is illustrated in Figures 12 and 13. These figures reveal that for c3=const. and very large values of c1 or for c1=const. and very small values of c3, the values of K˜I converge to the values of KI. In other words, the responses of the 3-PG Model approach the ones according to classical elasticity.

The effect of c1 and c3 on the nonclassical stress intensity factor L˜I is illustrated in Figures 14 and 15. The principal observation is that for smaller values of c3, the value of L˜I gets smaller as well. On the other hand, if c3=const., then the values of L˜I increase with increasing values of c1.

### 4.4 Results for mode-II

The graphs of the angular functions fαβII for both the analytical and the finite element solutions are shown in Figure 16. It can be seen that the graphs fit very well.

Since there are two stress intensity factors of Mode-II, L˜II,1, and L˜II,2, we find it convenient, in the vicinity of the crack tip, to consider the angular distributions of ΨαβΨ¯αβ and μαβγ themselves instead of the corresponding distributions of the angular functions hαβII,1, hαβII,2, gαβγII,1, and gαβγII,2. For r=const., Figures 1719 illustrate angular distributions of Ψαβ and μαβγ for both the analytical and the numerical solutions. Once more, we can recognize, that the analytical and the numerical results agree with great accuracy, which also verifies the assumed symmetry conditions of Mode-II for Ψαβ (see Section 4.8 “Symmetry conditions” in Part I). Further, the values of Ψ¯ij have been determined to be Ψ¯11=0, Ψ¯22=0, and Ψ¯12=0,00391 and this, in turn, verifies the existence of the constant terms Ψ¯ijΨ¯ijII.

The effect of c1 and c3 on K˜II is illustrated in Figures 20 and 21 and is similar to the effect of c1 and c3 on K˜I (cf. Figures 12 and 13).

The effect of c1 and c3 on the nonclassical stress intensity factors L˜II,1 and L˜II,2 is illustrated in Figures 2225. Again, this effect is quite similar to the effect of c1 and c3 on L˜I, but the stress intensity factors L˜II,1 and L˜II,2 both are negative in contrast to the positive stress intensity factor L˜I.

With regard to Mode-II crack problems, it is of interest to recall that the general solutions of Ψαβ0 depend on three constants of integration, that is, B0, E0, and F0. By virtue of the boundary conditions on the crack faces, we obtained the relation F0=E0, which is independent of the crack geometry, that is, the crack length. We do not know further conditions to relate B0 and E0 with each other. However, keeping in mind that, in general, stress intensity factors depend on the crack geometry and the applied loading conditions, one may ask if there exists a relation between B0 and E0, or equivalently between L˜II,1 and L˜II,2, which is independent of the crack length. To clarify this question, we consider the effect of the crack length a and the applied shear stress Σ on the ratio L˜II,1L˜II,2. Evidently, this ratio will depend on the material parameters c1 and c3. The effect of a and Σ is illustrated in Figure 26ad. In all figures, the width b and the length 2h of the specimen are the same and equal to the values given in Section 4.1. However, a and Σ are different for each figure. It can be recognized from Figure 26a and b that the ratio L˜II,1L˜II,2 does not depend on the applied loading, Σ, which might be expected, for the problem is linear. But the ratio L˜II,1L˜II,2 depends on the crack length, a, as can be seen by comparing Figure 26a and c or Figure 26b and d. This result may be seen as a justification of considering L˜II,1 and L˜II,2 as different stress intensity factors.

## 5. Concluding remarks

A finite element model for the 3-PG Model has been developed and verified with reference to the expected responses of a square section with a hole subjected to displacement-controlled tension loading. Then, the finite element model has been employed to discuss Mode-I and Mode-II crack problems. From these investigations, one can draw the following conclusions.

1. The finite element and the analytical solutions fit with good accuracy.

2. In particular, the finite element simulations confirm the order of singularity, r12, and the assumed symmetry conditions for the micro-deformation very well. Further, they confirm the existence of a constant term in the asymptotic expansion of the micro-deformation.

3. By using numerical simulations with the developed finite element model, the effect of the nonclassical material parameters on the classical and the nonclassical stress intensity factors has been investigated. Especially, limiting cases of material parameters leading to responses of classical elasticity have been analyzed.

4. For Mode-II loading conditions, two independent stress intensity factors have been assumed to be present in the near-tip fields of the double stress. The numerical investigation seem to verify this remarkable feature.

## Acknowledgments

The authors would like to thank the TU Darmstadt for support of publishing this work in open access.

## Authors’ contributions

Each one of the authors, Broese, Frischmann, and Tsakmakis, contributed the same amount of work for the three chapters “Mode-I and Mode-II Crack Tip Fields in Implicit Gradient Elasticity Based on Laplacians of Stress and Strain,” Parts I–III.

Dr. C. Broese: theory and numerical simulations.

Dr. J. Frischmann: analytical solution and numerical simulations.

Prof. Dr. Tsakmakis: theory and analytical solution.

Therefore, it is a joint work by all three authors.

## Funding

The first and second authors acknowledge with thanks the Deutsche Forschungsgemeinschaft (DFG) for partial support of this work under Grant TS 29/13-1.

## Competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this chapter.

## Abbreviations

 eq. equation

## References

1. 1. Williams ML. On the stress distribution at the base of a stationary crack. Journal of Applied Mechanics. 1957;24:109-114
2. 2. Diegele E, Elsässer R, Tsakmakis C. Linear micropolar elastic crack-tip fields under mixed mode loading conditions. International Journal of Fracture. 2004;129:309-339
3. 3. Hughes TJR. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Englewood Cliffs, New Jersey: Prentice-Hall; 1987
4. 4. Gould GL. Introduction to Linear Elasticity. 2nd ed. New York: Springer-Verlag; 1994
5. 5. Henshell RD, Shaw KG. Crack tip finite elements are unnecessary. International Journal for Numerical Methods in Engineering. 1975;9:495-507
6. 6. Barsoum RS. On the use of isoparametric finite elements in linear fracture mechanics. International Journal for Numerical Methods in Engineering. 1976;10:25-37

Written By

Carsten Broese, Jan Frischmann and Charalampos Tsakmakis

Submitted: 10 July 2020 Reviewed: 17 August 2020 Published: 21 September 2020