Modeling and Evaluation of a Multi-Stable Hybrid Energy Harvester

Abstract

This article develops a multi-stable hybrid energy harvester (MSHEH) which consists of a piezoelectric energy harvester (PEH) and an electromagnetic energy harvester (EMEH). By tuning two parameters, the MSHEH can achieve a mono-stable, bi-stable, and tri-stable state, respectively. A numerical procedure is developed to compute the EMEH’s transduction factor. The obtained result is validated experimentally. Using the equivalent magnetic 2-point dipole theory, the restoring force model of the magnetic spring is established. The obtained model is verified experimentally. The energy harvesting performances of the MSHEH under the four different configurations (linear, mono-stable, bi-stable and tri-stable) subjected to frequency sweep excitations are evaluated by simulation and validated by experiment. The comparative analysis focuses on power output, accumulated harvested energy, and effective energy-harvesting bandwidth. The optimum load resistances are investigated by Pareto front optimizations. The following key findings are obtained. When subjected to high-level frequency sweep excitation, the tri-stable configuration exhibits the widest frequency bandwidth and the highest total accumulated harvested energy. When subjected to low-level frequency sweep excitation, the bi-stable configuration is more efficient in energy harvesting. The best performance trade-off between the PEH and EMEH can be achieved by selecting the optimum load resistances properly.

1. Introduction

A vibration energy harvester (VEH) is a device that converts ambient mechanical energy into electrical energy. There are various ambient mechanical energies that can be captured, such as structural vibration [1], machinery vibration [2], and human motion [3]. The VEH provides a promising solution to a growing demand for self-sustainable power supply for wearable electronic devices and wireless sensor node networks, especially when deploying conventional power sources such as power lines or batteries is inconvenient or impractical [4].

In general, a traditional VEH consists of a linear oscillator that has a narrow operation frequency bandwidth. Over the last two decades, there has been a growing interest in enhancing the working bandwidth and energy harvesting efficiency of VEHs for different environments. Introducing nonlinearity is one of the promising solutions to broaden the working bandwidth of VEHs. Various nonlinear VEHs have been proposed [5]. According to the system stability state, the nonlinear VEHs can be classified as mono-stable and multi-stable, such as bi-stable or tri-stable. A mono-stable energy harvester reported in [6] consists of a piezoelectric cantilever beam with a tip magnet subjected to an external magnetic field generated by a pair of fixed magnets. Such a mono-stable system can exhibit softening or hardening behaviors when the magnetic interaction is adjusted. The energy harvesting performance of a mono-stable energy harvester was investigated in [7]. The study showed that the high-branch oscillation leads to a high energy harvesting efficiency. A bi-stable energy harvester can be used to improve energy harvesting performance by utilizing the snapping-through feature. As proved in [8], the inter-well oscillation of a bi-stable energy harvester can significantly enhance its power output performance. The study reported in [9] showed that a bi-stable energy harvester with an elastic magnifier can provide higher power output and wider working bandwidth. One of the drawbacks of bi-stable energy harvesters is the requirement of a sufficient excitation level in order to overcome the barrier of the potential wells. The tri-stable energy harvester was proposed to address this drawback. Based on the configuration of the bi-stable energy harvester in [8], a tri-stable energy harvester was achieved by tuning the orientation43s [10] and the positions [11,12] of the two fixed magnets. Moreover, the performance of an asymmetric tri-stable energy harvester was investigated in [12]. The studies showed that the proposed tri-stable energy harvester outperforms the bi-stable energy harvester under the low-level excitation in terms of the voltage output.

On the other hand, the concept of hybrid systems has been proposed to enhance both the energy density and the power output. There are two kinds of hybrid systems: the first one can convert multiple energy sources such as solar, thermal or vibration into electricity [13]; the second one can convert a single energy source such as vibration into electricity through multiple conversion mechanisms [14]. This study considers the second type. There are three main transduction mechanisms for the VEH, namely, piezoelectric [15], electromagnetic [16], and electrostatic [17]. Each of them has its own advantages and disadvantages. For example, the piezoelectric energy harvester (PEH) has high energy density and is easy to deploy, the electromagnetic energy harvester (EMEH) shows the benefits of high current output and ease of maintenance, and the electrostatic energy harvester has the advantages of compact design and wider working bandwidth. A VEH combined with two or more transduction mechanisms is referred to as hybrid energy harvester (HEH), yielding better efficiency and robustness [18,19]. This paper focuses on the HEH consisting of a PEH and an EMEH.

An HEH proposed in [20] consists of a cantilever beam patched with a PEH and attached with a tip magnet that moves inside a coil placed on the base. The study provided an approach of coupling the PEH and EMEH to increase the power output. A power management circuit was designed in [21] to overcome the impedance mismatching issue of the HEH. The HEHs proposed in [22] and [23] utilized a 2-degree-freedom structure to improve the power output. To enhance the performance of the HEH under ultra-low frequency excitation, the frequency up-conversion design of the HEH was proposed in [24]. In addition, a multi-modal HEH was developed in [25] to make the system able to operate at four different resonant modes, significantly widening the operation bandwidth. There have been conflicting views on the benefits of a linear HEH under harmonic excitation. For example, a recent study [26] showed that under harmonic excitation, an idealized two-port linear HEH with the electrical loss neglected offers little benefit in terms of the maximum output power. On the other hand, introducing nonlinearity to the HEH has been explored by some researchers. For example, a mono-stable HEHs proposed in [27,28] showed that the nonlinearity can significantly boost the energy density and widen the frequency bandwidth. A bi-stable HEH was proposed in [29] to improve the power output. In the study, an approximate method was used to simplify the modeling of the coupled system. A bi-stable HEH developed in [30] used the tunable stiffness design to achieve better adaptability for various environments. Further, the studies reported in [31,32] showed that the tri-stable HEH is beneficial for enhancing both operation bandwidth and output power compared with the mono-stable and bi-stable HEHs.

The above review indicates a need for a nonlinear HEH whose stability states can be adjusted in order to achieve better adaptability in terms of power output and frequency bandwidth. To address such a need, based on our previous study [33], a tunable multi-stable hybrid energy harvester (MSHEH) is proposed in this study. Differently from the existing designs that use two external magnets [34,35,36,37], the MSHEH employs a single external magnet, which makes the magnetic spring more compact and makes implementation of an EMEH easy. The EMEH is realized by placing one set of six coils above and one set of six coils below the two moving magnets. With this novel arrangement, the magnetic flux on both the moving magnets’ upper surface and the lower surface can be effectively utilized, and the space efficiency of the EMEH can be improved compared with the existing designs, such as the ones in [29,38,39].

