Magnetic Frequency Tuning of a Shape Memory Alloy Thermoelectric Vibration Energy Harvester

1. Introduction

Energy storage from oscillating sources is a promising solution for powering low-power electronic devices [1,2,3,4]. Among the many sources of untapped ambient energy, mechanical vibration, light and thermal energy are the most widely used [5,6,7]. The high response speed, on which the efficiency of the energy provided by the mechanism depends, is a challenge for actuators based on SMA elements. The time to reach the phase transitions of these harvesters can be reduced. In other words, the response rate is limited by the heating time and especially the cooling time of the SMA element actuator. Therefore, reducing the time to reach the phase transition temperature is essential to improve the reaction rate. One of the approaches to increase the actuation speed of SMA is to reduce the dimensions. Triggering rates of up to 1.6 kHz have been reported when the thickness of the SMA element is reduced to close to 1 μm. These are several orders of magnitude higher than the speeds typically associated with SMA actuators [8]. However, the size of the SMA is not the only condition that determines the actuation speed. For some actuators, the thermal response may be slower or faster due to changes in the cooling or heating behaviour of the actuator. This means that even if actuators have similar SMA thicknesses, the actuation speeds would be different depending on the equivalent mass of the entire actuation system or the amount of supplied energy required for actuation [9,10]. Changing the chemical composition of the SMA in the actuator allows cyclic transformation of the alloy over a relatively lower temperature range. The addition of traces of Ti to the Cu-Zn-Al alloy resulted in a significant reduction in transformation phase temperatures [11].

The efficiency of energy harvesting in piezoelectric cantilever (PEC) devices depends on the magnitude of the frequency at which the mechanism oscillates [12,13]. To increase the energy harvesting efficiency, many methods have been proposed to optimise the energy harvesting by adjusting the oscillation frequency [14,15]. One solution to controlling the oscillation frequency is to use multiple piezoelectric cantilevers whose mechanical resonant frequencies are distributed across the desired oscillation frequency in such a way that, at any given frequency, at least one circuit operates close to its resonant frequency. The most general approach to the design of such a mechanism is given in [16,17]. The frequency configuration of the PEC can also be achieved by varying the magnitude of the effective mass of the beam [18] or the stiffness of the PEC [19]. The stiffness of the PEC can be varied by changing the cantilever geometry or the material properties. The solution development includes the properties of the PECs, their elements and the properties of the electrical circuit in the calculations [20]. Methods for frequency scaling [21,22], implementation of non-linear methods [23,24,25] have also been proposed.

Bouhedma at al have presented a piezoelectric vibrational energy harvester whose frequency characteristics can be altered using magnets with up to 18% change [26]. To achieve high electrical power stability, a concept of piezoelectric energy harvester that mimics the switch mechanism has been proposed. The extracted energy is increased by a collision between two magnets [27]. The dynamics of a bistable magnetic oscillator is controlled by the strength and polarization of the added magnets [28]. The extracted energy voltage reaches the highest value in nonlinear oscillations, with the output power improved by up to 600% compared to linear oscillations. The functionality of vibrational energy harvesters using magnets is reviewed, and monostable and bistable configurations introduced by magnets with different excitation levels are discussed to improve the efficiency of the harvested energy output [29]. The dynamic characteristics of broadband tristable magnetic coupling harvesters have been investigated [30]. Compared with the bistable state, the results show that this harvester easily transitions through its different states, providing high output energy over a wider frequency range. A broadband nonlinear quin-stable energy harvester is proposed [31]. To control the generation of piezoelectric energy, magnets are usually deployed in repulsive or attractive configuration [32]. This hybrid approach increases the power generation capabilities of the harvester due to the combined piezoelectric and electromagnetic operation under the same external excitation [33]. Furthermore, a new generation of devices can be developed using materials and structures based on magnetoelectric effects. Magnetoelectric effects in ferromagnetic-ferroelectric layered composites arise due to magnetostriction and piezoelectric effect during the ferroelectric phases. The important advantage of magnetic field-based methods compared to mechanical methods is their remoteness and no contact. Ferrite composites have potential applications in sensors, autonomous power sources, electrically controlled microwave devices, and others [34]. Several reviews on piezoelectric energy harvesting have been published, most of them mainly focusing on biomedical applications [35], microelectromechanical systems (MEMS) [36], nanogenerators [37], materials, and energy harvesting from various sources [38].

