An investigation on the use of Galerkin and assumed mode methods for analyzing the dynamic behavior of piezoelectric energy harvester under Vortex induced vibration

Abstract

In this paper a new mathematical modeling for two configurations of piezoelectric energy harvesters under vortex induced vibration is presented. At the first (second) configuration, a cylinder under vortex flow is attached to the end of the cantilever beam as parallel (vertical) to the beam. In contrast to previous works that for applying the Galerkin method, the common cross section of the beam and the cylinder has been considered as the boundary of the configuration, here, a new mathematical modeling is presented which considers the boundary conditions of system at the end point of the configuration. It causes that the complexity of system is removed since the boundary conditions will be as geometric. Also, it causes that the notation of modeling for both directions of the cylinder that is, vertical parallel directions be identical. The coupling between the structure and wind flow are considered using wake oscillator Van der Pol model. The discretized equations of motion are solved using both numerical and perturbation method. The perturbation method presented in this paper is a new technique for the configuration with vertical cylinder. The results show that the perturbation method can predict the dynamic behavior of system correctly. For more verification the system is discretized according the assumed mode method presented in previous work. The comparison between Galerkin method and assumed mode method shows a good agreement.

Keywords

Vortex induced vibration piezoelectric energy harvester perturbation technique Galerkin method assumed mode method

1. Introduction

When a wind flow passes over a cylinder, the boundary layer separation often occurs. This separation causes a fluid rotation and vortex shedding behind the cylinder. Depending on the value of Reynolds number, the vortex pattern shedding might result periodic. The periodic pattern causes a periodic lift force on the cylinder. If the cylinder be attached to the end of a flexible beam, then the periodic lift force will oscillate the beam. This type of vibration is called vortex induced vibration or VIV (Liang et al., 2021). When the vortex shedding frequency approaches the natural frequency of the oscillator, a significant increase of the value of oscillations amplitude takes place. This phenomenon is called lock-in phenomena. Recently, the use of piezoelectric layer for harvesting the electric energy from a beam under vortex induced vibration has been noticed. In this configuration a piezoelectric layer coated on the beam is subjected to fluctuating stress due to the fluctuating bending of the beam, and it causes that the electricity is generated (Wang et al., 2020).

The position and direction of the attached cylinder to the beam may be categorized into four main groups. At the first group, the cylinder is attached to the middle part of a clamped-clamped beam as a dumbbell (Alimanesh and Zamanian, 2024; Su and Lin, 2020; Su and Wang, 2021). At the second group, the cylinder is attached to the end of a cantilever beam as vertical to the beam (Akaydin et al., 2012; Dai et al., 2014; Hou et al., 2020; Huang and Zhong, 2023; Javed and Abdelkefi, 2019; Lai et al., 2021), and at the third group it is attached to the cantilever beam as parallel that is, it is attached along the length of the beam (Dai et al., 2016; Jia et al., 2018; Sun et al., 2019a; Zhang et al., 2017a; Zhao et al., 2019). At the forth group the shape of the cylinder has been changed for improving the performance of aeroeleastic energy harvester, for example the shape of the cylinder has been assumed as a combination of the circular and other shapes such as cubic (Rezaei and Talebitooti, 2019; Sun et al., 2019b). In some works, the cylinder has been attached as asymmetric to the end of the beam (Shan et al., 2017; Sui et al., 2022), or an interface cylinder has been seated near the cylinder (Tang et al., 2022; Zhang et al., 2017b).

In previous works the equation of motion have been obtained in discretized form using two procedures. At the first procedure the energy relations of the system have been discretized using assumed modes method, and then the ordinary differential equations of motion have been obtained applying the Euler-Lagrange equation (Hou et al., 2020; Javed and Abdelkefi, 2019; Jia et al., 2018; Shan et al., 2017; Sui et al., 2022). At the second procedure, firstly the partial differential equation of system have been derived using Newton’s second law or Hamilton principle, and then the equations of motion have been obtained in discretized form using Galerkin method (Alimanesh and Zamanian, 2024; Lai et al., 2021; Rezaei and Talebitooti, 2019; Sun et al., 2019a).

The literature review of the previous works indicates that both procedure of assumed modes method and Galerkin method have been used for discretizing the motion equations of a beam under vortex induced vibration. It states that for applying the Galerkin method, the common cross section of the beam and the cylinder has been considered as the boundary of the configuration. It causes the complexity of the problem due to the dynamic nature of the boundary conditions. In this paper a new mathematical modelling is presented which considers the cylinders as a part of the main beam, and considers the boundary conditions of system at the end point of the configuration. It causes that the complexity of system is removed since the boundary conditions will be as geometric. Also, it causes that the notation of modelling for both directions of the cylinder that is, vertical and parallel directions be identical. In other words this paper is an attempt to present a deep understanding of the use of Galerkin method for both directions of the cylinder. For verification of the results, the equations of motion are also obtained in discretized form by applying the assumed modes procedure used in previous work. The discretized form of the equations are solved using Rkf45 numerical method, and the results of two procedure are compared together. For more verification the equations of motion are solved using multiple scale perturbation technique which has not been presented in previous works for vertical direction of the cylinder. It means that this paper is also an attempts to present an analytical perturbation solution when the cylinder is attached to the beam as vertical.

2. Modeling and formulation

The considered configurations consist of a cylinder attached to the free end of a flexible beam, oriented parallel (vertically) to the beam, as illustrated in Figure 1(a)–(c). These configurations are hereafter referred to as Configuration 1 and Configuration 2, respectively. It is assumed that a portion of the beam is coated with a piezoelectric layer. The cylinder is exposed to a steady wind flow with $300 < Re < 1.5 \times 10^{5}$ . If the XYZ coordinate system is defined as shown in Figure 1, the structure undergoes vibrations in the Y-direction due to the periodic nature of vortex shedding behind the cylinder. This periodic excitation induces mechanical strain in the beam, which in turn generates an electric potential difference across the two surfaces of the piezoelectric layer. Previous studies have demonstrated that when an electric voltage is applied across the surfaces of a rectangular piezoelectric layer bonded to the surface of a beam, the piezoelectric effect causes the layer to tend to change its length. However, since the layer is constrained by the underlying substrate, this tendency results in the generation of mechanical bending moments at the two ends of the piezoelectric layer, which are transferred to the substrate. Consequently, the voltage harvested by the piezoelectric layer due to the beam’s vibration acts as an actuating input to the beam itself. These induced bending moments at the beginning and end of the piezoelectric layer are denoted by $M_{1 P}$ and $M_{2 P}$ , respectively, as shown in Figure 1 (Hagedorn and DasGupta, 2007; Dai et al., 2014; Firouzi et al., 2016; Jia et al., 2018). Based on the findings of these previous works, the dynamic behavior of the system is governed by the following equation of motion:

