MODAL ANALYSIS OF NON-DIAGONALIZABLE CONTINUOUS SYSTEMS WITH APPLICATION TO WIND TURBINE BLADES By Xing Xing A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering – Doctor of Philosophy 2015 ABSTRACT MODAL ANALYSIS OF NON-DIAGONALIZABLE CONTINUOUS SYSTEMS WITH APPLICATION TO WIND TURBINE BLADES By Xing Xing This work represents an investigation of the modal analysis of distributed parameter systems whose stiffness or damping terms are non-diagonalizable with an undamped modalcoordinate transformation. The non-diagonalizability may be caused by nonmodal damping or stiffness that includes parametric excitation. The modal properties for these kinds of problems will be investigated. An approach for analyzing the complex modes of continuous systems with nonmodal damping is first developed. As an example, a cantilevered beam with damping at the free end is studied. Assumed modes are applied to discretize the eigenvalue problem in state-variable form, to then obtain estimates of the natural frequencies and state-variable modal vectors. The finite-element method is also used to get the mass, stiffness, and damping matrices for the state-variable eigenvalue problem. A comparison between the complex modes and eigenvalues obtained from the assumed-mode analysis and the finite-element analysis shows that the methods produce consistent results. The assumed-mode method is then used to study the effects of the end-damping coefficient on the estimated normal modes and modal damping. Most modes remain underdamped regardless of the end-damping coefficient. There is an optimal end-damping coefficient for vibration decay, which correlates with the maximum modal nonsynchronicity. As an experimental example of a non-modally damped continuous system, an end-damped cantilevered beam is studied for its complex modal behavior. An eddy-current damper is applied considering its noncontact and linear properties. The state-variable modal decomposition method (SVMD) is applied to extract the modes from impact responses. Characteristics of the mode shapes and modal damping are examined for various values of the damping coefficient. The eigenvalues and mode shapes obtained from the experiments are consistent with the numerical analysis of the model, although there is variation relative to sampling parameters. Over the range of damping coefficients studied in the experiments, we observe a maximum damping ratio in the lowest underdamped mode, which correlates with the maximum modal nonsynchronicity. The vibration model of a horizontal-axis wind turbine blade can be approximated as a rotating pretwisted nonsymmetric beam, with damping and gravitational and aeroelastic loading. The out-of-plane (flapwise) and in-plane (edgewise) motion of a wind turbine blade are examined with simple aeroelastic damping effects. Hamilton’s principle is applied to derive the in-plane and out-of-plane equations of motion, and the partial differential equation is linearized and then discretized by the assumed-mode method. A simple quasi-steady bladeelement airfoil theory is applied to obtain the aeroelastic damping. The analysis is performed on three blades of different size. The effects of nonproportional damping are not strong, but are seen to become more significant as the blade size increases. The results provide some experience for the validity of making modal damping assumptions in blade analyses. A perturbation approach is developed to analyze the perturbation effect of the parametric excitation on the unperturbed linear modes. The method is applied in the examples of a three-mass system and a wind turbine blade. In wind-turbine blades, the parametric excitation has a weak effect on the non-resonant unperturbed linear modal responses. Copyright by XING XING 2015 To someone who saves me from trouble and leads me in the right way, to someone who always has a positive life attitude v ACKNOWLEDGMENTS I sincerely thank Professor Brian Feeny for his guidance, assistance, and friendship during my graduate studies. I want to thank him for supporting me with the opportunity to enter and carry out research in the field of vibration and dynamics. He showed me the right way to do research work by his rigorous acdemic attitude, his enthusiasm and his profound knowledge. His great teaching in nonlinear vibration and chaos is also impressive. I thank Professor Steve Shaw, Seungik Baek, and Kirk Dolan for serving on my doctoral committee. I am grateful to have had the opportunity to learn from each of them. Professor Shaw’s excellent teaching in courses on theory of vibration and advanced dynamics gave me a solid foundation for my studies. His teaching is exceptional. I also have the opportunity to learn continuum mechanics from Professor Baek, which is very useful in my research. I would also like to thank Professor Dolan for helpful discussions about research and academis. I thank the National Science Foundation (CMMI-1335177), Permawick Company and the Carl Coppola Endowment Foundation to provide me with the funding to carry on with my research. I thank my parents for their help and support all through my life. Thank you for raising me up and teaching me a right life attitude. I would also like to thank my girlfriend for her advice in my life. vi TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Chapter 1 Introduction . . . . . . . . . . . . . . . . . . 1.1 Objective and Motivation . . . . . . . . . 1.2 Research Background . . . . . . . . . . . 1.2.1 Nonmodally Damped Distributed Systems . . 1.2.2 Modal Identification Methods of Experimental 1.2.3 Eddy Current Damping . . . . . . . . . . . . 1.2.4 Wind-Turbine Blade Modeling . . . . . . . . . 1.2.5 Aerodynamics of Wind Turbine Blades . . . . 1.2.6 This Work . . . . . . . . . . . . . . . . . . . . 1.3 Contributions . . . . . . . . . . . . . . . Chapter 2 2.1 2.2 2.2.1 2.3 2.3.1 2.3.2 2.3.3 2.3.4 2.4 Chapter 3 3.1 3.2 3.3 3.4 3.4.1 3.4.2 3.4.3 3.4.4 . . . . . . . . . . 1 1 3 3 4 5 6 8 10 12 Complex Modal Analysis of a Non-Modally Damped Continuous Beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter Introduction . . . . . . . . . . . . . . . . . . . . . . . . Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . Assumed-Mode Formulation of a Damped Beam . . . . . . . . . . . . Example : End-Damped Cantilevered Beam . . . . . . . . . . . . Assumed-Mode Formulation . . . . . . . . . . . . . . . . . . . . . . . Convergence Study of Different Assumed Mode Functions . . . . . . . Finite-Element Method . . . . . . . . . . . . . . . . . . . . . . . . . . Interpretation of Beam Modal Properties . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 13 14 14 21 21 21 24 28 33 Experimental Study on Complex Modes of an End Continuous Beam . . . . . . . . . . . . . . . . . . . . . . Chapter Introduction . . . . . . . . . . . . . . . . . . Eddy Current Damping Formulation . . . . . . . . . . State-Variable Modal Decomposition . . . . . . . . . . Cantilevered Beam Experiment . . . . . . . . . . . . . Experimental Setup Description . . . . . . . . . . . . . . . Analysis of the Model . . . . . . . . . . . . . . . . . . . . Modal Identification Results . . . . . . . . . . . . . . . . . Modal Coordinates . . . . . . . . . . . . . . . . . . . . . . 38 38 39 41 42 42 45 47 56 vii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Modal Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Damped . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 . . . . . . . . . . . . . . . . . . . . Chapter 4 Modal Analysis of a Wind Turbine Blade . . . . . . . . . 4.1 Chapter Introduction . . . . . . . . . . . . . . . . . . . . 4.2 Equation of Motion Formulation . . . . . . . . . . . . . . 4.2.1 Edgewise Bending . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Flapwise Bending . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Linearization and Modal Reduction . . . . . . . . . . . . 4.3.1 Edgewise Bending Modal Analysis . . . . . . . . . . . . . . . . 4.3.2 Flapwise Bending Modal Analysis . . . . . . . . . . . . . . . . 4.4 Study of the Effects of Rotational Position due to Gravity 4.5 Nonmodal Aerodynamic Damping . . . . . . . . . . . . . 4.5.1 Flapwise Bending . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.2 Edgewise Bending . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.3 Effect of Equilibrium on Stiffness . . . . . . . . . . . . . . . . 4.5.4 Modal Analysis Results . . . . . . . . . . . . . . . . . . . . . . 4.5.5 Quantification of the Nonmodal Damping . . . . . . . . . . . 4.5.6 A Review of Unsteady Dynamic Stall Lift Model . . . . . . . . 4.5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 . 63 . 64 . 65 . 67 . 68 . 70 . 74 . 79 . 84 . 86 . 91 . 93 . 95 . 103 . 106 . 108 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 110 110 114 116 117 118 Chapter 6 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . 6.1 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1 Complex Modal Analysis of a Nonmodally Damped Continuous Beam 6.1.2 Experimental Study on Complex Modes of an End Damped Continuous Beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.3 Modal Analysis of a Wind Turbine Blade . . . . . . . . . . . . . . . . 6.1.4 Modal Analysis of Parametrically Excited Systems . . . . . . . . . . 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 122 122 3.5.1 3.5.2 3.6 Complex Orthogonal Decomposition Verification Complex Orthogonal Decomposition Theory . . . . . Modal Identification Results . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 5 Modal Analysis of Parametrically Excited Systems 5.1 Chapter Introduction . . . . . . . . . . . . . . . . . 5.2 Perturbation Approach for Nonresonant Response . 5.3 Example: Three Mass System . . . . . . . . . . . . 5.4 Example : Wind Turbine Blade . . . . . . . . . . . . 5.4.1 Quantification of the Nonmodal Stiffness Matrix . . . . . 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 58 60 61 123 124 125 127 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 Appendix A Study of the Effect of the Number of Samples and Filter Frequency on SVMD . . . . . . . . . . . . . . . . . . . . . . . . . . 130 Appendix B Comparison for the axial stiffness and bending stiffness in a rotating beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 viii BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 ix LIST OF TABLES Table 2.1 Eigenvalue comparison from two kinds of assumed-mode functions and finite element method, c = 50 kg/s . . . . . . . . . . . . . . . . 24 MAC values comparing the assumed mode (uniform-beam modes) and FEM methods with varying damping coefficient. (*Mode 1 is a real mode pair when c = 50 kg/s and 1000 kg/s.) . . . . . . . . . . . 28 Damping ratio and modal frequencies (c = 50 kg/s) compared to undamped modal frequencies . . . . . . . . . . . . . . . . . . . . . 31 Traveling index computed from the complex modes with varying damping coefficient. (*Mode 1 is a real mode pair when c = 50 kg/s or 1000 kg/s.) . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Undamped modal frequencies comparison with and without added mass (Hz) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Table 3.2 Modal frequencies obtained from the model and experiments (Hz) . 50 Table 3.3 Damping ratio comparison for different values of damping coefficient (kg/s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Traveling index comparison for different values of damping coefficient (kg/s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 MAC values comparing the modes from the experimental results and from the model with various damping coefficient (kg/s) . . . . . . . 56 Comparison of edgewise modal frequencies (Hz) of 23-m turbine blade between the static case and the operational case, when Ω is 1.5 rad/s 72 Comparison of edgewise modal frequencies (Hz) of 63-m turbine blade between the static case and and the operational case, when Ω is 1.25 rad/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Comparison of edgewise modal frequencies (Hz) of 100-m turbine blade between the static case and and the operational case, when Ω is 0.78 rad/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Table 2.2 Table 2.3 Table 2.4 Table 3.1 Table 3.4 Table 3.5 Table 4.1 Table 4.2 Table 4.3 x Table 4.4 Table 4.5 Table 4.6 Table 4.7 Table 4.8 Table 4.9 Table 4.10 Table 4.11 Table 4.12 Table 4.13 Table 4.14 Table 4.15 Table 4.16 Comparison of flapwise modal frequencies (Hz) of 63-m turbine blade between the static case and the operational case when the rotation speed is 1.25 rad/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Comparison of flapwise modal frequencies (Hz) of 100-m turbine blade between the static case and the operational case when the rotation speed is 0.78 rad/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Comparison of edgewise modal frequencies (Hz) of the 23-m wind turbine blade at different rotation angles . . . . . . . . . . . . . . . 82 Comparison of flapwise modal frequencies (Hz) of the 23-m wind turbine blade at different rotation angles . . . . . . . . . . . . . . . . . 83 Comparison of edgewise modal frequencies (Hz) of the 63-m wind turbine blade at different rotation angles . . . . . . . . . . . . . . . 83 Comparison of flapwise modal frequencies (Hz) of the 63-m wind turbine blade at different rotation angles . . . . . . . . . . . . . . . . . 83 Comparison of edgewise modal frequencies (Hz) of the 100-m wind turbine blade at different rotation angles . . . . . . . . . . . . . . . 83 Comparison of flapwise modal frequencies (Hz) of the 100-m wind turbine blade at different rotation angles . . . . . . . . . . . . . . . 83 Damping ratio and modal frequencies compared to the undamped modal frequencies for flapwise bending of 23-m turbine blade at rated wind speed u = 16 m/s and Ω = 1.57 rad/s . . . . . . . . . . . . . . 96 Damping ratio and modal frequencies compared to the undamped modal frequencies for edgewise bending of 23-m turbine blade at rated wind speed u = 16 m/s and Ω = 1.57 rad/s . . . . . . . . . . . . . . 97 Damping ratio and modal frequencies compared to the undamped modal frequencies for flapwise bending of 63-m turbine blade at rated wind speed u = 11.4 m/s and Ω = 1.25 rad/s . . . . . . . . . . . . . 98 Damping ratio and modal frequencies compared to the undamped modal frequencies for edgewise bending of 63-m turbine blade at rated wind speed u = 11.4 m/s and Ω = 1.25 rad/s . . . . . . . . . . . . . 99 Damping ratio and modal frequencies compared to undamped modal frequencies for flapwise bending of 100-m turbine blade at rated wind speed u = 11.3 m/s and Ω = 0.78 rad/s . . . . . . . . . . . . . . . . 101 xi Table 4.17 Damping ratio and modal frequencies compared to undamped modal frequencies for edgewise bending of 100-m turbine blade at rated wind speed u = 11.3 m/s and Ω = 0.78 rad/s . . . . . . . . . . . . . . . . 104 Table 4.18 Traveling index for flapwise bending of the three blades at rated wind speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 Table 4.19 “Nonmodal ratio” for flapwise bending of the three blades at rated wind speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 Table 4.20 “Nondiagonality” for flapwise bending of the three blades at rated wind speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Table 5.1 Comparison of the Ω and the first modal frequency ω1 of edgewise and flapwise bending (Unit:rad/s) . . . . . . . . . . . . . . . . . . . 114 Table 5.2 “Nonmodal ratio” for edgewise bending of the three blades at rated wind speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Table 5.3 “Nonmodal ratio” for flapwise bending of the three blades at rated wind speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Table A.1 MAC values comparing the modes from the experimental results and from the model with varying number of samples and varying filter frequency (ff) when c = 10 kg/s . . . . . . . . . . . . . . . . . . . . 141 xii LIST OF FIGURES Figure 2.1 Cantilevered beam with damper at the free end . . . . . . . . . . . 15 Figure 2.2 Real (solid lines) and imaginary (dashed lines) parts of the sampled displacement modal functions generated from the assumed-mode method. In the top left, two real modes are grouped as a “real mode pair”. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Real and imaginary parts plotted against each other in the complex plane (assumed-mode method shown with solid lines, and finiteelement method shown with dots). The top left shows two modes of the real mode pair plotted against each other. . . . . . . . . . . . . 23 Variation of the eigenvalues with increasing number n of assumed modes (uniform-beam modes) . . . . . . . . . . . . . . . . . . . . . 25 Variation of the eigenvalues with increasing number n of assumed modes (uniform-beam modes) for the real mode pair . . . . . . . . 26 Variation of the eigenvalues with increasing number n of assumed modes defined as modified Duncan Polynomials . . . . . . . . . . . 27 Figure 2.7 Beam element with nodal displacements and rotations . . . . . . . . 28 Figure 2.8 Real (solid lines) and imaginary (dashed lines) parts of the displacement modal vectors generated from the finite-element method. In the top left, two real modes are grouped as a ”real mode pair”. . . . . . 29 Mode shapes from the assumed-mode method with varying damping coefficient. (When c = 1 kg/s, mode 1 is complex, so real mode 1 and 2 are real part and imaginary part of mode 1 respectively. ) (a) real parts (b) imaginary parts except for the cases of the real mode pair (c = 50 and 1000 kg/s), in which the second real modes are plotted. 35 Variation of the nonmodal frequencies with increasing damping coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Figure 2.11 Variation of the eigenvalues with increasing damping coefficient . . 36 Figure 2.12 Variation of the damping ratio with increasing damping coefficient . 37 Figure 2.3 Figure 2.4 Figure 2.5 Figure 2.6 Figure 2.9 Figure 2.10 xiii Figure 2.13 Variation of the nonsynchronicity index with increasing damping coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Figure 3.1 The magnetic field of the eddy-current damper . . . . . . . . . . . . 39 Figure 3.2 A schematic diagram of the experimental steup . . . . . . . . . . . 44 Figure 3.3 The experimental setup of the cantilevered beam with an eddy-current damper . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 The comparison for the first four modes with and without mass at undamped case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Comparison of the variation of the traveling index with and without added mass, with increasing damping coefficient . . . . . . . . . . . 48 Comparison of the variation of the damping ratio with and without added mass, with increasing damping coefficient . . . . . . . . . . . 49 The FFT plot of the cantilevered beam experiment with and without eddy-current damper . . . . . . . . . . . . . . . . . . . . . . . . . . 50 Figure 3.8 The acceleration signal shown with different number of samples . . . 51 Figure 3.9 The variation of the eigenvalues with increasing damping coefficient (circles indicate the experimental values) . . . . . . . . . . . . . . . 52 The modes comparison in complex plane when the damping coefficient is 10 kg/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Figure 3.11 The modes comparison when the damping coefficient is 10 kg/s . . . 54 Figure 3.12 The damping ratio comparison between the model and experiments at different damping coefficients . . . . . . . . . . . . . . . . . . . . 55 The traveling index comparison between the model and experiments at different damping coefficients . . . . . . . . . . . . . . . . . . . . 56 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 2000 sample points . . . . . . . . . . . . . . . . . . . . . . . . . 58 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 12500 sample points . . . . . . . . . . . . . . . . . . . . . . . . 59 Figure 3.4 Figure 3.5 Figure 3.6 Figure 3.7 Figure 3.10 Figure 3.13 Figure 3.14 Figure 3.15 xiv Figure 3.16 The modes comparison in complex plane when the damping coefficient is 10 kg/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 The real and imaginary parts of the modes when the damping coefficient is 10 kg/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Figure 4.1 The deflection of an in-plane inextensible beam . . . . . . . . . . . . 64 Figure 4.2 The three-dimensional deflection of an inextensible beam. The coordinate system depicted rotates with the hub. . . . . . . . . . . . . . 65 The modal frequencies for edgewise bending vibration with the variation of rotation speed . . . . . . . . . . . . . . . . . . . . . . . . . 71 Figure 4.4 The first four modes for edgewise bending vibration without rotation 72 Figure 4.5 The comparison between the edgewise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 23-m blade when the rotation speed is 1.5 rad/s . . . . . . . . . . . 73 The comparison between the edgewise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 23-m blade when rotation speed is 15 rad/s . . . . . . . . . . . . . . 74 The comparison between the edgewise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 63-m blade when rotation speed is 1.25 rad/s . . . . . . . . . . . . . 75 The comparison between the edgewise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 100-m blade when rotation speed is 0.78 rad/s . . . . . . . . . . . . 76 The modal frequencies for flapwise bending vibration of 23-m blade with the variation of rotation speed . . . . . . . . . . . . . . . . . . 77 The comparison between the flapwise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 23-m blade when the rotation speed is 1.5 rad/s . . . . . . . . . . . 78 The comparison between the flapwise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 23-m blade when the rotation speed is 15 rad/s . . . . . . . . . . . . 79 Figure 3.17 Figure 4.3 Figure 4.6 Figure 4.7 Figure 4.8 Figure 4.9 Figure 4.10 Figure 4.11 xv Figure 4.12 The comparison between the flapwise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 63-m blade when rotaton speed is 1.25 rad/s . . . . . . . . . . . . . 80 The comparison between the flapwise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 100-m blade when the rotation speed is 0.78 rad/s . . . . . . . . . . 81 The edgewise modes comparision for 23-m wind turbine blades when the blade rotates at horizontal and downward positions . . . . . . . 82 The flapwise modes comparision for 23-m wind turbine blades when the blade rotates at horizontal and downward positions . . . . . . . 84 The edgewise modes comparision for 63-m wind turbine blades when the blade is in horizontal and downward positions . . . . . . . . . . 85 The flapwise modes comparision for 63-m wind turbine blades when the blade is in horizontal and downward positions . . . . . . . . . . 86 The edgewise modes comparision for the 100-m wind turbine blade when the blade are at horizontal and downward positions . . . . . . 87 The flapwise modes comparision for the 100-m wind turbine blade when the blade are at horizontal and downward positions . . . . . . 88 A comparison plot of lift coefficient between experiment data and ONERA dynamic stall model . . . . . . . . . . . . . . . . . . . . . . 89 Figure 4.21 Translating airfoil with lift and drag forces acting . . . . . . . . . . 90 Figure 4.22 The comparision between the modes of the nonmodally damped case and those of the undamped case for 23-m blade . . . . . . . . . . . 97 The first four modes of the flapwise bending of 23-m wind turbine blade with nonmodal damping . . . . . . . . . . . . . . . . . . . . . 98 The first four modes of the edgewise bending of 23-m wind turbine blade with nonmodal damping . . . . . . . . . . . . . . . . . . . . . 99 Figure 4.13 Figure 4.14 Figure 4.15 Figure 4.16 Figure 4.17 Figure 4.18 Figure 4.19 Figure 4.20 Figure 4.23 Figure 4.24 Figure 4.25 The first four modes of the flapwise bending of 63-m wind turbine blade with nonmodal damping . . . . . . . . . . . . . . . . . . . . . 100 Figure 4.26 The first four modes of the flapwise bending of 63-m wind turbine blade with nonmodal damping . . . . . . . . . . . . . . . . . . . . . 101 xvi Figure 4.27 The first four modes of the flapwise bending of 100-m wind turbine blade with nonmodal damping . . . . . . . . . . . . . . . . . . . . . 102 Figure 4.28 The first four modes of the edgewise bending of 100-m wind turbine blade with nonmodal damping . . . . . . . . . . . . . . . . . . . . . 103 Figure 4.29 The variation of flapwise bending damping ratio with varying wind speed for 23-m blade . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Figure 4.30 The variation of flapwise bending damping ratio with varying wind speed for 100-m blade . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Figure 5.1 The example of three mass system . . . . . . . . . . . . . . . . . . . 114 Figure 5.2 The comparison between the condition without perturbation and with perturbation for three mass system . . . . . . . . . . . . . . . . . . 116 Figure 5.3 The comparison between the conditions without perturbation and with perturbation for 23-m blade (a) when t = 0 (b) when t = quarter period. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 Figure 5.4 The comparison between the conditions without perturbation and with perturbation for 63-m blade (a) when t = 0 (b) when t = quarter period. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 Figure 5.5 The comparison between the conditions without perturbation and with perturbation for 100-m blade (a) when t = 0 (b) when t = quarter period. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Figure A.1 The comparison of damping ratio, traveling index, and modal frequencies based on N = 12500 sample points . . . . . . . . . . . . . . 131 Figure A.2 The comparison of the real and imaginary parts of the mode shapes based on N = 12500 sample points when c = 10 kg/s . . . . . . . . 132 Figure A.3 The comparison of damping ratio, traveling index, and modal frequencies based on N = 6000 sample points . . . . . . . . . . . . . . 133 Figure A.4 The comparison of the real and imaginary parts of the mode shapes based on N = 6000 sample points when c = 10 kg/s . . . . . . . . . 133 Figure A.5 The comparison of damping ratio, traveling index, and modal frequencies based on N = 3000 sample points . . . . . . . . . . . . . . 134 xvii Figure A.6 The comparison of the real and imaginary parts of the mode shapes based on N = 3000 sample points when c = 10 kg/s . . . . . . . . . 134 Figure A.7 The comparison of damping ratio, traveling index, and modal frequencies based on N = 2000 sample points . . . . . . . . . . . . . . 135 Figure A.8 The comparison of the real and imaginary parts of the mode shapes based on N = 2000 sample points when c = 10 kg/s . . . . . . . . . 135 Figure A.9 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 2000 sample points . . . . . . . . . . . . . . . . . . . . . . . . . 136 Figure A.10 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 3000 sample points . . . . . . . . . . . . . . . . . . . . . . . . . 136 Figure A.11 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 6000 sample points . . . . . . . . . . . . . . . . . . . . . . . . . 137 Figure A.12 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 12500 sample points . . . . . . . . . . . . . . . . . . . . . . . . 137 Figure A.13 The comparison of the real and imaginary parts of the mode shapes based on N = 2000 sample points with filter frequency at 5 Hz when c = 10 kg/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 Figure A.14 The comparison of the real and imaginary parts of the mode shapes based on N = 2000 sample points with filter frequency at 10 Hz when c = 10 kg/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 Figure A.15 The comparison of the real and imaginary parts of the mode shapes based on N = 3000 sample points with filter frequency at 5 Hz when c = 10 kg/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 Figure A.16 The comparison of the real and imaginary parts of the mode shapes based on N = 3000 sample points with filter frequency at 10 Hz when c = 10 kg/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 Figure A.17 The comparison of the real and imaginary parts of the mode shapes based on N = 12500 sample points with filter frequency at 5 Hz when c = 10 kg/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 Figure A.18 The comparison of the real and imaginary parts of the mode shapes based on N = 12500 sample points with filter frequency at 10 Hz when c = 10 kg/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 xviii Chapter 1 Introduction 1.1 Objective and Motivation Self-adjoint, undamped distributed parameter systems have real normal modes and natural frequencies. We can use separation of variables to generate a differential eigenvalue problem. In some simple situations, we can then solve the equation of motion governing the free vibration without damping. Typically whether the geometry is simple or complicated, the eigenvalue problem often needs a numerical solution. For a damped distributed parameter system, the normal modes are the same as the undamped modes if the system is modally damped in the sense of Caughey [1]. When the system is nonmodally damped, the eigensolutions can have complex values and therefore involve mathematical difficulties in solving the complex differential eigenvalue problem [2–8]. Complex modes have some properties that are different from the real modes. When modes become complex, modal motions may be nonsynchronous [9]. In industry, when rotating parts exist in the system, then the modes could be complex and real normal modes may not be enough as a tool to study the vibration system. The objective of this work are listed as below. Our goal is to study vibration problems with the following two features: nonmodal damping and parametric excitation stiffness. In these two cases, the damping matrix or the stiffness would be non-diagonalizable with an undamped modal coordinate transformation. 1 First, we aim to develop an approach for studying the complex modes and the associated modal properties of a nonmodally damped continuous system. We expect to develop a straightforward, noniterative approach for solving the complex eigenvalue problem. Then this approach will be applied to a cantilevered beam example to examine the modal properties such as modal damping and modal nonsynchronicity. Second, we intend to conduct an experiment of a nonmodally damped beam to verify the complex modal properties which we obtain from the theoretical model based on the approach. With an adjustable damping coefficient, we aim to obtain the modal properties at different damping levels. We also will test the functionality of a pair of complex modal identification methods, state-variable modal decomposition and orthogonal modal decomposition. Third, in the operation of a wind turbine, the blade endures a variety of loads, such as gravitational forces, centrifugal forces and aerodynamic forces. The objective is to study the modes of a rotating wind-turbine blade under these loading conditions. The aerodynamic forces could cause nonmodal damping in the continuous beam model. Therefore, we will apply the approach that we developed to study the nonmodal damping in the vibration of a wind-turbine blade. Fourth, rotation can bring parametric excitation terms in the equation of motion of a wind-turbine blade, and make the modal analysis more difficult. We aim to develop a perturbation analysis approach to study the parametric excitation effect. 2 1.2 1.2.1 Research Background Nonmodally Damped Distributed Systems In distributed parameter systems, like strings, beams, shafts and rods, the mass, stiffness, and a certain class of damping operators are diagonalizable with the undamped modal coordinate transformation. The system may not be diagonalizable when the damping is nonmodal in the sense of Caughey [1], or the stiffness includes the parametric excitation. The partial differential equation of motion governing the free vibration of an undamped or proportionally damped beam can be solved by using the standard separation of variables method, which leads to a differential eigenvalue problem that often requires a numerical solution. When non-modal damping is introduced to the structure, the natural frequencies and modes of vibration become complex and the usual separation of variables solution becomes much more difficult. The general formulation for complex eigenvalue problems is still under development. Singh and others [10–12] developed the eigenproblem formulation and solution techniques which allow the determination of the complex eigenvalues and complex normal modes of transversely vibrating beams having a wide variety of non-proportional viscous damping configurations. An iterative Newton’s method was applied to compute the eigenvalues and further the modal parameters. Hull [13, 14] developed a closed-form series solution for the axial wave equation with a viscously damped boundary condition at the free end. Oliveto et al. [15–17] analyzed a simply supported beam with two rotational viscous dampers attached at its ends and developed the complex mode superposition method. A numerical procedure was also applied to obtain the complex frequencies and modes of vibration. Krenk [18–20] used a state-variable formulation to express the partial differential equation (PDE), which 3 is second-order in time derivatives, as a system of two PDEs, that are first-order in time derivatives. He then used separation of variables to establish an eigenvalue problem system, which was solved iteratively. Sirota and Halevi [21] presented the solution of the free response of damped boundary structures to general initial conditions in the form of an infinite sum of products of spatial and time functions by Laplace domain approach. Some other research studying complex modes can be found in [22–24]. In this work we describe a non-iterative method based on assumed modes, and apply into the complex modal analysis of an end-damped beam. 1.2.2 Modal Identification Methods of Experimental Modal Analysis Output-only modal analysis methods only use output signals such as acceleration, velocity and displacement to extract modal information, which is an advantage in some applications which the input signal is hard to capture. The output-only modal analysis can be time-based or frequency based. A few time-based methods include the eigensystem realization algorithm [25], Ibrahim time domain method [26], independent component analysis [27, 28], least squares complex exponential method [29], Prony’s method [30, 31] and stochastic subspace identification methods [32]. Some frequency-based approaches include orthogonal polynomial methods [33], complex mode indicator function [34] and frequency domain decomposition [35]. As variants of proper orthogonal decomposition methods, the smooth orthogonal decomposition (SOD) [36–38] and the state-variable modal decomposition (SVMD) [39–41] are new time-domain output-only methods developed in recent years and have shown good results for 4 modal analysis of free response simulations. Compared to the Ibrahim time domain method and Prony’s method, SOD or SVMD do not involve the possibility of oversized state matrices and their spurious modes. Compared to frequency domain decomposition, these methods do not involve spectral density functions. In the work of Farooq [41, 42], the state-variable modal decomposition method was applied to the modal parameter identification of structural systems. When applied to a beam experiment, the decomposition method showed good results compared to the theoretical results of an Euler beam. Quantitative assessment was conducted to investigate the effects of samples number and time record length on the accuracy of modal identification. A sensitivity analysis was also conducted to verify the limitations on the number of sensors in the real experiments regarding to the number of active modes. Experimental modal analysis on complex modes is an important topic for the analysis of the nonmodal damping, common in complicated engineering problems, and also wave motions. Ibrahim presented a technique for computing a set of normal modes from a set of measured complex modes [43]. Ewins proposed an approach for the modal identification of lightly damped structures [44]. Feeny presented a complex generalization of proper orthogonal decomposition [9, 45], and used it to decompose wave motion and to extract dominant complex modes. 1.2.3 Eddy Current Damping There are many ways to induce damping in mechanical systems. Eddy currents can provide a contact-free, nearly linear damping source. Eddy currents are generated when a nonmagnetic conductive metal moves in a magnetic field [46–48]. These eddy currents induce their own magnetic field with opposite polarity of the applied field, which causes a resistive 5 force. Eddy currents have been used for many interesting dynamics applications [47]. Eddy currents can be applied to brake system by using the rotational movement of a conductive medium between two oppositely poled magnets to induce an electro-motive force in the material [49–56]. Eddy current damping can be applied to rotational systems to suppress the lateral vibrations of rotor shafts [57–59]. Some other applications of magnetic damping can be seen in [60–62]. Inman et al. [46, 63] built an eddy current damper with permanent magnets and a conducting sheet, and then applied it to the vibration suppression of a beam. The theoretical model for the eddy current damper was derived using the electromagnetic theory, which can help to predict the damping characteristics and therefore the dynamic behavior of the structure. Two geometries of electromagnetic fields were investigated for the eddy-current damping, including the cases when the direction of motion of the conducting sheet is perpendicular and parallel to the magnet’s face. In this work, eddy-current damping will be used to provide linear, but nonmodal damping in an experiment. 1.2.4 Wind-Turbine Blade Modeling Research on rotating cantilevered beams has been popular for years. Such engineering problems include turbine blades and helicopter blades. Here, we summarize just some of the papers which have investigated the vibration of rotating beams. Yoo et al. [64–67] analyzed the coupled bending-axial vibration and uncoupled bending vibration using the RayleighRitz method. Coupled bending-bending-torsion vibration analysis was presented by Anghel et al. [68]. With the development of computer and software, mode shapes, and especially complex mode shapes, were investigated in an easier way [69–72] since the large load of cal6 culation is easier to handle. Cooley and Parker [73,74] generalized previous works and investigated coupled bending and torsional motion. Extended operators and the Galerkin method were applied to discretize the equations. Other related research can be found in [75, 76]. There has been a variety of work on modeling of rotating turbine blades. Hodges and Dowell [77] derived the nonlinear partial differential equations of bend-bend-torsion motion for a twisted helicopter rotor blade. Equations of motion for a rotating wind turbine blade were developed by Wendell [78, 79]. Then, in application to horizontal-axis wind-turbine blades, Kallosoe added gravity and pitch action into the equations of motion [80]. A nonlinear model was developed by Larsen and Nielsen [81, 82] for the coupled edgewise and flapwise bending vibrations of pre-twisted wind-turbine blades with a periodically moving support and applied a reduced form of the model to carry out stability analysis. Turhan and Bulut [83,84] investigated the in-plane bending vibration of a beam rotating with a periodically fluctuating speed. Nonlinear equations that model the flexural potential in non-uniform beam have been developed by Caruntu [85]. Ishida et al. [86] investigated the fundamental vibration characteristic of an elastic blade of the wind turbine. The nonlinear vibration analysis of the superharmonic resonance was performed. Ramakrishnan and Feeny [87] used the method of multiple scales to analyze resonances of the forced nonlinear Mathieu equation, as an approximate representative of blade motion. With the recent increase in the size of wind turbines, the stability and vibration of wind turbine blades has become more important. Sandia National Laboratories have some work on developing a large blade for horizontal-axis wind turbines. A 100-m blade for a 13.2 MW horizontal axis wind turbine was presented [88], and therefore new challenges come with the increase of the size, such as strength, fatigue, deflection and buckling. The research in this area includes the aerodynamic performance of the blades, resonance analysis, and modal 7 analysis. Bir and Oyague [89] showed the modal chracteristics of the 23-m GRC (Gearbox Reliability Collaborative) turbine blade. The modes of natural vibrations of a particular 63-m wind turbine blade were partly investigated in the work of Kallesoe [80, 90]. In time-periodic systems, such as wind turbine blades under steady operating speeds, the time-periodic stiffness matrix presents challenges to modal analysis. In experimental modal analysis, there are a few methods, such as Floquet modal analysis, for analyzing the time varying modes based on data measurements. Allen [91–96] proposed a system identification routine for estimating the modal parameters of a linear time-periodic system including the periodic mode shapes, and applied this to simulated measurements from a model of a rotating wind turbine. However, there is no mature theoretical method for analyzing the effect of the parametric excitation on the modal response of the stationary system. 1.2.5 Aerodynamics of Wind Turbine Blades Wind turbines have been used for centuries. However, research on the aerodynamics on wind turbines has become a popular topic within the last 30 years. Horizontal-axis wind turbines frequently experience high aerodynamic loads during normal service, and these aerodynamic loads bring excessively high stress on the turbine blades, which will finally reduce the service life of the wind turbine [97]. A lot of sources such as the stochastic nature of the wind, flexibility of the wind turbine structure, and the unsteady, three-dimensional character of the flow, combine to make the accurate estimation of the loading impossible [98]. Wind-turbine-blade aerodynamic phenomena can be broadly categorized by machine operating states into two primary cases. Rotational augmentation determines the blade aerodynamic response at zero and low rotor yaw angles. Dynamic stall dominates blade aerodynamics at moderate to high yaw angles [99]. Rotational augmentation includes research such 8 as rotationally induced cross flows and centrifugal forces in the presence of flow separation. Dynamic stall occurs when the effective angle of attack of the airfoil is above its normal static stall angle and accompanied with much larger phase variation in the unsteady airloads due to significant hysteresis in the flow [100]. A simple model for aeroelastic analsyis of a rotor was created by Eggleston and Stoddard [101], in which steady, linear aerodynamic assumptions were employed. Oregon State University developed a computer code FLAP and the analysis was constrained to allow only flapping motions for a cantilevered blade. FLAP uses a linearized blade element method (BEM) aerodynamic model with a frozen wake and can predict the blade load in steady or turbulent winds [98]. To reduce the undertainties in the calculation of aerodynamic loads, the computational fluid dynamics (CFD) methods are applied to the rotor and wake flow [102]. The flow field of wind turbine blades is complex due to a variety of factors, which is shown in previous research [103,104]. Eggers and Digumarthi formulated a model for characterizing the centrifugal and Coriolis effects on a deeply stalled wake residing on the blade suction surface [105]. The rotational and inertial effects were introduced to wind-turbine blade flows by the empirical stall delay model of Tangler and Selig [106]. Hansen [107] presented two experimental methods for estimating the modal damping of a wind turbine during operation. The first method is based on the assumption that a turbine mode can be resonated by a harmonic force at its natural frequency, and thus the decaying response after the end of excitation can be used to estimate the damping. The second method is operational modal analyisis based on stochastic subspace identification [32]. 9 1.2.6 This Work In this work we aim to analyze some continuous systems which cannot be decoupled using the undamped modal coordinate transformation. We focus first on nonmodal damping, where we present a method of analysis and apply it to a structure with concentrated damping elements and then to wind-turbine blades. Wind-turbine blades then prompt us to examine multi-degree-of-freedom and continuous systems with parametric excitation. We first describe an assumed-modes approach for the modal analysis of nonmodally damped continuous systems. Real assumed modes are applied to discretize the system. The discretized eigenvalue problem is then cast in state-variable form to obtain the natural frequencies and state-variable modal vectors based on the assumed modal coordinates. Displacement parts of real and complex state-variable modal vectors are recombined with the assumed modes to approximate characteristic shapes of the original system. The approach is applied to an end-damped cantilevered beam, and characteristics of the mode shapes and modal damping are examined for various values of the damping coefficient. As a comparision, the finite-element method is also applied to get the mass, stiffness, and damping matrices and which are then used to solve the eigenproblem in state-variable form. A cantilevered beam setup with an adjustable eddy-current damper at the end is built and experiments are conducted to verify the modal behavior with the variation of the damping coefficient. State-variable modal decomposition (SVMD) is applied to extract the complex modes from the measured response. The finite-element method is applied to build a theoretical model. The trends of modal properties such as damping ratio and modal nonsychronicity with varying damping coefficient are studied experimentally and compared with the prediction from the model. As verification, complex orthogonal decomposition (COD) is also 10 applied to obtain the modes and modal properties from the free response and compare with the results from SVMD. This is the first work in which SVMD is applied to an experiment with decidedly complex modes. We applied Hamilton’s principle to derive the in-plane and out-of-plane equations of motion of a horizontal-axis wind-turbine blade. The equations of motion are discretized by the assumed-mode approach with uniform cantilevered beam modes. We examine the modal properties with varying rotation speeds and compare the modal properties with and without rotation effects. Nonmodal aerodynamic damping is derived from a quasi-steady lift/drag aerodynamic model. The complex modes are analyzed while neglecting the effects of rotation to single out the effect that aerodynamic damping may have on the modes. Modal properties including damping ratio and modal nonsynchronicity in this nonmodally damped system are investigated. Data from 23-m, 63-m, and 100-m wind turbine blades are applied to the theoretical model. We compare the difference of the modal properties with the variance of the size of blades. The turbine-blade PDEs derived from Hamilton’s principle, under the assumption of steady wind and steady rotation, have parametric excitation. We develop a simple perturbation approach to analyze the perturbation effect of the parametric excitation on the unperturbed system. Since wind turbines are defined to operate well below the first modal frequency, in this initial work, we look at the nonresonant condition. The perturbation of the parametric excitation terms in the equation of motion of a wind turbine blade is analyzed, to investigate the perturbation effect of the parametric excitation terms on the modal response of stationary system. We also take a three mass system as an example to illustrate the modal parametric effect. 11 1.3 Contributions The contributions of this work are listed below ❼ We present a noniterative analytical and numerical approach to describe complex modes in nonmodally damped continuous systems. Most previous studies require iterative solutions or difficult numerical solutions. We sort out the complex and real modes in the original coordinates. The convergence with two choices of assumed-mode trial functions is studied. ❼ We study the behavior of an end-damped Euler-Bernoulli beam for various levels of damping, and show its modal damping and nonsynchronicity. We also reveal the behavior of the fundamental undamped mode as it becomes a real mode pair with increased damping. ❼ The behavior of the beam is examined in an experiment. Modes are extracted by SVMD in its first application to a system with complex modes. ❼ The modal analysis of the damped wind-turbine blade provides insight to the value of using static lab studies to represent rotating blades operating in the wind. ❼ The perturbation study of the parametrically excited systems gives some understanding of the perturbation of modal behavior in non-resonant conditions. 12 Chapter 2 Complex Modal Analysis of a Non-Modally Damped Continuous Beam 2.1 Chapter Introduction This chapter outlines a method in which real assumed-mode trial functions are applied to discretize the system. The discrete eigenvalue problem is then cast in state-variable form to obtain the natural frequencies and state-variable modal vectors. Displacement parts of real and complex state-variable modal vectors are recombined with the assumed-mode trial functions to approximate characteristic shapes of the original system. The approach is applied to study the behavior of the end-damped cantilevered beam. As a comparision, the finite-element method is also applied to get the mass, stiffness, and damping matrices, which are then used to solve the eigenproblem in state-variable form. Characteristics of the mode shapes and modal damping are examined for various values of the damping coefficient. The nonsynchronicity of the complex modes is quantified by a traveling index. The convergence with the variation of the number of assumed modes of different assumed-mode trial functions is examined. The real modal pair is defined and 13 studied when the damping coefficient is large. The optimal damping coefficient for modal properties such as damping ratio and traveling index are presented. 2.2 Problem Formulation In this section we outline a scheme for modal analysis of a nonmodally damped continuous system, and apply it to an end-damped Euler-Bernoulli beam. This approximately represents elastic structural components in bending, with concentrated dampers. The scheme involves two stages of coordinate transformations, including an assumed mode expansion and discretesystem diagonalization. This is similar to the approach used in [80, 108], except here the discrete system is examined in state-variable form. 2.2.1 Assumed-Mode Formulation of a Damped Beam The example of a cantilevered beam with a damper at the free end is shown in Figure 2.1. The lumped damper leads to non-proportional damping. Based on the Euler-Bernoulli beam theory, the equation of motion of the beam is obtained as ρ¨ u + cuδ(x ˙ − l) + ∂2 ∂x2 EI(x) ∂2 ∂x2 u=0 (2.1) where u(x, t) is the transverse displacement, ρ is the mass per unit length, δ(x − l) is the Dirac delta function, in m−1 , c is the damping coefficient in kg/s, E is Young’s modulus, I is the cross-sectional area of moment of inertia, and the dots indicate partial derivatives with respect to time, t. Since the damper is lumped directly into the PDE, the boundary conditions are 14 u(0, t) = 0 ∂ u(0, t) = 0 ∂x ∂2 u(l, t) = 0 ∂x2 ∂3 u(l, t) = 0. ∂x3 Figure 2.1 Cantilevered beam with damper at the free end An assumed modal expansion, which is based on the Galerkin projection [109], is used to reduce the order of the model. The transverse displacement is expanded as u(x, t) ≈ N i=1 ui (x)qi (t), where ui (x) are the assumed mode trial functions, which satisfy the geo- metric boundary conditions, and qi (t) are the generalized modal coordinates of the chosen trial function. We substitute this expansion into equation (2.1), multiply by uj (x), and integrate over the length of the beam to obtain a set of second-order ordinary differential equations, for j = 1, · · · , N . The reduced-order model then has the form M¨ q + Cq˙ + Kq = 0 15 (2.2) where q is a vector of N assumed modal coordinates, and M, K, and C are the mass, stiffness, and damping matrices. The elements of N × N matrices M, K, and C are L mij = ρ(x)ui (x)uj (x)dx (2.3) EIui (x)uj (x)dx (2.4) cui (x)uj (x)δ(x − l)dx (2.5) 0 L kij = 0 L cij = 0 for i, j = 1, · · · , N , where the primes indicate derivatives with respect to space, x. Since the assumed modes are not the true normal modes, M, C and K are generally non-diagonal. Furthermore, C is generally non-modal (non-Caughey [1]). Also, since equation (2.2) represents a self-adjoint system, matrices M, C and K are symmetric. In this case, M, C and K would also result if an assumed-mode (trial function) expansion for u(x, t) were applied to the energy expressions and virtual work in a direct variational derivation. Equation (2.2) can be put in state-variable form by letting yT = [q˙ T qT ], and introducing the equation Mq˙ − Mq˙ = 0 [110]. As such, Ay˙ + By = 0 (2.6)      0 M −M 0  where matrices A =  ,B=  are 2N × 2N . Seeking a solution of the M C 0 K form y(t) = φeαt leads to a general eigenvalue problem, αAφ + Bφ = 0 16 (2.7) Furthermore, the modal solution can be written as y = [q˙ T qT ]T = φeαt . The eigenvectors have the form φT = [vT wT ] = [αwT wT ], where w is the displacement partition of the modal vector [111]. In matrix form, the eigenvalue problem is AΦΛ + BΦ = 0, where Λ is a diagonal matrix of eigenvalues α, and Φ = [φ1 , φ2 , . . . , φ2N ] is the state-variable modal matrix. The state-variable modal matrix can be written in terms of its velocity and displacement partitions as      V   WΛ  Φ = [φ1 , φ2 , . . . , φ2N ] =  = , W W (2.8) where matrix W = [w1 , w2 , . . . , w2N ] is the displacement configuration modal matrix. Thus the multimodal solution is      q˙   WΛ    = Φz =   z, q W (2.9) where z is a vector of state-variable modal coordinates. Then a general response is q = Wz. For the complex-modal oscillations, wj and αj (t), and hence zj , come in complex conjugate pairs. The real form of the associated modal oscillation in q can thus be expressed as qj = (wj zj + wj∗ z∗j )/2, (2.10) where the symbol * denotes complex conjugation. So q is the vector of “assumed-mode” trial function coordinates, i.e. coordinates that modulate the trial functions, and z is a vector of 2N modal state variables, and is a linear combination of vectors q and q. ˙ Using y = Φz in equation (4.31) leads to equations that are decoupled in the modal state variables z. We can organize the state-variable modes to separate the real modes, wjr , j = 1, . . . , 2m, 17 from the complex modes, wk , k = 1, . . . , 2n, and put them in a modal matrix as W = [w1r , . . . , w2mr , w1 , w1∗ , w2 , w2∗ , . . . , wn , wn∗ ] = [Wr , Wc ]. (2.11) We then write the general response as q = Wz, where the modal state variables are now ordered as z = [z1r , . . . , z2mr , z1 , z1∗ , . . . , zn , zn∗ ]T = [zTr , zTc ]T , where the vector zr contains the real-valued state variables and the vector zc contains the complex-valued state variables. The displacement can thus be expressed as q = qr + qc = Wr zr + Wc zc = Wr zr + Hh + H∗ h∗ , (2.12) where the matrix Wc has been organized such that half of its modal vectors are listed in H, and the other half, the complex conjugates, are in H∗ , and likewise the corresponding halves of the complex state variables in vector zc are contained in h and h∗ . Thus, the real-valued displacement modal vectors in Wr characterize the part of the motion in the original M, C, and K system that has real eigenvalues and eigenvectors in the associated state-variable system (4.31), which are the over-damped part of the response, and the complex-valued modal vectors in H (or equivalently H∗ ) characterize the underdamped part of the M, C, and K system. If the damping were proportional (or more generally “modal” in the sense of Caughey [1], i.e. diagonalizable with the undamped modal coordinate transformation), then the real modes would come in pairs that are either identical or proportional, and the displacement partitions, Wc , of the complex modes could also be expressed as real (numerical software may produce complex eigenvectors with linearly dependent real and imaginary parts), although they would 18 correspond to the complex eigenvalues of the state-variable system (4.31). In the case of real eigenvalues, the pairs of identical (or proportional) real eigenvectors are associated with pairs of distinct eigenvalues (if not critically damped) which provide the two time constants of exponential decay of the composite over-damped displacement mode. When the damping is non-modal in the sense of Caughey, there are no longer pairs of identical (or proportional) real eigenvectors, and therefore no associated pairs of time constants that define a single composite displacement mode. However, if the damping matrix can be considered to be a small perturbation of a modal damping matrix, it is likely that “similar” pairs of real modal vectors can be associated. Now we relate the complex modes of the M, C, and K system to the original continuous system whence they came, for which N u(x, t) ≈ uj (x)qj (t) = U(x)q, (2.13) j where U(x) = [u1 (x), . . . , uN (x)] is a 1 × N matrix of assumed modal functions. Using q(t) = Wz(t), and particularly equation (2.12), with equation (2.13), we have u(x, t) ≈ U(x)q(t) = U(x)[Wr zr (t) + Hh(t) + H∗ h∗ (t)], (2.14) u(x, t) ≈ U(x)Wr zr (t) + U(x)Hh(t) + U(x)H∗ h∗ (t). (2.15) or In this expression, u(x, t) generally can consist of an over-damped contribution (should m = 0 and over-damped modes exist) of modal shapes, contained in Ur (x) = U(x)Wr , which are modulated by exponentially decaying state-variable modal coordinates. Since 19 U(x) is 1 × N and Wr is N × 2m, the matrix of real modal functions Ur (x) is 1 × 2m. More commonly for vibration systems, u(x, t) typically consists of an under-damped contribution (i.e. n = 0) of modal shapes, contained in matrix function Uc (x) = U(x)H, or its complex conjugate, which are modulated by exponentially decaying complex state-variable modal oscillations. Since U(x) is 1 × N and H is N × n, the matrix of complex modal functions Uc (x) is 1 × n. The n complex modal functions have real and imaginary parts, which collaborate to describe a non-synchronous oscillation (e.g. see [9]). In the analysis that follows, we will choose assumed modal functions organized in matrix U(x) to discretize an example system (a cantilevered beam with end dashpot) into the form (2.2), and then do a state-variable modal analysis on the state-variable description (4.31) of the assumed-mode discretized system (2.2) to obtain the displacement partitions of, generally, real and complex state-variable modal vectors, listed in Wr and H. We will then characterize the real and complex mode shapes of the non-modally damped cantilevered beam as Ur (x) = U(x)Wr (2.16) Uc (x) = U(x)H. (2.17) 20 2.3 2.3.1 Example : End-Damped Cantilevered Beam Assumed-Mode Formulation In this specific example, the beam properties are: elastic modulus E = 190 GPa, density ρ = 7500kg/m3 , width w = 0.051 m, thickness t = 0.003 m. These parameters imply a mass per unit length m = 1.1907 kg/m, and cross-sectional area moment of inertia I = 1.17 × 10−10 m4 . The damping coefficient is set to c = 50 kg/s. Both real and imaginary plots of the first four modes are shown in Figure 2.2, of which the solid lines are real parts, and the dashed lines are imaginary parts. A real mode pair is also depicted. It turns out that the real mode pair is associated with mode 1 of the undamped system, as we will see later, and the complex modal pairs are associated with other modes. The terminology used here is that “real mode 1” and “real mode 2” refer to the real mode pair associated with mode 1. “Mode 2”, “mode 3” and “mode 4” refer to the associated complex-conjugate modal pair. Figure 2.3 shows the complex domain plots, i.e. imaginary part versus real part parameterized in x, as solid lines. Thus x = 0 is the location of the clamp and is at the origin of the plots. The real mode pair, having no imaginary parts, is instead plotted as real mode 2 versus real mode 1. If a complex mode (or a real mode pair) were synchronous, the plots would form a straight line, since the real and imaginary parts (or a real mode pair) would be proportional. The deviation from the straight line indicates nonsynchronicity, which will be quantified later. 2.3.2 Convergence Study of Different Assumed Mode Functions The convergence of eigenvalues with varying assumed modes number was also investigated. The assumed mode functions applied in this work are the uniform beam modes, as shown in 21 real mode pair mode2 0.3 Real&Imaginary Parts Real part 0.2 0 −0.2 −0.4 0 0.2 0.4 0.6 0.8 Axial location (m) 1 mode3 0 −0.1 0 Real Imaginary 0.2 0.4 0.6 0.8 Axial location (m) 1 mode4 0.2 Real&ImaginaryParts Real&Imaginary Parts 0.1 −0.2 0.2 0.1 0 −0.1 −0.2 0.2 0 0.2 0.4 0.6 0.8 Axial location (m) 0.1 0 −0.1 −0.2 1 0 0.2 0.4 0.6 0.8 Axial location (m) 1 Figure 2.2 Real (solid lines) and imaginary (dashed lines) parts of the sampled displacement modal functions generated from the assumed-mode method. In the top left, two real modes are grouped as a “real mode pair”. equation (2.18). Φ(x)n = sin βn x − sinh βn x − ( sin βn l + sinh βn l )(cos βn x + cosh βn x) cos βn l + cosh βn l (2.18) Using the uniform-beam modes as assumed modes, as the number of assumed modes increases, the eigenvalues converge rapidly, such that the first three pairs of the complex eigenvalues have converged by n =6, and the real eigenvalues converge as n gets larger than about eight. But when the number is larger than 15, an inaccuracy caused by the hyperbolic function calculations emerges, and the eigenvalues diverge, as shown in Figure 22 real mode pair mode 2 0.3 Imaginary part Real mode 2 0.5 0 −0.5 −1 Imaginary part Imaginary part 0.2 0.1 0 −0.1 0 0.05 0.1 Real part 0.1 0 −0.1 −0.08 −0.06 −0.04 −0.02 Real part mode 4 0.2 −0.25−0.2−0.15−0.1−0.05 Real mode 1 mode 3 −0.2 −0.05 0.2 0.15 0 0.1 0 −0.1 −0.2 −0.1 0 0.1 0.2 Real part 0.3 Figure 2.3 Real and imaginary parts plotted against each other in the complex plane (assumed-mode method shown with solid lines, and finite-element method shown with dots). The top left shows two modes of the real mode pair plotted against each other. 2.4 and Figure 2.5. Therefore, the assumed mode number in this paper is set to be eight. A similar convergence behavior was shown in [64], in which the assumed modes functions were polynomial functions. A few modal analysis papers have applied different kinds of polynomial functions as the assumed-mode trial functions, such as Legendre, Chebyshev, integrated Legendre, orthogonalized Duncan, Bardell’s functions and Kwak’s functions [112–115]. Here we applied Karunamoorthy’s modified Duncan polynomials [116] as a comparison to the uniform beam characteristic modal functions. Duncan polynomials for bending are given as 1 1 1 φ(x)n = (n + 2)(n + 3)xn+1 − n(n + 3)xn+2 + n(n + 1)xn+3 6 3 6 (2.19) These equations satisfy the boundary conditions of a cantilevered beam and are linearly 23 Table 2.1 Eigenvalue comparison from two kinds of assumed-mode functions and finite element method, c = 50 kg/s Mode No. real mode pair mode 2 mode 3 mode 4 mode 5 mode 6 mode 7 mode 8 Uniform Beam Modes 1.6146 + 0i 759.12 + 0i 7.856 + 76.312i 25.309 + 253.11i 44.367 + 537.20i 57.906 + 925.38i 65.505 + 1414.9i 67.567 + 1992.3i 65.295 + 2667.4i Modified Duncan Polynomials 1.6145 + 0i 808.71 + 0i 7.860 + 76.286i 25.492 + 252.862i 45.498 + 536.642i 61.519 + 924.687i 68.250 + 1463.61i 288.414+2560.78i 158.606+3985.17i Finite-element 1.6145 + 0i 798.41 + 0i 7.858 + 76.288i 25.454 + 252.88i 45.285 + 536.68i 60.326 + 925.12i 69.26 + 1412.1i 74.297 + 1994.7i 77.201 + 2672.3i independent. By othogonalizing Duncan trinomials using the Gram-Schmidt process, we could obtain the modified Duncan polynomials. By applying the polynomial assumed-mode functions to our problem, the convergence result is shown in Figure 2.6. The convergence is worse for higher modes compared to the results using uniform beam characteristic modal functions. The eigenvalues obtained from the two kinds of assumed-mode functions are listed in Table 3.2. The number of assumed-mode functions are both set to be eight. It shows there is significant difference for the higher modes. 2.3.3 Finite-Element Method As a second way to solve this problem, and to compare the results obtained from the assumedmode method, the finite-element method (FEM) is used to find the mass and stiffness matrices, and further to form the discretized equations of motion. FEM discretizes the whole domain of the beam into smaller elements, and calculates the displacements at discrete nodes. Each beam element has two nodes, and at each node, the transverse nodal displacements di and rotations θi are considered. So four degrees of freedom exist in each element matrix: 24 Norm of Eigenvalue 500 mode2 mode3 400 mode4 300 200 100 0 6 8 10 12 14 16 n Figure 2.4 Variation of the eigenvalues with increasing number n of assumed modes (uniformbeam modes) FEM node point deflections d1 and d2 , and rotations θ1 and θ2 . The element matrices for an Euler-Bernoulli beam are [117, 118]  22h 54 −13h  156     2 2   22h 4h 13h −3h mh   Me =    420   54 13h 156 −22h     2 2 −13h −3h −22h 4h   (2.20)  6h −12 6h   12     2 −6h 2h2  6h 4h EI    Ke = 3     h −12 −6h 12 −6h     6h 2h2 −6h 4h2 (2.21) where h is the length of the beam element. Once the element mass and stiffness matrices have been computed for each element, 25 Real Eigenvalues 8 6 Real Eigenvlue 1 4 Real Eigenvalue 2 100 2 0 6 8 10 12 14 16 n Figure 2.5 Variation of the eigenvalues with increasing number n of assumed modes (uniformbeam modes) for the real mode pair they can be combined using a method of superposition to obtain global mass and stiffness matrices [117]. The global damping matrix was obtained from the boundary condition, as  0   .. .  C=  0   0  ··· 0 0  . . .. ..  . . .    · · · c 0   ··· 0 0 (2.22) This process leads to the equations of motion of the form M¨ x + Cx˙ + Kx = 0, where x = [d1 , θ1 , . . . , dn , θn ]. The vector of modal displacements, and the natural frequencies of the system, can be determined by solving the eigenvalue problem. Since C is nonproportional and non Caughey, we solve the eigenvalue problem in state-variable form, and keep the displacement partition of the eigenvectors. In this example, 50 elements were used in the calculation, as 50 falls into the convergence range. 26 Norm of Eigenvalue 500 mode2 400 mode3 mode4 300 200 100 0 6 7 8 9 10 11 12 13 n Figure 2.6 Variation of the eigenvalues with increasing number n of assumed modes defined as modified Duncan Polynomials For plotting, the transverse displacements terms d1 , d2 of the modal vectors were extracted from the eigenvectors, dropping the terms of rotations θ1 , θ2 . Both real and imaginary plots of the first four transverse-displacement modes are shown in Figure 2.8, of which the solid lines are real parts, and the dashed lines are imaginary parts. Dotted lines in Figure 2.3 show two modes of a real mode pair plotted against each other and the real parts and imaginary parts of the three complex modes plotted against each other in the complex plane. Mode plots obtained from the assumed mode method and the FEM are compared in Figure 2.3, and they match very well. The eigenvalues of the lower modes are shown in Table 3.2. The results show the natural frequencies generated from the assumed-mode method and finite-element method are highly consistent. The modal assurance criterion (MAC) values comparing the assumed-mode method and finite-element method are shown in Table 2.2. The function of the MAC [119] is to provide a measure of consistency between two estimates of a modal vector. This provides an additional confidence factor in the evaluation of modal vectors from different modal parameter 27 Figure 2.7 Beam element with nodal displacements and rotations Table 2.2 MAC values comparing the assumed mode (uniform-beam modes) and FEM methods with varying damping coefficient. (*Mode 1 is a real mode pair when c = 50 kg/s and 1000 kg/s.) damping coefficient c 0 1 mode 1∗ 0.9997 0.9996 mode 2 mode 3 mode 4 0.9997 0.9996 0.9886 0.9886 0.9709 0.9709 50 0.9998 0.9940 0.9956 0.9840 0.9648 1000 0.9998 0.9434 0.9955 0.9836 0.9659 estimation algorithms. The MAC modal scale factor is defined as M ACcdr = ∗ |2 |ψcr T ψdr ∗ ψ T ψ∗ ψcr T ψcr dr dr (2.23) The values near unity in Table 2.2 show that the modes from the assumed-mode method using undamped uniform beam modes and finite-element method are highly consistent. 2.3.4 Interpretation of Beam Modal Properties The real and imaginary parts of the modes plotted with various damping coefficients c are shown in Figure 2.9. The imaginary part dominates the modes when the damping coefficient is set to be 1 or 1000 kg/s. When either part, alone, dominates, the modal motion is nearly 28 real mode pair mode 2 0.3 Real and imaginary parts Real modes 0.2 0 −0.2 −0.4 0 0.2 0.4 0.6 0.8 Axial location (m) 1 mode 3 0.1 0 −0.1 −0.2 0 0.2 0.4 0.6 0.8 Axial location (m) 0.1 0 −0.1 −0.2 0 0.2 0.4 0.6 0.8 Axial location (m) real part imaginary part mode 4 0.2 Real and imaginary parts Real and imaginary parts 0.2 0.2 0.1 0 −0.1 −0.2 1 1 0 0.2 0.4 0.6 0.8 Axial location (m) 1 Figure 2.8 Real (solid lines) and imaginary (dashed lines) parts of the displacement modal vectors generated from the finite-element method. In the top left, two real modes are grouped as a ”real mode pair”. synchronous. As c approaches zero, the modes approach those of the cantilevered beam. When c is around 50 kg/s, the real and imaginary parts are at comparable scales for the modes shown. When the damping coefficient becomes sufficiently larger, the real parts of modes 2, 3 and 4 become small, making their modal motions become more synchronous. As the damping coefficient becomes very large, modes 2, 3 and 4 (and presumably higher modes) of the fixed-free boundary condition with the lumped damper approach the first, second, and third (and presumably higher) modes of a fixed-pinned boundary condition, as shown in Figure 2.10. This is because the damper is nearly rigid. As c increases, mode 29 1 of the cantilevered beam will have become a real mode pair which approaches a static deformation as c reaches infinity. A similar behavior was seen in the case of a torsional damper approaching a fixed-sliding-roller condition [120, 121]. Figure 2.11 provides the plot of the locus of eigenvalues in the complex plane for the first four modes with changing damping coefficient c. The real parts are near zero when c is small. For each eigenvalue except for mode 1, the real part reaches a minimum (most negative) value, and becomes less negative with further increase in the damping coefficient. Thus, each persistently complex mode’s modal damping decay constant has a local maximum with respect to the damping coefficient c for this dashpot arrangement. As the value of c increases, the damped modal frequency gradually decreases, with the most significant decrease occuring in the range in which the decay rate is largest. The decrease in frequency would be consistant with the damped frequency with modal damping. A similar plot was computed in [12] using an iterative numerical scheme. The results in [12], however, did not capture the real mode. The mode pair which is real for c = 50 kg/s comes from a mode that is a complex conjugate pair for low values of c. The associated eigenvalues are plotted with square symbols in Figure 2.11. Certainly as c approaches zero, the modal response is oscillatory, with complex conjugate eigenvalues. When increasing c hits a critical value, the eigenvalues become real and split apart very quickly with further increasing c. The undamped modal frequency cannot be faithfully compared to that computed from non-modally damped eigenvalues, since the non-modally damped eigensolution does not decouple the M, C, K system in the same way. However, considering the complex modal pair as a “mode”, the associated eigenvalues can be expressed in a form α1,2 = −ζωm ± iωm ζ 2 − 1, from which a “modal” frequency ωm and damping ratio ζ can be computed. These values are shown in Table 2.3 for the case of c = 50 kg/s. Compared to the undamped 30 Table 2.3 Damping ratio and modal frequencies (c = 50 kg/s) compared to undamped modal frequencies Modes No. real mode pair mode 2 mode 3 mode 4 ζ ωm ωn (c = 0) 10.864 35.01 17.14 0.1025 76.69 107.4 0.1002 254.16 300.7 0.0841 538.58 589.21 natural frequencies ωn in the undamped case (c = 0), the “modal” natural frequencies ωm are decreased by the effect of non-modal damping. Figure 2.12 shows the plot of damping ratios with changing damping coefficient c. Mode 1 is underdamped when c is small and becomes overdamped when c comes to the critical point around c = 9 kg/s. For each of modes 2, 3 and 4, the value of damping ratio increases with c, before it comes to a local maximum, and then decreases. A nonsynchronicity index or “traveling index” [9] indicates the degree of nonsynchronicity of the real and imaginary parts of the complex modes. The nonsynchronicity index is defined as the reciprocal of the condition number of the matrix whose two columns are the real and imaginary components of the complex mode. The condition number, and thereby the nonsynchronicity index, simultaneously capture the degree of linear independence, and the relative magnitudes between the two vectors, both features of which conspire to affect nonsychronicity. In the case of computing the nonsynchronicity index of a real-mode pair, both vectors are independently normalized, and so the index only captures the degree of linear independence. The relative strength of participation of the real mode pair depends on initial conditions. If the nonsynchronicity index has a value near zero, this is an indication that the modes are highly synchronous, while a value near unity indicates that the modes are highly nonsynchronous. Nonsynchronous modal motions can be used to describe traveling waves, and hence the term “traveling index” was used in [9], although “nonsynchronicity 31 Table 2.4 Traveling index computed from the complex modes with varying damping coefficient. (*Mode 1 is a real mode pair when c = 50 kg/s or 1000 kg/s.) damping coefficient c 0 1 50 1000 ∗ mode 1 0 0.005 0.510 0.646 mode 2 0 0.034 0.221 0.011 mode 3 0 0.018 0.342 0.017 mode 4 0 0.012 0.383 0.023 index” may be more general. In reference [12], effort was made to quantify nonsynchronicity by computing the Im/Re ratio, which is a ratio between the magnitudes of the imaginary and real parts of a complex mode. The Im/Re ratio, however, is dependent on an arbitrary complex scalar constant that can multiply the eigenvectors. It also does not regard the relative shapes of the real and imaginary parts. Table 2.4 provides the values of the nonsynchronicity index for different values of the damping coefficient c. Figure 2.13 shows the variation of the nonsynchronicity indices with increasing damping coefficient c. The nonsynchronicity index of a given mode starts from zero for the undamped case, then increases to a maximum value if the mode remains complex, and finally decreases with further increase in the damping coefficient. For the case of the mode that became a real-mode pair, the traveling index continued to increase with increasing c. Comparison between Figures 2.11 and 2.13 shows that the maximum nonsychronicity of each mode correlates with the mode’s maximum decay rate. As the value of c reaches values for which the dashpot effectively dissipates energy, the modes behave nonsynchronously for this dashpot arrangement. The large c value in Table 2.4 shows how the synchronicity returns for oscillations with large end damping. 32 2.4 Conclusion We have outlined a scheme for the complex modal analysis of non-modally damped distributed parameter systems, and applied this approach to the analysis of an end-damped cantilevered beam. The method involves the use of assumed modes of the undamped system, which are combined based on the state-variable eigensolution of the discretized system. The finite-element method was also utilized to get the mass, stiffness, and damping matrices and further to build and solve a state-variable eigenproblem. The eigenvalues and modal vectors obtained from the assumed-mode method were consistent with those from the finite-element method, as indicated by nearly unit MAC values. These results show that the distributed parameter assumed-mode based method is an efficient way to solve the complex mode problem when non-modal damping is included. The assumed-mode method involves computations of integrals for the low-order modeling of mass, damping and stiffness matrices, and subsequent computations are noniterative and involve lower-order matrices. The assumed-mode method should be equally applicable to structures with nonuniformities, without significant changes in the approach. We applied this method to study features of the end-damped cantilevered beam as a function of the damping coefficient c. With this damping arrangement, most modes are underdamped regardless of c. Each underdamped mode has optimal values of c for generating modal damping and modal nonsynchronicity. As c gets large, the damper becomes more rigid, and the beam behavior approaches that of a clamped-pinned beam. Thus, the contributions of this work are in the presentation of a noniterative method for approximating real and complex modes in continuous, generally damped systems, its verification through comparisons between finite-element discretizations and two choices of 33 assumed-mode trial functions with convergence studies, and in the study of the end-damped beam and the dependence of its modal damping and modal nonsynchronicity on the damping coefficient. 34 0 0.05 −0.1 −0.2 0 0.5 Axial location (m) mode 3 −0.1 1 0.1 0.05 0.05 Real part 0.1 0 −0.05 −0.1 c=1 c=50 c=1000 −0.05 0 Real part mode 2 0.1 Real part Real mode 1* real mode pair 0.1 0 0.5 Axial location (m) mode 4 1 0 0.5 Axial location (m) 1 0 −0.05 0 0.5 Axial location (m) −0.1 1 (a) real mode pair mode 2 0.4 Imaginary part Real mode 2* 1 0.5 0 −0.5 −1 0 0.5 Axial location (m) 0.2 0 −0.4 1 c=1 c=50 c=1000 −0.2 0 mode 3 0.5 Axial location (m) 1 mode 4 0.3 0.4 Imaginary part Imaginary part 0.2 0.1 0 −0.1 0.2 0 −0.2 −0.2 0 0.5 Axial location (m) −0.4 1 0 0.5 Axial location (m) 1 (b) Figure 2.9 Mode shapes from the assumed-mode method with varying damping coefficient. (When c = 1 kg/s, mode 1 is complex, so real mode 1 and 2 are real part and imaginary part of mode 1 respectively. ) (a) real parts (b) imaginary parts except for the cases of the real mode pair (c = 50 and 1000 kg/s), in which the second real modes are plotted. 35 600 Nonmodal Frequencies 500 mode1 mode2 mode3 mode4 clamped−pinned mode1 clamped−pinned mode2 clamped−pinned mode3 400 300 200 100 0 0 10 20 30 40 50 60 Damping coefficient 70 80 90 100 Figure 2.10 Variation of the nonmodal frequencies with increasing damping coefficient 600 500 Imaginary part 400 300 200 100 0 −50 mode1 mode2 mode3 mode4 −40 −30 −20 −10 0 Real part Figure 2.11 Variation of the eigenvalues with increasing damping coefficient 36 0.25 mode1 mode2 mode3 mode4 Damping ratio 0.2 0.15 0.1 0.05 0 0 20 40 60 Damping coefficient 80 100 Figure 2.12 Variation of the damping ratio with increasing damping coefficient 0.7 0.6 Traveling index 0.5 mode1 mode2 mode3 mode4 0.4 0.3 0.2 0.1 0 0 20 40 60 Damping coefficient 80 100 Figure 2.13 Variation of the nonsynchronicity index with increasing damping coefficient 37 Chapter 3 Experimental Study on Complex Modes of an End Damped Continuous Beam 3.1 Chapter Introduction In this chapter, a cantilevered beam experiment setup is built, with an eddy-current damper applied at the end as nonmodal damping. The state-variable modal decomposition (SVMD) is applied to extract the modes from the free impact response of an end-damped cantilevered beam and to study the modal properties such as mode shapes, modal frequencies, damping ratios and modal nonsynchronicity. The experimental results of the modal properties are compared with the model based on the theories in Chapter 2. State-variable modal coordinates and modal assurance criterion values are used to evaluate the quality of the modal decomposition. Complex orthogonal decomposition is applied in comparison to the modal identification results obtained from SVMD. 38 3.2 Eddy Current Damping Formulation An eddy-current damper was built and applied as the damping in this experiment, chosen based on its non-contact and linear properties. The theoretical model for the eddy current damper is derived using the electromagnetic theory in [46, 63], which can help to predict the damping characteristics and the dynamic behavior of the structure. Two arrangements of electromagnetic fields were investigated for the eddy-current damping, including the cases when the moving direction of the conducting sheet is perpendicular and parallel to the magnet’s face. In our experiment, we placed the conducting sheet parallel to the magnet’s face in order to obtain linear damping when the conducting sheet moves in the magnetic field. Figure 3.1 The magnetic field of the eddy-current damper When the conducting sheet moves in the magnetic field of a permanent magnet, as shown in Figure 3.1, eddy currents are induced, based on the induction effect discovered by Faraday [48]. These eddy currents induce their own magnetic field with opposite polarity of the applied field, which causes a resistive force in the opposite direction of the velocity. Thus the movement of the conducting sheet is resisted by forces induced by eddy currents, which can therefore be applied as a non-contact damping mechanism. 39 If the surface charges are ignored [47,122], the current density J induced in the conducting sheet is defined as J = σ(v × B) (3.1) where σ is the conductivity of the conducting sheet, v is the velocity of the conducting sheet, J is the current density, and B is the magnetic flux density. The eddy current force is defined as J × BdV F= (3.2) V where V is the volume between the magnet and the conducting sheet. Substituting equation (3.1) into equation (3.2), as an approximation, the force induced by the eddy currents can be expressed as F = −σδBz2 Se v (3.3) where δ is the thickness of the conductor plate, Bz2 is the magnitude of the magnetic flux in the z direction, Se is the equivalent contact area between the conducting sheet and the magnet, and v is the velocity of the conducting sheet. This force oppose the velocity, and can be interpreted as a damping force, F = −cv The damping coefficient c can be identified as c = σδBz2 Se (3.4) So the induced repulsive force is proportional to the velocity of the conductive metal, and the damping is linear as long as the magnetic flux, conductivity, and sheet thickness are 40 constant. This can be applied in a cantilevered beam [122] to provide nonmodal damping as described in Chapter 2. 3.3 State-Variable Modal Decomposition As a kind of output-only experimental modal analysis method, the SVMD method [41, 108] can be applied to extract the state-variable modes from free responses in the cantilevered beam experiments. The state-variable modal decomposition method enables frequency, damping, and mode-shape estimations from free-response ensemble data for damped linear systems. The eigenvalue problem is formed from correlations of state-variable ensembles. The SVMD eigenvalue problem can be related to the linear model’s state-variable eigenvalue problem and therefore can be used for determining modal parameter information. In SVMD, a generalized eigenvalue problem is derived from measured state-variable ensemble data. Sensed outputs from the free response of the system are incorporated into ensemble matrices and these ensemble matrices are used to create correlation matrices. Then the eigenvalue problem is generated based on the correlation matrices. The method handles both modal and nonmodal damping cases, regardless of whether the modes are real or complex. First, an ensemble matrix Y is built, which is constructed by velocity and displacement vectors Y = [y(t1 ), y(t2 ), · · · , y(tN )] (3.5) where y(t) = [x˙ 1 (t), · · · x˙ M (t); x1 (t), · · · xM (t)]T . T Then a correlation matrix is built as P = YY N and a second correlation matrix is built 41 T ˙ 1 ), · · · , y(t ˙ N )]. as Q = YW N , where W(t) = [y(t Then the SVMD eigenvalue problem is cast as αPφ = Qφ (3.6) Matrix Q is unsymmetric, which allows for complex eigensolutions. The solution of this eigenvalue problem produces an eigenvalue matrix and an eigenvector matrix which contain modal information. In principle, the real and imaginary parts of the eigenvalues indicate modal decay rates and frequencies. The inverse transpose of the eigenvector matrix contains the mode shapes. That is, if Φ = [φ1 , φ2 , · · · , φ2N ], then Ψ = Φ−T = [ψ1 , · · · , ψ2N ] contains estimates of the state variable modal vectors, ψ j , j = 1, · · · , 2N [41]. Since y = [x˙ T ; xT ]T as in equation 3.5, and the state variable modal vectors ψ represent characteristic shapes of y, the lower half of the ψ j indicate the displacement modal vectors. That is, if ψj = [z Tj , wTj ]T , then wj are the displacement modes. 3.4 3.4.1 Cantilevered Beam Experiment Experimental Setup Description This section describes the modal parameter estimation of a cantilevered beam experiment, under varying end-damped conditions. SVMD was applied to extract the modes from the free response of the system. The experimental setup is described below. A uniform aluminum beam as shown in Figure 3.3 was used in this experiment. 30 PCB model number 352B10 accelerometers were placed on the beam via beeswax, equally spaced from the clamp to the beam tip. The beam material properties were taken to be the 42 standard values: elastic modulus E = 68.9 GPa and density ρ = 2700 kg/m3 . The beam had a width w = 0.025 m and thickness t = 0.003 m. The mass per unit length with sensors was m = 0.29 kg/m (without the sensors mass m = 0.21 kg/m from the parameters, but mass of the sensors was accounted for in calculating the theoretical natural frequencies of the system), and the cross-sectional area moment of inertia was I = 5.625 × 10−11 m4 . As shown in Figure 3.2, a piece of copper plate was glued to the end of the cantilevered beam as the conducting sheet in the eddy-current damper. A permanent magnet was placed below the copper plate. The distance between the magnet and the copper was adjustable, and thus the damping coefficient was variable. The length and width of the copper plate was 2 inches and 1 inch. The composition of the permanent magent is NdFeB52 and the maximum magnetic flux of the magnet was measured to be 5000 Gauss. The conductivity of copper is σ = 5.8 × 107 Ω/m, the thickness of the copper was δ = 1.09 × 10−3 m, and the effective area was Se = 9.78 × 10−4 m2 . The magnitude of magnetic flux in the experiment was measured by a Gauss meter at the point A shown in Figure 3.2. The direction of B was assumed to be normal to the copper plate, such that equation 3.3 applies. In the experiment, the distance between the permanent magnet of the copper plate was adjusted between 1 mm and 10 mm, and the magnetic flux was measured to vary between 2800 Gauss and 4500 Gauss. Based on the eddy-current damping theory described above, the damping coefficient was calculated to vary between 5 and 12.5 kg/s. The beam was excited with an impact hammer, and the resulting multi-modal free response was monitored. Measurement signals from the accelerometers were then sampled using a National Instruments data acquisition system (PXI-1042Q). The data was sampled at a rate of 25 kHz. Acceleration signals were first obtained from the experiments and therefore the ensemble 43 Figure 3.2 A schematic diagram of the experimental steup Figure 3.3 The experimental setup of the cantilevered beam with an eddy-current damper matrix of acceleration was first built. The velocity V and displacement X ensembles were obtained by using the Matlab integration routine “cumtrapz”. All ensembles were then processed to remove the respective means. Data were high-pass filtered at a filter frequency of 2 Hz in order to remove a low frequency “drift” in the integrated signal ensembles. The frequency 2 Hz is about half of the theoretical first modal frequency, which is chosen based on the study in [42]. A second order high-pass filter was applied prior to and after each T YWT numerical integration of the signals. The correlation matrices P = YY N and Q = N were then built. Further mode extraction was done in Matlab by post processing. 44 3.4.2 Analysis of the Model Chapter 2 describes a modal-analysis approach for nonmodally damped continuous systems, in which the continuous system is discretized by using assumed modes or finite elements, and the resulting discretized system is analyzed in state-variable form. The state-variable eigenvalues approximate the modal parameters, and the eigenvectors are recombined with the basis functions of the discretization to describe the system modes. In this section, finite-element analysis (FEA) is used to generate the discrete model. The mass and stiffness matrices of the cantilevered beam were obtained from the FEA and the damper at the boundary was used to generate the damping matrix. Then a statevariable eigenvalue problem was generated and the modal frequencies and modal vectors were determined. The uniform beam model was modified to accomodate the added component in this experiment. Since a copper plate was added to the end of the beam, the mass of the copper was approximated and added to the relevant beam elements. 60 elements were used in the calculation. The damping force was distributed into the related beam elements (the last two) in the area of the copper plate. The mass per unit length parameter was modified based on the approximate mass of the 30 sensors. The undamped modal frequencies with and without the added mass were compared (the mass of sensors have already been approximated into the model in both cases), as shown in Table 3.1. It shows that the modal frequencies are significantly reduced when the mass was added to the cantilevered beam. Therefore, the modes have also been affected by the mass, and the comparison for the first four undamped modes is shown in Figure 3.4. A “traveling index” [9, 72] indicates the degree of nonsynchronicity of the real and imag- 45 Table 3.1 Undamped modal frequencies comparison with and without added mass (Hz) Parameter without mass ω1 5.95 ω2 37.3 ω3 104.4 ω4 204.5 with mass 4.8 32.5 94.6 189.9 inary parts of the complex modes. The comparison of the traveling index, for various values of the damping coefficient, with and without the added mass is shown in Figure 3.5. It shows that without the added mass, the peak value of the traveling index for mode 2 is 0.6, while this value drops to 0.32 with the added mass. With added mass, the peak traveling index decreases, and reaches its maximum at larger values of damping coefficient compared to the curve without added mass. Without added mass, it reaches the peak when the damping coefficient is 6 kg/s, while the peak point moves to 9 kg/s with the added mass. For mode 3 and mode 4, the peak value of the traveling index has also been reduced by the added mass. The peak points also shift towards a larger damping coefficient value with the added mass. The peak points for mode 3 and mode 4 shift from 11 kg/s to 22 kg/s and from 16 kg/s to 43 kg/s, respectively. The comparison of the variation of the damping ratio with and without the added mass is shown in Figure 3.6. It shows that without the added mass, the peak value of the damping ratio for mode 2 is 0.23, while this value drops to 0.12 with the added mass. With added mass, the damping ratio decreases and reaches the peak at larger damping coefficient values compared to the curve without added mass. Without added mass, it reaches the peak when the damping coefficient is 6 kg/s, while the peak point shifts right to 9 kg/s with the added mass. For mode 3 and mode 4, the peak value of the damping ratio has also been reduced by the added mass. The peak points also shift towards a larger damping coefficient value 46 mode 2 0.3 0.3 0.2 0.25 0.1 real part real part mode1 0.35 0.2 0.15 −0.1 0.1 −0.2 0.05 −0.3 0 0 0.2 0.4 Axial location (m) −0.4 0.6 0.3 0.2 0.2 0.1 0.1 real part 0.3 0 −0.1 −0.3 −0.3 0.2 0.4 Axial location (m) −0.4 0.6 0.6 mode 4 −0.1 −0.2 0 0.2 0.4 Axial location (m) 0 −0.2 −0.4 0 without mass with mass mode 3 real part 0 0 0.2 0.4 Axial location (m) 0.6 Figure 3.4 The comparison for the first four modes with and without mass at undamped case with the added mass. The peak point for mode 3 and mode 4 shift from 10 kg/s to 21 kg/s and from 15 kg/s to 43 kg/s, respectively. 3.4.3 Modal Identification Results We first verified the modal frequencies and the damping by the fast Fourier transform (FFT). The FFT plot was generated to show the amplitude across the modal frequencies. As shown in Figure 3.7, we can see that the peaks at the first three modal frequencies (modes 2, 3 and 4) have an obvious reduction, presumably, by the effect of damping. The first mode was not obtained in the experiments because the first modal frequency 4.8 Hz is below the limit of the accelerometers (5∼10 Hz). So we were unable to generate a meaningful first mode in our modal decomposition results. Furthermore, as the damping 47 0.7 mode2 no mass mode3 no mass mode4 no mass mode2 with mass mode3 with mass mode4 with mass 0.6 Traveling index 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 25 30 Damping coefficient 35 40 45 50 Figure 3.5 Comparison of the variation of the traveling index with and without added mass, with increasing damping coefficient coefficient increases, the first mode becomes a real mode pair, which is not a decaying oscillation but an exponential decay. Essentially, it is a zero-frequency mode and cannot be captured by the accelerometers. The first modal frequency could be raised by shortening the beam or increasing the thickness of the beam. However, the damping will become less significant in that case. The first N = 12500 points (0.5 seconds) were kept for data processing to minimize the contribution of the high frequency noise dominating the decayed signals in the later part of the ensembles [41]. For extracting modes 3 and 4, since the high-frequency modes were damping out faster, the data were further pared down to N = 3000 points (0.12 seconds) or N = 2000 points (0.08 seconds). The acceleration of sensor 1 with different time length is shown in Figure 3.8. When number of the samples are 12500 (0.5 seconds), the signal still has a small oscillation in the latter part, which we expect to include mostly the mode 48 0.25 mode2 no mass mode3 no mass mode4 no mass mode2 with mass mode3 with mass mode4 with mass Damping ratio 0.2 0.15 0.1 0.05 0 0 5 10 15 20 25 30 Damping coefficient 35 40 45 50 Figure 3.6 Comparison of the variation of the damping ratio with and without added mass, with increasing damping coefficient information for mode 2. Therefore, it is reasonable to use a longer sampling time like 12500 (0.5 seconds) to identify mode 2. When the number of samples are 3000 or 2000, there are more high frequency effects shown in the signal. So a shorter samples could give a better modal identification result for higher modes. A detailed comparison of results with different data sampling is shown in Appendix A. The theoretical frequency values, the frequency values from the FFT, and the SVMD identified modal frequencies are shown in Table 3.2. These values were obtained when the damping coefficient was 10 kg/s. The theoretical frequencies of modes 2, 3 and 4 are 32.5, 94.6 and 189.9 Hz, and the SVMD identified experimental modal frequencies are 32, 95.2 and 189.8 Hz. The modal frequencies between the finite-element model and the experimental results were consistent, although with a difference between 1% and 3%. 49 12000 No damping When damping coefficient is 10kg/s 96 10000 |FFT(x) | 8000 190 6000 4000 2000 32 0 50 100 150 200 250 300 Frequency [Hz] 350 400 450 500 Figure 3.7 The FFT plot of the cantilevered beam experiment with and without eddy-current damper Table 3.2 Modal frequencies obtained from the model and experiments (Hz) Parameter FEA model FFT SVMD ω2 32.5 32.9 33.6 ω3 94.6 96.0 93.8 ω4 189.9 190.0 187.9 The eigenvalues in comparison to the theoretical results from the model with varying damping coefficients are shown in Figure 3.9. The values were obtained using different data sampling and filter frequencies, and optimal values are presented. The model predicts that the real parts are nearly zero when c is small. When c increases, each real part reaches a minimum (most negative) value, and becomes less negative with further increase in the damping coefficient. The variation of the experimental results coincides with those results from the model. 50 40 acceleration acceleration 40 20 0 −20 −40 0 5000 samples −20 1000 2000 samples 3000 40 acceleration acceleration 0 −40 0 10000 40 20 0 −20 −40 0 20 500 1000 samples 1500 2000 20 0 −20 −40 0 200 400 600 samples 800 1000 Figure 3.8 The acceleration signal shown with different number of samples The experimental modes obtained from SVMD were compared with the theoretical modes from a computational model. Figure 3.10 shows the modes plotted in complex plane, with imaginary part against the real part. Figure 3.11 shows the modes plotted with imaginary part and real part separately. For modes 2 and 3, the modes were obtained using N = 12500, with the filter frequency set as 2 Hz. For mode 4, the mode were obtained using N = 3000, with the filter frequency set as 2 Hz. A detailed study and comparison of the mode shapes with the variation of data sampling and filter frequency can be found in Appendix A. We can see from both Figure 3.10 and Figure 3.11, that the modal identification of the lowest underdamped mode, i.e. mode 2, is good. The modes selected from the experimental result are highly consistent with the model. For mode 3 and mode 4, the experimental modes and model’s modes nearly coincide with each other, with a small deviation at the latter half. The consistency between the experimental results and the model can be quantified by the 51 1400 1200 imaginary part 1000 800 600 400 200 0 −120 mode 2 mode 3 mode 4 experiment −100 −80 −60 real part −40 −20 0 Figure 3.9 The variation of the eigenvalues with increasing damping coefficient (circles indicate the experimental values) values of MAC, which will be shown later. The quality of the modes could also be indicated by the modal coordinates and will be presented in next section. Considering the complex modal pair as a “mode” [72], the associated eigenvalues can be expressed in the form α1,2 = −ζωm ± iωm ζ 2 − 1, from which a “modal” frequency ωm and a “modal damping ratio” ζ can be computed. The values of damping ratios calculated with various damping coefficients are shown in Table 3.3. They are compared with the model and shown in Figure 3.12. As with the modal frequencies and mode shapes, the values were obtained using different data sampling and filter frequencies, and optimal values are presented. The model predicts that with the increase of the damping coefficient, the damping ratio of mode 2 increases to a maximum value when the damping coefficient is around 10 kg/s, and then decreases. For modes 3 and 4, the damping ratio keeps increasing in the range of the experimental damping coefficient. The trend of the experimental results 52 mode 2 0.3 Imaginary part 0.2 0.1 0 −0.1 Experiments Model −0.2 −0.2 −0.1 0 0.1 Real part 0.2 0.3 mode 4 0.3 0.2 0.2 Imaginary part Imaginary part mode 3 0.3 0.1 0 −0.1 −0.2 0.1 0 −0.1 −0.2 −0.2 −0.1 0 0.1 Real part 0.2 0.3 −0.1 −0.05 0 Real part 0.05 0.1 Figure 3.10 The modes comparison in complex plane when the damping coefficient is 10 kg/s matches with the model prediction. The values of traveling indices calculated with various damping coefficients are shown in Table 3.4. They are compared with the model results and shown in Figure 3.13. The model predicts that with the increase of the damping coefficient, the traveling index of mode 2 increases to a maximum value when the damping coefficient is around 10 kg/s, and then decreases. For modes 3 and 4, the traveling index keeps increasing in the range of the experimental damping coefficient. The trend of the presented experimental results obeys the model prediction. The nonsynchronicity is also apparent in the plot of the modes when the damping coefficient is 10 kg/s, since the phases of the real and imaginary parts have a significant difference, as shown in Figure 3.11, and the complex mode plot is not confined to 53 mode 2 Real & Imaginary part 0.4 0.2 0 −0.2 Experiments Model −0.4 0 0.1 0.2 0.3 0.4 location (m) 0.5 0.6 mode 3 mode 4 0.4 Real & Imaginary part Real & Imaginary part 0.4 0.2 0 −0.2 −0.4 0 0.1 0.2 0.3 0.4 location (m) 0.5 0.2 0 −0.2 −0.4 0.6 0 0.1 0.2 0.3 0.4 location (m) 0.5 0.6 Figure 3.11 The modes comparison when the damping coefficient is 10 kg/s a line, as in Figure 3.10. The modal assurance criterion (MAC) values comparing the selected modes from the experiment and the model are shown in Table 3.5. The function of the MAC [119] is to provide a measure of consistency between two estimates of a modal vector. This provides an additional confidence factor in the evaluation of modal vectors from different modal parameter estimation algorithms. A value near unity means the modes compared are highly Table 3.3 Damping ratio comparison for different values of damping coefficient (kg/s) mode number mode 2 mode 3 mode 4 c=5 0.078 0.026 0.007 c=7.5 c=10 c=12.5 0.108 0.123 0.113 0.031 0.05 0.055 0.009 0.012 0.018 54 Table 3.4 Traveling index comparison for different values of damping coefficient (kg/s) mode number mode 2 mode 3 mode 4 c=5 0.197 0.069 0.018 c=7.5 c=10 c=12.5 0.283 0.319 0.271 0.121 0.137 0.159 0.052 0.062 0.069 0.14 mode 2 mode 3 mode 4 experiment 0.12 Damping ratio 0.1 0.08 0.06 0.04 0.02 0 0 10 20 30 Damping coefficient 40 50 Figure 3.12 The damping ratio comparison between the model and experiments at different damping coefficients consistent. The MAC modal scale factor is defined as M ACcdr = ∗ |2 |ψcr T ψdr ∗ ψ T ψ∗ ψcr T ψcr dr dr (3.7) where the ∗ denotes complex conjugation. The values in Table 3.5 show that the modes from the experiments and those from the model are particularly consistent for the first mode. 55 0.35 mode 2 mode 3 mode 4 experiment 0.3 Traveling index 0.25 0.2 0.15 0.1 0.05 0 0 10 20 30 Damping coefficient 40 50 Figure 3.13 The traveling index comparison between the model and experiments at different damping coefficients Table 3.5 MAC values comparing the modes from the experimental results and from the model with various damping coefficient (kg/s) damping coefficient c 5 7.5 10 12.5 mode 2 0.9807 0.9780 0.9985 0.9892 mode 3 0.9556 0.9193 0.9767 0.9857 mode 4 0.9242 0.8497 0.9231 0.8126 3.4.4 Modal Coordinates Modal coordinates can help to distinguish the true modes from the spurious modes, and indicate the quality of the decomposition [41, 42]. The state-variable vector can be written as n y(t) = qi (t)ψ i=1 where qi (t) are the modal-coordinate state variables. 56 (3.8) In ensemble matrix form, Y = ΨQ (3.9) where elements of the ith row of Q are the ith sampled “state variable modal coordinates” (SVMC) or “modal state variables”. Thus the modal state variables are simply estimated by Q = Ψ−1 Y = ΦT Y (3.10) In vector form, y(t) = Ψq(t), and q(t) = ΦT y(t). The modal coordinates obtained from the SVMD are shown in Figure 3.14, for the case of N = 2000. A smooth periodic non-noisy modal coordinate history can indicate a true mode. From the plot, we can see mode 2 and mode 3 have smoothly shaped modal coordinates. For mode 4 and mode 5, the latter half of the modal coordinate history is distorted by noise. We can still consider the mode as a candidate for an actual mode if the first half is “good”. Quantification of “good” is discussed in detail in [42]. The quality of the modal frequencies can help further to determine if the mode is a true mode or a spurious mode. The modal coordinates obtained from the SVMD when the samples are 12500 are shown in Figure A.12. From the plot, we can see mode 2 has a smoothly shaped modal coordinate. For mode 3, 4 and 5, the latter half of the modal coordinate history is distorted by noise. Since higher modes generally dissipate quickly, when the samples are large, the latter part cannot capture the the higher modes. Therefore, it is better to use a small sample to identify the higher modes. The effect of the data record length N is examined in more detail in the appendix. 57 −6 6 −6 x 10 4 2 SVMC 3 SVMC 2 4 2 0 0 −2 −2 −4 0 500 1000 1500 Time Histoy −4 2000 −6 1 x 10 0 500 1000 1500 Time Histoy 2000 500 1000 1500 Time Histoy 2000 −7 x 10 5 x 10 SVMC 5 SVMC 4 0.5 0 0 −0.5 −1 0 500 1000 1500 Time Histoy −5 2000 0 Figure 3.14 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 2000 sample points 3.5 3.5.1 Complex Orthogonal Decomposition Verification Complex Orthogonal Decomposition Theory Here we would like to use another modal decomposition method, complex orthogonal decomposition (COD) [9], to check what we have obtained from SVMD. First, real measured signals need to be expressed in complex form. Suppose the real signal is yj = y(xj , t), j = 1, · · · , M where M is the number of sensors distributed on the structure. Applying the Hilbert transform of yj (t), i.e. yHj (t) = Im(zj (t)), then the complex analytic signal can be written as zj (t) = yj (t) + iyHj (t) [9, 123]. Given the signals zj (t), j = 1, · · · , M , we generate vectors zj = [zj (t1 ) · · · zj (tN )]T , by sampling at times t1 through tN . Then an M × N complex ensemble matrix is assembled as 58 −4 1 −5 x 10 4 2 SVMC 3 SVMC 2 0.5 x 10 0 −0.5 0 −2 −1 0 2000 4000 6000 −4 0 8000 10000 12000 2000 Time Histoy 2 x 10 1 0 SVMC 5 SVMC 4 8000 10000 12000 −5 x 10 −5 −10 0 6000 Time Histoy −5 5 4000 0 −1 2000 4000 6000 −2 0 8000 10000 12000 Time Histoy 2000 4000 6000 8000 10000 12000 Time Histoy Figure 3.15 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 12500 sample points ¯ T , where the Z = [z1 · · · zM ]T . Then we built a complex correlation matrix R = (1/N )ZZ bar indicates complex conjugation. ¯ T is complex Hermitian, it has real eigenvalues and complex eigenvectors. Since R = R ¯ Ti uj = δij The Hermitian nature of R implies that the normalized eigenvectors satisfy u where δij is the Kronecker δ, and uj are dimensionless normalized eigenvectors. The eigenvectors of a complex correlation matrix R are the complex orthogonal wave modes (COMs), and the eigenvalues (COVs) indicate the mean squared amplitudes of COD modal coordinates [9]. A theoretical connection between COMs and linear normal modes has not been established. However, if one mode dominates, we might expect it to be representative. 59 3.5.2 Modal Identification Results The experimental modes obtained from COD were compared with the theoretical modes from the model. Figure 3.16 shows the modes plotted in the complex plane, with the imaginary part against the real part. Figure 3.17 shows the modes plotted with imaginary part and real part separately. The COVs corresponding to modes 2, 3 and 4 were 5.95 × 103 , 4.13 × 103 , and 1.15 × 103 , respectively, indicating that modes 2 and 3 were dominant, but not overwhelmingly. mode 2 0 −0.2 Imaginary part −0.4 −0.1 Experiments Model −0.05 0 0.05 Real part mode 3 0.1 mode 4 0.3 0.3 0.2 0.2 Imaginary part Imaginary part 0.2 0.1 0 −0.1 −0.2 −0.2 −0.1 0 0.1 0.2 Real part 0.1 0 −0.1 −0.2 −0.2 −0.1 0 0.1 0.2 Real part Figure 3.16 The modes comparison in complex plane when the damping coefficient is 10 kg/s We can see from both Figure 3.16 and Figure 3.17, that the modal identification of the lowest underdamped mode, i.e. mode 2, is very good. The modes from both the experimental result and the model are highly consistent. For mode 3 and mode 4, the experimental modes and the model correlate with each other, but have small deviations at the latter half. The modes obtained from COD are supportive of the SVMD results. 60 0.2 0 −0.2 −0.4 Experiments Model 0 0.2 0.4 location (m) mode 3 0.6 mode 4 0.4 Real & Imaginary part Real & Imaginary part Real & Imaginary part mode 2 0.4 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 0.4 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 Figure 3.17 The real and imaginary parts of the modes when the damping coefficient is 10 kg/s 3.6 Conclusion We have conducted an experiment for investigating the complex modal behavior of an enddamped cantilevered beam. Eddy-current induced damping was applied at the end of a cantilevered beam, to generate the nonmodal damping. SVMD was applied to extract the complex modes from the free response of the cantilevered beam and to analyze the modal properties, including modal frequencies, mode shapes, damping ratio and modal nonsychronicity. The modal properties from the experimental results were compared with those from a numerical analysis of the model. The extracted modal parameters were quite sensitive to sampling parameters, the study of which was added to the appendix. The trends of the modal frequencies and mode shapes obtained from the experiments were consistent with those from the model. The variation of the damping ratio and traveling index with varying 61 damping coefficient also followed the predictions from the model. Over the range of damping coefficients studied in the experiments, we observed a maximum damping ratio in the lowest underdamped mode, which correlated with the maximum modal nonsynchronicity. The results of modal assurance criterion values showed that the SVMD modes were consistent with the modes from the model. This is the first experiment in which SVMD was used to extract structural normal modes that were anticipated to be complex. As verification, the COD method also produced similar results as SVMD method. The mode shapes obtained from COD correlated with the modes from the model very well. 62 Chapter 4 Modal Analysis of a Wind Turbine Blade 4.1 Chapter Introduction Modal analysis of wind-turbine blades involve complicated issues, such as rotation and aerodynamic effects. In this chapter, we will look at these effects separately. Chapter 5 will address the modes with dynamic rotation. A horizontal-axis wind-turbine blade can be considered as a rotating cantilevered beam, subjected to gravitational loading and aerodynamic loading. Industry tends to consider in-plane (edgewise) deflection separate from out-of-plane (flapwise) flection, and we will do the same in this work. The in-plane and out-of-plane motions are examined with simple aeroelastic damping effects. The aeroelastic model used is based on a simple quasi-steady blade-element airfoil theory. The complex modes are analyzed while neglecting the rotation effects and parametric excitation terms, and thus the nonmodal damping effect could be isolated. The assumptions in this work include: the beam is inextensible, and therefore, the axial displacement is not considered (See Appendix 2); the rotation speed of the blade was assumed to be constant; a quasi-steady model was applied for lift coefficient formulation; the airfoil type was assumed to be the same along the length of the blade for the nonmodal damping 63 modal analysis. 4.2 Equation of Motion Formulation The equation of motion of the rotating beam is obtained by applying the extended Hamilton’s principle [124]. Referring to Figure 4.1 and Figure 4.2, we sketch the process for an inextensible beam with large deflection (a detailed explanation of the inextensible assumption could be found in Appendix B), rotating in a horizontal axis with flexure in the plane of rotation. The potential energy is formulated to include gravitational loading, nonlinear curvature for large excitation, and the effect of the shortening of projection of the beam. Figure 4.1 The deflection of an in-plane inextensible beam In the structure of a horizontal-axis wind turbine, the turbine blade is attached to a hub which we assume to rotate at a constant rate. In reality, the speed will be slowly varying, considering the aerodynamic forces on the blades, the coupling of blades to the rotor, and 64 the factors brought by generator, gear box, friction and control. Here we start with the case of fixed angular speed Ω. We will first formulate the edgewise bending equation of motion, and then the flapwise bending equation. Figure 4.2 The three-dimensional deflection of an inextensible beam. The coordinate system depicted rotates with the hub. 4.2.1 Edgewise Bending In Figure 4.2, y(x, t) is the edgewise deflection. For the edgewise bending, the kinetic energy, the gravitational potential energy and the bending potential energy are obtained as L T = L1 1 ˙ 2 dx m(x)v · vdx + J(x)(y + φ) 2 2 0 0 L Vg = m(x)g[x(1 − cos φ) + s(x, t) cos φ + y(x, t) sin φ]dx 0 L Vb = 1 EI(x)y 2 (1 − 3y 2 )dx 2 0 65 (4.1) where the velocity is v = (−s˙ − Ωy)er + [Ω(x − s) + y]e ˙ φ , and the foreshortening at a position 2 4 x along the beam is s(x, t) = 0x ( y2 + y8 )dx. The virtual work due to nonconservative forces per unit length, f (x, t), is δW = 0L f (x, t)δydx. By applying extended Hamilton’s principle, 0L (δT − δV + δW )dt = 0, an integral-partial differential equation of motion for edgewise bending vibration is obtained as m(−¨ y + Ω2 y + 2ω s) ˙ + (Jz (x)¨ y ) − (EIy − 3EIy y 2 ) − (3EIy 2 y ) −m(¨ s + 2Ωy˙ + Ω2 (x − s) + g cos(Ωt))(y + y 3 /2) L + m(¨ s + 2Ωy˙ + Ω2 (w − s) + g cos(Ωt))dw(y + y 3 /2) x +f (x, t) = mg sin Ωt (4.2) with boundary conditions y(0, t) = y (0, t) = 0 at x = 0 and Jz (L)¨ y + (EIy − 3EIy y 2 + 3EIy 2 y ) = 0 EIy − 3EIy y 2 = 0 at x = L, where y(x, t) is the transverse beam displacement, w is an integration variable for the x domain, E is the Young’s modulus, I is the area moment of inertia of a cross-section of the beam and Jz is the mass moment of inertia per unit length, both about the neutral axis in the z direction. 66 4.2.2 Flapwise Bending In Figure 4.2, z(x, t) is the flapwise deflection. For flapwise bending, the kinetic energy, the gravitational potential energy and the bending potential energy are obtained as L 1 m(x)v · vdx 0 2 T = L m(x)g[x(1 − cos φ) + s(x, t) cos φ]dx Vg = (4.3) 0 L 1 EI(x)z 2 (1 − 3z 2 )dx 0 2 Vb = where the velocity is v = −se ˙ r +Ω(x−s)eφ +ze ˙ k , and the foreshortening at a position x along 2 4 the beam is s(x, t) = 0x ( z2 + z8 )dx. The virtual work of nonconservative aerodynamic and damping forces is δW = 0L f (x, t)δzdx. The equation of motion for flapwise bending vibration is obtained as −m¨ z + (Jy (x)¨ z ) − (EIz − 3EIz z 2 ) − (3EIz 2 z ) −m(¨ s + Ω2 (x − s) + g cos(Ωt))(z + z 3 /2) L + m(¨ s + Ω2 (w − s) + g cos(Ωt))dw(z + z 3 /2) + f (x, t) = 0 x with boundary conditions z(0, t) = z (0, t) = 0 at x = 0 and Jy (L)¨ z + (EIz − 3EIz z 2 + 3EIz 2 z ) = 0 EIz − 3EIz z 2 = 0 67 (4.4) at x = L, whereI is the area moment of inertia of a cross-section of the beam and Jy is the mass moment of inertia per unit length, both about the neutral axis in the y direction. In the equation of motion obtained, there are geometric nonlinear terms, linear terms with parametric excitation, i.e. sin Ωt, cos Ωt terms. These terms derive from the gravitational potential energy and possibly the aerodynamic force as well. 4.3 Linearization and Modal Reduction To linearize the equation of motion 4.2, we ignore the nonlinear cross terms with the velocity and acceleration via the s˙ and s¨ terms, and drop the aerodynamic force term, to get the equation for the edgewise model as m¨ y + (EIy ) − mΩ2 y + mΩ2 xy − ( L mΩ2 wdw)y x L mgdw cos Ωt = mg sin Ωt +mg cos Ωty − (4.5) x The flapwise bending vibration is governed by the equation of motion as shown in equation (4.4). When linearized, this equation reduces to m¨ z + (EIz ) + mΩ2 xz − ( L mΩ2 wdw)z x L +mg cos Ωtz − mgdw cos Ωt = 0 (4.6) x These equations are parametrically excited, which poses complications. For wind turbines 68 that are not too large, we can ignore this. We address it later. An assumed modal expansion is used to reduce the order of the linearized equation of motion. According to the assumed-modes method, the transverse displacement is expanded as y(x, t) ≈ N i=1 ui (x)qi (t), where ui (x) are the assumed modes derived from the uniform cantilevered beam, and qi (t) are the assumed modal coordinates. This expression is substituted into the equation (4.5). Multiplying by uj (x), and integrating over the length of the turbine blade, we get the second-order ODE as 0 m(x)ui uj dx + qj L +qj 0 0 m(x)gui uj dx cos Ωt − qj 0 EI(x)ui uj dx − qj m(x)Ω2 ui uj dx − qj L L +qj L L L q¨j L Ω2 ( L x 0 0 mwdw)ui uj dx L L mdw)ui uj dx cos Ωt = mg g( 0 mΩ2 ui uj dx x 0 ui dx sin Ωt (4.7) The second order ODE which governs the flapwise bending vibration is obtained in a similar way as L q¨j L +qj 0 0 L m(x)ui uj dx + qj L +qj 0 L m(x)Ω2 ui uj dx − qj L m(x)gui uj dx cos Ωt − qj 0 69 x EI(x)ui uj mwdwui uj dx L g( 0 Ω2 0 L x mdw)ui uj dx cos Ωt = 0 (4.8) 4.3.1 Edgewise Bending Modal Analysis Without considering the aerodynamic force and parametric excitation, the damping matrix is neglected, and then the mass and stiffness matrices are obtained as L mij = L kij = L + 0 0 m(x)ui uj dx 0 L EI(x)ui uj dx − m(x)Ω2 ui uj dx − L 0 Ω2 0 mΩ2 ui uj dx L x mwdwui uj dx (4.9) The dimensional parameters of a 23-m blade were obtained from the National Renewable Energy Laboratory (NREL) technical report of Bir and Oyague [89]. Figure 4.3 shows the plot of the modal frequencies of the first four modes with varying rotation speed. The range of the rotation speed varies from 1 rad/s to 10 rad/s in the figure. However, for most wind turbine blades, the speed is around 1-2 rad/s. The modal frequencies obtained show good agreement with the values published in the report of Bir and Oyague. It is observed that when the rotation speed increases, the modal frequencies increases slightly. This is because the centrifugal inertia forces increases as the angular speed increases. It obeys the prediction shown in reference [64]. Over the 1-2 rad/s of most wind turbines, however, the rotation speed has little effect on the modal frequencies. As such, we consider the edgewise bending mode shapes without rotation, which means the terms with rotation speed in the mass and stiffness matrices are neglected at this moment as 70 45 40 mode1 mode2 mode3 mode4 Modal Frequencies 35 30 25 20 15 10 5 0 1 2 3 4 5 6 Rotation Speed 7 8 9 10 Figure 4.3 The modal frequencies for edgewise bending vibration with the variation of rotation speed L mij = 0 m(x)ui uj dx L kij = 0 EI(x)ui uj dx (4.10) The first four vibration modes are shown in Figure 4.4. To check the effect of rotation on the mode shapes, we then include the terms with rotation in the calculation and compare the mode shapes under the two conditions, the static case which doesn’t include the rotation terms and the operational case which include the rotation terms. Figure 4.5 shows the mode shape variation with and without rotation. The dashed lines represent the mode shape when the wind turbine blade is stationary. The cross marked lines represent the mode shape when the wind turbine blade is rotating. The 71 Modal Displacement −0.1 −0.2 −0.3 −0.4 Modal Displacement edge mode 2 0 20 40 Axial location edge mode 3 0.2 0.1 0 −0.1 −0.2 0 20 40 Axial location 0.2 0 −0.2 −0.4 −0.6 60 Modal Displacement Modal Displacement edge mode 1 0 20 40 Axial location edge mode 4 60 0 20 40 Axial location 60 1 0.5 0 −0.5 60 0 Figure 4.4 The first four modes for edgewise bending vibration without rotation Table 4.1 Comparison of edgewise modal frequencies (Hz) of 23-m turbine blade between the static case and the operational case, when Ω is 1.5 rad/s Modes No. static case mode 1 3.14 mode 2 10.10 mode 3 22.95 mode 4 42.88 operational case 3.16 10.12 22.98 42.91 rotation speed here is 1.57 rad/s. We can see there is no observable difference between the modes shapes under the two conditions. But this will change when the rotation speed increases. Figure 4.6 shows the mode shapes when the rotation speed is set to 15 rad/s. We can see a small difference between the modes without rotation and those with rotation. The edgewise modal frequencies between the static case and the operational case are compared. The comparison of modal frequencies (Hz) of 23-m turbine blade with and without rotation effect is shown in Table 4.1. The comparison of modal frequencies (Hz) of 63-m turbine blade with and without rota72 edge mode 1 edge mode 2 0.4 Modal Displacement Modal Displacement 0.2 0 −0.2 −0.4 0 5 10 15 20 Axial location (m) edge mode 3 0.2 0 −0.2 −0.4 0 5 10 15 20 Axial location (m) 0 −0.2 −0.4 25 0 5 10 15 20 Axial location (m) With Rotation Without Rotation edge mode 4 0.4 Modal Displacement Modal Displacement 0.4 0.2 0.2 0 −0.2 −0.4 25 25 0 5 10 15 20 Axial location (m) 25 Figure 4.5 The comparison between the edgewise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 23-m blade when the rotation speed is 1.5 rad/s tion is shown in Table 4.2. The dimensional parameters of a 63-m blade were obtained from the NREL technical report of Jonkman [90]. Figure 4.7 shows the mode shape variation with and without rotation. The comparison of modal frequencies (Hz) of 100-m turbine blade with and without rotation is shown in Table 4.3. The dimensional parameters of a 100-m blade were obtained from the Sandia National Laboratories technical report of Griffith and Ashwill [88]. Figure Table 4.2 Comparison of edgewise modal frequencies (Hz) of 63-m turbine blade between the static case and and the operational case, when Ω is 1.25 rad/s Modes No. static case mode 1 1.06 mode 2 3.95 mode 3 9.17 mode 4 16.74 73 operational case 1.07 3.97 9.19 16.77 edge mode 1 edge mode 2 0.4 Modal Displacement Modal Displacement 0.2 0 −0.2 −0.4 0 5 10 15 20 Axial location (m) −0.2 0 5 With Rotation Without Rotation edge mode 3 10 15 20 Axial location (m) 25 edge mode 4 0.4 Modal Displacement Modal Displacement 0 −0.4 25 0.4 0.2 0 −0.2 −0.4 0.2 0 5 10 15 20 Axial location (m) 0.2 0 −0.2 −0.4 25 0 5 10 15 20 Axial location (m) 25 Figure 4.6 The comparison between the edgewise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 23-m blade when rotation speed is 15 rad/s 4.8 shows the mode shape variation with and without rotation. 4.3.2 Flapwise Bending Modal Analysis The flapwise bending vibration of the wind turbine blade is governed by the equation (4.8). Without considering the aerodynamic force and parametric excitation, the damping matrix is neglected, and then the mass and stiffness matrices are obtained as 74 edge mode 1 edge mode 2 0.4 Modal Displacement Modal Displacement 0.2 0 −0.2 −0.4 0 20 40 Axial location (m) 60 edge mode 3 0.2 0 −0.2 −0.4 0 20 40 Axial location (m) 0 −0.2 −0.4 0 With Rotation Without Rotation 0.4 Modal Displacement Modal Displacement 0.4 0.2 60 edge mode 4 0.2 0 −0.2 −0.4 60 20 40 Axial location (m) 0 20 40 Axial location (m) 60 Figure 4.7 The comparison between the edgewise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 63-m blade when rotation speed is 1.25 rad/s L mij = m(x)ui uj dx 0 L L kij = 0 EI(x)ui uj dx + 0 m(x)Ω2 ui uj dx − L 0 Ω2 (4.11) L x mwdwui uj dx Compared to the edgewise case in equation 4.9, the flapwise stiffness has one less term. Table 4.3 Comparison of edgewise modal frequencies (Hz) of 100-m turbine blade between the static case and and the operational case, when Ω is 0.78 rad/s Modes No. static case mode 1 0.64 mode 2 1.91 mode 3 4.67 mode 4 8.23 75 operational case 0.65 1.93 4.69 8.25 edge mode 1 edge mode 2 0.4 Modal Displacement Modal Displacement 0.2 0 −0.2 −0.4 0 20 40 60 80 Axial location (m) edge mode 3 0.2 0 −0.2 −0.4 0 20 40 60 80 Axial location (m) 0 −0.2 −0.4 100 0 20 40 60 80 Axial location (m) With Rotation Without Rotation edge mode 4 0.4 Modal Displacement Modal Displacement 0.4 0.2 0.2 0 −0.2 −0.4 100 100 0 20 40 60 80 Axial location (m) 100 Figure 4.8 The comparison between the edgewise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 100-m blade when rotation speed is 0.78 rad/s Figure 4.9 shows the variation of the modal frequencies with the increace of the rotation speed. The range of the rotation speed varies from 1 rad/s to 10 rad/s. It is observed that when the rotation speed increases, the modal frequencies increase slightly more than in the edgewise case. However, again the effect of rotation on the frequencies is negligible in the low Ω operating range. The modal frequencies obtained show good agreement with the values published in the report of Bir and Oyague. Then we could include the terms with rotation in the calculation and compare the mode shape under the two conditions. Figure 4.10 shows the mode shape variation with and without rotation. The dashed lines represent the mode shape when the wind turbine blade is stationary. The solid lines represent the mode shape when the wind turbine is operating at 76 25 Modal Frequencies 20 mode1 mode2 mode3 mode4 15 10 5 0 1 2 3 4 5 6 Rotation Speed 7 8 9 10 Figure 4.9 The modal frequencies for flapwise bending vibration of 23-m blade with the variation of rotation speed the rotation speed 1.57 rad/s. We can see there is no observable difference between the modes shapes under the two conditions, but this will change when the rotation speed increases. Figure 4.11 shows the mode shapes when the rotation speed is set to 15 rad/s. We can see from mode 1 that significant difference already shows between the modes without rotation and those with rotation. However, when the size of wind turbine blade increases, the difference between the modal frequencies of the static case and those of the operational case increases. The comparison of modal frequencies (Hz) of 63-m turbine blade with and without rotation at Ω = 1.25m/s is shown in Table 4.4. The dimensional parameters of a 63-m blade were obtained from the NREL technical report of Jonkman [90]. Figure 4.12 shows the mode shape variation with and without rotation. The comparison of modal frequencies (Hz) of 100-m turbine blade with and without 77 flap mode 1 flap mode 2 0.4 Modal Displacement Modal Displacement 0.2 0 −0.2 −0.4 0 5 10 15 20 Axial location (m) flap mode 3 0.2 0 −0.2 −0.4 0 5 10 15 20 Axial location (m) 0 −0.2 −0.4 25 0 5 10 15 20 Axial location (m) With Rotation Without Rotation flap mode 4 0.4 Modal Displacement Modal Displacement 0.4 0.2 0.2 0 −0.2 −0.4 25 25 0 5 10 15 20 Axial location (m) 25 Figure 4.10 The comparison between the flapwise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 23-m blade when the rotation speed is 1.5 rad/s rotation at Ω = 1.25m/s is shown in Table 4.5. The dimensional parameters of a 100-m blade were obtained from the Sandia National Laboratories technical report of Griffith and Ashwill [88]. Figure 4.13 shows the mode shape variation with and without rotation. Table 4.4 Comparison of flapwise modal frequencies (Hz) of 63-m turbine blade between the static case and the operational case when the rotation speed is 1.25 rad/s Modes No. static case mode 1 0.66 mode 2 1.90 mode 3 4.40 mode 4 7.92 78 operational case 0.71 1.96 4.46 7.98 flap mode 1 flap mode 2 0.4 Modal Displacement Modal Displacement 0.2 0 −0.2 −0.4 0 5 10 15 20 Axial location (m) 25 flap mode 3 −0.2 0 5 With Rotation Without Rotation 10 15 20 Axial location (m) 25 flap mode 4 0.4 Modal Displacement Modal Displacement 0 −0.4 0.4 0.2 0 −0.2 −0.4 0.2 0 5 10 15 20 Axial location (m) 0.2 0 −0.2 −0.4 25 0 5 10 15 20 Axial location (m) 25 Figure 4.11 The comparison between the flapwise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 23-m blade when the rotation speed is 15 rad/s 4.4 Study of the Effects of Rotational Position due to Gravity A horizontal-axis rotating blade is affected by its rotational position and the rotational motion. We studied the effects of rotation speed on the modes in section 4.3.1 and section 4.3.2. Relative to the angular position, the blade undergoes changing gravity forces, and so the modes of the blade at different rotation positions could be different. The stiffness of the blade will increase due to added tension in the downward position, and it will decrease due to added compression in the upward position. Therefore, it is insightful to study the effect of rotational effect on the modes of the blade at different rotation positions. 79 flap mode 1 flap mode 2 0.4 Modal Displacement Modal Displacement 0.4 0.2 0 −0.2 −0.4 0 20 40 Axial location (m) 60 flap mode 3 0.2 0 −0.2 −0.4 0 20 40 Axial location (m) 0 −0.2 −0.4 0 With Rotation Without Rotation 0.4 Modal Displacement Modal Displacement 0.4 0.2 60 flap mode 4 0.2 0 −0.2 −0.4 60 20 40 Axial location (m) 0 20 40 Axial location (m) 60 Figure 4.12 The comparison between the flapwise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 63-m blade when rotaton speed is 1.25 rad/s The modes of 23-m wind turbine blade at different rotation positions (downward position and horizontal position) were investigated and compared. As shown in Figure 4.14 and Figure 4.15, there is no significant difference for the edgewise and flapwise modes in horizontal and downward positions. The modal frequencies of the 23-m blade at the downward and horizontal positions have minor differences, as shown in Table 4.6 and Table 4.7. It can be seen that the modal Table 4.5 Comparison of flapwise modal frequencies (Hz) of 100-m turbine blade between the static case and the operational case when the rotation speed is 0.78 rad/s Modes No. static case mode 1 0.51 mode 2 1.54 mode 3 3.38 mode 4 5.87 80 operational case 0.54 1.58 3.46 5.94 flap mode 1 flap mode 2 0.4 Modal Displacement Modal Displacement 0.4 0.2 0 −0.2 −0.4 0 20 40 60 80 Axial location (m) flap mode 3 −0.2 0 20 With Rotation Without Rotation 40 60 80 Axial location (m) 100 flap mode 4 0.4 Modal Displacement Modal Displacement 0 −0.4 100 0.4 0.2 0 −0.2 −0.4 0.2 0 20 40 60 80 Axial location (m) 0.2 0 −0.2 −0.4 100 0 20 40 60 80 Axial location (m) 100 Figure 4.13 The comparison between the flapwise bending modes of the static case (dashed lines) and those of the operational case (solid lines) for 100-m blade when the rotation speed is 0.78 rad/s frequencies of the first four modes have minor differences, around 1%. Next, the modes of 63-m wind turbine blade at different rotation positions (downward position and horizontal position) were investigated and compared. As shown in Figure 4.16 and Figure 4.17, there is no significant difference for the edgewise and flapwise modes at horizontal and downward positions. However, although the differences of mode shapes of the 63-m wind turbine blade are indistinguishable, the modal frequencies of the blades in the downward and horizontal positions have minor differences, as shown in Table 4.8 and Table 4.9. It can be seen that the modal frequencies of the first flapwise mode have a difference around 6.4%. The vertical (downward) frequencies of mode 1 is larger than the horizontal frequencies, because the ef81 edge mode 1 edge mode 2 0.4 0.3 Modal Displacement Modal Displacement 0.4 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) edge mode 3 0.2 0 −0.2 −0.4 0 5 10 15 20 Axial location (m) 0 −0.2 −0.4 25 0 5 10 15 20 Axial location (m) horizontal Position Vertical Position edge mode 4 0.4 Modal Displacement Modal Displacement 0.4 0.2 0.2 0 −0.2 −0.4 25 25 0 5 10 15 20 Axial location (m) 25 Figure 4.14 The edgewise modes comparision for 23-m wind turbine blades when the blade rotates at horizontal and downward positions Table 4.6 Comparison of edgewise modal frequencies (Hz) of the 23-m wind turbine blade at different rotation angles Modes No. downward position horizontal position mode 1 3.189 3.161 mode 2 10.145 10.124 mode 3 23.013 22.978 mode 4 42.988 42.911 fect of the gravity forces increases the stiffiness at downward position and therefore raises the frequency. The 100-m blade has a visible difference between the mode shapes of horizontal position and downward position, as shown in Figure 4.18 and Figure 4.19. The modal frequencies of the blades in the downward and horizontal positions have a significant difference, as shown in Table 4.10 and Table 4.11. It can be seen that the modal frequencies of the first edgewise mode have a difference around 12%, while the modal frequencies of the first flapwise mode have a difference around 9.3%. 82 Table 4.7 Comparison of flapwise modal frequencies (Hz) of the 23-m wind turbine blade at different rotation angles Modes No. downward position horizontal position mode 1 2.017 1.995 mode 2 5.459 5.442 mode 3 11.767 11.728 mode 4 22.096 22.039 Table 4.8 Comparison of edgewise modal frequencies (Hz) of the 63-m wind turbine blade at different rotation angles Modes No. downward position horizontal position mode 1 1.112 1.075 mode 2 4.017 3.978 mode 3 9.293 9.201 mode 4 16.938 16.779 Table 4.9 Comparison of flapwise modal frequencies (Hz) of the 63-m wind turbine blade at different rotation angles Modes No. downward position horizontal position mode 1 0.735 0.688 mode 2 1.975 1.902 mode 3 4.489 4.328 mode 4 7.956 7.738 Table 4.10 Comparison of edgewise modal frequencies (Hz) of the 100-m wind turbine blade at different rotation angles Modes No. downward position horizontal position mode 1 0.688 0.645 mode 2 2.016 1.933 mode 3 4.835 4.694 mode 4 8.516 8.256 Table 4.11 Comparison of flapwise modal frequencies (Hz) of the 100-m wind turbine blade at different rotation angles Modes No. downward position horizontal position mode 1 0.594 0.539 mode 2 1.639 1.581 mode 3 3.604 3.463 mode 4 6.166 5.935 83 flap mode 1 flap mode 2 0.4 0.3 Modal Displacement Modal Displacement 0.4 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) flap mode 3 0.2 0 −0.2 −0.4 0 5 10 15 20 Axial location (m) 0 −0.2 −0.4 25 0 5 10 15 20 Axial location (m) horizontal Position Vertical Position flap mode 4 0.4 Modal Displacement Modal Displacement 0.4 0.2 0.2 0 −0.2 −0.4 25 25 0 5 10 15 20 Axial location (m) 25 Figure 4.15 The flapwise modes comparision for 23-m wind turbine blades when the blade rotates at horizontal and downward positions 4.5 Nonmodal Aerodynamic Damping Wind turbine blades are often tested in a lab, which seldom accounts for rotational and aerodynamic effects. A wind turbine blade in the field can have aerodynamic loads that, if expanded, will have terms proportional to deflection rate, which can classify as linear damping terms, although the proportionality may be x dependent. Thus wind turbines are likely to have non-modal damping. In this section, we derive for the aerodynamic force and introduce this into the EOM. The modal properties of the nonmodally damped system such as damping ratio and modal nonsychronicity are investigated. The goal of this work is to formulate the complex modes of flapwise motion and edgewise motion, due to the non-modal damping generated from a very simple quasi-steady aerodynamic model, see whether the complexity of modes is significant, and compare the results from a variety of different blades. Compared to blades in the lab, blades in the field have additional aeroelastic loads, and also centrifugal and cyclic effects due to rotation. This 84 edge mode 1 edge mode 2 0.4 0.3 Modal Displacement Modal Displacement 0.4 0.2 0.1 0 −0.1 −0.2 0 20 40 Axial location (m) −0.2 0 horizontal Position Vertical Position edge mode 3 20 40 Axial location (m) 60 edge mode 4 0.4 Modal Displacement Modal Displacement 0 −0.4 60 0.4 0.2 0 −0.2 −0.4 0.2 0 20 40 Axial location (m) 0.2 0 −0.2 −0.4 60 0 20 40 Axial location (m) 60 Figure 4.16 The edgewise modes comparision for 63-m wind turbine blades when the blade is in horizontal and downward positions work focuses on the damping due to aeroelasticity based on a simple model, and provides some experience for understanding whether laboratory modal testing can be representative of blades in the field. As the angle of attack of an airfoil increases, the lift coefficient increases until it comes to the stall condition [125]. Airfoils in dynamic stall exhibit large hysteresis loops in lift curve. However, in our case we are concerned with small motions and the stability of equilibrium, which does not involve stall. Measurements of lift coefficients with dynamically varying angles of attack involve hysteresis, such as in Figure 4.20 reproduced from [126]. Stall involves large hysteresis loops, as in the right column of the figure. As we are interested in small oscillations without nonlinear stall, we expect small hysteresis effects as in the left column of the figure. A quasi-steady model neglects this hysteresis. The neglected hysteresis omits some energy added to the system (reducing the effective damping) [125]. As such, the quasi-steady model is a simplifying approximation. Some applications or study of quasi85 flap mode 1 flap mode 2 0.4 0.3 Modal Displacement Modal Displacement 0.4 0.2 0.1 0 −0.1 −0.2 0 20 40 Axial location (m) −0.2 0 20 40 Axial location (m) horizontal Position Vertical Position flap mode 3 60 flap mode 4 0.4 Modal Displacement Modal Displacement 0 −0.4 60 0.4 0.2 0 −0.2 −0.4 0.2 0 20 40 Axial location (m) 0.2 0 −0.2 −0.4 60 0 20 40 Axial location (m) 60 Figure 4.17 The flapwise modes comparision for 63-m wind turbine blades when the blade is in horizontal and downward positions steady model can be found in reference [127–131]. 4.5.1 Flapwise Bending Adding aeroelastic forces to the linearized equation of motion 4.20 for flapwise bending vibration, and omitting the direct excitation term for modal analysis, yields m¨ z + F (z, z) ˙ + (EIz ) + mΩ2 xz − ( L mΩ2 wdw)z x L +mg cos Ωtz − mgdw cos Ωt = 0 (4.12) x where z is the out of plane flexural displacement, x is the location along the axis of the blade, t is the time, the primes indicate derivatives with respect to space, x. F (z, z) ˙ = F0 (x) + A(x)z˙ + B z˙ 2 are the aerodynamic forces per unit length from a quasistatic lift/drag model, as derived below. In this section we isolate the effect of nonmodal damping, i.e. we 86 edge mode 1 edge mode 2 0.4 0.3 Modal Displacement Modal Displacement 0.4 0.2 0.1 0 −0.1 −0.2 0 20 40 60 80 Axial location (m) −0.2 0 20 horizontal Position Vertical Position edge mode 3 40 60 80 Axial location (m) 100 edge mode 4 0.4 Modal Displacement Modal Displacement 0 −0.4 100 0.4 0.2 0 −0.2 −0.4 0.2 0 20 40 60 80 Axial location (m) 0.2 0 −0.2 −0.4 100 0 20 40 60 80 Axial location (m) 100 Figure 4.18 The edgewise modes comparision for the 100-m wind turbine blade when the blade are at horizontal and downward positions omit for now the effect of excitation and linearize about equilibrium. We take a cross section from the chord of the blade, and consider the aerodynamic forces per unit length which the blade endures dummy operation, as shown in Figure 4.21. L is the lift force, D is the drag force, vrel is the relative velocity of the wind with respect to the blade, β is the pre-twist angle and θ is relative wind angle. We will first consider the aerodynamic force projected in the flapwise bending direction. The relative velocity with flapwise deformation is vrel = vwind −vblade = −uk−(Ωxeφ + zk)), ˙ where Ω is the rotating speed of the hub, u is the wind speed, k is the unit vector of the coordinate in the z direction, and eφ is the unit vector of the coordinate in the x direction tangential to the circular path of rotation. 87 flap mode 1 flap mode 2 0.4 0.3 Modal Displacement Modal Displacement 0.4 0.2 0.1 0 −0.1 −0.2 0 20 40 60 80 Axial location (m) flap mode 3 −0.2 0 20 horizontal Position Vertical Position 40 60 80 Axial location (m) 100 flap mode 4 0.4 Modal Displacement Modal Displacement 0 −0.4 100 0.4 0.2 0 −0.2 −0.4 0.2 0 20 40 60 80 Axial location (m) 0.2 0 −0.2 −0.4 100 0 20 40 60 80 Axial location (m) 100 Figure 4.19 The flapwise modes comparision for the 100-m wind turbine blade when the blade are at horizontal and downward positions Then the magnitudes of lift and drag forces per unit length are derived as [101] 1 2 L = ρCL (α)cVrel 2 1 2 D = ρCD (α)cVrel 2 (4.13) where ρ is the air density, c is the chord length, CL (α) is the lift coefficient, and CD (α) is the drag coefficient. CL (α) = CL0 + CL1 α + · · · , CD (α) = CD0 + CD1 α2 + · · · , where CL0 , CL1 , CD0 and CD1 are constants generated from experiments and could be obtained in technical reports. Since the blade cross-section varies along its length, these quantities 88 Figure 4.20 A comparison plot of lift coefficient between experiment data and ONERA dynamic stall model are all functions of x. Considering the geometry, tan θ = u + z˙ xΩ (4.14) z˙ where β is the pretwist angle, and so the angle of attack is α = β + θ = β + tan−1 u+ xΩ . The aerodynamic force per unit length projected in the k direction for flapwise bending is expressed as Fk = L cos θ + D sin θ 1 2 cos θ + 1 ρC (α)cV 2 sin θ = ρCL (α)cVrel rel 2 2 D 1 1 = ( ρCL (α)cΩx + ρCD (α)c(u + z)) ˙ (u + z) ˙ 2 + (Ωx)2 2 2 (4.15) Expanding in Taylor series about z˙ and dropping the higher order terms, we obtain (u + z) ˙ 2 + (Ωx)2 ≈ u2 + (Ωx)2 + u u2 + (Ωx)2 z˙ (4.16) Substituting equation (4.16) back into equation (4.15), and taking the linear part, we 89 Figure 4.21 Translating airfoil with lift and drag forces acting obtain Fk = Fk0 + Fk1 z˙ + · · · , where 1 (Ωx)2 u 1 Fk1 = [( ρcCL1 + ρc(CD0 + CD1 (β + tan−1 ( ))2 ) u2 + (Ωx)2 2 2 2 2 xΩ (Ωx) + u 1 1 u + ( ρcCL0 Ωx + ρcCL1 Ωx(β + tan−1 ( )) 2 2 xΩ u u2 + (Ωx)2 (4.17) ]z˙ An assumed modal expansion is used to reduce the order of the linearized equation of motion. According to the assumed-modes method, the transverse displacement is expanded as z(x, t) ≈ N i=1 ui (x)qi (t), where ui (x) are the assumed modes derived from the uniform cantilevered beam, and qi (t) are the assumed modal coordinates. This expression is substituted into the equation of motion (4.12). Multiplying by uj (x), and integrating over the length of the turbine blade, we get a set of second-order ODEs as 90 M¨ q + Cq˙ + Kq = Q0 (4.18) where Q0 is due to the constant term in the aerodynamic force. Then we shift the coordinates qˆ = q − K−1 Q0 , drop the hat, and have equation 4.18 about the static deflection. Omitting the parametric excitation terms for now, such that we focus on the structural and aeroelastic contributions to the modes, the elements of the N × N matrices M, K, and C are L mij = m(x)ui uj dx, 0 L L kij = EI(x)ui uj + 0 0 (Ωx)2 L L m(x)Ω2 ui uj dx − 0 Ω2 L x mwdwui uj dx 1 1 u 2 [( ρcCL1 + ρc(C + C (β + arctan( )) ) u2 + (Ωx)2 D D 2 + u2 0 1 2 2 xΩ (Ωx) 0 1 1 u u + ( ρcCL0 Ωx + ρcCL1 Ωx(β + arctan( )) ]ui uj dx 2 2 xΩ u2 + (Ωx)2 cij ≈ (4.19) for i = 1 · · · , N and j = 1 · · · , N . 4.5.2 Edgewise Bending The equation of motion for edgewise bending vibration is obtained as m¨ y + F (y, y) ˙ y˙ + (EIy ) − mΩ2 y + mΩ2 xy L −( mΩ2 wdw)y + mg cos Ωty − x L mgdw cos Ωt = 0 x 91 (4.20) For aerodynamic loading in the edgewise direction under edgewise deflection, the relative ˙ φ )), where y(x, t) is the edgewise velocity is vrel = vwind − vblade = −uk − (Ωxeφ + ye deflection displacement. Considering the geometry tan θ = u xΩ + y˙ (4.21) Then the linear part of the aerodynamic force projected in the eφ direction for edgewise bending is derived as Feφ = L sin θ − D cos θ 1 2 sin θ − 1 ρC (α)cV 2 cos θ = ρCL (α)cVrel rel 2 2 D 1 1 ˙ = ( ρCL (α)cu − ρCD (α)c(Ωx + y)) (Ωx + y) ˙ 2 + u2 2 2 (4.22) Expanding the nonlinear terms in Taylor series in y, ˙ and dropping the higher order terms, we get 1 1 1 ≈ − y˙ Ωx + y˙ Ωx (Ωx)2 (u + y) ˙ 2 + (Ωx)2 ≈ u2 + (Ωx)2 + Ωx u2 + (Ωx)2 (4.23) y˙ Applying to equation 4.22, we get Feφ = Feφ + Feφ y˙ + · · · , where 0 1 92 1 u2 1 u Feφ = [( ρcCL1 + ρc(CD0 + CD1 (β + arctan( ))2 ) 2 2 1 2 2 Ωx (Ωx) + u u2 + (Ωx)2 u 1 1 − ( ρcu(CL0 + CL1 (β + arctan( ))) + ρcΩx 2 Ωx 2 Ωx u 2 ((CD0 + CD1 (β + arctan( )) ))) ]y˙ Ωx u2 + (Ωx)2 (4.24) Applying the assumed modal expansion described in the previous section, and casting the modal coordinates about the equilibrium due to the constant Feφ term, the equation of 0 motion about equlibrium are M¨ q + Cq˙ + Kq = 0, where the elements of the N × N matrices M, K, and C are obtained as L mij = m(x)ui uj dx, 0 L L kij = EI(x)ui uj + 0 0 m(x)Ω2 ui uj L dx − 0 Ω2 L x mwdwui uj dx L 1 1 u2 u [( ρcCL1 + ρc(CD0 + CD1 (β + arctan( ))2 ) u2 + (Ωx)2 2 2 2 2 Ωx (Ωx) + u 0 1 u 1 − ( ρcu(CL0 + CL1 (β + arctan( ))) + ρcΩx 2 Ωx 2 u 2 Ωx ((CD0 + CD1 (β + arctan( )) ))) ]ui uj dx Ωx u2 + (Ωx)2 cij = (4.25) for i = 1 · · · , N and j = 1 · · · , N . 4.5.3 Effect of Equilibrium on Stiffness In the previous derivation, we dropped the constant terms Fk0 and Feφ . Here we do a 0 detailed analysis to check the effect of the constant terms on the equilibrium of the system, 93 and in turns, the effect on the linearized systems. We start with the equation of motion, choosing edgewise bending as an example as m¨ y + L(y) + fs (y) + f0 + fa (y) ˙ =0 (4.26) where fs (y) represents the nonlinear stiffness terms in the equation of motion, where y is a list of y and its derivatives, as needed in the expression, and fa (y) ˙ represents the nonlinear terms from the aerodynamic damping. Letting y˙ = y¨ = 0, we obtain the equilibrium y0 (x) from L(y0 ) + fs (y0 ) + f0 + fa (0) = 0 (4.27) Considering small motions about equilibrium, we let u = y−y0 and insert it into equation (4.26) to produce m¨ u + L(y0 + u) + fs (y0 + u) + f0 + fa (u) ˙ =0 (4.28) Expanding the nonlinear terms in a Taylor series, we obtain m¨ u + L(y0 ) + L(u) + fs (y0 ) + Dfs |y0 u + o(2) + f0 + fa (0) + Dfa |0 u˙ + o(2) = 0 (4.29) Then with the equilibrium equation (4.27), we obtain m¨ u + L(u) + Dfs |y0 u + Dfa |u0 u˙ = 0 94 (4.30) The third term is a stiffness term, and the fourth term is a damping term. Assuming u and its derivatives are small, the terms Dfs |y0 u have a small effect on the system stiffness, and so they are neglected as we study the modal issues associated with the damping term, Dfa |u0 u. ˙ 4.5.4 Modal Analysis Results Equation (4.18) can be put in state-variable form by letting yT = [q˙ T qT ], and introducing the equation Mq˙ − Mq˙ = 0 [110]. As such, Ay˙ + By = 0 (4.31)      0 M −M 0  where matrices A =  , B =    are 2N × 2N . Seeking a solution of the M C 0 K form y(t) = φeαt leads to a general eigenvalue problem, αAφ + Bφ = 0 (4.32) Solving the eigenvalue problem and recombining the eigenvectors φj with the assumed modal functions ui (x) in the way described in Chapter 2, allows us to approximate the modes. Eight assumed modes were used in the assumed-modes method, as we applied it to three example blades: the 23-m NREL, the 63-m NREL and the 100-m Sandia blade. Each of these blades has a varying profile with x. For simplicity, we assumed a constant airfoil type, NACA-64 airfoil. The dimensional parameters of a 23-m blade were obtained from the NREL technical report of Bir and Oyague [89]. The parameters used in the calculation are: CL0 = 0.58, CL1 = 3.1, CD0 = 0.1, CD1 = 8.25, ρ = 1.18 kg/m3 , Ω = 1.57 rad/s 95 and u = 16m/s. With these parameters, the tip speed is ΩL = 36.6 m/s, and the tip speed ratio is ΩL/u = 2.29. Considering the complex modal pair as a “mode”, the associated eigenvalues can be expressed in the form α1,2 = −ζωm ±iωm ζ 2 − 1, from which a “modal” frequency ωm and a “nonmodal damping ratio” ζ were computed. The modal frequencies and damping ratio obtained under the aerodynamic damping are shown in Table 4.12, where ωn is undamped modal frequency. The first four modes including the real and imaginary parts of the complex modes of flapwise bending are shown in Figure 4.23. The damping ratios of the flapwise modes are small, and the modes are dominated by the real part, so the turbine blade can be well approximated by real modes of stationary beam. This can be verifed by the comparision between the modes of the nonmodally damped case and those of the undamped case, as shown in Figure 4.22. The dotted lines in the plot represent for the undamped mode, and match with the solid lines very well, which represent for the imaginary part of the complex mode in the nonmodally damped case. The first four modes of the edgewise bending for the 23-m wind turbine blade are shown in Figure 4.24, and the modal frequencies and damping ratio are shown in Table 4.13. The damping ratios of the edgewise modes are smaller compared to those of the flapwise modes. Table 4.12 Damping ratio and modal frequencies compared to the undamped modal frequencies for flapwise bending of 23-m turbine blade at rated wind speed u = 16 m/s and Ω = 1.57 rad/s Modes No. mode 1 mode 2 mode 3 mode 4 ζ 0.071 0.023 0.009 0.004 α 0.892±12.450i 0.805±34.104i 0.670±73.608i 0.606±138.40i ωm ωn 12.482 12.533 34.113 34.194 73.611 73.688 138.41 138.476 The effect of nonmodal damping is more significant on the 63-m wind-turbine blade. The 96 flap mode 1 flap mode 2 0.3 Real and imaginary parts Real and imaginary parts 0.4 0.3 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) 25 flap mode 3 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) 0.1 0 −0.1 −0.2 0 Real Imaginary undamped 0.3 Real and imaginary parts Real and imaginary parts 0.3 0.2 25 5 10 15 20 Axial location (m) 25 flap mode 4 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) 25 Figure 4.22 The comparision between the modes of the nonmodally damped case and those of the undamped case for 23-m blade dimensional parameters of a 63-m blade were obtained from the NREL technical report of Jonkman [90]. For an operating speed Ω = 1.25 rad/s and wind speed u = 11.4 m/s, the tip speed is ΩL = 78.7 m/s, and the tip speed ratio is ΩL/u = 6.9. The first four modes of flapwise bending are shown in Figure 4.25, and the modal frequencies and damping ratio are shown in Table 4.14. The modal frequencies and damping ratio of edgewise bending are listed in Table 4.15. Compared to the 23-m wind turbine blade, the damping ratio of the Table 4.13 Damping ratio and modal frequencies compared to the undamped modal frequencies for edgewise bending of 23-m turbine blade at rated wind speed u = 16 m/s and Ω = 1.57 rad/s Modes No. mode 1 mode 2 mode 3 mode 4 ζ 0.031 0.009 0.004 0.002 α 0.608±19.889i 0.617±63.596i 0.549±144.35i 0.475±269.59i 97 ωm ωn 19.898 19.860 63.599 63.734 144.35 144.372 269.59 269.913 flap mode 2 Real and imaginary parts Real and imaginary parts flap mode 1 0.4 0.3 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) 25 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) 0.2 0.1 0 −0.1 −0.2 0 Real Imaginary 0.3 Real and imaginary parts Real and imaginary parts flap mode 3 0.3 0.3 25 5 10 15 20 Axial location (m) 25 flap mode 4 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) 25 Figure 4.23 The first four modes of the flapwise bending of 23-m wind turbine blade with nonmodal damping flapwise modes increases. The aerodynamic loading exerts a significant damping on the large wind-turbine blade. The first four modes of edgewise bending are shown in Figure 4.26, and the modal frequencies and damping ratio are listed in Table 4.15. The damping ratios of the first four modes increase compared to the values of 23-m blade, but are much smaller than those of the flapwise modes. Table 4.14 Damping ratio and modal frequencies compared to the undamped modal frequencies for flapwise bending of 63-m turbine blade at rated wind speed u = 11.4 m/s and Ω = 1.25 rad/s Modes No. mode 1 mode 2 mode 3 mode 4 ζ α ωm ωn 0.173 0.732±4.165i 4.229 4.457 0.048 0.581±12.014i 12.028 12.325 0.019 0.536±27.757i 27.762 28.048 0.010 0.505±49.852i 49.855 50.141 98 edge mode 2 Real and imaginary parts Real and imaginary parts edge mode 1 0.4 0.3 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) 25 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) 0.2 0.1 0 −0.1 −0.2 0 Real Imaginary 0.3 Real and imaginary parts Real and imaginary parts edge mode 3 0.3 0.3 25 5 10 15 20 Axial location (m) 25 edge mode 4 0.2 0.1 0 −0.1 −0.2 0 5 10 15 20 Axial location (m) 25 Figure 4.24 The first four modes of the edgewise bending of 23-m wind turbine blade with nonmodal damping The trend continues with the 100-m wind-turbine blade. When the length of the blade increases, the damping ratio increases. The dimensional parameters of a 100-m blade were obtained from the Sandia National Laboratories technical report of Griffith and Ashwill [88]. For an operating speed Ω = 0.78 rad/s and wind speed u = 11.3 m/s, the tip speed is ΩL = 78 m/s, and the tip speed ratio is ΩL/u = 6.9. The first four modes are shown in Figure 4.27, and the modal frequencies and damping ratio are listed in Table 4.16. The modal Table 4.15 Damping ratio and modal frequencies compared to the undamped modal frequencies for edgewise bending of 63-m turbine blade at rated wind speed u = 11.4 m/s and Ω = 1.25 rad/s Modes No. mode 1 mode 2 mode 3 mode 4 ζ α ωm ωn 0.057 0.384±6.653i 6.664 6.748 0.015 0.374±24.860i 24.863 24.979 0.006 0.359±57.664i 57.665 57.783 0.003 0.351±105.25i 105.25 105.37 99 flap mode 1 flap mode 2 0.3 Real and imaginary parts Real and imaginary parts 0.4 0.3 0.2 0.1 0 −0.1 −0.2 0 20 40 Axial location (m) −0.1 −0.2 0 Real and imaginary parts Real and imaginary parts 0 Real Imaginary 0.3 0.3 0.2 0.1 0 −0.1 −0.2 20 40 Axial location (m) 0.1 60 flap mode 3 0 0.2 60 20 40 Axial location (m) 60 flap mode 4 0.2 0.1 0 −0.1 −0.2 0 20 40 Axial location (m) 60 Figure 4.25 The first four modes of the flapwise bending of 63-m wind turbine blade with nonmodal damping frequencies and damping ratio of edgewise bending are shown in Table 4.17. Compared to 63-m blade, the damping ratios increase and the value of the first flapwise mode reaches 0.25. Comparing the motion-dependent angle of attack in flapwise bending and edgewise bendz˙ −1 u + ing, in flapwise bending, α = β + θ = β + tan−1 u+ xΩ = tan xΩ xΩ z. ˙ (Ωx)2 +u2 In edgewise u u −1 u − bending, α = β + θ = β + tan−1 xΩ+ y. ˙ As a rough estimate, the y˙ = tan xΩ (Ωx)2 +u2 coefficient of z˙ compared to the coefficient of y˙ is increasing with the increase of the length of the blade, i.e. the tip speed ratio increases, and thus the angle of attack of the flapwise direction varies more than that of the edgewise direction. Therefore, this implies a greater variation in the force with changes of velocity in the flapwise motion, and so the nonmodal damping in flapwise direction is much stronger than that of the edgewise direction. 100 edge mode 1 edge mode 2 0.2 Real and imaginary parts Real and imaginary parts 0.4 0.3 0.2 0.1 0 −0.1 −0.2 0 20 40 Axial location (m) 60 edge mode 3 0.1 0 −0.1 −0.2 0 20 40 Axial location (m) 0 −0.1 −0.2 0 Real Imaginary 0.2 Real and imaginary parts Real and imaginary parts 0.2 0.1 60 edge mode 4 0.1 0 −0.1 −0.2 60 20 40 Axial location (m) 0 20 40 Axial location (m) 60 Figure 4.26 The first four modes of the flapwise bending of 63-m wind turbine blade with nonmodal damping A “traveling index” [9] indicates the degree of nonsynchronicity of the real and imaginary parts of the complex modes. The traveling index is defined as the reciprocal of the condition number of the matrix whose two columns are the real and imaginary components of the complex mode. The traveling index of the flapwise bending modes of the three blades are listed in Table 4.18. The values of mode 2 are higher compared to the values of the other Table 4.16 Damping ratio and modal frequencies compared to undamped modal frequencies for flapwise bending of 100-m turbine blade at rated wind speed u = 11.3 m/s and Ω = 0.78 rad/s Modes No. mode 1 mode 2 mode 3 mode 4 ζ α ωm ωn 0.256 0.845±3.197i 3.307 3.382 0.072 0.711±9.782i 9.808 9.931 0.029 0.633±21.624i 21.633 21.747 0.015 0.559±37.151i 37.156 37.272 101 flap mode 1 flap mode 2 0.3 Real and imaginary parts Real and imaginary parts 0.4 0.3 0.2 0.1 0 −0.1 −0.2 0 20 40 60 80 Axial location (m) 0.1 0 −0.1 −0.2 100 0 20 Real Imaginary flap mode 3 40 60 80 Axial location (m) 100 flap mode 4 0.3 Real and imaginary parts 0.3 Real and imaginary parts 0.2 0.2 0.1 0 −0.1 −0.2 0 20 40 60 80 Axial location (m) 100 0.2 0.1 0 −0.1 −0.2 0 20 40 60 80 Axial location (m) 100 Figure 4.27 The first four modes of the flapwise bending of 100-m wind turbine blade with nonmodal damping three modes, which means mode 2 is more nonsynchronous. With the increase of the size of the blade, the modal nonsynchronicity increases. However, in all cases, the modes are highly synchronous. The wind speed is varying in a wide range in the real operation of wind turbines. The typical cut-in and cut-out wind speeds are 3 m/s and 25 m/s [88–90]. Figure 4.29 shows the variation of the nonmodal damping ratios of flapwise bending for the 23-m wind turbine blade. Figure 4.30 shows the variation of the nonmodal damping ratios of flapwise bending for the 100-m wind turbine blade. We can see that the damping ratio increases with the increase of the wind speed for both of the two blades. The damping ratio of 23-m blade increases more with wind than that of the 100-m blade, since the Ωx term in the denominator increases 102 edge mode 1 edge mode 2 0.2 Real and imaginary parts Real and imaginary parts 0.4 0.3 0.2 0.1 0 −0.1 −0.2 0 20 40 60 80 Axial location (m) −0.1 0 20 Real Imaginary edge mode 3 40 60 80 Axial location (m) 100 edge mode 4 0.2 Real and imaginary parts Real and imaginary parts 0 −0.2 100 0.2 0.1 0 −0.1 −0.2 0.1 0 20 40 60 80 Axial location (m) 100 0.1 0 −0.1 −0.2 0 20 40 60 80 Axial location (m) 100 Figure 4.28 The first four modes of the edgewise bending of 100-m wind turbine blade with nonmodal damping from 23-m blade to 100-m blade. 4.5.5 Quantification of the Nonmodal Damping Nonsynchronous modes occur when the damping is significantly nonmodal, i.e. non-diagonalizable by the undamped modal coordinate transformation, but also when the damping is significant. Here we want to define criteria to quantify the degree to which the damping is nonmodal and significant. The undamped modal coordinate transformation gives UT MU = I, UT KU = Λ, UT CU = D, where U is the modal matrix deriving from the eigenvalue problem of the mass and stiffness matrices, K − ω 2 M = 0. If the damping is modal, then D is diagonal. 103 Table 4.17 Damping ratio and modal frequencies compared to undamped modal frequencies for edgewise bending of 100-m turbine blade at rated wind speed u = 11.3 m/s and Ω = 0.78 rad/s Modes No. mode 1 mode 2 mode 3 mode 4 ζ α ωm ωn 0.066 0.263±3.956i 3.965 4.053 0.021 0.256±12.045i 12.047 12.140 0.009 0.262±29.398i 29.399 29.478 0.005 0.264±51.768i 51.769 51.850 Table 4.18 Traveling index for flapwise bending of the three blades at rated wind speed Modes No. 23 − m 63 − m mode 1 0.009 0.015 mode 2 0.024 0.046 mode 3 0.014 0.023 mode 4 0.008 0.015 100 − m 0.021 0.057 0.030 0.021 To quantify the size of the damping terms relative to the stiffness terms, we define a “nonmodal ratio” as the summation of the singular values of matrix D over the summation of the singular values of matrix Λ, where the singular values quantify the amplitude of a matrix. The “nonmodal ratio” of the three blades are shown in Table 4.19. It shows that with the increase of the size of the blades, the amplitude of the damping matrix becomes larger compared to the stiffness matrix, which correlates with the variation of the damping ratio of the three blades. We then define a “nondiagonality” to quantify how nondiagonal the matrix D is. Each eigenvector uj of matrix D has angles with each basis vector. For each eigenvector, we find the minimum angle of the vectors uj to the basis vectors. Then define the “nondiagonality” as the maximum over the eigenvectors of these minimal angles. Max cosine corresponds to Table 4.19 “Nonmodal ratio” for flapwise bending of the three blades at rated wind speed Blade “Nonmodal ratio” 23 − m 63 − m 100 − m 1.46e-6 1.26e-5 2.93e-5 104 0.09 0.08 Flapwise Damping Ratio 0.07 0.06 mode1 mode2 mode3 mode4 0.05 0.04 0.03 0.02 0.01 0 0 5 10 15 Wind Speed (m/s) 20 25 Figure 4.29 The variation of flapwise bending damping ratio with varying wind speed for 23-m blade Table 4.20 “Nondiagonality” for flapwise bending of the three blades at rated wind speed Blade 23 − m “Nondiagonality” 1.052 63 − m 100 − m 1.055 1.081 minimum angle. For a given uj , we seek the max cosine to represent the minimum angle to a basis vector of D. We then take the maximum of these angles (minimum of cosines), over all eigenvectors uj . Mathematically, this is mini (maxi (abs(uj i))). If the matrix D is diagonal, then its eigenvectors line up with the vector basis of D. The “nondiagonality” of the three blades at rated wind speed are shown in Table 4.20. It shows that with the increase of the size of the blades, the value of the “nondiagonality” increases, i.e. the nondiagonality of the matrix D increases. 105 0.3 Flapwise Damping Ratio 0.25 0.2 mode1 mode2 mode3 mode4 0.15 0.1 0.05 0 0 5 10 15 Wind Speed (m/s) 20 25 Figure 4.30 The variation of flapwise bending damping ratio with varying wind speed for 100-m blade 4.5.6 A Review of Unsteady Dynamic Stall Lift Model There are many methods that could be used for the prediction of the aerodynamics of unsteady lift model of wind turbine blade, such as Theodorsen theory [132], ONERA method [133–135], Leishman-Beddoes method [100] and Larsen’s approach [126]. Theodorsen’s theory is a tool for modelling the sectional aerodynamics using unsteady thin aerofoil theory. It simulates the quasi-periodic first harmonic variations in angle of attack. Added mass and wake vorticity effects are included compared to quasi-steady thin airfoil model [136, 137]. Some comparison between Theodorsen’s unsteady theory and measurements of the unsteady lift on airfoils can be found in [136]. In [137], the Theodorsen’s unsteady aerodynamic model was cast into a low-dimensional, state-space representation. 106 1 ¨ + α˙ − a¨ CL = π[h α] + 2π[α + h˙ + α( ˙ − a)]C(k) 2 (4.33) where the first term in the equation represents for the added-mass effect, while the second term accounts for the lift attenuation by the wake vorticity. C(k) is the transfer function, which can be expressed in Hankel functions. The ONERA model describes the unsteady dynamic stall with a set of differential equations. Tran and Petot [133,134] first presented a system of differential equations relating the aerodynamic forces and the variables of the airfoil section to simulate the time delay effects of the flow. The parameters in the equations were obtained by identification of the test results (engineering tests in wind tunnel), so the ONERA model is a semi-empirical model. Peters [135] then presented a modification of the ONERA model, including the pitch-plunge distinction, the unsteady free-stream, and large angles of attack. For small angles of attack, the model reduces to Theodorsen theory for steady free stream. In the ONERA model, the relation between the lift coefficient of the airfoil and the angle of attack can be expressed by a set of equations as follows dθ d2 θ dCl1 + λCl1 = λaθ + (λs + δ) + s 2 dτ dτ d τ 2 d Cl2 dC d∆Cz dθ +α ¯ γ¯ l2 + γ¯ 2 (1 + α ¯ )2 Cl2 = −¯ γ 2 (1 + α ¯ )2 [∆Cz + C¯ ] 2 dτ dθ dτ d τ (4.34) Cl = Cl1 + Cl2 where Cl1 and Cl2 are the lift coefficients in the linear and non-linear regions of angle of attack, θ is the total angle of attack, and ∆Cl is the difference between the extended linear lift curve and the acual static lift curve [138]. The parameters λ, a, δ, α ¯ , γ¯ , C¯ are functions of 107 angle of attack and can be obtained from wind tunnel tests by parameter identification. A more recent model was presented by Larsen [126], which is based on a backbone curve in the form of the static lift as a function of angle of attack. The static lift is described by two parameters, the lift at fully attached flow and the degree of attachment. The nonstationary effects are formulated by three mechanisms, a delay of the lift coefficient of fully attached flow, a delay of the development of separation, and a lift contribution due to leading edge separation. The total lift coefficient can be expressed as cL (t) = cL,d (t) + cL,v (t) 1 = cos4 ( θd (t))[cL0 (α) − c1 (t) − c2 (t)] + cL,v (t) 4 (4.35) where the first term cL,d (t) represents for the time delay of the lift coefficient, and the second term cL,v (t) represents for the value of the induced lift after the initiation of the diminishing effect. The static coefficient to be determined is cL0 (α), which can be found from experimental static lift coefficient. We focus on linear modal analysis in the nonmodal damping. In this case, the variation of the hesteresis phenomena could be less significant. Therefore, we apply the quasi-steady model first for the formulation of the nonmodal damping, and the unsteady lift models could be brought into the formulation in the future. 4.5.7 Conclusion In this work, the quasi-steady aerodynamic loading in both flapwise and edgewise directions were formulated and the associated damping terms were derived for analysis. We isolated 108 the effect of nonmodal damping in the equation of motion for a rotating wind-turbine blade and further solved for the modes in flapwise and edgewise bending. The data from the NREL technical reports were applied in the computations and the samples of 23-m, 63-m and 100-m wind turbine blades were analyzed. The plots of the first four modes and the calculation of damping ratios show that the aerodynamic loading has a stronger damping effect on the flapwise bending compared to that of the edgewise bending. With the increase of the length of wind turbine blades, the damping ratio and nonsynchronicity increase. The aerodynamics loading exerts a more significant effect on large size wind turbine blades. The damping ratio increases with the increase of the wind speed. The high synchronicity of the modes suggests that the modes can be approximated well by real modes. Therefore, considering aerodynamics (but not rotation effects), modal studies in labs on fixed blades should be representative of modal shapes of aerodynamically loaded blades. The aerodynamic loading modeling of the wind turbine blade in this paper used a simple quasisteady aeroelastic model, and it could be improved by applying a dynamic model, incorporating the unsteady aerodynamic effects such as stall hysteresis. We also approximated the calculations by using a uniform airfoil type in each blade. In our next calculations, we will apply the spatially variable airfoil type as indicated in the technical reports for each blade studied here. 109 Chapter 5 Modal Analysis of Parametrically Excited Systems 5.1 Chapter Introduction In time-periodic systems, such as wind turbines, the time-periodic stiffness matrix will bring difficulties to modal analysis. In experimental modal analysis, there are a few methods, such as Floquet modal analysis, for analyzing the time varying modes based on data measurements, as mentioned in the introductory chapter. However, there is no mature theoretical method to analyze the effect of the parametric excitation on the modal response of the stationary system. So here we develop a perturbation approach to analyze the perturbation effect of the parametric excitation on the unperturbed system. 5.2 Perturbation Approach for Nonresonant Response Basically one obtains the response of the nonlinear system by perturbing the response of the linear system which is obtained by deleting all the nonlinear terms [139]. There are a few ways to implement the perturbation methods, such as straightforward expansion, multiple scales, harmonic balance and averaging method. Here we use straightforward expansion method to solve for non-resonant responses. Since wind turbines are defined to operate well 110 below the first modal frequency, in this intial work, we look at the nonresonant condition. A small, dimensionless parameter is introduced, which is the order of the amplitude of the motion and can be used as a crutch in obtaining the approximate solution. So we start with the equation of motion, which can be written as follows: M¨ q + (K0 + K1 cos(Ωt))q = 0 (5.1) where the dimensions of M, K0 , K1 , q are N × N , N × N , N × N and N × 1, respectively. We assume the solution of the response can be represented by an expansion having the form q = q 0 + q1 (5.2) Substituting equation (5.2) into equation (5.1), we get ˜ 1 (t)q0 = 0 ¨ 1 ) + K0 q0 + K0 q1 + K M(¨ q0 + q (5.3) ˜ 1 (t) = 1 K1 (eiΩt + e−iΩt ) with K 2 Equating the coefficients of , we obtain Order 0 M¨ q0 + K0 q0 = 0 (5.4) ˜ 1 (t)q0 M¨ q1 + K0 q1 = K (5.5) Order 1 Equation 5.4 can be solved assuming synchronous motion, leading to the familiar eigen- 111 value problem, (K0 − ω02 M)φ = 0, which has eigenvectors φj and eigenvalues ω0j , j = 1, · · · , N . Therefore, equation (5.4) has a modal solution of the form 1 i(ω t+β) q0j = a(e 0j + c.c.)φj 2 (5.6) 2 M)φ = 0, j = 1, · · · , N , where φj and ω0j are eigenvectors and eigenvalues of (K0 − ω0j j and a and β regulates the amplitude and phase of the mode. Substituting equation (5.6) into equation (5.5), we get 1 i(Ω+ω0j )t i(Ω−ω0j )t M¨ q1 + K0 q1 = AK1 φj (e +e + c.c.) 4 (5.7) where A = 21 aeiβ . We introduce ωij as ω1j = Ω + ω0j (5.8) ω2j = Ω − ω0j We seek a particular solution to equation (5.7) of the form 2 q1j = iω t χij e ij + c.c. (5.9) i=1 and we obtain 1 2 M)−1 K φ χij = A(K0 − ωij 1 j 4 (5.10) 2 M is invertible, which means ω is not equal to a natural frequency of provided K0 − ωij ij the unperturbed system. Substituting equation (5.10), (5.9) and (5.6) into equation (5.2), 112 we obtain the “modal” solution of the response iω0j t qj = A(e 1 + c.c.)φj + A 4 2 2 M)−1 K φ eiωij t + c.c. (K0 − ωij 1 j (5.11) i=1 2 = ω 2 . Under this condition, from equation (5.8), the analysis breaks down when where ωij 0k Ω = 0 and Ω = |ω0j ± ω0k |. In this straightforward perturbation scheme, a and β in the constant A = 12 aeiβ are determined by initial conditions. Thus, the j th perturbed nonresonant modal response consists of an unperturbed modal re2 M)−1 K φ sponse of shape φj plus a series of perturbations of the shapes χij = 14 A(K0 −ωij 1 j modulated at frequencies ωij = Ω ± ω0j . The net response is quasiperiodic of fluctuating shape. So at this point, we have obtained the non-resonant response of the time periodic system, in the form of a system of original stationary equation of motion solutions plus a series of corrected solutions from the parametric excitation perturbation. The result can show the perturbation effect of the parametric excitation on a specific linear mode. For each particular mode, we can plot the solution in equation (5.11) to visualize how the perturbed part affects the modes of the original system. In the application of wind turbine, the rotation speed Ω is much smaller compared to the modal frequencies of the blade ωi . A detailed comparison of the Ω and the first modal frequency ω1 of edgewise and flapwise bending can be obtained in Table 5.1. Therefore, we can focus on the non-resonant case of the perturbation response. 113 Table 5.1 Comparison of the Ω and the first modal frequency ω1 of edgewise and flapwise bending (Unit:rad/s) Blade Ω Edgewise ω1 Flapwise ω1 5.3 23 − m 1.57 19.86 12.53 63 − m 1.25 6.75 4.46 100 − m 0.78 4.05 3.38 Example: Three Mass System We start with a three mass system as an example of the perturbation analysis on a multidegree-of-freedom system. Figure 5.1 shows an example of the three mass system. Figure 5.1 The example of three mass system The properties of this three degree of freedom system are assigned as follows: m1 = m2 = m3 = m, k1 = k2 = k3 = k4 = k. The equations of motion are obtained as m1 x¨1 + k1 x1 + k2 (x1 − x2 ) = 0 (5.12) m2 x¨2 + k2 (x2 − x1 ) + k3 (x2 − x3 ) = 0 (5.13) m3 x¨3 + k3 (x3 − x2 ) + k4 x3 = 0 (5.14) Combining all three equations, we can obtain the matrix form of the equations of motion 114 for this discretized mass system:   m 0 0       0 m 0     0 0 m So  x¨1 x¨2 x¨3      2k −k 0  x1           + −k 2k −k  x2  = 0         0 −k 2k x3     m 0 0   2k −k 0          M =  0 m 0  K0 = −k 2k −k          0 0 m 0 −k 2k (5.15) (5.16) Now introducing the excitation matrix K1 , the equations of motion become M¨ x + (K0 + K1 cos Ωt)x = 0 (5.17) where K1 = K0 in this example. Now we can use the approach developed in the previous section to solve for the solutions of this discretized mass system. Substitute the matrices M, K0 and K1 into equation 5.11 and plot the response, then we can get the the plot of the three modes in Figure 5.2. In this example, the value of is set to 0.1. The plot in Figure 5.2 is a certain instant randomly selected from a video of the response. In Figure 5.2, the solid line is the modal response without perturbation, while the dashed line is the one with perturbation. We could see that the mode shapes with perturbation and without perturbation don’t completely coincide with each other, which means the perturbations did affect the modal response of the unperturbed system. 115 mode 2 mode 3 1.5 1 1 1 0.5 0 −0.5 Without Perturbation With Perturbation −1 −1.5 0.5 1 1.5 2 2.5 Mass location 3 Modal Displacement 1.5 Modal Displacement Modal Displacement mode 1 1.5 0.5 0 −0.5 −1 3.5 −1.5 0.5 0.5 0 −0.5 −1 1 1.5 2 2.5 Mass location 3 3.5 −1.5 0.5 1 1.5 2 2.5 Mass location 3 3.5 Figure 5.2 The comparison between the condition without perturbation and with perturbation for three mass system 5.4 Example : Wind Turbine Blade In the previous three mass example, the perturbation with the parametric excitation shows an effect on the modal response of the stationary part of the system. Now we want to apply the perturbation analysis to the wind turbine blades, and see what effect it has. The elements of matrices M, K0 and K1 are obtained from equation 4.7 in Chapter 4, as L mij = m(x)ui uj dx 0 L L k0ij = 0 L + 0 EI(x)ui uj dx − m(x)Ω2 ui uj 0 L dx − mΩ2 ui uj dx Ω2 0 x L L k1ij = 0 m(x)gui uj dx cos Ωt − (5.18) L mwdwui uj dx L g( 0 x mdw)ui uj dx cos Ωt The plots of wind turbine blade perturbed modal responses show a snapshot from the animated motion, as shown in Figure 5.3, Figure 5.4, and Figure 5.5, when the time is 0 and a quarter of the period of vibration. The solid line is the modal response without perturbation, while the dotted line is the one with perturbation. The data from three different blades, the NREL 23-m blade, the NREL 63-m blade, and the Sandia 100-m blade were applied. 116 The results are shown in Figure 5.3, Figure 5.4, and Figure 5.5. From these plots, we could see that the mode shapes with perturbation and without perturbation coincide with each other pretty well, which indicates the perturbed excitation doesn’t significantly affect the modal response. The perturbation part doesn’t show a significant effect on the wind turbine modal response, because the periodic excitation matrix is much smaller compared to the unperturbed modal matrix. Comparing between the three mass system and this wind turbine blade example, it may indicate the perturbation effect depends on the magnitude of the parametric excitation matrix over the orignial modal matrix. 5.4.1 Quantification of the Nonmodal Stiffness Matrix Similar to the “nonmodal ratio” defined in section 4.5.5, here we want to define criteria to quantify the degree to which the parametric excitation stiffness matrix is nonmodal and significant. The undamped modal coordinate transformation gives UT MU = I, UT KU = Λ, UT CU = D, where U is the modal matrix deriving from the eigenvalue problem of the mass and stiffness matrices, K − ω 2 M = 0. For parametric excitation stiffness matrix, we have UT K1 U = P We define a “nonmodal ratio” as the summation of the singular values of matrix P over the summation of the singular values of matrix Λ, where the singular values quantify the amplitude of a matrix. The “nonmodal ratio” of edgewise bending for the three blades are shown in Table 5.2 and those of the flapwise bending are shown in Table 5.3. It shows that with the increase of the size of the blades, the amplitude of the matrix K1 becomes larger compared to the matrix K0 . 117 Table 5.2 “Nonmodal ratio” for edgewise bending of the three blades at rated wind speed Blade “Nonmodal ratio” 23 − m 63 − m 100 − m 3.75e-5 1.27e-4 2.56e-4 Table 5.3 “Nonmodal ratio” for flapwise bending of the three blades at rated wind speed Blade “Nonmodal ratio” 5.5 23 − m 63 − m 100 − m 4.59e-5 1.58e-4 2.55e-4 Conclusion We developed an approach for the perturbation analysis of the parametric excitation in the wind turbine equation of motion. The example of three mass system shows the perturbed effect on the unperturbed system by the excitation terms when the stiffness matrix is of a different size. It is shown for wind turbine blades, the parametric excitation does not have a significant effect on the modal response of the blade. However, if nonlinear terms or resonances are considered, this could be different. On the other hand, this result shows that it is reasonable to neglect the parametric excitation when investigating the mode shapes of the turbine blades with a relative small size. 118 mode 1 mode 2 2 Modal Displacement Modal Displacement 2 1 0 −1 −2 0 5 10 15 20 Axial location (m) mode 3 0 −1 0 5 10 15 20 Axial location (m) −1 0 5 Without Perturbation With Perturbation 2 1 −2 0 −2 25 Modal Displacement Modal Displacement 2 1 25 mode 4 1 0 −1 −2 25 10 15 20 Axial location (m) 0 5 10 15 20 Axial location (m) 25 (a) mode 1 mode 2 2 Modal Displacement Modal Displacement 2 1 0 −1 −2 0 5 10 15 20 Axial location (m) mode 3 −1 10 15 20 Axial location (m) 25 mode 4 2 Modal Displacement Modal Displacement 0 −2 25 0 5 Without Perturbation With Perturbation 2 1 0 −1 −2 1 0 5 10 15 20 Axial location (m) 1 0 −1 −2 25 0 5 10 15 20 Axial location (m) 25 (b) Figure 5.3 The comparison between the conditions without perturbation and with perturbation for 23-m blade (a) when t = 0 (b) when t = quarter period. 119 mode 1 mode 2 2 Modal Displacement Modal Displacement 2 1 0 −1 −2 0 20 40 Axial location (m) 60 mode 3 0 −1 0 20 40 Axial location (m) −1 0 20 40 Without Perturbation Axial location (m) With Perturbation mode 4 2 1 −2 0 −2 Modal Displacement Modal Displacement 2 1 1 0 −1 −2 60 60 0 20 40 Axial location (m) 60 (a) mode 1 mode 2 2 Modal Displacement Modal Displacement 2 1 0 −1 −2 0 20 40 Axial location (m) 60 mode 3 1 0 −1 −2 0 20 40 Axial location (m) 0 −1 −2 0 20 40 Without Perturbation Axial location (m) With Perturbation mode 4 2 Modal Displacement Modal Displacement 2 1 1 0 −1 −2 60 60 0 20 40 Axial location (m) 60 (b) Figure 5.4 The comparison between the conditions without perturbation and with perturbation for 63-m blade (a) when t = 0 (b) when t = quarter period. 120 mode 1 mode 2 2 Modal Displacement Modal Displacement 2 1 0 −1 −2 0 20 40 60 80 Axial location (m) mode 3 1 0 −1 −2 0 20 40 60 80 Axial location (m) 0 −1 −2 100 0 20 40 60 80 Without Perturbation Axial location (m) With Perturbation mode 4 2 Modal Displacement Modal Displacement 2 1 1 0 −1 −2 100 100 0 20 40 60 80 Axial location (m) 100 (a) mode 1 mode 2 2 Modal Displacement Modal Displacement 2 1 0 −1 −2 0 20 40 60 80 Axial location (m) mode 3 1 0 −1 −2 0 20 40 60 80 Axial location (m) 0 −1 −2 100 0 20 40 60 80 Without Perturbation Axial location (m) With Perturbation mode 4 2 Modal Displacement Modal Displacement 2 1 100 100 1 0 −1 −2 0 20 40 60 80 Axial location (m) 100 (b) Figure 5.5 The comparison between the conditions without perturbation and with perturbation for 100-m blade (a) when t = 0 (b) when t = quarter period. 121 Chapter 6 Conclusions and Future Work 6.1 Concluding Remarks We have investigated the modal analysis of non-diagonalizable vibration systems, particularly continuous systems when the damping is nonmodal and there is parametric excitation in the system. 6.1.1 Complex Modal Analysis of a Nonmodally Damped Continuous Beam We have outlined a scheme for the complex modal analysis of non-modally damped distributed parameter systems, and applied this approach to the analysis of an end-damped cantilevered beam. The method involves the use of assumed modes of the undamped system, which are combined based on the state-variable eigensolution of the discretized system. The finite-element method was also utilized to get the mass, stiffness, and damping matrices and further to build and solve a state-variable eigenproblem. The eigenvalues and modal vectors obtained from the assumed-mode method were consistent with those from the finite-element method, as indicated by nearly unit MAC values. These results show that the distributed parameter assumed-mode based method is an efficient way to solve the complex mode problem when non-modal damping is included. The 122 assumed-mode method involves computations of integrals for the low-order modeling of mass, damping and stiffness matrices, and subsequent computations are noniterative and involve lower-order matrices. The assumed-mode method should be equally applicable to structures with nonuniformities, without significant changes in the approach. We applied this method to study features of the end-damped cantilevered beam as a function of the damping coefficient c. With this damping arrangement, most modes are underdamped regardless of c. Each underdamped mode has optimal values of c for generating vibration decay and modal nonsynchronicity. 6.1.2 Experimental Study on Complex Modes of an End Damped Continuous Beam We have conducted an experiment for investigating the complex modal behavior of an enddamped cantilevered beam. Eddy-current induced damping was applied at the end of a cantilevered beam to generate the nonmodal damping. SVMD was applied to extract the complex modes from the free response of the cantilevered beam and to analyze the modal properties including modal frequencies, mode shapes, damping ratio and modal nonsychronicity. The modal properties from the experimental results were compared with those from a numerical analysis of the model. The modal frequencies and mode shapes obtained from the experiments were consistent with those from the model. The variation of damping ratio and traveling index with varying damping coefficient also agreed with the predictions from the model. Over the range of damping coefficients studied in the experiments, we observed a maximum damping ratio in the lowest underdamped mode, which correlated with the maximum modal nonsynchronicity. 123 The results of modal assurance criterion values showed that the SVMD modes were consistent with the modes from the model. As verification, the COD method also produces similar results as SVMD. The mode shapes obtained from COD correlated with the modes from the model. 6.1.3 Modal Analysis of a Wind Turbine Blade We applied the assumed-mode complex modal analysis to study HAWT blades. Based on the equations of motion of a turbine blade, the edgewise and flapwise bending modes of turbine blades were analyzed. Comparing the natural frequencies obtained with the technical report of Bir and Oyague [89], it shows a good agreement. It was shown that when the rotating speed of the hub is slow, i.e. at the speed range of an rotating wind turbine, the mode shapes do not depend significantly on whether the conditions are with or without rotation. The modal frequencies of the blade in the vertical and horizontal positions have minor differences. This provides the theoretical foundation for the research work by building a small model of the turbine blade in the lab. However, the difference between the modal frequencies at horizontal and vertical positions becomes larger when the size of the turbine blade increases and it is verified on the 63-m and 100-m wind turbine blades. The effects of rotational position due to gravity of the three different blades, the NREL 23-m blade, the NREL 63-m blade and the Sandia 100-m blade, were compared. With the increase of the size of the blade, the difference of the modal frequencies at different rotation positions (horizontal position and vertical position) increases. The aerodynamic loading in both flapwise and edgewise directions were derived and the damping terms were derived for analysis. We isolated the effect of nonmodal damping in the equation of motion for a rotating wind turbine blade and further solved for the modes in 124 flapwise and edgewise bending. The data from the NREL and Sandia technical reports were applied in the computation and the samples of 23-m, 63-m and 100-m wind turbine blades were analyzed. From the plot of the first four modes and the calculation of damping ratios, it shows that the aerodynamic loading has a much stronger effect on the flapwise bending compared to the edgewise bending. With the increase of the length of wind turbine blades, the damping ratio and nonsynchronicity increased. The aerodynamic loading exerts a more significant effect on large wind turbine blades. The damping ratio increases with the increase of the wind speed. However, the damping ratios vary greatly for 23-m blade, and vary little for 100-m blade, which means the effect of the wind speed on the nonmodal damping of the wind turbine blades is different for different blades. The aerodynamic loading modeling of the wind turbine blade in this paper used simple quasisteady aeroelastic model, and it could be improved by applying a dynamic model, considering the unsteady aerodynamic effects such as stall hysterisis [100, 133, 134]. The rotational effects were not considered in this Chapter, and could be further added into modeling in the future. 6.1.4 Modal Analysis of Parametrically Excited Systems We developed an approach for the perturbation analysis of the parametric excitation in the wind turbine equation of motion. The example of three mass system shows the perturbed effect on the unperturbed system by the excitation terms when the stiffness matrix is of a different size. It is shown for wind turbine blades, the parametric excitation does not have a significant effect on the modal response of the blade. However, if nonlinear terms or resonances are considered, this could be different. On the other hand, this result shows that 125 it is reasonable to neglect the parametric excitation when investigating the mode shapes of the turbine blades with a relative small size. 126 6.2 Future Work Future work will include the following: ❼ The assumed-mode approach for complex modal analysis can be applied to more severely nonmodally damped systems, such as the flutter of a flag. It also may be worthwhile to consider the large hydropower turbine blade. In this case, a much larger density of water will bring a more significant nonmodal damping, and the associated modal properties could be investigated. ❼ The experiment in Chapter 3 did not capture the first mode. However, if we replace the accelerometers with strain gages and redo the experiment, the first mode should be obtained. The strain gages have a better ability of measuring low frequency signals. Then we can have a better knowledge of the real mode pair of the first mode. ❼ We can consider replacing the eddy-current damper with other kinds of concentrated dampers, such as VIBEX, which is a damping product produced by the Permawick Corporation. It would be worthy to experimentally quantify the damping properties and complex modes, and to see how they change when the damping material changes in the real experiment. ❼ The aerodynamic loading modeling of the wind turbine blade in this paper used simple quasi-steady aeroelastic model, and it could be improved by applying a dynamic model, such as the unsteady modeling of ONERA [133, 134]. The rotational effects were not considered in the nonmodal damping, and could be further added into the model in the future. ❼ The perturbation analysis of parametrically excited MDOF systems presented here 127 was explortory, and only considered nonresonant conditions. It would be important to futher interpret and improve the perturbation analysis of the parametric excitation in wind turbine application. A comparison between the numerical modal solution and the perturbation modal solution would be useful. Future work could more closely tie the concepts to nonlinear normal mode theory, or analyze resonant cases, e.g. with the multiple scale method. 128 APPENDICES 129 Appendix A Study of the Effect of the Number of Samples and Filter Frequency on SVMD The number of samples applied in SVMD can have a significant effect on the extracted modes, as in reference [42]. The filter frequency set in the high-pass filter could also have an effect on the extracted modes. Therefore, here we present a more detailed study of the dependence of the results on the number of samples N and the high-pass filter frequency. Number of Samples As mentioned in Chapter 3, a large data sampling could contribute a better mode for lower modes and a relatively small data sampling may be better for higher modes, considering higher modes decay faster. However, state-variable modal coordinates need to be examined because even if for higher modes, in the latter part of the modal coordinates, small oscillations may still contain modal information. Therefore, a detailed study of the dependency of number of samples is necessary. The damping ratio, traveling index, and modal frequencies obtained when the number of samples N = 12500 are plotted against the results from the model and shown in Figure A.1. 130 The lines represent the values from the model, the triangles represent the values of mode 2, the circles represent the values of mode 3, and the asterisks represent the values of mode 4. The mode shapes obtained when N = 12500 are shown in Figure A.2. The plots of the modal properties by using the number of samples N = 6000 are shown in Figure A.3. The mode shapes obtained by using N = 3000 are shown in Figure A.4. The plot of N = 3000 are shown in Figure A.5. The mode shapes obtained when N = 3000 are shown in Figure A.6. The plot of N = 2000 are shown in Figure A.7. The mode shapes obtained when N = 2000 are shown in Figure A.8. N=12500 N=12500 0.12 Damping ratio 0.1 0.08 0.06 N=12500 0.35 1400 0.3 1200 0.25 1000 imaginary part mode 2 (model) mode 3 (model) mode 4 (model) mode 2 (exp) mode 3 (exp) mode 4 (exp) Traveling index 0.14 0.2 0.15 800 600 0.04 0.1 400 0.02 0.05 200 0 0 10 20 30 40 Damping coefficient 50 0 0 10 20 30 40 Damping coefficient 50 0 −120 −100 −80 −60 −40 real part −20 0 Figure A.1 The comparison of damping ratio, traveling index, and modal frequencies based on N = 12500 sample points The state-variable modal coordinates for each different data sampling case when damping coefficient is 10 kg/s are shown in Figure A.9, Figure A.10, Figure A.11 and Figure A.12. The mode shapes obtained from N = 12500 were better compared to those from N = 2000 for mode 3 and mode 4. Compare the modal coordinates of N = 12500 and N = 2000, we can see in the latter half of N = 12500, there was still some oscillation for mode 3 and mode 4, which could contain some mode information. Therefore, in this experiment, N = 12500 produced a better mode. 131 mode 3 0.2 0.2 0.2 0.1 0 −0.1 −0.2 0 −0.2 0.2 0.4 location (m) 0 −0.1 −0.2 0.6 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.1 0 −0.1 −0.2 −0.1 Real & Imaginary part Real & Imaginary part 0.2 0 0.1 −0.2 −0.1 0 0.1 0.2 Experiments Real part Model mode 3 0.4 0.4 −0.4 Imaginary part 0.3 −0.2 −0.1 0 0.1 0.2 Real part mode 2 Real & Imaginary part mode 4 0.3 Imaginary part Imaginary part mode 2 0.3 0.6 −0.05 0 0.05 Real part mode 4 0.1 0.4 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 Figure A.2 The comparison of the real and imaginary parts of the mode shapes based on N = 12500 sample points when c = 10 kg/s Modal Identification with Different Filter Frequencies As mentioned earlier in Chapter 3, to remove a low frequency “drift” in the integrated signal ensembles, data were high-pass filtered. A higher filter frequency can remove the drift better, but brings in more similar time constant, perhaps distortion. A previous study in reference [42] showed that a maximum of one half of the fundamental frequency by FFT should be employed as the filter frequency. The first fundamental frequency cannot be captured in this experiment due to the low limit of the accelerometers and real mode pair of the first mode. We filtered the signal with different filter frequencies, and the results are presented below. 132 N=6000 N=6000 mode 2 (model) mode 3 (model) mode 4 (model) mode 2 (exp) mode 3 (exp) mode 4 (exp) 0.12 Traveling index Damping ratio 0.1 0.08 0.06 N=6000 0.35 1400 0.3 1200 0.25 1000 imaginary part 0.14 0.2 0.15 800 600 0.04 0.1 400 0.02 0.05 200 0 0 10 20 30 40 Damping coefficient 0 50 0 10 20 30 40 Damping coefficient 50 0 −120 −100 −80 −60 −40 real part −20 0 Figure A.3 The comparison of damping ratio, traveling index, and modal frequencies based on N = 6000 sample points mode 3 0.2 0.2 0.2 0.1 0 −0.1 0.1 0 −0.1 −0.2 0 Real part mode 2 −0.2 Experiments Model 0.4 Real & Imaginary part 0.4 0.2 0.2 0 −0.2 0 0.2 0.4 location (m) 0.6 0 Real part 0.2 0 −0.1 −0.1 mode 3 0.1 mode 4 0 −0.2 0 0 Real part 0.4 0.2 −0.4 0.1 −0.2 Real & Imaginary part −0.2 −0.4 Imaginary part 0.3 −0.2 Real & Imaginary part mode 4 0.3 Imaginary part Imaginary part mode 2 0.3 0.2 0.4 location (m) 0.6 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 Figure A.4 The comparison of the real and imaginary parts of the mode shapes based on N = 6000 sample points when c = 10 kg/s 133 N=3000 N=3000 mode 2 (model) mode 3 (model) mode 4 (model) mode 2 (exp) mode 3 (exp) mode 4 (exp) 0.12 Traveling index Damping ratio 0.1 0.08 0.06 N=3000 0.35 1400 0.3 1200 0.25 1000 imaginary part 0.14 0.2 0.15 800 600 0.04 0.1 400 0.02 0.05 200 0 0 10 20 30 40 Damping coefficient 0 50 0 10 20 30 40 Damping coefficient 0 −120 50 −100 −80 −60 −40 real part −20 0 Figure A.5 The comparison of damping ratio, traveling index, and modal frequencies based on N = 3000 sample points mode 3 0.2 0.2 0.2 0.1 0 −0.1 0.1 0 −0.1 −0.2 0 0.2 Real part −0.2 0 0.2 Experiments Real part Model mode 3 0.4 mode 2 Real & Imaginary part 0.4 0.2 0 −0.2 0 0.2 0.4 location (m) 0.6 0 −0.1 −0.1 −0.2 0.2 0.4 location (m) 0.1 mode 4 0 0 0 Real part 0.4 0.2 −0.4 0.1 −0.2 Real & Imaginary part −0.2 −0.4 Imaginary part 0.3 −0.2 Real & Imaginary part mode 4 0.3 Imaginary part Imaginary part mode 2 0.3 0.6 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 Figure A.6 The comparison of the real and imaginary parts of the mode shapes based on N = 3000 sample points when c = 10 kg/s 134 N=2000 N=2000 mode 2 (model) mode 3 (model) mode 4 (model) mode 2 (exp) mode 3 (exp) mode 4 (exp) 0.12 Traveling index Damping ratio 0.1 0.08 0.06 N=2000 0.35 1400 0.3 1200 0.25 1000 imaginary part 0.14 0.2 0.15 800 600 0.04 0.1 400 0.02 0.05 200 0 0 10 20 30 40 Damping coefficient 0 50 0 10 20 30 40 Damping coefficient 0 −120 50 −100 −80 −60 −40 real part −20 0 Figure A.7 The comparison of damping ratio, traveling index, and modal frequencies based on N = 2000 sample points mode 3 0.2 0.2 0.2 0.1 0 −0.1 0.1 0 −0.1 −0.2 0 0.2 Real part −0.2 0 0.2 Experiments Real part Model mode 3 0.4 mode 2 Real & Imaginary part 0.4 0.2 0 −0.2 0 0.2 0.4 location (m) 0.6 0 −0.1 −0.1 −0.2 0.2 0.4 location (m) 0.1 mode 4 0 0 0 Real part 0.4 0.2 −0.4 0.1 −0.2 Real & Imaginary part −0.2 −0.4 Imaginary part 0.3 −0.2 Real & Imaginary part mode 4 0.3 Imaginary part Imaginary part mode 2 0.3 0.6 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 Figure A.8 The comparison of the real and imaginary parts of the mode shapes based on N = 2000 sample points when c = 10 kg/s 135 −6 −6 x 10 6 4 2 SVMC 3 SVMC 2 4 2 0 0 −2 −2 −4 x 10 0 500 1000 1500 Time Histoy −4 2000 −6 500 1000 1500 Time Histoy 2000 500 1000 1500 Time Histoy 2000 −7 x 10 1 0 5 x 10 SVMC 5 SVMC 4 0.5 0 0 −0.5 −1 0 500 1000 1500 Time Histoy −5 2000 0 Figure A.9 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 2000 sample points −6 5 −6 x 10 4 x 10 SVMC 3 SVMC 2 2 0 −5 0 −2 −10 0 1000 2000 Time Histoy −4 3000 −6 1 1.5 3000 1000 2000 Time Histoy 3000 x 10 1 SVMC 5 SVMC 4 1000 2000 Time Histoy −6 x 10 0.5 0 −0.5 −1 0 0.5 0 −0.5 0 1000 2000 Time Histoy 3000 −1 0 Figure A.10 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 3000 sample points 136 −5 1 2 0.5 1 0 −0.5 −1 −1.5 x 10 3 SVMC 3 SVMC 2 −5 x 10 1.5 0 −1 −2 0 2000 4000 Time Histoy −3 6000 0 −5 6000 4 SVMC 5 SVMC 4 2000 4000 Time Histoy x 10 6 0.5 0 −0.5 −1 6000 −6 x 10 1 2000 4000 Time Histoy 2 0 −2 0 2000 4000 Time Histoy −4 6000 0 Figure A.11 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 6000 sample points −4 1 −5 x 10 4 2 SVMC 3 SVMC 2 0.5 x 10 0 −0.5 0 −2 −1 0 2000 4000 6000 −4 0 8000 10000 12000 2000 Time Histoy x 10 2 8000 10000 12000 x 10 1 0 −5 −10 0 6000 −5 SVMC 5 SVMC 4 5 4000 Time Histoy −5 0 −1 2000 4000 6000 −2 0 8000 10000 12000 Time Histoy 2000 4000 6000 8000 10000 12000 Time Histoy Figure A.12 SVMD modal coordinates are shown for mode 2 to mode 5 based on N = 12500 sample points 137 mode 3 0.2 0.2 0.2 0.1 0 −0.1 0.1 0 −0.1 −0.2 0 0.2 Real part −0.2 0 0.2 Experiments Real part Model mode 3 0.4 mode 2 Real & Imaginary part 0.4 0.2 0 −0.2 0 0.2 0.4 location (m) 0 −0.1 −0.1 −0.2 0.2 0.4 location (m) 0.1 mode 4 0 0 0 Real part 0.4 0.2 −0.4 0.6 0.1 −0.2 Real & Imaginary part −0.2 −0.4 Imaginary part 0.3 −0.2 Real & Imaginary part mode 4 0.3 Imaginary part Imaginary part mode 2 0.3 0.2 0 −0.2 −0.4 0.6 0 0.2 0.4 location (m) 0.6 Figure A.13 The comparison of the real and imaginary parts of the mode shapes based on N = 2000 sample points with filter frequency at 5 Hz when c = 10 kg/s mode 3 0.2 0.2 0.2 0.1 0 −0.1 0.1 0 −0.1 −0.2 0 0.2 Real part −0.2 0 0.2 Experiments Real part Model mode 3 0.4 mode 2 Real & Imaginary part 0.4 0.2 0 −0.2 0 0.2 0.4 location (m) 0.6 0 −0.1 −0.1 −0.2 0.2 0.4 location (m) 0.1 mode 4 0 0 0 Real part 0.4 0.2 −0.4 0.1 −0.2 Real & Imaginary part −0.2 −0.4 Imaginary part 0.3 −0.2 Real & Imaginary part mode 4 0.3 Imaginary part Imaginary part mode 2 0.3 0.6 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 Figure A.14 The comparison of the real and imaginary parts of the mode shapes based on N = 2000 sample points with filter frequency at 10 Hz when c = 10 kg/s 138 mode 4 0.2 0.2 0.2 0.1 0 −0.1 −0.2 Real & Imaginary part 0.2 0 −0.2 0 0.2 0.4 location (m) 0 −0.1 −0.2 −0.2 0 0.2 Experiments Real part Model mode 3 0.4 0.4 −0.4 0.1 0.2 0 −0.2 −0.4 0.6 0 0.2 0.4 location (m) 0.1 0 −0.1 −0.2 −0.1 Real & Imaginary part 0 0.2 Real part mode 2 Imaginary part 0.3 −0.2 Real & Imaginary part mode 3 0.3 Imaginary part Imaginary part mode 2 0.3 0.1 0.2 0.4 location (m) 0.6 0.4 0.2 0 −0.2 −0.4 0.6 0 Real part mode 4 0 Figure A.15 The comparison of the real and imaginary parts of the mode shapes based on N = 3000 sample points with filter frequency at 5 Hz when c = 10 kg/s mode 3 0.2 0.2 0.2 0.1 0 −0.1 0.1 0 −0.1 −0.2 0 0.2 Real part −0.2 0 0.2 Experiments Real part Model mode 3 0.4 mode 2 Real & Imaginary part 0.4 0.2 0 −0.2 0 0.2 0.4 location (m) 0.6 0 −0.1 −0.1 −0.2 0.2 0.4 location (m) 0.1 mode 4 0 0 0 Real part 0.4 0.2 −0.4 0.1 −0.2 Real & Imaginary part −0.2 −0.4 Imaginary part 0.3 −0.2 Real & Imaginary part mode 4 0.3 Imaginary part Imaginary part mode 2 0.3 0.6 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 Figure A.16 The comparison of the real and imaginary parts of the mode shapes based on N = 3000 sample points with filter frequency at 10 Hz when c = 10 kg/s 139 mode 3 0.2 0.2 0.2 0.1 0 −0.1 0.1 0 −0.1 −0.2 −0.2 0 Real part 0.2 mode 2 0 Real part 0.2 −0.2 0 0.2 0.4 location (m) 0 Real part 0.1 mode 4 0.4 0.2 0 −0.2 −0.4 0.6 −0.1 mode 3 Real & Imaginary part Real & Imaginary part 0 0 −0.1 0.4 0.2 0.1 −0.2 −0.2 Experiments Model 0.4 −0.4 Imaginary part 0.3 −0.2 Real & Imaginary part mode 4 0.3 Imaginary part Imaginary part mode 2 0.3 0 0.2 0.4 location (m) 0.6 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 Figure A.17 The comparison of the real and imaginary parts of the mode shapes based on N = 12500 sample points with filter frequency at 5 Hz when c = 10 kg/s mode 3 0.2 0.2 0.2 0.1 0 −0.1 0.1 0 −0.1 −0.2 0 0.2 Real part −0.2 0 0.2 Experiments Real part Model mode 3 0.4 mode 2 Real & Imaginary part 0.4 0.2 0 −0.2 0 0.2 0.4 location (m) 0.6 0 −0.1 −0.1 −0.2 0.2 0.4 location (m) 0.1 mode 4 0 0 0 Real part 0.4 0.2 −0.4 0.1 −0.2 Real & Imaginary part −0.2 −0.4 Imaginary part 0.3 −0.2 Real & Imaginary part mode 4 0.3 Imaginary part Imaginary part mode 2 0.3 0.6 0.2 0 −0.2 −0.4 0 0.2 0.4 location (m) 0.6 Figure A.18 The comparison of the real and imaginary parts of the mode shapes based on N = 12500 sample points with filter frequency at 10 Hz when c = 10 kg/s 140 Table A.1 MAC values comparing the modes from the experimental results and from the model with varying number of samples and varying filter frequency (ff) when c = 10 kg/s Mode Number mode 2 mode 3 mode 4 ff type ff=2 Hz ff=5 Hz ff=10 Hz ff=2 Hz ff=5 Hz ff=10 Hz ff=2 Hz ff=5 Hz ff=10 Hz N = 12500 N = 6000 N = 3000 N = 2000 0.9985 0.9149 0.9800 0.8737 0.8908 0.9675 0.8944 0.8897 0.9146 0.9124 0.9077 0.8702 0.9767 0.9765 0.9729 0.9666 0.9765 0.9791 0.9665 0.9662 0.9680 0.9790 0.9675 0.9664 0.9318 0.9309 0.9231 0.9240 0.7197 0.9297 0.9285 0.9271 0.9419 0.9341 0.9230 0.9296 A calculation of the MAC values for different filter frequency are listed in Table A.1. 141 Appendix B Comparison for the axial stiffness and bending stiffness in a rotating beam In reference [64], Yoo and Shin derived equations of motion (EOM) of a rotating beam, including the EOM in axial direction and bending direction. Here we did a rough comparison between the first modal frequency of the axial direction and the bending direction. Suppose the assumed modal functions for the bending and axial displacement are the same in this rough estimation, as φi = 1 − cos γx (B.1) Here we neglect Ω2 terms, as we see in Chapter 4, Ω has little effect on modal frequencies of wind-turbine blades. The first element of the stiffness term in the axial direction is obtained in a simplifying form as a = EA K11 L 0 φ1,x φ1,x dx = EA( 2π 2 L 2 sin (γx)dx ) L 0 (B.2) where E is Young’s modulus, A is the cross-sectional area, φ1,x is the first order derivative of the assumed modal function with respect to x. The first element of the stiffness term in transverse direction is obtained in a simplifying form as 142 b = EI K11 L 0 L φ1,xx φ1,xx dx = EI γ 4 cos2 (γx)dx (B.3) 0 where I is the area of moment of inertia of a cross-section of the beam. Then after intergration, we obtain EAπ 2 8L 4 EIπ b = K11 32L3 a = K11 Kb (B.4) 2 2 1 h π The ratio between these two stiffness terms are K11 a = 48 ( 2 ) L 11 Since in a wind turbine blade, h is much less than L, the ratio between the bending stiffness and the axial stiffness is very small. Therefore, the axial stiffness is very large, and we can assume the beam to be inextensible when we derive the equation of motion in bending. 143 BIBLIOGRAPHY 144 BIBLIOGRAPHY [1] T. K. Caughey and M. O’Kelly. Classical normal modes in damped linear systems. Journal of Applied Mechanics, 32(3):583–588, 1965. [2] F. E. Udwadia. A note on non-proportional damping. Journal of Engineering Mechanics, 135(11):1248–1256, 2009. [3] M. I. Friswell and U. Prells. Non-proportional damping and defective systems. In International Conference on Noise and Vibration Engineering, Katholieke University Leuven, Belgium, September 2000. [4] N. Blaise, T. Canor, and V. Denoel. Equivalent static wind loads for structures with non-proportional damping. In 11th International Conference RASD 2013, Pisa, July 2013. [5] U. Prells and M. I. Friswell. A measure of non-proportional damping. Mechanical Systems and Signal Processing, 14(2):125–137, 2000. [6] S. Adhikari. Damping modelling using generalized proportional damping. Journal of Sound and Vibration, 293:156–170, 2006. [7] J. Angeles and S. Ostrovskaya. The proportional-damping matrix of arbitrarily damped linear mechanical systems. Journal of Applied Mechanics, 69:649–656, 2002. [8] K. R. Chung and C. W. Lee. Dynamic reanalysis of weakly non-proportionally damped systems. Journal of Sound and Vibration, 111(1):37–50, 1986. [9] B. F. Feeny. A complex orthogonal decomposition for wave motion analysis. Journal of Sound and Vibration, 310:77–90, 2008. [10] G. Prater and S. Singh. Quantification of the extent of non-proportional viscous damping in discrete vibratory systems. Journal of Sound and Vibration, 104(1):109–125, 1986. [11] R. Singh and G. Prater. Complex eigensolution for longitudinally vibrating bars with a viscously damped boundary. Journal of Sound and Vibration, 133(2):364–367, 1989. 145 [12] G. Prater and R. Singh. Eigenproblem formulation, solution and interpretation for non-proportionally damped continuous beams. Journal of Sound and Vibration, 143(1):125–142, 1990. [13] Andrew J. Hull. A closed form solution of a longitudinal beam with a viscous boundary condition. Technical report, Naval Undersea Warfare Center Detachment, 1992. [14] A. J. Hull. A closed form solution of a longitudinal bar with a viscous boundary condition. Journal of Sound and Vibration, 169:19–28, 1994. [15] G. Oliveto and A. Santini. Complex modal analysis of a continuous model with radiation damping. Journal of Sound and Vibration, 192:15–33, 1996. [16] G. Oliveto and A. Santini. Time domain response of a one-dimensional soil-structure interacting model via complex analysis. Engineering Structures, 18:425–436, 1996. [17] G. Oliveto, A. Santini, and E. Tripodi. Complex modal analysis of a flexural vibrating beam with viscous end conditions. Journal of Sound and Vibration, 200(3):327–345, 1997. [18] S. Krenk. Vibrations of a taut cable with an external damper. Journal of Applied Mechanics, 67:772–776, 2000. [19] S. R. Nielsen and S. Krenk. Whirling motion of a shallow cable with viscous dampers. Journal of Sound and Vibration, 265:417–435, 2003. [20] Steen Krenk. Complex modes and frenquencies in damped structural vibrations. Journal of Sound and Vibration, 270(4-5):981–996, 2004. [21] L. Sirota and Y. Halevi. Modal representation of second order flexible structures with damped boundaries. Journal of Vibration and Acoustics, 135:064508(1–6), 2013. [22] J. Zarek and B. Gibbs. The derivation of eigenvalues and mode shapes for the bending motion of a damped beam with general end conditions. Journal of Sound and Vibration, 78:185–196, 1985. [23] B. Yang and X. Wu. Transient response of one-dimensional distributed systems: a closed form eigenfuction expansion realization. Journal of Sound and Vibration, 208:763–776, 1997. 146 [24] J. A. Main and N. P. Jones. Free vibrations of taut cable with attached damper. Journal of Engineering Mechanics, 128:1062–1071, 2002. [25] J. Juang and R. Pappa. An eigensystem realization algorithm for modal parameter identification and model reduction. Journal of Guidance, 8 (5):620–627, 1985. [26] S. Ibrahim and E. Mikulcik. A method for the direct identification of vibration parameters from the free response. Shock and Vibration Bulletin, 47(4):183–198, 1977. [27] G. Kerschen, F. Poncelet, and J. Golinval. Physical interpretation of independent component analysis in structural dynamics. Mechanical Systems and Signal Processing, 21(4):1561–1575, 2007. [28] F. Poncelet, G. Kerschen, J. Golinval, and D. Verhelst. Output-only modal analysis using blind source separation techniques. Mechanical Systems and Signal Processing, 21:2335–2358, 2007. [29] D. L. Brown, R. J. Allemang, R. D. Zimmerman, and M. Mergeay. Parameter estimation techniques for modal analysis. In SAE Technical Paper 790221, 1979. [30] S. Braun and Y. Ram. Determination of structural modes via the Prony model : system order and noise induced poles. Journal of the Acoustical Society of America, 81(5):1447–1459, 1987. [31] J. Arruda, P. Campos, and J. Golinval. Experimental determination of flexural power flow in beams using a modified Prony method. Journal of Sound and Vibration, 197(3):309–328, 1996. [32] P. V. Overschee and B. De Moor. Subspace Identification for Linear Systems : TheoryImplementation-Applications. Kluwer Academic Publishers, 1996. [33] M. Richardson and D. Formenti. Parameter estimation from frequency response measurements using rational fraction polynomials. In Proceedings of the International Modal Analysis Conference, 1982. [34] C. Shih, Y. Tsuei, R. Allemang, and D. Brown. Complex mode indication function and its application to spatial domain parameter estimation. Mechanical Systems and Signal Processing, 2:367–377, 1988. [35] L. Zhang R. Brincker and P. Andersen. Modal identification of output-only systems using frequency domain decomposition. Smart Materials and Structures, 10:441–445, 2001. 147 [36] D. Chelidze and W. Zhou. Smooth orthogonal decomposition-based vibration mode identification. Journal of Sound and Vibration, 292:461–473, 2006. [37] W. Zhou and D. Chelidze. Generalized eigenvalue decomposition in time domain modal parameter identification. Journal of Vibration and Acoustics, 130(1):011001, 2008. [38] U. Farooq and B. F. Feeny. Smooth orthogonal decomposition for randomly excited systems. Journal of Sound and Vibration, 316(1-5):137–146, 2008. [39] B. F. Feeny and U. Farooq. A nonsymmetric state-variable decomposition for modal analysis. Journal of Sound and Vibration, 310(4-5):792–800, 2008. [40] B. F. Feeny and U. Farooq. A state-variable decomposition method for estimating modal parameters. In ASME International Design Engineering Technical Conferences, Las Vegas, 2007. [41] U. Farooq and B. F. Feeny. An experimental investigation of state-variable modal decomposition for modal analysis. Journal of Vibration and Acoustics, 134:021017(1– 8), 2012. [42] Umar Farooq. Orthogonal Decomposition Methods for Modal Analysis. PhD thesis, Michigan State University, 2009. [43] S. R. Ibrahim. Computation of normal modes from identified complex modes. AIAA Journal, 21:446–451, 1983. [44] N. Maria and D. J. Ewins. A new approach for the modal identification of lightly damped structures. Mechanical Systems and Signal Processing, 3(2):173–193, 1989. [45] B. F. Feeny. Complex modal decomposition for estimating wave properties in onedimensional media. Journal of Vibration and Acoustics, 135(1-2), 2013. [46] H. A. Sodano, J. Bae, D. J. Inman, and W. K. Belvin. Concept and model of eddy current damper for vibration suppression of a beam. Journal of Sound and Vibration, 288:1177–1196, 2005. [47] H. A. Sodano. Development of Novel Eddy Current Dampers for the Suppression of Structural Vibrations. PhD thesis, Virginia Polytechnic Institute and State University, 2005. 148 [48] R. P. Feynman, R. B. Leighton, and M. Sands. The Feynman Lectures on Physics. Addison-Wesley Publishing Company, 1977. [49] L. Cadwell. Magnetic damping: analysis of an eddy current brake using an air track. Journal of Physics, 64:917–923, 1996. [50] L. Davis and J. Reitz. Eddy currents in finite conducting sheets. Journal of Applied Physics, 42:4119–4127, 1971. [51] M. Heald. Magnetic braking improved theory. American Journal of Physics, 56:521– 522, 1988. [52] K. Lee and K. Park. Optimal robust control of a contact-less brake system using an eddy current. Mechatronics, 9:615–631, 1999. [53] K. Nagaya, H. Kojima, Y. Karube, and Kibayashi. Braking force and damping coefficient of eddy current brakes consisting of cylindrical magnets and plate conductors of arbitrary shape. IEEE Transactions on Magnetics, MAG-20:2136–2145, 1984. [54] D. Schieber. Optimal dimensions of rectangular electromagnet for braking purposes. IEEE Transactions on Magnetics, 11(3):948–952, 1975. [55] E. Simeu and D. Georges. Modeling and control of an eddy current brake. Control Engineering Practice, 4:19–26, 1996. [56] H. Wiederick, N. Gauthier, D. Campbell, and P. Rochon. Magnetic braking: Simple theory and experiment. American Journal of Physics, 55:500–503, 1987. [57] R. Fredrick and M. Darlow. Operation of an electromagnetic eddy current damper with a supercritical shaft. ASME Journal of Vibration and Acoustics, 116:578–580, 1994. [58] R. Fung, J. Sun, and S. Hsu. Vibration control of the rotating flexible-shaft/multiflexible-disk system with the eddy-current damper. Journal of Vibrations and Acoustics, 124:519–526, 2002. [59] Y. Klingerman, O. Gottlieb, and M. Darlow. Analytic and experimental evaluation of instability in rotordynamic system with electromagnetic eddy-current damper. Journal of Vibration and Acoustics, 120:272–278, 1998. 149 [60] M. Karnopp. Permanent magnet linear motors used as variable mechanical damper for vehicle suspensions. Vehicle System Dynamics, 18:187–200, 1989. [61] G. Larose, A. Larsen, and E. Svensson. Modeling of tuned mass dampers for wind tunnel tests on a full-bridge aeroelastic model. Journal of Wind Engineering and Industrial Aerodynamics, 54:427–437, 1995. [62] Y. Matsuzaki, Y. Ishikubo, T. Kamita, and Ikeda. Vibration control system using electromagnetic forces. Journal of Intelligent Material Systems and Structures, 8:751– 756, 1997. [63] J. Bae, M. Kwak, and D. Inman. Vibration suppression of a cantilever beam using eddy current damper. Journal of Sound and Vibration, 284:805824, 2005. [64] H. H. Yoo and S. H. Shin. Vibration analysis of rotating cantilever beams. Journal of Sound and Vibration, 212:807–828, 1998. [65] H. H. Yoo, J. E. Cho, and J. T. Chung. Modal analysis and shape optimization of rotating cantilever beams. Journal of Sound and Vibration, 290(1):223–241, 2006. [66] H. S. Lim, J. Chung, and H. H. Yoo. Modal analysis of a rotating multi-packet blade system. Journal of Sound and Vibration, 325(3):513–531, 2009. [67] H. Kim, H. H. Yoo, and J. Chung. Dynamic model for free vibration and response analysis of rotating beams. Journal of Sound and Vibration, 332(22):5917–5928, 2013. [68] G. Surace, V. Anghel, and C. Mares. Coupled bending-bending-torsion vibration analysis of rotating pretwisted blades: an integral formulation and numerical examples. Journal of Sound and Vibration, 206(4):473–486, 1997. [69] T. Yokoyama. Free vibration characteristics of rotating Timoshenko beams. International Journal of Mechanical Sciences, 30:743–755, 1988. [70] S. Naguleswaran. Lateral vibration of a centrifugally tensioned Euler-Bernoulli beam. Journal of Sound and Vibration, 176:613–624, 1994. [71] S. Lee and S. Lin. Bending vibrations of rotating nonuniform Timoshenko beams with an elastically restrained root. Journal of Applied Mechanics, 61:949–955, 1994. [72] X. Xing and B. F. Feeny. Complex modal analysis of a non-modally damped continuous beam. In proceedings of the ASME IDETC, Portland, OR, August 2013. 150 [73] C. Cooley and R. Parker. Vibration of spinning cantilever beams undergoing coupled bending and torsional motion. In In proceedings of the ASME IDETC, Portland, OR, August 2013. [74] C. G. Cooley. High-Speed Dynamics and Vibration of Planetary Gears, Vibration of Spinning Cantilevered Beams, and an Efficient Computational Method for Gear Dynamics. PhD thesis, Ohio State University, 2012. [75] K. V. Avramov, C. Pierre, and N. Shyriaieva. Flexural-flexural-torsional nonlinear vibrations of pre-twisted rotating beams with asymmetric cross-sections. Journal of Vibration and Control, 13(4):329–364, 2007. [76] M. R. Silva. Non-linear flexural-flexural-torsional-extensional dynamics of beams-1. formulation. Int. J. Solids Structures, 24(12):1225–1234, 1988. [77] D. H. Hodges and E. H. Dowell. Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades. Technical report, NASA, 1974. [78] J. Wendell. Aeroelastic stability of wind turbine rotor blades. Technical report, U.S. Department of Energy, 1978. [79] J. Wendell. Simplified aeroelastic modeling of horizontal-axis wind turbines. Technical report, Massachusetts Inst. of Tech., Cambridge, USA, 1982. [80] B. S. Kallesoe. Equations of motion for a rotor blade, including gravity, pitch action and rotor speed variations. Wind Energy, 10:209–230, 2007. [81] J. W. Larsen and S. R. Nielsen. Non-linear dynamics of wind turbine wings. International Journal of Nonlinear Mechanics, 41:629–643, 2006. [82] J. W. Larsen and S. R. Nielsen. Non-linear parametric instability of wind turbine wings. Journal of Sound and Vibration, 299:64–82, 2007. [83] O. Turhan and G. Bulut. Dynamic stability of rotating blades (beams) eccentrically clamped to a shaft with fluctuating speed. Journal of Sound and Vibration, 280(35):945–964, 2005. [84] O. Turhan and G. Bulut. On nonlinear vibrations of a rotating beam. Journal of Sound and Vibration, 322:314–335, 2009. 151 [85] D. Caruntu. On nonlinear forced response of nonuniform blades. In Proceeding of the 2008 ASME Dynamic Systems and Control Conference, number DSCC2008-2157, 2008. [86] T. Inoue, Y. Ishida, and T. Kiyohara. Nonlinear vibration analysis of the wind turbine blade. Journal of Vibration and Acoustics, 134(3), 2012. [87] V. Ramakrishnan and B. F. Feeny. Resonances of a forced mathieu equation with reference to wind turbine blades. Journal of Vibration and Acoustics, 134:064501 (5 pages), 2012. [88] D. T. Griffith and T. D. Ashwill. The Sandia 100-meter all-glass baseline wind turbine blade: Snl 100-00. Technical report, Sandia National Laboratories, 2011. [89] G. S. Bir and F. Oyague. Estimation of blade and tower properties for the gearbox research collaborative wind turbine. Technical report, National Renewable Energy Laboratory, U.S. Department of Energy, 2007. [90] J. M. Jonkman. Dynamics modeling and loads analysis of an offshore floating wind turbine. Technical report, National Renewable Energy Laboratory, 2007. [91] M. S. Allen. A global, single-inputmulti-output (simo) implementation of the algorithm of mode isolation and application to analytical and experimental data. Mechanical Systems and Signal Processing, 20:1090–1111, 2006. [92] M. S. Allen and J. H. Ginsberg. Floquet modal analysis to detect cracks in a rotating shaft on anisotropic supports. In 24th International Modal Analysis Conference (IMAC XXIV), St.Louis, MO, 2006. [93] M. S. Allen. Floquet experimental modal analysis for systems identification of linear time-periodic systems. In ASME 2007 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, Las Vegas, Nevada, USA, 2007. [94] M. S. Allen. Frequency-domain identification of linear time-periodic systems using LTI techniques. Journal of Computational and Nonlinear Dynamics, 4, 2009. [95] M. S. Allen, M. W. Sracic, S. Chauhan, and M. H. Hansen. Output-only modal analysis of linear time-periodic systems with application to wind turbine simulation data. Mechanical Systems and Signal Processing, 25:1174–1191, 2011. 152 [96] M. S. Allen. Global and Multi-Input-Multi-Output (MIMO) Extensions of the Algorithm of mode isolation (AMI). PhD thesis, Georgia Institute of Technology, 2005. [97] S. Schreck and M. Robinson. Rotational augmentation of horizontal axis wind turbine blade aerodynamic response. Wind Energy, 5:133–150, 2002. [98] A. C. Hansen and C. P. Butterfield. Aerodynamics of horizontal axis wind turbines. Annual Review of Fluid Mechanics, 25:115–149, 1993. [99] S. Schreck and M. Robinson. Horizontal axis wind turbine blade aerodynamics in experiments and modeling. In IEEE Transactions on Energy Conversion, volume 22, 2007. [100] J. G. Leishman. Challenges in modelling the unsteady aerodynamics of wind turbines. Wind Energy, 5:85–132, 2002. [101] D. M. Eggleston and F. S. Stoddard. Wind Turbine Engineering Design. Van Nostrand Reinhold Company, New York, 1987. [102] H. Snel. Review of aerodynamics for wind turbines. Wind Energy, 6:203–211, 2003. [103] D. Shipley, M. Miller, M. Robinson, M. Luttges, and D. Simms. Evidence that aerodynamic effects, including dynamic stall, dictate HAWT structure loads and power generation in highly transient time frames. Technical report, National Renewable Energy Laboratory, Golden, CO, 1994. [104] M. Robinson, R. Galbraith, D. Shipley, and M. Miller. Unsteady aerodynamics of wind turbines. In AIAA Paper 95-0526, 1995. [105] A. Eggers and R. Digumarthi. Approximate scaling of rotational effects on mean aerodynamic moments and power generated by CER blades operating in deep-stalled flow. In 11th ASME Wind Energy Symposium, 1992. [106] J. Tangler and M. Selig. An evaluation of an empirical model for stall delay due to rotation for HAWTs. Technical report, National Renewable Energy Laboratory, Golden, CO, 1997. [107] M. H. Hansen, K. Thomsen, and P. Fuglsang. Two methods for estimating aeroelastic damping of operational wind turbine modes from experiments. Wind Energy, 9:179– 191, 2006. 153 [108] R. A. Caldwell.Jr. and B. F. Feeny. Output-only modal analysis of a nonuniform beam experiment by using decomposition methods. ASME IDETC’ 11, Washington, August, pages 28–31, 2011. [109] L. Meirovitch. Fundamentals of Vibrations. McGraw-Hill, 2001. [110] L. Meirovitch. Analytical Methods in Vibrations. Macmillan, 1967. [111] J. H. Ginsberg. Mechanical and Structural Vibrations, Theory and Applications. John Wiley and Sons, New York, 2001. [112] R. L. Bisplinghoff, H. Ashley, and R. L. Halfman. Aeroelasticity. Addison-Wesley Publishing Company, 1983. [113] Y. Y. Kim. Flexural-torsional coupled vibration of rotating beams using orthogonal polynomials. Master’s thesis, Virginia Polytechnic Institute and State University, 2000. [114] I. T. Minhinnick. The theoretical determination of normal modes and frequencies of vibration. Technical report, Aeronautical Research Council, London, 1956. [115] N. M. Auciello. Flapwise bending vibration of rotating Euler-Bernoulli beam with non-uniform tapers. In 1st International Virtual Scientific Conference, Volume 1, June2013, 2013. [116] S. P. Timoshenko. History of Strength of Materials. Dover Publications Inc, New York, 1953. [117] Dary L. Logan. A First Course in the Finite Element Method. PWS Publishing Company, 1997. [118] Jason Salmanoff. A finite element, reduced order, frequency dependent model of viscoelastic damping. Master’s thesis, Virginia Polytechnic Institute and State University, 1997. [119] Randall J. Allemang. The modal assurance criterion - twenty years of use and abuse. Sound and Vibration, 2003. [120] K. J. Liu. A theoretical study on distributed active vibration suppression of mechanical systems with integrated distributed piezoelectric sensor and actuator. Master’s thesis, Univ. of Kentucky, 1989. 154 [121] H. Tzou, D. Johnson, and K. Liu. Damping behavior of cantilevered structronic systems with boundary control. Journal of Vibration and Acoustics, 121:402–407, 1999. [122] J. S. Bae, M. K. Kwak, and D. J. Inman. Vibration suppression of a cantilever beam using eddy current damper. Journal of Sound and Vibration, 284:805–824, 2004. [123] A. V. Oppenheim and R. W. Schafer. Discrete-Time Signal Processing. Prentice Hall, Englewood Cliffs, NJ, 1989. [124] L. Meirovitch. Principles and Techniques of Vibrations. Upper Saddle River, N.J. : Prentice Hall, 1997. [125] T. S. Reddy and K. R. Kaza. A comparative study of some dynamic stall models. Technical report, NASA Technical Memorandum 88917, 1987. [126] J. Larsen, S. Nielsen, and S. Krenk. Dynamic stall model for wind turbine airfoils. Journal of Fluids and Structures, 23:959–982, 2007. [127] J. Manwell, J. MacGowan, W. Stein, A. Rogers, and M. Jones. Development of quasisteady analytical models for wind/diesel systems. Proc. Ninth ASME Wind Energy Symposium, ASME Publication SED, 9:223–230, 1990. [128] J. Tangler and G. Bir. Evaluation of RCAS Inflow Models for Wind Turbine Analysis. DIANE Publishing, 2004. [129] S. Jung, T. No, and K. Ryu. Aerodynamic performance prediction of a 30 KW counterrotating wind turbine system. Renewable Energy, 30:631–644, 2005. [130] W. Skrzypinski. Analysis and Modeling of Unsteady Aerodynamics with Application to Wind Turbine Blade Vibration at Standstill Conditions. PhD thesis, DTU Wind Energy, 2012. [131] X. Chen and R. Agarwal. Inclusion of a simple dynamic inflow model in the blade element momentum theory for wind turbine application. International Journal of Energy and Environment, 5:183–196, 2014. [132] T. Theodorsen. General theory of aerodynamic instability and the mechanism of flutter. Technical report, NACA Report 496, 1935. [133] C. T. Tran and D. Petot. Semi-empirical model for the dynamic stall of airfoils in view of the application to the calculation of responses of a helicopter blade in forward flight. 155 The Sixth European Rotorcraft and Powered Lift Aircraft Forum, Bristol, England, 1980. [134] K. McAlister, O. Lambert, and D. Petot. Application of the ONERA model of dynamic stall. Technical report, NASA Technical Paper 2399, AVSCOM Technical Report 84A-3, 1984. [135] D. A. Peters. Toward a unified lift model for use in rotor blade stability analysis. Journal of the American Helicopter Society, 30:32–42, 1985. [136] R. L. Halfman. Experimental aerodynamic derivatives of a sinusoidally oscillating airfoil in two-dimensional flow. Technical report, NACA-TR-1108, 1951. [137] S. L. Brunton and C. W. Rowley. Empirical state-space representations for Theodorsens lift model. Journal of Fluids and Structures, 38:174–186, 2013. [138] J. P. Rogers. Applications of an analytic stall model to time-history and eigenvalue analysis of rotor blades. Journal of the American Helicopter Society, 29:25–33, 1982. [139] A. H. Nayfeh and D. T. Mook. Nonlinear Oscillations. Wiley, NewYork, 1979. 156