Access the full text.

Sign up today, get DeepDyve free for 14 days.

Acoustics
, Volume 1 (4) – Sep 24, 2019

/lp/multidisciplinary-digital-publishing-institute/effect-of-tip-mass-length-ratio-on-low-amplitude-galloping-T0hBzXwr2C

- Publisher
- Multidisciplinary Digital Publishing Institute
- Copyright
- © 1996-2021 MDPI (Basel, Switzerland) unless otherwise stated Disclaimer The statements, opinions and data contained in the journals are solely those of the individual authors and contributors and not of the publisher and the editor(s). MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Terms and Conditions Privacy Policy
- ISSN
- 2624-599X
- DOI
- 10.3390/acoustics1040045
- Publisher site
- See Article on Publisher Site

acoustics Article Eect of Tip Mass Length Ratio on Low Amplitude Galloping Piezoelectric Energy Harvesting Mohammad Yaghoub Abdollahzadeh Jamalabadi Department of Mechanical, Robotics and Energy Engineering, Dongguk University-Seoul, Seoul 04620, Korea; abdollahzadeh@dongguk.edu Received: 21 August 2019; Accepted: 20 September 2019; Published: 24 September 2019 Abstract: Galloping beams were exposed to the wind free stream and is used for sustainable wind-power harnessing. In this paper, the eect of tip mass on the performance of a galloping energy harvester is investigated by simple modeling of the system, which is useful for broad engineering applications of these systems. Here, the piezoelectric layer attached to a cantilever beam with a tip mass exposed to the wind free stream is used as an energy harvester. A ﬂuid–solid interaction model is used to simulate the problem. The ﬂuid–solid interaction model is composed of the experimental data for aerodynamic loads and one-dimensional structural model of piezoelectric and beam material with Euler–Bernoulli beam theory. The governing partial dierential equations of the system are solved analytically by use of the approximation method. The resulting model is conﬁrmed by preceding experimental results. The eects of the tip mass length ratio on the onset of galloping, the level of the produced voltage, and the harvested power are determined analytically. As shown by increase of the length of tip mass for the constant beam and piezoelectric length, the inertia of the system increases while the tip displacement and onset of galloping decrease. Keywords: Wind energy; sustainable energy; galloping; piezoelectric; energy harvesting 1. Introduction The elasticity of piezoelectric materials makes them possible energy harvesting materials from environmental vibrations. A review of vibration-based micro-generator piezoelectric energy harvesters by Saadon and Sidek [1] showed the wide applicability of piezoelectric power harvesting devices. The piezoelectric material can convert large amounts of strain to electrical energy. Two modes are important in the mechanical energy conversion of piezoelectric materials. One mode is based on applying the external force parallel to the poling direction, and the second is perpendicular. As the perpendicular mode has the lower coupling coecient, it is used mostly in engineering applications. Saadon and Sidek performed a review of vibration-based micro-electro-mechanical systems (MEMS) piezoelectric energy harvesters [1]. To convert the environmental excitations to strain, usually the bimorph conﬁguration is used. An experimental investigation of the bimorph cantilever model for piezoelectric energy harvesting from base excitations was done by Erturk and Inman [2]. Many environmental sources such as raindrop impacts [3], sea waves [4,5], and galloping [6–9] could be used as a vibrational source of this device. Nonlinear modeling and analysis of piezoelectric energy harvesting from transverse galloping [10,11] and its coupling with electric loads [12] are performed in various studies. The eciency and optimal design of piezoelectric energy harvesting induced by galloping phenomena are typically explored mathematically. Although the models of References [13,14] are more accurate than the original model of Reference [9], they achieve unrealistic results though their simulations were not validated by proper experimental results. The use of the experiment results of Reference [8] as a benchmark is not a proper choice as the eect of tip mass is negligible in that Acoustics 2019, 1, 763–793; doi:10.3390/acoustics1040045 www.mdpi.com/journal/acoustics Acoustics 2019, 1 764 study. The models presented in References [9,13–15] caused a great error in the calculation of harvested power which could be misleading for young researchers in the design of energy harvesting devices by the aid of piezoelectric materials. Although Sirohi and Mahadik [9,16] showed that their results are in good agreement with their numerical calculations, it is suspicious that they reached identical analytical results with two dierent methods. They used the formula for vibrating base excitation which does not include the eect of tip mass inertia [16] and coupled formulation [9], but they reached the same analytical results. Based on the results of Reference [9,13,14], the order of magnitude of power estimated is watts while the author ’s experience [17,18] and benchmark [8] show maximum attainable power in milliWatts. Also, the eect of tip mass is studied in References [19,20]. The cases studied in Reference [19,20] are related to the cases where the blu body subjected to wind ﬂows or acceleration is concentrated at tip mass, and thus, the eect of the mass moment of inertia on dynamic calculation and the length of tip mass on the applied force and momentum is negligible. The degree of freedom of movement of the cantilever beam in References [19–39] is in a way that it cannot sense the eect of distributed tip mass and that the concentrated mass can lead to a sucient solution. Modeling and experimental veriﬁcation of proof mass eects on vibration energy harvester performance by Kim et al. [19] and a comparative study of tip cross sections for ecient galloping energy harvesting by Yang et al. [20] are about the geometrical optimization of the system. Parkinson et al. investigated the various aspects of ﬂow-induced vibration such as the following: square-section cylinder [21], the combined eects of vortex-induced vibration and galloping [22–24], turbulence eects on galloping of blu cylinders [25], wind-induced instability [26], aeroelastic instability [27–29], wake source model for blu body potential ﬂow [30], square prism [31], galloping response of towers [32], and combined eects of galloping and vortex resonance [33]. Since Den Hartog [34] presented the galloping criterion in transmission line vibration due to sleet, the aeroelastic galloping of prismatic bodies has been studied in many researches [35]. Recently, Jamalabadi and Kwak presented a dynamic model of a galloping structure equipped with piezoelectric wafers and energy harvesting [36]. Their model is used in broad engineering applications. Also, Abdelkeﬁ et al. [37], Ewere and Wang [38], and Rostami and Fernandes modeled the eect of inertia and ﬂap on autorotation applied for hydrokinetic energy harvesting [39]. Although their model is comprehensive, it fails to present a simple form to use as a rule of thumb by engineers. This paper is about the mass and the eect of geometry (tip mass length ratio) of the mass, which can aect the applied force and momentum. To compare with the results of Reference [8], the D-shaped tip mass has been considered. Although one case cannot conclude the eect of tip mass, as there was no experimental data for other shapes in the literature, just the modiﬁed versions of the previous method are presented. Dierent shapes of tip mass could be discussed in future researches. The manuscript has followed several studies on galloping energy harvesting, in particular, the conﬁguration which is similar to that presented in Reference [9]. In this paper, the problem is solved completely and is correctly addressed. Furthermore, the approximate solution method is used to give an explicit solution to the problem for engineers who do not prefer to involve the mathematical complexity of the common methods for partial dierential Equation (PDE) solutions. 2. Mathematical Modeling As presented in Figure 1, the system is composed of a cylinder tip mass (blu body) mounted on an elastic beam made of aluminum and covered by two piezoelectric sheets (sheets are bonded on both sides) in an open environment with the wind velocity of V . The two piezoelectric plates plotted in Figure 1 are connected in parallel with opposite polarity to an electrical impedance R to harvest higher galloping energy. The free-body diagram of the tip mass which is exposed to an incoming ﬂow in the z-direction with a magnitude of V is depicted in the top view of the problem for the D-shape cylindrical mass in Figure 2. The beam mounted blu body experiences galloping in the x-y plane in they direction when V is greater than the onset galloping velocity. The material properties of the system (cantilever beam, air, and the tip body) are presented in Table 1. As the system is exposed to Acoustics 2019, 2 FOR PEER REVIEW 3 system is exposed to the wind, it starts to vibrate (in the y-direction) to reach a stable oscillation. The strain produced in the base beam is converted to an electric charge in piezoelectric sheets. Acoustics 2019, 1 Table 1. Geometrical and material properties. 765 Symbol Description and Unit Value the wind, it starts to vibrate (in the y-direction) to reach a stable oscillation. The strain produced in the −3 ρa Air density (kg·m) 1.225 base beam is converted to an electric charge in piezoelectric sheets. mt Tip mass (g) 65 Lr Length of the tip body (mm) 235 Table 1. Geometrical and material properties. D Width of the tip body (mm) 30 Symbol Description and Unit Value V∞ Wind velocity (m/s) 4.02 Air density (kgm ) 1.225 L Length of the beam (mm) 90 m Tip mass (g) 65 L Length of the tip body (mm) 235 wb Width of the beam material layer (mm) 38 D Width of the tip body (mm) 30 tb Beam material layer thickness (mm) 0.635 V Wind velocity (m/s) 4.02 −2 Eb Beam material Young's modulus (GN·m) 70 L Length of the beam (mm) 90 w Width of the beam material layer (mm) 38 −3 ρb Beam material density (kg·m) 2700 t Beam material layer thickness (mm) 0.635 Lp Length of the piezoelectric sheets (mm) 72.2 E Beam material Young’s modulus (GNm ) 70 wp WidtBeam h of the piezoel material density ectric la(kg yer (mm) m ) 36.2 L Length of the piezoelectric sheets (mm) 72.2 tp Piezoelectric layer thickness (mm) 0.267 w Width of the piezoelectric layer (mm) 36.2 −2 Ep Piezoelectric material Young's modulus (GN·m) 62 t Piezoelectric layer thickness (mm) 0.267 E Piezoelectric material Young’s modulus−(GN 3 m ) 62 ρp Piezoelectric material density (kg·m) 7800 Piezoelectric material density (kgm ) 7800 −1 d31 Strain coefficient of the piezoelectric layer (pC N ) −320 d Strain coecient of the piezoelectric layer (pC N ) 320 −11 ε33 " Permittivity Permittivity cocomponent mponent at co at constant nstant strain (n strain (nF F·mm) 33.6 ) 33.6 R Load resistance (MW) 0.7 R Load resistance (MΩ) 0.7 Figure 1. Front view of the problem. Figure 1. Front view of the problem. Acoustics 2019, 1 766 Acoustics 2019, 2 FOR PEER REVIEW 4 Figure 2. Top view of the problem for a D-shape cylindrical mass. Figure 2. Top view of the problem for a D-shape cylindrical mass. 2.1. Beam Modeling 2.1. Beam Modeling The frequency of the vibration is equal to the fundamental mode of the beam (here, it is the The frequency of the vibration is equal to the fundamental mode of the beam (here, it is the bending mode and is a function of the Reynolds number as claimed in Reference [8]). The solid part of bending mode and is a function of the Reynolds number as claimed in Reference [8]). The solid part the system is modeled as follow: of the system is modeled as follow: The displacement of the beam is a function of time and is expressed using the modal expansion The displacement of the beam is a function of time and is expressed using the modal expansion form (exact mode shapes of the structure): form (exact mode shapes of the structure): w(x, t) = (x)q (t) (1) i i wx (,t) = φ(x)q(t) (1) ii where w(x,t) is the displacement of the beam; ' (t) is the ith mode shape; and q (t) is the time-dependent i i where w(x,t) is the displacement of the beam; φi(t) is the ith mode shape; and qi(t) is the part of the beam displacement, which is also mentioned in the literature as the mode coordinate. time-dependent part of the beam displacement, which is also mentioned in the literature as the Each of the displacement functions is named a mode, and the shape of the displacement curve is mode coordinate. Each of the displacement functions is named a mode, and the shape of the named the mode shape. The Euler–Bernoulli beam model is composed of the second derivate of the displacement curve is named the mode shape. The Euler–Bernoulli beam model is composed of the internal moment; piezoelectric coupling term; and internal damping of the structure (strain rate), second derivate of the internal moment; piezoelectric coupling term; and internal damping of the where the galloping aerodynamic moment and force are applied by the Dirac delta function at the structure (strain rate), where the galloping aerodynamic moment and force are applied by the Dirac end of the beam and the viscous air damping coecient is neglected. The exact mode shapes of the delta function at the end of the beam and the viscous air damping coefficient is neglected. The exact structure are obtained by Euler–Bernoulli beam theory, including forcing, damping, and piezoelectric mode shapes of the structure are obtained by Euler–Bernoulli beam theory, including forcing, coupling terms: damping, and piezoelectric coupling terms: 2 2 2 3 @ w @ @ w @ w A + I E + c = 2 2 2 2 @t22 @x @x 2@x @t 3 h i ∂∂ww ∂ ∂w x =L 1 2 ρAI++E c = (2) E d w (t + t ) (x x ) V(t)+ p 31 p22 p 2 1 2 dx x =L ∂∂tx ∂x ∂ 1x 1∂t M (x L) + F (x L) tip tip xL = dx d 0 (2) Ed w () t+− t δ(x x) V(t)+ where E is Young s modulus of material, pp 31 pI isb the mass moment 1 of inertia (calculated with respect to dx xL = the axis perpendicular to the applied loading and passes through the centroid of the cross-section), d 1 1 c is the internal damping coecient of structure (with the unit of kgm s ), is the density of MFδδ (x−+ L) (x− L) tip tip dx material, A is the cross-sectional area, E is the piezoelectric material’s Young s modulus, d is the p 31 strain coecient of piezoelectric, w is the width of the piezoelectric layer, t is the piezoelectric layer where E is Young′s modulus of mater p ial, I is the mass moment of inertia p (calculated with respect to thickness, L and L are respectively the starting and ending points of the piezoelectric material sheets the axis perp 1 endicu 2 lar to the applied loading and passes through the centroid of the cross-section), c −1 −1 on the beam, is the delta Dirac function, V is the produced voltage at the electrodes, M is the is the internal damping coefficient of structure (with the unit of kgm ·s ), ρ is the densi tip ty of eective moment applied at the end of the beam by the tip mass, and F is the eective force applied material, A is the cross-sectional area, Ep is the piezoelectric material’s Young′s modulus, d31 is the tip at the end of the beam by the tip mass. The left-hand side (LHS) of Equation (2) is the simpliﬁcation strain coefficient of piezoelectric, wp is the width of the piezoelectric layer, tp is the piezoelectric layer of the time-dependent linear theory of elasticity, known as engineer s beam theory (a special case thickness, L1 and L2 are respectively the starting and ending points of the piezoelectric material of Timoshenko beam theory or classical beam theory), which oers a tool for calculating the small sheets on the beam, δ is the delta Dirac function, V is the produced voltage at the electrodes, Mtip is the effective moment applied at the end of the beam by the tip mass, and Ftip is the effective force applied at the end of the beam by the tip mass. The left-hand side (LHS) of Equation (2) is the Acoustics 2019, 1 767 deﬂection characteristics of beams under lateral external loads. The ﬁrst term in the LHS characterizes the inertial eect derived from kinetic energy while the second one in the LHS signiﬁes the eective stiness derived from potential energy due to internal forces. The right-hand side of Equation (2) has units of force per length composed of two sources of a distributed load of piezoelectric sheet and tip mass eect at the point of tip mass-beam connection. The piezoelectric load considered here is the ideal case which will be discussed more in Section 2.5, and the tip mass eects are retrieved from the ﬁrst and second integrations of Equation (2) concerning the longitudinal direction which leads to shear force in the beam and to bending moment in the beam. Estimate the exact mode shapes of the structure by dropping forcing, damping, and piezoelectric coupling terms from Euler–Bernoulli beam theory. As in the absence of a transverse load, the free vibration of the beam leads to the sum of harmonic vibrations (based on the Fourier decomposition); here, a sinusoidal form for the time-dependent part is assumed, and then, the resulting eigenvalue equation for the exact mode shapes of for each value of frequency of the beam is obtained by a partial dierential equation: @ (x) E(x)I(x) = (x)A(x)! (x) (3) 2 2 @x @x where ! is the natural frequency of the ith mode. The left-hand side of Equation (3) presents the second derivate of the bending moment with respect to the longitudinal position, while the right-hand side of Equation (3) presents the inertial terms. The Rayleigh method oers the same method for computation of the fundamental frequency of the system. As the displacement is not unique and depends on the frequency, the time-dependent part is given by the following: h i .. . 2 0 0 q (t) + 2!q (t) + ! q (t) = (L ) (0) E d w (t + t ) V(t)+ i p p 31 p p i i i i h i 1 2 0 2 0 3 V Da q (t) (L)L + (L) (L)L + (L)L + a r 1 i 2 1 i i r r 2 3 4 3 0 2 (L)L + 2 (L) (L)L + (4) 6 r 7 6 i i i r 7 6 7 . 6 0 2 2 3 0 3 4 7 1 1 6 7 2 (L) (L)L + (L) (L)L + 6 7 V Da q (t) r r a 3 6 i i 7 i i 2 1 i 6 7 6 7 4 L 5 0 4 r (L) where is the air density, L is the length of the tip body, D is the width of the tip body, V is the a r 1 wind velocity, L is the length of the beam, t is the beam material layer thickness, dot is used to present the derivative with respect to time, and prime is used for the ﬁrst derivate versus location. The values of q (actually constants as there is one for each mode i), in general, are complex and are found by the initial conditions on the displacements and velocity of the beam. As the characteristic time scale of the ﬂow motion is less than the characteristic time scale of the oscillations, the quasi-steady hypothesis is used to evaluate the galloping aerodynamic terms in the above equations. 2.2. Aerodynamic Modeling The ﬂow-induced vibrations are divided into two parts as forced vibrations (vibrations in the along-wind direction) and self-induced vibrations (or aeroelasticity due to vibrations in the across-wind direction). Examples of forced vibrations are gust (turbulence eects), bueting, and vortex shedding (without “lock-in” eect), and examples of induced vibrations are vortex shedding (with “lock-in” eect), pipeline vibration nearby of the bed, galloping, and bridge ﬂutter. For the normal structures, the dynamic forces that must be considered based on engineering codes are the forces in the wind direction (usually turbulence and gusts), and vortex-induced oscillations such as ﬂutter and galloping are, as a whole, of no real importance. Although there is no detailed explanation of galloping in the American Society of Mechanical Engineers (ASME) Code, vibrations in the galloping eect were shortlisted as a possible cause of the Tacoma Narrows bridge failure in 1940 between the vibrations induced by ﬂuid ﬂow. The galloping, which in some references is referred to as “dancing vibrations”, also appears in the vibrations of group of tethers or risers on a tension-leg platform, vibration of Acoustics 2019, 1 768 asymmetric ice-coated power lines, and the vibration of a ﬂow line connected to the leg of an oshore tower, which are considered examples of galloping in engineering problems. In practice, structures always experience conditions of turbulent ﬂow. Bearman et al. [21] performed some experiments on the ﬂow-induced vibration of a square-section cylinder and measured the aerodynamic damping. They showed that galloping instability and vortex resonance (due to negative aerodynamic damping) in some proﬁles is inﬂuenced by the turbulence of the airﬂow (e.g., in a wind tunnel test or natural winds). The nonlinear aeroelastic behavior of the system causes the failure of the prediction of galloping oscillation scales. Also, suitably choosing an aerodynamically equivalent reference frame with the unsteady situation that symbolizes the steady condition is another problem. Models of the combined eects of vortex-induced oscillation and galloping are presented by Corless et al. [22–24]. They considered combined eects to capture the measured values of the amplitude of vibrations. Also, Laneville and Parkinson [25] showed that the diculty of calculating that value (by wind tunnel tests on galloping of blu cylinders) is due to the turbulence damping eects which increase the Scruton number. Wind-induced aeroelastic instability of towers and structures (as a nonlinear oscillator) are mathematically modeled through the works of Parkinson et al. [26–32]. Parkinson [26] discovered that the criterion for the quasi-steady assumption is that any wake impact experienced by the oscillating body through one period should not aect the next period of the motion of the body. It means that, if the body (with characteristic length of D and exposed to free stream velocity of V ) vibrates with the natural frequency of f (for a low frequency, it typically is 1 n about one Hertz) exposing a ﬂuid vortex impact at the beginning of a period of galloping, the vortex in ﬂuid (with the center velocity of V ) should be moved downstream adequately far away (at minimum, ten times the characteristic length of body) until the end of that period. Since, at one period later when the mass body returns to the location of the beginning of galloping, that vortex no longer disturbs the ﬂuid ﬂow around the body (f V /10D). Meanwhile, the frequency disturbance to the mean ﬂow n 1 generated by vortex shedding frequency (f = V /5D) should be, at minimum, two times greater than the oscillation frequency. Bearman et al. [21] have a more conservative limit for the natural frequency of the square prism structure (f V /30D). It means that the frequency disturbance to the mean ﬂow n 1 generated by the vortex shedding frequency should be, at minimum, six times greater than the natural oscillation frequency of the structure for the quasi-steady supposition to be appropriate. To obtain the force terms of the right-hand side of Equation (5) from Equation (2), the following is performed: Based on the free body diagram of the y-z plane, which is depicted in Figure 2, the force is as follows: [ ] F = V DL c cos + c sin ds (5) tip a r L D Similarly, momentum is as follows: M = V DL s[c cos + c sin]ds (6) tip a r L D which can be obtained by application by the wind on the tip mass. Table 2 presents the experimental coecients. As the lift and drag coecients are known, the galloping force coecient based on the angle of the wind normal to the surface (angle of attack) is obtained from the projection of lift and drag forces on the y-direction as follows: C = c cos + c sin (7) y L D in which the angle of wind normal to the surface (angle of attack) is deﬁned by w + sw L L = arctan (8) 1 Acoustics 2019, 1 769 In the literature for commonly employed blu bodies, the galloping force is given in the form of a third-order polynomial function as a function of the galloping position [13]. The ﬁrst three aerodynamic coecients are usually used to characterize the galloping force with zero coecients for even terms: 2 3 ! ! . . 6 7 1 6 y y 7 2 6 7 6 7 F = V DL a + a (9) y a r 6 1 3 7 1 4 5 2 V V 1 1 where a and a are aerodynamic empirical coecients, for which the D-shaped and other cylindrical 1 3 cross sections are found in Table 2. The force terms in Equation (2) are replaced by the terms appearing in Equation (4) by the integration of F over the tip mass: y = w + sw (10) L L tip where s is from 0 to L The eective wind momentum over tip mass is calculated around the beam r. free end. Table 2. Experimental coecients for dierent cylindrical cross sections. Isosceles 30 D-section Isosceles 53 Square a 2.9 0.79 1.9 2.3 a 6.2 0.19 6.7 18 2.3. Piezoelectric Modeling There are dierent ways of explaining the basic equations of the piezoelectric materials based on which variables are of interest. Matrix relationships are extensively used for ﬁnite element modelling. For analytical approaches, some of the relations are commonly valuable so the problem can be shortened. For example, matrix relationships describe strain in direction 3 as a function of stress and ﬁeld. To obtain the piezoelectric sensing voltage part of the formulation in the right-hand side of Equation (5) from Equation (2), the following is performed: First, the Gauss law over the piezoelectric sheets is considered: V(t) dQ = (11) R dt where R is a purely resistive load to harvest the produced electrical energy. In the current study, the internal resistance is neglected with respect to the load resistance. It is an appropriate AC circuit for initial evaluation of the electrical responses; therefore, it is repeatedly employed in the literature. The charge accumulated by the surface electrodes is found by surface integration (all over the two electrodes, i.e., 2w L ) from the electric displacement: p p Q = D.n dA (12) The vector of electric displacement appeared in Equation (12) in 3–1 mode is simpliﬁed to a scalar. The constitution equation of the piezoelectric is as follows: D = " E + E d " (13) 33 p 31 Acoustics 2019, 1 770 which shows the relation between the electric displacement vector and the electrical ﬁeld vector. In the current study, by use of the 3–1 mode, the matrix formulation is simpliﬁed to the scalar formulation. The vector of the electric ﬁeld is simpliﬁed to a scalar as follows: E = (14) Strain at the top surface of each piezoelectric sheet (electrode) is evaluated from the following: t + t p @ w " = ( ) (15) @x In some references, it is referred to as the membrane strain of piezoelectric materials and is the axial strain at the center of the piezoelectric material. By replacing deﬁnitions of strain at the top surface and electric ﬁeld through piezoelectric materials from Equations (14) and (15) in the deﬁnition of electric displacement (Equation (13)), the electric charge accumulated at the top surface of piezoelectric sheets covered by an electrode can be estimated. The time dierence of electric charges is used in the right-hand side of Equation (12) to ﬁnd the balance between the voltage of the electrode and beam displacement as follows: " # h i w L p p V(t) . 0 0 2" V(t) + + (L ) (0) E d w (t + t ) q (t) = 0 (16) p p p p 33 31 b i i i t R The ﬁrst term in the above equation is the electrical current through the two piezoelectric sheets w L p p with the capacitance of C = 2" ; the second term in the above equation is the electrical current p 33 through the impedance of electrical load connected to the piezoelectric sheets to harvest the energy; and the third term represents the eect of the rate of electrical charge accumulated on electrode surfaces. Detailed proof of Equation (16) is in Appendix A.4. 2.4. Approximate Method In this study, an approximate method is used to solve the problem: x (3L x) (x) = (L) (17) 2L The approximation used in Equation (17) is obtained from the static deﬂection of the clamped, free beam in the static condition under the unit tip force applied at the free end. It is good to notice that, in Reference [9], Sirohi and Mahadik used a single cubic shape function for beam deﬂection. Furthermore, the proposed function of Equation (17) is not completely satisfying Equation (3) as it is just an approximation to the exact solution. As shown in Figure 3, that function is in good agreement with the ﬁnal solution by the ﬁnite element method [13] with less complexity. Based on the single and parabola shape functions, which is assumed to represent the beam shape, its ﬁrst derivate versus location is derived as follows: 3x(2L x) (x) = (L) (18) 2L The eective values of system stiness, which is obtained from the potential energy (strain energy) of the device, is as follows: Z ! 1 @ w e f f P.E = EI dx = w (19) 2 2 @x 0 Acoustics 2019, 1 771 The eective mass, which is obtained from the kinetic energy of the system, is as follows: Z Z L L 1 . 1 . e f f 2 2 2 2 K.E = Aw dx + A y ds = w ! (20) tip tip tip 2 2 2 0 0 y can be evaluated from Equation (10). The explicit form of the eective mass is found from tip the following: h i 9I tip M = m [F (L) F (0)] + 2m F (L ) F (0) + m [1 + 2 ] + (21) p p e f f b M M M M tip 4L Detailed proof of Equation (21) is in Appendix A.2. The explicit form of the eective stiness is given by the following: h i K = k [F (L) F (0)] + 2k F (L ) F (0) (22) e f f b K K p K p K Detailed proof of Equation (22) is in Appendix A.1. The dimensionless function of mass eect versus location is as follows: Z ! 5 2 2 (x) x (5x 35Lx + 63L ) F (x) = d = (23) (L) 140L The dimensionless function of stiness eect versus location is as follows: h i (x) 3x 3 2 2 F (x) = L dx = x 3Lx + 3L (24) (L) L where the stiness constants of the components (beam and piezoelectric sheet) as a function of physical and geometrical characteristics are as follows: " # 3 3 E w t p p E w p p p t t b b k = + t + (25) 3 3 2 2 12L 3L b b E w t b b k = (26) 12L The mass constants of the components (tip, beam, and piezoelectric sheet) as a function of physical and geometrical characteristics are as follows: m = Adx = w t L (27) b b b b b m = A dx = w t L (28) p p p p p p p m L tip I = A s ds = (29) tip tip tip where the ratio is deﬁned by the following: L (L) 3L = = (30) (L) 2L Equation (21) contains all the dynamic mass eect caused by the tip mass and its mass moment of inertia, while the authors in Reference [9] (See Equation (26) in Reference [9].) missed the coupling of Acoustics 2019, 1 772 deﬂection and rotation term, i.e., m 2 . The ﬁrst term of Equation (21) presents the eect of the beam tip material. The second term of Equation (21) presents the eect of piezoelectric sheets on the eective mass. Here, they started from the clamped point until the length of L . The third term of Equation (21) presents the tip mass eect and the cross eect of rotation and displacement. The last term presents 9I tip w w L L the eect of rotation (i.e., K.E = = I ). The election of the dimensionless mode shape (in rot 2 tip 2 2 4L Equation (17)) will aect the dimensionless derivate of the mode shape (See Equation (18).), which is used to calculate piezoelectric coupling terms (See Equation (3).), and the eective mass and stiness distribution, which is used to evaluate the natural frequency of the system. Based on those functions which can be plotted as a function of the dimensionless variable (x/L), the calculations of the system can be performed. Figure 3 presents the dimensionless function for calculating eective mass, eective stiness, and the electromechanical coupling coecient. Figure 3 is a general ﬁgure that can be used to assess the composite beams where its elasticity and mass distributions are the functions of the position. As shown, as the location of the stiener element becomes closer to the clamped end, it is more eective (90 percent of the eectiveness of stiness is at the left half of the beam, and 50 percent of it is located at the ﬁrst 20 percent of the position) and, as the location of the mass element becomes closer to the free end, it is more eective (90 percent of the eectiveness of the mass is at the last 40 percent of the beam, and 50 percent of it is located at the last 20 percent of the position). For the cantilever beam with constant force at the free end considered here, the endpoint values can be calculated from Equation (18) which one can obtain the explicit expression as follows: (L) = (L) (31a) From Equations (23) and (24), they can be calculated as follows: F (L) = 0.2357 (31b) F (L) = 3 (31c) As mentioned before, the selected approximate solution (See Equation (17).) follows the case where it is the static solution of the cantilever beam with a concentrated tip force. This lumped parameter model has identical eciency with the distributed parameter model as it is degraded from the distributed parameters. It is good to notice that the fundamental frequency of the system could be found from the following: e f f ! = (32) e f f Although the natural frequency could be found from Equation (11) as (L) exists in both the nominator and denominator, it is needed for the next calculations. The coecient (L) is obtained by normalizing the eective mass to unity as follows: (L) = (33) e f f Acoustics 2019, 2 FOR PEER REVIEW 10 free end, it is more effective (90 percent of the effectiveness of the mass is at the last 40 percent of the beam, and 50 percent of it is located at the last 20 percent of the position). For the cantilever beam with constant force at the free end considered here, the endpoint values can be calculated from Equation (18) which one can obtain the explicit expression as follows: φφ (L) = (L) (31a) From Equations (23) and (24), they can be calculated as follows: Φ (L) = 0.2357 (31b) Φ= () L 3 (31c) As mentioned before, the selected approximate solution (See Equation (17).) follows the case where it is the static solution of the cantilever beam with a concentrated tip force. This lumped parameter model has identical efficiency with the distributed parameter model as it is degraded from the distributed parameters. It is good to notice that the fundamental frequency of the system could be found from the following: eff ω = (32) eff Although the natural frequency could be found from Equation (11) as exists in both the φ (L) nominator and denominator, it is needed for the next calculations. The coefficient is obtained φ (L) by normalizing the effective mass to unity as follows: φ (L) = (33) eff Acoustics 2019, 1 773 Figure 3. Dimensionless function for calculating of eective mass, eective stiness, Figure 3. Dimensionless function for calculating of effective mass, effective stiffness, and and electromechanical coupling coecient. electromechanical coupling coefficient. 2.5. Analytical Solution 2.5. Analytical Solution Substituting the approximate solution from Equation (17) in governing Equations (4) and (16) by Substituting the approximate solution from Equation (17) in governing Equations (4) and (16) considering only one mode, the equations of the system are simpliﬁed as follows: by considering only one mode, the equations of the system are simplified as follows: .. . L (6L3L ) p p w + 2!w + ! w = E d w (t + t ) V(t)+ L L L p 31 p p 3 b 2L M e f f . (34) . 3 2 4 Da w L V Da w L a 3 r a 1 L r 1 L 2 3 1 + + + 1 + 2 + 2 + + 2M 3 2M V 5 e f f e f f " # " # w L L (6L 3L ) p p V(t) p p . 2" V(t) + + E d w (t + t ) w (t) = 0 (35) 33 p 31 p p b L t R p 2L Detailed proof of Equation (34) is in Appendix A.3. The analytical solution of the above system of equations (i.e., Equations (34) and (35)) in steady-state operation is proposed as a trigonometric function of time with the same frequencies: x (3L x) w(x, t) = w sin(!t) (36) max 2L V(t) = V sin(!t + ') (37) max where t is the time, x is the location through the beam measured from the clamped end, L is the length of the beam, w is the maximum deﬂection of the beam which occurs at the free end, V is the max max maximum produced voltage by the piezoelectric sheets, ! is the angular velocity of the galloping motion, and ' is the phase dierence between the tip motion and voltage production. Replacing the assumptions of Equations (36) and (37) in Equation (35) and integrating both sides of the resulting equation from time instant 0 to time instant T/2, where the T is the period of motion, the relationship between w and V is obtained as follows: max max 3(2L L )L (t + t )w E d R! p p p b p p 31 V = r w (38) max max w L p p 2L 1 + 2" R! p Acoustics 2019, 1 774 Multiplying Equation (33) by w and Equation (34) by V, integrating the resulting equations concerning time from instant 0 to T/2 (time duration as half the period), and equating the coupling terms between two equations gives the following: u 2 u 3E d w (t +t )L (2LL ) p p p p p u 31 b u R u 2 3 V DL a u a r 2L 1 u 1 u 2!M + 1 + + u e f f 2 2 3 w L u p p t 1+ 2" R! w = (39) max DL a 3 a r 3 2 2 3 ! 1 + 2 + 2 + + 8 2V 5 It is good to notice that, in integrating from 0 to T/2, there is no variation in kinetic energy (ﬁrst term in Equation (33)), elastic energy (third term in Equation (33)), and the energy stored in capacitance (ﬁrst term in Equation (33)) since they will not appear in the ﬁnal equations (i.e., equation of energy balance). Also, the other coecient that aects the value open circuit obtained from Equation (38) in case of connecting to a load resistance is the phase dierence, which is deﬁned by a phase dierence of voltage and tip mass position as follows: 0 1 B C B C B C B C B C ' = arctan (40) B C B C w L p p @ A 2" R! Derivation of Equation (46) is straightforward from putting the values of voltage and displacement at the initial time (i.e., w = w !, V = V sin('), and V = !V cos(')) into Equation (34). max max max The eect of the phase dierence between voltage and displacement on the value of the load-connected circuit to the voltage of the open circuit is obtained from Equation (38) and is as follows: V V max max = cos(') (41) w w max max R!1 By applying the condition of a positive value under the square root, the onset of galloping (minimum value of the velocity where the expression under square root is positive) is derived from Equation (39) and is as follows: 3E d w (t +t )L (2LL ) p p p p p 31 b 2R 2L 4!M + e f f w L p p 1+ 2" R! V = (42) onset DL a 1 + + a r Based on the deﬁnition of the onset of galloping in Equation (42), the maximum tip deﬂection of the beam (See Equation (39).) is rearranged as a function of the velocity: u s u ! a 1 + + V V V onset 1 1 w = 1 (43) max ! V V 2 3 onset onset a 1 + 2 + 2 + + 4 5 Equation (43) in the current study is equal to Equation (18) in Reference [15] in the limit of a concentrated tip mass ( ! 0 ) and in considering the Taylor expansion of the arctangent function Acoustics 2019, 1 775 1 max where a = . The average harvested power over time (P = ) is obtained by integrating the ave 3 2R power during a period of motion: 3a 1 + + (2L L )L (t + t )w E d R p p p p p 3 V V b 31 2 1 1 P = ! V 1 (44) ave onset 2 4 V V w L p p 2 3 onset onset a 1 + 2 + 2 + + L 1 + 2" R! 2.6. A Note on Coupling Term Another value which is reported in the literature is the value : L (6L 3L ) p p = p E d w (t + t ) (45) p p p 31 b 2L M e f f where is the static coupling factor of the considered mode (i.e., k transversal coupling factor in the driving part). This parameter could appear as a simpler coecient of the voltage at the ﬁrst term in the right-hand side of Equation (34). Since the measured output of the piezoelectric sensing device is usually less than the expected ideal value, the above value needs modiﬁcation as follows: = (46) clamped theory piezo where is the performance of the piezoelectric at the clamped beam conﬁguration in comparison piezo with the ideal behavior. The performance of the piezoelectric at the clamped beam conﬁguration is E d 2 31 obtained in Reference [8] by use of the coecient k = . The modiﬁcation leads to the value of 0.2326 for the performance of the piezoelectric at the clamped beam conﬁguration in comparison with the ideal behavior (See Equation (28) in Reference [8].). If the correction of the value of the capacity of the piezoelectric sheet in comparison with the static piezoelectric sheet is considered, the performance of the piezoelectric at the clamped beam conﬁguration in comparison with the ideal behavior used in the coupling factor is obtained as follows: E d = 1 (47) piezo Based on that the deﬁnition of piezoelectric performance, Equation (4) for one mode and Equation (35) are restated as follows: .. . V Da qL a 1 r 2 1 q + 2!q + ! q = V(t) + 1 + + piezo 2M 3 e f f 3 (48) Da ' q L a r L 2 3 + 1 + 2 + 2 + + 2V 5 V . C V + + q = 0 (49) which is similar to the Equations (27) and (28) in Reference [14]. The only dierence with Reference [14] is the sign of , where in that reference, d is considered a positive value, while here, the negative sign is considered. The coupling factor has aected the performance of the system and is measured from the relation between displacement and produced voltage. The relation between the voltage and tip mass deﬂection in the open circuit case is proportional to the value. It can be found in Equation (38) that the voltage to displacement coecient of the system at open circuit condition is as follows: 3(2L L )(t + t )t E d V p p b p p 31 max = (50) 4L " max R!1 33 Acoustics 2019, 1 776 aR a where lim = . As shown in Equation (50) for the oered conﬁguration, the width of 2 2 b 1+b R R!1 the piezoelectric sheet does not aect the performance of the system since in the engineering design, the width of the piezoelectric sheet could be as thin as possible (regarding other engineering limitations) to save the amount of material. Also, for the arbitrary set of discrete piezoelectric sheets, the above formulation can be rephrased as follows: 6(t + t )t E d V p p p 31 max = (L x ) (51) ave w L " max R!1 33 where x is the average distance between the piezoelectric sheet (midpoint of the piezoelectric sheet) ave x +x e i from clamped point (x = ), x is the minimum distance between the piezoelectric sheet (left of ave i the piezoelectric sheet) from clamped point, and x is the maximum distance between the piezoelectric sheet (right of the piezoelectric sheet) from clamped point. Detailed proof of Equation (51) is in Appendix A.5. It is noticeable that Equation (50) is a speciﬁc case of Equation (51); by substitution of x = in Equation (51), Equation (50) is retrieved. Another point in Equation (50) is that the ave conﬁguration of two piezoelectrics cause the same voltage-displacement relation as that of the single piezoelectric case. The reason is that the factor 2 in the nominator of the two moments generated by the piezoelectric sheets is reduced by factor 2 in the denominator of the series of capacitance (twice capacitor value). The conclusion from Equation (50) is the evaluation of the theoretical maximum harvesting power in the energy harvesting of a cantilever with an attached prism under aeroelastic galloping. Replacing the generated voltage in the case of connecting load resistance (See Equation (41).) and maximum displacement of the case of no damping (internal structural damping and the electrical damping) in Equation (50), the maximum voltage in the energy harvesting case is obtained as follows: 2 r a 1+ + V 1 3(2LL )(t +t )t E d 3 V V p p b p p 31 onset 1 1 V = cos(') 1 (52) max 3 4 ! V V 4L " 33 3 2 3 onset onset a 1+2 +2 + + Optimal wind velocity could be obtained from Equation (49) as follows: opt V = 2V (53) 1 onset If the above velocity is replaced in Equation (49) and the onset velocity is approximated from Equation (42) by the assumption of negligible electrical damping, the following equation is derived for optimal voltage: u 2a 1 + + t 1 3(2L L )(t + t )t E d 3 4!M p p b p p 31 e f f V = cos(') (54) opt 3 4 2 4L " ! 33 2 3 a 1 + 2 + 2 + + DL a 1 + + 3 a r 4 5 3 At the phase dierence of ' = /4, the following value for optimal power is obtained: 0 1 B C B C B M C (2L L )(t + t )t E d p p p p e f f B b 31 C 1 B C B C P = 12 (55) opt B C B 3 4 C 2 B " C @ 2 3 A DL a 1 + 2 + 2 + + a 1 + + a r 3 1 5 3 3. Results and Discussion Various aspects of the energy harvester are analyzed next to assess the properties of dierent parameters on the level of harvested energy. As a ﬁrst stage, the numerical model is validated and compared with the experimental data performed by Sirohi and Mahadik [9]. Acoustics 2019, 1 777 3.1. Validation by Experimental Results Further, the evaluation of the damping ratio of the system is needed. The value = 0.003 is used in Reference [14] (See page 251, Equation (28).), but here, Figure 9 of Reference [9] is used in the simulations. Figure 4 presents the evaluation of the damping ratio of the system from the measured impulse response of the beam with electrodes in an open circuit [9]. The estimated value by using the Acoustics 2019, 2 FOR PEER REVIEW 16 tip values of voltage presented in Figure 4 are = 0.0109 and f = 4.17 (1/s). Figure 4 Figure 4. . Evaluation of Evaluation of the the damping ratio o damping ratio of f the the syst system em frfr om the om the measured impulse measured impulse response of response of th the e beam with electrodes in an open circuit [7]. beam with electrodes in an open circuit [7]. The most important features of the numerical results and its comparison with the experimental results of Sirohi and Mahadik [9] is presented in Table 3. The results of the simulation of the system (of which the physical property is given in Table 1) with the aerodynamic coecient given in Table 2 are summarized in Table 3. As shown in the calculation of natural frequency, the relative error of the used method is less than 2 percent while the relative error of the Galerkin procedure by Abdelkeﬁ et al. [14] is more than 45 percent (See Table 2 in Reference [14].). They found that the angular velocity for the ﬁrst mode is 38.176 which resulted in 6.076 Hz in frequency. Since the results of the approximate method have less error in comparison with the ﬁnite element method, this dierence comes from the application of the ﬁrst order Rayleigh–Ritz method instead of the higher-order modal method or Galerkin method used in Reference [14], which aect the dierence in natural frequency and eective mass of the system. Figure 5. Comparison of measured and predicted steady-state voltages as a function of incident wind velocity (R = 0.7 MΩ). Unfortunately, there is no published data other than References [9,16] for the problem of transverse galloping while the beam is normal to the wind direction and tip mass size in the direction of the beam (from clamped point to free end) is comparable to the beam length. More experiment results from Reference [17] are presented in Table 4 for validating the theoretical study and the numerical simulations. Those results are related to a tip mass of 23 g with length and width of 70 and 60 mm, respectively. Moreover, the range of wind velocity is 8–12 (m/s). Instead of a two-layer piezoelectric, Jamalabadi et al. [17] used one piezoelectric sheet partially covering the cantilever. The start point of the piezoelectric sheet was not the clamped point for some manufacturing considerations, and the piezoelectric sheet was located at the position with the highest strain values. Acoustics 2019, 1 778 Table 3. Numerical and experimental results of Sirohi and Mahadik [9]. Symbol Description and Unit Value m Beam mass (g) 5.9 m Piezoelectric sheet mass (g) 5.4 Acoustics 2019, 2 FOR PEER REVIEW 16 m Total mass (g) 654 k Beam stiness (N/m) 77.9 k Piezoelectric sheet stiness (N/m) 69.3 k Total stiness (N/m) 646.3 f Natural frequency (1/s) 5 f Experimental natural frequency (1/s) 4.167 exp C Total capacitance of the piezoelectric layers (nF) 658.7 1 1/2 Coupling factor (mmCm kg ) 12.8 m 1 Voltage to tip displacement of the system (Vmm ) 0.57 max R!1 R Load resistance (MW) 0.7 R Optimal resistance (kW) 57 max v Onset of the galloping velocity in open circuit Condition (ms ) 22.79 onset v Onset velocity in load resistance of 0.7 MW (ms ) 40.95 onset Also, it is revealed in Table 3 that the onset of galloping (with the aerodynamic coecients of Table 2) for the system based on the Equation (42) is predicted at 22.79 m/s. The value at load resistance (0.7 MW) is 40.9455 m/s, while for 9.5 mph (4.25 m/s), the value of 30 Volts (which is approximately a 2-mm tip deﬂection) is seen (See Figure 8 of Reference [9], and consider mph = 0.44704 m/s.). One conclusion of this paradigm could be that observed motion, in that case, is not a real galloping, and the coecients of Table 2 are not recommended. Although the authors of References [13,14] solved this paradigm by using artiﬁcial values for the system damping ratio, they exposed unrealistic results for tip displacement in the next step. The aerodynamic coecient values which ﬁt the experimental results of Reference [9] are a = 7.2 and a = 2.1757. Additionally, the ﬁtted results are present in 1 3 Figure 4. Evaluation of the damping ratio of the system from the measured impulse response of the Figure 5. The ﬁtted coecients predict the onset of galloping at 5.6 mph (2.5 m/s) for R = 0.7 MW, beam with electrodes in an open circuit [7]. and the maximum error is less than 5 Volts. Figure 5. Comparison of measured and predicted steady-state voltages as a function of incident wind Figure 5. Comparison of measured and predicted steady-state voltages as a function of incident wind velocity (R = 0.7 MW). velocity (R = 0.7 MΩ). Unfortunately, there is no published data other than References [9,16] for the problem of transverse Unfortunately, there is no published data other than References [9,16] for the problem of galloping while the beam is normal to the wind direction and tip mass size in the direction of the transverse galloping while the beam is normal to the wind direction and tip mass size in the direction of the beam (from clamped point to free end) is comparable to the beam length. More experiment results from Reference [17] are presented in Table 4 for validating the theoretical study and the numerical simulations. Those results are related to a tip mass of 23 g with length and width of 70 and 60 mm, respectively. Moreover, the range of wind velocity is 8–12 (m/s). Instead of a two-layer piezoelectric, Jamalabadi et al. [17] used one piezoelectric sheet partially covering the cantilever. The start point of the piezoelectric sheet was not the clamped point for some manufacturing considerations, and the piezoelectric sheet was located at the position with the highest strain values. Acoustics 2019, 1 779 beam (from clamped point to free end) is comparable to the beam length. More experiment results from Reference [17] are presented in Table 4 for validating the theoretical study and the numerical simulations. Those results are related to a tip mass of 23 g with length and width of 70 and 60 mm, respectively. Moreover, the range of wind velocity is 8–12 (m/s). Instead of a two-layer piezoelectric, Jamalabadi et al. [17] used one piezoelectric sheet partially covering the cantilever. The start point of the piezoelectric sheet was not the clamped point for some manufacturing considerations, and the piezoelectric sheet was located at the position with the highest strain values. Table 4. Geometrical and material properties [17]. Symbol Description and Unit Value Air density (kgm ) 1.225 m Tip mass (g) 31.6 L Length of the tip body (mm) 70 D Width of the tip body (mm) 60 V Wind velocity (m/s) 9–11 L Length of the beam (mm) 210 w Width of the beam material layer (mm) 50 t Beam material layer thickness (mm) 1 E 69 Beam material Young’s modulus (GNm ) Beam material density (kgm ) 3067 x Start of the piezoelectric sheets (mm) 5 L Length of the piezoelectric sheets (mm) 30 w Width of the piezoelectric layer (mm) 40 t Piezoelectric layer thickness (mm) 0.5 E 61 Piezoelectric material Young’s modulus (GNm ) Piezoelectric material density (kgm ) 8148 d 310 Strain coecient of the piezoelectric layer (pCN ) " Permittivity component at constant strain (nFm ) 38.3 Table 5 presents the numerical and experimental results of Jamalabadi et al. [17]. As shown, the natural frequency of the system is evaluated at 7.33 Hz while the experimental data (ﬁrst peak of frequency response) happens at 7.5 Hz (2 percent error). Table 5. Numerical and experimental results of Jamalabadi et al. [17]. Symbol Description and Unit Value m Beam mass (g) 32.2 m Piezoelectric sheet mass (g) 4.4 m Total mass (g) 57.6 k Beam stiness (N/m) 31.0 k Piezoelectric sheet stiness (N/m) 27.7 k Total stiness (N/m) 122.35 f Natural frequency (1/s) 7.33 f Experimental natural frequency (1/s) 7.5 exp C Total capacitance of the piezoelectric layers (nF) 92 1 1/2 Coupling factor (mmCm kg ) m 1 Voltage to tip displacement of the system (Vmm ) 0.467 max R!1 0.444 Voltage to tip displacement of the system (Vmm ) max exp v Onset of the galloping velocity in open circuit condition (ms ) 27.8 onset v 148.8 Onset velocity in load resistance of 0.7 MW (ms ) onset It is good to note that, in some references, to remove the dierence of the frequency gained by the numerical method and the measured frequency, a nonlinear torsional spring is considered at the clamped point for the eects of the incomplete clamp. In this study, the clamp point is considered a Acoustics 2019, 1 780 Acoustics 2019, 2 FOR PEER REVIEW 18 complete clamp condition. Also, the voltage to tip displacement of the system is evaluated with a good agreement (5 percent error), since the correction of Equation (47) is not necessary (if the correction of −1 Voltage to tip displacement of the system (Vmm) 0.444 Equation (47) is applied, the relative error increases to 12 percent). Onset of the galloping velocity in open circuit condition Figure 6 presents the variations of the voltage results of Reference [17] for dierent values of wind vonset 27.8 −1 velocity. Inspecting this ﬁgure, the theory used here ﬁts (ms the) experimental results for a = 2.416 and a = 1 3 −1 201.3095. vonset Onset velocity in load resistance of 0.7 MΩ (ms) 148.8 Figure 6. Comparison of measured and predicted steady-state voltages as a function of incident wind Figure 6. Comparison of measured and predicted steady-state voltages as a function of incident wind velocity at open circuit condition. velocity at open circuit condition. 3.2. Eects of the Load Resistance and Tip Mass Length Ratio on the Harvester s Response 3.2. Effects of the Load Resistance and Tip Mass Length Ratio on the Harvester′s Response The eect of tip mass length on eective mass is plotted in Figure 7. The curve is plotted by use of Equation The effect (21). oAs f tip ma shown, ss len if the gth on e dynamic ffective corr mas ection s is plot is not ted in applied Figu on re the 7. The c tip mass, urve iit s plot willtlead ed by tous ae higher of Equation error in (21). As sh the calculati owon n, if the dynami of eective mass, c correcti natural on i frsequency not applied on the ti , etc. After thepcorr maection ss, it wil oflFigur lead to e 8a in higher er calculating ror in the calculation of the eective mass, it ewould ffective mass be ready , nat toube ral used frequenc in design y, etc. After the calculations correction of provided for Fig the ure lumped 8 in calc method, ulating tas he effect presented ive mass in Refer , it w ence ould be r [7]. For eadexample, y to be used in the in desi casegn ca investigated lculations prov in thisided paper for , for the lumped the tip mass method, length to as p beam resented length in ratio Reference ( = ) [7]. of 2.61, For the example, contribution in the c of tip ase inve mass stig in the ated e in ective this mass is 13.95 times of its static mass. It means the total mass of the system is approximately 14 times paper, for the tip mass length to beam length ratio ( ) of 2.61, the contribution of tip mass in the γ = the static mass of the tip mass. If the formulation in Reference [9] is adopted, the value is 6.1 instead. Also, it is revealed in Figure 8 that, if the length of the rigid tip body is ten times that of the cantilever effective mass is 13.95 times of its static mass. It means the total mass of the system is approximately length, the eective mass should be considered in vibrational calculations and is one hundred times 14 times the static mass of the tip mass. If the formulation in Reference [9] is adopted, the value is 6.1 the mass of the tip body. instead. Also, it is revealed in Figure 8 that, if the length of the rigid tip body is ten times that of the cantilever length, the effective mass should be considered in vibrational calculations and is one hundred times the mass of the tip body. The theoretical relation between the open-circuit voltage and the tip displacement as predicted by Equation (43) is about 570.1643 volts per each meter deflection of the tip of the beam. These values will be corrected by the phase difference angle of 3.9477 degrees which leads to the correction of 568.8111 volts per meter. It means that the maximum voltage of 34 corresponds to the deflection of 6 centimeters of tip deflection. If the correction of Reference [9] is considered, this value decreases to Acoustics 2019, 2 FOR PEER REVIEW 20 Acoustics 2019, 2 FOR PEER REVIEW 20 could not exceed 100 milliwatts. The resultant harvested power of the current case with that of could not exceed 100 milliwatts. The resultant harvested power of the current case with that of References [13,14] based on the coefficient of Table 2 could be compared in Figure 10. References [13,14] based on the coefficient of Table 2 could be compared in Figure 10. In the current study, the analytical solution is compared with some previous cases and it is In the current study, the analytical solution is compared with some previous cases and it is expected that the model is useful for low-amplitude energy harvesting which is in the range of linear expected that the model is useful for low-amplitude energy harvesting which is in the range of linear elastic material strain. As shown, the experimental results demonstrate that the values in Table 2 are elastic material strain. As shown, the experimental results demonstrate that the values in Table 2 are not valuable in the calculation of harvested power for the low-amplitude galloping regime. The a1 not valuable in the calculation of harvested power for the low-amplitude galloping regime. The a1 value found in this study is considerably higher than the a1 value presented in Table 2. This value found in this study is considerably higher than the a1 value presented in Table 2. This estimates the onset of galloping velocity much lower than that in the experimental results (See estimates the onset of galloping velocity much lower than that in the experimental results (See Equation (42)). Also, the a3 coefficient is 1 order to 2 orders of magnitude higher than the a3 Equation (42)). Also, the a3 coefficient is 1 order to 2 orders of magnitude higher than the a3 presented in Table 2. Even the difference in estimation of onset of galloping velocity ignores the 1–2 presented in Table 2. Even the difference in estimation of onset of galloping velocity ignores the 1–2 order of magnitude difference in a3 and could cause a one order of magnitude difference in final order of magnitude difference in a3 and could cause a one order of magnitude difference in final deflection (See Equation (39).) and voltage (See Equation (38).) and two orders of magnitude deflection (See Equation (39).) and voltage (See Equation (38).) and two orders of magnitude difference in power results (See Equation (44), and compare Figure 10 with Figure 4 of Reference difference in power results (See Equation (44), and compare Figure 10 with Figure 4 of Reference [17]). The difference between fitted coefficient for Reference [9] (Re ≈ 10 ) and that of Reference [17] [17]). The difference between fitted coefficient for Reference [9] (Re ≈ 10 ) and that of Reference [17] (Re ≈ 10 ) could be attributed to the various characteristic lengths and ranges of velocity which lead (Re ≈ 10 ) could be attributed to the various characteristic lengths and ranges of velocity which lead to various Reynolds numbers. While the Reynolds number for Reference [9] is about 11,400 the to various Reynolds numbers. While the Reynolds number for Reference [9] is about 11,400 the Reynolds number for Reference [17] is about 67,260. It shows that, to globalize the results of the Reynolds number for Reference [17] is about 67,260. It shows that, to globalize the results of the Acoustics 2019, 1 781 aerodynamic coefficient to the problem of beam galloping perpendicular to the wind direction, more aerodynamic coefficient to the problem of beam galloping perpendicular to the wind direction, more experiments are needed. experiments are needed. Figure 7. Effect of tip mass length on effective mass. Figure Figure 7. 7. E Effe ect ct of t of tip ip m mass ass len length gth on effecti on eective ve m mass. ass. Figure 8. The ratio of the displacement of the distributed tip mass and the displacement of the point Figure 8. The ratio of the displacement of the distributed tip mass and the displacement of the point Figure 8. The ratio of the displacement of the distributed tip mass and the displacement of the point tip mass. tip mass. tip mass. The theoretical relation between the open-circuit voltage and the tip displacement as predicted by Equation (43) is about 570.1643 volts per each meter deﬂection of the tip of the beam. These values will be corrected by the phase dierence angle of 3.9477 degrees which leads to the correction of 568.8111 volts per meter. It means that the maximum voltage of 34 corresponds to the deﬂection of 6 centimeters of tip deﬂection. If the correction of Reference [9] is considered, this value decreases to 46 mm, while as the measurement of deﬂection is not reported in Reference [9], it could not be compared with experimental data. The added damping ratio by applying the electrical impedance is shown in Figure 9. The eective damping ratio as a function of load resistance is as follows: 3E d w (t +t )L (2LL ) p 31 p p p p 2L = ! (56) elec w L p p 2!M 1 + 2" R! e f f 33 p Acoustics 2019, 1 782 Acoustics 2019, 2 FOR PEER REVIEW 21 Figure Figure 9. 9. The The e equivalent quivalent dam damping ping ratio of ratio of the the ele electrical ctrical load load im impedance. pedance. Figure 9 presents the equivalent damping ratio of the electrical load impedance as a function of load impedance. As shown, the minimum added damping ratio to the system at a load resistance of 0.7 MW is 0.01. The artiﬁcial value of 0.003 for total damping ratio at that point could not be considered. As one can observe from Equation (39), the damping coecient should have a maximum value at a speciﬁc conﬁguration: R = (57) opt 2" w L ! 33 p p At this value of resistance, the phase dierence of voltage and displacement would be /4 and the electrical damping is maximized. Figure 9 presents the peak value when the value (R = 57 kW) is less than the initial range used in the experiment of Reference [8]. To compare the power harvested in Reference [13] with the current study, the electrical damping is deﬁned as follows: C = (58) piezo 1 + C R! where C is deﬁned as the electrical damping resulting from the electromechanical coupling by Tan and Yan under Equation (5) in Reference [13], (here is ) is adopted by comparing their Equation (1) p piezo Figure 10. Variation of the amplitude of the harvested power with the electrical damping C at w L p p with Equation (48), and C = 2" . By observation in Equation (56), it is clear that, if the electrical p 33 different wind speeds. p damping ratio deﬁned in Equation (56) is multiplied by angular velocity, the electrical damping deﬁned in Equation (58) is obtained (C = 2! ). If harvested power is considered during a period of motion 4. Conclusions V DL a In this study, galloping of two piezoelectric layers attached to a r a cantile ver beam with a tip mass 2( + )!M 1 + + e f f M 2 3 e f f exposed to the wind free stream was studied. The fluid–solid interaction Equations are solved P = (59) max 2 4 DL a !R a r 3 3 analytically by use of the approximation method. The result2ing model 3 is confirmed by preceding 1 + 2 + 2 + + 8 2V 5 experimental results. The effects of the tip mass length ratio on the onset of galloping, the level of the produced voltage, and harvested power are determined. In this research are the following: and rearranged by the deﬁnition of Equation (58), harvested power is M (4!+2C) e f f • In this study, a simple analytical model is used to encounter the effect of a tip mass which could V a 1 + + 4M CV DL 1 3 a r e f f be used by engineers for design of energy harvesting devices with piezoelectric materials. P = (60) max 2 2 4 3R ! 2 3 • The dimensionless functions for calculation of effective mass, effective stiffness, and a 1 + 2 + 2 + + electromechanical coupling coefficient were presented. • The effect of tip mass length on effective mass was developed. • The onset of galloping for the system was predicted analytically. • The current aerodynamic coefficient in the literature causes great numerical errors. Acoustics 2019, 2 FOR PEER REVIEW 21 Acoustics 2019, 1 783 As seen from the above equation, Equation (58) does not presents an explicit expression of harvested power as a function of C because, in Equation (58), R is still a function of C (See Equation (58).). The eects of electrical damping on the harvested power are presented in Figure 10. For the conﬁguration deﬁned in Table 1, the maximum value of C is 3.9722. Figure 8 reveals that even such odd wind velocities can be accessible for steady power harvesting; the amount of harvested power could not exceed 100 milliwatts. The resultant harvested power of the current case with that of References [13,14] based on the coecient of Figure 9. Table 2The e could qube ivalent dam compared ping in ratio of Figure the ele 10. ctrical load impedance. Figure 10. Variation of the amplitude of the harvested power with the electrical damping C at dierent Figure 10. Variation of the amplitude of the harvested power with the electrical damping C at wind speeds. different wind speeds. In the current study, the analytical solution is compared with some previous cases and it is 4. Conclusions expected that the model is useful for low-amplitude energy harvesting which is in the range of linear elastic material strain. As shown, the experimental results demonstrate that the values in Table 2 are not In this study, galloping of two piezoelectric layers attached to a cantilever beam with a tip mass valuable in the calculation of harvested power for the low-amplitude galloping regime. The a value exposed to the wind free stream was studied. The fluid–solid interaction Equations are solved found in this study is considerably higher than the a value presented in Table 2. This estimates the onset analytically by use of the approximation method. The resulting model is confirmed by preceding of galloping velocity much lower than that in the experimental results (See Equation (42)). Also, the a experimental results. The effects of the tip mass length ratio on the onset of galloping, the level of the coecient is 1 order to 2 orders of magnitude higher than the a presented in Table 2. Even the produced voltage, and harvested power are determined. In this research are the following: dierence in estimation of onset of galloping velocity ignores the 1–2 order of magnitude dierence in a and could cause a one order of magnitude dierence in ﬁnal deﬂection (See Equation (39).) and • In this study, a simple analytical model is used to encounter the effect of a tip mass which could voltage (See Equation (38).) and two orders of magnitude dierence in power results (See Equation (44), be used by engineers for design of energy harvesting devices with piezoelectric materials. and compare Figure 10 with Figure 4 of Reference [17]). The dierence between ﬁtted coecient • The dimensionless functions for calculation of effective mass, effective stiffness, and 4 5 for Reference [9] (Re 10 ) and that of Reference [17] (Re 10 ) could be attributed to the various electromechanical coupling coefficient were presented. characteristic lengths and ranges of velocity which lead to various Reynolds numbers. While the • The effect of tip mass length on effective mass was developed. Reynolds number for Reference [9] is about 11,400 the Reynolds number for Reference [17] is about • The onset of galloping for the system was predicted analytically. 67,260. It shows that, to globalize the results of the aerodynamic coecient to the problem of beam • The current aerodynamic coefficient in the literature causes great numerical errors. galloping perpendicular to the wind direction, more experiments are needed. 4. Conclusions In this study, galloping of two piezoelectric layers attached to a cantilever beam with a tip mass exposed to the wind free stream was studied. The ﬂuid–solid interaction Equations are solved analytically by use of the approximation method. The resulting model is conﬁrmed by preceding experimental results. The eects of the tip mass length ratio on the onset of galloping, the level of the produced voltage, and harvested power are determined. In this research are the following: Acoustics 2019, 1 784 In this study, a simple analytical model is used to encounter the eect of a tip mass which could be used by engineers for design of energy harvesting devices with piezoelectric materials. The dimensionless functions for calculation of eective mass, eective stiness, and electromechanical coupling coecient were presented. The eect of tip mass length on eective mass was developed. The onset of galloping for the system was predicted analytically. The current aerodynamic coecient in the literature causes great numerical errors. The equivalent damping ratio of the electrical load impedance as a function of load impedance was calculated. The eect of tip mass on the analytical solution of the concentrated mass was calculated. The eect of tip mass is a decrease in the tip displacement in comparison with the point mass. The results are ﬁtted on the analytical solution, and the new aerodynamics coecients included the eect of tip mass. As shown by increases of the length of tip mass for the constant beam and piezoelectric length, the inertia of the system increases while the tip displacement and onset of galloping decrease. Conﬂicts of Interest: The author declares no conﬂict of interest. Nomenclature Air density (kgm ) m Tip mass (g) L Length of the tip body (mm) D Width of the tip body (mm) V Wind velocity (m/s) L Length of the beam (mm) w Width of the beam material layer (mm) t Beam material layer thickness (mm) 0 2 E Beam material Young s modulus (GNm ) Beam material density (kg m ) x Start of the piezoelectric sheets (mm) L Length of the piezoelectric sheets (mm) w Width of the piezoelectric layer (mm) t Piezoelectric layer thickness (mm) E Piezoelectric material Young’s modulus (GN m ) Piezoelectric material density (kg m ) d Strain coecient of the piezoelectric layer (pC N ) " Permittivity component at constant strain (nF m ) m Beam mass (g) m Piezoelectric sheet mass (g) m Total mass (g) k Beam stiness (N/m) k Piezoelectric sheet stiness (N/m) k Total stiness (N/m) f Natural frequency (1/s) f Experimental natural frequency (1/s) exp C Total capacitance of the piezoelectric layers (nF) 1 1/2 Coupling factor (mmCm kg ) V Voltage of the piezoelectric layers (V) w Tip displacement of the beam (m) v Onset of the galloping velocity (ms ) onset Acoustics 2019, 1 785 Appendix A Appendix A.1. Find the Eective Stiness of the System As mentioned before, the eective values of the system’s stiness are obtained from the potential energy (strain energy) of the device (See Equation (19).). Assuming the approximate solution of Equation (17) for one mode, the second derivate of the beam displacement is as follows: @ w(x, t) L x = 3 q(t) (A1) 2 3 @x L where = (x = L). By substitution of the second derivate of the beam displacement into the potential energy (strain energy) of the system presented in Equation (19) and considering the deﬁnitions in Equations (24)–(26), the potential energy of the beam is obtained as follows: K K e f f e f f 2 P.E = w = ( q(t)) 2 2 R R R 2 2 2 2 L L L 1 2 1 @ w @ w = EI dx = E I dx + E I dx 2 p p 2 b b 2 0 2 0 @x 0 @x R R 2 2 L L 1 Lx Lx = E I 3 q(t) dx + E I 3 q(t) dx b b 3 L p p 3 L 2 0 L 0 L R R L E I L 2 E I 2 p p p 2 b b 9 9 = ( q) (L x) dx + (L x) dx L 3 3 3 3 2L 0 L L 0 L R R L L (A2) 2 k 2 p 2 b 9 9 = ( q) (L x) dx + k (L x) dx L 3 p 3 2 L L 0 0 R R L L 2 k 9 2 9 2 = ( q) (L x) dx + k (L x) dx L 3 3 0 L 0 L h h ii h h ii L L 2 k b 3x 2 2 3x 2 2 = ( q) x 3Lx + 3L + k x 3Lx + 3L L 3 p 3 L L 0 0 k [F (L)F (0)]+2k [F (L )F (0)] b K K p K p K 2 = ( q(t)) whereby comparing the ﬁrst and last rows, Equation (22) is obtained. Appendix A.2. Find the Eective Mass of the System As mentioned before, the eective values of the system’s mass are obtained from the kinetic energy of the device (See Equation (20).). Assuming the approximate solution of Equation (17) for one mode, the ﬁrst derivate of displacement of the beam versus time is obtained as follows: . . x (3L x) w(x, t) = q(t) (A3) 2L where = (x = L). Also, assuming the approximate solution of Equation (18) for one mode, the ﬁrst derivate of the slope of the beam versus time is obtained as follows: . . x (3L x) w(x, t) = q(t) (A4) 2L where = (x = L). By substitution of the above equations into the kinetic energy of the system presented in Equation (20) and considering the deﬁnitions in Equation (10) (for the y ), Equation (23) (for the dimensionless function of mass tip eect versus location), and Equations (27)–(30) (for mass constants of tip, beam, and piezoelectric sheet), the kinetic energy of the beam is obtained as follows: Acoustics 2019, 1 786 M . M . e f f 2 e f f K.E = w = q(t) 2 2 R R L . L . 2 r 2 1 1 = Aw dx + A y ds tip tip 2 2 tip 0 0 R R R . L . L . L . 2 2 p r 1 1 0 = A w dx + A w dx + A w + sw ds b b p p tip tip L L 2 2 0 0 0 R 2 R 2 2 2 . . L L x (3Lx) p x (3Lx) = A q(t) dx + A q(t) dx L 3 p p L 3 b b 2 0 2L 0 2L L . r 3 q(t) 1 L + A q(t) + s ds tip tip L 2 2L R 2 R 2 2 2 . 2 L x (3Lx) L x (3Lx) = q A dx + A dx L b b 3 p p 3 2 0 2L 0 2L 2 2 . L 1 3s + q A 1 + ds L tip tip 2 2L R 2 R 2 4 4 . L L x (3Lx) p x (3Lx) = q A dx + A dx L p p b b 6 6 0 4L 0 4L h i (A5) . 2 A 2 r tip tip 3s 9s + q 1 + + 2 L 4L R R 4 2 2 4 2 2 . 2 L x (9L 6Lx+x ) L x (9L 6Lx+x ) = q A dx + A dx L b b 6 p p 6 2 0 4L 0 4L 3L 3L 1 r r + q A L + + L tip tip r 2 2L 8L R R . L 2 4 5 6 L 2 4 5 6 A L p b b 9L x 6Lx +x 9L x 6Lx +x = q dx + A L dx p p L 7 7 0 4L 0 4L h i . 2 9I 3L tip 1 r + q A L + + L tip tip r 2 2 2L 8L L L 5 2 2 5 2 2 p m x (5x 35Lx+63L ) x (5x 35Lx+63L ) = q + m L 7 7 140L 140L 0 0 . 2 m 9I tip tip + q [1 + 2 ] + L 2 2 8L 9I tip m [F (L)F (0)]+2m F (L )F (0) +m [1+2 ]+ 2 b M M p[ M p M ] tip . 4L = q whereby comparing the ﬁrst and last rows of Equation (A5), the Equation (21) is obtained. Appendix A.3. Find the Second-Order Governing Equation of Motion of the System Neglecting the dynamic eects of the longitudinal forces and rotational inertia and considering the governing equation of transversal vibration of a beam as presented by Equation (2), one can use the equivalent method to solve the equation. If both sides of Equation (2) are multiplied by the essential mode shape and integrated all over the system, one can obtain for the essential mode the following: L+L L L R R R .. . 2 iv iv A qdx + EI qdx + cI qdx = 0 0 0 h i 2 (A6) E d w (t + t )V (x x ) dx+ p p p 31 b 1 dx R R L L M (x L)dx + F (x L)dx tip tip 0 dx 0 where the ﬁrst term on the left-hand side is the inertial term of which the eective mass is found in Equation (21) by use of the method of energy ! ! Z Z 2 0 2 L L ' ' M = A dx + A 1 + s ds (A7) e f f tip tip ' ' 0 L 0 L and where the second term is the eective stiness which is found in Equation (19) by use of the method of energy 00 2 K = EI' dx (A8) e f f while the integration is converted based on the following partial integral relation, R R L L iv 000 0 000 '' dx = '' ] ' ' dx 0 0 R (A9) L L 000 0 00 00 2 = '' ] ' ' ] + ' dx 0 0 0 Acoustics 2019, 1 787 on the clamped boundary condition on the left (x = 0), and the free end on the right (x = L). Since the left-hand side of Equation (2) is rewritten in the form of the second-order dierential equation of a single degree of freedom motion of a point mass system with damping, the right-hand side is presented by a generic outside force as follows: .. . M q + C q + K q = e f f e f f e f f .. . (A10) M q + 2!q + ! q = f e f f e f f Also, if both sides of Equation (3) are multiplied by the essential mode shape and integrated all over the system. one can obtain the Eigen-function equation for the essential mode: K = M ! (A11) e f f e f f which is equal to Equation (32). Similar to the above derivations, which are for one mode condition, we have the following: .. . q + 2!q + ! = 0 0 (L ) (0) E d w (t + t )V p p p p 31 b (A12) h i + M (x L) + F (x L) (x)dx tip tip dx The right-hand side of Equation (A12) is simpliﬁed for the force term as follows: h i F (x L) (x)dx = F (A13) tip tip L and for the momentum term as follows: Z " # M (x L) (x)dx = M (A14) tip tip dx Equation (A12) is rewritten as .. . q + 2!q + ! = 0 0 (L ) (0) E d w (t + t )V (A15) p p 31 p p +F + M tip L tip The eective tip force is obtained by integration of Equation (9) over the tip mass length. Through that integration, the deﬁnition of Equation (10) (for the y ) should be considered. F is found from the following: tip tip Acoustics 2019, 1 788 " # r . . R 3 y y 1 2 F = V D a + a ds tip a 1 3 2 1 V V 1 1 " # L . . r 3 R . . 0 0 w +sw w +sw 1 2 L L L L = V D a + a ds a 1 3 2 1 V V 1 1 " # L . . . . R 0 0 3 q+s q q+s q L L 1 2 L L L L = V D a + a ds a 1 3 2 1 V V 1 1 R h i . . 1 0 = V Da q + s q ds a 1 L L L 2 1 . . 3 Da a 3 + q + s q ds L L 2V R h i h i . . 1 0 = V Da q + s q ds a 1 L L L (A16) 2 1 . 3 . . 2 . . 2 . 3 Da a 3 0 2 0 3 0 + q + 3s q q + 3s q q + s q ds L L L L L L 2V . L . = V DL a q + q a r L 2 1 2 . 3 . . 2 . . 2 . 3 DL a L a r 3 r L L 2 0 r 0 + q + 3 q q + L q q + q L L L 2V 2 L 4 L 1 L = V DL a q 1 + a r 1 L 2 1 2 2 3 0 0 0 DL a q a r 3( L ) L L L r r r L L 1 L + 1 + 3 + + 2V 4 1 L L L = V DL a q 1 + a r 1 L 2 1 2 . 3 DL a ( q) a r 3 L + 1 + 3 + + 2V 2 4 Morover, the eective tip moment is obtained by integration of Equation (6) where the force term is presented by Equation (9). Over and done with that integration, the deﬁnition of Equation (10) (for the y ) should be tip considered. M is found from tip " # r . . R 3 y y 1 2 M = V D a + a sds tip a 1 3 2 1 V V 1 1 " # L . . r . . 3 0 0 w +sw w +sw 1 2 L L L L = V D a + a sds a 1 3 2 1 V V 1 1 " # L . . . . r 3 R 0 0 q+s q q+s q L L 1 2 L L L L = V D a + a sds a 1 3 2 1 V V 1 1 R h i . . 1 0 = V Da q + s q sds a 1 L 2 L L . . 3 Da a 3 + q + s q sds 2V L L R h i h i . . 1 2 0 = V Da s q + s q ds a 1 L 2 L L 1 (A17) . 3 . . 2 . . 2 . 3 Da a 3 2 0 3 0 4 0 + q s + 3s q q + 3s q q + s q ds L L L 2V L L L q L q L r 1 2 L = V DL a + a r 1 2 1 2 3 . 3 q . . 2 . . 2 3 . 3 DL a ( L ) a r 3 L 0 3 2 0 r 0 + + L q q + L q q + q r L r L L L L 2V 2 4 5 1 1 r 2 L = V DL a q + a r 1 L 2 1 2 3 . 3 2 2 3 0 0 0 DL a ( q) L L L a r 3 L r r r 1 3 1 L L L + + + + 2V 2 4 5 1 L L L 1 2 1 = V DL a q + a r 1 L 2 1 2 3 . 3 DL a ( q) a r 3 L 1 3 2 1 3 + + + + 2V 2 4 5 1 Acoustics 2019, 1 789 The eective force applied by the tip mass on the right-hand side of Equation (16) is as follows: " # . 3 . DL a q a r 3( L ) 0 1 2 F + M = V DL a q 1 + + 1 + 3 + + tip L tip a r 1 L L L 2 2 2V 2 1 4 " # 2 3 DL a q a r 3 L 1 2 1 1 3 2 1 3 0 + V DL a q + + + + + a r 1 L 2 1 2 3 2V 2 4 5 L h i 1 1 0 = V DL a q 1 + + L + a r 1 L L r 2 1 2 2 3 . 3 DL a ( q) a r 3 L 1 3 1 2 2 3 0 + 1 + 3 + + + L + + + 2V 2 4 2 4 5 L 1 2 L 1 = V DL a q 1 + + + a r 1 L 2 2 2 3 L (A18) . 3 4 3 0 DL a q r a r 3 L 2 L 1 3 2 1 3 + 1 + 3 + + + + + + 2V 2 4 2 4 5 1 L h i 1 1 = V DL a q 1 + + + a r 1 L 2 1 2 2 3 4 3 DL a q a r 3 L 1 3 1 2 2 3 + 1 + 3 + + + + + + 2V 2 4 2 4 5 1 2 = V DL a q 1 + + a r 1 L 2 3 . 3 4 4 DL a q a r 3 L 2 3 + 1 + 2 + 2 + + 2V 5 Put the eective values of the above equations together h i .. . 2 2 0 0 M q + 2!q + ! q = (L ) (0) E d w (t + t )V(t)+ e f f L L p p 31 p p b . . 2 4 4 (A19) V Da q L Da q L a 1 L r a 3 L r 1 2 3 1 + + + 1 + 2 + 2 + + 2 3 2V 5 by normalizing the modes with the eective mass, = (A20) e f f then, the ﬁnal equation is obtained .. . q + 2!q + ! q = h i 0 0 (L ) (0) E d w (t + t )V(t)+ p p 31 p p b V Da q L a 1 L r (A21) 1 + + + 2 3 4 4 Da q L a 3 L r 2 3 1 + 2 + 2 + + 2V 5 by dividing both sides of Equation (A21) by ' 0 0 .. . (L ) (0) M q + 2! q + ! q = E d w (t + t )V(t)+ L L L L p 31 p p e f f b . . (A22) 2 4 V Da qL Da q L a 1 L r a 3 L r 1 2 3 1 + + + 1 + 2 + 2 + + 2 3 2V 5 Recall the normalizing condition (Equation (33)); then, Equation (A22) is rewritten as follows: .. . L (6L3L ) p p M w + 2!w + ! w = E d w (t + t ) V(t)+ p p p e f f L L L 3 31 b 2L . 3 (A23) 2 4 Da w L V Da w L a 3 r a 1 L r 1 L 2 3 1 + + + 1 + 2 + 2 + + 2 3 2V 5 Which, by dividing both sides with M , Equation (34) is retrieved. e Acoustics 2019, 2 FOR PEER REVIEW 28 LL (6L − 3 ) pp M ww++2( ξω ωw = Edwt+t)V(t)+ () eff LL L 3 p 31pp b 2L (A23) ρρ VDawLL Da w γγ a 1rLa 3 r ∞ L 23 1++γγ + 1+22 +γ +γ + 22 3 V 5 Which, by dividing both sides with Meff, Equation (34) is retrieved. Appendix A4. Find the Piezoelectric Effect on the Governing Equation of Motion of the System To prove the piezoelectric coupling terms in the Euler–Bernoulli beam equation of the system, Acoustics 2019, 1 790 the constitution equation of strain-electrical displacement should be considered as follows: Appendix A.4. Find the Piezoelectric Eect on the Governing Equation of Motion of the System ε=+ dE 31 (A24) To prove the piezoelectric coupling terms in the Euler–Bernoulli beam equation of the system, the constitution equation of strain-electrical displacement should be considered as follows: where σ is the stress and E is the electrical field. For the pin force theory, the piezoelectric field is modeled as follows: " = d E + (A24) FA== σσwt (A25) piezo p p p where is the stress and E is the electrical ﬁeld. For the pin force theory, the piezoelectric ﬁeld is modeled as follows: Based on the dimension of the piezoelectric, the electric field is found as defined by Equation F = A = w t (A25) p p p piezo (14). For the no-strain condition assumed in the pin force model (balance of electrical and mechanical strain Based ), the fo on the llow dimension ing equat of ion i the s piezoelectric, obtained: the electric ﬁeld is found as deﬁned by Equation (14). For the no-strain condition assumed in the pin force model (balance of electrical and mechanical strain), the following equation is obtained: F =Ed w t (A26) piezo p 31 p p F = E d w t (A26) p p p piezo 31 pt The piezoelectric sheet eects as a moment at the ending positions of the sheet are evaluated based on the The piezoelectric sheet effects as a moment at the ending positions of the sheet are evaluated modiﬁed pin force model as presented in Figure A1. based on the modified pin force model as presented in Figure A1. Figure A1. Front view of the problem. Figure A1. Front view of the problem. As shown, the value of each moment (top and bottom) is produced by a constant force located at the half-height of the piezoelectric sheet as follows: As shown, the value of each moment (top and bottom) is produced by a constant force located at the half-height of the piezoelectric sheet as follows: t + t M = F (A27) piezo piezo tt + p b MF = (A27) piezo piezo The eect of the moment of two piezoelectrics in the governing equation of the beam is as follows: The effect of the moment of two piezoelectrics in the governing equation of the beam is as M = E d w (t + t )V(t) (A28) p p p piezo 31 b follows: M =+ Ed w() t t V(t) (A28) piezo p 31 p p b Also, the derivation of the piezoelectric charge equation is clearer in the following: Acoustics 2019, 1 791 Also, the derivation of the piezoelectric charge equation is clearer in the following: V(t) dQ R dt = D.n dA dt = " E + E d " dA 33 p 31 dt t +t 2 Lp b d V @ w = 2 " + E d ( ) w dx 33 p 31 2 p t 2 dt 0 p @x t +t Lp b 2 d V @ w = 2 " + E d ( ) w dx 33 p 31 2 p dt 0 t 2 @x R R Lp 2" Vw Lp 2 33 p d d @ w = dx + E d (t + t ) w dx p p (A29) 31 p b 2 dt t dt 0 p 0 @x h i Lp Lp 2" Vw 33 p d d @w = dx + E d (t + t ) p 31 p b dt 0 p dt @x R h i Lp 2" Vw 33 p d d 0 0 = (L ) (0) q (t) dx + E d (t + t ) p 31 b p i p i i dt 0 t dt h i w L p p d d 0 0 = 2" (V) + E d (t + t ) (L ) (0) q (t) 33 p 31 p b p i t i i p dt dt h i w L . p p 0 0 = 2" V + E d (t + t ) (L ) (0) q (t) p p 33 31 p b t i i i which is equal to Equation (16). Appendix A.5. Find the Relation of Piezoelectric Voltage and Tip Displacement For the arbitrary set of discrete piezoelectric sheets, the above formulation rewritten as follows: 3(t + t )t E d V p p p b 31 max = (x (2L x ) x (2L x )) (A30) e e i i w L L " max R!1 p 33 while the slope dierence is obtained from Equation (18). It rearranged as follows: 3(t +t )t E d p p p 31 V b max 2 2 = 2L(x x ) (x x ) e e 3 i i max L L " R!1 p 33 3(t +t )t E d p b p p 31 = (2L (x + x )) (x x ) (A31) 3 e i e i L L " p 33 6(t +t )t E d p b p p 31 x +x e i = L ( ) (x x ) 3 e i L L " 2 p 33 Substituting L = x x in Equation (A31), it reduces to Equation (51). p e i References 1. Saadon, S.; Sidek, O. A review of vibration-based MEMS piezoelectric energy harvesters. Energy Convers. Manag. 2011, 52, 500–504. [CrossRef] 2. Erturk, A.; Inman, D.J. An experimentally validated bimorph cantilever model for piezoelectric energy harvesting from base excitations. Smart Mater. Struct. 2009, 18, 025009. [CrossRef] 3. Ilyas, M.A.; Swingler, J. Piezoelectric energy harvesting from raindrop impacts. Energy 2015, 90, 796–806. [CrossRef] 4. Abdollahzadeh Jamalabadi, M.Y.; Ahmadi, E. Design, Construction and Testing of a Dragon Wave Energy Converter. Am. J. Nav. Archit. Mar. Eng. 2016, 1, 7–15. 5. Abdollahzadeh Jamalabadi, M.Y.; Ahmadi, E. Experimental study on overtopping performance of a circular ramp wave energy converter. Rev. Energy Technol. Policy Res. 2017, 3, 1–16. 6. Tan, T.; Yan, Z. Electromechanical decoupled model for cantilever-beam piezoelectric energy harvesters with inductive-resistive circuits and its application in galloping mode. Smart Mater. Struct. 2017, 26, 035062. [CrossRef] 7. Barrero-Gil, A.; Alonso, G.; Sanz-Andres, A. Energy harvesting from transverse galloping. J. Sound Vib. 2010, 329, 2873–2883. [CrossRef] 8. Sirohi, J.; Mahadik, R. Piezoelectric wind energy harvester for low-power sensors. J. Intell. Mater. Syst. Struct. 2011, 22, 2215–2228. [CrossRef] 9. Sirohi, J.; Mahadik, R. Harvesting wind energy using a galloping piezoelectric beam. J. Vib. Acoust. 2012, 134, 011009. [CrossRef] Acoustics 2019, 1 792 10. Abdelkeﬁ, A.; Yan, Z.; Hajj, M. Modeling and nonlinear analysis of piezoelectric energy harvesting from transverse galloping. Smart Mater. Struct. 2013, 22, 025016. [CrossRef] 11. Bibo, A.; Daqaq, M.F. On the optimal performance and universal design curves of galloping energy harvesters. Appl. Phys. Lett. 2014, 104, 023901. [CrossRef] 12. Tan, T.; Yan, Z.; Hajj, M. Electromechanical decoupled model for cantilever-beam piezoelectric energy harvesters. Appl. Phys. Lett. 2016, 109, 101908. [CrossRef] 13. Tan, T.; Yan, Z. Analytical solution and optimal design for galloping-based piezoelectric energy harvesters. Appl. Phys. Lett. 2016, 109, 253902. [CrossRef] 14. Abdelkeﬁ, A.; Yan, Z.; Hajj, M.R. Performance analysis of galloping-based piezoaeroelastic energy harvesters with dierent cross-section geometries. J. Intell. Mater. Syst. Struct. 2014, 25, 246–256. [CrossRef] 15. Kazakevich, M.I.; Vasilenko, A.G. Closed Analytical Solution for Galloping Aeroelastic Self-Oscillations. J. Wind Eng. Ind. Aerodyn. 1996, 65, 353–360. [CrossRef] 16. Sirohi, J.; Mahadik, R. Harvesting Wind Energy Using a Galloping Piezoelectric Beam. In Proceedings of the ASME 2009 Conference on Smart Materials, Adaptive Structures and Intelligent Systems, Oxnard, CA, USA, 21–23 September 2009; pp. 443–450. 17. Abdollahzadeh Jamalabadi, M.Y.; Kwak, K.M.; Hwan, S.J. Galloping modeling of a cantilever beam with piezoelectric wafers. In Proceedings of the Joint Conference by KSNVE, ASK and KSME(DC), Gwangju, Korea, 26–18 April 2017; Volume 4, p. 478. 18. Abdollahzadeh Jamalabadi, M.Y.; Kwak, K.M.; Hwan, S.J. Flow-induced vibration of Piezoelectric D-shape Beam Wind Energy Harvesters. Proc. KSNVE Annu. Autumn Conf. 2016, 10, 54. 19. Kim, M.; Hoegen, M.; Dugundji, J.; Wardle, B.L. Modeling and experimental veriﬁcation of proof mass eects on vibration energy harvester performance. Smart Mater. Struct. 2010, 19, 045023. [CrossRef] 20. Yang, Y.; Zhao, L.; Tang, L. Comparative study of tip cross-sections for ecient galloping energy harvesting. Appl. Phys. Lett. 2013, 102, 064105. [CrossRef] 21. Bearman, P.W.; Gartshore, I.S.; Maull, D.J.; Parkinson, G.V. Experiments on ﬂow-induced vibration of a square-section cylinder. J. Fluids Struct. 1987, 1, 19–34. [CrossRef] 22. Corless, R.M.; Parkinson, G.V. A model of the combined eects of vortex-induced oscillation and galloping. J. Fluids Struct. 1988, 3, 203–220. [CrossRef] 23. Corless, R.M.; Parkinson, G.V. Mathematical modelling of the combined eects of vortex-induced vibration and galloping. Part II. In Proceedings International Symposium Flow-Induced Vibration and Noise, Blu-Body/Fluid and Hydraulic Machine Interactions; Pa¨ıdoussis, M.P., Dalton, C., Weaver, D.S., Eds.; ASME: New York, NY, USA, 1992; Volume 6, pp. 39–63. 24. Corless, R.M.; Parkinson, G.V. Mathematical modelling of the combined eects of vortex-induced vibration and galloping. Part II. J. Fluids Struct. 1993, 7, 825–848. [CrossRef] 25. Laneville, A.; Parkinson, G.V. Eects of Turbulence on Galloping of Blu Cylinders. Ph.D Thesis, University of British Columbia, Vancouver, BC, Canada. 26. Parkinson, G.V. Wind-induced instability of structures. Philos. Trans. R. Soc. 1971, 269, 395–409. [CrossRef] 27. Parkinson, G.V. Mathematical models of ﬂow-induced vibrations of blu bodies. In Flow-Induced Structural Vibrations; Naudascher, E., Ed.; Springer: Berlin, Germany, 1974; pp. 81–127. 28. Parkinson, G.V. Phenomena and modelling of ﬂow-induced vibrations of blu bodies. Prog. Aerosp. Sci. 1989, 26, 169–224. [CrossRef] 29. Parkinson, G.V.; Brooks, N.P.H. On the aeroelastic instability of blu cylinders. J. Appl. Mech. 1961, 28, 252–258. [CrossRef] 30. Parkinson, G.V.; Jandali, T. A wake source model for blu body potential ﬂow. J. Fluid Mech. 1969, 40, 577–596. [CrossRef] 31. Parkinson, G.V.; Smith, J.D. The square prism as an aeroelastic non-linear oscillator. Q. J. Mech. Appl. Math. 1964, 17, 225–239. [CrossRef] 32. Parkinson, G.V.; Sullivan, P.P. Galloping response of towers. J. Ind. Aerodyn. 1979, 4, 253–260. [CrossRef] 33. Parkinson, G.V.; Wawzonek, M.A. Some considerations of the combined eects of galloping and vortex resonance. J. Wind Eng. Ind. Aerodyn. 1981, 8, 135–143. [CrossRef] 34. Den Hartog, J.P. Transmission line vibration due to sleet. Trans. Am. Soc. Electr. Eng. 1932, 51, 1075–1076. [CrossRef] 35. Novak, M. Aeroelastic galloping of prismatic bodies. ASCE J. Eng. Mech. 1969, 95, 115–142. Acoustics 2019, 1 793 36. Abdollahzadeh Jamalabadi, M.Y.; Kwak, M.K. Dynamic modeling of a galloping structure equipped with piezoelectric wafers and energy harvesting. Noise Control Eng. J. 2019, 67, 142–154. [CrossRef] 37. Abdelkeﬁ, A.; Scanlon, J.M.; McDowell, E.; Hajj, M.R. Performance enhancement of piezoelectric energy harvesters from wake galloping. Appl. Phys. Lett. 2013, 103, 033903. [CrossRef] 38. Ewere, F.; Wang, G. Performance of galloping piezoelectric energy harvesters. J. Intell. Mater. Syst. Struct. 2014, 25, 1693–1704. [CrossRef] 39. Rostami, A.B.; Fernandes, A.C. The eect of inertia and ﬂap on autorotation applied for hydrokinetic energy harvesting. Appl. Energy 2015, 143, 312–323. [CrossRef] © 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Acoustics – Multidisciplinary Digital Publishing Institute

**Published: ** Sep 24, 2019

**Keywords: **Wind energy; sustainable energy; galloping; piezoelectric; energy harvesting

Loading...

You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!

Read and print from thousands of top scholarly journals.

System error. Please try again!

Already have an account? Log in

Bookmark this article. You can see your Bookmarks on your DeepDyve Library.

To save an article, **log in** first, or **sign up** for a DeepDyve account if you don’t already have one.

Copy and paste the desired citation format or use the link below to download a file formatted for EndNote

Access the full text.

Sign up today, get DeepDyve free for 14 days.

All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.