In previous publications, the authors of this paper have developed and investigated a self-excited thermomechanical oscillator [39] and an energy harvester[40], in which the heat of a constant-energy heater is converted into vibrations of an SMA filament. Two masses are mounted on the SMA filament, which perform vibration two perpendicular directions, longitudinal and transverse. The longitudinal vibration is due to the contraction of the SMA filament when it is close to the heater. These contractions shorten the filament, and this causes it to move away from the heater and subsequently cool with elongation due to the inertial forces of the masses. The transverse vibrations are due to the relative movement of the filament from the SMA in the vertical direction and are caused by inertial forces and the weight of the masses in addition to the contractions. The variation of the Young’s modulus due to the phase changes of the SMA cause another oscillation of higher frequency which is summed with the longitudinal oscillations. The energy harvester [40]implemented on this principle transforms the resulting oscillations, through piezoelectric cantilevers mounted at the ends of the filament. Studies have shown that the vibrations have a pseudo-stochastic character due to the hysteresis of the relative martensitic fraction.

The aim of this paper is to investigate the possibilities of increasing the generated power of the above-described self-excited thermal energy harvester by tuning the natural frequency of the vibrating system using permanent magnets.

2. Design Concept of Energy Harvester and Natural Frequency Tuning System

The thermoelectric energy harvester consists of an SMA filament that is fixed between two fixed supports O and C, symmetrically arranged relative to a fixed heater. Two identical spheres are fixedly mounted on the filament near the heater, on which magnetic forces due to symmetrically placed NdFeB magnets (M₁ and M₂) act. Two piezoelectric cantilevers ${OP}_1$and ${OP}_1$are positioned at the fixed ends of the SMA filament, which are mounted, using a heat-resistant foil, to follow the bending of the filament. The whole structure is symmetric about the y ordinate of the fixed coordinate system O₁xy (Figure 1). The filament material is a shape memory alloy of nitinol known by the trademark Flexinol^® [41]. For the piezoelectric cantilevers, two vibration sensors type DT1-052K [42] were used, which are constructed of a thin PVDF piezoelectric film metallized with silver on both sides representing flexible electrodes covered with plastic protective layers and galvanically connected with solder terminals.

The principle of operation of the energy harvester, as already mentioned above, is based on temperature changes above the heater, which are high near its surface and cause contractions of the SMA filament and its simultaneous moving away from the warm zone. Once in the lower temperature zone, the SMA filament softens and stretches under the inertial forces of the masses, bringing it closer to the heater and the action zone of the magnets. These actions are repeated continuously, resulting in a cyclic vibration process of the SMA filament and mass system. Near the supports, the vibrational bends of the SMA filament are transmitted to the piezoelectric cantilevers and they generate electrical energy (see Video from Supplementary Material).

3. Composing the Dynamical Model by Combining the Interactions of the Different Physical Domains

The distance between the two supports O and C can be adjusted, which creates possibilities to mount the SMA filament with pre-tension or pre-sag. The geometric parameter designations are in accordance with Figure 1 and the approach for determining the elastic and temperature deformations of SMA filament and PVDF cantilevers is similar to the approach described in [40].

3.1. Determination of Elsatic and Temperature Parameters of the SMA Filament and Bending of the Piezoelectric Cantilevers

Since half of the distance between supports O and C is denoted by l, one can express the pre-strain Δls₀ of this parameter

$$\Delta l_{s0}=l-l_s.$$

(1)

where l_s is for half the length of the SMA filament before tensioning.

A relative length ${\lambda }_0$ is introduced as

$${\lambda }_0={\displaystyle \frac{l_{01}}{l}}={\displaystyle \frac{l_{03}}{l}}.$$

(2)

where $l_{10}=O{A}_0$, $l_{30}=B_{0}C$ are the lengths of the filament sections when the ferromagnetic spheres are on the line OC.

Considering the symmetry and uniformly distributed deformations along the entire length of the filament, then for the end sections OA₀ and B₀C for the pre-strains $\Delta l_{s10}$and $\Delta l_{s30}$, follows

$$\Delta l_{s01}=\Delta l_{s03}=\Delta l_{s0}{\lambda }_0,$$

(3)

In this case for the middle section A₀B₀ its deformation can be written

$$\Delta l_{s012}=2\Delta l_{s0}\left(1-{\lambda }_0\right).$$

(4)

