Open access peer-reviewed chapter - ONLINE FIRST

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

By Carsten Broese, Jan Frischmann and Charalampos Tsakmakis

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

DOI: 10.5772/intechopen.93619

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 − 1 2 . 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 IIfor 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 IIas 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 Iapplies.

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 uixand Ψijxas well as the so-called virtual fields δujxand δΨjkxbelong 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 δujand Eq. (2) by δΨjkand 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 VPto indicate the summation of single integrals over VPi, which generally do not coincide. The meaning of the integration over VTis 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 Vis approximated by nelfinite 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 uiand Ψijare approximated by

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

where uiAand ΨijAare the unknown values of uiand Ψijat node A. The so-called shape functions NAuand NAΨbelong to finite-dimensional function spaces Shand Th, which approximate the function spaces Sand T, respectively. Similarly, δujand δΨjkare approximated by

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

with δujBand δΨjkBbeing constants. The functions δujhand δΨjkhare elements of finite-dimensional function spaces Vhand Wh, which approximate the function spaces Vand 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 VePand VeTbeing 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 Σijand σjkas well as the double stresses μijkin 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 δujBand δΨjkBmay 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 xiAare 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=±b2are assumed to be free of classical and nonclassical tractions. At the boundary y=b2, the macro-displacement component uy, the component Pxof the classical traction, and the nonclassical traction components Tijare assumed to vanish. At the boundary y=b2, the macro-displacement in the y-direction is given by uy=0.1mm, while Px=0and Tij=0are 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 Σyyalong the section A-A(see Figure 1a), as a function of the local coordinate

axr,E37

with xr. A so-called stress concentration factor kis 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-Ahas 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 16stripes of quadratic elements with reduced integration, that is, elements with 8nodes, is used. Every stripe consists of 180elements 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 720such stripes is used, whereas the mesh for the simulations of the 3-PG Model only needs 72such stripes of elements.

Figure 2 illustrates the distribution of the dimensionless stress ΣyyΣ0along the section A-A. The value of ΣyyΣ0at a=0represents the stress concentration factor kand turns out to be k=3,14for 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 c1and

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>0applies. 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Σ0along the section A-Ato be close to the classical one whenever c3is 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 c3on the values of the stress concentration factor k. We expect that for c1=const.and c30, the values of kwill 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 c1on 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 c1on kis illustrated in Figure 4 by the graph referred to as specimen 1. A convenient way to illustrate the effect of c1is 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,20and 200, respectively. The corresponding graphs of kas a function of c1are displayed in Figure 4 and are referred to as specimens 4, 20, and 200, respectively. Keeping in mind, that c1is an internal material length, we rescale the abscissa c1by considering the graphs of kas 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=11mmand 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 16stripes 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 Eand ν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=100MPato 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 pon 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=12with great accuracy. Therefore, the value p=12will be fixed, to avoid errors owing to inaccurate numerical determination of pwhen discussing predicted responses. All stress intensity factors are determined from Eqs. (49)(53) for the fixed value p=12by 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αβIand gαβγIfor 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=5106mmby K˜I2πr, r2πL˜Ic2c12μr2πL˜Ic1c32μand L˜I2πr, respectively. The stress intensity factors K˜Iand L˜Ihave been determined by least square fitting as described in the last section. In the finite element computations, the values c1=104mm2and c3=1for the nonclassical material parameters are chosen. The values of Ψ¯ijare 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Ψ¯ijIverify the existence of this constant terms. The fact that the analytical and the numerical results of the angular functions hαβIagree 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˜Iand K˜IIof the 3-PG Model will be, in general, different from the stress intensity factors KIand KIIof 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, c1and 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 c1or for c1=const.and very small values of c3, the values of K˜Iconverge 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 c1and c3on the nonclassical stress intensity factor L˜Iis illustrated in Figures 14 and 15. The principal observation is that for smaller values of c3, the value of L˜Igets smaller as well. On the other hand, if c3=const., then the values of L˜Iincrease with increasing values of c1.

4.4 Results for mode-II

The graphs of the angular functions fαβIIfor 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 Ψ¯ijhave been determined to be Ψ¯11=0, Ψ¯22=0, and Ψ¯12=0,00391and this, in turn, verifies the existence of the constant terms Ψ¯ijΨ¯ijII.

The effect of c1and c3on K˜IIis illustrated in Figures 20 and 21 and is similar to the effect of c1and c3on K˜I(cf. Figures 12 and 13).

The effect of c1and c3on the nonclassical stress intensity factors L˜II,1and L˜II,2is illustrated in Figures 2225. Again, this effect is quite similar to the effect of c1and c3on L˜I, but the stress intensity factors L˜II,1and L˜II,2both 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 Ψαβ0depend 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 B0and E0with 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 B0and E0, or equivalently between L˜II,1and L˜II,2, which is independent of the crack length. To clarify this question, we consider the effect of the crack length aand the applied shear stress Σon the ratio L˜II,1L˜II,2. Evidently, this ratio will depend on the material parameters c1and c3. The effect of aand Σis illustrated in Figure 26a–d. In all figures, the width band the length 2hof the specimen are the same and equal to the values given in Section 4.1. However, aand Σare different for each figure. It can be recognized from Figure 26a and b that the ratio L˜II,1L˜II,2does not depend on the applied loading, Σ, which might be expected, for the problem is linear. But the ratio L˜II,1L˜II,2depends 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,1and L˜II,2as 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

How to cite and reference

Cite this chapter Copy to clipboard

Carsten Broese, Jan Frischmann and Charalampos Tsakmakis (September 21st 2020). Mode-I and Mode-II Crack Tip Fields in Implicit Gradient Elasticity Based on Laplacians of Stress and Strain—Part III: Numerical Simulations [Online First], IntechOpen, DOI: 10.5772/intechopen.93619. Available from: