Limit-cycle oscillations of a pre-tensed membrane strip

Journal of Fluids and Structures 60 (2016) 1–22

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Limit-cycle oscillations of a pre-tensed membrane strip Ariel Drachinsky n, Daniella E. Raveh Aerospace Engineering Department, Technion – Israel Institute of Technology, Haifa 32000, Israel

a r t i c l e i n f o

abstract

Article history: Received 31 January 2015 Accepted 14 June 2015 Available online 12 November 2015

The paper presents a computational and experimental study of the nonlinear aeroelastic response of a pre-tensed, high aspect-ratio, thin membrane strip. The goal of the study is to derive and validate a computational model that can be used for analysis and design of membrane strips, for the purpose of energy harvesting from flutter at low airspeeds. The mathematical model is based on a pre-tensed-beam model, accounting for bending and torsional stiffening effects due to pretension and large deformations. The aerodynamic model is a potential flow model. The equations of motion are written as a set of nonlinear ordinary differential equations, using Galerkin’s method, and are simulated numerically. The nonlinear aeroelastic model is used to study the oscillation characteristics of the membrane strip in the various stability regions. The effects of the initial pretension and non-linear stiffening on the energy-harvesting potential of the system are studied. The combined effect of the preload on the flutter onset speed, on the flutter frequency and amplitude, and on the loss of orbital stability, indicate that an optimal preload can be determined based on the intended airspeed range for energy harvesting. A series of wind tunnel tests are conducted, in which the flutter onset velocity, and the post-flutter frequencies and amplitudes are measured. Good agreement between the experimental data and computational results validate the computational model. & 2015 Elsevier Ltd. All rights reserved.

1. Introduction In recent years there has been considerable interest in renewable energy resources, and specifically in wind power. One of the suggested devices for energy harvesting in low-speed winds is the Windbelt generator, by Shawn Frayne and Jordan McRae with the Humdinger Wind Energy Company (Ward and Freyne, 2007). The Windbelt device is a taut, high aspectratio membrane that flutters in low-speed winds. Magnets that are attached to the fluttering membrane move in and out of coils, thereby generating electrical power. A somewhat similar device was proposed by Sundararajan et al. (2012), where the mechanical vibrations are transformed into electrical power via piezo-electric devices. These devices motivated the current study of the post-flutter characteristics of a high aspect-ratio membrane (i.e., a membrane strip) in low subsonic flows. Aerodynamic and aeroelastic analyses of 2D membrane wings are well documented in recent literature (e.g. Smith and Shyy, 1995; Smith and Shyy, 1995; Scott et al., 2007; Gordnier, 2009; Rojratsirikul et al., 2009; Tiomkin et al., 2011). These studies typically refer to membranes that are restrained at the leading and trailing edges, and deal with equilibrium of tension, aerodynamics, and inertial forces on the membrane wing. There is also a vast body of literature on aeroelastic stability and limit-cycle oscillations (LCO) of plates and shells, specifically on plates that are clamped on all sides, or on n

Corresponding author. E-mail address:

http://dx.doi.org/10.1016/j.jfluidstructs.2015.06.007 0889-9746/& 2015 Elsevier Ltd. All rights reserved.

2

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

Nomenclature

Qa

a A Ac B Bc Caij

Qe θm V0 vd vf W w Wi

Ci Cmat E fi G H I J Ki k Kaij La Ma μ Nc N0

span length cross section area coil area chord length magnetic field on the center of the coil aerodynamic modal damping, normalized by the airspeed structural modal damping material damping coefficient modulus of elasticity frequency of the ith mode shear modulus membrane thickness area moment of inertia torsional rigidity of the cross section structural stiffness non-dimensional frequency aerodynamic modal stiffness, normalized by the square of the airspeed aerodynamic lift force per unit length aerodynamic pitching moment per unit length magnetic dipole constant number of rings in coil initial membrane pre-tension per unit length

Xi ΔN εx ̄ ζ θ θi ν ρ ρa σ0 ωf ωi Zm

aerodynamic pressure on a membrane element elastic stiffness term angle between magnet and coil axis airspeed divergence velocity flutter velocity plunge displacement of a strip plunge displacement of a membrane element plunge displacement of a strip in the ith mode shape displacement of the ith mode additional tension due to membrane deflection mid-plane strain averaged over the span damping ratio pitch displacement of a strip pitch displacement of a strip in the ith mode shape Poisson's ratio membrane density air density pre-stress in the membrane flutter frequency frequency of the ith mode distance between magnet and coil

cantilevered wings (e.g. Lian et al., 2003; Levin and Karpel, 2011; Patil and Hodges, 2001; Tang and Dowel, 2004). These studies use different models with different levels of complexity, from a simple beam model, to an arbitrary shell structure modeled with finite elements (FE). The aerodynamic models may vary from linear potential strip theory to full computational fluid dynamics (CFD) analysis. A recent study, motivated by NASA’s mold-line link project for noise reduction of transport aircraft, presented numerical analysis and wind-tunnel testing on membranes of low aspect ratio, in various boundary conditions (Kuo et al., 1972). None of these studies consider the aeroelastic stability and characteristics of high aspect-ratio membranes, as in the case of the Windbelt device. A somewhat similar, yet different problem studied in the literature is the problem of panel flutter (Kuo et al., 1972; Dowell, 1970; Wang and Dowell, 2012; Dowell, 1967). These studies focus on the stability and LCO of thin shells of different aspect ratios, and use similar structural models. However, in the problem of panel flutter, only one side of the panel is exposed to airflow, whereas in the case discussed in the current study both sides of the membrane are subjected to airflow. The current study presents a nonlinear LCO analysis of a membrane strip, of high aspect ratio, in low subsonic flow. The structural model is nonlinear, employing Galerkin’s method on a pre-tensed beam model, accounting for the large deformations during oscillation. The aerodynamic model is a linear potential-flow model. The analysis focuses on the variation of the LCO oscillation amplitude and frequency as a function of airspeed, and the various stability regions. A parametric study from an energy harvesting point of view is then presented. Wind tunnel tests serve for validating the mathematical model and its assumptions.

2. Mathematical model Fig. 1 shows the membrane-strip geometry, and the problem setup. The strip dimensions are a
b
h (a- span, b-chord, h-thickness), where a4 4b 4 4h. The membrane strip is pre-tensed at the short ends (or, clamped at one end and tensioned at the other end), while the long ends, which are the membranes’ leading and trailing edges, are free. As customary in dynamic aeroelastic problems, a modal technique is employed, assuming that the displacements of the membrane-strip can be described by superposition of a few characteristic shapes (typically, structural natural modes). A Galerkin approximation is used, with the first four natural modes as shape functions, in order to obtain a set of nonlinear ordinary differential equations (ODEs) describing the membrane’s oscillations due to wind loads. The equations are then integrated numerically using a Runge-Kutta based algorithm.

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

3

N

ρ V b

E,h,ρν ,

a

N Fig. 1. Schematic illustration of the problem.

θ h W

dy

z y

x

dx

b

Fig. 2. Degrees of freedom of the membrane.

