The possibilities to model crack propagation for glass in MARC

The possibilities to model crack propagation for glass in

MARC

Citation for published version (APA):

Goijaerts, A. M. (1994). The possibilities to model crack propagation for glass in MARC. (DCT rapporten; Vol. 1994.010). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/1994 Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

providing details and we will investigate your claim.

WFW-report 94.010

The possibilities to model crack propagation for glass in MARC

Ad Goijaerts

WF W-94.

o

1

o

- A - 11 - t. - h

__

nr - -

Title

Ac1 Goijaerts

The possibilities to model crack propagation for glass in MARC

Abstract

As part of the on-going research to model crack propagation for glass, the possibilities to implement this in the commercial Finite Element package MARC were investigated. The basic theories of fracture mechanics and damage mechanics, the primary tools for crack propagation, are given. Then, the possibilities to apply fracture mechanics in MARC are extensively described. Also the possibilities to apply damage mechanics in MARC for crack propagation in glass are discussed.

The theories about the Hertz contact and a possible test case, the Hertzian cone crack, about which a small literature review is provided, are mentioned in the appendix.

The main conclusion is that, at the moment, MARC is not a very appropriate Finite Element package for modelling crack propagation in glass.

WFW-94.018

WFw.94

.

o

10

contents

1 Introduction B 2 Fracture Mechanics 2 2.1 Introduction

. . .

2 2.2 Theory

. . .

2 2.2.1 Stress-intensity factors

. . .

2 2.2.2 J -integral

. . .

3 2.2.3 Crack propagation

. . .

4

2.3 Fracture mechanics in MARC

. . .

6

2.3.2 Mesh adaptation

. . .

8

+ 2.3.1 Applied method of Parks to obtain the

f -integral

. . .

6

3 Continuum Damage Mechanics 10

. . .

3.1 Introduction 10 3.2 Theory

. . .

BO 3.2.1 Basic concepts

. . .

18

3.2.2 One-dimensional continuum damage mechanics

. . .

11

3.2.3 Pathological mesh dependency

. . .

13

3.3 Damage mechanics in MARC

. . .

14

3.3.1 GURSON-model for ductile metals

. . .

14

3.3.2 OGDEN-model for rubber materials

. . .

14

3.3.3 Model for progressive composite failure

. . .

15

3.3.4 The low-tension material model

. . .

15

3.3.5 The test case and the results

. . .

17

4 Conclusions 21 References 23 A Hertzian cone crack 25 A . l The Hertz contact

. . .

25

A.2 The cone crack’s influencing conditions

. . .

27

WFW-94.018

1

Introduction

At the Philips Research Laboratories (NatLab) the fracture behaviour of glass is investi- gated. The aim is to develop a finite element model for the prediction of crack propagation. There are roughly two ways to describe crack propagation: fracture mechanics and dam- age mechanics. In principle fracture mechanics is a more accurate way but requires more computing time in a finite element model, due to continuous mesh adaptation during crack growth. Damage mechanics requires less computing time but is less accurate at present. The purpose of this study is to examine the possibilities of the commercial finite element package MARC to model crack propagation and, when possible, to compare the results with the developed model mentioned above, with respect to accuracy and computing time. In chapter 2 fracture mechanics will be discussed and in chapter 3 continuum damage mechanics will be discussed. Both chapters start with an overview of the theory followed by the possibilities for implementation in MARC. In chapter 4 some conclusions will be made.

In the appendix, a small literature review is given on a possible test case with well-known results: the Hertzian cone crack. If the Hertzian cone crack should be modelled, it would be very interesting to look whether it would be possible to influence the cone angle by different boundary conditions. If that would be the case, then a promising production method might be found.

This WFW-report is the result of a training period under supervision of J. Horsten and Prof. F. Baaijens. It is for me the second training period during the study Biomechanical Engineering a t the Eindhoven University of Technology.

WFW-94.010

2 Fracture Mechanics

2.1

Introduction

First in this chapter the basic theory of fracture mechanics will be discussed. T h e definition and the interpretation of the stress-intensity factors and the 1-integral will be given, as well as some aspects of crack propagation. Then the possibilities to apply fracture machanics in MARC will be described. The applied method will be discussed and the problem of mesh adaptation will be explained.

2.2 Theory

2.2.1 Stress-intensity factors

Fracture mechanics focusses on the analysis of a geometry with a crack. In figure 1 is shown that the stress at the crack-tip Is infinite.

+

t

Figure I: The stress singularity

The result of the infinite stress is a stress singularity. for the stress:

K o = -

&

o is infinite for T = 0.

1

l F

at the crack-tip

This can also be seen in the equation

