Performance enhancement of galloping-based piezoelectric energy harvesting by exploiting 1:1 internal resonance of magnetically coupled oscillators

In this study, we introduced 1:1 internal resonance in a magnetically coupled 2-degree-of-freedom galloping-based piezoelectric energy harvester to improve the energy-harvesting efficiency. The governing equations for the proposed magnetically coupled aeroelectromechanical system considering the effect of oscillating wake are established using the extended Hamilton principle, and the method of multiple scales is exploited to obtain approximate analytical solutions. Parameter sweeping calculations are conducted to validate the analytical predictions through comparison with analytical solutions, which show good matching. In addition, the overall dynamic behaviors for 1:1 internal resonance under a pure oscillating wake are investigated. Typical nonlinear characteristics such as jump, hysteresis, and frequency synchronization appear in the present system with varying design parameters. In particular, cusp bifurcations are observed under changes in the wind flow proportionality constant and oscillating wake coefficient, while synchronization regimes are observed in the plan of the frequency detuning parameter and oscillating wake coefficient. Chaos phenomena are observed as the gap distance between the bluff bodies decreases to less than 5 mm, which reveals that the frequency lock-in phenomenon can be strengthened by adjusting the magnet distance. From the perspective of output performance, both the voltage and power output results show that the exploitation of 1:1 internal resonance can significantly improve the output performance under a suitable magnet distance.


Introduction
With the dramatic increase in energy consumption due to rapid industrial development, a series of environmental pollution problems has continued to arise over the past few decades. The trend of future human development of renewable and sustainable energy is becoming a worldwide consensus. In recent decades, increasing interest has been shown in energy harvesting from the surrounding environment as an alternative to renewable energy. In particular, the recent rapid developments in wireless sensor networks have created new demands on self-powering technologies, such as portability, miniaturization, high efficiency, and environmental adaptability. In response to these demands, numerous researchers have been devoting their efforts toward the development of various energy-harvesting technologies.
The operating environment of sensor nodes and certain electric devices includes various sources of wasted energy, such as wind [1], mechanical [2][3][4], thermal [5], solar [6], and chemical [7] sources. These can be selected and converted into electric energy for electric consumption devices. Wind energy, as one of the most common wasted energy sources, is abundant and free. A large amount of wind energy has been applied to the traditional rotating wind turbine [8] for commercial power supplies because of its high output and high efficiency in wind energy conversion. On the other hand, the wind turbine has difficulty meeting the requirements of self-powered small-scale wireless sensors and their target electrical devices that need compact structural integration of their power supply components. Therefore, the natural shortage of wind turbines applied to small-scale electric devices has accelerated the study of compact wind energy-harvesting systems. Many of these are based on flowinduced vibrations such as vortex-induced vibration (VIV), galloping, flutter, and wake galloping. Among them, the galloping phenomenon is a representative flow-induced vibration that occurs when the structural system loses its stability. The galloping phenomenon is broadly applied in small-scale wind energy-harvesting systems because it usually exhibits a large oscillation amplitude and wide operating wind speed range compared to VIV, flutter and wake galloping [9][10][11].
Wind energy can be converted into electric energy through three typical transduction mechanisms: electromagnetic, electrostatic, and piezoelectric. Among them, the piezoelectric mechanism is the most popular owing to its high energy density and ease of application [11][12][13]. In the last decade, a variety of piezoelectric wind energy harvesters (PEHs) have been studied from the basic mechanisms of flow-induced vibration to the technologies of improving energy conversion efficiency [14][15][16][17][18]. Among various aerodynamic phenomenon-based wind energy harvesters, galloping-based piezoelectric wind energy harvesters (GPEHs) exhibit excellent performance because of their larger amplitude response and wider operational wind range. A general GPEH is composed of a cantilever beam attached to a piezo-patch and a crosssectionally asymmetric bluff body fixed on its free end. To further improve energy-harvesting efficiency and reduce the critical wind speed, some researchers have been devoting their efforts to optimizing the cross section of the bluff body [19][20][21][22][23][24], exerting a synergistic effect by combining aerodynamic phenomena (VIV, galloping) [25][26][27][28] and increasing the degrees of freedom (DOF) [29,30]. Lan et al. [29] proposed lumped parameter models for two structural configurations with different methods of arranging the bluff bodies under aerodynamic forces. They analytically and numerically compared the output performances between the proposed configurations. The results showed that the second configuration with the wind flow loaded on the outer bluff body can easily and remarkably reduce the critical wind speed and improve the output power. In addition, Hu et al. [10] recently proposed a comb-like beam (CombBeam)based GPEH, which consists of several parasitic beams mounted on a conventional cantilever beam. This structure is regarded as a multiple-degree-offreedom (MDOF) system. They showed that the MDOF CombBeam-based GPEH was superior to the conventional beam GPEH in reducing the cut-in wind speed from 2.24 to 1.96 m/s and enhancing the power by approximately 171.2% with an optimal resistance under a specific wind speed of 3 m/s. Starting from such a linear 2-DOF GPEH, the exploitation of an external nonlinear force in a flow-induced vibration energy harvester was also considered in some works to further enhance the energy-harvesting performance. Wang et al. [31] proposed a double-beam piezomagneto-elastic nonlinear wind energy harvester, in which two piezoelectric cantilever beams are attached to prism-like bluff bodies employing repulsive magnets fixed vertically and oriented to face each other.
Most recently, Sun et al. [9] investigated a magnetically coupled GPEH in a tandem configuration to improve the energetic performance, which consists of two elastic structures with two bluff bodies and a pair of repulsive magnets installed for the generation of mutual interaction between them. This work showed that a maximum overall average output power of 8.09 mW and power improvement rate of 65% were obtained compared to the case without magnetic interaction under an optimal magnetic gap distance. Yao et al. [32] investigated several vibration energy-harvesting systems with optimized L-shaped piezoelectric cantilever beam structures, revealing that the inverted bistable L-shaped beam demonstrated the highest power generation efficiency because this structure can pass through the potential barrier to achieve high energy-harvesting efficiency at either a lower excitation frequency or a higher excitation frequency. It is worth noting that these works have reported excellent energy-harvesting efficiency through various techniques, but they are still insufficient to give a comprehensive interpretation of these nonlinear dynamic systems, which plays an important role in providing design guidance. Moreover, few studies have reported on the effects of the natural frequency ratio between the oscillators on the overall output performance in a two-or multiple-DOF GPEH.
To the best of our knowledge, few study studies have reported on the effects of the natural frequency ratio between the oscillators on the overall output performance in a two-or multiple-DOF GPEH. To date, some systematic studies on internal resonance have been carried out on a general 2-DOF PEH in previous works. In these previous works [33][34][35][36], the dynamical behaviors show dramatic changes when one of the natural frequencies is close enough to the other for internal resonance to occur. Internal resonance has been extensively reported in external forced vibration energy-harvesting systems for output performance enhancement [37], which usually appears when the ratio between two natural frequencies of a given system is close to an integer and various interesting phenomena, such as frequency or (and) lock-in (synchronization), hysteresis, and chaotic behaviors, can result. Inspired by these works on internal resonance in a 2-DOF PEH, this study attempted to exploit 1:1 internal resonance in the development of a magnetically coupled gallopingbased piezoelectric energy harvester (MCR-GPEH) for a significant improvement in output performance.
This study proposed a lumped parameter model for the MCR-GPEH considering the effect of oscillating wake using the Hamilton principle, and a series of comprehensive parametric studies have been analytically and numerically conducted to predict the evaluation of the symmetric dynamic behavior when internal resonance occurs. This work provides a general framework and design guidance for a highefficiency wind energy harvester by exploiting 1:1 internal resonance through controlling the distance. Both the numerical and analytical results showed that a 1:1 internal resonance wind energy-harvesting system is expected to be a promising alternative in the development of a high-efficiency galloping-based wind energy-harvesting system.

