Stochastic thermodynamics of a piezoelectric energy harvester model

We experimentally study a piezoelectric energy harvester driven by broadband random vibrations. We show that a linear model, consisting of an underdamped Langevin equation for the dynamics of the tip mass, electromechanically coupled with a capacitor and a load resistor, can accurately describe the experimental data. In particular, the theoretical model allows us to define fluctuating currents and to study the stochastic thermodynamics of the system, with focus on the distribution of the extracted work over different time intervals. Our analytical and numerical analysis of the linear model is successfully compared to the experiments


I. INTRODUCTION
From microscopic organisms in the biosphere, life in general and human activities in particular critically depend on the conversion of different forms of energy into useful work. Harvesting energy from the environment is therefore a central task in many applications, where random fluctuations possibly arising from disparate sources at different scales, from the microscopic thermal Brownian motion in a fluid, to the macroscopic vibrations in means of transport, can be converted into work.
The well established rules of thermodynamics for macroscopic systems become blurred when fluctuations are relevant and have to be taken into account [1]. From a theoretical perspective, the study of fluctuations is addressed within the theory of stochastic thermodynamics, where the standard concepts of energy, heat, work and entropy, are extended to non-equilibrium systems, driven by external forces or coupled to different reservoirs. In this framework, the interest is focused on the fluctuations of the above quantities defined along a single trajectory in the stochastic motion of the system and on their probability distributions. Indeed, general relations have pushed the range of validity of standard thermodynamics into the realm of non-equilibrium regimes [2][3][4]: from the Fluctuation Relations [5][6][7][8][9] and the generalized fluctuation-dissipation relations [10,11], to the general results ruling work and heat exchanged in non-equilibrium transformations, such as the Jarzinski relation [12], the Crooks fluctuation theorem [13], or the Hatano-Sasa relation [14]. Very recently, thermodynamic uncertainty relations bounding the signal to noise ratio of a measured current have been also discovered [15]. In particular, the study of work and heat fluctuations has been the object of focus in several systems, such as overdamped linear Langevin Equation [16], particle diffusion in time-dependent potentials [17][18][19][20][21][22], Brownian particles driven by correlated forces [23], general thermal systems [24], asymmetric processes [25], underdamped Langevin Equation [26], or in transient relaxation dynamics [27]. The interest in these quantities is motivated by the search for optimization protocols in models of stochastic engines or, from a more theoretical perspective, by the general symmetry properties or by singular behaviors that work and heat distributions can show [28]. Experimental studies confirming theoretical predictions have been reported for instance in [29][30][31][32].
Energy harvesting model systems represent an interesting context where the concepts of stochastic thermodynamics can be applied, due to the fundamental relevance of random fluctuations. Very well studied examples are the Brownian (or molecular) motors [33], also known as ratchet models, where a probe is in contact with a thermal bath and the presence of a spatial asymmetry coupled to some non-equilibrium source allows to rectify the motion of the probe, with extraction of useful work. These systems have been studied theoretically and experimentally for instance in the context of granular media [34,35], where the dissipative interactions among grains induce the non-equilibrium behavior, or in biological motors [36], where active internal forces are at play.
More application studies have been carried on in different kinds of energy harvesters, that are based on the piezoelectric properties of some materials. In this case, macroscopic vibrations of the system, as for instance in a car or in a train, can induce small currents, from which an output power can be extracted to feed sensors or small electrical devices [37]. Typically, piezoelectric harvesters are employed in resonant cantilever structures (see Figure 1). The mechanical to electrical energy conversion mechanism is based on the piezoelectric effect that is the ability of some materials (notably crystals and certain ceramics) to generate an electric voltage in response to an applied mechanical stress. Due to their resonant nature, piezoelectric harvesters are typically studied in steady state sinusoidal conditions at frequencies belonging to their resonance band [38,39]. In particular, the main focus is on their energetic performance, that is on the mechanical and power electronic architectures and on the control techniques leading to the maximization of the extracted power [40,41]. Less attention has been devoted to the case of resonant piezoelectric harvesters excited by non-sinusoidal vibrations or solicited by white noise vibrations [42][43][44][45]. In particular, the theoretical analysis of [42] provided a stochastic description of the output power from resonant energy harvesters driven by broadband vibrations and output power dependence on signal bandwidth was considered. Instead, Ref. [44] proposed a methodology for the probabilistic analysis of a cantilever piezoelectric harvester under white Gaussian noise, without experimental validation. In [45], an experimental analysis on piezoelectric harvesters is carried out in the presence of harmonic, random and sine on random vibrations with particular reference to the electrical power extraction. However, in all previous studies on these energy harvesters no attention was devoted to the analysis of the fluctuations and distributions of relevant quantities such as the extracted power, in the general framework of stochastic thermodynamics of non-equilibrium systems.
Here we consider a typical piezoelectric harvester in a resonant cantilever structure driven by random broadband vibrations. Despite the several nonlinearities present in the experimental system, we show that a linear model with effective parameters can well reproduce the observed dynamics. In particular, we consider a mass in the presence of a harmonic potential, in contact with a source of white noise. The mass is also electromechanically coupled with a capacitor, which allow for power extraction through a load resistance. First, we show that the characteristic response curve, output power vs. load resistance, obtained from experiments is very well fitted by the analytical formula derived for the theoretical model. Then, we define a fluctuating work along a system trajectory, according to the standard approach of stochastic thermodynamics, and we focus on the study of the work fluctuations. We find that also the distributions of the work measured experimentally over different time intervals, are in very good agreement with those computed in numerical simulations of the linear model, using effective parameters.
Our study presents an experimental characterization of the distributions of integrated currents (output power) for a system that is used in applications as a valid energy harvester device. Moreover, our analysis shows that a simplified linear model, which allows for analytical computations, is able to accurately reproduce the experimental results, even at the fine level of fluctuations.

