Analysis of Three-Dimensional Cracks in Submodels

(1)

DEPARTMENT OF MANAGEMENT AND ENGINEERING

Analysis of Three-Dimensional

Cracks in Submodels

Master thesis in Solid Mechanics performed at

Linköping University

for

Saab Aerosystems

by

David Karlsson

LIU-IEI-TEK-A--07/0016--SE

January 2007

(2)

(3)

Datum Date 2007-01-19

Titel Analysis of Three-Dimensional Cracks in Submodels Title

Författare David Karlsson Author

Sammanfattning Abstract

A common technique to evaluate load paths in complex structures is to perform FE-calculations with relative large elements. This procedure gives no information regarding stress concentrations at e.g. holes or radius but this phenomenon can later on be investigated in details with local individual submodels. Displacements is taken from the global model and used to analyse stress concentrations and crack driving parameters in the submodel.

Today, the crack controlling stress intensity factors are in general cases obtained from handbook solutions of elementary cases. This method requires engineering judgements in a conservative manner and one way to improve the solution is to model the crack in its correct surroundings in a local three- dimensional submodel.

This master thesis is focused on the development of an automated support for analysing of three- dimensional cracks in submodels. The results from a global Nastran model can be imported to Trinitas and used for a more accurate stress and fatigue life analysis in a local model. Here a three-dimensional crack tip subdomain can be generated inside an eight point brick volume. The crack tip subdomain is specially designed and adjusted for accurate determination of stress intensity factors along the crack front. For example, all points are adjusted with respect to the brick volume and the crack size, triangular wedge elements are applied around the crack tip, the midpoints for these elements are moved to quarter points and the crack front is curved. The crack tip subdomain is validated against several reference cases and shows sufficiently good results with respect to the stress intensity factor.

Finally, the automated crack tip subdomain generation is applied to a geometrically complex part of a main wing carry-through bulkhead of a fighter aircraft in order to show the applicability of the procedure in an industrial environment.

Nyckelord Crack, fracture, stress intensity factor, submodel, three-dimensional Keyword Språk Language Svenska/Swedish x Engelska/English Rapporttyp Report category Licentiatavhandling x Examensarbete C-uppsats D-uppsats Övrig rapport

Serietitel och serienummer Title of series, numbering

ISRN: LIU-IEI-TEK-A--07/0016--SE Avdelning, institution

Division, Department Div of Solid Mechanics

Dept of Management and Engineering SE-581 83 LINKÖPING

(4)

(5)

Abstract

A common technique to evaluate load paths in complex structures is to perform FE-calculations with relative large elements. This procedure gives no information regarding stress concentrations at e.g. holes or radius but this phenomenon can later on be investigated in details with local individual submodels. Displacements is taken from the global model and used to analyse stress concentrations and crack driving parameters in the submodel.

Today, the crack controlling stress intensity factors are in general cases obtained from handbook solutions of elementary cases. This method requires engineering judgements in a conservative manner and one way to improve the solution is to model the crack in its correct surroundings in a local three-dimensional submodel.

This master thesis is focused on the development of an automated support for analysing three-dimensional cracks in submodels. The results from a global Nastran model can be imported to Trinitas and used for a more accurate stress and fatigue life analysis in a local model. Here a three-dimensional crack tip subdomain can be generated inside an eight point brick volume. The crack tip subdomain is specially designed and adjusted for accurate determination of stress intensity factors along the crack front. For example, all points are adjusted with respect to the brick volume and the crack size, triangular wedge elements are applied around the crack tip, the midpoints for these elements are moved to quarter points and the crack front is curved. The crack tip subdomain is validated against several reference cases and shows sufficiently good results with respect to the stress intensity factor.

Finally, the automated crack tip subdomain generation is applied to a geometrically complex part of a main wing carry-through bulkhead of a fighter aircraft in order to show the applicability of the procedure in an industrial environment.

(6)

(7)

Preface

This master thesis is the final project in the programme for a Master of Science degree in Mechanical Engineering. The work has been carried out at the Division of Solid Mechanics at Linköping University Institute of Technology for Saab Aerosystems in Linköping during the autumn 2006.

First of all I would like to thank my supervisor Hans Ansell at Saab Aerosystems for starting this project, sharing of contacts in the industry and location of applications on the Gripen aircraft.

Then I will express gratitude to my supervisor Bo Torstenfelt at the Division of Solid Mechanics at Linköping University for offering me this project, support with programming and always be close to assistance for discussion, guidance and advice. Finally I will thank Börje Andersson at FOI for reference cases, Lars Djärv at Saab for global results from Nastran and Anna-Lena Lagerström at Saab for drawings and coordinates of the wing attachment frame.

David Karlsson

(8)

(9)

1

Introduction

1.1 Saab

The history of Saab begins in Linköping in 1930 when the Swedish government decided to strengthen the country’s defence capabilities of its air force. The company was formed in 1937 and its primary aim was to meet the need for a domestic military aircraft industry in Sweden [1]. They became the dominant supplier to the Swedish Air Force and have created many generations of military aircraft, for example Tunnan, Lansen, Draken, Viggen and in the nineties Gripen.

Saab is today one of the leading high-technology companies in the world with approximately 12 000 employees. The main operations are defence, aviation, space and civil security. Saab develops, manufacture and delivers products and services for the defence market and for commercial markets. Saab has the world as its market, but research, development and production are carried out principally in Sweden.

Gripen is Saab’s most recent fighter and the world’s first fourth-generation fighter aircraft in active service. Saab Aerosystems has the overall system responsibility for the development of the Gripen and the aircraft is the business unit’s most important product. Other areas for Saab Aerosystems are unmanned aerial vehicles, tactical systems and advanced pilot training systems.

1.2 Background

A common technique to evaluate different load paths and global stresses in complex structure is to perform FE-calculations with relative large elements. In the global model the stiffness of each element are adjusted to be correct, but the exact geometry are forced to be simplified. This procedure gives no information regarding stress concentrations at e.g. holes or radius but this phenomenon can later on be investigated in details with local individual submodels. For a small number of elements in the global model coordinates, displacements and rotations are extracted and used as input data to the submodel. Here can concentrations and crack driving parameters be further investigated.

Today the crack controlling stress intensity factors for part through cracks are in general cases obtained from handbooks. The handbook solutions of elementary cases are derived from FE-analysis and often expressed in analytical formulas. The application of the elementary solutions in the design and sizing process involves and requires engineering judgments in a conservative manner.

(12)

One way to improve the quality of the solution and to reduce the amount of conservatism is to model the crack in its correct surroundings in a three-dimensional model. This procedure can easily be fit into the ordinary sizing process just by replacing the handbook solution with the solution from the local submodel.

1.3 Objective