\begin{matrix} \frac{\partial^{2}}{\partial X^{2}} (I_{E} (X) \frac{\partial^{2} W (X, t)}{\partial X^{2}}) + c (X) \frac{\partial W (X, t)}{\partial t} \\ + m (X) \frac{\partial^{2} W (X, t)}{\partial t^{2}} - \frac{\partial}{\partial X} (I_{m} (X) \frac{\partial^{3} W (X, t)}{\partial t^{2} \partial X}) \\ = F_{L} (X, t) + \frac{\partial M_{P} (X, t)}{\partial X} \\ W (0, t) = 0, \frac{\partial W (0, t)}{\partial X} = 0, \frac{\partial^{3} W (L_{e}, t)}{\partial X^{3}} = 0, \\ \frac{\partial^{2} W (L_{e}, t)}{\partial X^{2}} = 0 \\ L_{e} = {\begin{matrix} L_{s} + L_{c} \to for configuration 1, \\ L_{s} + D \to for configuration 2 \end{matrix} \end{matrix}

(1)

Where $W$ indicates the transverse deflection of the system along the Y-direction, and $t$ represents the time variable. Moreover, $m (X)$ denotes the mass per unit length of the configuration, $I_{E} (X)$ is the bending stiffness, $I_{m} (X)$ is the mass moment of inertia per unit length of the configuration, and $c (X)$ is the equivalent damping coefficient per unit length of the configuration due to viscosity and wind flow. $F_{L} (X, t)$ is the external force per unit length acting on the configuration due to vortex shedding effects, as illustrated in Figure 1. The term $M_{P} (X, t)$ represents the mechanical torque per unit length of the beam induced by the piezoelectric effect. Additionally, $L_{s}$ is the length of the substrate layer of the beam, $L_{c}$ is the length of the cylinder, and $D$ is the diameter of the cylinder. $ρ_{c}$ indicates the mass density of the cylinder, and $L_{e}$ is the total length of the configurations.

The function of $I_{E} (X)$ and $m (X)$ may be described as below by the use of Heaviside function $H_{L_{i}}$ :

\begin{array}{l} I_{E} (X) = I_{s} (H (X) - H (X - L_{1 p})) + I_{s p} (H (X - L_{1 p}) - H (X - L_{2 p})) + I_{s} (H (X - L_{2 p}) - H (X - L_{s})) + \\ I_{c} (H (X - L_{s}) - H (X - L_{e})), \\ m (X) = m_{s} (H (X) - H (X - L_{1 p})) + m_{s p} (H (X - L_{1 p}) - H (X - L_{2 p})) + m_{s} (H (X - L_{2 p}) - H (X - L_{s})) + \\ m_{c} (H (X - L_{s}) - H (X - L_{e})), \\ I_{m} (X) = I_{m c} (H (X - L_{s}) - H (X - L_{e})) \\ H_{L_{i}} = H e a v i s i d e (X - L_{i}) = {\begin{cases} 0 X < L_{i}, \\ 1 X \geq L_{i} \end{cases} \end{array}

(2)

where $I_{s}$ and $I_{sp}$ are the bending stiffness coefficient of single layered part and double layered part of system, respectively. Moreover $I_{c}$ is the bending stiffness coefficient of the cylinder part. Moreover $H_{L_{i}}$ is the Heaviside function described as the above, $L_{1 p}$ and $L_{2 p}$ are the distance of the left and right sides of the piezoelectric layer from the left clamped side of the beam, respectively. In addition $m_{s}$ and $m_{sp}$ are the mass per unit length of the configuration at the single and double layered part of beam, respectively, and $m_{c}$ is the mass of the cylinder per unit the length of axial X. Considering Figure 2, the amount of $m_{c}$ and $m_{s}$ and $m_{sp}$ may be evaluated as below.

\begin{matrix} m_{s} = ρ_{s} h_{s} b_{s}, m_{sp} = ρ_{s} h_{s} b_{s} + ρ_{p} h_{p} b_{p}, \\ m_{c} = ρ_{c} \frac{π}{4} D^{2} \to for configuration 1, \\ m_{c} = ρ_{c} L_{c} h_{c} (X) \to for configuration 2, \\ h_{c} (X) = 2 \sqrt{{(\frac{D}{2})}^{2} - {(\frac{D}{2} + L_{s} - X)}^{2}} + h_{s} \end{matrix}

(3)

In the above equation, $ρ_{s}$ , $ρ_{p}$ , and $ρ_{c}$ are the mass density of the substrate layer, piezoelectric layer and cylinder, respectively. Moreover $h_{s}$ , $h_{p}$ , and $h_{c}$ are the thickness of the substrate layer, thickness of piezoelectric layer and the height in the cross section of the cylinder which is vertical to axis X as shown in Figure 1. Moreover $b_{s}$ and $b_{p}$ are the width of the substrate layer and the piezoelectric layer. Considering Figure 1 the amount of $I_{c}$ and $I_{s}$ and $I_{sp}$ may be evaluated as below:

\begin{matrix} I_{s} = \frac{1}{12} E_{s} b_{s} h_{s}^{3} \\ I_{sp} (X) = (\frac{1}{3} E_{s} b_{s} (y_{1}^{3} - y_{0}^{3}) + \frac{1}{3} E_{p} b_{p} (y_{2}^{3} - y_{1}^{3})), \\ I_{c} = E_{c} \frac{π}{4} {(\frac{D}{2})}^{4} \to for configuration 1, \\ I_{c} = E_{c} \frac{1}{12} L_{c} {h_{c}}^{3} \to for configuration 2, \\ I_{mc} = ρ_{c} \frac{π}{4} {(\frac{D}{2})}^{4} \to for configuration 1, \\ I_{mc} = ρ_{c} \frac{1}{12} L_{c} {h_{c}}^{3} \to for configuration 2 \end{matrix}

(4)

Where $E_{s}$ , $E_{p}$ , and $E_{c}$ are the elasticity modulus of the substrate layer, piezoelectric layer and the cylinder, respectively. Considering Figure 2, $y_{1}$ and $y_{2}$ are the distance of lower and upper surface of the substrate layer and piezoelectric layer from neutrals axis in double layered part of system which are evaluated as:

\begin{matrix} y_{0} = - {\bar{y}}_{n} \\ y_{1} = - {\bar{y}}_{n} + h_{s} \\ y_{2} = (h_{s} + h_{p}) - {\bar{y}}_{n} \end{matrix}

(5)

Where ${\bar{y}}_{n}$ is distance of neutral axis from low surface of the substrate layer, and is evaluated as below (Javed and Abdelkefi, 2019):

{\bar{y}}_{n} = \frac{(h_{p} + h_{s}) E_{p} b_{p} h_{p}}{2 (E_{p} b_{p} h_{p} + E_{s} b_{s} h_{s})} + \frac{h_{s}}{2}

(6)

It has been shown in previous works that the electrical effect of a rectangular piezoelectric layer is equivalent to a concentrated mechanical torque at the end of piezoelectric layer which is shown by $M_{1 P}$ and $M_{2 P}$ in Figure 1. Therefore the equivalent mechanical moment due to the electrical effect of piezoelectric layer may be written as follow:

\begin{matrix} M_{P} = M_{P 2} δ (X - L_{p 2}) - M_{P 1} δ (X - L_{p 1}) \\ = \frac{E_{p} d_{31} {\bar{A}}_{P}}{h_{p}} V (t) (δ (X - L_{p 2}) - δ (X - L_{p 1})), \\ {\bar{A}}_{P} = b_{p} \int_{y_{1}}^{y_{2}} ydy = \frac{b_{p}}{2} (y_{2}^{2} - y_{1}^{2}) \end{matrix}

(7)

In which V(t) is the harvested voltage by piezoelectric layer, and $d_{31}$ denotes the piezoelectric constant, where axes 3 and 1 correspond to axes $Y$ and $X$ , respectively .The function of $c (X)$ may be evaluated as below (Dai et al., 2014):

\begin{matrix} c (X) = {\begin{matrix} c_{v} (H (X) - H (X - L_{s})) + (2 π S_{T} \frac{U}{D} γ ρ_{f} D^{2}) H (X - L_{s}) \to for configuration 1, \\ c_{v} (H (X) - H (X - L_{s})) + (2 π S_{T} \frac{U}{D} γ ρ_{f} D^{2}) L_{c} δ (X - L_{s}) \to for configuration 2, \end{matrix} \\ γ = \frac{c_{D}}{4 π S_{T}} \end{matrix}

(8)

In above equations, $U$ is the wind speed, and $ρ_{f}$ is the mass density of the wind flow. Moreover $γ$ is a parameter depending on the drag coefficient, $c_{D}$ , and may be approximated as constant and equal to 0.8. $S_{T}$ is the Strouhal number which is approximately constant and equal to 0.2 in range of $300 < Re < 1.5 \times 10^{5}$ . $F_{L} (X, t)$ is related to wake oscillator variable $Q (X, t)$ as follow (Dai et al., 2014; Jia et al., 2018; Zamanian and Garibaldi, 2019).

F_{L} (X, t) = {\begin{matrix} \frac{1}{4} ρ_{f} U^{2} C_{L 0} DQ (t) H (X - L_{s}) \to for configuration 1, \\ \frac{1}{4} ρ_{f} U^{2} C_{L 0} DQ (t) L_{c} δ (X - L_{s} - \frac{D}{2}) \to for configuration 2 \end{matrix}

(9)

The wake oscillator variable, $Q (t)$ , is defined as $2 C_{L} (t) / C_{L 0}$ where $C_{L} (t)$ is the time variable lift coefficient, and $C_{L 0}$ is the lift coefficient amplitude on a fixed cylindrical structure experiencing vortex shedding and is equal to 0.3 over a large range of Reynolds number. The wake oscillator variable is coupled to the structure motion as follows (Dai et al., 2014; Jia et al., 2018; Zamanian and Garibaldi, 2019).

\begin{matrix} \frac{\partial^{2} Q (t)}{\partial t^{2}} + η (2 π S_{T} \frac{U}{D}) (Q^{2} (t) - 1) \frac{dQ (t)}{dt} \\ + {(2 π S_{T} \frac{U}{D})}^{2} Q (t) = \frac{\bar{G}}{D} a_{c}, \\ a_{c} = {\begin{matrix} \frac{\partial^{2} W (X, t)}{\partial t^{2}} | \begin{matrix} X = L_{s} + \frac{L_{C}}{2} \end{matrix} \to for configuration 1, \\ \frac{\partial^{2} W (X, t)}{\partial t^{2}} | \begin{matrix} X = L_{s} + \frac{D}{2} \end{matrix} \to for configuration 2 \end{matrix} \end{matrix}

(10)

Where $η$ and $\bar{G}$ are the damping Van der Pol parameter and the coupling coefficient, and are typically equal to $0.3$ and $12$ , respectively (Facchinetti et al., 2004). The differential equation governing the output voltage of the piezoelectric layer has been derived using Gauss’s law in Abdelkefi et al. (2011). Based on the results of Gauss’s law, the electric potential V(t), between upper and lower surface of piezoelectric may be related to dynamic deflection as follow:

\begin{matrix} e_{31} h_{p} b_{p} \int_{L_{1 p}}^{L_{2 p}} \frac{\partial^{3} W}{\partial X^{2} \partial t} dX + \frac{\in_{33} b_{p} (L_{2 p} - L_{1 p})}{h_{p}} \frac{dV (t)}{dt} \\ + \frac{V (t)}{R} = 0 \end{matrix}

(11)

Where $e_{31} = E_{p} d_{31}$ , $R$ is the electrical load resistance of the piezoelectric layer, $\in_{33}$ is the electrical permittivity factor (dielectric constant). For analytical convenience, the following changes in variables are applied to the equations of motion.

w = \frac{W}{D}, x = \frac{X}{L_{s}}, τ = \frac{t}{T}, T = \sqrt{\frac{m_{s} L_{s}^{4}}{E_{s} I_{s}}}

(12)

By applying the above variables change into equation (1) the dimensionless form of motion equations will be as:

\begin{matrix} \frac{\partial^{2}}{\partial x^{2}} ({\bar{I}}_{E} (x) \frac{\partial^{2} w (x, τ)}{\partial x^{2}}) + {\bar{c}}_{v} (x) \frac{\partial w (x, τ)}{\partial τ} \\ + {\bar{c}}_{d} (x) Ω_{f} \frac{\partial w (x, τ)}{\partial τ} + \bar{m} (x) \frac{\partial^{2} w (x, τ)}{\partial τ^{2}} - \\ \frac{\partial}{\partial x} ({\bar{I}}_{m} (x) \frac{\partial^{3} w (x, τ)}{\partial τ^{2} \partial x}) = {\bar{f}}_{L} (x) Ω_{f}^{2} Q + {\bar{f}}_{p} (x) V (τ), \\ \frac{d^{2} q (τ)}{\partial τ^{2}} + η Ω_{f} (q^{2} (τ) - 1) \frac{dq (τ)}{d τ} + Ω_{f}^{2} q (x, τ) = G {\bar{a}}_{c} \end{matrix}

(13)

Where

\begin{array}{l} {\bar{I}}_{E} (x) = (H (x) - H (x - \frac{L_{1 p}}{L_{s}})) + \frac{I_{s p}}{I_{s}} (H (x - \frac{L_{1 p}}{L_{s}}) - H (x - \frac{L_{2 p}}{L_{s}})) + (H (x - \frac{L_{2 p}}{L_{s}}) - H (x - 1)) + \\ \frac{I_{c}}{I_{s}} (H (x - 1) - H (x - \frac{L_{e}}{L_{s}})), \\ \bar{m} (x) = (H (x) - H (x - \frac{L_{1 p}}{L_{s}})) + \frac{m_{s p}}{m_{s}} (H (x - \frac{L_{1 p}}{L_{s}}) - H (x - \frac{L_{2 p}}{L_{s}})) + (H (x - \frac{L_{2 p}}{L_{s}}) - H (x - 1)) + \\ \frac{m_{c}}{m_{s}} (H (x - 1) - H (x - \frac{L_{e}}{L_{s}})), \\ q (τ) = Q (t) | \begin{array}{l} t = T τ \end{array}, \\ {\bar{a}}_{c} = {\begin{cases} \frac{\partial^{2} w (x, τ)}{\partial τ^{2}} | \begin{array}{l} x = 1 + \frac{L_{C}}{2 L_{s}} \end{array} \to for configuration 1 \\ \frac{\partial^{2} w (x, τ)}{\partial τ^{2}} | \begin{array}{l} x = 1 + \frac{D}{2 L_{s}} \end{array} \to for configuration 2 \end{cases} \\ {\bar{c}}_{v} (X) = \frac{c_{v} T}{m_{s}} (H (x) - H (x - 1)), \\ {\bar{I}}_{m c} = \frac{I_{m c}}{m_{s} L_{s}^{2}} (H (x) - H (x - \frac{L_{e}}{L_{s}})) \end{array}

(14)

and

\begin{array}{l} {\bar{c}}_{d} (x) = {\begin{cases} \frac{γ ρ_{f} D^{2}}{m_{s}} (H (x - 1) - H (x - \frac{L_{e}}{L_{s}})), \to for configuration 1, \\ \frac{γ ρ_{f} D^{2}}{m_{s} L_{s}} L_{c} δ (x - 1 - \frac{D}{2 L_{s}}) \to for configuration 2 \end{cases}, \\ {\bar{f}}_{L} (x) = {\begin{cases} \frac{ρ_{f} C_{L 0} D^{2}}{16 π^{2} m_{s} S_{T}^{2}} (H (x - 1) - H (x - \frac{L_{e}}{L_{s}})), \to for configuration 1, \\ \frac{ρ_{f} C_{L 0} D^{2}}{16 π^{2} m_{s} S_{T}^{2} L_{s}} L_{c} δ (x - 1 - \frac{D}{2 L_{s}})) \to for configuration 2, \end{cases} \\ {\bar{f}}_{p} (x) = \frac{E_{p} d_{31} {\bar{A}}_{P}}{h_{p}} \frac{T^{2}}{D L_{s}^{2} m_{s}} (\frac{d}{d x} δ (x - \frac{L_{2 p}}{L_{s}}) - \frac{d}{d x} δ (x - \frac{L_{1 p}}{L_{s}})), Ω_{f} = \frac{2 π S_{T} U T}{D} \end{array}

(15)

By applying the dimensionless quantities of equation (12) into equation (11) it can be written as:

\frac{e_{31} h_{p} b_{p} D}{L_{s} T} \int_{L_{1 p} / L_{s}}^{L_{2 p} / L_{s}} \frac{\partial^{3} w}{\partial x^{2} \partial τ} d x + \frac{(L_{2 p} - L_{1 p}) \in_{33} b_{p}}{h_{p} T} \frac{d V (τ)}{d τ} + \frac{V (τ)}{R} = 0

(16)

3. Discretization using Galerkin procedure

Now, the free vibration response of the system is assumed as

w (x, τ) = φ (x) e^{i ω τ}

(17)

Where $ω$ is the natural frequency of free vibration of the structure, and $φ (x)$ is the mode shape of the system correspond to natural frequency $ω$ . By substituting equation (17) into equation (13), the differential equation governing the mode shape of the system is written as follow

\frac{d^{2}}{d x^{2}} ({\bar{I}}_{E} (x) \frac{d^{2} φ (x)}{d x^{2}}) - \bar{m} (x) ω^{2} φ (x) = 0

(18)

Galerkin method is used to solve the above equation. Firstly, the free vibration mode shape of the beam coated by piezoelectric layer and attached to the cylinder at the end are obtained as comparison function. For configuration 2, the comparison function are obtained assuming the cylinder as a cubic cylinder as shown in Figure 1(c). If one assumes that $ψ_{sp} (x)$ , $ψ_{s 1} (x)$ , $ψ_{s 2} (x)$ are the functions of mode shape at different part of the assumed system, then the differential equation may be rewritten as below, where $ω_{c}$ is the corresponding natural frequency for the comparison function.

\begin{matrix} - ω_{c}^{2} ψ_{s 1} (x) + \frac{d^{4} ψ_{s 1} (x)}{d x^{4}} = 0, 0 < x < \frac{L_{1 p}}{L_{s}}, \\ - ω_{c}^{2} (\frac{m_{sp}}{m_{s}}) ψ_{sp} (x) + (\frac{I_{sp}}{I_{s}}) \frac{d^{4} ψ_{sp} (x)}{d x^{4}} = 0 \frac{L_{1 p}}{L_{s}} < x < \frac{L_{2 P}}{L_{s}}, \\ - ω_{c}^{2} ψ_{s 2} (x) + \frac{d^{4} ψ_{s 2} (x)}{d x^{4}} = 0, \frac{L_{2 P}}{L_{s}} < x < 1, \\ - ω_{c}^{2} (\frac{m_{cc}}{m_{s}}) ψ_{c} (x) + (\frac{I_{cc}}{I_{s}}) \frac{d^{4} ψ_{c} (x)}{d x^{4}} = 01 < x < \frac{L_{e}}{L_{s}}, \end{matrix}

(19)

Where $m_{cc}$ and $I_{cc}$ are the amount of mass per unit length and second moment of area for the cylinder for obtaining the comparison functions. These values in configuration 1 are identical to $m_{c}$ and $I_{c}$ , and are as below in configuration 2 considering the cylinder as a cubic cylinder as shown in Figure 1(b).

\begin{matrix} m_{cc} = ρ_{c} π {(\frac{D}{2})}^{2} \to for configuration 1, \\ m_{cc} = ρ_{c} L_{c} D \to for configuration 2, \\ I_{cc} = E_{c} \frac{π}{4} {(\frac{D}{2})}^{4} \to for configuration 1, \\ I_{cc} = E_{c} \frac{1}{12} L_{c} D^{3} \to for configuration 2 \end{matrix}

(20)

The following boundary conditions must be satisfied for differential equation (19)

\begin{array}{l} ψ_{s 1} |_{x = 0} = 0, \frac{d ψ_{s 1}}{d x} |_{x = 0} = 0, \\ ψ_{s 1} |_{x = \frac{L_{1 p}}{L_{s}}} = ψ_{s p} |_{x = \frac{L_{1 p}}{L_{s}}}, \frac{d ψ_{s p}}{d x} |_{x = \frac{L_{1 p}}{L_{s}}} = \frac{d ψ_{s 1}}{d x} |_{x = \frac{L_{1 p}}{L_{s}}}, \\ I_{s} \frac{d^{2} ψ_{s 1}}{d x^{2}} |_{x = \frac{L_{1 p}}{L_{s}}} = I_{s} \frac{d^{2} ψ_{s p}}{d x^{2}} |_{x = \frac{L_{1 p}}{L_{s}}}, I_{s} \frac{d^{3} ψ_{s 1}}{d x^{3}} |_{x = \frac{L_{1 p}}{L_{s}}} = I_{s p} \frac{d^{3} ψ_{s p}}{d x^{3}} |_{x = \frac{L_{1 p}}{L_{s}}}, \\ ψ_{s 2} |_{x = \frac{L_{2 p}}{L_{s}}} = ψ_{s p} |_{x = \frac{L_{2 p}}{L_{s}}}, \frac{d ψ_{s p}}{d x} |_{x = \frac{L_{2 p}}{L_{s}}} = \frac{d ψ_{s 2}}{d x} |_{x = \frac{L_{2 p}}{L_{s}}}, \\ I_{s} \frac{d^{2} ψ_{s 2}}{d x^{2}} |_{x = \frac{L_{2 p}}{L_{s}}} = I_{s p} \frac{d^{2} ψ_{s p}}{d x^{2}} |_{x = \frac{L_{2 p}}{L_{s}}}, I_{s} \frac{d^{3} ψ_{s 2}}{d x^{3}} |_{x = \frac{L_{2 p}}{L_{s}}} = I_{s p} \frac{d^{3} ψ_{s p}}{d x^{3}} |_{x = \frac{L_{2 p}}{L_{s}}}, \\ ψ_{s 2} |_{x = 1} = ψ_{c} |_{x = 1}, \frac{d ψ_{s 2}}{d x} |_{x = 1} = \frac{d ψ_{c}}{d x} |_{x = 1}, \\ I_{s} \frac{d^{2} ψ_{s 2}}{d x^{2}} |_{x = 1} = I_{c c} \frac{d^{2} ψ_{c}}{d x^{2}} |_{x = 1}, I_{s} \frac{d^{3} ψ_{s 2}}{d x^{3}} |_{x = 1} = I_{c c} \frac{d^{3} ψ_{c}}{d x^{3}} |_{x = 1}, \\ \frac{d^{3} ψ_{c}}{d x^{3}} |_{x = \frac{L_{e}}{L_{s}}} = 0, \frac{d^{2} ψ_{c}}{d x^{2}} |_{x = \frac{L_{e}}{L_{s}}} = 0, \end{array}

(21)

The general solution of equation (19) may be expressed as follows

\begin{matrix} ψ_{s 1} (x) = C_{1} \sin (\sqrt{β_{1}} x) + C_{2} \cos (\sqrt{β_{1}} x) \\ + C_{3} \cosh (\sqrt{β_{1}} x) + C_{4} \sinh (\sqrt{β_{1}} x), \\ ψ_{sp} (x) = C_{5} \sin (\sqrt{β_{2}} x) + C_{6} \cos (\sqrt{β_{2}} x) \\ + C_{7} \cosh (\sqrt{β_{2}} x) + C_{8} \sinh (\sqrt{β_{2}} x), \\ ψ_{s 2} (x) = C_{9} \sin (\sqrt{β_{1}} x) + C_{10} \cos (\sqrt{β_{1}} x) \\ + C_{11} \cosh (\sqrt{β_{1}} x) + C_{12} \sinh (\sqrt{β_{1}} x), \\ ψ_{c} (x) = C_{13} \sin (\sqrt{β_{3}} x) + C_{14} \cos (\sqrt{β_{3}} x) \\ + C_{15} \cosh (\sqrt{β_{3}} x) + C_{16} \sinh (\sqrt{β_{3}} x) \\ β_{1}^{2} = ω_{c}^{2}, β_{2}^{2} = \frac{m_{sp} I_{s} ω_{c}^{2}}{m_{s} I_{sp}}, β_{3}^{2} = \frac{m_{cc} I_{s} ω_{c}^{2}}{m_{s} I_{cc}} \end{matrix}

(22)

By applying the boundary conditions described in (21) into equation (22), and setting the determinant of coefficient matrix equal to zero and solving the characteristic equation, the natural frequency will be obtained. Therefore the mode shape correspond to the natural frequency $ω_{c}$ may be written as below:

\begin{array}{l} ψ (x) = (H (x) - H (x - \frac{L_{1 p}}{L_{s}})) ψ_{s 1} (x) + (H (x - \frac{L_{1 p}}{L_{s}}) - H (x - \frac{L_{2 p}}{L_{s}})) ψ_{s p} (x) + \\ (H (x - \frac{L_{2 p}}{L_{s}}) - H (x - 1)) ψ_{s 2} (x) + (H (x - 1) - H (x - \frac{L_{e}}{L_{s}})) ψ_{c} (x), \end{array}

(23)

Now the solution of equation (18) is assumed as

φ (x) = \sum_{i = 1}^{N} {\bar{C}}_{i} ψ_{i} (x)

(24)

Where ${\bar{C}}_{i}, i = 1 . . n$ are the constant coefficients that must be obtained after applying Galerkin procedure, and $ψ_{i} (x)$ is the ith comparison function correspond to ith characteristic value. So, substituting equation (24) into equation (18), and multiplying the outcome by $ψ_{j} (x)$ and integrating it at the interval $x = 0 . . \frac{L_{e}}{L_{s}}$ , it results that

\begin{matrix} \int_{0}^{\frac{L_{e}}{L_{s}}} \sum_{i = 1}^{n} {\bar{C}}_{i} \frac{d^{2}}{d x^{2}} ({\bar{I}}_{E} (x) \frac{d^{2} ψ_{i} (x)}{d x^{2}}) ψ_{j} (x) dx \\ - ω^{2} \int_{0}^{\frac{L_{e}}{L_{s}}} \sum_{i = 1}^{n} {\bar{C}}_{i} \bar{m} (x) ψ_{i} (x) ψ_{j} (x) dx \end{matrix}

(25)

The natural frequency of system is obtained setting the determinants of constant coefficients ${\bar{C}}_{i}$ equal to zero in equation (25). After that substituting the natural frequency in the algebraic system, the coefficients of ${\bar{C}}_{i}$ are obtained. Now, the structure oscillation is assumed as $w (x, t) = φ (x) p (t)$ where $p (t)$ is time coordinate function correspond to the structure. By substituting this assumptions into equation (13), multiplying the outcome by $φ (x)$ and integrating on the whole length of the structure, the differential equation governing the system oscillations will be as:

\begin{array}{l} \frac{d^{2} p (τ)}{d τ^{2}} + 2 ξ ω_{n} \frac{d p (τ)}{d τ} + c_{d e q} Ω_{f} \frac{d p (τ)}{d τ} + ω_{n}^{2} p (τ) = f_{L e q} Ω_{f}^{2} q (τ) + f_{p e q} V (τ), \\ \frac{d^{2} q (τ)}{d τ^{2}} + η Ω_{f} (q^{2} (τ) - 1) \frac{d q (τ)}{d τ} + Ω_{f}^{2} q (τ) = G \frac{d^{2} p (τ)}{d τ^{2}}, \\ k_{e q} = (\int_{0}^{\frac{L_{e}}{L_{s}}} φ (x) \frac{d^{2}}{d x^{2}} (\bar{I} (x) \frac{d^{2} φ (x)}{d x^{2}}) d x = \int_{0}^{\frac{L_{e}}{L_{s}}} \bar{I} (x) \frac{d^{2} φ (x)}{d x^{2}} \frac{d^{2} φ (x)}{d x^{2}} d x, \\ m_{e q} = \int_{0}^{\frac{L_{e}}{L_{s}}} (\bar{m} (x) φ^{2} (x) + \frac{d}{d x} ({\bar{I}}_{m} (x) \frac{d φ (x)}{d x})) d x, ω_{n} = \sqrt{\frac{k_{e q}}{m_{e q}}}, \\ ξ = \frac{\int_{0}^{\frac{L_{e}}{L_{s}}} {\bar{c}}_{v} (x) φ^{2} (x)}{2 m_{e q} ω_{n}}, c_{d e q} = \frac{\int_{0}^{\frac{L_{e}}{L_{s}}} {\bar{c}}_{d} (x) φ^{2} (x) d x}{m_{e q}}, f_{L e q} = \frac{\int_{0}^{\frac{L_{e}}{L_{s}}} {\bar{f}}_{L} (x) φ (x) d x}{m_{e q}}, \\ f_{p e q} = \frac{\int_{0}^{\frac{L_{e}}{L_{s}}} {\bar{f}}_{p} (x) φ (x) d x}{m_{e q}}, G = \bar{G} {\begin{cases} φ {(x)}_{| x = 1 + \frac{L_{c}}{2 L_{s}}} \to for configuration 1, \\ φ {(x)}_{| x = 1 + \frac{D}{2 L_{s}}} \to for configuration 2 \end{cases} \end{array}

(26)

Substituting the assumption of $w (x, t) = φ (x) p (t)$ into equation (16), it may be rewritten as follows:

\begin{matrix} \frac{e_{31} h_{p} b_{p} D}{L_{s} T} \frac{dp}{d τ} \int_{L_{1 p} / L_{s}}^{L_{2 p} / L_{s}} \frac{d^{2} φ (x)}{d x^{2}} dx + \frac{\in_{33} b_{p} (L_{2 p} - L_{1 p})}{h_{p} T} \frac{dV (τ)}{d τ} \\ + \frac{V (τ)}{R} = 0 \end{matrix}

(27)

By renaming the expression of equation (27), one obtains

χ \frac{dp}{d τ} + C_{p} \frac{dV (τ)}{d τ} + \frac{V (τ)}{R} = 0

(28)

In which

\begin{matrix} χ = \frac{e_{31} b_{p} h_{p} D}{L_{s} T} (\frac{d φ (x)}{dx} | \begin{matrix} x = \frac{L_{2 p}}{L_{s}} \end{matrix} - \frac{d φ (x)}{dx} | \begin{matrix} x = \frac{L_{1 p}}{L_{s}} \end{matrix}), \\ C_{p} = \frac{ϵ_{33} b_{p}}{h_{p} T} (L_{2 p} - L_{1 p}) \end{matrix}

(29)

4. Discretization using assumed mode method

In this procedure which has been used in previous works, firstly the free vibration equations of motion are obtained using Hamilton principle assuming the joint cross section of the beam and cylinder as boundary of system. Then free vibration mode shapes of system are obtained. After that the discretized equations of motion governing on the vortex induced vibration are obtained using assumed mode method by applying the Lagrange equation. For applying the Hamilton principle, the variation of kinematic energy, $Π_{k}$ , and potential energy may be written as below:

\begin{array}{l} Π_{k} = \frac{1}{2} \int_{0}^{L_{s}} m (X) \frac{\partial^{2} W}{\partial t^{2}} d X + (\frac{1}{2} M_{c} {(\frac{\partial W}{\partial t} | \begin{array}{l} X = L_{s} \end{array} + r_{c} \frac{\partial^{2} W}{\partial X \partial t})}^{2} + \frac{1}{2} I_{M c} {(\frac{\partial^{2} W}{\partial X \partial t})}^{2}) | \begin{array}{l} X = L_{s} \end{array}, \\ Π_{ε} = \frac{1}{2} \int_{0}^{L_{s}} I_{E} (X) (\frac{\partial^{2} W}{\partial X^{2}})^{2} d X - \frac{1}{2} \int_{0}^{L_{s}} (H (X - L_{1 p}) - H (X - L_{2 p})) e_{31} b_{p} V (t) \frac{\partial^{2} W}{\partial X^{2}} d X \\ r_{c} = {\begin{cases} \frac{L_{c}}{2} \to for configuration 1 \\ \frac{D}{2} \to for configuration 2 \end{cases}, I_{M c} = {\begin{cases} \frac{1}{12} M_{c} L_{c}^{2} \to for configuration 1 \\ \frac{1}{2} M_{c} {(\frac{D}{2})}^{2} \to for configuration 2 \end{cases} \end{array}

(30)

Where $M_{c}$ is the total mass of the cylinder, and $I_{Mc}$ is the total mass moment of inertia around mass center for the cylinder .The potential energy, and variation of external work due to non-conservative external force may be written as below:

\begin{matrix} δ Π_{F} = c_{v} \int_{0}^{L_{c}} δ \overset{\cdot}{W} dX + \frac{1}{4} ρ_{f} U^{2} C_{L 0} DQ (t) \int_{L_{s}}^{L_{e}} δ W_{c} dX \\ - (2 π S_{T} \frac{U}{D} γ ρ_{f} D^{2}) \int_{L_{s}}^{L_{c}} {\overset{\cdot}{W}}_{c} δ W_{c} dX \end{matrix}

(31)

It should be noted that the amount of $I_{E} (X)$ is evaluated considering equation (2) with this difference that here it does not include the term of $I_{c} (X)$ . In addition $W_{c}$ describes the function of transverse displacement along the length of the cylinder which may be described as below according that it has been assumed as a rigid body:

\begin{matrix} W_{c} = W | X = L_{s} + r_{c} \frac{\partial W}{\partial X} | \begin{matrix} X = L_{s} \end{matrix} \\ r_{c} = {\begin{matrix} \frac{L_{c}}{2} \to For configuration 1, \\ \frac{D}{2} \to For configuration 2 \end{matrix} \end{matrix}

(32)

By applying the Hamilton principle, and setting $V (t) = 0, δ Π_{F} = 0$ , the free motion equations will be as:

\begin{array}{l} \int_{t_{1}}^{t_{2}} δ (Π_{k} - Π_{ε}) + δ Π_{F}) d t = 0, \\ \frac{\partial^{2}}{\partial x^{2}} (I_{E} (X) \frac{\partial^{4} W (X, t)}{\partial X^{4}}) + m (X) \frac{\partial^{2} W (X, t)}{\partial t^{2}} = 0 \\ W | \begin{array}{l} X = 0 \end{array} = 0, \frac{\partial W}{\partial X} | \begin{array}{l} X = 0 \end{array} = 0, \\ I_{E} (X) \frac{\partial^{3} W}{\partial X^{3}} - M_{c} (\frac{\partial^{2} W}{\partial t^{2}} + r_{c} \frac{\partial^{3} W}{\partial t^{2} \partial X}) |_{X = L_{s}} = 0 \\ I_{E} (X) \frac{\partial^{2} W}{\partial X^{2}} + I_{M c} (\frac{\partial^{3} W}{\partial X \partial t^{2}}) + M_{c} r_{c} (\frac{\partial^{2} W}{\partial t^{2}} + r_{c} \frac{\partial^{3} W}{\partial t^{2} \partial X}) |_{X = L_{s}} = 0 \end{array}

(33)

The dimensionless form of the equation (33) will be as below after substituting the variables of equation (12) into equation (33):

\begin{matrix} \frac{\partial^{2}}{\partial x^{2}} ({\bar{I}}_{E} (x) \frac{\partial^{2} w}{\partial x^{2}}) + \bar{m} (x) \frac{\partial^{2} w}{\partial τ^{2}} = 0, \\ w |_{x = 0} = 0, \frac{\partial w}{\partial x} |_{x = 0} = 0, \\ ({\bar{I}}_{E} (x) \frac{\partial^{3} w}{\partial x^{3}} - \frac{L_{s}^{3} M_{c}}{T^{2}} (\frac{\partial^{2} w}{\partial τ^{2}} + \frac{r_{c}}{L_{s}} \frac{\partial^{3} w}{\partial τ^{2} \partial x})) |_{x = 1} = 0 \\ ({\bar{I}}_{E} (x) \frac{\partial^{2} w}{\partial x^{2}} + \frac{L_{s}^{2} r_{c} M_{c}}{T^{2}} (\frac{\partial^{2} w}{\partial τ^{2}} + \frac{r_{c}}{L_{s}} \frac{\partial^{3} w}{\partial τ^{2} \partial x}) + \frac{L_{s} I_{Mc}}{T^{2}} \frac{\partial^{3} w}{\partial τ^{2} \partial x}) |_{x = 1} = 0 \end{matrix}

(34)

Assumed mode method is used to discretize the equation of motion. Firstly, the exact mode shape of free vibration of system that is, the mode shape considering the non-uniformity of the cross section, is extracted to be used as comparison function. Therefore, the free vibration response of system is assumed as:

w (x, t) = {\begin{matrix} ψ_{s 1} (x) e^{i ω t}, 0 \leq x \leq \frac{L_{1 p}}{L_{s}}, \\ ψ_{sp} (x) e^{i ω t}, \frac{L_{1 p}}{L_{s}} \leq x \leq \frac{L_{2 p}}{L_{s}}, \\ ψ_{s 2} (x) e^{i ω t}, \frac{L_{2 p}}{L_{s}} \leq x \leq 1 \end{matrix}

(35)

Where as mentioned before $ψ_{s 1} (x)$ , $ψ_{sp} (x)$ , $ψ_{s 2} (x)$ are the functions of mode shape at different parts of the system which here are obtained considering the joint cross section of the beam and cylinder as boundary of system. Substituting equation (35) into equation (34), it results in

\begin{matrix} - ω^{2} ψ_{s 1} (x) + \frac{d^{4} ψ_{s 1} (x)}{d x^{4}} = 0, 0 < x < \frac{L_{1 p}}{L_{s}}, \\ - ω^{2} (\frac{m_{sp}}{m_{s}}) ψ_{sp} (x) + (\frac{I_{sp}}{I_{s}}) \frac{d^{4} ψ_{sp} (x)}{d x^{4}} = 0 \frac{L_{1 p}}{L_{s}} < x < \frac{L_{2 P}}{L_{s}}, \\ - ω^{2} ψ_{s 2} (x) + \frac{d^{4} ψ_{s 2} (x)}{d x^{4}} = 0, \frac{L_{2 P}}{L_{s}} < x < 1 \end{matrix}

(36)

The general solution of equations with the given boundary conditions could be expressed as follows:

\begin{matrix} ψ_{s 1} (x) = C_{1} \sin (\sqrt{β_{1}} x) + C_{2} \cos (\sqrt{β_{1}} x) \\ + C_{3} \cosh (\sqrt{β_{1}} x) + C_{4} \sinh (\sqrt{β_{1}} x), \\ ψ_{sp} (x) = C_{5} \sin (\sqrt{β_{2}} x) + C_{6} \cos (\sqrt{β_{2}} x) \\ + C_{7} \cosh (\sqrt{β_{2}} x) + C_{8} \sinh (\sqrt{β_{2}} x), \\ ψ_{s 2} (x) = C_{9} \sin (\sqrt{β_{1}} x) + C_{10} \cos (\sqrt{β_{1}} x) \\ + C_{11} \cosh (\sqrt{β_{1}} x) + C_{12} \sinh (\sqrt{β_{1}} x), \\ β_{1}^{2} = ω^{2}, β_{2}^{2} = \frac{m_{sp} I_{s} ω^{2}}{m_{b} I_{sp}} \end{matrix}

(37)

where $ω$ is obtained by inserting the following boundary condition equations into equation (37), setting the determinant of coefficients matrix to zero, and solving the characteristic equation.

\begin{array}{l} ψ_{s 1} |_{x = 0} = 0, \frac{d ψ_{s 1}}{d x} |_{x = 0} = 0, \\ ψ_{s 1} |_{x = \frac{L_{1 p}}{L_{s}}} = ψ_{s p} |_{x = \frac{L_{1 p}}{L_{s}}}, \frac{d ψ_{s p}}{d x} |_{x = \frac{L_{1 p}}{L_{s}}} = \frac{d ψ_{s 1}}{d x} |_{x = \frac{L_{1 p}}{L_{s}}}, \\ I_{s} \frac{d^{2} ψ_{s 1}}{d x^{2}} |_{x = \frac{L_{1 p}}{L_{s}}} = I_{s p} \frac{d^{2} ψ_{s p}}{d x^{2}} |_{x = \frac{L_{1 p}}{L_{s}}}, I_{s} \frac{d^{3} ψ_{s 1}}{d x^{3}} |_{x = \frac{L_{1 p}}{L_{s}}} = I_{s p} \frac{d^{3} ψ_{s p}}{d x^{3}} |_{x = \frac{L_{1 p}}{L_{s}}}, \\ ψ_{s 2} |_{x = \frac{L_{2 p}}{L_{s}}} = ψ_{s p} |_{x = \frac{L_{2 p}}{L_{s}}}, \frac{d ψ_{s p}}{d x} |_{x = \frac{L_{2 p}}{L_{s}}} = \frac{d ψ_{s 2}}{d x} |_{x = \frac{L_{2 p}}{L_{s}}}, \\ I_{s} \frac{d^{2} ψ_{s 2}}{d x^{2}} |_{x = \frac{L_{2 p}}{L_{s}}} = I_{s p} \frac{d^{2} ψ_{s p}}{d x^{2}} |_{x = \frac{L_{2 p}}{L_{s}}}, I_{s} \frac{d^{3} ψ_{s 2}}{d x^{3}} |_{x = \frac{L_{2 p}}{L_{s}}} = I_{s p} \frac{d^{3} ψ_{s p}}{d x^{3}} |_{x = \frac{L_{2 p}}{L_{s}}} \end{array}

(38)

and the following boundary condition must be added:

\begin{matrix} ({\bar{I}}_{E} (x) \frac{d^{3} ψ_{s 2}}{d x^{3}} - \frac{L_{s}^{3} M_{c}}{T^{2}} (- ω^{2} ψ_{s 2} - ω^{2} \frac{r_{c}}{L_{s}} \frac{d ψ_{s 2}}{dx})) |_{x = 1} = 0, \\ ({\bar{I}}_{E} (x) \frac{\partial^{2} ψ_{s 2}}{\partial x^{2}} + \frac{L_{s}^{2} r_{c} M_{c}}{T^{2}} (- ω^{2} ψ_{s 2} - ω^{2} \frac{r_{c}}{L_{s}} \frac{d ψ_{s 2}}{dx}) \\ - \frac{L_{s} I_{Mc}}{T^{2}} ω^{2} \frac{d ψ_{s 2}}{dx}) |_{x = 1} = 0 \end{matrix}

(39)

Therefore, system mode shape will be:

\begin{array}{l} φ (x) = (H (x) - H (x - \frac{L_{1 p}}{L_{s}})) ψ_{s 1} (x) + (H (x - \frac{L_{1 p}}{L_{s}}) - H (x - \frac{L_{2 p}}{L_{s}})) ψ_{s p} (x) + \\ (H (x - \frac{L_{2 p}}{L_{s}}) - H (x - 1)) ψ_{s 2} (x) \end{array}

(40)

Assuming the response of system under vortex induced vibration as $w (x, τ) = p (t) φ (x)$ , the dimensionless form of the kinematic energy will be as:

\begin{array}{l} Π_{k} = \frac{1}{2} m_{e q} {(\frac{d p (τ)}{d τ})}^{2}, \\ m_{e q} = \frac{D^{2} L_{s}}{2 T^{2}} \int_{0}^{1} (m (X) |_{X = L_{s} x} φ^{2} (x)) d x + \\ (\frac{1}{2} M_{c} \frac{D^{2}}{T^{2}} {(φ (x) + \frac{r_{c}}{L_{s}} \frac{d φ (x)}{d x})}^{2} |_{\begin{array}{l} x = 1 \end{array}} + \frac{1}{2} I_{M c} \frac{D^{2}}{L_{s}^{2} T^{2}} {(\frac{d φ (x)}{d x})}^{2} |_{\begin{array}{l} x = 1 \end{array}}) \end{array}

(41)

Moreover, the potential energy of the system may be written as below:

\begin{array}{l} Π_{ε} = \frac{1}{2} k_{e q 1} p^{2} (τ) + k_{e q 2} p (τ) V (τ), \\ k_{e q 1} = \frac{D^{2}}{L_{s}^{3}} \int_{0}^{1} I_{E} (X) |_{X = L_{s} x} {(\frac{d^{2} φ (x)}{d x^{2}})}^{2} d x, \\ k_{e q 2} = - \int_{0}^{1} e_{31} b_{p} \frac{D}{L_{s}} \frac{d^{2} φ (x)}{d x^{2}} (H (x - \frac{L_{1 p}}{L_{s}}) - H (x - \frac{L_{2 p}}{L_{s}})) d x \end{array}

(42)

and the variations of the external forces including lift force, drag damping force and viscosity damping force will be as:

\begin{matrix} δ Π_{F} = {\bar{f}}_{v} \frac{dp}{d τ} δ p (τ) + {\bar{f}}_{L} Ω_{f}^{2} q (t) δ p (τ) + {\bar{f}}_{d} Ω_{f} \frac{dp}{d τ} δ p (τ), \\ {\bar{f}}_{v} = \frac{{\bar{c}}_{v} D L_{s}}{T} \int_{0}^{1} φ^{2} dx, \\ {\bar{f}}_{L} = \frac{1}{4} ρ_{f} \frac{D^{4}}{{(2 S_{T} π T)}^{2}} C_{L 0} L_{s} \int_{1}^{L_{e} / L_{s}} φ_{c} dx \to for configuration 1, \\ {\bar{f}}_{L} = \frac{1}{4} ρ_{f} \frac{D^{4}}{{(2 S_{T} π T)}^{2}} C_{L 0} L_{c} φ_{c} \to for configuration 2, \\ {\bar{f}}_{d} = (\frac{γ ρ_{f} D^{4} L_{s}}{T^{2}}) \int_{1}^{L_{e} / L_{s}} φ_{c}^{2} dx \to for configuration 1, \\ {\bar{f}}_{d} = (\frac{γ ρ_{f} D^{4}}{T^{2}}) L_{c} φ_{c}^{2} \to for configuration 2 \end{matrix}

(43)

In which $φ_{c}$ is the function of mode shape at the cylinder part, and may be evaluated as below considering that it is assumed as a rigid body:

\begin{matrix} φ_{c} = φ |_{x = 1} + (x - 1) \frac{\partial φ}{\partial x} |_{x = 1}, \\ φ_{c} = φ |_{x = 1} + \frac{D}{2 L_{s}} \frac{\partial φ}{\partial x} |_{x = 1} \end{matrix}

(44)

The equation of motion governing on $p (t)$ may be obtained using Lagrange equation as:

\frac{d}{d τ} (\frac{\partial}{\partial \dot{p}} (Π_{k} - Π_{ε})) - \frac{\partial}{\partial p} (Π_{k} - Π_{ε})) = \frac{δ Π_{F}}{δ p}

(45)

It results in:

\begin{matrix} \frac{d^{2} p (τ)}{d τ^{2}} + 2 ξ ω_{n} \frac{dp (τ)}{d τ} + c_{deq} Ω_{f} \frac{dp (τ)}{d τ} + ω_{n}^{2} p (τ) \\ = f_{Leq} Ω^{2} q (τ) + f_{peq} V (τ), \\ ω_{n} = \sqrt{\frac{k_{eq 1}}{m_{eq}}}, c_{deq} = \frac{{\bar{f}}_{d}}{m_{eq}}, \\ f_{Leq} = \frac{{\bar{f}}_{L}}{m_{eq}}, f_{peq} = \frac{k_{eq 2}}{m_{eq}}, \\ ξ = \frac{1}{2 Q_{f}} = \frac{{\vec{f}}_{v}}{2 m_{eq} ω_{n}}, \end{matrix}

(46)

The Van der Pol equation governing to the lift coefficient may be written as below:

\begin{matrix} \frac{\partial^{2} Q (t)}{\partial t^{2}} + η (2 π S_{T} \frac{U}{D}) (Q^{2} (t) - 1) \frac{dQ (t)}{dt} \\ + {(2 π S_{T} \frac{U}{D})}^{2} Q (t) = \frac{G}{D} a_{c}, \\ a_{c} = {\begin{matrix} \frac{\partial^{2} W (X, t)}{\partial t^{2}} | \begin{matrix} X = L_{s} \end{matrix} + \frac{L_{c}}{2} \frac{\partial^{3} W (X, t)}{\partial t^{2} \partial X} \to for configuration 1 \\ \frac{\partial^{2} W (X, t)}{\partial t^{2}} | \begin{matrix} X = L_{s} \end{matrix} + \frac{D}{2} \frac{\partial^{3} W (X, t)}{\partial t^{2} \partial X} \to for configuration 2 \end{matrix} \end{matrix}

(47)

The dimensionless for of the equation will be as below:

\begin{matrix} \frac{d^{2} q (τ)}{d τ^{2}} + η Ω_{f} (q^{2} (τ) - 1) \frac{dq (τ)}{d τ} + Ω_{f}^{2} q (τ) = G \frac{d^{2} p (τ)}{d τ^{2}}, \\ G = \bar{G} {\begin{matrix} φ | \begin{matrix} x = 1 \end{matrix} + \frac{L_{c}}{2 L_{s}} \frac{d φ}{dx} | \begin{matrix} x = 1 \end{matrix} \to for configuration 1, \\ φ | \begin{matrix} x = 1 \end{matrix} + \frac{D}{2 L_{s}} \frac{d φ}{dx} | \begin{matrix} x = 1 \end{matrix} \to for configuration 2 \end{matrix} \end{matrix}

(48)

And finally the differential equation governing the output voltage will be identical to equation (28).

5. Perturbation solution

To strike a balance between the terms of equations (26) and (28), the bookkeeping parameter $ε = 1$ which represent the order of terms is used as below

\begin{matrix} \frac{d^{2} p (τ)}{d τ^{2}} + ε 2 ξ ω_{n} \frac{dp (τ)}{d τ} + ε c_{deq} Ω_{f} \frac{dp (τ)}{d τ} \\ + ω_{n}^{2} p (τ) = ε f_{Leq} Ω_{f}^{2} q (τ) + ε f_{peq} V (τ), \\ \frac{d^{2} q (τ)}{d τ^{2}} + ε η Ω_{f} (q^{2} (τ) - 1) \frac{dq (τ)}{d τ} + Ω_{f}^{2} q (τ) \\ = ε G \frac{d^{2} p (τ)}{d τ^{2}}, χ \frac{dp (τ)}{d τ} + C_{p} \frac{dV (τ)}{d τ} + \frac{V (τ)}{R} = 0 \end{matrix}

(49)

Now the response of equation (49) is assumed as follows using multiple scales method of perturbation theory:

\begin{matrix} P (τ) = p_{0} (T_{0}, T_{1}) + ε p_{1} (T_{0}, T_{1}) . \\ q (τ) = q_{0} (T_{0}, T_{1}) + ε q_{1} (T_{0}, T_{1}) . \\ V (τ) = V_{0} (T_{0}, T_{1}) + ε V_{1} (T_{0}, T_{1}) \end{matrix}

(50)

Where $T_{0} = τ$ and $T_{1} = ε τ$ are time scales. By substituting equation (50) into equation (49) and separating terms with equal order of $ε$ , the following relations are obtained:

Order $(1)$ :

\begin{matrix} \frac{d^{2}}{d {T_{0}}^{2}} p_{0} (T_{0}, T_{1}) + ω_{n}^{2} p_{0} (T_{0}, T_{1}) = 0, \\ \frac{d^{2}}{d {T_{0}}^{2}} q_{0} (T_{0}, T_{1}) + Ω_{f}^{2} q_{0} (T_{0}, T_{1}) = 0, \\ C_{p} \frac{\partial}{\partial T_{0}} V_{0} (T_{0}, T_{1}) + \frac{1}{R} V_{0} (T_{0}, T_{1}) = - χ \frac{\partial}{\partial T_{0}} p_{0} (T_{0}, T_{1}) \end{matrix}

(51)

Order $(ε)$ :

\begin{matrix} \frac{\partial^{2}}{\partial {T_{0}}^{2}} p_{1} (T_{0}, T_{1}) + ω_{n}^{2} p_{1} (T_{0}, T_{1}) = - 2 \frac{\partial^{2} p_{0} (T_{0}, T_{1})}{\partial T_{0} \partial T_{1}} \\ - c_{deq} Ω_{f} \frac{\partial p_{0} (T_{0}, T_{1})}{\partial T_{0}} - ε 2 ξ ω_{n} \frac{\partial p_{0} (T_{0}, T_{1})}{\partial T_{0}} \\ + f_{Leq} Ω_{f}^{2} q_{0} (T_{0}, T_{1}) + f_{peq} V_{0} (T_{0}, T_{1}), \\ \frac{\partial^{2}}{\partial {T_{0}}^{2}} q_{0} (T_{0}, T_{1}) + Ω_{f}^{2} q_{0} (T_{0}, T_{1}) = - 2 \frac{\partial^{2} q_{0} (T_{0}, T_{1})}{\partial T_{0} \partial T_{1}} \\ - η Ω_{f} (q_{0}^{2} (T_{0}, T_{1}) - 1) \frac{\partial q_{0} (T_{0}, T_{1})}{\partial T_{0}} + G \frac{\partial^{2} p_{0} (T_{0}, T_{1})}{\partial {T_{0}}^{2}} \\ C_{p} \frac{\partial}{\partial T_{0}} V_{1} (T_{0}, T_{1}) + \frac{1}{R} V_{1} (T_{0}, T_{1}) = - χ \frac{\partial}{\partial T_{1}} p_{1} (T_{0}, T_{1}) \end{matrix}

(52)

The solution of equation (51) will be as follows:

\begin{matrix} p_{0} (T_{0}, T_{1}) = A (T_{1}) e^{i ω_{n} T_{0}} + \bar{A} (T_{1}) e^{- i ω_{n} T_{0}} \\ q_{0} (T_{0}, T_{1}) = B (T_{1}) e^{i Ω_{f} T_{0}} + \bar{B} (T_{1}) e^{- i Ω_{f} T_{0}}, \\ V_{0} = \frac{- i ω_{n} χ \bar{A} (T_{1}) R e^{- i ω_{n} T_{0}}}{iR C_{p} ω_{n} - 1} - \frac{i ω_{n} χ A (T_{1}) R e^{i ω_{n} T_{0}}}{iR C_{p} ω_{n} + 1} + \bar{K} e^{- \frac{T_{0}}{R C_{p}}} \end{matrix}

(53)

Where $A (T_{1})$ and $B (T_{1})$ are the complex coefficient and $\bar{A} (T_{1})$ and $\bar{B} (T_{1})$ are their complex conjugates. It is clear that the term of $\bar{K} e^{- \frac{T_{0}}{R C_{p}}}$ will be equal to zero in steady state response. It means that obtaining the coefficient of $A (T_{1})$ and $B (T_{1})$ . is sufficient for steady state response of system. These coefficients are obtained by solvability conditions applied in the following process. By substituting equation (53) into equation (52):

\begin{matrix} \frac{\partial^{2}}{\partial {T_{0}}^{2}} p_{1} (T_{0}, T_{1}) + ω_{n}^{2} p_{1} (T_{0}, T_{1}) = \\ - 2 ω_{n} i (\frac{\partial A (T_{1})}{\partial T_{1}} e^{i ω_{n} T_{0}} - \frac{\partial \bar{A} (T_{1})}{\partial T_{1}} e^{i ω_{n} T_{0}}) - \\ c_{deq} Ω_{f} (A (T_{1}) i ω_{n} e^{i ω_{n} T_{0}} - \bar{A} (T_{1}) i ω_{n} e^{- i ω_{n} T_{0}}) \\ - 2 ζ ω_{n} (A (T_{1}) i ω_{n} e^{i ω_{n} T_{0}} - \bar{A} (T_{1}) i ω_{n} e^{- i ω_{n} T_{0}}) + \\ f_{Leq} Ω_{f}^{2} (B (T_{1}) e^{i Ω_{f} T_{0}} + \bar{B} (T_{1}) e^{- i Ω_{f} T_{0}}) \\ + f_{peq} (\frac{- i ω_{n} χ \bar{A} R e^{- i ω_{n} T_{0}}}{iR C_{p} ω_{n} - 1} - \frac{i ω_{n} χ AR e^{i ω_{n} T_{0}}}{iR C_{p} ω_{n} + 1} + \bar{K} e^{- \frac{T_{0}}{R C_{p}}}), \\ \frac{\partial^{2}}{\partial {T_{0}}^{2}} q_{1} (T_{0}, T_{1}) + Ω_{f}^{2} q_{1} (T_{0}, T_{1}) \\ = - 2 Ω_{f} i (\frac{\partial B (T_{1})}{\partial T_{1}} e^{i Ω_{f} T_{0}} - \frac{\partial \bar{B} (T_{1})}{\partial T_{1}}) e^{- i Ω_{f} T_{0}} - \\ η Ω_{f} ({(B (T_{1}) e^{i Ω_{f} T_{0}} + \bar{B} (T_{1}) e^{- i Ω_{f} T_{0}})}^{2} - 1) \\ (i Ω_{f} B (T_{1}) e^{i Ω_{f} T_{0}} - i Ω_{f} \bar{B} (T_{1}) e^{- i Ω_{f} T_{0}}) + \\ G (- A (T_{1}) ω_{n}^{2} e^{i ω_{n} T_{0}} - \bar{A} (T_{1}) ω_{n}^{2} e^{- i ω_{n} T_{0}}) \end{matrix}

(54)

By assuming $Ω_{f} = ω_{n} + ε σ$ , where $σ$ is the detuning parameter, and by rearranging the outcome

\begin{array}{l} \frac{\partial^{2}}{\partial T_{0}^{2}} p_{1} (T_{0}, T_{1}) + ω_{n}^{2} p_{1} (T_{0}, T_{1}) = (- 2 ω_{n} \frac{\partial A (T_{1})}{\partial T_{1}} i - c_{d e q} Ω_{f} A (T_{1}) i ω_{n} - 2 ζ ω_{n} A (T_{1}) i ω_{n} + \\ f_{L e q} Ω_{f}^{2} B (T_{1}) e^{i σ ε T_{0}} - f_{p e q} \frac{i ω_{n} χ A R}{i R C_{p} ω_{n} + 1}) e^{i ω_{n} T_{0}} + N S T + c c, \\ \frac{\partial^{2}}{\partial T_{0}^{2}} q_{1} (T_{0}, T_{1}) + Ω_{f}^{2} q_{1} (T_{0}, T_{1}) = (- 2 i Ω_{f} \frac{\partial B (T_{1})}{\partial T_{1}} - η Ω_{f}^{2} B^{2} (T_{1}) \bar{B} (T_{1}) i + η Ω_{f}^{2} B (T_{1}) i - \\ G ω_{n}^{2} A (T_{1}) e^{- ε σ T_{0} i}) e^{i Ω_{f} T_{0}} + N S T + c c \end{array}

(55)

CC indicates conjugate of complex term, and NST represents all terms which are not secular. For solvability condition the coefficient of secular term in equation (55), that is, the coefficient of $e^{i ω_{s} T_{0}}$ and $e^{i Ω_{f} T_{0}}$ must become equal to zero: By putting it equal to zero and substituting $A (T_{2})$ as polar, $A (T_{1}) = \frac{1}{2} a (T_{1}) e^{i α (T_{1})}$ , and $B (T_{1}) = \frac{1}{2} b (T_{1}) e^{i β (T_{1})}$ , one obtains

\begin{matrix} - \overset{\cdot}{a} ω_{n} e^{i α} i + a \overset{\cdot}{α} ω_{n} e^{i α} - \frac{c_{deq} Ω_{f}}{2} ia ω_{n} e^{i α} - ζ ia ω_{n}^{2} e^{i α} \\ + \frac{f_{Leq} Ω_{f}^{2}}{2} b e^{i β} e^{i σ T_{1}} - \frac{1}{2} f_{peq} \frac{i ω_{n} χ aR}{iR C_{p} ω_{n} + 1} e^{i α} = 0, \\ - Ω_{f} \overset{\cdot}{b} e^{i β} i + Ω_{f} \overset{\cdot}{β} b e^{i β} - \frac{η Ω_{f}^{2}}{8} b^{3} i e^{i β} + \frac{η Ω_{f}^{2}}{2} b i e^{i β} \\ - \frac{G ω_{n}^{2} a}{2} e^{i α} e^{- i σ T_{1}} = 0 \end{matrix}

(56)

By multiplying $e^{- i α}$ at the first equation and multiplying $e^{- i β}$ at the second equation

\begin{array}{l} - \dot{a} ω_{n} i + a \dot{α} ω_{n} - \frac{c_{d e q} Ω_{f}}{2} i a ω_{n} - ζ i a ω_{n}^{2} + \frac{f_{L e q} Ω_{f}^{2}}{2} b e^{i (β - α + σ T_{1})} - \\ \frac{f_{p e q} ω_{n}^{2} R^{2} χ C_{p}}{2 (R^{2} C_{p}^{2} ω_{n}^{2} + 1)} - \frac{f_{p e q} ω_{n} R χ a}{2 (R^{2} C_{p}^{2} ω_{n}^{2} + 1)} i = 0, \\ - Ω_{f} \dot{b} i + Ω_{f} \dot{β} b - \frac{η Ω_{f}^{2}}{8} b^{3} i + \frac{η Ω_{f}^{2}}{2} b i - \frac{G ω_{n}^{2} a}{2} e^{i (α - β - σ T_{1})} = 0 \end{array}

(57)

By separating the real part and imaginary part of equation (28) and equating the outcome equal to zero

\begin{matrix} \overset{\cdot}{α} = + \frac{f_{peq} ω_{n} R^{2} χ C_{p}}{2 (R^{2} {C_{p}}^{2} {ω_{n}}^{2} + 1)} - \frac{f_{Leq} Ω_{f}^{2}}{2 a ω_{s}} b \cos (β - α + σ T_{1}), \\ \overset{\cdot}{a} = - \frac{c_{deq} Ω_{f}}{2} a - ζ a ω_{n} - \frac{f_{peq} a χ R}{2 (R^{2} {C_{p}}^{2} {ω_{n}}^{2} + 1)} \\ + \frac{f_{Leq} Ω_{f}^{2}}{2 ω_{n}} b \sin (β - α + σ T_{1}), \\ \overset{\cdot}{β} = \frac{G ω_{n}^{2} a}{2 Ω_{f} b} \cos (α - β - σ T_{1}) \\ \overset{\cdot}{b} = - \frac{η Ω_{f} b}{8} (b^{2} - 4) - \frac{G ω_{n}^{2} a}{2 Ω_{f}} \sin (α - β - σ T_{1}) \end{matrix}

(58)

Now by assuming $ψ = β - α + σ T_{1}$

\begin{matrix} \overset{\cdot}{ψ} = (\frac{G ω_{n}^{2} a}{2 Ω_{f} b} + \frac{f_{Leq} Ω_{f}^{2} b}{2 a ω_{n}}) \cos (ψ) - \frac{f_{peq} ω_{n} R^{2} χ C_{p}}{2 (R^{2} {C_{p}}^{2} {ω_{n}}^{2} + 1)} + σ \\ \overset{\cdot}{a} = - \frac{c_{deq} Ω_{f}}{2} a - ζ a ω_{n} - \frac{f_{peq} a χ R}{2 (R^{2} {C_{p}}^{2} {ω_{n}}^{2} + 1)} + \frac{f_{Leq} Ω_{f}^{2}}{2 ω_{n}} b \sin (ψ), \\ \overset{\cdot}{b} = - \frac{η Ω_{f} b}{8} (b^{2} - 4) + \frac{G ω_{n}^{2} a}{2 Ω_{f}} \sin (ψ) \end{matrix}

(59)

By putting derivative of $a, b$ and $ψ$ equal to zero in equation (59), the equation governing equilibrium solution amplitude will be obtained

6. Results and discussion

The geometrical and mechanical properties of the configurations have been shown in Table 1.

Table 1.

Geometrical and mechanical properties of the harvesters.

$h_{s} (mm)$	$L_{s} (mm)$	$E_{s} (Gpa)$	$b_{s} (mm)$	$ρ_{s} (kg / m^{3})$	$D (mm)$
0.635	267	73	32.5	2730	39.6
$h_{p} (mm)$	$L_{1 p} (mm)$	$E_{p} (Gpa)$	$b_{p} (mm)$	$ρ_{p} (kg / m^{3})$	$d_{31} (pC / N)$
0.534	0	66	25.4	7800	−190
Q	$L_{2 p} (mm)$	$ϵ_{33} (nF / m)$	$e_{31} (C / m^{2})$	$ρ_{c} (kg / m^{3})$	$ρ_{f} (kg / m^{3})$
100	31.8	13.28	−12.54	64	1.2
R (MΩ)	$L_{c} (mm)$	$E_{c} (Mpa)$
2.46	203	300

A comparison between the natural frequency of the system evaluated using Galerkin method and evaluated value using assumed modes method has been presented in Table 2. It must be noted that the length of the substrate layer are identical in all cases. Considering this table if the mode shapes of a uniform beam that is, the beam without both piezoelectric layer and the cylindrical shape are used as comparison function, then the evaluated amount for natural frequency does not converge on the obtained value by assumed mode method. It shows that even using three mode shapes of the uniform beam as comparison function, the convergence has not been occurred. It shows that if one uses two mode shapes of the beam with piezoelectric layer and without the cylindrical shape as comparison function, then the results converge on the value obtained by assumed mode method. Moreover, it shows that the convergence occurs only using one free vibration mode shape of the beam with both piezoelectric layer and equivalent cylindrical shape. As mentioned in Section 2 the dimensions of the equivalent cylinder have been considered identical to the cylinder dimension of the main system for configuration 1. In addition, the cylinder whose side length is equal to the diameter of the main cylinder, and its length is equal to the length of the main cylinder has been considered as equivalent cylinder in configuration 2. Here, the mode shapes of the equivalent system are called semi-exact mode shape. The comparison between the mode shape obtained using Galekin method and using assumed mode method has been shown in Figures 3 and 4 for configuration 1 and 2, respectively. It shows that there is a good agreement between the mode shape obtained by assumed mode method and the Galerkin method using one mode shape of semi-exact mode shape as comparison function.

Table 2.

First natural frequency of the system using different comparison functions.

The type ofconfiguration	Galerkin method using one mode shape of the beam with considering the effect of piezoelectric layer and without considering the effect of the cylinder in comparison function	Galerkin method using three mode shape of the beam with considering the effect of piezoelectric layer and without considering the effect of the cylinder in comparison function	Galerkin method using one mode shape of the beam without considering the effect of both piezoelectric layer and cylinder in comparison function	Galerkin method using three mode shape of the beam without considering the effect of both piezoelectric layer and cylinder in comparison function	Galerkin method using one mode shape of the beam with considering the effect of both piezoelectric layer and cylinder in comparison function	Assumed mode method
Configuration 1, $E_{c} = 100 M . Pas$	3.237	1.181	2.783	1.625	1.101	1.106
Configuration 1. $E_{c} = 1000 M . Pas$	9.726	1.214	7.771	1.771	1.101	1.106
Configuration 2, $E_{c} = 100 M . Pas$	1.690	1.619	2.238	1.697	1.606	1.604
Configuration 2, $E_{c} = 1000 M . Pas$	2.265	1.622	2.476	1.702	1.606	1.604

Figure 5 compares the amplitude of the harvested voltage evaluated by Galerkin method with the evaluated value by assumed mode method in different values of wind speed for configuration 1. To this aim the discretized equations by Galerkin procedure that is, equations (26)–(28), and discretized equations by assumed mode procedure that is, equations of (46)–(48) are solved using Runge–Kutta–Fehlberg numerical method, which is called RKF45. Variations of steady state amplitudes of the structure and wake oscillations with respect to the variations of wind speed are extracted for different values of the wind speed. To this aim, the value of $Ω_{f}$ , which is closely related to the wind speed, is considered much less than the natural frequency of linear system. Then, solving the equations of motion and by obtaining time history of solution, the amplitude of the steady state is recorded. By slightly increasing the value of $Ω_{f}$ that is, by slowly increasing the wind speed the previous step is repeated. The steady-state amplitude of the previous step is considered as the initial conditions of displacement for the next step, whilst the initial condition of the velocity is settled to zero. This work is followed until $Ω_{f}$ largely overpass the system natural frequency. It shows that there is a main difference between the amount of evaluated voltage by Galerkin method and assumed mode method when three mode shapes of a straight uniform beam are used as comparison function. It demonstrates that the lock-in domain evaluated by Galerkin method occurs at more flow speed in compare to the assumed mode method. Moreover it shows that the amount of the harvested voltage is approximately double times of the obtained value by assumed mode method. It shows that the difference decreases when two mode shapes of the beam with piezoelectric layer and without equivalent cylindrical shape are used as comparison function. Also it indicates that there is a good agreement between the Galerkin method and assumed mode method, if one uses only one mode shape of the beam with both piezoelectric layer and equivalent cylindrical shape (semi-exact mode shape).

Comparing Figures 5 and 6 shows that an increase of the value of the elasticity modulus of the cylinder does not change the evaluated harvested voltage. It is due to the fact that considering Figure 2, the deflected beam at the end of the cantilever beam that is, at the cylinder part is as straight line without curvature, and therefore the bending stiffness does not have a considerable effect.

The comparison of the evaluated harvested voltage using Galerkin method with the evaluated harvested voltage using assumed mode method has been shown for configuration 2 in Figure 7. It demonstrates that there is no a main difference between evaluated lock-in domain by Galerkin method and assumed mode method when three(three)(one) free vibration mode shape of beam without (with)(with) piezoelectric layer and without (without)(with) the equivalent cylindrical shape are used as comparison function. However, it shows that the harvested voltage evaluated by using three mode shapes of the beam without both piezoelectric layer and equivalent cylinder mass are more than other cases. Comparing Figures 6 and 7 demonstrates that the difference between different curves of Figure 6 is larger and considerable than the difference between different curves of Figure 7. These variations is due to the fact that the amount of evaluated natural frequency using three mode shapes of the straight beam (two mode shapes of the beam coated by piezoelectric layer) (one mode shape of the beam coated by piezoelectric layer and attached to the equivalent cylindrical shape) converge on the value obtained by assumed mode method where it has not been occurred for configuration 1 using three mode shapes of straight beam. It is due to the fact that the mass distribution of the cylinder in configuration 1 is scattered in the axial direction more than configuration 2, where the cylinder is perpendicular to the beam and the displacement of different points of the cylinder is approximately near together.

In follow the numerical solution of discretized equations by Galerkin method are compared to the analytical solution obtained by perturbation method. The comparison solution has been shown in Figure 8 for configuration 1. It demonstrates a good agreement between the results. It means that the numerical solution is verified.

The comparison between the harvested voltage by configuration 1 with configuration 2 has been shown in Figure 9. It shows that the amount of harvested voltage by configuration 1 is more than the harvested voltage by configuration 2. However, this figure shows that the horizontal cylinders needs more speed of the wind flow to stay in lock in domain.

Figures 10 and 11 present a comparison of output power between two energy harvester configurations. According to the results shown in these figures, Configuration 2 yields higher output power than Configuration 1 within the lock-in domain. Moreover, the figures demonstrate that when the resistance value is taken as R₀ = 2.46 MΩ, an increase in electrical resistance from 0.1R₀ to 0.5R₀ and then to R₀ leads to a corresponding increase in output power, provided that all other parameters—such as the mechanical and geometrical properties of the system—are held constant. However, when the electrical resistance is further increased to higher values (e.g. 5R₀ and 10R₀), the output power decreases. In other words, maximum power output is achieved at an optimal resistance value, approximately equal to R₀.

If the wind speed values on the horizontal axis in Figures 5 –11 are multiplied by a specific factor, the wind speed will be expressed in a dimensionless form. In this paper, this multiplier is set to $D / (2 π S_{T} UT)$ for both Configuration 1 and Configuration 2. It is important to note that the same multiplier is used for both configurations because they have been non-dimensionalized based on the same parameters outlined in Table 1 and equation (12).

As shown in the previous figures, at low wind speeds, the harvested voltage amplitude remains very small. However, as the wind speed gradually increases, a critical value is reached at which the amplitude suddenly rises to a significantly higher level. Subsequently, as the wind speed continues to increase, the amplitude remains high until a second critical value is reached, beyond which the amplitude abruptly drops back to a very low level. Figure 12 demonstrates a similar behavior during the decreasing the amount of the wind speed. Notably, the wind speeds at which the sudden changes in amplitude occur differ from those in the increasing phase. In other words, a clear hysteresis phenomenon is observed between the forward and backward wind speed variations. It must be noted that previous experimental work have indicated that a cylinder under vortex induced vibration exhibits an upper and a lower response branches (Khalak and Williamson, 1996), the present investigation reveals that these branches cannot be accurately predicted using the classical Van der Pol oscillator model used in this paper. This suggests that further studies employing alternative modeling approaches may be necessary to better capture and analyze this behavior.

7. Conclusion

In this paper a new mathematical modelling for two configurations of piezoelectric energy harvester under vortex induced vibration has been presented. At the first (second) configuration, a cylinder under vortex flow is attached to the end of the cantilever beam as parallel (vertical) to the beam. In contrast to previous works that for applying the Galerkin method, the common cross section of the beam and the cylinder has been considered as the boundary position, here the end point of the configuration has been considered as the boundary position. It has been shown that using this mathematical modelling, the complexity of boundary conditions may be removed since the boundary conditions will be appeared as geometric. Also, it has been shown that the provided notation of modelling may be used for both directions of the cylinder.

The discretized equations of motion have been solved using both numerical and also perturbation methods, which is a new technique for the configuration with vertical cylinder. The results showed that the perturbation method can predict the dynamic behavior of system correctly. For more verification the system has been discretized according the assumed mode method presented in previous work. The comparison between Galerkin method and assumed mode method showed that there is a good convergence between the Galerkin method and assumed mode method when three(three)(one) free vibration mode shape of beam without (with)(with) piezoelectric layer and without (without)(with) the equivalent cylindrical shape are used as comparison function in the configuration with vertical cylinder to the beam.

It shows that the amount of harvested voltage by configuration with cylinder along the lengths of the beam is more than the harvested voltage by configuration with cylinder vertical to the beam. However, it has been shown that the former configuration needs more speed of the wind flow to stay in lock in domain than later configuration.

Footnotes

Declaration of conflicting interests

The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The study was performed in Department of Mechanical Engineering, Kharazmi University, under a grant presented by Vice Chancellor in Research, which should be acknowledged.

ORCID iD

Mehdi Zamanian

Data availability statement

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

References

Abdelkefi

Najar

Nayfeh

, et al. (2011) An energy harvester using piezoelectric cantilever beams undergoing coupled bending–torsion vibrations. Smart Materials and Structures 20: 115007.

Akaydin

Elvin

Andreopoulos

(2012) The performance of a self-excited fluidic energy harvester. Smart Materials and Structures 21: 025007.

Alimanesh

Zamanian

(2024) Analysis of clamped-clamped piezoelectric energy harvester under vortex induced vibration considering the stretching effect. Journal of Intelligent Material Systems and Structures 35: 333–351.

Dai

Abdelkefi

Wang

(2014) Piezoelectric energy harvesting from concurrent vortex-induced vibrations and base excitations. Nonlinear Dynamics 77: 967–981.

Dai

Abdelkefi

Yang

, et al. (2016) Orientation of bluff body for designing efficient energy harvesters from vortex-induced vibrations. Applied Physics Letters 108: 053902.

Facchinetti

de Langre

Biolley

(2004) Coupling of structure and wake oscillators in vortex-induced vibrations. Journal of Fluids and Structures 19: 123–140.

Firouzi

Zamanian

Hosseini

SAA

(2016) Static and dynamic responses of a microcantilever with a T-shaped tip mass to an electrostatic actuation. Acta Mechanica Sinica 32: 1104–1122.

Hagedorn

DasGupta

Vibrations and Waves in Continuous Mechanical Systems. Chichester: John Wiley & Sons Ltd, 2007.

Hou

Shan

Zhang

, et al. (2020) Design and modeling of a magnetic-coupling monostable piezoelectric energy harvester under vortex-induced vibration. IEEE Access 8: 108913–108927.

10.

Huang

Zhong

(2023) Hydrokinetic energy harvesting from flow-induced vibration of a hollow cylinder attached with a bi-stable energy harvester. Energy Conversion and Management 278: 116718.

11.

Javed

Abdelkefi

(2019) Characteristics and comparative analysis of piezoelectric-electromagnetic energy harvesters from vortex-induced oscillations. Nonlinear Dynamics 95: 3309–3333.

12.

Jia

Shan

Upadrashta

, et al. (2018) Modeling and analysis of upright piezoelectric energy harvester under aerodynamic vortex-induced vibration. Micromachines 9: 667.

13.

Khalak

Williamson

CHK

(1996) Dynamics of a hydroelastic cylinder with very low mass and damping. Journal of Fluids and Structures 10: 455–472.

14.

Lai

Wang

Zhu

, et al. (2021) A hybrid piezo-dielectric wind energy harvester for high-performance vortex-induced vibration energy harvesting. Mechanical Systems and Signal Processing 150: 107212.

15.

Liang

Hao

Olszewski

(2021) A review on vibration-based piezoelectric energy harvesting from the aspect of compliant mechanisms. Sensors and Actuators 331: 112743.

16.

Rezaei

Talebitooti

(2019) Wideband PZT energy harvesting from the wake of a bluff body in varying flow speeds. International Journal of Mechanical Sciences 163: 105135.

17.

Shan

Deng

Song

, et al. (2017) A piezoelectric energy harvester with bending–torsion vibration in low-speed water. Applied Sciences 7: 116.

18.

Sui

Zhang

Yang

, et al. (2022) Modeling and experimental investigation of magnetically coupling bending-torsion piezoelectric energy harvester based on vortex-induced vibration. Journal of Intelligent Material Systems and Structures 33: 1147–1160.

19.

Sun

Seok

(2019b) Development of the optimal bluff body for wind energy harvesting using the synergetic effect of coupled vortex induced vibration and galloping phenomena. International Journal of Mechanical Sciences 156: 435–445.

20.

Sun

Zhao

Tan

, et al. (2019a) Low velocity water flow energy harvesting using vortex induced vibration and galloping. Applied Energy 251: 113392.

21.

Lin

(2020) Design and analysis of a vortex-induced bi-directional piezoelectric energy harvester. International Journal of Mechanical Sciences 173: 105457.

22.

Wang

(2021) Development of a non-linear bidirectional vortex-induced piezoelectric energy harvester with magnetic interaction. Sensors 21: 2299.

23.

Tang

Fan

Wang

, et al. (2022) Energy harvesting from flow-induced vibrations enhanced by meta-surface structure under elastic interference. International Journal of Mechanical Sciences 236: 107749.

24.

Wang

Geng

Ding

, et al. (2020) The state-of-the-art review on energy harvesting from flow-induced vibrations. Applied Energy 267: 114902.

25.

Zamanian

Garibaldi

(2019) Vortex induced vibration analysis of a cylinder mounted on a flexible rod. Wind and Structures 29: 471–485.

26.

Zhang

Abdelkefi

Dai

, et al. (2017a) Design and experimental analysis of broadband energy harvesting from vortex-induced vibrations. Journal of Sound and Vibration 408: 210–219.

27.

Zhang

Dai

Abdelkefi

, et al. (2017b) Improving the performance of aeroelastic energy harvesters by an interference cylinder. Applied Physics Letters 111: 073904.

28.

Zhao

Zhang

Wang

(2019) Optimization of galloping piezoelectric energy harvester with V-shaped groove in low wind speed. Energies 12: 4619.