Due to the specific symmetrical design shown in Figure 1, it can be assumed that in the middle section AB the filament moves parallel to the horizontal surface of the heater and therefore heats uniformly along its entire length. Since the end sections are relatively far from the heater they are assumed not to heat up and their temperature is low i.e., they are cold sections with a temperature equal to that of the room. The deformation of the cold (or end) sections due to the displacement of the masses and the initial tension for the symmetrical motion variant has the form

$$\Delta l_1=\sqrt{{\left(x_A-x_O\right)}^2+{\left(y_A-y_O\right)}^2}-l_{10}+\Delta l_{s0}{\lambda }_0$$

(5)

and

$$\Delta l_3=\sqrt{{\left(x_B-x_C\right)}^2+{\left(y_B-y_C\right)}^2}-l_{30}+\Delta l_{s0}{\lambda }_0,$$

(6)

where $x_O$, $x_A$, $x_B$, $x_c$ $y_O$, $y_A$, $y_B$, $y_c$ are the coordinates of the points O, A, B and C.

Under the same conditions for the deformation of the middle section can be written

$$\Delta l_2=\sqrt{{\left(x_A-x_B\right)}^2+{\left(y_A-y_B\right)}^2}-l_{20}+2\Delta l_{s0}\left(1-{\lambda }_0\right).$$

(7)

Since the end sections of the SMA filament are always cold and it is assumed that the room temperature is always below the final martensitic temperature M_f of NiTi, the stiffness of these sections can be expressed by

$$k_l={\displaystyle \frac{E_{m}A}{l_{01}}},$$

(8)

Where $E_m$ is the NiTi Young’s modulus in the martensitic state, and A is the cross-sectional area of the SMA filament.

For the stiffness of the middle section, it can be written

$$k_{Т}={\displaystyle \frac{E_{02}A}{l_{20}}},$$

(9)

but since this section is subjected to varying temperatures that cause all phase transformations in NiTi the Young’s modulus$E_{02}$also transforms depending on the current phase. According to the model discussed in [43], this module can be represented by the piecewise function

$$E_{02}=\left\{\hspace{0.17em}\begin{array}{c}\begin{array}{l}{E}_a-\left(E_a-E_m\right)R_m\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le {\epsilon }_{s2}\le {\epsilon }_m^y\\ E_a-\left(E_a-E_T\right)R_m\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_m^y\le {\epsilon }_{s2}\le {\epsilon }_m^d\\ E_a-\left(E_a-E_d\right)R_m\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_m^d\le {\epsilon }_{s2}\end{array}\end{array}\right.$$

(10)

where $E_m$, $E_T$, $E_d$, and $E_a$ are Young’s modulus of the fully twined, partially twined, detwinned martensite and austenite respectively; ${\epsilon }_m^y$ is the yield strain of the twined martensite; ${\epsilon }_m^d$ is the minimum strain of the detwinned martensite, and

$${\epsilon }_{s2}={\displaystyle \frac{\Delta l_2}{l_{02}}}$$

(11)

is the strain of the middle section.

The relative martensite fraction $R_m$ is expressed in the form

$$R_m={\displaystyle \frac{S_R}{2}}\left(1+\tanh {\displaystyle \frac{T_m-T}{{\lambda }_T}}\right)+R_C{S}_R$$

(12)

where $T$ is the current temperature of the NiTi, the temperature $T_m$ depends on whether the thermal process is heating or cooling and is obtained by the relation

$$T_m={\displaystyle \frac{1}{2}}\left\{\begin{array}{c}{A}_S+A_f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{dT}{dt}}>0\\ M_s+M_f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{dT}{dt}}<0\end{array}\right.$$

(13)

The unit of coefficient ${\lambda }_T$ is temperature, calculated by the expression

$${\lambda }_T={\displaystyle \frac{1}{2}}\left\{\begin{array}{c}{\displaystyle \frac{A_S+A_f-b_a}{c_a}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{dT}{dt}}>0\\ {\displaystyle \frac{M_s+M_f-b_m}{c_m}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{dT}{dt}}<0\end{array}\right.$$

(14)

where A_s, A_f, M_s, and M_f are the start and finish temperatures for austenite and martensite phases of NiTi. The limiting values of the coefficients, $b_a<0.5({А}_f-A_s)$ and $b_m<0.5(M_s-M_f)$, have been established experimentally [44].

The method for hysteresis modeling is described in detail in [43]. It is based on the analysis of temperature fluctuation point $T_f$ at which gradient change of a transient temperature occurs.

If $T_f$ is in the transient periods with both martensite and austenite fraction in the NiTi exist. Let denote the derivatives of the transient temperature versus time before fluctuation point with ${\dot{T}}_{f-}={\displaystyle \frac{dT}{dt}}\left|T\to T_{f^{-}}\right.$ and after fluctuation point with ${\dot{T}}_{f+}={\displaystyle \frac{dT}{dt}}\left|T\to T_{f^{+}}\right.$.

If the condition

$$\hspace{0.17em}{\dot{T}}_{f-}>0\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}{\dot{T}}_{f+}<0\hspace{0.17em}\hspace{0.17em}$$