The objective of this work is to pick a limited number of elements from a solved global FE-model and use the results to study three-dimensional cracks in local submodels. From the global model, solved by the FE-program Nastran, coordinates, displacements and rotations are taken and used as input values in the submodel, which is analysed in the FE-program Trinitas. In the detailed analysis, the geometry is supposed to be more accurate and an automated support shall be used for creating three-dimensional crack tip subdomains from single eight point brick volumes. The size and elliptical shape of the crack should be defined by the user. The function shall handle three types of common cracks, which are corner, surface and embedded cracks. The aim of the analysis is both to locate stress concentrations and determine the crack controlling stress intensity factor along the crack front. The stress intensity factor can be used in a fracture criteria in order to predict static failure and in a crack growth law in order to predict the stable growth of cracks under cyclic fatigue loading.

1.4 Procedure

The main work is concentrated on the automated routine that convert an eight point brick volume to a crack tip subdomain specially designed for calculation of stress intensity factors for three-dimensional cracks. First, a study of the fracture mechanics theory has to be done to understand how this crack tip subdomains shall be created to obtain accurate results in the FE-program Trinitas. After that, several reference cases need to be analysed to ensure that the mesh works for different sizes and shapes of cracks and for different load cases. In the application part, a wing carry-through bulkhead of the Gripen aircraft should be studied especially with respect to stress intensity factors for a three-dimensional corner crack. In this part results from a global model must be transferred and applied to the submodel.

1.5 Restrictions

The automated support for the crack tip subdomain generation will only work for standard eight point brick volumes. This means that it can be difficult to introduce a crack in all locations and different geometries in a model. The ideal shape of the transformed brick volume is a cube but the method will work for different sizes and shapes of the volume, for example circular sides. The results from the Trinitas program will be the crack controlling parameter called stress intensity factor. The stress intensity factor can be computed for three different modes of loading. This work will only consider the first mode which is the most common and dominated in industrial crack applications.

(13)

2

Theoretical background

2.1 Fracture mechanics

Fracture in machines has been a problem in engineering history ever since and may be more important today since constructions become more optimized and complex. Consider loaded structure. As long as the load is small enough, the component will only deform elastically and work as planned. If however, the component is sensitive to cracking due to e.g. inadequate design, defects from manufacturing, handling or bad quality materials, it may finally break and fracture will occur.

When fracture takes place in a component it will lose the function it was initially designed for. Due to the components function in the entire construction the fracture will cause different damages. The worst case is that the whole structure is depending on this detail and will collapse if the component breaks.

So knowledge about what happens with a loaded structure with crack initiations is important in order to avoid a failure or understand a fracture. Fracture mechanics is used to answer some main questions such as: will the crack grow, will it grow fast and unstable or slow and stable, if it grows stable at what rate does it grow, to what size can the crack grow and how many cycles or how long time does it take to become an unstable crack.

2.2 Stress analysis of cracks

Close to a hole, a radius or a groove, the stresses are higher than in undisturbed parts far away.

Figure 2.1 Irregular geometries which result in stress concentrations

(14)

The factor between the stress far away from the stress raiser, the nominal stress, and the maximum stress at the stress raiser is called stress concentration factor and is denoted Kt. The stress concentration factor for an isolated hole in a large plate subjected to

uni-axial tension, is three [2] and for a general case it can be written as

(2.1)

At the stress concentration the probability for a crack to initiate and propagation is higher than in the rest of the material. So in the design process it is important to avoid small radii because sharper notches give higher stresses.

As introduction to stress concentrations at the crack tip, a large plate with an elliptical hole subjected to uni-axial tension, Figure 2.2 is studied. The stress concentration factor is according to [2]

(2.2)

Figure 2.2 Large plate with an elliptic hole

Note that a circular hole a = b gives Kt = 3 as mentioned earlier. If the minor axis, b,

vanish to zero, the ellipse will be like a crack and the stress concentration factor goes to infinity. A crack is obtained when b is equal to zero. This implies that a crack is characterized by two surfaces meeting at a sharp edge at the crack tip. The crack front is the line connecting all adjacent sites where separation may occur subsequently. During continued separation this line will move along a surface called the fracture surface. The area of this surface will of course increase as the crack growths.

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = b a K_t 1 2 nom t Kσ σ_max =

(15)

To determine the stress at the material in front of the crack tip it is necessary to define three different modes that the crack can be loaded in [3], demonstrated in Figure 2.3.

Figure 2.3 Three modes of loading a crack

Mode I The crack is opened so the crack surface moves directly apart.

Mode II The crack is sheared in the plane where the crack surface slides over one another.

Mode III The crack is sheared in the plane where the crack surface moves relative to one another.

In general for a crack growing in a body the separation can be described as some combination of these three modes. In this work only Mode I will be studied because crack growth usually takes place in Mode I or close to it and Mode I is the most severe loading of the crack. Mode I is also encountered in the majority of actual engineering situations involving cracked components [3]. So consequently attention has been given both analytical and experimental results to the Mode I loading.

For the notation in Figure 2.4 the stresses in Mode I can be written as follows according to [3]

(2.3a)

(2.3b)

(2.3c)

KI is the stress intensity factor and will be further explained in the next chapter. It is

obvious that these local stresses could rise to extremely high values if r approaches zero, but this situation is prevented by the entering of plastic deformation at the crack tip [3].

Mode II Mode III Mode I z y x z y x z y x ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − = 2 3 sin 2 sin 1 2 cos 2 θ θ θ π σ r K_I x ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = 2 3 sin 2 sin 1 2 cos 2 θ θ θ π σ r K_I y ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = 2 3 cos 2 cos 2 sin 2 θ θ θ π τ r K_I xy

(16)

Figure 2.4 Stress components close to the crack tip

In two-dimensional problems, the stress in the z-direction can be divided into two different cases. First, if the plate is very thin the stress component σz can be assumed to

be zero. This state is called plane stress. Second, if the plate is thick the strain component εz can be assumed to be zero and the stress becomes σz = ν(σx + σy). This

condition is named plane strain or plane deformation. The out of plane shear stresses τxz

and τyz is zero in both cases.

0 =

z

σ for plane stress (2.4)

(

x y

)

z ν σ σ

σ = + for plane strain (2.5)

In a three-dimensional plate the material on the free surface is in a state of plane stress, because there are no stresses normal to the free surface. In the interior of the plate there are plane strain conditions near the crack tip [4]. This means that it is a region near the free surface where the stress state is neither plane stress or plane strain.

2.3 Stress intensity factor

The stress components in Equations 2.3a-c are proportional to a single constant, KI. If

this constant is known the entire stress distribution at the crack tip can be calculated from the stress equations. The constant K is called the stress intensity factor and defines the strength of the crack tip singularity 1/√r [4]. It can be seen as the magnitude of the stress field at the crack tip and this implies that the stresses near the crack tip increase in proportion to K. The subscript on K denotes the Mode of loading.

For a fix r and θ in the stress Equations 2.3a-c the only unknown parameter is KI. This

means that the stress intensity factor must be a function of both the applied load and the geometry of the cracked structure. So the stress intensity factor can be written as follows according to [5] τxy σx σy θ x y r Crack

(17)

(2.6) where σ is a stress, a is a measure of the crack length and f is a non-dimensional function of the geometry and load case.

