Robustness of nonlinear electromagnetic vibration energy harvester subjected to random excitation

Vibration Energy Harvesters (VEHs) are devices used to collect mechanical energy from the surrounding environment to supply low power electronic systems such as Wireless Sensor Nodes. In this paper, we introduce an electromagnetic VEH model and a semi-analytical method called Moment Equation Copula Closure (MECC) that is compared to Monte Carlo simulations. Those methods are then used to derive the maximum power that can be extracted from random vibration before analyzing the effect of cubic stiffness nonlinearity on the VEH robustness against the variation of the excitation spectrum. Unlike bistable nonlinearity, it is shown that Dufﬁng nonlinearity can be used to enhance the VEH power density and robustness with a limited effect on the harvested power.


INTRODUCTION
The Internet of Things (IoT) is the term used to describe a cluster of technologies enabling machine to machine (M2M) communication and machine to human interactions through the Internet [1]. Among the technologies of the IoT, Wireless Sensor Network (WSN) senses and communicates to gather data from different locations. Apart from a gateway and a base station, WSN are usually composed of a bunch of wireless sensor nodes cooperating to transmit the sensed data to a common point. Eventually, those nodes can be scattered in an unfriendly environment, which may lead to communication, electronic or mechanical issues. One of these issues concerns the power needed to supply the node; due to their limited lifespan or their inability to withstand extreme temperatures, usual batteries are sometimes not able to satisfy the nodes specifications [2]. Therefore, alternatives to batteries and power cables are needed to supply the wireless sensor nodes. Recent advances in the energy efficiency of electronic devices and the work done to optimize power management [3] place energy harvesting as a viable solution. Energy harvesting aims to collect energy from the environment, which is usually considered as losses or constraints, and convert it to usable electrical power for electronic devices. The surrounding source of energy can be radiation (light, radio wave), thermal gradient or mechanical movement [4,5].
We focus on harvesting the energy from mechanical vibrations (i.e. Vibration Energy Harvesting). One of the first Vibration Energy Harvesting study dates back to the 90s [6] and introduces the widely used basic principle of vibration energy scavenging: the device is modelled as a single degree of freedom mechanical oscillator which concentrates the vibration energy in the form of kinetic energy that is then converted into electrical energy through transduction mechanism. Generally, the electromechanical transduction is provided by one or a combination of the following three phenomena [7] : piezoelectricity, electromagnetic induction (based on the Lenz-Faraday law of induction) and the electrostatic effect. We will focus on electromagnetic vibration energy harvester (VEH) because, unlike the others, this transduction technic is already widely known and used in the market and has virtually an infinite life time (no wear nor mechanical stress) while providing high electromechanical coupling coefficient.
Linear VEH has been widely studied and is well understood [8]. Its main drawback is its narrow frequency band, which results in low flexibility against changes in the excitations most energetic frequency. Therefore, efforts have been made to overcome this issue (i.e. to widen the frequency band); some of them proposed to use multi degree of freedom (d.o.f.) VEH [9][10][11] but those solutions increase the volume (and/or mass) of the VEH impairing power generated per unit of volume (mass). Another proposed technique is called self-tuning and consist in adjusting the resonant frequency of the harvester to the source either by closed-loop electronic circuit or by using specific dynamic properties; however those techniques prove to be difficult to implement experimentally. In 2009, Cottone et al. [12] and Mann and Sims [13] proposed to introduce on purpose nonlinearity to the stiffness of a linear VEH, respectively a bi-stable and a Duffing type nonlinearity. These promising ideas, mainly because of the effect of the nonlinearity, which tends to bend the frequency response of VEH widening its frequency band, were soon followed by other works on nonlinear vibration energy harvesting [14][15][16][17][18][19][20][21][22].
Most of the previously cited studies deal with harmonic vibrations sources. However, it is of common knowledge that real-life vibration sources are wideband random signals. From this last observation, many researches have been carried regarding vibration energy harvesting under random excitation [23][24][25][26]. In 2014, Daqaq et al. [27] proposes a review of nonlinear energy harvesting. They state that both feasible Duffing type nonlinearities (hardening and bi-stable) can, at best, provide the same performances as their linear counterpart in terms of harvested power under random excitation. Moreover, Daqaq et al. infers that Duffing type nonlinearities could be used to reduce the effects of small drifts in the center frequency of colored noise excitation hence increasing the robustness of the harvester. Though the interest of Duffing type nonlinearities in vibration energy harvesting seems mitigated, the above conclusions remains general and as concedes Daqaq et al. [27] by advising a careful characterization of the excitation source, a case by case analysis would be necessary.
In this paper, we propose to study the case of nonlinear vibration energy harvester that would be attached to a truck. In fact, a truck has a lifespan of around twenty years and no access to a direct power source. Therefore, a WSN coupled to the trailer of the truck (e.g. to monitor its position) would greatly benefit from a vibration energy harvester. First, the vibration random source is presented. Thereafter, the nonlinear vibration energy harvester model with base excitation is presented together with a semianalytical method [28] to approximate the statistical moments of interest. This method is then compared to direct Monte Carlo simulation while analyzing the behavior of the nonlinear VEH under vibrations iduced by the truck. Finally, the robustness of the nonlinear VEH against changes in the excitation frequency spectrum is studied.