(15)

is satisfied, then the so-called scale factor $S_R$ is

$$S_R=\left\{\begin{array}{c}1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_S>T_f>M_f\\ R_m^f=\hspace{0.17em}{\displaystyle \frac{1}{2}}\left(1+\tanh {\displaystyle \frac{T_m-T_f}{{\lambda }_T}}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{A}_f>T_f>A_S\end{array}\right.$$

(16)

and

$$R_C=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_S>T_f>M_f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}or\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{A}_f>T_f>A_S$$

(17)

If at the fluctuation temperature $T_f$, the NiTi temperature changes from cooling to warming, i.e.,

$$\hspace{0.17em}{\dot{T}}_{f-}<0\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}{\dot{T}}_{f+}>0\hspace{0.17em}\hspace{0.17em}$$

(18)

The scale factor is calculated by

$$S_R=\left\{\begin{array}{c}\hspace{0.17em}1\hspace{0.17em}-\hspace{0.17em}{R}_m^f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_S>T_f>M_f\\ R_m^f=\hspace{0.17em}{\displaystyle \frac{1}{2}}\left(1+\tanh {\displaystyle \frac{T_m-T_f}{{\lambda }_T}}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_s\ge T_f\ge M_f\end{array}\right.$$

(19)

and

$$R_C=\left\{\begin{array}{c}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_f>T_f>A_f\\ 1-R_m^f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_s>T_f>M_f\end{array}\right.$$

(20)

Similarly, sub-minor hystereses are described, which occur in the case that already a minor hysteresis appears, and a sub-minor fluctuation temperature occurs in it during a transient. For these cases, analogous approaches described in [43]are used.

To simplify the general formulation of the problem, it is assumed that the temperature field above the heater in which the NiTi filament moves is stationary. In addition, the simplifying assumption that the NiTi filament itself obtains a potential temperature from this field by dependence

$$T=\left\{\begin{array}{c}\left(T_h-T_0\right)e^{-{\left({\displaystyle \frac{y-y_e}{\theta }}\right)}^2}+T_0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}y<y_e\\ T_h\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}{y}_e\le y\le y_f\\ \left(T_h-T_0\right)e^{-{\left({\displaystyle \frac{y-y_f}{\theta }}\right)}^2}+T_0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}y>y_f\end{array}\right.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}$$

(21)

where coordinates $y_e,\hspace{0.17em}{y}_f$ are determined experimentally [45] and they depend on the heater position, $T_h$ is the maximum temperature received by the filament from the heater and $\theta$ is a conditional length, and$T_0$is the temperature of the room.

Considering (21) and the hysteresis theory of the relative martensitic fraction (12)- (20), it is seen that the Young’s modulus ${Е}_{02}$ for the middle section is a potential function of the coordinates of the ferromagnetic masses and the temperature of the SMA filament, which is here also assumed to be a potential function. These properties facilitate the application of Lagrange equations of the second kind to the solution of the dynamics problem.

The bending of the piezoelectric beams due to the movement of the SMA filament is assumed to be as in [40] where, since the bent shape for the left cantilever is assumed to be a parabola of the form

$$y_l=a_l{\left(x+l\right)}^2=a_l{x}^2+2a_{l}lx+a_l{l}^{2}x\in \left[-l,-l+l_p\right]$$

(22)

and for the right cantilever the bent shape is described with the parabola

$$y_r=a_r{\left(x-l\right)}^2=a_r{x}^2-2a_{r}lx+a_r{l}^{2}x\in \left[l-l_p,l\right]$$

(23)

where $l_p$ is the length of the piezoelectric cantilever, $a_l$ and $a_r$ are factors that are determined assuming that the angle of inclination at the end of the cantilevers coincides with that of the SMA filament at the same point (Figure 2).