Knowledge about the stress intensity factor is very important because this parameter will determine whether a crack grows or not and how fast it might grow. A lot of work has been done in determination of KI analytically for both two- and three-dimensional

problems. In the reference [6] many cases are solved, but in general three-dimensional problems are complex to solve analytically.

2.4 Energy consideration

An alternative approach in fracture mechanics to the stress intensity factor is the energy criterion. Griffith was the first to propose the energy criterion for fracture, but Irwin is primarily responsible for developing the present version of this approach [4]. First define the potential energy, П, of an elastic body [2]

(2.7)

where U is the elastic strain energy stored in the body and W is the potential energy or the work done by external forces. The energy release rate, G, is defined as the rate of change in potential energy with respect to the crack area [2].

(2.8)

The energy release rate is a measure of the energy available for an increment of a crack extension. To define the fracture criteria in an energy consideration the same discussion used for the stress intensity can be employed. An unstable crack growth occurs when the energy release rate reaches a critical value, G = Gc, where Gc is a measure of

fracture toughness and can be considered a material property.

2.4.1 Relationship between KI and G

So far two parameters have been introduced to describe the behaviour of cracks, first the stress intensity factor and then the energy release rate, and for a linear elastic material there is a connection between these two parameters. Lets consider an infinite plate with a through thickness crack of length 2a subjected to a uniform tensile stress. The energy release rate is given by Equation 2.9 and the same problem solved with the stress intensity factor is defined in Equation 2.10 [4].

(2.9) (2.10) f a KI =σ π W U− = Π dA d G=− Π E a G 2 πσ = a KI =σ π

(18)

where σ is the tensile stress, a is the crack length and E is Young’s modulus. Note that for a constant G and KI the stress varies with 1/√a in both cases. By combining Equation

2.9 and 2.10 the following relationship between KI and G can be shown for plane stress

and plane strain [4]

(2.11) E

E' = _{for plane stress} _(2.12)

2 '

1−ν

= E

E for plane strain (2.13)

where E is the Young’s modulus and ν is the Poisson’s ratio. It can also be proved that Equation 2.11 is a general relationship applied to all configurations in Mode I [4].

2.5 The concept of damage tolerance

The damage tolerance concept assumes that cracks exist in the structure but ensures that the complete structure will not fail prior to the time that the defect is detected and repaired. To realize this several methods are used, for example multiple load paths and crack arresters such as stiffeners. A crack arrester decrease the crack growth as the crack passes beneath each stiffener. The damage tolerance design approach requires periodic inspections in order to ensure detection of critical cracks before it can reach a critical length. The component may then be repaired or replaced.

Figure 2.5 Periodic inspections in order to ensure safe inspection interval. (From Saab)

' 2 E K G₌ I

(19)

There are two main tasks within the damage tolerance concept where fracture mechanics methods and the stress intensity factor play a central role. Fracture mechanics is used to predict the behaviour of macrocracks. A macrocrack is here defined as a crack from the size of about a millimeter which propagates mainly under opening displacements (Mode I). The tasks are to analyse how cracks propagates under static and cyclic loading.

2.5.1 Residual strength

Prediction of the largest crack a structure can sustain under specific residual strength requirements can be done using a critical value of the stress intensity factor when unstable crack propagation will occur.

The critical value KIc is called fracture toughness and is a material constant that is

independent of the size and geometry of the cracked body under certain conditions [4]. The fracture criteria can be expressed as:

Ic I K

K = (2.14)

2.5.2 Fatigue crack growth

Prediction of the number of cycles a fatigue crack is growing stable can be done using the stress intensity factor. It is assumed that the conditions at the crack tip are uniquely defined by the current KI and the crack growth rate is controlled by the range of the

stress intensity factor ΔKI = (KImax - KImin) where KImax and KImin are the stress intensities

at the peak and valley of a stress cycle.

The crack growth rate has a linear relation to the stress intensity factor range in a log-log representation and can be described by Paris’ law and can be written as [2]

(2.15)

where C and n are material constants determined by experiments. This relation can be used to calculate the number of constant amplitude loading cycles, N, required to reach a certain crack length, a.

(2.16)

2.6 The finite element method

Physical phenomena in engineering mechanics can be studied by math models which lead to differential equations. Usually the problem is too complicated to be solved by analytical methods so it is necessary to use a numerical technique for the solution. A widely used technique is the finite element method (FEM), where general differential equations can be solved in an approximate manner. Instead of seeking solutions that hold over the entire region this method divide the region into smaller parts, finite

(

)

n I K C dN da Δ =

(

)

∫

_Δ =ac ai n I K C da N

(20)

elements, and the approximation can be carried out over each element. This means that a variable can for example vary in a non-linear way over the entire region but assumed to vary linearly over each element. The behaviour for each element can then be determined because the approximation made over each element is fairly simple. The solution can then be improved by using more elements to represent the structure. The collection and arrangement of all elements is called a mesh.

In the two-dimensional case these elements are normally triangular, with three or six nodes, or quadrilateral, with eight or nine nodes. Each element is connected at nodal points to neighbour elements. At the nodal points compatibility conditions are satisfied and stress and displacement variations are restricted within each element. In stress analysis the solution process starts with determination of the stiffness matrix, which relates nodal forces to nodal displacements, for each element. Then construct an overall stiffness matrix and finally solve the problem by use nine overall matrices to obtain unknown forces and displacements.

Today when most analysis are performed by computers it is still important for the engineer to have knowledge about the concepts and theories used in the programs to understand the problem, how to model it, the behaviour of finite elements, assumptions and limitations built into the software and judge if the results are realistic or not.

2.7 Mesh design for crack problems

This is the most important theory part for discussions concerning crack tip subdomains in this thesis. The design of a finite element mesh is more or less individual to each problem to solve and some programs have automatic mesh generation capabilities, but fracture mechanics is one case that requires some special knowledge of the user. Common element types for crack applications are 8-node elements for two-dimensional problems and 15-node wedge elements or 20-node brick elements for three-dimensional problems, see Figure 2.6.

Figure 2.6 Common elements used for crack problems in three dimensions

2.7.1 Quarter points

An important question in crack analysis is how to split the region close to the crack tip into finite elements. As the theory for linear elastic conditions showed in Chapter 2.2,

20-node brick 15-node wedge

(21)

stresses near a crack tip are proportional to 1/√r, where r is distance from crack tip. To overcome this difficulty one might either refine the mesh in the crack tip region or construct elements that can represent this singularity. For two-dimensional problems a refined mesh is quite possible but in three dimensions it becomes more costly. It turns out, as will be shown later, that conventional elements become capable of representing a 1/√r stress field when the mid-side nodes are moved to quarter points around the crack tip, showed in Figure 2.7.

Figure 2.7 Mid-side nodes moved to quarter points around the crack tip

To illustrate this behaviour in a very simple way a three node bar element with the mid-side node moved to quarter point will be studied. An isoparametric formulation is applied and it use so called shape functions to map both coordinates and displacements between global coordinates and reference coordinates. If the same shape functions are used for both coordinates and displacements an element is called isoparametric.

