Abstract

Wind energy harvesters are emerging as a viable alternative to standard, large horizontal-axis wind turbines. This study continues a recent investigation on the operational features of a torsional-flutter-based apparatus, proposed by the author to extract wind energy. The apparatus is composed of a non-deformable, flapping blade-airfoil. A nonlinear torsional spring mechanism, either simulated by a Duffing model or a hybrid Duffing–van der Pol model, installed at equally spaced supports, enables limit-cycle, post-critical vibration. To enhance the output power, stochastic resonance principles are invoked through a novel, negative stiffness mechanism that is coupled to the eddy current device for energy conversion. The output power is explored by numerically solving the stochastic differential equation of the model, accounting for incoming flow turbulence. Three main harvester types with variable configuration are examined; the chord length of the blade-airfoil, used for energy harvesting, varies between 0.5 and 1 m; the spanwise-length-to-chord aspect ratio is four. The flapping frequency varies between 0.10 and 0.25 Hz. The study demonstrates that exploitation of negative stiffness mechanism can improve the performance of the harvester.

1 Introduction

Several harvesting devices [1] have been proposed to exploit wind energy in the low wind speed range and for small-scale applications. These wind energy technologies are triggered by aeroelastic instability [2,3]. Most existing harvesters are designed as micro-mechanical units and used for recharging miniature sensors [4]. By contrast, “meso-scale” harvesting technologies are still an active research field. For example, a hybrid device has been installed on a highway bridge and utilizes two different energy sources (natural wind flow and structural vibrations) [5]. Furthermore, vortex-induced vibrations of multi-unit circular cylinders in water flows [6] have been considered. Harvesting devices, triggered by water axial-flow instabilities [7], or by flutter of an inverted and flexible flag [8], or by vortex-induced vibration of cylinders inside ventilation ducts [9,10], have been studied. The exploitation of coupled, bending-torsion vibration of streamlined bodies, i.e., the “flutter mill,” has been proposed for energy harvesting in wind flows [2,3,6,9]; several flapping foil mechanisms have been derived from this concept [11]. All of the above share the common task of scavenging the kinetic energy of the flow at low wind speeds.

Along this line of research, this paper continues recent investigations on the use of a novel wind energy technology, which considers a purely torsional-flutter harvester [12]. Contrary to the previously mentioned devices, the apparatus has a simple operational mechanism that is ideal for implementation; it is also designed to supplement energy to one or few residential housing units.

The proposed apparatus is inspired by seminal works [13,14]. It is composed of a rigid, flapping blade-airfoil that rotates about a pivot axis, positioned upwind with respect to the airfoil’s aerodynamic center. The blade-airfoil has a chord length equal to about 1 m and a spanwise length between 4 and 10 times the chord (i.e., aspect ratio (AR)). Both a nonlinear torsional spring mechanism (Duffing model) and a hybrid one (Duffing and van-del-Pol model), installed at equally spaced spanwise supports close to the pivot position, have been considered in previous studies to convert the kinetic energy of the flow to electrical energy by eddy currents.

Stochastic modeling and analysis of the harvester are used to evaluate the output power [12]. The state vector of the model includes physical states (flapping angle and angular velocity) and aeroelastic torque unsteady flow states, corrected for three-dimensional flow effects. Turbulence is also included in the stochastic dynamic equation. A wide range of stationary turbulent winds can be simulated, typical of complex urban flow environments.

In this study, an enhancement is explored by considering the principles of stochastic resonance, i.e., the aptitude of an external noise (turbulence) to amplify the unstable flapping. This condition is achieved by adding a negative stiffness mechanism, embedded within the transmission shaft that is connected to the eddy current power system. The negative stiffness destabilizes the torsional restoring torque provided by a positive, cubic nonlinear Duffing model. The conceptual idea is also borrowed from the bi-modal Duffing oscillators excited by stationary white noise, for which several investigations and solutions are available [15].

The negative stiffness mechanism is installed at equally spaced supports; it is used to enhance the amplitude of transverse, post-critical vibration of a permanent magnet within the electro-magnetic coil of the power system. Conversion to electrical power is in fact warranted by an eddy current power system with multi-loop magnetic coil and a translating permanent magnet [16]. Energy is stored in a battery for future use.

In this initial exploration, example configurations are considered with an adjustable negative stiffness. The rotation, pivot axis position is placed at the windward edge of the blade-airfoil. Stationary wind flows are exclusively examined. Turbulence intensity equal to 10% is used to investigate realistic types of wind flow. Stability is studied by solving the stochastic differential equations in time domain, by ensemble averaging and Monte-Carlo sampling.

