Investigative Study of the Effect of Damping and Stiffness Nonlinearities on an Electromagnetic Energy Harvester at Low Frequency Excitations

1. Introduction

Several studies have been conducted towards broadening the operational bandwidth of a Vibration Energy Harvester (VEH), beyond the resonant region [1,2,3,4,5,6,7]. To extend the operational bandwidth of a VEH system, the authors in [1] investigated a broadband energy harvester whose function is based on a combination of nonlinear stiffening effect and multimodal energy harvesting in order to attain a high bandwidth over a broad range of excitations (0.1–2.0g). In [2], the author extended some previous studies on using movable masses to extend the bandwidth of VEHs. The author demonstrated a novel method that involved embedding liquid in the system’s mass, used to extend the bandwidth of the VEH or tune the frequency without a significant reduction in the power output. Ramlan et al., in [3], demonstrated the potential benefits of nonlinear stiffness characteristics in energy harvesters. Two implementations of nonlinear stiffness characteristics were considered in the study. For the first implementation, using a bi-stable snap-through mechanism, it was shown that more electrical power was harvested, compared to a tuned linear device, for a given input excitation. Likewise, for the second implementation, using a hardening spring, it was also demonstrated that the bandwidth could also be extended, in comparison with an equivalent linear device. It should be noted that an equivalent linear VEH device provides the same maximum throw, at resonance, as the nonlinear VEH device of interest. Meanwhile, the maximum throw is dependent on the size of the VEH system (typically the electromagnetic-type) and is defined as the maximum distance that the mass of the VEH system can travel. This constitutes a mass-displacement constraint for electromagnetic-type energy harvesters. Wang et al., in [4], presented a novel automated method that tracks the resonant frequency of a VEH. The authors achieved this by employing a pair of cylindrical movable magnetic sliders on a cantilever beam, which increased the bandwidth of the VEH, as the slider could track the resonant frequency on the cantilever beam without manual involvement or external energy input. In [5], the authors suggested tuning the resonant frequency of a VEH system to align with that of the excitation frequency and set the electrical damping to be equal to the parasitic damping. The system was implemented using Power electronics circuits, which enabled the adjustment of both the damping characteristics and the resonant frequency, thus improving its efficiency. Studies in [6] and [7] focused on the comparison of the bandwidth provided by a Duffing-type energy harvester, with that of an equivalent linear harvester. The results obtained also demonstrated that a nonlinear harvester provided a larger bandwidth compared to the linear equivalent. Most of the published works have compared the Duffing-type VEH with its equivalent linear device. In addition, several parameter optimisations have been suggested to achieve similar results. Such parameter optimisations include the use of a mechanical damping [6], the integration of an optimum electrical load [7,8], and the application of multi-stage harvester models [9,10].

It has been evidently established in literature [11,12,13,14,15] that a nonlinear stiffness behaviour, as depicted in Figure 1, can be implemented using several geometric arrangements of permanent magnets, whose lines of flux cut across defined arrangements of copper coils. In [11], a nonlinear dual-function multi-modal energy harvester and vibration absorber (EHVA) for harvesting energy and suppressing vibration in low-medium frequency band, was presented. Two different methods were employed to extend the operating bandwidth of the system. These include the design of the multi-modal shapes of the EHVA as well as the hysteresis property of nonlinear softening springs, implemented by a novel permanent magnets structure. The authors in [12] investigated the steady state response of a specific VEH system under the condition of external and internal resonance, with focus on the double jump phenomenon. The frequency response curve shows the existence of resonance peaks tilting to the left and right of the natural frequency of the system. Wang and Zhu in [13] coupled a magnetic multi-stable device to a pendulum VEH in order expand its bandwidth, specifically in low-frequency operation. In [14], an experimental and theoretical study for designing a nonlinear electromagnetic converter-based magnetic spring was conducted. In this study, a special emphasis was given to the magnetic force acting on the moving magnet, based on two parameters – the volume of the magnets and the geometry of the two fixed magnets (i.e., disc or ring). Meanwhile, the authors in [15] numerically analysed a magnetic-spring-based electromagnetic energy harvester with piecewise nonlinear stiffness. It was demonstrated that the piecewise nonlinear stiffness behaviour, developed due to the interactions of the moving magnet’s flux on the coil, facilitated the response of the system in a wider frequency range, enabling generate more electrical energy.

Similarly, system damping characteristics have been employed in optimising the energy harvested by vibration energy harvesters [16,17,18,19,20,21,22,23,24]. In [16], the authors presented an unconventional way of achieving a tuneable resonant frequency (from high frequency to ultralow frequency) and extending both the bandwidth and peak of the energy harvested, simultaneously, by utilising distinct structurally supported displacement-dependent nonlinear damping property. This work was further extended by the authors in [17], where a scissor-like energy harvesting system, with equivalent nonlinear damping and linear stiffness characteristics, was developed. It was demonstrated that the scissor-structure provided beneficial nonlinear damping effects, thereby significantly improving the magnitude of the power harvested, as well as the bandwidth over which it was harvested. The beneficial effect of antisymmetric nonlinear damping to energy harvesters, was demonstrated by the authors in [18], when the system is subject to ambient random vibrations. Meanwhile, in [19], a VEH which employs a complementary metal-oxide-semiconductor-compatible 3D micro-electromechanical system coils and a ferromagnetic core, was presented. In order to describe the nonlinear electromagnetic damping coefficient and nonlinear attraction between the magnet and the ferromagnetic core, a systematic model was proposed. Thereafter, a vibration model was developed by considering nonlinear stiffness and damping coefficient to derive the dynamic characteristics and output performance of the system. The authors in [20] proposed a novel H-bridge circuit-based electromagnetic damper which enables a bi-directional flow of electrical energy between the electromagnetic damper and the energy storage. It was also demonstrated that this process enables the realisation of diverse damping behaviour with salient self-powered feature.