2.1. Structural model The structural model is that of a continuous pre-tensed beam. Due to the high aspect ratio of the membrane we assume that there is no chordwise bending, and therefore each strip (of length dx) of the membrane can move in two degrees of freedom (DOF): plunge (W) and pitch (ϴ), as shown in Fig. 2. We verify this assumption in Section 3.1 by comparing modal analysis results computed with this model and with a finite-element model that does include chordwise bending. The mathematical model assumes large displacements, small rotation angles, and small strains. In-plane bending and axial deformations are irrelevant for the analysis and thus are not considered. The equilibrium equation is written for an element dy of the strip (an element dx–dy of the membrane), where the element’s displacements can be written in terms of the strip’s DOFs as: wðx; yÞ ¼ W ðxÞ þθðxÞ Uy:

ð1Þ

Summing the forces on an element dx–dy yields the membrane’s dynamic equation that resembles the string equation:
 ∂ ∂wðx; yÞ ∂2 wðx; yÞ N ðx; yÞ ¼ ρh þ Q a þ Q e; ∂x ∂x ∂t 2

ð2Þ

where N is the tension in the span direction, per unit length, ρ is the membrane material density, and h is the membrane thickness. Qa is an aerodynamic pressure term, and Qe is the elastic term, resulting from bending and torsional rigidity. In Eq. (2) the aerodynamic and elastic terms are considered as external forces. They will be derived later on in this section. The tension in the membrane is provided by: Nðx; yÞ ¼ N 0 þΔN; where N0 is the pre-tension, and ΔN is the additional tension due to the lengthening of the membrane: 0sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1
2  Z Z
Eh a @ ∂w Eh a ∂wðx; yÞ 2 1þ  1Adx  dx; ΔN ¼ Eh U ϵx ðyÞ ¼ a 0 ∂x 2a 0 ∂x

ð3Þ

ð4Þ

where εx is the averaged strain (over the thickness) in the span direction. Substituting Eqs. (1), (3), and (4) into Eq. (2), the equilibrium equation, in terms of the strip’s DOF, becomes: !

2   Z a
∂ W ∂2 θ Eh ∂WðxÞ ∂θðxÞ 2 ∂2 W ∂2 θ þy þ y dx þy 2 þ N0 2 2 2 2a 0 ∂x ∂x ∂x ∂x ∂x ∂x
 ∂2 W ∂2 θ  ρh 2 þ y 2  Q a Q e ¼ 0: ∂t ∂t

ð5Þ

4

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

The elastic terms are formulated for a strip, based on an Euler Bernoulli beam model and Saint Venant’s theorem, and are written in terms of the strip’s cross-section properties EI and GJ: Z b=2 Z b=2 ∂4 W ∂2 θ Q e dy ¼ EI ; Q e y dy ¼ GJ 2 : ð6Þ 4 ∂x ∂x  b=2  b=2 The aerodynamic terms for the lift and pitching moment on a strip, La and Ma respectively, are based on Theodorsen’s unsteady aerodynamic model (detailed in Section 2.2): Z b=2 Z b=2 Q a dy  La ; Q a y dy  Ma : ð7Þ  b=2

 b=2

Summing the forces over the chord and requiring overall force and moment equilibrium leads to the force and moment equations for the strip: ∂2 W ∂2 W ∂4 W þ La  σ 0 A þ EI 2 2 ∂x ∂x4 ∂t   2   ∂2 θ E ∂ W  A UI ww þI p U I θθ  þ 2I U I ¼ 0; p wθ 2a ∂x2 ∂x2

ρA

ρI p

ð8aÞ

 ∂2 θ ∂2 θ þ M a  σ 0 I p þ GJ ∂x2 ∂t 2 " ! # 4 ∂2 W E  b ∂2 θ ¼ 0; 2I p UI wθ  þ I p UI ww þ AU I θθ 2 2a 80 ∂x ∂x2

ð8bÞ

where: Z I ww 

a

0


 ∂W 2 dx; ∂x 3

A ¼ bh;

hb Ip ¼ ; 12

Z I θθ 

a 0




∂θ ∂x

2

Z dx;

I wθ 

a 0




∂θ ∂x


 ∂W dx ∂x

N0 σ0 ¼ : h

ð9Þ

Provided sufficient pre-tension, the elastic bending stiffness is negligible compared to the membrane stiffness, and the bending stiffness terms are therefore omitted. The elastic torsional stiffness, however, cannot be neglected. For solution of the non-linear equation we implement Galerkin's method (Meirovitch, 2003), assuming that the nonlinear solution can be described by a linear combination of shape functions: W ðx; t Þ  X 1 ðt ÞW 1 ðxÞ þ X 2 ðt ÞW 2 ðxÞ þ X 3 ðt ÞW 3 ðxÞ þ X 4 ðt ÞW 4 ðxÞ θðx; t Þ  X 1 ðt Þθ1 ðxÞ þX 2 ðt Þθ2 ðxÞ þX 3 ðt Þθ3 ðxÞ þX 4 ðt Þθ4 ðxÞ;

ð10Þ

where Wi and θi are the plunge and pitch displacements in the ith shape function respectively, and Xi are the generalized displacements. The shape functions are the structural, natural mode shapes presented in Appendix A, obtained from eigensolution of the linearized homogeneous equations: ∂2 W ∂2 W  ρA 2 ¼ 0; ∂x2 ∂t

ð11aÞ

 ∂2 θ ∂2 θ  ρI p 2 ¼ 0: σ 0 I p þ GJ 2 ∂x ∂t

ð11bÞ

σ0 A

Note that in each mode either θi or Wi is zero due to the structural uncoupling of the bending and torsion. Implementing Galerkin's method on the element’s equilibrium equation (Eq. (2)) leads to:   Z a Z b=2 
∂ ∂wðx; yÞ ∂2 wðx; yÞ N ðx; yÞ ρh Q  Q a e U wi dy dx ¼ 0; ∂x ∂t 2 0  b=2 ∂x where wi is the element’s displacement in the ith mode. Substituting Eq. (1) yields: Z a Z b=2 Z a Z b=2 ½r dy UW i ðxÞdx þ ½r  Uy dy U θi ðxÞdx ¼ 0; 0

 b=2

0

 b=2

ð12Þ

ð13Þ

where r

 
 ∂ ∂wðx; yÞ ∂2 wðx; yÞ N ðx; yÞ ρh Qa Qe : 2 ∂x ∂x ∂t

Integrating on dy yields: Z a Z a ðr 1 ÞW i dx þ ðr 2 Þθi dx ¼ 0; 0

0

ð14Þ

ð15Þ

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

5

where r1 and r2 are the left-hand side terms of Eq. (8a) and (8b) respectively. Integrating on dx, four differential equations are obtained: X 1 þK 1 X 1 þξ11 X 31 þ ξ12 X 1 X 22 þ ξ13 X 1 X 23 þ ξ14 X 1 X 24 þξ234 X 2 X 3 X 4 Z a Z a 2 2 þ La U u1 dx þ Ma U u1 dx ¼ 0; Aρa 0 Aρa 0 X 2 þK 2 X 2 þξ21 X 2 X 21 þ ξ22 X 32 þ ξ23 X 2 X 23 þ ξ24 X 4 X 24 þξ134 X 1 X 3 X 4 Z a Z a 2 2 þ La U u2 dx þ M a U u2 dx ¼ 0; I p ρa 0 I p ρa 0 X 3 þK 3 X 3 þξ31 X 3 X 21 þ ξ32 X 3 X 22 þ ξ33 X 33 þ ξ34 X 3 X 24 þξ124 X 1 X 2 X 4 Z a Z 2 2 a þ La U u3 dx þ M a Uu3 dx ¼ 0; Aρa 0 ρa 0

