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