# Multi-modal vibration energy harvesting approach based on nonlinear oscillator arrays under magnetic levitation

I. Abed, N. Kacem, Noureddine Bouhaddi, M.L. Bouazizi

► To cite this version:

I. Abed, N. Kacem, Noureddine Bouhaddi, M.L. Bouazizi. Multi-modal vibration energy harvesting approach based on nonlinear oscillator arrays under magnetic levitation. Smart Materials and Structures, 2016, 25, pp.25018 - 25018. 10.1088/0964-1726/25/2/025018. hal-01512531

HAL Id: hal-01512531

https://hal.science/hal-01512531

Submitted on 23 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Multi-modal vibration energy harvesting approach
based on nonlinear oscillator arrays under magnetic
levitation

I. Abed<sup>a,b</sup>, N. Kacem<sup>a</sup>, N. Bouhaddi<sup>a</sup> and M.L. Bouazizi<sup>c,d</sup>

**Abstract.** We propose a multi-modal vibration energy harvesting approach based on arrays of coupled levitated magnets. The equations of motion which include the magnetic nonlinearity and the electromagnetic damping are solved using the harmonic balance method coupled with the asymptotic numerical method. A multi-objective optimization procedure is introduced and performed using a Non-dominated Sorting Genetic Algorithm (NSGA) for the cases of small magnet arrays in order to select the optimal solutions in term of performances by bringing the eigenmodes close to each other in terms of frequencies and amplitudes. Thanks to the nonlinear coupling and the modal interactions even for only three coupled magnets, the proposed method enable harvesting the vibration energy in the operating frequency range of 4.6 – 14.5 Hz, with a bandwidth of 190 % and a normalized power of 20.2 mW cm<sup>-3</sup>g<sup>-2</sup>.

PACS numbers: 41.20.Gz, 05.45.-a, 46.40.Ff, 46.15.Ff, 45.10.Db

<footnote><sup>a</sup>FEMTO-ST Institute, UMR 6174, Applied Mechanics Department, University of
Franche-Comté, UBFC, 24 chemin de l'Epitaphe, 25000 Besançon, France</footnote>

<footnote><sup>b</sup>National Engineering School of Tunis (ENIT), Tunis El Manar University (UTM),
Tunisia</footnote>

<footnote><sup>c</sup>Preparatory Engineering Institute of Nabeul (IPEIN), 8000 M'Rezgua, Nabeul, Tunisia</footnote>

<footnote><sup>d</sup>Mechanical Department, College of Engineering, Prince Sattam Bin Abdulaziz University, KSA</footnote>
# 1. Introduction

The natural energy sources have been frequently taken to produce inexhaustible electric energy in their local environment. The purpose is related to the reduction of power requirement and to replace a battery that has a limited lifetime, requires maintenance and cannot be used in hostile environments. Then, the mechanical structures were exposed to ambient and internal vibrations to produce energy by using adequate transducers [1] which convert the mechanical force to electric power. The accumulated energy in these cases can be stored and used. Within the emerging trend of smart systems, Vibration Energy Harvesters (VEHs) can be used in many fields such as environmental monitors, wireless sensors and medical implants [2, 3, 4].

In order to make VEHs usable, there are several types of transduction, where the most common transduction modes are piezoelectric [5, 6, 7, 8, 9], and electromagnetic [10, 11, 11, 12]. The conventional linear VEHs are usually designed to be resonantly tuned to the ambient dominant frequency. They have a narrow operating bandwidth that limits their application in real-world environments where the ambient vibrations have their energy distributed over a wide spectrum of frequencies, with significant predominance of low frequency. In order to overcome this limitation, several approaches have been proposed such as resonant frequency tuning techniques, multi-modal energy harvesting and the exploitation of nonlinearities.

Resonant frequency tuning technique can be implemented in a manual way by applying a preload [13, 14, 15], adjusting the pre-deflection [16], regulating the distance between the magnets with using a spring-screw mechanism [17] and varying the vertical relative distance to modify the attractive magnetic force [18] and adjusting the gravity center of the tip mass [18] or in a self-tuning way [19], applying voltage to the transducer [20], switching the shunt electrical load [21] and implementing an automatic controller
[22, 18].

The multi-modal energy harvesting technique can be implemented by exploiting
multiple bending modes of a continuous beam or by exploiting a cantilever array
configuration. For instance, Roundy et al. [23] proposed multiple proof masses attached
on a clamped-clamped beam and the implementation of this idea was realized by Yang
et al. [24] with an electromagnetic generator based on three fixed magnets along a beam
vibrating at the top of three coils and operating in three modes of vibration. Sari et
al. [25] studied a structure composed of different assembly generators resonating at
different resonance frequencies. This structure contains 35 beams in different lengths
and each beam supports a coil that moves to the top of a stationary magnet. Yang
et al. [26] analyzed a VEH based on coupled flexural vibration of two elastically and
electrically connected piezoelectric beams, while Eurturk et al. [27] modeled an L-shaped
cantilever piezoelectric energy harvester. Tadesse et al. [28] proposed a multimodal
hybrid harvester consisting of piezoelectric crystals bonded to a cantilever beam and at
its tip, a permanent magnet is attached to oscillate within a stationary coil fixed to the
top of the package. Arafa et al. [29] presented a 2-DOFs cantilever piezoelectric VEH
with a dynamic magnifier consisting of a spring-mass system which is placed between
the fixed end of the piezoelectric beam and the vibrating base structure, while Kim
et al. [30] suggested a VEH composed of two piezoelectric cantilevers coupled with a
common proof mass and utilizing both translational and rotational degrees of freedom.

Despite the fact that multimodal techniques enable wide bandwidth energy
harvesting, they require more sophisticated interface circuits than that for a single-
mode harvester due to the phase difference between the output signals in array
configurations [31]. Consequently, several researches have been oriented towards the
study of nonlinear systems. For instance, Daqaq et al. [32] reported softening
frequency response characteristics in a parametrically forced piezoelectric device with
structural nonlinearities. Masana and Daqaq [33] proposed a nonlinear analysis of
electromechanical generator composed by piezoelectric clamped clamped beam. Mann
and Sims [34] showed analytically and experimentally how magnetic levitation could
be used to extend device bandwidth through a hardening response. Nevertheless, this
extension is limited by dry friction dissipation phenomenon. Mahmoudi et al. [35]
proposed an alternative to overcome this issue by guiding the moving magnet vertically
in an elastic way by means of sandwich beams, combined with an hybrid piezoelectric
and electromagnetic transductions.

In this paper, unlike classical VEHs either linear [23, 25, 27, 30] and multimodal or
nonlinear and mono-frequency [32, 34, 35], to our knowledge we propose for the first time
a vibration energy harvesting approach based on an array of coupled levitated magnets
combining the benefits of nonlinearities and modal interactions. These benefits consist
essentially in enlarging the bandwidth and increasing the harvested power. The set of
coupled nonlinear equations of motion is solved using the harmonic balance method
(HBM) coupled with the asymptotic numerical method (ANM) [36]. Then, in order to
improve the performances of the proposed method in terms of bandwidth and harvested
power, a multi-objective optimization procedure is introduced and performed using
NSGA-II algorithm [37] for the cases of two and three coupled magnets. Comparisons
have been made between uncoupled VEHs based on a single levitated magnet, two and
three coupled levitated magnets and state of the art VEHs in order to emphasize the
high performances of the proposed multi-modal vibration energy harvesting method.
## 2. System modeling

### 2.1. *Design*

