Secondary Resonance Energy Harvesting with Quadratic Nonlinearity

Piezoelectric energy harvesters can transform the mechanical strain into electrical energy. The microelectromechanical transformation device is often composed of piezoelectric cantilevers and has been largely experimented. Most resonances have been developed to harvest nonlinear vibratory energy except for combination resonances. This paper is to analyze several secondary resonances of a cantilever-type piezoelectric energy harvester with a tip magnet. The conventional Galerkin method is improved to truncate the continuous model, an integro-partial differential equation with time-dependent boundary conditions. Then, more resonances on higher-order vibration modes can be obtained. The stable steady-state response is formulated approximately but analytically for the first two subharmonic and combination resonances. The instability boundaries are discussed for these secondary resonances from quadratic nonlinearity. A small damping and a large excitation readily result in an unstable response, including the period-doubling and quasiperiodic motions that can be employed to enhance the voltage output around a wider band of working frequency. Runge–Kutta method is employed to numerically compute the time history for stable and unstable motions. The stable steady-state responses from two different methods agree well with each other. The outcome enriches structural dynamic theory on nonlinear vibration.


Introduction
Applying pressure to piezoelectric materials causes them to produce an electric charge. The piezoelectric material would contract or expand in an electric field. Since their discovery, the principles of the piezoelectric effect have caught many scientists' attention [1,2]. Owing to the piezoelectricity, piezoelectric energy harvesters are invented to transform mechanical vibration into electrical energy. This kind of electromechanical device can provide a potential solution to capture mechanical energy for operating self-sustained low-power electronics. The technology has not only monetary gains in the reduction of the maintenance costs by replacing batteries, but also ecological ramifications in the reduction of the chemical waste [3].
However, the voltages generated by the piezoelectricity and its inverse are usually small [4,5]. More voltage output is desirable in piezoelectric instruments, including ultrasonic medical transducers, ultrasonic cleaning components, ultrasonic welding components, and accelerometers. To this end, cantilevered piezoelectric energy harvesters have been experimentally investigated in a wide range [6]. Linear energy harvesters efficiently work merely at the natural frequency. The piezoelectric cantilever can be designed to match the resonant frequency for a large response amplitude [7]. Such a design has a single peak around a narrow band of frequency. The amount of generated voltage is quite low when the energy harvester works off resonance. The large voltage output around a wide band of operating frequency has been of great concern [8].
(3) The subharmonic resonance is another way to improve the bandwidth of working frequency. The in-depth analyses can better explain those results observed in experiments, and approximately forecast the working conditions for energy harvesters based on stability theories.
(4) Some interesting phenomena, such as the period-doubling and the quasiperiodic motions, are numerically predicted in the hope of more experimental computations. The hard piezo-materials with higher mechanical quality factors might be beneficial for experimental verification and real use.
The rest is structured as follows. The nonlinear model of a cantilevered piezoelectric energy harvester subjected to a tip magnetic force is established in Section 2. Its natural frequencies and modal functions are calculated and used in the application of the improved Galerkin method. Section 3 analyzes the combination resonances based on a truncation system. The method of two times scales is employed to deal with the governing equation and its stability condition. Two subharmonic resonances are analyzed in Section 4. Section 5 explains a previous experiment by others, and discusses the effects of damping coefficient, excitation amplitude, and maximal voltage on secondary resonance energy harvesting with quadratic nonlinearity. Section 6 ends this paper with some concluding remarks.

