Ansys Fracture Tutorial

Universitat de Girona

FRACTURE MECHANICS Computer lab sessions

D. Trias

October 2012 This document can be found at: ftp://amade.udg.edu/amade/mme/MecFrac/MecFrac.pdf

Contents 1 Singular stresses 1.1

1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.1.1

Effect of element type, size and shape

. . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.1.2

Finite element discretisation of stresses at a crack tip . . . . . . . . . . . . . . . . . .

4

1.2

Quarter point / crack tip elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.3

Creating quarter mid-nodes at crack tip with ANSYS . . . . . . . . . . . . . . . . . . . . . .

7

1.3.1

Meshing with usual tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.3.2

Meshing with special tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

1.4

Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

1.5

Suggested problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

1.6

Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2 Computational Fracture Mechanics I: Computation of G

17

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.2

Finite Crack Extension Method (FCEM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

2.3

Crack Closure Method (CCM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

2.4

Virtual Crack Closure Technique (VCCT) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

2.5

Suggested exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.6

Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

3 Computational Fracture Mechanics II: Computation of K

27

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.2

The stress intensity factor (K) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.2.1

Numerical estimation of the stresses at the crack tip . . . . . . . . . . . . . . . . . .

28

3.2.2

Computation of K by stress extrapolation . . . . . . . . . . . . . . . . . . . . . . . . .

30

3.2.3

Computation of K by displacement extrapolation . . . . . . . . . . . . . . . . . . . . .

30

3.2.4

Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

Displacement extrapolation with quarter node elements . . . . . . . . . . . . . . . . . . . . .

32

3.3.1

Formulae for the stress intensity factor . . . . . . . . . . . . . . . . . . . . . . . . . .

32

ANSYS commands for the computation of K . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

3.4.1

Crack opening displacement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

3.4.2

KCALC command . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

3.3 3.4

iii

3.5

Proposed exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

3.6

Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4 Computational Fracture Mechanics III: Computation of the J-integral

37

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

4.2

The J integral with ANSYS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4.3


39

4.4

Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

5 Computational Fracture Mechanics IV: Cohesive zone modeling

41

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

5.2

Cohesive laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

5.2.1

Bilinear law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

5.2.2

Exponential law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

Cohesive elements in ANSYS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

5.3.1

Cohesive zone modeling with interface elements . . . . . . . . . . . . . . . . . . . . .

44

5.3.2

Cohesive zone modeling with contact elements . . . . . . . . . . . . . . . . . . . . . .

48

5.4

Some remarks on element size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

5.5


51

5.6

Further reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

5.7

Aknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

5.3

Chapter 1

Singular stresses 1.1

Introduction

Linear Elastic Fracture Mechanics deals with cracked solids. This means the stress and strain fields within this loaded solids are assumed to be influenced somehow by the presence of cracks. In order to analyze which is the efect of these cracks let us model a simple cracked plate. Example 1.1 (Analysis of a cracked plate). Let us assume we consider the analysis of the cracked plate component of Figure 1.1 by the Finite Element Method.

Figure 1.1: Plate with a central sharp crack.

For the sake of simplicity let us assume the material of the plate is steel and the applied pressure is 1 MPa. Solution to Example 1.1. The ANSYST M command sequence for this example is listed below. You can either type these commands on the command window, or you can type them on a file, then, on the command window enter /input, file, ext or just use copy and paste.

FINISH /CLEAR /TITLE Stress singularity

1

2

Fracture Mechanics

!PRE-PROCESSOR ******************************* /PREP7 !Geometrical parameters L = 50 ! Half length of component a = 10 ! Crack half length b = 0 ! Crack width th = 30 ! Component thickness el_len = 5 ! Element length (mesh) el_shape=0 ! 0: quad, 1: triangle ET,1,PLANE42 !4 node solid element !ET,1,PLANE82 !8 node solid element !Material properties MP,EX,1,210000 !Young modulus MP,PRXY,1,0.3 !Poisson’s coef. KEYOPT,1,3,3 !Keyoption for thickness activation R,1,th !Thickness for element type #1 !Geometry definition by keypoints K,1,a,0 K,2,L,0 K,3,L,L K,4,0,L K,5,0,b !Lines from keypoints L,1,2,(L-a)/el_len L,2,3,L/el_len L,3,4,L/el_len L,4,5,(L-b)/el_len L,5,1,a/el_len !Areas from lines AL,1,2,3,4,5 !Mesh LCCAT,5,1 MSHAPE, el_shape, 2D MSHKEY,2 AMESH,ALL FINISH !SOLUTION *********************************** /SOLU

Chapter 1. Singular stresses

DL,1,1,SYMM DL,4,1,SYMM SFL,3,PRES,-1 !Pressure on top tip SOLVE FINISH !POST-PROCESSOR ****************************

/POST1

PLDISP,1 !Deformed shape PLNSOL,S,Y ! Stresses Y direction !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Path plot of stress components at theta=0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !PATH,0DEG,2,10,50 !PPATH,1,,a,0,0 !PPATH,2,,3*a,0,0 !PDEF,SX,S,X,NOAV !PDEF,SY,S,Y,NOAV !PDEF,UY,U,Y,NOAV !PDEF,SY2,S,Y,AVG !PLPATH,SY2 !PLPATH,SX !PLPATH,UY !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Path plot of stress components at theta=45 deg !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! PATH,45DEG,2,10,50 PPATH,1,,a,0,0 PPATH,2,,3*a,3*a-a,0 PDEF,SX,S,X,NOAV PDEF,SY,S,Y,NOAV PDEF,UY,U,Y,NOAV PDEF,SY2,S,Y,AVG PLPATH,UY PLPATH,SX PLPATH,SY PRPATH,SY FINISH

This file can be found at: ftp://amade.udg.edu/mme/FEmet/1_stress_sing.dat

3

4

Fracture Mechanics

1.1.1

Effect of element type, size and shape

Let us consider now the effect of element size, type and shape in the stress field of the analyzed plate. To do so, let us fill in the following table: Element size

Triangle

Triangle

Quad

Quad

(el _l en )

PLANE42

PLANE82

PLANE42

PLANE42

1 mm 2 mm 3 mm 4 mm

1.1.2

Finite element discretisation of stresses at a crack tip

Try to explain what happens when regular elements are employed to approximate the stress field at a crack tip. You may employ Fig 1.2 by drawing on it the stress field, different discretization sizes and the approximation obtained by each one. You may start thinking what happens when constant strain elements are used and then try to figure out which is the solution given by linear elements.

Figure 1.2: Finite element discretization of the stresses near the crack tip

As you may have noticed, traditional engineering maximum stress-based failure analysis become senseless when stress singularities appear, so other approaches must be used. These are basically stress intensity factor (K) based approaches and energy release rate (G) based approaches. While G-based and K-based approaches are the main topic of next chapters, the next section gives an introduction of quarter node tip elements, which are used to capture efficiently the stress singularity in a finite element discretization.

5


1.2

Quarter point / crack tip elements p

The aim of these elements is to introduce in the element formulation a stress singularity of the 1/ r type. This is useful for stress related approaches, such as the numerical computation of the stress intensity factor. To do so we may start by recalling the isoparametric formulation of a 1D quadratic Lagrangian element:

Figure 1.3: 8 node isoparametric quadrilateral element

N1 N2 N3

1 = − ξ(1 − ξ) 2 = (1 − ξ2 ) 1 = ξ(1 + ξ) 2

(1.1) (1.2) (1.3)

Since for an isoparametric element the same approximation for the geometry and for the displacements is used, the geometry of the 1-3 edge may be expressed: x=

n=3 X i

1 1 Ni x i = − ξ(1 − ξ)x 1 + (1 − ξ2 )x 2 + ξ(1 + ξ)x 3 2 2

Figure 1.4: 8 node isoparametric element with quarter-side located mid-nodes

(1.4)

6

Fracture Mechanics

Then if we locate the origin of a quad element of side length L at node 1 and locate the mid-node (node #2) at x 2 = L/4, as shown in Figure 1.4: 1 L x = ξ(1 + ξ)L + (1 − ξ2 ) 2 4

(1.5)

And solving for ξ: r ξ = −1 + 2

x L

(1.6)

As you remember, the displacement approximation is given by: u=

n=3 X i

1 1 Ni u i = − ξ(1 − ξ)u 1 + (1 − ξ2 )u 2 + ξ(1 + ξ)u 3 2 2

(1.7)

where u 1 , u 2 and u 3 are the displacements at nodes 1,2 and 3. Using equation 1.6 in the former equation we obtain the expression of the displacements as a function of the geometry (x): r ¶µ r ¶ r ¶µ r ¶ µ µr ¶ µ 1 x x x x x x 1 u = − −1 + 2 − 2−2 u1 + 4 u 2 + −1 + 2 2 u3 2 L L L L 2 L L

(1.8)

If we now compute the strain in the x direction: µ ¶ µ ¶ µ ¶ ∂u ∂ξ ∂u 1 3 4 2 4 4 1 1 = =− p − u1 + p − u1 + − p + u3 (1.9) ∂x ∂x ∂ξ 2 2 xL L xL L xL L p We may easily verify that the former expression is a function of (1/ x) and, consequently, the strain εx =

field presents a singularity. However, for a 2D element we might apply the same method for the 8-node serendipity element of Figure 1.4. This way, the singularity would be present only along the edge 1-3, and εr would not be singular. We may obtain a radial singular stress field by collapsing nodes 4, 5 and 6 of the quad of Figure 1.4, ash shown in Figure 1.5.

Figure 1.5: Quarter tip element from collapsed 8 node quadrilateral


1.3

7

Creating quarter mid-nodes at crack tip with ANSYS

As seen in the former sections, if we have to deal with the stresses when stress singularities are present it is useful to use some kind of special elements. We also showed how, actually, these special elements may be constructed from regular elements just by placing the mid-side nodes at a quarter of the element side length of the crack tip. We usually may do this by our means or by using some special commands provided by the commercial finite element software.

1.3.1

Meshing with usual tools

We may construct this kind of mesh at a crack tip manually, that is, placing the corresponding mid-side nodes at a distance of a quarter of the side from the crack tip. The following ANSYS log file analyzes the same problem of Figure 1.1, using N, NGEN and E commands to build the mesh. Example 1.2. Model the cracked plate of Figure 1.1 using the isoparametric quarter node elements presented in Section 1.2. Locate the quarter nodes by using ANSYS’ common node and element generation commands. Solution to Example 1.2. The ANSYST M command sequence for this example is listed below. You can either type these commands on the command window, or you can type them on a file, then, on the command window enter /input, file, ext or just use copy and paste.

FINISH /CLEAR /TITLE Stress singularity - Mesh #2 !=================================== !PRE-PROCESSOR !----------------------------------/PREP7 ! Geometrical parameters in mm L = 5 ! Plate length a = 1 ! Crack length b = 0 ! Crack heigth th = 1 ! Plate thickness el_len = 0.2 ! Element length pi= 3.1415926535897932384626433832795 fx = a/10 rdiv = 4 ! radial divisions at crack tip (must be 4) tdiv = 12 ! angular divisions at crack tip rdiv2 = 10 ET,1,PLANE42

8

Fracture Mechanics

!Material properties MP,EX,1,210000 !Young modulus MP,PRXY,1,0.3 !Poisson’s ratio KEYOPT,1,3,3 !Keyoption to introduce thickness in element 42 R,1,th ! Thickness !Keypoints which define the geometry K,1,a,0 K,2,L,0 K,3,L,L K,4,0,L K,5,0,b K,6,a,L K,7,2*a,0 K,8,2*a,a K,9,a,a K,10,0,a K,11,a+2*fx,0 K,12,a,2*fx K,13,a-2*fx,0 ! Lines L,1,11,100,1000 !1 L,11,7,rdiv2 !2 L,7,2,10 !3 L,2,3,tdiv/4 !4 L,3,6,tdiv/4 !5 L,4,10,10 !6 L,10,5,tdiv/4 !7 L,10,9,tdiv/4 !8 L,9,8,tdiv/4 !9 L,7,8,tdiv/4 !10 L,5,13,rdiv2 !11 L,9,12,rdiv2 !12 !L,1,13,85,20 !13 L,1,13,100,1000 !13 LARC, 11, 12, 1, 2*fx !14 LARC, 12, 13, 1, 2*fx !15 LESIZE,14,,,tdiv/2 LESIZE,15,,,tdiv/2 L,1,12,100,1000 !16 L,6,4,tdiv/4 !17 L,9,6,10 !18 ! Nodes

CLOCAL,11,1,a,0,0 CSYS,11


N,1,0,0,0 N,rdiv,2*fx,0,0 FILL,1,,1,2,1,,,1 NGEN,tdiv*2+1,10,2,rdiv,1,0,90/tdiv,0,0 csys,1 ! Generate 1/4 nodes *DO,I,1,tdiv+1 FILL,1,2+20*(I-1),1,1000+I,,,,3 *ENDDO

*DO,I,1,tdiv+1 FILL,2+20*(I-1),4+20*(I-1),1,2000+20*(I-1),,,,1 *ENDDO

csys,1 ET,2,PLANE82 KEYOPT,2,3,3 !Keyoption to introduce thickness in element 82 R,2,th !thickness TYPE,2 ! Elements csys,11 *DO,I,1,tdiv E,2+20*(I-1),2+20*I,1,1,12+(I-1)*20,1001+I,1,1000+I *ENDDO E,2,4,24,22,2000,14,2020,12 EGEN,tdiv,20,tdiv+1

AL,2,10,9,12,14 ! Area created from lines AL,11,15,12,8,7 AL,3,4,5,18,9,10 AL,18,17,6,8 AL,1,14,16 AL,16,15,13 ASBA,1,5,,DELETE,DELETE ASBA,2,6,,DELETE,DELETE LCCAT,4,5 !19 LCCAT,9,10 !20 LCCAT,7,8 !21 TYPE,2 !MSHKEY,2 AMESH,1,4 !Mesh csys,11

9

10

Fracture Mechanics

NSEL,S,LOC,X,2*fx-fx/50,2*fx+fx/50 NUMMRG,NODE,0.001 theta=pi/(2*tdiv) !*DO,I,1,tdiv !NMODIF,12+20*(I-1),fx*cos(theta) !NMODIF,14+20*(I-1),2*fx*cos(theta) !*ENDDO /PSYMB,ESYS,1 csys,0 ALLSEL FINISH

/SOLU ANTYPE,0 D,1,UY,0 D,1,ROTZ,0 D,2,UY,0 D,2,ROTZ,0 DL,1,5,SYMM DL,2,1,SYMM DL,3,3,SYMM DL,6,4,SYMM DL,7,2,SYMM SFL,5,PRES,-1 !Pressure applied to upper line SFL,17,PRES,-1 ALLSEL SOLVE FINISH !=================================== !POST-PROCESSOR !----------------------------------/POST1 PLDISP,1 !Deformed shape PLNSOL,S,Y ! Stresses Y direction angle=pi/6 PATH,30DEG,2,10,100 PPATH,1,,a,0,0 PPATH,2,,a+cos(angle),sin(angle),0 PDEF,SX,S,X,NOAV

11


PDEF,SY,S,Y PDEF,UY,U,Y,NOAV PLPATH,SY PLPATH,SX PRPATH FINISH This file can be found at: ftp://amade.udg.edu/mme/FEmet/1_quarter_mid_nodes.dat

Figure 1.6 shows the different path plots for different elements sizes when using ANSYS’ PLANE42 and PLANE82 elements compared to the solution obtained with 2 mm Quarter mid-side node elements which were constructed the way described above. Results show as these elements are able to provide a solution p

similar to that obtained by much smaller elements. Since these elements include the 1/ r singularity, they do not provide a value for r = 0, which is the singular point.

7

6

0.5 mm PLANE42 TRI 2 mm PLANE42 TRI 5 mm PLANE42 TRI Quarter mid−side nodes σyy (MPa)

σyy (MPa)

5 4 3

5 4 3

2 1 0

2 2

4

6

r0 (mm)

6

8

3 2

2

4

r0 (mm)

6

6

8

10

1 mm PLANE82 QUAD 2 mm PLANE82 QUAD 5 mm PLANE82 QUAD Quarter mid−side nodes

5 σyy (MPa)

4

1 0

1 0

10

0.5 mm PLANE42 QUAD 2 mm PLANE42 QUAD 5 mm PLANE42 QUAD Quarter mid−side nodes

5 σyy (MPa)

1 mm PLANE82 TRI 2 mm PLANE82 TRI 5 mm PLANE82 TRI Quarter mid−side nodes

6

4 3 2

2

4

r0 (mm)

6

8

10

1 0

2

4

r (mm)

6

8

10

0

Figure 1.6: σ y y at the crack tip region for the analyzed plate, as obtained with different element type/size

and with Quarter mid-side node elements (45o path)

1.3.2

Meshing with special tools

Since doing this may be tricky, commercial finite element software usually have special tools for meshing the crack tip. The next box shows ANSYS command for crack tip meshing.

12

Fracture Mechanics

ANSYS Command: Meshing the crack tip

KSCON, NPT, DELR, KCTIP, NTHET, RRAT NPT number of the node located at the crack tip DELR Radius of first row of elements about crack tip KCTIP Crack tip singularity key. For our purposes its value should be 1 NTHET Number of elements in circumferential direction. Default is one per 30o . RRAT Ratio of 2nd row element size to DELR. Default is 0.5

The following ANSYS log file performs the same analysis but now the command KSCON is employed to mesh the crack tip. Example 1.3. Model the plate of Figure 1.1 using quarter node elements for the crack tip obtained with the KSCON command. Solution to Example 1.3. The ANSYST M command sequence for this example is listed below. You can either type these commands on the command window, or you can type them on a file, then, on the command window enter /input, file, ext or just use copy and paste.

FINISH /CLEAR /TITLE Stress singularities - KSCON /PREP7 ! Geometrical parameters in mm L = 50 ! Plate length a = 10 ! Crack length b = 0 ! Crack heigth th = 30 ! Thickness el_len = 2 ! Element length ET,1,PLANE82 !Material properties MP,EX,1,210000 ! Young modulus MP,PRXY,1,0.3 ! Poisson’s ratio KEYOPT,1,3,3 ! Keyoption to introduce thickness


R,1,th ! Thickness K,1,a,0 K,2,L,0 K,3,L,L K,4,0,L K,5,0,b K,6,2*a,0 K,7,a,a K,8,0,a K,9,a+sqrt(2)*a/2,sqrt(2)*a/2 K,10,a,L L,1,6,(2*a)/el_len L,6,2,(L-2*a)/el_len L,2,3,L/el_len L,3,10,(L-a)/el_len L,10,4,a/el_len L,4,8,(L-a)/el_len L,8,5,a/el_len L,5,1,a/el_len L,9,3,sqrt((L-a)**2+L*L)/el_len L,1,7,a/el_len L,7,8,a/el_len LARC,7,9,1,a,6 LARC,9,6,1,a,6 L,7,10,(L-a)/el_len KSCON,1,0.15,1,12,0.25 AL,2,3,9,13 AL,12,9,4,14 AL,1,13,12,10 AL,8,10,11,7 AL,11,14,5,6 AMESH,ALL FINISH /SOLU

DL,6,5,SYMM DL,7,4,SYMM DL,2,1,SYMM DL,1,3,SYMM SFL,4,PRES,-1 SFL,5,PRES,-1 SOLVE FINISH /POST1

13

14

Fracture Mechanics

PLDISP,1 PLNSOL,S,Y FINISH This file can be found at: ftp://amade.udg.edu/mme/FEmet/1_stress_sing_KSCON.dat

Figure 1.7 shows pathplots for σ y y as obtained with quarter mid-side node elements. Two of the plots correspond to elements obtained with the KSCON command. The third one corresponds to the solution given by the elements obtained in section 1.3.1. 30 KSCON1 KSCON2 modified mesh

σyy (MPa)

25 20 15 10 5 0 0

0.5

1

1.5 r0 (mm)

2

2.5

3

Figure 1.7: σ y y at the crack tip region for the analyzed plate, as obtained with with Quarter mid-side node

elements: several mesh refinement with KSCON and by-hand mesh modification (45o path)

1.4

Summary and conclusions

• When stress singularities are present in some kind of problem, the finite element analysis and specially the mesh must be carefully considered. • If regular elements are employed, a very fine discretization is needed • Special crack tip elements may be used to estimate adequately the stress singularity with a lower number of elements. • However the stress solution should still be used carefully, knowing that any solution with a finer mesh at the crack tip will lead to a higher value of the maximum value of the stress (which analytical value is ∞) • Energy release rate (G) - based analysis or stress intensity factor (K) - based analysis are a much better option for linear elastic fracture mechanics regime.


1.5

15

Suggested problems

Problem 1.1. Let us assume we want to model the bi-metal shown in Figure 1.1. Consider that the bi-metal thickness is 30 mm. Apply the displacement on the right side of the steel block. Material data: E st eel = 210000 MPa, νst eel = 0.3, E Al = 70000 MPa, ν Al = 0.3.

• Solve the problem for three different meshes and analyze the stress in the vertical direction. • Identify the zone which reaches higher stresses. What happens with the maximum value (absolute value) of the stress when remeshing? Why? • Construct an adequate mesh using KSCON command and check that, at least for a range of some element size at the crack tip, solution is not mesh dependent.

Figure 1.8: Aluminium-Steel Bi-metal for Suggested Problem 1.1

1.6

Further reading

Henshell R.D and Shaw K.G. (1975) Crack tip elements are unnecessary. International Journal for Numerical Methods in Engineering 9(3): 495-507. Saouma V.E. and Schwemmer D. (1984) Numerical evaluation of the quarter-point crack tip element. International Journal for Numerical Methods in Engineering 20(9): 1629-1641. Gray L.J., Phan AV. et al (2003) Improved quarter point crack tip element.Engineering Fracture Mechanics 70:269-283

Chapter 2

Computational Fracture Mechanics I: Computation of G 2.1

Introduction

As seen in the former chapter, the traditional-materials strength stress analysis of a cracked component may be hardly tackled. Although the stress discretisation may be improved by using crack tip elements, the meaningful analysis is generally that performed using the energy release rate (G). The energy release rate may be defined as the rate at which energy is dissipated (∆Π) when a crack grows, under constant boundary conditions:

¯ ∆Π ¯¯ G =− ∆A ¯constant B.C

(2.1)

As seen, the crack opening is measured in terms of the created area (∆A ). This is possible because G is a state function which means it only depends on the updated geometry and the geometry but not on how they change in the fracture process. So no matter if the crack grows we can compute G using different crack lengths. In the following sections, differents ways for the numerical computation of G, according to its definition of equation 2.1 will be summarized. We will apply these methods to the classical example of a cracked plate shown in figure 3.2.1, using an initial crack length a=10mm, L=50mm, unit pressure and a plate thickness of 30 mm. Although the crack length should be much more longer to satisfy the main assumptions involved, you may compare the results we will obtain with the Classical Beam Theory approximation, in order to get a rough idea of the magnitude: GC BT =

P 2 a 2 12P 2 a 2 = 2 3 BE I B h E

17

(2.2)

18

Fracture Mechanics

Figure 2.1: Plate with a side sharp crack. G CBT =

N /mm

Recall that equation 2.2 is an approximation and the analyzed geometry does not satisfy the basic assumpttions of the Classic Beam Theory, so the result for sure includes large errors. You may use it only to check the order of magnitude of the computed energy release rate in the following sections.

2.2

Finite Crack Extension Method (FCEM)

In Linear Elastic Fracture Mechanics (LFEM), in quasi-static conditions (which means the kinetic energy involved in the process may be neglected), the elastic stored energy (Π) equals the difference between the strain energy (U ) and the external work (W ):

Π =U −W

(2.3)

where the strain energy U : U=

1 · σi j εi j 2

(2.4)

According to the definition of G, it may be computed:

G=

Π(a + δa) − Π(a) ∆A

(2.5)

where, if B is the thickness of the component, ∆A = δa · B . By these assumption, G may be computed using two different finite element models with two different crack lengths and the same boundary conditions: • Finite Element Model 1. BC and crack length a • Finite Element Model 2. Same BC but crack length a + δa

Chapter 2. Computational Fracture Mechanics I: Computation of G

19

ANSYS Command: Computation of Strain Energy 1. Using User GUI: • General Post-processor → Element Table → Define Table → Add • General Post-processor → Element Table → Sum of each item 2. Using ANSYS commands: • AVPRIN • ETABLE • SSUM

Example 2.1. Compute the energy release rate (G) for the side-cracked plate of Figure 3.2.1 by means of the Finite Crack Extension Method. You may use the parametrized ANSYS model given below. Solution to Example 2.1. The ANSYST M command sequence for this example is listed below. You can either type these commands on the command window, or you can type them on a file, then, on the command window enter /input, file, ext or just use copy and paste. FINISH /CLEAR /TITLE Computation of G - Crack length: a /PREP7 !PRE-PROCESSOR L = 50 !Length of component (mm) a = 10 ! Crack length (mm) b = 0 ! Crack heigth (mm) th = 30 !Thickness (mm) el_len = 2 ! Element length ET,1,PLANE42 KEYOPT,1,3,3 !Keyoption to activate thickness R,1,th !Thickness assignment MP,EX,1,210000 !Young modulus MP,PRXY,1,0.3 !Poisson’s ratio K,1,a,0 K,2,L,0 K,3,L,L K,4,0,L K,5,0,b L,1,2,(L-a)/el_len L,2,3,L/el_len L,3,4,L/el_len

20

Fracture Mechanics

L,4,5,(L-b)/el_len L,5,1,a/el_len AL,1,2,3,4,5 LCCAT,5,1 ARSYM,Y,1 LCCAT,7,8 MSHKEY,2 AMESH,ALL ALLSEL NSEL,S,LOC,X,a,L ! All nodes at y=0, but those of the crack NSEL,R,LOC,Y,0

! are selected

CPINTF,ALL FINISH /SOLU !SOLUTION DL,10,2,ALL,0 SFL,3,PRES,-1

!Pressure

nsel,s,loc,y,L CP,1,UY,ALL allsel SOLVE FINISH /POST1 !GENERAL POST-PROCESSOR PLDISP,1 PLNSOL,S,Y FINISH

This file can be found at: ftp://amade.udg.edu/mme/FEmet/2_comp_g.dat

You may fill in the following table to compute the energy release rate according to the Finite Crack Extension Method: Magnitude

FE model 1 (a)

FE model 2 (a + δa)

Strain energy Displacement Force W Π

G IFCEM =

N /mm

21


2.3

Crack Closure Method (CCM)

The Crack Closure Method assumes that the energy which is dissipated when a crack grows some δa is the same energy needed to open the crack some δa . So if, as in the former method, two different finite element models are employed, the forces and displacements needed to close the crack by some δa may be computed as follows. • Finite Element Model 1 (Crack length = a ). Used to obtain the force acting at the crack tip. Since this force is actually a reaction, some rigid link or Multi Point Constraint (MPC) should be employed at the crack tip to obtain this force value. • Finite Element Model 2 (Crack length = a + δa ). The longer crack length is achieved by removing the rigid link at the former crack tip and keeping a rigid link at the new one. This model will be used to read the displacements needed to close the crack by δa . The Energy Release Rate (G) may be computed: h i (1) (2) G = F x(1) · δ(2) + F · δ x y y

1 2 · δa · B

(2.6)

where the first term of the addition corresponds to G I I and the second term to G I , that is: GI =

(1) (2) 1 F y · δy 2 δa · B

(2.7)

GI I =

1 F x(1) · δ(2) x 2 δa · B

(2.8)

(2) 0(2) (2) (2) 0(2) and δ(2) and superscripts (1) and (2) denote the model where the variable x = u x −u x , δ y = u y −u y

is taken.

Figure 2.2: Closure Crack Method

ANSYS Command: Removing CPs • CPINTF • NUMMRG

22

Fracture Mechanics

ANSYS Command: Reading displacement and force results • PRNLD • PRRSOL

You may fill in the following table to compute G by using the nodal values of force and displacement, according to the CCM: Magnitude

FE model 1 (a)

Displacement (top)

-

Displacement (bottom)

-

δy

-

FE model 2 (a + δa)

Force (F y ) G ICCM =

2.4

N /mm

Virtual Crack Closure Technique (VCCT)

The Virtual Crack Closure Technique makes a further assumption: the crack grows in a self-similar manner. This means that if we only look at the nearby of the crack tip, from one growth step to the next one, we would see about the same crack shape -the same displacements- and about the same forces acting at the crack tip. Consequently, instead of using two different models to get the forces at the crack tip and the displacement needed to close the crack by some δa , we may use the same model, so the computational efforts are reduced.

Figure 2.3: Virtual Crack Closure Technique

So, now G may be computed as:

23


£ ¤ G = F x · δx + F y · δ y

1 2 · δa · B

(2.9)

where, again, the first term of the addition corresponds to G I and the second term to G I I , that is: GI =

1 F x · δx 2 δa · B

(2.10)

GI I =

1 F y · δy 2 δa · B

(2.11)

where and δx = u x − u x0 , δ y = u y − u 0y . Finally, you may compute G through the formula derived by the Virtual Crack Closure Technique: Magnitude

FE model (a + δa)

Displacement (top) Displacement (bottom) δy

Force (F y ) G IVCCT =

2.5

N /mm

Suggested exercises

Problem 2.1. Model the CT specimen of Figure 2.4 using ANSYS with W=5mm, A=25 mm, B (thickness)=10 mm, C=50 mm, D = 62.5 mm, E=5mm and F=60 mm.

Figure 2.4: CT specimen

With this model you are going to analyze G as a function of the crack length (a), using different methods and for two different load cases: constant force (P=1000 N) and constant displacement (0.053mm). To make this you have to use in the ANSYS model Constraint Equations in the line where the crack will grow.

24

Fracture Mechanics

The crack should grow between 2 mm and 16mm, so a element size of 2 mm may be a good choice. You may use the provided file CT.dat. In each model you have to keep the following data: • External load or applied displacement • In the constant force loadcase, the displacement of the node where the force is being applied. • In the constant displacement loadcase, the reaction at the node where the displacement is being applied. • Force at the crack tip, in the vertical direction. • Displacement of the nodes closer to the crack tip, in the vertical direction. • Strain energy Compute: 1. For both loadcases, get the compliance curve for the specimen as a function of the crack length (C = f (a)). Derive numerically the obtained curve and use the computed derivative to compute G=

P 2 dC 2B da

2. G(a) for both loadcase using the FCEM, CCM, VCCT methods. Plot in the same graph the obtained curves together with the curve in the former question. (The four methods should give similar results, for the following, use only the curve obtained with FCEM.) 3. Assume that the material R-curve is given by:   G £1 − (1 − δa/c )3 ¤ for c f R=  Gc for

δa ≤ c f δa > c f

where c f = 1.25 mm and G c = 0.3N /mm . Find what load P produces instability. 4. Assume you are performing a laboratory test using a CT speciment, trying to measure the R-curve of a material. Which loadcase would you use? Why? 5. Analyze the effect of the element size on FCEM, CCM, VCCT methods. Obtain G for three different meshes with different element sizes using the three methods.

FINISH /CLEAR /TITLE Computation of G - Crack length: a /PREP7 !PRE-PROCESSOR a0=8


A= 25 W= 5 B=10 C=50 D=62.5 E=5 F=60

el_len = 2 ! Element length ET,1,PLANE42 KEYOPT,1,3,3 !Keyoption to activate thickness R,1,B !Thickness assignment MP,EX,1,210000 !Young modulus MP,PRXY,1,0.3 !Poisson’s ratio

K,1,0,0 K,2,C-A,0 K,3,C-A+E,w/2 K,4,D,w/2 K,5,D,F/3 K,6,D,F/2 K,7,C,F/2 K,8,0,F/2 K,9,C,F/3

L,1,2,(D-(D-C)-A)/el_len L,2,3,sqrt(E*E+W*W/4)/el_len L,3,4,(D-C+A-E)/el_len L,4,5,(F/3-W/2)/el_len L,5,6,(F/2-F/3)/el_len L,6,7,(D-C)/el_len L,7,8,C/el_len L,8,1,(F/2)/el_len L,5,9,(D-C)/el_len L,9,7,(F/2-F/3)/el_len AL,1,2,3,4,9,10,7,8 AL,9,5,6,10 LSYMM,Y,ALL AL,11,12,13,14,19,20,17,18 AL,19,15,16,20 ACCAT,1,2 ACCAT,3,4 MSHKEY,2 AMESH,5 AMESH,6 ALLSEL NSEL,S,LOC,X,0,C-A-a0 ! All nodes at y=0, but those of the crack

25

26

Fracture Mechanics

NSEL,R,LOC,Y,0

! are selected

CPINTF,ALL FINISH /SOLU !SOLUTION KD,9,UX,0 KF,9,FY,1000 ! Comment for displacement loadcase !KD,9,UY,10 ! Uncomment for displacement loadcase KD,18,ALL

ALLSEL SBCTRAN SOLVE FINISH /POST1 !GENERAL POST-PROCESSOR PLDISP,1 PLNSOL,S,Y FINISH

This file can be found at: ftp: // amade. udg. edu/ mme/ FEmet/ CT. dat

2.6

Further reading

Krüger R. (2002) The Virtual Crack Closure Technique: History, Approach and Applications. NASA/CR2002-211628. ICASE. Report No. 2002-10.

Chapter 3

Computational Fracture Mechanics II: Computation of K 3.1

Introduction

Although the use of the energy release rate is normally preferred in advanced analysis and in crack propagation simulation, the stress intensity factor K is widely used for design and verification of structures. While G is a energy-based magnitude, K is a stress related value and so, any computational method used to compute it will have to deal somehow with the stress singularity and its related issues we introduced in Chapter 1. As we saw in the first chapter, the stress discretisation in a finite element mesh may be improved by using crack tip elements, to avoid the strong mesh dependence produced by the stress singularity. In this chapter we will show how to use quarter mid-side node elements to discretize the stress field and then, some method to compute K.

3.2

The stress intensity factor (K)

The plane stress field in the nearby of a crack tip of a crack loaded in mode I can be approximated by the following expressions:

σxI

=

σIy

=

σxI y

=

µ ¶µ µ ¶ µ ¶¶ θ θ 3θ KI cos 1 − sin sin p 2 2 2 2πr µ ¶µ µ ¶ µ ¶¶ KI θ θ 3θ cos 1 + sin sin p 2 2 2 2πr µ ¶ µ ¶ µ ¶ KI θ 3θ θ sin cos cos p 2 2 2 2πr

(3.1) (3.2) (3.3)

were superscript I denotes mode I and θ and r are the polar coordinates (angle and distance, respectively) in a polar coordinate system with center at the crack tip. Analogously, the stress field ahead the crack tip in a Mode II situation is given by: 27

28

Fracture Mechanics

σxI I σIyI σxI Iy

µ ¶µ µ ¶ µ ¶¶ θ θ 3θ KI I = −p sin 2 + cos cos 2 2 2 2πr µ ¶ µ ¶ µ ¶ θ KI I θ 3θ sin = p cos cos 2 2 2 2πr µ ¶ µ ¶¶ µ ¶µ KI I θ θ 3θ = p cos 1 − sin sin 2 2 2 2πr

(3.4) (3.5) (3.6)

The stress at the crack tip might be seen as the limit: ¯ KI ¯ lim σiI j = p f iIj (θ)¯ θ=c r →0 2πr

(3.7)

where c denotes a constant. So if we are able to somehow know the stress field at the nearby of the crack tip we are able to compute K I : K I = lim σiI j r →0

p 2πr ¯ ¯ I f i j (θ)¯

(3.8)

θ=c

since

f iIj (θ)

are known trigonometrical functions. Analogously for K I I : K I I = lim σiI Ij r →0

p 2πr ¯ ¯ f iIjI (θ)¯

(3.9)

θ=c

3.2.1

Numerical estimation of the stresses at the crack tip

Let us recall the cracked plate of Figure that we analyzed in the former chapter, using the file 1_quarter_mid_nodes.d Using the ANSYS commands of the next box, you may obtain the stresses at a given path (radius). ANSYS Command: Path Plots

Command

Comments

PATH,NAME,nPTS,nSETS,nDIV

Defines geometrically a path by nPTS

PPATH,POINT,NODE,X,Y,Z,CS

Defines one of the points of the path. POINT is the ID of the point. NODE is a node number if the point is located in a NODE. X,Y,Z may be used to define the location of the point

PDEF,LABEL,ITEM,COMP,AVGLAB

Defines the ITEM (for instance STRESS) and COMP (for instance X) to plot and give it a LABEL

PLPATH,NAME

Plots the path labelled with NAME

Chapter 3. Computational Fracture Mechanics II: Computation of K

29

Example 3.1. Obtain a plot for σ y for the side-cracked plate, using a path plot. Solution to Example 3.1. The ANSYST M command sequence for this example is listed below. You can either type these commands on the command window, or you can type them on a file, then, on the command window enter /input, file, ext or just use copy and paste.

FINISH /POST1 ! Path plot of stresses PATH,0DEG,2,6,100 PPATH,1,,a,0,0 PPATH,2,,a+5,0,0 PDEF,SX,S,X,NOAV PDEF,SY,S,Y,NOAV PLPATH,SY PRPATH,SX,SY !Path results in a text file This file can be found at: ftp://amade.udg.edu/mme/FEmet/3_path.dat

A plot similar to that in Figure 3.1 should be obtained.

Figure 3.1: σ y at the nearby of the crack for θ = 0o .

30

Fracture Mechanics

3.2.2

Computation of K by stress extrapolation

Since we know how to obtain the stresses in the nearby of the crack tip we are able to obtain K I . If the stresses when r → 0 would not tend to ∞ we could compute K I with the stresses when r = 0, but since they do tend to ∞, we have to somehow compute numerically the limit of Eq. 3.8. To do so we compute K for each value of the stress σi j using equations as a function of r and plot the pairs (K,r). Since the values of the stress for small r are affected by the stress singularity we will neglect them and fit the linear variation of σi j (r ). The extrapolation for r = 0 gives a good approximation of K I , as shown in Figure 3.2 12

KI (MPa ⋅ mm

1/2

)

11

y = 0.183*x + 6.5

10 9 8 7 6 0

5

10

15

20

25

r (mm)

Figure 3.2: Computation of K I by stress extrapolation (from σ y

3.2.3

Computation of K by displacement extrapolation

The former procedure is strongly affected by the stress singularity. In a finite element procedure, stresses are obtained from the displacements and so may contain larger errors, specially in cases like this one where large stress gradients are present. For this reason, a more precise option is to use the displacement solution for the computation of the stress intensity factor, K . In this case for the region near the crack tip, the relations between the displacement field and K I and K I I are:     r  cos θ (κ − cos θ)   uI  2π 2 KI = 2G  sin θ (κ − cos θ)  r  vI  2

(3.10)

    r  sin θ (2 + κ + cos θ)   uI I  2π 2 = 2G KI I  cos θ (2 − κ − cos θ)  r  vII  2

(3.11)

where G is the shear modulus: G=

E 2(1 + ν)

(3.12)

and κ is a parameter which allows the simultaneous consideration of plane stress and plane strain cases, with:

31


3−ν 1+ν κ = 3 − 4ν κ =

for plane stress

(3.13)

for plane strain

(3.14)

Since in a finite element solution, the displacement field is generally a better solution than the stress field, the value of K I obtained in this manner (see Figure 3.4) should provide a better approximation.

9

KI (MPa ⋅ mm1/2)

8.5 8 7.5 y = 0.449*x + 6.78

7 6.5 6 0

0.5

1

1.5

2 r (mm)

2.5

3

3.5

4

Figure 3.3: Displacement extrapolation technique for the computation of K I (using u y ).

You may compare the numerical result with the analytical one, which may be obtained using the handbook formula of Figure 3.4.

Figure 3.4: Handbook expression for the analyzed case.

3.2.4

Remarks

• The reviewed techniques for the computation of K are first approximations. Further developments exist and are still object of current research. • Both stress and displacement extrapolation need of fine meshes to converge to the correct value of K. Crack tip elements are strongly recommended for the stress extrapolation method.

32

Fracture Mechanics

3.3

Displacement extrapolation with quarter node elements

When the stress singularity is very well discretized, which practically means when quarter-point isoparametric elements are used, some simple formulae can be applied with surprisingly accurate results. This formuale are obtained making θ = π in expressions 3.10 and 3.11, since for this angle the error is minimum. If you recall the quarter node elements shown in Section 1.2, the approximation of the displacement along edge 1-3 is given by: r u = u 1 + [4u 2 − u 3 − 3u 1 ]

r r + [2u 3 + 2u 1 − 4u 2 ] L L

(3.15)

and the same expression is valid for the vertical displacement v . If we substitute this approximation of the displacement field in expressions 3.10 and 3.11 we obtain:     r  cos θ (κ − cos θ)   4u 2 − u 3 − 3u 1  2π 2 KI = 4G  sin θ (κ − cos θ)  L  4v 2 − v 3 − 3v 1  2

(3.16)

    r  sin θ (2 + κ + cos θ)   4u − u − 3u  2π 2 3 1 2 KI I = 4G  cos θ (2 − κ − cos θ)  r  4v 2 − v 3 − 3v 1  2

(3.17)

and for mode II:

So we can compute K I and K I I substituting any value of angle θ and using the displacements at nodes 1,2 and 3. If we particularize the former expressions for θ = π and since v 1 = 0: 2G KI = κ+1

KI I

r

2G = κ+1

2π (4v 2 − v 3 ) L

r

2π (4u 2 − u 3 ) L

(3.18)

(3.19)

If we denote node 2 with an A and node 3 with a B and make L=` the former expressions may be written in the form given by Guinea et al. (2000): E0 KI = 4

r

2π (4v A − v B ) `

(3.20)

where E 0 is the effective elastic modulus defined as equal to E for plane stress and E /(1 − ν2 ) for plane strain. v A is the vertical displacement of the quarter mid-side node and v B the vertical displacement of the outer vertex node (See Figure 3.5) .

3.3.1

Formulae for the stress intensity factor

With some similar approach the following expressions may also be obtained to compute K I (Guinea et al, 2000): E0 KI = 2

r

2π vA `

(3.21)


33

Figure 3.5: Quarter-point singular elements and coordinates for near crack-tip field description. Source:

Guinea et al, 2000

E0 KI = 12

r

2π (8v A − v B ) `

(3.22)

These methods provide objective ways of computing K. Similar expressions can be obtained for K I I , using 3.11. Again in mixed mode situations the superposition principle may be applied

3.4

ANSYS commands for the computation of K

3.4.1

Crack opening displacement

ANSYS offers a built in method for the computation of the Stress intensity factor (K). Although this method is related to displacement extrapolation it is actually based on the concept of Crack Opening Displacement (COD) and uses the formula obtained by Paris and Sih which describes the crack opening near the crack tip for linear elastic-plastic materials: V KI 1 + κ p = p r 2G 2π

(3.23)

This expression can be easily obtained from 3.10 for θ = π. Since the crack opening displacement V can be obtained from the displacement solution at the nodes which define the crack face, the parameters A and B can be obtained by a simple linear fit. V p = A + Br r

(3.24)

V lim p = A r →0 r

(3.25)

Then, since:

K can be computed from 3.23: KI =

p 2G A 2π 1+κ

(3.26)

34

Fracture Mechanics

3.4.2

KCALC command

The main steps needed to perform the computation of K in a two-dimensional model are: 1. Define a path with three nodes. Where NODE1 must be the crack tip and NODE2 and NODE3 two nodes in the same crack face. If quadratic elements are used, a choice which gives good results is to use the three nodes of the crack tip element. 2. Define a cartesian local coordinate system with origin at the crack tip. 3. Execute the KCALC command

ANSYS Command: KCALC, KPLAN, MAT, KCSYM, KLOCPR

KPLAN Key to convert plane stress results into plane strain stress intensity factors: 0 - Plane strain and axisymmetric cases (default) 1 - Plane stress MAT Material number used in the extrapolation (defaults to 1). KCSYM Symmetry key: 0 or 1 - Half-crack model with symmetry boundary conditions in the crack-tip coordinate system. KII = KIII = 0. Three nodes are required on the path. 2 - Like 1 except with antisymmetric boundary conditions (KI = 0). 3 - Full-crack model (both faces). Five nodes are required on the path (one at the tip and two on each face). KLOCPR Local displacements print key: 0 - Do not print local crack-tip displacements. 1 - Print local displacements used in the extrapolation technique.

Example 3.2. Compute the Stress Intensity Factor for the plate of Figure 3.2.1 using ANSYS’ KCALC command. Solution to Example 3.2. The ANSYST M command sequence for this example is listed below. You can either type these commands on the command window, or you can type them on a file, then, on the command window enter /input, file, ext or just use copy and paste.


35

SOLVE FINISH CS,12,0,1,4,124 !Define local coordinate system at crack tip CSYS,12 ! Activate local coordinate system PATH,K1,3,10,50 ! Define 3-node path PPATH,1,1 PPATH,2,1013 PPATH,3,242 KCALC,0,1,0,1

!Execute KCALC

This file can be found at: ftp://amade.udg.edu/mme/FEmet/3_kcalc.dat

3.5

Proposed exercises

Problem 3.1. Numerical validation of Irwin’s hypothesis Irwin’s hypothesis may be used when plastic strains appear in the region near the crack tip. It is based on defining an equivalent case in the elastic regime, with an equivalent crack length. Let us keep working with the model 1_stress_sing_KSCON.dat we used in Chapter 1.

1. Introduce a perfect plasticity model as the material model. You can do this by adding the following lines after the material properties definition: TB,BKIN,1,1 TBDATA,1,270,0 where the value 270 MPa is the yield stress and 0 the hardening modulus. 2. Now increase the applied stress to a value that ensures that plastic strains appear near the crack tip (representative results are obtained for about 40 MPa). 3. Obtain a curve of S y for r between 0 and 0.3 mm, approximately. 4. Compute the equivalent crack length according to the Irwin’s hypothesis. In the former plot you can obtain the crack length for S y =270 MPa (the yield strength). Compare both values of the equivalent crack length. 5. Using the analytical expression of the stress field in the nearby of a singularity, plot the S y curve for the equivalent crack length of the former point. Compare this curve with the one of the question 3 6. Observe the results and comment on about the validity of Irwin’s hypothesis.

36

Fracture Mechanics

Problem 3.2. Consider again the side-cracked plate of Figure 3.2.1. Compute the mode I stress intensity factor using equations 3.21, 3.22, 3.20. Compare the results with those obtained with the other methods seen in this chapter. Comment on the results. Problem 3.3. Superposition principle Proof the principle of superposition can be used as schematized in Figure 3.6.

Figure 3.6: Proposed exercise # 2.

Consider a cracked plate submitted to an stress σ (A). Consider the same plate with the same stress σ but also closing stresses which make the crack remain closed (B). Consider the plate submitted only to the closing stresses but in the oposite direction (C).

1. Considering that the superposition principle is applicable for a single opening mode, discuss how could you compute K I(C ) . 2. Proof that K I(A) = K I(B ) + K I(C )

3.6

Further reading

Guinea G.V, Planas J. and Elices M. (2000) K I evaluation by the displacement extrapolation technique. Engineering Fracture Mechanics 66:243-255. Tada H., Paris P.C., and Irwin G.R. (2000) The Stress Analysis of Cracks Handbook. ASME Press. 3rd Edition

Chapter 4

Computational Fracture Mechanics III: Computation of the J-integral 4.1

Introduction

To complete the review of computational analysis for Linear Elastic Fracture Mechanics we will summarize the concept of J-integral and we will use ANSYS to compute it. Let us consider a line integral going around the crack tip and starting in one side of the crack and ending at the other side of the crack, as shown in figure 4.1.

Figure 4.1: J integral.

It can be shown that the following integral is independent of the path for any curve which satisfies the former conditions: ¶ Z µ ∂~ u J= U dy −~t · · ds ∂x Γ

(4.1)

were U is the strain energy density (U = 21 σ : ε), ~t is the traction vector defined by the external normal ~ n, ~ u is the displacement field and ds is an infinitesimal in the direction of the curve.

37

38

Fracture Mechanics

The integral is actually an equilibrium, for any path not including the crack, that is starting and ending at the same point, J=0, so if the curve starts at one side of the crack and ends at the other side, its value equals the energy inverted on the crack. The J-integral is also useful in non-linear fracture mechanics but, since in LEFM its value equals the energy release rate G ,

4.2

The J integral with ANSYS

As usual, ANSYS help describes properly the procedure to compute the J-integral. Here we summarize this procedure for bidimensional cracks: 1. Start the new computation of the J-integral with: CINT,NEW,ID where ID is an integer identifying the path, for instance 1. 2. Define the node at the crack tip and the crack plane normal with: CINT,CTNC,CMNAME where CMNAME is the name of a node component1 CINT,NORMAL,par1,par2 where par1 is a coordinate system identifier and par2 is an axis of the coordinate system 3. Specify the number of contours n to compute with the command: CINT,NCONTOUR,n 4. Activate the option for symmetry conditions, if present: CINT,SYMM,ON 5. Specify the output controls: OUTRES,ALL or OUTRES, CINT 6. Finally, the results for the value of the J integral may be listed or plotted: PRCINT,ID PLCINT,PATH,ID where ID is the crack identifier. Example 4.1. Compute the J-integral for the cracked plate of Figure 3.2.1, by using the ANSYS built-in method. 1 The command CM,CMNAME,NODE stores the selected nodes under a node component of name CMNAME.

Chapter 4. Computational Fracture Mechanics III: Computation of the J-integral

39

Solution to Example 4.1. The following commands may be used in any of the parametrized models we used before, with the crack tip located at (a,0) to define the J-integral computation. The ANSYST M command sequence for this example is listed below. You can either type these commands on the command window, or you can type them on a file, then, on the command window enter /input, file, ext or just use copy and paste. FINISH \PREP7 CSYS,1 NSEL,S,LOC,X,a !Select the crack tip node, located at (a,0) NSEL,R,LOC,Y,0 CM,CRACK,NODE NSEL,ALL CINT,NEW,1 CINT,CTNC,CRACK CINT,NORMAL,0,2 CINT,NCONTOUR,20 CINT,SYMM,ON OUTRES,CINT This file can be found at: ftp://amade.udg.edu/mme/FEmet/4_j_int.dat

After solving the model and in the \POST1 module, results ofr the J integral may be obtained with the commands PRCINT,1 or PLCINT,PATH,1. It is important to set a sufficient number of contours in the command CINT,NCONTOUR,n, so the integral converges to a value.

4.3

Proposed exercises

Problem 4.1. Compute the J-integral for the model 1_stress_sing_KSCON.dat we used in Chapter 1. Compare the value of J, with that of G and K, obtained in the corresponding examples.

4.4

Further reading

• Rigby R.H. and Aliabadi M.H. Decomposition of the mixed-mode J-integral - Revisited. International Journal of Solids and Structures 35(1):2073-2099, 1998.

Chapter 5

Computational Fracture Mechanics IV: Cohesive zone modeling 5.1

Introduction

Whilst Linear Elastic Fracture Mechanics assumes the presence of a crack in a perfectly elastic brittle or quasi-brittle material this is an idealization. Generally, in the nearby of the crack tip there exists a zone where the material is damaged due the presence of microcracks. When the number of microcracks grow a larger crack is formed and crack growth takes place. This region of the material is called Failure Process Zone or, in the case of crack growth modeling, cohesive zone. Some modeling techniques treat the material in this more realistic manner: before crack growth the region at the crack tip follows a failure process. This chapter summarizes the different possibilities included in ANSYS for the cohesive zone modeling.

5.2

Cohesive laws

Cohesive laws describe mathematically the separation or debonding of two material surfaces. They are usually presented as τ - δ curves. τ is the stress acting to separate the surfaces and δ the relative displacement between them. The different cohesive laws have some similarities: • Some positive slope region in which, when an increase in τ implies an increase in δ. • Some inflexion point τm . Once this point is reached, the cohesive material starts the failure/damage process. • Some negative slope region. Since the material is damaged the stress to achieve larger δ decreases. • Some δm for which τ = 0, which means the total damage of the material 41

42

Fracture Mechanics

The behaviour of the cohesive material is sketched in Figure 5.1. If the material is loaded with τ < τm , the unload follows the same path since the material is not damaged. On the other hand, if the material is loaded producing some τ > τm , the material starts to damage and then the unload follows the secant.

Figure 5.1: Cohesive law

This material behaviour can be modeled with different laws. Usually the linear (sometimes called bilinear), linear-parabolic, exponential and trapezoidal are included in the commercial FE software. They are sketched in Figure 5.2. ANSYS includes only the bilinear and exponential laws.

Figure 5.2: Usual cohesive laws

5.2.1

Bilinear law

This is the cohesive law used for the contact elements so the names of the variables are slightly different • P normal contact stress (tension). Equals τ • K n : normal contact stiffness • u n : contact gap. Equals δ

43

Chapter 5. Computational Fracture Mechanics IV: Cohesive zone modeling

Figure 5.3: Bilinear Cohesive law

• u¯ n : contact gap at the maximum normal contact stress (tension) • u nc : contact gap at the completion of debonding • d n : debonding (damage) parameter. d n = 0 for the virgin material and d n = 1 for the totally damaged material For mode II or mixed mode, additional parameters are required.

5.2.2

Exponential law

This law is the only one available for interface elements. τ = exp σmax ∆n exp−∆n exp−∆t

2

(5.1)

with: • σmax : stress for which crack opening starts • δn : maximum normal displacement • δt : maximum tangential displacement Parameters must be given so: Z τ(δ) dδ = Gc

5.3

(5.2)

Cohesive elements in ANSYS

ANSYS offers two different possibilities for the cohesive zone modelling. A straight forward manner is the use of interface elements. A second approach, is the use of ANSYS’ contact elements together with a cohesive law.

44

Fracture Mechanics

5.3.1

Cohesive zone modeling with interface elements

Element type Cohesive elements are referred as interface elements in the literature because of their topology. That is, the element is located in the interface between two solid structural elements to simulate the debonding process between them.The different interface elements available in ANSYS are shown in the next Table: Element Characteristics Interface Element Structural Elements 2D, linear

INTER202

PLANE42, VISCO106, PLANE182

2D, quadratic

INTER203

PLANE2,

PLANE82,

VISCO88,

SOLID95,

SOLID186,

SOLID46,

SOLID64,

PLANE183 3D, quadratic

INTER204

SOLID92, SOLID187

3D, linear

INTER205

SOLID45,

SOLID65, SOLID185, SOLIDSH190 Material definition As mentioned before, when using interface elements, the only material model which can be used is the exponential. It needs of three parameters: • σmax : maximum stress • δ¯n : normal displacement at maximum stress • δ¯t : tangential displacement at maximum stress

ANSYS Command: Material definition for interface elements TB,CZM,MAT,NTEMP,NPTS,EXPO TBDATA,1,SMAX,DN,DT where SMAX is σmax , DN is δ¯n and DT is δ¯t .

Example 5.1. DCB test modeling with interface elements Since the DCB specimen is controlled to always be in crack-opening situation, interface elements may be successfully employed in the modeling of this test. Solution to Example 5.1. The following file reproduces an Example of ANSYS’ verification manual which aim is to test its cohesive modeling with a DCB test.


45

The ANSYST M command sequence for this example is listed below. You can either type these commands on the command window, or you can type them on a file, then, on the command window enter /input, file, ext or just use copy and paste.

FINISH /CLEAR /COM,ANSYS MEDIA REL. 11.0 (10/27/2006) REF. VERIF. MANUAL: REL. 11.0 /TITLE, VM248, DELAMINATION OF DOUBLE CANTILEVER BEAM - 2D PLANE STRAIN /COM,

REF: ALFANO, G. AND CRISFIELD, M. A.,

/COM,

"FINITE ELEMENT INTERFACE MODELS FOR THE DELAMINATION ANALYSIS

/COM,

OF LAMINATED COMPOSITES: MECHANICAL AND COMPUTATIONAL ISSUES"

/COM,

INT. J. NUMER. METH. ENGNG 2001, 50:1701-1736.

/PREP7 ET,1,182

!* 2D 4-NODE STRUCTURAL SOLID ELEMENT

KEYOPT,1,1,2

!* ENHANCE STRAIN FORMULATION

KEYOPT,1,3,2

!* PLANE STRAIN

ET,2,182 KEYOPT,2,1,2 KEYOPT,2,3,2 ET,3,202

!* 2D 4-NODE COHESIVE ZONE ELEMENT

KEYOPT,3,3,2

!* PLANE STRAIN

MP,EX,4,1.353E5

!* E11 = 135.3 GPA

MP,EY,4,9.0E3

!* E22 =

9.0 GPA

MP,EZ,4,9.0E3

!* E33 =

9.0 GPA

MP,GXY,4,5.2E3

!* G12 =

5.2 GPA

!MP,GYZ,4,5.2E3 !MP,GXZ,4,3.08E3 MP,PRXY,4,0.24 MP,PRXZ,4,0.24 MP,PRYZ,4,0.46 GMAX = 0.004 TNMAX = 25

!* TENSILE STRENGTH

TB,CZM,5,,,EXPO

!* COHESIVE ZONE MATERIAL

TBDATA,1,TNMAX,GMAX,1000.0 RECTNG,0,100,0,1.5

!* DEFINE AREAS

RECTNG,0,100,0,-1.5 LSEL,S,LINE,,2,8,2

!* DEFINE LINE DIVISION

LESIZE,ALL,0.75 LSEL,INVE LESIZE,ALL, , ,200 ALLSEL,ALL TYPE,1

!* MESH AREA 2

MAT,4 LOCAL,11,0,0,0,0 ESYS,11 AMESH,2 CSYS,0 TYPE,2 ESYS,11 AMESH,1 CSYS,0 NSEL,S,LOC,X,30,100 NUMMRG,NODES ESLN TYPE,3 MAT,5

!* MESH AREA 1

46

Fracture Mechanics

CZMESH,,,1,Y,0,

!* GENERATE INTERFACE ELEMENTS

ALLSEL,ALL NSEL,S,LOC,X,100

!* APPLY CONSTRAINTS

D,ALL,ALL NSEL,ALL FINISH /SOLU ESEL,S,TYPE,,2 NSLE,S NSEL,R,LOC,X NSEL,R,LOC,Y,1.5

!* APPLY DISPLACEMENT LOADING ON TOP

D,ALL,UY,10 NSEL,ALL ESEL,ALL ESEL,S,TYPE,,1 NSLE,S NSEL,R,LOC,X NSEL,R,LOC,Y,-1.5

!* APPLY DISPLACEMENT LOADING ON BOTTOM

D,ALL,UY,-10 NSEL,ALL ESEL,ALL NLGEOM,ON AUTOTS,ON TIME,1 NSUBST,40,40,40 OUTRES,ALL,ALL SOLVE

!* PERFORM SOLUTION

FINISH /POST26 NSEL,S,LOC,Y,1.5 NSEL,R,LOC,X,0 *GET,NTOP,NODE,0,NUM,MAX NSEL,ALL NSOL,2,NTOP,U,Y,UY RFORCE,3,NTOP,F,Y,FY PROD,4,3, , ,RF, , ,20 /TITLE,VM248, DCB: REACTION AT TOP NODE VERSES PRESCRIBED DISPLACEMENT /AXLAB,X,DISP U (mm) /AXLAB,Y,REACTION FORCE R (N) /YRANGE,0,60 XVAR,2 PLVAR,4 PRVAR,UY,RF *GET,TMAX,VARI,4,EXTREM,TMAX

!* TIME CORRESPONDING TO MAX RFORCE

FINISH /POST1 SET, , , , ,TMAX

!* RETRIEVE RESULTS AT TMAX

NSEL,S,NODE, ,NTOP

!* SELECT NODE NTOP

*GET,RF_NTOP,NODE,NTOP,RF,FY

!* FY RFORCE AT NODE NTOP

*GET,UY_NTOP,NODE,NTOP,U,Y

!* DISP AT NODE NTOP CORRESPONDING TO RFORCE

RF_MAX = RF_NTOP*20

!* PLANE STRAIN OPTION AND WIDTH = 20 mm

SET,LAST

!* RETRIEVE RESULTS AT LAST SUBSTEP

*GET,RF_END,NODE,NTOP,RF,FY

!* FY RFORCE AT NODE NTOP AT LAST SUBSTEP

*GET,UY_END,NODE,NTOP,U,Y

!* DISP AT NODE NTOP CORRESPONDING TO RFORCE

RF_END = RF_END*20

!* PLANE STRAIN OPTION AND WIDTH = 20 mm

*DIM,LABEL,CHAR,2,2 *DIM,VALUE,,2,3


47

*DIM,VALUE2,,2,3 LABEL(1,1) = ’RFORCE’,’DISP ’ LABEL(1,2) = ’FY (N)’,’UY (mm)’ *VFILL,VALUE(1,1),DATA,60.0,1.0 *VFILL,VALUE(1,2),DATA,RF_MAX,UY_NTOP *VFILL,VALUE(1,3),DATA,ABS(RF_MAX/60.0),ABS(UY_NTOP/1.0) *VFILL,VALUE2(1,1),DATA,24,10.0 *VFILL,VALUE2(1,2),DATA,RF_END,UY_END *VFILL,VALUE2(1,3),DATA,ABS(RF_END/24.0),ABS(UY_END/10.0) /COM /OUT,vm248,vrt /COM,------------------- VM248 RESULTS COMPARISON -------------/COM, /COM,

|

TARGET

|

ANSYS

|

RATIO

/COM, /COM,MAX RFORCE AND CORRESPONDING DISP USING INTER202: /COM, *VWRITE,LABEL(1,1),LABEL(1,2),VALUE(1,1),VALUE(1,2),VALUE(1,3) (1X,A8,A8,’

’,F10.3,’

’,1F10.3,’

’,1F5.3)

/COM, /COM,RFORCE CORRESPONDING TO DISP U = 10.0 USING INTER202: /COM, *VWRITE,LABEL(1,1),LABEL(1,2),VALUE2(1,1),VALUE2(1,2),VALUE2(1,3) (1X,A8,A8,’

’,F10.3,’

’,1F10.3,’

’,1F5.3)

/COM,----------------------------------------------------------/OUT FINISH *LIST,vm248,vrt This file can be found at: ftp://amade.udg.edu/mme/FEmet/5_VM248.dat

Figure 5.4 shows the force-displacement curve for the DCB specimen, as obtained from the model using cohesive elements.

Figure 5.4: Force-displacement curve

48

Fracture Mechanics

Figure 5.5: Force-displacement curves for different parameters of the cohesive law

5.3.2

Cohesive zone modeling with contact elements

Element type On the other hand, for complex boundary conditions which may not always tend to open the crack, contact elements may also be employed. In this case Element Formulation Usage

Target element

CONTA171

linear

2-D 2-Node Surface-to-Surface

TARGE169

CONTA172

quadratic


TARGE169

CONTA173

linear


TARGE170

CONTA174

quadratic


TARGE170

Material definition When using contact elements the only material model which can be employed is the bilinear. This can be defined in ANSYS in two different ways: by maximum traction and maximum separation (CBDD) or by maximum traction and critical energy release rate (CBDE). ANSYS Command: Material definition for cohesive zone modeling through contact elements TB,CZM,MAT,NTEMP,NPTS,CBDX (changing X by D or E) TBDATA,1,C1,C2,C3,C4

Example: DCB test modeling with contact elements The DCB test may also be modelled with contact elements. This option requires higher computational time.


49

Solution to Example 5.2. The ANSYST M command sequence for this example is listed below. You can either type these commands on the command window, or you can type them on a file, then, on the command window enter /input, file, ext or just use copy and paste.

finish /clear /prep7 et,1,182

! solid 4 node element

keyopt,1,3,2

! plane strain

et,2,182 keyopt,2,3,2 et,3,169

! target 2d element

et,4,171

! 2d contact element

keyopt,4,12,5

! bonded: cohesive law must be defined

et,5, 169

! target 2d element

et, 6, 171

! 2d contact element

keyopt,6,4,2

! Nodal point contact

keyopt,6,2,4

! Lagrange multiplier method

MP,EX,1,1.353E5

!* E11 = 135.3 GPa

MP,EY,1,9.0E3

!* E22 =

9.0 GPa

MP,EZ,1,9.0E3

!* E33 =

9.0 GPa

MP,GXY,1,5.2E3

!* G12 =

5.2 GPa

MP,GYZ,1,5.2E3 MP,GXZ,1,3.08E3 MP,PRXY,1,0.24 MP,PRXZ,1,0.24 MP,PRYZ,1,0.46 kopen

= 1.e6

smax=25

! Stiffness contact ! Definition of Cohesive law

gic=0.26 tb,czm,2,1,1,cbde tbdata,1,smax,gic,smax,gic,1.e-5,1e15 !Geometry a=35 length=100 h=1.5 l=length/2 rectng,0,length,0,h rectng,0,length,0,-h e_size=0.5 esize,e_size type,1 mat,1 local,11,0,0,0,0 esys,11 amesh,2 csys,0 type,2 esys,11

! Element size

50

Fracture Mechanics

amesh,1 ! Contact between specimen two arms real,3

! real set of contact between arms

r,3 rmodif,3,3,-kopen

! Normal contact stiffness

rmodif,3,12,-kopen

! Tangential contact stiffness

asel,s,area,,1 nsla,s,1 nsel,r,loc,y,0 type,3 mat,2 esurf asel,s,area,,2 nsla,s,1 nsel,r,loc,y,0 nsel,r,loc, x, 0, a+e_size/2 real,3 mat,3 type,6 esurf allsel,all nsel,all finish /solu dk,6,all dk,3,uy,30 nsel,all esel,all eqslv,front neqit,200 nropt, unsymm nlgeom,on autots,on time,1 deltime,0.0005,0.000005,0.1 outres,all,all solve finish

This file can be found at: ftp://amade.udg.edu/mme/FEmet/5_DCB_comp.dat


5.4

51

Some remarks on element size

The use of cohesive elements, in any of the available forms, implies that the cohesive zone (failure process zone) must be meshed with a sufficient number of elements. A rule of the thumb is to use at least three cohesive elements for the failure process zone. Some estimation for the length of the FPZ should be used to determine a critical element size. For instance, the length of the fracture process zone for delamination in a unidirectional test specimen loaded in mode I can be estimated as: `fpz =

9π E 3G I c 32 (τo3 )2

(5.3)

where E 3 and τo3 are respectively the Young modulus and the strength for the direction 3 of the composite and G I c is the critical energy release rate for mode I. Under mixed-mode loading, the length of the failure process zone is larger than the obtained using the former expression, so the given estimation is conservative. For typical CFRP the latter expression leads to element size between 0.1 and 0.5 mm. Obviously this is unsuitable for large structures. Then, some engineering methods may be applied which allow the use of larger element sizes (Turon et al, 2007).

5.5

Proposed exercises

Perform a mesh-size dependence analysis for the simulation of the DCB test, using interface elements. You may use the Example given in section 5.3.1. Compare the results with equation 5.3.

5.6

Further reading

Mi, Y., Crisfield, M. A., Davies, G. A. O. and Hellweg, H. (1998) Progressive delamination using interface elements. Journal of Composite Materials 32(14):1246 − 1272. Alfano, G. and Crisfield, M. A. (2001) Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. International Journal for Numerical Methods in Engineering 50(7):1701 − 1736. Turon A., Dávila C.G., Camanho P.P., Costa J. (2007) An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Engineering Fracture Mechanics 74(10):1665−1682.

5.7

Aknowledgements

The autor is very grateful to Iñigo Llanos at IKERLAN for decisive help with the modelling of the cohesive zone using contact elements.

Ansys Fracture Tutorial

Recommend Documents