Fracture failure prediction using a moving finite element mesh refinement scheme

FRACTURE FAILURE ~RE~ICTI~~ USING A ~~~1~~ FINITE ELEMENT MESH REFINEMENT SCHEME

(Received 7 March 1989) Abstract-A two-dimens~oaa~ fraCrure mechanics tin&e element computer program using an automated mesh refinement scheme in regions of high stress ooncentration is presented. Co&ant stress-strain triangutar elemerus are used to discretize the structural system so as to predict the initiation and propagation of fracture, yield patterns, load-displacement history and failure load under mormtanicaly increasing loads, The material stress-strain curve is idealized by selecting yield points on the nonlinear portion of the curve and joining them by linear segments. For each yield point taken, a minimum element area is specified which is gradually decreased as one moves from a lower yield point stress (or strain) to a higher yield point stress (or strait& As Iosds are i~cr~a~t~~ regions of I&h stress/strain ~n~~~ratjon

in ZLR inpur coarse En&eetement mesh am deteeed by using either the van M&es or the St. Yenant yiei& rc&ed mtii their areas reach the premtibed minimum area requirements. A f&e-noded trianguIar element is used to confine the mesh ~e~n~e~t locally. To r~~st~bute the energy of the fractured element into the unfractured media the zero modulus unload-reload method is used. The method is demonstrated for a center-cracked rectangular panel under tensile load.

criterion. Then the &snests in these fegicms are repe&f&

EWRODUCTiUN The

study of crack propagation and fracture problems is of great importance in many engineering fields. Frequently, not only the final fracture strength of a structural system is of interest to the design engineer but also the yield patterns, the crack path, the crack Iength and the ~oad~is~la~m~t history up to collapse. Various researchers have suggested finite element analysis procedures in which the structural system as a whole is analyzed, and crack ~ro~~~atio~ paths and ultimate fracture load (or collapse load) are predicted, MiIIer et ar! fl] presented the ~p~l~catjon of an incremental eIastic-plastic finite eiement solution procedure for a 2024-T3 centercracked aluminum rectangular panel. They used constant stress-strain triangular elements and a nodal release method to account for crack p~o~aga~on, However, the Unite element ~~~~t~o~s did not shovv agreement with e~~~rnent~ rest&s obtained for similar specimens. Belie and Reddy[Z] proposed the ‘zero modulus unload-reload method’ to redistribute the loads after an element fractures. Using the constant. stress-strain triangular elements and the St, Venant yield criterion, they d~rno~s~rat~ the a~~~~~tjon of their technique for pi.tension test specimen and the center-cracked panel studied earlier by Miller et al. [I], The numerical results they reported agreed very well with experimental results provided an extremely fine mesh (which was manually input in their study) was used 23

in the caiculations, Although the work reported by Belie and Reddy [2J has definite merits, the use of St. Venant or maximum principal strain theory for yielding as well as fracture, limits its use to only uniaxiai loading and to regions with very small stress gradients ~~he~ore, meta&c solids are observed to yield at values closely predicted by von Mises or maximum distortion strain energy theory. Kukreti et al. [3] developed a more general plane stress twodimensional computer program using the zero modulus unload-reload method. This program included a library of two-d~rne~s~o~a~ elements {three- and sixnoded triangular elements and four- and eight-noded isoparametric quadrilateral elements), both von Mises and St. Venant yield criteria and a mesh preprocessor to refine an input crude mesh to give the desired fine mesh before the stress analysis is started, This procedure was further extended by Kukreti et a6 @] to solve thr~d~m~~s~ona~ fracture mechanics problems. Besides presenting an a~~Iication to a three-dimensional problem, they also demonstrated that even for planeWstress fracture problems the use of three-dimensional elements gives results very close to experimental vdues with a much cruder mesh, thus requiring much fess ~m~u~tional time than that needed to conduct a two-dimensional analysis with an equal degree of accuracy. From the aforementioned studies it is apparent that to successfully model plane stress fracture probferns using two~~m~ns~onai finite elements, it is

A. R. KUKRETIand

24

essential to incorporate an automated procedure to refine the mesh to any desired degree of accuracy in the yielded and fractured regions. This will eliminate the guess work in preparing a fine mesh in such regions and a much cruder mesh in regions where stress does not vary significantly. In this paper a moving local finite element mesh refinement scheme is presented to solve plane-stress fracture problems. As in previous studies [2-4], an incremental loading technique is employed to monotonically load the structure and predict its response until failure occurs. The actual uniaxial stress-strain curve of the material is idealized by selecting points (called yield points) on the nonlinear portion of the stress-strain curve, and joining these points by linear segments. For each yield point taken, a desired minimum element area is specified. In order to prevent a rapid increase of the number of system degrees of freedom at the early stages of the analysis, the value of this area is gradually decreased as one moves from a lower yield point stress (or strain) to a higher yield point stress (or strain). As loads are incremented, regions of high stress/strain concentration in an input coarse mesh are detected by using either the von Mises or the St. Venant yield criterion. Then, the elements in those regions are repeatedly refined until their areas reach the prescribed minimum element area requirements. Five-noded triangular elements, called fringe elements, are developed to engulf the refined regions, thus confining the mesh refinement locally.

A. H. SHAMMAA

% Cl

E2

c3

64

et

Strain (a)

Actual

stress-strain

Qi

Q,-I

7

/’

curve

/’ 6,

/’

\

Linear approximation

dl-1

(b)

Linear approximation

of a segment

Fig. 1. Linear approximation of the stress-strain curve.

segment, which is given by ai+oi+I

Esi = ___. Ei+Li+I

FINITE ELEMENT MODELING USED FOR FRACTURE MECHANICS

In order to analyze nonlinear behavior using the finite element method for elastic-plastic materials, a certain number of points (called yield points) are selected on the stress-strain curve, as shown in Fig. la. Then, the stress-strain curve is linearized by a number of linear segments which connect the selected yield points (see Fig. lb). The greater the number of points selected, the closer the linearized stress-strain curve will approach the actual nonlinear stress-strain curve. The modulus of elasticity of any segment, shown in Fig. lb, can be determined from ui-Cl-1 Ei = ~ Cj- CCi _1

(1)

where oi and ei are the engineering stress and strain, respectively, at the ith yield point of the stress-strain curve. The value of Poisson’s ratio, vi, is assumed to vary in the plastic range according to the following equation given by Bert et al. [5] and Nadai [6]: vi = i - (; - v,)(E,/E,)

(2)

where v, is the elastic Poisson’s ratio, E, is the elastic modulus, and Esi is the secant modulus at the ith

Therefore, for each segment of the stress-strain diagram, the material property values are represented by Ei and vi. A three-noded planar triangular element is used to discretize the system as an assemblage of finite elements. The element formulation is based on small deformation theory and uses an incremental loading technique to monotonically load the structure so as to progress along the linearized stress-strain curve of the material. Two approaches for checking yielding of an element are considered, namely the von Mises and St. Venant yield criteria. The von Mises yield criterion is considered both in terms of principal strains and principal stresses, in order to move along the plastic region of the stress-strain curve for elastic-perfectly plastic (e.g. steel) and plastic (e.g. aluminum) materials, respectively. To start the finite element analysis, all elements are assumed to have elastic material properties and a small load increment is applied to the structural system and is increased linearly until the effective stress or strain, computed by using either the von Mises or the St. Venant yield criterion, in at least one element reaches the lower yield point of the next segment of the linearized stress-strain curve. This element will yield and its material properties are

