Stochastic Thermodynamics of an Electromagnetic Energy Harvester

Abstract

We study the power extracted by an electromagnetic energy harvester driven by broadband vibrations. We describe the system with a linear model, featuring an underdamped stochastic differential equation for an effective mass in a harmonic potential, coupled electromechanically with the current in the circuit. We compare the characteristic curve (power vs. load resistance) obtained in experiments for several values of the vibration amplitude with the analytical results computed from the model. Then, we focus on a more refined analysis, taking into account the temporal correlations of the current signal and the fluctuations of the extracted power over finite times. We find a very good agreement between the analytical predictions and the experimental data, showing that the linear model with effective parameters can describe the real system, even at the fine level of fluctuations. Our results could be useful in the framework of stochastic thermodynamics applied to energy harvesting systems.

1. Introduction

Harvesting energy from environmental noise is a fundamental task in many contexts, from the microscales, as in Brownian motors driven by thermal fluctuations [1,2,3,4], to the macroscopic world, as for random vibrations in ships, bridges, railroad cars, and aircraft [5]. In order to treat systems where fluctuations play a central role, in recent years, thermodynamics concepts such as heat, work and entropy have been extended to the realm of stochastic processes and non-equilibrium systems, paving the way to the development of stochastic thermodynamics [6,7]. The interest in the study of these fluctuating quantities, such as heat and work defined along a single stochastic trajectory of the system, is motivated by the search for optimization protocols in models of engines where fluctuations are dominant [8] and by the general symmetry properties that the distribution functions have been shown to satisfy, both theoretically [9,10,11] and experimentally [12,13].

From a more practical perspective, in Internet of Things applications, energy harvesting systems, converting ambient energy into electricity, can represent an alternative solution to disposable batteries [14,15]. That is why energy harvesters have found application in very different fields [16,17,18]; see also [19,20]. In particular, vibration energy harvesters, which convert the wide available mechanical energy of vibrations into electric energy, have proven particularly attractive [21,22]. These devices are mainly based on piezoelectric materials or electromagnetic induction [23,24]. Both piezoelectric and electromagnetic vibration harvesters are typically operated in resonant structures and can efficiently operate only near resonance, even if electromagnetic harvesters are characterized by higher powers with respect to piezoelectric ones. In the literature, the study of vibration harvesters forced by non-sinusoidal or random vibrations has been deeply carried out in case of piezoelectric technology [25,26,27], and in case of hybrid electromagnetic-piezoelectric systems [28,29]. On the other hand, electromagnetic vibration harvesters are typically studied under purely sinusoidal vibrations tuned to their resonance frequency. Analysis of electromagnetic harvesters in the presence of random vibrations can be found in [30,31], but the literature lacks a stochastic analysis of the fluctuations and distributions of their relevant quantities.

In this work, we consider a linear model for an electromagnetic energy harvester driven by white noise, amenable for analytical solution, and compare in detail the theoretical predictions from the model with real experiments. The model consists of an underdamped Langevin equation describing the stochastic motion of an effective mass in the presence of a harmonic potential and electromechanically coupled with an external variable, the current in the circuit, whose dynamics is not directly affected by the noise. The model is, by construction, out of equilibrium, since detailed balance is broken by the coupling with the current variable. This can also be interpreted as a feedback mechanism acting on the mass and introducing a memory effect in the system [32,33]; see also [27] where we performed a similar analysis on a piezoelectric energy harvester. We first fit the effective model parameters to the characteristic curve of the system, namely the extracted power versus the load resistance, and then analyze the time correlation functions of the current and the fluctuations of the extracted work over finite time intervals. We find that the linear model with effective parameters can reproduce very well both the time decay of the correlation function and the work fluctuations at different times.

2. Experimental Setup

A typical electromagnetic vibration energy harvester, schematically shown in Figure 1, is based on a permanent magnet, which is attached to the housing by means of a spring system and can move relative to a coil fixed to the housing when a vibration is applied. The coil is connected to an electrical load, and the relative displacement between the magnet and the coil leads to the conversion of the mechanical energy of vibrations into electrical energy. The maximization of the energetic performance can be carried out on the basis of mechanical tuning techniques, acting on the mechanical characteristics, or on the basis of maximum power point tracking techniques, acting on power electronic converters connected to the coil terminals [34,35,36,37].