Figure 2.8 shows a three node bar element with the mid-side node moved to quarter point, where (a) is the physical element with a global coordinate x and (b) is the mapped element with a reference coordinate ξ. The reference coordinate map the physical element into a reference element where node 1 ξ = -1, node 2 ξ = 0 and node 3 ξ = +1.

Figure 2.8 Three node bar element with global coordinate respective reference coordinate

For linear elastic conditions the stress is proportional to the strain, written in general form in 2.17, which means that the 1/√r behaviour also will be represented in the strain field. (2.17) Crack tip x 1 2 3 3L/4 L/4

(a) Physical element with global coordinate

1 2 3

ξ

ξ = -1 ξ = +1

(b) Mapped element with reference coordinate

(22)

where {σ} is the stress matrix, [D] is the stress-strain constitutive matrix and {ε} is the strain vector which is defined below.

(2.18) {d} contains the nodal displacements in global coordinates, for a three point bar {d} = {u1 u2 u3}T, and [B] is the strain-displacement matrix given by 2.19 for a bar element.

(2.19)

where Ni are shape functions corresponding to node i. The shape functions are

expressed in terms of ξ, so the chain rule must be used for the components in the [B] matrix.

(2.20)

Each shape function Ni have unit value at node i and is zero at every other node. The

shape functions for the three nodes on the bar are defined as follows according to [7].

(2.21a-c)

The first derivatives in the chain rule expression in 2.20 can now be calculated.

(2.22a-c)

The second term in the chain rule in 2.20 is not immediately available, so first define a global x-coordinate as a function of ξ [7].

(2.23)

where n is the number of nodes in the element and xi are coordinates in the global

system corresponding to node i. Use expression 2.23 with node coordinates x1 = 0, x2 =

L/4 and x3 = L. i n i i x N x

∑

= = 1

{ }

ε =

[ ]

B

{ }

d

[ ]

_⎢⎣⎡ _⎥⎦⎤ ∂ ∂ = ... ... x N B i dx d d dN dx dN_i _i ξ ξ =

(

)

(

)

(

2

)

3 2 2 2 1 2 1 1 2 1 ξ ξ ξ ξ ξ + = − = + − = N N N ξ ξ ξ ξ ξ ξ + = − = + − = 2 1 2 2 1 3 2 1 d dN d dN d dN

(23)

(

)

2 1 4 +ξ = L x or =2 −1 L x ξ (2.24)

Now calculate the second term in the chain rule expression, which have the same value at all nodes.

(2.25)

Finally use Equation 2.22, 2.24 and 2.25 to calculate the [B] matrix, which gives the axial strain of the bar as follows.

Equation 2.26 now shows that the strain in x direction varies with 1/√x. Since the stress is proportional to the strain for linear elastic conditions, elements with quarter points is capable of representing a 1/√x stress field. In the same way as described above for a bar element, it can be shown that a transfer of the mid-side nodes to quarter points give a 1/√x stress field for quadratic and triangular elements [7]. If 9-noded elements in two dimensions are used a new problem arise, namely where shall the extra node be placed to retain the stress singularity. There is no general answer to that question but [8] had proposed 3/8 or 11/32.

2.7.2 Wedge elements

In three dimensions the quarter point modification results in a 1/√r stress singularity on the element edges for brick elements and within the element and on the edges for wedge elements [4]. This means that wedge elements introduces the stress singularity on all rays emanating from the crack tip in the entire element. So the elements around the crack tip are degenerated to 15-node or 18-node wedge element and it is recommended that l<a/8 and 30°<α<40° according to [7], where a is crack length, l is length of wedge element and α is angle of wedge element, see Figure 2.9.

Figure 2.9 15-node wedge elements around the crack tip l/4 l a α Crack tip

[ ]

{ }

⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛₋ ₊ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − = = 3 2 1 2 1 2 2 4 2 3 2 u u u Lx L Lx L Lx L d B x ε (2.26) Lx dx d 1 = ξ

(24)

2.7.3 Curved crack front

Another aspect that has to be considered in three-dimensional cracks and specially in this work is the presence of circular or elliptic crack front. It has been shown in [8] that if the crack tip region is meshed with elements that have only straight lines, Figure 2.10a, the square root singular behaviour exists only on the two planes which form the outer surfaces to the crack. If instead a cylindrical element with four curved lines is applied, Figure 2.10b, the singular behaviour is present on all rays emanating from and within a small distance of the crack edge.

Figure 2.10 Element with straight lines compared with curved lines

2.7.4 Spider mesh

An efficient mesh design for the crack tip region and further away from the crack has been proved to be some kind of spider mesh [4], illustrated in Figure 2.11. This mesh consists of an inner ring with triangular elements around the crack tip and outer rings of four sided elements. Since the crack tip region contains very large stress and strain gradients it is obvious that the mesh refinement should be greater closer to the crack tip.

Figure 2.11 Efficient mesh design around the crack

(25)

2.8 FE-program Trinitas

The Trinitas program is an integrated graphical environment for finite element analysis, developed by Bo Torstenfelt at the division of solid mechanics at Linköping University. The program includes activities for geometry modelling, boundary condition definition, mesh generation, simulation and evaluation of the results [9].

2.8.1 Method for calculation of KI

The Trinitas program will be used for the crack analysis and in this chapter there will be a description of the technique for determination of the stress intensity factor used in Trinitas. The stress intensity factor is extracted from either the stress or displacement field or by any energy based method. There are three basic ways to calculate the stress intensity factor, namely the crack opening displacement method, the stiffness derivative technique and the J-integral method, where the latter two are energy based methods. The crack opening displacement method is simplest but least accurate and the other two methods have approximately equal accuracy [7]. The stiffness derivative method is used in the Trinitas program and can be applied to both two- and three-dimensional bodies. The stiffness derivative technique proceeds as follows for a cracked body subjected to Mode I loading. The potential energy П in terms of the finite element method is given by [7]

(2.27)

where {d} is the global nodal point displacement vector, [K] is the global stiffness matrix and {F} is the global load vector. The energy release rate G, defined in Equation 2.8 as the rate of change in potential energy with respect to the crack area may now be applied to Equation 2.27. The determination of G is evaluated for a two-dimensional body under fixed load conditions and the crack area is therefore replaced with the thickness times the crack length, ∂A = t∂a.

(2.28)

The first term must be zero since the stiffness equation is [K]{d} = {F}and the third term must also be zero because of a constant load. Finally, the energy release rate can be calculated and used in Equation 2.11, which defined the relationship between G and KI,

as follows.

(2.29)

where E’ is defined in Equation 2.12 and 2.13 for plane stress respectively plane strain. The stress intensity factor can now be computed numerically by first solve the stiffness equation [K]{d} = {F} for {d}, then increase the crack a small amount, determine a new stiffness matrix and finally approximate the stiffness derivative as ∂[K]/ ∂a ≈ Δ[K]/Δa. This is most efficiently done on the element level.

{ }

[ ]

{ }

d a t K d E K G I T ∂ ∂ − = = 2 1 ' 2

{ }

_{

_{[ ]}

_{{ } { }}

_{} { }}

