FINITE ELEMENT ANALYSIS OF ROTATING DISKS by GREGORY (B.A.Sc, LYNN NIGH The University of British Columbia, 1 978) A THESIS SUBMITTED IN PARTIAL FULFILMENT OF T H E REQUIREMENTS FOR T H E DEGREE OF MASTER OF APPLIED SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Civil Engineering We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA May 1980 © Gregory Lynn Nigh ,1980 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copy- ing of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It is under- stood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Gregory Lynn Nigh Department of Civil Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, V6T 1W5 Canada ii ABSTRACT The finite element method is applied to the analysis of rotating disks. The equation governing the transverse deflection of a rotating disk is presented in both body-fixed and space-fixed coordinate systems. The quadratic stress triangle is used to determine the i n - plane stress distribution which contributes to the out-of-plane stiffness. The out-of-plane problem is analyzed using elements based on the high precision 18 degree of freedom compatible plate bending triangle. Free vibration problems of both axisymmetric and non-axisymmetric disks are solved in the body-fixed coordinate system. The numerical results obtained for axisymmetric disks show good agreement with results from other analyses. In a space-fixed coordinate system, the response of a rotating axisymmetric disk to a space-fixed transverse load is examined in some detail. It is found that there exists a fundamental critical rotation speed at which a component of a natural vibration mode becomes stationary in space, leading to an instability. This critical speed may be related to a non-dimensional stiffness parameter which, together with the non-dimensional hub radius, completely defines the disk geometry and stiffness. T h u s , it appears that the critical value of the stiffness parameter is a function only of the non-dimensional hub radius. The response of the disk past this iii critical stiffness is also examined. The results provide evidence of the existence of the higher critical speeds and indicate that the response is very sensitive to any approximations made with respect to the geometry of the disk or the in-plane stress distributions. iv T A B L E OF C O N T E N T S Page ABSTRACT ii T A B L E OF C O N T E N T S iv LIST OF T A B L E S vii LIST OF FIGURES ix ACKNOWLEDGEMENTS xi Chapter I. II. III. INTRODUCTION 1 1.1 Background 1 1.2 Purpose and Scope 3 THEORETICAL BACKGROUND 5 2.1 Introduction 5 2.2 Body-Fixed Coordinate System 5 2.3 Space-Fixed Coordinate System 8 2.4 Energy Functionals and Boundary Conditions 12 FINITE ELEMENT REPRESENTATION 18 3.1 Introduction 18 3.2 In-Plane Problem 19 V Chapter Page 3.3 3.4 IV. 3.2.1 Introduction 19 3.2.2 The Finite Element Solution 19 3.2.3 The Forced Exact Solution 25 Out-of-Plane Problem 27 3.3.1 Introduction 27 3.3.2 Linear Plate Bending Stiffness 28 3.3.3 Geometric Stiffness 31 3.3.4 Inertia Stiffness 34 Convergence and Bounds 38 SAMPLE PROBLEMS 41 4.1 Introduction 41 4.2 Test Problems 43 4.2.1 Geometric Stiffness Matrix 43 4.2.2 Combined Geometric and Inertia Stiffness Matrices 4.3 44 Vibrations of Axisymmetric Disks 50 4.3.1 Vibrations of a Rotating Membrane 51 4.3.2 Vibrations of a Rotating Disk-Without Hub Vibrations of a Rotating Disk-With 58 Hub 63 4.3.3 4.4 Stiffness Analysis of an Axisymmetric Disk . . . . 67 4.5 Analysis of a Non-Axisymmetric Disk 82 4.5.1 Static Deflection Response 82 4.5.2 4.5.3 Free Vibration Response In-Plane Stresses 84 89 vi Chapter V. Page CONCLUSIONS REFERENCES 96 99 APPENDIX A - MATRICES A S S O C I A T E D WITH T H E SOLUTION OF T H E IN-PLANE PROBLEM 101 APPENDIX B - MATRICES A S S O C I A T E D WITH T H E SOLUTION OF T H E O U T - O F - P L A N E PROBLEM 107 vii LIST OF T A B L E S Table 1.1 Page Non-Dimensional Frequency Coefficients for Membrane Vibration 45 4.2 Thin Rotating Ring Deflection Results 48 4.3 Non-Dimensional Natural Frequencies of a 4.4 Rotating Membrane Natural Vibration Frequencies (Hz) of a Rotating Disk 4.5 4.6 Non-Dimensional Frequency Parameters of a Rotating Disk Critical Stiffness from Finite Element Solution K = 0.2 4.7 Natural Frequencies of a Disk Rotating @^£p-. 4.8 Natural Frequencies of Non-Axisymmetric Disk 57 62 66 7Z * 7 6 90 A.l Transformation Matrix [T] 104 A. 2 Rotation Matrix [R] 105 A. 3 Strain Distribution Matrix [AA] 106 A. 4 Transformation Matrix for Forced Exact Stress Distribution [T] 106 Table Page B.1 Transformation Matrix [T] 110 B.2 Rotation Matrix [R] 111 ix LIST OF FIGURES Figure 3.1 3.2 •4.1 4.2 Page Coordinate Systems and Body Forces for Quadratic Stress Triangle Element 21 Coordinate Systems for 18 DOF Triangular Plate Bending Element 29 Layout of Finite Elements for Membrane Vibration 45 Vibrating Square Membrane, Relative Error in Eigenvalues 46 4.3 Thin Rotating Ring 48 4.4 Finite Element Grids for Disk Without Hub 53 4.5(a) Non-Dimensional In-Plane Radial 4.5(b) Distribution K = 0.0 Non-Dimensional In-Plane Circumferential Stress Distributions K = 0.0 4.6(a) 4.6(b) 4.7 Stress 54 55 Vibrating Rotating Membrane, Relative Error in Eigenvalues 59 Vibrating Rotating Membrane, Relative Error in Eigenvalues 60 Finite Element Grid for Disk with Hub 65 X Figure Page 4.8 Finite Element Grids for Stiffness Analysis 4.9(a) Non-Dimensional In-Plane Radial Stress Distributions 4.9(b) 4.10 4.11(a) 69 K = 0.2 70 Non-Dimensional In-Plane Circumferential Stress Distributions K = 0.2 71 Typical Deflected Shape Near Fundamental Critical Speed 74 Deflection Contours, Grid N = 3, Stresses Forced Exact 79 4.11(b) Deflection Contours 4.12 Non-Axisymmetric Disk Geometry 4.13 Finite Element Grid for Non-Axisymmetric Disk 4.14(a) Static Deflection of a Non-Axisymmetric Disk 4.14(b) Static Deflection of a Non-Axisymmetric Disk 4.14(c) Static Deflection of a Non-Axisymmetric Disk 4.15 Finite Element Grid for Non-Axisymmetric Disk 4.16 Typical Vibration Modes of a Non-Axisymmetric 4.17 Disk Non-Dimensional In-Plane Stress Distributions of a Non-Axisymmetric Disk 92 Non-Dimensional Stresses Along Slot Diameter .. 94 4.18 a = 0.001 80 83 .. 83 .... 85 86 .... .. 87 88 91 xi ACKNOWLEDGEMENT The author wishes to express his gratitude to his advisor, Dr. M.D. Olson, for his interest and patient guidance and advice during the preparation of this thesis. Dr. S.G. He would also like to thank Hutton for reviewing this work. The financial support of the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. Finally, the author would like to thank his wife, Janice, for her encouragement and support throughout his graduate work. 1 CHAPTER I INTRODUCTION 1.1 Background Rotating disks are a fundamental component of many machines and mechanisms. Often these disks are subject to a loading or excitation in the transverse (out-of-plane) direction. It then becomes necessary to determine the natural frequencies of vibration in the transverse direction, or the response to a transverse load. The problem is complicated by the fact that the disk is made stiffer due to the effects of rotation. In what appears to be one of the earliest works in this field. Lamb and Southwell [1] treated the problem of a complete ( i . e . , no hub) constant thickness axisymmetric disk in which the flexural rigidity was negligible and the transverse stiffness was derived solely from the membrane stresses due to rotation. A simple, exact, expression was derived for the natural frequency of vibration that depended only on the rotational speed, Poisson's ratio, and the number of nodal diameters and circles which characterize the particular mode of interest. In addition, for the case when flexural rigidity is not negligible, upper and lower bounds were obtained for the natural frequencies of vibration corresponding to some modes. Southwell [2] extended these results to the case in which the disk has a central hub. 2 The in-plane stresses f o r a complete d i s k were u s e d , b u t the o u t - o f - plane d i s p l a c e m e n t a n d slope were c o n s t r a i n e d at the h u b A g a i n , approximate results were o b t a i n e d f o r the case in which b o t h f l e x u r a l e f f e c t s a n d membrane s t r e s s e s stiffness. E v e r s m a n a n d D o d s o n [3] c o n t r i b u t e d to the o u t - o f - p l a n e c o n s i d e r e d the f r e e v i b r a t i o n s of an a x i s y m m e t r i c c o n s t a n t t h i c k n e s s d i s k hub, radius. transverse with a c e n t r a l with b o t h i n - p l a n e a n d o u t - o f - p l a n e d i s p l a c e m e n t s f u l l y c o n - s t r a i n e d at the h u b . T h e i r a n a l y s i s was f r e e of a n y approximations to the t h e o r y , a l t h o u g h a n u m e r i c a l solution of the d i f f e r e n t i a l e q u a t i o n was required. N o n - d i m e n s i o n a l g r a p h i c a l r e s u l t s were p r e s e n t e d which e n a b l e d e t e r m i n a t i o n of the n a t u r a l v i b r a t i o n f r e q u e n c i e s c o r r e s p o n d i n g to modes c h a r a c t e r i z e d b y 0, 1, or 2 nodal diameters or c i r c l e s . Mote [4] a n d Mote a n d Nieh [5] h a v e p r e s e n t e d both theoretical a n d e x p e r i m e n t a l i n v e s t i g a t i o n s critical speed instability. into the phenomenon of T h i s i n s t a b i l i t y o c c u r s when the wave s p e e d of the b a c k w a r d t r a v e l l i n g component of one of the v i b r a t i o n modes is the same as the rotation s p e e d . v i b r a t i o n mode which is s t a t i o n a r y leads to a in s p a c e , a n d is an -factor in the d e s i g n of thin r o t a t i n g In This important disks. r e c e n t y e a r s the f i n i t e element method has been a p p l i e d to the problem of f r e e v i b r a t i o n s of r o t a t i n g d i s k s . a n d Wilson resonant [6] a n d K e n n e d y a n d Gorman [7] Kirkhope have used efficient a n n u l a r elements to examine the v i b r a t i o n s of t h i n a x i s y m m e t r i c disks. 3 The effects of variable thickness and temperature distributions are also included in these analyses. Mote [8] has considered the problem of vibrations of non-axisymmetric (slotted) disks by using a twelve degree of freedom compatible plate element with a constant strain element to determine the membrane stresses due to rotational and thermal effects. The response of a rotating disk subject to a stationary verse load has received relatively little attention. trans- Mote [4] has discussed it briefly in the context of critical speed instability. Benson [9] and Benson and Bogy [10] have presented an analysis of this problem relative to a space-fixed coordinate system. They used a finite differ- ence scheme (in the radial direction) to obtain the steady-state response of a rotating axisymmetric disk to a space-fixed constant point load, for various rotation speeds. However, they did not observe the phenomenon of a critical speed instability as would be expected. 1. 2 Purpose and Scope In this thesis, the finite element method is applied to the analysis of rotating disks. Since it is desirable to be able to consider non-axisymmetric problems, high precision triangular elements are used. Both the in-plane stresses due to rotation and the out-of-plane response are determined by this method. A geometric stiffness matrix is derived which incorporates the in-plane stresses in the out-of-plane problem, and numerical results for vibrations of axisymmetric disks are obtained 4 which confirm the validity of the formulation. In addition, a special stiffness matrix is derived for the analysis of the steady-state out-ofplane deflection of an axisymmetric disk due to a transverse load. This problem is examined in some detail. axisymmetric problems are also examined. Some results from non- 5 CHAPTER THEORETICAL 2.1 II BACKGROUND Introduction In this chapter, the governing differential equation for the transverse deflection of a rotating disk will be developed. The problem of free vibrations of a rotating disk is most easily solved in a bodyfixed coordinate system. However, to treat the problem of a rotating disk subject to a constant point load which is stationary in space, the governing differential equation is transformed to a space-fixed coordinate system. The associated energy functionals and boundary conditions are also examined. The mid-surface of the disk is assumed to lie in the x - y plane. Rotation is about the z axis, and the angular velocity, Q, is constant. The disk thickness is assumed to be small, so that the effect of the transverse shear stresses on the deflected shape may be neglected. The material is linearly elastic, homogeneous and isotropic, and may be characterized by a modulus of elasticity, E, a Poisson's ratio, v, and a constant mass density, p. 2. 2 Body-Fixed Coordinate System Consider a disk of constant thickness, rotating at a constant speed, which is not subject to any external forces. As a result of 6 the rotation, a state of stress will exist in the disk. If the thick- ness of the disk is small relative to the other dimensions, the stresses may be assumed to be constant through the thickness. T h e n , to solve for the state of stress within the disk, the equilibrium equations from plane stress elasticity may be used, with the body forces equal to the inertia force. T h u s , in polar coordinates; pQ S = 0 r 2 . . . . (2. 1) in which R and S are the radial and circumferential body force components per unit volume respectively, and r is the radius. If the disk is axisymmetric, the state of stress in the plane may be easily solved for in a polar coordinate system. However, the more general case of a non-axisymmetric disk is to be considered. For this reason, the problem is solved in an x-y coordinate system, so that the body forces are; X = pfi x Y 2 = pQ y 2 . . . .(2.2) and the differential equations for in-plane equilibrium are; 3a 3a * + 3x 3a 2£ 3x + X = 0 3a 2- + 3y Y = 0 ^1 3 y + . . . .(2.3) 7 in which a , a and o are the usual in-plane stress components, x y xy The response of the disk in the transverse direction, whether it be a free vibration response or the response to a transverse load, will be affected by these in-plane stresses and body forces. T h u s , the differential equation of the deflected surface of a disk of constant thickness, when there are both stresses and body forces acting in the plane of the disk, is, from Reference [ 1 1 ] ; 3w _ 3w 3w + 2a + a xy ' x 3 y y „3y , 2 DV w-h o x 3x ... 2 2 2 +h 3J5L 3x x + Y 3y = q w2 (2.4) in which; w transverse deflection h thickness D plate rigidity V biharmonic operator q transverse load per unit area. Eh" I2(1-v ) z The in-plane stresses and body forces must satisfy the equilibrium relations of Equation 2 . 3 . To consider the problem of transverse vibrations of the disk, a transverse inertia force must be added to Equation 2.4 to render 3 w 2 D V w-h ' x 3x 2 3w 3 w + 2a + a xy 3 x 3 y y 3y 2 2 2 . , 3 w ,, 3 w 3 w +h X +Y = q - ph 3t 3x 3y 2 v 2 (2.5) 8 where 9 2 w 3t is the acceleration in the transverse direction. 2 Equation 2.5 is the differential equation of the deflected surface of a rotating disk, relative to a body-fixed coordinate system. That is, the x and y coordinate axes rotate in space at a constant speed, ft, with the disk. T h e n , because the in-plane body forces are time independent, the in-plane stresses will also be time independent and hence, functions only of position. The load q and deflection w are the only time dependent variables in Equation 2.5. 2.3 Space-Fixed Coordinate System If the transverse load, q, is stationary in space, then in a body-fixed coordinate system, it must be a function of time. is time dependent. Equation 2.5 is difficult to solve. If q Furthermore, if the load is stationary in space, it is likely the response in space which is of interest. T h u s , even if the response, w, could be obtained from Equation 2.5, it would have to be expressed relative to a coordinate system fixed in space. The alternative, which is described here, is to express Equation 2.5 relative to a space-fixed coordinate system. In order to make the transformation from a body-fixed coordinate system to a space-fixed coordinate system, the material derivative operator is used to evaluate the time derivative terms of Equation 2.5. by: For two dimensional problems, this operator is defined where v x and are the velocities of the body in the x and y directions, respectively. v x = For the problem described herein; -fiy v y = (2.7) fix where x and y are the coordinates, in space, of a point on the body. Then using Equation 2.6 to evaluate the transverse inertia force gives D , , , . 9w —=• (phw) = ph 9t Dt' , r, + 2phfi 2 ( 2 9 w 3x' - phfi : 9w 9x 9w 9x 9t 2 T {' y 3w X 29w +x 9x9y 9y 2 2xy 3w 9y at J 2 + 9w 2 (2.8) The inertia term of Equation 2.5 is then replaced by Equation 2.8 and the relations of Equation 2.2 are used to yield the transformed Equation 2. 5 as 10 3w _ + 2a 2 DV w - h h j 3x X 2 X 2 3w ( 3w 2 + a V 3x3y y 3y i 2 3w 2 2 + phfi' 3w ' 2 - 2xy 3x' 2 3 w 2 + x 3x3y . 3 w - . _. = q - ph + 2phQ 3t 3y' 3w 3w 3x3t 3y3t 2 2 2 (2.9) This is the differential equation of the deflected surface of a rotating disk in a space-fixed coordinate system. Attention here is focused on a particular point in space which is momentarily occupied by the disk, rather than on a particular point on the disk, as it is in a bodyfixed coordinate system. In this coordinate system, only the steady state problem will be considered, so that the time derivative terms on the right hand side may be neglected, yielding; 3 w 4 * 3x 2 3w 2 2 DV w - h + 2a xy '2 3 w 3x3y : y 1? 3 w y 2 2 + phfi 3w 2 + a 2xy 3x3y 3y > 2 2 3w 2 + x 3y = q (2.10) 2 This equation is applicable only to disks which are axisymmetric with respect to both their geometry and in-plane stresses. The basis for 11 this restriction may be realized if a disk with a radial slot is considered. Then at any point in space within the domain of the problem, the boundary conditions are constantly changing. When the free edges associated with the slot are at the particular point in space, free edge stress boundary conditions apply, but for the remainder of the revolution, there are no prescribed boundary conditions at that point. In addition, in a non-axisymmetric disk, the in-plane stresses, expressed as functions of coordinates in space, must also be functions of time. Only in the case of an axisymmetric disk are the in-plane stresses and boundary conditions independent of the choice of coordinate system. A rigorous derivation of the governing equations of a non-axisymmetric disk in a space-fixed coordinate system is beyond the scope of this work. For the problems considered herein, it is sufficient to note that this restriction on Equation 2.10 exists. Equations 2.5 and 2.10 both describe the deflection surface of a rotating disk with in-plane stresses and body forces. Their difference is due only to the coordinate system to which they are related. However, each shows certain advantages in the various types of problems to be solved. For solving free vibration problems, in which q = 0, Equation 2.10 contains terms which have mixed time and space derivatives, whereas Equation 2.5, which is suited to this type of problem, has separate time and space derivative terms. To determine the steady-state response to a transverse load, in which the time derivatives are zero, the situation is reversed. In Equation 12 2.10, w h i c h is r e l a t i v e to a s p a c e - f i x e d c o o r d i n a t e s y s t e m , q is a f u n c t i o n o n l y of p o s i t i o n , w h e r e a s it is a l s o a f u n c t i o n of time in E q u a t i o n w h i c h makes the p r o b l e m more d i f f i c u l t . Thus, both equations 2.5, will be u s e d in t h i s t h e s i s , d e p e n d i n g on the p a r t i c u l a r p r o b l e m to be solved. 2.4 Energy Functionals and Boundary Conditions A s the f i n i t e element r e p r e s e n t a t i o n s of the p r o b l e m s d e s c r i b e d h e r e i n r e q u i r e the use of the e n e r g y functionals, they are presented h e r e . They a r e b a s e d on s t r a i n e n e r g y c o n s i d e r a t i o n s a n d the p r i n c i p l e of minimum p o t e n tial e n e r g y . T h e s e f u n c t i o n a l s a r e also u s e d to d e r i v e the b o u n d a r y d i t i o n s a s s o c i a t e d w i t h the d i f f e r e n t i a l e q u a t i o n s g o v e r n i n g con- the o u t - o f - plane problem. When the i n - p l a n e a n d o u t - o f - p l a n e s y s t e m s a r e u n c o u p l e d , a n d the r o t a t i o n s p e e d is c o n s t a n t so t h a t the i n - p l a n e b o d y f o r c e s a r e i n d e p e n d e n t of time, the i n - p l a n e p r o b l e m may be t r e a t e d as a s t a t i c problem. T h e n the total p o t e n t i a l e n e r g y of the i n - p l a n e s y s t e m , n, is the s t r a i n e n e r g y , U, plus the p o t e n t i a l e n e r g y of the b o d y f o r c e s , W. Thus n = (2.11) u fw where U Eh u 2(1-v ) 2 2 x + v 2 y + 2vu v x y + 1-v, 2 (u y +v x ,2 dxdy ) . . .(2.12) 13 Xu + Yv W (2.13) dxdy in which u and v are the in-plane displacements in the x and y directions, respectively, and subscripts denote partial derivatives. ln the case of dynamic systems such as described by Equations 2.5 and 2.9, the rigorous formulation of the governing functionals is based on Hamilton's principle. r*1 L dt The principle may be written as = (2.14) 0 where L is the Lagrangian function which consists of the strain energy of the body, the potential energy of the loads and the kinetic energy of the moving body. In addition to satisfying the displacement boundary conditions as in potential energy functionals, the variations of the exact displacements must also be zero at the two arbitrary instants of time, t^ and t^. The strain energy in the body consists of two terms, Ug and U.., which are the strain enerqies associated with the flexure of the plate and the transverse components of the in-plane stresses, respectively. They are derived in Reference [11] and are; U U B M D 2 2 2 +w + 2vw w + 2 ( 1 - v ) w ^ dxdy [ xx yy xx yy xy j w 2 2 aw +aw +2a w w , dxdy x x y y xy x y j wl (2.15) (2.16) 14 The potential energy of the lateral loads is given b y ; W - - q w dx dy . . . The kinetic energy associated with transverse motion of the disk .(2.17) is; 2 1 2 T ph \ ^ } dxdy where . . . is the transverse velocity of the disk t Thus, the Lagrangian function is; .(2.18) in body fixed coordinates, a L = U B + U M - T + W . . . .(2.19) and it may be shown that application of Hamilton's principle leads to Equation 2.5. When solving free vibration problems with Equation 2.5 the lateral loads are, by definition, zero. The corresponding Lagrangian then consists only of the strain energies and the kinetic energy. If a harmonic time variation of the deflections is then assumed, the Lagrangian may be integrated with respect to time to yield a functional form suitable for use with the finite element method. In order to derive the functional corresponding to Equations 2.9 or 2.10, the material derivative of Equation 2.6 is used to describe the transverse velocity in a space fixed coordinate system. Thus [ 3 w1 f Dw 1 •j > of Equation 2.18 is replaced by "j p~f~ f • The kinetic energy * then consists of two parts; T , which contains time derivative terms, and U,, which does not. Thus; 15 T (2.20) + U, Using the relations of Equation 2.7, the inertia integral, U., phfi U, 2 2 2 2 . y w + x w - 2xyw w x y x y : 1 7 is; dxdy . . . .(2.21) and ph where w o . t 7 x y dxdy .(2.22) The equation of motion. Equation 2.9, may then be derived from Hamilton's principle; (U B +U M - U, - T + W) dt = 0 (2.23) As only the steady state problem is to be considered here, it may be shown that Equation 2.10 may be derived from the first variation of the functional II, where 11 = U B + U M " U l + W .(2.24) 16 The boundary in the same way. integrals for Equation 2.10 may also be obtained The resulting boundary integrals (which must vanish), in x - y coordinates, are; 9 6w 2 2 phft y w x M D( w x D (w 9 6w 2 2 phQ x w D(w yyy yy + N w x x + N w xy y ) , dy = 0 yy J + v w ) , dy = 0 yy j xx H - D(w 7 + (2-v) w xxx 9 6w 9 6 w 2 + phft xyw y 7 y 2 + phft xyw x + N w y y + N w xyx + (2-v) w ) , dx = o **y J + vw xx ) (2.25) dx = 0 where N , N , and N are the in-plane stress resultants and 6w, x y xy ^ 6w , and 6w x y are the arbitrary variations of displacement and slope which are kinematically admissable functions of x and y . It is not 2 obvious that the phft terms corresponding to the displacement variation vanish on the boundary of the disk as would be expected. In order 17 to confirm that they do, the line integrals associated with 6w expressed in polar coordinates. must be As Equation 2 . 1 0 is applicable only to axisymmetric disks, this is not a restriction. T h u s , the boundary integrals become; 1 o 6 w |i N w + N „ r r 'V9 r J I M ^ d0 == r <> 6 w J r ^ where Q r w "6 Q 1 d C) ;r r + r Q M de 3 0 . . . .(2.26) 0 J is the transverse shear and M twisting moments, respectively. and M „ are the bending and r0 r a The vanishing of these boundary integrals then leads to the general boundary conditions along any edge as; either 1 Nw +N Q - w r r r 6 r 0 Q - Q 1 8M re + r r 36 = 0 or 6w = 0 and either M = r 0 or 6w r = 0 . . . . .(2.27) For the disks considered here, the inner radius is clamped, allowing no deflection or normal slope, so that outer radius, the boundary is free, so that N 6w = 6 w and N r r Q r = 0 . On the are both zero as a result of the natural boundary conditions in the in-plane problem. Then in the out-of-plane problem, the natural boundary conditions on the free edge are; 3M -Q + r p JQ^- = 0 and M r = 0 . . . . .(2.28) 18 CHAPTER FINITE ELEMENT 3.1 III REPRESENTATION Introduction The finite element method has proven to be a powerful technique for obtaining an approximate numerical solution to problems in which the governing equations, problem geometry or boundary conditions are so complex that an analytical solution cannot be extracted. The theoretical foundations of this method, and the requirements for convergence of the solution have been well developed. This chapter will describe the application of the finite element method to the problem of rotating disks. The elements used herein are high precision triangular elements. Their stiffness matrices are derived by using an approximate displacement field to evaluate the strain energy in each element. The global coordinate origin is assumed to coincide with the center of rotation. The use of triangular rather than annular or sectorial elements enables non-axisymmetric problems to be considered. Since the in-plane stresses appear in the out-of-plane differential equation and energy integrals, some distribution of these stresses within each element must be determined before the out-of-plane problem may be solved. T h u s , the total finite element problem is broken into two separate problems, the in-plane problem and the out-of-plane problem. Of these two, the in-plane problem must be solved first. 19 3. 2 In-Plane Problem 3.2.1 Introduction The in-plane stresses as used in the out-of-plane problem may be determined by one of two methods. If the geometry of the problem is axisymmetric, there exists an exact solution for the in-plane stresses cy, and based on plane stress elasticity theory. Within each element this stress distribution may be approximated by a quadratic function that is forced to match the exact solution at each node and mid-side. Alternatively, an approximate stress distribution may be obtained by solving the in-plane problem via the finite element method. This is the only alternative if the disk geometry is not axisymmetric. 3.2.2 The Finite Element Solution The in-plane finite element problem is solved using a quadratic stress (strain) triangle (QST) element, which is a plane stress triangle with a complete cubic variation of both u and v displacements. The cubic variation of displacements allows a quadratic variation of strains (and hence, stresses). There are eighteen degrees of freedom associated with this element, namely, u ' ' 3u 3x ' 3u dy ' v at each of the three nodes. ' 3v 3x ' 3v 3y 20 The element is derived in a local E, - n coordinate system as shown in Figure 3.1. The in-plane displacements relative to this local system may be expressed as series functions of these local coordinates and the generalized coordinates a.; 20 u = m. I a. E, 20 n. n v i=1 = I p. aj E, q. n . . . .(3.1) i=1 where m., n., p. and q., integer are given in Appendix A . values ranging from 0 to 3 inclusive, A matrix relationship between the nodal displacements (in £-n. coordinates) and the generalized coordinates, may then be developed; w' where {A} - \ = [T] {A [ . . . is the vector of generalized coordinates, [T] is a tion matrix which is a function only of the element .(3.2) transforma- local coordinates a, b, and c, and {w'"} is the nodal deflection vector, including the displacements of the centroid. A rotation matrix [R] may also be derived that transforms the displacements in local coordinates to the displacements in global coordinates; W where {W } C L 1 = [R] (w C \ . . . .(3.3) is the deflection vector in global coordinates, still, including the centroidal displacements. 21 x,u Figure 3.1 Coordinate Systems and Body Forces for Quadratic Stress Triangle Element 22 Using the series form of u and v , the strain energy of the element may be evaluated in closed form, using Equation 2.12 . It may be expressed in the required quadratic form of aj as U = 1 {A} [k] T {A} . . . . .(3.4) This stiffness matrix [k] is a function of the material properties and element dimensions, and is reproduced, with the transformation and rotation matrices, in Appendix A. The strain energy may be related to the global deflection vector through the transformation and rotation matrices; U = 1 {W } G [R] [T _ 1 ] [k] [T _ 1 ] [R] {W } G . . . .(3.5) so that the stiffness matrix relative to the global coordinates is; [K] = [R] [T '] [k] [T '] [R] . . . . .(3.6) In order to solve the in-plane problem of rotating disks, a consistent body force load vector also needs to be developed for each element. As mentioned in Chapter 2 , the only body force acting in the plane is R = p f t r , which may be expressed in x - y coordinates as 2 X = p f i x and Y = pQ y. 2 2 These x - y body forces may be transformed to element local coordinates to give 23 x = pfi (c 2 ? + O Y = p f i ( c ^ + n) 2 in which the constants (x,cos6 + y,sin8) c = (3.7) ( - x ^ i n S + y^osB) where x , y , x^, y^ are the x - y nodal coordinates of nodes 1 and 3 respectively and 9 is the inclination of the local coordinate axes as shown in Figure 3.1. The potential energy of the body forces may then be evaluated from Equation 2.13, using the assumed displacement fields u and v and the body forces of Equation 3.7, to give 1 f 0 W=p h f i 2 i=1 ' 20 Z a. c m. n. m.+l n. 'n + E n ') d^dri 1 1 u + p h f i T a. 2 i=1 (c ti (c^ p. q. p. q.+ l 'n + K 'n )d^dn 1 (3.8) 1 1 Recognizing Euler's Beta function; F(m,n) K n d^dn : n+1 ( m+1 , u^ | " " a ( b ) m + 1 l m!n! j (m n+2) (3.9) + the potential energy of the body forces may be expressed as the vector product 24 W = {F } .(3.10) {A } T in which the elements of the vector {F} are given by F. = phft 2 = ph fi 2 CpF(m.,n.) + F(m.+1,n.) i - 1,2, — , 1 0 c F(p.q.) + F(p.,q. + 1) i = 11,12, ,20 .(3.11) Applying the same rotation and transformation matrices as before gives the element consistent load vector {P} = [R] T [T ] ] {F} in global coordinates; . (3.12) . The stiffness matrix and load vector are normally statically condensed to remove the degrees of freedom associated with the centroid displacements. If these centroid displacements have no associated load components, this procedure is straightforward. However, when such loads exist, special care must be taken when performing the static condensation and particularly when evaluating the strain energy, to include all the required terms. Once the global problem has been solved, the centroidal displacements of each element may be determined from the element nodal displacements. With the full element displacement vector {W Q }, the 25 generalized coordinates within each element may another transformation; {A} = . (3.13) [T ][R]{W } _ 1 C Finally, the quadratic stress distribution, may be determined through with respect to the local axes, be calculated as a function of the generalized coordinates through 30 35 [E] 3v [E] [AA] 3n 3u 3.2.3 The n which is 2 [AA] reproduced in is a function of Appendix A. Forced Exact Solution If the problem geometry is axisymmetric, become ( 3 . 14) K n 3v 95 where [E] is the elastic constitutive matrix and the generalized coordinates e the in-plane stresses functions only of the radial dimension, and are independent of the angular coordinate. identically zero. displacement As a result, the in-plane shear stress is If the boundary conditions require at the inner radius, a, and no surface tractions on the outer radius, b, the in-plane stresses are described a r = C (3 + v) 1 that there be by; no 26 Cjp 1 - C - -C ,2 + v) ' 1 + 3 v 3+v where r 1-v 3+v (3+v)-(1+v)K { (1+V)+(1-V)K J 2 and 2 2 (3.15) and K and L, are the non-dimensional hub radius, ^ , and radius, j-- , respectively. If the disk is solid ( i . e . , no hub), then both C and K are identically zero and the equations reduce to: a. - C^{3 + v) C^3 [1 - c ] 1+3v 3+v + v) (3.16) In order to generate the stress distribution, the value of a and o a is calculated at each node and mid-side of each element. These stresses are then transformed to the local £ - ri coordinate system. Then assuming a quadratic stress variation; 6 m. n. a = I i=1 a ; 'n jC 1 . . . .(3.17) 27 the stresses at each node and mid-side can be related to the coefficients of the quadratic through a transformation matrix, {a} = [T] {a} . . . .(3.18) This transformation matrix is, as before, only a function of the local coordinates and is the same for each of the stresses considered. It is reproduced in Appendix [a ,o r ,o r ) to be A. This is a more efficient technique of generating the in-plane stresses if the condition for axisymmetry is met. In the case in which there is no hub (K = 0), it also provides an exact representation of the in-plane stresses, since the variation of c? and Og is only quadratic r and can be described exactly by quadratic variations of , and a ^ - However, if K is non-zero, the exact solution cannot be described by 2 just a quadratic variation (due to the 1/e, terms) and so this technique will generate a stress distribution that is exact only at the nodes and mid-sides. 3. 3 3.3.1 Out-of-Plane Problem Introduction The out-of-plane problem as described by Equation 2.10 has three distinct relationships which contribute to the out-of-plane stiffness. The term D V w is associated with the linear bending 4 stiffness of a thin plate; the terms containing the in-plane stresses T ] 28 are associated with the geometric stiffness, and the terms proportional to phfi 2 are associated with what will be called the inertia stiffness. The element stiffness matrix corresponding to each of these terms will be derived separately, although they each use the same assumed displacement field and transformation and rotation matrices. 3.3.2 Linear Plate Bending Stiffness The stiffness matrix associated with the plate bending terms has been derived in detail in Reference [12], and its performance well documented. Its formulation is very similar to that of the considered Q S T . previously However, since it forms the basis for the other out- of-plane stiffness terms, a brief review of the derivation is presented here. Like the Q S T , this element is derived in a local E-ri ordinate system shown in Figure 3.2. co- In this case, the assumed out- of-plane displacement, w, is given by an incomplete quintic polynomial; 20 w = I i=1 m. a.E 1 n. n . . . 1 .(3.19) 1 where rrij and n., integer values ranging from 0 to 5 inclusive, are given in Appendix B. The quintic polynomial is made incomplete by deleting the E^ri (rrij = 4, nj = 1) term. In addition, the transformation matrix [T] is derived to include two constraint equations. These two c o n - straints, plus the missing E^n, force the normal slope along any edge 29 Figure 3.2 Coordinate Systems for 18 DOF Triangular Plate Bending Element 30 to vary only as a cubic. This has the combined affect of transforming from twenty generalized coordinates to eighteen nodal variables, and more importantly, enforcing inter-element compatibility, both in displacements and in slope normal to the element edge. The strain energy (Equation 2.15) may now be evaluated on the basis of this assumed displacement field, and expressed as a quadratic function of the generalized coordinates; U = \ {A} T [k] {A} . . . . .(3.20) The stiffness matrix [k] is again, a function of the material properties and element dimensions, and is reproduced in Appendix the coordinate transformation and rotation matrices as B . Applying before gives the element strain energy as a quadratic function of the eighteen nodal variables; T U = 1 {W } C T [R] T [T _ 1 ] The vector of nodal displacements,{W W ' 3 9w ' x 3 2 9w 3w y ' 9x2 ' at each of the three nodes. [k][T ][R]{W } _ 1 Q C . . . . .(3.21) } consists of 2 9w 3x3y ' 3 y 2 2 9w A consistent mass matrix and uniform load vector, based on the same assumed displacement field, may also be 31 derived. All of the matrices associated with this element are presented in Reference [12]. 3.3.3 Geometric Stiffness The stiffness associated with the in-plane stresses represents a first order of non-linearity in the system. This stiffness incorporates the effects of the transverse components of the in-plane stresses in the outof-plane problem. To derive the geometric stiffness matrix, the assumed displacement field and the in-plane stresses are written in a series form; 20 w = m. n. £ a. E n ' i=1 ; 1 ....(3.22) 1 N I o> = ^ b. E i=1 m. 1 n N n. a 1 ' = I c m. E ' i=1 1 n n. 1 N ov = m. n. T d. E n where a., m., and n. are as before. i' i' i . . . .(3.23) The variables b., c , and d. are i' i' the coefficients which describe the in-plane stresses ^ ov, a , and E n respectively, with positive directions as shown in Figure 3.2. i a r En The upper index, N, reflects the fact that any degree of variation of the in-plane stresses may be accommodated. However, as the QST was used to determine these stresses, and will allow, at best, a quadratic variation, N is set to six. Using a higher degree of variation would present no 32 conceptual difficulty, only the practical problem of actually generating the coefficients. Differentiating and squaring the expression for the displacement, and combining with the in-plane stress distributions, each term in the energy integral (Equation 2.16) may be expressed as a triple sum; 20 20 6 k~ 7 T T a.a.b. m.m. £ • • • • i i i J J i=1 j=1 k = 1 ' m - 8w 195 J 1 20 y 9w 9n • i 20 y • k 6 m y a.a.c. n.n. E, i i 1 J k j-|<~ ' n n IJ i n 1 1 "ik 2 2 J 1 i=1 j =1 k = 1 ' 9w '' 9w ' '5n m... ijk [ 9n = J I I I i i k i i j=1 k=1 I ' a j=1 a d m n... ijk m. + m. + m. i j k n = K n n. + n. + n. . i j k (3.24) Using these relations and the requirement of symmetry, the strain energy may be expressed as; 20 20 U I I 6 I \ i=1 j=1 k=1 m... -2 n... a.a.b. m.m. E, ' r\ ' ' J k i j m-1 n,.,-l) m . " ,-2 j k i]k~ a a d ( m . n . m . n ) 4 ' n |d£d + a.a.c, n.n. £ ::I j j k n + i j k + i i j k n . .(3.25) 33 Finally, taking the integration inside the summation, and recognizing Euler's Beta function to give f m n . n+1 / E n dEdn =c ja F(m, n) the r m + 1 . . *m+l\ m!n! -(-b) | j ( m + n + 2 ) , _ . • • • -(3.26) geometric stiffness matrix relative to the generalized coordinates may be determined as; f 6 k.. = h T <b, m.m.F(m... -2,n... ) + c. n.n.F(m... ,n... -2) ij k-i I J 'J 'J k i j ijk' ijk k 1 k k + d, (m.n. +m.n.)F(m... - 1,n... - 1) > k i j j i ijk ijk 1 Applying the same transformation and rotation matrices as . . . .(3.27) were used in the plate bending element gives the geometric stiffness matrix relative to the same eighteen nodal variables; [K ] c = [R] [T T 1 ] [k][T T ^[R] . . . . .(3.28) If the in-plane stresses are compressive, this matrix has the effect of making the plate more flexible, to the point that it may buckle. In the case of the rotating disk, the in-plane stresses are all tensile, and the matrix has the effect of making the plate stiffer. 34 3.3.4 Inertia Stiffness Because of the similarity between the inertia and geometric functionals, the inertia stiffness matrix by analogy to the geometric stiffness. may be easily derived Since the geometric stiffness matrix was developed by considering the strain energy in a local coordinate system, the inertia stiffness matrix will be developed in a local coordinate system as well. However, the strain energy is an invariant and does not depend on the coordinate system used, so the form of the integral is the same whether in global x - y coordinates or local E,-r\ coordinates. The inertia stiffness matrix does not exhibit this characteristic and hence a more rigorous be used than just exchanging transformation must £ and n for x and y. Considering an arbitrary point within the element domain (Figure 3 . 2 ) , the global coordinates may be expressed as a function of the local coordinates and constants, x = x Q + E, cos 9 - n sin 6 y = y Q + £ sin 6 + n cos 0 ....(3.29) from which the inverse relations may be obtained; £ = -x n = x_sin8-y Q cos 8 - y sin 9 + x cos 6 + y sin 0 Q n cos 9 - x sin 9 + y cos 8 . . . . . (3. 30) 35 By repeated use of the chain rule, the second derivatives of w in global coordinates may be expressed as functions of the second derivatives of w in local coordinates, so that; ( y 2 9w 2 9x' _ 8w 2 3 w - 2xy + x 3x9y 9y' 2 2 9w 2 2 (y cos 0 - x sin 0) 9E< 9 w 9 £9n 2 + 2(x sin 6 - y cos 6) (x cos 8 + y sin 0) 2 9 w 2 + (y sin 0 + x cos 0) 3n Using the Equations 3.29 . (3.31) 2 to express the x - y coordinates in terms of £ and r\ gives; f 2 y 29w 2 Y ~ 9x 2 _ 9w 29w 2xy + x 9y j 9x9y 2 9 w 2 (y cos 0 - x 0 Q sin 0 + ri) 3^ 2 9w + 2 ( x s i n 0 - y cos 0 - n ) ( y s i n 0 + x c o s 0 + £) 2 Q Q o Q 3E3n 3w ^ 3 £ 2 + ( y s i n 0 + x c o s 0 + £) n n T h u s , the inertia terms of the differential equation correctly in local E-ri coordinates. (3.32) are expressed 36 Considering the geometric stiffness terms in local coordinates; 2 9 w n 2 _ + 2oV._ 9£ ? ^ 9 w 2 , 3 w + 0_ 9£9n . =;— 8n n , . . 1 3 . 33J . . . .(3.34) . f 2 the following analogies can be made; 2 ( y c o s 8 - x s i n 6 + n) Q (x Q is Q sin 9 - y c o s 6 - n) ( y Q Q sin 6 + x c o s 6 + £) Q 2 ( y sin 8 + x cos 8 + £) Q is Q On the basis of this analogy, a is . a the inertia stiffness matrix is evaluated in exactly the same way as was the geometric stiffness matrix. Furthermore, these analogous stress terms are only quadratic functions of £ and n and hence there is no approximation in this derivation as there was in assuming that the stresses varied only quadratically. The analogous stresses may be expressed in a series form (where the overhead bar denotes inertia terms); \ 6 ^ ov ^ - I i=1 = m. bj K 1 6 _ m. n. J d. E n i=1 1 6 n. ' n ' a n = I i=1 m. n. c. K ' n 1 . . . . (3.35) 37 so that the inertia stiffness matrix relative to the generalized coordinates is 2 k.. = phQ ij ^ r- 6 T \b. m.m. F(m... -2, n... ) +c. n.n. F(m... , n... -2) ( k i j ijk ijk k i j ijk ijk + d. (m.n. + m.n.) F(m... - 1,n... - 1 ) V k i j j i ijk ' IJk J where m... = m. + m. + m. ijk i j k The coefficients b\, c., and dj 2 = 1 b r Q 2 c i d d (y = 1 2 s 0 i n 9 + x o c o . . .(3.36) — (y cos 6 - x sin 8) Q n... = m. + m. + n. ijk i j k are zero except for the following; — b and s 9 - = 3 2 ( y c o s 6 - x s i n 0) Q b Q & =1 - ' c 2 = 2 (y s i n 9 Q + x o c o s 9 ^ c 4 = 1 = ( x sin 6 - y c o s 6) (x c o s 6 + y sin 0) Q = Q Q ( x s i n 0 - y c o s 0) Q Q Q d^ = - ( x c o s 0 + y sin 0) d Q Q . The transformation and rotation 5 . . = 1 . .(3.37) matrices are applied as before to define the inertia stiffness matrix relative to the same eighteen nodal variables. Whereas the geometric stiffness matrix has the effect of making the structure stiffer (for tensile in-plane matrix makes the structure more flexible. stresses), the inertia stiffness Even though both the 38 geometric and inertia stiffness components are directly proportional to pft , 2 it is possible for the inertia stiffness to dominate all the stiffness terms and create a dynamic instability problem. This will be explored in more detail in the next chapter. 3.4 Convergence and Bounds The strength of any approximate method of solution is enhanced if it can be rigorously proven that in some limit, the approximate method will render the exact solution. Furthermore, the approximate method is made even more attractive if some asymptotic rate of convergence to the exact solution can be determined, and if the approximate solution itself represents some sort of bound on the exact solution. Finite element theory provides a method for assessing all of these characteristics. It can be shown that both the QST and plate bending elements will give answers that converge to the correct solution, and that the convergence rate with respect to the strain energy will be of the order of (1/N^) where N is the number of elements along one dimension. Furthermore, it can also be shown that the strain energy obtained by using these elements will represent a lower bound on the exact strain energy. The geometric stiffness matrix presents a slightly different case. It can be proven that this matrix will give solutions that converge to the exact solution, and that the strain energy obtained 39 r e p r e s e n t s a lower b o u n d on the e x a c t s t r a i n e n e r g y . b e c a u s e only c o n t i n u i t y of d e f l e c t i o n s c a n be shown that the asymptotic (and Furthermore, not slopes) convergence is r e q u i r e d , it rate with r e s p e c t to the g strain energy will b e of the o r d e r of v e r g e n c e p r o o f imposes stresses satisfy (1/N ). Unfortunately, the r a t h e r s t r i c t r e q u i r e m e n t that the the c o n in-plane the i n - p l a n e e q u i l i b r i u m e q u a t i o n s a n d b o u n d a r y d i t i o n s , a n d also be c o n t i n u o u s the i n - p l a n e s t r e s s e s across element b o u n d a r i e s . do not meet these r e q u i r e m e n t s , c a n n o t be p r o v e n a n d t h e r e a r e no b o u n d s T h e inertia stiffness for the e r r o r in the a s s u m e d with the geometric stiffness matrix general, convergence solutions. is r e p r e s e n t e d e x a c t l y displacement matrix, on the so In con- field. Thus, the c o n v e r g e n c e except by analogy rate s h o u l d be of p the o r d e r of (1/N ). As that of the linear b e n d i n g convergence this stiffness, it s h o u l d rate of the total o u t - o f - p l a n e In a p p l y i n g the a s s u m e d problem g e o m e t r y . faster convergence not a d v e r s e l y problem. the rate of c o n v e r g e n c e displacement f i e l d as b y Modelling of s t r a i g h t constant line s e g m e n t s . throughout ), the problem domain, is p r o p o r t i o n a l to the e r r o r in the a r e a . density straight It can where N is now the If the s t r a i n e n e r g y is the f a c t o r be number assumed then the e r r o r in s t r a i n Thus, not the modelling of the line segments i n t r o d u c e s an e r r o r in area r e p r e s e n t a t i o n . 2 (1/N problem is g o v e r n e d a c u r v e d edge with a s e r i e s of shown that this e r r o r d e c r e a s e s as than a f f e c t the s u c h h i g h p r e c i s i o n finite elements to the of a n a l y s i s of c i r c u l a r d i s k s , so much b y rate r e p r e s e n t s energy governing 40 convergence of the problems examined herein is the representation of the problem geometry, which leads to an expected convergence rate of the strain energy of the order of about (1/N 2 ). 41 CHAPTER IV SAMPLE PROBLEMS 4.1 Introduction The finite element formulation developed in the previous chapter is now used to investigate several problems. The geometric stiffness matrix and consistent mass matrix are applied to the problem of free vibrations of a square membrane to show that the predicted rate of convergence can occur under ideal conditions. no fundamental errors in the devised As a check that there are program, a rather unique problem is which makes use of all the out-of-plane stiffness matrices and for which an analytic solution may be estimated for comparison purposes. Free vibration results for rotating axisymmetric disks will be compared to solutions obtained by other methods, and the problem of a rotating disk subject to a space fixed transverse load is examined in some detail. Finally, several problems involving a non-axisymmetric disk will be considered. If the geometric stiffness matrix is used in the out-of-plane problem, some in-plane stress distribution is required. this is obtained by solving the in-plane problem Typically, using the QST finite elements, although in some cases the forced exact distribution is used. In all the problems, the same finite element grid is used for both the in-plane and out-of-plane problems. T h u s , it should be understood 42 that hereinafter, when the geometric stiffness matrix is used, the in-plane problem was also solved in order to obtain the in-plane stress distribution. In all the examples considered herein, only kinematic (forced) boundary conditions are enforced in the finite element solution, as required by the theorem of minimum potential energy. When such boundary conditions are to be enforced on a curved boundary, the conditions are approximated by using a local rotation of coordinates so that the nodal variables are relative to the normal and tangential directions of the boundary at that point. For the in-plane problem, the clamped boundary condition of the hub requires that both the u and v displacements and both of the tangential derivatives be zero. For the out-of-plane problem, this boundary condition constrains all of the variables except for the normal curvature variable. All computer calculations were done on an IBM 370/168 system with an AMDAHL CPU unit. Double precision was used throughout in order to minimize the effects of rounding errors and so that the high rate of convergence could be observed in some test examples. Eigen- problems were solved using a power iteration method which returned a specified number of the lowest eigenvalues and associated eigenvectors. The vibration modes were determined in each case by plotting the eigenvectors. 43 4. 2 4.2.1 Test Problems Geometric Stiffness Matrix As mentioned in the previous chapter, when using the geometric stiffness matrix, convergence to the exact solution can only be guaranteed if the in-plane stresses are represented exactly. In order to see the predicted rate of convergence, another restriction must be met. potential energy The theorem requires only deflection continuity, whereas the assumed displacement field is continuous with respect to both deflections and slopes. Hence, any exact solution used for comparison must be continuous in both displacements and slopes. The problem of free vibration of a rectangular membrane meets all of these requirements. As the problem boundaries for a rectangle can also be represented exactly, any error in the finite element solution may be attributed solely to the assumed displacement field. T h e n , the relative error in the eigenvalues should converge to zero at o the rate of 1/N , where N is the number of elements along a side. The mode shapes associated with the free vibrations of a rectangular membrane subject to uniform in-plane stresses are defined by sine waves in two directions. For a square membrane, the non-dimensional frequency coefficient is p a to 2 A = a n 2 = (m , 2 2, + n ) TT a 2 where m and n are the number of half waves in each direction, p is 44 the mass density, a is the dimension of the square, a is the in-plane stress and co is the natural frequency. Using the geometric stiffness matrix and consistent mass matrix, the eigenproblem was solved at three different grid sizes. Only the frequencies corresponding to the first four modes are presented, as these were the only modes which could be identified for the coarsest grid. For each grid and (2,1) size, the frequencies associated with modes were identical. (1,2) Symmetry and anti-symmetry conditions were applied when solving the problem at the finest grid size. The grids used and the non-dimensional frequencies obtained are shown in Figure 4.1 and Table 4.1, respectively. A log-log plot of relative error in the eigenvalues versus the number of elements along a side is shown in Figure 4.2. for reference. A line of slope 1/N is plotted The numerical results confirm the theoretical expectations, and in fact the convergence rate for the lowest eigenvalues was indeed 1/N 4.2.2 g above N = 2. Combined Geometric and Inertia Stiffness Matrices The accuracy and correctness of any newly formulated numerical solution technique remain suspect checked against known results. until the results have been In the case of the inertia stiffness matrix, this presents a difficult problem, as no exact solutions utilizing such a stiffness could be found. In this section a unique problem is examined which incorporates all of the out-of-plane stiffness effects and yet allows an approximate analytic solution. N = 4 F i g u r e 4.1 L a y o u t of F i n i t e Elements for Membrane TABLE Non-Dimensional Mode 4.1 F r e q u e n c y C o e f f i c i e n t s f o r Membrane N = 1 Vibration N = 2 N = 4 Vibration Exact (1,1) 19.75938400 19.73949275 19.73921000 19.73920881 (1,2)(2,1) 49.72110508 49.51355991 49.34866133 49.34802203 (2,2) 79.85243116 79.48476475 78.95934266 79.95683525 2 4 Number of Elements Per Side Figure 4.2 Vibrating Square Membrane, Relative Error in Eigenvalues 47 Consider a thin annular ring, rotating in its own plane at a constant angular velocity fi, as shown in Figure 4 . 3 . The radius R is very large compared to the cross-sectional dimensions of the ring, thickness h and width b. Equation ( 2 . 1 0 ) The x and y axes are space fixed so that is applicable. If the ring is free of any surface tractions, and is subject only to the radial in-plane body force prfi 2 where r is the radius and p is the mass density, then there exists an exact solution for the in-plane stresses based on plane stress elasticity, This solution may be degenerated so as to neglect any variation of the stress through the width of the ring, or a strength of materials approach may be used to render = Og pR 2 2 fi o = r a r ( Now consider a small element of the ring located on the x axis as shown in Figure 4 . 3 . N xy = a J i = 0 , and r6 N At this particular point, N 2 y y = 0 and x = R, so that Equation ( 2 . 1 0 ) 4 2 2 2 = a„h = phR fi . 8 2 2 x - o h = 0, r Furthermore, at this point becomes DV w - phR fi w + phR fi w yy yy - q or 4 DV w T h u s , in the limit that the rotating disk becomes a thin rotating annular ring, the governing differential equation approaches the static plate 48 y bh uniform ring width b thickness h N see Detail 111]N Detail Element of Ring on x Axis Showing Finite Element Grid Figure 4.3 Thin Rotating Ring T A B L E 4.2 Thin Rotating Ring Deflection Results Central Deflection Finite Element Bending Stiffness Only I Bending & Geometric Stiffness (inches) (inches) 1.3976 1.3958 0.6643 0.6630 Bending, Geometric & Inertia Stiffness I Beam Theory 1.4049 1.3958 49 e q u a t i o n f o r small elements of the r i n g . So if the p r o g r a m is used to a n a l y z e a small element of the r i n g , the s o l u t i o n s h o u l d b e c o m p a r able to the solution o b t a i n e d u s i n g only the b e n d i n g stiffness. S u c h a p r o b l e m was a n a l y z e d u s i n g a c o a r s e f i n i t e element g r i d of f o u r elements a r r a n g e d as shown of dimensions dimensions h : b : L : R was 1:10:100:10000, c a n c e l l a t i o n of t e r m s , so that the element S i n c e the example d e p e n d s on the it is important that all terms r e l a t i n g to the be of the same o r d e r of m a g n i t u d e , so that the c a n c e l l a t i o n will h a v e a s i g n i f i c a n t e f f e c t on the s o l u t i o n . to be simply s u p p o r t e d on the e d g e s y - ±^ the d e f l e c t i o n may b e a p p r o x i m a t e d b y governing T h e ratio were small c o m p a r e d to the r a d i u s , e n a b l i n g the a p p r o x i - mations to be r e a s o n a b l y a c c u r a t e . stiffness in F i g u r e 4.3. If the plate is assumed a n d f r e e on the o t h e r two, w = w^cos so that the d i f f e r e n t i a l e q u a t i o n becomes r _ ,TT,4 D (j-) O n this b a s i s , D ^= L - , ,-.22.71,2 phR Q (j-) + phR 2 2. TT , 2 Q, (j-) Try W 0 C O S L the c o e f f i c i e n t s were a d j u s t e d so that = phR ft 2 2 a n d f o r s i m p l i c i t y , a u n i f o r m l y d i s t r i b u t e d t r a n s v e r s e load was used. T h e r e s u l t s of this i n v e s t i g a t i o n a r e d i s p l a y e d in T a b l e 4.2. T h e d e f l e c t i o n s shown a r e the a v e r a g e s of the two nodes on the x-axis 50 for each case. Three cases were tested; one with only the bending stiffness, one using both the bending and geometric stiffnesses, and one using all of the out-of-plane stiffness terms. was actually done four times, The last case once on each of the positive and negative x and y axes, as a confirmation that all transformations were being done correctly. Any difference was negligible, and was likely due to a reversal of the orientation of the elements. The solutions based on one dimensional beam theory are shown for each case. The finite element solutions agree well with the beam theory results and confirm the expected phenomenon regarding the cancellation of geometric and inertia stiffness effects. 4.3 Vibrations of Axisymmetric Disks The equation governing the transverse vibrations of rotating disks was given in Chapter 2, relative to a body fixed coordinate system. If the center of rotation and the coordinate origin coincide, and if the problem is axisymmetric, it can be shown that Equation (2.5), when transformed to polar coordinates, is exactly the same equation solved by Lamb and Southwell [1] and Eversman and Dodson [3]. 51 In this section, the finite element formulation presented in the previous chapter is applied to the problem of free transverse vibrations of rotating axisymmetric disks, and the results compared with those obtained in References [1] and [3]. a means of assessing the convergence in which the in-plane stresses are This provides characteristics not for cases represented exactly, so that convergence cannot be proven and bounds on the solution do not exist. 4.3.1 Vibrations of a Rotating Membrane Lamb and Southwell [1] obtained the solution natural frequencies of transverse vibration of a rotating case in which the flexural rigidity is negligible. for disk the for the When there is no central hub, the resulting formula for the natural frequencies depends only on Poisson's ratio, v , the rotation speed fi, and the number of nodal diameters and nodal circles which characterize the associated mode shape. by The non-dimensional frequency parameter X is given 52 to x 2 (s + 2n) (s + 2n + 2) = -n = ST { I 3+v^ 8 J s — 2 fl+3v ' { 8 j w h e r e n is the n u m b e r of n o d a l c i r c l e s , s is the n u m b e r of n o d a l d i a m e t e r s , and is the n a t u r a l frequency. T h e g e o m e t r i c s t i f f n e s s m a t r i x a n d t h e c o n s i s t e n t mass m a t r i x were u s e d to o b t a i n the n a t u r a l f r e q u e n c i e s a n d v i b r a t i o n modes t h e f i n i t e element g r i d s s h o w n in F i g u r e 4.4. were generated u s i n g The in-plane the q u a d r a t i c s t r e s s t r i a n g l e . The using stresses out-of-plane v i b r a t i o n p r o b l e m was s o l v e d twice f o r each c a s e , one f o r c i n g a n t i symmetry along both coordinates axes ( A - A ) a n d the o t h e r f o r c i n g a n t i - s y m m e t r y a l o n g one a x i s a n d s y m m e t r y a l o n g the o t h e r a x i s (A-S). T h i s a l l o w e d all modes to be g e n e r a t e d e x c e p t f o r those w i t h no nodal diameters. In the s o l u t i o n o b t a i n e d b y Lamb a n d S o u t h w e l l , these modes i n v o l v e d a d i s p l a c e m e n t of the c e n t r e of the d i s k a n d so a r e not i n c l u d e d f o r c o m p a r i s o n The expressions here. f o r the i n - p l a n e s t r e s s e s w h i c h t h e r e is no c e n t r a l h u b a r e g i v e n stress f o r the case in in C h a p t e r 3. The shear is i d e n t i c a l l y z e r o , a n d the r a d i a l a n d c i r c u m f e r e n t i a l s t r e s s e s a r e o n l y q u a d r a t i c f u n c t i o n s of the r a d i u s . S i n c e the Q u a d r a t i c S t r e s s T r i a n g l e is c a p a b l e of modelling e x a c t l y a q u a d r a t i c stress d i s t r i b u t i o n , a n y e r r o r may be a t t r i b u t e d to the d i f f i c u l t y of m o d e l l i n g a curved boundary with s t r a i g h t line segments. radial and circumferential stresses, -^-^—, pft b 2 The non-dimensional f o r each of the t h r e e 2 f i n i t e element g r i d s a n d the e x a c t v a l u e s , a r e s h o w n in F i g u r e s 4.5(a) N = 2 N = 1 Figure 4.4 Finite Element Grids for Disk Without Hub N F i g u r e 4. 5(a) N o n - D i m e n s i o n a l I n - P l a n e Radial Stress Distributions K = 0. 0 Figure 4.5(b) Non-Dimensional In-Plane Circumferential Stress Distributions K =0.0 56 and 4.5(b), respectively. With the coarse g r i d , the poor representa- tion of the problem boundary significantly affects both of the stresses throughout the entire domain. However, with the finest g r i d , the radial stresses are represented nearly exactly, and any error in the circumferential stress distribution is confined mainly to the boundary region. The non-dimensional vibration frequencies obtained for each finite element g r i d , together with the exact values, for the first modes are listed in Table 4.3. ratio of 0.3. vectors. twenty-t These values are for a Poisson's The modes were determined by inspection of the eigen- At N = 1, only the first fourteen modes could be identified from the results for both the A - A and A-S cases, but with N = 2 and N = 3 all twenty-two modes could be identified. For N = 1 and N = 2 only the first few modes are ordered correctly, but at N = 3, the finite element solution produced the twenty-two modes in correct order. This is a distinct advantage over the analytic or semi-analytic methods, as these methods usually involve finding the frequency corresponding to a particular mode. T h u s , to find the fundamental frequency or a certain number of the lowest frequencies involves inspection of many mode shapes This advantage becomes more apparent when solving more difficult problems. It may be observed from Table 4.3 that convergence is not monotonic for all modes. However, if the absolute value of the relative error in each case is considered, convergence is monotonic, except for the mode (0,8), which was one of the highest modes obtained for the TABLE Non-Dimensional N a t u r a l F r e q u e n c i e s of a Rotating A From F i n i t e N = 1 Mode* 4.3 Elements Membrane A N = 2 N = 3 Exact (0,1) 0.9991645 0.9997996 0.9999322 1. 00 (0,2) 2.3267 2.3474 2.3493 2. 35 (0,3) 3.9319 4.0411 4.0478 4.05 (1,1) 5.8591 5.9481 5.9499 5. 95 (0,4) 6.8488 6.0758 6.0946 6.10 (0,5) 8.0989 8.4378 8.4882 8. 50 (1,2) 8.0544 8.9325 8.9476 8.95 (0,6) 11.514 11.093 11.229 11.25 (1,3) 12.412 12.246 12.296 12. 30 (2,1) 11.853 14.147 14.202 14. 20 (0,7) 16.227 13.975 14.308 14. 35 (1,4) 19.286 15.854 15. 991 16. 00 (0,8) 18.777 19.182 17.731 17.80 (2,2) 17.150 18.644 18.839 18. 85 (1,5) 19.630 20.021 20.05 (0,9) 21.385 21.472 21.60 (2,3) 23.504 23.823 23.85 (1,6) 22.422 24.361 24.45 (0,10) 27.414 25.508 25. 75 (3,1) 25.388 25.712 25. 75 (2,4) 28.725 29. 190 29.20 (1,7) 29.373 29.121 29. 20 Mode (n,s) n - n u m b e r of nodal c i r c l e s s - n u m b e r of nodal diameters 58 grid N = 1. The rates of convergence may be determined from the log-log plots of relative error versus N, as shown in Figures 4.6(a) -2 and 4.6(b). Generally, the rates are better than N , which corres- ponds to the rate of convergence of error of area representation of the problem domain. This probably reflects the fact that convergence is due to two effects in this example, the modelling of the problem geometry and the modelling of the in-plane stresses. At the grid size N = 3, the relative error for each of the twenty-two modes is less than 1% and less than 0.1% for most of the lower modes. The numerical results confirm that convergence will in fact occur, although the bounds on the solution are lost. Furthermore, this example indicates that good results may be obtained even with a fairly coarse g r i d . 4.3.2 Vibrations of a Rotating Disk - Without Hub Lamb and Southwell h ] have also presented bounds on the natural frequencies of vibration for certain modes for the case in which both the flexural rigidity and the effects of rotation are significant, when there is no central hub. The upper bounds are derived only for the modes with no nodal circles, and involve finding the roots of a quadratic equation which has rather unwieldy coefficients. The lower bounds are obtained by calculating the square of the natural frequency for each of the limiting cases (no flexural rigidity, and no rotational effects) and adding them to obtain the square of the lower limit of the natural frequency for the combined effect. However, this method 59 Grid Figure 4.6(a) Size Vibrating Rotating Membrane, Relative Error in Eigenvalues 60 Grid Size Figure 4.6(b) Vibrating Rotating Membrane, Relative Error in Eigenvalues 61 requires knowing the non-dimensional frequency coefficient for a vibrating circular plate for each mode of interest. The finite element grids shown in Figure 4.4 were used to determine the natural frequencies of vibration for a rotating disk when both rotational and flexural effects are present. The disk analyzed had a radius of sixty inches, a thickness of one inch and was rotating at 1, 000 RPM. The material properties corresponded to steel, with a modulus of elasticity of 30, 000 KSI, a Poisson's ratio of 0.3 and a weight per unit volume of 490 P C F . As for the membrane vibration problem, each out-of-plane problem was solved twice in order to force different symmetry requirements. The values of the natural frequencies for each mode are listed in Table 4.4, along with the upper and lower bounds where applicable. In contrast to the membrane problem, the order of the modes did not vary when the finite element grids were refined. The values of the natural frequencies in Table 4.4 are not dimensionless as they were in Table 4.3. This is because when both flexural and rotational effects are operative, a non-dimensional f r e quency is difficult to define and does not make the results any less cumbersome. However, a useful non-dimensional parameter is avail- able for this case which gives an indication of the relative importance of these two effects. a = The non-dimensional parameter 2 2Eh 2 3 ( 3 + v ) ( 1-v ) p f i 2 b " 62 TABLE Natural Vibration Frequencies 4.4 (Hz) From F i n i t e of a R o t a t i n g Disk Elements Upper Bound Mode* Lower Bound (0,1) 16.67 16.67 16.67 16.67 16.67 (0,2) 29.41 30.17 29.61 29. 51 29.44 (0,3) 47. 55 49.24 48.30 47. 96 47. 85 (1,1) 68.81 71.36 70.71 70. 07 - (0,4) - 82.08 73. 55 72.86 72. 93 (0,5) - (1,2) 106.6 - (0,6) (1,3) 155.5 (2,1) 168.4 N = 1 N = 2 N = 3 104.8 111.2 105.3 104. 1 117.5 111.2 110. 0 155.7 143. 4 141.7 160.0 160.2 158. 2 - 180. 5 177.7 •188.6 - 143.4 (0,7) - - - 185. 3 (1,4) - - - 214.2 (0,8) - - - 235. 1 - - 244.9 - (2,2) 240.1 (1,5) - - - 277.6 (0,9) - - - 290.8 240.4 298.8 (2,3) 303.7 - - 320. 7 - (3,1) 335.2 - - 340.6 - - (1,6) - - - 348.6 (0,10) - - - 352. 5 Mode ( n , s ) n - n u m b e r of nodal c i r c l e s 363.7 s - n u m b e r of nodal diameters 63 v a r i e s from z e r o , f o r the case of a p u r e membrane in which only effects are present, rotational to i n f i n i t y , f o r the case of a p u r e plate in which only f l e x u r a l e f f e c t s a r e p r e s e n t . which r e p r e s e n t s a f a i r l y s t i f f If the d i s k In this p a r t i c u l a r e x a m p l e , a = 0.0639 disk. has no c e n t r a l h u b , rotation in a plane s l i g h t l y then the mode (0,1) represents i n c l i n e d to the r e f e r e n c e p l a n e , so that the natural f r e q u e n c y associated with this mode is the rotation f r e q u e n c y . T h i s c h a r a c t e r i s t i c was also e x h i b i t e d in the case of the membrane solved p r e v i o u s l y , in which X - 1 f o r this mode. a r e not available f o r the o t h e r modes, conclusions from these r e s u l t s . As exact solutions it is d i f f i c u l t to draw a n y H o w e v e r , f o r the modes problem which hard have b o u n d e d f r e q u e n c i e s , the finite element solutions a p p e a r to c o n v e r g e values within the b o u n d s , a g a i n i n d i c a t i n g that a reasonable to d e g r e e of a c c u r a c y may be o b t a i n e d e v e n with a f a i r l y c o a r s e f i n i t e element grid. 4.3.3 Vibrations of a Rotating D i s k - With T h e problem of f r e e t r a n s v e r s e v i b r a t i o n s of a r o t a t i n g is now examined f o r the case when the d i s k and Dodson such a disk, [3 ] frequency has a c e n t r a l h u b . disk Eversman h a v e p r e s e n t e d r e s u l t s f o r the n a t u r a l f r e q u e n c i e s of for various f o r all combinations Hub non-dimensional h u b ratios a n d stiffnesses, of 0, 1, o r 2 nodal diameters or c i r c l e s . is related to the o t h e r f a c t o r s t h r o u g h a frequency parameter, defined by The natural non-dimensional 64 1 8 (3+v)a f " V2 + I « J ' 1+3v ' 8 which is determined graphically. The combined linear bending and geometric stiffness matrix and the consistent mass matrix were used to approximate the natural frequencies and vibration modes for the finite element grid shown in Figure 4 . 7 . The non-dimensional hub radius is 0 . 4 4 . To model the complete clamp of the hub, all kinematic variables were constrained on the inner radius boundary as described in the introduction of this chapter. The out-of-plane problem was solved only once, forcing symmetry along the diameter, generated. which enabled all possible modes to be A non-dimensional stiffness of a - 0.001 and a Poisson's ratio of 0 . 3 were used. The values of the fourth root of the non-dimensional f r e quency parameter for each mode as obtained from the finite element analysis are listed in Table 4 . 5 . The values of the parameter corresponding to the modes considered in Reference [3 ] are also listed. Only the first twenty frequencies were calculated, and none of the modes characterized by two nodal circles were present. In the few cases where results are available for comparison, the finite element method shows good agreement, with a maximum error of about 2%. Hence, good engineering accuracy is again obtained with a fairly coarse g r i d . Fiqure 4.7 Finite Element Grid for Disk With Hub T A B L E 4. 5 Non-Dimensional Frequency Parameters of a Rotating Disk; K = 0.44, a = 0.001, v = 0.3 / 7 Mode (0,0) 8. 24 8.2 (0,1) 8.67 8.6 (0,2) 9.68 9.7 (0,3) 10.87 (0,4) 12.07 (0,5) 13.24 (0,6) 14. 36 (0,7) 15.43 (0,8) 16. 38 (1,0) 13. 76 13. 5 (1,1) 13. 90 13. 7 (1,2) 14.32 14.2 (1,3) 14.95 (0,9) 17. 50 (1,4) 15.68 (0,10) 18.47 (1,5) 16.60 (0,11) 19.49 (1,6) 17. 59 (0,12) 20.42 (2,0) not observed 18.8 (2,1) not observed 18.9 (2,2) not observed 19.0 Mode (n,s) n - number of nodal circles s - number of nodal diameters 1From Finite Elements "From Reference 3 (2) 67 This example also demonstrates an advantage of the finite element method that was mentioned previously. Finding the fundamen- tal frequency or a certain number of the lowest frequencies with the finite element method does not involve considering many different mode shapes, as in the analytic or semi-analytic techniques. Rather, the lowest frequencies and corresponding modes are approximated by the finite element solution so that when a reasonably fine element grid is used, the order of the lowest modes is obtained fairly accurately. This ability is displayed in Table 4.5, in which case the semi-analytic method could not predict which modes (if any) occurred between modes (0,2) and (1,0) and modes (1,2) and (2,0). 4.4 Stiffness Analysis of an Axisymmetric Disk The problem of a rotating disk subject to a space-fixed, constant, transverse point load is now considered. This problem is analyzed in a space-fixed, rather than a body-fixed, coordinate system and is thus governed by the differential equation derived in Chapter 2 (Equation 2.10), so that the inertia stiffness matrix now contributes to the total problem stiffness. It should be emphasized that this analysis is valid only for cases in which the disk geometry and the in-plane problem are axisymmetric, although the out-of-plane problem need not be axisymmetric. For all the out-of-plane problems c o n - sidered in this section, the non-dimensional hub radius, K , is 0.2 and the non-dimensional load radius is 0.8. 68 The finite element grids used in this analysis are shown in Figure 4.8. The out-of-plane solution was assumed to be symmetric about the load diameter, so only half the disk was a n a l y z e d . 1 The grid was altered near the load point so that the load could be applied at a node. For the out-of-plane problem, a full clamp boundary condition was applied at the hub boundary. The in-plane stresses for each grid were generated by two methods. In the analysis of Reference [ 9 ] , the in-plane stresses corresponding to a complete disk were used for simplicity. order to remain consistent with that particular analysis, T h u s , in the in-plane stresses were generated by forcing a quadratic function to match the exact values, for a complete disk (K = 0.0), at the nodes and mid-sides of each element, as described in Section 3.2.3. Since the exact stress distribution is only a quadratic function, a nearly exact distribution of stresses should be generated throughout the problem domain. The in-plane stresses were also derived by applying a clamped boundary condition on the hub radius, and using the QST finite element with the consistent body force load vector to generate the stress distribution. The exact stress distribution involves a 2 1/r term, which makes it more difficult for the finite element to generate a reasonable stress distribution. Plots of the non-dimensional radial and circumferential stresses for the two grids, and the exact solution, are shown in Figures 4.9 (a) and (b) . The stresses are ^ h i s assumption is verified further on in this work. 69 Figure 4.8 Finite Element Grids for Stiffness Analysis 70 Exact Figure 4.9(a) Non-Dimensional In-Plane Radial Stress Distributions K =0.2 71 Exact Figure 4.9(b) Non-Dimensional In-Plane Circumferential Stress Distributions K = 0.2 72 non-dimensionalized by multiplying by P pfi b 2 . The maximum values 2 are 4.92 and 2.45 for the radial and circumferential stresses, respectively. These plots indicate that the finite element method is providing a reasonable approximation to the stress distribution, particularly in the region nearer 2 the outer edge, in which the 1/r term does not dominate. Although the results for the circumferential stress do not appear as good as those for the radial stress, along any radial line there is a good overall approximation to the stress distribution. An attempt was made to reproduce the results of Reference [9], in which the response of a rotating disk to a transverse load is obtained for various values of a . This out-of-plane problem was solved using each of the two grids and each of the two methods to generate the in-plane stress distributions, in order to observe the effects of the respective approximations. Since the computer program uses dimensional, rather than non-dimensional, parameters, the different values of a were obtained by varying the rotation speed. As the particular technique first used to solve the system of linear equations required that all the diagonal terms be positive, no solutions could be obtained past a certain value of a . This particular value of a , for each of the cases considered, was determined by extrapolating from a plot of the determinant of the stiffness matrix versus a, to determine when the determinant became zero. T h u s , it was found that, for the particular non-dimensional hub radius of 0.2, there exists a fundamental critical speed fi r R or non-dimensional stiffness, a r D , 73 at which the rotating disk has no stability in the transverse direction, and a standing (in space) wave appears. The values of this fundamental critical stiffness are shown in Table 4 . 6 . There is a fairly large variation in the values, which indicates that both of the approximations have significant effects. The deflected shape corresponding to a value of a just slightly greater than c t ^ , ^ is shown in Figure 4 . 1 0 . four cases This deflected shape, which was the same in each of the considered, is similar to the vibration mode The c o n - (0,2). tour intervals in the deflection plot represent equal proportions of the maximum deflection. T h u s , in Figure 4 . 1 0 , which has ten contour intervals, the contour lines adjacent to the nodal lines represent an absolute value of the deflection of 1 0 % of the maximum. This critical speed instability has been examined by Mote [ 4 ] and may also be related to the work done by Lamb and Southwell [ 1 ] . ln the latter work, it is noted that for any mode, the associated wave velocity in space, to , is given by t o = ft ± is the natural frequency corresponding fixed coordinate system, ft is the , s to the rotation mode, speed, number of nodal diameters associated with the mode. that there may be some combination of ft, to zero. where in to n a body- and s is the It is obvious and s such that a j g is This would mean that one of the natural vibration frequencies is associated with a mode shape that is standing in space. Mote has discussed this phenomenon in terms of stability of a rotating disk, and 2 defines the criterion of instability as ( t o ^ ± sft) 2 = g , harmonic excitation frequency of the transverse load. where 3 is the If the transverse load is constant and stationary in space, 6 = 0 . The Figure 4.10 Typical Deflected Shape Near Fundamental Critical Speed Table 4.6 Critical Stiffness from Finite Element Solution K = 0.2 Grid In-Plane Stresses N = 2 Forced Exact ( K = 0.0) N = 2 QST Generated (K N = 3 Forced Exact (K N = 3 QST Generated (K = 0.2) a = 0.2) = 0.0) CR 0.0707 [+00.0002 0.0814 ± 0.0002 0.0805 ± 0.0002 0.0881 ± 0.0002 CO criterion for instability may be stated as ft ± -j- - 0, which, from the previous reference, is when the natural vibration mode becomes fixed in space. The fundamental critical speed is then defined by ^ CR ~ m m n s = 1, 2, 3, s In order to ensure that this is the phenomenon actually being observed in the finite element problem, the free vibration problem was solved to determine c o grid N = 2 at a rotation speed corresponding to n ft The was modified slightly so that it was symmetric, and the free vibration problem was solved as in Section 4.2. As the vibration proble is not easily described by non-dimensional parameters, dimensional quantities are used. analysis, The material constants of steel were used in the so that the modulus of elasticity was 30,000 KSI, Poisson's ratio was 0.25 and the weight per unit volume was 490 PCF. The outer radius of the disk was 10 inches, the hub radius was 2 inches and the thickness was 0.1 inch. With these material constants and disk dimen- sions, the critical speed for the N = 2 grid with the QST generated i n plane stresses was 332 radians per second. The vibration modes and natural frequencies for this particular disk rotating at ft ^ in Table 4.7. (0,2) are listed The values of ft, co and s corresponding to the mode satisfy Mote's criterion of instability for 3 = 0, and, according to the analysis of Lamb and Southwell, result in a wave velocity, in space, of zero. TABLE Natural 4.7 F r e q u e n c i e s of a D i s k R o t a t i n g Mode 1 di 2 n (0,0) 461 (0,1) 493 (0,2) 664 (0,3) 1044 (0,4) 1603 (1,0) 2265 (0,5) 2322 (1,1) 2433 (1,2) 2938 (0,6) 3185 (1,3) 3776 (0,7) 4189 (1,4) 4878 (0,8) 5268 (1,5) 6318 (2,0) 6369 (0,9) 6627 (2,1) 6646 (2,2) 7414 (1,6) 7831 Mode ( n , s ) n - n u m b e r of nodal c i r c l e s s - n u m b e r of nodal 2 For fi @ fi = 332 r a d / s e c diameters 77 These results confirm that the instability observed in a spacefixed coordinate system is directly related to the free vibration problem in a body-fixed coordinate system. It may also be concluded that the critical value of the non-dimensional stiffness is a function only of the non-dimensional hub radius, since these are the only two non-dimensional parameters which define the instability problem. T h u s , given a non- dimensional hub radius, there exists a single value of the critical nondimensional stiffness, so that if the disk material and geometry are known, the fundamental critical speed may be determined directly. The results of Tobias and Arnold [16] indicate that a phase shift occurs at the resonant condition, so that the maximum deflection is not directly under the load. In order to ensure that the forced symmetry conditions in the analysis presented here have not precluded this phenomenon, an analysis of an entire disk was done. The full grid equivalent of Figure 4.8, N = 2, with the stresses forced exact (K = 0.0) was run at exactly the same value of 0, as used previously for the half g r i d , which was just slightly less than the fundamental c r i t ical value. The results were identical, to seven significant with those of the half g r i d , thus confirming the symmetry figures, assumption. Hence the present analysis did not reproduce the unsymmetric results reported in Reference [16]. These results were obtained by using a Choleski type decomposition of the total stiffness matrix in order to solve the transverse deflection problem. the diagonal. This type of routine fails if negative terms occur on T h u s , such a solution technique will not work at 78 stiffnesses past the critical stiffness, a „ However, from either a . D mathematical or a physical point of view, there must be a finite response at stiffnesses smaller than c t ^ . In order to obtain this response, a more general solution technique must be used. The response, in space fixed coordinates, to a stationary transverse point load was thus determined for values of a of 0.1, 0.05, 0.01 , 0.005 and 0.001, using the N = 3 grid. A non-dimensional hub radius of 0.2 was used for the out-of-plane problem, and the in-plane stresses were generated by forcing the exact values, based on a non-dimensional hub radius of 0.0. Plots of the deflection contours are shown in Figures 4.11(a) and 4.11(b). at a = 0.1, All of the plots, except display a single nodal diameter at approximately 90° from the load diameter. This characteristic is also exhibited for values of a significantly greater than a £R ( -9-' e a = would suggest that these particular values of n o t shown), which a are not very close to any of the higher critical stiffnesses. The location of the maximum deflections have been indicated on each of the plots in Figure 4.11(a). The non-dimensional maximum deflection, w, is defined by (3 + v ) b h 2 w = 2 (pQ ) w 8p where w is the dimensional deflection and p is the magnitude of the point load. The values of the non-dimensional maximums are shown below along with those from Reference [9]. 78a a Present Analysis Reference 9 4.40 3.51 -3.06 ±3.25 0.66 0.1 0.05 0.01 0.005 -4.01 -6.21 Neither the magnitudes or the deflected shapes of this analysis compare well with those of Reference [9]. No explanation for this could be found, although it is possible that either one or both solution techniques may be suffering from convergence problems. The results obtained by Mote [4] typically show a fairly narrow range, just above the fundamental critical speed, in which several of the other modes also become stationary in space, leading to an instability. Such an instability is characterized by a zero determinant of the stiffness matrix. Since the determinant at both a=0.1 and a=0.05 was positive, there is either no value of a in this range for which the determinant is zero,or an even number of critical values of a in this range. As the fundamental critical value lies within this range, there must be at least one other critical value, and possibly more. displayed here. T h u s , the trend in Mote's results is also It should be noted that the critical values of a could be determined explicitly by solving the eigenproblem Ax = ABx, in which A represents the linear bending stiffness matrix and B combined geometric and inertia stiffness matrix. values of a could then be related to the represents the The critical eigenvalues, X, and the eigenvectors, x, would define the mode shape (fixed in space) which was excited at each corresponding value of X . The results shown in Figure 4.11 do not show many of the dominant features of the plots in Reference [9], (with the complete 79 F i g u r e 4.11(a) D e f l e c t i o n C o n t o u r s , G r i d N = 3, Stresses Forced Exact N = 3 Stresses Forced Exact Figure 4.11(b) N = 2 Stresses Forced Exact Deflection Contours a = 0.001 81 disk in-plane stress distribution), which represent solutions to exactly the same problem considered here different technique. except that it is solved using a For this reason, the effects of the grid refine- ment and the in-plane stress distribution were examined for a = 0.001. The stiffness problem was solved at this value of a for the two different grids, using both the QST generated and forced exact inplane stresses based on non-dimensional hub radii of 0.2 and 0.0 respectively. The plots of these results, shown in Figure 4.11 (b) , illustrate the dramatic change in response obtained by refining the g r i d . The change due to the different in-plane stress distributions is not quite as obvious, but is still significant. T h u s , even though the grid N = 3 provides a close representation of the problem boundaries, and allows a reasonable approximation of the in-plane stress distribution, the significant change in response from N - 2 to N = 3 is disturbing, and raises doubts about the validity of the other results. However, this particularly low value of a corresponds to a very flexible disk ( e . g . , a ten-inch diameter steel disk rotating at about 30, 000 RPM) which the flexural stiffness is almost negligible. in Since the stiffness is then so dependent on the in-plane stresses, it would be expected that the response depends significantly on any approximations for these stresses. Furthermore, since the (destabilizing) used inertia stiffness matrix will have a dominant effect on the stiffness at very low values of a , it is reasonable to expect the response to be very unstable, and hence be very sensitive to any changes in the approximations used. 82 4. 5 Analysis of a Non-Axisymmetric Disk As the strength of the particular method formulated here lies in its ability to accept non-axisymmetric problem boundaries, some examples of such problems are now considered. The geometry of the disk modelled in this section, shown in Figure 4.12, is based on a prototype of a steel saw blade. The ten inch outer radius used in the model neglects the effect of the cutting teeth. In this section, the response of the static disk to a point load is compared to some experimental results, the free vibration problem of the non-rotating disk is examined, and the membrane stresses due to rotation are also considered. 4.5.1 Static Deflection Response, Using the finite element grid shown in Figure 4.13 with the linear bending stiffness element, the response of a static disk to a point load on its perimeter was analyzed. Experimental results were obtained for this particular example in which the load was applied at 15° intervals around the perimeter, starting at a point 5° off the centerline of the slot. finite element g r i d . This accounts for the lack of symmetry in the The load used in the experiment applied at a radius of 9.55 inches. was 26.1 pounds, For the finite element solution, the load was applied at the nodes on the outer radius. load points are indicated on the grid in Figure 4.13. The twelve Figure 4.12 Non-Axisymmetric Disk Geometry 6 Figure 4.13 Finite Element Grid for Non-Axisymmetric Disk 84 The results obtained from the finite element analysis and from the experiments are plotted in Figure 4.14, for the load at nodes 1, 3, 5, 8, 10, and 12. Although both experimental and finite element results were obtained for the load at all 12 nodes, the other plots show no unusual or unexpected trends or characteristics so are not presented. The finite element solution consistently underestimated the deflection response, although the deflected shape was predicted well. Because only the linear bending stiffness matrix was used in this analysis (for which convergence may be proven and bounds on the strain energy exist) it was expected that the solution would be stiff, but a higher order element such as used here should give better results than this. The discrepancy is likely due to the inability to obtain a truly clamped boundary condition when doing the experiment, whereas the clamped boundary condition is accurately modelled in the finite element solution. 4.5.2 Free Vibration Response The free vibration response (with no rotation) of this particular non-axisymmetric disk was analyzed using the finite element grid shown in Figure 4.15 (a). The grid is refined slightly near the tip of the slot to allow for the finite slot width. Four separate analyses were required in order to obtain all possible combinations of symmetry and anti-symmetry about each axis. F i g u r e 4. 14(a) Static Deflection of a N o n - A x i s y m m e t r i c Disk CO Figure 4.14(b) S t a t i c Deflection of a N o n - A x i s y m m e t r i c Disk T 2 3 4 5 6 7 8 9 10 II 12 Position 00 F i g u r e 4.14(c) S t a t i c Deflection of a N o n - A x i s y m m e t r i c Disk ^ (a) Figure 4.15 Finite Element Grid for Non-Axisymmetric Disk 89 The mode shapes and corresponding natural frequencies are listed in Table 4 . 8 , analysis. in the order in which they occurred in the The interesting feature of the vibration mode shapes is the appearance of a ^ nodal diameter along the slot diameter. This diameter has been labelled as ^ since it does not actually change the mode shape, but merely causes a phase shift between the two halves of the disk. The modes (0,3), (0,3j) and (0,4) are reproduced in Figure 4.16 to illustrate this feature. Although this effect occurs for all the modes listed in Table 4.8, it is interesting to note that as the number of nodal diameters increases, the difference in the natural frequencies for the whole mode and the half mode becomes negligible. Another feature that is illustrated by Figure 4.16 is that the maximum deflection occurs at the corner where the slot intersects the outer radius. For the particular modes shown here, this maximum deflection is at least 30% greater than the other local maximums. This feature was exhibited by all the modes, but the difference between the overall maximum and the local maximums was seldom this large. 4.5.3 In-Plane Stresses As there is no exact solution for the distribution of the i n - plane stresses in the non-axisymmetric problem, even the finite element solution to this problem is of interest. 8 - pQ b 2 a 2 The non-dimensional distributions are shown in Figure 4.17. stress These distributions - were obtained using the quadratic stress triangles with the grid shown 90 T A B L E 4.8 Natural Frequencies of Non-Axisymmetric Disk Mode * co n (Hz) Mode * co (Hz) n (0,0) 159 (1,2) 1186 (0,*) 160 (1,2*) 1261 (0,1) 163 (1,3) 1363 (0,1*) 166 (0,8) 1374 (0,2) 194 (0,8*) 1374 (0,2i) 206 (1,3*) 1443 (0,3) 267 (1,4) 1 597 (0,3*) 286 (1,4*) 1677 (0,4) 389 (0,9) 1721 (0,4*) 410 (0,9*) 1723 (0,5) 563 (1,5) 1915 (0,5|) 581 (1,5*) 1986 (0,6) 784 (0,10) 2111 (0,6*) 800 (0,10*) 2112 (0,7) 1030 (1,6) 2302 (1,0) 1063 (1,6*) 2363 (0,7*) 1063 (0,11*) 2541 (1,*) 1087 (0,11) 2541 (1,1) 1115 (1,7) 2732 (1,U) 1150 (1,7*) 2777 Mode (n,s) n - number of nodal circles s - number of nodal diameters 91 Mode (0,4) 389 Hz Figure 4.16 Typical Vibration Modes of a Non-Axisymmetric Disk 92 Shear Figure 4.17 Stress Non-Dimensional In-Plane Stress Distributions of a Non-Axisymmetric Disk 93 in Figure 4.15 (b) . Similar results were obtained using the grid shown in Figure 4.15 (a). The maximum values for the axisymmetric problem would be about 4.2, 1.6 and zero for the radial, circumferential and shear stresses, respectively. As was expected, the results indicated a stress concentration near the tip of the slot, so for clarity, the contours were omitted in this area. Except around the slot tip, there is no significant change in the radial stress distribution or the shear stress distribution, which is generally small, although not zero. stress, The circumferential though, shows quite a significant overall reduction in stress as a result of introducing the slot, and in fact, is reduced to zero, or nearly zero, around the whole perimeter, quite unlike the a x i symmetric distribution. This reduction in the circumferential stress is the motivation for putting the slot in the disk, since it will also reduce these stresses due to thermal effects. The stress distribution along the slot radius, from the hub to the tip of the slot, is plotted in Figure 4.18. This particular plot was obtained from the results using the grid in Figure 4.15 (a). As the tip of the slot is a free surface, it was expected that the radial stress and shear stress would approach zero at this point. To see if a more realistic stress distribution could be obtained, the grid in Figure 4.15 (b) was also used, including double-noding at the slot tip. The results represented no significant over those shown in Figure 4.18. improvement The problem lies in the finite Figure 4.18 Non-Dimensional Stresses Along Slot Diameter 95 element model of the prototype. Modelling the rounded slot tip with a few straight line segments introduces a stress singularity at the tip, which the finite element method is unable to accept, at least as formulated here. Further examination of this problem is beyond the scope of this work, and would likely have little, if any, effect on the overall problem. 96 CHAPTER V CONCLUSIONS A finite element method to determine the natural vibration modes and frequencies of rotating disks has been presented. The transverse vibration frequencies depend on both flexural and geometric stiffness. The finite element method has been applied to generate an approximate in-plane stress distribution, which is incorporated in the out-of-plane problem, which in turn is also solved using the finite element technique. Comparison with exact and approximate results obtained from other methods shows that, even though convergence cannot be proven and bounds on the solution do not exist, the results obtained from this finite element formulation do converge to the correct results. Furthermore, the finite element solutions indicate that good engineering accuracy may be obtained even with a fairly coarse grid of elements. As the finite elements used in this analysis were a triangular shape, no restrictions on the geometry of the problem had to be made. T h u s , the formulation presented herein, which applies a consistent load vector to generate the in-plane stress distribution due to rotation, may be used for the vibration analysis of non-axisymmetric rotating disks, bladed rotating disks; or other rotating machine elements which may be described by thin plate theory. 97 T h e finite element method has also been a p p l i e d to the problem of a r o t a t i n g a x i s y m m e t r i c d i s k subject to a s t a t i o n a r y transverse load. It was f o u n d that t h e r e e x i s t s a f u n d a m e n t a l c r i t i c a l rotation s p e e d at which the d i s k has no t r a n s v e r s e stiffness. material p r o p e r t i e s a n d d i s k geometry dimensional p a r a m e t e r , As a . completely d e s c r i b e the d i s k , dimensional h u b r a d i u s , corresponds a T h e rotation s p e e d , may all be d e f i n e d b y a single a n d the n o n - d i m e n s i o n a l hub radius it must be c o n c l u d e d that for a n y t h e r e e x i s t s a single c r i t i c a l v a l u e of to the fundamental c r i t i c a l rotation s p e e d . non- a Since non- which the c r i t i c a l s p e e d phenomenon is resonance related, higher critical were e x p e c t e d , a n d a l t h o u g h t h e y were not d e t e r m i n e d e x p l i c i t l y , e v i d e n c e of t h e i r e x i s t e n c e was o b s e r v e d . h i g h e r than and the the fundamental results critical speed was also The governing at speeds examined, i n d i c a t e d that the r e s p o n s e at h i g h s p e e d s s e n s i t i v e to a n y a p p r o x i m a t i o n s a n d as the d i s k T h e r e s p o n s e of the d i s k speeds is very made. e q u a t i o n s p r e s e n t e d h e r e a r e all geometry a n d loading are also s y m m e t r i c , symmetric, it seems l i k e l y that the u n s y m m e t r i c t h e o r e t i c a l a n d e x p e r i m e n t a l r e s u l t s R e f e r e n c e [16] un- of may be r e p r o d u c e d without modifications to the g o v e r n - ing e q u a t i o n s p r e s e n t e d in C h a p t e r 2. T h e lack of symmetry e x p e r i m e n t a l r e s u l t s may be d u e to s e v e r a l e f f e c t s , s u c h as non-linear effects, and i m p e r f e c t i o n s in the d i s k , been a c c o u n t e d for in this a n a l y s i s . solution has been c o n s i d e r e d in this damping, none of which F u r t h e r , o n l y the s t e a d y analysis. in the have state 98 For non-axisymmetric disks, the in-plane stress results from the QST solution show that the radial and shear stresses are not greatly affected by the introduction of slots in the disk. The circumferential stresses, though, are significantly reduced throughout the disk except at the tip of the slot. The free vibration results for a non-axisymmetric disk show the presence of a "j" nodal diameter occurring along the slot diameter, which does not alter the mode shape but causes a phase shift across the diameter. The difference in the natural vibration frequency due to this phase shift appears to become negligible in the modes characterized by a large number of nodal diameters. Only a brief treatment of the problem of non-axisymmetric disks was presented here. The topic presents many interesting areas for further detailed investigations, such as the effects of varying the number or length of slots and the effects of rotation on the natural vibration frequencies. The response of a rotating disk to a stationary load would also be of interest. non-axisymmetric As any sort of analytical solution of a non-axisymmetric problem would be extremely difficult, if at all possible, this problem best displays the capabilities and advantages of the finite element formulation presented here. To include the effects of a non-uniform thickness is a trivial extension of the method described. Furthermore, if the membrane stresses due to a thermal gradient were to be included, the finite element method also provides a consistent method for assessing these. 99 REFERENCES 1. Lamb, H. and Southwell, R.V., "The Vibrations of a Spinning Disk," Proceedings of the Royal Society (London), Vol. 99, July 1921, pp. 272-280. 2. Southwell, R.V., "On the Free Transverse Vibrations of a Uniform Circular Disc Clamped at its Centre; and on the Effects of Rotation," Proceedings of the Royal Society (London), Vol. 101, May 1 922, pp. 133-153. 3. Eversman, W. and Dodson, R.O., "Free Vibrations of a Centrally Clamped Spinning Disk," American Institute of Aeronautics and Astronautics Journal, Vol. 7, No. 10, October 1969, pp. 2010-2012. 4. Mote, C D . , J r . , "Stability of Circular Plates Subjected to Moving Loads," Journal of The Franklin Institute , Vol. 290, No. 4, October 1970, pp. 329-344. 5. Mote, C D . , J r . and Nieh, L . T . , "Control of Circular-disk Stability with Membrane Stresses," Experimental Mechanics, Vol. 11, No. 11, November 1971, pp. 490-498. 6. Kirkhope, J . and Wilson, C . J . , "Vibration and Stress Analysis of Thin Rotating Discs Using Annular Finite Elements," Journal of Sound and Vibration, Vol. 44, No. 4, 1976, pp. 461-474. 7. Kennedy, W. and Gorman, D., "Vibration Analysis of Variable Thickness Discs Subjected to Centrifugal and Thermal Stresses," Journal of Sound and Vibration, V o l . 53, No. 1 , 1977, pp. 83101. 8. Mote, C D . , J r . , "Stability Control Analysis of Rotating Plates by Finite Element: Emphasis on Slots and Holes," Journal of Dynamic Systems, Measurement, and Control, Vol. 94, No. 1, March 1972, pp. 64-70. 100 9. Benson, R.C., "Deflection of a Transversely Loaded Spinning Disk," Ph.D. Thesis, University of California, Berkeley, 1977. 10. Benson, R.C. and Bogy, D.B., "Deflection of a Very Flexible Spinning Disk Due to a Stationary Transverse Load," Journal of Applied Mechanics, Vol. 45, No. 3, September 1978, pp. 636-642. 11. Timoshenko, S. and Woinowsky-Krieger, S., Theory of Plates and Shells, McGraw-Hill Book Company, New York, 1959. 12. Cowper, G.R., Kosko, E., Lindberg, G.M. and Olson, M.D., "A High Precision Triangular Plate Bending Element," Aeronautical Report LR-514, December 1968, National Research Council of Canada. 13. Felippa, C . A . , "Refined Finite Element Analysis of Linear and Non-Linear Two-Dimensional Structures, " Ph.D. Thesis, University of California, Berkeley, 1966. 14. Cowper, G.R., Lindberg, G.M. and Olson, M.D., "A Shallow Shell Finite Element of Triangular Shape," International Journal of Solids and Structures, Vol. 6, No. 8, 1 970, pp. 1 133-1 156. 15. Cowper, G.R., Kosko, E., Lindberg, G.M. and Olson, M.D., "Static and Dynamic Applications of a High-Precision Triangular Plate Bending Element," American Institute of Aeronautics and Astronautics Journal, Vol. 7, No. 10, October 1969, pp. 1957-1965. 16. Tobias, S.A. and Arnold, R.N., "The Influence of Dynamical Imperfection on the Vibration of Rotating Disks," The Institution of Mechanical Engineers Proceedings, V o l . 171 , 1957, pp. 669-690. 101 APPENDIX A MATRICES A S S O C I A T E D WITH T H E SOLUTION OF T H E IN-PLANE PROBLEM The matrices described in Section 3.2 of this thesis are presented here in full for the sake of completeness. No readily available complete documentation of the quadratic stress triangle could be found. by Felippa [13]. One derivation has been presented The derivation described in this thesis is based on a paper by Cowper et^ a k [14], in which a shallow shell triangular finite element is reduced to a plane stress element. The geometry of this element and the body force components resulting from a rotation about the global coordinate origin are shown in Figure 3.1. The variables at each of the three nodes are U ' 3u 3x ' 3u 3y ' v ' 3v 3x ' 3v 3y * The displacement vector also includes the u and v displacements of the centroid. The integer exponents of Equation 3.1 are; m. = o. 1, o. 2, 1, o, 3, 2, 1, 0, o, o, 0, 0, o. o, 0, o. 0, 0 n. = o, o. 1, 0, 1, 2, 0, 1, 2, 3, 0, 0, 0, o. 0, 0, 0, 0, o, 0 Pj = o. o. o, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 2, 1, 0, 3, 2, 1, 0 = o. o, o. 0, o. o. o. 0, o, o, o, 0, 1, q i 1, 2, 0, 1, 2, 3. . . . .(A.l) 102 The components of the stiffness matrix, [ k ] , and consistent load vector for the body force due to rotation, {F} coordinates, are given k - E k.. ij K f r i j = m.m. F(m. + m. - 2,n. + n.) + q.q. F(p. + p.,q. +q. - 2) i j i j ' i j M^j i ' J ' ^ I j M , n.n. F(m. + m.,n. + n. - 2) + p.p. F(p. + p. - 2 , q . + q.) , ^ i j i j' i j ^rj "i *j ' M ] ) 1-v £ k h 2 1-v F. by; 2(1-v ) ^ relative to the generalized n.p. + vm.q. 1-v T i P j n + V i j J m q F(m. + p. - 1 ,n. + q. - 1) j i ' j M K F(m. + p. - 1 ,n . + q. - 1) i j i j (A.2) = phfi' C£ F(m.,n.) + F ( m . + 1,n.) i = = phfi' C F(p. q.} + F(p.,q. + 1) i = 11, • " , 2 0 = x_ cos 9 + y~ sin 9 F(m,n) n 3 # 7 c 3 n+1 I m+1 la c n , . >m+1 - (-b) 1, • • • , 1 0 = - x , sin 9 + y , cos 1 m!I nn l (m+n+2) ' 1 (A.3) (A. 4) 103 The stiffness matrix and load vector are transformed to global coordinates through the transformation and rotation matrices which are reproduced in Tables A . l and A . 2, respectively. T [K] = [R] {P} = [R] [T [T -1 '] T [k] [T Hence; -1 '] [R] ] {F } . . . .(A.5) . . . .(A.6) To obtain the in-plane stress distribution from the deflection solution vector, the transformation described by Equation 3.13 is used to obtain the generalized coordinates, a ^ which are then related to the strain (and hence, stress) distribution within each element. The matrix [AA] from Equation 3.14, which is used for this purpose, is presented in Table A . 3 . When the in-plane stress distribution is obtained by forcing it to match the exact solution at the nodes and mid-sides as described in Section 3.2.3, another transformation matrix is required. This matrix, which is used in Equation 3.18, is presented in Table A. 4. 104 TABLE A.l TRANSFORMATION MATRIX [T] 0 1 0 T T [T] 1 0 T 2 0 T 3 0 T 2 0 T T 3 0 4 0 where [T ] 2 [T ] 3 = 0 0 -b -2b 0 0 3b 1 0 -b 0 0 a 0 a 0 0 a 1 0 2a 0 0 3a 1 -b 0 b 0 1 0 0 0 1 0 2 2 _0 0 1 0 a 0 1 0 c 0 0 c 0 1 0 0 c 0 0 1 0 0 n a-b 3 c 3 (a 9 (a-b)c 9 0 3 0 | 0 0 0 0 0 0 0 0 0 0_ 0 2 b 3 2 0 0 2 0 2 0 a 2 0 0 0 0 0 0 2c 0 0 2 c 0 2 c 9 (a-b) 27 3 (a-b)c' 27 (a-- b ) c 27 2 c 27 j ~3 c 0 105 TABLE ROTATION R A.2 MATRIX 0 l cos 0 0 [ R, ] -sin 0 0 0 2 -sin 8cos 0 0 0 0 0 R 0 0 0 0 sin 8 1 0 0 cos 8 - s i n 8cos cos 8 . 2 -sin 8 . 2 sin 8 sin Bcos 0 cos sin 8cos ( 0 -sin bcos 2 cos 8 sin 2 6 R, 0 [R ] -sin 0 l sin 8cos 2. cos 8 . 2 -sin 9 -sin 8cos . 2 sin 6 0 R si n 8cos f cos 6 0 0 [R] where [R] cos 106 T A B L E A.3 STRAIN a a 2a 2 l3 a |(a a ) 3 + 1 2 DISTRIBUTION MATRIX [AA] a 4 15 2 a (a +2a ) 5 14 3a 5 l6 a (2a a ) 6 + 1 5 2a ? l8 2 a (a +3a ) 8 17 a g i9 ^20 (2a +2a ) g 18 (3a T A B L E A.4 TRANSFORMATION MATRIX FOR FORCED E X A C T S T R E S S DISTRIBUTION [T] 0 0 a 0 0 0 0 0 c (a-b) 2 a 2 (a-b) 4 0 0 2 c 4 c2 4 0 0 2 c 2 c 2 4 bj 4 : ac 4 -be g 2 1 0 + a )j 1 9 107 APPENDIX B M A T R I C E S A S S O C I A T E D WITH T H E S O L U T I O N OF OUT-OF-PLANE THE PROBLEM T h e m a t r i c e s d e s c r i b e d in S e c t i o n 3.3 of t h i s t h e s i s a r e p r e s e n t e d h e r e i n f u l l f o r the s a k e of c o m p l e t e n e s s . F o r a complete d e s c r i p t i o n of the d e r i v a t i o n of the 18 d e g r e e of freedom plate b e n d i n g e l e m e n t , see R e f e r e n c e s metry a n d c o o r d i n a t e s y s t e m s F i g u r e 3.2. w [12,15]. 3w 3y geo- u s e d f o r t h i s element a r e s h o w n in T h e v a r i a b l e s at each of the t h r e e nodes 3_w 3x The 3 w 3x3y ' 2 are 3 w 2 3y 2 T h e i n t e g e r e x p o n e n t s of the d i s p l a c e m e n t f i e l d as d e f i n e d b y Equation 3.19 a r e ; m. = 0, 1, 0, 2, 1, 0, 3, 2, 1, 0, 4, 3, 2, 1, 0, 5, 3, 2, 1, 0 n. = 0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4, 0, 2, 3, 4, 5 .(B.I) T h e components of t h e l i n e a r b e n d i n g s t i f f n e s s m a t r i x [k ] , the c o n s i s t e n t mass m a t r i x [ m ] , a n d t h e c o n s i s t e n t u n i f o r m d i s t r i b u t e d load v e c t o r {F}, r e l a t i v e to the g e n e r a l i z e d c o o r d i n a t e s , a r e g i v e n by; 108 Eh IJ = k.. 'I 12(1-v ) ,J k.. - 2 m.m.(m.-l) (m.-l) F(m. + m.-4, n. + n.) i j i J i J i j + n.n.(n.-1) (n. - 1) F(m. + m., n. + n. - 4) i j i j i j i j ' 2( 1-v)m.m.n.n. + vm.n.(m.-l) (n.-l) + vm.n.(m.-l) (n.-l) i j i j i j i j j i j i F(m. + m. - 2, n. + n. - 2) J ' J (B.2) 1 m.. = ph F. i = q F(m., n.) ^ i' i IJ F(m. + m., n. + n.) ' J .(B.3) J 1 ABA) The components of the geometric stiffness matrix, which was derived in Section 3.3.3 based on a quadratic stress distribution, a r e given by G k.. 'J = h T , b, m.m. F(m. + m. + m, - 2, n. + n. + n. ) j k ' i j k k=1 1 J 1 + c. n.n. F(m.+m.+m. , n.+n.+n. - 2 ) k i j i j k' i j k ' + d. (m.n. + m.n.) F(m. + m. + m. - 1 , n.+n.+n. - 1 ) k i j j i i j k i j k ' (B.5) 109 in which the coefficients fc> , c , and d k distribution. k k describe the in-plane stress The components of the inertia stiffness matrix as derived in Section 3.3.4 are given by 6 k!. = phfi I 2 {b m.m. F(m. + m. + m k R - 2, n ; + n. + n ) R k=1 + c, n.n. F(m. + m. + m. , n. + n. + n. - 2) k i j i j k' i j k . (B.6) . + d. (m.n.+m.n.)F(m.+m. +m. - 1 ,n. + n. + n. -1) k i j j i i j k ' i j k where the coefficients b , c , and d k k are defined by Equation 3.37. k In all of these expressions; F(m,n) = m + 1 , , ,m+1 a - (-b) c n l m!I n (B.7) (m+n+2) These matrices and vector are related to the global coordinate system through transformation and rotation matrices which are reproduced in Tables B.I and B.2, respectively. Hence, the various matrices in global coordinates; [K ] = [R] [T~ ] [k ][T" ][R] [K ] = [R] [T L C [K ] 1 [M] = = T 1 T _ 1 T L 1 ] [k ][T" ][R] T C 1 [R] [T' ] [k'][T" ][R] T 1 T 1 [R] [T" ] [m][T" ][R] T {P} = [ R ] [ T T 1 _ 1 T ] T 1 {F} . . . . .(B.8) o o o o o o o o O o o O £ O o O . J3 o o o o o o o o o o o o o o o o o o o o o o O O O O . 3rO . Q . O to rp o D O . o j O O O o o o o o o o o o o o O ro Q O O < N O T 5f M (N (D B J O ( if O o o o o o o o o O O O O O . Q O O O O O _Q l _Q ro o _Q o KO O O O O _ Q . O O _Q O Q o O O * — o r M O O o o N O t— O O o o f o o o O O O rxjo j O O O t (TJ O — o o O o O o O o O o u (J o o o o o o o o o o O O O o o u u Q u o O (J O ^ O O O O O O y o ^ O O O ( T J O r N o o o t » — o o j O O O o u o o O O U o o o o » — o O O o O u r N O o o u o o r ° o o o O » s O o o u u o o o ^ o o O o O j r- - O o J o o * j u n j O f N i o u o o O o O to O o to ro O O O nj D ( N C o O C N o _QOOO t ( O o o o r o a=r u u < > o .0 i n in o ( D n o _n o o o o o o o O j C M o y tn rM ( Q O . O O O a- 3 " ° o o O a u f O j in o o I N ro O o o o o o o o Q o o O o o O .• o o o o o o o o O O N — i o O o o O o o o O o O o o o o o f o o o o o o o o o 110 T A B L E B.2 ROTATION MATRIX cost R 1 0 0 0 R 1 sin 6 0 -sinB [R] [R] 0 [R,] cos( 2 cos 6 2sin0cos0 -sm0cos0 2 . 2 cos 0-sin 0 2 s in ( 0 .2 sin 0 -2sm0cos0 sin0cos0 2 cos (
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Finite element analysis of rotating disks
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Finite element analysis of rotating disks Nigh, Gregory Lynn 1980
pdf
Page Metadata
Item Metadata
Title | Finite element analysis of rotating disks |
Creator |
Nigh, Gregory Lynn |
Date Issued | 1980 |
Description | The finite element method is applied to the analysis of rotating disks. The equation governing the transverse deflection of a rotating disk is presented in both body-fixed and space-fixed coordinate systems. The quadratic stress triangle is used to determine the in-plane stress distribution which contributes to the out-of-plane stiffness. The out-of-plane problem is analyzed using elements based on the high precision 18 degree of freedom compatible plate bending triangle. Free vibration problems of both axisymmetric and non-axisymmetric disks are solved in the body-fixed coordinate system. The numerical results obtained for axisymmetric disks show good agreement with results from other analyses. In a space-fixed coordinate system, the response of a rotating axisymmetric disk to a space-fixed transverse load is examined in some detail. It is found that there exists a fundamental critical rotation speed at which a component of a natural vibration mode becomes stationary in space, leading to an instability. This critical speed may be related to a non-dimensional stiffness parameter which, together with the non-dimensional hub radius, completely defines the disk geometry and stiffness. Thus, it appears that the critical value of the stiffness parameter is a function only of the non-dimensional hub radius. The response of the disk past this critical stiffness is also examined. The results provide evidence of the existence of the higher critical speeds and indicate that the response is very sensitive to any approximations made with respect to the geometry of the disk or the in-plane stress distributions. |
Subject |
Disks, Rotating |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-03-22 |
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.0062857 |
URI | http://hdl.handle.net/2429/22267 |
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 |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1980_A7 N54.pdf [ 4.38MB ]
- Metadata
- JSON: 831-1.0062857.json
- JSON-LD: 831-1.0062857-ld.json
- RDF/XML (Pretty): 831-1.0062857-rdf.xml
- RDF/JSON: 831-1.0062857-rdf.json
- Turtle: 831-1.0062857-turtle.txt
- N-Triples: 831-1.0062857-rdf-ntriples.txt
- Original Record: 831-1.0062857-source.json
- Full Text
- 831-1.0062857-fulltext.txt
- Citation
- 831-1.0062857.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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062857/manifest