(I>

The factor K is called the stress-intensity factor and varies with the different types of fracture. There are three different types of fracture, with three typical loading modes, which are shown in figure 2.

WFW- 94.

o

1

o

1

MODEI MODE li MODE 111

OPENING MODE SLIDING MODE TEARING MODE

Figure 2: The three independent fracture modes [12]

For every loading mode the stress-intensity factor E(is given below:

with

p

a function of the relevant geometry parameters, o is the normal stress, 711 and 7111 are the shear stresses for mode I1 and mode I11 in different, perpendicular directions as in figure 2 and a is the crack length. If these three stress-intensity factors are known, a crack propagation criterion can be made: Crack propagation will occur if

with Ki, the critical stress-intensity factors (or fracture toughnesses), which are experi- mental material parameters. These stress-intensity factors are related to the components of the :-integral, which will be discussed in the next paragraph.

2.2.2 J-integral

In the preceding paragraph the stress situation was considered to obtain a criterion for crack propagation. One can also look at the energy state with use of the :-integral. The 1-integral is defined as follows:

Where

I’

is an arbitrary, closed curve, with interior region R, see figure 3.

WFW-94.010

/

+

\

Figure 3: Definition of the $integral

The vectors

G,

z a n d U, are the unit normal vector, the stress vector and the displacement vector on the contour

r.

W is the specific elastic energy in that point of F. T h e :-integral is equal to

o'

for all

r

enclosing a region without singularities or discontinuities. The

y-

integral is not equal to

o'

if

r

contains a singularity or a discontinuity. A crack is such a discontinuity with a singularity at its tip, and it can be shown that the y-integral does not depend on the choise of

I?

[la].

The components of the y-integral are related to the stress-intensity factors:

with p. = G = - _{2 ( l + v )} E ( E is the Young's modulus and v is the Poisson's ratio) and with

cv =

Ilfil

can be interpreted as the amount of energy, available for cracking, the surface energy

needed for the surfaces to be cracked. The 1-integral can also be interpreted as the generalised force working in the point of discontinuity within

I

'

.

for plane stress and a = 3 - 4.v for plane strain.

2.2.3 Crack propagation

The criterion to determine whether a crack will propagate or not is often equation 5 from paragraph 2.2.1. To determine the direction of the crack propagation, the following possibilities are given in literature according to Wolbert [15] (see for explanation figure 4):

P.

Maximum Circumferential Stress : ( o ~ o ) ~ ~ ~

2. Minimum Principal Shear Stress : ( ~ $ ) ~ i ~ 3. Minimum Strain Energy-density : (S),i,

: KI1 = 0 4. Local Symmetry

5. Direction of $integral

: 1

WF

w-

94 ~

o

1

o

with: (1 - v2)(0L2.,

+

0iY) - 2.741

+

S = 2E U x x X

Figure 4: Coordinates and stresses near the crack-tip [15]

The principle of local symmetry

( K I T

= O) is also called "pure opening" because it as- sumes that an already existing crack will grow in the direction where it will be torn apart by only mode I (opening mode) loading. According t o Nurse [16] the criteria, mentioned above, have similar results for small values of K I I / K I (< 0.5).

Considering the last interpretation of the 1-integral in the preceding paragraph the fifth possibility t o determine the direction of crack propagation is provided: If a crack prop- agates it propagates in the direction of the 1-integral. This is accepted to be one of the best propagation direction criterions along with the determination of the K-factors. To find out whether a crack will propagate or not, the length of the xintegral vector should be compared with a critical value of the :-integral vector; The crack will propagate if:

in accordance with equation 5. If calculated:

with

u,

the Raleigh wave velocity,

the crack propagates, the propagation velocity can be Jc

'u = (1 -

-)u,

l

l

1

l

l

which is a material parameter.

WFW- 94

.o

1

o

2.3

Fracture

mechanics

in

MARC

2.3.1 Applied method of Parks to obtain the j-integral

Parks [5] presented a finite element technique, based on the energy release rate G. It will briefly be explained with reference to figure 5.

Figure 5: The method sEParks [5]

There is a rigid and a deformed region needed for this method. The rigid region is the area enclosed by

î o .

The deformed region is the region between

ro

and

ri.

In MARC a virtual displacement (AZ) of the rigid region and the crack-tip can be prescribed, when the problem is solved for the upper plot of figure 5. Then the strain energy difference A U , which exists between the two plots of figure 5 , will be calculated. This strain energy difference contains the internal elastic energy and the work, done by the external loads, if the crack would propagate over Al. The energy release rate G is the amount of energy dissipated per unit crack advance [la]. Therefore, dividing AU by Al we obtain [5]:

WFW-94. 0 1 O

In the case of linear elastic material behaviour with only mode I loading, G is equivalent to Jl [P2].

Jl = G

04)

MARC apparently assumes that there is only mode I loading because equation 14 is stated in the manual [18] as being true. There is only mode

I

loading if K I I = O. This is the case if

AZ

is chosen in the direction of crack propagation because the crack propagates in the direction for which

