Aeroelastic analysis and nonlinear dynamics of an elastically mounted wing

Journal of Sound and Vibration 331 (2012) 5774–5787

Contents lists available at SciVerse ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Aeroelastic analysis and nonlinear dynamics of an elastically mounted wing M. Ghommem a, A. Abdelkefi a, A.O. Nuhait b, M.R. Hajj a,n a b

Department of Engineering Science and Mechanics, MC 0219, Virginia Polytechnic Institute and State University, Blacksburg, Virginia 24061, USA Department of Mechanical Engineering, King Saud University, Riyadh 11421, Saudi Arabia

a r t i c l e in f o

abstract

Article history: Received 13 February 2012 Accepted 19 July 2012 Handling Editor: L.N. Virgin Available online 18 August 2012

The response of an elastically mounted wing that is free to plunge and pitch, supported by nonlinear translational and torsional springs, and interacting with an incoming stream is analyzed. A tightly coupled model of the wing flow interaction is developed. A three-dimensional code based on the unsteady vortex lattice method is used for the prediction of the unsteady aerodynamic loads. The response of the wing shows a sequence of static and dynamic bifurcations and chaotic motions when increasing the flow speed. Pairs of stable solutions are observed over the different response regimes. The effects of the gust and structural nonlinearity on the wing’s response are also investigated. The results show that gust may lead to jumps between the pairs of solutions for static and dynamic equilibrium responses without impacting the boundaries of the different response regimes. As for the effect of the structural nonlinearity, increasing the nonlinear coefficient of the stiffness of the torsional spring yields lower static deflections and amplitudes of the limit cycle oscillations. & 2012 Elsevier Ltd. All rights reserved.

1. Introduction The aeroelastic response of structures and aircraft is of significant interest for many reasons. Advances in aircraft performance will only be met by using light-weight materials that may be very flexible. The generation of more power from wind turbines will require the use of larger blades which can only be realized through using lighter materials. These requirements have their own penalties in terms of aeroelastic constraints. Such constraints also play a role in the design and dynamic stabilization of long spanned bridges. In other applications, the aeroelastic response of wings has been proposed as means for energy harvesting to power micro sensors or actuators [1–3]. In all of these applications, competing nonlinear mechanisms (geometric, inertia, free-play, damping, aerodynamic,y) within aeroelastic systems may lead to different nonlinear dynamic response characteristics such as limit cycle oscillations (LCO), bifurcation, internal resonances, and chaotic motions [4–8]. More importantly, the transition from one state to another is dependent on the system parameters and operational conditions. In general, these responses are not easily characterized. Yet, the ability to predict them is very important for implementation of control strategies, optimization of the system’s performance or quantification of uncertainties associated with variations in its parameters and operational conditions. Numerical simulations of the response of aeroelastic systems can be carried out with different levels of fidelity. Although high fidelity simulations are implementable, their extensive computational requirements and cost limit the

n

Corresponding author. Tel.: þ 1 5402314190; fax: þ1 5402314574. E-mail address: [email protected] (M.R. Hajj).

0022-460X/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jsv.2012.07.040

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

5775