2 Theoretical Derivations

The apparatus (Fig. 1) is composed of a flapping blade-airfoil with dimensions b (half-chord width) and (spanwise, longitudinal length). The main degree-of-freedom (DOF) of the model is the torsional rotation α about pivot “O.” Aeroelastic torque is modeled by standard aeroelastic formulation for small angles of attack, corrected for three-dimensional flow effects that depend on the aspect ratio of the blade.

Fig. 1
Schematics of the torsional-flutter harvester—cross-sectional view on the XY horizontal plane
Fig. 1
Schematics of the torsional-flutter harvester—cross-sectional view on the XY horizontal plane
Close modal

The load model also replicates the main features of synoptic, stationary turbulent winds and combines mean outflow wind speed U with along-wind stationary turbulence u. Since the blade longitudinal plane is oriented on a vertical plane (XZ), gravity is not relevant to energy conversion.

The external torque is about pivot O in Fig. 1. Noting the distance ab from the mid-chord point, the pivot axis position can vary between the leading edge (a=1) and the one-quarter chord position (a=0.75) [12].

2.1 Non-Rigid Transmission to Electrical Power Circuit

2.1.1 Foreword.

Contrary to previous studies, a non-rigid transmission system is employed to transfer the flow energy to electrical energy. A schematic diagram of the energy conversion system and secondary power circuit is depicted in Fig. 2. The foremost leeward edge of the flapping airfoil (node A) undergoes rotation about the windward-edge pivot “O” (Fig. 1). A rigid linkage AB is connected by a hinge at B to a massless shaft BC. The BC shaft is equipped with a permanent magnet at the upper end (toward node C in Fig. 2) that slides inside a multi-loop coil. This mechanism is used to convert the rotation to transverse translation along the y axis. The conversion is based on eddy current, induced by the permanent magnet that slides inside the multi-loop coil. A frictionless collar is used to guide the magnet along a trajectory that, for small flapping amplitudes α, is approximately directed along the y axis in the figure. This type of multi-loop coil system is inspired by an electro-magnetic harvester, triggered by random vibration sources [16].

Fig. 2
Schematics of the torsional-flutter harvester—cross-sectional view on the XY horizontal plane
Fig. 2
Schematics of the torsional-flutter harvester—cross-sectional view on the XY horizontal plane
Close modal

2.1.2 Dynamic Equilibrium at Node “C”.

The equilibrium equation at node C of the transmission shaft BC is examined in detail (Fig. 2). In Fig. 2, p is the absolute translation at the massless linkage AB; pC,r is the relative displacement at “C.” A non-rigid connector is inserted between nodes “B” and “C” of the transmission shaft that connects the linkage AB to the electro-magnetic induction converter, i.e., the magnet that slides through a frictionless collar into the multi-loop coil. The moving magnet introduces magnetic induction and interacts with the “BC” shaft, translating in the y direction inside the winding coil (Fig. 2). The quantity I is the eddy-induced output current of the power system and Φ(e.m.c.) is the electro-mechanical coupling coefficient in units of Newton/Ampère [12].

The dynamic equilibrium equation of the translating mass, in the y direction and in time domain (t), reads
(1)
where mC denotes the mass of the permanent magnet and shaft BC; kBC<0 is the negative stiffness of the non-rigid transmission; pC,r is the relative transverse displacement between B and C, and f(e.m.)(t)=Φ(e.m.c.)I(t) is the time-dependent electromotive force that obeys the Faraday’s Law [1,16,17].

If the mass mC0 is negligible, the previous equation can be simplified as kCpC,r(t)=f(e.m.)(t)=FAB, which coincides with the internal constraint force FAB transferred to the primary system through linkage AB in the y direction. For small flapping angles about O, we note that p(t)(1a)bα(t).

2.1.3 Electrical Power System Equation.

The eddy current equation of the power circuit with electro-mechanical coupling is found by Faraday’s law of magnetic induction [1,12,16,17]:
(2)
In Eq. (2)RC is the resistance and LC is the inductance of the eddy current electrical circuit. The previous equation can be simplified by noting that, since kCpC,r(t)=f(e.m.)(t)=Φ(e.m.c.)I(t) and mC0, static condensation of the variable pC,r(t) is carried out. With p(t)(1a)bα(t), Eq. (2) becomes
(3)
Equation (3) is converted to dimensionless form by designating λRL=RC/(ωαLC) as the generalized impedance of the circuit. After normalizing the linear stiffness coefficient of the coupling link as kBC=χBCkαb2 (with χBC<0), the “soft” coupling parameter ΔBC is defined, with
(4)
In Eq. (4), the dimensionless electro-mechanical coupling coefficient of the “rigid shaft” is [12]
(5)
where I0α is the total polar mass moment of inertia of the flapping blade-airfoil about pivot O.
Finally, Eq. (6) is found by setting: τ=tωα as the dimensionless time, ωα as the angular frequency of the linear mechanical system and ι=I(t)/i0 as the dimensionless output current with reference current intensity i0=Φ(e.m.c.)(1a)ωα/RC (A or Ampère). This yields
(6)