VIBRATION SOURCE
Using seven measurements and analysis studies of trailer truck vibration made all around the world [29][30][31][32][33][34][35], we have been able to characterize the vibration source which we are interested in.
The considered vibration source is a stationary random excitation signalz(t) characterized by its Power Spectral Density (PSD), which the Fourier transform of the autocorrelation function. As can be seen on Fig. 1, all of the gathered PSD contains three main peaks corresponding to different vibration modes of the trailer: suspension, tires and structural modes. The most energetic mode is the one linked to suspension of the trailer and can be approached by a fitting function Szz given by equation 1 where ω is the angular frequency and A, ω t , ξ t are function parameters. This function is build taking into account the road roughness and the response of the suspension mode to it, a fitting example is also given on Fig. 1.
By fitting all the gathered PSD, we observe that the most energetic frequency of the vibration source f 0z varies from 1.7 to 4.2 Hz with a mean of 2.5 Hz. Harvesting this source hence is a real change for linear vibration energy harvesters dur to their narrow frequency band.

NONLINEAR GENERIC MODEL
In this paragraph, we describe the model of an electromagnetic nonlinear VEH depicted in Fig. 2.
The base of the VEH is submitted to a displacement z(t) (the vibration source) which causes a relative displacement x(t) of the mass m suspended to the base by a spring of stiffness k and a nonlinear cubic stiffness k 3 that can be implemented using magnets or a beam in large displacement [13,36]. We consider the mobile mass as a source of magnetic flux (magnet) and a copper coil of internal resistance R int and inductance L is wrapped around it creating an electromechanical coupling with coefficient α. The energy of the moving mass is then dissipated as mechanical losses (damping c m ) and electrical current i(t) in the copper coil. This electrical energy is then lost as heat in the internal resistance of the coil and harvested in the load modelled here by a resistance R l .
The mechanical movement x of the mass m can be modelled by equation (2).
The electric current i flowing through the coil by equation (3).
A usual and trustworthy assumption is to consider the role of the coil inductance L as negligible. Consequently, equation (2) and (3) decouple and only one differential equation maintains. Also, the instantaneous power P(t) harvested on the load is given by equation (4).
Thus, to obtain the harvested power, one has to solve equation (5) and use equation (6) knowing the speedẋ of the moving mass.ẍ Because we are working with stationary random excitations, x and P will have a random nature and exhibit stationarity after an initial transient response. To have meaningful performances metrics, we consider only the stationary regime and apply statistical averaging to the quantities of interest. Hence, we focus on the expected harvested power P(t) that is to be obtained by computing ẋ(t) 2 , the second order spectral momentum of x(t) since x(t) is zero.
The Duffing-type nonlinearity added to our electromagnetic VEH leads its mechanical oscillator to several types of behavior depending of the value of k 1 and k 3 (see Fig.3a).
If k 1 > 0 and k 3 = 0, the VEH has a linear behavior with a resonance angular frequency ω 0 equals to (k 1 /m).
If k 1 ≥ 0 and k 3 > 0, the VEH behaves as a Duffing oscillator characterized by a stiffness hardening. At the equilibrium position (x = 0), the stiffness is equal to k 1 ; then we consider the linear resonant frequency ω 0 as in the linear case.
If k 1 < 0 and k 3 > 0, the VEH has two stable equilibrium positions in ± (−k 1 /k 3 ) and can oscillate around or between those two positions: it presents a bi-stable behavior. At the stable equilibrium positions, the stiffness is equal to 2k 1 and we consider ω 0 = −2k 1 /m. Any other arrangement is not feasible in the vibration energy harvesting context. In the lexicon of oscillators, we refer to potential well. Those potential wells correspond to stable equilibrium positions that can be graphically observed as minima of the potential elastic energy, plotted on Fig. 3b. For bi-stable oscillators, we therefore speak of intra or inter-well motion when the mass oscillate around and respectively between its two stable equilibrium positions.

