On the Optimization of a Multimodal Electromagnetic Vibration Energy Harvester Using Mode Localization and Nonlinear Dynamics

: In this paper we study a generic model of a nonlinear quasiperiodic vibration energy harvester (VEH) based on electromagnetic transduction. The proposed device consists of multiple moving magnets guided by elastic beams and coupled by repulsive magnetic forces. A system of two degrees-of-freedom (DOFs) with tunable nonlinearity and mode localization is experimentally validated. The validated 2-DOFs harvester is optimized using a multiobjective optimization procedure to improve its harvested power and frequency bandwidth. An efﬁcient criterion using the modal kinetic energy of the ﬁnite element model is proposed to quantify the energy localized in the structure perturbed zones. Afterward, this concept has been generalized to a 5-DOFs VEH with two perturbed DOFs oscillators and the optimal performances are derived using a multiobjective optimization. This proposed model enables a signiﬁcant increase in the harvested power and frequency bandwidth by 101% and 79%, respectively, compared to that of the 2-DOFs device. Moreover, it has been shown that harvesting energy from two perturbed magnets among ﬁve provides almost the same amount of harvested energy and enhances the frequency bandwidth by 18% compared to those of the periodic system. Consequently, the harvester can be improved by reducing the transduction circuits number and the manufacturing cost.


Introduction
Over the last years, wearable devices and embedded systems have been increasingly introduced in diverse applications. Most portable devices depend on batteries as the main energy source. However, they have a limited lifespan even with continuous progress in producing high-energy-density batteries [1]. Therefore, collecting energy from the environment ambient sources and converting it into usable electricity is challenging [2][3][4]. Harvesting energy techniques from different sources have received much attention [5]. Energy harvesting from vibrations, as one of the most studied methods [6,7], has a major challenge in converting vibration energy into electrical energy. Multiple studies proposed different designs of vibration energy harvesters [8,9] used in wireless sensors which are considered as self powered systems. The energy storage is possible through different conversion technologies namely electrostriction [10], magnetostriction [11], piezoelectricity [12,13] and electromagnetism [14].
In spite of considerable progress, major limitations need to be overcome in the field of vibration energy harvesting [15]. For instance, most vibration energy harvesting devices effectively operate in a narrow bandwidth near to their resonance frequency. Consequently, their application is limited to specific domains and they can not be used where energy prevails over a larger bandwidth. To overcome this problem, several approaches concerning the vibration energy harvesters have been proposed and have proved reliable results. Among them, one can mention the introduction of the nonlinearity [16,17], the adoption of multimodal configurations [18,19], the combination of two or more of the techniques mentioned, etc.
Adopting a multimodal approach helps to cover a certain range of frequency to achieve a broader bandwidth. For that reason, multiple researchers have developed energy harvesters with multiple DOFs [20][21][22][23]. Benefits of the multimodal method with the functionalization of the energy localization phenomenon have been, recently, investigated [24,25]. The energy localization phenomenon, first developed by Anderson [26], is exploited to increase the amplitude of vibration and also to enhance the harvested power. When an irregularity is introduced in the periodic system, the energy will be localized in regions near to this imperfection instead of being propagated in an equal manner to all system's regions [27][28][29]. Although the fact that multimodal techniques enable a wider bandwidth for energy harvesters, their implementation leads to severe technological constraints due to the interface required to the electrical circuits and the high costs of the latter. Consequently, the nonlinearity is introduced to address this limitation and also to overcome the restriction of the conventional linear VEHs having a narrow operating bandwidth [30,31]. The nonlinearity of the device is presented in many works in different ways. The researches were oriented to change the natural frequency of a system by controlling its geometrical characteristics [32,33], by introducing external forces [34], via the interaction of the oscillator with the magnetic field [35][36][37][38] or by imposing high displacements [39], etc. These methods proved that the introduction of nonlinear dynamics increases significantly the frequency bandwidth of the devices. In recent works published by the authors of this paper, the benefits of the combination of the nonlinearity and the energy localization phenomenon have been investigated in a two degrees-of-freedom (2-DOFs) system [40,41]. It has been proved that the introduction of nonlinear dynamics enhances the frequency bandwidth and offers a higher energy localization robustness compared to the linear systems. Moreover, it has been confirmed that the functionnalization of the energy localization allows enhancing the harvested power. Furthermore, to produce a high-efficiency output harvester, the optimization of the vibration energy harvesters has gained a lot of interest. For instance, many studies have been conducted to optimize the output power and the structure aspects of the harvesters [42,43]. Several optimization approaches have been used such as the Non-Dominated Sorting Genetic Algorithm II (NSGA-II) [44], the new optimization methodology based on Artificial Intelligence (AI) [45], etc.
In this work, a multimodal electromagnetic vibration energy harvester with tuning nonlinear dynamics and energy localization is studied. The results of [40,41], allow the experimental validation of the 2-DOFs model. The validated 2-DOFs system is optimized using a multiobjective optimization procedure to improve the VEH harvested power and frequency bandwidth. An efficient criterion using the modal kinetic energy of the finite element model is proposed to quantify the energy localized in the perturbed zones of the quasiperiodic structure. This concept is then generalized to a 5-DOFs VEH where two DOFs oscillators among five are perturbed. A complete optimization procedure using the NSGA-II has been conducted to search for the optimal position of the introduced perturbations. Subsequently, a multiobjective optimization with tuning nonlinearity and mode localization is performed to enhance the harvesting output performance. The obtained results prove that the proposed method provides an efficient tool to design a large harvester with promising performances in terms of harvested power and frequency bandwidth compared to the current state of the art. Figure 1 illustrates the generic proposed vibration energy harvester. It consists of N neodymium magnets guided by elastic steel beams and weakly coupled by repulsive magnetic forces. These moving magnets are placed between top and bottom fixed magnets. The coupling is tuned by varying the gap between magnets thanks to the changement of position of the beams inserted into threaded rods. Wire-wound copper coils are wrapped around the moving magnets. When the device is subjected to a harmonic base excitation Y = Y 0 cos(ωt) where Y 0 is the imposed acceleration amplitude, each moving magnet oscillates around its equilibrium position. Consequently, a current is induced in each coil according to Lorentz' law. v n (n = 1, . . . , N) quantifies the displacement of the magnets. The two magnets at the ends of the array are considered to be fixed so that v 0 = v N = 0. The periodicity of the structure is broken by mistuning few magnet masses.