Conceptual description
A conventional piezoelectric wind energy harvester based on galloping or VIV aerodynamic phenomena is generally designed with a single cantilever beam and a bluff body fixed on its free end. The wind energy can be converted into structural vibration energy and then into electric energy through the piezoelectric transducer. Various aspects of this conversion have been thoroughly investigated in previous studies, such as fundamental theory construction, novel structure design and aerodynamic energy conversion mechanism.
In this study, we proposed a magnetically coupled 2-DOF galloping-based piezoelectric energy harvester to increase the energy-harvesting efficiency through double-step energy extraction from the incoming wind flow using the tandem arrangement of the bluff bodies, as shown in Fig. 1. The proposed structure is configured with two presumed arbitrary bluff bodies and individually connected to a fixture using an elastic component. Considering the merits of the galloping phenomenon in wind energy harvesting, as introduced in the Introduction section, the cross section shape of Fig. 1 Schematic of the magnetically coupled 2-DOF galloping-based piezoelectric energy harvester the bluff bodies is designed as a circular asymmetric shape for aerodynamic instability. The primary bluff body is arranged in front of the secondary bluff body with a tandem configuration along the wind incoming direction. The wind flow passing the primary bluff body generates wake flow that impacts the secondary bluff body and further gives rise to the variation in its dynamic behaviors. However, the oscillation of the secondary bluff body has a negligible influence on the primary bluff body due to the directionality of the wake flow, as investigated in [38]. In this study, we exploited a pair of mutually repelling magnets attached to the inner side of each bluff body to strengthen the coupling effect of the dynamic behaviors of the two bluff bodies, thus hoping to achieve a synergistic effect.
This study focuses on the investigation of the overall output performance of the two individual elastic structures under various natural frequency differences. The nonlinear characteristics induced by the wind flow and magnetic force lead to various nonlinear dynamic behaviors, such as the jump phenomenon, limit cycle, frequency synchronization and catastrophe bifurcation. The variation in the natural frequency difference between two weakly coupled mass-spring-damper systems results in synchronization and internal resonance, which benefits the wind energy-harvesting performance.