In accordance with the assumptions on the shape of the bending, the factors are found

$$a_l={\displaystyle \frac{y_A}{l_p\left(2x_A-l_p+2l\right)}}$$

(24)

$$a_r={\displaystyle \frac{y_B}{l_p\left(-2x_B-l_p+2l\right)}}$$

(25)

3.2. Determination of Magnetic Forces and Their Potential Energy

It is assumed that the magnetic force between a spherical ferromagnetic mass and a permanent magnet as depicted in Figure 3(a) is approximated by the function

$$F_m={\displaystyle \frac{\kappa }{{\left({\delta }_0+\delta \right)}^2}}$$

(26)

where $\kappa$ is called the imaginary magnetic mass, and ${\delta }_0$ is the imaginary initial gap and ${\delta }_0$ is the distance between the surface of the magnet and the center of the sphere. The same force function was determined in [46] based on numerical approximations and experiments.

The force function (26), also called the static magnetic force characteristic, has the graphical interpretation shown in Figure 3(b). The magnetic force at a conditional zero gap (it is called conditional because the distance between the sphere and the magnet can practically never be zero because of the roughness and deviations from the shapes of the two surfaces) at ${\delta }_0$=0 is the initial magnetic force

$$F_{m0}={\displaystyle \frac{\kappa }{{\delta }_0^2}}$$

(27)

from where the unit of the imaginary magnetic mass $\kappa$ [Nm²] is evident.

The schematic experimental setup shown in Figure 4 was used to determine the force exerted by the M₁ magnet on one of the spheres and its potential energy. A current position of the ferromagnetic mass at point A, located at position $y_A$. The influence of the displacement of the masses along the axis on the magnetic forces is neglected.

Taking into account that

$${\delta }_1=y_1+y_A$$

(28)

according to (26) the left magnet M₁ will exert the force

$$F_{ml}={\displaystyle \frac{-\kappa }{{({\delta }_0+y_1+y_A)}^2}}$$

(29)

The right magnet M2 will exert the same force because the coordinates of the two ferromagnetic masses are equal.

The potential energies of the two magnet forces are determined after integration

$$E_m=2{\displaystyle \int \frac{-\kappa }{{({\delta }_0+y_1+y_A)}^2}d}{y}_A={\displaystyle \frac{2\kappa }{{\delta }_0+y_1+y_A}}$$

(30)

3.3. Derivation of the Lagrange Equations of the Second Kind

In the energy harvester investigated here, the two ferromagnetic masses move along the $x$ and $y$ axes and, in addition, the two piezoelectric cantilever beams generate two different electric charges $q_l$of the left cantilever and $q_r$ of the right cantilever. This means that the generalized coordinates are six in total. Four of them the coordinates$x_{А}$, $x_B$, $y_{А}$, and $y_B$are mechanical and the two charges$q_l$ and $q_r$ are electrical.

Since the energy harvester is assumed to be in its symmetrical configuration, the following notations can be introduced $x_A=x$, $x_B=-x$, $y_A=y$, $y_B=y$, $q_l=q$, and$q_r=q$. Note that symmetry reduces the generalised coordinates from six to three $x$, $y$ and $q$. The symmetry also allows the following notations to be introduced $x_{О}=-l$, $x_C=l$, $y_{О}=0$, $y_C=0$, $l_1=l_3$, $\Delta l_1=\Delta l_3$, and $l_2=2(l-l_{10})$.

The kinetic energy $E_K$ of the symmetric mechanical system is obtained in the form

$$E_K=m({\dot{x}}^2+{\dot{y}}^2)$$

(31)

where the dot denotes the derivatives of the coordinates with respect to time, and $m$ is the mass of a ferromagnetic sphere.

The potential energy $E_p$ (except the piezoelectric one) consists of the mechanical potential energy and the magnetic potential energy, therefore it is expressed by

$$E_p=k_l\Delta l_1^2+{\displaystyle \frac{k_T}{2}}\Delta l_2^2+{\displaystyle \frac{2\kappa }{{\delta }_0+y_1+y}}$$

(32)

where

$$\Delta l_1^2=\sqrt{{(x+l)}^2+y^2}+{\lambda }_0{l}_{s0}-l_1$$

(33)

and

$$\Delta l_2=2\left[x+({\lambda }_0+1)l_{s0}-l+l_1\right]$$

(34)

The piezoelectric energy $E_{pe}$, according to [47], is derived in the form