The authors in [21,22] also employed nonlinear damping in extending the energy harvested by a vibration energy harvester. Based on the findings in [21] and [22], an analysis, design, and optimisation of a nonlinear VEH system was conducted in [23] and [24]. While no mass-displacement constraint was considered in [23], this was considered in [24]. In these studies, an optimum cubic damping parameter was designed for a desired power level, using the Output Frequency Response Function (OFRF) method. The OFRF of a nonlinear system is determined based on the class of nonlinear differential equation (NDE) the system belongs to and it shows the relationship between the output spectrum of the system and its nonlinear parameters [25,26]. It, thus, describes the system characteristics. The OFRF representation of the system studied in [25,26,27,28,29], were determined using the Least-Squares (LS) approach. However, this method requires several numerical simulations, using a (training) set of values of the system design parameters, to obtain the respective system output responses [30].

Vazquez et al. in [31], based on the characteristics of the nth-order Volterra operator being a multi-linear function of a combination of input signals, modelled the behaviour of the Volterra operators using their Associated Linear Equation (ALE) decompositions. These ALE decompositions, as discussed in [31,32], can be used as an analytical tool for analysing the Volterra class of nonlinear systems such as the Duffing equation. Based on this, it was further revealed in [33,34] that the ALE decompositions, for a Volterra class of nonlinear systems, can be used to determine a more accurate OFRF representation of the system, using a significantly lesser number of numerical simulations, compared to the LS approach. These methods were further extended in [35] to the Generalised Associated Linear Equation (GALE) decompositions, which considered a general class of nonlinear damping.

In this study, an analysis and design of a nonlinear VEH system is conducted, using the OFRF representations of the system output spectra, which are determined from the GALE decompositions of the nonlinear VEH model. In addition to using a nonlinear damping component to extend the average power of the VEH, a stiffness nonlinearity is also integrated to widen the operational frequency range of the harvesting device. It should be noted that the current study is an extension of the initial work by the authors in [34]. To the best of the authors’ knowledge, this is the first time the OFRF method, derived using the GALE decompositions, is employed in the design and optimisation of a nonlinear vibration energy harvester. Using the OFRF model, derived from the GALE decompositions, simplifies the design process since a polynomial of the system’s performance metrics, is derived in terms of the design parameters. A sequel study involving an experiemental validation of the design, is currently ongoing.

Subsequent sections of this paper are organised as follows—Section 2 presents the model formulation of the system of interest. In Section 3, the OFRF method is introduced, while Section 4 describes the determination of the OFRF structure. Section 5 discusses the evaluation of the GALE decompositions and in Section 6, the determination of the OFRF model, using the decomposed GALE contributions is demonstrated. Section 7 provides the results obtained and their corresponding discussions while the research findings are concluded in Section 7.

2. Model Formulation

A single degree-of-freedom (SDOF) vibration-based energy harvester, as illustrated in Error! Reference source not found. 2, having a suspended mass, $m$ and an oscillating support-base with displacement $y(t)$, is given. The mass is separated from the base using an isolation system modelled as a nonlinear damping system, connected parallel to a nonlinear spring. The damping system comprises a mechanical viscous damping, $c_1$ and an electrical damping, $c_3$. The electrical damping arises from the electromagnetic force generated by virtue of the non-Ohmic load resistance connected across the EM damper. The linear and cubic stiffness coefficients are $k_1$ and $k_3$ respectively.

Figure 2. SDOF of a vibration energy harvester with nonlinear stiffness and damping.

Figure 2. SDOF of a vibration energy harvester with nonlinear stiffness and damping.

The model of the SDOF VEH is a class of NDE and the equation of motion of the mass with respect to the relative displacement $z=x-y$ is given as

$$m\ddot{z}+c_1\dot{z}+c_3{\dot{z}}^3+k_{1}z+k_3{z}^3=-m\ddot{y}$$

(1)

For a harmonic base displacement with amplitude,$Y$, frequency, $\omega$and zero phase shift, the base displacement is given as $y=Y\sin (\omega t)$. Therefore Equation (1) becomes

$$m\ddot{z}+c_1\dot{z}+c_3{\dot{z}}^3+k_{1}z+k_3{z}^3=m{\omega }^{2}Y\sin (\omega t)$$

(2)

The nonlinear damping device absorbs an instantaneous power equal to the product of the instantaneous damper force and relative velocity of the VEH. Therefore, it yields an average power given as

$$P_{av}=\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}\left(c_3{\dot{z}}^3\right)}\cdot \dot{z}dt$$

(3)

For a single-frequency harmonic oscillation, where $z=Z\sin (\omega t)$, yields

$$P_{av}=\frac{3}{8}{c}_3{\omega }^4{Z}^4$$

(4)

In addition, it can be deduced that the output frequency response $Z$, of system (2), is a function of $\omega$, as well as the nonlinear parameters,$c_3$ and $k_3$. This implies that $P_{av}$, derived in Equation (4), is also a function of $c_3$, $k_3$ and $\omega$. Note that the resonant frequency is the frequency of interest, as the maximum power is extracted at this frequency.

3. The Output Frequency Response Function (OFRF)

