Reduced order nonlinear aeroelasticity of swept composite wings using compressible indicial unsteady aerodynamics

Journal of Fluids and Structures 92 (2020) 102812

Contents lists available at ScienceDirect

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

Reduced order nonlinear aeroelasticity of swept composite wings using compressible indicial unsteady aerodynamics ∗

Touraj Farsadi a,b , Mohammad Rahmanian c , , Altan Kayran a a

METUWind, Department of Aerospace Engineering, Middle East Technical University, 06800, Ankara, Turkey Department of Aerospace Engineering, Adana Alparslan Türkeş Science and Technology University, Adana, Turkey c Department of Mechanical Engineering, Jahrom University, P.O. Box: 74135-111, Jahrom, Iran b

article

info

Article history: Received 19 August 2019 Received in revised form 7 October 2019 Accepted 14 November 2019 Available online xxxx Keywords: Sweep angle Indicial aerodynamics Thin walled beams Aeroelastic instability Limit cycle oscillations

a b s t r a c t Nonlinear dynamic aeroelasticity of composite wings in compressible flows is investigated. To provide a reasonable model for the problem, the composite wing is modeled as a thin walled beam (TWB) with circumferentially asymmetric stiffness layup configuration. The structural model considers nonlinear strain displacement relations and a number of non-classical effects, such as transverse shear and warping inhibition. Geometrically nonlinear terms of up to third order are retained in the formulation. Unsteady aerodynamic loads are calculated according to a compressible model, described by indicial function approximations in the time domain. The aeroelastic system of equations is augmented by the differential equations governing the aerodynamics lag states to derive the final explicit form of the coupled fluid-structure equations of motion. The final nonlinear governing aeroelastic system of equations is solved using the eigenvectors of the linear structural equations of motion to approximate the spatial variation of the corresponding degrees of freedom in the Ritz solution method. Direct time integrations of the nonlinear equations of motion representing the full aeroelastic system are conducted using the well-known Runge–Kutta method. A comprehensive insight is provided over the effect of parameters such as the lamination fiber angle and the sweep angle on the stability margins and the limit cycle oscillation behavior of the system. Integration of the interpolation method employed for the evaluation of compressible indicial functions at any Mach number in the subsonic compressible range to the derivation process of the third order nonlinear aeroelastic system of equations based on TWB theory is done for the first time. Results show that flutter speeds obtained by the incompressible unsteady aerodynamics are not conservative and as the backward sweep angle of the wing is increased, post-flutter aeroelastic response of the wing becomes more well-behaved. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Numerous applications of composite thin walled beams (TWBs) are reported in a variety of engineering applications. The advanced theory of beams with different materials (functionally graded and composites) is investigated in several studies (Chaabane et al., 2019; Bourada et al., 2019; Atmane et al., 2017). The main reason for such a pervasive utilization of composites in the fields of aerospace, industrial and mechanical engineerings is their high specific strength and specific ∗ Corresponding author. E-mail address: [email protected] (M. Rahmanian). https://doi.org/10.1016/j.jfluidstructs.2019.102812 0889-9746/© 2019 Elsevier Ltd. All rights reserved.

2

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

stiffness ratios (Chikh et al., 2017). Introduction of composite materials in the past century has made it possible to perform tailoring operations on the structural characteristics of thin walled beams to maximize their efficiency with minimal weight penalties (Weisshaar, 1981; Shirk et al., 1986). Employing such modified thin walled structures in the construction of aerospace vehicles as well as the utilization of small margins of safety in the design process yields minimal weights with an accompanied major penalty of larger structural deformations. In the context of thin walled beams, which are used to model advanced wing configurations, the most appreciable publications in the last three decades are those of Librescu and his colleagues (Librescu and Song, 1992; Qin, 2001; Qin et al., 2002; Na et al., 2007; Librescu et al., 2008; Choo et al., 2012; Na et al., 2011). The advanced refined theory of thin-walled beams of arbitrary closed cross section, incorporating non-classical effects is presented by Librescu and Song (1991). It should be noted that the non-classical effects such as transverse shear effect (Zine et al., 2018; Kaci et al., 2018) could be important for thicker sections and the restrained warping could be important in elastically coupled blades such as in bend-twist coupled blades. The same authors have then used the theory of thin walled beams to investigate the subcritical aeroelastic response of swept anisotropic wing structures (Librescu and Song, 1992). Librescu and Song (2006) and Suresh and Nagaraj (1996) have extended the basic theory for TWBs to include several different complexities and verified the outcomes for different benchmark problems. In a set of two papers, Jung et al. (1999a, 2002b) have provided a thorough assessment and refined the structural model for composite rotor blades with elastic couplings. Thin and thick walled composite rotors are studied as well as either open-section and closed-section blades with arbitrary section shape, stacking sequence and restraint effects. Other complexities such as out-of-plane warping, warping restraint, transverse shear and large deformation are also taken into account. In one of their studies (Jung et al., 2002b) a semi-complementary energy functional is used to derive the beam force–displacement relations. Wall thickness effects are shown to become significant when the thicknessto-depth ratio of the beam reaches around 20%. Other factors such as the slenderness ratio and the layup angle are also shown to have non-negligible effects on the transverse shear behavior of the beam. Nonlinear formulation of TWBs has been studied in the literature for isotropic (Moore, 1986; Attard, 1986; van Erp et al., 1988) and composite structures (Bhaskar and Librescu, 1995). Bhaskar and Librescu (1995) utilized the Lagrangian description and Hamilton’s principle to obtain the geometrically nonlinear formulation of the problem. Such structural models can handle geometrical non-uniformities such as, taper ratio, pretwist and sweep (Farsadi et al., 2019; Farsadi and Hasbestan, 2019). The external store and engine thrust effects may be readily included in the governing equations. Wang and Qin (2016) investigated the nonlinear dynamical behavior of composite TWBs for the special case of simultaneous 1:2 internal and 1:1 external resonances. Their final equations of motion are then solved by the extended Galerkin approach and the method of multiple scales. In general, unsteady aerodynamic loads can be derived using different types of aerodynamic theories such as the strip theory, unsteady vortex-lattice method, indicial response theory, Euler/Navier–Stokes CFD techniques, and the classical and widespread doublet-lattice method (DLM). Although, Unsteady Vortex-Lattice Method and CFD methods are very precise and accurate in calculating aeroelastic response of lifting surfaces, these methods are intrinsically complex and high computational power is needed in aeroelastic analysis. Therefore, time domain unsteady indicial aerodynamic models including compressibility effects offer fast and reasonably accurate solutions to unsteady aeroelastic problems provided that linear aerodynamics suffices. In addition, the present aeroelastic computer algorithm (with indicial aerodynamics rather than CFD or UVLM) may be readily used in a vibration control program or can also work perfectly in conjunction with optimization programs where swiftness is an important factor to make simulations affordable. Ronch et al. (2018) calculated indicial aerodynamic functions in a theoretical manner by modifying the indicial functions of the incompressible flow obtained by Queijo et al. (1978) and combining their results with those of Leishman (2006). Their formulation extends the well-known aerodynamic theories, which are limited to thin airfoils in incompressible flow, to generic trapezoidal wing planforms. The incompressible indicial function concept could be extended to the compressible indicial function concept using the effective Mach number considering the sweep angle and the Prandtl–Glauert compressibility factor in the linear and nonlinear transonic regime while the latter does not fall in the scope of the present study. Aeroelastic analysis of long and slender composite wings is presented by Cesnik et al. (1996). They used an asymptotically correct cross sectional formulation for the structural part and a new 2D unsteady inflow finite state theory which is originally developed by Peters et al. (1995) . Some studies available in the literature (Jung et al., 1999b, 2002a) focused on modeling techniques of composite rotor blades and take all the restraint effects, out of plane warping and other non-classical effects into account. Carrera unified formulation which is known as CUF has also been used in aeroelastic analysis of composite wings/structures (Carrera et al., 2014). CUF is capable of modeling the structure with minimal degrees of freedom while the aerodynamic loading may be obtained by any of the available models. Petrolo exploited 1D structural model obtained by the CUF and coupled it with the doublet lattice method to investigate the flutter behavior of composite lifting surfaces (Petrolo, 2013). In a comprehensive study by Qin (2001), the refined analytical anisotropic TWB model is utilized and aeroelastic instability, dynamic aeroelastic response and active/passive aeroelastic control of advanced aircraft wings are systematically addressed. Moreover, the indicial aerodynamic concept is used to obtain the unsteady aerodynamic model. It is shown that for a fixed wing, compressibility has a significant influence on the determination of the classical flutter speed. Aeroelastic stability margins of aircraft wings at subsonic regime is studied via a unified aeroelastic model by Qin et al. (2002). Fazelzadeh and Hosseini (2007) applied Librescu’s theory to the problem of aerothermoelastic behavior of TWBs in supersonic gas flows. They obtained the quasi steady aerodynamic loading based on the first order piston theory. Shadmehri et al. (2007) discussed the flexural–torsional vibrations of laminated

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

3

TWBs with single-cell closed cross section. They have also investigated the aeroelastic stability of composite wings in an incompressible flow in a separate study (Haddadpour et al., 2008) where the aerodynamic model is constructed based on the Wagner’s function. Na et al. (2011) performed aeroelastic analysis of composite TWBs in compressible fluid regime. Effects of geometrical and material coupling and flight Mach number on the aeroelastic stability margins are studied in depth for laminated TWBs in subsonic regime. Sodja et al. (2016) conducted a dynamic aeroelastic wind tunnel experiment on the aeroelasticity of composite wings. The aeroelastic behavior is studied in terms of the critical flutter speed, LCO frequency and the tip deflection and torsion. In a study performed by Rajpal and De Breuker (2017), a reduced order aeroelastic model is employed to perform the aeroelastic and optimization analysis of composite wings subjected to gust loads. Stodieck et al. (2013, 2014) used tow steered composites to tailor the aeroelastic behavior and to improve the flutter characteristics of composite wings. Nonlinearities in wing structures are mainly due to geometric and material nonlinearities, nonlinear contacts and freeplays, while other factors such as unsteady aerodynamic loadings may also introduce nonlinearities to the problem (Kim, 2004). Dowell et al. (2003) provided a valuable reference on the effect of various nonlinearities in both structural and aerodynamic models on the aeroelastic response. Price et al. (1994) studied the nonlinear aeroelasticity of a twodimensional airfoil subjected to incompressible flow with a structural free-play nonlinearity in pitch degree of freedom, using power spectral density, phase portrait and Poincare plots. Tang and Dowell (2004) coupled the nonlinear beam and the nonlinear aerodynamic stall models to study the effects of nonlinear terms on the flutter and the post flutter (LCO) behavior of high aspect ratio wings. Patil et al. (2001) considered nonlinearity terms both in the aerodynamic and structural models in an aeroelastic model to investigate the aeroelastic response of high-aspect ratio wings. The doubletlattice aerodynamic theory as well as the nonlinear finite element solution methodology is used in their study. Eken and Kaya (2015) had a close look at limit cycle oscillations of cantilevered swept back wings with cubic nonlinearities in an incompressible fluid flow. In the current study, nonlinear aeroelasticity of composite wings is studied using a geometrically nonlinear formulation for the TWB. Geometrically nonlinear terms of up to third order are retained in the formulation. Geometrically nonlinear TWB theory for modeling the wing is a higher fidelity structural beam-wing model to be used in the aeroelastic system compared to the nonlinear Euler–Bernoulli or Timoshenko beam models, which are mostly used in the nonlinear aeroelastic studies in the literature. The structural equation of motion includes several non-classical effects such as warping inhibition and transverse shear strain. Using the aerodynamic states and the indicial compressible aerodynamics theory, explicit aerodynamic loadings involving lift and pitching moment are obtained in the time domain. The main aspect of this formulation is the Mach dependent exponential approximation of the indicial functions which makes it possible to perform direct stability analysis at any subsonic Mach number. The aerodynamic model is valid for thin airfoils. The ideal form of linear, attached, inviscid, compressible and unsteady aerodynamic flow based on the indicial theory is proposed by the present work. As a limitation for the present aeroelastic model, transonic regime including shock effects cannot be covered. The aerodynamic model is linear and the flow is assumed to be attached to the airfoil. The coupled field equations of motion are obtained by utilizing the Hamilton’s principle. Nonlinear aeroelastic responses are then obtained for composite wings with circumferentially asymmetric stiffness (CAS) configuration TWB structural model. In the solution procedure, the linear eigenfunctions are first determined by solving the linear structural equations of motion and then eigenfunctions are utilized as the base spatial solutions in the Ritz method for the solution of the nonlinear aeroelastic system of equations. According to the results obtained, qualitative and quantitative effects of the wing sweep angle, the fiber layup of the CAS configuration, on the LCO behavior are determined. It should be noted that to the best of authors’ knowledge, nonlinear aeroelasticity of wings based on geometrically nonlinear thin walled beam model by retaining nonlinear terms of up to third order and compressible unsteady aerodynamics based on indicial functions for the investigation of the post-critical aeroelastic response of composite wings has not been studied before. It is believed that present study provides results which can be used for comparison purposes in similar studies which utilize different geometrically nonlinear structural wing and compressible unsteady aerodynamic models in an aeroelastic framework. 2. Governing equations 2.1. Basic assumptions and kinematic relations Consider a swept aircraft wing composed of a single-cell thin walled beam, as shown in Fig. 1. Note that the center of twist and aerodynamic axes are not coincident and separated by the distance ab, with a being the dimensionless coefficient. Mathematically speaking, the wing box structure of the present study may be modeled as a thin walled beam which is described in Fig. 2 for a circumferentially asymmetric stiffness (CAS) configuration. In Fig. 2 L, l, h, d and b are the length, width, thickness, height and the semi-chord, respectively. The width and the height of the TWB at the root and the tip are denoted by lroot , droot and ltip , dtip , respectively. It should be noted that, Fig. 2 is a general model of TWB, however, in the present study, untapered TWB model is investigated. As depicted in Figs. 1 and 3, there is a need to define two different coordinate systems. A Cartesian fixed coordinate system (x, y, z) is located at the root as the main reference while the local coordinate system (n, s, z) is at the mid plane of the cross section of the TWB. Parameter s defines the local tangent to the middle surface and n is the axis perpendicular to the s coordinate.

4

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

Fig. 1. Schematic description of a swept aircraft wing structurally modeled as a doubly tapered TWB (a) Planform view (b) Wing’s detailed view (c) Cross-section view.

Fig. 2. Schematic demonstration of the CAS layup configuration and the associated parameter definitions.

The derivation of the governing equations of motion, basic assumptions considered in the present study follow those of Bhaskar and Librescu (1995). The considered structural model is similar to that developed by Farsadi et al. (2017, 2018). For a detailed description of the original structural model, the reader is referred to these references. Motion kinematics may be defined according to the six degrees of freedom available in the presented model. The displacement field for the TWB may be defined in terms of u(x, y, z , t), v (x, y, z , t), w (x, y, z , t) which are then described in terms of translations of the shear center of the TWB along the x, y and z axes (u0 (z , t), v0 (z , t), w0 (z , t)) and the three rotational degrees of freedom (θx (z , t), θy (z , t), φ (z , t)) about the x, y and z axes, Librescu and Song (2006), u(x, y, z , t) = u0 (z , t) − (y − n dx ) sin φ (z , t) − (x + n ds )(1 − cos φ (z , t)) ds dy

v (x, y, z , t) = v0 (z , t) + (x + n dy ) sin φ (z , t) − (y − n dx )(1 − cos φ (z , t)) ds ds w(x, y, z , t) = w0 (z , t) + (x +

θ

dy n ds ) y (z

(1)

, t) + (y − n dx )θ (z , t) − [Fw (s) + na(s)]φ ′ (z , t) ds x

where the over-prime symbol represents differentiation with respect to the z coordinate. Moreover, Fw (s), na(s) are the primary and secondary warping functions (Farsadi et al., 2017).

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

5

Fig. 3. Cross sectional view of the TWB and the associated variable definitions.

Strain–displacement relations are delineated according to the Green–Lagrange strain field definitions. Excluding the zero strain components, the remaining components in the local (n, s, z) coordinate system are then expressed as,

∂w 1 ∂u ∂v + [( )2 + ( )2 ] ∂z 2 ∂z ∂z

εzz =

γsz = γxz γnz = γxz

dx

+ γzy

ds dy

− γzy

ds

dy ds dx

(2)

+ Ψ (s)φ ′ + 2nφ ′

(3) (4)

ds

where definition for Ψ (s) is given in Farsadi et al. (2017). 2.2. Energy expressions Governing equations of motion are determined using the Hamilton’s principle. Based on the Hamiltonian dynamics, Hamilton’s principle is given in Eq. (5),

δH =

t2

∫

δ (T − V + Wn.c ) dt = 0

(5)

t1

where T , V are the total kinetic and strain potential energies, respectively and Wn.c stands for the work done by the external non-conservative loads. The energy stored in the structure as the strain potential is also defined in terms of non-zero strain components as, V =