Equation of Motion
The fourth order partial differential equations of the continuum system are derived using the Hamilton principle as detailed in Appendix A. The following equation of motion is then obtained: where v n stands for the transverse displacement of the nth beam (n = 1, .., N). and . denote respectively the derivatives with the spatial variable and the time. L, L c , ρ, I, E , c m , c e , F m , S and M are respectively the beam half-length, the magnet half-diameter, the steel density, the beam quadratic moment, the steel's Young modulus, the mechanical damping, the electrical damping, the magnetic force, the beam section and the magnet mass. For reasons of symmetry of the 1st mode, each beam is fixed at x = 0 and guided at x = L, so the associated boundary conditions are: To transform the continuous multiphysics problem into a system of discrete ordinary differential equations in the time domain, the Galerkin modal decomposition is used and is detailed in Appendix B. Consequently, dividing Equation (A11) of Appendix B by the equivalent mass Meq, the equation of motion of the nth DOF is written in terms of generalized coordinates as follows: ä n + cȧ n + ω 2 0 [(1 + 2β) a n − β (a n−1 + a n+1 )] + f nl a 3 where c = c eq /M eq is the equivalent viscous damping, β = k l mg /k l mec is the coupling factor, f nl = k nl mec /M eq is the mechanical nonlinear term while assuming that the magnetic nonlinear term is neglected compared to the mechanical nonlinear term (k nl mg /k nl mec = 0.28%), p stands for a mass ratio and ω 0 is the eigenfrequency of the decoupled 1-DOF oscillator. Y 1 and t 1 are defined in Equation (A6). M eq , c eq , k l mg , k l mec and k nl mec are defined in Equation (A12).
The mass mistuning coefficient α n is introduced in Equation (3). Thus, the equation of motion of the nth magnet can be written as follows: where α n = 1 for the mistuned DOF n α n = 1 for the non-mistuned DOF n with n = 1, .., N, as we consider two additional fixed magnets at the ends of the array: a 0 = a N+1 = 0.

