HEALTH MONITORING FOR FLYWHEEL ROTORS SUPPORTED BY ACTIVE MAGNETIC BEARINGS Except where reference is made to the work of others, the work described in this thesis is my own work or was done in collaboration with my advisory committee. This thesis does not include proprietary or classified information. _____________________________________________ Kelly Lowell Barber ______________________________ Subhash C. Sinha Professor Mechanical Engineering ______________________________ Jeffrey C. Suhling Quina Distinguished Professor Mechanical Engineering ______________________________ Jerry Fausz Program Manager Air Force Research Laboratory Kirtland AFB Albuquerque, NM ______________________________ George T. Flowers, Chair Professor Mechanical Engineering ______________________________ Robert L. Jackson Assistant Professor Mechanical Engineering ______________________________ Stephen L. McFarland Dean Graduate School HEALTH MONITORING FOR FLYWHEEL ROTORS SUPPORTED BY ACTIVE MAGNETIC BEARINGS Kelly Lowell Barber A Thesis Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Master of Science Auburn, Alabama August 7, 2006 iii HEALTH MONITORING FOR FLYWHEEL ROTORS SUPPORTED BY ACTIVE MAGNETIC BEARINGS Kelly Lowell Barber Permission is granted to Auburn University to make copies of this thesis at its discretion, upon the request of the individuals or institutions and at their expense. The author reserves all publication rights. ______________________________ Signature of Author ______________________________ Date of Graduation iv VITA Kelly Lowell Barber, son of Greg and Betty Barber, was born on August 10, 1981 in Union City, Tennessee. Graduating with honors from Union City High School in the spring of 1999, Kelly went on to pursue a Bachelor of Science degree in Engineering from the University of Tennessee at Martin. In the spring of 2004, he graduated Magna cum Laude and was inducted into the Order of the Engineer. During the summer of 2004 Kelly married Jennie A. Aden, daughter of Bill and Jo Margaret Aden, and later in the fall of 2004, they entered graduate school together at Auburn University. Kelly spent the summer of 2005 interning with the Air Force Research Laboratory at Kirtland Air Force Base in Albuquerque, New Mexico. v THESIS ABSTRACT HEALTH MONITORING FOR FLYWHEEL ROTORS SUPPORTED BY ACTIVE MAGNETIC BEARINGS Kelly Lowell Barber Master of Science, August 7, 2006 (B.S.E., University of Tennessee at Martin, 2004) 217 Typed Pages Directed by George T. Flowers Magnetic bearings are not a new technology in themselves, yet the control and implementation of such devices is still a budding science. Magnetic bearings are particularly attractive for some applications because of the low energy loss associated with their frictionless operation. By using electromagnetic forces, a rotor can be levitated without mechanical contact between the rotor and its supports. One major application of magnetic bearings is mechanical energy storage devices using flywheels, which potentially have a substantially higher energy storage density than more standard devices such as chemical batteries. Flywheels equipped with magnetic bearings can not only be used to store energy, but can also provide actuation for attitude control in satellite applications. The conceptual designs for such flywheel systems vi generally employ a disk consisting of a hub (either metal or composite) and a high strength composite rim spinning at high rotation speeds, 50kRPM or higher. This causes substantial stresses to be applied to the rim and hub. In addition, as energy is added or withdrawn from the system, the rotor speed changes over a wide range, resulting in cyclic stresses on the disk and possible fatigue induced cracks. The initiation and growth of such cracks has potentially disastrous implications, possibly causing the entire structure to be destroyed. Accordingly, health monitoring is of critical importance in maintaining the integrity of such devices. In this thesis, a health monitoring strategy based upon the acquisition and analysis of vibration measurements is described and evaluated. A common technique in this regard is to track changes in the synchronous vibration due to imbalance. However, such an approach must consider the controller strategy used with the magnetic bearings. Herein, a simulation model is developed that consists of a flywheel system supported by magnetic bearings, which are controlled using an adaptive strategy that suppresses synchronous vibration. The interaction between the rotor vibration and the controller responses are evaluated in order to provide insight into indicators of crack initiation and growth. The results and conclusions are also validated using an experimental test rig. Some insights and guidelines as to appropriate strategies for crack detection in rotor systems interacting with active bearing controllers are presented and discussed. vii ACKNOWLEDGEMENTS The author would like to acknowledge and thank George T. Flowers for his assistance and direction in the process of this thesis and for the endless presentation of opportunity. The author would also like to acknowledge his committee which includes Subhash Sinha, Jeffrey Suhling, and Robert Jackson from Auburn University, and Jerry Fausz from the Air Force Research Laboratory. A special acknowledgement and thank you is made to Nathan Martz and Brian Wilson at Kirtland Air Force Base and Alex Matras for all of their invaluable help and support. Finally, the author would like to thank his wife, Jennie, and the rest of his family for their unwavering support. viii Style manual or journal used: Journal of ASME_________________________________ Computer software used: Microsoft Word XP___________________________________ ix TABLE OF CONTENTS LIST OF FIGURES............................................................................................... LIST OF TABLES................................................................................................ NOMENCLATURE.............................................................................................. CHAPTER 1 ? INTRODUCTION........................................................................ CHAPTER 2 ? LITERATURE REVIEW............................................................. CHAPTER 3 ? EQUATIONS OF MOTION FOR ACTIVE MAGNETIC BEARINGS AND FLYWHEEL WITH MASS IMBALANCE.............................................................................. 3.1 ? General Structure.............................................................................. 3.2 ? Current, Position, and Force Correlation.......................................... 3.3 ? Linearization Using Taylor Series Expansion.................................. 3.4 ? Equations of Motion Including Rotor-Dynamics............................. 3.4.1 ? Generic Flywheel and Bearing Model Introduction and Assumptions...................................................................... 3.4.2 ? Derivation of the Homogeneous System............................................................................... 3.4.3 ? Inclusion of Mass Imbalance and Generalized Forces................................................................................ 3.4.4 ? Finalizing the Equations of Motion for the Active Magnetic Bearing Rig....................................................... CHAPTER 4 ? CONTROLLING THE ACTIVE MAGNETIC BEARING........ xii xviii xx 1 5 12 12 15 18 20 20 22 25 27 29 x 4.1 ? Amplifier Dynamics......................................................................... 4.2 ? Linear State Variable Model............................................................ 4.3 ? State Estimation................................................................................ 4.4 ? Observer Feedback Control.............................................................. 4.5 ? Perturbation and Integral Control..................................................... 4.5.1 ? Gravity and Mass Imbalance............................................. 4.5.2 ? Integral Control................................................................. 4.6 ? Adaptive Disturbance Rejection....................................................... 4.7 ? Control Voltage Splitting................................................................. CHAPTER 5 ? NUMERICAL SIMULATIONS AND RESULTS...................... 5.1 ? Adaptive Gain, Hp............................................................................ 5.2 ? Simulation Results............................................................................ 5.2.1 ? Ideal Conditions................................................................ 5.2.2 ? Noisy Conditions............................................................... 5.2.3 ? Pole Variation.................................................................... 5.2.4 ? ?H Variation....................................................................... 5.2.5 ? Discrepancies between Running Speed and ADR Frequency.......................................................................... 5.2.6 ? Introduction of Other Frequencies.................................... 5.3 ? Response Time................................................................................. CHAPTER 6 ? EXPERIMENTAL APPARATUS AND RESULTS................... 6.1 ? Hardware Configuration................................................................... 6.2 ? Controller Software and Code.......................................................... 29 31 33 36 38 38 39 39 41 46 47 48 48 51 59 64 67 75 88 90 91 96 xi 6.3 ? Experimental Results........................................................................ 6.3.1 ? Simulated Imbalance......................................................... 6.3.2 ? ?H Variation.................................................................... 6.3.3 ? Frequency Discrepancies................................................... 6.3.4 ? Additional Imbalance by Spinning the Rotor.................... 6.4 ? AFRL Rotor Test Rig....................................................................... 6.4.1 ? Test Wheel......................................................................... 6.4.2 ? Experimental Results......................................................... 6.5 ? Response Time................................................................................. CHAPTER 7 ? CONCLUSIONS AND RECOMMENDATIONS...................... 7.1 ? Conclusions...................................................................................... 7.2 ? Recommendations............................................................................ BIBLIOGRAPHY................................................................................................. APPENDICES..................................................................................................... APPENDIX A ? COMPUTER SIMULATION CODE........................................ APPENDIX B ? EXPERIMENTAL CONTROL CODE..................................... APPENDIX C ? MAGNETIC BEARING AND TEST WHEEL DRAWINGS.. APPENDIX D ? EQUIPMENT AND MANUFACTURER INFORMATION... 100 100 107 109 112 120 121 126 131 133 133 135 136 139 140 159 181 192 xii LIST OF FIGURES Figure 1-1 Radial Active Magnetic Bearings and Rotor................................. Figure 3-1 Heteropolar Configuration............................................................ Figure 3-2 Flux Path in a Pole Pair................................................................. Figure 3-3 Flux Paths throughout a Non-Split Flux, Heteropolar Stator........ Figure 3-4 Complementary Pole Pairs and Rotor for One Axis Actuation..... Figure 3-5 Generic Flywheel and Bearing Model with Mass Imbalance....... Figure 3-6 Displacement of the Rotor............................................................. Figure 3-7 Orthogonal Components of a Mass Imbalance............................. Figure 4-1 Block Diagram of an Open Loop Linear State Variable Model.... Figure 4-2 Block Diagram of State Estimator................................................. Figure 4-3 Block Diagram of Closed Loop Plant with Observer Feedback... Figure 4-4 Block Diagram of Complete Closed Loop System....................... Figure 5-1 Time Trace for Ideal Conditions................................................... Figure 5-2 H*p for Ideal Conditions................................................................ Figure 5-3 Time Trace for Noisy Conditions (Case I).................................... Figure 5-4 H*p for Noisy Conditions (Case I)................................................. Figure 5-5 Comparison of H*p for Ideal and Noisy Conditions (Case I)........ Figure 5-6 Time Trace for Noisy Conditions (Case II)................................... 2 13 14 15 17 21 24 26 32 35 37 45 49 50 51 52 53 55 xiii Figure 5-7 H*p for Noisy Conditions (Case II)................................................ Figure 5-8 Percentage Increase in H*p vs. Signal to Noise Ratio for ? = 20 Hz...................................................................................... Figure 5-9 Time Trace for Noisy Conditions (Case III)................................. Figure 5-10 H*p for Noisy Conditions (Case III).............................................. Figure 5-11 Time Trace for System Poles = -500 + 250i................................. Figure 5-12 H*p for System Poles = -500 + 250i.............................................. Figure 5-13 Settling Time of H*p as a Function of Pole Placement.................. Figure 5-14 Amplitude of H*p as a Function of Pole Placement....................... Figure 5-15 Time Trace for ?H = 100,000........................................................ Figure 5-16 H*p for ?H = 100,000..................................................................... Figure 5-17 Settling Time of H*p as a Function of Increase in ?H.................... Figure 5-18 H*p for Discrepancy between Running Speed and ADR Frequency (Case I)........................................................................ Figure 5-19 FFT of H*p between 10 - 15 Seconds for Discrepancy in Running Speed and ADR Frequency (Case I)............................... Figure 5-20 FFT of H*p between 15 - 20 Seconds for Discrepancy in Running Speed and ADR Frequency (Case I)............................... Figure 5-21 H*p for Discrepancy between Running speed and ADR Frequency (Case II)....................................................................... Figure 5-22 FFT of H*p between 10 - 15 Seconds for Discrepancy between Running Speed and ADR Frequency (Case II)............................. Figure 5-23 FFT of H*p between 15 - 20 Seconds for Discrepancy between Running Speed and ADR Frequency (Case II)............................. Figure 5-24 Percentage Change in H*p as a Function of Discrepancy between Running Speed and ADR Frequency.............................. 56 57 58 59 60 61 62 63 64 65 66 68 69 70 71 72 73 74 xiv Figure 5-25 H*p for Two Excitation Frequencies (Case I)................................ Figure 5-26 FFT of H*p for Two Excitation Frequencies (Case I).................... Figure 5-27 for Two Excitation Frequencies (Case II).............................. Figure 5-28 FFT of H*p for Two Excitation Frequencies (Case II).................. Figure 5-29 H*p for Two Excitation Frequencies (Case III)............................. Figure 5-30 H*p for Two Excitation Frequencies (Case IV)............................. Figure 5-31 H*p Curve for Secondary ADR Frequency.................................... Figure 5-32 Time Trace for Increase in Running Speed and ADR Frequency. Figure 5-33 H*p for Increase in Running Speed and ADR Frequency.............. Figure 6-1 Active Magnetic Bearing Rig with Motor and Adjustable Driver Figure 6-2 Close Up of Active Magnetic Bearing Rig................................... Figure 6-3 Diagram of Op-Amp Circuit......................................................... Figure 6-4 Project Box for Active Magnetic Bearing Rig.............................. Figure 6-5 Diagram of Collocation Issues...................................................... Figure 6-6 Time Trace for Simulated Imbalance at 20 Hz (Case I)................ Figure 6-7 H*p for Simulated Imbalance at 20 Hz (Case I)............................. Figure 6-8 Time Trace for Simulated Imbalance at 20 Hz (Case II).............. Figure 6-9 H*p for Simulated Imbalance at 20 Hz (Case II)........................... Figure 6-10 Time Trace for ?? = 80,000.......................................................... Figure 6-11 H*p for ?? = 80,000....................................................................... Figure 6-12 H*p for Frequency Discrepancy..................................................... Figure 6-13 FFT of H*p between 12 - 15 Seconds for Frequency Discrepancy 76 77 78 79 80 82 84 86 87 91 93 94 96 98 102 103 105 106 107 108 110 111 ? pH xv Figure 6-14 FFT of H*p between 15 - 30 Seconds for Frequency Discrepancy................................................................................... Figure 6-15 FFT of Time Trace for Rotor Spinning at 20 Hz.......................... Figure 6-16 FFT of H*p for Rotor Spinning at 20 Hz....................................... Figure 6-17 H*p for Simulated Imbalance with Spinning Rotor....................... Figure 6-18 FFT of Time Trace for Simulated Imbalance and Spinning Rotor............................................................................................. Figure 6-19 FFT of H*p for Simulated Imbalance and Spinning Rotor............ Figure 6-20 H*p for Simulated Imbalance and Rotor Speed at 20 Hz............... Figure 6-21 AFRL Magnetic Bearing Test Rig................................................ Figure 6-22 Cutaway Model of Test Wheel...................................................... Figure 6-23 Photograph of Delrin Test Wheel.................................................. Figure 6-24 Initial Wishbone Bracket Actuator Design................................... Figure 6-25 AFRL Magnetic Bearing Test Rig with Test Wheel and Wishbone Bracket Actuator.......................................................... Figure 6-26 Time Trace for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz............................................................................ Figure 6-27 H*p for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz........................................................................................... Figure 6-28 Time Trace of Rotor Speed for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz..................................................... Figure 6-29 FFT of H*p between 0-8 Seconds for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz............................................. Figure A-1 Simulation Model Block Diagram................................................ Figure A-2 Simulation Controller Block Diagram.......................................... Figure A-3 Simulation Estimator Block Diagram........................................... 112 113 114 115 116 117 118 121 122 123 124 125 127 128 129 130 141 142 143 xvi Figure A-4 Simulation ADR Block Diagram.................................................. Figure A-5 Simulation Gp Subsystem Block Diagram.................................... Figure A-6 Simulation Disturbance Vector Block Diagram............................ Figure A-7 Simulation Hp Subsystem Block Diagram.................................... Figure A-8 Simulation H*p Generator Block Diagram................ Figure A-9 Simulation Imbalance and Gravity Generator Block Diagram..... Figure A-10 Simulation Integral Controller Block Diagram............................. Figure A-11 Simulation Control Voltage Splitter Block Diagram.................... Figure B-1 Experimental Model Block Diagram............................................ Figure B-2 Experimental A/D Converter Block Diagram............................... Figure B-3 Experimental Position Adjustment Block Diagram...................... Figure B-4 Experimental Controller Block Diagram...................................... Figure B-5 Experimental Estimator Block Diagram....................................... Figure B-6 Experimental ADR Block Diagram.............................................. Figure B-7 Experimental Disturbance Vector Block Diagram........................ Figure B-8 Experimental Gp Subsystem Block Diagram................................ Figure B-9 Experimental Hp Subsystem Block Diagram................................ Figure B-10 Experimental H*p Generator Block Diagram................................ Figure B-11 Experimental Integral Controller Block Diagram......................... Figure B-12 Experimental Imbalance Generator Block Diagram..................... Figure B-13 Experimental Ramp Generator for Translational Imbalance Block Diagram.............................................................................. Figure B-14 Experimental Ramp Generator for Rotational Imbalance Block Diagram......................................................................................... 144 145 146 147 148 149 150 151 160 161 162 163 164 165 166 167 168 169 170 171 172 173 xvii Figure B-15 Experimental Control Voltage Splitter Block Diagram................ Figure B-16 Experimental D/A Converter Block Diagram............................... Figure C-1 Drawing of Magnetic Bearing Base Housing.............................. Figure C-2 Drawing of Magnetic Bearing Probe Housing............................. Figure C-3 Drawing of Isometric Views for Magnetic Bearing Probe Housing......................................................................................... Figure C-4 Drawing of Magnetic Bearing Stator Housing.............................. Figure C-5 Drawing of Isometric View for Magnetic Bearing Stator Housing......................................................................................... Figure C-6 Drawing of Magnetic Bearing Stator Housing Face Plate............ Figure C-7 Drawing of Magnetic Bearing Stator............................................ Figure C-8 Drawing of Magnetic Bearing Assembly Order............................ Figure C-9 Drawing of Magnetic Bearing Rotor............................................. Figure C-10 Drawing of Test Wheel................................................................. 174 175 182 183 184 185 186 187 188 189 190 191 xviii LIST OF TABLES Table 3-1 Description of Terms Used for Generic Flywheel and Bearings.. Table 5-1 Parameters for Ideal Conditions.................................................... Table 5-2 Parameters for Noisy Conditions (Case I)..................................... Table 5-3 Parameters for Noisy Conditions (Case II)................................... Table 5-4 Parameters for Noisy Conditions (Case III).................................. Table 5-5 Parameters for System Poles = -500 + 250i.................................. Table 5-6 Parameters for ?H = 100,000......................................................... Table 5-7 Parameters for Discrepancy between Running Speed and ADR Frequency (Case I)........................................................................ Table 5-8 Parameters for Discrepancy between Running Speed and ADR Frequency (Case II)....................................................................... Table 5-9 Parameters for Two Excitation Frequencies (Case I).................... Table 5-10 Parameters for Two Excitation Frequencies (Case II).................. Table 5-11 Parameters for Two Excitation Frequencies (Case III)................. Table 5-12 Parameters for Two Excitation Frequencies (Case IV)................. Table 5-13 Parameters for Secondary ADR Frequency.................................. Table 5-14 Parameters for Increase in Running Speed and ADR Frequency. Table 6-1 Magnetic Bearing and Rotor Properties........................................ Table 6-2 Parameters for Simulated Imbalance at 20 Hz (Case I)................ 22 49 52 55 58 60 65 68 72 76 79 81 82 84 87 92 103 xix Table 6-3 Parameters for Simulated Imbalance at 20 Hz (Case II)............... Table 6-4 Parameters for ?H = 80,000........................................................... Table 6-5 Parameters for Frequency Discrepancy......................................... Table 6-6 Parameters for Simulated Imbalance and Spinning Rotor............ Table 6-7 Parameters for Simulated Imbalance and Rotor Speed at 20 Hz.. Table 6-8 Parameters for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz............................................................................ Table 6-9 Time Responses of Experimental Plots.................................. 105 108 110 116 119 128 132 ? pH xx NOMENCLATURE A State matrix Ag Area of pole at air gap (m2) B Input matrix C Output matrix D Input matrix d Distance of mass imbalance from center of flywheel (meters ? m) e error vector F Force (Newtons ? N) Fn Net force (N) Gp Adaptive gain g Nominal air gap (m) Hp Adaptive gain ? pH Adaptive imbalance detection gain IP Polar moment of inertia (kilogram-meter2 ? kg-m2) IT Transverse moment of inertia (kg-m2) i Unit vector associated with body fixed X axis, Electrical current (Amperes ? A) ib Bias current (A) ic Control current (A) xxi J General mass moment of inertia (kg-m2) j Unit vector associated with body fixed Y axis K Feedback gain matrix, spring stiffness (N/m) Ka Amplifier gain (A/V) Kg Integral gain Ki Current stiffness (N/A) Ks Sensor gain (V/m) Kp Position stiffness (N/m) k Unit vector associated with body fixed Z axis L Observer gain matrix L1 Length from center of flywheel to left bearing (m) L2 Length from center of flywheel to right bearing (m) M Mass of flywheel and rotor (kg) Mo Moment (N-m) m Mass of imbalance (kg) N Number of coil turns n Number of states or variable Q Generalized force (N) q Degree of freedom T Kinetic energy (Joules ? J), transfer matrix t time (seconds ? s) ud Disturbance vector up Input vector xxii V Potential energy (J) v Voltage (volts ? V) X Orthogonal axis, translation (m) x State vector, rotor perturbation Y Orthogonal axis, translation (m) y Output vector Z Geometric axis of rotation Z' Principal axis of inertia ? Rotation about Y axis (radians ? rad) ? Rotation about X axis (rad) ? Disturbance matrix ? Weighting matrix ?f Frequency Discrepancy (Hz) ? Angle between geometric axis of rotation and principal axis of inertia (rad) ?o Permeability of free space ? Amplifier time constant (s) ? Disturbance function ?f Phase (rad) ? Overall angular velocity (rad/s) ? Frequency, running speed (hertz ? Hz, rad/sec) xxiii Subscripts b Bias c Control d Disturbance f Frequency G Adaptive gain Gp H Adaptive gain Hp i Current, arbitrary number p Position real Real plant Sin Sine component Cos Cosine component X X axis Y Y axis 1 Left bearing 2 Right bearing I Primary pole pair II Complementary pole pair xxiv Other Nomenclature ^ (circumflex) Estimated state, unit vector ? (dot) Time derivative _ (underscore) Vector 1 CHAPTER I INTRODUCTION Since the invention of the wheel as well as other rotating machinery, the world has seen the need for support bearings. One common factor among most bearings is the fact that they require mechanical contact between a rotating surface and a non-rotating surface. This contact induces energy loss, mechanical wear, and eventually failure of the bearing and sometimes even the rotor that it supports. Most modern applications of bearings address this issue with preventative methods. However, as technology has progressed over the decades, other means of supporting rotors have become available. Some of the more modern bearings incorporate elaborate systems using either air or a hydrodynamic fluid to reduce friction. In addition, by harnessing the power of electromagnetism, it has become possible to levitate a rotor, thus alleviating any energy loss and wear due to mechanical contact as well as the need for a complex lubrication system. Although the idea of magnetic bearings and magnetic levitation has been around since the mid-1800s, it wasn't until the advent of modern sensing equipment that practical designs have been developed. Magnetic bearings operate using electromagnets that generate only attractive forces and are therefore inherently unstable. Thus for a horizontally positioned rotor to maintain stability, two radial bearings with electromagnetic coils placed around the rotor are 2 required, as well as a thrust bearing acting in the axial direction, resulting in a total of five axes of support. The electromagnetic coils inside each bearing must be actively controlled using a digital processor which can take sensor input and generate the appropriate output to the coils. Figure 1-1 below shows a three-dimensional drawing of two radial active magnetic bearings with the rotor pulled out for illustrative purposes. Figure 1-1 ? Radial Active Magnetic Bearings and Rotor Active magnetic bearings are well suited for many applications in industry. Their ability to operate without contact does away with the need for lubrication systems making them particularly suitable for operation in a vacuum, at extreme (high or low) temperatures, and in corrosive fluids. This unique advantage also eliminates the high cost associated with complex cooling systems, regularly scheduled maintenance and spare parts. Since magnetic bearings require active control, inherent rotational vibration 3 as well as other perturbations can be compensated for, resulting in low energy loss and a much more energy efficient machine. Possible applications include gas turbine engines in aircraft and turbo-pumps in rocket engines. Although the use of magnetic bearings can certainly reduce the cost of operations over time, the initial investment of magnetic bearings has traditionally been high and is one of the main factors that have kept them from wider use and application. Active magnetic bearings have received much recent attention for their potential use in satellites and spacecraft. Specifically designed flywheels used in conjunction with magnetic bearings have the potential to not only surpass chemical batteries in their ability to store mechanical energy, but they can also be used for attitude control of the satellite itself. Moreover, this type of system could reduce the overall weight of the satellite resulting in a major reduction in cost of launch. Conceptually, these flywheels will be constructed using high strength composite materials, such as multi-direction composites (MDC) or filament-wound composites, resulting in high inertial properties allowing for a high energy storage density to weight ratio. These disks will typically spin at very high speeds with a possible operating range from 10kRPM to 100kRPM as energy is added via a motor and removed via a generator. As explained by Hai et al. [10 and 11], variable spin rates will produce significant cyclical stresses in the disk, resulting in possible delamination of the composite material, debonding of the hub and rim, and/or fatigue induced cracks. If any failures in the flywheel are not promptly detected, the integrity of the disk can be compromised and the entire structure face catastrophic failure. Fortunately, such flaws produce changes in the balance state of the flywheel. By monitoring certain adaptive control gains during operation, these types of failures can be 4 identified in their infancy. Appropriate methods can then be taken to contain the system in as safe a way as possible. Through a proposed vibration-suppression method called Adaptive Disturbance Rejection (ADR), compensation for sudden imbalances and perturbations introduced to the system begins almost instantaneously. This allows for the characteristics of crack growth to be possibly detected within a single revolution of the rotor. More importantly, this method can be used to detect slight changes in imbalance without referencing a time trace of the system and is robust with regard to signal noise. Some of the current methods of monitoring a flywheel rotor supported by magnetic bearings use algorithms in conjunction with a time trace of the rotor vibration. Should the rotor position exceed some boundary as defined by the algorithm, the system is shut down. This particular approach provides no insight to the actual cause of the boundary crossing and does not necessarily monitor the "health" of the flywheel itself. It is the purpose of this thesis to evaluate through both computer simulation and two experimental rigs the effectiveness of the ADR health monitoring approach. A derivation of the equations of motion, including rotor dynamics, will be presented as well as a state space controller with appended integral and adaptive controls. Next, the development and integration of a crack simulation test wheel will be discussed. Finally, a detailed discussion of the results of the simulations and experiments, and suggestions for future work will be presented. 5 CHAPTER II LITERATURE REVIEW Significant research and improvements have been made to the concept of levitation and magnetic bearings since the ideas were first conceived. Samuel Earnshaw [5] was the first to publish work about such ideas back in 1842. However, the realization of such a technically complex, yet simplistic and novel approach to solving the problems of mechanical contact in rotor systems has only been around since the 1960's. Since the magnetic bearing's inception, it has lent itself to new and innovative ideas in the area of controls engineering. With each new breakthrough, applications using magnetic bearings are being dreamed up and implemented that were beforehand not thought possible. Current major industrial uses of magnetic bearings consist of compressors, power turbines, overhung compressors, turbo expanders and turbo pumps. The most common applications are being used in centrifugal compressors and turbo expanders [27]. Recently however, magnetic bearings have been applied to other industries. Because of their ability to operate in sealed environments without lubrication, in addition to the ability to function submerged without regard for temperature restrictions or corrosive fluids, magnetic bearings are now being considered and used in the beverage and food industries as well as the semi-conductor equipment industry [20]. Even more fascinating are their use in biological and pharmaceutical applications involving life cell processing 6 [20]. In a paper by Shen et al. [24], it was reported that a new breed of artificial hearts using a permanent magnet synchronous motor rotor and pump impeller configuration levitated by magnetic bearings was being investigated. All of this new technology involving magnetic bearings is the result of years of dynamics and controls research. In 1981, Matsumura et al. [18] derived the equations of motion for the magnetic bearing and rotor including rotor dynamics and succeeded in showing how to get rid of steady state error using an integral type controller. Flowers et al. [7] utilized this integral control approach by augmenting and appending it to a state feedback controller. This allowed for a state space approach to controlling the system while still evading the traditional steady state error typical in state feedback systems. Wilson [29] developed and verified zero bias and low bias control law methods that globally asymptotically stabilize active magnetic bearing systems. He went on to construct the largest domain of definition possible while minimizing energy losses by reducing the total square flux required for regulation through the control laws. Although it is observed that zero bias laws require larger voltage inputs to the system, he notes that amplifier bandwidth is more of a limiting factor on actuation than voltage saturation. Knopse et al. [11] suggested the use of an adaptive open loop controller using recursive gain scheduling. This approach was proven to cancel out vibration as the rotor was spun up through the critical speed. Typically, the rotor will see the highest amplitudes of vibration as it approaches the speed congruent with the rotor's first natural frequency, also known as the critical speed. Mohamed et al. [19] developed a method of canceling out these high vibratory amplitudes using a nonlinear approach. He found that the gyroscopic motion of the rotor effected stability at high speeds and caused Hopf 7 bifurcation to periodic solutions as speeds varied through the critical speed. By deriving a nonlinear feedback controller, he was able to control the Hopf bifurcation tendency through the first natural frequency. This brings up a fundamental difference in control strategies. While many researchers tend to linearize the plant model for simplicity, magnetic bearings are in fact very nonlinear due to the electromagnetic physics involved. Thus, the linear versus nonlinear approaches to control algorithms have both their advantages and disadvantages, depending on the desired application. Although a rotor used in conjunction with magnetic bearings can be stabilized using many different control methods, vibration due to imbalance is inherent in all rotating machinery. There have been many methods developed to compensate and compress this natural tendency to oscillate. Ahmed et al. [1] discovered that by using a technique called Adaptive Force Balancing, a planar rotor could be asymptotically stabilized without any knowledge of the position of the center of mass. Shafai et al. [23] also used this technique along with LQG/LTR, H? and QFT algorithms to solve the imbalance problem. Although the Adaptive Force Balancing method presented by Ahmed et al. mainly addressed a single input, single output (SISO) system, a very brief treatment of a multi-input, multi-output (MIMO) system was given as well. However, another method of rejecting imbalance and general disturbances was presented by Fuentes et al. [8]. The paper developed a rigorous mathematical proof which set up a type of control law called model reference. By including a model reference which contains no disturbances and a control law that takes into account expected disturbances, the controller will track the output of the plant and force it to behave like the reference. 8 One of the advantages of this type of controller is that the disturbances to be rejected only have to be bounded, continuous functions. This type of controller lends itself very well to MIMO systems such as a horizontal magnetic bearing setup. It could however, be applied to any kind of system. Although the work done by Fuentes et al. [8] only contained numerical simulations, Matras et al. [15 and 16] implemented the controller on an active magnetic bearing rig showing that this new technique, entitled Adaptive Disturbance Rejection (ADR), could in fact suppress persistent excitations at known frequencies within a certain range. Work done later by Matras [17] showed heuristically that modifications to the control law allowed for the suppression of higher frequencies. By allowing one adaptive gain to act on the estimator output while the other adaptive gain remained connected to the sensor readings, he was able to achieve near infinite gain margins allowing for frequencies extending through and beyond the critical speed to be rejected. This thesis incorporates the different variations of the technique known as ADR and is discussed in general in later sections. Research involving magnetic bearings has gone beyond the field of controls and is now being focused by some at modifying the components of the magnetic bearings themselves. Maslen [14] has a compendium giving insights to different approaches to the structure of magnetic bearings. A new breed of magnetic bearings called superconducting magnetic bearings (SMB) is currently being investigated. Coombs et al. [4] describes the dynamics of SMB using the compound YBa2Cu3O7 as an active component and showed that although they have a high load capacity, they also have very little inherent damping and low stiffness properties. It is also stated that cyclical loads at or near the natural frequency can have catastrophic effects. Ichihara et al. [12] developed 9 an experimental radial SMB system using YBCO superconductors and an outer rotor of permanent magnets. This system was addressed as a 10kWh-class energy storage system and was successfully implemented storing 2.24kWh while operating at 7,500 RPM. Both Ichihara et al. [12] and Wilson [29] have produced results using magnetic bearings for a particular type of combined application that is being approached very seriously by the Air Force and NASA. This application consists of using composite flywheels supported by magnetic bearings for a combination of energy storage and attitude control. NASA [21] has already published numerical simulations pertaining to this application being used on the International Space Station called the Attitude Control and Energy Storage Experiment (ACESE). The setup employed a configuration of two flywheels, each spinning in opposite directions. The primary objective of the experiment was to demonstrate the ability to store the same amount of energy as two on-board batteries, whilst the second objective was to exert torque on the space station when storing or dissipating energy. Contingency plans such as the loss of one rotor were also simulated and the experiment overall gave good results. Fausz et al. [6] announced plans with the Air Force Research Laboratory to begin the FACETS (Flywheel Attitude Control & Energy Transmission and Storage) Grand Challenge. This work is currently in progress and combines energy storage, attitude control, and Power Management and Distribution (PMAD) duties for the flywheel system to achieve. By using four or more wheels in a reaction mode and utilizing the null subspace of the angular momentum dynamics of the wheels, simultaneous momentum management and power tracking can be accomplished efficiently. This work is being done largely as an attempt to combine 10 these primary systems into a single device to be used on satellites thus making them more efficient. The work in this thesis is part of the ongoing research for this project. In order for energy storage and attitude control to be implemented using flywheels, the flywheels must be spun up to very high speeds. The work done by NASA [21] reported a range of 15-50kRPM, while the operating range of the FACETS Grand Challenge requires 30-100kRPM. When operating at speeds in a range this high, the flywheels in question (which are typically constructed using composites) will see high stress loads. Thus, health monitoring and an understanding of crack dynamics in composite materials is an integral part of a working energy storage system. Recent research in the area of crack initiation, propagation, and composite failure shows promise. Zhu et al. [30] stated that the effect of cracks on rotors and flywheels should be considered in designing active magnetic bearing controllers. He observed that monitoring the super-harmonic components, specifically the 2X and 3X revolutions, in sub-critical speed regions could be used as an index to detect cracks. He also states that it is impossible to use these super-harmonics when dealing with super-critical speeds. Moreover, Amati et al. [2] showed that vibration monitoring of rotors on active magnetic bearings is significantly affected by the unbalanced magnetic pull applied by induction motors throughout sub-critical speeds prompting the need for a different strategy of condition monitoring. Hai [10] explains that the maximum stresses in composite rotors occur in the hub area. In a subsequent paper, Hai et al. [11] observed that cracks reduce the effective natural frequency and this effect is magnified as the speed increases. She also represents that the centrifugal effect of operating speeds acts as a negative stiffness and thus reduces 11 the effective stiffness of the overall system until "catastrophic failure occurs." Shiue et al. [25] developed a speed control approach for virtual containment, maintaining some level of substantial energy storage regardless of flywheel failure. In another paper by Shiue et al. [26], it is stated that insight into the size and severity of a flywheel flaw can be provided by monitoring imbalance. In the same paper, an experimental rig was set up where a mass was added to a flywheel with adhesive tape. The rotor was spun up until the mass was ejected thus causing an effective change in the balance state and simulating a crack initiation. In this thesis, an active magnetic bearing will be modeled and controlled using feedback and adaptive disturbance techniques. Via an experimental crack initiation simulation, results will be found such that imbalance changes can be detected on-line. These crack simulations will be done through a couple different methods. One method incorporates a specifically designed test wheel consisting of a movable mass that can be triggered by magnets through a wide range of speed. The other method uses the magnetic bearing controller to perturb the system using sinusoidal functions similar to that of an imbalance. By using these techniques with the aforementioned ADR control strategy implemented on active magnetic bearings, it will be shown that the characteristics of flywheel failure and imbalance changes can be identified by observing the automatic on- line adjustments of the adaptive gains. 12 CHAPTER III EQUATIONS OF MOTION FOR ACTIVE MAGNETIC BEARINGS AND FLYWHEEL WITH MASS IMBALANCE A horizontally positioned active magnetic bearing rig usually consists of two radial bearings, each one controlling the two translational modes of the rotor within the plane of the bearing. Consequently, the two non-axial rotational modes are controlled as well. Typically a thrust bearing is used to constrain any translational movement in the axial direction. This allows for the rotor to be completely controlled and centered within the bearings for a total of five degrees of freedom leaving only the ability to rotate axially remaining. Obviously, this is accomplished using electromagnets contained within each bearing. 3.1 ? General Structure Magnetic bearings are generally assembled using groups of electromagnetic pole pairs. Since the forces generated by the pole pairs are attractive in nature, the pole pairs require complementary pole pairs opposite of the rotor. The simplest and most common scheme uses four pole pairs such that the rotor can be actuated in both planar orthogonal directions. One general layout of these pole pairs for a radial actuator design is known as a heteropolar configuration. This type of configuration consists of the poles being placed 13 circumferentially around the rotor. A diagram depicting this arrangement can be seen below in figure 3-1. Figure 3-1 ? Heteropolar Configuration A pole pair, usually in the shape of a horseshoe, is constructed using stacked laminates of iron forming a core, and a conductive wire is coiled in opposite directions around each leg to form two poles of opposite polarity. The pole pairs extend out very close to the rotor such that only a very small air gap, typically a few thousandths of an inch, separates the poles from the rotor. As current is passed through the coils, a magnetic flux is generated. This flux, much like current in an electrical circuit, must follow a path. Beginning in the coil of one leg or pole, the flux passes through an air gap into the rotor. The part of the rotor that reacts with the magnetic flux, usually fabricated from electrochemically cut laminates stacked over a solid shaft, propagates the flux through to the second air gap and into the other leg of the pole pair of opposing polarity. In diagram 3-2, the flux path in a pole pair and rotor can be seen. As the current is 14 increased and decreased, the magnetic flux fluctuates thus varying the attractive forces applied. Figure 3-3 depicts the flux paths throughout a heteropolar stator with a non-split flux design. Figure 3-2 ? Flux Path in a Pole Pair Figures 3-1 through 3-3 represent only one type of radial actuator for a magnetic bearing. Although the heteropolar, non-split flux configuration is the simplest and most common, other designs and arrangements of pole pairs, flux paths, and stator types exist. Maslen [14] covers a variety of these in detail. However, the experimental rigs used for this research are of this type, so the focus will remain on this particular kind. Coils Flux Path Rotor Laminates Solid Rotor Core Pole Pair Heteropolar Stator Air Gap 15 Figure 3-3 ? Flux Paths throughout a Non-Split Flux, Heteropolar Stator 3.2 ? Current, Position, and Force Correlation The magnetic forces produced by the pole pairs of the stator manifest themselves as an attraction of the rotor face to the surface of each pole end. These attractive forces generated by the stator are dependent on the amount of electrical current provided to the coils and the distance or air gap between the face of the rotor and the surfaces of the pole pair. When the rotor is away from center position, a difference in the distance from the rotor face to each leg of the pole pair can be observed. Although this does tend to couple the two orthogonal directions to be controlled, the air gap is typically small enough that the translational movement of the rotor can be decoupled into the two separate axes of actuation thus simplifying the modeling procedure. Also, it should be noted that each 16 pole of the pole pair is at an angle to the general axis of actuation. However, since the pole pair is symmetric about this axis, the side forces produced by the poles cancel each other out leaving only the resultant force in the direction of the axis. Although the current (i) supplied to the coils and the air gap (g) between the rotor and the poles are considered to be variable inputs, other factors based on stator design also apply in determining the force provided by each pole. These other factors consist of the number of windings on each pole (N), the permeability of free space (?o), and the cross-sectional area at the end of the pole (Ag). Both coils of the pole pair are connected via the conductive wire such that the current passing through each coil is the same. Moreover, taking into account that there are two poles per pole pair, thus two coils and two air gaps, the general structure of the force equation below allows for the doubling of the windings and air gap to cancel each other out, hence producing equation 3-1. (3-1) Once again, the attractive nature of the pole pair and rotor requires that a second pole pair be present on the opposite side of the rotor such that the rotor can be actuated in either direction of the axis as needed. Hence, the net force provided along the axis of actuation is extended to include both pole pairs and can be seen in equation 3-2 below. (3-2) The subscripts I and II in the equation above denote the primary and complementary pole pairs respectively and a diagram detailing the geometry and terms of the equation can be found in figure 3-4. ??? ? ??? ? ? ?= 22 22 2 0 III III gn gg iiANF ? 2 2 2 0 g iANF g?= 17 Figure 3-4 ? Complementary Pole Pairs and Rotor for One Axis Actuation As seen in the figure above, x is the vertical axial perturbation from center position. Hence, if the term g is the uniform air gap, or bias position when the rotor is centered between the pole pairs, then gI and gII can be represented by equations 3-3 and 3-4 below. (3-3) (3-4) By using the same logic as above, the coil currents can be treated the same way. Such that if a bias current (ib) is applied to both pole pairs at all times and a control Y xggI ?= xggII += iI iII gI gII X x Ag ?? N 18 current (ic) is applied as needed, then the currents expressed earlier as iI and iII can now be expressed in a similar way to the position equations. It should be noted that when the rotor is being actuated upward in a positive x direction, the top current is increased and the bottom current is decreased, such that the currents iI and iII can be now be expressed as given in equations 3-5 and 3-6. (3-5) (3-6) The subscripts remain on the bias currents so that if necessary, separate bias currents can be assigned to the opposing pole pairs. One advantage to doing this is the ability to factor in and compensate for constant forces, such as gravity. By substituting equations 3-3 through 3-6 into equation 3-2, the magnetic force can now be represented as equation 3-7 below. (3-7) 3.3 ? Linearization Using Taylor Series Expansion In order to simplify the derivation and set up the model such that classical control methods can be used, equation 3-7 can be linearized with respect to x and ic using the first terms of a Taylor series expansion. By linearizing about x = 0 and ic = 0, the net force produced by the two pole pairs can be approximated by equation 3-8 below. (3-8) cbI iii I += cbII iii II ?= ( ) ( ) ( ) ( ) ??? ? ? ? ? ? + ?? ? += 2 2 2 2 2 0 xg ii xg iiANF cbcb gn III? c c nn n ii Fx x FF ? ?+ ? ?? 19 The partial derivatives of Fn from equation 3-7 can be thought of as stiffness terms as they act much like a spring stiffness. Thus, the partial derivative of the net force with respect to the position x will be called the position stiffness (Kp) and the partial derivative of the net force with respect to the control current ic will be called the current stiffness (Ki). These stiffness terms, evaluated at x = 0 and ic = 0, are respectively derived below, and end up as equations 3-9 and 3-10. (3-9) As before, the subscripts referring to the opposing pole pairs remain intact as the bias currents may differ between them. (3-10) Now that the magnetic force equation has been linearized about the inputs x and ic, it can be represented as equation 3-11 below. (3-11) ( ) ( ) ( ) ( ) 0,0 2 2 2 2 2 0 == ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ?? ? + ? ?== ? ? c III ix cbcb gp n xg ii xg iiAN xKx F ? ( ) ( ) ( ) ( ) 0,0 3 2 3 2 2 02 == ? ? ? ? ? ? ? ? + ?+ ? += c III ix cbcb gp xg ii xg iiANK ? ??? ? ??? ? += 3 22 2 02 g iiANK III bb gp ? ( ) ( ) ( ) ( ) 0,0 2 2 2 2 2 0 == ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? + ?? ? + ? ?== ? ? c III ix cbcb g c i c n xg ii xg iiAN iKi F ? ( ) ( ) ( ) ( ) 0,022 2 02 == ??? ? ??? ? + ?+ ? += c III ix cbcb gi xg ii xg iiANK ? ??? ? ??? ?= ??? ? ??? ?= 2 2 02 2 0 22 g iANK g iANK II II I I b gi b gi ?? IIIIII cicipn iKiKxKF ++= 20 3.4 ? Equations of Motion Including Rotor-Dynamics This section will build on some basic and advanced methods of dynamic analysis to derive the equations of motion for an active magnetic bearing rig supporting a flywheel with a mass imbalance. First, a generic flywheel and bearing model will be introduced and assumptions about the system will be stated. Next, a method of derivation using LaGrange's equations will be presented to form a set of equations of motion for the generic flywheel model. From this point, mass imbalance and generalized forces will be discussed and included into the equation set. Finally, the stiffness terms for the magnetic forces derived in the above section will be substituted in for the generic stiffness terms and generalized forces leading lastly to a completed list of the modeled equations for the active magnetic bearing rig. 3.4.1 ? Generic Flywheel and Bearing Model Introduction and Assumptions A generic homogeneous flywheel and rotor supported by two bearings is represented in figure 3-5. The flywheel includes a mass imbalance and a rotor passing through the center. Each bearing will be characterized at this point by a spring and two generalized forces for both orthogonal directions, X and Y. Also, the two bearings will be denoted by subscripts 1 and 2 such that all forces and values associated with the left bearings will have the subscript 1, and all forces and values associated with the right bearings will have the subscript 2. All six degrees of freedom will be represented in the model and will be labeled as X, Y, Z, ?, ?, and ?. Table 3-1 gives a list of the terms used and their relationship to this model. 21 Z X ? ? ? K1Y, F1Y-I K1X, F1X-I K2Y, F2Y-I K2X, F2X-I Y L1 L2 d (M, IT, IP) m F1Y-II F1X-II F2Y-II F2X-II Figure 3-5 ? Generic Flywheel and Bearing Model with Mass Imbalance Many assumptions are made regarding the generic model and are presented in greater detail by Vance [28]. The term ? is a constant running speed and is the time derivative of the angular degree of freedom ?. The angles ? and ? are considered to be small, such that the trigonometric identities of ? are defined by Cos(?) = 1, Sin(?) = ?, Cos2(?) = 1, and Sin2(?) = ?2. The same is true of ?. Also, ? and ? can be regarded as rotations about the space-fixed Y and X axes respectively and therefore cease to be proper Euler angles such that ? about Z-axis is the only rotation that sees coupling forces from ? and ?. For the derivation using LaGrange's equations, only the first and second order terms are retained while third order and higher are discarded. This allows for all important dynamics to be kept, yet allowing the final model to remain linear. Finally, as 22 stated earlier, the "small air gap" assumption is made such that the forces in the X and Y directions are decoupled from one another. Table 3-1 ? Description of Terms Used for Generic Flywheel and Bearings 3.4.2 ? Derivation of the Homogeneous System LaGrange's equations are a powerful tool used by many to perform dynamic analyses of systems and will be covered briefly in this section. First, the overall angular velocity (?) of the system must be taken into account. Equation 3-12 shows how this is set up. (3-12) Symbol Description Symbol Description X,Y,Z Cartesian Coordinates with origin in center of flywheel d Distance of imbalance mass from center of flywheel ?,? Rotations about the Y and X axes respectively L1 Distance from center of flywheel to left bearing ? Constant angular speed about the Z axis L2 Distance from center of flywheel to right bearing M Mass of flywheel and rotor K1X,K2X Stiffness values associated with spring forces along X axis for left and right bearings respectively IT Transverse moment of inertia of the flywheel about X and Y axes K1Y,K2Y Stiffness values associated with spring forces along Y axis for left and right bearings respectively IP Polar mass moment of inertial of the flywheel about Z axis F1X,F2X Generalized forces along X axis for left and right bearings respectively m Mass of imbalance F1Y,F2Y Generalized forces along Y axis for left and right bearings respectively ????? ++=? kJi ??? ' 23 The term ?J refers to the space fixed axis Y, the term ?'i refers to an intermediate rotational axis x', and the terms ?i , ?j , and ?k refer to a set of axes that move with the rotor but do not spin about the rotation axis. After considering the angular rotations of ? and ? about the Y and X axes respectively, the overall angular speed can now be expressed in terms of body fixed coordinates with respect to the flywheel and rotor as seen below in equation 3-13. (3-13) Next, LaGrange's equation will be defined in equation (3-14). (3-14) In the equation above, T represents the equation for the kinetic energy of the system, V signifies the potential energy of the system, Q stands for any generalized forces acting on the system, and qi represents the various degrees of freedom that the energies and forces act upon such that q1 = X, q2 = Y, q3 = ?, and q4 = ?. Below in equation 3-15, the formulation for the kinetic energy term of the system is given where J is the general mass moment of inertia of the flywheel and rotor and M is as defined in table 3-1. (3-15) For the term (???) in the equation above, equation 3-13 has a dot product taken with respect to itself and the small angle assumptions are applied. Afterward, all third order and higher terms are discarded resulting in the equation shown below. ?????? ?????? ?+??????+??????=? kSinjCosi )()( ?????? i iii QqVqT q T dt d = ? ?+ ? ?? ?? ? ? ? ?? ? ? ? ? ? ? ( )???+?? ? ? ??? ? += ?? JYXMT 21 22 21 24 x + L1? X Z ? L1 x L2 x - L2? (3-16) At this point, equation 3-16 can be substituted into equation 3-15 and the term J is split between its transverse and polar parts as seen in equation 3-17. (3-17) Next, consider the potential energy relation. While mathematically modeling the energy stored by the springs as shown in figure 3-5, translation as well as rotation of the rotor must be taken into account. In figure 3-6, an exaggerated depiction of the effective displacement of each end of the rotor with respect to the X axis is represented. The same applies to the displacement of the rotor along the Y axis. Equation 3-18 shows the cumulative potential energy stored by the system. (3-18) Figure 3-6 ? Displacement of the Rotor ( ) ?????? ?+?? ? ? ??? ?+ ??? ? ??? ?=??? ??? ?????? 2222 ?????? ?+?? ? ? ??? ? ++ ??? ? ??? ? += ????? ?????? 22 21 22 21 22 21 PT IIYXMT ( ) ( ) ( ) ( )22221211212222121121 ???? LXKLXKLYKLYKV XXYY ?+++++?= 25 Now that the kinetic and potential energy equations have been developed, the left hand side of equation 3-14 can be completed for all four degrees of freedom. The different terms of equation 3-14 with respect to q1 are calculated in equations 3-19. Equation 3-20 shows the addition of the terms in equations 3-19 via equation 3-14 giving a nearly completed equation of motion for the degree of freedom X. The same procedure can now be done for Y, ?, and ?. (3-19a) (3-19b) (3-19c) (3-20) 3.4.3 ? Inclusion of Mass Imbalance and Generalized Forces The mass imbalance will be modeled here as a point mass (m) at a defined distance (d) from the center of the flywheel. The force applied by the imbalance is dependent on not only the mass and distance from center, but also the rotational speed of the flywheel as seen below in equation 3-21. (3-21) ?? ? =? ? ? ? ? ? ? ? ? ? XM X T dt d 0=??XT ( ) ( )?? 2211 LXKLXKXV XX ?++=?? ( ) ( ) 12211 QLXKLXKXM XX =?+++?? ?? dmFimbalance 2?= 26 Y X m ? d d Cos(? t) d Sin(? t) The force, as applied by the imbalance must be separated into its respective X and Y components. By doing this, the contributions from this force can be easily included into their respective equations of motion. Figure 3-7 depicts the mass imbalance in the plane of the flywheel and its respective components. Figure 3-7 ? Orthogonal Components of a Mass Imbalance If the imbalance mass is in an X-Y plane that is not axially centered on the rotor, it will cause a moment to be applied to the rotor as it spins. When this happens, the principal axis of inertia becomes misaligned. Therefore, if the term ? represents the constant angle between the principal axis of inertia (Z') and the geometric axis of rotation (Z), the moment affecting the rotor can be described as shown in equation 3-22. (3-22) ( ) ?? 2 PTimbalance IIMo ?= 27 Much like the imbalance force above, the moment applied to the rotor must be separated into its angular components such that it can be included into the equations of motion for ? and ?. Now that the imbalance forces and moments have been described, the inclusion of these along with other forces and moments can be included into the overall generalized force terms Qi to complete the derivation of the equations of motion for the generic flywheel and bearing system. These generalized terms will be expressed as Q1 = ?FX, Q2 = ?FY, Q3 = ?MoY, and Q4 = ?MoX such that the generalized force with respect to X can be represented as seen below in equation 3-23. The completed equation of motion for q1 can be determined by substituting equation 3-23 into equation 3-20 as seen below in equation 3-24. (3-23) (3-24) 3.4.4 ? Finalizing the Equations of Motion for the Active Magnetic Bearing Rig The equations of motion as presented for the generic flywheel and bearing system describe the bearings as having actual springs with a positive stiffness. The position stiffness terms for the magnetic force as derived in section 3.3 can be modeled as springs also, but have a negative stiffness effect on the system due to the instability of the actual system. With this in mind, the position stiffness term from equation 3-9 can replace the spring stiffness terms in the generic equations with the minor addition of a leading ( ) IIXIXIIXIXX FFFFtCosumFQ ???? ?+?+=?= 221121 ?? ( ) ( ) ( ) IIXIXIIXIX XX FFFF tCosdmLXKLXKXM ???? ?? ?+?+ =?+++ 2211 2 2211 ???? 28 negative sign. The current stiffness terms for each pole pair from equation 3-10 can be multiplied by the appropriate control current and substituted for the generalized forces for each bearing. The effects of gravity (G) along the Y axis can also be included and the current stiffness terms moved to the left side of the equation. After factoring the bias currents out of the current stiffness terms and rearranging, the equations of motion for a homogeneous flywheel with a mass imbalance as supported on horizontal active magnetic bearings can be seen below in equations 3-26. (3-26a) (3-26b) (3-26c) (3-26d) ( ) ( ) ( ) ( ) ( )tCosdmiiiiK iiiiKKLKLXKKXM IIXIIXIXIXX IIXIIXIXIXXXXXX cbcbi cbcbipppp ?? ? 2 21 22222 211112121 =?? ????+? ???? ???? ?? ( ) ( ) ( ) ( ) ( ) ( )GmMtSindmiiiiK iiiiKKLKLYKKYM IIYIIYIYIYY IIYIIYIYIYYYYYY cbcbi cbcbipppp ++=?? ???++? ???? ???? ?? ?? ? 2 21 22222 111112121 ( ) ( ) ( ) ( ) ( ) ( )tCosIIiiiiKL iiiiKLKLKLXKLKLII PTcbcbi cbcbippppPT IIXIIXIXIXX IIXIIXIXIXXXXXX ??? ???? 2 2 1 2 2 2 121 22221 111112121 ?=?+ ??+???? ???? ???? ??? ( ) ( ) ( ) ( ) ( ) ( )tSinIIiiiiKL iiiiKLKLKLYKLKLII PTcbcbi cbcbippppPT IIYIIYIYIYY IIYIIYIYIYYYYYY ??? ???? 2 2 1 2 2 2 121 22222 111112121 ?=?? ?++??++ ???? ???? ??? 29 CHAPTER IV CONTROLLING THE ACTIVE MAGNETIC BEARING RIG Many different control laws and algorithms are used to control mechanical systems. Since active magnetic bearings are inherently unstable, a controller is required such that the actuation of the system can be directed to achieve stability of the rotor and maintain this stability during operation. Consequently, the ability to monitor the "health" or performance of the system is readily available through the magnetic bearings. As discussed earlier in Chapter 2, many complex methods of control have been researched and implemented on active magnetic bearings. In this chapter, a linear state space approach with an appended integral control will be covered. First, amplifier dynamics will be derived and the overall state space of the plant will be discussed. A state estimator and observer feedback controller will then be designed using classical controls methods. Next, modeled perturbations to the system and integral control will be discussed. Adaptive disturbance rejection will be introduced and appended to the controller and control voltage splitting will be addressed, completing the model. 4.1 ? Amplifier Dynamics The equations of motion, as derived in Chapter 3, model the system as a self- contained plant that is excited via electrical current. However, most controller hardware 30 including the hardware used in this research only outputs control voltages. Amplifiers are needed to transform the control voltages into currents that will operate the bearings. These amplifiers, having their own dynamics, need to be modeled and appended to the system equations to complete the plant for numerical simulations. Amplifiers may be nonlinear in nature. However, for the purpose of simplicity, they can be approximated as first order filters. Equation 4-1 below shows a first order differential equation representing the electrical current-voltage dynamics of an amplifier. (4-1) In the above equation, ? refers to the amplifier bandwidth, Ka is the amplifier gain, up represents the controlling input voltage, and ic remains the control current. Each pole pair of the system, eight in all, has an amplifier assigned to it. Moreover, each pole pair receives its own control current, and recalling equation 3-10, also has its own assigned bias current. It is beneficial at this point to define a set of complementary pole pairs as receiving equal but opposite control voltages. Keeping the respective subscripts such that -I represents the primary pole pair and -II represents the complementary pole pair, equation 4-2 describes the control voltage for the left bearing along the X axis with respect to the individual voltages per pole pair. (4-2) The pole pair subscript is dropped in the first term of the above equation since it represents a single control voltage pertaining to both pole pairs. Each amplifier equation can now be multiplied through by its respective bias current, as shown in equations 4-3. pacc uKii =+ ?? IIXIXX ppp uuu ?? ?== 111 31 (4-3a) (4-3b) These two complementary amplifier equations are subtracted from each other, as seen in equation 4-4, allowing for a single dynamic equation per actuation axis. This effectively couples both amplifiers together. (4-4) Once again, the equation above is with respect to the control currents in the left bearing along the X axis. The same process can be done with respect to the 2X, 1Y, and 2Y axes. A similarity in the format of the control and bias currents can now be seen between the equations of motion derived in Chapter 3 and the above amplifier equations. With that in mind, a state space formulation can be presented that includes the orthogonal translations, rotations, velocities, and currents. 4.2 ? Linear State Variable Model The equations of motion of a flywheel supported by active magnetic bearings with amplifier dynamics are of a formulation that can be combined into a state space. However, the imbalance and gravity terms as appended to the end of equations 3-26 will be neglected for now allowing for the modeled plant to represent a linearized time invariant homogeneous system. These imbalance terms will be added back later as perturbations to the system. A classical state space controls approach will be used here IXIXIXIXIXIX pbacbcb uiKiiii ?????? =+ ? 111111 ? IIXIIXIIXIIXIIXIIX pbacbcb uiKiiii ?????? ?=+ ? 111111 ? XIIXIXIIXIIXIXIXIIXIIXIXIX pbbacbcbcbcb uiiKiiiiiiii 11111111111 ?? ?? ? ? +=? ? ?? ? ? ?+? ? ?? ? ? ? ?????????? ??? 32 such that the linear state variable model can be expressed as a state equation and an output equation. These equations are shown below as equations 4-5 and 4-6 respectively. (4-5) (4-6) In the above equations, x represents the state space or vector field with respect to the linearly independent variables of the system. The term, up describes a vector of inputs into the system, and y defines the measured outputs of the system. Underscored variables designate vectors as such. The terms A, B, C, and D represent matrices of coefficients that linearly operate on the vector fields. In the case of this particular dynamic system, the term Dup will be neglected due to the fact that the inputs play no part in the output of the system. A block diagram depicting the open loop linear state variable model is shown in figure 4-1. Figure 4-1 ? Block Diagram of an Open Loop Linear State Variable Model The state space and input vectors for the dynamic system as described above are shown in equations 4-7. ? x B up y A C x puBxAx += ? puDxCy += 33 (4-7) The structure and values for the linear operators from equations 4-5 and 4-6 with respect to the state space and input can be found in Appendix A of this thesis. 4.3 ? State Estimation In order to model the active magnetic bearings and rotor as realistically as possible, it should be noted that only the positions can be measured as outputs of the system, hence the need for the output equation. Since the full state can not be recovered, a state estimator will be required to estimate those states that are not measured. The state estimator will be designed such that it is a dynamic system whose purpose is to estimate the full state of the system given the measured outputs of the system and the inputs as supplied by the controller. The estimator is configured such that its output and the actual measured outputs converge. However, before the estimator can be built and implemented, the state variable model must first undergo an observability test proving ? ? ? ??? ? ? ? ? ??? ? = ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? = ???? ???? ???? ???? ? ? ? ? Y Y X X IIYIIYIYIY IIYIIYIYIY IIXIIXIXIX IIXIIXIXIX p p p p p cbcb cbcb cbcb cbcb u u u u u iiii iiii iiii iiii Y X Y X x 2 1 2 1 2222 1111 2222 1111 ? ? ? ? 34 that all states are indeed observable allowing for the eigenvalues of the estimator to be chosen and the estimation to converge with the plant. The observability matrix checks the structure and conditions of the linear operators A and C. For a state space containing n states or variables, the rank of the observability matrix must be equal to n in order for the model to be deemed observable. Equation 4-8 below shows the construction and criteria of the observability matrix. (4-8) After the model is determined to be observable, the state estimator can be built and implemented into the model. The design of the estimator as presented below in equation 4-9 is much like the state equation with the addition of the error term . (4-9) The variables designated with the circumflex (^) are defined as the estimated states. By substituting equation 4-6 sans the input term Dup into the above equation for the estimated output, equation 4-9 can be rearranged into a more usable form as given in equation 4-10. A block diagram depicting the state estimator can be seen in figure 4-2. (4-10) n CA CA CA C rank n = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?1 2 M ?????? ?++= ?? ? ? yyLuBxAx p ( ) yLuBxLCAx p ++?= ? ? ? ?????? ? ? yyL 35 Figure 4-2 ? Block Diagram of State Estimator The error (e) between the estimator output and the measured states is defined as shown in equation 4-11 below. By noting the closed loop behavior of the estimator, the error dynamics, or instantaneous change in error, are determined as shown in equation 4- 12. (4-11) (4-12) The term L is a matrix of coefficients chosen such that the eigenvalues of the matrix A-LC contain all negative real parts. Since the system was deemed observable, all eigenvalues of this matrix can indeed be placed anywhere, thus allowing the error dynamics to be driven to zero. Hence, a convergence between the estimated and measured output will be attained. With the state estimator built, the full state can now be recovered and an observer feedback control law defined. A-LC B ? xu y L ? ? x ? ?= xxe ( )eLCAe ?=? 36 4.4 ? Observer Feedback Control In this section, a control law will be defined such that the output of the state estimator or observer will be operated on and fed back into the plant via the input. Normally the input vector to the plant tracks some desired reference. In this particular case however, the desired reference is a zero vector and thus will not be included as it would have no effect on the closed loop behavior of the plant. Equation 4-13 below shows the form of the control law. (4-13) The term K as seen above represents a matrix of coefficients that operate on the estimated state. Before the control law can be implemented into the system, effectively closing the plant loop, the structure of the plant must first be checked for controllability. Much like the observability of the plant as discussed earlier, the controllability of all states are checked by discerning the rank of a matrix containing a combination of the state matrix (A) and input matrix (B) called the controllability matrix. For a plant containing n states, the controllability matrix must be full rank where the rank is equal to n. The structure of the controllability matrix and its criteria are shown below in equation 4-14. (4-14) After the plant is proven controllable, the control law can now be substituted into the model plant from equation 4-5 resulting in equation 4-15. (4-15) ? ?= xKu ?? ?= xBKxAz [ ] nBABAABBrank n =?12 L 37 Keeping in mind the influence of estimation error and equation 4-11, the estimated state can be represented as the difference between the actual state and the error. Equation 4-16 shows the substitution of the estimation error, and equation 4-17 gives the overall dynamics of the closed loop system. (4-16) (4-17) By examining the triangular structure of the overall closed loop system matrix above, the separation principal can be instituted such that the feedback gain (K) and the observer gain (L) may be designed separately. Thus, similarly to the placement of the eigenvalues for the state estimator, the eigenvalues or poles of the closed loop plant can now be chosen such that stability of the plant can be easily achieved. This of course is done by picking appropriate values for the feedback gain matrix K. Figure 4-3 depicts the closed loop system thus far. Figure 4-3 ? Block Diagram of Closed Loop Plant with Observer Feedback puBxAx += ? xCy = ( ) yLuBxLCAx p ++?= ? ? ? y u ? x??= xKu p ( ) ( ) eBKxBKAexBKxAx +?=??=? ( ) ( ) ?? ? ?? ? ?? ? ?? ? ? ?= ?? ??? ?? ??? ? ? e x LCA BKBKA e x 0 38 4.5 ? Perturbation and Integral Control The linear state variable model up to this point has utilized the equations of motion for a homogenous flywheel supported by active magnetic bearings and the amplifiers required to actuate the system. In this section, the imbalance and gravity terms that were dropped earlier will be added back to the modeled plant as input disturbances and any steady state error as seen by the controller will be addressed via integral control. 4.5.1 ? Gravity and Mass Imbalance As discussed earlier in Chapter 3, the forces due to imbalance can be modeled as sinusoidal disturbances to the system. The effects of gravity can be modeled as a constant function. These perturbations to the system are added to the model through a disturbance vector. Thus, by adding this disturbance vector to the state equation, the imbalance and gravity terms will be reinstated and equations 3-26 represented in their full form. Equations 4-18 and 4-19 will now serve as the linear variable plant model. (4-18) (4-19) The model as presented above is the same as equations 4-5 and 4-6 with the addition of the term ud representing an input disturbance vector and ? symbolizing a real valued matrix that maps the disturbance vector into the plant. dp uuBxAx ?++= ? xCy = 39 4.5.2 ? Integral Control Now that gravity has been included into the plant, a steady state offset from the "zero" state will occur. The observer based feedback controller as derived in section 4.4 acts like a proportional-derivative controller in that it only operates on the positions, velocities, and currents of the observer output. Thus, any steady disturbance will impact the steady-state dynamics of the system. One approach to alleviating this problem was addressed by Flowers et al. [7] and is utilized in this research. This is done by appending an integral control block after the feedback gain matrix K. Equation 4-20 shows the structure of the integral control where Kg represents the integral gain. (4-20) Thus, for a balanced levitating rotor, the position can now be stabilized and controlled such that it remains in a centered position with respect to the bearings. 4.6 ? Adaptive Disturbance Rejection The control strategy covered in this section is an adaptive controller specifically developed to reject persistent excitations at known frequencies. Originally developed by Fuentes et al. [8], this strategy was implemented successfully on active magnetic bearings and covered in great detail by Matras et al. [15 ,16, and 17]. Hence, the control laws introduced here will only be briefly explained. Given a disturbance vector ud as seen in equation 4-18, which contains some linear combination of scalar functions, amplitudes, and unit vectors, a control law (up) can be defined using adaptive techniques such that the effects of the disturbances will be dtuKuu pgpp ? += 40 suppressed. The scalar functions represented in the disturbance vector can be constant functions, sinusoidal functions, or any other wave form with a known phase. However, for sinusoidal disturbances, like imbalance in this case, the phase does not have to be known since it can be represented by two sinusoids that are 90 degrees out of phase with each other. The amplitudes of these disturbances do not have to be known and can change over time since the control technique adapts itself to match the amplitude of the excitations. Equation 4-21 below defines a vector of known or expected disturbances. (4-21) Given that the linear model is at least output stabilizable and is almost strictly positive real (A.S.P.R.), two positive definite weighting matrices can be defined, ?G and ?H. Since the model in this case is indeed controllable, the first condition is satisfied. The second condition requires feedback of the position, velocity, and current states if the amplifiers are first order systems and the sensor dynamics negligible. This condition as pertaining to active magnetic bearings is shown in Matras [15] and extended further in [17]. A control law can now be defined as seen in equation 4-22 with the adaptive terms defined in equations 4-23 and 4-24. (4-22) (4-23) (4-24) ? ? ? ? ? ? ? ? ? ? ? ? = i d ? ? ? ? M2 1 dppp HyGu ?+= G T p yyG ??= ? H T dp yH ??= ? ? 41 The terms Gp and Hp are adaptive gains that operate on the output of the plant. The gain Gp can be thought of as a stabilizing gain such that any nonzero output from the plant will cause this gain to increase. The adaptive gain Hp scales the disturbance vector ?d to the necessary amplitudes such that the disturbances as introduced to the plant are effectively canceled out. The weighting matrices, ?G and ?H, control how quickly these gains can adapt. Typically, the gain Gp is not necessary unless Hp becomes bounded before the perturbations are rejected. In this case, Gp will adapt until Hp reaches sufficient amplitude to suppress the excitation. When used in parallel with the observer feedback control law, this technique will produce asymptotic output tracking with bounded adaptive gains thereby canceling out any known disturbances during operation. Work done by Matras [17] showed that the adaptive disturbance rejection as presented above proved effective for frequencies below the critical speed. However, he heuristically demonstrated that by using the observer output as an input to , frequencies above the critical speed could be suppressed as well. Both techniques are used in this research. The numerical simulations in the next chapter and the first experimental rig use the approach developed herein. The latter experimental rig uses this extended approach. 4.7 ? Control Voltage Splitting It is beneficial at this point to address some specifics of the active magnetic bearing rig. The rigs used in this research, as stated previously, only detect position at the ends of the rotor. This is done using proximitor probes along the orthogonal axes of the bearings. Thus, four position readings are available with respect to the 1X, 2X, 1Y, pG ? 42 and 2Y axes and serve as outputs of the plant. Since there are eight pole pairs, four per bearing, used to actuate the rotor, there are eight control voltages going to the appropriate amplifiers and serving as inputs. Recalling equation 4-7, the state space representing the plant has modeled the system up to this point as having four control inputs and four outputs with twelve states. Modeling the plant in this way allows for the system to be deemed observable and controllable. However, since there is a discrepancy in the defined number of control inputs to the modeled plant and real magnetic bearing, one more block must be appended to the closed loop. This block will split the four control voltages as output from the control laws into eight which will then be input into the modeled amplifiers. This allows for a more realistic plant model. Thus, the state space as defined in equation 4-7 can be thought of as the "modeled plant", and the state space as defined below in equation 4-25 can be thought of as the "real plant". The vector ud will only apply to the "real plant". (4-25) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ?? ? ? ?? ? ? ? ?? ? ? ?? ? ? ? + ? ?= ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? = ?? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? GmM tII tII tm tm u u u u u u u u u u i i i i i i i i Y X Y X x PT PTd p p p p p p p p realp c c c c c c c c real IIY IY IIY IY IIX IX IIX IX IIY IY IIY IY IIX IX IIX IX ? ? ? ? ? ? ? ? cos sin cos sin 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 43 Although the general state space structure of the "model plant" and "real plant" differ, the overall internal dynamics of both plants will remain the same. The main difference between the two plants other than the addition of the disturbance vector is the placement of the bias currents. The bias currents in the "model plant" were contained within the state space. For the "real plant" however, the bias currents with regard to the position and velocity equations of the magnetic bearing and rotor will remain a part of the current stiffness terms and reside in the "real" state matrix, Areal. Moreover, the amplifier equations of the "real plant" will not be multiplied by the bias currents as shown in equation 4-4. The amplifier equations will retain the form of equation 4-1 and will be effectively decoupled with respect to the complementary pole pairs. This will allow for eight individual control inputs. The "modeled plant" is still useful such that it will be used in the state estimator and the control laws will operate upon it as described earlier. However, each of the four control voltages coming out of the integral controller will now be split into two signals. The first signal will be added to a bias voltage such that the voltage to the amplifier of the primary pole pair will remain at a constant value unless excited by the controller. Recalling equation 4-2, the second signal will be subtracted from a bias voltage and sent to the amplifier of the secondary pole pair. For example, if the rotor needs to be moved up along an axis, the upper pole pair will see a boost in current and the lower pole pair will see a drop in current. Bias voltages are necessary since a negative voltage will have the same effect as a positive voltage with respect to the amplifiers and the magnetic force produced by the pole pair. Hence, the control voltages should not drop below zero nor should they reach 44 so high of a value that the amplifiers overload. Equations 4-26 show two general control inputs to the "real plant" as they are split from a single control signal. This will apply to all four controller outputs. (4-26a) (4-26b) In the equation 4-26, vb represents the bias voltage. Keeping in mind that the currents output by the amplifiers are dependent on the voltages sent to them, the bias currents can now be defined as a function of the bias voltages. Equation 4-27 shows the relationship between the bias current and voltage per amplifier. (4-27) A block diagram depicting the completed linear state space model encompassing the "model plant" with observer feedback, adaptive disturbance rejection control, and integral control as fed into a "real plant" can be seen in figure 4-4. pbp uvu II += pbp uvu IIII ?= bab vKi = 45 Figure 4-4 ? Block Diagram of Complete Closed Loop System In the above figure, the dotted line connecting the output of the observer to the adaptive disturbance rejection block is a simplified demonstration of the extension as performed by Matras [17]. As stated previously, this connection will only apply to the experimental results of the second test rig. With the plant loop now closed and completed, numerical simulations can take place and experimental rigs utilized. ( ) yLuBxLCAx p ++?= ? ? ? dprealreal uuBxAx ?++= ? )()( xCy real )(= ? = xKu p y ud ? x dtuKuu pgpp ? += dppp HyGu ?+= pbp pbp uvu uvu IIII II ?= += up Real Plant Model Plant d? 46 CHAPTER V NUMERICAL SIMULATIONS AND RESULTS The numerical simulations and results found in this chapter are a product of the active magnetic bearing equations of motion as derived in Chapter 3 and the control algorithms developed in Chapter 4. State space matrices and other parameters are defined and setup via MatLab m-files and block diagrams representing the control algorithms are created in Simulink. Once the system is setup within the program, initial conditions are defined and numerical integration can commence, hence producing the simulation. The purpose of the simulation is to represent the active magnetic bearing system as accurately as possible. The results from these simulations can give a large insight into the characteristics of the real system and help to decide optimal values for certain parameters essential to the operation of the real bearings. Although one should keep in mind that the simulations are indeed derived from a linearized plant, they are a very good indicator of what to expect from the actual magnetic bearings. In this chapter, graphical results and conclusions will be presented pertaining to the imbalance of a flywheel and the characteristics of the adaptive disturbance rejection controller. These simulated results will be validated by comparing them to the experimental results as produced by two different test rigs. The m-files and block diagrams utilized in this chapter can be found in Appendix A. 47 5.1 ? Adaptive Gain, Hp Recalling section 4.6 in the previous chapter, there are two adaptive gains, Gp and Hp, which constitute the adaptive disturbance rejection control. Since the imbalance amplitudes herein will be generally small, the gain Gp will not be used and its weighting matrix, ?G, set to zero. Thus, focus will remain on Hp. As seen in equation 4-26, the derivative of this adaptive gain is comprised of the disturbance vector ?d, the output of the plant, and a weighting matrix. Since ?d is constructed using two sinusoidal functions that are 90o apart, the gain Hp has a Sine and Cosine term for each of the four output vectors of the plant. It will be assumed at this point that the imbalance will only affect the translational components of the plant, X and Y. Although an imbalance can most certainly affect the rotational components, ? and ?, focus will remain on the translational imbalance since the results from that can be easily extended to a rotational case. Recalling section 3.4.3, a translational imbalance will affect both orthogonal translations of the rotor. Hence, it will only be necessary to examine one output vector, in this case X. Moreover, since the phase of the imbalance will not be known, it is imperative that the magnitude of the sum of the Sine and Cosine terms of Hp be observed. Equation 5-1 shows how this is set up. 5-1 Changes in different parameters can affect the outcome of the adaptive imbalance detection gain, H*p. Some of these parameters include the placed poles of the closed loop 22 CospSinpp HHH += ? 48 system, the running speed, imbalance mass, distance the mass travels, ?H value, noise level in the system, and discrepancy between the running speed and the frequency to be rejected, as well as many others. Although the magnitude, frequency, and the general shape of H*p can differ, the overall result is generally the same such that any increase in imbalance will cause either an increase in H*p or a discontinuity in the slope of H*p. 5.2 ? Simulation Results In this section, the parameters previously mentioned will be varied and results presented to give an overall idea of how they affect the system and the detection of an imbalance. Initially, a mass of one gram at a distance of two millimeters away from the center of the flywheel with a given running speed will make up the imbalance. At some time, t, the mass will ramp to a distance of four millimeters away from the center of the flywheel in 0.1 seconds providing an increase in imbalance and thus, simulating a crack growth in a flywheel. These values, with exception of the running speed, will remain constant throughout all of the simulations. First, ideal conditions should be considered. 5.2.1 ? Ideal Conditions Ideal conditions for this type of simulation can be defined as no noise and an exact match in the running speed and the frequency to be rejected by the ADR. Beginning with a running speed of 20 Hz, a time trace and H*p plot can be found on the next page along with a table of given parameters. The parameters defined in table 5-1 will be the baseline parameters for the simulation results shown in this chapter. 49 0 2 4 6 8 10 12 14 16 18 20-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 10 -7 Time (s) Tr an sla tio n i n X (m ) Figure 5-1 ? Time Trace for Ideal Conditions Running Speed (?) 20 Hz ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Plant Poles -1000 + 250i Table 5-1 ? Parameters for Ideal Conditions 50 0 2 4 6 8 10 12 14 16 18 200 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 -3 Time (s) Hp * Figure 5-2 ? H*p for Ideal Conditions It should be noted at this point that the simulations begin with the rotor already spinning, hence the observed transient response for the first 2 seconds or so. Regardless, at 10 seconds, it can be plainly observed in the time trace the effect of the increase in the imbalance. The adaptive imbalance detection gain also shows a dramatic increase in amplitude almost instantaneously. Since the running speed and the frequency to be rejected match exactly, this curve becomes a straight line. It will be shown in later results that as the two frequencies drift apart, this curve will become sinusoidal. However, the dramatic increase in amplitude will still be clearly noticeable. 51 5.2.2 ? Noisy Conditions In figure 5-1 above, the imbalance is clearly visible. However, the changes in imbalance are typically small and may be completely hidden by the noise so that the increase in imbalance is indiscernible as shown in figure 5-3. Even for such conditions, the noise in the time trace does indeed show up in H*p. However, the change in imbalance is still observable. Figures 5-3 and 5-4 depict this condition using the same parameters as before with white noise added. The amount of white noise superimposed in the system is defined by the signal to noise ratio of H*p shown in table 5-2. Figure 5-5 shows figures 5-2 and 5-4 superimposed on each other to further illustrate how the added noise affects the adaptive imbalance detection gain, H*p. 0 2 4 6 8 10 12 14 16 18 20-3 -2 -1 0 1 2 3 x 10 -7 Time (s) Tr an sla tio n i n X (m ) Figure 5-3 ? Time Trace for Noisy Conditions (Case I) 52 Running Speed (?) 20 Hz ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p 6.3303 Plant Poles -1000 + 250i Table 5-2 ? Parameters for Noisy Conditions (Case I) 0 2 4 6 8 10 12 14 16 18 200 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 -3 Time (s) Hp * Figure 5-4 ? H*p for Noisy Conditions (Case I) 53 0 2 4 6 8 10 12 14 16 18 200 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 -3 Time (s) Hp * Ideal Noisy Figure 5-5 ? Comparison of H*p for Ideal and Noisy Conditions (Case I) Since there are a variety of ways to discern the signal to noise ratio, it is important that it is explained how the value is derived. It is a simple procedure in which the average value of H*p is determined after a time when the initial transients have died down. This time is found by examining the case with no noise. Here, the average is taken between 6 and 10 seconds, so that not only the initial transients are neglected but so are the effects of the increase in imbalance. Then, the largest deviation is found and subtracted from that average. This value is considered the noise. The average amplitude of H*p is then divided by this noise to give a signal to noise ratio. 54 Of course, as the signal to noise ratio decreases, so does the visual value of the gain. The noise eventually saturates this gain until no useful information can be extracted from it. Moreover, using an ideal condition with no noise, the increase in H*p is 100 % assuming that the distance away from center the mass moves to is double the original distance. This can be seen in the previous plots. However, as the noise increases, this percentage decreases. In fact, after a certain signal to noise ratio has been reached, the only way to discern any information is by taking the averages of the gain before and after the imbalance occurs, neglecting the transients of course. Next, plots for the time trace and H*p with increased noise will be given. A plot of the percentage increase in H*p versus signal to noise ratio for a simulated running speed of 20 Hz will also be presented. These graphs depict how the increase in noise affects the curve. It is interesting that until a low signal to noise ratio is reached, a significant increase can still be ascertained. The percentage increase is found by subtracting the average of the curve before the increase in imbalance occurs from the average after the increase and dividing by the initial average. This value is then multiplied by 100. 55 0 2 4 6 8 10 12 14 16 18 20 -5 -4 -3 -2 -1 0 1 2 3 4 5 x 10 -6 Time (s) Tr an sla tio n i n X (m ) Figure 5-6 ? Time Trace for Noisy Conditions (Case II) Running Speed (?) 20 Hz ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p 0.5272 Plant Poles -1000 + 250i Table 5-3 ? Parameters for Noisy Conditions (Case II) 56 0 2 4 6 8 10 12 14 16 18 200 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 -3 Time (s) Hp * Figure 5-7 ? H*p for Noisy Conditions (Case II) Inspection of figure 5-7 shows that one can barely discern that the minimum values after the imbalance increase are indeed higher than before the imbalance. Although this graph is indeed extremely noisy, the percentage increase is still nearly 56%. Of course, the signal to noise ratio is very low, 0.5272 in this particular case. The next plot shows the percentage increase as a function of the signal to noise ratio. Many simulations were performed in which the noise level was varied. As the noise decreases and the signal to noise ratio increases, the curve goes to 100%. However, the percentage increase as well as the visual importance of H*p drops off significantly 57 after a signal to noise ratio of 2, and plummets after a ratio of 1. This is certainly an expected trend. 0 1 2 3 4 5 6 710 20 30 40 50 60 70 80 90 100 Signal to Noise Ratio Pe rce nta ge In cre as e i n H p* Figure 5-8 ? Percentage Increase in H*p vs. Signal to Noise Ratio for ? = 20 Hz All of the plots given thus far are for a running speed of 20 Hz. As would be expected, an increase in running speed below the first critical speed causes an increase in imbalance and thus an increase in the values of the adaptive imbalance detection gain. This can be seen for the case where the running speed and ADR frequency are set to 40 Hz. To accentuate the point of detecting imbalance in a noisy system, white noise will be added such that the imbalance will not be visible in the time trace. However, the H*p plot will show how the increase of imbalance is plainly visible and detectable. 58 0 2 4 6 8 10 12 14 16 18 20-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 10 -6 Time (s) Tr an sla tio n i n X (m ) Figure 5-9 ? Time Trace for Noisy Conditions (Case III) Running Speed (?) 40 Hz ?H 40,000 ADR Frequency 40 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p 8.8787 Plant Poles -1000 + 250i Table 5-4 ? Parameters for Noisy Conditions (Case III) 59 0 2 4 6 8 10 12 14 16 18 200 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 -3 Time (s) Hp * Figure 5-10 ? H*p for Noisy Conditions (Case III) 5.2.3 ? Pole Variation Changing the placement of the controller poles will of course change the system response. Thus, the behavior of H*p is changed as well. For the cases described above, the poles were placed at -1000+250i. This set of poles gives a quick response to the system. Decreasing the magnitude of the real part of the poles will decrease the time it takes for the transients to die down. The natural frequency of the system is also altered since the real part of the poles defines the critical speed in radians per second. To show how changing the poles affect the system, the following plots show the system running at 60 20 Hz with poles placed at -500+250i with no noise. This produces a visible difference from the case presented in figure 5-1. 0 2 4 6 8 10 12 14 16 18 20 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 10-7 Time (s) Tr an sla tio n i n X (m ) Figure 5-11 ? Time Trace for System Poles = -500 + 250i Running Speed (?) 20 Hz ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Plant Poles -500 + 250i Table 5-5 ? Parameters for System Poles = -500 + 250i 61 0 2 4 6 8 10 12 14 16 18 200 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 -3 Time (s) Hp * Figure 5-12 ? H*p for System Poles = -500 + 250i Comparing the time traces from figures 5-1 and 5-11, one can see that with the decreased magnitude of the real parts of the poles, the response is quicker yet the amplitude of the transients and imbalance is higher. This extends to H*p in figures 5-2 and 5-12. Likewise, decreasing the magnitudes of the real parts of the poles causes the ADR to respond and settle quicker to the imbalance change. Numerous simulations were done at 20 Hz to determine the settling time of the imbalance change as a function of different pole placements. The settling time is calculated by finding the time at which H*p reaches 0.1% of the average value after the imbalance has changed and the transients die down. This value is then subtracted from the time of the initial imbalance change to 62 give a settling time. As can be seen, the trend does show that as the magnitude of the real parts of the poles decrease, so does the settling time. -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -2000 1 2 3 4 5 6 7 8 9 10 Real Part of the System Poles Se ttli ng T im e ( s) Figure 5-13 ? Settling Time of H*p as a Function of Pole Placement Another interesting effect of the variation in pole placement is the amplitude at which H*p settles. Since the poles affect settling time and the steady state response of the plant, it is expected that not only should the transients of H*p change, but also the steady state value. As the magnitude of the real parts of the poles is decreased, the amplitude of the steady state oscillations in the time traces decrease. This allows H*p to reach a higher value, effectively suppressing more of the vibration as permitted by the placement of the poles. It is until the placed poles and resulting critical speed reach a value closer to the 63 running speed that the amplitude of H*p starts to decrease. This is expected because once the critical speed matches the running speed, the system becomes unstable. These results give insight into the benefits of pole placement with regard to imbalance detection. However, particular attention should be paid to other aspects as well since decreasing the magnitude of the real parts of the poles does effectively soften the system. Once again, figures 5-12 and 5-13 are for a running speed and ADR frequency of 20 Hz and a weighting matrix value of 40,000. -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -2001 1.5 2 2.5 3 3.5 4 4.5 5 x 10 -4 Real Part of the System Poles Av era ge A mp litu de of H p* be for e I mb ala nc e Figure 5-13 ? Amplitude of H*p as a Function of Pole Placement 64 5.2.4 ? ?H Variation Changing the poles is not the only way to affect the ADR response. As previously stated in section 4.6, altering the weighting matrix, ?H, on which the adaptive gain, Hp, operates can also affect how quickly H*p responds to an imbalance. In the cases presented above, the weighting matrix has been defined with a baseline of 40,000. By increasing this value, not only will the response be quicker, but the general initial transient amplitude and frequency will be higher until the imbalance is suppressed. Because the adaptive response is quicker, the imbalance seen in the time trace will be slightly decreased as can be seen in the figures below. For this case, the weighting matrix was set to 100,000. 0 2 4 6 8 10 12 14 16 18 20-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 10 -7 Time (s) Tr an sla tio n i n X (m ) Figure 5-15 ? Time Trace for ?H = 100,000 65 Running Speed (?) 20 Hz ?H 100,000 ADR Frequency 20 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Plant Poles -1000 + 250i Table 5-6 ? Parameters for ?H = 100,000 0 2 4 6 8 10 12 14 16 18 200 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 -3 Time (s) Hp * Figure 5-16 ? H*p for ?H = 100,000 66 As seen in figures 5-15 and 5-16 as compared with the alternate pole placement in figures 5-11 and 5-12, the time traces show a slight decrease in imbalance amplitude. Also, H*p responds similarly with respect to how quickly the imbalance is suppressed. Figure 5-16 also depicts the increase in frequency and amplitude of H*p during the initial transient stage. This gives insight to the dynamics of the controller as it works harder and faster to suppress the imbalance. The magnitude of the weighting matrix substantially influences the settling time of H*p as shown in figure 5-17. Numerous simulations were done where the magnitude of the weighting matrix was varied. The settling time was determined as explained in section 5.2.3. 0 1 2 3 4 5 6 7 8 9 10 x 105 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Delta H Se ttli ng T im e ( s) Figure 5-17 ? Settling Time of H*p as a Function of Increase in ?H 67 As seen in figure 5-17, an increase in the weighting matrix can drastically improve the response and settling time of H*p. As the weighting matrix value approaches infinity, the settling time approaches zero. However, the settling time can never truly be zero since the imbalance must occur before the ADR can take effect. Once again, all simulations contributing to figure 5-17 were performed for a running speed and ADR frequency of 20 Hz. 5.2.5 ? Discrepancies between Running Speed and ADR Frequency It is important to note the effects of discrepancies between the running speed and the frequency to be rejected. In this research as well many other cases, it is extremely difficult to synchronize the two frequencies. Although motor speeds are usually known within some range, they are not exactly known most of the time. For all previously presented simulations, the running speed and ADR frequency have been identical. However, realistic discrepancies between the two will now be investigated. Detecting imbalances using the adaptive imbalance detection gain, H*p, is a very robust procedure with regards to differing frequencies. Although the shape of H*p is altered, the imbalance is still obvious and very detectable within a certain range. In the previous sections, H*p settles to what appears to be a straight line. However, when there is a difference between the running speed and the frequency to be rejected, H*p will contain several frequencies as it suppresses the imbalance. The overall sinusoidal oscillation of H*p will exponentially decrease in amplitude until it rejects as much of the imbalance as it can. It will then settle to a bounded oscillation. This is shown in the 68 figure below where the running speed is 25 Hz and the ADR frequency is set at 20 Hz. Table 5-7 gives the parameters. 0 2 4 6 8 10 12 14 16 18 200 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 -4 Time (s) Hp * Figure 5-18 ? H*p for Discrepancy between Running Speed and ADR Frequency (Case I) Running Speed (?) 25 Hz ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Plant Poles -1000 + 250i Table 5-7 ? Parameters for Discrepancy between Running Speed and ADR Frequency (Case I) 69 0 10 20 30 40 50 60 70 80 90 100 110 1200 0.05 0.1 0.15 Frequency (Hz) Am pli tud e Figure 5-19 ? FFT of H*p between 10 - 15 Seconds for Discrepancy in Running Speed and ADR Frequency (Case I) In figure 5-19 above, a FFT analysis of H*p between 10 and 15 seconds is calculated to see which frequencies contribute to the sinusoidal nature of the transients. The first peak is at 5 Hz. This is the difference between the running speed and ADR frequency which makes up part of the exponentially decaying component or transient of H*p. The next peak is at 50 Hz, which is the 2X component of the running speed and finally there is a peak at 100 Hz, which is the 4X component. In the next plot, a FFT analysis is performed for H*p for the time between 15 and 20 seconds. This shows that the transient frequency has indeed died out leaving only the 2X and 4X components of the running speed, the steady state frequencies of the adaptive imbalance detection gain. 70 0 10 20 30 40 50 60 70 80 90 100 110 1200 0.05 0.1 0.15 Frequency (Hz) Am pli tud e Figure 5-20 ? FFT of H*p between 15 - 20 Seconds for Discrepancy in Running Speed and ADR Frequency (Case I) As can be observed, the 5 Hz frequency had died out completely. The mathematical reason for the appearance of these frequencies is that the ADR feeds back into itself through the output vector of the plant. The adaptive gain, Hp, can be divided into its sine and cosine parts with respect to ?d. Assuming a single X vector output by the plant can be characterized as a cosine imbalance term with ADR feedback and the weighting matrix is neglected, the equations governing the ADR would appear as shown below. (5-2) ( )tyH XSinp ?sin= ? 71 (5-3) (5-4) The term ?f in equation 5-4 represents the frequency difference between the running speed and the ADR as added to or subtracted from the running speed. As can be seen, the equations become very complex, especially when keeping equation 5-1 in mind as the final operation to derive H*p. Another simulation at a running speed of 50 Hz supports these results further. 0 2 4 6 8 10 12 14 16 18 200 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 -4 Time (s) Hp * Figure 5-21 ? ?pH for Discrepancy between Running speed and ADR Frequency (Case II) ( )tyH XCosp ?cos=? ( )( ) ( ) ( )( )tHtHty CospSinpfX ??? cossincos +???= 72 Running Speed (?) 50 Hz ?H 40,000 ADR Frequency 60 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Plant Poles -1000 + 250i Table 5-8 ? Parameters for Discrepancy between Running Speed and ADR Frequency (Case II) 0 20 40 60 80 100 120 140 160 180 200 2200 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Frequency (Hz) Am pli tud e Figure 5-22 ? FFT of H*p between 10 - 15 Seconds for Discrepancy between Running Speed and ADR Frequency (Case II) 73 0 20 40 60 80 100 120 140 160 180 200 2200 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Frequency (Hz) Am pli tud e Figure 5-23 ? FFT of H*p between 15 - 20 Seconds for Discrepancy between Running Speed and ADR Frequency (Case II) Once again, the difference in running speed and ADR frequency shows itself as a transient frequency of H*p which dies out over time leaving only the steady state frequencies, 2X and 4X of the running speed. Of course, as the ADR frequency moves away from the running speed or vice versa, the imbalance will be affected less. Hence the average steady-state value of H*p will be less as well. In figure 5-24, the plot shows how the discrepancy in frequencies affects the steady-state value of H*p as a percentage of the case with matching frequencies. 74 0 2 4 6 8 10 12 14 16 18 200 10 20 30 40 50 60 70 80 90 100 Difference between Running Speed and ADR Frequency (Hz) Pe rce nt of Hp * S tea dy S tat e A mp litu de Figure 5-24 ? Percentage Change in H*p as a Function of Discrepancy between Running Speed and ADR Frequency For figure 5-24, the value of the weighting matrix is 40,000. The ADR frequency was set to 20 Hz with the poles at -1000 + 250i. The running speed was then increased from 20 Hz to 40 Hz. It can be seen that as the difference between the two frequencies becomes greater, the steady-state amplitude of H*p exponentially decreases. However, even though this value decreases, the change in imbalance still causes H*p to increase 100% of it initial value assuming the imbalance distance has doubled. This can be seen in figures 5-18 and 5-21 above. It should also be mentioned that in the case where the running speed becomes less than the ADR frequency, the steady-state amplitude of H*p 75 still decreases. In fact, it decreases even faster due to the lesser amplitudes of imbalance with regard to the running speed. Noise has not been included into these plots for the purpose of clarity. However, as discussed in section 5.2.2, noise will certainly affect the outcome of H*p. Regardless, an imbalance should still be detectable within a certain range of noise and frequency discrepancy. 5.2.6 ? Introduction of Other Frequencies Another case to be considered is when other frequencies close to the running speed are introduced into the system. The ADR is inevitably affected by any and all frequencies within range of the running speed. When this happens, H*p takes on sinusoidal shapes. If one were to try to reject an imbalance at a known frequency and another sinusoidal excitation of a different frequency were introduced, the adaptive imbalance detection gain would most certainly react to both. One realistic example of this case is when using a controller to simulate an imbalance on a test rig and spinning the rotor as well. Since the rotor speed in some instances can only be approximated, there will be two excitation frequencies, the spinning rotor and the controller simulated imbalance. This particular example will be supported by results from an experimental rig in the next chapter. The results of a simulation where the running speed and ADR frequencies are set to 20 Hz are shown in the following pages. A secondary sinusoidal excitation of one third the amplitude is introduced at a frequency of 25 Hz. 76 0 2 4 6 8 10 12 14 16 18 200 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 -3 Time (s) Hp * Figure 5-25 ? H*p for Two Excitation Frequencies (Case I) Running Speed (?) 20 Hz Initial Imbalance Level 0.0316 N ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Secondary Excitation Frequency, Amplitude 25 Hz 0.01 N Plant Poles -1000 + 250i Table 5-9 ? Parameters for Two Excitation Frequencies (Case I) 77 0 10 20 30 40 50 60 70 80 90 1000 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 Frequency (Hz) Am pli tud e Figure 5-26 ? FFT of H*p for Two Excitation Frequencies (Case I) By examining the FFT results shown above, a peak at 5 Hz is readily observable. This is once again the difference between the rejected frequency and the introduced frequency. However, the smaller peak at 45 Hz is not the 2X component of the running speed or the secondary excitation. In fact, the results are even less expected. This smaller peak is the difference between the rejected and secondary frequencies subtracted from the 2X component of the secondary excitation (2 * 25 Hz - 5 Hz = 45 Hz). By performing a number of simulations, it has been deduced that where the frequency at which that peak will occur depends on the difference between the respective frequencies. Should the secondary frequency be larger than the rejected frequency, the difference will be subtracted from the 2X component of the secondary frequency. However, if the 78 frequency of the secondary excitation is less than the running speed, then the difference will be added to the 2X component of the secondary frequency. Below is an additional case at different speeds to support this supposition. The running speed and ADR frequency are set at 60 Hz and the secondary excitation has a frequency of 50 Hz. 0 2 4 6 8 10 12 14 16 18 200 1 2 3 4 5 6 7 8 9 x 10 -3 Time (s) Hp * Figure 5-27 ? H*p for Two Excitation Frequencies (Case II) H*p as shown above appears to settle to a nearly straight line. However, it should be noted that as the secondary frequency moves away from the running speed, the less contribution it will have on H*p. Table 5-10 and figure 5-28 show the parameters and FFT analysis results for the adaptive imbalance detection gain shown in figure 5-27 above. 79 Running Speed (?) 60 Hz Initial Imbalance Level 0.0316 N ?H 40,000 ADR Frequency 60 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Secondary Excitation Frequency, Amplitude 50 Hz 0.01 N Plant Poles -1000 + 250i Table 5-10 ? Parameters for Two Excitation Frequencies (Case II) 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 1500 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 Frequency (Hz) Am pli tud e Figure 5-28 ? FFT of H*p for Two Excitation Frequencies (Case II) 80 For case II presented above, the difference between the two frequencies is 10 Hz, which is exactly where the first peak appears. Also, since the secondary frequency was less than the rejected frequency, the difference is added to the 2X component of the secondary imbalance to produce the second peak (2 * 50 Hz + 10 Hz = 110 Hz). For secondary excitations with relatively low amplitude and frequencies that are a good distance away from the frequency to be rejected, imbalance detection is easily observable. However, as the two frequencies become closer, the frequency difference becomes the dominating part of the adaptive imbalance detection gain. Figure 5-29 shows a frequency difference of 1 Hz. 0 2 4 6 8 10 12 14 16 18 200 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10 -3 Time (s) Hp * Figure 5-29 ? H*p for Two Excitation Frequencies (Case III) 81 Running Speed (?) 20 Hz Initial Imbalance Level 0.0316 N ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Secondary Excitation Frequency, Amplitude 21 Hz 0.01 N Plant Poles -1000 + 250i Table 5-11 ? Parameters for Two Excitation Frequencies (Case III) Even with a difference of only 1 Hz, the imbalance can be detected. However, as the two frequencies approach each other, the frequency difference saturates H*p. Also, in a case where the secondary imbalance amplitude is significantly larger than the primary imbalance, detection can become even more difficult. When this happens, the only means of detection may be by checking for discontinuities in the slope of H*p. In figure 5-30, the frequency difference is 0.08 Hz. This means that H*p will oscillate at this frequency with the larger frequency components superimposed on this signal. Depending on when the imbalance occurs during the cycle of H*p, the discontinuity may be easily visible or it may not. There is also a large dependency on the amplitudes of the primary and secondary imbalances. In the case presented in figure 5-30, the secondary excitation has a magnitude of nearly three times the primary imbalance, yet detection is still possible. 82 0 2 4 6 8 10 12 14 16 18 200 0.5 1 1.5 2 2.5 3 x 10 -3 Time (s) Hp * Figure 5-30 ? H*p for Two Excitation Frequencies (Case IV) Running Speed (?) 20 Hz Initial Imbalance Level 0.0316 N ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Secondary Excitation Frequency, Amplitude 20.08 Hz 0.1 N Plant Poles -1000 + 250i Table 5-12 ? Parameters for Two Excitation Frequencies (Case IV) 83 The equations governing the ADR under these circumstances will be slightly different from those presented in section 5.2.5. However, the assumptions made will remain the same. In this case equations 5-2 and 5-3 will remain the same, but equation 5- 4 will have an appended sinusoidal term. This is shown in equations 5-5, 5-6, and 5-7 below. (5-5) (5-6) (5-7) A similar result can be obtained by adding more frequencies for the ADR to reject. In instances where motor speed is not known exactly, one might simply add more frequencies within the approximate range of the motor speed. This, however, proves more harmful since the additional ADR frequencies act similarly to secondary imbalance frequencies. The more frequencies added, the worse H*p can become. Each added frequency within a close range of another frequency works to suppress part of the vibration as seen by the output. Hence, those added frequencies also become part of the closed loop and work on each other. This can be seen in figure 5-31 below where a secondary ADR frequency was added. The parameters can be found in table 5-13. ( ) ( )( )( ) ( ) ( )( )tHtHtty CospSinpffX ????? cossincoscos +?+??+= ( )tyH XSinp ?sin=? ( )tyH XCosp ?cos=? 84 0 2 4 6 8 10 12 14 16 18 200 0.5 1 1.5 2 2.5 3 x 10 -3 Time (s) Hp * Figure 5-31 ? H*p for Secondary ADR Frequency Running Speed (?) 19.8 Hz ?H 40,000 ADR Frequencies 20 Hz 20.2 Hz Time of Imbalance Increase 10 s Signal to Noise Ratio for H*p Approaching ? Plant Poles -1000 + 250i Table 5-13 ? Parameters for Secondary ADR Frequency 85 In the instance of having two or more ADR frequencies, the governing equations become even more complex since there are now two ADR components for each frequency to be rejected. These are shown below where the adaptive terms Hp Sin(?) and Hp Cos(?) refer to the primary frequency to be rejected and the terms Hp Sin(?+?) and Hp Cos(?+?) refer to the second frequency to be rejected . (5-8) (5-9) (5-10) Looking back at figure 5-31, one can see again that the difference in the frequencies (in this case 0.2 Hz) dominates H*p. Noticing the discontinuity in the slope of the curve at 10 seconds is the only way one can discern that the imbalance occurred. It is stated in the basic theory of the ADR method that the frequency of the sinusoidal disturbance must be fully known for the controller to have full effect. However, as shown, H*p is very useful in detecting increases in imbalance when used correctly. Typically, if the running speed is not fully known, it is best to set only one ADR frequency within some range. Although the imbalance will not be suppressed fully, one will be able to detect an increase in amplitude of that frequency. It should be noted that all of these results are for constant rotor speeds. It is obvious that since H*p is very responsive to a wide range of frequencies as well as its own, speed shifts will also affect the ADR dramatically. Depending on how close the ( ) ( )tyH XSinp ?? sin= ? ( ) ( )tyH XCosp ?? cos= ? ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( )( )?? ? ? ? ?? ? ? ? ??++ ??+ ?= ?? ?? tHtH tHtH ty fCospCosp fSinpSinp X ?? ?? ? ?? ?? coscos sinsin cos 86 ADR frequency and running speed are, one can expect the dominate frequency of H*p to decrease or increase as the speed shifts. One way to overcome the sinusoidal tendency of H*p is to feed the fully known running speed back into the controller. Should the speed increase, the ADR will compensate for this type of increase in imbalance. However, when a mass shifts or a crack develops, there should be a significant change in H*p. In the figures below, the ADR frequency and running speed is a perfect match. The running speed is initially 20 Hz and then ramps up to 30 Hz between 10 and 15 seconds. At 12.5 seconds the imbalance distance is increased. Both the time trace and H*p are given as well as the parameters for the simulation. 0 2 4 6 8 10 12 14 16 18 20-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 10 -7 Time (s) Tr an sla tio n i n X (m ) Figure 5-32 ? Time Trace for Increase in Running Speed and ADR Frequency 87 Initial Running Speed (?) 20 Hz Final Running Speed (?) 30 Hz Ramp Time 10 - 15 s ?H 40,000 ADR Frequency Matches ? Time of Imbalance Increase 12.5 s Signal to Noise Ratio for H*p Approaching ? Plant Poles -1000 + 250i Table 5-14 ? Parameters for Increase in Running Speed and ADR Frequency 0 2 4 6 8 10 12 14 16 18 200 0.5 1 1.5 2 2.5 3 x 10 -3 Time (s) Hp * Figure 5-33 ? H*p for Increase in Running Speed and ADR Frequency 88 As can be observed in the above figures, H*p compensates for the increase in imbalance due to an increase in speed. However, when the imbalance mass is actuated, the curve responds dramatically until the speed reaches a constant value again. This result is certainly expected. By knowing the running speed exactly, one can foresee the change in H*p. Moreover, if a crack was to occur, then it should be detectable. 5.3 ? Response Time Generally, the response of the ADR to a change in imbalance is instantaneous. This, of course, is practically dependent on the specific situation of the system. For ideal conditions, H*p will begin to increase as soon as a change in the amplitude of some sinusoidal frequency is detected in the output of the plant. However, as discussed in previous sections, noise, frequency discrepancy, and added frequencies can have a large effect on how H*p responds. Ideally, one would like for the imbalance to be detectable within a single revolution of the rotor and, for ideal conditions, this is the case. By inspection of figure 5-29, it can be observed that the oscillation of H*p can obscure the immediate visibility of the increase although it is still for the most part instantaneous. Visibly identifying an increase in imbalance under these sinusoidal circumstances is largely dependent on when the imbalance happens in the cycle of the adaptive imbalance detection gain. However, noticing a discontinuity in the slope of H*p is a sure sign that conditions have indeed changed. Figure 5-31 is an example of this result. H*p in this figure is not bounded and thus, detection by crossing a bound is not possible. At 10 seconds however, a discontinuity in the slope of the curve is plainly visible as well as an increase in the amplitude of oscillation. 89 When noise is factored in, the results of H*p can be very irregular. Depending on the signal to noise ratio of H*p, detection may be noticeable within a single revolution or it may not. A variety of methods and algorithms could be used to detect the changes in H*p. For ideal conditions, simple boundary crossing methods would work fine. However, depending on the state of operations, other approaches might have to be considered. Regardless, sudden changes in H*p within a reasonable noise limit are certainly detectable. Moreover, the less prevalent that the noise is in H*p, the faster imbalances will be detected. 90 CHAPTER VI EXPERIMENTAL APPARATUS AND RESULTS The equations of motion and control algorithms were developed in Chapters 3 and 4. Chapter 5 presents the culmination of these ideas via computer simulations using the Matlab and Simulink. In this chapter, the results described in Chapter 5 will be validated using an experimental rig. This rig consists of two radial magnetic bearings, a one horsepower motor, and a rotor supported by the bearings. The bearings are actuated by a total of eight amplifiers to control four degrees of freedom. Also, six proximity probes, three per bearing, are used to sense the position of the rotor at all times. These are the basic components of the experimental rig. In the following sections, the hardware setup will be described in greater detail, the controller software and code will be discussed, and finally, experimental results will be presented and compared to the findings in Chapter 5. It should be noted that the experimental rig described in this chapter is the same rig used by Matras [15 and 16]. Hence, all information found in the next section is described previously by Matras in earlier works. Figure 6-1 shows a photograph of the active magnetic bearing rig along with the motor and adjustable drive. The lightly shaded cables coming from the circumference of the bearings are connected to the proximity probes and the darker wires towards the back power the bearings. 91 Figure 6-1 ? Active Magnetic Bearing Rig with Motor and Adjustable Driver 6.1 ? Hardware Configuration Originally used as a test model for a magnetic bearing manufacturer, the test rig described here eventually became the property of Auburn University and has since been used for research pertaining to the fields of dynamics and controls. Although little technical information was supplied by the manufacturer, steps were taken by Matras [15] to disassemble the magnetic bearing components to provide measurements and a general idea of its construction. The original drawings by Matras are included in Appendix C of this thesis. Also, basic magnetic properties of the rig as found by Matras can be seen in table 6-1. 92 Variable Description Value N Number of windings on each pole 90 Ag Cross sectional area at the end of each pole 0.00118 m2 g Air gap 0.00025 m M Mass of rotor 5.084 kg Table 6-1 ? Magnetic Bearing and Rotor Properties Each radial bearing stator contains eight poles in a heteropolar configuration which are paired such that each bearing has four pole pairs controlling two degrees of freedom. The axes associated with the two degrees of freedom are perpendicular to each other and rotated 45 degrees from a vertical and horizontal orientation. The stator is constructed using 54 magnetic-annealed iron laminates, each with oxide sheathing on the surface. This coating serves to reduce eddy currents which typically cause energy loss and generate considerable heat in the stator. The sections of the rotor that propagate the magnetic flux generated by the stator are made of the same type of laminates as the stator and are press fit along with solid sections onto a detachable sleeve. Two of these sleeves, one for each bearing, along with a spacer, slide onto a main shaft and are held together by two pieces that screw on to each end of the shaft via tapped sections of the rotor. Finally, a 1 horsepower motor aligned with the bearings is attached to the rotor with a flexible coupling. This motor can be run at speeds of up to 7000 rpm as controlled by an adjustable speed drive. It should be mentioned that having the motor plugged into the same electrical circuit as the bearing components causes considerable interference and audible noise, 93 much like a grinding sound. Originally noted by Matras [15], this problem can be alleviated by having the power to the motor supplied by a completely different electrical circuit than the circuit powering the bearing components. It is also worth mentioning that all metal to metal contact between the motor and the magnetic bearing rig must be eliminated for noise and interference to be at minimal levels. Because of this, the motor is mounted on an acrylic base and attached to the rotor via a rubber coupling. Figure 6-2 ? Close Up of Active Magnetic Bearing Rig A set of proximity probes, used to detect the position of the rotor, are needed for feedback purposes. Each probe uses eddy currents to detect the air gap between the tip of the probe and a solid section of the metal rotor. These signals are then sent to proximitor 94 which in turn outputs a voltage proportional to the detected air gap. The gain used by these probes is 7700 volts per meter. The voltage output by the proximitors is then fed into a signal processing board which amplifies the signal and offsets it such that the maximum resolution prior to analog to digital conversion can be achieved. The analog to digital converter has a maximum voltage range of +10 volts. Since the rotor is only capable of moving 0.25 mm from center position, the voltage range is 2 volts resulting in a voltage bound of -5 to -9 volts. The processing board amplifies this signal by a factor of -3 volts using an analog op-amp circuit and then sums it with another voltage such that when the rotor is in center position, the final output is 0 volts. The added voltage can be adjusted with a potentiometer allowing the output to be calibrated as needed. Since there are six of these proximity probes in all, the signal processing board has six of these particular circuits. The op-amps are powered using a +15 volt power supply. After the signal leaves the processing board, the overall sensor gain becomes Ks = 23610 V/m. Figure 6-3 shows a schematic for the op-amp circuit. Figure 6-3 ? Diagram of Op-Amp Circuit ? + ? + Vin +15 V +15 V -15 V 10 k? 100 k? 10 k? 10 k? 30 k? Vout 95 Since there are four pole pairs per bearing, eight amplifiers are needed to provide the appropriate currents. It was found by Matras [15] that these amplifiers had a gain of 1 amp per volt. It was also documented by Matras that the original bandwidth of the amplifiers was very low. By replacing the adjustment resistor, RH15, with a 3.3 M? resistor and removing the adjustment capacitor, CH17, the bandwidth was significantly raised yielding an approximate value of 1200 Hz per amplifier. Finally, two DC power supplies are required powering four amplifiers each. They were configured as directed by the manufacturer to provide the necessary voltage and encased in aluminum boxes for safety reasons. The amplifiers and signal processing board were then placed in a project box with 8 BNC connectors for control voltage inputs. These signals route through the amplifiers and are output to the bearings via 16 terminals on the project box. The outputs of the proximitors are plugged into the project box by a 25-pin connector, routed through the signal processing board, and are output to the controller via 6 BNC connectors. There is also a switch connected to a +5 volt terminal on the power supply that is output to the controller by a single BNC connector. This switch allows for the bearing to be turned on and off easily without the need of resetting the controller or turning off power supplies. Finally, the system that interfaces the project box to a personal computer is manufactured by dSPACE. This system is a state of the art external processor that provides analog to digital and digital to analog conversion and real-time processing via input and output boards with 32 BNC connectors each. In this particular case, 8 outputs are used for control voltages and 7 inputs are used for the 6 position signals and the 96 switch. The system interfaces with Simulink which allows for a control code to be created, compiled, uploaded to the dSPACE processor, and then ran in real-time. Figure 6-4 ? Project Box for Active Magnetic Bearing Rig 6.2 ? Controller Software and Code The dSPACE software interfaces with a Matlab module called Real-Time Workshop. This module takes a Simulink model, creates a computer code from it, and then compiles it to be used by the dSPACE processor. All values used in the Simulink code can be viewed and changed in real-time through separate software called Control Desk while the control code is on-line and being executed. This software allows the user to create a Graphical User Interface using sliders, knobs, and input textboxes to adjust controller values such as gains and actuation variables. Plotters, gauges, and displays allow the user to monitor signals and other important controller details. Control Desk 97 also provides the ability to save time traces in formats that can be read directly into Matlab for signal processing purposes. The Simulink model used to control the experimental rig is similar in structure to that of the simulations of Chapter 5 with a few significant differences. Obviously the plant is replaced with an output block for digital to analog conversion and an input block for analog to digital conversion. Moreover, the six position signals are sent through a subsystem that manipulates them such that the location of the rotor on each axis at the poles can be determined as four meaningful position signals. Keeping in mind that there are three proximity probes per bearing, two of these probes lie on the top and bottom of the same axis while the third is on the perpendicular axis. Because of this, the bottom reading of the two probes on the same axis is subtracted from the top reading and the result is then divided by two. This allows for an average position reading along that particular axis. The signal fed from the third probe is used directly. Also, because the probes are at a different axial location than the poles, collocation issues are a concern. Figure 6-5 on the next page presents a cutaway of the bearings and rotor showing the respective locations of the sensors and actuators. Equation 6-1 utilizes the values shown in figure 6-5 to adjust the position signals to the positions at the poles. (6-1) ( )31 1 2 11 xxl lxx ?+=? 98 Figure 6-5 ? Diagram of Collocation Issues There are also changes in the state space representation of the estimated plant. Since the position states of the rotor are reported by the proximitors and the amplifier dynamics do not need to be estimated, the velocities of the rotor are the only states that require estimation. Also, it was discovered that including amplifier dynamics in the state estimator caused instability in the system. This is due in large part to the fact that power to the bearings and amplifiers is typically not turned on until after the control code has been loaded and is already running. Hence, amplifier dynamics are neglected in the state space representation of the experimental controller. However, remembering section 4.1, the control voltages are input through the amplifier dynamics. Here the control voltages will be included into the equations of motion for the magnetic bearing and rotor via the current stiffness terms. Instead of these terms being in the form of Ki multiplied by a x1 x3 x1' l1 l2 Poles Proximitor Probes Magnetic Bearing Housings Rotor 99 control current in the state matrix A, Ki will now reside in the input matrix B and will be multiplied by the amplifier gain Ka and the input control voltage up. The estimator will be for the sole purpose of estimating the velocities since the actual positions as opposed to the estimated ones are fed to the control matrix, K. Also, due to the fact that the equations of motion for the state space representation are configured as two translations and two rotations, and the position output for the actual magnetic bearing is based on four translations, it should be noted that the vector of position outputs is multiplied by an invertible transformation matrix, T. This matrix and the format of transformation can be seen below in equations 6-2 and 6-3. (6-2) (6-3) The controller for the experimental rig also incorporates the previously mentioned switch such that when the input from the switch is greater than 1.5 volts, the control voltages output by the controller are multiplied by one. On the other hand, if the switch voltage becomes less than 1.5 volts, the control voltages are multiplied by zero effectively shutting the bearings off. This allows the user to switch the power to the bearings on and off as desired. ? ? ? ? ? ? ? ? ? ? ? ? ? ?= 2 1 2 1 010 010 001 001 L L L L T ( ) ( )??YXYYXX XTX =2121 100 Finally it should be mentioned that the signals read from the analog to digital converter are scaled from a range of +10 volts to +1 volt. Because of this, the inputs are multiplied by 10 as the first step of the control code. Also, all control voltages output from the digital to analog converter are amplified by a factor of 10 requiring the signals to be multiplied by 0.1 prior to being output from the controller. The Simulink model for the experimental rig and the Matlab m-file used to initiate the variables can be found in Appendix B. Manufacturer information regarding all hardware is given in Apendix D. 6.3 ? Experimental Results In this section, plots from the experimental apparatus will be presented as a means of validating the results from Chapter 5. Following the format of Chapter 5, first a time trace of the rotor vibration and a plot of H*p will be shown with parameters similar to those found in section 5.2.1. Next, plots will be shown where the parameter ?H has been varied. Finally, results with frequency discrepancy and additional imbalance frequencies will be presented. 6.3.1 ? Simulated Imbalance Experiments using the previously described experimental rig will be performed where the imbalance is simulated with the controller. Much like the simulations in Chapter 5, since the equations of motion for the system are known as well as the equations that model the imbalance, it is easy to simulate an imbalance with sine and cosine terms. These terms will be multiplied by constants representing mass, rotor speed, and distance of an imbalance mass from center. By simply superimposing these terms 101 with the control outputs, the rotor can be manipulated to follow the same type of path that a natural imbalance would cause. This also allows one to set the exact frequency at which the imbalance will oscillate. By knowing this frequency, the ADR can be accurately programmed to suppress it, allowing for a more controlled experiment. First, the rotor will be excited at 20 Hz with a ?H value of 40,000. The imbalance mass will be simulated as one gram, at an initial distance of one millimeter away from the center of the rotor. When activated, the imbalance mass will be virtually moved away from the center of the rotor an additional millimeter over a tenth of a second for a total distance of two millimeters. The parameters defined in table 6-2 will serve as a baseline for the series of experimental results found in this chapter. Of course, noise will be included into the plots due to the inherent noise in the system. Some of the sources of the noise are sensor runout which will be explained later, geometrical imperfections in the rotor, and non-uniform electrical properties of the rotor. Although steps have been taken to eliminate the noise element, some of these effects are always present. Similar to the simulations shown in Chapter 5 and because of the nature of imbalance, the results from only one axis on one bearing will be studied. This includes all time traces and H*p plots. Also, rotational imbalances will be neglected since the results from the translational imbalances can be easily extended to a rotational case. Again, only the Hp element of the ADR controller will be utilized here. The Gp component is effectively zeroed out and neglected. Figure 6-6 shows a time trace of the effects of a simulated translational imbalance as observed from the X1 axis. 102 0 5 10 15 20 25 30-4 -3 -2 -1 0 1 2 3 4 x 10 -6 Time (s) Tr an sla tio n i n X 1 ( m) Figure 6-6 ? Time Trace for Simulated Imbalance at 20 Hz (Case I) In figure 6-6, the results of a 30 second experiment are presented. Table 6-2 shows the parameters of the experiment and figure 6-7 shows H*p. As can be seen, H*p increases dramatically as soon as the imbalance is detected. By closer inspection (zooming in on the peaks in the time trace), it is seen that although the imbalance increase is actuated at a time of around 16.5 seconds, it is not visible in the time trace until about 16.54 seconds. The dramatic increase in H*p begins some time around 16.52 seconds. This means that the increase in imbalance can be detected through H*p approximately two hundredths of a second before it can be via the time trace. 103 Simulated Running Speed (?) 20 Hz Simulated Initial Imbalance Level 0.0158 N ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 16.4948 s Signal to Noise Ratio for H*p 185.4933 Plant Poles -750 + 250i Table 6-2 ? Parameters for Simulated Imbalance at 20 Hz (Case I) 0 5 10 15 20 25 300.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.065 Time (s) Hp * Figure 6-7 ? H*p for Simulated Imbalance at 20 Hz (Case I) 104 Once again, it can also be seen that since the offset of the simulated imbalance mass doubles (effectively doubling the imbalance force), the steady-state value of H*p is twice the previous value. Also, H*p does indeed settle to an approximately straight line. As discussed in the previous chapter, the noise seen in the time trace is filtered through and shows up in H*p. However, the effects of the noise on H*p do not hinder detection for this magnitude of imbalance. This is certainly expected and supports the results found in section 5.2.1 and 5.2.2. Next, a level of imbalance will be used such that it will not be detectable in the time trace. Instead of superimposing a noise element in the controller, the magnitude of the imbalance will be decreased. This, of course, will mean that more noise will be present in H*p and the signal to noise ratio will be much lower. However, imbalance detection will still be readily observable. In the next set of plots, the imbalance mass and simulated rotor speed will remain the same. However, the offset of the imbalance mass away from the center of the rotor will be decreased to 0.1 millimeters and will then be increased to 0.2 millimeters over a tenth of a second. Figures 6-8 and 6-9 will show the time trace and H*p plot respectively and table 6-3 will show the parameters of the experiment. 105 0 5 10 15 20 25 30-4 -3 -2 -1 0 1 2 3 x 10 -6 Time (s) Tr an sla tio n i n X 1 ( m) Figure 6-8 ? Time Trace for Simulated Imbalance at 20 Hz (Case II) Simulated Running Speed (?) 20 Hz Simulated Initial Imbalance Level 0.0016 N ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 20.1204 s Signal to Noise Ratio for H*p 27.5403 Plant Poles -750 + 250i Table 6-3 ? Parameters for Simulated Imbalance at 20 Hz (Case II) 106 0 5 10 15 20 25 302.5 3 3.5 4 4.5 5 5.5 6 6.5 x 10 -3 Time (s) Hp * Figure 6-9 ? H*p for Simulated Imbalance at 20 Hz (Case II) As expected, although the signal to noise ratio is considerably lower for this case, H*p responds dramatically to the slight increase in imbalance. This series of graphs clearly demonstrates that imbalances are not always readily detectable with a time trace. Noise will certainly be an issue with any electronic system and an active magnetic bearing is no different. However, changes in imbalances can still be detected by examining H*p. In the following sections, the imbalance distance will be set back to the initial value of one millimeter and ?H parameter will be varied. These results can then be compared to figures 6-6 and 6-7 and finally the findings in Chapter 5. 107 6.3.2 ? ?H Variation In Chapter 5, different pole placements were simulated showing the overall effects on H*p. For this experimental rig however, the system response is largely affected by pole placement. Good system poles have been found at -750 + 250i and will not be varied. Instead, only the weighting matrix, ?H, will be changed to show how the response of H*p is affected. Similar to the simulations in Chapter 5, increasing the magnitude of the weighting matrix serves to reduce the settling time of H*p when an increase in imbalance occurs. In the plots shown below, the magnitude of the weighting matrix has been doubled. 0 5 10 15 20 25 30-5 -4 -3 -2 -1 0 1 2 3 4 x 10 -6 Time (s) Tr an sla tio n i n X 1 ( m) Figure 6-10 ? Time Trace for ?H = 80,000 108 Simulated Running Speed (?) 20 Hz Simulated Initial Imbalance Level 0.0158 N ?H 80,000 ADR Frequency 20 Hz Time of Imbalance Increase 14.7500 s Signal to Noise Ratio for H*p 185.4933 Plant Poles -750 + 250i Table 6-4 ? Parameters for ?H = 80,000 0 5 10 15 20 25 300.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.065 Time (s) Hp * Figure 6-11 ? H*p for ?H = 80,000 109 Comparing the time traces, figures 6-6 and 6-10, it can be seen that the increase in imbalance is suppressed quicker with an increased value of ?H. Also by inspection of figures 6-7 and 6-11, it is clear that H*p settles faster with respect to the increase in the weighting matrix. These results are expected and thus support the simulation results found in the previous chapter. Once again, by examining figure 5-16 from Chapter 5, it is obvious that increasing the weighting matrix will dramatically decrease the settling time of H*p until a certain value is reached. However much like an exponential curve, the effects on H*p begin to lessen with increasingly higher values of ?H until no effects can be observed at all. In addition, high gain values tend to exacerbate measurement noise, which is undesirable. 6.3.3 ? Frequency Discrepancies Next, differences between the simulated imbalance frequency and ADR frequency will be explored on the experimental rig. It is expected that H*p will take on a sinusoidal shape as opposed to a straight line and that the steady state frequency of H*p will be twice the simulated running speed or the 2X component. Also, the difference between the simulated imbalance frequency and the ADR frequency is expected to make up part of the transient response of H*p after the increase in imbalance. In the following plots, the simulated running speed will be set to a value of 25 Hz while the ADR frequency will remain at 20 Hz. First, H*p will be shown followed by two plots of an FFT analysis. The first FFT plot shows the first few seconds after the increase in imbalance and the second will show the remaining time. This will illustrate how the frequency content changes during and after the transients of H*p have died down. 110 0 5 10 15 20 25 301 1.5 2 2.5 3 3.5 4 4.5 5 x 10 -3 Time (s) Hp * Figure 6-12 ? H*p for Frequency Discrepancy Simulated Running Speed (?) 25 Hz Simulated Initial Imbalance Level 0.0158 N ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 11.9412 s Signal to Noise Ratio for H*p 185.4933 Plant Poles -750 + 250i Table 6-5 ? Parameters for Frequency Discrepancy 111 0 10 20 30 40 50 60 70 80 90 100 110 1200 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Frequency (Hz) Am pli tud e Figure 6-13 ? FFT of H*p between 12 - 15 Seconds for Frequency Discrepancy Since the increase in imbalance occurs just before 12 seconds, an FFT is taken from 12 to 15 seconds. In figure 6-13 above, it can be seen that there is a major frequency component at 5 Hz and a second major frequency component at 50 Hz. The first peak accounts for the 5 Hz difference in simulated running speed and ADR frequency. The second peak is twice the simulated running speed. Looking next at figure 6-14, it can be observed that the 5 Hz frequency component has died out leaving only the 2X and 4X components. These results certainly support the findings in section 5.2.5 and further prove the validity of the simulations. In the next section, additional imbalance frequencies and their effects on H*p will be examined. 112 0 20 40 60 80 100 1200 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Frequency (Hz) Am pli tud e Figure 6-14 ? FFT of H*p between 15 - 30 Seconds for Frequency Discrepancy 6.3.4 ? Additional Imbalance by Spinning the Rotor In this section, the primary imbalance will be simulated as before, but with the additional imbalance of the rotor as it is spun up to a different speed. First however, a phenomenon called sensor runout must first be introduced. First described by Setiawan et al. [22] and again by Matras [15], this phenomenon is caused by the non uniform electrical properties of the rotor. When the shaft is rotated, this sensor runout occurs at frequencies which are integer values of the running speed. This means that if an FFT analysis of the rotor vibration was performed for the rotor spinning at a particular speed, there would be several peaks at not only the running speed, but also at integer multiples 113 of the rotor speed. To prove this here, the motor will be spun up to a speed of 20 Hz with the ADR enabled to reject at 20 Hz. An FFT analysis will be performed for the time trace. This is shown in figure 6-15 below. 0 20 40 60 80 100 120 140 160 180 2000 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Frequency (Hz) Am pli tud e ( m) Figure 6-15 ? FFT of Time Trace for Rotor Spinning at 20 Hz As can be seen in figure 6-15, the peak at 20 Hz is mostly suppressed by the ADR while frequency components at integer values of the running speed remain high. Because of this, H*p will certainly react to all of these frequencies. Of course, the extra frequency components can be reduced by simply adding in more ADR frequencies to reject. However, as discussed in section 5.2.6, adding in more ADR frequencies can have adverse effects on H*p. This will be investigated experimentally later. In the next plot, 114 an FFT analysis of H*p is shown with the ADR frequency set to 20 Hz and a weighting matrix value of 40,000. 0 20 40 60 80 100 120 140 160 180 2000 0.5 1 1.5 2 2.5 3 3.5 4 Frequency (Hz) Am pli tud e Figure 6-16 ? FFT of H*p for Rotor Spinning at 20 Hz By examination of figure 6-16, one can observe peaks at the running speed and integer values of the running speed. This can be explained by remembering section 5.2.6 where extra excitation frequencies were simulated. Hence, for every frequency component not equal to 20 Hz, the frequency difference between that frequency and the rejected frequency will show up in H*p. Moreover, since there is a major peak in the time trace at 40 Hz, a major peak in H*p at 20 Hz will occur which is the difference between the 40 Hz frequency and the rejected frequency of 20 Hz. The other peaks can be 115 explained with similar reasoning. As for the secondary peaks (also explained in section 5.2.6), they are part of the peaks present at integer multiples of the running speed. For example, considering the 40 Hz peak in the time trace, it can be observed that a secondary peak in the FFT analysis of H*p would lie at 60 Hz. This would be the 20 Hz difference subtracted from the 2X component of 40 Hz. For the next set of plots, the rotor speed will again be set to 20 Hz. However, a simulated imbalance at 25 Hz will be added and the ADR will be set to reject at 25 Hz. 0 5 10 15 20 25 300.04 0.05 0.06 0.07 0.08 0.09 0.1 Time (s) Hp * Figure 6-17 ? H*p for Simulated Imbalance and Spinning Rotor 116 Actual Running Speed ~20 Hz Simulated Imbalance Speed 25 Hz Simulated Initial Imbalance Level 0.0158 N ?H 40,000 ADR Frequency 25 Hz Time of Imbalance Increase 11.3486 s Signal to Noise Ratio for H*p 185.4933 Plant Poles -750 + 250i Table 6-6 ? Parameters for Simulated Imbalance and Spinning Rotor 0 10 20 30 40 50 60 70 80 90 1000 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 Frequency (Hz) Am pli tud e ( m) Figure 6-18 ? FFT of Time Trace for Simulated Imbalance and Spinning Rotor 117 0 10 20 30 40 50 60 70 80 90 1000 2 4 6 8 10 12 14 Frequency (Hz) Am pli tud e Figure 6-19 ? FFT of H*p for Simulated Imbalance and Spinning Rotor Examining figures 6-18 and 6-19, the cause of the extra frequency components can be explained. In figure 6-18, once again frequency components of the actual running speed and its multiples are readily seen. Also, there is a peak at 50 Hz which is the 2X component of the simulated imbalance. Inspection of figure 6-19 shows that there is a major peak at 5 Hz, the frequency difference between the rejected frequency and actual running speed. The next peak is at 15 Hz, which is the difference between 40 Hz component and rejected frequency. The 35 Hz peak is due to the 60 Hz component of the time trace. The peak at 45 Hz however, could be explained two different ways. The first explanation could be that the 5 Hz difference is added to the 2X component of the 20 Hz 118 frequency from the time trace and the second explanation could be the discrepancy between the 50 Hz component and the rejected frequency. The rest of the peaks can be established similarly. Regardless, figure 6-17 shows how H*p responds dramatically to the increase in imbalance. However, as shown in Chapter 5, when two imbalance frequencies are very close together, the difference between the secondary frequency and the rejected frequency will dominate H*p. This can be shown by setting the motor speed, simulated imbalance, and ADR frequencies to 20 Hz. Because the rotor speed is not exactly 20 Hz, the difference will appear in H*p. Figure 6-20 shows this effect. 0 5 10 15 20 25 300.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 Time (s) Hp * Figure 6-20 ? H*p for Simulated Imbalance and Rotor Speed at 20 Hz 119 Actual Running Speed ~20 Hz Simulated Imbalance Speed 20 Hz Simulated Initial Imbalance Level 0.0158 N ?H 40,000 ADR Frequency 20 Hz Time of Imbalance Increase 13.3213 s Signal to Noise Ratio for H*p 185.4933 Plant Poles -750 + 250i Table 6-7 ? Parameters for Simulated Imbalance and Rotor Speed at 20 Hz Figure 6-20 shows how the motor speed is very close to 20 Hz, but is not 20 Hz exactly. This induces a dominating low frequency, high magnitude sinusoidal behavior in H*p. The imbalance detection in this case must be performed by observing the discontinuity in the slope of H*p. Although the change in imbalance shows up in a different form, it is still readily observable. As previously explained, the extra frequency content is due to sensor runout and the frequency characteristics of H*p. It should be noted that the times of the increases in imbalance are known exactly in all the previous plots. This is due to the fact that the simulated imbalance is actuated with the controller via the program Control Desk. Although all plots using the experimental rig are performed for a rotor speed of 20 Hz, the presented results are meant to support the findings in Chapter 5. It is expected 120 that higher frequencies would yield similar results. In addition, the frequencies and rotor speed are kept low for safety reasons. 6.4 ? AFRL Rotor Test Rig For the next set of experiments, a different experimental rig was used. This rig was located and built on site at the Air Force Research Laboratory at Kirtland Air Force Base in Albuquerque, New Mexico for research in the FACETS program [6]. The FACETS (Flywheel Attitude Control-Energy Storage System) program, as described in chapter 2, used this rig as a means to test the adaptive control and other technologies on the ground first. Built and configured similarly to the one described in sections 6.1 and 6.2, the major difference with this active magnetic bearing rig was that the air gap is much smaller. This was a design choice such that it would yield greater stiffness with less current making it more efficient. Another difference was that the rotor was powered with an air turbine as opposed to an electric motor. This again was a design choice that allowed the rotor to be spun to higher speeds and decreased the electrical noise introduced to the system. One of the down sides to using an air turbine was the fact that it was powered with compressed air as supplied by the building. Hence, the speed of the rotor could not be fully controlled since the air pressure was varied using a simple mechanical valve. Moreover, any time another lab used compressed air, the supply to the air turbine was compromised. Much like an electric motor, the speed of the rotor was never exactly known. The controller was designed and implemented by Alex Matras [17] using PID methods and the aforementioned extension of the ADR control algorithm as well as 121 several error protocols. The system characteristics and controller design are explained in great detail in Matras [17] and will not be included here since this author's contribution to the controller on that particular rig was minimal. Figure 6-21 ? AFRL Magnetic Bearing Test Rig 6.4.1 ? Test Wheel The original intent of the AFRL rig for this research was to have a test wheel designed and implemented on the rotor to physically simulate an increase in imbalance. The test wheel was designed such that a small metal ball would sit in a cavity enclosed within the wheel. At a chosen time, strong actuation magnets would be moved close to the side of the test wheel to magnetically attract the ball toward a chute that extended 122 radially out from the center of the wheel. When the ball was moved in line with the enclosed chute, centrifugal forces caused by the spinning rotor would then throw the ball radially out and down the chute until it reached the outermost bound. This in turn would cause an increase in mass imbalance by changing the distance of the mass from the center of the wheel. The distance the ball traveled could be manipulated by simply screwing in a threaded end piece that served as the outer most bound for the chute. Also, the design incorporated several threaded holes around the circumference of the test wheel with the intent of balancing the wheel. Figure 6-22 presents a cutaway model of the test wheel which shows how the ball is actuated and designed to travel. Figure 6-22 ? Cutaway Model of Test Wheel Actuation Magnets Test Wheel Steel Imbalance Ball Chute Cavity Direction of Travel Threaded End Pieces 123 The test wheel was originally made out of aluminum because it was thought that the material possessed no magnetic qualities. However, it was discovered that when the aluminum wheel was spun up, it actually caused a repelling force when approached with a magnet. This was found to be due to eddy currents formed in the spinning aluminum by the magnet. Hence, the wheel was redesigned and built using a plastic material called Delrin for its high strength, low weight, and complete lack of magnetic properties. This alleviated the problem completely. A drawing for the construction of the test wheel can be found in Appendix C of this thesis. Figure 6-23 ? Photograph of Delrin Test Wheel For testing, two steel imbalance balls were used with each being placed in perpendicular chutes. The steel balls were actuated using two wishbone shaped brackets that were attached to an ACME screw. When the balls were in their initial positions in 124 the cavities, one bracket holding four actuation magnets was positioned so to attract the balls toward the flat side of the cavities. When the balls were to be magnetically moved, the ACME screw was turned by hand which moved the first bracket away from the side of the wheel and positioned the second bracket on the opposite side of the wheel closer. The second bracket was designed to hold four actuation magnets with opposite poles facing the wheel. The reasoning behind this was that while the first bracket was close to the wheel, the balls would be gradually magnetized. When the second bracket was moved close to the wheel, the opposite polarity would help cause the balls to be attracted toward the chutes. Figure 6-24 shows the initial wishbone bracket actuator. Figure 6-24 ? Initial Wishbone Bracket Actuator Design 125 It was found after some initial tests that the second bracket was not quite magnetically strong enough to pull the balls towards the chutes. The second bracket was then redesigned and built such that it would hold eight magnets. The magnets used in these brackets were chosen to be neodymium magnets because of their high magnetic properties and relatively small size. Figure 6-25 shows a photograph of the AFRL test rig with the test wheel and wishbone bracket actuator in place. In this photograph, the rotor is spinning. Figure 6-25 ? ARFL Magnetic Bearing Test Rig with Test Wheel and Wishbone Bracket Actuator After several tests, it was found that the test wheel could only operate correctly in a range of speeds between approximately 40 to 60 Hz. At slower speeds, the centrifugal 126 forces were not great enough to throw the balls down the chutes because of the strength of the magnets. At higher speeds, the magnets were unable to pull the balls to the chutes because the centrifugal forces caused increased rolling resistance between the steel balls and the wheel. Another observation was that the speed of the rotor changed when the second bracket approached the wheel due to the magnetic interaction between the magnets and the steel balls. Regardless, the test wheel was considered successful and useful with in this range of speed. 6.4.2 ? Experimental Results The experiments performed using this test rig along with the test wheel predates many of the simulations shown in Chapter 5. Also, since the rig was located in Albuquerque, New Mexico, it had limited availability to the Auburn University research team. The speed of the rotor was measured using a tachometer which can be seen at the end of the rotor in figure 6-25. Although the tachometer gave an approximate speed, it was not accurate enough to determine an adequate ADR frequency. Due to these factors, several ADR frequencies were added to the controller hoping to at least cover a range of speed at which the rotor was spinning. Different values were used for different experiments. Some experiments used frequencies that covered a range at the approximated running speed and others experiments used frequencies that covered ranges at the running speed and its 2X component. At the time of these experiments it was unknown how H*p reacted to additional ADR frequencies. Although the imbalance is only slightly detectable in the given plots by noticing the discontinuity in the slope of H*p, the results do indeed support the simulations in the later part of section 5.2.6. 127 It should be noted that since the controller was designed and implemented by someone other than this author, the values for the weighting matrices, ?G and ?H, will not correspond with all the previous experiments. Also, because of the nature of the PID controller, the poles were not placed in a state space sense and PID gain values were not recorded. Moreover, the time of the increase in imbalance will not be known exactly since there was no way to document this with the controller other than visually detecting it with H*p. The plots below are for an approximate running speed of 60 Hz and the H*p data is taken at the ADR frequency of 60 Hz. 0 2 4 6 8 10 12 14 16 18 20-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 Time (s) Tr an sla tio n i n X 1 ( m) Figure 6-26 ? Time Trace for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz 128 Actual Running Speed ~60 Hz ?G 50 Gp Saturation 8 ?H 30 ADR Frequencies 59.7 Hz 59.8 Hz 59.9 Hz 60.0 Hz 60.1 Hz 60.2 Hz Table 6-8 ? Parameters for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz 0 2 4 6 8 10 12 14 16 18 200 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Time (s) Hp * Figure 6-27 ? H*p for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz 129 Examining figures 6-26 and 6-27, it can be seen that the time trace gives no indicator of the time when the balls were actuated. However, the H*p plot shows the characteristic discontinuity in the slope of the curve between 17 and 18 seconds. It is also observed that the dominant sinusoidal frequency of H*p becomes smaller as the time progresses. By examining a plot of the rotor speed via the tachometer, it can be seen that the speed does indeed increase causing the frequency difference between the 60 Hz ADR frequency and the running speed to decrease. Figure 6-28 shows a plot of the tachometer signal. Once again, this data corresponds very well with the simulations found in Chapter 5. 0 2 4 6 8 10 12 14 16 18 2058.6 58.8 59 59.2 59.4 59.6 59.8 60 60.2 60.4 60.6 Time (s) Ro tor S pe ed (H z) Figure 6-28 ? Time Trace of Rotor Speed for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz 130 The tachometer plot doesn't specifically prove that the rotor speed changes because of the magnetic interaction between the balls and the magnets, but it does give evidence that supports it. Since any and all pressure changes in the compressed air supply to air turbine affects its speed, this could certainly be a cause of the change. However, it is probably not simply coincidental that the rotor sped up as the bracket with eight magnets was driven closer to the test wheel. 0 1 2 3 4 5 6 7 8 9 100 50 100 150 200 250 300 350 400 450 500 Frequency (Hz) Am pli tud e Figure 6-29 ? FFT of H*p between 0 - 8 Seconds for AFRL Magnetic Bearing Test Rig with Test Wheel at ~60 Hz An FFT analysis was performed for H*p for the time between 0 and 8 seconds since that was the span of time that the rotor speed stayed somewhat constant. The plot of the FFT results was scaled so that the major peak at about 0.7 Hz could be seen. This 131 peak does indeed account for the difference between the average speed of 59.3 Hz and the 60 Hz rejected frequency. Unfortunately, the only usable data with the test wheel is the above experiment. Three distinct experiments were run using the test wheel at speeds of 40, 50, and 60 Hz. The balls did not actuate properly with the 40 Hz experiment, and the 50 Hz experiment did not show any characteristics of increases in imbalance by examining H*p. Other experiments were ran using controller simulated imbalances but will not be included here as the results were similar to those from the other rig. As previously mentioned, the affects of multiple ADR frequencies had not been thoroughly investigated at that point and the disadvantages were unknown. Also, the time available for performing the experiments was very limited and there was no chance to examine the results in detail and perform any follow up experimentation. 6.5 ? Response Time The time it takes for H*p to respond to an increase in imbalance can be found using a variety of methods. For the sake of simplicity, the times reported here will be found by simply detecting when H*p exceeds some bound. This bound will be set as the maximum value of H*p before the simulated imbalance is activated. Times will not be found for H*p in figures 6-20 and 6-27 due to their sinusoidal nature. In Chapter 1 of this thesis, it was proposed that the imbalance could possibly be detected within a single revolution of the rotor. Table 6-9 on the following page shows the time of the actuated imbalance, the time H*p crosses the set bound, and the difference between them. The 132 table will also show the simulated running speed, the time it takes for one revolution at that speed, and the percentage of one revolution it takes for the imbalance to be detected. ? pH Curves Simulated Running Speed Time of One Revolution Time of Imbalance Increase Time of Imbalance Detection Detection Time Percentage of One Revolution Figure 6-7 20 Hz 0.05 s 16.4948 s 16.5154 s 0.0206 s 41 % Figure 6-9 20 Hz 0.05 s 20.1204 s 20.2078 s 0.0874 s 175 % Figure 6-11 20 Hz 0.05 s 14.7500 s 14.7783 s 0.0283 s 57 % Figure 6-12 25 Hz 0.04 s 11.9412 s 11.9793 s 0.0381 s 95 % Figure 6-17 25 Hz 0.04 s 11.3486 s 11.4430 s 0.0944 s 236 % Table 6-9 ? Time Responses of Experimental H*p Plots Examination of the table above shows that for high signal to noise ratios, H*p can indeed respond in less than one revolution. Even with decreased signal to noise ratios or additional excitations with magnitudes much higher than that of the increasing simulated imbalance, the response time of H*p is still within 3 revolutions. These results are of course for low speeds, but the results can certainly be extended to higher speeds. Better methods for detecting the increase in imbalance with H*p would surely decrease the detection time. However, it is beyond the scope of this thesis to thoroughly examine such methods. The methods used here are meant to provide only a basic idea of detection time. 133 CHAPTER VII CONCLUSIONS AND RECOMMENDATIONS 7.1 ? Conclusions There are different methods of attempting to detect imbalance increases associated with flywheels on active magnetic bearings. One such method uses algorithms that detect boundary crossing in time traces. Although this method may be effective to some degree, the results discussed in both Chapters 5 and 6 shows that small imbalance changes can not always be detected through the time trace, especially with noisy plots. However, it has been presented that the technique known as adaptive disturbance rejection can have many positive effects on vibration suppression. Not only can the technique effectively suppress vibrations at specified frequencies as presented in Matras et al. [15, 16, and 17], but increases in the magnitude of vibration, such as an imbalance change, can be detected easier and more efficiently with H*p than by examining time traces of rotor vibration. The method of detecting imbalance changes with H*p has also been shown to be fairly robust, especially when specifying only one ADR frequency to reject. In many cases, the rotor speed may not be known exactly. This does impede the vibration suppression element of the ADR, but it does not change the ability to detect changes in the magnitude of the imbalance. Many scenarios that can affect the response of H*p have 134 been researched. It was found that additional imbalance frequencies and additional ADR frequencies have the most adverse effect on H*p, especially when the additional frequencies are very close to the specified ADR frequency and/or primary imbalance frequency. However, it was also shown that with these types of scenarios, the change in imbalance was still detectable to some degree. It was concluded that the ADR was most effective in suppressing vibrations and responding to imbalance changes when only one frequency was specified to be rejected. Moreover, that frequency should be the exact speed of the rotor. When this is the case, not only can the ADR give the maximum vibration suppression, but H*p can react to perturbations in the imbalance level with the utmost efficiency. It should be mentioned that the speed at which H*p can respond to changes in imbalance levels with high signal to noise ratios and low speeds was found to be within one revolution of the rotor. Although noise does indeed filter into H*p and affect the detection time, imbalance changes with low signal to noise ratios were still detected within a few revolutions. It is certainly possible that even at very high rotor speeds, detection time via H*p can be short enough to enact safety procedures that can shut down a flywheel system safely. Overall, this research effort was successful. The experimental results found in Chapter 6 match the simulations presented in Chapter 5 very closely. This proves the validity of the linearized model as derived in Chapter 3. This also shows how the controller presented in Chapter 4 effectively floats the bearing, ensuring stability and allowing for changes in simulated or actual rotor imbalance to be detected quickly and efficiently. 135 7.2 ? Recommendations Although much time and effort has gone into this research, there are aspects and characteristics of the adaptive imbalance detection gain that need further research consideration. The systems of equations shown in chapter 5, equations 5-2 through 5-4, 5-5 through 5-7, and 5-8 through 5-10, should be examined and closed form solutions found. Not only would this analytically show the resulting frequencies in H*p, but would also allow one to analytically prove stability of the ADR for those particular scenarios. It would also be useful and insightful to conduct more experimentation with the test wheel described in Chapter 6. Due to a lack of time, limited experimentation with the AFRL facility was done. It is suggested that since the test wheel gives a more physical representation of changes in mass imbalance, future experiments with the flywheel should include changing the parameter, ?H, and examining the effects of multiple ADR frequencies further as well as using only a single ADR frequency. Other future work could also include researching the extensions made to the ADR technique as found in Matras [17] and determining if and how they might affect the outcome of H*p. All simulations and experimentations with exception of the AFRL rig neglect the Gp component of the ADR. The controller used on the AFRL rig utilized the adaptive gain, Gp. Finally, research should be done in the area of developing better algorithms to detect discontinuities in the slope of H*p. Results have shown that H*p can take on fairly straight lines or high amplitude, low frequency sinusoidal paths. Developing better algorithms would not only allow for the characteristics of a change in imbalance to be tracked easier, but also decrease detection time even more. 136 BIBLIOGRAPHY 1. Ahmed, J., Bernstein, D. S., 1999, "Adaptive Force Balancing of an Unbalanced Rotor," Proceedings of the 38th IEEE Conference on Decision and Control, pp. 773-778, Phoenix, Arizona. 2. Amati, N., Brusa, E., 2001, "Vibration Condition Monitoring of Rotors on AMB Fed by Induction Motors," IEEE/ASME International Conference on Advanced Intelligent Mechatronics Proc., pp. 750-756, Como, Italy. 3. Brogan, W. L., 1991, Modern Control Theory, 3rd Edition, New Jersey, Prentice Hall. 4. Coombs, T. A., Cardwell, D. A., Campbell, A. M., 1997, "Dynamic Properties of Superconducting Magnetic Bearings," IEEE Transactions on Applied Superconductivity, Vol. 7, No. 2, pp. 924-927. 5. Earnshaw, S., 1842, "On the Nature of Molecular Forces Which Regulate the Constitution of the Luminiferous Ether," Transactions of the Cambridge Philosophical Society, Vol. 7, pp. 97-112. 6. Fausz, J. L., Richie, D. J., 2000, "Flywheel Simultaneous Attitude Control and Energy Storage Using a VSCMG Configuration," Proceedings of the 2000 IEEE International Conference on Control Applications, pp. 991-995, Anchorage, Alaska. 7. Flowers, G. T., Sz?sz, G., Trent, V.S., Greene, M. E., 2001 "A Study of Integrally Augmented State Feedback Control for an Active Magnetic Bearing Supported Rotor System," Journal of Engineering for Gas Turbines and Power, Vol. 123, pp. 377-382. 8. Fuentes, R.J., Balas, M. J., 2000, "Direct Adaptive Rejection of Persistent Disturbances," Journal of Mathematical Analysis and Applications, Vol. 251, pp. 28-39. 9. Ginsberg, J. H., 1998, Advanced Engineering Dynamics, 2nd Edition, New York, Cambridge University Press. 137 10. Hai, X., 2002, "Dynamic Analysis of High-Speed Multi-Direction Composite Flywheel Rotors," M. S. Thesis, Auburn University, Auburn, Alabama. 11. Hai, X., Flowers, G., Fausz, J., Horvath, R., 2005, "Vibration Analysis and Health Monitoring of Cracks in Composite Disk Rotor Systems," Proceedings of ASME 2005 International Design Engineering Technical Conference & Computers and Information in Engineering Conference, pp. 1-8, Long Beach, California. 12. Ichihara, T., Matsunaga, K., Kita, M., Hirabayashi, I., Isono, M., Hirose, M., Yoshii, K., Kurihara, K., Saito, O., Saito, S., Murakami, M., Takabayashi, H., Natsumeda, M., Koshizuka, N., 2005, "Application of Superconducting Magnetic Bearings to a 10 kWh-Class Flywheel Energy Storage System," IEEE Transactions on Applied Superconductivity, Vol. 15, No. 2, pp. 2245-2248. 13. Knospe, C. R., Hope, R. W., Fedigan, S. J., Williams, R. D., 1995, "Experiments in the Control of Unbalance Response Using Magnetic Bearings," Mechatronics, Vol. 5, No. 4, pp. 385-400. 14. Maslen, E., 2000, Magnetic Bearings, University of Virginia, Charlottesville, Virginia, http://www.people.virginia.edu/~ehm7s/MagneticBearings/mag_brgs.pdf. 15. Matras, A. L., 2003, "Implementation of Adaptive Disturbance Rejection Control in an Active Magnetic Bearing System," M. S. Thesis, Auburn University, Auburn, Alabama. 16. Matras, A. L., Flowers, G. T., Fuentes, R., Balas, M., Fausz, J., 2003, "Suppression of Persistent Rotor Vibrations Using Adaptive Techniques," Proceedings of the 2003 ASME International Mechanical Engineering Conference and Exposition, Vol. 2, Washington, DC. 17. Matras, A. L., 2005, "Applied Adaptive Disturbance Rejection Using Output Redefinition on Magnetic Bearings," Ph. D. Thesis, University of Colorado, Boulder, Colorado. 18. Matsumura, F., Kobayashi, H., 1981, "Fundamental Equation for Horizontal Shaft Magnetic Bearing and Its Control System Design," Electrical Engineering in Japan, Vol. 101, No. 3, pp. 123-130. 19. Mohamed, A. M., Emad, F. P., 1993, "Nonlinear Oscillations in Magnetic Bearing Systems," IEEE Transactions on Automatic Control, Vol. 38, No. 8, pp. 1242- 1245. 138 20. "The Multiple 'Attractions' of a Magnetic Bearing," Processing Magazine, 1/10/2006, http://www.processingmagazine.com/productarticle.asp?ArticleID=827. 21. Roithmayr, C. M., 1999, "International Space Station Attitude Control and Energy Storage Experiment: Effects of Flywheel Torque," NASA Technical Memorandum 209100. 22. Setiawan, J. D., Mukherjee, R., and Maslen, E. H., 2002, "Synchronous Sensor Runout and Unbalance Compensation in Active Magnetic Bearings Using Bias Current Excitation", ASME Journal of Dynamic Systems, Measurement and Control, Vol. 124, pp. 14-24. 23. Shafai, B., Beale, S., LaRocca, P., Cusson, E., 1994, "Magnetic Bearing Control Systems and Adaptive Forced Balancing," IEEE Control Systems, Vol. 14, No. 2, pp. 4-13. 24. Shen, J. X., Vilathgamuwa, D. M., Tseng, K. J., Chan, W. K., 1999, "A Novel Compact PMSM with Magnetic Bearing for Artificial Heart Application," Proceedings of the 1999 IEEE Industry Applications Conference - 34th IAS Annual Meeting, Vol. 2, pp. 1201-1207, Phoenix, Arizona. 25. Shiue, F. W., Lesieutre, G. A., Bakis, C. E., 2001, "A Virtual Containment Strategy for Filament-Wound Composite Flywheel Rotors with Damage Growth," Journal of Composite Materials, Vol. 36, No. 9, pp. 1103-1120. 26. Shiue, F. W., Lesieutre, G. A., Bakis, C. E., 2002, "Condition Monitoring of Filament-Wound Composite Flywheels Having Circumferential Cracks," Journal of Spacecraft and Rockets, Vol. 39, No. 2, pp. 306-313. 27. Tessier, L., 1998, "Magnetic Bearing Systems - Tutorial," ASME/IGTI 1998 TurboExpo, Tutorial Handout. 28. Vance, J. M., 1988, Rotordynamics of Turbomachinery, New York, Wiley- Interscience. 29. Wilson, B. C. D., 2004, "Control Designs for Low-Loss Active Magnetic Bearings: Theory and Implementation," Ph. D. Thesis, Georgia Institute of Technology, Atlanta, Georgia. 30. Zhu, C., Robb, D. A., Ewins, D. J., 2003, "The Dynamics of a Cracked Rotor with an Active Magnetic Bearing," Journal of Sound and Vibration, Vol. 265, pp. 469- 487. 139 APPENDICES 140 APPENDIX A COMPUTER SIMULATION CODE Appendix A contains the block diagrams created in the program Simulink used to produce the simulations described in Chapter 5. Before each simulation, the model runs an m-file called state_space_setup_LRSR_magbear_real which defines all of the parameters and variables used in the model. This m-file is included as well as another m- file called imbalance_fft_analysis which was used to run FFT analyses on the various time traces and H*p plots. Descriptions of each of the blocks can be found in the back of this appendix. The blocks in the Simulink model with dropped shadows are subsystems which are shown on subsequent pages. All parameter and gain blocks are labeled and their values can be found in the m-file state_space_setup_LRSR_magbear_real. The outputs of the model go to scopes, which can view the outputs during the simulations, and output blocks. These output blocks allow for the data to be written to files and processed by Matlab. The lines connecting the various blocks have numbers associated with them. These numbers indicate the size of the vectors running to and from each block. Filename: imbalance_plant_real.mdl Simulation Parameters Start Time: 0.0 Stop Time: 20.0 Solver Type: Fixed-Step, ode4 (Runge-Kutta) Fixed Step Size: 0.0001 Model Callbacks Model Initialization Function: state_space_setup_LRSR_magbear_real Simulation Start Function: state_space_setup_LRSR_magbear_real 14 1 Fig ur e A -1 ? S im ula tio n M od el Bl oc k D iag ra m For the State-Space(Real Plant) block A = Areal B = [Breal Gamma] C = Creal D = Dreal Initial Conditions = initial Y X K*u Transition Matrix, T time To Workspace1 simout To Workspace x' = Ax+Bu y = Cx+Du State-Space(Real Plant) Uc In Uc Out Integral Control Out1 Imbalance-Gravity Generator Y U Uc Controller-ADR(Model Plant) Clock Uc In Uc Out CV Splitter Beta Band-Limited White Noise Alpha K*u 1/Ks[5x1] [5x1] [4x1] 4 [4x1] [13x1] [4x1] [4x1] [4x1] [4x1] [4x1] [4x1] [8x1][4x1] [4x1] 4 [4x1] imbalance_plant_real 14 2 Fig ur e A -2 ? S im ula tio n C on tro lle r B loc k D iag ra m 1 Uc K*u Transition Matrix, Inv(T) Y U x(hat) Estimator K*u Control Matrix, K Y ADR Out1 ADR 2 U 1 Y 4 4 4 4 4 12 [4x1] 4 4 [4x1] [4x1] imbalance_plant_real/Controller-ADR(Model Plant) 143 Figure A-3 ? Simulation Estimator Block Diagram z'(hat) z(hat) 1 x(hat) K*u L 1 s Integrator K*u B K*u A-LC 2 U 1 Y 4 12 12 12 1212 12 12 4 12 12 imbalance_plant_real/Controller-ADR(Model Plant)/Estimator 144 Figure A-4 ? Simulation ADR Block Diagram 1 ADR Out1 Phi Transpose Phi Phi Y Phi Transpose Phi Hp*Phi Hp Subsystem1 Y Gp*y Gp Subsystem1 1 Y 4 4 4 4 4 [2x1] [1x2] [4x1] [4x1] [4x1] imbalance_plant_real/Controller-ADR(Model Plant)/ADR 14 5 Fig ur e A -5 ? S im ula tio n G p S ub sys tem B loc k D iag ra m 1 Gp*y Product9 Product8 Product7 Product6 Product4 Product3 Product1 Product 1 s Integrator3 1 s Integrator2 1 s Integrator1 1 s Integrator -1 Gain5 -1 Gain4 -1 Gain2 -1 Gain1 0 Delta G(y) 0 Delta G(x) 0 Delta G(b) 0 Delta G(a) K*u 1/Ks 1 Y 4 4 4 imbalance_plant_real/Controller-ADR(Model Plant)/ADR/Gp Subsystem1 146 Figure A-6 ? Simulation Disturbance Vector Block Diagram 2 Phi 1 Phi Transpose Sine Wave uT Math Function1 uT Math Function Cos Wave [1x2] [1x2] [1x2]2 [2x1] imbalance_plant_real/Controller-ADR(Model Plant)/ADR/Phi 14 7 Fig ur e A -7 ? S im ula tio n H p S ub sys tem B loc k D iag ra m HpHp' 1 Hp*Phi Matrix Multiply Product1 Matrix Multiply Product 1 s Integrator In1 Hp Veiwer -1 Gain1 -C- Delta H -K- 1/Ks 3 Phi 2 Phi Transpose 1 Y 4 [4x1] 4 [2x2] [4x2] [4x2][1x2] [2x1] [4x2] [4x2] [4x2] imbalance_plant_real/Controller-ADR(Model Plant)/ADR/Hp Subsystem1 14 8 Fig ur e A -8 ? S im ula tio n H *p Ge ne ra tor B loc k D iag ra m adaptive curve Adaptive_Curve To Workspace1 Hp_Sin To Workspace Terminator Matrix Multiply Product4 uT Math Function3 u2 Math Function2 u2 Math Function1 sqrt Math Function Hp Sin X Hp Cos X-C- Constant1 Add 1 In1 [1x2] [1x2] [1x2] [4x2] 4 [1x4] [1x4] imbalance_plant_real/Controller-ADR(Model Plant)/ADR/Hp Subsystem1/Hp Veiwer 14 9 Fig ur e A -9 ? S im ula tio n I mb ala nc e a nd G ra vit y G en era tor B loc k D iag ra m 1 Out1 -C- omega d distance cos Trigonometric Function1 sin Trigonometric Function Sine Wave theta SaturationRamp Product4 Product2 Product1 u2 Math Function2 uT Math Function1 uT Math Function -C- Gravity (M+m)*G Cos Wave theta Clock [1x5] [5x1] 2 2 5 2 imbalance_plant_real/Imbalance-Gravity Generator 150 Figure A-10 ? Simulation Integral Controller Block Diagram 1 Uc Out Kg Kg 1 s Integrator 1 Uc In [4x1] [4x1] [4x1] [4x1] [4x1][4x1] [4x1] imbalance_plant_real/Integral Control 151 Figure A-11 ? Simulation Control Voltage Splitter Block Diagram 1 Uc Out vb8 vb8 vb7 vb7 vb6 vb6 vb5 vb5 vb4 vb4 vb3 vb3 vb2 vb2 vb1 vb1 uT Math Function1 uT Math Function 1 Uc In [4x1] 8 [1x8] [8x1] imbalance_plant_real/CV Splitter 152 Description of Simulink System and Subsystem Blocks ? imbalance_plant_real o This block diagram creates a state space with the vectors as defined in the m-file state_space_setup_LRSR_magbear_real. o The outputs are branched off. One branch sends the outputs through a block that divides the readings by the sensor gain so the four output vectors can be read in meters. o The second branch sends the outputs through a transformation matrix that transforms the 2-translational, 2-rotational coordinate system into a 4- translational coordinate system which are then sent to the control block. o The control voltages output by the controller are split so that one branch sends the control voltages back to the estimator and the other branch proceeds to the integral control block. o The output of the control voltage splitter is muxed with the imbalance/gravity vectors and sent back to the state space block as inputs ? imbalance_plant_real/Controller-ADR(Model Plant) o This subsystem sends the inputs to the state estimator block and the ADR block o The output of the state estimator is multiplied by the control matrix, K and is then added to the output of the ADR block. ? imbalance_plant_real/Controller-ADR(Model Plant)/Estimator o This subsystem is the state estimator as described in Chapter 4. ? imbalance_plant_real/Controller-ADR(Model Plant)/ADR o This subsystem feeds the inputs through the Gp and Hp subsystems along with the outputs of the subsystem Phi. o The outputs of Gp and Hp are added together and sent out. ? imbalance_plant_real/Controller-ADR(Model Plant)/ADR/Gp Subsystem1 o This subsystem divides the inputs by the sensor gain and then manipulates them as described in Chapter 4. ? imbalance_plant_real/Controller-ADR(Model Plant)/ADR/Phi o This subsystem creates the disturbance vectors as described in chapter 4. ? imbalance_plant_real/Controller-ADR(Model Plant)/ADR/Hp Subsystem1 o This subsystem takes the input vectors from the positions and Phi and manipulates them as described in Chapter 4. o The output of the integrator is sent to another block to obtain H*p as described in Chapter 5. 153 ? imbalance_plant_real/Controller-ADR(Model Plant)/ADR/Hp Subsystem1/Hp Viewer o This subsystem manipulates the outputs of the integrator from the Hp Subsystem1 block to obtain H*p as explained in Chapter 5. ? imbalance_plant_real/Imbalance-Gravity Generator o This subsystem creates the imbalance and gravity vectors derived in Chapter 3. ? imbalance_plant_real/Integral Control o This subsystem performs the integral control on the control voltages as described in Chapter 4. ? imbalance_plant_real/CV Splitter o This subsystem splits up the four control voltages and sends out eight outputs as described in Chapter 4. 154 % state_space_setup_LRSR_magbear_real.m % Simulink Plant variables % 6/6/05 % Kelly Barber clear all close all %------------------variable intitialization------------------------------ N=90; %number of turns per pole mu=pi*4e-7; %permeability of free space Ag=1.176e-3/2; %CS area per pole (m^2) g=0.25e-3; %air gap (m) M=5.414637; %mass of wheel and rotor (kg)(from solid edge) m=0.001; %mass of imbalance (kg)(estimated) Ka=1.2; %amplifier gain (A/V) Ks=23610; %sensor gain (V/m) Ip=0.004903; %moment of inertia (principal)-kg-m^2 (from solid edge) It=0.044950; %moment of inertia (transverse)-kg-m^2 (from solid edge) L1=0.090805; %length of bearing one to flywheel (m) L2=0.090805; %length of bearing two to flywheel (m) G=9.81; %gravity (m/s^2) vb1=1.2; %bias voltage 1 (v) vb2=1.2; %bias voltage 2 (v) vb3=1.2; %bias voltage 3 (v) vb4=1.2; %bias voltage 4 (v) vb5=1.2; %bias voltage 5 (v) vb6=1.2; %bias voltage 6 (v) vb7=1.2; %bias voltage 7 (v) vb8=1.2; %bias voltage 8 (v) ib1=vb1*Ka; %bias current 1 (A) ib2=vb2*Ka; %bias current 2 (A) ib3=vb3*Ka; %bias current 3 (A) ib4=vb4*Ka; %bias current 4 (A) ib5=vb5*Ka; %bias current 5 (A) ib6=vb6*Ka; %bias current 6 (A) ib7=vb7*Ka; %bias current 7 (A) ib8=vb8*Ka; %bias current 8 (A) z=2*mu*Ag*N^2; %from EM force equation Kp1x=z*(ib1^2+ib2^2)/g^3; %position stiffness 155 Kp2x=z*(ib3^2+ib4^2)/g^3; Kp1y=z*(ib5^2+ib6^2)/g^3; Kp2y=z*(ib7^2+ib8^2)/g^3; Ki=z/g^2; %current stiffness Kg=1.5; %integral gain tau=1/3000; %time constant w=50*(2*pi); %omega((Hz)*2*pi=(rad/sec)) theta=0; %imbalance term theta d=0.002; %imbalance term d(m) %----------------------end variable initialization----------------------- %---------------------------modeled plant-------------------------------- % states follow as such % [[X] % [Y] % [alpha] % [beta] % [X(dot)] % [Y(dot)] % [alpha(dot)] % [beta(dot)] % [Ib1*I1-Ib2*I2] x1 % [Ib3*I3-Ib4*I4] x2 % [Ib5*I5-Ib6*I6] y1 % [Ib7*I7-Ib8*I8]] y2 A=[zeros(4,4),eye(4,4),zeros(4,4); (Kz1x+Kz2x)/M,0,(L1*Kz1x-L2*Kz2x)/M,0,0,0,0,0,Ki/M,Ki/M,0,0; 0,(Kz1y+Kz2y)/M,0,-(L1*Kz1y-L2*Kz2y)/M,0,0,0,0,0,0,Ki/M,Ki/M; (L1*Kz1x-L2*Kz2x)/It,0,(L1^2*Kz1x+L2^2*Kz2x)/It,0,0,0,0,w*Ip/It, L1*Ki/It,-L2*Ki/It,0,0; 0,-(L1*Kz1y-L2*Kz2y)/It,0,(L1^2*Kz1y+L2^2*Kz2y)/It,0,0,-w*Ip/It,0,0,0, -L1*Ki/It,L2*Ki/It; 0,0,0,0,0,0,0,0,-1/tau,0,0,0; 0,0,0,0,0,0,0,0,0,-1/tau,0,0; 0,0,0,0,0,0,0,0,0,0,-1/tau,0; 0,0,0,0,0,0,0,0,0,0,0,-1/tau]; B=[zeros(8,4); (ib1+ib2)*Ka/tau,0,0,0; 0,(ib3+ib4)*Ka/tau,0,0; 156 0,0,(ib5+ib6)*Ka/tau,0; 0,0,0,(ib7+ib8)*Ka/tau]; C= [Ks*eye(4,4),zeros(4,8)]; D= zeros(4,8); %---------------------end model set up-------------------------------- %-------------------begin model control code-------------------------- sys=ss(A,B,C,0); P1=-1000+250i; P2=-1000-250i; P3=-6000; pole=[P1, P2, P1, P2, P1, P2, P1, P2, P3, P3, P3, P3]; %close loop poles (repeated) K=place(sys.a,sys.b,pole); %controller gain L=place(sys.a',sys.c',4*pole)'; %observer/estimator gain T=[1 0 L1 0; %transition matrix [x1,x2,y1,y2]'=T*[x,y,alpha,beta]' 1 0 -L2 0; 0 1 0 -L1; 0 1 0 L2]; %-------------------end model control code-------------------------------- %----------------------end modeled plant---------------------------------- %----------------------real plant----------------------------------------- % states follow as such % [[X] % [Y] % [alpha] % [beta] % [X(dot)] % [Y(dot)] % [alpha(dot)] % [beta(dot)] % [I1] % [I2] 157 % [I3] % [I4] % [I5] % [I6] % [I7] % [I8]] Areal=[zeros(4,4),eye(4,4),zeros(4,8); (Kz1x+Kz2x)/M,0,(L1*Kz1x-L2*Kz2x)/M,0,0,0,0,0,ib1*Ki/M,ib2*Ki/M,ib3*Ki/M, ib4*Ki/M,0,0,0,0; 0,(Kz1y+Kz2y)/M,0,-(L1*Kz1y-L2*Kz2y)/M,0,0,0,0,0,0,0,0,ib5*Ki/M,ib6*Ki/M, ib7*Ki/M,ib8*Ki/M; (L1*Kz1x-L2*Kz2x)/It,0,(L1^2*Kz1x+L2^2*Kz2x)/It,0,0,0,0,w*Ip/It,ib1*L1*Ki/It, ib2*L1*Ki/It,ib3*-L2*Ki/It,ib4*-L2*Ki/It,0,0,0,0; 0,-(L1*Kz1y-L2*Kz2y)/It,0,(L1^2*Kz1y+L2^2*Kz2y)/It,0,0,-w*Ip/It,0,0,0,0,0, ib5*-L1*Ki/It,ib6*-L1*Ki/It,ib7*L2*Ki/It,ib8*L2*Ki/It; zeros(8,8),-1/tau*eye(8,8)]; Breal=[zeros(8,8); Ka/tau,0,0,0,0,0,0,0; 0,-Ka/tau,0,0,0,0,0,0; 0,0,Ka/tau,0,0,0,0,0; 0,0,0,-Ka/tau,0,0,0,0; 0,0,0,0,Ka/tau,0,0,0; 0,0,0,0,0,-Ka/tau,0,0; 0,0,0,0,0,0,Ka/tau,0; 0,0,0,0,0,0,0,-Ka/tau]; Gamma= [zeros(4,5); 0,m,0,1,-sin(pi/4); m,0,1,0,-cos(pi/4); 0,0,0,0*0,0; 0,0,0*0,0,0; zeros(8,5)]; Creal= [Ks*eye(4,4),zeros(4,12)]; Dreal= zeros(4,13); %initial=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];%center position(It-Ip)*w^2 %resting position initial=[-0.25e-3*sin(pi/4),-0.25e-3*cos(pi/4),0,0,0,0,0,0,0,0,0,0,0,0,0,0]; %-------------------------end real plant---------------------------------- 158 % imbalance_fft_analysis.m % 6/7/05 % Kelly Barber ss1=3; %sample start time ss2=20; %sample end time % p=2^12; %fft sample - power of 2 dT=0.0001; %sample rate n1=(ss1/dT+1); %sample start point n2=(ss2/dT+1); %sample end point z=fft(SumHp(1,n1:n2),p); %fft of simout(:,#) f=(1/dT)*(0:(p/2))/p; %frequency in Hz. m=abs(z); %magnitude figure, plot(f,m(1,(1:(p/2)+1))), title('fft - X'), %plots freq vs. mag 159 APPENDIX B EXPERIMENTAL CONTROL CODE Appendix B contains the block diagrams created in Simulink to control the experimental rig as described in sections 6.1 and 6.2. This model is initialized by the m- file Omnicode_setup2, compiled through the real time workshop in Matlab, and then uploaded to the external processor dSPACE. Software called Control Desk allows for the manipulation of variables in the Simulink model and allows the user to save time traces and other signals. Included in the back of this Appendix is the m-file called omnicode_setup2.m and a description of the various Simulink subsystems. The FFT plots as seen in Chapter 6 are done using the FFT analysis m-file given in the previous appendix. Filename: The_Omnicode2.mdl Simulation Parameters Start Time: 0.0 Stop Time: 10800.0 Solver Type: Fixed-Step, Discrete(no continuous states) Fixed Step Size: 0.0001 Model Callbacks Model Initialization Function: Omnicode_setup2 Simulation Start Function: Omnicode_setup2 160 Figure B-1 ? Experimental Model Block Diagram Th e_ Om nic od e2 161 Figure B-2 ? Experimental A/D Converter Block Diagram The_Omnicode2/AD Conv 162 Figure B-3 ? Experimental Position Adjustment Block Diagram Th e_ Om nic od e2 /P os itio n 163 Figure B-4 ? Experimental Controller Block Diagram Th e_ Om nic od e2 /C on tro lle r-A DR (M od el Pla nt) 164 Figure B-5 ? Experimental Estimator Block Diagram Th e_ Om nic od e2 /C on tro lle r-A DR (M od el Pla nt) /E stim ato r 165 Figure B-6 ? Experimental ADR Block Diagram Th e_ Om nic od e2 /C on tro lle r-A DR (M od el Pla nt) /A DR 166 Figure B-7 ? Experimental Disturbance Vector Block Diagram Th e_ Om nic od e2 /C on tro lle r-A DR (M od el Pla nt) /A DR /P hi 167 Figure B-8 ? Experimental Gp Subsystem Block Diagram Th e_ Om nic od e2 /C on tro lle r-A DR (M od el Pla nt) /A DR /G p S ub sy ste m1 168 Figure B-9 ? Experimental Hp Subsystem Block Diagram Th e_ Om nic od e2 /C on tro lle r-A DR (M od el Pla nt) /A DR /H p S ub sy ste m1 169 Figure B-10 ? Experimental H *p Generator Block Diagram Th e_ Om nic od e2 /C on tro lle r-A DR (M od el Pla nt) /A DR /H p S ub sy ste m1 /S um of S qu are s 170 Figure B-11 ? Experimental Integral Controller Block Diagram Th e_ Om nic od e2 /In teg ral Co ntr ol 171 Figure B-12 ? Experimental Imbalance Generator Block Diagram Th e_ Om nic od e2 /Im ba lan ce G en era tor 172 Figure B-13 ? Experimental Ramp Generator for Translational Imbalance Block Diagram Th e_ Om nic od e2 /Im ba lan ce G en era tor /d Ra mp 173 Figure B-14 ? Experimental Ramp Generator for Rotational Imbalance Block Diagram Th e_ Om nic od e2 /Im ba lan ce G en era tor /Th eta R am p 174 Figure B-15 ? Experimental Control Voltage Splitter Block Diagram The_Omnicode2/Imbalance Generator/CV Splitter 175 Figure B-16 ? Experimental D/A Converter Block Diagram The_Omnicode2/Imbalance Generator/DA Conv 176 Description of Simulink System and Subsystem Blocks ? The_Omnicode2 o This block diagram serves as the controller for the active magnetic bearing described in sections 6.1 and 6.2. o The input is fed from the dSPACE analog to digital converter. o The output is sent out through the dSPACE digital to analog converter. o The switch voltage is sent through a logical operator to determine if the control voltage will be sent out or not. ? The_Omnicode2/AD Conv o This subsystem takes in the inputs from dSPACE, multiplies them by 10, and then splits off the signal from the switch. ? The_Omnicode2/Position o This subsystem takes the six position readings from the proximitors and multiplies them by the Adjust matrix to obtain four meaning position readings. o There is a calibration block for calibrating the position signals as needed. o The four position vectors are then divided by the sensor gain and multiplied by the Collocation matrix to adjust for collocation issues. ? The_Omnicode2/Controller-ADR(Model Plant) o This subsystem takes the inputs and feds them into the state estimator and the ADR block. o The position inputs are combined with the estimated velocities and then multiplied by the control matrix, K. o The output of the ADR block and the control matrix are added together. ? The_Omnicode2/Controller-ADR(Model Plant)/Estimator o This subsystem estimates the states as described in Chapter 4. ? The_Omnicode2/Controller-ADR(Model Plant)/ADR o This subsystem takes the position inputs and the output of the subsystem Phi and feeds them into the Gp and Hp subsystems. o The output of the Gp and Hp subsystems is added together and sent out. ? The_Omnicode2/Controller-ADR(Model Plant)/ADR/Phi o This subsystem creates the disturbance vectors for the ADR as described in Chapter 4. ? The_Omnicode2/Controller-ADR(Model Plant)/ADR/Gp Subsystem1 o This subsystem performs the mathematical operations as described in Chapter 4. 177 ? The_Omnicode2/Controller-ADR(Model Plant)/ADR/Hp Subsystem1 o This subsystem takes the positions and disturbance vectors and performs the mathematical procedures as described in Chapter 4. o The output from the integrator is fed into the Sum of Squares block to obtain H*p. ? The_Omnicode2/Controller-ADR(Model Plant)/ADR/Hp Subsystem1/Sum of Squares o This subsystem performs the necessary operations to obtain H*p as described in Chapter 5. ? The_Omnicode2/Integral Control o This subsystem performs the integral control as explained in Chapter 4. o The switch value resets the integrator when the switch is off ? The_Omnicode2/Imbalance Generator o This subsystem creates the imbalance vectors as described in Chapter 3. o The subsystem is setup so that the imbalance parameters such as frequency, mass, initial distance and angle, and imbalance activations can be manipulated on-line by Control Desk. ? The_Omnicode2/Imbalance Generator/d Ramp o This subsystem allows the user to utilize Control Desk to manipulate translational imbalance parameters such as final mass distance and the ramp time of the imbalance increase. ? The_Omnicode2/Imbalance Generator/Theta Ramp o This subsystem allows the user to utilize Control Desk to manipulate rotational imbalance parameters such as final imbalance angle and the ramp time of the angle increase. ? The_Omnicode2/CV Splitter o This subsystem takes the final four control voltages and splits them into eight separate outputs as described in Chapter 4. ? The_Omnicode2/DA Conv o This subsystem takes the eight final control voltages and passes them through a saturation block. o The voltages are then either multiplied by 1 or 0 depending on the switch value. o The voltages are finally fed through separate gains for the option of canceling them individually, multiplied by 0.1, and then sent to the dSPACE digital to analog converter. 178 % Omnicode_setup2 % Plant variables % 6/6/05 % Kelly Barber clear all close all %------------------variable intitialization------------------------------ N=90; %number of turns per pole mu=pi*4e-7; %permeability of free space Ag=1.176e-3/2; %CS area per pole (m^2) g=0.25e-3; %air gap (m) M=5.084; %mass of rotor (kg) m=0.001; %mass of simulated imbalance Ka=1.2; %amplifier gain (A/V) Ks=23610; %sensor gain (V/m) Ip=0.00000021085; %moment of inertia (principal)-kg-m^2 (from solid edge) It=0.0000054546; %moment of inertia (transverse)-kg-m^2 (from solid edge) L1=0.090805; %length of bearing one to flywheel (m) L2=0.090805; %length of bearing two to flywheel (m) Adjust=[0.5 -0.5 0 0 0 0; %transforms 6 prox reading into 4 positions 0 0 1 0 0 0; 0 0 0 0.5 -0.5 0; 0 0 0 0 0 1]; Collocation=[1.188 0 -0.188 0; %deals with collocation issues 0 1.188 0 -0.188; -0.188 0 1.188 0; 0 -0.188 0 1.188]; vb1=1; %bias voltage 1 (v) vb2=1; %bias voltage 2 (v) vb3=1; %bias voltage 3 (v) vb4=1; %bias voltage 4 (v) vb5=1; %bias voltage 5 (v) vb6=1; %bias voltage 6 (v) vb7=1; %bias voltage 7 (v) vb8=1; %bias voltage 8 (v) ib1=vb1*Ka; %bias current 1 (A) ib2=vb2*Ka; %bias current 2 (A) ib3=vb3*Ka; %bias current 3 (A) 179 ib4=vb4*Ka; %bias current 4 (A) ib5=vb5*Ka; %bias current 5 (A) ib6=vb6*Ka; %bias current 6 (A) ib7=vb7*Ka; %bias current 7 (A) ib8=vb8*Ka; %bias current 8 (A) z=2*mu*Ag*N^2; %from EM force equation Kp1x=z*(ib1^2+ib2^2)/g^3; %position stiffness Kp2x=z*(ib3^2+ib4^2)/g^3; Kp1y=z*(ib5^2+ib6^2)/g^3; Kp2y=z*(ib7^2+ib8^2)/g^3; Ki=z/g^2; %current stiffness Kg=1.5; %integral gain tau=1/3000; %time constant w=0; %omega %----------------------end variable initialization----------------------- %---------------------------modeled plant-------------------------------- % states follow as such % [[X] % [Y] % [alpha] % [beta] % [X(dot)] % [Y(dot)] % [alpha(dot)] % [beta(dot)] A=[zeros(4,4),eye(4,4); (Kz1x+Kz2x)/M,0,(L1*Kz1x-L2*Kz2x)/M,0,0,0,0,0; 0,(Kz1y+Kz2y)/M,0,-(L1*Kz1y-L2*Kz2y)/M,0,0,0,0; (L1*Kz1x-L2*Kz2x)/It,0,(L1^2*Kz1x+L2^2*Kz2x)/It,0,0,0,0,w*Ip/It; 0,-(L1*Kz1y-L2*Kz2y)/It,0,(L1^2*Kz1y+L2^2*Kz2y)/It,0,0,-w*Ip/It,0]; B=[zeros(4,4); 0,Ki*Ka/M,0,Ki*Ka/M; Ki*Ka/M,0,Ki*Ka/M,0; 0,L1*Ki*Ka/It,0,-L2*Ki*Ka/It; -L1*Ki*Ka/It,0,L2*Ki*Ka/It,0]; 180 C=[eye(4,4),zeros(4,4)]; D=zeros(4,4); %---------------------end model set up-------------------------------- %-------------------begin model control code-------------------------- sys=ss(A,B,C,0); P1=-750+250i; P2=-750-250i; pole=[P1, P2, P1, P2, P1, P2, P1, P2]; % close loop poles (repeated) K=place(sys.a,sys.b,pole); % controller gain L=place(sys.a',sys.c',4*pole)'; % observer/estimator gain T=[0 1 0 -L1; %transition matrix [y1,x1,y2,x2]'=T*[x,y,alpha,beta]' 1 0 L1 0; 0 1 0 L2; 1 0 -L2 0]; %-------------------end model control code-------------------------------- %----------------------end modeled plant---------------------------------- 181 APPENDIX C MAGNETIC BEARING AND TEST WHEEL DRAWINGS The following drawings of the magnetic bearing as described in section 6.1 were made by Alex Matras and can be found in Matras [15]. Those drawings are intended to provide a general idea of the dimensions and construction of the Auburn University magnetic bearing. The drawing for the test wheel was made by this author and is intended as a construction drawing. 182 Figure C-1 ? Drawing of Magnetic Bearing Base Housing 183 Figure C-2 ? Drawing of Magnetic Bearing Probe Housing 184 Figure C-3 ? Drawing of Isometric Views for Magnetic Bearing Probe Housing 185 Figure C-4 ? Drawing of Magnetic Bearing Stator Housing 186 Figure C-5 ? Drawing of Isometric View for Magnetic Bearing Stator Housing 187 Figure C-6 ? Drawing of Magnetic Bearing Stator Housing Face Plate 188 Figure C-7 ? Drawing of Magnetic Bearing Stator 189 Figure C-8 ? Drawing of Magnetic Bearing Assembly Order 190 Figure C-9 ? Drawing of Magnetic Bearing Rotor 191 Figure C-10 ? Drawing of Test Wheel 192 APPENDIX D EQUIPMENT AND MANUFACTURER INFORMATION Amplifiers (8) Copley Controls Corp. DC Brush Servo Amplifiers Model No: 4212Z Power Supplies (2) Copley Controls Corp. Unregulated DC Power Supplies Model No: 666 Proximitors (6) Bentley Nevada, Inc. Series 3300XL 5/8 mm Proximitor Part No: 330180-50-00 Proximity Probes (6) Bentley Nevada, Inc. Series 3300XL Proximity Probes Part No: 330101-08-16-05-02-00 Motor (1) Baldor Electric Company Super-E Industrial Motor Cat. No: EM3545 410 University Ave. Westwood, MA 02090 Tel: (781) 329-8200 Fax: (781) 329-4055 http://www.copleycontrols.com 1631 Bentley Parkway S. Minden, NV 89423 Tel: (775) 782-3611 http://www.bentleynevada.com P.O. Box 2400 Fort Smith, AR 72901 Tel: (800) 828-4920 http://www.baldor.com 193 Speed Control (1) Baldor Electric Company Adjustable Speed Drive Cat. No: ID15J101-ER Controller Hardware (1) dSPACE Inc. Processor: DS1005 D/A Board: DS2103 A/D Board: DS2003 28700 Cabot Dr., Suite 1100 Novi, MI 48377 Tel: (248) 567-1300 Fax: (248) 567-0130 http://www.dspaceinc.com