K I ~

= O [15] (No attention is payed at

Krrr

because the 2- dimensional case is considered). This means that this method cannot be used t o estimate the direction of crack propagation. Apparently MARC only uses this implementation for checking whether or not the crack will propagate in a given direction by comparing J1 to a certain J,,.

If a simple example with only mode I loading is considered, the accuracy of MARC, concerning the calculation of J1, can be checked. The example of figure 6 will be used, with a = 10, b = 20, h = 10 and cr = 5.0.104. The Young's modulus is chosen E = 7.0-1010 and Poisson's ratio v = 0.265, which are appropriate values for glass (all in SI-units). The used mesh is shown in figure 6. 8-node quadrilateral plane strain elements were used and at the crack-tip quarter point elements are used.

Quarter point elements are elements that can handle the stress singularity because two nodes are not situated at the middle of the sides but

i

from the crack-tip node and

2

from the other corner node. In this case they are also 8-node quadrilateral elements with three nodes at the crack-tip so that the element seems a triangular element. (See also figure 7.)

Figure 6: Example to evaluate the accuracy of the Jl calculation and the mesh

Also in figure 6 in the mesh, a typical stress pattern for a mesh with a crack-tip is drawn. According to Rooke and Cartwright [6] KI = 3.00

.

06

in the case of figure 6. With the use of equation 7 for plane strain, it follows that J1 = 9,39 theoretically. The cal- culation of MARC gives Jl = 9,60. This is a deviation of about 2%, which is not very large.

Although MARC calculates only

SI

and not the :--vector, there is a cumbersome way to estimate the direction of crack propagation. MARC calculates the strain energy dif- ference between the states before and after virtual movement of the rigid region. MARC can calculate AU for a certain

AZ