For a rigid shaft, it is noted that χBC+ and ΔBC1.

2.2 Dynamic Equilibrium Equation of the Mechanical System.

The physical states of the primary, mechanical system are torsional rotation α and derivative dα/dτ with respect to dimensionless time τ. The dynamic equilibrium equation of the flapping rotation is
(7)

In Eq. (7), the external aeroelastic torque M0z is about the pivot O, and I0α is the total polar mass moment of inertia about O. Structural damping is simulated through a linear term in Eq. (7) with damping ratio ζα; HNL(α,dα/dτ) is a function of α and dα/dτ that encloses the nonlinear restoring and damping torque mechanisms, and that enables power conversion and limit-cycle oscillation [12,18].

The moment is also composed of the electromotive torque about O
(8)

In Eq. (8), ι(τ) is the dimensionless current, i0 is the reference current intensity, and Φ(e.m.c.) is the electro-mechanical coupling coefficient, as derived in Sec. 2.1.3.

Duffing model: HNL(α,dα/dτ)=κα3, i.e., the nonlinear restoring force effect is simulated by the term κα3, with κ>0 being a suitable dimensionless constant.

Hybrid Duffing–van Der Pol model: The equilibrium incorporates both Duffing and van der Pol models [19] to improve energy conversion. The modified nonlinear restoring torque function reads HNL(α,dα/dτ)=[κα32ζαγα2dα/dτ]. Besides parameter κ, there is a nonlinear damping coefficient γ>0, which models the self-limiting feature of the negative damping mechanism proportionally to the linear damping term 2ζα.

3 Stochastic Differential Equations

3.1 Random Turbulence Wind Flows.

The aeroelastic torque is proportional to the mean flow speed U, contaminated by along-wind stationary turbulence u. The flapping foil of the apparatus has a reference diagonal dimension (2b) with 2b being the chord length, the spanwise width, and AR=/b the geometric aspect ratio. Therefore, the main transitory feature of the flow that influences the stability is the instantaneous magnitude of the flow speed.

The along-wind turbulence component u(t) depends on dimensional time t=τωα or, equivalently, τ=t/ωα. The turbulence influences the dynamic pressure and load. The dynamic pressure is also proportional to constant flow speed U. The total, instantaneous dynamic pressure is proportional to (m/s)2:
(9)

In Eq. (9), turbulence is random, stationary, and normalized as u^(τ)=u(τ)/U. The properties of the stationary process are represented, without any loss of generality, as a Gaussian white noise process, simulated in Sec. 3.4 by a Wiener process [20] of independent Gaussian increments with standard deviation σu^.

Finally, since the diagonal reference dimension of the blade-airfoil (2b) is small compared to the integral turbulence length scales, the gust load is fully coherent across the apparatus (Fig. 1).

3.2 Aeroelastic Torque With Random Perturbations.

Mean aerodynamic forces are zero since static lift force is zero at α=0 due to blade symmetry. The aeroelastic torque M0z depends on the constant mean speed U [12] with mean direction (orientation of the blade-airfoil) parallel to the x axis in Fig. 1. Aeroelastic torque is modeled in the time domain by Wagner’s indicial function formulation [21,22], later corrected for lift and torque reduction due to three-dimensional flow effects. The Wagner’s indicial function of the load is [23]
(10)
with reduced frequency kα=ωαb/U and suitable load parameters: c1=0.165, d1=0.0455, c2=0.335, d2=0.3 [23] for an idealized, symmetric NACA0012 airfoil section.

Furthermore, random load perturbation is introduced by replacing the deterministic quantity d2=0.3 with a random, time-dependent d2=(d2,m+δ2(τ)); d2,m=0.3 and δ2(τ) is a zero-mean, Gaussian perturbation; δ2(τ) also accounts for load measurement error and modeling simplifications. It is assumed that the unsteady torque features c1,d1,c2 are not affected either by turbulence or any other perturbation.

3.3 Dynamic Equilibrium With Leading Edge Flapping Pivot.

If the rotation axis is at the leading edge (a=1), the dynamic equation of the flapping angle from Eq. (7) is [12]
(11)