The Equation of Absolute Motion
The considered electromechanical device is primarily composed of a cantilevered structure ( Figure 1a) with a general electric circuit (Figure 1b). The slender and uniform cantilever, clamped at a framework, is assumed to be of a rectangular cross-section. A tip-proof mass is attached to adjust the natural frequency for real use. A cylindrical permanent magnet is fixed at the right end of the beam-type structure. The magnet, modeled as a portion of the tip lumped mass, is perpendicularly repelled by a similar cylindrical magnet which is attached to the framework. The structural parameters are the same compared to [3], as described in Table 1. In addition, an extremely thin piezoelectric ceramic patch covers the whole surface of the cantilever. Two electrodes cover the opposite transversal faces of the piezoceramic (PZT) layer. The two electrodes will be connected to a resistive load. When the piezoelectric device is subjected to mechanical stress, an electric charge is generated. As illustrated in Figure 1b, the circuitry includes a resistive load as well as an internal capacitance because the electrode pair brackets the PZT layer. According to Kirchhoff laws, the circuit equation can be expressed as [25,26] ( ) ( ) ( ) 3 eq 2 0 eq d , where 0 < X < L, the piezoelectric energy harvester offers the current source, Ceq represents the equivalent capacitance, Req represents the equivalent resistance, Θ represents the electromechanical coupling coefficient, V(T) represents the voltage output across the resistive load at the time of T, X is the spatial coordinate and U(X, T) represents the transversal displacement from the horizontal axis where the clamped end of the cantilever is the original point, as depicted in Figure 1.
The Kelvin-Voigt damping model is applied to contribute to the material damping effect of the system of viscosity coefficient Λ. A tip proof of mass mt (together with the tip magnet) is attached to  It is assumed that the two applied cylindrical magnets are both with height H and radius R. The initial distance between them is symbolically denoted by D while the entire structure is just equipped. ρ is the residual magnetic flux density. S is the common area between the two magnets. µ is the permeability of air. The magnetic force between the two cylindrical magnets was formulated in the form of [3,7,32] In addition, an extremely thin piezoelectric ceramic patch covers the whole surface of the cantilever. Two electrodes cover the opposite transversal faces of the piezoceramic (PZT) layer. The two electrodes will be connected to a resistive load. When the piezoelectric device is subjected to mechanical stress, an electric charge is generated. As illustrated in Figure 1b, the circuitry includes a resistive load as well as an internal capacitance because the electrode pair brackets the PZT layer. According to Kirchhoff laws, the circuit equation can be expressed as [25,26] where 0 < X < L, the piezoelectric energy harvester offers the current source, C eq represents the equivalent capacitance, R eq represents the equivalent resistance, Θ represents the electromechanical coupling coefficient, V(T) represents the voltage output across the resistive load at the time of T, X is the spatial coordinate and U(X, T) represents the transversal displacement from the horizontal axis where the clamped end of the cantilever is the original point, as depicted in Figure 1.
The Kelvin-Voigt damping model is applied to contribute to the material damping effect of the system of viscosity coefficient Λ. A tip proof of mass m t (together with the tip magnet) is attached to the cantilever of mass m b , length L, Young's modulus E, and area moment of inertial I. The internal bending moment of the cantilever can be represented as [25,26] M(X, T) = EI ∂ 2 U(X, T) where H( ) is the Heaviside step function, which is defined as the integral of the Dirac delta function. If 0 < X < L, the equation of absolute motion is where g is the acceleration of gravity, and the Dirac delta function satisfies The following is the boundary conditions at X = 0 or L. For a slender cantilever, the small tip magnet can be viewed as a point mass, and the magnetic force serves the nonlinear boundary condition at the right end of the cantilever. At its left end, the transversal displacement of the base is assumed to be Bsin(ΩT). Then, the boundary conditions are [25][26][27][28] U(X, T) = B sin(ΩT); ∂U(X, T) Materials 2020, 13, 3389 5 of 24 at X = 0, and at X = L. Equation (7) offers the nonlinearity, which is explicitly dependent on time.
Introduce dimensionless parameters as follows.
where e is the unit voltage (one volt). The dimensionless value η represents the ratio of the tip mass to the mass of the cantilever. From Equations (2), (4), and (8), the application of the Euler-Bernoulli beam theory yields the dimensionless governing equation, which is represented algebraically as ..
where 0 < x < 1, the prime denotes the (partial) derivative with respect to x, and the overdot denotes the (partial) derivative with respect to t. The corresponding boundary conditions are simplified as where the dimensionless magnetic force function is

The Equation of Relative Motion
The standard governing equation of relative motion round its equilibrium will be established. The equilibrium position can be analytically derived from Equation (9). To this end, introduce a new coordinate transformation as follows [27].
where A e satisfies Substituting Equation (12) into Equation (9) yields the standard governing equation. ..
Materials 2020, 13, 3389 6 of 24 The Taylor series expansion method is valid to deal with the nonlinear terms in the boundary condition (15).