in many directions. The direction for which AU has its maximum is the direction of crack propagation. (A crack will propagate in the direction where the greatest amount of energy can be won.) One can imagine an iterative way to come to the maximum AU and that can be impiemented in a subroutine in MARC. Besides this all, Parks states in his articie [Sj that the rigid region can be shrunk t o the single node a t the crack-tip. In the implementation in MARC, however, this gives bad results. Also, Parks wrote that no special elements are needed for this implementation but MARC suggests to use quarter-point elements to achieve a higher accuracy.

2.3.2 Mesh adaptation

If the crack propagation direction is found with a certain accuracy in one increment, then the crack-tip must be propagated over a certain distance in the foliowing increment. Then, again, the crack propagation direction is estimated in the next increment and in the following increment the crack will propagate, and so on.

The crack propagation velocity can be calculated as shown in paragraph 2.2.3. Either one takes a certain time step for the following increment and calculates the crack-tip displacement or one takes a certain crack-tip displacement and calculates the time step for the following increment.

The crack-tip or a region around the crack-tip must be moved. Eventually, this gives severe problems concerning the mesh as can be seen in figure 7. Four figures are given: a is the starting mesh, b is the mesh after one move of the crack-tip region and c and d are the meshes after respectively three and five moving increments. After several increments some elements of the mesh will get inside-out (which will happen in the first moving step after figure 7d, the node pointed at by the arrow will come into another element), and this terminates the calculation of MARC. Not only the nodes in the crack-tip region should move but almost all the nodes of the entire mesh ought to be moved, only the nodes at the edge of the mesh should, of course, not be moved. Between the crack-tip region and the nodes at the edge, the displacement can be interpolated: nodes close to the crack-tip region are moved by almost the same displacement as the crack-tip, nodes close to the edge should hardly be displaced.

This method only gives satisfying results if the crack propagation is small compared with the measures of the mesh. If the crack propagation is large compared with the measures of the mesh, then elements and nodes should be added behind the crack-tip and should be removed in front of the crack. In MARC there is no option available that can do this or rezone the complete mesh if the edges are given. However, there is a possibility to write a subroutine in MARC to define the mesh changes and rezone the complete mesh, by defining every node and element change. But this is very difficult to do for a n arbitrary mesh. The writers of MARC are also working on that. Nevertheless, if they succeed, the computing time will be large.

Another possibility is t o calculate the crack-tip displacement and then go to another Finite Element Method program which can make a mesh if the edges are given. SEPRAN is such a program and for other purposes it has already been achieved to let SEPRAN make a

WFW-94.

o

1

o

mesh for each increment in MARC.

Figure 7: Mesh problem by moving crack-tip region

WFW-94.

o

P

o

3

Continuum Damage

3.1

Introduction

Firtly, the theory of the continuum damage mechanics will be presented, starting with the basic concepts and one-dimensional continuum damage mechanics. Secondly, the possibilities t o model continuum damage mechanics in MARC will be discussed. Some continuum damage models that can be implemented will be mentioned, as well as the model that has been used for a test case, and the results of that test case.

3.2

Theory

3.2.1 Basic concepts

In a realistic material there are discrete flaws which have a very small length scale compared to the object size. In continuum damage mechanics these discrete flaws are averaged or integrated over a finite volume. An internal damage variable is introduced and then it is possible to relate this damage variable to the change in mechanical properties, like Young’s modulus and stresses.

s

Figure 8: Volume with damage

For defining the damage variable look at figure 8. This volume is considered to be large enough compared to the size of the local defects in order to justify integration of the local defects over the element and small compared to the global dimensions. At the cross section, with normal

Z

and area S, there exist micro-cracks and cavities. The effective area

3,

which is the stress-transmitting surface, is smaller than or equal t o

S.

The area of the cavities is associated with the damage and is defined as:

WFW-94.810

The damage variable related to the plane with normal 5 is now defined as:

Because the effective area can not be larger than the original area and because a n area can not be negative: O

L

,!?

I: S.

This means that O

L

SD

I:

S and as a result 8

5 D,

5

1. If no damage has occured the effective surface is equal to the original surface and

D,

= O. If the siirface is entirely damaged than the effective area is e q d to 0 and

D,

= 1. It, s h d d be notified that such damage variable can be different in every load direction (it is related t o 6).

3.2.2 One-dimensional continuum damage mechanics

When one-dimensional damage mechanics are considered only one damage-variable D =

D,

is needed. This damage-variable can be used to get an effective stress and a n effective Young's modulus:

o o

o=--- -

o = x a = -

-

s

S

s

s - s ,

1 - %

_s

1 - D

With 6 the effective stress.

x a a = - f o r : 0 2 0 1 - D x a a=--- f o r ; o

<

O 1 - h D (17)

with h the closing-coefficient, O

5

h

I:

1, characteristic for the closing of micro-cracks and cavities (apparent damage reduction).

The principle of strain-equivalence [13] states that the effective stress working on the undamaged material will have the same deformation ( E =

$)

as the original stress working

on damaged material ( E =

z),

with the effective Young's modulus (see figure 9).

/

D=O

.

same

strain E F

\ D+O

Figure 9: The principle of strain-equivalence

WFW-94.010

Since =

&,

the principle of equivalence implies:

E

= (1 - D ) E ₍₂₀₎

Now the damage-variable will be drawn as a function of strain, where E,, is the critical

strain at which the first damage will appear and E, is the ultimate strain at which material

fracture will occur with an ultimate damage D , 1 (see figure 10).

I

Figure 10: Damage as a function of strain

The relation between damage and strain can be linear (1), it can be like (2), or it can be exponential (3 and 4). The influence on the stress-strain relation can be seen in fig- ure 11.

Figure 11: Stress-strain relation with damage

The curves 1 t o

4

correspond approximately with the curves in figure 10. The area under a curve, multiplied by a characteristic crack width (see paragraph 3.2.3), is the work done

WFW-94.010

by the stress in the fracture process, also called the fracture energy, G j .

Brekelmans [13] has mentioned phenomenological damage laws for three characteristic mechanical load cases: large plastic deformation, creep and fatigue. Damage laws for brittle materials, like glass, were not mentioned because they give some problems and those problems will be mentioned in paragraph 3.3. In that paragraph also the strategy to get to two- and three-dimensional damage models by means of the one-dimensional description given above will be discussed.

3.2.3 Pathological mesh dependency

In many cases with a continuum damage mechanics application in the Finite Element Method a pathological mesh dependency is observed. At present it is an important research topic to make a continuum damage model which is mesh independent. The cause of the mesh dependency can be found in the physical characteristic crack width. In a finite element model the crack width is related to the characteristic element length. I t has not yet been achieved to make a fully mesh independent continuum damage model.

WF W- 94.

o

1

o

3.3 Damage mechanics

in MARC

3.3.1 GURSON-model for ductile metals

In ductile materials, given the appropriate loading conditions, voids will form in the ma- terial. They will grow and coalesce, leading to crack formation and potential failure. Experimental studies have shown that these processes are strongly influenced by the pres- ence of hydrostatic stress. Gurson studied microscopic voids in materiais and derived a set of modified constitutive equations for eiastic-piastic materials considering growing yield surfaces [8]. Tvergaard and Needleman modified the model of Gurson with respect to the behaviour at small void volume fractions and at void coalescence by introducing two extra parameters in the constitutive equations of Gurson. The resulting model has been implemented in MARC (see [18] volume A page 6-106).

Glass is not a ductile, elastic-plastic material and failure of glass has nothing to do with the formation of voids where the model has been based upon. So this model can not be used to model crack propagation in glass.

3.3.2 OGDEN-model for rubber materials

This model is intended for filled rubber-like materials under repeated loading (see fig- ure 12).

4

E

Figure 12: Hysteresis loop in damaged rubber material

During this repeated loading softening can occur. Softening is the decrease of the stiffness by increasing stress and softening is often associated with damage. In MARC one standard scalar function for representing the current amount of damage is implemented, which is a modification of the one proposed by Sirno [li], but it is possible to define your own function. To capture the experimental observation that the stress softens under cyclic deformation the instantaneous strain energy is now modified by multiplying with this scalar function.

This option can also be used in MARC with a viscoelastic material, but it is not suitable to model crack propagation in glass.

WFW-94.

o

10

3.3.3 Model for progressive c o m p o s i t e failure

A

model has been put into MARC to allow the progressive failure of certain types of composite materials. The aspects of this model are defined below.

o Failure occurs when any one of the failure criteria is satisfied. Failure criteria are mostly a combination of critical stresses and strains and can be defined by the user.

o Upon failure, the material moduli for orthotropic materials at the integration points are changed such that all moduli have the lowest original modulus.

o Upon failure for isotropic materials, the failed moduli are taken as 10% of the original moduli.

o If there is only one modulus, such as in a beam or truss problem, the failed modulus is taken as 10% of the original one.

o There is no healing of the material.

This model could perhaps be used to model crack propagation in glass but it is rather inflexible. The model can not be defined or changed by the user, such as that 10% for instance. Another reason for choosing the next model is that in this model the Young’s modulus is changed in all directions (in the isotropic case), which is not the case in the next model.

3.3.4 The low-tension m a t e r i a l m o d e l

In MARC an option is implemented to handle low-tension material, like concrete. These materials are characterized by the stress-strain relation in figure 13.

I

_ _ _ _ -

A

I

/ I C I II

Figure 13: Uniaxial stress-strain relation for low-tension material

W F W- 9 4. 0 1 O

With o,, the critical stress, E,, the critical strain, at which the damage formation starts, E the Young's modulus, E, the softening rnoduIus, Eu the virtual modulus a t a certain strain E ~ , with the ultimate cracking strain, at which the

material is fully broken. When the strain is increasing the stress will increase too, until the stress reaches oer. Then the damage variable

D

begins to increase (figure 10 and figure 11) and increases till D = 1 at E,.

This stress-strain relationship is uniaxial, but that does not mean that it can only be used for one-dimensional problems. Áiso two- and three-dimensional problems can be handled by this model. In that case the principal stresses must be calculated. i f the crack propa- gates, it will propagate in the direction perpendicular to the largest principal tensile stress. Then, in the direction of the largest principal stress the uniaxial stress-strain relation is applied and the new virtual stiffness, valid for an integration point in a strain increment in the direction of the largest principal stress, will be calculated. The stiffness perpendicular t o the largest principal stress will not be altered. That is the main difference with the model of paragraph 3.3.3.

This expansion of the uniaxial model to the two- and three-dimensional case is imple- mented in MARC and will be used for the modelling of crack propagation in glass. To implement that model the different values for a,,, E and E, must be defined.

8 N

E = 7.0. 10"

5

Now

E,

can be found because the fracture energy for glass is known to be approximately 8.0

-$-

= 8.0 and the fracture energy divided by the characteristic element length rep- resents the area under the stress-strain curve, w. So now the softening modulus can be calculated.

With this softening modulus, the stress-strain relation for glass should be like that in figure 14. It is noted that

between E,, and E,, and

oCT 3.0. 10

<

E , ~ , which is in contrast with figure 13.

Figure 14: Uniaxial stress-strain relation for glass

WFW-94.010

E, depends on the characteristic element size, but this element size should be smaller than

1.25

.

rn

to produce an E, larger than &CT and that is not very realistic. So in this

case the stress-strain relation is non-unique and if it is used in MARC it will be difficult to Iet this problem converge. For that reason we choose to model the stress-strain relation given in figure 15.

Figure 15: The modelled stress-strain relation

This is a relatively simple model: when the stress reaches the value of the critical stress, the material breaks at once, the stiffness is made zero in one increment. But for glass this simple model is closer to reality than the model of figure 13. With this model the test case in the next paragraph is performed.

3.3.5 The test case and t h e results

As discussed in paragraph 3.2.3, an important aspect of the used model that must be checked is the mesh dependency. If the method is mesh dependent, then for the same problem different meshes would result in different, solutions. This is of course undesirable. The used test case to see whether or not the implemented model is mesh dependent is S ~ Q W ~ in figure 16.

WFW-94.01

o

Figure 16: The used test case

This is a square 2-dimensional model. The square will be divided into 400 8-node quadri- lateral plane-strain elements. All elements will have the same stiffness ( E = 7.Q.1010

3)

and the same critical stress (ocr = 3.0

.

lo8

3)

except one element at the left side of the square (see figure 16). That element will have a stiffness which is nearly zero. This weak element will cause stress concentration when o is raised. Due to the stress concentration the element a t the right side of the weak element will, with increasing o, reach the critical stress first and will crack first. Following this way of thinking, the crack ought to propa- gate horizontally to the right and this corresponds to experimental facts. While the crack is propagating, the stress transmitting area gets smaller and then o should also decrease, so that the maximum stress in the mesh should only be a little larger than the critical

stress, so that not too many elements would crack in the same increment, because this should cause numerical problems. This is done with the AUTO INCREMENT option in MARC. This option controls that the maximum stress in the mesh is only a little larger than the critical stress.

I I I I I I I I I I I I I I I I I I I

i,

I

h

Figure 17: Two different meshes

WFW-94.

o

B

o

Now the problem of figure 16 is solved with 2 different meshes, shown in figure 17. The two different meshes give two different solutions as can be seen in figure 18

L

I "" . ':

Figure 18: Two different solutions after 20 increments

The plotted contour lines in these plots represent the stress in y-direction. If there is a stress concentration of the y-component of stress then the element where the concentration is located will crack first.

For obtaining a better view on the process of fracture during the 20 increments, plots are included in appendix

B

for the regular square mesh and the irregular oblique square mesh. The solution of the regular mesh is in agreement with experimental results but the solution of the irregular mesh is not the correct solution for this problem. The crack follows the line of elements instead of the horizontal direction. The irregularity of the stress contour lines also indicates that the solution is questionable.

However, the solution is not as bad as expected because the crack passes a line of elements in the horizontal direction twice. (For this mesh the weak element is chosen beneath the middle, so that it is situated on the line of elements that has the largest angle with the x-axis. This way, the chance that the solution crosses a line of elements is the greatest. The correct solution should of course still go horizontally.)

When the solution crosses such a line there are severe convergence problems. The reason for that is not perfectly clear but one can regard that point as a bifurcation point: there is a choice between two almost 'good' solutions, and MARC has problems to pick the best. Although the solution is not as bad as expected, an employee of the MARC Analysis Research Company [18] claims that the model of figure 13 should be mesh independent. This means that there should be a difference in mesh dependency between the model of figure 13 and the model of figure 15. One could think that there should be a difference in Convergence because in one model the stiffness Is decreased radically and in the other it

is decreased gradually, step by step, but the model with the radically decreasing stiffness also converges. Still it is not clear why there should be a difference in mesh dependency between the two models. The greatest concentration of the second component of stress will still, most of the time, be at the element in the same line, at the right side of the cracked element. Future research is recommended to investigate the mesh dependency of the model of figure 13.

Additional numerical techniques, such as proposed e.g. in [2O], should be implemented in MARC. That could propably help in the battle against mesh dependency.

WF W-94. O

P

O

4

Conclusions

The following conclusions are made:

o In its present form the fracture mechanics procedures in MARC can not be used to model realistic crack propagation. The crack propagation criterion is not very sophisticated and an automatic rezoning facility is lacking. If there comes a more advanced rezoning option in MARC, the modelling of crack propagation will be possilole but the method will use very much computing time.

o In its present form, the damage mechanics in MARC can not model crack propaga- tion for glass properly, the method exhibits mesh dependency. The mesh dependency for other materials is not investigated but there are no indications that the method would be mesh independent for other materials. Additional numerical techniques should be implemented.

m Since both methods already fail in the simple case of a straight uniformly loaded crack, it can not be expected that a more complex case, such as a Hertzian cone crack (appendix A ) , can be treated properly.

WFW-94.010

WFW-94.010

I

References

F.C. Frank and

B.R.

Lawn

On the theory of Hertzian fracture.

Proceedings of the Royal Society, Vol. 299, pp. 291-306, 1967 K.L. Johnson, J.J. O'Connor and A.C. Woodward

The effect of the indenter elasticity on the Hertzian fracture of brittle materials. Proceedings of the Royal Society, Vol. 334, pp. 95-117, i973

I. Finnie and S. Vaidyanathan Hertzian ring cracks.

"Fracture Mechanics of Ceramics", Vol. 1, edited by R.C. Bradt, D.P.H. Hassel- man and F.F. Lange, Plenum Press New York, 1974

B.R.Lawn, T.R. Wilshaw and N.E.W. Hartley

A

computer simulation study of Hertzian cone crack growth. International Journal of Fracture, Vol. 10, pp. 1-16, 1974 D.M. Parks

A

stiffness derivative finite element technique for determination of crack tip stress intensity factors.

International Journal of Fracture, Vol. 10, pp.487-502, 1974 D.P. Rooke and D.J. Cartwright

Compendium of Stress intensity factors. Hillingdon Press, Uxbridge, UK, 1976

C.G. Knight, M.V. Swain and M.M. Chaudhri Impact of small steel spheres on glass surfaces.

Journal of Materials Science, Vol. 12, pp. 1573-1586, 1977 A.L. Gurson

Continuum theory of ductile rupture by void nucleation and growth. Journal of Engineering Materials and Technology, Vol. 99, pp. 2-15, 1977 Dov Bahat and Martin R. Sharpe

Dependence of Hertzian fracture angle on Poisson's ratio and indentation tech- niques.

Journal of Materials Science, Vol. 17, pp. 1167-1170, 1982

K.L.

Johnson Contact mechanics.

Cambridge University Press, 1985 J.C. Sim0

On a fully three-dimensional finite strain viscoelastic damage model: Formulation and computational aspects.

Computer Methods in Applied Mechanics and Engineering, Vol. 60, pp. 153-173, 1987

LI

WFW-94.010

P. Schreurs

Dictaat breukmechanica, nr. 4694.

Eindhoven University of Technology, Mechanical engeneering, 1990 M. Brekelmans

Dictaat Continuum Damage Mechanica, nr. 4693.

Eindhoven University of Technology, Mechanical engeneering, 1990 Y. Akimune

Hertzian cone crack in S i c caused by spherical particle impact. Journal of Materials Science Letters, Vol. 9, pp. 659-662, 1990 P.M.M. Wolbert

Fracture mechanics.

Technical Note 039/90NE, Philips Centre for Manufacturing Technology, 1990 A.D. Nurse and E.A. Patterson

A technique to predict the direction of edge cracks.

International Journal of Mechanical Science, pp. 253-264, 1990 Li Yingzhi and D.A. Hills

The Hertzian cone crack.

Transactions of the ASME, Journal of applied Mechanics, Vol. 58, pp. 120-127, 1991

MARC Analysis Research Corporation

User information manual, Volume A-E), Revision K.5. MARC Zoetermeer, 1992