[ ]

_{{ } { } { }}

a t F d d a t K d F d K a t d A G T T T ∂ ∂ + ∂ ∂ − − ∂ ∂ − = ∂ Π ∂ − = 2 1

{ }

d T

[ ]

K

{ } { } { }

d ₋ d T F = Π 2 1

(26)

3

The crack tip subdomain

This chapter begins with introducing two definitions that will be used in the continuation work, illustrated in Figure 3.1. The eight point brick volume will be generated to a mapped topology adjusted for accurate calculation of stress intensity factors for curved three-dimensional cracks, called the crack tip subdomain.

Figure 3.1 Two separated definitions

This part has been focused on how the eight point brick volume shall be decomposed to the crack tip subdomain in order to give accurate results of stress intensity factor calculations and have the ability to be applied at real applications. Later on in the thesis the complete workflow will be described and the crack tip subdomain will be introduced in a submodel.

First it has to be decided which type of geometry it should be on the replaceable volume in the submodel. The automated generation of the crack tip subdomain will be implemented in Trinitas which has three types of volumes available for 3D models, namely brick, wedge and tetra. The best geometry in this application is a brick volume, because it is common in 3D analysis, easy to introduce a crack topology and simple to cut with respect to symmetry.

(27)

To give accurate results and have the ability to be applied in real applications the crack tip subdomain has to fulfil a list of requirements that will be further discussed below.

• Replace an eight point brick volume in the submodel.

• Have smaller elements near the crack tip and wedge elements with quarter points close to the crack tip.

• Have three of the bricks six surfaces free from elements, since these surfaces may contain boundary conditions.

• Be built up with brick, wedge and tetra volumes, since this three volume types can be handled in Trinitas for 3D analysis.

• Adjust its topology with respect to the size of the eight point brick volume and the crack.

• Handle both circular and elliptic crack fronts.

• Manage to work in presence of radius, hole or groove.

3.1 Eight point brick volume

First of all the crack tip subdomain shall replace an eight point brick volume in the existing submodel. The area around the crack must therefore be built up with brick volumes. To analyse a corner crack two bricks must be replaced with the crack definition in between, for a surface crack four bricks shall be replaced and for an inner crack eight bricks have to be replaced. This means that the total number of elements easily can be decreased in the case of symmetry.

3.2 Elements close to the crack tip

As stated in the theoretical part in Chapter 2.7 the mid points on the elements close to the crack tip should be moved to quarter points to represent the stress singularity. It is also preferable if these elements are of wedge type for better accuracy. Furthermore it was declared that the entire mesh should be formed as some kind of spider mesh, with larger elements further away from the crack tip, see Figure 2.11.

3.3 Three surfaces free

An eight point brick volume has six surfaces and three of this shall preferably be free from nodes, lines and surfaces which results from refined elements. The reason for this is that the crack tip subdomain can be built in with other elements and have boundary conditions on these surfaces. It is also very good if the mesh can, as much as possible, be refined with Trinitas existing automated model mesh routine without spread elements over the three free surfaces. In Figure 3.2 a simple 2D case is illustrated, where (a) shows an arbitrary geometry there two elements shall be replaced with crack tip subdomains. Picture (b) shows a topology that fulfils this criterion and (c) shows an example that spread lines over all sides which are not preferred.

(28)

3.4 Three types of volumes

All these inner volumes that build up the whole crack tip subdomain is required to be made of maximum three different types of volumes. The ground for this statement is that the automated crack tip subdomain generation will be implemented in the Trinitas program which has three types of volumes available for 3D models, namely brick, wedge and tetra, showed in Figure 3.3.

Figure 3.3 Three volume types with corresponding refinement affect

The figures also demonstrate how a refinement of one line will affect the volume in Trinitas. It is important to be aware of how the volume type will be affected by automatic refinement, since this function should be available for the user. Note that the tetra will spread points over all lines, which means that the refinement is hard to control. So the conclusion is that the brick type should be used as much as possible to retain the control and ability to direct the refinement in the topology.

3.5 Adjust topology

Next point in the requirement list is the ability for the topology to adjust its points with respect to the eight point brick volume and the size of the crack. The best geometry for

(b) (c)

(a)

(a) Brick (b) Wedge (c) Tetra

(29)

the crack tip subdomain is a cube, but the automated function has to work for other geometries as long as they are built up with eight points. Of course there will be a limit for how large the differences can be between the width, depth and height. The second thing that will influence the crack tip subdomain is the crack size specified by the user. In this case there also will be a limit for how small or big the crack can be in relation to the eight point brick volume. Figure 3.4 below shows an example of a warped brick volume that should be handled by the automated crack tip subdomain generation.

Figure 3.4 A warped volume with dashed crack front

3.6 Curved crack fronts

The crack front is always curved which means that the crack can either be circular or elliptic. A circular crack is often used in the presence of a corner crack and an elliptic crack is employed in both surface and embedded cracks. The user shall specify the lengths of the two axes that form the crack front, but the relation between the axes can not be too large. This is normally not a problem, since in the case where elliptic cracks are applied the assumed relation of the axes is often two before the crack start to propagate.

3.7 Radius, hole and groove

The last demand is about curved lines on the replaced eight point brick volume, since stress concentrations will appear close to radii, holes and grooves. In the area with a stress concentration the probability for a crack initiation is large, so it is necessary for the crack tip subdomain to have the ability to handle this kind of geometries. Figure 3.5 shows an example with a groove where two volumes must be replaced to introduce a corner crack. The result of a cylindrical surface on a brick volume is that every line in the crack tip subdomain that lies on that surface must be curved.

(30)

3.8 Result

All these requirements resulted in three different types of basic topologies which can be mapped to a crack tip subdomain. The first one is the most simple and the third topology fulfils all requirements mentioned earlier. This chapter will give a description of all topologies and a comparison of the accuracy to calculate stress intensity factors will follow in Chapter 5.1.

3.8.1 Topology 1

The first topology has straight lines describing the crack front and brick element close to the crack tip, illustrated in Figure 3.6. The topology is built with 30 brick and 4 wedge elements. The refined view in (c) has still three surfaces free from elements.

Figure 3.6 Three different views of topology 1

(a) Over (b) Under

(31)

In this topology the crack front is curved and two extra lines around the crack front is applied, but the elements close to the crack tip are still brick elements, see Figure 3.7. The topology contains 52 brick and 4 wedge elements.

Figure 3.7 Three different views of topology 2

(a) Over (b) Under

(32)

In the last topology the elements around the crack front are replaced with wedge elements, Figure 3.8c, and this topology fulfil all requirements stated before. It is built with at the minimum 55 brick and 20 wedge elements.

Figure 3.8 Four different views of topology 3

(a) Over (b) Under

(33)

4

The crack tip subdomain implementation

