- Email: [email protected]
Nonlinear aeroelastic analysis of high-aspect-ratio wings in low subsonic flow
Acta Astronautica 70 (2012) 6–22
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro
Nonlinear aeroelastic analysis of high-aspect-ratio wings in low subsonic flow K. Eskandary, M. Dardel n, M.H. Pashaei, A.K. Moosavi Department of Mechanical Engineering, Babol University of Technology, P.O. Box 484, Shariati Avenue, Babol, Iran
a r t i c l e in f o
abstract
Article history: Received 9 December 2010 Received in revised form 13 June 2011 Accepted 11 July 2011 Available online 5 August 2011
In this study, aeroelastic characteristics of high-aspect-ratio wing models with structural nonlinearities in quasi-steady aerodynamics flows are investigated. The studied wing model is a cantilever wing with double bending and torsional vibrations and with large deflection ability in accordance with Hodges–Dowell wing model. This wing model is valid for long, straight and thin homogeneous isotropic beams. The aerodynamics model is based on quasi-steady aerodynamic which is valid for aerodynamic flows without wake, viscosity and compressibility effects. The effect of different parameters such as mass ratios and stiffness ratios on flutter and divergence velocities and limit cycle oscillation amplitudes are carefully studied. & 2011 Elsevier Ltd. All rights reserved.
Keywords: High-aspect-ratio wing Flap-lag and torsion Flutter Limit cycle
1. Introduction Helicopter blades and some models of aircrafts have long wings so that flexibility and nonlinear deformations have important roles in them. Structural flexibility of this type of structures is very high, and structures are suffering from large deflections and deformations. Large deflection of these structures is describable by structural nonlinear modeling. Presences of nonlinear structural factors with aerodynamic forces lead to undesirable aeroelastic phenomena. Therefore, careful review and understanding of these phenomena in aeroelastic systems is vital. Modeling of long wings has been done by several researchers which some of them have been mentioned here. Hodges and Dowell developed the equations of motion by two complementary methods, Hamilton’s principle and the Newtonian method [1]. Tang and Dowell constructed an experimental high-aspect-ratio wing aeroelastic model with a slender body at the tip, and measured the response due to flutter and limit-cycle oscillations
n
Corresponding author. E-mail address: [email protected] (M. Dardel).
0094-5765/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.actaastro.2011.07.017
(LCO) in a wind-tunnel test. They also developed a theoretical model and made calculation to correlate with the experimental data [2]. Tang and Dowell presented a nonlinear response analysis of a high-aspect-ratio wing aeroelastic model excited by gust loads theoretically and examined it in a wind-tunnel [3]. Patil, Hodges and Cesnik described a formulation for aeroelastic analysis of aircraft with high-aspect-ratio wings. The analysis was a combination of a geometrically exact mixed formulation for dynamics of moving beams and finite-state unsteady aerodynamics with the ability to model dynamic stall [4]. Malatkar explored the impact of kinematic structural nonlinearities on the dynamics of a highly deformable cantilevered wing. Two different theoretical formulations were presented and analyzed for nonlinear behavior [5]. Nichkawde presented a study of the nonlinear vibrations of metallic cantilever beams and plates subjected to transverse harmonic excitations. Both experimental and theoretical results were presented [6]. In this work, aeroelastic behavior of Hodges–Dowell wing model with incompressible quasi-steady aerodynamic model is investigated. A brief description of the Dowell–Hodges beam is given here. The theory is intended for application to long, straight, slender, homogeneous,
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Nomenclatures ah b Cv, Cw Cy c, c dD, dL dFv, dFw dMx dMxf E e G g h I1, I2 Ia Iy J L
dimensionless distance from wing section mid-chord to elastic axis wing section semi-chord section lag and vertical bending damping section twist damping wing chord and dimensionless chord, c/L section drag and lift forces section lag and vertical component forces section pitch moment about elastic axis section pitch moment about one-quarter chord modulus of elasticity section mass center from elastic axis shear modulus gravitational constant flap displacement vertical, lag area moments mass moment of inertia about elastic axis mass moment of inertia per unit length about center of mass torsional stiffness constant wing span
isotropic beams with large displacements and is accurate to second order based on the restriction that squares of bending slopes, twist, h/L, and c/L are small with respect to unity. Radial nonuniformities (mass, stiffness, twist, etc.), lag offsets of the mass centroid and tension axes from the elastic axis, and warping of the cross section are included. Other more specialized details are not considered, such as blade root feathering flexibility, torque offset, blade sweep, and curvature; nor are configurations considered in which the feathering bearing is replaced with a torsionally flexible strap [1]. The effect of wakes, compressibility and viscosity of incompressible quasi-steady aerodynamic model are ignored. In current study, determination of flutter and divergence velocities and clarifying of limit cycle amplitudes are investigated. Effect of changing mass ratio and distance between elastic axis and center of mass and other parameters are carefully examined in aeroelastic properties. This study is done for the different wing section models to evaluate the effect of mass and inertia on the limit cycle amplitude and flutter velocity. The effects of structural damping are considered in structural modeling.
m ra t U, Un UL Vi, Wi v w x xf xa
mass per unit length of the wing radius of gyration about elastic axis time freestream velocity and dimensionless velocity linear flutter velocity generalized coordinates for bending lag or edgewise bending deflection vertical or flapwise bending deflection position coordinate along wing span position of flexural axis (pivot) dimensionless distance from elastic axis to center of mass y pitch angle of wing section yi generalized coordinates for torsion y0 steady angle of attack at root section m mass ratio z, Z dimensionless flap displacement r air density t dimensionless time o frequency of limit cycle oscillations oz, oZ, oa natural frequencies in flap, lag and pitch on1 , on2 frequency ratio ð Þ0 dð Þ=dx ð: Þ dð Þ=dt
beams. In this study, the wing has no initial twist and initial warping. In current wing model, the elastic and mass axes are not coincided, and structural damping is included in equations of motions. The considered wing model is a cantilevered wing model with bending displacements of w and v and torsion of y, which are shown in Fig. 1. The flap deflection is denoted by w, positive downward direction, lag deflection is denoted by v, positive in stream-wise direction and the torsion angle y, positive nose up. The sketch of wing section is shown in Fig. 2, where c is the chord, b is the semi-chord length, ahb is the distance from wing section mid-chord to elastic axis and xab is the distance from elastic axis to the center of mass. The length of the wing is L. According to Hodges–Dowell equations for a uniform elastic wing and with ignoring the wing section warping,
2. Aeroelastic equations 2.1. Structural equations of motions In this section, a brief description of structural equations of motion of wing model based on Hodges–Dowell wing model is presented. A Hodges–Dowell wing model presents double bending and torsional vibrations. Hodges– Dowell wing model is second-order equations which are valid for long, straight and thin homogeneous isotropic
7
Fig. 1. Flexible deflected wing.
8
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
section. According to quasi-steady aerodynamic theory, aerodynamic moments and forces are as follows [8]: c € Lðx,tÞ ¼ rpb2 w xf y€ þ rpb2 U y_ 2 _! _ 3 y w 2 cxf þ rpU c y þ þ ð7Þ 4 U U c c _ w xf y_ Mxf ðx,tÞ ¼ rpb2 xf 2 2 rpb4 _ 3 y cxf rpb2 U y_ 4 8 _! _ 3 y 1 w 2 2 cxf þ rpU ec y þ þ rpUc3 y_ 4 16 U U ð8Þ
Fig. 2. Schematic figure of the wing section.
where the aerodynamics center of e and distance between elastic axis to the wing mid-chord are defined by
equations of motions will be as follows [7]: dF € _ ¼ w € EI1 wð4Þ þ ðEI2 EI1 Þðyv00 Þ00 þmwmx a by þ Cw w dx
ð1Þ
e¼
dMx dx
ð2Þ
In addition, by ignoring the drag force for this low speed, inviscid aerodynamics model, we can write
€ þ Cy y_ ¼ GJ y þ ðEI2 EI1 Þw00 v00 þmra2 y€ mxa bw 00
EI2 vð4Þ þðEI2 EI1 Þðyw00 Þ00 þmv€ þCv v_ ¼
dFv dx
xf 1 , c 4
xf ¼
b þ ah b 2
dD ¼ 0 ð3Þ
where m, Iy, EI1, EI2, GJ, Cw,Cy and Cv are the mass per unit length of the wing, mass moment of inertia per unit length about mass center, flap and lag bending stiffness, torsional stiffness, the flap, torsional and lag structural damping coefficients, respectively. e is the distance between mass and elastic axes and x is the position coordinate along the wing span. dFw/dx, dMx/dx and dFv/dx are the aeroedynamics force and moments, which will be described as follows. Important structural nonlinear terms have been included in Eqs. (1)–(3) and third- or higher-order nonlinear terms have been ignored. 2.2. Quasi-steady aerodynamics model According to the work presented in [2,3,7], aerodynamic forces and moment are as follows: dFw ¼ dL þ y0 dD
ð4Þ
dMx ¼ dMxf
ð5Þ
dFv ¼ dD þ y0 dL
ð6Þ
where y0 is a constant angle of attack (which is assumed to be zero), L(t) is the aerodynamic lift force, M(t) is the aerodynamic moment, D(t) is the aerodynamic drag force. In works done by Dowell et.al. aerodynamic force has been made in accordance with Onera aerodynamics model. In present work, quasi-steady aerodynamic model has been used. Also, in this aerodynamics model, the flow is assumed to be incompressible, inviscid and irrotational. Quasi-steady models assume that there are only four contributions to the aerodynamic forces, horizontal air velocity U, at angle y(t) to cross section, wing flap velocity, _ wðtÞ, normal component of torsion velocity, ðxf xÞy_ ðtÞ, and local velocity induced by the vorticity around the airfoil. xf is the elastic axis position of the wing cross
ð9Þ
ð10Þ
In current study, the angle of attack is assumed to be zero. 2.3. Discretizing aeroelastic equations using assumed mode method With substituting aerodynamics forces and moment according to Eqs. (7)–(9) in structural equations of motions of (1)–(3), the complete aeroelastic equations will be obtained. These equations can be descritized in accordance with Galerkin’s method. In this method, at first the w, v and y displacements are expressed in terms of to the assumed mode method. Then with substituting these displacements in equations of motions and multiplying each equation with an admissible mode function and integrating along the whole area of the wing, the equations of motions will be descritized. In assumed mode method, displacements are written based on the products of time dependent generalized coordinates and mode functions which satisfy the geometric boundary conditions. In accordance with the assumed mode method the bending and torsional displacements can be described as follows: X wðx,tÞ ¼ Wi ðtÞci ðxÞ ð11Þ
yðx,tÞ ¼ vðx,tÞ ¼
X X
yi ðtÞfi ðxÞ
ð12Þ
Vi ðtÞai ðxÞ
ð13Þ
where Wi(t), yi(t) and Vi(t) are the generalized coordinates and ci(x), fi(x) and ai(x) are the mode functions that satisfy the geometric boundary conditions. By combining structural Eqs. (1)–(3) with aerodynamics Eqs. (7)–(9) and multiplying both sides of the equation in a desired mode shape and integrating along the wing, the original partial differential equations are converted to the ordinary differential equations. For further simplification of the equations, reducing the amount of the necessary calculations and selecting the
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
appropriate mode shape, high-order derivative terms are converted to lower-order derivative terms by integrating by part method. The resulting equations are as follows: Z L X XX Wj EI1 ci cð4Þ Vj yk j dx þ ðEI2 EI1 Þ 0
Z
L
0
þm
X
þCw
ci a00j f00k dx þ2 € j W
X
¼ rpb2
Z Z
þ rpb2 U
X
Z
yj
L 0
Z
þCy
X
y_ j
y€ j Z
Z
¼ rpb2 xf xf
0
ci fj dx þ y_ j
Z
X
L
ci fj dx
0
Z
L
ci cj dx
0
ð14Þ
XX X
L
fi fj dxme
0
Z
_ j W
Wj Vk € j W
Z
Z
L 0
L
fi cj dx
0
fi fj dx
2 c X
y€ j
Z
Z
L 0
L
fi cj dx
fi fj dx
0
rpb4 X €
yj
8
X Z L 3 cxf rpb2 U y_ j fi fj dx 4 0 X Z L Z X _ j þ rpUec2 U yj fi fj dx þ W þ
X
Vj
X
y_ j
X 1 rpUc3 y_ j 16
EI2
3 cxf 4
Z
L 0
L
0
þm
Z
L
fi fj dx
0
X
00
V€ j
¼ y0 rpb2
þ rpb2 U X
_ j W
Z 0
00
X Z 0
L
Z
L
0
L
ai aj dx þ Cv
X
€ j W
y_ j
Z
Z
0
L 0
L
0
fi cj dx
ð15Þ
XX 000
Wj yk
0
ai cj fk dx þ X
V_ j
Z 0
Z 0
L
ð4Þ
ai cj fk dx
L
The aeroelastic Eqs. (14)–(16) are written in non-dimensional form by defining the following parameters: rffiffiffiffiffiffiffiffiffi U t w v Ia m t ¼ 1 , z ¼ , Z ¼ , ra ¼ , , m¼ b b b mb2 rpb2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u R u GJ 1 ðf 00 Þ2 dx U , oa ¼ t 2 R01 1 2 , Un ¼ boa Iy L 0 ðf1 Þ dx vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u R 1 00 2 u EI ðc Þ dx oz ¼ t 14 R01 1 2 mL 0 ðc1 Þ dx vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u R 1 00 2 u EI ða Þ dx o o x oZ ¼ t 24 R01 1 2 , on1 ¼ z , on2 ¼ Z , x ¼ L oa oa mL ða1 Þ dx ð19Þ where z and Z are the dimensionless flap and lag deflections, oz, oZ and oa are the first natural frequencies of uncoupled flap, lag and torsional modes, U*, oni , m and ra are the dimensionless velocity, bending to torsion frequency ratios, mass ratio and the radius of gyration about the elastic axis, respectively, and x is the dimensionless coordinate along the span of the wing. These dimensionless parameters are defined in order to select the necessary parameters in accordance with the values given in [9], to validate the presented model in accordance the flutter velocities in this reference. With substituting these dimensionless parameters in Eqs. (14)–(16) and multiplying these equations in ci, fi and ai, respectively, the following equations of motion will be obtained:
on12
ai aj dx
Z L c X € ai cj dx xf yj ai fj dx 2 0
X Z ai fj dx þ rpUc U yj
ð18Þ
fi fj dx
fi fj dx
0
an ¼ ðcos bn cosh bn Þ=ðsin bn sinh bn Þ
0
L
L
0
Z
L
ai cj fk dx þ2
þ
0
Z
ai að4Þ j dx þðEI2 EI1 Þ
Z
fi c00j a00k dx
L
€ j W
bn ¼ 1:8751, 4:6941, 7:8547, 10:9955, 10:9955þ p,
And for torsional modes: 2n1 x fn ¼ sin p 2 L
ci fj dxrpUc
c X
2
ci fj dx
L
L
0
L
0
Z L c X € ci cj dx xf yj ci fj dx 2 0
L 0
X
X
y€ j
0
Appropriate mode shape for description of the displacements of w, v, and y is obtained from beams and bar with equivalent boundary conditions to the wing model. For descriptions of w and v displacements, Euler–Bernoulli’s cantilevered beam mode functions are used and for description of torsional displacement of y, mode functions of fixed– free bar are used. These mode functions are as follows: For bending modes: x x x x ð17Þ cn ¼ sin bn sinh bn þ an cos bn cosh bn L L L L That
fi f00j dx þðEI2 EI1 Þ
2 þmKm
Z
ci að4Þ fk dx j
ci cj dx
0
3 cxf þ 4 GJ
X
y_ j
X Z U yj
0
L
L
€ j W
X
Z
ci a000j f0k dx þ
ci cj dxme
0
X
L
L 0
_ j W
Z
9
U
R1 0
n2 R 1
N c21 dx X
002 0 c1 dx j ¼ 1
m
ai fj dx
X Z L 3 cxf ai cj dx þ y_ j ai fj dx 4 0 ð16Þ
Z
1 0
N Z a X xa h y€ j
L 0
zj
Z
0
¼
j¼1
1
ci cj dx
c00i c00j dx þ 1 1
0
ci fj dx þ
N ð22ah Þ X y_
m
j¼1
Z
1
X N
m
0
Z
j¼1
1
j
z€ j
1 0
ci cj dx
N Cw 2 X z_ n moa U m j¼1 j
ci fj dx
n2 R 1 2 N X N o1 0 c1 dx X EI2 1 Zy EI1 U n2 R 1 c’’2 dx j ¼ 1 k ¼ 1 j k 0
1
N 2X
mj¼1
yj
Z 0
1
ci fj dx
10
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Z
1
0
R1
2
0
1
!
N X
1 0 f1 dx 4e R U n2 1 f02 dx ra2 m N X
Z
y€ j
þ
!
1
c0i a000j fk dx
0
yj
Z
1
0
j¼1
1
fi fj dx
0
j¼1
Z
c0i a00j f0k dxþ
ð21Þ ða2 þ ð1=8ÞÞ 1þ h 2 ra m
f0i f0j dx þ
ah xa þ ra2 m ra2
X N
z€ j
Z
1
fi cj dx
0
j¼1
!
N Z 1 Cy ðð1=2Þ þ ð4e1Þðð1=2Þah ÞÞ X y_ j fi fj dx Iy oa U n ra2 m 0 j¼1
Z 1 N 4e X EI ¼ 2 z_ j fi cj dx 2 1 EI1 ra m j ¼ 1 0 R1
o1n2
r 2 U n2 a
o2n2 U n2
R1
0
N X N c21 dx X
0
c002 1 dx j ¼ 1 k ¼ 1
R1
0
N a21 dx X
0
a002 1 dx j ¼ 1
R1
þ
þ
Z
Zj
moa U n j ¼ 1
m
yj
j¼1
Z 0
Z 0
1
fi c00j a00k dx
N X
Z€ j
Z
j¼1
Z 0
1
0
ai aj dx Z
m j¼1
ai fj dx
1
ai fj dx
ð22Þ
1
N y0 X ai aj dx z€ j
1
0
j¼1
N 2y0 X
a00i a00j dx þ
Z_ j
Z
N y0 ah X y€ j
m
1 0
N X
Cv
z j Zk
0
j¼1
Z N 2y0 X z_
m
j
j¼1
Property
Value
Wing
Case 1
Case 2
Span (L) (m) Chord (c) (m) Dimensionless distance from wing section mid-chord to elastic axis (ah) Dimensionless distance from elastic axis to center of mass (xa) Flap structural modal damping (Cw) Lag structural modal damping (Cv) Torsional structural modal damping (Cy) Mass ratio (m) Radius of gyration about elastic axis (ra) Natural frequencies in flap bending (oz) (rad/s) Natural frequencies in lag bending (oZ) (rad/s) Natural frequencies in torsion (oa) (rad/s)
0.6299 0.1018 0.3
0.6299 0.1018 0.19
0.22
0.02
0.02 0.025 0.031 48 0.558 26 p
0.02 0.025 0.031 36.8 0.558 18 p
172.6 p 254.6 p 210 p 180 p
Flight conditions Density of air (r) (kg/m3)
1.1
1.1
1
ai cj dx
N y0 ð22ah Þ X y_ j
m
Table 1 Wing model data.
Z
1 0
ai fj dx
1 0
ai cj dx
R1 N X N EI1 on22 0 a21 dx X ¼ 1 zy R EI2 U n2 1 a002 dx j ¼ 1 k ¼ 1 j k 0 1 ! Z Z 1
0
00
1
0
a0i cj fk dx þ
0
000
a0i cj fk dx
ð23Þ
The interval of the integral is [0,1]. Then the whole aeroelastic Eqs.(21)–(23) can be written briefly as follows: € þ½Ceq fqg _ þ ½Keq fqg ¼ fFN g ½Meq fqg
ð24Þ
where [Meq], [Ceq], and [Keq] are mass, damping and stiffness matrices, respectively, and {FN} is a vector of structural nonlinearity. The vector {q} is defined as T fqg ¼ fzyZg . Now after presenting the whole aeroelastic equations, the validity and results of these equations will be examined.
Fig. 3. Damping of eigenvalues of the linear aeroelastic system vs. free stream velocity, case 1.
3. Aeroelastic results In this section, the aeroelastic results of the presented wing model for different cases will be studied. For this purpose, at first the wing structural and aerodynamics parameters will be selected. Here, two different wing section parameters are selected, which are shown in Table 1. These parameters are selected in compatibility with those parameters presented in [2,9]. These parameters are not the same to those given in these references, since in [2] the parameters are for a wing which represent stall phenomenon and parameters given in [9] are for a simple wing model with flap and pitch motions and incompressible, irrotational aerodynamics flow, which are different with those presented in current work. Hence
combinations of the parameters given by these references are given in Table 1. 3.1. Eigenvalue solution of the linear aeroelastic system At first eigenvalue solution of the linear aeroelastic system which determine the stability of the linear part of the system, will be studied. Changes of eigenvalues of the linear aeroelastic equation with free stream velocity for the case 1 of Table 1, have been shown in Figs. 3–5. Damping and frequency of eigenvalues are given by the real and imaginary parts of eigenvalues l of the linear system. Different combinations of mode functions are
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 4. Frequency vs. damping of the eigenvalues of the linear aeroelastic system, case 1.
11
Fig. 6. Damping of eigenvalues of the linear aeroelastic system vs. free stream velocity, case 2.
Fig. 5. Frequency of eigenvalues of the linear aeroelastic system vs. free stream velocity, case 1.
Fig. 7. Frequency vs. damping of the eigenvalues of the linear aeroelastic system, case 2.
selected, and the convergences of the eigenvalue solution are studied. For this figures, the number of mode functions for flap, lag and torsion displacements are 4, 2 and 2, respectively. With increase in the number of mode functions above these mentioned mode numbers, minor changes will occur in the predictions. A sample convergence of the solution based on the number of the mode functions will be presented in the next section. As seen in Fig. 3 with increase in the free stream velocity, one of the branches of the real part of the eigenvalues of the linear system will change its sign and become positive. In this case, the initially stable system becomes unstable. In flutter velocity a branch of
eigenvalues has intersection with the real axis. For this intersection, the real part of corresponded eigenvalue is zero. For the case 1 of Table 1 the flutter velocity is UF ¼102.271 m/s. From Fig. 4 the frequency of the flutter velocity is about 83:75 Hz. From Fig. 4 it is clear that two branches of eigenvalues are approaching toward each other gradually and near the flutter velocity they have the less distance from each other that show the bending–torsion flutter mode of the presented wing model. Between these branches, an internal resonance occurs. Also from Fig. 4 with increase in the free stream velocity, at first all of the branches move to the left side of the imaginary axis that
12
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
means more stability of the system, but with increase in velocity, two of these branches are moving toward imaginary axis that means system tends to the instable. After getting the branches to the imaginary boundaries and passing them, linear system became unstable. Similar forms of the behaviors are presented in [9,10]. Similar calculations for the case 2 of Table 1 are shown in Figs. 6–8. In this case the flutter velocity is UF ¼ 100:924 m=s. Also, the system has a divergence velocity of UD ¼ 123:666m=s. In divergence velocity the frequency of oscillation is zero (a pure real eigenvalue), in which the static aeroelastic phenomenon will occur and aerodynamics forces will dominate to the structural forces. 3.2. Solution convergence study
Fig. 8. Frequency of eigenvalues of the linear aeroelastic system vs. free stream velocity, case 2.
Table 2 Convergence solution analysis of flutter and divergence velocity versus the number of structural mode functions. Number of mode functions
Flutter velocity
Divergence velocity
w
h
v
Case 1
Case 2
Case 1
Case 2
2 2 3 4 4 5
1 2 2 2 3 3
1 1 2 2 2 3
100.922 102.246 102.268 102.271 102.378 102.378
100.947 100.902 100.921 100.924 100.926 100.926
205.144 205.144 205.144 205.144 205.144 205.144
123.666 123.666 123.666 123.666 123.666 123.666
Here the effect of different number of mode functions on the convergence of the solution will be studied. This convergence study is briefly described in Table 2 and Figs. 9–11. In Table 2, the variation of the flutter and divergence velocities are shown. As seen from Table 2, with increase in the number of mode functions, there are not any remarkable changes in the flutter velocity predictions from row 3 of this table. By increasing the numbers of bending and torsional modes form case 2-11 (flap bending, torsional, lag bending) to case 2-2-1, a remarkable change has been occurred in the flutter velocity. But after this, by increasing the mode numbers from 2-2-1 to 3-2-2 and then 4-2-2 there is not many changes in these predictions. The divergence study of the limit cycle amplitudes based on the number of mode functions are shown in Figs. 9–11. These figures show the displacement at the tip of the wing. As seen from these figures, there are minor changes of the predictions based on the 3-2-2 and 4-2-2 number of mode functions. Hence selection of the 4, 2 and 2 mode functions for flap, torsion and lag displacements, respectively, are very justified.
Fig. 9. Flap displacement of the nonlinear aeroelastic wing model in UN ¼102.271 m/s for different combinations of mode functions (tip of the wing) case -1.
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
13
Fig. 10. Torsional displacement of the nonlinear aeroelastic wing model in UN ¼ 102.271 m/s for different combinations of mode functions (tip of the wing) case -1.
Fig. 11. Lag displacement of the nonlinear aeroelastic wing model in UN ¼ 102.271 m/s for different combinations of mode functions (tip of the wing) case -1.
Fig. 12. Wing tip (a) vertical bending deflection, (b) torsion angle and (c) lag bending deflection vs. time for nonlinear solution in UN ¼0.9UF—case 1 (tip of the wing) case -1.
14
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 13. Wing tip (a) vertical bending deflection, (b) torsion angle and (c) lag bending deflection vs. time for nonlinear solution in UN ¼1.1UF—case 1 (tip of the wing) case -1.
Fig. 14. Flap displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -1 (tip of the wing).
Fig. 15. Torsion displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -1 (tip of the wing).
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 16. Lag displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -1 (tip of the wing).
Fig. 17. Flap displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -2 (tip of the wing).
Fig. 18. Torsion displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -2 (tip of the wing).
15
16
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 19. Lag displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -2 (tip of the wing).
Fig. 20. The flap displacement amplitude along the span of the wing for different velocities, case 1.
Fig. 21. The torsion displacement amplitude along the span of the wing for different velocities, case 1.
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
3.3. The amplitude of the limit cycle oscillation After examination of the linear stability of aeroelastic model and selection of the correct number of mode functions, time response analysis of nonlinear aeroelastic model in different velocities will be examined. For example, the time response analysis of two different velocities of the case 1 wing is shown in Figs. 12 and 13. As seen from Fig. 12, in velocity less than the linear flutter velocity, vibration amplitudes decreases with time and finally damps. In this velocity, all of the eigenvalues of the linear system have negative real part, and the linear system is stable. Hence, the amplitude of the oscillation gradually moves to zero. This is by increasing the free
17
stream velocity above the flutter boundary, limit cycle oscillations will occur. In Fig. 13 the reason is that in a velocity greater than linear flutter velocity, the linear system has a complex conjugate positive real part eigenvalues which tend to increase the amplitude of the vibrations. In the absence of the structurally nonlinear term, the amplitude of the oscillations tend to infinity, while with increase in the amplitude of the oscillations, the nonlinear structural terms sufficiently become large and prevent further increase of the vibration amplitude. Hence, finally the oscillation will settle down in a nonlinear absorber which is a limit cycle in this case. Effects of the different free stream velocity in the time response of nonlinear aeroelastic model are shown
Fig. 22. The lag displacement amplitude along the span of the wing for different velocities, case 1.
Fig. 23. The lag static displacement along the span of the wing for different velocities, case 1.
18
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
in Figs. 14–16 for the case 1. These displacements are shown for the tip of the wing. From Fig. 14, with increasing velocity beyond linear flutter velocity, the amplitude of the tip flap displacements increases at first and then decreases. Similar forms of the behaviors are presented for torsion and lag displacements. From Figs. 14 and 15 the flap and torsion displacements are symmetric, but the lag displacement of Fig. 16 is not symmetric. The lag displacement has a static deformation in the direction of the flow and the whole
of the wing will oscillate in the vicinity of this static displacement. According to this figure in velocities greater than flutter velocity, the aerodynamics force is divided in two sections. Some part of the aerodynamics and structural forces counteract each other without inclusion of the inertial forces and produce a static displacement of the wing, and remaining parts of the aerodynamics forces will counteract with structural and inertial forces and produce oscillations around this static state. Hence the limit cycle amplitudes decrease with the increase in
Fig. 24. Wing cross section oscillations at the nodal line of w displacement of 0.5 m for velocity of 1.01UF m/s, case 1.
Fig. 26. Wing cross section oscillations at the nodal line of w displacement of 0.58 m for velocity of 1.01UF m/s, case 1.
Fig. 25. Wing cross section oscillations at the nodal line of w displacement of 0.54 m for velocity of 1.01UF m/s, case 1.
Fig. 27. Wing cross section oscillations at the nodal line of w displacement of 0.54 m for velocity of 1.06UF m/s, case 1.
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
19
Fig. 28. Wing cross section oscillations at the nodal line of w displacement of 0.54 m. for velocity of 1.09UF m/s, case 1.
velocity, the flap displacement at the tip of the wing is nearly constant. With the increasing free stream velocity up to 1.01UF, the torsion displacement decreases, Fig. 21. For the lag displacement according to Fig. 22, the maximum amplitude at the tip of the wing decreases with the increasing free stream velocity up to 1.12UF and then increases. But at the mid of the wing the amplitude of the limit cycle oscillation decreases for velocities greater than the flutter velocity. As it is shown in Fig. 16, there is a certain static displacement in lag displacements. This static displacement of different velocities in along the span of the wing is shown in Fig. 23. According to Fig. 23 to UN ¼1.05UF with the increasing velocity, displacement amplitude increases. But maximum amplitude at the tip of the wing decreases with the increasing velocity. The cross section variations at velocity of the 1.01UF of the case 1 wing are shown in Figs. 24–26. For velocities greater than flutter velocity, there is a nodal line for the w dislacement and the cross section has rotation around the node point of 0:54 m. The cross section variations for points before and after the nodal point of the w displacement are shown in Figs. 24 and 26. The variations of the wing cross section at different velocities are shown in Figs. 27–29. According to these figures, in velocity of 1.06UF, nodal line gets to the end of the wing section and for velocities greater than it, the nodal line is made in the out of cross section. The position of the wing nodal line at different velocities is shown in Fig. 30. Also, it should be mentioned that constitution of nodal line is due to counter reaction of effects of w and torsional displacements, where that the sum of these parameters is zero. The configuration of the wing in velocity of 1.1UF is shown in Fig. 31. As expected, there is a nodal line for the w displacement.
Fig. 29. Wing cross section oscillations at the nodal line of w displacement of 0.54 m. for velocity of 1.15UF m/s, case 1.
Fig. 30. Wing cross section oscillations at the nodal line of w displacement of 0.54 m for velocity of 1.06UF m/s, case 1.
velocity. Similar results for the case 2 of Table 1 are shown in Figs. 17–19. As it is clear in Fig. 17 for lag displacement, with the increasing free stream velocity greater than flutter velocity, limit cycle amplitude increases. Also, as it is shown in Figs. 18 and 19, there is a similar situation for torsion and lag displacements. The maximum amplitudes of flap, torsion and lag displacements along the span of the wing or the configuration of the wing are shown in Figs. 20–22, respectively. According to Fig. 20, for velocities greater than flutter
20
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 31. Wing deflections at velocity of 1.1UF m/s, case 1.
Fig. 32. Phase portrait for the displacement at the tip of the wing tip (a) flap, (b) torsion and (c) lag displacements in UN ¼1.1UF—case 2.
A typical phase portrait for the displacements at the tip of the wing at the velocity of U1 ¼ 1:1UF is shown in Fig. 32. As seen from these figures, the motions gradually settle down in the limit cycle.
3.4. The effects of different aeroelastic parameters on the flutter velocity
inertial term will decrease and hence, the flutter velocity will increase. In case of the increase in xa from 0.5 to 0.5 which is shown in Fig. 33(b), it should be mentioned that with increase in xa, the linear flutter velocity increases. Fig. 33(c) shows the same results for the linear flutter velocity with different values of xa and ah. 4. Conclusion
Finally, the effect of different aeroelastic parameters on the linear flutter velocity of wing will be studied. For this purpose the effect of mass ratio m, dimensionless distance from wing section mid-chord to elastic axis ah and dimensionless distance from elastic axis to the mass center xa will be studied. The variation of the flutter velocity with mass ratio m and ah is shown in Fig. 33(a). As it is clear from Fig. 33(a), for a constant m with the increase in ah from 0.5 to 0.5, the linear flutter velocity decreases, and for constant ah, with the increase in m from 20 to 60, the linear flutter velocity increases. With the increase in mass ratio, the
In this study the aeroelastic characteristics of highaspect-ratio wing in low subsonic flow is considered. The high-aspect-ratio wing is modeled by Hodges–Dowell wing model which exhibit flap—lag and torsion displacement and is a structurally nonlinear wing model. The aerodynamics model was linear. Flutter, divergence and limit cycle amplitudes of different wing sections are studied and the presence of the nodal line in limit cycle oscillation of the wing is clarified. The effect of different aeroelastic properties such as mass ratio, etc. are examined in this study.
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
21
Fig. 33. The effect of different wing parameters on aeroelastic characteristics of wing, case 1.
References [1] D.H. Hodges, E.H. Dowell, Nonlinear Equations of Motion for the Elastic Bending and Torsion of Twisted Nonuniform Rotor Blades, NASA Technical Note, December 1974. [2] D. Tang, E.H. Dowell, Experimental and theoretical study on aeroelastic response of high-aspect-ratio wings, AIAA Journal 39 (8) (2001).
[3] D. Tang, E.H. Dowell, Experimental and theoretical study of gust response for high-aspect-ratio wing, AIAA Journal 40 (3) (2002). [4] M.J. Patil, D.H. Hodges, C.E.S. Cesnik, Nonlinear Aeroelastic Analysis of Aircraft with High Aspect Ratio Wings, AIAA-98-1955. [5] P. Malatkar, Nonlinear Vibrations of Cantilever Beams and Plates, Ph.D. Dissertation, Virginia Polytechnic Institute and State University, 2003.
22
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
[6] C. Nichkawde, 2006, Nonlinear Aeroelastic Analysis of High Aspect Ratio Wings using the Method of Numerical Continuation, M.Sc. Dissertation, Texas A&M University. [7] D.M. Tang, E.H. Dowell, Effects of geometric structural nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings, Journal of Fluids and Structures 19 (2004) 291–306. [8] Y.C. Fung, An Introduction to the Theory of Aeroelasticity, John Wiley and Sons, Inc., 1995.
[9] B.H.K. Lee, L. Liu., K.W. Chung, Airfoil motion in subsonic flow with strong cubic nonlinear restoring forces, Journal of Sound and Vibration 281 (2005) 699–717. [10] B. Ghadiri, M. Razi, Limit cycle oscillations of rectangular cantilever wings containing cubic nonlinearity in an incompressible flow, Journal of Fluids and Structures 23 (2007) 665–680.
Contents lists available at ScienceDirect
Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro
Nonlinear aeroelastic analysis of high-aspect-ratio wings in low subsonic flow K. Eskandary, M. Dardel n, M.H. Pashaei, A.K. Moosavi Department of Mechanical Engineering, Babol University of Technology, P.O. Box 484, Shariati Avenue, Babol, Iran
a r t i c l e in f o
abstract
Article history: Received 9 December 2010 Received in revised form 13 June 2011 Accepted 11 July 2011 Available online 5 August 2011
In this study, aeroelastic characteristics of high-aspect-ratio wing models with structural nonlinearities in quasi-steady aerodynamics flows are investigated. The studied wing model is a cantilever wing with double bending and torsional vibrations and with large deflection ability in accordance with Hodges–Dowell wing model. This wing model is valid for long, straight and thin homogeneous isotropic beams. The aerodynamics model is based on quasi-steady aerodynamic which is valid for aerodynamic flows without wake, viscosity and compressibility effects. The effect of different parameters such as mass ratios and stiffness ratios on flutter and divergence velocities and limit cycle oscillation amplitudes are carefully studied. & 2011 Elsevier Ltd. All rights reserved.
Keywords: High-aspect-ratio wing Flap-lag and torsion Flutter Limit cycle
1. Introduction Helicopter blades and some models of aircrafts have long wings so that flexibility and nonlinear deformations have important roles in them. Structural flexibility of this type of structures is very high, and structures are suffering from large deflections and deformations. Large deflection of these structures is describable by structural nonlinear modeling. Presences of nonlinear structural factors with aerodynamic forces lead to undesirable aeroelastic phenomena. Therefore, careful review and understanding of these phenomena in aeroelastic systems is vital. Modeling of long wings has been done by several researchers which some of them have been mentioned here. Hodges and Dowell developed the equations of motion by two complementary methods, Hamilton’s principle and the Newtonian method [1]. Tang and Dowell constructed an experimental high-aspect-ratio wing aeroelastic model with a slender body at the tip, and measured the response due to flutter and limit-cycle oscillations
n
Corresponding author. E-mail address: [email protected] (M. Dardel).
0094-5765/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.actaastro.2011.07.017
(LCO) in a wind-tunnel test. They also developed a theoretical model and made calculation to correlate with the experimental data [2]. Tang and Dowell presented a nonlinear response analysis of a high-aspect-ratio wing aeroelastic model excited by gust loads theoretically and examined it in a wind-tunnel [3]. Patil, Hodges and Cesnik described a formulation for aeroelastic analysis of aircraft with high-aspect-ratio wings. The analysis was a combination of a geometrically exact mixed formulation for dynamics of moving beams and finite-state unsteady aerodynamics with the ability to model dynamic stall [4]. Malatkar explored the impact of kinematic structural nonlinearities on the dynamics of a highly deformable cantilevered wing. Two different theoretical formulations were presented and analyzed for nonlinear behavior [5]. Nichkawde presented a study of the nonlinear vibrations of metallic cantilever beams and plates subjected to transverse harmonic excitations. Both experimental and theoretical results were presented [6]. In this work, aeroelastic behavior of Hodges–Dowell wing model with incompressible quasi-steady aerodynamic model is investigated. A brief description of the Dowell–Hodges beam is given here. The theory is intended for application to long, straight, slender, homogeneous,
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Nomenclatures ah b Cv, Cw Cy c, c dD, dL dFv, dFw dMx dMxf E e G g h I1, I2 Ia Iy J L
dimensionless distance from wing section mid-chord to elastic axis wing section semi-chord section lag and vertical bending damping section twist damping wing chord and dimensionless chord, c/L section drag and lift forces section lag and vertical component forces section pitch moment about elastic axis section pitch moment about one-quarter chord modulus of elasticity section mass center from elastic axis shear modulus gravitational constant flap displacement vertical, lag area moments mass moment of inertia about elastic axis mass moment of inertia per unit length about center of mass torsional stiffness constant wing span
isotropic beams with large displacements and is accurate to second order based on the restriction that squares of bending slopes, twist, h/L, and c/L are small with respect to unity. Radial nonuniformities (mass, stiffness, twist, etc.), lag offsets of the mass centroid and tension axes from the elastic axis, and warping of the cross section are included. Other more specialized details are not considered, such as blade root feathering flexibility, torque offset, blade sweep, and curvature; nor are configurations considered in which the feathering bearing is replaced with a torsionally flexible strap [1]. The effect of wakes, compressibility and viscosity of incompressible quasi-steady aerodynamic model are ignored. In current study, determination of flutter and divergence velocities and clarifying of limit cycle amplitudes are investigated. Effect of changing mass ratio and distance between elastic axis and center of mass and other parameters are carefully examined in aeroelastic properties. This study is done for the different wing section models to evaluate the effect of mass and inertia on the limit cycle amplitude and flutter velocity. The effects of structural damping are considered in structural modeling.
m ra t U, Un UL Vi, Wi v w x xf xa
mass per unit length of the wing radius of gyration about elastic axis time freestream velocity and dimensionless velocity linear flutter velocity generalized coordinates for bending lag or edgewise bending deflection vertical or flapwise bending deflection position coordinate along wing span position of flexural axis (pivot) dimensionless distance from elastic axis to center of mass y pitch angle of wing section yi generalized coordinates for torsion y0 steady angle of attack at root section m mass ratio z, Z dimensionless flap displacement r air density t dimensionless time o frequency of limit cycle oscillations oz, oZ, oa natural frequencies in flap, lag and pitch on1 , on2 frequency ratio ð Þ0 dð Þ=dx ð: Þ dð Þ=dt
beams. In this study, the wing has no initial twist and initial warping. In current wing model, the elastic and mass axes are not coincided, and structural damping is included in equations of motions. The considered wing model is a cantilevered wing model with bending displacements of w and v and torsion of y, which are shown in Fig. 1. The flap deflection is denoted by w, positive downward direction, lag deflection is denoted by v, positive in stream-wise direction and the torsion angle y, positive nose up. The sketch of wing section is shown in Fig. 2, where c is the chord, b is the semi-chord length, ahb is the distance from wing section mid-chord to elastic axis and xab is the distance from elastic axis to the center of mass. The length of the wing is L. According to Hodges–Dowell equations for a uniform elastic wing and with ignoring the wing section warping,
2. Aeroelastic equations 2.1. Structural equations of motions In this section, a brief description of structural equations of motion of wing model based on Hodges–Dowell wing model is presented. A Hodges–Dowell wing model presents double bending and torsional vibrations. Hodges– Dowell wing model is second-order equations which are valid for long, straight and thin homogeneous isotropic
7
Fig. 1. Flexible deflected wing.
8
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
section. According to quasi-steady aerodynamic theory, aerodynamic moments and forces are as follows [8]: c € Lðx,tÞ ¼ rpb2 w xf y€ þ rpb2 U y_ 2 _! _ 3 y w 2 cxf þ rpU c y þ þ ð7Þ 4 U U c c _ w xf y_ Mxf ðx,tÞ ¼ rpb2 xf 2 2 rpb4 _ 3 y cxf rpb2 U y_ 4 8 _! _ 3 y 1 w 2 2 cxf þ rpU ec y þ þ rpUc3 y_ 4 16 U U ð8Þ
Fig. 2. Schematic figure of the wing section.
where the aerodynamics center of e and distance between elastic axis to the wing mid-chord are defined by
equations of motions will be as follows [7]: dF € _ ¼ w € EI1 wð4Þ þ ðEI2 EI1 Þðyv00 Þ00 þmwmx a by þ Cw w dx
ð1Þ
e¼
dMx dx
ð2Þ
In addition, by ignoring the drag force for this low speed, inviscid aerodynamics model, we can write
€ þ Cy y_ ¼ GJ y þ ðEI2 EI1 Þw00 v00 þmra2 y€ mxa bw 00
EI2 vð4Þ þðEI2 EI1 Þðyw00 Þ00 þmv€ þCv v_ ¼
dFv dx
xf 1 , c 4
xf ¼
b þ ah b 2
dD ¼ 0 ð3Þ
where m, Iy, EI1, EI2, GJ, Cw,Cy and Cv are the mass per unit length of the wing, mass moment of inertia per unit length about mass center, flap and lag bending stiffness, torsional stiffness, the flap, torsional and lag structural damping coefficients, respectively. e is the distance between mass and elastic axes and x is the position coordinate along the wing span. dFw/dx, dMx/dx and dFv/dx are the aeroedynamics force and moments, which will be described as follows. Important structural nonlinear terms have been included in Eqs. (1)–(3) and third- or higher-order nonlinear terms have been ignored. 2.2. Quasi-steady aerodynamics model According to the work presented in [2,3,7], aerodynamic forces and moment are as follows: dFw ¼ dL þ y0 dD
ð4Þ
dMx ¼ dMxf
ð5Þ
dFv ¼ dD þ y0 dL
ð6Þ
where y0 is a constant angle of attack (which is assumed to be zero), L(t) is the aerodynamic lift force, M(t) is the aerodynamic moment, D(t) is the aerodynamic drag force. In works done by Dowell et.al. aerodynamic force has been made in accordance with Onera aerodynamics model. In present work, quasi-steady aerodynamic model has been used. Also, in this aerodynamics model, the flow is assumed to be incompressible, inviscid and irrotational. Quasi-steady models assume that there are only four contributions to the aerodynamic forces, horizontal air velocity U, at angle y(t) to cross section, wing flap velocity, _ wðtÞ, normal component of torsion velocity, ðxf xÞy_ ðtÞ, and local velocity induced by the vorticity around the airfoil. xf is the elastic axis position of the wing cross
ð9Þ
ð10Þ
In current study, the angle of attack is assumed to be zero. 2.3. Discretizing aeroelastic equations using assumed mode method With substituting aerodynamics forces and moment according to Eqs. (7)–(9) in structural equations of motions of (1)–(3), the complete aeroelastic equations will be obtained. These equations can be descritized in accordance with Galerkin’s method. In this method, at first the w, v and y displacements are expressed in terms of to the assumed mode method. Then with substituting these displacements in equations of motions and multiplying each equation with an admissible mode function and integrating along the whole area of the wing, the equations of motions will be descritized. In assumed mode method, displacements are written based on the products of time dependent generalized coordinates and mode functions which satisfy the geometric boundary conditions. In accordance with the assumed mode method the bending and torsional displacements can be described as follows: X wðx,tÞ ¼ Wi ðtÞci ðxÞ ð11Þ
yðx,tÞ ¼ vðx,tÞ ¼
X X
yi ðtÞfi ðxÞ
ð12Þ
Vi ðtÞai ðxÞ
ð13Þ
where Wi(t), yi(t) and Vi(t) are the generalized coordinates and ci(x), fi(x) and ai(x) are the mode functions that satisfy the geometric boundary conditions. By combining structural Eqs. (1)–(3) with aerodynamics Eqs. (7)–(9) and multiplying both sides of the equation in a desired mode shape and integrating along the wing, the original partial differential equations are converted to the ordinary differential equations. For further simplification of the equations, reducing the amount of the necessary calculations and selecting the
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
appropriate mode shape, high-order derivative terms are converted to lower-order derivative terms by integrating by part method. The resulting equations are as follows: Z L X XX Wj EI1 ci cð4Þ Vj yk j dx þ ðEI2 EI1 Þ 0
Z
L
0
þm
X
þCw
ci a00j f00k dx þ2 € j W
X
¼ rpb2
Z Z
þ rpb2 U
X
Z
yj
L 0
Z
þCy
X
y_ j
y€ j Z
Z
¼ rpb2 xf xf
0
ci fj dx þ y_ j
Z
X
L
ci fj dx
0
Z
L
ci cj dx
0
ð14Þ
XX X
L
fi fj dxme
0
Z
_ j W
Wj Vk € j W
Z
Z
L 0
L
fi cj dx
0
fi fj dx
2 c X
y€ j
Z
Z
L 0
L
fi cj dx
fi fj dx
0
rpb4 X €
yj
8
X Z L 3 cxf rpb2 U y_ j fi fj dx 4 0 X Z L Z X _ j þ rpUec2 U yj fi fj dx þ W þ
X
Vj
X
y_ j
X 1 rpUc3 y_ j 16
EI2
3 cxf 4
Z
L 0
L
0
þm
Z
L
fi fj dx
0
X
00
V€ j
¼ y0 rpb2
þ rpb2 U X
_ j W
Z 0
00
X Z 0
L
Z
L
0
L
ai aj dx þ Cv
X
€ j W
y_ j
Z
Z
0
L 0
L
0
fi cj dx
ð15Þ
XX 000
Wj yk
0
ai cj fk dx þ X
V_ j
Z 0
Z 0
L
ð4Þ
ai cj fk dx
L
The aeroelastic Eqs. (14)–(16) are written in non-dimensional form by defining the following parameters: rffiffiffiffiffiffiffiffiffi U t w v Ia m t ¼ 1 , z ¼ , Z ¼ , ra ¼ , , m¼ b b b mb2 rpb2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u R u GJ 1 ðf 00 Þ2 dx U , oa ¼ t 2 R01 1 2 , Un ¼ boa Iy L 0 ðf1 Þ dx vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u R 1 00 2 u EI ðc Þ dx oz ¼ t 14 R01 1 2 mL 0 ðc1 Þ dx vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u R 1 00 2 u EI ða Þ dx o o x oZ ¼ t 24 R01 1 2 , on1 ¼ z , on2 ¼ Z , x ¼ L oa oa mL ða1 Þ dx ð19Þ where z and Z are the dimensionless flap and lag deflections, oz, oZ and oa are the first natural frequencies of uncoupled flap, lag and torsional modes, U*, oni , m and ra are the dimensionless velocity, bending to torsion frequency ratios, mass ratio and the radius of gyration about the elastic axis, respectively, and x is the dimensionless coordinate along the span of the wing. These dimensionless parameters are defined in order to select the necessary parameters in accordance with the values given in [9], to validate the presented model in accordance the flutter velocities in this reference. With substituting these dimensionless parameters in Eqs. (14)–(16) and multiplying these equations in ci, fi and ai, respectively, the following equations of motion will be obtained:
on12
ai aj dx
Z L c X € ai cj dx xf yj ai fj dx 2 0
X Z ai fj dx þ rpUc U yj
ð18Þ
fi fj dx
fi fj dx
0
an ¼ ðcos bn cosh bn Þ=ðsin bn sinh bn Þ
0
L
L
0
Z
L
ai cj fk dx þ2
þ
0
Z
ai að4Þ j dx þðEI2 EI1 Þ
Z
fi c00j a00k dx
L
€ j W
bn ¼ 1:8751, 4:6941, 7:8547, 10:9955, 10:9955þ p,
And for torsional modes: 2n1 x fn ¼ sin p 2 L
ci fj dxrpUc
c X
2
ci fj dx
L
L
0
L
0
Z L c X € ci cj dx xf yj ci fj dx 2 0
L 0
X
X
y€ j
0
Appropriate mode shape for description of the displacements of w, v, and y is obtained from beams and bar with equivalent boundary conditions to the wing model. For descriptions of w and v displacements, Euler–Bernoulli’s cantilevered beam mode functions are used and for description of torsional displacement of y, mode functions of fixed– free bar are used. These mode functions are as follows: For bending modes: x x x x ð17Þ cn ¼ sin bn sinh bn þ an cos bn cosh bn L L L L That
fi f00j dx þðEI2 EI1 Þ
2 þmKm
Z
ci að4Þ fk dx j
ci cj dx
0
3 cxf þ 4 GJ
X
y_ j
X Z U yj
0
L
L
€ j W
X
Z
ci a000j f0k dx þ
ci cj dxme
0
X
L
L 0
_ j W
Z
9
U
R1 0
n2 R 1
N c21 dx X
002 0 c1 dx j ¼ 1
m
ai fj dx
X Z L 3 cxf ai cj dx þ y_ j ai fj dx 4 0 ð16Þ
Z
1 0
N Z a X xa h y€ j
L 0
zj
Z
0
¼
j¼1
1
ci cj dx
c00i c00j dx þ 1 1
0
ci fj dx þ
N ð22ah Þ X y_
m
j¼1
Z
1
X N
m
0
Z
j¼1
1
j
z€ j
1 0
ci cj dx
N Cw 2 X z_ n moa U m j¼1 j
ci fj dx
n2 R 1 2 N X N o1 0 c1 dx X EI2 1 Zy EI1 U n2 R 1 c’’2 dx j ¼ 1 k ¼ 1 j k 0
1
N 2X
mj¼1
yj
Z 0
1
ci fj dx
10
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Z
1
0
R1
2
0
1
!
N X
1 0 f1 dx 4e R U n2 1 f02 dx ra2 m N X
Z
y€ j
þ
!
1
c0i a000j fk dx
0
yj
Z
1
0
j¼1
1
fi fj dx
0
j¼1
Z
c0i a00j f0k dxþ
ð21Þ ða2 þ ð1=8ÞÞ 1þ h 2 ra m
f0i f0j dx þ
ah xa þ ra2 m ra2
X N
z€ j
Z
1
fi cj dx
0
j¼1
!
N Z 1 Cy ðð1=2Þ þ ð4e1Þðð1=2Þah ÞÞ X y_ j fi fj dx Iy oa U n ra2 m 0 j¼1
Z 1 N 4e X EI ¼ 2 z_ j fi cj dx 2 1 EI1 ra m j ¼ 1 0 R1
o1n2
r 2 U n2 a
o2n2 U n2
R1
0
N X N c21 dx X
0
c002 1 dx j ¼ 1 k ¼ 1
R1
0
N a21 dx X
0
a002 1 dx j ¼ 1
R1
þ
þ
Z
Zj
moa U n j ¼ 1
m
yj
j¼1
Z 0
Z 0
1
fi c00j a00k dx
N X
Z€ j
Z
j¼1
Z 0
1
0
ai aj dx Z
m j¼1
ai fj dx
1
ai fj dx
ð22Þ
1
N y0 X ai aj dx z€ j
1
0
j¼1
N 2y0 X
a00i a00j dx þ
Z_ j
Z
N y0 ah X y€ j
m
1 0
N X
Cv
z j Zk
0
j¼1
Z N 2y0 X z_
m
j
j¼1
Property
Value
Wing
Case 1
Case 2
Span (L) (m) Chord (c) (m) Dimensionless distance from wing section mid-chord to elastic axis (ah) Dimensionless distance from elastic axis to center of mass (xa) Flap structural modal damping (Cw) Lag structural modal damping (Cv) Torsional structural modal damping (Cy) Mass ratio (m) Radius of gyration about elastic axis (ra) Natural frequencies in flap bending (oz) (rad/s) Natural frequencies in lag bending (oZ) (rad/s) Natural frequencies in torsion (oa) (rad/s)
0.6299 0.1018 0.3
0.6299 0.1018 0.19
0.22
0.02
0.02 0.025 0.031 48 0.558 26 p
0.02 0.025 0.031 36.8 0.558 18 p
172.6 p 254.6 p 210 p 180 p
Flight conditions Density of air (r) (kg/m3)
1.1
1.1
1
ai cj dx
N y0 ð22ah Þ X y_ j
m
Table 1 Wing model data.
Z
1 0
ai fj dx
1 0
ai cj dx
R1 N X N EI1 on22 0 a21 dx X ¼ 1 zy R EI2 U n2 1 a002 dx j ¼ 1 k ¼ 1 j k 0 1 ! Z Z 1
0
00
1
0
a0i cj fk dx þ
0
000
a0i cj fk dx
ð23Þ
The interval of the integral is [0,1]. Then the whole aeroelastic Eqs.(21)–(23) can be written briefly as follows: € þ½Ceq fqg _ þ ½Keq fqg ¼ fFN g ½Meq fqg
ð24Þ
where [Meq], [Ceq], and [Keq] are mass, damping and stiffness matrices, respectively, and {FN} is a vector of structural nonlinearity. The vector {q} is defined as T fqg ¼ fzyZg . Now after presenting the whole aeroelastic equations, the validity and results of these equations will be examined.
Fig. 3. Damping of eigenvalues of the linear aeroelastic system vs. free stream velocity, case 1.
3. Aeroelastic results In this section, the aeroelastic results of the presented wing model for different cases will be studied. For this purpose, at first the wing structural and aerodynamics parameters will be selected. Here, two different wing section parameters are selected, which are shown in Table 1. These parameters are selected in compatibility with those parameters presented in [2,9]. These parameters are not the same to those given in these references, since in [2] the parameters are for a wing which represent stall phenomenon and parameters given in [9] are for a simple wing model with flap and pitch motions and incompressible, irrotational aerodynamics flow, which are different with those presented in current work. Hence
combinations of the parameters given by these references are given in Table 1. 3.1. Eigenvalue solution of the linear aeroelastic system At first eigenvalue solution of the linear aeroelastic system which determine the stability of the linear part of the system, will be studied. Changes of eigenvalues of the linear aeroelastic equation with free stream velocity for the case 1 of Table 1, have been shown in Figs. 3–5. Damping and frequency of eigenvalues are given by the real and imaginary parts of eigenvalues l of the linear system. Different combinations of mode functions are
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 4. Frequency vs. damping of the eigenvalues of the linear aeroelastic system, case 1.
11
Fig. 6. Damping of eigenvalues of the linear aeroelastic system vs. free stream velocity, case 2.
Fig. 5. Frequency of eigenvalues of the linear aeroelastic system vs. free stream velocity, case 1.
Fig. 7. Frequency vs. damping of the eigenvalues of the linear aeroelastic system, case 2.
selected, and the convergences of the eigenvalue solution are studied. For this figures, the number of mode functions for flap, lag and torsion displacements are 4, 2 and 2, respectively. With increase in the number of mode functions above these mentioned mode numbers, minor changes will occur in the predictions. A sample convergence of the solution based on the number of the mode functions will be presented in the next section. As seen in Fig. 3 with increase in the free stream velocity, one of the branches of the real part of the eigenvalues of the linear system will change its sign and become positive. In this case, the initially stable system becomes unstable. In flutter velocity a branch of
eigenvalues has intersection with the real axis. For this intersection, the real part of corresponded eigenvalue is zero. For the case 1 of Table 1 the flutter velocity is UF ¼102.271 m/s. From Fig. 4 the frequency of the flutter velocity is about 83:75 Hz. From Fig. 4 it is clear that two branches of eigenvalues are approaching toward each other gradually and near the flutter velocity they have the less distance from each other that show the bending–torsion flutter mode of the presented wing model. Between these branches, an internal resonance occurs. Also from Fig. 4 with increase in the free stream velocity, at first all of the branches move to the left side of the imaginary axis that
12
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
means more stability of the system, but with increase in velocity, two of these branches are moving toward imaginary axis that means system tends to the instable. After getting the branches to the imaginary boundaries and passing them, linear system became unstable. Similar forms of the behaviors are presented in [9,10]. Similar calculations for the case 2 of Table 1 are shown in Figs. 6–8. In this case the flutter velocity is UF ¼ 100:924 m=s. Also, the system has a divergence velocity of UD ¼ 123:666m=s. In divergence velocity the frequency of oscillation is zero (a pure real eigenvalue), in which the static aeroelastic phenomenon will occur and aerodynamics forces will dominate to the structural forces. 3.2. Solution convergence study
Fig. 8. Frequency of eigenvalues of the linear aeroelastic system vs. free stream velocity, case 2.
Table 2 Convergence solution analysis of flutter and divergence velocity versus the number of structural mode functions. Number of mode functions
Flutter velocity
Divergence velocity
w
h
v
Case 1
Case 2
Case 1
Case 2
2 2 3 4 4 5
1 2 2 2 3 3
1 1 2 2 2 3
100.922 102.246 102.268 102.271 102.378 102.378
100.947 100.902 100.921 100.924 100.926 100.926
205.144 205.144 205.144 205.144 205.144 205.144
123.666 123.666 123.666 123.666 123.666 123.666
Here the effect of different number of mode functions on the convergence of the solution will be studied. This convergence study is briefly described in Table 2 and Figs. 9–11. In Table 2, the variation of the flutter and divergence velocities are shown. As seen from Table 2, with increase in the number of mode functions, there are not any remarkable changes in the flutter velocity predictions from row 3 of this table. By increasing the numbers of bending and torsional modes form case 2-11 (flap bending, torsional, lag bending) to case 2-2-1, a remarkable change has been occurred in the flutter velocity. But after this, by increasing the mode numbers from 2-2-1 to 3-2-2 and then 4-2-2 there is not many changes in these predictions. The divergence study of the limit cycle amplitudes based on the number of mode functions are shown in Figs. 9–11. These figures show the displacement at the tip of the wing. As seen from these figures, there are minor changes of the predictions based on the 3-2-2 and 4-2-2 number of mode functions. Hence selection of the 4, 2 and 2 mode functions for flap, torsion and lag displacements, respectively, are very justified.
Fig. 9. Flap displacement of the nonlinear aeroelastic wing model in UN ¼102.271 m/s for different combinations of mode functions (tip of the wing) case -1.
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
13
Fig. 10. Torsional displacement of the nonlinear aeroelastic wing model in UN ¼ 102.271 m/s for different combinations of mode functions (tip of the wing) case -1.
Fig. 11. Lag displacement of the nonlinear aeroelastic wing model in UN ¼ 102.271 m/s for different combinations of mode functions (tip of the wing) case -1.
Fig. 12. Wing tip (a) vertical bending deflection, (b) torsion angle and (c) lag bending deflection vs. time for nonlinear solution in UN ¼0.9UF—case 1 (tip of the wing) case -1.
14
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 13. Wing tip (a) vertical bending deflection, (b) torsion angle and (c) lag bending deflection vs. time for nonlinear solution in UN ¼1.1UF—case 1 (tip of the wing) case -1.
Fig. 14. Flap displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -1 (tip of the wing).
Fig. 15. Torsion displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -1 (tip of the wing).
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 16. Lag displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -1 (tip of the wing).
Fig. 17. Flap displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -2 (tip of the wing).
Fig. 18. Torsion displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -2 (tip of the wing).
15
16
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 19. Lag displacement of the nonlinear aeroelastic wing model in different free stream velocity for the case -2 (tip of the wing).
Fig. 20. The flap displacement amplitude along the span of the wing for different velocities, case 1.
Fig. 21. The torsion displacement amplitude along the span of the wing for different velocities, case 1.
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
3.3. The amplitude of the limit cycle oscillation After examination of the linear stability of aeroelastic model and selection of the correct number of mode functions, time response analysis of nonlinear aeroelastic model in different velocities will be examined. For example, the time response analysis of two different velocities of the case 1 wing is shown in Figs. 12 and 13. As seen from Fig. 12, in velocity less than the linear flutter velocity, vibration amplitudes decreases with time and finally damps. In this velocity, all of the eigenvalues of the linear system have negative real part, and the linear system is stable. Hence, the amplitude of the oscillation gradually moves to zero. This is by increasing the free
17
stream velocity above the flutter boundary, limit cycle oscillations will occur. In Fig. 13 the reason is that in a velocity greater than linear flutter velocity, the linear system has a complex conjugate positive real part eigenvalues which tend to increase the amplitude of the vibrations. In the absence of the structurally nonlinear term, the amplitude of the oscillations tend to infinity, while with increase in the amplitude of the oscillations, the nonlinear structural terms sufficiently become large and prevent further increase of the vibration amplitude. Hence, finally the oscillation will settle down in a nonlinear absorber which is a limit cycle in this case. Effects of the different free stream velocity in the time response of nonlinear aeroelastic model are shown
Fig. 22. The lag displacement amplitude along the span of the wing for different velocities, case 1.
Fig. 23. The lag static displacement along the span of the wing for different velocities, case 1.
18
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
in Figs. 14–16 for the case 1. These displacements are shown for the tip of the wing. From Fig. 14, with increasing velocity beyond linear flutter velocity, the amplitude of the tip flap displacements increases at first and then decreases. Similar forms of the behaviors are presented for torsion and lag displacements. From Figs. 14 and 15 the flap and torsion displacements are symmetric, but the lag displacement of Fig. 16 is not symmetric. The lag displacement has a static deformation in the direction of the flow and the whole
of the wing will oscillate in the vicinity of this static displacement. According to this figure in velocities greater than flutter velocity, the aerodynamics force is divided in two sections. Some part of the aerodynamics and structural forces counteract each other without inclusion of the inertial forces and produce a static displacement of the wing, and remaining parts of the aerodynamics forces will counteract with structural and inertial forces and produce oscillations around this static state. Hence the limit cycle amplitudes decrease with the increase in
Fig. 24. Wing cross section oscillations at the nodal line of w displacement of 0.5 m for velocity of 1.01UF m/s, case 1.
Fig. 26. Wing cross section oscillations at the nodal line of w displacement of 0.58 m for velocity of 1.01UF m/s, case 1.
Fig. 25. Wing cross section oscillations at the nodal line of w displacement of 0.54 m for velocity of 1.01UF m/s, case 1.
Fig. 27. Wing cross section oscillations at the nodal line of w displacement of 0.54 m for velocity of 1.06UF m/s, case 1.
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
19
Fig. 28. Wing cross section oscillations at the nodal line of w displacement of 0.54 m. for velocity of 1.09UF m/s, case 1.
velocity, the flap displacement at the tip of the wing is nearly constant. With the increasing free stream velocity up to 1.01UF, the torsion displacement decreases, Fig. 21. For the lag displacement according to Fig. 22, the maximum amplitude at the tip of the wing decreases with the increasing free stream velocity up to 1.12UF and then increases. But at the mid of the wing the amplitude of the limit cycle oscillation decreases for velocities greater than the flutter velocity. As it is shown in Fig. 16, there is a certain static displacement in lag displacements. This static displacement of different velocities in along the span of the wing is shown in Fig. 23. According to Fig. 23 to UN ¼1.05UF with the increasing velocity, displacement amplitude increases. But maximum amplitude at the tip of the wing decreases with the increasing velocity. The cross section variations at velocity of the 1.01UF of the case 1 wing are shown in Figs. 24–26. For velocities greater than flutter velocity, there is a nodal line for the w dislacement and the cross section has rotation around the node point of 0:54 m. The cross section variations for points before and after the nodal point of the w displacement are shown in Figs. 24 and 26. The variations of the wing cross section at different velocities are shown in Figs. 27–29. According to these figures, in velocity of 1.06UF, nodal line gets to the end of the wing section and for velocities greater than it, the nodal line is made in the out of cross section. The position of the wing nodal line at different velocities is shown in Fig. 30. Also, it should be mentioned that constitution of nodal line is due to counter reaction of effects of w and torsional displacements, where that the sum of these parameters is zero. The configuration of the wing in velocity of 1.1UF is shown in Fig. 31. As expected, there is a nodal line for the w displacement.
Fig. 29. Wing cross section oscillations at the nodal line of w displacement of 0.54 m. for velocity of 1.15UF m/s, case 1.
Fig. 30. Wing cross section oscillations at the nodal line of w displacement of 0.54 m for velocity of 1.06UF m/s, case 1.
velocity. Similar results for the case 2 of Table 1 are shown in Figs. 17–19. As it is clear in Fig. 17 for lag displacement, with the increasing free stream velocity greater than flutter velocity, limit cycle amplitude increases. Also, as it is shown in Figs. 18 and 19, there is a similar situation for torsion and lag displacements. The maximum amplitudes of flap, torsion and lag displacements along the span of the wing or the configuration of the wing are shown in Figs. 20–22, respectively. According to Fig. 20, for velocities greater than flutter
20
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
Fig. 31. Wing deflections at velocity of 1.1UF m/s, case 1.
Fig. 32. Phase portrait for the displacement at the tip of the wing tip (a) flap, (b) torsion and (c) lag displacements in UN ¼1.1UF—case 2.
A typical phase portrait for the displacements at the tip of the wing at the velocity of U1 ¼ 1:1UF is shown in Fig. 32. As seen from these figures, the motions gradually settle down in the limit cycle.
3.4. The effects of different aeroelastic parameters on the flutter velocity
inertial term will decrease and hence, the flutter velocity will increase. In case of the increase in xa from 0.5 to 0.5 which is shown in Fig. 33(b), it should be mentioned that with increase in xa, the linear flutter velocity increases. Fig. 33(c) shows the same results for the linear flutter velocity with different values of xa and ah. 4. Conclusion
Finally, the effect of different aeroelastic parameters on the linear flutter velocity of wing will be studied. For this purpose the effect of mass ratio m, dimensionless distance from wing section mid-chord to elastic axis ah and dimensionless distance from elastic axis to the mass center xa will be studied. The variation of the flutter velocity with mass ratio m and ah is shown in Fig. 33(a). As it is clear from Fig. 33(a), for a constant m with the increase in ah from 0.5 to 0.5, the linear flutter velocity decreases, and for constant ah, with the increase in m from 20 to 60, the linear flutter velocity increases. With the increase in mass ratio, the
In this study the aeroelastic characteristics of highaspect-ratio wing in low subsonic flow is considered. The high-aspect-ratio wing is modeled by Hodges–Dowell wing model which exhibit flap—lag and torsion displacement and is a structurally nonlinear wing model. The aerodynamics model was linear. Flutter, divergence and limit cycle amplitudes of different wing sections are studied and the presence of the nodal line in limit cycle oscillation of the wing is clarified. The effect of different aeroelastic properties such as mass ratio, etc. are examined in this study.
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
21
Fig. 33. The effect of different wing parameters on aeroelastic characteristics of wing, case 1.
References [1] D.H. Hodges, E.H. Dowell, Nonlinear Equations of Motion for the Elastic Bending and Torsion of Twisted Nonuniform Rotor Blades, NASA Technical Note, December 1974. [2] D. Tang, E.H. Dowell, Experimental and theoretical study on aeroelastic response of high-aspect-ratio wings, AIAA Journal 39 (8) (2001).
[3] D. Tang, E.H. Dowell, Experimental and theoretical study of gust response for high-aspect-ratio wing, AIAA Journal 40 (3) (2002). [4] M.J. Patil, D.H. Hodges, C.E.S. Cesnik, Nonlinear Aeroelastic Analysis of Aircraft with High Aspect Ratio Wings, AIAA-98-1955. [5] P. Malatkar, Nonlinear Vibrations of Cantilever Beams and Plates, Ph.D. Dissertation, Virginia Polytechnic Institute and State University, 2003.
22
K. Eskandary et al. / Acta Astronautica 70 (2012) 6–22
[6] C. Nichkawde, 2006, Nonlinear Aeroelastic Analysis of High Aspect Ratio Wings using the Method of Numerical Continuation, M.Sc. Dissertation, Texas A&M University. [7] D.M. Tang, E.H. Dowell, Effects of geometric structural nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings, Journal of Fluids and Structures 19 (2004) 291–306. [8] Y.C. Fung, An Introduction to the Theory of Aeroelasticity, John Wiley and Sons, Inc., 1995.
[9] B.H.K. Lee, L. Liu., K.W. Chung, Airfoil motion in subsonic flow with strong cubic nonlinear restoring forces, Journal of Sound and Vibration 281 (2005) 699–717. [10] B. Ghadiri, M. Razi, Limit cycle oscillations of rectangular cantilever wings containing cubic nonlinearity in an incompressible flow, Journal of Fluids and Structures 23 (2007) 665–680.