P hone-service: f3 1- 79- 5 1841 i

K. Zeng, K. Breder and D.J. Rowcliffe

The Hertzian stress field and formation of cone cracks - I. Theoretical approach.

Acta Metallica Materia, Vol. 40, pp. 2595-2600, 1992 E.J. Sluys

Wave propagation, localisation and dispersion in softening solids.

PhD-thesis, Delft University of Technology, W.D. Meinema B.V. Delft, 1992

WFW-94.010

A Hertzian

cone

crack

In this appendix the findings of examining the literature on the Hertzian cone crack will be described. First, the Hertz theory of elastic contact will be discussed briefly. Then the cone crack and all the conditions that can influence it will be mentioned.

A.1

The Hertz

caritact

The Hertz theory describes the contact between two elastic solids with dissimilar ellipsoid surfaces. Nearly all contact forms between two elastic solids can be described with this theory because a plane surface can be regarded as an ellipsoid with infinite principal radius. On the other hand, looked at micro scale, there will always be spherical contact because a surface will never be an exactly straight plane. The restrictions and assumptions made in the Hertz theory are the following: The profiles are parabolic, the surfaces frictionless and the elastic half-space theory can be used. Also the displacements should be small compared to the characteristic radii. The basics of the Hertz theory will be discussed with reference to figure 19.