capability to perform detailed analyses that are based on multiple simulations as would be required for a full system characterization. On the other hand, simulations that make use of simplified models of the aerodynamic flow and structural responses have been shown to be very useful in the analysis of aeroelastic systems. Of particular interest is the combination of the unsteady vortex lattice method (UVLM) to determine the aerodynamic loads with different fidelity levels of structural models. Nuhait and Mook [9] used the two-dimensional unsteady vortex lattice method (UVLM) to compute the unsteady aerodynamic loads on a wing and its oscillations when placed in a uniform flow speed. Furthermore, they used the method of images to simulate ground effects and found that the wing may be significantly destabilized when placed in the proximity of the ground. Wang et al. [10] considered a nonlinear aeroelastic model that combined a nonlinear finite element implementation of a geometrically nonlinear beam model and a three-dimensional version of UVLM. The governing equations of the coupled dynamical system were integrated simultaneously to predict accurately the wing response. The aeroelastic framework was implemented for a rigid wing model, Goland wing, and a Hale type of aircraft. Their results were compared and validated against reported numerical computations. In this work, the aeroelastic response of an elastically mounted rigid wing that is supported by nonlinear translational and torsional springs and interacting with an incoming uniform stream is characterized. The objective is to discern the transition in the wing’s response between different states and equilibrium responses and to determine how the nonlinear response of the wing is impacted by variations in gust loads and in the nonlinear coefficients of the supporting springs. In these simulations, the flow and structure are treated as elements of a single dynamical system. The unsteady aerodynamic loads are predicted with a three-dimensional version of the unsteady vortex lattice method. A strong coupling scheme that is based on Hamming’s fourth-order predictor corrector method [11] is used to determine the structural response. Details of implementation of both aerodynamic simulation and predictor corrector are presented in the following sections. Bifurcations and transitions in the response of the wing as the wind speed is increased are discussed in the results sections. The effects of introducing a gust in the incoming flow and varying the nonlinear coefficients of the supporting springs are also discussed. A summary and conclusions are presented in the last section. 2. Model of the aeroelastic system 2.1. Structural model The elastically mounted wing considered in this work (see Fig. 1) is free to rotate (pitch motion) about an elastic axis and translate vertically (plunge motion). The system’s parameters and their values are presented in Table 1. Denoting by h and a the plunge deflection and pitch angle, respectively, one can write the governing equations of this system as ! ! ! ! !    0 kh ðhÞ ch 0 mT mW xa b L h h_ h€ ¼ þ þ , (1) 0 ka ðaÞ 0 ca Ia mW xa b M a a_ a€ where mT is the total mass of the wing and its supporting structure; mW is the wing mass alone; Ia is the mass moment of inertia about the elastic axis; b is half-chord length; xa is the nondimensional distance between the center of mass and the elastic axis; ch and ca are the plunge and pitch viscous damping coefficients, respectively; L and M are the aerodynamic lift and moment about the elastic axis, and kh and ka are the structural stiffness for the plunge and pitch motions, respectively. The representative parameters of this stiffness are approximated in polynomial form by kh ðhÞ ¼ kh0 , ka ðaÞ ¼ ka0 þ ka2 a2 ,

Fig. 1. Schematic of the elastically mounted wing.

(2)

5776

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

Table 1 Parameters of the aeroelastic system. Span c ¼2b xa

rp

Wing span (m) Wing chord (m) Nondimensional distance between center of gravity and elastic axis Air density (kg/m3)

1.3 0.27 0.331 1.225

mw mT Ia ca ch ka0 kh0 ka2

Mass of the wing (kg) Mass of wing and supports (kg) Moment of inertia about the elastic axis (kg m2) Pitch damping coefficient (kg m2/s) Plunge damping coefficient (kg/s) Stiffness in pitch (N m) Stiffness in plunge (N/m) Nonlinear stiffness in pitch (N m)

2.049 12.387 0.0558 0.036 27.43 6.833 2844.4 500

Fig. 2. Wing discretization and an example of a typical vortex element. The points P1, P2, P3, and P4 form the corners of the element. The control point (centroid of the element) Pc is determined by taking the intersection of the vectors linking the corner points. Gs and Gc denote the circulations around a vortex segment along the spanwise and chordwise directions, respectively.

where kh0 and ka0 are respectively the linear spring coefficients of the translational and torsional springs and ka2 is used to account for cubic stiffness of the torsional spring. 2.2. Aerodynamic model The unsteady flow over the elastically mounted wing is modeled using a three-dimensional unsteady vortex lattice method (UVLM). This method has been used extensively to determine aerodynamic loads and aeroelastic responses [9,12–16]. The UVLM applies only to incompressible and inviscid flows where the separation lines are known a priori. This tool is capable of simulating flows past moving thin wings and capturing the unsteady effects of the wake, but not the viscous effects, flow separation at the leading-edge, and extreme situations of strong wing–wake interactions. In UVLM, it is assumed that the flow field is inviscid everywhere except in the boundary layers and the wake. A set of discrete vortex rings is placed on the wing to represent a viscous boundary layer in the limit of infinite Reynolds number. Each vortex ring consists of four short straight vortex segments. The outcome is a lattice of vortex lines that divide the body surface into elements (see Fig. 2). The shape of each element is chosen to be rectangular. In previous studies [12,17], it was pointed out that such elements are superior to other shapes in terms of accuracy and rate of convergence. In the UVLM simulation, the position of, and circulation around, the vortex segments in the wake, and the circulations around the vortex segments on the wing are unknowns. The velocity field associated with these vortices is calculated according to Biot–Savart law. In the real flow, the velocity associated with the vorticity in the boundary layer and wake interacts with the incoming stream in a way that the no-penetration condition is satisfied on the surface of the wing. Thus, imitating the real situation, we solve for the circulations of the bound vortex segments by imposing the no-penetration condition at the centroid of each element, the so-called control point. The vorticity is introduced into the wake by shedding vortex segments from the trailing edge. To render the wake forcefree, these vortices are moved with the fluid particle velocity and their individual circulations remain constant. These vortices influence the flowfield around the oscillating wing. Because they were generated previously in the motion, the flow around, and hence the aerodynamic loads on the wing, depend on the history of the motion, unlike other aerodynamic

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