ð17aÞ

ð17bÞ

ð17cÞ

X 4 þK 4 X 4 þξ41 X 4 X 21 þξ42 X 4 X 22 þ ξ43 X 4 X 23 þ ξ44 X 34 þξ123 X 1 X 2 X 3 Z a Z a 2 2 þ La U u4 dx þ M a U u4 dx ¼ 0; I p ρa 0 I p ρa 0 (17d)where the K and ξ terms are defined in Appendix A. Structural, viscous damping can be added in several ways. One approach is to assume a linear relation between the displacement rate and the force:  Z a Z b=2  ∂w U wi dy dx C mat ð18aÞ Ci ¼ ∂t 0  b=2 and adding Eq. (18a) to Eq. (12). Another approach is to assume a linear relation between the strain rate and the damping force, thus between the rate of curvature change and the resultant force  Z a Z b=2  ∂3 w Uwi dy dx Ci ¼ C mat 2 ð18bÞ ∂x ∂t 0  b=2 and adding Eq. (18b) to Eq. (12). Note that Cmat has different dimensions in Eqs. (18a) and (18b). Since the coefficients Cmat are difficult to establish, it is customary to assume a modal damping coefficient, ζi, based on experience or vibration test: C i ¼ 2ζ i ωi :

ð18cÞ

In the following analyses we assume no structural damping, ζi ¼Cmat ¼0. 2.2. Aerodynamic model The aerodynamic lift and moment, per unit length, are based on Theodorsen’s strip theory, assuming potential flow, a flat, rigid, and thin section that oscillates in two degrees of freedom of pitch and plunge, in a constant frequency, with no aerodynamic interaction between sections. The aerodynamic force in each DOF consists of two terms: one due to modal displacement and one due to modal velocity. Z a Z a



_ ; θ; θ_ U W i U dxþ _ ; θ; θ_ U θi U dx La W; W M a W; W 0

0 _ ; θ_ : ð19Þ ¼ FX i ðW; θÞ þF X_ i W Substitution of the oscillatory derivatives found in Wright and Cooper (2007) and integration yields: a a FX 1 ¼ ρa V 20 Lz ðkÞX 1 þρa V 20 Lθ ðkÞX 2 2 2 a a FX 2 ¼ ρa V 20 b M z ðkÞX 1 þ ρa V 20 b M θ ðkÞX 2 2 2 a a FX 3 ¼ ρa V 20 Lz ðkÞX 3 þρa V 20 Lθ ðkÞX 4 2 2 a a FX 4 ¼ ρa V 20 b M z ðkÞX 3 þ ρa V 20 b M θ ðkÞX 4 ; 2 2 ab ab L _ ðkÞX_ 1 þ ρa V 0 Lθ_ ðkÞX_ 2 4 z 4 2 2 ab ab _ _ M z_ ðkÞX 1 þρa V 0 M _ ðkÞX_ 2 F X 2 ¼  ρa V 0 4 4 θ

F X_ 1 ¼  ρa V 0

ð20aÞ

6

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

ab ab L _ ðkÞX_ 3 þ ρa V 0 Lθ_ ðkÞX_ 4 4 z 4 2 2 ab ab Mz_ ðkÞX_ 3 þ ρa V 0 M _ ðkÞX_ 4 ; F X_ 4 ¼ ρa V 0 4 4 θ F X_ 3 ¼ ρa V 0

ð20bÞ

where k is the non-dimensional frequency, defined as: k¼

ωb : 2v

ð21Þ

After normalization of the generalized aerodynamic forces, similarly to the structural terms, and addition of the damping terms, Eq. (22) is obtained. The coefficients are provided in Appendix A. Note that there is no linear coupling between DOFs 1,2 and DOFs 3,4, neither in the structural nor in the aerodynamic terms.   X€ 1 þ ðC 1  vC a11 ÞX_ 1  vC a12 X_ 2 þ K 1  v2 K a11 X 1  v2 K a12 X 2 þ ξ11 X 31 þ ξ12 X 1 X 22 þ ξ13 X 1 X 23 þξ14 X 1 X 24 þ ξ234 X 2 X 3 X 4 ¼ 0;

ð22aÞ

  X€ 2 þ ðC 2  vC a22 ÞX_ 2  vC a21 X_ 1 þ K 2  v2 K a22 X 2  v2 K a21 X 2 þ ξ21 X 2 X 21 þ ξ22 X 32 þ ξ23 X 2 X 23 þξ24 X 4 X 24 þ ξ134 X 1 X 3 X 4 ¼ 0;

ð22bÞ

  X€ 3 þ ðC 3  vC a33 ÞX_ 3  vC a34 X_ 4 þ K 3  v2 K a33 X 3  v2 K a34 X 4 þ ξ31 X 3 X 21 þ ξ32 X 3 X 22 þ ξ33 X 33 þξ34 X 3 X 24 þ ξ124 X 1 X 2 X 4 ¼ 0;

ð22cÞ

  X€ 4 þ ðC 4  vC a44 ÞX_ 3  vC a43 X_ 3 þ K 4  v2 K a44 X 4  v2 K a43 X 3 þ ξ41 X 4 X 21 þ ξ42 X 4 X 22 þ ξ43 X 4 X 23 þξ44 X 34 þ ξ123 X 1 X 2 X 3 ¼ 0:

ð22dÞ

2.3. Solution of the aeroelastic equations Eq. (22) is solved numerically using Matlab’s ODE45 solver, which is an adaptive solver that integrates a system of ODEs using a 4th order Runge-Kutta formula (with a 5th order correction). Since the aerodynamic terms are functions of the oscillation frequency, the solution is followed by spectral analysis. After obtaining the oscillation frequency, the aerodynamic terms are updated, and the problem is simulated again. This procedure is repeated until convergence. 2.4. Power calculation For the sake of comparison and parametric study, an approximate estimation of the energy harvested by the fluttering membrane is derived. The calculation is aimed to be comparative rather than yield accurate values. The energy harvesting system is composed of a magnet that is attached close to the base of the membrane, and two coils attached to a rigid frame, as illustrated in Fig. 3. The calculation is based on the assumption that the displacement of the magnet is small in comparison with the distance between the magnet and the coils. The electromotive force (EMF) is given by:
 ∂ðΦBc Þ ∂ðΦBc Þ ∂ðΦBc Þ Ε ¼ Nc ¼ Nc Zm þ θm ; ð24Þ ∂t ∂Z ∂θ where Nc is the number of rings in the coil, Bc is the electromagnetic field, Zm is the location of the coil, perpendicular to the magnet (and to the membrane), and θm is the angle between the magnet and the coil. For the sake of simplicity, the

Fig. 3. Illustration of the energy harvesting system.

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

electromagnetic field is approximated by that of a dipole: ! ! ! μ μ θ2 μ 3θ2 Bc ¼ 3 ð3 cos θm  1Þ  3 3 1  m  1 ¼ 3 2  m ; 2 2 Zm Zm Zm

7

ð25Þ

