Theoretical model and analysis of an electromagnetic vibration energy harvester with nonlinear damping and stiffness based on 3D MEMS coils

Based on the structure and fabrication characteristics of 3D coils published in [8], an innovative MEMS VEH is designed (figure 1(a)). This device includes a permanent magnet, four pairs of planar springs, a pair of 3D coils, and a silicon steel sheet. The overall dimension of the device is 20 mm × 20 mm × 1 mm, which is made by bonding two silicon wafers. The springs are processed by wire cutting and are glued on both sides of the bonded substrate. The permanent magnet is placed on the center plane. By installing a shim at the bottom of the permanent magnet, the initial position of magnet can be offset from the center of silicon steel sheet, which is shown as x₀ in figure 1(b).

Zoom In Zoom Out Reset image size

3D MEMS coils are first applied in EM-VEHs, and a silicon steel sheet is integrated in this structure to constrain the magnetic field. The magnetising direction of the permanent magnet is parallel to the coils axis direction and perpendicular to the vibration direction. Compared with the traditional MEMS VEHs using planer coils, our model has a much higher magnetic flux gradient.

Figure 1(b) shows the symbols and geometrical parameters of the magnetic path. The air gap δ and initial magnet offset (IMO) x₀ can significantly affect the dynamic response and output performance. Therefore, the influences of these two parameters are analyzed and optimised in this study.

The proposed EM-VEH is a typical 'spring-mass-damper' vibration system. Figures 1(c) and (d) show the vibration diagrams, where m is the mass of the vibration block. The stiffness terms include the mechanical stiffness K_m (the resilience of the springs) and the electromagnetic stiffness K_e (the attraction between the magnet and the silicon steel sheet). When the eddy current loss in the silicon steel core is ignored, the electromagnetic damping d_e is related to the interaction between the magnet and coils, which leads to the conversion of vibration energy to electricity. d_m is the mechanical damping which mainly includes the air damping effect. Considering that the air damping is significantly small compared with the electromagnetic damping, and that the proposed device can apply vacuum packaging technology to further reduce the energy dissipation, the mechanical damping d_m is set to be zero in this study.

The equation of motion for the system is established by summing the dynamic forces on the block and this is described as follows:

$$\begin{equation}\ddot{x}+ \frac{{{d_e} + {d_m}}}{m}\dot{x} + \frac{{{K_m}(x) + {K_e}(x)}}{m} = - \ddot{y}. \end{equation} \tag{ 1 }$$

The following sections describe the methods for modeling the stiffness and damping coefficient to predict the dynamic response and output performance accurately.

According to Faraday's law of electromagnetic induction, the electromotive force V can be derived as shown in equation (2):

$$\begin{equation}V{\text{ = }}\frac{d}{{dt}}\left(\sum {{\varphi _i}} \right) = \frac{{Nd\varphi }}{{dt}}.\end{equation} \tag{ 2 }$$

Here, ϕ_i is the magnetic flux through the ith loop, N is the total number of coil turns, and ϕ is the average magnetic flux along the axis of silicon steel. As the silicon steel can constraint the magnetic field and reduce the magnetic leakage, the magnetic flux along the coil axis can be assumed constant. Then the power consumed by the coils and the load circuit can be calculated using equation (3).

$$\begin{equation}P = \frac{{{V^2}}}{{{R_c} + {R_L}}} = \frac{{{{(N\frac{{d\varphi }}{{dx}})}^2}}}{{{R_c} + {R_L}}}{\left( {\frac{{dx}}{{dt}}} \right)^2}.\end{equation} \tag{ 3 }$$

Here, R_c is the resistance of coil windings, which can be calculated using the classical direct-current resistance formula, and R_L is the load resistance. Simultaneously, the power extracted by coils can also be expressed by the product of damping force and speed, as shown in equation (4), which is equal to the power consumed by the circuit.

