Modeling the coupled electro-mechanical response of a torsional-flutter-based wind harvester with a focus on energy efficiency examination

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

Contents lists available at ScienceDirect

Journal of Wind Engineering & Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia

Modeling the coupled electro-mechanical response of a torsional-flutter-based wind harvester with a focus on energy efficiency examination Luca Caracoglia Department of Civil and Environmental Engineering, Northeastern University, 400 Snell Engineering Center, 360 Huntington Avenue, Boston, MA 02115, USA

A R T I C L E I N F O

A B S T R A C T

Keywords: Wind-based energy harvesting Aeroelastic vibration Torsional flutter Post-critical dynamics Electro-mechanical coupling

Wind energy is a rapidly evolving field of research because of the need for clean energy resources. Large horizontal axis wind turbines (HAWT) are often employed to increase output power and energy production. On the other hand, “specialized” wind-based energy systems have been proposed to capture the wind energy resource in the low wind speed range and for intermediate-scale applications, e.g., one or few residential housing units. Wind harvesters, triggered by various aeroelastic instability regimes, have been studied in recent years [e.g., Matsumoto et al. (2006)]. Along this line of research, the writer has examined a torsional-flutter-based apparatus for extracting energy from the wind flow. This paper presents some recent advancements, a new fully-coupled electro-mechanical model and the numerical results of an ongoing investigation.

1. Introduction The exploitation of aeroelastic phenomena for energy production is a relatively new technological area in comparison with other research activities in the field of aeroelasticity. Nevertheless, several new contributions and studies have appeared in the literature. One of the first interesting examples refers to the concept of the flutter mill (Tang et al., 2009) for extracting energy from the wind. This conceptual idea has captivated the interest of the research community in the last decades. Various solutions have been proposed and examined by researchers (Ahmadi, 1978; Farthing, 2009; Kwon et al., 2011; Matsumoto et al., 2006; Shimizu et al., 2008; Tang et al., 2009). Among the most interesting examples, the flapping wing power generator (Shimizu et al., 2008) is based on the principle of the classical flutter of airfoils, involving a combined pitching-heaving motion. Another apparatus, which exploits the two-degree-of-freedom coupled flutter of a streamlined body in plunging and pitching motion, has been investigated both analytically and experimentally for energy extraction (Zhu, 2011). The coupled-mode flutter-mill concept has also been proposed and used for water pumping and irrigation purposes (Farthing, 2009). The configuration of this kind of aeroelastic harvesters is composed of either a rigid blade, mounted on a mobile support (Abdelkefi et al., 2012a) or several flexible blades (Abdelkefi et al., 2012b; Dunnmon et al., 2011; Kwon et al., 2011). If the blade or airfoil is

exposed to an air flow speed beyond the critical flutter speed, large vibration can be triggered. The resulting limit-cycle oscillation regime can be exploited and converted to electrical energy (Dowell et al., 2004; Dunnmon et al., 2011). The limit-cycle is influenced by nonlinear aeroelastic effects. Oscillation amplitudes must be controlled both to optimize energy extraction from the flow and facilitate conversion to electrical power. Along the same line, the use of airfoil-like plates equipped with porous screens has been recently proposed as a new development of the flutter mill concept (Pigolotti et al., 2016, 2017). In most successful applications of the flutter mill concept (Abdelkefi et al., 2012a; Dunnmon et al., 2011) the use of piezo-electric materials has been indicated as a practical technology that converts elastic deformation energy to electricity. Piezo-electric actuators have also been recently used to harvest energy from transverse galloping of square prisms in fluid flow (Abdelkefi et al., 2013, 2014), a mechanism conceptually simpler than coupled flutter, and from flutter-induced vibration of miniature beams, predominantly employed for sensor design (Casadei and Bertoldi, 2014). Capitalizing from some recent advances in the exploitation of piezo-electric materials, a number of design configurations of galloping-based harvesters have also been recently examined and experimentally tested, either at small scale (Tomasini and Giappino, 2016) or miniature scale [e.g., for installations inside air ventilation systems (Biscarini et al., 2016; Gkoumas et al., 2017; Petrini et al., 2014)]. The piezo-electric technology is typically efficient if either

E-mail address: [email protected]. https://doi.org/10.1016/j.jweia.2017.10.017 Received 16 February 2017; Received in revised form 15 October 2017; Accepted 16 October 2017 Available online xxxx 0167-6105/© 2017 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Caracoglia, L., Modeling the coupled electro-mechanical response of a torsional-flutter-based wind harvester with a focus on energy efficiency examination, Journal of Wind Engineering & Industrial Aerodynamics (2017), https://doi.org/10.1016/ j.jweia.2017.10.017

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

Nomenclature qem ðwem Þ R1 , R2 RC s t U U* wem α αI αmax

The following symbols are used in this study B average magnetic flux density of the winding coil [teslas] b half-width (half-chord) of the blade-airfoil harvester [m] CðkÞ complex Theodorsen function CP dimensionless power coefficient of the harvester cj ; dj parameters of the Wagner function with j ¼ f1; 2g e rotational-axis eccentricity, measured from structural support (Fig. 1) FðkÞ real part of Theodorsen function fe:m: electro-motive force induced by winding coil on the harvester GðkÞ imaginary part of Theodorsen function I time-varying induced current in the winding coil (electrical system) [amperes] I0α total mass moment of inertia of the rotating bladeairfoil [kg∙m2] Imax Maximum current during limit-cycle post-critical stages (one period) [amperes] bi imaginary unit k k* kα LC ℓ ℓC M0z b 0z;ae;u M Me:m: m N Pin , Pout p

γ* ε ζp ζα ι ιmax κ λRL μae;j ; νae;j ρ τ ΦðsÞ Φe:m:c:

reduced frequency reduced frequency of the harvester at flutter still-air reduced angular frequency of the bladeairfoil mechanism inductance of the coil (electrical system) [henries] longitudinal length (span or height) of the blade-airfoil [m] effective coil length (electrical system) total dimensional torsional moment about pole O [N∙m]

Ψ ω* ωp

unsteady dimensionless aeroelastic torque about pole O total dimensional electromotive torque about pole O [N∙m] mass of the 1dof model illustrating electro-mechanical coupling [kg] number of winding coils (electrical system) Maximum (instantaneous) input and output powers [watts] sliding dof of the magnet inside the winding coil

ωα Operators ð⋅Þ' ; ð⋅Þ''