where μ is the force of the dipole. Since the magnet is located close to the base of the membrane, where the displacements are small, the modal forces that will result from its presence will be relatively small. Therefore the aeroelastic system is assumed to be unaffected by the electromagnetic system.1 The oscillation of the magnet is therefore defined by the movement of the membrane. The calculation is conducted for a magnet located at 0.1 span, assuming that the coils are aligned with the magnet. Z m ¼ Z 0 þX 1 sin ð0:1π Þ ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    2 θm ¼ θ2 xx þ θ2 yy ¼ ðX 2 sin ð0:1π ÞÞ2 þ  X 1 π=a U cos ð0:1π Þ :

ð26Þ

The energy harvesting system is assumed to be connected to an energy absorber that is modeled by a component with large electrical resistance. The effective resistance resultant from the coil’s resistance and self-inductance is assumed negligible. Therefore, the power is given by; 0 !2 Ε ðt Þ2 9θ2m =2  6 3θ2m ðψNAC μÞ2 @ ð27Þ Ψ ; Z m  3 θm ; Ψ  P ðt Þ  4 R R Zm Zm where Ac is the area of the coil. The power is averaged over a cycle. Assuming harmonic oscillations, and assuming that the displacements are small enough, allowing for realization, the averaging of the power expression over a cycle yields power estimation that is proportional to the square of the velocity amplitude (or to the square of the frequency and the oscillation amplitude).  2 ð28Þ P  p ω UX 1_amp : This approximation does not yield an exact solution for the power values, but can be used to explain the trends presented. 2.5. Strain calculation The analysis is validated with an experiment presented in Section 5. The main measurement in the experiment is the strain in the longitudinal direction, on a pre-selected point on the membrane. Hence, in order to get a better basis for comparison with the experimental data, the strain in the membrane is evaluated numerically. The strain in the membrane is the result of three components: ϵ ¼ ϵ0 þ ϵb þ ϵT :

ð29Þ

The first component is strain due to pre-tension: ϵ0 ¼

N0 : EA

ð30Þ

The second component is the strain due to the bending of the membrane: σ Mh ∂2 w h ¼ εb ¼ ¼ E 2EI ∂x2 2 and according to Eqs. (1) and (10)
2  4 X ∂2 w ∂ Wi ∂2 θi ð x; t Þ  X þ y : i ∂x2 ∂x2 ∂x2 i¼1 Note that though the bending contribution to the stiffness was negligible, the strain due to bending is not. The third component results from the additional tension in the membrane, due to the large deformation (ΔN):  !2 Z a
Z a X 4
ΔN 1 ∂wðx; yÞ 2 X 2i ∂W i ∂θi ¼ þy dx ¼ dx: ϵT ¼ Eh 2a 0 ∂x 2a 0 ∂x ∂x i¼1

ð31Þ

ð32Þ

ð33Þ

Note that while the first component is static, the second is dependent on X(t) and therefore oscillates in the LCO frequency. The third component is dependent on X(t)2 and therefore oscillates in a frequency that is double the LCO frequency, 1 The actual effect of the coupling would lead to additional damping in the system that would slightly increase flutter speeds and reduce the amplitudes. However, though this is expected to decrease the harvested power, the trends and general behavior of the system would not change.

8

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

and not about zero: X  X 0 cos ðωt Þ-X 2 

X0 ð1 þ cos ð2ωt ÞÞ: 2

ð34Þ

3. Numerical test case The membrane in the test case is made of Mylar, with the geometrical parameters and properties detailed in Table 1. 3.1. Linear analysis Modal analysis and linear flutter analysis were conducted for the membrane using three different models as follows: 1. Modal analysis using the pre-tensed beam model derived above. Since flutter analysis is linear, linearization was made about zero deflections. A k-method flutter analysis (Wright and Cooper, 2007) was performed using the resulting fourmode system, with strip-method aerodynamics. 2. Modal analysis in Nastran using shell elements. The resulting 10 modes were used in ZAERO g-method flutter analysis with unsteady panel aerodynamics (Zona 6) (ZONA Technologies, 2011). 3. Modal analysis for a FE beam model, augmented with the geometric stiffness matrix (due to pretension), in a Matlab code. The resulting 10 modes were used in a k-method flutter analysis, with strip theory aerodynamics. The objective of analysis 2 was to verify the assumption that there is no chordwise bending in the first modes, and thus that the pre-tensed beam model is appropriate. In addition, it provided comparison between strip and panel aerodynamic models. The objective of analysis 3 was to support the exclusion of the bending stiffness terms from the pre-tensed beam model, and verify that four modes are sufficient for flutter analysis. A comparison between the results of the three analyses is presented in Table 2, showing good agreement. The first bending frequency computed by model 1 (without the bending stiffness), is slightly lower (by 1.1%), which leads to a 3% reduction in flutter speed. Fig. 4 shows the V-g plot of the pre-tensed-beam model. Flutter occurs at 6.2 m/s for the first bending and first torsion modes, and divergence, of the torsion mode, at 8.6 m/s. Similar plots were obtained with the other two models. The comparison shows good agreement between the three models, both in the structural modal properties and in the flutter and divergence velocities. This validates the beam and strip assumptions, the neglect of the bending term, and the adequacy of the number of modes used in the analysis. 3.2. Nonlinear analysis Eq. (22) was simulated at velocities up to five times the flutter velocity. In order to assess the sensitivity of the problem to initial conditions (IC) the response in each airspeed was simulated twice, with two sets of IC, of large and small displacements, where a large displacement is an order of magnitude larger than the computed LCO amplitude, and a small displacement is an order of magnitude lower. Fig. 5 shows a plot of the maximum displacement of the mid-span leading edge (LE) point at steady oscillations versus the ratio of airspeed to flutter speed (where flutter speed is the linearly computed Table 1 Structural parameters. Property

a

b

h

σ0

E

ρ

ζ

ρa

ν

Value Units

596 mm

25 mm

0.25 mm

3.89 MPa

6980 MPa

1430 kg/m3

0.0 %

1.225 kg/m3

0.39 –

Table 2 Modal and flutter analysis comparison between different models.

1st Bending [Hz] 1st Torsion [Hz] 2nd Bending [Hz] 2nd Torsion [Hz] Vf [m/s] Vd [m/s]

Four modes pre-tensed beam model

10 Modes shell þ panel model

10 Modes FE beam þstrip model

43.9 49.1 87.9 98.2 6.2 8.6

44.2 48.7 88.4 98.3 6.4 8.9

44.4 49.2 88.8 98.5 6.4 8.6

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

0.05

100

0

90 80

mode 1 mode 2 mode 3 mode 4

-0.1 -0.15 -0.2 -0.25

60 50 40 30

-0.3

20

-0.35

10

-0.4

mode 1 mode 2 mode 3 mode 4

70 frequency, Hz

g (added damping)

-0.05

9

0

2

4

6

8

10

12

14

16

18

0

20

0

2

4

6

velocity, m/s

8

10

12

14

16

18

20

velocity, m/s

Fig. 4. Vg plot obtained from the pre-tensed beam model. Left: frequency vs. velocity, right: damping vs. velocity. 0.06

0.05 large IC small IC

(wmax )/a

0.04

0.03

0.02

0.01

0 0.5

1

1.5

2

2.5 v/vf

3

3.5

4

4.5

Fig. 5. LE maximum displacement vs. airspeed.

flutter speed). At high velocities, the oscillation is chaotic and no steady state is achieved. The values plotted at these velocities are the maximum values of the last 100 out of 3000 periods. Fig. 6 shows the average displacement (over time) of the mid-span leading edge point at steady oscillations. Three stability thresholds are seen in Figs. 5 and 6, dividing the membrane oscillations into four regions, as a function of airspeed, as described in Table 3. Fig. 7 shows the modal displacements as a function of airspeed (where the modal bending deformations are normalized by span length, and the torsion deformations by span-over-half-chord). It is seen that the third and fourth modes remain stable (oscillations decay). This is unlike in the linear stability analysis where there is also flutter of the third and fourth modes. This is due to the stiffening, resulting from large deformations of the first and second modes. Fig. 8 shows the membrane’s mid-span leading edge oscillations in each of the four regions. While in region 2 the oscillations are about zero-displacement, in region 3 there is a static offset displacement. It is due to this static offset that the maximum displacement in Fig. 4 appears to depend on the IC. If this static offset is subtracted from the solution, both branches coincide. Spectral analysis was conducted at each speed, using the FFT algorithm. Fig. 9 shows the variation of the dominant frequency as a function of airspeed. It is seen that the frequency increases bi-linearly with the airspeed, independent of the IC. The increase in frequency results from aerodynamic stiffness, and also from the nonlinear elastic term, which is governed by the cube of the displacement magnitude. The change in frequency growth rate seen at 2.2 times flutter speed corresponds to the static offset in the membrane oscillation. Sample results of the spectral analysis from each of the stability regions are shown in Fig. 10. Though the oscillation is non-harmonic, a dominant frequency appears in all cases, for both modes. In regions one and two, the oscillation is very close to harmonic, and only one frequency is observed. Note that in region one, the amplitude is very small. This is due to