Fracture failure prediction using a moving FE mesh modified using eqns (2) and (3). A second cycle is performed in which loads are incremented linearly until, as before, the effective stress or strain in at least one element reaches its next segment’s lower yield point. Consequently, this element’s material properties are modified. This process is repeated in cycles until at least the effective stress or strain in one element reaches the last yield point on the linearized stress-strain curve (representing the fracture point). Then, the element is said to have fractured and the initial fracture load is the sum of all load increments accrued in all previous cycles. Stresses, strain and displacements are proportioned in the same manner. When an element fractures, the contribution of that element to the structural system must be eliminated (i.e. the element must be made ‘inert’) by redistributing the load carried by that element to the remaining unfractured system. This is accomplished by using the procedure proposed by Belie and Reddy [2], as follows. 1. The fractured element is made ‘inert’ by changing its modulus of elasticity to zero and Poisson’s ratio to 0.5. 2. The whole structure is unloaded following the elastic response of the unfractured specimen while retaining each element’s progress along the linearized stress-strain curve. 3. The system is reloaded back to the original fracture load level with the fractured element removed. The unloading-reloading procedure is repeated until another element fractures and is made ‘inert’. If all elements surrounding a node fracture, then to continue the analysis the degrees of freedom of that node are condensed out. If, during the reloading cycle, conducted to distribute the strain energy of the fractured element into the unfractured specimen, another element fractures before reaching the previous fracture load level, then unstable crack growth is said to occur and the procedure is terminated. The total load at this stage is defined as the ultimate fracture load of the specimen.

25

A

point

t

5 0-l

Strain

Fig. 2. Minimum area specified for each yield point taken on the linearized stress-strain curve.

based on either the von Mises or the St. Venant yield criterion, should be equal to or less than the value specified for AMIN (1). If not, then this element should be repeatedly refined until in the refined mesh the area of the yielded element satisfies this condition. In order to prevent a rapid increase in the number of system degrees of freedom at the early stages of the analysis, the value of this area is gradually decreased as one moves from a lower yield point stress (or strain) to a higher yield point stress (or strain), as shown in Fig. 2. One refinement method that can be adopted is to divide the element into two by bisecting its largest angle. At the same time, the element joining the element bisected must also be divided into two, as

MOVING MESH REFINEMENT SCHEME USED

Once an element is found to be highly stressed, it should be subdivided (refined) into smaller triangular elements. These smaller elements must be compatible with the existing triangular elements. This means that there can be no sagging between the interfacing edges of the new elements and the neighboring elements. In this paper a moving local finite element mesh refinement scheme is presented, in which for each yield point selected on the linearized stress-strain material curve, there is a prescribed minimum area requirement specified, as shown in Fig. 2. For example, the value specified for AMIN (1) means that the area of the element in which either the effective stress or effective strain becomes equal to 0, or c, , respectively,

l

Original

0 Added

nodes nodes

Fig. 3. Refinement method by bisecting an angle.

A. R. KUKRETI and A. H.

26

0

Original nodes

0 Added nodes

~/~/ Y

Y

(a) Somph starting mesh

(b)