MECC METHOD
To solve equation (5) whenz(t) is a colored noise, several methods are available [37]. The most accurate uses Monte Carlo analysis of direct numerical computation of the oscillator response, however, this method is very time consuming since a large amount of simulation must be done on a sufficient time period. For efficiency reasons, we decide to propose a less accurate but faster method for a first analysis: the Moment Equation Copula Closure (MECC) method introduced by Joo et al. [28]. This method is pratically described below.
To present this method, we consider the stochastic differential equation (7) that can be easily linked to equation (5).
Assuming that bothz and x have zero mean and are stationary, one can obtain the moment equations given by (8) and (9).
Where s and t are two time values and τ = t − s. Also, the function C is the covariance function; e.g. C xz (τ) is the covariance function of x and z at time τ.
The nonlinear term (k 3 ) implies the presence of a forth order statistical moment term in equations (8) and (9). This excludes a direct resolution of these moment equations. So, a closure technique is necessary to approach the solution of the system. The Gaussian closure technique [38] is the reference when the oscillator is unimodal (i.e k 1 > 0) and supposes that x has a Gaussian probability density function (pdf). This last assumption doesnt hold for bi-stable oscillator. Thus, we consider a pdf f based on the analytical solution of the Fokker-Planck equation under white noise excitation, it is defined in equation (10). With: Where γ is a free parameter related to the energy level of the system.
As indicated by its name, the main step of Joo et al.s MECC method is to provide a closure of the moment equations (8) and (9) by using Gaussian copula densities to find the linear relations given in equations (12) and (13). Those relations are approximation obtained from Taylor expansions.
In our case, it turns out that ρ xz and ρ xx have the same expression which is given by equation (14).
Where er f (·) is the Gauss error function and F(x) the cumulative distribution function obtained from the pdf given by equation (10). (8) and (9) and using the relation given in (12) and (13), we are able to use the Wiener-Khinchin theorem to obtain the PSD equations (15) and (16).

From the moment equations
Then, we obtain the expression of the PSD of the displacement response S xx as a function of Szz.
From equation (17), we can derive the expression of the variance of the response x(t) 2 : The main step of the MECC method is to find the triplet (γ,ρ xz ,ρ xx ) which allows to respect at best equations (14) and (18). A cost function J of (γ,ρ xz ,ρ xx ) equals to the sum of the squared difference of the left and right terms of prcited equations.
Finally, the minimization of J give us a triplet (γ,ρ xz ,ρ xx ) that is the key to compute approximate values of most of the statistical quantities of the oscillators response. Among those quantities, ẋ(t) 2 is the most interesting for us and is obtained from equation (19) that comes from the linear property of time derivation in the frequency domain.
To obtain reference analytical values, we analyze the linear VEH response under truck trailer vibration as characterized before. Linearity implies that x(t) will be Gaussian as long asz(t) is Gaussian and equation (19) turns to equation (20).
For sake of simplicity, we consider that the harvester is tuned to the most energetic frequency of the excitation f 0z (i.e k 1 = (2π f 0z ) 2 ). Therefore, a white noise approximation as defined by Preumont [39] is acceptable and we obtain equation (21). S max is the maximum value of Szz.
Taking back the notations of equations (5) and (6), we are able to write the expression of the expected power (22).
Impedance matching (i.e find the load impedance R l that maximizes the harvested power) is done for load resistance equals to R int (1 + α 2 /(R int c m )) . Moreover, we consider an optimized VEH that is the electromechanical coupling coefficient is maximized and the losses (electrical and mechanical) are minimized (α 2 /(R int c m ) → +∞). In this case, the expected harvested power reaches its maximum P(t) max value given in equation (23).
The result expressed in equation (23) corroborated state of the art findings [26]. Also, it is worth noting that it is highly valuable, as an industrial point of view because it allows one to instantaneously judge the feasibility of a vibration energy harvesting solution by comparing the power needs to the extractable power from a known vibration source.
From this development, we can also retain some design rules for a linear VEH under random vibration, that is to say the parameters to minimize (R int ,c m ) and to maximize (α); furthermore, it is interesting to note that the expression α 2 /(R int c m ) can be used to quantify and compare the performances of linear electromagnetic VEHs.
As our study focuses on the VEH nonlinearities, we consider in the following that the VEHs intrinsic parameters (m, c m , α, R int ) are fixed to realistic values extracted from one of our prototype and are listed in Table 1.

