On Efficiency of Two-Degree-of-Freedom Galloping Energy Harvesters with Two Transducers

1. Introduction

In the contemporary world, energy plays a pivotal role across all facets of human existence. Powering electronic gadgets, illuminating and heating homes, and propelling machinery all hinge on energy. However, as the population grow and technologies advance, we increasingly encounter challenges concerning energy production and distribution.

The energy crisis has become a pervasive issue demanding attention. Concurrently, with technological advancements, the concept of the Internet of Things (IoT) is burgeoning [1]. IoT delineates a network of interconnected devices that communicate and share data to enhance the daily lives of users. Smart home gadgets, vehicles, environmental sensors, and myriad other devices can interconnect, analyze data, and autonomously make decisions. Yet, as IoT devices proliferate, challenges in providing them with sustainable energy sources emerge. Conventional power methods such as batteries or power cords may prove inconvenient, have limited durability, or be impractical for distributed IoT networks.

Energy harvesting technologies, such as flow-induced vibration recovery, offer potential solutions. Notably, Bladeless Vortex power plants, activated by von Karman vortices, promise a cost reduction of up to 70% compared to traditional wind farms [3]. However, their dependence on specific flow conditions presents limitations.

Galloping emerges as an alternative vibration excitation method. Early studies explore its potential for energy recovery [4], emphasizing the crucial role of aerodynamic coefficients and the device's quasi-stationary nature. Subsequent research [5] delves into analytical, numerical, and experimental analyses of galloping energy harvesters (GEH), examining structural damping's impact and developing electromechanical models.

Due to the fundamental influence of flow body geometry [6], research focuses on optimizing shapes to enhance efficiency. While rhomboid geometries are prone to galloping instability, adding a stream splitter on the opposite side favors power generation [7,8]. Additionally, devices with irregular shapes demonstrate benefits [9,10].

Nonlinear elasticity's influence on GEH performance is extensively studied [11,12]. It is found to enhance efficiency at low flow rates, albeit potentially reducing operating speeds [13,14]. Various types of nonlinearities are assessed for their impact on device power [15], with duffing-type non-linear GEH variants analyzed for effectiveness [16].

Expanding GEH degrees of freedom shows promise [17], with series-connected mass systems and devices extended to multiple degrees explored [18,19]. Analytical solutions for mathematical models of such devices are developed [20], and experimental research on innovative designs, like magnetically coupled beams, shows increased power generation potential [21].

However, research on two-degree-of-freedom systems lacks exploration of utilizing two electromechanical transducers. Addressing this gap, this study aims to assess the efficiency of a GEH with two degrees of freedom. It considers transducer placement between the stationary base and lower mass and between masses, aiming to contribute to enhancing energy harvesting technologies.

2. Efficiency of Reference Variant

While the primary focus of this study lies in analyzing the efficiency of systems with two degrees of freedom (GEH2D), understanding the characteristics of this class of devices can be facilitated by comparing them with a system possessing well-established, as discussed in [3,4], properties - a linear device with one degree of freedom (GEH1D), hereafter referred to as the reference device. Therefore, it is justified to conduct an analysis of the efficiency of this system, particularly as we aim to highlight a certain intriguing feature of this variant, which has not been previously addressed in the literature. According the cited articles, the mathematical model of GEH1D takes the form (1) with parameters detailed in Table 1.

$$\widehat{m}\ddot{z}+\widehat{c}\dot{z}+\widehat{k}z-\theta \xi =F_L=-\frac{1}{2}\widehat{\rho }{\widehat{u}}^2\widehat{h}\left(a_1\frac{\dot{z}}{\widehat{u}}+a_3{\left(\frac{\dot{z}}{\widehat{u}}\right)}^3\right),$$

(1a)

$$\widehat{C_p}\dot{\xi }+\frac{\xi }{\widehat{R}}+\widehat{\theta }\dot{x}=0.$$

(1b)

The efficiency $\eta$ of such a system is defined as the ratio of the average power generated by the device $\widehat{P_g}$ over a time period equal to the period $\widehat{T}$ of the voltage function ξ(τ) to the surface flow power density $\widehat{P_f}$:

$${\eta }^L=\frac{\widehat{P_g}}{\widehat{P_f}}=\frac{\frac{1}{\widehat{T}}\underset{0}{\overset{\widehat{T}}{\int }}\frac{{\mathsf{\xi }}^2\left(\tau \right)}{\widehat{R}}d\tau }{\frac{1}{2}\widehat{\rho }\widehat{h}{\widehat{u}}^3},$$

(2)

