A comparative study of nonlinear aeroelastic models for high aspect ratio wings

Journal of Fluids and Structures 85 (2019) 249–274

Contents lists available at ScienceDirect

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

A comparative study of nonlinear aeroelastic models for high aspect ratio wings ∗

Amir Hossein Modaress-Aval a , Firooz Bakhtiari-Nejad a , , Earl H. Dowell b , David A. Peters c , Hossein Shahverdi d a

Department of Mechanical Engineering, Amirkabir University of Technology, Tehran, Iran Department of Mechanical Engineering, Duke University, Durham, NC, United States c Department of Mechanical Engineering, Washington University St. Louis, Mo, United States d Department of Aerospace Engineering, Amirkabir University of Technology, Tehran, Iran b

article

info

Article history: Received 8 August 2018 Received in revised form 1 December 2018 Accepted 7 January 2019 Available online xxxx Keywords: Aeroelastic Nonlinear High aspect ratio wing Subsonic 2D and 3D aerodynamic models

a b s t r a c t Linear and nonlinear aeroelastic analyses and results for a cantilevered beam subjected to an unsteady subsonic airflow are presented as a model of a high aspect ratio wing. A third-order nonlinear beam model is used as the structural model. This model includes nonlinear inertial and damping terms as well as nonlinear stiffness terms. To model aerodynamic loads a Wagner model based on a state-space model along with 2D and 3D Peters’ aerodynamic models are used and the results are compared. The Galerkin method is used to transform partial differential equations to ordinary differential equations. For aeroelastic analysis of a high aspect ratio wing, using a state-space Wagner model reduces the number of degrees of freedom. The distribution of the lift and inflow derived using 2D and 3D Peters aerodynamic models are compared and the range of aspect ratios has been determined for which a 2D aerodynamic model is sufficient. There is a good agreement between the flutter speed obtained in this study and those reported in the literature. Beside nonlinear stiffness terms, the effects of nonlinear inertial and damping terms may be significant for flow velocities above the flutter speed. As the flow velocity increases further the oscillations lose their periodicity and become chaotic. Using the different number of structural modes, the convergence of the results is shown. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction High-aspect-ratio wings are becoming more important in the design of modern flying vehicles, e.g. high-altitude long-endurance (HALE) aircraft. These aircraft are being considered for military reconnaissance, long-term surveillance, environmental sensing, and communications relay missions. As a result of weight restrictions, the wings are very flexible and this may lead to large deformations. So the structural characteristics and thus aeroelastic response of these wings may be nonlinear. Thus nonlinear structural as well as aerodynamic modeling is important for the analysis and design of these flexible structures. Additionally, a simple yet realistic model is desired that can be used in designing a control system. Aeroelastic analysis of high-aspect-ratio wings has been studied by many researchers. Most studies have used linear beam theory to simplify the wing structural model and did not consider the geometrical nonlinearities due to large deflections. However, during the last two decades, there has been a growing interest in the nonlinear aeroelastic analysis of high aspect ∗ Corresponding author. E-mail address: [email protected] (F. Bakhtiari-Nejad). https://doi.org/10.1016/j.jfluidstructs.2019.01.003 0889-9746/© 2019 Elsevier Ltd. All rights reserved.

250

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Nomenclature a b bi EIi e GJ ix iy iz i1 i2 i3 M qi R U∞ u

v T XYZ x y, z Z

αi βi γi δ W nc ϵij ξ ηζ θi λi π ρ ρi ρ∞ σij ψ, θ , φ

Nondimensional distance between mid-chord (M.C) and elastic axis (E.A) (the normalizing factor is b) Reference semi-chord of the lifting surface Aerodynamic coefficients (Peters’ aerodynamic model) Bending rigidity, i = 1, 2 Axial strain Torsional rigidity Unit vectors along the x, y and z axes Unit vectors along ξ , η and ζ axes Number of state variable in Peters’ model Distributed external loads, i = 2,4 Deformed position vector of an arbitrary point on the representative cross-section Freestream air velocity Axial deformation Displacement of the reference point on the observe cross-section along y Kinetic energy Cartesian coordinate system describing the undeformed geometry Spanwise coordinate centered at root of the wing, dimensionless on L Coordinate of an arbitrary point on the representative cross-section The number of spanwise elements (in 2D Peters’ model) Wagner function coefficients Wagner function coefficients Aerodynamic coefficients (Wagner aerodynamic model) Virtual work of the nonconservative forces Local strains Local Cartesian coordinate system associated with curvilinear coordinates Rotation about ii , i = 1, 2, 3 Aerodynamic states Elastic energy Density Curvature about ii Freestream air density Local stress Angle of rotation about Z , Y1 and ξ axes respectively (3-2-1)

ratio wings. Geometric structural nonlinearities may arise from the coupling among elastic flap bending, chordwise bending, and torsion. Patil et al. (1999) studied the effects of geometrical nonlinearities on aeroelastic behavior of high aspect ratio wings. In their study, a complete nonlinear mixed variational formulation is used for the structural dynamics of the wing and the two-dimensional (2D) Peters aerodynamic model is used for the loading. They characterized the effects of steadystate lift and drag. Two years later, they (Patil et al., 2001a) studied the nonlinear aeroelastic analysis of a high aspect ratio wing using same structural and aerodynamic models considering the effects of stall. In 2004 (Patil and Hodges, 2004), they used a geometrically exact structural model based on a complete nonlinear mixed variational formulation, with a nonplanar aerodynamic theory for investigating the effects of geometrical nonlinearities on aeroelastic behavior of high-aspect-ratio wings. In all of their studies, finite elements in space and time were used for time integration. Tang and Dowell (2001) did an experimental and theoretical study on the aeroelastic response of a high aspect ratio wing. They used structural equations of motion based on nonlinear beam theory which is combined with the ONERA aerodynamic stall model to study the effects of geometric structural nonlinearity and steady angle of attack on flutter and limit cycle oscillation (LCO) of high-aspect-ratio wings. In 2002, they (Tang and Dowell, 2002) did an experimental and theoretical study of gust response of high aspect ratio wing using the same structural and aerodynamic models. Tang and Dowell (2004) considered geometrically nonlinear effects using second-order geometric nonlinear terms with the ONERA aerodynamic stall model. They used a simpler nonlinear beam theory in comparison to the model used by Patil and Hodges, but consider the effects of dynamic stall on the aeroelastic response of high aspect ratio wing. Shams et al. (2008) considered a wing in vertical and torsional motion using the second-order form of nonlinear general flexible Euler–Bernoulli beam equations for structural modeling. Unsteady linear aerodynamic theory based on an integral

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

251