In Eq. (11), Φ0=Φ(0)=0.5; νae,1(τ), μae,1(τ), νae,2(τ), and μae,2(τ) are four time-dependent aeroelastic states; ϵ=πρb4/I0α is the inertia parameter [12], and ρ is the air density. Parameter η3D=AR/(AR+2) accounts for three-dimensional flow effects on lift and torque [22].

On the left-hand side of Eq. (11), the nonlinear function HNL(α,dα/dτ) is again used to differentiate between Duffing model and Hybrid Duffing–van der Pol one, as in Eq. (7).

3.4 State-Space Vector Equation.

Equation (11) is combined with Eq. (6) to form a state-space model, composed of seven nonlinear, coupled electro-mechanical equations. Since the turbulence variable u^ is a Gaussian process, u^ dependency is represented by a Wiener noise process Bu^(τ) of independent unit-variance Gaussian increments. Similarly, aeroelastic load perturbation δ2(τ) depends on a second scalar Wiener noise BΔ2(τ). The final vector equation is a stochastic differential equation of the Itô-type [24]:
(12)

Wem=[α,dα/dτ,νae,1,νae,2,μae,1,μae,2,ι]T in Eq. (12) is the state vector; Wem includes both physical, aeroelastic states and dimensionless output current ι.

In Eq. (12), Wiener noise Bu^(τ) separately addresses turbulence perturbation from BΔ2(τ), used for load perturbation. Quantity qem,NL is a nonlinear vector-function; tNL,u^ is a nonlinear turbulence diffusion vector-function that depends on turbulence intensity σu^. QL,Δ2 is a constant, diffusion matrix that controls the load perturbation and depends on the standard deviation of δ2(τ), σd2. Both qem,NL and tNL,u^ are nonlinear; they are two 7×1 vector-functions.The non-zero elements of the 7×7QL,Δ2 matrix are two only. The Wong-Zakai [25] correction terms are introduced in both QL,Δ2 and tNL,u^. Derivation of qem,NL, QL,Δ2, and tNL,u^(Wem) is omitted for the sake of brevity. A detailed description of these equations may be found in previous works [12,18].

Finally, Eq. (12) must be solved numerically [26] with appropriate initial conditions, imposed on the random Wem at initial time τ=0.

3.5 Mean-Square Stability.

The stability analysis is studied numerically. It relies on the use of moment Lyapunov exponent (MLE) [27], usually considered in the context of wind-excited nonlinear systems [28,29]. The MLE measures the propensity of the system’s slow dynamics to asymptotically exhibit a diverging, unstable oscillatory trend. First, Eq. (12) is solved in weak form, i.e., by numerical integration that is repeated several times by Monte-Carlo sampling [27]. Second, mean-square stability requires the evaluation of the second MLE of the sub-vector Υ=[α,dα/dτ,ι]T:
(13)
where E[.] is the expectation operator applied to the Euclidean vector norm; Υ(τl) is evaluated at time τl, which is sufficiently large to approximate the true MLE value, i.e., limτl+ΛΥ(2) [27]. The vector Υ in Eq. (13) also includes the output current ι to evaluate the effect on operational conditions and energy conversion.

As explained by Xie [27,30],limτl+ΛΥ(2) is a measure of the slow dynamics of a nonlinear system; it can be interpreted as the role of damping in the linear deterministic case. Specifically, ΛΥ(2)<0 indicates that the variance of the 2-norm of zero-mean vector Υ in Eq. (13) vanishes as time tends to infinity, i.e., the system is asymptotically stable in the mean squares.

From a practical perspective, E[.] in Eq. (13) is approximated by arithmetic averages after collecting a sample of 200 repeated, numerical solutions of Eq. (12). Furthermore, Eq. (13) is numerically found and examined at τ>300, i.e., after about 47 equivalent periods of aeroelastic flapping and, if the Υ<0 is noted, the harvester is labeled as stable, i.e., not engaged in energy conversion.

3.6 Output Power Estimation.

As described in Sec. 3.5, Eq. (12) is numerically solved in weak form over time τ0 with appropriate initial conditions, imposed on the probability distribution of Wem. If MLE analysis through Eq. (13) detects an unstable regime and post-critical operational conditions, the output power at sustained flapping can also be evaluated numerically. The output power is a random process; its expectation is [18]
(14)
with Ψ found from Eq. (5); E[ι2(τ)] is the variance of ι.

4 Description of the Modeled Units and Flow Properties