5777

modeling tools like blade-element theory. The flow is symmetric about the center chord plane of the wing. At t ¼ 0, there is no wake vorticity yet. At t ¼ Dt, the vorticity starts to shed from the trailing-edge. The same shedding mechanism is carried out for the next time steps. This process is repeated for any desired number of steps. Further details on the UVLM implementation are provided in Ref. [18,19]. The wakes created behind the elastically mounted wing for two different flow configurations are plotted in Fig. 3. In Fig. 3(a), the wing undergoes a static deflection (i.e., fixed angle of attack). It is noted that the roll-up of the wake vortices does not change its orientation. In Fig. 3(b), the wing oscillates freely (i.e., the angle of attack changes sign). In this case, the roll of the wake vortices alternates between positive and negative orientations. To compute the aerodynamic loads, we follow the approach of Katz and Plotkin [13,18]. This approach accounts for the leading-edge suction force. First, we compute the local lift and drag contributions of each element. Then, we transform the elemental lift and drag to the inertial frame. We then sum the elemental forces to compute the total lift and moment. To demonstrate the capability of UVLM of appropriately capturing the flow around the elastically mounted wing and accurately predicting the associated aerodynamic loads, we perform simulations over different configurations for fixed and plunging wings. Fig. 4 compares the lift coefficient as a function of the angle of attack as obtained using the 3-D UVLM code with values measured by Spearman and Feigh [20] and Magill [21]. The wing has a rectangular planform with an aspect ratio of 1.7 and a NACA 23012 airfoil. The numerical results obtained from UVLM are slightly lower than the experimental results. However, the slopes are about the same. Fig. 5 shows a comparison of the predicted lift from simulations using the 3-D UVLM on a plunging rectangular wing that has an aspect ratio of four and is placed at a constant angle of attack, a ¼ 51, with coefficients obtained by Katz and Plotkin [18]. The plunging amplitude is set equal to 0.1 c. Three different reduced frequencies (k ¼ o  c=2=U where o is the plunging frequency, c is the wing chord, and U is the freestream velocity) are considered. Four chordwise and 13 spanwise panels have been used in both simulations. As expected, the variations in the aerodynamic lift using the current simulations agree very well with those obtained by Katz and Plotkin [18].

Z

Z

X

X

Y

Y

Fig. 3. UVLM simulation of the wake of the elastically mounted wing for different configurations (a) fixed wing at a constant angle of attack and (b) oscillating wing.

Experimental results: Spearman and Feigh (1998)

0.6

Experimental results: Magill (2000) Current UVLM simulations

0.5

CL

0.4 0.3 0.2 0.1 0 −0.1

0

2

4

6

8

α (deg) Fig. 4. Experimentally (Spearman [20] and Magill [21]) and numerically (current UVLM simulations) obtained lift coefficient, C L , as function of the angle of attack a.

5778

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

0.1 0

κ = 0.5

−0.1

κ = 0.3

CL

−0.2

U∞

α

κ = 0.1

−0.3 −0.4 −0.5 −0.6 −0.7 0

1

2

3

ω.t

4

5

6

7

Fig. 5. Variations of the lift coefficient for a plunging wing in a uniform stream. Comparison of current results – solid line with numerical results of Katz and Plotkin [18] (UVLM) – squares.