10

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

4

x 10

-3

large IC small IC

3 2

waverage /a

1 0 -1 -2 -3 -4 0.5

1

1.5

2

2.5 v/vf

3

3.5

4

4.5

Fig. 6. LE average displacement vs. airspeed. Table 3 Different stability regions. Region v/vf 1 2 3 4

0–1 Stable. All oscillations decay. 1–2.2 LCO 1st and 2nd modes. 2.2–3.7 LCO 1st and 2nd modes with a static offset (the oscillation is about non-zero deflection). The onset of this phenomenon does not occur at the linearly computed divergence speed. 4 3.7 Gradual loss of periodicity, turning into chaotic oscillations.

0.045 mode 1 large IC mode 2 large IC mode 3 large IC mode 4 large IC mode 1 small IC mode 2 small IC mode 3 small IC mode 4 small IC

0.04

Xmax (normalized)

0.035 0.03 0.025 0.02 0.015 0.01 0.005 0 0.5

1

1.5

2

2.5 v/vf

3

3.5

4

4.5

Fig. 7. Modal amplitude vs. airspeed.

the decay of the oscillations in this region. Region two is characterized by a classical LCO. In region three, a large peak at zero frequency corresponds to the static offset of the vibration. In addition, a content of higher harmonics can be seen, especially in the torsion (second) mode. In region four, a dominant frequency, accompanied by a broad spectrum of low magnitudes is observed. We recall that the aerodynamic model assumes vibration in a single frequency. The oscillations in regions three and four do not adhere to this assumption, thus the results there are phenomenological and not necessarily accurate. Fig. 11 shows a time history of each of the strain components discussed in Section 2.5, and the total strain, at a sample airspeed. It shows that although the bending term is negligible in the stiffness calculation, it affects the measured strain.

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

2

x 10

-4

region-1, u/uf=0.9

6

x 10

-3

11

region-2, u/uf=1.37

X LE (normalized)

X LE (normalized)

4 1

0

-1

2 0 -2 -4

-2 2.5

2.6

2.7

-6 2.5

2.8

2.6

t [sec] region-3, u/uf=3.46

0.03

2.7

2.8

t [sec] region-4, u/uf=4.09

0.05

X LE (normalized)

XLE (normalized)

0.02 0.01 0 -0.01

0

-0.02 -0.03 2.5

2.6

2.7

-0.05 2.5

2.8

2.6

t [sec]

2.7

2.8

t [sec]

Fig. 8. Sample time history plots at the different stability regions.

95 large IC small IC

oscillation frequency/flutter frequency

90 85 80 75 70 65 60 55 50 45 0.5

1

1.5

2

2.5 v/vf

3

3.5

4

4.5

Fig. 9. Oscillation frequency vs. airspeed.

3.3. Sensitivity to number of modes Sensitivity analysis was conducted to study the number of modes required for the analysis. The analysis was conducted using two, four and six modes. The results were similar and are presented in Fig. 12. At speeds lower than the speed of the second bifurcation, all higher modes decayed with time. In the chaotic region, following the second bifurcation, the higher modes do not necessarily decay. However, since the model is not valid in this region, the actual modal response is not of significance, only the bifurcation point is. Moreover, in the experiment the oscillations were observed with a stroboscope. The results show a clear bending-torsion-flutter oscillation. If higher modes were present, their contribution was small.

12

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

region-1, u/uf=0.99

0

X1 X2

-2

10

-4

10

-6

10

region-2, u/uf=1.37

0

10

amplitude

amplitude

10

-2

10

-4

10

-6

0

2

4

10

6

0

2

4

w/wf region-3, u/uf=3.46

0

-2

10

-4

10

-6

10

region-4, u/uf=4.09

0

10

amplitude

amplitude

10

6

w/wf

-2

10

-4

10

-6

0

2

4

6

10

0

2

4

w/wf

6

w/wf

Fig. 10. Frequency content at the different stability regions.

9

x 10

-4

u/uf=1.55

8 7

mid chord strain [-]

6 5

preload term tenssion term bending term total

4 3 2 1 0 -1 2.5

2.505

2.51

2.515

2.52

2.525 2.53 t [sec]

2.535

2.54

2.545

2.55

Fig. 11. Strain time history.

4. Parametric study The effects of two parameters were examined in the parametric study:

 The preload.  The non-linearity-parameter λ. 4.1. Preload A plot presenting the flutter/divergence onset, and the second bifurcation point, is presented in Fig. 13. The second bifurcation point is the upper limit of the validity of the analysis model, since above this threshold the oscillation is in

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

13

0.045 mode 1 four modes mode 2 four modes mode 3 four modes mode 4 four modes mode 1 two modes mode 2 two modes mode 1 six modes mode 2 six modes mode 3 six modes mode 4 six modes mode 5 six modes mode 6 six modes

0.04 0.035

(wmax )/a

0.03 0.025 0.02 0.015 0.01 0.005 0 0.5

1

1.5

2

2.5 v/vf

3

3.5

4

4.5

Fig. 12. Modal results for sensitivity analysis to mode number.

wind velocity [m/s]

15

10

5

calculated flutter velocity calculated divergence velocity non linear second bifurcation 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

preload [kg] Fig. 13. Flutter, divergence, and second bifurcation onset velocities as a function of preload.