Natural Frequencies and Modal Functions
In general, the Galerkin method involves the modal functions of the linearized equation. In the following, the natural frequencies and modal functions are both to be determined. Some of their properties are to be shown. To this end, the standard governing equation goes into a linear form that is independent of the viscoelastic damping, electrical coupling, and external excitation.
This linear partial differential Equation (18) has been solved [27]. It can be expressed as where n = 1, 2, 3, . . . and the modal functions satisfy the following equation.
where the natural frequencies are obtained from For convenience, the nth modal function with condition φ n (1) = 1 is in the form of Here, the integral of the product of two different modal functions satisfies Besides, two other integrals are following.
Equations (23)- (25) is useful in the use of the Galerkin method. where Substituting Equations (24) and (25) into Equation (26) gives According to the boundary conditions (15), the substitution of z (1, t) + α The following is the second step, choosing appropriate trial functions. As a general rule, it is convenient to assume expansions for the relative displacement in the form of the weighted sum of linear modes (19). To perform calculations, N-term truncation that has a finite number of terms is used.
The truncation with the use of Equation (27) makes the electrical component of Equation (14) become Meanwhile, the mechanical counterpart (29) becomes ..
where n = 1,2,3, . . . N, and Equations (31) and (32) are the discrete truncation system that is being sought and to be analyzed.

The Multiscale Analyses
To investigate combination resonances, at least two-term truncation is in need. When N = 2, the governing equations become .. ..
As for a weakly nonlinear energy harvester, the solution to Equation (34) is sought in the form of where ε is a small parameter, the new time variable T 0 is equal to t, and T 1 = εt, which is a slower time scale than T 0 . T 1 characterizes the modulation of the amplitudes and phases due to viscoelastic damping and possible resonances. The derivatives with respect to time t become where When the frequency of the external excitation is away from the natural frequency, the effect of the excitation is small unless its amplitude is hard. To examine subharmonic and combination resonances, the excitation should appear in the lowest-order perturbation equation. To this end, it is necessary to reorder the electromechanical coupling coefficient, viscoelastic damping, and external excitation as follows.
Substituting Equations (35)-(38) into Equation (34) and equating the coefficient of ε on both sides yield Equating the coefficient of the square of epsilon on both sides can give The general solution of Equation (39a) is where Y 0 (T 1 ) is a function controlled by initial conditions. For the steady-state motion, the general solution of Equation (39a,b) is taken in the form of where cc is the complex conjugate of all the preceding terms, and Then, the second-order Equation (40a) becomes whose solution is In the next sections, two subharmonic resonances and combination resonances will be analyzed on the basis of Equation (40b).

The Stability Analyses
In the beginning, the summed harmonic resonance is discussed. Introduce the detuning parameter according to Hence, The second-order Equation (40b) becomes where NST stands for terms that do not produce secular terms. The corresponding solvability condition can be obtained by eliminating secular terms.
Its solution is in the form of where r is the real characteristic value. Putting Equation (50) into Equation (49), and separating real and imaginary parts yield and further from Equation (51a,b), (52) As we know, the system is stable if and only if all eigenvalues are negative real numbers [31]. Otherwise, the electromechanical device is unstable. The critical situation is that the characteristic value is equal to zero. Consequently, the boundary of the stable region satisfies Separating real and imaginary parts of Equation (53) yields The stability boundary can be simplified from Equations (54b) and (55).
In general, the combination resonance of summed type occurs when which is obtained on the basis of Equation (43)  Near the considered combination resonance, Figure 2 plots the stability boundaries (56) for different excitation and viscoelastic damping. The dotted, solid, dashed and dotted-dashed lines represent α = 0, α = 0.001, α = 0.002, and α = 0.003, respectively. The considered parameters are c 0 = 0.5084, d = 0.15, h = 0.004167, η = 0.02, c = κ = 0.000001, and θ = 0.0002. It is revealed that the combination resonance might theoretically occur with a small damping, a large enough excitation amplitude, and a suitable excitation frequency close to the sum of the first two natural frequencies.
combination resonance might theoretically occur with a small damping, a large enough excitation amplitude, and a suitable excitation frequency close to the sum of the first two natural frequencies. In contrast, the forced motion is usually stable when the excitation amplitude is small enough and the excitation frequency is off the combination resonance, such as the blue point S0 (depicted in Figure 2) where β = 0.001 and σ0 ≈ −0.4047 (ω = 25).
According to Equation (50), Y1 and Y2 decay with time in the stable, off-resonance region Therefore to the first-order approximation, the steady-state voltage output is approximately given by Equation (45).
In addition, the differential harmonic oscillations are also deduced when ω ≈ ω2 − ω1. The multiscale analysis is quite similar to that of the summed harmonic oscillations. Nevertheless, the potential resonance cannot take place when the excitation frequency goes around the difference between the first two natural frequencies. The stable steady-state response possesses the same solution (59).