The contributions of the present paper lie in four aspects. Firstly, the proposed MSHEH is novel in terms of stability tuning and the EMEH design. Secondly, a numerical modeling procedure is developed to determine the transduction factor of the EMEH. Thirdly, a comparative study is conducted to evaluate the energy harvesting performances of four different configurations subjected to the frequency sweep excitation. Fourthly, a Pareto front optimization is conducted to maximize the power output of both EMEH and PEH under harmonic excitations with various exciting frequencies. In addition, further optimization is conducted to maximize the accumulated harvested energy for both EMEH and PEH under high-level frequency up-sweep excitation.

2. Apparatus and Modeling

Figure 1a shows a solid modeling drawing of the proposed MSHEH. As shown in the figure, a thin stainless-steel beam is clamped to a platform which is fastened to a base by using four aluminum extrusions. Each side of the upper end of the beam is attached by a piezoelectric transducer or PZT (S128-J1FR-1808YB, Midé), while its lower end is fixed with a small cylindrical magnet B and attached with a holder for an assembly of two identical cylindrical magnets, A and C. The holder can be fixed on any position along the beam by sliding. A large cylindrical magnet D is fixed in a holder that can slide vertically in a stand on the base. When the cantilever beam is at its equilibrium position or undeflected, the four magnets are situated on the same vertical plane, and magnets B and D are collinear. By sliding the holder for magnet D, the distance between magnet B and magnet D can be adjusted. By sliding the holder for magnets A and C along the beam, the distance between magnets A/C and magnet B, and the distance between magnets A/C and magnet D can be adjusted.

To add an EMEH to the system, 12 coils are placed symmetrically between magnets A and C, i.e., 6 coils above and 6 coils below. Each of the coils is held in a holder that allows individual adjustment of the coil’s position and orientation. Through adjustment, the end surfaces of the coils are approximately parallel to the oscillation trajectory of magnets A and C. Figure 1b illustrates the spatial positions of the coils: those on the side of magnet C are labelled as 1 to 6 while those on the side of magnet A are labelled as 1′ to 6′. Figure 1c shows the polarities of the four magnets where ${\mathit{m}}_A$, ${\mathit{m}}_B$, ${\mathit{m}}_C$, ${\mathit{m}}_D$ are the magnetic moment vectors, and $A_0$, $B_0$, $C_0$ and A, B, C represent the center positions of magnets A, B, and C when the beam is in undeformed and deformed states, respectively. Note that the origin of the coordinate system is fixed at $B_0$.

Figure 2a,b show the front view and side view of Figure 1a, respectively, where d is the distance between magnet B and magnet D when the beam is undeformed, $h$ is the distance between magnet B and magnets A, C, l is the length of the cantilever beam, and $d_t$ is the distance between the axis of magnet B and that of magnets A and C. As shown in Figure 2b, $x$ and $z$ represent the transverse and longitudinal displacements of the center of magnet B relative to $B_0$, respectively. $\alpha$ is the angle between ${\mathit{m}}_B$ and ${\mathit{m}}_D.$ Since the slope of the beam’s tip is relatively small, it is assumed that $\angle {BOB}_0\approx \alpha$.

Figure 3 shows a lumped parameter model for the simplified system. In the figure, $m=0.09 \mathrm{Kg}$ represents the equivalent mass at the tip of the beam, $w_b$ and $x$ are the displacement of the base and the equivalent mass relative to the base, respectively, $c_m$ = 0.0058 N/m is the mechanical damping coefficient of the system, and $k_n$ is the nonlinear stiffness including the effects of the cantilever beam and the magnetic interaction. The PEH’s circuit is given on the right side of the figure, where $\theta =8.515\times {10}^{-3}$ N/V is the electromechanical coefficient of the PEH, which is identified by the experimental method proposed in [40], and $R_{lp}$ is the resistance of a load resistor connected to the output of the PEH. The EMEH’s circuit is given on the left side of the figure, where $K_t$ is the total transduction factor of the EMEH, $v_{em}$ is the inductive voltage or so-called electromotive force (EMF) of the EMEH, $R_c$ and $L_c$ are resistance and inductance of one coil, respectively, and $R_{le}$ is the resistance of a load resistor connected to the output of the EMEH. Note that as the 12 coils are connected in series, their total resistance and inductance are 12$R_c$ and 12$L_c$, respectively.

Based on Newton’s second law and Kirchhoff’s current law, the governing equations of the system can be derived as follows:

$$m\ddot{x}+c_m\dot{x}+f_e+f_n+\theta v_p=-m{\ddot{w}}_b$$

(1)

$$c_p{\dot{v}}_p+\frac{v_p}{R_{lp}}-\theta \dot{x}=0$$

(2)

where $v_p$ is the voltage over the load resistor of the PEH, $c_p=50\times {10}^{-9}$ F is the capacitance of the PEH, $f_n$ is the total restoring force, and $f_e$ is the electromagnetic force caused by the changes in the magnetic flux through the coils. Based on Lenz’s law, the electromagnetic force can be expressed as follows:

$$f_e=K_{t}I$$

(3)

where $K_t=2{\sum }_{i=1}^6{K}_{ti}$ is the total transduction factor with $K_{ti}$ as the transduction factor for the ith coil and $I$ is the current in the EMEH’s circuit. Note that the values of the transduction factors of coils 1 to 6 and coils 1′ to 6′ are equal, since they have identical configurations at the upper and lower sides of magnet C and magnet A, respectively. Applying Kirchoff’s law to the circuit of the EMEH yields:

$$12L_c\frac{dI}{dt}+\left(12R_c+R_{le}\right)I=v_{em}$$

(4)

By using a multimeter, it is found that $R_c=0.9$ $\mathsf{\Omega }$. By using an inductance meter, it is found that $L_c=0.454$ mH. Since the frequency of vibration considered in this study does not exceed 20 Hz, the inductive impedance of the coil is negligible compared with $R_c$. Thus, the current can be written in the following form:

$$I=\frac{v_{em}}{12R_{coil}+R_{le}}$$

(5)

3. Determination of the EMEH’s Transduction Factor

Due to the unique design of the EMEH, the determination of its total transduction factor is not straightforward. In what follows, a numerical method is employed for this purpose. According to Faraday’s law, the EMF of the EMEH can be expressed as:

$$v_{em}=-2{\displaystyle \sum _{i=1}^6\frac{d{\Phi }_i}{dt}}=-2{\displaystyle \sum _{i=1}^6\frac{\partial {\Phi }_i}{\partial x}}\dot{x}$$

(6)

where ${\Phi }_i$ is the total magnetic flux through the ith coil. In fact, the magnetic flux is not evenly distributed throughout the whole coil due to the complex orientation of the magnets. Thus, each coil is sliced into n layers and the magnetic flux in the jth layer is assumed to be uniformly distributed and denoted as ${\varphi }_{ij}$.

As shown in Figure 4, the layer closest to magnet A or C is labelled as layer 1, which means the bottom layer for the upper side coil and the top layer for the lower side coil are the first layer. Thus, the total magnetic flux in the ith coil can be expressed as:

$${\Phi }_i=\frac{N}{n}{\displaystyle \sum _{j=1}^n{\varphi }_{ij}}$$