$$\begin{equation}P = {d_e}{\left( {\frac{{dx}}{{dt}}} \right)^2}.\end{equation} \tag{ 4 }$$

Hence the electromagnetic damping coefficient can be derived based on the conversion of energy, as shown in equation (5). Table 1 lists the coil parameters and impedance data.

$$\begin{equation}{d_e} = \frac{{{{(N\frac{{d\varphi }}{{dx}})}^2}}}{{{R_c} + {R_L}}}.\end{equation} \tag{ 5 }$$

The magnetic flux through one cross section can be calculated by:

$$\begin{equation}\varphi {\text{ = }}\int {{B_{axis}}dS} = {B_{ave}}(WT).\end{equation} \tag{ 6 }$$

Here, B_ave is the average magnetic flux density normal to the section, and WT is the area of cross section. The magnetic field is generated by the magnet and the magnetised silicon steel. Hence, the magnetic flux density is a complicated nonlinear function of magnet size, silicon steel dimension, air gap, and magnet position. It is significantly difficult to derive a precise theoretical model to calculate the magnetic field distribution. Engel-Herbert et al presented an analytic expression for magnetic field of a single cubic permanent magnet in three dimensions [14]. The magnetic field strength along the magnetisation direction can be expressed by summing several arctangent terms of position coordinates. The magnetic field distribution of these two cases have the following common characteristics: (1) the magnetic induction intensity is symmetrically distributed, implying that the function is an even function: B(x) = B(−x) and dB/dx = 0 at x = 0 ; and (2) at infinity location the magnetic induction intensity decays to zero: B(±∞) = 0 . Inspired by this expression, in this study, the magnetic flux density is described by the arctangent and polynomial functions, as shown in equation (7).

$$\begin{equation}\frac{\varphi }{{WT}}{\text{ = }}{F_1}\arctan \left(\frac{1}{{\sqrt {{F_2} + {F_3}{{(\frac{x}{T})}^2} + {F_4}{{(\frac{x}{T})}^4}} }}\right).\end{equation} \tag{ 7 }$$

The ϕ/WT is the average magnetic flux density, which remains unchanged when the overall structure is scaled up or down, such that this function can be generalized to any similar geometrical structure. The coefficients F_i are related to the geometrical parameters, including the size of the air gap, magnet, and silicon steel sheet, and are also related to the coercive force of magnet. Equation (7) gives a functional form of the magnetic flux distribution, and the specific coefficient values need to be obtained by experiments or numerical simulations. In this study, magnetic flux data considering different positions and air gap sizes are simulated using MAXWELL 3D (figure 2(a)). The coercivity of magnet is 838 000 A m⁻¹, and the relative permeability is 1.045. The B-H curve data of silicon steel sheet is listed in table 2. Table 3 presents the coefficients values of different air gaps, which show good agreement with the simulation results and also proved the practicality and validity of equation (7).

Zoom In Zoom Out Reset image size

Combining equations (5) and (7), the damping coefficient can be derived as follows.

$$\begin{align} {d_e}(x) & = \frac{{{N^2}}}{{{R_c} + {R_L}}}\frac{{{{(WT{F_1})}^2}}}{4}\frac{1}{{({F_{\text{2}}} + {F_{\text{3}}}{{(\frac{x}{T})}^2} + {F_{\text{4}}}{{(\frac{x}{T})}^4})}} \nonumber \\ & \quad\times \mathop {}\limits^{} \mathop {}\limits^{} \mathop {}\limits^{} \frac{{\text{1}}}{{{{({F_{\text{2}}} + {F_{\text{3}}}{{(\frac{x}{T})}^2} + {F_{\text{4}}}{{(\frac{x}{T})}^4} + 1)}^2}}}{\left(\frac{{2{F_{\text{3}}}}}{{{T^2}}}x + \frac{{4{F_{\text{4}}}}}{{{T^4}}}{x^3}\right)^2}. \end{align} \tag{ 8 }$$