Equation (2) describes the dynamic model of the VEH system, which is a typical base-excited duffing-equation, with an integrated nonlinear damping. A similar system has been studied using methods such as, the Harmonic balance method [36], and multiple scales [37]. However, the aforementioned methods only facilitate the analytical study of nonlinear systems. On the other hand, the OFRF, which is employed in this study, does not only facilitate the analytical study of nonlinear systems but also enables the design and optimization of such systems. However, for the OFRF method to perform appropriately, the system of interest must operate within a stable regime. The main benefit of the OFRF method is that it can provide an explicit analytical relationship between the design objective and the system nonlinear parameters. This can significantly facilitate the system’s design and optimization process. A comprehensive explanation of the OFRF concept can be found in [25,26,27,28,29,30].

Let us examine the differential equation in Equation (5), which describes a class of Volterra systems

$${\displaystyle \sum _{\overline{m}=1}^{\overline{M}}{\displaystyle \sum _{p=0}^{\overline{m}}{\displaystyle \sum _{{\mathrm{k}}_1,\mathrm{...}{\mathrm{k}}_{\overline{m}}=0}^{\mathrm{K}}{c}_{p,\overline{m}-p}}}}({\mathrm{k}}_1,\cdots {\mathrm{k}}_{\overline{m}}){\displaystyle \prod _{i=1}^p\frac{d^{{\mathrm{k}}_i}z(t)}{d{t}^{{\mathrm{k}}_i}}}{\displaystyle \prod _{i=p+1}^{\overline{m}}\frac{d^{{\mathrm{k}}_i}y(t)}{d{t}^{{\mathrm{k}}_i}}}=0$$

(5)

where $\overline{M}$ is the maximum degree of nonlinearity, in terms of the system’s input, $y(t)$ and output, $z(t)$, and $\mathrm{K}$ is the order of the derivative. According to the OFRF method, as described in [26], the output frequency response of Equation (5) can be described by a polynomial function in terms of the system nonlinear parameters as

$$Z(\mathrm{j}\omega )={\displaystyle \sum _{{\delta }_1=0}^{n_1}\cdots {\displaystyle \sum _{{\delta }_{S_N}=0}^{n_{S_N}}{\Psi }_{({\delta }_1,\dots ,{\delta }_{S_N})}}}(\mathrm{j}\omega ){\kappa }_1^{{\delta }_1}\dots {\kappa }_{S_N}^{{\delta }_{S_N}}$$

(6)

where $n_i$ are the maximum order of ${\kappa }_i$, for $i=1,\dots ,S_N$ in the polynomial expression of the output spectrum, $Z(\mathrm{j}\omega )$ of Equation (6). The OFRF coefficients, ${\Psi }_{({\delta }_1,\dots ,{\delta }_{S_N})}(\mathrm{j}\omega )$ are frequency functions with complex values. They are also dependent on the system linear parameters and system input, where ${\delta }_i=0,\dots ,n_i$ and $i=1,\dots ,S_N$. In addition, ${\kappa }_1^{{\delta }_1}\dots {\kappa }_{S_N}^{{\delta }_{S_N}}$ is known as the OFRF structure. They are a set of monomials in terms of the system nonlinear characteristic parameters. If the set of monomials in the OFRF polynomial, of the nth-order output spectrum, is denoted as $\mathfrak{M}$ and the vector of the frequency function, is denoted as $\Theta (\mathrm{j}\omega )$, the OFRF can be then be described as

$$Z(\mathrm{j}\omega )=\mathfrak{M}\cdot \Theta {(\mathrm{j}\omega )}^{\mathrm{T}}$$

(7)

Where

$$\mathfrak{M}={\displaystyle \underset{n=1}{\overset{s_N}{\cup }}{\mathrm{{\rm E}}}_n}$$

(8)

Here, $S_N$ is the maximum order of nonlinearity considered for this study and the set of monomials ${\mathrm{{\rm E}}}_n$ can be derived using the method in [15] as

$$\begin{array}{l}{\mathrm{{\rm E}}}_n=\left[{\displaystyle \underset{{\mathrm{k}}_1,\dots ,{\mathrm{k}}_{\mathrm{n}}=0}{\overset{\mathrm{K}}{\cup }}\left[c_{0,n}({\mathrm{k}}_1,\dots ,{\mathrm{k}}_{\mathrm{n}})\right]}\right]\cup \left[{\displaystyle \underset{\overline{m}-p=1}{\overset{\mathrm{n}-1}{\cup }}{\displaystyle \underset{p=1}{\overset{\mathrm{n}-(\overline{m}-p)}{\cup }}{\displaystyle \underset{{\mathrm{k}}_1,\dots ,{\mathrm{k}}_{\mathrm{n}}=0}{\overset{\mathrm{K}}{\cup }}\left(\left[c_{p,(\overline{m}-p)}({\mathrm{k}}_1,\dots ,{\mathrm{k}}_{\overline{m}})\right]\otimes {\mathrm{{\rm E}}}_n{}_{-(\overline{m}-p),p}\right)}}}\right]\\ \cup \left[{\displaystyle \underset{p=2}{\overset{\mathrm{n}}{\cup }}{\displaystyle \underset{{\mathrm{k}}_1,\dots ,{\mathrm{k}}_{\mathrm{n}}=0}{\overset{\mathrm{K}}{\cup }}\left(\left[c_{p,0}({\mathrm{k}}_1,\dots ,{\mathrm{k}}_{\overline{m}})\right]\otimes {\mathrm{{\rm E}}}_{n,p}\right)}}\right]\end{array}$$

(9)

Note that the character ‘$\otimes$ ’ is the Kronecker product and given by