II. EXPERIMENTAL SETUP
The schematic representation of the experimental setup for the piezoelectric harvester is shown in Figure 2. It consists of a cantilever structure with a tip mass, whose displacement in time x(t) due to the vibrations provided by the shaker, induces a voltage v p (t) across the electrical load.
In all the experimental tests that we have carried out, we have used the commercial piezoelectric harvester MIDE PPA-4011 loaded by different electrical load resistances. This product incorporates four piezoelectric wafers resulting in significant performance improvements with respect to other models by MIDE. The harvester has been mounted on a shaker by using a support providing the possibility of different clamping positions. The addition of tip masses in order to define mechanical properties of the resonant structure is also possible. A picture of the whole experimental setup is shown in Figure 2. The shaker VT-500 by Sentek (500 N rated force and 450 m/s 2 maximum acceleration) has been used to get the desired input vibrations. The controller Spider 81 allowed application of the desired voltage signal to the shaker amplifier and to carry out the recording, by means of its built-in acquisition board, of the voltage across the load resistance. An accelerometer 3055D2 by Dytran (sensitivity 100 mV/g on the range 50 g) has been used to monitor the applied acceleration in order to implement a closed-loop feedback vibration control. In order to study the response of the system to broadband vibrations, we fed the shaker with a Gaussian signal generated with standard software (MATLAB), with a sampling rate f = 5 kHz. In Figure 3 typical waveforms of the input acceleration and of the voltage across the load resistance (for R = 2200 Ω) recorded during experimental tests are shown.
The first quantity we analyzed is the extracted power P harv that can be obtained from the average dissipated heat on the load resistance for unit time, P harv = v 2 p /R. We show in Figure 4 the characteristic curves P harv vs. R, for different values of the input acceleration. From our analysis, we find an optimal resistance value R * ∼ 2000 Ω, which is independent of the shaking amplitude.