The experimental setup that has been employed in all the tests is shown in Figure 2. In particular, the considered harvester is the Model-D by Revibe. It is an electromagnetic vibration energy harvester characterized by a mechanical resonance frequency of about 100 Hz. The harvester is placed on a shaker, the VT-500 by Sentek, used as the source of desired vibrations. In order to generate the shaker vibrations, its driving current is provided by an LA-800 power amplifier. The Spider-81 controller by Crystal Instruments implements a closed loop vibration control by providing the driving signal to the power amplifier on the basis of the acceleration of the shaker vibrations measured by the Dytran 3055D2 accelerometer. It is worth noting that, as shown in Figure 2 (right), the considered electromagnetic harvester comes as a black box that, apart from the forcing by the input vibration, is accessible only at the electric terminals. Therefore, even if its detailed mechanical structure cannot be accessed, it will be shown that a linear model with effective parameters can well reproduce the measured dynamics.

The harvester Revibe Model-D was forced by broadband vibrations of Gaussian type with a sampling rate $f=5$ kHz and different standard deviations ($0.25$ g, $0.5$ g and $0.75$ g), which were generated with MATLAB software and provided to the shaker controller. In Figure 3, the waveforms of the input acceleration, of the voltage and of the power across the load resistance recorded during one of the experimental tests are shown as an example.

3. Theoretical Model

In order to reproduce the dynamics observed in the experimental systems, we consider the following underdamped linear Langevin equation

$$\begin{array}{ccc}\hfill \dot{x}& =& v\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill M\dot{v}& =& -k_{s}x-\gamma v-\theta I+M\xi \hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill L_c\dot{I}& =& \theta v-(R_c+R)I,\hfill \end{array}$$

(3)

where $\xi$ is white noise with zero mean and variance $\langle \xi \left(t\right)\xi \left(t^{\prime }\right)\rangle =2D_0\delta (t-t^{\prime })$. In the above equations, x represents the position of the magnet with respect to the coil, v represents its velocity, M represents the mass, $\gamma$ represents the viscous damping due to the air friction, $k_s$ is the elastic constant of the spring effective system, I is the current at the electrical terminals, $L_c$ and $R_c$ are the coil inductance and resistance, respectively, R is the load resistance and $\theta$ is the effective electromechanical coupling factor. In particular, during its movement, the permanent magnet, apart from the viscous damping force of the medium and the elastic force of the spring, is subject to the electromagnetic (Lorentz) force due to the interaction between the current I flowing into the coil and the induction field of the magnet. As shown in detail in [38], $\theta$ takes into account the coil geometrical properties, its number of turns and the magnetic field strength. It is the coefficient that, in the linear model expressed by Equation (2), links the electromagnetic force and the current at the electrical terminals. This current is due to the electromotive force induced in the coil that, as shown in Equation (3), is proportional to the magnet speed through $\theta$.

In our experimental setup, noise is introduced through the vibrations applied by the shaker to the energy harvester. The role of the shaker mimics the effect of real broad-band vibrations that can be measured for instance in a car or in an aircraft. In our model, such fluctuations are treated as white noise, affecting the dynamical equation for the mass velocity.

The linear model can be solved (see below and the Appendix A for details) to obtain the mean extracted power as a function of the parameters

$$P_{harv}=\langle I^2\rangle R=\frac{D_0{\theta }^2{M}^{2}R}{[\gamma (R+R_c)+{\theta }^2][\gamma L_c+M(R+R_c)]+\gamma k_s{L}_c^2},$$

(4)

where $\langle \dots \rangle$ represents an average over the stationary state.

Fitting Parameters

Using Equation (4), we can fit the parameters of the linear model to the experimental data. Because of the large number of parameters, we also exploit the temporal decay of the correlation function of the current signal (see next section) in order to better estimate some of them (in particular the effective viscosity $\gamma$). We also impose a constraint on the parameters $k_s$ and M, such that $k_s={(2\pi \ast 100)}^{2}M$, dictated by the mechanical resonance of the harvester. We obtain the following values: $\gamma =1.80\pm 0.05$ Kg/s, $M=0.048\pm 0.005$ Kg, $k_s=18811\pm 50$ Kg/s${}^2$, $\theta =29.9\pm 0.5$ N/A, $R_c=227.6\pm 0.5$ $\mathsf{\Omega }$, $L_c=0.124\pm 0.005$ H, which give the analytical curve reported in Figure 4.

4. Temporal Autocorrelation Function

Beyond the study of the mean extracted power, we are interested in exploring the possibility that the model could also well reproduce other features observed in the experimental system. Therefore, in this section, we focus on the autocorrelation function of the measured current $I\left(t\right)$, namely the quantity $C_{II}\left(t\right)=\langle I\left(t\right)I\left(0\right)\rangle$ in the stationary state. This provides information on the characteristic time-scales in the system.

An analytical expression for $C_{II}\left(t\right)$ can be obtained explicitly due to the linear nature of the model. In particular, we introduce the column vector $X={(x,v,I)}^T$ and define the coupling matrix