1

∫ L∮ ∑ Nl ∫

2

0

C k=1

nk

[σzz εzz + σsz εsz + σnz εnz ](k) dn ds dz

(6)

nk−1

where the linear constitutive equation relating stresses to strains is given in many research articles, e.g. Farsadi et al. (2017, 2018) or Librescu and Song (2006). Substituting the strain expressions given by Eqs. (2)–(4) into Eq. (6) in the same manner as done in Farsadi et al. (2017, 2018), integrating over the thickness of the TWB and exploiting the force/moment resultant definitions given in Farsadi et al. (2017, 2018), the total strain energy is expressed as,

⎡ V =

1 2

L

∫ 0

⎤

Tz (w ′ 0 + 12 (u′ 0 + v ′ 0 )) + Qy (θx + v ′ 0 cos φ − u′ 0 sin φ ) 2

2

⎢ ⎥ ⎣+Qx (θy + u′ 0 cos φ + v ′ 0 sin φ ) + My (θ ′ y − u′ 0 φ ′ sin φ + v ′ 0 φ ′ cos φ )⎦ dz +Mx (θ x − u 0 φ cos φ − v 0 φ sin φ ) + Mz φ − Bw φ + ′

′

′

′

′

′

′′

1 Λ 2 z

φ

(7)

′2

where Tz (axial force), Qx (chord-wise shear force), Qy (flap-wise shear force), Mx (flap-wise bending moment), My (chordwise bending moment), Mz (Saint-Venant twist moment) and Bw (bi-moment or warping torque) are explicitly defined in the previous study of the authors (Farsadi et al., 2017).

6

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

Force and moment resultants may be determined explicitly in terms of the six degree of freedom as given in Eq. (8) which are also reported by Farsadi et al. (2017, 2018).

⎧ ⎫ Tz ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ Qx ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ Q ⎪ ⎪ ⎨ y⎪ ⎬

⎡

a11

⎢ ⎢ ⎢ ⎢ My =⎢ ⎢ Mx ⎪ ⎪ ⎪ ⎪ ⎢ ⎪ ⎪ ⎪ ⎢ Mz ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎣ ⎪ ⎩ Bw ⎪ ⎭ Λz sym

a12 a22

a13 a23 a33

a14 a24 a34 a44

a15 a25 a35 a45 a55

a16 a26 a36 a46 a56 a66

⎤⎧

⎫

a18 ⎪ w′ + 0.5(u′ + v ′ 0 ) ⎪ ⎪ ⎪ ⎪ θ +0 u′ cos φ 0+ v ′ sin ⎪ a28 ⎥ ⎪ φ ⎪ ⎪ ⎪ y 0 0 ⎪ ⎪ θ + v ′ cos φ − u′ sin φ ⎪ ⎪ ⎪ a38 ⎥ x 0 0 ⎥⎪ ⎪ ⎪ ⎨ ⎬ ′ ′ ′ ′ ′ ⎥ a48 ⎥ θ y − u 0 φ sin φ + v 0 φ cos φ ′ ′ ′ ′ ′ a58 ⎥ ⎪θ x − u 0 φ cos φ − v 0 φ sin φ ⎪ ⎥⎪ ⎪ ⎪ ⎪ ⎪ a68 ⎥ ⎪ φ′ ⎪ ⎪ ⎪ ⎦⎪ ′′ ⎪ a78 ⎪ −φ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ 2 a88 0.5φ ′

a17 a27 a37 a47 a57 a67 a77

2

2

(8)

The geometric nonlinearity incorporated in the current study makes the stiffness matrix in Eq. (8), fully populated. Explicit definition for the (aij ) coefficients in the stiffness matrix are reported in the previous study of the authors (Farsadi et al., 2017). The next important energy responsible for the inertial effects is the kinetic energy of the TWB (T ) which is defined as,

T =

1

∫ L∮ ∑ Nl ∫

2

0

C k=1

nk

ρs R˙ 2 dndsdz

(9)

nk−1

where over dot represents differentiation with respect to the temporal coordinate, ρs is the mass density of the TWB and, referring to Figs. 1–2, the position vector R of an arbitrary point on the TWB is defined by, R = (x + u)i + (y + v )j + (z + w )k

(10)

Using the deformed position of an arbitrary point (u, v, w ) in the kinetic energy definition, Eq. (9), its variation may be deduced as the following,

δT = −

L

∫

[I1 δ u0 + I2 δv0 + I3 δw0 + I4 δθx + I5 δθy + I6 δφ + I7 δφ ′ ]dz

(11)

0

where I1 = b1 u¨ 0 − b2 (φ¨ cos φ − φ˙ 2 sin φ ) − b3 (φ¨ sin φ + φ˙ 2 cos φ ) I2 = b1 v¨ 0 − b2 (φ¨ sin φ + φ˙ 2 cos φ ) + b3 (φ¨ cos φ − φ˙ 2 sin φ ) I3 = b 1 w ¨ 0 + b2 θ¨x + b3 θ¨y − b7 φ¨ ′ I4 = (b4 + b12 )θ¨x + (b6 − b13 )θ¨y + b2 w ¨ 0 − (b8 − b14 )φ¨ ′

(12)

I5 = b 3 w ¨ 0 + (b5 + b11 )θ¨y + (b6 − b13 )θ¨x − (b9 + b15 )φ¨ ′ I6 = b3 (v¨ 0 cos φ − u¨ 0 sin φ ) − b2 (u¨ 0 cos φ + v¨ 0 sin φ ) + (b4 + b5 + b12 + b11 )φ¨ I7 = b 7 w ¨ 0 + (b8 − b14 )θ¨x + (b9 + b15 )θ¨y − (b10 + b16 )φ¨ ′ and

{b1 , b2 , b3 , b4 , b5 , b6 , b7 , b8 , b9 , b10 } =

∮

{b11 , b12 , b13 , b14 , b15 , b16 } =

∮

C

m0 {1, y, x, y2 , x2 , xy, Fw (s), yFw (s), xFw (s), Fw2 (s)}ds m2 {(

C Nl

{m0 , m2 } =

∑∫ k=1

dy ds

nk

2

) ,(

dx ds

2

) ,

dy dx ds ds

, a(s)

ρs(k) {1, n2 }dn

dx ds

, a(s)

dy ds

, a2 (s)}ds (13)

nk−1

2.3. Unsteady aerodynamic loading in the subsonic compressible flow regime Work done by the non-conservative external aerodynamic loading is the last portion required to complete the Hamilton’s relation, Eq. (5). Variation of the non-conservative work done on the TWB is given as,

δ Wn.c =

∫

L

[Lae (z , t)δv0 (z) + Mae (z , t)δφ (z)] dz

(14)

0

where Lae and Mae are the unsteady aerodynamic lift and pitching moment about the reference axis, respectively. In the present study, explicit expressions for the unsteady aerodynamic loading in the subsonic compressible flow regime in time domain are obtained using indicial aerodynamics which is well accepted by the scientific community (Tobak, 1954; Lomax, 1961; Leishman, 2006; Bisplinghoff et al., 2013). Indicial response is a mathematical concept which

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

7

cannot be obtained directly from experiments. For incompressible and inviscid flows, closed form solutions are available for indicial responses (Bisplinghoff et al., 2013) while for compressible subsonic flows, there is no direct analytical solution, even though approximate solutions may be obtained from inversions of periodic aerodynamic responses in the frequency domain (Leishman, 1993) or analytical generalization of the incompressible indicial functions to obtain compressible indicial functions (Ronch et al., 2018; Berci and Righi, 2017). In the incompressible unsteady aerodynamics, the aerodynamic loadings originate from the two main sources which are circulatory and non-circulatory. Non-circulatory lift and moment depend on the instantaneous accelerations and velocities of the wing. However, in a compressible medium where the speed of sound is finite, the non-circulatory flow does not adapt itself to the changes of the boundaries instantly, hence non-circulatory lift and moment do not only depend on the instantaneous values of acceleration and velocity but also on the history of motion. For compressible flows, it is not customary to make distinction between the circulatory and the non-circulatory flow since the history of the motion becomes important and apparent mass is no longer a meaningful concept (Bisplinghoff et al., 2013; Marzocca et al., 2002). In the present study, indicial functions are employed to express the unsteady aerodynamic loading in the compressible flow regime. For arbitrary small motions of the thin airfoil in subsonic flow, with respect to the reference axis placed at the leading edge of the airfoil, downwash velocity corresponding to pitching and plunging motions can be expressed as Bisplinghoff et al. (2013),

wa (x, z , t) = w ˆ a (z , t) − xφˆ a (z , t)

(15)

where the plunging and the pitching motions w ˆ a (z , t) and φˆ a (z , t) are given by Eqs. (16)–(17).

w ˆ a (z , t) = v˙ 0 − U cos Λφ + φˆ a (z , t) = φ˙ +

∂φ U sin Λ ∂z

∂v0 U sin Λ ∂z

(16) (17)

It should be noted that Eqs. (16)–(17) are written for a swept wing configuration with the wing sweep angle Λ. Also, for the sake of brevity the free stream fluid velocity U∞ is represented by U hereinafter. Eqs. (15)–(17) show that w ˆ a (z , t) term gives constant downwash, whereas the pitching motion φˆ a (z , t) induces a linear variation of downwash with x measured from the leading edge. Following the detailed derivations presented in Appendix A, unsteady aerodynamic lift Lae (z , t) and pitching moment Mae (z , t) about the pitch axis located at a distance b(a + 1) behind the leading edge are expressed as,

[

Lae (z , t) = 2CLφ ρ U cos Λb2 φˆ a (z , 0)φ¯ cq (t) + D2 (z , t)

]

[ ] − CLφ ρ U cos Λb w ˆ a (z , 0)φ¯ c (t) + D1 (z , t) [ ] Mae (z , t) = 4CLφ ρ U cos Λb3 φˆ a (z , 0)φ¯ cMq (t) + D4 (z , t) [ ] ˆ a (z , 0)φ¯ cM (t) + D3 (z , t) − 2CLφ ρ U cos Λb2 w

(18)

(19)

where Di ; i = 1..4 are defined as, D1 (z , t) =

t

∫

dw ˆ a (z , τ )

D2 (z , t) =

t

∫

dφˆ a (z , τ ) dτ

0

D3 (z , t) =

t

∫

D4 (z , t) =

(20)

φ¯ cq (t − τ ) dτ

(21)

dw ˆ a (z , τ ) dτ

0 t

∫

φ¯ c (t − τ ) dτ

dτ

0

dφˆ a (z , τ ) dτ

0

φ¯ cM (t − τ )dτ

(22)

φ¯ cMq (t − τ ) dτ

(23)

In order to include the 3-D effects of the finite span wing, lift curve slope CLφ is given by Diederich (1951), 2π cos Λe

CLφ = √ 1−

M2

√

1+

Λe 2 1 ( 2 cos ) AR 1−M 2

+

2 cos Λe AR

,

(24)

where AR is the aspect ratio and

√ cos Λe = √

1 − M2

1 − M 2 cos2 Λ

cos Λ.

(25)

8

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

2.4. Parametric approximation of indicial functions in subsonic potential flow Exponential representations of the Mach number dependent indicial functions, defined with respect to the leading edge of the airfoil, may be utilized to obtain the state space form of the aerodynamic loads (Eq. (26)). In the literature, coefficients of the exponential representations of the indicial functions are obtained numerically for a limited number of Mach numbers (M = 0.5, 0.6, 0.7, 0.8) (Bisplinghoff et al., 2013; Mazelsky and Drischler, 1952; Mazelsky, 1952; Leishman, 1988). However, in an aeroelastic analysis usually iterative methods are used and for an airspeed, damping of the system is traced. To take the compressibility effects into account properly, compressible indicial functions have to be evaluated for the Mach number which corresponds to the input airspeed. However, since the compressible indicial functions are known for discrete Mach numbers in the literature, an appropriate interpolation method is required to evaluate the compressible indicial functions for any Mach number.

[ φc ,cq,cM ,cMq (M , τ ) = α0c ,0cq,0cM ,0cMq (M) −

3 ∑

] −βi τ

αic ,icq,icM ,icMq (M)e

(26)

i=1

In the following, a methodology is presented to extract the coefficients of the exponential representations of the indicial functions for any Mach number less than 0.8 in the compressible flow regime. As shown in Eq. (26), it is assumed that all plunging and pitching indicial functions are established in terms of summation of exponential functions with the same Mach-independent power coefficients (βi , i = 1, 2, 3). It is further assumed that for all exponential representations of the indicial functions, Mach-independent power coefficients βi are equal to their counterparts in the exponential approximation of the plunging indicial function for the lift for M = 0.5 and they are given by β1 = 0.0754 , β2 = 0.372 , β3 = 1.89 (Mazelsky and Drischler, 1952). Thus, for each indicial function in Eq. (26), there are only four Mach-dependent base coefficients αic ,icM ,icq,icMq (M) , i = 0, 1, 2, 3 to be determined. The asymptotic values of the indicial functions (Eq. (26)) should recover their counterparts in the steady flow, and √ they are obtained by multiplying the steady state incompressible indicial functions by the Prandtl–Glauert factor 1/ 1 − M 2 (Bisplinghoff et al., 2013)

φc (∞) = √ φcM (∞) = φcq (∞) = φcMq (∞) =

1 1 − M2 −1

,

(27)

√

,

(28)

√

,

(29)

√

.

(30)

4 1 − M2 3 4 1 − M2 −1 4 1 − M2

Utilizing Eqs. (27)–(30), the first base coefficients of the plunging and pitching indicial functions for the lift and the pitching moment (α0c ,0cq,0cM ,0cMq (M)) are found. It should be noted that asymptotic values of the indicial functions, which are equal to the first base coefficients, are independent of the power coefficients (βi , i = 1, 2, 3). For the determination of the remaining three base coefficients for each indicial function (αic ,icM ,icq,icMq (M), i = 1, 2, 3), three conditions have to be imposed. In the time interval 0 ≤ τ ≤ M2M , Lomax et al. (1952) have obtained exact expressions for the dimensionless +1 indicial functions given by, 2 [ 1− πM [ −1 φcM (τ ) = 1− πM

φc (τ ) =

φcq (τ ) =

1

[ 1−

τ 2M

τ 2M

τ

]

(1 − M) , (1 − M) + (1 − M) +

(31)

τ2 8M

τ2

]

(M − 2) , (1 −

M

(32)

]

) ,