Numerical Simulations when ω = ω2 + ω1
For the sake of numerical computations, the applied two-term truncation system (34) is rewritten in the form of ( ) ( ) ( ) where pn (n = 1, 2) are introduced and regarded as velocity.
Runge-Kutta method is one powerful numerical technique for solving differential equations [29]. The C programming language has been used to implement the numerical program of Equation (60). In the following, the potential steady-state responses are numerically computed by the use of the classical fourth-order Runge-Kutta method. The dimensionless time step is 2π/(10,000ω) in the C program source code. All the initial conditions are zero as the final behavior of the forced system is independent of them. In contrast, the forced motion is usually stable when the excitation amplitude is small enough and the excitation frequency is off the combination resonance, such as the blue point S 0 (depicted in Figure 2) where β = 0.001 and σ 0 ≈ −0.4047 (ω = 25).
According to Equation (50), Y 1 and Y 2 decay with time in the stable, off-resonance region Therefore to the first-order approximation, the steady-state voltage output is approximately given by Equation (45).
In addition, the differential harmonic oscillations are also deduced when ω ≈ ω 2 − ω 1 . The multiscale analysis is quite similar to that of the summed harmonic oscillations. Nevertheless, the potential resonance cannot take place when the excitation frequency goes around the difference between the first two natural frequencies. The stable steady-state response possesses the same solution (59).

Numerical Simulations When
For the sake of numerical computations, the applied two-term truncation system (34) is rewritten in the form of . q n = p n ; .
where p n (n = 1, 2) are introduced and regarded as velocity.
Runge-Kutta method is one powerful numerical technique for solving differential equations [29]. The C programming language has been used to implement the numerical program of Equation (60). In the following, the potential steady-state responses are numerically computed by the use of the classical fourth-order Runge-Kutta method. The dimensionless time step is 2π/(10,000ω) in the C program source code. All the initial conditions are zero as the final behavior of the forced system is independent of them.  Figure 3c compares the steady-state phase diagram from two different methods. The solid line stands for numerical results (60), and the solid blue dots for analytical approximation (59). These figures support that the steady-state response is periodic, although there is a quantitative difference between analytical prediction and numerical simulation. In Figure 4d, the Poincare map theory is also employed to verify the periodic motion. The multiscale analysis predicts that only a point emerges in the phase diagram. The solid blue dot is the analytical result. The hollowed dots are based on numerical computation. Figure 4d reveals that the steady-state voltage output is a periodic response whose period is the same as the external excitation. Figure 3 discusses the time history of voltage output at Point S0 in Figure 2. The initial history of output voltage is illustrated in Figure 3a where the considered parameters are α = 0.001, β = 0.001, c0 = 0.5084, d = 0.15, h = 0.004167, η = 0.02, c = κ = 0.000001, θ = 0.0002, and ω = 25 (σ0 ≈ −0.4047). The transient motion is short and the steady-state motion occurs in the terminal interval [600, 602]. Figure 3b compares the analytical voltage output (59) with the steady-state response (60) based on the Runge-Kutta method. Figure 3c compares the steady-state phase diagram from two different methods. The solid line stands for numerical results (60), and the solid blue dots for analytical approximation (59). These figures support that the steady-state response is periodic, although there is a quantitative difference between analytical prediction and numerical simulation. In Figure 4d, the Poincare map theory is also employed to verify the periodic motion. The multiscale analysis predicts that only a point emerges in the phase diagram. The solid blue dot is the analytical result. The hollowed dots are based on numerical computation. Figure 4d reveals that the steady-state voltage output is a periodic response whose period is the same as the external excitation.      Figure 4b. The phase diagram in Figure 4c is not a simple closed curve, which supports the conclusion that the forced vibration is   Figure 4b. The phase diagram in Figure 4c is not a simple closed curve, which supports the conclusion that the forced vibration is aperiodic. The Poincare map in Figure 4d illustrates the quasiperiodic response, as a closed curve appears in the phase diagram.

Numerical Simulations When
To discuss the potential combination resonance of difference type, two numerical computations are carried out here. To numerically support the analytical prediction (59) It can be observed that the response remains stable in Figure 6 with a large excitation amplitude and small damping. The external excitation determines the period of the steady-state response. The numerical result quantitatively explains the theoretical prediction (59). It is very evident that the combination resonance of difference type does not take place in the absence of internal resonance.  Besides, the stable voltage output is independent of viscoelastic damping according to Equations (43) and (59). It can be theoretically concluded that the damping has no effect on the stable Besides, the stable voltage output is independent of viscoelastic damping according to Equations (43) and (59). It can be theoretically concluded that the damping has no effect on the stable steady-state response. To some extent, the numerical computations with different damping almost result in an identical outcome.