$$A=\left(\begin{array}{ccc}0& -1& 0\\ \frac{k_s}{M}& \frac{c}{M}& \frac{\theta }{M}\\ 0& -\frac{\theta }{L_c}& \frac{R_c+R}{L_c}\end{array}\right),$$

(5)

so that the Equations (1)–(3) can be rewritten in vectorial form as

$$\dot{X}=-AX+\eta ,$$

(6)

where $\eta ={(0,\xi ,0)}^T$. Then, defining the covariance matrix $\sigma =\langle X^{T}X\rangle$ as

$$\sigma =\left(\begin{array}{ccc}{\sigma }_{xx}& {\sigma }_{xv}& {\sigma }_{xI}\\ {\sigma }_{vx}& {\sigma }_{vv}& {\sigma }_{vI}\\ {\sigma }_{Ix}& {\sigma }_{Iv}& {\sigma }_{II}\end{array}\right),$$

(7)

at stationary, one has the constraint [39]

$$D=\frac{A\sigma +\sigma A^T}{2},$$

(8)

where D is the noise matrix

$$D=\left(\begin{array}{ccc}0& 0& 0\\ 0& D_0& 0\\ 0& 0& 0\end{array}\right).$$

(9)

The solution of Equation (8) provides the covariance matrix as a function of the model parameters (see Appendix A). In order to compute the correlation function $C_{II}\left(t\right)$, we first introduce the response matrix

$$G\left(t\right)=e^{-At}.$$

(10)

Then, the correlation $C_{II}\left(t\right)$ can be expressed in terms of the matrix elements of $G\left(t\right)$ and of the covariance matrix as follows [39]:

$$C_{II}\left(t\right)=G_{Ix}\left(t\right){\sigma }_{xI}+G_{Iv}\left(t\right){\sigma }_{vI}+G_{II}\left(t\right){\sigma }_{II},$$

(11)

where the symbols $G_{ij}\left(t\right)$ denote the response of the variable $i\in \{x,v,I\}$ to a perturbation on the variable $j\in \{x,v,I\}$. This algebraic computations have been handled with Mathematica. The general form of $C_{II}\left(t\right)$ is an exponential decay, modulated by sinusoidal oscillations.

The autocorrelation from the experimental data is obtained by averaging the whole recorded temporal signal (excluding initial and final transient regimes) on a sliding window of $0.1$ s. We then compare the results from experiments with the correlations obtained analytically and with numerical simulations of the model Equations (1)–(3), as reported in Figure 5, for several values of the load resistance R and for acceleration $a=0.25$ g (other values of a show the same behaviors). First, we observe that analytical predictions and numerical simulations perfectly agree, confirming the correctness of the analytical approach. More interestingly, the comparison with experimental data shows that the model very well reproduces the characteristic frequency of the signal, and the amplitude decay of the oscillations, with small discrepancies for larger values of the load resistance R.

5. Stochastic Energetics

We now analyze the system from a stochastic thermodynamics [7] approach. In particular, we focus on the work fluctuations at finite times, comparing the model predictions obtained from numerical simulations, with experimental data. According to the prescription of stochastic energetics [40], we define the heat exchanged along a trajectory in a time interval $[0,\tau ]$ with the surrounding medium as

$$Q_{ex}\left(\tau \right)=-{\int }_0^{\tau }\gamma v{\left(t\right)}^{2}dt,$$

(12)

and the energy fed into the system from the external driving as the integral of the injected power $P_{inj}=M\xi \left(t\right)v\left(t\right)$

$$E_{inj}\left(\tau \right)=M{\int }_0^{\tau }\xi \left(t\right)v\left(t\right)dt,$$

(13)

where the product is meant according to the Stratonovich definition. Here, we consider these two contributions separately, at variance with standard Brownian systems where the two terms arise from the coupling of the particle with the fluid. Indeed, in our system, thermal fluctuations are completely negligible, whereas $Q_{ex}$ is the heat dissipated to the environment due to the viscosity $\gamma$. The injected power is instead due to the white noise introduced by the (macroscopic) shaker.

Next, using the Langevin Equation (2), we can rewrite these two terms as follows

$$\begin{array}{ccc}\hfill Q_{ex}\left(\tau \right)+E_{inj}\left(\tau \right)& =& {\int }_0^{\tau }\left[-\gamma v{\left(t\right)}^2+Mv\left(t\right)\xi \left(t\right)\right]dt\hfill \\ & =& \frac{1}{2}M[v{\left(\tau \right)}^2-v{\left(0\right)}^2]+\frac{1}{2}{k}_s[x{\left(\tau \right)}^2-x{\left(0\right)}^2]+\theta {\int }_0^{\tau }v\left(t\right)I\left(t\right)dt\hfill \\ & =& \Delta E+W^{\prime }\left(\tau \right),\hfill \end{array}$$