Modal Localization Phenomenon
The energy localization phenomenon occurs under conditions of internal weak coupling in nearly periodic structures. In fact, when a small disorder is introduced, the symmetry of the periodic system is broken leading to energy confinement in the perturbed region. To illustrate the energy localization phenomenon, a multi-degree-of-freedom system consisting of an array of 10 weakly coupled beams considered to be perfectly periodic, shown in Figure 2, is proposed. A finite element method is, then, developed to study this system. The elastic beams are bi-clamped and the coupling between the beams is insured via spring connexions where their stiffnesses are fixed so that the coupling is weak ( 1%). The assumption of the weak coupling leads to the creation of closed modes where all normal frequencies ω n may be expressed with respect to the reference frequency ω 0 as follows: where i = 1, 2, ..., N, ∆ = 2 β (1 − cos( iπ N+1 )) << 1 and β 1% is the coupling factor. The Taylor expansion of cos( iπ N+1 ) gives the following approximate solution of all normal frequencies: This assumption permits to create closed modes in order to study the energy localization phenomenon. The natural eigenfrequencies of the 10-beams system are reported in Table 1. Table 1. Natural eigenfrequencies of the 10-beams system in Hz.
ω 10 120.8924 120.8996 120.9109 120.9253 120.9418 120.9590 120.9756 120.9900 121.0013 121.0084 The first mode shapes and their sensitivity to the introduced mistuning are shown in Figure 2. The vibrations of the periodic system, observed in Figure 2a, are uniformly distributed. Then, the mistuning was introduced to the 10-beams structure by varying the density of three arbitrary beams among 10 by 6% (the 3rd, 6th and 9th beams and the 2nd, 5th and 8th beams counting from the bottom in Figure 2b,c, respectively). The perturbed beams have notably more displacements than the others and this is modeled in terms of kinetic energy confined in the mistuned regions shown in Figure 2b,c. Comparing the results of the configurations illustrated in Figure 2c,d where the irregularities are placed in the same positions but with different amount, the energy localization is more accentuated and kinetic energies of the perturbed beams are significantly higher then the others in the case of the density mistuning by 10%. Figure 2. Bending vibration mode of the coupled beam system subjected to variations in mass density by (a) 0% (periodic system), (b) 10% of the 2nd, 5th and 8th beams counting from the bottom, (c) 6% and (d) 10% of the 3rd, 6th and 9th beams counting from the bottom.
Consequently, to quantify the mode localization and predict its occurrence, a criterion is proposed. The ratio between the modal kinetic energies of the concerned local DOFs and the global structure will be calculated, for each mode, as follows: where M n is the mass matrix of the area made up of n elements, x n is the eigenvector restriction of this area, X and M are respectively the eigenvector of the considered mode and the mass matrix of the full model. According to that, E 12 applied to the 1st and 2nd DOFs in periodic ( Figure 2a) and quasiperiodic ( Figure 2b) configurations are respectively of 20% and 50%. It can be concluded through this example that the localized energy in 2-DOFs of the quasiperiodic system is equal to the energy of 5-DOFs in the periodic system.
In the particular case of discrete system with N DOFs, this criterion can be written as follows for the DOF n: where m n and x n are respectively the mass and the displacement of the considered DOF n.