The First Subharmonic Resonance
In this section, introduce the detuning parameter according to and then Equation (40b) becomes Eliminating those terms that produce secular terms has The solution to Equation (64b) is The solution to Equation (64a) is in the form of where r, A 1r , and A 1i are real numbers. Equation (64a) can be rewritten as For a nontrivial solution, two characteristic values are Here, Y 1 oscillates and decays if However, Y 1 decays without oscillating if Therefore to the first-order approximation, the steady-state voltage output is approximately given by Otherwise the forced motion is unstable. The stability boundary is When the exciting frequency is close to twice the first natural frequency, the stability strongly depends on the viscosity damping. If the damping coefficient is very small, the output voltage might grow without bound according to Equation (66). In addition, when the excitation amplitude is small enough, the forced motion is stable off the first subharmonic resonance, such as the blue point S 1 (depicted in Figure 7) where β = 0.001 and σ 1 ≈ −0.09674 (ω = 8).
Otherwise the forced motion is unstable. The stability boundary is ( ) ( ) When the exciting frequency is close to twice the first natural frequency, the stability strongly depends on the viscosity damping. If the damping coefficient is very small, the output voltage might grow without bound according to Equation (66). In addition, when the excitation amplitude is small enough, the forced motion is stable off the first subharmonic resonance, such as the blue point S1 (depicted in Figure 7) where β = 0.001 and σ1 ≈ −0.09674 (ω = 8).

Numerical Simulations when ω = 2ω1
At the stable point, Point S1 in Figure 7, the numerical simulation (60) is carried out and compared with the analytical steady-state voltage output (71).

Numerical Simulations When ω = 2ω 1
At the stable point, Point S 1 in Figure 7, the numerical simulation (60) is carried out and compared with the analytical steady-state voltage output (71). Figure 8a is the initial history of output voltage with α = 0.001, β = 0.001, c 0 = 0.5084, d = 0.15, h = 0.004167, η = 0.02, c = κ = 0.000001, θ = 0.0002, and ω = 8 (σ 1 ≈ −0.09674). In the terminal interval [600, 604], the steady-state response, the phase diagram, and the Poincare map are respectively compared in Figure 8b-d. The solid black lines stand for numerical results (60), and the solid blue dots for analytical prediction (71). In Figure 8d, the hollowed dots are numerically obtained by Poincare mapping. Two different methods almost yield the same mapping results in Figure 8d. Obviously, the system generates a steady-state voltage in the analytically stable region, and the oscillating frequency is the same as the excitation frequency.  In the unstable region, the excitation amplitude is required to be particularly large.  Figure 9a describes an irregular initial motion, and Figure 9b illustrates a steady-state voltage output, which is also observed from the phase diagram in Figure 9c and the Poincare map in Figure 9d. According to the Poincare map theory, the steady-state response is period-doubling because there are two points in the phase plane.

The Second Subharmonic Resonance
In this case, introduce the detuning parameter according to

The Second Subharmonic Resonance
In this case, introduce the detuning parameter according to and then Equation (40b) becomes Eliminating those terms that produce secular terms gives The solution to Equation (75a) is The solution to Equation (75b) can be expressed in the form of where r, A 2r , and A 2i are real numbers. Equation (75b) can be rewritten as For a nontrivial solution, two characteristic values are The electromechanical system is stable if and only if all roots (79) of the characteristic equation possess a negative real part. Then, the potential steady-state motion is Here, the stability boundary is Near the second subharmonic resonance, Figure 10 plots the stability boundaries (81) for different combinations of external excitation and viscoelastic damping. The dotted, solid, dashed and dotted-dashed lines represent α = 0, α = 0.0001, α = 0.0002, and α = 0.0003, respectively. The considered parameters are c 0 = 0.5084, d = 0.15, h = 0.004167, η = 0.02, c = κ = 0.000001, and θ = 0.0002. It is revealed that the output voltage might grow without bound according to Equation (77). The necessary conditions are theoretically an extremely small damping, a large enough excitation amplitude, and the suitable excitation frequency close to twice the second natural frequency. As the method of multiple time scales requires that the response amplitude is small, the forced motion is usually stable off the subharmonic resonance, such as the blue point S 2 (depicted in Figure 10) where β = 0.001 and σ 2 ≈ 0.2874 (ω = 43).