Governing equations
In this section, a mathematical model of the proposed magnetically coupled 2-DOF system is used to predict the dynamic behaviors and output voltage performance. The proposed system consists of two discrete mass-spring-damper components that are coupled with a weak nonlinear magnetic force. The governing equations can be formulated using Hamilton's principle as follows: where t represents time, d T t and d P t represent the variations of the total kinetic and potential energy contained in this system, respectively, and dW nc is the virtual work done by the nonconservative forces. The kinetic energy in the proposed system includes the translational and rotational kinetic energies of the outer and inner bluff bodies. Thus, the total kinetic energy T t in the present system is given by where the overdot of a symbol denotes its derivative with respect to time t; m 1 and m 2 denote the total masses of the bluff bodies (primary and secondary, respectively) and their corresponding attached permanent magnets, respectively; and y 1 and y 2 denote the displacement of the primary and secondary bluff bodies. The potential energy of the proposed system includes several components, such as the electrical, electromechanical, mechanical, and magnetic potential energy. Subsequently, the total potential energy P t of the proposed system is expressed as follows: where k 1 and k 2 represent the spring constants of the primary and secondary springs; c p1 and c p2 represent the capacitance of the primary and secondary piezo components; h e1 and h e2 represent the piezoelectric coupling coefficients; _ k 1 and _ k 2 represent the flux linkage of the equivalent circuit of the piezoelectric components; and U mag represents the potential energy from the magnetic field generated by the two repulsive magnets, which can be described by a dipole-dipole model due to their small size compared to those of the bluff bodies [39].
where l 0 is the vacuum permeability; r BA is the vector from the center of magnet A to the center of magnet B; and m A and m B are the magnetic moment vectors of the outer and inner magnets, respectively. For a permanent magnet, the magnitude of the magnetization vector can be evaluated using the magnet's residual flux density B r , that is, m ¼ B r =l 0 . Under a series of mathematical manipulations, the magnetic potential energy can be expressed as where y 12 ¼ y 1 À y 2 , referring to [9], m 1 and m 2 are the magnitudes of the magnetization vectors, v 1 and v 2 are the volumes of the primary and secondary magnets, and D m is the minimum gap distance between the center points of the two magnets. Taking the derivative with respect to y 12 of the magnetic potential energy and then expanding the resulting equation in a Taylor series at y 12 ¼ 0, the magnetic force can be obtained as The virtual work done by the nonconservative forces is given by where R L1 and R L2 represent the load resistances for the primary and secondary external circuits, respectively; c 1 and c 2 represent the damping constants; h w i denote the wake force coefficients, assuming both are equal to H w ; and F g1 and F g2 are the aerodynamic transverse forces acting on the bluff bodies, which can be expressed as [26] where the indices i = 1 and 2 represent the corresponding symbols for the primary and secondary bluff bodies, respectively; q a represents the air density; H i and U i are the height of the bluff bodies and the wind speed applied to each bluff body; C gi represents the dimensionless coefficients of the aerodynamic transverse force, which is a function of the attack angle according to the assumption of a quasi-steady state of the bluff bodies in the wind flow field according to previous studies [22,40]. C gi is an inherent characteristic related to the cross-sectional shape of the bluff body and the Reynolds number, which can be obtained through series expansion of the transverse aerodynamic force at different rotation angles of the bluff body. It can be expressed as , and A 5i are the aerodynamic coefficients. In this study, a third-order expansion is adopted to describe the nonlinear dynamic behaviors resulting from the aerodynamic force.
Subsequently, we performed several manipulations of Eqs. (2), (3) and (7) using variational calculus, substituting the results into Eq. (1) to gather the coefficients of dy 1 , dy 2 , dk 1 ðtÞ and dk 2 ðtÞ. Furthermore, we conducted a normalization manipulation for calculation convenience. Thus, a set of governing equations can be obtained as where the definitions of the normalized parameters in Eq. (9) are listed in Table 1. The set of governing equations given in Eq. (9) indicates that the two primary and secondary oscillators are mechanically coupled to each other through a nonlinear magnetic force. Additionally, the oscillating aerodynamic wake force generated by the adjacent bluff body during oscillation is considered to avoid the loss of generality of the mathematical modeling, although it produces only a tiny disturbance to the dynamic behaviors of the bluff body compared to that of the magnetic force. It should be noted that the oscillating wake force acting on the bluff bodies is not equal due to the unilaterality of the incoming wind flow. Accordingly, an order selection considering the wake force scale should be conducted to find an analytical solution using the perturbation method.