Numerical-Experimental Confrontation
The concept combining the benefits of geometric nonlinearities and energy localization is proposed to enhance the performances of a periodic weakly coupled electromagnetic VEH device. These effects on the frequency bandwidth as well as on the harvested power have been investigated in [40]. Following this study, it has been concluded that the energy can be harvested from one perturbed magnet instead of two and only one electric circuit can be used while the maximum harvested energies are comparable and the bandwidth is enhanced compared to the results of the periodic structure. Subsequent work included the manufacturing and experimental characterization of the proposed device under harmonic excitations [41]. The nonlinearity level has been tuned by varying the critical resistances and the energy localization has been controlled by the mass mistuning. Through this study, it has been confirmed that the experimental tuning of these phenomena allows the enhancement of the VEH performance in terms of frequency bandwidth and harvested energy and the compromise solution that maximizes simultaneously the main objectives has been experimentally determined. Moreover, it has been proved that the nonlinear dynamics offer a higher robust energy localization compared to the linear system. Despite the significant improvement in the performances of the 2-DOFs device, a complete multiobjective optimization needs to be performed for further enhancement of the generic harvester. First, a system of 2-DOFs is studied. The power is harvested only from the oscillations of the 2nd perturbed magnet. Therefore, the harvested power is expressed as follows: where A 2max , R load , R int and δ stand respectively for the maximum amplitude of the frequency response of the second perturbed DOF oscillator, the transduction circuit load resistance, the coil internal resistance and the electromagnetic coefficient.
The system of equations of a 2-DOFs model where the second DOF oscillator is perturbed, is generated according to Equation (4). It is solved using ode45 method. The amplitudes and the harvested power are calculated. All the numerical simulations are performed with a basis acceleration equal to a rms = 1 g, a gap d = 50 mm and a coupling coefficient β = 0.11%. The kinetic energies based on the criterion defined in Equation (8) while varying the value of the mass mistuning is plotted in Figure 3a. The largest difference between the kinetic energies is obtained when α = 1.06 as observed in this Figure. Based on the test bench described in [41], several experiments have been performed in order to validate the numerical results. To do that, the magnet masses are perturbed by adding small masses to the magnets. In order to accurately measure voltage peak, up and down frequency sweeps for all the following experimental tests, which allow the capture of the bifurcation points of the nonlinear frequency response, have been done. The amplitudes and masses of the 2-DOFs oscillators are being experimentally obtained while varying the mass mistuning coefficents, the corresponding kinetic energies are calculated through Equation (8) and plotted in Figure 3a. As shown, the optimal mass mistuning is α = 1.06 and is in good agreement with that of the numerical simulations.
The harvested power with load resistances is numerically calculated for the optimal mistuning value α = 1.06 based on Equation (9). While varying the load resistances during the experimental tests, the current flowing in each load resistance R load provides an electric power calculated as P = max(V) 2 /R load where V is the voltage generated by the coil. As shown in Figure 3b, the optimal load resistance R * load which gives the maximum power P max is 6 Ω for both numerical and experimental results.  A quantitative comparison between theoretical and numerical results is achieved. In fact, the optimal mass mistuning and the optimal load resistance are in good agreement with the experimental results (Error < 1%). The model predicts the values of the maximum harvested powers with error of 2%.

Optimization of the Validated Model
During simulations, it was possible to find a maximum power but a minimum bandwidth. While the performance of the harvester in terms of harvested power and frequency bandwidth has been improved by tuning nonlinearity and mode localization, a multiobjective optimization needs to be done to further enhance the harvester and obtain its optimal parameters. This approach has been applied in many fields where optimal decisions need to be taken in the presence of trade-offs between several objectives.
In the following, this procedure is introduced using an extension of the Genetic Algorithm (GA) [46] for multiple objectives which is the Non-dominated Sorting Genetic Algorithm II (NSGA-II) [47]. The latter is related to evolutionary multiobjective algorithm that aims at improving the adaptive fit of a population to a Pareto front composed of a set of compromise solutions between the objectives. The formulation of the problem is defined in Table 2.

Objective Function
Maximize: (P, BW) = f (α, β, R load ) Constraints Where R c stands for the critical load resistance used to tune the nonlinearity level. It has been obtained from experiments performed in [41] and it is equal to 15 Ω. The load resistances used in the following optimization procedures should be higher than R c to guarantee a nonlinear behavior of the device.
The simulations performed generate the frequency responses of all the candidate solutions using the ode45 function of MATLAB. While assuming that the bifurcation point of the nonlinear response curves coincides with the maximum value of the frequency response, the bandwidth and the average harvested power are calculated.
After running the multiobjective algorithm convenient to this problem, the following Pareto front, illustrated in Figure 4, is generated.  The compromise solution that suits for the present problem is the one which maximizes simultaneously the two objective functions. Therefore, the selected Pareto solution and its corresponding parameters are defined in Table 3. Table 3. Results of the multiobjective optimization of the 2-DOFs harvester.

