The finite element method vs the edge function method for linear fracture analysis

Engineering Fmcfm Mechanics Vol. 13, pp. 4345 Pergamon Press Ltd., 1960 Printed in Great Britain

THE FINITE ELEMENT METHOD VS THE EDGE FUNCTION METHOD FOR LINEAR FRACTURE ANALYSIS J. F. FLEMING, J. J. GUYDISH, J. R. PENTZ and C. E. RIJNNION School of Engineering, University of Pittsburgh, Pittsburgh, PA 15261,U.S.A. and G. P. ANDERSON Thiokoi Corporation, Wasatch Division, Brigham City, UT 84301,U.S.A. Abstract-An investigation was performed to compare the Finite Element Method and the Edge Function Method for the solution of two dimensional linear elastic fracture problems. The comparison was made to determine whether one of the procedures had a distinct advantage over the other for fracture analysis. Although each method has minor advantages for the computation of specific quantities or the analysis of specific types of problems, there is no major advantage for either procedure which is immediately apparent.

Because of the considerable investment which many organizations have already expended on developing a Finite Element capability, and the additional investment which would be required to imple~nt the Edge Function Method the application of the Edge Function Method will probably be very limited.

IF STRUCTURES, in which fracture considerations can be a controlling factor, are to be designed for both safety and economy the design engineer must have information concerning the effect of cracks upon the strength of the structure. For very simple geometries, analytical solutions have been formulated for the elastic stress distribution around a crack. In fact, for many situations, design tables have been developed. For more complex geometries, however no tabulated solutions are available, thus the designer is required to resort to numerical calculation procedures to obtain an approximate solution with an accuracy acceptable for engineering design. The purpose of this discussion is to compare two such numerical procedures, the Finite Element Method and the Edge Function Method, for application to two dimensional problems of linear elastic fracture mechanics. The Finite Element Method is a well known numerical technique, suitable for both elastic and elastic-plastic problems, which relies upon simulating the structure by a group of small elements with the accuracy of the solution, in general, increasing as the size of the elements are decreased. The Edge Function Method is a technique, which is restricted to elastic problems at the present time, in which large ma’cro, as opposed to micro, elements can be used because the governing elastic equations are completely satisfied for the domain under analysis. Both procedures are computer oriented and require the solution of a large set of linear simultaneous equations. The actual computation process for either procedure can be completely automated, thus the main concern of the analyst is: developing a mathematical model to represent the problem to be analyzed; preparing input data for a computer program; and interpreting the final computer output. The present comp~son has been made for a series of test cases to obtain an indication of situations in which the Edge Function Method might possess an advantage, in spite of its limitation to linear elastic solutions, recognizing the present large investments many organizations have made in developing their Finite Element capabilities. The parameters considered in the comparison of the two procedures include: the ease with which a problem can be set up for analysis on the computer; the information available from the computer analysis; the amount of post processing of output required to arrive at useful engineering design information; the cost of the analysis in terms of both computer resources and man hours required for data prep~ation and output inte~retation; and the accuracy of the final results. Both single material and b&material problems are considered.

APPLICATION OF THE FINITE ELEMENT METHOD