(electrical system) nonlinear vector-function for post-critical simulations aerodynamic functions needed to solve for flutter resistance of the winding coil [ohms] dimensionless time, s ¼ Ut=b continuous time [s] wind speed [m/s] minimum operational wind speed [m/s] State vector of the coupled electro-mechanical system angular rotation of the blade-airfoil about pole O initial rotation at time τ ¼ 0 (system's release) maximum rotation amplitude of the blade-airfoil during limit-cycle post-critical stages frequency ratio at flutter γ * ¼ ωα =ω*

generalized inertia parameter, ε ¼ πρb4 ℓðI0α Þ1 damping ratio of the 1dof model (illustrating coupling) structural damping ratio of the harvester time-varying dimensionless induced current of the harvester (electrical system) maximum value (magnitude) of the dimensionless induced current during limit-cycle post-critical stages nonlinear stiffness coefficient of the torsional spring model generalized impedance coefficient aeroelastic states with j ¼ f1; 2g air density [kg/m3] dimensionless time of the state-space model, τ ¼ kα s aeroelastic load function (indicial function) dimensional electro-mechanical coupling coefficient [N/ampere] Dimensionless electro-mechanical coupling coefficient angular frequency of the harvester at flutter [rad/s] angular frequency of the 1dof model (illustrating coupling) [rad/s] Still-air rotational frequency of the blade-airfoil [rad/s] first and second derivative with respect to s

post-critical flutter state, stable beyond the critical flutter speed. This phenomenon is a prerogative of long-span bridges with bluff deck girders (Scanlan and Tomko, 1971). It is also possible in the case of streamlined airfoils if the center (axis) of rotation of the airfoil's cross-section coincides with the leading edge of the profile (Bisplinghoff et al., 1955; Kakkavas, 1998). It is observed that the pitching rotation mechanism of an H-shaped cross section, similar to a bluff bridge deck and prone to torsional flutter, was proposed in the recent past to convert wind power to electric energy (Ahmadi, 1979). An H-shaped cross section is believed to be adequate for its propensity to exhibit torsional flutter. However, it also has several disadvantages. For example, sensitivity to incoming turbulence may lead to non-negligible buffeting vibration prior to the flutter onset, which needs to be avoided to prolong the lifetime of the apparatus. Furthermore, a bluff body would tend to produce larger drag loads, which must be considered while designing the structural support system and torsional-restoring force rotating mechanism. Additionally, larger wakes generated by a bluff body would reduce efficiency if an array of devices were used (for example, staggered from each other). The blade-airfoil cross-section and apparatus proposed herein are simple, versatile, and compact; they are possibly attractive in comparison

vibration frequencies are greater than 10 Hz (Priya and Inman, 2009) or the energy transfer is enhanced by repeated impacts (Zhu and Zhang, 2015). Nevertheless, the high cost of piezo-electric materials restricts applications to the scale of small devices. For applications at a larger scale, the use of vortex-induced vibration of “suspended” circular cylinders, connected by springs to a fixed-base support has been proposed (Bernitsas et al., 2008, 2009; Mehmood et al., 2013). For example, the VIVACE apparatus, designed for underwater applications (Bernitsas et al., 2008, 2009), exploits the lock-in harmonic vibration of the cylinder in the direction orthogonal to the flow. This device is ideal if the anticipated flow speed does not vary considerably. Nevertheless, the overall efficiency can significantly diminish in the case of irregular and turbulent flow conditions, typical of wind flow in the atmospheric boundary layer. Consequently, applications are possibly restricted to a narrow (lock-in) wind speed range. This paper builds on the results of a previous study that examined the technical feasibility of a wind-based energy harvester, taking advantage of the leading-edge torsional flutter instability of a blade-airfoil (Caracoglia, 2010) for wind energy conversion. Torsional flutter is a single-mode aeroelastic instability phenomenon, which triggers a diverging vibration of a flexible body. Sustained vibration is possible in a

2

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

Fig. 1. Conceptual model and schematics of the apparatus: (a) top view of a typical cross-section on horizontal plane, (b) close-up view illustrating the system for energy conversion and (c) 3D rendering.

with other harvesters, either using flapping airfoils or based on torsional flutter for power generation as described in a recent literature review study (Young et al., 2014). Since the energy conversion mechanism of the proposed apparatus has not been fully explained in the previous study by Caracoglia (2010), further investigation is carried out in this paper. The aim of this paper is to further expand the previous model, by accounting for the effects of electro-mechanical coupling in an effort to predict the electrical output current and evaluate the output during active energy conversion. In contrast with other wind-based harvesters the technology developed herein is less sensitive to wind turbulence. The operational mechanism of the apparatus is suitable for smaller-scale applications in the range of moderate wind speeds, where HAWT systems (Manwell et al., 2010) are either less efficient or impractical. A NACA 0012 airfoil section is employed to trigger the torsional flutter. Magnetic induction is employed for energy conversion (Kwon et al., 2013).

The novel extended model, presented in this work, is composed of seven nonlinear fully-coupled electro-mechanical dimensionless equations; it is written in state-space form. The aeroelastic forces, needed to trigger the harvesting mechanism, are simulated through unsteady lift theory for incompressible flow and indicial functions. Numerical integration is employed to solve the differential equations. Investigations reveal that the “aptitude” to energy conversion is controlled by Ψ, a dimensionless electro-mechanical coupling coefficient, and by λRL , a coefficient related to the generalized impedance of the power circuit. 2. Description of the apparatus The previous study identified the preliminary configuration of the proposed apparatus (Caracoglia, 2010). The main results are briefly summarized in this section. A schematic view of the apparatus is reported in Fig. 1. Fig. 1a presents the generic cross section of the harvester on a

3

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

reference horizontal plane; the z axis is oriented out of plane in this figure. Fig. 1b schematically illustrates the main principles of the electro-mechanical system for energy storage. Fig. 1c shows a 3D rendering of a hypothetical configuration of the apparatus, oriented with vertical rotation axis (z). This orientation is considered to avoid dependence on gravity loads. In the figure panels, the main body of the blade is rigid, has width (chord length) 2b and longitudinal length ℓ. The apparatus is composed of a blade-airfoil, a NACA0012 rigid blade. It is vertically connected to a support structure by means of a flexible torsional mechanism that enables the rotation α about the vertical axis (z) in the proximity of the leading edge (pole “O”, Fig. 1a). The reference zero-angle axis in the figure corresponds to the mean wind direction; U is the mean wind speed. The figure also illustrates the location of the shaft (transmission) mechanism connecting the apparatus to the electrical converter (Fig. 1b). The energy conversion principle has been inspired by an electromagnetic generator, exploiting traffic-induced vibration on bridge structures (Kwon et al., 2013). The shaft “B” moves (motion coordinate denoted by p in Fig. 1b) and is equipped with a permanent magnet (Kwon et al., 2013), which generates a magnetic field moving through a coil. The motion of the magnet through the coil induces an electromotive force (and voltage) in the coil. If the coil is connected to an external electrical circuit, the current I is induced. The induced current in turn produces a magnetic field with opposite polarity, which acts by retarding the transverse translation of the magnet (Donoso et al., 2009). This effect leads to electro-mechanical coupling within the system (Stephen, 2006). The produced energy must be stored in a battery for later use.

are derivatives with respect to dimensionless time s ¼ tU=b; CðkÞ ¼ FðkÞ þ biGðkÞ is the complex Theodorsen function (Theodorsen, 1935), bi is the imaginary unit and k is the reduced frequency. Three-dimensional load effect, typical of wings with finite ratios ℓ=b (Bisplinghoff et al., 1955), has not been included. Combination of Eqs. (1) and (2) yields (Caracoglia, 2010):

  α'' ð1 þ 9ε=8Þ þ α kα2 þ εk 2 R1 ðkÞ þ 1:5α' εkR2 ðkÞ ¼ 0:

πρb4 ℓðI0α Þ1 is a normalized inertia parameter, which identifies various operational regimes. Dissipation proportional to α' as a function of U and R2 is evident in Eq. (3). An increment of U causes a decrement of k and the destabilizing effect on α' . The unstable regime depends on the sign of R2 , which depends on k and may (or may not) become negative. In the case of leading-edge torsional flutter of airfoils some combinations of frequency kα and inertia parameter ε are necessary to produce aeroelastic instability and, consequently, the triggering of the apparatus. More detailed discussion is provided in a subsequent section. 3.2. Aeroelastic model and post-critical flutter vibration: discussion of the main findings Beyond flutter onset, a benign post-critical flutter regime is established by controlling the angles of attack during vibration. This is achieved, in this application, by introducing a nonlinear restoring force mechanism through suitable modeling of the torsional connection element at O in Fig. 1a. This mechanism is modeled as a nonlinear torsional spring element, which operates in parallel with the linear spring of the linear dynamic model. The post-critical 1dof dynamic equilibrium equation is found from Eqs. (1)–(3) by introducing two new terms: a Duffing oscillator spring model [e.g. Nayfeh and Mook (1995)] proportional to α3 with κ > 0 being a dimensionless parameter, and a linear term 2ζα kα simulating internal friction with structural damping ratio ζα . This transformation yields:

3. Reduced-order electro-mechanical coupled model 3.1. Aeroelastic model and critical flutter condition: discussion of the main findings The apparatus relies on vibration triggered by one degree-of-freedom (1dof) flutter, which must be appropriately controlled and calibrated to achieve adequate efficiency. The dynamics of the apparatus is modeled as a 1dof pitching (“flapping”) blade-airfoil of chord length 2b and height (longitudinal span length) ℓ, with longitudinal axis oriented in the z direction in Fig. 1. Rotation angles α and moments about pole O are taken positive if clockwise. The element, connecting the blade with the structural support, is composed of a flexible short plate (Fig. 1a); it is simulated in the model as a concentrated linear torsional spring placed at O. The dynamic equation (rotation α) without damping is:

d2 ½αðtÞ M0z þ ω2α ¼ : dt 2 I0α

 α'' ðsÞð1 þ 9ε=8Þ þ 1:5εα' ðsÞ þ 2ζ α kα α' ðsÞ þ kα2 αðsÞ þ κα3 ðsÞ ¼  εCðkÞ αðsÞ þ 1:5α' ðsÞ :

(4)

A more accurate model, capable of describing a wider-range of postcritical large-amplitude vibration regimes, should possibly be used instead of Eq. (4), for example accounting for nonlinear fluid-structure interaction in the damping term α' of Eq. (4), similar to the concept of van der Pol oscillator and based, for example, on the ONERA empirical model (Dunn and Dugundji, 1992; Tran and Petot, 1981). This more sophisticated model has not been considered for the sake of simplicity and because the present analysis primarily focuses on the coupled system. Eq. (4) describes, in a mixed time-frequency formulation, the dynamics induced by unsteady aeroelastic torque about O. In order to solve b 0z;ae;u ðsÞ ¼ εCðkÞ½αðsÞ þ 1:5α' ðsÞ Eq. (4) in s time domain, the torque M on the right-hand side must be expressed through indicial functions (Bisplinghoff et al., 1955):

(1)

The torsional aeroelastic moment M0z in Eq. (1) is referenced to O in the figure. The torque depends on the uniform wind field of mean speed U with direction parallel to the x axis in the figure. I0α is the total polar mass moment of inertia of the moving apparatus about O; ωα is the angular frequency and t is the time variable. The effect of structural damping will be later included. Mean aerodynamic forces are approximately zero; static lift force is negligible due to blade symmetry if α  0. Other aerodynamic loads on the structural support and components are not considered. The unsteady moment M0z is derived from classical aerodynamic theory (Bisplinghoff et al., 1955; Dowell et al., 2004) by combining the effect of transverse lift force L(ae)h and torsional moment M(ae)α, schematically illustrated in Fig. 1a. The moment becomes

      9 M0z ¼ πρb2 U 2 ℓ  1:5α' þ α''  α þ 1:5α' CðkÞ 8

(3)

The two functions R1 ðkÞ ¼ k ½FðkÞ þ biGðkÞ and R2 ðkÞ ¼ k1 ½1 þ FðkÞ þ biGðkÞ in Eq. (3) are complex and depend on k. The coefficient ε ¼ 2

2 b 0z;ae;u ðsÞ ¼ M

(2)

3

6 1:5∫ s Φ' ðσÞα' ðs  σÞdσ þ 1:5Φ0 α' ðsÞ 7 0 7 ε6 4 þ∫ s Φ' ðσÞαðs  σÞdσ þ Φ0 αðsÞ 5: 0

(5)

Eq. (5) is written in a non-standard form, compared to Bisplinghoff et al. (1955). This equivalent form is valid for s  0 only; it is still acceptable in the present context and investigation. The function ΦðsÞ ¼ ½1  c1 expðd1 sÞ  c2 expðd2 sÞ is the R.T. Jones' model of the Wagner

In Eq. (2) k ¼ ωb=U is a generic reduced frequency, kα ¼ ωα b=U is the reduced rotational frequency of the oscillator, α' ¼ dα=ds and ¼ d2 α=ds2

4

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

Φe:m:c: ¼ NBℓC (Kazmierski and Beeby, 2011) with N number of coil turns, B average magnetic flux density inside the collar of the winding coil and ℓC effective coil length (ℓC ≠ℓ). Equations (8) and (9) can be adapted to study the present harvester, after observing that the dimensional electro-motive torque Me:m: with respect to pole O in Fig. 1b (Shaft A) can be approximately determined as Me:m:  2bfe:m: for moderate vibration amplitudes and small α. This assumption implies that a concentrated electro-motive force [Eq. (8)] is applied at the leeward edge of the blade-airfoil. Therefore, the transverse translational velocity (dp=dt) in Eq. (8) can be interpreted as the vibration of the shaft B in Fig. 1b with p  2bα for small pitch angles. The comparison with the translating configuration described by Eq. (8) suggests that an electro-motive torque Me:m: ¼ 2Φe:m:c: IðtÞ [Nm] can formally be added to the unsteady aeroelastic torque on the right-hand side of the dimensional dynamic equation [equivalent to Eq. (7)]. If Eq. (9) is adapted to the rotating dynamics of the proposed apparatus with p  2bα and dimensionless time τ is utilized, the following equation is obtained:

function (Bisplinghoff et al., 1955; Jones, 1939) with Φ0 ¼ Φðs ¼ 0Þ ¼ ½1  c1  c2 . The theoretical form of the R.T. Jones’ formula is initially employed in the simulations discussed in Sections 5.1–5.3; the parameters are c1 ¼ 0:165, d1 ¼ 0:0455, c2 ¼ 0:335, d2 ¼ 0:300. Since the aerodynamic loads of the blade-airfoil may differ from the theoretical values in a realistic implementation (Fig. 1), a supplementary investigation is presented in Section 5.4 to examine the influence of a variation in the load model parameters on the system dynamics. Irrespective of the numerical values of the parameters, it is convenient to rewrite the integrals in Eq. (5) as aeroelastic states νae;1 , νae;2 in Eq. (6a), and aeroelastic states μae;1 , μae;2 in Eq. (6b) below. It must be noted that the following equations are valid provided that a lag-two aeroelastic model is employed to describe the unsteady loads. s

∫ 0 Φ' ðσÞα' ðs  σÞdσ ¼ ½νae;1 ðsÞ   þ ½νae;2 ðsÞ  ¼   s d1 s s c1 d1 ∫ 0 e α'ðs  σÞdσ þ c2 d2 ∫ 0 ed2 s α'ðs  σÞdσ ;   s ∫ 0 Φ' ðσÞαðs  σÞdσ ¼ μae;1 ðsÞ  þ μae;2 ðsÞ ¼   s s c1 d1 ∫ 0 ed1 s αðs  σÞdσ þ c2 d2 ∫ 0 ed2 s αðs  σÞdσ :

(6a)

d½IðτÞ 2bΦe:m:c: d½αðτÞ RC ¼  IðτÞ: dτ dt LC ωα LC

(6b)

It is also convenient to define the dimensionless induced current ιðtÞ from the relationship IðτÞ ¼ ιðtÞ½Φe:m:c: 2bωα =RC . The quantity ½Φe:m:c: 2bωα =RC  has in fact the dimensions of a reference current [A]. If the angular velocity is constant and equal to one (rad/s), the left-hand side of Eq. (10) must vanish, and the induction effect produces a constant-amplitude current equal to this reference current. Consequently, Eq. (10) can be simplified and written in dimensionless form as:

For further implementation, dimensionless time τ ¼ ωα t ¼ kα s is used in Eqs. (4)–(6); subsequent transformations lead to a new dynamic equilibrium equation:

    9 d2 α 1:5ε dα  þ α þ κα3 ¼ 1þ ε þ 2ζα 2 þ 8 kα dτ dτ   2 3 dα þ Φ α þ 1:5k 0 α 7 ε6 dτ 7:  26 5 kα 4 þ1:5ðνae;1 þ νae;2 Þ þ μae;1 þ μae;2 

    dι RC dα dα ¼  ι ¼ λRL ι : dτ ωα LC dt dt

(7)

    9 d2 α 1:5ε dα  þ α þ κα3 ¼ 1þ ε þ þ 2ζ α 2 8 k dτ dτ α     ε 3 dα 3 þ ðνae;1 þ νae;2 Þ þ μae;1 þ μae;2  Ψι  2 Φ0 α þ kα 2 dτ 2 kα 

The model is based on the findings documented by Kwon et al. (2013), which have been adapted to describe the behavior of the proposed apparatus. Kwon et al. (2013) present a sensor device, designed to harvest energy from the vertical vibration of a 1dof mass-spring-dashpot system through eddy currents. Inspired by Kwon et al. (2013) and the derivations in Priya and Inman (2009), this study first examines the oscillation of a 1dof mechanical system in the p vertical direction, after noticing the similarity with the vibration inside the magnetic coil shown in Fig. 1b. This system configuration leads to electro-mechanical coupling. The equations of this coupled system are, in their simplest form,

d2 ½pðtÞ d½pðtÞ þ mω2p pðtÞ ¼ fext þ fe:m: ; þ 2mζ p ωp dt 2 dt

d½IðtÞ Φe:m:c: d½pðtÞ RC ¼  IðtÞ: dt dt LC LC

(11)

In the previous equation λRL ¼ RC =ðωα LC Þ is a generalized impedance coefficient of the power circuit. If the effect of Me.m. is included on the right-hand side of Eq. (7) with IðτÞ ¼ ιðtÞ½Φe:m:c: 2bωα =RC , the following dynamic equation is found:

3.3. Modeling the electrical system and output current equation

m

(10)

(12)

in which Ψ ¼ 4b2 ðΦe:m:c: Þ2 =ðωα I0α RC Þ is designated as the dimensionless electro-mechanical coupling coefficient of the apparatus. 3.4. Resultant model of the coupled electro-mechanical system Differentiating the right-hand side of Eq. (6) with respect to τ ¼ ωα t ¼ kα s leads to four first-order differential equations of the aeroelastic states (with j ¼ f1; 2g)

(8)

dνae;j dα dj ¼ dj cj  νae;j ; dτ kα dτ

(13)

 dμae;j dj  ¼ cj α  μae;j : dτ kα

(14)

(9)

In the dynamic equation [Eq. (8)] the mass, structural angular frequency and damping ratio are denoted by m, ωp and ζp . Eq. (9) describes the power output relationship of the electrical circuit, secondary system or “coil C”, with IðtÞ induced current (in amperes, A), LC inductance of the winding coil (in henries) and RC resistance (in ohms). The external loads are an external force fext and the electro-motive force fe:m: . The fe:m: is found in accordance with Faraday's law of magnetic induction (Kwon et al., 2013; Stephen, 2006) as fe:m: ¼ Φe:m:c: IðtÞ. The quantity Φe:m:c: , which couples the two equations above, is the electro-mechanical coupling coefficient in the units of newton/ampere [N/A]. If the motion of the magnet in Fig. 1b is completely immersed in the coil, we find

Eqs. (13) and (14) are combined with Eqs. (11) and (12) to find the nonlinear coupled electro-mechanical (em) state-space model. The state vector and dynamic model equations are defined as follows, with T denoting transpose operator:

 T dα wem ðτÞ ¼ αðτÞ; ; νae;1 ðτÞ; νae;2 ðτÞ; μae;1 ðτÞ; μae;2 ðτÞ; ιðτÞ dτ ¼ ½wem;1 ; …; wem;7 T

5

(15)

L. Caracoglia

qem

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

wem;2 8 > >  8 9 > 1:5ε  > > > > 3 > > >  w þ κw þ 2ζ  w > > > em;1 em;2 α em;1 > > > > kα > > > > > > > > > > > > > > > > > > > > > < = 2 3 >   > 8 > 3 3 > > Φ k ðw w þ þ w þ w Þ > 0 em;1 α em;2 em;3 em;4 > > 8 þ 9 6 7 > > > > 2 2 > > 7 > ε6 > > > > 6 7 > > >  Ψwem;7 >  26 > > > 7 < > > k > > 5 α4 > > > > ¼ : ; > þwem;5 þ wem;6 > > > > > > > > > c1 d1 wem;2  d1 kα1 wem;3 > > > > > > c2 d2 wem;2  d2 kα1 wem;4 > > > > > > kα1 d1 ðc1 wem;1  wem;5 Þ > > > > > > k1 d2 ðc2 wem;1  wem;6 Þ : α

9 > > > > > > > > > > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > > > > > > > > > > ;

;

(16)

λRL ðwem;2  wem;7 Þ

dwem ¼ qem ðwem Þ: dτ

apparatus (ωα ) and angular frequency of the harmonic aeroelastic vibration. Eq. (19) is solved first to obtain the reduced frequency at incipient flutter. Inspection reveals that there is only one root (k ¼ k*  0:04034), which is independent of the harvester geometry and configuration. This remark indicates that, when structural damping is negligible (ζα ¼ 0), alteration of the flutter condition by modification of specific physical properties, at least theoretically, is not possible for improving the performance of the harvester. This is a limitation of the proposed apparatus in comparison with other flutter mills, which exploit the coupled plunging-pitching instability and can be more suitably fine-tuned. In addition, even though unsteady aerodynamics is necessary for triggering flutter, the solution seems approximately quasi-steady in terms of the loads. The first equation [Eq. (18)] is sequentially solved after Eq. (19) and substitution of k ¼ k* to get the frequency ratio at incipient flutter. This substitution leads to the derivation of γ * ¼ ωα =ω* where ω* becomes the dimensional flutter vibration angular frequency. Examination of Eq. (18) suggests that the squared frequency ratio (Caracoglia, 2010) at flutter

(17)

The vector function qem in Eq. (16) is seven-dimensional; the top two terms describe the dynamics of the mechanical oscillator [Eq. (12)], the subsequent four scalar equations pertain to the unsteady fluid-oscillator model [Eqs. (13) and (14)] and the last one is the induced current equation [Eq. (11)]. The dynamic model in state-space form is concisely written as in Eq. (17) using qem ; it is nonlinear because of the effect of the Duffing-type spring model acting on wem;1 with parameter κ. 4. Practical considerations for the preliminary design of the harvester A previous study (Caracoglia, 2010) has shown that Eq. (4) can be solved by postulating an incipient flutter condition exhibiting simple harmonic motion. The linear flutter problem (κ ¼ 0) is solved after transformation of this equation into two algebraic equations that need to be simultaneously solved to determine the two unknown variables: the reduced frequency at incipient instability k* and the angular frequency ω* at developing flutter. This step is traditionally encountered in aeroelasticity (Scanlan and Rosenbaum, 1968). Even though derivation and solution of the two equations has been described in a previous study (Caracoglia, 2010), new developments and further examination are discussed in this section in an attempt to determine minimum operational regimes and potential limitations of the harvester, useful for future implementation. The two equations, i.e., the real and imaginary parts of the complex-valued flutter algebraic equation, are reported below in the case of zero structural damping (ζα ¼ 0):

  FðkÞ 3 GðkÞ 9 þ ; γ2 ¼ 1 þ ε  2 þ k 2 k 8

(18)

2GðkÞ þ 3ð1 þ FðkÞÞk ¼ 0:

(19)

ðγ * Þ2 ¼ ðωα =ω* Þ2 is a linear function of the generalized inertia parameter ε only, after noticing that the quantity enclosed within brackets in the equation is constant and independent of other quantities. In contrast with the above-described limitation, this aspect is perceived as an advantage since design of the apparatus is essentially controlled by a single parameter ε and is consequently straightforward. Investigation on the relationship between frequency ratio γ * ¼ ωα =ω* and generalized inertia parameter ε at incipient flutter is presented in Fig. 2a. In the figure panel the relationship between ðγ * Þ2 ¼ ðωα =ω* Þ2 and ε > 0, under the condition k ¼ k* , can be seen (thin continuous line). The relationship between ðγ * Þ2 and ε is linear and monotonically decreases toward zero. No physical solution is possible for ε > εmax  1:748E  3; this value becomes an upper limit beyond which torsional flutter is impossible. The design of the harvester is thus restricted to a narrow combination of variables b, I0α and ℓ. This behavior is not unexpected; it has been explained in the reference literature (Bisplinghoff et al., 1955) that some combinations of the system parameters lead to a non-flutter condition, which must be avoided from an energy harvesting viewpoint. Furthermore, the relationship between the frequency ratio γ * and ε in the range 0  ε  εmax (thick continuous line) becomes approximately pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi γ *  1  572:2ε (Fig. 2a), suggesting a weak dependence of the flutter frequency on the generalized inertia parameter ε only and, consequently,

The form of the equations above is slightly different but equivalent to the original formulae derived in Caracoglia (2010); this particular form is meant to highlight some specific properties of the solution, which are relevant to the efficiency of the harvester. In Eq. (18) the quantity γ * ¼ ωα =ω* is a “frequency ratio” between intrinsic angular frequency of the

6

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

concerns with a finite-length airfoil would be the reduction of lift coefficient slope (2π) with the decrement of the aspect ratio ℓ=b and, consequently, the slope of the aerodynamic torque. Examination of more specific formulas for airfoils (Argentina and Mahadevan, 2005; Bisplinghoff et al., 1955) suggests that the two-dimensional solution of the aerodynamic lift is essentially attained (80%) with ℓ=b > 8.

a limited variation of ω* from the still-air angular frequency of the mechanical apparatus (ωα ). This remark leads to a simple and approximate expression for the critical flutter speed U * at which the device becomes operational:

U ¼

    ω b 1 24:788 p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi : ¼ ω b b  ω α α k k γ  1  572:2ε

(20)

5. Post-critical operational stages and energy conversion

This expression is only approximate in the presence of structural damping and electro-mechanical coupling but provides some useful indications. It suggests that the minimum operational speed is linearly proportional to both b and ωα and, therefore, small compact devices working at low frequencies are preferable in order to maintain U * as low as possible. The quantity within brackets in Eq. (20) is independent of the system parameters apart from ε. It is therefore labeled as “universal curve” of the normalized flutter speed 1=ðk* γ * Þ for this harvester. The quantity 1=ðk* γ * Þ is plotted in Fig. 2b as a function of ε > 0. The curve is a “square-root hyperbola”, which has a moderately flat graph, insensitive to the main system parameters or configuration for a relatively wide interval of ε. The function approaches the vertical asymptote at ε ¼ εmax  1:748E  3 where the flutter speed becomes infinity. This remark confirms that this apparatus may have a good application range. Moreover, the flutter speed curve in Fig. 2a is a unique and general result. This function may be used to represent any physical configuration through suitable rescaling [Eq. (20)] and does not require any additional item to identify the various operational regimes. It provides a physical clarification for the partially unexplained similarities observed in Caracoglia (2010) among several examined cases. Finally, the inspection of Eq. (20) and Fig. 2 also indicates that, if the frequency of the apparatus ½ωα =ð2πÞ and the inertia per unit length ðI0α =ℓÞ can be preselected, the longitudinal length ℓ does not influence the flutter solution. Thus, simulations and calculations may be referred to a unit-length apparatus (ℓ ¼ 1), as implicitly assumed in the remainder of the paper. Adequacy of this assumption implies negligible threedimensional flow effects, typical of finite-length airfoils. One of the

5.1. Examination of post-critical dynamics by numerical integration After noticing the proportionality of operational (or flutter) wind speed U * with angular frequency ωα and reference chord b through Eq. (20), a representative set of harvester configurations is selected. The selection is based on the constraint that the apparatus should be ideally operational in the range of wind speeds between 10 m/s and 15 m/s, implying that the angular frequency should approximately correspond to ωα variable between 3 and 5 s for a “compact” apparatus of moderate dimensions (0 < b < 1 m). Thus, three harvester configurations are considered and are investigated in detail: Type 0, Type 1 and Type 2. The fundamental properties of each configuration along with the minimum operational U * speed in the absence of structural damping (ζα ¼ 0), determined from the linear theory [Section 4 and Eq. (18)], can be found in Table 1. The free vibrations at various U and with a given set of initial conditions [Eq. (17)] are solved by numerical integration. A fifth-order Runge-Kutta-Fehlberg method (Sewell, 2005) is used for this purpose with time step suitably selected. A non-zero initial rotation angle αI is necessary to trigger the apparatus. The initial trigger condition is set to αI ¼ 0:01 radians in all the simulations. Eq. (17) is consequently solved with state vector at initial time wem ðτ ¼ 0Þ ¼ ½αI ; 0; …; 0T . As explained in Section 4, in the absence of structural damping (ζα ¼ 0) and nonlinearity (κ ¼ 0) the minimum flutter threshold is controlled by the inertia parameter ε (a function of both I0α =ℓ and b) and is proportional to b and ωα (Fig. 2). The configurations presented in Table 1 are therefore selected from this remark; they represent a lowmass, low-frequency apparatus. An example of dynamic rotation, numerically solved, is presented in Fig. 3a for a Type-0 apparatus and in Fig. 3b for a Type-1 one; two of the basic configurations from Table 1 are studied. In both cases a nonlinear torsional-spring restoring mechanism with κ ¼ 100 is chosen from previous studies (Caracoglia, 2010). Type-0 device is designed with ε  7:7E  4 whereas Type-1 device has ε  1:6E  4. Vertical scales are not the same in the panels of Fig. 3 to examine stable vs. unstable dynamic conditions. One unit of time τ=ð2πÞ on the horizontal axis approximately coincides with one vibration period. Fig. 3a and b illustrate both a stable case in the top two panels (wind speed U ¼ 10 m/s) and one unstable case in the bottom two panels (wind speed U ¼ 15 m/s) for the nonlinear harvester. Time histories of the dynamic rotation (α) and dimensionless induced current (ι) are respectively presented on the left-side and right-side panels. In the stable preflutter scenario (top panels) the system dissipates energy and vibration decays from the initial angle αI ¼ 0:01. The electro-mechanical coupling effect, exhibited as an indirect induced current ι (top right panels in both Fig. 3a and b), acts in this case as an additional source of damping. In the post-flutter scenario, the model predicts that the pitching of the bladeairfoil is fully developed for Type-0 after about approximately 300 cycles (bottom panels of Fig. 3a), whereas it is not fully sustained for the second device (Fig. 3b). In the second case, a longer exposure to wind would possibly be needed to achieve self-sustained harmonic motion. More detailed investigations have been carried out to examine the post-flutter vibrations in terms of energy efficiency and dimensionless induced current (ι) for other combinations of the design parameters, starting from the three basic configurations summarized in Table 1. The results of the parametric study are presented in the next sub-sections to

Fig. 2. Investigation on the relationship between frequency ratio γ * ¼ ωα =ω* and generalized inertia parameter ε at incipient flutter (with critical flutter reduced frequency k ¼ k*  0:0403): (a) γ * and ðγ * Þ2 vs. ε, (b) “universal curve” of the normalized leadingedge torsional flutter speed 1=ðk* γ * Þ vs. ε. 7

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

Table 1 Properties of the nonlinear harvester (basic reference configurations). Basic config.

ωα =ð2πÞ [Hz]

b [m]

I0α =ℓ [kg-m2/m]

ζ α [%]

κ []

Ψ []

λRL []

U * [m/s]

Type 0 Type 1 Type 2

0.25 0.20 0.10

0.25 0.20 0.50

20 40 300

0.25(y) 0.30(y) 0.30(y)

100 100 100

0.01(#) 0.01 (#) 0.01 (#)

0.75(y) 0.75(y) 0.75(y)

13.0(z) 12.7(z) 10.7(z)

Notes: (#) Ψ will be later varied in Fig. 4, Fig. 8 and Fig. 9; (y) ζ α and λRL will be later varied in Fig. 5 and Fig. 6; (z) estimated by linear flutter analysis with ζ α ¼ 0 and Eq. (20).

Fig. 3. Electro-Dynamic post-critical analysis for Type-0 (a) and Type-1 (b) basic harvester configuration (Table 1) at wind speeds U ¼ f10; 15g [m/s]. Dynamic rotation α about pivot axis (pole O) and dimensionless current ι vs. dimensionless time τ=ð2πÞ.

The quantity ιmax is defined as the maximum (i.e., the magnitude or absolute value) of the output current ι observed during limit-cycle postcritical vibration. This reference quantity is determined by numerical simulation over an interval between 0 and 400 cycles [τ=ð2πÞ]. The normalization of this quantity, according to ιðtÞ ¼ IðτÞ=½Φe:m:c: 2bωα =RC , also ensures independence on the geometry and dimensions of the apparatus and, therefore, enables the cross-comparison among the various configurations in Table 1. The duration of the simulations ensures that an asymptotic value ιmax is approximately attained because of the sustained limit-cycle oscillation. Further inspection of the results (not shown) has confirmed that most data points, presented in Fig. 4, belong

evaluate the efficiency of the harvester through the study of a representative value of induced current. 5.2. Influence of the dimensionless electro-mechanical coupling coefficient (Ψ ) Fig. 4 illustrates, on the vertical axis of each panel, the maximum dimensionless output current of the apparatus (ιmax ) at U ¼ 15 m/s. Results are obtained for all three basic configurations by varying the coupling coefficient Ψ from its reference values. The reference “starshaped” points correspond to the configurations described in Table 1. 8

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

between Ψ and ιmax around Ψ ¼ 1E  4 but are not shown in the figure. They seem to be influenced in part either by the numerical integration or the computation of ιmax . Practically speaking, numerical simulations are interrupted below this small Ψ value. 5.3. Influence of the generalized impedance coefficient (λRL) and structural damping ratio (ζ α ) Fig. 5 examines the combined effect induced by a variation of generalized impedance coefficient (λRL ) and damping ratio (ζα ) for all three basic devices equipped with reference eddy-current system. The basic cases with electro-mechanical coupling coefficient Ψ ¼ 0:01, described in Table 1, are considered. The projected values of ιmax at U ¼ 15 m/s are compared, in this parametric study, to the reference values for Type-0 device (Fig. 5a), Type-1 device (Fig. 5b) and Type-2 device (Fig. 5c). The vertical-scale axis is not the same in the figure panels to point out the differences. In the case of Type-0 (Fig. 5a) and Type-1 (Fig. 5b) devices a smaller scale is used to discern small output current

(a)

(b)

(c) Fig. 4. Maximum output current (ιmax ) of the harvester vs. electro-mechanical coupling coefficient (Ψ) at U ¼ 15 m/s: (a) Type 0, (b) Type 1, (c) Type 2. Results for the three basic configurations (Table 1) are indicated by an asterisk symbol; the horizontal-axis scale is not the same in the panels.

to this category, and that very few data points exhibit partly-attained post-critical pitching. Analysis of Fig. 4 reveals that efficiency is inversely proportional to Ψ and that the relationship is nonlinear. No power can be extracted if Ψ exceeds 0.03 for Type-0 and Type-1 devices, whereas the Type-2 configuration (Fig. 4c) appears to be more adequate over a wider range of Ψ. The horizontal-axis scale has been expanded in Fig. 4c in comparison with the upper panels to examine the whole range of Ψ with nonzero ιmax . Nevertheless, even though the range of Ψ exhibiting active energy conversion reduces for Type-1 device (Fig. 4b), the largest value of ιmax is associated with, for this type of apparatus, a relative increment of output current greater than 40%. It must be noted that the output current curve must decrease to zero as Ψ tends to zero. This first parametric investigation suggests that a low-frequency large-dimension apparatus (such as Type 2), which can convert energy at low wind speeds [according to Eq. (20)] over a wider interval of Ψ, may not always be efficient compared to a Type-1 configuration. Small variations and some discontinuities are also noted in the relationship

Fig. 5. Maximum output current (ιmax ) vs. generalized impedance coefficient (λRL ) at U ¼ 15 m/s for a harvester with Ψ ¼ 0:01 and variable damping ratio (ζ α ): (a) Type 0, (b) Type 1, (c) Type 2. Results for the three basic configurations (Table 1) are indicated by an asterisk symbol; the vertical-scale axis is not the same in the panels. 9

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

a comparable interval of λRL . Contrary to Type-0 device, the sensitivity to structural damping is less pronounced in a larger apparatus with higher frequency and more difficult to start [Eq. (20)]. For example, the performance still seems acceptable with ζ α ¼ 0:7% (dashed line in Fig. 6).

values at larger damping ratios. The damping ratios are gradually increased from the original values reported in Table 1. The onset of torsional flutter is mainly controlled by damping, either induced by fluidstructure interaction or influenced by simulated eddy-current through electro-mechanical coupling. Therefore, a relative increment of any damping source tends to suppress the vibration and reduces the performance of the harvester. Examination of Fig. 5 confirms that low structural damping ratio is preferable because it ensures active operational conditions over a wider range of λRL . The ιmax curve found for Type-1 device with ζα ¼ 0:3% (dashed-dotted line in Fig. 5b) is qualitatively similar to the curve obtained for Type-0 device with ζα ¼ 0:25% (dashed-dotted line in Fig. 5a). Both the range 0 < λRL < 0:2 and the range λRL > 1:0 must be avoided in these two configurations. A small increment of structural damping (from 0.25% to 0.40%) can reduce by more than 50% the induced output current in the case of Type-0 device (Fig. 5a). Furthermore, no energy transfer is already observed at moderate damping above ζα ¼ 0:3% for Type-1 device. Further investigations appear to indicate that the aptitude to energy conversion is controlled by a combination of the parameters λRL and ζα . Small variations in the parameters noticeably influence the performance and reduce the output current ιmax . In general, ζα values larger than 0.5% are not acceptable since the onset of active vibration is shifted toward higher U. The properties of Type-0 configuration, such as b, ωα and I0α , appear to be more suitable since the potential operational range of the device is warranted for a broader set of Ψ and λRL . A preliminary recommendation for the parameters of Type-0 and Type-1 devices in the range of the examined wind speeds is: 0:4 < λRL < 1:0 and 0:01 < Ψ < 0:03. In contrast, Type-2 device (Fig. 5c) exhibits a remarkably good behavior over a wide range of damping ratios between ζ α ¼ 0:3% and ζ α ¼ 0:9%. The three figure panels also suggest that the reference configuration (star-shaped point on the graph) does not coincide with the best performance (contrary to the points used for Type-0 and Type-1 devices). The interval of λRL , corresponding to an active device, appears to double for λRL > 2:0 with significant growth in the values of ιmax . Therefore, a supplementary investigation on output current (ιmax ) vs. generalized impedance coefficient (λRL ) is carried out for Type-2 harvester at U ¼ 15 m/s with variable ζα by increasing the electromechanical coupling coefficient from Ψ ¼ 0:01 to Ψ ¼ 0:025. The results of this investigation are summarized in Fig. 6. In Fig. 6 the device remains operational for 0 < λRL < 1:6 with the largest ιmax values above 0.08, obtained with ζ α ¼ 0:3%. As expected, reduction of the interval of λRL is accompanied by a diminution of efficiency. For example, a 20% diminution of ιmax is noted at λRL ¼ 1:0 with ζ α ¼ 0:3% between Fig. 5c and Fig. 6 (dashed-dotted lines in both figures). Interestingly, the ιmax vs. λRL curves are qualitatively similar to the numerical solutions found with Type-0 device (Fig. 5a), i.e., defined over

Dynamic response of the harvester is undoubtedly influenced by the shape of the profile undergoing fluid-structure interaction. The choice of a basic blade-airfoil profile (NACA0012) is good because the apparatus is simple. Nevertheless, an assumption has been made (small rotations) to justify the choice of the unsteady aeroelastic torque loads by Eqs. (5) and (6). Inspection of the dynamic solution substantially reveals that vibration amplitudes (e.g., Fig. 3) are of the order of 5–10 , which are still an acceptable value to ensure approximate validity of the initial assumption. Nevertheless, a more appropriate model, incorporating nonlinear features of the lift curve as a function of the attack angle α may be necessary (Dunn and Dugundji, 1992; Tran and Petot, 1981), especially at larger α. A second important remark is related to the general features of the loads. In fact, local wind pressures may not exclusively be controlled by the geometry of the moving blade but potentially affected either by interference or shadow effect, caused by the support mast (Fig. 1). Furthermore, it has been recently observed that interference with the ground may influence the aeroelastic behavior of flexible flags (Dessi and Mazzocconi, 2015), which are similar to the proposed apparatus (at least in principle). The configuration of Fig. 1, mounted on a building roof, may consequently have an impact on the flow-harvester dynamic interaction. Therefore, it is pertinent to investigate potential variations in the aeroelastic loads and deviations from the theoretical formula of the airfoil's indicial functions postulated by Jones (1939). This problem is often encountered in the field of bridge aerodynamics and bluff bodies (Scanlan et al., 1974). It is usually resolved by wind tunnel experimentation to find the loads. Since the main objective of this study is the derivation of the coupled electro-mechanical interaction model, the problem of the aeroelastic load is examined by parametric numerical investigation only. Supplementary experimental studies are beyond the scope of this paper. The set of alternative load cases includes cases “D” “E” and “C”, which are presented in Fig. 7. Three aeroelastic load cases are constructed from the reference model by R.T. Jones, ΦðsÞ ¼ ½1  c1 expðd1 sÞ  c2 expðd2 sÞ with s ¼ τ=kα , by synthetic modification of the parameters c1 ¼ 0:165, d1 ¼ 0:0455, c2 ¼ 0:335, d2 ¼ 0:300. The altered coefficients are found by proportionally changing all the values of the parameters through multiplication by a constant. The variations in the coefficients of the function ΦðsÞ are moderate in load cases D (1.2 times increment) and E (1.1 times increment), whereas a reduction (0.5 times) is examined with load case C. The

Fig. 6. Additional investigation on output current (ιmax ) vs. generalized impedance coefficient (λRL ) for Type-2 harvester at U ¼ 15 m/s with Ψ ¼ 0:025 and variable damping ratio (ζ α ).

Fig. 7. Unsteady torque load indicial functions Φ vs. dimensionless time (s ¼ τ=kα ) associated with the aeroelastic load cases D, E, C as well as Jones (1939), used in the parametric study.

5.4. Examining the effects of a variation in the unsteady aeroelastic loads

10

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

that, despite the relative small variations in the unsteadiness of the aeroelastic loads, remarkable differences are exhibited by the output current curves. In the case of Type-0 (Fig. 8a) and Type-2 (Fig. 8c) devices, the range of Ψ values that correspond to an actively engaged device doubles. This observation indicates a nonlinear relationship between the load variations and the output of the electrical circuit (particularly for cases D and E). The relative increment of the output ιmax from the reference R.T. Jones’ configuration (dashed-dotted line in all the figure panels) approximately varies between þ20% and to þ40% with case D. The gain of ιmax is more moderate for Type-2 device (Fig. 8b). Even though the output behavior has distinct features for each apparatus, the shape of the curves is very similar for load cases D and E in all the panels. This remark suggests that the physical mechanism driving the aeroelastic interaction and the consequent energy transfer have not been altered by the load variations (at least at U ¼ 15 m/s). The flutter condition and the post-critical dynamics are approximately the same as before. In contrast, a distinct diminution of the load unsteadiness (aeroelastic case C) induces an abrupt change in the output. No energy can be extracted from the system irrespective of the type of device. Fluttering is consequently suppressed. This result is not unexpected since torsional flutter for such a streamlined profile requires the load model to be unsteady. Furthermore, it also confirms the relevance of accurate examination of the loads, which may potentially lead to an inefficient output. It must be noted that the condition Ψ→0 leads to ιmax →0 since the coupling among the system equations [Eq. (16)] is excluded. This aspect is perhaps not immediately perceivable from the graphs of ιmax vs. Ψ (Fig. 4 and Fig. 8) due to the choice of the scale and because a small but non-zero lower limit of the order of 1E-4 has been selected for Ψ > 0 in the numerical simulations. A more detailed investigation should possibly include analysis of smaller values of Ψ to better identify this special condition. Finally, since evaluation of ιmax is carried out by data point inspection using the results of the numerical solver, small numerical issues are also possible; these aspects warrant special consideration.

indicial functions of the various cases are plotted in Fig. 7 as a function of dimensionless time s and compared to the reference solution by R.T. Jones. An increase in the model coefficients coincides with a more pronounced unsteadiness in the loads; the function ΦðsÞ in Fig. 7 for cases D and E asymptotically tends to the steady aerodynamic torque [ΦðsÞ→1] at a slower rate. On the contrary, case C evaluates a situation in which the transition time to a steady aerodynamic load condition is very rapid. Even though selection of the ΦðsÞ model parameters may conceivably require more evidence and verification, it seems appropriate for this parametric analysis. Examination of Fig. 7 also suggests that the overall trend of the synthetic unsteady indicial functions is plausible with the behavior of a predominantly streamlined control surface [e.g., absence of “overshooting effect”, as explained in Scanlan et al. (1974)]. Fig. 8 illustrates the main results of this parametric investigation. The figure examines the maximum output current coefficient vs. the coupling coefficient Ψ for the various aeroelastic loads cases. All the three device types are considered at the wind speed U ¼ 15 m/s and with constant generalized impedance λRL ¼ 0:75. Inspection of cases D and E reveals

6. Preliminary evaluation of output power Even though the induced current ιmax may be employed for estimation of the output power [for example approximately as λRL ðιmax Þ2 in dimensionless units], the dimensionless variables ι or ιmax may not be fully adequate to study the performance. The selected normalization of the physical and dimensional variables IðtÞ and RC possibly leads to the conclusion that larger output current values can be achieved for smallscale devices in comparison with large-dimension blade-airfoils and or other design conditions. Therefore, it is useful to more carefully analyze the output power of the harvester. Examination of output power and harvester efficiency is presented in this section. Approximately, a harvester configuration with a unit-length apparatus (ℓ ¼ 1 m) could provide a maximum of 300–400 W of power at sustained vibration with rotation amplitude of about 8 and in the range of wind speeds between 10 m/s and 15 m/s. This power value is based on the maximum power that can be extracted from the blowing wind, if three-dimensional flow and other effects are not considered. The power indicated above is compatible with the power of small-scale wind turbines of various configurations and rotor axis orientations, in which either the reference rotor diameter or the scale of the mechanism employed for power generation is about one or few meters (Wood, 2011). More precisely, the input power (Pin ) can be found by combining U  U * beyond the critical flutter speed threshold with the reference dynamic pressure 1=2ρU 2 and the swept area ½2bℓð2αmax Þ, corresponding to the maximum amplitude during one limit-cycle transverse vibration of the airfoil's leeward edge. The leeward edge motion depends on the rotation amplitude, which is a function of time t. Under the hypothesis of small sustained harmonic rotations about O over one limit-cycle period, the expression for the maximum input power is:

Fig. 8. Influence of aeroelastic loads (Fig. 7) and electro-mechanical coupling coefficient (Ψ) on the maximum output current (ιmax) at U ¼ 15 m/s, for a harvester with λRL ¼ 0:75: (a) Type 0, (b) Type 1, (c) Type 2. Results for the basic configurations are indicated by an asterisk symbol; the horizontal-scale axis is not the same in the panels.

11

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

Pin ¼ 2ρU 3 bℓαmax :

(21)

The previous expression imply that the limit-cycle periodic rotations are approximately harmonic with αmax being the maximum rotation or rotation amplitude. The maximum output power (Pout ) is compatibly calculated by examining the relationship between RC and the maximum value of IðtÞ, during limit-cycle post-critical vibration of the blade-airfoil. This leads to the following approximate expression for the maximum instantaneous output power as a function of Imax , the maximum (or the magnitude) of IðtÞ:

Pout  RC ðImax Þ2 ¼

ð2bωα Φe:m:c: ιmax Þ2 : RC

(22)

(a)

As before, the instantaneous output current IðtÞ is represented in a harmonic form, during limit-cycle energy conversion. The quantities in Eqs. (21) and (22), particularly αmax and ιmax , must be determined at a wind speed U  U * . The ratio between Eq. (22) and Eq. (21) leads to the dimensionless maximum output/input power coefficient CP , which is a measure of efficiency in dimensionless form (0 < CP < 1). The quantity CP is obtained, after some simplifications, as:

Pin πΨðιmax Þ kα3 ¼ : Pout 2εαmax 2

CP ¼

(23)

In the previous equation, both αmax and ιmax are evaluated numerically by integration of the reduced-order dynamical system, as described in the previous section. Consequently, the inertia parameter ε and the dimensionless electro-mechanical coupling coefficient Ψ must be initially calculated, e.g., from Table 1. The power coefficient in Eq. (23) is quadratically proportional to the magnitude of the output current ιmax and linearly depends on Ψ (at least explicitly). The results in Eq. (23) are also a function of the mean flow speed U, structural, mechanical and electrical properties of the apparatus. Similarly, the reduced rotational frequency of the blade-airfoil apparatus kα must be evaluated at the specific U value of the actual post-critical scenario. The assessment of the quantities must accordingly be repeated at various U beyond the critical flutter speed threshold to characterize efficiency and operational regimes. Even though an explicit relationship is not evident from Eq. (23), results are also indirectly influenced by impedance λRL (as in Fig. 6) and by the aeroelastic loads (as in Fig. 8). Fig. 9 illustrates the preliminary examination of CP as a function of Ψ for the harvester configurations analyzed in Section 5.2, (a) Type 0, (b) Type 1 and (c) Type 2, at U ¼ 15 m/s with constant impedance λRL ¼ 0:75 and ε. The results are expanded from the three basic configurations, listed in Table 1, which are indicated by an asterisk (star-shaped) symbol. The three panels in Fig. 9 are found from the ιmax vs. Ψ empirical curves shown in Fig. 4 and the corresponding rotation amplitudes αmax (specific values not shown for the sake of brevity). Clearly, low values of CP are found by Eq. (23). One of the reasons is associated with the fact that CP depends on the cubic value of the reduced frequency kα ; the kα values of the three basic configurations, determined at this post-critical speed U, are quite small (between 0.017 for Type-1 device and 0.026 for Type0 device at U ¼ 15 m/s). Undoubtedly, using larger values of kα would improve the performance. For example, a beneficial increase of kα can be achieved by incrementing ωα . This effect must, however, be balanced with a negative influence on critical flutter speed threshold (Section 4). In any case, the CP predictions can be considered as satisfactory for the present study but they are a first approximation of actual power efficiency. More specifically, the CP results are highly sensitive to the selected configuration, with local relative variations about two to three times (for example between the Type-0 and Type-2 cases). Type-2 configuration is more efficient at this wind speed U since the harvester is active for a wider range of Ψ and reaches the largest CP because of ιmax . For comparison purposes, rotor power coefficients, typically found for

(b)

(c) Fig. 9. Preliminary examination of the power coefficient (CP ) vs. electro-mechanical coupling coefficient (Ψ) for the harvester configurations examined in Fig. 4 at U ¼ 15 m/s and with λRL ¼ 0:75: (a) Type 0, (b) Type 1, (c) Type 2. Results for the three basic configurations (Table 1) are indicated by an asterisk symbol; the horizontal-axis scale is not the same in the panels.

rotating wind energy machines, are considered. For example, actual values for large-scale horizontal-axis wind turbines are at most between 0.4 and 0.5 [e.g., Hau (2013)] depending on the number of blades and tip-speed ratios, with the ideal Lanchester–Betz–Joukowsky (van Kuik, 2007) limit of a turbine being equal to 16/27. In the case of vertical-axis wind turbines efficiency is usually lower but achieved with a simpler operational mechanism; power coefficients may be located between 0.15 and 0.40 [e.g., Hau (2013)] depending on the rotor configuration (Darreius, Savonius, etc.) and tip-speed ratio of the blades. The discrepancy between these power coefficients and the CP values in Eq. (23) is apparent. The results in Fig. 9 and the use of Eq. (23) are, however, not completely conclusive since several simplifications have been introduced to enable the calculation of the CP . For example, Eq. (23) tends to suggest that efficiency is inversely proportional to the rotational angle (the pitching amplitude), which is perhaps counterintuitive. More careful examination of the output power is possibly needed through variation of 12

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

other quantities (such as U) or properties (λRL , ε, etc). Numerical integration and automated (i.e., empirical) extraction of the ιmax may also have had an impact on the current results. Similar conclusions, predicting a low output power, were previously drawn (Caracoglia, 2010) and were primarily attributed to the use of linear aeroelastic loads during post-critical motion in Eqs. (5) and (6).

Abdelkefi, A., Yan, Z., Hajj, M.R., 2013. Modeling and nonlinear analysis of piezoelectric energy harvesting from transverse galloping. Smart Mater. Struct. 22 (2), 025016 (1–9). Abdelkefi, A., Yan, Z., Hajj, M.R., 2014. Performance analysis of galloping-based piezoaeroelastic energy harvesters with different cross-section geometries. J. Intell. Mater. Syst. Struct. 25 (2), 246–256. Ahmadi, G., 1978. Aeroelastic wind energy converter. Energy Convers. 18, 115–120. Ahmadi, G., 1979. An oscillatory wind energy converter. Wind Eng. 3 (3), 207–215. Argentina, M., Mahadevan, L., 2005. Fluid-flow-induced flutter of a flag. PNAS - Proc. Natl. Acad. Sci. U. S. A. 102 (6), 1829–1834. Bernitsas, M.M., Ben-Simon, Y., Raghavan, K., Garcia, E.M.H., 2009. The VIVACE converter: model tests at high damping and Reynolds number around 105. J. Offshore Mech. Arct. Eng. 131 (1), 011102. Bernitsas, M.M., Raghavan, K., Ben-Simon, Y., Garcia, E.M.H., 2008. VIVACE (Vortex Induced Vibration Aquatic Clean Energy): a new concept in generation of clean and renewable energy from fluid flow. J. Offshore Mech. Arct. Eng. 130 (4), 041101. Biscarini, G., Petrini, F., Gkoumas, K., Bontempi, F., 2016. Piezoelectric EH from flowinduced structural vibrations. In: Proceedings of the XIV Conference of the Italian Association for Wind Engineering (IN-VENTO-2016), Terni, Italy. Bisplinghoff, R.L., Ashley, H., Halfman, R.L., 1955. Aeroelasticity. Dover Publications Inc., Mineola, NY, USA. Caracoglia, L., 2010. Feasibility assessment of a leading-edge-flutter wind power generator. J. Wind Eng. Ind. Aerod. 98 (10–11), 679–686. Caracoglia, L., 2017. Numerical investigations on the operational regimes of a torsionalflutter-based wind harvester. Procedia Eng. 199, 3434–3439. Casadei, F., Bertoldi, K., 2014. Wave propagation in beams with periodic arrays of airfoilshaped resonating units. J. Sound Vib. 333, 6532–6547. Dessi, D., Mazzocconi, S., 2015. Aeroelastic behavior of a flag in ground effect. J. Fluids Struct. 55, 303–323. Donoso, G., Ladera, C.L., Martin, P., 2009. Magnet fall inside a conductive pipe: motion and the role of the pipe wall thickness. Eur. J. Phys. 30, 855–869. Dowell, E.H., Clark, R., Cox, D., Curtiss Jr., H.C., Edwards, J.W., Hall, K.,C., Peters, D.A., Scanlan, R.H., Simiu, E., Sisto, F., Strganac, T.W., 2004. A Modern Course in Aeroelasticity - Fourth Revised and Enlarged Edition. Kluwer Academic Publishers, Dordrecht, The Netherlands. Dunn, P., Dugundji, J., 1992. Nonlinear stall flutter and divergence analysis of cantilevered graphite/epoxy wings. AIAA J. 30 (1), 153–162. Dunnmon, J.A., Stanton, S.C., Mann, B.P., Dowell, E.H., 2011. Power extraction from aeroelastic limit cycle oscillations. J. Fluids Struct. 27 (8), 1182–1198. Farthing, S.P., 2009. Vertical axis wind turbine induced velocity vector theory. Proc. Institution Mech. Eng. - Part A J. Power Energy 223 (2), 103–114. Gkoumas, K., Petrini, F., Bontempi, F., 2017. Piezoelectric vibration energy harvesting from airflow in HVAC (Heating Ventilation and Air Conditioning) systems. Procedia Eng. 199, 3444–3449. Hau, E., 2013. Wind Turbines: Fundamentals, Technologies, Applications, Economics, third ed. Springer-Verlag, Berlin/Heidelberg, Germany. Jones, R.T., 1939. The Unsteady Lift of a Finite Wing. Technical Note 682. NACA. Kakkavas, C., 1998. Computational Investigation of Subsonic Torsional Airfoil Flutter. MS Thesis. Naval Postgraduate School, Monterrey, California, USA. Kazmierski, T.J., Beeby, S., 2011. Energy Harvesting Systems: Principles, Modeling and Applications. Springer, New York, New York, USA. Kwon, S.-D., Lee, H., Lee, S., 2011. Wind energy harvesting from flutter. In: Proceedings of the 13th International Conference on Wind Engineering (ICWE-13), Amsterdam, NL, Paper No. 071 (CD-ROM). Kwon, S.-D., Park, J., Law, K., 2013. Electromagnetic energy harvester with repulsively stacked multilayer magnets for low frequency vibrations. Smart Mater. Struct. 22 (5), 055007. Manwell, J.F., McGowan, J.G., Rogers, A.L., 2010. Wind Energy Explained. Theory, Design and Application, second ed. John Wiley and Sons, Chichester, UK. Matsumoto, M., Okubo, K., Keisuke, M., Ito, Y., 2006. Fundamental study on the efficiency of power generation system by use of the flutter instability. In: Proceedings of the 2006 ASME Pressure Vessels and Piping Division Conference, Vancouver, British Columbia, Canada. ASME Paper No. PVP2006-ICPVT-11–93773. Mehmood, A., Abdelkefi, A., Hajj, M.R., Nayfeh, A.H., Akhtar, I., Nuhait, A.O., 2013. Piezoelectric energy harvesting from vortex-induced vibrations of circular cylinder. J. Sound Vib. 332 (19), 4656–4667. Nayfeh, A.H., Mook, D.T., 1995. Nonlinear Oscillations. John Wiley and Sons, Hoboken, New Jersey, USA. Petrini, F., Gkoumas, K., Bontempi, F., 2014. Piezoelectric energy harvesting under air flow excitation. In: Carassale, L., Repetto, M.P. (Eds.), Wind Engineering in Italy Proceedings of the 13th Conference of the Italian Association for Wind Engineering (In-Vento 2014), June 22-25, 2014. Genova University Press (electronic proceedings), Genoa, Italy. Pigolotti, L., Mannini, C., Bartoli, G., 2016. Experimental investigations on a flat plate equipped with porous screens undergoing classical flutter oscillations. In: Proceedings of the XIV Conference of the Italian Association for Wind Engineering (IN-VENTO-2016), Terni, Italy. Pigolotti, L., Mannini, C., Bartoli, G., Thiele, K., 2017. Critical and post-critical behaviour of two-degree-of-freedom flutter-based generators. J. Sound Vib. 404, 116–140. Priya, S., Inman, D.J., 2009. Energy Harvesting Technologies. Springer Science, New York, NY, USA. Scanlan, R.H., Beliveau, J.G., Budlong, K.S., 1974. Indicial aerodynamic functions for bridge decks. J. Eng. Mech. ASCE 100 (EM4), 657–673. Scanlan, R.H., Rosenbaum, R., 1968. Introduction to the Study of Aircraft Vibration and Flutter. Dover Publications, New York, New York, USA.

7. Conclusions and outlook The analysis of variability in the aeroelastic loads has revealed a nonnegligible influence of the model used for describing the unsteady aerodynamics on the output current. Variations of the order of thirty to forty percent, depending on the type of apparatus, have been noted despite the small modifications in the parameters of the load functions. This remark reaffirms the crucial role of unsteady aerodynamics in relation to this specific type of harvester. Moreover, it is observed that these variations have been found without introducing any modification to the values of steady aerodynamic torque (and lift) coefficients. It is therefore possible that a variation in the critical flutter speed of the linear model may be responsible for this apparently remarkable result. Unfortunately, Eqs. (18) and (19) can only provide an approximate indication; they are not fully appropriate to predict flutter under the syntheticallymodified load cases. The use of a nonlinear aeroelastic load model is needed to overcome these issues. This model could more accurately simulate fluid-induced “damping” effects (briefly discussed in the paper only). In addition, studies should be carried out to evaluate the importance of three-dimensional flow effects for apparatuses with small length/chord ratios; preliminary and partial results (Caracoglia, 2017) are only available. Limitations of the proposed apparatus are also present. They include, for example, a potentially large cut-in wind speed for devices with large chords b and operational frequencies (above 1 Hz). This problem is confirmed by the examination of Eq. (20), which recognizes some restrictions dictated by the flutter type and, on occasion, possibly large and unsuitable wind speed threshold. This is a supplementary constraint that would make the apparatus less attractive under variable operational conditions, if the goal is to generate energy daily with sporadic occurrences of relevant sustained winds. The issues of output current, power and efficiency have also been analyzed. Numerical results indicate that, although the performance coefficient is lower than the one typically attributed to conventional wind energy systems, preliminary and promising values have been found. The power coefficient has been calculated theoretically using basic conversion laws and physical principles. Although this approach can be used in the case of wind turbines, it is perhaps not fully applicable to the proposed harvester since the solution of the post-critical stage of conversion has been carried out numerically. More studies are needed (especially experimental) to fully address the calculation of both induced current and output power, which need to be evaluated for a wider range of wind speeds and system properties. The use of other dimensionless groups or variables, different from the ones used in Eq. (23) and inspired by the form of Eqs. (5) and (6), may be more adequate. Finally, in-situ deployment of the apparatus requires further investigation, after experimental examination along with prior verification in wind tunnel to improve the estimation of efficiency. Despite several technological challenges, the present study suggests that the exploitation of wind energy through the harvesting of flowinduced vibration is an interesting practical subject that deserves careful consideration. References Abdelkefi, A., Nayfeh, A.H., Hajj, M.R., 2012a. Design of piezoaeroelastic energy harvesters. Nonlinear Dyn. 68 (4), 519–530. Abdelkefi, A., Nayfeh, A.H., Hajj, M.R., Najar, F., 2012b. Energy harvesting from a multifrequency response of a tuned bending-torsion system. Smart Mater. Struct. 21 (7), 075029. 13

L. Caracoglia

Journal of Wind Engineering & Industrial Aerodynamics xxx (2017) 1–14

Scanlan, R.H., Tomko, J.J., 1971. Airfoil and bridge deck flutter derivatives. J. Eng. Mech. ASCE 97 (EM6), 1717–1737. Sewell, G., 2005. The Numerical Solution of Ordinary and Partial Differential Equations, second ed. John Wiley and Sons, Hoboken, New Jersey, USA. Shimizu, E., Isogai, K., Obayashi, S., 2008. Multiobjective design study of a flapping wind power generator. J. Fluids Eng. ASME 130 (2), 021104. Stephen, N.G., 2006. On energy harvesting from ambient vibration. J. Sound Vib. 293 (1–2), 409–425. Tang, L., Païdoussis, M.P., Jiang, J., 2009. Cantilevered flexible plates in axial flow: energy transfer and the concept of flutter-mill. J. Sound Vib. 326 (1–2), 263–276. Theodorsen, T., 1935. General Theory of Aerodynamic Instability and the Mechanism of Flutter. Technical Report 496. National Advisory Committee for Aeronautics, Washington, DC, USA. Tomasini, G., Giappino, S., 2016. Galloping-based piezo-aeroelastic energy harvester: numerical models and wind tunnel experimental tests. In: Proceedings of the XIV

Conference of the Italian Association for Wind Engineering (IN-VENTO-2016), Terni, Italy. Tran, C.T., Petot, D., 1981. Semi-empirical model for the dynamic stall of airfoils in view of the application to the calculation of responses of a helicopter blade in forward flight. Vertica 5 (1), 35–53. van Kuik, G.A.M., 2007. The Lanchester–Betz–Joukowsky limit. Wind Energy 10 (3), 289–291. Wood, D., 2011. Small Wind Turbines: Analysis, Design, and Application. SpringerVerlag, London, UK. Young, J., Lai, J.C.S., Platzer, M.F., 2014. A review of progress and challenges in flapping foil power generation. Prog. Aerosp. Sci. 67, 2–28. Zhu, J., Zhang, W., 2015. Coupled analysis of multi-impact energy harvesting from lowfrequency wind induced vibrations. Smart Mater. Struct. 24 (4), 045007. Zhu, Q., 2011. Optimal frequency for flow energy harvesting of a flapping foil. J. Fluid Mech. 675, 495–517.

14