multiple frequencies. In the experiment, above these speeds the oscillations lost orbital stability. The large amplitudes nonperiodic oscillations make this threshold a limit on energy harvesting. Fig. 13 shows that as the preload is increased, the “working region”, namely the range of wind velocities that result in periodic oscillations decreases. Higher preload means smaller working region for energy harvesting. A plot of the average power at different wind velocities, for different values of preload, is presented in Fig. 14. The plot shows that the input energy increases with airspeed for all preload values. Lower preload values tend to yield more power. This can be explained by the dependency of the power on the oscillation amplitude and the frequency. The increase of frequency with the preload is much smaller than the decrease in the amplitude (not shown). At a preload of 0.5–1 kg, the differences become small, and optimum lower bound for the preload is reached. An optimal pretension cannot be determined for the whole airspeed range. A specific pretension would have to be determined based on the target airspeed range.

4.2. Divergence before flutter Fig. 13 shows that for high values of preload linear flutter analysis predicts divergence before flutter. In the non-linear analysis, as in the experiment, there is no static divergence, and the membrane always oscillates upon stability loss. This is due to the nonlinear stiffening.

14

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

The effect of the nonlinear stiffening on divergence can be explained using a simplified, static 2D model shown in Fig. 15. A plate is mounted on a nonlinear spring having linear and cubic components. Moment equilibrium on the plate yields: ρV 0 2 SC Lα eα=2 ¼ K 0 α þ ξα3 :

ð35Þ

The non-trivial solution of the equation is: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r

α¼ 7 ρV 0 2 SeC Lα  2K 0 =2ξ:

ð36Þ

As long as the numerator is negative, only the trivial solution exists. As the numerator becomes positive, a new stable solution is created. This is the static divergence. Unlike in the linear case, the deflection in this “divergence” is finite. As a result of the static deflection the stiffness is increased, thus increasing the frequencies of the structure, and changing the flutter velocity. In the membrane case, as a result of the stiffness change, the bending and torsion frequencies get closer, thus reducing the flutter frequency. In order to verify this assumption the following analysis was conducted for the reduced two mode system: A static nonlinear analysis was conducted to obtain the post-divergence fixed-point equilibrium point, by solving Eq. (37):   ð37aÞ K 1 v2 K a11 X 1  v2 K a12 X 2 þξ11 X 31 þ ξ12 X 1 X 22 ¼ 0;   K 2 v2 K a22 X 2  v2 K a21 X 2 þξ21 X 2 X 21 þ ξ22 X 32 ¼ 0: The structural stiffness was updated based on linearization about the new equilibrium point: " # 2ξ12 X 1 X 2 K 1 þ 3ξ11 X 21 þ ξ12 X 22 : GK post_divergence ¼ 2ξ21 X 1 X 2 K 2 þ 3ξ22 X 22 þ ξ21 X 21

ð37bÞ

ð38Þ

Finally, a new flutter solution was computed with the updated stiffness matrix. An example for this procedure is given for a case of 4.3 kg pretension. The divergence velocity is 10.8 m/s. The linear flutter velocity is 15 m/s. The static equilibrium solutions of the two-mode system for different velocities is given in Fig. 16. The corresponding updated stiffness matrix is presented in Fig. 17. The variation of the flutter velocity and flutter frequency calculated with the post divergence stiffness matrix is presented in Fig. 18

Fig. 14. Average power as a function of airspeed, for different values of preload.

Fig. 15. : 2D strip on a nonlinear spring model.

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

15

Fig. 18 shows the flutter speed computed with the updated, post-divergence stiffness matrix, presenting a steep drop in the flutter speed, and an increase in the frequency. The flutter velocity drops below the divergence velocity almost immediately, then gradually rises as the stiffness of the structure is increased. Fig. 19 shows the time history of the membrane motion, showing the onset of divergence and flutter with 4.3 kg preload, for airspeeds above the divergence and below flutter speed. The displacement starts as non-oscillating divergence (0 oto0.1), and then turns into LCO. A few limitations on the discussion above are noted: The oscillations shown in Fig. 19 start at zero frequency and develop into LCO of a certain frequency. In between it switches through a range of frequencies. The aerodynamic model used to compute this response assumes a single frequency – in this case, the LCO frequency. Therefore, the results shown in Fig. 19 are not accurate, but rather serve to demonstrate the mechanism of the phenomenon 0.015 X1static X2static *b/2

0.01

X/a

0.005

0

-0.005

-0.01

-0.015 0.95

1

1.05

1.1

1.15

1.2

1.25

1.3

v/vd

Fig. 16. Static post-divergence solution equilibrium point.

x 10

5

2.5

800

600 k12

k11

2 400

1.5 200

1 0.95

1

1.05

1.1

1.15 v/vd

1.2

1.25

0 0.95

1.3

7

1

1.05

1.1

1

1.05

1.1

1.15 v/vd

1.2

1.25

1.3

1.15

1.2

1.25

1.3

5

x 10

x 10

2

3

1.5 k22

k21

2.5 1

2 0.5

0 0.95

1

1.05

1.1

1.15 v/vd

1.2

1.25

1.3

1.5 0.95

v/vd

Fig. 17. Linearized stiffness matrix about post divergence equilibrium point.

16

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

1.6

100 90

1.4

80

ff(post divergence ) [Hz]

vf(post divergence )/vd

1.2

1

0.8

0.6

70 60 50 40 30 20

0.4 10 0.2 0.9

0.95

1

1.05

1.1 v/vd

1.15

1.2

1.25

1.3

0 0.9

0.95

1

1.05

1.1 v/vd

1.15

1.2

1.25

1.3

Fig. 18. post-divergence flutter velocity (right) and flutter frequency (left).

Fig. 19. Time history of the oscillation onset preload¼4.3 kg, velocity¼ 12.5 m/s.

of LCO at divergence speed. For the same reason, it is very difficult to accurately predict the onset of flutter in the “divergence before flutter” preload regions with the aerodynamic model used in this study. 4.3. Non-linearity parameter λ The nonlinearity parameter λ is defined in Eq. (A2). It multiplies all of the nonlinear terms in the membrane’s equation of motion (A2). The λ parameter represents the magnitude of the nonlinear stiffness terms, and is linearly dependent on the modulus of elasticity. Figs. 20 and 21 show the power and oscillation-amplitude at different airspeeds and for different λ values. The oscillation frequency is independent of λ. The relation between the oscillation amplitude and λ and between the power and λ is given in Eq. (34). It is noted that although λ is dependent on the modulus of elasticity, changing the modulus of elasticity has a wider effect on the problem due to its effect on the shear modulus, and thus on the torsion frequencies and flutter velocity. 1 Amplitude p pffiffiffi; λ

1 Pp : λ

ð39Þ

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

17

300

λ=1.62 λ=2.7 λ=3.77, nominal λ=4.86 λ=5.94 λ=7.56 λ=10.8

250

power [watt]

200

150

100

50

0

5

6

7

8 9 velocity [m/sec]

10

11

12

Fig. 20. Variation of average power with airspeed for different values of λ, preload of 2.5 kg. 0.01

λ=1.62 λ=2.7 λ=3.77, nominal λ=4.86 λ=5.94 λ=7.56 λ=10.8

0.009

amplitude (normalized)

0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0

5

6

7

8 9 velocity [m/sec]

10

11

12

Fig. 21. Variation of oscillation amplitude with airspeed for different values of λ, preload of 2.5 kg.