III. THEORETICAL MODEL
In order to describe the observed experimental results and to extend the study of the system to fluctuating quantities relevant in the stochastic thermodynamics framework, we consider the following linear model where ξ is white noise with zero mean and correlation ξ(t)ξ(t ′ ) = 2D 0 δ(t − t ′ ). In the above equations, x represents the displacement of the tip mass M , v its velocity, γ the viscous damping due to the air friction, K s the stiffness of the cantilever in the harmonic approximation, v p the voltage across the load resistance R, C p the effective capacitance in the circuit, and θ the effective electromechanical coupling factor of the transducer. In Equation (2) we have neglected the thermal fluctuations on the tip mass, which are too small to affect its motion.
In the system we can identify several characteristic times: τ 1 = M/γ, τ 2 = M/K s , τ 3 = C p γ/θ 2 , τ 4 = C p R. τ 1 is the relaxation time of the tip mass due to the viscous damping, τ 2 is the relaxation time within the harmonic potential, τ 3 represents the typical timescale of the coupling between the proof mass and the capacitor, τ 4 is the characteristic time of the RC circuit. Among the various characteristic times, τ 3 is the only one depending on quantities belonging to both the electrical and the mechanical subdomain of the whole harvesting system. Hence, its physical meaning is not as intuitive as in the case of the other characteristic times. In any case, in order to better highlight its role, it is possible to show [41] that there is a link between the amplitude of the speed of the tip mass and the amplitude of the external acceleration, when the harvester operates in sinusoidal conditions at the resonance frequency and in open circuit (no load, R → ∞). In particular, one has the relation [41] v max = a max θ M Therefore, the speed amplitude, in the above operating conditions, depends on all the three characteristic times, τ 1 , τ 2 , and τ 3 , that assume finite values. It does not depend on τ 4 , since it is unbounded in open circuit (R → ∞).
where the memory kernel Γ(t) has the simple exponential form Let us note that here the friction memory kernel Γ(t) is not associated with any noise term. This puts the system by construction in a non-equilibrium state, because the fluctuation-dissipation relation of the second kind does not hold.

A. Average Values
The static properties of the linear model can be obtained by standard methods [46]. We introduce the column vector X = (x, v, v p ) T and the coupling matrix so that the Equations (1)-(3) can be rewritten in vectorial form aṡ where η = (0, ξ, 0) T . Defining the covariance matrix σ = X T X as at stationarity one has the constraint where D is the noise matrix From Equation (10) one gets the following relations for the covariance matrix elements The stationary distribution is a multivariate Gaussian where σ −1 denotes the inverse matrix of σ. The explicit expressions of the elements of σ are reported in Appendix V. The average output power is the heat dissipated into the resistance per unit time and its explicit expression as a function of the parameters is

B. Fitting the Model to Experimental Data
The linear model described by the Equations (1)-(3) contains several physical parameters that are directly controlled in the experiments and others that can be fitted to match the measured data. In particular, the parameter D 0 that quantifies the amplitude of the white noise is related to the acceleration a provided by the shaker and to the sampling rate 1/∆t of the input signal, D 0 = a 2 ∆t/2, where ∆t = 1/f = 0.0002 s. Furthermore, the parameters K s and M are related to the characteristic frequency of the system, which for the experimental apparatus is K s /M = 2π × 140 Hz. Finally, the capacitor C p is estimated as C p ∼490 nF. The other parameters can be fitted to the experimental data through the analytical expression (20) as a function of the load resistance R. For the case a = 9.81 m/s 2 , we obtain the following values for the model parameters: M = 0.0083 ± 0.0002 Kg, θ = 0.0195 ± 0.002 N/V, γ = 0.359 ± 0.05 Kg/s. This set of parameters is used also for other values of the shaker accelerations used in the experiments, a = 0.8 × 9.81 m/s 2 and a = 1.2 × 9.81 m/s 2 , providing a very good agreement between analytical predictions and experimental data, as shown in Figure 4.