2.3. Coupling formulation and scheme The most significant complication in solving for the aeroelastic response of a wing such as the one considered here is the fact that the prediction of the aerodynamic loads requires the wing’s motion and vice versa. The approach followed in this paper is to treat the flow and the wing as elements of a single dynamical system and to integrate all of the governing equations numerically, simultaneously, and interactively in the time domain. There is a fundamental complication with this approach: to predict the aerodynamic loads one must know the motion of the wing, and to predict the motion of the wing one must know the aerodynamic loads. To overcome this complication, an iterative scheme that accounts for the interaction between the aerodynamic loads and the motion of the wing is employed. The procedure is based on Hamming’s fourth-order predictor–corrector method [11]; and hence the added mass, the added damping, and/or added stiffness will not appear explicitly in the equations of motion. The Hamming’s method requires the solution to be known at three previous time steps and the current one. As a result, different methods are used to generate the solution for the first three time steps. For the first time step, Euler and modified Euler methods are used as predictor–corrector schemes. For the second time step, Adams–Bashforth two-step predictor and Adams–Moulton two-step corrector schemes are used. For the third time step, Adams–Bashforth three-step predictor and Adams–Moulton three-step corrector schemes are used. All of the above predictor–corrector schemes solve a system of first-order differential equations. For the sake of subsequent analyses, we define the state variables 2 3 2 3 X1 h 6 7 6 X2 7 6 h_ 7 6 7 7 X¼6 (3) 7 6 X3 7 ¼ 6 4 5 4a5 X4 a_ and rewrite the equations of motions in the first-order differential form as X_1 ¼ X 2 Ia kh0 Ia c mw xa bka0 mw xa bca mw xa bka2 3 Ia mw xa b M X1 h X2 þ X3 þ X4 þ X 3  L X_2 ¼  d d d d d d d X_3 ¼ X 4 m x bk m x bc m k mT c a mT ka2 3 mw xa b mT w a w a T a0 h0 h Lþ X_4 ¼ X1 þ X2  X3 X4 X3 þ M, d d d d d d d

(4)

where d ¼ mT Ia ðmw xa bÞ2 . In order to describe the numerical procedure for solving the equations of motion, we put them in the following set of system of first-order differential equations: y_ j ðtÞ ¼ f j ðt,y1 ðtÞ,y2 ðtÞ,y3 ðtÞ,y4 ðtÞÞ

for j ¼ 1,2,3, and 4,

(5)

with initial conditions yj ð0Þ ¼ aj

for j ¼ 1,2,3, and 4,

(6)

where yj, fj, and aj are the j-components of column vectors with four components each. Assuming that the solutions at the three previous time steps are computed using the special algorithms mentioned earlier, the right-hand side of Eq. (5) is

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

5779

evaluated at the present time step ti and the solution is predicted using Hamming’s predictor ypj ðt i þ DtÞ ¼ yj,4 þ

4Dt ½2f j,1 f j,2 þ 2f j,3  3

for j ¼ 1,2,3, and 4,

(7)

where yj,4 ¼ yj ðt i 3DtÞ f j,1 ¼ f j ðt i ,y1 ðt i Þ,y2 ðt i Þ, . . .Þ f j,2 ¼ f j ðt i Dt,y1 ðt i DtÞ,y2 ðt i DtÞ, . . .Þ f j,3 ¼ f j ðt i 2Dt,y1 ðt i 2DtÞ,y2 ðt i 2DtÞ, . . .Þ:

Fig. 6. Time histories of (a) the pitching angle and (b) the plunging displacement for different flow speeds.

Fig. 7. Bifurcation diagrams: variations of the static solution of (a) the pitching angle and (b) the plunging displacement with the flow speed.

Fig. 8. Bifurcation diagrams: variations of the dynamic solution (peak amplitudes of oscillations) of (a) the pitching angle and (b) the plunging displacement with the flow speed.

5780

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

The values of y’s and f’s are updated as yj,4 ¼ yj,3 ;

yj,3 ¼ yj,2 ;

f j,3 ¼ f j,2 ;

and

yj,2 ¼ yj,1 ; f j,2 ¼ f j,1 :

The predicted solutions, ypj ðt i þ DtÞ, are modified by using the truncation-error estimates from the previous time step (ti), simply ynj ðt i þ DtÞ ¼ ypj ðt i þ DtÞ þ

112 tej ðt i Þ 9

for j ¼ 1,2,3, and 4:

(8)

All of the components of tej ðt i Þ are set equal to zero for the first time the Hamming scheme is used. The solutions of Eq. (5) are set equal to the modified solutions computed by Eq. (7); namely 0

yj ðt i þ DtÞ ¼ ynj ðt i þ DtÞ

for j ¼ 1,2,3, and 4:

(9)

Fig. 9. Time histories (a, e, i, m), phase portraits (b, f, j, n), power spectra (c, g, k, o), and Poincare´ sections (d, h, l, p) of the pitching angle when U¼ 24 m/s, U ¼25 m/s, U ¼ 26 m/s and U ¼27 m/s, respectively.

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