Trinitas shall contain an automated support for the implementation of the crack tip subdomain. First the user has to specify the length of the two axes describing the curved crack front. Then pick an eight point brick volume and define the two lines describing the crackposition. The program will then create the crack tip subdomain built with points, lines, surfaces and volumes. All points will adjust its coordinates with respect to the crack size and the geometry of the picked volume. After the first crack tip subdomain is created the user will have the ability to copy this to nearby brick volumes. Note that the copied crack tip subdomain will use the same crack size but adjust its points to the new volume geometry. The crack front will automatically be a group and defined as a crack in Trinitas. The code is written as general as possible for the ability to add new topologies in the program, but still use the same subroutines. Here follows a description of some important subroutines used to implement the crack tip subdomain.

4.1 Get geometry mask

In this routine all points are defined with an 8-node 3D isoparametric formulation which uses shape functions to map the coordinates. All points are described in reference coordinates (-1 ≤ r, s, t ≤ +1) and will later on be mapped to global coordinates (x, y, z). Here are also all lines, surfaces and volumes defined for the topology, which must be done in order to mesh a model in Trinitas. Points build lines, lines build surfaces and surfaces build volumes. In the case of cylindrical surfaces on the eight point picked volume this subroutine will change all lines on cylindrical surface to curved lines. Finally, this function will call for modify shape coordinates.

4.2 Modify shape coordinates

This function will adjust all points with respect to the specified crack size and the picked volume geometry. The points close to the crack tip are very dependent on the crack size and points further away are more dependent on the brick volume geometry. In Appendix A an example is attached of a short piece from the code.

4.3 Define new points

This subroutine map all points from the reference coordinates to the global coordinates with 8-node 3D isoparametric shape functions. In the case of cylindrical surfaces the function use 20-node 3D isoparametric shape functions to describe the curved lines. After this routine all points, lines, surfaces and volumes can be created and the crack tip subdomain can finally be visualized on the screen.

(34)

5

Accuracy and verification

5.1 Influence of geometrical parameters

This part will investigate how different geometrical parameters will influence the accuracy of the stress intensity factor KI. The reference case used is a corner crack edge

bar, which are a common crack in 3D applications, showed in Figure 5.1.

Figure 5.1 Corner crack reference case

There does not exist any exact analytic solution to this reference problem, so lets study two different approximate solutions found in Appendix B and C. The first approximation found in Appendix B comes from Tada, Paris and Irwin [6] and is a very simple expression for the calculation of KI compared with the second solution in

Appendix C from Newman and Raju [4]. The first case is solved for an infinite body which can be simulated with the following input parameters: a = 0.1 m, c = 0.1 m, W = 1 m, t = 1 m and σ = 1 Pa. The result will be symmetric around 45° and can be seen in Table 5.1 with percent differences.

Angle Analytic KI

Newman, Raju Tada, Paris, Irwin Analytic KI %Diff.

10° 0.394 0.415 +5.3

20° 0.385 0.405 +5.2

30° 0.380 0.400 +5.3

45° 0.378 0.397 +5.0

Table 5.1 KI [N/m3/2] of two approximate analytic solutions of a corner crack a

c

(35)

The difference of the results is surprising big. In the following tests the Newman, Raju analytic approximation has been used, because it is more complex and comes from a newer reference.

In the following tests the crack tip subdomain will be refined until the solution of KI

converges. The height is not used in the analytic solution so it will be high enough to not affect the result and topology 3 is used if no other is mentioned. The KI values will

increase closer to the free surfaces, but on the two surfaces they must be zero. The analytic solution is therefore not valuable on the free surfaces and the Trinitas program makes an approximation of the elements on the free surface. So in the following analysis the KI values on the free surfaces (angle 0° and 90°) are neglected. In all these

tests no bending are applied, which means that KI will be symmetric around 45°.

5.1.1 Type of topology

Three different types of topologies have been introduced, but how much will the differences affect the result of the stress intensity factor. The input values of the parameters are: a = 0.1 m, c = 0.1 m, W = 1 m, t = 1 m, h = 2 m and σ = 1 Pa. The result as function of the angle is presented in Table 5.2 with percent differences.

Angle Analytic KI Top. 1 %Diff. Top. 2 %Diff. Top. 3 %Diff.

10° 0.394 0.405 +2.8 0.405 +2.8 0.407 +3.3

20° 0.385 0.390 +1.3 0.392 +1.8 0.394 +2.3

30° 0.380 0.348 -8.4 0.383 +0.8 0.386 +1.6

45° 0.378 0.358 -5.3 0.380 +0.5 0.383 +1.3

Table 5.2 KI [N/m3/2] of the analytic solution compared with the three topologies

The large differences in topology 1 are caused by larger elements close to the crack tip and that the crack front is built with straight lines. Both topology 2 and 3 shows good results and both have better accuracy closer to the centre. It is hard to say why number two shows smaller errors in this case, but keep in mind that the analytic solution is an approximation. It can also be mentioned that the elements close to the crack tip with quarter points must be square in the normal plane to the crack front, for topology 1 and 2. This phenomenon must the user have in mind if the mesh is refined around the crack tip. In topology 3 a refined mesh will automatically improve the elements around the crack tip, because all these elements are of wedge type.

5.1.2 Geometry of the eight point brick volume

This part will look into the relations between the width and height of the eight point brick volume, see Figure 5.2. The influence of the depth will not be analysed because of symmetry. The input values of the parameters are: a = 0.1 m, c = 0.1 m, W = 1 m, t = 1 m, h = 2 m and σ = 1 Pa. The result is showed in Table 5.3.

(36)

Figure 5.2 Eight point brick volume with dashed crack front

Angle Analytic KI H = 0.2 _{L = 0.2} %Diff. H = 0.1 _{L = 0.2} %Diff. _{L = 0.16}H = 0.2 %Diff.

10° 0.394 0.407 +3.3 0.405 +2.8 0.402 +2.0 20° 0.385 0.394 +2.3 0.391 +1.6 0.390 +1.3 30° 0.380 0.386 +1.6 0.382 +0.5 0.384 +1.1 45° 0.378 0.383 +1.3 0.379 +0.3 0.382 +1.1 60° 0.380 0.386 +1.6 0.382 +0.5 0.385 +1.3 70° 0.385 0.394 +2.3 0.391 +1.6 0.394 +2.3 80° 0.394 0.407 +3.3 0.405 +2.8 0.407 +3.3

Table 5.3 KI [N/m3/2] of the analytic solution compared with different geometries

It can be seen that lower height on the eight point brick volume will give lower KI and

smaller length will result in lower KI for points affected of the decreasing length, but the

results are still very good. For larger differences in height and length the topology could not be drawn for this size of crack.

5.1.3 Relation between crack size and the eight point brick volume

The size of the brick volume will be changed, meanwhile the crack size is holding constant. All three sides of the brick volume have the same length. This will give a relation of the crack size with respect to the brick volume. The input values of the constant parameters are: a = 0.1 m, c = 0.1 m, W = 1 m, t = 1 m, h = 2 m and σ = 1 Pa. Figure 5.3 illustrates three different relations where the crack length divided by the brick volume length is presented in percent. The result follows in Table 5.4.

Figure 5.3 Relations between crack length and brick volume

L H

(37)

Angle Analytic KI 10 % %Diff. 50 % %Diff. 75 % %Diff.