$${\mathrm{{\rm E}}}_{n,p}={\displaystyle \underset{i=1}{\overset{\mathrm{n}-p+1}{\cup }}\left({\mathrm{{\rm E}}}_i\otimes {\mathrm{{\rm E}}}_{\mathrm{n}-i,p-1}\right),}{\mathrm{{\rm E}}}_{\mathrm{n},1}={\mathrm{{\rm E}}}_{\mathrm{n}},{\mathrm{{\rm E}}}_1=[1]$$

(10)

Then the set of monomials can be obtained as $\mathfrak{M}={\displaystyle \underset{n=1}{\overset{s_N}{\cup }}{\mathrm{{\rm E}}}_n}$

4. Determination of the OFRF Structure

Firstly, the OFRF representations of the output spectrum of system (2), $Z(\mathrm{j}\omega )$ and the harvestable average power of Equation (4), $P_{av}$, are obtained in terms of the design parameters, $c_3$ and $k_3$. It is observed that Equation (2) belongs to the class of Volterra system of Equation (4) in [25], with $\overline{M}=3$ and $\mathrm{K}=2$. The system parameters are obtained as $c_{10}(2)=m$, $c_{10}(1)=c_1$, $c_{10}(0)=k_1$, $c_{30}(000)=k_3$, $c_{30}(111)=c_3$, and $c_{01}(0)=-m{\omega }^{2}Y$.

If the set of monomials in the OFRF representation of the nth-order output spectrum of system (2), is denoted by $\mathfrak{M}$, and the complex-valued OFRF coefficients is denoted by $\Theta (\mathrm{j}\omega )$, then the OFRF representation can be written as

$$Z(\mathrm{j}\omega )=\mathfrak{M}\cdot \Theta (\mathrm{j}\omega )$$

(11)

Applying the algorithm, as presented in [38], to obtain the OFRF monomials, $\mathfrak{M}$, up to the 7th order, yields

$$\begin{array}{l}{\mathrm{{\rm E}}}_1=[1]\\ {\mathrm{{\rm E}}}_3=[c_3{k}_3]\\ {\mathrm{{\rm E}}}_5=[c_3^2{c}_3{k}_3{k}_3^2]\\ {\mathrm{{\rm E}}}_7=[c_3^3{c}_3^2{k}_3{c}_3{k}_3^2{k}_3^3]\end{array}$$

(12)

Therefore, $\mathfrak{M}={\displaystyle \underset{n=1}{\overset{s_N}{\cup }}{\mathrm{{\rm E}}}_n}=[1,c_3,k_3,c_3^2,c_3{k}_3,k_3^2,c_3^3,c_3^2{k}_3,c_3{k}_3^2,k_3^3]$

It should be noted that for improved accuracy, higher orders can be considered. Furthermore, the OFRF representation, as presented in Equation (11), which comprises the monomials obtained, as presented in Equation (12), and its respective OFRF coefficients, ${\Theta }_{n|r}(\mathrm{j}\omega )$ (yet to be determined), can be represented in the form

$$Z_{\mathrm{OFRF}}(\mathrm{j}\omega )={\displaystyle \sum _{n=1}^{S_N}{\mathrm{{\rm E}}}_{n|r}\cdot {\Theta }_{n|r}(\mathrm{j}\omega )}$$

(13)

where $r=0,2,\dots ,h$ and $h$ is the maximum number of elements in ${\mathrm{{\rm E}}}_n$. Rewriting Equation (13) yields

$$\begin{array}{l}{Z}_{\mathrm{OFRF}}(\mathrm{j}\omega )={\Theta }_{1|0}(\mathrm{j}\omega )+c_3\cdot {\Theta }_{3|1}(\mathrm{j}\omega )+k_3\cdot {\Theta }_{3|2}(\mathrm{j}\omega )\\ +c_3^2\cdot {\Theta }_{5|1}(\mathrm{j}\omega )+c_3{k}_3\cdot {\Theta }_{5|2}(\mathrm{j}\omega )+k_3^2\cdot {\Theta }_{5|3}(\mathrm{j}\omega )+c_3^3\cdot {\Theta }_{7|1}(\mathrm{j}\omega )\\ +c_3^2{k}_3\cdot {\Theta }_{7|2}(\mathrm{j}\omega )+c_3{k}_3^2\cdot {\Theta }_{7|3}(\mathrm{j}\omega )+k_3^3\cdot {\Theta }_{7|4}(\mathrm{j}\omega )\end{array}$$

(14)

Thus, to determine the OFRF coefficients, ${\Theta }_{n|r}(\mathrm{j}\omega )$, the GALE decompositions of system (2) are first computed up to the 7th-order. In the next section, the evaluation of the GALE decompositions, for the NDE system of interest, is demonstrated.

5. Evaluation of the Generalised Associated Linear Equations (GALEs)

For a nonlinear system of the Volterra class given in system (2), the following substitutions can be made

$$z(t)={\displaystyle \sum _{n=1}^{\infty }{z}_n(t)}$$

(15)

Rewriting system (2) in a general form, by leaving all the linear elements on the LHS and substituting Equation (15) yields

$$\begin{array}{l}m{\displaystyle \sum _{n=1}^{\infty }{\ddot{z}}_n}+c_1{\displaystyle \sum _{n=1}^{\infty }{\dot{z}}_n}+k_1{\displaystyle \sum _{n=1}^{\infty }{z}_n}=\\ m{\omega }^{2}Y\sin (\omega t)-{\displaystyle \sum _{\underset{\_}{j}=3}^{\underset{\_}{L}}{c}_{\underset{\_}{j}}{\left({\displaystyle \sum _{n=1}^{\infty }{\dot{z}}_n}\right)}^{\underset{\_}{j}}}-{\displaystyle \sum _{\underset{˜}{j}=3}^{\underset{˜}{L}}{k}_{\underset{˜}{j}}{\left({\displaystyle \sum _{n=1}^{\infty }{z}_n}\right)}^{\underset{˜}{j}}}\end{array}$$