The Finite Element Method (FEMJ is a numerical procedure in which a continuum is represented by a mathematical model consisting of a number of individual elements joined together at pre-determined nodes. The displacements at the nodes are computed by solving a set of linear simultaneous equations obtained from equilib~um requirements. In these equations 43

44

.I. F. FLEMING et al.

the displacements at the nodes are related to the external applied loads by the combined stiffnesses of the individual elements connected to the nodes. In matrix form the equilibrium equations are

where [K] is the total stiffness matrix for the system, which is obtained from the sum of the stiffnesses of the individual elements; (0) is the vector of displacements at the element nodes; and {W} is the vector of external applied load at the nodes. This approach has been applied to a wide variety of problems in solid mechanics including: plane stress and plane strain analysis; plate bending; thin and thick shell analysis; and three dimensional stress analysis. In this discussion the interest is solely in the application of the FEM to the elastic analysis of two Dimensions structural com~nents confining cracks and in particular to the determination of the stress dist~bution around the crack. in the application of the FEM to any problem is solid mechanics the accuracy of the analysis is highly dependent upon the mathematical model which is used to represent the physical situation. The solution of eqn (1) will result in the exact nodal displacements for this approximate mathematical model. Therefore, it is important that the mathematical model be chosen with physical properties that cIosely represent those of the structure being analyzed and that the array of elements be chosen to accomodate any anticipated deformation pattern or stress dist~bution. The irn~~nt variables in the mathematical model are: the size, shape and distribution of the individual finite elements; the assumed displacement distribution within the individual elements; the manner in which the restraints for the continuum are simulated at discrete nodal points; and the manner in which the external loads are applied. In the analysis of plane stress or plane strain problems, which includes the problems of interest here, a number of different types of elements have been developed. The most commonly used elements are of triangular or quadrilateral shape with either a linear or quadratic displacement distribution within the element. Higher order elements have been used in some special situations. The FEM has been used in fracture analysis in a number of previous investigations using a wide variety of elements. The analyses reported here were limited to the use of four element types: three node constant strain triangular elements (CST); six node linear strain triangular elements (LST); four node quadrilateral elements formed from four CST elements by condensing out the center node (CSQ); and eight node linear strain isoparametric quadrilateral elements (LSQ). The latter element is often refered to in the literature as 8-NIQ. These elements are available in a wide variety of commercial computer programs, all of which perform essentially the same computations during the anaiysis. The difference between the programs is prima~ly in the input required and the form of the output. Several speciaiized crack tip finite elements have been developed which include the effect of the elastic stress singularity at the tip of the crack[l, 21. Although these specialized elements are a very important contribution to fracture analysis and can decrease the FEM computational effort for a given accuracy, they were not used in any of the analyses performed in this investigation because of their lack of general availability in commercial finite element computer programs. The well known technique, discussed by Henshell and Shaw[3], of moving the mid-side node of the LSQ element nearest the crack tip, to the quarter point of the side, was also not used, For the analyses performed in this investigation automatic mesh generators were used to generate the finite element. arrays, from a minimum amount of input information, thus greatly reducing the input data preparation requirements. The use of the mesh generators permitted the analysis of many more individual problems than would have been possible without them, although the finite element arrays obtained might not have resulted in an absolute optimization of computer run time. The li~tations of the mesh generators required the insertion of nodes which might not have been present if an optimum manual generation procedure were used. It was felt that any trade off between data preparation time and computer run time was well justified. Since this investigation consisted of a comparison of the application of the FEM and the EFM in linear fracture mechanics, it was important for a realistic comparison that the FEM be

The finite element method vs the edge function method

45

used as eficiently as practicabIe in the solution of crack problems. Two questions which immediately were encountered where: (1) What should be the properties of the finite element mesh for each of the elements being used in the investi~tion?; and, (2) IIow should the computed results by compared? As to the first question, it was decided that one of the first steps in the investigation would be a study of the efficiency of various linite element meshes with an aim toward mesh optimization, within the limitations of the mesh generators being used. Therefore, an extensive study was performed to investigate mesh parameters for several different crack problems for which solutions had been previously reported. With respect to the second question, it was realized that vacuity would be incountered in attempti~ to compare directly the stresses in the region of the crack tip for the various problems invest~ated due to the singularity in stresses at that point. Therefore, it was decided that two computed quantities would be used for comparison, the stress intensity factor and the total strain energy. Both of these quantities are easily computed from the results of both the FEM and EFM analyses and both are also quantities which have been reported in the literature for a number of standard problems. In addition, the minimum energy principle furnishes a criterion for a comparison of overall accuracy. FEM mesh o~timi~at~o~ for crack ~rob~ern~ As stated previously, the FIB analyses in this investigation were limited to the use of elements with displacement distributions which were-either linear (constant strain) or quadratic (linear strain). Studies have shown that a high degree of mesh refinement is required to obtain acceptable solutions with these elements near singular points, particularly the linear displacement element. Further, the accuracy of the solution varies greatly with the size and shape of the elements and with the general configuration of the mesh. For example, it is known that very small elements are required near the crack tip, however, in order to reduce the number of nodes in the mesh, and the resulting number of degrees of freedom in the mathematics model, it is common practice to increase the size of the elements away from the crack tip. There are many methods of doing this, and consequently, individual analysts, working independently on the same crack problem, are likely to prepare a variety of meshes, all different in appearance. Unfortunately, these different meshes can also result in very different solutions. Evidence shows that mesh co~guration can have a significant effect on the accuracy of the computations. A study was performed as part of this invest~ation~for several simple single material crack problems, in order to establish guidelines for the creation of efficientmeshes using the CSQ and LSQ elements. Quadrilateralelements were used rather than triangular elements solely because of the availabilityof an automatic mesh generator for these elements during the initial phase of the investigation. At a later phase in the investigation a mesh generator was also developed for triangular elements, after which both quad~late~~ and triangular elements were used. The study was, by necessity, somewhat limits in scope and only a few elemen~y problems~with well documented analytical solutions, were considered. The accuracy of the solution for each of the meshes analyzed was judged by comparing the Mode I stress intensity factor, KI, obtained from the finite element solutions, to that tabulated by Tada et aL[4]. Since the accuracy obtained for the stress intensity factor also varies with the manner in which it is computed from the results of the finite element solution, a secondary aim of the study was a comparison of the two most commonly used methods of performing these calculations. These methods, which have been well documented in the literature, are: the displacement method[5], in which the stress intensity factor is found by ext~~lati~ the value from the computed displacements along the crack back to the crack tip; and the energy release rate method%], in which the rate of change in strain energy with respect to crack length, which is mathematicahy reiated to the stress intensity factor, is determined from the computed strain energies for two slightly different lengths of crack. The majority of the cases analyzed in deducing optimum mesh parameters were for a single edge crack in a long strip of linear elastic material with various crack length to specimen width ratios, subjected to a u~orm tensile end stress, as shown in Fig. 1. Due to symmetry, only half of the specimen was modeled by the finite element mesh. A typical mathematics model demonst~ti~ the general form of the quad~late~ finite element mesh used for the analysis of

J. F. FLEMING et al.

46

b Fig. l.Typical crackedspecimen,

these cases is shown in Fig. 2. Mathematical models with the number of nodes ranging from less than one hundred to over 1500 were considered. A few cases were also analyzed for a single edge crack specimen subjected to bending and a double edge crack specimen subjected to tensile loads. For the double edge crack problems, both axes of symmetry were utilized thus requiring only one quarter of the specimen to be modeled. In alt, over 200 finite element analyses were performed in this phase of the investigation.

crack

hp

-c1

Fig. 2. Typical single material finite element mesh.

The finite element methodvs the edge function method

41

A complete description of this study is given in a recently published paper by Guydish and Flemmingf’l]. The description includes: a discussion of the work of previous investigators; the manner in which the finite element mesh plasters were varied; the procedures used for computing the stress intensity factors from the finite element solutions; and the conclusions reached from the study. only the conclusions will be discussed here. The first inportant conclusion is that, at least for the elementary singularity problems studied, the finite element method with linear displacement elements can be used to obtain a solution that is for all practical purposes correct. There is a certain minimum error, however, that appears to be due mainly to difficulties with the indirect methods of computing Kl, rather than any weakness in the finite element method itself. Values within 1% of the analytics solution were obtained, using the displace~nt method to compute KI, with meshes of 500 nodes or less for each specimen considered. Note that due to the use of symmetry this corresponds to a much smaller number of nodes than would be required if the entire specimen were considered in the mathematical model. Further, 4% error or less was consistently obtained with meshes of less than 150 nodes. This seems to be about the minimum usable number of nodes for these problems, as the accuracy falls off rapidly below this point. Nevertheless, the evidence shows that the displacement method can be relied upon to give good estimates of KI as long as there are enough nodes along the crack to result in a reasonable exhalation of the value of ZU to the crack tip. If too few nodes are present along the crack considerable error can be introduced. The energy release rate method gave somewhat less accurate results in this study, however, it is believed that this is due to inconsistencies in applying the method rather than any weakness in the method. The differences in the results obtained by the two methods were usually well within an error level consistent with engineering design requirements. The second conclusion is that the use of quadratic displacement elements produces a considerable improvement in the computed value of KI. With these elements, the error was consistently 4% or less for a mesh with only 75 nodes; however, this advantage may be negated by the increase in computer time required. For the quadrilateral element programs used in this study there was approximately a linear relationship between the required computer execution time and the number of degrees of freedom, for the problem size ranges considered, for both the linear displacement and quadratic displacement elements. The quadratic displacement element program, however, required more than twice as much time as the linear displacement element program for a given number of degrees of freedom. This might or might not be true for other computer programs used by other analysts since program execution times can vary greatly depending upon the efficiency of the computational algorithms used, the bandwidth of the stiffness matrix, and the size of the core available in the computer. It was evident throughout the study that the accuracy of the solution was highly sensitive to mesh co&uration. The tinal conclusion, then, is the establishment of several parameters which are important to efficient mesh design. They are: (a) minimum element size; (b) node spacing in the near tip region; and (c) the slenderness ration of individual elements. In primula, the following sidelines were formulated. (a) It is necessary that the elements immediately adjacent to the crack tip be very small in proportion to the crack length, regardless of the manner by which KI is computed from the finite element results. There appears to be an optimum size of about 0.005 of the crack length for the linear dimensions of these elements. (b) The most efficient mesh results from a smooth progression of node spacing, starting with the minimum element size at the crack tip, and extending over a distance of at least one crack length in all directions. The progression may start within a few nodes of the crack tip since a large number of small elements is not needed as long as the profession is smooth. Beyond the distance of one crack length there is no apparent li~tation on element size. A geometric progression was used in this study, however, an arithmetic progression, as suggested by Anderson[6], may be somewhat more effective. In the arithmetic progression, the location of the nodes, relative to the crack tip, vary in the same order as the square root of I function which describes the displacement field, where r is the radial distance from the crack tip. (c) Very lbng slender elements anywhere in the mesh produce noticeable errors in the results. For best results, the length to width ratio of individu~ elements should not exceed approximately 50.