10° 0.394 0.409 +3.8 0.407 +3.3 0.399 +1.3

20° 0.385 0.395 +2.6 0.394 +2.3 0.388 +0.8

30° 0.380 0.386 +1.6 0.386 +1.6 0.381 +0.3

45° 0.378 0.383 +1.3 0.383 +1.3 0.379 +0.3

Table 5.4 KI [N/m3/2] of the analytic solution compared with different crack/brick relations

The conclusion is that smaller crack size in relation to the converted brick volume will give higher KI values and with that worse result. It is surprising that the last case with

75 % crack/brick relation shows so small errors even though many elements outside the crack front are really badly shaped.

5.1.4 Refinement of the crack tip subdomain

Here will the existing automated mesh refinement function in Trinitas be used. This part will investigate how the total numbers of elements in the crack tip subdomain will affect the result for both topology 2 and 3. The first test has the smallest number of elements possible for the required points along the crack front. The third test is an example where the three free surfaces are divided into smaller elements. Input values are: a = 0.1 m, c = 0.1 m, W = 1 m, t = 1 m, h = 2 m and σ = 1 Pa. The result for topology 2 is presented in Table 5.5 and for topology 3 in Table 5.6.

Angle Analytic KI 118 el. %Diff. 284 el. %Diff. _{div. sur.}644 el. %Diff.

10° 0.394 0.402 +2.0 0.405 +2.8 0.407 +3.3

20° 0.385 0.389 +1.0 0.392 +1.8 0.395 +2.6

30° 0.380 0.381 +0.3 0.383 +0.8 0.388 +2.1

45° 0.378 0.378 ±0.0 0.380 +0.5 0.385 +1.9

Table 5.5 KI [N/m3/2] of the analytic solution compared with number of elements for topology 2

Angle Analytic KI 157 el. %Diff. 407 el. %Diff. _{div. sur.}824 el. %Diff.

10° 0.394 0.408 +3.6 0.407 +3.3 0.409 +3.8

20° 0.385 0.395 +2.6 0.394 +2.3 0.397 +3.1

30° 0.380 0.387 +1.8 0.386 +1.6 0.391 +2.9

45° 0.378 0.385 +1.8 0.383 +1.3 0.388 +2.6

Table 5.6 KI [N/m3/2] of the analytic solution compared with number of elements for topology 3

The conclusion for topology 2 is that more elements give higher values of KI.

Additional elements should give a more accurate result which not seems to correspond in this test, but have in mind that the real solution probably lies a little above this

(38)

analytic approximation. In the other case, when topology 3 was applied, the results are more equal and independent of the number of elements. This shows that topology 3 is more stable and can for sure be used with fewer elements without affect the result. Finally, it can be concluded that keeping three of the six surfaces in the crack tip subdomain free from divided elements not affect the result so much, according to the third test for both topologies.

5.1.5 Serendipity or Lagrange type model mesh

There are two major types of finite elements available in Trinitas for 3D analysis. First, there is the Serendipity type with a 20-node brick element and a 15-node wedge element. The second is the Lagrange type with the higher order elements, the 27-node brick and the 18-node wedge. One advantage with the Lagrange type is that less elements are needed and one disadvantage is that the analysis time is longer. So therefore this test will be carried out with two different numbers of elements. Input parameters are: a = 0.1 m, c = 0.1 m, W = 1 m, t = 1 m, h = 2 m and σ = 1 Pa. The result is showed in Table 5.7, where the last row is the CPU analysis time in Trinitas.

Angle Analytic _K

I

157 el.

Ser. %Diff. 157 el. Lag. %Diff. 407 el. Ser. %Diff. 407 el. Lag. %Diff. 10° 0.394 0.408 +3.6 0.409 +3.8 0.407 +3.3 0.407 +3.3 20° 0.385 0.395 +2.6 0.395 +2.6 0.394 +2.3 0.394 +2.3 30° 0.380 0.387 +1.8 0.388 +2.1 0.386 +1.6 0.386 +1.6 45° 0.378 0.385 +1.8 0.385 +1.8 0.383 +1.3 0.383 +1.3

Time 14 s 55 s 25 s 100 s

Table 5.7 KI [N/m3/2] of the analytic solution compared with Serendipity and Lagrange type

The results of these tests are that the Serendipity finite elements can be used without loosing in accuracy. The use of the Serendipity type instead of the Lagrange will reduce the analysis time, which can be very useful when larger models are analysed.

5.1.6 Conclusions

This geometric parameter influence analysis has showed that both topology 2 and 3 give accurate results, but topology 3 is more stable with respect to number of elements in the crack tip subdomain. It has also appeared that the Serendipity element is recommended, because it gives good results in ¼ of the analysis time compared with Lagrange type. It can be noted that all results lies above the analytic solution. One reason for this might be that there not exists any exact analytic solution for the corner edge crack problem.

(39)

5.2 Reference cases

In this chapter, it will be investigated how the crack tip subdomain will behave in four different reference cases, namely corner crack, surface crack, embedded crack and finally a corner crack close to a hole. These kinds of tests are very important to confirm that the crack tip subdomain can be used and trusted in real applications. Both topology 2 and 3 will be studied and in all tests the mesh will be refined until the solution converges. As said before the KI values will increase closer to the free surfaces, but on

free surfaces they must be zero. The analytical solution is therefore not valuable on the free surfaces and Trinitas makes an approximation over the last element close to the surface, so KI values on free surfaces are neglected.

5.2.1 Corner crack

The first case is a corner crack edge bar which has been used in Chapter 5.1. This test will include both tensile stress and bending stress, so the Newman, Raju reference case in Appendix C has to be used. Two different sets of input variables will be studied, first with tensile stress and then with added bending stress. The input values for the first case are: a = 0.1 m, c = 0.1 m, W = 1 m, t = 1 m, σm = 1 Pa, σb = 0 Pa and the eight point

brick volume is a cube with length 0.2 m. The result can be seen in Table 5.8 for topology 2 and 3. In Table 5.9 the bending stress is applied, σb = 1 Pa.

Angle Analytic KI Top. 2 %Diff. Top. 3 %Diff.

10° 0.394 0.405 +2.8 0.407 +3.3 20° 0.385 0.392 +1.8 0.394 +2.3 30° 0.380 0.383 +0.8 0.386 +1.6 45° 0.378 0.380 +0.5 0.383 +1.3 60° 0.380 0.383 +0.8 0.386 +1.6 70° 0.385 0.392 +1.8 0.394 +2.3 80° 0.394 0.405 +2.8 0.407 +3.3

Table 5.8 KI [N/m3/2] for reference case with a corner crack

10° 0.767 0.790 +3.0 0.795 +3.6 20° 0.745 0.759 +1.9 0.763 +2.4 30° 0.730 0.734 +0.5 0.739 +1.2 45° 0.718 0.720 +0.3 0.726 +1.1 60° 0.716 0.720 +0.6 0.725 +1.3 70° 0.722 0.732 +1.4 0.736 +1.9 80° 0.736 0.754 +2.4 0.758 +3.0