Numerical Simulations when ω = 2ω2
At the stable point S2 in Figure 10, the numerical simulation is used to verify the analytical steady-state voltage output (80). Figure 11a is the initial history of output voltage with α = 0.001, β = 0.001, c0 = 0.5084, d = 0.15, h = 0.004167, η = 0.02, c = κ = 0.000001, θ = 0.0002, and ω = 43 (σ2 ≈ 0.2874). The period is much smaller than those of previous examples if a periodic response exists. In the terminal interval [600, 601], the steady-state response, the phase diagram and the Poincare map are respectively compared in Figure 11b-d. The solid black lines stand for numerical results (60), and the solid blue dots for analytical prediction (80). In Figure 11d, the hollowed dots are numerically obtained by Poincare mapping, and approach the analytical result. In the stable region for the second subharmonic resonance, the system might generate a steady-state voltage output, and the oscillating frequency is the same as the excitation frequency.
It is numerically convinced that the second subharmonic resonance exists when ω = 2ω2 ≈ 42.7 (σ2 = 0). Figure 12 gives a growing response without bound. During the unstable motion, the considered parameters are α = 0.0001, β = 0.01, c0 = 0.508394, d = 0.15, h = 0.004167, η = 0.02, c = κ = 0.000001, and θ = 0.0002. The initial history and its phase diagram are plotted in Figure 12a,b, respectively. It can be seen that the output voltage is far away from the original point. This conclusion is very useful to capture vibratory energy. However, the large-amplitude swaying

Numerical Simulations When ω = 2ω 2
At the stable point S 2 in Figure 10, the numerical simulation is used to verify the analytical steady-state voltage output (80). Figure 11a is the initial history of output voltage with α = 0.001, β = 0.001, c 0 = 0.5084, d = 0.15, h = 0.004167, η = 0.02, c = κ = 0.000001, θ = 0.0002, and ω = 43 (σ 2 ≈ 0.2874). The period is much smaller than those of previous examples if a periodic response exists. In the terminal interval [600, 601], the steady-state response, the phase diagram and the Poincare map are respectively compared in Figure 11b-d. The solid black lines stand for numerical results (60), and the solid blue dots for analytical prediction (80). In Figure 11d, the hollowed dots are numerically obtained by Poincare mapping, and approach the analytical result. In the stable region for the second subharmonic resonance, the system might generate a steady-state voltage output, and the oscillating frequency is the same as the excitation frequency.

The Experimental Explanation
So far, we have not carried out the resonances experimentally. Nevertheless, Lin et al. tested a nonlinear energy harvester around a subharmonic resonance [11,17]. This kind of secondary resonance does not exist in any other linear system.
In [11], a piezoelectric cantilever with magnetic coupling is used to enhance the width of the frequency band. Their nonlinear system is simplified as a single-degree-of-freedom model. The natural frequency is approximately 10 Hz and the resonant frequency is a bit larger. When the excitation frequency is about 24 Hz, an obvious peak appeared in Figure 2 of [11]. The vibration amplification phenomenon means at least two conclusions: (a) A potential subharmonic resonance exists in the piezoelectric energy harvester; (b) the magnetic force can generate quadratic nonlinearity. However, the subharmonic resonance did not catch the authors' attention at once. It is numerically convinced that the second subharmonic resonance exists when ω = 2ω 2 ≈ 42.7 (σ 2 = 0). Figure 12 gives a growing response without bound. During the unstable motion, the considered parameters are α = 0.0001, β = 0.01, c 0 = 0.508394, d = 0.15, h = 0.004167, η = 0.02, c = κ = 0.000001, and θ = 0.0002. The initial history and its phase diagram are plotted in Figure 12a,b, respectively. It can be seen that the output voltage is far away from the original point. This conclusion is very useful to capture vibratory energy. However, the large-amplitude swaying motion is capable of causing structural damage. A combination of security and performance should be considered in detail for real use.

The Experimental Explanation
So far, we have not carried out the resonances experimentally. Nevertheless, Lin et al. tested a nonlinear energy harvester around a subharmonic resonance [11,17]. This kind of secondary resonance does not exist in any other linear system.
In [11], a piezoelectric cantilever with magnetic coupling is used to enhance the width of the frequency band. Their nonlinear system is simplified as a single-degree-of-freedom model. The natural frequency is approximately 10 Hz and the resonant frequency is a bit larger. When the excitation frequency is about 24 Hz, an obvious peak appeared in Figure 2 of [11]. The vibration