Owing to large deformation, the planar spring exhibits nonlinear stiffness characteristics along the thickness direction that can be approximated as a cubic polynomial. The spring parameters of this study are shown in figure 3, and the thickness is 60 μm. Every VEH has two planar springs arranged on the top and bottom sides of the silicon substrate.

Zoom In Zoom Out Reset image size

The nonlinear resilience of single planar copper spring is simulated using ANSYS, and is obtained by equation (9).

$$\begin{equation}\begin{gathered} {F_{ms}} = \frac{{\text{1}}}{{\text{2}}}({P_{m1}}x + {P_{m3}}{x^3}{\text{)}} \hfill \\ \mathop {}\limits^{} \mathop {}\limits^{} {\text{ = }}0.09748x + 3.44 \times {10^{ - 7}}{x^3} \hfill \\ \end{gathered}. \end{equation} \tag{ 9 }$$

The magnet is attracted by the magnetised silicon steel core as shown in figure 1(c). In [11], the interaction force between two permanent magnets are derived by equating magnets to magnetic dipoles, as shown in equation (10), which is a rational function of position x.

$$\begin{equation}{F_e} = \frac{{3\mu {m_1}{m_2}}}{{2\pi }}\frac{x}{{{{({x^2} + {d^2})}^{5/2}}}}.\end{equation} \tag{ 10 }$$

However, in this study, the size of the silicon steel core is too large to be assumed as magnetic dipoles. Additionally, the magnetization intensity in the silicon steel core varies with the magnet position. Therefore, the function form of equation (10) needs to be adjusted to fit the structure of this paper. Based on this, an empirical rational function, as in equation (11), is proposed to describe the magnetic force between the magnet and silicon steel core.

$$\begin{equation}\frac{F}{{WT}} = \frac{{{K_{11}}(\frac{x}{T})}}{{{K_{20}} + {K_{22}}{{(\frac{x}{T})}^2} + {K_{24}}{{(\frac{x}{T})}^4}}}.\end{equation} \tag{ 11 }$$

Here, F/WT is the attraction per unit facing area of magnet and silicon steel, which remains constant when the overall structure is scaled up or down. The coefficients K_ij can be calculated for different geometrical designs.

Similar to equation (7), the discrete values of attraction with different positions and air gaps are simulated (figure 2(c)). The fitting coefficients are then derived from these values (table 4). The fitting results agree well with the simulation data.

The total stiffness is the sum of the electromagnetic and mechanical stiffness. The IMO x₀ can be controlled by installing an adjust shim on the bottom of the permanent magnet. It is worth noting that, due to the attraction of silicon steel core, an original deformation of planar spring is unavoidable when IMO is not zero, and this is shown in figure 4. The original deformation H₀ is the real solution of equation (12).

$$\begin{equation}{P_{m1}}{H_0} + {P_{m3}}H_0^3 = {F_e}({x_0}).\end{equation} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

The relative position of the magnet, silicon steel core and planar spring at position coordinate x is shown in figure 5. For ease of reading, only one planar spring is plotted in figure 5. Then the restoring force produced by silicon steel core and planar spring at any position can be calculated using equation (13). The potential curve is derived by the integral of restoring force, and the equilibrium position is set to be the zero potential energy point.

$$\begin{align} F(x) & = {F_e}(x) + {F_m}(x - {H_0} - {x_0}) \nonumber \\ & \mathop {}\limits^{} \mathop {}\limits^{} \mathop {}\limits^{} = \frac{{(WT){K_{11}}(\frac{x}{T})}}{{{K_{20}} + {K_{22}}{{(\frac{x}{T})}^2} + {K_{24}}{{(\frac{x}{T})}^4}}} \nonumber\\ & \quad+ \mathop {}\limits^{} \mathop {}\limits^{} \mathop {}\limits^{} \mathop {}\limits^{} {P_{m1}}(x - {H_0} - {x_0}) + {P_{m3}}{(x - {H_0} - {x_0})^3} \end{align} \tag{ 13 }$$