Inspired by vibration energy harvesting based on a single levitated magnet [34], we propose an extension of magnetic levitation for a multi-degree of freedom (MDOF) VEH. Using this approach, one can take advantage of modal interactions and magnetic nonlinearities in order to obtain high performances in terms of bandwidth and harvested power. As shown in Figure 1, the considered device is composed of $n + 2$ magnets $M_i$ where $i \in [0, n + 1]$. $M_0$ and $M_{n+1}$ are fixed with respect to a Teflon tube inside which $M_1, M_2...M_{n-1}$ and $M_n$ are subjected to magnetic levitation forces. All magnets are placed vertically in such a way that all opposed surfaces have the same pole and wire-wound copper coils are wrapped horizontally around the separation distance between each two adjacent magnets.

<figure><img src="image_2.png" /><figcaption>Figure 1: A schematic diagram of the MDOF VEH based on magnetic levitation.</figcaption></figure>
## 2.2. Magnetic forces

In order to describe the electromagnetic induction model represented in Figure 1, $n + 1$ reference frames have been applied. The first reference frame is fixed in space and is used to describe the motion amplitude $Y_0$ and excitation frequency $\Omega$ of the outer housing. The other reference frames (designed as $x_1, x_2...x_{n-1}$ and $x_n$) describe the motion of the moving magnets which are subjected to gravitational forces $\vec{P}_j = M_j\vec{g}$ and magnetic forces (Figure 2) expressed as follows [10]:

$$ \vec{F}_{j,j-1} = -\frac{\mu_0 Q_{M_{j-1}} Q_{M_j}}{4\pi (d - x_j + x_{j-1})^2} \vec{x} \qquad (1) $$

$$ \vec{F}_{j,j+1} = \frac{\mu_0 Q_{M_{j+1}} Q_{M_j}}{4\pi (d + x_j - x_{j+1})^2} \vec{x} \qquad (2) $$

<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic of a moving magnet subject to magnetic and gravitational forces.</figcaption></figure>

where $Q_{M_j}$ is the magnetic intensity of the $j$th magnet, $\mu_0$ is the magnetic permeability, $d$ is the gap separating each two adjacent magnets, $x_0 = x_{n+1} = Y$ is the excitation amplitude and $g = 9.81m/s^2$ is the gravitational constant. A change of variable $v_j = x_j - Y$ yields the resulting restoring force applied on each moving magnet
$M_j$:

$$ 
\begin{aligned}
\vec{F}_j^m &= \vec{F}_{j,j-1} + \vec{F}_{j,j+1} + \vec{P}_j \\
&= \left( \frac{\mu_0 Q_{M_j}}{4\pi} \left( \frac{Q_{M_{j-1}}}{(d+v_{j-1}-v_j)^2} - \frac{Q_{M_{j+1}}}{(d+v_j-v_{j+1})^2} \right) - M_j g \right) \vec{x}
\end{aligned} 
\quad (3) $$

where the magnetic intensity of each magnet $M_j$ is written as $Q_{M_j} = S_j H_{cj}$, with $S_j$ and $H_{cj}$ are respectively the cross section and coercive force of the magnet $M_j$.

At static equilibrium, the resulting restoring force verifies $\vec{F}_{(j)}^m = \vec{0}$. Assuming that the equilibrium position of each moving magnet is $v_j = 0$, we obtain the following equation:

$$ Q_{M_j} = \frac{4\pi M_j g d^2}{\mu_0 (Q_{M_{j-1}} - Q_{M_{j+1}})} \quad (4) $$

Substituting Equation (4) into Equation (3), yields:

$$ \vec{F}_j^m = \left( \frac{M_j g d^2}{(Q_{M_{j-1}} - Q_{M_{j+1}})} \left( \frac{Q_{M_{j-1}}}{(d+v_{j-1}-v_j)^2} - \frac{Q_{M_{j+1}}}{(d+v_j-v_{j+1})^2} \right) - M_j g \right) \vec{x} \quad (5) $$

### 2.3. Electromagnetic damping forces

When the device is subjected to an external mechanical vibration, each moving magnet $M_j$ oscillates around its equilibrium position and a current is induced in each coil as shown in Figure 3, resulting in the following electrical damping forces:

$$ F_j^e = \alpha_{j-1} i_{j-1} - \alpha_j i_j \quad (6) $$

where the electromechanical coupling coefficient $\alpha_j = N_j B l$, with $B$ the average magnetic field strength, $N_j$ is the number of coil turns and $l$ is the coil length.
<figure><img src="image_4.png" /><figcaption>Figure 3: Equivalent electro-mechanical system with n moving magnets based on magnetic levitation system.</figcaption></figure>

In order to overcome the issue of phase difference between the output currents in array configurations, an independent interface circuit is proposed for each current (Figure 3). By applying Kirchhoff's theory to the electrical circuits, we obtain the following equation:

$$R_j i_j = \alpha_j (\dot{v}_{j+1} - \dot{v}_j) - r_{int} i_j \quad (7)$$

Equations (6) and (7) give the following expression:

$$F_j^e = c e_{j-1} (\dot{v}_j - \dot{v}_{j-1}) - c e_j (\dot{v}_{j+1} - \dot{v}_j) \quad (8)$$

The electrical damping can be expressed as a function of the internal resistance of the coil, the resistance of the harvesting circuit and the electromechanical coupling coefficient $\alpha_j$ [9].

$$c_{ej} = \frac{\alpha_j^2}{R_j + r_{int}} \quad (9)$$

## 2.4. Equations of motion

The application of Newton's first law to each moving magnet $M_j$ leads to the following equation of motion:

$$M_j \ddot{v}_j + c \dot{v}_j + F_j^e + F_j^m = -M_j \ddot{Y}; \quad j = 1, 2, \dots, n \quad (10)$$
Expanding the non-linear magnetic forces written in Equation (3) in Taylor series up to the third order, a system of coupled nonlinear equations is obtained as follows:

$$ 
\begin{aligned}
M_j \ddot{v}_j + (c + c_{ej-1} + c_{ej}) \dot{v}_j - c_{ej} \dot{v}_{j+1} - c_{ej-1} \dot{v}_{j-1} & \\
+ k_{j,j+1} (v_j - v_{j+1}) + k_{j,j-1} (v_j - v_{j-1}) & \\
+ \alpha_{j,j+1} (v_j - v_{j+1})^2 - \alpha_{j,j-1} (v_j - v_{j-1})^2 & \\
+ \lambda_{j,j+1} (v_j - v_{j+1})^3 + \lambda_{j,j-1} (v_j - v_{j-1})^3 & \\
= -M_j \ddot{Y}; \quad j = 1, 2, \dots, n &
\end{aligned} 
\quad (11) $$

This expansion is valid for magnet displacements below 50% of the gap [35]. The change of variable $Q_{M_j} = Q'_{M_j} M_j$ is substituted into Equation (5), and assuming that the magnetic intensities are equal ($Q'_{M_{j-1}} = Q'_{M_{j+1}}$), the linear and nonlinear stiffnesses can be written as:

$$ 
\begin{aligned}
k_{j,j+1} &= \frac{2g}{d} \frac{M_j M_{j+1}}{(M_{j-1}-M_{j+1})} & k_{j,j-1} &= \frac{2g}{d} \frac{M_j M_{j-1}}{(M_{j-1}-M_{j-1})} \\
\alpha_{j,j+1} &= \frac{3}{2d} k_{j,j+1} & \alpha_{j,j-1} &= \frac{3}{2d} k_{j,j-1} \\
\lambda_{j,j+1} &= \frac{2}{d^2} k_{j,j+1} & \lambda_{j,j-1} &= \frac{2}{d^2} k_{j,j-1}
\end{aligned} 
\quad (12) $$