5781

The left superscript 0 denote the number of iteration. At this stage (end of predictor section), the number of iteration is set equal to zero. At this level, the responses of the wing are found and supplied to the aerodynamic model in order to compute the aerodynamic loads at the new time step t i þ Dt. The right-hand sides of Eq. (5) are re-evaluated at the new time step t i þ Dt. The number of iterations k, at the beginning of the corrector section, is set equal to one. Hamming’s corrector is used to compute the corrected solutions, k yj,1 ðt i þ DtÞ; namely k

yj,1 ðt i þ DtÞ ¼ 18½9yj,2 yj,4 þ 3Dtðf j,1 þ 2f j,2 f j,3 Þ

for j ¼ 1,2,3, and 4,

(10)

where yj,2 ¼ yj ðt i Þ yj,4 ¼ yj ðt i 2DtÞ f j,1 ¼ f j ðt i þ Dt, k1 y1 ðt i þ DtÞ, k2 y1 ðt i þ DtÞ, . . .Þ f j,2 ¼ f j ðt i ,y1 ðt i Þ,y2 ðt i Þ, . . .Þ f j,3 ¼ f j ðt i Dt,y1 ðt i DtÞ,y2 ðt i DtÞ, . . .Þ f j,4 ¼ f j ðt i 2Dt,y1 ðt i 2DtÞ,y2 ðt i 2DtÞ, . . .Þ: The absolute value of the difference between the solutions at the present and previous iteration, k ej , is estimated as follows: k

ej ¼ 9k yj k1 yj 9

for j ¼ 1,2,3, and 4:

(11)

k

The values of ej are checked to determine if all its components are less than a specified error tolerance or not. If the condition is not satisfied, the computation process is repeated starting at the top of the corrector section of the Hamming method. The iteration number, k, is increased by one (kþ1). The present iteration number is k while the previous iteration is set equal to k1. The right-hand sides of Eq. (5) are re-evaluated at time step t i þ Dt using the corrected solution, k1 yj,1 ðt i þ DtÞ.

10−2

d

10−4

10−6

10−8

0

5

10

15

Time (s) Fig. 10. Sensitivity to initial conditions for the aperiodic attractor at U ¼27 m/s.

0

0.18

2

4

6

8

10

12

14

16

18

0.16

20

0.02 0.01

0.14 0 −0.01

0.1 0.08

−0.02

0.06

−0.03

h(m)

α (rad)

0.12

0.04 −0.04 0.02 0

−0.05 0

2

4

6

8

10 12 Time (s)

14

16

18

20

Fig. 11. Time histories of the pitching (blue line) and plunging (black line) motions at U ¼ 32 m/s. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

5782

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

Once more, Hamming’s corrector is used to compute the corrected solutions, k yj,1 ðt i þ DtÞ using Eq. (11). The absolute value of the difference between the solutions at the present and previous iteration, k ej , is estimated using Eq. (11) and its values are checked to determine if all its components are less than the error tolerance. The computation processes are repeated for several iterations until the solutions converge or the number of iterations k reach certain value at which the solution is considered to be diverged. After the solutions convergence, the truncation error, tej ðt i þ DtÞ is estimated as follows: tej ðt i þ DtÞ ¼

9 k ½ y ðt þ DtÞypj ðt i þ DtÞ 121 j i

for j ¼ 1,2,3, and 4:

(12)

The corrected solutions, k yj ðt i þ DtÞ are modified by using the new truncation-error estimate, tej ðt i þ DtÞ, k

yj ðt i þ DtÞ ¼ k yj ðt i þ DtÞtej ðt i þ DtÞ

t1

Ugust

t2

tn

for j ¼ 1,2,3, and 4:

plate chord

U

xG Fig. 12. Deterministic gust representation.

Fig. 13. Effects of the gust on the pitch and plunge motions for U ¼19 m/s and when (a) U gust ¼ 5%U, (b) U gust ¼ 8%U and (c) U gust ¼ 10%U.

(13)

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

5783