Simulations examine various configurations of the apparatus, with rotation axis pivot at the apex of the blade-airfoil (a=1). Three “types” are selected from Refs. [12,18], with the main properties described in Table 1; in this table, quantities such as angular frequency and damping ratio must be interpreted as the properties of the linearized dynamic equation of the harvester.

Table 1

Main harvester and flow properties

I0αωα2πζαARσu^σd2
Typeb (m)((kg×m2)/m)(Hz)(%)(–)(%)(–)
00.25200.250.254100.07
10.25400.250.304100.07
20.503000.100.304100.07
I0αωα2πζαARσu^σd2
Typeb (m)((kg×m2)/m)(Hz)(%)(–)(%)(–)
00.25200.250.254100.07
10.25400.250.304100.07
20.503000.100.304100.07

The nonlinear stiffness parameter in the Duffing model is constant and set to κ=100, irrespective of the type. The nonlinear damping parameter γ=0.5 is used in the hybrid Duffing–van der Pol model. Coupling with the power circuit is achieved by setting Ψ=0.01 and λRL=0.75 [12].

Initial conditions are imposed on the random vector Wem at τ=0 by considering an initial flapping angular motion at τ=0 that triggers the instability [12]. The random initial amplitude of α is set to zero mean and standard deviation equal to 2 deg, coincident with small, realistic angular deviations from the static equilibrium.

In all the simulations, the wind flow is turbulent with intensity σu^=10% and load perturbation with standard deviation σd2=0.07.

5 Numerical Results

Numerical simulations are executed in the range 0<Umax20m/s only, since the efficacy of the harvester at low wind speeds is desirable.

In the next four figures, the legend labels “Ty.0,” “Ty.1,” and “Ty.2,” respectively, refer to type-0, type-1, and type-2 apparatuses, described in Table 1.

5.1 Operational Regimes.

Figure 3 illustrates the reference simulation result for the linear model (κ=1E7,γ=0) found for turbulent wind flows of intensity 10% and moderate load perturbation. The flexible transmission link has a moderate negative stiffness with χBC=0.1.

Fig. 3
ΛΥ(2)versus time τ for the linear model and apparatus at mean flow speeds: (a)U=12.0m/s and (b)U=16.6m/s
Fig. 3
ΛΥ(2)versus time τ for the linear model and apparatus at mean flow speeds: (a)U=12.0m/s and (b)U=16.6m/s
Close modal

Examining the linear dynamic model is useful to evaluate the incipient instability threshold, in a linear system, divergent vibrations are only possible and no limit cycle exists. The figure demonstrates that both type-2 and type-0 apparatuses, with transmission shaft equipped with negative stiffness mechanism, tend to be unstable at low wind speeds, e.g., U=12m/s.

By contrast, Fig. 4 examines the Duffing case with nonlinear stiffness coefficient κ=100 and negative stiffness apparatus. The system is stable at U=12.0m/s for all the wind speeds examined. The type-2 case diverges due to numerical integration instability [18]. Type-2 apparatus is unstable and therefore operational at U>12.0. Stability, however, may still correspond to operational wind conditions, since limit-cycle vibration entails a periodic response (zero damping).

Fig. 4
ΛΥ(2)versus time τ for the nonlinear Duffing model and apparatus at mean flow speeds: (a)U=12.0m/s and (b)U=14.8m/s
Fig. 4
ΛΥ(2)versus time τ for the nonlinear Duffing model and apparatus at mean flow speeds: (a)U=12.0m/s and (b)U=14.8m/s
Close modal

Finally, Fig. 5 examines the hybrid model at U=19.8m/s with both nonlinear damping and nonlinear stiffness in the primary system (κ=100, γ=0.5); the negative stiffness effect through the passive connector is again set to χBC=0.1. Turbulence of intensity σu^=10% and load perturbation σd2=0.07 are the same as before.

Fig. 5
ΛΥ(2)versus time τ for the hybrid Duffing–van der Pol model at a mean flow speed U=19.8m/s
Fig. 5
ΛΥ(2)versus time τ for the hybrid Duffing–van der Pol model at a mean flow speed U=19.8m/s
Close modal

Again, an apparently stable behavior is observed at high wind speeds U<20m/s, suggesting that the apparatus may require re-designing the transmission mechanism between flapping blade tip and permanent magnet (Fig. 2), if the hybrid Duffing–van der Pol apparatus is employed.

5.2 Output Power.