form of Wagner function is used for determination of aerodynamic loading on the wing. They showed there is a need for at least three structural modes for nonlinear aeroelastic response calculation beyond the flutter boundary. Jian and Jinwu (2009) did an aeroelastic analysis of a high aspect ratio wing taking into account structural geometrical nonlinearities and dynamic stall combining a geometrically exact, nonlinear anisotropic beam model with nonlinear ONERA (Edlin) dynamic stall model. Numerical results indicate that structural geometrical nonlinearities could lead to LCO without stall occurring. Harmin and Cooper (2011) presented an aeroelastic Reduced Order Model (ROM) for wing structures that experience geometric nonlinearities. The methodology consists of a combined modal/finite element approach coupled with unsteady doublet lattice method to generate a reduced order model capable of predicting the aeroelastic response of a flexible wing structure. In 2012 Ovesy et al. (2012) studied geometrically nonlinear effects in aeroelastic analysis of a high aspect ratio wing using nonlinear inextensional beam model with semi-empirical ONERA aerodynamic model. However, their model did not include nonlinear inertial and damping terms. In 2012 Xie et al. (2012) published an experimental study similar to that of Tang and Dowell. In this study, they used a finite element model coupled with the lifting line method for the static aeroelasticity and the strip theory considering curved surface effects for the flutter analysis. They investigated LCO behavior using time-marching solutions. In 2013, Arena et al. (2013) presented the nonlinear aeroelastic modeling and behavior of HALE wings, undergoing large deformations and exhibiting dynamic stall. A parametric structural model of wings based on an exact kinematic approach is formulated and coupled with an incompressible unsteady aerodynamic model that is obtained via an indicial formulation accounting for viscous effects, including dynamic stall and flow separation (a modified Beddoes–Leishman model is employed). The nonlinear aeroelastic equations are then implemented and solved by a Finite Element code. The results from the aeroelastic analysis confirm the trends for non-linear geometry effects on flutter speed. They analyzed the post-flutter behavior of HALE wings using their computational aeroelastic model. Bakhtiari-Nejad et al. (2017) used state-space Wagner model to model aerodynamic forces to do the aeroelastic analysis of a high aspect ratio wing. The large wing deformations observed in high aspect-ratio wings caused by the flexibility increase also change the flight mechanics of the aircraft and may become an issue. Most of these applications consider HALE (Hallissy and Cesnik, 2011). Patil et al. (2001b) published a study concerning HALE aircraft’s wings with high aspect ratios. Their theory was based on a mixed variational formulation based on the exact intrinsic equations for dynamics of beams in moving frames and a finitestate air load model for deformable airfoils on fixed and rotating wings. Results were obtained from such an analysis for aeroelastic behavior as well as overall fight dynamic characteristics of a complete aircraft model representative of HALE aircraft. A computational code, called NATASHA (Nonlinear Aeroelastic Trim and Stability of HALE Aircraft), was developed by Patil and Hodges (2006) to investigate the HALE aircraft. They presented a theory for flight-dynamic analysis of highly flexible flying-wing configurations. In this framework, the structure is modeled using the non-linear geometrically exact beam equations (Hodges, 2003) and the aerodynamic model consists of 2D Peters’ aerodynamic theory. To solve the set of final equations, the beam is discretized into finite elements. They showed that the flight-dynamic characteristics of the deformed aircraft are completely different as compared with a rigid aircraft. The authors of Castellani et al. (2017), Zhao and Ren (2011) and Krüger (2007) used a multibody dynamic approach for investigating the flight dynamic and nonlinear aeroelasticity of flexible aircraft. Chang et al. (2008) presented an analysis and parametric study of the flight dynamics of highly flexible aircraft to predict the atypical flight dynamic characteristics of highly flexible flying wings with one or more fuselages, wings, and/or tails. They used geometrically exact, intrinsic beam elements as the structural model which is coupled with an aerodynamic model consisting of 2D, large angle-of-attack, unsteady theory for the lifting surfaces. A central-difference discretization scheme is applied to the intrinsic governing equations in space to obtain a numerical solution. Sotoudeh et al. (2010) perform some validation studies (non-linear structural dynamics, structural dynamics and aeroelasticity) to evaluate the accuracy of employing the computation code of NATASHA (nonlinear aeroelastic trim and stability for HALE aircraft) which is based on the published geometrically exact, fully intrinsic beam equations to study HALE aircraft at a conceptual and preliminary stage. This study presented a wide range of results from NATASHA, comparing them with well-known solutions of beam stability and vibration problems, experimental data, or results from well-established computer codes. Mardanpour et al. (2013) investigated the effects of engine placement and sweep on the flutter characteristics of an aftswept flying wing using the NATASHA framework. Two years later, in another study (Mardanpour et al., 2014), Mardanpour and et al. did a nonlinear aeroelastic analysis of a high aspect ratio wing considering the effects of time-dependent thrust using the code NATASHA. Time-marching analyses were conducted using time-dependent thrust to simulate the non-linear aeroelastic response of a high aspect-ratio wing. They observed a stable LCO below flutter speed for an impulse engine excitation if the engine is placed at the wing tip and behind the elastic axis. The initial flight tests of an X-56A system, performed by Lockheed and AFRL in the latter half of 2013 and extending into early 2014, collected flight data on highly flexible structures and flutter suppression control technology. The X-56A is a small, fixed-wing aircraft that was designed to exhibit body freedom flutter, symmetric wing bending–twisting, and asymmetric wing bending–twisting modes within its flight envelope. The goal of The X-56A to advance aeroservoelastic technology through flight research using a low-cost, modular, remotely piloted aircraft. Palacios and his co-workers developed a multidisciplinary framework named SHARP (Simulation of High-Aspect-Ratio Planes) at the Imperial College London (Hesse and Palacios, 2014; Murua et al., 2012). The main goal of this framework is to enable the study of maneuvering flexible aircraft dynamic behavior (including aircraft and offshore wind turbines

252

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

applications) undergoing large deformations. The aeroelastic response of the vehicle, which may be subject to large wing deformations at trimmed flight, is captured by coupling a displacement-based flexible-body dynamics formulation with an aerodynamic model based on the unsteady vortex lattice method. This framework allows for computing: trim conditions; gust and maneuver loads; stability analysis; dynamic response; flutter suppression; and gust load alleviation. A non-linear aeroservoelastic framework was developed by Wang et al. (2016) to study flight simulations of flexible aircraft in modal coordinates. The study relied on the projection of the intrinsic equations of geometrically nonlinear composite beams on the linear normal modes at a reference condition. They made the use of a linear three-dimensional (3D) FE model to create such an intrinsic beam description of the structure, which is then coupled with a 2D unsteady aerodynamic model. Linear H∞ control synthesis and closed-loop nonlinear simulations were explored on a highly flexible flying wing under large-amplitude discrete gusts. The Wagner aerodynamic model has rarely been used to study aeroelastic analysis of high aspect ratio wings. One of the main objectives of the present study is to introduce a method using Wagner aerodynamic model which can reduce the number of degree of freedom compared with other aerodynamic models. This could be an asset in the design of control systems. To show the advantages of the current formulation, a careful comparison of two different 2D aerodynamic models has been done. In addition, the lift distribution obtained using 2D and 3D aerodynamic models have been presented and the range of aspect ratios has been determined for which a 2D aerodynamic model is sufficient to do an aeroelastic analysis. It has also been shown that besides nonlinear stiffness terms, the effects of nonlinear inertial and damping terms are significant for flow velocities above the flutter speed. In the present study, linear and nonlinear aeroelastic analyses of high aspect ratio wing using different aerodynamic models are presented. A third-order nonlinear beam model is used as the structural model and to model aerodynamic loads different incompressible aerodynamic models have been adopted: Wagner and Peters’ aerodynamic models. The study shows the number of degrees of freedom using the state-space Wagner model is less than the number of degrees of freedom using 2D and 3D Peters’ aerodynamic model. The equations are discretized using Galerkin’s method. First, the small perturbation flutter boundary is determined. Subsequently, LCO is considered when nonlinear effects on the aeroelastic behavior are significant. Using different number of structural modes, the convergence of the results is shown. 2. Aeroelastic model 2.1. Structural model Various methods have been used to study the aeroelastic behavior of a high aspect ratio wing. All of them have used a beam model to study the aeroelastic behavior of these wings. Here the governing equations of vibrations considering geometric nonlinearities have been derived. To fully account for geometric nonlinearities, local stress and strain measures have been used. The nonlinear equations of motion have been expanded into polynomials about the static equilibrium position. Since in this study the initial angle of attack is zero and also the cross section of the wing is rectangular (a symmetric airfoil), there is no static force. So it is assumed that the equilibrium position is close to the undeformed position. The nonlinear terms are kept up to third order. For a beam with one end at s = 0 being fixed or hinged, and the other end at s = L being free or sliding with no external load acting along the x direction, one may assume the beam is inextensional (i.e. the axial strain is zero but the axial deformation is not zero) (Nayfeh, 2004). It should be noted because the drag force is neglected, only bending and torsional vibrations have been considered and chordwise vibration is ignored. Although there are several works (Patil et al., 2001a; Tang and Dowell, 2004; Ovesy et al., 2012) in which the drag force and in-plane bending have been considered, for the cases studied in the present paper, there is no significant difference in flutter speed found in this study and that in Patil et al. (2001a), Tang and Dowell (2004) and Ovesy et al. (2012). In addition, in this study, the initial angle of attack is zero and consequently, the drag force is neglected as is chordwise bending. With these two assumptions (inextensional beam and ignoring chordwise bending), the governing equations of motions reduce to two equations, one for flapwise motion and another for twisting motion. The governing equations of vibrations can be written as follows (for more details see Appendix A and chapter 4 of Nayfeh, 2004) mv¨ + D33 v