πM 2M 4M 2 [ ]} { −2 3τ 3τ 2 τ3 1 2 3 φcMq (τ ) = 1− (1 − M) + (1 − M) + M + (1 − M) . 3π M 4M 32M 2 16M 3 4

(33) (34)

where τ = Ut /b. Further descriptions related to determining the aerodynamic indicial functions for compressible flow in the Theodorsen’s coordinate are given in Appendix A. Four conditions are imposed to obtain the four unknown coefficients in Eq. (26) for each indicial function. The first condition is based on the fact that asymptotic values of the indicial functions should recover their counterparts in the steady flow. Two of the conditions are obtained by evaluating the indicial functions at the non-dimensional times of τ = 0 and τ = M2M . The last condition is obtained by fitting polynomials to the already known exponential representations of +1

the indicial functions at Mach numbers 0.5, 0.6, 0.7 and 0.8 (Mazelsky and Drischler, 1952; Mazelsky, 1952; Leishman,

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

9

Fig. 4. Comparison of the indicial lift functions for the plunging airfoil at (a) M = 0.5 (b) M = 0.6.

1988; Farsadi and Javanshir, 2012; Bisplinghoff et al., 2013). It should be noted that for the plunging indicial function for the lift (φc ), exponential representation of the indicial function is known at Mach numbers 0.5, 0.6, 0.7 and 0.8 whereas for the remaining indicial functions exponential representations are known at Mach numbers 0.5, 0.6 and 0.7. In this respect, plunging indicial function for lift is approximated by a third order polynomial given by Eq. (35).

φ˜ c (τ ) = ac3 τ 3 + ac2 τ 2 + ac1 τ + ac0

(35)

The known exponential representations of the plunging indicial function for the lift (φc ) for Mach numbers 0.5, 0.6, 0.7 and 0.8 are evaluated at four different non-dimensional times (Ut1 /b = 4 (M = 0.5); Ut2 /b = 4.8(M = 0.6); Ut3 /b = 5.6(M = 0.7); Ut4 /b = 6.4(M = 0.8)) and the four unknown coefficients (ac0 , ac1 , ac2 , ac3 ) of the polynomial approximation of the indicial function φc are determined. The same process is repeated for the other indicial functions but second order polynomials are fitted since the exponential representations for the indicial functions φcq,cM ,cMq (M , τ ) are available for Mach numbers 0.5, 0.6 and 0.7. Third equation is then obtained by equating the polynomial approximation of the indicial function Eq. (35) evaluated at the non-dimensional time 40 to the relevant exponential approximation of the indicial function given by Eq. (26) at the same non-dimensional time. As shown later in the article, all indicial functions level out with the increase of the non-dimensional time. Hence, for all compressible indicial functions non-dimensional time 40 is selected to represent the behavior at a higher non-dimensional time. With this approach, unsteady compressible indicial functions can be evaluated at any Mach number. This novel method of evaluating the indicial functions at any Mach number and dimensionless time is very crucial for the aeroelastic stability analysis in the compressible flow regime which is described in the next section. Fig. 4 compares the indicial lift functions for the plunging airfoil approximated by the methodology described and the already available exponential representation of the indicial function. Fig. 4a and b are obtained at the Mach numbers M = 0.5 and M = 0.6, respectively. Once the Mach dependent base coefficients (αic ,icM ,icq,icMq (M) , i = 0, 1, 2, 3) of the lift and the moment compressible indicial functions for the plunging and pitching motion with respect to the coordinate system established at the leading edge of the airfoil are determined for any Mach number, base coefficients with respect to the axis located at b(a + 1) distance behind the leading edge can be determined utilizing the relations given by Eqs. (36)–(39). It should be noted that these relations are the same as those given by Eqs. (100)–(103) for the indicial functions themselves owing to the fact that, the approximate exponential representation of the four indicial functions are defined with the same Mach-independent power coefficients (βi , i = 1, 2, 3).

α¯ ic = αic

(36) a

1

α¯ icM = αicM + ( + )αic

(37) 2 1 α¯ icq = αicq − ( + )αic (38) 2 2 a 1 a 1 α¯ icMq = αicMq + ( + )(αicq − αicM ) − ( + )2 αic (39) 2 2 2 2 where i = 1, 2, 3. Having determined the base coefficients (αic ,icM ,icq,icMq (M) , i = 0, 1, 2, 3) with respect to the axis located at b(a + 1) distance behind the leading edge, indicial functions with respect to the Theodorsen’s coordinate are defined by Eq. (40). 2 a

[ φ¯ c ,cq,cM ,cMq (M , τ ) = α¯ 0c ,0cq,0cM ,0cMq (M) −

3 ∑ i=1

] −βi τ

α¯ ic ,icq,icM ,icMq (M)e

(40)

10

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

The integral expressions Eqs. (20)–(23), which appear in the lift and moment expressions in Eqs. (18)–(19), are expressed again by substituting the exponential representations of the indicial functions (Eq. (40)) into Eqs. (20)–(23). For instance, the integral D1 (z , t) in Eq. (20) is expressed as, D1 (z , t) =

t

∫

dw ˆ a (z , τ )

[ α¯ oc −

dτ

0

3 ∑

] α¯ ic (M)e−βi τ dτ

(41)

i=1

The first term in Eq. (41) gives the downwash multiplied by the asymptotic value of the plunging indicial function for the lift. By defining the integrals involving the exponential terms in Eq. (41) by the aerodynamic lag terms Bic (z , t) and making use of the Leibniz integral rule, Eq. (41) can be re-written as, D1 (z , t) = α¯ 0c w ˆ a (z , t) −

3 ∑

α¯ ic (M)Bic (z , t)

(42)

i=1

Similarly, integrals for D2 (z , t), D3 (z , t) and D4 (z , t) in Eqs. (21)–(23) are expressed as, D2 (z , t) = α¯ 0cq φˆ a (z , t) −

3 ∑

α¯ icq (M)Bicq (z , t)

(43)

i=1

D3 (z , t) = α¯ 0cM w ˆ a (z , t) −

3 ∑

α¯ icM (M)BicM (z , t)

(44)

α¯ icMq (M)BicMq (z , t)

(45)

i=1

D4 (z , t) = α¯ 0cMq φˆ a (z , t) −

3 ∑ i=1

It should be noted that in transforming Eq. (41) into Eq. (42), it is assumed that the wing starts the motion at rest, hence w ˆ a (z , 0) = φˆ a (z , 0) = 0. When the Leibniz integral rule is applied to the integrals involving the exponential terms in Eq. (41) , it can be shown that the aerodynamic lag terms Bic (z , t) are defined by Eq. (46). B˙ ic (z , t) + βi

U cos Λ b

Bic (z , t) =

dw ˆ a (z , t)

; i = 1, 2, 3

dt

(46)

Similarly, the other aerodynamic lag terms are defined by Eqs. (47)–(49). B˙ icq (z , t) + βi B˙ icM (z , t) + βi B˙ icMq (z , t) + βi

U cos Λ

b U cos Λ

b U cos Λ b

Bicq (z , t) =

BicM (z , t) =

BicMq (z , t) =

dφˆ a (z , t) dt dw ˆ a (z , t) dt dφˆ a (z , t) dt

; i = 1, 2, 3

(47)

; i = 1, 2, 3

(48)

; i = 1, 2, 3

(49)

Finally, assuming that the wing is at rest initially, the explicit expressions of the unsteady lift and moment in compressible flow regime are given by Eqs. (50) and (51).

[ Lae (z , t) = 2CLφ ρ U cos Λb

2

α¯ 0cq φˆ a (z , t) −

3 ∑

] α¯ icq Bicq (z , t)

(50)

i=1

[ − CLφ ρ U cos Λb α¯ 0c w ˆ a (z , t) −

3 ∑

] α¯ ic Bic (z , t)

i=1

[ Mae (z , t) = 4CLφ ρ U cos Λb

3

α¯ 0cMq φˆ a (z , t) −

3 ∑

] α¯ icMq BicMq (z , t)

(51)

i=1

[ − 2CLφ ρ U cos Λb

2

α¯ 0cM w ˆ a (z , t) −

3 ∑

] α¯ icM BicM (z , t)

i=1

2.5. Nonlinear aeroelastic integral equation of motion The final governing equations of motion are acquired in their integral representation after the substitution of the potential and kinetic energy variations as well as the non-nonconservative work done by the unsteady aerodynamic loading into the Hamilton’s principle expression in Eq. (5). Explicit definitions for the unsteady lift and moment loadings

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

11

are given by Eqs. (50) and (51). Accordingly, to fully define the time domain equations of motion, the aerodynamic lag states Bic ,icM ,icq,icMq (z , t) have to be also defined by augmenting the system of equations by Eqs. (46)–(49). To keep the formulation more concise, the integral equation of motion is reported here for the case of circumferentially asymmetric stiffness (CAS) layup configuration, and the respective nonlinear solution for the aeroelastic problem is presented. As partly depicted in Fig. 2 for the case of CAS layup configuration, the fiber angles on the outermost and innermost layers on the top and bottom flanges have opposite signs (θ (y) = −θ (−y); θ (x) = −θ (−x)). In the CAS layup configuration, due to the flapwise bending-torsion and extension-chordwise bending couplings some of the stiffness coefficients in Eq. (8) vanish and consequently the general nonlinear system of equations simplifies to some extent. The non-vanishing stiffness coefficients in Eq. (8) include a12 , a18 , a28 , a56 and a37 . It is worth mentioning that, the stiffness coefficient a56 is the flapwise bending-torsion coupling stiffness which defines the amount of load mitigation in flexible wing structures. Furthermore, the non-zero mass/inertia terms in the CAS layup configuration include b1 , b4 , b5 , b10 , b11 and b12 , where their definitions are given by Farsadi et al. (2018). Retaining both square and cubic nonlinear terms in the equation of motion, an explicit expression is obtained after a rather cumbersome manipulation as given in Eq. (52).

⎡ L

∫ 0

⎤

f1 (∆) δ u0 + f2 (∆) δ u′ 0 + f3 (∆) δv0 + f4 (∆) δv ′ 0 +

⎢ ⎣f5 (∆) δw0 + f6 (∆) δw′ 0 + f7 (∆) δθx + f8 (∆) δθ ′ x + f9 (∆) δθy + f10 (∆) δθ

′

y

+ f11 (∆) δφ + f12 (∆) δφ + f13 (∆) δφ ′

⎥ ⎦ dz = 0

(52)

′′

where ∆ denotes each one of the degrees of freedom u0 , v0 , w0 , θx , θy , φ . The functions fi in Eq. (52) are given in Appendix B and the test functions in Eq. (52) are variations of the degrees of freedom (δ u0 , δv0 , δw0 , δθx , δθy , δφ ). 3. Method of solution Nonlinear solution algorithms are mainly based on iterative schemes which can be practically very computationally expensive especially when complicating effects such as fluid-structure interactions are also included (Dimitriadis, 2017). In order to minimize the computational costs, approximations of the initial solution to the nonlinear problem must be improved. One of the most efficient methods to have an efficient initial guess to the solution of the nonlinear problem is to use the linear eigensolutions. Specifically speaking, we may utilize the solution to the problem of linear TWBs with CAS layup configuration as the initial approximation. To achieve this objective, using a Ritz based solution methodology, solution of the linear structural equations are acquired and the mode shapes for the TWB are determined. These mode shapes are then utilized in the nonlinear aeroelastic analysis of the composite wing to construct a modal basis solution for the problem and minimize the computational costs. 3.1. Mode shape determination for the Cantilevered TWBs A Ritz based solution method is adopted to obtain the eigenfunctions of linear composite cantilevered TWBs with CAS layup. For cantilevered wing structures, it is necessary to employ trial functions capable of satisfying the essential boundary conditions in order to discretize the spatial coordinate via the expansions given in Eq. (53), N ∑

u0 (z , t) =

ψiu (z)ηui (t)

i=1 N ∑

v0 (z , t) =

ψiv (z)ηvi (t)

i=1 N

w0 (z , t) =

∑

ψiw (z)ηwi (t)

i=1

θx (z , t) =

N ∑

(53)

ψ

η

x i (z) xi (t)

i=1 N

θy (z , t) =

∑

ψiy (z)ηyi (t)

i=1

φ (z , t) =

N ∑

φ

ψi (z)ηφi (t)

i=1

where N is the required number of terms for an acceptable convergence. Spatial functions for all degrees of freedom are taken as identical as shown in Eq. (54). u,v,w,θx ,θy ,φ

ψi

(z) = z i ; i = 1...N

(54)

12

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

The first seven expansion terms (N = 7) are included in the expansion relations in Eq. (53) to obtain fully converged results. To obtain the mode shapes of the TWB, all nonlinear terms are first omitted in Eq. (52) and the expansions given in Eq. (53) are substituted back into the modified governing equation of motion. All test functions in Eq. (52) are chosen from the spatial functions given by Eq. (54) which vanish at the root of the wing (z = 0) and their respective derivatives are determined accordingly. Integrating Eq. (52) over the axis of the TWB, one obtains the following system of equations, Ms η¨ + Ks η = 0

(55)

where Ks , Ms are respectively the structural stiffness and mass matrices which are expressed in Farsadi et al. (2018) and the vector of state variables is defined as,

{ }T η6N ×1 = ηuT ηvT ηwT ηxT ηyT ηφT { } where ηuT = ηu1 (t) ηu2 (t) · · · ηuN (t) , etc.

(56)

The system of equations given by Eq. (55) can be recast into the state space representation using the following state variables as shown in Eq. (57). G12N ×1 = ηT

}T

η˙ T

{

(57)

Utilizing Eq. (57), the state space form is obtained as shown in Eq. (58)

˙ = VG G [ V=

0 1 −M− s Ks

]

I 0

(58)

¯ λt in Eq. (58) yielding, Standard form of an eigenvalue problem is obtained by introducing G = Ge ¯ = 0. (V − Iλ)G

(59)

¯ is the corresponding eigenvector. where λ is the eigenvalue and G 3.2. Nonlinear aeroelastic analysis of composite wings Modal solution is computationally efficient in nonlinear aeroelastic analysis of composite wings structurally modeled as TWBs. To construct the modal basis solutions, one may utilize the linear eigenfunctions determined in the foregoing section as the spatial approximation to all degrees of freedom. Outcome of the eigenvalue analysis for the dominant (first) m right eigenvectors corresponding to the translational and rotational deformations of the linear system is the reduced modal matrix R as given in Eq. (60). Note that the m dominant eigenvectors are those corresponding to the lowest m eigenvalues.

[

T

R6N ×m = Ru m×N

Rvm×N T

Rw m×N T

T

Rxm×N

yT

Rm×N

φT

Rm×N

]

(60)

A reduced order model is consequently established for any of the six degrees of freedom ∆(u0 , v0 , w0 , θx , θy , φ ) in terms of the relevant trial function ψ ∆ and R∆ which is composed of the dominant m right eigenvectors and the modal coordinate as, T

∆(z , t) = ψ ∆ R∆ ϑ (t)

(61) ∆

where ϑ (t) is the generalized modal coordinate vector of dimension m × 1, R is the reduced modal matrix of dimension N × m composed of dominant m right eigenvectors corresponding to any of the degrees of freedom ∆ and ψ ∆ is the vector of trial functions, of dimension N × 1, corresponding to any of the degrees of freedom ∆. In the nonlinear system, the test functions for all degrees of freedom are defined by pre-multiplying the trial functions T ψ ∆ by the reduced modal matrix L∆ which is composed of the dominant m left eigenvectors of the linear system. T

δ ∆(z) = L∆ ψ ∆

(62)

To acquire the reduced order system of nonlinear equations, the modal expansion in Eq. (61) and the test functions in Eq. (62) are substituted back into the governing equations of motion given in Eq. (52), resulting in the following nonlinear system of equations,

¯ NL (t) = 0 Mt ϑ¨ (t) + Ct ϑ˙ (t) + Kt ϑ (t) − Z(t) + H

(63)

where Mt , Ct and Kt are the reduced order mass, aerodynamic damping and stiffness matrices of dimension m × m, and ¯ NL (t) are the reduced order vectors of dimension m × 1 including the aerodynamic lag states and the structural Z(t) and H nonlinear terms, respectively.

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

13

In terms of the left and right eigenvectors (LT , R), reduced order mass, damping and stiffness matrices in Eq. (63) are defined as, Mt = LT (Ms + Mae )R Ct = LT Cae R

(64)

Kt = LT (Ks + Kae )R where the structural mass Ms and the stiffness Ks matrices are given in Farsadi et al. (2018) while the aerodynamic mass Mae , stiffness Kae and the damping Cae matrices of dimension 6N × 6N are defined in Appendix C. ¯ NL (t), which includes the nonlinear terms, originates from the integral equation of motion in In Eq. (63), the vector H Eq. (52). Following the substitution of Eq. (104) into Eq. (52), nonlinear terms are grouped together as shown in Eq. (65).

⎡ HNL (t) =

L

∫ 0

1

1

(Nu2 + Nu3 )δ u′ 0 + (Nv2 + Nv3 )δv ′ 0 + (Nw2 + Nw3 )δw ′ 0 + (Nx2 + Nx3 )δθx

⎢ ⎢+(N 22 + N 32 )δθ ′ + (N 21 + N 31 )δθ + (N 22 + N 32 )δθ ′ x y y x x y y y y ⎣ 31 32 33 21 22 23 ′ ′′ +(Nφ + Nφ )δφ + (Nφ + Nφ )δφ + (Nφ + Nφ )δφ

⎤ ⎥ ⎥ dz ⎦

(65)

¯ NL (t) in Eq. (63), square and cubic nonlinear terms N 2 , N 3 To determine the reduced order vector of nonlinear terms H defined by Eqs. (105), modal definitions of the degrees of freedom given by Eq. (61) and the variations of the degrees of freedom given by Eq. (62) are substituted into Eq. (65) resulting in the expression given by Eq. (66) after a rather complex manipulation.

⎡

v

w

Lu ψ ′ (N¯ u2 + N¯ u3 ) + Lv ψ ′ (N¯ v2 + N¯ v3 ) + Lw ψ ′ (N¯ w2 + N¯ w3 )

¯ NL (t) = H

T

T

u

T

⎤

⎥ ∫ L ⎢ xT x 2 1 ⎢+L ψ (N¯ x + N¯ x31 ) + LxT ψ ′ x (N¯ x22 + N¯ x32 ) + LyT ψ y (N¯ y21 + N¯ y31 ) ⎥ ⎥ dz ⎢ ⎢ yT ′ y ¯ 22 32 ⎥ 31 φ T ′ φ ¯ 22 32 φ T φ ¯ 21 ¯ ¯ ¯ ) ) + L ψ ( N + N + N ) + L ψ ( N + N 0 ⎣+L ψ (N y y φ ⎦ φ φ φ 3 3 φ T ′′ φ ¯ 23 +L ψ (Nφ + N¯ φ ) ¯j

(66)

T

In Eq. (66), the nonlinear coefficients Ni are defined in terms of the trial functions (ψ ∆ ), derivatives of the trial functions, reduced modal matrices (R∆ ) and the modal coordinates ϑi (t) and they are given explicitly in Appendix D. 3.3. Final state space form of aeroelastic equations The aerodynamic damping matrix Cae in Eq. (64) does not include the aerodynamic lag terms defined in Eqs. (50)–(51). The lag terms Bic , Bicq , BicM , BicMq are collected in the aerodynamic lag vector Z which is derived from the virtual work done by the lift and the moment due to the aerodynamic lag states Bic ,icM ,icq,icMq (z , t). The virtual work done by the unsteady lift and moment due to the aerodynamic lag states in the integral equation of motion (52) is given by,

∫ L{

lag

lag

}

lag ( Lae |lag c (z , t) + Lae |cq (z , t))δv0 + ( Mae |cM (z , t) + Mae |cMq (z , t))δφ dz

(67)

0

where the unsteady aerodynamic loads through the lag states are given by Eq. (68). Lae |lag ¯ 1c B1c (z , t) + α¯ 2c B2c (z , t) + α¯ 3c B3c (z , t) c (z , t) = −CLφ ρ U cos Λb α

[

]

2 Lae |lag α¯ 1cq B1cq (z , t) + α¯ 2cq B2cq (z , t) + α¯ 3cq B3cq (z , t) cq (z , t) = 2CLφ ρ U cos Λb

[

]

lag

Mae |cM (z , t) = −2CLφ ρ U cos Λb2 α¯ 1cM B1cM (z , t) + α¯ 2cM B2cM (z , t) + α¯ 3cM B3cM (z , t)

[

lag

]

Mae |cMq (z , t) = 4CLφ ρ U cos Λb3 α¯ 1cMq B1cMq (z , t) + α¯ 2cMq B2cMq (z , t) + α¯ 3cMq B3cMq (z , t)

[

(68)

]

Substituting the Eq. (68) into Eq. (67) and employing the definition of the variations δv0 and δφ given by Eq. (62) and performing manipulations, the aerodynamic vector of lag states Z is obtained as,

Z(t) = α¯ 1c I

[

α¯ 2c I

···

α¯ 2cMq I

⎧ ⎫ ˆ 1c (t) ⎪ ⎪ B ⎪ ⎪ ⎪ ⎪ ⎪ ˆ 2c (t) ⎪ ⎪ ⎪ B ⎪ ⎪ ⎨ ⎬ ] . .. α¯ 3cMq I ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪Bˆ 2cMq (t)⎪ ⎪ ⎪ ⎪ ⎩Bˆ ⎭ (t) 3

(69)

cMq

ˆ ic (t), Bˆ icq (t), Bˆ i (t), Bˆ i (t); i = 1, 2, 3 are the transformed aerodynamic where I is the identity matrix of order m×m and B cM cMq lag state vectors of order m × 1 defined by, ∫ L T ˆ ic (t) = B CLφ ρ U cos ΛbLv ψ v Bic (z , t)dz 0

14

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812 L

∫

ˆ icq (t) = B

2CLφ ρ U cos Λb2 Lv ψ v Bicq (z , t)dz , T

0

ˆ i (t) = B cM

L

∫ 0

ˆ i (t) = B cMq

2CLφ ρ U cos Λb2 Lφ ψ φ BicM (z , t)dz ,

L

∫ 0

T

(70)

4CLφ ρ U cos Λb3 Lφ ψ φ BicMq (z , t)dz , T

Moreover, including the aerodynamic lag states, there would be 12m number of additional unknowns in Eq. (63), and to perform the nonlinear aeroelastic analysis an additional 12m number of equations is required. Therefore, the nonlinear system of equations (Eq. (63)) is augmented by the differential equations which define the aerodynamic lag states (Eqs. (46)–(49)). In order to incorporate the auxiliary equations given by Eqs. (46)–(49) into the nonlinear system of equations in Eq. (63), the following manipulations are made.

• Both sides of Eq. (46), which defines the aerodynamic lag states in the lift expression due to the plunging motion, are multiplied by CLφ ρ U cos Λb and by the test function δv0 , • Both sides of the second equation in Eq. (47), which defines the aerodynamic lag states in the lift expression due to the pitching motion, are multiplied by 2CLφ ρ U cos Λb2 and by the test function δv0 , • Both sides of the third equation in Eq. (48), which defines the aerodynamic lag states in the moment expression due to the plunging motion, are multiplied by 2CLφ ρ U cos Λb2 and by the test function δφ , • Both sides of the fourth equation in Eq. (49), which defines the aerodynamic lag states in the moment expression due to the pitching motion, are multiplied by 4CLφ ρ U cos Λb3 and by the test function δφ , and for each case the resulting expressions on both sides are integrated along the length of the wing. As the outcome of the above manipulations, Eqs. (46)–(49) are transformed into Eq. (71)–(74). L

∫ 0

L

∫ 0

L

∫ 0

L

∫ 0

CLφ ρ U cos ΛbLv ψ v [B˙ ic (z , t) + βi T

U cos Λ b

Bic (z , t)]dz = D1c ϑ˙ (t) + D2c ϑ¨ (t)

U cos Λ

2CLφ ρ U cos Λb2 Lv ψ v [B˙ icq (z , t) + βi T

b

Bicq (z , t)]dz = D1cq ϑ˙ (t) + D2cq ϑ¨ (t)

U cos Λ

2CLφ ρ U cos Λb2 Lφ ψ φ [B˙ icM (z , t) + βi T

b

4CLφ ρ U cos Λb3 Lφ ψ φ [B˙ icMq (z , t) + βi T

BicM (z , t)]dz = D1cM ϑ˙ (t) + D2cM ϑ¨ (t)

U cos Λ b

BicMq (z , t)]dz = D1cMq ϑ˙ (t) + D2cMq ϑ¨ (t)

(71) (72) (73) (74)

where i = 1, 2, 3 and D1c , D2c , D1cq , D2cq , D1cM , D2cM , D1cMq , D2cMq are the m × m coefficient matrices defined by, L

∫ D1c = 0

L

∫ D2c = − 0 L

∫ D1cq = 0

L

∫ D2cq = 0

L

∫ D1cM = 0

0 L

∫ D2cMq = 0

T

T

φ

T

φT

Rφ ) dz

φ

T

T

Rφ ) dz