Compromise Solution
(P * , BW * ) = (69.84 µW·cm −3 ·g −2 , 1.038 Hz) With α * = 1.03 The corresponding optimal parameters α * , β * and R * load are reproduced in the experimental tests to compare the numerical results to the experimental ones. The measured performances are illustrated in Figure 5. The error between numerical and experimental results in terms of maximal harvested energy and maximal frequency bandwidth are respectively 5% and 0.5%. Thus, we can consider that the optimized model is in good agreement with the experimental results.

Optimization of a Multiple-Degree-of-Freedom: Five-Coupled-Beams
While a proof-of-concept of the 2-DOFs structure was designed, modeled, fabricated and characterized, demonstrating improved power and bandwidth performance, the model is generalized to a quasi-periodic 5-DOFs and its analytical solution is derived.

Optimal Position of the Introduced Mistuning
The governing equations of the generalized periodic 5-DOFs system is written as follows: Based on the characterization of the previous configuration and the benefits of the combined techniques, two moving magnets' masses are mistuned among the 5-DOFs. Each perturbed magnet mass will have its own mass mistuning coefficient and its own corresponding load resistance. As a choice, we set the same gap between the magnets. The possible combinations of the positions of the two mistuned DOF oscillators are: (a,b), (a,c), (a,d), (a,e), (b,c), (b,d) and (b,e) where (a) for example stands for the DOF (1) as illustrated in Equation (10). The decision of the position of the two DOFs to be perturbed will be made according to the case that provides the maximum energy harvested as formulated in Table 4. Table 4. Objective function and constraints applied to search for the position of the mistuning.

Objective Function
Maximize: P = f (α 1 , α 2 , β, R 1 load , R 2 load ) Constraints To do that, a continuous optimization procedure involving discrete variables is used [48]. The set of discrete variables contains α 1 and α 2 which stand for the mistuning coefficients of the two perturbed oscillators. For these two values, the multiple possible combinations are treated. A continuous optimization algorithm, implemented in MAT-LAB, is called according to each combination of discrete variables. During the continuous optimization, those discrete variables denote the position and the value of the mistuning simultaneously. For that, they are modeled as a first step as discrete and then as continuous. In the case of two successive perturbed magnets, the combination is eliminated and not treated. This choice is made in order to have distributed perturbations over the network so that the energy localization is more efficient. During the conducted simulations, this system of equations is solved by the ode45 MATLAB solver. The optimization returns the position of the mass mistuning and the maximum harvested energy of the best configuration. Results indicate that the 5-DOFs system with mistuning of the 2nd and 4th DOFs, illustrated by the equivalent model in Figure 6, give the maximum harvested energy. Thus, the harvested power will be calculated through this expression in the following study: P = P 2 + P 4 , where P 2 and P 4 are the harvested powers from the DOFs 2 and 4, respectively.

Multiobjective Optimization of the Five-Coupled-Beams Harvester
Aiming to simultaneously enhance the harvested energy and the frequency bandwidth of the 5-DOFs harvester with two mistuning, a multiobjective optimization is carried based on the genetic algorithm. The Pareto front of the kinetic energies of the mistuned magnets is plotted, as illustrated in Figure 7a, in order to visualize their trending variations. From one solution to another, the kinetic energies of the two perturbed dofs are conflicting. To improve the harvester performance, the kinetic energies of the perturbed DOFs should be maximized so their sum will be maximized in the following. In order to take advantage of the mistuned DOFs, both of them should vibrate in close proportions so they should have close energies. For that, during this procedure, the energy rates of the mistuned DOFs are controlled. This constraint is considered as the additional subjective preference information to choose the optimal solution enhancing the harvester. To do that, the energy localization rate between the 2nd and the 4th DOFs is defined as follows: Then, the mutiobjective optimization problem is formulated in Table 5 as follows: Table 5. Objective function and constraints applied to the optimization problem for the 5-DOFs harvester.