Analytical solutions based on the method of multiple scales
In an attempt to understand the nonlinear dynamic behaviors of the present system, some semianalytical techniques, such as the ratio approach [41] and parameterized Adomian decomposition method [42], can be employed to obtain the overall analytical solutions. In particular, the internal resonance between the individual oscillators is investigated in this study; therefore, the perturbation method of multiple scales [43] is utilized to obtain the approximate analytical solutions for Eq. (9). The time dependence can be expanded into two time scales in the form T 0 ¼ t (fast time scale) and , where e is a scaling parameter. The new time derivatives with respect to t can be expressed as: where D n ¼ o oT n and O e 2 ð Þ denotes the terms higher than the second power of e. The asymptotic series solutions are expanded up to the first order, such that where the symbols ð Þ 0 and ð Þ 1 represent the e 0 -and e 1order solutions for the dynamic responses of the primary and secondary oscillators. The oscillators in the present system are assumed to be weakly coupled by the external magnetic force and oscillating wake force. An ordering for the terms in Eq. (9) can be determined by scaling the viscous damping, magnetic coupling, electromechanical coupling, and aerodynamic force in the same order. In other words, the scaled relevant coefficients or terms can be expressed as It should be noted that the oscillating wake coupling coefficient generated by the secondary oscillator is defined at a higher order of e than e 2 considering the fact that the effect of the wake generated from the secondary oscillator to the primary bluff body is much weaker than that generated from the primary oscillator to the secondary oscillator. Substituting the scaled terms into Eq. (10) and expanding the governing equations in terms of Eqs. (10) and (12), we can obtain the ordered equations after collecting the coefficients of e 0 and e 1 , which yields the following results: e 0 order equations:

Mass normalized electromechanical coupling coefficients
Normalized wake coupling coefficientsh w i ¼h w i =m i Normalized magnetic force coefficients for primary oscillator Normalized magnetic force coefficients for secondary oscillators 1 e 1 order equations: The solutions of the zeroth-order Eqs. (13)-(16) can be written as: where cc denotes the complex conjugate of the preceding term; A 1 and B 1 denote the unknown complex-valued functions for primary and secondary oscillators, respectively; the overbar denotes the complex conjugate; substituting the zeroth-order solutions into the right-hand side of the first-order Eqs. (17)(18)(19)(20) and expanding the equations, it can be observed that frequency modes with a frequency ratio of 1:1 are coupled via the oscillating wake force and the linear terms of the magnetic force. Additionally, it can also be observed that the frequency modes with a frequency ratio of 1:3 are coupled via cubic nonlinearities of the magnetic force. In this study, the 1:1 internal resonance plays a decisive role in the improvement of dynamic behaviors and the output performance of the present system, referring to the results of a preliminary numerical analysis. After eliminating the undesired nonsecular terms and collecting coefficients of harmonics, the following can be obtained: where NST stands for the terms that do not produce secular terms. Then, it is necessary to introduce the 1:1 internal resonance condition as where r is the internal resonance detuning parameter, which represents the closeness of the two natural frequencies of the primary and secondary oscillators. Substituting x 2 and x 1 in Eq. (27) into Eqs. (25)- (26) and after mathematical manipulation, the simplified equations can be expressed as Eliminating the terms that lead to secular terms in the above equations and expressing the complexvalued function in polar form where a and b represent the unknown amplitudes for the primary and secondary oscillators, respectively, and b 1 and b 2 represent the phases of the corresponding dynamic responses. Next, separating the real and imaginary parts in Eqs. (28)-(29) yields the modulation equations as where the time derivatives of b 1 and b 2 in Eq. (32) and (34) can be vanished in terms of the relationship of c 0 ¼ r þ b 0 2 À b 0 1 , where the prime means the time derivative of T 1 , and thus the modulation equations with 4-dimensional phase space (a; b; b 1 ; b 2 ) can be reduced to a 3-dimensional phase space (a; b; c). Additionally, the steady-state responses can be determined by setting the time derivatives in the modulation equations equal to zero. The resulting equations are a set of transcendental equations that are difficult to solve by an analytical method; however, it is fairly straightforward to solve them numerically [44] for semianalytical solutions. Subsequently, we can obtain the output voltage generated by the piezoelectric element using Eqs. (23) and (24), where both homogeneous and nonhomogeneous solutions are contained in them. It is known that the homogeneous solutions would converge to zero on account of the electric energy dissipation in the external resistor over time. Consequently, we can find the output voltage amplitude through the nonhomogeneous terms, and the resultants can be expressed as a function of the amplitude of the corresponding oscillators as Subsequently, the output power can be described as According to the expressions of output voltage and power, we can compare the energy-harvesting performance between the internal resonance case and the noninternal resonance case.