Finally, the time is set equal to t i þ Dt. The process may be repeated for the next time step t i þ 2Dt and the new truncation error, tej ðt i þ DtÞ, which is computed at this time level, t i þ Dt, will be used in the next time level, t i þ 2Dt. 3. Results and discussion Time histories of the pitching and plunging motions of the wing for different flow speeds are presented in Fig. 6. For U r 19 m=s, the wing undergoes initial motions and settles back to its original position of zero pitch and plunge. When the flow speed becomes larger than U¼19 m/s, the wing moves to a static equilibrium position where both pitch and plunge have constant values that are different than zero. The amplitudes of both pitch and plunge deflections increase as the flow speed is increased. At U ¼22 m/s, the wing starts oscillating around as ¼ 0:057 rad and hs ¼ 0:0058 m. The amplitude of the oscillations as well as the mean values increase as the flow speed is increased. It is obvious that the wing exhibits different equilibrium conditions as the speed of the incoming flow is varied. To gain more insight into the static and dynamic behaviors of the wing when varying the flow speed U, the bifurcation diagrams showing variations of the static and dynamic solutions of the pitching and plunging motions are plotted in Figs. 7 and 8, respectively. In Fig. 7, the solid and dashed lines correspond to the stable and unstable equilibrium points, respectively. As U is increased, the zero solution remains stable until a critical value of U ¼ U s ¼ 19 m=s is reached. At this speed, a static bifurcation (supercritical pitchfork bifurcation) occurs. For U 4U s , two symmetric stable solutions are noted. These two nontrivial solutions lose stability at U¼21.4 m/s. At this speed, the wing undergoes a Hopf bifurcation leading to limit cycle oscillations (LCO). It is important to note that as the speed is increased beyond 21.4 m/s, the mean value of the pitch motion decreases while that of the plunge motion increases. This indicates that the energy distribution among the two different motions is dependent on the speed. It should also be noted that the upper and lower branches of the bifurcation diagrams are obtained by considering different initial conditions. Fig. 8 shows the variations of the peak amplitudes of the pitching and plunging oscillations at speeds between 18 m/s

Fig. 14. Effects of the gust on the pitch and plunge motions for U¼ 21 m/s and when (a) U gust ¼ 5%U, (b) U gust ¼ 10%U and (c) U gust ¼ 13%U.

5784

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

and 25 m/s. Clearly, the wing undergoes a supercritical Hopf bifurcation near U¼21.4 m/s that is exhibited through the smooth transition from static to oscillating solutions. Fig. 9 shows time histories, phase portraits, power spectra, and Poincare´ sections of the pitch motion for U¼24, 25, 26, and 27 m/s. At U ¼ 24 m/s, the time history, phase portrait, and Poincare´ section presented in Fig. 9(d) shows that the response is periodic. Furthermore, we observe that the power spectrum of the pitch angle exhibits a peak near f ¼ f N  2:2 Hz and smaller peaks at its second and third harmonics. These peaks originate from the cubic nonlinearity present in the structural model and the static deflection that gives rise to a quadratic nonlinearity in the presence of a cubic nonlinearity. As the flow speed is increased to U¼25 m/s, the motion is characterized by a period doubling. This is shown in the two loops of the phase portrait and the appearance of a subharmonic in the power spectrum. The Poincare´ section, Fig. 9(h), shows a finite number of discrete points indicating that the motion remains periodic. At U¼ 26 m/s, another period doubling takes place. This is exhibited in the four loops of the phase portrait and by the second subharmonic of the fundamental frequency. Again, the finite number of discrete points in the Poincare´ section in Fig. 9(l) indicates that the flow is periodic. Beyond this sequence of period-doubling bifurcations, the motions are characterized as chaotic. This characterization is based on the aperiodic attractor at U¼27 m/s and is exhibited by the unorganized set of points in the Poincare´ section indicating a chaotic motion. To check further if this attractor is chaotic, its sensitivity to a small perturbation in the initial conditions is investigated. As such, two initial points in the plunge motion, separated by d0 ¼ 106 , are considered and their time evolutions are recorded. The variations of the separation distance d with time are plotted in Fig. 10. Clearly, the separation grows exponentially in the range 7 s o time o15 s. Thus, there is a positive Lyaponov exponent g associated with this attractor. It should be noted here that when the flow speed is increased to U ¼32 m/s the wing experiences initial oscillations and then settles down to a static equilibrium (as ¼ 0:157 rad and hs ¼ 0:032 m) as shown in Fig. 11. As such, the region of chaotic response is followed by a region of static equilibrium.

Fig. 15. Effects of the gust on the pitch and plunge motions for U ¼22 m/s and when (a) U gust ¼ 5%U, (b) U gust ¼ 8%U and (c) U gust ¼ 10%U.

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