Figure 6 presents preliminary output power results at various mean flow speeds U for an apparatus equipped with nonlinear Duffing restoring torque mechanism and negative stiffness transmission shaft using a soft coupling parameter χBC=0.1. This figure transfers the information, originally included in Fig. 4 and discussed in the previous sub-section, to energy conversion. In the two figure panels, the mean output power is estimated numerically by Eq. (14) in the time interval 0τ200. The power estimation for U=14.8m/s in the panel of Fig. 6(b) is interrupted at τ70, since numerical integration issues arise (also refer to Fig. 4(b)).

Fig. 6
Mean output power E[Pout(τ)]versus time τ for the nonlinear Duffing model and apparatus at mean flow speeds: (a)U=12.0m/s and (b)U=14.8m/s
Fig. 6
Mean output power E[Pout(τ)]versus time τ for the nonlinear Duffing model and apparatus at mean flow speeds: (a)U=12.0m/s and (b)U=14.8m/s
Close modal

This figure demonstrates that type-2 apparatus is active and produces energy that might be stored in a battery for future re-use [12]. Nevertheless, the available power is rather small and possibly unsuitable to practical design and implementation. Additional investigations are necessary to optimize the negative stiffness transmission shaft with soft coupling parameter.

6 Concluding Remarks

An initial inspection of the numerical results suggested that the use of a hybrid model with a nonlinear torsional stiffness term and a nonlinear damping term may produce a stable system at all wind speeds. Therefore, any non-linearity, installed in the primary mechanism that enables the airfoil flapping, may neutralize the torsional flutter and, consequently, be detrimental to the energy conversion.

To enhance the harvester vibrations, a flexible link was proposed, exploiting the properties of a negative stiffness mechanism connected to the eddy current power device. The study demonstrates that the output power can be increased, approximately by a factor of two for a type-2 apparatus, and that the performance can be remarkably improved.

In spite of the promising outcome, future studies are still necessary to evaluate the output power at other flow regimes, induced by the flapping blade-airfoil and enhanced by the negative stiffness mechanism. Energy conversion may also require further improvements. For instance, since the added mass of the secondary system connected to the negative spring mechanism has been neglected, enhancements should be possible by specifically including this quantity. Consequently, the ratio between mass moments of inertia of the flapping foil and the secondary power device system could possibly be adjusted to achieve a better performance.

Footnote

2

This work was completed while the author was on sabbatical leave as a visiting professor (January–April 2024) in the Department of Civil, Environmental and Mechanical Engineering, University of Trento (Italy), and a visiting professor (June 2024) in the Department of Civil and Environmental Engineering, University of Perugia (Italy).

Acknowledgment

This work was completed while L. Caracoglia was on sabbatical leave in early 2024; L. Caracoglia gratefully acknowledges the support of the Universities of Trento and Perugia (Italy) during his sabbatical leave in 2024. This material is based in part upon work supported by the National Science Foundation (NSF) of the United States of America, Award CMMI-2020063. Any opinions, findings, and conclusions or recommendations are those of the author and do not necessarily reflect the views of the NSF.

Conflict of Interest

The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data Availability Statement

Data and models that support the findings of this study are available from the corresponding author upon reasonable request.

Nomenclature

a =

dimensionless position of the pivot “O”

b =

half-chord width of the blade-airfoil

=

spanwise length of the blade-airfoil

p =

absolute transverse displacement of DOF B

t =

time (s)

u =

along-wind stationary turbulence (m/s)

z =

vertical axis coordinate

u^ =

normalized along-wind stationary turbulence

I =

output current of the power system (A)

U =

mean outflow wind speed (m/s)

ci =

parameters of the Φ function (i=1,2)

d2,m =

mean value of random exponent parameter d2

di =

parameters of Φ function (i=1,2)

kBC =

dimensional axial stiffness, flexible BC shaft

kα =

reduced frequency, one-DOF flapping foil

mC =

lumped mass of the transmission shaft (kg)

pC,r =

relative transverse displacement BC (transmission shaft)

qem,NL =

nonlinear drift vector-function

tNL,u^ =

nonlinear turbulence diffusion function

Bu^ =

scalar Wiener noise for turbulence

BΔ2 =

scalar Wiener noise for load perturbation

FAB =

internal constraint force of linkage AB

I0z =

polar mass moment of inertia of flapping foil about “O” (kgm2)

LC =

impedance of output power circuit (Henries)

M0z =

aeroelastic torque (N m)

M(e.m.) =

electromotive torque (N m)

RC =

resistance of output power circuit (Ohms)

HNL =

nonlinear restoring torque function

QL,Δ2 =

diffusion matrix of load perturbation

Wem =

random state vector

AR =

aspect ratio of the blade-airfoil

α =

flapping angle of the blade-airfoil, about “O”

γ =

nonlinear dimensionless damping (van der Pol)

ΔBC =