$$\begin{gathered} y=\frac{z}{\widehat{h}},v=\xi \frac{\widehat{\theta }}{\widehat{m}\widehat{h}{{\widehat{\omega }}_n}^2},{\widehat{\omega }}_n=\sqrt{\frac{\widehat{k}}{\widehat{m}}},c=\frac{\widehat{c}}{\widehat{m}{\widehat{\omega }}_n},u=\frac{\widehat{u}}{\widehat{h}{\widehat{\omega }}_n},\rho =\widehat{\rho }\frac{{\widehat{h}}^2}{2\widehat{m}},r=\widehat{C_p}\widehat{R}{\widehat{\omega }}_n, \\ \kappa =\frac{{\widehat{\theta }}^2}{\widehat{C_p}\widehat{m}{{\widehat{\omega }}_n}^2},t=\tau {\widehat{\omega }}_n, \end{gathered}$$

(3)

the device dimensionless mathematical model can be written as:

$$\ddot{y}+c\dot{y}+y-v=\rho \left(a_{1}u\dot{y}+a_3\frac{{\dot{y}}^3}{u}\right),$$

(4a)

$$\dot{v}+\frac{v}{r}+\kappa \dot{y}=0.$$

(4b)

Assuming that $T_c$ is the period of the function $v\left(t\right)$, efficiency in terms of dimensionless parameters can be expressed as:

$${\eta }^L=\frac{\frac{1}{T_c}{\int }_0^{T_c}\frac{v^2\left(t\right)}{\kappa r}dt}{\rho u^3}.$$

(5)

To characterize the efficiency ${\eta }^L$ of the system, it is essential to derive the voltage function $v\left(t\right)$ generated by it, which can be accomplished using the harmonic balance method. Consequently, it was assumed that the limit cycle of the system would be described by a set of solutions in the form:

$$y=A_y\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t\right),$$

(6a)

$$v=A_v\cos \left(\omega t+\phi \right),$$

(6b)

where, the four unknown quantities $A_y$, $A_v$, $\omega =\frac{\widehat{\omega }}{{\widehat{\omega }}_n}$, $\phi$ represent the dimensionless vibration amplitude, dimensionless voltage amplitude, dimensionless vibration frequency, and the phase shift between the oscillator vibrations and voltage oscillations, respectively. The parameter $\widehat{\omega }$ denotes the unknown dimensional frequency of the system's vibrations. Substituting solutions of the form (6) into the model (1) leads to the transformation of the system of differential equations into a system of algebraic equations:

$$A_y\left(1-{\omega }^2\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t\right)+A_y^3\frac{\rho {\omega }^3{a}_3}{u}{\mathrm{s}\mathrm{i}\mathrm{n}}^3\left(\omega t\right)+A_y\omega \left(u\rho a_1-c\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega t\right)-A_v\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t+\phi \right)=0,$$

(7a)

$$A_v\left(\frac{\cos \left(\omega t+\phi \right)}{r}-\omega \sin \left(\omega t+\phi \right)\right)-A_y\kappa \omega \mathrm{s}\mathrm{i}\mathrm{n}\left(\omega t\right)=0.$$

(7b)

This condition will be satisfied for every moment of time $t$ if and only if the sum of the coefficients with corresponding time functions equals zero. From equation (7b), the following was deduced:

$$\frac{A_v\left(\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi \right)-r\omega \mathrm{s}\mathrm{i}\mathrm{n}\left(\phi \right)\right)}{r}=0,$$

(8a)

$$\frac{A_v\left(r\omega \mathrm{c}\mathrm{o}\mathrm{s}\left(\phi \right)+\mathrm{s}\mathrm{i}\mathrm{n}\left(\phi \right)\right)+A_y\kappa r\omega }{r}=0.$$

(8b)

Equation (8a) shows that:

$$\mathrm{t}\mathrm{g}\left(\phi \right)=\frac{1}{\omega r},\sin \left(\phi \right)=\frac{1}{\sqrt{{\left(\omega r\right)}^2+1}},\cos \left(\phi \right)=\frac{\omega r}{\sqrt{{\left(\omega r\right)}^2+1}}.$$

(9)

Based on equations (8b) and (9), the expression describing the relationship between the vibration amplitude and the voltage amplitude can be derived:

$$A_v={-A}_y\frac{r\kappa \omega }{\sqrt{1+r^2{\omega }^2}}.$$

(10)

By balancing harmonics of equation (7a), the remaining two algebraic equations, necessary to determine an approximate solution of the device model are obtained. After taking into account the previously derived relations (9) and (10), they take the form:

$$A_y{k}_e^L\omega -A_y\left({\omega }^2-1\right)=0,$$

(11a)

$$A_{y}u\rho \omega a_1+A_y^3\frac{3\rho {\omega }^3{a}_3}{4u}-A_{y}c\omega -A_y{e}^L\omega =0,$$

(11b)

where $k_e^L=\kappa \frac{r^2{\omega }^2}{1+r^2{\omega }^2}$ and $e^L=\kappa \frac{r}{1+r^2{\omega }^2}$ are the piezoelectric stiffness and electrical damping of the linear system. Based on the above system of equations, it can be shown that:

$$\begin{array}{c}{\omega }_1^2\\ {\omega }_2^2\end{array}=\frac{r^2\left(1+\kappa \right)-1\pm \sqrt{4r^2+{\left(r^2\left(1+\kappa \right)-1\right)}^2}}{2r^2},$$

(12a)

$${A_y}^2=\frac{4u\left(c+e^L-u\rho a_1\right)}{3\rho {\omega }^2{a}_3}.$$

(12b)

Equation (12a) shows that ${\omega }_2^2<0$, regardless of the system parameters. Therefore, in the following part of the work notation ${\omega }_1=\omega$ has been adopted. Returning now to the general definition of efficiency (5) and substituting the expressions (6b), (10), (12) into it, we obtain:

$${\eta }^L=\frac{2e^L\left(c+e^L-u\rho a_1\right)}{3u^2{\rho }^2{a}_3}.$$

(13)

Figure 1 illustrates the efficiency characteristics of the system under discussion. The efficiency is depicted as the ratio ${\eta }^L/{\eta }_p^L$ and such representation maintained throughout the study. The figure also includes the analogous characteristic obtained numerically for initial conditions $y\left(0\right)=0.1$, $\dot{y}\left(0\right)=0$, $n\left(0\right)=0$, utilizing the fourth-order Runge-Kutta method. In this figure it can be observed that the system undergoes Hopf bifurcation at certain speed, from now on referred to as the critical flow speed $u_{cr}$,. Its value results directly from (13):

$${\eta }^L=\frac{2e^L\left(c+e^L-u\rho a_1\right)}{3u^2{\rho }^2{a}_3}=0\Rightarrow u=u_{cr}=\frac{c+e^L}{\rho a_1},$$

(14)

Another noteworthy observation from the same set of graphs is that irrespective of the parameter set, there exists a specific nominal flow speed $u=u_p^L$ for which the efficiency ${\eta }^L$ attains a maximum value, denoted as the peak efficiency ${\eta }_p^L$ . These quantities are given by the expressions:

$$\frac{\delta {\eta }^L}{\delta u}=0\Rightarrow u=u_p^L=2\frac{c+e^L}{\rho a_1}=2u_{cr},$$

(15)

$${\eta }_p^L={\eta }^L\left(u_p^L\right)=\frac{-a_1^2{e}^L}{6a_3\left(c+e^L\right)}.$$

(16)

The aforementioned identities were previously derived in [3,4]. However, the remainder of the article discusses entirely original content. An unexplored property of GEH1D, which can be inferred from Figure 1 or deduced from the identity $u_p^L=2u_{cr}$ , is noteworthy. A system with a low critical speed $u_{cr}$ will experience a more pronounced decline in efficiency due to the deviation of the flow speed $u$ at which it operates from the nominal speed $u_p^L$ . Let the measure of this phenomenon be the flow speed bandwidth in which the system efficiency ${\eta }^L$ does not fall below 90% of the maximum efficiency ${\eta }_p^L$, hereafter referred to as the high efficiency band $B^L$. According to the definition, $B^L$ is given as:

$${\eta }^L=0.9{\eta }_p^L\Rightarrow \begin{array}{c}{u}_2\\ u_1\end{array}=\frac{\rho a_1\left(c+e^L\right)\pm \frac{\sqrt{0.4}}{0.9}\rho a_1\left(c+e^L\right)}{\frac{\mathrm{0,9}}{2}{\left(\rho a_1\right)}^2},$$

(17)

$$B^L=u_2-u_1=\sqrt{\frac{160}{81}}\frac{\left(c+e^L\right)}{\rho a_1}\approx 1.4u_{cr}\approx 0.7u_p^L.$$

(18)

The quantities: critical speed ${\mathit{u}}_{\mathit{c}\mathit{r}}$, nominal speed ${\mathit{u}}_{\mathit{p}}^{\mathit{L}}$, peak efficiency ${\mathit{\eta }}_{\mathit{p}}^{\mathit{L}}$ and high efficiency bandwidth ${\mathit{B}}^{\mathit{L}}$ will be further referred to as the criterion parameters. They have been labeled on Fig. 1.

3. Efficiency of Two Degree of Freedom System

Based on the information provided in section 1, the presence of additional mass in systems with two degrees of freedom implies the necessity to consider the placement of the electromechanical transducer. It may be positioned between the main mass and the stationary base (Fig. 2a), between the masses (Fig. 2b), or in both of these locations (Fig. 2c).

Figure 2. Subvariants of two degree of freedom system.

Figure 2. Subvariants of two degree of freedom system.

To evaluate the influence of piezoelectric placement on efficiency, the efficiency characteristics of the variant with both transducers will be derived. Subsequently, the ramifications of removing one of them will be examined.

The dynamics of devices within the discussed category can be elucidated by a general dimensionless mathematical model of the form (20), where, in addition to the identities (3), the following also holds:

$$\begin{gathered} y_i=\frac{z_i}{\widehat{h}},v_i=\xi \frac{\widehat{{\theta }_i}}{\widehat{m_i}\widehat{h}{{\widehat{\omega }}_n}^2},c_i=\frac{{\widehat{c}}_i}{{\widehat{m}}_1{\widehat{\omega }}_n},k_2=\frac{{\widehat{k}}_2}{{\widehat{k}}_1},\kappa =\frac{{\widehat{{\theta }_i}}^2}{\widehat{C_{pi}}\widehat{m_1}{{\widehat{\omega }}_n}^2}, \\ \theta =\frac{{\widehat{\theta }}_1}{{\widehat{\theta }}_2},M=\frac{{\widehat{m}}_2}{{\widehat{m}}_1},i=\mathrm{1,2} \end{gathered}$$

(19)

$${\ddot{y}}_1+y-k_2\left(y_2-y_1\right)+c_1{\dot{y}}_1-c_2\left({\dot{y}}_2-{\dot{y}}_1\right)-v_1+v_2=\rho \left(a_{1}u{\dot{y}}_1+a_3\frac{{\dot{y}}_1^3}{u}\right),$$

(20a)

$${\kappa }_1{\dot{v}}_1+\frac{v_1}{r_1}+{\dot{y}}_1,$$

(20b)

$${M\ddot{y}}_2+k_2\left(y_2-y_1\right)+c_2\left({\dot{y}}_2-{\dot{y}}_1\right)-v_2=0,$$

(20c)

$${\kappa }_2{\dot{v}}_2+\frac{v_2}{r_2}+{\dot{y}}_2-{\dot{y}}_1=0.$$

(20d)

To derive the solutions of the GEH2D mathematical model, the procedure outlined in [20] was adapted. It was assumed that the approximate solution of the model (20) will have the following form:

$$y_1=A_1\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{1}t\right)+B_1\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{1}t\right)+G_1\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{2}t\right)+H_1\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{2}t\right),$$

(21a)

$$v_1=n_1\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{1}t\right)+n_2\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{1}t\right)+n_3\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{2}t\right)+n_4\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{2}t\right),$$

(21b)

$$y_2=A_2\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{1}t\right)+B_2\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{1}t\right)+G_2\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{2}t\right)+H_2\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{2}t\right),$$

(21c)

$$v_2={\vartheta }_1\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{1}t\right)+{\vartheta }_2\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{1}t\right)+{\vartheta }_3\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{2}t\right)+{\vartheta }_4\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{2}t\right),$$

(21d)

Substituting the above identities into equations (20b) and (20d), and then balancing the harmonics $\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{1}t\right)$, $\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{1}t\right)$, $\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{2}t\right)$, $\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{2}t\right)$, allows for the derivation of the relationships between the amplitudes of voltage and the amplitudes of vibration.

$$\begin{gathered} n_1=B_1{\omega }_1{\epsilon }_{D11}-A_1{\delta }_{D11},n_2=-A_1{\omega }_1{\epsilon }_{D11}-B_1{\delta }_{D11}, \\ {n}_3=H_1{\omega }_2{\epsilon }_{D12}-G_1{\delta }_{D12},n_4=-G_1{\omega }_2{\epsilon }_{D12}-H_1{\delta }_{D12}, \end{gathered}$$

(22a)

$$\begin{gathered} {\vartheta }_1=\theta \left({\epsilon }_{D21}\left(\theta B_2-B_1\right){\omega }_1+{\delta }_{D21}\left(A_1-\theta A_2\right)\right), \\ {\vartheta }_2=\theta \left({\epsilon }_{D21}\left(A_1-\theta A_2\right){\omega }_1+{\delta }_{D21}\left(B_1-\theta B_2\right)\right), \\ {\vartheta }_3=\theta \left({\epsilon }_{D22}\left(\theta H_2-H_1\right){\omega }_2+{\delta }_{D22}\left(G_1-\theta G_2\right)\right), \\ {\vartheta }_4=\theta \left({\epsilon }_{D22}\left(G_1-\theta G_2\right){\omega }_2+{\delta }_{D22}\left(H_1-\theta H_2\right)\right). \end{gathered}$$

(22b)

where: ${\delta }_{Dij}={\kappa }_i\frac{r^2{\omega }_{\mathrm{j}}^2}{1+r^2{\omega }_{\mathrm{j}}^2}$ and ${\epsilon }_{Dij}={\kappa }_i\frac{r}{1+r^2{\omega }_{\mathrm{j}}^2}$ represent the piezoelectric stiffness and electric damping of the $i$-th piezoelectric for vibration with the $j$-th frequency. Harmonic balancing of $\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{1}t\right)$ and $\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{1}t\right)$ of the algebraic equation obtained by substituting the solutions (21) into equation (20c) and using the identities (22) leads to the system of equations:

$$A_2=B_1{p}_1+A_1{q}_1,$$

(23a)

$$B_2=B_1{q}_1-A_1{p}_1,$$

(23b)

where:

$$p_1=\frac{(k_2+\theta {\delta }_{D21})(c_2+{\theta }^2{\epsilon }_{D21}){\omega }_1-(c_2+\theta {\epsilon }_{D21})(k_2+{\theta }^2{\delta }_{D21}-M{\omega }_1^2){\omega }_1}{{(c_2+{\theta }^2{\epsilon }_{D21})}^2{\omega }_1^2+{(k_2+{\theta }^2{\delta }_{D21}-M{\omega }_1^2)}^2},$$

(23c)

$$q_1=\frac{(c_2+\theta {\epsilon }_{D21})(c_2+{\theta }^2{\epsilon }_{D21}){\omega }_1^2+(k_2+\theta {\delta }_{D21})(k_2+{\theta }^2{\delta }_{T21}-M{\omega }_1^2)}{{(c_2+{\theta }^2{\epsilon }_{D21})}^2{\omega }_1^2+{(k_2+{\theta }^2{\delta }_{D21}-M{\omega }_1^2)}^2}.$$

(23d)

By substituting the above identities into the harmonic balance equations for $\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{1}t\right)$ and $\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{1}t\right)$ of expression (20a), the following identities can be formulated:

$$A_1{\alpha }_1+B_1{\beta }_1-A_1\frac{3\rho a_3{\omega }_1\left(\right(A_2^2+B_2^2){\omega }_1^2+2(G_2^2+H_2^2\left){\omega }_2^2\right)}{4u}=0,$$

(24a)

$$A_1{\beta }_1-B_1{\alpha }_1+B_1\frac{3\rho a_3{\omega }_1\left(\left(A_2^2+B_2^2\right){\omega }_1^2+2\left(G_2^2+H_2^2\right){\omega }_2^2\right)}{4u}=0,$$

(24b)

where:

$${\alpha }_1=\left(p_1\right(k_2+{\theta }^2{\delta }_{D21})+(c_1+{\epsilon }_{T11}+\left(c_2+{\epsilon }_{D21}\right)(1-\theta q_1)){\omega }_1-u\rho {\omega }_1{a}_1),$$

(25a)

$${\beta }_1=1+{\delta }_{T11}+\left({\delta }_{D21}\theta +k_2\right)\left(1-\theta q_1\right)-{\omega }_1\left(p_1\left(c_2+{\theta }^2{\epsilon }_{D21}\right)+{\omega }_1\right).$$

(25b)

Adapting the above procedure to balance the harmonics $\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{2}t\right)$ and $\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{2}t\right)$ of equations (20c) and (20a) allows for the derivation of the following identities:

$$G_2=H_1{p}_2+G_1{q}_2,$$

(26a)

$$H_2=H_1{q}_2-G_1{p}_2,$$

(26b)

$$G_1{\alpha }_2+H_1{\beta }_2-G_1\frac{3\rho a_3{\omega }_2\left(2\right(A_1^2+B_1^2){\omega }_1^2+(G_1^2+H_1^2\left){\omega }_2^2\right)}{4u}=0$$

(26c)

$$G_1{\beta }_2-H_1{\alpha }_2+H_1\frac{3\rho a_3{\omega }_2\left(2\right(A_1^2+B_1^2){\omega }_2^2+(G_1^2+H_1^2\left){\omega }_2^2\right)}{4u}=0$$

(26d)

where:

$$p_2=\frac{(k_2+\theta {\delta }_{D22})(c_2+{\theta }^2{\epsilon }_{D22}){\omega }_2-(c_2+\theta {\epsilon }_{D22})(k_2+{\theta }^2{\delta }_{D22}-M{\omega }_2^2){\omega }_2}{{(c_2+{\theta }^2{\epsilon }_{D22})}^2{\omega }_2^2+{(k_2+{\theta }^2{\delta }_{D22}-M{\omega }_2^2)}^2},$$

(27a)

$$q_2=\frac{(c_2+\theta {\epsilon }_{T22})(c_2+{\theta }^2{\epsilon }_{D22}){\omega }_2^2+(k_2+\theta {\delta }_{D22})(k_2+{\theta }^2{\delta }_{D22}-M{\omega }_2^2)}{{(c_2+{\theta }^2{\epsilon }_{D22})}^2{\omega }_2^2+{(k_2+{\theta }^2{\delta }_{D22}-M{\omega }_2^2)}^2},$$

(27b)

$${\alpha }_2=\left(p_2\right(k_2+{\theta }^2{\delta }_{D22})+(c_1+{\epsilon }_{T12}+\left(c_2+{\epsilon }_{T22}\right)(1-\theta q_2)){\omega }_2-u\rho {\omega }_2{a}_1),$$

(27c)

$${\beta }_2=1+{\delta }_{D12}+\left({\delta }_{D22}\theta +k_2\right)\left(1-\theta q_2\right)-{\omega }_2\left(p_2\left(c_2+{\theta }^2{\epsilon }_{D22}\right)+{\omega }_2\right).$$