The Experimental Explanation
So far, we have not carried out the resonances experimentally. Nevertheless, Lin et al. tested a nonlinear energy harvester around a subharmonic resonance [11,17]. This kind of secondary resonance does not exist in any other linear system.
In [11], a piezoelectric cantilever with magnetic coupling is used to enhance the width of the frequency band. Their nonlinear system is simplified as a single-degree-of-freedom model. The natural frequency is approximately 10 Hz and the resonant frequency is a bit larger. When the excitation frequency is about 24 Hz, an obvious peak appeared in Figure 2 of [11]. The vibration amplification phenomenon means at least two conclusions: (a) A potential subharmonic resonance exists in the piezoelectric energy harvester; (b) the magnetic force can generate quadratic nonlinearity. However, the subharmonic resonance did not catch the authors' attention at once. They only discussed the output voltage when the excitation frequency is 6.5 Hz, 10 Hz, and 13 Hz. Further study is carried out in [17].
An improved experiment is set up by adjusting the distance between two magnets [17]. More time histories of the output voltage are observed and compared. At the excitation frequency of 20 Hz, the output voltage is periodic. The nonlinear cantilever shows the response amplitude five times as large as that of a linear cantilever without magnetic coupling. When the excitation frequency is between 23 Hz and 30 Hz, the amplitude of output voltage is almost the same as the linear cantilevered system. These experiments have indicated that magnetic nonlinearity is capable of improving the performance of energy harvesting. Compared with a linear harvester, the nonlinear system can indeed provide more electrical energy because of the subharmonic resonances or other nonlinear characteristics.
As theoretical compensation, we intend to apply our theory to the experiment in [17]. Although the steady-state response of voltage is periodic, its period is not that of the external excitation any longer. Instead, the period of the voltage output is exactly twice the period of the excitation at any subharmonic resonance. The period-doubling motion can also easily be identified from Figure 3 of [17].
Moreover, not all the voltage output is period-doubling. According to our theoretical results, the steady-state response might be a periodic motion of single frequency. The more complex response might occur at the combination resonance of the summed type. In the absence of internal resonances, there is no amplitude amplification phenomenon when the excitation frequency approaches the difference between two natural frequencies. These theoretical analyses not only forecast the stability boundary, but also demonstrate different types of the output voltage. Obviously, the stability boundary is quite difficult to be analyzed by experiments. The theoretical significance is that it can be used to more reasonably design energy harvesters.

The Effects of Material Damping
Piezoelectric ceramics are ferroelectric synthetic materials. In general, it can be manufactured from Pb (lead), Zr (zirconium), and Ti (titanium). This compound class shows higher piezo-mechanical efficiency than quartz. Meanwhile, a different combination of Pb, Zr, and Ti is capable of producing diversity in the damping of piezoelectric materials. Table 2 displays the different mechanical quality factors and dielectric losses for the piezo-materials manufactured by APC International, Ltd., formerly American Piezo Ceramics, Inc., Mackeyville, PA, USA [33].
In practice, piezoelectric materials also show changing parameter values owing to their thickness, actual shape, surface finish, shaping process, and postprocessing. In general, a hard piezo-material, which is harder to depolarize than a soft piezo-material, exhibits higher mechanical quality factors and lower dielectric losses. This combination of piezoelectric and mechanical properties makes hard piezo-materials ideal for piezoelectric energy harvesting. In order to design an energy harvester working at the subharmonic or combination resonance, APC 841, APC 844, APC 880, and APC 881 are preferable. Here is what we have to say. The conclusions in [17] are mainly based on an experimental test. The saying, "the stochastic, subharmonic, and ultraharmonic responses produce an average of threefold to fivefold increase in voltage production", is not supported by our in-depth analyses. The comparison should be on the basis of the stability theory proposed in this article. For example, the nonlinear energy harvester has no superiority if the nonlinearity is weak enough, the excitation amplitude is small enough, or the external excitation is off-resonance.

The Effects of Excitation Amplitude
The excitation amplitude is quite important in the stability analyses. The stability of the output voltage is directly affected by the amplitude of the external excitation. As is described in Section 3.1, these secondary resonances occur when the excitation amplitude is large enough.
In the stable region, the output voltage is periodic and the period of the steady-state response is exactly the same as that of the external excitation. From Equation (43), the response amplitude of output voltage is almost proportional to the amplitude of the external excitation. Within the proportional limit, the response amplitude linearly grows with the increased excitation amplitude. When the excitation amplitude is large enough, a secondary resonance might occur.
In the unstable region, an increased excitation amplitude can cause a great rise in the width of the frequency band. This phenomenon is just desirable, as depicted in Figure 2, Figure 7, and Figure 10. As we all know, the linear energy harvester only works at a single frequency. According to the stability analyses, the nonlinear energy harvester can greatly improve the range of working frequency. The enhanced performance enables the nonlinear device to work in a random environment.
Besides, we disagree that the output voltage has the probability of growing without bound, even though the theoretical analysis forecasts the growing response and the numerical simulation follows. The conception is based on the following two reasons at least: (a) The response amplitude of the piezoelectric cantilever must be less than its length; (b) The methodology requires a weak response, which is unsuitable for the large-amplitude motion.