(77)

Rφ ) dz

T

T

L

(78) T

T

T

2CLφ ρ U cos Λb2 (Lφ ψ φ ψ v Rv ) dz T

2CLφ ρ U 2 b3 sin 2Λ (Lφ ψ φ ψ ′ T

4CLφ ρ U cos Λb3 (Lφ ψ φ ψ ′

φ

T

φ

T

Rφ ) dz

Rφ ) dz

(75) (76)

vT

2CLφ ρ (U cos Λ)2 b2 (Lφ ψ φ ψ φ Rφ − tan Λ Lφ ψ φ ψ ′ Rv ) dz

0 L

T

T

CLφ ρ U cos Λb (Lv ψ v ψ v Rv ) dz

2CLφ ρ U cos Λb2 (Lv ψ v ψ ′

∫

D1cMq =

T

CLφ ρ U 2 b2 sin 2Λ(Lv ψ v ψ ′

D2cM = −

∫

T

CLφ ρ (U cos Λ)2 b(Lv ψ v ψ φ Rφ − tan ΛLv ψ v ψ ′

(79) (80) (81) (82)

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

15

Comparing Eq. (70) and Eqs. (71)–(74) and noting that the coefficients multiplying the aerodynamic lag terms on the left hand side of Eqs. (71)–(74) are independent of time, one may recast Eqs. (71)–(74) as,

ˆ ic (t) dB

U cos Λ

ˆ ic (t) = D1c ϑ˙ (t) + D2c ϑ¨ (t) B b U cos Λ ˆ icq (t) = D1cq ϑ˙ (t) + D2cq ϑ¨ (t) + βi B dt b ˆ icM (t) dB U cos Λ ˆ icM (t) = D1cM ϑ˙ (t) + D2cM ϑ¨ (t) + βi B dt b ˆ icMq (t) U cos Λ dB ˆ icMq (t) = D1cMq ϑ˙ (t) + D2cMq ϑ¨ (t) + βi B dt b dt ˆ icq (t) dB

+ βi

(83) (84) (85) (86)

ˆ ic , Bˆ icq , Bˆ icM , Bˆ icMq ; (i = 1, 2, 3) are the transformed aerodynamic lag state vectors of dimension m × 1 defined where B by Eq. (70). It should be noted that to obtain Eqs. (83)–(86) from Eqs. (71)(74), it is assumed that semi-chord length b is constant along the wing. For tapered wings, one may use the mean semi-chord length only in the U cos Λ/b term in Eqs. (83)–(86). Gathering Eq. (63) and Eqs. (83)–(86) together, the final consistent system of nonlinear aeroelastic equations for the composite swept wings structurally modeled as TWB with the CAS layup is obtained. Eq. (87) gives the state space representation of the final expression of the nonlinear aeroelastic system of equations,

[

DCM Dd

]

0 −I

{ } ϑ˜ ˜ dt B

[

d

14m×14m

14m×1

D + K 0

Dα Dβ

] 14m×14m

{ } ϑ˜ ˜ B

14m×1

{

¯ NL H + 0

} =0

(87)

14m×1

where

⎧ ⎫ ˆ 1c ⎪ B ⎪ { } [ ] [ ] ⎨ ⎬ { } { } ϑ Ct Mt Kt 0 . ˜ ˜ . ϑ = ˙ ; B = ; DCM = ; DK = . ⎪ I 0 0 −I ϑ 2m×1 ⎪ ⎩ˆ ⎭ 2m×2m 2m×2m B3cMq 12m×1 ⎤ ⎡ ¯β 0 ⎤ ⎡ D 0 0 β1 Im×m 0 0 ⎥ ⎢ 0 D¯ β 0 0 Λ⎣ ⎥ ; D¯ β = − U cos Dβ = ⎢ 0 β2 Im×m 0 ⎦ b ⎣0 ¯β 0 ⎦ 0 D 0 0 β3 Im×m 3m×3m ¯β 0 0 0 D 12m×12m ] [[ ] α1c Im×m · · · α3cMq Im×m m×12m Dα = 0m×12m 2m×12m ⎤ ⎡⎡ ⎤ Im×m

. ⎥ ⎢⎢ ⎢⎣ .. ⎦ ⎢ ⎢ Im×m 3m×m ⎢⎡ ⎢ Im×m ⎤ ⎢ ⎢⎢ . ⎥ ⎢⎣ .. ⎦ ⎢ ⎢ I ⎢ m×m ⎤3m×m Dd = ⎢⎡ ⎢ Im×m ⎢⎢ . ⎥ ⎢⎣ . ⎦ ⎢ . ⎢ ⎢ Im×m 3m×m ⎤ ⎢⎡ ⎢ Im×m ⎢ ⎢⎢ .. ⎥ ⎣⎣ . ⎦ Im×m

.

[

.

[

.

[

.

[

D1c

D1cq

D1cM

D1cMq

3m×m

D2c

]

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ] ⎥ D2cq m×2m ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ] D2cM m×2m ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ] ⎥ D2cMq m×2m ⎦

(88)

m×2m

12m×2m

Aeroelastic response of the wing is obtained in the time domain for a prescribed set of initial conditions by direct integration of Eq. (87) via the Runge–Kutta method. The time history plots are then utilized for nonlinear stability analysis.

4. Numerical results and discussions As the first case, numerical verification studies are conducted and then linear and nonlinear results are presented for different TWB configurations.

16

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

Table 1 Geometric and material properties of the wings used for the validation of the aeroelastic system based on geometrically linear structural model. Model

EI (N m2 )

GJ (N m2 )

L (m)

b (m)

a

Iα (Kg m2 )

ml (Kg/m)

Λ (deg)

ρ (Kg/m3 )

c (m/s)

Freon (%)

40A 40A 40D 73 72 25A 15 12 15 13 14 24

15.00 15.00 14.50 16.80 22.40 18.50 54.70 156.8 54.70 152.4 94.62 30.96

10.15 10.15 9.55 15.6 11.7 5.60 13.0 41.3 13.0 31.82 26.49 8.06

0.630 0.630 0.630 0.660 0.660 0.813 0.813 0.422 0.813 0.462 0.574 0.553

0.0509 0.0509 0.0509 0.0509 0.0509 0.0509 0.0509 0.0978 0.0509 0.0713 0.0713 0.0719

−0.20 −0.20 −0.21 −0.12 −0.12 −0.20 −0.08 −0.07 −0.08 −0.08 −0.08 −0.02

0.00024 0.00024 0.00023 0.00022 0.00024 0.00013 0.00017 0.00014 0.00017 0.00051 0.00054 0.00038

0.339 0.338 0.322 0.476 0.524 0.138 0.289 0.637 0.289 0.439 0.461 0.238

0 0 15 30 15 60 60 15 60 30 45 45

1.108 1.721 1.092 1.020 1.731 0.491 3.914 1.891 2.540 3.844 3.710 1.747

153.78 153.98 155.14 151.59 149.25 153.91 156.90 156.18 160.06 132.79 156.45 149.00

89 90 89 95 94 88 92 97 90 99 85 93

Table 2 Critical Mach numbers and frequencies of different wing models. Model

40A 40A 40D 73 72 25A 15 12 15 13 14 24

Flutter frequency (ωf , rad/s)

Flutter Mach number (Mf )

Flutter speed (m/s)

Experimental (Barmby et al., 1951)

Present

Experimental (Barmby et al., 1951)

Present

Experimental (Barmby et al., 1951)

Present

0.50 0.45 0.51 0.69 0.59 0.79 0.51 0.79 0.62 0.68 0.56 0.54

0.526 0.475 0.482 0.708 0.610 0.766 0.537 0.812 0.637 0.704 0.577 0.556

76.89 69.29 79.12 104.60 88.06 121.59 80.02 123.38 99.24 90.3 87.61 80.46

80.89 73.14 74.78 107.33 91.04 117.90 84.26 126.82 101.96 93.49 90.27 82.84

61 56 62 24 30 29 37 55 36 61 54 49

63 59 59 27 34 26 41 59 40 63 57 53

4.1. Verification studies The linear TWB model has already been verified in a previous study of the authors by comparing the natural frequencies of the stationary CAS configuration TWB with those determined by the finite element solution performed by MSC Nastran (Farsadi et al., 2017). Unsteady aerodynamics model based on the representation of the lift and the moment by the compressible indicial functions is verified by the experimental study of Barmby et al. (1951). The experimental flutter speed and frequency results obtained for the aluminum alloy wing in the subsonic compressible flow regime are compared with the predictions of the present study. In the model used in the verification study, structural model is based on the geometrically linear equations and exactly the same Ritz based methodology introduced for the nonlinear aeroelastic analysis of the composite wing is applied. In this respect, seven trial functions (N = 7) are used in the series given in Eq. (53) for the bending (v0 ) and torsional (φ ) deformations and first seven (m = 7) eigenfunctions are used for the modal solution of the nonlinear system. Table 1 gives the geometric and material properties of the wings used for the validation of the aeroelastic system based on geometrically linear structural model. In order to adapt the present TWB structural model to the wing models of Barmby et al. degrees of freedom are limited to the flapwise bending deflection (v0 ) and the elastic twist (φ ) and the assignments given by Eq. (89) are made for the stiffness and the mass properties. It should be noted that the a parameter in Table 1 stands for the non-dimensional offset between the shear center and the mid-chord, as shown in Fig. 1. A negative value of the a parameter indicates that the shear center is forward of the mid chord. a66 = GJ , a55 = EI , b4 + b5 = Iα , b1 = ml

(89)

Table 2 compares the wind tunnel results of flutter Mach numbers and flutter frequencies with the results obtained by the proposed compressible indicial aerodynamics formulation. Results given in Table 2 show that present linear aeroelastic analysis results are in good agreement with the experimental results of Barmby et al. (1951) for all variants of the TWB models including the effect of the wing sweep angle Λ. For further confidence on the accuracy of the presented model, another verification is also presented. To this aim, the flutter speed of the Goland cantilevered wing with bending displacement and torsional rotational degrees of freedom in compressible flow reported by Lin and Iliff (2000) is compared with the flutter speed obtained in the present study. Lin and Iliff (2000) have presented a closed-form formula for the lift and moment coefficients of a lifting surface in twodimensional, unsteady, compressible, subsonic flow utilizing an explicit analytical solution of the Possio equation. The

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

17

Table 3 Geometric and material properties of the cantilevered Goland wing used in the verification study (Lin and Iliff, 2000). Parameters

Values

ml Iα EI GJ a b L

0.746 slug/ft 1.943 slug/ft2 23.6 × 106 lb ft 2.39 × 106 lb ft −0.33 3 ft 20 ft

Table 4 Comparison of the flutter speed and Mach number to those of Lin and Ilif (Lin and Iliff, 2000) at the sea level. Lin and Iliff (Lin and Iliff, 2000) Present study

Flutter speed (m/s)

Flutter Mach number

Reduced frequency

136 139

0.40 0.409

0.46 0.44

Table 5 Geometric and material properties of composite CAS configuration wings used in aeroelastic numerical studies. Material properties E1 (GPa) E2 (GPa) E3 (GPa) G12 (GPa) G13 (GPa) G23 (GPa)

ν21 = ν32 = ν31

Geometric properties 206.8 5.17 5.17 3.10 3.10 2.55 0.25

L (m) l (m) d (m) h (m) θ (deg) ρs (kg/m3 ) a Λ (deg) b (m)

10 0.757 0.12 0.03 −45 to −90 1528 −0.3 −40 to +60 0.8

geometric configuration for the test case is reported in Table 3 and comparison of the flutter speeds at the sea level are given in Table 4 which shows that flutter speeds agree considerably well. 4.2. Linear aeroelastic analysis of composite wings Prior to the nonlinear aeroelastic simulations of composite wings, linear aeroelastic analyses are implemented to estimate the stability margins for a better evaluation of the nonlinear aeroelastic simulation results. The Ritz-based methodology is followed in the linear aeroelastic analysis of CAS configuration TWB. Nine trial functions (N = 9) are used in the expansions given in Eq. (53) for all six degrees of freedom. For the modal reduction, the number of mode shapes to be included is found to be 7(m = 7) for an acceptable convergence. Setting the nonlinear vector in Eq. (87) equal to zero, linear aeroelastic equations of motion is obtained. Considering ˜ }T , one obtains the standard q = q0 eλt where λ is the complex valued eigenvalue of the aeroelastic system and q = {ϑ˜ B form of the generalized eigenvalue problem. Geometric and material properties of composite wings structurally modeled as TWBs, which are studied in this section, are given in Table 5. 4.2.1. Fiber angle effect To see the effect of compressibility on the aeroelastic instability margins, linear aeroelastic analysis of the composite wings has been performed by incompressible indicial unsteady aerodynamics utilizing approximation of the Wagner’s function by two aerodynamic lag terms and also by compressible indicial unsteady aerodynamics utilizing approximation of the indicial functions by a total of twelve aerodynamic lag terms. Table 6 presents the critical Mach numbers and the corresponding critical frequencies for unswept wings (Λ = 0) for four different fiber angles of the CAS configuration TWB. From the results given in Table 6, it is evident that compressibility accounts for reduction in aeroelastic instability speeds. Actually, instability speeds obtained by the incompressible unsteady aerodynamics are not conservative. From Table 6, it is also seen that fiber angle has a significant effect on the aeroelastic instability speed. For fiber angles other than −90◦ , the aeroelastic instability type is flutter and as the absolute value of the fiber angle increases, flutter Mach number decreases and for the −90◦ fiber angle case, critical frequency becomes zero indicating divergence instability. For the wing with no sweep, fiber angle of −45◦ results in the highest critical Mach number. It should be noted that for the incompressible case and the fiber angle of −45◦ , the flutter Mach number is high enough for a possible shock formation. However, we still present this result to show how the flutter speed changes with the fiber angle bearing in mind that our model should not be used for speeds above M = 0.8. It should be noted that as Fig. 1 shows, fiber angle is measured from

18

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812 Table 6 Critical Mach number and frequencies of composite CAS configuration wings obtained by compressible and incompressible indicial unsteady aerodynamics. Fiber angle (θ )

−45 −60 −75 −90

Critical frequency (rad/s)

Critical Mach No. Compressible

Incompressible

Difference (%)

Compressible

Incompressible

Difference (%)

0.73 0.67 0.67 0.57

0.87 0.80 0.75 0.65

16.27 16.54 10.24 11.76

34.56 41.82 66.31 0.00

38.5 46.5 69.3 0.00

10.23 10.06 4.31 0.00

the chord-wise axis and when the fiber angle is ±90◦ , fibers are aligned along the wing axis. For negative fiber angles, fibers are aligned towards the leading edge, hence for negative fiber angles other than −90◦ , flutter instability occurs at higher airspeeds compared to the TWBs with positive fiber angles. As a matter of fact, for positive fiber angles, critical Mach numbers are less than 0.3, consequently incompressible unsteady aerodynamics model is sufficient to study the aeroelastic instability characteristics. 4.2.2. Sweep angle effect Fig. 5 shows the effect of the sweep angle on the critical Mach number and critical aeroelastic instability frequency for four different fiber angles of the CAS configuration TWB and for the range of the sweep angle −45◦ ≤ Λ ≤ 45◦ . In Fig. 5, positive sweep angle implies backward-swept wing and negative sweep angle implies forward-swept wing. Critical Mach number results are presented in the range 0.3 ≤ M ≤ 0.8 since the unsteady compressible aerodynamic model based on indicial functions is valid in this range. For the fiber angle of −90◦ , bending-twisting coupling effect due to the material anisotropy is non-existent. For this configuration, Fig. 5 shows that for forward-swept wings, divergence instability occurs and as the sweep angle is increased such that it becomes positive and the wing becomes backwardswept, critical Mach number increases and aeroelastic instability switches to flutter type. This is the typical effect of the wing sweep on the aeroelastic stability characteristics of wings which have no material coupling. It is seen that for fiber angles other than −90◦ , extreme forward-swept wings also have low critical Mach numbers with divergence instability. However, for moderate forward sweep angles, critical Mach number can be increased significantly with the use of fiber angles other than −90◦ with the −45◦ fiber angle being the most effective in increasing the critical Mach number. For the wing with zero sweep angle, the highest critical Mach number occurs for the fiber angle of −45◦ . For moderate backwardswept wings, off-axis fiber angles cause reduction in the critical Mach number. It should be noted that modeling of the compressible unsteady aerodynamics via indicial functions allows frequency domain solution of the flutter problem for the geometrically linear case. For the flutter analysis of composite wings in the compressible flow regime utilizing the geometrically linear TWB model, the present frequency domain approach is deemed to be very advantageous compared to other approaches which generally require time domain solutions, such as coupled CFD and structural finite element methods. Especially, in the initial phase of the wing design, the effect of frequent design changes on the flutter speed and frequency of composite swept wings can be assessed very effectively in terms of model preparation and solution time by the frequency domain solution utilizing the compressible unsteady aerodynamics via indicial functions. 4.3. Nonlinear aeroelastic analysis of composite wings As previously mentioned, nonlinear aeroelasticity of composite wings in the present study is intended to be conducted in time domain by means of the fourth order Runge–Kutta direct integration technique. In order to handle this integration, state space representation given in Eq. (87) is employed as the input for the Runge–Kutta method. It is also worth mentioning that the wing is initially at rest but in order to provide a small disturbance for the transient solution, an initial value of 1.0e−06 is assigned to all modal coordinates. Due to the intrinsic character of aeroelastic systems which fall in the category of self-excited vibrations, imposing a small value of initial condition may ultimately lead to attenuation or magnification of the response. In cases when the steady state aeroelastic response damps out, the system is stable while on the contrary, for unstable conditions, limit cycle oscillations appear. In what follows, time history plots and the associated phase plots and Power Spectrum Density (PSD) diagrams are given for the response of the wing tip (z = L). Initially, the nonlinear aeroelastic system of equations is confirmed by comparing the time domain solution calculated with the linear aeroelastic solution for a rectangular wing with the fiber angle of −75◦ in the CAS configuration. Referring to Table 6, it is seen that the predicted flutter Mach number is 0.67 according to the linear structural analysis. Fig. 6 gives the time history plots of the flapwise wing deflection at the subcritical, critical and supercritical Mach numbers. It should be noted that the low amplitude of oscillations is mainly due to the very stiff TWB structure that is used in the present study. As it is shown in Fig. 6a, in nonlinear analysis for Mach values less than 0.67, disturbance generated by the initial conditions imposed, attenuate due to the aerodynamic damping. With the increase in the Mach number, limit cycle oscillation starts at the bifurcation speed which corresponds to the flutter Mach number obtained by the linear aeroelastic analysis. As shown in Fig. 6b, beyond the critical Mach number, in the supercritical region, amplitude of the limit cycle oscillation increases substantially and as seen in Fig. 6c, the magnitude of oscillations of the flapwise deflection

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

19

Fig. 5. Effect of sweep and lamination angles on the stability margins of composite wings using the compressible indicial aerodynamic model (a) Critical Mach No. vs. the sweep angle (b) Critical frequency vs. sweep angle.

is much higher than the magnitude of the oscillations at the critical condition. For small initial conditions, time domain flutter solution of the geometrically nonlinear problem gives the same flutter speed as the geometrically linear problem. Closeness of the flutter speeds determined by the linear and nonlinear aeroelastic analysis has also been reported by Tang and Dowell (2001). In the present study, the intrinsic character of the nonlinear aeroelastic system has been investigated by exciting the wing by a very small initial disturbance and utilizing linear unsteady aerodynamics. It is to be noted that the due to the very stiff TWB structure that is used in the present study, LCO oscillation amplitudes are very low. For comparison, the study of Shams et al. (2008) is taken as reference. 4.3.1. Compressibility effect For a wing with a fiber angle of θ = −75◦ , Fig. 7 shows the bifurcation diagram of the torsional (φ ) and the flapwise (v0 ) degrees of freedom at the wing tip calculated both with compressible and incompressible indicial aerodynamic theories. The incompressible model is the one proposed by Farsadi et al. (2018) while the compressible model is the indicial formulation presented in the current study. A delay of the onset of flutter (Hopf bifurcation) is found in the incompressible model. Accordingly, flutter speed determined by the incompressible model is overestimated, which can be dangerous if the flight vehicle is designed to operate near its flutter boundaries. This shows the significance of compressible aerodynamic modeling for a more precise aeroelastic design. Moreover, Fig. 7 also shows that the slopes of the bifurcation curves in the compressible model are more significant compared to incompressible curves. This implies that the sensitivity of the LCO amplitude to the increase of the Mach number is more dramatic in the compressible aerodynamic model in the post critical zone. 4.3.2. Fiber angle effect Bifurcation diagrams of composite wings with the CAS configuration TWB structural model are presented in Fig. 8 for fiber angles of −75◦ and −45◦ . For the CAS configuration TWB with fiber angles of −75◦ and −45◦ and zero sweep angle (Λ = 0◦ ) , bifurcation points are at Mach numbers 0.67 and 0.736, respectively. Bifurcation points identified by the nonlinear aeroelastic analysis coincide with the flutter Mach numbers determined by the linear aeroelastic analysis given in Table 6. Fig. 8a shows that for the −45◦ fiber angle case, not only the instability occurs at a higher Mach number, but post-flutter limit cycle oscillation amplitudes of the flapwise and torsional deformation of the wing tip are also lower than the −75◦ fiber angle case. This observation clearly shows that post-flutter response of the composite wing with fiber angle of −45◦ is also more well behaved compared to the −75◦ fiber angle case. It should be noted at Mach numbers sufficiently higher than the critical Mach number, nonlinear aeroelastic response becomes either quasi-periodic or chaotic, and for these responses, maximum amplitude of the flapwise wing tip deflection and the maximum wing twist angle are plotted in the bifurcation diagram shown in Fig. 8a. In the presented results, limit cycle oscillations are initiated through a Hopf bifurcation and this continues to the point where quasi-periodic oscillations emerge. The initiation of quasi-periodic oscillations corresponds to the period-doubling (PD) bifurcation. The PD bifurcation then continues and provides a route to chaos (C). These two types of bifurcations are now marked on the bifurcation plots by using PD and C symbols. As the sweep angle increases from zero to a larger value, such as Λ = 45◦ , significant qualitative and quantitative differences appear in the nonlinear response. For the 45◦ sweep angle, results are shown in Fig. 8b. In this case, for the fiber angles

20

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

Fig. 6. Time histories of the flapwise wing tip deflection, Fiber angle (θ ) = −75◦ (a) Subcritical (M = 0.99Mcritical = 0.6633) (b) Subcritical (M = Mcritical = 0.67, ωf = 66.31 rad/s) (c) Subcritical (M = 1.01Mcritical = 0.6767).

Fig. 7. Bifurcation diagram for flapwise and torsion degrees of freedom using compressible and incompressible aerodynamic models θ = −75◦ .

of −75◦ and −45◦ , bifurcation occurs at Mach numbers 0.779 and 0.529, respectively. Comparing bifurcation diagrams of the unswept and swept wing configurations in Fig. 8, one can see that for the fiber angle of −75◦ , the critical Mach number increases when the backward sweep angle increases and the reverse behavior is observed for the −45◦ fiber angle case. It should be noted that the critical speeds determined by the nonlinear aeroelastic analysis through the bifurcation diagrams are in perfect agreement with the flutter speeds, which are obtained by the linear aeroelastic analysis, given in

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

21

Fig. 8. The effect of lamination angle on the bifurcation diagram of the composite swept and unswept wings.

Fig. 5a. Moreover, although the overall trend of the response for a specific fiber angle remains more or less tact, the LCO amplitude decreases at nonzero sweep angles compared to the unswept wing configuration. For the unswept wing with the CAS configuration fiber angle of −75◦ , Figs. 9–11 give the time histories of the flapwise wing tip displacement and wing tip twist and their respective power spectral density functions (PSDs) at postflutter Mach numbers, M /Mcritical = 1.01, 1.05, 1.09. From Fig. 9, it is seen that for the low post-critical Mach number (M /Mcritical = 1.01), aeroelastic response is purely periodic. At the post-flutter Mach ratio of 1.05, Fig. 10 shows that oscillation amplitudes increase and the response is quasi-periodic and responses are composed of more than one dominant frequency. Finally, as seen in Fig. 11, at the post-flutter Mach ratio of 1.09, nonlinear aeroelastic response becomes chaotic and PSD plots show that broadband range of dominant frequencies exist in the responses. 4.3.3. Sweep angle effect For the composite wing with the CAS configuration fiber angle of −75◦ , Fig. 12 gives the bifurcation diagram which shows the effect of the sweep angle on the nonlinear aeroelastic response of the wing. Fig. 12 compares the oscillation amplitudes of the flapwise tip deflection and tip twist of the unswept (Λ = 0◦ ) and the backward-swept wing (Λ = 30◦ , 45◦ ). For the backward-swept wings, critical Mach numbers determined by the nonlinear aeroelastic analysis at Λ = 30◦ , 45◦ are 0.692 and 0.779, respectively. These results agree with the linear aeroelastic analysis as depicted in Fig. 5. As seen in Fig. 12, for the −75◦ fiber angle configuration, with the introduction of the wing sweep, critical Mach number increases and amplitudes of oscillations reduce compared to the unswept wing at the same post-critical Mach number.

22

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

Fig. 9. Response qualification of the composite wing with θ = −75◦ and zero sweep angle at M /Mcritical = 1.01 (a) Flapwise deformation (b) Torsional deformation.

Fig. 10. Response qualification of the composite wing with θ = −75◦ and zero sweep angle at M /Mcritical = 1.05, (a) Flapwise deformation (b) Torsional deformation (c) PSD of flapwise deformation (d) PSD of torsional deformation.

To aid the interpretation of the nonlinear aeroelastic response of the composite wing with the −75◦ fiber angle configuration, phase plots are prepared for the unswept wing and backward-swept wings with sweep angles Λ = 15◦ , 30◦ , 45◦ . Figs. 13–15 present the phase portraits for the flapwise bending motion of the wing tip, for the Mach ratios of 1.01, 1.05 and 1.09, respectively. As Fig. 13 shows, at the low post-critical Mach ratio (M /Mcritical = 1.01), nonlinear aeroelastic response of the wing is purely periodic with narrow closed orbits in the phase plot for all sweep angles. For the post-critical Mach ratio of 1.05, Fig. 14 shows that nonlinear aeroelastic responses of wings with low sweep angles show quasi-periodic behavior which appears as closed thicker circuits in the phase plane. For the unswept wing, in the phase plane, there are two main thicker loops and as the sweep angle is increased, response of the wing becomes periodic

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

23

Fig. 11. Response qualification of the composite wing with θ = −75◦ and zero sweep angle at M /Mcritical = 1.09, (a) Flapwise deformation (b) Torsional deformation (c) PSD of flapwise deformation (d) PSD of torsional deformation.

Fig. 12. The effect of sweep angle on the bifurcation diagram of the composite wing with θ = −75◦ .

again. Finally, for the higher post-critical Mach ratio of 1.09, phase plots given in Fig. 15 show that nonlinear aeroelastic responses of the unswept wing and the wing with the sweep angle of 15◦ are chaotic which is evident from the numerous intermittent circuits. However, as the sweep angle is increased, nonlinear aeroelastic response of the wing tip becomes quasi-periodic again. Phase plots given in Figs. 13–15 all show that for the −75◦ fiber angle configuration, by increasing the backward sweep angle, nonlinear response of the wing becomes more well-behaved and the amplitude of oscillations of the flapwise bending motion of the wing tip reduce. Due to the combined effect of the lowering of the effective angle

24

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

Fig. 13. Phase plane plots of the flapwise bending motion for different sweep angles at M /Mcritical = 1.01 and θ = −75◦ , (a) Λ = 0◦ (b) Λ = 15◦ (c) Λ = 30◦ (d) Λ = 45◦ .

of attack and the dynamic pressure, as the backward sweep angle is increased, nonlinear response of the wing becomes more well-behaved. Figs. 16–18 show the Poincare plots which are constructed for the flapwise bending degree of freedom of the wing tip (v0 ) at the section φ = 0 by plotting crossings in both directions. At the low post-critical Mach ratio (M /Mcritical = 1.01), Fig. 16 shows that for all sweep angles, there are single dots which are indicative of the periodic response. For the post-critical Mach ratio of 1.05, Fig. 17 shows that for the unswept wing, there are finite number of dots indicating quasi-periodic response. As the sweep angle is increased, number of dots in the Poincare maps decreases, and for the sweep angle of 45◦ , there is again single dot which indicates purely periodic response which agrees with the last phase plot given in Fig. 13. For the post-critical Mach ratio of 1.09, Poincare maps, given by the first two plots in Fig. 18, for the unswept wing and the wing with the sweep angle of 15◦ have numerous dots indicating a chaotic response. It is seen that for the thicker orbits in the phase plane, the number of dots emerging in the Poincare plots is also high. So in the post-critical region, there are routes to chaos and the aeroelastic behavior may be chaotic depending on the configuration. 5. Conclusion In the present study, geometrically nonlinear aeroelastic analysis of composite wing structures is performed using a thin walled beam theory and taking the non-classical effects into account. The structural equations of motion are obtained for the CAS configuration TWB by including up to third order nonlinear terms. The unsteady compressible aerodynamics model is constructed in the time domain by using a rather involved interpolation method which allows the calculation of indicial functions at any Mach number up to 0.8, by incorporating twelve aerodynamic lag states. Calculation of the indicial functions at any Mach number allows frequency domain solution of the linear aeroelastic problem for the very fast determination of the flutter speed for the geometrically linear structural model. To construct the nonlinear coupled-field system of equations, twelve auxiliary equations governing the aerodynamic lag states are included in the formulation. Nonlinear aeroelastic system of equations for the composite wing is solved using the linear base functions to approximate the spatial variation of the degrees of freedom of the TWB. Time response of the nonlinear aeroelastic system is obtained via the Runge–Kutta direct integration algorithm. Effects of fiber angle and the sweep angle of the CAS configuration TWB on the post-critical response of the composite wing are studied by providing bifurcation diagrams, phase portraits

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

25

Fig. 14. Phase plane plots of the flapwise bending motion for different sweep angles at M /Mcritical = 1.05 and θ = −75◦ , (a) Λ = 0◦ (b) Λ = 15◦ (c) Λ = 30◦ (d) Λ = 45◦ .

and one side Power Spectral Density (PSD) plots. It is concluded that for the flutter analysis of composite wings in the compressible flow regime utilizing the geometrically linear TWB model, the present frequency domain approach which is based on modeling of the compressible unsteady aerodynamics through the compressible indicial functions, is deemed to be very advantageous. Frequency domain solution allows very fast calculation of the flutter speed and the associated flutter frequency in the compressible flow regime by performing a Mach sweep. Results show that critical speeds determined by the linear aeroelastic analysis of the composite wings via the frequency domain solution agree very well with the bifurcation speeds obtained by the time domain solutions of the nonlinear aeroelastic system of equations of the composite wing structurally modeled as TWB. Moreover, it is shown that compressibility accounts for reduction in aeroelastic instability speeds and instability speeds obtained by the incompressible unsteady aerodynamics are not conservative. For the composite wing with the CAS configuration TWB structural model, fiber angle is seen to be a very influential parameter on the instability speed of the composite wing. For negative fiber angles, the form of the instability is of flutter type because fibers are oriented towards the leading edge and divergence instability is deferred and flutter instability occurs at higher airspeeds and the compressibility effect is not negligible compared to the TWBs with positive fiber angles. Furthermore, it is shown that combination of the fiber angle of the CAS configuration and the sweep angle of the wing has significant effect on the critical speeds as well as on the post-flutter nonlinear aeroelastic response of the composite wing. For the fiber angle of −75◦ , it is shown that as the backward sweep angle of the wing is increased, post-flutter aeroelastic response of the wing becomes more well behaved compared to wings with low backward sweep angle with lower LCO amplitude. This study presents nonlinear aeroelasticity of wings based on geometrically nonlinear thin walled beam model by retaining nonlinear terms of up to third order and compressible unsteady aerodynamics based on indicial functions for the investigation of the post-critical aeroelastic response of composite wings with sweep. The novel aspects of the study is the explicit derivation that ends up with final form of nonlinear aeroelastic system of equations given by Eq. (87) and the interpolation method employed for the evaluation of the indicial functions at any Mach number in the subsonic compressible range. It is believed that present study provides results which can be used for comparison purposes in similar studies which utilize different geometrically nonlinear structural wing and compressible unsteady aerodynamic models in an aeroelastic framework.

26

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

Fig. 15. Phase plane plots of the flapwise bending motion for different sweep angles at M /Mcritical = 1.09 and θ = −75◦ , (a) Λ = 0◦ (b) Λ = 15◦ (c) Λ = 30◦ (d) Λ = 45◦ .

Nomenclature The following is a list of symbols used in the article along with their respective definitions. AR : Aspect ratio CAS : Circumferentially Asymmetric Stiffness lay-up configuration TWB : Thin Walled Beam u(x, y, z , t); v (x, y, z , t); w (x, y, z , t) : Displacement field of the thin walled beam u0 , v0 , w0 : Translations of the shear center along x, y, z axes, respectively θx , θy , φ : Rotations of the cross section about x, y, z axes, respectively b, l, d, h, L : Half chord, width, height, thickness and length of the thin walled beam, respectively x, y, z : Cartesian fixed coordinate system located at the root of the thin walled beam s, z , n : Local coordinate system located at the mid plane of the cross section of the thin walled beam rn : Geometric quantity Fw (s) : Primary warping function na(s) : Secondary warping function a : Dimensionless coefficient c : Speed of sound Ψ : Torsional function T , V : Total kinetic and strain potential energies, respectively Wn.c : Work done by external non-conservative loads Nl : Number of composite layers k : Index of the layer number εij , σij : Strain and stress components, respectively γzy , γzx : Transverse shear measures of the cross section Tz : Generalized beam axial force

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

27

Fig. 16. Poincare plots for a variety of sweep angles at the plane of φ = 0, M /Mcritical = 1.01 and θ = −75◦ , (a) Λ = 0◦ (b) Λ = 15◦ (c) Λ = 30◦ (d) Λ = 45◦ .

Qx , Qy : Chord-wise and flap-wise shear forces, respectively Mx , My , Mz : Flap-wise bending moment, chord-wise bending moment, Saint-Venant twist moment, respectively Bw : bi-moment or warping torque Λz : Torque due to the higher order stress couple aij : Stiffness matrix coefficients ∂X ∂2X , ∂t ∂t2 ′ ′′ ∂ X ∂ 2 X X , X : ∂ z , ∂ z2

X˙ , X¨ :

R : Position vector of an arbitrary point ρs : Mass density of the thin walled beam ρs(k) : Mass density of the thin walled beam at kth layer ml : Mass per unit length of the thin walled beam Iα : Mass moment of inertia Lae , Tae : Unsteady aerodynamic lift and moment, respectively ϑ (t) : Generalized modal coordinates ψ ∆ : Vector of trial functions ∆ : Any of the six degrees of freedom u0 , v0 , w0 , θx , θy , φ EI , GJ : Bending and torsion rigidities, respectively λ : Complex values eigenvalue of the aeroelastic system ωf : Frequency corresponding to critical flutter condition w ˆ a (z , t), φˆ a (z , t) : Plunging and pitching velocities of the airfoil, respectively wa (x, z , t) : Downwash velocity corresponding to plunging and pitching motions Λ : Sweep angle U∞ , U : Free stream fluid velocity ρ∞ , ρ : Mass density of the free stream CLφ : Local lift curve slope

28

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

Fig. 17. Poincare plots for a variety of sweep angles at the plane of φ = 0, M /Mcritical = 1.05 and θ = −75◦ , (a) Λ = 0◦ (b) Λ = 15◦ (c) Λ = 30◦ (d) Λ = 45◦ .

t , t0 , τ : Dimensional; non-dimensional time variables M : Mach number φc (t) : Indicial lift function due to unit step change of the plunging motion φcM (t) : Indicial moment function about the leading edge due to unit step change of the plunging motion L′T (z , t), MT′ (z , t) : Lift and moment due to plunging motion evaluated at the reference coordinate at the leading edge φcq (t) : Indicial lift function due to unit step change of the pitching rate φcMq (t) : Indicial moment function about the leading edge due to unit step change of the pitching rate L′q (z , t), Mq′ (z , t) : Lift and moment due to pitching motion evaluated at the reference coordinate at the leading edge L′ (z , t), M ′ (z , t) : Total lift (positive upwards) and moment (about the leading edge and positive nose down) LT (z , t), MT (z , t) : Lift and moment due to plunging motion evaluated at the Theodorsen’s coordinate (behind the leading edge) Lq (z , t), Mq (z , t) : Lift and moment due to pitching motion evaluated at the Theodorsen’s coordinate (behind the leading edge) φ¯ c (t), φ¯ cM (t), φ¯ cq (t), φ¯ cMq (t) : Indicial functions for compressible flow in Theodorsen’s coordinate αic , αicq , αicM , αicMq : Indicial Mach dependent coefficients evaluated at the leading edge φ˜ c (t) : Approximated Indicial function for plunging motion α¯ ic , α¯ icq , α¯ icM , α¯ icMq : Indicial Mach dependent coefficients evaluated at the Theodorsen’s coordinate (behind the leading edge) Bic (z , t) : Aerodynamic lag terms θ : Fiber angle ψ u , ψ v , ψ w , ψ x , ψ y , ψ z : Trial functions ηu , ηv , ηw , ηx , ηy , ηz : Time dependent variables G : Vector of state variables Ms , Cs , Ks : Structural mass, damping and stiffness matrices Mae , Cae , Kae : Aerodynamic mass, damping and stiffness matrices

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

29

Fig. 18. Poincare plots for a variety of sweep angles at the plane of φ = 0, M /Mcritical = 1.09 and θ = −75◦ , (a) Λ = 0◦ (b) Λ = 15◦ (c) Λ = 30◦ (d) Λ = 45◦ .

Mt , Ct , Kt : Reduced order mass, damping and stiffness matrices L, R : Left and right eigenvectors of the linear equations, respectively Z : Vector of lag states ˆ ic (t), Bˆ icq (t), Bˆ i (t), Bˆ i (t) : Transformed aerodynamic lag state vectors B cM cMq ¯ NL : Nonlinear vector and reduced order nonlinear vector, respectively HNL , H ν12 , ν13 , ν23 : Poisson’s ratio Un : Normal component of the free stream velocity E1 , E2 , E3 : Young’s modulus in principle directions G12 , G13 , G23 : Shear modulus in principle directions N : Number of trial functions in the solution algorithm m : Number of reduced order modes taken into account Nu2,v,w,x,y,φ ,Nu3,v,w,x,y,φ : Second and third order nonlinear terms βi : Mach-independent power coefficients Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Appendix A. Compressible indicial functions The indicial lift function, (φc (t)), and the indicial moment function about the reference axis located at the leading edge (φcM (t)) are defined as the responses to the unit step change of the plunging motion w ˆ a at the leading edge. Consequently, in analogy with the incompressible flow, the lift and moment due to the plunging motion w ˆ a (z , t) are defined as,

[

L′ T (z , t) = −CLφ ρ U cos Λb w ˆ a (z , 0)φc (t) +

t

∫ 0

dw ˆ a (z , τ ) dτ

φc (t − τ ) dτ

] (90)

30

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

[

t

∫

M ′ T (z , t) = −2CLφ ρ U cos Λb2 w ˆ a (z , 0)φcM (t) +

dw ˆ a (z , τ ) dτ

0

φcM (t − τ )dτ

] (91)

In a similar fashion, φcq (t) and φcMq (t) are introduced as the indicial lift and the moment functions about the leading edge due to the unit step change of the pitching rate (φˆ a ) at the leading edge. As a result, the lift and moment due to the pitching motion φˆ a (z , t) are defined by Eqs. (92) and (93).

[ L q (z , t) = 2CLφ ρ U cos Λb

2

M q (z , t) = 4CLφ ρ U cos Λb

3

′

φˆ a (z , 0)φcq (t) +

dτ

0

[ ′

dφˆ a (z , τ )

t

∫

φˆ a (z , 0)φcMq (t) +

φcq (t − τ ) dτ

dφˆ a (z , τ )

t

∫

]

dτ

0

(92)

] φcMq (t − τ ) dτ

(93)

Total aerodynamic lift (positive upwards) and moment about the leading edge (positive nose-up) are expressed as, L′ (z , t) = L′ T (z , t) + L′ q (z , t)

(94)

M ′ (z , t) = M ′ T (z , t) + M ′ q (z , t)

(95)

In Eqs. (94)–(95), the prime symbol indicates that the aerodynamic loads are evaluated in the reference coordinate located at the leading edge of the wing. Aerodynamic loads in the so called Theodorsen’s coordinate, which is located at a distance b(a + 1) behind the leading edge, are then defined by Eqs. (96)–(99) which are obtained by employing the general law for transferring the axis of a moment (Bisplinghoff et al., 2013).

[

ˆ a (z , 0)φ¯ c (t) + LT (z , t) = L′T (z , t) = −CLφ ρ U cos Λb w

t

∫

dw ˆ a (z , τ ) dτ

0

φ¯ c (t − τ ) dτ

] (96)

MT (z , t) = M ′ T (z , t) + b(a + 1)L′ T (z , t)

(97)

[

ˆ a (z , 0)φ¯ cM (t) + = −2CLφ ρ U cos Λb2 w [

t

∫

dw ˆ a (z , τ ) dτ

0

Lq (z , t) = Lq (z , t) − CLφ ρ U cos Λb (a + 1) φˆ a (z , 0)φc (t) + 2

′

[ = 2CLφ ρ U cos Λb

2

φˆ a (z , 0)φ¯ cq (t) +

t

∫ 0

dφˆ a (z , τ ) dτ

φ¯ cM (t − τ )dτ t

∫

dφˆ a (z , τ )

Mq (z , t) = M q (z , t) − 2CLφ ρ U cos Λb (a + 1) φˆ a (z , 0)φcM (t) +

[ + Lq (z , t)b(a + 1) = 4CLφ ρ U cos Λb

φc (t − τ ) dτ

(98)

φ¯ cq (t − τ ) dτ

[

3

]

]

3

′

dτ

0

]

φˆ a (z , 0)φ¯ cMq (t) +

t

∫

dτ

0 t

∫ 0

dφˆ a (z , 0)

] φcM (t − τ )dτ

dφˆ a (z , τ ) dτ

(99)

] φ¯ cMq (t − τ ) dτ

In Eqs. (96)–(99) a new set of aerodynamic indicial functions for compressible flow are defined in the Theodorsen’s coordinate as

φ¯ c (t) = φc (t)

(100) a

1

φ¯ cM (t) = φcM (t) + ( + )φc (t) 2 a

2 1 φ¯ cq (t) = φcq (t) − ( + )φc (t) 2 2 a 1 a 1 φ¯ cMq (t) = φcMq (t) + ( + )[φcq (t) − φcM (t)] − ( + )2 φc (t) 2 2 2 2

(101) (102) (103)

It should be noted that the indicial functions corresponding to compressible aerodynamics are actually Mach number dependent, but this dependency is not explicitly shown in Eqs. (100)–(103) to simplify the writing.

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

31

Appendix B. Linear/nonlinear parameter definitions in aeroelastic equations of motion f1 (∆) = b1 u¨ 0 f2 (∆) = a22 (θy + u′ 0 ) + a12 w ′ 0 + Nu2 + Nu3 f3 (∆) = b1 v¨ 0 − Lae f4 (∆) = a33 (θx + v ′ 0 ) − a37 φ ′′ + Nv2 + Nv3 f 5 (∆ ) = b1 w ¨0 f6 (∆) = a11 w ′ 0 + a12 θy + Nw2 + Nw3 1

f7 (∆) = a33 (v ′ 0 + θx ) − a37 φ ′′ + (b4 + b12 )θ¨x + Nx2 + Nx3 f8 (∆) = a55 θ

′

f9 (∆) = a12 w

x

′

+ a56 φ + ′

0

2 Nx2

1

(104)

2 Nx3

+

1

1

+ a22 (u 0 + θy ) + (b5 + b11 )θ¨y + Ny2 + Ny3 ′

2

f10 (∆) = a44 θ ′ y + Ny2 + Ny3

2 1

1

f11 (∆) = (b4 + b5 )φ¨ + Nφ2 + Nφ3 − Mae 2

2

f12 (∆) = a66 φ ′ + a56 θ ′ x + b10 φ¨ ′ + Nφ2 + Nφ3 3

f13 (∆) = a37 (v ′ 0 + θx ) − a77 φ ′′ + Nφ2 + Nφ3

3

where the quadratic N 2 and cubic N 3 nonlinear terms are defined as, Nu2 = a22 v ′ 0 φ − a33 (v ′ 0 φ + θx φ ) − a55 θ ′ x φ ′ + ( 21 a28 − a56 )φ ′

2

+a12 ( 32 u′ 0 2 + 12 v ′ 0 2 + θy u′ 0 ) + a11 u′ 0 w′ 0 + a37 φ ′′ φ Nv2 = a12 (u′ 0 v ′ 0 + v ′ 0 θy + w ′ 0 φ ) + a22 (u′ 0 φ + θy φ ) − a33 u′ 0 φ + a44 φ ′ θ ′ y + a11 v ′ 0 w ′ 0 Nw2 =

a11 2 (u′ 0 2

a18 2

+ v ′ 0 2 ) + a12 v ′ 0 φ +

φ′2

1

Nx2 = −a33 u′ 0 φ 2

Nx2 = −a55 u′ 0 φ ′ 1

Ny2 =

a12 2 (u′ 0 2

+ v ′ 0 2 ) + a22 v ′ 0 φ +

a28 2

φ′2

2

Ny2 = a44 v ′ 0 φ ′ 1

Nφ2 = a22 (u′ 0 v ′ 0 + v ′ 0 θy ) − a33 (u′ 0 v ′ 0 + u′ 0 θx ) + a12 v ′ 0 w ′ 0 + a37 u′ 0 φ ′′ 2

Nφ2 = a44 v ′ 0 θ ′ y + a28 (u′ 0 φ ′ + θy φ ′ ) − 2a56 u′ 0 φ ′ − a55 u′ 0 θ ′ x + a18 w ′ 0 φ ′ 3

Nφ2 = −a37 u′ 0 φ Nu3 = −a22 (u′ 0 φ 2 + 21 θy φ 2 ) − a12 ( 21 w ′ 0 φ 2 + u′ 0 v ′ 0 φ ) +

a18 ′ 2 u0 ′ 2 a12 2 2 2 ) ( u′ 0 3 ′0 ) 2 2 a44 ′ 0 ′ a55 ′ x ′ a56

+a55 u′ 0 φ ′ 2 + a33 u′ 0 φ 2 − a44 θ ′ y φ ′ φ + Nv3 = a22 v ′ 0 φ 2 − a33 ( 12 θx φ 2 + v ′ 0 φ

v

+v

=−

φ

a11 2 ′ 3 ( ′ 0 u′ 0 0 ) 2 a12 ′ 3 2 Nw u0 2 a33 ′ 31 2 Nx 0 2

+

=−

2 Nx3

31

Ny

32

Ny

+

a18 2

v 0φ + ′2

′

+ u′ 0 v ′ 0 2 )

φ

φ

−

+ φv

v φ −

θ φφ−

φφ ′ 2

v φ

= −a55 v ′ 0 φ ′ φ =

a22 ′ u0 2

φ2

= a44 u′ 0 φ ′ φ

31

Nφ = −a55 v ′ 0 φ ′ φ − a56 v ′ 0 φ ′ − 2

−a44 (u′ 0 φ ′ θ ′ y ) − a22 (u′ 0 θy φ + φ 32

Nφ = a55 (φ

+

a11 3 (u′ 0 2

a18 ( 2

φ

′ ′ 2 u0

′ ′ 2 u0

v

− φv

− v ′ 0 θ ′ x φ ) + a44 (φ v

+φ v

′ ′ 2 0 )

′ ′ 2 0

+

a88 2

φ

′3

+ v ′ 0 3 + u′ 0 w ′ 0 φ )

a12 2 ( ′ 0 u′ 0 2 2 ′ 2 u′ 0 0 )

+ a33 (v ′ 0 θx φ + φv ′ 0 2 − φ u′ 0 2 )

− u′ 0 θ ′ y φ ) − 2a56 (v ′ 0 φ ′ φ )

(105)

32

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

Appendix C. Aerodynamic mass, damping and stiffness matrices

⎡

⎤

⎢0 ⎢ ⎢0 ⎢ ⎢ ⎢ ⎢0 Kae = ⎢ ⎢ ⎢0 ⎢ ⎢ ⎢0 ⎣ 0 Ka22 (i, j) =

0

0

0

0

0 ⎥

Ka22

0

0

0

Ka26 ⎥

0

0

0

0

0

0

0

0

0

0

0

0

Ka62

0

0

0

L

∫

⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ ⎦ (106)

Ka66 vT

CLφ ρ U 2 bα0c cos Λ sin Λψiv ψ ′ j dz

0

L

∫

φ

∫

L

φT

CLφ ρ U 2 b2 α0cq sin 2Λψiv ψ ′ j dz CLφ ρ (U cos Λ) bα0c ψi ψj dz − , j) = − 0 ∫ L0 vT φ CLφ ρ U 2 b2 α0cM sin 2Λψi ψ ′ j dz Ka62 (i, j) = 0 ∫ ∫ L L φT φ φ φT CLφ ρ U 2 b3 α0cMq sin 2Λψi ψ ′ j dz CLφ ρ (U cos Λ)2 b2 α0cM ψi ψj dz − 2 Ka66 (i, j) = −2 0 0 ⎤ ⎡ Ka26 (i

⎢0 ⎢ ⎢0 ⎢ ⎢ ⎢ ⎢0 Cae = ⎢ ⎢ ⎢0 ⎢ ⎢ ⎢0 ⎣ 0 Ca22 (i, j) =

0

0

0

0

0 ⎥

Ca22

0

0

0

Ca26 ⎥

0

0

0

0

0

0

0

0

0

0

0

0

Ca62

0

0

0

L

∫ 0

Ca26 (i, j) = −2

, j) = 2

Ca62 (i

2

⎢0 ⎢ ⎢0 ⎢ ⎢ ⎢ ⎢0 Mae = ⎢ ⎢ ⎢0 ⎢ ⎢ ⎢0 ⎣ 0

L

∫ L

0

φ

CLφ ρ U cos Λb2 α0cq ψi ψjv dz T

φT

CLφ ρ U cos Λb2 α0cM ψiv ψj dz L

∫ 0

φ

φT

CLφ ρ U cos Λb3 α0cMq ψi ψj dz

⎤ 0

0

0

0

0 ⎥

Ma22

0

0

0

Ma26 ⎥

0

0

0

0

0

0

0

0

0

0

0

0

Ma62

0

0

0

, j) = π ρ∞ U cos Λ

∫L

a (i, j) = π ρ∞ U cos Λ M26

∫L

a M22 (i

(107)

Ca66

CLφ ρ U cos Λbα0c ψiv ψj dz

∫

, j) = −4 ⎡

⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ ⎦ vT

0

Ca66 (i

vT

a M62 (i, j) = π ρ∞

∫L

a M66 (i, j) = π ρ∞

∫L

0 0

0 0

⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ ⎦

Ma66 v

(108) 6N ×6N

vT

b ψi ψ j dz 2

φ

ab3 ψi ψ vj T dz φT

ab3 ψiv ψ j dz φ

φT

b4 ( 81 + a2 )ψi ψ j dz

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

33

j

Appendix D. Nonlinear coefficients N¯ i

N¯ u2 = a22

m m ∑ ∑

vT

m m ∑ ∑

(ψ ′ Rv )r (ψ φ Rφ )s ϑr ϑs − a33 T

r =1 s=1 m m

r =1 s=1 m m

− a33

∑∑

(ψ x Rx )r (ψ φ Rφ )s ϑr ϑs − a55 T

T

+ ( a28 − a56 ) 2

+

2

m m ∑ ∑

(ψ ′

φT

Rφ )r (ψ ′

φT

Rφ )s ϑr ϑs +

r =1 s=1

a12

+ a11

∑∑

m m ∑ ∑

vT

vT

(ψ ′ Rv )r (ψ ′ Rv )s ϑr ϑs + a12

(ψ

T ′u

(ψ

T ′u

R u ) r (ψ

T ′w

Ru )r (ψ

T ′v

Rw )s ϑr ϑs + a37

r =1 s=1

N¯ v2 = a12

m m ∑ ∑

∑∑

(ψ ′

wT

Rv )s ϑr ϑs + a12

+ a22 + a44

T

(ψ y Ry )r (ψ φ Rφ )s ϑr ϑs − a33 T

(ψ ′

a11 ∑ ∑ 2

+ a12

T

φT

1

2

vT

∑∑

∑∑ r =1 s=1 m m

1

N¯ y2 =

a12 ∑ ∑ 2

+ a22

r =1 s=1 m m

∑∑

yT

uT

T

2

N¯ y2 = a44

m

∑∑ r =1 s=1 m m

1

N¯ φ2 = a22

∑∑ r =1 s=1 m m

+ a22

∑∑

uT

T

(ψ ′′

φT

Rφ )r (ψ φ Rφ )s ϑr ϑs T

vT

(ψ ′ Rv )r (ψ y Ry )s ϑr ϑs T

uT

T

uT

(ψ ′ Ru )r (ψ φ Rφ )s ϑr ϑs vT

(ψ ′ Rv )r (ψ ′

∑∑

uT

T

uT

φT

a11 ∑ ∑ 2

wT

vT

Rw )s ϑr ϑs

vT

(ψ ′ Rv )r (ψ ′ Rv )s ϑr ϑs

r =1 s=1 m m

a18 ∑ ∑ 2

T

(ψ ′

φT

Rφ )r (ψ ′

φT

(ψ ′ Ru )r (ψ ′

r =1 s=1

uT

vT

uT

(ψ ′ Rv )r (ψ φ Rφ )s ϑr ϑs + T

vT

φT

uT

vT

vT

T

uT

vT

(ψ ′ Rv )r (ψ ′

(112)

Rφ )s ϑr ϑs

(ψ ′ Ru )r (ψ ′ Ru )s ϑr ϑs +

− a33

r =1 s=1

(111)

Rφ )s ϑr ϑs

(113) m m a12 ∑ ∑

2

r =1 s=1 m m

a28 ∑ ∑ 2

vT

vT

(ψ ′ Rv )r (ψ ′ Rv )s ϑr ϑs

(ψ ′

φT

Rφ )r (ψ ′

φT

(114)

Rφ )s ϑr ϑs

r =1 s=1

Rφ )s ϑr ϑs

(115)

(ψ ′ Ru )r (ψ ′ Rv )s ϑr ϑs

(116)

(ψ ′ Rv )r (ψ y Ry )s ϑr ϑs

r =1 s=1 m m ∑ ∑

(110)

(ψ ′ Ru )r (ψ φ Rφ )s ϑr ϑs

∑∑

(ψ ′ Ru )r (ψ φ Rφ )s ϑr ϑs

r =1 s=1 m

uT

(ψ y Ry )r (ψ ′ Ru )s ϑr ϑs

r =1 s=1 m m

(ψ ′ Rv )r (ψ φ Rφ )s ϑr ϑs +

r =1 s=1 m m

N¯ x2 = −a55

m m ∑ ∑

m m ∑ ∑

Rφ )r (ψ ′ Ry )s ϑr ϑs + a11

uT

r =1 s=1 m m

N¯ x2 = −a33

uT

(ψ ′ Ru )r (ψ ′ Ru )s ϑr ϑs

r =1 s=1

m m ∑ ∑

r =1 s=1 m m

(ψ ′ Ru )r (ψ ′ Ru )s ϑr ϑs +

r =1 s=1 m m

∑∑

Rφ )s ϑr ϑs

r =1 s=1

r =1 s=1 m m

N¯ w2 =

m m ∑ ∑

a12

2

m m ∑ ∑

Rw )r (ψ φ Rφ )s ϑr ϑs + a22

r =1 s=1 m m ∑ ∑

3

r =1 s=1 m m

r =1 s=1 m m ∑ ∑

φT

r =1 s=1

r =1 s=1 m m

+ a12

xT

(ψ ′ Rx )r (ψ ′

(109)

r =1 s=1

r =1 s=1 m m ∑ ∑

T

r =1 s=1

r =1 s=1

1

1

vT

(ψ ′ Rv )r (ψ φ Rφ )s ϑr ϑs

(ψ ′ Ru )r (ψ ′ Rv )s ϑr ϑs − a33

m m ∑ ∑ r =1 s=1

uT

T

(ψ ′ Ru )r (ψ x Rx )s ϑr ϑs

34

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

+ a12

m m ∑ ∑

vT

(ψ ′ R v ) r (ψ ′

2

vT

yT

∑∑

uT

φT

(ψ ′ R u ) r (ψ ′

r =1 s=1 m m

+ a28

T

∑∑

(ψ y Ry )r (ψ ′

Rφ )s ϑr ϑs

φT

Rφ )s ϑr ϑs

Rφ )s ϑr ϑs − 2a56

xT

uT

(ψ ′ Ru )r (ψ ′ Rx )s ϑr ϑs + a18

m m ∑ ∑

uT

φT

Rφ )s ϑr ϑs

φT

Rφ )s ϑr ϑs

(ψ ′ Ru )r (ψ ′

m m ∑ ∑

(ψ ′

wT

Rw )r (ψ ′

r =1 s=1

r =1 s=1 m m 3

φT

r =1 s=1

m m ∑ ∑

N¯ φ2 = −a37

uT

(ψ ′ Ru )r (ψ ′′

(117)

r =1 s=1

− a55

m m ∑ ∑

(ψ ′ Rv )r (ψ ′ Ry )s ϑr ϑs

∑∑ r =1 s=1 m m

+ a28

Rw )s ϑr ϑs + a37

r =1 s=1

r =1 s=1 m m

N¯ φ2 = a44

wT

∑∑

uT

(ψ ′ Ru )r (ψ φ Rφ )s ϑr ϑs T

(118)

r =1 s=1

N¯ u3 = (a33 − a22 )

m m m ∑ ∑ ∑

uT

(ψ ′ Ru )r (ψ φ Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

T

(119)

r =1 s=1 p=1

−

−

m m m a22 ∑ ∑ ∑

2

(ψ y Ry )r (ψ φ Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

T

T

r =1 s=1 p=1

m m m a12 ∑ ∑ ∑

2

(ψ ′

wT

Rw )r (ψ φ Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

T

r =1 s=1 p=1

− a12

m m m ∑ ∑ ∑

vT

uT

(ψ ′ Ru )r (ψ ′ Rv )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1

+

+

m m m a11 ∑ ∑ ∑

2

uT

uT

uT

uT

vT

vT

(ψ ′ Ru )r (ψ ′ Ru )s (ψ ′ Ru )p ϑr ϑs ϑp

r =1 s=1 p=1

m m m a11 ∑ ∑ ∑

2

(ψ ′ Ru )r (ψ ′ Rv )s (ψ ′ Rv )p ϑr ϑs ϑp

r =1 s=1 p=1

− a44

m m m ∑ ∑ ∑

yT

(ψ ′ R y ) r (ψ ′

φT

Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1

+ (a55 +

m m m a18 ∑ ∑ ∑ ′ uT u φT φT ) (ψ R )r (ψ ′ Rφ )s (ψ ′ Rφ )p ϑr ϑs ϑp 2 r =1 s=1 p=1

N¯ v3 =

−

m m m a11 ∑ ∑ ∑

2

vT

uT

r =1 s=1 p=1

m m m 3a12 ∑ ∑ ∑

2

uT

(ψ ′ Rv )r (ψ ′ Ru )s (ψ ′ Ru )p ϑr ϑs ϑp vT

vT

(ψ ′ Rv )r (ψ ′ Rv )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1

+ (a22 − a33 )

m m m ∑ ∑ ∑

vT

(ψ ′ Rv )r (ψ φ Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

T

r =1 s=1 p=1

−

a33

+(

2 a18 2

m m m ∑ ∑ ∑

(ψ x Rx )r (ψ φ Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

T

T

r =1 s=1 p=1

+ a44 )

m m m ∑ ∑ ∑ r =1 s=1 p=1

vT

(ψ ′ Rv )r (ψ ′

φT

Rφ )s (ψ ′

φT

Rφ )p ϑr ϑs ϑp

(120)

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

−

+

m m m a12 ∑ ∑ ∑

2

uT

uT

T

vT

vT

vT

35

(ψ ′ Ru )r (ψ ′ Ru )s (ψ φ Rφ )p ϑr ϑs ϑp

r =1 s=1 p=1

m m m a11 ∑ ∑ ∑

2

(ψ ′ Rv )r (ψ ′ Rv )s (ψ ′ Rv )p ϑr ϑs ϑp

r =1 s=1 p=1 m m m ∑ ∑ ∑

− a55

xT

(ψ ′ Rx )r (ψ ′

φT

Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1 m m m ∑ ∑ ∑

− a56

(ψ ′

φT

Rφ )r (ψ ′

φT

Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1

N¯ w3 = − 1

N¯ x3 = −

m m m a12 ∑ ∑ ∑

2

2

N¯ x3 = −a55

T

vT

T

T

(121)

r =1 s=1 p=1

m m m a33 ∑ ∑ ∑

2

uT

(ψ ′ Ru )r (ψ φ Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp (ψ ′ Rv )r (ψ φ Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

(122)

r =1 s=1 p=1 m m m ∑ ∑ ∑

vT

(ψ ′ Rv )r (ψ ′

φT

Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

(123)

r =1 s=1 p=1 1

N¯ y3 =

m m m a22 ∑ ∑ ∑

2

2

N¯ y3 = a44

uT

(ψ ′ Ru )r (ψ φ Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

T

(124)

r =1 s=1 p=1 m m m ∑ ∑ ∑

uT

(ψ ′ Ru )r (ψ ′

φT

Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

(125)

r =1 s=1 p=1 1

N¯ φ3 = −a55

m m m ∑ ∑ ∑

vT

(ψ ′ Rv )r (ψ ′

φT

Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1

− a56

m m m ∑ ∑ ∑

vT

(ψ ′ Rv )r (ψ ′

φT

φT

R φ ) s (ψ ′

Rφ )p ϑr ϑs ϑp

r =1 s=1 p=1

−

−

−

m m m a12 ∑ ∑ ∑

2 a12 2

− a44

uT

uT

vT

vT

vT

r =1 s=1 p=1 m m m ∑ ∑ ∑

(ψ ′ Rv )r (ψ ′ Rv )s (ψ ′ Rv )p ϑr ϑs ϑp

r =1 s=1 p=1

m m m a12 ∑ ∑ ∑

2

vT

(ψ ′ Rv )r (ψ ′ Ru )s (ψ ′ Ru )p ϑr ϑs ϑp

(ψ ′

wT

uT

Rw )r (ψ ′ Ru )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1 m m m ∑ ∑ ∑

uT

φT

uT

T

uT

uT

T

vT

vT

T

vT

T

(ψ ′ Ru )r (ψ ′

yT

Rφ )s (ψ ′ Ry )p ϑr ϑs ϑp

r =1 s=1 p=1

− a22

m m m ∑ ∑ ∑

(ψ ′ Ru )r (ψ y Ry )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1

− a22

m m m ∑ ∑ ∑

(ψ ′ Ru )r (ψ ′ Ru )s (ψ φ Rφ )p ϑr ϑs ϑp

r =1 s=1 p=1

+ a22

m m m ∑ ∑ ∑

(ψ ′ Rv )r (ψ ′ Rv )s (ψ φ Rφ )p ϑr ϑs ϑp

r =1 s=1 p=1

+ a33

m m m ∑ ∑ ∑ r =1 s=1 p=1

(ψ ′ Rv )r (ψ x Rx )s (ψ φ Rφ )p ϑr ϑs ϑp T

(126)

36

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

+ a33

m m m ∑ ∑ ∑

vT

vT

T

uT

uT

T

(ψ ′ Rv )r (ψ ′ Rv )s (ψ φ Rφ )p ϑr ϑs ϑp

r =1 s=1 p=1

− a33

m m m ∑ ∑ ∑

(ψ ′ Ru )r (ψ ′ Ru )s (ψ φ Rφ )p ϑr ϑs ϑp

r =1 s=1 p=1 2

N¯ φ3 =

m m m a88 ∑ ∑ ∑

2

(ψ ′

φT

Rφ )r (ψ ′

φT

φT

Rφ )s (ψ ′

φT

Rφ )s (ψ φ Rφ )p ϑr ϑs ϑp

Rφ )p ϑr ϑs ϑp

(127)

r =1 s=1 p=1

− 2a56

m m m ∑ ∑ ∑

vT

(ψ ′ Rv )r (ψ ′

T

r =1 s=1 p=1

+ (a55 +

m m m a18 ∑ ∑ ∑ ′ φ T φ uT uT ) (ψ R )r (ψ ′ Ru )s (ψ ′ Ru )p ϑr ϑs ϑp 2 r =1 s=1 p=1

− a55

m m m ∑ ∑ ∑

xT

vT

(ψ ′ Rx )r (ψ ′ Rv )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1

+ (a44 +

m m m a18 ∑ ∑ ∑ ′ φ T φ vT vT ) (ψ R )r (ψ ′ Rv )s (ψ ′ Rv )p ϑr ϑs ϑp 2 r =1 s=1 p=1

m

− a44

m

m ∑∑∑

yT

uT

(ψ ′ Ry )r (ψ ′ Ru )s (ψ φ Rφ )p ϑr ϑs ϑp T

r =1 s=1 p=1

References Atmane, H.A., Tounsi, A., Bernard, F., 2017. Effect of thickness stretching and porosity on mechanical response of a functionally graded beams resting on elastic foundations. Int. J. Mech Mater. Des. 13 (1), 71–84. Attard, M.M., 1986. Nonlinear theory of non-uniform torsion of thin-walled open beams. Thin-Walled Struct. 4 (2), 101–134. Barmby, J.G., Cunningham, H.J., Garrick, I., 1951. Study of effects of sweep on the flutter of cantilever wings. Tech. Rep. NACA-TR-2121, National Advisory Committee for Aeronautics. Langley Aeronautical Lab., Langley Field, VA, United States. Berci, M., Righi, M., 2017. An enhanced analytical method for the subsonic indicial lift of two-dimensional aerofoils – with numerical cross-validation. Aerosp. Sci. Technol. 67, 354–365. Bhaskar, K., Librescu, L., 1995. A geometrically non-linear theory for laminated anisotropic thin-walled beams. Internat. J. Engrg. Sci. 33 (9), 1331–1344. Bisplinghoff, R., Ashley, H., Halfman, R., 2013. Aeroelasticity. Dover Books on Aeronautical Engineering. Dover Publications. Bourada, F., Bousahla, A.A., Bourada, M., Azzaz, A., Zinata, A., Tounsi, A., 2019. Dynamic investigation of porous functionally graded beam using a sinusoidal shear deformation theory. Wind Struct. 28 (1), 19–30. Carrera, E., Cinefra, M., Petrolo, M., Zappino, E., 2014. Finite Element Analysis of Structures through Unified Formulation. Wiley. Cesnik, C., Hodges, D., Patil, M., 1996. Aeroelastic analysis of composite wings. In: 37th Structure, Structural Dynamics and Materials Conference. Salt Lake City, UT, USA, pp. 1113–1123. Chaabane, L.A., Bourada, F., Sekkal, M., Zerouati, S., Zaoui, F.Z., Tounsi, A., Derras, A., Bousahla, A.A., Tounsi, A., 2019. Analytical study of bending and free vibration responses of functionally graded beams resting on elastic foundation. Struct. Eng. Mech. 71 (2), 185–196. Chikh, A., Tounsi, A., Hebali, H., Mahmoud, S., 2017. Thermal buckling analysis of cross-ply laminated plates using a simplified hsdt. Smart Struct. Syst. 19 (3), 289–297. Choo, J., Yoon, G., Song, J.-S., Kwon, T., Na, S., Qin, Z., 2012. Robust aeroelastic control of a thin-walled wing structure with model uncertainty. J. Aerosp. Eng. 25 (2), 320–333. Diederich, F., 1951. A plan-form parameter for correlating certain aerodynamic characteristics of swept wings. Tech. Rep. NACA-TN-2335, National Advisory Committee for Aeronautics. Langley Aeronautical Lab., Langley Field, VA, United States. Dimitriadis, G., 2017. Introduction to Nonlinear Aeroelasticity. In: Aerospace Series, Wiley. Dowell, E., Edwards, J., Strganac, T., 2003. Nonlinear aeroelasticity. J. Aircr. 40 (5), 857–874. Eken, S., Kaya, M.O., 2015. The effect of sweep angle on the limit cycle oscillations of aircraft wings. Adv. Aircr. Spacecr. Sci. 2 (2), 199–215. van Erp, G., Menken, C., Veldpaus, F., 1988. The non-linear flexural-torsional behaviour of straight slender elastic beams with arbitrary cross sections. Thin-Walled Struct. 6 (5), 385–404. Farsadi, T., Asadi, D., Kurtaran, H., 2019. Flutter improvement of a thin walled wing-engine system by applying curvilinear fiber path. Aerosp. Sci. Technol. 93, 105353. Farsadi, T., Hasbestan, J., 2019. Calculation of flutter and dynamic behavior of advanced composite swept wings with tapered cross section in unsteady incompressible flow. Mech. Adv. Mater. Struct. 26 (4), 314–332. Farsadi, T., Javanshir, J., 2012. Expansion of indicial function approximations for 2-d subsonic compressible aerodynamic loads. In: ASME 2012 International Mechanical Engineering Congress and Exposition. American Society of Mechanical Engineers Digital Collection, pp. 519–528. Farsadi, T., Rahmanian, M., Kayran, A., 2018. Geometrically nonlinear aeroelastic behavior of pretwisted composite wings modeled as thin walled beams. J. Fluids Struct. 83, 259–292. Farsadi, T., Sener, O., Kayran, A., 2017. Free vibration analysis of uniform and asymmetric composite pretwisted rotating thin walled beam. In: Proceedings of the ASME International Mechanical Engineering Congress and Exposition, Florida, USA, IMECE2017-70531. Fazelzadeh, S.A., Hosseini, M., 2007. Aerothermoelastic behavior of supersonic rotating thin-walled beams made of functionally graded materials. J. Fluids Struct. 23 (8), 1251–1264.

T. Farsadi, M. Rahmanian and A. Kayran / Journal of Fluids and Structures 92 (2020) 102812

37

Haddadpour, H., Kouchakzadeh, M., Shadmehri, F., 2008. Aeroelastic instability of aircraft composite wings in an incompressible flow. Compos. Struct. 83 (1), 93–99. Jung, S.N., Nagaraj, V., Chopra, I., 1999a. Assessment of composite rotor blade modeling techniques. J. Am. Helicopter Soc. 44 (3), 188–205. Jung, S.N., Nagaraj, V.T., Chopra, I., 1999b. Assessment of composite rotor blade modeling techniques. J. Am. Helicopter Soc. 44 (3), 188–205. Jung, S.N., Nagaraj, V.T., Chopra, I., 2002a. Refined structural model for thin- and thick-walled composite rotor blades. AIAA J. 40 (1), 105–116. Jung, S.N., Nagaraj, V., Chopra, I., 2002b. Refined structural model for thin-and thick-walled composite rotor blades. AIAA J. 40 (1), 105–116. Kaci, A., Houari, M.S.A., Bousahla, A.A., Tounsi, A., Mahmoud, S., 2018. Post-buckling analysis of shear-deformable composite beams using a novel simple two-unknown beam theory. Struct. Eng. Mech. 65 (5), 621–631. Kim, K., 2004. Nonlinear aeroelastic analysis of aircraft wing-with-store configurations (Ph.D. thesis). Texas A&M University. Leishman, J., 1988. Validation of approximate indicial aerodynamic functions for two- dimensional subsonic flow. J. Aircr. 25 (10), 914–922. Leishman, J., 1993. Indicial lift approximations for two-dimensional subsonic flow as obtained from oscillatory measurements. J. Aircr. 30 (3), 340–351. Leishman, G.J., 2006. Principles of Helicopter Aerodynamics with CD Extra. Cambridge university press. Librescu, L., Na, S., Qin, Z., Lee, B., 2008. Active aeroelastic control of aircraft composite wings impacted by explosive blasts. J. Sound Vib. 318 (1–2), 74–92. Librescu, L., Song, O., 1991. Behavior of thin-walled beams made of advanced composite materials and incorporating non-classical effects. Appl. Mech. Rev. 44 (11 Part 2), 174–180. Librescu, L., Song, O., 1992. On the static aeroelastic tailoring of composite aircraft swept wings modelled as thin-walled beam structures. Compos. Eng. 2 (5–7), 497–512. Librescu, L., Song, O., 2006. Thin-walled Composite Beams: Theory and Application. Springer Science & Business Media. Lin, J., Iliff, K.W., 2000. Aerodynamic lift and moment calculations using a closed-form solution of the possio equation. Tech. Rep. NASA Technical Report, NASA/TM-2000-209019. Lomax, H., 1961. Indicial aerodynamics. Tech. Rep. Technical Report, 62N10912. Lomax, H., Heaslet, M.A., Fuller, F.B., Sluder, L., 1952. Two-and three-dimensional unsteady lift problems in high-speed flight. Tech. Rep. NACA-TR-1077. Marzocca, P., Librescu, L., Chiocchia, G., 2002. Aeroelastic response of a 2-d airfoil in a compressible flow field and exposed to blast loading. Aerosp. Sci. Technol. 6 (4), 259–272. Mazelsky, B., 1952. Determination of indicial lift and moment of a two-dimensional pitching airfoil at subsonic mach numbers from oscillatory coefficients with numerical calculations for a mach number of 0.7. Tech. rep., National Aeronautics and Space Administration, Washington DC. Mazelsky, B., Drischler, J.A., 1952. Numerical determination of indicial lift and moment functions for a two-dimensional sinking and pitching airfoil at mach numbers 0.5 and 0.6. Tech. rep., National Aeronautics and Space Administration, Washington DC. Moore, D.B., 1986. A non-linear theory for the behaviour of thin-walled sections subject to combined bending and torsion. Thin-Walled Struct. 4 (6), 449–466. Na, S., Librescu, L., Marzocca, P., Yoon, G., Rubillo, C., Bong, K., 2007. Robust aeroelastic control of two-dimensional supersonic flapped wing systems. Acta Mech. 192 (1–4), 37–47. Na, S., Song, J.-S., Choo, J.-H., Qin, Z., 2011. Dynamic aeroelastic response and active control of composite thin-walled beam structures in compressible flow. J. Sound Vib. 330 (21), 4998–5013. Patil, M.J., Hodges, D.H., Cesnik, C.E., 2001. Limit-cycle oscillations in high-aspect-ratio wings. J. Fluids Struct. 15 (1), 107–132. 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. Petrolo, M., 2013. Flutter analysis of composite lifting surfaces by the 1d carrera unified formulation and the doublet lattice method. Compos. Struct. 95, 539–546. Price, S., Lee, B., Alighanbari, H., 1994. Postinstability behavior of a two-dimensional airfoil with a structural nonlinearity. J. Aircr. 31 (6), 1395–1401. Qin, Z., 2001. Vibration and aeroelasticity of advanced aircraft wings modeled as thin-walled beams–dynamics, stability and control (Ph.D. thesis). Qin, Z., Marzocca, P., Librescu, L., 2002. Aeroelastic instability and response of advanced aircraft wings at subsonic flight speeds. Aerosp. Sci. Technol. 6 (3), 195–208. Queijo, M., Wells, W.R., Keskar, D.A., 1978. Approximate indicial lift function for tapered, swept wings in incompressible flow. Tech. Rep. NASA Technical Report, Report No.: NASA TP 1241.9. Rajpal, D., De Breuker, R., 2017. Preliminary aeroelastic design framework for composite wings subjected to gust loads. In: International Forum on Aeroelasticity and Structural Dynamics, Como, Italy. Ronch, A.D., Ventura, A., Righi, M., Franciolini, M., Berci, M., Kharlamov, D., 2018. Extension of analytical indicial aerodynamics to generic trapezoidal wings in subsonic flow. Chin. J. Aeronaut. 31 (4), 617–631. Shadmehri, F., Haddadpour, H., Kouchakzadeh, M., 2007. Flexural–torsional behavior of thin-walled composite beams with closed cross-section. Thin-Walled Struct. 45 (7–8), 699–705. 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. Shirk, M.H., Hertz, T.J., Weisshaar, T.A., 1986. Aeroelastic tailoring-theory, practice, and promise. J. Aircr. 23 (1), 6–18. Sodja, J., Werter, N., Dillinger, J.K., De Breuker, R., 2016. Dynamic response of aeroelastically tailored composite wing: Analysis and experiment. In: 57th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. p. 0469. Stodieck, O., Cooper, J.E., Weaver, P.M., Kealy, P., 2013. Improved aeroelastic tailoring using tow-steered composites. Compos. Struct. 106, 703–715. Stodieck, O., Cooper, J.E., Weaver, P., Kealy, P., 2014. Optimisation of tow-steered composite wing laminates for aeroelastic tailoring. In: 55th AIAA/ASMe/ASCE/AHS/SC Structures, Structural Dynamics, and Materials Conference. p. 0343. Suresh, J., Nagaraj, V., 1996. Higher-order shear deformation theory for thin-walled composite beams. J. Aircr. 33 (5), 978–986. Tang, D., Dowell, E., 2001. Experimental and theoretical study on aeroelastic response of high-aspect-ratio wings. AIAA J. 39 (8), 1430–1441. Tang, D., Dowell, E., 2004. Effects of geometric structural nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings. J. Fluids Struct. 19 (3), 291–306. Tobak, M., 1954. On the use of the indicial-function concept in the analysis of unsteady motions of wings and wing-tail combinations. Tech. Rep. NACA-TR-1188. Wang, X., Qin, Z., 2016. Nonlinear modal interactions in composite thin-walled beam structures with simultaneous 1: 2 internal and 1: 1 external resonances. Nonlinear Dynam. 86 (2), 1381–1405. Weisshaar, T.A., 1981. Aeroelastic tailoring of forward swept composite wings. J. Aircr. 18 (8), 669–676. Zine, A., Tounsi, A., Draiche, K., Sekkal, M., Mahmoud, S., 2018. A novel higher-order shear deformation theory for bending and free vibration analysis of isotropic and multilayered plates and shells. Steel Compos. Struct. 26 (2), 125–137.