(27d)

By adding equation (24a) multiplied by $B_1$ to equation (24b) multiplied by $A_1$, an expression was obtained that allows the frequency ${\omega }_1$ to be explicitly formulated:

$$1+k_2+{\delta }_{D11}-k_2{q}_1+{\delta }_{D21}\theta \left(1-\theta q_1\right)-p_1\left(c_2+{\theta }^2{\epsilon }_{D21}\right){\omega }_1-{\omega }_1^2=0.$$

(28)

Similarly, by adding equation (26c) multiplied by $H_1$ to equation (26d) multiplied by $G_1$, an equation was obtained from which the frequency ${\omega }_2$ can be derived:

$$1+k_2+{\delta }_{D12}-k_2{q}_2+{\delta }_{D22}\theta \left(1-\theta q_2\right)-p_2\left(c_2+{\theta }^2{\epsilon }_{D22}\right){\omega }_2-{\omega }_2^2=0.$$

(29)

Note that in each pair of parameters (($p_1$, $p_2$), ($q_1$, $q_2$), (${\epsilon }_{D21},{\epsilon }_{D22}$), (${\delta }_{D11},{\delta }_{D12}$), (${\delta }_{D21},{\delta }_{D22}$), the only difference is the frequency in their definition – ${\omega }_1$ or ${\omega }_2$. Considering the similarity between equations (28) and (29), one can conclude that the frequencies ${\omega }_1$ and ${\omega }_2$ must be equal, thereby excluding the possibility of polymodal vibrations in the system. The expression describing the vibration frequency ${\omega }_2={\omega }_1$ can therefore be derived by solving only one of the above equations.

The relationship between the vibration amplitude and the system frequency remains unknown. The first of the equations necessary to determine this relationship was obtained by subtracting equation (24b) multiplied by $B_1$ from equation (24a) multiplied by $A_1$. The second one is the difference of equation (26c) multiplied by $G_1$ and equation (26d) multiplied by $H_1$:

$$\begin{array}{l}{p}_1\left(k_2+{\theta }^2{\delta }_{D21}\right)+\left(c_1+c_2\left(1-q_1\right)+{\epsilon }_{D11}+{\epsilon }_{D21}\theta \left(1-\theta q_1\right)\right){\omega }_1\\ -u\rho {\omega }_1{a}_1-\frac{3\rho {\omega }_1\left(A_y^2{\omega }_1^2+2G_y^2{\omega }_2^2\right)a_3}{4u}=0,\end{array}$$

(30a)

$$\begin{array}{l}{p}_2\left(k_2+{\theta }^2{\delta }_{D22}\right)+\left(c_1+c_2\left(1-q_1\right)+{\epsilon }_{D12}+{\epsilon }_{D22}\theta \left(1-\theta q_2\right)\right){\omega }_2\\ -u\rho {\omega }_2{a}_1-\frac{3\rho {\omega }_2\left(2A_y^2{\omega }_1^2+G_y^2{\omega }_2^2\right)a_3}{4u}=0,\end{array}$$

(30b)

where $A_y^2=A_1^2+B_1^2$ and $G_y^2=G_1^2+H_1^2$. Set of equations (30) has three non-trivial solutions in terms of $A_y^2$ and $G_y^2$, which after recalling the ${\omega }_2={\omega }_1$ identity, take the form::

$$\begin{gathered} A_y^2=\frac{4u\left(p_1\left(k_2+{\theta }^2{\delta }_{D21}\right)+\left(c_1+c_2\left(1-q_1\right)+{\epsilon }_{D11}+{\epsilon }_{D21}\theta \left(1-\theta q_1\right)\right)-u\rho a_1\right)}{3\rho {{\omega }_1}^3{a}_3}, \\ {G}_y^2=0 \end{gathered}$$

(31a)

$$\begin{gathered} A_y^2=0, \\ {G}_y^2=\frac{4u\left(p_1\left(k_2+{\theta }^2{\delta }_{D21}\right){\omega }_1+\left(c_1+c_2\left(1-q_1\right)+{\epsilon }_{D11}+{\epsilon }_{D21}\theta \left(1-\theta q_1\right)\right)-u\rho a_1\right)}{3\rho {{\omega }_1}^3{a}_3}, \end{gathered}$$

(31b)

$$\begin{gathered} A_y^2=\frac{4u\left(3p_1\left(k_2+{\theta }^2{\delta }_{D21}\right){\omega }_1+c_1+c_2\left(1-q_1\right)+{\epsilon }_{D11}+{\epsilon }_{D21}\theta \left(1-\theta q_1\right)\right)-u\rho a_1}{9\rho {{\omega }_1}^2{a}_3}, \\ {G}_y^2=\frac{4u\left(-3p_1\left(k_2+{\theta }^2{\delta }_{D21}\right){\omega }_1+c_1+c_2\left(1-q_1\right)+{\epsilon }_{T11}+{\epsilon }_{D21}\theta \left(1-\theta q_1\right)\right)-u\rho a_1}{9\rho {{\omega }_1}^2{a}_3}. \end{gathered}$$

(31c)

These expressions, along with the previously derived identities, enable the explicit formulation of solutions for the system (20) in the form (21). It should be noted that the identities given by equations (31a) and (31b) correspond to the same solution in the form (21); thus, only one of them, namely (31a), will be further analyzed. Moreover, no set of system parameters and initial conditions has been found that would lead to the excitation of vibrations with amplitude (31c). Therefore, this solution was considered unstable, and the efficiency has been derived based only on (31a):

$${\eta }_D=\frac{\frac{1}{T_c}\underset{0}{\overset{T_c}{\int }}\left(\frac{{v_1}^2\left(t\right)}{r_1{\kappa }_1}+\frac{{v_2}^2\left(t\right)}{r_2{\kappa }_2}\right)dt}{\rho u^3}=\frac{2\lambda (\mu -u\rho a_1)}{3u^2{\rho }^2{a}_3},$$

(32)

where:

$$\lambda =({\epsilon }_{D11}+{\epsilon }_{D12}{\theta }^2(1-2\theta q_1+{\theta }^2({p_1}^2+{q_1}^2)\left)\right),$$

(33a)

$$\mu =\left(c_1+\left(c_2+\theta {\epsilon }_{D11}\right)\left(1-\theta q_1\right)\right)+\frac{p_1\left(k_2+{\theta }^2{\delta }_{TD21}\right)}{{\text{ω}}_1}.$$

(33b)

Critical speed $u_{cr}^D$ and nominal speed $u_p^D$ of GEH2D are given by the expressions:

$${\eta }_D={\mathsf{\gamma }}_1\frac{2\lambda (\mu -u\rho a_1)}{3u^2{\rho }^2{a}_3}=0\Rightarrow u=u_{cr}^D=\frac{\mu }{\rho a_1},$$

(34)

$$\frac{\delta {\eta }_D}{\delta u}=0\Rightarrow u=u_p^D=2\frac{\mu }{\rho a_1}=2u_{cr}^D.$$

(35)

Now it is possible to derive peak efficiency of the system:

$${\eta }_p^D={\eta }^D\left(u_p^D\right)=\frac{-\lambda a_1^2}{6{\mathsf{\gamma }}_2{a}_3}.$$

(36)

To fully determine the features of the analyzed variant, it is necessary to define its high efficiency bandwidth $B^D$. It is:

$${\eta }^D=\mathrm{0,9}{\eta }_p^D\Rightarrow \begin{array}{c}{u}_2\\ u_1\end{array}=\frac{\frac{2}{0.9}\rho a_1\mu \pm \frac{\sqrt{0.4}}{0.9}\rho a_1\mu }{{\left(\rho a_1\right)}^2},$$

(37)

$$B^D=u_2-u_1=\sqrt{\frac{160}{81}}\frac{\mu }{\rho a_1}\approx 1.4u_{cr}^D\approx 0.7u_p^D.$$

(38)

Figure 3 depicts the relationship between efficiency ${\eta }^D$ and flow speed $u$, represented by the function (32), compared with an analogous relationship obtained through numerical integration of the model (20) for initial conditions: $y_1\left(0\right)=0.1$, ${\dot{y}}_1\left(0\right)=0$, $n_1\left(0\right)=0$, $y_2\left(0\right)=0$, ${\dot{y}}_2\left(0\right)=0$, $n_2\left(0\right)=0$. Moreover, the figure shows the values characterizing the efficiency of the variant - $u_{cr}^D$, $u_p^D$, ${\eta }_p^D$ and $B^I$.

4. Comparison of Systems

As suggested by equations (34), (35), and (38), critical speeds, nominal speed, and high-efficiency bandwidths are interrelated in the same manner in both analyzed GEH variants. Consequently, comparing the efficiency of these systems can only be achieved by comparing their peak efficiencies. Since the efficiency of both devices depends on their nominal speeds, it is necessary to compare systems with the same nominal speeds. The electrical damping of the linear system $e^L$ will therefore be selected in such a way that this condition is satisfied. In the following derivation, it is assumed that the structural damping of the linear system c is equal to the structural damping $c_1$ of the damper connecting the lower mass of the GEH2D with the base, i.e.: $c=c_1$. The shape of the flowing body, represented by the coefficients $a_1$ and $a_3$, and the density of the fluid remain the same for all variants.

$$u_p^I=u_p^L,$$

(39a)

$$2\frac{\mu +{\mathsf{\gamma }}_2}{\rho a_1}=2\frac{e^L+c_1}{\rho a_1},$$

(39b)

$$e^L={\mathsf{\gamma }}_2+\mu -c_1.$$

(39c)

According to equations (16) and (39c), the efficiency of a system with one degree of freedom with a nominal speed equal to the nominal speed of the GEH2D is equal to:

$${\eta }^L\left(u_p^D\right)=-\frac{\left(\mu -c_1\right)a_1^2}{6\mu a_3}.$$

(40)

The ratio of peak efficiencies of compared systems with the same nominal speed is therefore given by:

$$\frac{{\eta }^L\left(u_p^D\right)}{{\eta }^D\left(u_p^D\right)}=\frac{s_1+s_2}{s_1+s_3},$$

(41)

where:

$$s_1={\epsilon }_{T11}\left(\right(k_2+{\theta }^2{\delta }_{D21}\left)\right(k_2+{\theta }^2{\delta }_{T21}-2M{{\omega }_1}^2)+{(c_2+{\theta }^2{\epsilon }_{D21})}^2{{\omega }_1}^2+M^2{{\omega }_1}^4),$$

(42a)

$$s_2=M{{\omega }_1}^2\left((k_2{\epsilon }_{D21}-c_2{\delta }_{D21})(\theta -1)\theta +M{{\omega }_1}^2(c_2+\theta {\epsilon }_{D21})\right),$$

(42b)

$$s_3={\theta }^2{\epsilon }_{D21}\left(M{{\omega }_1}^2\left(2k_2\right(\theta -1)+M{{\omega }_1}^2)+({c_2}^2{{\omega }_1}^2+{k_2}^2){(\theta -1)}^2\right).$$

(42c)

Clearly, according to (41), discussed system with two degree of freedom, will have greater efficiency than reference system if and only if $s_3>s_2$. For variant with only one, lower transducer (Fig. 2a), where ${\epsilon }_{D21}={\delta }_{D21}=0$, this inequality takes the form:

$$0\ngtr c_2{M}^2{{\omega }_1}^4,$$

(43)

what implies that such a device cannot be more efficient than the reference system, regardless of its parameters. This conclusion contradicts the results presented in [17]. The reason for this discrepancy lies in the fact that in the cited work, the operating conditions of the compared systems were not standardized - the devices had different critical speeds and, consequently, different nominal speeds. Despite the unquestionable value of this article, the conclusion stated therein can be subject to questioning. In the case of the variant with only the upper transducer (Fig. 2b), where ${\epsilon }_{T11}={\delta }_{T11}=0$, the condition for efficiency improvement takes the following form:

$${\epsilon }_{D21}\ngtr {\epsilon }_{D21}+c_2,$$

(44)

what indicates that this variant does not offer efficiency enhancement. Similarly, an identical and impossible-to-satisfy requirement is associated with the special case of the third variant, where there are two identical transducers, meaning: ${\epsilon }_{D11}={\epsilon }_{T21}$, ${\delta }_{D11}={\delta }_{D21}$ and $\theta =1$. However, the inequality $s_3>s_2$ can be satisfied for the most general variant – the one equipped with two different transducers. In this case, it can be reduced to the condition:

$$\begin{gathered} \left(\theta -1\right)\theta \left({\epsilon }_{D21}{k_2}^2\left(\theta -1\right)\theta +c_2{\delta }_{D21}M{{\omega }_1}^2+{c_2}^2{\epsilon }_{D21}{{\omega }_1}^2\left(\theta -1\right)\theta +{\epsilon }_{D21}{k}_{2}M{{\omega }_1}^2\left(2\theta -1\right)\right) \\ >M^2{{\omega }_1}^4(c_2-{\epsilon }_{D21}(\theta -1\left)\theta \right), \end{gathered}$$

(45)

which, after substituting the solution of (27), can be solved numerically for a chosen parameter. For the parameters presented in Table 2, solving the above inequality for $\theta$ results in obtaining a threshold value of approximately $\theta \approx 1.13$. An example of the efficiency characteristics of the GEH1D with $e^L$ given by (39c) and GEH2D is shown in Figure 4, where condition (45) has been met by adopting parameters according to Table 2, with $\theta =1.3>1.13$.

5. Conclusions

The aim of the study was to investigate the efficiency of a galloping energy harvester with two degrees of freedom, considering its three variants, which differ in the number and location of electromechanical transducers. The realization of this objective commenced with the analysis of the reference variant with one degree of freedom. Utilizing the harmonic balance method, an approximate solution of the mathematical model of the system was derived, followed by the formulation of an expression describing the efficiency of the variant. Based on this, key criteria parameters were defined, providing comprehensive information about the variant's efficiency: peak efficiency, high efficiency bandwidth, critical speed, and nominal speed.

Subsequent sections of the study delineate the multitude of configurations that a system with two degrees of freedom can adopt: three different sub-variants of the device were characterized, differing from each other by the location and number of the electromechanical transducers. Parameters characterizing the efficiency of all subvariants were derived, demonstrating that the critical speed, nominal speed, and high efficiency bandwidth are related to each other in the same manner as in the case of the reference system. It was then demonstrated that among the three indicated subvariants of the two-degree-of-freedom system, only the one with two transducers can be more efficient than the reference system.