(cl

Mesh after element 0

Mesh after element

@

is refined

is refined

Fig. 4. Mesh re~ne~ent using four-noded and five-neded fringe elements.

shown in Fig, 3. Referring to this figure, it can be seen that by this procedure a new node 5 is introduced alongside 2-3. The advantage of this method is that it is easy to automate. Furthermore, no new special element is introduced as a result of the refinement process to preserve compatibility. However, one grave disadvantage of this method is that after a number of refinements, the mesh is bound to end up having slender trianguiar elements (i.e. triangles with very acute angles), which may lead to large computational errors as shown by Carey and Humphrey [7] and Zienkiewicz [8]. Another method for element refinement, which is adopted in this study, is based on the principles suggested by Carey [3]. This method is explained in detail in the remainder of this section. It

SHAMMAA

overcomes the disadvantages of the aforementioned method and adds to its advantages. Consider a mesh as shown in Fig. 4a, which has 14 elements. Say element 1 is the most highly stressed (or strained) element and has to be refined. Therefore, element 1 is sub~vided into four congruent subelements numbered as 1, 15, 16 and 17 (see Fig. 4b) by introducing three new midside nodes. Furthermore, to continue the refinement of the highly stressed (or strained) region, ail elements that have a common edge with element 1 are also subdivided (i.e. elements 2, 5 and 6 in this example) using the same procedure as used for element 1. Thus, the refined mesh will have 26 elements, as shown in Fig. 4b. In this figure, the nodes denoted by a solid circle are the nodes that existed before the current mesh re~nement was started, and those denoted by an open circle are the nodes added to obtain the current refined mesh. Note that in Fig. 4b the elements 7-9 and 1l-13 are four-noded elements and these engulf the locally refined region. Such elements will be called ‘fringe elements’. Now, say in the next cycle, element 3 has to be refined. So, now, element 3 and the element that has a common edge with it, i.e. element 4, are subdivided using the same procedure as used in the previous refinement cycle. The refined mesh so obtained is shown in Fig. 4~. In this figure, note that elements 9 and Ii-13 are four-noded fringe elements, whereas elements 7 and 8 are five-noded fringe elements. Hence, for the same problem, there are three types of triangular elements present in the mesh, namely three-noded, four-noded and five-noded triangular elements. To obtain the system stiffness matrix for this mesh, three different types of element stiffness matrices would have to be derived and coded. Moreover, the locations of these three types of elements, with respect to other elements, have to be known at all times in the computer program, especially for the purposes of further re~nement. For example, to refine an element its type, location, orientation with respect to neighboring elements and the types of all elements that have a common node with it have to be known. Since at any refinement state three types of elements may exist in different positions and orientations, it may require a very complex automated computer logic to keep the bookkeeping in an ever growing mesh. Therefore, for programming ease, only two types of eiement are used, namely three-noded triangular parent elements and five-noded fringe triangular elements, Whenever during retinement a four-noded element is created, a fifth node will be added on one of the other edges of the element to make it a five-noded element. To a~ompIish this, the following rule is followed: all other element edges that come into any of the three nodes of a refined element must have midside nodes. Following this rule, the nodal configuration that would be obtained after refining element 1 of Fig. 4a would be as shown in Fig. Sa instead of Fig. 4b. Thus, all fringe elements (4 and 7-14) are five-noded trian-

Fracture failure prediction using a moving FE mesh Y

I

27

3

ILLL 4

5

2

x

(a)

Nodal and element after

element

0

configuration (a)

is refined

Fringe element nodes

Y

1

x

I

(b) (b)

Nodal and element

configuration

after

is refined

element

@

Element numbers and nodal coordinates of the congruent element obtained by adding node 6 ( l original nodes, 0 added nodes and 0 fictitious nodes)

Fig. 5. Mesh refinement using five-noded fringe elements only.

gular elements. Similarly, if element 3 has to be further refined, then the mesh obtained will be as shown in Fig. 5b instead of Fig. 4c. It should be noted that in the five-noded fringe element, the displacement variation along an edge which contains the midside node would be bilinear. Since this element has five nodes, the resulting element stiffness matrix will be 10 x 10 in size. This is formulated next. Consider a fringe element with corner nodes 1, 2 and 3 and midside nodes 4 and 5 added along edges 2-3 and 3-1, respectively, during refinement, as shown in Fig. 6a. The stiffness matrix of the fringe element can be derived by assembling the stiffness matrices of the four congruent triangular elements obtained by adding a sixth node on the third remaining side, i.e. node 6 along edge 1-2, as shown in Fig. 6b. Note that the solid nodes are original nodes before refinement, the empty nodes 4 and 5 physically exist due to the refinement and node 6 is a fictitious node for the purpose of driving the fringe element stiffness matrix. The assembling is done, based on the constraint that the displacement vector at midside node 6 is equal to the average of the displacement vectors of the corner nodes 1 and 2 (see Fig. 6c), i.e. (4)

x I-u,

(c)

Element nodal displacement numbers

Fig. 6. Fringe element configuration.

The 10 x 10 stiffness matrix, [K],, of the fringe element can be obtained by first assembling the element stiffness matrices of the four congruent triangular elements shown in Fig. 7, thus giving a 12 x 12 element stiffness matrix, [K], which corresponds to the six-noded triangular element (nodes denoted as 1, 2, 3,4, 5 and 6). Then, node 6 is condensed out using

A. R. KUKRETIand A. H. (

STApT

SHAMMAA

)

,Iw( and setting-up

t Stoge II : Nodal optimization-i_m/ t Stoge III : Assembly and solution

t Stage JX : Determination of yielded elements I Stage P : Checking for element refinement 1 Stage XI : Mesh plotting and modifications

Fig. 7. General program organization.

the constraint given by eqn (4). This can be done by using the following transformation, [T]: r

4

where

1000000000 0100000000

,

0010000000 0001000000 0000100000

(5)

PI=

0000010000

(12 x 10)

0000001000 0000000100 0000000010 0000000001 ,Of0000000 1 I d-0-000000 1 2

(6)

Fracture failure prediction using a moving FE mesh shown in Fig. 6. These strain components puted from

In short form, eqn (5) can be written as

ACID

Ed,

PI

CE]n = Pin

(7)

(12 x 1) = (12 x 10) (10 x 1)

are com-

WI,

(11)

where, for the congruent element n, [B]” is the strain-displacement transformation matrix and {q’), is the element nodal displacement vector with node 6 being condensed out using eqn (4). In view of the nodal coordinates defined in Fig. 6b and the nodal displa~ments defined in Fig. 6c, the nodal displacement vectors for the congruent elements are defined as

where (q} is the nodal displacement vector of the six-noded triangular element and {q}, is the nodal displacement vector of the fringe element (nodes are 1,2,3,4 and 5), as shown in Fig. 6. Thus the stiffness matrix, [K],, of the fringe element is obtained from

iTI

P-IT wl Klr (10x10)=(10X12)(12x12)(12x10)

29

(@

where superscript T denotes transpose of a matrix. The strain components vector, (t-j,, of the fringe element can be obtained from

14lr

{Elf _ Plr

(9)

(3 x 1) - (3 X 10) (10 X 1)

where [BJ, is the strain-displacement transformation matrix of the fringe element. To obtain the elements of this matrix, it is assumed that the strain components of the fringe element are equal to the average of the strains of the four congruent elements

@2- iY3) [WI

2

0

0

(&-X2)

($x3---2)

c_Y*-fY3)

(3 x 6) = z

0

(3 x 6) = G PI2

2

Pa

2

:

$3

j

cvz-Y3)

(3 x 6) = z

r

1

0

0

-i lx 3

0

x,

X2

-Y2

$ys

0

g_Y2+y3)

0

-y2

0

(x3 -x2)

0

x3

0

x2

cy2 -Y3)

-x3

y3

x2

-y2

4Y2

0

0

0

-5 ‘X 2

0

-i ‘X 2

--ix2

;y,

--d--X*

fY2

(see Fig. 6), i.e. @)n @}I ,’ i (3 x 1) 4,=, (3 x 1)

(10)

where {L),, for n = 1, 2, 3 and 4 represents the strain component vectors of the four congruent elements

-Y2

--f(X2+Xj) 0

-4(x2+x3>

iy*

1

0

y3

@3-x2)

PI4 2 (3 x 8j=&

-Y2

-x2)

0

0

0

-5X3

ul"2v2u4u4~~~~)

Wd)

and the strain-displacement transformation matrices for the congruent elements are obtained as

$52+y,f

fCY*--Y3)

&3-x2)

iq’fr =(", (8 x 1)

fY3

0

fCY2-Yd

(124

-672

-Y3)

0 --(x3 -x2)

x2

(134

0

x2 -y2

1

Wb)

(13c)

0

-Y3

0

-(X3-X2)

0

X3

-cy*--Y3)

x3

-Y,

I