Objective Function
Maximize: (P, BW, τ 24  Pareto optimal solutions of the multiobjective optimization are depicted in Figure 7b and the two-dimensional projection is reported in Figure 8a. Among the multiple Pareto solutions, one will be chosen. The selection criterion is based on the solution that minimizes the difference between E 2 and E 4 . To do that, for each solution located inside the compromised interval, the corresponding difference between E 2 and E 4 is calculated and is illustrated in 'red' in Figure 8a. Consequently, taking into consideration this additional subjective preference information, the compromised solution that maximizes the harvested power and the frequency bandwidth simultaneously is reported in Table 6. Table 6. Results of the optimization problem of the 5-DOFs harvester.

Compromised Solution
(P * , BW * ,τ * 24 ) = (140.923 µW·cm −3 ·g −2 , 1.87 Hz, 64.45) Introducing these optimal values in the model, the frequency response in terms of harvested power is calculated and plotted in Figure 8b.  The simultaneous functionalization of the nonlinearity, the energy localization phenomenon and the multimodal configuration show an improvement up to 101% of the harvested power and 79% of the frequency bandwidth compared to the performances given by the 2-DOFs nonlinear system with perturbing only one DOF. The 5-DOFs periodic system is also studied in order to compare its performance to the one of the optimized quasiperiodic 5-DOFs system. In this case, the power is harvested from all the vibrating 5-DOFs oscillators (P = ∑ 5 n=1 P n ). Following the same strategy, the Pareto front of the periodic system maximizing the harvested power is plotted in Figure 9. The optimal parameters that maximize simultaneously the harvested power and the frequency bandwidth for the periodic 5-DOFs system are β = 1.1%, R load = 22 Ω. The obtained results show that the proposed quasiperiodic model provides a larger bandwidth and comparable harvested power compared to the periodic one. In fact, the difference between the harvested powers in the two cases is of 6.15% and the frequency bandwidth is higher by 18% in the periodic system. From these results, it can be concluded that the functionnalization of nonlinearity and mode localization in a multimodal device overcome the challenge of increasing the harvested energy and the frequency bandwidth. Moreover, the combination of these phenomena allows harvesting energy from only two DOFs oscillators instead of five while keeping comparable performance. The benefits of this property consist of reducing the cost of the electrical circuits to be implemented and the technological constraints of the structure.
To compare the optimized proposed harvester performance with the current stateof-the-art, the volume figure of merit (FoM v ) proposed by [49] has been chosen among diverse performance metrics in the literature for being the most general criterion. The volume figure of merit is the ratio of the harvester useful power output transferred to the load to the maximum theoretical power flowing into an equivalent device. This equivalent device has a cubic geometry with the same volume as the original harvester one but with a proof mass having the density of gold ρ Au occupying this volume half, while the other half is destined to oscillation [49]. Hence, this figure of merit is calculated according to the following expression: where Y 0 is the imposed acceleration amplitude, ω 0 is the resonant frequency, ρ Au is the density of gold and Vol is the harvester volume. According to this figure of merit, the performance of the current work harvester and the ones of other harvesters are illustrated in Figure 10. As shown, the optimized proposed harvester provides competitive performance compared to the harvesters based on electromagnetic transduction as well as the ones based on hybrid piezoelectric-electromagnetic transduction.  [50] [51] [52] [53] [54] [55] [56] [57] [Current work] [35] [23]

Conclusions
In this work, a generic model of a nonlinear quasiperiodic VEH based on electromagnetic transduction has been studied. A 2-DOFs VEH with tunable nonlinearity and mode localization is experimentally validated. Then, the validated 2-DOFs structure is optimized using a multiobjective optimization procedure with respect of the harvested power and the frequency bandwidth. To quantify the energy localized in the perturbed zones of the quasiperiodic structure, an efficient criterion based on the modal kinetic energy of the finite element model is presented. Afterward, the concept is generalized to a 5-DOFs structure and the optimal performances are derived using a multiobjective optimization procedure. It has been proved that the optimal parameters that improve the quasiperiodic 5-DOFs device performance enable an enhancement up to 101% and 79% in terms of harvested power and frequency bandwidth, respectively, compared to the 2-DOFs harvester. Moreover, through the performance comparison between the quasiperiodic 5-DOFs system with well-chosen mistuning positions and the periodic one, it has been proved that the functionalization of the nonlinearity and the energy localization allows more efficient frequency bandwidth and comparable harvested powers. Hence, the proposed quasiperiodic harvester is improved by reducing the transduction circuits number and the manufacturing cost. In addition, it has been shown that the optimized harvester has competitive performance compared to the current state-of-the-art. Finally, despite the fact that the generic harvester performances are optimized, an efficient energy harvester for small-scale realistic applications is challenging. Future work will include the miniaturization and experimental validation of the device with more than 2-DOFs as well as the optimization of its performance while investigating the nonlinear energy localization phenomenon and its stability [58].

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Dynamic Equations of Motion of the Continuous Structure
Based on the problem formulation detailed in Mahmoudi et al. [35], the continuum multiphysics system, including the equation of motion of the present structure and the magnetic transduction equation, is developed. To simplify the development of equations, the magnet mass has been considered to be a localized point mass for the mode of interest (1st). FEM simulations under ANSYS have been run to support this hypothesis. The illustrations of the distributed magnet mass and of the localized point mass are shown in Figure A1a,c, respectively. As reported in Figure A1b,d, the frequency of the beam and its mode shape are the same for the two configurations (Errors of 0.005%, 0.003% and 0.09%, respectively in terms of displacement in the center, slope and frequency).

Appendix A.2. Equation of Motion
The fourth order partial differential equations of the continuum system are derived using the Hamilton principle. Applying the Hamilton's variational approach to the dynamical system, we have the following equation: where W n ci , W n p , W n ext are respectively the works of the non-conservative forces, damping forces and exterior forces of the nth beam (n = 1, .., N).
The magnets and the beams are identical and the distance d between the magnets is the same. When the magnets oscillate, magnetic and electromagnetic forces are created as a result of the variation of the magnetic field. The magnetic force depends on the magnetic intensities Q and the gap d is defined as follows: where Q = H c S m , H c is the coercive force, S m is the section of the magnetic pole and v n denotes the transverse displacements of the nth beam (n = 1, ..., N).
For δH = 0, the following equation is obtained: Figure A1. L, ρ, I, E , c m , c e , F m , S and M are respectively the beam half-length, the steel density, the beam quadratic moment, the steel's Young modulus, the mechanical damping, the electrical damping, the magnetic force, the beam section and the magnet mass.
For reasons of symmetry, the beam is fixed at x = 0 and guided at x = L, so the associated boundary conditions are: v n (0, t) = v n (0, t) = 0 EIv n (L, t) = Mv n (L, t) and . denote respectively the derivatives with respect to the spatial and the time variable.
The magnetic transduction is provided by a coil wound around the separation distance between two consecutive magnets. The following mechanical-magnetic coupling equation is defined as follows: i n (t) = δ R load + r intv n | x=L (A5) where δ is the electromagnetic coefficient and R load and r int are respectively the load and the internal resistances.
To simplify the equations, the following dimensionless variables are introduced: The coupled continuum multiphysics problem is then equivalent to the following equation:

Appendix B. Reduced Model by Galerkin Modal Decomposition
A reduced model is generated by Galerkin modal decomposition, transforming the continuous multiphysics problem into a system of discrete ordinary differential equations in the time domain. In order to solve the previous system, the displacements are projected on a single mode basis. The displacement of each DOF oscillator n is written as follows: u n (X, t 1 ) = φ(X)a n (t 1 ) where φ(x) is the projection base and a n (t) are the generalized coordinates.
To simplify the modal projection of the electromagnetic force, the latter is developed in Taylor series at order 3 and as already proved by Mahmoudi et al. [35], the magnetic force can be expressed as a combination of the linear stiffness k l mg and the cubic stiffness k nl mg as follows: F m = k l mg u n + k nl mg u 3 n (A9) It is assumed that the magnetic nonlinearity is neglected compared to the mechanical nonlinearity in the case of a weak coupling between the two beams.