soft coupling coefficient, Eq. (4)

δ2 =

random perturbation to parameter d2,m

ϵ =

normalized inertia parameter

ζα =

structural damping ratio of the flapping foil

η3D =

three-dimensional flow effect parameter

ι =

dimensionless induced current, power circuit

κ =

dimensionless cubic torsional stiffness

ΛΥ(2) =

second moment Lyapunov exponent

λRL =

generalized impedance of the power circuit

μae,i =

aeroelastic state (i=1,2)

νae,i =

aeroelastic state (i=1,2)

σd2 =

standard deviation of δ2(τ)

σu^ =

standard deviation of u^

τ =

dimensionless time

τl =

discrete time instant used to find MLE

Υ =

sub-vector of random vector Wem

Φ =

unsteady aeroelastic forcing function

Φ0 =

unsteady aeroelastic forcing function at τ=0

Φ(e.m.c.) =

dimensional electro-mechanical coupling coefficient (N/A)

χBC =

dimensionless stiffness of the flexible shaft

Ψ =

electro-mechanical coupling coefficient

ω =

angular vibration frequency (rad/s)

ωα =

pulsation of the one-DOF flapping foil (rad/s)

Abbreviations and Operators

DOF =

degree-of-freedom

E[.] =

expectation operator

MLE =

moment Lyapunov exponent

T =

transpose operator

References

