The, UNIVERSITY _of NEWCASTLE ‘AUSTRALIA = Discipline of Civil, Surveying and Environmental Engineering CIVL3830 Stress and Finite Element Analysis Part 2: Finite Element Methods Lecture Notes 2005 Lecturer: Dr Andrew Abbo © Copyright. Professor S.W. Sloan All rights reserved.L INTRODUCTION . 2 ONE-DIMENSIONALLAGRANGE INTERPOLATION 2.1 Linea Interpolation 222QuadntleInerpstsion 223 General Formula for Lagrange interpoation 3 THE FINITE ELEMENT CONCEPT....... 13.1 Bnotgy Stored in an Eatc Rod Soiete to» Uniform Bod Force 53.2Deritioa of Sites Equation or Rod Element using Minimum Pott Baer) 133 Finite Element Ans ofan Elatio Rod Sujet oa Uniform Body Fore. 3:4 Comparinon of Fate Element Solution wth Exact Sostion 4 NUMERICAL INTEGRATION AND COORDINATE TRANSFORMATION . “41 One-dimensional Gris Integration 42Numerclly Integrated Blast Rod Elen 43 Worked Etnpe of Nunedcall Intepated Red Blemert {44 Derivation of Stiness Equations or Rod Hkment sng Vital Work ‘$ THE GLOBAL STIFFNESS EQUATIONS 5.1 Proper ofthe Global tlinss Mase 52Rendwidt $53 Storage ofa Symmes Banded Mati 51 Gauss Binion 55 Choesl Decomposition 56 Worked Example of Canis Decompsison $5.3 Worked Example Incorporating Boundary Condos 6 FINITE ELEMENT PROGRAMMING {61 Spectestion of Finite Element Grid 62 Boundary Conions 63 Assmbly of Element Euston ving Sicring Veaore ‘7 TRUSS ELEMENT . 712A Liner Ths Element 112 Formation of Element ifs wing Virol Werk 73 Namerclletgetion of Element Sines Mae 7.4 Formulation of Element Force Vectr 15 Worked Example of Numerically Integrited Tris erent 716A Worked Brame of Tus: Azali wig Fine Elements 8 BEAMELEMENT « 8:1 A Linear Beam Elenea with Shear Deformation £2 Formation of Linear Beam Element Using Vial Work 8.3 Nomerzl Integration of FlementSiffess Mari, £4 Formation of Element Free Vector £5 Worked Ermpleof Numerically Intepsted Beam Element 9 TWO-DIMENSIONAL INTERPOLATION AND INTEGRATION 91 Linear Interpolation Over Hangle 92 Qusdratilaterpolation Over Tange 93 The Jacobian Matric for Tangles 944 Numeric Integration Formule for angles 95 Linear Interpolation Over Quarters{96 Quadratic Interpolation Over Quadetral 9.7 The Jacobian Mate for Quadratrals 948 Numeral Integration Forma or Quadra AO PLANE STRESS ELEMENTS 10 Srese Strain Relations for Plane Sess 102 Plne Stress Linear Dangle 103 Formulation of Plan Stes Liner Tingle Using Viral Work 104 Numerical Integration of Stites Mats or Pane Stes Linear ange 105 Formation of Fore Vetor or lane Suess Linear Tune 106 Worked Exacple of Nate Integrated Plane Stress Linear Tangle 10.7 Pane Suess Quadeatio Tange 1018 Numerial Integration of Stifaes Matis for Pane tess Quadratic Tangle 1019 Formulation of Force Vecr fr Plane Sues Quadratic angle PLANE STRAIN AND AXISYMMETRIC ELEMENTS .. 11.4 StreseStrnn Relations for Plae Stain 11.2 Plne Strain nest angle 113 Stes Stan Relations for Axsyamcty 314 Assgmmeti Linew Tangle 11.5 Formation of Axiymmeti Linea Tange Using Vital Work 11 Nomerieal otegratin of Sess Mat forAsmeti Linear Triangle 11.7 Formulation of Force Vector for Asiymmeti Linear Tangle 11.8 Worked Fximple of Force Vector fr Aslgpmmetc Linear ingleCHAPTER 1 INTRODUCTIONL INTRODUCTION The finite clement method is now used widely in all branches of engineering. It is especially useful in civil engineering for analysis of structural and geotechnical problems. (One of the biggest advantages of the finite element technique is its generality, many different types of problems may be tackled using essentially the same set of steps. Some ‘examples of finite elements are shown below: mt mes linear displacement model quadratic displacement mode ‘wo-dimensional two-dimensional degrees of freedom 6 degrees of freedom Figure 11: Two-noded truss element Figure 1.2: Three-noded truss element mL al linear displacement model quadratic displacement model two-dimensional ‘two-dimensional G degrees of freedom LD degrees of freedom Figure 13: Three-noded triangular element Figure 1.4: Sicnoded triangular element ys aa linear displacement model ‘quadratic displacement model two-dimensional two-dimensional B degrees of freedom 16 degrees of freedom Figure 15: Four-noded quadrilateral ‘Figure 1.6: Eight-noded quadrilateral clement clementlinear displacement model ‘quadratic displacement model ‘three-dimensional three-dimensianal 12 degrees of freedom 30 degrees of freedom Figure 1.7: Four-noded tetrahedral Figure 18: Ten-noded tetrahedral element ‘element linear displacement model three-dimensional 2A degrees of freedom Figure 1.9: ight-noded brick element Figure 1.10: Twenty-noded brick element ‘The two-noded truss element is one ofthe simplest elements and may be used to predict the deflections and forces in pin-jointed frames. P FFigore 1.11: Finite element analysis of pin jointed truss‘The two-dimensional solid elements may be used to model awide variety of structures. Advantages of the finite element technique include the ability to model arbitrary ‘geometries and complicated loadings. Figure 1.13: Plane strain tunnel subjected to gravity loading In this course, we shall be concentrating on the formulation of the displacement finite clement method for analysis of static linear elastic problems, The technique isnot limited to this area, however, and is widely used for analysis of nonlinear materials and vibration. ‘The steps involved in applying the displacement finite element method to linear elastic problems are essentially a follows: ‘+ Select displacement model. This model wil typically be based on a polynomial and ‘employ Lagrange interpolation ‘+ Discretse the problem into nodes and elements. ‘+ Derive the element stiffness equations using the principle of minimum potential energy or virtual work, ‘+ Assemble the element sifiness equations to form the global stiffness equations. ‘+ Apply the boundary conditions and solve the global stiffness equations for the ‘unknown displacements, ‘+ Compute the stresses and strains using the approximate nodal displacements, Each ofthese phases willbe discussed in considerable detailCHAPTER 2 LAGRANGE INTERPOLATION2 ONE-DIMENSIONAL LAGRANGE INTERPOLATION Lagrange interpolation is used almost exclusively in the derivation of displacement based finite elements. The use of polynomial shape functions permits high order elements to be formulated in which the element stfinesses are evaluated using numerical integration. This feature is especially important for three-dimensional and curved-sided elements where analytic integration of the stifinesses is impossible, 2.1 Linear Interpolation Consider a function u(x) which we wish to approximate byalinear function f(z) over the interval x, : Linea one-dimensional shape functions ‘We may check the validity of (24) by computing the value of f(z) at selected values of x Ate eywie) «22 Dy, BOD, oy HO) Gay Goat Axe 11 —a)), +26, = Another useful check on the shape functions is that they must sum to unity according to 1 = Bay +) NG) + Na) = PNG) = 1 ‘This condition follows from substituting 1 arbitay constant in equation (2.1), where c is an 2.2 Quadratic Interpolation (Consider a function 1) which we wish to approximate by a quadratic function fts)over the interval x, = x = 25. The approximation i) is defined to coincide with u(x) exactly at ey ‘Figure 2.3: Quadratic one-dimensional Lagrange interpolation the end points so that ix) = u(s,) = my, Ce) = wa) = up, and des) = ws) = wy asshown in Figure 23. Iti, of couse, necessary to match i) with u(x) at three distinct Points in order to define the quadratic variation uniquely. The approximation i(e) is writen in terms ofthe interpolation functions N(x), NG) and Ns) according to GG) = Nye), + Neus + NG @s) ‘At node one, we must have ig) = i = NyGe)ey + Walesa + Neds = 11 and hence Ny(e,) = 1 with Nee;) = Na(e;) = 0. Simiay at nodes two and thre, ‘Nofea) = Nats) = Lith NyCe3) = NGs)-= Notes) = NG) = 0 Since the equation ofthe lines passing through nodestwo and three are x = x,andx = x3, the condition N(x) = Nj(,) = Os satisied by guessing that Nj (2) is ofthe form. iQ) = Hema) x3) here dis an unknown multiplier. Substituting the condition that N,(x;) = 1, thisimplios —— “Gp and hence On EEE e5 Using the same procedure at nodes two and three we obtain Nyy = LED) aye) en = Ye = *2) Ne) = EE es Substituting (2.6), (27) and (28) into equation (25), we see that the quadratic ‘approximation i) is defined by @-e-a) , , @—aye- Grea) ee ‘The shape functions N(), N2)and Na(x) vary quadratically between x, and x3as shown in Figure 24. figs) = N ‘Figure 2.4: Quadratic one-dimensional shape functions ‘The correctness of (2.9) may be verified by checking that d(x) coincides with u(x) at xy. and x3 Atx =a,2 ==), LANA, MDs, oy = eae a) Galea) eae a Manny i = A=, 5 =H), , e—aN=2D, Gres ee mse) 2 Atrexg iif = @s—2),, Gs De—a), . @—aes—ad, © ° Ge sane =H)" * GN =99)"? * Ga Vey HD By dircet substitution, we may also readily verify that te NG) + Ne) + Ns) = YN) 2.3 General Formula for Lagrange Interpolation Instead of deriving the shape functions from frst principles it is posible to use the general formula f@) = Y Neu, (2.10) & where (uy land m ~ 1 is the degree ofthe polynomial approximation required. Applying equation (2.11) with m = 3, we obtain the same shape functions that were derived in the previous Section acording io aS oe 2 eee) m= Te 5- eBay ie iss F@-s) _ @-me-s NO = TE" eer jet = Bae -») By dinton te shape ution aay ply the lowing properties (= for ji (a) = 0 for j= 1.2, andj atDMG) =1CHAPTER 3 { THE FINITE ELEMENT CONCEPT16 ‘3 THE FINITE ELEMENT CONCEPT In this section the finite clement concept is introduced by considering a simple one- HE a Ny = SD 2 HD yy (628) ax) E2-H ‘Te compute the displacement at the midpoint ofthe element, where x = L/4, we have ~2xk axe and ly, = 3812 thu, = ‘The above equation indicates thatthe displacement at the midpoint of the element is ‘merely equal tothe average ofthe displacements at its ends, which swhat we would expect intuitively ‘Similar calealations may be performed for element two, where = 1/2) . erp ~ In computing the above shape functions, iti crucial to note that the local node numbers tone and two of element two correspond to the global node numbers two and three. Accordingly, the displacement at any pint in element two may now be written as = Ny + Noid = NU, +NUs ‘To compute the displacement a the midpoint of the element, where SL/A , we have wake Nya 2-2x8 ted £72and edu, +h undu, +10, iL? 16E ‘The finite clement solution for the displacement at various point in the rod is shown in Tble A, x 0 Lis rR 3Li4 L : 0 ama a2 ie ae 16E 8 16E 26 ‘Table 3.1 Finite clement solution fortwo linear elements ‘To compute the strains at any point inside an element, we differentiate the displacements. In element one for example, equation (3.10) gives Ay My eta bre and hence [Note that the strainis constant over each element. This is because the strain isa frstorder differential of the linearlyvarying displacement. Sinilai for element two, we obtain, eg ew eg ae My, 4 = Gb id + Gd = G+ GU, G2) with an, ect G30) aN, 2 ae? 63 and hence‘Comparing the strains between element one and element two, we see that they are different. This implies that the strains are discontinuous between elements. The displacements, on the other hand, which are the primary variables, are continuous throughout the mesh. Generally speaking, the jump in strain between adjacent elements {is reduced as the number of elements is increased, (Once the strains have been found, the stresses may be computed using the stress-strain (constitutive) law. For uniaxial elastic conditions, the axial stress is related to the axial strain via Young's modulus F according to on Ee 6.32) Since the strains are constant across each element, so are the stresses. The stresses may jump between adjacent elements. For element one we have op = Ex SEL. EL EX gg ‘Similarly, for element two we xX ML Om BX GEG ‘To assess the performance of our finite element mode, it i instructive to compare the results from the two-element mesh with the exact solution, ‘The accuracy of the nodal displacements and element stresses is of particular interest, 3.4 Comparison of Finite Element Solution with Exact Solution ‘The equilibrium equation of the rod is given by (3.1). Substituting (3.2) and (3.32), the governing differential equation becomes ait ye £ab +0 "This may be integrated to give where cis an integration constant, ‘The boundary conditions at the free end of the rod, where x = L, are oy = Eilujdx = du/dx = 0. These define the integration constant as = EL/E. Thus the exact solution for the strain distribution ia the rod is ag t= Fu -2 833) Using (6.32), thestrosss are given by ox Bes = RL 63 Integrating (3.33) again, we obtain the displacements asoo iu -¥) G35) where the integration constants er since u = 0 when x = 0 ‘The exact displacements are compared with the two-element finite element solution in ‘Table 32 and Figure 35. x ° us 12 sue L 7 5 ae | ame w@ 16 ae 2 ans 0 pray ma t RE a %exor 0 =14 o ~7 0 ‘Table 3.2; Errors infinite element displacements For this problem, the finite element displacements match the exact displacements at the nodes, butnot atpointsinside the elements, In general, thenodal displacements from finite clement analysis do not match the exact displacements, and itis necessary to refine the ‘mesh by trial and error until the displacement errors are small. Tt also interesting to note that the displacements are never greater than the exact displacements. This is always tric for displacement finite clement analysis of linea elastic materials where it can be shown that the finite clement stifiness matrix is invariably too “stiff”, The fact that the nodal isplacements approach the exact displacements from below is very useful in finding suitable mesh by tral and error, 06 —_____ 05: oa. aE 931 finite clement mm 0 035 050 075 100 ‘Figure 3.5: Comparison of finite element solution with exact solution0 ‘The errors n the stresses are generally much larger than the errors inthe displacements, as shown in Table 3.3 and Figure 26. element 1 clement 2 x 0 La Le 12 | 3u4 L o ama | aera | axes | mua | xu | xrys @aar| %E | ae | mp | me | ma | o Yeenor | 25 0 50 50 0 = ‘Table 3.3: Errors in finite element stresses Unlike the displacements, the computed stresses are most accurate atthe centre of each clement, This is because the finite element procedure attempts to model the exact stress field in a least squares sense. Because the stres is constant within each clement, and no continuity conditions are imposed, itis inevitable that stress jumps, such as those shown in Figure 3.6, occur between elements, This introduces ambiguity ifthe stress is sampled at the nodes. In Table 33, for example, node two basa stress of 3X1/4in element one bat stress of XL/4 in element two. Quite cleat, as the number of clementsis increased, the size of the stress jumps, and the corresponding ambiguity, will iminish. In displacement finite element analysis, its usual to sample the stestes at points inside the elenment and the displacements atthe nodes. The precise location ofthe points fr sampling the stresses dependson the typeof element and willbe discussed in subsequent letures. For our simple problem, however, the finite element stresses match the exact stresses ifthe former are ‘sampled at the centre of each element finite element 0 035050075 1.00 i Figure 3.6: Comparison of finite element solution with exact solutionCHAPTER 4 NUMERICAL INTEGRATION AND COORDINATE TRANSFORMATION30 4 NUMERICAL INTEGRATION AND COORDINATE TRANSFORMATION InChapter 3, the principle ofminimum potential energy was used to derive the element stiffness equations for a simple rod element. ‘To deduce these equations it is necessary to perform integrations over the volume of atypical element. Even forsimple elements these {integrations become quite complex. In this Chapter we describe how analytic integration may be avoided by the use of coordinate transformation and numerical integration. [Numerical integration is now standard in most finite element codes and will be used throughout the remaining Chapters. The advantages of using numerical integration, as ‘opposed to analytic integration, are as follows: ‘+ All integrations are carried out by the computer. For many types of elements, especially those involving two or three dimensions, analytic integration is both tedious and impractical + Elements with curved boundaries can be developed. Analytic integration for these ‘ypes of elements is usually intractable in two or three dimensions, + Properties such as Young's modulus and cross-sectional area may be allowed to vary throughout each element. Because numerical integration permits the material ‘properties to vary from point to point inside an element, itisalso possible to deal with nonlinear problems ina straightforward manner. Although there area number of methods avilable for evaluatingintegrals numerically, the ‘one most often used in finite element analysis is due to Gauss. This method will be used exclusively in the remaining Chapters. 4.1 One-dimensional Gauss Integration Gauss’s formulae are writen for finite integrals over the interval (—1, 1] and specify @ set of coordinates and weights. In one dimension these rules are ofthe form. | Fea = Speeym, (aay ‘where nis the numberof integration (or Gauss) points ,isthe-coordinate for integration point i, and ws the weight for integration point “umber of points coordinate weight exact for polynomial a x om of degree T 0 z T We 1 a Ue 1 0 59 3 59 5 6 5/9 ‘Table 4.1 Data for Gaussian integration The rules ae designed for cases where f(a) is polynomial. Indeed, a rule with n points is exact fora polynomial integrand of degree 2n~1 of less. For example, a one-point rulea gives no error for a constant or linear function, Whereas a two-point rule gives no error for «polynomial up to and including a cubic. Gauss rules may also be used for cases where the integrand isnot spolynomialbut the result wil, n general, be inexact. The sampling points and weights for various Gauss formulae are shown in Table 4.1. Other formulae for higher order polynomial integrands are also available and are given in many numerical analysis tems ‘illustrate the use of these rules, consider the third degree polynomial $0) = 48 +32 +21 ‘The exact analytical solution forthe integral of f(x) Over [—1, 1s ["Wosatrnetans Since f(a) is of degree three, a two-point Gauss rule is needed to evaluate this integral without error. Using the two-point rule we see that J fede = > Fea)w; = 0/8) x 1+ FCA 3) x1 («07 y+ 3.8) + 21/08) + 1) (4-218) + 347.8) + 264.8) + 1) “[F)o (Fe) land the Gauss rule performs the integration exactly, 42 Numerically Integrated Elastic Rod Element To evaluate the integrals that are necessary to form the element stiffness equations it is convenient to introduce a coordinate transformation, 1n this procedure we map the “real” clement onto a “model” element using a standard transformation. ‘The model clement is chosen so that all of the required integrations can be done nuerically in a SMraightforward manner. The shape functions and integrations, which were previously expressed in terms of the real coordinates, are now evaluated in terms of the model coordinates. To illustrate the process, consider the linear two-noded rod element that was, considered in Section 3.2, real coordinates model coordinates Figure 4.1: Coordinate transformation for numerical integration32 ‘The model coordinate, is chosen to range from —1 19 +1 so that Gaussian integration ‘may be applied easly. To derive the shape functions in terms ofthe model coordinate, we again use Lagrangian interpolation. For the shape function associated with node one we we know N, = Oat = 1 andthus guess that My = Ny) = AG -1) ‘where 4 isa constant to be determined. Using the condition that Ny = 1 when & gives 4 = -1/2, and thus 46-1 (42) Similarly, forthe shape funetion associated with node wo, weknow Ny = Oat § = ~1 and sess that Ny = NG) = 1G +1) In this case Ny = Lwhen § = 1 and hence 2 = 1/2, This gives N,=3E+1) (43) The displacement at any point inside each rod clement is again given by the equation an + Nat = [Me safe] Nu= INT “s where u, and ware the nodal displacements, N, and Nare respectively given by (8.2) and (43), and no [Mi A] {in addition to using the shape functions to define the displacements, we also use them to define the coordinates. The 3-coordinate at any point within the element can be written in a similar form tothe displacement according to NT (as) =n n= [Mi maf} x where xTe fs *] and and 3p ae the nodal cordnates. Substituting for the shape fictions and rearranging, we see that the model coordinate & is related to real coordinate x according.8 Bonn by my fe gee ‘When the same shape functions are used for both the coordinates andthe displacements, as in equations (4,4) and (4.5), the element is said to be isoparametric. Numerically integrated elements are often called isoparametric elements in the literature and the two terms are used synonymously. 1t should be pointed here that itis not necessary to use the same shape functions to model the displacements and the coordinates. For example, with ' quadratic rod element, itis possible to use quadratic shape functions to model the displacements bu linear shape functions to model the coordinates. Such an element is said tobe subparamettic. Thistype of element i often used where the geometry ofthe problem is simple and a high order is not required to model the coordinates, Its also possible to use superparametric elements, where the expansion used to model the coordinates is of higher order than the expansion used to model the displacements. Both ofthese types of elements, however, ae much less common than isoparametric elements. Differentiating (4.4) with respect tox, the axil strain i the od i given by (46) an where axial strain in rod cy 2 8 nip ma ut = (1; u;) = nodal displacements for rod Inordertocomputethestran-dsplacement matrix andhence te strain at any point inside am clement, we need to evaluate partial derivatives of the shape functions wth respect to the real coordinate x. Since the shape funetios are now expressed in terms ofthe model coordinate these derivatives ere found using the chain rule according 0 ANE _ d (1, de Ba ebes)g - aN, NodE | d (1, us Ze Ea ~ Aa +o)s 2 BE In the above, the quantity dx/d is found by differentiating (45) to give a, ay, fay, avy] f) an, em sty = Me, -[# ehh) 3s as) = Ee be "Note that dx/a is nown as the determinant of the Jacobian matrix and is usualy denoted 4s detJ. For a one-dimensional case the Jacobian is a 1x1 matrix and hence its {determinant isequal:o the matriitset. Substituting forthe partial derivatives ofthe shapeu functions in (4.8), we see that det J isin fact equal to L/2, where Lis the length of the rod element. (Once the strain has been computed, the ail stress inthe rod is given by the stress-strain relation o.= Ee where F is Young’s modulus. This equation is often written in the more general form of o=De a9) where stress-strain mati, oT = (on) axial sures vector ‘Te derive the stiffness equations for a numerically integrated element, we again invoke the ‘principle of minimum potential energy. Using the Jacobian transformation to change the Integration variable the strain energy stored in atypical rod element may be written as SHAY onde 4f ecteccersat Substtingequstions (7) and earning es 54] cuyepe ons a 4[ erwrene cers “ “ 4o'[ won aaa a0) ‘The external work done by the body force on the element is “ waa [Rude af aNee af xnecua , a 4 Defining b? (aan) ‘The element stiffness equations are again found using the principle of minimum potential energy, This states that35 a a aly Er =0 Where a is the potential energy for the element. Inserting (4.10) and (4.11) and differentiating we obtain the element stiffness equations as where k 4f BDB dets df an) J NT b dots af (a3) and fe (44) Equations (4.12) and (4.13) may be integrated numerically using the Gauss rule to give Sere, oets,v, (as) SINT b dots, (416) k f ‘here is the number of integration points, the value of the weight for integration PointiB; = BCE, isthe strsin- [ee] xk kx? (423) Lm) ay fio] a2 au ef} 425) ‘This is the same as the element force vector given in equation (3.17), which was derived by exact integration. 44 Derivation of Stitfness Equations for Rod Element using Virtual Work ‘The governing equations for finite element analysis may also be derived from the principle of virtual Work. Indeed, for many problems in structural mechanics, the virtual work principle is simpler to apply than the principle of minimum potential energy and is also more general ‘The principle of virtual work states that if@ and d--dd are any two displacements fickis satisfying the displacement boundary conditions, then equilibrium is automatically satisfied provided the following equation holds [steer [ora [sas 25 {In the above we note that e isthe strain vector whose length equals the number of strain ‘components, isthe stress vector winose length equals the numberof stress components, bis the imposed body force vector whose length equals the numberof foree componeats, 1 is the applied surface pressure vector whose length equals the number of force components, Vis the volume ofthe body, and Sis the area of the body over which suface pressures are applied. It should be stressed that, provided the changes in geometry ae infinitesimal, equation (426) is applicable to materials which have either linear or nonlinear stress-strain ‘elations. This means it can be used to derive the governing finite element equations for avery general class of materials. The principle of minimum potential eneray, on the other ‘hand, is only applicable to linear elestc or nonlinear elastic materials which obey Hooke's law. ‘ollustrat how the principle of virtual work canbe used to derive the element stiffness {equations in structural mechanics, we again consider the rod element discussed in Sections 3:2and 42. In this ease the virtual work equation becomes A [ ee (62n38 ‘where the strain vector is ¢ = {e,}, the stress vectors 0 ‘vector is 6d = {u), and the body force vector is b = [R}. (og), the virtual displacement Now from equation (4.7) it follows that T= Bul = WaT and hence eT = duTBT (428) ‘Substituting (47), (4.9) and (428) into the lft hand side of (4.27) we see that * et af area =f de BB ude owtaf BIDBdets dw (629) After isertng (4) ito the right hand sie of (427), we lo obtain ‘| Combining (4.27), (4.29) and (430) gives saToac =f aetwtnds = aota [" nr vaesae (420 asta f BTDBdet J dE w= owt [ NT b det dg ‘Since the variation duin the nodal displacementsis arbitrary andthe above equation must hhold forall such variations, it follows that af wopaudu =a [ NT bets df ‘This expression may be written in the usual form ose ‘where the element stifiness matrix and force vector are defined as, xeaf" oon aeaae toa f Wrause and the nodal displacement vectors These stiffness equations are identical to those shown in (4.12), (413) and (4.14), which were obtained using the principle of minimum potential energy.CHAPTER 5 THE GLOBAL STIFFNESS EQUATIONS0 ‘THE GLOBAL STIFFNESS EQUATIONS This section discusses various methods for solving the NXN set of global stiffness equations KUSF Where K isthe global sifness matrix, Fis the global force vector, and U are the unknown lisplacements. In many applications ofthe finite element method, the global stifiness ‘matrix has a number of key characteristics which enable vet ecient solution provedures to be developed, These procedures are necessary because practical applications of the finite element method technique may involve very large sets of equations, It is not ‘uncommon for atypical finite element mesh to involve 10,000 or more unknowns 1 Properties of the Global Stiffness Matrix For many problems in strctural mechanics, particularly those involving linear elastic ‘material behaviour, the global stfiness matrix has the following properties: + Itissparse and banded ‘This means that many ofthe entries in K are zero and that ‘the nonzeros are clustered ina “band” about the diagonal, Ignoring al of the zeros ‘outside the band leads o significant savings and isa property which should always be exploited, + Itissymmetric. This means thatthe enteies in K are symmetric about the diagonal s0 that Ky = Ky. Exploiting symmetry typically seduces the number of operations required to solve linear equations by a factor of four and also reduces the storage demands by a factor of two, ‘* It is positive definite This means that the scalar quantity UTKU > 0 for all N-mensional vectors U xO. For lina elastic nays the quantity JUTKU actualy corresponds othe stain energy othe gntem IF Ks posite deta thi impli thes an inverse andthat the unknowns Ucanalvaysbe found using Gauen clinaton witout row interchanges le. without pong. Sine ne Svein needet no entra nonzrosoutsie th diagonal band ar creed an it epee o vero the origina K matrix wi acre form, Note hata postive deans matrix aio has the following properties +The lagonal ets ral posi, otha > 0 fr f= 1,208 +The absolute wae af he largest entry onthe diagonals notes tha the abso ae ofthe rgesteney nye clin tbe mar Mathematical this condi expressed os magi | Dal © Foreach i= j,k < Ky Ky Because of the above properties, most finite element codes solve the global stiffness equations using some form of decomposition procedure which does not need row interchanges. Common algorithms for solving finite element equations include banded, profile and frontal solution schemes. These may all be regarded as variants of Gause ¢limination and differ only in the way in which the sparsity and structure of the linear ‘equations are exploited, In these notes we will concentrate on banded solution methodsa ‘which use Choleski decomposition. The steps involved i this process willbe deseribed in ‘greater detail later in this section. 52 Bandwidth ‘Consider the mesh of three linear rod elements with four nodes as shown below. Figure 5: Three element mesh numbered for low bandwidth ite that the global stiffness matrix: + element 1 conte to equation 1 and2 + clement 2 contsbutes to equtions 2 and3 + slement 3 contributes to quitions3 and 4 ‘Ths the gob stitfaes mati fo this mesh Wile ofthe form shown in Figure 2 [XS symmetric Kar ON BN 0 Nk Re bod igure 5.2: Structure of bal stitbess matin wih low bandwidth numbering ‘We note that this is symmetric, sparse and banded. For a matrix with N rows the semi-bandwidth, B, may be expressed mathematically as B= max |b) (ay where bythe row semi-bandwidth for row i the dtference between +1 and the column index ofthe frst nonzero entry in row. For the stifnes matrix shown ia Figure 52, the senivhandwith i given as B = max (2) = max(1,2.2.2) = 2 ‘The semi-bandwidth is dependent on the manner in which the nodes inthe finite element grid are numbered. To illustrate this point, let the previous mesh be renumbered as shown below.2 [Figure 5.3: Three element mesh numbered for high bandwidth For this case we note that: * element 1 contributes to equations 1 and 4 ‘+ element 2 contributes to equations 4 and2 + element 3 contributes to equations 2 and 3 ‘and the global stiffness matrix is ofthe form shown in Figure 5.4 be PRs \.symmenic 1 0 Ka 1 0 Kn Ke 2 N[Ku Ke 0 Kul R—,—_} Figure $4: Strocture of lob sitness mat with high bandwidth numbering With this new numbering the semi-bandwidth is B doubled. ax(1,1,2,4) = 4 and has thus {In order tobe efficient, the nodes should be numbered so that the bandwidth is kept as ssnall as possible. For a banded equation solver this not only reduces the amourt of memory that is required to store the stiffness matrix, but also reduces the solution time. Fora given finite element grid, the semi-bandwidth may be computed using the formula BedD+1) 62) where ds the number of equations (degrees of freedom) per node and D is the maximum Positive difference between node numbers on any element in the grid. For the mesh numbering in Figure 5.3, we haved=1 (one equation pernode) and D=3 (for clement 1 which has nodes 1 and 4). Thus B=1(3+1)=4, which is the same 28 what we observe from Figure 5.4. Sinced remains fixed for the given ype of problem, equation (5:2) suggests that the semi-bandwidth is determined solely by the quantity D. Thus, to ensure ‘that the semi-bandwidth is kept small, che grid should be numbered so thatthe maximums lifference between the node numbers for each clement is kept to a tinimumm, 53 Storage of a Symmetric Banded Matrix Because of its special structure, a symmetric banded matrix can be stored very ctficienty, ll entries outside the semi-bandwigth may be discarded (as they ee all 2er0) ‘and only the upper (or lower) triangle inside the band is stored. ‘The total storage“ requirement for this bandwidth solution scheme is typically equal to the product of the semi-bandwidth and the number of equations. To ilustrate this point, the global stiffness ‘matrix for the low bandwidth numbering of Figure 5.1 maybe storedas chown n Figure 55. Since the semi-bandwidth for this cases, the total memory requirement is only Sstorage units [StS symmetric * Ky ar Kas. Ko Kn Oke EK, D> lon ky 0 ON Ke Ky Kes Ku bs be dummy element which is not used Figure 55: Storage of global stiffness matrix with low bandwidth numbering For the high bandwidth numbering of Figure 53, the memory requirement is doubled to 16 storage units as shown below. nS. symmenic xox Kn Ka 0 Ku — s+ fumay element which is not used ‘Besemi-bandwd storage required=NxB=4xé=16 Figure 56; Storage of global stifness matrix wth high bandwidth numbering Because a number of zeros are stored, the poor numbering clearly results ina waste of computer memory. It will also result in a les efficient solution, since the number of operations (as measured by the number of multiplications) involved in a bandwidth solver is soughly proportional to NB $4 Gaussian Elimination {In practical finite element codes, the global stiffness equations are most often solved using some form of Gaussian elimination. This isa direct solution technique which does ‘not require iteration. Note tha the inverse of the global stiffness matrieisnever computed explicitly, since ths destroys the banded nature of the equations and leads to large storagedemands. The key idea behind the Gauss procedure isto convert the set of equations to "upper triangular form. Once this is done, the equations may aze solved by backward substitution, ‘Tp illustrate the Gaussian elimination process, consider the simple set of matrix equations shown below 1 oy; 4 Ku=|1 3 20} =P 423} 63) jo 2 24 (Us, 10} ‘To convert (5.3) to upper triangular form, we need to create zeros below the diagonal in columns one and two. Starting with column one, a zero inthe frst column of row two can be obtained by computing Ky for j= 3,23 32-22 19) = [0 25 2} ‘Thus we subtract an appropriate multiple ofrow one from row two to form anew row tw, ‘Note that the multiplier of 1/28 obtained by dividing the first entry in row two by the fst entryinrow one. The denominator ofthis multipliers called the pivot which, quite clearly, ‘can never equal zero ifthe equations have a solution. Moving now to row three, there is already a zero incolumn one and thisrow does not need tobe modified. This result is also obtained if we follow the above procedure to obtain Ryn ky, for jaa 4 = 2 9)-92 19 = 23) ‘The right-hand side vectors also need to be modified simultaneously using the same multipliers to give yy dxaeu =0-0x4 0 “After these operations, the orginal set of equations have now become 21 oY) fa lo 25 ava = dur lo 2 2}[Ys] 10) Note that row one is unchanged, whereas rows two and three would, in general, be modified. To createazerobelow the diagonal in column two, rows one and two now remain ‘unchanged with row three being modified according to KG =KY- 2eky for 1,2,345 (0 2 2}~os)0 25 2} 00 os} ‘Thus we now subtract a multiple of the (modified) second row from the (modified) third row to obtainyet another modified third row. Taemultplierin this case isfoundby dividing the second entry in row three by the second entry in row two, Similarly forthe right-hand side we obtain FS =F} ~08F; = 10-08 x 11 = 12 ‘This completes the conversion ofthe equations to upper tiangularform, as they now read 21 oyu) 4 lo 25 2 /4U2p = dur 64) oo oles} [a ‘The solution can be obtained by backward substitution, As the name implies this process starts withthe last equation and works baciowards as shown below 04U;=12 > Uy=3 25U,+2U;=11 > 2U,+ Inspection of (54) reveals that Gaussian elimination preserves the bended structure of the equations, since no nonzeros are created in locations that were intially ascupied by zeros. This property permits the original matrix tobe over-vritten by the upper triangular form during the elimination process without the need for additional storage Ifthe solution t0(5.3) i found by computing the inverse of K directly, we obtain A rapa) py ue {Ubewip=[1 2 -2Hish = 42 65) U5] 1-2 2silio) |, Clearly K~1 isnot banded. Since new nonzeros have been created at locations (1,3) and Gu) itrequires more storage than the upper triangle formed by Gaussian elimination, 5.5 Choleski Decomposition ‘A widely used method for solving the global stiffness equations is known as Choleski decomposition. This procedure uses the result that asymmetric postive definite matrix K can always be written ia the form KLE where Lis lower triangular matrix Ty 0. 0 lear La ~ 0 Legs Lg ~~ Lon46 ‘The Key to Choleski method is thus computing I, Once this is found, the solution to the system KU = Fis equivalent to the solution of UUs 66) ‘This equation may be solved in two stages. Inthe frst stage we et Y = L™U and solve LY=F 6n for ¥using forward substitution. The process of forward substitution, which isthe opposite of backward substitution, can be used inthis equation because Lis ower triangular, Ta the second stage, the solution U is found by using backward substitution to solve ue 68) ‘This.to stage approach may appear unnecessarily complicated compared with Gauss climination, Note, however, that once we have performed the decomposition fo obtain L, we may solve for various right-hand sides merely by using forward and backward substitution. This feature is convenient in many engineering applications. For example, itis often necessary to analyse the behaviour of a structure, whose global stfiness matrix remains constant, under many different load conditions, ‘The Choleski procedure may be summarised a follows ‘Decomposition Phase for Computing 1 Set ty = VR Fori= 2m set by = Ky/Ly For i = 2j...n-t do steps 4 and s 2 Fores Quam set ‘Backward Substitution 1 SOU, = Ya/bon"7 2 Pele mtet mt Ye fn 340] i In the above algorithm, the banded nature of the matrix K has not been exploited. ‘Moreover, we have assured that additional storage is avaiable for holding L. In practical finite codes, L would overavrite K and advantage would also be taken of the banded structure to reduce the numberof arithmetic operations. ‘5.6 Worked Fxample of Choleski Decomposition ‘Tillustrate the Choleski procedure, consider the system of equations (5.3) which were solved previously using Gaussian elimination, 1 oyu, xu=|1 3 210, o 2 2] || La = Kn/ly 0/8 =0 2 4 [s- Za BR = HR S$. j=3 Ly [es > bls |- 78 (2-0) = 2/275 Ars 6 ts [ro- Sa] SRB = iB ‘Thus the L matrix is Ro check : LL? as required - Yay] — 1a‘Thus the solution vectoris UF = [1 2 3] 8.7 Boundary Conditions Infinite clement analysis of structural mechanies problems, there are two common types of prescribed boundary conditions + A zero prescribed displacement. This occurs st 2 fixed support where the corresponding displacement U; + Annonzero prescribed displacement. This occurs when a support is allowed to move bby prescribed amount, c, so that U; = c. An example where nonzero prescribed displacements occur is when a rigid footing is pushed into a soft soil mass ‘There are a variety of methods available for dealing wit prescribed boundary conditions in the solution ofthe global stiffness equations. In these notes, we use a very effective deletion method which is not only simple to implement but also very efficient since all prescribed degrees of freedom are removed from the solution phase. To illustrate the process involved, consider the simple set of stiffness equations shown below Fin Kz Kis] (Ux) (Fy Kn, Kaz Kas|{Ua} = {Fe [Ks1 Kaz Kas}[Us] [Fs Let Ua be prescribed to equal a known value e. Writing the equations in full we have KyUy + Kye + KyUs KyUy + Kye + KyUy = Fy + Ry 69 KU, + Kye + Kyy = Fy Note that the reaction R, must be added to the right hand side of the second equation ‘because Up is prescribed. Rearranging these equations we obtainKn Ke ][Ui] | fF ~Kne Ki Ke fe fe 6.10) KyU, + Kno + Kyl and (1) = Fan Sgt - Tb solve the equations (5.9), we first solve (5.10) inthe usual way which yields U7, and Uy Since U; is already fixed we can then use equation (5.11) to compute the nodal reaction Rp. Thus, in this procedure, we simply delete all ows and columns associated with the prescribed value and, at the same time, modify the right-hand side of the remaining equations. Ifthe nodal reactions are needed, each row with a prescribed value must be stored prior to deletion, It's interesting to note that if the value for U, has been prescribed to be zero, which is ‘the most common case in practice, then the right-hand side force vector does not need t0 bbe modified at all and we simplysolve Ky Ky] [] _ [Fr Ky Ky |] Us * ] Fs 6.2) ‘The reaction due tothe fixed value of Ucan, ofcourse, again be computed using equation oan). In cases where more than one displacement is prescribed, the same procedure is followed. For example,let U, = cand U, = c,, Writing out the equations in fullwe have Kut, + Kae + Kya = Fy +R, Kae, + Keg + Kaa = Fy + Ry 13) Kacy + Kygey + KyUy. = Fy ‘The equation tobe solved fr Uys [Bs] = fs ues ~ Kae) nd the reactions corresponding to U; and U,are given by Kuti + Kez + Kus ~ Ry = Kuty + Key + KnUy — Overall, the deletion procedure is efficient since the number of equations to be solved is minimised. Moreover, when dealing with fixed supports where the correspondingpreseribed values are ero, no modification ofthe force vector is necessary and the method is relatively straightforward to implement. 5.8 Worked Example Incorporating Boundary Conditions ‘An example of an elastic rod which is fixed at one end and subject to a uniform Jongitudinal body force Xisshown in Figure 5.7. Two linear rod elements of lengthL/2 are ‘used to model the rod, Figure 5.7: Elastic rod subjected to a uniform body force ‘This example was considered previously in Chapter 3 where itwas shown in equation (3.24) thatthe global stiffness equations for the finite element model are given by’ 1a ayy; 1 reals 2 alfusb. auth, a 0-1 1) [05] 1 [AC the fixed support we have the boundary condition that U7, = 0. To incorporate this condition we delete row one and column one of the matrix equations to give [2 ty Asmentioned in the previous section, no modification ofthe force vectors necessary since the prescribed value is ro. Solving these equations we obtain r= (U, vu {0 me ¥| 2E Because ofthe fixed support, the first equation of (5.14) is modified to Uy ~ Avy + ox vy AE +n, 3ELE sives the unknown force Ry 28 ~ HA, A ALR L 4 XAL ‘The applied forces and reactions for this problem are shown in Figure 5.8. Note the net reaction at node one is ~ KAL + ZAL/4 = -3¥AL/4 and that, most importantly, SAppted Forces + Reactionsst ‘This equation isa fundamental expression of equilibrium and must always be satisfied in any structural finite element analy oe Figure 5.8: Applied forces and reactions for elastic rod under uniform body force ‘Suppose that instead of being fixed, the lft support is moved by an amount U; = cto the right. For the boundary conditions shown, this implies that the whole structure is subject toa rigid body translation, as wel asthe applied body forces. To incorporate this condition, columa one of (5.14) is transferred to the right-hand side and used tomodify the {oree veetor. Simultaneously, row one is deleted to give 2-7/0) 2) (_2e4) [ALE , BAe) 28 = 4ek | en L 4 « “) ae 1 es 1 ° Solving these equations yields the sokution valu uo-[o Bre H+ ‘As expected, the rigid body translation U, = esimply adds the same constant amount to all the displacements that were computed previously forthe case U, = 0. Usingrow one ‘ofthe global stiffness equations the reaction a the left support remains unchanged ands given byCHAPTER 6 FINITE ELEMENT PROGRAMMING6 FINITE ELEMENT PROGRAMMING This Chapter discusses how finite element analysis is implemented. Ta perform a finite clement analysis, the user needs to specify the element topology (geometry) of the problem, the boundary conditions, the applied loading, andthe material properties. Each of these steps will be discussed in detail. It should be noted at this stage that there is no ‘uniform standard for preparing finite element data, and each package typically ses its wa format, In this Chapter, we describe the format that is used for the Newcastle Finite Element Library, 641 Specification ofa Finite Element Grid ‘Tobe able to solve a finite element problem on the computer, we need to specify the finite element grid precisely. Typically, we must supply the following information ‘+The total number of nodes ‘©The total number of elements, ‘+ The geometry ofthe problem. This is done by specifying nodal coordinates ‘+The topology ofthe finite element grid. This is done by specifying the number of rhodes for each element and alist of nodes that define it ‘+The material properties. This is done by specifying the numberof material types and the material properties foreach material type. This information isinked othe finite clement grid by assigning a material type (typically an integer) to each element. ‘+The number of degrees of freedom (vatiable) at each node, "The boundary conditions + The applied loads, Different types of problems, of course, require different types of input but the above is usually abareminimum. To illustrate how this information is specified in practice, consider ‘a mesh of linear elements as shown below. ®@ = Uy slobal freedoms Figure 6.1: Mesh of two-noded rod elements ‘The rod is subject to a uniform body force X andis of uniform cross-section. For this grid Wwe have the basic control data + total number of elements=7O7ELS=3 ‘+ number of degrees of freedom/nod: + total degrees of freedom: + umber of material ypes=NUMMAT=1 (E andA are constant)55 + mumber of nodesielemer -NODEL=2 For this example, the geometry of the problem is specified by the s-oordinats of the nodes. A general scheme, which s applicable to ny number of dimensions, stores the ‘coordinates in an array COORD which has size (TOTNOD,DIMEN) where DIMEN isthe ‘number of phyicel dimensions inthe problem. For the above one-dimensional case, DIMEN equals one and all thex-coorénate ae stored inthe fist column of COORD. ‘With reference 1 Figure 6.1 we ths have FCOORD(,1)) fF1] frsvordinate of node 1 Jcoorn@,1)| _|¥2] _ |s-coorainate of node 2| * |COORDG, 1)| ~ fr] ~ |s-coordinate of node 3} [coord (,1)] |e] [p-coordinate of node 4| coorD ‘The topology ofthe grid is specified by listing the nodes which define each element. These nodesare usually specified ina consistent manner andin the above example itis convenient to arrange the nodes in ascending sequence of their 2-coordinate. ‘The nodes for each ‘element are stored ine singe row of the array ELTOP as shawn below ELTOP = |ELTOP(2,1) ELTOP@,2)| = |1 4]—nodes for element 2 [ELTOP(A,1) ELTOP(,2) 1] nodes for element 1 leL TOP.) ELTOP(3,2) 14 2] nodes for element 3 To hold the nodes forall elements the array ELTOP must be of sie (JOTELS,NOD! ‘The nodes defining an element are found in ELTOPU, J), where J=1,2..0.NODEL. ‘The array ELTOP may also be used to store additional information about each element. {In particular, it may be used to store the number of nodes for each element as well ats ‘material type. The size of ELTOP is then (TOTELS, NODEL +2) and for our mesh of linear rod elements it would look ike 12.3 1]—element 1 121 4/“clement 2 124 2|~element 3 ELTOP ‘Columns one and two now contain, respectively, the material type and number of nodes for each element. Columns three and four old the element node numbers and aze again stored row by row. With this storage scheme, the material type for element Iis found in ELTOP(,1), whilst the number of nodes is found in ELTOPU2). “The nodes defining clement fare found in ELTOP(, J), where J=3,4...LTOP(T2)+2. For the one-dimensional rod element, the material properties that need to be specified are Young's modulus, E, and cross-sectional area. A simple method for storing these is to create two vectors, say E and A, of size NUMMAT. ‘The Young's modulus and cross-sectional area for element /, with material ype M==ELTOP(,1), are then found in E(if) and A(M) respectively 6.2 Boundary Conditions In order to specify the boundary conditions.we need to define which of the global degrees of freedom are prescribed, At this point tis again wseful to distinguish between56 the local freedom numbers, which pertain to an element, andthe global freedom nambers, Which pertain to the overall mesh. For a mesh of three linear rod elements shown in Fig 5.1, each element has the local and global freedoms shown below u how A048 <2 soba 0S) 0G 0 0GEEEO Tr 7 OT 3 T T tocat Figure 62: Local and global freedoms for meth of linear rod elements Note that ‘+ element one has local freedoms 1 and 2 and global freedoms 3 and 1 ‘+ clement two has local freedoms 1 and? and global freedoms 1 and 4 ‘+ clement three has local freedoms 1 and 2 and global freedoms 4 and 2 1m general, the convention is thatthe global freedom numbers are related to the node ‘number Iva the formulae ‘+ for onesdimensioa: global freedom =(U)= I ensions: global freedoms =(U) + fortwo: 2 ‘Thasforour grid of linear rod elements, the global freedom numbers correspond precisely to the node numbers. ‘There are a number of different strategies for incorporating the boundary conditions in a finite element program, but here we will use the approach based ona nodal freedom array. The nodal freedom array, NF, is of size (TOTNOD, DOFNOD) and indicates the equation numbers, thats the final row and column numbers in the global stiffness equations, associated with each global freedom, If NF([J)=0, this indicates that the Jth degree of freedom at node Fs suppressed, IfNF(LJ)=K. this implies that the Jth ‘degree of freedom at node is assembled into the Kth row/column in the global stiffness ‘equations. 4 2% globat Freedoms ® ® + F NFamy en Figure 6.3: Nodal freedom numbers incorporating boundary conditions For the mesh of linear rod elements show in Figure 63, the nodal freedom array is given by INF. 1) INF, 1) |NFG,1) LF.)Ea By convention, the entries in the nodal freedom array are formed by running through the nodes in ascending sequence and allocating numbers to each degree of freedom that is tunrestrained, All fixed degrees of freedom are given a value of zero. The nodal freedom, array holds the nodal freedom numbers aad is used to form a steering vector for each clement. Allowing forthe fact that prescribed degrees of freedom are never assembled, the steering vector tells us where the clement stiffness matrix entries should be inserted in the global stiffness matrix. 63 Assembly of Element Equations using Steering Vectors ‘illustrate the construction and application of steering vectors and the nodal freedom, array, consider the mesh of linear 40d elements shown in Figure 6.3. Element one, for txampl, is defined by the nodes three and one in ELTOP. Thus the corresponding nodal freedom numbers are held in [NF3.1) NF 1] and the stering vector is 0 4 STEER! = {NF(,1) NFU} Elements two and three, on the other hand, are defined by nodes one and four and four and two, respectively, in ELTOP and hence STEER? = (VF(1) NF(4,1)} STEER? = (VF) N21) ‘Tosee how the steering vectors can be used to assemble the global equations, we note that the element sifess matrices are ofthe general form e-fha] ef] fi] ‘where the superscripts denote the element number. For the ease of element one, the clement stifiness entries are mapped into the global sifiness matrix according to 123 o1 if, 0 0) 1 Of Ha) ige : ¥ “Tee " xeao a ol 3f0 0 o [Note thatthe entries above and to the left ofthe element stiffness matrix are simply the steering vector which indicate that ‘+ Ff; ismot assembled into the global stiffness matrix sine its row and column indices ‘rom the steering vector are both zer0 + Elz is not assembled into the global stiffness matrix since its row index from the steering vector is zero ‘+ Ej, snot asembled into the global stiffness matsx ince its column index from the steering vector is zero58 +H lsassembled into X;, of the global stiffness matrixsince its row and column indices fom the steering vector are (1) Similarly, for elements two and three, we have where + HyisascembledintoK,, ofthe plobal stifinessmatixsinc its row and cola indices fiom the steering vetor are (11) + Hlsassembled into Kj ofthe global stifinessmatixsine its ow andcolumn indices from the steering vetar ae (1.3) + Byisassembledinto Ky, ofthe plobalsifnessmatixsineitsrow and column indices fom the steering vetr ae @.1) + Ehisassembled into Ky ofthe global stifiness matrixsince itsrow and column indices ‘fom the steering vector are (3,3) snd 12 f+, 0 Skee) 0 3, ae where +) jnassembled into K ofthe global stfiness matrixsinceits tow and column indives {rom the steering vector are (3,3) + hpisassembled into K, ofthe global stfiness matrixsince its ow and column indices from the steering vector are (32) + Hisassembledinto Ky ofthe global stifinessmatrxsincetsrow and column indices fom the steering estoraze (23) + ,isassembled into Kof the global stiffness matrix since its row and column indices ‘fom the steering vector are (2,2) ‘The stecting vector i also used to assemble the element force vectors into the global force vector. Noting that the element force vectors are of the general form off9 for element one we have aff off insert, fe) SS p-adto slo where ‘© flisnot assembled into the global force vector since its row index from the steering vector is zer0 ‘+ ffisassembled into F, ofthe global force vector since its row index from the steering vector is (1) Similarly for elements two and three we have 1fA+8 insert St p=2} 0 34 where ‘+ ffisassembled into F ofthe global force vector since its row index from the steering ‘vector is (1) + Alisssembedinto ofthe global force vector since isrow ndexfrom theterng ‘vector is (3) and , ifar8 e a peal 8 sleek mere + Alisasembled into F, ofthe global force vector since is rowindes fom the steering vector is) © flisassembled into vector is 2) 9 the global force vector since its row index from the steering ‘Thus by creating the nodal freedom array, and hence the element steering vectors, we have incorporated the boundary condition U = into the global stiffness equations which are ‘of the form60 1 fm fa+A Ky | 403 a (2) +] [yy arf ‘These equations may be solved to give U, Usand Uy Note that 6.1) could also bederived Dbyassembling all ofthe element stiffness and force Vector entries into the global stfiness equations and then deleting the rows and column associated with the fixed displacements, This procedure, howeves, is awkward to implement ina finite clement code (since rows and columns need tobe shifted to remove mull entries) and isalso wasteful of computerstoraze.TRUSS ELEMENT CHAPTER 7@ ‘TRUSS ELEMENT ‘The behaviour of pin-jointed trusses, in which the members are assumed not to transmit ‘moments, is ideally suited to analysis by the finite element method, Although this type of structure may be analysed using a variety of manual techniques, the finite clement approach has the advantage that it gives both the joint deflections and the member forces Moreover, trusses of variable cross-section present no special difficulty and the effets of -. loading by self-weight can aso be incorporated, 7.1 A Linear Truss Element ‘The coordinate system fora two-dimensional trus elements shown in Figure 7.1. The J matric and hence its determinant is equal tothe matrix itself. Note that both the Jacobian matrix and the ttrain-displacement matrix B are constant throughout the element since the shape functions which describe the geometry and displacements ofthe truss clement are assumed tobe linear. (Once the strain is known, the axial stress in the trus is given by the strest-strain relation Oe = Eee (7s) ‘where £ is Young's modulus. This equation is often written in the more general form of =De ‘where the axial stress vector o and the stress-strain matrix Dare defined as 716) 172 Formula ‘Toderive the element stifiness equations we invoke the principle of virtual work which, forthe two-dimensional truss element of Figure 7.1, may be written as 4 [ Segond = PB + P60, + Peps + Pysbrn aan Inthe above, A denotes the cross-sectional area ofthe truss, Seqisthe virtual exal strain, {is the axial stress, (bx ,dy,) ae the virtual displacements at the nodes, and (P,P) are the point loads applied at the nodes. Note that we have assumed thatthe truss is loaded only by point loads applied at the nodes, and not by distributed body forces such as6 selfweight. The latter may be included by considering the appropriate contsbution to the right hand side of the virtual work equation, but will not be considered here for the sake of simplicity, yuation (7.17) may be writen in the usual matrix form as Af oeece = Pon 8) where € = (¢) ithe strain veetor, © (Pa Ba Pa Ba) bul = (bu 6%, ) is the stress vector and P’ Now and De ‘where Bis the strain-displacement matrix and Dis the stress-strain matrix, Substitutingin (7.18) and rearranging gives ota [wince at = 507? 019) Since the variation éu in the nodal displacements is arbitrary, and (7.19) must hold for all such variations, it follows that af BYDBdery dé u = P ‘This expression may be written inthe usual form kunt where a J BTDB dots dé (720) Pp and the nodal displacement vector uis given by (7). Having established the governing equations for the trussclement, we ean now proceed with the details of their evaluationo ‘73 Numerical Integration of Blement Stiffness Matrix ‘Th derive the stiffness equations for a mumercally integrated element, we substitute (7.12), (7.14) and (7.16) ito equation (720) and apply numerical integration. This gives k= 4S BrDR, set3,w, 21) ‘where isthe sumber of integration points, ws the value of the weight for integration Point, B= BE; )isthestrain-splacement matrixevaluatedat = &, dst J, = detJ(@,) isthe determinant of the Jacobian evaluated at # = &, and the value ofthe model coordinate for integration point In te above itis assumed thatthe area and sfiness propertics are constant across cach element. If this isnot the ese, then (7.21) would be ‘placed by an expression ofthe form k= Sai, ders, Aw, (2) where 4, = AG;) and D, = D(,) are the area and streststrain matrix, respectively, computed at § = £, The generality of (7-22) ilustrates the ease with which complicated properties can be incorporated into an element when numerical integration is used, 74 Formal tion of Element Force Vector In the analysis of trusses, the loads are usually applied as concentrated forces at the ‘nodes. This permits the global force vector to be assembled dizecty, without having to ‘compute the element force vector. One common exception to this rule is when the self ‘weight ofthe truss elements are included. In this case, which i not covered here, the nodal {orces arising from the distributed body force must be computed element by element and ‘then inserted into the global force vector, 17 Worked Example of Numerically Integrated Truss Element ‘To choose an appropriate rule for evalusting (7.21) we need to determine the order of the highest degree polynomial that occurs inside the summation sign, For the sifiness ‘matrick, the terms in the matrix B are one order lower than the order of the expansion that is used to model the element displacements. In the case ofthe linear truss element, the B ‘matritis made up of constant terms and thus the order of its highest degree polyno! is zero. The order of the Jacobian matrix, and hence its determinant, can be determined by inspection of equation (7.14). Substituting the derivatives of the shape functions we see that det Ji given by detd = ee x)Jeosa + 4 ~yy)sina Noting that cosa = (x, ~ x))/Land sina +» — Y)/L, this ean be simplified to L(costa + sine) = & det = E(costa + sinta) = & (723) toprove that det Jis constant throughoutthe clement. Thus the highest degree term inthe ‘element stifiness matri, for a linear trusselement with constant ares and constant Young's ‘modulus, is equal to the product of terms which are all constant. The highest degree term6 inside the summation for kis therefore also a constant and, with reference to Table 4.1, it is necessary to use only a one-point Gaus rule. ‘To perform the numerical integration by hand, we frst note that from (7.12) fewsa -sina cosa sina] 24) Suistiuting (7.23) and (7.24) into (7.21), and using the one-point rule in Table 4.1, we Obtain the element stiffness matrix as sine foot pone] x & x feos Jone] § x2 costa cosasina -cos*a -cosasina] over -seasoa cata couine| 09 ‘Asthe clement stifinest is again symmetric and positive definite this implies thatthe global stiiness matrix will also be symmetric and positive definite. 76 A Worked Example of Truss Analysis using Finite Elements Consider the truss shown in Figure 72, which has three nodes and three elements and is subject toa vertical downward force at node two. All elements have Young's modulus E and cross-sectional area.A. Elements one and two are of length /2L, whilst element three is of length 22. clement and nodal freedom numbers U ‘Figure 7.2: Truss geometry and nodal freedom numbers Asin previous examples, the topology ofthe grid is specified by listing the nodes which ‘define each element. These nodes ae usualy specified in aconsstent manner, but foreachtruss clement the nodes can be specified in any order. The nodes for each element aze stored in a single row ofthe array ELTOP as shown below [ELTOP(1,1) ELTOP(1,2)] [i 2]+-nodes for element 1 ELTOP = |ELTOPQ,1) ELTOPQ,2)| = |3 2|—nodes for element 2 [ELTOP(,1) ELTOPG,2)| [3 1]~nodes for element 3 "The nodal freedom NF is constructed using the nodal freedom numbers shown. As in previous examples, the entries inthe nodal freedom array are formed by running through ‘the nodes in ascending sequence and allocating numbers to each degree of freedom that isunrestrained. By convention, theu-displacement isnumbered before the -displacement at each node and all fixed degrees of freedom are given a value of 2er0, For the mesh of linear rod elements shovin in Figure 7.2, the nodal freedom array is given wy [NFL 1) NFQ,2)] [0 oO} node 1 NF = |NF.1) NF@,2)] = [1 2]— node 2 INFG.1) NFG,2)] [3 o]— node 3 [Note that after the boundary conditions have been applied, there are only three unknown displacements, This implies the global stifiness matrix will be 3x3 and the global force vector will be 331. Because the load is applied to the truss directly at a node, the global {force vector can be assembled without forming element force vectors. Indeed, i this case te element force vectors can only be found after the solution forthe displacements has been obtained. For our tru problem, the point load P is applied downward at nodal ao: os ‘The nodal freedom array NF holds the nodal freedom umbers and, in conjunction with the element node numbers listed in the element topology array ELTOF, i used to form @ Steering vector for each element. Allowing forthe fact that presribed degrees of freedom are never assembled, te stcering vector tells us where the element stiffness matrixentries should be inserted inthe global stiness matrix, nodal freedom numbers Figure 7.3: Local and global freedom numbers for truss element one‘The local freedom numbers and nodal freedom numbers fo element one aze as shown in Figure 7.3. Note that local freedoms one and two are associated withthe node thats listed first in the element topology array (in tis case node one). The loal freedoms three and four are then associated with the other node. ‘This order also determines the value of the angle @ which, by definition, it measured potitive anticlockwise from the x-axis about the first node. Using equation (7.25), te stiffness matrix for element one is given by ° (7.27) ‘With reference to Figure 73, the steering vector for element one is, STHER! « [NFQ1)NFQ,2)_NFQ1) NFQQ2) Which indicates that ‘+ rowone and cokumn one of the element stiffness matrix are not assembled into the slobal stifiness matrix ‘+ row two and column two of the element stiffness matrix are not assembled into the slobal stiffness matrix ‘+E jisassembled into K;, of the global sifinessmatrisince its row and column indices fom the steering vector are (1,1) + Kisassembied into X,zof the global stifinessmatrixsince its row and column indices fom the steering vector are (1,2) + Hhisassembled into X,, of the global sitfness matrixsince its row and column indices from the steering vector are (21) ‘+ isassembled into Xz ofthe global sifiness matrixsine its row and column indices ‘rom the steering vector are (2.2) “fe ol: whe 1 oo; toca freedom numbers nodal freedom numbers ‘Figure 74: Local and global freedom numibers for truss element two For clement two, the local freedom numbers and nodal freedom numbers are shown in Figure 74, Thus its element stifness matrix and steering vector aren raay raf 1 14 Wagar 114 728) tat {NFG.1)_-NFG,2) NFQ,1) NFQ,2)} STEER’ B 0 1 a The steering vector indicates that ‘+ row two and column two of the element stifiness matrix are not assembled into the ‘lobalstifiness matrix + jisassembled into Kas ofthe plobalstifinessmatrixsinceisrow and columningices from the steering vector are 33) + isassembled into Ky ofthe global stiiness matrixsinceitsrow and column indices from the steering veetor are (3,1) + Hgisessembledinto Kp ofthe pobalstiffessmetrix since itsrow and column indices srom the steering vector are (3,2) +B isessembledinto X,,ofthe global stifiness matrix since itsrowandcolumnindices from the steering vector are (1,3) + Hyisassembled into K;, ofthe global stiness matrixsince its ow andcolumn indices from the steering vector ae (11) + Hsassembled into X,.ofthe global stifiness matrixsineitsrowand column indices from the steering vector are (1.2) + yisassembledinto Kx ofthe global sifiess matrix since itsrowandcofumn indices fom the stering vector are (23) + Kis assembled into Ky ofthe global stifinessmatrxsinceitsrow and column indices from the stering vector are (2,1) + HLjisassembled into Kz ofthe global stitiness matrixsinceitsowand column indices from the steering vector are (22) a= 180" cosa sina = 0 ays 2 vito vyh0 3 1 9 3 o4 @o«4 ou @ % local freedom numbers nodal freedom numbers Figure 75: Local and global freedom numbers for truss element three For element three, the local freedom numbers and nodal freedom numbers are shown in Figure 7.5. The corresponding element sifiness matrix and steering vector aren fo -a g 00 00 . a9 2hr|-6 0 A o ) 00 00 STEER? = [VFG,1) NFG,2) NFC.3) NFUL2]}= 0 0 9) “The steering vectorindcstes that ‘row two and column two of the element sitfness matrix are not assembled into the ‘lobal stiftess matrix + rowthree and column shree of the element stifiness matrix arenot assembled into the slobal stiffness matrix ‘+ row four and column four ofthe element stifiness matrix are not assembled into the slobal stiffness matrix ‘+ Kj jisassembled into Xs ofthe global stiffness matrixsince itsrow andcolumn indices fom the steering vector are (33) ‘Adding the stiffness conteibutions from each element gives the global stiffness matrix According to K= we + we + we 19 tact 0 0 Sel 1 of Bt]s 1 ifs tefo o | HL oof Ls a a] lo o a 200 4 -felo 2 1 030) 11 B41 ‘Combining (7:30) with the global force vector, given by (7.26),we obtain the global stiffness equations as| nf) 9 Alo 2 ipral = |p| 22r| : oa 1 1 Bally} Lo) Solving this system furnishes the unknown nodal displacements as, 1pL Vz = (i+ 1/2) BL u-%Together with the known boundary conditions U, specifies the solution for the displacements. as this completely (Once the displacements are known, the strainsin each trusselement can be computedusing the strain-displacement relations (7.11). In element one, for example, and Busted PL EA Wary = (E+ yee this implies thatthe axial strsin for element one is Similarly for elements two and three (ox 0 AE: a 3 ~ 20 Once the strains are known, the axial stresses in the truss can be computed from the stress-strain relation (7.15) according to = BE eomresin = i™ (compression) (734) tension) (738) Fix old = - 2p (compression) 738) (compression) 37 Fheold= 4P (tenon) 038) ‘These forces canbe verified by checking equilibrium ateach nodeinthe truss, Atnode two, for example, we assume the forces shown in Figure 7.6. P @ <> FP A Figure 7.6: Equilibrium for node two in truss From equilibrium Sh = ~ Fleoses + Foss YB = = Flsings ~ Fleins = P ‘These equations give the axial forces in elements one end two as -P- (compression) :P_ (compression) ‘hich are identical to those obtained from the finite element analysis,BEAM ELEMENT CHAPTER 816 8 BEAM ELEMENT In the analysis of beam deformations, two different types of behaviour need to be distinguished. For slender beams we negleet shea effects and assume thatthe plane which isinitialy normal othe neutral ais remains normal o this axis after deformation has taken place. This assumption islustrated in Figure 8. Figure 8.1: Beam with no shear deformation In this discussion v denotes the vertical movement of ofthe beam (positive upwards) and 8 denotes the rotation of the normal to the neutral axis (positive anticlockwise about yo). For the case of thick beams, shear deformations cause warping across the beam section, ‘To model this analytically we still assume tha plane sections remain plane during losding, ‘but allow the plane normal to the neutral axis to rotate as shown in Figure 82. neutral axis Figure 82: Beam with shear deformationn ‘The rotation ¢ defines a constant shearing strain through the depth ofthe beam, Since the actual shearing stress and strain are not constant (they vary parabolically throughout a rectangular cross-section), ¢ is more precisely defined as an equivalent shearing strain acting on ¢ shear area 4y. ‘Thus we can write 1) ‘where ris the shear stress, G = £/2(1 + »)is the shear modulus, Q isthe shear force, A is the cross-ectional area, A, is the effective shear ares, and k= A,/A ~ 5/6 for ‘rectangular sections. From elementary beam theory, the moment at any section is given by 8, m= ard 62) ‘where isthe moment of inertia and £/is the flexural stiffness. Equations (8.1) and (8.2) define the “strese-sirain” relations for the beam and may be combined in the matrix forma {a [> cal a ‘This can also be written inthe more general form onDe (84) ‘where the generalised “strain” and “stress” vectors are given by 65) 66) HO Plo Gak en Note that in this cas, the stress vector is actually comprised of stress resultants (the ‘moment and shear force) whilst the strain vectoriscomprised ofa curvature andan average shear strain. For this reason these vectors are, respectively, often referred to as the generalised stress vector and the generalised strain vector. (Se (oe) Hl ][o5 postive shear force negative shear force and the “stress-strain” matrix is Figure 83: Sign conventions for moment and shear force% ‘The sign conventions forthe moments and shear forces are shown in Figure 83. 8.1 A Linear Beam Element with Shear Deformation ‘To derive the stiffness equations for a general beam element, ve allow forthe effect of shear deformations and use the principle of virtual work. Consider the two-noded linear ‘beam clement shown in Figure 8.4 where each node has two degrees of freedom, a vertical isplacement v (positive upwards) anda rotation 6 (postive anticlockwise) yy Figure 8.4: Two-noded beam with shear deformation ‘The linear shape functions sre expressed in terms of the model coordinate £ according to N= -E-1) 1 m- de+n Once their nodal values are known, the displacements and rotations at any point inside a ‘beam element can be found from Nyy + Non = [Mi afr} 68) A we} 69) ‘These equations can be expressed more compactly in the matrix form of| Y) fm ox, oe lof “] ON, 0 NA] y¥2p = NO (810) e: = NO; + N30y where N, ON, 0 N=low, ow, fr % v2 &} (au)Bp In order to compute the sifiness matrix for the beam element itis necessary to derive the strain-displacement matrix B. This matrix relates the strain vector at any point inside the ‘beam, as defined by equation (8.5), to the nodal displacement vector of equation (8.11), Differentiating (8), we obtain (0; It} ee Differentiating (8.8) and using (8.9), this quantity may also be expressed as o-[ f} -[K ofr as) Aisplacement relations as ay, ("4 “ |I0, Ny lg ce le, Bu eu) Where the strain vector eis defined in (8.5), the element displacement vector wis defined in (8.11) and the strain-displacement matrix Bis an, aN Cao =e Bel, (15) ae For an isoparametric formulation, the shape functions that define the displacements are also used to define the coordinates. Thus thox-coordinate at any point within the clement can be written in similar form according to TEN, +My, = [My wf} 616) Differentiating (8.16) with respect to the model coordinate & gives the determinant of the Jacobian as80 ‘To compute the B matrix, we require the partial derivatives of the shape functions with ‘espectto the cartesian coordinates andy, These derivatives cannotbe evaluated directly 8s the shape functions are written in terms ofthe model coordinates and 7, s0 we invoke the chain rule according ‘The partial derivatives of the shape functions with respect to the model coordinates, together with the det J, enable the required values to be found. 8.2 Formulation of Linear Beam Element Using Virtual Work For the beam clement, the virtual work equation may be written as 4 (elf +560) a= [ances my +P C0, +6380, (18 Where A denotes the cross-sectional are, 8 2is the vistwal curvature, Mis the moment, 8 he vital equivalent ser sai, Q's the sea oc, sth ital verte aiplacenent te dstrutedou por ot eng (ry) ae he ial verte dipdements athe nodes, (P,P, he vera! pon loud applied athe nde (28; .68,) are the virtual rotationsat thenodes, and (C, ,;) are the couples applied a the nodes. a ee affsteu [ars aw where for 40, 40 604) Now, from (810) and (8.4) it follows that eat outNtaL oe = ouTBT ‘Substituting these equations in (8.19), together withthe stress-strain relations (8.4), and rearranging gives, ota womae sare = af NTeaerae eo a0 Ly Since the variation Buin the nodal displacements is arbitrary, and (8.20) must hol forall such variations it follows that a ‘This expression defines the clement sifiness relations winced w= f Ntecouae +P ka=t where the element stiffness matrix may be identified as ke 4f BtDp det 34f (62 and the element force vertor is [ wrccuae oe 2 The nodal displacement vector for the beam element is given by (811). 8.3 Numerical Integration of Element Stiffness Matrix ‘To derive the stiffness equations for a numericallyintegrated element, we substitute 8.15), (827) and (8.17) into equation (8.21) and apply numerical integration. This gives k= AS BIDS, det Jp, ey) ‘where nis the number of integration points, w; the value ofthe weigh for integration Point i, B, = B(,) is the strain-isplacement matrix evaluated at integration point i, et J, = 4etJ,) is the determinant of the Jacobian evaluated at integration point i, end is the value ofthe model coordinate for integration point 8.4 Formulation of Element Force Veetor In beam analysis, the loads may be applied ina distributed manner or as concentrated forces/couples atthe nodes. These types of loading are shown, respectively, in Figure 85 and Figure 86. Note thatthe positive direction ofthese actions is as shown. Since the ‘wo-noded beam element uses linear shape functions to model the displacements and rotations, its posible to compute consistent nodal forces fr a lineal varying distributed2 lod. For problems involving more complicated distributed load, it would be necessary to.se an element of higher order yy a) Figure 86: Beam clement with forces and couples applied at nodes In the most general case the element force vector comprises both forces and moments. If the mesh ofbeam elements islozded onlybypoint forces and couples at the nodes, as show in Figure 8.6, then the global force vector can be assembled ditecly (inthe same way a3 for the truss element discussed in Chapter 7). For example, consider the mesh of two linear beam elements shown in Figure 87. Py a Ne, aa ae ® ® ® Figure 8.7: Mesh of beam elements with forces and couples applied at nodes Ignoring boundary con: (624) For the case of distributed loading, the first term on the let hand side of (622) may be integrated numerically using the Gauss rule to give the element force vector as DyNT det (625)where isthe numberof integration points, wis the value of the weight fr integration pointi, N, = NG) isthe shape funtion matrix evaluated a integration point isthe ttaction vector evaluated atintegraton point det J, = detJ(,)'s the determinant ofthe Jacobian evaluated at integration point and §, the value ofthe model coordinate for integration point: Once each element force vector i formed, itis then inserted into the slobal force vector using the element tering Veto. 8.5 Worked Example of Numerically Integrated Beam Element To illustrate the steps involved in forming the stiffness relations by numerical integration, 2 worked example is presented for the twomnoded beam element. In performing the numerical integration by hand, we first note that aNidé _ ad (1 rs he - 8-4 ag dt Ba 00+ 9)e- 3 ago t Sain (826) Using (8.15) ths gives the strain-displacement matrix as 4 1 a | ae ma 1 1 Hi -, by THE all 2 j thy ow -n,| ‘To choose an appropriate rule for evaluating (8.28, itis necessary to determine the highest order polynomial that occurs inside the summation sign. Inthe ese ofthe linear clement, the B matrix comprises some linear terms and some constant terms and thus corresponds toa polynomial of degree one. Provided the beam element isnot curved, the Jacobian matrix, and hence its determinant, is always constant. Thus we see that the highest 4egrec term in the element stiffness matrix is equal to the product of terms which are of ‘degree one. The highest order term inside the summation fork is therefore of degree two and itis necessary to use a two-point Gauss rule. For , equation (8.28) becomesry TH Oph e x -eme where, from Table 41, the coordinates and weights are Ae UB w= ey wat Substituting ME) = M0) 3) MED) = NCU) NG) = NAGI) = Nab) = NCA) and the integration point weights gives Gak As expected, the element stiffness matrix is symmetric, For the element force vector, consider the case shown in Figure 88 where the beam clement i subject to linearly varying distributed load defined by 9 = Na, + Node (830)Figure 88: Beam clement with linear distributed loading Using equation (8.2) the element force vectors given by [MGd 0 Shey fee") ae 0 NG) PED) 0 Inseeg @@2] * % xm, ° “8 rs Dd |wstep [Ga + Mae d93} x E xm, aM sG0(uEda + NAED9R) o 2 WENGE Dan + NaEi)aa) o Inside the summation sign, the highest order terms are products of the shape functions. Since these functions are of degree one, thisimplies thatthe highest order term is quadratic ‘and a two-point rue is again sufficient. Using the two-point rule scheme in Table 4.1 we obtain fegp(mche + wide) ° t= dy (x Avg, + nxAy Nat (side weds) 0 xke ie8% (Ay weg, + nc4 IN (me Bn + NL +p) 0 xh ) + Which simplifies to es + do a 5lhe +e 5 tek ‘As an immite check on this lement force vector we nts hat Sonoda foves = 4(qy + 43) = toa load apple to beam ‘Asaspecal cat, the lament force vector fora unifrm dtbuted load cane obtained tysuisinaingg, = 42 = atone ft alo tohCHAPTER 9 TWO-DIMENSIONAL INTERPOLATION AND INTEGRATION9 TWO-DIMENSIONAL INTERPOLATION AND INTEGRATION Lagrangian interpolation is not only applicable in one-dimension, it ean also’ be extended to two- and three-dimensions. Two- and three-dimensional finite clements are ‘used widely in engineering analysis to predict the behaviour of plates and other continuum problems. This Chapter describes how to derive the shape functions for a number of ‘two-dimensional triangular and quadrilateral finite elements. Indeed, using the method described, it is possible to derive the shape functions for any two- or three-dimensional finite element. 9.1 Linear Interpolation Over a Triangle ‘The linear triangle was one of the earliest clements used in two-dimensional applications. Although higher-order formulations, such as the quadratic triangle, are generally preferable to this element, itis still widely used in practice. The linear triangle hnas three nodes, and when used for analysis of linear elastic solids, has two degrees of freedom at each node, These degrees of freedom correspond to the displacements in the ‘x-and yediections as shown in Figure 9.1, 1 on ® mapping ye @ bd ® (0) a0 ao ¢ Try ttl coordinates ‘model coordinates Figure 9.1: The linear triangle ‘The mapping for triangular elements is based on the unit triangle, which has model coordinates (7) and vertices located at (1,0), (00) and (Q,1). To ensure that the determinant of the Jacobian is positive, the nodes for each clement are always listed {anticlockwise (this point will be discussed further in a later Section). For an isoparametiic formulation, the same shape functions are used to describe both the coordinates and the 4isplacements throughout each element. The coordinates a any point inside the clement are given by Noxy + Nyx3 = Nx, (91) Y= My + Nave + Nays (92)‘where the shape functions Aare linear functions ofthe model coordinates (67). Similarly, the displacements are defined by w= Nye, + Ney + Nyy = Nay 3) YN + Navy + Nyy = Nyy 4) ‘The precte form of the shape functions canbe found using the same technique that was ‘sed forthe one-dimensional case. ‘Atmode one, where (1) = (Bj) = (1.0), We must have Ey) = NyG pm) + NaE yada + Naya) = uy and hence Mn) =1 es) Navn) = 0 06) NaGyn) 0 on Similarly for node two, where (€.9) = 73) = 0,1), we must have Eats) = NyGortadty + Ney tah + NaEona)ty = ty hich gives NiExm) = 0 (98) NxM) = (94) Ns) = 0 03) Atnode three, Gm) = (75) = (0.0) and HEa9) = NyCorta)y + NgCons Ma + N3G Therefore MGon) =0 20) NyGyn3) = 0 (94) Noy) #1 ean ‘The conditions imposed by equations (9.5)~(9.11) may be concisely summarised as NGn) = NG) where (f.m) are the model coordinates for node of an element with m nodes. for j =i (9.12) for j= tj.m and j ei 13)0 ‘Toderive the shape function Ny, we fits note thatitmust be 2610 at nodes two and three, ‘This condition follows irom equation (9.13) and may be writen as Ni Gy) = N(Gyns) = 0. Since the equation oftheline passing theoughnodes to and threes = 0, thete requirements are automaticaly satisfied by guessing that, ofthe form ™ ag) Where isan unknown multiplier. Substtting the condition that Ny(Ey.m)) = ty which follows from equation (9.12), this implies that 1-4) asl and henee ™ ‘The derivation of N follows a similar pattem. Since the equation ofthe line passing through nodes one and three is = 0, the conditions (G41) = Ny(Eyns) automatically satisfied by guessing that Nis ofthe form N= 4) Substituting the requirement that Na(@,) sives a and hence Maa Finally, for NV, we note thatthe equation of the line passing through nodes one and two is &+—1= 0 It follows that the conditions N,E,m,) = NyGqn,) = 0 are ‘automatically satisfied by guessing that Nis ofthe form (+n-1) Substituting the requirement that Ny(b53) = 1 gives & Nya1-b=9 Jn summary, the shape functions for the three-noded triangle are Moe (24) Men 13) Nye 1-Eon (016) ‘These functions may be checked to ensure that they satisfy all of the required conditions, Atnode one where (fm) = (1,0) y,a At node two where (&n) = (0,1) At node three where (67) = (0,0) Another useful check on the shape functions is that they must sum to unity according to Mam tm= Suet a ‘As discase in Chaptr2, thi condition follows from the ned forthe element tobe able ‘to mode! a uniform displacement field. Summing (9.14, (9.15) and (9.16)-we see that this requirement is satisfied 9.2 Quadratic Interpotat Over a Triangle ‘The quadratictriangle for continuum analysis, shown in Figure 92, has hassixnodes and twelve degrees of freedom. Unlike the linear element described inthe previous section, ‘the quadratic element can have curved sides and is thus suited to madelling domains with curved boundaries. ” ®. uy (0,1)! %y rapping © ager @ 12,12) i < L (0,172) ean © @ OO: ny @ (00) 20) G0) real coordinates ‘model coordinates igure 92: The quadratic triangle ‘The coordinates at any point inside the element are given by2 Nyry + Nyey + Nyy + Neg + Nets + Net Sun an YN + Nava + Nays * Nave + Nov + Ney = DN (0.18) Where the shape functions Nj are quatiratic functions of the model coordinates (7). ‘Similarly, the displacements are defined by. We Nay Nats Nyy tN Nyy + Neus = Say 009) =t PEN EN Ey + Neve Noo tei SN (020) ‘The precise form of the shape functions can again be found from first principles using the ‘same arguments that were developed for the linear triangle. ‘Tocompute N;,wenote thatthe equations ofthe two lines which passthrough ll ofthe nodes except node one are & = Oand £ ~ 1/2 = 0, Ifwe guess that My = AGE 1/2) ‘his cutomatically ensures that 1 is zero at all nodes except node one and hence satisfies ‘equation 9.13). The multiplier is found using equation (9.12) which states that MGon) = AIG, ~ 172) =1 ‘Substituting (€).9;) = (1,0) then gives ACYL = 1/2) = 1 and and hence B= 26NE - 1/2 omy ‘Te deduce Np, we note that the equations of the lines that pass through all nodes except ‘node two are § = Oand y = 0. A suitable guess for N', which satisfies (9.13) is thus Na = 260) Using (0.12), Ny must aso satisfy NiGom) = 469 =1 where Ea) = (0/2, 1/2). This implies AQ/2\/2)= 1Which gives A = 4 and §= 4G) 2) ‘Similar calculations can be performed to yield the remaining shape functions as y= ayy (23) 4@ye #90 2a Ns= UE+n-12)E 49-1 925) 4G\E +1 — 1) (26> 9.3 The Jacobian Matrix for Triangles ‘To compute the element sifiness matrix in terms ofthe model coordinates & and nit is necessary to define the Jacobian matrix for two dimensions. The two-dimensional Jacobian matrix is fax 2) a ae -|\n 27 ae The determinant and inverse of J follow as ad _ me seu = BEB ea) x» my 3 wala 3 (029) an aE Consider the Jacobian matrix forthe three-noded triangle, Differentiating equations (9.1) and (9.2) gives z aXe (030) z (931) 2 eo 2 ro) Further diferentiation of equations (9.14), (9.15) and (9.16) yields the shape function derivatives asBo 0 ay mh Meo Me, e Substituting these derivatives into equations (9.30)+(9.33), we obtain the terms of the Jacobian matrix according to Bean-y (036) Benen 3 Benny (38) Hrs (39) Inserting these values in (9.28) gives the important result that St = Gy — 202-9) — G29) = 2A (04g Where 4 isthe element area. This result is valid for any triangular element with straight sides and gives a physical interpretation of the Jacobian matrix. It also serves as a useful data check in compater codes since the arce ofeach element, and hence det J, must always be positive If this isnot the case, then there isan errr in the nodal coordinate data or the list of element nodes. Unfortunately this simple geometrical interpretation does aot hold {or triangular clementswith curvedsides,since the determinant ofthe Jacobian no longer constant throughout the clement. The condition that det J is positive, however, is always neoessary, 9.4 Numerical Integration Formulae for Tiangles ‘Tocompute the stiffness matrifor a triangular clement, it isnecessary to evaluate integrals of the form Hi Jf penser [ [ ands (a1) ‘Numerical integration rules for triangles replace these integrals by summations according to 1 pee fee I J J 1G.mdands = Y FE,.ndm, (9.42) where » is the number of integration (or Gauss) points, (fn) are the coordinates for imtegration point i, and w, i the weight for integration point i The sampling points and eights for various Gauss formulae are shown in Table 91. ‘As in the one-dimensional case, the rules are designed for cases where f(y) is @ polynomial. With reference to Table 91, we see that a one-point rule gives no ervesfoy e95 constant or linear function, whereas a three-point rule gives no error for @ polynomial up ‘oand including a quadratic. Other formulae for higher order polynomial integrands are ao available and are given in many numerical analysis texts, Like their one-dimensional relatives, the triangular integration rales be used for cases where the integrand is not a polynomial, but the result is no longer exact. “We. poins coordinates weight degree of” ” & u polynomial Z i B wm r zB 175 76 3 6 28 Ys 6 1/6 us ‘0.65902 76223 74002 028193 33685 5301/12 0.65902 76223 74002 010503 90080 72877 1/12 0.23193 33685 53031 0.65902 76223 14052 1/12 0.23193 33685 53031 0.10803 90090 7287 1/12 0.10903 90090 72877 0.65902 76223 14092 1/12 0.10903 90000 72877028193 33685 53081 1/12 ‘Table 9.1 Data for Gaussian integration over triangles ‘illustrate the use of these rules, consider the integral ef [Perens mane ‘Evaluating this exactly gives raf @- 3 +84 Nag = 3/8 ‘Using the one-point rule fom Table 9. gives 1= fmm G+ 8m + adm, Letxte bl GAgX5t 92 7/8 ‘This isnot exact because we need a quadratic rule to avoid any error. For the three-point rule we obtain | L=fGs.my)vy + FCarmve + FEs.nsMs=3/8 As expected this results exact as we have used a rle with quadratic accuracy to integrate ‘ quadratic integrand, ‘95 Linear Interpolation Over a Quadrilateral ‘The linear quadrilateral has four nodes, and when used for analysis of linear elastic solids, has two degrees of freedom at each node, This element and its degrees of freedom are shown in Figure 93, val ,29® Gn G1) Zu eal coordinates ‘model coordinates Figure 9.3: The linear quadbilateral ‘The mapping for quailateral clementsisbased on asguare whose side length two units ‘The four vertices have the model coordinates of (1,1), (—1,1),(—1,—1)and (1,~1) and the nodes are again listed anticlockwise to ensure that the determinant ofthe Jacobian i Positive. The model coordinates are chesen 10 range from ~1 to +1 a8 this proves ‘onvenientin applying the one-dimensional Gaus integration formula to twordimensionsl problems. For an isoparametric formulation, the coordinates and the displacements throughout each element are given by Nye + Naty + Nya + Nate = SN (043) tes Y= Ma + Naya + Nos + Nag = D ND (9.48) = Nyy + Ngly + Nyy + Nyy = 5 Nit & = Nv, + Nave + Nays + Navy = Ney a ‘The precise form ofthe shape functions can again be found from first principles using the ‘conditions imposed by equations (9.12) and (9.13).7 ‘Tocompute Nj, we note that the equations ofthe two lines which passthrough a of the nodes except node one are € + 1 = Oand + 1 = 0, we guess that Ny =46+ D+) {his eutomatically ensures that Nis zero at all nodes except node one and hence satisfies equation (9.13). The multiplier is found using equation (9.12) which states that NyGyn) = AG, + Dn + D Substituting @,,,) = (1,1) then gives a+ a+ 1-1 and hence Nya 4E+ DG +0 (045) Applying an identical procedure atthe remaining nodes gives 46 -Da+1 (9.46) 3e6-na-0 (47) -}e+D@-1 (9.48) 9.6 Quadratic Interpolation Over a Quadriateral ‘The most common form of quadratic quadrilateral has eight nodes, and when used for analysis of near elastic solids, has to degrees of freedom at each node. This clement and its degrees of freedom are shown in Figure 9.4, Asin the case ofthe quadratic triangle, the quadratic quadrilateral can have curved sides, (4 Od) Gy) QoQ @ oO © —- (-19) Go) o @ Ci) @=1) G1) x,u real coordinates ‘model coordinates Figure 94: The quadratic quadrilateral8 For nisoparamettic formulation, the coordinates and the displacements throughout each element are given by x= DM y= My ie w= SN, ee v= Ny, & ‘To compate 1, we note that the equations ofthe tree lies which passthrough all of thenodes excepinodeoneare +1= 0,7 +1=Oandé +n ~ 10, Kweguessthat Ny = AG + 1M + WE ty - 1) this automatically ensures that Nis zero at al odes except node one and hence satisfies ‘equation (9.13). The multiplier is found using equation (9.12) which stipulates that NiGyn) AG. + Dlr + DE tm ~ Y= 1 Substituting (&).m)) = (1,1) then gives AQ ++ Hesi-yer or aes and hence M=4@4+ 00+ DE +n-1) Applying an identical procedure atthe remaining nodes gives N= RE +E -DO+D Ny = -HE- Do + Dey -E~1) Ne= 4E-Do+ Do- 0 4E- Ne ne + +1) Now E+ INE ~ I - Ny= -}E+@- WE =7-1)8 Ny=-hE+D0 + D@-D 9:7 The Jacobian Mates for Quadrilaterals Consider the Jacobian matrix forthe four-noded quadrilateral, Differentiating equations (943) and (9:44) gives z Stan 49) # 050) 2 es) 4 Sa,, (952) Further differentiation of equations (9.48), (9.46), (9.47) and (9.48) yields the shape function derivatives as wy, er dot 053) 58) 55) (56) Substituting these derivatives into equations (9.49)-(9.52), we obtain the terms of the Jacobian matrix according to Baje- a -y tet }e- ate sn Balet ander ntn = Be dor- monty +EO- yt —yon (059) Be hort ns -70 + 301- vy ts 208 (060) Inserting these values in (9.28) we see that, for the general case of an arbitrary quadrilateral, detJ is a polynomial of degree two. This implies thatthe value of the Jacobian determinant varies from point to point inside the quadrilateral, unlike the ‘mapping forthe linear triangle. There i, however, ane important special case where the100 detJ becomes constant and a simple physical interpretation of its meaning becomes possible. This is when the sides of the element form 4 rectangle (or a parellelogram) as shown in Figure 9.5. igure 9.5: The linear rectangle For thiscase, we see that yomea Hom =0 YrwaO yeh une and hence =}O+fOE=<0 1 ©+}On=0 gig ale Se ae +t e 8 50H) +508 = 4 yom a0 nomad “Inserting these values in (928) gives the impostant result that detF = (@/2N(6/2) — (0)(0) = 4/4 where isthe element area. This results valid for any quadrilateral element whose sides form a rectangle or a parallelogram, including the higher order elements, For (quadrilateral of arbitrary shape, the condition that det J is positive must stil be satisfied. (61) (9.62) (983) st) (965)1 9.8 Numerical Integration Formulae for Quadrilaterals T compute the stiffness matrix for « quadrilateral element, itis necessary t0 evaluate Sntegrals of the form aH I J f SG n)ddy [ SG.n)dé (0.56) ‘Numerical integration rules for quadrilaterals are similar to those for triangles andreplace these integrals by summations according to Ts { | HG natin = 7G, (sn where nis the numberof integration (or Gauss) points, (Fy 1) ate the coordinates for integration point nd wi the weight for ntgraton pola The sampling points and weighs fr vorious Gaus formalae over a quailatra are shown in Table 9. “These rules area straightforward generalisation of the one-dimensional case and are thus designed for integrands where (E,7) is polynomial. With reference to Table 9.2, we see ‘that a one-point rule gives no error fora constant or linear function, whereas 2 four-point rule gives no error for apolynomial up to and including acubie. Other formulae for higher ‘order polynomial integrands may be constructed by applying the one-dimensional rules in two-dimensions. opens ——~—~S~Crdates ‘welght degree of n & 1 w;__ polynomial Z 0 v 2 T ‘Tle 9.2 Data for Gaussian integration over quedrilaterals102 ‘Consider the integral 4 ia + bn + WP hdbdy = 8/3 ‘Applying the one-point rule from Table 92 gives T= FG nds = P+ an + ape = (0+0x0+02 =o ‘This's not exact because we need a quadratic rule to avoid any error. For the four-point rule we obtain T= fEombea + FG anda +f Esotadvs + FEtedme = [BF + A C8) + BY) + [Ay + ey Aas) + a/ By + [as BY + a/ Bye 8) + CB) + [ory + a Bay 8) + a/8Ja) = 5/5 As expected, this result is exact as we have used a rule with cubic accuracy to integrate a (quadratic integrand.CHAPTER 10 PLANE STRESS ELEMENTS108 10 PLANE STRESS ELEMENTS. ‘This chapter illustrates the development of two-dimensional elements for analysis of linear elastic solids under conditions of plane strain. The governing equations are derived {rom the principle of virtual work and the consistent element force vectors are established for a vatity of loading conditions. Inthree dimensions, the stateof stress at pointis defined by the nine stress components shown in Figure 10.1. Tensile normal stresses are taken as positive, The first subscript on ‘the shear stress denotes the axis which is normal to the shear stress plane whilst the second subscript denotes the coordinate direction in which the shear stress acts. For faces where the normal tensile stress is in a positive coordinate direction, the shear stesss positive if it points in a postive coordinate direction. For the faces where the normal tensile tress 's ia a negative coordinate direction, the shear stress is positive if it points in a negative coordinate direction. The positive senses fr the three normal stresses andthe six shear stresses are shown in Figure 10.1 Figure 10.1: Three-dimensional stress state By taking moments of the various forces acting on each face of Figure 10.1, it may be shown that and thus the three-dimensional stress vector comprises three normal stresses and three shear stresses. For anisotropic linear elastic solid, the strains ae given by Hooke's law according 1 | [ox - (ay + 04)) [o, - Ho. + 00) qoay (10.2)10s Defining the Lamé constant 2 T= (103) ‘equations (10.1) may be inverted to give the stresses in terms of the strains ae G+ 2G) +A + 02) ty = Gy 0) = C+ 2G + Ales +e) be = Ge ana) 05 = G+ 2G) + Hex + 6) ta = Ga ‘These relations define the elastic constitutive law which sof the general form o=De 05) where OT = {0 05 05 Tay Te Tord aT = (en by fe Yay Ye Yor) and the stress-strain matrix is R+2G 4 4 0 0 0 p +m i 0 0 0 p-| 2 2 442G 0 0 0 0 0 o Goo 0 o o o Go 0 0 o 0 0G 10.1 Stress-Strain Relations for Plane Stress Plane stress loading occurs commonly inthe analysis of beams and thin plate, where it is assumed the stress in the out-of-plane direction 9, is zero. The three nonzero stress ‘components 0,0) and are shovinin Figure 10.2. Thebodyis assumed tobe of thickness ‘hin the 2-dtection } ‘Figure 10.2: Plane stress loading105 ‘Under conditions of plane stress, we have o;= rye = ta =O and the {hreedimensional stress and strain vectors reduce to oT = {0,05 ty) (10.6) {e & Yy) 07) Defining 2G 477406 the corresponcing stress-strain matrix becomes 4G A oO p-| 4 A426 0 08) 0 0 4 [Note that the condition o, = 0 implies thatthe out-of-plane strain can be expressed in terms of the remaining strains according to a= e ta) 102 Plane Stress Linear Triangle ‘To derive the stiffness equations forthe plane stress linear triangle we again use the principle of virtual work. The element has three nodes and two degrees of freedom ateach node. These degrees of freedom correspond tothe displacements in thex- andy-directions as shown in Figure 103, , ayy +O nays > 20a) Bs Oo a ® (0.0) 0) € Te esos dc Figure 10.3: The plane stress linear triangle ‘The interpolation and mapping for triangular elements has already been discussed in Chapter 9 where it was noted that, for an isoparametric formulation, the sume shape functions are used to describe both the coordinates and the displacements. ‘The displacements may be computed from107 Muy + Mate + Ny 209) VE Nyy + Navy + Nyy (20.10) where the shape functions are given by (2.14), (9.15) and (9.16) according to Ni =e (10.11) iw) [Ny 0 Ny 0 Ny 0] fut Where the shape function matrix N and the element nodal displacement vector w are Sefined as MON, 0; 0 0M, ON, ON; (aos) we fk ma oa sm) nas) Inorderto compute the stiffness matrix forthe triangular element itis necessary to derive the strain-displacement matrix B, This matrix relates the strain vector at any point inside the triangle, as defined by equation (10.7), to the nodal displacement vector of equation (10.16). Differentiating equations (10.9) and (10.10), we obtain an ao.a7) a= (20.18) ang) fan, ra Collecting terms into a single matrix equation gives the strain-displacement relations10 ee) [ig N ay +, Vf yar 8 at OO lly, ay, aw, anglf slelo 2h go Meo BM 3 y lye aN aN aN, Ny aN ANG] ey be) [ot at oe a Bly These are of the standard form e=Bu 029 ‘whore the strain vector eis defined in (10.7), the element displacement vector w is defined {in (10.16) and the strain-displacement matrix Bis [2M (10.21) a For an isoparametric formulation, the shape fonctions that define the displacements are also used to define the coordinates. Thus thex-andy-coordinates at any point within the clement can be written inthe form x Nyx + Nat + Nyt (10.22) y= Ny + Naya + Novo (20.23) As discussed in Chapter 9, the determinant of the Jacobian for two-dimensions is ao oe ag oy Differentiating (10.22) and (10.23) with respect to the model coordinates £ and » gives dees aoa) eM May ag ag aE Sie sie als and henceact 2025) ‘The strain-displacement matrix (10: 2 contains partial derivatives ofthe shape functions with respect tothe real coordinates andy. As the shape functions have been expressed in terms of the model coordinates § and 7, these derivatives cannot be found directly and must be computed using the chain rue of differentiation. This rae gives (any) ae fae | foe |} ae anit “Jae 2) 2M tar) [ala (20.28) aN where J is the two-dimensional Jacobian matrix defined in equation (927). Inverting (10.26), we can obtain the derivatives of the ith shape function with respect to x and ) according to oan where J? is the inverse of the Jacobian matrix defined by (9.29). These derivatives complete the formation ofthe B matrix. 103 Formolaon of Plane Stres Linear Tangle Using Virtual Work Fora plane stress element of thickness the viral work equation may be writen as af (Ge.04 + 8665 + bya = af (Gge + dug ]at + af (uk + oFea + Pdi Pybry + sae + Pygbts + Pyar0 where A denotes te element area, (5, de, 67) ae the virtual strains, (0, dy ty) are the stresses, (bu dy) are the virtual displacements, (g, q,)are the surface pressuresin thex and ‘yiections applied over the edge of length L, (X ) are the body forces in the x and ‘yediections, (6u,,;) are the virtual displacements at node i, and (P,P) are the point loads in the x and y-directions applied at node ‘This equation may be written in matex form as if setoas = af oarsan +f saTbaa + PTu (1028) where {Be 86) br} bat = fou 6 T= 9] ey du? = (5, 87, By bry buy drx} ‘Now, from (10.18) and (1020) it follows that dal = ou NT eT = ua? ‘These equations, together wth the planestressform ofthe constitutive relations (10.5) can be substituted in (10.28) to give ota wPonaer sad = oof Nieat cant Nica © do" (029) Since the variation dun the nodal displacements is arbitrary, and (10.29) must hold forall such variations, i follows that Af winacr rae w anf nea on Nthae sr ly Nh ‘Thus the element stifness relations may be written inthe standard form ust ‘where the element stiffness matrix and force vector arem kaaf smase oo + fuss wou se ans ‘and the nodal displacement vector is defined by (10.16) 104 Numerical Integration of Stiffness Matrix for Plane Stress Linear Triangle ‘To drive the stiffness equations for amumerically integrated element, we substitute (10.8), (20.21) and (10.25) into equation (10.30) and apply numerical integration. This gives x BIDS, det; (1032) ‘where n is the number of integration points, w;is the value ofthe weight for integration point i, B, = B(E,.n,) is the strain-isplacement matrix evaluated at integration point i, det J, = det 1G,.7, is the determinant ofthe Jacobian evaluated at integration point, and G7) are the values of the mode! coordinates for integration point i 10.5 Formulation of Force Vector for Plane Stress Linear Triangle In plane stress the loads may be applied as point forces atthe nodes, pressures on the clement edges, or distributed body forces over the element volume. The latter case is ‘common when gravity loading is imposed. These types of loading are shown, respectively, in Figure 10.4, Figure 1055 and Figure 10.8. Note that the positive direction of these actions is as shown. Since the three-noded triangle uses linear shape functions to model the displacements, itis possible to compute consistent nodal forces fora linearly varying ‘edge pressure or body force oad, Problems involving more complicated distributed loads ‘would require an clement of higher order. 6 ‘Figure 10.4: Plane stress linear triangle with point foree loadingm ™ normal tration shear traction Figure 10.5: Plane stress linear triangle with pressure loading Figure 10.6: Plane stress linear triangle with body force loading the elementis loaded only by point forces atthe nodes, as shown in Figure 104 then the lobal force vector can be assembled directly (inthe same way a8 forthe truss clement ‘scussedin Chapter), For example, consider the mesh of two linear plane stress elements shown in Figure 107. Figure 10.7: Mesh of plane stress linear triangles with point forces applied at nodes ‘Ignoring boundary conditions, the global force vector for this case is simply413 FT a [Pa 0 Pep yg 0 Py 0 Py) For the case of pressure loading onan element ede, the traetions are usually specified in {ermsof their normal and shear components as shown in Figure 10.5. Consider the loaded ‘edge shown in Figure 108, where the model coordinate F ranges from E = 1 atnode one toF = + Latnode to. Figure 108: Tractions on edge of plane stress linear triangle ‘Along the edge, the normal ang shear tractions may be expressed as n= Nid + Nada = Nyda + Meda where M=-}E-0) (00.33) R= 3€+) 0.34) sche unalone ier sap eons (0) athe psc a Smet oe mtngennsbuton tom he oplaiadg tonsotctnad Cae eae tom cnc tgcnadeatosheoniee i m0 0% ea Neat = | ofan thn, 0 iid (2038) of In this equation, the shape Functions are given by (10.11), (10.12) and (10.13) and are evaluated along an edge. As the three edges ofthe triangle are described by the equations £ = 0,9 = Gor + — 1 = O,thenoneofthese conditions will be true along the loaded16 ‘edge. This implies that one of the three shape functions will also be zero. For the case shown Figure 10.8, we have £ +1 ~ 1 = Oand hence Ny = 0, Since the pressures and slsplacements vary linearly in terms of the model edge coordinate & , equation (10.35) can be replaced by the equivalent expression (20.36) This is more convenient to evaluate numerically as the limits range from ~1 to +1 when the integral isexpressed in terms of. With reference to Figure 10.8, we see that resolving inthe horizontal and vertical directions gives “ de= joa + guna = aft + on$h oan 4 ay = aysina ~ gycosa = 9/5 ~ an SE (20.38) Now along the element edge we also have Mix, + Nox y= M+ May and hence a N e (2040) (aot)us where WW, oN, 00g N=/o N, 0 N,00 2 4 9,2) (act + gn] at een (40.42) Equation (10.41) may be integrated numerically using the Gauss rule to give the surface [pressure contribution tothe element force vector as AYN yw (2043) ‘here ms the numberof integration points, wis the value ofthe weight for integration pointi, 8, = NCE; )is the edge shape function matrix evaluated at integration point, ia evaluated at integration point, and £} is the value of the model edge coordinate & for integration point, Once each element force vector is formed, i then inserted into the slobal force vector using the element steering e=tor ‘The final type of loading to be considered, shown in Figure 10.6, i caused by body forces. ‘Neglecting contributions from other types of loading, the consistent nodal forces are given by the second term on the right-hand side of (10:31) according to ‘This equation may be integrated numerically aoording to Ta AS NTH, ders, aos) ‘where isthe number of integration points, isthe value ofthe weight for integration point i, N, = N(E,,1,) is the shape function matrix evaluated at integration point i, et J, = det JG,,7, )isthe determinantof the Jacobian evaluated at integration pointi, and (i,m) ate the vals of the mode coordinates for integration point 1046 Worked Example of Numerically Integrated Plane Stress L ar Triangle ‘Th illustrate the steps involved in forming the stiffness relations by numerical integration, a worked example is presented for the plane stress linear triangle. In performing the numerical integration by hand, we first note from (9.40) that the determinant of the Jacobian is constant and given by116 dety = 24 (1045) Substituting equations (9:36)(2.39) in (929), we also see that » » | % no oH wala #|"afy-n 4-9 a” & Using (10.27), the partial derivatives ofthe shape functions can then be found using aN, aN, aN |e =e x me &!_ [%- ys Ym My] avy aN, aN3)2a]x,— x) 4-H 2 —¥Y (20.46) Tasertng these derivatives in (10.21 we obtain the strsin-isplacement matrix Ds 0 yn 0 ye 0 Bed}O mm 0 xy 0 mm os) be Yn ta Jn Xn Ja where mato I aos set AT ae Math Meee In the plane stress constitutive relations defined by equation (10.8), let us define the D matrix Pu Pn o p=|P2 Pa 0 a048) 0 0 Dy with Ds nn = A + 2G, Diy = Aand Dy = 6. Substituting (20.45), (10.47 and (10.48) into (10.32) yields the element stiffness according to[Yas 0 a0] 0 xn Yale, e 0 x2 ene aps 2 A] wend aa] Ona yal ¢ Yi. 0 xy Ox Jy To evaluat this stifiness, 7 Dy 0 2Dn 0 0 Dsl fs 9 yn O yy 0 0 xm 0 x3 0 ra)e2dxw, Fan Yon Fas Yar Xt Yi is necessary to determine the highest order polynomial that ‘occurs inside the summation sign. Inthe case of the linear triangle the B matrix is made up of constant terms and thus corresponds toa polynomial of degree zero. Similarly, the Jacobian matrixis also constant. Thus we see thatthe highest degree entry inthe clement Stifinessmatrixisequal othe product ofconstant terms, Thus thehighestorderterminside the summation for kis therefore also a constant and, from Table 9.1, itis necessary to use only a one-point Gauss rule For m1, the stiffness matrix summation becomes ¥23 0 25] 0 aa Yn ke Ps © Aallyt De 9 [a] 2A) Ox yall g 24 Yig 0 ty 0 xa Yn Dyyavn ? Patna Puvatn Duratn Dy Yat Dsvatas Pyain Pvtatn Pun + * + Dsstais DesYors Dstt x- 4 U4) Pwo Pntmts Params + + DasYn%2 Daa DasYnts Pura Pinte Pwarn # 4 ? Pskatn Dy Yotn Py tata Duyatn Dutota Davarry + + # Pyknta Psion Padirs PP Diz 0] fr 0 yn 0 2 0 O xq 0 xy 0 xyheaa xd 0 Ds) 32 Yas *13 Yo1 Ya a2, 2 Panry + Dayiaya, Padsta + DasYtn Pret + Pada gmmevic: Puri # Dasa, Daviakn Das Ya%n Duran Divaus Toillustrate the formation of the element force vector, consider the edge loading shown ln Figure 10-5 for the special case where the normal edge traction is uniform and the shear traction is zero. Substituting q, = 0, equation (10.42) simplifies to i 2 where, from equations (10.39) and (10.40), the derivatives are given by X, ®, -le- Fn Bae} am, VX, A =1 BrP Pr oy) ‘To integrate equation (10.41) numerically using a Gauss rule, we need to determine the highest degree of the terms inside the suramation, From (10.33) and (10.34), the shape function matric i linear inthe edge coordinate . Since the terms in tare all constants, this implies thatthe highest order term in the summation is linear and, hence, a one-point rule is sufficient to integrate the force vector exactly. For the one-point rue in Table 4.1, equation (10.41) becomes, mE) 0 0 M@|r ong, realMGy 0 ff 20°91 0 WE] Fe - x) o 0 0 Substituing F, = Oand w, = 2we obtain 1/2 O-¥)] 0 1216 ay Hea — x) can p2 ff 22-1 ak) or-yp 0 1/2/]-$6@, ~ 2) He ~20) oo 0| oo Q Sitice x» ~ xj = Leos@ and yz ~ y, = Lsina, where L isthe length of the edge and ais its inclination to the x-axis, this force vector may also be written as ait. | sina 2 |ncosa ° 049 ‘As an immediate check on these values we note that S nodal xforces = gy X hk x sina = total road applied to edge S nodal y-frces = gq X AL x cosa = total load applied to edge i 10: Plane Stress Quadratic Tangle ' ‘As shown in Figure 10.9, the plane stress quadratic triangle has six nodes and two “degrees of freedom at each node. Using sn isoparametric formulation, quadratic functions fare used to model both the displacements and the geometry and the element may have ‘curved sides @ Qeasny ® @ OO: © ar) Go) ‘model coordinates Figure 109: The quadratic tiangle ‘The interpolation and mapping for the quadratic triangle elements has already been { discussed in Chapter 9. The displacements are computed using = Nit, + Nyy + Nyy + Nyt + Nous + Notts (049) 8 Nyy + Nyy, + Ngby + Navy + Navy + Neve 4030) ‘where the shape functions aze given by equations (8.21)(9.26) according to I N= 2OE- 17> M= 4 Ns= 2) ~ 1/2) y= AONE +=, Ny= 2E+n~1/2VE +91) Ne=-4GVE +9 It is convenient for later work to express equations (10.49) and (10.50) in the single equation120 yf =] 0M, ON, ON, OM ONS ON, {] [eeteeeaeaenral where the shape function matrix N and the element nodal displacement vector w are defined as Ny ON, ON, 0 Ny ON ON, 0 N=1 0 N, 0 .N, 0s 0.Ny 0 Ng O Ng Gost) we fla Hs ta vey Hs a Hs Hy Hy ve) (aos) ‘The strain-displacement matrix B is obtained by differentiating equations (1049) and (20.50) according to al INy aN, aNg i} ee an, fay ew OTN, 5 ° aNg [b a Bll. These are ofthe standard form where the strain vector eis defined in (10.7), the element displacement vector wis defined in (1052), and the strain-displacement matrix B is,2 (ios) For an isoparametric formulation, thex-andy.coordinates at any point within the element can be writen inthe form 2 Nyx + Naty + Nyy + Nate + Noxs + Nex (1054) y= Naya + Nova + Nays + Nave + Neds + Neds (20355) Differentiating (10.54) and (10.55) with respect to the model coordinates & and gives (1056) ‘The strain-displacement matrix B (10.53) contains partial derivatives of the shape functions with respect tothe real coordinates andy. These are computed inthe same way as for the linear clement by using equation (10.27). 10.8 Numerical Integration of Stiffuess Matrix for Plane Stress Quadratic Triangle ‘The virtual work equation for the plane stress quadratic triangle is identical to that for the Iinear plane stress triangle, except that the point force and virtual displacement vectors become Pra (Py By Pa By Pa Pa Pa Pa Ba Pa Pu Pos Ps Ps By} Ba = my by tp Bry bis dry buy dg dug Svs du, br Apart from this minor change the derivation ofthe stiffness equations forthe two types of elements are the same. To form the stiffness equations for the quadratic element, we substitute (108), (1053) and (1056) into. equation (10.30) and apply numericalm integration. The numberof integration points required depends on whether the triangle hhas straight orcurved sides. Inthe former case, che srain-Sisplacement matrix Bisa linear function ofthe model coordinates and and the determinant ofthe Jacobian is constant. ‘Thus the highest degree term inside the summation sgn is quadratic and, from Table 9.1, a three-point integration rule is necessary to evaluate the element stiffness exactly, When the element has curved sides, de is no longer constant and the terms in the B matrix ‘become very complicated. Indeed, from (10.56), we see that the determinant of the Jacobian may have terms which are of quadratic degree. Using (10.27), we also see that the derivatives in the B matrix re given by aN, = wy aN y © dety | ae am ~ on a which are ratios of wo quadratic polynomials. In equation (10.30), one ofthe 1/detJ terms cancels with the det in the summation but there is still one left over. Once the matrix ‘product is multiplied out, we conclude that the highest order term in the stiffness matrix is the ratio ofa fourth order polynomial and a quadratic polynomial. For an element which has highly curved sides, asuitable strategy would be to use a rule which integrates a fourth ‘order polynomial exactly. More commonly, where only one edge ofthe element is curved, ‘a three-point rule is usually of sufficient accuracy. 10.9 Formulation of Force Vector for Plane Stress Quadratic Triangle ‘Under point force loading, the global frce vector canbe assembled directly as described for the plane stress linear triangle. For edge pressure and body force loading itis generally necessary to use numerical integration to compute the consistent nodal forces. Since this clement uses a quadratic expansion to model the displacements, it can accommodate 2 ‘quadratically varying edge pressure or body fore load, ae Figure 10.10: Tractions on edge of plane stress quadratic triangle Consider the loaded edge shown in Figure 10.10. Along this edge the model coordinate has values of F = -1 at node one, F = 0 at node two, and & = + 1 atnode three, The ‘normal and shear tractions may be expressed asFez} 0 = Mian + Nodua + NsQus a= Nida + Noga + Mage, where x f 28E-) DE+1) me Le+y are quadratic one-dimensional shape functions and (q,;,q,) are the prescribed nodal values of (gq, 42). Neglecting contributions from other ypes ofloading, the consistent nodal forces forthe element are again given by equation (10.41) in which w_[ 0% 0M 0 000000 N=lo m0 N, 0 N000000 and Finally, we note that the consistent nodal forces caused by body force loading are again defined by equation (10.44), except that det J is now given by (10.56) and N is given by (033).CHAPTER 11 PLANE STRAIN AND AXISYMMETRIC ELEMENTS125 1 PLANE STRAIN AND AXISYMMETRIC ELEMENTS, ‘This chapter illustrates the development of two-dimensional elements for analysis of linear elastic solids under conditions of plane strain and axial symmetry. The governing equations are derived from the prineiple of virtual work and the consistent element force vectors ae established fora variety of loading conditions {L.A Stress-Strain Relations for Plane Strain Plane strain conditions occur wherever the loading inthe out-of-plane direction is ong. ‘Typical examples of plane strain loading are embankments and strip footings subject 10 uniform pressure. Away from the ends of the body we assume that the strain in the _zairection is 2r0, and therefore consider a two-dimensional slice of unit thickness. The three independent nonzero stess components are again Ga, dy, and Ty as shown in Figure 11.1. a0 Figure 11.1: Plane strain loading Sulsttuting the conditions e, = t3¢ = tay = Othe stress and strain vectors become the same as (10.6) nd (10.7). The corresponding stress-strain matrix is n47G 4a} of whe 4 au where G and. ate given, respectively by equations (10.2) and (10.3) The condition é, = 0 implies that the out-of-plane stress can be expressed in terms of the remaining stresses according to| 2, = (a, +93) 11.2 Plane Strain Linear Triangle ‘The derivation of the stiffness equations fr the plane strain linear triangle is identical {0 that forthe plane stress linear triangle except that equation (11.1) replaces equation (20.8 the stress-strain matrix. Also the thickness of the element, fis usually assumed {oe unity inthe stiffness matrix and the force vectors.aw 113 Stress-Strain Relations for Axisymmetry _Axisymmetric deformation occurs when the body and the loading is symmetric in the cylindrical coordinate estem rz and 9, It occurs, for example, in the analysis of circular footings and piles subject to a vertical load. The four nonzero stress components are dy, 5 opand ras shown in Figure 11.2. Figure 112: Axisymmetric loading Under conditions of axisymmetry, the three-dimensional stress and strain vectors reduce 10 OF & (0, 05 te Oph «ay T= (ere Ye a3) ‘The corresponding stress-strain matrix becomes a+2G 2 0 2 a at7G 0 4 0 0) GOO aay a a 0 2+26 Note that the entries inthe first three rows and columns of the axisymmetric stress-strain matrix are identical to those inthe plan strain case. Thus the axisymmetric constitutive relations are obtained by simply augmenting the plane strain constitutive relations with an additional row and column. 11.4 Asisymmetric Linear Triangle ‘Axisymmetric elements, although usually drawn in two-dimensions, are actually ring shaped as shown in Figure 11.3 The interpolation andmapping for the axisymmetric linear triangle is identical to that forthe plane sires inear triangle. The displacements are again Nyy + Noita + Nguy aus) Nye + Nava + Ng5 a6)28 ‘where the shape functions are given by equations (9.14), (9.15) and (9.6). = Pin x9® gue 113: The aaymaett liner tanele 1tis gain convenient to expres these equations inthe matrix form fi iu) _ [v0 0m OTe -esataflee on eed 3) ith the shape function matrixN and the element nodal displacement vector u being given by (10.15) and (10.16). ‘Under axisymmetric loading, the strains are defined according to as) ais) i | 1.10) ‘This leads to the plane stress strain-displacement relations e=Bu au) ‘with the strain vector e defined in (11.3), the element displacement vector u defined in (20.16), andany aN, bar Oe o Bo : (ana) Jax, aN, aM eee x N, foo F For an isoparamettc axisymmetric triangle, the radial and vertical coordinates for any poiat inside the element may be expressed inthe usual form Nyry + Nata + Nats «uy Ney + Nata + Naty au and the determinant of the Jacobian becomes ae dete ab on ~ O70 Differentiating (11.13) and (11.14) with respect to the model coordinates & and 7 and substituting inthe above gives os (Se [SHe}-E with respect tthe coordinates and. Asbefore, these derivatives cannotbe founddiectly 1nd must be computed using the chain rule of differentiation according to et = a (11s) a] (aN) fant alla bar| az})anf ~ "any ate} (ae (aN; (aN a = my, (a1.16) where iy ale,130 ‘Equation (1.16) enables the derivative terms in the B matricto be computed for any value of the model coordinates & and. The Nir terms, which arise because ofthe hoop strain, ‘an easily be found by using the shape functions definitions together with (11.13) 11.5 Formulation of Axisymmetric Linear Triangle Using Virtual Work ‘The virtual work equation for an axisymmetric element may be written as (6010; + d8202 + Sate + degagha¥ = | (bugr + dvajas + f f f 2 Phy tn + Pb + PBs + where ¥ denotes the element volume, (Ber 8; 37m de) ate the vistual strains, (0, 0 tn ay ate the streses, (By) are the virtual displacements, (gq) are the surface presszes inthe rand2-irections applied over the area, (Jf Z) ae the body forces inthe F ands-drections, (b4,d,) ae the virtual displacements at node i, and (P,P) are the Poin (fing) load in ther and»-direcions applied at node i. r Y= rdrded = rdAds Figure 11.4: Axigmmetric volume, dS = rdLdo ‘Figure 11.5: Axisymmetric surface area [Noting the geometry shown in Figure 11.4 and Figure 11.5, the equation may be written in matrix form as symmetric virtual work31 a {tara = 2 | satura san aatordaeeton cuay | where deh = (ber Bee Oye Bey} - oT = [oy o: te 04) da? = ou oy C= lo a Ry [ (Pa Pa Pa Pa Pa Pa} aul = {04 87, dy Ory buy by) ‘Now from (11.7) and (1.1) itfollows that ea = utNT f det = buat ‘These equations together with the axisymmetric form ofthe stress-strain matrix (11.4), can be substituted in (11.17) to give stax | wtbpact 1rda e an | Waray + d0t20 thr + oa'h Since the variation du is arbitrary, and the virtual work principle must hold forall such variations it follows that ax wore arian ax wterat anf Nioraa <9 “Thus the element stiffness and force vector are k= 25 { BtDBéets rat 18) af strats af Nodes cuss) 11.6 Numerical Integration of Stif mess Matrix for Axisymmetric Linear Triangle ‘To derive the stiffness equations for numerically integrated element, we substitute (10-25), (11.4) and (11.12) into equation (11.18) and apply numerical integration. This ives132 K= 20) BPDB det; w; 120) ‘where mis the numberof integration points, isthe weight for iategration point, isthe r-coordinate for integration point i, By = B(S,m) i the staindisplacement matrix evaluate at integration point, detJ, = detJ.7,) isthe determinant of the Jacobian evaluated at integration point i and (,.7,) ate the values ofthe model coordinates for Integration point “Unlike the plane ses and plane strain linear triangle, the stiffness smatraxfor the axisymmetric linear triangle cannot be computed exactly de to the Ir term {nthe Bmatrx. Allhoughone of thelr terms ances with ther term inside the summation, ‘hore isstill another left which means that some of the integrand entries are in fact ratios of polynomials. Tn order to obtain an aocurate, Dut not exact, sffness atx for axisymmetric conditions it isusuallysuficent to choose a ule which ignores te effect of ‘the emainingIrterm. Expanding the matrixproduct inequaton(1i.20),thebighest order term conteins products ofthe shape functions divided bythe reoordinate. Ignoring the effect ofthe later, we would choose @ 3-point Gauss integration rule which integrates a quadratie polynomial exact. 11.7 Formulation of Force Vector for Axisymmetric Linear Triangle ‘Under point force loading, the global force vectorcanbe assembled direcly as described for the plane stress case. It is important to note, however, that a point force under suisymmetric conditions actually corresponds to a ring line load. For edge pressure and body force loading itis again convenient to use numerical integration to compute the consistent nodal forees, ar Figure 11.6: Tractions on edge of axisymmetric linear triangle ‘Consider the loaded edge shown in Figure 11.6, which sweeps out an annular surface about ‘the vertical axis. The normal and shear tractions are again written as 90 = Bian + X44, a= Midu + Nida where the one-dimensional shape functions N, are given by (10.33) and (10.34) and @nis gu) 86 the prescribed nodal values of (gn, 4,). Following the same argument that was ‘developed for the plane stress element in Section 10, the consistent nodal forces, which are given by the first term on the right-hand side of (11.19), may be expressed as1 0 0M I, 0 |(27 wan] IR ol iran 1.21 | of, 4s] (ait) ray 00 “The tractons shown n Figure 11.6,whenresolvedin the radial and vertical directions give - sina = q,$e + gy ay = qycosa + qgsina = 4,5F + ane 122) ge gsina ~ qyoosa = 4/$2 ~ ane (1123) [Now along the element edge we aso have = Mn+ Mr (4124) My, + yz and hence 125) 11265) Combining equations (11.22)-(11.26), and inserting in (11.21), gives the consistent nodal forces as (quan where,(aia) Equation (11.27) may be integrated numerically using the Gauss rule to give the surface pressure contribution to the element force veetor as fam YN try 129) ‘where mis the number of integration points, w; is the value of the weight for integration point, 1 isther-coordinste for integration point N, = N(G,)is the edge shape function matrix evaluated at integration point /, value of the model edge coordinate & for integration point i iis evaluated at integration point i, and E, isthe For body force loading, the consistent nodal fores are given by the second term on the right-hand side of (1119) according to IN; 0 rt] 0 Mr Iw. 0 |r} NTbrd4 = 24 0 Nalfaf dott rand ojo Iw, 0 on, ‘This equation may be integrated numerically according to t= an SNP, det, cum tore isthe mmberfintegration pots, whe vale of the weight for integration point, isthe rcoordiate fo integration plat, N, = NG) isthe shape nton fart ctaloated at intgration point dx, ~ dG) te detertinact ofthe Tooobian evaluated at integration poat fend (Em) ae the values of the model coordinates for integration ponte 11.8 Worked Example of Force Vector for Axisymmetric Linear Triangle Consider the edge loading shown in Figure 11.6 forthe special case where the normal edge traction is uniform and the shear traction is zero. Substituting g, = 0, equation (11.28) becomes 22) anil ‘where, from equations (11.25) and (11.26), the derivatives are given by15 ale = a) {In equation (11.29), the shape function matrix is linear in the edge coordinate &. From (11.24), we see that this also true forthe radius?. Since the terms int are all constants for ‘uniform applied traction and a straight edge, this implies thatthe highest order term in the summation is quadratic. For the tvo-point rule in Table 4.1, which is exact for a cubic polynomial, equation (11.29) becomes INE) 0 0 FGD|r ae lee tan HE) _o | tin 0 ME)||-Fe- o 0 oO MG, 0 o NE) a) + 2n[ ME) 0 Yb (Gan + RE) «2 0 N€)| 7 o 0 ‘0 Substituting mat sy met 132) we obtain te $e, - 20) L la-4y,+4 x(Ja-4y, +40 a ( Bt 20 ty136 +2 0-2) . | (o+dr ete which reduces to @-a)Gn +4n) Ja — rr + $12) a, +2, gery @2-20GR + 3p) He. — Gri + 3r2)| o ° Since ry ~ ry = Leosa and vector may also be written a 2; = Lsina, where Li the length of the edge, this force Gr, +4r)sing |-Gr, + 41a eosa f= guald Gr + 3pa)sing Gr: + Sr oosar 0 0 ‘Asan immediate check on these values we note that Snot rfores = gy x 20d, + Lx sna = ttl road applied to edge Sonoda forces = = 4 x 2dr, + rE x cosa total elond applied to edgeseo