Exploiting Stochastic Resonance Principles to Influence the Efficiency of a Torsional-Flutter Energy Harvester in Turbulent Winds¹

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.

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.

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].

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].

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].

If the mass $mC≈0$ 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)≈(1−a)bα(t)$⁠.

For a rigid shaft, it is noted that $χBC→+∞$ and $ΔBC→1$⁠.

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].

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τ)=[κα3−2ζαγα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ζα$⁠.

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.

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).

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.

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).

$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×7$$QL,Δ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$⁠.

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.

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.

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$⁠.

Numerical simulations are executed in the range $0<Umax≤20m/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.

Figure 3 illustrates the reference simulation result for the linear model (⁠$κ=1E−7$⁠,$γ=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$⁠.

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).

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.

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.

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)).

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.

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.

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.

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 and models that support the findings of this study are available from the corresponding author upon reasonable request.

$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$⁠)

DOF =

degree-of-freedom

$E[.]$ =

expectation operator

MLE =

moment Lyapunov exponent

$T$ =

transpose operator