DESIGN OPTIMIZATION AND HEALTH MONITORING FOR MULTI-DIRECTIONAL
COMPOSITE FLYWHEEL ENERGY STORAGE SYSTEMS
Except where reference is made to the work of others, the work described in this dissertation
is my own or was done in collaboration with my advisory committee. This dissertation
does not include proprietary or classified information.
_____________________________________
Xin Hai
Certificate of Approval:
__________________________ _____________________________
Yasser A. Gowayed George T. Flowers, Chair
Professor Professor
Polymer and Fiber Engineering Mechanical Engineering
_________________________ _____________________________
Subhash C. Sinha Malcolm J. Crocker
Professor Distinguished University Professor
Mechanical Engineering Mechanical Engineering
_________________________ _____________________________
Amnon J. Meir Stephen L. McFarland
Professor Acting Dean
Mathematics Graduate School
DESIGN OPTIMIZATION AND HEALTH MONITORING FOR MULTI-DIRECTIONAL
COMPOSITE FLYWHEEL ENERGY STORAGE SYSTEMS
Xin Hai
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of
the Requirements for the
Degree of
Doctor of Philosophy
Auburn, Alabama
August 7, 2006
iii
DESIGN OPTIMIZATION AND HEALTH MONITORING FOR MULTI-DIRECTIONAL
COMPOSITE FLYWHEEL ENERGY STORAGE SYSTEMS
Xin Hai
Permission is granted to Auburn University to make copies of this dissertation at its
description, upon the request of individuals or institutions and at their expense.
The author reserves all publication rights.
___________________________
Signature of Author
___________________________
Date of Graduation
iv
VITA
Xin Hai, daughter of Lingxi Zheng and Fei Hai, was born on December 10, 1970,
in Zhejiang province, China. She entered Zhejiang University ? in Hangzhou, China ?
in September 1989. After four years of study, she graduated with a Bachelor of Science
degree in Mechanical Engineering in July 1993. Then, she worked as a Mechanical
Engineer in Hangzhou Angel Limit Co. in China from August 1993 to February 1998. As
a senior Mechanical Engineer she worked in the Hangzhou Wahaha Group Co. in China
from February 1998 to April 1999. She entered the Graduate School of Auburn
University in September 1999 and completed a Master of Science in Mechanical
Engineering degree in May 2002. She began to pursue her doctoral program in September
2002. She married Xundan Shi, son of Dingbiao Shi and Guoying Shao, on June 19,
1996. Xin Hai and Xundan Shi have two sons, Kevin H. Shi and Alex H. Shi who are five
and one years old, respectively.
v
DISSERTATION ABSTRACT
DESIGN OPTIMIZATION AND HEALTH MONITORING FOR MULTI-DIRECTIONAL
COMPOSITE FLYWHEEL ENERGY STORAGE SYSTEMS
Xin Hai
Doctor of Philosophy, August 7, 2006
(Master of Science, Auburn University, 2002)
(B.S., Zhejiang University, 1993)
248 Typed Pages
Directed by George T. Flowers
Flywheels are used to store kinetic energy through the rotation of an inertial mass,
which usually rotates at very high speeds. Generally, that energy can be extracted (using a
motor-generator system) in the form of electricity. For space applications, where most
advanced designs are used, weight is a major consideration. Storing as much energy as
possible in as little mass as possible is a fundamental requirement for such systems.
Rotational speeds for certain designs may exceed 60,000 RPM, which introduce very high
stress levels into the hub and rim. Failure of the structure can result in an almost explosive
disintegration, which can prove deadly to anyone in the vicinity of such a disaster.
Accordingly, two major and often contradictory challenges in flywheel design are improved
energy storage density and operational safety. In order to address these challenges, a new
vi
approach to strengthen rotating structures with additional reinforcement in the radial
direction, along with the typical hoop direction reinforcement, which is called multi-
direction composites (MDC), has been developed. This work explores basic issues of safety
and energy storage density with specific application to the MDC approach. Specifically, the
goals of this research effort were to:
? Develop a methodology to optimize the design, based upon maximizing the
energy storing density (ESD) of the complete flywheel system (including
shaft, hub, and rim).
? Design and evaluate a prototype configuration consisting of a steel shaft, a
multiple-layer domed hub and a multi-directional composite rim.
? Develop and evaluate a health monitoring methodology for such systems
based upon variations in imbalance level.
? Investigate the application of the Gabor analysis technique, a joint time-
frequency method, as a vibration analysis tool for crack identification.
? Investigate the natural frequency variations due to the presence of the cracks
in a composite disk rotor system, both experimentally and numerically.
vii
ACKNOWLEDGEMENTS
The author would like to express her deep sincere thanks to her advisor, Dr.
George T. Flowers, Professor, Department of Mechanical Engineering, for his
instruction, guidance, encouragement and patience in the completion of the research and
dissertation. His suggestions, criticism and general exchange of ideas contributed much
to this dissertation.
Thanks are also due to Dr. Yasser Gowayed, Professor, Department of Polymer
and Fiber Engineering, Dr. Subhash Sinha, Professor, Dr. Malcolm Crocker,
Distinguished University Professor, Department of Mechanical Engineering, and Dr.
Amnon. J. Meir, Professor, Department of Mathematics, for their academic guidance and
kind help.
Furthermore, my sincere thanks go to my husband, Xundan Shi, who gave me
great support and courage to help me finish my doctoral program. Also I want to thank
my family members, parents, parents-in-law and two lovely sons. They provided me with
inspiration and endless support during the completion of my Ph.D. degree in Mechanical
Engineering, at Auburn University.
viii
Style manual or journal used: Guide to Preparation and Submission of Thesis and
Dissertation
Computer software used: MS Word2000, ANSYS5.7
ix
TABLE OF CONTENTS
LIST OF FIGURES. . . . . . . . . . . . . . . . . . . . . . . . . . . xv
LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . .xxi
NOMENCLATURE . . . . . . . . . . . . . . . . . . . . . . . . . xxiii
CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW. . . . . . . . . 1
1.1 History of the Flywheel. . . . . . . . . . . . . . . . . . . . 3
1.2 Use of Flywheel Technology for Mechanical Energy Storage. . . . . 3
1.3 The Use of Composite Materials for Flywheels. . . . . . . . . . . 5
1.4 Structure Optimization of a Flywheel Energy Storage System. . . . .7
1.4.1 Terms of Flywheel Energy Storage System. . . . . . . . . .8
1.4.2 Flywheel Energy Storage System Optimal Design. . . . . . . 9
1.4.3 Flywheel Energy Storage System Optimal Procedure. . . . . 13
1.4.4 Design Optimization using Finite Element Analysis. . . . . 14
1.5 Safe Operation of a Flywheel Energy Storage System. . . . . . . . 15
1.6 Experimental Crack Investigations of the Flywheel Energy Storage
System. . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.7 Numerical Investigations of the Influence of a Crack on the Vibration
x
Flywheel Energy Storage System. . . . . . . . . . . . . . . . . 18
1.8 The Structure of the Dissertation. . . . . . . . . . . . . . . 20
1.9 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . 20
CHAPTER 2 DESIGN OPTIMIZATION METHODOLOGY FOR MULTI-
DIRECTIONAL COMPOSITE FLYWHEEL HUBS. . . . . . . . 26
2.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . .26
2.2 Analysis. . . . . . . . . . . . . . . . . . . . . . . . . .29
2.2.1 Configuration. . . . . . . . . . . . . . . . . . . . .30
2.2.2 Coordinate Systems. . . . . . . . . . . . . . . . . . 30
2.3 Optimization Design Approach. . . . . . . . . . . . . . . . .31
2.3.1 Design Variables. . . . . . . . . . . . . . . . . . . .32
2.3.2 State Variables. . . . . . . . . . . . . . . . . . . . 32
2.3.3 Constraint Requirements. . . . . . . . . . . . . . . . 32
2.3.4 Objective Function. . . . . . . . . . . . . . . . . . .33
2.3.5 Governing Equations for Optimization Processing. . . . . . 33
2.3.6 Design Optimization Cycles. . . . . . . . . . . . . . .35
2.4 Results and discussion. . . . . . . . . . . . . . . . . . . . 36
2.4.1 Domed-hub Model Results. . . . . . . . . . . . . . . 37
2.4.2 Optimization Results of 2-D Domed Hub-Shaft-Rim Model with
the Thickness of the Hub Port at 4.8 mm and 6.0 mm. . . . . 38
2.4.3 Discussion. . . . . . . . . . . . . . . . . . . . . . 39
xi
2.5 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . .40
CHAPTER 3. HEALTH MONITORING METHODOLOGY FOR MDC FLYWHEEL
SYSTEMS BASED UPON VARIATIONS IN IMBALANCE
ECCENTRICITY LEVEL . . . . . . . . . . . . . . . . . . .64
3.1 Computational procedure . . . . . . . . . . . . . . . . . . 64
3.2 The calculation of imbalance eccentricity of the mass center (U
m
) . . . 67
3.3 Results and discussion. . . . . . . . . . . . . . . . . . . .73
3.3.1 Relationship between U
m
and rotation speed . . . . . . . . 73
3.3.2 Relationship between U
m
and various material properties . . .74
3.3.3 Relationship between U
m
and crack size . . . . . . . . . 75
3.3.4 Relationship between U
m
and ? . . . . . . . . . . . . . 75
3.4 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . .76
CHAPTER 4. GABOR TECHNIQUES AS A TOOL FOR THE VIBRATION
ANALYSIS AND HEALTH MONITORING OF CRACKS IN ROTOR
SYSTEMS. . . . . . . . . . . . . . . . . . . . . . . . .92
4.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . 92
4.2 Gabor Analysis. . . . . . . . . . . . . . . . . . . . . . . 96
4.2.1 FFT Drawbacks. . . . . . . . . . . . . . . . . . . . 96
4.2.2 Discrete Gabor Expansion. . . . . . . . . . . . . . . .97
xii
4.2.3 Comparison of Wavelet and Gabor Analysis. . . . . . . . 99
4.3 Experimental Results. . . . . . . . . . . . . . . . . . . . 101
4.3.1 Natural Frequency Measurements. . . . . . . . . . . . 102
4.3.2 Rotor Vibration Analysis. . . . . . . . . . . . . . . .103
4.4 Conclusions. . . . . . . . . . . . . . . . . . . . . . . .104
CHAPTER 5. VIBRATION ANALYSIS AND HEALTH MONITORING OF CRACKS
IN COMPOSITE DISK ROTOR SYSTEMS. . . . . . . . . . . 115
5.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . 115
5.2 Experimental Investigation. . . . . . . . . . . . . . . . . . 116
5.2.1 Testing Sample. . . . . . . . . . . . . . . . . . . .116
5.2.2 Testing Rig. . . . . . . . . . . . . . . . . . . . . 117
5.2.3 Experimental Results and Discussions. . . . . . . . . . 117
5.3 Numerical Simulation Using Finite Element Analysis. . . . . . . 118
5.3.1 3-D Model. . . . . . . . . . . . . . . . . . . . . .118
5.3.2 Grid Independence Study. . . . . . . . . . . . . . . .119
5.3.3 Comparison of the Results between the Experimental and
Simulation. . . . . . . . . . . . . . . . . . . . . 120
5.3.4 Parameters of Numerical Simulations. . . . . . . . . . .121
5.3.5 Numerical Simulation Results and Discussions. . . . . . . 122
5.3.5.1 Natural Frequencies Variation with a Variable Length
Crack along the Disk Hoop Direction. . . . . . . . 122
xiii
5.3.5.2 Natural Frequencies Variation with a Constant Length
Crack at Different Locations along the Disk Radial
Direction. . . . . . . . . . . . . . . . . . . . 123
5.3.5.3 Natural Frequencies Variation with a Variable Length
Crack along the Disk Radial Direction. . . . . . . 125
5.4 Simplified Rotor Model and Simulation. . . . . . . . . . . . . 127
5.4.1 Positions of the Disk Masses. . . . . . . . . . . . . . 127
5.4.2 Velocities of the Disk Masses. . . . . . . . . . . . . . 127
5.4.3 Kinetic and Potential Energy of the Rotor. . . . . . . . . 128
5.4.4 Equations of Motion of the Rotor. . . . . . . . . . . . 129
5.4.5 Study of Numerical Simulation. . . . . . . . . . . . . 131
5.4.6 Results and Discussions. . . . . . . . . . . . . . . . 131
5.5 Conclusions. . . . . . . . . . . . . . . . . . . . . . . .133
CHAPTER 6. CONCLUSIONS AND RECOMMENDATIONS. . . . . . . . . 164
REFERENCES. . . . . . . . . . . . . . . . . . . . . . . . . . . . .168
APPENDICES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
A. C++ CODE FOR ?CALCULATIONS OF HUB GEOMETRY AND
MATERIAL PROPERTIES?. . . . . . . . . . . . . . . . . . 181
B. ANSYS CODE FOR ?CONSTRUCTION AND OPTIMIZATION OF
COMPOSITE SHAFT-HUB-RIM ANSYS MODEL?. . . . . . . . 194
xiv
C. FORTRAN CODE FOR ?CALCULATION OF MAXIMUM ENERGY
STORING DENSITIY?. . . . . . . . . . . . . . . . . . . . 200
D. FORTRAN CODE FOR ?IMBALANCE ECCENTRICITY OF THE MASS
CENTER CALCULATIONS?. . . . . . . . . . . . . . . . . . 204
E. MATLAB CODE FOR ?SIMULATION OF CRACKED DISK MODEL?
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212
xv
LIST OF FIGURES
Figure 1.1 MDC flywheel with assembled layers . . . . . . . . . . . . . .23
Figure 1.2 MDC flywheel stress distribution direction . . . . . . . . . . . 24
Figure 1.3 Geometrical dimensions of the MDC flywheel system . . . . . . . 25
Figure 2.1 2-D configuration of the domed-hub for the flywheel system. . . . . .42
Figure 2.2 3-d configuration of the shaft/domed-hub/rim for the flywheel system. .43
Figure 2.3 The scheme of global cylinder coordinate system and local element
coordinate systems. . . . . . . . . . . . . . . . . . . . . . 44
Figure 2.4 Optimization design procedure of a composite domed hub-shaft-rim
flywheel system. . . . . . . . . . . . . . . . . . . . . . .45
Figure 2.5 Flowchart of ANSYS optimization procedure. . . . . . . . . . . .46
Figure 2.6 ANSYS mesh grid systems of the 2-D domed composite hub. . . . . 47
Figure 2.7 S
X
stress distribution of the 2-D domed hub model. . . . . . . . . 48
Figure 2.8 S
XY
stress distribution of the 2-D domed hub model. . . . . . . . .49
Figure 2.9 S
Y
stress distribution of the 2-D domed hub model. . . . . . . . . 50
Figure 2.10 U
X
deformation of the 2-D domed hub model. . . . . . . . . . . 51
Figure 2.11 U
X
deformation distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 4.8 mm and the ratio of the elliptical radii
a/b = 2.01. . . . . . . . . . . . . . . . . . . . . . . . . 52
xvi
Figure 2.12 S
X
stress distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 4.8 mm and the ratio of the elliptical radii
a/b = 2.01. . . . . . . . . . . . . . . . . . . . . . . . . 53
Figure 2.13 S
XY
stress distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 4.8 mm and the ratio of the elliptical radii
a/b = 2.01. . . . . . . . . . . . . . . . . . . . . . . . . 54
Figure 2.14 S
Y
stress distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 4.8 mm and the ratio of the elliptical radii
a/b = 2.01. . . . . . . . . . . . . . . . . . . . . . . . . 55
Figure 2.15 U
SUM
deformation distribution of 2-D domed hub-shaft-rim model with
the thickness of the hub port t = 6.0 mm and the ratio of the elliptical radii
a/b = 1.99. . . . . . . . . . . . . . . . . . . . . . . . . 56
Figure 2.16 U
X
deformation distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 6.0 mm and the ratio of the elliptical radii
a/b = 1.99. . . . . . . . . . . . . . . . . . . . . . . . . 57
Figure 2.17 S
XY
stress distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 6.0 mm and the ratio of the elliptical radii
a/b = 1.99. . . . . . . . . . . . . . . . . . . . . . . . . 58
Figure 2.18 Distribution of the allowable tip speed vs. the ratio of a/b (The maximum
shear stress is equal to the criterion shear stress for the flywheel rotor)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
Figure 2.19 Distribution of the maximum shear stress vs. the ratio of a/b (The rotating
speed ?= 60,000 rpm) . . . . . . . . . . . . . . . . . . . . .60
xvii
Figure 2.20 Distribution of the maximum ESD vs. the ratio of a/b. . . . . . . .61
Figure 3.1 Illustration of hoop-direction crack. . . . . . . . . . . . . . .78
Figure 3.2 Illustration of radial-direction crack. . . . . . . . . . . . . . . 79
Figure 3.3 Top view of crack model. . . . . . . . . . . . . . . . . . . 80
Figure 3.4 3-D model. . . . . . . . . . . . . . . . . . . . . . . . . 81
Figure 3.5 3-D model. . . . . . . . . . . . . . . . . . . . . . . . . 81
Figure 3.6 Configuration of node mass distribution. . . . . . . . . . . . . .82
Figure 3.7 Relationship between U
m
and rotation speed. . . . . . . . . . . .83
Figure 3.8 Relationship between U
m
and different material properties
(Rotation speed: ? = 60,000 rpm). . . . . . . . . . . . . . . . 84
Figure 3.9 Relationship between U
m
and crack length
(l: length of crack, c: circumference, rotation speed: ? = 60,000 rpm)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
Figure 3.10 Relationship between U
m
and ? in hoop layer radial crack
(Rotation speed: ? = 60,000 rpm). . . . . . . . . . . . . . . 86
Figure 4.1 Drawbacks of the Fourier Transform ---- lose time information. . . .106
Figure 4.2 Drawbacks of the Fourier Transform ---- can?t track frequency change 107
Figure 4.3 Gabor analysis used for frequencies distribution. . . . . . . . . .108
Figure 4.4 Comparisons of the wavelet analysis and Gabor analysis
(a) A typical run-up signal of a cracked rotating shaft;
(b) The wavelet scalogram; (c) The Gabor spectrogram. . . . . . . 109
Figure 4.5 Experimental setup for crack detection. . . . . . . . . . . . .110
Figure 4.6 Comparison of the run-up signals of intact (a) and cracked (c) shafts and
xviii
their corresponding spectrograms (b,d) . . . . . . . . . . . . .111
Figure 4.7 The Gabor coefficients of an intact shaft (a) overall frequency range; (b)
the fundamental frequency; (c) the second harmonic; (d) the third
harmonic; and (e) the FFT-based power spectrum. . . . . . . .112
Figure 4.8 The Gabor coefficients of a shaft with a 5 mm deep crack in the middle.
(a) overall frequency range; (b) the fundamental frequency and the
sideband components; (c) the second harmonic of the running speed; and
(d) the FFT-based power spectrum. . . . . . . . . . .113
Figure 5.1 Schematic diagram of experimental test article. . . . . . . . . . .135
Figure 5.2 Experimental setup for the frequency measurement of the composite
crack disk with heavy rim (a) Out-of-plane direction test setup
(b) In-plane direction test setup. . . . . . . . . . . . . . . . .136
Figure 5.3 FEA model (a) Drawing and grid system; (b) Fourth mode shape with the
associated frequency 52.5 Hz for the crack case; (c) Fifth mode shape with
the associated frequency 315 Hz (in-crack direction); (d) Sixth mode
shape with the associated frequency 321 Hz (off-crack direction). . .137
Figure 5.4 Top view of crack disk model. . . . . . . . . . . . . . . . . .138
Figure 5.5 The mesh grid changes along the disk hoop direction (at n
z
= 4, n
R
= 24, a
hoop crack lies in the middle radius with the angle 60?). . . . . . .139
Figure 5.6 The mesh grid changes in the disk thickness direction (at n
?
= 84, n
R
= 24,
a hoop crack lies in the middle radius with the angle 60?). . . . . .140
Figure 5.7 The mesh grid changes along the disk radial direction (at n
z
= 6, n
?
= 84, a
hoop crack lies in the middle radius with the angle 60?) . . . . . .141
xix
Figure 5.8 The path for crack growing in hoop direction. The crack start point P
1
is
fixed, the other end extends along the middle radius ( r = ? (Ri+Ro) )
circumference C
r
, till point P
2
. . . . . . . . . . . . . . . . . .142
Figure 5.9 Natural Frequencies variation with a variable crack length along the disk
hoop direction. . . . . . . . . . . . . . . . . . . . . . . .143
Figure 5.10 Natural Frequency variation with the variable crack length along the disk
hoop direction (all modes) . . . . . . . . . . . . . . . . . . .144
Figure 5.11 The path for crack moving in radial direction. The crack length is fixed,
the crack center position starts from point P
2
, moves to P
3
. . . . . .145
Figure 5.12 Natural Frequency variation with a constant length crack at different
locations along the disk radial direction. . . . . . . . . . . . . 146
Figure 5.13 Natural Frequency variation with a constant length crack at different
locations along the disk radial direction (all modes). . . . . . . . 147
Figure 5.14 The path for crack growing in radial direction. The crack center is fixed at
point P
1
, the two ends of crack extend to points P
2
and P
3
. . . . . . .148
Figure 5.15 Natural Frequency variation with a variable length crack along the disk
radial direction. . . . . . . . . . . . . . . . . . . . . . . .149
Figure 5.16 Natural Frequency variation with a variable length crack along the disk
radial direction (all modes). . . . . . . . . . . . . . . . . . 150
Figure 5.17 Simplified rotor model with flexible hub. . . . . . . . . . . . .151
Figure 5.18 Coordinate systems for support frame and rotating shaft frame. . . . 152
Figure 5.19 Influence of a crack on forward and backward whirling natural
frequencies . . . . . . . . . . . . . . . . . . . . . . . . .153
xx
Figure 5.20 Natural frequency variation of the primary mode with stiffness ratio for
selected operating speeds (a) Forward; (b) Backward. . . . . . . .154
xxi
LIST OF TABLES
Table 2.1 Optimization of the flywheel rotor with the thickness of the hub port
t = 4.8 mm (rotating speed ?= 60,000 rpm). . . . . . . . . . . . 62
Table 2.2 Optimization of the flywheel rotor with the thickness of the hub port
t = 6.0 mm (rotating speed ?= 60,000 rpm). . . . . . . . . . . . 63
Table 3.1 Material properties of the 3-D rim model . . . . . . . . . . . . 87
Table 3.2 The U
m
at different rotational speeds in different layers. . . . . . . .88
Table 3.3 The U
m
due to different material properties (Rotation speed: ? = 60,000
rpm). . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Table 3.4 The U
m
observed for different crack sizes of hoop layer hoop cracks
(Rotation speed: ? = 60,000 rpm). . . . . . . . . . . . . . . .90
Table 3.5 The U
m
due to different value of ? of a hoop layer radial crack
(Rotation speed: ? = 60,000 rpm). . . . . . . . . . . . . . . . 91
Table 4.1 Measured natural frequencies for the intact and cracked shafts. . . . 114
Table 5.1 Material properties of the composite disk. . . . . . . . . . . . 155
Table 5.2 Measured Natural Frequencies for Out-of-Plane 1
st
Bending Mode. . 156
Table 5.3 Measured Natural Frequencies for Out-of-Plane 2
nd
Bending Mode. . 157
Table 5.4 Measured Natural Frequencies for 1
st
In-Plane Mode. . . . . . . .158
Table 5.5 Natural frequencies and mode shapes from finite element analysis. . 159
xxii
Table 5.6 Comparison of natural frequencies from finite element analysis and
experimental testing. . . . . . . . . . . . . . . . . . . . .160
Table 5.7 Explanation of terms used to describe modes. . . . . . . . . . . 161
Table 5.8 Parameters of the simple rotor disk model. . . . . . . . . . . .162
Table 5.9 Comparison of natural frequencies from finite element analysis and
theoretical model simulation. . . . . . . . . . . . . . . . . . 163
xxiii
NOMENCLATURE
Numbers
1,2,3,4 test points for in-plane direction tests
English Symbols
a,b,c,d test points for out-of-plane direction tests
a the long radius of the elliptical edge of the hub for the shaft-hub-rim flywheel, m
b the short radius of the elliptical edge of the hub for the shaft-hub-rim flywheel, m
C
r
the circumference length in the disk middle radius, m
C
x
viscous damping element in the X direction
C
y
viscous damping element in the Y direction
E modulus of elasticity, Pa
E
x
, modulus of elasticity for composite material in the X direction, Pa
E
y
, modulus of elasticity for composite material in the Y direction, Pa
E
z
modulus of elasticity for composite material in the Z direction, Pa
ex the shape factor for the hub design
ESD the energy storing density for the flywheel system, kWhr/kg
F
cj
the natural frequency values for the crack model, Hz, j=1,2,3,4,5,6
xxiv
F
nj
the natural frequency values for the no-crack model, Hz, j=1,2,3,4,5,6
Gxy torsion modulus of elasticity for composite material in the X-Y plane, Pa
Gyz torsion modulus of elasticity for composite material in the Y-Z plane, Pa
Gxz torsion modulus of elasticity for composite material in the X-Z plane, Pa
h the height of the hub or rim, m
I the mass moment of inertia of the rotor around the spin axis, kgm
2
K
i
spring stiffness element connected the discrete i
th
mass element mass and shaft,
N/m, i = 1,2,3,4
K
c
spring stiffness with the crack, N/m
K
x
support spring stiffness in the X direction, N/m
K
y
support spring stiffness in the Y direction, N/m
L
c
the length of the crack, m
m
i
discrete i
th
mass element , kg, i = 1,2,3,4
m
h
the mass of the hub, kg
Nuxy Poisson?s ratio for the composite material in X-Y plane
Nuyz Poisson?s ratio for the composite material in Y-Z plane
Nuxz Poisson?s ratio for the composite material in X-Z plane
n
z
the element numbers divided in the disk thickness direction
n
?
the numbers of the mesh grids along the disk hoop circumference direction
n
R
the numbers of the mesh grids along the disk radial direction
p
i
the coordinate of the discrete i
th
mass element, m , i = 1,2,3,4
Ppot the potential energy of the rotor system, J
R the radius of the hub, m
xxv
R
cc
the radius at the crack center position of the disk sample for the radial direction,
R
i
? R
cc
? R
o
, m
R
i
the inner radius of the disk sample, m
R
o
the outer radius of the disk sample, m
r
o
the outer radius of the hub/ the outer radius of the rim for the shaft-hub-rim
flywheel, m
r
s
the radius of the shaft for the shaft-hub-rim flywheel, m
r
i
the inner radius of the rim for the shaft-hub-rim flywheel , m
r the middle radius of the disk sample, ?(R
o
+R
i
), m
?R the difference between the outer radius and the inner radius, R
o
- R
i
, m
T the kinetic energy of the rotor system
t the thickness of the disk sample /the thickness of the hub flat part for the shaft-
hub-rim flywheel, m
U
m
the imbalance eccentricity center of mass for the disk sample, m
V the total volume of the hub
V
hub
the volume of the ellipse hub for the shaft-hub-rim flywheel, m
3
V
s
the volume of the steel shaft for the shaft-hub-rim flywheel, m
3
V
rim
the volume of the multiple directional rim for the shaft-hub-rim flywheel, m
3
V
i
the velocity component of the discrete i
th
mass element, m/s, i = 1,2,3,4
v
i
the deformation of the discrete i
th
mass element, m, i = 1,2,3,4
X the displacement in the X direction for the support coordinate system, m
Y the displacement in the Y direction for the support coordinate system, m
Z the displacement in the Z direction for the support coordinate system, m
X the displacement in the X direction for the rotating frame coordinate system, m
Y the displacement in the Y direction for the rotating frame coordinate system, m
Z the displacement in the Z direction for the rotating frame coordinate system, m
Greek Symbols
? the rotational speed of the rotor system, rpm or Hz
? the rotational speed of the flywheel system, rpm or Hz
?
i
the phase angle of the discrete i
th
mass element, rad, i = 1,2,3,4
? the density of material, kg/m
3
? Poisson?s ratio
? the ratio of the inner radius/outer radius for the hub, that is a/b
Subscripts
c crack
cc position of the crack center
h hub
i inner or i
th
mass, i = 1,2,3,4
j number of the mode, j=1,2,3,4,5,6
o outer
pot potential energy
r disk middle radius
X X direction
xxvi
Y Y direction
xxvii
Z Z direction
xy X-Y plane
yz Y-Z plane
zx Z-X plane
CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW
Flywheels have been used for thousands of years in many fields, but until the
1970?s flywheels had not been used mostly as energy storage devices. The energy stored
in flywheel usually is the kinetic energy and is proportional to the square of the rotational
speed, so flywheels can be used as energy storage sources for space vehicles, satellite
devices, etc., typically rotate at relatively high speeds in order to store as much energy as
possible. A major milestone occurred during the 1960?s and 70?s when NASA sponsored
programs, which proposed energy storage flywheels as possible primary sources for
space missions.
With the development of composite material manufacturing technology, multi-
directional composites (MDC) with a layered structure can be used in flywheels for space
vehicles. The use of composite materials provides numerous advantages over metals,
such as reduced weight and increased strength, as reported by Gowayed and Flowers
(2002). There are two principles about the light weight for a high speed rotational
flywheel: one advantage is that the bearing support system gets ultra-low friction, and
another advantage is that the inertial inference and thermal loading which causes stress
concentration in the material at high rotational speeds are minimized. High strength is
needed to achieve maximum rotational speed. Figure 1.1 shows the structure of a typical
layered MDC flywheel. The thin layers are made of radial reinforcement composite
1
materials whereas the thick layers are made of hoop reinforcement composite materials.
The hoop reinforcement layer and radial reinforcement layer are assembled alternately in
order to form a composite structure reinforced in both hoop and radial directions. This
approach helps remedy the classical problem of radial crack propagation in thick,
filament wound flywheels and eliminates kinks that can drastically reduce hoop strength.
The MDC flywheel also allows the tip speed to be increased, thus enhancing the energy
capacity by allowing more radial tensile stress distributions.
In order to maximize both the energy storage capacity and the energy storage
density (ESD), it is necessary to optimize the design by applying our understanding of the
static and dynamic behaviors of MDC flywheel systems rotating at high speeds. In this
study, understanding the distribution of the stress, the deformation of the geometry, is
critical in order to decide the optimum structure, geometry, parameters and rotational
speeds for a flywheel system.
Cracks and voids are common defects in high-speed rotating flywheel systems
and are a precursor to fatigue-induced failure. Another goal of this study is to investigate
the initiation and propagation of the cracks both numerically and experimentally. Cracks
are often initiated in or near the regions where have very high local stresses distribution.
When stresses increase to some threshold level or failure point, the leading edge of the
crack (crack tip) starts to propagate. Cracks could cause localized flexibility in the
structure of the rotor and thereby, affect the dynamic behavior to some extent. There are
two basic kinds of cracks: transverse cracks and torsion cracks, which depend on the
types of stress fields. In most rotor systems, bending stress tends to be dominant.
Accordingly, the present study considers only transverse cracks.
2
1.1 History of the Flywheel
The flywheel is one of the mankind?s earliest inventions and has been widely used
for a long time. In ancient China, the four wooden flywheels made the earliest fight
vehicles. In Egyptians, about 2,500 years ago, people created the earliest chariot by using
wooden flywheels. Flywheel systems were widely used in everyday life (Genta 1985):
water wheels for pumping water and wind driven flywheels that generated power. The
historical development of flywheels and their uses has largely been dependent on
advances in both materials and machine technology, coupled with opportunity or
necessity (Horner, et al., 1996). The first great advances were made in the 18
th
century
during the Industrial Revolution. The use of flywheels in stream engines and the
widespread use of metal instead of wood in the construction of machines were two
important developments. Flywheels built of cast iron, with the greater density of material,
could combine a greater mass and moment of inertia in the same space and structure.
There are numerous applications of flywheels; one modern application is the use of
automotive flywheels to keep electric cars running smoothly.
1.2 Use of Flywheel Technology for Mechanical Energy Storage
The origins and use of flywheel technology for mechanical energy storage began
several hundred years ago, and was developed throughout the whole Industrial
Revolution. A flywheel is a unique rotating machine in that the design purpose is to store
energy by the combined effects of speed and mass (Bowler, 1997). It is universally
appreciated as a kinetic form of stored energy. Because the flywheel rotor rotating at a
high-speed stores a large mount of inertia, the flywheel system can be used in many
3
applications that use stored inertia in the flywheel rotor to balance fluctuations in the
velocity of a rotational machine caused by inertia and external forces (Genta, 1985). One
example of this application is the flywheel rotor in steam engines.
A compactly sized and lightweight component with a large amount of energy
storage storing density is very important for a spacecraft. Space vehicle programs
consistently effort to reduce satellite rim mass to increase payload capacity and/or reduce
launch and fabrication costs. At the same time, performance demands on satellite systems
or spaceship continue to increase, creating a challenge for space vehicle technology
development. Although chemical batteries have been used long time in spacecraft design,
mechanical batteries such as flywheels offer significant weight and life benefits.
Flywheel battery systems have the potential to store more amounts of energy with very
little weight compared to chemical batteries and they can be used as attitude control
actuators to replace reaction flywheel assemblies and control moment gyros. Flywheel-
based systems provide both high energy storing density and attitude control functionality,
and solve both of these issues.
Until recent times, the development of flywheels has been limited by a lack of
suitable materials and useful technologies. Steel has been used for a long time in the past
but has limitations on speed/weight ratio and safety for the spaceship project. With the
development of modern fibers and resins, a new range of exciting materials has become
available for the engineer and designer. This has led to many studies in the 1970?s and
80?s on the general concept of a flywheel energy storage system.
In the 1990?s, with the developments in strong lightweight composite materials,
magnetics bearings and solid state electronic devices, new designs are now significantly
4
performing all earlier facts in several important respects, which can be summarized as
follows (Horner (1996)):
? Quick energy recovery
? High efficiency
? Low maintenance and long service life
? High stored energy vs. small volume
? High output power levels
? High stored energy vs. light weight
? Low product and operational costs
In 2000, Truong et al., introduced the Flywheel Energy Storage Demonstration
Project, initiated at the NASA Glenn Research Center, as a possible replacement for the
Battery Energy Storage System on the International Space Station. In 2000, Fausz et al.,
reported that the Flywheel Attitude Control, Energy Transmission and Storage
(FACETS) system should combine all or parts of the energy storage, attitude control, and
power management and distribution (PMAD) subsystems into a single system, thus
significantly decreasing flywheel mass (and volume).
1.3 The Use of Composite Materials for Flywheels
When a flywheel energy storage system is rotating at high speed, the stress
distribution in thin rims is mainly concentric (or hoop), so composite materials with the
high unidirectional strength of their fibers can be used to advantage. A composite
5
flywheel rotor was designed and developed by Potter and Medicott (1986), as part of a
program to develop a flywheel energy storage system, which used in vehicle applications.
In the following decade, advanced composite flywheels have been developed in order to
meet the demands for energy storage systems with both a high energy density and long
life. In the study of Curtiss, et al., (1995), it was shown that the composite Carbon fiber-
Epoxy disk rotor is capable of a 38% higher rim speed or 91% greater rotor energy
density, than a comparable rotor built of high strength weight ratio isotropic Titanium or
steel alloys.
Generally, the flywheel energy storage system consists of a flywheel rotor, a
supporting device, a charge/discharge device, and a safety device. Composite flywheel
energy storage technologies currently compete with advanced electro-chemical batteries
in applications, which require high specific energy and power.
Many investigators in the field have advanced new designs and manufacturing
technologies to increase the storing density and ensure the safety operation of the
flywheel. Among them, technological advances in high strength, low-density composite
materials have been interested in increasing the energy storing density.
Kojima, et al., (1997), proposed a carbon fiber reinforced plastic (CFRP)
flywheel. The study results show that a flywheel, with high-modulus graphite/epoxy
filament wound composites that are designed for better stress distribution, is able to rotate
at a higher speed. In 1999, Huang tested a polar woven flywheel design; this dual
function capability provides weight savings and an improvement in life and reliability of
the total spacecraft system. The life expectancy of a satellite flywheel system is at least
10 years, during which it is expected to undergo 50,000 cycles (charge/discharge) while
6
supplying 3 to 5 kW-h of usable energy. In 2002, a new approach to strengthen flywheels
with additional reinforcement in the radial direction along with the typical hoop direction
reinforcement?Multi-Direction Composite (MDC) flywheel systems?was reported by
Gowayed and Flowers.
1.4 Structure Optimization of a Flywheel Energy Storage System
The efficiency of flywheel energy storage not only depends on the wheel
materials, but also on the connecting and driving devices. Another essential issue in the
design of the flywheel rotor is how the permanent magnet for the motor/generator is
attached to the hub that is connected to the flywheel rotor. The optimum structural design
of rotating flywheels has fascinated many investigators. Published reports of such
investigations date back to as early as 1912. During the 1960?s and 1970?s, with the
development of computers, a number of computational schemes appeared in the literature
for the optimum design of rotating flywheels. However, it was not until the 1990's when
microelectronics, magnetic bearing systems and high power density motor-generators
became available, that flywheel energy storage systems became really practical. During
the past few decades, the structure optimal design and control for the flywheel energy
storage system has been widely studied and many articles on the optimal design and
structure development of the flywheel energy storage system provide an overview of the
current technologies [Lashley et al., 1989; Niemeyer et al., 1989; Jee and Kang 2000;
Sahin et al., 2001; Zheng and Fabien 2002; Shen et al., 2003; Swett et al., 2005; Zheng et
al., 2005].
7
1.4.1 Terms of Flywheel Energy Storage System
Through a series of workshops and meetings, the aerospace flywheel community
has defined some common terms as they appear to high performance (i.e., high specific
energy) flywheels. From the NASA space flywheel program, the following definitions
have been agreed upon:
? Flywheel System (a.k.a.: flywheel unit, or simply "flywheel): All
mechanical and electrical components required to charge or discharge
power.
? Rotor: The complete rotating assembly, include rim, hub, shaft and other
device parts.
? Rim: The main rotating mass of the flywheel rotor system. Rim is
contained in the most of the stored energy
? Hub: The part of the rotor that connects the rim to the shaft. Because the
hub mass contributes minimally to the total stored energy for the flywheel
system, it should be as light as possible.
? Shaft: The shaft is rotating axial part that lies on the spin axis of the rotor
and is the rotor interface connection to the portions of the magnetic
bearings and motor generators. Some flywheel rotor systems may not
require a shaft.
? Flywheel Energy Storage (FES): A flywheel-based system used only for
8
charge and discharge of energy. For spaceship or satellite applications, this
implies the use of at least two flywheels in a counter-rotating
configuration.
? Total Stored Energy (TSE): The TSE of the flywheel system is the total
energy stored in the rotor of the flywheel system (as defined above).
? Energy Storing Density (ESD): Energy stored divided by mass.
1.4.2 Flywheel Energy Storage System Optimal Design
The optimization design of a flywheel energy storage system is a very
complicated and time costing process, which must include consideration of shape, stress
distribution, rotor dynamic characteristics, safety operation and cost, etc. The important
issues to be considered in the optimization design of a flywheel energy storage system are
listed by Horner (1996):
? Power transfer and discharge
? Flywheel shape and dimension
? Energy losses
? Suspension system connection and operation
? Manufacture and product
? Failure management and detection
A reasonable approach to the design is essential to make sure that an efficient
9
system is realized. Many investigations have been conducted to optimize the flywheel
system design. Especially with the development of high strength composite materials and
simulation software skills, a number of computational schemes and methodologies for the
design of flywheel system have been reported in the literature.
The shape or dimension of an energy storage flywheel is generally chosen in a
way to maximize the energy storing density (ESD). Structural shape optimization is a
very important subject in engineering today. Difficulties persist in the effective use of
optimization techniques; some caused by the complexities of structure, some caused by
the conflicting criteria that needs to be taken into account between the high cost and low
life.
Among them, Grudkowski et al., (1996), described the major design
considerations for a composite flywheel. For composite rotor design, it is critical to note
that the stress allowable for the boundary conditions of operation must be used. The
purpose to maximize flywheel energy storing density is affected by the hoop and radial
design values obtained. Therefore, high-strength, low-density fibers are preferred, and
multiple ring structural designs have been used as one approach to overcome the low
strength in the radial direction.
In 1998, Ha et al., (1999), developed one type of technological advancement in
high-strength, low-density composite material that has been influential in increasing the
energy storage capability. Ha?s method of increasing the total energy storage is to
fabricate the hybrid rotor in multiple rims of different composite materials. This design of
the hybrid composite rotor can reduce the tensile stress in the radial direction, thus
overcoming the low strength in that direction without reduction of the hoop directional
10
strength. Ha et al., (1999), also used optimal design methods to select the dimensions of a
hybrid composite flywheel rotor and permanent magnet. That flywheel rotor design
consists of multiple rims of advanced composite materials with a permanent magnetic
rotor attached inside the flywheel rim. In this way, the centrifugal forces of the magnet
rotating together with the composite flywheel rotor act to reduce the tensile stress in the
radial direction. As a result, designs using the hybrid composite rims have significantly
reduced the stresses in rotors and flywheels that have incorporated this design. They have
attained about twice the total storage energy of those that use either material on its own.
The present optimization approach can be used to develop an efficient flywheel rotor for
various applications. Ha also used this multiple rings structure for his thick-walled
composite ring design (1998), and his hybrid-rim-type composite flywheel rotor design
(2001).
To maximize the specific energy storage of a flywheel rotor system, Emerson and
Bakis (2000), suggested that it is desirable to maximize the top rotational speed of the
?charged? rotor. A maximum energy storage flywheel design is reached when the
stresses, due to the various loads (rotational, interference fit, thermal load), match the
rotor material strength. Today, most composite flywheel rotor designs have been based
on anisotropic elasticity equations (Lekhnitskii, 1968), and experimentally measured
material stress strengths. Most design variables are chosen from rotor geometries and
composite material properties for maximizing the storing energy. Such designs may
depend on the different application requirements, but almost always include a rotor with
one or multiple concentric rings containing hoop-oriented or radial-oriented fibers.
Arnold et al. (2002), and Portnov et al. (2005), presented their studies with multiple rims
11
or disks to meet the flywheel operation requirements.
Another successful example was presented by Gowayed and Flowers (2002).
With the use of high performance magnetic bearings, high strength composite materials,
and improved power and electronics technologies, there has evolved enabling
technologies of flywheel energy storage systems in spacecraft applications. Recent years
have seen a rapid growth of heterogeneous product fabrication techniques, many as a
result of the development of layer manufacturing processes or prototyping.
In 2002, Gowayed and Flowers presented a new approach to MDC flywheel
optimum design. The rotor is assembled from layers of pure radial or hoop
reinforcements, as shown in Figure 1.2. This approach helps to solve the problems due to
radial crack initiation and growth in filament wound flywheels and eliminates crimp that
can reduce hoop strength. It also allows the designer to increase the tip speed by safely
allowing an even distribution of radial tensile stress, significantly increasing the service
lifetime. Figure 1.3 presents the geometrical dimensions and material properties of a
typical MDC flywheel system.
Thielman and Fabien (2000) studied a similar approach. They proposed a stacked-
ply composite flywheel with alternating plies of tangential reinforcement and radial
reinforcement ply. It was shown that the performance of the flywheel could be improved
by optimizing the orientation of the fibers in the radial reinforcement ply.
Huang and Fadel (2000) demonstrated how Kumar and Dutta?s modeling
techniques could be applied to two different kinds of heterogeneous flywheel rotor
systems. One approach consists of using a finite number of distinct materials and the
other of using two or more materials with continuous volume fraction variation.
12
1.4.3 Flywheel Energy Storage System Optimal Procedure
The process of finding an optimum design for composite flywheel rotors in
general is very complex and time consuming. It involves in optimizing a large number of
parameters, including rotational speed, critical speed, dynamic mode shapes, flywheel
geometry, material characteristics, material lay-up, and stress distribution, etc.
The optimization procedures are dependent on objective functions, constraints and
model types. To judge whether a flywheel design is good, criteria or state constraints first
need to be established and results critiqued for performance satisfaction. Common
optimization design cycles of the flywheel energy storage system can be formulated as a
minimization (or maximization) problem in the following form (Ebrahimi, 1988):
Min (Max)
F(x) (Objective function)
Subject to:
I
i
(x) ? 0; i=1,?,m (Inequality constraints)
E
k
(x) = 0; k=1,?,n (equality constraints)
where:
; j=1,?,l (side constraints, x is the n-dimensional design
vector)
k
jj
i
j
xxx ??
Design variables are usually geometric parameters such as length of the shaft,
thickness of the hub, radius of the rim, material property orientations, and rotating speed
13
or even node coordinates. Limits on design variables are called side constraints, the full
set of design variables is called as the design vector.
Stress distributions, displacement deformations, temperature variations and
natural frequency changes are typical state variables. Limits placed on state variables are
termed behavior constraints and act to limit design variables response and define design
model feasibility. An optimal design problem subjected to constraints is known as a
constraints problem. If no limits are imposed, it is known as an unconstraint problem.
When a design result satisfies the constraints (limits), it is a feasible design model;
otherwise, it is an infeasible design.
Typically, weight, volume, cost or energy storage density can be the design
characteristic desired as an objective function. Of course, when all of state variables are
chosen to below the allowable limits, this design is acceptable, and all of the design
variables are simultaneously optimized to achieve the ?best? and most satisfied design.
1.4.4 Design Optimization using Finite Element Analysis
With the development of FEA technology, numerical investigations of flywheel
energy storage systems became possible and have been widely used since the early
1970?s. Muller et al., (1996) introduced some basic optimal design procedure of ANSYS
software. The model of ANSYS optimization procedure is suitable for any analysis and
simulation types with its finite element approximation. ANSYS model can accept the
predefined variables with analytical constraints. The method is well fitted to shape
optimization because ANSYS optimization model regards the geometric dimension, load
types, and the boundary conditions as parameterized design variables and objective
14
functions.
As long as the analytical sensitivities are not big and necessary problem for the
ANSYS approximate optimization model, this method can be used in a wide variety of
finite element analysis options and response quantities, such as material property non-
linearities, geometry or dimension non-linearities, large strain and displacement, user or
customer predefined constraints and dynamic load response are solvable with the ANSYS
approximate optimization method, detail reports could be found in Muller et al., (1996)?s
paper.
Hauser (1997), gave an example of using the ANSYS model to do the shape
optimization of a flywheel rotor. He used the finite element program (ANSYS) to
conduct a flywheel system. With that model, the calculation and optimization of the
flywheel were presented, and the results verified by comparing them with experimental
results. In his study, an existing flywheel system was calculated and the first suggestions
for its geometric optimization were made by the FEA. The aim of his development was to
obtain a powerful tool for the calculation and optimization of the flywheel design. This
procedure was used in part in our own study. Abdi et al. (1994), tried a shape
optimization of a complex shell under complex criteria associated with the cost, the
fatigue life and the weight of the structure. All finite element and shape optimization
computations were achieved with the ANSYS program.
1.5 Safe Operation of a Flywheel Energy Storage System
The safe operation of a flywheel energy storage system is the main bottleneck in
energy storage technology, especially if the flywheels are being considered for serious
15
space applications. Usually, the rotors are exposed to extreme loading conditions by high
centrifugal forces and partially superposed elevated temperatures due to an integrated
electrical machine, bearings or friction. A flywheel energy storage battery usually
operates at very high speed.
A detailed composite rotor failure analysis of a sophisticated fail-safe design has
been reported experimentally by Fischer et al., (2002). The tests allowed the rotor to fail
catastrophically. The different flywheel tests were made with catastrophic rotor burst, and
the significant harm to the environment is explained in their study.
As in most high speed rotating machinery, there are several flaw types that could
occur in flywheel systems, such as worn bearings and bowed/cracked shafts. Cracks are
always hidden risks in rotational machinery, especially in a rotational flywheel rotor.
Under fatigue load operating conditions, the crack may be produced as a result of the
inherent flaws or manufacturing defects in the materials. As a consequence, many
practical rotordynamic systems contain shaft/rotor elements that are highly susceptible to
transverse cross-sectional cracks caused by fatigue. For composite flywheels,
progressive damage, like delamination of composite, deboned between rim and hub, etc.,
is another common type of fault. For this reason, methods of making early detection and
diagnosis of cracks have been the subject of flywheel energy storage system
investigations. Critical questions for a cracked rotor investigation are: (1) ?How does
crack initiation and propagation affect the dynamic characteristics of a high-speed
rotating rotor system?? and (2) ?What is a good way to detect and monitor the crack
behavior in a timely fashion??.
There are two main methods used for the crack investigation of the flywheel
16
energy storage system: one is experimentation and the other is numerical analysis and
finite element simulation. These two methods can be used individually or combined.
1.6 Experimental Crack Investigations of the Flywheel Energy Storage System
Vibration measurements have been used widely in machinery condition
monitoring. In the 1990?s, several researchers reported experimental investigations of
incipient failure and dynamic behaviors. Simpson et al., (1990), utilized ultrasound to test
the energy storage flywheels that were fabricated with S2 glass-epoxy composites. The
experiment was design to study changes in the ultrasonic properties as a function of strain
history and identified possible predictors of incipient failure. Tension specimens of the
flywheel material were loaded uniaxially, and the ultrasonic properties (the shear and
longitudinal wave velocities and the attenuation) were measured as a function of strain.
The velocities were found to be excellent indicators of the maximum strain incurred at
each point of the flywheel, and the attenuation delineates the region in which the stress is
high enough to initiate micro cracking in the matrix.
Zheng et al. (1997), tested the input-output stability of the cracked rotor system. It
was suggested that fatigue crack is an important reason for the existence of low frequency
vibration components and a small fatigue crack may drive a system to instability in a very
short time. Wu and Huang (1998), evaluated crack flexibility for a shaft-disk rotor with
crack energy released. From the FFT analysis of the displacement responses, the second
harmonic component was extracted and served as a good index to detect the crack
properties. Bachschmid et al. (2000), carried a model-based diagnostic approach to
identify a transverse cracked rotor vibrations with 1X, 2X and 3X rev. components.
17
Darpe et al. (2003), investigated the steady state response of a cracked rotor
experimentally. They proposed that the response of the rotor to axial impulse excitation
could be used for diagnosis of rotor cracks, because spectral response of the cracked rotor
with and without axial excitation is found to be distinctively different. Also, in their
(2002) study paper, it was shown that the spectral responses to periodic multiple axial
impulses were present in the bending natural frequency as well as in the side bands
around impulse excitation frequency and its harmonics.
Green and Casey (2003) investigated the dynamics of a rotating shaft with
transverse cracks extensively. The second harmonic component of the system response
was shown to be the primary response characteristic resulting from the introduction of a
crack.
1.7 Numerical Investigations of the Influence of a Crack on the Vibration Flywheel
Energy Storage System
To predict accurately the response of a system to the presence of a transverse
crack, an appropriate numerical model is essential. Once the crack is included in the
system model, unique characteristics of the system response can be identified and
attributed directly to the presence of the crack. These predicted indicators then serve as
target observations for monitoring systems.
So far, many researchers have studied numerical and simulation investigations of
the dynamic behaviors of a cracked shaft and crack detection. Lee and Ng (1994), studied
the natural frequencies and mode shapes of cracked, simply supported beam based on the
Rayleigh-Ritz principle. El-Dannanh and Farghaly (1994), applied Dimarogonas? crack
18
model to compute the modal frequency parameters of stationary shafts for different crack
locations. Chen and Chen (1995), adopted a finite element model of a Timishenko beam
to study the stability of a rotating cracked shaft subjected to an end load. It was shown
that the influences of the crack are significant for the dynamic behavior of the system; as
the depth of the crack reaches certain values, it makes the forward whirling frequencies
of a shaft to become less than the backward ones.
Tsai et al., (1996, 1998), studied the shafts with a transverse open crack or multi-
crack (1997). The crack was modeled as a joint of a local spring and the transfer matrix
method was employed. Based on their work, the position of crack can be predicted by
comparing the fundamental mode shapes of the shaft with or without a crack. Takahashi
(1998), also presented an analysis for the vibration and stability of a non-uniform shaft
with a crack by use of the transfer matrix approach. Information about the effects of
cracks was provided in their paper.
Sekhar et al., (2000), discussed the effects of crack on rotor system instability.
The FEM method was used and it has been observed that the instability speed was
reduced considerably with increase in crack depth. He et al. (2001), proposed and
described a genetic algorithm? based method for shaft crack detection. Using genetic
algorithms avoids some of the weaknesses of traditional gradient-based analytical search
methods, and according to the results, it is possible to predict the shaft crack location and
configuration.
Sinou et al., (2003), presented both numerical and experimental tools to predict
the first critical speed of an engineering machine. Altstadt et al., used strain based
damage model to simulate the crack propagation. Guo et al. (2003), explored the
19
influence of one or more cracks on the natural frequencies and different models of a
rotating shaft with transverse cracks by FEM methods.
1.8 The Structure of the Dissertation
This dissertation consists of six Chapters.
A review of the literature on flywheel structure optimal design and rotor crack
detection is presented in Chapter1. Chapter 2 discusses a design optimization
methodology for multi-directional composite (MDC) flywheel hubs using Finite Element
Analysis Method. Chapter 3 presents a health monitoring methodology for MDC
flywheel systems based upon variations in imbalance eccentricity level, proposing a
relationship between the crack propagation and the imbalance eccentricity of the mass
center U
m
. Chapter 4 presents a Gabor technique as a tool for the vibration analysis and
health monitoring of cracks in rotor systems. Chapter 5 includes the investigations on the
natural frequency variations caused by the presence of the cracks in a composite disk
rotor system, made both experimentally and numerically. Chapter 6 covers the
conclusions of this study and also recommendations for future research.
1.9 Conclusions
For all of the structural optimal design and development in the previous
mentioned flywheel energy storage system studies, the materials are mostly homogenous
and only a few are made to optimize the design through selectively placing different
materials at different locations. One of major contributions of this study is the
development of a more comprehensive approach and methodology to obtain optimized
20
structures for MDC flywheel systems. The shape and material distribution were
considered and a specific optimized design for a MDC flywheel system was
demonstrated.
MDC flywheel systems to be used in spacecraft must be designed to operate
safely. The existence and propagation of cracks are one of the most common failures
found in high rotation speed flywheel systems. Understanding how cracks propagate is
also very important. The development and evaluation of a health monitoring
methodology for such systems based upon variations in imbalance eccentricity level is
another major challenge to the present research. Three post-processing software products
were developed; two are used to calculate the critical speeds of the flywheel systems
from the ANSYS FEA model, while the third is used to calculate the mass center
imbalance eccentricity due to the data from the FEA model.
The above investigations into the change of dynamics devoted to diagnosing the
cracked shafts are almost always based on a qualitative knowledge of cracked rotor
behavior. All the efforts contribute positively to detection work on cracked shafts, but
these works are still not developed enough for practical applications. Another major
contribution of this study is that developed an experimental method to detect the dynamic
behaviors of cracks using Gabor analysis techniques. A study of the application of Gabor
analysis techniques in the detection and monitoring of lateral cracks in rotor shafts has
been conducted.
A more sophisticated finite element (FEA) model should be built to further
explore the relationship between natural frequencies and the crack properties and to
verify the results from the experiment study. Following a simplified theoretical model
21
developed for the cracked disk, and the stiffness variations associated with crack depths
and locations of the disk are the investigations into a series of parametric studies using
Matlab.
22
Figure 1.1 MDC flywheel with assembled layers
23
2-D reinforcement:
? Hoop
Radial?
Axis of rotation
tresses
Hoop
stresses
Radial s
Figure 1.2 MDC flywheel stress distribution direction
24
40% S-glass in hoop
20% S-glass in radial
5
5
86 49
10
60
10
100
19
165
52% T1000 in hoop
8% T1000 in radial
T1000 hoop, vf=60%
T 1000 radial, vf=60%
S-glass hoop, vf=60%
S-glass radial, vf=60%
10
10
Figure 1.3 Geometrical dimensions of the MDC flywheel system
(Unit: mm)
25
CHAPTER 2 DESIGN OPTIMIZATION METHODOLOGY FOR MULTI-
DIRECTIONAL COMPOSITE FLYWHEEL HUBS
An optimum design has been performed to maximize the energy storing density
(ESD) of a particular type of multi-direction composite flywheel rotor with a domed-hub
attached inside the rim for an energy storage system. The flywheel rotor consists of a
steel shaft, a multiple layered domed hub and a multiple direction composite rim. The
optimal design obtained in this study shows that the ESD of the flywheel rotor can be
significantly increased by properly fabricating the hub shape using the optimal thickness t
and ratio of a/b.
2.1 Introduction
Flywheel energy storage systems have presented the potential to store very large
amounts of energy with minimum weight when compared to the traditional chemical
battery systems. However, high energy storing density (ESD) flywheel is needed to
obtain a lightweight design that is competitive with a chemical battery system. Also,
there are major concerns regarding the safety of operation and the structural integrity of
some critical components, such as the rim, hub and shaft. The stresses at high speeds (up
to 60,000 RPM) can be very high and can cause fiber breakage, micro-cracking,
debonding, and delamination, which can lead to premature failure of the flywheel design,
26
as studied by Fisher, C., et al., (1999). The consequences of failures of the flywheel
systems are a major concern that could affect the flywheel energy storing systems of a
spacecraft system or satellite. Therefore, an optimal design approach of the flywheel
system is very important in order to ensure failure-safe operation.
In order to increase the energy storing density (ESD) of the flywheel, appropriate
rotor material and reasonable structure are needed. Due to their high specific strength and
the possibility they provide for creating a load-adapted property profile, composite
materials are ideally suited for such purposes.
Compared with isotropic materials, the flywheel with the composite orthotropic
material has evident advantages, such as high stress strength, lightweight, long lifetime
and high energy storing density, as reported by Gowayed and Flowers (2002). They
developed a composite flywheel made with reinforcement in the hoop and radial
directions (Multiple Direction Composites MDC) to investigate the optimal structure
design. The rotor is assembled from layers of pure hoop or radial reinforcements. That
approach could reduce the classical problem of radial crack propagation in thick, filament
wound wheels and eliminate kinks or crimp and let the designer increase the tip speed,
thereby enhancing the specific energy of the wheel by safely allowing radial tensile stress
distributions. The multiple layer structure breaks through the restriction of the ratio of
inner radius and outer radius, which could evidently improve the distribution of the shear
stress and increase the flywheel energy storing density (ESD), as indicated by Li, W., et
al., (2001).
Usually, engineering design is an iterative process that subject to obtain a best or
optimal design. The optimal design is one that satisfies the design requirements with a
27
minimum expense of the certain factors such as cost and weight. The usual method
towards the optimal design is a traditional one. The desired objective function and
boundary condition performance of the design are first defined. Then a model is
developed with the constraints of requiring the objective function and performance
requirements. Next, an analysis of the optimal procedure is obtained and the results
evaluated by applying the design requirements. The design configuration is then created
in a way to better archive the design needs. Within the optimal design procedure, the
simulation model is built, analyzed, evaluated and recreated through a series of
improving design loops until the criteria were satisfied.
Optimum shape design of rotating disks is a concept that has attracted many
investigators. With the development of computers, a number of computational schemes
for design of rotating flywheel systems appeared in literature. Sandgren, E., et al., (1983)
developed the disk contour by a general thickness functional form. The resulting stress
analysis problem can be solved numerically. Ebrahimi, N., (1987) reported how to design
a rotating flywheel with a continuously differentiable thickness function and introduced a
defined objective function for the flywheel design problem.
Fully integrated optimization requires software features such as an analysis and
optimization database with solid and parametric modeling. Solid models are prerequisites
for efficient shape optimization so that even 3-D models can be described with the
minimum of variables. ANSYS (Version 7.0) optimization capability encompasses these
features. For the structural analysis, a four nodes axi-symmetric finite element is used to
build the numerical model. The maximum shear stress criterion for the stress analysis is
used to determine the failure for each element. For an elliptical shaped hub, the thickness
28
of the hub port t and the ratio of the long radius to the short radius, a/b, are chosen as the
design variables.
2.2 Analysis
In this study, the main objective of the flywheel energy storing system
optimization procedure is to develop and test a domed elliptical shaped hub for the
composite MDC flywheel based on the work done by Gowayed and Flowers (2002). A
multi-directional composite (MDC) flywheel rotor with a box-section hub has been
introduced, with a hub that is flexible for in-plane deformation. The hub is flexible where
needed, and rigid for out-of-plane shear load that is generated from the out-of-plane
vibration mode shapes. The goal for the hub design is to modify the shape of the hub,
increase the maximum shear stress and deformation of the hub in the limited space,
estimate the maximum energy storing density, and practice a reasonable optimal
approach for this specific problem.
As the hub lies in the heart of the flywheel structure, the shear stress distribution
for the hub was essential when the system is in high-speed rotational mode. The domed
hub connects the shaft and rim, and is required for the high tensile and compressive inter-
laminar strength is lightweight and the shape specific to protect if a crack occurred, with
maximum energy density and minimum deflection. From this study, they tried to test a
new hub design with an unusual compressive-braced domed hub. An ingenious
combination of dome hub designs, with lightweight, clear span advantage of a tensegrity
dome. So the hub design is the critical issue in this study.
29
2.2.1 Configuration
The configuration of the domed hub is shown in Figure 2.1. This figure only
shows the 2-D half domed-hub structure since the hub is symmetric about the rotational
axis. The half domed-hub includes one port and two elliptical radii. The distance from the
center coordinate to the outside of the port hub is constant, which is equal to 115.0 mm.
The thickness of the port is denoted as t. The long radius of the elliptical radius is denoted
as a, whereas the short radius is denoted as b. The calculation of the hub parameters,
such as thickness t, and the ratio of a/b, will be discussed in details.
The volume of the hub V
hub
changes with the variations of the thickness t and the
ratio of a/b. Hence, the energy storing density (ESD) varies with the design variables.
The aim of this study was to get a calculation tool for the calculation and optimization of
the domed hub-rim-shaft flywheel systems. This study has tried to identify the
remarkable advantages of the flywheel systems that can be gained by using an optimized
elliptical shaped hub in comparison to a constant thickness box-section shaped hub.
Figure 2.2 shows a 3-D configuration of the shaft, domed-hub and rim for the
flywheel system. The inner part is a steel shaft and the outer part is a multiple directional
composite (MDC) rim. A domed-hub connects the shaft and rim. The radius of the shaft
is fixed (r
s
= 19.05mm) and the inner and outer radii of the rim r
i
and r
o
are equal to
115mm and 165mm, respectively. So the volume of shaft (V
s
) and rim (V
rim
) are constant.
2.2.2 Coordinate Systems
A global cylindrical (R, ?, Z components) coordinate system shown in Figure 2.3
is introduced. R or X is the radial direction, ? or Y is the hoop direction, and Z is the
30
axial direction. The geometry items (keypoints, lines, areas, volumes, etc.) are built in the
space of global coordinate system. The components of the domed hub are modeled as a
series of blocks with different composite material properties. Since the structural aspects
of the hub are made from composite orthotropic material, a different coordinate system is
required to assign the material properties for every block of the hub. Another local shell
element coordinate system is then used to establish the orientation of material properties.
Right-handed orthogonal systems were applied to all element coordinate systems, as
illustrated in Figure 2.3.
In this study, a C++ code (listed in Appendix A) was used to create the composite
material properties, thickness and length of the each domed hub segment. The ANSYS
software was used to create the geometric model, apply the material properties, boundary
conditions and rotational speed, obtain stress distribution and produce the optimization
design, which was presented in Appendix B.
2.3 Optimization Design Approach
This study deals with the optimization of a domed-hub shape for shaft-hub-rim
flywheel system case under high rotating conditions. The optimal structure is that with
minimum weight and maximum storing energy (in other words, maximum energy storing
density ESD) and without any rotor failure behavior such as fiber breakage, debonding,
etc. The design variables are the various thicknesses of the hub flat t and the ratio of the
two-elliptical radii a/b. The formulation of the above criteria is the maximum shear stress
limit of the composite materials. All finite element and optimization computations have
been made using the ANSYS program; the shape optimization of the domed-hub has
31
been achieved using an axi-symmetric version of this structure.
2.3.1 Design Variables
Design variables represent those aspects of the design that can be varied to obtain
a maximum/minimum objective function. Design variables should be independent
variables. The thickness t of the port of the hub and the ratio a/b of two elliptical radii are
considered design variables.
2.3.2 State Variables
State variables represent the response of the design to loading and boundary
conditions, and to changes in geometry. Every state variable is a function of one or more
design variables. The maximum shear stress of the hub ?
max
(in plane ? 0.45 GPa) is
imposed to prevent the design variables from reaching unacceptable values reported by
Gowayed and Flowers (2002). Maximum stresses are compared with experimentally
measured stresses after fatigue data using a Hill polynomial criterion and employing a
factor of safety. Hill?s failure criterion is a generalization of the Von Mises-Hencky
maximum distortion energy theory and can be applied to anisotropic materials.
2.3.3 Constraint Requirements
If the hub expansion is less than that of the rim at any rotational speed level, it
results in a hub/rim separation. As a result, the hub expansion is set to be higher than that
of the rim inner radius with a constant fraction to be determined by the designer. In this
study, the maximum deformation (expansion) is greater than the preset interference
32
(5.08e-3 mm) between the hub and rim. The rim and hub can be connected with
CONTACT element through a preset interference.
As the flywheel rotor is applicable in the space ship and for the benefit of the
maximum ESD, the high rotational speed is required. In the design process, the rotating
speed of 60,000 RPM is used in all optimal designs.
2.3.4 Objective Function
The objective function is a single variable that characterizes the aspect to be
maximum/minimum. It is a function of one or more of the design variables. In this study,
the maximum energy storing density (ESD) is taken into the optimization design.
2.3.5 Governing Equations for Optimization Processing
The total stored energy (TSE) in the flywheel is given as:
2
2
1
?ITSE = , (2-1)
where:
I = the mass moment of inertia of the rotor about the spin axis.
? = the angular speed of the rotor system.
For the shaft-hub-rim flywheel system, the overall total stored energy (TSE) can
be written in terms of the inner radius r
i
, outer radius r
o
and ratio of r
i
/r
o
of each
33
component, shaft, hub, and rim:
)
)4(
))((
(
24
4
ex
rrrh
TSE
ex
o
ex
o
ex
o
?
???????
=
?
?
?
????
,
(2-2)
where:
? = density of components for the shaft, hub or rim, to be calculated from the rule of
mixtures and the knowledge of fiber and matrix densities, total and relative fiber
volume fractions (Gowayed and Flowers (2002)).
? = the rotational speed of the flywheel system.
r
o
= the outer radius of the shaft, hub or rim.
? = the ratio of the inner radius/outer radius for shaft, hub or rim.
ex = the shape factor, ex=0 means a constant section shaft, rim or hub, any other shape
can be embedded within a mathematical formulation (Gowayed and Flowers
(2002)).
h = the height of the shaft, hub or rim.
The energy storing density (ESD) is defined as the total stored energy (TSE)
devided by the total mass of the flywheel system:
34
TSE
ESD
mass
==
442
(())
()
(4 )
ex ex ex
oo o
hr r r
V
ex
?? ? ?
?
??
??? ? ? ? ?
/()?
?
??,
(2-3)
where:
V = the volume of the flywheel system, V
shaft,
V
hub
or V
rim.
2.3.6 Design Optimization Cycles
The optimization for a domed hub-shaft-rim flywheel is thus formulated as
follows:
Maximize ESD
Subjected to ?
max
? ?
crit
Find ESD, thickness t, and ratio a/b
(2-4)
In order to solve the nonlinear optimization problem formulated in equation (2-4),
a parametric model was built and optimization analysis was performed. The procedure of
optimum design, together with the structural analysis, is shown in Figure 2.4. Most
importantly, as in the procedure, the effect of the optimization function flowchart is
shown in Figure 2.5.
To start, the composite material properties and initial geometric parameters were
calculated from the C++ code, then these parameters were input to ANSYS Parametric
Design Language, using these data and ANSYS optimization operation system, a model
35
was generated parametrical, and the optimization procedure was run to test how well the
fitted ESD matched the state variables with the design variables varying. Finally, a
FORTRAN code (Appendix C) was used to calculate the maximum energy storing
energy using the output file of ANSYS as the input.
According to Muller et al., (1996)?s paper and ANSYS?s optimization design
program, the procedure is introduced as follow: For each design loop, the results from the
last loop were used as the starting point, and the new design variables are obtained from
the modified old design model. As the new deign model built, another finite element
analysis is performed. An improved approximate sub-problem is performed which allows
a new design to be obtain. This loop or procedure continues until either convergence
reaches or the problem stops. Convergence is defined to achieve when all constraints on
design variables and state variables are satisfied and the changes in all design variables
and the objective function is within the criterion between al loops.
Although the above numerical methods, which are used in the maximum of the
approximate sub-problem, can?t be considered as the real state, they present a reliable
way to reach the maximum of the energy storing density (ESD) approximation.
2.4 Results and Discussion
For this specific domed-hub optimal problem, several optimization techniques
were applied. A two-dimensional, axi-symmetric model was used in the FE, and due to
the symmetry, only one-quarter of the flywheel rotor was presented in the following
discussion. The design variables are interactive and that the one parameter that can?t be
optimized independent of the other parameters, usually starts with the primary variable;
36
then the second one, then next, until gets all parameters fit the optimal design(s).
In this work, for the sake of the minimum weight of the flywheel rotor system, the
thickness of the hub port t is the primary variable within the two design variables; as long
as the hub can support the rim and the whole rotor system; a smaller thickness value is
recommended. The first step of the hub optimal design is to find out the logical range of
the port thickness of the hub, then based on the reasonable thickness values, decide the
associated best fitting a/b variable.
2.4.1 Domed-hub Model Results
To illustrate the procedure of deciding the thickness t of the port of the hub,
Figure 2.6 shows a 2-D domed hub model built with ANSYS. According to Gowayed and
Flowers (2002) and other constraint requirements, the manufactured range of the hub port
thickness t is from 3.0 mm to 9.0 mm with 0.6mm increment change, the ratio of the two
elliptical radii a/b is from 1.95 to 2.05 with 0.1 increment change.
Starting with the thickness of hub port as t = 3.0 mm, and the ratio of two
elliptical radii a/b = 2.0 as the initial values, Figure 2.6 is the element mesh plot view.
Different colors of each segment represent different composite material properties. All
the material properties come from the C++ code [Gowayed and Flowers (2002)].
For the boundary conditions, the two ends of hub are fixed, and the rotational
speed 60,000 RPM is applied. The expansion-matching criterion of the rotor is imposed
by comparing the expansion within the hub to the expansion at the inner radius of the
rim. After running static cases, it is found that the thickness of the hub port t must be
equal to or greater than 4.8 mm to reach the constraint requirement.
37
The stress and deformation distributions of the thickness t = 4.8 mm case are
illustrated in Figures 2.7~2.10. From these plots, it is shown that the maximum
deformation lies in the hub port in x-direction, in which direction the rim is added, since
the maximum deformation (1.3685 mm) is greater than the preset interference (5.08e-3
mm) between the hub and the rim. Through the preset interference, the rim and the hub
can be connected with CONTACT element. The maximum shear stress lies in the outer
side of the domed edge; the value is 0.22 GPa, which is less than the criterion shear stress
?
max
0.45 GPa. So, in this simple hub case, the model is satisfied with the basic
requirements of optimization. It can be continued for the following optimization design to
find the maximum EDS by increasing the rotational speed till the maximum stress values
reach the constrain conditions.
2.4.2 Optimization Results of 2-D Domed Hub-Shaft-Rim Model with the Thickness of
the Hub Port at 4.8 mm and 6.0 mm
Figures 2.11 ~ 2.14 and Figures 2.15 ~ 2.17 show the results of the optimization
design when the thickness t of the hub port is 4.8 mm and 6.0 mm, respectively. After
adding the steel shaft and the composite rim, the whole model was developed. The
CONTACT elements were applied to the connection of shaft/hub and hub/rim,
respectively. The maximum shear stresses in both cases are in the inner side of the domed
hub edge. The results of the two different thickness of the hub port 4.8 mm, 6.0 mm
cycles are listed in Tables 2.1 and 2.2, respectively. In all cases, dynamical stability was
monitored as design variable changed.
From the Figures 2.11 ~ 2.14 and Table 2.1, it is shown that the shaft/hub/rim
38
flywheel model with the thickness t of the hub port equal to 4.8 mm can?t satisfy the
criterion shear stress (0.45GPa) as long as it rotating at 60,000 RPM within the a/b value
range studied. The maximum value of the shear stress for every case exceeds the criterion
shear stress. It implies the thickness t of the hub port should be increased properly in
order to meet all constrains we have.
Figures 2.15 ~ 2.17 and Table 2.2 show the optimal design results of another t
value of the hub port. When the thickness t equals to 6.0 mm, as the a/b value changes,
there are two cases that can be run in the 60,000 RPM safely while the maximum shear
stress is less than the criterion stress value. When the ratio of the elliptical radii a/b
equals to 1.99 or 2.0, the allowable tip speeds reach 61,115 RPM and 61,306 RPM,
respectively. The energy storing density (ESD) is 0.7312059 (kWhr/kg) for case a/b =
1.99, which is a little bit greater than that of case a/b = 2.0 with ESD = 0.7297543
(kWhr/kg). The two possible solutions were found for the optimal design task, and
significantly, the tip speeds and the energy storing density could exceed the design
requirements.
2.4.3 Discussion
A plot of the allowable maximum rotating speeds for each case is shown in Figure
2.18. The dotted line represents the all tip speeds results for the cases of the thickness t =
4.8 mm, the solid line represents the top speeds for the cases of the thickness t = 6.0 mm,
and the straight line represent the rotating speed of 60,000 RPM that is the requirement
speed for the shaft/hub/rim flywheel system of this study. One can see that the maximum
rotating speeds of two cases (a/b = 1.99 or 2.0 with the thickness t = 6.0 mm) are above
39
60,000 RPM.
Figure 2.19 shows the distribution of the maximum shear stress vs. the ratio of a/b
while the rotating speed is 60,000 RPM. Like the previous plot, the dotted line represents
the maximum shear stress results for all cases with the thickness t = 4.8 mm, the solid
line represents the maximum shear stress for all cases of the thickness t = 6.0 mm, and
the straight line represents the shear stress of 0.45 GPa that is the criterion stress for the
shaft/hub/rim flywheel system in this study. The shear stress of two cases (a/b = 1.99 or
2.0 with the thickness t = 6.0 mm) is less than the criterion stress, which means the fiber
will not fail at the rotating speed of 60,000 RPM.
Figure 2.20 shows the energy storing density (ESD) versus the ratio of a/b. The
max ESD is obtained at a/b = 2.01 and 1.99 when the thickness is 4.8 mm and 6.0 mm,
respectively.
Finally, the optimization design consequence is the case with the thickness of the
hub port t equals to 6.0 mm, and the ratio of the elliptical radii of a/b = 1.99. For this
case, the maximum shear stress is less than the criterion shear stress (0.45GPa) and the
ESD is 0.7312059 (kWhr/kg), which is listed in the Table 2.2.
From the above discussion, it is implied that increasing the hub port thickness and
adjusting the a/b value, the energy storing energy ESD could increase and the maximum
shear stress of the domed hub may decrease.
2.5 Conclusions
By taking the advantage of the symmetric structural, a quarter model can be
implemented. The finite element method is applied to the study by using a commercial
40
finite element code. ANSYS is employed to obtain the static and dynamic behavior of the
flywheel. The shaft, domed-hub and rim were modeled by solid brick four nodes
axisymmetric element and a special contact-target element was used to account for the
connection between the separate bodies, such as shaft and hub, hub and rim. The proper
interference was applied to each pair contact-target element. Static and transient analyses
were conducted to obtain the deformation and the stress distribution of the flywheel. By
changing the thickness of hub t and the values of a/b, a parametric study was performed
to investigate the shear stress response when the flywheel is rotating at a very high speed
(up to 60,000 RPM). The results of this study confirmed the proposed developed model
might be used for the optimization of the future design.
41
115.00 mm
Thickness t
b
a
Figure 2.1 2-D configuration of the domed-hub for the flywheel system
42
1
X
Y
Z
3d-hub,tz=4.8
APR 29 2003
10:43:56
ELEMENTS
MAT NUM
Figure 2.2 3-D configuration of the shaft/domed-hub/rim for the flywheel system
43
?, Y
Z
R
1
, X
1
R
2
, X
2
R
3
, X
3
?
1
, Y
1
?
2
, Y
2
?
3
, Y
3
Z
3
Z
2
Z
1
?
R, X
Figure 2.3 The scheme of global cylinder coordinate system and local element
coordinate systems
44
Create ANSYS log file
Build the ANSYS simulated model
Run ANSYS optimization procedure
?
max
? ?
crit
Calculate the energy storing density
Save usable or feasible design
1. a/b; 2. ?
max;
3. ESD
Yes
No
Find Max. ESD
withMin..?ma
Plot a/b vs. ?
max
Plot a/b vs. ESD
Update t,
a/b
Finish
Run the C++ code to get the
material properties and the initial
geometric parameters (t, a/b)
Figure 2.4 Optimization design procedure for a composite domed hub-shaft-rim
flywheel system
45
Figure 2.5 Flowchart of ANSYS optimization procedure
46
NOV 14 2002
21:31:46
optimization of the hubdesign a_ober_b=2
Y
Z X
1
ELEMENTS
MAT NUM
Figure 2.6 ANSYS mesh grid systems of the 2-D domed composite hub
47
.703E+09
1
NODAL SOLUTION
STEP=1
SUB =1
TIME=1
SX (AVG)
DMX =.045103
SMN =-.819E+0 9
SMX =.703E+09
MN
MX
X
Y
Z
-.819E+09
-.650E+09
-.481E+09
-.312E+09
-.143E+09
.263E+08
.195E+09
.365E+09
.534E+09
NOV 14 2002
21:40:25
Figure 2.7 S
X
stress distribution of the 2-D domed hub model
48
.220E+09
1
NODAL SOLUTION
STEP=1
SUB =1
TIME=1
SXY (AVG)
DMX =.045103
SMN =-.220E+0 9
SMX =.220E+09
MN
MX
X
Y
Z
-.220E+09
-.171E+09
-.122E+09
-.734E+08
-.245E+08
.245E+08
.734E+08
.122E+09
.171E+09
NOV 14 2002
21:41:16
Figure 2.8 S
XY
stress distribution of the 2-D domed hub model
49
.147E+09
1
NODAL SOLUTION
STEP=1
SUB =1
TIME=1
SY (AVG)
DMX =.045103
SMN =-.155E+0 9
SMX =.147E+09
MNMX
X
Y
Z
-.155E+09
-.122E+09
-.879E+08
-.544E+08
-.208E+08
.128E+08
.463E+08
.799E+08
.113E+09
NOV 14 2002
21:40:49
Figure 2.9 S
Y
stress distribution of the 2-D domed hub model
50
.013685
1
NODAL SOLUTION
STEP=1
SUB =1
TIME=1
UX (AVG)
DMX =.045103
SMN =-.117E-03
SMX =.013685
MN
MX
X
Y
Z
-.117E-03
.001417
.00295
.004484
.006017
.007551
.009084
.010618
.012151
NOV 14 2002
21:42:10
Figure 2.10 U
X
deformation of the 2-D domed hub model
51
1
DISPLACEMENT
STEP=10
SUB =2
TIME=10
DMX =.007897
MAR 11 2003
21:33:30
Y
Z X
Figure 2.11 U
X
deformation distribution of 2-D domed hub-shaft-rim model with
the thickness of the hub port t = 4.8 mm and the ratio of the elliptical radii a/b = 2.01
52
.165E+10
1
NODAL SOLUTION
STEP=10
SUB =2
TIME=10
SX (AVG)
DMX =.007897
SMN =-.644E+0 9
SMX =.165E+10
MN
MX
X
Y
Z
-.644E+09
-.389E+09
-.134E+09
.121E+09
.375E+09
.630E+09
.885E+09
.114E+10
.140E+10
MAR 11 2003
21:34:01
Figure 2.12 S
X
stress distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 4.8 mm and the ratio of the elliptical radii a/b = 2.01
53
54
1
MN
MX
X
Y
Z
optimization of the hubdesign a_ober_b=2.03256109
-.422E+09
-.323E+09
-.223E+09
-.123E+09
-.231E+08
.767E+08
.177E+09
.276E+09
.376E+09
.476E+09
MAR 11 2003
21:34:54
NODAL SOLUTION
STEP=10
SUB =2
TIME=10
SXY (AVG)
DMX =.007897
SMN =-.422E+09
SMX =.476E+09
Figure 2.13 S
XY
stress distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 4.8 mm and the ratio of the elliptical radii a/b = 2.01
1
MNMX
X
Y
Z
optimization of the hubdesign a_ober_b=2.03256109
-.192E+10
-.159E+10
-.126E+10
-.923E+09
-.588E+09
-.254E+09
.797E+08
.414E+09
.748E+09
.108E+10
MAR 11 2003
21:34:26
NODAL SOLUTION
STEP=10
SUB =2
TIME=10
SY (AVG)
DMX =.007897
SMN =-.192E+10
SMX =.108E+10
Figure 2.14 S
Y
stress distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 4.8 mm and the ratio of the elliptical radii a/b = 2.01
55
1
X
Y
Z
optimization of the hubdesign a_ober_b=2
MAR 27 2003
11:07:47
DISPLACEMENT
STEP=10
SUB =2
TIME=10
DMX =.006702
1.99
Figure 2.15 U
SUM
deformation distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 6.0 mm and the ratio of the elliptical radii a/b = 1.99
56
1
MN
MX
X
Y
Z
optimization of the hubdesign a_ober_b=2
0
.745E-03
.001489
.002234
.002979
.003724
.004468
.005213
.005958
.006702
MAR 27 2003
11:08:19
NODAL SOLUTION
STEP=10
SUB =2
TIME=10
USUM (AVG)
DMX =.006702
SMX =.006702
1.99
Figure 2.16 U
X
deformation distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 6.0 mm and the ratio of the elliptical radii a/b = 1.99
57
58
1
MN
MX
X
Y
Z
optimization of the hubdesign a_ober_b=2
-.377E+09
-.286E+09
-.195E+09
-.104E+09
-.137E+08
.770E+08
.168E+09
.258E+09
.349E+09
.440E+09
MAR 27 2003
11:07:15
NODAL SOLUTION
STEP=10
SUB =2
TIME=10
SXY (AVG)
DMX =.006702
SMN =-.377E+09
SMX =.440E+09
1.99
Figure 2.17 S
XY
stress distribution of 2-D domed hub-shaft-rim model with the
thickness of the hub port t = 6.0 mm and the ratio of the elliptical radii a/b = 1.99
59
Maximum Rotating Speed
54000
55000
56000
57000
58000
59000
60000
61000
62000
1.98 1.99 2.00 2.01 2.02 2.03 2.04 2.05
a/b
R
o
ta
tin
g
S
p
e
e
d
(R
P
M
)
__thickness =6.0 mm- - -thickness =4.8 mm
Figure 2.18 Distribution of the allowable tip speed vs. the ratio of a/b
(The maximum shear stress is equal to the criterion shear stress for the flywheel rotor)
Maximum Shear Stress
4.0E+08
4.2E+08
4.4E+08
4.6E+08
4.8E+08
5.0E+08
5.2E+08
5.4E+08
5.6E+08
5.8E+08
6.0E+08
1.98 1.99 2.00 2.01 2.02 2.03 2.04 2.05
a/b
Shear Stress (GPa)
- - -thickness = 4.8 mm __thickness =6.0 mm
Figure 2.19 Distribution of the maximum shear stress vs. the ratio of a/b
(The rotating speed ?= 60,000 RPM)
60
Maximum ESD
0.40
0.45
0.50
0.55
0.60
0.65
0.70
0.75
1.98 2.00 2.02 2.04
a/b
ESD (KWh/Kg)
- - - -thickness = 4.8 mm
__thickness =6.0 mm
Figure 2.20 Distribution of the maximum ESD vs. the ratio of a/b
61
Table 2.1 Optimization of the flywheel rotor with the thickness of the hub port
t = 4.8 mm (rotating speed ?= 60,000 RPM)
(a/b) ESD
(KWhrs/Kg)
V
tot
(m
3
) ?
max
(GPa)
1.98 0.607058 0.003294787 0.524e9
1.99 0.5998307 0.003303771 0.528e9
2.0 0.6699522 0.003210075 0.480e9
2.01 0.7071015 0.003151479 0.462e9
2.02 0.6947668 0.003174478 0.467e9
2.03 0.6746069 0.003208308 0.476e9
2.04 0.5735387 0.003382578 0.588e9
2.05 0.6726946 0.003223103 0.477e9
62
Table 2.2 Optimization of the flywheel rotor with the thickness of the hub port
t = 6.0 mm (rotating speed ?= 60,000 RPM)
(a/b) ESD
(KWhrs/Kg)
V
tot
(m
3
)
?
max
(GPa)
1.98 0.6730617 0.00345878 0.462e9
1.99 0.7312059 0.003372043 0.434e9
2.0 0.7297543 0.003416068 0.440e9
2.01 0.5517544 0.003767504 0.543e9
2.02 0.604307 0.003620698 0.502e9
2.03 0.6890825 0.003472563 0.452e9
2.04 0.6563657 0.003547441 0.469e9
2.05 0.6193641 0.003616594 0.466e9
63
CHAPTER 3 HEALTH MONITORING METHODOLOGY FOR MDC FLYWHEEL
SYSTEMS BASED UPON VARIATIONS IN IMBALANCE ECCENTRICITY LEVEL
MDC flywheel systems to be used in spacecraft must be designed to operate
safely. A certain level of safety should be assured with a sufficiently high probability.
The occurrence of cracks is one of the most common failures found in high rotation speed
flywheel systems. Understanding what causes cracks to occur and how they propagate is
also very important, as well as static analysis and modal dynamic analysis. Fractures of
composite laminates and multi-layer structures have been widely studied by Chudnovsky,
A. and Li, R. S., (1995). The analysis of the crack tip field is a fundamental topic that has
received significant attention.
An investigation was undertaken to study the damage initiation behavior in a
layered structure composite rim subjected to a radial or hoop direction crack in different
layers. The focus of this research effort was to develop an understanding of the crack
propagation behavior and the calculation procedures of the mass center imbalance
eccentricity related to the existence of cracks.
3.1 Computational Procedure
3-D models for an MDC layered flywheel with cracks in various directions were
developed. The displacement at every node was computed via a series of steady-state
64
cases.
The crack is closed when the speed of flywheel is zero, but starts to open when
the flywheel starts to rotate. In order to simulate the crack in the FEA model, special
treatment is introduced for the nodes within the crack area. There are two nodes at the
points on the crack and these two nodes can separate when the flywheel is rotating.
Figures 3.1 and 3.2 illustrate hoop-direction and radial-direction cracks, respectively.
Figure 3.1 (a) is the initial section view for the hoop-direction crack area. The crack area
is marked with a small rectangle, showing where the crack was created. Within this area
there are two grid points superimposed at each of the points marked with a filled circle.
In order to simulate a hoop direction crack, these two points are not connected to each
other. One grid point belongs to the left cell, while the other grid point belongs to the
right cell. Thus, these two points can move in opposite directions when the flywheel
starts to rotate at high speed. This is shown in Figure 3.1 (b). Similarly, the illustration of
the occurrence of radial direction crack is shown on Figures 3.2 (a) and (b). Again,
initially there are actually two grid points superimposed at each point in the crack area
marked, one of which belongs to the top cell and the other to the bottom cell.
As mentioned in Chapter 2, the flywheel has a sandwich-like structure. One layer
is made of hoop-reinforcement composite material and the other layer is made of radial-
reinforcement composite materials. Thus, the crack can grow in two directions: the hoop
direction and the radial direction, as shown in Figure 3.3. This produces four different
types of crack:
1) Radial-direction crack that develops in a hoop reinforcement material layer. This is
refereed to (hoop-radial) crack, where the first word represents the attribution of the
65
layer, and the second word indicates the direction of the crack.
2) Radial-direction crack develops in a radial reinforcement material layer (radial-
radial);
3) Hoop-direction crack that develops in a hoop reinforcement material layer (hoop-
hoop);
4) Hoop-direction crack that develops in a radial reinforcement material layer (radial-
hoop);
Figures 3.4 and 3.5 present the 3-D layered structure model developed for crack
analysis. In order to simplify the computation, only the rim was considered and the shaft
and hub were replaced by 32 evenly-spaced springs which were connected with the rim.
The mass of these springs was taken to be negligible, and the spring constants of these
springs were assumed to be very large.
The rim was made of two different kinds of composite materials. The material
properties are listed in Table 3.1. Different colors represent different materials. The
STATIC mode was set for all the cases and the mesh was generated by a sweep method.
The computation domain was divided into 48 parts evenly in the hoop direction. A
20
z
?19
x
grid system was adopted in every section. Since the number of grid points
increased for the 3-D simulation, the computational time also increased dramatically.
The boundary conditions were as follows: At the center of every spring, a DOF
Ux, Uy, Uz=0 constraint is applied. Where the end of each spring was connected with
the rim, Uz=0 constraints were applied. The rotational speed was varied from 10,000
RPM to 60,000 RPM.
66
3.2 The Calculation of Imbalance Eccentricity of the Mass Center (U
m
)
Due to the imperfect geometry and the occurrence of the crack, the rotating
flywheel generates a mass imbalance force. This force is a function of the square of the
rotating speed, and therefore grows rapidly with large increases in speed. This makes
balancing the rotor a critical element in the development of the flywheel. The imbalance
eccentricity of the mass center calculation is one of the most important steps in the entire
rotor balancing procedure. In this chapter, it is assumed that the material is consistent and
the shape is perfect. Any mass imbalance force present is produced only by the
occurrence of the crack. The numerical simulation results for a flywheel with a crack
provide the displacement at every grid point. Based on this information, the imbalance
eccentricity
Um
of the center of mass can be calculated.
Since FEA analysis only provides the displacement at every grid point, the
equivalent mass of each node cell has to be calculated in order to get the imbalance
eccentricity of the center of mass. Usually the mass center of every cell is located in the
center of the cell. The best way to weight the equivalent mass for every grid point
depends on the weighting criteria. The calculation of the imbalance eccentricity in the
center of mass involves the first moment of the inertia about the rotational axis. The
weighting criteria are defined as follows:
i
m
The first moment of inertia should remain consistent regardless of whether the
center of mass is assumed to be located in the center of the cell or the mass of the cell is
divided between all its grid points. Notice that the first moment of inertia is defined as
follows:
67
?
?= dmrS . (3-1)
It is necessary to introduce the node mass distribution factor (NMDF) for every
node. The equivalent mass at each node from one cell can be calculated from the product
of the NMDF and the total mass of the cell. The NMDF is defined as the ratio of the
equivalent mass of the node and the total mass of the cell. The derivation of the NMDF
can be obtained through the weighting criteria mentioned above.
Referring to Figure 3.6, the mass m? of the small area between inner radius
and the random radius
i
r
r is given as:
hdrrm ????=? ??2 . (3-2)
The total mass of the area between two concentric circles and is given as:
t
M
i
r
o
r
hrrM
iot
????= ?? )(
22
. (3-3)
It is assumed that the total mass of the cell is centralized on two nodes, node 1 and
node 2. The ratio of mass centralized on node 1 and the total mass, NMDF of node 1, is x
whereas the ratio of mass centered on node 2 and the total mass, NMDF of node 2, is 1-x.
Thus we have:
68
hrrxrrxhrrM
ioioiot
??????+??=????= ????? )]()1()([)(
2222
22
. (3-4)
The first moment of inertia of the area between the two centric circles is:
?
?????=
o
i
r
r
rhdrrS )2(
1
?? . (3-5)
The first moment of inertia, assuming the mass is only centered at node 1 and node 2, is
given as:
hrrrxrrrxS
oioiio
???????+???= ??? ])()1()([
222
2
2
. (3-6)
The weighting criterion requires:
21
SS = . (3-7)
So we have:
hrrrxrrrxrhdrr
oioi
r
r
io
o
i
???????+???=?????
?
????? ])()1()([)2(
2222
. (3-8)
Eliminating constants? , , h? ? from both sides of the equation, we have:
69
))1()((2
222
oi
r
r
io
rxrxrrdrr
o
i
??+??=
?
. (3-9)
Evaluating the integration on the left-hand side of Equation (3-9), we get:
))()(()(
3
2
33
ooiioioio
xrrxrrrrrrr ?++?=? . (3-10)
Solving for x, we get:
22
22
3
1
3
1
3
2
oi
oioi
rr
rrrr
x
?
??
= , (3-11)
where , . Thus, given the radii of the inner edge and the outer edge of the area,
the NMDF of node 1 (at the inner edge) can be calculated from Equation (3-11) and the
NMDF of node 2 (at the outer edge) is simply equal to 1-x.
oi
rr < 0>x
With the aid of the above process, the mass that is centered at the cell center can
be transformed into a series of masses centered at every grid point, consist with the first
moment of inertia. Combined with the displacement at every grid point, the imbalance
eccentricity of the center of mass thus can be calculated.
The original coordinates of grid point i are , and the displacements of grid
point i are , when the flywheel rotating at high speed. Then the new coordinates,
, , are given as:
i
x
i
y
i
ux
i
uy
'
i
x
'
i
y
70
',
'
,
.
ii i
ii i
x xux
yyuy
=+
=+
(3-12)
The first moment of inertia of the flywheel when it is stationary is given as:
,
.
xi
yi
Sm
Sm
=
=
?
? i
x
y
'
i
i
x
y
(3-13)
If we assume that the geometry of the flywheel is perfect, we have:
0,
0.
xi
yii
Smx
Smy
==
==
?
?
(3-14)
It is obvious that the cells will deform if the flywheel is rotating at a high speed
and the geometry of the flywheel will no longer be perfect. We assume that the
equivalent mass at every grid point is not going to change with the deformation of the
cell. Thus, the first moment of inertia of the deformed flywheel is given as:
''
'
,
.
x
yi
Sm
Sm
=
=
?
?
(3-15)
Substituting Equation (3-12) and Equation (3-14) into Equation (3-15), we have:
71
72
x
uy
''
''
,
.
xiii i i
yiii i i
Smxmxmuxmu
Smymymuym
==+=
==+=
????
????
(3-16)
The total mass of the flywheel can then be calculated as:
?
=
iT
mM . (3-17)
Thus, the imbalance eccentricity of the center of mass in the x and y directions are given
as:
t
x
mc
M
S
x
'
= , (3-18)
T
y
mc
M
S
y
'
= , (3-19)
where : the imbalance eccentricity of the center of mass of the rim in the X direction,
mc
x
mc
y : the imbalance eccentricity of the center of mass of the rim in the Y direction.
Thus, the overall imbalance eccentricity of the center of mass is given as:
22
mcmcm
yxU += . (3-20)
A FORTRAN code based on the procedure above is presented in Appendix D.
this can be used to calculate the imbalance eccentricity of the center of mass using the
output file of ANSYS as the input.
3.3 Results and Discussion
3.3.1 Relationship between U
m
and Rotation Speed
Figure 3.7 presents the relationship between U
m
and rotation speeds. For all the
cases, the distance between the middle of a hoop crack and the axis of rotation
equals . The distance between a radial crack center and the axis of rotation
also equals .
0.2/)(
inout
rr +
0.2/)(
inout
rr +
Rotational speed can have a significant effect on the crack behavior. A series of
cases were tested for four different of cracks. The rotation speeds are varied from 10,000
RPM to 60,000 RPM. For every case, the imbalance eccentricity in the center of mass
was calculated using the new post-processing software developed. The results are
summarized in Table 3.2 and plotted in Figure 3.7. Referring to Figure 3.7, the imbalance
eccentricity of the center of mass U
m
grows as the speed of rotation increases. The
imbalance eccentricity of a flywheel with a crack developing in the hoop reinforcement
material is larger than that of a flywheel with a crack developing in the radial
reinforcement material, regardless of the direction the crack is developing in, both in the
hoop and radial directions. Notice that the layer of hoop reinforcement material is thicker
than the layer of radial reinforcement material. This implies the multi-layered structure
can be used to successfully impede the propagation of the crack. For the direction in
which the crack is developing, the imbalance eccentricity of a flywheel with a crack
developing in the radial direction is greater than with a crack developing in the hoop
73
direction, even if the length of the crack developing in the hoop direction is greater than
the length of the crack developing in the radial direction in our models. This implies that
a crack developing in the radial direction generates a bigger imbalance mass force than
does a crack developing in the hoop direction. This phenomenon can be explained as
follows:
For a crack developing in the hoop direction, the displacement of the points
between the crack and the outer radius increases the mass imbalance eccentricity.
However, the displacement of the points between the crack and inner radius decreases the
mass imbalance. For a crack developing in the radial direction, the displacements of the
points on both sides of the crack increase the mass imbalance eccentricity. Thus, a crack
developing in the radial direction is more harmful than a similar crack developing in the
hoop direction.
3.3.2 Relationship between U
m
and Various Material Properties
Figure 3.8 presents the relationship between imbalance eccentricity U
m
and the
ratio of modulus of elasticity values. Material properties are another factor that can affect
imbalance, so a series of cases with different ratios of the modulus of elasticity value of
the material to the modulus of elasticity value of a standard material were tested for radial
and hoop direction cracks in the hoop reinforced layer. The ratios of the modulus of
elasticity values were 10, 8, 5, 2, 0.5, 0.2, 0.125, and 0.1. The rotational speed was
60,000 RPM. For every case, the imbalance eccentricity U
m
in the center of mass was
calculated using the new post-processing software that had been developed. The results
are summarized in Table 3.3 and plotted in Figure 3.8. Referring Figure 3.8, the flywheel
74
made of a material with a small value for the modulus of elasticity produced a greater
imbalance U
m
, whereas the flywheel made of a material with a large value for the
modulus of elasticity produced a smaller imbalance eccentricity U
m
. The linearity
exhibited in Figure 3.8 implies that the relationship between the imbalance eccentricity
U
m
and the ratio of the modulus of elasticity and the standard modulus of elasticity is an
exponential function.
3.3.3 Relationship between U
m
and Crack Size
A series of cases modeling a flywheel with different lengths of cracks developing
in the hoop direction were tested. Table 3.4 lists the results showing the amount of
imbalance eccentricity U
m
observed for various sizes of cracks. Figure 3.9 shows the
relationship between imbalance eccentricity U
m
and crack length. The cracks were
located in the hoop-reinforced layers (hoop-hoop). The rotation speed of the flywheel
was 60,000 RPM. As the crack length increased, the imbalance eccentricity U
m
increased
when the crack was shorter than half of the circumference, then Um decreased when the
crack was longer than half of the circumference. The imbalance eccentricity Um, reached
a maximum value when the crack was half as long as the circumference. The curve
presented in Figure 3.9 also exhibits a symmetric about l/c =0.5.
3.3.4 Relationship between U
m
and ?
? is defined as the ratio of the inner diameter to the outer diameter. If the outer
radius remains constant then ? increases as the inner radius increases. Only the flywheels
75
with a crack developing in the radial direction were considered. The cracks were in the
hoop reinforcement material layer. The ? values were 0.2, 0.4 and 0.6. The height of the
flywheel remains constant. The effect of the crack length was also investigated. The
distance between the middle of the crack and the rotational axis is . Table
3.5 lists the imbalance eccentricity for every case. Figure 3.10 shows the relationship
between the imbalance eccentricity U
0.2/)(
inout
rr +
m
and the ratio of the crack length and the
difference between the outer radius and the inner radius. Figure 3.10 shows that with the
same ? value, the imbalance eccentricity U
m
also increases as the ratio of the crack length
and the difference between the outer radius and the inner radius increases. With the same
length of crack, the imbalance eccentricity U
m
increases as ? increases. This implies that
a flywheel with small ? can reduce the imbalance eccentricity in the center of mass. On
the other hand, a flywheel with small ? also lowers the storage efficiency. The design of a
flywheel in a spacecraft usually requires high-energy storage density. Thus, a flywheel
with small ? is not a good choice for a flywheel to be used in a spacecraft.
3.4 Conclusions
The assumptions used in the calculations and the principles governing the
imbalance eccentricity in the center of mass were presented in this chapter. Post-
processing software for the calculation of imbalance eccentricity U
m
based on the output
of ANSYS was developed and used to study the relationship between the imbalance
eccentricity of the center of mass U
m
and rotational speed, variations in the material
properties, variations in the crack size and variations in ?. The results can be
76
summarized as follows:
(1) The imbalance eccentricity of the center of mass U
m
increases as the speed of
rotation increases.
(2) A crack developing in the radial direction is more harmful for a flywheel.
(3) For a flywheel with a crack developing in the hoop direction, as the crack
length increases, the imbalance eccentricity U
m
increases up to the point
where the crack is shorter than half the circumference, and then decreases
once the crack is longer than half the circumference. It reaches its maximum
value when the crack length equals half the circumference.
(4) A flywheel with small ? can reduce the imbalance eccentricity in the center of
mass, but this lowers the energy storage efficiency.
77
(a)
(b)
Figure 3.1 Illustration of hoop-direction crack.
78
(a)
(b)
Figure 3.2 Illustration of radial-direction crack.
79
(b) Radial direction crack
(a) Hoop direction crack
Figure 3.3 Top view of crack model
80
Figure 3.4 3-D model
Figure 3.5 3-D model
81
dr
r ri
Node 1: NMDF(x)
Node 2: NMDF(1-x)
ro
h
Figure 3.6 Configuration of nodal mass distribution
82
RotatingSpeed(rpm)
U
nba
l
a
nc
e
o
f
M
a
s
s
C
e
n
t
e
r
U
m
(
?
m)
0 10000 20000 30000 40000 50000 60000 70000
10
-4
10
-3
10
-2
10
-1
10
0
10
1
10
2
Hoop-Rad
Rad-Rad
Hoop-Hoop
Rad-Hoop
Figure 3.7 Relationship between U
m
and rotation speed
83
E/E
s
U
m
(
?
m)
10
-1
10
0
10
1
10
-1
10
0
10
1
10
2
10
3
Hoop-Rad
Hoop-Hoop
Figure 3.8 Relationship between U
m
and different material properties
(Rotation speed: ? = 60,000 RPM)
84
l/c
U
m
(
?
m)
0.20.40.60.8 1
2
4
6
8
10
12
14
16
18
Figure 3.9 Relationship between U
m
and crack length
(l: length of crack, c: circumference, rotation speed: ? = 60,000 RPM)
85
Crack Length (%)
U
m
(
?
m)
0 25507510
10
-3
10
-2
10
-1
10
0
10
1
?=0.6
025507
?=0.2
?=0.4
Figure 3.10 Relationship between U
m
and ? in hoop layer radial crack
(Rotation speed: ? = 60,000 RPM)
86
Table 3.1 Material properties of the 3-D rim model
Young?s Modulus Density Poisson?s Ratio
PaeE
x
1217802.0= 288013.0 ?= eU
xy
PaeE
y
1060262.0= 42648.0=
yz
U
Rim radial
layer
PaeE
z
1060262.0=
3
1520.0 /kg m? =
288013.0 ?= eU
xz
PaeE
x
1060262.0= 42648.0=
xy
U
PaeE
y
1060262.0= 26000.0=
yz
U
Rim hoop
layer
PaeE
z
1217802.0=
3
1520.0 /kg m? =
26000.0=
xz
U
87
Table 3.2 The U
m
at different rotational speeds in different layers
Rotation speed U
m
(?m)
(RPM) Hoop-radial Radial-radial Hoop-hoop Radial-hoop
10,000.00 7.752445E-01 3.079648E-03 1.763526E-01 4.335275E-04
20,000.00 3.402593E-00 1.358098E-02 7.721130E-01 2.388601E-03
30,000.00 7.654059E-00 3.137149E-02 1.740034E-00 4.425180E-03
40,000.00 1.496168E+01 6.014105E-02 3.395000E-00 5.964242E-03
50,000.00 2.125401E+01 8.624395E-02 4.817300E-00 1.059468E-02
60,000.00 3.060827E+01 1.229544E-01 6.948030E-00 1.946457E-02
88
Table 3.3 The U
m
due to different material properties
(Rotation speed: ? = 60,000 RPM)
U
m
(?m)
Ratio of modulus
of elasticity
(E/Es)
Hoop-radial Hoop-hoop
10.0 2.729340E-00 6.667237E-01
8.0 3.423327E-00 8.124267E-01
5.0 6.245662E-00 1.360854E-00
2.0 1.438898E+01 3.545602E-00
1.0 3.060827E+01 6.948030E-00
? 6.824654E+01 1.443283E+01
1/5 2.178154E+02 3.255985E+01
1/8 4.111557E+02 5.537692E+01
1/10 5.588994E+02 6.419068E+01
89
Table 3.4 The U
m
observed for different crack sizes of hoop layer hoop
cracks (Rotation speed: ? = 60,000 RPM)
Crack Nodes Crack Sizes
U
m
(m) (?m)
1-node 0.015118914 1.469620E-00
3-node 0.045356743 3.501807E-00
5-node 0.075594572 5.505914E-00
7-node 0.105832401 7.387419E-00
9-node 0.136070229 9.195531E-00
11-node 0.166308058 1.080063E+01
13-node 0.196545887 1.229428E+01
15-node 0.226783716 1.339638E+01
17-node 0.257021545 1.446539E+01
19-node 0.287259373 1.534829E+01
21-node 0.317497202 1.570897E+01
23-node 0.347735031 1.597986E+01
25-node 0.37797286 1.603758E+01
27-node 0.408210688 1.567029E+01
29-node 0.438448517 1.507164E+01
31-node 0.468686346 1.425701E+01
33-node 0.498924175 1.320116E+01
35-node 0.529162004 1.187100E+01
37-node 0.559399832 1.048380E+01
39-node 0.589637661 8.588857E-00
41-node 0.61987549 6.841281E-00
43-node 0.650113319 4.903017E-00
48-node 0.725707891 8.009648E-03
90
Table 3.5 The U
m
due to different values of ? for a hoop layer radial crack
(Rotation speed: ? = 60,000 RPM)
Crack
ratio (%)
U
m
(?m)
?=0.6
Crack
ratio
(%)
U
m
(?m)
?=0.2
Crack ratio
(%)
U
m
(?m)
?=0.4
0 9.869918E-03 0 2.755247E-03 0 5.642903E-03
16.67 4.613105E-02 8.33 5.113679E-03 22.22 2.397524E-02
33.33 2.647962E-01 25 1.608824E-02 44.44 1.327298E-01
50 7.711860E-01 41.67 3.970804E-02 55.56 3.097230E-01
66.67 1.941827E-00 58.33 8.969170E-02 66.67 5.963567E-01
83.33 3.882489E-00 75 2.982346E-01 77.78 9.657845E-01
100.0 5.848550E-00 91.67 6.583967E-01 88.89 1.490675E-00
91
CHAPTER 4 GABOR TECHNIQUES AS A TOOL FOR THE VIBRATION
ANALYSIS AND HEALTH MONITORING OF CRACKS IN ROTOR SYSTEMS
4.1 Introduction
The dynamics and diagnostics of cracked rotor have been gaining importance in
recent years. Although rotors can be carefully designed for fatigue loading and high level
of safety by using high-quality composite materials and precise manufacturing
techniques, failures of rotors as a result of cracks occur often, particularly in high-speed
rotating machines. Cracks often initiate in or near the regions with high local stresses.
When stresses increase to some threshold level, the leading edge of the crack (crack tip)
starts to propagate. Cracks cause localized flexibility in the structure and therefore affect
the dynamic behavior to some extent. There are two basic kinds of cracks: transverse
cracks and torsion cracks due to different stress types. In most rotor systems, bending
stress tends to be dominant. Accordingly, the present study considers only transverse
cracks.
With the growth of a crack, the bending stiffness of the shaft decreases, the
unbalance level changes, producing changes in both the amplitudes and the phases of the
1X (first fundamental harmonic component) and higher harmonic components. The
natural frequencies are also affected by the change of the stiffness, depending on the rotor
mode shape and the location of the crack. Some modes are more strongly affected than
92
others.
Harmonic analysis is one of the most popular and useful analysis methods or tools
for engineers and scientists. First, what is the meaning of ?harmonic?? Usually when one
uses the term ?harmonic?, one means that the measured frequency values are integer (or
fractional) multiples of a fundamental frequency or rotating speed, as explained by Shao
H., et al., (2003). By measuring and analyzing of the different harmonic components of
the signal sets (plotting or calculating the amplitudes and phases of each component),
engineers can determine whether or not the rotor system is operating normally. Typically,
the traditional harmonic analysis methods, such as the Fourier transform, are used for
constant (or near constant) rotor speed operation, so that the fundamental frequency does
not change appreciably with time. However, when the rotor speed is increasing or
decreasing significantly, the harmonic components may overlap. For such cases, the
traditional power spectrum methods are not able to distinguish and extract the vibration
signals caused by differing excitation sources.
In the present study, Gabor analysis, or a joint time-frequency method, is
evaluated as a signal analysis technique for crack detection and health monitoring. Gabor
analysis is a special set of time-varying harmonic analysis. From the concept of
expansion and series, Dennis Gabor (1900-1979), introduced the idea of that expanding a
signal into a set of weighted frequency modulated Gaussian functions, the details of this
method could be found in Qian (2002).
Historically, accomplishing crack detection using conventional analysis methods
is a very challenging task. Fourier methods are generally often not sufficiently sensitive
to allow for observing subtle changes in primary frequency components; so many
93
researchers have modeled and studied the dynamic response of cracked rotors using
different approaches. Bicego, V., et al., (1999), Bertholet, Y., Qu, J. and Jacobs, L., et al.,
(2000), Harmon, L. and Baaklini, G. (2000), used ultrasonic equipment to detect the
crack effects. Kim and Kim (2003), used the Newmark time integration method and the
fast Fourier transform (FFT) to study the effects of various parameters, such as the crack
depths and crack positions. Adewusi and Al-Bedoor (2001), presented the dynamic
response of an overhang rotor with a transverse crack using the discrete wavelet
transform (DWT) analysis technique. Shekhar and Prabhu (1998), studied the transient
response of a cracked rotor passing through its critical speed. In their extended study
(2001) and (2003), they reported that the continuous wavelet transform (CWT) could be
applied to extract feathers of transverse cracks from the time-domain signals of rotor-
bearing systems. They concluded that crack propagating produces changes in vibration
amplitudes of frequency scale levels corresponding to 1X, 2X and 4X harmonics. Gentile
and Messina (2003), also used the continuous wavelet transform (CWT) to detect cracks.
Green, I., et al., (2003), developed an asymmetry crack model to show that a 2X
harmonic component is a primary response characteristic resulting from the introduction
of a crack.
The archived literature also includes applications of the Gabor transform for
detecting crack characteristics. Liu, G., et al., (1998), studied crack edge-reproduced
wave. In their work, crack edge-reproduced wave (ultrasonic signals) had been detected
by using the Gabor transform of time-frequency analysis. Quek, S.K., et al., (2001),
examined the sensitivity of wavelet technique in the detection of cracks in beam
structures. The results show that the Gabor wavelet transform is a useful tool in detection
94
of cracks in beam structures. Li, Zheng and Sia, S., et al., (2004), considered the
continuous Gabor wavelet transform to be applied to analyzing flexural waves in a
cantilever beam with an edge crack. The results show that it is an effective damage
detection method and can be extended to detect the damage of complicated structures.
Peng, G., et al., (2004), presented an active monitoring method, based on Lamb wave and
wavelet transform, to determine damage locations. It is concluded that the wavelet
transform, using the Gabor wavelet, effectively decomposes the differential signal into its
time-frequency components, and the peaks of the time-frequency distribution near the
center frequency of the exciting signal indicate the arrival times of waves. By the
calculation of time delays between the arrival of the differential and exciting signals, the
damage localization points can be obtained. However, the author did not find in the
literature that Gabor analysis techniques have been used for rotating shafts or flywheel
disk systems.
In this work, Gabor analysis is then used to examine the characteristics of the
respective frequency components over time in order to develop a methodology for the
identification of crack properties based upon this information. First, a discussion of the
basics of Gabor analysis is presented. Experimental data obtained from a laboratory test
rig is used for this purpose. The characteristics of a vibration signal from a rotating shaft
with a transverse crack are examined. From the study, one can see that the Gabor analysis
is not only more robust and easier to implement, but also offers more insight into the
physical process.
95
4.2 Gabor Analysis
4.2.1 FFT Drawbacks
The Fourier transform (FFT) is the most widely used signal processing method.
However, it is only strictly applicable for stationary signals. For signals involving run-up
or run-down, the FFT-based spectrum cannot precisely distinguish the frequency
components.
There are three main drawbacks of the Fourier transform. First, FFT may lose
time information. For example, the basis functions of the Fourier transform, pure sinusoid
functions, last from negative infinite to positive infinite in time. Thus, for a specific
frequency, the spectrum describes the total energy of this frequency in the time domain.
Also, the spectrum cannot indicate how these frequencies evolve over time. Figure 4.1
shows that the two different time domain signals (a) linear chip signals; (b) reversed
linear chip signals; get the same spectrum distribution in Figure 4.1(c) plot by using FFT
transform.
Second, the traditional FFT-based spectrum can only be used to analyze stationary
signals and cannot track the change of frequency. For example, in the spectrum of the
chirp signal in Figure 4.1, the frequency bandwidth of signals is proportional to the rate
of the frequency change. In Figure 4.2, the spectrums of the two-chirp signals overlap,
while the Gabor transform can distinguish them.
Third, another problem of the Fourier transform is that the length of the data
limits the resolution of the FFT. The frequency resolution of FFT is defined as:
96
dataofnumber
spanfrequency
f =?
.
(4-1)
As a result, long time signals have to be obtained for better resolution.
4.2.2 Discrete Gabor Expansion
In the present study, Gabor analysis and the modal-based spectrum are used to
overcome these drawbacks.
For a given discrete time sequence , the corresponding discrete Gabor
transform can be computed with the basis functions, which are the modulated and shifted
version of
][kx
][k? :
?
?
=
?
?=
1
0
2
,
][*][
~
L
k
N
nkj
nm
emTkkxC
?
? , (4-2)
where
C
m,n
= the Gabor coefficients
T = the discrete time sampling interval
N = the total number of frequency lines
][
~
kx is the preprocessed data of . The reason for the use of ][kx ][
~
kx originates
from the artificial jumps at the beginning and at the end of the input signal, because the
analysis window ?[k] oversteps the boundaries and causes some information to be missed.
97
An empirical solution is to extend the data samples at the beginning and the end
according to the length of the analysis window. Although there are several methods to
achieve this, such as zero padding, symmetric, periodic and spline extensions, our results
show that symmetric extension and spline extension are more suitable because the signals
extended by the two methods are smooth at the boundaries.
In Equation (4-3), since a basis function is concentrated in a region on the time-
plane, the Gabor coefficients can describe how the frequency contents change in the time
domain. The Gabor spectrogram can thus be calculated based on the Gabor coefficients
as
],[],[
'
,','
''
,
kiWVDCCkiGS
hhnm
Dnnmm
nmD ?
??+?
= , (4-3)
where WVD
h,h'
[i,k] is the cross Wigner-Ville distribution of the frequency-modulated
Gaussian functions. The order of the Gabor spectrogram, D, controls the degree of
smoothing from Qian (2002). Figures 4.6 to 4.8 illustrate the run-up and steady-state
signals, respectively, and their Gabor spectrograms are where all the components can be
clearly seen. The abscissa is time, the ordinate is frequency and the intensities of the
color represent the magnitudes of the coefficients.
After computing the Gabor coefficients by using Eq. (4-2), one can reconstruct
the components using a time-varying filter, which is actually a two-dimensional binary
mask function, to modify the Gabor coefficient as
98
nmnmnm
CMC
,,,
?
= . (4-4)
The component of interest can then be extracted.
As long as some requirements are satisfied in Qian (2002), the component can be
reconstructed as
??
?
=
?
=
?=
1
0
1
0
2
,
][
?
][
M
M
N
n
N
nkj
nm
emTkhCkx
?
, (4-5)
where h[k] is called the synthesis function and the integer M represents the total number
of time points give by . In reference Qian (2002), it is shown that if the
functions h[k] and ?[k] are identical, the Gabor coefficients of the resulting waveform in
the time domain will be optimally close to , in the sense of least square error. This
process is called orthogonal-like Gabor transformation, by Qian (2002).
TLM /=
nm
C
,
?
4.2.3 Comparison of Wavelet and Gabor Analysis
Some previous works done by wavelet analysis have been introduced. Wavelet
and Gabor analyses are two sorts of joint time-frequency analyses. The continuous
Wavelet and Gabor transforms of a time signal x(t) are defined as
?
+?
??
?
>==< dt
a
bt
tx
a
xbaCWT
ba
)()(
1
,),(
*
,
?? , (4-4)
99
and
dtebtgtxgxbCGT
tj
b
?
?
?
?
+?
??
?
?>==< )(*)(,),(
,
, (4-5)
respectively, where
)(
1
,
a
bt
a
ba
?
= ?? , (4-6)
tj
b
ebtgg
?
?
)(
,
?= , (4-7)
and * denotes the complex conjugate.
Comparison of Equations (4-4) and (4-5) show that they are similar, but that the
elementary functions are different. The wavelet transform uses the dilated and translated
version of the mother wavelet ?(t) to decompose the signal, while the Gabor transform
uses the modulated and shifted copy of g(t). Because the basis functions are concentrated
at a region (b, ?) on the time-frequency plane, the Gabor transform can indicate how the
frequency changes over time.
It is shown that the bandwidth of the elementary functions in the wavelet
transform varies with its center frequency, while the bandwidth of the elementary
functions in the Gabor transform does not. This difference leads to the fact that the time-
frequency resolutions of these two transforms are different. The wavelet transform has
high time-resolution but low frequency-resolution at high frequency, and high frequency-
resolution but low time-resolution at low frequency.
100
Figure 4.4 illustrates a typical run-up signal of a cracked rotating shaft, and its
wavelet scalogram. In Figure 4.4 (b), although the magnitude is displayed in decibels,
only the component associated with 1X rotating speed can be seen. All the other
components become blurred in the scalogram. However, as Qian (2002) explained that
the Gabor transform has a constant time-frequency resolution once the analysis window
is determined.
Although both wavelet and Gabor transforms can describe how the frequency
components change with time, the wavelet transform is more suitable for transient
signals, such as engine knocks, because the basis functions in the wavelet transform have
finite-duration. On the other hand, with the Gabor coefficients, the basis functions have to
be long enough to cover a certain number of oscillations. Since the signals obtained from
the cracked rotors are still relatively stationary in a short time, the Gabor transform is
more applicable.
4.3 Experimental Results
An experimental facility is used to obtain data for analysis in this study. It
consists of a ? inch diameter shaft held vertically by two ball bearings and a motor
coupling, as shown in Figures 4.5(a) and 4.5(b). Proximity probes affixed to the bearing
supports are used to measure the vibration signals. A National Instruments PCMCIA
6036E card and LABVIEW software are used for data acquisition and analysis. Cracks
are deliberately introduced into the shaft for various depths and axial locations using a
metal cutting saw.
For a steel shaft, E=200*10
9
Pa, ?=7850 kg/m
3
. The diameter of the shaft used is
101
0.015875m, the length of the shaft between the two bearings is 0.365m. Then the
fundamental natural frequency is:
0
4
236
2
EI
f
mL
?
==Hz. (4-8)
4.3.1 Natural Frequency Measurements
First, the natural frequencies are measured for shafts with a variety of cracks
depths and locations. For each case, the shaft was inserted into the support bearings and
vibratory responses were excited using an impulse hammer. The resulting data was
analyzed using the covariance method to extract the natural frequencies. Table 4.1 shows
sample results for an intact rotor and for two rotor shafts with 5 mm depth cracks at
separate locations. For each sample, 20 hammer tests were conducted and the results
averaged to obtain the resulting value. The measured natural frequency of intact shaft is
226 Hz, about 4% less than the fundamental natural frequency, which means the
measured value is reliable.
As expected, the natural frequency of the intact shaft is higher than either of those
with transverse cracks. Interestingly, the cracks acted to lower the effective stiffness of
the shaft. The sample with the mid-span crack has a lower natural frequency than that of
the shaft with the crack near to its end, which is also an expected result. Similar results
were obtained for configurations with other crack depths. It should be noted that there
are two basic parameters that influence natural frequency changes that occur as a result of
102
transverse crack development, crack depth and crack axial location.
4.3.2 Rotor Vibration Analysis
Next, a series of tests was conducted for the purpose of identifying differences
between the vibration signatures of intact and cracked rotor shafts. Tests were conducted
for a variety of rotor speeds and conditions, including run-up and steady-state operation.
Figure 4.6 shows the time responses and Gabor spectrograms for the vibration
measurements from an intact shaft and from a shaft with a 5 mm depth crack at its center.
As mentioned earlier, the shafts were mounted vertically so as to remove the effects of
gravity from the test rotors. The rotor speed was increased in a linear fashion up to 4500
RPM over 10 seconds and then held steady for an additional 8 seconds.
Figure 4.6(b) shows the Gabor spectrogram for the intact shaft. The synchronous
vibration associated with the rotor imbalance eccentricity is readily observed. There are
also a variety of multi-synchronous components present in this spectrogram. The 2X, 3X,
and 4X components are quite distinct for all running speeds and the 5X and 6X can also
be observed for the higher running speeds. After the running speed reaches 4500 RPM
(75 Hz), there is a noticeable widening of the spectrogram associated with the 3X
component (225 Hz). This is due to the presence of the rotor natural frequency close to
this component (at about 226 Hz).
The Gabor spectrogram for the cracked shaft can be seen in Figure 4.6(d). Again,
the synchronous component (1X) is readily discernible for all running speeds, as is the
2X component. However, the higher harmonics are not nearly as clearly defined as for
the intact shaft. In addition, the natural frequency of the cracked shaft appears when the
103
running speed reaches about 35 Hz and is visible for the remainder of the test. Sideband
frequencies for several of the harmonics of the running speed can also be seen.
In order to obtain more detailed insights into the basic features of the
spectrograms for the two shafts, additional testing was conducted at a fixed running speed
of 4500 RPM. Figure 4.7 shows the FFT-based power spectrum and Gabor spectrograms
for the intact shaft. The order of the Gabor spectrogram is four. Figure 4.7 (a) shows the
frequency range up to 256 Hz. The first three harmonics of the running speed can be seen
clearly, with the magnitudes of the three harmonics decreasing with their orders. This
indicates that the energy is concentrated in the first harmonic of the running speed. Figure
4.7 (b) through (d) show the Gabor spectrograms of the first three harmonics separately.
There is no indication of any sidebands associated with any of these harmonics.
The FFT-based power spectrum and the Gabor spectrograms for the cracked shaft
are shown in Figure 4.8. In contrast to the spectrograms for the intact shaft, here the
magnitudes of the 2X and 3X vary periodically in the time domain. In addition, the
magnitude of the fundamental bending natural frequency is even higher than the second
and third harmonics of the running speed. This is probably due to the fact that it is close
to the 3X component. Twenty Hz sideband components also appear for the first and
second harmonics (1X and 2X) of the running speed.
4.4 Conclusions
A study of the application of Gabor analysis techniques to the detection and
monitoring of lateral cracks in rotor shafts has been conducted. Experimental vibration
data generated from a laboratory test rig was analyzed, with comparisons made between
104
results for an intact shaft and for a sample shaft with a machined crack. As expected, the
stiffness of the shaft decreases as the crack depth increases and when the location of the
crack is more toward the shaft center. The rotor-speed run-up results also show that the
crack reduces the critical speed of the rotor systems. Comparisons of the spectrograms for
the intact and cracked shafts show distinctive differences between the relative magnitudes
of the various major frequency components. These include the prominence of the rotor
natural frequency and the presence of sidebands in the cracked shaft results. From this
study, the Gabor expansion based analysis is shown to be a good and achievable tool for
the current work.
105
Linear chirp
Reversed linear chirp
FFT Spectrum
(c)
(b)
(a)
Figure 4.1 Drawbacks of the Fourier Transform ---- lose time information
106
Figure 4.2 Drawbacks of the Fourier Transform ---- can?t track frequency
change
107
Figure 4.3 Gabor analysis used for frequencies distribution.
108
(c)
(b)
(a)
Figure 4.4 Comparisons of the wavelet analysis and Gabor analysis
(a) A typical run-up signal of a cracked rotating shaft; (b) The wavelet
scalogram; (c) The Gabor spectrogram.
109
(a) (b)
Figure 4.5 Experimental setup for crack detection.
110
Figure 4.6 Comparison of the run-up signals of intact (a) and cracked (c)
shafts and their corresponding spectrograms (b,d)
111
(e)
Figure 4.7 The Gabor coefficients of an intact shaft
(a) Overall frequency range; (b) the fundamental frequency; (c) the second
harmonic; (d) the third harmonic; and (e) the FFT-based power spectrum
112
sidebands
(d)
sidebands
2X, 151 Hz
1X, 75.5 Hz
the fundamental natural frequency, 207 Hz
3X, 227 Hz
Figure 4.8 The Gabor coefficients of a shaft with a 5 mm deep crack in the middle.
(a) Overall frequency range; (b) the fundamental frequency and the sideband components;
(c) the second harmonic of the running speed; (d) the FFT-based power spectrum.
113
Table 4.1 Measured natural frequencies for the intact and cracked shafts
5 mm crack at the mid-span 5 mm crack at the end
Intact
shaft
In-crack
direction
Off-crack
direction
In-crack
direction
Off-crack
direction
Natural frequency
(Hz)
226.28
209.619 214.12 212.989 215.792
Standard deviation
(Hz)
1.824
0.9172 1.1226 0.8941 1.1301
114
CHAPTER 5 VIBRATION ANALYSIS AND HEALTH MONITORING OF CRACKS
IN COMPOSITE DISK ROTOR SYSTEMS
Cracks are general defects in rotating systems and are a precursor to fatigue-
induced failure. Identifying the presence and growth of cracks is a critical concept for the
health monitoring and diagnostics of such systems. A combined computational and
experimental study of the vibration characteristics of a composite hub flywheel rotor
system with a cracked hub disk is presented.
5.1 Introduction
High-speed flywheel energy storage systems offer the potential for substantially
improved energy storage densities as compared to conventional chemical batteries and
are seriously being considered for advanced satellite applications. To improve energy
storage densities, the rotor hub (providing the shaft/rim interface in a flywheel rotor) is an
attractive candidate for reducing rotor mass. The hub tends to lie close to the shaft and
contributes little to the overall energy storage capacity while adding to the system mass.
Some candidate designs use composite materials and allow significant flexibility.
However, to date there has been little work specifically directed at understanding the
influence of cracks in composite hub flywheel systems. It is critical that reliable health
monitoring strategies be developed for these systems if they are to be a successful
115
technology.
Accordingly, the present work investigates the influence of cracks on the dynamic
characteristics of composite-disk hub flywheel systems. An experimental model was
tested and the results used to validate a finite element analysis model, which was then
used to examine the basic characteristics of such a system. A simplified rotor model with
a flexible hub was developed based upon insights from these investigations. Using this
model, a detailed parametric study was conducted in order to investigate the relationship
between natural frequency and rotational speeds with and without a crack.
5.2 Experimental Investigation
In order to gain some initial insights into the influence of cracks on the dynamic
characteristics of a flexible hub rotor system and to provide validation data for a finite
element analysis (FEA) study, the natural frequencies of a basic test article were
experimentally determined. This includes measurement in the ?in-crack? and ?off crack?
directions for both in-plane and out-of-plane vibration modes.
5.2.1 Testing Sample
A schematic diagram of the basic test article used in this study is shown in Figure
5.1 and photographs are shown in Figure 5.2. The test article consists of a steel shaft
supporting a thin composite disk with an attached steel rim. This test article is intended to
simulate a flywheel system with a fairly flexible composite hub and a relatively massive
rim (for energy storage). The outer radius is R
o
= 76.2 m and the inner radius is R
i
= 7
mm. The hub disk thickness is t = 3.53 mm. The two sides of the rim are attached to the
116
sample disk as shown in the diagram. The total thickness of the rim is 60 mm. The
difference between the inner radius and the outer radius of the rim is 32 mm. The testing
points for the out-of-plane and in-plane vibration measurements are marked with the
symbols:
? Representing the testing point in the plane of the disk for out-of-plane direction testing
? Representing the testing point in the plane of the disk for in-plane direction testing
The material properties for the composite disk are summarized in Table 5.1. The
steel rim has a mass of 9.09Kg.
5.2.2 Testing Rig
The experimental set-up used in this study is shown in Figure 5.2. It consists of
the test article clamped in block supports attached to the slip table of a vibratory shaker.
A data acquisition system and two accelerometers affixed to the rotor system are used to
obtain and record the vibration signals. A hoop direction crack is deliberately introduced
into the composite disk during the lay-up of the disk. There are two testing
configurations, one for the out-of-plane direction and another for the in-plane direction,
as shown in Figures 5.2(a) and 5.2(b), respectively.
5.2.3 Experimental Results and Discussions
Using this setup, a series of experimental tests were conducted in order to identify
differences between the vibration signatures of the ?in-crack? and ?off-crack? directions
and to determine the natural frequencies in both directions. Tables 5.2 and 5.3 show some
117
sample results for the out-of-plane direction tests. Four points located in the plane of the
disk were picked as the testing points, as shown in Figure 5.1. Point ?a? lies in ?in-crack?
direction, points ?b? and ?d? are each located 90
o
from ?a? in the ?off-crack? directions,
and point ?c? is located 180
o
degrees from ?a?.
Similarly, sample results for the in-plane direction tests are summarized in Table
5.4. Locations for four points are all on the outer edge of the rim. Point ?1? lies in the ?in-
crack? direction, points ?2? and ?4? are each located 90
o
degrees from point ?1? in the
?off-crack? directions, and point ?3? is located 180
o
degrees from point ?1?.
As expected, the natural frequencies associated with the ?off-crack? directions are
higher than those associated with the ?in-crack? direction, both for the out-of-plane and
the in-plane tests. This is due to the decrease in effective stiffness of the rotor in the ?in-
crack? directions associated with the presence of the crack.
5.3 Numerical Simulation Using Finite Element Analysis
5.3.1 3-D Model
For the next stage of this effort, a finite element model was developed using
ANSYS. The same geometry and material properties of the experimental test article are
used for this model. Figure 5.3(a) shows a drawing of the model rotor system.
Cylindrical coordinates are used, with x representing the radial direction, y representing
the hoop direction, and z is the axial direction. The hoop direction is also the fiber
direction. Figure 5.4 is the top view of hoop direction crack for disk model.
118
5.3.2 Grid Independence Study
The finite element method is a numerical analysis technique for obtaining
approximate solutions to a wide variety of engineering problems. It envisions the solution
region as built up of many small, interconnected subregions or elements. A finite element
model gives a piecewise approximation to the governing equations, and a solution region
can be analytically modeled or approximated by replacing it with an assemblage of
discrete elements. The accuracy of a numerical solution heavily depends on the grid size,
which is used for the given problem. In order to determine the proper grid size for this
study, the grid independency test was conducted to see how the choice of grid affects the
value of the natural frequency for this model. The test was mainly focused on grid size
along three directions: the disk hoop direction, the disk thickness direction and the radial
direction. Graphs n
?,
n
z
and n
R
represent the total element numbers divided along the
hoop circumference direction, disk thickness direction, and disk radial direction,
respectively. The natural frequencies of first, second, fourth, fifth and sixth mode are
commonly used as a sensitivity measure of the accuracy of the solution.
For a given model with n
z
= 4, n
R
= 24 and a hoop crack lies in the middle radius
with an angle of 60?, six different grid sizes are tested: n
?
= 36, 48, 60, 72, 84 and 120.
Figure 5.5 shows the dependence of natural frequencies of five modes on the grid size n
?
.
Comparison of the natural frequencies of second mode among six different simulation
cases reveals that finer grid densities give rise to smaller predicted values in an
asymptotic fashion, while the natural frequencies of the rest modes monitored (first,
fourth, fifth and sixth mode) increases as the grid becomes finer and finer. For the natural
119
frequencies of all the modes monitored, the two grid distributions n
?
= 84 and 120 gave
nearly identical results.
Similar tests were performed in order to decide the proper grid sizes along disk
thickness and radial directions. Figure 5.6 shows the variation of natural frequencies with
the grid size in disk thickness direction, n
z
. According to the results among the five cases
we tested (n
z
= 1, 2, 4, 6 and 8), one can see that the results from n
z
= 6 and 8 are almost
identical. Results from Figure 5.7 imply that natural frequencies are not sensitive to n
R
when n
R
is greater than 24.
Considering both the accuracy and the computational time, the following
calculations of this study were all performed with n
?
= 84, n
z
= 6 and n
R
= 24.
5.3.3 Comparison of the Results between the Experimental and Simulation
Configurations with and without a crack in the composite hub were investigated
using this model. The location of the crack is the same as in the experimental test article
in order to do a close comparison between numerical simulation results and experimental
results. Some example results are shown in Table 5.5.
For the cracked case, the first two modes are axis-symmetric bending modes, with
two distinct frequencies associated with the ?in-crack? and ?off-crack? directions. The
third mode is a torsional mode, which we are not interested in for this study. The fourth
mode is an umbrella-type bending mode, with the associated natural frequency of 52.5
Hz, as shown in Figure 3(b). The fifth and sixth modes are in-plane stretching modes,
with two distinct frequencies associated with the ?in-crack? and ?off-crack? directions,
120
respectively, as shown in Figures 5.3(c) and 5.3(d).
A comparison of the results from the experimental testing and finite element
simulation is shown in Table 5.5. The natural frequencies associated with all the modes
of the cracked case are smaller than those for the cases without a crack. This implies that
the presence of a crack noticeably influences the natural frequencies for the primary
modes. It can also be observed that the frequency values for the ?in-crack? direction are
consistently smaller than those for the ?off-crack? direction.
5.3.4 Parameters of Numerical Simulations
Furthermore, a parameter variation study was conducted in order to deeply
examine the influence of crack length and crack position on the different natural
frequency. To simplify description, physical modes are replaced by number-mode and
are listed in Table 5.7. For example, the first mode is the mode for out-of-plane 1
st
bending axis-symmetric mode in-crack direction, the second is the mode for out-of-plane
1
st
bending axis-symmetric mode off-crack direction; the third is the 1
st
torsional mode,
and is not of concern in this study and therefore ignored. The fourth is the out-of-plane
2
nd
bending umbrella mode, the fifth is the 1
st
in-plane stretching mode in-crack direction,
and the sixth is the 1
st
in-plane stretching mode off-crack direction.
The parameter studies are as follows:
? Crack length growing along the disk hoop direction
? Crack length moving along the disk radial direction
? Crack length growing along the disk radial direction
121
5.3.5 Numerical Simulation Results and Discussions
5.3.5.1 Natural Frequencies Variation with a Variable Length Crack along the Disk Hoop
Direction
Figure 5.8 shows the crack?s location and its growing direction which is
investigated in this study. The crack is located along the hoop direction at the middle
radius of the disk sample, the starting point P
1
is unchanged and the crack length varies
from 0% to 97% (0?~330?) of the disk circumference. The case with 0% of the disk
circumference means no crack, while the case with 97% of the disk circumference
implies that there is a long crack along the circumference.
Figure 5.9 shows the variation of the natural frequencies associated with the five
different modes, respectively. The y-axis is the natural frequency F
cj
normalized by F
nj
,
where F
cj
is the natural frequency for the model with the crack, and F
nj
is the natural
frequency for the case without the crack. The x-axis is the length of crack L
c
normalized
by C
r
, where C
r
is the circumference of the middle radius r of the disk model.
Based on these figures, it is observed that all the natural frequencies decrease
when the crack length increases. Figures 5.9 (a) and (b) show the variation of natural
frequencies associated with first and second modes, respectively. Upon comparing (a)
and (b), one can see that the natural frequencies of out-of-plan 1
st
bending mode decrease
faster in ?in-crack? direction than in ?off-crack? direction. The variation of natural
frequencies of in-plane 1
st
stretching mode in both ?in-crack? and ?off-crack? directions
with the increasing of crack length are shown in Figures 5.9 (d) and (e), respectively.
Similar trend also can be observed in Figures 5.9 (d) and (e). For a given crack length,
natural frequencies in ?in-crack? direction are always smaller than that in ?off-crack?
122
direction, both for out-of-plane and in-plane modes. It implies that the existence of crack
has more remarkable effects on the stiffness in the ?in-crack? direction than in the ?off-
crack? direction.
The variation of the natural frequencies (shown in Figure 5.9(c)) associated with
the fourth mode exhibit a different trend compared with the others in Figure 5.9. The
natural frequency decreases slowly as the ratio of L
c
/ C
r
reaches 0.5. It corresponds to
the sector angle of crack length equals to ?. It indicates a phase change for the mode
shape. This is because the fourth mode is the 2
nd
out-of-plane bending mode, also referred
to as the umbrella type. A crack has more influence on axial direction or z-axis than in-
crack or off-crack direction. As the crack length increases, the imbalance eccentricity
center of mass U
m
moves down towards off-crack direction when the crack length is less
than half of the circumference. When the crack length is greater than half of the
circumference, U
m
moves up back to the disk center.
Figure 5.10 shows normalized natural frequency variations with the crack length
growing along the disk hoop direction for all five modes. Again, with the increase of
crack length, the relative change of natural frequencies associated with the first mode is
greater than any other modes. The min value of F
cj
/ F
nj
is about 0.2. It implies the crack
length (size) along the hoop direction in a disk has more noticeable effects on the low
modes.
5.3.5.2 Natural Frequencies Variation with a Constant Length Crack at Different
Locations along the Disk Radial Direction
In order to evaluate how the location of crack affects the natural frequencies, a
123
series of cases are generated in this study. Figure 5.11 shows the schematic diagram of
the parameter change. A crack located along the radial direction has a fixed length of
1/8
th
of ?R. The center of the crack locates seven points that are evenly spaced between
the inner and outer radius of the disk. Natural frequencies associated with the first six
modes are obtained from the 3-D finite element model for the above seven cases.
Figure 5.12 shows the variation of the natural frequencies associated with the five
different modes, respectively. Similarly, the y-axis is F
cj
normalized by F
nj
while the x-
axis is R
cc
normalized by R
o
.
Based on these figures, it is observed that the natural frequencies change when the
crack position changes. Examining curves for the first mode (shown in Figure5.12a),
second mode (shown in Figure 5.12b), fifth mode (shown in Figure5.12d) and sixth mode
(shown in Figure5.12e), one can conclude that the natural frequencies associated with the
?in-crack? direction are more sensitive to the location of crack than those for the ?off-
crack? direction. As shown in Figure 5.12a, the natural frequencies associated with the
first mode increase, as the crack is closer to the outer side of the disk. In other words, a
crack closer to the inner side of the disk has more remarkable effects on natural
frequencies associated with the first mode. Also, one can see that the whole curve is
below 1.0. This means that the natural frequencies associated with the first mode with a
crack are always smaller than that without a crack. Figure 5.12d shows the natural
frequencies associated with the fifth mode. The natural frequency reaches its minimum
when R
cc
/ R
o
= 0.5.
Although the imbalance eccentricity center of mass U
m
doesn?t change as long as
the crack length is constant, the moment of the centrifuge about the axial direction
124
changes when the crack location moves from the inner radius to outer radius. Figure
5.12(c) shows the variation of the natural frequencies associated with the fourth mode. It
is very different from other modes; the maximum decrease of the frequency is at the outer
disk edge, and then increases gradually as the crack moves toward the inner disk side.
This implies that a crack close to the outer disk side exhibits its most notable effects on
the natural frequency associated with the fourth mode.
The natural frequencies variation with location of cracks for all five modes is
plotted in Figure 5.13. For the case studied, all the curves are within a band with 0.96 <
Fc/Fn < 1.0. This implies that the crack position that changes along the radial direction in
a disk does not have important effects on the natural frequency.
From the above results, it can be noted that the change in the different modes
associated with the natural frequencies can be used to pinpoint the crack location.
5.3.5.3 Natural Frequencies Variation with a Variable Length Crack along the Disk
Radial Direction
The parametric study was repeated for the same disk rotor systems by considering
the crack to be lengthened in the disk radial direction. Figure 5.14 shows the crack?s
positions and lengths that we studied. The crack?s center is kept at point P
1,
which lies in
the middle of radius r. the two ends of crack extend to the disk inner side and outer side,
respectively. Six different lengths of crack are studied: 17.14%, 34.29%, 51.43%, 68.57%
and 85.7% of ?R, where ?R is the difference between the outer and the inner radius of
the disk).
Figure 5.15 shows the variation of the natural frequencies associated with the five
125
different modes, respectively. The y-axis is F
cj
normalized by F
nj
, where F
cj
is the natural
frequency for the model with the crack, and F
nj
is the natural frequency for the case
without the crack. The x-axis is L
c
normalized by ?R, which is the ratio of the crack
length and the difference between the outer and the inner radius of the disk model.
It is observed that all the natural frequencies decrease as the crack length
increases along the disk radial direction. After comparing the first mode (shown in
Figure5.15a) with the second mode (shown in Figure 5.15b), the fifth mode (shown in
Figure5.15d) with the sixth mode (shown in Figure5.15e), and one can see that the
natural frequencies associated with the ?in-crack? direction decrease faster than the
corresponding values for the ?off-crack? direction. It is also observed that the change of
the natural frequencies associated with the lower modes is greater than that of the higher
modes as the crack length increases. Figure 5.15 (c) shows the variation of the natural
frequencies associated with the fourth mode, the natural frequencies monotonously
decrease as the crack length increases. This is because the imbalance eccentricity center
of mass U
m
moves down towards the off-crack direction with the increasing of the crack
length. This is the noticeable difference as compared with Figure 5.9 (c).
Figure 5.16 shows the natural frequency variation as the crack length increases
along the disk radial direction for all five modes. Similarly, with the increase of the crack
length, the natural frequencies associated with the first mode decrease faster than any
other modes. The minimum value of F
cj
/ F
nj
is about 0.68. This implies that the natural
frequencies associated with the lower modes are more sensitive to the crack length along
the radial direction than at the higher modes. This is very similar to what we have
observed for the cases with crack length changing in the hoop direction.
126
5.4 Simplified Rotor Model and Simulation
Based upon these results, the analytical study was expanded to consider the
combined effects of rotor speed and crack size on natural frequency characteristics. A
schematic diagram of a simplified model developed for this purpose is shown in Figure
5.17. This model is meant to capture the characteristics of the primary in-plane modes (5
and 6). It consists of a rigid massless shaft and a flexible hub attached to discrete masses
representing the rigid rim. The shaft is symmetrically supported by springs of stiffness,
k
x
, k
y
and viscous damping elements, C
x
, and C
y
at both ends. The flexible disk and rim
are simulated by four identical mass elements m
1
, m
2
, m
3
, and m
4,
with the associated
deformations, ?
1
, ?
2
, ?
3
, and ?
4
. The spring elements that connect the mass elements to the
shaft are k
1
, k
2
k
3
, and k
4
. ? is the rotational speed and is assumed to be constant.
5.4.1 Positions of the Disk Masses
The equations of motion are developed in the following manner. First, the
location of the i
th
disk mass is denoted as:
( )[ ] [ ]JY
i
t
i
vRIX
i
t
i
vR
i
p
))
v
+++++++= )sin()()cos( ????
,
1, 2,3, 4,
(1).
2
i
i
i
?
?
=
=?
(5-1)
5.4.2 Velocities of the Disk Masses
The velocity of each disk mass is then obtained by differentiating the
127
corresponding position vector with respect to time, yielding:
[ ]
()[]
JY
i
t
i
vR
IX
i
t
i
vRK
JY
i
t
i
vR
i
t
i
v
IX
i
t
i
vR
i
t
i
v
i
V
)
))
)
&
&
v
&
&
v
++++
+++?+
++++++
+++?+=
)sin()(
)cos(
)cos()()sin(
)sin()()cos(
??
???
?????
?????
. (5-2)
5.4.3 Kinetic and Potential Energy of the Rotor
The kinetic energy of the rotor consists of two parts, the kinetic energy of the disk
and the kinetic energy of the shaft, and is represented as:
?
=
++?=
4
1
22
)(
2
1
)(
2
1
i
hiii
YXmVVmT
&&
vv
. (5-3)
Assume all the discrete masses are identical and equal to m, substitution of the
equations (2) into equations (3) gives:
)(
2
1
)]}(2)(2)[{sin(
)]}(2)(2)[{cos(
)](2
)([2)(2
2
1
22
4
1
2
4
1
2
22222
4
1
22
4
1
2
YXm
vRYvRXvXvYtm
vRXvRYvYvXtm
YXYX
YXYXmvRmvmT
h
i
iiiii
i
iiiii
i
i
i
i
&&
&
&&
&
&
&&
&
&&
&&
&
++
+++?+++
++++?++
?+
++++++=
?
?
??
=
=
==
?????
?????
?
??
.
(5-4)
The potential energy of the complete rotor system is the total energy stored in all
128
of the spring elements and is represented as:
22
4
1
2
2
1
2
1
2
1
YkXkvkP
yx
i
iipot
++=
?
=
.
(5-5)
5.4.4 Equations of Motion of the Rotor
Substitution of equations (4) and (5) into Lagrange?s equations, produces the
equations of motion.
44
11
4
22
1
48 cos( )()4 sin( )
4[cos()()] 4
ii i i
i
iih x
i
mX m Y m t v m v t
mtRvmXmXcXk
??????
??? ?
==
=
?+ +? +
?++?+
??
?
&& &
&& &
&& &
0X=
0Y
, (5-6.a)
44
11
4
22
1
4 8 sin( )( ) 4 cos( )
4[sin()()] 4
ii i i
iih y
i
mY m X m t v m v t
mtRvmYmYcYk
? ?? ? ??
??? ?
==
=
++ ++ +
?++?+
??
?
&& &
&& &
&& &
=
, (5-6.b)
22
2
cos( )( ) 2 cos( ) 2 sin( )
sin( ) sin( )( ) 4 ( )
cos( ) 0
ii i i i
ii
iii
mv m t X m Y t mX t
mY t m t Y m Rv
mt Xkv
? ??????
??? ?? ?
???
++?+++
?+ +?+
?++=
&& & &
&&
&&
?
. (5-6.c)
Only symmetric supports are considered, k
x
and k
y
are set with the same value,
EOM are converted to a constant coefficient form by transforming the support
129
coordinates (X, Y and Z) to the rotating frame ( X , Y and Z ), Eigen-analysis is then used
to analyze the resulting set of equations.
The transform relationship between space-fixed coordinate and shaft-fixed
coordinate system is shown Figure 5.18. The transform equations are:
cos( ) sin( )XX tY t? ?=?, (5-7.a)
sin( ) cos( )YX tY t? ?=+, (5-7.b)
where X, Y are coordinates for the space-fixed frame, and X ,Y are coordinates
for the shaft-fixed frame, Z and Z are in the coincident direction.
Substitution of equations (5-7.a , b) into equations (5-6), results in the EOM based
on shaft-fixed frame:
44
11
4
222
1
(4 ) (16 2 ) cos( )( ) 4 sin( )
4[cos()()](16 )
hh i
i
ii hx
i
mmX m m Y m v m v
mRvmmkXcc0Y
? ??
?? ?? ?
==
=
+?+ + ?
?+?+?+?
??
?
&& &
&& &
&
?
=
, (5-8.a)
44
11
4
222
1
(4 ) (16 2 ) sin( )( ) 4 cos( )
4[sin()()](16 )
hh i i
i
ii hy
i
mmY m m X m v m v
m Rv m m kYcYcX0
i
? ??
?? ?? ?
==
=
?+ ? + ? ?
++++??
??
?
&& &
&& &
&
?
=
, (5-8.b)
2
22
(4 4 )cos()
(4 4)sin()4()
ii i
ii
mv mX m Y m X
mY m X m Y m R v k v
?? ?
????
+? ?
++ ? ? ++=
&& &
&&
&& &
0
i
. (5-8.c)
130
5.4.5 Study of Numerical Simulation
The model represented by the equations of motion 5-8(a-c) was used to simulate
and analyze the system behavior for a variety of parametric configurations. Due to the
additional flexibility of the crack, the stiffness in the crack direction is softer than that of
off-crack direction. Thus, the presence of a crack is modeled as a smaller stiffness
coefficient, k
i
, associated with the quadrant in which the crack is located. By using same
parameters as the test sample, the Matlab codes (Appendix E) are run to simulate the
spectrum distribution of the cracked disk model. The parameters for the rotor model are
summarized in Table 5.8.
5.4.6 Results and Discussions
Figure 5.19 shows the rotational frequencies of the same rotor system both with
and without a crack. Forward and backward whirling frequencies for cases with and
without a crack are shown as dotted and solid lines, respectively. One can see that both
the solid and the dotted lines converge to a point when the rotational speed is zero. As
expected for both cases with and without a crack, the forward whirling frequency
increases as the rotational speed increases, while the backward whirling frequency
decreases as the rotational speed increases. However, there is a noticeable difference in
both the forward and backward frequencies between cases with a crack and without a
crack as the rotational speed increases. This implies that the presence of a crack has a
more remarkable influence for greater rotational speeds.
In addition, the relationship between the stiffness ratio and natural frequency is
studied. The coefficient k
c
stands for the spring stiffness simulating the rotor quadrant
131
with the crack. The case without a crack corresponds to k
1
/k
c
= 1. For a case with a crack,
k
1
/k
c
is always greater than one, as the crack results in a smaller effective k
c
. The
magnitude of k
c
decreases as the crack size increases. Inspection of Figure 5.20 shows
that:
(1) As the stiffness ratio increases (the crack size increases) both the forward and
backward natural frequencies decrease;
(2) When the rotational speed is zero, the fundamental frequency does not show a marked
change as the stiffness ratio increases;
(3) As the rotational speed increases, the frequency drops substantially as the stiffness
ratio increases, which is in agreement with the results shown in Figure 5.18. From
Figure 5.20, when the crack size grows, the possible tip rotating speed of the model
will decrease quickly. This is because the rotational speed can have a significant
effect on the crack behavior. The imbalance eccentricity of the center of mass U
m
grows as the speed of rotation increases even the crack size doesn?t change. This part
is discussed in details at Chapter 3.3.1 and Figure.3.7.
(4) The maximum stiffness ratio is smaller for higher rotational speeds. Obviously, this
indicates that the allowable crack size for a stable system decreases as the rotational
speed increases. The centrifugal effect of the operating speed acts as a negative
stiffness, reducing the effective stiffness of the overall system until catastrophic
failure occurs. Both the forward natural frequency (Figure 5.20(a)) and backward
natural frequency (Figure 5.20(b)) show similar tendencies.
A comparison of the results from the finite element analysis and theoretical model
simulation is shown in Table 5.9. One set data is from the ANSYS simulation parameter
132
study case (variations of the mode 5 in-plane 1
st
stretching mode as the crack grows along
the disk radial direction) and one from the simplified model case (variations of the
fundamental natural frequency as the stiffness ratio changes when the rotating speed is
zero).
From these two data sets, it is indicated that there is no big change for the
frequency value as the crack grows along the radial direction; both sets shows the same
tendency as the rotating speed is zero.
5.5 Conclusions
It is well known that the natural frequencies decrease due to the presence of a
crack. An investigation of the effects of a hoop-direction crack on the dynamic
characteristics of a flexible-composite-hub/rigid-rim rotor system has been presented.
This work has consisted of an experimental study, which was used to benchmark and
validate a finite element model. Using that model, a parameter variation study was
extended to examine the effects of crack length on the natural frequencies of the in-plane
and the out-of-plane vibration modes, as well as differences between similar modes for
the ?in-crack? and ?off-crack? directions. Finally, a simplified model for the primary in-
plane mode (a stretching mode) was developed and used to evaluate the effect of
operating speed on the frequency characteristics. Based upon the results of these studies,
insights into vibration characteristics that are potentially useful for the purposes of rotor
health monitoring were obtained. Specifically, it was observed that the presence of a
crack serves to reduce the effective natural frequency of the rotor system and this effect is
magnified as operating speeds increase. In addition, there is a directional dependence of
133
these frequency characteristics due to this softening effect arising from a crack. The
frequencies associated with the ?in-crack? direction are generally smaller than those
associated with the ?off-crack? direction.
Based on these results, it suggested that the changes of the natural frequencies
could be used for crack detection and the monitoring of rotor systems.
134
t
t
1
Crack
R
o
R
i
1
2
4
a
c
3
t
3
t
2
b
d
Figure 5.1 Schematic diagram of the experimental test article
135
(b)
(a)
Figure 5.2 Experimental setup for the frequency measurement of the composite crack
disk with heavy rim. (a) Out-of-plane direction test setup; (b) In-plane direction test setup.
136
1
MN
MX
X
Y
Z
0
.035909
.071819
.107728
.143637
.179547
.215456
.251365
.287275
.323184
DEC 1 2004
11:05:41
NODAL SOLUTION
STEP=1
SUB =6
FREQ=321.184
USUM (AVG)
RSYS=0
DMX =.323184
SMX =.323184
1
MN
MX
X
Y
Z
0
.036748
.073495
.110243
.14699
.183738
.220485
.257233
.293981
.330728
DEC 1 2004
11:03:54
NODAL SOLUTION
STEP=1
SUB =4
FREQ=52.45
USUM (AVG)
RSYS=0
DMX =.330728
SMX =.330728
1
X
Y
Z
NOV 23 2004
13:50:50
ELEMENTS
MAT NUM
1
MN
MX
X
Y
Z
0
.035964
.071927
.107891
.143854
.179818
.215781
.251745
.287708
.323672
DEC 1 2004
11:04:52
NODAL SOLUTION
STEP=1
SUB =5
FREQ=314.745
USUM (AVG)
RSYS=0
DMX =.323672
SMX =.323672
(a) (b)
(c) (d)
Figure 5.3 FEA model (a) Drawing and grid system; (b) Fourth mode shape with the
associated frequency 52.5 Hz for the crack case; (c) Fifth mode shape with the associated
frequency 315 Hz (in-crack direction); (d) Sixth mode shape with the associated frequency 321
Hz (off-crack direction).
137
Crack in hoop direction
Figure 5.4 Top view of disk crack model
138
17
17 . 2
17 . 4
17 . 6
17 . 8
18
18 . 2
18 . 4
18 . 6
18 . 8
19
36 48 60 72 84 96 108 120
Frequency (Hz)
Second mode
57.0
57.2
57.4
57.6
57.8
58.0
58.2
58.4
58.6
58.8
59.0
36 48 60 72 84 96 108 120
Frequency (Hz)
Fourth mode
310.0
311.0
312.0
313.0
314.0
315.0
316.0
317.0
318.0
319.0
320.0
36 48 60 72 84 96 108 120
Frequency (Hz)
Fifth mode
310.0
312.0
314.0
316.0
318.0
320.0
322.0
36 48 60 72 84 96 108 120
Frequency (Hz)
Sixth mode
17
17.2
17.4
17.6
17.8
18
18.2
18.4
18.6
18.8
19
36 48 60 72 84 96 108 120
Frequency (Hz)
First mode
n
?
(e)
n
?
(d)
n
?
(c)
n
?
(b)
n
?
(a)
Figure 5.5 The mesh grid changes along the disk hoop direction (at n
z
= 4, n
R
=
24, a hoop crack lies in the middle radius with the angle 60?)
139
17 . 0
17 . 2
17 . 4
17 . 6
17 . 8
18 . 0
18 . 2
18 . 4
18 . 6
18 . 8
19 . 0
0123456789
Frequency (Hz)
First mode
17 . 0
17 . 2
17 . 4
17 . 6
17 . 8
18 . 0
18 . 2
18 . 4
18 . 6
18 . 8
19 . 0
0123456789
Frequency (Hz)
Second mode
57.0
57.2
57.4
57.6
57.8
58.0
58.2
58.4
58.6
58.8
59.0
0123456789
Frequency (H
z)
Fourth mode
310.0
312.0
314.0
316.0
318.0
320.0
322.0
324.0
326.0
328.0
330.0
0123456789
Fre
que
nc
y
(Hz)
Fifth mode
310.0
312.0
314.0
316.0
318.0
320.0
322.0
324.0
326.0
328.0
330.0
0123456789
Frequency (Hz)
Sixth mode
(b)
(c)
(a)
n
z
n
z
(d) (e)
n
z
n
z
n
z
Figure 5.6 The mesh grid changes in the disk thickness direction (at n
?
= 84, n
R
= 24, a hoop crack lies in the middle radius with the angle 60?)
140
(a)
16.0
16.5
17.0
17.5
18.0
18.5
19.0
8 1216202428
NR
Frequency (Hz)
First mode
16.0
16.5
17.0
17.5
18.0
18.5
19.0
8 1216202428
Frequency (Hz)
Second mode
n
R
NR
(b)
n
R
40.0
44.0
48.0
52.0
56.0
60.0
64.0
68.0
8 1216202428
Frequency (Hz)
Fourth mode
NR
n
R
(c)
(d)
300.0
305.0
310.0
315.0
320.0
325.0
330.0
8 1216202428
Frequency (Hz)
Fifth mode
NR
n
R
(e)
300.0
305.0
310.0
315.0
320.0
325.0
330.0
8 1216202428
Frequency (Hz)
Sixth mode
NR
n
R
Figure 5.7 The mesh grid changes along the disk radial direction (at n
z
= 6, n
?
=
84, a hoop crack lies in the middle radius with the angle 60?)
141
R
i
R
o
r
P
2
P
1
L
c
Figure 5.8 The path of a crack growing in the hoop direction.
(The crack start point P
1
is fixed, the other end extends along the middle radius
( r = ? (Ri+Ro) ) circumference C
r
, until reaches point P
2.
)
142
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Fc1 / Fn1
First mode
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Fc2 / Fn2
Second mode
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Fc4 / Fn4
Fourth mode
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Fc5 / Fn5
Fifth mode
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Fc6 / Fn6
Sixth mode
Lc / Cr
Lc /
(b)
(a)
Lc / Cr
(c)
Lc / Cr Lc / Cr
(d) (e)
Figure 5.9 Natural Frequencies variation with a variable crack length along the
disk hoop direction
143
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Fc / Fn
First mode Second mode Fourth mode
Fifth mode Sixth mode
Lc / Cr
Figure 5.10 Natural Frequency variation with the variable crack length
along the disk hoop direction (all modes)
144
P
3
P
2
R
o
R
cc
Figure 5.11 The path for crack moving in radial direction.
(The crack length is fixed, the crack center position starts from point P
2
, moves to P
3.
)
145
0.95
0.96
0.97
0.98
0.99
1. 0 0
1. 0 1
0 0.2 0.4 0.6 0.8 1
Rcc / Ro
Fc2 / Fn2
Second mode
(b)
0.95
0.96
0.97
0.98
0.99
1. 0 0
1. 0 1
0 0.2 0.4 0.6 0.8 1
Rcc / Ro
Fc1 / Fn1
First mode
(a)
0.95
0.96
0.97
0.98
0.99
1. 0 0
1. 0 1
0 0.2 0.4 0.6 0.8 1
Rcc / Ro
Fc4 / Fn4 Fourth mode
0.95
0.96
0.97
0.98
0.99
1. 0 0
1. 0 1
0 0.2 0.4 0.6 0.8 1
Rcc / Ro
Fc5 / Fn5
Fifth mode
(d)
(c)
0.95
0.96
0.97
0.98
0.99
1. 0 0
1. 0 1
0 0.2 0.4 0.6 0.8 1
Rcc / Ro
Fc
6
/ Fn6
Sixth mode
(e)
Figure 5.12 Natural Frequency variation with a constant length crack at different
locations along the disk radial direction
146
0.95
0.96
0.97
0.98
0.99
1.00
1.01
0 0.2 0.4 0.6 0.8 1
Rcc / Ro
Fc / Fn
First mode Second mode Fourth mode
Fifth mode Sixth mode
Figure 5.13 Natural Frequency variation with a constant length crack at different
locations along the disk radial direction (all modes)
147
?R
P
3
P
1
P
2
L
c
Figure 5.14 The path for crack growing in radial direction.
(The crack center is fixed at point P
1
, the two ends of crack extend to P
2
and P
3
.)
148
149
0.6
0.7
0.8
0.9
1
1. 1
0 0.2 0.4 0.6 0.8 1
Fc6 / Fn6
Sixth mode
Lc / RLc / ?R
(e)
0.6
0.7
0.8
0.9
1
1. 1
0 0.2 0.4 0.6 0.8 1
L
Fc1 / Fn1
First mode
c / DRLc / ?R
(a)
0.6
0.7
0.8
0.9
1
1. 1
0 0.2 0.4 0.6 0.8 1
Fc2 / Fn2
Second mode
Lc / R
Lc / ?R
(b)
0.6
0.7
0.8
0.9
1
1. 1
0 0.2 0.4 0.6 0.8 1
L
Fc4 / Fn4
Fourth mode
c / RLc / ?R
(c)
0.6
0.7
0.8
0.9
1
1. 1
0 0.2 0.4 0.6 0.8 1
Fc5 / Fn5
Fifth mode
Lc / RLc / ?R
(d)
Figure 5.15 Natural Frequency variation with a variable length crack along the
disk radial direction
0.6
0.7
0.8
0.9
1
1.1
0 0.2 0.4 0.6 0.8 1
Lc / R
Fc / Fn
First mode Second mode Fourth mode
Fifth mode Sixth mode
Lc / ?R
Figure 5.16 Natural Frequency variation with a variable length crack along
the disk radial direction (all modes)
150
?
k
1
k
2
k
3
k
4
?
4
?
3
?
2
?
1
4
3
2
1
Z
Y
X
Figure 5.17 Simplified rotor model with flexible hub
151
Z, Z
?
X
Y
X
Y
?
?
Figure 5.18 Coordinate systems for support frame and rotating shaft frame
152
100
150
200
250
300
350
400
0 100 200 300 400 500 600
Rotational speed (Hz)
Natural frequency (Hz)
balanced system (no crack)
unbalanced system (with a crack, K1/Kc=10)
Figure 5.19 Influence of a crack on forward and backward whirling
natural frequencies
153
100
150
200
250
300
350
0123456
Stiffness ratio Log(k1/kc)
Natural frequence (Hz)
w=0 Hz w=100 Hz w=300 Hz w=500 Hz
300
310
320
330
340
350
360
370
380
390
0123456
Stiffness ratio Log(K1/Kc)
Natural Frequence(Hz)
w=0 Hz w=100 Hz w=300 Hz w=500 Hz
(a)
(b)
Figure 5.20 Natural frequency variation of the primary mode with stiffness ratio for selected
operating speeds (a) Forward; (b) Backward
154
Table 5.1 Material properties of the composite disk
Ex 6.67 GPa
Ey 140 GPa
Ez 6.67 GPa
Gxy 3.20 GPa
Gxz 2.39 GPa
Gyz 3.19 GPa
Nuxy 0.260
Nuyz 0.260
Nuxz 0.395
?
1520 kg/m
3
155
Table 5.2 Measured Natural Frequencies for Out-of-Plane 1
st
Bending Mode
Test Point Natural Frequency (Hz)
a 16.3
b 17.0
c 16.3
d 17.0
156
Table 5.3 Measured Natural Frequencies for Out-of-Plane 2
nd
Bending Mode
Test Point Natural Frequency (Hz)
a 53.0
b 53.3
c 53.0
d 53.3
157
Table 5.4 Measured Natural Frequencies for the 1
st
In-Plane Mode
Test Point Natural Frequency (Hz)
1 300
2 316
3 312
4 316
158
Table 5.5 Natural frequencies and mode shapes from finite element analysis
Natural Frequency
(Hz)
Mode
Without Crack With Crack
16.5
(in-crack direction)
1
st
Bending
axis-symmetric
mode
16.9
16.8
(off-crack direction)
Torsional mode
43.7 43.6
2
nd
Bending
umbrella mode
53.6 52.5
315
(in-crack direction)
1
st
In-plane
stretching mode
325
321
(off-crack direction)
159
Table 5.6 Comparison of natural frequencies from finite element
analysis and experimental testing
Mode Simulation
(Hz)
Experiment
(Hz)
Error (%)
(In-crack
direction)
16.5 16.3 1.2
Out-of-
plane
1
st
Bending
(Off-crack
direction)
16.8 17.0 1.2
Out-of-plane
2
nd
Bending
52.5 53.2 1.4
(In-crack
direction)
315 300 4.6 1
st
In-plane
(Off-crack
direction)
321
314 2.2
160
Table 5.7 Explanation of terms used to describe modes
Number-Mode in FEA model
Physical Modes
In-crack direction Off-crack direction
Out-of-plane
1
st
bending
axis-symmetric mode
First mode Second mode
1
st
Torsional mode
Third mode (not interested)
Out-of-plane
2
nd
bending
Umbrella mode
Fourth mode
In-plane
1
st
Stretching mode
Fifth mode Sixth mode
161
Table 5.8 Parameters of the simple rotor disk model
Balanced case Unbalanced Case
Damping factor for the support
C
x
, Cy
0.8 0.8
Damping factor for the mass
c
1
c
2
c
3
c
4
0.5 0.5
System support spring
k
x
k
y
(Kg/m
3
)
5.0e7 5.0e7
The element spring of Mass
k
1
k
2
k
3
k
4
(Kg/m
3
)
5.0e7 k
2
k
3
k
4
=5.0e7
k
1
changes from 5.0e7 to 5.0e2
Mass of the disk
m
1
m
2
m
3
m
4
(Kg)
2.45 2.45
Mass of the hub
m
h
(Kg)
0.2 0.2
Radius of the hub
R (m)
0.076 0.076
162
Table 5.9 Comparison of natural frequencies from finite element analysis and
theoretical model simulation
ANSYS Model
Rates(Lc / ?R) 0 0.1714 0.3429 0.5143 0.6857 0.857
Fc/Fn 1 0.998543 0.994824 0.988718 0.979543 0.96482
Crack length grows along the disk radial direction (Data for the mode 5 in-plane 1
st
stretching mode)
Matlab Model
Kc/K1 0 0.1 0.01 0.001 0.0001 0.00001 0.000001
Fc/Fn 1 0.995944 0.97965 0.978678 0.978703 0.977838 0.977838
When the rotating speed is zero, the variations of natural frequency
163
CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS
A study aimed at developing (i) an effective optimal design methodology for
advanced composite rotors, and (ii) evaluating dynamic characteristics and developing
insights for health monitoring has been presented. The major contributions of this work
are:
1. In order to obtain the optimization design of a particular MDC flywheel with
the maximum energy storing density (ESD), a flywheel optimization algorithm has been
proposed and implemented. First, a methodology was developed to present an optimal
design procedure for the structure of MDC flywheel systems. To solve the nonlinear
optimization problem, a parametric model was built and optimization analysis was
performed by using finite element software ANSYS. Second, a specific optimal structure
design was successfully demonstrated of a domed hub-shaft-rim flywheel rotor. A
parametric study was performed by changing the thickness, t, of the hub port, the value of
the ratio a/b, and the maximum energy storing density (ESD) of the flywheel rotors was
calculated while the shear stress response of the model was recorded for each speed of
vibration. A specific optimized design for a MDC flywheel system was demonstrated and
a prototype fabricated. The results of this study confirmed the proposed model developed
might be useful for the optimization of future designs.
2. For safety issues of the MDC flywheel systems, the initiation and propagation
164
of a crack is one of the most important failures found in high rotation speed flywheel
systems. To understand the crack propagation behavior, general observations of the crack
propagation mechanism are presented in detail from the numerical investigations. The
relationships between the imbalance eccentricity of the mass center (U
m
) and the
variations in rotation speeds, crack growth size, material properties of the modulus of
elasticity, and the radii ratio ? (ratio of inner radius and outer radius) were studied
extensively using finite element analysis. The imbalance eccentricity of the center of
mass U
m
increases when the rotating speed is increasing. A crack lying in the radial
direction is more harmful than that that in the hoop direction for a flywheel. As a crack
length varies along the hoop direction, the imbalance eccentricity U
m
increases as the
crack length increases, where the crack is shorter than half the circumference. Then the
imbalance eccentricity U
m
decreases once the crack is longer than half the circumference.
It reaches its maximum value when the crack length equals half the circumference. A
flywheel with a small value of ? can reduce the imbalance eccentricity of the center of
mass, but this will lower the energy storage efficiency.
3. To obtain further insight into the influence of cracks in a flywheel system, more
experimental and simulation work were conducted.
First, a methodology to detect the dynamic behavior of a shaft rotor with cracks
using Gabor?s analysis techniques was developed. The application of a joint time-
frequency method, or Gabor analysis, as a tool for crack identification and health
monitoring in rotating machinery was investigated. Experimental vibration data was
generated for a set of sample shafts with different crack depths and locations. Gabor
165
analysis was then used to examine the characteristics of the respective frequency
components over time in order to develop a methodology for the identification of crack
properties based upon this information. From the study, one can see that the Gabor
analysis is not only robust and easy to implement, but also offers more insight into the
physical process. There is a significant agreement between the experimental results and
the theoretical calculations.
Next, the presence and growth of crack characteristics of a composite hub
flywheel rotor system with a cracked hub disk were studied. Experimental testing of both
in-plane and out-of-plane vibration characteristics, using a rotor with a composite disk
hub supporting a relatively massive rim, was conducted. A crack was deliberately
introduced into the hub disk during fabrication. Then, a finite element (FEA) model was
developed to explore further the relationship between natural frequencies and crack
properties. Finally, a simplified theoretical model for the primary in-plane vibration mode
was developed and used in a series of parametric studies. The analytical model was
validated by experimental measurements. It was observed that the presence of a crack
tends to affect both the magnitudes and distribution of the rotor natural frequencies.
Certain primary frequencies for rotors with a crack are smaller than for those without a
crack. In addition, the frequency values associated with the ?in-crack? direction are
generally smaller than those associated with the ?off-crack? direction, introducing non-
symmetry into the rotordynamics that can serve as an indicator for rotor health
monitoring.
The results obtained from this study should be extremely useful for optimizing the
design and safe operation of a MDC flywheel system. Suggestions for further work are
166
summarized as follows:
1. For the experimental study of crack detection of the shaft rotor, only the effects
of two cases (the crack at the mid-span and at the end) have been tested. Conducting
additional experimental tests and examining various crack depths and axial positions will
be very helpful to understand the crack behaviors further and will extend the conclusions
of this present study.
2. Due to the limitation of our experimental capability (no spin pit), rotating tests
of the cracked disk rotor were not performed study. Extending the experimental rotating
test of rotor disk and applying Gabor techniques for crack detection would be a useful
extension of my Gabor analysis work.
3. Consider the influence of non-symmetric supports on the rotor dynamic
behavior and overall stability of the simplified rotor model. This will produce time-
periodic coefficients and require Floquet analysis.
167
REFERENCES
Abdi, R. and Touratier, M., ?Shape Optimization of a Pressure Vessel under Plastic
Flow, Plastic Instability, Weight and Fatigue Life Criteria?, Computers and Structures,
Vol. 51, pp. 559-570, 1994.
Adewusi, S. A. and AL-Bedoor, B. O., ?Wavelet Analysis of Vibration Signals of an
Overhang Rotor with a Propagating Transverse Crack,? Journal of Sound and Vibration,
Vol. 246, No.5, pp. 777-793, 2001.
Altstadt, E. E., Werner, M. and Dzugan, J., ?FEM Simulation of Crack Propagation in
3PB and CT Specimens,? www.fzrossendorf.de/FWS/publikat/JB02/JB_02_R13.pdf.
Arnold, S. M., Saleeb, A. F. and Al-Zoubi, N. R., ?Deformation and Life Analysis of
Composite Flywheel Disk Systems?, Composites, Part B, Vol. 33, pp. 433-459, 2002.
Bachschmid, N., Pennacchchi, P., Tanzi, E., and Vania, A., ?Identification of Transverse
Crack Position and Depth in Rotor Systems,? Meccanica, Vol. 35, pp. 563-582, 2000.
Bicego, V., Lucon, E., Rinaldi, C., and Crudeli, R., ?Failure Analysis of a Generator
168
Rotor with a Deep Crack Detected During Operation: Fractorgraphic and Fracture
Mechanics Approach,? Nuclear Engineering and Design, Vol. 188, pp. 173-183, 1999.
Bitterly, J. G., ?Flywheel technology past, present, and 21
st
century projections?, IEEE
Aerospace and Electronic Systems Magazine, Vol. 13, No. 8, Aug, IEEE, pp 13-16, 1998.
Bowler, E. M., ?Flywheel Energy Systems: Current Status and Future Prospects?,
Magnetics Material Producers Association Joint Users Conference, Sept. 22-23, 1997.
Chen, L. W. and Chen, H. K., ?Stability Analysis of a Cracked Shaft Subjected to the
End Load,? Journal of Sound and Vibration, Vol. 188, No. 4, pp. 497-513, 1995.
Chudnovsky, A. and Li, R. S., ?Asymptotic analysis of the crack tip field in the vicinity
of a biomaterial interface?, American Society of Mechanical Engineers, Manufacturing
Engineering Division, MED, Proceedings of the 1995 ASME international Mechanical
Engineering Congress and Exposition, San Francisco, CA, Nov. 12-17, 1995.
Curtiss, D. H., Mongeau, P. P. and Puterbaugh, R. L., ?Advanced Composite Flywheel
Structural Design for a Pulsed Disk Alternator?, IEEE Transactions on Magnetics, Vol.
31, No. 1, 1995.
Darpe, A. K., and Gupta, K., ?Analysis of the Response of a Cracked Jeffcott Rotor to
Axial Excitation,? Journal of Sound and Vibration, Vol. 249, No. 3, pp. 429-445, 2002.
169
Darpe, A. K., Gupta, K. and Chawla, A., ?Experimental Investigations of the Response of
a Cracked Rotor to Periodic Axial Excitation,? Journal of Sound and Vibration, Vol. 260,
pp. 265-286, 2003.
Ebrahimi, N., ?Optimum design of flywheels?, Composites and Structures, Vol. 29, No.
4, pp. 681-686, 1988.
Emerson, R. P. and Bakis, C. E. ?Viscoelastic behavior of composite flywheels?,
International SAMPE Symposium and Exhibition (Proceedings), Vol. 45 No. II, SAMPE,
pp 2001-2012, 2000.
Fausz, J. L. and Richie, D. J., ?Flywheel simultaneous attitude control and energy storage
using a VSCMG configuration?, IEEE Conference on Control Applications -
Proceedings, Proceedings of the 2000 IEEE International Conference on Control
Applications, Anchorage, USA, Sep 25-Sep 27, 2000.
Fischer, G., Grothaus, R., and Hufenbach, W., ?Experimental Failure Effect Analysis of
High-speed Flywheel Rotors?, ASTM Special Technical Publication, Mar. 11-12 2002,
Pittsburgh, PA.
Fisher, C. and Lesieutre, G., ?Health Monitoring of High Energy Density Rotor for Use
in Spacecraft Applications,? International SAMPE symposium and exhibition
170
(proceeding), Vol. 44, No. 2, pp. 2145-2154, 1999.
Genta G., ?Kinetic energy storage ? Theory and practice of advanced flywheel systems?,
Butterworths, British, 1985.
Gentile, A., and Messina, A., ?On the Continuous Wavelet Transforms Applied to
Discrete Vibrational Data for Detecting Open Cracks in Damaged Beams,? Solids and
Structures, Vol. 40, pp. 295-315, 2003.
Gowayed, Y. and Flowers, G. T., ?Optimal Design of Multidirection Composite
Flywheel Rotors,? Polymer Composites, Vol. 23, No. 3, pp. 433-441, 2002.
Green, I., and Casey, C., ?Crack Detection in a Rotor Dynamic System by Vibration
Monitoring ? Part I: Analysis,? Proceedings of ASME Turbo Expro 2003, Atlanta, June
16-19, 2003.
Grudkowski, T. W., Dennis, A. J., Meyer, T. G. and Wawrzonek, P.H., ?Flywheels for
energy storage?, SAMPE Journal, Vol. 32, No. 1, Jan-Feb, pp 65-69, 1996.
Guo, D., Chu, F. and He, Y., ?Vibration Analysis of the Rotating Shafts with Transverse
Surface Cracks,? Proceedings of ASME 2003 Design Engineering Technical Conference,
Chicago, September 2-6, 2003.
171
Ha, S. K. and Jeong, H. M., ?Optimum Design of Thick-walled Composite Rings for an
Energy Storage System,? Journal of composite materials, Vol. 32, No.9, 1998.
Ha, S. K. and Kim, D. J., ?Optimal Design of a Hybrid Composite Flywheel Rotor
using Finite Element Methods,? International SAMPE symposium and exhibition
(proceeding), Vol. 44, No. 2, pp. 2119-2132, 1999.
Ha, S. K., Yang, H. and Kim, D. J., ?Optimal design of a hybrid composite flywheel with
a permanent magnet rotor?, Journal of Composite Materials, Vol. 33, No. 16, pp 1544-
1575, 1999.
Ha, S. K., Kim, D. J. and Sung, T. H., ?Optimum design of multi-ring composite
flywheel rotor using a modified generalized plane strain assumption?, International
Journal of Mechanical Sciences, Vol. 43, pp 993-1007, 2001.
Harmon, L. M. and Baaklini, G. Y., ?Ultrasonic Resonance Spectroscopy of Composite
Rings for Flywheel?, Nondestructive and Evaluation of Materials and Composites V,
Mar. 7-8, 2001, Newport Beach, CA.
Hauser, A.O., ?Calculation of Superconducting Magnetic Bearings using a Commercial
FE-program (ANSYS),? IEEE Transactions on Magnetics, Vol. 33, No 2, pp. 1572-1575,
1997.
172
He, Y., Guo, D. and F. Chu, ?Using Genetic Algorithms and Finite Element Methods to
Detect Shaft Crack for Rotor-bearing System,? Mathematics and Computers in
Simulation, Vol. 57, pp.95-108, 2001.
Horner, R. E. and Pround, N. J., ?Key Factors in the Design and Construction of
Advanced Flywheel Energy Storage Systems and Their Application to Improve
Telecommunication Power Back-up,? Proceedings of the 1996 18th International
Telecommunications Energy Conference, INTELEC, Oct 6-10 1996, Boston, MA, USA.
Huang B. C., ?Polar woven flywheel rim design?, International SAMPE Symposium and
Exhibition (Proceedings), Vol. 44, No. II, SAMPE, pp 2133-2144, 1999.
Huang, J. and Fadel G. M., ?Heterogeneous Flywheel Modeling and Optimization?,
Materials and Design, Vol. 21, pp. 111-125, 2000.
Jee, D. H. and Kang, K. J., ?A Method for Optimal Material Selection Aided with
Decision Making Theory?, Materials and Design, Vol. 21, pp. 199-206, 2000.
Kim, S. S. and Kim, J. H., ?Rotating Composite Beam with a Breathing Crack,?
Composite Structures, Vol. 60, pp. 83-90, 2003.
Kojima, A., Ikeda, T., Seki, S. and Itoh, T, ?Flywheel Energy Storage System?, 5th
JAPAN International SAMPE Symposium, pp.249-252, 1997.
173
Lashley, C. M., Ries, D. M., Zmood, R. B., Kirk, J. A. and Anand, D.K., ?Dynamics
Considerations for a Magnetically Suspended Flywheel,? Proceedings of the 24th
Intersociety Energy Conversion Engineering Conference - IECEC-89. Part 3,
Electrochemical Storage and Conversion, Aug 6-11, 1989, Washington, DC, USA.
Lee, H. P., and Ng, T. Y., ?Natural Frequencies and Modes for the Flexural Vibration of
a Cracked Beam,? Applied Acoustics, Vol. 42, pp. 151-163, 1994.
Lekhnitskii, S. G., ?Anisotropic Plates?, English translation, Gordon and Breach Science
Publishers, New York, NY, 1968.
Li, W. and Shen, Z., ?Composite Material Flywheel Structure and Energy Storing
Energy,? Taiyangneng Xuebao, Vol. 22, No. 1, pp. 96-101, 2001.
Li, Z., Xia, S., Fan, J., and Su, X., ?Damage Detection of Beams Based on Experimental
Wavelet Analysis of Flexible Waves,? Key Engineering Materials, v 261-263, No. II,
Advances in Fracture and Failure Prevention: Proceedings of the Fifth International
Conference on Fracture and Strength of Solids (FEOFS2003): Second International
Conference on Physics and Chemistry, 2004.
Liu, G., Wu, G.; Jia, D., Shangguan, J., Wang, Z. and Zhou, H., ?Application of Gabor
Transform for Detecting Crack Edge-reproduced Wave from Reactor Pressure Vessel,?
Hedongli Gongcheng/Nuclear Power Engineering, Vol. 19, No. 6, pp. 514-518, 1998.
174
Muller, G. and Tiefenthaler, P., ?Design Optimization with the Finite Element Program
ANSYS,? International Journal of Computer Applications in Technology, Vol. 7, No. 3,
pp. 271-277, 1996.
Niemeyer, W. L., Studer, P., Kirk, J. A., Anand, D. K. and Zmood, R. B., ? A High
efficiency motor/generator for magnetically suspended flywheel energy storage system,?
Proceedings of the 24th Intersociety Energy Conversion Engineering Conference -
IECEC-89. Part 3: Electrochemical Storage and Conversion, Aug. 6-11, 1989,
Washington, DC.
Peng, G., Yuan, S. F., and Xu, Y. D., ?Damage Location of Two-dimensional Structure
Based on Wavelet Transform and Active Monitoring Technology of Lamb Wave,?
Proceedings of SPIE - The International Society for Optical Engineering, Vol. 5393,
Nondestructive Evaluation and Health Monitoring of Aerospace Materials and
Composites III, pp. 151-160, 2004.
Portnov, G., Uthe, A. N., Cruz, I. Fiffe, R. P. and Arias, F., ?Design of Steel-Composite
Multirim Cylindrical Flywheels Manufactured by Winding with High Tensioning and In
Situ Curing 2. Numerical Analysis,? Mechanics of Composite Materials, Vol. 41, No. 3,
2005.
Potter, K.D. and Medlicott, P.A.C., ?Development of a composite flywheel rotor for
vehicle applications ? a study of the interaction between design, materials and
175
fabrication?, Materials Science Monographs, High Tech - The Way Into the Nineties,
Proceedings of the 7th International Conference of the Society for the Advancement of
Material and Process Engineering European Chapter, 1986.
Prabhakar, S., Sekhar, A. S. and Mohanty, A. R., ?Detection and Monitoring of Cracks in
a Rotor-Bearing System Using Wavelet Transforms,? Mechanical Systems and Signal
Processing, Vol. 15, No. 2, pp. 447-450, 2001.
Qian, S., ?Introduction to Time-Frequency and Wavelet Transforms,? Prentice Hall PTR.
Upper Saddle River, 2002.
Quek, S. T., Wang, Q., Zhang, L. and Ang, K.-K., ?Sensitivity Analysis of Crack
Detection in Beams by Wavelet Technique,? International Journal of Mechanical
Sciences, Vol. 43, No. 12, pp. 2899-2910, 2001.
Sahin, F. , Tuckey, A. M. and Vandenput, A. J. A., ?Design, Development and Testing of
a High-speed Axial-flux Permanent-magnet Machine?, 36th IAS Annual Meeting -
Conference Record of the 2001 Industry Applications, Sep. 30-Oct. 4, 2001, Chicago, IL.
Sandgren, E. and Ragsdell K. M., ?Optimal Flywheel Design with a General Thickness
Form Representation,? ASME Journal of Mechanisms, Transmission, and Automation in
Design, Vol. 105, pp. 425-433, 1983.
176
Sekhar, A. S., and Prabhu, B. S., ?Condition Monitoring of Cracked Rotors through
Transient Response,? Mechanism and Machine Theory, Vol. 33, No. 8, pp. 1167-1175
1998.
Sekhar, A. S, ?Crack Detection and Monitoring in Rotor Supported on Fluid Film
Bearings: Start-up vs. Run-Down,? Mechanical Systems and Signal Processing, Vol. 17,
No. 4, pp. 897-901, 2003.
Shao, H., Jin, W. and Qian, S., ?Order Tracking by Discrete Gabor Expansion,? IEEE
Transaction on Instrumentation and Measurement, Vol. 52, No. 3, pp. 754-761, 2003.
Shen, S. and Veldpaus, F. E., ?Analysis and Control of a Flywheel Hybrid Vehicular
Powertrain?, IEEE Transaction on Control Systems Technology, Vol. 12, No. 5, 2004.
Shen, J. Y. and Fabien, B. C., ?Optimal Control of a Flywheel Energy Storage System
with a Radial Flux Hybrid Magnetic Bearing?, Journal of the Franklin Institute, Vol.
339, pp. 189-210, 2002.
Simpson, W. A. Jr. and McClung, R. W., ?Ultrasonic detection of fatigue damage in
glass-epoxy composite flywheels?, ASTM Special Technical Publication, International
Symposium on Damage Detection and Quality Assurance in Composite Materials, Nov
13-14, 1992.
177
Sinou, J., Villa, C., and Demailly, D., ?Rotordynamics Analysis Experimental and
Numerical Investigations,? Proceedings of ASME 2003 Design Engineering Technical
Conference, Chicago, September 2-6, 2003.
Swett, D. W. and Blanche IV, J. G., ?Flywheel Chargong Module for Energy Storage
Used in Electromagnetic Aircraft Launch System?, IEEE Transactions on Magnetics,
Vol. 41, No. 1, 2005.
Takahashi, I., ?Vibration and Stability of a Cracked Shaft Simultaneously Subjected to a
Follower Force with an Axial Force,? International Journal of Solids Structures, Vol. 35,
No. 23, pp. 3071-3080, 1998.
Thielman, S., and Fabien, B. C., ?An Optimal Control Approach to the Design of
Stacked-ply Composite Flywheels?, Engineering Computations, Vol. 17, No. 5, pp. 541-
555, 2000.
Truong, L. V., Wolff, F. J. and Dravid, N. V., ?Simulation of flywheel-electrical system
for aerospace applications?, Proceedings of the Intersociety Energy Conversion
Engineering Conference, 35th Intersociety Energy Converson Engineering Conference,
Las Vegas, NV, USA, Jul 24-Jul 28, 2000.
Tsai, T. C., and Wang, Y. Z., ?Vibration Analysis and Diagnosis of a Cracked Shaft,?
Journal of Sound and Vibration,? Vol. 192 No. 3, pp. 607-620, 1996.
178
Tsai, T. C., and Wang, Y. Z., ?The Vibration of a Multi-crack Rotor,? International
Journal of Mechanical Science,? Vol. 39, No. 9, pp. 1037-1053, 1997.
Tsai, T. C. and Wang, Y. Z., ?The Vibration of a Rotor with a Transverse Open Crack,?
Proceedings of the National Science Council, ROC (A), Vol. 22, No. 3, pp. 372-384,
1998.
Wu, M. C. and Huang, S. C., ?Vibration and Crack Detection of a Rotor with Speed-
dependent Bearing?, International Journal of Mechanical Science, Vol. 40, No. 6, pp.
545-555, 1998.
Zhang, J., Banks, S. P. and Alleyne, H., ?Optimal Attitude Control for Three-axis
Stabilized Flexible Spacecraft?, Acta Astronautica, Vol. 56, pp. 519-528, 2005.
Zheng, G.T., and Leung, A. Y., ?On the Stability of Cracked Rotor Systems,? Journal of
the Chinese Society of Mechanical Engineering, Vol. 19, No. 1, pp. 125-133, 1997.
179
APPENDICES
180
c***********************************************************************
c APPENDIX A
c
c C++ CODE FOR ?CALCULATIONS OF HUB GEOMETRY AND MATERIAL PROPERTIES?
c
c***********************************************************************
// the file used for parameters output
// Nov. 11 2002
#include
#include
#include
#include
#define MY_PI 3.141592654
#define NUM_ELEMENTS 40
#define NUM_LAYERS 7
void calculate_geometry(float);
void material_properties(int, double );
void uni();
void matinv(double *);
void LeastSquare(double x[500],double y[500],double a[20],int n,int m);
void CurveFitting(double x[500],double y[500],int n,int m);
double
alpha_d[NUM_LAYERS][NUM_ELEMENTS],dz[NUM_LAYERS][NUM_ELEMENTS],tz[NUM_LA
YERS][NUM_ELEMENTS],location[NUM_LAYERS][NUM_ELEMENTS];
double ef1,ef2,gf,nuf,em,gm,num,vf,vm,e1,e2,g12,nu12,nu21,alpha,
a_over_b;
double
k,Cl[3][3],Cg[NUM_LAYERS][3][3],c,c2,c3,c4,s,s2,s3,s4,Ccomposite[3][3];
//--------------------H BEGIN---------------------------------------
int H_size = 0;
double H_len[500], H_tz[500], H_e1[500], H_e2[500];
double H_g12[500], H_nu12[500];
double H_tc,H_dc,H_tcn[100];
double H_xi[500],H_yi[500],H_xo[500],H_yo[500];
double XO[500], tmpXO[500],tmptz[500];
float aob[100],M_Stress[100],M_Uy[100];
float speed[100];
int number = 0;
int tnum;
float tc=0;
//--------------------H END---------------------------------------
FILE *outfile;
int main(int argc, char *argv[])
{
char comm[50];
181
char filename[50];
//--------------------H BEGIN-------------------------------------
//open output file
outfile = fopen("out.txt","w+");
// THE INTERVAL OF a_over_b
for (int aa =1;aa<=8;aa++)
{
double initial=1.99;
a_over_b=initial+0.01*(aa-1);
//a_over_b=1.8;
tnum=(int)(a_over_b*100.0+0.00000000000001)+1;
//tnum=(int)(a_over_b*100.0);
printf("DDDDDD TNUM=%d a_over_b=%f ",tnum,a_over_b);
//--------------------H END--------------------------------------
uni();
tc=0.8;
calculate_geometry(tc);
//printf("Before Curve Fitting\n");
//int cf=(int)H_size*0.6;
for(int ia=1;ia<=H_size;ia++)
{
//fprintf(outfile,"XO[%d]=%d H_tz[%d]=%f\n",ia,ia,ia,H_tz[ia-1]);
XO[ia]=(double)ia;
}
//Curve Fitting Out Side
CurveFitting(XO,H_tz,H_size,2);
//fprintf(outfile,"\nAfter Curve Fitting\n");
fprintf(outfile,"\n");
for(int i=1;i<=H_size;i++)
{
// fprintf(outfile,"\n");
fprintf(outfile,"tz%d(%d)=%f\n",tnum,i,H_tz[i-1]/1000.);
//fprintf(outfile,"\n");
}
sprintf(filename,"%d",tnum);
182
strcpy(comm, "copy out.txt data");
strcat(comm, filename);
strcat(comm, ".txt");
system(comm);
}
//close output file
fclose(outfile);
return 0;
}
//------------------------H FUNCTION------------
// NAME : LeastSquare |
//DESCREPTION : LeastSquare Method |
// DATE : 07/13/2002 |
//------------------------------------------------
void LeastSquare(double x[500],double y[500],double a[20],int n,int m)
{
double s[20],t[20],b[20];
double dt1,dt2,dt3,z,d1,p,c,d2,g,q,dt;
int i,j,k;
for(i=1;i<=m;i++)
{
a[i]=0.0;
}
if(m>n) m=n;
if(m>20) m=20;
z=0.0;
for(i=1;i<=n;i++)
{
z=z+x[i]/n;
}
b[1]=1.0;
d1=n;
p=0.0;
c=0.0;
for(i=1;i<=n;i++)
{
p=p+(x[i]-z);
c=c+y[i];
}
c=c/d1;
p=p/d1;
a[1]=c*b[1];
if(m>1)
{
183
t[2]=1.0;
t[1]=-p;
d2=0.0;
c=0.0;
g=0.0;
for(i=1;i<=n;i++)
{
q=x[i]-z-p;
d2=d2+q*q;
c=y[i]*q+c;
g=(x[i]-z)*q*q+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[2]=c*t[2];
a[1]=c*t[1]+a[1];
}
for(j=3;j<=m;j++)
{
s[j]=t[j-1];
s[j-1]=-p*t[j-1]+t[j-2];
if(j>=4)
{
for(k=j-2;k>=2;k--)
{
s[k]=-p*t[k]+t[k-1]-q*b[k];
}
}
s[1]=-p*t[1]-q*b[1];
d2=0.0;
c=0.0;
g=0.0;
for(i=1;i<=n;i++)
{
q=s[j];
for (k=j-1;k>=1;k--)
{
q=q*(x[i]-z)+s[k];
}
d2=d2+q*q;
c=y[i]*q+c;
g=(x[i]-z)*q*q+g;
184
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[j]=c*s[j];
t[j]=s[j];
for (k=j-1;k>=1;k--)
{
a[k]=c*s[k]+a[k];
b[k]=t[k];
t[k]=s[k];
}
}
dt1=0.0;
dt2=0.0;
dt3=0.0;
for (i=1;i<=n;i++)
{
q=a[m];
for (k=m-1;k>=1;k--)
{
q=q*(x[i]-z)+a[k];
}
dt=q-y[i];
if (fabs(dt)>dt3) dt3=fabs(dt);
dt1=dt1+dt*dt;
dt2=dt2+fabs(dt);
}
return;
}
//------------------------H FUNCTION------------
// NAME : CurveFitting |
//DESCREPTION : Curve Fitting LeastSquare Method |
// DATE : 07/13/2002 |
//------------------------------------------------
void CurveFitting(double x[500],double y[500],int n,int m)
{
double z;
int mi,ti,i,j,k;
double tx=0.4;
double tmp,df,min;
double a[20];
double xbk[500],ybk[500],xx[500],yy[500],yyy[500];
185
printf("Enter Curvefitting");
//Make a copy
for(i=1;i<=n;i++)
{
xbk[i]=x[i];
ybk[i]=y[i];
}
for(i=1;i<=n;i++)
{
xx[i]=x[i]*x[i];
yy[i]=y[i]*y[i];
}
z=0.0;
for(i=1;i<=n;i++)
{
z=z+xx[i]/(double)n;
}
LeastSquare(xx,yy,a,n,m);
for(i=1;i<=m;i++)
{
printf(" A[%d] = %f\n",i,a[i]);
}
for(k=1;k<=n;k++)
{
tx=xx[k];
df=tx-z;
tmp=a[1];
for (j=1; j<=m-1; j++)
{
tmp=tmp+a[j+1]*pow(df,j);
}
yyy[k]=sqrt(fabs(tmp));
//H_yo[k]=sqrt(fabs(tmp));
}
ti=(int)(0.75*n);
min=10000000.0;
mi=0;
for(i=1;i<=ti;i++)
186
{
tmp=fabs(yyy[i]-ybk[i]);
//fprintf(outfile,"tmp=%f\n",tmp);
if(tmp 1)
return( 0)*/
;
if (fabs(max) < fabs(mtxout[j][k]))
{
row = j;
colum = k;
max = mtxout[j][k];
}
}
}
/* Row intercahnge */
ipivot[colum] += 1;
if (row != colum) {
for (l = 0; l < 3; l++) {
max = mtxout[row][l];
mtxout[row][l] = mtxout[colum][l];
mtxout[colum][l] = max;
}
}
index[i][0] = row;
192
index[i][1] = colum;
pivot[i] = mtxout[colum][colum];
mtxout[colum][colum] = 1.0;
for (l = 0; l < 3; l++)
mtxout[colum][l] /= pivot[i];
for (j =0; j < 3; j++)
if (j != colum) {
max = mtxout[j][colum];
mtxout[j][colum] = 0.0;
for (l=0; l<3; l++)
mtxout[j][l] -= mtxout[colum][l] * max;
}
}
for (i = 0; i < 3; i++) {
l = 3 - 1 - i;
if (index[l][0] != index[l][1]) {
row = index[l][0];
colum = index[l][1];
for (k=0; k<3; k++) {
max = mtxout[k][row];
mtxout[k][row] = mtxout[k][colum];
mtxout[k][colum] = max;
}
}
}
for (i=0;i<3;i++)
for (j=0;j<3;j++)
*(ptr+i*3+j) = mtxout[i][j];
}
193
!********************************************************************
!* APPENDIX B
!*
!* ANSYS CODE FOR ?CONSTRUCTION AND OPTIMIZATION OF
!* COMPOSITE SHAFT-HUB-RIM ANSYS MODEL?
!*
!*********************************************************************
!*
!*
/filnam, hubdes
be=2.000
nub=31
*dim,boo,array,nub
*do,i,1,nub,1
boo(i)=1.80+(i-1)*0.01
*enddo
*do,i,1,nub,
*if,be,ge,boo(i),then
*if,be,lt,boo(i+1),then
*if,(be-boo(i)),gt,(boo(i+1)-be),then
tempi=i+1
*else
tempi=i
*endif
tepb=boo(tempi)
*endif
*endif
*enddo
be=tepb
/title, optimization of the hubdesign a_ober_b=%be%
*afun,rad
tnum=100*tepb
be%tnum%=tepb
/input,output,txt
/prep7
w=4000
pi=3.1415926
den=1.52e3
exo=0.
const=3600000.0
CON=1/10
CON1=20
!*element type
!*
194
ET,1,PLANE42
KEYOPT,1,1,1
KEYOPT,1,2,0
KEYOPT,1,3,1
KEYOPT,1,5,0
KEYOPT,1,6,0
!*
!*material properties
!*
*do,i,1,np%tnum%,1
uimp,i,ex,,,mex%tnum%(i),
uimp,i,ey,,,mey%tnum%(i),
uimp,i,ez,,,mey%tnum%(i),
uimp,i,prxy,,,mnuxy%tnum%(i),
uimp,i,prxz,,,0.2,
uimp,i,pryz,,,0.2,
uimp,i,gxy,,,mgxy%tnum%(i),
UIMP,i,GYZ,,,16.0317e9,
UIMP,i,GXZ,,,16.0317e9,
uimp,i,dens,,,den,
*enddo
!*
UIMP,99,EX,,,40.4e9,
UIMP,99,EY,,,143.6e9,
UIMP,99,EZ,,,143.6e9,
UIMP,99,PRXY,,,0.26,
UIMP,99,PRYZ,,,0.26,
UIMP,99,PRXZ,,,0.26,
UIMP,99,GXY,,,56984126984.126984,
UIMP,99,GXZ,,,16.0317e9,
UIMP,99,GXZ,,,16.0317e9,
UIMP,99,DENS,,,1.52e3,
!=============================================================
UIMP,100,EX,,,200e9,
UIMP,100,PRXY,,,0.3,
UIMP,100,GXY,,,76923076923.076920,
UIMP,100,DENS,,,7860,
!=============================================================
!*area properties
!*
*do,i,1,np%tnum%,1
am%tnum%(i)=i
*enddo
!*
!*create the gemeotry of the model
!*
rs=0.01905 !radius of the shaft
ho=0.0250 !height of half flat part
tz0=tz%tnum%(1) !thickness of the flat part
li=0.115-rs-tz0 !distance
195
hi=li/be%tnum%
!*create the coordinates of the keypoints
!*
xi%tnum%(1)=li
yi%tnum%(1)=0
temp1=hi*hi
*do,i,1,np%tnum%
d1=xi%tnum%(i)*xi%tnum%(i)
d2=yi%tnum%(i)*yi%tnum%(i)
ra%tnum%(i)=sqrt(d1+d2)
d3=(ra%tnum%(i)+tz%tnum%(i))/ra%tnum%(i)
xo%tnum%(i)=d3*xi%tnum%(i)
yo%tnum%(i)=d3*yi%tnum%(i)
yi%tnum%(i+1)=yi%tnum%(i)+len%tnum%(i+1)-len%tnum%(i)
temp2=yi%tnum%(i+1)*yi%tnum%(i+1)
xi%tnum%(i+1)=sqrt(abs(1.0-temp2/temp1))*li
*enddo
d1n=xi%tnum%(np%tnum%+1)*xi%tnum%(np%tnum%+1)
d2n=yi%tnum%(np%tnum%+1)*yi%tnum%(np%tnum%+1)
ra%tnum%(np%tnum%+1)=sqrt(d1n+d2n)
d3n=(ra%tnum%(np%tnum%+1)+tz%tnum%(np%tnum%))/ra%tnum%(np%tnum%+1)
xo%tnum%(np%tnum%+1)=d3*xi%tnum%(np%tnum%+1)
yo%tnum%(np%tnum%+1)=d3*yi%tnum%(np%tnum%+1)
!*coordinates transfer
!*
*do,i,1,np%tnum%+1,1
xir%tnum%(i)=xi%tnum%(i)+rs
xor%tnum%(i)=xo%tnum%(i)+rs
yir%tnum%(i)=yi%tnum%(i)+ho
yor%tnum%(i)=yo%tnum%(i)+ho
*enddo
!*create the keypoints
!*
*do,i,1,np%tnum%+1,1
k,i,xir%tnum%(i),yir%tnum%(i)
k,i+100,xor%tnum%(i),yor%tnum%(i)
*enddo
!*creat area
!*
J=1
*do,i,1,np%tnum%,1
196
A,I+1+(J-1)*100,I+(J-1)*100,100*J+I,100*J+I+1
*enddo
!*
!*assign the material properties to area
*do,j,1,np%tnum%,1
ASEL,ALL
VMIN=(J-1)+1
VMAX=J
ASEL,S,,,VMIN,VMAX,1
AATT, am%tnum%(J), , 1, 0
*enddo
!*
ASEL,ALL
MSHAPE,0,2d
MSHKEY,1
AMESH,ALL
!*
/NUMBER,1
!*
!*create the flat part
!*
k,998,li+rs,0
k,999,0.115,0
k,997,0.165,0
k,996,0.165,ho
*do,i,np%tnum%+1,np%tnum%+1
a,998,999,101,1
asel,all
aatt,am%tnum%(1),,1,0
MSHAPE,0,2d
MSHKEY,1
AMESH,ALL
*enddo
!*
!*creat the rim
a,999,997,996,101
asel,all
aatt,am%tnum%(1),,1,0
MSHAPE,0,2d
MSHKEY,1
AMESH,ALL
!*
!
!*creat the shaft
k,1000,0,0
k,1001,rs,0
k,1002,0,0.08
k,1003,rs,0.08
k,1004,0,yir%tnum%(np%tnum%+1)
k,1005,0,yor%tnum%(np%tnum%+1)
!*
a,1000,1001,np%tnum%+1,1004
197
a,np%tnum%+1,1004,1005,100+np%tnum%+1
a,1005,100+np%tnum%+1,1003,1002
asel,all
aatt,99,,1,0
MSHAPE,0,2d
MSHKEY,1
AMESH,ALL
EPLOT
!
PIVCHECK,OFF
finish
!*solution
!*
/SOLU
!*
lsel,s,loc,x,0
dl,all,,ux,0
dl,all,,uy,0
allsel,all
lsel,s,loc,y,0
dl,all,,symm
!dl,all,,ux,0
allsel,all
OMEGA,0,w,0,0
/STATUS,SOLU
SOLVE
FINISH
!*postprocessing
!*
/POST1
SET,FIRST
PLDISP,2
PLNSOL,U,X,0,1
PLNSOL,S,EQV
nsort,s,eqv
*get,smax,SORT,,MAX
etable,evol,volu
!etable,engd,engrd
ssum
*get,vtot,ssum,,item,evol
!vtot=CON-vtot
!enged=CON1-enged
finish
SAVE
198
!*******************************************************
!*
!* OPTIMIZATION OF SHAFT-HUB-RIM ANSYS MODEL
!*
!*****************************************************
!*
!*
!**********************************************************
!* Second Pass: Create optimization input.
!**********************************************************
anfile='hubde02'
anext='inp'
!
!resume,hubdesign.db
! ENTER OPT AND IDENTIFY ANALYSIS FILE
/opt
opanl,anfile,anext
!
! IDENTIFY OPTIMIZATION VARIABLES
opvar,be,dv,1.95,2.05 ! DVs: outter radius
opvar,smax,sv,,0.45e9 ! SV: Maximum equivalent stress
!opvar,vtot,obj,,,1 ! OBJ: Total volume, tolerance = 1.0
opvar,engv,obj,,,1
!
! RUN THE OPTIMIZATION
opkeep,on ! Save best design
optype,subp ! Subproblem approximation method
opsave,anfile,opt0 ! Save the current opt database
opexe
!
! REVIEW RESULTS
oplist,all,,,1 ! List all design sets
plvaropt,smax ! SV smax vs. set number
!plvaropt,enged ! OBJ vtot vs. set number
plvaropt,vtot
xvaropt,be
plvaropt,engd
finish
199
C**********************************************************************
C APPENDIX C
C
C FORTRAN CODE FOR ?CALCULATION OF MAXIMUM ENERGY SYORING DENSITIY?
C XIN HAI
C
C 11/23/2003
C**********************************************************************
REAL XI(1000), YI(1000), XO(1000), YO(1000)
C DOUBLE PRECISION H_TOP, H_BOT, TMP_IN, TMP_OUT
C DOUBLE PRECISION X_CUR, W, STEP,X_SHAFT_OUTER
C DOUBLE PRECISION X_RIM_INNER, X_RIM_OUTER
C DOUBLE PRECISION V_TOT, V_HUB, V_SHAFT, V_RIM
C DOUBLE PRECISION ENG_TOT, ENG_HUB, ENG_SHAFT, ENG_RIM, ENG_DESITY
C FILE FOR X VTALUE FOR INNER SIDE
OPEN(62,FILE='XI.TXT',STATUS='UNKNOWN')
C FILE FOR Y VALUE FOR INNER SIDE
OPEN(65,FILE='YI.TXT',STATUS='UNKNOWN')
C=============================================
C FILE FOR X VALUE FOR OUTER SIDE
OPEN(72,FILE='XO.TXT',STATUS='UNKNOWN')
C FILE FOR Y VALUE FOR OUTER SIDE
OPEN(75,FILE='YO.TXT',STATUS='UNKNOWN')
C FILE FOR ROTATING SPEED
OPEN(80,FILE='W.TXT',STATUS='UNKNOWN')
C=============================================
OPEN(85,FILE='RESULT.DAT',STATUS='UNKNOWN')
C=============================================
C READ ROTATING SPEED FIRST
READ(80,*)W
C READ XI, YI, XO, YO HERE
C=============================================
C READ XI
C=============================================
NPOINT=0
I=0
262 READ(62,*,END=162)TEMP
I = I + 1
NPOINT = NPOINT + 1
XI(I) = TEMP/1000.0
GO TO 262
162 CONTINUE
200
C=============================================
C READ YI
C=============================================
I = 0
265 READ(65,*,END=165)TEMP
I = I + 1
YI(I) = TEMP/1000.0
GO TO 265
165 CONTINUE
IF (I .NE. NPOINT) THEN
WRITE(85,*)'DATA IS WRONG!!!'
WRITE(85,*)'NUMBER OF POINT IN XI.DAT IS NOT EQUAL TO YI.DAT'
END IF
C=============================================
C READ XO
C=============================================
I = 0
272 READ(72,*,END=172)TEMP
I = I + 1
XO(I) = TEMP/1000.0
GO TO 272
172 CONTINUE
IF (I .NE. NPOINT) THEN
WRITE(85,*)'DATA IS WRONG!!!'
WRITE(85,*)'NUMBER OF POINT IN XI.DAT IS NOT EQUAL TO XO.DAT'
END IF
C=============================================
C READ YO
C=============================================
I = 0
275 READ(75,*,END=175)TEMP
I = I + 1
YO(I) = TEMP/1000.0
GO TO 275
175 CONTINUE
IF (I .NE. NPOINT) THEN
WRITE(85,*)'DATA IS WRONG!!!'
WRITE(85,*)'NUMBER OF POINT IN XI.DAT IS NOT EQUAL TO YO.DAT'
END IF
C DEFINE CONSTANTS
PI = 3.1415926
FLAT_THICKNESS=6.0/1000.0
201
C CALCULATE THE ENERGY DENSITY FOR HUB
NSTEP = 200
X_SHAFT_OUTER=19.05/1000.0
X_RIM_INNER=115.0/1000.0
X_RIM_OUTER=165.0/1000.0
STEP = (X_RIM_INNER-X_SHAFT_OUTER)/NSTEP
C MATERIAL PROPERTIES
DEN_SHAFT=7860.0
DEN_RIM=1520.0
DEN_HUB=1520.0
C ASSIGN INITIAL VALUES
V_HUB = 0.0
ENG_HUB = 0.0
DO I = 1, NSTEP
TMP_IN = STEP*(I-1)
X_CUR=X_RIM_INNER-TMP_IN
CALL INTERPOLATE(NPOINT,XI(1),YI(1),X_CUR,H_BOT)
CALL INTERPOLATE(NPOINT,XO(1),YO(1),X_CUR,H_TOP)
V_HUB =V_HUB+2.0*PI*X_CUR**2.0*(H_TOP-H_BOT)-
& 2.0*PI*(X_CUR-STEP)**2.0*(H_TOP-H_BOT)
TMP_IN = (X_CUR-STEP)**4.0
TMP_OUT = X_CUR**4.0
ENG_HUB = ENG_HUB + 2.0*PI*DEN_HUB*(H_TOP-H_BOT)*
& (TMP_OUT-TMP_IN)*W*W/(4.0*3600000.0)
END DO
C ADD FLAT PART OF HUB
V_HUB = V_HUB + 2.0*PI*(X_RIM_INNER**2.0 -
& (X_RIM_INNER-FLAT_THICKNESS)**2.0)*25.0/1000.0
TMP_IN = (X_RIM_INNER-FLAT_THICKNESS)**4.0
TMP_OUT = X_RIM_INNER**4.0
ENG_ADD = 2.0*PI*DEN_HUB*25.0/1000.0*
& (TMP_OUT-TMP_IN)*W*W/(4.0*3600000.0)
ENG_HUB = ENG_HUB + 2.0*PI*DEN_HUB*25.0/1000.0*
& (TMP_OUT-TMP_IN)*W*W/(4.0*3600000.0)
C CALCULATE ENERGY AND VOLUME FOR SHAFT
V_SHAFT = PI*X_SHAFT_OUTER**2.0*100.0/1000.0
ENG_SHAFT = PI*DEN_SHAFT*100.0/1000.0*X_SHAFT_OUTER**4.0*W*W
& /(4.0*3600000.0)
C CALCULATE ENERGY AND VOLUME FOR RIM
V_RIM = PI*(X_RIM_OUTER**2.0-X_RIM_INNER**2.0)*50.0/1000.0
202
ENG_RIM = PI*DEN_RIM*50.0/1000.0*(X_RIM_OUTER**4.0-
& X_RIM_INNER**4.0)*W*W/(4.0*3600000.0)
V_TOT = V_HUB + V_SHAFT + V_RIM
ENG_TOT = ENG_HUB + ENG_SHAFT + ENG_RIM
ENG_DESITY = ENG_TOT/(9.8*DEN_SHAFT*V_SHAFT + 9.8*DEN_HUB*V_HUB +
& 9.8*DEN_RIM*V_RIM)
WRITE(85,*)'VOLUME OF SHAFT = ',V_SHAFT,' m^3'
WRITE(85,*)'VOLUME OF HUB = ',V_HUB,' m^3'
WRITE(85,*)'VOLUME OF RIM = ',V_RIM,' m^3'
WRITE(85,*)'TOTAL VOLUME = ',V_TOT,' m^3'
WRITE(85,*)'ENERGY OF SHAFT = ',ENG_SHAFT,' Kwhr'
WRITE(85,*)'ENERGY OF HUB = ',ENG_HUB,' Kwhr'
WRITE(85,*)'ENERGY OF RIM = ',ENG_RIM,' Kwhr'
WRITE(85,*)'ENG_DESITY = ',ENG_DESITY
WRITE(85,*)'PROGRAM FINSHED NORMALLY'
STOP
END
C====================================================
C LINEAR INTERPOLATION OF Y ACCORDING TO GIVEN X
C====================================================
SUBROUTINE INTERPOLATE(NPOINT,X,Y,X_IN_ORG,Y_OUT)
REAL X(NPOINT),Y(NPOINT)
C DOUBLE PRECISION X_IN,Y_OUT
C FIND THE POSITION OF X_IN IN THE TABLE
X_IN = X_IN_ORG
IF (X_IN_ORG .LE. X(NPOINT)) X_IN = X(NPOINT)
IF (X_IN_ORG .GE. X(1)) X_IN = X(1)
INDEX=0
DO I=1, NPOINT-1
IF((X_IN .LE. X(I)) .AND. (X_IN .GE. X(I+1)) ) THEN
INDEX = I
GOTO 88
END IF
END DO
88 CONTINUE
IF(INDEX .EQ. 0) THEN
WRITE(*,*) 'PROGRAM FIND WRONG INDEX!!! STOP'
STOP
ELSE
Y_OUT= (Y(I+1)-Y(I))/(X(I+1)-X(I))*(X_IN-X(I))+Y(I)
END IF
END
203
c***********************************************************************
c APPENDIX D
c
c FORTRAN CODE FOR ?IMBALANCE ECCENTRICITY OF THE MASS CENTER
c CALCULATIONS?
c
c XIN HAI
c
c***********************************************************************
c
c program to calculate the imbalance eccentricity of mass
c center for crack model of the rotational rim
c
c***********************************************************************
program main
c n=number of nodes,nn=number of hoop-crack nodes
c nrn=number of radial-crack nodes
parameter (n=22002,nn=0,nrn=36,n3p=0)
common /var/ ldown(nn),lup(nn),id(n),ldnr(nrn)
common /vv/ h(100),Idin(n3p,3),Idout(n3p,3)
c*****************************************************
c
integer itn,tid
real x(n),y(n),z(n)
real mass(n),hold
real ux(n),uy(n),uz(n),tux,tuy,tuz
double precision am,tm
real itx,ity
real itz
open(unit=11,file='rad-radnode.lis',status='OLD')
open(unit=12,file='rr6283.lis',status='OLD')
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
open(unit=113,file='laydownhh.lis',status='OLD')
open(unit=114,file='layuphh.lis',status='OLD')
open(unit=115,file='laynr.lis',status='OLD')
open(unit=116,file='Inpoints.lis',status='OLD')
open(unit=117,file='Outpoints.lis',status='OLD')
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
open(unit=21,file='nmass.dat',status='UNKNOWN')
open(unit=22,file='unbalance.dat',status='UNKNOWN')
c read the coordinates of node
itn=0
ih=0
i=1
22 read(11,*,end=23)id(i),x(i),y(i),z(i)
itn=itn+1
204
c-----save all different height levels in array h
ihflag=0
do j=1,ih
if(abs((z(i)-h(j))).le.0.0000000001) then
ihflag=1
endif
enddo
c-----save new height level
if(ihflag.eq.0) then
ih=ih+1
h(ih)=z(i)
endif
i=i+1
go to 22
c write(6,*)itn
23 continue
c-----sort array h in ascendent order
do 288 i=1,ih
do 288 j=1,ih-1
if(h(j).gt.h(j+1)) then
hold=h(j)
h(j)=h(j+1)
h(j+1)=hold
end if
288 continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c get the hoop-crack node number
do j=1,nn
read(113,*) ldown(j)
read(114,*) lup(j)
end do
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c get the radial-crack node number
do j=1,nrn
read(115,*) ldnr(j)
end do
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c----get id number of 3 points group
do j=1,n3p
205
read(116,*) Idin(j,1),Idin(j,2),Idin(j,3)
read(117,*) Idout(j,1),Idout(j,2),Idout(j,3)
enddo
c calculate the mass distribution of each node
tm=0.0
do i=1,itn
xx=x(i)
yy=y(i)
zz=z(i)
am=0.0
idd=id(i)
call massdt(xx,yy,zz,am,idd,ih)
mass(i)=am
tm=tm+am
write(21,*) id(i),mass(i)
end do
write(6,*) tm
c read the deformations of node
25 read(12,*,end=26) tid,tux,tuy,tuz
iflag=0
do j=1,itn
if(id(j).eq.tid) then
ux(j)=tux
uy(j)=tuy
uz(j)=tuz
iflag=1
endif
enddo
if(iflag.eq.0) write(6,*)tid,' WARNING: POINTS NOT FIT!'
go to 25
26 continue
c
c calculate the imbalance eccentricity
itx=0.0
ity=0.0
itz=0.0
xcm=0.0
ycm=0.0
zcm=0.0
rcm=0.0
do i=1,itn
uxx=ux(i)
uyy=uy(i)
uzz=uz(i)
um=mass(i)
itx=itx+uxx*um
ity=ity+uyy*um
206
itz=itz+uzz*um
end do
xcm=itx/tm
ycm=ity/tm
zcm=itz/tm
rcm=sqrt(xcm**2+ycm**2+zcm**2)
write(6,*) xcm,ycm,zcm,rcm
write(22,*) xcm,ycm,zcm,rcm
stop
end
c*****************************************************
c
c
c*****************************************************
c
c subroutine to calculate the mass
c
c*****************************************************
subroutine massdt(xx,yy,zz,am,idd,ih)
c-----------------------------------------------------
parameter (n=22002,nn=0,nrn=36,n3p=0)
common /var/ ldown(nn),lup(nn),id(n),ldnr(nrn)
common /vv/ h(100),Idin(n3p,3),Idout(n3p,3)
double precision dhUp,dhDown,massIn,massOut
double precision rr,dr,x,p,r1,r2,am,small
c the number of node in each element layer
parameter (nr=48)
small=0.000001
c ri= inside radius of rim
c ro= outdise radius of rim
ri=0.066
ro=0.165
c zbot= bottom edge of rim
c ztop= top edge of rim
zbot=-0.0255
ztop=0.0255
c increment of r
dr=0.0055
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c get the increment of h
do jj=2,ih-1
207
if(abs(zz-h(jj)).le.small) then
dhUp=(h(jj+1)-h(jj))/2.0
dhDown=(h(jj)-h(jj-1))/2.0
endif
enddo
c-----Modify dhUp and Dhdown value when point is on top or bottom
if(abs(zz-zbot).le.small) then
dhDown=0.0
dhUp=(h(2)-h(1))/2.0
endif
if(abs(zz-ztop).le.small) then
dhUp=0.0
dhDown=(h(ih)-h(ih-1))/2.0
endif
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c material properties
den=1520.0
pi=3.1415926
c find the place of each node
rr=sqrt(xx**2+yy**2)
c-----addjust rr
do i=0,20
if(abs(rr-(ri+i*dr)).le.(0.01*dr)) rr=ri+i*dr
enddo
c-----------------------------------------------
c calculate mass
massIn=den*pi*(dhUp+dhDown)*(rr**2-(rr-dr)**2)/nr
massOut=den*pi*(dhUp+dhDown)*((rr+dr)**2-rr**2)/nr
c--------------------------------------------
c inside node
if (abs(rr-ri).le.(0.01*dr)) then
r1=rr
r2=rr+dr
call CalculateX(r1,r2,p)
x=p
am=x*massOut
c inside edge node
if ((abs(zz-zbot).le.small).or.(abs(zz-ztop).le.small)) then
am=am/1.0
end if
208
c inside non-edge node
return
end if
c----------------------------------------------
c outside node
if (abs(rr-ro).le.(0.01*dr)) then
r1=rr-dr
r2=rr
p=0
call CalculateX(r1,r2,p)
x=p
am=(1.0-x)*massIn
c outside edge node
if ((abs(zz-zbot).le.small).or.(abs(zz-ztop).le.small)) then
am=am/1.0
end if
c outside non-edge node
return
end if
c-------------------------------------------------
c non-side node
if ((abs(rr-ri).gt.(0.01*dr)).and.(abs(rr-ro).gt.(0.01*dr))) then
r1=rr-dr
r2=rr
call CalculateX(r1,r2,p)
x=p
am1=(1.0-x)*massIn
r1=rr
r2=rr+dr
call CalculateX(r1,r2,p)
x=p
am2=(x)*massOut
am=am1+am2
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c modify mass on hoop-crack nodes (laydown, layup)
do j=1,nn
if (idd.eq.ldown(j)) then
am=am1
end if
if (idd.eq.lup(j)) then
am=am2
end if
enddo
209
c-----modify mass on 3 points cases
c-----take the points in Inside points list
c-----treat type 1 points
do j=1,n3p
if(idd.eq.Idin(j,1)) then
am=am1
endif
c-----treat type 2 points
if(idd.eq.Idin(j,2)) then
am=am2/2.0
endif
c-----treat type 3 points
if(idd.eq.Idin(j,3)) then
am=am2/2.0
endif
enddo
c-----take the points in Outside points list
c-----treat type 1 points
do j=1,n3p
if(idd.eq.Idout(j,1)) then
am=am2
endif
c-----treat type 2 points
if(idd.eq.Idout(j,2)) then
am=am1/2.0
endif
c-----treat type 3 points
if(idd.eq.Idout(j,3)) then
am=am1/2.0
endif
enddo
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c modify mass on radial-crack nodes
do j=1,nrn
if (idd.eq.ldnr(j)) then
am=am/2.0
end if
enddo
210
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c non-side and edge node
if ((abs(zz-zbot).le.small).or.(abs(zz-ztop).le.small)) then
am=am/1.0
endif
c non-side and non-edge node
return
end if
end
c************************************************************
c
c
c************************************************************
c subroutine to calculate the portion of distribution
c
subroutine CalculateX(r1,r2,x)
double precision r1,r2,x
x=(2./3.*r1**2-1./3.*r1*r2-1./3.*r2**2)/(r1**2-r2**2)
return
end
c************************************************************
211
%***********************************************************************
% APPENDIX E
%
% MATLAB CODE FOR ?SIMULATION OF CRACKED DISK MODEL?
%
% XIN HAI
%
%***********************************************************************
%balanced case based on support frame equations
function yprime=disk02(t,y)
%damping factor for the support
c=0.8;
%damping factor for the mass
c1=0.5;
c2=0.5;
c3=0.5;
c4=0.5;
%support spring
kx=5.0e7;
ky=5.0e7;
%mass spring
k1=5.0e7;
k2=5.0e7;
k3=5.0e7;
k4=5.0e7;
%mass of the disk and hub
m=2.45;
mh=0.2;
%rotating speed w=0,10,30,100,200,350,500
w=10;
%square of the rotating speed
ww=w*w;
%radius of the hub
R=0.076;
%variables
x=y(1);
yy=y(3);
v(1)=y(5);
v(2)=y(7);
v(3)=y(9);
v(4)=y(11);
xd=y(2);
yd=y(4);
vd(1)=y(6);
vd(2)=y(8);
vd(3)=y(10);
212
vd(4)=y(12);
for i=1:4,
phi(i)=3.14159/2*(i-1);
end
sum1=0;
for i=1:4,
sum1=sum1+vd(i)*sin(w*t+phi(i));
end
sum2=0;
for i=1:4,
sum2=sum2+vd(i)*cos(w*t+phi(i));
end
sum3=0;
for i=1:4,
sum3=sum3+(v(i)+R)*cos(w*t+phi(i));
end
sum4=0;
for i=1:4,
sum4=sum4+(v(i)+R)*sin(w*t+phi(i));
end
yprime = zeros(12,1)';
%variable dX/dt
yprime(1)=y(2);
%variable dY/dt
yprime(3)=y(4);
%variable dV1/dt
yprime(5)=y(6);
%variable dV2/dt
yprime(7)=y(8);
%variable dV3/dt
yprime(9)=y(10);
%variable dV4/dt
yprime(11)=y(12);
%mass matrix
mass=[mh+4*m,0,m*cos(w*t+phi(1)),m*cos(w*t+phi(2)),m*cos(w*t+phi(3)),m*c
os(w*t+phi(4))];
mass=[mass;0,mh+4*m,m*sin(w*t+phi(1)),m*sin(w*t+phi(2)),m*sin(w*t+phi(3)
),m*sin(w*t+phi(4))];
mass=[mass;m*cos(w*t+phi(1)),m*sin(w*t+phi(1)),m,0,0,0];
mass=[mass;m*cos(w*t+phi(2)),m*sin(w*t+phi(2)),0,m,0,0];
mass=[mass;m*cos(w*t+phi(3)),m*sin(w*t+phi(3)),0,0,m,0];
mass=[mass;m*cos(w*t+phi(4)),m*sin(w*t+phi(4)),0,0,0,m];
213
ydd(1)=8*m*w*yd+4*m*ww*x-c*xd+4*m*w*sum1+4*m*ww*sum3-kx*x;
ydd(2)=-8*m*w*xd+4*m*ww*yy-c*yd-4*m*w*sum2+4*m*ww*sum4-ky*yy;
ydd(3)=2*m*w*yd*cos(w*t+phi(1))-
2*m*w*xd*sin(w*t+phi(1))+m*ww*yy*sin(w*t+phi(1))+m*ww*x*cos(w*t+phi(1))+
4*m*ww*(R+v(1))-c1*vd(1)-k1*v(1);
ydd(4)=2*m*w*yd*cos(w*t+phi(2))-
2*m*w*xd*sin(w*t+phi(2))+m*ww*yy*sin(w*t+phi(2))+m*ww*x*cos(w*t+phi(2))+
4*m*ww*(R+v(2))-c2*vd(2)-k2*v(2);
ydd(5)=2*m*w*yd*cos(w*t+phi(3))-
2*m*w*xd*sin(w*t+phi(3))+m*ww*yy*sin(w*t+phi(3))+m*ww*x*cos(w*t+phi(3))+
4*m*ww*(R+v(3))-c3*vd(3)-k3*v(3);
ydd(6)=2*m*w*yd*cos(w*t+phi(4))-
2*m*w*xd*sin(w*t+phi(4))+m*ww*yy*sin(w*t+phi(4))+m*ww*x*cos(w*t+phi(4))+
4*m*ww*(R+v(4))-c4*vd(4)-k4*v(4);
ypp=inv(mass)*ydd';
yprime(2)=ypp(1);
yprime(4)=ypp(2);
yprime(6)=ypp(3);
yprime(8)=ypp(4);
yprime(10)=ypp(5);
yprime(12)=ypp(6);
yprime=yprime';
214
%clear all;
%close all;
tspan=[0,0.5];
for i=1:12,
xo(i)=0.0;
end
%variables
%x=xo(1);
%yy=xo(3);
%v(1)=xo(5);
%v(2)=xo(7);
%v(3)=xo(9);
%v(4)=xo(11);
%xd=xo(2);
%yd=xo(4);
%vd(1)=xo(6);
%vd(2)=xo(8);
%vd(3)=xo(10);
%vd(4)=xo(12);
xo(1)=0.0000001;
xo(2)=0.000000;
xo(3)=0.0000001;
xo(4)=0.000000;
%---------------------------
[t1,y1]=ode45('disk02',tspan,xo);
figure (2)
grid on;
subplot(3,3,1);
plot(t1,y1(:,1));
title('Time Vs. X');
grid on;
subplot(3,3,2);
plot(t1,y1(:,3));
title('Time Vs. Y');
grid on;
subplot(3,3,3);
plot(t1,y1(:,5));
title('Time Vs. V1');
grid on;
subplot(3,3,4);
plot(t1,y1(:,7));
title('Time Vs. V2');
grid on;
215
subplot(3,3,5);
plot(t1,y1(:,9));
title('Time Vs. V3');
grid on;
subplot(3,3,6);
plot(t1,y1(:,11));
title('Time Vs. V4');
grid on;
216
%unbalanced system program based on support frame equations
function yprime=disk03(t,y)
%damping factor for the support
c=0.8;
%damping factor for the mass
c1=0.5;
c2=0.5;
c3=0.5;
c4=0.5;
%support spring
kx=5.0e7;
ky=5.0e7;
% unbalanced mass spring
% change k1 to simulate the crack
% when w=0, k1 can change from 5.0e7 to 5.0e2
% when w=100, k1 can change from 5.0e7 to 2.0e5
% when w=300, k1 can change from 5.0e7 to 2.0e6
% when w=500, k1 can change from 5.0e7 to 4.0e6
%k1=5.0e6;
k1=5.0e5;
%k1=5.0e4;
%k1=5.0e3;
%k1=5.0e2;
%k1=2.0e5;
%k1=2.0e6;
%k1=4.0e6;
k2=5.0e7;
k3=5.0e7;
k4=5.0e7;
%mass of the disk and hub
m=2.45;
mh=0.2;
%rotating speed
% for unbalanced cases, w=0,100,300,500 Hz
%w=0;
w=100;
%w=300;
%w=500;
%square of the rotating speed
ww=w*w;
%radius of the hub
R=0.076;
%variables
x=y(1);
yy=y(3);
v(1)=y(5);
217
v(2)=y(7);
v(3)=y(9);
v(4)=y(11);
xd=y(2);
yd=y(4);
vd(1)=y(6);
vd(2)=y(8);
vd(3)=y(10);
vd(4)=y(12);
for i=1:4,
phi(i)=3.14159/2*(i-1);
end
sum1=0;
for i=1:4,
sum1=sum1+vd(i)*sin(w*t+phi(i));
end
sum2=0;
for i=1:4,
sum2=sum2+vd(i)*cos(w*t+phi(i));
end
sum3=0;
for i=1:4,
sum3=sum3+(v(i)+R)*cos(w*t+phi(i));
end
sum4=0;
for i=1:4,
sum4=sum4+(v(i)+R)*sin(w*t+phi(i));
end
yprime = zeros(12,1)';
%variable dX/dt
yprime(1)=y(2);
%variable dY/dt
yprime(3)=y(4);
%variable dV1/dt
yprime(5)=y(6);
%variable dV2/dt
yprime(7)=y(8);
%variable dV3/dt
yprime(9)=y(10);
%variable dV4/dt
yprime(11)=y(12);
%mass matrix
218
mass=[mh+4*m,0,m*cos(w*t+phi(1)),m*cos(w*t+phi(2)),m*cos(w*t+phi(3)),m*c
os(w*t+phi(4))];
mass=[mass;0,mh+4*m,m*sin(w*t+phi(1)),m*sin(w*t+phi(2)),m*sin(w*t+phi(3)
),m*sin(w*t+phi(4))];
mass=[mass;m*cos(w*t+phi(1)),m*sin(w*t+phi(1)),m,0,0,0];
mass=[mass;m*cos(w*t+phi(2)),m*sin(w*t+phi(2)),0,m,0,0];
mass=[mass;m*cos(w*t+phi(3)),m*sin(w*t+phi(3)),0,0,m,0];
mass=[mass;m*cos(w*t+phi(4)),m*sin(w*t+phi(4)),0,0,0,m];
ydd(1)=8*m*w*yd+4*m*ww*x-c*xd+4*m*w*sum1+4*m*ww*sum3-kx*x;
ydd(2)=-8*m*w*xd+4*m*ww*yy-c*yd-4*m*w*sum2+4*m*ww*sum4-ky*yy;
ydd(3)=2*m*w*yd*cos(w*t+phi(1))-
2*m*w*xd*sin(w*t+phi(1))+m*ww*yy*sin(w*t+phi(1))+m*ww*x*cos(w*t+phi(1))+
4*m*ww*(R+v(1))-c1*vd(1)-k1*v(1);
ydd(4)=2*m*w*yd*cos(w*t+phi(2))-
2*m*w*xd*sin(w*t+phi(2))+m*ww*yy*sin(w*t+phi(2))+m*ww*x*cos(w*t+phi(2))+
4*m*ww*(R+v(2))-c2*vd(2)-k2*v(2);
ydd(5)=2*m*w*yd*cos(w*t+phi(3))-
2*m*w*xd*sin(w*t+phi(3))+m*ww*yy*sin(w*t+phi(3))+m*ww*x*cos(w*t+phi(3))+
4*m*ww*(R+v(3))-c3*vd(3)-k3*v(3);
ydd(6)=2*m*w*yd*cos(w*t+phi(4))-
2*m*w*xd*sin(w*t+phi(4))+m*ww*yy*sin(w*t+phi(4))+m*ww*x*cos(w*t+phi(4))+
4*m*ww*(R+v(4))-c4*vd(4)-k4*v(4);
ypp=inv(mass)*ydd';
yprime(2)=ypp(1);
yprime(4)=ypp(2);
yprime(6)=ypp(3);
yprime(8)=ypp(4);
yprime(10)=ypp(5);
yprime(12)=ypp(6);
yprime=yprime';
219
%clear all;
%close all;
tspan=[0,1];
for i=1:12,
xo(i)=0.0;
end
%variables
%x=xo(1);
%yy=xo(3);
%v(1)=xo(5);
%v(2)=xo(7);
%v(3)=xo(9);
%v(4)=xo(11);
%xd=xo(2);
%yd=xo(4);
%vd(1)=xo(6);
%vd(2)=xo(8);
%vd(3)=xo(10);
%vd(4)=xo(12);
xo(1)=0.0000001;
xo(2)=0.000000;
xo(3)=0.0000001;
xo(4)=0.000000;
%---------------------------
[t1,y1]=ode45('disk03',tspan,xo);
figure (1)
grid on;
subplot(2,2,1);
plot(t1,y1(:,5));
title('Time Vs. V1');
grid on;
subplot(2,2,2);
plot(t1,y1(:,7));
title('Time Vs. V2');
grid on;
subplot(2,2,3);
plot(t1,y1(:,9));
title('Time Vs. V3');
grid on;
subplot(2,2,4);
plot(t1,y1(:,11));
title('Time Vs. V4');
grid on;
220
wol=y1(:,1);
yol=y1(:,3);
v1=y1(:,5);
v2=y1(:,7);
v3=y1(:,9);
v4=y1(:,11);
save xau.txt wol -ascii;
save yau.txt yol -ascii;
save v1u.txt v1 -ascii;
save v2u.txt v2 -ascii;
save v3u.txt v3 -ascii;
save v4u.txt v4 -ascii;
221