$$E_{pe}={\displaystyle \frac{q^2}{C(1-k_{33}^2)}}-{\displaystyle \frac{2d_{33}{K}_a}{C(1-k_{33}^2)}}q\Delta l_p+{\displaystyle \frac{K_a\Delta l_p^2}{(1-k_{33}^2)}}$$

(35)

where

$$K_a=3I/s^E{l}_p^3$$

(36)

$$\Delta l_p={\displaystyle \frac{y{h}_p}{l_p-2\left(x+l\right)}}$$

(37)

is the lengthen of the piezoelectric cantilever

$$C={\displaystyle \frac{{\epsilon }_{33}w(l_p+\Delta l_p)}{h_p}}$$

(38)

is equivalent capacitance of the piezoelectric transducer, $d_{33}$,${\epsilon }_{33}$ are piezoelectric constant of PVDF and the dielectric permittivity respectively, $s^E$is the compliance of the piezoelectric cantilever under constant electric field, $k_{33}^2$ is the electromechanical coupling coefficient of PVDF, $I$, $w$, and $h_p$ are the second moment of area width, and the thickness of the PVDF cantilever respectively.

Based on the above results, the Lagrangian $\Lambda$ of the magnetoelectromechanical system is written in the form

$$\Lambda =m({\dot{x}}^2+{\dot{y}}^2)-k_l\Delta l_1^2-{\displaystyle \frac{k_T}{2}}\Delta l_2^2-{\displaystyle \frac{2\kappa }{{\delta }_0+y_1+y}}-{\displaystyle \frac{q^2-2d_{33}{K}_{a}q\Delta l_p+C{K}_a\Delta l_p^2}{C(1-k_{33}^2)}}$$

(39)

The system of three coupled equations describing the behavior of the multiphysics energy harvester is derived according to the Lagrange equations, which has the form

$${\displaystyle \frac{d}{dt}}{\displaystyle \frac{\partial \Lambda }{\partial {\dot{q}}_i}}-{\displaystyle \frac{\partial \Lambda }{\partial q_i}}=Q_i,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,2,3$$

(40)

in which the generalized coordinates are $q_1=x$, $q_2=y$ and, $q_3=q$ and the generalized forces $Q_i$ have the form

$$\begin{array}{l}{Q}_1=-{\beta }_x\dot{x}\\ Q_2=-2mg-{\beta }_y\dot{y}\\ Q_3=-R_L\dot{q}\end{array}$$

(41)

Where ${\beta }_x$, ${\beta }_y$ are the viscous resistance coefficients along the axis $x$and $y$, respectively, $g$ is the ground acceleration, $R_L$is the load resistance in the circuit of the piezoelectric cantilevers when connected in parallel.

Applying the equations (40) to (39), the following system is derived:

$$\left|\begin{array}{c}\begin{array}{l}m\ddot{x}+k_l{\displaystyle \frac{\Delta l_1\left(x,y\right)\left(x+l\right)}{\Delta l_1\left(x,y\right)-D_1}}+2k_T\left(y\right)\left(D2+x\right)+{\displaystyle \frac{h_p\left(K_a{d}_{33}{l}_p+q\right)yq}{C_p{K}_p{\left[h_{p}y-2l{l}_p+l_p^2-2l_{p}x\right]}^2}}-\\ -{\displaystyle \frac{2K_a{h}_p^2{y}^2}{K_p{\left(l_p-2x-2l\right)}^3}}=-{\beta }_x\dot{x}\\ \end{array}\\ \begin{array}{l}m\ddot{y}+k_l{\displaystyle \frac{\Delta l_1\left(x,y\right)y}{\Delta l_1\left(x,y\right)-D_1}}+{\displaystyle \frac{d{k}_T\left(y\right)}{dy}}{\left(D2+x\right)}^2-{\displaystyle \frac{\left(2K_a{d}_{33}{l}_p+q\right)\left(l_p-2x-2l\right)h_{p}q}{2C_p{K}_p{\left[h_{p}y-2l{l}_p+l_p^2-2l_{p}x\right]}^2}}+\\ -{\displaystyle \frac{K_a{h}_p^{2}y}{K_p{\left(l_p-2x-2l\right)}^2}}+{\displaystyle \frac{2\kappa }{{\left({\delta }_0+y_1+y\right)}^2}}=-mg-{\beta }_x\dot{y}\\ \\ {\displaystyle \frac{2d_{33}{K}_a{h}_{p}y+2\left(2l-l_p\right)q+4qx}{C_p{K}_p\left(h_{p}y-2l_{p}l+l_p^2-2l_{p}x\right)}}=-R_L\dot{q}\end{array}\end{array}\right.$$