Oscillating wake-induced internal resonance
In this section, we first investigated the electric characteristics of the external circuit and some fundamental aerodynamic features as a preliminary investigation for the further examination of coupled dynamic behaviors. Additionally, the oscillating wake generated from the primary bluff body has an effect on the dynamic behaviors of the secondary bluff body. To validate the proposed analytical model and give a comprehensive understanding of the oscillating wakecoupled system, we examine the dynamic behaviors considering the pure oscillating wake coupling effect by assuming the magnet gap distance to be infinite. Subsequently, some bifurcation analyses are carried out to clarify the inherent characteristics of the internal resonance. The relevant parameters and the corresponding values used in this study are listed in Table 2. In a wind energy-harvesting system, the critical wind speed is one of the most important factors in evaluating the output performance. Here, we first investigate and clarify the influence of the external load on the critical wind speed in terms of secular equations. The energy-harvesting system would lose stability when the linear mechanical damping changes to be negative. Nevertheless, the resultant damping term is determined by the combined effect of mechanical damping and electric damping in an electromechanically coupled system. From the modulation equations of Eqs. (31) and (33), the difference in the phase angle between the two oscillators can be zero or 180°before the oscillators start to vibrate, which cause the triangular functions to be zero. Therefore, the criteria equation for wind speed can be derived from Eqs. (31) and (33) as Furthermore, the amplitudes a and b can be calculated as In this equation, A i3 is always negative, so it is clear that the above equation is valid only if Therefore, the critical wind speed can be expressed as According to the expression of critical wind speed, we can investigate the relationship between the external resistance and the critical wind speed. As shown in Fig. 2, the critical wind speed increases as the external resistance increases and then decreases with a further increase in the external resistance, which matches the results in [46]. The up-and-down phenomenon of critical wind speed versus external resistance results from the electromechanical damping terms of E i2 , which cause a maximum value damping at a resistance value of 0.176 MX for the primary oscillator. The resistance value is also regarded as an optimal value to obtain a maximum output power in the external circuit. As shown in Fig. 2a, the electromechanical coupling strength also shows a significant influence on the critical wind speed. It can be observed that the strong coupling effect between the electric domain and mechanical domain results in an (a) (b) Fig. 2 Variation in critical wind speed versus external resistance with different a electromechanical coupling coefficients and b detuning parameters increase in the critical wind speed. On the other hand, the detuning parameter r describes the difference between the natural frequencies of the primary and secondary oscillators, which causes a change in the secondary natural frequency, as shown in Fig. 2b. This further indicates that an increase in natural frequency is able to pull down the critical wind speed. In a general energy-harvesting system, impedance matching between mechanical and electric systems is essential to obtain a maximum output power. Mechanical energy would be dissipated during the energyharvesting operation because the mechanical energy is transformed into electric energy. However, the increase in electric damping results in a dramatic decrease in the dynamic response of the oscillator, as shown in Fig. 3, especially in the low wind speed range. This result gives a pertinent interpretation of the fact that reaching the critical wind speed is postponed as the external resistance approaches the optimal value, as shown in Fig. 2. For a consistent evaluation of the effect of internal resonance and other key parameters, the external resistance is treated as infinite, in other words, an open circuit in the next section.