(134

where Al-A4 are the areas of the congruent elements l-4, respectively. In eqns (12a)-(12d), ( > denote a row matrix. Substituting eqn (11) into eqn (10) gives (14)

A. R. KUKRETIand A. H. SHAMMAA

30

3081 Model D computer. The computer program locates and refines proper regions of a starting coarse mesh, locates crack initiation and propagation path, predicts yielding patterns and provides a loaddeflection history up to failure. The program analyzes two-~mensional structures using three-noded and five-noded triangular (fringe) elements based on small displacement theory. Only structures having isotropic materials can be analyzed using this program. The program developed is divided into five main stages. These stages are explained in the following subsections and are shown in Fig. 7.

Any congruent element nodal displacement vector, (q‘)#, when node 6 is condensed out, can be related to the fringe element nodal displacement vector, (q],, by the following transformation: WI” Wf” &71/ (m x l)=(m x lO)(lO x 1)

(15)

where m equals 6, 6, 6 and 8 for congruent elements 1,2, 3 and 4, respectively, and (qjf has been defined earlier in eqns (5) and (7). The elements of the matrices [T’], will contain zeros and ones and can be easily written. Substituting eqn (15) into eqn (14) gives

Stage I: input and setting-up A coarse mesh generated either by hand or by an automated mesh preprocessor is required to start the analysis. This stage involves the following steps. 1. Basic info~ation that controls the program operations is first read. This information includes the following input data: total number of nodes (NUMNP), total number of elements (NUMEL) and a flag IFAIL indicating which yield criterion is to be used (IFAIL = 1, 2 and 3 for St. Venant yield criterion expressed in terms of principal strains, von Mises yield criterion expressed in terms of principal

Now, comparing eqns (9) and (16) gives the required strain-displacement matrix for the fringe element as (17) which, upon substitution yields:

of eqns (13a-d) and [T’],,

(Yz+Yd

(X,-2%) (%--5X2)

0

(2Y,--Y3)

--(x2+.%)

0

-(x3+x*) (Yz-f-Yd

Y3

0

0

-xj

0

Y3

x2

-x3

-Y2

0

x* -Yz

-&2+y3)

0 (2x,+x3)

0

Px,+-%)

-(2Y,+Yd

.

j (18)

Finally, the stress components vector, (o>J, of the fringe element is determined from

rq

PI hllf @If (3 x 1) = (3 x 3) (3 x IO) (10 X 1) where [D] is the standard elasticity matrix of the material comprising the fringe element, which for an isotropic material plane stress problem is given by

11 0

Iv

p]_EI -v2

v 1 O

0

0

(20)

y

where E = modulus of elasticity and v = Poisson’s ratio of the material. DESCRIPTIONOF THE COMPUTER

PROGRAM

DEVELOPED

A computer program was developed using the procedures described in the previous two sections and was executed on the University of Oklahoma’s IBM

stresses and von Mises yield criterion expressed in terms of principal strains, respectively). 2. The x and Y coordinates of each node are then input and stored, respectively, in arrays X(J) and Y(f), where J = 1,2,. . . , NUMNP. In addition, element connectivity (i.e. the system nodes belonging to the element) is provided and stored in an array called NI(I, J), where I = 1,2,. . . , NUMEL and NEM (=number of element nodes). J=l,2,..., Also, since the structural system may consist of elements involving different materials, a material code is provided for each element and stored in an array MAT(I), where I = 1,2, . . . , NUMEL. 3. The constrained and loaded nodal numbers and the value of each nodal load is input. 4. Using the input node numbers, the nodal displacement numbers are computed and stored in an array calied IDF(J), where J = 1,2,. . . ,NDOF (= total number of system degrees of freedom). 5. Subroutine PROP is now called. This subroutine sets up the material property arrays. These arrays are needed to input the material data for each segment of the stress-strain curve joining the yield points. Therefore, the following arrays are set up. (a) ET(f) = an array containing the value of

Fracture failure prediction using a moving FE mesh modulus of elasticity for each segment of the linearized stress-strain curve, where I = 1,2, . . . , NOPC (=number of yield points taken on the stress-strain curve). (b) XNU{I) = an array containing the value of Poisson’s ratio at each yield point, where Z = 1,2, . . . , NOPC. (c) STRAIN(Z) = an array containing strain value at each yield point, where Z = 1,2, . . . , NOPC. (d) STRESS(Z) = an array containing stress value at each yield point, where Z = 1,2, . . . , NOPC. (e) AMIN = array containing the minimum area to be achieved at each yield point, where Z = 1,2,. . . , NOPC. The (NOPC + I)th element of these arrays is set equal to 999.0, which denotes fracture. 6. The system load vector, B(Z), where Z=l,2,. . . , NODF, is formulated by assembling the nodal loads in the appropriate rows using the loaded nodal numbers. The values of the elements in the system load vector, (B), for the nodal directions which are not loaded are kept as zeroes. Stage II: nodal optimization

This stage consists of subroutines JSET, SETUP and OPTNUM. These subroutines are explained in the subsequent paragraphs. Subroutine JSET: this subroutine finds the current ~nimum half bandwid~, IBD. Subroutine SETUP: this subroutine first formulates the array {JMEM}, in which the .Zth row contains the number of element edges coming into the Jth system node. In the present version of the program, it is fixed that the number of element edges coming into any node should not exceed 12. The subroutine then constructs a nodal connectivity array called {MEMJT) in which 12 rows are reserved for each system node. So the size of this array is (12 * NUMNP). For the Jth system node, the rows reserved in the array (MEMJT} start from the (12*&Z- l)+ I)th row and end at the (12*J)th row. How many element edges or how many nodes are connected to node 3 is known from the array (JMEM}. Say this is equal to K. The node numbers of such nodes connected to node J are stored in rows (12 * (J - 1) - 1) to (12 *J) of the array {MEMJT). If K < 12, then in the (12 * (.Z - 1))th row to the (12 * (.Z - 1) + K)th row, the numbers of the K nodes associated with the Jth system nodes are entered and the remaining rows from the (12 * (J - 1)th row to the (12 * .Z)th row will have zero entries. Subroutines OPTNUM: this subroutine optimizes the nodal numbering so that the minimum bandwidth is obtained. The node nurn~~ng optimi~tion is based on the procedure suggested by Collins [lo]. Stage III: assembly and solution

This stage consists of the following subroutines: CYCLE, CYCLEl, SOLVE, SOLVEI, and TRIAN. In this stage, first the elements of the upper half

31

bandwidth of the system stiffness matrix are constructed and stored in a rectangular array, [A]. Then, the system stiffness equilibrium equations are solved for the unknown nodal ~spla~ments using the Gauss elimination method. Subroutines CYCLE and SOLVE are used when the number of degrees of freedom is less than or equal to 2000. Then, the computer core (for IBM 3081) is sufficient to perform an in-core solution in a reasonable time and cost. When the number of degrees of freedom is greater than 2000, subroutines CYCLE1 and SOLVE1 are used and an out-of-core solution scheme is formulated so that large problems can also be analyzed. The subroutines of this stage are described in the subsequent paragraphs. Subroutine CYCLE involves the following steps. 1. The system stiffness matrix (rectangular), [A], is initialized. 2. Element loop starts. 3. For the first cycle, the number of element nodal displacements, MDF, and the element nodal displacement number vector, MID(Z) (where Z = 1,2,. . . , MDF), are first stored in a direct access file number (say 1). Then, depending on the yield criterion to be used, either the element strain-displacement transformation matrix, [ZJB], or the element stressdisplacement transformation matrix, [DB], and the element stiffness matrix, [EK], are formulated for each element and stored as the next two records, respectively, in the direct access file number 1. The formulation of matrix [BB] (or [DB] and [EK]) is accomplished by calling element subroutine TRIAN (described later). 4. The element stiffness matrix, [EK], is assembled at the proper locations of the system stiffness matrix, [A 1, using the element nodal displacement numbers stored in the array (MID). 5. Steps 2-4 are repeated until the element loop ends. 6. The diagonal elements of the system stiffness matrix are checked to find whether all the elements su~ounding a node have fractured. If any diagonal element of matrix [A] is equal to zero, then the corresponding nodal degree of freedom is constrained to continue the analysis. 7. If at certain system nodes, the nodal displacements are prescribed as a result of given boundary conditions, the appropriate rows and columns of the system stiffness matrix, [A], and the system load vector, {Z?), are modified so as to describe the particular boundary condition. Subroutine CYCLE1 involves the same steps as previously presented for subroutine CYCLE, except that the system stiffness matrix is subdivided into a number of blocks. Each block is a square matrix of the size (NBD x NBD), where NBD = 3 * (IDIFF + 1). The bandedness

and symmetry

properties

(21) of the

32

A. R. KUKRETIand A. H. SHAMMAA

system stiffness matrix are utilized by recognizing that only (2 * NBD) equations are needed at one time in the computer core to perform the Gauss elimination procedure for one block (i.e. convert it to an upper triangular matrix if a complete square system stiffness matrix was utilized). Therefore, at any time, space for two blocks is created in-core and the stiffness matrix is formulated by sequentially assembling two blocks at a time. The first block of the system stiffness matrix is formulated by assembling the element stiffness matrices of those elements whose nodal displacement numbers fall between 1 and NBD. Some elements should share degree of freedom numbers belonging both to the first and second blocks. After the first block has been formulated, it is stored in a sequential file (say file number 2) and the second block is shifted up in its place, thus reserving the remaining bottom space for the third block. The procedure is repeated and now the element stiffness matrices of those elements whose nodal displacement numbers fall between (NBD + 1) and (2 * NBD) are assembled. Then, the second block is stored in sequential file number 2, following the previous record already stored. In this manner, the system stiffness matrix is formulated in blocks. Thus, the Nth block of the system stiffness matrix is built up by assembling the element stiffness matrices of those elements whose nodal displacement numbers fall between ((N - 1) * NBD + 1) and (N * NBD). The number of blocks, NBL, required is computed from NBL = 1 + (NDOF - l)/NBD

(22)

where NDOF is the total number of degrees of freedom. Subroutine SOLVE uses the Gauss elimination procedure to solve the system stiffness equilibrium equation for nodal displacements when an in-core solution is performed. The conventional Gauss elimination method is modified to account for the fact that only the upper half bandwidth elements of the system stiffness matrix are stored in a rectangular array. Subroutine SOLVE1 is similar to subroutine SOLVE, except it is called when the system stiffness matrix is stored out-of-core. In this process, two blocks of system stiffness matrix are read in-core and Gauss elimination method is applied. Then the first block will be triangularized and is then stored out-ofcore in sequential file number 2. The second block, which has been partially triangularized, is shifted up and the third block is read into the in-core of the second block. This process is repeated until all the blocks have been triangularized. Nodal displacements are computed by the back substitution process by bringing the last two blocks in-core to compute the

Load at first yield point =

nodal displacements of the last block. Then, to compute the nodal displacements of the block before last, it is shifted down in place of the last block and the next highest block is brought in-core. This process is repeated until all the unknown nodal displacements are computed. Subroutine TRIAN is called from subroutines CYCLE and CYCLE1 for the first cycle, and from subroutine MODFY (explained later) for subsequent cycles. This subroutine generates and stores the required element matrices for the three-noded and the five-noded fringe triangular elements. It should be recalled that before subroutine TRIAN is called for an element, its total number of nodal displacements (MDF) and nodal displacement numbers (stored in {MID}) are computed and stored as the first and second records, respectively, out of a set of four direct access records generated for each element. Then subroutine TRIAN is called, which performs the following two steps for each element. 1. If a strain yield criterion is selected (i.e. IFAIL = 1 or 3 for St. Venant or von Mises yield criterion using effective strain, respectively), the strain-displacement transformation matrix, [BB], is computed and stored as the third record of the direct access file. But, if a stress yield criterion is selected (i.e. IFAIL = 2 for von Mises criterion using effective stress), then the stress-displacement transformation matrix, [DB], is computed and stored instead of [BB]. 2. The stiffness matrix, [EK], is computed and stored as the last record of the direct access file. Four such records are generated for each element and stored in direct access file number 1. Stage IV: determination of yielded and fractured elements

In this stage the yield criterion is applied to predict element yielding or fracture. As mentioned before, a flag IFAIL, which can take a value equal to 1, 2 or 3, determines which yield criterion is to be used. If IFAIL = 1, then the St. Venant yield criterion expressed in terms of principal strains is used by calling subroutine VENANT. If IFAIL = 2, then the von Mises yield criterion expressed in terms of principal stresses is used by calling subroutine VONMIS. If IFAIL = 3, then the von Mises yield criterion expressed in terms of principal strains is used by calling subroutine VONMIS. Subroutines VENANT and VONMIS are similar and so only subroutine VONMIS is described in the subsequent paragraph. Subroutine VONMIS involves the following steps. 1, For each element principal stresses or strain are evaluated. Depending whether IFAIL = 2 or 3, effective stress or strain, respectively, is computed using the von Mises yield criterion. 2. In the first cycle, the first yield load is calculated from

Stress (or strain) at first yield point of the stress-strain Maximum element effective stress (or strain)

curve .

(23)

Fracture failure~r~~tioa usinga moving FE mesh For later cycles, a load inurement is applied to predict further yielding. The Ioad increment is evaluated as follows, (a) The cumulative maximum total effective stress (or strain) of each element is calculated at the present cycle. fb) Then it is sub~a~ted from the stress (or strain) value of the next yield point selected on the material stress-strain curve. (c) The rn~~irn~ of this offered when divided by the ~creme~tai effective stress (or strain) due to unit toad gives the load i~~ement. 3. The total effective stress (or strain’) of each element is compared with the stress (or strain) value of the next yield point selected on the material stress-strain curve to check if the element has yielded. Xfthe total effective stress (or strain) becomes equal to the next yield point stress for strain) value, then it is said to have yielded. All yielded element numbers are stored in an array called IEFQV), where N=l,2,..., NEF (= number of yielded and fractured elements at the current cycle), 4. Displa~me~~ at loaded nodes and the total appli~ load at the current cycle are computed and printed out. 5. If an element fr~tu~s
All elements whose numbers are stored irr the array (IEF) are candidates for refinement. This stage uses subrautines AREACK, REFINE, ADJ and ADE. These subroutines are explained in the subsequent para~aphs. Subroutine AREACK checks the areas of the element num~rs stored in the array (IEF) against the ~Rimurn areas read in the array ~A~~~ provided in subroutine PROP. Now, according to the stress or strain level of an element, whose number is stored in the array (IEF] (say the element stress level is somewhere in the second seventy, the element area is checked with the ~~espondi~g specified minimum area (shown as AMIN (2) in Ftg, 2). If the element area is greater than the corresponding specieed rn~~murn area, then this element number is stored in an array called IEFR flu), where K-1,2,. . . 3NEFR (= total number of elements to be refined~, and is a candidate for re~nement. The ~finement is done by calling subroutine REFINE,

33

which is dewibed in the next paragraph, The aforementioned check is performed on all efements whose numbers are stored in the array (IEF]. Subroutine REFINE refines the elements whose numbers are stored in the array (IEFR). However, there are situatjons where the content of an eles ment is not possible (these s~t~ationsdo not exist in the star&g mesh). For example, consider the starting mesh shown in Fig. 8a. An array to store ~nfo~a~~n on how many times any element has been refined is generated and called LEVEL (M), where M=1,2,..., BETEL. So for the starting mesh the array (LEVEL) is i~i~ali~ed to zeroes denoting that no element has yet been refined. Say element 2 is refined, as shown in Fig. 8b. Then, LEVEL (1) = LEVEL (2) = - . . = LEVEL (16) = 1, which indicates that elements l-16 have all been refined once. Now, say that element 16 has to be refined next. It can be seen from Fig, 8b that adding three midside nodes on the edges of element 16 creates incumpatibility between element 17 and 18,each of which may end up with four nodes on one of its edges. This situation is detected by using the array ~LE~L~ in the following manner” {a) ff LEVEL (16) - LEVEL (17) = 0, and if LEVEL (lo] - LEVEL (18) = 0, then no ~~ompa~bility problem exists in the re~nement procedure adopted for eIement 16, (b) But, if LEVEL (14) -LEVEL (17) =J, or if LEVEL (16) - LEVEL (18) = N, where J and 1v are any scalar numbers, then this indicates that in~m~atibil~t~ exists between element 16 and its neighbo~n~ efements. In order to over~me this problems both f and N have to be made equal to zero by refining element, 17 J number of times and element f 8 N number of times. For the mesh shown in Fig. 8b, both f and N are equal to 1. As a result, both elements 17 and 18 have to be refined once more, as shown in Fig. 8~. So finally, element 16 is refined as shown in Fig. 8d. It can be seen from this figure that this re~~erne~t process overcomes the problem of “refinement incompatibil~ty’~with the nei~hbo~ng elements. This procedure of first refining adjoining elements (i.e, eIements 17 and 18 in this example) and then the element concerned (i.e* element I6 in the e~~ple~ will be called ‘swappings. Subrout~~ REFINE ~~voivesthe follo~~g steps. 1. Cheek if any element whose number is stored in the array (EEFR) has a problem of tenement incompa~bi~t~~as explained in the previous ~ragraph. If incompatibility exists, then remove the element number from the respective row in the array {ZEFRf and insert the proper neighbo~~g element numbers, In other words, perform swapping. Repeat this step, if needed. 2. Loeate element numbers of all elements which have two nodes common with any efement whose number is stored in the array {IEFR). All such elements found are called tioining’ elements and

A. R.

34

(a) A sample

and A. H.

refining

refining

element

element

@

SHAMMAA

in both arrays {IEFR} and {JELEM} are the elements which are to be subdivided, check refinement incompatibility with the neighboring elements for the elements whose numbers are stored in the array {JELEM}. 4. Find the fringe elements for the elements whose numbers are stored in the array {IEFR} and store their numbers in an array called JEFELEM (N), , NFJE (= total number of fringe whereN=1,2,... joining elements). 5. Call subroutine ADJ. This subroutine adds three new midside nodes on the elements whose numbers are stored in the arrays {IEFR} and {JELEM}, as shown in Fig. 9. 6. Add two new midside nodes (if not already added from the previous step) on the proper sides of the elements whose numbers are stored in the array {JFELEM}. 7. Next call subroutine ADE. This subroutine adjusts the node connectivity indices for the elements which have been refined. For example, for element 4 of Fig. lOa, the node connectivity is changed from l-2-3 to 18-19-20. As a result, there are three new elements (elements 20, 21 and 22 of Fig. lob) which have to be added to the current mesh. For this example, this is accomplished by the following steps. (4 Node connectivity is established for elements 20, 21 and 22 as 2&10-l& 18-3-19 and 19-5-20, respectively, as shown in Fig. lob, and then stored in an array called {NI}. (b) The new elements are assigned the same material code as well as the number of refinements (stored in array {LEVEL}) of the parent element 4.

mesh

(b) The mesh after

Cc) The mesh after and @

KUKRETI

@

, @

IO

5 (a)

An element

to be refined

3

(d) The mesh after to

Fig.

8.

refining

element

@

and

A

@

19

@

0

An illustration of refinement incompatibility.

the numbers of these elements are stored in an array called JELEM (L), where L = 1,2, . . . , NJE (= number of joining elements). 3. Since the elements whose numbers are stored

5L--d--AO (b)

New mldsides

18

20

nodes added

Fig. 9. Addition of new nodes in subroutine ADJ.

35

Fracture failure prediction using a moving FE mesh

(a) Element

@

records generated for each element. This step is performed by calling subroutine EDOF. 2. For the new added and yielded elements, the strain-displacement transformation matrix or the stress-displacement transformation matrix, depending on the type and form of the yield criterion selected, is computed and stored as the third record of the direct access file. 3. For the new added and yielded elements, the element stiffness matrix is computed and stored as the last record on the direct access file. Four such records are generated for each element and stored in direct access file number 1. Steps 2 and 3 are executed by calling subroutine TRIAN. Stages V and VI are repeated in cycles until crack instability is detected, as explained in the previous section.

nodal connectivity

NUMERWAL EXAMPLE: ANALYSIS OF A CKNTER-CRACKED PANEL

5

20

IO

(b) New element odded

Fig. 10. New elements added in subroutine AK%

(c) For the new elements, the total stress and strain components are kept the same as for the parent element 4, up to the current cycle. (d) The progress of the new elements along the material stress-strain curve is assumed to be the same as for the parent element. 8. Change the node ~nn~tivity of the fringe elements, whose numbers are stored in the array (JFELEM), from those pertaining to a three-noded element to one pertaining to the five-noded triangular element. 9. Adjust the total number of nodes, NUMNP, and the total number of elements, NUMEL, to their new values after the current refinement. Stage V is repeated until all elements whose numbers are stored in the array {IEF} reach their respective prescribed minimum area requirement.

A common specimen tested in the laboratory to predict fracture and related stress intensity factors is the center-cracked panel. Miller er al. [l] have presented experimental results for a center-cracked panel specimen made of 2024-T3 alumina, the configuration of which is shown in Fig. Il. The pre-crack length in the specimen is 5 in. Belie and Reddy [2] and Kukreti et al. [4] have also compared finite element analysis results for this center-cracked panel with the experimental results reported by Miller et al. [l]. In the present study, this problem is solved using the ~mputer program developed and the results are

Pre-crack length

30in

In this stage the mesh can be plotted and yielded elements marked using some standard plotting software package. In this study, the Basic Plotting Package [11] available at the University of Oklahoma was used. In addition, displacement~train transformation or displacement-stress transformation and stiffness matrices for each yielded element and for each new added element are updated in subroutine MODFY, which involves the following steps. 1. For the new added elements, the total number of element nodal displacements, MDF, and the element nodal displa~ment num~rs, MID (J), MDF, are computed and stored as the J-l,&..., first and second records of a set of four direct access

Thickness m0. I in

k

Ibin.

-I

Fig. 1I. The cracked panel configuration.

36

A.R. KUKRETI and A. H. SHAMMAA

compared with the aforementioned experimental and finite element results. The cracked panel shown in Fig. 11 is symmetrical about an axis passing through the crack length and also about an axis passing through the midcrack length and perpendicular to it. Therefore, only a quarter of the panel needs to be considered in the finite element model by choosing appropriate boundary conditions along the nodes on both axes of symmetry. The starting crude finite element mesh used for this quarter panel is shown in Fig. 12. Five yield points are selected on the stress-strain curve for 2024-T3 aluminum, and the material property values for these yield points and the different linearized regions are shown in Table 1, As discussed in the preceding section, a minimum element area is specified for each yield point selected on the stress-strain curve. The finite element mesh is automatically refined such that all elements which have reached a certain yield point stress o,r strain have the minimum area specified. To investigate the effect of gradually decreasing the magnitude of the area of the yielded elements as their state of effective stress or strain approaches the fracture point, six sets of minimum element areas for the various yield points are selected. The values of these areas are shown in Table 2. This results in six finite element meshes to be run using the computer program developed. These are tabulated as Run 1, Run 2, . . . , Run 6 in Table 2. A value of 999.0 for the minimum element area is adopted if a mesh is not to be refined any further. Thus, in Run 1 the input crude mesh shown in Fig. 12 is not refined. Note that the last column of Table 2 contains the actual minimum element area achieved at the crack tip. Only the von Mises yield criterion is tried for this example because Kukreti et al. [4] had concluded for this example problem that the von Mises yield criterion produces better results than the St, Venant yield criterion. The yield loads and fracture loads obtained from the computer program presented in this paper for the center-cracked panel example are presented in Table 3. The last column of this table contains the yield loads and fracture load obtained by using a three-dimensional solid finite element, as reported by Kukreti et al. [4] for their finest mesh. For this mesh, the minimum efement volume was equal to 0.002 in’

‘-

I-

Crack length Fig. 12. The finite element model.

(which corresponds to a minimum element area 0.02 in’). Comparing yield load results in Run 1 and Run 2, it can be seen that by decreasing the minimum element areas of Run 1 by 75% (which gives Run 2), the value of the first yield load predicted is 23% less. Also, in Run 2 the fracture toad predicted is 10% less than that predicted by Run 1. Similarly, compa~ng the fracture load predicted by Run 1 with that of Run 6, it is seen thai by decreasing the minimum element area by 99.9% at the fifth yield point (fracture point on the material stress-strain cruve), a 19% lower fracture load is predicted in Run 6 than that in Run 1. It can also be seen from Tables 2 and 3 that comparing the results of the finest mesh reported by Kukreti et al. [4] with those obtained from Run 3 (for which the minimum element area is about twice that of [4]), the first, second, third and fourth yield loads predicted by the computer program presented in this

Table 1. Material properties at yield points se&ted on the stress-strain curve Yield point no.

Strain (in/in)

Stress (ksi)

0

0.0

:

0.02175 0.0045

46.923 56.923

3 4 5 6

0.06375 0.10650 0.19500 999.0

65.384 69.615 13.462 999.0

0.0

x

-i

Yield region number?

Poisson’s ratio

Modulus of elasticity (ksi)

Elastic 1 2 3 4 Fracture

0.313 0.429 0.474 0.486 0.491 0.500

10427.0 579.7 201.5 98.9 43.5 0.0

f The ith yield region connects linearty the ith yield point to the (i + l)th yield point on the stress-strain curve.

Fracture failure prediction using a moving FE mesh

37

Table 2. Minimum element area required at the yield points selected on the stress-strain Yield point no. 0

1 2 3 4 5 6

curve Actual area at crack

Input minimum areas (in*) Run 1

Run 2

Run 3

Run 4

Run 5

Run 6

tiu (in21

999.0 999.0 999.0 999.0 999.0 999.0 999.0

0.781 0.700 0.700 0.700 0.700 0.700 999.0

0.781 0.700 0.100 0.100 0.100 0.100 999.0

0.781 0.700 0.100 0.030 0.030 0.030 999.0

0.781 0.700 0.100 0.030 0.010 0.010 999.0

0.781 0.700 0.100 0.030 0.010 0.003 999.0

0.781t 0.19501 0.04888 0.0122]/ 0.0031~ 0.0008tt 0.0008tt

t For all runs. $ For all runs except Run 1. 5 For all runs except Runs 1 and 2. 11 For Runs 4, 5 and 6. q For Runs 5 and 6. tt For Run 6.

paper are 26%, lo%, 0.9% and 0.6%, respectively, higher than those predicted by them. However, the fracture load predicted by this study is 1.9% lower than that predicted by them. The fracture loads predicted by Belie and Reddy [2] (using a two-dimensional triangular three-noded linear finite element), Kukreti et al. [4] (using a threedimensional solid finite element) and the computer program presented in this paper are listed in Table 4. In addition, the experimental fracture load obtained by Miller et al. [1] is 36.0 kips. The last column of Table 4 shows the precentage difference between the finite element analysis predicted fracture load reported by various researchers and the experimental value. It can be seen that for the mesh which produces fracture results closest to the experimental value, the computer program presented in this paper predicts a fracture load which is 0.8% (lower) away from the experimental value. On the other hand, the best result reported by Belie and Reddy [2] and Kukreti et al. [4] for the fracture load is 20.8% (higher) and 1.4% (higher), respectively, away from the experimental value. The computer program developed was found to be computationally more efficient than that reported by Kukreti et al. [4]. For the different minimum element areas, the Central Processor Unit (CPU) computer time required using the computer programs developed by Kukreti et al. [4] and the present study are shown in Fig. 13. From the figure, it can be seen that

for the range of minimum element areas tried by the two computer programs, the computer program developed in this study is more efficient than that developed by Kukreti et al. [4]. Thus by using an automated mesh refinement procedure, the computer program can be made to predict results more accurately and efficiently. The various yield and fracture loads versus the minimum element area in the crack vicinity are plotted for the six meshes tried in Fig. 14. The fracture loads found by Belie and Reddy [2], Kukreti et al. [4] and Miller et al. [1] are also overlaid in this figure. For the computer program developed, it can be seen from this figure that as the minimum element area is decreased, the different yield and fracture loads converge to lower values. Also, as the area is decreased, the fracture load predicted by the program converges from a higher value towards the experimental value. Finally, for very fine element areas, the predicted fracture load ultimately crosses the experimental value to a value which is 6.7% lower than the experimental one. It is possible that this may be attributed to the accumulation of truncation and round-off errors. A value of the minimum element area that should correspond to the exact fracture load can be interpolated between Runs 2 and 3. Such a value is computed to be equal to 0.078 in2. It can also be seen from Fig. 14 that the fracture load curve predicted by this study corresponds to a more flexible system than the systems which correspond to the

Table 3. Predicted vield loads at different load levelst Predicted yield loads (kips) Yield point

1st 2nd 3rd 4th Fracture

Run 1

Run 2

Run 3

Run 4

Run 5

Run 6

141

19.5 31.2 35.9 38.3 41.5

15.1 29.0 32.8 35.0 37.2

15.1 28.5 31.9 33.8 35.7

15.1 28.5 31.1 32.7 34.6

15.1 28.5 31.1 32.4 33.6

15.1 28.5 31.1 32.4 33.6

12.0 25.9 31.6 33.6 36.4

t For the six cases selected, the minimum element area of yielded elements at the various yield loads (including fracture) is different from the one reported in [4].

38

A. R. KUKRETIand A. H. SHAMMAA Table 4. Fracture loads at different minimum element areas predicted by different researchers and the present study Fracture load (kips) Present study Experimental [I]

Run no.

Load

% Difference?

41.5

+86.1 + 18.1 +15.3

2

37.2

+80.6 + 8.3 -I- 3.3

3

35.7

+ 55.6 + 1.4 - 0.8

4

34.6

+47.2 - 3.9

48.5

5

33.9

+34.7 - 5.8

43.5

6

33.6

f20.8 - 6.7

r21

[4]

67.0

42.5

65.0

39.0

56.0

36.5

36.0 53.0

-

t Negative percentage difference means below experimental value and positive percentage difference means above experimental value. For each fracture load, the first value listed is that obtained in [2], the second is from [4] and the third is from the present study.

fracture curves predicted by Belie and Reddy [2] and Kukreti et al. [4]. An approximate value of the fracture load can be computed using the following equation [12]: P=

Pmhn- Pnh,

hn- ha

(24)

where P = P,,, = P, = h, = h, =

the extrapolated fracture load, the fracture load predicted in Run the fracture load predicted in Run minimum element area in Run m, the minimum element area in Run

m, n, and n.

25 -

m-

0

Two-dlmenrlonal

0

Three-dimenaonal

(present study1 (Ref t41)

CPU computer time versus minimum element element area.

0

1st yield

l

2nd yield

.

4th yield

A

Fracture

.

Fracture 143 *

o ----

3rd yteld Experimental fracture Cl1

Fracture t21

14. Yield and fracture loads versus minimum area.

Fracture failure prediction using a moving FE mesh Table 5 shows the results obtained by applying eqn (24) taking m = i for Run i (where i = 1,. . . ,6). It can be seen from this table that the precentage difference between predicted and extrapolated fracture load decreases as the mesh is made finer. Note that the extrapolated load, P,would always be lower than the two predicted loads, P,,, and P,,, involved in its evaluation (i.e. P,, < P,,,). When P is equal to P,, then eqn (24) is no longer used. The yielded loads predicted versus the total number of elements in the finite element mesh for the different cases studied are presented in Fig. 15. It can be seen from this figure that the number of elements does not necessarily increase as the minimum area of the elements is decreased gradually as the fracture point on the material stress-strain curve is ap proached. This is because smaller elements result in a more precise prediction of stresses and strains in the panel, thus producing a better localized refined mesh which does not extend needlessly in regions where it is not desired. The minimum element area achieved by the fractured elements was decreased further than that reported in Table 2. For most of the cases, it was found that the results did not vary appreciably (maximum 5% variation for predicted yield loads and fracture load). However, in cases when a very small minimum element area was specified at the fracture point, it was found that the predicted fracture load suddenly increased to a very high value in comparison to the yield load predicted at the previous yield point. Table 6 shows the minimum area requirements and the results of such a run using the von Mises yield criterion. Comparing the results in Table 6 with the results of Run 6 in Table 3, it can be seen that the solution of the run shown in Table 6 diverged only at at the fifth yield point (i.e. the fracture point). The cause of such a sudden divergent trend can be attributed to computer truncation and round-off errors due to solving a set of very large numbers of linear simultaneous equations [S, 121. As the minimum element area at any yield point is signi~~tIy decreased, the system degrees of freedom of the resulting refined finite element mesh would increase. Table 7 lists the number of the resulting finite element mesh system degrees of freedom at each yield point for the six cases obtained by performing the six runs shown in Table 2. The number of system degrees of freedom for the case presented in Table 6 are also presented in Table 7. It can be seen that the number

39

0

1st yield

l

2ndyield

0

3rd ymld

9

4th ymld

1

FWJCtUW

Fig. IS. Yield and fracture loads versus total number of elements in the finite element mesh.

of system degrees of freedom for the case reported in Table 6 is much larger (2086% larger than Run 1 mesh and 43% larger than Run 4 mesh) than for the other six causes studied. This gives large solution errors.

Table 5. Extrapolated

fracture loads from the computer kU”S

Run no.

Predicted fracture load (kips)

Extrapolated fracture load (kips)

% Difference

2 3 4 5 6

37.2 35.7 34.6 33.9 33.6

35.8 35.3 34.5 33.9 33.6

3.8 1.1 0.3 0.0 0.0

Table 6. A run with diverging fracture load Yield point

Input minimum area (in2)

Actual area at crack tip (inz)

Yield load (kips)

I

0.700 0.200 0.010 0.003 0.0007

0.1950 0.1950 oBO31 0.0008 0.0002

15.1 29.0 30.9 31.9 88.1

2 3 4 5

Table 7. Number of degrees of freedom at different yield levels for all the runs Yield point tI0.

1 2 3 4 5

Number of degrees of freedom Run 1

Run 2

112 112 112 112 112

140 220 266 292 362

Run3 140 612 738 788 1112

Run 4

Run 5

140 612 902 1106 1706

140 612 902 1144 1686

Run6 140 612 902 1144 1634

Run shown in Table 6 140 432 1268 1796 2448

A. R. KUKRETIand A. H. SHAMMM

40 CONCLUSIONS

From the numerical problem presented in this paper, the following conclusions can be made. 1. The automated mesh refinement results in faster convergence of the solution than a manual mesh refinement in which guess work has to be used. 2. The automated mesh refinement improves the prediction accuracy of the results obtained. As the minimum element area requirements are decreased, the solution improvement converges towards the true value. However, if extremely small area requirements are imposed for fractured elements, the solution may diverge due to computer truncation and round-off errors resulting by solving a very large system of linear simultaneous equations. 3. For problems that can be modeled as twodimensional problems, the two-dimensional computer program developed predicts results which are comparable with those predicted with the threedimensional computer program developed by Kukreti et al. [4]. Furthermore, the two-dimensional computer program developed requires less computer CPU time than that required by the three-dimensional computer program.

2. 3.

4. 5.

advances in computerized aerospace structural analysis and design. Comput. Struct. 7, 315-326 (1977). R. G. Belie and J. N. Reddy, Direct prediction of fracture for two-dimensional plane stress structures. Compur. Sfrucr. 11, 49-53 (1980). A. R. Kukreti, A. S. Khan and A. Iranmanesh, Fracture prediction in elasto-plastic problems using finite element method. Conference, Application of Fracture Mechanics to Materials and Structures, Frieburg, West Germany (June 1983). A. R. Kukreti, A. S. Khan and A. Kumar, A threedimensional finite element program to predict ultimate fracture load. Comput. Srrucf. 32, 1325-1345 (1989). C. W. Bert, E. J. Mills and W. S. Hyler, Effects of variation in Poisson’s ratio on plastic tensile instability. J. Basic Engng, Trans. ASME 89D, 35-39

(1967). 6. A. Nadai,

Theory of Flow and Fracture

of Solids,

2nd Bdn. Vol. 1. McGraw-Hill. New York (1950). 7. G. F. Carey and D. L. Humphrey, Mesh refinement and interactive solution methods for finite element computations. Int. J. Numer. Meth. Engng 17, 1717-1734 (1981). 8. 0. C. Zienkiewicz, The Finite Element Method, 3rd Bdn. McGraw-Hill, New York (1977). 9. G. F. Carey, A mesh-refinement scheme for finite element computations. Comput. Meth. appl. Mech. Engng 7, 93-105 (1976).

10. R. J. Collins, Bandwidth reduction by automatic renumbering. Int. J. Numer. Meth. Engng 6, 345-356 (1973).

REFERENCES 1. R. E. Miller Jr, B. F. Backman, H. B. Hansteen, C. M. Lewis, R. A. Samuel and S. R. Varanasi, Recent

11. R. J. Collins, Basic Plotting Guide. University of Oklahoma, Norman, OK (1980). 12. R. D. Cook, Concepts and Applications of Finite Elemeni Analysis, 2nd Bdn. John Wiley, New York (1981).