′′′′

= −mv

′′

[∫ s ∫ L

− mv ′

s

(

v˙ + v v¨ dsds 2

′ ′

)

]

0

s

∫

(

) v˙ ′2 + v ′ v¨ ′ ds

0

[ − D33

4v ′ v ′′ v ′′′ + v ′2 v ′′′′ + v ′′3 − 2φφ ′′ v ′′ −2φ ′2 v ′′ − 4φφ ′ v ′′ − φ 2 v ′′′′

]

(1)

[ ] − D22 2φ ′2 v ′′ + 2φφ ′′ v ′′ + 4φφ ′ v ′′ + φ 2 v ′′′′ + q2 + j3 v¨ ′′ j1 φ¨ − D11 φ = (D33 − D22 ) v ′′2 φ + q4 ′′

(

)

(2)

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

253

Also since the beam is inextensional we have 1 u′ = − v ′ 2 2

(3)

It should be noted that the first two integral terms in Eq. (1) are nonlinear inertial and damping terms which in most previous studies have been ignored. The remaining terms are nonlinear stiffness terms. All the parameters are defined in Appendix. In the next section, aerodynamic models used in this study are described. 2.2. Aerodynamic models Wagner and Peters’ aerodynamic models have been used. These models are based on the assumption of incompressible potential flow and consequently the effects of viscosity have been neglected. For the rest of the paper, the state-space Wagner aerodynamic model and the Peters aerodynamic model are simply called the Wagner model and the Peters’ model respectively. 2.2.1. Wagner model In the Wagner aerodynamic model, the lift and moment , q2 and q4 , are as follows

¨ q2 = L = ρ∞ π b2 [−¨v + U∞ φ˙ − baφ]

[

+ π bρ∞ U∞ −˙v + U∞ φ + b

(

1 2

) ] − a φ˙

(4)

− 2π b2 ρ∞ U∞ [γ1 λ1 + γ2 λ2 ] q4 = Me.c =

( [ ( ) ) ] 1 1 ρ∞ π b3 −av¨ + a − U∞ φ˙ − b a2 + φ¨ 2 8 ( )[ ( ) ] 1 1 2 + π b ρ∞ U ∞ a + −˙v + U∞ φ + b − a φ˙ 2 2 ( ) 1 [γ1 λ1 + γ2 λ2 ] − 2π b2 ρ∞ U∞ a +

(5)

2

where ρ∞ and U∞ are freestream air density and freestream air velocity, respectively. b is a reference semi-chord of the lifting surface and a is a nondimensional distance between the mid-chord (M.C) and the elastic axis (E.A), as can be seen in Fig. 1. It should be noted that if E.A. is located between the leading edge and the M.C, it is negative and if it is located on the opposite side (between trailing edge and M.C), it is positive. λi and γi are aerodynamic states and aerodynamic coefficients respectively which are as follows

λi =

∫

t

e−

βi U b

(t −τ )

w (τ ) dτ

(6)

0

γi = αi βi

U

(7)

b

where αi and βi are Wagner function coefficients and w (t) is

w (t ) = −˙v + φ U∞ + b

(

1 2

) −a φ

(8)

For a beam where φ and v varies along the beam, we have

w = w(t , x)

(9)

It can be shown that λi s are governed by differential equations as follows:

λ˙ i +

βi U b

λi = w (t)

(10)

From Eqs. (6) and (9), it can be concluded that

λi = λi (t , x)

(11)

The dependence of λi on x is an important consideration when Galerkin method is used to solve the final equations as will be discussed later. In the Galerkin method, the equations of motion will be multiplied by mode shapes and then will be integrated over the wing span.

254

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 1. Schematic showing geometry of the wing section.

When Eq. (10) is multiplied by the mode shapes and integrated from 0 to L, because of the dependence of λi to x, new variables Yi and Xi will be introduced respectively as shown in the Eqs. (12) and (13).