Parametric analyses under the wake effect
In this subsection, we investigated the dynamic responses of both oscillators under no magnetic force coupling by setting the distance between the two magnets to infinity. The oscillating wake generated by the primary bluff body periodically acts on the secondary bluff body; however, it is not easy for the oscillating wake of the secondary bluff body to inversely affect that of the primary bluff body. This phenomenon is analytically and numerically investigated, and the results are shown in Fig. 4a. In this numerical calculation, we adopted a means of wind speed sweeping to illustrate the evolution of dynamic behavior with time. By comparing the numerical results to the analytical results, good matching can be seen in a wind speed range from 0 to 10 m/s. The strength of the incoming wind exerted on the secondary bluff body tends to be weaker due to the wind flow obstruction of the primary bluff body In other words, the actual wind speed loaded to the secondary bluff body is smaller than that to the primary bluff body. This fact necessitates a function relationship between the wind speeds of U 1 and U 2 . We assume that U 2 is linearly related to U 1 and expressed as U 2 ¼ gU 1 . Here, g is a wind speed proportionality constant to quantitatively express the reduction in wind speed, which is always smaller than one. Subsequently, we analytically investigated the amplitude responses of the secondary bluff body versus wind speed under different g with H w = 0 Ns/m and r = 0 rad/s, as shown in Fig. 4b. The wind speed range in this study ranges from 0 to 10 m/s based on the fact that the usual wind speed in an urban environment is assumed to be within this range. The amplitudes of the bluff body increase along the axis of the wind speed for all the various g cases. However, the critical wind speed moves to the higher range as g decreases. The subplot in Fig. 4b shows the variation in the dimensionless critical wind speed U cr (U cr =U 0 ) over the wind speed proportionality constant g, where U 0 denotes the critical wind speed when g is zero. It further shows that U cr is very sensitive to the variation of g at a low range of g, which demonstrates a dramatic increase as g approaches zero. A natural difference in the dynamic responses occurs in the secondary bluff body when the wake coupling effect is accounted for in the dynamical analysis. Response bifurcation and multiple solution phenomena start to appear when H w is given a small value of 0.002 Ns/m, as shown in Fig. 5. Here, two stable branches and one unstable branch can be determined through the stability evaluation based on the eigenvalues of the Jacobian matrix of the righthand sides of Eqs. (31)- (34), as stated in [47]. The cross point between the bifurcation curve and the dynamic response curve is defined as the saddle node. It can be clearly observed that the critical wind speed is pulled down to the lower wind range referring to that case without consideration of wake oscillation. In addition, it is also observed that the saddle node can be significantly reduced from 9.2 to 3.8 m/s as g increases from 0.08 to 0.16, and the second stable branch tends to be close to the first stable branch as g increases, as illustrated in Fig. 5a and b.
The influence of the detuning parameter r on the dynamic response of the secondary bluff body is also investigated in the presence of oscillating wake coupling, as shown in Fig. 6. The graphs depicted similar dynamic behaviors with the cases shown in Fig. 5; however, we can see that the variation of r exhibits few effects on the saddle node position, which is much different from the varying case of g. To further evaluate the analytical results, a series of numerical simulations were conducted by sweeping the parameter over time. A jump phenomenon occurs at the saddle node point when the wind speed increases from zero to a higher wind speed range in both the r = 0.15 and r = 5 cases, as shown in Fig. 6a and b. It is worth noting that the maximum amplitude appears close to the first stable branch at a smaller r and tends to approach the second stable branch when a larger r is introduced. On the other hand, the amplitude responses decrease along the routine of the first stable branch when the wind speed decreases backward. From the above comparisons, good matching between the analytical and numerical results is achieved, which further validates the proposed prediction model.

Perspective in saddle-node bifurcation
As previously mentioned, the saddle node point appears as we vary the associated parameters, and the saddle node plays an important role in the dynamic response evaluation for a nonlinear coupled system. Accordingly, a systematic analysis of the overall bifurcation curve is necessary to reveal the mechanism contained in the wake-coupled oscillators. Response bifurcation usually occurs at positions where the derivative of the response to the varying parameter is zero. A 3-dimensional graph for the amplitude of the secondary oscillator in terms of two varying parameters can be identified based on the resultant polynomial equation. It is expected that a surface fold phenomenon of the dynamic response would occur (a) (b) Fig. 7 Saddle node bifurcation of the oscillating wake-coupled nonlinear oscillator in the a g À H w plane and r À H w plane due to the nonlinearity of the coupled parameter and aerodynamic force. For a coherent exhibition of the saddle-node bifurcation curve, we reduce the polynomial equation including 3 varying parameters to that including 2 parameters by vanishing the amplitude of the secondary oscillator using the original polynomial equation and its derivative with the amplitude of the secondary oscillator. As shown in Fig. 7, a series of saddle-node bifurcation curves are generated in the g À H w and r À H w plane, where the stable and unstable fixed points collide. Additionally, a cusp bifurcation occurs when two saddle-node bifurcation curves meet, as shown in Fig. 7a and b. The folds degenerate at these cusp points; that is, the frequency synchronization tends to be difficult when the varying parameter passes the cusp point. This shows that the cusp point moves to higher values of g and H w as r increases, as illustrated in Fig. 7a. In Fig. 7b, the shadow regions are bounded by the saddle node curves and saddle-node bifurcation curve for a saddle and two saddle cycles (SSN). We can see from this graph that the shadow region is symmetric with the axis of r = 0, which indicates that the bounded synchronization region is strongly dependent on the natural frequency difference of the two oscillators. Furthermore, g has a significant influence on the area of the bounded regime, namely, the increase in g causes area expansion of the bounded regime.