1.
Priya
,
S.
, and
Inman
,
D. J.
,
2009
,
Energy Harvesting Technologies
,
Springer Science
,
New York
.
2.
Matsumoto
,
M.
,
Mizuno
,
K.
,
Okubo
,
K.
, and
Itô
,
Y.
,
2006
, “
Fundamental Study on the Efficiency of Power Generation System by Use of the Flutter Instability
,” Proc. 2006 ASME Pressure Vessels and Piping Division Conference, Vol. 9,
American Society of Mechanical Engineers
, pp.
277
286
, ASME Paper PVP2006-ICPVT11-93773.
3.
Pigolotti
,
L.
,
Mannini
,
C.
,
Bartoli
,
G.
, and
Thiele
,
K.
,
2017
, “
Critical and Post-critical Behaviour of Two-Degree-of-Freedom Flutter-Based Generators
,”
J. Sound Vib.
,
404
, pp.
116
140
.
4.
Abdelkefi
,
A.
,
Nayfeh
,
A. H.
, and
Hajj
,
M. R.
,
2012
, “
Design of Piezoaeroelastic Energy Harvesters
,”
Nonlinear Dyn.
,
68
(
4
), pp.
519
530
.
5.
Le
,
H.
,
Kwon
,
S.-D.
, and
Law
,
K.
,
2023
, “
Hybrid Energy Harvesting From Wind and Bridge Vibrations
,”
16th Int. Conference on Wind Engineering (ICWE2023)
,
Florence, Italy
,
Aug. 27–31
,
Abstract No. 34
.
6.
Bernitsas
,
M. M.
,
Raghavan
,
K.
,
Ben-Simon
,
Y.
, and
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
,”
ASME J. Offshore Mech. Arct. Eng.
,
130
(
4
), p.
041101
.
7.
Singh
,
K.
,
Michelin
,
S.
, and
de Langre
,
E.
,
2012
, “
Energy Harvesting From Axial Fluid-Elastic Instabilities of a Cylinder
,”
J. Fluids Struct.
,
30
, pp.
159
172
.
8.
Eugeni
,
M.
,
Elahi
,
H.
,
Fune
,
F.
,
Lampani
,
L.
,
Mastroddi
,
F.
,
Romano
,
G. P.
, and
Gaudenzi
,
P.
,
2020
, “
Numerical and Experimental Investigation of Piezoelectric Energy Harvester Based on Flag-Flutter
,”
Aerosp. Sci. Technol.
,
97
, p.
105634
.
9.
Gkoumas
,
K.
,
Petrini
,
F.
, and
Bontempi
,
F.
,
2017
, “
Piezoelectric Vibration Energy Harvesting From Airflow in HVAC (Heating Ventilation and Air Conditioning) Systems
,”
Procedia Eng.
,
199
, pp.
3444
3449
.
10.
Petrini
,
F.
, and
Gkoumas
,
K.
,
2018
, “
Piezoelectric Energy Harvesting From Vortex Shedding and Galloping Induced Vibrations Inside HVAC Ducts
,”
Energy Build.
,
158
, pp.
371
383
.
11.
Young
,
J.
,
Lai
,
J. C. S.
, and
Platzer
,
M. F.
,
2014
, “
A Review of Progress and Challenges in Flapping Foil Power Generation
,”
Prog. Aerosp. Sci.
,
67
, pp.
2
28
.
12.
Caracoglia
,
L.
,
2018
, “
Modeling the Coupled Electro-mechanical Response of a Torsional-Flutter-Based Wind Harvester With a Focus on Energy Efficiency Examination
,”
J. Wind Eng. Ind. Aerodyn.
,
174
, pp.
437
450
.
13.
Ahmadi
,
G.
,
1979
, “
An Oscillatory Wind Energy Converter
,”
Wind Eng.
,
3
, pp.
207
215
.
14.
Roohi
,
R.
,
Hosseini
,
R.
, and
Ahmadi
,
G.
,
2023
, “
Parametric Study of an H-Section Oscillatory Wind Energy Converter
,”
J. Ocean Eng.
,
270
, p.
113652
.
15.
Meimaris
,
A. T.
,
Kougioumtzoglou
,
I. A.
, and
Pantelous
,
A. A.
,
2019
, “
Approximate Analytical Solutions for a Class of Nonlinear Sstochastic Differential Equations
,”
Eur. J. Appl. Math.
,
30
(
5
), pp.
928
944
.
16.
Kwon
,
S.-D.
,
Park
,
J.
, and
Law
,
K.
,
2013
, “
Electromagnetic Energy Harvester With Repulsively Stacked Multilayer Magnets for Low Frequency Vibrations
,”
Smart Mater. Struct.
,
22
(
5
), p.
055007
.
17.
Stephen
,
N. G.
,
2006
, “
On Energy Harvesting From Ambient Vibration
,”
J. Sound Vib.
,
293
(
1–2
), pp.
409
425
.
18.
Caracoglia
,
L.
,
2024
, “
Stochastic Performance of a Torsional-Flutter Harvester in Non-stationary, Turbulent Thunderstorm Outflows
,”
J. Fluids Struct.
,
124
, p.
104050
.
19.
Nayfeh
,
A. H.
, and
Mook
,
D. T.
,
1995
,
Nonlinear Oscillations
,
Wiley–VCH Verlag GmbH & Co. KGaA
,
Weinhein
.
20.
Grigoriu
,
M.
,
2002
,
Stochastic Calculus. Applications in Science and Engineering
,
Birkhäuser
,
Boston, MA
.
21.
Bisplinghoff
,
R. L.
,
Ashley
,
H.
, and
Halfman
,
R. L.
,
1955
,
Aeroelasticity
,
Dover Publications Inc.
,
Mineola, NY
.
22.
Argentina
,
M.
, and
Mahadevan
,
L.
,
2005
, “
Fluid-Flow-Induced Flutter of a Flag
,”
Proc. Natl. Acad. Sci. USA
,
102
(
6
), pp.
1829
1834
.
23.
Jones
,
R. T.
,
1939
, “The Unsteady Lift of a Finite Wing,” Tech. Rep. NACA-TN-682, National Advisory Committee for Aeronautics, Washington, DC.
24.
Itô
,
K.
,
1951
,
On Stochastic Differential Equations
(
Memoirs of the American Mathematical Society
); Vol.
4
,
American Mathematical Society
,
Providence, RI
.
25.
Wong
,
E.
, and
Zakai
,
M.
,
1965
, “
On the Relation Between Ordinary and Stochastic Differential Equations
,”
Int. J. Eng. Sci.
,
3
, pp.
213
229
.
26.
Kloeden
,
P. E.
,
Platen
,
E.
, and
Schurz
,
H.
,
1994
,
Numerical Solution of Stochastic Differential Equations Through Computer Experiments
,
Springer-Verlag
,
Berlin-Heidelberg
.
27.
Xie
,
W.-C.
,
2005
, “
Monte Carlo Simulation of Moment Lyapunov Exponents
,”
ASME J. Appl. Mech.
,
72
(
2
), pp.
269
275
.
28.
Náprstek
,
J.
,
2001
, “
Stability Domains of Wind-Excited Random Nonlinear Systems Through Lyapunov Function
,”
J. Wind Eng. Ind. Aerodyn.
,
89
(
14–15
), pp.
1499
1512
.
29.
Pospíšil
,
S.
,
Náprstek
,
J.
, and
Hračov
,
S.
,
2006
, “
Stability Domains in Flow-Structure Interaction and Influence of Random Noises
,”
J. Wind Eng. Ind. Aerodyn.
,
94
(
11
), pp.
883
893
.
30.
Xie
,
W.-C.
,
2006
,
Dynamic Stability of Structures
,
Cambridge University Press
,
New York
.