(16)

The GALE decompositions of system (16) can be obtained for the nth order, for $n=1,2,\dots ,S_N$, where $S_N$ is the maximum order of the system nonlinearity considered, as demonstrated in [34,35] thus,

$$\begin{array}{l}m{\displaystyle \sum _{n=1}^{S_N}{\ddot{z}}_n}+c_1{\displaystyle \sum _{n=1}^{S_N}{\dot{z}}_n}+k_1{\displaystyle \sum _{n=1}^{S_N}{z}_n}=\\ m{\omega }^{2}Y\sin (\omega t)\\ -{\displaystyle \sum _{n=1}^{S_N}{\displaystyle \sum _{\underset{\_}{j}=3}^{n=3}{c}_{\underset{\_}{j}}{\displaystyle \sum _{{\underset{\_}{j}}_1=1}^{n-\underset{\_}{j}+1}\cdots }}}{\displaystyle \sum _{{\underset{\_}{j}}_i=1}^{n-\underset{\_}{l}+i-{\underset{\_}{j}}_1-\cdots {\underset{\_}{j}}_{i-1}}\cdots }{\displaystyle \sum _{{\underset{\_}{j}}_l=0}^{n-{\underset{\_}{j}}_1-\cdots {\underset{\_}{j}}_i\cdots -{\underset{\_}{j}}_{l-1}}{\dot{z}}_{{\underset{\_}{j}}_1}{\dot{z}}_{{\underset{\_}{j}}_i}\dots {\dot{z}}_{{\underset{\_}{j}}_l}}\\ -{\displaystyle \sum _{n=1}^{S_N}{\displaystyle \sum _{\underset{˜}{j}=3}^{n=3}{k}_{\underset{˜}{j}}{\displaystyle \sum _{{\underset{˜}{j}}_1=1}^{n-\underset{˜}{j}+1}\cdots }}}{\displaystyle \sum _{{\underset{˜}{j}}_i=1}^{n-\underset{˜}{j}+i-{\underset{˜}{j}}_1-\cdots {\underset{˜}{j}}_{i-1}}\cdots }{\displaystyle \sum _{{\underset{˜}{j}}_l=0}^{n-{\underset{˜}{j}}_1-\cdots {\underset{˜}{j}}_i\cdots -{\underset{˜}{j}}_{i-1}}{z}_{{\underset{˜}{j}}_1}{z}_{{\underset{˜}{j}}_i}\dots z_{{\underset{˜}{j}}_l}}\end{array}$$

(17)

where the summation of all the sub-indices of ${\dot{z}}_{\underset{\_}{j}}$ and $z_{\underset{˜}{j}}$on the RHS has to be equal to $n$i.e.,$({\underset{\_}{j}}_1+\dots +{\underset{\_}{j}}_l=n)$ and $({\underset{˜}{j}}_1+\dots +{\underset{˜}{j}}_l=n)$. In computing the GALE decompositions, the low-order output responses contribute to the immediate higher-order responses up to the maximum order considered. For an estimation of the total output responses up to the $S_N\mathrm{th}-\mathrm{order}$ and its corresponding output spectrums,