Magnetically force-induced internal resonance
The fundamental dynamic behaviors of the oscillators undergoing oscillating wake coupling have been parametrically investigated in the previous section without consideration of actual physical conditions. In galloping-based energy harvesting, an asymmetric bluff body is usually adopted to generate the transverse aerodynamic force, which leads to vibration of the oscillator as the system loses its stability. In the galloping phenomenon, the wake shed from the asymmetric bluff body is easily broken as wind flow is strongly separated for the edge of the bluff body, in which the wake intensity would be significantly reduced. As an alternative, we consider utilizing a weak magnetic force effect to bring out the internal resonance to improve the energy-harvesting performance. In this section, we concentrate on the clarification of the synchronization mechanism and evaluation of the improved performance of the magnetically coupled wind energy harvester under internal resonance.

Synchronization analysis
To further validate the analytical solution in the magnetically coupled case, we compared the analytical predictions and numerical solutions by examining the effects of some key parameters on the dynamic responses of the primary and secondary oscillators. Herein, the numerical solutions were obtained by sweeping the parameters with varying time. Figure 8 shows the dynamic responses for each oscillator when the detuning parameter changes with fixed associated parameters, in which a very good matching between the analytical predictions and numerical solutions can be observed throughout the concerned range of detuning parameters. This graph shows that a slight amplitude decrease in the primary oscillator (Fig. 8a) corresponds to a significant amplitude increase in the secondary oscillator (Fig. 8b) as r approaches zero, where the 1:1 internal resonance occurs. This phenomenon indicates that an energy exchange occurs between the two modes of oscillation during their oscillation. It is worth noting that not only stable solutions but also unstable branches of the secondary bluff body were observed in the analytical prediction, as shown in Fig. 8b. The unstable branch is strongly dependent on the wind speed proportionality constant g; in other words, this unstable branch is directly induced by the wind flow applied to the bluff body.
To represent the dynamic strength distribution over time at various frequencies and examine the frequency or amplitude synchronizing regime, spectrograms via fast Fourier transformation are adopted, as shown in the second row of Fig. 8. It can be observed that a concentrated and consistent energy distribution of the primary oscillator is located near 20 Hz, and a similar phenomenon can be found in the spectrogram of the secondary oscillator even though some energy distribution appears in other frequency regimes. This spectrogram indicates that the oscillating frequency of the secondary oscillator is locked on that of the primary oscillator. However, the secondary natural frequency varies by adjusting the detuning parameter r, which gives rise to some subharmonic or superharmonic responses in the secondary oscillators, as shown in the spectrogram of the secondary oscillator. Obviously, the energy distribution level is higher as r approaches zero, which reflects the advantages of using 1:1 internal resonance to improve the energy conversion efficiency from wind energy to structural vibration energy and further enhancing the wind energy-harvesting performance.
As the wind speed proportionality constant g increases to 0.16, the wind flow-induced response curve changes from an unstable branch to a stable branch, as shown in Fig. 9b. In this case, the dynamic response of the secondary oscillator exhibits a combination of the two stable branches. It also shows that the secondary oscillator amplitude is enlarged as r approaches zero owing to the frequency lock-in phenomenon, as shown in the spectrogram of the secondary oscillator. Compared to the case with g = 0.08, the increase in the strength of the wind flow applied to the secondary oscillator promotes the occurrence of galloping-induced vibration and further becomes the dominant response. It can be observed in the secondary spectrogram that the main energy concentration varies over time following the increase Fig. 9 Dynamic resonance comparisons between the analytical results and numerical results versus varying natural frequency differences and spectrogram analysis for the a primary oscillator, b secondary oscillator and spectrograms for each oscillator at D m = 20 mm, H w = 0.08, and U 1 = 10 m/s at g = 0.16 in r with time variation, which means the secondary oscillator oscillates with its own ascendant natural frequency. Subsequently, some slight changes in the oscillating frequency of the primary oscillator occur when g increases to 0.16, as shown in the spectrogram of the primary oscillator.
As mentioned in the previous section, g plays an important role in the dynamic response of the secondary oscillator. In a practical situation, g is usually a very small value due to the strength deduction of wind flow, as experimentally investigated in [38]. Accordingly, it becomes necessary to improve the energy-harvesting performance by exploiting the magnetic force for the synchronization phenomenon. However, an appropriate magnetic force is crucial to cause the synchronization phenomenon to happen, since a strong mode coupling results in a structural binding phenomenon, which causes a 2-DOF system to collapse to a 1-DOF system. On the other hand, an extremely weak coupling results in invalidation of mode coupling and thus causes the occurrence of synchronization to be difficult. In this study, we investigated the influence of magnet distance on the responses of both oscillators. Subsequently, we investigated the dynamic responses of both oscillators with a small magnet distance for strong mechanical coupling. The results are illustrated in Fig. 10, in which it can be observed that both the primary and secondary oscillators are totally bound to each other, and they have almost the same dynamic behaviors over the varying frequency range. Based on the prediction accuracy point, the analytical prediction results shown in Fig. 10a and b do not give an accurate amplitude matching that of the numerical results. The difference between the numerical results and analytical prediction occurs because the chaos phenomenon occurs in this strongly coupled system as the detuning parameter r varies, and chaos itself is difficult to predict. In the spectrogram of each oscillator in the second row of Fig. 10, it can be identified that a frequency evolution occurs over time, especially when the frequency evolution is irregular and infinite in the chaos regions. The other regime, except the chaotic regime, also shows a quasi-chaos state even though the dominant frequencies exhibit the same trend of gradually increasing for both of the oscillators.
To further clarify the influence of the magnetic force strength on the system dynamic responses, we conducted an analytical investigation and compared the dynamic behaviors of the oscillators according to various magnet distances. The results are shown in Fig. 11, in which it can be observed that the energy bound phenomenon occurs at a magnet distance of 20 mm, where an ascending amplitude of the primary bluff body corresponds to a descending amplitude of the secondary oscillator. Additionally, in this case, the wind flow-induced response is unstable. As the magnet distance increases to 40 mm, the energy bound effect tends to recede, and a combined dynamic response occurs. The wind flow-induced response becomes dominant as the magnet distance further increases to 60 mm. In summary, an appropriate magnet distance is significant in improving the effect of energy transformation between dynamic modes.