5. Experimental study 5.1. Setup To verify the computational model a series of wind tunnel tests were conducted, in which a 0.025 by 0.6 m Mylar membrane was clamped to the tunnel ceiling and pre-loaded by various weights, ranging from 0.5 to 4.5 kg. Fig. 22 shows a schematic sketch and a photo of the test setup. By gradually increasing the flow speed divergence and flutter speeds were measured for each loading case. At speeds above flutter speed the frequency and amplitude of the membrane vibrations were recorded. The measurement setup consists of two strain gauges, an accelerometer, a stroboscope, and a simple camera. The airspeed is measured by the internal wind speed gauges that are part of the wind tunnel. All of the measurement devices except for the stroboscope are connected to the wind-tunnel data acquisition system, and are thus synchronized. The two unidirectional longitudinal strain gauges are mounted on two opposite sides of the membrane middle axis. The gauges are located  50 mm above the bottom clamping device (the exact distance depends on the pretension). These gauges are used for two purposes: Measure the oscillation’s amplitude and frequency, and provided online measurement of the pre-tension. The latter is highly important since the preload is somewhat reduced during the experiment- First due to the clamping, and then during the oscillations. The accelerometer is attached to the lower clamping device, close to the base of the membrane. Though the clamping device is relatively rigid in comparison with the membrane, the membrane-induced oscillations were large enough to yield sufficient readings for spectral analysis. In addition to the strain gauges and the

18

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

accelerometer, in certain parts of the experiment a stroboscope was used to visualize the oscillation’s mode shape. This was done by setting the scope to a frequency that is close to the flutter frequency. Since the stroboscope operates in a single frequency, when the oscillation is harmonic (or relatively close to harmonic), a clean mode shape is seen. When the oscillation is non-harmonic the picture obtained by the stroboscope is smeared, giving a good indication on the loss of orbital stability. The experiments took place in the Aerospace Engineering Department’s subsonic wind tunnel. The wind tunnel has an open stream circulation. It allows work at wind speeds in the range of 0–90 m/s. The inlet has an area ratio of 1:25. Eight nets are installed in it in order to reduce turbulence and the effect of outside wind gusts. The experiment cell is built from four frames that, together, create a square section with height and width of one meter and length of three meters. The wind tunnel enables Reynolds numbers of up to six million, and low turbulence level of less than half percent. The stagnation pressure is atmospheric pressure, and the air temperature is the same as the outside temperature.

Wind Top clamping

tunnel

device Tightening Membrane

tablets Long plates

Drilled String

stiffening

attachment

tablets

element

String

Weight

Fig. 22. Test setup, left: schematic sketch, right: photo.

5 4.5

average tension [kg]

4 3.5 3 2.5 2 1.5 1

free, before test clamped, before test clamped, after test

0.5 0

0

1

2

3

4

5

weight [kg] Fig. 23. Measured vs. actual preload at different stages of each experiment.

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

19

5.2. Test procedure and results The test was performed with various preloading weights, ranging from 0.5 to 4.5 kg. The test procedure for each loading case was as follows: First the bottom side of the membrane was unrestrained, allowing it to lengthen freely. Then a weight was suspended on the membrane bottom end, and a measurement of the strain was taken. The bottom clamping device was tightened on the membrane and another strain measurement was taken. The airspeed in the tunnel was gradually raised until flutter onset. After the flutter velocity was established, the speed was slightly decreased and an airspeed sweep through the flutter velocity was recorded. The LCOs at various airflow velocities above flutter were then recorded, in increments of  0.5 m/s. The length of each recording time window was about three seconds. The airspeed was then raised until loss of orbital stability. Finally, the airspeed was gradually lowered to zero while recording the oscillations. Fig. 23 shows the preload as measured on the free (unclamped) membrane and on the clamped membrane, before and after each test, showing that the preload was reduced due to the clamping and also during the oscillations. For comparison with results from numerical analysis, the values of the clamped load before test were used. Fig. 24 shows a comparison of the flutter (and divergence) velocities and frequencies for each loading case, as computed numerically and measured in the wind tunnel. The flutter frequencies are in good agreement with the numerical results. The flutter speed is in good agreement up to preload of  3.5, where the analysis yields divergence before flutter. In the

15 60

frequency [Hz]

10

5

0

0.5

1

1.5

2 2.5 3 measured preload [kg]

3.5

4

40

30

20

calculated flutter velocity calculated divergence velocity experimental flutter velocity analytical second bifurcation experimantal second bifurcation

calculated flutter frequency calculated first frequency calculated second frequency experimental flutter frequency

10

4.5

0

0.5

1

1.5

2 2.5 3 measured preload [kg]

3.5

4

Fig. 24. Flutter velocities (left) and frequencies (right) as a function of preload – Comparison of experimental and numerical results.

65 60 55

frequency [Hz]

wind velocity [m/s]

50

50 45

N0=1.50, analysis N0=2.25, analysis N0=3.00, analysis N0=3.50, analysis N0=4.30, analysis N0=1.48, experiment N0=2.23, experiment N0=2.98, experiment N0=3.79, experiment N0=4.18, experiment

40 35 30 25

4

6

8

10

12

14

velocity [m/sec] Fig. 25. Variation of LCO frequency with airspeed for different preload values: analysis vs. experiment.

4.5

20

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

7

x 10

-4

6

strain [-]

5

4

3

2

experiment calculated

1

0 0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

time [sec] Fig. 26. Time history comparison of strains, pretension¼2.25 kg, air velocity¼ 8.1 m/s.

experiment there is no divergence, but the LCO onset velocities in this airspeed range correspond to the computed divergence velocities, as discussed in Section 4.2. The second bifurcation velocities, measured are slightly lower than the computed but the trend is similar. Fig. 25 shows the change in LCO frequency with increased airspeed, for different values of preload, as obtained in the experiment and in the analysis. The resolution of the numerical results is 0.61 Hz. The numerical results are in good agreement with the experiment. The loss of orbital stability in the experiment, in the low pretension values, occurred in a lower speed than that predicted in the analysis (not shown). Fig. 26 shows a comparison of sample time history plot of the strains from the analysis and the experiment. The plot is taken from the results with 2.25 kg preload at airspeed of 8.1 m/s. A slight static offset appears between the two sets of results. This offset probably results from the preload reduction during the oscillation that was presented in Fig. 23. Otherwise, the results are in good agreement. Both plots show an oscillation in two frequencies, the higher of which is dominant, as discussed in Section 2.5 and shown in Section 3.2, Fig. 10. The plot also shows good correspondence between the calculated amplitude and the amplitude recorded in the experiment.

6. Summary and conclusions The paper presented a numerical and experimental study of the phenomenon of LCO in a pre-tensed membrane strip, for the purpose of energy harvesting. The mathematical model used in the numerical analysis was of a pre-tensed beam, with nonlinear, coupled bending and torsional stiffening due to large deformations. The aerodynamic model was a potential strip theory model. The experimental study was based on a series of wind tunnel tests, in which the flutter onset velocity was detected, and the post-flutter LCO oscillations were recorded via strain gauge readings. The numerical results were in good agreement with the experimental results, in terms of flutter velocity and oscillation frequencies, for most preload values. In the region of high preload, linear stability analysis predicts divergence before flutter. In the experiment, and in the nonlinear analysis, LCO occurred at the (linearly predicted) divergence speed. This was explained to be a result of the nonlinear stiffening, due to large deformations. From an energy harvesting point of view, the study showed that smaller values of preload generally yield more power due to higher oscillation amplitudes. However a lower bound is reached where the effect of increase in oscillation amplitude and decrease in frequency are canceled out. It was also shown that for each preload value there is a working region of airspeeds suitable for energy harvesting. This region is bounded from below by the flutter velocity and from above by the second stability threshold, after which the oscillations lose their orbital stability and become non-periodic. This working region becomes smaller with higher preload. Further research is recommended in which the coupling between the electromagnetic and aeroelastic systems is taken into account. The current study does not accurately predict the LCO beyond the loss of orbital stability and the onset of LCO at the divergence velocities. This could be improved by use of time domain, instead of frequency domain, aerodynamics.

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