The Effects of Maximal Voltage
In addition, there is another reason why the output voltage must be limited. Piezoelectric ceramic materials have the ability to both generate and respond to voltages. The output voltage controls the applied electric field. An ultra-high electric field can change the piezoelectric properties of piezoelectric ceramics after polarization. Accordingly, the piezoelectric energy harvester should be operated in safe working environments.
The alternating current field is illustrated in this article. The applied maximal voltage (per unit thickness) of some APC materials is provided for reference by APC International, Ltd. 5-7 VAC/mil for APC 850, APC 851, APC 854, APC 855 where one inch = 1000 mil. 9-11 VAC/mil for APC 840, APC 841, APC 842, APC 844, APC 880, APC 881 [33]. From a physical point of view, the mechanical properties of piezoelectric materials are also linked to the time duration of the applied electric field. Ultrahigh electric fields are often easier to cause changes in system parameters, enhanced nonlinear responses, material fatigue damages, and even structural fracture behaviors.
All in all, the theoretical analyses provide a novel solution to enhance electrical performance as well as a wider range of working frequencies to some extent. However, the output voltage or excitation amplitude must be less than its maximum permitted. The efficiency to harvest mechanical energy is worthy of being experimentally studying. A combination of security and performance should be considered in detail for real use.

Conclusions
For the first time, the summed combination resonance is developed to harvest vibratory energy. A cantilever-type piezoelectric energy harvester with tip magnets is investigated. The mathematical model is an integro-partial differential equation with time-dependent boundary conditions. A newly-built methodology is proposed to deal with the issue of continuum mechanics. Several dynamic phenomena are found to be of interest and novelty, such as the quasiperiodic motion and the theoretically growing response without bound.
The subharmonic resonances are also approximately analyzed by the improved Galerkin method. Based on the two-term truncation system with quadratic nonlinearity, the multiscale method is developed to solve these secondary resonances. The numerical simulation obtained by Runge-Kutta method well matches those stable results obtained analytically. The stability theory is quite useful to design a piezoelectric energy harvester. To some extent, these secondary resonances may provide both a much larger voltage output and a wider bandwidth of working frequency, as partially experimented in [11,17].
The differential combination resonance is discussed analytically and numerically, too. The main conclusions are as follows.

•
When the excitation frequency approaches the sum of two natural frequencies, the combination resonance of the summed type exists. The stability boundary is provided approximately but analytically and sensitive to the viscoelastic damping. In the stable region, the one-order approximation of steady-state response is formulated. In the unstable region, the quasiperiodic motion may occur. The unstable response is readily found with small damping and a large enough excitation amplitude.

•
When the excitation frequency approaches the difference between two natural frequencies, the differential harmonic oscillation does not exist. Their long-term response is simply harmonic and their frequency output is the same as the excitation frequency. To some extent, viscoelastic damping has no effect on the steady-state response. This outcome is drawn in the absence of internal resonance.

•
When the excitation frequency approaches twice the natural frequency, the subharmonic resonance exists. The stability boundaries are provided for the first two vibration modes. Increasing damping decreases the unstable region, where there might be a period-doubling motion, or a theoretically growing response without bound. On the contrary, simply-harmonic voltage output is produced off the subharmonic resonances. The steady-state response in the stable region is explicitly expressed and numerically verified.
In summary, the summed combination resonances and the subharmonic resonances can be exploited to harvest more mechanical energy. Different from the primary resonance, these secondary resonances might improve the bandwidth towards two different directions. In general, the primary resonance merely enhances the bandwidth towards one direction, which is well known as the soft-spring or hard-spring behavior, respectively. The proposed stability analyses can be used to design the piezoelectric energy harvester. The stability boundary is very sensitive to material damping. In the stable region, the nonlinear energy harvester cannot produce more voltage output than the linear device, which is independent of damping. On the contrary, the secondary resonance energy harvester must work in an unstable region. With the help of further experiments, the low-cost piezoelectric energy harvester can support a wide range of industries, including structural health monitoring sensors and wireless networks. To this end, APC 841, APC 844, APC 880, and APC 881 are preferable piezoelectric ceramic materials to manufacture these electronic applications.