J. F. FLEMING

48

ef

al.

Using the above as guidelines for mesh generation, excellent accuracy and efficiency have been achieved in calculating the stress intensity factor. It must be emphasized, however, that these were were very simple single material problems, nevertheless, there is no apparent reason why these guidelines should not be valid for most other cases. FE&f u~~ljcati5~s to bi-mate~al ~r5blems There are many physical situations in which two or more materials are bonded together to form a single structural element. An important problem is that of adhesive debonding along a bi-material interface, especially in the vicinity of the crack tip. There is added mathematical complexity for this case compared to the single material case since the order of the singularity depends upon the physical properties of the bonded materials and an oscillatory factor of the form cos(lnr) is introduced[8]. A few analytical solutions for this type of problem have been reported in the literature[9]. Since the purpose of this discussion is a comparison of the application of the FEM and the EFM to fracture problems, the main interest here is how the accuracy of the two methods might be compared for bi-material situations. It was decided that, for this investigation, the stress intensity factor would again be used, however, special problems were encountered for the bi-material problem. In cracked homogeneous bodies the stress intensification near the crack tip arises solely from a geometric discontinuity. For pure symmetric or pure antisymmetric loading a single stress intensity factor can be defined, however, for bi-material bodies with a crack along the material interface, the stress intensi~c~tion occurs due to both a geometric discontinuity and a material discontinuity. It has been shown, for example, by Erdogan[9] that, for this situation, Mode I loading will result in both KI and KU stress intensity factors, thus for bi-material problems it is necessary to extract both stress intensity factors from the results of a finite element analysis. The technique used was a semi-graphical procedure which relates the complex stress intensity factor, KI + iK& to the complex displacements along the crack[ 101.As long as a sufficient number of nodes in included along the crack in the finite element mesh this procedure will give accurate values. The finite element analyses for the bi-mate~al crack problems were performed using CST and LST elements. The meshes were developed using the mesh generation guidelines discussed previously for single material problems. Due to the lack of symmetry on either side of the

interface _

material and crack

Fig. 3. Typical

bi-material

finite element

mesh.

The finite element method vs the edge junction method

49

crack, resulting from the two different materials, the finite element meshes used for the bi-material problems required more modes than those used for the single material problems. A typical finite element mesh for a b&material single edge crack problem is shown in Fig. 3. APPLICATION OF THE EDGE FUNCTION METHOD The Edge Function Method (EFM) is a numerical procedure which can be applied to a wide variety of linear Sunday value problems. The basic philosophy of the procedure is that an analytical function corresponding to the solution of a boundary value problem is composed of a linear combination of one or more of a group of functions which individually satisfy the governing differential equation for the problem. The EFM consists of obtaining an approximate solution by requiring the combination of these functions to satisfy the boundary conditions in an average sense or by multi-lint discrete matching. The pioneer in the application of the EFM to problems in solid mechanics is Quinlanll I]. The procedure has been applied by Quinlan to plane elasticity and plate bending problems while several other investigators have extended the procedure to other ~undary value problems such as shell analysis and plate vibrations[lZ151. The procedure differs from the FEM in that the governing differential equation is satisfied exactly at all points while the boundary conditions are only approximately satisfied. In the FEM the soverning differential equation is only approximated and the boundary conditions are satisfied exactly at a set of discrete points. In addition, the EFM results in a continuous functional solution which may be evaluated at any desired point to obtain numerical values. On the other hand the FEM results in the solution for displacements and stresses at discrete locations, de~nding upon the particular finite element mesh used in the mathematical model. To obtain numerical values at other locations requires some type of numerical interpolation. There are four basic types of functions which can be used to form the solution in the EFM. These are: edge functions, which are used to represent the behavior on and near the edges of the region; corner functions, which contain the terms necessary to describe the behavior near vertices in single material regions and particularly to represent the singularity in the elastic stress at a crack tip; bi-wedge functions, which describe the behavior near vertices consisting of angular wedges of two dissi~liar materials; and polar functions, which provide global stress fields and rigid body movements. The combination of these four basic types of functions is known as the function mix, A detailed discussion of the theory for the EFM, which describes the algebraic form of the functions used and the manner in which they are combined to satisfy the boundary conditions for crack problems, is given in Ref.[lO]. The discussion which follows will be limited to the application of the method. Computer uppiicatjon uj EFIU

Due to the amount of computation involved, manual application of the EFM is impractical for even the simplest problems. During the application of the procedure it is necessary to solve a large system of linear simultaneous equations in which the unknowns are the coefficients for the various functions which are included in the function mix. The number of simultaneous equations varies with the number of individual sides of the region being analyzed and with the number of harmonics being considered in the sinusoidally varying edge functions in order to satisfy the presc~~ Sunday conditions. For even very simple problems several hundred simultaneous equations may be encountered. A significant part of the investigation being discussed here consisted of developing a computer program for applying the EFM to fracture problems. The resulting program analyzes a plane elasticity problem having a polygonal boundary with straight sides, loaded with any combination of dist~buted or concentrate loads on the edges. The region may contain one or more cracks. Both single material regions and bi-material regions, with the crack extending along the interface between the two materials, may be considered. Any edge or point common to two edges, i.e. a vertex, may have prescribed displacements. The input data to the program consists of: a description of the geometry of the polygonal region being analyzed; the material properties; the loads acting on the edges; any prescribed displacements of edges or vertices; and the function mix desired. The output consists of: the EFhi Vol. 13, No. 1-D

J. F. FLEMING ef al.

50

displacements and stresses at all vertices and at a series of points along each edge; the stress intensity factors; and the stresses and displacements at any specific points or along any straight lines or circular arcs specified by the user. The FORTRAN source listing for the program, a description of it’s operation and a User’s Manual are given in Ref. 1101. The program developed in this investigation for the EFM has one distinct advantage over those used for the FEM in that an indication of the accuracy of the final solution is available. As described previously the EFM satisfies the governing differential equation exactly while approximating the boundary conditions. Since the program prints out the stresses and displacements on the boundaries, the accuracy to which the boundary conditions are satisfied can be easily seen. If poor accuracy is obtained it is possible to re-analyze the problem with a different function mix to attempt to improve the solution. From the experience gained in running the program during this investigation, it has been found that if the boundary conditions are satisfied within the desired engineering accuracy, the computed behavior at the interior of the region will also be within acceptable engineering accuracy. COMPARISON OF FEM AND EPM FOR CRACK PROBLEMS As stated previously, the purpose of this investigation was to compare the FEM and the EFM as applied to linear fracture problems. To accomplish this comparison a number of in~vidu~ crack problems were analyzed using each procedure. The four basic physical situations which were considered are shown in Fig. 4. They consisted of a long strip with: a single edge crack (SEC); a double edge crack (DEC); a center crack (CC); and an infinite series of cracks (ISC). For the majority of the individual problems a uniform tensile stress was considered to be acting on the end of the strip, although a few problems ‘were analyzed with other loadings. All computations were performed on the University of Pittsburgh DEC PDP-10

c

I_

4

DEC

b

;

cc

b

1 / j / /b b _!_

1

ISC

Fig. 4. Crack geometriesconsideredin investigation.

51

The finite element method vs the edge function method

Computer System. Symmetry was considered for each problem where possible in order to reduce the amount of computational effort required. Single material

problems

Table 1 summarizes the results of a group of single material problems which were analyzed by both the FEM and the EFM. Each of these individual problems was a long narrow strip in plane strain, corresponding to one of the physical situations shown in Fig. 4, with various ratios of the crack length to strip width, u/b, subjected to a uniform tensile stress, a, at its ends. The material was assumed to have properties similar to steel, i.e. a modulus of elasticity, E, of 30000 ksi and a Poisson’s ratio, V,of 0.30. The values which are tabulated for each problem are the nondimensionalized Mode I stress intensity factor ZU/a~(a), obtained from the analyses, and the execution time required on the DEC PDP-10. In addition, the per cent differences between the computed stress intensity factor and the reference value tabulated by Tada et al. [4] is shown. It can be seen from the few cases tabulated here that an accurate value for the stress intensity factor can be obtained by both the FEM and EFM. In fact, the percentage differences from the reference values which are given are somewhat misleading since Tada et al. acknowledge that their tabulated values are, in most cases, accurate only within approximately 0.5%. Any difference between the computed and reference values of less than 0.5% is meaningless and for all practical purposes the computed value can be considered to be exact. As stated previously, the results obtained from the FEM are highly dependent upon the mesh used in the mathematical model. For all of the problems presented in Table 1 the accuracy can be determined by comparing KI to the reference value; however, for most practical applications of the FEM to fracture problems a reference solution will not be available. The Table 1. Comparison of FEM and EFM single material problems

Crack

SEC

SEC

SEC

DEC

DEC

DEC

cc

cc

cc

Geometry h b Inches

3.0

3.0

3.0

3.0

3.0

3.0

3.0

3.0

3.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

a Method

KIlov/(a) Numerical Reference

EFM

EE 3:7303

FEM/CSQ

4.9122

0.4

0.5

3.7284

Percent difference

CPU Times (min : set)

-0.80 -0.10 0.05

6:21.00 7:03.00 4 : 32.36

- 1.94

3:51.56

0.13

10:37.81

- 1.42

3:50.19

0.46

18:59.37

- 1.42

2 :47.72

-0.18

4: 16.09

-0.77

3 :42.20

5.0096 EFM

5.0160

FEM/CSQ

7.0348

0.6

7.1367 EFM

7.16%

FEM/CSQ

1.9771

0.4

2.0056 EFh4

2.0019

FEM/CSQ

2.0451

0.5

2.0610 EFM

2.0719

FEM/CSQ

2.1648

0.6

0.53

5 : 03.24

- 0.38

3:46.03

0.91

6 : 47.50

- 0.57

3:52,46

- 0.03

5 : 29.33

- 0.91

3 : 54.35

2.1731 EFM

2.1928

FEM/CSQ

1.9544

0.4

I.%57 EFM

1.9651

FEMKSQ

2.0841

0.5

2.1034 EFM

2.1034

FEM/CSQ

2.2848

0.6

0.00

9: 18.24

- 1.05

4:00.57

0.03

6:48.93

2.3091 EFM

Note: All solutions are plane strain.

2.3098

J. F. FLEMING et al.

S?

difficulty which faces the analyst is how to determine, for these cases, whether an accurate solution has been obtained. It has been demonstrated by both Guydish[16] and Pentz[ 171 that the computed strain energy, which is available from most FEM computer programs, may be used as an indicator of solution accuracy. If the computed strain energy, for a given physical situation, is plotted vs the per cent error in the value of KI for various finite element meshes, for the same element type, it will be found that the strain energy will approach a maximum value (i.e. potential energy approaches a minimum value) as the error decreases. Therefore, if several different solutions are available for the same physical situation where the difference is due to the finite element mesh which is used, the most accurate overall solution will be that with the greatest strain energy. The value of the strain energy for any individual problem is in itself not of any significance for accuracy determination; however, by running several problems for different mesh parameters the strain energies can be used as a guide to relative solution accuracies. It should be recalled that the minimum energy theorem guarantees only an overall accurate solution, not necessarily a point by point or locally accurate solution. For example, near the crack tip the local accuracy may not rank with the overall accuracy, although, one is inclined to make a positive correlation. The same approach cannot be relied upon for different individual solutions, obtained for the same physical situation, using the EFM with various function mixes. Any solution obtained by the EFM, using a particular function mix, is a solution for some physical situation, however, the boundary conditions might not be those which are desired. Since there is no guarantee that the strain energy will vary in any definite pattern with a change in the boundary conditions, the strain energy cannot be used as an indicator of solution accuracy. The accuracy of the solution can be determined, however, by checking stresses and displacements on the edges for compliance with the desired boundary conditions. If the boundary conditions are satisfied within an acceptable tolerance the total solution should also be acceptable. Therefore, the function mix which best satisfies the boundary conditions will also give an accurate solution in the interior of the region. Bi-material

problems

Table 2 shows a summary of several of a group of bi-material problems which were analyzed by both the FEM and the EFM. There are no reference solutions available in the literature for these problems; therefore, a percentage error for the solutions cannot be given. The comparison must be made directly between the values attained for the non-dimensionahzed stress intensity factors using the two methods. Several different finite element meshes were considered to obtain the FEM solutions which are presented while several function mixes were used for the EFM solution. The best solutions were chosen on the basis of the change in the

Table 2. Comparison of FEM and EFM bi-material problems Crack

SEC

Geometry h b a (inches)

3.0

2.0

Method

Kl/ad(a)

KH/ad/(a)

Ko/&(a)

CPU times (min: set)

FEM/CST

4.9710

0.6319

5.0110

18:48.35

EFM

5.0334

0.4649

5.0548

38: 55.74

1.0

DEC

3.0

2.0

1.0

FEM/CST FEM/LST EFM

2.0148 2.0512 2.0789

0.2637 0.2251 -0.1873

2.0320 2.0635 2.0873

20:46.52 66:35.01 31:S8.21

cc

3.0

2.0

1.0

FEM/CST FEM/LST EFM

1.9943 2.0347 2.0724

0.1784 0.1346 0.0943

2.0023 2.0391 2.0745

20:55.19 70:23.36 32: 20.28

FEM/CST

I .8898

0.1933

1.8997

22 : 39.44

ISC

3.0

2.0

1.0 EFM

1.9919

1.9961

31:2264

Note: All solutions are plane strain.

-0.1301

The finite element method vs the edge function method

53

strain energy for the FEM and the computed boundary conditions for the EFM as described previously. In the single material problems there was only one stress intensity factor to be considered for Mode f loading. In bi-material problems, both Mode I and Mode II stress intensity factors are presented. It can be seen from Table 2 that the values of H, obtained for any of the reported physical situations, agree very well between the FEM and EFM solution. The value of KG however, shows a wide variance. Since the complex stress intensity factor, Ko, is a combination of both KI and KZ1 K=KI+iKII

(2)

or K = Ku eiB

(3)

where Ko =

(KY + Kuy2

(4)

it might be more appropriate to compare the values obtained for Ko for each method. The tabulated values for Ko are in good agreement, therefore, it can be concluded that essentially the same solutions (i.e. differences of l-S%) can be obtained by the two methods. As a final demonstration of the accuracy of the EFM a solution was obtained for the case of a center crack in an infinite plate which was solved exactly by Rice and Sih[I8]. Table 3 shows a comparison between the stress intensity factors for the solution of Rice and Sih, the solution obtained by the EFM and a finite element solution obtained by Lin and Marl191 using a special 17 node bi-material crack tip element. All these solutions agree within 1.4%, which verifies the accuracy of the EFM solution. Computer~ec~tio~ time An indication of the efficiency of the FEM and the EFM can be obtained by comparing the computer execution times listed in Tables 1 and 2. For the single material problems the execution times for the EFM were greater for most individual problems than required for the FEM. For the bi-material problems the reverse is true. The reason for the change becomes apparent when one considers the difference in the two types of problems. In the single material problems the crack corresponds to a line of symmetry, therefore, it was possible to use a mathematical model consisting of the material on only one side of the crack. In the bi-material problems the sym~try condition did not exist and the mathematical model included material on both sides of the crack. Therefore the total number of nodes in the mathematical model, for the FEM, and the number of sides for the ERM, approximately doubled. In the FEM the execution time increases with the number of nodes being considered while in the EFM it increases with the number of sides. For the size of problems being considered here the increase in execution time was greater for the FEM than for the EFM. No general conclusions can be reached as to which method requires the greater execution time since it is obviously dependent upon problem size. Although an attempt was made to use efficient pro~mming methods for all of the programs used in this study, a great deal of time or effort was not spent to obtain an optimum computer code. It is quite possible that the execution times for any of the programs could be considerably improved if some time were spent in code optimization, particularly in disk input and output operations.

Table 3. Center crack in infinite bi-material plate Solution Rice and Sih Lin and Mar EFM

KI

KII

I.7512 0.1283 1.7671 0.1290 1.7754 0.1294

Ko 1.7559 1.77lJ 1.7801

J. F. FLEMING et al.

54

CONCLUSIONS

The investigation has demonstrated that both the FEM and the EFM can be used to obtain accurate solutions to linear elastic fracture problems. The following conclusions can be stated: (1) Both methods require the use of a digital computer. The computer resources required by the FEM and EFM programs used in this investigation are comparable. As the size of the problem increases the computer execution time increases for both methods. (2) For either method it might be necessary to perform more than one computer analysis to obtain an accurate solution for any given physical situation. In the FEM analyses the type of elements and the size and distribution of elements near the crack tip can be varied. By comparing the computed strain energies the relative accuracy of various FEM analyses can be determined. For the EFM the accuracy of any individual analysis can be determined from the aggreement between the computed stresses and displacements on the edges and the desired boundary conditions. The quality of the solution can be varied by changing the function mix. (3) The ease with which either method can be applied is highly dependent upon the input requirements and output data for the computer program. If a mesh generator is available for the FEM programs the input requirements for the two methods are comparable. The absence of a mesh generator gives a major advantage to the EFM since considerable time and effort is required to describe a large finite element mesh manually, while it is only necessary to describe the edges of the region for the EFM. One additional advantage of the EFM is in the output data available from the computer. Since a continuous functional solution has been generated it is possible to compute stresses and displacements at any desired points in the region. To obtain the same information from a FEM analysis requires interpolation which can lead to unacceptable errors. (4) The stress intensity factor is obtained directly in the EFM since it is one of the coefficients computed in the function mix. In the FEM the stress intensity factor must be inferred by a data reduction process using the computed displacements. (5) The FEM is easily adaptable to other types of problems than those considered here, such as three dimensional solids, non-linear material properties and inelastic behavior. The extensions of the FEM to these additional problem types in restricted by the lack of the proper functions for use in the function mix. Considerable additional analytical work would be necessary to find such sets of functions before extension of the EFM could be contemplated. (6) In view of the rather widespread investment in and availability of commercial FEM programs and their applicability to non-linear and inelastic material behavior, it seems that the use of the EFM will be restricted to very special situations. The investigation has demonstrated that the EFM can be a useful tool in the analysis of linear elastic fracture problems. The EF’M computer program which has been developed can be applied to a wide variety of plane elasticity problems, Its accuracy is comparable with that obtained from any commercially available FEM computer programs. Acknowledgemenfs-The research reported here was sponsored by National Science Foundation Grant No. ENG 7621640 to the University of Pittsburgh with M. L. Williams as Principal Investigator. Acknowledgement is given to J. L. Swedlow of Carnegie Mellon University for his numerous technical suggestions and contributions and to C. C. Yates of the University of Pittsburgh for acting as Project Administrator. The CSQ and LSQ computer programs, which were used in the investigation, were provided by Thiokol Corporation, Wasatch Division, specifically for use in this project. The data reduction procedure for extracting the bi-material stress intensity factors from the results of the FEM analyses was provided by R. E. Smelser. Acknowledgement is also given to P. M. Quinlan who provided an excellent insight into the EFM through a short course he presented at Georgia Institute of Technology.

REFERENCES [I] E. Byskov, The calculation of stress intensity factors using the finite element method with cracked elements. Int. 1 Fracture 6(2) (1970). [2] W. K. Wilson, Finite element methods for elastic bodies containing cracks. Mechanics of Fracture. Noordhoff, Leaden. [3] R. D. Henshell and K. G. Shaw, Crack tip elements are unnecessary. Int. J. Numerical Methods in Engng 10, (1975). [4] H. Tada P. C. Paris and G. R. Irwin, The Stress Analysis of Cracks Handbook, Del Research Corp., Hellertown, Pennsylvania (1973). [5] S. K. Chan, I. S. Tuba and W. K. Wilson, On the finite element method in linear fracture mechanics. Engng Fracture Mech. 2 (1970). [6] G. P. Anderson V. L. Ruggles and G. S. Stibor, Use of finite element computer programs in fracture mechanics. Int. J. Fracture 7(l) (1971).

The finite element method vs the edge function method

55

[7] J. J. Guydish and J. F. Fleming, Optimization of the finite element mesh for the solution of fracture problems. Engng Fracture Mech. 10 (1). (1978). 181 _ _ M. L. Williams, The stresses around a fault or crack in dissimilar media. Bulletin of the Seismological So&v - ofI America, Vol. 49, No. 2, (1959). [9l F. Erdogan, Stress distrtbution in a non-homogeneous elastic plane with cracks. f. Appi. Me& Trons ASME 30 (1%3). [IO] J. F. Fleming, J. J. Guydis, J. R. Pentz, C. E. Runnion and J. L. Swedlow. A comparison of the finite element and edae function methods for- linear fracture analysis. Rer Rep. SE’JEC DO 77-994, School of En~nee~~, U~versity if Pittsburah. National Science Foundation Proiect ENG 74-2164tl.(1977). [Ill P. M, Q&an, The A method for Skew plates. Pmt. 4th U.S. !&r. Co&. Appt. Mech. (l%Z). WI M. J. G’Callaghan, The edge function method for linear boundary value problems with discontinuous boundary conditions. Ph.D. Dissertation, University College Cork, Ireland (1972). 1131G. V. Kelly, The A method for boundary value problems. OAR Res. Rep. (1%7). [I41 I. H. Tai and W. A. Nash, Vibrations of thin plates-a new approach. AFOSR Scientific Rep. AFOSR-TR-744789 (19731. [I51 hf. J. ‘G’CaIlaghan,W. A. Nash and P. M. Quinlan, Vibrations of thin elastic shells-a new approach. AFOSR Interim Rep. AFOSR-TR-764731 (197.5). of Civil u4 J. J. Guydish, Application of the finite element method in linear fracture mechanics. M. S. Thesis, blent En~n~~ng, Unive~ity of Pittsbu~h (1976). WI J. R. Pentz, The finite element and edge function methods in linear fracture mechanics. N. S. Thesis, Department of Civil Engineering, University of Pittsburgh (1977). WI J. R. Rice and G. C. Sih, Plane problems of cracks in dissimilar media, 1. Appl Mech. Tnvls ASME, Ser. E 32 (1%5). [I91 K. Y. Lin and J. W. Mar, Finite element analysis of stress intensity factors for cracks at a bi-material interface. Inr. J. Fracture 12(4)(1976). (Received 5 March 197%received for pnbficafion 21 June 1979)