UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Finite element analysis of rotating disks Nigh, Gregory Lynn 1980

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

FINITE ELEMENT ANALYSIS OF ROTATING DISKS by GREGORY LYNN NIGH ( B . A . S c , The University of British Columbia, 1 978) A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Civil Engineering We accept this thesis as conforming to the required standard THE 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 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. Department of Civil Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, V6T 1W5 Canada and study. I further agree that permission for extensive copy-Gregory Lynn Nigh 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 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-axi-symmetric 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 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 TABLE OF CONTENTS Page ABSTRACT i i TABLE OF CONTENTS iv LIST OF TABLES vii LIST OF FIGURES ix ACKNOWLEDGEMENTS xi Chapter I. INTRODUCTION 1 1.1 Background 1 1.2 Purpose and Scope 3 II. 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 III. FINITE ELEMENT REPRESENTATION 18 3.1 Introduction 18 3.2 In-Plane Problem 19 V Chapter Page 3.2.1 Introduction 19 3.2.2 The Finite Element Solution 19 3.2.3 The Forced Exact Solution 25 3.3 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 3.4 Convergence and Bounds 38 IV. 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 44 4.3 Vibrations of Axisymmetric Disks 50 4.3.1 Vibrations of a Rotating Membrane 51 4.3.2 Vibrations of a Rotating Disk-Without Hub 58 4.3.3 Vibrations of a Rotating Disk-With Hub 63 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 Free Vibration Response 84 4.5.3 In-Plane Stresses 89 vi Chapter Page V. CONCLUSIONS 96 REFERENCES 99 APPENDIX A - MATRICES ASSOCIATED WITH THE SOLUTION OF THE IN-PLANE PROBLEM 101 APPENDIX B - MATRICES ASSOCIATED WITH THE SOLUTION OF THE OUT-OF-PLANE PROBLEM 107 vii LIST OF TABLES Table Page 1.1 Non-Dimensional Frequency Coefficients for Membrane Vibration 45 4.2 Thin Rotating Ring Deflection Results 48 4.3 Non-Dimensional Natural Frequencies of a Rotating Membrane 57 4.4 Natural Vibration Frequencies (Hz) of a Rotating Disk 62 4.5 Non-Dimensional Frequency Parameters of a Rotating Disk 66 4.6 Critical Stiffness from Finite Element Solution K = 0.2 7 Z* 4.7 Natural Frequencies of a Disk Rotating @^£p-. 7 6 4.8 Natural Frequencies of Non-Axisymmetric Disk 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 Page 3.1 Coordinate Systems and Body Forces for Quadratic Stress Triangle Element 21 3.2 Coordinate Systems for 18 DOF Triangular Plate Bending Element 29 •4.1 Layout of Finite Elements for Membrane Vibration 45 4.2 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 Stress Distribution K = 0.0 54 4.5(b) Non-Dimensional In-Plane Circumferential Stress Distributions K = 0.0 55 4.6(a) Vibrating Rotating Membrane, Relative Error in Eigenvalues 59 4.6(b) Vibrating Rotating Membrane, Relative Error in Eigenvalues 60 4.7 Finite Element Grid for Disk with Hub 65 X Figure Page 4.8 Finite Element Grids for Stiffness Analysis 69 4.9(a) Non-Dimensional In-Plane Radial Stress Distributions K = 0.2 70 4.9(b) Non-Dimensional In-Plane Circumferential Stress Distributions K = 0.2 71 4.10 Typical Deflected Shape Near Fundamental Critical Speed 74 4.11(a) Deflection Contours, Grid N = 3, Stresses Forced Exact 79 4.11(b) Deflection Contours a = 0.001 80 4.12 Non-Axisymmetric Disk Geometry 83 4.13 Finite Element Grid for Non-Axisymmetric Disk .. 83 4.14(a) Static Deflection of a Non-Axisymmetric Disk . . . . 85 4.14(b) Static Deflection of a Non-Axisymmetric Disk 86 4.14(c) Static Deflection of a Non-Axisymmetric Disk . . . . 87 4.15 Finite Element Grid for Non-Axisymmetric Disk .. 88 4.16 Typical Vibration Modes of a Non-Axisymmetric Disk 91 4.17 Non-Dimensional In-Plane Stress Distributions of a Non-Axisymmetric Disk 92 4.18 Non-Dimensional Stresses Along Slot Diameter .. 94 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. He would also like to thank Dr. S.G. 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 for a complete disk were used, but the out -of -plane displacement and slope were constra ined at the hub radius. Aga in , approximate results were obtained for the case in which both f lexural effects and membrane stresses contr ibuted to the out-of -p lane st i f fness. Eversman and Dodson [3] cons idered the free transverse vibrations of an axisymmetric constant thickness disk with a central hub, with both in-plane and out-of -p lane displacements fully c o n -strained at the hub. The i r analysis was free of any approximations to the theory, although a numerical solution of the differential equation was requ i red . Non-dimensional graphical results were presented which enable determination of the natural v ibrat ion frequencies corresponding to modes character ized by 0, 1, or 2 nodal diameters or c i rc les . Mote [4] and Mote and Nieh [5] have presented both theoretical and experimental investigations into the phenomenon of cr it ical speed instabi l i ty. Th i s instabil ity occurs when the wave speed of the backward travel l ing component of one of the vibrat ion modes is the same as the rotation speed. Th i s leads to a resonant vibrat ion mode which is stationary in space, and is an important -factor in the design of thin rotating d i sks . In recent years the finite element method has been appl ied to the problem of free v ibrat ions of rotating d i sk s . Kirkhope and Wilson [6] and Kennedy and Gorman [7] have used efficient annular elements to examine the vibrat ions of thin axisymmetric d i sk s . 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 trans-verse load has received relatively little attention. 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-of-plane deflection of an axisymmetric disk due to a transverse load. This problem is examined in some detail. Some results from non-axisymmetric problems are also examined. 5 CHAPTER II THEORETICAL BACKGROUND 2.1 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 body-fixed 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 coordin-ate 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. Then, 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. Thus, in polar coordinates; pQ2 r S = 0 . . . . (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 = p f i 2 x Y = pQ2y . . . .(2.2) and the differential equations for in-plane equilibrium are; 3 a 3 a * + ^1 + X = 0 3 x 3 y 3 a 3 a 2£ + 2- + Y = 0 . . . .(2.3) 3x 3y 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. Thus, 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 ] ; DV w-h 3 2w _ 3 2w 3 2w o + 2a + a x ... 2 xy y „ , w 2 in which; w h D V q 3x ' x3y 3y transverse deflection thickness plate rigidity Eh" + h x 3J5L + Y 3x 3y I 2 ( 1 - v z ) biharmonic operator transverse load per unit area. = q (2.4) 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 D V w-h 3 2w ' x 3 x 2 + 2a xy 32w 3 x 3 y + a 3 2w y 3 y 2 . , v 3 w ,, 3 w + h X + Y 3x 3y = q - ph 3 2w 3t 2 (2.5) 8 9 2w where is the acceleration in the transverse direction. 3 t 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. Then, 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. If q is time dependent. Equation 2.5 is difficult to solve. Furthermore, if the load is stationary in space, it is likely the response in space which is of interest. Thus, 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. For two dimensional problems, this operator is defined by: where v x and are the velocities of the body in the x and y directions, respectively. For the problem described herein; v x = - f i y v y = fix (2.7) 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 , , , . 9 2w T , r, —=• (phw) = ph + 2phfi Dt' 9t 9 2w 3 2w {' y 9x 9t + X 9y at J ( 2 9 w 3x' 2xy 3 2w 2 9 2w + x 9x9y 9 y - phfi : 9 w 9 w 9x (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 DV  hw - h 3 2w _ 3 2w 3 2w ' j + 2a + a X 3x 2 X y 3x3y V 3 y 2 i + phfi' ( 2 3 2w 3x' 3 2w 2 3 2w - 2xy + x 3x3y 3y' . 3 w - . _. = q - ph + 2phQ 3 t 2 3 2w 3 2w 3x3t 3y3t (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 body-fixed 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; DV 4w - h 3 2w * 3x 2 + 2a xy 3 2w 3 2w + a 3x3y y 3y 2 > + phfi : '2 3 2w y 1? 2xy 3 2w + x 2 3 2w 3x3y 3y 2 = q (2.10) 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 co-ordinate 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, which is re la t i ve to a space - f i xed coord inate sy s tem, q is a funct ion only of po s i t i on , whereas it is a lso a func t i on of time in Equat ion 2.5, which makes the problem more d i f f i c u l t . T h u s , both equat ions wi l l be used in th is thes i s , depend ing on the pa r t i c u l a r problem to be s o l ved . 2.4 Ene rgy Funct ionals and Bounda r y Cond i t i ons A s the f in i te element representat ions of the problems de sc r i bed here in r equ i re the use of the energy funct iona l s , they are p re sented h e r e . They are based on s t ra in energy cons iderat ions and the p r i n c i p l e of minimum po ten -t ia l e ne r g y . These funct iona l s are also used to de r i ve the bounda ry c o n -d i t ions assoc iated with the d i f f e ren t i a l equat ions g o v e r n i n g the ou t - o f -plane prob lem. When the in -p lane and ou t -o f -p l ane systems are uncoup led , and the rotat ion speed is constant so that the in -p lane body forces are independent of time, the i n -p lane problem may be t reated as a stat ic prob lem. Then the total potent ia l energy of the i n -p lane sy s tem, n, is the s t r a i n ene r g y , U, p l u s the potent ia l energy of the body f o r ce s , W. Thu s n = u fw (2.11) where U Eh 2 ( 1 - v 2 ) 2 2 1 - v , ,2 u + v + 2vu v + (u + v ) x y x y 2 y x d xdy . . .(2.12) 13 W Xu + Yv dxdy (2.13) 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. The principle may be written as r*1 L dt = 0 (2.14) 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 B D 2 2 2 w +w + 2vw w +2(1-v)w^ w l [ xx yy xx yy xy j dxdy (2.15) U M 2 2 a w + a w +2a w w , x x y y xy x y j dxdy (2.16) 14 The potential energy of the lateral loads is given by; W - - q w dx dy . . . .(2.17) The kinetic energy associated with transverse motion of the disk is; 2 1 T 2 ph \ ^ } d x d y . . . .(2.18) where is the transverse velocity of the disk in body fixed coordinates, a t Thus, the Lagrangian function is; 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; T + U, 15 (2.20) Using the relations of Equation 2.7, the inertia integral, U., is; U, p h f i : 2 2 2 2 . y w + x w - 2xyw w 1 x y 7 x y dxdy . . . .(2.21) and p h 7 x y dxdy .(2.22) where w - . The equation of motion. Equation 2.9, may then be o t 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 1 1 = U B + U M " U l + W .(2.24) 16 The boundary integrals for Equation 2.10 may also be obtained in the same way. The resulting boundary integrals (which must vanish), in x-y coordinates, are; 9 6 w 2 2 2 phft y w + phft xyw + N w + N w M 7 x 7 y x x xy y D( w + (2-v) w ) , xxx x y y J dy = 0 9 6 w D (w + v w ) , xx yy j dy = 0 9 6 w 2 2 2 phQ x w + phft xyw + N w + N w H y x y y x y x - D(w + (2-v) w ) , yyy * *y J dx = o 9 6 w D(w + vw ) yy xx dx = 0 (2.25) where N , N , and N are the in-plane stress resultants and 6w, x y xy ^ 6w , and 6w are the arbitrary variations of displacement and slope x y 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 must be expressed in polar coordinates. As Equation 2.10 is applicable only to axisymmetric disks, this is not a restriction. Thus, the boundary integrals become; 1 1 d M r Q 6 w i N w + N „ - w Q - C ) + ; r o |J I r r 'V9 r " 6 r 3 0 de <> 6 w M ^ d0 = = 0 . . . .(2.26) J r r ^ J where Q is the transverse shear and M and M „ are the bending and r r r0 a twisting moments, respectively. The vanishing of these boundary integrals then leads to the general boundary conditions along any edge as; either 1 1 8 M re N w + N Q - w Q - Q + - = 0 or 6 w = 0 r r r 6 r 0 r r 3 6 and either M = 0 or 6w = 0 . . . . . ( 2 .27 ) r r For the disks considered here, the inner radius is clamped, allowing no deflection or normal slope, so that 6w = 6w r = 0 . On the outer radius, the boundary is free, so that N r and N r Q 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 III FINITE ELEMENT REPRESENTATION 3.1 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 con-ditions 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 displace-ment 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 differen-tial equation and energy integrals, some distribution of these stresses within each element must be determined before the out-of-plane problem may be solved. Thus, 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, ' 3u 3u 3v 3 v u ' 3x ' dy ' v ' 3x ' 3y at each of the three nodes. 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 co-ordinates and the generalized coordinates a.; 20 20 m. n. p. q. u = I a. E, n v = I aj E, n . . . .(3.1) i=1 i=1 where m., n., p. and q., integer values ranging from 0 to 3 inclusive, are given in Appendix A . A matrix relationship between the nodal displacements (in £-n. coordinates) and the generalized coordinates, may then be developed; w' - \ = [T] {A [ . . . .(3.2) where {A} is the vector of generalized coordinates, [T] is a transforma-tion matrix which is a function only of the element 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 displace-ments in global coordinates; W L 1 = [R] ( w C \ . . . .(3.3) where {W C} 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} T[k] {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 {WG} [R] [ T _ 1 ] [k] [ T _ 1 ] [R] {WG} . . . . ( 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 2 r , which may be expressed in x-y coordinates as X = p f i 2 x and Y = pQ2y. These x-y body forces may be transformed to element local coordinates to give 23 x = p f i 2 ( c ? + O Y = p f i 2 ( c ^ + n) in which the constants (x,cos6 + y,sin8) c = (-x^inS + y^osB) (3.7) 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 0 f W = p h f i 2 Z a. i=1 ' u m. n. m.+l n. (ccti 'n 1 + E 1 n ') d^dri 20 + p h f i 2 T a. i=1 1 p. q. p. q.+ l ( c ^ 'n 1 + K 'n 1 )d^dn (3.8) Recognizing Euler's Beta function; F(m,n) K n d^dn n+1 ( m+1 , u^ m + 1 l m!n! : | a " ( " b ) j (m+n+2) (3.9) the potential energy of the body forces may be expressed as the vector product 24 W = { F }T { A } .(3.10) in which the elements of the vector {F} are given by F. = phft2 CpF(m.,n.) + F(m.+1,n.) i - 1,2, — , 1 0 = ph fi2 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 in global coordinates; 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. displacements of each element may be determined from the element nodal Q displacements. With the full element displacement vector {W }, the {P} = [ R ] T [T ] ] {F} . . (3.12) The stiffness matrix and load vector are normally statically Once the global problem has been solved, the centroidal 25 generalized coordinates within each element may be determined through another transformation; {A} = [T _ 1][R]{W C} . (3.13) Finally, the quadratic stress distribution, with respect to the local axes, may be calculated as a function of the generalized coordinates through [E] 30 3 5 3v 3n 3u 3v 95 [E] [AA] e K n 2 n (3 . 14) where [E] is the elastic constitutive matrix and [AA] is a function of the generalized coordinates which is reproduced in Appendix A. 3 .2 .3 The Forced Exact Solution If the problem geometry is axisymmetric, the in-plane stresses become functions only of the radial dimension, and are independent of the angular coordinate. As a result, the in-plane shear stress is identically zero. If the boundary conditions require that there be no displacement at the inner radius, a, and no surface tractions on the outer radius, b, the in-plane stresses are described by; a = C (3 + v) r 1 26 C j p + v) 1 - C - -C ' 1 + 3 v ,2 3+v where and r 1-v 3+v (3+v)-(1+v)K { (1+V)+(1-V)K 2 J 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) [1 - c ] C^3 + v) 1+3v 3+v (3.16) In order to generate the stress distribution, the value of a and oa 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 a j C ; 'n 1 . . . .(3.17) i=1 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 co-ordinates and is the same for each of the stresses [ar,o ,o r ) to be considered. It is reproduced in Appendix 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?r and Og is only quadratic and can be described exactly by quadratic variations of , and a ^ T ] -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 Out-of-Plane Problem 3.3.1 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 DV 4w is associated with the linear bending stiffness of a thin plate; the terms containing the in-plane stresses 28 are associated with the geometric stiffness, and the terms proportional to phfi2 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 dis-placement 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 previously considered QST. 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 QST, this element is derived in a local E-ri co-ordinate system shown in Figure 3.2. In this case, the assumed out-of-plane displacement, w, is given by an incomplete quintic polynomial; 20 m. n. w = I a.E 1 n 1 . . . .(3.19) i=1 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 con-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 dis-placements 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 B . Applying the coordinate transformation and rotation matrices as before gives the element strain energy as a quadratic function of the eighteen nodal variables; T T T U = 1 {WC} [R] [ T _ 1 ] [k ] [T _ 1 ] [R]{W C } . . . . .(3.21) Q The vector of nodal displacements,{W } consists of 2 2 2 9 w 9 w 3 w 9 w 9 w W ' 3 x ' 3 y ' 9x2 ' 3x3y ' 3 y 2 at each of the three nodes. 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 out-of-plane problem. To derive the geometric stiffness matrix, the assumed dis-placement field and the in-plane stresses are written in a series form; 20 m. n. w = £ a. E 1 n ' ; . . . . ( 3 . 2 2 ) i = 1 1 N N m. n. m. n. o> = I b. E 1 n 1 a = I c E ' n 1 ^ i=1 ' 1 i=1 N m. n. ov = T d. E n . . . .(3.23) where a., m., and n. are as before. The variables b., c , and d. are i' i' i i' i' i the coefficients which describe the in-plane stresses ov, a , and ar ^ E n En respectively, with positive directions as shown in Figure 3.2. 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 displace-ment, and combining with the in-plane stress distributions, each term in the energy integral (Equation 2.16) may be expressed as a triple sum; 8 w 195 J 20 20 6 m - k ~ 2 " i k 7 T T a.a.b. m.m. £ n • • • • i i i 1 J k 1 J i=1 j=1 k = 1 ' 1 9 w 9n 20 20 6 m nj-|<~2 y y y a.a.c. n.n. E, I J n ' • i • i i i 1 J k 1 J i=1 j =1 k = 1 '5n ' 9 w ' ' 9 w ' [ 9n J I I I a i a i d k m i n i K n j=1 j=1 k=1 I ' m... = m. + m. + m. ijk i j k n... = n. + n. + n. . ijk i j k (3.24) Using these relations and the requirement of symmetry, the strain energy may be expressed as; U 20 20 6 I I I \ i=1 j=1 k=1 m... -2 n... a.a.b. m.m. E, ' r\ ' ' J k i j m : : I. " ,-2 m-1 n,.,-l) + a.a.c, n.n. £ j j k n i]k~ + a i a j d k (m.n . + m.n i ) 4 ' j k n i j k | d £ d n . .(3.25) 33 Finally, taking the integration inside the summation, and recognizing Euler's Beta function to give F(m, n) f rm n . n + 1 / m + 1 . . *m+l\ m!n! , _ . E n dEdn =c ja - (-b) | ( m + n + 2 ) j • • • -(3.26) the geometric stiffness matrix relative to the generalized coordinates may be determined as; 6 f k.. = h T <b, m.m.F(m... -2,n... ) + c. n.n.F(m... ,n... -2) ij k-i I k 1 J ' J k ' J k k i j ijk' ijk + d, (m.n. +m.n.)F(m... - 1,n... - 1) > . . . .(3.27) k i j j i ijk ijk 1 Applying the same transformation and rotation matrices as 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 ] T [ k ] [ 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 may be easily derived by analogy to the geometric stiffness. Since the geometric stiffness matrix was developed by considering the strain energy in a local co-ordinate 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 transformation must be used than just exchanging £ 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 . 2 9 ) from which the inverse relations may be obtained; £ = - x Q cos 8 - y Q sin 9 + x cos 6 + y sin 0 n = x _ s i n 8 - y 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 deriva-tives of w in local coordinates, so that; ( 2 92w _ 8 2w 2 3 2w y - 2xy + x -9x' 3x9y 9y' (y cos 0 - x sin 0) 2 92w 9E< + 2(x sin 6 - y cos 6) (x cos 8 + y sin 0) 9 2w 9 £9n + (y sin 0 + x cos 0) 2 9 2w 3 n 2 . (3.31) Using the Equations 3.29 to express the x-y coordinates in terms of £ and r\ gives; f 2 2 2 2 9 w _ 9 w 2 9 w y Y ~ 2xy + x 9x 9x9y 9y 2 j (y 0 cos 0 - x Q sin 0 + ri) 2 9 2w 3 ^ + 2(x Qsin 0 - y Qcos 0 -n ) ( y o s i n0 + x Qcos 0 + £) 92w 3 E 3 n 32w + (y n s in 0 + x n cos 0 + £) ^ 3 £ (3.32) Thus, the inertia terms of the differential equation are expressed correctly in local E-ri coordinates. 36 Considering the geometric stiffness terms in local coordinates; 2 2 2 9 w _ 9 w , 3 w f , n + 2oV._ + 0_ =;— . . . . 1 3 . 33J ? 9 £ ^ 9 £ 9 n n 8 n 2 the following analogies can be made; 2 (y Q cos 8 - x Q s i n 6 + n) is a (x Q sin 9 - y Q cos 6 - n) (y Q sin 6 + x Q cos 6 + £) is a 2 (y Q sin 8 + x Q cos 8 + £) is . . . . .(3.34) On the basis of this analogy, 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 m. n. 6 m. n. \ - I bj K ' n ' a = I c. K ' n 1 ^ i=1 1 n i=1 6 _ m. n. ov = J d. E n . . . . (3.35) ^ i = 1 1 37 so that the inertia stiffness matrix relative to the generalized coordinates is 2 6 r -k.. = phQ T \b. m.m. F(m... -2, n... ) +c. n.n. F(m... , n... -2) ij ^ ( 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. and n... = m. + m. + n. . . . (3 .36) ijk i j k ijk i j k The coefficients b\, c., and dj are zero except for the following; — 2 — -b 1 = (y Q cos 6 - x Qsin r8) b 3 = 2(y Qcos 6 - x Q s i n 0) b & = 1 2 -c i = ( y 0 s i n 9 + x o c o s 9 ' c 2 = 2 ( y Q s i n 9 + x o c o s 9^ c 4 = 1 d 1 = (x Q sin 6 - y Q cos 6) (x Q cos 6 + y Q sin 0) d 2 = ( x Q s i n0 - y Q cos 0) d^ = - (x Q cos 0 + y Q sin 0) d 5 = 1 . . . . . (3 .37) The transformation and rotation 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 stresses), the inertia stiffness matrix makes the structure more flexible. Even though both the 38 geometric and inertia stiffness components are directly proportional to p f t 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 represents a lower bound on the exact strain energy . Furthermore, because only continuity of deflections (and not slopes) is r equ i red , it can be shown that the asymptotic convergence rate with respect to the g strain energy will be of the order of (1/N ). Unfortunate ly, the c o n -vergence proof imposes the rather str ict requirement that the in-plane stresses satisfy the in-plane equi l ibrium equations and boundary c o n -dit ions, and also be continuous across element boundar ies. In genera l , the in-plane stresses do not meet these requirements, so convergence cannot be proven and there are no bounds on the solutions. The inertia st i ffness matrix is represented exactly except for the er ror in the assumed displacement f ie ld. T h u s , by analogy with the geometric stiffness matrix, the convergence rate should be of p the order of (1/N ). As this rate represents faster convergence than that of the linear bending st i f fness, it should not adversely affect the convergence rate of the total out -of -p lane problem. In app ly ing such high precis ion finite elements to the problem of analysis of c i rcu lar d i sks , the rate of convergence is governed not so much by the assumed displacement field as by the modelling of the problem geometry. Modelling a c u r v e d edge with a series of stra ight line segments introduces an error in area representat ion. It can be 2 shown that this error decreases as (1/N ), where N is now the number of straight line segments. If the strain energy density is assumed constant throughout the problem domain, then the error in strain energy is proportional to the er ror in the area. T h u s , the factor govern ing 40 convergence of the problems examined herein is the representation of the problem geometry, which leads to an expected convergence rate 2 of the strain energy of the order of about (1/N ). 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 con-vergence can occur under ideal conditions. As a check that there are no fundamental errors in the program, a rather unique problem is devised 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. Typically, this is obtained by solving the in-plane problem 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. Thus, 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 con-ditions 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 Test Problems 4.2.1 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. The potential energy 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. Then, 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 2 2 , 2 a 2, 2 A = to = (m + n ) TT a n 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 g r i d . For each grid size, the frequencies associated with modes (1,2) and (2,1) were identical. 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. A line of slope 1/N is plotted for reference. The numerical results confirm the theoretical expectations, and in fact the convergence rate for the lowest eigen-g values was indeed 1/N above N = 2. 4.2.2 Combined Geometric and Inertia Stiffness Matrices The accuracy and correctness of any newly formulated numer-ical solution technique remain suspect until the results have been checked against known results. In the case of the inertia stiffness matrix, this presents a difficult problem, as no exact solutions utiliz-ing 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 Figure 4.1 Layout of Finite Elements for Membrane Vibrat ion T A B L E 4.1 Non-Dimensional Frequency Coeff ic ients for Membrane Vibrat ion Mode N = 1 N = 2 N = 4 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 4 7 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. The x and y axes are space fixed so that Equation ( 2 . 1 0 ) 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 2 2 Og = pR fi or = a r ( Now consider a small element of the ring located on the x axis as shown in Figure 4 . 3 . At this particular point, N x - o r h = 0 , 2 2 N = a J i = 0 , and N = a„h = phR fi . Furthermore, at this point xy r6 y 8 y = 0 and x = R, so that Equation ( 2 . 1 0 ) becomes 4 2 2 2 2 DV w - phR fi w + phR fi w - q yy yy or 4 DV w Thus, in the limit that the rotating disk becomes a thin rotating annular ring, the governing differential equation approaches the static plate y 48 uniform ring width b thickness h see Detail b h N 111] N Detail Element of Ring on x Axis Showing Finite Element Grid Figure 4.3 Thin Rotating Ring TABLE 4.2 Thin Rotating Ring Deflection Results Central Deflection Finite Element (inches) Bending Stiffness Only I 1.3976 Bending & Geometric Stiffness 0.6643 Bending, Geometric & Inertia Stiffness I 1.4049 Beam Theory (inches) 1.3958 0.6630 1.3958 49 equation for small elements of the r ing . So if the program is used to analyze a small element of the r i n g , the solution should be compar-able to the solution obtained us ing only the bending st i f fness. Such a problem was analyzed using a coarse finite element gr id of four elements arranged as shown in F igure 4.3. The ratio of dimensions h : b : L : R was 1:10:100:10000, so that the element dimensions were small compared to the radius, enabl ing the app rox i -mations to be reasonably accurate. Since the example depends on the cancellation of terms, it is important that all terms relating to the st iffness be of the same order of magnitude, so that the cancellation will have a significant effect on the solution. If the plate is assumed to be simply supported on the edges y - ± ^ and free on the other two, the deflection may be approximated by w = w^cos so that the govern ing dif ferential equation becomes r _ , T T , 4 , , - . 2 2 . 7 1 , 2 2 2. TT , 2 D (j-) - phR Q (j-) + phR Q, (j-) Try W 0 C O S L On this basis, the coefficients were adjusted so that D ^= = phR 2 f t 2 L and for s implicity, a uniformly d i s t r ibuted t ransverse load was used. The results of this investigation are d isp layed in Table 4.2. The deflections shown are the averages of the two nodes on the x-ax is 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. The last case was actually done four times, 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]. This provides a means of assessing the convergence characteristics for cases in which the in-plane stresses are not 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 for the natural frequencies of transverse vibration of a rotating disk for the case in which the flexural rigidity is negligible. 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. The non-dimensional frequency parameter X is given by 52 2 to x = -n = ST where n is the number of nodal c i r c l e s , s is the number of nodal d iameters , and is the natura l f r e q u e n c y . The geometric s t i f fnes s matr ix and the cons i s tent mass matr ix were used to obtain the natu ra l f requenc ie s and v i b ra t i on modes u s ing the f i n i te element g r i d s shown in F i gu re 4.4. The in -p lane s t resses were generated u s i ng the quadra t i c s t re s s t r i ang l e . The ou t - o f - p l ane v i b r a t i on problem was so lved twice for each case, one f o r c i n g a n t i -symmetry along both coord inates axes ( A - A ) and the other f o r c i n g an t i - s ymmet ry a long one ax i s and symmetry a long the other ax i s ( A - S ) . Th i s al lowed all modes to be generated except for those with no nodal d iameters . In the so lut ion obta ined by Lamb and Southwe l l , these modes i nvo l ved a d isplacement of the cent re of the d i sk and so are not inc luded for comparison he re . The exp re s s i on s fo r the i n -p l ane s t resses for the case in which there is no cen t r a l hub are g i ven in Chap te r 3. The shear s t re s s is i dent i ca l l y ze ro , and the radia l and c i r cumfe ren t i a l s t res ses a re only quadra t i c funct ions of the r ad i u s . S ince the Quadra t i c S t re s s T r i ang l e is capable of modell ing exac t l y a quadra t i c s t re s s d i s t r i b u t i o n , any 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 model l ing a c u r v e d bounda r y wi th s t r a i gh t l ine segments. The non-d imens iona l rad ia l and c i r cumfe ren t i a l s t re s se s , - ^ - ^ — , for each of the three p f t 2 b 2 f i n i te element g r i d s and the exact va lues , are shown in F i gu re s 4.5(a) (s + 2n) (s + 2n + 2) { 3+v^ 2 fl+3v ' I 8 J — s { 8 j N = 1 N = 2 N Figure 4.4 Finite Element Grids for Disk Without Hub Figure 4. 5(a) Non-Dimensional In-Plane Radial Stress Distr ibut ions 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 grid, the poor representa-tion of the problem boundary significantly affects both of the stresses throughout the entire domain. However, with the finest grid, 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 grid, together with the exact values, for the first twenty-t modes are listed in Table 4.3. These values are for a Poisson's ratio of 0.3. The modes were determined by inspection of the eigen-vectors. 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. Thus, 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 T A B L E 4.3 Non-Dimensional Natural Frequencies of a Rotating Membrane Mode* A From Finite Elements A Exact N = 1 N = 2 N = 3 (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 - number of nodal c irc les s - number 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 grid. 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 Size Figure 4.6(a) 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 vibrat-ing 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 PCF. 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 fre-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. The non-dimensional parameter 2 2Eh a = 2 3(3+v) ( 1-v ) p f i 2 b " 62 T A B L E 4.4 Natural V ibrat ion Frequencies (Hz) of a Rotating Disk Mode* Lower Bound From Finite Elements Upper Bound N = 1 N = 2 N = 3 (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) - 111.2 105.3 104. 1 104.8 (1,2) 106.6 117.5 111.2 110. 0 -(0,6) - 155.7 143. 4 141.7 143.4 (1,3) 155.5 160.0 160.2 158. 2 -(2,1) 168.4 - 180. 5 177.7 •-(0,7) - - - 185. 3 188.6 (1,4) - - - 214.2 -(0,8) - - - 235. 1 240.4 (2,2) 240.1 - - 244.9 -(1,5) - - - 277.6 -(0,9) - - - 290.8 298.8 (2,3) 303.7 - - 320. 7 -(3,1) 335.2 - - 340.6 -(1,6) - - - 348.6 -(0,10) - - - 352. 5 363.7 Mode (n,s) n - number of nodal c irc les s - number of nodal diameters 63 varies from zero, for the case of a pure membrane in which only rotational effects are present, to inf in i ty, for the case of a pure plate in which only f lexural effects are present. In this part icular example, a = 0.0639 which represents a fair ly st iff d i sk. If the disk has no central hub , then the mode (0,1) represents rotation in a plane sl ightly incl ined to the reference plane, so that the natural f requency associated with this mode is the rotation f requency . This character ist ic was also exhibited in the case of the membrane problem solved prev ious ly , in which X - 1 for this mode. As exact solutions are not available for the other modes, it is di f f icult to draw any hard conclusions from these results . However, for the modes which have bounded f requenc ies , the finite element solutions appear to converge to values within the bounds , again indicating that a reasonable degree of accuracy may be obtained even with a fair ly coarse finite element g r i d . 4.3.3 Vibrat ions of a Rotating Disk - With Hub The problem of free t ransverse vibrations of a rotating disk is now examined for the case when the disk has a central hub . Eversman and Dodson [3 ] have presented results for the natural frequencies of such a d isk, for various non-dimensional hub ratios and st i f fnesses, for all combinations of 0, 1, or 2 nodal diameters or c i rc les . The natural f requency is related to the other factors through a non-dimensional f requency parameter, defined by 64 1 8 f " V 2 + ' 1+3v ' (3+v)a I « J 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.44. 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, which enabled all possible modes to be generated. 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 fre-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 grid. Fiqure 4.7 Finite Element Grid for Disk With Hub TABLE 4. 5 Non-Dimensional Frequency Parameters of a Rotating Disk; K = 0.44, a = 0.001, v = 0.3 Mode / 7 (2) (0,0) (0,1) (0,2) (0,3) (0,4) (0,5) (0,6) (0,7) (0,8) (1,0) (1,1) (1,2) (1,3) (0,9) (1,4) (0,10) (1,5) (0,11) (1,6) (0,12) (2,0) (2,1) (2,2) 8. 24 8.67 9.68 10.87 12.07 13.24 14. 36 15.43 16. 38 13. 76 13. 90 14.32 14.95 17. 50 15.68 18.47 16.60 19.49 17. 59 20.42 not observed not observed not observed 8.2 8.6 9.7 13. 5 13. 7 14.2 18.8 18.9 19.0 Mode (n,s) n - number of nodal circles s - number of nodal diameters 1 From Finite Elements "From Reference 3 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 con-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 analyzed. 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. Thus, in order to remain consistent with that particular analysis, 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 Stress Distributions K Circumferential = 0.2 72 P non-dimensionalized by multiplying by . The maximum values p f i 2 b 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 approxi-mation 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 extra-polating from a plot of the determinant of the stiffness matrix versus a, to determine when the determinant became zero. Thus, 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 . This deflected shape, which was the same in each of the four cases considered, is similar to the vibration mode ( 0 , 2 ) . The con-tour intervals in the deflection plot represent equal proportions of the maximum deflection. Thus, 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 10% 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 s = ft ± , where t o n is the natural frequency corresponding to the mode, in a body-fixed coordinate system, ft is the rotation speed, and s is the number of nodal diameters associated with the mode. It is obvious that there may be some combination of ft, to and s such that a j g is zero. 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 2 defines the criterion of instability as ( t o ^ ± sft) = g , where 3 is the harmonic excitation frequency of the transverse load. 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 a CR N = 2 Forced Exact ( K = 0.0) 0.0707 [+00.0002 N = 2 QST Generated (K = 0.2) 0.0814 ± 0.0002 N = 3 Forced Exact (K = 0.0) 0.0805 ± 0.0002 N = 3 QST Generated (K = 0.2) 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 n at a rotation speed corresponding to ft The grid N = 2 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. The material constants of steel were used in the analysis, 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 in-plane stresses was 332 radians per second. The vibration modes and natural frequencies for this particular disk rotating at ft ^ are listed in Table 4.7. The values of ft, co and s corresponding to the mode (0,2) 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. T A B L E 4.7 Natural Frequencies of a Disk Rotating @ fi Mode 1 2 di 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 - number of nodal c irc les s - number of nodal diameters 2 F o r fi = 332 rad/sec 77 These results confirm that the instability observed in a space-fixed 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. Thus, given a non-dimensional hub radius, there exists a single value of the critical non-dimensional 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 grid, which was just slightly less than the fundamental crit-ical value. The results were identical, to seven significant figures, with those of the half grid, thus confirming the symmetry assumption. Hence the present analysis did not reproduce the unsymmetric results reported in Reference [16]. These results were obtained by using a Choleski type de-composition of the total stiffness matrix in order to solve the transverse deflection problem. This type of routine fails if negative terms occur on the diagonal. Thus, such a solution technique will not work at 78 stiffnesses past the critical stiffness, a „ D . However, from either a 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). All of the plots, except at a = 0.1, 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 ( e -9 - ' a = n o t shown), which would suggest that these particular values of 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 2 h 2 w = ( p Q ) 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 0.1 4.40 0.66 0.05 3.51 0.01 -3.06 -4.01 0.005 ±3.25 -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. Thus, the trend in Mote's results is also displayed here. 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 represents the combined geometric and inertia stiffness matrix. The critical values of a could then be related to the 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 igu re 4.11(a) Def lec t ion C o n t o u r s , G r i d N = 3, S t resses Forced Exac t N = 3 Stresses Forced Exact N = 2 Stresses Forced Exact Figure 4.11(b) Deflection Contours a = 0.001 81 disk in-plane stress distribution), which represent solutions to exactly the same problem considered here except that it is solved using a different technique. 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 in-plane 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 grid. The change due to the different in-plane stress distributions is not quite as obvious, but is still significant. Thus, 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) in which the flexural stiffness is almost negligible. Since the stiffness is then so dependent on the in-plane stresses, it would be expected that the response depends significantly on any approximations used for these stresses. Furthermore, since the (destabilizing) 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 un-stable, 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. This accounts for the lack of symmetry in the finite element grid. The load used in the experiment was 26.1 pounds, applied at a radius of 9.55 inches. For the finite element solution, the load was applied at the nodes on the outer radius. The twelve load points are indicated on the grid in Figure 4.13. 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. Figure 4. 14(a) Static Deflection of a Non-Axisymmetric Disk CO Figure 4.14(b) Static Deflection of a Non-Axisymmetric Disk T 2 3 4 5 6 7 8 9 1 0 II 1 2 Position 0 0 Figure 4.14(c) Static Deflection of a Non-Axisymmetric 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, in the order in which they occurred in the analysis. 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 in-plane stresses in the non-axisymmetric problem, even the finite element solution to this problem is of interest. The non-dimensional stress distributions are shown in Figure 4.17. These distributions 8 a - pQ2b2 -were obtained using the quadratic stress triangles with the grid shown 90 TABLE 4.8 Natural Frequencies of Non-Axisymmetric Disk * Mode co (Hz) n * Mode co n (Hz) (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 Stress Figure 4.17 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, circumfer-ential 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. The circumferential stress, 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 axi-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 improvement over those shown in Figure 4.18. 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. Thus, 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 The finite element method has also been applied to the problem of a rotating axisymmetric disk subject to a stationary t ransverse load. It was found that there exists a fundamental cr it ical rotation speed at which the disk has no t ransverse st i f fness. The rotation speed, material propert ies and disk geometry may all be defined by a single non-dimensional parameter, a . As a and the non-dimensional hub radius completely descr ibe the d i sk , it must be concluded that for any non-dimensional hub rad ius , there exists a single crit ical value of a which corresponds to the fundamental cr it ical rotation speed. Since the crit ical speed phenomenon is resonance related, higher cr it ical speeds were expected, and although they were not determined expl ic i t ly , evidence of their existence was observed . The response of the disk at speeds higher than the fundamental cr it ical speed was also examined, and the results indicated that the response at high speeds is very sensitive to any approximations made. The govern ing equations presented here are all symmetric, and as the disk geometry and loading are also symmetric, it seems u n -likely that the unsymmetric theoretical and experimental results of Reference [16] may be reproduced without modifications to the g o v e r n -ing equations presented in Chapter 2. The lack of symmetry in the experimental results may be due to several ef fects, such as damping, non- l inear ef fects, and imperfections in the d i sk , none of which have been accounted for in this analys is . Fu r ther , only the steady state solution has been cons idered in this analys is . 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 character-ized 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 non-axisymmetric disk to a stationary load would also be of interest. 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 . , Jr. 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, Vol. 53, No. 1 , 1977, pp. 83-101. 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," Aero- nautical 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, Vol. 171 , 1957, pp. 669-690. 101 APPENDIX A MATRICES ASSOCIATED WITH THE SOLUTION OF THE 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. One derivation has been presented by Felippa [13]. The derivation described in this thesis is based on a paper by Cowper et^  ak [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 3u 3 u 3v 3 v U ' 3x ' 3y ' v ' 3x ' 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 q i = o. o, o. 0, o. o. o. 0, o, o, o, 0, 1, 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} relative to the generalized coordinates, are given by; k - E h k ^ 2 ( 1 - v 2 ) i j k.. = m.m. F(m. + m. - 2,n. + n.) + q.q. F(p. + p.,q. +q. - 2) ij i j i j ' i j M ^ j K i ' J ' ^ I M j 1-v , n.n. F(m. + m.,n. + n. - 2) + p.p. F(p. + p. - 2,q. + q.) , ^ i j i j ' i j ^ r j " i * j ' M ] ) 1 -v n.p. + vm.q. F(m. + p. - 1 ,n. + q. - 1) j K i ' j M f 1 -v T n i P j + V m i q j J F(m. + p. - 1 ,n . + q. - 1) i j i j (A.2) F. = phfi' C£ F(m.,n.) + F (m. + 1,n.) i = 1, • • • , 1 0 = phfi' CnF(p.#q.} + F(p.,q. + 1) i = 11, • " , 2 0 r = x_ cos 9 + y~ sin 9 c = -x, sin 9 + y, cos £ 3 7 3 n 1 ' 1 (A.3) F(m,n) n+1 I m+1 , . >m+1 c la - (-b) I n l m! n (m+n+2) (A. 4) 103 The stiffness matrix and load vector are transformed to global coordin-ates through the transformation and rotation matrices which are reproduced in Tables A . l and A. 2, respectively. Hence; T -1 T -1 [K] = [R] [T '] [k] [T '] [R] . . . .(A.5) {P} = [R] [T ] {F } . . . .(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 T A B L E A . l T R A N S F O R M A T I O N MATR IX [T] where [ T 2 ] = [ T 3 ] T 1 0 0 T 1 T 2 0 [T] 0 T 3 0 T 4 0 T 2 0 T 3 0 1 -b 0 b 2 0 0 - b 3 0 0 0 | 0 1 0 -2b 0 0 3b 2 0 0 0 0 0 1 0 -b 0 0 b 2 0 0 j 1 a 0 2 a 0 0 3 a 0 0 0 0 1 0 2a 0 0 3a 2 0 0 0 _0 0 1 0 a 0 0 2 a 0 0_ 1 0 c 0 0 2 c 0 0 0 ~3 c 0 1 0 0 c 0 0 0 2 c 0 0 0 1 0 0 2c 0 0 0 n a-b 3 c 3 (a 9 (a-b)c 9 2 c 9 ( a - b ) 3 (a-27 - b ) 2 c 27 (a -b )c ' 27 c 27 105 T A B L E A .2 R O T A T I O N MATR IX [R] where [ R, ] cos 0 0 -sin 0 0 [R] 0 2 cos 6 -sin 8cos 0 -sin 8cos . 2 sin 6 R l 0 0 0 0 R 1 0 0 0 0 R l 0 0 0 0 R, 0 si n 8cos f 2. cos 8 0 . 2 - s in 9 - s in 8cos sin 8 0 cos 8 0 0 sin 8cos . 2 -s in 8 cos -sin bcos . 2 sin 8 sin Bcos 0 sin 8cos ( 2 cos 8 [ R 2 ] cos 8 -sin 6 sin cos 106 TABLE A.3 STRAIN DISTRIBUTION MATRIX [AA] a 2 2a 4 a 5 3a ? 2a g a g a l 3 a15 2 a l 6 a l 8 2 a i 9 ^20 |(a 3 +a 1 2) (a 5+2a 1 4) (2a 6 +a 1 5) (a 8+3a 1 7) (2a g+2a 1 8) ( 3 a 1 0 + a 1 9 ) j TABLE A.4 TRANSFORMATION MATRIX FOR FORCED EXACT STRESS DISTRIBUTION [T] 0 (a-b) 2 a 2 2 0 0 c 2 c 2 a 0 (a-b) : 4 4 bj 4 0 0 0 0 ac 4 - b e 0 0 2 c 0 2 c 4 c2 4 107 A P P E N D I X B M A T R I C E S A S S O C I A T E D WITH T H E SOLUT ION OF THE O U T - O F - P L A N E P R O B L E M The matrices de s c r i bed in Sect ion 3.3 of th i s thes i s a re p r e -sented he re in f u l l fo r the sake of completeness. For a complete de sc r i p t i on of the de r i va t i on of the 18 degree of freedom plate bend i ng element, see References [12,15]. The geo-metry and coord inate systems used for th i s element are shown in F i gu re 3.2. The va r i ab le s at each of the three nodes are The in teger exponents of the d isp lacement f ie ld as de f ined by Equat ion w 3_w 3x 3 w 3y 3 2 w 3x3y ' 3 2 w 3y 2 3.19 a re ; 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 ) The components of the l inear bend i ng s t i f fness matr ix [k ] , the con s i s -tent mass matr ix [m] , and the cons i s tent un i form d i s t r i b u t e d load vecto r {F}, re lat ive to the genera l i zed coord ina tes , a re g i ven b y ; 108 Eh - k.. , J 12(1-v2) 'I k.. = m.m.(m.-l) (m.-l) F(m. + m.-4, n. + n.) IJ 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) 1 J ' J (B.2) m.. = p h F(m. + m., n. + n.) IJ ' J 1 J .(B.3) F. = q F(m., n.) i ^ i ' i 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.. = h T , b, m.m. F(m. 'J k = 1 1 J 1 + m. + m, - 2, n. + n. + n. ) j k ' i j k + 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>k, c k , and d k describe the in-plane stress distribution. The components of the inertia stiffness matrix as derived in Section 3.3.4 are given by 6 k!. = phfi 2 I {b km.m. F(m. + m. + m 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 + 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 . 6 ) where the coefficients b k , c k , and d k are defined by Equation 3.37. In all of these expressions; F(m,n) = c m + 1 , , ,m+1 a - (-b) m! n I n l (m+n+2) (B.7) These matrices and vector are related to the global coordin-ate 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 L ] = [ R ] T [ T ~ 1 ] T [ k L ] [ T " 1 ] [ R ] [ K C ] = [ R ] T [ T _ 1 ] T [ k C ] [ T " 1 ] [ R ] [K 1] = [ R ] T [ T ' 1 ] T [ k ' ] [ T " 1 ] [ R ] [M] = [ R ] T [ T " 1 ] T [ m ] [ T " 1 ] [ R ] {P} = [ R ] T [ T _ 1 ] T {F} . . . . .(B.8) o o o o o o o o o i n 3" r o a- =r o o u ° y o o u u <-> t n o CD .0 (N i n i n 110 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 u o o O O O O O . Q O O O O O f t j O O O O O O . • £ o J3 o o to rp o o j o o o o o o o u (J (D -Q o o o O O O O O O O ( J O (J O O u o o o o o o o o o o o o o o o u o o u o o o o o o o o j a o o o o o n j o o o u o o o o ro IN ro rM O O . D O . Q O O O ( T j O ^ O O O O O O O O o 3- rO <N 5f M (N . Q . Q O . O O O (D B J O ( O O O O O O O O O O o a- CM if CN o o o o o o o o o o o o y o ^ j o o ^ u o o O O O O O . Q O O O O O ( T J O u o o u o o o o o _n o _QOOO rxjo r o o o o o u o o o o _Q _Q o _Q o o nj to l r o KO r o o o o o o o o o o o O O O O O f N O O O O O t N t j O r N O O f N O O O _Q O t— O O O (TJ O »— O O u o o » — o o O _ Q . Q o r M O O to n j O f N o o o r s i o o o o O O * — O O O O O i — O O O U ° O O O O o o o o o r- o o o o o o o o o o o * - o o o o o » — o o o o o o o TABLE B.2 ROTATION MATRIX [R] [R] R 1 0 0 0 R 1 0 0 0 [R,] cost sin 6 -sinB cos( 2 2 cos 6 2sin0cos0 s in ( 2 . 2 -sm0cos0 cos 0-sin 0 sin0cos0 .2 2 sin 0 -2sm0cos0 cos ( 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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

Comment

Related Items