Zoom In Zoom Out Reset image size

$$\begin{equation}P(x) = \int\limits_{{x_0}}^x {F(x)dx}. \end{equation} \tag{ 14 }$$

The air gap can influence the magnetic force, whereas the IMO can affect the mechanical stiffness. Therefore, the influences of the air gap and IMO on the stiffness characteristics are analysed in this section. Figures 6(a) and (b) show the stiffness and potential energy curves for different air gaps when the IMO is 0 μm. With the decrease in the air gap, the stiffness of the system is improved, and because of this a steep potential energy curve is generated. Figures 6(c) and (d) delineate the stiffness and potential energy for different offsets when the air gap is 0.7 mm. The curves for other air gaps are similar. For a nonzero IMO, the curves demonstrate asymmetrical characteristics. The stiffness is considerably flatter on the right side of the initial position than on the left side, which indicates a lower potential energy on the right side. An increase in the IMO leads to a rapid increase and decrease in the stiffness on the left and right sides, respectively.

Zoom In Zoom Out Reset image size

Based on the proposed model, the vibration dynamic of this system can be expressed using equation (15). The output power and other related performance parameters can be derived using equation (16). The average output power during steady-state operation is calculated by equation (17).

$$\begin{align} \ddot{x} & + \frac{1}{m}\Bigg[\frac{{{N^2}}}{{{R_c} + {R_L}}}\Bigg(WT\frac{d}{{dx}}\Bigg({F_1} \arctan \nonumber\\ & \quad \frac{1}{{\sqrt {{F_2} + {F_3}{{(\frac{x}{T})}^2} + {F_4}{{(\frac{x}{T})}^4}} }}\Bigg)\Bigg)^2\Bigg]\dot{x} + \frac{1}{m}({P_{m1}}x + {P_{m3}}{x^3}) \nonumber\\ & + \frac{{WT}}{m}\frac{{{K_{11}}(\frac{x}{T})}}{{{K_{20}} + {K_{22}}{{(\frac{x}{T})}^2} + {K_{24}}{{(\frac{x}{T})}^4}}} = - \ddot{y} \end{align} \tag{ 15 }$$

$$\begin{equation}{P_e} = {I^2}{R_L} = {\left(\frac{{N\frac{{d\varphi }}{{dx}}}}{{{R_c} + {R_L}}}\right)^2}{\left(\frac{{dx}}{{dt}}\right)^2}{R_L}\end{equation} \tag{ 16 }$$

$$\begin{equation}{P_{ea}} = \frac{1}{T}\int {{P_e}dt}. \end{equation} \tag{ 17 }$$

It is difficult to obtain the analytical solution of equation (15). Hence, in this section, the vibration equation is solved numerically. The excitation ÿ is given as a fixed frequency sine function, a₀sin(ωt). The equation is discretised using the Runge–Kutta method. Then, the steady vibration response for a given initial state at a certain excitation condition can be numerically solved. In this study, the initial velocity is set to be zero, and the initial position is the lowest point of the potential energy curve, which is also the IMO. By adjusting the frequency and extracting the vibration amplitude and output power during steady-state at each frequency, the amplitude–frequency curves of system vibration and output performance can be obtained. The frequency corresponding to the maximum amplitude is the resonance frequency. The amplitude– frequency characteristics are related to the air gaps, IMO, excitation intensity and load resistance. Hence, for each design parameter or operating condition, it is necessary to calculate the amplitude– frequency curve separately, and then compare its performance.

Equation (15) is a nonlinear second-order ordinary differential equation, which is very complicated. Hence it is difficult to derive the complete analytical solution. However, when the excitation acceleration is small enough, the vibration block only vibrates in a limited range, and then the system can be regarded as equivalent to a linear system which can be solved analytically. Therefore, in this study, the numerical method is verified by comparing it with the analytical solution under small excitation.