NONLINEAR VEH UNDER TRUCK VIBRATION
In this section, we focus on the response of the nonlinear VEH submitted to the vibration source characterized during truck transportation. First the variability of the most enegetic frequency is not taken into account and we have: S max equals to This approach has two purposes, the first one is to analyze the behavior of the nonlinear VEH under truck vibrations. The second one is to compare both computation methods and identify the weakness of the MECC method in terms of accuracy for future uses.
Note that for any results presented, impedance matching is satisfied, that is to say: R l is set to the value that maximizes P(t) .

Hardening nonlinearity
On Fig. 4, we plotted the expected harvested power P(t) for different value of ω 0 /ω max (ω max is the angular frequency corresponding to f 0z ) and k 3 under excitation with the nominal PSD. As it is the most accurate, we first consider only the results obtained through the Monte Carlo simulations (Fig. 4a): 100 samples. For values of the nonlinear coefficient k 3 under a given threshold, the VEH behave linearly: the harvested power is maximized (42 mW) when the VEHs linear resonance frequency is tuned to the most energetic frequency of the excitation ω max . When nonlinearity is sufficient, the VEHs apparent stiffness increases and lower values of linear stiffness are necessary to guarantee tuning to the excitation; a slight decrease of the harvested energy is then observed. For even higher values of k 3 , low linear stiffness cannot maintain tuning and the harvested power drops significantly. In the present case, we note that the hardening-type nonlinearity has no positive effect on the harvested power. Fig. 4 also helps us to appraise the accuracy of the MECC method for hardening-type nonlinear VEH submitted to colored noise. Results obtained using the MECC method (Fig. 4b) shows good agreements with those obtained with Monte Carlo simulations (Fig. 4a). Only the decrease of the power, when the detuning due nonlinearity is balanced by low linear stiffness, is badly captured. This weakness is mainly caused by lack of accordance of the chosen pdf of equation (10) with true displacement pdf. Even if the MECC method shows good performances, it cannot be used to draw conclusions about expected harvested power of a hardening-type nonlinear VEH under truck vibration. However, one as to mention that MECC method is much faster than Monte Carlo simulation which needs approximately 60 hours to compute Fig. 4 against 1.3 hours for its counterpart. Therefore, to fasten Monte Carlo simulation, MECC method could, for example, safely be used to quickly identify the value of R l that satisfies impedance matching.
To illustrate this ability of the MECC method, we plotted P(t) against R l using Monte Carlo and MECC methods on Fig.  5. The chosen point (ω 0 = ω max /2 and k 3 =6.25·10 4 N.m −3 ) being one where MECC exhibits its poorest performances. Observing Fig. 5, we note that despite an obvious difference between the results, those obtained through MECC method follows the same trend as those found using Monte Carlo. Hence impedance matching can be identify using the MECC method.