Performance evaluation
The 1:1 internal resonance is adopted for the mode synchronization between the oscillators, and thus the output performances are expected to be effectively improved. In this study, a maximum power based on an optimal resistance is assumed to evaluate the output performance.
The total RMS voltage is calculated by combining the output voltages from the primary and secondary oscillators according to Eq. (35) with an optimal resistance of 0.18 MX. Subsequently, the total maximum average power can be calculated according to Eq. (36). As shown in Fig. 12, the output voltage and average power undergo a dramatic decrease of 66% and 175%, respectively, as r approaches zero (that is, the difference between the two natural frequencies is close to zero) at an extremely small magnet distance of 15 mm because a collapse of the number of degrees of freedom occurs when the coupling force exceeds a throttle value. The total RMS voltage with a maximum value of approximately 60 V and an average power of 0.021 W with two peaks occurs near the position when r is zero as the magnet distance increases to 20 mm. Subsequently, as we continually increase the magnet distance to 25 mm, the two electric output peaks change to one single peak near the position where r equals zero, as shown in Fig. 12.
In summary, a significant enhancement of both the total RMS voltage and average power can be achieved by exploiting nonlinear magnet force coupling. The 1:1 internal resonance of two oscillators leads to frequency synchronization as the natural frequencies of the two oscillators are close to each other, further promoting enhancement of the output performance.

Conclusion
In this study, we proposed a 1:1 internal resonance technique to improve the output performance in a magnetically coupled wind energy harvester with two tandem oscillators arranged along the wind flow direction. The proposed technique was numerically and analytically verified to be advantageous, while a significant improvement in the total RMS voltage and average power is readily achieved as frequency synchronization occurs when the two natural frequencies equal each other. Based on the results obtained, several important conclusions can be deduced as follows: 1. The magnetically coupled governing equations representing the transverse dynamics of the bluff bodies and electric circuitry are derived based on the extended Hamilton principle. Comparisons between the analytical prediction solutions and numerical results are conducted and exhibit very good matching between them. 2. Several saddle-node bifurcation analyses are conducted to reveal the nature of the system with the g À H w and r À H w plane. The cusp bifurcation and synchronization region are identified based on the bifurcation analyses. 3. The frequency synchronization under magnetic coupling can be triggered as the natural frequencies of the oscillators approach each other under a suitable magnet distance, as verified both in analytical prediction solutions and numerical results. 4. The total RMS voltages and average power generated by both of the oscillators show a significant improvement when the magnet distance is larger than 15 mm, which benefits energy-harvesting performance.
This study provides an important guideline for the design of a magnetically coupled high-efficiency 2-DOF wind energy harvester using 1:1 internal mode resonance. It is expected that the proposed method has potential in the future design of compact and efficient wind energy harvesters, with broad and practical industrial applications, such as wireless sensor networks.