A NEW EIGHTEEN PARAMETER TRIANGULAR ELEMENT FOR GENERAL PLATE AND SHELL ANALYSIS by TERRANCE WILLIAM BEARDEN B.A.Sc. (1974 ) The Un ivers i ty of Manitoba A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in the Department of CIVIL ENGINEERING We accept th is thes i s as conforming to the required standard The Un iver s i t y of B r i t i s h Columbia A p r i l , 1976 (c^ T e r r a n c e W i l l i a m B e a r d e n , 1976 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t he r e q u i r e m e n t s f o r an advanced deg ree a t t he U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s tudy . . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y pu rpo se s may be g r a n t e d by t he Head o f my Department o r by h i s r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i t h o u t m y w r i t t e n p e r m i s s i o n . T.W. Bearden Department o f C i v i l E n g i n e e r i n g The U n i v e r s i t y o f B r i t i s h Co l umb i a 2075 Wesbrook P l a c e V a n c o u v e r , Canada V 6 T 1W5 A p r i l , 1976 - - i i -ABSTRACT The purpose of th i s invest igat ion was to develop an eighteen parameter f l a t t r iangular f i n i t e element for analyzing plate and shel l structures. The development of the element was accomplished by combining a plate bending element with a new . plane stress element. The well known nine parameter t r i ang le using the normal displacement and two slopes at each vertex was used for the plate bending element. This element contains an incomplete cubic for the normal displacement. For the in-plane element, complete cubics were used i n i t i a l l y for the displacements^and then various constraints were imposed to reduce the number of generalized co-ordinates to nine, namely the two in-plane displacements and an in-plane rotat ion at each vertex. One of the constra ints , namely that the included angle at each vertex was invar iant , destroyed the completeness of the element. However, the element was compatible in the plane. A patch-type test of the in-plane element showed that i t could not represent a l l constant s t ra in states exactly. However, the errors were small. The complete element was then tested on a plane stress cant i lever beam, a square plate subjected to membrane stresses only, a c y l i n d r i c a l s h e l l , a spherical shel l and a non-prismatic folded plate structure. In a l l cases, reasonable engineering accuracy was achieved with modest grids of elements. Thus i t was concluded that the incompleteness of the in-plane element was not too important. F i n a l l y , a compatible beam element was formulated and tested to supplement the t r iangular element. The beam element formulation included unsymmetric crosssections. - i i i -TABLE OF CONTENTS Page Abstract i i List of Tables v List of Figures v i i Symbols z Acknowledgements x i i i Chapter 1: I n t r o d u c t i o n 1 Chapter 2: General Information 3 2.1 F i n i t e Element Technique 3 2.1.1 D e s c r i p t i o n of Method 3 Chapter 3: D e r i v a t i o n of the Element's P r o p e r t i e s 5> 3.1 General Information 5 3.2 In-Plane Formulation 6 3.2.1 I n t e g r a t i n g the S t i f f n e s s M a t r i x 23 3.2.2 S t a t i c Condensation of C e n t r o i d a l Degrees of Freedom 2U 3.2.3 C h a r a c t e r i s t i c s of the Plane Stress Element 25 3.3 Bending Element Formulation 25 3.4 Assembling the In-Plane and Bending Element S t i f f n e s s e s 37 3.5 Co-ordinate Transformations hO 3.6 Element Dimensions U6 3.7 Summary of the Combined Element Chapter 4: Stress Computations 4.1 In General 1+8 4.2 In-Plane Stresses 5*0 - i v -4.2.1 Consistent Formulation 5o 4.2.2 Constant S t r a i n Formulation 52 4.2.3 L i n e a r S t r a i n Formulation 57 4.3 Bending Stresses 70 Chapter 5: Beam S t i f f e n e r Element 71 5.1 Symmetric Bending 72 5.2 Unsymmetric Bending 83 Chapter 6: Numerical A p p l i c a t i o n s 89 6.1 Constant Stress A p p l i c a t i o n s 89 6.2 C a n t i l e v e r Beam Problem 9U 6.3 P a r a b o l i c a l l y Loaded Square P l a t e 99 6.4 C y l i n d r i c a l S h e l l Roof 112 6.5 Poi n t Loaded S p h e r i c a l S h e l l 121 6.6 Non-Prismatic Folded P l a t e S t r u c t u r e 1 3 0 6.7 Beam S t i f f e n e r A p p l i c a t i o n 1u2 Conclusions 1U5 Appendix A: Computer Program Information 1u8 A . l D i s c u s s i o n of Program 1u8 A.2 Input Data 151 A.3 Flow Chart 160 B i b l i o g r a p h y - V -LIST OF TABLES Table 5.1 S t i f f n e s s M a t r i x f o r Symmetric Beam S t i f f e n e r Element C a n t i l e v e r Beam Problem: 6.4 Tip D e f l e c t i o n and Normal Stress 6.6 Stresses (C.S.T., L.S.T., Consistent Formulation and Exact) Page 3.1 Trigonometric R e l a t i o n s f o r the Element 10 3.2 Transformation M a t r i x f o r Plane Stress Element 18 3.3 Area Co-ordinates f o r Numerical I n t e g r a t i o n of the Bending Element 36 3.4 Strain-Displacement M a t r i x f o r P l a t e Bending 3k 4.1 Strain-Displacement M a t r i x f o r L.S.T. Stresses 61+ 82 Constant Stress A p p l i c a t i o n s : 6.1 D e f l e c t i o n s f o r Constant Shear Stress 92 6.2 D e f l e c t i o n s f o r Constant Normal Stress 92 6.3 D e f l e c t i o n s f o r Constant Moment 93 97 P a r a b o l i c a l l y Loaded Square P l a t e : 6.5 D e f l e c t i o n s and S t r a i n Energy f o r Various G r i d Sizes 103 105 C y l i n d r i c a l S h e l l Roof: 6.7 D e f l e c t i o n s f o r Various Gridworks 115* 6.8 Stresses (L.S.T.) Poi n t Loaded S p h e r i c a l S h e l l : 6.9 D e f l e c t i o n s f o r Various G r i d Sizes 121+ 6.10 Stresses (L.S.T.) 121+ 6.11 Stress Comparison (C.S.T., L.S.T., Consistent Formulation, and Exact) ^ ^ - v i Non-Prismatic Folded P l a t e : 6.12 D e f l e c t i o n s vs G r i d Sizes 6.13 L o n g i t u d i n a l Stresses Appendix A A.2.1 Format of Input Data Cards - v i i -LIST OF FIGURES Figure Page 3.1 Co-ordinate Systems 8 3.2 Rotations for an Element 12 3.3 Degrees of Freedom of Bending Element 2 7 3.4 Area Co-ordinates 28 3.5 Co-ordinate Systems Ui 4.1 Constant Strain Triangle 5 2 4.2 Area Co-ordinate of L.S.T. 57 4.3 Tangential Displacements (for L.S.T. Stresses) 65 4.4 Normal Displacements (for L.S.T. Stresses) 67 5.1 Beam Stiffener Element 71 5.2 Beam Stiffener Geometry (Strong Direction) 72 5.3 Beam Stiffener (Degree of Freedom Strong Direction) 78 5.4 Beam Stiffener Geometry (Weak Direction) 7 9 5.5 Beam Stiffener (Degree of Freedom Weak Direction) 80 5.6 Beam Stiffener (Torsion) 80 5.7 Resultant Beam Stiffener (12 Degrees of Freedom) 81 5.8 Beam Subjected to Couples 83 5.9 Deflected Beam Under Pure Bending 86 Constant Stress Applications: 6.1 Constant Shear Stress 91 6.2 Constant Normal . Stress 91 6.3 Constant Moment 91 Cantilever Beam Problem 6.4 Cantilever Beam (Grids and Loading) 9 6 6.5 Tip Deflection vs_ No. of Degrees of Freedom 9 8 P a r a b o l i c a l l y Loaded Square P l a t e : 6.6 General Layout and Loading 102 6.7 S t r a i n Energy vs T o t a l Number of Degrees of Freedom 10lt 6.8 10 vs T o t a l Number of Degrees of Freedom 10U 6.9 N - and N „ vs F i n i t e Element G r i d S i z e xD yB 108 6.10 6.11 6.12 N . vs F i n i t e Element G r i d S i z e yA — N vs F i n i t e Element G r i d S i z e yc — N . vs F i n i t e Element G r i d S i z e yD — 109 110 111 C y l i n d r i c a l S h e l l Roof: 6.13 General Layout 11U 6.14 W - D e f l e c t i o n Along Edge B - C 116 6.15 Wg vs T o t a l Number of Degrees of Freedom 117 6.16 N x Along Edge A - B 119 6.17 M y Along Edge D - C 120 Point Loaded Sphere: 6.18 General Layout and Loading 123 6.19 D e f l e c t i o n at Pole vs F i n i t e Element Gr i d 125 6.20 Normal Displacement vs Angle $ Near Pole 126 6.21 Membrane Stresses vs $ Angle Near Pole 128 6.22 Membrane Stresses vs $ Angle Remote from Pole 129 Non-Prismatic Folded P l a t e S t r u c t u r e : 6.23 General Layout and Loading 133 6.24 P l a t e Geometry 134 6.25 Model and F i n i t e Element Mesh Patterns 135 6.26 V e r t i c a l D e f l e c t i o n Along Fold Line c 137 6.27 L o n g i t u d i n a l Stress Along Fold L i n e c 139 - i x -6.28 L o n g i t u d i n a l S t r e s s A l o n g F o l d L i n e E (c.L.^ 6.29 T r a n s v e r s e Moment a t Midspan Beam S t i f f e n e r Problem: 6.30 G e n e r a l Layout 6.31 Load Case 2 A.1.1 Beam S t i f f e n e r S e c t i o n P r o p e r t i e s LIST OF SYMBOLS D e f i n i t i o n Area of t r i a n g l e Column vec t o r of polynomial c o e f f i c i e n t s Element dimension, Figure 3.5 C o e f f i c i e n t s of the displacement polynomials, Equation 3.2 Strain-displacement matrix Cosine & s i n e of the angle a, Figure 3.1 Centre l i n e Constant s t r a i n t r i a n g l e f o r m u l a t i o n E l a s t i c i t y matrix Modulus of e l a s t i c i t y of the m a t e r i a l being modelled & that of the f i n i t e elements E c c e n t r i c i t y , Figure 5.2 Equation Modified E u l e r ' s beta f u n c t i o n , Eqn. 3-33 Load ve c t o r Figure Shear modulus of e l a s t i c i t y Inches Moment of inertiaf'is P o l a r moment of i n e r t i a S t i f f n e s s m a t r i x Area co-ordinates used i n the l i n e a r s t r a i n t r i a n g l e Length of beam s t i f f e n e r element Area co-ordinates used i n p l a t e bending f o r m u l a t i o n - x i -L.S.T. Li n e a r s t r a i n t r i a n g l e N Number of s u b - d i v i s i o n s ( g r i d refinement) N_^ Shape functions used i n the p l a t e bending f o r m u l a t i o n N" N N Membrane s t r e s s e s x" y ' xy N e N* % {P} Load vector T^,T^,T^ Length of element's sides r ,r R a d i i of curvature y' z [T] Transformation m a t r i x ex. Equation 3-28 U S t r a i n Energy u, v Displacements i n the x & y d i r e c t i o n , r e s p e c t i v e l y un ^ j Normal displacement at node i to node j i i Normal displacement at node k i i j _ _ Tangential displacement at node i to node j ii k Tan g e n t i a l displacement at node k w Normal out of x,y plane displacement W. ,W I n t e r n a l & e x t e r n a l work i ? 'e x,y,z G l o b a l c a r t e s i a n co-ordinates a Angle tangent to element s i d e & a x i s , Figure 3.1 {6} D e f l e c t i o n vector {e:} S t r a i n v e c t o r L o c a l co-ordinates, Figure 3.5 {a} Stress v e c t o r 6,cf> Angle, Figure 6.12 [X] D i r e c t i o n cosine m a t r i x •Q^ Shear s t r a i n at node i - x i i -v Poisson's r a t i o to R o t a t i o n ^ Corresponds to Subscripts A With respect to polynomial c o e f f i c i e n t s b F l e x u r a l a c t i o n (bending) c Centroid of element G Glob a l co-ordinate system L L o c a l co-ordinate system P Membrane a c t i o n (plane s t r e s s ) x,y,£,£ Denote d e r i v a t i v e s of displacement w i t h respect to x,y,S,C 6 With respect to the a c t u a l degree of freedom Superscript T Denotes t r a n s p o s i t i o n of rows and columns of a m a t r i x S p e c i a l Symbols [ ] Denotes a matr i x { } Denotes a column v e c t o r - x i i i -The author would l i k e to express his appreciation to his advisor. Dr. M. D. Olson for his assistance, encouragement and for reading and checking the entire thesis. Also, the author would l i k e to thank Dr. N. D. Nathan for his assistance. The author would l i k e to extend his thanks to Sarah Dahabieh for typing the manuscript. I am also very grateful to the Defence Research Board and the National Research Council for their financial assistance. - 1 -CHfiPTER 1 IJS1TR0DUCTIDN The method of f i n i t e elements originated about twenty years ago i n the f i e l d of engineering and has since developed immensely. The basic idea behind the method i s that a solution region can be approximated by replacing i t with an assemblage of descrete elements. The f i n i t e element procedure reduces the problem to one of a f i n i t e number of unknowns by d i v i -ding the solution region into elements and by expressing the unknown f i e l d i n terms of assumed approximating functions within each element. The appro-ximating or interpolation functions used herein are defined i n terms of the values of the displacement f i e l d variables at specific points called nodes. These nodal variables are the unknowns which are solved for. The interpo-lation functions cannot be chosen a r b i t r a r i l y because certain compatibility conditions have to be satisfied. The accuracy of the solution depends not only on the size of elements used but also on the interpolation functions incorporated. One major advantage of the f i n i t e element method i s that the force - displacement or stiffness characteristics of each element can be computed and then the elements assembled to represent the stiffness of the overall structure. When choosing the interpolation functions that are to be incor-porated i n deriving an element's stiffness characteristics, discretion has to be used. The higher the order of the functions used the more complex the formulation becomes and the problem size increases greatly demanding more computer memory to be u t i l i z e d . However i f the polynomials are very low i n order then accuracy can be lost even though a great many elements can be -2-used in the. analysis. It is. desired to develop, a relatively low order trian-gular finite element .in .plane stress: which.when cxrf>lned..with. a triangular plate bending element, can be used .to model shell and folded plate structures with, reasonable accuracy:, and ecohamy. Starting, with., complete cubic polyno-mials to represent the two in-^plane displacements and. an in-plane rotation at each node of the plane stress element, various constraints are then introduced to reduce the number of degrees of freedom to nine for the element. This element was then combined with, the well known Zienkiewicz nine parameter plate bending triangular element which uses a cubic polynomial for the normal displacement. A computer program employing the new eighteen degree of freedom triangular finite element was developed. Various shell structures and a folded plate one were analyzed and results were compared with analytical solutions. Subsequently a twelve degree of freedom unsymmetrical beam stiffener element was formulated so that stiffened plate and stiffened shell structures could also be modeled. - 3 -CHAFIER 2 .. GENERAL IJ^RMATION 2.1 Finite Element Technique: The assumptions used in the finite element method herein are: 1] The element' s thickness is uniform. 2) The.material Is elastic, isotropic and homogeneous. 3) Elements are assumed to be connected onlyvnode points. 4) Relation between forces and deformations is linear. 5) Small deflection theory is assumed from plate theory; therefore there isn't any coupling of the in-plane and bending actions. 2.1.1 Description of Method - Displacement Approach: Using the potential energy (P.E.) principle, we assume a displacement field within the element. For equilibrium, the P.E. is a minimum and the internal work (strain energy) is equivalent to the work done by the external forces acting on the element. From this approach the stiffness characteristics of the element can be defined. This is illustrated below: The stresses in a continuum are expressed in terms of strains { a } = I D ] { e } 2 - 1 where { a } = Stress Vector I D J = Elasticity Matrix { e } - Strain Vector The strains at any point within an element can.be described in terms of the nodal displacements as { e } = I B ] { S' } 2 - 2 -k-where { 6 } = nodal displacement vector I B ] = strain-displacement matrix. * assuming a v i r t u a l displacement {6} at the nodes, the external work V/e done by the nodal loads { P} i s : W = { <S }*T { P } 2-3 e Similarly the internal work done by the element when subjected to the v i r t u a l displacement i s : W± = / { e }*T { a } dvol. 2-4 vol. substituting equations 2-1 and 2-2 into 2-4 yields W± =;'{« J* 1 I B ] T [ D ] [ B ] ' { 5 } dvol. 2-5 equating the internal work with the external work yields { 5 } { P } = { 6 } I B ] [ D ] [ B ] dvol { 6 } 2-6 then for an arbitrary v i r t u a l displacement {• 6 } r v o l and { P } = [K] { 6 } 2-8 { p > = A ^ i t B ] T t D ] [ B ] dvol { 6 } 2-7 So I K ] = [ B ] T [ D ] [ B ] dvol. 2-9 where [ K ] = element stiffness matrix -5"-CHAFTER 3 DERIVATION OF. THE ELEMENT'S PROPERTIES 3.1 General Information; A triangular element i s used because i t s shape affords easy applica-tion to many types of problems where rectangular elements could not be used. For example modelling odd shaped objects and desiring sub-sequent grid refinements i n regions of high stress gradients. This i s i llustrated later with a non-prismatic folded plate roof and various shell roofs. I t i s assumed that the behavior of a continuously curved surface can be adequately represented by the behavior of a surface b u i l t up of small, f l a t elements. From plate theory small deflections are. assumed so that the in-plane and bending actions are assumed uncoupled within each f l a t element. It i s desired to make the f i n i t e element as near to being compatible as possible. A compatible element i s one which satisfies sufficient inter-element continuity requirements that the total potential energy i n the structure converges monotonically towards a minimum as the (7) mesh of finxte elements i s progressively refined '. The potential energy i s a iidnimum when; among a l l the kinematically admissible dis-placements, those satisfying the equilibrium conditions make the potential energy stationary. The definition of compatibility may also be expressed as follows; i f a dependent variable i n a structure enters the energy expression with highest derivative of order q > 6, then the (q - 1) derivative of that variable must be continuous (7) between adjacent compatible elements . For plate bending, the - 6 -highest derivative i s two so the f i r s t derivative of the normal dis-placement or the slope must be continuous. The element to be compa-ti b l e must have continuous slopes (rotations) and displacements for modelling plate and shell structures. For i n plane or membrane action, the highest derivative i s one, so only the displacements and not the slopes have to be continuous for compatibility. Convergence to the correct ndnimum. potential energy i s obtained i f the polynomials are complete to order P. Where P i s the maximum derivative i n the energy expression. Only completeness to order P i s necessary for convergence. The f i n i t e element described herein i s the result of combining an in-plane and a plate bending element. For the in-plane portion the highest derivative i n the energy expres-sion i s one, so only complete f i r s t order polynomials i n u and V are required to ensure convergence of the potential energy. The energy expression for plate bending has highest derivatives of order two. Then at least a complete quadratic polynomial must be used for the normal displacement to ensure convergence. 3.2 In-Plane Element Formulation: As mentioned earlier i t i s desired to combine a 9 degree of freedom triangular plane stress element with the well known Zienkiewicz 9 parameter plate bending triangular element to represent folded plate and shell structures. So the two displacements u and V and an i n -plane rotation are used at each node to define the plane stress f i n i t e element, (refer to f i g . 3.5) - 7 -Beginning with complete cubic polynomials for each of the two i n -plane displacements, constraints are introduced to force the displace-ment p a r a l l e l to an edge to vary only l i n e a r l y along that edge and to force the included angle at each vertex to remain f ixed. Condensa-t ion of the remaining two degrees of freedom then yields the 9 parameter element. Then pro ceding as mentioned above, the 9 x 9 st iffness matrix i n l o c a l co-ord i s developed-. Starting with complete cubic polynomials: 2 2 3 2 2 u = a-^ + a 2£u+ a 3c + a ^ + a^z, + a gC + a^K + aQ£, t, + a ^ + 3 a 1 0 ? 3-1 y= a l l + a 12 5 + a 1 3 ? + a 1 4 ? 2 + a15^ + a 1 6 ^ + a 1 7 5 3 + a 1 8 ^ + a i g 5 c 2 + a2Qt3 3_2 But constraints are to be introduced to force the displacement p a r a l l e l to an edge to vary l inear ly along i t , so, u can be rewritten omitting the squared and cubic terms i n £ only, (-for s i de one.) . ^ r ^ J . r ^ 2 ^ 2 ^ 2 ^ 31 3 - la u = a 1 + a 2 £ + a 3 ? + a^ + a^. + a g C e + a^t ; + a g? In series notation: u = i a. ^ ? p i 3-lb i = 1 1 1 0 l i n i V = S a . -> ' 5 1 1 C 1 3-2a i = 1 1 + 8 -8-where { m } T = ^ 0 1 0 1 0 2 1 0 ^ ) ' i P = <0 0 1 1 2 1 2 3 ) { 1 } T = < 0 1 0 2 1 0 3 2 1 0 ) { n }T = <0 0 1 0 1 2 0 1 2 3) So i n i t i a l l y we begin with 18 parameters and wish to reduce these to 9. F i r s t Reduction: Force a displacement p a r a l l e l to an edge to vary linearly along i t . (refer to f i g . 3.1) i i l e t s = Sin a c = cos a then u = uc + Vs I I 3-3 V = - us + v*c Fig. 3.1 Co-ordinate Systems also £ = 1£0 + xc - ys X, = X S + yc 3-4 -9-Referring to equations lb and 2a and substituting equations 4 then; 8 * u = s= a. ( to + xc - y s ) n i ( xs + yc ) p i 3-lc 1 = 1 and V = £ a... - ( £ 0 + xc - ys ) 1 1 ( xs + yc ) n i 3-2b i = 1 1 + 8 but u = uc + Vs from equations 3 so the tangential displacement along an edge i s u u = c [A- a. ( So + xc - ys ) n i ( xs + yc ) p i ] + i = 1 10 _ _ l i _ _ n i + s [ ^ a^*( Co + xc - ys) ( xs + yc ) ] 3-5 i = 1 1 8 and we are interested i n u along an edge, where y = 0 therefore u = c [ J a. ( Co + x c ) m ( xs ) p i ] + s [ £ a. . 0 (. Co + xc f 1 i = 1 1 i = 1 1 + 8 ( x s ) n i ] 3-6 For u to vary linearly along an edge, we want the squared and cubic terms of x to vanish: Squared terms: sc a 1 2 + cs (ca 4 + sa 1 3) + s^ Cca5 + sa 1 4) + 2 + s3c Co a 1 5 + cs2£ 0 (ca 6 + sa 1 6) + , 2 + s Co Cca ? + sa 1 7) = 0 3_7 - 1 0 -Cubic Terms: 3 2 2 sc a, c + c s (ca, + sa n r) + s c (car> + sa, _) + 15 6 16 5 17 3 + s (ca g + sa 1 8) = 0 3-8 Note: Equations 7 and 8 are constraint equations for sides 2 and 3 of the element. Therefore actually 4 constraints are applied, leaving 5 parameters to be removed to ' ;y'ield the 9 desired. The cos a and s i n a ( c and s ) should actually be subscripted where j = 1, 2, 3 (side number). Table 3.1 Trigonometric relations for the Element Side j Cj Sj 1 1 0 0 Z -a/r 2 c/r 2 a 3 b/r 3 c / r 3 -b where r± = a + b, r 2 = 7 a 2 + c 2 , r 3 = 7 b 2 + c 2 the lengths of the 3 sides of the elements. are In-Plane Rotations: Define the rotation of one side of the element to be where '3 V = 3_J_ 3_V 3_C 9_V 3 X 3 X 3 ? . 3 X 3 ? from equations 4 and 3 i - 1 = c = s 3 x 3 x therefore 3 V _ „ 3 V/ , B 3 V — — = c + s — 3 x 3 5 3 £ c ( - u^ s c ) + s ( - u s +:.V c) 3-9 where u = - ^ i , etc. % 9 £ from equations lb and 2 a * 1 = 1 1 ^ 5 i = 1 1 + 8 8 i = 1 a. p. 1 *x mi pi-1 10 i = 1 ai+8. n i g 1 1 C 1 1 1 " 1 then r 8 n i pi-1 10 l i n i - l i n„ + s. [ - s. g? a.p.? + c. a. _ n . r r ] 3-10 J i = 1 J l = 1 for the j t h side of the element Define; The rotation at a node to be the average of the 2 side rotations at the node (refer to f i g . 3.2) Fig. 3.2 Rotations for an Element UJl = U J 1 2 + U J13 2 ^2 = OJ 2 1 •+'*23 3-11 2 w3 = "32 + % 1 -13-Rotation at node '.(1) • Co-ordinates (E,^, ^ ) = C -b, 0 1 U J 1 2 = a 1 Q - 2a 1 2b + 3a 1 5b 2 3-12 2 2 UJ 1 3 = - bca 2 + b C a 1 Q - 2a^o•+ 3a 1 5b ) + 2 ~2 r3 r3 2 2 2 - c C a 3 - a^. ;b + a gb ) + be ( a l ; L - a 1 3 b + a l gb ) 3-13 2 ~2 r3 r3 then -be bf. b 2, o r i + + " l = r 2 a2 + a10 ( r 2 + 1 ) " 2 b a i 2 ( 1 + r 2 > + 3 a 1 5 b r 2 - ^ C a 3 - a*:1 b + a gb 2) + ^ | ( - a 1 3b + a l gb 2)0* | r3 r3 3-14 Rotation at node (2): ( K2i C 2 ) = U , 0 } 2 ^21 = a10 + 2 a 1 2 a + 3 a 1 5 a 3 - 1 5 W i 3 = ^ a 2 + ^ C a10 + 2 a12.. a + 3 a 1 5 a 2 ) + r2 r2 c 2 2 ca 2 " " 2 ( a3 + a 4 a + a 6 a } " ~2 ( a l l + a 1 3 a + a 1 6 a } r2 r2 3-16 -14 -t h e n 2 2 2 u)2 = I ^ a e + a 1 Q C 1 + \ 1 + 2aa 1 2 ( 1 + ^ ) + 3a 2a l f j C 1 + ^ ) + r2 r2 r2 r2 2 - - j - ( a 3 + a 4a + a ga ) - <| ( a n + a^a + a ^ ) ] * 1 r2 r2 2 3-17 Rotations at node (3): C C 3/ K3) = C o, c) — 2 ^- 2 W32 = r 2 ( a2 + a H - c + a 7 C 1 + r 2 ( a10 + a13° + a 1 7 C } + 2 2 2 — 2 22: 9 " r 2 C a3 + 2 a 5 C + ^ 8 ° } " r 2 ( a l l + 2 a 1 4 C + 3 a 1 8 C > 3-18 b e 2 b ^ ^31 = - r 2 ( a2 + a 4 C + a 7 C } + r 2 ( a10 + a 1 3 C + a 1 7 C } + 2 c b e " r f ( a3 + 2 a 5 c + 3 a 8 c 2 ) + r 2 ( a l l + 2 a 1 4 c + 3 a 1 8 c 2 ) 3-19 t h e n % = * c C r f - r f J C a 2 + a 4 c + a ? c 2 ) + C ^ - A ) C a10 + a 1 3 C + a 1 7 ° 2 ) "' ~ °2 ( r|" + } ( a3 + 2 a 5 C + 3 a 8 c 2 ) + •; b ^ „ a _ + C ' ( r 2 " r 2 } ( a l l + 2 a i 4 c + 3 a i s c 2 ) 1 * J 3-20 - 1 5 -Shear Strains: Shear strain is normally defined as the change in angle from a right angle but since our element's sides are not initially at right angles to one another, we have to redefine the shear strains as: Define Y H The difference of the side rotations at a node. So Y l = U»12 3-21a 2 = ^ 3 -"21 3-21b 2 = W31 - Ui„„ 32 3-21c Shear strain at node (I) is : From equations 12, 13 substituted into equation 21a, yields 2 2 Y l 5 [ T T a a2 + a10 ( 1 - T" > - 2 a12 b ( 1 " T" ) + r3 r3 r3 + 3a1 5b2 ( 1 - 4 ) + ^ ( a 3 - a b + agb2) -r3 r3 r3 C a u - a 1 3b + a1 6b2) J \ 3-22 Shear strain at node (2): From equations 15 and 16 substituted into equation 21b gives 2 2 T2 = [ " ^T a 2 + a10 C 1 " V ± - 2 a a 1 2 ( 1 " h ] + r2 r2 r2 -16-+ 3a a. 2 2 1 5 C 1 - ^ ) + V C a 3 + a 4 a + a g a 2 ! + 2 | r2 r2 r2 C a l l + a 1 3 a + a 1 6 a 2 ) ] * I 3-23 Shear Strain at node (3); Substitute equations 18 and 19 into 21c yields: 2 2 = I c ( a 7 + b T ) ( a 2 + a 4c + a 7c 2) + ( \ - \ ) r2 r3 r2 r3 ( a l n + a,-.c + a 1 7 c 2 ) + c 2 C K- - K- ) ( a. + 2a Kc + 3a Qc 2) + 2 ~T r 3 r 2 C ( T + T » ( a l l + 2 a 1 4 c + ^ I S ^ ' ]4 r2 r3 3-24 Summarizing the generalized displacement vector i s : { 6 } = < u r V 1, ktU 1, ^,...^3 u,, V C , 0 9 desired degrees of freedom 2 fshear strains set to 0 J square & cubic side 3 centroidal degrees of freedom to be statically| condensed later. square & cubic side 2 Note: The nodal shear strains are a l l set to zero later, Y l = Y2 = Y3 = ° -17-Transformatdon matrix relating the degree of freedom to the polynomial coefficients i s : { « •} = I T ] { A } 18 x 1 18 x 18 18 x 1 3-25 where [ T ] = Transformation matrix { A } = a l ft The transformation matrix i s written out i n f u l l on next page (Table 3.2) where - c^ s^ are sine and cosine of angle.a for side i - B, C and A are dimensions of the element - r.^, r 2 , r ^ are the lengths of the 3 sides of the element. - 1 8 -1 CN N CD cn cn cn N CD -19-i CO U rijIcNCN CN CO ml u CO U CO | CN CVplcN CO ml n + CN IcN CN 6 ' H CN . U l CM u CO CN CO U CN co 6 i CN CN >H CN i co m i m d , CN CQ ro CN ro M CN CN |CN CN rtjl U + H >< COl CN < CN u >—»,CN CN . rill U 1 ,CN CO ! e 1 u ; CN U ! I ] C N ° CN CO U CN 1 CN CN rj CN C\P|cN CO ! ml u + CN IcN CN 5 l 54 U l CN CN PQ CNP|CN CO ml M + rH CN^ICN CN rtjl n + H rt! I j I .CN CO ml CN i t I CN. <! 6 + rH CN CN >H CN U <f|cN CN U 1 . CN CO S 1 M U l CN i I m i H I CN H + CN |CN ro ml u < . ... — CN CN CN u CN~"|CN CO ml u + CN |CN CN H | CN TABLE 3.2: TRANSFORMATION. MATRIX FOR PLANE STRESS ELEMENT (CONT'D) (3) u U >' Ift.oi 1 — 1 ! 1 1 i (A - B) 3 C 3 Z (A - B) 3 3 c 2 9 CA - B? C 3 3 3 3 c 3 ! 1 t i i l 1 — I BC 2 R 3 c 2 -BC 2 B 2C 2 ! 24 \ \ - AC c 2 - 2 2 ACI ^ 2 2 2 A C 2 r 2 C (A + B ) ;c2 (1 - 1 ) ~ 2 2 1 2 2 2 r 2 r 3 ;2 r 3 r 2 • C 2 (A + B ) C 3 Cl - 1 0 2 2 ; 2 2 2 r 2 r 3 ; r 3 r 2 I :3CA + B ) 0 2 2 2 r 2 R3 3C 4(1 - 1 ) 2 2 2 r r 3 2 1 I c 2 s 1 C 2 b 2 1 c s 2 C 2 S 2 « A C 2 S 2 A ' a — 1 — . i ! 3 1 c s j 2 2 1 2 2 c 2 s 2 c 2 s 3 c 2 s C 3 B 3 c 3 s 2 - c 2 s 3 2 B "°3 S 3 2 & 1 - 1 ' i 1 C 3 S 3 2 2 C 5 S3 c 3 s 3 .1 O 1 TABLE 3.2: TRANSFORMATION MATRIX FOR PLANE STRESS ELEMENT (CONT'D) (4) A - B h — 2 — 1(1 - B ) 1 : ^ 2. l l (1 i2 A 2) ~2~ r2 2 2 ( A W ) 1 2 2 2 r2 r3 C 3 -BC 2 r 3 CA 2r? ft Q - A ' l I e_(A_ + B_) 2 2 2 r2 r3 (A - B) CA - B) C 3 3 - B Cl - B |) B ft. S2 C2 S 3 C 2 2rt CA^ 2r? 2 2 C CA~_- Bf) „ 2 2 2 r„ h C^CA + B ) c 2 S 2 C 3 S 2 9 CA - B)' 2 2 3 B^Cl - B Z) 2 2 3 A ZC1 - AH 3S 2C-2 ( A S2 C2 -3S 3C- B S3 C3 2 2 c 2 s 2 i-c 3s 3 2 6 2 2 C3 S3 S2 C2 •S 3B S 3 C b3 C3 -22-Stiffness Matrix; The elemental stiffness matrix can be obtained from the strain energy. In plane stress the energy expression i s : U = Et / /[ U| 2 +Y1?-2 + 2 v u ^ r + 1 - V ( u ? +-V^.)2 ] d£ d? 2C 1 -v f ) ' A 2 area of element 3-26 where E i s Young's modulus, t i s the plate thickness and V" i s Poisson's ratio. Equations lb and 2a^are substituted into equation 26 and the integra-tions are carried out to yie l d the quadratic strain energy form u f = Et ^ { A } T [ K j { A } 3-27 =— A 2 C 1 -V/ ) know { <5 } = I T ] { A } then , { A } = I T 2 { 5 } 3-28 putting equation 28 i n 27 yields T - I I T ] - J - { 6 } | (1 -vf) U* = \ E ^ [[ J " 1 A [ K J [ T ] _ 1 { S } 3-29 Equate the strain energy to the external work done by the loads { p }: 2 a -yr) L * and ' { P } = i y { 5 } { 6 } 3-30 - 2 3 -then T I.KS] = I s ] - 1 I kA] I s ] " 1 3-31 11 x 11 11 x 18 18 x 18 18XJ1 where [ s ] i s the f i r s t 11 columns of [ T ] since yi = Y? 18 X 11. 18 x 18 Y 3 = 0 and the square and cubic terms of sides 2 and 3 are set to zero. 3.2.1 Integrating the Stiffness Matrix The entries of the stiffness matrix [ K &] may be determined i n closed form. F i r s t a formula for the integral / / £ m ? n de dz = F ( mf n ) 3-32 A taken over the area of the element i s obtained (4) where _ n + 1 f m + 1 . m + •! I mi i i i 3-33 2)1 When equations lb and 2a are substituted into equation 26 and we incorporate the symmetry requirement, the result i n closed form i s ; K „ = nu m. F ( nu + nu - 2, p. + p. ) + n ± n. F ( l± + l y n ± + n.. -2) + I P ± Pj F ( iru + nu, vL + p.. - 2) + 1 ± 1.. F l . . ( 1. + 1. - 2, n. + n.) ] + [ 1 - v p. 1. + vm.n. ] F ( m. + 1. - 1, p. + n. - 1 ) + [ 1 - v p. 1. + vm.n. ] F 3 i 3 i — ^ — 1 3 1 J C mi + 1. - 1, p ± + n. - 1 ) 3-34 - 2 4 -where inland p^run from 1: to 8 and l;,and n^run from 1 to 10> as defined following equation 2a. 3.2.2 Condensation o f Centroidal Degrees of Freedom: Since the centroidal displacements Uc and Vc l i e inside the element, these displacements w i l l be unaffected when the elements are joined together to represent the structure. Therefore we may solve for them before the elements are joined together, without affecting the f i n a l result. lYLuxindzing the potential energy i n one element: [ K 5 ] [ « ; ' ] = 11 x 11 11 x 1 9 x 9 3-35 hi 2 x 9 2 x 1 Evaluating: K. 11 61 + K12 82 = P. 1 3-36 hi 61 + K22S2 = P2 3-37 Solving for S i n equation 37 62 = hi ( P2 " W Equation 36 becomes 11*1 + ^ 2 ( P2 'hl^ = P l -25"-9E h C *11 " K12K22- V ^ P l - ^ 2 ^ 2 P2 62 3 " 3 8 and P = K 6^ 3-39 Therefore { P }* = ? 1 - ^ 2*22 P2 and * -1 I K J = K n - K 1 2K2 2 3 _ 4 0 9 x 9 * * where { P } and I K ] are the f i n a l load vector and stiffness matrix for the nine degree of freedom plane stress element. 3.2.3 Characteristics of the Plane Stress Element The tangential displacement along an edge i s continuous for a linear variation. The other in-plane displacement normal to each edge varies cubically along the edge. At the nodes of the element the rotation i s continuous but i t i s not continuous along the element's sides. Because of the restriction ^ = = = 0, the element i s not complete but the approximation affects the element's performance only sl i g h t l y as w i l l be illustrated later i n sameiftUYtverical appli-es cations. Inter-element compatibility (C ) i s easily achieved. 3.3 Bending Element Formulation: The Zienkiewicz nine parameter plate bending triangular element i s used with the nine parameter plane stress element. Nine degree of freedom would imply that a complete cubic be used for out of plane displacements W. However a d i f f i c u l t y arises as the f u l l cubic -26-expansion contains ten terms and any omission has to be made arbitra-r i l y . To retain a certain symmetry of appearance Cisotropy) a l l ten terms could be retained and two coefficients made equal to li m i t the number of unknowns to nine. Several such p o s s i b i l i t i e s have been investigated but a much more serious problem occurs. The transformation matrix becomes singular for certain orientations of the triangle sides. This happens for instance, when two sides of the triangle are paral l e l to the x and y axes. 0. C. Zienkiewicz pointed out that d i f f i c u l t i e s of such asymmetry can be avoided by the use of area co-ordinates (9) . R.D. Cook also pointed out that invariance could be achieved by the use of area co-ordinates (2) . The nine terms of a cubic expression are formed iby. the products of a l l possible cubic term combinations (9) i n area co-ord.; 9 W = £ •wi 3-41 1 = 1 where : ISL = shape functions and are defined as follows: N i = L l + % # 1 2 l 3 - L1 L2 " L1 L3 2 = - b x ( L 2 ! ^ + | L I L 2 3 ) + b 2 ( L 3 L 2 + 1 L L L 2 L 3 ) N 3 = - C3 ( L1 L2 + I L 1 L 2 L 3 } + °2 ( L 3 L 1 + \ L1 L2 L3 ) N4 = L2 + L2 L3 + L 2 L 1 " L2 L3 " L 2 L 1 B " b l C + i + b 3 C L1 L2 + I L1 L2 L3 } N5 N 6 = " c l ( L2 L3 + \ L1 L2 L3 } + a3 C L1 L2 + \ L1 L2 L3 } N 7 = L 3 + + - - L 3 L 2 -27-N„ = - b„ C + I 1 + b± C L 2 L 2 + j L X L 2 L 3 ) '8 ~2 N Q = - c 0 C L?Jjn + i '1 2 1 9 ~ ^2 ^ ^ 1 ' 2 L1 L2 L3 1 + c l C L2 L3 + 2 L1 L2 L3 ] 3-42 = Triangular or area co-ordinates Effectively then an incomplete cubic i n w i s used. where b i = ^2 .!-y3 b2 = ^3 - *1 h3 = y±- y2 C n = X„ -3 ~ *2 °2 x l _ x3 C3 ~ X2 *1 3-42A ** X,U Fig. 3.3 Degrees of Freedom of the Bending Element As shown i n Fig. 3.3 the nine parameters chosen to represent the element's configuration are: r 9,, CW. where x ;8W 97 and .- 8w 9x -28-Area co-ordinates C 1^, L 2 , ) are used since the formulation i s more direct and easier. C refer to f i g . 3.4) (x 3,y 3) ( x r y i ) 0 where A1 = area of triangle 3P2,etc. .© (x«,y?) t node Fig. 3.4 Area Co-ordinates ' 5 « L ^ T O T A L = h + A2 + A3 SO L X + L 2 + L 3 = 1 L3 = A 3-43 3-44 Area of Element: ( A ) 1 I 1 2 A = det X l X2 X 3 *2 ^3 evaluating and using equations 42A yields = c 2 b x - Clb2 3-45 -29-Relationship between Cartesian arid Area Co-ordinates: We know: 1 * x 1 i 1 x l X2 *3 y2 ^3 L L 3 J evaluating equation 46 yields X = X 1 L 1 + X2 L2 + X3 L-y = y l L l + y 2 L 2 + y 3 L 3 but from equation 44 L 3 = 1 - ^ - 1 ^ so equations 46b and 46c using equation 44a become: w e won't: 3-46 3-46b 3-46c 3-44a 3-4-70 3 - 4 1 b 3 - 4 8 Where { L } = second derivatives of the area co-ordinates • (L^) { X } = second derivatives of cartesian co-ordinates [ J ]• = Jacobian matrix Note: second derivatives are used since strain operator has the same derivatives. so from equation 3-48: 2 3L: 3L: '3 ^ 3 L X 3 L , = [ J ] 3-49 -30-Using the chain rule to relate the 2 co-ordinate systems 9 3 x 3 ' 3 y 3 9 L 1 3 L 1 9 x 3 L 1 9 i r 2 2 3x Z 3y Similarly: 3 = 3_x _3_ 3_y_ _J_ 3L 2 3L 2 3x 3L 2 3y - " C i ~ + b l ~ 3x x 3y then: 2 2 3 2 3 2 9 a 2 - = < £ - 4 - - 2c 0b„ + b? 8 9L^ 3x^ 2 2 3x3? 9y " c 2 9 ' 2 c b 9 2 + * 2 9 2 —2 - , — j - 2 ,b, + b, — T 9L 2 X 3x^ 1 1 3x 3y 1 9 y 2 and 3 2 a 2 ^ i B 9 + b i c o * 2 .2 _ i _ = - c c _ ^ + (^2 ) 2 1 _ . b b 3 3L 13L 2 X Z 3xZ 2 3x3y 1 2 3y 2 In matrix form, expressing equation 3-49 we obtain: r 3tf 2 3 3L: 3L L 9L 2 - C 1 C 2 " b l b 2 " C 2 b 2 - c,b, 1 1 C l b 2 + b l c 2 - 1 r 2 h 3X 2 9/ 9x9y j -31-Stxain-Displacement Relationship:. [ B J We know: { e } = { X } { W } 3-55 Where { e } = Strain vector < Bxy^ 9 { W } = < : w.fcj. = displacement (normal to plane) i = 1 1 1 •L ; = { w } { N } { X } = [ J J - 1 * - { L } from equation 48 Then equation 3-55 can be expressed as: { e } =.. I J ] { L } { w } { N;} = [ J ] _ 1 [ B ] { H } 3-55a so [ B ] = { L } { H } 3-56 Knowing L 3 = 1 - - L 2 , equations 42 can be rewritten, eliminating their dependence on : N, = 3L 2 - 2L 3 - 2 L ^ _ 2 L 2 L 2 + ^ " b 3 ( I L ? L 2 + i L l V J L 1 L 2 1 + b2 ( L 1 1 1 2 + 2 L 1 L 2 _ 2 L 1 L 2 ) - -1 2_ — C 0 ( — L7L„ + L,L„- ^ L,L? ) + cu ( CL? - L? -3 v 2 " l ^ ' 2" AT J2~' 2 "1"2 ' " ~2 v "1 1 1 2 + 2hL2 " 2 L1 L2 } I L 1 L 2 " + ! l 2 l 2 + -32-N4 2 3 = 3L* - 2L* — ZL^L^ — 2 L?^^ 2 L^L 2 N5 = - b l C L 2 3 2 3 i , i o 1 ~ 2 L 2 L 1 " L2 + 2 L 1 L 2 " 2 L?*2 1 + B 3 C 2 + 2" L 1 L 2 " - \ L^L 2) N6 - - C l T L 2 3 2 3 1 2 1 " 2 L 2 L 1 ~ L2 + 2 L 1 L 2 " L 1 L 2 } + C 3 ( 2" L 1 L + ? L1 L2 ' - \ ) N 7 = 1 - 3L2 - + 4L 2L 2 + 2 L 3 + 4 1 ^ - 3L 2 + 2L 3 N8 = - - 2 I 2 - ^ 2 ^ + 1^2 + ^ + 1 ^ ) + + b l ( L 2 " 2 V 1 " ^ 2 + K L 2 + \ L 1 L 2 + L2 3 > N 9 =- - - ^ - I L ^ + I L ^ + L ^ I L ^ ) + + C l ( L 2 - | L l L 2 - 2 L 2 + l L l L 2 + i 4 L 2 + L 2 3 ) 3-42a Now equation 3-56 can be evaluated 3l4 [ B ] = 3L; 9 L 1 9 L 2 J <C N l ' N2 ' N 3 ' N 9 ^ > o r -33-2 3 L7 3 L ? 3 1^ 8 I S 9 L: 9 L L 3 L 2 9 L L 9 L 2 3 3-56a 3L 19L 2 This yields the [ B ] matrix printed on the following page. (Table 3.4). TABLE 3.4 STR&INHDISPLACEMENT MATRIX (BENDING ELEMENT) [ B ] = 6 - 12LX + - 4 L 2 -b 3 L 2 + b 2 (2 -6L X - 3L2) - c 3 L 2 + c 2 (2-6L± - 3L2) - 4L„ 2 b l L 2 " b 3 L 2 C1 L2 " c 3 L 2 - 6 + 12L1 + 8L 2 *2 { ~ 4 + 6L X + 3L2) + b l L 2 ^ 2 C-4 + 6L X + 3L2)+ + C ; L L 2 " 4 L 1 b 3 L l " b 2 L x C 3 L 1 " C 2 L 1 6 - 4L X -12L 2 -b-L (2 - 3Lj - 6L 2 I + b 3 L l 3L X + -6L ) + -6 + 8L± + 12L2 " b 2 L l + b± (-4 + 6L, + 3^) " V l + C l (-4 + 6L 2 + + 3LX) 2 - 4L X -4 L2 •*3 ( L 1 " L 2 + i) + b 2 c - 3 L l -L 2 + J) - c 3 ( L r L 2 •* * C2 ( - 3 ^ - ^ + 2 - 4L 1 -4L 2 -b x (-3L 2 -t - L , + 1) +b. x 2" J (L 2 - ^ + 1) --3—± - c r % 3 L 2 - L i ' + 2^ + + C 3 ( L 2" L 1 + 2^ -4 + ?8Ly. v f 8 L 2 V • - ' -- b 2 ( 3 L 1 + L 2 + + 3 L 2 -c 2 ( 3 L 1 + L 2 + C 1 ( L 1 + 3 L 2 - 3 ) . -35-Stiffness Matrix From the vir t u a l work, approach, as discussed i n section 2.1.1, the stiffness matrix i s derived from equation 2-9 I K ] = j I B 7 J D J J B J d Area 2-9 area For plate bending, the e l a s t i c i t y matrix [ D ] i s I D 2 •a 1 V . 0 E t 3 12 ( 1 - v 2 ) V 1 0 3-56b 0 0 1 - v 2 Then equation 2-9 can be rewritten as: [ K ] = 3 = ^ [ R J * weight * area 2-9b i = 1 A numerical integration procedure i s used which i s exact for the element. Since strain varies linearly so does stress. This yields a.quadratic order for the stiffness. Refer to Table 3.3 where i = side no. of the element. -36-TABLE 3.3: AREA CD-ORDINATES FOR NUMERICAL INTEGRATION Side L l L2 L 3 wt 1 1 1 0 1 2 2 3 2 0 1 1 1 2 2 3 3 1 0 1 1 2 2 3 Note: Refer to figure 3.4. -37-3.4 . : Assemblingthe ,In-Elane and Bending ..Element ,.Stif finesses—; Shell, and folded plate, structures. support .their., applied loadings by a coupling of,in-plane resistance and bending resistance. Thus structural action may be represented by cabining the in-plane s t i f f -ness with, the bending stiffness. The resulting stiffness matrix i n local co-ordinates treats the in-plane and bending actions as being independent of each other (uncoupled). However when the stiffness matrix i s referenced to the global system coupling does result between the membrane and bending actions. When the elements modelling a body are assembled, coupling also exists between adjacent elements. The degrees of freedom chosen to describe the in-plane action are: { l P x } 9 = ^ 1 ^ 1 ' "2 where = 6 ( i n plane rotation) Similarly those degrees of freedom used for bending action are: { 6 b } T =/wx , e x l , e y l , w2 , ...e y 3\ 1 x 9 ' For the combined element ( in-plane and bending), the displacement vector i s to be arranged as follows: { 6 h } T = \ u i ' V i ' wx , e x l , e y l , e z l , ^ , . . . e A 1 x 18 / In terms of the force - displacement relationships: In-Plane { F } = I K ] { 5 } P P P -38-Where I J = in-plane stiffness .matrix i n local co-ordinate. 9 x 9 bending: { F b } = 1 ^ 3 { sb } where = plate bending stiff-matrix i n local co-ordinate. Combined: 18 x 1 K 1 1 2 Px 2 K 1 2 ? px ? K 1 3 * b l l 3 x 3 ^12 *bl3 1 X ] K 2 P 1 V 1 K 3 P K 2 1 P 2 x 2 K P 2 x 2 K 2 3 P 2 x 2 *b21 ^22 *b23 K 4 P l x l K 5 P L x 1 K 6 P l x l K 3 1 P 2 x 2 32 *b 2 x 2 33 K J J P 2 x 2 p *b31 *b32 *b33 K 7 P 1 x ll P L x 1 P l x l 18 x 18 -39-To further c l a r i f y what has happened, here i s how the bending stiffness matrix has been paritioned and addressed for use i n \ K 1 9 x 1 *b 11 3 x 3 *bl2 ^13 *b21 *b22 *b23 *b31 *b32 ^33 9 x 9 r W -W 1 'xl 'yi 2 x2 y2 x3 I y3 9 x 1 And the in-plane stiffness matrix i s addressed as follows { F p} = K n p l l K p l 2 K p l 3 P K J P Kp21 Kp22 Kp23 K 4 P K I P p Kp31 Kp32 Kp33 K 7 P K 8 P u 3 I K ] P -Uo-Now the resultant stiffness matrix can be used. Fran, the stiffness method technique, the stiffness matrix for a structure Is developed by surnrung the element stiffness matrices at the appropriate nodes. However prior to surtrning the element matrices, they must a l l be referenced to a ccmmon co-ordinate system (global or structure co-ordinate). 3.5 Co-ordinate Transformations Each element i n this study has associated with i t , i t s own local co-ordinate system. Each system has the same orientation with respect to the element, regardless of how the element may be orientated i n global co-ordinates the element displacement f i e l d i s expressed i n terms of local co-ordinates and as long as the displace-ment function used has a balanced representation of terms, then invarianoe w i l l be achieved even though incomplete polynomials are used (2). Fig. 3;5 illu s t r a t e s the use of local and global axes. The local degree of freedom must be related to the global co-ordinate system so that the t o t a l structure stiffness matrix can be computed by simply sunning the element stiffnesses. Let some matrix say [ X ] relate the 2 co-ordinate systems. F i g . 3.5 Co-ordinate Systems -k2-Then w I A ] V w 3-57 and where 9Y ^ J 3-57a u , v , w , e x , e y , U j 5 are the degrees of freedom i n terms of global co-ordinates. <^ir V, w, e x, • e^r 6 ^ are the degrees of freedom i n terms of the local axes for each node of the element. The [ X ] i s merely a matrix of direction cosines: [ U = XxX XxY AxZ XyX V V z XzX AzY XzZ Ebr the whole element, the transformation from global to local i s : -1*3-e . yi z l "2 18 x 1 3 x 3 <iiagbn"a L mo-t r i ces M 18 x 1 v, e. 'xi Zl 18 x 1 Or { 6 L } [ T ] { 6 G } 3-58 The elemental stiffness matrix i n local co-ordinate i s transformed to the global co-ordinate system: I K_ ] = [ T ] I K ] [ T ] 3-59 18 x 18 18 x 18 18 x 18 To find the direction cosine matrix { X ], consider the equation of the plane which passes through nodes 1, 2 and 3 of figure 3.5. The global co-ordinates of the three nodes are X. Y. Z. for i = 1, 2, 3. -kk-X - X± X - Y x Z - Z± det x2 - X l y 2 - Y L X3 " X l Y3 " Y l Z2 " Z l Z 3 " Z 1 = 0 3-60 This yields: C x - x 1 ) I c y 2 - Y1 ) c z 3 - z x ) - ( y 3 " Y l } ( z2 " Z l } ] + - ( y - y 1 ) I C x 2 - x 1 ) ( z 3 - z 1 ) - ( x 3 - x x ) ( z 2 - z x ) ] + + C z - z 1 ) I ( x 2 - x x ) ( y 3 - y x ) - ( x 3 - ) ( y 2 - y x ) = 0 Or; x I ( y 2 - Y l ] C z 3 - z± ) - ( y 3 - y x ) ( z 2 - Z ±) ] + + y I ( x 3 " x l } C z2 " z l } " ( x 2 " x i > ( z 3 - z x ) ] + + z I C x 2 - x x ) ( y 3 - Y l ) - (x 3 - x x ) t y 2 - Y l ) ] = constant 3-61 Or A X + B Y + C Z = constant 3-62 Now relate to the X - Y axes vthe direction cosines of the normal to the X - Y plane riitto axe* A = c o n P o n e n t projection (general) vector length vector length is E = " \ / A 2 + B 2 + C 2 then: AzX = A E zY B E vzZ C E 3^63 - U 5 -The clirection cosines for the x axis from node 1 to 2 are: Vector length. = a + b = 1^ Or \ = V ^ 2 - xx )2 + ( y2 - yx )2 + ( z2 - Z l ]2' Then h h h Theadirection cosines for the y axis are: The y axis i s perpendicular to both x and z so use the dot product. This yields: Xyx AxX + AyY AxY + AyZ XxZ ~ 0 Xyx Xzx + V AzY + AyZ XzZ " 0 v + y + y = 1 and the vector length i s : L = j[ ( Y 2 - y± ) c - ( z 2 - z± ) B ] 2 + [ ( x 2 - x x ) c - ( z 2 - z x) A ] 2 + [ .( x 2 - ^ ) B - ( y 2 - Y l ) A ] 2 then = _-C y, - y, .). c - ( z„ - z ) B xyx ^ ± £ ± _ ( x_ - x ) c - ( z„ - z, ) A AyY ± ± ± 1 - U 6 -and yz Cx2 - X j j B - (y 2 - Yj) 3-66 3.6 Element Dimensions (in terms of global co-ordinate) refer to f i g . 3.5. The length of side (1) has already been defined as 1^. The length of side (2), between nodes (2) and (3) i s L,: 1 2 =V(X3 " X 2 ) 2 + ^3 " Y 2 ) 2 + C z3 ' Z 2 ) 2 ' a =1 (x 2 - X 3 ) (x 2 - X x ) + (y 2 - y 3) (y 2 - y^) + ( z 2 - Z 3 ) ( z 2 - z^) ] h So b = 1 - a c = V l 2 - a 3.7 Summary of the Ccaribined Element: As a result of combining the nine parameter plane stress element with the nine parameter plate bending element, the triangular element can model the six possible movements at a node i n space, namely three translations and three rotations. The displacement tangential to each side of the element varies linearly, but the normal (membrane) displacement and plate bending vary cubically along the edges. Therefore, a l l displacements are continuous between elements. However the bending slopes normal to the edge are not compatible except at the nodes. As mentioned i n section 3.1 i f the element i s to model shells and plate structures, both displacements and slopes must be continuous for the element and from element to adjacent element. However only -47-slope continuity i s sati s f i e d at the nodes and not along the sides of the element. For slope continuity, along the-sides of the element, both the translations and the rotations have to be continuous. These sacrifices, did not hinder the elements performance to a great extent as w i l l be ill u s t r a t e d later i n some numerical examples. - 4 8 -CHAPTER 4 STRESS (XMPUTATICNS 4.1 In General: When a structure i s modelled using f i n i t e elements, the deflections of the nodes are solved from the force-deformation equation: The master stiffness matrix i s decomposed and the nodal displacements are easily computed. Before the stresses can be computed, this deflection vector should be transformed to the local system for each element and the in-plane movements and the bending movements separated. These have to be separated because the e l a s t i c i t y matrices [ D ] and the strain -displacement matrices [ B ] are different for each "type of action. Even though the maximum stress at a point i n a body i s total l y independent of any co-ordinate system used, i t i s convenient here to work with the local system for each element. Using equation 3-58 { 5 } = . [ T ] ' { fip. } where I T ] = transformation 4-1 where tX] = master stiffness matrix i n global co-ordinate master deflection vector i n global co-ordinate 18 x 1 18 x 18 18 x 1 matrix. Now the local solution vector of deflections can be broken down for each element as follows: w. Or 18 x 1 ' 9 x 1 9 x 1 L p b 4-3 Now the stresses can be computed. In general, the strains are computed from equation 2-2 where { e } = [ B ] ' { 6 } I B ] = strain-displacement - matrix 4-3a and then the stresses frcm equation 2-1 are: { a } = I D j { e } The r e s u l t a n t s t r e s s at a node i s computed by c a l c u l a t i n g the average s t r e s s of a l l the surrounding element c o n d i t i o n s . -50-4.2 In Plane S t r e s s e s : In the f o l l o w i n g , three d i f f e r e n t methods f o r approximating the s t r e s s e s from the c a l c u l a t e d displacements are presented. The f i r s t method (c o n s i s t e n t ) .uses the same stra i n - d i s p l a c e m e n t matrix as i n the s t i f f n e s s c a l c u l a t i o n . The second method (C.S.T) uses the strain-displacement matrix from the constant s t r e s s t r i a n g l e and j u s t ignores the r o t a t i o n a l degree of freedom at each node. The t h i r d method (L.S.T.) uses the l i n e a r s t r a i n t r i a n g l e s t r a i n - d i s p l a c e m e n t matrix and i n v o l v e s c a l c u l a -t i n g e f f e c t i v e mid-side node displacements, thus making use of the nodal r o t a t i o n s . 4.2.1 Consistent Formulation: The word c o n s i s t e n t i m p l i e s s o l v i n g f o r the s t r e s s e s i n the usual manner described i n the previous s e c t i o n , d e r i v i n g the strain-displacement matrix [B] that i s c o n s i s t e n t w i t h the element f o r m u l a t i o n of Sec. 3.2. Having the s o l u t i o n vector of the in-plane displacements {^p^ i n l o c a l co-ordinate, we can proceed to s o l v e f o r the s t r e s s e s anywhere w i t h i n the element. {<5 } P [D ] {e } P P 4-4 where ty {a } p {e } P 1 V 0 E V 1 0 1 - v z 0 0 1-v 2 -* fax plane s t r e s s v e c t o r ) oy L Txy the s t r a i n s {£ } P [B ] {6 } P P 4-5 - 5 1 -where I B_ J = strain-<Hsplacement matrix. P So equation 4-4 can be expressed as: = I D J I B ]"{$'} 4-6 P p p p D 3 x 1 3 x 3 3 x 3 3 x 1 The strain-displacement matrix i s formulated from the original displa-cement polynomials: u - . J a, ^ C P 1 3-Ib i = 1 V = < ai+6 5 C 3 2 a x = 1 where: { m } T = ^ 0 1 0 1 0 2 1 o) T { p } = (0 0 1 1 2 1 2 3) T { 1 } = <(0 1 0 2 1 0 3 2 1 0) { n } T = (0 0 1 0 1 2 0 1 2 3) knowing <3u Equations 4-7 can be evaluated anywhere i n the element and i n particular at the nodes resulting i n nine strains per element. -52-An alternate approach., of computing, the plane stresses was t r i e d i n an effor t to improve on the previously mentioned method. This i s based on the Linear Strain. Triangle approach. Also the constant Strain Triangle was computed to. serve as a comparison to the v a l i d i t y of the results. 4.2.2 Constant Strain Formulation; (C.S.T.) This triangular element shown i n Figure 4.1 has six degrees of freedom, two per node. The element i s rather limited because of i t s simplicity. It can only represent a constant state of strain ( and stress) across the entire element. However, this element does y i e l d stresses which can serve as an approximate check on higher order elements. I . ^a,^ Fig. 4.1 Constant Strain-Triangle Degree of Freedom -53-Again we start with the solution vector i n local co-ordinates {& } P where ' V -9 x 1 u. u. but for the constant strain triangle we only need a linear variation for the displacement polynomials since the strain-displacement operator i s only of f i r s t order, giving a constant strain variation across the element. u = a-^ + a2? + a-jC; y = a 4 + a 55 + a 6 ? 4-8 4-9 Then only six parameters are required so we w i l l use only u and V at each of the nodes. Writing equations 4-8 and 4-9 i n matrix form yields: 1 ? e 0 0 0 0 0 0 1 l n u — I V -5U-Or { U } = I a : { 9 } 4-10 where { 8 } = vector of prescribed degree of freedom I a ] = matrix of prescribed coefficients We want to find a relation which relates the u and v displacements directly to the triangle's degree of freedom. { U } = { A j { s } 4-11 where { U } = assumed displacement f i e l d { s } = actual degree of freedom of element So assume that { 6 } i s related t o { s } b y { g } : { s } = I B ]'{ 8 } 4-12 Building the [ g ] matrix : @ node CD : £ = ?=?-,_ s l = a l + V l + a3h S2 = a4 + V l + V] § node (2) K = %2 ? = ?2 s 3 = a x + a 2 ? 2 + a 3 ? 2 s 4 = a 4 + a 5 ? 2 + a 6 ? 2 -55-@ node 13} : S = 5 3 ? = e ' s 5 = a x + a ^ + a^Cg s 6 =. a 4 + a 5 ? 3 + a 6 ? 3 4-13 Putting equations 4-13 i n matrix form of equation 4-12 results i n the following [ g ] matrix. [ 6 ] 6 x 6 1 0 1 0 1 0 ? i 0 0 ?2 0 5, 0 1 0 1 0 1 0 0 c l 0 0 From equation 4-12, { 0 } can be solved for: -1 { 6 } = [ B J A ' { s } 4-14 But we want { U.} = I A J { s } 4-11 and know { U } = I a ] { 6 } 4-10 then substituting equation 4-14 into equation 4-10 yields { U } = I s J I 6 ] 1 { s } 4-15 So [ A J = . I a J I 6 J -56--1 Now we can proceed to derive the strain displacement matrix. Strains are expressed as { £ } = [ L ] { U } 4-16 where L = strain displacement operator for plane stress a n d 1 O 95 0 k 3. 3_ L 'BC 95' e — strain matrix 1 ex ey yxy Substituting equation 4-15 into equation 4-16 yields { e > = [ L J I a J J g ] 1 { s } ,-1 4-17 Then [ B ] - [ L ] [ a ] [ & ] = the strain-displacement matrix. Once the strains are computed the stresses follow from the equation. { cr } = I D 1 { e } 4-18 where D = elasticity matrix for plane stress, defined in section 4.2.1. a = stress matrix $ a X ~\ ay Txy J -57-4.2.3 Linear Strain Formulation: (L.S.T.) This element has mid-side nodes as- well as the nodes at the vertices, as shown i n figure 4.2. The element i s one higher order than the constant strain element (section 4.2.2) because not only can i t represent a constant state of strain (and stress), i t can also model a linear variation. Implicit from the t i t l e , the strains over the element are to vary linearly, so since the strain-operator for plane stress i s f i r s t order, the assumed displacement f i e l d must vary quadratically. Here area co-ordinates are used since the computation i s more direct, and ef f i c i e n t . Noting that complete quadratics i n two space require: 2 2 \ _ 4 (3) _ (V) 6 parameters for each displacement 2 ' 2 function. Then i n total 12 values of displacement must be computed for each element. A logical choice would be to use u and V at the mid-side nodes with the 3 existing nodes (refer to figure 4.2) . j j, 1,£ 2' J l3 = a r e a co-ordinate A , + A 2 + A 3 = A T Fig. 4.2 Linear-Strain Triangle -58-The area co-ordinates are defined as follows: 1 2^ ~' ^2 3^ 3 (Same as section 3.3} Noting that each shape function should be unity at one node and zero at a l l others (since the shape functions are actually interpola-tion formulae), we obtain the following shape functions: Then N l : = Z1 (2^ - 1) N 2 = = % 2 (2£2 - 1) N 3 = = &3 (2SL3 - 1) N4 = - **i*2 N 5 = - u2i3 N6 = U = 12 : ^ N. i = 1 1 s. l where Or { U J * - ^ , a,, V2, u 3, { U } = [ A ] { s } where O3} = s o l u t l o n vector of displacements 4-19 4-20 4-20a - 5 9 -where [ A r 12 x 2 %2 (2i2 l 3 ^2%3 « 2£ 3 4& S. 3 1 - 1 1 - I I - 1} l2 (2£2 l 3 ®l2 4 £1^2 0 0 4<3\ -1] -1) -1) Area of the Element: A 1 1 1 A 1 **• l l *2 X 3 V l V 3 V3 let | det X l " X 3 X2 " X3 Yl " Y 3 Y2 y 3 a^ = x„ - x, a« — x, - x 3 2 2 "1 "3 1 a3 " *2 ~ X Bl = V2 b2 = V3 b 3 = y x - y. - 6 0 -then 2 I a 2 ^ i ~ b 2 a l ^ Relating area ccHordinates to Cartesian (36-ordinates: 1 1 1 xl X2 X ^2 y 4-22 Expanding equation 4-22 and making use of the fact that £1 + l 2 + l 3 = 1 Or a3 l SL1 - i2 4-23 gives: x = x 1 £ 1 + x2l2 + x 3 (1 - £ X - A ). y = y x ^ + Y 2 * 2 + ^3 ( 1 - £ i _ V 4-24 We want { L }= [ J J { X } 4-25 where { L } = { X } = [ J ] = f i r s t derivative of area co-ordinates ( I. ) 1 f i r s t derivative of Cartesian co-ordinates Jacobian matrix Note : F i r s t derivatives are used because the strain operator i s f i r s t order (contains only f i r s t derivatives). -61-From equation 4-25 3 A, 3A, J r J ] _3_ 3x _3_ [ ^3y 4-26 Using the chain rule to relate the two co-ordinates systems: 3 _ 3x 3 • • _3y_ _3_ 9&2 9 A 1 3 X 9 Y 9x and: b2 ~ 9y 4-2 7a 9 A, 3x _3_ 3y 2 3 £ 2 3x 3£ „ 3y - a — + b — 9x ± 3y 4-2 7b then expressing equation 4-2 7a and 4-2 7b i n the form of equation 4-26, yields the following Jacobian matrix. I J ] = a2 " b 2 - a, .Ss inverse i s : u2 *2 where A = area of the element - 6 2 -Stxain-Displaoement Relationship: We know:. { e } - { X } I U } 4-28 where U = displacement (in-plane), { e } = vector of strains then from equation 4-25: So { X } = [ J J 1 { L } ,-1 { , £ • - } = I J J { L } { U } ,-1 4-29 4-30 = I • J J ' ( L } I A ] { s } 4-30a Therefore [ B J - { L } [ A J i s the strain-displacement matrix, where L = strain operator = Evaluating equation 4-30 a [ J ] 1 I L ] = i 2 A 1 0 3 3£. _3_ 0 3x 0 3 3y _JL _9_ 3y 3x •d • + b2 3£ 2 "I 1 - s "1 ~ K 9 ^ 3£ 2/ ( b i — + b2 — \ 9A-, yields: -63-^ L2t± - 11 A 2 C2^ 2 - II £2 .(2£2 - 1) Cl - A2 -Aj.) [2 (1 - £ 2 - A 1 ) - 1] [ A 12 x 2 4 V 2 a - A 2 - J ^ I 1 2 a - ^ - j » 2 ) - 1 ] o 4*1*2 4A2 (1 - * 2 - ^ 4 C 1 - ^ - A 2 ) ^ 4£ 2 ( 1 - * 2- i±) 4 ( 1 - l± - £2) £x The resulting [ B ] premultiplied by [ J ] for use i n equation 4 - 30a i s shown i n table 4.1 on the next page. TABLE 4.1: /graMNrDISElraQMMT. MATRIX FOR L . S . T . ^ ( 4 ^ - 1 ] 0 a l ^ 4 j i l " 1 5 0 ^ ( 4^ - 1 J b]_ C 4A 1 - 1 ) fc>2 ( 4 i 2 - 1 } 0 C 4£ 2 - 1 ) 0 a2 ^ 4 i l2 ~ 1 } b 2 C 4A2 - 1 ) b ± ( 4A x + 4£ 2 - 3 ) + +fc>2 ( 4£ 2 + 4 ^ - 3 ) 0 a ± ( 4 ^ + 4A2 - 3 ) + +3^ ( 4*2 + 4H± -.3 ) 0 a x C 4^ + 4£ 2 - 3) + t a 2 ( 4£ 2 + 4 ^ - 3 ) b x ( 4A1 + 4«-2 - 3 ) + fb 2 ( 4£ 2 + 4^ - 3 ) 4 ( + b ^ ) 0 4 ( a ^ 2 + a ^ ) 0 4 ( a ^ 2 + ^ \ ) 4 ( b ^ 2 + b 2 ^ ) b x ( - 4 ^ ) + +b2 ( 4 - 4^ - 8^ } 0 a x C - 4^ ) + a 2 ( 4 - 4 ^ - 8 ^ ) 0 a x ( - 4 ^ ) + + a 2 C 4 - 4 - 8 ^ ) b x (-4^ 2 ) + + b 2 ( 4 - 4 « x - 8 £ 2 ) b x ( 4 - 8 ^ - 4 ^ ) + + b 2 ( -4 ^ ) 0 ( 4 - 8 ^ - 4 ^ ) + + . ( - 4 ^ ) 0 a x ( 4 - 8^ - 4^) + + a 2 ( - 4 ^ ) b l ( 4 " 8 \ " 4 ^ ) + + b 2 ( - 4 ^ ) Where A = Area of Element - 6 5 -Now the stresses can be computed: { cr } = I D ] { e } 4-31 Or { a } = I D ] I B ] * [ S ) 4-32 where I B J * = J J ] " 1 I B ] and {, s 3 = the u and y displacements at the corner and midside nodes. Mid-Side Node Displacements: We are not through yet because u and V of the mid-^ side nodes have not been defined. From the u and V displacements of the three vertices, we must somehow derive reasonable mid-side displacements. First resolve the cartesian u and V displacements into tangential and normal displacements at each vertex. The tangential displacements are shown in figure 4.3. Fig. 4.3 Tangential Displacements =66-The relations used to resolve the cartesian displacements into tangential components at the vertices are; @ Node CD : u t l 2 = U ; L 4-33 u t l 3 = - S^sin g 1 - cos ^ 4-34 @ Node (2) : Ut21 = ^2 . ' 4-35 Ut23 = V2 s i n 2^ ~ U2 0 0 6 B2 4 - 3 6 @ Node (3) : Ut32 = W3 s i n 62 " "3 0 0 3 B2 4 - 3 7 ut31 = - s i n % ~ u 3 ° ° s &2 4-37. Now define the tangential displacement of a mid-side node to be the average of i t s end node tangential displacements. u t 4 = u t l 2 + ut21 4-39 u ^ = ut23 + ut32 4-40 u t 6 = Ut31 + U t l 3 4-41 The normal displacements are shown on figure 4.4. 1SZ3 stoe NO. Fig. 4.4 Normal Displacements The relations used to resolve the cartesian displacements into normal components at the vertices are: @ Node (1) u ,_ = ViL nl2 1 U n l 3 = U l S 2 J 1 H ~ V l 0 0 3 4-42 4-43 @ Node (2) Un21 = V2 Un23 =~ U2 s i j l 62 " W2 0 0 3 62 4-44 4-45 @ Node (3) u h31 u 3 s i n g i ~ V ' 3 ° ° s 6 X Un32 = " "3 S i n g2 " V3 0 0 3 4-46 4-47 -68-Define the normal displacement of the mid-side node by an interpol--ation of Her*ni+iafi cubic CshapeJ functions. Un4 = unl2 C 1 " 3 ^ + 2 ^ 1 + s l * i & " 2 ^ + ^ + un21 ^ ~ + + s ^ ( g 3 -l2\ 4-48 = in-plane rotation at node i s^ = length of side i 5 = a running dimensionless parameter, varying linearly along an edge from 0 at starting node to 1 at end node. Therefore = ^ of the mid-point of a side. As can be seen, the Hermit i-an polynomials are a set of shape functions for an element's side at the ends of which the slopes and values of the normal displacements are used as variables. Simplifying equation 4-48 yields: Un4 = U n l 2 + ^ l + V l " S1 W2 4 ~ 4 8 Similarly: U n 5 = u n 2 3 + + ^32 - 4 ~ 4 9 Un6 = ^ J l + f3^3 + u n l 3 " S 3 W 1 4 " 5 0 - 6 9 -The f i n a l step before these mid-side displacements can be used to compute the strains i s to transform them back, to the cartesian co-ordinate system. The following equations perform the task: @ Node (1) uL x = u t l 2 4-51 VL X = V2 4-52 @ Node (2): uL 2 = U t 2 1 4-53 VL 2 = u n 2 1 4-54 @ Node (3) ^ 3 = " Ut32 0 0 3 S2 " Ut31 S l n B31 4 " 5 5 V L 3 = Ut32 s i s i g2 " Ut31 0 0 3 631 4 " 5 6 @ Node (4): uL 4 = u t 4 4_ 5 7 ^ 4 = u n 4 4-58 @ Node C5) ^ 5 = " Un5 3 i n h ~ u t 5 0 0 3 62 4 " 5 9 YL 5 = - u^cos g 2 + u t 5 sin 6 2 4-60 -70-§ Node C6) ^ 6 = Un6 s 5 s i g l " u t 6 0 0 3 §1 V L 6 =~ un6 0 0 3 h ~ u t 6 s h l h 4-61 4-62 4.3 Bending Stresses: The bending stresses are computed by applying the equation: { % } = I \ 11 % J ^ b } where { a, } = b i m ^ L xy = in.- k. iP_ in. of length Cassumed to be valid over the whole element) Note and [ ] i s defined by equation 3-56b [ ] i s defined i n Table 3.4 { 6^ } = solution vector of displacement s. O e W -to £ 1 3 . 3 . 3 ) V / * 3 J A^atA. the stresses M M < J be e v a l u a t e d a/*jinhere. uMW -the, ele-»v\eT>t but IA paT--t»cuU»r ot i t e nodes as used herein. CHAPTER 5 BEAM STIFFENER ELEMENT A beam stiffener element i s used i n conjunction with a plate. The beam element strengthens the plate, increasing the flexural r i g i d i t y of the system. Deformations caused by bending are considered and as i n section 3.2 for f i n i t e element formulation the assumed displacement fields are substituted into the strain energy and an element stiffness matrix i s obtained. An unsymmetrical section implies that the centroid of the cross section and shear centre do not coincide. Consider an " L " shaped section which acts as a stiffener for the plate shown below y,y Note ; Right-handed system i s used throughout. ONE FINITE ELEMENT CF A PLATE > x,u node at mid-ht. of sldb. ^ L-SECTION BEAM STIFFENER ELEMENT Fig. 5.1 Beam Stiffener Element Define: - bending about y-axis i n the x-z plane to be the strong action of the stiffener. - bending about the z-axis i n the x-y plane to be the weak action of the stiffener. -72-Let e = vertical distance from the shear centre of the beam to the neutral plane of the plate section. If possible we want the beam stiffener element to be compatible with the adjoining plate i t i s stiffening. Consider, f i r s t the bending of a symmetric section. The bending of a doubly unsymmetric section i s merely a change i n the stiffness matrix of the symmetric bending case. 5.1 Symmetric Bending: Strong Direction: - consider bending i n the ver t i c a l (z - x ) plane about the y - axis. Fig. 5.2 Beam S t i f f ener Geometry (Strong Direction) -73-To describe the beam's displaced position u, w and 6 y are used at each end of the element. From geometry: dw UQ = u c - e _ a 5-1 dx Wg = w c 5-2 For a compatible element, since w of plate i s a cubic variation, then w of beam must be cubic. W c = a l + aZ* + a 3 ^ + a4-# 5 " 3 l e t u, = a 5 + agx + a ? / 5-4 Then from equation 5-1 ^ = a 5 + a g ^ + a y - e ( a 2 + a 32x + a 43x 2) = (a 5 - ea 2) + (a g - e2a3) x + (s^ - e3a4) x 2 5-5 but i f u of beam i s to be continuous with u of plate, i t must have a linear variation. 2 Then the x term of equation 5-5 must vanish: (a ? - e3a4) = 0 Or a ? = 3a 4 e 5-6 -7U-.Relate the Degree.of Freedom teethe polynanial coefficients;. { 5 } = I T ] ' [ A } 5-7 where { <5 } = < °yi & { A } = w„ 2 *2 *3 < V d5 and I T J = transformation matrix Knowing: dw dx a^ - 2a^x - 3a^xz 5-8 = (a 5 - ea 2) + (a g - 2ea3) x 2 3 wB = a 1 + a 2x + 3-^. + a 4x 5-9 and @ node (1) x = 0 ; @ node (2) x = £ -75-Then: I T ] = 0 -e Q 0 1 Q 1 0 0 0 0 o 0 -1 0 0 0 o Q -e -2e£ 0 1 t 1 £ £ 2 £ 3 0 0 0 -1 -2£ 2 -3£ 0 0 NcwAthe stiffness matrix i n terms of the polynomial coefficients can be developed from the strain energy (U): U= ^ / t w " )2dx + ^ A u ' )2 & 2 ° 2 ° 5-11 Where the f i r s t term of equation 5-11 i s the strain enrgy stored in the beam due to pure bending arid v . ' ; ' the second term i s due to axial deformation Using expression 5-3 and 5-4 and 5-6, equation 5-11 becomes: U = — [ 4a 2£ + 12a 3a 4£ 2 + 12a 2£ 3 J + = 2a 2£ EI + 6a 3a 4£ 2 EI + 6a 2£ 3 EI + a 2£ — + 3a ga 4e£ 2 EA + 4- aJ- «2„3 ^ 2 + 6a. e £ — 2 = E I 2a 2£ I + 6a 3a 4£ 2I + a 2 £ 3 C6I + 6e2A } + a 2£ ^ . + 6~ 2 + 3ea ga 4£ A ] 5-12 - 7 6 -But we also know that (quadratic form of Ul U = j { A }T I 1^ 3 { A } 5-13 So writing equation 5-12 i n the form of equation 5-13 and making use of symmetry yields: t K A ] = 2E 3Z21 3« 2I £ 3C6l+6e 2A) 3 „2 A 2 -3 2, y e n A HA 9 — Let I. = I + e A , [ Kft ] becomes: 2 6£ 2I 12£ 3I 0 3e£ 2A 1 A, £A -77-The stiffness; matrix I K^ J i n terms, of the polynomial coefficients can be transformed to be expressed In terms of the nodal degree of freedom as follows: Know { A } = [ T J " 1 { 6 } 5-10 and U = j { A }T I K^ J { A } 5-13 then substituting equation 5-10 into equation 5-13 U = \ { <5 }T I T J " 1 T [ K^] [ T ] _ 1 { S } = |' { fi }T I K 5 J { 6 } 5-14 where [ K g ] = [ T ] - 1 ?: [ J [ T ] _ 1 i s the local stiffness matrix of the beam s t i f f ener element. The resulting K matrix i s : L 2 0 12(r 24e 2) Symn et r i c - e L 2 2 2 -6L(r +&\ 2 2 2 4IT Cr +e ) - L 2 0 eL 2 L 2 0 -12(r 2+e 2l 6LCr 2+e 2l 0 12(r 2+e 2l eL 2 -6LCr 2-^ 2] 2 2 2 2 i r ( r +e^] - e L 2 6LCr 2+e 2l lL 2Cr 2-te 2l -78-' A A ~ X-SECTIONAL A R E A . f MA node. X,u Fig. 5.3 BEAM STIFFENER CDEGREE OF FREEDOM) STRONG DIRECTION -79-Vfeak direction bending i n the horizontal (x - y) plane about z -axis. The formulation i s analogous to that of the strong direction formula-tion except the x and y axes are used instead of the x and z axes. The degree of freedom used to describe the deformed position of the beam are u, y and © at each end of the stiffener. Refer to Fig. 5.4 z 3 and Fig. 5.5. B r . c , N A 2. . . ( »- XtUL Fig. 5.4 Beam Stiffener Geometry ( Weak Direction) From Geometry: dV % = u c " e g - H 5 " 1 5 dx ^ B ~ * c A complete cubic i n y i s used since the variation along the plate i s cubic. 2 3 V = a nx + a„x + aoc + a.x 5-16 C 1 2 3 4 - 8 0 -and so on;, as-, i n the strong direction derivation. v, X. u Fig. 5.5 Beam S t i f f ener (Degree of Freedom Weak Direction) The same stiffness matrix i s obtained as given for the strong direc-tion formulation, but velo-hve + h e »° -r»g-5. 5. Torsion: So far we have not considered twisting Crotation) of the section. Refer to figure 5.6. Where $ = Fig. 5.6 Beam S t i f f ener (Torsion) -81-Slnce we do not have continuity of © between nodes Cl) and (2), there i s no point i n striving for a compatible, element i n torsion. But the twist at the nodes w i l l be compatible If we use a linear variation for tf>. * = C l - £ l +*2 C - £ l Knowing, the strain energy for Torsion (U) i s 5-17 U = G J e f f f\$')2 dx 5-18 The stiffness i s : I Krp ] G J eff 1 -1 -1 1 Where J - polar moment of inertia 1 3 = -=%hb for thin sections 3A Now we can combine the stiffnesses together to form the overall XU e 4 2 12 x 12 matrix for one beam element Fig. 5.7 Resultant Beam Stiffener C 12 degrees of freedom) I CO ro I -83-5.2 Unsymmetric.Bending: (reference 81 Consider Bending of a beam by couples ity and mz acting i n two ar b i t r a r i l y chosen perpendicular axial planes zx and yx. Refer to figure below: Fig. 5.8 Beam Stiffener Subjected to Couples Bending i n zx plane: Assume that the magnitudes of the couples are such that bending occurs i n the zx plane, so that the neutral axis i n each cross section i s para l l e l to the y axis. The radius of curvature due to the bending i s r^, and the bending stresses w i l l be: crx = E e 5-19 and e = r - 5-20 ^z Ez Therefore cr = — 5-21 z -8U-Then the bending couples, can be expressed as: M - = /zcr dA = ^ j r 5-22 Y A x r z M = / ycr dA = E I y z 5-23 A r z Bending i n xy plane: If the magnitudes of the couples cause bending i n the xy plane, then you get the analogous equations: Bending stresses: av = E e and e = Y Therefore °x = F- 5 " 2 4 Y Bending couples: M = / Y dA = E I z 5-25 A r Y M = / z f f dA = E I y z 5-26 x A Note: I = / y 2 dA rY I , = / z 2 dA 5-28 / yz dA A - 8 5 -and I r 5-28 Bending i n Both xy and xz plane: (coupled action) In the general case, the beam deflects.:in both planes. The relations between the bending moments and curvatures are obtained by combining Mie equations for the uncoupled cases. VT "RT My = Z x + yz 5-29 r r z y M z = ^ + E I y z 5-30 r r y z Since the beam stiffener w i l l always be attached to a plate, then the bending (couples) can always be chosen to act i n the xz plane, so M =0. z Therefore from equation 5-30 Or I I _z yz r r y z rY='— r z 5 " 3 1 I yz then equation 5-29 becomes: - 8 6 -M « E, y . i - , . i _Y_ + _vz_ E_ r 5-32 For pure bending, the total strain energy i s : U = V z 5-33 Note: 2 •> _ 1 . J q _ dx d z „ , r. 3Z - — ; de • = - w " dx 5-34 z r dx z Also e = ^ 5-35 EI y Then equation 5-33 can be expressed as: u = _Y_ 5-36 2EI y - 8 7 -If the strain energy due to axial deformations i s also considered, the strain energy for an incremental length, dx i s : U = 2EI 2 „ 2 , . EA /M^ dx + - / (u' ) dx 5-37 Incorporating equation 5-32 and 5-34 U = E 21 H I , EA ^ , , x 2 . 21 -,2 I I - I y z yz *rCw")2 ^ + — r (u 1 ) 2 dx v c 2 ^ 5-38 Note: I f a symmetric section i s evaluated using equation 5-38, I = 0 , then the equation becomes yz a z U= ' ( w " ) 2 d x + — ^ (u 1 ) 2 dx o C O o C 21 y which i s the same as equation 5-11 i i z ¥ o ) We know W c ~ WB = a l + ?2X + a 3 x + a 4 x ~ 2 u c = a r + a^x + a, *5 + a 6 x + a 7 x So proceeding as was. done for the symmetric case, the stiffness matrix can be developed. -88-The resulting matrix i s found.to be analogous to the one given for a symmetric section i n table 5.1 where the values Re and rE become: Re rE (e 2 + U 2! 2 2 (Ei + U>) and U 2 = 2 2 r , AT - I z Y z yz 2 ,2 * * * u = 2 2 r AI - I ^ z z yz 2 .2 z E = distance from z axis to shear centre (e y -> horizontal distance) e = distance from center of gravity of slab to shear centre (e z -* ve r t i c a l distance) -89-CHAPTER 6 NUMERICAL APPLICATIONS 6.L • Constant Stress A p p l i c a t i o n s ; The nine degree of freedom plane s t r e s s element developed i n s e c t i o n 3.2 i s f i r s t t e s t e d under simple s t r e s s c o n d i t i o n s . This i s to see how much the element's incompleteness hinders i t s performance. Because of the nodal shear s t r a i n c o n s t r a i n t s , the element w i l l not be able to model the true s t r e s s s t a t e e x a c t l y but perhaps i t w i l l be able to make a good or reasonable approximation to i t . A square p l a t e supported as shown i n f i g u r e 6.1 i s subjected to a constant shear s t r e s s , a constant normal s t r e s s and a l i n e a r l y v a r y i n g normal s t r e s s (constant moment). The p l a t e i s modelled by two f i n i t e elements and dimensionless u n i t s are used throughout. R o t a t i o n a l degrees of freedom are allowed at a l l nodes but because the nodal shear s t r a i n s ( r o t a t i o n s ) are constrained f o r each element ( s e c t i o n 3.2) there may be some d i s -crepancy here. The constant shear s t r e s s s t a t e i s simulated by lo a d i n g the p l a t e as shown i n f i g u r e 6.1. Table 6.1 presents the r e s u l t i n g d e f l e c t i o n s and the exact values are als o t a b u l a t e d d i r e c t l y under these values. The u displacements are the same as the exact and the v displacements are only about 7% i n e r r o r of the exact values. At the f r e e end of the p l a t e , (nodes 3 and 4) the r o t a t i o n a l r e s u l t s are reasonable; only about 7% e r r o r . The c a n t i l e v e r e d p l a t e using the same g r i d and boundary c o n d i t i o n s as i n the constant shear s t r e s s l o a d i n g ( f i g u r e 611) i s used to model constant normal s t r e s s . The lo a d i n g i s i l l u s t r a t e d i n f i g u r e 6.2. Table 6.2 -90-compares the r e s u l t i n g d e f l e c t i o n s from the f i n i t e element i d e a l i z a t i o n to the exact ones. The u displacements are 16% i n e r r o r at the fr e e end but the v displacements and the r o t a t i o n s are much greater i n e r r o r . A l i n e a r l y v a r y i n g normal s t r e s s (constant moment) i s simulated by the loading shown i n f i g u r e 6.3. Table 6.3 presents the displacements using the f i n i t e elements and a l s o the exact displacements. The u and v displacements are about 26% i n e r r o r at the fr e e end. The r e l a t i v e e r r o r i n the r o t a t i o n s at the f r e e ends are 21% (node '3) and 4% (node 4 ) . In g e n e r a l , the element i s unable to model these s t r e s s s t a t e s e x a c t l y as was expected due to i t s incompleteness which i s due to c o n s t r a i n i n g the nodal shear s t r a i n s ( r o t a t i o n s ) . The displacements and r o t a t i o n s i n the constant shear s t r e s s case were only s l i g h t l y i n e r r o r of the exact values. However, i n the other two cases w i t h the exception of the u displacements, the p r e d i c t i o n s from the element were r e l a t i v e l y poor. I t i s ^ i n t e r e s t i n g to note the s t r a i n energy r e s u l t s from these t e s t s . That i s , i n the first.-.two c o n s t a n t . s t r e s s t e s t s , the s t r a i n energy e r r o r was only 7.5 and 3.1%, r e s p e c t i v e l y , whereas f o r the l a s t l i n e a r s t r e s s case, i t was much higher at 32%. In examples to f o l l o w , we s h a l l see to what extent the element's incompleteness hinders i t s performance. - u i _ CONSTANT S T R E S S APPLICATIONS d» MENS I ON LESS UNITS: | e= i.o L - i.o t = 10 -y=0.3 o .v FREE AT ALL NODES j SIMULATING 0.5 FIGURE er- CONSTANT SHEAR STRESS *-o.5 SIMULPO-IIHG or * - 1.0 12. 12. FIGURE 6,2: CONSTANT NORMAL STRESS SIMULATING M= LO FIGURE 3 ' CONSTANT SENDING MOMENT •92-CONSTANT STRESS APPLICATION TABLE 6.1: DEFLECTIONS FOR CONSTANT SHEAR STRESS POINT - ^ - ^ D I S P L . POINT (1) POINT (2) POINT (3) POINT (4) U 0.0000 0.0000 0.0000 0.0000 V 0.0000 0.0000 2.4044 2.4044 U) 1.2022 1.2022 1.2022 1.2022 ^EX 0.0000 0.0000 0.0000 0.0000 V E X 0.0000 0.0000 2.6000 2.6000 EX 1.3000 1.3000 1.3000 1.3000 NOTE: EX = EXACT VALUE S t r a i n Energy, U = 1.2022, U ' = 1.3000 TABLE 6.2: DEFLECTIONS FOR CONSTANT NORMAL STRESS P O I N T ^ ^ ^ ^ ^ - ^ D I S P L . POINT (1) POINT (2) POINT (3) POINT (4) U 0.0000 0.0000 0.8353., 0.8353 V 0.1353 0.0000 0.1353' 0.0000 : Si 0.1690 -0.6355 -0.6355 0.1690 U E X 0.0000 0.0000 1.0000 1.0000 V E X 0.3000 0.0000 0.3000 0.0000 EX 0.0000 0.0000 0.0000 0.0000 U = 0.4847, U_ v = 0.5000 -93-TABLE 6.3: DEFLECTION FOR CONSTANT MOMENT POINT ^ ^ ^ ^ ^ ^ ^ D I S P L . POINT (1) POINT (2) POINT (3) POINT (4) U 0.0000 0.0000 4.4241 -4.4241 V -0.1911 0.0000 4.4241 4.2329 0) -2.6446 -0.6625 9.5108 11.4929 U E X 0.0000 0.0000 6.0000 -6.0000 V E X 0.0000 0.0000 6.0000 6.0000 " E X 0.0000 0.0000 12.0000 12.0000 U = 4.0934, U „ v . = 6.0000 -9h-6.2 C a n t i l e v e r Beam Problem: Although the element developed i n Chapter 3 was meant to model plate"and s h e l l type s t r u c t u r e s , i t was f e l t that i t would be b e n e f i c i a l to compare the nine degree, of freedom plane s t r e s s element ( S e c t i o n 3.2) to other common elements i n a f a m i l i a r plane s t r e s s a p p l i c a t i o n . The w e l l known c a n t i l e v e r beam was s e l e c t e d to be modelled. The beam has u n i t thickness and i s loaded by lumping the p a r a b o l i c a l l y v a r y i n g shear s t r e s s at the end nodes as loads. The m a t e r i a l p r o p e r t i e s and the various gridworks used i n the a n a l y s i s are shown i n f i g u r e 6.4. The boundary c o n d i t i o n s at the c a n t i l e v e r e d end are f i x e d e n t i r e l y . Since the nodal r o t a t i o n and IX-displacement here are f i x e d , then the ^d i s p l a c e m e n t between the nodes (along the elements', side) i s constrained to be zero a l s o . The r e s u l t s are compared w i t h the constant s t r a i n t r i a n g l e (C.S.T.) and the l i n e a r s t r a i n t r i a n g l e (L.S.T.). Table 6.4 presents the t i p (end) d e f l e c t i o n obtained from the C.S.T., L.S.T., as w e l l as the nine d.o.f. element ( s e c t i o n 3.2). From the t a b l e , i t appears that the C.S.T. has a higher convergence r a t e but the nine d.o.f. element i s more accurate f o r a given g r i d of elements. The L.S.T. d e f l e c t i o n s are s u p e r i o r to both the other two elements. For a g r i d of four the three element types y i e l d reasonably accurate d e f l e c t i o n s . The exact d e f l e c t i o n i s computed from f l e x u r a l theory., Figure 6.4 i l l u s t r a t e s the performance of the three elements as more g r i d refinements are used. A l l three types of elements appear to be converging at a reasonable r a t e to the exact t i p d e f l e c t i o n . However, f o r r e l a t i v e l y coarse g r i d s , the nine d.o.f. element -95-i s f a r more accurate than the C.S.T. element. R e f e r r i n g back to t a b l e 6.4, the s t r e s s e s obtained from the three types of elements are presented f o r the various g r i d , s i z e s . A l l the s t r e s s e s appear to be reasonably accurate and are converging i t appears to the exact value. The nine d.o.f. element s t r e s s e s are more accurate than the C.S.T. s t r e s s e s . The L.S.T. s t r e s s e s are b e t t e r than the other two element types, however, i t r e q u i r e s f a r more d.o.f. f o r a given gridwork than the other two element types. In general, the nine d.o.f. element performed b e t t e r than the constant s t r a i n element and not q u i t e as good as the l i n e a r s t r a i n element. -96-L = 48" -*-X E = 30,000 kS | V = O.E5 t = l.O «K-PARA&OL ICALLY VARYING END SHEAR TOTAL LOAD P= 4o K. 20 K 20 K 6-667 K 26467 K 6.647K GRID2. GRIQ4 F\GURE <o-A-l CANTILEVER BEAM ( L0ADIN6 < GR»Os) CANTILEVER BEAM TABLE 6.4: TIP DEFLECTION AND NORMAL STRESS TIP DEFLECTION (IN) ***NORMAL STRESS (K/IN 2) FINITE ELEMENT GRID CONSTANT STRAIN TRIANGLE (C.S.T.) NINE D.O.F. PLANE STRESS ELEMENT* LINEAR STRAIN TRIANGLE (L.S.T.) C.S.T. PRESENT;-. TRIANGLE* L.S.T. 1 0.0909 0.2302 2 0.1988 0.2983 0.3550 43.28 47.855 59.145 4 0.3115 0.3291 0.3556 53.51 55.774 60.024 EXACT ** ! • j 0.3558 i j i 60 FINITE ELEMENT ELEMENT NO. OF DEGREES TYPE GRID'"-\ OF FREEDOM C.S.T. ,1 16 9 D.O.F. 24 C.S.T. ? 50 9 D.O.F. 74 L.S.T. 160 C.S.T. 162 9 D.O.F. \\ 242 L.S.T. 576 NOTE: *''fRefe-rs tg the -9 D.O.F. p l a n e str.ess t r i a n g l e d e r i v e d i n s e c t i o n 3.2 h e r e i n ( u s i n g C S T " s t r e s s " c a l c u l a t i o n s ) . ** E x a c t s o l u t i o n o b t a i n e d from f l e x u r a l t h e o r y . *** Normal s t r e s s a x (a) x = 12" and y = 6.0" C A N T I L E V E R BEAM, P R O B L E M * EXACT TIP DEFLECTION = 0 . 3 5 5 8 IN. L S.T-5 % ERROR N O T E > * FROM F L E X U R A L THEORY h. FOR C.S.T. ELEMENT. 0 FOR 3 D.O.F. E L E M E N T ( R E F E R TO SEC. 3.2J o FOR L..ST- E L E M E N T . I CO I 4 h loo 200 30O 40© TOTAL NO. OF DEGREES OF FREEDOM soo 575 FIGURE 6.5 '• T I P D E F L E C T I O N VS NO. OF DEGREES OF FREEDOM -99-6.3 Parabolically Loaded Square plate: This example serves to test the new plane stress element derived i n section 3.2. Here the elements only act i n plane stress because they are loaded i n their plane. No bending stresses are induced. The problem i s a square plate loaded on two opposite sides by a paraboli-c a l l y distributed normal stress. The other two sides are free. The loading and a typical grid layout are shown i n figure 6.6, making use of symmetry only one quarter of the plate i s modelled. Various gridworks are used and the results are compared to an exact solution (3) . The load vector used i s a consistent one based on the v i r t u a l work of the parabolic distribution times the cubic distribution for the edge displacement. Table 6.5 ill u s t r a t e s a comparison of the various deflections and strain energy with the exact solution for various gridworks. With each refinement i n grid, the deflections u^, u c, Vc, V D appear to be converging monotonically to some values s l i g h t -l y i n error of the exact values. The reason for this apparent error i s due to the fact that the element i s incomplete and the nodal shear strains are constrained, making the element s t i f f e r . The strain energy i s also converging i n a similar manner, monotonically to a value 13% i n error of the exact. Figure 6.7/. il l u s t r a t e s the manner in which the strain energy converges with each refined gridwork. Similarly i n figure 6.'8^ the end deflection i n the direction of -100-the applied, load,ia„,cx)nyerging-but...to.a-value„rouc^ i n error of.the exact. The resultant, stress at. a. node was; computed by calculating the average stress of all...the surrounding element contributions. Some characteris-t i c or typical stresses are compared i n table 6.6 with the exact for various, grid refinements.. The linear strain (L.S.T.), constant strain CC.S.T.) and consistent formulation (section 4.2.1) are computed. As illustrated a l l three values compare to the exact with only slight error. For a l l gridworks N^g , N^, and are exactly the same for the three stress computations. In general the C.S.T. values were better where the three stress values differed. Figure 6 .9 shows the rapid convergence of and N for refined gridworks to values only s l i g h t l y i n error (2%) of the exact. The reason again i s due to the element's incompleteness and the nodal shear strain constraints making the element s t i f f e r , therefore inhibiting i t from absorbing as much strain energy as i t would i f i t were complete and no constraints introduced. The variation of with grid refinements i s il l u s t r a t e d i n figure 6.10. Here the consistent formulation (section 4.2.1) stresses are very poor and converge slowly to a value approximately 30% i n error of the exact. The C.S.T. stresses however converge rapidly and are only s l i g h t l y i n error (for 10 x 10 gridwork 4%). The L.S.T. values on the other hand converge more slowly than the C.S.T. values and for a gridwork of ten are .\%ti.M error with the exact. This i s 101-; probably •. due to:. the. displacement...f i e l d ...1 imitations, put ..on. the mid-side nodes (section 4.2.31. . Figure .6.11 again i l l u s t r a t e s the superiority of the C.S.T.. stresses! for.convergence and relatively small error for N ^ vs grid.size. Here the consistent formulation stresses appear better but they have converged to a value i n error of the exact, whereas the C.S.T. are s t i l l converging. The variation of N with grid size i s plotted on figure 6.1Z The L.S.T., C.S.T. and consistent formulation yielded identical results for each gridsize. Convergence again i s rapid and i s only about 3% i n error of the exact value. In general the deflection values compare closely with the exact ones and the stresses (C.S.T.) converge quickly to values only s l i g h t l y i n error of the exact solution. The L.S.T. stresses are not quite as good as the C.S.T. but i n some instances are the same. In general the consistent formulation stresses converge slowly and i n some instances are i n great error with the correct values. F i g . 6 . 6 : GENERAL LAYOUT & LOADING ( 3 X 3 GRID) -103-• TABLE.6.5': DEFLECTIONS:. AND. STRAIN ENERGY: ...... , PARABOLICALLY LOADED SQUARE P L A T E F I N I T E ELEMENT .GRID 10 E t Ug 1 0 2 E t u ( i - v ^ i c 10E V 10 E t V n STRAIN ENERGY U 2 * • l O E t r U ( l - v 2 ) L V 2 IXI - 0 . 9 7 2 6 2 . 4 2 3 5 1 .49079 3 . 9 0 6 1 5 . 2 . 3 3 0 5 5 2 2X2 - 1 . 0 4 9 7 1 1 .57514 1 .32862 4 . 2 1 9 8 6 2 . 3 7 4 3 0 8 3 X 3 - 1 . 1 1 6 2 1 1 .79154 1 .25831 4 . 2 8 9 8 2 2 . 3 9 5 8 0 0 4X4- - 1 .15116 2 . 0 1 2 0 3 1 . 2 1 5 3 3 4 . 3 2 7 6 6 2 . 4 0 7 9 7 2 5X5 - 1 . 1 7 1 5 2 . 1 5 5 5 6 1 .1889 4 . 3 5 0 6 0 2 . 4 1 5 3 2 2 ; 6 X 6 - 1 .1844 2 . 2 4 6 9 7 1 .17198 4 . 3 6 5 5 7 2 . 4 2 0 0 8 3 IOXIO - 1 .2076 2 . 4 0 2 3 1 .14202 4 . 3 9 3 6 6 2 . 4 2 8 8 8 2 EXACT - 1 .519928 1 . 7 8 3 7 1 . 2 7 7 2 7 5 . 0 7 3 4 7 8 2 . 7 9 3 5 6 9 5 * NOTE: U FOR WHOLE P L A T E . PARABOLICALLY LOADED SQUARE PLATE 2.9 + 2.8 2.7 D l i 2 ' 6 N •f Ui o r>2.5 2.4 2.3 2.2 EXACT STRAIN ENERGY = 2.7935 "1 1 h 1 0 % E R R O R H 1 1 i —I 1 I h 20 40 60 80 100 " T O T A L . * O.O.F. 200 320 i STRAIN ENERGY VS TOTAL NO. OF DEGREES, OF FREEDOM Figure 6.7".A 5.0 4.5 0 > UJ Z 4.0 + OU 4,; 3.5 + I EXACT 10'V/ = 5.0734 I07o E i ? R O R 3.0 + 20 40 60 80 100 140 200 320 10 Vj VS TOTAL NO. OF DEGREES OF FREEDOM Fig u r e 6i8 -105' TABLE 6.6: STRESSES PARABOLICALLY LOADED SQUARE PLATE FINITE STRESS C.S.T. L.S.T. G0NSISTEIS1T EACT ELEMENT GRID (K/INI FORMULATION FORMULATION FOPMULATION N » xA 0.51201 0.145689 - 0.14095 N -yA 0.98334 0.52098 0.85904 NxB - 0.11547 - 0.11547 - 0.115468 0.0 yB 0.26352 0.26352 0.263517 0.41067 1 N xc 0.097672 0.19957 0.0 N yc 0.56899 0.295608 0.0 NxD 0.31081 0.31081 0.31081 0.41067 V 0.87447 0.87447 0.87447 1.0 N . xA - 0.047936 - 0.011036 0.12776 - 0.14095 V 0.73788 0.77478 0.56219 0.85904 - 0.01905 - 0.019095 - 9.019095 0.0 yB 0.37108 0.37108 0.371079 0.41067 2 N xc 0.039165 0.12447 0.073391 0.0 N yc 0.42188 0.62819 0.144118 0.0 NXD 0.3809 0.3809 0.380904 0.41067 V 0.96216 0.96216 0.96216 1.0 N » xA - 0.978739 - 0.10797 0.12403 - 0.14095 V 0.78049 0.75126 0.57773 0.85904 - 0.003557 - 0.003557 - 0.003557 0.0 ys 0.38307 0.38307 0.383068 0.41067 3 N xc 0.02494 0.097614 0.04086 0.0 N yc 0.3098 0.4980 0.09228 0.0 N ^ XD 0.39493 0.39493 0.39493 0.41067 N 0.97066 0.97066 0.97066 1.0 -106-TABLE 6.6: STRESSES (CONT'D) PARABOLICALLY LOADED SQUARE PLATE FINITE STRESS C.S.T. L.S.T. CONSISTENT EXACT ELEMENT GRID CK/IN1 FORMULATION FORMULATION FORMULATION NXA N . yA - 0.09141 - 0.14005 0.12215 - 0.14095 0.79883 0.75019 0.58526 0.85904 Nx3 N .„ N XC N yc 0.004855 0.004855 0.004855 0.00 0.39143 0.39143 0.39143 0.41067 4 0.017319 0.24129 0.077801 0.39898 0.02837 0.047519 0.0 0.0 V 0.40544 0.40544 0.405441 0.41067 0.97369 0.97369 0.973689 1.0 xA N -yA - 0.098101 - 0.15529 - 0.12113 - 0.14095 0.80855 0.75135 0.589336 0.85904 xB V 0.0098762 0.009875 0.009876 0.0 0.39638 0.39638 0.396379 0.41067 5 N xc N „ yc 0.012509 0.196605 0.064284 0.32961 0.02199 0.05460 0.0 0.0 NxD N ^ yD 0.41144 0.41144 0.41144 0.41067 0.97471 0.97471 0.97471 1.0 N , xA N -yA - 0.102087 - 0.16399 0.12048 - 0.14095 0.81430 0.75239 0.59174 0.85904 NxB 0.013098 0.01309 0.013097 0.0 0.39940 0.39940 0.39940 0.41067 6 N XC N yc 0.009304 0.165477 0.05470 0.27958 0.01816 0.04574 0.0 0.0 NxD 0.41483 0.41483 0.41483 0.41067 yD 0.97489 0.97489 0.97489 1.0 . - 1 0 7 -TABLE 6.6 CONT'D: STRESSES PARABOLICALLY LOADED SQUARE PLATE FINITE ELEMENT GRID STRESS (K/IN) C.S.T. FORMULATION L.S.T. FORMULATION CONSISTENT FORMULATION EXACT NXA -0.112275 -0.17772 0.11934 -0.14095 NYA 0.823605 0.7545 0.59566 0.85904 IO NXB NYB 0.018819 0.40412 0.018819 0.40411 0.018819 0.40411 0.0 0.41067 NXC 0.003467 0.034398 0.0 NYC 0.100574 0.17194 0.0 NXD 0.41922 0.41922 0.41067 N Y H 0.97364 0.97364 1.0 PARABOLICALLY LOADED SU_UARE~PLATE 0.26 ! 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1-1 2 3 4 5 6 8 10 FINITE ELEMENT GRID SIZE (N) -109-PARABOLICALLY LOADED SQUARE PLATE 0 . 8 6 0 . 8 2 EXACT N = 0 . 8 5 9 0 4 YA C.S.T. 0 . 7 8 L.S.T. 0 . 7 4 0 . 7 0 0 . 6 6 0 . 6 2 f 0 . 5 8 1 CONSISTENT FORMULATION 0 . 5 4 + 0 . 5 0 H 1 1 1 1 1 h —I 1 1 1 —I 1 h 7 8 9 1 0 4 5 6 FINITE ELEMENT GRID SIZE (N) N . VS FINITE ELEMENT GRID SIZE YA FIGURE 6 . JO, -110-PARABOLICALLY LOADED SQUARE PLATE CONSISTENT FORMULATION ^EXACT Nye = O'O 4 5 6 7 FINITE ELEMENT GRID SIZE (N) 10 N y c VS FINITE.ELEMENT..GRID SIZE FIGURE 6.1:1 -111-PARABOLICALLY LOADED SQUARE PLATE .960 I — | 1 1 1 1 1 1 1 1 1 1 1 1 1 r 2 3 4 5 6 7 8 9 10 FINITE ELEMENT GRID SIZE (N) NUT, VS FINITE ELEMENT GRID SIZE FIGURE 6.12. - 1 1 2 -C y l i n d r i c a l Shel l Roof A c y l i n d r i c a l she l l roof is modelled using the f l a t t r i angu la r element derived in Chapter 3- The she l l shown in f i gure 6.13, because of i t s conf igurat ion is between what is normally termed a shallow she l l and what is defined as a deep s h e l l . For th i s ceason, the exact ana l y t i ca l resu l t s obtained by shallow she l l theory and those obtained from deep she l l theory are both presented. The shel l is loaded by i t s own weight, and the loads are lumped at the nodes as v e r t i c a l fo rces . The geometry and a typ ica l gridwork is shown in f i gure 6.13- Only one quarter of the she l l is modelled using symmetry. The de f l ec t i on s are tabulated in table 6.7 and are compared to exact values from reference 3- A l l de f l ec t i ons converge very rap id ly to values only s l i g h t l y in error of the exact, even for r e l a t i v e l y coarse g r id s . The reason for th i s convergence to a value s l i g h t l y in er ror of the exact is because the element is <i:ncomplete and shear s t r a i n constra ints are imposed at the nodes. Figure 6.14 i l l u s t r a t e s g raph ica l l y the va r i a t i on of d e f l e c t i o n w along edge B - c for various gr id refinements. Note the rapid convergence of w to a value s l i g h t l y b o f f the ana l y t i ca l one, even for coarse gridworks. Figure 6.15 p lots wD vs the tota l number of degrees of freedom. Convergence again is rap id , to a value only \% in error of the exact for a 10 x 10 g r i d . Results from a f i f t e e n degree of freedom t r iangu lar element which combjimes the constant s t r a i n t r i ang l e (six degrees of freedom) for the membrane act ion and the Zienkiewicz nine parameter p la te bending element (Ref. 10) is a l so presented. Note the larger error for the f i f t e e n d.o. f . element, even when more d .o. f . are used. For very course -113-gridworks, the f i f t e e n d.o. f . element is fa r too s t i f f compared to the eighteen d.o. f . element used here in. The stress at a node is the average of the surrounding element s t ress cont r ibut ions . The C.S.T. stresses were the same as the L.S.T. s t resses . These values were used instead of the cons is tent formulation stresses (sect ion 4.2.1.) because of the improved accuracy and convergence c h a r a c t e r i s t i c s as was i l l u s t r a t e d in sect ion 6.1. Table 6.8 compares the various membrane and bending stresses (sect ion 4.3) and the s t r a in energy with the exact values (3) obtained from shallow she l l theory. The stresses and s t r a in energy converged to values only s l i g h t l y in error of the exact. However the bending s t ress M appears to be f luc tua t ing cons iderably. The va r i a t i on of N x along edge A - B is p lo t ted for the d i f f e r e n t gridworks in f i gure 6.16. As shown for success ive ly f i ne r g r i d s , the N g^ is rap id ly approaching a value s l i g h t l y in er ror of the exact. In f i gure 6.17 the d i s t r i b u t i o n of along edge D - c is shown for the various g r i d s i ze s . Again i t is seen that even for the extremely coarse gr id the error is smal l . In general the de f lec t ions and stresses compare except iona l ly well with the exact values but appear to be converging to values s l i g h t l y in error of the pred ic ted. The reason as was mentioned e a r l i e r is due to the fact that the element is incomplete and shear s t r a in constra ints are imposed at the elements' nodes. NOTE: MODEL 1/4 OF STRUCTURE DUE TO SYMMETRY CONDITIONS. F i g . 6.13: CYLINDRICAL SHELL GEOMETRY ( 3 x 3 G R I D ) -115-TABLE 6.7 : Iff LECTIONS CYLINDRICAL SHELL FINITE ELEMENT GRID NET NO. CF EQUATIONS iou A ( IN) A WB (IN) V B (IN) 10 W . (,IN) ' 2 x 2 30 -0.735 -4.571 2.375 6.010 3 x 3 63 -1.049 -3.629 1.912 5.281 4 x 4 108 -1.201 -3.530 1.861 5.234 5 x 5 165 -1.285 -3.527 1.860 5.275 10 x 10 630 -1.417 -3.564 1.881 5.414 EXACT * ? -3.607 ? EXACT * * -1.5133 -3.7033 1.9637 5.2494 NOTE : * FROM DEEP SHELL THEORy * * FROM SHALLOW SHELL THEORY +2.0 W - DEFLECTION ALONG EDGE B-C FIGURE 6.14-C Y L I N D R I C A L S H E L L ROOF - 3.G + \A/s E X A C T = 3.601 I N . I ERROR 3.4 + 3.Z + 3.0 + 2 8 N O T E : * F R O M D E E P S H E L L THEORY F L A T 15 D.O.F. E L E M E N T ( C S I + 9 P L A T E . &EKIOIK16) R E F tO. O FLAT 18 D.O.F. E L E M E N T ( PRESENT) I 2 .6 2.4-+ 4- 1* 200 -4o0 GOO 80O iooo I ZOO I400 T O T A L NO. O F O E G R E E S O F F R E E D O M oo Wa vs T O T A L N O . O F D E G R E E S O P F R E E D O M F I G U R E 6 . 1 5 -118-CVUNDRICAL SHELL ROOF TABLE 6.3•: STRESSES . (C.S.T.) GRID SIZE NXB (K/IN) MXC (K-IN/IN) MYC (K-IN/IN) STRAIN ENERGY (K-IN) # 2 x 2 3.065 -0.099 -2.010 70.93 3 x 3 4.536 -0.075 - 1.822 56.62 4 x 4 5.223 0.005 - 1.788 55.56 5 x 5 5.591 0.085 - 1.744 55.79 10 x 10 6.113 0.254 -1.650 56.78 EXACT * ? ? ? EXACT * * 6.4124 0.0927 -2.056 58.828 NOTE : # STRAIN ENERGY FOR TOTAL STRUCTURE * FROM DEEP SHELL THEORY * * FROM SHALLOW SHELL THEORY K^B-KFREE EDGE AT LINE OF SYMMETRY X DISTANCE FROM PT. A (INCHES) Nx ALONG EDGE A-B FIGURE 6.16 My - ALONG EDGE D-C FIGURE 6.17 -121-6.5 Point. Loaded.Spherical .Shell ... A spherical, shell. Is modelled by. using, the. f l a t eighteen degree of . freedom triangular, f i n i t e element derived.in Chapter 3. The spheri-c a l s h e l l subjected to a point load creates regions of large bending stresses, a region where there are mainly membrane stresses and a region of high stress concentration. Because the problem i s axisym-metric and rectangular co-ordinates are used throughout the analysis, one quarter of half the sphere i s modelled. A non-uniform grid spacing i s used where the ratios of successive elements are taken as 1 : 2 : 3 : 4, .. .N, to provide a better representation of results i n the region of high stress gradients near the pole. Figure 6.18 illustrates the general layout and the loading. The deflections resulting from tteepoint loading at the pole are tabulated i n table 6 . 3 . The exact values obtained analytically reference 6 are also given. Convergence i s rapid and for a relative-l y coarse grid, the deflections at the pole and equator compare extremely well. Figure 6.19 illustrates graphically the rapid convergence to a value sl i g h t l y i n error of the exact of the z di s -placement at the pole with the number of elements used. The normal displacement near the pole vs the colatitude direction along the sphere i s shown i n figure 6 . 20 with the exact values. The displace-ments compare very well with the exact ones. The stresses based on the C.S.T. and L.S.T. formulation are presented in table 6.IH with the exact ones. The L.S.T. stresses appear to be better where they d i f f e r radically from the C.S.T. values. In many -122-instances the C.S.T. and L.S.T. stresses..are aliiost;.the same. In general the stresses, are only: slightly in error, .of. the exact values. Table 6.) 0 compares the stresses -N. . N_, M e and M at the pole and equator.with the exact.ones based on shallow shell theory. These stresses are frcm the L.S.T. formulation and are reasonably close to the exact values for the finer grids. The distribution of the membrane stresses near the pole in the oolatitude direction are shown in figure 6.21:. Again both N© and are close to the exact solution (FLUGGE - reference 6) . More remote from the pole at oola-titude angles of twenty to ninety degrees, the membrane stresses are compared to the exact ones in figure 6.22L. The stresses based on the finite element solution follow very closely the distribution of the stresses based on the analytical results. In general the displacements away from the pole region compare very closely with the exact solution. Away from the pole the bending stresses die out and the membrane stresses dominate. These membrane stresses remote from the pole follow the analytical values very closely. Again i t is seen (figure 6.I9») as in the previous examples that the displacements appear to be converging rapidly but to values only slightly in error of the exact. This is due to the fact that the element is incomplete and also because of the nodal shear constraints introduced. - 1/4 OF SPHERE IS MODELLED USING A NONUNIFORM GRID FOR COLATITUDE DIRECTION. - RATIO OF SIDES IS 1:2:3...N N = No. OF ELEMENTS HIGH. - N = 4 IS SHOWN. X FIG. 6.15: SPHERICAL SHELL - 1 2U-TABLE 6.3; DEFLECTIONS POINT, LOADED SPHERE,... F INITE DEFLECTION DEFLECTION I LEMENT RID @ POLE CIN..1 (Et W/P) @ EQUATOR (3N.1 CEt W/PI 2 6.1054 - 0.2908 NOTE: 4 8.0346 - 0.2235 * Value for deep shell 8 20.1638 - 0.1993 Theory. 10 21.8660 - 0.19901 ** Value for Shallow Shell 12 22.3918 - 0.1984 Theory. 14 22.4478 - 0.1979 * EXACT 21.200 ? STRESS UNITS: ** EXACT i 21.093 - 0.2069 N = K. / IN-ip M = K. - IN/ TN. TABLE 6 .10 : STRESSES (L.S.T.) FINITE AT POLE AT EQUATOR ELEMENT GRID. N e \ M Qx 10 - 1 M^x 10 1 N e M„ x 10~ 3 % x 1 0 " 3 2 - 0.9589 2.0684 - 0.0788 - 0.069 0.4996 0.3043 0.5376 0.2162 4. 4.7133 5.099 0.501 0.287 -0.1861 0.1329 - 6.319 - 5.987 8 12.388 12.331 1.307 - 0.601 -0.1334 0.1473 0.0267 0.00310 10 12.660 12.637 2.210 -0.913 -0.1318 0.1462 0.0279 0.00347 12 12.465 12.455 3.012 - 1.209 -0.1301 0.1458 0.0264 0.00217 14 12.180 12.177 3.685 - 1.470 -0.1287 0.1458 0.0278 0.00134 E ** XACT 10.313 10.313 oo oo -0.1592 + 0.1592 = 0 * 0 i -125-POINT LOADED SPHERE 24.0 + 22.0 m 20.0 « S5 w o 18.0 Pn fn O H W § 16.0 Pn! CO M Q N 14.0 12.0 + 10.0 t 8.0 5 2 ERROR EXACT W = 21.2 * NOTE: * FROM DEEP SHELL THEORY 1 1 1 1 h 6 8 10 12 14 16 NUMBER OF ELEMENTS IN BOTH DIRECTIONS (N) DEFLECTION AT POLE VS FINITE ELEMENT GRID FIGURE 6.13-. •126-<J! DEGREES (EROM PEAK) NORMAL DISPLACEMENT VS ANGLE Q NEAR POLE FIGURE 6.2Q -127-TABLE 6.11 STRESSES POINT LOADED SPHERE FINITE ELEMENT GRID STRESSES (K/IN) C.S.T. FORMULATION L.S.T. FORMULATION EXACT 2 N0P NQP QE -0.5974 2.1221 -0.6921 0.0746 -0.9589 2.0684 0.4996 0.3043 10.313 10.313 -0.1592 0.1592 4 N 6 P N QE 3.3125 4.7298 -0.18113 0.15643 4.7133 5.099 -0.1861 0.1329 10.313 10.313 -0.1592 0.1592 8 N0P N INQP QE 10.546 11.7947 -0.1603 0.1498 12.388 12.331 -0.1334 0.1473 10.313 10.313 -0.1592 0.1592 10 N e P N INQP N 9 E QE 11.5366 12.305 -0.1564 0.14658 12.660 12.637 -0.1318 0.1462 10.313 10.313 -0.1592 0.1592 12 N0P N INQP QE 11.7286 12.2365 0.1542 0.1457 12.465 12.455 0.1301 0.1458 10.313 10.313 -0.1592 0.1592 NOTE: - SUBSCRIPT P ^ P O L E - SUBSCRIPT Ecp EQUATOR MEMBRANE STRESSES VS ANGLE ($ NEAR POLE FIGURE 6.Z\9 MEMBRANE STRESSES VS (ft ANGLE REMOTE FROM POLE FIGURE 6.22 -130-6.6 Non-pri smatic. .Folded-Plate Structure The next application i s to a non-prismatic folded plate structure. This structure w i l l deform i n an ^symmetrical manner when subjected to loads. The structure's geometry and loading are shown i n figure 6.23. A uniform line loading i s applied at +ojp -Folcl lines as indicated. The basic plate units which make up the structure are trapezoidal i n shape (shown in figure 6.24). The eighteen degree of freedom f i n i t e element developed i n Chapter -3 i s used. The various gridworks employed are shown i n figure 6.25. The results are compared with: (1) Experimental (reference 5) - Tests performed on a scale model. (2) Analytical (reference 5) - A theory for long non-prismatic folded plates i s presented and applied. (3) High Order Finite Element (Beavers - reference 1) - A f i n i t e element representation using a high order f i n i t e element i s presented. A complete quintic polynomials -isa used for bending and complete cubics are u t i l i z e d for the two in-plane displace-ments. An eighteen degree of freedom in-plane element i s combined with an eighteen degree of freedom plate bending element, resulting i n a thirty-six degree of freedom triangular element. Sprecial constraint equations are also introduced for the skewed boundaries. The stresses are computed by averaging the stresses of the surrounding element contributions that are a l l coplanar to each other. The stresses presented include the membrane stresses which are constant „ •' -131-over the thickness of. the element,.and the small bending stresses that. are. assumed.to be. constant-across the widttuof, the element but are extremely small. . The stresses, from.the L.S.T. formulation Csection .4.2.31'. are presented for ..the. membrane portion and the bending .stresses based on section 4.3 are presented. Table 6.represents the ve r t i c a l deflection along.fold lines c and E for the various f i n i t e element grids. The deflections obtained from the element derived i n Chapter 3 herein, compare reasonably well with the Beavers Cl), experimental (5) and analytical (5) results. This i s shown graphically on figure 6.26- the vertical deflection along fold line c i s plotted for each of the gridworks used. Notice the steady convergence toward the analytical result for each subse-quent grid refinement. The longitudinal stresses along fold line c and E are tabulated i n table 6.13 for each gridwork. Again the values appear to be steadily converging to the experimental and analytical results wAh each grid refinement. The stresses presented are based on the C.S.T. and L.S.T. stresses (section 4.2.2 and 4.2.3) . Figures 6.27 and 6.28 i l l u s t r a t e graphically the variation of longitudinal stresses along fold lines C and E respectively for the different gridworks. In both figures one can see the rapid convergence toward the analytical values. The transverse moment at midspan i s ill u s t r a t e d i n figure 6.29 for the gridwork N = 128. The moments are compared with Beaver's high order f i n i t e element results. The results are not too far apart -132-and i n some ..instances ,differ_ only- slightly. In general the. f i n i t e ..element... representation, of ..the non-prismatic folded plate structure yielded reasonably, good results. i -133-X NOTE; ALL ANGLES y = 40° LOAD = 2.334 #/IN. OF HORIZONTAL PROJECTION ALUMINUM MATERIAL E = 10.4 x 103 KSL t = 0.063 ir = 0.33 8 PLATES IN ALL NONPRISMATIC FOLDED PLATE FIGURE 6.23 - 1 3 U -NONPRISMATIC FOLDED PLATE ACTUAL PLATE DIMENSIONS EXAGGERATED PLATE GEOMETRY PLATE GEOMETRY FIGURE 6.24 NONPRISMATIC FOLDED PLATE -135-MODEL PLAN VIEW 32 ELEMENT MESH 64 ELEMENT MESH 128 ELEMENT MESH FIGURE 6.25: MODEL & MESH PATTERNS TABLE 6.J2: DEFLECTIONS NONPRISMATIC FOLDED PLATE Fir ELE ITE MENT DISTANCE ALONG FOLD VERTICAL DEFLECTION ALONG FOLD LINE C X 1 Q 3 IN VERTICAL DEFLECTION* ALONG FOLD LINE E X l n - 3 GRID LINE OF L F.E. BEAVERS EXPT. ANALYTIC F.E. EXPT. ANALYTIC 1/4 L 1.402 2.2 2.3 2.1 2.027 3.1 3.1 32 I 1/2 L 2,168 3.5 3.1 3.2 2.502 3.9 3.8 3/4 L 1.800 2.9 2.6 2.7 1.652 2.4 2.5 1/4 L 1.532 2.2 2.3 2.1 2.264 3.1 3.1 64 1/2 L 2.384 3.5 3.1 3.2 2.755 3.9 3.8 3/4 L 2.003 2.9 2.6 2.7 1.814 2.4 2.5 I 1/4 L 1.816 2.2 2.3 2.1 2.693 3.1 3.1 128 1 I 1/2 L 2.832 3.5 3.1 3.2 3.304 3.9 3.8 i i 3/4 L 2.365 2.9 2.6 2.7 2.165 2.4 2.5 i I I NOTE: F.E. = RESULTS FROM FINITE ELEMENT DERIVED IN CHAPTER 3 i NCNPRISMATIC FOLD PLATE NOTE: o For N= 128 a For N = 64 ^ • For N = 32 A For EXPTAL — For ANALYTIC FIGURE 6.26: VERTICAL DEFLECTION ALONG FOLD LINE C -138-TABLE 6.I3 LONGITUDINAL STRESSES (PSI) NONPRISMATIC FOLDED PLATE j FINITE ELEMENT GRID DISTANCE ALONG FOLD LINE C DISTANCE ALONG FOLD LINE E (C.L.) 1/4 L 1/2 L 3/4 L L 1/4 L 1/2 L 3/4 L L 32 -94.14 •* -178.90 -271.69 -386.66 -190.95 -324.32 -237.79 -350.47 -196.44 -334.67 -245.49 -396.11 -87.00 -201.09 -85.60 -162.48 64 -94.95 -203.70 -303.97 -438.39 -259.71 -415.70 -187.64 -270.51 -264.30 -427.83 -271.63 -455.02 -84.73 -240.86 -41.48 -135.80 12 8 -226.47 -279.25 -466.69 -533.34 -465.92 -546.30 -245.90 -294.44 -488.65 -574.35 -488.20 -575.68 -260.57 -331.90 -77.55 -99.71 EXP ANA TAL. LYTIC -356. -344. -586 -573 -666. -635. -710. -676. -702 -655 -418 -428 NONPRISMATIC FOLDED PLATE °FOR N = 64 o FOR N = 128 A FOR EXPTAL FOR ANALYTIC FIGURE 6.27: LONGITUDINAL STRESS ALONG FOLD LINE C (L.S.T. STRESSES) NPM -PRISMATIC-. FOLDED PLAT E L 1 ? NOTE: • FOR N = 32 a FOR N = 64 o FOR N = 128 A FOR EXPTAL FOR ANALYTIC FIGURE 6.28: IOSfGIIUDINAL STRESS ALONG FOLD LINE E CL. (L.S.T. STRESSES) NONPRISMATIC FOLDED PLATE — A REPRESENT BEAVER'S RESULT FOR N = 32 UNITS FT-LBS/FT FIGURE 6.29: TRANSVERSE MOMENT AT MIDSPAN For N = 128 -11*2-Beam .Stiff ener ..Application; The twelve degree- of., freedcm beam., s t i f f ener. element developed i n section 5.1 i s tested using two., different, load cases. The section in each case i s syirtmetric. The general layout i s shown i n figure 6.30. The beam elements support a thin flexurally weak plate which i s modelled with the f i n i t e element developed i n Chapter 3. For load case one, the beam i s simply supported and a ver t i c a l load i s applied at midspan. From flexural theory, the maximum deflection i s computed as a check. The beam elements yielded an answer less than two per cent i n error. Load case two i s a moment applied at 30 degrees to the major principal plane of the section (refer to figure 6 31?) at each end of the simply supported beam. The deflection was again computed from ref-erence 8 as a check. The result using the beam elements was less than one per cent i n error. I t appears that the stiffness matrix derived i n section 5.1 for the beam element, using the strain energy approach i s an accurate representation. -11*3-S Sxio E = 30,000 K S I FIGURE S30'': BEAM STIFFENER PROBLEM -I--BEAM STIFFENER PROBLEM LOAD CASE (1) VERTICAL LOAD APPLIED AT MIDSPAN. RESULTS: P = -l.OK FROM FLEXURAL THEORY A C L . = PL_ = -0.021" MAX PROGRAM YIELDED % ERROR IN A = 1.9% LOAD CASE ( 2 ) : 4 8 E I — A C L . = -0.0206" MAX. Please refer to figure 6.31" - M i (—' V 8118.4 F i g u r e . 6.31 M = 50"K ONE SUPPORT: s.s . = 43.301"K M z = - 2 5 . 0 0 " K OTHER SUPPORT S . S My = - 4 3 . 3 0 1 " K M z = 2 5 . 0 0 " K A = 5 . 3 4 , I y = 5 6 . 9 , I z = 3 . 8 J = 6 0 . 7 , L = 1 4 4 " R z = 0 . 8 4 4 Ry = 3 .264 RESULTS: PREDICTED A MAX = 0 . 5 7 2 2 ' FROM TLMOSHENKO STR. OF MAT'LS. PG. 232 A MAX 2 = (ML 2 cose ( L 2 \ + (ML2 s iNe \ k 8 E i Y ; ^ 8 E I Z ; PROGRAM RESULTS: A MAX •Y » 2 , 2 . 2 X Y + A Z 1 = 0 . 5 6 7 7 " % ERROR IN A = 0.786% - 1 U 5 -CCNGLUSIONS Presented herein has been a shallow shell element of ar b i -trary triangular shape. The element was; developed by caribining a nine degree of freedom plate bending element with a nine degree of freedom in-plane element. An incomplete cubic polynomial was used to describe the normal out of plane displacement and cubic polynomials were used to describe the two in-plane displacements. Constraints and st a t i c condensation were used to reduce the number of generalized co-ordinates for the in-plane displace-ments. The eighteen degree of freedom triangular f i n i t e element was developed with the intent of modelling plate and shell structures. I t i s assumed that the behavior of a continuously curved surface can be adequately represented by the behavior of a surface b u i l t up of small f l a t elements. The stresses are computed three different ways. The consistent formulation Cstxain-displacement matrix, etc.) i s compared with the constant strain triangle stresses. A technique i s developed to compute the midside node displacements from the vertex nodes and the element configuration. Then the linear strain triangle stresses are computed and compared to the other two stress results. To assess the new nine parameter plane stress element, a parabolically loaded square plate was modelled. The plate, due to i t s in-plane loading, had only membrane stresses. The deflections and strain energy converged rapidly to values only marginally i n error of the exact solution. The consistent stresses were very poor but the constant strain -1 US-triangle and linear strain.triangle-stresses-converged xapidly to values that compared closely, with, the exact values:. A cylindrical shell roof was represented next. Loaded only by i t s own weight, the load was lumped at the various nodes as v e r t i c a l forces. In general the deflections and stresses (C.S.T. and L.S.T.) con-verged rapidly to values only s l i g h t l y off the analytical results. Even for relatively coarse grids, the results obtained were reasonable. Vfe wanted to investigate how the element might perform i n regions of large bending stresses, regions of large membrane stresses and f i n a l l y i n regions of high stress concentration. So a point loaded spherical shell was modelled. The results again indicated relatively rapid conver-gence and reasonable accuracy with the analytical values for both deflections and stresses. In each case the deflections, stresses and strain energy appeared to converge f a i r l y rapidly toward values sli g h t l y i n error of the analytically predicted ones. This characteristic i s attributed to the fact that shear strain constraints were used at the nodes and the f i n i t e element i s incomplete. A non-prismatic folded plate structure was studied next. We were not sure how the element would act for this type of unsymmetrical bending and whether the fold lines might introduce errors. However, the results were quite encouraging. The deflections and stresses were compared to experimental, analytical and a f i n i t e element analysis using a high order element. -11*7-A twelve degree of freedom beam stiffener element was formulated using the .strain.energy, expression, with, the intent of combining i t with, the f i n i t e element.. At f i r s t the formulation was performed for a symmetric crossecfdon. Then two numerical, examples were tested. The def-lections were only marginally i n error with those predicted, from flexural theory even when the beam stiffeners were loaded unsymmetrically. Later the formulation was generalized to include beam stiffeners with unsymmetri-ca l crossections. -1 US -APPENDIX A . l , DKCUSSICN OF PROGRAM A computer program using Fortran IV language was developed for the analysis of folded plate and shell structures. The program u t i l i z e s the eighteen degree of freedom f i n i t e element and the twelve degree of free-dom beam s t i f f ener element based on the theory discussed ea r l i e r . A general flow chart of the program i s given i n Appendix A.3.- " Given a structure, a geometrical model i s constructed from i t . The model i s divided up into a suitable gridwork of elements. These triangular elements should have relatively low aspect ratios although i t i s not essential. Next the apexes of these elements are numbered but care should be taken so as to rninimize the band width of the master- stiffness matrix. With the nodal points numbered, the degrees of freedom are determined next by surmiing the constraint numbers. For each node i t must be determined i f some nodal movements are inhibited from motion or not. This vector of nodal move-ments (constraints) represents the boundary and syrnrnetry conditions of the structure. The appropriate node numbers are then associated with each element; The beam stiffeners are treated the same way. Note that each beam stiffener element only extends over the region of one f i n i t e element. This way the band width from the f i n i t e elements i s not destroyed. -11+9-The.main features of the program can be .considered to be divided into the following procedure: 1) Number nodal degrees of freedom, establish, band width, check problem size, and read i n Finit e Element data. 2) If beam stiffeners are used read i n the pertinent data. 3) Compute the bending element stiffness matrix. 4) Compute the in-plane element stiffness matrix. 5) Combine the bending and in-plane stiffness matrices and build the structure (master) stiffness matrix. 6) If beam stiffeners are used compute each beam stiffener's stiffness matrix and add to the structure stiffness matrix. 7) Build the master load vector. 8) Solve for the unknown degrees of freedom '(nodal displacements). 9) Compute the membrane stresses and bending stresses for each element then find the resultant values at each node by averaging a l l surrounding element contributions. Of course co-ordinate transformations and other steps have been omitted but these represent the core to the whole procedure. The program i s set up to handle 2,000,000 bytes. One million of a r e these v . set aside for the master stiffness matrix. This means that the -15"0-Master stiffness matrix can handle 125,000 double precision, words Cor two f u l l words.!.., The othea: 1,.000,000 bytes...are used by .the remainder of the program. .. The.examples presented herein, did. not. u t i l i z e a l l of. this availa-ble core area. Note: A l l units are expressed In kips and inches. A l l real numbers are double precision and a l l integers are single f u l l words. - 1 5 1 -. - APPENDIX. A. 2 INPUT DATA A description of Input items:Is discussed, following Table A.2.1. TABLE A.2.1: FORMAT OF INPUT DATA CARDS IDENTIFIER DESCRIPTION FORTRAN FORMAT CARD COLUMNS 1 NLC NSTRT NDOF ir T E NBEAM NOELEM ITER 3 NE NNODES NVAR TOTAL NO. OF LOAD CASES 15 STRUCTURE IDENTIFICATION NO. 15 CONTROL FOR DUPLICATING DEGREE 15 OF FREEDOM NO. ( NO. = NUMBER) POISSON'S RATIO FOR F.E. CFINITE ELEMENT) THICKNESS OF F.E. YOUNGE'S MODULUS OF ELASTICITY FOR F.E. TOTAL NO. OF BEAM STTFFENERS USED (XNTROL FOR WHETHER PROBLEM IS TO BE SOLVED WITH OR WITHOUT F.E. NUMBER OF ITERATIONS REQUESTED FOR VARIABLE BANDWIDTH MATRIX DECOMPOSITION ROUTINE TOTAL NO. OF FINITE ELEMENTS IN PROBLEM 15 TOTAL NO. OF NODES IN PROBLEM 15 NO. OF VARIABLES (DEGREE OF FREEDOM) PER NODE 15 F5.3 F5.3 F15.2 15 15 15 I - 5 6- 10 II - 15 I - 5 6- 10 II - 25 26-30 31-35 36-40 I - 5 6-10 I I - 15 - 1 5 2 -...TABLE-A. 2.1 (CONT'D) IDENTIFIER.' rjESCKLPTION FORTRAN CARD FORMAT COLUMNS NNODEL NO. OF NODES PER ELEMENT ICOFW FIELD WIDTH: USED FOR READING IN EACH ELEMENT'S NODE NUMBERS. (EXPLAINED FOLLOWING TABLE) NODAL DATA (FOR EACH NODE) X, Y, AND Z aX3FC)INATES AND NODAL CONSTRAINTS (IX VECTOR) : IF NDOF = OR TF NDOF = OR IF NDOF = 0 1 2 FINITE ELEMENT DATA (FOR EACH ELEMENT) ICO (I, J) , J NODE NO.'S FOR THE I'TH ELEMENT : IF ICOFW = 2 OR IF ICOFW = 3 BEAM STTFFENER DATA (FOR EACH STIFFENER) JNL (LOWER NODE NO.) JNG (GREATER NODE NO.) JNP (ORIENTATION NODE USED TO DEFINE ORIENTATION OF STIFFENER'S WEAK PLANE) X (JNP) Y (JNP) Z (JNP) ) J GLOBAL (XORDINATES OF JNP NOT DEFINED IF JNP IS NOT SPECIFIED - EXPLAINED POUJOWING THIS TABLE. 15 15 3F10.0 612 614 615 312 313 15 15 15 F10.1 F10.1 F10.1 16-20 21.25 1 -30 31-42 31-54 31-60 1 - 6 1 - 9 I - 5 6- 11 I I - 15 I- 10 I I - 20 21-30 -153-TABLE A.2.1 CCONT'D] IDENTIFIER DESCRIPTION FORTRAN CARD FORMAT COLUMNS FOR EACH BEAM STIFFENER ** RG RS A ECG ECS PJ E IZ IYZ RY RADIUS OF GYRATION CGREATERj RZ RADIUS OF GYRATION (SMAUEFt) TOTAL XSECT. AREA OF STIFFENER EY ECCFJSITRICITY (GREATER) EZ EOCENTRICITY C S M A U » E R ) POLAR MOMENT OF INERTIA YOUNGE'S MODULUS OF ELASTICITY FOR BEAM STIFFENER MATERIAL SHEAR MODULUS OF ELASTICITY MOMENT OF INERTIA W.R.T. Z AXIS PRODUCT OF MOMENT OF INERTIA. W.R.T. Y AND Z AXES. F7.3 F7.3 F7.3 F7.3 F7.3 F7.3 F7.3 F7.3 F7.3 F7.3 1 - 7 8- 14 15-21 22-28 29-35 36-42 43-49 50-56 57-63 64-70 Note: * I f the beam stiffener i s symmetric then the value of IZ can be any value other than zero, but i t must be entered. (Stiffener w i l l bend with only R ) y ** I f a l l beam stiffeners are the same shape, enter a 0.0 for RG on sub-sequent cards and the values on the previous card are assumed. ** If a l l beam s t i f feners are of the same material, enter a 0.0 for E on subsequent cards and the values on the previous card are assumed. Refer to Figure A.1.1. -15V TABLE A.2.1 CGOSIT1D) IDENTIFIER DESCRIPTION FORTRAN CARD FORMAT COLUMNS LOAD INFORMATION JNODES TOTAL NO. OF LOADED NODES IVERT CONTROL USED TO INDICATE IF LOADS NEED TO BE TRANSFORMED TO THE GLOBAL SYSTEM. 15 15 1 - 5 6- 10 FOR EACH LOADED NODE (K AND INCHES} KNODE LOADED NODE NO. FX LOAD APPLIED IN X-DIRECTION FY LOAD APPLIED IN Y-DIRECTION FZ LOAD APPLIED IN Z-DIRECTION MX MOMENT APPLIED ABOUT X-AXIS MY MOMENT APPLIED ABOUT Y-AXIS MZ MOMENT APPLIED ABOUT Z-AXIS IEL IF LOADS ARE TO BE TRANSFORMED TO THE GLOBAL SYSTEM, THIS IS THE ELEMENT WHICH IS NORMAL TO FZ AND PARALLEL TO FX AND FY. (EXPLAINED FOLLOWING THIS TABLE) 15 F10.2 F10.2 F10.2 F10.2 F10.2 F10.2 15 1 - 5 6- 15 16-25 26-35 36-45 46-55 56-65 1 - 5 Detailed Description: The first card of the program allows the user to assign the structure an identification number so that he may easily refer to i t at some future date. -155-If more than one load case.is to be applied to the structure,, then the solu-tion routine saves, the decomposed ..structure , s t ^ the subse-quent displacements and. stresses, are computed..very, quickly, without having to.decompose,the .structure stiffness matrix each time.. The NDOF i s used to f a c i l i t a t e where one wants to ..assign, duplicate, degree of freedom numbers to various nodes. Here i s how i t Is used: - I f no duplicate degree of freedom numbering i s desired, leave NDOF blank. - If you wish to use duplicate degree of freedom numbering, then - for reading i n actual degree of freedom number i n fields of 4, enter 1 for NDOF - for reading i n actual degree of freedom number i n fields of 5, enter 2 for NDOF Example Want node 13's degree of freedom to be same as node 4's degree of freedom, then enter - 4 - 4 - 4 - 4 - 4 - 4 for constraints of node 13. Example Want w of node 16 to be same as U of node 5, then enter 0 1 19 0 1 1 for constraints of node 16,, where +he ac+uoL d o f no. '3 C o + h e ac+uoL <A.o:f- no. -for u. duspl. o"f n o d e 5. The second card defines the material properties of the f i n i t e elements. A l l fi n i t e elements are assumed to have the same thickness T. If beam stiffeners are used, then enter the total number ,'(NBEAM) . I f no beam s t i f f eners are used, then leave NBEAM blank. Note: Each beam s t i f f ener element extends over the length of one f i n i t e element only. I t may be desired to run a - 1 5 6 -s t r u c t u r e t h a t i s composed., o f beam., s t i f feners. only:, o r you may. wish, t o n e g l e c t the strength, o f t h e .adjoining. s l a b of. f i n i t e . elements. ... If.. the. user e n t e r s a 1 f o r NOELEM, t h e n only, the stiffness;.of..the,beam.stxffeners. i s considered and the d e f l e c t i o n s due. t o the a p p l i e d l o a d s are. computed. Normally one would wish.to i n c l u d e the e f f e c t . o f t h e . p l a t e of. f i n i t e elements so NOELEM i s l e f t blank.. ITER i s the number, o f . i t e r a t i o n s used by the s o l u t i o n r o u t i n e f o r decomposing the v a r i a b l e bandwidth, master s t i f f n e s s m a t r i x . F o r no i t e r a t i o n s , t h i s v a l u e i s l e f t blank. The t h i r d c a r d i s used t o i n d i c a t e the t o t a l number o f nodes and elements i n a problem. F or the element used, NV7AR, the number o f degree o f freedom per node i s s i x . The number o f nodes per element, NNODEL i s t h r e e . The v a r i a b l e ICOFW i s used t o i n d i c a t e the w i d t h o f the f i e l d s f o r r e a d i n g the node numbers o f each element. - I f the t o t a l number o f nodes i s l e s s than o r equal t o 99 then s e t ICOFW = 2 and the nodes are read i n f i e l d s o f 2. - I f the t o t a l number o f nodes i s g r e a t e r than 99 then s e t ICOFW =3 and the nodes w i l l be read i n f i e l d s o f 3. The f o u r t h i tem regards s p e c i f y i n g the g l o b a l x, y, and z co-ordinates and the s i x c o n s t r a i n t values f o r each node. The s i x c o n s t r a i n t s correspond t o u, V, w, e , 0 , 6 movements. E i t h e r a 1 (free) o r a 0 (fixed) i s entered x y z f o r each o f the c o n s t r a i n t s . Note: a l l l ' s do not heed t o be entered s i n c e a b l a n k here represents a 1 Cfree ele:, unconstrained motion). The f i f t h , item e n t a i l s denoting the th r e e node numbers which correspond t o each element. These values are entered t h r e e p e r ca r d (each element) and -157-are i n f i e l d widths according to the value of ICOFW. If beam s t i f f eners are not. used i n a particular problem, then, items six and seven can be disregarded. Item six regards inputting the lower node number (JNLl, the greater node number (OTSfG) and a third node number (JNPl of the beam s t i f f ener. The JNP node number's co-ordinates are used to define the orientation of the weak plane of the s t i f f ener. There are three cases which could exist: (1) The weak axis of the s t i f f ener i s i n the x - y plane (horizontal) i . e . the stiffener i s v e r t i c a l . Then JNP does not need to be entered. The X (JNP), Y (JNP) and Z (JNP) does not need to be entered either. (2) The weak plane of the stiffener i s not i n the horizontal plane but i t s orientation can be defined by using the co-ordinates of a known node. Then the node number i s entered for JNP. On the following card enter only - 0.0 for X (JNP) and leave Y (JNP) and Z (JNP) blank. ( 3 1 The weak plane of the stiffener i s not i n the horizontal plane and i t s orientation has to be described by introducing the co-ordinates of a new constrained node. Give JNP a number i n the range [ NNODES + 200 < NNODES + 400] and on the next card, enter the values of X (JNP), Y (JNP) , Z (JNP). Item six i s done for each beam stiffener. Item seven i s also done for each beam stiffener. On each card, one per stiffener, the section and material properties are entered (noted i n Table A.2.1) . -15*8-I t e m e i g h t t h e t o t a l number o f l o a d e d n o d e s . QNODESl i s . , e n t e r e d . I f t h e l o a d s a r e a c t i n g i n t h e Z - d i r e c t i o n . ( v e r t i c a l ) , , t h e n t h e r e i s n o n e e d t o t r a n s f o r m t h e m , s o a b l a n k , o r 0 i s e n t e r e d f o r I V E R T . I f t h e y a r e a c t i n g i n a d i f f e r e n t d i r e c t i o n , t h e n t h e y s h o u l d b e t r a n s f o r m e d . t o t h e g l o b a l s y s t e m b e f o r e t h e . m a s t e r l o a d v e c t o r i s b u i l t , s o a 1 i s e n t e r e d f o r I V E R T . I t e m n i n e ; f o r e a c h l o a d e d n o d e , i t s number i s e n t e r e d a n d t h e n i t s m a g n i -t u d e ( F ,. F , V M , M , M 1. I f t h e l o a d h a s t o b e t r a n s f o r m e d ( IVERT = 1) , x y x y z t h e n t h e n e x t c a r d s h o u l d i n d i c a t e t h e number o f t h e e l e m e n t ( IEL ) f o r w h i c h F z i s n o r m a l a n d F x a n d F ^ a r e a c t i n g i n t h e same p l a n e ( w . r . t . l o c a l a x e s o f t h e e l e m e n t ) . I f IVERT = 0 t h e n IEL i s n o t e n t e r e d . N o t e : I t i s i m p o s s i b l e t o l o a d i n a d i r e c t i o n w h i c h i s c o n s t r a i n e d f r o m m o t i o n . - 1 5 9 -The x-axis runs of the stiffener runs along i t s length ONE BEAM ELEMENT FIGURE A. 1.1 BEAM STIFFENER SECTION PROPERTIES -160-APPENDIX A .3 FLOW CHART FOR COMPUTER PROGRAM START READ & WRITE -STRUCTURE I.D. NO., NO. OF LOAD CASES -F.E.- MATERIAL PROPERTIES, NO. OF NODES, ELEMENTS & BEAM STIFFENERS. READ & WRITE - ALL FINITE ELEMENT (F.E.) DATA, i.e. NODAL CONSTRAINTS & GLOBAL CO-ORDINATES, NODE NO. FOR EACH ELEMENT. - COMPUTE THE ACTUAL D.O.F. NO. FOR EACH ELEMENT & THE HALF BAND WIDTH FOR MASTER STIFFNESS MATRIX (K). YES -161-READ & WRITE - ALL BEAM STIFFENER DATA: i . e . MATERIAL; & SECTION PROPERTIES & STIFFENERS' ORIENTATION FOR EACH FINITE ELEMENT COMPUTE THE 18X18 ELEMENT STIFF. MATRIX [R] IN GLOBAL CO-ORD. TRANSFORM THE NODAL GLOBAL CO-ORD. TO THE LOCAL SYSTEM COMPUTE THE 9X9 PLATE BENDING STIFFNESS MATRIX (je] BUILD CK] BY ADDING [RJ TO IT © COMPUTE THE IN-PLANE , STIFFNESS MATRIX [k ] I) COMBINE THE TWO MATRICES (kpj" & C^ b] TO GET 18X18 [k,] IN LOCAL CO-ORDINATE. TRANSFORM THE [ f e j MATRIX TO GLOBAL CO-ORD - 1 6 2 -rC CONVERT THE CONSTANT BAND WIDTH [K] TO A VAR-IABLE BAND WIDTH [Kj INITIALIZE MASTER LOAD VECTOR TO ZERO. READ & WRITE -LOADS & IF NECESSARY TRANSFORM THEM TO GLOBAL CO-ORD-FOR EACH BEAM STIFFENER. TRANSFORM THE NODAL CO-ORDS. DEFINING THE BEAM'S ORIENTATION TO THE LOCAL SYS & COMPUTE ITS LENGTH. COMPUTE THE 12X12 BEAM STIFFENER STIFFNESS MATRIX (R b s) & TRANSFORM TO GLOBAL SYSTEM. BUILD [Kj BY ADDING ( k b s ) TO IT. DECOMPOSE (INVERT) THE MASTER STIFFNESS MATRIX [K] & IF MORE THAN ONE LOAD CASE SAVE. COMPUTE THE DEF-LECTIONS OF THE NODES FROM PK-1 = 8 © FOR EACH FINITE ELEMENT TRANSFORM THE GLOBAL ^ \ (CO-ORDINATE OF THE 3 NODES ] V TO THE LOCAL SYSTEM. J ( FROM THE SOLUTION VECTOR OF DISP. IN OBAL CO-ORDINATE, COMPUTE THE STRESSES TRANSFORM SOLUTION VECTOR (<5) TO LOCAL SYSTEM SEPARATE (6) INTO ITS 9 INPLANE (6 ) & 9 BENDING COMPONENTS COMPUTE THE BEND-ING STRESSES ( 0 , ) D /^T COMPUTE THE AVERAGE >v ' / BENDING STRESSES, PLANE \ I STRESSES (C.S.T.) & PLANE J V STRESSES-(LST) AT EACH / NODE >/ COMPUTE THE PLANE STRES-SES USING ONLY u & V AT EACH NODE (C.S.T.) COMPUTE THE DISPL. OF MID-SIDE NODES BY USING THE END NODE DISPLS & HERMITIAN POLYNOMIALS FOR A STDF.-COMPUTE THE PLANE STRESSES USING L.S.T. _ FORMULATION STOP -164-BIBLIOGRAPHY Beavers, J.E., " T h e o r e t i c a l and Experimental I n v e s t i g a t i o n s of Non-P r i s m a t i c Folded P l a t e S t r u c t u r e s " , Ph.D. Thesis presented to V a n d e r b i l t U n i v e r s i t y , N a s h v i l l e , Tennessee, 1974. Cook, R.D., "Concepts and A p p l i c a t i o n s of F i n i t e Element A n a l y s i s " , John Wiley and Sons, I n c . , Toronto, 1974. Cowper, G.R. , L i n d b e r t , G.M., Olson, M.D., "A Shallow S h e l l F i n i t e Element of T r i a n g u l a r Shape", I n t e r n a t i o n a l J o u r n a l of S o l i d s and S t r u c t u r e s , V o l . 6, 1970, pp. 1133-1156. Cowper, G.R., Kosko, E., Lindberg, G.M., Olson, M.D., "A High P r e c i s i o n T r i a n g u l a r P l a t e Bending Element", NRC No. 10621, Ottawa, 1968. Johnson, CD., and Lee, T.T. , "Experimental Study of Non-Prismatic Folded P l a t e s " , J o u r n a l of the S t r u c t u r a l D i v i s i o n , ASCE, No. ST6, June 1968, pp. 1441-1455. Olson, M.D., " A n a l y s i s of A r b i t r a r y S h e l l s Using Shallow S h e l l F i n i t e Elements", T h i n - S h e l l S t r u c t u r e s : Theory, Experiment and Design, Fung, Y.C and Sechler, E.E. , P r e n t i c e - H a l l , Inc., New Jersey, 1974. Olson, M.D., " C o m p a t i b i l i t y of F i n i t e Elements i n S t r u c t u r a l Mechanics" Proceedings of World Congress on F i n i t e Element Methods i n S t r u c t u r a l Mechanics", Bournemouth, England, Oct. 12-17, 1975. Timoshenko, S., "Strength of M a t e r i a l s " , D. Van Nostrand Co., New York, 1958. Z i e n k i e w i c z , O.C, "The F i n i t e Element Method i n Engineering Science", Second E d i t i o n , McGraw-Hill Book Co., London 1971. Clough, R.W., Johnson, R.J., "A F i n i t e Element Approximation f o r the A n a l y s i s of Thin S h e l l s " , I n t . J . S o l i d s S t r u c . , No. 4, pp. 43-60, 1968.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A new eighteen parameter triangular element for general...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A new eighteen parameter triangular element for general plate and shell analysis Bearden, Terrance William 1976
pdf
Page Metadata
Item Metadata
Title | A new eighteen parameter triangular element for general plate and shell analysis |
Creator |
Bearden, Terrance William |
Date Issued | 1976 |
Description | The purpose of this investigation was to develop an eighteen parameter flat triangular finite element for analyzing plate and shell structures. The development of the element was accomplished by combining a plate bending element with a new plane stress element. The well known nine parameter triangle using the normal displacement and two slopes at each vertex was used for the plate bending element. This element contains an incomplete cubic for the normal displacement. For the in-plane element, complete cubics were used initially for the displacements and then various constraints were imposed to reduce the number of generalized co-ordinates to nine, namely the two in-plane displacements and an in-plane rotation at each vertex. One of the constraints, namely that the included angle at each vertex was invariant, destroyed the completeness of the element. However, the element was compatible in the plane. A patch-type test of the in-plane element showed that it could not represent all constant strain states exactly. However, the errors were small. The complete element was then tested on a plane stress cantilever beam, a square plate subjected to membrane stresses only, a cylindrical shell, a spherical shell and a non-prismatic folded plate structure. In all cases, reasonable engineering accuracy was achieved with modest grids of elements. Thus it was concluded that the incompleteness of the in-plane element was not too important. Finally, a compatible beam element was formulated and tested to supplement the triangular element. The beam element formulation included unsymmetric crosssections. |
Subject |
Plates (Engineering) Shells (Engineering) |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-02-09 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0062480 |
URI | http://hdl.handle.net/2429/19867 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1976_A7 B42.pdf [ 6.77MB ]
- Metadata
- JSON: 831-1.0062480.json
- JSON-LD: 831-1.0062480-ld.json
- RDF/XML (Pretty): 831-1.0062480-rdf.xml
- RDF/JSON: 831-1.0062480-rdf.json
- Turtle: 831-1.0062480-turtle.txt
- N-Triples: 831-1.0062480-rdf-ntriples.txt
- Original Record: 831-1.0062480-source.json
- Full Text
- 831-1.0062480-fulltext.txt
- Citation
- 831-1.0062480.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062480/manifest