(7)

where N is the turns of the coil. The magnetic flux in the jth layer for the ith coil is given by

$${\varphi }_{ij}={\displaystyle \underset{A}{\iint }\left(B_{zij}\cos \beta +B_{xij}\sin \left|\beta \right|\right)}dA$$

(8)

where $B_{xij}$ and $B_{zij}$ is the magnetic flux density in the x and z direction at that layer, respectively, β is the angle of the coil from the horizontal, and A is the area of the end surface of the coil. Substituting Equation (7) into Equation (6) gives:

$$v_{em}=-2\frac{N}{n}{\displaystyle \sum _{i=1}^6{\displaystyle \sum _{j=1}^n\frac{\partial {\varphi }_{ij}}{\partial x}}}\dot{x}=K_t\dot{x}$$

(9)

According to Equation (9), the transduction factor $K_{ti}$ for each coil can be defined as follows:

$$K_{ti}=-\frac{N}{n}{\displaystyle \sum _{j=1}^n\frac{\partial {\varphi }_{ij}}{\partial x}}$$

(10)

which is a function of the change rate of the magnetic flux with respective to the displacement x. In this study, a finite element analysis software Comsol Multiphysics is utilized to compute the change rates of the magnetic flux of the six coils 1 to 6 when magnets C and B are oscillating through them. For the sake of simplicity, the influence of magnet A on coils 1 to 6 is ignored. The geometry of the model built in Comsol is shown in Figure 5a. It should be noted that each of the coils is modelled as n disks to represent the n layers and meshed individually. As shown in Figure 5b, $D_{coil}$ and $h_{coil}$ are the diameter and height of the coil, respectively, $d_g$ is the air gap between the end surfaces of magnet C and the coils 2, 5, and $d_s$ is the lateral distance between the bottom center of the coils 4 and 6 and the center of magnet D. All the values of the parameters of the coils and magnets used in the simulation are listed in

Table 1. Parameters of the coils.

Symbol	Name	Coils 1, 4	Coils 2, 5	Coils 3, 6
$D_{coil}$ (m)	Diameter of the coil	0.029	0.029	0.029
$h_{coil}$ (m)	Height of the coil	0.015	0.015	0.015
$\beta$ (degree)	Angle of the coil from the horizontal	−10	0	10
N	Turn number	245	245	245
$d_s$ (m)	Distance from the center of magnet D	−0.0353	0	0.0353
$d_g$ (m)	Air gap between magnet C and coils 2, 5	0.002	0.002	0.002

and

Table 2. Parameters of the magnets.

Symbol	Name	Magnet A & C	Magnet B	Magnet D
$D_{mag}$ (m)	Diameter of magnet	2.54 × 10−2	7.94 × 10−3	2.54 × 10−2
$l_A, l_C, l_B$, $l_D$ (m)	Hight of magnets	9.525 × 10−3	7.94 × 10−3	2.54 × 10−2
$B_r$ (T)	Residual flux density	1.44	1.28	1.28
	Material	N50	N42	N42

, respectively.

In the simulation, the number n of layers for each coil is set to 12, and magnets B and C oscillate from $x=-0.065$ m to $x=0.065$ m. In order to simulate the trajectory of magnets B and C, the displacement of the center of magnets B and C in the z-axis is modelled as $z\left(x\right)$ which can be derived from the trigonometric relationship in the triangle $O{R}_0{B}_0$ in Figure 2b, as follows:

$$z(x)=l-\sqrt{l^2-x^2},$$

(11)

and α is the angle of magnets B and C from the horizontal and is approximated as ${\sin }^{-1}\left(x/l\right)$, based on the triangle $O{R}_0{B}_0$.

Figure 6 shows the magnetic flux distributions when magnets B and C move from the farthest left position to the farthest right position. In particular, Figure 6a–c illustrate the situations when magnet C is concentric with coils 1, 2 and 3, respectively. During the simulation, the change rate of the magnetic flux through each layer of the coils with the different displacements is recorded.

Based on Equation (10), the transduction factors for all six coils can be computed, and the results are shown as solid lines in Figure 7a,b. Then, these results are curve-fitted using piecewise functions, which are the sum of three sine functions in the specific displacement ranges and can be defined as follows:

$$K_{ti}(x)=\left\{\begin{array}{ll}{\displaystyle \sum _{j=1}^3}{(-1)}^i{a}_{uj}\sin (b_{uj}(x+c_i))& d_i^{\min }\le x\le d_i^{\max }(i=1,2,3)\\ 0& x<d_i^{\min } \& x>d_i^{\max } (i=1,2,3)\end{array}\right.$$

(12)

$$K_{ti}(x)=\left\{\begin{array}{ll}{\displaystyle \sum _{j=1}^3}{(-1)}^{i-1}{a}_{lj}\sin (b_{lj}(x+c_i))& d_i^{\min }\le x\le d_i^{\max }(i=4,5,6)\\ 0& x<d_i^{\min } \& x>d_i^{\max } (i=4,5,6)\end{array}\right.$$

(13)

where $a_{uj}$ and $b_{uj}$ are the curve-fitting constants for coils 1, 2, 3 and $a_{lj}$, $b_{lj}$ are the curve-fitting constants for coils 4, 5, 6; $d_i^{\min }$ and $d_i^{\max }$ are the equation limits for the ith coil, and $c_i$ is the coordinate translation for the ith coil. The equation limits are $d_1^{\min }=-0.0706$ m, $d_4^{\min }=-0.0690$ m, $d_1^{\max }=d_4^{\max }=0$, $d_2^{\min }=-0.0353$ m, $d_5^{\min }=-0.0345$ m, $d_2^{\max }=0.0353$ m, $d_5^{\max }=0.0345$ m, $d_3^{\min }=d_6^{\min }=0$, $d_3^{\max }=0.0706$ m, and $d_6^{\max }=0.0690$ m. The coordinate translations are $c_1=0.0353 \mathrm{m}, c_4=0.0345$ m, $c_2=c_5=0$, $c_3=-0.0353 \mathrm{m}, \mathrm{and} c_6=-0.0345$ m, and

Table 3. Curve-fitting constants of Equations (12) and (13).

Using the original turn number N	Constants	${\mathit{a}}_{\mathit{u}\mathbf{1}}$	${\mathit{b}}_{\mathit{u}\mathbf{1}}$	${\mathit{a}}_{\mathit{u}\mathbf{2}}$	${\mathit{b}}_{\mathit{u}\mathbf{2}}$	${\mathit{a}}_{\mathit{u}\mathbf{3}}$	${\mathit{b}}_{\mathit{u}\mathbf{3}}$
	Value	−0.94	88.80	−0.29	179.47	−0.06	344.84
	Constants	${\mathit{a}}_{\mathit{l}\mathbf{1}}$	${\mathit{b}}_{\mathit{l}\mathbf{1}}$	${\mathit{a}}_{\mathit{l}\mathbf{2}}$	${\mathit{b}}_{\mathit{l}\mathbf{2}}$	${\mathit{a}}_{\mathit{l}\mathbf{3}}$	${\mathit{b}}_{\mathit{l}\mathbf{3}}$
	Value	−1.06	99.08	−0.28	207.30	−0.084	376.82
Using equivalent turn number Ne	Constants	${\mathit{a}}_{\mathit{u}\mathbf{1}}$	${\mathit{b}}_{\mathit{u}\mathbf{1}}$	${\mathit{a}}_{\mathit{u}\mathbf{2}}$	${\mathit{b}}_{\mathit{u}\mathbf{2}}$	${\mathit{a}}_{\mathit{u}\mathbf{3}}$	${\mathit{b}}_{\mathit{u}\mathbf{3}}$
	Value	−0.68	87.93	−0.22	177.48	−0.05	343.65
	Constants	${\mathit{a}}_{\mathit{l}\mathbf{1}}$	${\mathit{b}}_{\mathit{l}\mathbf{1}}$	${\mathit{a}}_{\mathit{l}\mathbf{2}}$	${\mathit{b}}_{\mathit{l}\mathbf{2}}$	${\mathit{a}}_{\mathit{l}\mathbf{3}}$	${\mathit{b}}_{\mathit{l}\mathbf{3}}$
	Value	−0.78	99.08	−0.20	207.30	−0.06	376.81

lists the obtained curve-fitting constants. The curve-fitting results are shown in dotted lines in Figure 7a,b. Overall, the curve-fitting results agree well with the numerical results. Figure 7c shows the total transduction factor, which is a strongly nonlinear function of $x$, reaching the maximum values around $\left|x\right|=0.017$ m and $\left|x\right|=0.045$ m, respectively.

According to the verification method proposed in [41], an experiment is carried out to verify the computed transduction factors for the six coils. Since the functions of the transduction factor of coils 3 and 6 are symmetric with coils 1 and 4, in this case only the transduction factors of coils 1, 2, 4 and 5 need to be measured. The range of the initial displacement position x has been chosen from $-0.04 \mathrm{m}$ to 0 m with an interval of 0.005 m. The measurement results are shown in Figure 8. It can be seen that the experimental results (blue stars) are slightly lower than the original simulation results (red lines). The deviation arises from the assumption that all the turns are located at the outer shell of the coil. The assumption lead to an overestimation of the coil turn number $N$ in Equation (10). To enhance the accuracy of the transduction factors, the equivalent turn number $N_e$ for the coils need to be estimated. Based on the obtained experimental results, an approximate equivalent turns number is found to be 180 through trial and error. Based on the simulation results for coils 2 and 5, the constants of the curve fitting functions shown in Equations (12) and (13) can be obtained, and are and listed in Table 3. As shown in Figure 8, the measured data, the simulation results based on the equivalent turn number $N_e$ (green lines) and the curve fitting functions (black lines) match well.

4. Determination of the Nonlinear Restoring Force

The total restoring force $f_n$ of the system in the x-direction consists of an equivalent force $f_g$ due to the gravity, a restoring force $f_b$ due to the beam’s elasticity, an attractive magnetic force $f_{DBx}$ between magnet D and magnet B and two repulsive magnetic forces: $f_{DAx}$ between magnet D and magnet A, and $f_{DCx}$ between magnet D and magnet C. Since magnets A and C are identical and symmetrical about the central line of the beam, the values of $f_{DAx}$ and $f_{DCx}$ are equal. Then, the total restoring force can be expressed as:

$$f_n=f_g+f_b+f_{DBx}+f_{DAx}+f_{DCx}=mg\cdot \tan \left(\alpha \right)+k_{b}x+f_{DBx}+2f_{DAx}$$

(14)

where $k_b$ = 90.1 N/m is the stiffness of the beam, which can be determined experimentally. In what follows, the analytical restoring forces $f_{DBx}$ and $f_{DAx}$ will be found using the equivalent magnetic 2-point dipole model proposed in [42]. To have a better understanding of the magnetic force model, Figure 9a,b show the front view of the apparatus when the beam is undeformed and deformed.

As shown in Figure 9a,b, ${\mathit{r}}_{15}, {\mathit{r}}_{25}, {\mathit{r}}_{35}$ and ${\mathit{r}}_{45}$ are the vectors from $Q_5$ to $Q_1, Q_2, Q_3$ and $Q_4$, respectively, and ${\mathit{r}}_{16}, {\mathit{r}}_{26}, {\mathit{r}}_{36}$ and ${\mathit{r}}_{46}$ are the vectors from $Q_6$ to $Q_1, Q_2, Q_3$ and $Q_4$ respectively, where $Q_i$, $i=1,\dots ,6$ are the total surface charges of the magnets defined by:

$$\begin{array}{l}{Q}_1=-M{S}_B Q_2=M{S}_B Q_3=M{S}_A Q_4=-M{S}_A \\ Q_5=-M{S}_D Q_6=M{S}_D \end{array}$$

(15)

where $S_A=5.07\times {10}^{-4} {\mathrm{m}}^2$, $S_B=4.95\times {10}^{-5} {\mathrm{m}}^2$ and $S_D=5.07\times {10}^{-4} {\mathrm{m}}^2$ are the surface area of magnets B, A and D, respectively, $M=B_r/\mu$ is the magnetization of magnets A, B and D, where $B_r$ is the magnetic residual flux density; their values are listed in Table 2, and $\mu =4\pi \times {10}^{-7}\mathrm{H}/\mathrm{m}$ is the vacuum permeability.

The magnetic force between magnet B and magnet D is considered first. Based on the Boit–Savart law, the magnetic force exerted by magnet B on magnet D is the combination of the magnetic force exerted from $Q_1$ and $Q_2$ to $Q_5$ and $Q_6$, which is given in the following equation:

$${\mathit{f}}_{DB}=Q_1\frac{{\mu }_0}{4\pi }\left(Q_5\frac{{\mathit{r}}_{15}}{{\left|{\mathit{r}}_{15}\right|}^3}+Q_6\frac{{\mathit{r}}_{16}}{{\left|{\mathit{r}}_{16}\right|}^3}\right)+Q_2\frac{{\mu }_0}{4\pi }\left(Q_5\frac{{\mathit{r}}_{25}}{{\left|{\mathit{r}}_{25}\right|}^3}+Q_6\frac{{\mathit{r}}_{26}}{{\left|{\mathit{r}}_{26}\right|}^3}\right)$$

(16)

where ${\mathit{r}}_{15}$, ${\mathit{r}}_{16}$, ${\mathit{r}}_{25}$ and ${\mathit{r}}_{26}$ can be derived from the position vectors of $Q_1$, $Q_2$, $Q_5$ and $Q_6$, respectively. According to Equation (14), to obtain the total restoring force, only the $f_{DBx}$ is considered, which can be expressed as follows [33]:

$$f_{DBx}=-\frac{{\mu }_0}{4\pi }\left[Q_1\left(x-\frac{l_B}{2}\sin \alpha \right)\left(\frac{Q_5}{{\gamma }_1}+\frac{Q_6}{{\gamma }_2}\right)+Q_2\left(x+\frac{l_B}{2}\sin \alpha \right)\left(\frac{Q_5}{{\gamma }_3}+\frac{Q_6}{{\gamma }_3}\right)\right]$$

(17)

where the expressions ${\gamma }_1$, ${\gamma }_2$, ${\gamma }_3$ and ${\gamma }_4$ are given in [33]. Further, the magnetic force between magnet A and D in the x-direction can also be obtained as:

$$f_{DAx}=-\frac{{\mu }_0}{4\pi }\left[Q_3\left(x-h\sin \alpha -\frac{l_A}{2}\right)\left(\frac{Q_5}{{\gamma }_5}+\frac{Q_6}{{\gamma }_6}\right)+Q_4\left(x-h\sin \alpha +\frac{l_A}{2}\right)\left(\frac{Q_5}{{\gamma }_7}+\frac{Q_6}{{\gamma }_8}\right)\right]$$

(18)

where ${\gamma }_5$, ${\gamma }_6$, ${\gamma }_7$ and ${\gamma }_8$ are also defined in [33]. By substituting Equations (17) and (18) into Equation (14), the total restoring force can be obtained.

To validate the model, the five different configurations are considered: Case (I) d = 0.0605 m, h = 0.0035 m; Case (II) d = 0.0496 m, h = 0.0058 m; Case (III) d = 0.0452 m, h = 0.0058 m; Case (IV) d = 0.0339 m, h = 0.0092 m; and Case (Ⅴ) d = 0.0330 m, h = 0.0079 m. Among them, the first case is the mono-stable configuration, and the second and third cases are the bi-stable configurations. By applying the original values of the total charges that are listed in the first column of

Table 4. Values of the total charges on the surfaces of different magnets.

		Original Values	Optimum Values (6 Parameters)	Optimum Values (5 Parameters)
Magnet A	$Q_3$	$580.64$	$1150$	$642.36$
	$Q_4$	$580.64$	$1516.7$	$903.53$
Magnet B	$Q_1$	$50.4$	$0.16$	$0$
	$Q_2$	$50.4$	$58.83$	$35.62$
Magnet D	$Q_5$	$516.12$	$207.12$	$353.39$
	$Q_6$	$516.12$	$321.6$	$577.56$
$J_1$		0.3420	0.0786	0.0784

, the simulation results are plotted as red lines in Figure 10. To verify the accuracy of the model, the total restoring forces of the system under various configurations are measured by using the restoring force surface method [33]. The results corresponding to the five chosen cases are plotted as black dots in Figure 10. By comparing the measured data and the simulation results using the original values of $Q_1$ to $Q_6$, it can be found that the model fails to predict the magnitudes of Cases III, IV, and V or the bi-stable and tri-stable cases.

To improve the accuracy of the model, a genetic algorithm-based identification approach proposed in [33] is applied. In this approach, $Q_1$ to $Q_6$ are treated as six independent parameters to be identified by minimizing an objective or fitness function, defined below:

$$J_1(Q_1,Q_2,Q_3,Q_4,Q_5,Q_6)=\sqrt{\frac{1}{5N_d}{\displaystyle \sum _{j=1}^5{\displaystyle \sum _{i=1}^{N_d}{(f_{jm}(x_i)-f_{ja}(x_i))}^2}}}$$

(19)

where $f_{jm}\left(x_i\right)$ is the measured restoring forces that are smoothened by a spline fitting, $f_a\left(x_i\right)$ is the analytical restoring forces based on Equation (14), and $N_d=101$ is the number of training data for each case, according to [33]. Once the six parameters have been identified, the neglectable parameter (with an almost zero value) can be set to zero, then an optimization for the five independent parameters can be conducted. All the identified values of the total charges and their corresponding fitness values are listed in Table 4. As shown in the table, the five-parameter optimization has the lowest fitness value. With the results, the recalculated restoring forces are plotted as blue lines in Figure 10.

In what follows, the optimal values with the five-parameter optimization are used. By integrating the total restoring forces with respect to $x$, the potential energies of the five cases can be found and plotted in Figure 11a. By varying the tuning parameters $d$ and $h$, the stability state region can be generated and plotted in Figure 11b. This figure reveals the tunability of the system. For both the lower limit and upper limit of the parameter $d$, the system is a mono-stable one, regardless of the value of $h$. To have a bi-stable system, the distance $d$ should be around the middle of the tuning range so that the repelling force between magnets A, C and magnet D is strong enough. And by maintaining a decrease in the parameter $d$, a tri-stable system can be achieved. Figure 12a shows the potential energies vs. $x$ and $d$ by fixing h at 0.02 m where C1, C2 and C3 represent the crossing points of the line C and the borderline between the strong mono-stable and tri-stable, the tri-stable and bi-stable, and the bi-stable and week mono-stable, respectively, while Figure 12b shows the potential energies vs. $x$ and $h$ by fixing d at 0.035 m, where D1 and D2 represent the crossing points of the line D and the borderline between the medium mono-stable and tri-stable and the tri-stable and bi-stable, respectively. It can be found that the region for the tri-stable is the narrowest one, which indicates that the tri-stable system has the highest sensitivity when changing the parameters.

5. Evaluation of the Performances of the MSHEH

In this section, the performances of the MSHEH are evaluated under harmonic frequency sweeping excitation. Both up-sweep and down-sweep excitations are conducted numerically and experimentally. For this purpose, four different configurations are considered: linear, mono-stable, bi-stable and tri-stable. The restoring force of the linear system is defined by setting $f_{DBx}$, $f_{DAx}$, $f_{DCx}$ to zero in Equation (14), while the last three configurations are Case (I), Case (III) and Case (IV), defined in the previous section. The load resistances for the EMEH and PEH are set to the boldfaced values in Table 7.

5.1. High-Level Acceleration

In the simulation, the acceleration amplitude is set to $1.6$ $\mathrm{m}/{\mathrm{s}}^2$. The frequency of the harmonic excitation is varied by

$$f=f_i+\frac{f_e-f_i}{T}t$$

(20)

where for the up-sweep, $f_i=2$ Hz and $f_e=8$ Hz, and for the down-sweep, $f_i=8$ Hz and $f_e=2$ Hz, and $T=360$ s. The initial conditions are set as $x\left(0\right)=\dot{x}\left(0\right)=0$ for the linear, mono-stable and tri-stable system, and $x\left(0\right)=-0.019$ m and $\dot{x}\left(0\right)=0$ for the bi-stable system. The numerical simulation is conducted by solving Equations (1), (2), (4) and (9) using Matlab function ode45. The instant power output is used to measure the energy harvesting performance of the system, which can be calculated from the instant voltage over the load resistor and the corresponding load resistance value. Considering the instant power-output results for different configurations may overlap each other, to better represent all the results, the amplitudes of the instantaneous power outputs are obtained by picking the upper envelopes of the instant power-output signals. The simulation results are shown in Figure 13, where ${\widehat{P}}_{em}$ and ${\widehat{P}}_p$ are the instantaneous power of the EMEH and PEH, respectively.

Figure 13a clearly shows that for the EMEH, the linear and tri-stable configurations outperform the mono-stable and bi-stable configurations in terms of the peak output powers. The EMEH with the mono-stable and tri-stable configurations show obvious hardening behaviors, which leads to a wider effective energy harvesting bandwidth. With the frequency rising, the bi-stable system first switches between the intra-well oscillation and chaotic inter-well oscillation at 3.7 Hz and then resumes the intra-well motion after 4.4 Hz. And, owing to its lower potential barriers, the tri-stable system starts with the periodic inter-well oscillation at 2 Hz and then switches to the intra-well oscillation at 4.3 Hz. It can be seen that the periodic inter-well motion of the tri-stable system generates more power compared to the chaotic inter-well motion of the bi-stable system, and the intra-well motion has the lowest energy harvesting efficiency among the three motion modes.

As shown in Figure 13b, the trends for the power outputs of the PEH of the four configurations are similar to those of the EMEH. In addition, it should be noted that the value of the instant power output of the EMEH under the mono-stable and tri-stable configurations are close to each other from 2 Hz to 3.4 Hz. However, the value of the power output of the PEH under the mono-stable configuration is lower than that under the tri-stable configuration. The main reason is that the voltage output of the EMEH is only related to the transduction factor and the velocity of the moving magnets according to Equation (9). Since both the mono-stable system and tri-stable system perform the large amplitude oscillation under the low-frequency excitation (lower than 3.4 Hz), the velocity of the moving magnets of the two systems are close when passing the high-power output regions ($x=\pm 0.017$ m), which explains the similar power output level. On the other hand, the power output of the PEH mainly depends on the displacement of the cantilever beam’s tip. The two side potential wells of the tri-stable system lead to a larger amplitude response at the inter-well oscillation mode than that of the mono-stable system. Thus, the PEH with the tri-stable configuration shows higher power output than the PEH with the mono-stable configuration.

To validate the above simulation, an experiment is conducted. As shown in Figure 14, the apparatus is fixed on a slip table that is driven by a shaker (2809, Brüel and Kjær) through a stinger. The shaker is driven by an amplifier (2718, Brüel and Kjær). Two laser reflex sensors (CP24MHT80, Wenglor) are used to measure the displacements of the beam’s tip and the base, respectively. A computer equipped with a dSPACE dS1104 data acquisition board is used to collect the signals from the laser sensors and the voltage signals of the EMEH’s load resistor and the PEH’s load resistor, and to output the excitation voltage signal to the power amplifier. To control the experiment, a program is developed by using the MATLAB Simulink which is interfaced with dSPACE Controldesk Desktop software.

The experimental results are shown in Figure 15. It can be seen that the trends in the results agree with the simulation ones for the higher-frequency region. However, for the lower-frequency region, the experimental results for the four systems are lower than their simulation counterparts. Such discrepancy can be attributed to the limit of the shaker, because 4 Hz exceeds the lower limit of the ideal working range of the shaker, which causes the actual acceleration of the excitation to be much lower than 1.6 $\mathrm{m}/{\mathrm{s}}^2$. Nevertheless, the experiment results indicate that the model used in the simulation is valid. In what follows, more simulation is carried out to further evaluate the performances of the MSHEH.

Figure 16 shows the simulation results of the power outputs for the EMEH and PEH under the frequency down-sweep (8 to 2 Hz) excitation. The overall trends in the power outputs of the linear, mono-stable and bi-stable configurations are similar to those in the frequency up-sweep excitation. The power output of the tri-stable system is not as high as in the up-sweep test, since it mainly performs the intra-well oscillation. The bi-stable and tri-stable systems start to jump when the frequency decreases to 4.3 Hz and 3.5 Hz, respectively, which is the situation when the two systems just overcome the threshold of their local potential well and switch to the chaotic inter-well oscillation mode. As shown in Figure 17, the experimental results generally agree with the simulation ones.

To better measure the bandwidth of the MSHEH, the accumulated harvested energy $E_{em}$ of the EMEH and the accumulated harvested energy $E_p$ of the PEH are defined as:

$$E_{em}(t)={\displaystyle {\int }_0^t{P}_{em}(\tau )d\tau }$$

(21)

$$E_p(t)={\displaystyle {\int }_0^t{P}_p(\tau )d\tau }$$

(22)

where $P_{em}$ and $P_p$ are the instantaneous power of the EMEH and PEH, respectively. Figure 18a,b show the total accumulated harvested energy $E\left(t\right)=E_{em}\left(t\right)+E_p\left(t\right)$ under the frequency up-sweep excitation and down-sweep excitation, respectively. Note that the relationship between the time $t$ and the frequency $f$ is defined by Equation (20). The effective frequency range of energy harvesting can be defined as the region where the increase rate of $E\left(t\right)$ is equal to or greater than 0.1 $\mathrm{J}/\mathrm{Hz}$. The total bandwidth can be obtained by taking the sum of the frequency range of the up-sweep and down-sweep tests. For example, the effective energy harvesting bandwidth for the tri-stable system is 2.36 Hz (ranging from 2 Hz to 4.36 Hz) and 0.1 Hz (ranging from 2 Hz to 2.1 Hz) for the frequency up-sweep and down-sweep excitations, respectively, and the total bandwidth will be 2.36 Hz. In addition, the total accumulated harvested energy for each configuration is represented by $E_{em}\left(T\right)$ under the up-sweep excitation. All the results are listed in

Table 5. The effective frequency bandwidths of the systems under the high-level frequency sweep excitation.

	Up-Sweep Frequency Range (Hz)	Down-Sweep Frequency Range (Hz)	Total Bandwidth (Hz)	$\mathit{E}\left(\mathit{T}\right)$ (J)
Linear	4.53–5.23	4.54–5.24	0.71	1.78
Mono-stable	2.00–3.73	2–2.88	1.37	1.62
Bi-stable	3.68–4.83	3.48–5.26	1.78	0.72
Tri-stable	1.66–4.36	2.00–2.10	2.36	3.86

, and the results clearly show all three nonlinear configurations show wider bandwidth compared to the linear one. The tri-stable system has the largest total bandwidth and total accumulation of harvested energy among the four configurations, as it is able to enter the periodic inter-well oscillation mode.

5.2. Low-Level Acceleration

To investigate the energy harvesting performances of the system under the excitation with low-level acceleration, a series of simulations are conducted. In the simulation, the acceleration amplitude is set to $0.3$ $\mathrm{m}/{\mathrm{s}}^2$ and the frequency varies from 2 to 8 Hz (up-sweep) and 8 to 2 Hz (down-sweep) in a duration of 360 s, and the same initial conditions as those in the previous section are used. In this time, the bi-stable and tri-stable will perform the low-amplitude intra-well oscillation at their side potential wells and middle potential well, respectively.

By following the same procedure discussed previously, the $E\left(t\right)$ for the system under frequency up-sweep and down-sweep excitation can be obtained, and is shown in Figure 19. As the values of $E\left(t\right)$ under low-level excitation are much lower than those under the high-level excitation, the threshold of the increase slope is chosen as $2\times {10}^{-3}$ J/Hz when the effective frequency range is identified. The effective bandwidth of the system and the $E\left(T\right)$ of the different configurations are both summarized in

Table 6. The effective frequency bandwidths of the systems under low-level frequency sweep excitation.

	Up-Sweep Frequency Range (Hz)	Down-Sweep Frequency Range (Hz)	Total Bandwidth (Hz)	$\mathit{E}\left(\mathit{T}\right)$ (J)
Linear	4.60–5.10	4.62–5.12	0.52	$3.82\times {10}^{-2}$
Mono-stable	2.00–2.37	2.00–2.08	0.37	$1.69\times {10}^{-2}$
Bi-stable	3.07–5.34	4.14–5.72	2.65	$9.35\times {10}^{-2}$
Tri-stable	3.07–3.44	2.00–3.66	1.66	$1.32\times {10}^{-2}$

. The results indicate that the bi-stable and tri-stable possess a wider efficient energy harvesting bandwidth, and the bi-stable system, in particular, can harvest more energy compared to other configurations. The higher energy harvesting efficiency is attributed to the high-power output regions of EMEH located around the equilibrium points of the bi-stable system (at $\left|x\right|=0.02 \mathrm{m}$). And this can guarantee the high efficiency of energy harvesting even when the bi-stable system performs the low-amplitude intra-well oscillation.

6. Pareto Front Optimization

To maximize the power output of the system, it is crucial to determine the optimum resistance value. Traditionally, this involves applying impedance matching to each component within a hybrid energy harvester [19]. However, the complex coupling effect between the PEH and the EMEH warrants further consideration. A traditional impedance matching may not be sufficient to ensure the optimum overall performance of the system. The explanation is shown as follows: in this apparatus, the deployment of a large number of coils results in a significantly high peak value for $K_t$, leading to substantial electromagnetic damping forces from the EMEH. The force will significantly impact the dynamics of the system, particularly in multi-stable configuration cases. Although higher currents in the EMEH can increase its power output, the resulting large damping force may hinder the system from performing the inter-well oscillations. In other words, increasing the power output of the EMEH may scarify the power output of the PEH. Therefore, a proper compromise between the power output of EMEH and PEH needs to be considered when one chooses the optimum $R_{le}$ and $R_{lp}$.

In this study, the MATLAB Global Optimization Toolbox based on the genetic algorithm is employed to solve such a multiple-objective optimization problem. The average power outputs of the EMEH and PEH are defined as follows:

$${\overline{P}}_{em}=I_{rms}^2{R}_{le},$$

(23)

and

$${\overline{P}}_p=\frac{v_{rms}^2}{R_{lp}},$$

(24)

respectively, where $I_{rms}$ and $v_{rms}$ are the root mean square value of the output current of EMEH and voltage of the PEH, respectively. The search range is from 0.1 $\mathsf{\Omega }$ to 300 $\mathsf{\Omega }$ for $R_{le}$ and from 0.1 MΩ to 5 MΩ for $R_{lp}$. Since the program is based on the minimization of the objective functions, the two objective functions are set to ${-\overline{P}}_{em}$ and ${-\overline{P}}_p.$ The population size and the maximum number of the generation are set to 500 and 50, respectively; and the same initial conditions as those in Section 5 are used, and the amplitude of the acceleration of the harmonic excitation is set to $2 \mathrm{m}/{\mathrm{s}}^2$.

After implementing the optimization program for the four systems under the excitation with six different frequencies (2.5 Hz, 3 Hz, 3.5 Hz, 4 Hz, 4.5 Hz and 5 Hz), the best results of the so-called Pareto front are shown as black dots in Figure 20. To find the best trade-off point, the distance between the origin of the plot and each best result is evaluated. The point with the shortest distance is considered to have the best trade-off between ${-\overline{P}}_{em}$ and ${-\overline{P}}_p$, shown as red dots in Figure 20. Then, the total power output ${\overline{P}}_t$ is the sum of ${-\overline{P}}_{em}$ and ${-\overline{P}}_p$, corresponding to this point. It is important to note that the results presented on the Pareto front offer decision support for configuring the system to meet diverse application requirements. In practical scenarios, the priority may lean towards either the EMEH or the PEH, dictating that the optimal point could be selected from either the left or right side of the best trade-off point identified in this study.

Table 7. The optimum load-resistance values and total power output for the four configurations under the harmonic excitation.

Configuration	Frequency (Hz)	${\mathit{R}}_{\mathit{le}}^{\ast }$$(\mathsf{\Omega }$)	${\mathit{R}}_{\mathit{lp}}^{\ast }$$(\mathbf{M}\mathsf{\Omega }$)	${\overline{\mathit{P}}}_{\mathit{em}}$ (mW)	${\overline{\mathit{P}}}_{\mathit{p}}$ (mW)	${\overline{\mathit{P}}}_{\mathit{t}}$ (mW)
Linear	2.5	82.612	1.49	0.03	0.004	0.03
	3	19.25	1.11	0.44	0.04	0.48
	3.5	13.76	0.92	0.10	0.02	0.12
	4	44.69	0.77	0.44	0.04	0.48
	4.5	49.30	0.74	6.45	0.16	6.61
	5	204.18	0.60	47.10	1.54	48.64
Mono-stable	2.5	43.00	1.23	10.34	0.23	10.56
	3	54.36	1.04	21.36	0.39	21.76
	3.5	24.42	0.93	0.22	0.02	0.24
	4	14.92	0.87	0.08	0.01	0.09
	4.5	15.26	0.77	0.03	0.01	0.04
	5	15.97	0.67	0.02	0.01	0.02
Bi-stable	2.5	69.19	1.03	0.54	0.01	0.54
	3	256.11	0.73	12.76	0.75	13.51
	3.5	282.41	1.37	19.07	1.14	20.21
	4	43.90	0.79	12.94	0.12	13.06
	4.5	123.43	1.23	5.75	0.16	5.90
	5	75.73	0.63	6.90	0.06	6.97
Tri-stable	2.5	94.56	1.41	17.76	0.55	18.31
	3	164.22	1.04	19.09	0.85	19.94
	3.5	54.67	1.58	4.81	0.12	4.93
	4	46.17	2.06	4.35	0.10	4.46
	4.5	34.74	0.61	7.26	0.05	7.31
	5	19.07	0.62	0.10	0.02	0.11

liststhe optimum results $R_{le}^{\ast }$ and $R_{lp}^{\ast }$, and the corresponding ${\overline{P}}_t$ of the four configurations when reaching the maximum at a different exciting frequency: 5 Hz (linear); 3 Hz (mono-stable); 3.5 Hz (bi-stable); and 3 Hz (tri-stable). As summarized in [19], the power output of HEH consisting of a PEH and EMEH generally ranges from 1 $\mathsf{\mu }W$ to 100 mW; therefore, the power output level for the proposed apparatus in this study is considered reasonable.

In the above optimization, the maximum harvested powers of the EMEH and PEH are chosen as the objective functions and the harmonic excitation with a constant frequency is considered. The result shows that the linear configuration outperforms the other three ones, confirming the well-known knowledge that the linear energy harvester is the best choice if the ambient vibration is harmonic, with a fixed frequency. In the previous frequency sweep-excitation simulation, the best compromised values in Table 7 were used in order to compare the four configurations, based on the benchmark of the linear configuration. The results have shown that the nonlinear configurations outperform the linear one in terms of the accumulated harvested energy and the frequency bandwidth. A natural question arises as to what the best load resistances are if the MSHEH is subjected to a frequency sweep excitation and the accumulated harvested energies are chosen to be the objective functions.

To answer this question, a further optimization is conducted. The MSHEH is subjected to the high-level frequency up-sweep excitation. The two objectives are set as $-E_{em}\left(T\right)$ and $-E_p\left(T\right)$, respectively. By following the same simulation procedure as outlined for the high-level frequency up-sweep tests in the previous section, the $E\left(t\right)$ for each configuration can be obtained. The setting of the optimization is the same as above, and the same initial conditions are used as those in Section 5. Considering the computational cost, the duration of the excitation signal is chosen as $T=100$ s. The obtained Pareto fronts for the four configurations are shown in Figure 21, where the best trade-off points are identified by red circles. The optimum resistance values $R_{le}^{\ast }$ and $R_{lp}^{\ast }$, and the corresponding accumulated harvested energy for the EMEH and PEH $E_{em}^{\ast }\left(T\right)$ and $E_P^{\ast }\left(T\right)$, and total accumulated harvested energy $E^{\ast }\left(T\right)$ are listed in

Table 8. The optimum load-resistance values and total power output for the four configurations under the frequency up-sweep excitation.

Configuration	${\mathit{R}}_{\mathit{le}}^{\ast }$$(\mathsf{\Omega }$)	${\mathit{R}}_{\mathit{lp}}^{\ast }$$(\mathbf{M}\mathsf{\Omega }$)	${\mathit{E}}_{\mathit{em}}^{\ast }\left(\mathit{T}\right)$ (J)	${\mathit{E}}_{\mathit{p}}^{\ast }\left(\mathit{T}\right)$ (J)	${\mathit{E}}^{\ast }\left(\mathit{T}\right)$ (J)
Linear	194.80	0.54	0.24	$1.0\times {10}^{-2}$	0.25
Mono-stable	256.86	0.58	0.95	$5.0\times {10}^{-2}$	1.00
Bi-stable	66.28	0.91	0.26	$0.3\times {10}^{-2}$	0.26
Tri-stable	258.21	1.09	1.26	$7.6\times {10}^{-2}$	1.34

. It can be seen that $E^{\ast }\left(T\right)$ for nonlinear configurations outperforms that for the linear configuration. Figure 22 compares the $E\left(t\right)$ (solid lines) of the MSHEH with the optimum load resistances from Table 7, referred to as Opt 1 and those (dashed lines) with the optimum load resistances from Table 8, referred to as Opt 2. Several observations can be made. Here, $E\left(T\right)$ represents the $E\left(t\right)$s’ value at 8 Hz. Firstly, the $E\left(T\right)$s from the linear configuration remain almost unchanged for both cases. Secondly, the $E\left(T\right)$ from the bi-stable configuration for Opt 2 sees an increase compared with that for Opt 1. Thirdly, the $E\left(T\right)$ from the mono-stable configuration for Opt 2 increases significantly. Fourthly, the tri-stable configuration for Opt 2 still exhibits the best performance when compared to the other three.

7. Conclusions

In this study, we present the development and evaluation of a multi-stable hybrid energy harvester (MSHEH). The system is equipped with both an electromagnetic energy harvester (EMEH) and a piezoelectric energy harvester (PEH), offering two tuning variables (h and d) for selecting the different stability states. A novel arrangement of coils in the EMEH has been implemented to enhance energy harvesting efficiency across various oscillation modes. A numerical approach is employed to determine the transduction factor for the EMEH. The obtained results are validated experimentally. The magnetic restoring force model is established based on the equivalent magnetic 2-point dipole model and is validated experimentally. The accuracy of the model is further improved by the genetic-algorithm identification approach. This refined model was used to map the stability state region. Four different configurations of the MSHEH, namely linear, mono-stable, bi-stable, and tri-stable, were chosen to evaluate the energy harvesting performances of the MSHEH through both simulation and experiment.

In the performance evaluation, the MSHEH’s four configurations are subjected to frequency up-sweep or down-sweep base excitation with high-level acceleration and low-level acceleration, respectively. The results revealed that under the high-level excitation, the mono-stable and multi-stable configurations exhibit a wider working bandwidth than the linear one. Particularly, owing to the shallower barrier of the potential wells, the tri-stable system is able to perform the large-amplitude periodic inter-well oscillation, which means that it has the widest frequency bandwidth (2.36 Hz) and highest total accumulated harvested energy (3.86 J) among the four configurations. When the system is under low-level excitation, both bi-stable and tri-stable harvesters perform the low-amplitude intra-well oscillation around the side potential wells and the middle potential well, respectively. The results show that the bi-stable system outperforms the others in terms of effective bandwidth (2.65 Hz) and total accumulated harvested energy ($9.35\times {10}^{-2}$ J). Due to the high-power output regions of the EMEH being located around the two side equilibriums of the bi-stable configuration, the EMEH’s power output remains sufficiently high, even though the system only performs low-amplitude intra-well oscillations.

In the end, a Pareto front optimization is employed to find the optimum values for $R_{le}$ and $R_{lp}$ by balancing the power output for the EMEH and PEH when the system is under harmonic excitation with various frequencies. The results demonstrate that the value of the optimum $R_{le}$ is higher when the amplitude of the oscillation is larger, and the values of the optimum $R_{lp}$ are inversely proportional to the frequency of the excitation. In addition, another Pareto optimization is conducted to further improve the accumulated harvested energy for both EMEH and PEH under the high-level frequency up-sweep excitation. The results demonstrate that the total accumulated harvested energies of the nonlinear configurations outperform the linear one.