Figure 19: The Hertz theory of elastic contact [IO]

WFW-94.010

There are two surfaces: surface

P

has principal radius R1 and surface 2 has principal radius R2. (The surfaces are chosen convex for convenience; if a surface is concave the principal radius will be negative.) Now the deformation, due to a normal load

P,

is considered. During compression distant points in the two bodies 2-1 and Tz move towards O, parallel

to the z-axis, by displacements

61

and 62 respectively. If the solids would not deform, their profiles would overlap as shown by the dotted lines in figure 19. Due t o the contact pressure, the surface of each body is displaced parallel to O z by an amount U,, and Üz2 relative to the distant points TI and T2. Then there is a contact circle with radius a if axisymmetry is assumed and in the 2-dimensionai case the contact area is an infinite strip with width 2a. The pressure p ( ~ ) as a function of T

(1.1

5

a), a, and the relative

displacement

S

=

Si

+

62 can be calculated for the axisymmetric problem as well as for the 2-dimensional case [lo], introducing the following parameters:

-1 E * =

(T+-)

1 - v s 1 - v ;

E2 and:

The solution for the axisymmetric case is given by:

113

a2 9 P

a = - =

R 16RE*2

With the maximum pressure po:

While the solution for the two dimensional case (two cylinders with parallel axes) is:

a2 4~ & = - = -

R

TE* 2P 11' p ( r ) = -7j

(2

- 7-21 ra

The pressure, the contact area radius (or the contact area half width) and the normal displacements can be used to calculate the Hertz stress distribution, given in figure 20 in comparison with a stress distribution caused by a uniform pressure.

WFW-94.0 10

Uniform pressure p 1 _Hertz 1

-1.0 -0.5 -0.5 O O -0.5 -0.5 -1.0 -1.0 -1.5 -1.5

Figure 20: The stress distribution at the surface and along the axis of symmetry [lo]

For the axisymmetric case the stress distribution, at the surface (upper part of figure 20) and along the axis of symmetry (bottom of figure 20) caused by (left) uniform pressure and (right) Hertz pressure acting on a circular area radius a, is given. Here, g T , 0 0 and

gz are the principal stresses in cylindrical coordinates and 71 = ;loT - g e l is the principal shear stress.

A.2 The cone crack’s influencing conditions

Already in 1967 a stress field beneath a Hertzian contact area was described by means of the formulae given in section A.I. That stress field is given in figure 21.

WFW-94.010

0-005 0.005

Figure 21: An axisymmetric stress field under Hertzian contact

Frank and Lawn [i] also stated for the first time that the crack propagation should, t o a first approzimation, be perpendicular to the most tensile principal stress, and thus cor- respond, in a torsion-free stress field, to a surface delineated by the trajectories of the other two principal stresses. These trajectories are given in figure 2 1 by the dotted lines, starting at 0.8a, a and 1 . 2 ~ . This was in accordance with the known experimental facts. Later Lawn e.a. [4] and Johnson e.a. 621 found that a Hertzian crack starts at the surface at a position of respectively 1 . 0 5 ~ and 1 . 2 5 ~ . This range will be probably even wider. The crack radius can also be influenced by the size of the indenter, according to Johnson e.a. [2]. The cone angle, the angle between the crack and the vector normal to the surface, is a very interesting item that one would like to control. This is the case if for instance the Hertzian breaking of glass is seen as a producing method to make holes in a plate of glass. In this application one would like the cone angle to be very small. Actually the cone angle can vary from 35 deg. (Knight e.a. [7]) to 72 deg. (Yingzhi and Hills [17]) and there is probably even a wider range.

According to Akimune [14] and Johnson e.a. [2] the cone angle can become smaller when the Young’s modulus of the indenter becomes larger for both static and dynamic impact. Johnson e.a. [2] also states that the fracture radius depends on the ratio of the indenter Young’s modulus to the specimen Young’s modulus.

According to Akimune [14] and Knight e.a. [7] the cone angle can get smaller in a dynamic experiment by a higher impact velocity of the indenter.

Zeng e.a. [19], Bahat and Sharp [9] and Lawn e.a. [4] state that the cone angle gets smaller when the Poisson’s ratio gets smaller. According to Finnie and Vaidyanathan [3] it is even possible t o get a linear correlation between the cone angle and the Poisson’s ratio, like: No literature has been found about the possibilities to influence the cone angle by stress conditions; bending stress as well as tensile stress. One could imagine that the cone angle will be different with another value of the bending stress q, and the tensile stress ot in figure 22.

6 ( V ) = Q V Jr

p.

WFW-94.010

Figure 22: The cone angle can probably be influenced by bending and tensile stresses

In literature cone crack propagation has been simulated but a coupled simulation was not found. With a coupled simulation the following procedure is meant: In every time step the crack criteria are checked according t o the last configuration, considering the propagated crack. And then the crack is propagated over a small finite displacement. In literature often the stress field is calculated, and then, assuming that the propagating crack does not influence the stress field, the stress criteria are applied. This is done for instance by Zeng e.a. [19]. This is not a very accurate way of calculating crack propagation. The method, which is in development at the Philips Research Laboratories (NatLab), applies a coupled method to calculate crack propagation.

WFW-94.

o

1 0

B M A R C - P ~ O ~ S

Regular

square

mesh

>

1

W O al al al m al m al al m al W t- Q hu W m 7 v) O 8 $ t g g O