(14)

where

$$\Delta E=\frac{1}{2}M[v{\left(\tau \right)}^2-v{\left(0\right)}^2]+\frac{1}{2}{k}_s[x{\left(\tau \right)}^2-x{\left(0\right)}^2]$$

(15)

is the mechanical energy variation, and

$$W^{\prime }\left(\tau \right)=\theta {\int }_0^{\tau }v\left(t\right)I\left(t\right)dt$$

(16)

can be interpreted as the electromechanical work, which allows us to read Equation (14) as the first law for stochastic thermodynamic quantities. The effective work extracted by the harvester is defined as the difference between $W^{\prime }\left(\tau \right)$ and the heat dissipated in the parasite resistance $R_C$

$$W\left(\tau \right)=W^{\prime }\left(\tau \right)-R_C{\int }_0^{\tau }I{\left(t\right)}^{2}dt.$$

(17)

Next, exploiting Equation (3), the work W can be related to the dissipated heat in the load resistance R in the time interval $[0,\tau ]$

$$Q_{diss}\left(\tau \right)=R{\int }_0^{\tau }I{\left(t\right)}^{2}dt,$$

(18)

via the relation

$$W\left(\tau \right)=Q_{diss}\left(\tau \right)+\frac{L_c}{2}[I{\left(\tau \right)}^2-I{\left(0\right)}^2].$$

(19)

The difference between the two quantities is given by a term non-extensive in time, that can be neglected for long trajectories.

In our experimental setup, we have direct access to the current I in the load circuit and, thus, to the work W, according to Equations (18) and (19). Therefore, we can compare the behavior of the fluctuations of these quantities with those obtained from the numerical simulations of the linear model, where the parameters have been fixed through the fitting procedure described above. Let us note that here, the extracted power is directly related to the dissipated heat by the Joule effect in the load resistance. Therefore, by definition, since it involves a square of the electrical current, the work is always positive. This rules out the possibility of observing negative fluctuations, even at small times, and thus the validity of a work fluctuation relation cannot be explored.

In Figure 6 we compare the power fluctuations measured in the experiments for different times and several parameters, with the results of numerical simulations of the linear model. A very good agreement is found in all cases. This shows that the effective linear model provides an accurate description of the real system, not only for average values, but even at the fine level of fluctuations.

Let us also comment on the functional forms of the work distributions at finite times. We observe that, before converging towards a Gaussian distribution, the probability density function shows very long tails at large values of $W/\tau$. This means that, over small time intervals, there is a finite probability to extract a work much larger than the average value even if the most probable value is smaller than the average.

Finally, we observe that the efficiency of the system can be defined as the ratio $\eta =P_{harv}/\langle P_{inj}\rangle$, where, from Equation (13), $\langle P_{inj}\rangle =M{D}_0$. Due to the linearity of the model, the dependence on $D_0$ appears as a prefactor in both the expressions and thus cancels out in the ratio. This results in a behavior of $\eta$ as a function of the load resistance R similar to the one reported in Figure 4, with a peak around $R=500\mathsf{\Omega }$. As a function of the other parameters, as expected, we find that the efficiency increases with the electromechanical coupling $\theta$ and decreases with the parasite resistance $R_c$.

6. Conclusions

We have shown that a linear model of an underdamped Langevin equation with an electromechanical coupling can very well reproduce the behavior of a real electromagnetic energy harvester, not only with respect to the average values of the measured current, but also with respect to its temporal correlations and fluctuations. The analysis of these quantities and the comparison between analytical prediction and experimental data allow one to obtain a better estimate of the effective parameters for the model.

The study of the stochastic energetics of the system could be extended to other relevant quantities, such as entropy production. In order to address this issue from an experimental perspective, however, a more complete access to the dynamical variables is necessary. In particular, this requires the measurements of the position and velocity of the mass, which are not directly accessible with the actual experimental setup. Indeed, for instance, as illustrated in Refs. [32,33] for a similar linear model, the contribution of the medium entropy production takes the form $\Delta s_m=-1/T{\int }_0^{\tau }dt[M\dot{v}\left(t\right)+\theta \tilde{I}(x,v)+k_{s}x\left(t\right)]v\left(t\right)$, where T is related to the noise amplitude $D_0$ and $\tilde{I}(x,v)$ is obtained from a projection over the current distribution probability [33] and explicitly involves the variables x and v. We plan to address the study of entropy production in future works, in particular for a piezoelectric energy harvester [27], where direct access to all the dynamical variables is possible, and therefore, an experimental measurement of entropy production can be compared to theory and simulations.