5785

4. Effects of parameter variations on wing response The different regions of aeroelastic response and transition from one state to another must depend on the system’s parameters which include the flow speed, as shown above, and the wing’s structural parameters. Because the different bifurcations and response amplitudes are associated with these parameters, variations in any of them could lead to significant variations in the wing’s response. In this section, we investigate the effects of introducing a gust in the incident flow and varying the nonlinear stiffness coefficient on the system’s response. 4.1. Gust effects To determine the gust effects, we add a gust that is represented by a one cosine profile to the mean freestream to yield U fst ¼ U þU G , where U G is the gust velocity and is given by [22]    8 U gust 2px0 > > 1cos < 2 xG UG ¼ > > :0

xG , U xG for t 4 , U

for t o

x0 o 0 (14) x0 4 0,

where x0 ¼ xUt and xG denotes the spatial wavelength which is set equal to 0.4 c. Fig. 12 illustrates the progression of a selected gust in space and its position in relation to the wing. The gust effects are determined by considering different values for U gust , as a percentage of U, over each of the response states as observed above.

Fig. 16. Effects of the gust on the pitch and plunge motions for U¼ 24 m/s and when (a) U gust ¼ 3%U, (b) U gust ¼ 5%U and (c) U gust ¼ 10%U.

5786

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

Time histories of the pitching and plunging motions for different gust amplitudes when U¼ 19 m/s, i.e. in the stable region prior to the static pitchfork bifurcation, are shown in Fig. 13. Clearly, the gust does not affect the response of the system as both motions are damped for all gust values. At U¼21 m/s, which is in the region of static equilibrium pitchfork bifurcation, we note that for both the pitch and plunge motions, the gust affects the system’s response only when its value is large, about U gust ¼ 13%U. Furthermore, we note that the static equilibrium points of the pitch and plunge responses change from positive to negative value and inversely. Beyond the point of the Hopf bifurcation, two different cases are considered. At U¼22 m/s, the gust changes the system’s response from oscillating about positive mean pitch and plunge values to one that oscillates about negative mean values for gust amplitudes of about 10%. At a higher speed U¼24 m/s, similar variations in the response are observed but at lower gust amplitudes. Clearly, the gust impacts the response of the wing by changing the mean values of the pitch and plunge motions without affecting the type of equilibrium over different regions. Furthermore, Figs. 14–16 shows that the gust amplitude which impacts the response are smaller as the flow speed is increased.

Fig. 17. Time histories of (a) the pitching angle and (b) the plunging displacement for different values of the nonlinear pitching stiffness coefficient ka2 .

Fig. 18. Variations of the mean values of (a) the pitching angle and (b) the plunging displacement with the nonlinear pitching stiffness coefficient ka2 for U ¼ 22 and 24 m/s.

Fig. 19. Variations of the LCO amplitudes of (a) the pitching angle and (b) the plunging displacement with the nonlinear pitching stiffness coefficient ka2 for U ¼22 and 24 m/s.

M. Ghommem et al. / Journal of Sound and Vibration 331 (2012) 5774–5787

5787

4.2. Effects of structural nonlinearities Time histories of the pitching and plunging motions at U¼22 m/s and for different values of ka2 are plotted in Fig. 17. Clearly, increasing ka2 reduces both of the static deflections and amplitudes of limit cycle oscillations (LCOs). In Fig. 18, the variations of the mean value solutions of the pitching and plunging motions with ka2 for U ¼22 and 24 m/s are plotted. Increasing the value of ka2 decreases the static pitch angle and plunge displacement and the slope of their variations as well. Similar observations can be made for the variations of the LCO amplitudes of the pitching and plunging motions with ka2 as shown in Fig. 19. We note also that the change is much more pronounced at the larger speeds. 5. Conclusions In this work, an elastically mounted wing free to plunge and pitch, supported by nonlinear translational and torsional springs, and interacting with an incoming uniform stream has been analyzed. The results show that the wing undergoes a sequence of static and dynamic bifurcations as the wing speed is increased. The dynamic bifurcations include a sequence of period doubling that develop into a chaotic motion as the speed increased. The effects of the gust and structural nonlinearity on the wing aeroelastic response have also been investigated. The results show that gust may lead to a jump between pairs of stable solutions over the different response regimes. This is a manifestation of the system’s nonlinearities. Of interest are two notions. The first is that the gust does not impact the boundaries of the different response regimes. The second is that the gust amplitudes at which the wing’s response is smaller at the higher wind speeds. As for the effect of the structural nonlinearity, increasing the nonlinear coefficient of the stiffness of the torsional spring have yielded lower static deflections and LCO amplitudes. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