Bistable nonlinearity
On Fig. 6a, we plotted the expected harvested power P(t) for different value of ω 0 /ω max (note that ω 0 is here equal to −2k 1 /m) and stable equilibrium positions ( −k 1 /k 3 ) under excitation with the nominal PSD. Focusing on results obtained from Monte Carlo simulations, we first note that the maximum harvested power (37 mW) is lower than that of unimodal VEH in the given conditions. Then, it can be seen that when the equilibrium positions are close to each other( −k 1 /k 3 =0.1 cm), the harvested power is very low; this can be interpreted saying that the apparent stiffness of the harvester (which exhibits fully inter-well movement) cant tune the excitation frequency spectrum (the VEH would have to be unimodal). When the equilibrium positions deviate from one another, intra-well movement tend to be more frequent until almost no jumps occurs between the equilibrium positions. At this point ( −k 1 /k 3 =1.5 cm), the harvester can be seen as a unimodal oscillator (the mass moves around only one position) and we see that it gives the most power when ω 0 = ω max . Once again, in the above situation, the nonlinearity brings no benefit to energy harvesting.
As done with hardening type nonlinearity and Fig. 4, we can have hints about the performances of the MECC method when comparing Fig. 6a (Monte Carlo) and Fig. 6b (MECC). For small stable equilibrium positions ( −k 1 /k 3 <1 cm), the two methods give very similar trends assessing the good performances of the MECC methods. However, this last method completely fails to capture the behavior of the bi-stable VEH when the equilibrium positions are too far from another; in these conditions, the excitation is not sufficient to enable intra-well are more frequent than inter-well movement and minimization of the MECC cost function doesnt adequately satisfyneither equations (14) nor (18).  This weak point of the MECC method has been pointed out by Joo et al. [28]. Once again, one note that the MECC method is much faster than Monte Carlo simulations.

Variation of the excitation frequency
In this section, we aim to test, in our application, a property of the hardening nonlinearity which tends to robustify the performances of the harvester against variation of the excitation frequency. Our identified vibration source now has its most energetic frequency f 0z (ω max ) varying from 1.7 to 4.2 Hz.
To cut down simulation time, we select three configurations which characteristics are listed in Table 2. Configuration 1 corresponds to the optimal linear vibration energy harvester; Configuration 2 is the nonlinear configuration tuned to 1.7 Hz which gives the best results and Configuration 3 is Configuration 2 with a sligthly stronger nonlinearity.
When submitting those configurations to our nominal configuration, we first note that the nonlinearity has an impact on the maximum displacement of the mobile mass. In fact, the maximum displacement with Configuration 1 is 72.4 mm for a 42 mW power. With Configuration 2, the power drops to 39.5 mW (-6 %) and the maximal displacement to 50.6 mm (-30 %). Finally, with Configuration 3, the harvested power is 38.3 mW (-9 %) while the maximal displacement is 32.6 mm (-55 %). Therefore, we note that nonlinearity decreases maximal displacement with a minor effect on harvested power: it increases the power density of the vibration energy harvester.
Considering now the variation of the excitation frequency, Fig. 7 shows the harvested power against the most energetic frequency of the excitation f 0z . This plot shows that the cubic nonlinearity add robustness to the harvester against variation of the excitation spectrum. For example, we show that Config-  uration 1 produces more that 30 mW between 2.1 and 3.2 Hz while this frequency band is from 2.1 and 3.9 Hz (+64 %) for Configuration 3.
We can then conclude that in our application, the introduction of a Duffing hardening nonlinearity offers two benefits : increase of the power density (power per unit of volume) and gain of robustness against variation of the excitation spectrum.

CONCLUSION
In this paper, we have drawn several conclusions on nonlinear vibration energy harvesting under stochastic excitation. First, using the proposed model, we have shown that the maximal power that can be extracted from a random vibration source is equal to πS max /m where S max is the maximal value of the source acceleration power spectral density and m the mass of the harvester. Then, the MECC method [28] used to solve nonlinear stochastic ordinary differential equations has been presented. It has been shown that for moderate nonlinearity, MECC has a reasonable computational cost with a good accuracy compared to Monte Carlo simulations and can therefore be used as an alternative tool to reduce the computational time.
Finally, we showed that the Duffing hardening nonlinearity offers several benefits in our application which concerns the scavenging the vibrations induced by a truck. In fact, designing a VEH with this kind of nonlinearity increases the robustness of the device against variation of the excitation spectrum and enhances the power density by lowering the maximal displacement with a limited impact on the harvested power.