(40)

First it can be noted that the result is not symmetric around 45° when the bending is applied in the second case. Furthermore the result of the analysis show that the analytic solution lies below the tests for all points. This can partly be explained by variations in solutions for corner crack reference cases. As shown in Chapter 5.1 a different reference gave about 5 % higher values than Newman, Raju. It can also be seen that topology 3 always gives a higher KI compared to topology 2. The results are worse closer to the

free surfaces, but it is still in reasonable limits.

5.2.2 Surface crack

Here a surface crack in an infinite body will be studied, see Appendix D for the Newman, Raju analytic solution. The input data in the first test are a = 0.1 m, c = 0.1 m, W = 1 m, t = 1 m, σm = 1 Pa, σb = 0 Pa and the cubic brick volume length is 0.2. The

stress intensity factor will be symmetric around 90°, so the result for topology 2 and 3 in Table 5.10 goes from 10° to 90°. In the second test a bending stress is added, σb = 1 Pa,

see result in Table 5.11.

10° 0.399 0.403 +1.0 0.406 +1.8 20° 0.389 0.388 -0.3 0.390 +0.3 30° 0.382 0.377 -1.3 0.379 -0.8 45° 0.375 0.369 -1.6 0.372 -0.8 60° 0.373 0.367 -1.6 0.369 -1.1 70° 0.372 0.367 -1.3 0.369 -0.8 80° 0.372 0.368 -1.1 0.370 -0.5 90° 0.372 0.368 -1.1 0.370 -0.5

Table 5.10 KI [N/m3/2] for reference case with a surface crack

10° 0.775 0.787 +1.5 0.792 +2.2 20° 0.752 0.751 -0.1 0.756 +0.5 30° 0.732 0.722 -1.4 0.727 -0.7 45° 0.712 0.700 -1.7 0.706 -0.8 60° 0.701 0.689 -1.7 0.694 -1.0 70° 0.697 0.685 -1.7 0.689 -1.1 80° 0.695 0.684 -1.6 0.688 -1.0 90° 0.695 0.684 -1.6 0.687 -1.2

Table 5.11 KI [N/m3/2] for reference case with a surface crack including bending

The results differs in both too big and too small values and topology 2 still shows a little smaller values in all points compared to topology 3. In general topology 3 shows improved results with respect to topology 2. When the bending stress is applied the errors becomes a little bigger, but all results from both topologies are still sufficiently good.

(41)

5.2.3 Embedded crack

The reference case used in this section is an embedded elliptical crack. For an infinite body this case can be solved exact analytically, see the Tada, Paris, Irwin exact solution in Appendix E and the Newman, Raju approximate solution in Appendix F. This case will also include a comparison between these two analytic solutions. The example is not so common to analyse in real applications, but it is one case that can be solved exact in 3D. The input values of the parameters are: a = 0.05 m, c = 0.1 m, W = 1 m, t = 1 m and

σ = 1 Pa. The angle for an elliptic crack is measured as showed in Figure 5.4 and in an embedded crack there will not exist any free surface, so all angles are possible to use. The result is presented in Table 5.12 and it is symmetric in every 90°.

Figure 5.4 Embedded elliptical crack

Angle Analytic KI

Tada (exact) Analytic KNewman I % Diff. Top. 2 % Diff. Top. 3 % Diff.

0.0° 0.231 0.232 +0.4 0.232 +0.4 0.233 +0.9 19.0° 0.248 0.248 ±0.0 0.241 -2.8 0.243 -2.0 38.0° 0.280 0.280 ±0.0 0.263 -6.0 0.266 -5.0 54.6° 0.304 0.305 +0.3 0.285 -6.2 0.288 -5.3 67.7° 0.318 0.318 ±0.0 0.303 -4.7 0.307 -3.5 78.9° 0.325 0.325 ±0.0 0.316 -2.8 0.321 -1.2 90.0° 0.327 0.328 +0.3 0.321 -1.8 0.326 -0.3

Table 5.12 KI [N/m3/2] for reference case with an elliptic embedded crack

First it can be concluded that the two analytic solutions correspond very well, which was not the case for the corner crack used before. Now the differences become bigger, for the analysed model than before, and it is probably caused by the introduced elliptic crack front. The end points correspond very well to the analytic solution but the inner points are worse. It can also be seen that topology 3 shows better results than topology 2. Note that the stress intensity factor is greater at the minor axle which means that the crack will grow in that direction.

a

c θ

(42)

5.2.4 Corner crack close to a hole

The last reference case is obtained from Börje Andersson at FOI (Swedish Defence Research Agency). The relative error in the achieved results is about 0.1 %. The example contains a large plate with a small hole in the centre and applied to uni-axial tension stress, see Figure 5.5. The hole will lead to a stress concentration and close to the hole two corner cracks are applied, illustrated in Figure 5.6a. This reference case is the most complicated, since it introduces both stress concentrations and the need for elements with one cylindrical side.

Figure 5.5 Large plate with a small hole applied with stress

Figure 5.6 Reference case from FOI

In this case only topology 3 will be used, since it is the best topology for elliptical cracks. Three different sets of input values will be analysed see Table 5.13 below and for comparison of results the coordinate system in Figure 5.6b is used. All results including percent differences can be found in Table 5.14a-c.

σ σ H W t R a c Crack c x y a

(43)

a c R t W H σ

Input 1 200 100 400 400 40000 40000 1

Input 2 100 100 200 200 20000 20000 1

Input 3 50 100 100 100 10000 10000 1

Table 5.13 Three sets of inputs

x y KI FOI KI input 1 % Diff.

99.2 25.4 35.3 34.8 -1.4 96.7 50.7 34.7 34.0 -2.0 92.5 75.8 34.1 33.0 -3.2 81.5 115.8 33.3 31.7 -4.8 61.5 157.6 32.2 30.5 -5.3 45.7 177.9 31.7 30.5 -3.8 25.1 193.6 31.7 31.3 -1.3

x y KI FOI KI input 2 % Diff.

98.4 17.6 24.3 23.8 -2.1 93.8 34.7 23.6 23.3 -1.3 86.2 50.7 23.6 23.2 -1.7 70.7 70.7 24.2 23.9 -1.2 50.7 86.2 25.7 25.3 -1.6 34.7 93.8 27.6 27.5 -0.4 17.6 98.4 30.4 30.5 +0.3

X Y KI FOI KI input 3 % Diff.

96.8 12.6 14.9 14.2 -4.7 88.9 22.8 15.5 14.6 -5.8 78.8 30.8 16.7 15.7 -6.0 57.9 40.8 19.0 17.9 -5.8 37.9 46.3 21.6 20.5 -5.1 25.3 48.4 23.7 23.2 -2.1 12.7 49.6 26.5 26.2 -1.1

Table 5.14 KI for a reference case with a crack close to a hole

(a) Result for input 1

(c) Result for input 3 (b) Result for input 2

Analysis of Three-Dimensional Cracks in Submodels

DEPARTMENT OF MANAGEMENT AND ENGINEERING