, 0 8 x

p- m

2

I ' w y F W O W O m O m u? W O O a W O W OD Q r

I=

W O W 0 0 o 7 v

5

2

2

e

2

r Q

..

M e T

.. ..

30 Philips Electronics N.V. 1994

WFW-94. OP O >

J

W

o

O al

t

o

t

al

t

al o 8 m W O al In N h W O UJ W O W c.> v> O O eo W O IC O UJ h O In IC eo o h $ 8 6 O

* o g z

5

lD F p: $! ul O 0 0

..

v> eo

.. ..

-

Philips Electronics N.V. 1994 31

WFW- 94. O 1 O

WFW-94.010

WFW-94.010 W O O W O al W O (0 O m W al al al al al al al u! W O W IC W O al In m W O U rn ₈ W 9 6 IC

$ 2

IC O O 0 0 0:

..

co 8 g o s 2 In l-

i &

b n! r

2

? T P

f

*

9 T u! "!

.. ..

.- Z Y > b l L 34 Philips Electronics N.V. 1994

WFW-94.010

Irregular oblique square mesh

WFW-94.010

1

WFW-94.010

WFW-94

.o

1

o

>

1i

WFW-94.018 >

A

Q E 1

_-

i I- C c z h e r

c

Philips Electronics N.V. 1994 39