AN INFINITE FINITE ELEMENT by RONALD FREDERICK DNGLESS B.A.Sc. (1971) The U n i v e r s i t y of B r i t i s h Columbia A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE i n the Department of CIVIL ENGINEERING We accept t h i s t h e s i s as conforming to the r e q u i r e d standard The U n i v e r s i t y of B r i t i s h Columbia September 1973 In presenting t h i s thesis in p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It i s understood that copying or publication of t h i s thesis f o r f i n a n c i a l gain s h a l l not be allowed without my written permission. R. P. Ungless Department of C i v i l Engineering The University of B r i t i s h Columbia Vancouver 8, Canada September, 1973. AS I N F I N I T E F I N I T E ELEMENT ABSTRACT T h i s t h e s i s i s c o n c e r n e d w i t h d e v e l o p i n g a f i n i t e e l e m e n t model o f i n f i n i t e s i z e t o f a c i l i t a t e t h e s t r e s s a n a l y s i s o f t h r e e d i m e n s i o n a l , s e m i - i n f i n i t e b o d i e s g o v e r n e d by t i m e i n d e p e n d e n t l i n e a r e q u a t i o n s . An e l e m e n t o f t r i a n g u l a r p l a n f o r m e x t e n d i n g t o i n f i n i t y i n one d i r e c t i o n i s d e v i s e d . T h e r e a r e t h r e e nodes p e r e l e m e n t , e a c h node h a v i n g t h r e e d i s p l a c e m e n t d e g r e e s o f f r e e d o m . The i n f i n i t e e l e m e n t i s u s e d i n c o n j u n c t i o n w i t h r e g u l a r f i n i t e e l e m e n t s t o r e p r e s e n t t h e s t i f f n e s s o f t h e s e m i - i n f i n i t e s o l i d w h i c h has p r e v i o u s l y been assumed t o be z e r o o r i n f i n i t e i n r e g u l a r f i n i t e e l e m e n t m o d e l s . The i n f i n i t e e l e m e n t r e l i e v e s t h e c o m p u t a t i o n a l p r o b l e m c a u s e d by l a r g e n umbers o f e l e m e n t s w h i c h h a s l i m i t e d t h e use o f t h e f i n i t e e l e m e n t method i n t h r e e d i m e n s i o n a l h a l f s p a c e p r o b l e m s . A l a r g e r e d u c t i o n i n t h e number o f d e g r e e s o f f r e e d o m i s p o s s i b l e w i t h t h e u s e o f t h e i n f i n i t e e l e m e n t b e c a u s e t h e a r t i f i c i a l b o u n d a r y o f p r e v i o u s m o d e l s i s e l i m i n a t e d . The a c c u r a c y o f t h e e l e m e n t i s t e s t e d on two e x a m p l e s whose e x a c t s o l u t i o n s a r e known. E o t h i n v o l v e a s e m i - i n f i n i t e s o l i d , one l o a d e d w i t h a s u r f a c e i i p e r p e n d i c u l a r p o i n t l o a d , and the other with a s u r f a c e p a r a l l e l p o i n t l o a d . The r e s u l t s compare f a v o u r a b l y with the theory of e l a s t i c i t y s o l u t i o n s . i i i TABLE OF CONTENTS P a g e ABSTRACT i i TABLE OF CONTENTS i v LIST OF FIGURES v l ACKNOWLEDGEMENTS v i i CHAPTER 1. INTRODUCTION 1 1.1 Background 1 1.2 Purpose and Scope 2 2. ELEMENT GEOMETRY 4 2.1 Topology. 4 2.2 Transformation of Axes 9 2.3 Transformation of Global Side Vectors to Local Side Vectors 13 3. AREA COORDINATES » 15 4. ELEMENT DISPLACEMENT FUNCTION 20 5. ELEMENT STIFFNESS MATRIX 26 5.1 Formulation of Element S t i f f n e s s Matrix 26 5.2 Anal y t i c Integration of the S t i f f n e s s C o e f f i c i e n t s 28 i v CHAPTER Page 6 . SEMI-INFINITE NUMERICAL INTEGRATION 38 6.1 Goodwin-Moran Method 38 6 . 2 A n a l y t i c a l v e r s u s N u m e r i c a l I n t e g r a t i o n 39 7 . NUMERICAL EXAMPLES 47 7 .1 Example 1 - S e m i - I n f i n i t e S o l i d Loaded by a S u r f a c e P o i n t Load P e r p e n d i c u l a r to t h e S u r f a c e 49 7 . 2 Example 2 - S e m i - I n f i n i t e S o l i d Loaded by a S u r f a c e P o i n t Load P a r a l l e l t o t h e S u r f a c e 5 5 8. DISCUSSION 60 BIBLIOGRAPHY 62 APPENDIX A. ELEMENT STIFFNESS MATRIX FOR LOCAL COORDINATE SYSTEM 63 B. GOODWIN-MORAN METHOD OF SEMI-INFINITE INTEGRATION 66 C. LISTING OF FORTRAN SUBROUTINE OF THE GOODWIN-MORAN METHOD 70 D. GENERALIZED LOAD VECTOR 72 v LIST OF FIGURES F i g u r e Page 2.1 A Three Dimensional I i i f i n i t e F i n i t e Element. . . . 5 2.2 P l a n View o f the Base Plane 6 2.3 L o c a l and G l o b a l Coordinates f o r an I n f i n i t e F i n i t e Element 10 3.1 D e f i n i t i o n o f Area Coordinates 16 4.1 Adjacent L o c a l Coordinate Systems 24 7.1 S e m i - I n f i n i t e Solid. Loaded by a S u r f a c e P o i n t Load P e r p e n d i c u l a r to the S u r f a c e 50 7.2 F i n i t e Element Mesh o f H e m i s p h e r i c a l Bowl w i t h One Quarter Symmetry 53 7.3 Displacement of H e m i s p h e r i c a l Boundary under a V e r t i c a l U n i t Load (X-Z Plane) 54 7.4 S e m i - I n f i n i t e S o l i d Loaded by a S u r f a c e P o i n t Load P a r a l l e l to the S u r f a c e 55 7.5 Displacement o f H e m i s p h e r i c a l Boundary Loaded P a r a l l e l to S u r f a c e (X-Z Plane) 58 7.6 Displacement of H e m i s p h e r i c a l Boundary Loaded P a r a l l e l to S u r f a c e (X-Y Plane) 5& v i ACKNOWLEDGEMENTS I wish to express my a p p r e c i a t i o n f o r the guidance and adv i c e given by my a d v i s o r . Dr. D. L. Anderson, throughout the p r e p a r a t i o n of t h i s t h e s i s . Thanks a l s o are due to Ian M i l l e r and B r i a n Badomske f o r many h e l p f u l d i s c u s s i o n s . F i n a l l y , I would l i k e to thank the R a t i o n a l Research C o u n c i l of Canada f o r t h e i r f i n a n c i a l support. September, 1973 Vancouver, B r i t i s h Columbia Vl l 1 CHAPTER 1 IIIHODUCTION }.± Background Much use of the f i n i t e element method has been made in s o i l and rock mechanics problems dealing with i n f i n i t e halfspaces. I t has been normal practice to replace the i n f i n i t e halfspace under study by a f i n i t e body with an a r t i f i c i a l boundary at some ar b i t r a r y distance from a point of interest (e.g., point of loading, location of surface structure). Since a f i n i t e number of elements i s a necessity, the reason for this boundary i s quite apparent. Because of computer storage l i m i t a t i o n s , three dimensional f i n i t e element analyses cannot afford many elements and therefore, approximating with a fixed or free boundary near the point of i n t e r e s t may be i n gross error. The -number of normal f i n i t e elements required to achieve a certain accuracy may be quite large. This w i l l r e s u l t in very large numbers of simultaneous equations, which may place a severe l i m i t a t i o n on ) 2 the usefulness of the method for p r a c t i c a l problems. Also the band width of the r e s u l t i n g equation system becomes large leading to large computer storage requirements. A f i n i t e element of i n f i n i t e s i z e , with which this thesis concerns i t s e l f , could minimize some of these problems. 1.2 Purpose and Scope The purpose of th i s thesis i s to construct a model with which the stress analysis of three dimensional, semi-i n f i n i t e bodies may be f a c i l a t e d . A f i n i t e element model i n which the element i s of i n f i n i t e size i s proposed to replace the present a r t i f i c i a l boundary conditions. A reduction of the degrees of freedom and band width of an i n f i n i t e halfspace problem may be possible with the use of the i n f i n i t e f i n i t e element. The standard method requires c a r e f u l placement of an a r t i f i c i a l boundary where stresses or displacements can be assumed negli g i b l e while the i n f i n i t e element method makes no such requirement. As a r e s u l t , fewer elements are needed leading to fewer degrees of freedom and a lesser band width than the standard method. Attention w i l l be confined to the consideration of l i n e a r e l a s t o s t a t i c problems characterized by time independent l i n e a r equations. This thesis considers an element of 3 triangular plan form extending to i n f i n i t y i n one d i r e c t i o n . Each node of t h i s t r i a n g l e has three displacement degrees of freedom and these nodes form the usual boundary. application of the element developed i s i l l u s t r a t e d by a f i n i t e element solution of the Boussinesg problem using the i n f i n i t e elements. The r e s u l t s are compared to exact a n a l y t i c a l solutions of the problem. A computer program was developed for the f i n i t e element analysis employing the i n f i n i t e element. 4 CHAPTER 2 ELEMENT GEOMETRY 2 . J To£olocjy_ The element model, shown i n F i g . 2 . 1 , i s referred to as an " i n f i n i t e f i n i t e element." This paradoxical term has been coined because the element extends to i n f i n i t y i n the d i r e c t i o n normal to the f i n i t e plane tr i a n g l e formed by the nodes. Fig. 2 . 1 shows the geometry of the element and the two right hand cartesian coordinate systems - global( X,Y,Z ) and l o c a l (x,y,z). The o r i g i n of the l o c a l coordinate system i s shown as being at node I and the x axis p a r a l l e l to the I - 3 side. This choice of axis permits integration over the area of the t r i a n g l e without dividing the in t e g r a l into two parts. The element has three nodes located at the corners. These nodes a l l l i e in the l o c a l z=0 plane. This plane w i l l be referred to as the "base plane". The l o c a l z axis i s directed into the element. 5 Fig. 2.1 A Three Dimensional I n f i n i t e F i n i t e Element 6 F i g . 2.2 Plan View of the Base Plane Three generalized displacements (u, v, w) are assigned to each corner node of the element shown in F i g . 2.2, r e s u l t i n g i n nine degrees of freedom per element. The in-plane 7 displacements u and v act in the x and y di r e c t i o n respectively and the out of plane displacement w acts in the z d i r e c t i o n . The element edges or "rays" st a r t from the nodes and t r a v e l off to i n f i n i t y in the z d i r e c t i o n . The edges can be represented by vectors V_ ( L = , , 3 ) with origins at the nodes. The vectors are shown on F i g . 2 . 1 and can be expressed as since the vector i s straight. In the above eguation the Qk 1 b % are constants and the I ,J ,TT vectors are unit vectors in the l o c a l (x,y,z) coordinate system. The superscript denotes a guantity expressed in the l o c a l coordinate system. If (*-_;y- ) are the nodal coordinates of the element in the z=0 flane and i f the directions of the edges of the element are V. = Q-L + b~j +• k then on any plane given by z=constant the nodal coordinates are y . = y ° + b x J «• J 1 1 ( 2 . 2 ) I n i t i a l l y the "rays" of each element w i l l be specified i n terms of global coordinates by the vectors 8 V * = A t I - B,J + Q K Ci-1,3) (2.3) where are unit vectors i n the global ( X,Y,Z ) system. The superscript <y denotes a quantity expressed in the global coordinate system. In Section 2.3 the transformation of the global side vectors ( V ) to l o c a l side vectors ( V ) i s carried out. The constants A[,rB;,C- may form part of the element input data or may be automatically generated. One method of doing this i s by adding the unit vectors normal to the base planes of a l l elements common to a node. The resultant vector becomes the "ray" originating at that node. The l a t t e r approach has been adopted in a computer program written for i n f i n i t e f i n i t e element problems. Some r e s t r i c t i o n s on the direction of the rays must be observed i f the method i s to produce r e a l i s t i c r e s u l t s . An element which has a l l three rays perpendicular to the base plane ( i . e . , CX^ t b. - O ) w i l l produce erroneous re s u l t s due to the chosen displacement function (see Chapter U). The s t i f f n e s s of the exact solution goes to zero i n t h i s case yielding i n f i n i t e displacements whereas the f i n i t e element s t i f f n e s s i s f i n i t e . At the other extreme, an element in which any one of i t s rays l i e s i n the same plane as the base plane, or almost so (i . e . , , )o •-*•<») w i l l experience numerical d i f f i c u l t i e s and 9 w i l l give poor r e s u l t s . A t h i r d r e s t r i c t i o n requires that the cross sectional area increase with increasing z. I f this i s not s a t i s f i e d , a zero cross se c t i o n a l area at some value of i w i l l l i k e l y occur producing numerical d i f f i c u l t i e s . Also rays must not be allowed to cross as t h i s also leads to a zero cross sectional area. 2.. 2 Transformation of Axes To transform the s t i f f n e s s matrix of an element expressed i n the l o c a l element coordinate system to a global system, and then assemble into the structure s t i f f n e s s matrix, requires a transformation matrix between the l o c a l coordinate system (x,y,z) and the global coordinate system ( X,Y,Z ). The base plane and the two coordinate systems of an i n f i n i t e f i n i t e element are shown in Fig. 2.3. There are several choices for directions of the l o c a l coordinates. For t h i s thesis the l o c a l coordinate system s h a l l be chosen as having the x axis directed along the I k side of the t r i a n g l e as shown i n Fig. 2.3. 1 0 Fig. 2.3 Local and Global Coordinates for an I n f i n i t e F i n i t e Element — fr The vector defines the ik side of the triangle which i n global coordinates i s 11 r v T.-Y, (2.4) The d i r e c t i o n cosines are found by dividing the components of t h i s vector by i t s length giving a vector of unit length. r \ 7\ x X J x... Y k(. 7 . (2.5) with 2 + Y.1 + I* k<. The l o c a l z d i r e c t i o n w i l l be normal to the "base plane". The vector product of two sides of the triangle defines t h i s normal d i r e c t i o n , thus 12 r V = V. x V. - •< Y.Z.. X..Z J L k l X« Y.. Y7, JL K t j i Kc ( 2 . 6 ) the magnitude of Vz- i s equal to twice the area of the t r i a n g l e . L-The d i r e c t i o n cosines of the z axis are obtained simply by dividing the components of V z by i t s length 2.A . This gives the unit vector f r k = * Y > - 2A Y Z - Y- Z . X k . Y -• - X . . Y K I J L J l k i (2.7) The y axis i s perpendicular to the x and z axes, and so i s found by the cross product of vectors i n the x and z d i r e c t i o n s . If vectors of unit length are taken in each of these directions as defined by Eqs. 2.5 and 2.7, we have 13 r J r A k X L - (2.8) A transformation matrix can be written containing the di r e c t i o n cosines in the global directions of the L iJ , K vectors of the l o c a l system. This matrix w i l l be denoted by T and i s given by T = \ x \ Y \ Z V V V \* \t \ z : J : (2.9) The element s t i f f n e s s matrix in global coordinates i s derived from the l o c a l element s t i f f n e s s matrix by the following r e l a t i o n [4] 9 - r T i * (2-10) k'= T k T 2.3__Transformatign_gf_Glcbal To complete the positioning of an element in space, the element edges or rays must be described. Their d i r e c t i o n 14 can be described by a vector whose form i s given in Eg. 2.1 in terms of l o c a l coordinates. I n i t i a l l y the rays w i l l be described by global coordinates thereby necessitating a transformation to l o c a l coordinates. If the ray vectors i n the l o c a l and global systems are written as local-: = OL T 4 lo[ J + C . ' ^ (.= 1,3 global: V.? = A; I + B t J + C . K L = l ' 3 then the unknown c o e f f i c i e n t s CL^ , i b i r C _ must be determined from the known c o e f f i c i e n t s A : . B i . C : . The scalar product of the global ray vector with a unit, l o c a l vector produces the component of the ray vector in the unit l o c a l vector d i r e c t i o n . Thus a- - V*- t - A ; \ „ + B ; Ji, Y + C L A > z c;.V . T - IT - < \ . \ x * B ; A , Y * C ; A l Z then following the form of Eg. 2.1; the constants are redefined as b ^ = b i . 7 C l ' ,2.11, c . = c / / , < , I 15 CHAPTER 3 MM_COORDINATES In the development of a f i n i t e element one of the basic necessities i s to establish a r e l a t i o n between the continuous displacements u,'v,w and the nodal displacements. To aid in thi s matter, i t i s convenient to introduce a special set of normalized coordinates for a triangle ca l l e d area or natural coordinates. Let A, . A , and r\ be the areas of the three subtriangles subtended by P and the corners (see F i g . 3.1), the index designating the opposite corner number. The area coordinates 0- of the point P are defined by (3.1) where A = / \ 1 + A 3 L + A 3 = area of the element. Thus the area coordinates of the nodes of the element are: 16 node I node 2. node 3 also, a useful i d e n t i t y i s proved using the d e f i n i t i o n p i - f . - o ( 1 , 0 , 0 ) p t - i , f , = o ( 0 , 1 , 0 ) = 0 , ?3= I ( 0 , 0 , 1 ) = ^ w n l c h c a n D e e a s i l y of P • and A . 3 (0,0,1) (l,0.0) F i g . 3 .1 D e f i n i t i o n of Area C o o r d i n a t e s . 1 7 The relationship between the x-y coordinates and the area coordinates i s l i n e a r and can be shown to be given by the expression _ P-x y x , x 2 x 3 Pi (3.2) Inverting this equation and using the identity + ^3 = ' ?, I ' ? x > - I 2A (3.3) J - y where 2 A X. y. y* 18 - ( v ^ / ) ( y ; ^ 3 i ) + ( < ^ v ) ( y > b z ) = A + B z + Cr2-(3.4) Rewriting Eg. 3;3 i n index form gives r Pi ( v v) (y;* i ° k 0 - ( x ; + a k Z ) ( y ; + b j Z ) (3.5) 19 Or J ' A + B ^ + C z 1 where i n Eg. 3.5 the indices take on the permutation J K 1 Z 3 Z 3 I 3 I 1. In a plane p a r a l l e l to the x,y plane, the area coordinates are functions of the x,y nodal coordinates of the tr i a n g l e and as we proceed i n the z di r e c t i o n these same nodal coordinates become functions of z. Thus, the area coordinates become functions of z. Area coordinates of a tr i a n g l e i n a plane p a r a l l e l to top surface tr i a n g l e or "base plane" 20 CHAPTER 4 ELEMENT DISPLACEMENT FUNCTION The role of the displacement function of an element i s to establish a r e l a t i o n between the nodal displacements and the displacements of a l l points i n the element. There are two i n t e r r e l a t e d factors which influence the se l e c t i o n of a displacement function. F i r s t , the p a r t i c u l a r displacement degrees of freedom that describe the model must be selected. In Chapter 2 these were chosen as the displacements of the nodes i n the three coordinate directions. Second, to ensure that the numerical r e s u l t s approach the correct solution, the function should s a t i s f y certain conditions such as providing compatible displacements between adjacent elements and the a b i l i t y to describe r i g i d body displacements of the element. A fixed boundary at i n f i n i t y w i l l prohibit a r i g i d body displacement of the whole element. However a r i g i d body displacement of each cross section of the element (z=constant) i s a possible condition which the displacement function must be 21 able to describe. This describes a state in which only shear stresses i n the plane z=constant e x i s t . The shape function which relates the continuous displacements u,v,w to the nine nodal displacements takes the form U = FKy) Pfz) s = A s (4.D where r \ U U = \ v > w V J c o n l i n u o u s dis pla ce me n l s in ocal coorolinale d i r e c t i o n s and 22 s = < > = nodal d i s p l a c e m e n t s in loca c o o r d i n a t e d i r e c t i o n s The variation of the continuous displacements u,v,w in the x-y plane i s given by F(X y) • T h e simplest function for F ( x ^ ) , which w i l l s a t i s f y compatibility of displacements when the base planes of adjacent elements l i e i n the same plane, i s a lin e a r r e l a t i o n between nodal and continuous displacements, Since a li n e a r r e l a t i o n between the (x,y) coordinates and the area coordinates e x i s t s , a lin e a r function of area coordinates can be used for interpolation over the triangle lying in the x-y plane ( i . e . , z=constant). Therefore, F (x,y0 w i l l be expressed as follows ?, o o f , o O f , o o Rx,y) - O f, o O o o p. ° (<*.2) o O f, O 0 f t 0 ° p. In order for the displacements to vanish at i n f i n i t y 23 u,v and w must approach zero as z tends to i n f i n i t y , must be a decay function that modifies the displacements u,v,w by decreasing them monotonically with increasing z u n t i l they vanish at i n f i n i t y . Also must be unity at z=0. A function for POO. which s a t i s f i e s the above requirements and yet remains simple enough to carry out analytic integration in some cases, i s P ( z ) = 0 + Z ) " (4.3) and w i l l be used throughout the thesis. The chosen displacement function jf(*,_y) w i l l s a t i s f y continuity between elements whose base triangles l i e i n the same plane. However this i s , i n general, not the case and d i s c o n t i n u i t i e s at the interface can occur. To i l l u s t r a t e how th i s discontinuity may occur, an example i n Fig. 4.1 i s given. « 24 Fig. 4.1 Adjacent l o c a l Coordinate Systems A discontinuity arises because a point on the boundary between adjacent elements would have di f f e r e n t values of the l o c a l coordinate z, as shown i n Fig. 4.1. This i n turn w i l l produce d i f f e r e n t values of displacement for the same point on the boundary. This discontinuity can be minimized i f the angle of the ray vector to each base plane i s approximately the same. However the f i n i t e element derived i n t h i s thesis i s based on a non-conforming displacement function. In summary a very elementary form of displacement function has been chosen to represent the displacements of the element. This yields a simple f i n i t e element but one which, i t w i l l be shown, gives good r e s u l t s for the deformations of semi-i n f i n i t e s o l i d s . 2 6 CHAPTER 1 5 ELEMENT STIFFNESS MATRIX 5^ 1^ Formulation of Elejment_Stiffness The element s t i f f n e s s matrix i s developed for a homogeneous, lin e a r e l a s t i c material. C l a s s i c a l e l a s t i c i t y theory i s used in the development with only the lin e a r portions of the strain-displacement r e l a t i o n s being retained, thereby l i m i t i n g the analysis to small displacements. The strain-displacement equations are given by: (5.1) 27 or ( : L U where 6 x , £^, £z are the normal s t r a i n compoments and X X 0 are the components of shear s t r a i n . If the assumed form of the displacement function given in Eg. A.3 i s inserted into the strain-displacement equation (Eq. 5.1), the strains may be expressed in matrix notation i n terms of the unknown nodal displacements as ( : L u = L A s = B s (5.2) It can be shown [ 1 ] that the element s t i f f n e s s matrix in l o c a l coordinates, K , i s given by B T E B dV _ _ _ (5.3) V where the integration i s carried out over the volume of the element. E i s the constitutive matrix r e l a t i n g stress to s t r a i n , and i s given by 28 - (lrVXl"2v) I v/|-v V/-v O O O I v/-v O O O 5 Y M O ,-2v 20-v) O O l-2v 2 ( l -v ) O o O 1-2V (5.4) E. and V are Y,oungfs modulus and Poisson's r a t i o where respec t i v e l y . Terms of the l o c a l s t i f f n e s s matrix k given i c f u l l in appendix A. are 5_.2 ftPgl£tic_ Integration,,,of_the_Stiffness_Coefficients The s t i f f n e s s c o e f f i c i e n t s as given by Eg. 5.3 and shown i n Appendix A must be integrated over the volume of the element. The f i r s t two integrations, carried out in the x-y plane, cover the cross-sectional area of the element for any value of z. The th i r d integration i n the z-direction completes the volume i n t e g r a l by ranging from the base plane ( i . e . , z=0) to i n f i n i t y . 29 As the integration proceeds the calculations become more and more laborious, thus making numerical integration desirable. However analytic integration can be carried out for the f i r s t two (x,y) integrations with reasonable ease, and t h i s was done in a l l cases to reduce computing time. The t h i r d integration in the z dir e c t i o n proved to be too d i f f i c u l t for analytic evaluation in a l l but a single case. Therefore resort was made to numerical integration i n every case. If cross-sections of the element lying in the x-y plane are examined at various values of z, the resulting t r i a n g l e s w i l l , in general, be found to have di f f e r e n t orientations. For t h i s reason a rotation and tra n s l a t i o n of the l o c a l coordinate system i s convenient i n carrying out the analytic integration. Two cross-sections for d i f f e r e n t values of z are shown below. Consider two rectangular coordinate systems x-y and x'-y* (see figure below) with axes rotated and translated with respect to each other. 01 a point P in space has coordinates (x,y) cr (x' #y*) r e l a t i v e to these coordinate systems. The eguations of transformation between coordinates are given by 3 1 or x '= ( x - a ) c o s c x - (y-b) s m ^ y' = (x-a)sm^ +(y-b)cos^ (5-5) x r x ' c o s ^ + y ' s i n © < + a y - -x'smoc + y' coS°< + b ( 5 ' 6 ) where the angle c*- i s the counterclockwise rotation of the x* and j' axes with respect to the x and y axes. Also the o r i g i n of the x'-y' coordinate system i s located at (o^ b) r e l a t i v e to the x-y coordinate system. The x*-y» coordinate system w i l l be attached to a cross-sectional area a distance z from the base plane and therefore w i l l be orientated d i f f e r e n t l y for each new cross section. The x-y coordinate system w i l l be attached to the base plane thereby f i x i n g i t s orientation. In pa r t i c u l a r , the x dir e c t i o n of each coordinate system w i l l be aligned along the 1 - 3 edge with the origi n at the 1 node. These axes are shown below. 32 Zl= O Using the above diagrams, the sine and cosine of the angle ^ of Eqs. 5 . 5 and 5 . 6 are found to be S i n oc, = y . - y 3 COS <x = X 3 - K ! where for example X 3 i s the X coordinate of the point 3 t and X 3 i s the X coordinate of the point 3 . Also a s h i f t i n g of axes takes place, such that the translation parameters are Q - * , , ^ = jf • The next step i s to evaluate the l i m i t s of the integration. The upper l i m i t on the integration with respect to x 1 (see above diagram on the right) i s x = x , f - ( *3 ~ x* | y and the lower l i m i t i s 3 3 X = • y The range of integration with respect to y 1 i s from zero to The i n d i v i d u a l terms of the s t i f f n e s s matrix K have the form r I -UL O LL ( 5 . 7 ) where UL and L L refer to the upper and lower l i m i t s of the x' variable. The term inside the brackets w i l l be integrated a n a l y t i c a l l y . On examination of the integrand of Eg. 5 . 7 for a l l of the terms of the s t i f f n e s s matrix, i t i s found that there are six d i s t i n c t types which can be expressed i n the following manner. Type 1 : Type 2 : Type 3 : Type H : Type 5 : f ( x . y , z ) ' 34 Type 6 : ^ (x ,y z ) = ^ ( z ) Since in each case Oj (z.^ i s not a function of x* or y* we can remove th i s term from within the bracket cf Eg. 5.7. Let I. - dx,' d y K C T ) The evaluation of X, reduces to finding the area i n t e g r a l of a triangular cross section lying i n the x-y plane, however this area has already been derived in Eg. 3 . 4 . Thus X . becomes I = (A+ BZ + C Z z (5.8) Let I J dx'd y Using Eg. 5.6 re s u l t s i n I. = Jx' oly' x ' c o S <* + y S i n c < + X ( The X , term may be integrated by inspection by recognizing i t as a type 1 i n t e g r a l since X, i s a function of z only. Carrying out the integration with respect to x* and y' and 35 substituting the l i m i t s yields (x 2'+x 3')co5o< + y I 'sm<*)+ x, (A + Bz + Cz*) With the substitution of the expressions for 5"in << and COS «< there r e s u l t s , + X.CA+BZ + CZ7-) (5.9) In an i d e n t i c a l manner the four remaining inte g r a l s may be evaluated to give the following r e s u l t s . I y c J x ' d y ' ^ J (5. 10) I 4 -d x ' d y ' 1 (x3-x.) xt' (x / + X,') + X , 36 y* + 2x, x, z (A + BZ + C z z ) \ y l ^ y ' (5.11) 1 >i 12 x ; 1 - ( y , - y * ) ( x 3 - * . ) [ y * Y x » ' + 2 x ; ) j + (x 3 - x , ) |_yz'j + 2y, [ I , ] - y ' (A+BZ + C Z 2 ) (5.12) x y d x' ol y' 37 + i 2 ~J r (VO -(y.-Ja) >i (x,'+ X2') + (vO^-y^Ty; I. + X. t l . ] - vy, (A + B Z + C Z ^ ) 2 (5.13) The primed c o o r d i n a t e s may be e l i m i n a t e d by using Eg.- 5.5 and the i d e n t i t y so t h a t the i n t e g r a l s c o n t a i n only unprimed c o o r d i n a t e s . The i n t e g r a t i o n with r e s p e c t to 2 remains to be c a r r i e d out. The X values of Eqs. 5.8 to 5.13 which are f u n c t i o n s of z only are m u l t i p l i e d by ^"^0 t o f ° r n i the i n t e g r a n d which w i l l be n u m e r i c a l l y i n t e g r a t e d as d e s c r i b e d i n Chapter 6. 38 CHAPTER 6 SEMI-INFINITE NUMERICAL INTEGRATION 6.1 Goodwin-Moxan^Method The Goodwin-Moran method of numerical integration,* which i s a form of the f a m i l i a r t r a p e z o i d a l r u l e method, i s used f o r the s e m i - i n f i n i t e i n t e g r a t i o n . I t has been noted t h a t the t r a p e z o i d a l r u l e given by the e x p r e s s i o n fa)dt = h 2 ^ i M J (6.1) o f t e n g i v e s s u r p r i s i n g l y a c curate r e s u l t s f o r i n t e g r a l s over a doubly i n f i n i t e range. To take advantage of t h i s phenomenon, a t r a n s f o r m a t i o n t h a t transforms i n t e g r a l s over other ranges to the doubly i n f i n i t e range i s made. For the i n t e g r a l of i n t e r e s t i n t h i s t h e s i s , namely \ (-(x)cJx # the t r a n s f o r m a t i o n * = £ i s used. T h i s transforms the i n t e g r a l to a doubly i n f i n i t e range whereupon Eg. 6.1 can be used. The r e s u l t i n g quadrature r u l e i s expressed as f o l l o w s 39 r J = N = i 4- p J ' J (6.2) The complete d e r i v a t i o n of Eq. 6.2 appears i n Appendix B. The parameter S i s a s c a l e f a c t o r to cente r the summation on the area which makes a s i g n i f i c a n t c o n t r i b u t i o n while a v o i d i n g r e g i o n s which do not c o n t r i b u t e a p p r e c i a b l y to the i n t e g r a l . The parameter P s e r v e s as a node p o s i t i o n e r or i n t e r v a l s pacing parameter analogous to h i n the t r a p e z o i d a l r u l e . In t h i s f o r m u l a t i o n the i n t e r v a l spacing between nodes forms a geometric p r o g r e s s i o n which i s mere l o g i c a l than uniformly spaced nodes f o r an i n f i n i t e i n t e g r a l . 6.2 An a j i v t i c a l _ y e r s u s _ N u m e r i c a l n I n t e g r a t i o n In order to t e s t the accuracy of the numerical i n t e g r a t i o n , a comparison with the exact i n t e g r a l can be made I/* f o r some of the e n t r i e s i n the l o c a l s t i f f n e s s matrix j\ T h i s comparison i s made only f o r one i n t e g r a l whose i n t e g r a n d does not c o n t a i n a d e r i v a t i v e with r e s p e c t to z and which can be i n t e g r a t e d a n a l y t i c a l l y with some ease. The comparison w i l l be made f o r element kj-, of k ^ 40 which i s given i n Appendix A as E (i-v)(o<+/3) k 0+vXl - 2 v ) (6.3) xn 6.2^1 An a l y t i c a l . Integratign_.of k t l Pi The area coordinate ^ i s given i n Chapter 3 by At+Blz+Clz*+ (DI + E L Z ) X + (F^Q.z)y A + B z + Cz ' (3.6) From Chapter 4, i s given as (4-3) D i f f e r e n t i a t i n g Eg. 3.6 as indicated i n Eg. 6.3 and l e t t i n g E(l-v)(* + ;fl) _ f yields (l+v)0"2v) ^ (D.+E.ZXFJ + S.Z) Jxdyolz J (\*zy ( A + B z - C z A X 41 (D, + E i z ) ( F l + G ,z ) [Ve a ( z ) ]o l z o ( l . z ) 1 ( A + B z + Cz. 1 )* Since the integrand i s not a function of x or y, the f i r s t two integrations y i e l d the area of the triangular prism as a function of z. This expression for the area has already been obtained i n Chapter 3 as ZA = A+Bz+Cz. Therefore A.+ B z + C z ' Area(z) = Z The integration with respect to z can now be carried out ( D , + E , z Y F > G , z ) d ; K - ^ " z \ ( l + z ) V A + Bz + C z J ) - 1 2 r «o Oz)* (A + Bz + Cz*) 42 L e t (| +z) ~ "t 1 z J , t 2 ( A * + B*t+C*t* ) Z 9 olt + t T + J t T I 43 and edge d i r e c t i o n parameters (see Eqs. 2.1 and 3.4) a = O b, = -I a» =' O b 2 = | a 3 - | b 3 = 0 A= *'Xy;-y;) + *:(y;-y;) + < ( y ; - y : ) o o + I ( i - o ) B - ^ ( b , - b 1 ) ^ - ( b , - b 1 ) + x;(b,-b 1) + y,'(a,-a,) + y 1 ,(q,-a 1) + y ; ( a , - a 2 ) = 0 + 0 + 1 = 3 [ i - ( - o ; + o + i -o ] + o C = a , ( b 3 - k ) + flt(b,-b,) + a . f k - b , ) 0 + O = 2 U5 A* = A - B + C = O B* =' B- 2C = - I C * = C - 2 D, = yl-y: = - i F, = x ; - x ; = - i E, = b 3 - - i ( i . ' V O . - - i 9, - D,F, = ( - 0 ( - l ) = I <Jt = D,G, +E,F, = (HY-I) + (-O(-i) = z Sa = E.Gt, = ( - O ( - i ) - | g ; » 3 , - 2 9 3 = 2 - 2 ( i ) = o ^ = B * l - 4 A * C * - ( - i ) ( -0 - ( 4 ) (o ) ( 2 ) = ! p 00 — g,' i t = ^ 1 <-> L J , T J z 46 = 1 2 - \ Lz - ? L r r - y ( . 3 4 4 5 7 ) Using the Goodwin-Moran numerical integration scheme with N = number of nodes = 15,; S = scale factor = 2 and P -geometric progression factor = '/^ and using double precision output, gave a value of K^, = (. 34£>5 5 ^ ^ • This guadrature i s only approximate of course but the accuracy attained i n the previous example shows that the numerical res u l t i s an excellent replacement for the exact r e s u l t . I 2 ( 2 ) + ( - l ) + | = 1 2 ( 2 ) + ( - | ) - | 2 2 47 CHAPTER 7 NUMERICAL EXAMPLES Two examples were chosen against which the res u l t s of an analysis using the i n f i n i t e f i n i t e elements could be compared to exact a n a l y t i c a l solutions. The two examples b a s i c a l l y derive from the solutions for a point load on an e l a s t i c halfspace, in one case applied perpendicular and in the other tangential, to the surface of the halfspace. The actual examples considered the halfspace with a hemispherical region removed from the top surface, and with loads applied to the surface of t h i s hemispherical bowl egual to the stresses that would arise from a concentrated lead on the entire half space. The stresses and displacements in the actual examples are therefore i d e n t i c a l to the stresses and displacements i n the halfspace problem in the region outside of the hemispherical cutout. U8 Point L o a d on a rlal£space. F i n i f e E l e m e n t E x a m p l e There were several reasons for choosing such examples. F i r s t was the desire to test the element in an example roughly s i m i l a r to a situ a t i o n where i t was considered to be useful, thus the half space problem. Second was the desire to test the i n f i n i t e f i n i t e elements only and not in combination with regular three dimensional f i n i t e elements, thus the hemispherical cutout. Third the desire to test the element, or at least a reasonable combination cf elements, for loads primarily applied i n the i n f i n i t e (2) d i r e c t i o n and for loads applied perpendicular to the i n f i n i t e d i r e c t i o n , thus the two examples with loads normal and tangential to the surface. 1»9 11... -, fi ? § ,B P1 §,., 3 _ Z „,. ? § J tz I *? I i B .1 f .§,_ ,j? ,9 1 I ,£L £ °. § £? e 3 _ I? Y. 3. Surface Point Load mPerpendiculat The equations for the deflections of points in a semi-i n f i n i t e s o l i d loaded with a point load perpendicular to the surface can be found in [1] as U 1 IT P E + r V _ p lit r •f (7.1) P ZTT i n which 50 Z W X \ R Z x,u Fig. 7.1 Semi-Infinite Solid Loaded by a Surface Point Load Perpendicular to the Surface The stresses due to a perpendicular point load on the surface of a s e m i - i n f i n i t e s o l i d are given in [1] as 0~ = — — ** 2TT R 2 P I CT =7 — — XY lir R z cr x z c r YY R 3 ^ y \ R R + Z R(R+Z)' 3XYZ _ 3P X Z Z 2TT R P J_ 2ir R* 3ZY R 3 ( l - 2 v ) ( 2 R + Z ) X Y R (R + Z ) 2 (7.2) - ( l - Z v ) / Y ^ Z R + Z U V \ R R+Z R(R+Z)7 5 1 cj- - 3P XZ YZ 2.TT R 5 _ _3P YZ_* ZZ o D 5 Z7T R These stresses are used in the generation of the nodal stress vector j? (see Appendix D) . Due to a x i a l symmetry, any r a d i a l section may represent the s t i f f n e s s of the halfspace. However the boundary conditions are most e a s i l y inserted i f one-quarter of the semi-i n f i n i t e s o l i d i s considered. The quarter surface of the hemispherical bowl forms the base plane for the i n f i n i t e f i n i t e elements. This surface i s subdivided into a mesh of f i f t y - f o u r elements and thirty-seven nodes as shown in F i g . 7.2. The V displacements are set equal to zero along the X axis and the U displacements are assigned to be zero along the Y axis. Therefore points d i r e c t l y beneath the load are allowed W displacements only. Appendix D describes the load vector used to approximate the stresses distributed over the base plane. These stresses, given by Eqs. 7.2, vary ncn-linearly from point to point. However to simplify the calculations to determine the 52 l o a d v e c t o r , the value of the s t r e s s vector i s c a l c u l a t e d at each node and then assumed to vary l i n e a r l y between nodes. F i g . 7.3 compares the r e s u l t s obtained from the i n f i n i t e f i n i t e element a n a l y s i s using the mesh shown i n F i g . 7.2 with the exact e l a s t i c i t y s o l u t i o n given by Eqs..7,2. The i n f i n i t e f i n i t e element r e s u l t s are found to vary by 4 ° / o on the average from the exact s o l u t i o n . D i r e c t l y beneath the l o a d , agreement wit h i n 2. / o i s achieved. A l l t e s t s are f o r a u n i t load ( i . e . , P = I ). Example 1 r e s u l t s i n a problem with ninety-seven degrees of freedom and a s t i f f n e s s matrix with a h a l f band width of t w e n t y - f i v e . The c e n t r a l p r o c e s s i n g u n i t (CPU) time on a IBM 360-67 was 96 seconds f o r the e n t i r e problem i n c l u d i n g the numerical i n t e g r a t i o n of the s t i f f n e s s matrix. The numerical i n t e g r a t i o n parameters used i n t h i s example are as f o l l o w s number of nodes N| = I 5 s c a l e f a c t o r S = I.O geometric p r o g r e s s i o n f a c t o r Y" = 0 . 5 53 F i g . 7.2 F i n i t e Element Mesh of Hemispherical Bowl with One Quarter Symmetry F I G . 7 .3 D E F L E C T I O N S OF H E M I S P H E R I C A L BOUNDARY L O A D E D P E R P E N D I C U L A R TO S U R F A C E ( X - Z P L A N E ) 55 7.2 Example 2 - _ Seffli-Ipf i n i t e Solid Loaded by., a_Surf acg_Point Load P a r a l l e l to the Surface The equations for the deflections of points in a semi-i n f i n i t e s o l i d loaded by a point load p a r a l l e l to the surface i n the X d i r e c t i o n as shown in F i g . 7 . 4 are [ 1 ] u = Q (uv) i ZTT E R R z (R+Z) Y = Q (W) X Y _ ZTT E R 3 I -(l-2v) R (R + Z ) L w _ Q (l+v) X ZTT E R* Z , (l-2v)R R X R(R+Z) ( 7 . 3 ) * x.u F i g . 7 . 4 Semi-Infinite Solid Loaded by a Surface Point Load P a r a l l e l to the Surface 56 The s t r e s s e s due t o a p o i n t l o a d i n t h e X d i r e c t i o n p a r a l l e d t o t h e s u r f a c e o f s e m i - , i n f i n i t e s o l i d a r e [ 1 ] o X cr = — 27T R 3 XY Z7T R 3 3 X R2 _(l. 2 vWl-2v)RY3- X ^ R + Z ) ' ( R + Z ) 2 \ R 2 ( R + Z ) 3 X * + ( l - Z v ) R 1 R 2 ( R + Z ) 2 3 Q X*Z 2TT R 5 X Z ( 3 R + Z ) R Z ( R + Z ) ( 7 . 4 ) Q X cr = — — YY R 3 Y Z 1 1 R ' 3 Q X Y Z 2TI R 5 -(l - 2 v ) + ( l - 2 v ) R : ( R + Z ) : I - Y2 ( 3 R + z V R 2 ( R + Z ) _ 3 Q X Z ( T , = — -zz z a T R A g a i n s y m m e t r y r e g u i r e s o n l y one g u a r t e r o f t h e s e m i -i n f i n i t e s o l i d need be c o n s i d e r e d . T h e r e f o r e t h e same mesh a s E x a m p l e 1, F i g . 7.2, may be u s e d . The \/ d i s p l a c e m e n t s a r e s e t e q u a l t o z e r o a l o n g t h e X a x i s . I t c a n be shown f r o m a symmetry a r g u m e n t t h a t t h e V and w d i s p l a c e m e n t s a l o n g t h e 57 I axis are also zero. In Figs. 7.5 and 7.6 a test comparing the results obtained from the i n f i n i t e f i n i t e element mesh shown in Fi g . 7^2 with the exact e l a s t i c i t y solution given i n Egs. 7.3 i s shown. The i n f i n i t e f i n i t e element r e s u l t s are found to vary in the order of 3 A from the exact solution. Example 2 re s u l t s i n a problem with 91 degrees of freedom and a s t i f f n e s s matrix with a half band width of 24. The CPD time on a IBM 360-67 was 99 seconds for the entire problem including the numerical integration of the s t i f f n e s s matrix. The numerical integration parameters are the same as those used i n Example 1. 58 F I G . 7.5 D E F L E C T I O N S O F H E M I S P H E R I C A L B O UNDARY L O A D E D P A R A L L E L TO S U R F A C E ( X - Z P L A N E ) FIG. 7.6 DEFLECTIONS OF HEMISPHERICAL BOUNDARY LOADED PARALLEL TO SURFACE (X-Y PLANE) 6 0 C H A P T E R 8 D I S C U S S I O N Instead of dealing with an imposed r i g i d or free boundary, as has been commonplace in s o i l mechanics problems using f i n i t e elements, a f l e x i b l e boundary formed by the i n f i n i t e f i n i t e elements has been introduced. The purpose of using such i n f i n i t e f i n i t e elements i s not to evaluate stresses at some great distance from the surface, but to r e a l i s t i c a l l y represent the s t i f f n e s s of the sem i - i n f i n i t e s o l i d which has previously been assumed to be zero or i n f i n i t e . Close to the surface and near points of load, regular f i n i t e elements would be used to evaluate the stresses and displacements. The solution of a halfspace problem using regular f i n i t e elements requires numerous elements, many of which are necessary onl:y to reach the a r t i f i c i a l boundary where stresses or displacements can be assumed n e g l i g i b l e . This boundary can be reduced with the introduction of the i n f i n i t e f i n i t e elements since they represent the s t i f f n e s s of the surrounding region. 61 Thus, with the shrinking of the boundary, a large reduction i n the number of degrees of freedom i s possible. In l i g h t of the accuracy attained using a very simple displacement function, a higher order element i s not warranted for halfspace problems. 62 BIBLIOGRAPHY 1. S c o t t , R. F., Fundamentals of S o i l Mechanics, Addison-Wesley, Reading, Mass., 1963, pp. 498-500. 2. Squir e , W i l l i a m , I n t e g r a t i o n f o r Engineers and S c i e n t i s t s , American E l s e v i e r , Hew York, 1970. 3. Squire, W i l l i a m , "Numerical E v a l u a t i o n of I n t e g r a l s using Moran Tran s f o r m a t i o n , " West V i r g i n i a U n i v e r s i t y , Department of Aerospace E n g i n e e r i n g , TR-14, 1969. 4. Z i e n k i e w i c z , 0. C., The F i n i t e Element Method i n E n g i n e e r i n g Science, McGraw-Hill, London, 1971. 63 APPENDIX A ELEMENT STIFFNESS MATRIX FOB LOCAL COORDINATE SYSTEM The element s t i f f n e s s matrix f\ may conveniently be partitioned as follows: SYM : : <A.1) The submatrices are a l l (3x3) arrays and the matrix i s a (9x9) symmetric matrix. The partitioned matrix i s shown in f u l l i n Eq. A.2 in which Zv Z(\-v) 6 4 0~ p s E ( i - v ) r ( l+v ) ( i - 2 v ) where E and 'V are Young*s modulus and Poisson's r a t i o . Each term of the submatrix k Lj must be integrated over the volume of the element. X —> < cn + . _» II 1 >* • —> + • —> K r > " T 5 X M • *J _ i _ • ** Q — > 1 % - r Kl >^ + + a - . X — » • —> • -» o + X • .> X X 1 ^ + N 4 L U II 66 APPENDIX B GOODWIN—MOBAN METHOD OF S E M I - I N F I N I T E INTEGRATION The a p p l i c a t i o n o f t h e t r a p e z o i d a l r u l e t o t h e r a n g e -<?o to + oo i s t h e b a s i s o f t h e Goodwin-Moran method. The s i m p l e t r a p e z o i d . a l r u l e + oo j = +N f(x)dx = h X f ( j ^ ) (B.1) - oo J = "N g i v e s s u r p r i s i n g l y good r e s u l t s f o r i n t e g r a n d s whose d e r i v a t i v e s v a n i s h q u i c k l y a s x —> 0 0 . S i n c e + 00 r +00 (B.2) - 00 - 00 Eq. B.1 c a n be r e w r i t t e n a s 67 + 00 J J = - N to center the summation and avoid adding negligible terms for either large positive or negative arguments. When integrals over ranges other than -oo to +oo are encountered, a transformation to the doubly i n f i n i t e range can be made to make use of the trapezoidal rule's accuracy. I f the in t e g r a l to be transformed i s taken as I = food then a t r a n s f o r m a t i o n of X = 6 ^ + Ol w i l l g i ve I - •f ( e 1 + a ) e 1 o l t - oo S h i f t i n g the variable as in Eg. B.2 gives 6 8 I = \ f ( e ^ + a ) e t t S < J t oo and then approximating the i n t e g r a l in the same manner as Eg. B.1 gives j a + N I = K > e j h * J a ) J Now l e t r s e ^ - Then Also l e t O = c . These substitutions produce the result poo J=+* A f ( x ) J x - S | L r | ^ r i f ( a ^ S r J ) . o j - - N The approximation becomes exact as [S]-> oo and P —*• . The absolute value of the term i s taken so that r can be less than I.O and the res u l t i s unchanged as the + J and - J term are simply interchanged from the case where the r e c i p r o c a l of T i s used. The quadrature rule then becomes 69 + N (B.3) 70 APPENDIX C LISTING OF FORTRAN SUBROUTINE OF THE GOODWIN-KORAN METHOD The t f o l l o w i n g i s a l i s t i n g of a F o r t r a n subprogram which uses the Goodwin-Moran method to evaluate an i n t e g r a l over a s e m i - i n f i n i t e range. Eg. B.3 of Appendix E i s the b a s i s of the subprogram. The input to the subprogram i s e x p l a i n e d i n the l i s t i n g while the output, QGEOM, the value of the i n t e g r a l , i s ret u r n e d to the c a l l i n g program. 71 FUNCTION QGEOM (GRAND) C C EVALUATION OF INTEGRAL OVER SEMI-INFINITE RANGE BY GOODWTN-C MORAN METHOD USING POINTS IN A GEOMETRIC PROGRESSION C A=LOWER LIMIT OF INTEGRAL,IT CAN BE POSITIVE,NEGATIVE OR C ZERO C GHAND=EXTERNAL FUNCTION SPECIFYING INTEGRAND C S=POSITIVE NUMBER SERVING AS SCALE FACTOR C RAT=POSITIVE NUMBER (NOT 1.) USED AS FACTOR IN GEOMETRIC C PROGRESSION C N=NUMBER OF NODES (ACTUALLY 2*N+1 NODES ARE USED) C THE FUNCTION GRAND IS EVALUATED AT POINTS BETWEEN C A+S*RAT**N AND A+S/RAT**N. SELECT VALUES TO COVER C SIGNIFICANT RANGE. DEFINE GRAND TO GUARD AGAINST C UNDERFLOWS. C COMMON/INTPAR/A,S,RAT,N RP=S RM=S QGEOM=S*GRAND (A+S) DO 30 J=1,N RP=RP*RAT RM=RM/RAT 30 QGEOM=QGEOM+RP*GRAND(A + RP)+RM*GRAND (A*RM) QGEOM= QGEOM*ABS (ALOG (RAT)') RETURN END 72 APPENDIX D GEOHAL1ZED_LOAD_VECTOB In the f i n i t e element examples of Chapter 7, the nodal f o r c e s a c t i n g a t the nodes of the h e m i s p h e r i c a l boundary are taken t o be s t a t i c a l l y e q u i v a l e n t to the f o r c e s that would be produced by the s t r e s s e s , from a p o i n t load on the s u r f a c e of the h a l f s p a c e , a c t i n g on t h i s h e m i s p h e r i c a l s u r f a c e . I t i s the purpose of t h i s appendix to convert these given boundary s t r e s s e s to nodal p o i n t l o a d s . The i n and out of plane s t r e s s e s C|_(^B) a c t i n g on the base plane are o b t a i n e d from the Bcussinesq s o l u t i o n of the theory of e l a s t i c i t y [ 1 ] , The area c o o r d i n a t e (3 B i s d e f i n e d f o r the base plane only (i.e.,z=0). Let R be the node f o r c e s i n g l o b a l c o o r d i n a t e s that are s t a t i c a l l y e q u i v a l e n t to . When v i r t u a l node displacements C) P t a l s o i n g l o b a l c o o r d i n a t e s , are produced, the work done by these node f o r c e s must be equal to t h a t done by the s u r f a c e s t r e s s e s • That i s . 73 dA (D. 1) where ^ . ( ^ & ) = tl— = continuous displacements i n the g l o b a l system and A i s an i n t e r p o l a t i n g f u n c t i o n d e f i n e d by Eg. 4.2. The d i s t r i b u t e d s t r e s s e s ^ ( ^ g ) v a r y n o n - l i n e a r l y from p o i n t t o p o i n t but f o r convenience they w i l l be r e s t r i c t e d to vary l i n e a r l y . With the d i s t r i b u t e d s t r e s s e s 0^ v a r y i n g l i n e a r l y from node to node, only two node values of 0|_ , one at each end along the boundary, are necessary to s p e c i f y 0|_ . I f the nodal v a l u e s of C|_ are s p e c i f i e d i t w i l l be necessary to i n t e r p o l a t e l i n e a r l y between nodes i n order to d e f i n e CJ_ . A l i n e a r i n t e r p o l a t i o n f u n c t i o n has been used i n Chapter 4 to r e l a t e the nodal displacements to the displacements of a l l p o i n t s i n an element. T h i s same f u n c t i o n can be used to r e l a t e the nodal a p p l i e d s t r e s s e s to the s t r e s s e s a p p l i e d on the base plane. Then, •4 U A o i n which A ^ i s a l i n e a r i n t e r p o l a t i o n f u n c t i o n and t h e r e f o r e 74 i s i d e n t i c a l to Eg. 4.2. The v e c t o r _P c o n t a i n s the nodal values of s t r e s s e s C|_(^ B)and i s 9 * 1 . Rewriting Eg. D. 1 g i v e s That i s . (D.3) Uo c o n f u s i o n should a r i s e over the d i f f e r e n t A' S . The A'S are i n t e r p o l a t i n g f u n c t i o n s d e f i n e d f o r the base plane only while the M p e r t a i n i n g to i s the area of the base plane. The i n t e g r a t i o n can be c a r r i e d out e x p l i c i t l y and w i l l c o n t a i n c r o s s product terms. I t can be shown t h a t the i n t e g r a l s y i e l d the f o l l o w i n g r e s u l t s ' 1.3 75 F i n a l l y , upon i n t e g r a t i o n of Eq. D.3, the node f o r c e s are given by IX R R R IY IZ 2X < R 2Z R R R 3X 3Y 3Z A iz 2 r + r IX 2.X + r 3 X + r IY 2Y + r 3 Y 2 r + r IZ 2Z + >"3I + Z r + r 3 X r IY + P 3 V + 2. r 2 z (0.») where the terms R^x , R^y and and p are •ix ' P iV a n a P i Z the X # Y and Z components at node i of the g e n e r a l i z e d nodal f o r c e and the d i s t r i b u t e d s t r e s s r e s p e c t i v e l y . The v e c t o r ^ roust be found from the s i x s t r e s s components CTT^ . a c t i n g a t a p o i n t . In tensor n o t a t i o n t h i s i s given as T; 76 where h j i s a component of a u n i t normal t o t h e base p l a n e i n the j lh d i r e c t i o n . , when e v a l u a t e d a t the nodes, become terms i n the _|0 m a t r i x i . e. P z t = TY = « C T ( n, + C r f i x + CC 3 n3 e v a l u a t e d at node 2. The s i x unique terms o f t h e s t r e s s m a t r i x C~ a r e g i v e n f o r each problem i n S e c t i o n s 7 . 1 and 7.2. A v a l u a b l e check cn the g e n e r a l i z e d l e a d v e c t o r may be o b t a i n e d by summing the g e n e r a l i z e d f o r c e s i n any p a r t i c u l a r d i r e c t i o n and comparing i t t o the pr o d u c t of the s t r e s s i n the same d i r e c t i o n and the area over which i t a c t s . From e q u l i b r i u m , t h e s e two r e s u l t s s h o u l d be e g u a l .
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Infinite finite element
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Infinite finite element Ungless, Ronald Frederick 1973
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Infinite finite element |
Creator |
Ungless, Ronald Frederick |
Publisher | University of British Columbia |
Date Issued | 1973 |
Description | This thesis is concerned with developing a finite element model of infinite size to facilitate the stress analysis of three dimensional, semi-infinite bodies governed by time independent linear equations. An element of triangular plan form extending to infinity in one direction is devised. There are three nodes per element, each node having three displacement degrees of freedom. The infinite element is used in conjunction with regular finite elements to represent the stiffness of the semi-infinite solid which has previously been assumed to be zero or infinite in regular finite element models. The infinite element relieves the computational problem caused by large numbers of elements which has limited the use of the finite element method in three dimensional halfspace problems. A large reduction in the number of degrees of freedom is possible with the use of the infinite element because the artificial boundary of previous models is eliminated. The accuracy of the element is tested on two examples whose exact solutions are known. Both involve a semi-infinite solid, one loaded with a surface perpendicular point load, and the other with a surface parallel point load. The results compare favourably with the theory of elasticity solutions. |
Subject |
Structural analysis (Engineering) |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-03-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0050530 |
URI | http://hdl.handle.net/2429/32603 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1973_A7 U54_5.pdf [ 2.76MB ]
- Metadata
- JSON: 831-1.0050530.json
- JSON-LD: 831-1.0050530-ld.json
- RDF/XML (Pretty): 831-1.0050530-rdf.xml
- RDF/JSON: 831-1.0050530-rdf.json
- Turtle: 831-1.0050530-turtle.txt
- N-Triples: 831-1.0050530-rdf-ntriples.txt
- Original Record: 831-1.0050530-source.json
- Full Text
- 831-1.0050530-fulltext.txt
- Citation
- 831-1.0050530.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0050530/manifest