( ) βi U βi U ψi λ˙ i + λi dx = Y˙i + Yi b b 0 ) ∫ L ( βi U βi U θi λ˙ i + λi dx = X˙ i + Xi

∫

L

b

0

(12) (13)

b

Similarly, in the structural bending and pitching equations new variables Yi and Xi will be introduced. In the 2D Peters’ model, inflow state variables are not a function of x. This feature is the main advantage of the Wagner model with respect to 2D Peters’ model when the Galerkin method is used to solve the equations. As a result of this feature, for an aeroelastic analysis of a high aspect ratio wing using Wagner model, the division of the wing into a number of smaller spanwise sections called aerodynamic elements is not needed. So the number of degree of freedom using Wagner model is less than the number of degree of freedom using 2D Peters’ model. Accordingly, the number of inflow state variables in Wagner model depends on the number of mode shapes used in each direction as follows T = 2P + 2Q

(14)

where P and Q are the number of mode shapes used in pitching and bending directions respectively. For the simplest case when M and N are one which is enough for linear analyses, T = 4. 2.2.2. Peters’ model In Peters’ aerodynamic model, the lift and moment , q2 and q4, are as follows (Peters et al., 2007) q2 = L = ρ∞ π b2 −¨v + U∞ φ˙ − baφ¨

[

]

[ ( ] ) 1 + 2π bρ∞ U∞ −˙v + U∞ φ + b − a φ˙ − λ0

(15)

2

q4 = Me.c =

[ ( ) ( ) ] 1 1 ρ∞ π b3 −av¨ + a − U∞ φ˙ − b a2 + φ¨ 2 8 [ ( ) ] 1 + 2π b2 ρ∞ U∞ −˙v + U∞ φ + b − a φ˙ − λ0

(16)

2

where ρ∞ , U∞ , b are same as before. λ0 can be derived using 2D or 3D inflow models. In 2D inflow model, λ0 is defined as follows ∞

1∑ λ0 ∼ bn λ n = 2

(17)

n

bn = (−1)n−1

(N + n − 1)! 1 1≤n≤N −1 (N − n − 1)! (n!)2

bN = (−1)N +1

(18) (19)

For an airfoil, 2D Peters’ state variables, λi s, are governed by differential equations as follows: [A] λ˙ +

{ }

U b

( ( ) ) 1 {λ} = {C } h¨ + U∞ θ˙ + b − a θ¨ 2

(20)

where θ and h describe airfoil motion and λs are state variables which are only function of time. For more details see Peters et al. (2007) and Peters et al. (1995).

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

255

Table 1 Number of total inflow states with respect to M and N. Highest power of x

Harmonic number, m

0 1 2 3 4 5 6 7 8

0

1

1 1 2 2 3 3 4 4 5

1 1 2 2 3 3 4 4

Number of total inflow states

2

3

1 1 2 2 3 3 4

1 1 2 2 3 3

4

1 1 2 2 3

5

1 1 2 2

6

1 1 2

7

1 1

8

1

1 3 6 10 15 21 28 36 45

When there is a wing instead of an airfoil, θ and h are functions of x, where x varies from 0 to L along the wing or beam. Thus we divide the whole beam into the number of spanwise elements and for each element Eq. (20) is used in which θ and h describe the middle of each element’s motion. In the simplest case when the whole beam is considered as one element, θ and h describe the motion of x = 2L . When there are two elements, θ and h describe the motion of x1 = 4L , x2 = 3L . 4 In 3D inflow model, λ0 is defined as follows

λ0 =

∞ ∑

φjr (x) J0

(

rb L

r ,j

)

] [ r αj cos (r ψ) + βjr sin (r ψ)

(21)

The radial expansion functions φ˜ jr (x) have the following form which are related to Legendre function. In addition, r ,j

... =

∑M ∑N

r =0 j=r +1,r +3,... . . . in which M is the total number of harmonics and N determine the number of radial function for each harmonic (N − 1 = highest pow er of x ). Consequently, r is the harmonic number and j is the polynomial number. ψ is the azimuth of the blade which is π2 . The choice of inflow states is based on the relationships in Table 1. The table shows the number of radial shape functions for each harmonic (m) in order to have radial terms up to a given power of x .

∑

r

φjr (x) =

P j (ν)

(22)

ν

r

√

where P j (ν) are the associated Legendre function of the first kind and ν = 1 − x and x = x/L. The αjr s and βjr s are inflow state variables which should be calculated from a set of ordinary differential equations below:

{ }∗

[M] αjr

{ }∗

[M] βjr

[ ]−1 { r } 1 mc + Lc αj = {τn } [ ]−1 { r } + Ls βj =

2

(23)

2 1

{τ ms } (24) 2 n mc ms τn and τn are the pressure coefficients which serve as the inflow generalized forces, or the inflow forcing functions. [M] , [Lc ] , and [Ls ] represent coefficient matrices that relate the pressure coefficients to the inflow state variables. The pressure coefficients offer the interface between the wake and lift models and are defined as follows τ

me n

≡

1( 2

τ

mc n

+ iτ

ms n

)

=

1 2π

1

(∫

)

fm Lφ (x) dx e−imψ m n

(25)

0

Lq is dimensionless circulatory lift on ρ∞ V2 R. The functions fm which appear in the τnm equations, can be found for any given chordwise pressure distribution. For more details the reader is referred to Peters and He (1995) and Peters and He (1991). As was noted, for aeroelastic analysis of a high aspect ratio wing, when using 2D Peters’ model, the beam needs to be divided into a number of spanwise aerodynamic elements, but using Wagner model the division of the beam into aerodynamic elements is not needed. So the number of degree of freedom using Wagner model is less than the number of degree of freedom using Peters’ model. Although in 3D Peters’ model which combines the lift and moment in 2D Peters’ model with a 3D wake model, the sectional division of the wing is not needed, as it was shown in Table 1, the number of total inflow states and consequently the number of degree of freedom is much higher than the number of degree of freedom using Wagner model. However, 3D Peters’ model is much more suitable for aeroelastic analysis of low aspect ratio wings rather than 2D models (Peters and He, 1995, 1991). In addition, 2D and 3D Peters’ models can be used when the bending deflection along the chord of the wing is considered (Peters et al., 2007). Also, for aeroelastic analysis of a rotating beam, the 3D Peters’ model is more accurate as well. 2.3. Aeroelastic model Substituting q2 and q4 from the desired aerodynamic model (Wagner or 2D Peters’ models) into Eqs. (1) and (2) respectively, nonlinear aeroelastic equations will be derived. In the next step to determine the boundary of dynamic

256

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

instability or flutter, a linear aeroelastic analysis has been studied. For this the nonlinear terms in the governing equations have been omitted. Next using Galerkin approach the equations are written in matrix form as follows [M] {x¨ } + [C ] {x˙ } + [K ] {x} = 0

(26)

where [M ], [C ] and [K ] are aeroelastic inertia matrix, aeroelastic damping matrix, and aeroelastic stiffness matrix respectively and {x} are generalized coordinates of the system which are a function of time. Then using eigenvalue method, the boundary of dynamic instability has been determined. 3. Results The results of the study are presented in four main parts: 3.1. Linear aeroelastic (flutter) analysis of an airfoil In this part, linear aeroelastic analysis of a 2D airfoil using Wagner and 2D Peters’ models has been done. The flutter speeds computed by these two aerodynamic models have been compared in details. 3.2. Unsteady lift results for a wing with prescribed simple harmonic motion In order to study the accuracy of the 2D aerodynamic theories in predicting the unsteady lift on a high aspect ratio wing a comparison between 2D and 3D aerodynamic models on low and high aspect ratio wings has been done. This comparison is based on a dimensionless analysis. In this case,(first) a cantilever wing is given a prescribed simple harmonic elastic pitching motion about its quarter-chord. This motion θ x, t , is specified through a single radial shape function. The pitching motion and resulting lift are assumed to be simple harmonic

( ) θ x, t = θˆ Θ (x) eiωt

(27)

where ω = ωVR and ˆ denotes a complex quantity (magnitude and phase). Also, x is the spanwise coordinate centered at root of the wing, dimensionless on L. In next step, the heaving motion is added. Similarly, heaving motion is specified through a single radial shape function and assumed to be simple harmonic

ˆ (x) eiωt w x, t = hW

(

)

(28)

In this case, the structural equations are not needed, and the distribution of the lift and inflow have been presented. 3.3. Linear aeroelastic (flutter) analysis of a high aspect ratio wing In this part, the flutter analysis using 2D and 3D aerodynamic models has been done. Considering linear governing equations and using the Galerkin approach, ordinary differential equations have been derived and using eigenvalue solutions, the flutter speed has been determined. 3.4. Nonlinear aeroelastic analysis of a high aspect ratio wing This part includes the nonlinear aeroelastic analysis of a high aspect ratio wing using 2D aerodynamic models considering all nonlinear terms up to third order. Similar to the previous part, using Galerkin approach, the nonlinear governing equations have been solved and the time response of the system above the speed of flutter has been presented. To provide nonlinear equations for Galerkin method, v and ϕ are substituted with their approximate modal expansions. In this study, the first five natural modes in each motion are considered. As mentioned above, this study has been divided into four main parts. In the first part, the results of aeroelastic analysis of an airfoil using Wagner and 2D Peters’ aerodynamic models is shown. Then, a comparison between 2D and 3D aerodynamic models on the low and high aspect ratio wings with simple harmonic motion has been done. In the third part, the results of linear aeroelastic analysis for a high aspect ratio wing using these 2D aerodynamic models represented and finally, in the last part, the results of nonlinear aeroelastic analysis of a high aspect ratio wing are shown. In each part, the results obtained from using different aerodynamic models are compared with each other. Also, the convergence of the results is shown using different numbers of mode shapes, state variables, and elements. Tables 2 and 3 give the structural and planform data for the airfoil and wing models under investigation.

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

257

Table 2 Parameters of the 2D airfoil (NACA 0012). Chord Mass Moment of inertia (C.G) a xA a Compression spring stiffness (Kh ) Torsional spring stiffness (Kθ ) Altitude Density of air a

1.83 m (6 ft) 19.62 kg (1.3447 Slugs) 4.1021 kg − m2 (3.0256 Slugs − ft2 ) −0.2 0.1 200 kg/m (134.47 lb/ft) 26.1306 kg − m/rad (1891 ft − lb/rad) 20 km 0.0889 kg/m3

Distance between elastic axis and center of gravity (positive in flow direction). Table 3 Parameters of the wing. Half-span Chord Mass per unit length Moment of inertia (50% chord) Spanwise elastic axis Center of gravity Bending rigidity Torsional rigidity Altitude Density of air

16 m 1m 0.75 kg/m 0.1 kg m 50% chord 50% chord 2 × 104 N m2 1 × 104 N m2 20 km 0.0889 kg/m3

Table 4 Speed of flutter using 2D Peters’ model with the different number of state variables (M). M U (m/s) Different with Wagner model (%)

2 45.96 8.20

3 52.42 −4.69

4 50.58 −1.02

5 50.34 −0.54

6 49.68 0.78

7 49.89 0.34

8 50.04 0.06

9 50.19 −0.24

10 50.02 0.1

11 50.55 −0.96

12 48.90 2.34

3.5. Linear aeroelastic (flutter) analysis of an airfoil In this part, the results of an aeroelastic analysis of an airfoil using Peters’ and Wagner models have been presented. At first, the flutter speed was determined using Peters’ model with the different number of state variables (M). Table 4 shows the speed of flutter for M = 2 to M = 12. Fig. 2 shows the variation of the speed of flutter respect to the number of state variables in the aerodynamic model. As can be seen in Fig. 2 when five to ten state variables are used in Peters’ model, there is a little change in the flutter speed. As the number of state variables is increased to more than ten, there is a numerical divergence. So, it is concluded that using five to seven state variable in Peters’ model is enough and using more than 10 state variables should be avoided. Second, the flutter speed using Wagner model is determined to be 50.07 m/s. Fig. 3 shows damping with respect to velocity using Peters’ model with the different number of state variables from five to nine and compares with the result obtained from using the Wagner model. As can be seen, there is a good agreement between the results. 3.6. Unsteady lift results for a wing with simple harmonic motion The unsteady lift results on a cantilever fixed wing oscillating with simple harmonic pitching motion are presented first. For comparison with the results of Nibbelink (1992), the angle of attack has a linear distribution, θ (x) = x and the wing has an aspect ratio of 6. Also, it is suggested in Nibbelink (1992) and NIibbelink and Peters (1993) that for 3D Peters’ model it is enough to consider M = 8 and consequently, from Table 1, N = 9. A wing with a high aspect ratio of 32 is considered. Secondly, the unsteady lift results (on a cantilever fixed wing oscillating with simple harmonic pitching and heaving motion ) 3 2 4 6 1 with the distribution of θ (x) = 21 3x − x and W (x) = 11 (15x − 5x + x ) are provided here. In the second part, a wing with a high aspect ratio of 32 is only considered. Finally, the distribution of the lift for some aspect ratios between 6 and 32 are presented that shows above what aspect ratio the 3D effects of the aerodynamic model can be ignored. Fig. 4 shows the distribution of sectional lift’s real part and inflow using three different aerodynamic models. The pitching motion is only considered and the wing has an aspect ratio (AR) of 6. k is reduced frequency. The result using 3D Peters’ model for the wing of AR = 6 has been reported by Nibbelink (1992). As can be seen, for low aspect ratio wings there is a significant difference between the distribution of the lift and inflow for 2D and 3D aerodynamic models. Fig. 5 shows the same results when AR = 32. From Fig. 5 we can see there is a good agreement between the distribution of the lift and inflow obtained from 2D and 3D aerodynamic models for a high aspect ratio wing. When both heaving and pitching motions are considered, Fig. 6 presents the distribution of sectional lift’s real part and inflow when AR = 32. As can be seen from Fig. 6, for the majority part of high aspect ratio wings the distribution of the lift

258

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 2. Variation of flutter speed with respect to the number of state variables.

Fig. 3. Damping as a function of freestream velocity.

and inflow for 2D and 3D aerodynamic models are nearly the same and there is only a small difference at the tip section of the wing which is caused by the tip loss effect modeled in the 3D aerodynamic model. As will be shown in the next part, this small difference in the distribution of the lift between 2D and 3D aerodynamic models for a high aspect ratio wing has a minor effect on the flutter speed. Fig. 7 shows the distribution of the lift for aspect ratios of 10, 15 and 20. As can be seen above the aspect ratio of 15, the 3D effects of the aerodynamic model can be ignored and a 2D aerodynamic model is sufficient to do an aeroelastic analysis. 3.7. Linear aeroelastic (flutter) analysis of a high aspect ratio wing In this part, a linear aeroelastic analysis of a beam, using Wagner, 2D and 3D Peters’ aerodynamic models is done and the results are compared. Fig. 8 shows damping with respect to velocity when Wagner aerodynamic model is used. As can be seen, the flutter speed is 32.99 m/s. This is close to the values determined by earlier investigators (Patil et al., 2001a; Tang and Dowell, 2004; Ovesy et al., 2012). Table 5 illustrates the speed of flutter based on 2D Peters’ model using a different number of state variables and elements. M is the number of state variable and Z is the number of spanwise aerodynamic elements. As can be seen, there is excellent convergence as M and Z increase. It can be seen from Table 5, that by considering five spanwise elements and using five state variables for each element, there is a good agreement with the results obtained from Wagner model. So it can be concluded that in order to do the aeroelastic analysis of a high aspect ratio wing, in this case, using 5 elements along the beam is enough to reach convergence. Fig. 9 presents the variation in flutter speed using a different number of state variables when five spanwise elements are considered. It is apparent from Figs. 2 and 9 that the variation in flutter speed for a 2D airfoil and a 3D wing is very similar

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

259

Fig. 4. The distribution of sectional lift’s real part and inflow (pitching motion, AR = 6). Table 5 Flutter speed using different number of state variables and elements. M M M M M M M

=2 =3 =4 =5 =6 =7 =8

Z =1

Z =2

Z =3

Z =4

Z =5

Z =6

30.40 34.91 30.95 30.66 31.12 31.38 31.46

30.86 36.32 33.12 32.69 32.73 33.00 33.17

30.92 36.45 33.41 32.98 32.95 33.22 33.39

30.95 36.49 33.50 33.07 33.02 33.29 33.46

30.96 36.51 33.55 33.12 33.06 33.32 33.49

30.96 36.52 33.57 33.14 33.07 33.34 33.51

to each other. In addition, as was noted above, considering five to seven state variables for each spanwise element is enough to do the aeroelastic analysis of a high aspect ratio wing. Table 6 compares the flutter speed derived using 2D and 3D aerodynamic models in the current study and those reported by previous studies. It can be seen from Table 6 that there is a good agreement between the results obtained using 2D and 3D aerodynamic models for a high aspect ratio wing. Table 7 compares the number of sufficient aerodynamic states used in each aerodynamic model for doing flutter analysis. Using the Wagner aerodynamic model requires the smallest number of aerodynamic states and thus can significantly decrease the computational time.

260

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 5. The distribution of sectional lift’s real part and inflow (pitching motion, AR = 32). Table 6 Speed of flutter using different number of state variables and elements. Aerodynamic model

Previous studies

Wagner 2D Peters (M = 8, Z = 5) 3D Peters (M = 8, N = 9) Patil et al. (2001a) Tang and Dowell (2004) Ovesy et al. (2012) Speed of flutter (m/s) 32.99

33.49

32.78

32.21

32.5

33.8

Table 7 Number of sufficient aerodynamic states used in each aerodynamic model for doing flutter analysis. Aerodynamic model

Number of aerodynamic states

Wagner (P = 1, Q = 1)

2D Peters (M = 8, Z = 5)

3D Peters (M = 8, N = 9)

4

40

45

3.8. Nonlinear aeroelastic analysis of a high aspect ratio wing As it was noted before, by considering all nonlinear terms up to third order, a nonlinear aeroelastic analysis has been done using a Galerkin method a time-marching response solution. Different numbers of mode shapes are used to show the convergence of the results. To show how many mode shapes in Galerkin method is needed to reach convergence in final results, up to five mode shapes are taken in each motion variable and the results obtained from using the different number of mode shapes are compared with each other. It should be noted that the initial angle of attack is zero and Wagner aerodynamic model is only

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

261

Fig. 6. The distribution of sectional lift’s real part and inflow (pitching and heaving motions, AR = 32).

used to show the convergence of the results obtained from using the different number of mode shapes. Fig. 10(a) shows time response of the system above flutter speed in detail when five mode shapes are taken into account in each direction. As can be observed in Fig. 10(a), taking first three modes into account is sufficient. Fig. 10(b) and (c) illustrate the amplitude of the mode shapes when four mode shapes are taken into account in each direction. Ei is used to show the amplitude of ith mode shape. From Fig. 10, it can be seen that the first and the second modes are dominant in response and the effect of the third mode is weak and the fourth mode is almost negligible. Also, as shown in Fig. 10(c), there is period doubling and this is the most common of several routes to chaos for a nonlinear dynamical system. However, the nonlinear effect is not strong enough to make the overall response chaotic. Fig. 11 compares the time response of the system above the flutter speed obtained from Wagner and 2D Peters’ aerodynamic models when three first modes considered in each direction. There is a small difference in frequencies which is caused by the small differences in the stiffness and damping matrix obtained from these two aerodynamic models. In spite of small changes in frequencies, there is almost no difference in the response amplitude obtained from the two different aerodynamic models. Fig. 12 shows the time response of the system at different freestream velocities comparing results obtained from using the different number of mode shapes. As can be seen, with increasing speed, the difference between the results is increased. So it can be concluded that as the speed increases above the flutter speed, more mode shapes are needed to capture the response of the system. Also, it is observed that system’s response beyond the flutter speed goes into a limit cycle oscillations which becomes increasingly complex with increasing speed. As the speed increases above flutter speed, the oscillations loses periodicity and become chaotic. Fig. 13 shows Fast Fourier Transform (FFT) of the response at different freestream velocities. As can be seen from Fig. 13, as the speed increases, the frequency band of the response becomes wider. So, as the speed increases, the oscillations become chaotic and more mode shapes are needed to capture the response of the system with a broader range of frequency. A similar observation has been also reported by Patil et al. (2001a). To show that considering five mode shapes are enough for a good approximation of the response at the speed of 40 m/s, the amplitude of the mode shapes

262

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 7. The distribution of sectional lift’s real part (pitching motion, different aspect ratios).

are presented in Fig. 14 when five mode shapes are taken into account in each direction. It should be noted that when speed increases above 40.5 m/s, such as 41.5 m/s, the response becomes sensitive to the parameters of the solution method such as time step and tolerance, while for speeds below 40.5 m/s the response is not sensitive to the parameters of the solution method. A possible explanation for this sensitivity may be the lack of accuracy of the response as the speed increase above 40

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

263

Fig. 8. Damping with respect to velocity for beam (Wagner aerodynamic model).

Fig. 9. Variation of flutter speed with respect to the number of state variables with 5 elements.

m/s when only five mode shapes are taken in each direction. So it may be that for speeds above 40 m/s using more than five mode shapes is necessary to determine the response. This needs further study. The root mean square (RMS) of the response with respect to freestream velocity is shown in Fig. 15. It can be seen from Fig. 15, as speed increases, RMS increases and after the speed of 39 m/s there is a sudden increase in RMS and then there is no change in RMS. To show the effects of nonlinear inertial and damping terms, the response of the system has been obtained with and without these terms at different freestream velocities and the results are compared in Fig. 16. As can be seen, the effects of nonlinear inertial and damping terms can be ignored and nonlinear stiffness terms are predominant near the speed of flutter, but as speed increases above the speed of flutter, nonlinear inertial and damping terms have considerable effect on the response. Figs. 17 and 18 show time history response of tip displacement and tip twist of the wing using first five modes at speeds of 35 m/s and 40 m/s respectively. As can be seen in Fig. 17 the amplitude of tip twist at the speed of 35 m/s is less than 8 degree which shows that using this aerodynamic model is suitable for this problem. Also as shown in Fig. 18 the amplitude of tip twist in the speed of 40 m/s is less than 15 degrees and nonlinear aerodynamic models which include flow separation is still not needed. But it should be noted that in this study the steady angle of attack was zero. So when there is a steady angle of attack which causes an increase in tip twist angle, a nonlinear aerodynamic model is needed to consider the effects of dynamic stall. 4. Conclusions In this paper, analysis and computational results of a linear and nonlinear aeroelastic model of a high aspect ratio wing subjected to subsonic airflow have been presented. A third-order nonlinear beam model has been used as a structural model. This model includes nonlinear inertial and damping terms as well as stiffness terms. To model aerodynamic loads, different aerodynamic models, Wagner and 2D and 3D Peters’ models, have been used. The comparison between these aerodynamic models is useful for determining a suitable dynamic model with both accuracy and computational efficiency which can be used to design a control system. Using different numbers of structural mode shapes, the convergence of the results has been

264

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 10. Time response above flutter speed (35 m/s) using different number of mode shapes and the amplitude of each mode shapes when 4 mode shapes are considered.

shown. There is a good agreement between results obtained using different aerodynamic models for high aspect ratio wings. A summary of the results is as follows:

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

265

Fig. 11. Time response above flutter speed (35 m/s) using three mode shapes.

• In 2D Peters’ model, using five to seven state variables is enough and there is a good agreement between the results obtained using the 2D Peters’ and Wagner models.

• For aeroelastic analysis of a high aspect ratio wing using 2D Peters’ model, in this case, considering 5 spanwise elements is enough to reach convergence. In addition, considering five to seven state variables for each element is sufficient.

• Comparing Wagner and Peters’ models which both are state-space aerodynamic models with the Theodorsen (frequency domain) or Wagner(time domain) model, these two models do not include the Theodorsen function, C (k),

• • • • • •

where k is the reduced frequency. So, it is easier to use the state space or time domain models in forming an eigenvalue problem using the k or p–k method. Also, they are more convenient to derive a dynamic model which can be used in the design of control systems. Of course, the Wagner and Theodorsen models are fundamentally equivalent being essential Fourier Transforms of each other. For wings with the aspect ratio of 15 or more, a 2D aerodynamic model is accurate enough to conduct an aeroelastic analysis. For aeroelastic analysis of a high aspect ratio wing with constant chord using state-space Wagner model decreases the number of degrees of freedom compared to the 2D and 3D Peters’ models. This can significantly decrease the computational time. The effects of nonlinear inertial and damping terms can be ignored near the flutter speed, but as the airspeed increases above the flutter speed, they have considerable effect on the limit cycle oscillation response. As the speed increases, more structural shape functions are needed to capture the limit cycle oscillation response of the wing. It is observed that the system’s response above the flutter speed which goes into limit cycle oscillations becomes increasingly complex with increasing speed. As the speed increases the oscillations lose periodicity and become chaotic. For the cases studied, the amplitude of the tip twist above the flutter speed is less than 15 degrees and nonlinear aerodynamic models which include flow separation is not needed. However, when there is a steady angle of attack, the

266

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 12. Time response in different freestream velocities using different number of mode shapes.

amplitude of the tip twist increases and consequently a nonlinear aerodynamic model may be needed to consider the effects of dynamic stall. This could well lead to an increase in the drag force and one may need to consider chordwise

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 13. FFT of the response in different freestream velocities.

267

268

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 14. The amplitude of the mode shapes at a speed of 40 m/s.

Fig. 15. The root mean square (RMS) of the response with respect to freestream velocity.

motion in addition to out of plane motion. Clearly considering the effect of the initial angle of attack on post flutter behavior and chaotic motion of the wing is an issue of interest and is proposed for the future study.

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 16. Time response with and without nonlinear inertial and damping.

269

270

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. 17. Tip displacement and tip twist of the wing using first five modes at a speed of 35 m/s.

Appendix. Derivation of the structural model A schematic model of the wing is shown in Fig. A.19. Two coordinate systems have been introduced to describe the undeformed and deformed geometries of an initially straight beam where the XYZ axes are in Cartesian coordinate system describing the undeformed geometry and the ξ ηζ axes are orthogonal curvilinear coordinate system describing the deformed geometry. Moreover, u, v, and w denote the displacements of the reference point on the representative crosssection (as shown in Fig. A.19) along the X , Y , and Z directions, respectively; s denotes the undeformed length from the root of the beam to the reference point. The effects of gravity have been neglected in the structural model. Furthermore, X and ξ represent the neutral axis of the beam before and after deformation, respectively; ix , iy and iz denote the unit vectors along the X , Y , and Z axes; and i1 , i2 , and i3 denote the unit vectors along ξ , η and ζ axes. v denotes the displacement of the reference point on each cross-section along y and φ specifies rotation about final position of ξ -axis. We use three consecutive Euler angles ψ, θ, and φ to describe the rotation from the undeformed position to the deformed one. These angles are executed in this order that first, the XYZ system is rotated by an angle ψ about the Z axis to an intermediate coordinate system X1 Y1 Z . Second, the X1 Y1 Z system is rotated by an angle θ about the Y1 axis to an intermediate coordinate system ξ Y1 Z2 . Finally, the ξ Y1 Z2 system is rotated by an angle φ about the ξ axis to an intermediate coordinate system ξ ηζ . In this case, w and θ are zero, as a result of ignoring chord-wise motion. The unit vectors of two coordinate systems XYZ and ξ ηζ are related by: i1 i2 i3

{ }

ix

{ } = [T ] iy iz

(A.1)

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

271

Fig. 18. Tip displacement and tip twist of the wing using first five modes at a speed of 40 m/s.

where

1 0 0

0 cos φ − sin φ

[ [T ] =

0 sin φ cos φ

][

cos ψ − sin ψ 1

sin ψ cos ψ 0

0 0 1

] (A.2)

is a rotational transform matrix. ψ and φ are angles of rotation about Z and ξ axes respectively. In this paper, only a summary of the derivation of equations will be presented and for more detail, readers are referred to Nayfeh (2004). In order to derive governing equations, we use Hamilton’s principle which states

∫

t

(δ T − δπ + δ Wnc )dt = 0

(A.3)

t0

where T is the kinetic energy, π is the elastic energy, and δ Wnc is the variation of the neoconservative energy due to external loads and damping. It can be shown that

δT =

b

∫ a

∫

ρ R¨ .δ RdAds a

(A.4)

272

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

Fig. A.19. Schematic diagram of the wing model showing local and global coordinate.

and

δπ =

∫ L∫ 0

(σ11 δϵ11 + σ12 δϵ12 + σ13 δϵ13 + σ22 δϵ22 + σ23 δϵ23 + σ11 δϵ33 )dAds

(A.5)

It follows from Fig. A.19 that the deformed position vector R of an arbitrary point on the representative cross-section is given by: R = (s + u) ix + v iy + yi2 + zi3

(A.6)

As it was mentioned before, to fully account for geometric nonlinearities, local stress and strain measures have been used. Hence, the following local displacement field can be introduced as: U = u1 i1 + u2 i2 + u3 i3

(A.7)

where u1 = u01 − yθ3 + z θ2 ,

(A.8)

u2 = u02 − z θ1

(A.9)

u3 = u03 + yθ1

(A.10)

In these equations, θ1 , θ2 and θ3 are the rotation angles of each cross-section along the beam with respect to the ξ , η and ζ axes, respectively. ξ ηζ as defined earlier is the local Cartesian coordinate system associated with curvilinear coordinates and attached to each cross-section along the beam. So, the rotation angles of the cross-section are infinitesimal with respect to ξ ηζ . Similarly, ui s are infinitesimal with respect to ξ ηζ . So we have u01 = u02 = u03 = 0,

θ1 = θ2 = θ3 = 0

(A.11)

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

∂ u0 ∂ u02 = 3 = 0, ∂s ∂s

∂ u01 ≡e ∂s

∂θ2 ≡ ρ2 , ∂s

∂θ1 ≡ ρ1 , ∂s

273

(A.12)

∂θ3 ≡ ρ3 ∂s

(A.13)

It follows from Eqs. (A.7) and (A.13) that

∂U = (e − yρ3 + z ρ2 ) i1 − z ρ1 i2 + yρ1 i3 ∂s

(A.14)

∂U ∂U = =0 ∂y ∂z

(A.15)

Consequently, the local strains are given by

ϵ11 =

∂U .i1 = e − yρ3 + z ρ2 ∂s

(A.16)

ϵ12 =

∂U ∂U .i2 + .i1 = −z ρ1 ∂s ∂y

(A.17)

∂U ∂U .i3 + .i1 = yρ1 ∂s ∂z = ϵ33 = ϵ23 = 0

ϵ13 =

(A.18)

ϵ22

(A.19)

The nonconservative work term can be defined as:

δ Wnc =

L

∫

(q2 δv + q4 δθ1 ) ds

(A.20)

0

where qi are distributed external loads. After some calculations and expanding all the nonlinear transcendental terms into polynomials about the static equilibrium position and keeping the nonlinear terms up to third order, the governing equations of vibrations can be written as follows (for more detail see chapter 4 of Nayfeh, 2004) mv¨ + D33 v ′′′′ = −mv ′′

− mv

′

s

[∫ s ∫ ∫

( 2 ) v˙ + v ′ v¨ ′ dsds 0

L s

( ′2 ) v˙ + v ′ v¨ ′ ds

[0 − D33

]

4v ′ v ′′ v ′′′ + v ′2 v ′′′′ + v ′′3 − 2φφ ′′ v ′′ −2φ ′2 v ′′ − 4φφ ′ v ′′ − φ 2 v ′′′′

]

(A.21)

[ ] − D22 2φ ′2 v ′′ + 2φφ ′′ v ′′ + 4φφ ′ v ′′ + φ 2 v ′′′′ + q2 + j3 v¨ ′′ j1 φ¨ − D11 φ ′′ = (D33 − D22 ) v ′2 φ + q4

(

)

(A.22)

Also since the beam is inextensional we have 1 u′ = − v ′ 2 (A.23) 2 In Eqs. (A.21) to (A.23), q2 and q4 are distributed external loads, u is axial deformation and m, j1 , j3 , D11 , D22 and D33 are defined as

∫

ρ dA

m= A

(A.24)

274

A.H. Modaress-Aval, F. Bakhtiari-Nejad, E.H. Dowell et al. / Journal of Fluids and Structures 85 (2019) 249–274

∫ j1 =

ρ (y2 + z 2 )dA

(A.25)

ρ y2 dA

(A.26)

A

∫ j3 = A

D11 = GJ , D22 = EI2 , D33 = EI3

(A.27)

where ρ , GJ, EI2 , and EI3 are density, torsional rigidity, bending rigidity in the chordwise direction and bending rigidity in the flapwise direction, respectively. References Arena, A., Lacarbonara, W., Marzocca, P., 2013. Nonlinear aeroelastic formulation and postflutter analysis of flexible high-aspect-ratio wings. J. Aircraf. Bakhtiari-Nejad, Firooz, Aval, A.H.M., Dowell, EarlH, Shahverdi, Hossein, 2017. Linear and nonlinear aeroelastic analysis of a high aspect ratio wing subjected to incompressible air flow. In: ASME 2017 IDETC/CIE International Design Engineering Technical Conferences & Computers and Iformation in Engineering Conference. Castellani, M., Cooper, J.E., Lemmens, Y., 2017. Nonlinear static aeroelasticity of high-aspect-ratio-wing aircraft by finite element and multibody methods. J. Aircr. 54 (2), 548–560. Chang, C.-S., Hodges, D.H., Patil, M.J., 2008. Flight dynamics of highly flexible aircraft. J. Aircr. 45 (2), 538–545. Hallissy, B., Cesnik, C., 2011. High-fidelity aeroelastic analysis of very flexible aircraft. Harmin, M., Cooper, J., 2011. Aeroelastic behaviour of a wing including geometric nonlinearities. Aeronaut. J. 115 (1174), 767–777. Hesse, H., Palacios, R., 2014. Reduced-order aeroelastic models for dynamics of maneuvering flexible aircraft. AIAA J.. Hodges, D.H., 2003. Geometrically exact, intrinsic theory for dynamics of curved and twisted anisotropic beams. AIAA J. 41 (6), 1131–1137. Jian, Z., Jinwu, X., 2009. Nonlinear Aeroelastic Response of High-aspect-ratio Flexible Wings. Chin. J. Aeronaut. 22 (4), 355–363. Krüger, W., 2007. Multibody dynamics for the coupling of aeroelasticity and flight mechanics of highly flexible structures. In: Proceedings of IFASD. Mardanpour, P., Hodges, D.H., Neuhart, R., Graybeal, N., 2013. Engine placement effect on nonlinear trim and stability of flying wing aircraft. J. Aircr.. Mardanpour, P., Hodges, D.H., Rezvani, R., 2014. Nonlinear aeroelasticity of high-aspect-ratio wings excited by time-dependent thrust. Nonlinear Dynam. 75 (3), 475–500. Murua, J., Palacios, R., Graham, J.M.R., 2012. Applications of the unsteady vortex-lattice method in aircraft aeroelasticity and flight dynamics. Prog. Aerosp. Sci. 55, 46–72. Nayfeh, Ali H.P.P., 2004. Linear and Nonlinear Structural Mechanics. Wiley-VCH, Germany. Nibbelink, B., 1992. Finite-state inflow applied to aeroelastic flutter of fixed and rotating wings(Ph. D. Thesis). NIibbelink, B., Peters, D., 1993. flutter calculations for fixed and rotating wings with State-space inflow dynamic. In: 34th Structures, Structural Dynamics and Materials Conference. Ovesy, H.R., Nikou, A., Shahverdi, H., 2012. Flutter analysis of high-aspect-ratio wings based on a third-order nonlinear beam model. Proc. Inst. Mech. Eng. G 227 (7), 1090–1100. Patil, M.J., Hodges, D.H., 2004. On the importance of aerodynamic and structural geometrical nonlinearities in aeroelastic behavior of high-aspect-ratio wings. J. Fluids Struct. 19 (7), 905–915. Patil, M.J., Hodges, D.H., 2006. Flight Dynamics of Highly Flexible Flying Wings. J. Aircr. 43 (6), 1790–1799. Patil, M.J., Hodges, D.H., Cesnik, C.E., 1999. Characterizing the effects of geometrical nonlinearities on aeroelastic behavior of high-aspect ratio wings. In: NASA Conference Publication. NASA. Patil, M.J., Hodges, D.H., Cesnik, C.E., 2001a. Limit-cycle oscillations in high-aspect-ratio wings. J. Fluids Struct. 15 (1), 107–132. Patil, M.J., Hodges, D.H., Cesnik, C.E.S., 2001b. Nonlinear aeroelasticity and flight dynamics of high-altitude long-endurance aircraft. J. Aircr. 38 (1), 88–94. Peters, D.A., He, C.J., 1991. Correlation of measured induced velocities with a finite-state wake model. J. Amer. Helicopter Soc. 36 (3), 59–70. Peters, D.A., He, C.J., 1995. Finite state induced flow models. II - Three-dimensional rotor disk. J. Aircr. 32 (2), 323–333. Peters, D.A., Hsieh, M.C.A., Torrero, A., 2007. A state-space airloads theory for flexible airfoils. J. Amer. Helicopter Soc. 52 (4), 329–342. Peters, D.A., Karunamoorthy, S., Cao, W.M., 1995. Finite state induced flow models. Part I: Two-dimensional thin airfoil. J. Aircr. 32 (2), 313–322. Shams, S., Lahidjani, M.S., Haddadpour, H., 2008. Nonlinear aeroelastic response of slender wings based on Wagner function. Thin-Walled Struct. 46 (11), 1192–1203. Sotoudeh, Z., Hodges, D.H., Chang, C.-S., 2010. Validation studies for aeroelastic trim and stability analysis of highly flexible aircraft. J. Aircr. 47 (4), 1240– 1247. Tang, D., Dowell, E.H., 2001. Experimental and theoretical study on aeroelastic response of high-aspect-ratio wings. AIAA J. 39 (8), 1430–1441. Tang, D., Dowell, E.H., 2002. Experimental and Theoretical Study of Gust Response for High-Aspect-Ratio Wing. AIAA J. 40 (3), 419–429. Tang, D.M., Dowell, E.H., 2004. Effects of geometric structural nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings. J. Fluids Struct. 19 (3), 291–306. Wang, Y., Wynn, A., Palacios, R., 2016. Nonlinear Modal Aeroservoelastic Analysis Framework for Flexible Aircraft. AIAA J.. Xie, C., Liu, Y., Yang, C., 2012. Theoretic analysis and experiment on aeroelasticity of very flexible wing. Sci. China Technol. Sci. 1–12. Zhao, Z., Ren, G., 2011. Multibody dynamic approach of flight dynamics and nonlinear aeroelasticity of flexible aircraft. AIAA J. 49 (1), 41–54.