$$\left\{\begin{array}{l}z(t)={\displaystyle \sum _{n=1}^{S_N}{z}_n(t)}\\ Z(\mathrm{j}\omega )={\displaystyle \sum _{n=1}^{S_N}{Z}_n(\mathrm{j}\omega )}\end{array}\right.$$

(18)

For $S_N=7$, the following GALE decompositions are obtained

$$\left\{\begin{array}{l}m{\ddot{z}}_1+c{\dot{z}}_1+k{z}_1=m{\omega }^{2}Y\sin (\omega t)\\ m{\ddot{z}}_3+c{\dot{z}}_3+k{z}_3=-k_3{z}_1^3-c_3{\dot{z}}_1^3\\ m{\ddot{z}}_5+c{\dot{z}}_5+k{z}_5=-3k_3{z}_1^2{z}_3-3c_3{\dot{z}}_1^2{\dot{z}}_3\\ m{\ddot{z}}_7+c{\dot{z}}_7+k{z}_7=-3k_3(z_1{z}_3^2+z_1^2{z}_5)-3c_3({\dot{z}}_1{\dot{z}}_3^2+{\dot{z}}_1^2{\dot{z}}_5)\end{array}\right.$$

(19)

The continuous-time output response of system (2) and its corresponding output spectrum, where $Z_n(\mathrm{j}\omega )=\mathrm{fft}(z_n(t))$ , are respectively expressed as

$$\left\{\begin{array}{l}z(t)=z_1(t)+z_3(t)+z_5(t)+z_7(t)\\ Z(\mathrm{j}\omega )=Z_1(\mathrm{j}\omega )+Z_3(\mathrm{j}\omega )+Z_5(\mathrm{j}\omega )+Z_7(\mathrm{j}\omega )\end{array}\right.$$

(20)

The cumulative structure of the individual nth-order GALE contributions, up to the 7th -order, is presented in Figure 3. Similarly, Figure 4 shows the output spectrum for each nth-order contribution of the GALE decompositions, up to the 7th-order. It is observed that at resonance there is a significant contribution by the individual decompositions. Meanwhile, Figure 5 demonstrates the nth-order contributions of the GALE decompositions, up to the 9th-order for a range of nonlinear damping values.

6. Determination of the OFRF Model using the GALE Contributions

Equation (20) shows the output spectrum of system (2) determined from the Fast Fourier Transform of the GALE contributions obtained. The nth-order output spectrum of each GALE contribution is equal to the corresponding nth-order component of the OFRF representation thus

$$\left\{\begin{array}{l}{Z}_{\mathrm{ALEs}}(\mathrm{j}\omega )=Z_{\mathrm{OFRF}}(\mathrm{j}\omega )\\ {\displaystyle \sum _{n=1}^{S_N}{Z}_n(\mathrm{j}\omega )}={\displaystyle \sum _{n=1}^{S_N}{\mathrm{{\rm E}}}_{n|r}\cdot {\Theta }_{n|r}(\mathrm{j}\omega )}\end{array}\right.$$

(21)

From Equation (21), it can be deduced that

$$\begin{array}{l}{Z}_1(\mathrm{j}\omega )={\Theta }_{1|0}(\mathrm{j}\omega )\\ Z_3(\mathrm{j}\omega )=c_3\cdot {\Theta }_{3|1}(\mathrm{j}\omega )+k_3\cdot {\Theta }_{3|2}(\mathrm{j}\omega )\\ Z_5(\mathrm{j}\omega )=c_3^2\cdot {\Theta }_{5|1}(\mathrm{j}\omega )+c_3{k}_3\cdot {\Theta }_{5|2}(\mathrm{j}\omega )+k_3^2\cdot {\Theta }_{5|3}(\mathrm{j}\omega )\\ Z_7(\mathrm{j}\omega )=c_3^3\cdot {\Theta }_{7|1}(\mathrm{j}\omega )+c_3^2{k}_3\cdot {\Theta }_{7|2}(\mathrm{j}\omega )+c_3{k}_3^2\cdot {\Theta }_{7|3}(\mathrm{j}\omega )\\ +k_3^3\cdot {\Theta }_{7|4}(\mathrm{j}\omega )\end{array}$$

(22)

Subsequent analysis in this study has been done using the following system parameter values presented in Table 1 with $\Omega =\omega /{\omega }_n$.

To obtain the OFRF coefficients up to the 7th-order, five simulations are required using five different values of $c_3(c_{3r})$and $k_3(k_{3r})$, where $r=1,2,3,4,5.$, as given in Table 2.

The OFRF coefficients can be determined for any frequency of interest using Equation (23) given as

$$\begin{array}{l}{\Theta }_{1|0}(\mathrm{j}\omega )=Z_1(\mathrm{j}\omega )\\ \left[\begin{array}{c}{\Theta }_{3|1}(\mathrm{j}\omega )\\ {\Theta }_{3|2}(\mathrm{j}\omega )\end{array}\right]={\left[\begin{array}{cc}{c}_{31}& k_{31}\\ c_{32}& k_{32}\end{array}\right]}^{-1}\left[\begin{array}{c}{Z}_{31}(\mathrm{j}\omega )\\ Z_{32}(\mathrm{j}\omega )\end{array}\right]\\ \left[\begin{array}{c}{\Theta }_{5|1}(\mathrm{j}\omega )\\ {\Theta }_{5|2}(\mathrm{j}\omega )\\ {\Theta }_{5|3}(\mathrm{j}\omega )\end{array}\right]={\left[\begin{array}{ccc}{c}_{31}^2& c_{31}{k}_{31}& k_{31}^2\\ c_{32}^2& c_{32}{k}_{32}& k_{32}^2\\ c_{33}^2& c_{33}{k}_{33}& k_{33}^2\end{array}\right]}^{-1}\left[\begin{array}{c}{Z}_{51}(\mathrm{j}\omega )\\ Z_{52}(\mathrm{j}\omega )\\ Z_{53}(\mathrm{j}\omega )\end{array}\right]\\ \left[\begin{array}{c}{\Theta }_{7|1}(\mathrm{j}\omega )\\ {\Theta }_{7|2}(\mathrm{j}\omega )\\ {\Theta }_{7|3}(\mathrm{j}\omega )\\ {\Theta }_{7|4}(\mathrm{j}\omega )\end{array}\right]={\left[\begin{array}{cccc}{c}_{31}^3& c_{31}^2{k}_{31}& c_{31}{k}_{31}^2& k_{31}^3\\ c_{32}^3& c_{32}^2{k}_{32}& c_{32}{k}_{32}^2& k_{32}^3\\ c_{33}^3& c_{33}^2{k}_{33}& c_{33}{k}_{33}^2& k_{33}^3\\ c_{34}^3& c_{34}^2{k}_{34}& c_{34}{k}_{34}^2& k_{34}^3\end{array}\right]}^{-1}\left[\begin{array}{c}{Z}_{71}(\mathrm{j}\omega )\\ Z_{72}(\mathrm{j}\omega )\\ Z_{73}(\mathrm{j}\omega )\\ Z_{74}(\mathrm{j}\omega )\end{array}\right]\end{array}$$

(23)

Therefore, the GALE-generated OFRF representation of the output spectrum of system (2) can be expressed as

$$\begin{array}{l}Z(\mathrm{j}\omega ;c_3,k_3)={\Theta }_{1|0}(\mathrm{j}\omega )+c_3\cdot {\Theta }_{3|1}(\mathrm{j}\omega )+k_3\cdot {\Theta }_{3|2}(\mathrm{j}\omega )+c_3^2\cdot {\Theta }_{5|1}(\mathrm{j}\omega )\\ +c_3{k}_3\cdot {\Theta }_{5|2}(\mathrm{j}\omega )+k_3^2\cdot {\Theta }_{5|3}(\mathrm{j}\omega )+c_3^3\cdot {\Theta }_{7|1}(\mathrm{j}\omega )+c_3^2{k}_3\cdot {\Theta }_{7|2}(\mathrm{j}\omega )\\ +c_3{k}_3^2\cdot {\Theta }_{7|3}(\mathrm{j}\omega )+k_3^3\cdot {\Theta }_{7|4}(\mathrm{j}\omega )\end{array}$$

(24)

The OFRF coefficients of Equation (24) have been determined using the GALE approach. The benefit of using the GALE approach is that the number of numerical simulations required to determine the OFRF of the system is significantly reduced [33]. To obtain the respective OFRF representation of the average power harvestable by the VEH system, via the nonlinear damping system, the OFRF representation of the output spectrum in Equation (24) is substituted in Equation (4) to yield

$$P_{av}(\omega ,c_3,k_3)=\frac{3}{8}{c}_3{\omega }^4{\left|Z(\mathrm{j}\omega ;c_3,k_3)\right|}^4$$

(25)

The OFRF representation of the quartic magnitude of the output spectrum, ${\left|Z(\mathrm{j}\omega ;c_3,k_3)\right|}^4$ can also be derived and represented in a polynomial form [20] as

$${\left|Z(\mathrm{j}\omega ;c_3,k_3)\right|}^4={\displaystyle \sum _{\overline{n}=0}^7{\displaystyle \sum _{\overline{m}=0}^{\overline{n}}{\upsilon }_{\overline{m},\overline{n}-\overline{m}}(\omega )c_3^{\overline{m}}{k}_3^{\overline{n}-\overline{m}}}}$$

(26)

where $\underset{\_}{N}=7$ is the maximum nth order nonlinearity while ${\upsilon }_{\overline{m},\overline{n}-\overline{m}}$ are functions of frequency and represent the OFRF coefficients of ${\left|Z(\mathrm{j}\omega ;c_3,k_3)\right|}^4$. Therefore, the average power can be represented as

$$P_{av}(\omega ,c_3,k_3)=\frac{3}{8}{\omega }^4\cdot {\displaystyle \sum _{\overline{n}=0}^7{\displaystyle \sum _{\overline{m}=0}^{\overline{n}}{\upsilon }_{\overline{m},\overline{n}-\overline{m}}(\omega )c_3^{\overline{m}+1}{k}_3^{\overline{n}-\overline{m}}}}$$

(27)

The OFRF-based results are obtained and compared with that determined using the Runge-Kutta 4 algorithm (ODE45 in MATLAB) for both the output spectrum of system (2) and the average power harvested by the VEH system. The comparisons were conducted for different combinations of parameter values beyond the training set of the nonlinear parameters $c_3$ and $k_3$. The results are presented in Figure 6 and Figure 7 for the pair of parameters,$c_3$= 0.45 Ns³m^-3, $k_3$= 250 Nm^-3and $c_3$= 0.25 Ns³m^-3, $k_3$= 270Nm^-3, respectively. In the next section, results obtained from the GALE-generated OFRF representations, presented in Equations (24) and (27), are provided and their implications discussed extensively.

7. Results and Discussion

The OFRF representation of the output spectrum, $Z$ of system (2) was derived using the GALE decompositions evaluated from the same system. The OFRF representation of the output spectrum was subsequently used to estimate the average power generated by the electromagnetic damper. This was performed for a pair of nonlinear parameter values, $c_3$ and $k_3$ beyond the range over which the OFRF representation was determined i.e., $c_3\in [0.3,0.4]$Ns³m^-3 and $k_3\in [0,220]$Nm^-3. As observed in Figure 6 and Figure 7, the OFRF representation accurately represents the actual output spectrum and average power of the VEH, respectively. These results clearly demonstrate a good match between the OFRF representation and the more accurate results from direct numerical simulations. This demonstrates the benefits of the OFRF methodology as it evidently describes the system dynamics over the entire spectrum. Note that the wobbles observed around the resonant regions in Figure 6 and Figure 7 are due to the use of parameters beyond the design (training) range. Using parameters further beyond the design range will cause the system to approach instability.

In the implementation of such a nonlinear VEH system, the cubic damping nonlinearity can be contributed by an electromagnetic damper whose characteristics is dependent on that of the load resistance of the energy harvesting circuit. The hardening stiffness nonlinearity can be realised by the application of magnetic springs, wherein a mass of permanent magnet is levitated between two stationary magnets [39]. However, such magnetic springs can also contribute a damping component typically referred to as magnetic damping [40].

7.1. Effect of a hardening spring and cubic damping on the VEH system

Nonlinearities contributed by a hardening-type spring stiffness is integrated into a standard linear harvesting device, to expand the bandwidth over which power can be harvested. The bandwidth expansion is caused by a shift in the resonant frequency to higher frequencies. To demonstrate the effect of integrating a hardening-type stiffness characteristic to the dynamics of a VEH system with cubic damping nonlinearity, numerical studies are conducted. The integration of a hardening spring-effect can be implemented by employing magnetic springs [39]. Using the OFRF representations provided in Equations (24) and (25), the effect of the stiffness parameter, $k_3$can be observed in Figure 8. It is clearly seen to extend the operational bandwidth of the nonlinear VEH system due to the shift in the resonant frequency of the output spectra. The operational bandwidth of the VEH system, with and without stiffness characteristic, are denoted as $△{\Omega }_2$ and $△{\Omega }_1$ respectively. In addition to this, an apparent increase in the dynamic range of the nonlinear VEH system can also be observed. This is logical as the average power of the nonlinear VEH system, given in Equation (25), is a function of the excitation frequency, $\omega$. Several studies in literature focused on increasing the bandwidth of linear VEH devices with the integration of a hardening-type stiffness and compared the duffing-type harvester with a standard linear harvester. Figure 9 shows the effect of varying the nonlinear stiffness characteristics on the output spectrum and average power of the nonlinear VEH system. In Figure 9a, while the relative displacement of the VEH system remains relatively constant as $k_3$ increases, it is evident in Figure 9b that the average power of the VEH device increases, within the resonant region, $\Omega =1.2$, as $k_3$ increases. However, the power remains relatively constant within the low and high frequency regions, where $\Omega \ll 1$ and $\Omega \gg 1$ respectively.

Furthermore, the effect of a variation of the nonlinear cubic damping $c_3$, on the VEH system, was also investigated. This was demonstrated by varying $c_3$ while fixing $k_3$, as shown in Figure 10. It is revealed that while the maximum span of the VEH system reduces by an insignificant amount, the average power harvested increases significantly around and beyond the resonant region. It should be noted that to implement the nonlinear cubic damping force characteristics, the current flowing through the nonlinear load (Energy harvesting circuit) should be proportional to the cube of the voltage across it [21].

In the current study, a hardening-type stiffness is integrated in a VEH system with cubic damping nonlinearity. However, a comparison of the VEH with damping and stiffness nonlinearities, against its linear equivalent, has not been considered here. This is due to the unavailability of a basis for such comparison to be made. Comparisons of this nature have never been reported for Electromagnetic-type VEH devices, with damping and stiffness nonlinearities, to the best of the authors’ knowledge.

Furthermore, most studies in literature considered the Duffing-type harvesters that exhibit the jump phenomenon, and they were majorly designed to operate within the larger stable branch. Nevertheless, it is imperative to note that if the VEH model experiences a jump phenomenon, the sum of the GALE decompositions will not converge to the actual output spectrum around the jump region. Therefore, the OFRF representations will poorly describe the actual output spectra of the system and, consequently, will be an inappropriate method to conduct the system analysis and design. However, this problem is resolved here with the integration of both linear and nonlinear damping characteristics.

7.2. Optimisation of an unconstrained nonlinear VEH system

Using the OFRF representation of the average power as determined in Equation (27), an optimisation problem can be formulated as

$$\begin{array}{l}\underset{\left[c_3,k_3\right]}{\max }{P}_{av}({\omega }_r,c_3,k_3)\\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{l}0.2\le c_3\le 0.5\\ 0\le k_3\le 400\end{array}\right.\end{array}$$

(28)

The solution to the unconstrained optimisation problem is simple and can be determined using the MATLAB fminsearch or fmincon function. Moreover, using the OFRF representations of Equations (24) and (27), the relationships between the design parameters,$c_3$, $k_3$, the output spectrum, and the average power, at the resonant frequency, can be established. These are presented using surface plots in Figure 11 and Figure 12. It can be deduced from Figure 12 that the average power is significantly sensitive to the nonlinear stiffness characteristic, $k_3$ but less sensitive to the nonlinear damping characteristic, $c_3$, at lower values of $k_3$. However, at higher values of $k_3$, the average power increases to an optimal value, as $c_3$ increases to 0.46, then it starts to decline. It should be noted that this design is only valid within the design range of the parameters of interest.

In Figure 13, the output spectrum and average power of the VEH system, for a variation of $c_3$ and $k_3$, are provided at $\Omega \ll 1$ and $\Omega \gg 1$.

8. Conclusions

In this study, the analysis and design of a vibration energy harvester with damping and stiffness nonlinearities were considered. Nonlinear stiffness was introduced into the mechanical subsystem of the vibration energy harvester, to improve the operational frequency range (bandwidth) of the VEH system. However, the nonlinear damping was introduced to extend the average power of the VEH system. While the nonlinear stiffness can be realised using magnetic springs or geometrically horizontal springs, the nonlinear damping can be realised using an energy harvesting circuit with nonlinear characteristics. Such an energy harvesting circuit will provide a current which varies proportionally with the cubic power of the voltage.

A polynomial representation of the system model, i.e., the Output Frequency Response Function (OFRF), was derived. However, the OFRF representation was determined using the Generalised Associated Linear Equation (GALE) decompositions of the system model. It was observed that by using the GALE method, the number of numerical simulations needed to determine the OFRF of the actual system was considerably reduced. Subsequently, the GALE-generated OFRF representation of the actual system was verified to ensure it clearly represented the system output spectra.

The effect of a hardening stiffness nonlinearity, as well as a nonlinear cubic damping, on the VEH system, was investigated. The application of low-level excitation for this study ensured no jump phenomenon was exhibited. This is because the GALE and OFRF concepts are applicable only to a class of nonlinear systems stable at zero equilibrium and which can be described by a Volterra series model. The exhibition of a jump, by the system model, will nullify the suitability of the methods employed in this study.

The results obtained in this study show that the nonlinear stiffness characteristic extends the operational bandwidth as well as the harvested power of the VEH system, hence improving its performance. It was also revealed that employing a nonlinear cubic damping characteristic, in the presence of a nonlinear stiffness characteristic, improved the performance of the VEH system. Using the OFRF representation, optimal values of the VEH design parameters can be determined for any desired power level within and beyond the design range. Future studies will focus on the design of VEH systems with damping and stiffness nonlinearities subject to mass-displacement constraints inherent in practical VEH systems.