THE INERTER PENDULUM VIBRATION ABSORBER: WITH APPLICATIONS IN OCEAN WAVE ENERGY CONVERSION AND HYDRODYNAMIC RESPONSE SUPPRESSION By Aakash Gupta A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering—Doctor of Philosophy 2023 ABSTRACT The annual power incident on the ocean-facing coastlines of North America is over 400 GW. Capturing a small fraction of this energy can significantly contribute to meeting energy demands. Therefore, there is a renewed research interest in converting energy from ocean waves. Typically, ocean wave energy capturing devices, known as wave energy converters (WECs), are placed in deep water as the wave energy is higher in the deep water compared to shallow water. To reduce the cost of installing and maintaining WECs in deep water, they can be integrated with existing offshore floating platforms in the ocean. For such integration, traditional WECs, operating on the principle of linear resonance, have a natural period in heave close to a typical wave period to generate a large heave resonant response and hence high-efficiency wave power production, which causes large platform motions. In other words, wave power production and hydrodynamic stability of the platform are conflicting objectives in traditional linear WECs. Therefore, simultaneous wave energy conversion and response suppression of the platform is necessary. To address this issue, in this work, a device known as an inerter pendulum vibration absorber (IPVA) is proposed combining the inerter with a parametrically excited centrifugal pendulum. Two system variations are studied: the IPVA and IPVA-PTO, marking the absence and presence of an electromagnetic power take-off (PTO) system. Both the IPVA and the IPVA-PTO are integrated with a single-degree-of-freedom (sdof) structure: a primary mass, and a spar, respectively. The efficacy in suppressing vibrations is studied in the case of the sdof IPVA system, whereas wave energy conversion and response suppression are analyzed for the spar IPVA-PTO. For both systems, a nonlinear energy transfer phenomenon in which the energy is transferred between the primary mass (or spar) and the pendulum vibration absorber. For the sdof IPVA system, it is shown that the energy transfer is associated with the 1:2 internal resonance of the pendulum induced by a period-doubling bifurcation. A perturbation analysis shows that a pitchfork bifurcation and a period-doubling bifurcation are necessary and sufficient conditions for this internal resonance to occur. Harmonic balance analysis, in conjunction with Floquet theory, along with the arc-length continuation scheme, is used to predict the boundary of internal resonance in the parameter space and verify the perturbation analysis. Furthermore, the effects of various system parameters on the boundary are examined. Next, the sdof IPVA is compared with a linear benchmark and an autoparametric vibration absorber and shows more efficacious vibration suppression. For the spar IPVA-PTO system, a similar analysis shows the nonlinear energy transfer, which is used to convert the vibrations of the spar into electricity while reducing its hydrodynamic response. Similar to the IPVA, a period-doubling bifurcation results in 1:2 internal resonance, which is necessary and sufficient for nonlinear energy transfer to occur. The hydrodynamic coefficients of the spar are computed using a commercial boundary element method code. The period-doubling bifurcation is studied using the harmonic balance method. A modified alternating frequency/time (AFT) approach is developed to compute the Jacobian matrix involving nonlinear inertial effects of the IPVA-PTO system. The response amplitude operator (RAO) in heave and the capture width of the spar IPVA-PTO are compared with its linear counterpart, and the spar IPVA-PTO system outperforms the linear energy harvester with lower RAO and higher capture width. Experiments containing integration of the IPVA and the IPVA-PTO system with an sdof system (or “dry" spar in the case of IPVA-PTO) are performed in order to verify the analysis. Next, both the IPVA and the IPVA-PTO systems are integrated with a spar-floater combination and analyzed for their performance. Near the first resonance frequency, the spar-floater IPVA system shows a period-doubling bifurcation and energy transfer similar to the sdof IPVA system and outperforms the linear benchmark for hydrodynamic response suppression. On the other hand, the spar-floater integrated IPVA-PTO system is analyzed for its performance near both resonance frequencies. It is shown that near the first resonance, the spar-floater IPVA-PTO system’s response undergoes a period-doubling bifurcation, and for small electrical damping, shows energy transfer. However, near the second resonance, secondary Hopf bifurcation is observed. A rich set of pendulum responses, such as primary and secondary harmonics, quasi-periodic, non-periodic, and rotation, are observed. Rotation is shown to provide the best energy conversion among all the identified responses. Finally, the electrical damping of the system is varied to find the optimal values for which the largest energy conversion occurs in the system, and it is found that the optimal electrical damping for energy transfer is associated with pendulum’s rotation. Copyright by AAKASH GUPTA 2023 To Akshat, in loving memory v ACKNOWLEDGEMENTS I would like to express my sincere gratitude toward my advisor Dr. Wei-Che Tai, whose mentorship and support helped me throughout my doctoral studies. His enormous knowledge and enthusiasm have greatly helped me throughout my studies. I sincerely thank my committee members, Dr. Zhaojian Li and Dr. Xiaobo Tan. I would also like to thank Dr. Brian Feeny for his contribution to my academic and personal development. I want to mention VibeLab members Dr. Daniel Segalman, Fatemeh, Jamal, Gaurav, Ayse, Jun, Joel, Guoxin, and Majid, for cultivating a great work environment and for the time spent together in biking, hiking, camping, and movie nights. I would also like to thank all my friends, especially Madhumitha, Parth, Anshu, Prashant, Saharan, Rahul, Santhosh, Mrudul, Yoganandh, Saru, Varun, and, Prashasti, for their support throughout my Ph.D. studies. I thank my parents, Meena and Purshottam, and my siblings and cousins, Aadarsh, Anuj, Chinki, Muskaan, and Anmol, for their constant love and support. Finally, I would like to extend my special thanks to Varun Thumbe, Anirudh Suresh, and Saurabh Pargal, whose mathematical and intuitive discussions helped me realize my love for research and motivated me to keep learning new things. vi LIST OF TABLES . . . LIST OF FIGURES . . . . . . . . . . . . . TABLE OF CONTENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Overview of the work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Background . . 1.3 Motivation and literature survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Outline of the document . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 1 1 2 2 9 CHAPTER 2 . . . . ANALYSIS OF THE IPVA SYSTEM: A STUDY IN VIBRATION . 10 SUPPRESSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Overview . . 10 2.2 Design of the Inerter Pendulum Vibration Absorber . . . . . . . . . . . . . . 2.3 Internal resonance of IPVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 Numerical demonstration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5 Parametric studies . . 2.6 Discussion . . 25 2.7 Experimental analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . OCEAN WAVE ENERGY CONVERSION USING THE IPVA-PTO SYSTEM WITH A SPAR . . . . . . . . . . . . . . . . . . . . . . . . . 40 . 40 3.1 Overview . 3.2 The IPVA-PTO system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3 . 47 Internal resonance and bifurcation boundaries . . . . . . . . . . . . . . . . . 3.4 Numerical demonstration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 . 54 3.5 Discussion . 3.6 Experimental analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . INTEGRATION OF THE IPVA WITH A SPAR-FLOATER SYSTEM: A STUDY IN VIBRATION SUPPRESSION . . . . . . . . . . . . . . . 67 4.1 Overview . . 67 4.2 Design of the system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.3 Equations of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.4 Period-doubling bifurcation in the IPVA system . . . . . . . . . . . . . . . . . 72 4.5 Nonlinear energy transfer and energy harvesting potential . . . . . . . . . . . . 73 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.6 Parametric studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 . 4.7 Conclusion . . . . CHAPTER 5 INTEGRATION OF THE IPVA-PTO WITH A SPAR-FLOATER SYSTEM: STUDY IN OCEAN WAVE ENERGY CONVERSION . . . . . . . . . . 80 . 80 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Overview . vii 5.2 Modeling of the system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.3 Bifurcation analysis of the IPVA-PTO system . . . . . . . . . . . . . . . . . . 84 5.4 Performance analysis near first resonance frequency . . . . . . . . . . . . . . . 86 5.5 Performance analysis near second resonance frequency . . . . . . . . . . . . . 90 5.6 Parametric studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.7 Optimization of electrical damping value . . . . . . . . . . . . . . . . . . . . . 95 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.8 Conclusion . . . . . . CHAPTER 6 CONCLUSION AND PROPOSAL FOR THE FUTURE WORK . . . . 6.1 Concluding remarks . 6.2 Suggested future direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 . 103 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 APPENDIX A: BIFURCATION TRACKING ALGORITHM . . . . . . . . . . . . . . . . . 114 APPENDIX B: JACOBI-ANGER EXPANSION . . . . . . . . . . . . . . . . . . . . . . . . 116 viii LIST OF TABLES Table 2.1 Labels depicting the different parts of the experimental setup . . . . . . . . . . . 33 Table 2.2 Experimental parameters for the IPVA system . . . . . . . . . . . . . . . . . . . 35 Table 3.1 Parameters for the experimental IPVA-PTO system . . . . . . . . . . . . . . . . 61 Table 4.1 Parameters for spar-floater system (all data in SI units) . . . . . . . . . . . . . . 70 ix LIST OF FIGURES Figure 1.1 Annual mean power density and annual mean best direction (→) for waves around the world [4] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 1.2 A point absorber wave energy converter which generates energy using the heave motion [7] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 Figure 2.1 Rack-pinion and ball-screw IPVA systems: schematics and computer-aided (CAD) realizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Figure 2.2 Pitchfork bifurcation of 𝜙0 for the following parameter values 𝑓 = 0.007, 𝜇𝑟 = 0.2, 𝜔 = 0.9, 𝜉 = 𝜉 𝑝 = 0.005 . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Figure 2.3 Parametric instability boundary for the following parameter values 𝜂 = 0.3, 𝜇𝑟 = 0.25, 𝜉 = 0.005, 𝜉𝑝 = 0.005 . . . . . . . . . . . . . . . . . . . . . . 20 Figure 2.4 FFT and time series of periodic solutions at point 1 in Fig. 2.3 . . . . . . . . . . 21 Figure 2.5 FFT and time series of periodic solutions at point 2 in Fig. 2.3 . . . . . . . . . . 22 Figure 2.6 FFT and time series of non-periodic solutions at point 3 in Fig. 2.3 . . . . . . . 22 Figure 2.7 Poincare section at points 1, 2, and 3 as marked in Fig. 2.3. The value of 𝜙 is normalized between 0 and 2𝜋 . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Figure 2.8 Parametric instability boundary for 𝜇𝑟 = 0.4, 𝜉 = 0.005, 𝜉 𝑝 = 0.01 for various values of 𝜂 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Figure 2.9 Parametric instability boundary for 𝜂 = 0.4, 𝜉 = 0.005, 𝜉𝑝 = 0.01 for various values of 𝜇𝑟 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Figure 2.10 Parametric instability boundary for 𝜂 = 0.3, 𝜇𝑟 = 0.2, 𝜉 𝑝 = 0.005 for various values of 𝜉 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Figure 2.11 Parametric instability boundary for 𝜂 = 0.4, 𝜇𝑟 = 0.4, 𝜉 = 0.005 for various values of 𝜉 𝑝 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Figure 2.12 Comparison of IPVA with the linear system for various parameters . . . . . . . 26 Figure 2.13 The autoparametric vibration absorber [57] . . . . . . . . . . . . . . . . . . . . 27 Figure 2.14 Comparison of IPVA and the autoparametric vibration absorber for a set of parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Figure 2.15 Comparison of IPVA and the autoparametric vibration absorber for a set of parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 x Figure 2.16 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 2.17 Numerical stability boundary of energy transfer corresponding to the experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 Figure 2.18 Time series and measured FFT of the pendulum motion . . . . . . . . . . . . . 36 Figure 2.19 Frequency response of 𝜃 and 𝜙 for the IPVA and the linear system (where applicable) at 0.2 𝑔 base excitation acceleration . . . . . . . . . . . . . . . . . 38 Figure 2.20 Frequency response of 𝜃 and 𝜙 for the IPVA and the linear system (where applicable) at 0.21 𝑔 base excitation acceleration . . . . . . . . . . . . . . . . . 38 Figure 3.1 The IPVA-PTO system integrated with a spar: schematic and CAD realization, and the benchmark linear system . . . . . . . . . . . . . . . . . . . . . . . . . 41 Figure 3.2 Comparison of added mass and radiation damping between Yeung et al. [77] and Ansys simulation for 𝑎 = 0.2m, ℎ = 1m, 𝑑 = 0.25m . . . . . . . . . . . . . 46 Figure 3.3 Ansys AQWA model for calculation of hydrodynamic coefficients . . . . . . . . 46 Figure 3.4 Added mass, radiation damping, and diffraction and Froude-Krylov force for sparD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Figure 3.5 Parametric instability boundary for 𝜂 = 0.4, 𝜇𝑟 = 0.4, 𝜇𝑔 = 0.03, 𝜉 = 0.05, 𝜉𝑝 = 0.01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Figure 3.6 FFT and time series of periodic solutions at point 1 in Fig. 3.5 . . . . . . . . . . 54 Figure 3.7 FFT and time series of periodic solutions at point 2 in Fig. 3.5 . . . . . . . . . . 55 Figure 3.8 FFT and time series of periodic solutions at point 3 in Fig. 3.5 . . . . . . . . . . 55 Figure 3.9 Comparison of response amplitude operator (RAO) and capture width between the IPVA-PTO and linear system . . . . . . . . . . . . . . . . . . . . . . . . . 57 Figure 3.10 Experimental setup. The labels denote the following — 1: Shaker, 2: Spectral analyzer and controller, 3: Accelerometer signal conditioner, 4: Shaker signal amplifier, 𝐴: Top plate, 𝐵: Ball-screw system, 𝐶: Ball-screw mounting plate , 𝐷: Sun and planet gear, 𝐸: Middle plate, 𝐹: Base plate, 𝐺: Primary mass, 𝐻: Connecting springs, 𝐼: Pendulum, 𝐽: Carrier, 𝐾: Generator, 𝐿: Load resistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 Figure 3.11 PD boundary of energy transfer corresponding to the experimental setup for two values of load resistance . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Figure 3.12 Measured time series and FFT of the output voltage at 0.42 𝑔, at 3.6 Hz . . . . . 62 xi Figure 3.13 Frequency response of relative motion between the top plate and the base plate, and power for the IPVA-PTO and the linear system at 0.6 𝑔 base excitation acceleration with 10 ohm resistance. 1:2 stands for secondary harmonic solutions, NP stands for non-periodic solution . . . . . . . . . . . . . 64 Figure 3.14 Frequency response of relative motion between the top plate and the base plate, and power for the IPVA-PTO and the linear system at 0.9 𝑔 base excitation acceleration with 5 ohm resistance and 3 ohm resistance for the linear system and 5 ohm for the IPVA system . . . . . . . . . . . . . . . . . . . 64 Figure 4.1 Ocean wave energy harvesting design . . . . . . . . . . . . . . . . . . . . . . . 67 Figure 4.2 Ansys AQWA model for calculation of hydrodynamic coefficients . . . . . . . . 70 Figure 4.3 Added mass, radiation damping, and diffraction and Froude-Krylov force for spar and floater . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 Figure 4.4 Stability boundary for the primary harmonics of the system for 𝜇𝑟 = 0.2, 𝜂 = 0.3, 𝜉 = 0.05, 𝜉𝑝 = 0.02, 𝜇𝑝 = 0.039, 𝜇𝑏𝑠𝑐 = 0.1 . . . . . . . . . . . . . . . . 74 Figure 4.5 Time series and FFT of the pendulum motion at point ×1 . . . . . . . . . . . . 74 Figure 4.6 Time series and FFT of the pendulum motion at point ×2 . . . . . . . . . . . . 75 Figure 4.7 Time series and FFT of the pendulum motion at point ×3 . . . . . . . . . . . . 75 Figure 4.8 The frequency responses of the linear and the IPVA system . . . . . . . . . . . 76 Figure 4.9 Effect of various system parameters on the stability boundary, all the values are from Fig. 4.4 except for the parameter value specified in the figure . . . . . . 79 Figure 5.1 The IPVA-PTO system integrated in a spar-floater setup for wave energy conversion where: a. Ocean wave energy conversion setup, b. the IPVA-PTO system, and c. equivalent mathematical model . . . . . . . . . . . . . . . . . . 80 Figure 5.2 Ansys Aqwa mesh model, added mass, radiation damping, and diffraction and Froude-Krylov force for spar and floater . . . . . . . . . . . . . . . . . . . 82 Figure 5.3 Bifurcation boundary for the primary harmonics of the system for 𝜉 = 0.01, 𝜉𝑒 = 0.015, 𝜉𝑒=0.02, 𝜇𝑟=0.4, 𝜂=0.4, 𝜇 𝐴∞,1=0.0554, 𝜇 𝐴∞,2=0.6381 , 𝜔𝑟=6.928, 𝜔0=3.20 rad/s, 𝜇𝑔=0.01, 𝜇 𝑝=0.02, 𝜇𝑏𝑠𝑐=0.04, 𝜇𝐹=0.2097. Markers represent the following: Primary — ‘+’, Chaotic — ‘*’, One-fifth subharmonics — ‘★’, Secondary — ‘□’, Rotation — ‘◦’, Intermittent rotation — ‘△’ . . . . . 86 Figure 5.4 Real and imaginary parts of eigenvalues whose magnitude equals one with respect to frequency, other eigenvalues not shown . . . . . . . . . . . . . . . . 87 xii Figure 5.5 Poincare sections corresponding to points 1 and 2 in Fig. 5.3 show that the secondary Hopf bifurcation results in quasi-periodic motion . . . . . . . . . . . 87 Figure 5.6 Period-doubling bifurcation boundary for two values of electrical damping with mechanical damping 𝜉 = 0.02. All other parameters are from Fig. 5.3 . . . 88 Figure 5.7 Comparison between RAO and normalized capture width of the IPVA-PTO and the linear system for different wave heights for 𝜉 = 0.02, 𝜉𝑒,𝑙 = 0.0216 and 𝜉𝑒 = 0.01, all other parameters from caption of Fig. 5.3 . . . . . . . . . . . 90 Figure 5.8 Comparison between RAO and normalized capture width of the IPVA-PTO and the linear system for different wave heights for 𝜉 = 0.02, 𝜉𝑒,𝑙 = 0.0216 and 𝜉𝑒 = 0.02, all other parameters from caption of Fig. 5.3 . . . . . . . . . . . 90 Figure 5.9 Comparison between RAO and normalized capture width of the IPVA-PTO and the linear system for different wave heights for 𝜉𝑒,𝑙 = 0.07, and the boundary corresponding to Fig. 5.3. All marked wave heights in cm . . . . . . 91 Figure 5.10 Comparison between RAO and normalized capture width of the IPVA-PTO and the linear system for different wave heights for 𝜉 = 0.01, 𝜉𝑒,𝑙 = 0.07 and 𝜉𝑒 = 0.02, all other parameters from caption of Fig. 5.3. The corresponding bifurcation boundary can be found in Fig. 5.11d . . . . . . . . . . . . . . . . . 92 Figure 5.11 Effect of various system parameters on the bifurcation boundary, all the values are from Fig. 5.3 except for the parameter value specified in the figure. Readers are referred to the web version of this article for clarity . . . . . . . . . 95 Figure 5.12 ‘◦’ represents the Pareto frontal (non-dominated) solutions, ‘∗’ represents all the solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Figure 5.13 maxRAO and H2CW as a function of the electrical damping value for wave height = 5 cm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Figure 5.14 maxRAO and H2CW as a function of the electrical damping value for wave height = 4 cm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Figure 5.15 RAO and normalized capture width for two different values of 𝜉𝑒 for wave height 5 cm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 Figure 5.16 RAO and normalized capture width for two different values of 𝜉𝑒 for wave height 4 cm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 Figure 6.1 Comparison 𝜃, 𝜓 and 𝜙 responses between harmonic balance method and direct numerical simulation for wave height 3.7 cm and excitation frequency of 1.12 Hz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 . . . . . . xiii Figure 6.2 Stability boundaries for primary rotations and the response of the pendulum for wave height = 3.7 cm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 xiv CHAPTER 1 INTRODUCTION 1.1 Overview of the work This work proposes a novel system known as inerter-based vibration energy converter, called the inerter pendulum vibration absorber (IPVA), consisting of an inerter and a pendulum vibration absorber. The main objective of this study is to analyze the IPVA system for its effectiveness in simultaneous vibration suppression and energy harvesting. The harmonic balance method, along with a modified alternating-frequency time (AFT) method, is used to analyze the system. The stability of the periodic response of the system is determined using the Floquet theory, and it is observed that the primary harmonic response of the system can bifurcate either via a period- doubling bifurcation or secondary Hopf bifurcation depending on the type of integration between the IPVA (with or without generator) and linear system (single or two degree-of-freedom). This work presents two main systems: the IPVA system and the IPVA-PTO system. The significant difference between these two systems is that the IPVA-PTO system has a generator integrated with the pendulum and is studied for vibration suppression and energy harvesting. In contrast, the IPVA system does not have a generator and is primarily concerned with vibration suppression. First, the IPVA system is integrated with a single-degree-of-freedom (sdof) system to understand its effect on the performance of the sdof system. After this, the IPVA-PTO is integrated with an offshore floating platform (spar) to understand the impact of added mass and radiation damping on the system. For the IPVA system, it has been observed that it outperforms an autoparametric system in terms of the motion suppression of the primary system and energy harvesting potential. Furthermore, for the IPVA-PTO system, the response amplitude operator (RAO) in the heave of the spar and the capture width (hydrodynamic energy) of the IPVA-PTO- integrated spar are better than a benchmark linear PTO as the former has a lower RAO and higher capture width. Experiments are performed on the sdof “dry” IPVA and “dry” IPVA-PTO systems to observe their vibration suppression and energy conversion effect. After this, both the IPVA and the IPVA-PTO is integrated with a two-degree-of-freedom wave energy converter consisting of a 1 spar and a floater to analyze their efficacy in hydrodynamic response mitigation and wave energy conversion compared to a linear benchmark. 1.2 Background Since the last decade, renewable energy has received a lot of interest in generating electricity due to the depletion of fossil fuels and the threat of climate change. Of all the most popular renewable energy techniques, solar, wind and wave, wave energy is the most spatially concentrated, with an average power of 2 − 3 kW/m2 just below the ocean surface on the area perpendicular to the direction of the wave propagation [1]. Many researchers have worked to quantify the wave energy resources present in the oceans [2, 3, 4] and obtained similar results. One such quantification is shown in Fig. 1.1. As seen from this figure, ocean wave energy resources are enormous over the coastline of North America. It is estimated that the annual average wave power incident on the ocean-facing coastlines of North America is over 400 GW (about 80% electricity consumption for the entire continent [4]). Despite enormous resources, the cost of using existing wave energy converters (WECs) to generate electricity from ocean waves is higher than that of solar and wind energy conversion technologies due to the installation, mooring/foundation costs, operation, and maintenance costs, which account for 40%–50% of the life costs of the wave energy project [5]. A promising way to reduce these costs is to integrate WECs with offshore floating platforms because they can share infrastructure, equipment, mooring and anchoring systems, and survey and monitoring methods [6]. One of the ways to do this is to integrate WECs with oil rigs so that they can generate energy and directly supply electricity to the rigs. 1.3 Motivation and literature survey Despite efforts in reducing the cost of wave energy conversion technology, the LCOE (Levelized Cost of Energy) of wave energy (around 570 $/MWh) is significantly larger than that of onshore wind or solar photovoltaic energy (≈30 $/MWh) [8, 9]. To reduce the LCOE, researchers are investigating ways to integrate wave energy converters with existing offshore structures in the ocean [6]. Structures like floating wind turbines [10] and spar platforms used in offshore oil and gas industry [11, 12, 13, 14] are feasible candidates for such integration. Specifically, spar platforms 2 Figure 1.1 Annual mean power density and annual mean best direction (→) for waves around the world [4] Figure 1.2 A point absorber wave energy converter which generates energy using the heave motion [7] 3 establish the buoyancy and stability on a long and slender cylinder that goes deep below the water surface, thereby having good hydrodynamic response/stability and large water depths (600∼2500 meters for spar platforms in the Gulf of Mexico [15]). On the other hand, heaving WECs, like the one shown in Fig. 1.2, convert the relative heave motion between oscillating bodies into electricity and have a high wave energy conversion efficiency when operating at resonance [16]. As wave energy resources are more abundant in deep water than in shallow water, it is viable to integrate WECs and spar platforms to lower the LCOE. The integration of the spar and different types of heaving WECs [11, 12, 13, 17] has been getting attention. Specifically, the integration of spar with an annular floater has been widely studied to assess its feasibility in reducing the cost while converting power; see [11, 13, 18, 19], for example. Existing numerical studies [11, 17] suggest that such integration can lead to a 7%–30% capture width ratio (hydrodynamic efficiency) of wave energy production, which is comparable to existing heaving WECs [20]. According to the scaling law in [20], heaving WECs of a larger diameter would have a higher capture width ratio. A typical spar platform in the Gulf of Mexico, e.g., the Horn Mountain, has a diameter of 30 meters (BSEE data [21]). If the spar-WEC integration in [11, 17] were scaled up to this diameter, the peak mean wave power in operational conditions would be 2.4–10 MW (current floating wind turbines have 5 MW wind power). Although it shows promising wave energy conversion, such integration does not ensure a good hydrodynamic response of the platform. Past studies have shown that the integration with heaving WECs amplifies the platform heave and pitch motion [11, 12, 13], and even causes Mathieu instability [19], which would aggravate fatigue of the mooring and riser systems and even lead to failure of the whole system [22, 23]. This deterioration of hydrodynamic response and stability can be explained as follows. Generally speaking, a spar platform has a 20–30 s heave natural period [24, 25], which is far away from typical incident wave periods (5–10 s [26]) to avoid large heave resonant response. On the other hand, traditional heaving WECs operate based on the basic principle of linear resonance, thereby having a natural period in heave close to a typical wave period to generate a large heave resonant response and hence high-efficiency wave power 4 production. When a heaving WEC is integrated with a platform, this significant heave resonant response can give rise to large platform heave/pitch motions. In other words, increasing wave power production is at the cost of deteriorating the hydrodynamic response of the platform, which is the fundamental problem this research aims to resolve: wave power production and hydrodynamic response reduction are conflicting objectives. Moreover, another issue with traditional WECs is that the half-power frequency bandwidth of the WECs is low, which is an issue since the dominant wave period typically changes with time. Therefore, researchers are looking at ways to convert wave energy over a broadband frequency range without sacrificing the stability of the spar. Fantai et al. [27] used active control in co-located offshore wind-wave systems. They showed that actively altering the wave field with a WEC array with model predictive control before being incident on a floating offshore wind turbine (FOWT) can result in both motion reduction and reliable energy conversion. Zhang et al. [28] demonstrated that liquid column dampers can be used to maintain the stability of FOWTs. Another direction being investigated to widen the frequency bandwidth is the use of nonlinearity in the system. For example, nonlinear stiffness mechanisms [29, 30, 31] have been investigated to enhance the frequency bandwidth of point absorbers. Methods including increasing the number of degrees of freedom in the point absorber to increase wave energy conversion efficiencies while reducing the resonance frequency for spherical submerged bodies [32] have also been proposed. However, for two-body wave energy converters, no passive means exist that can achieve simultaneous hydrodynamic response suppression and wave energy conversion at the same time. Much research has been published considering nonlinearities to broaden the bandwidth of energy converters. For example, energy converters with intentionally introduced nonlinear stiffness have been widely studied [33]. Although this approach is proven effective, it is hard to realize the required nonlinear stiffness in large-scale structures. On the other hand, when rotary electromagnetic energy converters are considered, inertia nonlinearity can be implemented via creative transmission mechanisms. For example, Zuo et al. [34] developed the mechanical transmission concept known as mechanical motion rectification (MMR) for energy conversion. The MMR mechanism changes 5 the inertia of the energy converter in a piece-wise manner. Several studies [35, 36, 37] have demonstrated that such piece-wise inertia can broaden the bandwidth. Because wave energy conversion is generally designed for large-scale energy harvesting, other types of inertia nonlinearity deserve more investigation. One of the ways for large-scale energy conversion is using inertial nonlinearity by making the use of a device called inerter coupled with nonlinear vibration absorbers. An inerter is a mechanical device with two terminals, each of which exerts an equal and opposite inertial force proportional to the relative acceleration between the terminals [38]. The inerter amplifies the inertial effects of a small mass by using motion transmission mechanisms, fluids, and levers [38]. By virtue of its mass amplification effect, the inerter has been studied to improve the performance of various passive vibration mitigation techniques in the last decade. Ikago et al. [39] developed the tuned viscous mass damper (TVMD), which consisted of a tuning spring in series with the inerter and a viscous damper in parallel. It was shown that the TVMD outperformed the viscous damper alone when applied to a seismically excited single-degree-of- freedom (sdof) structure. Furthermore, Lazar et al. [40] proposed the tuned inerter damper (TID), wherein the inerter was substituted for the oscillating mass of a tuned mass damper (TMD). The TID and TMD were compared in seismically excited multiple-degree-of-freedom (MDOF) structures and demonstrated similar effectiveness. Later, Lazar et al. [41] considered the TID in suppressing the midspan vibration of cables and showed that the TID outperformed the optimal viscous damper. Moreover, Qian et al. [42] studied serial and parallel connections between the TID and a base-isolation system and concluded that the serial TID outperformed the parallel TID for practical structures. The inerter has also been applied to enhance the inertial effects of dynamic vibration absorbers (DVAs). Marian and Giaralis [43] proposed the tuned mass damper inerter (TMDI), which consisted of a TMD and the inerter in series. In a 3DOF structure simulation, they showed that for achieving similar vibration control performance, the weight of the TMDI was four times lighter than the TMD. Furthermore, De Domenico and Ricciardi [44] incorporated the TMDI in a base-isolation system and demonstrated that the displacement demand of the base-isolated structure could be 6 significantly reduced. Moreover, Joubaneh and Barry [45] studied the performance of four models of electromagnetic resonant shunt TMDI (ERS-TMDI) on both vibration suppression and energy harvesting and identified the best model. Their parametric studies showed that increasing the inertance enhances the performance of the best model in terms of both vibration mitigation and energy harvesting. On the other hand, Tai [46] proposed the tuned inerter-torsional-mass damper (TITMD), which integrated the inerter and a torsional mass damper. Compared with the TMDI, the TITMD achieved 20 − 70% improvement with identical weights. For the case of two-body wave energy converters (like spar-floater systems), Asai et al. [47] found that the addition of an inerter to the system enhances the energy conversion efficiency and broadens the frequency bandwidth. In recent years, the inerter has been integrated with nonlinear vibration absorbers. Qian and Zuo [48] considered the effects of adding an inerter to a nonlinear vibration absorber (NVA). The nonlinear vibration absorber consisted of a tuned mass damper with a nonlinear spring containing both linear and cubic stiffness. They observed that the spring-inerter-damper (SID) system added to the beam outperformed the nonlinear vibration absorber without the inerter. Furthermore, Kakou and Barry [49] added a nonlinear spring to the electromagnetic resonant shunt tuned mass damper- inerters (ERS-TMDI) [45] to analyze the implications of coupling a nonlinear spring to the system. Two configurations of the system, one with the energy harvester between the tuned mass and ground (Configuration-1) and the other with the energy harvester between primary structure mass and the tuned mass (Configuration-2), were compared for their efficacy in vibration suppression and energy harvesting of the system. It was found that Configuration-1 exhibits a higher range of feasible forcing without degrading the performance compared to Configuration-2. It was also observed that for optimal Configuration-1, higher nonlinear stiffness, inerter magnitude and resistance, and lower capacitance and inductance improved the energy harvesting performance. Yang et. al [50] proposed the nonlinear inertance mechanism (NIM) created by combining oblique inerters with one common hinged terminal and the other terminals fixed. It was shown that the addition of NIM could enhance the vibration isolation capabilities of a system. The NIM was combined with two different isolators: a spring-damper isolator and a nonlinear quasi-zero-stiffness (QZS) isolator. 7 After adding the NIM, the linear isolator showed bending of the frequency response curve towards the low-frequency range and reduction of the original peak values in dynamic response. For QZS systems, the larger frequency range of small dynamic response amplitude and lower kinetic energy of the mass was observed after the addition of NIM. In this work, we integrate the inerter, and a centrifugal pendulum with a generator attached such that the centrifugal pendulum is parametrically excited by the inerter. The integrated system is referred to as the inerter pendulum vibration absorber power take-off (IPVA-PTO). Various systems wherein a pendulum vibration absorber is parametrically excited by a primary structure have been studied extensively. One class of such systems are the autoparametric vibration absorbers, which give rise to interesting nonlinear responses, such as internal resonance [51, 52], amplitude- modulated response [53], and chaos [51]. The operation principle of AVAs is briefly explained as follows. When a primary structure is excited by force containing a frequency close to a natural frequency 𝜔1 of the structure, the structure will undergo resonance, thereby having large vibration. One way to mitigate the resonant response is to couple an AVA to the primary structure via quadratic nonlinearity and have the AVA’s natural frequency 𝜔2 tuned around half the natural frequency of the primary structure, i.e., 𝜔2 ≈ 𝜔1/2. When the force amplitude increases beyond a critical value, the primary structure’s vibration amplitude stops increasing or becomes saturated, and its vibration energy is transferred to the AVA. As a result, the AVA will have a resonance- like response, and then an electromechanical transducer is used to convert the resonance-like response into electricity. As such, structural vibration mitigation and vibration energy harvesting are achieved simultaneously. Because the AVA is internally excited by the primary structure and yet has a resonance-like response, the resonance is known as internal resonance. Note that when the natural frequencies satisfy 𝜔2 ≈ 𝜔1/2, they are known to be (nearly) commensurable in the literature [54], which is a necessary condition for internal resonances. If a heaving WEC employed an AVA as a power take-off (PTO) unit, it would have the potential to overcome the fundamental problem faced by integrating a floating platform and WECs. However, AVAs are not applicable to wave energy conversion because of two major shortcomings. First, because the natural frequencies 8 of AVA and heaving WEC must be commensurable (𝜔2 ≈ 𝜔1/2 ), AVA tends to be too long and heavy. For example, the natural frequency in the heave of the spar-WEC integration in [11] is about 0.1 Hz, so it would require an AVA of a natural frequency around 0.05 Hz to mitigate the resonant heave response. For the pendulum-type AVAs in [55], this suggests a pendulum length of 100 meters. Second, the saturation and energy transfer phenomenon of AVAs generally occur within a narrow range of excitation frequencies [52]. Ocean waves typically come in a spectrum of frequencies; that is, the dominant excitation frequency of ocean waves may change from time to time. Therefore, if this dominant frequency changed, AVAs would become detuned and therefore perform with low efficiency. Consequently, these two shortcomings call for a novel type of internal resonance vibration absorber that (a) can create a similar saturation phenomenon without having commensurable natural frequencies so that it can be lightweight and compact and (b) is responsive to a broad range of excitation frequencies or is broadband. Therefore, motivated by the autoparametric vibration absorbers and enormous wave energy potential promised by the successful integration of a floating platform and WECs, this work aims at studying the hydrodynamics and wave energy conversion capabilities of a spar-WEC integration where the WEC incorporates the inerter pendulum vibration absorber power-take off system (IPVA- PTO). The IPVA-PTO does not require commensurable natural frequencies and, therefore, can have a light and compact design. 1.4 Outline of the document The introductory chapter provides background, motivation, and the literature survey regarding the research statement. Chapter 2 talks about the integration of the IPVA with a single-degree- of-freedom system, including experimental analysis. Chapter 3 delves into the application of the IPVA-PTO system in ocean wave energy conversion along with experiments on the “dry” system and talks about the effect of the IPVA-PTO on the stability of the spar. Chapter 4 and 5 discuss the integration of IPVA and IPVA-PTO system with a spar-floater system for the purpose of hydrodynamic response suppression and energy conversion, respectively. Chapter 6 concludes the primary research findings, and avenues for future research are discussed. 9 CHAPTER 2 ANALYSIS OF THE IPVA SYSTEM: A STUDY IN VIBRATION SUPPRESSION 2.1 Overview To understand the behavior of the inerter pendulum vibration absorber (IPVA) system for ocean wave energy harvesting, we first analyze and perform experiments on the dry single-degree-of- freedom system by integrating the IPVA system into it. By dry, it means that the system is excited in a dry state, without any ocean wave effects on it. This is done to get an idea of how the dynamics of an oscillating system change when the IPVA is integrated with it and the energy harvesting potential of the pendulum. 2.2 Design of the Inerter Pendulum Vibration Absorber Both systems have two degrees of freedom, one associated with the pinion’s angular displacement (𝜃), and the other with the pendulum’s angular displacement relative to the pinion (𝜙). 2.2.1 Equations of motion Although different mechanisms are used, the working principle of both systems is identical. Therefore, their equations of motion are identical. Lagrange’s equations are used to derive the equations of motion. We derive the equations of motion for the rack-pinion system. First, the kinetic energy of the system is derived. We have the total kinetic energy 𝑇 = 𝑇𝑀 + 𝑇𝑐 + 𝑇𝑝 (2.1) where 𝑇𝑀, 𝑇𝑐 and 𝑇𝑝 is the kinetic energy of the primary mass 𝑀, carrier-pinion composite and the pendulum respectively. The coordinate system is defined in Fig. Fig. 2.1a. Therefore, the position vector of the pendulum r𝑝 can be obtained as 𝒓 𝒑 = (cid:0)𝑅𝑝 sin(𝜃) + 𝑟 sin(𝜃 + 𝜙)(cid:1) ˆ𝑰 + (cid:0)𝑅𝑝 cos(𝜃) + 𝑟 cos(𝜃 + 𝜙)(cid:1) ˆ𝑲. (2.2) where ˆ𝑰 is the direction of the motion of the primary mass, and ˆ𝑲 is the vertical direction. Differentiating it, we obtain (cid:164)𝒓 𝒑 = (cid:0)𝑅𝑝 cos(𝜃) (cid:164)𝜃 + 𝑟 cos(𝜃 + 𝜙) (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1)(cid:1) ˆ𝑰 + (cid:0)−𝑅𝑝 sin(𝜃) (cid:164)𝜃 − 𝑟 sin(𝜃 + 𝜙) (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1)(cid:1) ˆ𝑲.(2.3) 10 a SDOF system incorporating rack- pinion IPVA b CAD realization of rack-pinion IPVA system c SDOF system containing ball- screw IPVA system d CAD realization of the ball-screw IPVA system Figure 2.1 Rack-pinion and ball-screw IPVA systems: schematics and computer-aided (CAD) realizations Therefore, the kinetic energy of the pendulum is 𝑇𝑃 = 1 2 𝑚 (cid:164)𝒓 𝒑 · (cid:164)𝒓 𝒑 + 1 2 𝐽𝑝 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1) 2 , which simplifies to 𝑇𝑃 = 1 2 𝑚 (cid:16) 𝑝 (cid:164)𝜃2 + 𝑟 2 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1) 2 + 2𝑅𝑝𝑟 cos (𝜙) (cid:164)𝜃 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1) (cid:17) 𝑅2 + 𝐽𝑝 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1) 2 . 1 2 (2.4) (2.5) 𝑇𝑀 = 𝑇𝑝 = 1 2 1 2 𝑀 (cid:0)𝑅 (cid:164)𝜃(cid:1) 2 , 𝑇𝑐 = 𝐽 (cid:164)𝜃2, 𝐽𝑝 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1) 2 + 1 2 1 2 (cid:16) 𝑚 𝑝 (cid:164)𝜃2 + 𝑟 2 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1) 2 + 2𝑅𝑝𝑟 cos (𝜙) (cid:164)𝜃 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1)(cid:17) 𝑅2 (2.6) 11 are the kinetic energy of the structure, carrier, and pendulum. Here, 𝐽 is the moment of inertia of the carrier-pinion composite, 𝐽𝑝 is the moment of inertia of the pendulum with respect to its center of mass 𝑉 = 1 2 𝑘𝑥2 = 1 2 𝑘 𝑅2𝜃2. (2.7) To account for energy loss at the pivot point of the pendulum, a torsional viscous damping coefficient 𝑐 𝑝 is introduced. The virtual work done by the force 𝐹 = 𝐹0 sin (Ω𝑡), the damping torque in the pendulum, and the damping force in the primary mass can be derived as 𝐹𝛿𝑥, −𝑐 𝑝 (cid:164)𝜙𝛿 (cid:164)𝜙 and −𝑐 (cid:164)𝑥𝛿 (cid:164)𝑥, respectively, where 𝑐 𝑝 and 𝑐 are the torque damping coefficient in the pendulum and damping coefficient of the viscous damper between primary mass and ground respectively. Then the virtual work done by the force 𝐹, the damping torque (due to 𝑐 𝑝 and viscous damping 𝑐) are derived as 𝛿𝑊 = 𝐹 𝑅𝛿𝜃 − 𝑐 𝑝 (cid:164)𝜙𝛿𝜙 − 𝑐𝑅2 (cid:164)𝜃𝛿𝜃. (2.8) Therefore, the equations of motion of the system obtained using the Euler-Lagrange formulation are written as (cid:16) 𝑀 𝑅2 + 𝐽 + 𝑚𝑅2 𝑝 + 𝑚𝑟 2 + 2𝑚𝑅𝑝𝑟 cos 𝜙 (cid:17) (cid:165)𝜃 + (cid:16) 𝑚𝑟 2 + 𝑚𝑅𝑝𝑟 cos 𝜙 + 𝐽𝑝 (cid:17) (cid:165)𝜙 + 𝑐𝑅2 (cid:164)𝜃 +𝑘 𝑅2𝜃 − 2𝑚𝑅𝑝𝑟 (cid:164)𝜙 (cid:164)𝜃 sin 𝜙 − 𝑚𝑅𝑝𝑟 (cid:164)𝜙2 sin 𝜙 = 𝐹0 sin (Ω𝑡) 𝑅, (cid:16) 𝑚𝑟 2 + 𝐽𝑝 (cid:17) (cid:165)𝜙 + 𝑚 (cid:16) 𝑟 2 + 𝑅𝑝𝑟 cos 𝜙 (cid:17) (cid:165)𝜃 + 𝑐 𝑝 (cid:164)𝜙 + 𝑚𝑅𝑝𝑟 (cid:164)𝜃2 sin 𝜙 = 0 (2.9) For initial analysis, it is assumed that the pendulum is made of a point mass such that its moment of inertia with respect to the pivot point is much larger than the pendulum’s moment of inertia with respect to its center of mass, i.e., 𝑚𝑟 2 >> 𝐽𝑝. Furthermore, 𝐽 = (𝑚 𝑝 + 𝑚𝑐)𝑅2 𝑔, where 𝑚 𝑝 and 𝑚𝑐 are the pinion mass and carrier mass, respectively, and 𝑅𝑔 is the radius of gyration. As the primary mass 𝑀 is much larger than the sum of 𝑚 𝑝 and 𝑚𝑐 and as the pinion radius 𝑅 and the radius of gyration 𝑅𝑔 have the same order of magnitude, it is assumed that 𝑀 𝑅2 >> 𝐽. Without loss of generality, 𝐽 and 𝐽𝑝 are neglected. We rescale the time and convert (2.9) into a dimensionless form for further analysis using the following parameters, 12 𝜇𝑟 = 𝜉 = 𝑚𝑅2 𝑝 𝑀 𝑅2 𝑐 2𝜔0𝑀 , 𝜔0 = √︂ 𝑘 𝑀 , 𝜔 = Ω 𝜔0 , 𝜏 = 𝜔0𝑡, 𝜂 = , 𝜉𝑝 = 𝑐 𝑝 2𝜔0𝑀 𝑅2 , 𝑓 = 𝐹0 𝑀 𝑅𝜔2 0 , 𝑟 𝑅𝑝 𝑑 () 𝑑𝜏 , ()′ = . (2.10) Denote x = [𝜃, 𝜙]𝑇 and f = [ 𝑓 sin 𝜔𝜏, 0]𝑇 . The dimensionless equations of motion are obtained as Mx′′ + Cx′ + Kx + g (x, x′, x′′) = f (2.11) where M =        1 + 𝜇𝑟 (cid:0)1 + 𝜂2(cid:1) 𝜇𝑟𝜂2 𝜇𝑟𝜂2 0 2𝜉        , C =        (2𝜃′′ + 𝜙′′) cos 𝜙 − 𝜙′ (2𝜃′ + 𝜙′) sin 𝜙 , K = 𝜇𝑟𝜂2 2𝜉 𝑝        0 𝜃′′ cos 𝜙 + 𝜃′2 sin 𝜙 g (x, x′, x′′) = 𝜇𝑟𝜂        ,               1 0 0 0 .        (2.12) It is worth noting that the strength of the nonlinear inertial terms g (x, x′, x′′) is proportional to 𝜇𝑟 and 𝜂. The moment of inertia ratio 𝜇𝑟 can be readily magnified by adjusting the ratio 𝑅 𝑝 𝑅 , thereby creating strong nonlinear inertial effects with a small pendulum mass. For example, for a mass ratio 𝑚 10 leads to 𝜇𝑟 = 0.3, indicating that the inertia effect is magnified 𝑀 = 3%, a ratio 𝑅 𝑝 𝑅 = √ by a factor of ten. Furthermore, the pendulum length ratio 𝜂 is proportional to the length of the pendulum. Therefore, a long pendulum leads to strong nonlinear inertial effects. 2.3 Internal resonance of IPVA According to the studies on autoparametric resonance, internal resonance plays an essential role in transferring the kinetic energy of a primary structure to the pendulum vibration absorber [56, 57, 58]. As will be demonstrated in Sec. 2.6, when internal resonance occurs to the IPVA, a similar energy transfer phenomenon is observed, resulting in vibration mitigation of the primary structure. In this section, we will determine the conditions for which internal resonance will occur to the IPVA. To this end, we use the harmonic balance method to determine the parametric instability of the system. 13 2.3.1 Harmonic balance method By virtue of the harmonic balance method, periodic solutions of the system are assumed to take the following form 𝜃 𝑝 (𝜏) = 𝑃 ∑︁ (cid:16) 𝑝=1 Θ𝑐 𝑝 cos (cid:17) (cid:16) 𝑝𝜔𝜏 𝜈 + Θ𝑠 𝑝 sin (cid:16) 𝑝𝜔𝜏 𝜈 (cid:17)(cid:17) , 𝜙 𝑝 (𝜏) = Φ0 + 𝑃 ∑︁ (cid:16) 𝑝=1 Φ𝑐 𝑝 cos (cid:17) (cid:16) 𝑝𝜔𝜏 𝜈 + Φ𝑠 𝑝 sin (cid:17)(cid:17) (cid:16) 𝑝𝜔𝜏 𝜈 (2.13) where Θ𝑝, Φ𝑝 and Φ0 are unknown Fourier coefficients to be determined. Note that 𝜈 ∈ N accounts for subharmonics. Furthermore, a constant Φ0 is included to consider asymmetric oscillation of the centrifugal pendulum [59]. Denote by x𝑝 = (cid:2)𝜃 𝑝, 𝜙𝑝(cid:3)𝑇 the vector of assumed periodic solutions. After substituting (2.13) into the equations of motion (2.11), we obtain the following residue term R(𝜏) = Mx′′ 𝑝 + Cx′ 𝑝 + Kx𝑝 − g (cid:16) x𝑝, x′ 𝑝, x′′ 𝑝 (cid:17) − f. (2.14) To obtain an expression relating the Fourier coefficients, a Galerkin procedure [60] is used to project (2.14) on the orthogonal trigonometric basis, yielding 2𝑃 + 1 nonlinear algebraic equations ℎ0 ( ˆx) = ℎ𝑐 𝑝 ( ˆx) = ∫ 2 𝜋𝜈 𝜔 0 ∫ 2 𝜋𝜈 𝜔 0 R(𝜏)𝑑𝜏 = 0, ℎ𝑠 𝑝 ( ˆx) = ∫ 2 𝜋𝜈 𝜔 0 R(𝜏) sin (cid:17) (cid:16) 𝑝𝜔𝜏 𝜈 𝑑𝜏 = 0, R(𝜏) cos (cid:17) (cid:16) 𝑝𝜔𝜏 𝜈 𝑑𝜏 = 0 (2.15) with ˆx = (cid:2)𝚯𝑇 , 𝚽𝑇 , Φ0 (cid:16) 𝑝, x′′ 𝑝 Note that g x𝑝, x′ (cid:3)𝑇 (cid:17) , Θ = (cid:2)Θ𝑐 1 , · · · , Θ𝑐 𝑃, Θ𝑠 1 , · · · , Θ𝑠 𝑃 (cid:3)𝑇 and 𝚽 = (cid:2)Φ𝑐 1 will result in composite trigonometric terms like cos , · · · , Φ𝑐 (cid:16) , · · · , Φ𝑠 𝑃 𝑃, Φ𝑠 1 Φ𝑠 𝑝 sin ( 𝑝𝜔𝜏/𝜈) . . (cid:17) (cid:3)𝑇 These terms can be expanded using the Jacobi–Anger expansion, namely, an infinite series of products of Bessel functions and trigonometric functions [61] ( see Appendix B for the expansion formulas). For the current study, the Jacobi-Anger expansion is truncated at Bessel functions of order up to a third to capture the necessary nonlinear effects. We solve (2.15) for the Fourier coefficients using Newton-Raphson method. Substitution of the Fourier coefficients into (2.13) will lead to the periodic solutions. The stability of the periodic solutions will be determined in the next section. 14 2.3.2 Stability To determine the stability of the periodic solutions, small perturbations are introduced to (2.13) as follows 𝜃 (𝜏) = 𝜃 𝑝 (𝜏) + 𝛿𝜃 (𝜏) and 𝜙(𝜏) = 𝜙 𝑝 (𝜏) + 𝛿𝜙 (𝜏) (2.16) where |𝛿𝜃 (𝜏)| << 1 and |𝛿𝜙 (𝜏)| << 1. Denote the vector of small perturbations by 𝜹 = (cid:2)𝛿𝜃, 𝛿𝜙(cid:3)𝑇 . Substitution of (2.16) into (2.11) and linearization with respect to 𝜃 𝑝 (𝜏) and 𝜙 𝑝 (𝜏) yield (cid:18) M + (cid:19) 𝜕g 𝜕x′′ (cid:18) C + 𝜹′′ + (cid:19) 𝜕g 𝜕x′ (cid:18) K + 𝜹′ + (cid:19) 𝜕g 𝜕x 𝜹 = 0 (2.17) where the Jacobian matrices 𝜕g/𝜕x′′, 𝜕g/𝜕x′, and 𝜕g/𝜕x are evaluated at x = x𝑝, x′ = x′ 𝑝, and x′′ = x′′ 𝑝 wherever appropriate. Note that the Jacobian matrices are periodic functions of period 𝑇 = 2𝜋𝜈/𝜔, e.g., 𝜕g/𝜕x (𝜏) = 𝜕g/𝜕x (𝜏 + 𝑇). Because (2.17) have periodic coefficients, one can use Floquet theory to determine the stability [62]. To this end, (2.17) is transformed into the state- space form and numerically integrated using Matlab’s ode45 (based on an explicit Runge-Kutta integration method) over one period 𝑇 to obtain the fundamental matrix.] To ensure accuracy, the absolute and relative tolerance was taken to be 10−9. If any eigenvalue of the fundamental matrix has a magnitude greater than unity, the periodic solutions are unstable. Although (2.17) can be used to determine the stability of arbitrary periodic solutions, it is a numerical approach; hence, it is hard to understand how internal resonance occurs to the IPVA. To gain physical insights, we also use a semi-analytical approach to determine the stability. To that end, we apply a multiple-scale approach to (2.17) as follows. Because a compact and lightweight design of the IPVA system is preferred in practical applications, we consider 𝜇𝑟 << 1. Assuming that the parameters 𝜇𝑟, 𝜉, 𝜉 𝑝 are small quantities, we set 𝜇𝑟 = 𝜖 ˆ𝜇𝑟, 𝜉 = 𝜖 ˆ𝜉, 𝜉𝑝 = 𝜖 ˆ𝜉 𝑝 and introduce the following asymptotic expansions 𝛿𝜃,𝜙 (𝜏) = 𝛿(0) 𝜃,𝜙 (𝜏0, 𝜏1, · · · ) + 𝜖 𝛿(1) 𝜃,𝜙 (𝜏0, 𝜏1, · · · ) + · · · 𝜏𝑘 = 𝜖 𝑘 𝜏, 𝑘 = 0, 1, · · · 𝑑 𝑑𝜏 + · · · + 𝜖 = 𝜕 𝜕𝜏1 𝜕 𝜕𝜏0 15 (2.18) where |𝜖 | << 1 is a small bookkeeping parameter. After substituting (2.18) into (2.17) and collecting terms that will lead to parametric instabilities, the equation obtained in order 𝑂 (𝜖 0) is 𝜕2𝛿(0) 𝜃 𝜕𝜏2 0 + 𝛿(0) 𝜃 = 0, 𝜕2𝛿(0) 𝜙 𝜕𝜏2 0 2 ˆ𝜉 𝑝 ˆ𝜇𝑟𝜂2 𝜕𝛿(0) 𝜙 𝜕𝜏0 + 𝐴 (cid:0)x𝑝(cid:1) 𝜂 + 𝛿(0) 𝜙 = 0 where 𝐴 (cid:0)x𝑝(cid:1) = cos(𝜙 𝑝) (cid:17) 2 (cid:16) 𝜃′ 𝑝 − sin(𝜙 𝑝)𝜃′′ 𝑝 (2.19) (2.20) is a periodic coefficient of period 𝑇. It is worth noting that the first equation in (2.19) shows that 𝛿(0) 𝜃 are stable harmonic functions. Therefore, the stability of the periodic solutions is determined by the second equation in (2.19). When the nonlinearity is weak, periodic solutions (2.13) are dominated by primary harmonics, i.e., 𝑃 = 1 and 𝜈 = 1. As (2.19) is derived by the multiple-scale approach, it is accurate when the nonlinearity is weak. Therefore, in addition to (2.17) (Floquet theory), (2.19) is used to determine the boundary of parametric instability for periodic solutions of primary harmonics, which will explain how internal resonance occurs to the IPVA. Thus, 𝑃 = 1 and 𝜈 = 1 are substituted in (2.13) to obtain 𝜃 𝑝 = Θ𝑐 1 sin (𝜔𝜏). After substituting these expressions of Θ and Φ into (2.19) and expanding in terms of Bessel functions 1 sin (𝜔𝜏) and 𝜙 𝑝 = Φ0 + Φ𝑐 1 cos (𝜔𝜏) + Φ𝑠 1 cos (𝜔𝜏) + Θ𝑠 up to the third order, we arrive at a damped Mathieu equation as follows 𝜕2𝛿(0) 𝜙 𝜕𝜏2 0 2𝜉 𝑝 𝜇𝑟𝜂2 𝜕𝛿(0) 𝜙 𝜕𝜏0 + + 𝑢 ( ˆx) 𝜂 𝛿(0) 𝜙 + 𝑣 ( ˆx) 𝜂 cos(𝜔𝜏 − 𝛾)𝛿(0) 𝜙 = 0 (2.21) 16 where 𝑢 ( ˆx) = 𝜔2Θ1c (Φ0) 2 {Θ1 [𝐽0(Φ1) + 𝐽2(Φ1)s (2𝛼)] + 2𝐽1(Φ1)s (𝛼)}, 𝑣 ( ˆx)2 = [𝑏0𝑏1 + 2𝑏2 (𝑏1 + 𝑏3)] c (𝛼) + (cid:16) 2𝑏0𝑏2 + 2𝑏1𝑏3 − 𝑏2 1 (cid:17) c (2𝛼) + (2𝑏0𝑏3 − 𝑏1𝑏2) c (3𝛼) − 𝑏1𝑏3c (4𝛼) + 𝑏2 0 + 5𝑏2 1 4 + 𝑏2 2 + 𝑏2 3 , 𝑏0 = s(Φ0)Θ1𝜔2𝐽0(Φ1), 𝑏1 = −s(Φ0)Θ2 1 𝜔2𝐽1(Φ1), 𝑏2 = −s(Φ0)Θ1𝜔2𝐽2(Φ1), 𝑏3 = −s(Φ0)Θ2 1 √︂ (cid:16) (cid:17) 2 𝜔2𝐽3(Φ1), (cid:17) 2 (cid:16) Θ1 = Θ𝑐 + 1 1Θ𝑐 (cid:18) Φ𝑠 1Θ𝑐 Φ𝑐 Θ𝑠 1 − Φ𝑐 + Φ𝑠 , Φ1 = 1Θ𝑠 (cid:19) 1Θ𝑠 1 1 1 1 𝛼 = tan−1 √︂ (cid:16) (cid:17) 2 + (cid:16) Φ𝑠 1 (cid:17) 2 , Φ𝑐 1 (2.22) where c (·) = cos (·), s (·) = sin (·), and 𝐽𝑛 (·) denotes Bessel functions of the first kind of order 𝑛. Note that the detail of phase angle 𝛾 is not provided because 𝛾 is irrelevant to stability. There are two things worth noting in (2.21) and (2.22). First, 𝑏0 = 𝑏1 = 𝑏2 = 𝑏3 = 0 or 𝑣 ( ˆx) = 0 when Φ0 = 0. Because 𝑣 ( ˆx) is the magnitude of parametric excitation, no parametric instabilities can occur when 𝑣 ( ˆx) = 0. In other words, nonzero asymmetric oscillation, i.e., Φ0 ≠ 0, is a necessary condition for parametric instabilities. Second, the linear stiffness term 𝑢 ( ˆx) is composed of nonlinear inertial coupling induced by the carrier motion Θ1. As the nonlinear inertial coupling results in linear stiffness per se, the pendulum can have internal resonance without having any linear stiffness. Compared to the autoparametric vibration absorbers, which would need low linear stiffness to tune their natural frequency around half the natural frequency of the primary structure [56, 57, 58], the nonlinear inertial coupling of the IPVA enables compact designs. To determine the boundary of parametric instability, (2.21) is transformed to the standard form of Mathieu equation [51] 𝜕2Ψ 𝜕𝑤2 + 𝑝 ( ˆx) Ψ − 2𝑞 ( ˆx) cos(2𝑤)Ψ = 0 (2.23) 17 where (cid:32) Ψ = 𝛿(0) 𝜙 exp 𝑝 ( ˆx) = 4𝑢 ( ˆx) 𝜂𝜔2 − (cid:33) 2 ˆ𝜉 𝑝𝑤 ˆ𝜇𝑟𝜂2𝜔 4 ˆ𝜉2 𝑝 ˆ𝜇2 𝑟 𝜂4𝜔2 , 2𝑤 = 𝜔𝜏 − 𝛾 , 𝑞 ( ˆx) = − 2𝑣 ( ˆx) 𝜂𝜔2 (2.24) The boundary of parametric instability for (2.23) corresponds to the transition curves in the 𝑝 − 𝑞 plane [63]. Because we seek the boundary that occurs with low force magnitudes, we compute the transition curve with the lowest 𝑝 and 𝑞 values. Mathematically, this transition curve is expressed as 𝑝 = A1 (𝑞), where A1 (𝑞) are the characteristic values for even Mathieu functions with characteristic exponent 1 and parameter 𝑞 [63]. In this paper, A1 (𝑞) is computed by the “MathieuCharacteristicA” function of Wolfram Mathematica 11.3. Note that 𝑝, 𝑞 are functions of 𝑓 and 𝜔. Therefore, 𝑝 = A1 (𝑞) is solved with (2.15) simultaneously to yield the transition curves in the 𝑓 − 𝜔 plane. 2.3.3 Pitchfork bifurcation As mentioned in Sec. 2.3.2, nonzero asymmetric oscillation, i.e., Φ0 ≠ 0, is necessary to induce parametric instabilities. To determine when it occurs, 𝜃 𝑝 = Θ𝑐 1 sin (𝜔𝜏) and 𝜙 𝑝 = Φ0 + Φ𝑐 1 sin (𝜔𝜏) are substituted into (2.17) to solve for stable periodic solutions with Φ0 ≠ 0. We use the pendulum length ratio 𝜂 as the bifurcation parameter to obtain 1 cos (𝜔𝜏) + Θ𝑠 1 cos (𝜔𝜏) + Φ𝑠 a bifurcation diagram that shows the parameter space wherein Φ0 ≠ 0 will occur. To track the bifurcation points with varying 𝜂, a bifurcation tracking algorithm based on arc-length continuation is used with (2.17); see Appendix A for the detail. Figure 2.2 shows a bifurcation diagram of Φ0 with varying 𝜂. Three bifurcation branches were obtained using three sets of initial conditions (one corresponding to each branch, namely, the lower, middle, and upper). It can be observed that Φ0 undergoes a supercritical pitchfork bifurcation at a critical value of 𝜂. After this critical value of 𝜂, Φ0 ≠ 0 and parametric instabilities become possible. For the rest of this paper, we will only explore parametric instabilities with the parameters that lead to Φ0 ≠ 0. 18 2.3.4 Period-doubling bifurcation Within the parameter space wherein Φ0 ≠ 0 exist, the boundary of parametric instability is computed in the 𝑓 −𝜔 plane. To find an initial bifurcation point for the bifurcation tracking algorithm described in Appendix A, 𝜔 = 0.8 is set, and (2.17) is repeatedly used to compute the Floquet multipliers as 𝑓 decreases until the maximum magnitude of the Floquet multipliers becomes unity. Afterward, the bifurcation tracking algorithm will generate the boundary as described in Appendix A. To verify the boundary is indeed of parametric instability, the Mathieu equation (2.23) is used to generate the transition curve as described in Sec. 2.3.2. The boundary and transition curve for a set of parameters are shown in Fig. 2.3. Although the transition curve underestimates the boundary, they are in qualitative agreement. Specifically, the discrepancy between the two curves increases as the force magnitude 𝑓 increases. Because the perturbation method predicted the transition curve, it is reasonable that it is more accurate for small force magnitudes. Thus, the comparison verifies the claim that the boundary indicates parametric instability. To gain more insight, the Floquet exponents corresponding to a few points on the boundary are computed and found equal to ±𝑖𝜋/𝑇, where √ −1. According to [64], this indicates period-doubling bifurcation. Since periodic-doubling 𝑖 = bifurcation is a co-dimension one bifurcation, it is a curve in a parameter plane [65]. Therefore, the parametric instability boundary is, in fact, a boundary of period-doubling bifurcation. When this bifurcation occurs, the pendulum oscillation will have subharmonics of 𝜔/2, i.e., 𝜈 = 2 in (2.13). It is worth noting that the autoparametric vibration absorbers also have a similar bifurcation behavior; that is, subharmonics of half excitation frequency induced by parametric instabilities [56, 57, 58]. Within the parameter space wherein subharmonics of 𝜔/2 exist, the stability of the subharmonics can further be investigated. Preliminary investigations indicate the presence of another period- doubling bifurcation, implying that subharmonics of quarter frequency will appear. Therefore, it is hypothesized that there exists a cascade of period-doubling bifurcations in the 𝑓 − 𝜔 space, which eventually leads to chaotic motions of the system. However, determining the boundary of this additional period-doubling bifurcation is out of this paper’s scope. Figure 2.3 is a bifurcation diagram that shows the parameter space for qualitatively different 19 Figure 2.2 Pitchfork bifurcation of 𝜙0 for the following parameter values 𝑓 = 0.007, 𝜇𝑟 = 0.2, 𝜔 = 0.9, 𝜉 = 𝜉 𝑝 = 0.005 Figure 2.3 Parametric instability boundary for the following parameter values 𝜂 = 0.3, 𝜇𝑟 = 0.25, 𝜉 = 0.005, 𝜉𝑝 = 0.005 solutions defined by the instability boundary. By locating the parameters in Fig. 2.3, the qualitative behavior of the corresponding solutions can be predicted. For example, ×2 resides in the parameter space just above the boundary. Accordingly, periodic solutions of primary harmonics and subharmonics of excitation frequency are predicted at ×2. Next, by direct numerical integration, we verify the predictions by Fig. 2.3. 2.4 Numerical demonstration To verify the bifurcation analysis in Sec. 2.3.4, numerical integration (Matlab’s ode45) is used to obtain the solutions of (2.11) at three representative points in Fig. 2.3 (denoted by markers “×” followed by numbers, e.g., ×2). Among these three points, points ×1 and ×2 lead to periodic solutions, whereas point ×3 leads to non-periodic solutions. The fast Fourier transform (FFT) 20 0.10.20.30.40.50.60.70.8-1.5-1-0.500.511.50stableunstable(2.17)(cid:10)(2.23)(cid:10) of the periodic solutions is computed to reveal the frequency components, which are shown in Fig. 2.4a, Fig. 2.5a and Fig. 2.6a. On the other hand, a time series of the solutions are presented to show the dynamical behaviors, shown in Fig. 2.4b, Fig. 2.5b, Fig. 2.6b, Fig. 2.4c, Fig. 2.5c and Fig. 2.6c. Note that the frequencies ˆ𝜔 of the FFT are normalized with respect to the excitation frequency. It follows that primary harmonics correspond to components at ˆ𝜔 = 1, subharmonics of half excitation frequency correspond to components at ˆ𝜔 = 0.5, etc. There are several things worth noting in Fig. 2.3, Fig. 2.4, Fig. 2.5 and Fig. 2.6. First, the prediction at point ×1 is in good agreement with the numerical solutions. As shown in Fig. 2.3, point ×1 is below the instability boundary. It is expected that primary harmonics dominate the periodic solutions. This prediction is verified by Fig. 2.4, which shows that the periodic solutions at point ×1 have the largest components at ˆ𝜔 = 1, corresponding to primary harmonics. Furthermore, in Fig. 2.3, as we increase the value of 𝑓 and reach point ×2, the primary harmonics undergo a period-doubling bifurcation. As a result, subharmonics of half excitation frequency should arise. As shown in Fig. 2.5, subharmonics of half excitation frequency exist, which verifies the prediction in Fig. 2.3. Second, the parameters at ×3 lead to strong non-periodic solutions composed of both oscillation and intermittent rotations of the pendulum, as shown in Fig. 2.6c and Fig. 2.7c. Similar non-periodic solutions are also observed in autoparametric resonance vibration absorbers [58]. a FFT at point 1 b Time series of 𝜃 at point 1 c Time series of 𝜙 at point 1 Figure 2.4 FFT and time series of periodic solutions at point 1 in Fig. 2.3 In addition to FFT, the Poincaré sections demonstrate the period-doubling bifurcations predicted 21 00.511.5210-510-410-310-210-1100101102FFTFFT of FFT of 22002250230023502400-0.2-0.15-0.1-0.0500.050.10.150.22820284028602880290029202940296029803000-23.85-23.8-23.75-23.7-23.65-23.6-23.55-23.5-23.45 a FFT at point 2 b Time series of 𝜃 at point 2 c Time series of 𝜙 at point 2 Figure 2.5 FFT and time series of periodic solutions at point 2 in Fig. 2.3 a FFT at point 3 b Time series of 𝜃 at point 3 c Time series of 𝜙 at point 3 Figure 2.6 FFT and time series of non-periodic solutions at point 3 in Fig. 2.3 by Fig. 2.3. The Poincaré sections are computed by the Hénon trick [66], which are defined as 𝑃𝑛 (x0) = x𝑛 (𝜏0 + 2𝑛𝜋/𝜔; x𝑛−1, 𝜏0) , 𝑛 = 1, 2, · · · (2.25) where x𝑛−1 and x𝑛 are the solutions of the system (2.11) which pass through the Poincaré section at time 𝜏 = 𝜏0 + 2(𝑛 − 1)𝜋/𝜔 and 𝜏 = 𝜏0 + 2𝑛𝜋/𝜔, respectively. Successively, the points x0, x1 = 𝑃 (x0), x2 = 𝑃2 (x0) , · · · correspond to the intersection of the trajectory x (𝜏; x0, 𝜏0) with the sections at 𝜏 = 𝜏0, 𝜏0 + 2𝜋/𝜔, 𝜏0 + 4𝜋/𝜔, · · · , respectively. To demonstrate the period-doubling bifurcation, Poincaré sections corresponding to point ×1 and ×2 are plotted in Fig. 2.7. As shown in Fig. 2.7a, ×1 leads to a fixed point on the Poincaré section, corresponding to period-1 solutions. Figure 2.7b, on the other hand, shows two fixed points corresponding to period-2 solutions. Therefore, it is clear that the system has undergone a period-doubling bifurcation when moving from ×1 to ×2. 22 00.511.5210-510-410-310-210-1100101102FFTFFT of FFT of 33503400345035003550-0.2-0.15-0.1-0.0500.050.10.150.232003250330033503400345035003550-18-17.8-17.6-17.4-17.2-1700.511.5210-510-410-310-210-1100101FFTFFT of FFT of 1950200020502100215022002250230023502400-0.3-0.2-0.100.10.20.3500100015002000250030003500-8-6-4-202468 a Poincaré section at point 1 b Poincaré section at point 2 c Poincaré section at point 3 Figure 2.7 Poincare section at points 1, 2, and 3 as marked in Fig. 2.3. The value of 𝜙 is normalized between 0 and 2𝜋 Figure 2.8 Parametric instability boundary for 𝜇𝑟 = 0.4, 𝜉 = 0.005, 𝜉 𝑝 = 0.01 for various values of 𝜂 Figure 2.9 Parametric instability boundary for 𝜂 = 0.4, 𝜉 = 0.005, 𝜉 𝑝 = 0.01 for various values of 𝜇𝑟 23 44.555.5-1.5-1-0.500.511.50.911.11.21.31.41.51.61.7-0.09-0.085-0.08-0.075-0.07-0.06501234567-0.6-0.4-0.200.20.40.60.80.80.850.90.9510.010.020.030.040.050.060.070.080.090.1f = 0.3 = 0.4 = 0.50.80.820.840.860.880.90.920.940.960.9810.010.020.030.040.050.060.070.080.090.1fr = 0.2r = 0.3r = 0.4 Figure 2.10 Parametric instability boundary for 𝜂 = 0.3, 𝜇𝑟 = 0.2, 𝜉𝑝 = 0.005 for various values of 𝜉 Figure 2.11 Parametric instability boundary for 𝜂 = 0.4, 𝜇𝑟 = 0.4, 𝜉 = 0.005 for various values of 𝜉 𝑝 2.5 Parametric studies In this section, we analyze the effect of the parameters on the instability boundary. We consider four parameters in (2.10), namely 𝜇𝑟, 𝜂, 𝜉, and 𝜉 𝑝. It can be seen that these parameters can be varied independently of each other. Therefore, we will observe the effect by varying one parameter while keeping the other parameters constant. We start by increasing the value of 𝜂 while keeping the others constant. From Fig. 2.8, it can be observed that increasing 𝜂 does not make any significant change in the lowest 𝑓 value for parametric instability to occur, which corresponds to the vertex of the boundaries. That means that value of 𝜂 should not influence the energy transfer capabilities of the system by a lot. However, a minimum threshold value of 𝜂 is required for period-doubling bifurcation, as discussed in subsection 2.3.3. We next vary 𝜇𝑟. From Fig. 2.9, it can be observed 24 0.80.850.90.9510.010.020.030.040.050.06f = 0.01 = 0.03 = 0.050.810.830.850.870.890.910.930.950.970.990.0050.0150.0250.0350.0450.0550.0650.0750.0850.0950.105fp = 0.005p = 0.01p = 0.015 that the required values of 𝑓 to attain parametric instability decrease as 𝜇𝑟 increases. This can be attributed to the fact that the inertia supplied by the pendulum vibration absorber increases as the mass amplification factor 𝜇𝑟 increases. The value of 𝜇𝑟 can be controlled by changing the ratio 𝑅𝑝/𝑅, which can be adjusted by changing the carrier radius (𝑅𝑝). While keeping the other parameters constant, we vary the viscous damping ratio 𝜉 and observe its effects. From Fig. 2.10, it can be seen that the requirement of 𝑓 to achieve nonlinear energy transfer increases with an increase in the viscous damping. Similarly, we vary the 𝜉 𝑝 while keeping the other parameters constant. We see that the values of 𝑓 required to achieve parametric instability increase as the viscous damping increases. The observations on the effects of both viscous damping match well with the effect of viscous damping on parametric instability — the larger the viscous damping, the larger force it takes to cause parametric instability [63]. Last but not least, the parameter 𝜇𝑟 significantly influences the instability boundary, as demonstrated in Sec. 2.5. As seen in Fig. 2.9, a larger 𝜇𝑟 not only leads to lower force magnitudes required to cause internal resonance but a wider frequency bandwidth of internal resonance, which is beneficial in terms of vibration mitigation. The parameter 𝜇𝑟 can be readily increased by changing the ratio of 𝑅𝑝/𝑅 without incurring considerable weight to the system, which is attributed to the mass amplification effect of the inerter. 2.6 Discussion At the beginning of this study, it was proposed that a nonlinear energy transfer phenomenon similar to autoparametric resonance occurs when parametric instability occurs. To demonstrate this, we compare the proposed system with two systems, a linear benchmark and an autoparametric vibration absorber with a parametrically excited pendulum [57]. The linear system here is characterized by locking the pendulum at its initial position (𝜙 = 0), effectively removing all the nonlinearities in the system. By setting 𝜙 = 𝜙′ = 𝜙′′ = 0 in (2.11), the equation of motion of 25 the linear system is written as (cid:104) 1 + 𝜇𝑟 (cid:16) 1 + 𝜂2 + 2𝜂 (cid:17)(cid:105) 𝜃 ′′ + 2𝜉𝜃 ′ + 2𝜉 𝑝𝜃 ′ + 𝜃 = 𝑓 sin 𝜔𝜏. Using 𝜃 = 1 2Θ𝑒𝑖𝜔𝜏 + 𝑐.𝑐., the equation can be solved to obtain Θ = (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 𝑓 2𝑖𝜔 (cid:0)𝜉 + 𝜉 𝑝(cid:1) + 1 − 𝜔2 (cid:2)1 + 𝜇𝑟 (cid:0)1 + 𝜂2 + 2𝜂(cid:1) (cid:3) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) . (2.26) (2.27) We computed the root-mean-square (RMS) of the IPVA system and compared it with the linear system. The comparison is shown in Fig. 2.12. The response from the 2401-th to 3000-th cycle was used to compute the RMS to eliminate transient effects. The IPVA parameters used in Fig. 2.12a and Fig. 2.12b correspond to Fig. 2.3 and Fig. 2.11, respectively. a Comparison of linear system ( 𝑓0 = 𝑓0 0.007) with IPVA for different values for 𝜂 = 0.3, 𝜇𝑟 = 0.25, 𝜉 = 0.005, 𝜉𝑝 = 0.005 b Comparison of linear system ( 𝑓0 = 0.025) with IPVA for different 𝑓0 values for 𝜂 = 0.4, 𝜇𝑟 = 0.4, 𝜉 = 0.005, 𝜉𝑝 = 0.015 Figure 2.12 Comparison of IPVA with the linear system for various parameters Several things are worth noting in Fig. 2.12. First, it is shown that the response of the primary structure flattens for a range of excitation frequencies. In comparison with the response of the linear system, the IPVA shows significant vibration suppression with the flattening region. For example, as shown in Fig. 2.12b, 𝜃 for 𝑓0 = 0.025 flattens for 𝜔 ∈ [0.81, 0.87]. In comparison with Fig. 2.3 and Fig. 2.11, it is clear that the flattening occurs when the system is within the parametric 26 0.850.90.9500.10.20.30.40.50.6 = x/R (rad)f0 = 0.007 (IPVA)f0 = 0.011 (IPVA) f0 = 0.007 (Linear)0.810.830.850.870.890.910.930.950.970.9900.511.522.5 = x/R (rad)f0 = 0.025 (IPVA)f0 = 0.035 (IPVA)f0 = 0.025 (Linear)0.820.840.860.880.230.2350.240.2450.25 instability boundary. For example, as shown in Fig. 2.11, when 𝜉 𝑝 = 0.015 and 𝑓0 = 0.025, the system is within the parametric instability boundary for 𝜔 ∈ [0.81, 0.87]. This observation agrees with Fig. 2.12b. Second, within the flattening region, the response of the primary structure barely increases despite an increase in the force magnitude, suggesting a saturation phenomenon similar to autoparametric vibration absorbers [67] and nonlinear vibration absorbers with quadratic nonlinearities [68]. Note that the response of the IPVA system in Fig. 2.12b is non-periodic for 𝑓 = 0.035 within a range of 𝜔; thus, different initial conditions may lead to different rms responses. To examine the effect of initial conditions, ten RMS responses were computed and plotted at 6 discrete 𝜔 values, each corresponding to a different initial condition vector (cid:2)𝜃, 𝜙, (cid:164)𝜃, (cid:164)𝜙(cid:3)𝑇 that was randomly chosen from a standard normal distribution with zero mean and a unit standard deviation; see the inset in Fig. 2.12b. Figure 2.13 The autoparametric vibration absorber [57] Next, we compare IPVA with the autoparametric vibration absorber shown in Fig. 2.13. Because the autoparametric system oscillates in the vertical direction, the ball-screw IPVA is considered hereinafter. The equation of motion of the autoparametric system is written as [57] (𝑀 + 𝑚) (cid:165)𝑥 + 𝑐 (cid:164)𝑥 + 𝑘𝑥 + 𝑚𝑙 (cid:16) (cid:165)𝜙 sin 𝜙 + (cid:164)𝜙2 cos 𝜙 (cid:17) = 𝐹0 sin (Ω𝑡) 𝑚𝑙2 (cid:165)𝜙 + 𝑐𝑎 (cid:164)𝜙 + (𝑚𝑔𝑙 + 𝑚𝑙 (cid:165)𝑥) sin 𝜙 = 0 (2.28) where 𝑥 and 𝜙 represent the primary structure and pendulum angular displacement. 𝑀, 𝑘, and 𝑐 are the mass, stiffness, and viscous damping coefficient of the primary structure, respectively, and 𝑚 and 𝑙 are the pendulum mass and length, respectively. Furthermore, a viscous damping coefficient 𝑐𝑎 is introduced to account for energy loss at the pivot point of the pendulum. 27 According to [57], when the pendulum’s natural frequency is tuned around half of the natural frequency of the primary structure, the system shows autoparametric resonance for a certain set of parameters when excited harmonically. This autoparametric resonance results in energy transfer from the primary structure to the pendulum, thereby achieving vibration suppression of the primary structure. Because the IPVA system and the autoparametric system achieve vibration suppression in a similar way, the latter is an ideal benchmark system for comparison. For a fair comparison, the primary structure parameters, excitation force magnitude, and pendulum mass are kept identical in both systems. Specifically, primary mass 𝑀 = 5 kg, natural frequency of the primary structure 𝜔0 = √︁𝑘/𝑀 = 4𝜋rad/s, pendulum mass 𝑚 = 0.5 kg, force magnitude 𝐹0 = 0.491 N, and 𝜉 = 0.005 or 𝜉 = 0.01. Note that two values of 𝜉 are considered to examine the performance of the IPVA when the damping ratio of the primary structure changes. The remaining parameters pertaining to the autoparametric system are 𝑐𝑎 = 8.68 × 10−5 N·m·s and 𝑙 = 12.42cm, which were taken from ref.[57]. Specifically, the pendulum length was chosen to achieve autoparametric resonance, and the pendulum damping coefficient was determined from 𝜉𝑎 = 𝑐𝑎/(2𝑚𝑙2𝜔 𝑝) = 0.05, where 𝜔 𝑝 = √︁𝑔/𝑙 is the natural frequency of the pendulum. The remaining parameters pertaining to the IPVA system are 𝑅, 𝑅𝑝, and 𝑟. Three sets of 𝑅, 𝑅𝑝, and 𝑟 were chosen: (a) 𝑅 = 2.49 cm, 𝑅𝑝 = 4.97 cm, and 𝑟 = 1.99 cm; (b) 𝑅 = 1.78 cm, 𝑅𝑝 = 3.55 cm, and 𝑟 = 1.42 cm; and (c) 𝑅 = 2.07 cm, 𝑅𝑝 = 4.14 cm, and 𝑟 = 1.66 cm. These three sets are labeled as IPVA (a), IPVA (b), and IPVA (c), respectively, in Fig. 2.14 and Fig. 2.15. These three sets all lead to 𝜇𝑟 = 0.4 and 𝜂 = 0.4, and lead to 𝑓 = 0.025, 𝑓 = 0.035, and 𝑓 = 0.030, respectively. This way, the dependence of the IPVA on different values of 𝑓 will be examined. The RMS response of the IPVA system and autoparametric system were computed using the same direct numerical integration scheme with the same settings that were used to generate Fig. 2.12. The effects of initial conditions were examined for both systems when their responses were observed to be non-periodic using the same method used to obtain Fig. 2.12b. There are several things worth noting in Fig. 2.14. First, it can be seen that the response curve of the primary structure displacement 𝑥 for both IPVA (a) and IPVA (b) flatten for a range of 𝜔, 28 of for system frequency a Comparison 𝜃 IPVA and response of autoparametric with 𝜂 = 0.4, 𝜇𝑟 = 0.4, 𝜉 = 0.005, and for IPVA 𝜉 𝑝 = 0.015, whereas 𝜉𝑎 = 0.05 for autoparametric system b Comparison of the frequency response of pendulum’s velocity for IPVA and autoparametric system with 𝜂 = 0.4, 𝜇𝑟 = 0.4, 𝜉 = 0.005, and for IPVA 𝜉 𝑝 = 0.015, whereas 𝜉𝑎 = 0.05 for autoparametric system Figure 2.14 Comparison of IPVA and the autoparametric vibration absorber for a set of parameters of for frequency a Comparison 𝑥 response of IPVA and autoparametric vibration absorber with 𝜂 = 0.4, 𝜇𝑟 = 0.4, 𝜉 = 0.01. For IPVA, 𝜉 𝑝 = 0.015, and for autoparametric vibration absorber, 𝜉𝑎 = 0.05 b Comparison of the velocity of the pendulum for IPVA and autoparametric vibration absorber with 𝜂 = 0.4, 𝜇𝑟 = 0.4, 𝜉 = 0.01. For IPVA, 𝜉 𝑝 = 0.015, and for autoparametric vibration absorber, 𝜉𝑎 = 0.05 Figure 2.15 Comparison of IPVA and the autoparametric vibration absorber for a set of parameters demonstrating the energy transfer phenomenon for two different sets of parameters. Specifically, IPVA (b) has a more compact design (smaller 𝑅, 𝑅𝑝, and 𝑟) and shows better performance. 29 0.80.850.90.95100.0050.010.0150.020.0250.03x (m)IPVA (a)IPVA (b)Autopara. in [57]0.830.850.877.588.5910-30.80.850.90.95100.511.522.533.544.55IPVA (a)IPVA (b)Autopara. in [57]0.80.850.90.95100.0050.010.015x (m)IPVA (c)Autopara. in [57]0.80.850.90.95100.511.522.533.544.5IPVAAutopara. in [57] Although the autoparametric system shows similar vibration suppression, both IPVA (a) and IPVA (b) outperform it. Second, let us examine the pendulum response in Fig. 2.14b. As seen, for both IPVA (a) and IPVA (b), the pendulum response significantly increases within the 𝜔 range of parametric instability (𝜔 ∈ [0.81, 0.87] for IPVA (a) and 𝜔 ∈ [0.80, 0.89] for IPVA (b)), indicating that the kinetic energy of the primary structure transfers to the pendulum, resulting in the response flattening observed in Fig. 2.14a. It is noteworthy that both IPVA (a) and IPVA (b) have a larger pendulum angular velocity than the autoparametric system. Similarly, Fig. 2.15a and Fig. 2.15b show the comparison of IPVA (c) with the autoparametric system for 𝜉 = 0.01. As can be observed, IPVA (c) outperforms the autoparametric system in terms of vibration suppression. Furthermore, it also leads to a larger pendulum angular velocity. In addition to better vibration suppression, the IPVA system has two other advantages compared to the autoparametric system. First, it generates higher pendulum angular velocities, as shown in Fig. 2.14b and Fig. 2.15b. Kecik and Boroweic [69] proposed an energy harvesting autoparametric system where they installed an electromagnetic generator at the pendulum pivot point to convert the pendulum angular motion into electricity. As the larger angular velocity, the larger electricity can be generated, and the larger angular velocity in the IPVA system may lead to better performance in terms of energy harvesting, which remains to be explored in the future. Second, the IPVA system leads to a more compact design. The largest length in the IPVA system is the sum 𝑅𝑝 + 𝑟 of the carrier radius and pendulum length. IPVA (a) has 𝑅𝑝 + 𝑟 = 6.69 cm, the maximum among the three. On the contrary, the autoparametric system requires a long pendulum (𝑙 = 12.42 cm) as it needs this length to tune the natural frequency. Next, the experimental setup is discussed to verify the anaylsis performed. 2.7 Experimental analysis To verify the analysis performed, experiments consisting of a IPVA integrated with a single- degree-of-freedom system are conducted. As has been demonstrated in previous sections (Sec. 2.3.4 and Sec. 2.4), the system shows an internal resonance phenomenon. Thus, the aim of the experiments is to verify the energy transfer between the primary system and the pendulum vibration 30 absorber and the secondary resonance phenomenon of the IPVA (harmonics of frequency 𝜔/2). An experimental setup is created by integrating a single-degree-of-freedom system with the IPVA system as shown in Fig. 2.16. Table 2.1 shows the description of labels for various components of the experimental setup. The top plate, marked by plate 𝐴, supports the primary mass. The base plate system contains three plates, plate 𝐵, 𝐶 and 𝐷. Plate 𝐷 is directly connected to the shaker, whereas the nut of the ball-screw system connects plate attached to the plate 𝐵. Note that the plates 𝐵, 𝐶 and 𝐷 are rigidly connected. A coupler connects the ball-screw system with the carrier which hosts the pendulum. The eight springs are connected between the primary mass system and the base plate system. Therefore, due to the excitation of the base plate system, the top plate system starts moving. The relative motion between the top plate and the base plate system drives the screw (as the nut is fixed to the base plate system), i.e. 𝜃 = 𝑥 − 𝑦 𝑅 where 𝑥 is the motion of plate 𝐴 and 𝑦 is the motion of plate 𝐷. The pendulum is free to rotate and is connected to the carrier by a ball bearing. A shaker (APS 113) excites the base plate by controlling the motion of the plate 𝐷 and is driven by a spectral analyzer and controller (Spider 80x) using an amplifier. There is an accelerometer connected to the base plate for closed-loop control of the shaker, and an accelerometer is connected to the top plate to monitor the motion of the top plate system. An encoder is mounted on the pendulum shaft to measure the motion of the pendulum using a microcontroller (US Digital). 2.7.1 Equations of motion for the experimental setup The equations of motion of the experimental system can be directly derived from (Sec. 2.2.1), converting the forced excitation into base excitation. To obtain this equation of motion, we assume perfect transmission by the ball-screw system. Thus, the resultant equation of motion is given by ¯M (cid:165)X + ¯C (cid:164)X + ¯KX + ¯G(X, (cid:164)X, (cid:165)X) = ¯F (2.29) 31 Figure 2.16 Experimental setup 1 + 𝜇𝑟 + 𝜇𝑟𝜂2 + 2𝜂𝜇𝑟 cos(𝜙) + 𝜇𝑏𝑠𝑐 + 𝜇 𝑝 𝜇 𝑝 + 𝜇𝑟𝜂2 + 𝜇𝑟𝜂 cos(𝜙) 𝜇 𝑝 + 𝜇𝑟𝜂2 + 𝜇𝑟𝜂 cos(𝜙) 𝜇 𝑝 + 𝜇𝑟𝜂2 ,        2𝜉 0 0 2𝜉 𝑝               , ¯G = 𝜇𝑟𝜂 (cid:169) (cid:173) (cid:173) (cid:171) −2𝜃′𝜙′ sin(𝜙) − 𝜙′2 sin(𝜙) 𝜃′2 sin(𝜙) , (cid:170) (cid:174) (cid:174) (cid:172) (2.30) where ¯M = ¯K =               ¯F = (cid:169) (cid:173) (cid:173) (cid:171) , ¯C = 1 0 0 0        −𝑦′′ −𝑡 𝑓 𝑝sign( (cid:164)𝜙) (cid:170) (cid:174) (cid:174) (cid:172) with X = [𝜃, 𝜙]𝑇 and 𝑡 𝑓 𝑝 = 𝑇 𝑓 𝑀 𝑅2 , where 𝑇 𝑓 is the friction between the shaft of the pendulum and its bearing. Here, 𝑥 is the motion of the top plate, 𝑦 is the motion of the shaker, 𝜃 is the rotation of the carrier, 𝑅𝑝 is the location of the shaft of the pendulum with respect to the center of the 𝑐 𝑝 𝑀𝜔0𝑅2 , 𝑐 𝑝 being the damping between the pendulum’s shaft and ball bearing. 𝑀 is the mass of the top-plate (the ball-screw system, 𝜉 𝑝 is the damping coefficient of the pendulum defined by 𝜉 𝑝 = primary structure), 𝐽𝑏𝑠𝑐 includes the inertia of carrier, ball-screw, coupler and encoder and 𝐽𝑝 is the 32 ABCDGEFHIJ213 Label 1 2 3 A B C D E F G H I J Description Accelerometer signal conditioner Spectral analyzer and controller Shaker signal amplifier Top plate Ball-screw mounting plate Middle plate Lower plate Top added mass Ball-screw system Pendulum Carrier Encoder data acquisition box Shaker Table 2.1 Labels depicting the different parts of the experimental setup moment of inertia of the pendulum system with respect to its center of mass. We assume theoretical values for all the inertia, masses, transmission ratio, springs, and lengths in the system. Therefore, the values of 𝑚, 𝑅𝑝, 𝑀, 𝑅, 𝐽𝑝, 𝐽𝑏𝑠𝑐, 𝑘, 𝑟 and 𝜔0 are assumed to be known; see Table 2.2. This leaves 𝑐 and 𝑐 𝑝 to be experimentally identified, which is discussed next. 2.7.2 Obtaining the mechanical damping 𝑐 The mechanical damping value 𝑐 is obtained by removing the pendulum from the experimental setup and obtaining its frequency response function (FRF). The measured FRF and theoretical FRF are correlated with each other to obtain the mechanical damping value. The linear system without the pendulum has the following equation of motion where ˆ𝑀 (cid:165)𝑥 + 𝑘𝑥 + 𝑐 (cid:164)𝑥 = 𝑀 (cid:165)𝑦 + 𝑐 (cid:164)𝑦 + 𝑘 𝑦 ˆ𝑀 = 𝑀 (1 + 𝜇𝑏𝑠𝑐) , 𝜇𝑏𝑠𝑐 = 𝐽𝑏𝑠𝑐 𝑀 𝑅2 (2.31) (2.32) Using 𝑥 = 𝑋 𝑒𝑖𝜔𝑡 and 𝑦 = 𝑌 𝜔2 𝜔2 plate and the top plate respectively, and defining 𝑟 = 𝜔 𝜔0 , we obtain 𝑒𝑖𝜔𝑡 where 𝑋 and 𝑌 denote the displacement amplitude of the base 33 𝑋 𝑌 = 𝑏𝑠𝑐𝑟 4 − 2𝜇𝑏𝑠𝑐𝑟 2 + 4𝜉2𝑟 2 + 1 𝜇2 (𝜇𝑏𝑠𝑐 + 1)2 𝑟 4 − 2𝑟 2 (cid:0)𝜇𝑏𝑠𝑐 − 2𝜉2 + 1(cid:1) + 1 (2.33) Matlab’s “fit” function is used to fit the data and obtain the value of 𝜇𝑏𝑠𝑐, 𝜔0, and 𝑐. We found that the value of 𝜇𝑏𝑠𝑐 and 𝜔0 were close to the theoretical values, and therefore theoretical values of 𝜇𝑏𝑠𝑐 and 𝜔0 are chosen for the analysis of experiments. 2.7.3 Obtaining the pendulum damping 𝑐 𝑝 To obtain the pendulum damping, the pendulum and carrier were removed from the ball-screw system and fixed it rigidly to the ground. The pendulum was given some initial rotational velocity and the free rotational behavior of the pendulum was observed as it came to a stop. It was observed that damping alone is insufficient to capture the effects of the energy dissipation in the pendulum, so a friction term was considered in the analysis of the interaction between the shaft of the pendulum and its bearing. The rotation of the pendulum can be modeled using the following equation of motion (cid:16) 𝐽𝑝 + 𝑚𝑟 2(cid:17) (cid:165)𝜙 + 𝑐 𝑝 (cid:164)𝜙 + 𝑇 𝑓 = 0. This equation (2.34) is integrated with respect to time to obtain (cid:16) 𝐽𝑝 + 𝑚𝑟 2(cid:17) (cid:0) (cid:164)𝜙 𝑓 − (cid:164)𝜙𝑖(cid:1) + 𝑐 𝑝 (cid:0)𝜙 𝑓 − 𝜙𝑖(cid:1) + 𝑇 𝑓 Δ𝑡 = 0. (2.34) (2.35) where Δ𝑡 is the time it took for pendulum to stop, subscripts 𝑖 and 𝑓 on 𝜙 and (cid:164)𝜙 denote the initial and final values of 𝜙 and (cid:164)𝜙 respectively. Clearly, (cid:164)𝜙 𝑓 = 0 as the pendulum comes to a full stop. The equation can be solved to obtain 𝑐 𝑝. The experimental parameters are tabulated in Table 2.2 for ready reference. 2.7.4 Experimental verification In this section, the simulation model is verified with the experiments. First, the stability boundary is obtained and shown in Fig. 2.17 using the theoretical and fitted parameters shown in Table 2.2. Since it is known that there should be vibration suppression above the stability 34 Parameter 𝜉 𝑝 𝜉 𝜇𝑟 𝜂 𝜔0 𝜇𝑔 𝑇 𝑓 𝜇 𝑝 𝜇𝑏𝑠𝑐 Value 0.00418 0.242 4.072 0.371 96.98 rad/s 0.0029 0.0014N 0.388 8.294 Table 2.2 Experimental parameters for the IPVA system Figure 2.17 Numerical stability boundary of energy transfer corresponding to the experimental setup boundary, an acceleration value of 0.2 𝑔 is used to run the experiments, which is well above the stability boundary for a range of frequencies. Before observing the experimental results and comparing them to the linear system, a few points need to be stated: 1. The ball-screw system is assumed to be 100% efficient, with no loss in transmission 2. The effects of accelerometer cables on the system are neglected, and the encoder cable is assumed to only contribute to mechanical damping and moment of inertia of the ball-screw system 3. The springs are modeled to be linear for the analytical analysis 35 3.544.5Frequency (Hz)0.150.20.250.30.350.4Acceleration (g) a Time series of the pendulum motion b Measured FFT of the pendulum motion Figure 2.18 Time series and measured FFT of the pendulum motion 2.7.5 Verification of the internal resonance Previously it was defined that the secondary solutions are those solutions which contain harmonics of frequency 𝜔 2 when the excitation frequency is 𝜔. As has been reported in [70], internal resonance is required for the energy transfer to occur between the primary system and the pendulum vibration absorber and thus suppressing the vibrations of the primary system. To observe the secondary resonance of the system, experiments are performed on the IPVA system by fixing a frequency and increasing base acceleration value such that the secondary resonance is found. For this particular case, an excitation acceleration amplitude of 0.2 𝑔 and an excitation frequency of 4.16 Hz led to the secondary resonance, which means that the response of the system at this excitation frequency and acceleration amplitude will have harmonics of frequency half the excitation frequency (2.08 Hz), along with the excitation frequency (4.16 Hz). As can be seen from Fig. 2.18a, the motion of the pendulum measured by the encoder shows harmonics of frequency 𝜔 2 (≈ 2.08Hz) as evident from the time series of the pendulum’s motion and the FFT shown in Fig. 2.18b. Moreover, it can be observed that harmonics of frequency 3𝜔 2 and 2𝜔 also show up in the motion of the pendulum. 36 47.647.84848.248.448.648.84949.249.4t (s)-0.2-0.100.10.20.30.40.5012345678Frequency (Hz)10-210-1100101FFT of X 6.18954Y 1.18961X 8.2659Y 1.31328X 4.13295Y 15.3356X 2.05659Y 8.70668 2.7.6 Comparison with the linear system Next, the IPVA system is compared with the linear system to benchmark its vibration suppression capabilities. For this, the base excitation acceleration is fixed at 0.2 𝑔. The linear system is chosen by fixing the system’s pendulum at an angle such that the natural frequency of the linear system is equal to the resonant frequency of the IPVA system. Ten frequency points are chosen for which the experiments are run on the IPVA system to get the root-mean-squared (RMS) values of the amplitude of ball-screw rotation (𝜃) and pendulum’s motion (𝜙), and obtain the mean and standard deviations for both the rotations. Similar experiments are performed for the linear system by converting the frequency response function between the acceleration of the top plate 𝐴 and the input base acceleration from the shaker’s to RMS data for 𝜃 for the same ten frequencies as the IPVA system. Sixteen frequency response function were experimentally calculated to obtain these RMS 𝜃 values for the linear system. These RMS values are plotted against each other for comparison in Fig. 2.19a. The “Linear LB” and “Linear UB” labels show lower bound and the upper bound of the RMS values of 𝜃 for the linear system as calculated from the frequency response function, and error bars are plotted for the ball-screw motion (𝜃) for the IPVA system. Fig. 2.19b shows the motion of the pendulum observed experimentally. It is clear that the RMS values for 𝜃 do not follow a resonance like behavior for the range of frequencies shows, whereas the linear system does. This verifies that there is nonlinear energy transfer between the ball-screw (the primary mass motion), and the pendulum as cross-verified from the energy pumping in the pendulum, see Fig. 2.19b. Next, the case when the base excitation acceleration is 0.21 𝑔 is taken. Shown in Fig. 2.20a and Fig. 2.20b, we see results similar to that of 0.2 𝑔 base acceleration. However, the pendulum has more energy in this case compared to 0.2 𝑔 excitation due to a saturation-like phenomenon, where the motion of the primary system (𝜃) does not increase significantly even though the excitation acceleration is increased. This has also been observed in a numerical study by the authors [70]. 2.8 Conclusion This study analyzes the IPVA system proposed in [71] with a focus on vibration suppression of a linear oscillator subject to single harmonic excitation. It is shown that for a given excitation 37 Comparison between a experimental 𝜃 RMS value for linear and the IPVA system b Experimental frequency response RMS of the pendulum for the IPVA system Figure 2.19 Frequency response of 𝜃 and 𝜙 for the IPVA and the linear system (where applicable) at 0.2 𝑔 base excitation acceleration Comparison a between experimental 𝜃 RMS value for linear and the IPVA system b Experimental frequency response RMS of the pendulum for the IPVA system Figure 2.20 Frequency response of 𝜃 and 𝜙 for the IPVA and the linear system (where applicable) at 0.21 𝑔 base excitation acceleration 38 44.14.24.34.44.5Frequency (Hz)0.20.250.30.350.40.450.50.55 (rad)Linear LBLinear UBExp. IPVA44.054.14.154.24.254.34.354.44.454.5Frequency (Hz)3.544.555.566.5744.14.24.34.44.50.20.250.30.350.40.450.50.55Linear LBLinear UBExp. IPVA44.14.24.34.44.54.6Frequency [Hz]2468101214Mean phidot [rad/s] force magnitude, the pendulum parameters can be chosen such that internal resonance occurs to the pendulum vibration absorber for a specific range of excitation frequencies. It is also shown that a pitchfork bifurcation and period-doubling bifurcation of the pendulum response are necessary and sufficient conditions for internal resonance. Furthermore, when internal resonance occurs, the kinetic energy of the linear oscillator transfers to the pendulum, resulting in vibration suppression of the linear oscillator. A saturation phenomenon similar to autoparametric vibration absorbers and nonlinear vibration absorbers is observed in the IPVA system; that is, the response of the linear oscillator saturates despite the increase in the force magnitude. Meanwhile, the increased energy due to the increase in the force magnitude seems to transfer to the pendulum, resulting in an increased pendulum response. Furthermore, the system is compared to the autoparametric vibration absorber. It outperforms the autoparametric vibration absorber in vibration absorption and energy harvesting capabilities. Finally, experiments integrating the IPVA system with a single-degree-of-freedom spring mass system are performed. The secondary resonance which corresponds to the crossing of 1:2 internal resonance boundary is experimentally shown in the IPVA system. Furthermore, it is observed for two different acceleration values, that above the period-doubling bifurcation boundary (1:2 internal resonance boundary), the motion of the primary system was suppressed, and its energy transferred to the pendulum vibration absorber. 39 CHAPTER 3 OCEAN WAVE ENERGY CONVERSION USING THE IPVA-PTO SYSTEM WITH A SPAR 3.1 Overview The inerter pendulum vibration absorber (IPVA) system with an electromagnetic power take-off, known as the IPVA-PTO system, is integrated with a spar to study its ocean wave energy conversion efficacy. We analyze the system with hydrodynamic effects and perform experiments on the “dry” single-degree-of-freedom system by integrating the IPVA system into it. By dry, it means that the system is excited in a dry state, without any ocean wave effects on it. Both numerical and experimental analysis is compared with a linear benchmark. 3.2 The IPVA-PTO system Figure 3.1 shows the IPVA-PTO system mounted between a spar and a fixed frame. The spar (primary structure) is floating in the water, with hydrodynamic stiffness 𝑘 in the direction of heaving (𝑥). The system consists of a lead screw and nut mounted between a fixed reference and the spar such that the heaving displacement is converted into the angular displacement 𝜃 through 𝑥 = 𝑅𝜃, where 𝑅 = 𝐿/2𝜋, where 𝐿 is the screw lead. The carrier is fixed to the screw and has the same angular displacement. The pendulum pivots on the point of the carrier, which is located at a distance of 𝑅𝑝 from the carrier center. The pendulum length is 𝑟 and has the angular displacement 𝜙 with respect to the screw. The planetary gear system combines the pendulum and screw motion into one angular motion (cid:164)𝜃 − (cid:164)𝜙 input to the generator for electricity. A generalized force 𝐹𝑔 is the force on the spar assumed to contain a single harmonic for the current analysis due to the wave-structure interaction. Its derivation is discussed further in Sec. 3.2.2. The inductance of the generator is neglected because typical wave frequencies are too low to be significant [72]. This means the 𝑐𝑒 (cid:0) (cid:164)𝜃 − (cid:164)𝜙(cid:1) 2, where 𝑐𝑒 is the electrical electrical power generated by the generator can be written as 1 2 damping in the system. For this study, the 1:100 sparD model in [73] is considered. Note that the heaving plate and the mooring lines are ignored, and the spar is assumed to have only heaving motion. The mass of the 40 a Schematic of the system b CAD realization of the IPVA-PTO system c The benchmark linear system Figure 3.1 The IPVA-PTO system integrated with a spar: schematic and CAD realization, and the benchmark linear system 41 ∇∇𝑥𝑐Ocean bed𝑀db2azxIPVA-PTOpendulum mass𝜃planet gear𝜙𝜓𝜃generator𝜓=θ−𝜙: sun gear & generator angular motionsun gear drives generatorscrewyzcarrierscrewnutgeneratorplanet gearsun geargear box𝜃𝜙𝜃: screw & carrier angular motion𝜙: pendulum angular motion Anut connected to the fixed referencehousing fixed to spar centrifugal pendulumfixed referenceNutScrewGenerator GearsCarrierPendulumTo ball screw spar 𝑀 is given as 18.83 kg, the draft 𝑑 is taken to be 0.96 m, and the diameter of the spar 2𝑎 is 0.16m. The depth of the water is ℎ = 3 m, and the height of the spar 𝑑 + 𝑏 is 1.06 m, which makes the freeboard 𝑏 = 0.1 m. The wave amplitude is 6 cm per [73]. The density of water is assumed to be 𝜌 = 1025 kg/m3 and acceleration due to gravity 𝑔 = 9.81 m/s2. 3.2.1 Equations of motion We use Lagrange’s equations to derive the equations of motion for the IPVA-PTO. The procedure outlined below is similar to the procedure explained in [70]. First, the total kinetic energy of the system is derived as where 𝑇 = 𝑇𝑀 + 𝑇𝑐 + 𝑇𝑝 + 𝑇𝑔 (3.1) 𝑇𝑀 = 𝑇𝑝 = + 𝑇𝑟 = 1 2 1 2 1 2 1 2 𝑀 (cid:0)𝑅 (cid:164)𝜃(cid:1) 2 , 𝑇𝑐 = 𝐽 (cid:164)𝜃2, 1 2 𝐽𝑝 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1) 2 𝑚 (cid:16) 𝑝 (cid:164)𝜃2 + 𝑟 2 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1) 2 + 2𝑅𝑝𝑟 cos (𝜙) (cid:164)𝜃 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1)(cid:17) 𝑅2 𝐽𝑟 (cid:0) (cid:164)𝜃 − (cid:164)𝜙(cid:1) 2 (3.2) are the kinetic energy of the spar, carrier, pendulum, and generator, respectively. Here, 𝐽 is the moment of inertia of the screw and carrier, 𝐽𝑝 is the moment of inertia of the pendulum with respect to its center of mass, and 𝐽𝑟 is the moment of inertia of the rotor of the generator. Furthermore, the potential energy is written as 𝑉 = 1 2 𝑘𝑥2 = 1 2 𝑘 𝑅2𝜃2. (3.3) To account for energy loss due to the screw motion, a viscous damping coefficient 𝑐 is introduced. The virtual work done by the force 𝐹𝑔, the electrical damping torque in the rotor, and the viscous damping force can be derived as 𝐹𝑔𝛿𝑥, −𝑐𝑒 (cid:0) (cid:164)𝜃 − (cid:164)𝜙(cid:1) 𝛿 (𝜃 − 𝜙) and −𝑐 (cid:164)𝑥𝛿𝑥, respectively, where 𝑐𝑒 and 𝑐 are the torsional electrical damping coefficient of the generator [72] and damping coefficient of the viscous damper, respectively. Then the virtual work done by the force 𝐹𝑔, the damping torque 42 (due to 𝑐𝑒 and viscous damping 𝑐) are derived as 𝛿𝑊 = 𝐹𝑔 𝑅𝛿𝜃 − 𝑐𝑒 (cid:0) (cid:164)𝜃 − (cid:164)𝜙(cid:1) 𝛿 (𝜃 − 𝜙) − 𝑐𝑅2 (cid:164)𝜃𝛿𝜃. (3.4) The pendulum is assumed to be a point mass; thus, 𝐽𝑝 is neglected. Furthermore, the moment of inertia of the screw and carrier is assumed to be small compared to the moment of inertia of other objects in the system; thus, 𝐽 is neglected. Therefore, the equations of motion of the system obtained using Lagrange’s equations are written as (cid:16) 𝑀 𝑅2 + 𝐽𝑟 + 𝑚𝑅2 (cid:16) + 𝑚𝑟 2 + 𝑚𝑅𝑝𝑟 cos 𝜙 𝑝 + 𝑚𝑟 2 + 2𝑚𝑅𝑝𝑟 cos 𝜙 (cid:17) (cid:165)𝜃 (cid:17) (cid:165)𝜙 + 𝑐𝑅2 (cid:164)𝜃 + 𝑐𝑒 (cid:0) (cid:164)𝜃 − (cid:164)𝜙(cid:1) +𝑘 𝑅2𝜃 − 2𝑚𝑅𝑝𝑟 (cid:164)𝜙 (cid:164)𝜃 sin 𝜙 − 𝑚𝑅𝑝𝑟 (cid:164)𝜙2 sin 𝜙 = 𝐹𝑔 𝑅, 𝑚𝑟 2 (cid:165)𝜙 + 𝑚 (cid:16) 𝑟 2 + 𝑅𝑝𝑟 cos 𝜙 (cid:17) (cid:165)𝜃 + 𝑐𝑒 (cid:0) (cid:164)𝜙 − (cid:164)𝜃(cid:1) +𝑚𝑅𝑝𝑟 (cid:164)𝜃2 sin 𝜙 = 0 (3.5) We rescale the time and convert (3.5) into a dimensionless form for further analysis using the following parameters, , 𝜔0 = √︂ 𝑘 𝑀 , 𝜔 = Ω 𝜔0 , 𝜏 = 𝜔0𝑡, 𝜂 = , 𝜉 = 𝑐 2𝜔0𝑀 , 𝜉𝑒 = 𝑐𝑒 2𝜔0𝑀 𝑅2 , 𝑓𝑔 = , 𝑟 𝑅𝑝 𝐹𝑔 𝑀 𝑅𝜔2 0 𝜇𝑟 = 𝜇𝑔 = ()′ = 𝑚𝑅2 𝑝 𝑀 𝑅2 𝐽𝑟 𝑀 𝑅2 𝑑 () 𝑑𝜏 . , (3.6) Let us define x = [𝜃, 𝜙]𝑇 and f𝑔 = (cid:2) 𝑓𝑔, 0(cid:3)𝑇 . The dimensionless equations of motion are obtained as Mx′′ + Cx′ + Kx + g (x, x′, x′′) = f𝑔 (3.7) 43 where M = C =               𝜇𝑔 + 1 + 𝜇𝑟 (cid:0)1 + 𝜂2(cid:1) 𝜇𝑟𝜂2 − 𝜇𝑔 𝜇𝑟𝜂2 − 𝜇𝑔 𝜇𝑔 + 𝜇𝑟𝜂2 2 (𝜉 + 𝜉𝑒) −2𝜉𝑒 −2𝜉𝑒 2𝜉𝑒 , K =        1 0 0 0        , ,               g (x, x′, x′′) = 𝜇𝑟𝜂 (2𝜃′′ + 𝜙′′) cos 𝜙 − 𝜙′ (2𝜃′ + 𝜙′) sin 𝜙 𝜃′′ cos 𝜙 + 𝜃′2 sin 𝜙               (3.8) 3.2.2 Hydrodynamic coefficients We use the linear wave theory to determine the hydrodynamic coefficients of the spar subject to regular incident waves, which assumes that the fluid is inviscid and irrotational. Based on the theory, the hydrodynamic force on the spar consists of three components: Froude-Krylov force, diffraction force, and radiation force. The first two correspond to the incident wave field without and with the spar, respectively, and the last one is due to the oscillation of the spar. The Froude-Krylov and diffraction force together give rise to the excitation force while the radiation force gives rise to the added mass and radiation damping [74, 26], which can be represented by the well-known Cummins’ equation [75]: 𝐹𝑔 = −𝐴∞ (cid:165)𝑥 − ∫ ∞ 𝜎=0 𝑘 𝑅 (𝜎) (cid:164)𝑥 (𝑡 − 𝜎) 𝑑𝜎 + 𝛾𝐹𝑒 (𝑡) (3.9) where 𝛾 is the wave amplitude, 𝐹𝑒 is the excitation force per wave amplitude, and the radiation impulse response kernel, 𝑘 𝑅 (𝜎), and the radiation infinite-frequency added mas, 𝐴∞, are related to the radiation frequency-dependent damping and added mass 𝐵𝑅 (Ω) and 𝐴𝑅 (Ω), through Ogilvie’s relations [76] 𝐵𝑅 (Ω) = 𝐴𝑅 (Ω) = 𝐴∞ − 1 Ω ∫ ∞ 𝜎=0 ∫ ∞ 𝜎=0 𝑘 𝑅 (𝜎) cos (Ω𝜎) 𝑑𝜎 𝑘 𝑅 (𝜎) sin (Ω𝜎) 𝑑𝜎 (3.10) 44 and 𝐴∞ = lim Ω→∞ 𝐴𝑅 (Ω) (3.11) After substituting (3.10) and (3.11) into (3.7) and using (3.6), the equation of motion now becomes (M + A∞) x′′ + Cx′ + +g (x, x′, x′′) = f ∫ ∞ 𝑠=0 K𝑅 (𝑠) x′ (𝜏 − 𝑠) 𝑑𝑠 + Kx (3.12) (3.13) In (3.12), A∞ =        𝐴∞ 𝑀 0 0 0        , K𝑅 (𝑠) =        𝜅𝑅 (𝑠) 0 0      0   , f =    𝑓𝑒 (𝜏) 0    (cid:16) (cid:17) where 𝑠 = 𝜔0𝜎 is the normalized time, 𝜅𝑅 (𝑠) = 𝑘 𝑅 (𝑠/𝜔0) / impulse response kernel, and 𝑓𝑒 = 𝛾𝐹𝑒/(𝑀 𝑅𝜔2 𝑀𝜔2 0 0) is the dimensionless excitation force. is the normalized radiation 3.2.3 Determination of hydrodynamic coefficients The hydrodynamic coefficients, Froude-Krylov and diffraction forces are determined using Ansys AQWA. A convergence test is conducted to match the published results in [77]. A particular case with 𝑎 = 0.2 m, ℎ = 1 m, and 𝑑 = 0.25 m is chosen. Fig. 3.2 shows the comparison between added mass and radiation damping between Ansys and the published results. Note that 𝑚0 is the solution to the equation 𝑚0 tanh (𝑚0) = Ω2ℎ 𝑔 (3.14) From Fig. 3.2a and Fig. 3.2b, it can be seen that the Ansys solutions and the published results are in close agreement with each other. After verifying the Ansys model, we adopt the same setting to simulate the sparD system with the parameters provided in Sec. 3.2. Fig. 3.3 shows the mesh of the system used for analysis and Fig. 3.4 shows the added mass, radiation damping, and excitation force obtained using Ansys for sparD. 45 g a Added mass comparison b Radiation damping comparison Figure 3.2 Comparison of added mass and radiation damping between Yeung et al. [77] and Ansys simulation for 𝑎 = 0.2m, ℎ = 1m, 𝑑 = 0.25m Figure 3.3 Ansys AQWA model for calculation of hydrodynamic coefficients a Added mass b Radiation damping c Diffraction and Froude- Krylov force Figure 3.4 Added mass, radiation damping, and diffraction and Froude-Krylov force for sparD 46 012345m0a33.544.555.5AnsysYeung012345m0a00.20.40.60.8AnsysYeung051015201.061.081.11.121.141.161.18AR()0510152000.020.040.060.080.1BR()024681012141618020406080100120140160180200Fe 3.3 Internal resonance and bifurcation boundaries Internal resonance is essential for nonlinear vibration absorbers to absorb the vibration energy of a primary structure they are attached to [51, 57, 58, 70]. To investigate the internal resonance of the IPVA-PTO system, we will first determine the conditions required to have internal resonance using the harmonic balance method in conjunction with the modified alternating frequency/time approach discussed in Section 3.3.2. The convolution terms are handled separately, as discussed in Section 3.3.3. Specifically, 1:2 internal resonance will be the focus of the analysis as it has been shown to achieve energy absorption from an oscillating structure seeking vibration mitigation in the authors’ past study [70]. This study uses the harmonic balance method to investigate internal resonance. 3.3.1 Harmonic balance analysis In a fashion similar to [70], we assume the periodic solutions of the system take the form 𝜃 𝑝 (𝜏) = Θ0 + 𝑃 ∑︁ (cid:104) 𝜙 𝑝 (𝜏) = Φ0 + 𝑝=1 𝑃 ∑︁ 𝑝=1 Θ𝑐 𝑝 cos (cid:17) (cid:16) 𝑝𝜔𝜏 𝜈 + Θ𝑠 𝑝 sin (cid:16) 𝑝𝜔𝜏 𝜈 (cid:17)(cid:105) , (cid:104) Φ𝑐 𝑝 cos (cid:17) (cid:16) 𝑝𝜔𝜏 𝜈 + Φ𝑠 𝑝 sin (cid:17)(cid:105) (cid:16) 𝑝𝜔𝜏 𝜈 (3.15) where Θ𝑝, Φ𝑝, Θ0 and Φ0 are unknown Fourier coefficients to be determined. As has been reported previously [70], Φ0 is included to consider asymmetric pendulum oscillations [59], whereas Θ0 is introduced to expedite the derivation. It is found that Θ0 = 0 as the spar oscillates symmetrically with respect to the free surface. On the other hand, a nonzero Φ0 is necessary for 1:2 internal resonance [70]. Next, we substitute (3.15) into (3.12) to obtain the residue term R(𝜏) = Mx′′ 𝑝 + 𝑝 + Cx′ (cid:16) + Kx𝑝 − g x𝑝, x′ ∫ ∞ 𝑠=0 𝑝, x′′ 𝑝 K𝑅 (𝑠) x′ (𝜏 − 𝑠) 𝑑𝑠 (cid:17) − f. (3.16) 47 where x𝑝 = (cid:2)𝜃 𝑝, 𝜙𝑝(cid:3)𝑇 is the assumed solution vector. Using Galerkin’s projection on the orthogonal trigonometric basis, we can obtain ℎ0 (X) = ℎ𝑠 𝑝 (X) = ℎ𝑐 𝑝 (X) = ∫ 2 𝜋𝜈 𝜔 0 ∫ 2 𝜋𝜈 𝜔 0 ∫ 2 𝜋𝜈 𝜔 0 R(𝜏)𝑑𝜏 = 0, R(𝜏) sin R(𝜏) cos (cid:17) (cid:16) 𝑝𝜔𝜏 𝜈 (cid:17) (cid:16) 𝑝𝜔𝜏 𝜈 𝑑𝜏 = 0, 𝑑𝜏 = 0 where X = [X0, X𝑐 1 , X𝑠 1 , . . . , X𝑐 𝑝, X𝑠 𝑝]𝑇 , where X0 = [Θ0, Φ0], X𝑐 1 = (cid:2)Θ𝑐 1 (3.17) , Φ𝑠 1 (cid:3), and , Φ𝑐 1 (cid:3), X𝑠 1 = (cid:2)Θ𝑠 1 so on. Eqn. (3.17) is a set of nonlinear algebraic equations and can be solved by iterative methods like the Newton-Raphson scheme. However, since the term g (cid:16) x𝑝, x′ 𝑝, x′′ 𝑝 (cid:17) will result in composite trigonometric terms such as cos (cid:16) Φ𝑠 𝑝 sin ( 𝑝𝜔𝜏/𝜈) (cid:17) , one needs to use special expansions (example, using Bessel functions) to extract the Fourier coefficients of g (denoted by G). Moreover, for stability analysis, one would require to compute the Jacobian matrix 𝜕G 𝜕X , which is computationally expensive in the frequency domain. To that end, we use a modification of the alternating frequency/time domain (AFT) method [78, 79, 64] to compute the value of G and 𝜕G 𝜕X . Furthermore, to apply the harmonic balance method, we would require to obtain the Fourier coefficients of the convolution term. This calculation is explained in Section 3.3.3. 3.3.2 Modified alternating frequency/time approach We follow the AFT algorithm as implemented in [64] with a modification that takes care of nonlinear inertia. Let us define ¯𝜔 = 𝜔 𝜈 and expand the nonlinear function by Fourier series such that (cid:16) g x𝑝, x′ 𝑝, x′′ 𝑝 (cid:17) = G0 + 𝑃 ∑︁ (cid:16) 𝑝=1 which can be rewritten as G𝑐 𝑝 cos ( 𝑝 ¯𝜔𝜏) + G𝑠 𝑝 sin ( 𝑝 ¯𝜔𝜏) (cid:17) (cid:16) g x𝑝, x′ 𝑝, x′′ 𝑝 (cid:17) = (T ( ¯𝜔𝜏) ⊗ 𝐼𝑛) G 48 where T( ¯𝜔𝜏) = [1 cos ( ¯𝜔𝜏) sin ( ¯𝜔𝜏) . . . cos (𝑃 ¯𝜔𝜏) sin (𝑃 ¯𝜔𝜏)] G = [G𝑇 0 , (cid:0)G𝑐 1 (cid:1)𝑇 , (cid:0)G𝑠 1 (cid:1)𝑇 , . . . , (cid:0)G𝑐 𝑃 (cid:1)𝑇 , (cid:0)G𝑠 𝑃 (cid:1)𝑇 ]𝑇 (3.18) where ⊗ is the Kronecker tensor product and 𝐼𝑛 is the 𝑛 × 𝑛 identity matrix. Similarly, x = (T ( ¯𝜔𝜏) ⊗ 𝐼𝑛) X where X is the vector consisting of all retaining Fourier coefficients; see (3.17). The AFT method finds the time domain value of the nonlinear function and transforms it back to the frequency domain, i.e. 𝐷𝐹𝑇 −1 −−−−−→ x𝑝 (𝑡), x𝑝 (𝑡)′, x𝑝 (𝑡)′′ −→ g X (cid:16) x𝑝, x′ 𝑝, x′′ 𝑝 (cid:17) 𝐷𝐹𝑇 −−−−→ G(X) (3.19) where 𝐷𝐹𝑇 represents the discrete Fourier transform. Taking 𝑁 uniformly separated time points, 𝑡𝑖 = 2𝜋𝑖 𝜔𝑁 , and defining the vectors 𝑝 (𝑡1), x𝑇 𝑝 (𝑡2) . . . , x𝑇 𝑝 (𝑡𝑁 )(cid:3)𝑇 ¯x = (cid:2)x𝑇 ¯g = (cid:2)g𝑇 (𝑡1), g𝑇 (𝑡2) . . . , g𝑇 (𝑡𝑁 )(cid:3)𝑇 . and (cid:17) (cid:17) sin cos 1 ...           the Fourier coefficients can be computed using (cid:16) 2𝜋 𝑁 ... (cid:16) 2𝜋𝑁 𝑁 (cid:16) 2𝜋 𝑁 ... (cid:16) 2𝜋𝑁 𝑁 1 cos U = sin (cid:17) (cid:17) . . . (cid:17) (cid:16) 2𝜋𝑃 𝑁 cos (cid:17) (cid:16) 2𝜋𝑃 𝑁 sin . . . cos (cid:17) (cid:16) 2𝜋𝑃𝑁 𝑁 (cid:17) (cid:16) 2𝜋𝑃𝑁 𝑁 sin ,           (cid:16) U−1 ⊗ I𝑛 G = (cid:16) U−1 ⊗ I𝑛 X = (cid:17) (cid:17) ¯g ¯x 49 (3.20) (3.21) (3.22) (3.23) where U−1 = 1/𝑁 1 2 cos( 2𝜋 𝑁 ) 2 sin( 2𝜋 𝑁 ) ... . . . . . . . . . ... 1 2 cos( 2𝜋𝑁 𝑁 ) 2 sin( 2𝜋𝑁 𝑁 ) ... 2 cos( 2𝜋𝑃 𝑁 ) 2 sin( 2𝜋𝑃 𝑁 ) . . . 2 cos( 2𝑃𝜋𝑁 𝑁 ) . . . 2 sin( 2𝑃𝜋𝑁 𝑁 )                                         is the Moore-Penrose pseudo-inverse of U. Also, where ¯x′ = 𝜔 (U∇ ⊗ I𝑛) 𝑋, ¯x′′ = 𝜔2 (cid:16) U∇2 ⊗ I𝑛 (cid:17) 𝑋, . . . ∇𝑖 ∇ =  0                                . . . ∇𝑃 , ∇𝑖 = Thus, the Jacobian matrix of G can be obtained through 𝑖 0     −𝑖 0    .        𝜕G 𝜕 𝑋 (cid:16) = U−1 ⊗ I𝑛 (cid:16) + U−1 ⊗ I𝑛 where (cid:17) 𝜕 ¯g 𝜕 ¯x (cid:17) 𝜕 ¯g 𝜕 ¯x′′ (U ⊗ I𝑛) + 𝜔2 (cid:104) (cid:16) U∇2(cid:17) (cid:16) U−1 ⊗ I𝑛 (cid:17) 𝜕 ¯g 𝜕 ¯x′ (cid:105) , ⊗ I𝑛 𝜔 [(U∇) ⊗ I𝑛] 𝜕 ¯g 𝜕r = 𝜕 ¯g 𝜕r (cid:12) (cid:12) (cid:12) (cid:12)𝑡=𝑡1 . . . 𝜕 ¯g 𝜕r (cid:12) (cid:12) (cid:12) (cid:12)𝑡=𝑡 𝑁                         are the Jacobian of the nonlinear function g in the time domain and r can be ¯x, ¯x′ and ¯x′′. 50 (3.24) (3.25) (3.26) (3.27) (3.28) 3.3.3 Fourier coefficients of the convolution term We use the method described in [74] to obtain the Fourier coefficients of the convolution term. After substituting (3.15) into the convolution term, one arrives at 𝐼 = ∫ ∞ 0 𝜅𝑅 (𝑠) (cid:104) 𝑝 ¯𝜔 𝑃 ∑︁ 𝑝=1 − Θ𝑐 𝑝 sin ( 𝑝 ¯𝜔 (𝜏 − 𝑠)) + Θ𝑠 𝑝 cos ( 𝑝 ¯𝜔 (𝜏 − 𝑠)) (cid:105) 𝑑𝑠 (3.29) where 𝐼 = ∫ ∞ 0 𝜅𝑅 (𝑠) 𝜃′ 𝑝 (𝜏 − 𝑠) 𝑑𝑠. Expanding the trigonometric functions and using (3.10) lead to 𝑃 ∑︁ (cid:104) ¯𝜔 𝑝 𝐼 = −Θ𝑐2𝜉 ( 𝑝) 𝐵 (cid:0)𝜔0 ¯𝜔 𝑝(cid:1) + Θ𝑠 ¯𝜔 𝑝 𝜇( 𝑝) 𝐴 (cid:0)𝜔0 ¯𝜔 𝑝(cid:1)(cid:105) sin (cid:0) ¯𝜔 𝑝𝜏(cid:1) 𝑝=1 + ¯𝜔 𝑝 (cid:104) Θ𝑐 ¯𝜔 𝑝 𝜇( 𝑝) 𝐴 (cid:0)𝜔0 ¯𝜔 𝑝(cid:1) + Θ𝑠2𝜉 ( 𝑝) 𝐵 (cid:0)𝜔0 ¯𝜔 𝑝(cid:1) (cid:105) cos (cid:0) ¯𝜔 𝑝𝜏(cid:1) (3.30) where ¯𝜔 𝑝 = 𝑝 ¯𝜔, 𝜇( 𝑝) 𝐴 (cid:0) ¯𝜔 𝑝(cid:1) = 𝐴∞ − 𝐴𝑅 (cid:0)𝜔0 ¯𝜔 𝑝(cid:1) 𝑀 𝐵𝑅 (cid:0)𝜔0 ¯𝜔 𝑝(cid:1) 𝑀𝜔0 . 2𝜉 ( 𝑝) 𝐵 (cid:0) ¯𝜔 𝑝(cid:1) = , (3.31) Substituting (3.30) into (3.17) leads to the Fourier coefficients of the convolution term. 3.3.4 Stability of the periodic solutions When the nonlinearity in the system is weak, the periodic solutions (3.15) are dominated by the primary harmonics. Therefore, 𝜈 = 1, 𝑝 = 1 are chosen in this study. We determine the stability of the periodic solutions by introducing small perturbations in (3.15) as follows: 𝜃 (𝜏) = 𝜃 𝑝 (𝜏) + 𝛿𝜃 (𝜏) and 𝜙(𝜏) = 𝜙 𝑝 (𝜏) + 𝛿𝜙 (𝜏) (3.32) where |𝛿𝜃 (𝜏)| << 1 and |𝛿𝜙 (𝜏)| << 1. Let us define 𝜹 = (cid:2)𝛿𝜃, 𝛿𝜙(cid:3)𝑇 . Substitution of (3.32) into (3.12) and linearization with respect to 𝜃 𝑝 (𝜏) and 𝜙 𝑝 (𝜏) yield (cid:18) M + A∞ + (cid:19) 𝜕g 𝜕x′′ (cid:18) C + 𝜹′′ + (cid:19) 𝜕g 𝜕x′ 𝜹′ + ∫ ∞ 𝜏=0 K𝑅 (𝑠) 𝜹′ (𝜏 − 𝑠) 𝑑𝑠 + (cid:18) K + (cid:19) 𝜕g 𝜕x 𝜹 = 0 51 (3.33) where the Jacobian matrices 𝜕g/𝜕x′′, 𝜕g/𝜕x′, and 𝜕g/𝜕x are evaluated at x = x𝑝, x′ = x′ 𝑝, and x′′ = x′′ 𝑝 using the AFT method described in Section 3.3.2. Eqn. (3.33) is a set of linear ordinary differential equations of periodic coefficients. One can use the Floquet theory to determine the stability. To this end, (3.33) is transformed into the state-space form and numerically integrated over one period 𝑇 to obtain the fundamental matrix. It is, however, difficult to numerically integrate the convolution term. To overcome this difficulty, we use the method illustrated in [80] to obtain a state-space representation that governs the dynamics of the convolution term, which is briefly explained as follows. Denote by 𝑦(𝜏) the convolution term in (3.33), i.e. 𝑦(𝜏) = ∫ ∞ 0 𝜅𝑅 (𝑠) 𝛿′ 𝜃 (𝜏 − 𝑠) 𝑑𝑠 = ∫ 𝑡 −∞ 𝜅𝑅 (𝜏 − 𝑠) 𝛿′ 𝜃 (𝑠) 𝑑𝑠 which, after Laplace transform, gives, ˜𝑦 ( ˜𝑠) = ˜𝐻 ( ˜𝑠) ˜𝛿′ 𝜃 ( ˜𝑠) where ˜𝑦 and ˜𝜹 ′ are the Laplace transform of 𝑦(𝜏) and 𝛿′(𝜏) respectively, and ˜𝐻 ( ˜𝑠) = ∫ ∞ 0 𝜅𝑅 (𝜏) 𝑒 ˜𝑠𝜏𝑑𝜏. Using this transfer function, we find a state-space realization of order 𝑛 with state variable 𝛿𝜃 (𝜏) and matrices A ∈ R𝑛×𝑛, B ∈ R𝑛×1, C ∈ R1×𝑛 and 𝐷 ∈ R, such that w′ (𝜏) = Aw (𝜏) + B𝛿′ 𝜃 (𝜏) 𝑦 (𝜏) = Cw (𝜏) + 𝐷𝛿′ 𝜃 (𝜏) (3.34) The state space conversion mentioned above can be achieved using imp2ss function in Matlab. To use imp2ss, 𝜅𝑅 (𝜏) can be calculated using (3.10). We can further reduce the above state space model using balmr function in Matlab. After getting a reduced state space model for the convolution integral, the convolution term can be characterized by the state space model (3.34). As stated previously, when the nonlinearity is weak, the solutions are dominated by the primary harmonics. Therefore, we determine the stability of the primary harmonics. Specifically, we will determine when the primary harmonics undergo period-doubling bifurcation. When this occurs, the system will have 1:2 internal resonance [51, 57, 58, 70]. 52 Figure 3.5 Parametric instability boundary for 𝜂 = 0.4, 𝜇𝑟 = 0.4, 𝜇𝑔 = 0.03, 𝜉 = 0.05, 𝜉𝑝 = 0.01 3.3.5 Period doubling bifurcation We use the bifurcation tracking algorithm developed in [70] to obtain the bifurcation boundary for period-doubling bifurcation in the 𝑓𝑒 −𝜔 plane. The bifurcation boundary for a set of parameters is shown in Fig. 3.5. When this bifurcation occurs, the pendulum oscillations will have harmonics of 𝜔/2, i.e. 1:2 internal resonance. It is worth noting that autoparametric vibration absorbers also use 1:2 internal resonance to achieve vibration energy absorption [51, 57, 58]. 3.4 Numerical demonstration We now perform direct numerical simulation on the system (using Matlab’s ODE45) to verify the bifurcation boundary obtained in Section 3.3.5. We solve (3.12) at three points marked in Fig. 3.5 (×1, ×2 and ×3). To solve the system numerically, we obtain the state space representation of the convolution term in (3.12) in a way similar to one described in Section 3.3.4, thus giving ˜w′ (𝜏) = ˜A ˜w (𝜏) + ˜B𝜃′ (𝜏) 𝑧 (𝜏) = ˜C ˜w′ (𝜏) + ˜𝐷𝜃′ (𝜏) (3.35) where 𝑧 (𝜏) ≈ ∫ ∞ 0 𝜅𝑅 (𝑠)𝜃′ (𝜏 − 𝑠) 𝑑𝑠. As can be observed from Fig. 3.6a, Fig. 3.6b, Fig. 3.7a, Fig. 3.7b, Fig. 3.8a, and Fig. 3.8b, points ×1 and ×2 lead to periodic solutions whereas point ×3 leads to non-periodic solutions. This can be validated with the fast Fourier transform (FFT) of the solutions at these three points, which 53 0.50.60.70.80.911.10.050.10.150.20.250.3fe123Applied force is shown in Fig. 3.6c, Fig. 3.7c and Fig. 3.8c, respectively. Note that the frequencies ˆ𝜔 of the FFT are normalized with respect to the excitation frequency. Thus, primary harmonics correspond to components at ˆ𝜔 = 1, harmonics of half excitation frequency correspond to components at ˆ𝜔 = 0.5, and so on. In Fig. 3.5, it can be seen that the point ×1 is below the bifurcation boundary. We expect the solutions to be periodic and dominated by the primary harmonic, which is verified by Fig. 3.6. Now, as we increase the value of 𝑓𝑒 to reach point ×2, a period-doubling bifurcation occurs and we expect harmonics of frequency 𝜔/2. This claim is readily verified from Fig. 3.7, where subharmonics of half excitation frequency indeed exist. As such, it is found that the pendulum has 1:2 internal resonance at point ×2. Finally, the parameters at ×3 lead to strong non-periodic solutions composed of both oscillation and intermittent rotations of the pendulum, as shown in Fig. 3.8a. This is similar to the non-periodic solutions observed in autoparametric resonance vibration absorbers [58]. a Time series of 𝜙 at point 1 b Time series of 𝜃 at point 1 c FFT at point 1 Figure 3.6 FFT and time series of periodic solutions at point 1 in Fig. 3.5 3.5 Discussion In this study, we investigate how the 1:2 internal resonance of the IPVA-PTO can be exploited for wave energy conversion. To show the efficacy of the system, we compare it with a linear wave energy converter. The linear system is characterized by a linear wave energy converter that is obtained by removing the pendulum from the IPVA, as shown in Fig. 3.1c. We tune the linear inerter to the same natural frequency as the IPVA for a more fairer comparison. This is achieved by increasing the inertance of the generator. Considering the other physical parameters same as the 54 2.1052.112.1152.12t1040.70.80.911.11.21.31.42.1152.122.1252.13t104-0.3-0.2-0.100.10.20.300.511.522.5310-610-410-2100FFTFFT of FFT of a Time series of 𝜙 at point 2 b Time series of 𝜃 at point 2 c FFT at point 2 Figure 3.7 FFT and time series of periodic solutions at point 2 in Fig. 3.5 a Time series of 𝜙 at point 3 b Time series of 𝜃 at point 3 c FFT at point 3 Figure 3.8 FFT and time series of periodic solutions at point 3 in Fig. 3.5 IPVA-PTO, the equation of motion of the linear system will be 𝑀𝑙𝜃′′ + 𝐶𝑙𝜃′ + ∫ ∞ 𝜏=0 = 𝑓𝑒 𝜅𝑅 (𝑠) 𝜃′ (𝜏 − 𝑠) 𝑑𝑠 + 𝐾𝑙𝜃 where 𝑀𝑙 = 1 + 𝜇𝑔 + 𝜇 𝐴∞ + M, 𝜇 𝐴∞ = 𝐴∞ 𝑀 , 𝐶𝑙 = 2 (𝜉 + 𝜉𝑒) , 𝐾𝑙 = 1. (3.36) (3.37) and M is the additional inertia used for resonant frequency matching. We solve this system with the method described in Section 3.3.4 and compare it with the IPVA-PTO system. Note that the excitation force is also frequency dependent. For the current study, we fix the value of 𝑅 to 28 cm. Fig. 3.5 shows the force as a variation of the normalized angular frequency of the system. When the excitation force magnitude is above the bifurcation boundary, we expect that the vibration energy 55 2.0852.092.0952.12.1052.11t1040.511.52.1242.1262.1282.132.1322.134t104-0.3-0.2-0.100.10.20.300.511.522.5310-610-410-2100FFTFFT of FFT of 11.051.11.151.21.251.3t104-20-10010201.9051.911.9151.921.9251.93t104-0.4-0.200.20.400.511.522.5310-610-410-2100102FFTFFT of FFT of of the spar transfers to the IPVA. To see that, we compare the response amplitude operator (RAO) in heave and the normalized power of the IPVA-PTO and the linear system. The RAO is defined as and the capture width is defined as 𝑅 𝐴𝑂 = 𝑅𝜃 𝛾 , 𝑃 = 𝑃𝐼 𝑃𝑉 𝐴−𝑃𝑇𝑂/𝑃𝑙𝑖𝑛𝑒𝑎𝑟 (3.38) (3.39) which is the ratio of the power that the IPVA-PTO converts over the power that the linear benchmark converts. The instantaneous converted wave power for the IPVA and linear system are prospectively written as 𝑃𝐼 𝑃𝑉 𝐴−𝑃𝑇𝑂 = 𝑐𝑒 ( (cid:164)𝜃 − (cid:164)𝜙)2 𝑃𝑙𝑖𝑛𝑒𝑎𝑟 = 𝑐𝑒 ( (cid:164)𝜃)2 (3.40) In this work, the root-mean-square (RMS) of the converted power for 600 wave periods is computed for both linear and nonlinear systems. As can be seen from the flattening of the oscillations in Fig. 3.9a (saturation like phenomenon as previously reported in [70]) of the spar and the increase in the capture width within the same frequency range in Fig. 3.9b, there is an energy transfer from the spar to the pendulum. Note that the capture width is normalized by the maximum of the linear system for easy reference. Furthermore, it can be observed from Fig. 3.9b that the IPVA-PTO system significantly outperforms the linear system. Note that in Fig. 3.9a, approximately between Ω ∈ [2.3, 2.8], there is a presence of chaotic-like motions. The figures represent one such instance of the initial conditions. 3.6 Experimental analysis To verify the analysis performed, experiments consisting of the IPVA-PTO integrated with a “dry” single-degree-of-freedom (sdof) system (without any hydrodynamic effects on the system) are conducted. The “dry” system should still show the internal resonance phenomenon and it can be explained as follows. The hydrodynamic effects on the system add a frequency-dependent mass, 56 a Comparison of the response amplitude operator of the IPVA- PTO system and the linear system to different angular with respect frequencies b Normalized capture width of the IPVA-PTO system and the linear system with respect to different angular frequencies Figure 3.9 Comparison of response amplitude operator (RAO) and capture width between the IPVA-PTO and linear system damping and a force to the system. If we fix a frequency value - then the mass, the stiffness, and the damping matrix can be calculated at that frequency value, and a force can be obtained corresponding to the boundary point of the internal resonance boundary. This boundary can be obtained for all the frequency values as explained. Now, if the mass, stiffness, and damping matrix do not change with frequency, we can still obtain a force corresponding to internal resonance boundary at a specified frequency, and do the procedure for all the frequencies to obtain an internal resonance boundary. Thus, in principle, the internal resonance boundary for the dry system can be obtained. Although it may be different than the system with hydrodynamic effects, the boundary will still exist. The aim of the experiments is to characterize and compare the energy conversion and response suppression capabilities of the IPVA-PTO system with a linear benchmark. The experimental setup is shown in Fig. 3.10. The caption of Fig. 3.10 shows the description of labels for various components of the experimental setup. The top plate, marked by plate 𝐴, supports the primary mass. The base plate system contains three plates, plates 𝐶, 𝐸, and 𝐹. Note that the plates 𝐶, 𝐸 and 𝐹 are rigidly connected while Plate 𝐹 is directly connected to the shaker. The nut of the ball-screw system connects the plate 𝐴 while the screw is supported by a thrust bearing attached to 57 22.533.5012345Response Amplitude Operator (RAO)IPVA-PTOLinear22.533.500.511.5Normalized Capture WidthIPVA-PTOLinear Figure 3.10 Experimental setup. The labels denote the following — 1: Shaker, 2: Spectral analyzer and controller, 3: Accelerometer signal conditioner, 4: Shaker signal amplifier, 𝐴: Top plate, 𝐵: Ball-screw system, 𝐶: Ball-screw mounting plate , 𝐷: Sun and planet gear, 𝐸: Middle plate, 𝐹: Base plate, 𝐺: Primary mass, 𝐻: Connecting springs, 𝐼: Pendulum, 𝐽: Carrier, 𝐾: Generator, 𝐿: Load resistance plate 𝐶. A coupler connects the ball-screw system with the carrier (𝐽), which hosts the pendulum. Eight springs are connected between the primary mass system and the base plate system so that the two systems can move relative to each other. The relative motion drives the screw resulting in 𝜃 = 𝑥 − 𝑦 𝑅 , where 𝑥 is the motion of plate 𝐴 and 𝑦 is the motion of plate 𝐶. The pendulum (𝐼) is free to rotate, and its shaft is supported by the carrier through a ball bearing. A planetary gear system is hosted on the pendulum and the rotor of the generator (𝐷), whose housing is fixed to plate 𝐸. The sun gear and planet gear are mounted to the generator shaft and the pendulum shaft, respectively. A load resistance (𝑅𝐿) is connected to the generator terminals to measure the harvested electric power. Finally, a shaker (APS 113) excites the base plate system by controlling the motion of the plate 𝐹 and is driven by a controller (Spider 80x) using an amplifier. There is an accelerometer connected to the base plate for closed-loop control of the shaker, and another accelerometer is connected to the top plate to monitor the motion of the top plate system. The generator terminals are connected to a spectral analyzer (Spider 80x) to measure its voltage. 58 is given by where ¯M = ¯K =               3.6.1 Equation of motion for the experimental setup The equation of motion of the experimental system can be directly derived from (3.12) by converting the forced excitation into base excitation (of the shaker), adding friction to the pendulum shaft and the generator shaft, and removing the hydrodynamic effects. To obtain this equation of motion, we assume perfect transmission by the ball-screw system. The resultant equation of motion ¯M (cid:165)X + ¯C (cid:164)X + ¯KX + ¯G(X, (cid:164)X, (cid:165)X) = ¯F (3.41) 1 + 𝜇𝑟 + 𝜇𝑟𝜂2 + 2𝜂𝜇𝑟 cos(𝜙) + 𝜇𝑏𝑠𝑐 + 𝜇 𝑝 + 𝜇𝑔 𝜇 𝑝 + 𝜇𝑟𝜂2 + 𝜇𝑟𝜂 cos(𝜙) − 𝜇𝑔 −𝜉𝑒 2 (cid:0)𝜉𝑒 + 𝜉 𝑝(cid:1)        , ¯G = 𝜇𝑟𝜂 (cid:169) (cid:173) (cid:173) (cid:171)        1 0 0 0 , ¯C = 2 (𝜉 + 𝜉𝑒)        −𝑦′′ − 𝑡𝑔sign (cid:0) (cid:164)𝜃 − (cid:164)𝜙(cid:1) −𝑡 𝑝sign (cid:0) (cid:164)𝜙(cid:1) + 𝑡𝑔sign (cid:0) (cid:164)𝜃 − (cid:164)𝜙(cid:1) −𝜉𝑒 (cid:170) (cid:174) (cid:174) (cid:172) ¯F = (cid:169) (cid:173) (cid:173) (cid:171) with X = [𝜃, 𝜙]𝑇 , 𝑡 𝑝 = 𝜇 𝑝 + 𝜇𝑟𝜂2 + 𝜇𝑟𝜂 cos(𝜙) − 𝜇𝑔        −2𝜃′𝜙′ sin(𝜙) − 𝜙′2 sin(𝜙) 𝜇 𝑝 + 𝜇𝑟𝜂2 + 𝜇𝑔 𝜃′2 sin(𝜙) , , (cid:170) (cid:174) (cid:174) (cid:172) (3.42) 𝑇𝑔 𝑀 𝑅2 , where 𝑇𝑝 is the frictional torque between the pendulum shaft and its bearing and 𝑇𝑔 is that of the generator rotor. Here, 𝑥 is the motion of the top plate, 𝑦 𝑇𝑝 𝑀 𝑅2 and 𝑡𝑔 = is the motion of the shaker, 𝜃 is the rotation of the carrier, 𝑅𝑝 is the location of the pendulum shaft shaft defined by 𝜉 𝑝 = with respect to the center of the ball-screw system and 𝜉 𝑝 is the damping ratio of the pendulum 𝑐 𝑝 𝑀𝜔0𝑅2 , 𝑐 𝑝 being the damping coefficient between the pendulum shaft and 𝑐𝑒 𝑀𝜔0𝑅2 , where 𝑐𝑒 is the electrical damping in the generator due to its internal resistance and well as the load resistance and can be written as 𝑐𝑒 = 𝜅2/(𝑅𝑖𝑛𝑡 + 𝑅𝐿), ball bearing. Furthermore, 𝜉𝑒 = where 𝜅 is the torque constant of the generator, and 𝑅𝑖𝑛𝑡 and 𝑅𝐿 are internal and load resistance respectively. 𝑀 is the mass of the top plate and the primary mass (i.e., the primary structure), 𝐽𝑏𝑠𝑐 includes the inertia of the carrier, ball-screw, and coupler, and 𝐽𝑝 is the moment of inertia of the pendulum system with respect to its center of mass. We assume theoretical values for all the inertia (except 𝐽𝑏𝑠𝑐), masses, transmission ratio, springs, geometric lengths, and electrical damping 59 in the system. Therefore, the values of 𝑚, 𝑅𝑝, 𝑀, 𝑅, 𝐽𝑝, 𝑘, 𝑟, 𝑐𝑒 and 𝜔0 are assumed to be known; see Table 3.1. This leaves 𝐽𝑏𝑠𝑐, 𝑐 and 𝑐 𝑝 to be experimentally identified. The calculation of pendulum damping 𝑐 𝑝 has already been discussed in [81], and the same procedure was used to obtain the pendulum damping in this study. Therefore, in the next sub-section, we will discuss the identification of 𝑐 and 𝐽𝑏𝑠𝑐. 3.6.2 Obtaining the mechanical damping 𝑐 and ball-screw-carrier inertia 𝐽𝑏𝑠𝑐 The mechanical damping value 𝑐 is obtained by removing the pendulum head and the horizontal shaft (attached to the pendulum). This leaves three inertias on the linear system — the primary mass (which includes the top plate 𝐴), the ball-screw-carrier assembly, and the vertical shaft (where the pendulum’s horizontal shaft was attached originally) with the ball-bearing supporting it. From this experimental setup, the frequency response function (FRF) was obtained. The measured FRF and theoretical FRF are correlated with each other to obtain the mechanical damping value as follows. The linear system without the pendulum has the following equation of motion where ˆ𝑀 (cid:165)𝑧 + 𝐶 (cid:164)𝑧 + 𝐾𝑧 = −𝑀𝛼 (cid:165)𝑦 ˆ𝑀 = 𝑀 (1 + 𝛼) , 𝛼 = 𝜇𝑏𝑠𝑐 + 𝐽𝑝𝑠 𝑀 𝑅2 (3.43) (3.44) with 𝐽𝑝𝑠 being the inertia value of the vertical pendulum shaft. The FRF is obtained by substituting 𝑦 = 𝑌 𝜔2 𝑒𝑖𝜔𝑡 and 𝑧 = 𝑍 𝜔2 𝑒𝑖𝜔𝑡 into (3.43), where 𝑌 and 𝑍 denote the displacement amplitude of the base plate and the amplitude of the relative displacement between the top plate and base plate, respectively. By using the peak amplitude method [82], one can find 𝛼 (and therefore 𝐽𝑏𝑠𝑐, assuming theoretical value for 𝐽𝑝𝑠) and 𝑐. Ten different FRF experiments were performed and 𝑐 and 𝐽𝑏𝑠𝑐 values were obtained by averaging their individual experimental values. The experimental parameters are tabulated in Table 3.1 for ready reference. 3.6.3 Verification of the internal resonance In this section, the 1:2 internal resonance phenomenon is verified against the experiments. First, the PD boundary is obtained and shown in Fig. 3.11 using the theoretical and fitted parameters 60 Parameter 𝜉 𝑝 𝜉 𝜇𝑟 𝜂 𝜔0 𝜇𝑔 𝑇𝑝 𝑇𝑔 𝜇 𝑝 𝜇𝑏𝑠𝑐 Value 0.00418 0.1266 6.571 0.55 85.74 rad/s 0.442 0.0014N-m 4.26 × 10−5N-m 0.756 4.5117 Table 3.1 Parameters for the experimental IPVA-PTO system shown in Table 3.1 for two values of load resistance — 5 ohm and 10 ohm. Since it is known that there should be response suppression above the PD boundary, various acceleration values are used to run the experiments, including an acceleration of 0.6 𝑔 for Fig. 3.11a and an acceleration of 0.9 𝑔 for Fig. 3.11b, marked by horizontal line (—) on both the figures. Before observing the experimental results and comparing them to the equivalent linear benchmark, a few points need to be stated: 1. The ball-screw system is assumed to be 100% efficient, with no loss in transmission 2. The effects of accelerometer cables on the system are neglected 3. The springs are modeled to be linear for the analytical analysis As reported in [70], for a sdof system, 1:2 internal resonance is accompanied by an energy transfer between the primary system and the pendulum vibration absorber and thus suppressing the response of the primary system. To observe the secondary harmonic solutions (defined by the response of the system immediately after the period-doubling bifurcation, containing harmonics of half the excitation frequency along with the excitation frequency) of the system, the PD boundary shown in Fig. 3.11a is referred to. As shown, a marker × slightly above the boundary corresponds to an excitation acceleration amplitude of 0.42 𝑔 and an excitation frequency of 3.6 Hz. This set of acceleration and frequency are chosen as they are just above the boundary, where we expect 61 a 10 ohm (𝜉𝑒 = 0.0737) b 5 ohm (𝜉𝑒 = 0.0889) Figure 3.11 PD boundary of energy transfer corresponding to the experimental setup for two values of load resistance a Time series of the pendulum motion b Measured FFT of the pendulum motion Figure 3.12 Measured time series and FFT of the output voltage at 0.42 𝑔, at 3.6 Hz to observe secondary harmonic solutions. As can be seen from Fig. 3.12, the voltage measured by the spectral analyzer shows harmonics of frequency 𝜔 2 as evident from the time series of the pendulum’s motion shown in Fig. 3.12a and the fast Fourier transform (FFT) shown in Fig. 3.12b. This verifies the phenomenon of internal resonance. 3.6.4 Comparison with the linear system Next, the IPVA-PTO system is compared with a linear system to benchmark its response suppression and energy harvesting capabilities. For this analysis, the linear system is chosen by removing the pendulum and adding inertia to the ball-screw-carrier system such that the natural 62 33.23.43.63.844.2Frequency (Hz)0.40.450.50.550.60.65Acceleration (g)33.23.43.63.844.2Frequency (Hz)0.90.9511.051.11.151.21.251.3Acceleration (g)35.83636.236.436.636.83737.237.437.6t (s)-2-1.5-1-0.500.511.52Voltage (V)12345678f (Hz)10-210-1100FFT of VoltageX 3.60806Y 1.21278X 1.79487Y 0.123049X 5.40293Y 0.0597345X 7.1978Y 0.240968 frequency of the linear system is equal to the resonance frequency of the IPVA-PTO system. The internal resistance of the generator is 3.17 ohms. A few resistance values (3, 5, and 10 ohm, with optimal resistance value for maximum energy of the linear system close to 5 ohm) were chosen to perform experiments for the linear system for a better comparison with the IPVA-PTO system. Note that the response of the IPVA-PTO system can be categorized into three types. First, the pendulum can have primary harmonic oscillations where the oscillation frequency of the pendulum is equal to the excitation frequency. Second, the pendulum can have a secondary harmonic response due to 1:2 internal resonance where the pendulum’s response contains a frequency of half the excitation frequency in tandem with the primary harmonic response. Lastly, the pendulum can go into non- periodic motion, consisting of a rich set of frequencies, or can go into rotation. The methodology applied to generate the results for the IPVA-PTO system can be explained as follows: For a given excitation and load resistance value, frequency response functions (FRFs) of the voltage and the top plate acceleration are generated with respect to the shaker’s acceleration using sine sweep input. Then, the frequency range within which the response consisting of only primary harmonics is identified. For the identified frequency range, the root-mean-square (RMS) voltage and relative displacement between the top plate and base plate response are calculated from the mean of the RMS responses obtained for ten independent sine sweeps. For the frequency range within which the response is secondary/non-periodic/rotating, the range is discretized into frequency points with a step-size of 0.1 Hz, and experiments are run on the IPVA-PTO system to get the RMS values of the amplitude of the relative displacement between the top plate and the shaker, and the generated power using a data of approximately 120 seconds per experiment (total data of 1200 second for a frequency, with transients removed). Ten frequency response functions were calculated for the linear system to obtain the RMS power and relative displacement between the top plate and the shaker. All the RMS values are plotted against each other for comparison. First, let us look at the case when the load resistance is 10 ohms. Also, let us choose a large acceleration to observe non-periodic solutions in the pendulum’s response. As can be seen from Fig. 3.11a, the acceleration of 0.6 𝑔 is well above the PD boundary from 3.2 Hz to 4 Hz, 63 a Comparison between measured RMS of relative motion between the top and the base plate for linear and IPVA-PTO system b Comparison between measured RMS power value for the linear and the IPVA-PTO system Figure 3.13 Frequency response of relative motion between the top plate and the base plate, and power for the IPVA-PTO and the linear system at 0.6 𝑔 base excitation acceleration with 10 ohm resistance. 1:2 stands for secondary harmonic solutions, NP stands for non-periodic solution Comparison between a experimental relative motion between the top plate and the base plate RMS value for linear and the IPVA system Voltage Comparison between b experimental RMS value for the linear and the IPVA system. The legend is same as part a of the figure Figure 3.14 Frequency response of relative motion between the top plate and the base plate, and power for the IPVA-PTO and the linear system at 0.9 𝑔 base excitation acceleration with 5 ohm resistance and 3 ohm resistance for the linear system and 5 ohm for the IPVA system and the response goes into non-periodic motion, with secondary harmonic at some instances (as observed from the time series of the response, not shown here). The system shows excellent response suppression compared to the linear system as evident from Fig. 3.13a. Although the 64 33.23.43.63.844.2Frequency (Hz)246810121416Relative motion b/w top and base plate (m)10-3IPVA-PTOLinearNP1:233.23.43.63.84Frequency (Hz)00.050.10.150.20.250.30.350.40.45Power (W)IPVA-PTOLinearNP1:233.23.43.63.844.2Frequency (Hz)12345678Relative motion b/w top and base plate (m)10-3IPVA-PTO 5 ohmIPVA-PTO 5 ohmLinear 5 ohmLinear 3 ohm33.23.43.63.844.2Frequency (Hz)00.10.20.30.40.50.60.70.80.91Power (W) response suppression is significant, the power generated is in question. There can be instances where the power generated outperforms the linear system — when the pendulum goes into 1:2 internal resonance (secondary harmonic response), but there are cases when the power generated is less than the linear system, for non-periodic motion. For the next case, the load resistance for the IPVA-PTO system is 5 ohm, and thus the electrical damping in the IPVA-PTO system is higher as a consequence. Two resistance values are used for the linear system, 3 and 5 ohm, for the sake of comparison. Due to larger electrical damping, a higher input acceleration is required to cross the PD boundary. This time, we choose 0.9 𝑔 as the acceleration magnitude as it is close to the boundary, and therefore we can observe secondary harmonic solutions. In Fig. 3.14, it is observed that the system shows primary harmonic response up to around 3.4 Hz, and then the system starts showing secondary harmonic response up to approximately 3.8 Hz (with some transient non-periodic responses). After 3.8 Hz, the primary harmonic solution appears again and stays stable for the rest of the frequencies (as observed from the time series of the pendulum’s response). For this set of parameters, it is clear that the IPVA system not only outperforms the linear benchmark in terms of response suppression but also has higher energy conversion compared with the linear system. This is due to the absence of non- periodic response in the IPVA-PTO system. Thus it should be noted that although non-periodic responses are good at response suppression, it is the secondary harmonic response and rotation that can perform simultaneous response suppression and energy conversion. Next, we perform the numerical analysis of the system near the second resonance frequency. 3.7 Conclusion This study analyzes the modification of the IPVA system [70] with an added generator (referred to as IPVA-PTO) for wave energy conversion. The IPVA-PTO is integrated with a spar, and the dynamics of the system is analyzed. It is observed that a nonlinear energy transfer phenomenon similar to that observed in [70] exists in the IPVA-PTO system too. It is also shown that because of this energy transfer phenomenon, the IPVA-PTO system is shown to outperform its linear counterpart in terms of wave capture width. This analysis is verified by performing experiments 65 on the single-degree-of-freedom “dry” IPVA-PTO system, where the primary harmonic solution is shown to bifurcate via a period-doubling bifurcation. A set of pendulum responses (primary harmonic, secondary harmonic, non-periodic) is observed in the sdof “dry” IPVA-PTO system. Two different load resistance values are used to compare the IPVA-PTO and the linear system, and it is shown that although the non-periodic response of the IPVA-PTO system is only suitable for response suppression, the presence of secondary harmonics can outperform the linear benchmark both in response suppression and energy conversion. 66 CHAPTER 4 INTEGRATION OF THE IPVA WITH A SPAR-FLOATER SYSTEM: A STUDY IN VIBRATION SUPPRESSION 4.1 Overview This chapter talks about the modeling of the spar-floater integrated IPVA system. The system is compared with a linear benchmark for efficacy in wave energy conversion potential and hydrodynamic response suppression. 4.2 Design of the system Figure 4.1 Ocean wave energy harvesting design Figure 4.1(a) shows a spar and an annular floater floating in water. For simiplicity, the spar and floater are constrained such that they can only move in the heaving (𝑥) direction relative to the waterline. Figure 4.1(b) shows the IPVA system consisting of a lead screw, a carrier, and a pendulum vibration absorber. The nut of the lead screw is fixed to the floater while the screw is supported by a thrust bearing that is fixed to the spar through a housing. As a result, the relative heaving displacement 𝑥1 − 𝑥2 is converted into the angular displacement 𝜃 through 𝑥1 − 𝑥2 = 𝑅𝜃, where 𝑅 = 𝐿/2𝜋, and 𝐿 is the screw lead. The carrier is fixed to the screw such that they have the same angular displacement (𝜃). The pendulum pivots on a point of the carrier which is located at a 67 carrierscrewnutpendulum !(b)nut fixed to floaterhousing fixed to spar (c)spar heave "#CGfloater heave "$IPVA(a)wavesparnut motion "#%"$"waterline &’#"$"#IPVAhydrostatic stiffness ($)*,$)*,#’$hydrostatic stiffness (#+-. distance of 𝑅𝑝 from the carrier center. The pendulum has length 𝑟 and an angular displacement 𝜙 with respect to the screw. Figure 4.1(c) shows the mathematical model of the system where 𝑘1 and 𝑘2 denote the hydrostatic stiffness in heave of the spar and the floater, respectively, and 𝑀1 and 𝑀2 denote the mass of the spar and the floater, respectively. Furthermore, the wave motion generates hydrodynamic forces 𝐹𝑔,1 and 𝐹𝑔,2, exciting the spar and the floater, respectively. Linear wave theory is assumed for the current analysis, and the derivation of the hydrodynamic coefficients on the system, along with the equation of motion is discussed in the next section. 4.3 Equations of motion To facilitate deriving the equations of motion, the following coordinate transformation is employed: 𝑥2 = 𝑅𝜓 𝑥1 = 𝑅 (𝜃 + 𝜓) (4.1) Euler-Lagrange mechanics is used to derive the equations of motion. The kinetic and potential energy of the system are determined, followed by the hydrodynamic coefficients of the system. Finally, the virtual work due to the forces applied to the system is derived and the equations of motion are normalized. First, the sum of the kinetic energy of the spar, floater and the IPVA system, and the sum of the potential energy of the spar and floater are given by 𝑇 = + 𝑉 = 1 2 1 2 1 2 (cid:104) (cid:164)𝜃2 (cid:16) 𝐽𝑝 (cid:0) (cid:164)𝜃 + (cid:164)𝜙(cid:1) 2 + 1 2 𝑀1 (cid:0)𝑅 (cid:164)𝜃 + 𝑅 (cid:164)𝜓(cid:1) 2 + 𝑘1(𝑅𝜃 + 𝑅𝜓)2 + 𝐽𝑏𝑠𝑐 (cid:164)𝜃2 + 𝑚 1 2 𝑀2𝑅2 (cid:164)𝜓2, 1 2 1 𝑘2𝑅2𝜓2, 2 𝑟 2 + 2𝑟 𝑅𝑝 cos 𝜙 + 𝑅2 𝑝 (cid:17) + 𝑟 2 (cid:164)𝜙2 + 2𝑟 (cid:164)𝜃 (cid:164)𝜙(𝑟 + 𝑅𝑝 cos 𝜙) (cid:105) (4.2) where 𝑇 and 𝑉 are kinetic and potential energy respectively. Here 𝐽𝑏𝑠𝑐 and 𝐽𝑝 are the moments of inertia of the ball-screw-carrier assembly and of the pendulum respectively. Further, 𝑚 denotes the mass of the pendulum. Next, the hydrodynamic coefficients of the spar and the floater subject to incident regular waves are determined by the linear wave theory, which assumes that the fluid is 68 inviscid and irrotational [74]. Based on the linear wave theory, the hydrodynamic force on the spar and the floater consists of three components: Froude-Krylov force, diffraction force, and radiation force. The first term corresponds to the undisturbed incident wave field without the present of spar-floater system, whereas the diffracted force is the result of modification in the incident wave field due to presence of the spar-floater system, and the radiation force results from the oscillations in the spar-floater system. The Froude-Krylov and diffraction force together give rise to the excitation force while the radiation force gives rise to the added mass and radiation damping [74, 26], which can be represented by the well-known Cummins’ equation [75] 𝐹𝑔,𝑖 = −𝐴∞,𝑖 (cid:165)𝑥𝑖 − ∫ ∞ 𝜎=0 𝑘 𝑅,𝑖 (𝜎) (cid:164)𝑥𝑖 (𝑡 − 𝜎) 𝑑𝜎 + 𝛾𝐹𝑖 (𝑡), 𝑖 = 1, 2, (4.3) where 𝛾 is the wave amplitude, 𝐹𝑔,𝑖 is the incoming wave force, 𝑓𝑖 = 𝛾𝐹𝑖 is the excitation (Froude- Krylov and diffraction) force, and 𝐹𝑖 is the excitation force per wave amplitude. Here 𝐹𝑔,𝑖, and as a consequence 𝐹𝑖, are assumed to be sinusoidal with angular frequency Ω equal to the angular frequency of the incoming wave. The radiation impulse response kernel, 𝑘 𝑅,𝑖 (𝜎) and the radiation infinite-frequency added mass, 𝐴∞,𝑖, are related to the radiation frequency-dependent damping and added mass 𝐵𝑅,𝑖 (Ω) and 𝐴𝑅,𝑖 (Ω), through Ogilvie’s relations [76] 𝐴𝑅,𝑖 (Ω) = 𝐴∞,𝑖 − 1 Ω and ∫ ∞ 𝜎=0 ∫ ∞ 𝜎=0 𝐵𝑅,𝑖 (Ω) = 𝑘 𝑅,𝑖 (𝜎) cos (Ω𝜎) 𝑑𝜎, 𝑖 = 1, 2, 𝑘 𝑅,𝑖 (𝜎) sin (Ω𝜎) 𝑑𝜎, 𝑖 = 1, 2, (4.4) 𝐴∞,𝑖 = lim Ω→∞ 𝐴𝑅,𝑖 (Ω) . (4.5) The hydrodynamic coefficients 𝐴𝑅,𝑖 (Ω) and 𝐵𝑅,𝑖 (Ω), Froude-Krylov and diffraction force 𝑓𝑖 (𝑡) are determined using Ansys AQWA. For the calculation of the hydrodynamic coefficients, the spar-floater system was modeled together but was not coupled physically, though they influenced each other’s hydrodynamics. A convergence test was performed like previously illustrated in [83] 69 Parameter 𝑀1 𝑀2 𝑘1 𝑘2 𝑅 𝐴∞,1 𝐴∞,2 Value 19.46 4.08 198.79 2001.04 0.0672 1.078 12.418 Table 4.1 Parameters for spar-floater system (all data in SI units) Figure 4.2 Ansys AQWA model for calculation of hydrodynamic coefficients to match the published results in [84]. After verifying the Ansys model, the same setting is adopted to simulate the spar-floater system for the rest of the paper. The height of the spar is taken as 1.06 m with a draft of 0.96 m (the height of the spar below the surface of the water), and a diameter of 0.16 m, which is the 1:100 reduced sparD model with water depth of 3 m described in [85]. The floater, on the other hand, is an annulus with a depth of 0.02 m with inner diameter of 0.317 m and outer diameter of 0.595 m. Figure 4.2 shows the mesh of the system used for analysis and Fig. 4.3 shows the added mass, radiation damping, and excitation force obtained using Ansys for both the spar and the floater. Table 4.1 shows the physical properties of both the spar-floater system (obtained via Ansys Aqwa) and the ball-screw transmission chosen for this study. All the data is in SI units. Next, let’s look at the virtual work due to various forces on the system. The total virtual work in the system is 𝛿𝑊 = 𝛿𝑊𝐴𝑀 + 𝛿𝑊𝑅𝐷 + 𝛿𝑊𝐹 + 𝛿𝑊𝐷 (4.6) where 𝛿𝑊𝐴𝑀, 𝛿𝑊𝑅𝐷, 𝛿𝑊𝐹 and 𝛿𝑊𝐷 are the virtual works due to the added mass, radiation damping, 70 a Added mass b Radiation damping c Diffraction Froude-Krylov force and Figure 4.3 Added mass, radiation damping, and diffraction and Froude-Krylov force for spar and floater excitation force and mechanical damping, respectively, and are calculated as follows 𝛿𝑊𝐴𝑀 = −𝐴∞,1 (cid:20)∫ ∞ 𝛿𝑊𝑅𝐷 = − (cid:0) (cid:165)𝜃 + (cid:165)𝜓(cid:1) 𝑅2 (𝛿𝜃 + 𝛿𝜓) − 𝐴∞,2𝑅2 (cid:165)𝜓𝛿𝜓 𝑅2𝜅1 (𝜎) (cid:0) (cid:164)𝜃 (𝑡 − 𝜎) + (cid:164)𝜓 (𝑡 − 𝜎)(cid:1) 𝑑𝜎 (cid:21) (𝛿𝜃 + 𝛿𝜓) 0 (cid:20)∫ ∞ − 𝑅2𝜅2 (𝜎) (cid:164)𝜓 (𝑡 − 𝜎) 𝑑𝜎 (cid:21) 𝛿𝜓 0 𝛿𝑊𝐹 = 𝑅 𝑓1 (𝛿𝜃 + 𝛿𝜓) + 𝑅 𝑓2𝛿𝜓 𝛿𝑊𝐷 = −𝑐𝑅2 (cid:164)𝜃𝛿𝜃 − (cid:0)𝑐 𝑝(cid:1) (cid:164)𝜙𝛿𝜙 (4.7) Using these equations, the 𝜃, 𝜙 and 𝜓 contributions due to virtual work into the equations of motion can be obtained. After substituting the kinetic and potential energy and the virtual work, and using the following normalization parameters 𝜂𝑔 = 𝐽𝑟 𝑀1𝑅2 , 𝜇𝐹 = 𝑀1 𝑀2 , 𝜂 = 𝜏 = 𝜔0𝑡, 𝜉 = 𝑐 2𝜔0𝑀1 , 𝜉𝑝 = 𝑟 𝑅𝑝 , 𝜔0 = 𝑐 𝑝 2𝜔0𝑀1𝑅2 √︂ 𝑘1 𝑀1 , 𝜔2 = √︂ 𝑘2 𝑀2 , 𝜇𝑟 = 𝑚𝑅2 𝑝 𝑀1𝑅2 , 𝜔 = , ()′ = 𝑑 () 𝑑𝜏 , 𝑓𝑒,1 = 𝑓1 𝑀1𝑅 , 𝜔2 0 𝑓𝑒,2 = 𝜔2 0 𝜇 𝐴∞,1 = 𝐴∞,1 𝑀1 , 𝜇 𝐴∞,2 = 𝐴∞,2 𝑀1 the following equations of motion are obtained 𝑀𝑥′′ + 𝐶𝑥′ + 𝐾𝑥 + R (𝑥′) = 𝐹 71 , 𝜔𝑟 = Ω 𝜔0 𝑓2 𝑀1𝑅 , 𝜔2 𝜔0 , (4.8) (4.9) 00.10.20.30.40.50.60.70.80.91Frequency (Hz)051015202530354045Added mass (kg)SparFloater00.10.20.30.40.50.60.70.80.91Frequency (Hz)020406080100120Radiation damping (N-s/m)SparFloater00.10.20.30.40.50.60.70.80.91Frequency (Hz)05001000150020002500Excitation force per wave height (N/m)SparFloater 𝑀 = (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) 2𝜉 0 where 𝐶 = (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) R = 1 𝑀1𝜔2 0 (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) where 𝑎11 1 + 𝜇 𝐴∞,1 𝜇 𝑝 + 𝜂2𝜇𝑟 + 𝜂𝜇𝑟 cos 𝜙 1 + 𝜇 𝐴∞,1 𝜇𝐹 + 1 + 𝜇 𝐴∞,1 + 𝜇 𝐴∞,2 𝜇 𝑝 + 𝜂2𝜇𝑟 + 𝜂𝜇𝑟 cos 𝜙 0 0 𝜇 𝑝 + 𝜂2𝜇𝑟 (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) 0 0 0 0 0 0 2𝜉 𝑝 , 𝐾 = (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) 1 1 0 (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) ∫ ∞ 0 0 1 1 + 𝜇𝐹𝜔2 𝑟 0 , 𝐹 = (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) 𝜅1 (𝜏) [𝜃′ (𝑡 − 𝜏) + 𝜓′ (𝑡 − 𝜏)] 𝑑𝜏 (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) 0 0 𝑓𝑒,1 + 𝑓𝑒,2 −𝜂𝜇𝑟 𝜃′2 sin 𝜙 𝑓𝑒,1 + 2𝜂𝜇𝑟 𝜃′𝜙′ sin 𝜙 + 𝜂𝜇𝑟 𝜙′2 sin 𝜙 (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) ∫ ∞ 0 𝜅1 (𝜏) [𝜃′ (𝑡 − 𝜏) + 𝜓′ (𝑡 − 𝜏)] 𝑑𝜏 + ∫ ∞ 0 𝜅2 (𝜏) 𝜓′ (𝑡 − 𝜏) 𝑑𝜏 0 𝑎11 = 𝜂2𝜇𝑟 + 𝜇𝑏𝑠𝑐 + 𝜇𝑟 + 𝜇 𝑝 + 2𝜂𝜇𝑟 cos 𝜙 + 1 + 𝜇 𝐴∞,1 , 𝑥 = (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) 𝜃 𝜓 𝜙 (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) 4.4 Period-doubling bifurcation in the IPVA system Following the idea laid out in the previous work by the authors [70, 83], the boundary of parametric instability where the primary harmonic solution of the system becomes unstable can be obtained. To that end, the harmonic balance analysis with the modified alternating frequency time (AFT) method and Floquet theory are used as elaborated in [83]. As has been reported previously in [83, 70], the instability in the primary solutions for the system due to period-doubling bifurcations gives rise to additional terms in the ball-screw (𝜃) and pendulum’s response (𝜙), referred to as secondary solutions (with harmonics of frequency 𝜔 2 where 𝜔 denotes the wave excitation frequency) in the system. This period-doubling bifurcation gives rise to 1:2 internal resonance in the system. We use the bifurcation tracking algorithm developed in [70] with modified AFT method proposed in [83] to obtain the period-doubling bifurcation boundary. The boundary in Fig. 4.4 shows the wave height on the y-axis and wave frequency on the x-axis. Like the results in [83, 70], below this boundary, the primary harmonic solution is stable and above this boundary, the period of the solution doubles. To verify this behavior, the hydrodynamic 72 response with three wave heights (×1, ×2, and ×3, ×3 not shown), corresponding to 4.5 cm, 5.4 cm and 12 cm at a frequency of 0.44 Hz, is simulated by using the explicit Runge-Kutta (4,5) integration method. The integration kernel is evaluated using an impulse to state-space converter function described in [80, 83]. Fig. 4.4, Fig. 4.4 and Fig. 4.4 show the response of the pendulum (𝜙) and the fast Fourier transform (FFT) of the response of the pendulum, and the ball-screw system (𝜃). Note that ˆ𝜔 represents the frequency of oscillation normalized with respect to the wave frequency, which means that ˆ𝜔=1 corresponds to the primary component of the solution, and ˆ𝜔 = 1 oscillation of angular frequency 𝜔 2 represents 2 , 𝜔 being the wave excitation frequency. As can be observed from these figures, below the stability boundary the solution is strictly primary, whereas at point ×2, the motion of the pendulum (along with the ball-screw) contains harmonics of frequency 𝜔 2 , in addition to the primary harmonics. If we further increase the wave height to point ×3, the time series of the pendulum’s motion shows intermittent rotation and oscillation, whereas the FFT shows the motion of the pendulum (and the ball-screw) consisting of many frequencies, showing non-periodic behavior. This observation has been recorded previously by the authors in [70], and the response can be attributed to a cascade of period-doubling bifurcation which eventually can lead to chaotic-like motions in the system, as evident from response at point 3. It has been noted previously in this work that the crossing of the period-doubling bifurcation boundary results in energy transfer between the primary (spar-floater) system and the pendulum vibration absorber. To demonstrate this, we fix a wave height, say 5.5 cm (which is above the instability boundary for a frequency range from approximately 0.42 Hz to 0.465 Hz, as marked by vertical dashed lines in Fig. 4.4), and look at the frequency response of the proposed system. 4.5 Nonlinear energy transfer and energy harvesting potential As mentioned in the previous section, for the wave height of 5.5 cm, we expect an energy transfer between the spar-floater system and the pendulum for the specified frequency range. To benchmark this system, we use a linear system defined as the system with the pendulum locked such that its first natural frequency matches the resonant frequency of the IPVA system (given by the frequency corresponding to the lowest wave height in the stability boundary; see Fig. 4.4). We 73 Figure 4.4 Stability boundary for the primary harmonics of the system for 𝜇𝑟 = 0.2, 𝜂 = 0.3, 𝜉 = 0.05, 𝜉𝑝 = 0.02, 𝜇𝑝 = 0.039, 𝜇𝑏𝑠𝑐 = 0.1 a Time series of the pendulum motion at point ×1 b FFT of the pendulum motion at point ×1 Figure 4.5 Time series and FFT of the pendulum motion at point ×1 perform direct numerical simulation for both the IPVA and the linear system using the explicit Runge-Kutta (4,5) integration method (implemented in Matlab’s ode45), with integration kernel evaluated using an impulse to state-space converter [80, 83]. Fixing the wave height to 5.5 cm, we calculate the response amplitude operator (RAO) of the spar for both the IPVA and the linear system, defined by the response of the spar divided by the wave height. Figure 4.8a shows the comparison of the RAO between the IPVA and the linear system. The range for the secondary resonance is marked by vertical dash-dot lines. Figure 4.8b shows the values of motion of the 74 0.40.420.440.460.48Frequency (Hz)0.0450.050.0550.060.0650.070.0750.08Wave height (m)1256405660568057005720574057605780-1-0.500.511.522.500.511.522.5310-610-410-2100FFTFFT of FFT of a Time series of the pendulum motion at point ×2 b FFT of the pendulum motion at point ×2 Figure 4.6 Time series and FFT of the pendulum motion at point ×2 a Time series of the pendulum motion at point ×3 b FFT of the pendulum motion at point ×3 Figure 4.7 Time series and FFT of the pendulum motion at point ×3 ball-screw ( (cid:164)𝜃2) and the pendulum ( (cid:164)𝜙2) of the IPVA system, and the motion of the ball-screw ( (cid:164)𝜃2) for the linear system. In the literature, ball-screw and pendulum angular motions are used to drive electrical generators for wave energy production and the electrical power is proportional to the angular velocities squared, [86, 83, 87], for example. Therefore, the angular velocities squared are used to examine the energy conversion potential in this work. As can be seen in Fig. 4.8a and Fig. 4.8b, the linear system sacrifices the spar’s response for the resonant ball-screw angular motion. The IPVA system, on the other hand, has a significant smaller RAO in comparison with the linear system. It can be readily seen that till around the frequency 0.4 Hz, the ball-screw motion of the linear system and the IPVA system closely follow each other 75 56405660568057005720574057605780-1012300.511.522.5310-410-2100FFTFFT of FFT of 0.80.911.11.2104-200-150-100-5000.511.522.5310-610-410-2100102104FFTFFT of FFT of a Comparison between motion of spar for IPVA and linear system for the parameters mentioned in Fig. 4.4 with 𝜉 𝑝 = 0.02 b Comparison of energies in the IPVA system for the parameters mentioned in Fig. 4.4 with 𝜉 𝑝 = 0.02 c Comparison between motion of spar for IPVA and linear system for the parameters mentioned in Fig. 4.4 with 𝜉 𝑝 = 0.01 d Comparison of energies in the IPVA system for the parameters mentioned in Fig. 4.4 with 𝜉 𝑝 = 0.01 Figure 4.8 The frequency responses of the linear and the IPVA system in amplitude. However, the solution starts to deviate around 0.4 Hz due to increase in motion of the pendulum, and soon around 0.42 Hz, in the range corresponding to the secondary solutions (marked by dash-dot lines in the figure) in the IPVA system, the energy transfer from the primary system to the pendulum begins. Hence, it is evident that the energy in the ball-screw transfers to the pendulum due to internal resonance, resulting in suppression of the spar motion and pumping of the pendulum’s kinetic energy. From Fig. 4.8b, it can be clearly seen that the motion of the linear system is better than the IPVA system if the we consider the motion of the pendulum alone for 76 0.350.40.450.5Frequency (Hz)00.511.522.5RAO of the sparIPVA-PTOLinear0.350.40.450.5Frequency (Hz)00.511.522.53Comparison of different angular motions0.350.40.450.5Frequency (Hz)00.511.522.5RAO of the spar IPVA-PTOLinear0.350.40.450.5Frequency (Hz)00.20.40.60.811.2Comparison of different angular motions the IPVA system. However, if we use the difference between the angular motion of the pendulum and the ball-screw, the motion of the IPVA system is significantly higher than the linear system. This is because it has been observed that the ball-screw and the pendulum motion are generally out of phase for the IPVA system [83, 70]. A mechanism that converts the difference of angular motions into one angular motion to drive a generator can be readily achieved using a planetary gear setup; see [83]. This mechanism justifies the use of the angular motion difference to examine the energy conversion potential. Note that since the pendulum is locked for the linear system, there is effectively no pendulum damping in the linear system’s equation of motion. Therefore, it is worth investigating how the pendulum damping can affect the angular motions in the IPVA system (while the linear system stays the same). To study the effect of pendulum’s damping on the performance of the IPVA system, we calculate the 1:2 internal resonance or period-doubling bifurcation boundary for the IPVA system with a different damping value, taken as 𝜉 𝑝 = 0.01. This boundary is shown in Fig. 4.9d as a part of parametric studies (to be discussed later). It can be seen that the stability boundary for the case of 𝜉 𝑝 = 0.01 is significantly lower than 𝜉 𝑝 = 0.02 shown in Fig. 4.4. Therefore, we choose a wave height of 2.5cm to perform the numerical frequency response analysis for the case of 𝜉 𝑝 = 0.01. The frequency response is shown in Fig. 4.8c and Fig. 4.8d. We can observe a similar energy transfer phenomenon and vibration suppression in the spar’s RAO. However in this case, motion value of the IPVA defined in terms of (cid:0) (cid:164)𝜃 − (cid:164)𝜙(cid:1) 2 is around 3 times higher than the linear system compared to the case of 𝜉 𝑝 = 0.02, where the same motion is around 1.5 times higher. Another thing worth mentioning is that the RAO of the spar for 𝜉 𝑝 = 0.01 is higher than that of 𝜉 𝑝 = 0.02, though it still outperforms the linear system. Next, let us look at the effects of system parameters on the stability boundary. 4.6 Parametric studies To analyze the effects of various system parameters on the stability boundary of internal resonance (period-doubling bifurcation), we vary the following parameters 𝜇𝑟, 𝜂, 𝜉, and 𝜉 𝑝 while keeping the other parameters fixed at their values mentioned in Fig. 4.4. First, we see the effect 77 of 𝜇𝑟 on the stability boundary. Recall that 𝜇𝑟 is the mass amplification factor in the system given by 𝜇𝑟 = 𝑚𝑅2 𝑝 𝑀 𝑅2 . As the 𝜇𝑟 value increases, we see the wave height required to cross the internal boundary decreases as evident from Fig. 4.9a. Therefore to control the wave height required to cross the boundary, one can readily change the value of 𝑅 𝑝 𝑅 , which is the ratio of the distance of the pendulum pivot point with respect to the carrier over the effective radius of the ball-screw system. Next, the effect of 𝜂, defined by 𝑟 𝑅 𝑝 where 𝑟 is the length of the pendulum, on stability boundary is analyzed. As observed from Fig. 4.9b, for larger values of 𝜂, the wave height required to cross the boundary for a given frequency decreases. Fig. 4.9c and Fig. 4.9d show the effect of mechanical and pendulum damping respectively on the stability boundary. It can be observed that the wave height required to reach 1:2 internal resonance increases with increase in both the mechanical and pendulum damping value. The next section discusses the experimental verification for the IPVA system. 4.7 Conclusion This study analyzes the incorporation of the IPVA system [70] into a heaving spar-floater system to study the energy transfer between the spar-floater system and the pendulum vibration absorber. The hydrodynamic response and wave energy conversion potential of the integrated system, when the wave frequency is near the spar resonance frequency, are investigated using numerical frequency response simulations. A harmonic balance method is used to determine the boundary of period- doubling bifurcation in the parameter plane of wave height and wave frequency. According to the boundary, one can determine a combination of the wave height and frequency such that 1:2 internal resonance occurs to the IPVA system. It is observed that this 1:2 internal resonance is associated with a nonlinear energy transfer phenomenon similar to that observed in [70, 83]. It is also shown that because of this energy transfer phenomenon, the IPVA-PTO system achieves a lower maximum RAO, and a higher energy transfer potential compared to the linear benchmark, when the relative angular motion between the ball-screw and the pendulum is used as a measure of energy conversion potential. The effect of the pendulum damping is also characterized in the IPVA system, since it is the only parameter which is missing from the linear benchmark. It is found that for lower 78 a Stability boundaries for different values of 𝜇𝑟 b Stability boundaries for different values of 𝜂 c Stability boundaries for different values of 𝜉 d Stability boundaries for different values of 𝜉 𝑝 Figure 4.9 Effect of various system parameters on the stability boundary, all the values are from Fig. 4.4 except for the parameter value specified in the figure pendulum damping, a lower wave height is required to achieve the energy transfer phenomenon and the angular motions of the IPVA system are significantly higher than the linear system, although the RAO value also increases compared to the higher pendulum damping case. Finally, parametric studies showing the effects of various system parameters on the 1:2 internal resonance boundary were also shown in this study. 79 0.380.40.420.440.460.480.5Frequency0.040.060.080.10.120.140.160.18Wave height (m)r=0.15r=0.2r=0.250.380.40.420.440.460.480.5Frequency0.040.060.080.10.120.140.16Wave height (m)=0.25=0.3=0.350.380.40.420.440.460.480.5Frequency0.020.040.060.080.10.120.140.16Wave height (m)=0.03=0.05=0.070.380.40.420.440.460.480.5Frequency0.020.040.060.080.1Wave height (m)p=0.005p=0.01p=0.02 CHAPTER 5 INTEGRATION OF THE IPVA-PTO WITH A SPAR-FLOATER SYSTEM: STUDY IN OCEAN WAVE ENERGY CONVERSION 5.1 Overview This chapter talks about the modeling of the spar-floater integrated IPVA-PTO system. The system is compared with a linear benchmark for efficacy in wave energy conversion and hydrodynamic response suppression near both its resonance frequencies. A discussion on optimization of electrical damping for maximizing the wave energy conversion is discussed. 5.2 Modeling of the system In this section, first, the operation principle of the system is discussed, and then the equations of motion of the system are derived, taking into account the hydrodynamic effects. 5.2.1 Operation principle of the system Figure 5.1 The IPVA-PTO system integrated in a spar-floater setup for wave energy conversion where: a. Ocean wave energy conversion setup, b. the IPVA-PTO system, and c. equivalent mathematical model Figure 5.1a shows a spar and an annular floater floating in water connected by the IPVA-PTO system. For simplicity, the spar and floater are constrained such that they can only move in the heaving (𝑥) direction relative to the waterline. Figure 5.1b shows the IPVA system consisting of a lead screw, a carrier, a generator, and a pendulum vibration absorber. The nut of the lead screw is fixed to the floater while the screw is supported by a thrust bearing that is fixed to the spar through a housing. As a result, the relative heaving displacement 𝑥1 − 𝑥2 is converted into the angular 80 displacement 𝜃 through 𝑥1 − 𝑥2 = 𝑅𝜃, where 𝑅 = 𝐿/2𝜋 and 𝐿 is the screw lead. The carrier is fixed to the screw such that they have the same angular displacement (𝜃). The pendulum pivots at a distance of 𝑅𝑝 from the carrier center, and has a mass of 𝑚 and a moment of inertia 𝐽𝑝 with respect to the center of mass. Furthermore, it has a length of 𝑟 and an angular displacement 𝜙 with respect to the screw. A sun gear (mounted on the carrier via bearings) hosts the generator rotor while the generator housing is fixed to the spar. A planet gear is mounted to the pendulum shaft and meshes with the sun gear. It can be shown that the rotational motion of the generator rotor with respect to its housing is given by 𝜃 − 𝜙. Therefore, the instantaneous energy converted by the generator will be 𝑐𝑒 (cid:0) (cid:164)𝜃 − (cid:164)𝜙(cid:1) 2, where 𝑐𝑒 denotes the electrical damping of the generator. Figure 5.1c shows the mathematical model of the system, where 𝑘1 and 𝑘2 denote the hydrostatic stiffness of the spar and floater in the heaving direction, respectively, and 𝑀1 and 𝑀2 denote the mass of the spar and the floater, respectively. Furthermore, the wave motion generates hydrodynamic forces 𝑓1 and 𝑓2, exciting the spar and floater, respectively. For the analysis in this study, linear wave theory is assumed, and the derivation of the equation of motion is discussed next. 𝜓 = 𝑥2/𝑅, 𝜃 = (𝑥1 − 𝑥2) /𝑅 (5.1) 5.2.2 Equations of motion The equations of motion for a similar system (without the generator, sun gear and planet gears) have been previously derived in [88]. Linear wave theory is considered in the derivation of the equation of motion, where the fluid is assumed to be inviscid and irrotational [75]. Several important relations from [88] are provided below for reference. First, the following coordinate transformation is employed to facilitate the derivation: Next, the hydrodynamic forces on the system consist of three components: the radiation force, the incident (Froude-Krylov) force, and the diffraction force. The Froude-Krylov and diffraction forces determine the excitation force [74, 26] on the spar and the floater, and together with the radiation force, they can be represented by the well-known Cummins’ 81 a Ansys Aqwa model mesh b Added mass c Radiation damping d Diffraction and Froude-Krylov force Figure 5.2 Ansys Aqwa mesh model, added mass, radiation damping, and diffraction and Froude- Krylov force for spar and floater equation [75]: 𝐹𝑔,𝑖 = −𝐴∞,𝑖 (cid:165)𝑥𝑖 − ∫ ∞ 𝜎=0 𝜅𝑖 (𝜎) (cid:164)𝑥𝑖 (𝑡 − 𝜎) 𝑑𝜎 + 𝛾𝐹𝑖 (𝑡), 𝑖 = 1, 2 (5.2) where 𝐹𝑔,𝑖 is the incoming wave force, 𝛾𝐹𝑖 = 𝑓𝑖 is the excitation (Froude-Krylov and diffraction) force, 𝛾 is the wave amplitude, 𝜅𝑖 (𝜎) is the radiation impulse response kernel, and 𝐴∞,𝑖 represents the radiation infinite-frequency added mass. These hydrodynamic forces and coefficients of the system are simulated using Ansys Aqwa. The details can be found in [88]. Fig. 5.2 shows the mesh of the system, the simulated added mass, radiation damping, and Froude-Krylov and diffraction force on the spar and the floater. Note that in Fig. 5.2c, the spar’s radiation damping is insignificant compared to the floater’s, and therefore is not visible in the figure. By virtual of linear wave theory, 82 00.511.5/2 = Frequency (Hz)01020304050Added mass (N-s/m)SparFloater00.511.5/2 = Frequency (Hz)050100150200250300350Radiation damping (N-s/m)SparFloater00.511.5/2 = Frequency (Hz)05001000150020002500Excitation force per wave height (N/m)SparFloater 𝐹𝑔,𝑖 is assumed to be sinusoidal. As a consequence, 𝐹𝑖 is also sinusoidal with angular frequency Ω, i.e., the frequency of the incoming wave. Finally, 𝜅𝑖 (𝜎) and 𝐴∞,𝑖 are related to the radiation frequency-dependent damping and added mass 𝐵𝑅,𝑖 (Ω) and 𝐴𝑅,𝑖 (Ω) through Ogilvie’s relations [76, 74]. The latter are shown in Fig. 5.2b and Fig. 5.2c. After considering the contributions of kinetic energy, potential energy, virtual work, and hydrodynamic effects on the system, the following equation of motion is obtained Mx′′ + Cx′ + Kx + R(x′) = F, (5.3) with where M = K = F = (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) 2 (𝜉 + 𝜉𝑒) 0 −2𝜉𝑒 𝑎11 𝑎12 𝑎13 𝑎12 𝑎22 𝑎13 0 0 𝑎33 , C = (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) 0 0 −2𝜉𝑒 1 1 𝑟 0 1 1 + 𝜇𝐹𝜔2 , R(x′) = (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) 𝑓𝑒,1 + 2𝜂𝜇𝑟 𝜃′𝜙′ sin 𝜙 + 𝜂𝜇𝑟 𝜙′2 sin 𝜙 1 𝑀1𝜔2 0 0 0 0 𝑓𝑒,1 + 𝑓𝑒,2 −𝜂𝜇𝑟 𝜃′2 sin 𝜙 , (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) 0 2𝜉𝑒 𝑟1 𝑟1 + 𝑟2 0 0 0 (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) , x = (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) 𝜃 𝜓 𝜙 , (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) (5.4) 𝑎11 = 1 + (cid:16) 1 + 𝜂2(cid:17) 𝜇𝑟 + 𝜇𝑏𝑠𝑐 + 𝜇 𝑝 + 2𝜂𝜇𝑟 cos 𝜙 + 𝜇𝑔 + 𝜇 𝐴∞,1 𝑎12 = 1 + 𝜇 𝐴∞,1, 𝑎13 = −𝜇𝑔 + 𝜇 𝑝 + 𝜂2𝜇𝑟 + 𝜂𝜇𝑟 cos 𝜙, 𝑎22 = 1 + 𝜇𝐹 + 𝜇 𝐴∞,1 + 𝜇 𝐴∞,2, 𝑎33 = 𝜇 𝑝 + 𝜂2𝜇𝑟 + 𝜇𝑔, 𝑟1 = 𝑟2 = ∫ ∞ 0 ∫ ∞ 0 𝜅1 𝜅2 (cid:19) (cid:19) (cid:18) 𝑠 𝜔0 (cid:18) 𝑠 𝜔0 [𝜃′ (𝜏 − 𝑠) + 𝜓′ (𝜏 − 𝑠)] 𝑑𝑠, 𝜓′ (𝜏 − 𝑠) 𝑑𝑠, (5.5) 83 and 𝜇𝑔 = 𝜔2 = , 𝜇𝐹 = , 𝜇𝑟 = , 𝜇𝑏𝑠𝑐 = , 𝜂 = 𝑟 𝑅𝑝 𝑀1 𝑀2 𝑚𝑅2 𝑝 𝑀1𝑅2 , 𝜏 = 𝜔0𝑡, 𝜉 = 𝑑 () 𝑑𝜏 , 𝑓𝑒,1 = 𝜔2 𝜔0 , ()′ = 𝐽𝑟 𝑀1𝑅2 √︂ 𝑘2 𝑀2 , 𝜔𝑟 = Ω 𝜔0 𝑐𝑒 2𝜔0𝑀1𝑅2 𝑓2 𝑀1𝑅 𝜔2 0 𝐴∞,1 𝑀1 , 𝜔0 = √︂ 𝑘1 𝑀1 , 𝐽𝑏𝑠𝑐 𝑀1𝑅2 𝑐 2𝜔0𝑀1 𝑓1 𝑀1𝑅 𝜔2 0 , , 𝑓𝑒,2 = , 𝜇 𝐴∞,1 = , 𝜇 𝐴∞,2 = 𝜔 = 𝜉𝑒 = , 𝜇𝑝 = 𝐽𝑝 𝑀1𝑅2 , 𝐴∞,2 𝑀1 , (5.6) In (5.6), 𝐽𝑟 and 𝐽𝑏𝑠𝑐 are the moment of inertia of the generator’s rotor and ball-screw-carrier assembly, respectively. Also, 𝑐𝑒 is the electrical damping due to the generator, and 𝑐 is the mechanical damping between the spar and the floater. The damping in the pendulum is assumed to be significantly less than the electrical damping and is ignored in the current study. 5.3 Bifurcation analysis of the IPVA-PTO system From our previous studies [70, 89, 83, 81], it has been shown that a period-doubling bifurcation of the primary harmonic solution (corresponding to 1:1 resonance) results in a nonlinear energy transfer between the primary system and the pendulum vibration absorber. Therefore, we look for the bifurcation boundary that determines when a bifurcation occurs to the primary harmonic solution of the current system. To that end, the harmonic balance analysis, in conjunction with the alternating frequency time method (AFT) and Floquet theory are used to find the bifurcation boundary as explained in [83]. Since the spar-floater system (without considering the IPVA-PTO) has two degrees of freedom, we expect two resonance frequencies for the system (as the pendulum does not have any natural frequency of its own), where the first and second resonance mode are dominated by the spar and floater, respectively. Thus, we will determine the bifurcation boundary within a frequency range that includes the two resonance frequencies. It was observed in [83, 70, 81, 89] that when the period-doubling bifurcation occurs, the secondary harmonic solution (with harmonics of frequency 𝜔 2 , where 𝜔 is the excitation frequency) emerges in the response. The secondary harmonic solution is attributed to 1:2 internal resonance 84 in the system and is accompanied by nonlinear energy transfer between the primary system and the pendulum. However, for the current system, we see two types of bifurcations — period- doubling and secondary Hopf bifurcation. Therefore, the bifurcation tracking algorithm developed in [70, 83], with modifications for tracking the secondary Hopf bifurcation, is developed in this work to obtain the bifurcation boundary of the current system (which contains the period-doubling and the secondary Hopf bifurcation boundary). Fig. 5.3 shows the bifurcation boundary with wave height on the ordinate and wave frequency on the abscissa for parameters listed in its caption. Several things are worth noting in Fig. 5.3. First, the two resonance frequencies each correspond to the two local minima of the wave height marked by blue bullet marker (•) in the figure. Next, there are two types of bifurcations via which the primary solution becomes unstable: the solid line represents the period-doubling bifurcation whereas the dashed line represents the secondary Hopf bifurcation. This is verified by examining the real and imaginary parts of the eigenvalues of the fundamental solution matrix along the boundary. The eigenvalues whose magnitude equals one are plotted in Fig. 5.4. As shown, till around the frequency of 0.8 Hz, it is found that the eigenvalues are real and equal to −1. This corresponds to the period-doubling bifurcation. Between 0.8 Hz and 1.1 Hz, the imaginary part of the eigenvalues become nonzero and the eigenvalues become complex conjugates, corresponding to the secondary Hopf bifurcation. There is a region between 1.1 Hz and 1.35 Hz in which both the period-doubling and secondary Hopf bifurcation boundaries co-exist, after which the secondary Hopf bifurcation boundary disappears. Furthermore, we investigate two points ×1 and ×2 in Fig. 5.3 above the calculated secondary Hopf bifurcation boundary, which should show a quasi-periodic response. Fig. 5.5 shows that the projection of the Poincare section in the plane of the pendulum’s angular displacement (𝜙) and pendulum’s velocity (𝜙′) is a closed curve for both ×1 and ×2. This means that the system response is quasi-periodic; hence the secondary Hopf bifurcation boundary is verified (period-doubling bifurcation boundary is verified previously in [89] and not shown for the sake of brevity). Now from Fig. 5.3, we see different types of responses above the bifurcation boundary, such as secondary harmonic solutions, non-periodic (chaotic-like) oscillations, rotation, quasi-periodic 85 Figure 5.3 Bifurcation boundary for the primary harmonics of the system for 𝜉 = 0.01, 𝜉𝑒 = 0.015, 𝜉𝑒=0.02, 𝜇𝑟=0.4, 𝜂=0.4, 𝜇 𝐴∞,1=0.0554, 𝜇 𝐴∞,2=0.6381 , 𝜔𝑟=6.928, 𝜔0=3.20 rad/s, 𝜇𝑔=0.01, 𝜇 𝑝=0.02, 𝜇𝑏𝑠𝑐=0.04, 𝜇𝐹=0.2097. Markers represent the following: Primary — ‘+’, Chaotic — ‘*’, One-fifth subharmonics — ‘★’, Secondary — ‘□’, Rotation — ‘◦’, Intermittent rotation — ‘△’ response as indicated by the markers (see the caption). For the current analysis, it is our assumption that the period-doubling bifurcation will result in energy transfer between the primary system and the pendulum as evident from various studies on autoparametric systems [51]. On the other hand, how the secondary Hopf bifurcation influences the system performance remains unexplored. To verify the assumption and characterize the secondary Hopf bifurcation, we investigate the response of the system near the first and second resonance frequencies for various wave heights, and mechanical and electrical damping values in the following sections. 5.4 Performance analysis near first resonance frequency Near the first resonance frequency (defined as the interval between 0.25 Hz and 0.5 Hz), we see that the primary harmonic solution undergoes the period-doubling bifurcation (giving rise to secondary harmonic solutions). To show the effect of different electrical damping on the energy transfer phenomenon, we choose two electrical damping values: 𝜉𝑒 = 0.01 and 𝜉𝑒 = 0.02, and the mechanical damping value of 𝜉 = 0.02, while keeping all other parameters the same as mentioned in Fig. 5.3. The period-doubling bifurcation boundary (henceforth referred to as the PD boundary) 86 0.511.5Frequency (Hz)0.020.030.040.050.060.070.080.090.10.110.12Wave height (m)PDQP12 a Real part of eigenvalues with magnitude 1 b Imaginary part of eigenvalues with magnitude 1 Figure 5.4 Real and imaginary parts of eigenvalues whose magnitude equals one with respect to frequency, other eigenvalues not shown a Poincare section at point 1 b Poincare section at point 2 Figure 5.5 Poincare sections corresponding to points 1 and 2 in Fig. 5.3 show that the secondary Hopf bifurcation results in quasi-periodic motion for both cases is shown in Fig. 5.6. For the frequency range within which the wave height is above the PD boundary, nonlinear energy transfer is expected between the spar-floater system and the pendulum as previously reported by the authors in [70, 89, 83, 81]. Thus, to benchmark the current system, we define a linear system and compare two measures: the response amplitude operator (RAO) of the spar and the normalized capture width with respect to frequency. The RAO of the spar is defined as the heaving amplitude of the spar per unit wave height. On the other hand, the normalized capture width is defined as the value of the energy converted by the IPVA-PTO or the linear system divided by the maximum energy converted by the linear system. For a fair 87 0.20.40.60.811.21.41.6Frequency (Hz)-1.1-1-0.9-0.8-0.7-0.6-0.5-0.4-0.3Real part of eigenvalues0.20.40.60.811.21.41.6Frequency (Hz)-1-0.8-0.6-0.4-0.200.20.40.60.81Imaginary part of eigenvalues1.581.61.621.641.661.681.71.721.74-0.1-0.0500.050.10.150.20.250.30.354.84.955.15.25.35.45.50.250.30.350.40.450.50.550.6 a Bifurcation boundary for 𝜉 = 0.02, 𝜉𝑒 = 0.01 b Bifurcation boundary for 𝜉 = 0.02, 𝜉𝑒 = 0.02 Figure 5.6 Period-doubling bifurcation boundary for two values of electrical damping with mechanical damping 𝜉 = 0.02. All other parameters are from Fig. 5.3 comparison, we define the linear system as the system with the pendulum locked at a position such that the frequency value corresponding to the maximum capture width of the linear system matches the first resonance frequency of the IPVA-PTO system. Numerical integration of the equations of motion for both the linear and the IPVA-PTO system is performed using Matlab’s “ode45” function. The integration kernel in the equation of motion is evaluated using an impulse to state-space converter function described in [83, 80]. Finally, optimal electrical damping in the linear system, denoted by 𝜉𝑒,𝑙, is calculated in order to maximize the power converted from the linear benchmark. To calculate the optimal linear electrical damping, a constraint is enforced on the electrical damping, which limits its maximum value to 0.07. This constraint reflects the fact that commercially available generators whose rated power is close to our simulation results have a similar limit on the electrical damping, as confirmed in our repeated market surveys. Now, for the first simulation, we use the parameters used to obtain Fig. 5.6a, which gives 𝜉𝑒,𝑙 = 0.0216 for the linear system (calculated using the nonlinear optimization solver in Matlab). We simulate the system for various wave heights of 0.5, 1.5, and 2.5 cm. As can be seen from Fig. 5.6a, at the resonance frequency (the frequency corresponding to the minimum wave height of the boundary), the wave height 0.5 cm is below, 1.5 cm is slightly above, and 2.5 cm is significantly above the PD boundary. Fig. 5.7 shows the RAO and normalized capture width obtained. Several 88 0.320.340.360.380.40.420.440.460.480.5/2 = Frequency (Hz)0.010.020.030.040.050.060.070.080.09Wave height (m)0.380.40.420.440.460.480.5/2 = Frequency (Hz)0.040.0450.050.0550.060.0650.07Wave height (m) things can be observed in Fig. 5.7. First, for the wave height of 0.5 cm, the maximum RAO of the IPVA-PTO system is worse than the linear system, although it outperforms the linear system in terms of the normalized capture width. Next, at the wave height of 1.5 cm, we start to see energy transfer between the primary system (spar-floater combination) and the IPVA-PTO system, as evidenced by the widening of the half-power bandwidth in the normalized capture width of the IPVA-PTO (while the maximum power is almost equal). The RAO suppression is better than the wave height 0.5 cm but is still worse than the linear system. This can be attributed to the fact that the optimal electrical damping of the linear system is 𝜉𝑒,𝑙 = 0.0216, which is almost twice the electrical damping in the IPVA-PTO system (𝜉𝑒 = 0.01). Finally, at the wave height of 2.5 cm, we see the pendulum and the spar in the IPVA-PTO system have non-periodic motion for frequencies greater than 0.37 Hz to 0.47 Hz (from time series analysis, not shown). In this case, the IPVA-PTO system outperforms the linear system in terms of both RAO suppression and normalized capture width. We also see the half-power bandwidth of the normalized capture width increase for the IPVA-PTO system compared with the previous wave heights, showing the nonlinear energy transfer phenomenon in effect. For the second simulation, the parameters used to obtain Fig. 5.6b are used for analysis, which again gives optimal 𝜉𝑒,𝑙 = 0.0216 for the linear system. Note that due to the increase in the electrical damping, a larger wave height is required to cross the PD boundary. Therefore, for this case, the system is simulated for wave heights of 3, 4.5, and 6 cm, which is below, slightly above, and significantly above the PD boundary, respectively at the resonance frequency. Fig. 5.8 shows the RAO and normalized capture width obtained. Here it can be observed that the IPVA-PTO system outperforms the linear system for all three wave heights in terms of the maximum RAO and normalized captured width. This can be contrasted with the case of 𝜉𝑒 = 0.01 as follows. First, the electrical damping of the IPVA-PTO system (𝜉𝑒 = 0.02) is larger than the previous case (𝜉𝑒 = 0.01), and therefore has better RAO response compared to the case for 𝜉𝑒 = 0.01. Next, the maximum normalized capture width for the IPVA-PTO system with 𝜉𝑒 = 0.02 is smaller than the case when 𝜉𝑒 = 0.01. Furthermore, the IPVA-PTO system’s response is little influenced by crossing the PD 89 bifurcation boundary, which suggests that the energy transfer phenomenon is insignificant when the electrical damping is too large. a RAO. Readers are referred to the web version of this article for clarity b Normalized capture width Figure 5.7 Comparison between RAO and normalized capture width of the IPVA-PTO and the linear system for different wave heights for 𝜉 = 0.02, 𝜉𝑒,𝑙 = 0.0216 and 𝜉𝑒 = 0.01, all other parameters from caption of Fig. 5.3 a RAO. Readers are referred to the web version of this article for clarity b Normalized capture width Figure 5.8 Comparison between RAO and normalized capture width of the IPVA-PTO and the linear system for different wave heights for 𝜉 = 0.02, 𝜉𝑒,𝑙 = 0.0216 and 𝜉𝑒 = 0.02, all other parameters from caption of Fig. 5.3 5.5 Performance analysis near second resonance frequency From Fig. 5.3, it is expected that for the frequency range near the second resonance (defined as the range from 0.8 Hz to 1.1 Hz for the referred figure), we see that the primary harmonic response undergoes secondary-Hopf bifurcation, giving rise to quasi-periodic solutions. The quasi-periodic 90 0.250.30.350.40.450.5/2 = Frequency (Hz)0.20.40.60.811.21.41.6RAOIPVA-PTO 2.5IPVA-PTO 1.5IPVA-PTO 0.5Linear0.250.30.350.40.450.5/2 = Frequency (Hz)00.511.522.53Normalized capture widthIPVA-PTO 2.5IPVA-PTO 1.5IPVA-PTO 0.5Linear0.250.30.350.40.450.5/2 = Frequency (Hz)0.20.40.60.811.21.41.6RAOIPVA-PTO 6IPVA-PTO 4.5IPVA-PTO 3Linear0.250.30.350.40.450.5/2 = Frequency (Hz)00.20.40.60.811.21.41.61.8Normalized capture widthIPVA-PTO 6IPVA-PTO 4.5IPVA-PTO 3Linear response then becomes non-periodic and eventually rotational as wave height increases, as evident from the markers in Fig. 5.3. To benchmark the IPVA-PTO system for this case too, the linear system is defined as the system with the pendulum locked at a position such that the frequency value corresponding to the maximum capture width for the linear system matches the second resonance frequency of the IPVA-PTO system. The optimal electrical damping for the linear system is determined (denoted by 𝜉𝑒,𝑙), to maximize the power converted from the linear benchmark with the constraint 𝜉𝑒 ≤ 0.07, similar to the analysis near the first resonance frequency. For this case, we choose two different electrical damping values 𝜉𝑒 = 0.015 and 𝜉𝑒 = 0.02 while keeping 𝜉 = 0.01. The corresponding bifurcation boundaries can be found in Fig. 5.3 and Fig. 5.11d. a RAO. Readers are referred to the web version of this article for clarity b Normalized capture width Figure 5.9 Comparison between RAO and normalized capture width of the IPVA-PTO and the linear system for different wave heights for 𝜉𝑒,𝑙 = 0.07, and the boundary corresponding to Fig. 5.3. All marked wave heights in cm Now, for the first simulation, we take 𝜉 = 0.01, 𝜉𝑒 = 0.015 for the IPVA-PTO system, with the rest of the parameters taken from the caption of Fig. 5.3, which gives 𝜉𝑒,𝑙 = 0.07 for the linear system. Note that the optimal electrical damping would be 1.6583 without the constraint, which would be unattainable by commercially available generators. We simulate the system for various wave heights of 3, 3.5, and 5.5 cm. The results are shown in Fig. 5.9. As can be seen from Fig. 5.3, all the chosen wave heights are above the PD boundary near the first resonance frequency. Therefore, we can observe response suppression at the first resonance peak in Fig. 5.9a (as evident by the RAO). Next, near the second resonance frequency, we can see three types of behavior 91 0.40.60.811.21.4/2 = Frequency (Hz)0.20.40.60.811.2RAOIPVA-PTO 5.5IPVA-PTO 3.5IPVA-PTO 3.0Linear0.40.60.811.21.4/2 = Frequency (Hz)0.511.522.53Normalized capture widthIPVA-PTO 5.5IPVA-PTO 3.5IPVA-PTO 3.0LinearRotationQP to NPPrimary a RAO. Readers are referred to the web version of this article for clarity b Normalized capture width Figure 5.10 Comparison between RAO and normalized capture width of the IPVA-PTO and the linear system for different wave heights for 𝜉 = 0.01, 𝜉𝑒,𝑙 = 0.07 and 𝜉𝑒 = 0.02, all other parameters from caption of Fig. 5.3. The corresponding bifurcation boundary can be found in Fig. 5.11d depending on the wave height. For the wave height of 3.0 cm, which is below the bifurcation boundary near the second resonance frequency, the pendulum’s response is dominated by primary harmonics, which means there is no energy transfer between the primary system and the pendulum, and the RAO closely follows the linear system for the frequency range 0.7 Hz to 1.5 Hz. If we move to 3.5 cm, we see the primary system and the pendulum giving rise to a quasi-periodic response (which eventually devolves into a non-periodic response for this wave height) due to the secondary Hopf bifurcation. In this case, the RAO of the system also rises around 1.1 Hz, and the IPVA-PTO system is worse in terms of RAO than the linear system in the neighborhood of 1.1 Hz. Furthermore, when we increase the wave height to 5.5 cm, we start seeing pendulum rotation (as expected from the markers in Fig. 5.3) for a large range of frequencies (0.5 Hz to 1.35 Hz). This is in tandem with the increase in the RAO value in the frequency range of rotation of the pendulum, and the performance is worse than the linear system for this frequency range. However, it should be noted that the maximum RAO of the linear system for the whole frequency range is worse than the IPVA-PTO system. In conclusion, for this set of wave heights, we can see the response of the pendulum going from primary to non-periodic and eventually to rotation, with rotation showing the maximum wave energy conversion. 92 0.40.60.811.21.4/2 = Frequency (Hz)0.20.40.60.811.2RAOIPVA-PTO 4.5IPVA-PTO 3.4IPVA-PTO 1.5Linear0.20.40.60.811.21.41.6/2 = Frequency (Hz)0123456Normalized capture widthIPVA-PTO 4.5IPVA-PTO 3.4IPVA-PTO 1.5LinearRotationRotationPrimary Next, we change the electrical damping coefficient to 𝜉𝑒 = 0.02 for the IPVA-PTO system while keeping the other values the same as Fig. 5.9 for both the IPVA-PTO and the linear system. Fig. 5.10 shows the RAO and normalized capture width for both the IPVA-PTO and the linear system. Similar to the case of Fig. 5.9, we can observe that after a certain wave height (which corresponds to crossing the bifurcation boundary), the IPVA-PTO system outperforms the linear benchmark in terms of the normalized capture width. Furthermore, for this set of parameters, the authors could not find a quasi-periodic or non-periodic solution; the pendulum went directly into rotation. Further analysis will be required to determine the attractor for the rotational solution, which is outside the scope of this study. Having established the effectiveness of the IPVA-PTO system compared to the linear benchmark near both resonance frequencies, we now study the effects of various system parameters on the bifurcation boundary. 5.6 Parametric studies To analyze the effects of various system parameters on the bifurcation boundary (i.e., period- doubling bifurcation and secondary Hopf bifurcation boundary), we vary the following parameters: 𝜇𝑟, 𝜂, 𝜉, and 𝜉𝑒, while keeping the other parameters fixed at their values mentioned in Fig. 5.3. First, we see the effect of 𝜇𝑟 on the bifurcation boundary. Recall that 𝜇𝑟 is the mass amplification factor in the system given by 𝜇𝑟 = 𝑚𝑅2 𝑝 𝑀 𝑅2 . As the 𝜇𝑟 value increases, we see the wave height required to cross the bifurcation boundary decrease as evident from Fig. 5.11a. Therefore, to control the wave height required to cross the boundary, one can readily change the value of 𝑅 𝑝 𝑅 , which is the ratio of the distance of the pendulum pivot point with respect to the carrier over the effective radius of the ball-screw system. Next, the effect of 𝜂, defined by 𝑟 𝑅 𝑝 , where 𝑟 is the pendulum’s length, on the bifurcation boundary is analyzed. As observed from Fig. 5.11b, for larger values of 𝜂, the wave height required to cross the boundary for a given frequency decreases. Fig. 5.11c and Fig. 5.11d show the effect of mechanical and electrical damping respectively on the bifurcation boundary. It can be observed that the wave height required to cross the bifurcation boundary increases with an increase in either the mechanical or electrical damping value. From these parametric variation figures, it can be observed that the two local minima of the 93 bifurcation boundary can be positioned higher or lower relative to each other by changing the system parameters. If we are to install the IPVA-PTO system for energy conversion, the values of 𝜇𝑟, 𝜂, and 𝜉 cannot be changed once the system is manufactured. However, we can easily control the electrical damping by manipulating the system’s electrical components. Therefore, if the aim is to convert energy and suppress the response of the spar simultaneously, one has to minimize the RAO near the first resonance frequency (RAO at the first resonance frequency is typically larger than at the second resonance) and maximize the capture width near the second resonance frequency (as evident from Fig. 5.9 and Fig. 5.10). The reason is as follows: the power spectral density of wave excitation, although narrow band, is essentially random, for example, JONSWAP [90]. Typically, spars are designed such that their natural frequency (around 0.33 Hz to 0.5 Hz) [24, 25] is around half to one-sixth of typical incident wave frequency (around 0.1 to 0.2 Hz) [26]. Although small, there is a nonzero probability that there is an excitation near the first resonance frequency of the system, which can compromise the stability of the spar. On the other hand, the second resonance frequency of the system can be made close to the peak frequency of the wave’s power spectral density by design; thus, it is better to convert energy near the second resonant peak, while maintaining the stability of the spar near the first resonance frequency. From the observations in Sec. 5.5, it can be said that the maximum RAO will be reduced when the wave height is above the bifurcation boundary for the first resonance frequency, and the maximum capture width is maximized. Looking at Fig. 5.11d, it can be observed that when the electrical damping is varied, the wave height corresponding to the first resonance frequency of the bifurcation boundary changes significantly, whereas the height of the second resonance frequency does not change. Let us assume that we have an incoming wave height of less than 3 cm with frequency varying over time. It may be possible to generate larger energy by keeping the electrical damping 𝜉𝑒 = 0.02, but to save the spar from damage (by keeping the RAO as low as possible), the electrical damping value of 0.01 is preferred so that the spar is stable. Therefore, one can control the electrical damping value to keep the spar safe while generating significant energy compared to a linear benchmark. 94 a Bifurcation boundary for different values of 𝜇𝑟 b Bifurcation boundary for different values of 𝜂 c Bifurcation boundary for different values of 𝜉 d Bifurcation boundary for different values of 𝜉𝑒 Figure 5.11 Effect of various system parameters on the bifurcation boundary, all the values are from Fig. 5.3 except for the parameter value specified in the figure. Readers are referred to the web version of this article for clarity 5.7 Optimization of electrical damping value As can be seen from Fig. 5.11d, the electrical damping in the system can affect the wave heights corresponding to both the resonance frequencies in the system. For example, it is clear that the wave height corresponding to the first resonance frequency is higher than the wave height at the second resonance frequency at 𝜉𝑒 = 0.025, whereas, for 𝜉𝑒 = 0.015, the wave height at the first resonance frequency is lower than the second. It can also be observed that the wave height at the second resonance frequency is not as sensitive compared to the first resonance frequency. Therefore, it can be assumed that there is a scope to optimize the electrical damping for best energy 95 0.511.5/2 = Frequency (Hz)0.020.030.040.050.060.070.080.090.10.11Wave height (m)r=0.3r=0.35r=0.40.511.5/2 = Frequency (Hz)0.020.030.040.050.060.070.080.090.10.110.12Wave height (m)=0.35=0.4=0.450.511.5/2 = Frequency (Hz)0.020.030.040.050.060.070.080.090.10.11Wave height (m)=0.01=0.03=0.050.511.5Frequency0.020.040.060.080.10.120.14Wave height (m)e=0.01e=0.015e=0.02 conversion and response suppression. To that end, after fixing a wave height, different electrical damping values are taken and the maximum RAO (maxRAO) and the 𝐻2-norm of non-dimensional capture width (H2CW), obtained by normalizing the maximum value of capture width to 1, over the frequency range [0.25, 1.5] Hz are obtained corresponding to each 𝜉𝑒 value. The reason for choosing maxRAO and H2CW as measures of efficacy of wave energy converter is due to the fact the ocean waves are essentially random. Therefore, one would like to avoid the worst RAO response (hence the maxRAO), and maximize the energy capture throughout the frequency range around both resonances (H2CW). Next, a Pareto front is obtained in the space of maxRAO and H2CW by varying the electrical damping. A point in the maxRAO-H2CW space is said to dominate the other when its H2CW is higher and maxRAO is lower than the other solution. For this study, two wave heights are considered: 5 cm and 4 cm. The parameters from Fig. 5.3 are used, except the electrical damping which is varied from the range 𝜉𝑒 ∈ [0.015, 0.035] for the wave height 5 cm and 𝜉𝑒 ∈ [0.005, 0.025] for the wave height 4 cm, with a discretization of 0.0002 used for simulations. The electrical damping range was chosen such that the system response becomes strictly primary harmonic at the maximum value in the range. Fig. 5.12a shows the Pareto front for the wave height 5 cm. From Fig. 5.13, the variation of maxRAO and H2CW can be observed with respect to the electrical damping. The maximum value of H2CW occurs at 𝜉𝑒 = 0.0224, where the maxRAO value is around 0.9621. Now from Fig. 5.3, it is known that the non-periodic response typically occurs before rotation response of the pendulum. As electrical damping increases, the system response should go from rotation to non-periodic as evident from the increase in the wave height corresponding to bifurcation boundary in Fig. 5.3. In general, as the electrical damping increases, the response of the pendulum changes like this: Intermittent Rotation → Rotation → Non-periodic → Secondary harmonic/Quasi-periodic (depending on the frequency value) → Primary harmonic. Now, from Fig. 5.14b at 𝜉𝑒 = 0.0224, the H2CW value is 1 and it drops quickly to around 0.62 at 𝜉𝑒 = 0.0262. The hypothesis is that when this drop happens, the solution goes from rotation to non-periodic (as electrical damping is increasing), and the energy conversion capability is reduced greatly. This can be verified by the frequency response of RAO and normalized capture width 96 a Pareto front for wave height = 5 cm b Pareto front for wave height = 4 cm Figure 5.12 ‘◦’ represents the Pareto frontal (non-dominated) solutions, ‘∗’ represents all the solutions for 𝜉𝑒 = 0.0224 and 𝜉𝑒 = 0.0262 in Fig. 5.13. As can be seen for the case of 𝜉𝑒 = 0.0224, the normalized capture width is significantly higher and the response is rotational (as verified from the time series, not shown here) as opposed to lower normalized capture width and non-periodic response of the pendulum in the case of 𝜉𝑒 = 0.0262. The RAO near the first peak is worse for 𝜉𝑒 = 0.0224 compared to 𝜉𝑒 = 0.0262, however the global RAO is worse for the case of 𝜉𝑒 = 0.0262. Now, let’s look at the wave height 4 cm. The Pareto front is shown in Fig. 5.14. Similar to the case of wave height 5 cm, Fig. 5.14 shows the variation of maxRAO and H2CW with respect to the electrical damping. Again, we suspect that as the solution goes from rotation to non-periodic with increasing electrical damping, and the energy conversion capability is reduced greatly. From the time series of the solution, it can be verified that the peak of H2CW (for electrical damping 𝜉𝑒 = 0.0162) has rotational pendulum response. Now, let us again take two electrical damping values: 𝜉𝑒 = 0.0162 and 𝜉𝑒 = 0.0202, and simulate the frequency response of RAO and normalized capture width. As can be seen from Fig. 5.13, from the case of 𝜉𝑒 = 0.0162, the normalized capture width is significantly higher when the response is rotational and reduces significantly when the pendulum’s response becomes non-periodic for 𝜉𝑒 = 0.0202. However, in this case, the RAO is globally worse for 𝜉𝑒 = 0.0162 compared to 𝜉𝑒 = 0.0202. 97 0.550.60.650.70.750.80.850.90.951H2CW0.511.522.53maxRAO0.40.50.60.70.80.91H2CW0.511.522.533.5maxRAO a Maximum RAO as a function of electrical damping b 𝐻2 norm of energy as a function of electrical damping Figure 5.13 maxRAO and H2CW as a function of the electrical damping value for wave height = 5 cm a Maximum RAO as a function of electrical damping b 𝐻2 norm of energy as a function of electrical damping Figure 5.14 maxRAO and H2CW as a function of the electrical damping value for wave height = 4 cm 98 0.0150.020.0250.030.035Electrical Damping11.21.41.61.822.22.42.62.8Max RAO0.0150.020.0250.030.035Electrical Damping0.60.650.70.750.80.850.90.951H2CW0.0050.010.0150.020.025Electrical Damping0.511.522.533.5Max RAO0.0050.010.0150.020.025Electrical Damping0.40.50.60.70.80.91H2CW a RAO b Normalized capture width Figure 5.15 RAO and normalized capture width for two different values of 𝜉𝑒 for wave height 5 cm a RAO b Normalized capture width Figure 5.16 RAO and normalized capture width for two different values of 𝜉𝑒 for wave height 4 cm 5.8 Conclusion This study analyzes the incorporation of the IPVA system [70] with a power take-off into a heaving spar-floater system to study the integration for hydrodynamic response suppression and wave energy conversion. The hydrodynamic response and wave energy conversion are investigated using numerical frequency response simulations. A harmonic balance method determines the wave height for which the primary harmonic solution bifurcates. Two types of bifurcation of the primary harmonic solution are observed: period-doubling bifurcation near the first resonance and secondary Hopf bifurcation near the second resonance in the parameter plane of wave height and 99 0.20.40.60.811.21.41.6Frequency (Hz)00.511.522.5RAOe=0.0224e=0.02620.20.40.60.811.21.41.6Frequency (Hz)0123456Normalized capture widthe=0.0224e=0.02620.20.40.60.811.21.41.6Frequency (Hz)00.20.40.60.811.2RAOe=0.0162e=0.02020.20.40.60.811.21.41.6Frequency (Hz)00.511.522.533.544.55Energye=0.0162e=0.0202 wave frequency. According to the boundary, one can determine a combination of the wave height and frequency such that the primary solution becomes unstable and either 1:2 internal resonance or quasi-periodic response occurs to the IPVA-PTO system. Near the first resonance frequency, above the wave height required for bifurcation, a rich set of pendulum responses are observed, like secondary harmonics, non-periodic solutions, and rotation. Furthermore, it is observed that the bifurcation of primary harmonic near the first resonance frequency is associated with a nonlinear energy transfer phenomenon similar to that observed in [70, 83] when the electrical damping is sufficiently small. It is also shown that because of this energy transfer phenomenon, the IPVA- PTO system achieves a lower maximum RAO and a higher energy conversion than the linear benchmark. For higher values of electrical damping, the energy transfer phenomenon is not very clear, but the IPVA-PTO system still outperforms the linear system in terms of both maximum RAO and normalized capture width. Near the second resonance frequency, the pendulum’s response can involve quasi-periodic solutions, non-periodic oscillations, and rotation. It is shown that the rotation of the pendulum is excellent for wave energy conversion compared to the linear benchmark. Even though the linear system has a better maximum RAO near the second resonance frequency, the IPVA-PTO system has a globally lower maximum RAO (the overall maximum RAO typically occurs near the first resonance frequency). The effect of system parameters is studied on the bifurcation boundary of the system, with a particular focus on electrical damping since it can be actively controlled in practice. By controlling the electrical damping, the wave height required for crossing the bifurcation boundary at first resonance frequency can be changed significantly, which can be exploited to stabilize the spar. Finally, studies on electrical damping show that the normalized capture width of the system is maximum when the pendulum’s response is rotational and quickly worsens if the electrical damping is increased to the point that pendulum’s response becomes non-periodic. 100 CHAPTER 6 CONCLUSION AND PROPOSAL FOR THE FUTURE WORK 6.1 Concluding remarks In this work, a new system described as an inerter pendulum vibration absorber was proposed (IPVA). The IPVA system works on the principle of inertial nonlinearity to broaden the frequency bandwidth of the response and thereby allowing vibration suppression and energy harvesting over a larger range of frequencies. The IPVA system was analyzed in two forms - the first form is the dry IPVA system integrated with a single-degree-of-freedom system to observe the vibration suppression and energy harvesting, and the second form where the IPVA system was integrated to the spar for the application of ocean wave energy harvesting. In both cases, the pendulum’s parameters can be chosen such that 1:2 internal resonance happens for a range of frequencies. It has been shown that a pitchfork bifurcation and a subsequent period-doubling bifurcation are necessary and sufficient conditions for this 1:2 internal resonance to occur. When this internal resonance occurs, it leads to energy transfer between the primary system and the pendulum vibration absorber. This means the vibration of the primary system is suppressed, whereas the energy can be harvested from the pendulum. A saturation phenomenon similar to that observed in autoparametric vibration absorbers and other nonlinear vibration absorbers is observed in the IPVA system; that is, the response of the linear oscillator saturates despite the increase in the force magnitude. Meanwhile, the increased energy due to the increase in the force magnitude seems to transfer to the pendulum, resulting in an increased pendulum response. Furthermore, the dry system is compared to the autoparametric vibration absorber. It is shown to outperform the autoparametric vibration absorber in terms of vibration absorption and energy harvesting capabilities. Similarly, the wet system (consisting of IPVA integrated with spar) is shown to outperform the linear benchmark system in terms of energy harvesting and vibration suppression. This analysis is verified by performing experiments on the single-degree-of-freedom IPVA and “dry” IPVA-PTO system, where the primary harmonic solution is shown to bifurcate via a period-doubling bifurcation. The secondary resonance which corresponds to the crossing of 1:2 101 internal resonance boundary is experimentally shown for both the experimental systems. For the case of the IPVA system, that above the period-doubling bifurcation boundary (1:2 internal resonance boundary), the motion of the primary system was suppressed, and its energy transferred to the pendulum vibration absorber for two different acceleration values. In the case of the “dry” IPVA-PTO system, a rich set of pendulum responses (primary harmonic, secondary harmonic, non-periodic, rotation) is observed. Two different load resistance values are used to compare the IPVA-PTO and the linear system, and it is shown that although the non-periodic response of the IPVA-PTO system is only suitable for response suppression, the presence of secondary harmonics can outperform the linear benchmark both in response suppression and energy conversion. Finally, the IPVA and the IPVA-PTO system were integrated with a heaving spar-floater system to study the integration for hydrodynamic response suppression and wave energy conversion. The hydrodynamic response and wave energy conversion are investigated using numerical frequency response simulations. A harmonic balance method determines the wave height for which the primary harmonic solution bifurcates. For the case of the IPVA system, a period-doubling bifurcation boundary can be determined, which is associated with a nonlinear energy transfer phenomenon, resulting in a lower maximum RAO, and a higher energy transfer potential compared to the linear benchmark, when the relative angular motion between the ball-screw and the pendulum is used as a measure of energy conversion potential. For the case of the IPVA-PTO system, two types of bifurcation of the primary harmonic solution are observed: period-doubling bifurcation near the first resonance and secondary Hopf bifurcation near the second resonance in the parameter plane of wave height and wave frequency. Near the first resonance frequency, above the wave height required for bifurcation, a rich set of pendulum responses are observed, like secondary harmonics, non-periodic solutions, and rotation. Furthermore, it is observed that the bifurcation of primary harmonic near the first resonance frequency is associated with a nonlinear energy transfer phenomenon similar to that observed in [70, 83] when the electrical damping is sufficiently small. It is also shown that because of this energy transfer phenomenon, the IPVA-PTO system achieves a lower maximum RAO and a higher energy conversion than the linear benchmark. For higher 102 values of electrical damping, the energy transfer phenomenon is not very clear, but the IPVA-PTO system still outperforms the linear system in terms of both maximum RAO and normalized capture width. Near the second resonance frequency, the pendulum’s response can involve quasi-periodic solutions, non-periodic oscillations, and rotation. It is shown that the rotation of the pendulum is excellent for wave energy conversion compared to the linear benchmark. Even though the linear system has a better maximum RAO near the second resonance frequency, the IPVA-PTO system has a globally lower maximum RAO (the overall maximum RAO typically occurs near the first resonance frequency). The effect of system parameters is studied on the bifurcation boundary of both the IPVA and IPVA-PTO system, with a particular focus on pendulum damping for the IPVA system and electrical damping for the IPVA-PTO system since it can be actively controlled in practice. In the case of the IPVA-PTO system, by controlling the electrical damping, the wave height required for crossing the bifurcation boundary at first resonance frequency can be changed significantly, which can be exploited to stabilize the spar. Finally, for the IPVA-PTO system, studies on electrical damping show that the normalized capture width of the system is maximum when the pendulum’s response is rotational and quickly worsens if the electrical damping is increased to the point that the pendulum’s response becomes non-periodic. 6.2 Suggested future direction As has been concluded in this work, the spar-floater IPVA-PTO system shows the highest capture width when the response of the pendulum is rotational. The rotations of the pendulum are primary in nature, in the sense that if we assume 𝜙 to be the pendulum angle of rotation, 𝜙 = 𝑆𝜔𝜏 + ¯𝜙 (6.1) where 𝑆 can be +1 or −1 depending on the rotation direction of the pendulum, with ¯𝜙 = 𝜙0 + 𝑝ℎ𝑖𝑐 cos (𝜔𝜏) + 𝜙𝑠 sin (𝜔𝜏) in this case, where 𝜔 is the normalized excitation angular frequency, 𝜏 is normalized time, see (5.6), and 𝜙𝑐 and 𝜙𝑠 are constants. We call this solution the primary rotation. In the plane of frequency and wave height, the solution below the rotational response is non-periodic in nature and has been shown to be unfavorable for energy conversion in this 103 work. Therefore, future work can involve the prediction of rotational solutions of the pendulum. To that end, one can substitute (6.1) in (5.3) and obtain an equation in terms of new variables 𝜃, ¯𝜙, and 𝜓. Then the bifurcation boundary where the primary harmonic solution of the system becomes unstable and turns into a non-periodic solution can be obtained using the harmonic balance approach coupled with the modified AFT method. Therefore, in practice, with the knowledge of the bifurcation boundary above which primary rotations of the spar-floater IPVA-PTO system are stable, the electrical damping can be controlled such that the pendulum’s response is always rotating. Similarly, as the wave height increases, it has been observed that the primary rotations bifurcate either into secondary harmonic rotations or quasi-periodic rotations (defined as a combination of a secondary harmonic component with rotation or a quasi-periodic component with rotation, respectively). If the wave height is increased further, it can be seen that the solutions further bifurcate into non-periodic solutions characterized by intermittent rotations. One would ideally like to avoid non-periodic intermittent rotations for electricity generation. Therefore, a bifurcation boundary that predicts the bifurcation of primary rotation into either secondary harmonic rotation or quasi-periodic rotation can be useful to ensure the electrical damping is not too low that the solution is not primary rotational anymore. 6.2.1 Preliminary analysis For the first study to analyze the primary rotations of the pendulum, we do not consider the hydrodynamic effects on the system as a simplification. Furthermore, we make a small angle approximation for the harmonic part of the pendulum’s motion, that is, ¯𝜙 = 𝜙0 + 𝑁 ∑︁ 𝑘=1 𝜙𝑘𝑐 cos (𝑘𝜔𝜏) + 𝜙𝑘 𝑠 sin (𝑘𝜔𝜏) (6.2) with 𝜙ℎ = ¯𝜙 − 𝜙0, sin (𝜙ℎ) ≈ 𝜙ℎ − 𝜙3 ℎ/6, cos (𝜙ℎ) ≈ 1 − 𝜙2 ℎ/2 104 (6.3) a 𝜃 b 𝜓 Figure 6.1 Comparison 𝜃, 𝜓 and 𝜙 responses between harmonic balance method and direct numerical simulation for wave height 3.7 cm and excitation frequency of 1.12 Hz c 𝜙 Similar to 𝜙, 𝜃 and 𝜓 can be written as 𝜃 = 𝜃0 + 𝜓 = 𝜓0 + 𝑁 ∑︁ 𝑘=1 𝑁 ∑︁ 𝑘=1 𝜃 𝑘𝑐 cos (𝑘𝜔𝜏) + 𝜃 𝑘 𝑠 sin (𝑘𝜔𝜏) , 𝜓𝑘𝑐 cos (𝑘𝜔𝜏) + 𝜓𝑘 𝑠 sin (𝑘𝜔𝜏) . For the preliminary study, 𝑁 = 2 is chosen to apply the harmonic balance analysis. The results of the harmonic balance method (HBM) are compared with direct numerical simulation (DNS) using Matlab’s ode45 and are shown in Fig. 6.1. As can be observed, the harmonic balance analysis predicts the solution with excellent accuracy. Therefore, we can move forward with the bifurcation analysis to identify the stability zones for the primary rotation of the pendulum. Using harmonic balance with Floquet theory, two stability boundaries are obtained in the plane of wave height and wave frequency, which bound the primary rotation responses of the pendulum. In Fig. 6.2a, the region below the top boundary and above the bottom boundary corresponds to the primary rotations of the pendulum. To verify this, we choose a wave height = 3.7 cm and plot the frequency response of the pendulum’s motion. In Fig. 6.2b, the response of the pendulum is primary rotations from a frequency range of 1.03 Hz to 1.22 Hz (as verified from the time series, not shown here), whereas the stability boundaries predict a primary rotation between the range [1.2, 1.28] Hz. The error in the prediction can be attributed to the small angle approximation of 𝜙ℎ shown in (6.3). Therefore, more analysis needs to be performed to accurately predict the primary rotations in the system. 105 2174217621782180218221842186-0.6-0.4-0.200.20.40.6DNSHBM2174217621782180218221842186-1-0.500.51DNSHBM2174217621782180218221842186477547804785479047954800DNSHBM a Stability boundaries for primary rotations b Response of the pendulum for wave height = 3.7 cm Figure 6.2 Stability boundaries for primary rotations and the response of the pendulum for wave height = 3.7 cm After obtaining an accurate prediction of primary rotations in the system with no added mass and radiation damping, the next logical step would be to integrate the added mass and radiation damping in the system and predict the occurrence of primary rotations. These predictions can be used to tune the electrical damping of the system such that the pendulum stays in the primary rotation region, given a wave frequency and a wave height, resulting in maximum energy conversion of the system. 106 1.11.21.31.4Frequency (Hz)0.010.020.030.040.050.060.070.08Wave height (m)TopBottom11.051.11.151.21.251.31.35Frequency (Hz)11.522.533.5 BIBLIOGRAPHY [1] Johannes Falnes. A review of wave-energy extraction. Marine structures, 20(4):185–201, 2007. [2] Stephen Barstow, Gunnar Mørk, Lasse Lønseth, and Jan Petter Mathisen. Worldwaves wave energy resource assessments from the deep ocean to the coast. Journal of Energy and Power Engineering, 5(8):730–742, 2011. [3] Andrew M Cornett. A global wave energy resource assessment. In The Eighteenth international offshore and polar engineering conference. OnePetro, 2008. [4] Kester Gunn and Clym Stock-Williams. Quantifying the global wave power resource. Renewable Energy, 44:296–304, 2012. [5] Roger Bedard, George Hagerman, Mirko Previsic, Omar Siddiqui, Robert Thresher, and Bonnie Ram. Final summary report, project definition study, offshore wave power feasibility demonstration project. Electric Power Research Institute Inc, 2005. [6] HP Nguyen, CM Wang, ZY Tay, and VH Luong. Wave energy converter and large floating platform integration: A review. Ocean Engineering, 213:107768, 2020. [7] Shangyan Zou, Ossama Abdelkhalik, Rush Robinett, Giorgio Bacelli, and David Wilson. Optimal control of wave energy converters. Renewable energy, 103:217–225, 2017. [8] Elena Baca, Ritu Treisa Philip, David Greene, and Hoyt Battey. Expert elicitation for wave energy lcoe futures. Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States), 2022. [9] US Energy Information Association. Levelized costs of new generation resources in the annual energy outlook 2022. US Department of Energy, March, 2022. [10] Carlos Pérez-Collazo, D Greaves, and G Iglesias. A review of combined wave and offshore wind energy. Renewable and Sustainable Energy Reviews, 42:141–153, 2015. [11] Made Jaya Muliawan, Madjid Karimirad, and Torgeir Moan. Dynamic response and power performance of a combined spar-type floating wind turbine and coaxial floating wave energy converter. Renewable Energy, 50:47–57, 2013. [12] Madjid Karimirad and Kourosh Koushan. Windwec: Combining wind and wave energy In 2016 IEEE International Conference on Renewable inspired by hywind and wavestar. Energy Research and Applications (ICRERA), pages 96–101. IEEE, 2016. [13] Zhengshun Cheng, Ting Rui Wen, Muk Chen Ong, and Kai Wang. Power performance and dynamic responses of a combined floating vertical axis wind turbine and wave energy converter concept. Energy, 171:190–204, 2019. [14] Constantine Michailides. Hydrodynamic response and produced power of a combined structure consisting of a spar and heaving type wave energy converters. Energies, 14(1):225, 2021. 107 [15] Bureau of Safety and Environmental Enforcement. “platform structures online query,”, FEB, 2021. Accessed: February 1, 2021. [16] Aurélien Babarit, Jorgen Hals, Made Jaya Muliawan, Adi Kurniawan, Torgeir Moan, and Jorgen Krokstad. Numerical benchmarking study of a selection of wave energy converters. Renewable energy, 41:44–63, 2012. [17] C Michailides. Hydrodynamic response and produced power of a combined structure consisting of a spar and heaving type wave energy converters. energies 2021, 14, 225, 2021. [18] Made Jaya Muliawan, Madjid Karimirad, Torgeir Moan, and Zhen Gao. Stc (spar-torus combination): a combined spar-type floating wind turbine and large point absorber floating wave energy converter—promising and challenging. In International Conference on Offshore Mechanics and Arctic Engineering, volume 44946, pages 667–676. American Society of Mechanical Engineers, 2012. [19] Ling Wan, Zhen Gao, and Torgeir Moan. Experimental and numerical study of hydrodynamic responses of a combined wind and wave energy converter concept in survival modes. Coastal Engineering, 104:151–169, 2015. [20] Aurélien Babarit. A database of capture width ratio of wave energy converters. Renewable Energy, 80:610–628, 2015. [21] Bureau of safety and environmental enforcement. (feb, 2021) “platform structures online query,”. Accessed on Feb, 2021. [22] Hideyuki Suzuki and Akira Sato. Load on turbine blade induced by motion of floating platform and design requirement for the platform. In International Conference on Offshore Mechanics and Arctic Engineering, volume 42711, pages 519–525, 2007. [23] Minnan Yue, Qingsong Liu, Chun Li, Qinwei Ding, Shanshan Cheng, and Haitian Zhu. Effects of heave plate on dynamic response of floating wind turbine spar platform under the coupling effect of wind and wave. Ocean Engineering, 201:107103, 2020. [24] BJ Koo, MH Kim, and RE Randall. Mathieu instability of a spar platform with mooring and risers. Ocean engineering, 31(17-18):2175–2208, 2004. [25] A Subbulakshmi and R Sundaravadivelu. Heave damping of spar platform for offshore wind turbine with heave plate. Ocean Engineering, 121:24–36, 2016. [26] Changwei Liang. On the Dynamics and Design of a Wave Energy Converter with Mechanical Motion Rectifier. PhD thesis, State University of New York at Stony Brook, 2016. [27] Fantai Meng, Nataliia Sergiienko, Boyin Ding, Binzhen Zhou, Leandro Souza Pinheiro Da Silva, Benjamin Cazzolato, and Ye Li. Co-located offshore wind–wave energy systems: Can motion suppression and reliable power generation be achieved simultaneously? Applied Energy, 331:120373, 2023. 108 [28] Zili Zhang and Christian Høeg. Vibration control of floating offshore wind turbines using liquid column dampers. In Journal of physics: conference series, volume 1037, page 032002. IOP Publishing, 2018. [29] Xiantao Zhang, Haicheng Zhang, Xiao Zhou, and Ze Sun. Recent advances in wave energy converters based on nonlinear stiffness mechanisms. Applied Mathematics and Mechanics, 43(7):1081–1108, 2022. [30] Benjamin W Schubert, Nataliia Y Sergiienko, Benjamin S Cazzolato, William SP Robertson, and Mergen H Ghayesh. The true potential of nonlinear stiffness for point absorbing wave energy converters. Ocean Engineering, 245:110342, 2022. [31] Zhijia Wu, Carlos Levi, and Segen F Estefen. Wave energy harvesting using nonlinear stiffness system. Applied Ocean Research, 74:102–116, 2018. [32] Elie Al Shami, Xu Wang, and Xueyu Ji. A study of the effects of increasing the degrees of freedom of a point-absorber wave energy converter on its harvesting performance. Mechanical Systems and Signal Processing, 133:106281, 2019. [33] Mohammed F Daqaq, Ravindra Masana, Alper Erturk, and D Dane Quinn. On the role of nonlinearities in vibratory energy harvesting: a critical review and discussion. Applied Mechanics Reviews, 66(4):040801, 2014. [34] Lei Zuo, Gopinath Reddy Penamalli, John Wang, Rui He Zheng, Xiao Hui Lei, Jorge F Lam- Ki, Zhongjie Li, and Teng Lin. High-efficiency energy generator for harnessing mechanical vibration power, July 19 2016. US Patent 9,394,876. [35] Changwei Liang, You Wu, and Lei Zuo. Broadband pendulum energy harvester. Smart Materials and Structures, 25(9):095042, 2016. [36] Wei-Che Tai, Mingyi Liu, Yue Yuan, and Lei Zuo. On improvement of the frequency bandwidth of nonlinear vibration energy harvesters using a mechanical motion rectifier. Journal of Vibration and Acoustics, 140(5):051008, 2018. [37] Mingyi Liu, Wei-Che Tai, and Lei Zuo. Toward broadband vibration energy harvesting via mechanical motion-rectification induced inertia nonlinearity. Smart Materials and Structures, 27(7):075022, 2018. [38] Malcolm C Smith. The inerter: A retrospective. Annual Review of Control, Robotics, and Autonomous Systems, 3:361–391, 2020. [39] Kohju Ikago, Kenji Saito, and Norio Inoue. Seismic control of single-degree-of-freedom structure using tuned viscous mass damper. Earthquake Engineering & Structural Dynamics, 41(3):453–474, 2012. [40] IF Lazar, SA Neild, and DJ Wagg. Using an inerter-based device for structural vibration suppression. Earthquake Engineering & Structural Dynamics, 43(8):1129–1147, 2014. [41] IF Lazar, SA Neild, and DJ Wagg. Vibration suppression of cables using tuned inerter dampers. Engineering Structures, 122:62–71, 2016. 109 [42] Feng Qian, Yifan Luo, Hongxin Sun, Wei Che Tai, and Lei Zuo. Optimal tuned inerter dampers for performance enhancement of vibration isolation. Engineering Structures, 198:109464, 2019. [43] Laurentiu Marian and Agathoklis Giaralis. Optimal design of a novel tuned mass-damper– inerter (tmdi) passive vibration control configuration for stochastically support-excited structural systems. Probabilistic Engineering Mechanics, 38:156–164, 2014. [44] Dario De Domenico and Giuseppe Ricciardi. An enhanced base isolation system equipped with optimal tuned mass damper inerter (tmdi). Earthquake Engineering & Structural Dynamics, 47(5):1169–1192, 2018. [45] Eshagh Farzaneh Joubaneh and Oumar Rafiou Barry. On the improvement of vibration mitigation and energy harvesting using electromagnetic vibration absorber-inerter: Exact h2 optimization. Journal of Vibration and Acoustics, 141(6), 2019. [46] Wei-Che Tai. Optimum design of a new tuned inerter-torsional-mass-damper passive vibration control for stochastically motion-excited structures. Journal of Vibration and Acoustics, 142(1), 2020. [47] Takehiko Asai and Keita Sugiura. Numerical evaluation of a two-body point absorber wave energy converter with a tuned inerter. Renewable energy, 171:217–226, 2021. [48] Feng Qian and Lei Zuo. Tuned nonlinear spring-inerter-damper vibration absorber for beam vibration reduction based on the exact nonlinear dynamics model. Journal of Sound and Vibration, 509:116246, 2021. [49] Paul Kakou and Oumar Barry. Simultaneous vibration reduction and energy harvesting of a nonlinear oscillator using a nonlinear electromagnetic vibration absorber-inerter. Mechanical Systems and Signal Processing, 156:107607, 2021. [50] Dynamic analysis and performance evaluation of nonlinear inerter-based vibration isolators. Nonlinear Dyn, 99:1823–1839, 2020. [51] H Hatwal, AK Mallik, and A Ghosh. Forced nonlinear oscillations of an autoparametric system—part 1: periodic responses. 1983. [52] A Vyas and AK Bajaj. Dynamics of autoparametric vibration absorbers using multiple pendulums. Journal of Sound and Vibration, 246(1):115–135, 2001. [53] AK Bajaj, SI Chang, and JM Johnson. Amplitude modulated dynamics of a resonantly excited autoparametric two degree-of-freedom system. Nonlinear dynamics, 5(4):433–457, 1994. [54] Ali H Nayfeh and Dean T Mook. Nonlinear oscillations. John Wiley & Sons, 2008. [55] Krzysztof Kecik. Assessment of energy harvesting and vibration mitigation of a pendulum dynamic absorber. Mechanical Systems and Signal Processing, 106:198–209, 2018. [56] H Hatwal, AK Mallik, and A Ghosh. Non-linear vibrations of a harmonically excited autoparametric system. Journal of Sound and Vibration, 81(2):153–164, 1982. 110 [57] Y Song, H Sato, Y Iwata, and T Komatsuzaki. The response of a dynamic vibration absorber system with a parametrically excited pendulum. Journal of Sound and Vibration, 259(4):747– 759, 2003. [58] Jerzy Warminski and Krzysztof Kecik. Instabilities in the main parametric resonance area of a mechanical system with a pendulum. Journal of Sound and Vibration, 322(3):612–628, 2009. Special issue from the Second International Conference on Nonlinear Dynamics, Kharkov, Ukraine, 25-28 September 2007, held in honour of the 150th anniversary of Alexander Mikhailovich Lyapunov. [59] M Sharif-Bakhtiar and SW Shaw. Effects of nonlinearities and damping on the dynamic response of a centrifugal pendulum vibration absorber. 1992. [60] Thibaut Detroux, Ludovic Renson, Luc Masset, and Gaëtan Kerschen. The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems. Computer Methods in Applied Mechanics and Engineering, 296:18–38, 2015. [61] David Edward Newland. Nonlinear aspects of the performance of centrifugal pendulum vibration absorbers. 1964. [62] MN Hamdan and TD Burton. On the steady state response and stability of non-linear oscillators using harmonic balance. Journal of sound and vibration, 166(2):255–266, 1993. [63] Ivana Kovacic, Richard Rand, and Si Mohamed Sah. Mathieu’s equation and its generalizations: overview of stability charts and their features. Applied Mechanics Reviews, 70(2), 2018. [64] Lihan Xie, Sébastien Baguet, Benoit Prabel, and Régis Dufour. Bifurcation tracking by harmonic balance method for performance tuning of nonlinear dynamical systems. Mechanical Systems and Signal Processing, 88:445–461, 2017. [65] Thomas L Carroll and Louis M Pecora. Nonlinear dynamics in circuits. World Scientific, 1995. [66] Emmanuel Gourdon, Nicholas A Alexander, Colin A Taylor, Claude-Henri Lamarque, and Stéphane Pernot. Nonlinear energy pumping under transient forcing with strongly nonlinear coupling: Theoretical and experimental results. Journal of sound and vibration, 300(3- 5):522–551, 2007. [67] Ting Tan, Zhimiao Yan, Yajian Zou, and Wenming Zhang. Optimal dual-functional design for a piezoelectric autoparametric vibration absorber. Mechanical Systems and Signal Processing, 123:513–532, 2019. [68] Shafic S Oueini, Ali H Nayfeh, and Jon R Pratt. A nonlinear vibration absorber for flexible structures. Nonlinear Dynamics, 15(3):259–282, 1998. [69] K Kecik and M Borowiec. An autoparametric energy harvester. The European Physical Journal Special Topics, 222(7):1597–1605, 2013. 111 [70] Aakash Gupta and Wei-Che Tai. The response of an inerter-based dynamic vibration absorber with a parametrically excited centrifugal pendulum. Journal of Vibration and Acoustics, pages 1–39, 2022. [71] Aakash Gupta and Wei-Che Tai. Broadband and enhanced energy harvesting using inerter pendulum vibration absorber. In International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, volume 83969, page V007T07A007. American Society of Mechanical Engineers, 2020. [72] Changwei Liang, Junxiao Ai, and Lei Zuo. Design, fabrication, simulation and testing of an ocean wave energy converter with mechanical motion rectifier. Ocean Engineering, 136:190– 200, 2017. [73] T Seebai and R Sundaravadivelu. Response analysis of spar platform with wind turbine. Ships and Offshore Structures, 8(1):94–101, 2013. [74] Alexis Mérigaud. A harmonic balance framework for the numerical simulation of non-linear wave energy converter models in random seas. PhD thesis, National University of Ireland, Maynooth (Ireland), 2018. [75] WE Cummins, Wwllfil Iiuhl, and Alu Uinm. The impulse response function and ship motions. 1962. [76] T Francis Ogilvie. Recent progress toward the understanding and prediction of ship motions. In 5th ONR Symp. on Naval Hydrodynamics, 1964. [77] Ronald W Yeung. Added mass and damping of a vertical cylinder in finite-depth waters. Applied Ocean Research, 3(3):119–133, 1981. [78] TM Cameron and Jerry H Griffin. An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems. 1989. [79] Thibaut Detroux, Ludovic Renson, Luc Masset, and Gaëtan Kerschen. The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems. Computer Methods in Applied Mechanics and Engineering, 296:18–38, 2015. [80] Erlend Kristiansen, Åsmund Hjulstad, and Olav Egeland. State-space representation of radiation forces in time-domain vessel models. Ocean Engineering, 32(17-18):2195–2216, 2005. [81] Aakash Gupta, Van Tuan Kiet Duong, and Wei-Che Tai. Nonlinear Energy Transfer of a Spar-Floater System Using the Inerter Pendulum Vibration Absorber. Journal of Vibration and Acoustics, 145(5):051005, 09 2023. [82] David J Ewins. Modal testing: theory, practice and application. John Wiley & Sons, 2009. [83] Aakash Gupta and Wei-Che Tai. Ocean wave energy conversion of a spar platform using a nonlinear inerter pendulum vibration absorber. In International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, volume 86311, page V010T10A013. American Society of Mechanical Engineers, 2022. 112 [84] Scott J Beatty, Matthew Hall, Bradley J Buckham, Peter Wild, and Bryce Bocking. Experimental and numerical comparisons of self-reacting point absorber wave energy converters in regular waves. Ocean Engineering, 104:370–386, 2015. [85] T Seebai and R Sundaravadivelu. Dynamic analysis of slack moored spar platform with 5 mw wind turbine. Ocean Systems Engineering, 1(4):285–296, 2011. [86] Xiaofan Li, ChienAn Chen, Qiaofeng Li, Lin Xu, Changwei Liang, Khai Ngo, Robert G Parker, and Lei Zuo. A compact mechanical power take-off for wave energy converters: Design, analysis, and test verification. Applied Energy, 278:115459, 2020. [87] Daniil Yurchenko and Panagiotis Alevras. Parametric pendulum based wave energy converter. Mechanical Systems and Signal Processing, 99:504–515, 2018. [88] Duong Van Tuan Kiet Gupta, Aakash and Wei-Che Tai. Nonlinear energy transfer of a spar-floater system using the inerter pendulum vibration absorber. Journal of Vibration and Acoustics, page in print, 2023. [89] Aakash Gupta and Wei-Che Tai. Ocean wave energy conversion with a spar-floater system using a nonlinear inerter pendulum vibration absorber. In International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. American Society of Mechanical Engineers, 2023. [90] Klaus F Hasselmann, Tim P Barnett, E Bouws, H Carlson, David E Cartwright, K Eake, JA Euring, A Gicnapp, DE Hasselmann, P Kruseman, et al. Measurements of wind-wave growth and swell decay during the joint north sea wave project (jonswap). Ergaenzungsheft zur Deutschen Hydrographischen Zeitschrift, Reihe A, 1973. [91] Amol Marathe and Anindya Chatterjee. Asymmetric Mathieu equations. Proceedings of the Royal Society A, 462(2070):1643–1659, 2006. * 113 APPENDIX A: BIFURCATION TRACKING ALGORITHM A.1 Algorithm formulation A bifurcation tracking algorithm based on [91] is formulated as follows. First, define a vector function consisting of all the algebraic equations in (2.15) as follows. f ( ˆx; 𝝀) = ℎ0 ( ˆx; 𝝀) ℎ𝑐 1 ( ˆx; 𝝀) ... ℎ𝑐 𝑃 ( ˆx; 𝝀) ℎ𝑠 1 ( ˆx; 𝝀) ... ℎ𝑠 𝑃 ( ˆx; 𝝀) (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) (A.1.1) where ˆx is the vector consisting of all the Fourier coefficients in (2.15) and 𝝀 are bifurcation parameters of interest, e.g., 𝝀 = 𝜂 in Sec. 2.3.3 and 𝝀 = [ 𝑓 , 𝜔]𝑇 in Sec. 2.3.4. Second, define a scalar function 𝑔( ˆx; 𝝀) that outputs the maximum magnitude of eigenvalues of the fundamental matrix of (2.17). A bifurcation point will satisfy the following equations h( ˆx; 𝝀) = (cid:169) (cid:173) (cid:173) (cid:171) f ( ˆx; 𝝀) 𝑔( ˆx; 𝝀) − 1 = 0, (cid:170) (cid:174) (cid:174) (cid:172) (A.1.2) Now suppose a solution ˆx0 to (A.1.2) is found for a bifurcation point 𝝀 = 𝝀0, e.g., by numerical integration. To obtain a neighboring bifurcation point, we use an arc length continuation method. Assume there exists a neighboring solution ˆx1 = ˆx0 + 𝜹 ˆx0 for a neighboring bifurcation point √︁||𝜹 ˆx0||2 + ||𝜹𝝀0||2 << 1, where ||.|| depicting the ℓ2 norm. 𝝀1 = 𝝀0 + 𝜹𝝀0 satisfying (A.1.2) with Imposing the constancy of arc length constraint, we obtain ||𝜹 ˆx0||2 + ||𝜹𝝀0||2 = || ˆx1 − ˆx0||2 + ||𝝀1 − 𝝀0||2 = 𝑠2, (A.1.3) where 𝑠 is a sufficiently small arc length. Therefore, solving (A.1.2) and (A.1.3) together will give a neighboring bifurcation point. 114 A.2 Numerical implementation The bifurcation tracking algorithm is implemented numerically using Newton-Raphson iterations as follows. Suppose a solution ˆx0 for a bifurcation point 𝝀0 is known. Define a new vector function p( ˆx, 𝝀; ˆx0, 𝝀0) = h( ˆx; 𝝀) || ˆx − ˆx0||2 + ||𝝀 − 𝝀0||2 − 𝑠2       (A.2.1) which needs to equal 0 to obtain a neighboring bifurcation point. Suppose two neighboring solutions ˆx𝑘−1 and ˆx𝑘 are found for two bifurcation points 𝝀𝑘−1 and 𝝀𝑘 , respectively. Denote by y = (cid:2) ˆx𝑇 , 𝝀𝑇 (cid:3)𝑇 the solution to (A.2.1). A neighboring solution is initially guessed to be y 𝑔 𝑘+1 = y𝑘 + (y𝑘 − y𝑘−1) = 2y𝑘 − y𝑘−1. (A.2.2) Next, (A.2.1) is solved using Newton-Raphson iterations with the initial guess. The correction to the initial guess is given by y(𝑛+1) 𝑘+1 = y(𝑛) 𝑘+1 − J−1p (cid:16) y(𝑛) 𝑘+1; ˆx𝑘 , 𝝀𝑘 (cid:17) , 𝑛 = 0, 1, · · · with y(0) 𝑘+1 = y 𝑔 𝑘+1 being the initial guess and J is the Jacobian matrix. To calculate J, we use the forward difference method, i.e., the 𝑣𝑡ℎ column of J is J(:, 𝑣) = (cid:16) y(𝑛) 𝑘+1 p + 𝜖e𝑣; ˆx𝑘 , 𝝀𝑘 𝜖 (cid:17) (cid:16) − p y(𝑛) 𝑘+1; ˆx𝑘 , 𝝀𝑘 (cid:17) . (A.2.3) where e𝑣 is the 𝑣𝑡ℎ column in the 𝑁 × 𝑁 identity matrix, 𝑁 being the dimension of y and 0 < 𝜖 << 1. For the calcuations in this study, 𝜖 = 10−5 was taken and the convergence criterion ||𝑦 (𝑛+1) 𝑘+1 || < 10−10 was used. Newton-Raphson iterations will give the value of y𝑘+1, this − 𝑦 (𝑛) 𝑘+1 (along with y𝑘 ) can be further used to calculate y𝑘+2, thereby tracking the bifurcation points. 115 APPENDIX B: JACOBI-ANGER EXPANSION The Jacobi-Anger formulas used for expansion are as follows [61] cos(𝜙0 sin 𝜓) = 𝐽0(𝜙0) + 2𝐽2(𝜙0) cos(2𝜓) + 2𝐽4(𝜙0) cos(4𝜓) + . . . , sin(𝜙0 sin 𝜓) = 2𝐽1(𝜙0) sin(𝜓) + 2𝐽3(𝜙0) sin(3𝜓) + 2𝐽5(𝜙0) sin(5𝜓) + . . . (B.1) (B.2) where 𝐽𝑚 (𝜙0) is a Bessel function of the first kind of order 𝑚. For 𝜙0 = 2.0 radians (115 deg), 𝐽4 (𝜙0) = 0.034 and 𝐽5 (𝜙0) = 0.007. Therefore, only the first two terms are required in the expansion for good accuracy over a wide range of pendulum oscillations. 116