CHARACTERIZATION AND DYNAMIC ANALYSIS OF DAMPING EFFECTS IN COMPOSITE MATERIALS FOR HIGH-SPEED FLYWHEEL APPLICATIONS Except where reference is made to the work of others, the work described in this dissertation is my own or was done in collaboration with my advisory committee. This dissertation does not include proprietary or classified information. Alfonso Moreira Certificate of Approval: Malcolm J. Crocker Distinguished University Professor Mechanical Engineering George T. Flowers, Chair Professor Mechanical Engineering A. Scottedward Hodel Associate Professor Electrical and Computer Engineering Subhash C. Sinha Professor Mechanical Engineering Joe F. Pittman Interim Dean Graduate School CHARACTERIZATION AND DYNAMIC ANALYSIS OF DAMPING EFFECTS IN COMPOSITE MATERIALS FOR HIGH-SPEED FLYWHEEL APPLICATIONS Alfonso Moreira A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Auburn, Alabama May 10, 2007 iii CHARACTERIZATION AND DYNAMIC ANALYSIS OF DAMPING EFFECTS IN COMPOSITE MATERIALS FOR HIGH-SPEED FLYWHEEL APPLICATIONS Alfonso Moreira Permission is granted to Auburn University to make copies of this dissertation at its discretion, upon request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iv VITA Alfonso Moreira, son of R. Alfonso Moreira and Irma Cejudo, was born on May 21, 1975, in Santiago, Chile. He obtained his Bachelor of Science degree in Acoustical Engineering and the title of Acoustical Engineer from Universidad Austral de Chile in November 2000, with his thesis work on acoustics laboratories. He joined the Mechanical Engineering Department at Auburn University as a Research Assistant in January 2002, and became a Doctor of Philosophy candidate in May 2006. His scientific interests cover vibration analysis, nonlinear systems, and acoustics. v DISSERTATION ABSTRACT CHARACTERIZATION AND DYNAMIC ANALYSIS OF DAMPING EFFECTS IN COMPOSITE MATERIALS FOR HIGH-SPEED FLYWHEEL APPLICATIONS Alfonso Moreira Doctor of Philosophy, May 10, 2007 (B.Sc. Universidad Austral de Chile, November 2000) 133 Typed Pages Directed by George T. Flowers The directional mechanical properties of carbon fiber reinforced composite materials make them suitable for components of flywheel energy storage systems. Particularly the hub-rim interface is a component where fiber reinforced composite materials can be applied to reduce rotor mass to achieve high energy densities. However, these materials can introduce significant flexibility and damping into the system, that raise stability issues. This research work consisted of an investigation of the material damping of carbon fiber reinforced epoxy composites and a study of the effect of the material damping on the stability of composite high speed flywheel rotors. In order to characterize the damping of the composite material, a number of beam samples, cut from laminate plates in various configurations, were tested under several boundary conditions. Different methods were used for the extraction of the desired characteristics. The results vi are presented, described and detailed in this dissertation. A prototype of a flywheel rotor was also examined to determine the amount of damping of its composite hub-rim interface and compare these results with the ones of the tests on laminate beams. In addition, a model that captures the main features of flywheel systems was developed, and different configurations were simulated to determine the main factors governing stable ranges of operation. It was observed that some inherent features of flywheel systems allow assumptions that greatly simplify the analysis of the model. Parameter variation studies are presented and discussed in detail. Substantial insight into factors that govern the stability of this kind of high speed rotor system was obtained. vii ACKNOWLEDGMENTS Throughout this research work Dr. George Flowers? experience and vision continuously encouraged me to confront the challenges presented with a practical and open-minded approach. The guidance and advise provided by my other dissertation committee members: Dr. Crocker, Dr. Sinha and Dr. Hodel and the external reader, Dr. Gowayed, were extremely valuable to address the technical content and writing style of this dissertation. The support of my parents Irma and Alfonso and my future wife Tiphaine was crucial to focus on my work and maintain a positive attitude to achieve the completion of this research. viii Style manual or journal used: International Journal of Acoustics and Vibration Computer software used: Microsoft Office 2003, Matlab R2006a, LabView 7.0, Ansys 10 ix TABLE OF CONTENTS LIST OF FIGURES .......................................................................................................... xii LIST OF TABLES.............................................................................................................xv CHAPTER 1 INTRODUCTION .........................................................................................1 CHAPTER 2 BACKGROUND AND LITERATURE REVIEW .......................................3 2.1 Vibration Damping ........................................................................................................3 2.1.1 Damping Models.........................................................................................................5 2.1.2 Measurement of Vibration Damping ..........................................................................6 2.2 Fiber Reinforced Composite Polymers........................................................................12 2.3 Damping in Fiber Reinforced Composite Materials....................................................14 2.4 Modeling of Rotor Systems .........................................................................................16 2.4.1 Rotordynamic analysis..............................................................................................19 2.5 Rotordynamic Instability .............................................................................................20 2.6 Rotordynamic Instability caused by Internal Friction Damping..................................21 2.7 Flywheel as an energy storage system.........................................................................23 CHAPTER 3 DAMPING IN FIBER REINFORCED COMPOSITE MATERIALS .......25 3.1 Experiments .................................................................................................................25 3.2 Beam Supported on Bonded Stud with Random Excitation........................................26 3.3 Excitation at the Center of the Sample ........................................................................31 x 3.3.1 Comparison with Analytical Model..........................................................................32 3.3.2 Modal Damping of Samples Mounted with Stud in the Center................................35 3.4 Cantilever Beams with Swept Sine Excitation of Base ...............................................37 3.4.1 Relation between input acceleration and output displacement.................................41 3.4.2 Discussion of results .................................................................................................43 3.4.3 Finite Element Model of Cantilever Beam Configuration........................................44 3.5 Axially Loaded Beams.................................................................................................46 3.5.1 Observed behavior ....................................................................................................49 3.5.2 The Method of Free Damped Vibrations by Time Blocks .......................................52 3.5.3 Results.......................................................................................................................54 3.6 Natural Frequencies and Damping of a Sample Rotor ................................................56 3.6.1 Measurement Setup...................................................................................................57 3.6.2 Data Acquisition and Analysis..................................................................................58 CHAPTER 4 MODELING AND ANALYSIS OF FLYWHEEL SYSTEMS..................60 4.1 Model of a Flywheel System .......................................................................................61 4.1.1 Model parameters......................................................................................................70 4.1.2 Translational Model of a Flywheel System ..............................................................75 4.1.3 Experiment of Rotor with Flexible Hub-Rim Interface............................................83 CHAPTER 5 CONCLUSIONS .........................................................................................86 REFERENCES ..................................................................................................................92 APPENDIX........................................................................................................................98 Computer Codes.................................................................................................................98 xi Code to generate Figure 3.18 and Figure 3.19...................................................................98 Frequency and amplitude dependence of damping in samples loaded axially ................100 Eigenvalue analysis of matrix A from Eq. (4.33) ............................................................105 Stable thresholds for ranges of values of k H and c H .........................................................107 Stable thresholds for ranges of values of k B and c B .........................................................109 Derivation of Li?nard-Chipart conditions........................................................................111 xii LIST OF FIGURES Figure 2.1. Mass-spring-damper system ...........................................................................4 Figure 2.2. Decay of vibrations of a viscously damped single degree of freedom system .............................................................................................................9 Figure 2.3. Magnitude of frequency response function for a viscously damped system ...........................................................................................................11 Figure 2.4. Randomly oriented units of material (From Hyer [10] )..................................12 Figure 2.5. Poor transverse properties (From Hyer [10] ) ..................................................14 Figure 2.6. Rankine?s model ...........................................................................................17 Figure 3.1. Composite beam with bonded stud for mounting on the shaker...................27 Figure 3.2. Measurement setup .......................................................................................28 Figure 3.3. Damping ratio at first natural frequencies for three fiber alignments...........30 Figure 3.4. Loss factor as function of fiber alignment. (From Suarez et al. [17] ).............30 Figure 3.5. Beam excited at center point.........................................................................31 Figure 3.6. Transfer function of beam attached to the shaker at midpoint......................32 Figure 3.7. Modal damping at four first natural frequencies...........................................36 Figure 3.8. Dog-bone shaped end of the sample .............................................................37 Figure 3.9. Measurement setup .......................................................................................38 Figure 3.10. Experimental transfer function between base and tip of beam.....................39 xiii Figure 3.11. Natural frequency and damping ratio vs. input amplitude acceleration for the three first natural frequencies........................................40 Figure 3.12. Natural frequency and damping ratio vs. displacement of the tip of the beam for the three first natural frequencies ............................................42 Figure 3.13. Modal damping vs. displacement of the beam end.......................................43 Figure 3.14. Mode shapes of vibration at the three first natural frequencies, obtained from finite element model..............................................................45 Figure 3.15. Magnitude of the frequency response obtained from finite element model ............................................................................................................46 Figure 3.16. Test rig ..........................................................................................................47 Figure 3.17. Experimental setup........................................................................................49 Figure 3.18. Two fixed-exponential fittings of the top envelope of free vibration decay of the first natural frequency (593 Hz). Damping ratios differ by 72.8% .......................................................................................................51 Figure 3.19. The same data as Figure 3.18 with y axis shown on a logarithmic scale ..............................................................................................................51 Figure 3.20. Experimental free decay of vibration and the envelopes of the fittings with constant damping ratio and linearly changing damping ratio ...............................................................................................................54 Figure 3.21. Damping ratio vs. vibration amplitude for different 1st natural frequencies....................................................................................................55 xiv Figure 3.22. Damping ratio vs. frequency for various amplitudes of sinusoidal free vibration................................................................................55 Figure 3.23. Measurement setup .......................................................................................57 Figure 4.1. Eight degree of freedom model of a flywheel system ..................................62 Figure 4.2. Relations between coordinate systems..........................................................64 Figure 4.3. Imaginary part of eigenvalues for a range of running speeds.......................73 Figure 4.4. Real part of eigenvalues for a range of running speeds................................73 Figure 4.5. Imaginary parts of eigenvalues for a long symmetric rotor..........................74 Figure 4.6. Translational model for Flywheel with flexible hub-rim interface (Solid Edge drawing by Alex Matras [48] ).....................................................76 Figure 4.7. Maximum stable running speed for hub damping ratio ? H = 0.002 to 0.02 and hub stiffness k H = 0 to 100000 kg/s 2 ..............................................80 Figure 4.8. Out-of-phase mode destabilizes first for k H = 20000 N/m (? = 0.02)...........81 Figure 4.9. In-phase mode destabilizes first for k H = 50000 N/m (? = 0.02) ..................81 Figure 4.10. Maximum stable running speed for bearing damping ratio ? B = 0.015 to 0.1 and bearing stiffness k B = 0 to 100000 kg/s 2 ......................................82 Figure 4.11. Experimental Set-Up.....................................................................................83 Figure 4.12. Rotor rig experiencing unstable behavior.....................................................85 xv LIST OF TABLES Table 3.1. Properties of TORAYCA T300 ........................................................................26 Table 3.2. Relation between first three natural frequencies ..............................................35 Table 3.3. Modal damping at four first natural frequencies ..............................................36 Table 3.4. Natural frequencies and modal damping values...............................................58 Table 4.1. Imaginary parts of eigenvectors (x 10 -3 )...........................................................75 Table 4.2. Model Parameters ............................................................................................ 79 Table 4.3. Average values of experimental data............................................................... 84 xvi NOMENCLATURE X 0 : Vibration amplitude ?: Damping ratio ? : Phase angle ? n : Radian natural frequency ? d : Damped radian natural frequency f n : Natural frequency (Hz) k: Stiffness m: Mass t: time c: Damping coefficient c c : Critical damping coefficient ?: Logarithmic decrement F 0 : Force amplitude G: Frequency response function i: 1? ? 0 : Resonance frequency ? 1 , ? 2 : Half power frequencies xvii T: Transmissibility X b : Displacement of the base X m : Displacement of the mass ?: Rotational speed r: Radial displacement M: Mass imbalance v: Rotating whirl vector ?: Angle between imbalance and rotating whirl vectors V f : Volume fraction of fibers W f : Weight of fibers W m : Weight of matrix ? f : Density of fibers ? m : Density of matrix f i : Natural frequencies Y(x): Beam deflection E: Young?s modulus I: Moment of inertia L: Length of the beam, Lagrangian ? i : Radian natural frequencies a in : Acceleration of base a out : Response point acceleration d out i : Response point displacement of mode i xviii ?: Damping capacity c ? : Bending damping coefficient of shaft-rim interface k ? : Bending stiffness of shaft-rim interface c x : Extensional damping coefficient of shaft-rim interface k x : Extensional stiffness of shaft-rim interface c B? : Equivalent bending damping coefficient of bearings k B? : Equivalent bending stiffness of bearings c Bx : Equivalent extensional damping coefficient of bearings k Bx : Equivalent extensional stiffness of bearings x R , y R : Translational coordinates of the rim x H , y H: Translational coordinates of the hub ? R , ? R : Rotation of the rim about y and x, respectively ? H , ? H : Rotation of the hub about y and x, respectively ??? x yz: Space fixed coordinate system x yz: Body fixed coordinate system x' y' z', x'' y'' z'': Auxiliary coordinate systems ?: Total angular velocity v G : Velocity It: Transversal moment of inertia Ip: Polar moment of inertia T: Kinetic energy V: Potential energy xix x F G : Translational damping force in hub-rim interface F ? G : Rotational damping force in hub-rim interface k x , k H : Translational stiffness coefficient in hub-rim interface k ? : Rotational stiffness coefficient in hub-rim interface c x , c H : Translational damping coefficient in hub-rim interface c ? : Rotational damping coefficient in hub-rim interface rot R x : x coordinate of rim in rotational system st y : y coordinate in the space fixed system Bx F G : Translational damping force in bearings B F ? G : Rotational damping force in bearings c Bx , c B : Translational damping coefficient in bearings c B? : Rotational damping coefficient in bearings k Bx , k B : Translational stiffness coefficient in bearings k B? : Rotational stiffness coefficient in bearings R x? : Virtual displacement of the rim in the x direction H ?? : Virtual rotation of the hub with respect to the y axis q i : Generalized coordinate Q i : Generalized force [M]: Mass matrix [C]: Damping matrix [G]: Gyroscopic matrix xx [K]: Stiffness matrix [I]: Identity matrix ro R : Outer radius of rim ri R : Inner radius of rim ? comp : Volumetric density of carbon-epoxy w R : Width of rim m R : Mass of rim It R : Transversal moment of inertia of rim Ip R : Polar moment of inertia of rim ro H : Outer radius of hub ri H : Inner radius of hub w H : Width of hub ? Al : Volumetric density of aluminum m Hub : Mass of hub It Hub : Transversal moment of inertia of hub Ip Hub : Polar moment of inertia of hub ? B? : Rotational damping ratio in bearings m H : Mass of hub m R : mass of rim m T : Total mass a: Mass ratio 1 CHAPTER 1 INTRODUCTION High-speed flywheel energy storage systems offer the potential for substantially improved energy storage densities as compared to conventional chemical batteries. In the past years, they have been seriously considered for advanced satellite and vehicle applications. A major concern for such components is the energy/weight ratio or energy storage density. The hub-rim interface, which connects the hub mounted on a shaft to a massive rim, is an attractive candidate for reducing rotor mass. The rim is intended to concentrate mass as far from the shaft axis as possible, but the hub-rim interface lies close to the shaft and contributes little to the overall energy storage capacity while adding to the system mass. Some candidate designs use composite materials that can be tailored to withstand the stresses although possessing low mass. In addition, fiber reinforced composite rotors are regarded as safer than metallic rotors, since their failure modes are normally less destructive [1] . However, composite materials allow significant flexibility and tend to have relatively high internal damping, which may produce stability problems. Material damping is, in general, a very complex phenomenon and it is difficult to characterize its properties for a broad range of conditions. There are a variety of methods that are available to evaluate the damping for small amplitude vibrations and relatively narrow frequency ranges. Some methods focus on the natural frequencies of vibration, such as the free damped vibration method and the resonance curve (or half-power 2 bandwidth) method. Other approaches involve the response characteristics at frequencies somewhat removed from resonance, such as the hysteresis method. The former methods provide information that is more directly applicable to assessing stability characteristics and are the methods of choice for this work. From measurements on a set of composite material samples obtained from the Polymer and Fiber Engineering Department, the magnitude of the vibration damping of the material under several conditions was determined and the dominant mechanisms on the dissipation of energy in these composite materials were identified as well. Further, the effects of potential factors affecting the dynamic characteristics of the material were ascertained within ranges determined by applications on which the use of the composite material would offer significant benefits. A tentative design of a flywheel energy storage system was considered as a baseline to establish these factors. In addition, a model for analysis of the dynamics and stability of a flywheel system was developed and a series of parameter variation studies are discussed here in detail. A number of useful conclusions and insights for design of such systems were obtained and are presented. 3 CHAPTER 2 BACKGROUND AND LITERATURE REVIEW 2.1 Vibration Damping Along with mass and stiffness, damping determines the essential dynamic characteristics of a structure. While mass and stiffness are associated with energy storage, damping relates to the conversion of mechanical energy into other forms of energy, such as heat or sound. Damping in general affects only vibrational motions around the resonance frequencies of a system. If a classical mass-spring-damper system is considered (See Figure 2.1), for excitation frequencies that are considerably lower than the natural frequency of the system, the motion is mainly determined by the spring force, and is known as stiffness controlled. If, on the other hand, the frequency of the excitation force is considerably above the natural frequency of the system the inertia of the mass will have a greater effect on the response. This region is usually called mass controlled. However if the excitation frequency matches the natural frequency, that is, at resonance, the spring and inertia effects cancel each other. The excitation force provides energy to the system. The energy increases until a steady state is reached, in which the energy supplied per cycle is equal to the energy lost per cycle due to damping [2] . As a result of an increase of vibration damping in a system one finds that unforced and transient vibrations decay faster, and amplitudes of vibration of structures at resonances are reduced. However, some damping mechanisms can be detrimental to the 4 performance of a system. Damping forces present in moving parts, or frictional forces produced between moving parts, can generate self excited vibrations. Figure 2.1. Mass-spring-damper system A considerable amount of literature is available on the subject of vibration damping and in particular on material damping. Almost every modern vibration book has a section dedicated to it. Linacre [3] [4] in his publications on Iron & Steel (1950) provides some of the earliest reviews on damping research to date. Crandall [5] investigated the nature of damping, pointing out some amplitude and frequency dependence of damping and the limitations of some idealized models. Some authors such as Lazan [6] investigated the characteristics of vibration damping in more depth. Lazan?s text contains not only a thorough description of the most common models for damping characterization, but also a comprehensive compendium of levels of damping for different materials, specifying in most cases the testing conditions used such as vibration modes used (torsion, axial, bending, etc.) and environmental conditions. Another important text dedicated to the c m k 5 topic, although centered on vibration control by means of viscoelastic materials, is the one by Nashif et al. [7] . 2.1.1 Damping Models Damping is present in systems from several different disciplines, so there is a variety of damping mechanisms as well as approaches to interpret and describe them. Three major models are used to describe damping in mechanical vibrations: Coulomb, viscous and hysteretic damping. Each of these models describes a different phenomenon producing dissipation of vibration energy. Coulomb damping is caused by kinetic friction between sliding dry surfaces. Viscous damping is a form of fluid damping in which the damping force is proportional to velocity. Hysteretic damping, also referred to as solid damping, is caused by the internal friction or hysteresis when a solid is deformed [8] . Viscous damping is the most common of these three mechanisms. Strictly speaking, viscous damping only describes damping produced by laminar flow or by fluid passing through a slot, as in a shock absorber [8] , but it is frequently used to describe other types of energy dissipation without incurring great errors, when the dissipative forces are small. For the specific case of internal friction, the theory of elastic hysteresis is the most widely accepted. This model is based on the fact that the relation between stress and strain is nonlinear and different for the loading and unloading. However, a few more detailed theories have been developed that provide other explanations of the phenomenon of vibration damping and more detailed or versatile models that in turn add complexity to the analysis. Of these, the most relevant ones are the theory of linear hereditary elasticity or viscoelasticity, based on integral relations between stresses and strains, the theory of 6 microplastic deformation where dissipation is attributed to motion of dislocations in micro volumes, and Zener?s thermodynamic theory, in which dissipation is considered to be a consequence of the heat fluxes between parts with different stresses [9] . 2.1.2 Measurement of Vibration Damping The methods used for experimental investigation of energy dissipation response are classified into two groups. The first group consists of the so called direct methods, based on direct measurements of energy dissipation. The second group is the indirect methods, in which changes in other parameters such as amplitude and frequency are related to the amount of energy dissipation [9] . The energy method is a direct method in which the electrical or mechanical excitation required to maintain steady-state vibrations in a sample provide a direct measure of the energy being dissipated. The thermal method is a direct method that relies on the hypothesis that the majority of the energy dissipated is transformed into heat, and thus it uses a measure of the heat generated by the vibrational motion as a direct measure of the energy dissipated. It is apparent that the difficulties encountered when trying to accurately quantify the heat generated from the vibration process make this method hard to apply. Moreover, heat is not the only mechanism of energy dissipation, since irreversible changes in the structure of the material such as dislocation movements and cracks growth also take part of the effective vibrating energy [9] . The method of the hysteresis loop is perhaps the most popular among the direct methods. It uses the area of the hysteresis loop formed by the stress-strain curve during cyclic loading and unloading of the sample as a measure of the energy loss. Since the 7 relative energy dissipation for materials is very small, the area enclosed in the hysteresis curve is very small, so a high accuracy is required for the measurement of the strain [9] . Another difficulty found in the application of this method is that it needs a very precise tracking of the phase of each, the stress and strain measurements, a factor that becomes more significant as the excitation frequency is increased. The indirect methods include the method of free damped vibrations and the resonance curve or half-power bandwidth method. To explain the former let us consider the simple ideal linear mass-spring-damper system shown in Figure 2.1. The forces involved in the motion of this system are: -kx, produced by the spring; and - cx , produced by the damper and the equation of motion, when no external excitation is applied, is 0.mx cx kx+ +=  (2.1) If the system is released from a position X 0 respect to its equilibrium position, the displacement in the following instants follows the expression - 0 () cos( ) n t d xt X e t ?? ? ?=+, (2.2) for ? < 1, where ? is referred to as damping ratio [2] , and is defined as c c c ? = , (2.3) and c c is known as the critical damping coefficient and is defined by 2 2 cn ckmm?==. (2.4) ? is a phase angle that depends on the initial velocity, and ? n and ? d represent the undamped and damped radian natural frequencies of the system, respectively. They relate to the other parameters by 8 2 nn k f m ? ?== , and (2.5) 2 1 - dn ? ? ?= , (2.6) and f n is called the (undamped) natural frequency. The damping ratio ? usually has a very small value for structural materials, which means that ? d and ? n are sufficiently close to each other to allow the approximation ? d = ? n . The right hand side of Eq. (2.2) contains a cosine function with amplitude X 0 e -?? n t that decreases with a rate of ? ? n as time t increases. The time trace representing this free decay of the oscillations of the system after an excitation has ceased, provides a clear graphical way of seeing the effect of damping, as shown in Figure 2.2. The method of free damped vibrations uses this trace to obtain a measure of the energy dissipation from the decay in the amplitude of the vibration on one or more cycles of vibration. A common measure used in this method is the logarithmic decrement, ?. For n cycles of the free vibration decay, it is defined as: 1 ln i in X nX ? + ?? = ?? ?? , (2.7) where X i and X i+n represent the values of x at two peaks separated by n cycles. Under the assumption that a system is viscously damped and ? <<1, it follows from Eq. (2.2) that ? = 2 ? ? , or 1 ln 2 i in X nX ? ? + ?? = ?? ?? . (2.8) 9 Figure 2.2. Decay of vibrations of a viscously damped single degree of freedom system The second indirect method to consider is the resonance curve or half-power bandwidth method. If a vertical sinusoidal force defined by F(t) = F 0 cos(?t) is applied to the mass of the system described above, then the motion of the mass after transients have vanished, is also sinusoidal. The ratio between the resulting displacement and the force applied is called the frequency response function, G, and for this system in particular it has the form: 2 0 () 1 () 1 - ( / ) 2 / nn Xi Gi Fi ? ? ? ??? == + . (2.9) The magnitude of this complex expression is the real expression [] 2 2 2 1 () 1 - ( / ) 2 / nn Gi? ?? ??? = ??+ ?? . (2.10) It is possible to obtain experimentally a curve for the frequency response versus excitation frequency. Using Eq. (2.10) it can be shown that for ? <<1, the damping ratio is given by D i s p l a c e m ent X i X i+4 X 0 ? n t 2? 4? 6? 8? 10? 12? 14? ? n t X 0 10 21 0 - 2 ? ? ? ? = , (2.11) where ? 0 is the frequency at which the peak of the curve is obtained and ? 1 and ? 2 are the two frequencies, one below and one above, for which the frequency response is () 1 2 ? times the one at resonance. These frequencies are often called the half-power points because at these the energy stored in the system (and that dissipated by it), which is proportional to the square of the amplitude, is half of the maximum value. Figure 2.3 shows the magnitude of the frequency response of a viscously damped system for several values of damping ratio ?. Another expression can be used to identify the value of damping for a single degree of freedom system in a similar way. If on the mass-spring-system of Figure 2.1 the force is applied at the base instead, then the sinusoidal motion of the base produces a corresponding motion of the mass. The transmissibility of a system is a measure of how much the motion of the base or foundation influences the motion of the mass for a range of frequencies. For a system with viscous damping under sinusoidal excitation, the transmissibility, T, follows [] [] 2 2 2 2 1 2 / () 1- ( / ) 2 / n m b nn X T X ??? ? ?? ??? + == ??+ ?? , (2.12) where X b represents the displacement of the base and X m the displacement of the mass 11 Figure 2.3. Magnitude of frequency response function for a viscously damped system As can be seen, this expression is slightly different from that of Eq. (2.10). However for ?<<1 neglecting the term [2 ? ?/? n ] 2 on the top introduces little error on the estimation of the transmissibility, and thus the expression can be approximated as [] 2 2 2 1 () 1- ( / ) 2 / m b nn X T X ? ?? ??? =? ??+ ?? . (2.13) From this expression the damping ratio can be related to the displacement of the mass and that of the point where the force is applied. The advantage of this expression over the one for the magnitude of the frequency response is that the displacement of the point where the force is applied can be tracked with non-contact techniques that don?t interfere with the system. ? 1 ? 20.4 0.6 0.8 1 1.2 1.4 1.6 0 1 2 3 4 5 6 7 8 ? /? n |G (i?)| ? = 0.2 ? = 0.14 ? = 0.09 ? = 0.07 ? = 0.05 12 2.2 Fiber Reinforced Composite Polymers The search for lighter materials that allow further tailoring of the designs of structural components has found an answer in common materials such as iron, copper, nickel, carbon, and boron. To varying degrees, these materials have directionally dependent mechanical properties, with the directional dependence being due to the strength of the interatomic and intermolecular bonds [10] . Some directions exhibit stronger bonds than others and a material unit (which can range from the molecular to the macroscopic level) in which these bonds are aligned in certain directions is very stiff and considerably stronger in those directions. However, in the other directions usually the material is much softer and weaker. If a material is fabricated in bulk form, it will contain randomly oriented units of material, as shown in Figure 2.4, and the bulk material will have the same mechanical properties in all directions. These properties will reflect in general the properties of the weakest link of the unit. Figure 2.4. Randomly oriented units of material (From Hyer [10] ) 13 If, on the other hand, a material is processed in a manner that permits the alignment of the strong and stiff directions of all the basic units, some of the high strength and stiffness properties of all the basic material units can be preserved, along selected directions. Long and thin elements of material referred to as whiskers where units are aligned can be formed. Their mechanical properties can be close to those of a single unit if enough care is taken in processing. However, the process of enlarging a whisker by adding more basic units inevitably causes imperfections that significantly affect the strength and stiffness of the whisker and become the weak link in the material. Nevertheless, the units formed by adding to the length of whiskers, called fibers, have significant lengths, so they can be easily aligned in one direction to provide directional reinforcement to a structure. At the same time fibers can be aligned and grouped in what is called a tow, which further improves the handling of the fibers, especially when their diameters are small as is the case of most forms of carbon fibers. Fiber tows are embedded and bonded to another material in order to make use of them. This material is often called the matrix, and is usually softer and weaker than the reinforcement material [10] . A fiber reinforced composite material is formed by the embedding of a parallel array of strong, stiff fibers or tows in a matrix. Loads applied along the direction of the fibers will be transmitted to the fibers, which will assume most of the resistance to the load, as in Figure 2.5. However, if the load is applied perpendicular to the alignment of the fibers as in Figure 2.5b-c, a great part of it will be entirely in the matrix material. 14 Figure 2.5. Poor transverse properties (From Hyer [10] ) 2.3 Damping in Fiber Reinforced Composite Materials Conventional metallic materials exhibit very low values of damping. It is customary to assume that most of the energy dissipation in metallic structures occurs at the joints or in added damping treatments. Polymer composites, on the other hand, exhibit large values of damping. This has often been regarded as a positive characteristic, since damping is desirable for many applications where persistent oscillations are detrimental to performance. However, there are applications where excessive damping can cause severe problems, and thus a proper characterization of their dynamic behavior becomes critical to generate optimal designs. A considerable amount of work has been done in the field of dynamic characterization of composite materials. A comprehensive review of the research in this area is given in two publications by Gibson [12] [13] . Bert [14] also reviewed the early contributions to the field of dynamic behavior of composite materials and structures, a (a) (b) (c) Fiber direction Transverse directions 15 more experimental approach covering dynamic stiffness and damping, vibration of structural elements, and low-velocity transverse impact of laminated panels. The textbook by Zinoviev and Ermakov [9] covers the basics of damping analysis in composites and provides some measurement data. Bert [15] reviewed the theory of damping in fiber-reinforced composites for perfectly-bonded viscoelastic composites. Chaturvedi [16] provided an overview of the analytical and experimental characterization of damping in polymer composites for discontinuous and continuous fiber reinforcements. Suarez et. al. [17] investigated the influence of fiber length and fiber orientation on damping of polymer composite materials. Chia [18] published a review of the geometrical nonlinear static and dynamic behavior of composite laminates. Plunkett [19] reviewed the damping mechanisms believed to be present in fiber composite laminates. Yen and Cunningham [20] studied the effect of anisotropy in mode shapes and frequency distribution on graphite-epoxy plates, finding that the behavior is quite different to that of isotropic plates. Damping in composite materials is attributed to a number of sources, namely: a) The viscoelastic nature of the matrix and/or fiber materials. In composites with a polymeric matrix this effect is more pronounced [21] . b) Thermoelastic damping due to heat flow. It is assumed that the heat flows between areas at different stress states and consequently at different temperatures [9] . c) Coulomb friction generated from the slip in the matrix/fiber interface at unbonded regions. 16 d) Energy dissipation at cracks and delaminations, also related to Coulomb damping produced at damaged locations [21] . e) Viscoplastic damping, non-linear damping at large amplitudes of vibration, due to high levels of stress and strain. Adams and Maheri [22] have determined that the non-linearity in damping can be attributed to plastic deformation beyond certain critical stress level. Kenny and Marchetti achieved to correlate the load level, the high damping of plastic origin, and its thermal effects for carbon and graphite fiber reinforced polymers [23] . f) Hwang [24] concluded that the effects of transverse shear on the damping of laminated beams in flexural vibration and of interlaminar stresses on the damping of laminates under extensional vibration are most important in thick laminates. The data available for damping in polymer composite materials is very dissimilar [6] . The types of matrix and fiber materials, fiber length, curing temperature, laminate configuration, etc. are all factors that can greatly affect the energy dissipation properties of the material. Values for damping are often found in literature with poor reference to the method used for measurement, the environmental conditions, and the characteristics of the material selected. This complicates the task of comparing and validating experimental results. 2.4 Modeling of Rotor Systems The study of rotating structures is a field that has developed more as an experimental science than a theoretical one. The first analysis of a spinning shaft was 17 presented by W. J. Rankine in 1869. Rankine chose a model, shown in Figure 2.6, to examine the equilibrium conditions of a frictionless, uniform shaft disturbed from its initial position. In his analysis he neglected the Coriolis acceleration in the second equation of motion [25] and thus predicted incorrectly that rotating machines were not able to exceed their critical speed. Figure 2.6. Rankine?s model Rankine?s assertion was contradicted by contemporaries such as Foppl, whose demonstration of the existence of a stable supercritical running speed was not widely recognized, and De Laval, who in 1889 was able to run a single stage steam turbine at a supercritical speed. It was after almost 50 years that Henry H. Jeffcott performed the task of clarifying the issue and satisfactorily explained the phenomenon using a model that consists of a massive unbalanced disk mounted half way between rigid bearing supports on a flexible shaft of negligible mass, and where viscous damping opposes absolute Y m m? 2 r ? r kr X 18 motions of the disk. He was able to explain how the rotor whirl amplitude is maximized at the critical speed, ? = ? c , but diminishes as ? > ? c [26] . Further details can found in the article by Nelson [27] . Figure 2.7 shows the Jeffcott rotor model in whirling motion. The shaded square M represents an unbalanced mass. The whirl speed, ? ?=  , is the time rate of change of the angle ? . If the angle ? remains constant relative to the rotating whirl vector v, the whirl speed and the shaft speed are the same, thus the whirl is called synchronous. If, on the other hand, the angle ? has a rate of change 0? ?  , the whirling motion is referred to as non-synchronous. Figure 2.7. The Jeffcott rotor in whirling motion ? X Y ? ? M v 19 2.4.1 Rotordynamic analysis Modern high speed rotating machines are encountered in several applications where extreme production or storage of energy is desired. Their ability to achieve high shaft speeds allows them to deliver high energy densities and flow rates. This comes at the expense of high inertial loads and potential problems like vibration, shaft whirl and rotordynamic instability [26] . Rotordynamic analysis deals with the planning, design and adjustments to the designs of rotating machinery. Some of its main objectives are [26] : ? Predicting critical speeds, defined as the angular rates ? at which vibration due to imbalance of the rotor (the assembly of rotating parts) is a maximum. These speeds can be calculated from design data so that they are avoided when setting operational speeds. Rotordynamic analysis also offers methods to evaluate how modifications of the parameters will affect a design when critical speeds must be distanced from a given operational speed. ? Calculate the locations and masses adequate to achieve balancing of rotors, in order to reduce the amplitude of synchronous vibration. ? Predict threshold speeds at which dynamic instability occurs and determining suitable modifications in the design so as to suppress dynamic instabilities. This can be a challenging task, since destabilizing forces are hard to identify qualitatively and quantitatively, and thus it becomes difficult to represent them accurately in mathematical models. 20 2.5 Rotordynamic Instability Within the problems found in rotating machinery, synchronous whirl (produced by imbalance) is the most common. However, although nonsynchronous whirl is less frequent, it can severely damage a machine. Within these nonsynchronous phenomena is the rotor whirling that becomes unstable when a certain speed (called the threshold speed of instability) is reached, which has proven to have devastating effects on rotor systems. It is produced by tangential forces that act in the direction of the instantaneous motion. They are usually referred to as following or destabilizing forces and its magnitude can be proportional to the whirl velocity, in which case they are considered as negative damping, or proportional to rotor displacement, classified as a cross-coupled stiffness. Several mechanisms have been identified or at least are believed to produce rotordynamic instability. Oil whip is probably the most common source of instability in hydrodynamic bearings. It occurs when the shaft in the bearing is disturbed from equilibrium and the oil film starts to drive it in a whirling motion. This can occur until a point when the oil frequency matches a natural frequency of the system and remain unchanged as the running speed continues to increase. This is the phenomenon known as oil whip, which may cause destructive vibration [28] . Other less common sources of rotor instability are fluid ring seals, similar in nature to oil whip; internal friction in or between rotating parts; Alford?s forces, produced by irregular circumferential blade-tip clearances in an eccentric rotor [29] ; trapped liquids inside a hollow shaft or rotor; and dry friction whip, produced by rubbing friction between the rotor and stator, which originates a backward whirl motion [26] . Of all these, rotor instability caused by internal friction is the central interest of this work. 21 2.6 Rotordynamic Instability caused by Internal Friction Damping As early as 1924, observations of rotor instability were reported by Newkirk [30] , a phenomenon he referred to as ?whip?. In order to determine if unbalance was the cause for the observed phenomenon, he conducted a study using a test rotor to simulate a compressor unit, and drew a series of important conclusions [32] : ? Refinement in rotor balance does not affect the onset speed of whirling or whirl amplitude. ? Whirling always occurrs above the first critical speed. ? The whirl speed is constant regardless of the rotational speed. ? Misalignment of the bearings increases stability. ? Introducing damping into the foundation increases the whirl threshold speed. ? In a well balanced rotor, a disturbance is sometimes required to initiate the whirl motion. Newkirk realized that this phenomenon could not be attributed to critical-speed resonance, since the high vibrations encountered always occurred super-critically, i.e. above the first critical speed, and refinement of balance had no effect upon diminishing the whirl amplitudes. It was Kimball [31] (1924) who suggested that internal shaft friction can be responsible for shaft whirling. He postulated that below the rotor critical speed the internal friction damps out the whirling motion, while above the critical speed the internal friction sustains the whirl [32] . He attributed this effect to the hysteresis of the metal undergoing alternate stress reversal cycles. This led Newkirk to extend the Jeffcott model by adding a force normal to the deflected rotor, with which he could demonstrate that the rotor is unstable above the first critical speed. However, since he did not include the 22 effect of the flexibility and damping of the supports, he could not explain theoretically several key points of his experimental observations. Research work in the area developed with the years and some landmark papers and texts were published that covered the topic of instability produced by internal friction damping and other mechanisms. Ehrich [33] (1964) was able to determine that the ?consideration of the stabilizing effects of external friction leads to the more general conclusion that shaft whirl may occur at any natural mode?. He established that the rotational speed at which instability occurs is governed by the ratio of external friction to internal friction. By modeling a flexible rotor on elastic supports, Gunter [32] (1967) was able to come up with an analytical expression to predict the onset speed of instability and provided a theoretical explanation to Newkirk?s findings. He also proved that the threshold speed of whirl instability can be increased by decreasing the foundation stiffness. Then in 1969 he and Trumpler [34] showed that in the absence of bearing damping a symmetric flexible foundation reduces the rotor critical speed and also the whirl threshold speed. They also concluded that addition of internal damping greatly improves the threshold speed. They extended the investigation to consider an asymmetric foundation finding that, even with no damping added, the onset speed of instability is largely increased. Lund is widely recognized for his fundamental contributions to rotordynamic analysis, and this is also the case with internal friction instability. In his paper [35] (1974), he extended the Myklestad-Prohl method for calculating critical speeds, to calculate the damped natural frequencies of a general flexible rotor supported in fluid film journal 23 bearings. His method is significant because of its versatility to simulate virtually any practical shaft geometry and support configuration. Bently and Muszynska [36] (1982) determined that the effective rotor damping was reduced due to internal damping during sub-synchronous and backward precessional vibrations produced by other instability mechanisms, and verified that internal damping is indeed a source of rotor instability. Some authors have treated the topic of composite materials used for rotor systems. The work of Wettergren [37] (1998), dealt with the characterization of high- modulus carbon fibers in an epoxy matrix, to be used in shafts. Previously, the work by Chen [38] (1978) consisted in modeling an overhung flywheel rotor system with a flexible shaft, in which the rim was attached to the hub by elastic bands of unidirectional Aramid- Epoxy. This work was valuable in establishing some analytical tools for analyzing a flywheel with flexible hub-rim interface, but it did not address the characterization of the level of damping of the composite material used and its direct effect on the system?s stability. 2.7 Flywheel as an energy storage system Flywheels as energy storage systems have a long history. However only in the past decades they have been considered for more serious applications, and thus further research has been put into developing more efficient designs [39] . It was in the early 70?s that the idea of using reinforced plastics as a way to increase the energy/weight ratio started to be developed. Recently, interest has been shown in incorporating composite flywheels in aerospace applications as energy storage and combined systems for energy 24 storage and attitude control. This has generated a series of research efforts by governmental agencies together with the academic community, to create safe and reliable flywheel systems. Gowayed et al. [41] (2002) established some criteria for the optimal design of composite flywheel rotors, using both closed nonlinear and finite element analysis optimization. They maximized the total energy of the rotor as a function of geometrical and physical characteristics of the composite rim and the rotational speed. They also analyzed the potential of using closed form analyses to give initial estimates of optimal designs, and finite element analysis for more accuracy and a better insight on manufacturing approaches. Jansen et al. [42] (2002) described some changes in the design of the flywheel module at NASA Glenn Research Center. They incorporated a composite rim and magnetic bearings, among other improvements. They were able to meet the safety margins at the certification speed of 66,000 RPM. 25 CHAPTER 3 DAMPING IN FIBER REINFORCED COMPOSITE MATERIALS 3.1 Experiments A set of carbon fiber reinforced composite polymer plates was prepared by the Polymer and Fiber Engineering Department at Auburn University, in order to characterize their damping characteristics. These plates were all fabricated using the prepreg method, in which several sheets containing aligned carbon fibers are bonded using an epoxy matrix in a high temperature press. Different numbers of layers and relative alignments of the sheets were used to span a variety of configurations. In order to determine the extent to which the mounting conditions interfere with the proper determination of the material damping, a series of experiments were conducted using different mounting conditions. Some of the configurations aimed at characterizing the material itself, and some were designed in order to include the boundary effects that would be present in a real application. This would allow separating the contributions of the material, the configuration, and the type of mounting to the damping, providing additional information useful for the application of the results in the design of more complex structures. The samples were prepared from Carbon-Epoxy prepreg sheets using high- strength carbon fibers of the type TORAYCA ? T300. This kind of fiber has 3000 26 monofilaments in a tow, where each filament has an approximate diameter of 7 ?m. Some characteristics of this kind of material are given in Table 3.1. Tensile Strength Tensile Modulus Elongation Density Fiber Type ksi MPa Msi GPa % g/cm 3 T300 512 3530 33.4 230 1.5 1.76 Table 3.1. Properties of TORAYCA T300 Each sheet has a thickness of 0.12 mm and the volume fraction of the final plates is 62%. The fiber volume fraction, V f , can be obtained using the equation mf f f mmf W V WW ? ?? = + , (3.1) where W f = weight fraction of fibers, W m = weight fraction of matrix, ? f = density of fibers [g/cm 3 ], and ? m = density of matrix [g/cm 3 ]. 3.2 Beam Supported on Bonded Stud with Random Excitation First, samples supported with a bonded stud were tested. For this test the plate of composite material that was used had an alignment configuration [0?,0?,0?,90?,0?,0?,0?] (seven layers), which is very close to a unidirectional laminate plate. It had a thickness of 1.08 mm and an area of about 300 x 300 mm. The layers of the plate showed mismatched 27 edges, causing the epoxy bonding at the borders to be uneven. From borders in such condition, delamination can develop at the edges and spread towards the center of the plate. In order to prevent such delamination, stripes of about 25 mm wide were cut around the border of the plate. Special care was taken to maintain the alignment of the main axes of the plate with the fiber directions. Then the shapes of samples with different fiber alignments were drawn on the plate with some added margin for each sample, to allow polishing of the edges to get straight samples. The plate was cut using a ceramic tile saw and the edges of the beams were smoothed and polished using a file and sand paper. Figure 3.1. Composite beam with bonded stud for mounting on the shaker The dimensions of the beam samples were chosen to maximize the use of the area of the plate that was in good conditions and to have a width significantly smaller than the length, so that the assumptions of the Euler-Bernoulli beam theory could be satisfied for analysis. The test samples that were used consisted of beams with a width of 12 mm and lengths ranging from 120 to 250 mm. A threaded stud was bonded to one side of each 12 mm 1.08 mm 28 sample, close to the edge, and then this stud was fastened to the head of the shaker. This system provides the characteristics of a cantilever clamping, but the friction that is produced at the fixed end is significantly reduced. This allows a better isolation of the internal damping of the sample from the damping provided by the friction in the border. Figure 3.2. Measurement setup The damping value at the first natural frequency of each sample is quantified using an equivalent viscous damping ratio which is obtained from the transfer function magnitude plots using the half-power bandwidth method. Vibration signals were measured at the point where the sample is connected to the head of the shaker (as the input) and at a point near the end of the sample (as the output), where the highest amplitudes of the first mode of vibration are achieved. Samples were excited with a random noise signal of limited bandwidth around the center frequency of vibration of the first natural mode and the responses between the input and output were averaged over multiple frames. The first natural mode of vibration and the damping ratio were determined from each resulting Bode diagram. Head of Shaker Measurement Points Sample 29 The procedure was repeated using beams of different fiber alignment and different lengths to obtain the values of damping ratio for a range of natural frequencies. The fiber alignments tested were 0?, 90? and 45?, where each alignment represents the angle between the the 6 laminae with 0? and the longitudinal axis of the beam. The results of these measurements are shown in Figure 3.3. Examination of this figure provides some very interesting insights. First, it is very important to note that the damping ratios for each of the three sample fiber directions are approximately constant in the range of frequency studied. From a modeling perspective, this result indicates a linear (viscous type) characteristic over the frequency ranges tested, which serves to greatly simplify the basic analyses. Also, as expected, the damping levels change dramatically as a function of fiber alignment. The lowest values were observed for the 0? configuration (where all layers are aligned at 0? except the center one) at about 0.2%. Somewhat higher values were seen for the 90? configuration at about 0.25%. Substantially higher values were noted for the 45? configuration, at about 0.4%. This is in agreement with the results of similar studies that assessed the damping as a function of fiber alignment. Figure 3.4 shows the result of a study by Suarez et al. [17] where damping (represented by the loss factor) is shown as a function of fiber alignment. Since flexible hub designs will probably be constructed by winding the prepreg material around a mold, the alignment angles will vary considerably, but damping values will certainly fall within the ranges obtained. 30 20 30 40 50 60 70 80 90 100 110 120 130 0 0.002 0.004 0.006 0.008 0.01 1st natural frequency (Hz) dam pi ng r a t i o ? sample: 0deg sample: 90deg sample: 45deg Figure 3.3. Damping ratio at first natural frequencies for three fiber alignments Figure 3.4. Loss factor as function of fiber alignment. (From Suarez et al. [17] ) GRAPHITE EPOXY SPECIMENS Continuous fiber Fiber volume fraction 0.675 Fiber loss factor 0.0015 LEGEND ? Corrected force-balance approach ? Mean experimental value I Experimental scatter Fiber direction (degrees) 0.01 0.02 0.03 0.04 0.05 0.06 0 0 15 30 45 60 75 90 31 3.3 Excitation at the Center of the Sample An attempt was made to use a configuration in which the excitation was applied at the center of the beam. This connection was again made with a threaded stud connected to the sample by means of high strength epoxy. The intention was to achieve a system that would behave as a free-free beam, at least for the odd modes of vibration, since the connection point would be located at a node of those vibration modes. Figure 3.5. Beam excited at center point The sample used for this set of measurements came from a plate with 22 layers in a [0?,90?,0?,90?,0?,90?,0?,90?,0?,90?,0?] S configuration (where S stands for symmetric) that better simulates the conditions found on a component for a real application, as compared to the samples used to obtain the dependence of damping on fiber alignment, in the previous section. This sample plate has a thickness of 3.16 mm and the sample beam has 281 mm of length and 15 mm of width. The sample was cut using the same provisions as described above to avoid the delaminated sections that are present on the edges of the plate. A broadband frequency transfer function was obtained between measurements of random noise using accelerometers mounted at the shaker?s head and at a point close to Head of shaker Threaded Stud Sample 32 the tip of one of the sides of the beam in bending, which resulted in the response shown in Figure 3.6. The accelerometers used were miniature accelerometers of approximately 1 g of mass, including the effect of the attached cable. Figure 3.6. Transfer function of beam attached to the shaker at midpoint The natural frequencies for this system are identified to be at f 1 = 181 Hz, f 2 = 1078.75 = 5.959 (181) Hz, (3.2) f 3 = 2937.75 = 16.231 (181) Hz, 3.3.1 Comparison with Analytical Model Let us consider the Euler-Bernoulli beam equation of motion 4 4 4 () - ( ) 0 dY x Yx dx ? = , (3.3) where 2 4 n m EI ? ? = . (3.4) 33 Y(x) is the deflection of the beam at x. The form for the solution of this kind of equation is known to be: ( ) sin( ) cos( ) sinh( ) cosh( )Yx A x B x C x D x? ???=++ + (3.5) The boundary conditions for a beam of length L in a free-free configuration are 2 2 0 () 0 x dYx dx = = , 2 2 () 0 xL dYx dx = = , (3.6) 3 3 0 () 0 x dY x dx = = , and 3 3 () 0 xL dY x dx = = . If the second and third derivatives of the assumed solution are evaluated at the boundaries and the conditions given in (3.6) are applied, it is concluded that for a free-free configuration the following equation defines the natural frequencies of vibration: cosh( ) cos( ) 1LL? ? = . (3.7) There are an infinite number of values of ?L that satisfy this equation, the first ones being: ?L = (0, 4.73, 7.853, 10.996, 14.137, 17.279, ?) (3.8) and from Eq. (3.4), the natural frequencies for a free-free beam are 4 (0, 22.373, 61.67, 120.903, 199.86, 298.56, ...) n EI mL ? = , (3.9) or 4 (0, 1, 2.75, 5.404, 8.933, 13.34, ...) 22.373 n EI mL ? =??, (3.10) where ? n = 0 corresponds to a rigid body displacement. The remaining values of ? n are the predicted natural frequencies for a free-free beam. Since in this case the beam is excited in the center, the odd modes, in which the halves of the beam oscillate out of 34 phase, are not present (except ? 0 = 0). This occurs because the center point is a node for these vibration modes. Thus, the spacing between the first three (even) natural frequencies for this situation should be given by the even indexed values in Eq. (3.10), i.e. ? 2 = 5.404 ? 1 , (3.11) ? 3 = 13.34 ? 1 , etc. This theoretical relation between the natural frequencies for the free-free beam does not match the results of the experiment very closely. A possible cause is the type of connection between the stud and the beam. Since the system was being mechanically excited, a solid connection between the stud and beam was necessary and, consequently, the area of contact could not be kept too small. This means that the connection was not on a ?point? but rather a ?small area? at the center of the beam, so points around the center of the beam were constrained and could not deflect freely. This would imply that the system would have features of a double cantilever beam instead. In order to investigate further, an analysis similar to that done above for a free-free beam was performed. The boundary conditions for a cantilever beam are 0 () 0 x Yx = = , 0 () 0 x dY x dx = = , (3.12) 2 2 () 0 xL dYx dx = = , and 3 3 x = L () = 0 dY x dx . Evaluating the first, second, and third derivatives of (3.5) in the boundaries and using these conditions, the equation that defines the natural frequencies of vibration of a cantilever beam is obtained: 35 cosh( ) cos( ) -1LL? ? = . (3.13) The values of ?L that satisfy this equation are ?L = (1.8751, 4.694, 7.85, 10.99, 14.137, 17.28, ?), (3.14) where the rigid body mode of vibration, ? 0 = 0 on Eq. (3.8), is not present as expected. Again, using Eq. (3.4) for a cantilever beam, the natural frequencies are 4 (3.5160, 22.0336, 61.6225, 120.7801, 199.8548, 298.5984, ...) n EI mL ? = , (3.15) or 4 (1, 6.27, 17.53, 34.35, 56.84, 84.93, ...) 3.516 n EI mL ? =??, (3.16) and the spacing between the first three natural frequencies for this situation is given by ? 2 = 6.27 ? 1 , and (3.17) ? 3 = 17.53 ? 1 . A comparison of these results can be observed in Table 3.2. Experimental Free-Free B. C. Cantilever B. C. ? 2 = ? 1 ? 5.959 5.404 6.27 ? 3 = ? 1 ? 16.231 13.34 17.53 Table 3.2. Ratios between first three natural frequencies 3.3.2 Modal Damping of Samples Mounted with Stud in the Center The same procedure applied in Section 3.2 was used in this case to extract the modal damping values of the beam mounted at the center for the first four natural frequencies by means of the half-power bandwidth method. The value of the damping 36 ratio is around 0.2% for all of the vibration modes. The values are shown in Table 3.3 and Figure 3.7. Frequency (Hz) Damping Ratio 181 0.001409 1078.75 0.002665 2937.75 0.001612 5560 0.002316 Table 3.3. Modal damping at four first natural frequencies Figure 3.7. Modal damping at four first natural frequencies Results from this set of experiments show no clear functional relation of the values of damping with the modes of vibration. They all lie in the same range, which is also the range found for samples with 0? and 90? in the testing of beams taken from the 37 plate where the laminae are aligned. It can be concluded that the damping ratio obtained with this setup is basically a constant value around 0.2%. 3.4 Cantilever Beams with Swept Sine Excitation of Base A series of experiments using beams in a cantilever configuration were conducted to compare the results with those from the measurements using bonded studs. It was hypothesized that the epoxy connection with the stud could be providing significant dissipation, so further investigation was required. Another concern was that for that experiment the frequency response measurement was being made between the head of the shaker and a point close to the tip of the beam, so any dissipation occurring: (a) in the connection of the stud and the shaker, (b) the stud itself, or (c) its connection to the beam, would be included in the measurement. Figure 3.8. Dog-bone shaped end of the sample For these experiments, the sample plate was the same as the one used previously for the measurement with a center attachment, with 22 layers in a [0,90,0,90,0,90,0,90,0,90,0] S configuration, a thickness of 3.16 mm, 240 mm of length 38 and 10 mm of width. An aluminum clamping base was used, in which special care was taken in matching the edges for an even boundary of the cantilever attachment of the beam. The section of the sample that connected to the base was shaped like the end of a dog-bone, in which the width of the beam was kept larger at the clamping end, as shown in Figure 3.8. This tended to place the point of maximum bending stress away from the connection to the base, in such a way that the damping in the connection would not substantially influence the measurement of the material damping. Figure 3.9. Measurement setup The base was mounted on an electromagnetic shaker, which provided a narrow band sine sweep transversal excitation of fixed acceleration passing through one of the three first natural frequencies of vibration. Miniature accelerometers were placed on the top of the base and at a point close to the tip of the beam, serving to measure the input for feedback control and the output response, respectively. Both accelerometers were LDS Dactron LASER Shaker Control System LDS V408 Electromagnetic Shaker Desktop computer ENDEVCO 22 - Miniature Piezoelectric accelerometers LDS PA 500L Power Amplifier 39 connected to a two channel charge conditioning amplifier, and the output signals were routed to a controller system to provide a measure of the vibration of the shaker for the feedback control of the excitation, and to obtain the transfer function between the two measurement points. The test setup is shown in Figure 3.9. The result was a series of curves of the form of that shown in Figure 3.10. From these curves the damping ratio could be calculated using the half power bandwidth Method, explained in Section 2.1.2. Figure 3.10. Experimental transfer function between base and tip of beam The results of the measurements are shown in Figure 3.11 as a function of the input acceleration. However, it proved of interest to examine the results as a function of the output displacement. This magnitude was not monitored directly, but it was obtained in the manner shown in Section 3.4.1. 320 325 330 335 340 345 350 355 360 0 10 20 30 40 50 60 70 80 frequency (Hz) m a gn i t u de ( g / g ) 40 0 0.2 0.4 0.6 50 50.2 50.4 50.6 50.8 51 input acceleration (g) nat u r al f r e que nc y (H z ) 0 0.2 0.4 0.6 1 1.5 2 2.5 x 10 -3 input acceleration (g) dam pi n g rat i o ? 0 0.2 0.4 0.6 339 339.2 339.4 339.6 339.8 340 input acceleration (g) n a t u ral f r equ enc y (H z ) 0 0.2 0.4 0.6 1.1 1.2 1.3 1.4 1.5 x 10 -3 input acceleration (g) da m p i n g r a t i o ? 0 1 2 3 937.8 938 938.2 938.4 938.6 938.8 939 input acceleration (g) na t u ral f r equ enc y (H z ) 0 1 2 3 3 4 5 x 10 -4 input acceleration (g) da m p i ng ra t i o ? Figure 3.11. Natural frequency and damping ratio vs. input amplitude acceleration for the three first natural frequencies 41 3.4.1 Relation between input acceleration and output displacement Taking the spans of the accelerations measured at the excitation (input) and response (output) points, the output accelerations can be represented as a linear function of the input accelerations as a out = A a in + B (3.18) For the first mode of vibration of the carbon epoxy sample, the input and output accelerations went from 0.1 g to 0.57 g and 11.35 g to 40 g, respectively, where g is the gravitational acceleration, 9.81 m/s 2 . Replacing these values adequately in Eq. (3.18), it is obtained that, for this mode of vibration, A = 60.957 and B = 5.25. Considering this result and the relation between acceleration and displacement for a pure sinusoidal motion, a = ? 2 d, the displacements of the tip of the beam at the first natural frequency (output displacement out I d ) can be obtained from the input accelerations using the expression 22 ( 60.9574 5.25) 9.81 4 50.5 in I out I a d ? ? + = , (3.19) where f 1 = 50.5 Hz and a in is in g?s (1 g = 9.81 m/s 2 ). In the same way, the displacements of the tip of the beam at the second and third vibration modes can be obtained. These are given by 22 ( 63.54 1.575) 9.81 4 339.5 in II out II a d ? ? + = , (3.20) and 22 ( 77.72 3.74) 9.81 4 938.5 in III out III a d ? ? + = . (3.21) Using Eqs. (3.19), (3.20), and (3.21) another representation of the natural frequency and damping plots presented in Figure 3.11 as a function of the displacement of the tip of the 42 beam can be obtained, as shown in Figure 3.12. Another way of displaying the results is to place the plots of modal damping vs. beam end displacement alongside with a fixed scale for the y axis of damping ratio, in order to compare their magnitudes, as shown in Figure 3.13. 1 2 3 4 x 10 -3 50 50.2 50.4 50.6 50.8 51 output peak displacement (m) 1 st nat ur al f r equ enc y (H z ) 1 2 3 4 x 10 -3 1 1.5 2 2.5 x 10 -3 output peak displacement (m) dam pi ng r a t i o ? 2 4 6 8 x 10 -5 339 339.2 339.4 339.6 339.8 340 output peak displacement (m) 2 nd nat ural f r equenc y (H z ) 2 4 6 8 x 10 -5 1.1 1.2 1.3 1.4 1.5 x 10 -3 output peak displacement (m) dam pi ng rat i o ? 0 2 4 6 8 x 10 -5 938 938.5 939 output peak displacement (m) 3 rd nat ural f r equenc y ( H z ) 0 2 4 6 8 x 10 -5 3 4 5 x 10 -4 output peak displacement (m) dam pi ng rat i o ? Figure 3.12. Natural frequency and damping ratio vs. displacement of the tip of the beam for the three first natural frequencies 43 1 2 3 4 x 10 -3 0 0.5 1 1.5 2 2.5 x 10 -3 da m p i n g r a t i o ? FIRST NATURAL FREQUENCY 50.5 Hz 2 4 6 8 x 10 -5 output peak displacement (m) SECOND NA TURA L FREQUENCY 339.5 Hz 0 2 4 6 8 x 10 -5 THIRD NA TURA L FREQUENCY 938.5 Hz Figure 3.13. Modal damping vs. displacement of the beam end 3.4.2 Discussion of results The set of measurements gave as a result that for different amplitudes of vibration, the natural frequencies of vibration (strictly speaking, the damped natural frequency) had different values. Also, the damping ratio showed a clear dependence on the amplitude of vibration. The slight decrease of the natural frequency with the amplitude of vibration corresponds to a clear geometrical nonlinearity. The effective stiffness of the beam decreases as the amplitude of vibration increases, so it behaves like a spring with a softening effect. The increase in damping ratio with increasing amplitude of vibration can be attributed to the intensification of the friction between the layers of the composite beam. Shear friction appears to be one of the predominant mechanisms of energy dissipation in composite materials, as concluded before from the study of the dependence of damping on fiber alignment. 44 3.4.3 Finite Element Model of Cantilever Beam Configuration In order to verify the results obtained from the experiments in Section 3.4.1, in terms of the spacing between natural frequencies, a finite element model of the beam mounted in a cantilever configuration with base displacement was developed. The value for the Young?s Modulus was adjusted in such a way that the first natural frequency of the model closely matched the first natural frequency obtained experimentally. The value for the Young?s modulus that resulted in a satisfactory agreement between the model and the experiment is E = 50.4 GPa, which is around 35% less than what is expected for an ideal plate, and is considered a reasonable deviation. The ratios between the natural frequencies of the three first modes of vibration closely match the analytical and experimental results. Another valuable result of this simulation is the series of plots shown in Figure 3.14, where the mode shapes of vibration at the first three natural frequencies are shown. Figure 3.15 shows the magnitude of the frequency response obtained from the simulation. 45 Figure 3.14. Mode shapes of vibration at the three first natural frequencies, obtained from finite element model 46 Figure 3.15. Magnitude of the frequency response obtained from finite element model 3.5 Axially Loaded Beams In order to apply the mechanical characteristics of the material in question to the modeling of flywheel systems, it was necessary that responses be observed for a variety of vibration amplitudes and natural frequencies. A wider range of frequencies than that considered in Section 3.2 had to be considered to approach the range of natural frequencies associated with a high speed flywheel system (on the order of 1 kHz). The shortest samples available could practically not have a natural frequency greater than 100 Hz. Shorter beams that would achieve higher frequencies yielded unreliable damping measurements due to end clamping effects that are difficult to control. In order to extend the measurement range, the samples were subjected to a tensile load so as to increase the effective natural frequency and, at the same time, include the effect of preload and the high levels of stress present in flywheel components. Attempts were made to use a tensile Frequency (Hz) Magnitude 47 testing machine in this regard, but the clamps used to fasten the samples allowed some lateral displacement that complicated the measurement of vibration while the samples were stretched. A specially designed test rig was developed and constructed to allow stretching of the samples with a tight attachment of the clamps. A photograph of the test rig is shown in Figure 3.16. The left side is fixed to the base by two large bolts and the right side can slide smoothly within the limits of the clearance between the fastening bolt and the associated hole. The desired tension is set by means of the fine pitched stretch control bolt on the far right, which pulls the sliding clamp towards a fixed block. Once the desired natural frequency for measurement is obtained, the vertical bolt is fastened, fixing the right end of the sample in that position. The test samples are made of carbon-epoxy in a [0?,90?,0?,90?,0?,90?,0?] configuration. They have a thickness of 1.05 mm, a width of 10 mm, and an effective length of 110 mm, measured between the innermost sides of the clamp fillets. These fillets machined at each end (a dog-bone configuration), were added to minimize the effect of the friction between the sample ends and the clamps in the overall vibration decay. Figure 3.16. Test rig Stretch control bolt Fastening bolt Fixed attachment 48 A strain gage was placed on the surface of the sample, aligned with its longitudinal axis, to determine the strain in the sample. This allowed for the calculation of the applied load (given knowledge of the Young?s Modulus of the specimen) and to relate the load applied with the first natural frequency of bending vibration. A dummy strain gage was bonded to a slab of the same material as the sample, to complete a half bridge configuration, which accounts for any thermal stresses occurring in the strain gage mounted on the sample. The measurement procedure consisted of setting the tension of the sample, applying an initial displacement and measuring the free decay of the amplitude of vibration of the first natural frequency, using a laser vibrometer focused at the center of the sample. The measured signals were recorded by a computer equipped with a data acquisition system, where further signal filtering was performed to isolate the vibration at the first natural frequency from small effects coming from resonance frequencies of the rig and other natural frequencies of vibration of the beam. The test setup is shown in Figure 3.17. For the analysis of the signal, the method of free damped vibrations was applied to blocks of data in the time domain. The length of each block was chosen to be of 30 cycles, i.e., including 31 peaks, after studying the correlation of results at different block sizes. This matched the criterion used in a similar study, in which 20 was determined as the minimum number of cycles to be considered for each block [43] . Vibration signals were acquired using an NI-4552 computer based data acquisition system. This system has an excellent amplitude flatness and very low total harmonic distortion. The sampling frequency used to register the vibration decays was 132,300 Hz, in order to obtain 49 properly shaped discretized representations of the sinusoidal decays in which the amplitudes (peak points) were as close to the actual values as possible. Larger sampling frequencies could not be handled by the buffer of the computer system used. The sensitivity of the laser vibrometer was set to 80 ?m/V, at which the full scale input limit is 1.3 mm and the resolution is 0.32 ?m. Figure 3.17. Experimental setup 3.5.1 Observed behavior Although the dependence of damping ratio on vibration displacement is widely recognized, there has been little work that has employed the method of free vibrations to assess such dependence [43] . Using the clamped-clamped configuration it has been observed that for vibration at significantly different amplitudes the damping ratio of the material greatly changes, as was observed in the measurement with a cantilever configuration as well, in Section 3.4 . This means that the vibration decay of a simulated simple degree of freedom mass-spring-viscous damper system does not adequately model Clamped Clamped LabVIEW VI Displ. ? Freq. 0 100 200 300 400 500 600-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 60 70 80 90 1000 50 100 150 200 250 300 LabVIEW VI Strain Laser Vibrometer Strain Gage Strain Gage Conditioner NI PCI Data Acquisition Vibrometer Controller Load Load Dummy Gage 50 the actual shape of a vibration decay obtained experimentally. In Figure 3.18 the black curve is formed by connecting the local peaks of the displacement of the center of a clamped-clamped beam with negligible loading in free vibration decay. The blue line was generated using a value of damping ratio ? = 0.00167, and it can be seen that it fits the decay rate observed at the beginning of the decay. However, as the displacement amplitude decreases, this line cannot follow the decay rate of the curve obtained experimentally. In the same way the red line, generated using a damping ratio value of ? = 0.000966, provides a good fit for the low amplitude region, but cannot follow the graph at high amplitudes. Figure 3.19, which is just another representation of the data in Figure 3.18, where the vertical axis is represented by a logarithmic scale, further illustrates the dramatic differences between the damping ratios at low amplitudes. As the displacement amplitude of the vibration increases, the value of the damping ratio increases as well. In dynamic systems, where stability can be highly dependent on internal damping, such increase may shift the effective stability thresholds considerably for some designs. 51 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 1 2 x 10 -4 time (s) d i sp l a ce m e n t ( m ) experimental ? =0.00167 ? =0.000966 Figure 3.18. Two fixed-exponential fittings of the top envelope of free vibration decay of the first natural frequency (593 Hz). Damping ratios differ by 72.8% 0 0.2 0.4 0.6 0.8 1 1.2 1.4 10 -7 10 -6 10 -5 10 -4 10 -3 time (s) di sp l a ce m e n t ( m ) experimental ? =0.00167 ? =0.000966 Figure 3.19. The same data as Figure 3.18 with y axis shown on a logarithmic scale 52 3.5.2 The Method of Free Damped Vibrations by Time Blocks The rate of reduction of free vibrations is typically quantified using the logarithmic decrement of vibration, ?, or the dissipation factor, ?, the corresponding relative energy dissipation [9] . The logarithmic decrement can be related to the damping ratio, ?, and also to the energy dissipation or damping capacity, ?. It is determined over several (n) cycles of the decay of vibration of a single degree of freedom system from the displacement amplitudes, using 1 ln 2 2 i in A nA ? ??? + ==?, (3.22) where A i and A i+n are the amplitudes corresponding to the i th and the (i+n) th cycles of the vibrations, respectively. The damping ratio describes the decay in the time response of a linear damped single-degree-of-freedom system subjected to an initial displacement, A, as shown in Eq. (3.23) -2 () cos( ) n d t xt Ae t ?? ? ?=+, (3.23) where ? d is the frequency of damped free vibration, ? n is the natural frequency and ? is the phase. The value of the damping ratio is the averaged characteristic of the energy dissipation in ?n? cycles of the vibration. For an amplitude-independent damping, the value of the damping ratio, ?, is unique and the classical spring mass damper system shown in Eq. (3.23) can model the vibration decay. However if the damping is amplitude-dependent, the value changes, and the damping can be related to the average amplitude in the range considered, (A i + A i+n )/2 [9] . For such cases, this amplitude dependency must be incorporated into the damping function if the dynamic behavior of the dynamic system being considered is to be 53 accurately modeled. One approach is to modify the damping function 2 n x??  to directly account for amplitude dependence. This can be done by adding a linear dependence on instantaneous displacement, in the form 0 2( ) n ax x? ?+  . Thus, the original equation of motion 2 2 0 nn xxx?? ?+ +=  , (3.24) becomes the modified equation 2 0 2( ) 0 nn xaxxx???++ +=  . (3.25) Please note that an assumption is made that the linear dependence on vibration displacement amplitude (described above and observed experimentally) will be preserved if a function of instantaneous displacement is used instead. This assumption allows generality and ease of implementation of our model. Figure 3.20 shows an example of the experimental decay of the vibration and the lines formed by the peaks of the decays of the two single degree of freedom models of Eqs. (3.24) and (3.25) (using a parametric ?best? fit to the experimental data). The model of Eq. (3.24) (using a constant damping ratio ? of 0.0011) provides a good fit to the decay of vibrations but cannot follow it properly, particularly at higher amplitudes of vibration. However, the model from Eq. (3.25) that takes into account the dependence of the effective damping ratio on vibration displacement can be seen to follow the envelope much more precisely. 54 Figure 3.20. Experimental free decay of vibration and the envelopes of the fittings with constant damping ratio and linearly changing damping ratio 3.5.3 Results Similar reference amplitudes were used to measure local damping ratios for all the time traces that were recorded. Four vibration decays were registered, at 593, 677, 735 and 744 Hz. The vibration amplitudes selected were in a range between 40 and 75 ?m, which was present in all of the measured signals. The results obtained for the damping ratio show clear trends with regard to dependence on frequency and vibration. Linear functions to describe the change in damping ratio, both with respect to frequency and with respect to displacement, are a logical first candidate. As seen in Figure 3.21, the linear functional form describes the dependence on sinusoidal vibration amplitude in a proper way, but the same cannot be said for the dependence on natural frequency in Figure 3.22. In these figures the dots represent the damping ratio values at some chosen 0 0.2 0.4 0.6 0.8 1 1.2 1.4 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x 10 -4 time (s) d i sp l a ce m e n t ( m ) experimental linear: ? =0.00094+8.4 x constant: ? =0.0011 ? = 0.0022 ? = 0.00096 ?(x)=.00094+8.4x ? = 0.0011 55 amplitudes of sinusoidal free vibration (in mm) and natural frequencies (in Hz) and the lines are the linear fits of the data. 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 x 10 -5 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 x 10 -3 Vibration amplitude (m) Dam p i n g ra t i o ? (x , f ) 593Hz 677Hz 735Hz 774Hz Figure 3.21. Damping ratio vs. vibration amplitude for different 1 st natural frequencies 550 600 650 700 750 800 850 1 1.1 1.2 1.3 1.4 1.5 1.6 x 10 -3 Frequency (Hz) D a m p i n g r a t i o ? 0.04 mm 0.044 mm 0.048 mm 0.052 mm 0.056 mm 0.06 mm 0.064 mm Figure 3.22. Damping ratio vs. frequency for various amplitudes of sinusoidal free vibration 56 The results shown in Figure 3.21 reflect what was observed before, that an increase in the amplitude of the sinusoidal vibration causes the damping to increase. Figure 3.22, indicates that an increase in the axial load applied, with the consequent increase in the natural frequency of vibration of the specimen, produces a decrease in the damping. It was not possible to separate the effect that changing the natural frequency had on damping from the effects on damping of the conditions causing the change of frequency, such as the load applied or contact with added masses, and thus the overall effect on damping of vibration frequency alone could not be asserted. 3.6 Natural Frequencies and Damping of a Sample Rotor The experimental study described above was extended to consider the characteristics of a sample flywheel hub design developed at Auburn University. Measurements of damping at natural frequencies were performed over a 236 mm diameter carbon fiber hub-rim interface, built using an epoxy matrix, with the fibers woven in the shape of two side domes with a center ring and mounted on a steel shaft. Measurements were also conducted over a complete rotor, including the rim mounted on the interface. In order to mount the rim, an axial tensile load was applied at each end of the hub, which serves to stretch the hub axially and (correspondingly) compress the hub diameter. Then, the rim could be slid into position and the hub released. An epoxy adhesive was placed at the connection surface between the hub-rim interface and the rim and allowed to set before any testing of the complete system was performed. 57 3.6.1 Measurement Setup An aluminum structure was constructed to mount the specimen rigidly. It consisted of a base holding two massive towers to which the shaft was clamped. The base was fixed to a plate connected to a high power electromagnetic shaker. The shaker was driven by random noise in the bands of interest (explained in Section 3.6.2) using compatible software. The input vibration and the response of the specimen to this stimulus were measured using two laser vibrometers. These signals were routed to a computer through a PCI data acquisition card and recorded. A photograph of the experimental is shown in Figure 3.23. Special care was taken in measuring the response parallel to the direction of the excitation axis. Figure 3.23. Measurement setup 58 3.6.2 Data Acquisition and Analysis The frequency response of the hub-rim interface was measured at the hub center between 0 and 5000 Hz. Using this broadband frequency response the main natural frequencies were first identified. Then each mode was excited separately using a narrow band random excitation, and the frequency response functions averaged until discernible smooth response curves were obtained. The frequency response function at the natural frequencies was curve fitted using polynomials and the damping ratio, ?, for each mode was then extracted from the resulting curve using the half-power bandwidth method. Natural Frequency (Hz) Damping Ratio 1589 0.0113 Shaft- Interface 3816 0.0125 Shaft-Interface- Rim 100 0.0279 Table 3.4. Natural frequencies and modal damping values Table 3.4 shows the values of damping obtained for each case (hub-rim interface with and without the rim mounted) and natural frequency. The hub-shaft system showed two relatively high natural frequencies (at 1589 Hz and 3816 Hz, respectively). The damping ratios were somewhat greater than 1%, dramatically higher than those observed for the coupon samples of similar material. Most likely, frictional interaction between the hub and the shaft accounts for this result. As expected, the fundamental radial natural frequency of the complete rotor (with the rim attached) was substantially lower, at 59 100 Hz, than for the hub-shaft alone. It also proved to have a damping level even higher than that of the shaft-hub system, at about 2.8%. So, it appears that a main source of internal damping for such systems is the internal friction between the various components of the rotor rather than the material damping associated with each individual component. This could become the dominant effect in the instability of a rotor, and thus special care must be taken in the mounting of the components. It is important to remark that after the rim was mounted, the hub-rim interface was in a state of compression, by the action of the rim. 60 CHAPTER 4 MODELING AND ANALYSIS OF FLYWHEEL SYSTEMS Rotating machinery has opened a range of possibilities in applications like power generation and energy storage, and is one of the most widely used elements in advanced mechanical systems. However, characterizing the vibration and instability phenomena that are associated with the operation of such devices can be a challenging task. It becomes necessary, when performing such analyses, to make some assumptions, and thus it is critical to recognize the restrictions of a model and the suitability when trying to adapt it to other similar analyses. Typical rotordynamic studies consider a flexible shaft and rigid disks and/or blades attached to it, an analysis that has been very useful for the characterization of systems such as steam turbines. Some studies have considered effects of disk flexibility as well, but it is common practice to neglect them. However, the dynamics of a flywheel system for energy storage introduces further complexity to the problem and demands another approach for its analysis. The fact that the energy is more effectively accumulated farther from the center of rotation makes it desirable to concentrate as much of the mass of the system in that location and to reduce the mass of other components that lay closer to the shaft. The search for new materials and construction methods for flywheel energy storage systems is a continuous process of testing and development. An auspicious 61 opportunity has been found in fiber reinforced composite materials. Their strength to weight ratio and the ability to tailor designs by aligning the fibers in the directions where the maximum stresses are expected have originated several studies of the feasibility and problems that could arise from its use. The main problems observed have to do with the fact that light and thin structures built with fiber reinforced polymers can withstand the stresses, but introduce flexibility, which can be detrimental to the stability. Consequently, components built with composite materials may possess such a degree of flexibility that modeling them as rigid would yield erroneous conclusions from the analysis, and thus the flexibility must be accounted for in the modeling stage. Moreover, the high damping levels of polymeric materials and friction arising from the interaction of different components of a composite are usually desirable, but in rotordynamics they can have dramatically harmful effects. The interface element between the hub and the outer rim is an attractive component to optimize, with the objective of reducing mass in mind and being aware of the long studied problems arising from shaft flexibility. A feasible design consists of winded carbon fibers in a polymeric matrix. A prototype of a rotor including this characteristic was previously presented in Figure 3.23. An analysis of this kind of rotor system is provided below. 4.1 Model of a Flywheel System In order to predict the critical speeds of a flywheel system with relatively flexible components that introduce damping forces acting between moving parts, a model is presented which takes into account translational as well as rotational degrees of freedom. 62 The model assumes a rigid shaft to which a hub disk is attached. This hub disk is connected to a rigid rim by means of a massless and flexible hub-rim interface, with associated stiffness and damping properties. The masses of the system are concentrated on the hub and on the rim, since their masses are significantly greater than those of the other components, as depicted in Figure 4.1. However, if desired, the shaft can be incorporated in this model as well, in the mass and inertia terms of the hub. A noteworthy observation is that this kind of system tends to have a short shaft span in between bearings, and thus the flexibility of the shaft becomes less of an issue as compared to the case of other kinds of turbomachinery, such as multi-stage centrifugal compressors or turbo-pumps. The illustration in Figure 4.1 solely shows a non-proportional representation of the parts involved for a clear understanding of their degrees of freedom and interaction. Figure 4.1. Eight degree of freedom model of a flywheel system x ? y ? z ? 63 Motion of the system is described by using the 8 generalized coordinates ? R , ? R , ? H , ? H , x R , y R , x H , and y H , representing the angular and translational motions of the rim and the hub. Damping and stiffness parameters for the bearings and shaft-rim interface are considered to be symmetric. The spin speed ? is considered constant. The parameters involved in the analysis are: c ? : bending damping coefficient of shaft-rim interface k ? : bending stiffness of shaft-rim interface c x : extensional damping coefficient of shaft-rim interface k x : extensional stiffness of shaft-rim interface c B? : equivalent bending damping coefficient of bearings k B? : equivalent bending stiffness of bearings c Bx : equivalent extensional damping coefficient of bearings k Bx : equivalent extensional stiffness of bearings Rotations of the shaft-hub and rim, are considered to occur in the order: ? about ?y , ? about x', and ? about z'', where ? ?=  , as shown in Figure 4.2. The total angular velocities of the hub and the rim, have the form '' 'zx y? ???= + +   . (4.1) Then, considering the relation between the coordinate systems of Figure 4.2, the angular velocity of each component can be expressed in terms of the body fixed coordinate system x yz: 64 Figure 4.2. Relations between coordinate systems (cos( ) - sin( ) ) (cos( ) '' - sin( ) ''), (cos( ) - sin( ) ) cos( )(sin( ) cos( ) ) - sin( ) ), ( cos( ) cos( )sin( )) ( cos( )c zxy yz zxy xyz x ??? ? ?? ? ??? ? ??? ??? ????? ?? ?= + + =+ + + =+ +      os( ) - sin( )) ( - sin( )) .y z?? ? ?? ?+   (4.2) Also, using a space fixed coordinate system, the linear velocity of each component can be expressed as: 22 ?? , and . GG GG vxxyy vv x y =+ ?= + G  GG  (4.3) x', x'' y' ? z'' z' y'' ?z ?y , y' ? ?x x' z' y'' ? x'' z'', z y x 65 The expression for the kinetic energy, T, of each component is: {}[]{} 22 22 1111 ( ) ( ) 2222 T GG G GG G T mxy H mxy I=++??=++?? G G   , (4.4) where {} cos( ) cos( )sin( ) cos( ) cos( ) - sin( ) - sin( ) ? ???? ? ???? ?? ? ??+ ?? ?= ?? ??      , and [] 00 00 00 G It IIt Ip ? ? ? ? = ? ? ? ? ? ? . (4.5) Then Eq. (4.5) can be approximated as 22 22 2 111 ( ) ( ) ( -2 ) 22 GG T m x y It Ip? ?????=++++   , (4.6) considering that the angular displacements are small. The total kinetic energy of the system is 22 2 2 2 22 2 2 2 11 1 ( ) ( ) ( - 2 ) 22 2 11 1 ( ) ( - 2 ) ( ) . 22 2 RR R TH H R RR RR R T HH TH H Tmxy mxy Ip It Ip It ???? ?? ? ??? ?? =++++ + +++ ++    (4.7) The potential energy, V, consists only of the elastic energy on the bearings and on the hub-rim interface, and can be written as 2222 1111 (- ) (- ) 222 ( - ) ( - ) . 2 xR H xR H BxH BxH RH R H BH BH Vkxx kyy k ky kk ?? ?? ?? ?? ? ? =++++ ++ (4.8) The dissipation forces acting on the hub-rim interface in the translational degrees of freedom, x F G , expressed in the rotating (rot), body fixed reference frame x yz are ( ) - ( - ) ( - ) rot rot rot rot xxRH RH F cxxxyyy=+ G  . (4.9) 66 The velocity components must be transformed into the space fixed rotating frame. The components of position in the rotating frame, and rot rot x y , expressed in terms of the space fixed or stationary reference system, and stst x y , are cos( ) sin( ) , and cos( ) - sin( ) . rot st st rot st st xx ty t yy tx t ?? ?? =+ = (4.10) Eqs. (4.10) are differentiated with respect to time, arriving at: cos( ) - sin( ) sin( ) cos( ) , and cos( ) - sin( ) - sin( ) - cos( ) . rot st st st st rot st st st st x x tx ty ty t yy ty tx tx t ??? ? ?? ??? ??? =++ =     (4.11) The unit vectors and x y of Figure 4.2 must also be expressed in terms of ?? and x y , as () '' cos( ) ''sin( ) ''cos( ) sin( ) ( 'cos( ) 'sin( )) ??? cos( )( cos( ) - sin( )) sin( ) cos( ) sin( )( cos( ) sin( )) xx ty t xt ty z tx z t y z x ? ? ?? ? ? ? ??? ??? ? =+ =+ + ++ (4.12) () '' cos( ) - ''sin( ) cos( )('cos() 'sin())- 'sin( ) ????? cos( ) cos( ) sin( )( cos sin( )) - sin( )( cos - sin )() () yy tx t ty z x t ty z x tx z ? ? ??? ? ??????? ? = =+ + (4.13) which, for small angles ? and ? become () () () () ?? cos sin , and cos - sin . ttxxy yyx ??= + = (4.14) So the velocity components for each part, hub and rim, are ( ) ( ) ?? - rot rot st st st st x xyy x y x yx y??+=+ +  , (4.15) and the translational dissipation forces in the hub-rim interface have the form: 67 ()( )?? - ( - ) ( - ) ( - ) - ( - ) x x RH RH RH RH Fcxx yyx yy xxy??=++?? ?? G   . (4.16) Analogously, the components of angular position in the rotating frame, and rot rot ? ? , in terms of the space fixed reference system, and stst ? ? , are () () () () cos sin , and cos - sin , rot st st rot st st tt tt ?? ?? ?? ? ?? ? =+ = (4.17) which after differentiation provide: cos( ) - sin( ) sin ( ) cos( ) , and cos( ) - sin ( ) - sin ( ) - cos( ) . rot st st st st rot st st st st tt t ttt t ?? ????????? ?? ????????? =++ =     (4.18) So the angular velocity components for each part, hub and rim, are ( ) ( ) ?? - rot rot st st st st x yxy?? ??????+=+ +  , (4.19) and the expression for the rotational dissipation forces, F ? G , in the hub-rim interface is ()() ??- ( - ) ( - ) + ( - ) - ( - ) . RH R H R H RH Fc x y ?? ?? ??? ?? ??? ?? =+ ?? G   (4.20) The translational dissipative forces due to the bearing flexibility, Bx F G , are ( )??- Bx Bx H H Fcxxyy=+ G , (4.21) and the rotational dissipative forces, B F ? G , are: ( ) ??- BBH H Fc x y ?? ??=+ G   . (4.22) The virtual work done by these forces is the product of the forces by the virtual displacements in the corresponding directions: ()( ) ( ) ()() ( ) ???? - - - - . xRH RH BxH H RH R H B H H WF x xx y yy F xx yy Fx ?? ????? ?? ?? ?? ?? ?? ?? ?? =? + + ? +?? ?? +? + + ? + GG (4.23) 68 By means of Lagrange?s equations, which state: - i ii dL L Q dt q q ???? = ?? ?? ??  , (4.24) where L is the Lagrangian, L = T ? V, and Q i represents the generalized force for the coordinate q i , a system of eight equations describing the dynamics of the model is obtained, which can be written as RR R R ? R ? H ? R ? R ? H ? H RR ? RRR? H ? R ? R ? H ? H HH ? RHh? B? H ? R ? R ? H ? It ? + Ip ? ? + c ? - c ? + ? c ? + k ? - ? c ? - k ? = 0 It ? + c ? - Ip ? ? - c ? + k ? - ? c ? - k ? + ? c ? = 0 It ? - c ? + Ip ? ? + (c + c ) ? - ? c ? - k ? + ? c ? + (k +             B? H HH ? R ? B? HHH? R ? R ? B? H ? H RR xR xH xR xR xH xH RR xR xH xR xR xH xH k) ? = 0 It ? - c ? + (c + c ) ? - Ip ? ? - k ? + ? c ? + (k + k ) ? - ? c ? = 0 mx + cx - cx + kx + ? cy - kx - ? cy = 0 my + cy - cy - ? cx + ky + ? cx - ky =           HH xR x Bx H xR xR x Bx H xH HH xR x Bx H xR xR xH x Bx H 0 m x - c x + (c + c ) x - k x - ? c y + (k + k ) x + ? cy = 0 m y - c y + (c + c ) y + ?c x - k y - ? c x + (k + k ) y = 0.       (4.25) Terms can be grouped to express the equations in the matrix form [ ] [ ] [ ] 0My C Gy Ky++ + =    , (4.26) where R R H H R R H H y x y x y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ??  , R R H H R R H H y x y x y ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ??           , (4.27) 69 [] 0000000 0000000 00 0000 0 000 000 0 0000 00 0 00000 0 0 000000 0 0000000 R R H H R R H H It It It It M m m m m ?? ?? = ?? ?? , (4.28) [] -- 0000 -- 0 -- 0 0 0 0 -- 00 0 0 00 0 0 - - 00 0 0- - 00 0 0 -- 00 0 0 - - B B xx x x xx x x xxxBxx x xxxBx kck c ck c k kckk c ck ckk K kck c ck c k kckk c ck ckk ?? ? ? ?? ? ? ????? ?? ??? ?? ?? ?? ?? ?? ?? ?? ?? ?? ??+ + = + + ?? , and (4.29) [] --00000 0-000 0 -0 - 000 0 0 - 00 0 0 00 0 0 0- 0 00 0 00 0 - 00 0 0-0 0 00 0 00- 0 R R BH HB xx xx xxBx x xBx cIpc Ip c c cccIp cIpcc CG cc cc ccc ccc ?? ?? ??? ??? ? ? ? ? ?? ?? + + += + + ?? ?? . (4.30) Then, in order to express this system in the state space form [ ] x Ax=  , (4.31) (4.32) 81 81 () () () x x yt xt yt ? ? ? ? = ? ? ? ? ??     and 70 , (4.33) are defined. 4.1.1 Model parameters Suitable values for the parameters involved in the model of Eq. (4.31) were determined. Some of these parameters will be varied, with the goal of obtaining useful information for the optimization of designs. However, if not otherwise stated, the following values will be the ones used for analysis. The properties assigned to the rim are: ro R = Outer radius of rim = 0.16 m, ri R = Inner radius of rim = 0.08 m, ? comp = Volumetric density of Carbon-Epoxy = 1400 kg/m 3 , w R = Width of rim = 0.06 m, and m R = Mass of rim = 22 (- ) OR IR R C rrw? ? = 5.1 kg. From these physical properties the polar and transversal moments of inertia of the rim, Ip R and It R, can be calculated using Ip R = 22 1 ( ) 2 RR R mro ri+ = 0.081 kg?m 2 . It R = 222 1 (3( ) ) 12 RRRR mroriw++ = 0.042 kg?m 2 , and Similarly, the properties for the hub are [] [][ ] [][] [] [] -1 -1 88 88 16 16 - - 0 xx x MCG MK A I ?? + = ?? ?? 71 ro H = Outer radius of hub = 0.06 m, ri H = Inner radius of hub = 0 m, w H = Width of hub = 0.07 m, ? Al = Volumetric density of aluminum = 2700 kg/m 3 , m Hub = Mass of hub = 22 (- ) H HHAl ro ri w? ? = 2.14 kg, It Hub = Transversal moment of inertia of hub = 222 1 (3( ) ) 12 HHHH mroriw++ = 0.0028 kg?m 2 , Ip Hub = Polar moment of inertia of hub = 22 HH H 1 m (ro + ri ) 2 = 0.0038 kg?m 2 . Stiffness and damping parameters must be defined in the translational and rotational degrees of freedom for the hub-rim interface and for the bearings. The values chosen are k B? = 12300 N/rad, k Bx = 1.9e6 N/m, k ? = 37000 N/rad, k x = 6 x 10 6 N/m, ? B? = Bearing rotational damping ratio = 0.02, ? Bx = Bearing extensional damping ratio = 0.02, ? ? = Hub-rim interface rotational damping ratio = 0.0015, ? x = Hub-rim interface extensional damping ratio = 0.0015, c B? = 2 B BH H k It It ? ? ? = 0.23 N s /m, c Bx = 2 Bx Bx H H k m m ? = 82.2 N s /m, 72 c ? = 2 R R k It It ? ? ? = 0.12 N s /m, and c x = 2 x x R R k m m ? = 16.4 N s /m. Extracting the eigenvalues of A for the parameters given above and the running speed, ?, varying from 0 to 80,000 RPM, plots of the imaginary and real parts of the eigenvalues of the state space system are obtained. These plots are shown in Figure 4.3 and Figure 4.4. Examination of the results of the simulation using the model of the complete rotor, with consideration to the shape and structure of the hub rim interface, allow to argument that rotational modes of vibration have natural frequencies that are considerably larger than those of translational modes. The forward whirling modes that represent critical speeds, i.e., the ones that intersect with the line describing the frequencies equal to the running speeds in Figure 4.3 are all translational. The rotational modes do not intersect that line but in their backward directions (negative slopes). Moreover, since the hub-rim interface is built in a way that renders it relatively stiff in those directions and the ratio IpH/ItH is close to 2 (because the hub-interface-rim system is close to a thin disk), these natural frequencies increase rapidly with running speed and do not represent critical speeds, thus stability problems are not expected to occur in these modes. 73 Figure 4.3. Imaginary part of eigenvalues for a range of running speeds Figure 4.4. Real part of eigenvalues for a range of running speeds The advantages of having a rotor with a thin disk, in terms of stable regions, have been pointed out above. The following analysis shows how the dynamics of the system 2 4 6 8 10 12 14 x 10 4 -70 -60 -50 -40 -30 -20 -10 0 Running speed (RPM) R eal par t of ei ge nv al ues ( H z ) 6.5 7 7.5 8 8.5 9 x 10 4 -3 -2 -1 0 1 *** Rotational modes *** Translational modes *** Translational modes 0 1 2 3 4 5 6 x 10 4 0 200 400 600 800 1000 1200 Running speed (RPM) F r equ enc y ( H z ) Critical speeds *** Rotational modes *** Translational modes *** Translational modes *** Running speed = frequency 74 would change if the rim was enlarged in the axial direction. A rotor with the same characteristics as before has been modeled, with the difference that the mass of the rim has been distributed in a thinner but longer rim. The new rim is 40 cm long (compared to 6 cm of the original). The radii are roR = 12.6 cm and riR = 11.4 cm, maintaining in this way the center radius of 12 cm. In Figure 4.5 the imaginary parts of the eigenvalues of this rotor, for a range of running speeds, are shown. Figure 4.5. Imaginary parts of eigenvalues for a long symmetric rotor What can be seen from this figure is that the forward (positive slope) rotational mode with a higher frequency has increased its ?root? frequency (the frequency at 0 RPM) to over 1 kHz, due to the change in the moments of inertia. At the same time the rotational mode with lower frequency has maintained its root frequency, while the increase in the frequency of the forward mode (the positive slope) became smaller. This occurs because of the decrease in the ratio of the moments of inertia, Ip R /It R , from almost 0 1 2 3 4 5 6 x 10 4 0 200 400 600 800 1000 1200 Running speed (RPM) F r equ enc y ( H z ) Out-of-phase rotational mode In-phase rotational mode *** Rotational modes *** Translational modes *** Translational modes *** Running speed = frequency 75 2 (for a thin disk) to about 0.7 for this case. The main consequence of this occurrence is that this mode now intersects the line describing the frequencies equal to the running speeds and thus becomes a critical speed, with the potential of becoming an unstable mode for a certain running speed. An analysis of the eigenvectors, which describe the mode shapes of the vibrations at the frequencies corresponding to the eigenvalues in Figure 4.5, indicates that the rotational mode with a higher frequency corresponds to the out-of-phase mode, i.e. the mode in which the shaft and hub tilt to one side while the rim tilts to the other. The eigenvectors used for this analysis are shown in Table 4.1. It can be seen that for 1983.4 Hz, the rotational mode with higher frequency, the real parts of the eigenvectors have ??s and ??s (of the rim and hub) with opposite signs, while for 751.7 Hz, their signs are equal. Imaginary part of eigenvalue (Hz) -1983.4 1983.4 -751.7 751.7 ? R 0.0000 + 0.0005i 0.0000 - 0.0005i 0.1337 - 0.0000 0.1337 + 0.0000i ? R 0.0005 - 0.0000i 0.0005 + 0.0000i -0.0000 - 0.1337i -0.0000 + 0.1337i ? H -0.0006 - 0.0567i -0.0006 + 0.0567i 0.0673 - 0.0009i 0.0673 + 0.0009i Eigenvectors ? H -0.0567 + 0.0006i -0.0567 - 0.0006i -0.0009 - 0.0673i -0.0009 + 0.0673i Table 4.1. Imaginary parts of eigenvectors (x 10 -3 ), rotational displacement components 4.1.2 Translational Model of a Flywheel System The observation made in the previous section about the stability of the rotational modes and the fact that rotational and translational modes appear decoupled in the 76 previous analysis, allows us to further simplify our model to a purely translational model, in which rotational effects are neglected. So our original system now can be transformed into that shown in Figure 4.6, where the elastic bands with translational and rotational degrees of freedom have been replaced by springs and dampers in the x and y directions. This simplification greatly reduces computational time. Figure 4.6. Translational model for Flywheel with flexible hub-rim interface (Solid Edge drawing by Alex Matras [48] ) The reduced model can be represented by the system of equations - ( ) - - ( ) 0 - ( ) - - ( ) 0 - - - 0 HH xR x BxH xR xR x BxH xH HH xR x Bx H xR xR xH x Bx H RR xR xH xR xR xH xH RR mx cx c c x kx cy k k x cy my cy c c y cx ky cx k k y mx cx cx kx cy kx cy my ? ? ?? ++ ++ + = ++ + ++ = ++ =           - - - 0. xR xH xR xR xH xH cy cy cx ky cx ky?? (4.34) 77 Complex variable notation z = x + iy can be adopted to simplify the analysis, allowing us to express Eq. (4.34) as (-) (-) (-)0 (-) (-)- (-)0. HH BH BxH x H R x H R x R H RR xR H xR H xR H mz cz k z c z z k z z ic z z mz czz kzz iczz ? ? ++ + + + = =        (4.35) The terms are reorganized and the parameters THR mmm= + , (4.36) H T m a m = , (1 - ) R T m a m = , (4.37) 2 (1 - ) (1 - ) xx x x R TH kk ak mam m ? == = , (4.38) and 2 BxBx B TH kk mam ? == (4.39) are defined for expedite analysis. Then Eq. (4.35) can be expressed as 22 2 (1 - ) (-) (-) (-)0 (- ) (- ) - (- )0, (1 - ) (1 - ) xxxBB HHH HR HR RH TTT xx RRHxRH RH ac icc zzz zz zz zz am a a am am cic zzzz z am am ?? ? ?? +++ + + = =        (4.40) in a way that is more comfortable for design purposes. Eq. (4.40) can then be expressed in the matrix form [ ]{ } [ ]{ } [ ]{ } 0Mz Dz Kz+ +=   (4.41) in order to introduce the equations of motion in a computer program and perform eigenvalue analyses, where H R z z z ? ? = ? ? ?? , H R z z z ? ? = ? ? ??    , H R z z z ? ? = ? ? ??    , (4.42) 78 [] 10 01 M ?? = ?? ?? , [] ( )- 1 - (1 - ) (1 - ) Bx x xx T cc c aa D cc m aa + ? ? ? ? ? ?= ? ? ? ? ? ? , and (4.43) [] 22 2 22 (1- ) - - (1- ) - - (1 - ) (1 - ) xx Bx x TT xx xx TT cc ai a i mm K aa cc ii am am ?? ? ? ? ?? ?? ???? ++ ???? ???? = ?? ???? + ???? ???? ?? . (4.44) In state space form, the equations are 22 2 22 0100 (1 - ) - (1 - ) () -- 0001 --- (1 - ) (1-) (1 - ) (1 - ) xx Bx x H H TTBx x H H TT R R R xx xx TT TT cc ai a i z z mmcc c z z aam am z z z cc cc ii am m am am ?? ? ? ? ?? ?? ?? ?? ???? ++ ?? ???? + ???? ?? ?? = ?? ?? ?? ?? ??? ? + ??? ? ??? ? ??       R z ? ? ? ? ? ? ? ? ? ? ? ? ?? . The model parameters for the CFRC flywheel system are design variables that were chosen for convenience. The set of model parameters used for this part of the study, unless otherwise noted, are shown in Table 4.2. 79 Parameter Value m T 10.15 kg a 0.3 ? x 628 rad/s ? B 53.7 rad/s c x 165.8 kg/s c B 84 kg/s Table 4.2. Model Parameters Variation studies were conducted to assess the influence of damping and stiffness of the hub and bearings, and running speed on rotor dynamic stability. Figure 4.7 shows a simultaneous variation of the stiffness value for the hub-rim interface and of its damping ratio. Each point on the lines (or surface) represents the threshold running speed above which the rotor becomes unstable for a certain parameter configuration. This means that there is a distinctive operating speed below which the system is always stable for a given parametric configuration. As the level of internal damping in the hub-rim interface is increased, this threshold speed steadily decreases. Also, it should be noted that there are two distinctive vibratory modes for this model, which are those one would expect for any two mass system connected by springs. One consists of an in-phase mode in which the rim and shaft-hub move essentially in the same direction, but at generally different amplitudes. The other is an out-of-phase mode in which the rim and shaft-hub move in opposition to each other. 80 As the hub stiffness increases, there is a breakpoint in each of the curves plotted in Figure 4.7. These breakpoints are associated with a transition from the dominant mode (destabilized at a lower rotor speed) being the mode where hub and rim move out-of- phase to it being the one where they move in-phase. Figure 4.8 shows how a rotor of these characteristics and with a low stiffness will be destabilized by the out-of-phase mode, while a higher stiffness will yield the in-phase mode unstable at a lower running speed, as shown on Figure 4.9. Figure 4.7. Maximum stable running speed for hub damping ratio ? H = 0.002 to 0.02 and hub stiffness k H = 0 to 100000 kg/s 2 81 Figure 4.8. Out-of-phase mode destabilizes first for k H = 20000 N/m (? = 0.02) Figure 4.9. In-phase mode destabilizes first for k H = 50000 N/m (? = 0.02) 0 2000 4000 6000 8000 10000 12000 14000 -25 -20 -15 -10 -5 0 5 10 15 Running Speed (RPM) R eal par t o f ei ge nv al ues *** In-phase mode *** Out-of-phase mode ? Stability threshold 0 2000 4000 6000 8000 10000 12000 14000 -25 -20 -15 -10 -5 0 5 10 15 Running Speed (RPM) R eal part of ei g env al u e s *** In-phase mode *** Out-of-phase mode ? Stability threshold 82 Another simultaneous parameter variation was made for the parameters that characterize the bearings, keeping all the other parameters fixed. For this study the stiffness of the hub-rim interface was 2,760,688 N/m, equivalent to ? H = 628 rad/s, as specified above. The bearing damping ratio was varied from 0.015 to 0.1 (1.5 to 10 %) and the stiffness went from 0 to 100,000. The results of this study are shown in Figure 4.10, where it can be seen that an increase in bearing stiffness is beneficial in the sense of increasing the maximum predicted running speed, effect that is more evident for higher damping ratio values. At the same time, it can be observed that a higher bearing damping ratio is in general beneficial, while this benefit gets more significant as the bearing stiffness increases. Figure 4.10. Maximum stable running speed for bearing damping ratio ? B = 0.015 to 0.1 and bearing stiffness k B = 0 to 100000 kg/s 2 83 4.1.3 Experiment of Rotor with Flexible Hub-Rim Interface In order to gain greater insight on the qualitative behavior of this kind of system where the hub-rim interface flexibility is significant, a rig was built to operate and exhibit instability at relatively low running speeds, for safety issues. This was achieved by fixing a rubber hub-rim interface between two aluminum rims and two aluminum hubs, and mounting the structure on an aluminum shaft, so that the flexibility of the rubber interface was much more significant than that of the shaft in bending. The setup that was used is shown in Figure 4.11. Figure 4.11. Experimental Set-Up The mass of the rim was significant, but the deflection of the stretched rubber due to gravity was minor. The shaft was mounted on journal bearings and connected to a bolts shaft rim rim rubber bearing interfac bearing 84 servo motor. Running speeds were monitored by use of a proximity probe positioned close to the shaft. In order to determine the system parameters before run-up a setup consisting of a proximity probe and a signal analyzer was used. The shaft was made as rigid as physically possible, by placing the bearings close to the center disc. Impulsive force was applied and the decay of vibrations was recorded for the shaft on the bearings and for the rim on the flexible interface. Hub Shaft/Bearing Damping ratio (?) 0.0073 0.115 Spring constant, k (N/m) 26994 129009 Damping constant, c (N?s/m) 1.84 51.129 Natural Frequency (rad/s) 221.8 742.3 Table 4.3. Average values of experimental data The effective masses of the shaft and rim were measured to be m S = 0.238 kg and m R = 0.567.kg, respectively. The natural frequency of radial vibration was obtained assuming the damped frequency to be the same as the undamped natural frequency (which is a good assumption considering the low values of damping that were measured) and, from knowledge of the mass, the hub stiffness was calculated. By assuming the damping to be viscous, the damping ratio was determined using the ratio of amplitudes of Eq. (2.8). The same process was performed on the shaft alone, mounted on the bearings, 85 to determine its natural frequency and the damping provided by the contact with the bearings. Table 4.3 shows the calculated stiffness and damping values. Figure 4.12 shows frames from a video that was taken during the passage of the rotor (using a quasi static speed increase) through the critical speed. The first picture depicts the rotor while spinning at about 4200 rpm. Then at the critical speed (4980 rpm), the forward whirl amplitude increased abruptly, as shown in the two other pictures. As a result, the rubber was torn into pieces and the steel shaft bent to about 20 degrees. Figure 4.12. Rotor rig experiencing unstable behavior What was clear from the observation of the video can still be seen on the pictures: the whirl motion in the radial direction is significant compared to that on the transverse direction, i.e., rotation with respect to an axis perpendicular to the shaft. This reaffirms what was concluded in the eigenvalue analysis for a rotor with a thin rim. 86 CHAPTER 5 CONCLUSIONS A detailed study of damping in carbon-fiber epoxy composite structures has been conducted. The work has consisted of a series of experimental and simulation studies aimed at assessing the magnitude of damping and the influence of vibration amplitude and frequency on damping amount. This work has particularly considered the effects of damping in carbon fiber epoxy composite materials for application to flywheel energy storage systems. In the modeling of such systems, this material acts as an interface between the shaft-hub and the rim of such a system. First, a number of different configurations of fiber reinforced epoxy composites were experimentally evaluated. These were (1) sample beams mounted on a bonded stud attached to a shaker, on and off center, (2) beam samples supported in a cantilever configuration where the cantilever side was excited, (3) clamped-clamped beam samples with an axial load applied to the attachments, and (4) a prototype of a rotor with a steel shaft and hub, and fiber reinforced composite polymer hub-rim interface and rim. Results of all these experimental quantification of vibration damping studies have been documented and the results discussed in detail in this dissertation. Some particularly interesting results are: ? A vibration damping analysis of quasi-unidirectional composite material beam samples shows a clear dependence on the alignment of the fibers. The highest 87 values of vibration damping were registered for an alignment of the fibers of 45? with respect to the longitudinal axis of the beam. This effect has been observed in previous studies and it is attributed to the fact that at off-axis angles the shear stress is higher, having its maximum at an angle close to 45?. ? The measured damping ratios are in the range of 0.1 to 0.4% for all beam sample experiments, and for cross ply composites, which are the ones of more interest for this study since they better simulate a real component, most results show values around 0.2% or less. ? Experiments on beam samples mounted on bonded studs show slightly higher values of damping than those mounted in cantilever or clamped-clamped configurations. This suggests that the benefit of contacting the sample in only a small area with a stud comes with the disadvantage of exposing it to contact with the epoxy bonding, which appears to substantially increase (and undesirable) energy dissipation and the resulting damping ratios. ? The dog-bone configuration utilized in some of the experimental studies is an excellent way to reduce the dissipative effect of contact with the experimental sample and the testing apparatus. The values for vibration damping from that configuration are noticeably lower than those obtained using the bonded stud configuration and the values obtained from the cantilever and clamped-clamped configurations are in good agreement. ? The values of vibration damping obtained from the measurements of frequency response of the rotor prototype, with and without the rim mounted, are much 88 higher (about 10 times) than the values obtained from the measurements of the beam samples. This indicates that although material damping may play a significant role, friction between moving elements is the critical factor in the overall internal damping. ? Nonlinear effects have been considered in the measurements using the cantilever configuration and the clamped-clamped configuration. An increase of the vibration damping with the amplitude of the vibrations is observed in all of the measurement results. However, the relative significance of this effect may be dependent on the sample configuration and be substantially different for configurations other than those that were tested. ? The vibration damping does not depend on the frequency of vibrations in the experiments made with the quasi-unidirectional samples. On the other hand, in the clamped-clamped experiment, some variations in damping with natural frequency (achieved by loading the samples) are observed. However, it was not possible to isolate the influence of other parameters and establish a clear functional dependence on the frequency alone. In the modeling of flywheel systems in this work a design was considered in which the flexibility of the system is assigned entirely to the bearings and the hub-rim interface, including as well the associated damping effects. Some concluding remarks on this respect are given below: ? The configuration of the prototype tested makes it considerably stiffer in the axial direction than the radial one. However, regardless of this consideration, the translational modes were shown with this study to be more susceptible to 89 instability problems than the rotational even for similar stiffness values. This is further exaggerated by the fact that the large ratio between the moments of inertia, Ip/It, for a thin disk (or short rotor) produces rotational modes of vibration at frequencies which do not constitute critical speeds. ? It has been shown that for each practical configuration a safe range of operation below the threshold speed of instability can be determined. This threshold depends on the overall parametric configuration, but some clear observations can be made about the relative importance of certain parameters. o An increase in the amount of damping in the hub-rim interface causes the operation range to be reduced steadily. o For a hub-rim interface with very low stiffness, the mode that in general becomes unstable at a lower running speed is the in-phase one, in which the rim and shaft-hub move essentially as a whole. As the stiffness of the hub-rim interface is increased, there is a clear breakpoint after which the mode that becomes unstable at a lower running speed is the out-of-phase, which shows a motion characterized by the rim and hub moving in opposite directions. o A study involving the variation of the stiffness and damping of the bearings shows that an increase in the bearing stiffness is beneficial for the stable range of operation, especially if the value of damping is high. At the same time, a higher damping ratio is in general beneficial, while this benefit gets more significant as the stiffness increases. Some specific fundamental contributions of this work are: 90 ? A clear nonlinear functional behavior of carbon fiber reinforced composite polymers has been identified. ? The potential instability resulting from the flexibility of the hub-rim interface, for which there is scarce treatment in the literature on rotordynamic stability, has been examined in some detail. ? A method of analysis for flywheel rotordynamic stability, considering the particular characteristics of this kind of systems and possible simplifications, has been developed and presented. ? Key factors that determine stability of flywheel systems and their interactive roles have been identified and analyzed. There is certainly a great deal of fertile ground for further investigations on the topic of material damping in composite materials and structures. Some suggestions for future work in this area are: ? Further analysis of the dependence of damping in fiber reinforced composites on natural frequency and their response to forced harmonic excitation out of resonance. ? An investigation of the effects of voids and damage on damping and natural frequency of composite materials. ? A detailed investigation of the influence of non-symmetrical and nonlinear damping effects on the dynamic behavior and stability of rotating composite structures. 91 ? More detailed modeling and analysis of flywheel systems in which the flexibility of the hub-rim interface and the rotor shaft are both considered significant effects. 92 REFERENCES [1] Gabrys, C. W. and Bakis, C. E. Composite Materials: Testing and Design, 13th Vol., STP 1242, S. J. Hooper, Ed., American Society for Testing and Materials, Conshohocken, PA, Design and Testing of Composite Flywheel Rotors, 3-22, (1997). [2] Beranek, L. L. and Ver, I. L. Noise and Vibration Control Engineering, John Wiley and Sons, (1992). [3] Linacre, E. ?Damping Capacity, Part I: Introduction and Technique of Measurement?, Iron and Steel, May, 153-156, (1950). [4] Linacre, E. ?Damping Capacity, Part II-Effects of Conditions of Measurement?, Iron and Steel, June, 285-286, (1950). [5] Crandall, S. H. ?The Role of Vibration Damping in Vibration Theory?, Journal of Sound and Vibration, 11 (1), 3-18, (1970). [6] Lazan, B. J. Damping of Materials and Members in Structural Mechanics, Pergamon Press, (1968). [7] Nashif, A. D., Jones, D. I. G. and Henderson, J. P. Vibration Damping, John Wiley and Sons, (1985). [8] Steidel, R. F. An Introduction to Mechanical Vibrations, John Wiley & Sons, (1989). [9] Zinoviev, P. A. and Ermakov, Y. N. Energy Dissipation in Composite Materials, Technomic Publishing Company Inc., (1994). 93 [10] Hyer W. M. Stress Analysis of Fiber-Reinforced Composite Materials, WCB McGraw-Hill, (1998). [11] DeGarmo, E. P., Black, J. T. and Kohser, R. A. Materials and Processes in Manufacturing, Prentice Hall, (1997). [12] Gibson, R. F. "Dynamic Mechanical Properties of Advanced Composite Materials and Structures: A Review of Recent Research", The Shock and Vibration Digest, 22 (8), 3-12, (Aug 1990). [13] Gibson, R. F. "Dynamic Mechanical Properties of Advanced Composite Materials and Structures: A Review", The Shock and Vibration Digest, 19 (7), 13-22, (July 1987). [14] Bert, C. W. Proc. 1986 SEM Spring Conf. on Experimental Mechanics, New Orleans, LA, Dynamic Behavior Of Composites: An Overview, 747-751, (1986). [15] Bert, C. W. Appl. Mech. Div. Symp. Ser. ASME, 38, Composite Materials: A Survey of the Damping Capacity of Fiber Reinforced Composites, 53-63, (1980). [16] Chaturvedi, S. K. Encyclopedia of Composites, Lee, S. (Ed.), VCH Publishing Co., NY, Damping of Polymer Matrix Composite Materials, (1989). [17] Suarez, S. A., Gibson, R. F., Sun, C. T. and Chaturvedi, S. K. ?The Influence of Fiber Length and Fiber Orientation on Damping and Stiffness of Polymer Composite Materials?, Experimental Mechanics, 26 (2), 175-184, (1986). [18] Chia, C. Y. ?Geometrically Nonlinear Behavior of Composite Plates?, Applied Mechanics Reviews, 41 (12), 439-451, (1988). [19] Plunkett, R. Proceedings of the IUTAM Symposium, Blacksburg, VA, Damping Mechanisms in Fiber Reinforced Laminates, 93-104, (1983). 94 [20] Yen, S. C. and Cunningham, F. M. Proceedings 1985 SEM Spring Conf. on Experimental Mechanics, Las Vegas, NV, Vibration Characteristics Of Graphite-Epoxy Composite Plates, 60-67 (1985). [21] Gibson, R. F. Principles of Composite Material Mechanics, (McGraw Hill Inc., 1994). [22] Adams, R. D. and Maheri, M. R. ?Dynamic Flexural Properties of Anisotropic Fibrous Composite Beams?. Comp. Sci. Tech., 50 (4), 497-514, (1994). [23] Kenny, J. M. and Marchetti, M. ?Elasto-Plastic Behavior of Thermoplastic Composite Laminates Under Cyclic Loading?, Composite Structures, 32 (1), 375-382, (1995). [24] Hwang, S. J. Ph.D. Thesis, University of Idaho, Characterization of the Effects of Three Dimensional States of Stress on Damping in Composite Laminates, (1988). [25] Gunter, E. J. ?Lund?s Contribution to Rotor Stability: The Indispensable and Fundamental Basis of Modern Compressor Design?, Journal of Vibrations and Acoustics, 125, 462-470, (2003). [26] Vance, J. M. Rotordynamics of Turbomachinery, John Wiley and Sons, Inc., (1988). [27] Nelson, F. C. Sound and Vibration magazine, A brief history of early rotor dynamics, (Jun 2003). [28] Berry, J. E., Technical Associates of Charlotte, PC, Machinery Lubrication Magazine, Oil Whirl and Whip Instabilities - Within Journal Bearings, (May 2005) [29] Alford, J. S. ?Protecting Turbomachinery from Self-Excited Rotor Whirl?, Journal of Engineering for Power, 33 (5), 333?344 (1965). [30] Newkirk, B. L. "Shaft Whipping", General Electric Review, 27, 169, (1924). 95 [31] Kimball, A. L. Jr. ?Internal Friction Theory of Shaft Whirling?, General Electric Review, 27, 244-251, (1924). [32] Gunter, E. J. NASA SP-113, Dynamic Stability of Rotor-Bearing Systems, (1966). [33] Ehrich, F. F. ?Shaft Whirl Induced by Rotor Internal Friction?, ASME Journal of Applied Mechanics, 31, 279-282, (1964). [34] Gunter, E. J. Jr. and Trumpler, P. R. ?The Influence of Internal Friction on the Stability of High Speed Rotors with Anisotropic Supports?, Journal of Engineering for Industry, Nov, 1105-1113, (1969). [35] Lund, J. W. ?Stability and Damped Critical Speeds of a Flexible Rotor in Fluid-Film Bearings?, ASME Journal of Engineering for Industry, 96 (2), 509-517, (1974). [36] Bently, D. E., Muzsynska, A. NASA Conference Publication, Instability in Rotating Machinery, NAS 1-55-2409, Rotor Internal Friction Instability, 337-348, (1982). [37] Wettergren, H. L. ?Material Damping in Composite Rotors?, Journal of Composite Materials, 32 (7), 652-663 (1998). [38] Chen, L.T. Doctoral Thesis, University of Oklahoma, Whirling Response and Stability of Flexibly Mounted, Ring-Type Flywheel Systems, (1979). [39] Genta, G. Kinetic Energy Storage, Butterworths, London, (1985). [40] Columbia University Press, Columbia Electronic Encyclopedia, flywheel (in mechanics). [41] Gowayed, Y., Abel-Hady, F., Flowers, G. T., and Trudell, J. J. ?Optimal Design of Multi-Direction Composite Flywheel Rotors?, Polymer Composites, 23 (3), 433-441, (2002). 96 [42] Jansen, R. H., Hervol, D. S., Dever, T. P., Anzalone, S. M., Trudell, J. J., and Kenny, A. NASA/TM-2002-211788, Redesign of Glenn Research Center D1 Flywheel Module, (2002). [43] Kinra, V. K., Ray, S., Zhu, C., Friend, R. D. and Rawal S. P. ?Measurement of Nonlinear Damping in a Gr/Al Metal-Matrix Composite?, Experimental Mechanics, 37 (1), 5-10, (Mar 1997). [44] Ewins, D. J. Modal Testing: Theory, Practice and Application, Research Studies Press LTD, Baldock, England, (2000). [45] Sun, C. T., Wu, J. K. and Gibson R. F. ?Prediction of material damping of laminated polymer matrix composites?, Journal of Materials Science, 22 (3), 1006-1012, (1987). [46] Ungar, E. E and Kerwin, E. M Jr. ?Loss Factors of Viscoelastic Systems in Terms of Energy Concepts?, J. Acoust. Soc. Am, 34 (7), 954-957, (July 1962). [47] Chen, J., Gowayed, Y., Moreira, A. and Flowers, G. ?Damping of Polymer Composite Materials for Flywheel Applications?, Polymer composites, 26, (2005). [48] Moreira, A., Flowers G. T., Matras, A., Balas, M. and Fausz J. IDETC/CIE 2005, The Influence of Internal Damping on the Rotordynamic Stability of a Flywheel Rotor with Flexible Hub, (2005). [49] Moreira, A., Flowers, G. T. and Gowayed, Y. ICSV 10, Effects of Internal Damping on the Dynamic Stability of High-Speed Composite Rotor Systems, (2003). [50] Wettergren, H. L. Doctoral Thesis, Linkoping University, Rotordynamic Analysis with Special Reference to Composite Rotors and Internal Damping, (1996). 97 [51] Wallace, M. and Bert, C. Proceedings of the Oklahoma Academy of Science, 59, Experimental Determination of Dynamic Young's Modulus and Damping of an Aramid- Fabric/Polyester Composite Material, 98-101, (1979). [52] Cloud, C. H., Maslen, E. H. and Barret, L. E. Proc. IMechE, Evaluation of Damping Ratio Estimation Techniques for Rotordynamic Stability Measurements, 541-550, (2004). [53] Kinra, V. K. and Wolfenden, A. ASTM/STP 1169, American Society for Testing and Materials, Philadelphia, PA, M3D: Mechanics and Mechanisms of Material Damping, (1992). [54] Gunter, E. J. ?The Influence of Internal Friction on the Stability of High Speed Rotors?, Journal of Engineering for Industry, Nov (8), 683-688, (1967). 98 APPENDIX Computer Codes Computer codes used for this work are printed below. All of them were developed in Matlab 7.01. Code to generate Figure 3.18 and Figure 3.19 % This program shows how the damping at high amplitudes is much higher than % that at low ones. In this case for d2, peaks of HP filtered d1, the % damping ratio at low frequencies is .241%, while at high frequencies it % is .335%. This represents an increase of 38.8%. clear i j pktime peak fs fn s pkt expfitdecay expfitdecay2 % d2 is d1 hp filtered and restarted(from peak1) load ss2d593; sig=data593; sigl=length(sig); fs=132300; fn=593; frstpk=15; zh=.00167; zl=.000966; PkSCAL=.4; % High amplitude fit i=1; for j=4: sigl-2 if sig(j)>sig(j-1)&sig(j)>sig(j- 2)&sig(j)>sig(j+1)&sig(j)>sig(j+2)&sig(j)>sig(j-3)&sig(j)>sig(j- 4)&sig(j)>sig(j+3)&sig(j)>sig(j+4)... sig(j)>sig(j-5)&sig(j)>sig(j- 6)&sig(j)>sig(j+5)&sig(j)>sig(j+6)&sig(j)>sig(j-7)&sig(j)>sig(j- 8)&sig(j)>sig(j+7)&sig(j)>sig(j+8); pktime(i)=j/fs; 99 peak(i)=sig(j); i=i+1; end end pktime=pktime-pktime(1); pkl=length(peak); expfitdecayH=[peak(frstpk)*exp(-zh*2*pi*fn.*(pktime-pktime(frstpk)))]; % Low amplitude fit expfitdecay2=.69*peak1*exp(-.00242)... clear pktime peak i j i=1; for j=4: sigl-2 if sig(j)>sig(j-1)&sig(j)>sig(j- 2)&sig(j)>sig(j+1)&sig(j)>sig(j+2)&sig(j)>sig(j-3)&sig(j)>sig(j- 4)&sig(j)>sig(j+3)&sig(j)>sig(j+4)... sig(j)>sig(j-5)&sig(j)>sig(j- 6)&sig(j)>sig(j+5)&sig(j)>sig(j+6)&sig(j)>sig(j-7)&sig(j)>sig(j- 8)&sig(j)>sig(j+7)&sig(j)>sig(j+8); pktime(i)=j/fs; peak(i)=sig(j); i=i+1; end end pktime=pktime-pktime(frstpk); figure(67) expfitdecayL=[(peak(frstpk)*PkSCAL)*exp(-zl*2*pi*fn.*pktime)]; plot(pktime(frstpk:pkl),peak(frstpk:pkl),'k.-') hold on plot(pktime(frstpk:pkl),expfitdecayH(frstpk:pkl),'b') plot(pktime(frstpk:pkl),expfitdecayL(frstpk:pkl),'r','LineStyle','--') xlabel('time (s)','FontSize',12) ylabel('displacement (m)','FontSize',12) legend('experimental',strcat('\zeta =',num2str(zh)),strcat('\zeta = ',num2str(zl))) hold off figure(68) semilogy(pktime(frstpk:pkl),peak(frstpk:pkl),'k.-') hold on semilogy(pktime(frstpk:pkl),expfitdecayH(frstpk:pkl),'b') semilogy(pktime(frstpk:pkl),expfitdecayL(frstpk:pkl),'r','LineStyle','- -') % axis([0 1.4 7e-9 1.1e-2]) xlabel('time (s)','FontSize',12) ylabel('displacement (m)','FontSize',12) legend('experimental',strcat('\zeta =',num2str(zh)),strcat('\zeta = ',num2str(zl)),3) hold off 100 Frequency and amplitude dependence of damping in samples loaded axially % Frequency and amplitude dependence of damping in samples loaded axially. clc clear p leg format long ii=1; figure(24) hold all range=30; for amp=(4:.4:6.4)*1e-5; load ss2d524 sig=data524; fn=524; fs=132300; jump=floor(fs/fn); sigl=length(sig); sigtime=0:1/fs:(sigl-1)/fs'; i=1; j=8; while jsig(j-1)&&sig(j)>sig(j- 2)&&sig(j)>sig(j+1)&&sig(j)>sig(j+2)... &&sig(j)>sig(j-3)&&sig(j)>sig(j- 4)&&sig(j)>sig(j+3)&&sig(j)>sig(j+4)... &&sig(j)>sig(j-5)&&sig(j)>sig(j- 6)&&sig(j)>sig(j+5)&&sig(j)>sig(j+6)... &&sig(j)>sig(j-7)&&sig(j)>sig(j- 8)&&sig(j)>sig(j+7)&&sig(j)>sig(j+8); pktime(i)=j/fs; peak(i)=sig(j); if abs(peak(i)-amp)<3e-7 ttt=i; end i=i+1; j=j+jump-20; continue end j=j+1; end pktime=pktime'; peak=peak'; pkl=length(peak); amp524(ii)=peak(ttt-10); zeta524(ii)=log(peak(ttt-range)/peak(ttt))/2/pi/fn/(pktime(ttt)- pktime(ttt-range)); clear pktime peak ttt 101 load ss2d593 sig=data593; fn=593; fs=132300; jump=floor(fs/fn); sigl=length(sig); sigtime=0:1/fs:(sigl-1)/fs'; i=1; j=8; while jsig(j-1)&&sig(j)>sig(j- 2)&&sig(j)>sig(j+1)&&sig(j)>sig(j+2)... &&sig(j)>sig(j-3)&&sig(j)>sig(j- 4)&&sig(j)>sig(j+3)&&sig(j)>sig(j+4)... &&sig(j)>sig(j-5)&&sig(j)>sig(j- 6)&&sig(j)>sig(j+5)&&sig(j)>sig(j+6)... &&sig(j)>sig(j-7)&&sig(j)>sig(j- 8)&&sig(j)>sig(j+7)&&sig(j)>sig(j+8); pktime(i)=j/fs; peak(i)=sig(j); if abs(peak(i)-amp)<3e-7 ttt=i; end i=i+1; j=j+jump-20; continue end j=j+1; end pktime=pktime'; peak=peak'; pkl=length(peak); amp593(ii)=peak(ttt-10); zeta593(ii)=log(peak(ttt-range)/peak(ttt))/2/pi/fn/(pktime(ttt)- pktime(ttt-range)); clear pktime peak ttt load ss2d677 sig=data677; fn=677; fs=132300; jump=floor(fs/fn); sigl=length(sig); sigtime=0:1/fs:(sigl-1)/fs'; i=1; j=8; while jsig(j-1)&&sig(j)>sig(j- 2)&&sig(j)>sig(j+1)&&sig(j)>sig(j+2)... 102 &&sig(j)>sig(j-3)&&sig(j)>sig(j- 4)&&sig(j)>sig(j+3)&&sig(j)>sig(j+4)... &&sig(j)>sig(j-5)&&sig(j)>sig(j- 6)&&sig(j)>sig(j+5)&&sig(j)>sig(j+6)... &&sig(j)>sig(j-7)&&sig(j)>sig(j- 8)&&sig(j)>sig(j+7)&&sig(j)>sig(j+8); pktime(i)=j/fs; peak(i)=sig(j); if abs(peak(i)-amp)<3e-7 ttt=i; end i=i+1; j=j+jump-20; continue end j=j+1; end pktime=pktime'; peak=peak'; pkl=length(peak); amp677(ii)=peak(ttt-floor(range/2)); zeta677(ii)=log(peak(ttt-range)/peak(ttt))/2/pi/fn/(pktime(ttt)- pktime(ttt-range)); clear pktime peak ttt load ss2d735 sig=data735; fn=735; fs=132300; jump=floor(fs/fn); sigl=length(sig); sigtime=0:1/fs:(sigl-1)/fs'; i=1; j=8; while jsig(j-1)&&sig(j)>sig(j- 2)&&sig(j)>sig(j+1)&&sig(j)>sig(j+2)... &&sig(j)>sig(j-3)&&sig(j)>sig(j- 4)&&sig(j)>sig(j+3)&&sig(j)>sig(j+4)... &&sig(j)>sig(j-5)&&sig(j)>sig(j- 6)&&sig(j)>sig(j+5)&&sig(j)>sig(j+6)... &&sig(j)>sig(j-7)&&sig(j)>sig(j- 8)&&sig(j)>sig(j+7)&&sig(j)>sig(j+8); pktime(i)=j/fs; peak(i)=sig(j); if abs(peak(i)-amp)<3e-7 ttt=i; end i=i+1; j=j+jump-20; continue 103 end j=j+1; end pktime=pktime'; peak=peak'; pkl=length(peak); amp735(ii)=peak(ttt-floor(range/2)); zeta735(ii)=log(peak(ttt-range)/peak(ttt))/2/pi/fn/(pktime(ttt)- pktime(ttt-range)); clear pktime peak ttt load ss2d774 sig=data774; fn=774; fs=132300; jump=floor(fs/fn); sigl=length(sig); sigtime=0:1/fs:(sigl-1)/fs'; i=1; j=8; while jsig(j-1)&&sig(j)>sig(j- 2)&&sig(j)>sig(j+1)&&sig(j)>sig(j+2)... &&sig(j)>sig(j-3)&&sig(j)>sig(j- 4)&&sig(j)>sig(j+3)&&sig(j)>sig(j+4)... &&sig(j)>sig(j-5)&&sig(j)>sig(j- 6)&&sig(j)>sig(j+5)&&sig(j)>sig(j+6)... &&sig(j)>sig(j-7)&&sig(j)>sig(j- 8)&&sig(j)>sig(j+7)&&sig(j)>sig(j+8); pktime(i)=j/fs; peak(i)=sig(j); if abs(peak(i)-amp)<3e-7 ttt=i; end i=i+1; j=j+jump-20; continue end j=j+1; end pktime=pktime'; peak=peak'; pkl=length(peak); amp774(ii)=peak(ttt-floor(range/2)); zeta774(ii)=log(peak(ttt-range)/peak(ttt))/2/pi/fn/(pktime(ttt)- pktime(ttt-range)); clear pktime peak ttt plot([593 677 735 774],[zeta593(ii),zeta677(ii),zeta735(ii),zeta774(ii)],'.') 104 p(ii,:)=polyfit([593 677 735 774],[zeta593(ii),zeta677(ii),zeta735(ii),zeta774(ii)],1); leg(ii)=amp*1e3; ii=ii+1; end pfreq=[sum(p(:,1))/7,sum(p(:,2))/7]; xlabel('Frequency (Hz)') ylabel('Damping ratio \zeta') title(strcat('Damping ratio \zeta(x,f) vs frequency for various amplitudes of vibration')) legend(strcat(num2str(leg(1)),' mm'),strcat(num2str(leg(2)),' mm'),strcat(num2str(leg(3)),' mm'),strcat(num2str(leg(4)),' mm'),strcat(num2str(leg(5)),' mm'),strcat(num2str(leg(6)),' mm'),strcat(num2str(leg(7)),' mm')) for ii=1:length(p) plot((550:850),polyval(p(ii,:),(550:850))) end % AXIS([560 810 8e-4 1.05e-3]) hold off pamp593=polyfit(amp593,zeta593,1) pamp677=polyfit(amp677,zeta677,1) pamp735=polyfit(amp735,zeta735,1) pamp774=polyfit(amp774,zeta774,1) zeta=[zeta593' zeta677' zeta735' zeta774']; %experimental zetas %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% pamp1fit=polyfit(2*pi*[593 677 735 774],[pamp593(1) pamp677(1) pamp735(1) pamp774(1)],1) pamp2fit=polyfit(2*pi*[593 677 735 774],[pamp593(2) pamp677(2) pamp735(2) pamp774(2)],1) figure(98) plot(amp593,zeta593,'.g')%,'Color',[.6 .6 .6]) hold on plot(amp677,zeta677,'.b')%,'Color',[.4 .4 .4]) plot(amp735,zeta735,'.r')%,'Color',[.2 .2 .2]) plot(amp774,zeta774,'.k')%,'Color',[0 0 0]) legend('593Hz','677Hz','735Hz','774Hz',4) plot(3e-5:.1e-6:7.9e-5,polyval(pamp593,3e-5:.1e-6:7.9e- 5),'g')%,'Color',[.5 .5 .5]) plot(3e-5:.1e-6:7.9e-5,polyval(pamp677,3e-5:.1e-6:7.9e- 5),'b')%,'Color',[.3 .3 .3]) plot(3e-5:.1e-6:7.9e-5,polyval(pamp735,3e-5:.1e-6:7.9e- 5),'r')%,'Color',[.2 .2 .2]) plot(3e-5:.1e-6:7.9e-5,polyval(pamp774,3e-5:.1e-6:7.9e- 5),'k')%,'Color',[0 0 0]) xlabel('Vibration amplitude (m)') ylabel('Damping ratio \zeta(x,f)') title('Damping ratio vs amplitude for different 1st natural frequencies') hold off 105 Eigenvalue analysis of matrix A from Eq. (4.33) In this Matlab code, the parameters for the flywheel model with 8 degrees of freedom (2 rotational and 2 translational for each, hub and rim) are specified and then the mass, stiffness, damping and gyroscopic matrices are formed. These are all expressed in the state space form in matrix A (AA in the program), for which an eigenvalue analysis is performed using a range of running speeds, to obtain the stable bounds of operation at each speed under steady state. The code allows the user to incorporate or not the physical characteristics of a rigid shaft, which add up to the ones of the hub. %BASIC FLYWHEEL MODEL WITH RIGID SHAFT (OR SIMPLE MASSLESS SHAFT) %THIS MODEL ASSUMES NO MASS IN THE INTERCONNECTION BETWEEN THE HUB %AND THE RIM %THIS MODEL ASSUMES LINEAR DAMPING clc tcpu=cputime; %--RIM PROPERTIES roR=.16; %m riR=.08; %m wR=.06; %m rhoC=1400; %kg/m^3 comp (6.8e-3/.01165/.0033/.12367) mR=pi*(roR^2-riR^2)*wR*rhoC; %kg ItR=1/12*mR*(3*(roR^2+riR^2)+wR^2); %kg m^2 IpR=1/2*mR*(roR^2+riR^2); %kg m^2 IRratio=IpR/ItR; %--HUB PROPERTIES roH=.06; %m riH=0; %m wH=.07; %m rhoAl=2700; %kg/m^3 Al at 20 deg C %wo/shaft mH=pi*(roH^2-riH^2)*wH*rhoAl; %kg mH=2.3 ItH=1/12*mH*(3*(roH^2+riH^2)+wH^2); %kg m^2 IpH=1/2*mH*(roH^2+riH^2); %kg m^2 %w/shaft 106 %--SHAFT PROPERTIES % roS=.02; %m % riS=0; %m % wS=.4; %m % mS=pi*(roS^2)*wS*rhoAl; %kg mH=2.3 % ItS=1/12*mS*(3*(roS^2)+wS^2); % IpS=1/2*mS*(roS^2); % % mH=pi*(roH^2-riH^2)*wH*rhoAl; %kg mH=2.3 % ItH=1/12*mH*(3*(roH^2+riH^2)+wH^2)+ItS; %kg m^2 % IpH=1/2*mH*(roH^2+riH^2)+IpS; %kg m^2 % mH=mH+mS; IHratio=IpH/ItH; mT=mH+mR; ItT=ItH+ItR; %--BEARING SUPPORT STIFFNESS kBTHETA=ItT*(5000*2*pi/60)^2; %wBTHETA=523 Hz kBX=mT*(5000*2*pi/60)^2; %--INTERCONNECTION STIFFNESS kTHETA=3*kBTHETA; kX=3*kBX; %--BEARING SUPPORT DAMPING %----DAMPING RATIOS zBTHETA=0.02; zBX=0.02; %----DAMPING CONSTANTS cBTHETA=ItH*(2*zBTHETA*sqrt(kBTHETA/ItH)); cBX=mH*(2*zBX*sqrt(kBX/mH)); %--INTERCONNECTION DAMPING %----DAMPING RATIOS zTHETA=.0015; %55*1e-5; zX=.0015; %55*1e-5; %----DAMPING CONSTANTS cTHETA=ItR*(2*zTHETA*sqrt(kTHETA/ItR)); cX=2*zX*sqrt(kX/mR)*mR; %--FORM MASS, STIFFNESS, DAMPING, AND GYROSCOPIC MATRICES %----MASS MROT=[ItR,0,0,0;0,ItR,0,0;0,0,ItH,0;0,0,0,ItH]; MTRAN=[mR,0,0,0;0,mR,0,0;0,0,mH,0;0,0,0,mH]; %----DAMPING CROT=[cTHETA,0,-cTHETA,0;0,cTHETA,0,-cTHETA;- cTHETA,0,cTHETA+cBTHETA,0;... 0,-cTHETA,0,cTHETA+cBTHETA]; CTRAN=[cX,0,-cX,0;0,cX,0,-cX;-cX,0,cX+cBX,0;0,-cX,0,cX+cBX]; %-----STIFFNESS KROT=[kTHETA,0,-kTHETA,0;0,kTHETA,0,-kTHETA;- kTHETA,0,kTHETA+kBTHETA,0;... 0,-kTHETA,0,kTHETA+kBTHETA]; KTRAN=[kX,0,-kX,0;0,kX,0,-kX;-kX,0,kX+kBX,0;0,-kX,0,kX+kBX]; %----FORM SPEED DEPENDENT TERMS 107 GYRO=[0,IpR,0,0;-IpR,0,0,0;0,0,0,IpH;0,0,-IpH,0]; KRCROSS=[0,cTHETA,0,-cTHETA;-cTHETA,0,cTHETA,0;0,-cTHETA,0,cTHETA;... cTHETA,0,-cTHETA,0]; KTCROSS=[0,cX,0,-cX;-cX,0,cX,0;0,-cX,0,cX;cX,0,-cX,0]; %--RUNNING SPEED w0=20000/60*2*pi; %20000 RPM dw=200/60*2*pi; ww=zeros(1,2000); for i=1:4000 w=0+dw*(i-1); %--FORM TOTAL MASS, STIFFNESS, AND DAMPING MATRICES MTOT=[MROT,0*eye(4,4);0*eye(4,4),MTRAN]; CTOT=[CROT+w*GYRO,0*eye(4,4);0*eye(4,4),CTRAN]; KTOT=[KROT+w*KRCROSS,0*eye(4,4);0*eye(4,4),KTRAN+w*KTCROSS]; %--FORM STATE MATRIX MI=inv(MTOT); AA=[-MI*CTOT,-MI*KTOT;eye(8,8),0*eye(8,8)]; [v,d]=eig(AA); ww(i)=w; dr(i,1:16)=real(diag(d)'); di(i,1:16)=abs(imag(diag(d)')); end figure(1) plot(ww/2/pi*60,dr(:,:),'.') % axis([0 8e4 -10 10]); title('Real part') xlabel('Running speed (RPM)') figure(2) plot(ww/2/pi*60,di(:,:)/2/pi,'.') % axis([0 3e4 0 700]); hold on plot(ww/2/pi*60,ww/2/pi,'k.') hold off xlabel('Running speed (RPM)') ylabel('Frequency (Hz)') title('Imaginary part') tcpu=cputime-tcpu; Stable thresholds for ranges of values of k H and c H This code is the one used to vary the parameters k H and c H and generate Figure 4.7: 108 % Finds the running speed and hub stiffness at which % the destabilizing eigenvalue changes for each % damping ratio value % zeta: damping ratio % kh: hub stiffness % Weigchange: w value at which eigenvalue change occurs % Keigchange: kh value at which eigenvalue change occurs % Change: (kk x 2) matrix containing the eigenvalue change % |zsd | |zs | % |zsdd| = [A] |zsd| zs = xs + i ys z shaft % |zrd | |zr | % |zrdd| |zrd| zr = xr + i yr z rim % mt=ms+mr;a=ms/mt;wb=sqrt(kb/mt); clc clear time=cputime; wb=53; mt=10;a=.3; ch=165.8; cb=84; %rad/s , kg %for kk=1:29 kk=19; zeta(kk)=kk*.001+.001; for jj=1:201 kh(jj)=1+500*(jj-1); wh=sqrt(kh(jj)/((1-a)*mt)); ch=2*zeta(kk)*((1-a)*mt)*wh; for ii=1:80000; w=ii/2; W(jj,kk)=w*60/2/pi; A=[0 1 0 0;-(wb^2+wh^2*(1-a)-i*w*ch/mt)/a -(cb+ch)/a/mt - (-wh^2*(1-a)+i*w*ch/mt)/a ch/a/mt; 0 0 0 1;wh^2-i*w*ch/mt/(1-a) ch/(1-a)/mt -(wh^2- i*w*ch/mt/(1-a)) -ch/(1-a)/mt]; L(ii,1:4)=eig(A)'; if real(L(ii,1))>.0000001 W1(jj,kk)=w*60/2/pi; %RPM elseif real(L(ii,2))>0.0000001 W2(jj,kk)=w*60/2/pi; %RPM elseif real(L(ii,3))>0.0000001 W3(jj,kk)=w*60/2/pi; %RPM elseif real(L(ii,4))>0.0000001 W4(jj,kk)=w*60/2/pi; %RPM end end end %end figure(1) %set(gcf,'DefaultAxesColorOrder',CO) hold on plot(kh,W1,'k') plot(kh,W4,'k') 109 plot(kh,W3,'k') plot(kh,W2,'k') title('max stable running speed for zeta: .002 -> .02') xlabel('k_h (Kg/s^2)') ylabel('running speed (RPM)') figure(3) mesh(kh,zeta,W') xlabel('k_h (kg/s^2)') ylabel('zeta') zlabel('running speed (RPM)') figure(2) plot(zeta,Weigchange,'k') title('Maximum stable running speed for each zeta') xlabel('\zeta') ylabel('running speed (RPM)') Weigchange=Weigchange'; Keigchange=Keigchange'; zeta=zeta'; Stable thresholds for ranges of values of k B and c B This code is the one used to vary the parameters k B and c B and generate Figure 4.10: % cb changes \omega % Finds the running speed and hub stiffness at which the destabilizing % eigenvalue changes for each damping ratio value % zeta: damping ratio % kh: hub stiffness % Weigchange: w value at which eigenvalue change occurs % Keigchange: kh value at which eigenvalue change occurs % Change: (kk x 2) matrix containing the eigenvalue change % |zsd | |zs | % |zsdd| = [A] |zsd| zs = xs + i ys z shaft % |zrd | |zr | % |zrdd| |zrd| zr = xr + i yr z rim % mt=ms+mr;a=ms/mt;wb=sqrt(kb/mt); clc clear tiempo=cputime; % cb= 2 zeta wb mt 110 % zeta= cb / (2 wb mt) % cb=84; % wb= 53; wh= 628; ch= 165.8; mt= 10; a= .3; %rad/s , kg zeta=zeros(1,18); W=zeros(100,18); for kk=1:18 %18 zeta(kk)=kk*.005+.01; %(.015 ---> .1) kb=zeros(1,100); for jj=1:100 %100 wb=jj; kb(jj)=wb^2*mt; cb=2*zeta(kk)*mt*wb; for ii=1:80000; %80000 w=ii/2; W(jj,kk)=w*60/2/pi; A=[0 1 0 0;-(wb^2+wh^2*(1-a)-i*w*ch/mt)/a -(cb+ch)/a/mt -(- wh^2*(1-a)+i*w*ch/mt)/a ch/a/mt; 0 0 0 1;wh^2-i*w*ch/mt/(1-a) ch/(1-a)/mt -(wh^2- i*w*ch/mt/(1-a)) -ch/(1-a)/mt]; L(ii,1:4)=eig(A)'; if real(L(ii,1))>1e-8 % W1(jj,kk)=w*60/2/pi; %RPM break elseif real(L(ii,2))>1e-8 % W2(jj,kk)=w*60/2/pi; %RPM break elseif real(L(ii,3))>1e-8 % W3(jj,kk)=w*60/2/pi; %RPM break elseif real(L(ii,4))>1e-8 % W4(jj,kk)=w*60/2/pi; %RPM break end end end end % kh=wh^2*((1-a)*mt) % figure(1) % %set(gcf,'DefaultAxesColorOrder',CO) % hold on % plot(kh,W1,'k') % plot(kh,W4,'k') % plot(kh,W3,'k') % plot(kh,W2,'k') % title('max stable running speed for zeta: .002 -> .02') % xlabel('k_h (Kg/s^2)') % ylabel('running speed (RPM)') figure(3) mesh(kb,zeta,W') xlabel('bearing stiffness (Kg/s^2)') 111 ylabel('bearing damping ratio, \zeta_ _b') zlabel('running speed (RPM)') time=cputime-time; Derivation of Li?nard-Chipart conditions For the derivation of the Li?nard-Chipart conditions for the stability of a purely translational model with 4 degrees of freedom of the rotor with flexible hub-rim interface on rigid shaft mounted on flexible bearings, it was necessary to form the matrices, extract the determinant of the matrix ( )?I - A , and then select and group terms of equal degree in ? from the resulting polynomial into a0, a1, ?, a8. Then it was simple to form the matrices required by the procedure to obtain the Li?nard-Chipart criteria. These criteria are a simplification (which reduce computation time) of the probably more popular Routh-Hurwitz criteria. The expressions for the output conditions 1 to 8 alone use around 70 pages, so they will not be listed here. Although they are so extensive, the calculations of these algebraic expressions for each run with a new set of parameters (or at least one or two parameters changing) is less demanding computationally than extracting eigenvalues each time. The only drawback is the lack of an output of eigenvectors, which give direct information about the modes being involved in the instability, but it is always possible to return to the eigenvalue analysis to study interesting phenomena that may appear to be occurring. The program to perform what has been described above is given below: % Form Li?nard-Chipart criteria for stability of 4 % translational degrees of freedom model of rotor 112 % with flexible hub-rim interface on rigid shaft % mounted on flexible bearings. clear clc syms mR mH cX cBX kX kBX w L %--FORM MASS, STIFFNESS, DAMPING, AND GYROSCOPIC MATRICES %----MASS MTRAN=[mR,0,0,0;0,mR,0,0;0,0,mH,0;0,0,0,mH]; %----DAMPING CTRAN=[cX,0,-cX,0;0,cX,0,-cX;-cX,0,cX+cBX,0;0,-cX,0,cX+cBX]; %-----STIFFNESS KTRAN=[kX,0,-kX,0;0,kX,0,-kX;-kX,0,kX+kBX,0;0,-kX,0,kX+kBX]; %----FORM SPEED DEPENDENT TERMS KTCROSS=[0,cX,0,-cX;-cX,0,cX,0;0,-cX,0,cX;cX,0,-cX,0]; %--FORM TOTAL MASS, STIFFNESS, AND DAMPING MATRICES MTOT=MTRAN; CTOT=CTRAN; KTOT=KTRAN+w*KTCROSS; %--FORM STATE MATRIX MI=inv(MTOT); AA=[-MI*CTOT,-MI*KTOT;eye(4,4),0*eye(4,4)]; pL=simplify(det(L*eye(8,8)-AA)) a0=1 a1=simplify((2*mR*mH^2*cX+2*mR^2*mH*cBX+2*mR^2*mH*cX)/mR^2/mH^2) a2=simplify((2*mR*mH^2*kX+2*mR*cX^2*mH+4*mR*mH*cX*cBX+2*mR^2*cBX*cX+2*m R^2*mH*kBX+2*mR^2*mH*kX+mR^2*cBX^2+mR^2*cX^2+cX^2*mH^2)/mR^2/mH^2) a3=simplify((2*mR*cX^2*cBX+2*mR*cBX^2*cX+4*mR*cX*mH*kX+4*mR*mH*cX*kBX+4 *mR*cBX*kX*mH+2*cX^2*mH*cBX+2*mH^2*kX*cX+2*mR^2*cX*kBX+2*mR^2*cX*kX+2*m R^2*cBX*kBX+2*mR^2*cBX*kX)/mR^2/mH^2) a4=simplify((2*w^2*cX^2*mR*mH+2*mR*mH*kX^2+2*mR*cBX^2*kX+w^2*cX^2*mR^2+ 2*mR*cX^2*kBX+cX^2*cBX^2+mR^2*kBX^2+mH^2*kX^2+4*mR*cX*cBX*kX+4*mR*cBX*c X*kBX+4*mR*mH*kX*kBX+cX^2*mH^2*w^2+2*cX^2*mH*kBX+4*cBX*kX*cX*mH+2*mR^2* kBX*kX+mR^2*kX^2)/mR^2/mH^2) a5=simplify((2*w^2*cX^2*mR*cBX+2*mR*kBX^2*cX+2*cX^2*cBX*kBX+4*mR*cX*kBX *kX+4*mR*cBX*kX*kBX+2*mR*cBX*kX^2+2*mH*kX^2*cBX+2*cX^2*mH*cBX*w^2+4*mH* kX*cX*kBX+2*cBX^2*kX*cX)/mR^2/mH^2) a6=simplify((2*w^2*cX^2*mR*kBX+cBX^2*kX^2+2*mR*kBX^2*kX+2*mR*kBX*kX^2+w ^2*cX^2*cBX^2+2*kBX*kX^2*mH+cX^2*kBX^2+2*w^2*cX^2*mH*kBX+4*cBX*kX*cX*kB X)/mR^2/mH^2) a7=simplify((2*cX^2*cBX*kBX*w^2+2*kBX^2*kX*cX+2*kBX*kX^2*cBX)/mR^2/mH^2) a8=simplify((cX^2*kBX^2*w^2+kBX^2*kX^2)/mR^2/mH^2) %Lienard-Chipart Criteria Cond8=a8 %>0 Cond7=simplify(det([a1 a3 a5 a7 0 0 0;a0 a2 a4 a6 0 0 0;0 a1 a3 a5 a7 0 0;0 a0 a2 a4 a6 0 0;0 0 a1 a3 a5 a7 0;0 0 a0 a2 a4 a6 0;0 0 0 a1 a3 a5 a7])) 113 Cond6=a6 %>0 Cond5=simplify(det([a1 a3 a5 0 0;a0 a2 a4 0 0;0 a1 a3 a5 0;0 a0 a2 a4 0;0 0 a1 a3 a5])) Cond4=a4 %>0 Cond3=simplify(det([a1 a3 0;a0 a2 0;0 a1 a3])) Cond2=a2 %>0 Cond1=a1 clear mR mH cX cBX kX kBX w L a0 a1 a2 a3 a4 a5 a6 a7 a8