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

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.

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___________________________________

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

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..........................................................

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)...................................

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..............................

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

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)................

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

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.

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.

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)

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

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

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)

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.

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

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.

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