Equation (11) can be written in its matrix form. Then, in order to ensure the symmetry of the rigidity matrix, we assume that $k_{j,j+1} = k_{j+1,j}$ which results in the following relation between the magnets in term of mass

$$ M_{j-1} - M_{j+1} = M_j - M_{j+2}; \quad j = 1, 2, \dots, n-1 \quad (13) $$

Moreover, the positivity of the linear rigidity terms ($k_{j,j+1} \ge 0$) is conditioned by the following inequality:

$$ M_{j-1} - M_{j+1} \ge 0; \quad j = 1, 2, \dots, n \quad (14) $$
## 2.5. Power transferred to the electrical circuit

The magnetic transduction is ensured by $n - 1$ coils. The oscillations of the movable magnets cause magnetic field variations in the separation zones, which provides an induced current (Lenz's Law). The induced current can be expressed as a vibration velocity function $\dot{v}_j(t) = V_j \Omega \sin(\Omega t)$. The instantaneous electrical power $P(t)$ can be written as follows:

$$P(t) = \frac{\Omega^2}{2} \left( \sum_{j=1}^{n+1} (c_{ej}(V_{j-1} - V_j)^2) \right) (1 - \cos(2\Omega t)) \quad (15)$$

When the system reaches the steady state regime, the average power delivered to the electrical load takes the following form:

$$P_m = \frac{\Omega}{2\pi} \int_0^{\frac{2\pi}{\Omega}} P(t)dt = \frac{\Omega^2}{2} \left( \sum_{j=1}^{n+1} (c_{ej}(V_{j-1} - V_j)^2) \right) \quad (16)$$

Assuming that $r_{int} \ll R_j$, the mean power $P_m$ is a good approximation of the harvested power. The proposed configuration has the advantage of avoiding the phase difference between the induced currents. However, it requires sophisticated circuits in order to deliver the energy separately.

## 2.6. Solving procedure

The solving procedure is purely computational and based on the harmonic balance method (HBM) coupled with the asymptotic numerical method (ANM) [36]. This technique gives the periodic solutions of a dynamical system when a control parameter is varied. But its limitation is that the equation to solve should be quadratic and it is not always easy to recast any system to a polynomial quadratic form. Once the system is transformed into a quadratic form, the HBM is used to obtain the associated nonlinear algebraic system which is solved using the ANM. The application of this solving technique to Equations (11) and (16) is detailed in the Appendix.
### 3. Results and discussion

#### 3.1. Specifications of a nonlinear SDOF VEH

We start by calculating the response of a single nonlinear levitated magnet under harmonic base excitation shown in Figure 4. Although, previous calculations of this problem exist in the literature [34], they are either limited to slightly nonlinear oscillators and/or the dependence of the excitation amplitude on the drive frequency is not taken into account. We solve it here while addressing both issues and as a precursor to the many-oscillator case, treated in the following sections.

<figure><img src="image_5.png" /><figcaption>Figure 4: A schematic diagram of a SDOF VEH based on magnetic levitation.</figcaption></figure>

For n=1, the equation of motion of a Single Degree-of-Freedom (SDOF) VEH based on magnetic levitation can be written as follows:

$$ M_1 \ddot{v}_1 + (c + c_{e1} + c_{e2}) \dot{v}_1 + (k_{12} + k_{10}) v_1 + (\alpha_{12} - \alpha_{10}) v_1^2 + (\lambda_{12} + \lambda_{10}) v_1^3 = -M_1 \ddot{Y} \quad (17) $$

The natural frequency of the associated linear oscillator is:

$$ f_r = \frac{1}{2\pi} \sqrt{\frac{k_{12} + k_{10}}{M_1}} \quad (18) $$
Substituting Equation (12) into Equation (18), the natural frequency becomes:

$$f_r = \frac{1}{\pi} \sqrt{\frac{g(M_0 + M_2)}{2d(M_0 - M_2)}} \qquad (19)$$

where $M_0 \ge M_2$ as given by the inequality (14). Remarkably, the natural frequency depends only on the masses of the two fixed magnets. This is a consequence of assuming equal magnetic intensities and surfaces for the three magnets. Obviously, this is difficult to obtain in practice and the frequency may slightly depend on the mass of the moving magnet. Introducing the frequency $f_0 = \frac{1}{\pi}\sqrt{g/2d}$ and the dimensionless parameter $\mu = M_2/M_0$, Equation (19) becomes:

$$f_r = f_0 \sqrt{\frac{1 + \mu}{1 - \mu}} \qquad (20)$$

<figure><img src="image_6.png" /><figcaption>Figure 5: Variation of the natural frequency of a VEH based on a single moving magnet with respect to the gap $d$ and the mass ratio $\mu$.</figcaption></figure>

Figure 5 displays the variation of the natural frequency of a single moving magnet with respect to the dimensionless parameter $\mu$ and the gap $d$ separating the magnets. Notably, $f_r$ has a lower bound which is $f_0$ when $\mu = 0$, i.e. $M_2 = 0$. It corresponds to the case of a single repulsive force provided by the fixed magnet on the bottom which is indispensable to counter gravity. Moreover, Equation (19) is incompatible with the configuration $M_2 = M_0$ for which the static equilibrium position of the moving
magnet changes due to multistability and hence the resonance frequency expression will be modified.

Substituting Equation (12) into Equation (17) and dividing by $M_1$, we obtain the following equation:

$$ \ddot{v}_1 + 2\xi\omega_r\dot{v}_1 + \omega_r^2 v_1 - \frac{3g}{2d^2}v_1^2 + \frac{2}{d^2}\omega_r^2 v_1^3 = \kappa \cos(\Omega t) \quad (21) $$

where $\omega_r = 2\pi f_r$, $\xi = \frac{c + c_{e1} + c_{e2}}{2M_1\omega_r}$ and $\kappa = Y_0\Omega^2$. Equation (21) is a Duffing equation with a quadratic nonlinearity under harmonic excitation. To analyze this equation of motion, we use perturbation techniques which are well adapted to small excitation and damping, typically valid in levitated magnet-based resonators [34]. To facilitate the perturbation approach, in this case the method of averaging [38] for its ease of use, a standard constrained coordinate transformation is introduced, as given by:

$$ v_1 = A(t) \cos[\Omega t + \beta(t)] \quad (22) $$

In addition, since near-resonant behavior is the principal operating regime of the proposed system, a detuning parameter, $\sigma$, is introduced, as given by:

$$ \Omega = \omega_r + \epsilon\sigma \quad (23) $$

where $\epsilon$ is the small nondimensional bookkeeping parameter. Separating the resulting equations and averaging them over the period in the t-domain results in the system's averaged equations, in terms of amplitude and phase, which are given by:

$$ \dot{A} = \frac{1}{2\omega_r} [\kappa \sin \beta - 2\xi\omega_r^2 A] \quad (24) $$

$$ \dot{\beta} = -\sigma + \frac{1}{2A\omega_r} \left[ \frac{3\omega_r^2}{2d^2} A^3 + \kappa \cos \beta \right] \quad (25) $$

The steady-state motions occur when $\dot{A} = \dot{\beta} = 0$, which corresponds to the singular points of Equations (24) and (25). Thus, the frequency response equation can be written
in its parametric form with respect to the phase $\beta$ as follows:

$$ \sigma = \frac{3\kappa^2}{16d^2\xi^2\omega_r^3} \sin^2 \beta + \xi\omega_r \cot \beta \quad (26) $$

$$ A = \frac{\kappa}{2\xi\omega_r^2} \sin \beta \quad (27) $$

The critical amplitude is the oscillation amplitude $A_c$ above which bistability occurs [39, 40, 41]. Thus, $A_c$ is the transition amplitude from the linear to the nonlinear behavior. At the critical drive, the resonance curve exhibits a point of infinite slope, called the critical point. Moreover, at the same point, the phase curve also exhibits an infinite slope at the same detuning as the resonance curve itself. Mathematically, $A_c$ is defined as the oscillation amplitude for which the equation $\frac{d\sigma}{d\beta} = 0$ has a unique solution $\beta_c = \frac{\pi}{3}$. Thus, the critical force is deduced:

$$ \kappa_c = \frac{8\sqrt{2}}{3\sqrt[4]{3}} d\xi^{3/2} \omega_r^2 \quad (28) $$

The critical amplitude $A_c$ is obtained by substituting Equation (29) into Equation (27) at the point $\beta = \pi/2$.

$$ A_c = \frac{4\sqrt{2}}{3\sqrt[4]{3}} d\sqrt{\xi} \quad (29) $$

<figure><img src="image_7.png" /><figcaption>Figure 6: Forced frequency responses of the SDOF VEH with design parameters listed in Table 1 for three configurations: linear ($Y_0 = 0.7 \text{ mm}$), critical ($Y_0 = 2 \text{ mm}$), nonlinear ($Y_0 = 4.8 \text{ mm}$). Solid lines denote stable periodic solutions and dashed lines represent unstable periodic solutions.</figcaption></figure>
<table><caption>Table 1: Design parameters of the considered SDOF VEH.</caption><thead><tr><th>Parameter</th><th>M<sub>0</sub><br>(g)</th><th>M<sub>1</sub><br>(g)</th><th>M<sub>2</sub><br>(g)</th><th>c<sub>e1</sub><br>(Ns/m)</th><th>c<sub>e2</sub><br>(Ns/m)</th><th>c<br>(Ns/m)</th><th>d<br>(mm)</th></tr></thead><tbody><tr><td>Value</td><td>80</td><td>80</td><td>10</td><td>0.303</td><td>0.303</td><td>0.119</td><td>15</td></tr></tbody></table>

Figure 6 shows the evolution of the harvester frequency response (Equations (26) and (27)) with respect to the excitation amplitude for design parameters listed in Table 1. These parameters have been chosen in order to obtain a SDOF VEH resonating at a frequency below 10 Hz ($f_r = 6.5$ Hz) with a low electromagnetic damping. Beyond the critical amplitude, the bandwidth definition is reconsidered and determined with respect to the highest bifurcation point. Notably, the harvester bandwidth is significantly enlarged when the moving magnet is driven in the nonlinear regime.

## 3.2. Multi-objective optimization of a 2-DOFs VEH

For n=2, Equation (11) leads to the following coupled equations of motion:

$$ 
\left\{ 
\begin{array}{l}
\left( 
\begin{array}{l}
M_1 \ddot{v}_1 + (c + c_{e1} + c_{e2}) \dot{v}_1 - c_{e2} \dot{v}_2 \\
\quad + k_{12} (v_1 - v_2) + k_{10} v_1 + \alpha_{12} (v_1 - v_2)^3 \\
\quad - \alpha_{10} v_1^2 + \lambda_{12} (v_1 - v_2)^3 + \lambda_{10} v_1^3
\end{array}
\right) = -M_1 \ddot{Y} \\[2ex]
\left( 
\begin{array}{l}
M_2 \ddot{v}_2 + (c + c_{e2} + c_{e3}) \dot{v}_2 - c_{e2} \dot{v}_1 \\
\quad + k_{23} v_2 + k_{21} (v_2 - v_1) + \alpha_{23} v_2^2 \\
\quad - \alpha_{21} (v_2 - v_1)^2 + \lambda_{23} v_2^3 + \lambda_{21} (v_2 - v_1)^3
\end{array}
\right) = -M_2 \ddot{Y}
\end{array}
\right.
\quad (30)
$$

In this case, we seek to take advantage of the nonlinear magnetic coupling between the oscillators in order to enhance the bandwidth and the harvested power of the considered device. To do so, two objective functions $g_1$ and $g_2$ have been defined: the first one is related to the distance between the two natural frequencies of the system $f_1$ and $f_2$ while the second one concerns the intermodal distance between the harvested powers
<figure><img src="image_8.png" /><figcaption>Figure 7: Application of the adapted Morris method to the objective functions $g_1$ and $g_2$.</figcaption></figure>

taken at the eigenfrequencies.

$$g_1 = \frac{f_2 - f_1}{f_1} \qquad (31)$$

$$g_2 = \frac{|P_1 - P_2|}{P_{n1}} \qquad (32)$$

Where $P_1$ and $P_{n1}$ are the harvested and normalized powers at the frequency $f_1$ (Equation (16) for $\Omega = \Omega_1 = 2\pi f_1$) and $P_2$ is the harvested power at the frequency $f_2$ (Equation (16) for $\Omega = \Omega_2 = 2\pi f_2$).

Since $g_1$ and $g_2$ depend in several design parameters, as a first step, a Sensitivity Analysis (SA) has been performed in order to select unessential parameters that can be held constant. Among several SA methods, the Morris method is a specialized randomized one-factor-at-time SA design [42] and its efficiency was improved by using a Latin hypercube sampling instead of a random sampling and radial points instead of trajectories [43]. This method was used to determine the boundary between high and low influence factors among the parameters $c_{e1}$, $c_{e2}$, $c_{e3}$, $M_1$, $M_2$, $M_3$ and $d_0$. This limit is defined as the first variation is equal or greater to $d_{max}/\gamma$, where $d_{max}$ is the largest distance between two consecutive maximum factors and $\gamma = 10$.
Figure 7 shows that the first objective function $g_1$ depends only on the magnet masses and the parameters $c_{e1}$ and $c_{e2}$ have a low influence on the variation of the second objective function $g_2$. Consequently, $c_{e1}$ and $c_{e2}$ are held constant in the optimization procedure discussed below.

In order to improve the performances of the two degrees-of-freedom nonlinear VEH, the resonance peaks must be very close to each other in terms of natural frequencies and modal harvested power amplitudes. To do so, the objective functions $g_1$ and $g_2$ are used in an optimization procedure (NSGA-II) [37]. The two cost functions are simultaneously minimized with respect to the five retained design parameters ($M_1, M_2, M_3, d, c_{e3}$).

$$ 
\begin{array}{l}
\text{Minimize} \quad [g_1(x), g_2(x)], \\
\\
x = [M_1, M_2, M_3, d, c_{e3}]^T, \\
\\
\text{subject to} \quad M_3 + M_0 = M_1 + M_2, \\
\\
M_0 \geq M_2, M_1 \geq M_3
\end{array}
\qquad (33)
$$

<figure><img src="image_9.png" /><figcaption>Figure 8: Preferred solutions of the multi-objective optimization problem written in Equation (33).</figcaption></figure>

The sensitivity analysis was followed by a global optimization in order to localize local and global optima for each objective function separately using Monte Carlo simulations in the design space. Then, the multi-objective optimization was performed
<figure><img src="image_10.png" /><figcaption>Figure 9: Forced frequency responses of a 2-DOFs VEH in terms of (a) amplitude and phase (b) harvested power. The design parameters correspond to those of solution 1 in Table 2 and the excitation amplitude is $Y_0 = 2 \text{ mm}$.</figcaption></figure>

around the global optima. The Pareto-optimal solutions obtained by NSGA-II after 200 generations are shown in Figure 8. The Pareto optimal front is convex and its density and range depend on the design parameter constraints. Three solutions have been chosen on the Pareto front and their corresponding design parameters are listed in Table 2.

<table><caption>Table 2: Design parameters of the three selected optimal solutions in Figure 8.</caption><thead><tr><th></th><th>Solution 1</th><th>Solution 2</th><th>Solution 3</th><th>Range of variation</th></tr></thead><tbody><tr><td>$g_1(cm^3g^2)$</td><td>2.5 * 10<sup>-7</sup></td><td>3.6 * 10<sup>-3</sup></td><td>5.2 * 10<sup>-3</sup></td><td>[ ]</td></tr><tr><td>$g_2$</td><td>0.72</td><td>0.63</td><td>0.62</td><td></td></tr><tr><td>$M_1$ (g)</td><td>80</td><td>80</td><td>80</td><td>[10-250]</td></tr><tr><td>$M_2$ (g)</td><td>32</td><td>31.6</td><td>32</td><td>[10-250]</td></tr><tr><td>$M_3$ (g)</td><td>10</td><td>26.3</td><td>32</td><td>[10-250]</td></tr><tr><td>$c_{e3}$ (Ns/m)</td><td>0.27</td><td>0.71</td><td>1.1</td><td>[0.1-5]</td></tr><tr><td>$d_0$ (mm)</td><td>15</td><td>15</td><td>15</td><td>[5-15]</td></tr><tr><td>$c_{e1}$ (Ns/m)</td><td>1.19</td><td>1.19</td><td>1.19</td><td>-</td></tr><tr><td>$c_{e2}$ (Ns/m)</td><td>0.13</td><td>0.13</td><td>0.13</td><td>-</td></tr><tr><td>$c$ (Ns/m)</td><td>0.119</td><td>0.119</td><td>0.119</td><td>-</td></tr></tbody></table>
The frequency responses in terms of amplitude, phase and harvested power are displayed in Figure 9 for design parameters of solution 1 and an excitation amplitude $Y_0 = 2$ mm. It is shown that the two moving magnets vibrate in phase up to an excitation frequency of 5.7 Hz. Then, a phase between the two signals rises and reaches a maximum of $\pi/3$ for an excitation frequency equal to 9.5 Hz. Unlike the case of a single moving magnet, the bandwidth of a multi-degree of freedom VEH is measured on the power frequency response between an attenuation of $-3\,\text{dB}$ on the first peak and the horizontal slope preceding the last resonance peak as shown in Figure 9(b). Remarkably, the bandwidth is significantly high; it reaches 110% with respect to the first natural frequency.

<figure><img src="image_11.png" /><figcaption>Figure 10: Forced frequency responses of the three optimal solutions for a 2-DOFs VEH in term of normalized harvested power.</figcaption></figure>

The frequency responses of the three selected solutions are compared in term of normalized harvested power as shown in Figure 10. The first retained solution has the highest normalized power close to $12 \, \text{mW} \, \text{cm}^{-3} \, \text{g}^{-2}$. The specifications of the three design solutions are given in Table 3, approving the high performances of the first optimal solution.
<table><caption>Table 3: Performances of 2-DOFs VEHs for the three selected optimal solutions in Figure 11.</caption><thead><tr><th>Solution</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>f<sub>1</sub>(Hz)</td><td>5.2</td><td>6.1</td><td>6.5</td></tr><tr><td>f<sub>2</sub>(Hz)</td><td>8.9</td><td>10</td><td>10.6</td></tr><tr><td>P<sub>n</sub>(mWcm<sup>-3</sup>g<sup>-2</sup>)</td><td>10.4</td><td>8.8</td><td>9.2</td></tr><tr><td>BW(%)</td><td>120</td><td>85</td><td>84</td></tr></tbody></table>

### 3.3. Multi-objective optimization of a 3-DOFs VEH

Similarly to the case of two moving magnets, an optimization procedure is performed using NSGA-II for $n = 3$. It involves four objective functions: two of them concern the distance between the eigenfrequencies while the two others represent the distance between the resonance peaks in term of harvested power.

$$g_1 = \frac{|P_3 - P_1|}{g_1^* P_{n1}} \qquad (34)$$

$$g_2 = \frac{|P_2 - P_1|}{g_2^* P_{n1}} \qquad (35)$$

$$g_3 = \frac{f_3 - f_1}{g_3^* f_1} \qquad (36)$$

$$g_4 = \frac{f_2 - f_1}{g_4^* f_1} \qquad (37)$$

where $g_i^*$ is the extrema of the function $g_i$. The adapted Morris method was used to eliminate parameters with a low influence on the objective functions leading to four retained design parameters which are $M_2, M_3, M_4$ and $d$. The considered optimization
<figure><img src="image_12.png" /><figcaption>Figure 11: Preferred solutions of the multi-objective optimization problem written in Equation (38).</figcaption></figure>

problem is written as follows:

$$
\begin{array}{ll}
\text{Minimize} & [g_1(x), g_2(x), g_3(x), g_4(x)], \\
& x = [M_2, M_3, M_4, d]^T, \\
\text{subject to} & M_1 - M_3 = M_2 - M_4, \\
& M_0 - M_2 = M_1 - M_2, \\
& M_1 - M_3 = M_2 - M_4, \\
& M_0 \geq M_2, M_1 \geq M_3, M_2 \geq M_4
\end{array}
\qquad (38)
$$

The Preferred solutions of the obtained distribution by 25 population size for 1000
generation are shown in Figure 11(a). Three solutions displayed in Figure 11(b) have
been chosen and their corresponding design parameters are listed in Table 4. The
frequency responses in terms of amplitude, phase and harvested power are displayed in
Figure 12 for design parameters of solution 1 and an excitation amplitude Y₀ = 2.3 mm.
It is shown that the two moving magnets M₂ and M₃ vibrate either in phase or with a
low phase shift in the frequency range [0 20]Hz. Remarkably, the displacement of M₁
has a low phase shift with respect to M₂ and M₃ displacements, up to an excitation
frequency of 6 Hz and reaches a maximum of 2π/5 for an excitation frequency equal to
11 Hz. Finally, M₁ and M₃ vibrate in phase for f = 15 Hz.
<table><caption>Table 4: Design parameters of the three selected optimal solutions in Figure 11.</caption><thead><tr><th></th><th>Solution 1</th><th>Solution 2</th><th>Solution 3</th><th>Range of variation</th></tr></thead><tbody><tr><td>g<sub>1</sub></td><td>0.27</td><td>0.48</td><td>0.62</td><td>[0 - 1]</td></tr><tr><td>g<sub>2</sub></td><td>0.01</td><td>1</td><td>0.05</td><td>[0 - 1]</td></tr><tr><td>g<sub>3</sub></td><td>0.86</td><td>0.85</td><td>0.85</td><td>[0 - 1]</td></tr><tr><td>g<sub>4</sub></td><td>0.63</td><td>0.92</td><td>0.92</td><td>[0 - 1]</td></tr><tr><td>M<sub>2</sub> (g)</td><td>90</td><td>110</td><td>128</td><td>[10 - 250]</td></tr><tr><td>M<sub>3</sub> (g)</td><td>10</td><td>67</td><td>63</td><td>[10 - 250]</td></tr><tr><td>M<sub>4</sub> (g)</td><td>10</td><td>43</td><td>50</td><td>[10 - 250]</td></tr><tr><td>d<sub>0</sub> (mm)</td><td>15</td><td>15</td><td>15</td><td>[5 - 15]</td></tr><tr><td>c<sub>e1</sub> (Ns/m)</td><td>1.02</td><td>1.02</td><td>1.02</td><td>-</td></tr><tr><td>c<sub>e2</sub> (Ns/m)</td><td>0.76</td><td>0.76</td><td>0.76</td><td>-</td></tr><tr><td>c<sub>e3</sub> (Ns/m)</td><td>1.24</td><td>1.24</td><td>1.24</td><td>-</td></tr><tr><td>c<sub>e4</sub> (Ns/m)</td><td>1.7</td><td>1.7</td><td>1.7</td><td>-</td></tr><tr><td>c (Ns/m)</td><td>0.119</td><td>0.119</td><td>0.119</td><td>-</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Figure 12: Forced frequency responses of a 3-DOFs VEH in terms of (a) amplitude and phase (b) harvested power. The design parameters correspond to those of solution 1 in Table 4 and the excitation amplitude is $Y_0 = 2.3 \text{ mm}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Forced frequency responses of the three optimal solutions for a 3-DOFs VEH in terms of normalized harvested power.</figcaption></figure>
The frequency responses of the three selected solutions are compared in term of normalized harvested power as shown in Figure 13. The first retained solution has the highest normalized power close to $25 mWcm^{-3}g^{-2}$. The specifications of the three design solutions are given in Table 5, approving the high performances of the first optimal solution.

<table><caption>Table 5: Performances of 3-DOFs VEHs for the three selected optimal solutions in Figure 11.</caption><thead><tr><th></th><th>Solution 1</th><th>Solution 2</th><th>Solution 3</th></tr></thead><tbody><tr><td>f<sub>1</sub>(Hz)</td><td>4.3</td><td>5.4</td><td>5.3</td></tr><tr><td>f<sub>2</sub>(Hz)</td><td>6.8</td><td>9.9</td><td>9.7</td></tr><tr><td>f<sub>3</sub>(Hz)</td><td>11</td><td>13.8</td><td>13.5</td></tr><tr><td>P<sub>n</sub>(mWcm<sup>-3</sup>g<sup>-2</sup>)</td><td>20.22</td><td>19.42</td><td>18.68</td></tr><tr><td>BW(%)</td><td>190</td><td>168</td><td>173</td></tr></tbody></table>

### 3.4. Comparative study

In this section, comparisons between several magnetic levitation based-VEHs will be presented in order to highlight the advantages of devices having an array of coupled nonlinear oscillators in terms of bandwidth, harvested power and operating frequency range.

#### 3.4.1. Nonlinear coupling benefits
The performances of two SDOF VEHs denoted $VEH_1^1$ and $VEH_2^1$ are compared to those of a 2-DOFs VEH denoted $VEH^2$. The design parameters of the 2-DOFs VEH correspond to those of the first optimal solution listed in Table 2. A significant comparison is obtained if the design parameters of $VEH_1^1$ and $VEH_2^1$ are chosen in such a way that their natural frequencies match perfectly those of $VEH^2$. The latter has its first natural frequency equal to 5.2/ Hz which is higher than the lower bound limit of the natural frequency of a SDOF VEH given by $f_0 = \frac{1}{\pi}\sqrt{g/2d}$ (see section 3.1) if we maintain the same gap value. Consequently, the gap between
<figure><img src="image_15.png" /><figcaption>Figure 14: Forced frequency responses of two SDOF VEHs in terms of (a) amplitude and phase (b) overall harvested power. The design parameters are listed in Table 6 and the excitation amplitude is $Y_0 = 3$ mm.</figcaption></figure>

The magnets is chosen $d = 24$ mm for $VEH_1^1$ and $VEH_2^1$ while the values of the other design parameters are taken from those of $VEH^2$ except the top magnet masses which are adjusted to ensure a perfect matching between the natural frequencies (Table 6).

<table><caption>Table 6: Design parameters of $VEH_1^1$ and $VEH_2^1$.</caption><thead><tr><th>Parameter</th><th>$M_0$<br>(g)</th><th>$M_1$<br>(g)</th><th>$M_2$<br>(g)</th><th>$c_{e1}$<br>(Ns/m)</th><th>$c_{e2}$<br>(Ns/m)</th><th>$c$<br>(Ns/m)</th><th>$d$<br>(mm)</th></tr></thead><tbody><tr><td>$VEH_1^1$</td><td>80</td><td>80</td><td>10</td><td>1.19</td><td>0.13</td><td>0.119</td><td>24</td></tr><tr><td>$VEH_2^1$</td><td>80</td><td>32</td><td>46</td><td>1.19</td><td>0.27</td><td>0.119</td><td>24</td></tr></tbody></table>

The frequency responses of $VEH_1^1$ and $VEH_2^1$ in terms of amplitude, phase and harvested power are displayed in Figure 14 for design parameters listed in Table 6 and an excitation amplitude $Y_0 = 3$ mm. It is shown that the dynamic behavior of $VEH_1^1$ is linear while $VEH_2^1$ is characterized by a hardening frequency response below the critical amplitude. The maximum of phase between the two VEHs exceeds $\pi/2$ and the overall normalized harvested power of the two SDOF VEHs is about $6.4 \text{ mW cm}^{-3} \text{g}^{-2}$ which is lower than the one of the $VEH^2$ as shown in Figure 15.
<figure><img src="image_16.png" /><figcaption>Figure 15: Forced frequency responses of $VEH_1^1 + VEH_2^1$ and $VEH_2^2$ in term of normalized harvested power.</figcaption></figure>

The performances of $VEH_2^2$ are compared to the overall performances of $VEH_1^1$ and $VEH_1^2$ in Table 7, which emphasize the advantages of the nonlinear coupling between the moving magnets. Indeed, the magnetic coupling induces nonlinearities leading to an enhancement of the frequency bandwidth via a nonlinear spring hardening effect. Moreover, unlike the case of two SDOF harvesters, the magnetic parameters can be tuned to modify the coupling terms in such a way that the natural frequencies of the 2-DOFs VEH become closer and consequently the modal interactions increase leading to high vibration amplitudes, i.e. the harvested power increases.

<table><caption>Table 7: A comparison of performances between $VEH_1^1 + VEH_2^1$ and $VEH_2^2$.</caption><thead><tr><td>Parameter</td><td>Frequency (Hz)</td><td>P<sub>n</sub>(mWcm<sup>-3</sup>g<sup>-2</sup>)</td><td>BW(%)</td></tr></thead><tbody><tr><td>VEH<sub>1</sub><sup>1</sup> + VEH<sub>2</sub><sup>1</sup></td><td>5.3-12.6</td><td>6.4</td><td>115</td></tr><tr><td>VEH<sub>2</sub><sup>2</sup></td><td>5.1-12.1</td><td>10.4</td><td>120</td></tr></tbody></table>

### 3.4.2. Modal interaction benefits
The frequency responses of $VEH_2^2$ and a 3-DOFs VEH denoted $VEH_3^3$ are compared in term of normalized harvested power as shown in Figure 16. The design parameters of the 3-DOFs VEH correspond to those of the first optimal solution listed in Table 4. The 3-DOFs VEH has the highest normalized power close to $25 \text{ mWcm}^{-3}\text{g}^{-2}$. The performances of both devices are given in Table
8, approving the high performances of $VEH^3$ compared to $VEH^2$ and giving credits to the modal interactions and their benefits in such systems. The remarkable specifications of the proposed devices based on arrays of levitated magnets are highlighted in Table 9 with respect to existing magnetic harvesters.

<figure><img src="image_17.png" /><figcaption>Figure 16: Forced frequency responses of $VEH^2$ and $VEH^3$ in term of normalized harvested power.</figcaption></figure>

<table><caption>Table 8: A comparison of performances between $VEH^2$ and $VEH^3$.</caption><thead><tr><th>Parameter</th><th>Frequency (Hz)</th><th>P<sub>n</sub>(mWcm<sup>-3</sup>g<sup>-2</sup>)</th><th>BW(%)</th></tr></thead><tbody><tr><td>VEH<sup>2</sup></td><td>5.1 - 12.6</td><td>12.1</td><td>120</td></tr><tr><td>VEH<sup>3</sup></td><td>4.6 - 14.5</td><td>20.2</td><td>190</td></tr></tbody></table>

<table><caption>Table 9: A comparison of performances between several energy harvesters based on magnetic levitation and the proposed VEHs.</caption><thead><tr><th>Refs</th><th>Frequency (Hz)</th><th>Acceleration (g)</th><th>Power (mW)</th><th>Normalized Power (mWcm<sup>-3</sup>g<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Abu Riduan et al. [10]</td><td>7-10</td><td>0.5</td><td>2.090</td><td>0.1</td></tr><tr><td>Daniel et al. [44]</td><td>13</td><td>0.25</td><td>5.9</td><td>15.33</td></tr><tr><td>Marin et al. [45]</td><td>47.75-51.75</td><td>0.2</td><td>25.5</td><td>0.54</td></tr><tr><td>Cepnik et al. [46]</td><td>50</td><td>1</td><td>20.6</td><td>1.01</td></tr><tr><td>Byung-chul et al. [47]</td><td>16</td><td>0.2</td><td>1.52</td><td>1.07</td></tr><tr><td>VEH<sup>2</sup></td><td>5.1-12</td><td>0.22</td><td>335</td><td>10.4</td></tr><tr><td>VEH<sup>3</sup></td><td>4.6-14.5</td><td>0.17</td><td>795</td><td>20.2</td></tr></tbody></table>
4. Conclusion

In this paper, we propose a multi-modal vibration energy harvesting method based
on arrays of coupled levitated magnets. The equations of motion have been derived
taking into account the magnetic nonlinearity and electromagnetic damping. They have
been solved using the harmonic balance method coupled with the asymptotic numerical
method while considering frequency dependent excitation amplitude. The case of a
SDOF VEH has been analytically investigated and close form expression of the critical
amplitude was extracted. Multi-objective optimization procedures are introduced and
performed using NSGA-II algorithm for the cases of two and three moving magnets
in order to select the ideal solutions in term of performances. The advantages of
the nonlinear coupling have been pointed out by comparing the joint performances
of two uncoupled VEHs and those of a 2-DOFs VEH. The latter enables harvesting
the vibration energy in the frequency range 5.1 – 12 Hz with a bandwidth of 120 %
and a normalized harvested power of 10.4 mWcm⁻³g⁻². Moreover, a 3-DOFs VEH
proves to be even more performant with an operating frequency range of 4.6 – 14.5 Hz,
a bandwidth of 190 % and a normalized harvested power of 20.2 mWcm⁻³g⁻² which
highlights the benefits of modal interactions within the proposed harvesting energy
method.

With these promising results, the work in progress concerns the modeling of
multimodal interactions for large oscillator arrays. Indeed, when the number of coupled
magnets increases, the number of excited modes increases over a wide frequency range
leading to ultra-wide bandwidth devices and a high harvested power. Finally, the
nonlinearity can be functionalized in such periodic structures in order to generate
particular collective dynamics [48] with a large number of stable multi-mode solution
branches suitable for energy harvesting.
Finally, the experimental validation of the proposed multi-modal energy harvesting method is an important task that we will attempt to accomplish in the near future. Even though the implementation of an experimental protocol is expensive and takes time, it remains one of our principal objectives.

## Acknowledgments

This project has been performed in cooperation with the Labex ACTION program (contract ANR-11-LABX-01-01).

## Appendix

Equations (11) and (16) can be written in the following form:

$$ \dot{w} = f(t, w, \Omega) \qquad (39) $$

Where $w$ is a vector of unknowns, $f$ is periodic in $t$ and $(\Omega)$ is a real parameter.

The key point of this method lies in the quadratic recast of Equation (39) by introducing the following set of auxiliary variables:

$$ \begin{aligned} u_j &= \dot{v}_j, w_j = v_j^2; \forall j \in [1, n-2] \\ w_{j,j+1} &= v_j v_{j+1}; \forall j \in [1, n-3] \\ r_1 &= \Omega, r_2 = Y_0 \cos(\Omega t), r_3 = \Omega^2 \end{aligned} \qquad (40) $$

Then, Equation (39) can be written in matrix form as follows :
$$
\begin{bmatrix} \dot{v}_j \\ \dot{u}_j \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ Y_0 \cos(\Omega t) \\ 0 \\ 0 \end{bmatrix} + \begin{bmatrix} u_j \\ Ul_j \\ -w_j \\ -w_{j,j+1} \\ -r_1 \\ -r_2 \\ -r_3 \\ P \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \Omega \\ 0 \\ 0 \\ 0 \end{bmatrix} + \begin{bmatrix} 0 \\ Uq_j \\ v_j^2 \\ v_j v_{j+1} \\ 0 \\ 0 \\ r_1^2 \\ Pq \end{bmatrix} \\
m(\dot{Z}) = c(t, \Omega) + l(Z) + \Omega l_{c1} + q(Z, Z) \tag{41}
$$

with :

$$
\begin{align}
U l_j &= \begin{pmatrix}
-\frac{c+c_{ej-1}+c_{ej}}{M_j} u_j + \frac{c_{ej}}{M_j} u_{j+1} + \frac{c_{ej-1}}{M_j} u_{j-1} \\
-\frac{k_{j,j+1}}{M_j} (v_j - v_{j+1}) - \frac{k_{j,j-1}}{M_j} (v_j - v_{j-1}) \\
-\frac{\alpha_{j,j+1}}{M_j} (w_j + w_{j+1} - 2w_{j,j+1}) \\
+\frac{\alpha_{j,j-1}}{M_j} (w_j + w_{j-1} - 2w_{j,j-1})
\end{pmatrix} \tag{42} \\
U q_j &= \begin{pmatrix}
-\frac{\lambda_{j,j+1}}{M_j} (v_j w_j - v_{j+1} w_{j+1} + 3v_j w_{j+1} - 3w_j v_{j+1}) \\
-\frac{\lambda_{j,j-1}}{M_j} (v_j w_j - v_{j-1} w_{j-1} + 3v_j w_{j-1} - 3w_j v_{j-1}) + r_3 r_2
\end{pmatrix} \\
P q &= -\frac{r_3}{2} \sum_{j=1}^{n-1} (c_{ej} (w_j + w_{j+1} - 2w_{j,j+1}))
\end{align}
$$

where:

$Z = [v_1...v_n \ u_1...u_n \ w_1...w_n \ w_{12}...w_{n-1,n} \ r_1 \ r_2 \ r_3 \ P]^T$ is the vector of unknowns with size $N_{eq} = 4(n-2)+3$, $c(t, \Omega)$ and $l_{c1}$ are constant vectors, $m(\dot{Z})$ and $l(Z)$ are linear vectors and $q(Z, Z)$ is a quadratic vector.

The harmonic balance method is applied to Equations (41) to decompose the
solution $Z(t)$ into a truncated Fourier series.

$$
Z(t) = Z_0 + \sum_{k=1}^{H} Z_{c,k} \cos(k\Omega t) + \sum_{k=1}^{H} Z_{s,k} \sin(k\Omega t) \quad (43)
$$
The column vector $U$, with size $(2H+1)Neq$, where $Neq$ is the number of equations in Equation (41), collects the components of the Fourier series as:

$$U = [Z_0^T + Z_{c,1}^T + Z_{s,1}^T + Z_{c,2}^T + Z_{s,2}^T + \dots + Z_{c,H}^T + Z_{s,H}^T]^T \quad (44)$$

Substituting Equation (44) into Equation (41), assembling the terms of the same harmonic and neglecting the highest order, we obtain the following system:

$$\Omega M(U) = C + L(U) + Q(U, U) \quad (45)$$

where $M(\cdot), C, L(\cdot)$ and $Q(\cdot, \cdot)$ are operators which depend on $m(\cdot), c, l(\cdot)$ and $q(\cdot, \cdot)$ as it is defined by Cochelin et al. [36]. Once the algebraic system is obtained, we solve it using a continuation technique. In order to apply the ANM, equation (45) is transformed into:

$$R(U, \omega) = C + L(U) + Q(U, U) - \omega M(U) = 0 \quad (46)$$

## References

[1] Williams C and Yates R 1996 *Sensors and Actuators A: Physical* **52** 8–11

[2] Donelan J M, Li Q, Naing V, Hoffer J A, Weber D J and Kuo A D 2008 *Science* **319** 807–810

[3] Ding Z, Perlaza S, Esnaola I and Poor H 2014 *Wireless Communications, IEEE Transactions on* **13** 846–860

[4] Wang W, Wang N, Jafer E, Hayes M, O’Flynn B and O’Mathuna C 2010 *Environmental Science and Information Application Technology (ESIAT), 2010 International Conference on* **3** 367–372

[5] Ottman G, Hofmann H, Bhatt A and Lesieutre G 2002 *Power Electronics, IEEE Transactions on* **17** 669–676

[6] Roundy S and Wright P K 2004 *Smart Materials and Structures* **13** 1131

[7] Dutoit N E, Wardle B L and Kim S G 2005 *Integrated Ferroelectrics* **71** 121–160

[8] Sodano H A, Inman D J and Park G 2005 *Journal of Intelligent Material Systems and Structures* **16** 799–807

[9] Stephen N 2006 *Journal of Sound and Vibration* **293** 409–425

[10] Foisal A R M, Hong C and Chung G S 2012 *Sensors and Actuators A: Physical* **182** 106–113
[11] Glynne-Jones P, Tudor M, Beeby S and White N 2004 *Sensors and Actuators A: Physical* **110** 344–349

[12] Beeby S P, Torah R N, Tudor M J, Glynne-Jones P, O'Donnell T, Saha C R and Roy S 2007
*Journal of Micromechanics and Microengineering* **17** 1257

[13] Leland E S and Wright P K 2006 *Smart Materials and Structures* **15** 1413

[14] EichhornC G and WoiasP 2008 *Proc. PowerMEMS (Sendai)* 309312

[15] Hu Y, Xue H and Hu H 2007 *Smart Materials and Structures* **16** 1961

[16] Loverich J, Geiger R and Frank J 2008 *Proc. SPIE* **6928** 692805

[17] Challa V R, Prasad M G, Shi Y and Fisher F T 2008 *Smart Materials and Structures* **17** 015035

[18] Reissman T, Wolff E M and Garcia E 2009 *Proc. SPIE* **7288** 72880G

[19] Gu L and Livermore C 2010 *Applied Physics Letters* **97** 081904

[20] Youngsman J M, Luedeman T, Morris D J, Anderson M J and Bahr D F 2010 *Journal of Sound and Vibration* **329** 277 – 288

[21] Wu X, Lin J, Kato S, Zhang K, Ren T and Liu L 2008 *Proc. PowerMEMS* 245–248

[22] ZhuD, RobertsS, TudorMJ and BeebyS 2008 *Proc. PowerMEMS* 229–232

[23] Roundy S, Leland E, Baker J, Carleton E, Reilly E, Lai E, Otis B, Rabaey J, Wright P and Sundararajan V 2005 *Pervasive Comput. IEEE* **4** 28–36

[24] Yang B, Lee C, Xiang W, Xie J, He J H, Kotlanka R K, Low S P and Feng H 2009 *Journal of Micromechanics and Microengineering* **19** 035001

[25] Sari I, Balkan T and Kulah H 2008 *Sensors and Actuators A: Physical* **145146** 405 – 413

[26] Yang Z and Yang J 2009 *Journal of Intelligent Material Systems and Structures* **20** 569–574

[27] Erturk A, Renno J M and Inman D J 2008 *Journal of Intelligent Material Systems and Structures*
**20** 1–16

[28] Tadesse Y, Zhang S and Priya S 2009 *Journal of Intelligent Material Systems and Structures* **20** 625–632

[29] Arafa M, Akl W, Aladwani A, Aldraihem O and Baz A 2011 *Proc. SPIE* **7977** 79770Q

[30] Kim I H, Jung H J, Lee B M and Jang S J 2011 *Applied Physics Letters* **98** 214102

[31] Tang L, Yang Y and Soh C 2013 Broadband vibration energy harvesting techniques *Advances in Energy Harvesting Methods* ed Elvin N and Erturk A (Springer New York) pp 17–61 ISBN 978-1-4614-5704-6

[32] Daqaq M F, Stabler C, Qaroush Y and Seuaciuc-Osrio T 2009 *Journal of Intelligent Material*
[33] Masana R and Daqaq M F 2011 *Journal of Vibration and Acoustics: Transactions of the ASME* **16** 60366052

[34] Mann B and Sims N 2009 *Journal of Sound and Vibration* **319** 515 – 530

[35] Mahmoudi S, Kacem N and Bouhaddi N 2014 *Smart Materials and Structures* **23** 075024

[36] Cochelin B and Vergez C 2009 *Journal of Sound and Vibration* **324** 243 – 262

[37] Deb K, Pratap A, Agarwal S and Meyarivan T 2002 *Evolutionary Computation, IEEE Transactions on* **6** 182–197

[38] Nayfeh A H 1981 *Introduction to Perturbation Techniques* (Wiley, New York, USA)

[39] Kacem N, Hentz S, Pinto D, Reig B and Nguyen V 2009 *Nanotechnology* **20** 275501

[40] Kacem N, Arcamone J, Perez-Murano F and Hentz S 2010 *Journal of Micromechanics and Microengineering* **20** 045023

[41] Kacem N and Hentz S 2009 *Applied Physics Letters* **95** 183104

[42] Morris M D 1991 *Technometrics* **33** pp. 161–174

[43] Sohier H, Piet-Lahanier H and Farges J L 2015 *Acta Astronautica* **108** 18 – 29

[44] Apo D J and Priya S 2014 *Energy Harvesting and Systems* **1** 1 – 2

[45] Marin A, Turner J, Ha D S and Priya S 2013 *Smart Materials and Structures* **22** 075008

[46] Cepnik C, Radler O, Rosenbaum S, Strhla T and Wallrabe U 2011 *Sensors and Actuators A: Physical* **167** 416 – 421

[47] Lee B C, Rahman M A, Hyun S H and Chung G S 2012 *Smart Materials and Structures* **21** 125024

[48] Bitar D, Kacem N, Bouhaddi N and Collet M 2015 *Nonlinear Dynamics* **82** 749–766