C. Stochastic Energetics
The theoretical model allows us to study fluctuations and distributions of thermodynamic quantities defined at the level of the single trajectory. In particular, according to Sekimoto [47], we define the heat exchanged along a trajectory in a time interval τ with the surrounding medium as and the energy fed into the system from the external driving as the integral of the injected power P inj = M ξ(t)v(t) The product in the above equation is meant according to the Stratonovich prescription. Note that the heat in Equation (21) is released toward the medium. Using the Langevin Equation (2), we can rewrite these two terms as follows where is the mechanical energy variation and can be interpreted as the work performed by the harvester. Equation (23) represents the first principle for stochastic thermodynamic quantities. The average output power in the stationary state is then where vξ = D 0 and the last equality follows from Equation (23) using the stationary result ∆E = 0. The output power can be also expressed as the dissipated heat in the load resistance in the time interval τ which, exploiting Equation (3), can be related to the fluctuating work W by This shows that the difference between the fluctuations of W and Q diss is a term non-extensive in time. In our system, we have numerically checked that, for the studied time intervals, these border terms can be neglected. However, more generally, there are systems where such terms can be relevant for fluctuations, see for instance [17,48].
We have studied numerically and experimentally the work W (τ ) for different values of τ . The experimental evaluation of this quantity has been obtained integrating the time series of the output signal for the voltage v p (t) (the whole recorded time series were 600 s long). The results are reported in Figure 5. First, we stress the good agreement between the work distributions obtained from experiments and numerical simulations. This shows that the linear model with white noise is able to accurately describe the experimental system, even at the level of fluctuations, in the range of explored parameters. The work distributions present a pronounced asymmetry for small time intervals τ , characterized by an exponential tail for large values of W/τ , as also observed for the functional form of the injected power distribution in the overdamped Langevin equation, obtained analytically in [16]. At large times, the distributions seem to converge towards a Gaussian form, symmetric around the mean. In order to analyze the behavior of the work distributions measured in experiments as a function of the load resistance we have fitted the tail of the curves at short times (τ = 0.01 s) with an exponential function f (x) ∼ exp(−x/α), extracting the parameter α. For the curves at large times (τ = 1 s) we have fitted the data with a Gaussian function g(x) ∼ exp(−(x − µ) 2 /2σ 2 ) to obtain the variance σ 2 as a function of R. The results are shown in Figure 6. We observe that both α and σ 2 have a non-monotonic behavior, with a maximum appearing around the value which maximizes the mean extracted power, showing that fluctuations are larger in proximity of the optimal working point.

IV. CONCLUSIONS
We have studied experimentally a piezoelectric energy harvester driven by random broadband vibrations, focusing on the behavior of the extracted power as a function of the load resistance and of the vibration amplitude. We have shown that a linear model, consisting of an underdamped Langevin equation for the tip mass coupled with voltage dynamics reproduces very well the experimental data.
Moreover, the theoretical model allowed us to address the issue related to the behavior of the fluctuations of nonequilibrium currents, such as the extracted work in a time interval. This analysis plays a central role in the context of stochastic thermodynamics, where the properties of the system are scrutinized at the level of single trajectories. The comparison between the results of numerical simulations of the model and the experimental data showed that the linear system of equations provides a good approximation of the real system, reproducing the same behavior of the work distributions as a function of the parameters.
Our findings represent the first experimental study of the work fluctuations in a piezoelectric energy harvester and show that the observed behaviors can be consistently rationalized within a simple model, paving the way to future analyses. In particular, from the theoretical perspective, due to the linear nature of the model, the analytical computation of the work distribution and its large deviations function could be obtained with a path integral approach, as for instance described in [16]. Moreover, it could be interesting to modify the model by adding a noise source also in the equation for the voltage, obtaining a system of two coupled Langevin equations, or considering a bistable potential for the tip mass, as proposed in [49,50]. Analyses of other quantities such as heat or entropy production could be also carried out along the same lines. On the experimental side, it could be interesting to perform a similar study of current fluctuations in a system driven by realistic vibration sources, like cars or trains, taken from available databases, or where the simple linear resistance load is replaced by a diode bridge, as often considered in applications.