A. Abdelkefi, A.H. Nayfeh, M.R. Hajj, Modeling and analysis of piezoaeroelastic energy harvesters, Nonlinear Dynamics 67 (2012) 925–939. A. Abdelkefi, A.H. Nayfeh, M.R. Hajj, Design of piezoaeroelastic energy harvesters, Nonlinear Dynamics 68 (2012) 519–530. A. Abdelkefi, A.H. Nayfeh, M.R. Hajj, Enhancement of power harvesting from piezoaeroelastic systems, Nonlinear Dynamics 68 (2012) 531–541. E.H. Dowell, D. Tang, Nonlinear aeroelasticity and unsteady aerodynamics, AIAA Journal 40 (2002) 1697–1707. A. Abdelkefi, R. Vasconcellos, F.D. Marques, M.R. Hajj, Bifurcation analysis of an aeroelastic system with concentrated nonlinearities, Nonlinear Dynamics 69 (2011) 57–70. C. Chabalko, M.R. Hajj, W. Silva, Interrogative testing for nonlinear identification of aeroelastic systems, AIAA Journal 46 (2008) 2657–2658. A. Abdelkefi, R. Vasconcellos, F.D. Marques, M.R. Hajj, Modeling and identification of freeplay nonlinearity, Journal of Sound and Vibration 331 (2012) 1898–1907. H.C. Gilliatt, T.W. Strganac, A.J. Kurdila, An investigation of internal resonance in aeroelastic systems, Nonlinear Dynamics 31 (2003) 1–22. A.O. Nuhait, D.T. Mook, Aeroelastic behavior of flat plates moving near the ground, Journal of Aircraft 47 (2010) 464–474. Z. Wang, P.C. Chen, D.D. Liu, D.T. Mook, M.J. Patil, Time domain nonlinear aeroelastic analysis for hale wings, Proceedings of the 47th AIAA/ASME/ASCE/ AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA Paper No. 2006-1640, 2006. H.A.L.B. Carnahan, J.O. Wilkes, Applied Numerical Methods, Wiley, New York, 1969. Z. Wang, Time-Domain Simulations of Aerodynamic Forces on Three-Dimensional Configurations, Unstable Aeroelastic Responses, and Control by Network Systems, PhD Thesis, Virginia Tech, Blacksburg, VA, 2004. B.K. Stanford, P.S. Beran, Analytical sensitivity analysis of an unsteady vortex-lattice method for flapping-wing optimization, Journal of Aircraft 47 (2010) 647–662. M.S. Vest, J. Katz, Unsteady aerodynamic model of flapping wings, AIAA Journal 34 (1996) 1435–1440. T.E. Fritz, L.N. Long, Object-oriented unsteady vortex lattice method for flapping flight, Journal of Aircraft 41 (2004) 1275–1290. R. Palacios, J. Murua, R. Cook, Structural and aerodynamic models in nonlinear flight dynamics of very flexible aircraft, AIAA Journal 48 (2010) 2648–2658. S. Preidikman, Numerical Simulations of Interactions among Aerodynamics, Structural Dynamics, and Control Systems, PhD Thesis, Virginia Tech, Blacksburg, VA, 1998. J. Katz, A. Plotkin, Low-Speed Aerodynamics, Cambridge University Press, Cambridge, MA, 2001. M. Ghommem, M.R. Hajj, D.T. Mook, B.K. Stanford, P.S. Beran, R.D. Snyder, L.T. Watson, Global optimization of actively morphing flapping wings, Journal of Fluids and Structures 33 (2012) 210–228. M.L. Spearman, K.M. Feigh, An airplane configuration with an inboard wing mounted between twin fuselages, AIAA Paper No. 1998-440-695, 1998. S. Magill, Progress report on stability tunnel and computer panel method of large civil transport, Technical Report, Department of Aerospace and Ocean Engineering, Virginia Tech, 2000. E.A. Dowell, H.C. Curtiss, R.H. Scanlan, F. Sisto, A Modern Course in Aeroelasticity, Kluwer Academic Publishers, The Netherlands, 1989, p. 119.