(42)

where

$$C_p={\displaystyle \frac{{\epsilon }_{33}w}{h_p}}$$

(43)

and

$${\displaystyle \frac{d{k}_T\left(y\right)}{dy}}={\displaystyle \frac{d}{dy}}{\displaystyle \frac{E_{2}A}{l_2}}={\displaystyle \frac{A}{l_2}}{\displaystyle \frac{d{E}_{20}}{dy}},$$

(44)

The Young’s modulus derivative ${\displaystyle \frac{d{E}_{20}}{dy}}$for the middle section of SMA filament here is calculated according to the relation

$${\displaystyle \frac{d{E}_{20}}{dy}}=\left\{\hspace{0.17em}\begin{array}{c}\begin{array}{l}\left(E_m-E_a\right){\displaystyle \frac{d{R}_m}{dy}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0<{\epsilon }_{s2}\le {\epsilon }_m^y\\ \left(E_T-E_a\right){\displaystyle \frac{d{R}_m}{dy}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_m^y<{\epsilon }_{s2}\le {\epsilon }_m^d\\ \left(E_d-E_a\right)\hspace{0.17em}{\displaystyle \frac{d{R}_m}{dy}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_m^d<{\epsilon }_{s2}\end{array}\end{array}\right.,$$

(45)

where the derivative of the relative martensite fraction considering (21) is

$${\displaystyle \frac{d{R}_m}{dy}}=\left\{\begin{array}{c}{\displaystyle \frac{S_R{T}_h}{2{\lambda }_T{\theta }^2}}{e}^{-{\displaystyle \frac{1}{{\theta }^2}}\left(y-y_e\right)}\left(1-\tanh {\displaystyle \frac{T_m-T_h{e}^{-\frac{1}{{\theta }^2}\left(y-y_e\right)}+T_0}{{\lambda }_T}}\right)\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}y\le y_e\\ 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}{y}_e<y<y_f\\ {\displaystyle \frac{S_R{T}_h}{2{\lambda }_T{\theta }^2}}{e}^{-{\displaystyle \frac{1}{{\theta }^2}}\left(y-y_f\right)}\left(1-\tanh {\displaystyle \frac{T_m-T_h{e}^{-\frac{1}{{\theta }^2}\left(y-y_f\right)}+T_0}{{\lambda }_T}}\right)\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}y>y_f\end{array}\right.,$$

(46)

From the second equation of the system (42), it can be concluded that the magnetic system has a contribution which is expressed by a potential force of the form

$$F_{m\Sigma }={\displaystyle \frac{2\kappa }{{\left({\delta }_0+y_1+y\right)}^2}}$$

(47)

$$F_{m\Sigma }=-{\displaystyle \frac{2\kappa }{{\left({\delta }_0+y_1+y\right)}^2}}$$

(48)

This force presented in Maclaurin’s series leads to

$$F_{m\Sigma }\approx -{\displaystyle \frac{2\kappa }{{\left({\delta }_0+y_1\right)}^2}}+{\displaystyle \frac{4\kappa }{{\left({\delta }_0+y_1\right)}^3}}у-O\left({у}^2\right)$$

(49)

The factor in front of $у$ indicates that the magnetic system will increase the natural mechanical frequency of vibration along the $у$ coordinate.

The system of ordinary differential equations (42) was solved numerically using the data in Table 1.

4. Experimental Studies on the Influence of Magnetic Tuning

The experimental investigations of the influence of the magnetic system on the natural frequency tuning were carried out using a dedicated experimental setup. The schematic concept of this experimental setup is shown in Figure 5(a). Figure 5 (b) shows a top view of the whole experimental setup. Figure 5 (c) is a close-up view of the left side of the energy harvester with the PVDF cantilever 2, horizontal position sensor 3, ferromagnetic mass 4, heater 5 and magnets 10. Figure 5(d) is a picture of the right side of the energy harvester showing the sphere 4, heater 5, vertical position sensor 6, SMA filament 7, right support 8 and magnets 10.

In the experimental setup to the energy harvester comprising the heater 5 located in the middle strand of the CMA 7, the ferromagnetic spheres 4 on which the NdFeB magnets 10 act directly. Sensors are added to determine the horizontal 3 and vertical 6 displacement of the masses. These sensors are of the type of APO-075-002-000 from TT Elastronics Ltd. UK. Variable load resistors 1 are connected in parallel to the outputs of the two piezoelectric beams 2, the resistance of which is adjusted by an ohmmeter 9. The voltages from the piezoelectric beams and the two position sensors are fed into a Data Acquisition System DAQ USB - 6003 from National Instruments, then processed by a laptop using a LabVIEW 11 program to save to an Excel file.

In order to investigate the influence of magnets on the output voltage, two types of experiments were conducted. In the first, the results of the displacement along the two axes and the generated output voltage with the magnets removed were recorded. In the second type of experiments, the same signals were recorded but in the presence of magnets. In these experiments, the location and number of magnets was chosen to obtain the maximum output voltage.

A comparison between the output voltages of the energy harvester without and with permanent magnets is shown in Figure 6. The plot of the resulting voltage without magnets is shown in Figure 6(a), and the same output voltage but with magnets added is given in Figure 6(b). It can be seen that the magnet increases the vibration frequency of the piezoelectric cantilevers and furthermore increases the magnitude of the generated output voltage.

From the graphs, it is found that the magnets resulted in an increase in the output peak-to-peak voltage from 6.12 to 14,14 V, which is about 130% increase in output voltage of the energy harvester with magnets. The magnets increased the output voltage frequency from 8,3 Hz to 9,2 Hz, which in this case is about an 11% increase in frequency.

Maximum power P=1.9 µW for energy harvester without magnet is obtained at load resistance of R_L=7 MOhms. For energy harvester with magnets the obtained maximum power P=8.18 µW at a load resistance R_L = 5 MOhms. This result shows that the magnets contributed to increasing the output optimum power by 330%.

5. Comparison of Theoretical and Experimental Data to Validate the Dynamic Model

The dynamical model consisting of the system of differential equations (42) differs from the model solutions in [40] in the added potential magnetic energy. For the solution of the system the same approach based on implicit fifth order Euler method and the same initial conditions according to Table 1 was used.

Figure 8 (a), (b) and (c) show the experimental results for the output voltage, the displacement of the spheres along the y-axis and x-axis, respectively. For comparison, Figure 8 (d), (e) and (f) show the results obtained from the solution of system (42).

The numerical simulation results show a reasonably good confirmation of the experimental data. In analogous studies in [40], it was found that due to the multiple hystereses and the variable nature of the SMA filament heating, the vibrations have a pseudo-stochastic character that affects both the amplitudes and the frequencies and phases. The variation of the temperature modulus with the Young’s modulus leads to a change in the frequencies of three types of oscillations one thermal along the x-axis and two mechanicals along both axes x, and y. The x-axis displacements sum, which is the cause of the modulated vibrations along this axis.

The system of ordinary differential equations [42] is highly non-linear due to discontinuities on the right-hand sides of the functions caused by hysteresis. This makes them stiff and difficult to solve numerically. Attempts to use standard methods such as Runge-Kutta-Fehlberg 7(8) were unsuccessful, and a solution was only obtained by controlling the time step and number of iterations using the implicit Euler method [49] (see the Maple code form Supplementary Material).

6. Conclusions

In a newly developed thermal energy harvester with an SMA filament, a system of permanent magnets was added to influence the displacement of ferromagnetic masses attached to the SMA filaments. Since the heat source is at a constant temperature, thermal energy conversion is achieved by self-excited vibrations of the SMA filament and subsequent mechanoelectrical conversion by piezoelectric PVDF cantilevers. It has been shown experimentally and theoretically that permanent magnets increase the natural vibration frequency of the system SMA filament and ferromagnetic masses.

As a result of the action of the permanent magnets, the output voltage power of the thermal energy harvester increased from 1.9 µW to 8.18 µW, which is an increase of about 330%. An increase in the output peak-to-peak voltage from 6.12 V to 14.14 V was also observed.

The experimental setup developed is capable of measuring vibration processes in real time, where no integration is required to determine the displacement. This is because the signal for the displacement along both axes is obtained from a non-contact inductive sensor which only needs to be calibrated. If an inertial sensor such as an accelerometer were used, the double integration required to calculate the displacement would introduce large errors.

An interesting feature of the thermal energy harvester under consideration is that it can operate at cryogenic temperatures due to the wide temperature range over which the phase transformation of NiTi can be varied.