21

Appendix A. Mode shapes and coefficients of the final aeroelastic equation A.1. Mode shapes and frequencies:

W 1 ðxÞ ¼ sin

πx

; a

θ1 ðxÞ ¼ 0;

ω1 ¼

rffiffiffiffiffi σ0 ρ sffiffiffiffiffiffiffiffiffiffiffiffiffi GJ π σ 0 þ Ip π a

πx

; ω2 ¼ θ2 ðxÞ ¼ sin a a ρ
 rffiffiffiffiffi 2πx 2π σ 0 ; θ3 ðxÞ ¼ 0; ω3 ¼ W 3 ðxÞ ¼ sin a a ρ sffiffiffiffiffiffiffiffiffiffiffiffiffi
 GJ 2πx 2π σ 0 þ Ip ; ω4 ¼ W 4 ðxÞ ¼ 0; θ4 ðxÞ ¼ sin a a ρ W 2 ðxÞ ¼ 0;

ðA1Þ

A.2. Coefficients of the final aeroelastic equation: Lz ðkÞ K a1 1 ¼ K a33 ¼  ρa ρA ;

K a21 ¼ K a43 ¼

 bρa M z ðkÞ ; 2ρI p

C a11 ¼ C a33 ¼  bρa Lz ðkÞ ; 2ρA

a Lα ðkÞ K a1 2 ¼ K a34 ¼ bρ2ρA

2

K a22 ¼ K a44 ¼

b ρa M α ðkÞ 4ρI p

2

C a1 2 ¼ C a34 ¼ b ρa Lα ðkÞ

2

4ρA

3

C a21 ¼ C a43 ¼  b ρa M z ðkÞ C a22 ¼ C a44 ¼ b ρa M α ðkÞ 4ρIp

K 1 ¼ ω21 ;

K 2 ¼ ω22 ; 2

K 3 ¼ 4ω21 ;

8ρIp

K 4 ¼ 4ω22 2

2

b λ b b λ ; ξ13 ¼ λ; ξ14 ¼ λ; ξ234 ¼ ; 16 6 12 2 2 3b λ 3b λ ξ22 ¼ ; ξ23 ¼ λ; ξ24 ¼ ; ξ134 ¼ 2λ; 80 20 2 2 b b 2 ξ31 ¼ λ; ξ32 ¼ λ; ξ33 ¼ 4λ; ξ34 ¼ b λ; ξ124 ¼ λ; 12 6 2 2 3b λ 3b λ ξ41 ¼ 1λ; ξ42 ¼ ; ξ43 ¼ 12λ; ξ44 ¼ ; ξ123 ¼ 2λ; 20 5 C 1 ¼ 2ζ1 ω1 ; C 2 ¼ 2ζ 2 ω2 ; C 3 ¼ 2ζ 3 ω3 ; C 4 ¼ 2ζ 4 ω4 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffi Eπ 4 π σ0 π σ 0 þGJ=I p ; ω2  λ  4 ; ω1  a ρ a ρ ρa 1 ξ11 ¼ λ; 4 3λ ξ21 ¼ ; 4

ξ12 ¼

ðA2Þ

References Dowell, E.H., 1970. Panel flutter: a review of the aeroelastic stability of plates and shells. AIAA Journal 8 (3), 385–399. Dowell, E.H., 1967. Nonlinear oscillations of a fluttering plate. II. AIAA Journal 5 (10), 1856–1862. Gordnier, R.E., 2009. High-fidelity computational simulation of a membrane wing airfoil. Journal of Fluids and Structures 25, 897–917. Kuo, C.C., Morino, L., Dugundji, J., 1972. Perturbation and harmonic balance methods for nonlinear panel flutter. AIAA Journal 10 (11), 1479–1484. Lian, Y., Shyy, W., Viieru, D., Zhang, B., 2003. Membrane wing aerodynamics for micro air vehicles. Progress in Aerospace Sciences 39, 425–465. Levin, D., Karpel, M., 2011. Limit cycle oscillations of plate-type fins using increased-order models. Journal of Aircraft 48 (2), 715–719. Meirovitch, L., 2003. Methods of Analytical Dynamics. McGraw-Hill Book Company, New York (Sixteenth Printing, 1993), Dover Publications, Mineola, NY. Patil, M.J., Hodges, D.H., 2001. Limit cycle oscillations in high-aspect-ratio wings. Journal of Fluids and Structures 15, 107–132. Rojratsirikul, P., Wang, Z., Gursul, I., 2009. Unsteady fluid–structure interactions of membrane airfoils at low reynolds numbers. Experiments in Fluids 46, 859–872. Sundararajan, V., Romero, E., Bonilla, N., Martinez, C., 2012. Energy Harvesting from Fluttering Membranes. In: Proceedings of the 10th Latin American and Caribbean Conference for Engineering and Technology, Panama City, Panama, pp. 1–6. Smith, R.W., Shyy, W., 1995. Computational model of flexible membrane wings in steady laminar flow. AIAA Journal 33 (10), 1769–1777. Smith, R.W., Shyy, W., 1995. Computation of unsteady laminar flow over a flexible two-dimensional membrane wing. Physics of Fluids 7 (9), 2175–2184. Scott, R.C., Bartels, R.E., Kandil, O.A., 2007. An aeroelastic analysis of a thin flexible membrane. In: Proceedings of the 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA, Honolulu, Hawaii, pp. 7459–7475, AIAA-2007-2316. Tiomkin, S., Raveh, D.E., Arieli, R., 2011. Parametric study of a two dimensional membrane wing in viscous laminar flow. In: Proceedings of the 29th AIAA Applied Aerodynamics Conference, AIAA, Honolulu, Hawaii, AIAA-2011-3023. Tang, D.M., Dowel, E.H., 2004. Effects of geometric structural nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings. Journal of Fluids and Structures 19, 291–306. Ward, L., Freyne, S., 2007. Windbelt, Cheap Generator Alternative, Set to Power Third World, June. 〈http://www.popularmechanics.com/science/energy/ solarwind/4224763?series ¼ 37an〉.

22

A. Drachinsky, D.E. Raveh / Journal of Fluids and Structures 60 (2016) 1–22

Wang, I., Dowell, E.H., 2012. Flutter of rectangular plates in three dimensional incompressible flow with various boundary conditions: theory and experiment. In: Proceedings of the ASME International Design Engineering Technical Conferences & Computers and Engineering Conference, ASME, Chicago, IL, pp. 1–13. Wright, J.R., Cooper, J.E., 2007. Introduction to Aircraft Aeroelasticity and Loads. American Institute of Aeronautics and Astronautics160. ZONA Technologies, 2011. ZAERO Theoretical Manual, Version 8.5, ZONA Technologies. pp.3-16–3-31.