Equation (18) is the quasi-linear equation of equation (15). The constant electromagnetic damping coefficient d_e can be calculated using the gradient of magnetic flux at the bottom of potential well. The equivalent stiffness k is the first-order component of stiffness curve at the potential well bottom. The damping coefficient and stiffness can be calculated using equations (19) and (20).

$$\begin{equation}m \ddot{x} + d \dot{x} + kx = m{a_0}\sin \omega t\end{equation} \tag{ 18 }$$

$$\begin{equation}{d_e} = \frac{{{{\left( {N{{\left. {\frac{{d\varphi }}{{dx}}} \right|}_{x = {x_0}}}} \right)}^{\text{2}}}}}{{{R_c} + {R_L}}}\end{equation} \tag{ 19 }$$

$$\begin{equation}k = {\left. {\frac{{d({F_m} + {F_e})}}{{dx}}} \right|_{x = {x_0}}}.\end{equation} \tag{ 20 }$$

The displacement amplitude during the steady state can be derived from equation (18), as:

$$\begin{equation}X = \frac{{m{a_0}}}{{({d_e} + {d_m}){\omega _n}}}.\end{equation} \tag{ 21 }$$

Here, ω_n is the resonant frequency, which is $\sqrt {\frac{k}{m}}$. The average electrical energy extracted by coils is,

$$\begin{equation}{P_{e\_ave}} = \frac{{{m^2}a_0^2}}{2}\frac{{{d_e}}}{{{{\left( {{d_e} + {d_m}} \right)}^2}}}.\end{equation} \tag{ 22 }$$

Then, the power consumed by the load can be expressed as:

$$\begin{equation}{P_{L\_ave}} = \frac{{{m^2}a_0^2}}{2}\frac{{{d_e}}}{{{{\left( {{d_e} + {d_m}} \right)}^2}}}\frac{{{R_L}}}{{{R_c} + {R_L}}}.\end{equation} \tag{ 23 }$$

In this study, the output performances under 0.1 mg (1 g = 9.8 m s⁻²) excitation acceleration have been calculated numerically and compared with the linear analytical solution. Table 5 lists the structure parameters and the corresponding damping and stiffness coefficients.

The vibration amplitude and output power at resonant frequency calculated by equations (21 )–( 23) and numerical method are listed in table 6. It can be observed that when the excitation is small enough, the numerical results are in good agreement with the analytical solution of the linear system. Hence, the reliability of numerical method is verified.

The influences of load resistance to output power are analysed in this section. In this study, η is defined as the ratio of load resistance R_L and coil resistance R_c: η = R_L/R_c . When η is close to zero, most of the power extracted by coils is consumed by the coil resistance; therefore, the output power on the load is significantly low. Moreover, with the increase in load resistance, the system damping decreases, resulting in the reduction of the total power extracted by coils. Therefore, there is an optimal load resistance matching the coil resistance, which could maximise the output power.

Figure 7 shows the output power varying with resistance ratio for different parameters (δ–IMOs are 0.5 mm–25 μm, 0.5 mm–650 μm, 0.7 mm–150 μm, 1.0 mm–750 μm, separately). It can be observed that for different air gaps, the optimal resistance ratio is also different, which is 1.4, 1.6, 1.0 and 1.0 for 0.5 mm–25 μm, 0.5 mm–650 μm, 0.7 mm–150 μm and 1.0 mm–750 μm, respectively. For the linear system, the optimal load resistance can be derived based on equation (23), as:

$$\begin{equation}{R_L} = \frac{{{{\left( {N\frac{{d\varphi }}{{dx}}} \right)}^2}}}{{{d_m}}} + {R_c}.\end{equation} \tag{ 24 }$$

Zoom In Zoom Out Reset image size

The matching load resistance increases with increasing the gradient of the magnetic flux. However, due to the nonlinear characteristics, the specific value and changing rules of optimal load are significantly different from that given by equation (24). This is because in a nonlinear system, the relationship between the vibration amplitude and damping coefficient is not a simple inverse proportional function described by equation (21). However, it is affected by high-order nonlinear components of damping and stiffness curve.

The numerical results also indicate that the sensitivity between optimal load and structural parameters is reduced in this nonlinear system. For example, in the air gap of 0.5 mm, when the resistance ratio is 1, the output power is approximately 5% different from the maximum power. Hence in this study, the resistance ratio of 1.0 is used to perform a unified analysis of the output performance for different air gaps, which could make the load resistance match with the coil resistance in most cases.

Figure 8 shows that the output power under 1 g excitation acceleration (1 g = 9.8 m s⁻²) varies with frequency for different initial magnet offsets, when the air gap ranges from 0.5–1 mm. For the same air gap, the IMO affects the electromagnetic stiffness and flux gradient, which consequently influence the output performance of the VEH. Figure 9 demonstrates the variation in resonant frequency and corresponding output power density with the IMO. Owing to the nonlinear stiffness characteristics of planar springs, the resonant frequency increases slightly with increasing offset near the centre. Away from the centre, the electromagnetic stiffness is rapidly reduced, which gradual decreases the resonant frequency. As for the output power density, a peak is observed at positions of 25, 150, 250, and 750 μm for air gaps of 0.5, 0.7, 0.8, and 1.0 mm, respectively. Reducing the air gap can improve the output power density and simultaneously move the peak position closer to the centre. The maximum output power density is 608.8 μW cm⁻³ with an air gap of 0.5 mm and initial magnet offset of 25 μm, which is considerably higher than that of the previously published data on VEHs [13, 15–18] (table 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Contraposed to a 0.7 mm air gap, the influence of excitation acceleration is investigated, including 0.5 g, 1 g, 2 g, 3 g, 4 g and 5 g (1 g = 9.8 m s^-2). The power density and normalized power density (NPD) for 0–250 μm IMO under different excitation accelerations are simulated (figure 10(a)). When the excitation is less than 2 g, the maximum output power is observed at an IMO of 150 μm. For 3 g and 4 g, the optimal IMO increases to 200 μm. However, at 5 g excitation, the output power constantly increases with an increase in the IMO. To normalise the acceleration level, the NPD is calculated (figure 10(b)), which indicates that although the output power increases with excitation acceleration, the NPD reduces when the excitation is enhanced. For a 0.7 mm air gap, the NPD has a peak value of 521.14 μW cm⁻³ g⁻² with a 150 μm IMO and 0.5 g acceleration. The resonant frequency also slightly varies with acceleration, differing from the linear vibration mode. However, when the excitation exceeds 2 g, this variation trend changes at the 150 μm offset, also indicating that the dynamic vibration characteristic is closely related to the excitation level.

Zoom In Zoom Out Reset image size

The influence of air gap, IMO, and excitation can be intuitively determined from the vibration phase portrait. Figure 11 shows the phase locus of different initial magnet offsets under 1 g and 5 g excitations for a 0.7 mm air gap. The cloud chart describes the output power map for each vibration state (x, ). The cloud chart, calculated using equation (16), is only related to the damping characteristic, d_e(x). The excitation and IMO simply change the vibration locus in this chart.

Zoom In Zoom Out Reset image size

For 1 g acceleration (figure 11(a)), the maximum instantaneous power (∼ 200 μW) appears at a 150 μm IMO. For IMOs greater than 150 μm, although the magnetic gradient is larger, the phase locus is restricted in a low power region owing to the higher system stiffness, producing lesser electricity. For 5 g acceleration, the vibration phase locus continually increases with an increase in IMO (figure 11(b)), with the maximum instantaneous power occurring at 250 μm.